portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 1 | ! This file is part of MOM6, the Modular Ocean Model version 6. | |
| 2 | ! See the LICENSE file for licensing information. | |
| 3 | ! SPDX-License-Identifier: Apache-2.0 | |
| 4 | ||
| 5 | !> Provides the ocean stochastic equation of state | |
| 6 | module MOM_stoch_eos | |
| 7 | ||
| 8 | use MOM_diag_mediator, only : register_diag_field, post_data, diag_ctrl | |
| 9 | use MOM_error_handler, only : MOM_error, FATAL | |
| 10 | use MOM_file_parser, only : get_param, param_file_type | |
| 11 | use MOM_grid, only : ocean_grid_type | |
| 12 | use MOM_hor_index, only : hor_index_type | |
| 13 | use MOM_isopycnal_slopes, only : vert_fill_TS | |
| 14 | use MOM_random, only : PRNG, random_2d_constructor, random_2d_norm | |
| 15 | use MOM_restart, only : MOM_restart_CS, register_restart_field, is_new_run, query_initialized | |
| 16 | use MOM_time_manager, only : time_type | |
| 17 | use MOM_unit_scaling, only : unit_scale_type | |
| 18 | use MOM_variables, only : thermo_var_ptrs | |
| 19 | use MOM_verticalGrid, only : verticalGrid_type | |
| 20 | !use random_numbers_mod, only : getRandomNumbers, initializeRandomNumberStream, randomNumberStream | |
| 21 | ||
| 22 | implicit none ; private | |
| 23 | #include <MOM_memory.h> | |
| 24 | ||
| 25 | public MOM_stoch_eos_init | |
| 26 | public MOM_stoch_eos_run | |
| 27 | public stoch_EOS_register_restarts | |
| 28 | public post_stoch_EOS_diags | |
| 29 | public MOM_calc_varT | |
| 30 | ||
| 31 | !> Describes parameters of the stochastic component of the EOS | |
| 32 | !! correction, described in Stanley et al. JAMES 2020. | |
| 33 | type, public :: MOM_stoch_eos_CS ; private | |
| 34 | real, allocatable :: l2_inv(:,:) !< One over sum of the T cell side side lengths squared [L-2 ~> m-2] | |
| 35 | real, allocatable :: rgauss(:,:) !< nondimensional random Gaussian [nondim] | |
| 36 | real :: tfac = 0.27 !< Nondimensional decorrelation time factor, ~1/3.7 [nondim] | |
| 37 | real :: amplitude = 0.624499 !< Nondimensional standard deviation of Gaussian [nondim] | |
| 38 | integer :: seed !< PRNG seed | |
| 39 | type(PRNG) :: rn_CS !< PRNG control structure | |
| 40 | real, allocatable :: pattern(:,:) !< Random pattern for stochastic EOS [nondim] | |
| 41 | real, allocatable :: phi(:,:) !< temporal correlation stochastic EOS [nondim] | |
| 42 | logical :: use_stoch_eos!< If true, use the stochastic equation of state (Stanley et al. 2020) | |
| 43 | real :: stanley_coeff !< Coefficient correlating the temperature gradient | |
| 44 | !! and SGS T variance [nondim]; if <0, turn off scheme in all codes | |
| 45 | real :: stanley_a !< a in exp(aX) in stochastic coefficient [nondim] | |
| 46 | real :: kappa_smooth !< A diffusivity for smoothing T/S in vanished layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1] | |
| 47 | ||
| 48 | !>@{ Diagnostic IDs | |
| 49 | integer :: id_stoch_eos = -1, id_stoch_phi = -1, id_tvar_sgs = -1 | |
| 50 | !>@} | |
| 51 | ||
| 52 | end type MOM_stoch_eos_CS | |
| 53 | ||
| 54 | contains | |
| 55 | ||
| 56 | !> Initializes MOM_stoch_eos module, returning a logical indicating whether this module will be used. | |
| 57 | 1 | logical function MOM_stoch_eos_init(Time, G, GV, US, param_file, diag, CS, restart_CS) |
| 58 | type(time_type), intent(in) :: Time !< Time for stochastic process | |
| 59 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 60 | type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure | |
| 61 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 62 | type(param_file_type), intent(in) :: param_file !< structure indicating parameter file to parse | |
| 63 | type(diag_ctrl), target, intent(inout) :: diag !< Structure used to control diagnostics | |
| 64 | type(MOM_stoch_eos_CS), intent(inout) :: CS !< Stochastic control structure | |
| 65 | type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure. | |
| 66 | ||
| 67 | ! local variables | |
| 68 | integer :: i,j | |
| 69 | ||
| 70 | 1 | MOM_stoch_eos_init = .false. |
| 71 | ||
| 72 | 1 | CS%seed = 0 |
| 73 | ||
| 74 | call get_param(param_file, "MOM_stoch_eos", "STOCH_EOS", CS%use_stoch_eos, & | |
| 75 | "If true, stochastic perturbations are applied "//& | |
| 76 | 1 | "to the EOS in the PGF.", default=.false.) |
| 77 | call get_param(param_file, "MOM_stoch_eos", "STANLEY_COEFF", CS%stanley_coeff, & | |
| 78 | "Coefficient correlating the temperature gradient "//& | |
| 79 | 1 | "and SGS T variance.", units="nondim", default=-1.0) |
| 80 | call get_param(param_file, "MOM_stoch_eos", "STANLEY_A", CS%stanley_a, & | |
| 81 | "Coefficient a which scales chi in stochastic perturbation of the "//& | |
| 82 | "SGS T variance.", units="nondim", default=1.0, & | |
| 83 | 1 | do_not_log=((CS%stanley_coeff<0.0) .or. .not.CS%use_stoch_eos)) |
| 84 | call get_param(param_file, "MOM_stoch_eos", "KD_SMOOTH", CS%kappa_smooth, & | |
| 85 | "A diapycnal diffusivity that is used to interpolate "//& | |
| 86 | "more sensible values of T & S into thin layers.", & | |
| 87 | units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T, & | |
| 88 | 1 | do_not_log=(CS%stanley_coeff<0.0)) |
| 89 | ||
| 90 | ! Don't run anything if STANLEY_COEFF < 0 | |
| 91 | 1 | if (CS%stanley_coeff >= 0.0) then |
| 92 | 0 | if (.not.allocated(CS%pattern)) call MOM_error(FATAL, & |
| 93 | "MOM_stoch_eos_CS%pattern is not allocated when it should be, suggesting that "//& | |
| 94 | 0 | "stoch_EOS_register_restarts() has not been called before MOM_stoch_eos_init().") |
| 95 | ||
| 96 | 0 | allocate(CS%phi(G%isd:G%ied,G%jsd:G%jed), source=0.0) |
| 97 | 0 | allocate(CS%l2_inv(G%isd:G%ied,G%jsd:G%jed), source=0.0) |
| 98 | 0 | allocate(CS%rgauss(G%isd:G%ied,G%jsd:G%jed), source=0.0) |
| 99 | call get_param(param_file, "MOM_stoch_eos", "SEED_STOCH_EOS", CS%seed, & | |
| 100 | 0 | "Specfied seed for random number sequence ", default=0) |
| 101 | 0 | call random_2d_constructor(CS%rn_CS, G%HI, Time, CS%seed) |
| 102 | 0 | call random_2d_norm(CS%rn_CS, G%HI, CS%rgauss) |
| 103 | ! fill array with approximation of grid area needed for decorrelation time-scale calculation | |
| 104 | 0 | do j=G%jsc,G%jec |
| 105 | 0 | do i=G%isc,G%iec |
| 106 | 0 | CS%l2_inv(i,j) = 1.0 / ( (G%dxT(i,j)**2) + (G%dyT(i,j)**2) ) |
| 107 | enddo | |
| 108 | enddo | |
| 109 | ||
| 110 | 0 | if (.not.query_initialized(CS%pattern, "stoch_eos_pattern", restart_CS) .or. & |
| 111 | is_new_run(restart_CS)) then | |
| 112 | 0 | do j=G%jsc,G%jec ; do i=G%isc,G%iec |
| 113 | 0 | CS%pattern(i,j) = CS%amplitude*CS%rgauss(i,j) |
| 114 | enddo ; enddo | |
| 115 | endif | |
| 116 | ||
| 117 | !register diagnostics | |
| 118 | CS%id_tvar_sgs = register_diag_field('ocean_model', 'tvar_sgs', diag%axesTL, Time, & | |
| 119 | 0 | 'Parameterized SGS Temperature Variance ', 'None') |
| 120 | 0 | if (CS%use_stoch_eos) then |
| 121 | CS%id_stoch_eos = register_diag_field('ocean_model', 'stoch_eos', diag%axesT1, Time, & | |
| 122 | 0 | 'random pattern for EOS', 'None') |
| 123 | CS%id_stoch_phi = register_diag_field('ocean_model', 'stoch_phi', diag%axesT1, Time, & | |
| 124 | 0 | 'phi for EOS', 'None') |
| 125 | endif | |
| 126 | endif | |
| 127 | ||
| 128 | ! This module is only used if explicitly enabled or a positive correlation coefficient is set. | |
| 129 | 1 | MOM_stoch_eos_init = CS%use_stoch_eos .or. (CS%stanley_coeff >= 0.0) |
| 130 | ||
| 131 | 1 | end function MOM_stoch_eos_init |
| 132 | ||
| 133 | !> Register fields related to the stoch_EOS module for resarts | |
| 134 | 1 | subroutine stoch_EOS_register_restarts(HI, param_file, CS, restart_CS) |
| 135 | type(hor_index_type), intent(in) :: HI !< Horizontal index structure | |
| 136 | type(param_file_type), intent(in) :: param_file !< structure indicating parameter file to parse | |
| 137 | type(MOM_stoch_eos_CS), intent(inout) :: CS !< Stochastic control structure | |
| 138 | type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure. | |
| 139 | ||
| 140 | call get_param(param_file, "MOM_stoch_eos", "STANLEY_COEFF", CS%stanley_coeff, & | |
| 141 | "Coefficient correlating the temperature gradient "//& | |
| 142 | 1 | "and SGS T variance.", units="nondim", default=-1.0, do_not_log=.true.) |
| 143 | ||
| 144 | 1 | if (CS%stanley_coeff >= 0.0) then |
| 145 | 0 | allocate(CS%pattern(HI%isd:HI%ied,HI%jsd:HI%jed), source=0.0) |
| 146 | call register_restart_field(CS%pattern, "stoch_eos_pattern", .false., restart_CS, & | |
| 147 | 0 | "Random pattern for stoch EOS", "nondim") |
| 148 | endif | |
| 149 | ||
| 150 | 1 | end subroutine stoch_EOS_register_restarts |
| 151 | ||
| 152 | !> Generates a pattern in space and time for the ocean stochastic equation of state | |
| 153 | 0 | subroutine MOM_stoch_eos_run(G, u, v, delt, Time, CS) |
| 154 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 155 | real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), & | |
| 156 | intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1]. | |
| 157 | real, dimension(SZI_(G),SZJB_(G),SZK_(G)), & | |
| 158 | intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]. | |
| 159 | real, intent(in) :: delt !< Time step size for AR1 process [T ~> s]. | |
| 160 | type(time_type), intent(in) :: Time !< Time for stochastic process | |
| 161 | type(MOM_stoch_eos_CS), intent(inout) :: CS !< Stochastic control structure | |
| 162 | ||
| 163 | ! local variables | |
| 164 | real :: ubar, vbar ! Averaged velocities [L T-1 ~> m s-1] | |
| 165 | real :: phi ! A temporal correlation factor [nondim] | |
| 166 | integer :: i, j | |
| 167 | ||
| 168 | ! Return without doing anything if this capability is not enabled. | |
| 169 | 0 | if (.not.CS%use_stoch_eos) return |
| 170 | ||
| 171 | 0 | call random_2d_constructor(CS%rn_CS, G%HI, Time, CS%seed) |
| 172 | 0 | call random_2d_norm(CS%rn_CS, G%HI, CS%rgauss) |
| 173 | ||
| 174 | ! advance AR(1) | |
| 175 | 0 | do j=G%jsc,G%jec |
| 176 | 0 | do i=G%isc,G%iec |
| 177 | 0 | ubar = 0.5*(u(I,j,1)*G%mask2dCu(I,j)+u(I-1,j,1)*G%mask2dCu(I-1,j)) |
| 178 | 0 | vbar = 0.5*(v(i,J,1)*G%mask2dCv(i,J)+v(i,J-1,1)*G%mask2dCv(i,J-1)) |
| 179 | 0 | phi = exp(-delt*CS%tfac * sqrt(((ubar**2) + (vbar**2))*CS%l2_inv(i,j))) |
| 180 | 0 | CS%pattern(i,j) = phi*CS%pattern(i,j) + CS%amplitude*sqrt(1-phi**2)*CS%rgauss(i,j) |
| 181 | 0 | CS%phi(i,j) = phi |
| 182 | enddo | |
| 183 | enddo | |
| 184 | ||
| 185 | end subroutine MOM_stoch_eos_run | |
| 186 | ||
| 187 | !> Write out any diagnostics related to this module. | |
| 188 | 0 | subroutine post_stoch_EOS_diags(CS, tv, diag) |
| 189 | type(MOM_stoch_eos_CS), intent(in) :: CS !< Stochastic control structure | |
| 190 | type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure | |
| 191 | type(diag_ctrl), intent(inout) :: diag !< Structure to control diagnostics | |
| 192 | ||
| 193 | 0 | if (CS%id_stoch_eos > 0) call post_data(CS%id_stoch_eos, CS%pattern, diag) |
| 194 | 0 | if (CS%id_stoch_phi > 0) call post_data(CS%id_stoch_phi, CS%phi, diag) |
| 195 | 0 | if (CS%id_tvar_sgs > 0) call post_data(CS%id_tvar_sgs, tv%varT, diag) |
| 196 | ||
| 197 | 0 | end subroutine post_stoch_EOS_diags |
| 198 | ||
| 199 | !> Computes a parameterization of the SGS temperature variance | |
| 200 | 0 | subroutine MOM_calc_varT(G, GV, US, h, tv, CS, dt) |
| 201 | type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. | |
| 202 | type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure | |
| 203 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 204 | real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & | |
| 205 | intent(in) :: h !< Layer thickness [H ~> m] | |
| 206 | type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure | |
| 207 | type(MOM_stoch_eos_CS), intent(inout) :: CS !< Stochastic control structure | |
| 208 | real, intent(in) :: dt !< Time increment [T ~> s] | |
| 209 | ||
| 210 | ! local variables | |
| 211 | real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: & | |
| 212 | 0 | T, & !> The temperature (or density) [C ~> degC], with the values in |
| 213 | !! in massless layers filled vertically by diffusion. | |
| 214 | 0 | S !> The filled salinity [S ~> ppt], with the values in |
| 215 | !! in massless layers filled vertically by diffusion. | |
| 216 | real :: hl(5) !> Copy of local stencil of H [H ~> m] | |
| 217 | real :: dTdi2, dTdj2 !> Differences in T variance [C2 ~> degC2] | |
| 218 | integer :: i, j, k | |
| 219 | ||
| 220 | ! Nothing happens if a negative correlation coefficient is set. | |
| 221 | 0 | if (CS%stanley_coeff < 0.0) return |
| 222 | ||
| 223 | ! This block does a thickness weighted variance calculation and helps control for | |
| 224 | ! extreme gradients along layers which are vanished against topography. It is | |
| 225 | ! still a poor approximation in the interior when coordinates are strongly tilted. | |
| 226 | 0 | if (.not. associated(tv%varT)) allocate(tv%varT(G%isd:G%ied, G%jsd:G%jed, GV%ke), source=0.0) |
| 227 | !$omp target enter data map(to: h, tv%T, tv%S) | |
| 228 | !$omp target enter data map(alloc: T, S) | |
| 229 | 0 | call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth*dt, T, S, G, GV, US, halo_here=1, larger_h_denom=.true.) |
| 230 | !$omp target exit data map(from: T, S) | |
| 231 | !$omp target exit data map(release: h, tv%T, tv%S) | |
| 232 | ||
| 233 | 0 | do k=1,G%ke |
| 234 | 0 | do j=G%jsc,G%jec |
| 235 | 0 | do i=G%isc,G%iec |
| 236 | 0 | hl(1) = h(i,j,k) * G%mask2dT(i,j) |
| 237 | 0 | hl(2) = h(i-1,j,k) * G%mask2dCu(I-1,j) |
| 238 | 0 | hl(3) = h(i+1,j,k) * G%mask2dCu(I,j) |
| 239 | 0 | hl(4) = h(i,j-1,k) * G%mask2dCv(i,J-1) |
| 240 | 0 | hl(5) = h(i,j+1,k) * G%mask2dCv(i,J) |
| 241 | ||
| 242 | ! SGS variance in i-direction [C2 ~> degC2] | |
| 243 | dTdi2 = ( ( G%mask2dCu(I ,j) * (G%IdxCu(I ,j) * ( T(i+1,j,k) - T(i,j,k) )) & | |
| 244 | + G%mask2dCu(I-1,j) * (G%IdxCu(I-1,j) * ( T(i,j,k) - T(i-1,j,k) )) & | |
| 245 | 0 | ) * G%dxT(i,j) * 0.5 )**2 |
| 246 | ! SGS variance in j-direction [C2 ~> degC2] | |
| 247 | dTdj2 = ( ( G%mask2dCv(i,J ) * (G%IdyCv(i,J ) * ( T(i,j+1,k) - T(i,j,k) )) & | |
| 248 | + G%mask2dCv(i,J-1) * (G%IdyCv(i,J-1) * ( T(i,j,k) - T(i,j-1,k) )) & | |
| 249 | 0 | ) * G%dyT(i,j) * 0.5 )**2 |
| 250 | 0 | tv%varT(i,j,k) = CS%stanley_coeff * ( dTdi2 + dTdj2 ) |
| 251 | ! Turn off scheme near land | |
| 252 | 0 | tv%varT(i,j,k) = tv%varT(i,j,k) * (minval(hl) / (maxval(hl) + GV%H_subroundoff)) |
| 253 | enddo | |
| 254 | enddo | |
| 255 | enddo | |
| 256 | ! if stochastic, perturb | |
| 257 | 0 | if (CS%use_stoch_eos) then |
| 258 | 0 | do k=1,G%ke |
| 259 | 0 | do j=G%jsc,G%jec |
| 260 | 0 | do i=G%isc,G%iec |
| 261 | 0 | tv%varT(i,j,k) = exp(CS%stanley_a * CS%pattern(i,j)) * tv%varT(i,j,k) |
| 262 | enddo | |
| 263 | enddo | |
| 264 | enddo | |
| 265 | endif | |
| 266 | end subroutine MOM_calc_varT | |
| 267 | ||
| 268 | 0 | end module MOM_stoch_eos |