MOM_stoch_eos.F90

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
6module mom_stoch_eos
7
9use mom_error_handler, only : mom_error, fatal
10use mom_file_parser, only : get_param, log_version, param_file_type
11use mom_grid, only : ocean_grid_type
15use mom_restart, only : mom_restart_cs, register_restart_field, is_new_run, query_initialized
16use mom_time_manager, only : time_type
20!use random_numbers_mod, only : getRandomNumbers, initializeRandomNumberStream, randomNumberStream
21
22implicit none ; private
23#include <MOM_memory.h>
24
29public mom_calc_vart
30
31!> Describes parameters of the stochastic component of the EOS
32!! correction, described in Stanley et al. JAMES 2020.
33type, 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
52end type mom_stoch_eos_cs
53
54contains
55
56!> Initializes MOM_stoch_eos module, returning a logical indicating whether this module will be used.
57logical 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 ! This include declares and sets the variable "version".
69# include "version_variable.h"
70 integer :: i,j
71
72 mom_stoch_eos_init = .false.
73
74 cs%seed = 0
75
76 call log_version(param_file, "MOM_stoch_eos", version, "")
77 call get_param(param_file, "MOM_stoch_eos", "STOCH_EOS", cs%use_stoch_eos, &
78 "If true, computes stochastic perturbations that can be applied "//&
79 "to the EOS in various places.", default=.false.)
80 call get_param(param_file, "MOM_stoch_eos", "STANLEY_COEFF", cs%stanley_coeff, &
81 "Coefficient correlating the temperature gradient "//&
82 "and SGS T variance.", units="nondim", default=-1.0)
83 if ((cs%stanley_coeff < 0.0) .and. cs%use_stoch_eos) call mom_error(fatal, &
84 "STANLEY_COEFF must be set >= 0 if STOCH_EOS is true.")
85 call get_param(param_file, "MOM_stoch_eos", "STANLEY_A", cs%stanley_a, &
86 "Coefficient a which scales chi in stochastic perturbation of the "//&
87 "SGS T variance.", units="nondim", default=1.0, &
88 do_not_log=.not.cs%use_stoch_eos)
89 call get_param(param_file, "MOM_stoch_eos", "KD_SMOOTH", cs%kappa_smooth, &
90 "A diapycnal diffusivity that is used to interpolate "//&
91 "more sensible values of T & S into thin layers.", &
92 units="m2 s-1", default=1.0e-6, scale=gv%m2_s_to_HZ_T, &
93 do_not_log=.not.cs%use_stoch_eos)
94
95 ! Don't run anything if STANLEY_COEFF < 0
96 if (cs%stanley_coeff >= 0.0) then
97 if (.not.allocated(cs%pattern)) call mom_error(fatal, &
98 "MOM_stoch_eos_CS%pattern is not allocated when it should be, suggesting that "//&
99 "stoch_EOS_register_restarts() has not been called before MOM_stoch_eos_init().")
100
101 allocate(cs%phi(g%isd:g%ied,g%jsd:g%jed), source=0.0)
102 allocate(cs%l2_inv(g%isd:g%ied,g%jsd:g%jed), source=0.0)
103 allocate(cs%rgauss(g%isd:g%ied,g%jsd:g%jed), source=0.0)
104 call get_param(param_file, "MOM_stoch_eos", "SEED_STOCH_EOS", cs%seed, &
105 "Specfied seed for random number sequence ", default=0)
106 call random_2d_constructor(cs%rn_CS, g%HI, time, cs%seed)
107 call random_2d_norm(cs%rn_CS, g%HI, cs%rgauss)
108 ! fill array with approximation of grid area needed for decorrelation time-scale calculation
109 do j=g%jsc,g%jec
110 do i=g%isc,g%iec
111 cs%l2_inv(i,j) = 1.0 / ( (g%dxT(i,j)**2) + (g%dyT(i,j)**2) )
112 enddo
113 enddo
114
115 if (.not.query_initialized(cs%pattern, "stoch_eos_pattern", restart_cs) .or. &
116 is_new_run(restart_cs)) then
117 do j=g%jsc,g%jec ; do i=g%isc,g%iec
118 cs%pattern(i,j) = cs%amplitude*cs%rgauss(i,j)
119 enddo ; enddo
120 endif
121
122 !register diagnostics
123 cs%id_tvar_sgs = register_diag_field('ocean_model', 'tvar_sgs', diag%axesTL, time, &
124 'Parameterized SGS Temperature Variance ', 'None')
125 if (cs%use_stoch_eos) then
126 cs%id_stoch_eos = register_diag_field('ocean_model', 'stoch_eos', diag%axesT1, time, &
127 'random pattern for EOS', 'None')
128 cs%id_stoch_phi = register_diag_field('ocean_model', 'stoch_phi', diag%axesT1, time, &
129 'phi for EOS', 'None')
130 endif
131 endif
132
133 ! This module is only used if explicitly enabled or a positive correlation coefficient is set.
134 mom_stoch_eos_init = cs%use_stoch_eos .or. (cs%stanley_coeff >= 0.0)
135
136end function mom_stoch_eos_init
137
138!> Register fields related to the stoch_EOS module for resarts
139subroutine stoch_eos_register_restarts(HI, param_file, CS, restart_CS)
140 type(hor_index_type), intent(in) :: hi !< Horizontal index structure
141 type(param_file_type), intent(in) :: param_file !< structure indicating parameter file to parse
142 type(mom_stoch_eos_cs), intent(inout) :: cs !< Stochastic control structure
143 type(mom_restart_cs), pointer :: restart_cs !< A pointer to the restart control structure.
144
145 call get_param(param_file, "MOM_stoch_eos", "STANLEY_COEFF", cs%stanley_coeff, &
146 "Coefficient correlating the temperature gradient "//&
147 "and SGS T variance.", units="nondim", default=-1.0, do_not_log=.true.)
148
149 if (cs%stanley_coeff >= 0.0) then
150 allocate(cs%pattern(hi%isd:hi%ied,hi%jsd:hi%jed), source=0.0)
151 call register_restart_field(cs%pattern, "stoch_eos_pattern", .false., restart_cs, &
152 "Random pattern for stoch EOS", "nondim")
153 endif
154
155end subroutine stoch_eos_register_restarts
156
157!> Generates a pattern in space and time for the ocean stochastic equation of state
158subroutine mom_stoch_eos_run(G, u, v, delt, Time, CS)
159 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
160 real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
161 intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
162 real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
163 intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
164 real, intent(in) :: delt !< Time step size for AR1 process [T ~> s].
165 type(time_type), intent(in) :: time !< Time for stochastic process
166 type(mom_stoch_eos_cs), intent(inout) :: cs !< Stochastic control structure
167
168 ! local variables
169 real :: ubar, vbar ! Averaged velocities [L T-1 ~> m s-1]
170 real :: phi ! A temporal correlation factor [nondim]
171 integer :: i, j
172
173 ! Return without doing anything if this capability is not enabled.
174 if (.not.cs%use_stoch_eos) return
175
176 call random_2d_constructor(cs%rn_CS, g%HI, time, cs%seed)
177 call random_2d_norm(cs%rn_CS, g%HI, cs%rgauss)
178
179 ! advance AR(1)
180 do j=g%jsc,g%jec
181 do i=g%isc,g%iec
182 ubar = 0.5*(u(i,j,1)*g%mask2dCu(i,j)+u(i-1,j,1)*g%mask2dCu(i-1,j))
183 vbar = 0.5*(v(i,j,1)*g%mask2dCv(i,j)+v(i,j-1,1)*g%mask2dCv(i,j-1))
184 phi = exp(-delt*cs%tfac * sqrt(((ubar**2) + (vbar**2))*cs%l2_inv(i,j)))
185 cs%pattern(i,j) = phi*cs%pattern(i,j) + cs%amplitude*sqrt(1-phi**2)*cs%rgauss(i,j)
186 cs%phi(i,j) = phi
187 enddo
188 enddo
189
190end subroutine mom_stoch_eos_run
191
192!> Write out any diagnostics related to this module.
193subroutine post_stoch_eos_diags(CS, tv, diag)
194 type(mom_stoch_eos_cs), intent(in) :: cs !< Stochastic control structure
195 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
196 type(diag_ctrl), intent(inout) :: diag !< Structure to control diagnostics
197
198 if (cs%id_stoch_eos > 0) call post_data(cs%id_stoch_eos, cs%pattern, diag)
199 if (cs%id_stoch_phi > 0) call post_data(cs%id_stoch_phi, cs%phi, diag)
200 if (cs%id_tvar_sgs > 0) call post_data(cs%id_tvar_sgs, tv%varT, diag)
201
202end subroutine post_stoch_eos_diags
203
204!> Computes a parameterization of the SGS temperature variance
205subroutine mom_calc_vart(G, GV, US, h, tv, CS, dt)
206 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
207 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
208 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
209 real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
210 intent(in) :: h !< Layer thickness [H ~> m]
211 type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure
212 type(mom_stoch_eos_cs), intent(inout) :: cs !< Stochastic control structure
213 real, intent(in) :: dt !< Time increment [T ~> s]
214
215 ! local variables
216 real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: &
217 t, & !> The temperature (or density) [C ~> degC], with the values in
218 !! in massless layers filled vertically by diffusion.
219 s !> The filled salinity [S ~> ppt], with the values in
220 !! in massless layers filled vertically by diffusion.
221 real :: hl(5) !> Copy of local stencil of H [H ~> m]
222 real :: dtdi2, dtdj2 !> Differences in T variance [C2 ~> degC2]
223 integer :: i, j, k
224
225 ! Nothing happens if a negative correlation coefficient is set.
226 if (cs%stanley_coeff < 0.0) return
227
228 ! This block does a thickness weighted variance calculation and helps control for
229 ! extreme gradients along layers which are vanished against topography. It is
230 ! still a poor approximation in the interior when coordinates are strongly tilted.
231 if (.not. associated(tv%varT)) allocate(tv%varT(g%isd:g%ied, g%jsd:g%jed, gv%ke), source=0.0)
232 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.)
233
234 do k=1,g%ke
235 do j=g%jsc,g%jec
236 do i=g%isc,g%iec
237 hl(1) = h(i,j,k) * g%mask2dT(i,j)
238 hl(2) = h(i-1,j,k) * g%mask2dCu(i-1,j)
239 hl(3) = h(i+1,j,k) * g%mask2dCu(i,j)
240 hl(4) = h(i,j-1,k) * g%mask2dCv(i,j-1)
241 hl(5) = h(i,j+1,k) * g%mask2dCv(i,j)
242
243 ! SGS variance in i-direction [C2 ~> degC2]
244 dtdi2 = ( ( g%mask2dCu(i ,j) * (g%IdxCu(i ,j) * ( t(i+1,j,k) - t(i,j,k) )) &
245 + g%mask2dCu(i-1,j) * (g%IdxCu(i-1,j) * ( t(i,j,k) - t(i-1,j,k) )) &
246 ) * g%dxT(i,j) * 0.5 )**2
247 ! SGS variance in j-direction [C2 ~> degC2]
248 dtdj2 = ( ( g%mask2dCv(i,j ) * (g%IdyCv(i,j ) * ( t(i,j+1,k) - t(i,j,k) )) &
249 + g%mask2dCv(i,j-1) * (g%IdyCv(i,j-1) * ( t(i,j,k) - t(i,j-1,k) )) &
250 ) * g%dyT(i,j) * 0.5 )**2
251 tv%varT(i,j,k) = cs%stanley_coeff * ( dtdi2 + dtdj2 )
252 ! Turn off scheme near land
253 tv%varT(i,j,k) = tv%varT(i,j,k) * (minval(hl) / (maxval(hl) + gv%H_subroundoff))
254 enddo
255 enddo
256 enddo
257 ! if stochastic, perturb
258 if (cs%use_stoch_eos) then
259 do k=1,g%ke
260 do j=g%jsc,g%jec
261 do i=g%isc,g%iec
262 tv%varT(i,j,k) = exp(cs%stanley_a * cs%pattern(i,j)) * tv%varT(i,j,k)
263 enddo
264 enddo
265 enddo
266 endif
267end subroutine mom_calc_vart
268
269!> \namespace mom_stoch_eos
270!!
271!! This module provides the foundation of the Stanley parameterization (\cite stanley2020) for correcting the
272!! computation of density. Density is not a prognostic variable in MOM6; it is computed for various purposes
273!! in various places. The correction to this calculation provided by this module has been implemented
274!! in some places where density is used, but not all.
275!!
276!! To use the correction, first set <code>STOCH_EOS=True</code>. Then, choose the constant c from (25) of
277!! \cite stanley2020. This is controlled using <code>STANLEY_COEFF</code>. Setting a negative value will
278!! result in an error. \cite stanley2020 found a value of 0.2 offline, coarsening from 0.1 to 1 degree
279!! resolution. \cite kenigson2022 proposed a value of 0.5 in a 2/3 degree resolution model.
280!!
281!! Whether the correction is deterministic or stochastic can be controlled using the variable
282!! <code>STANLEY_A</code>. Setting this to 0.0 uses the deterministic version, while a value of 1.0 produces
283!! the stochastic version. Reducing from 1 to 0 smoothly transitions from stochastic to deterministic.
284!!
285!! To turn the correction on in various parts of the code, use
286!! - <code>USE_STANLEY_PGF=True</code> for the pressure gradient force (cf. \cite kenigson2022)
287!! - <code>USE_STANLEY_ISO=True</code> to correct the computation of isopycnal slopes (used in many places)
288!! - <code>USE_STANLEY_GM=True</code> to use the parameterization within GM (cf. \cite agarwal2023)
289!! - <code>USE_STANLEY_ML=True</code> to use the parameterization within the mixed-layer restratification
290!! parameterization. It applies to both the OM4 and Bodner schemes. (cf. \cite agarwal2023)
291!!
292!! For ensemble simulations, the random number generator seed can be controlled using the parameter
293!! <code>SEED_STOCH_EOS</code>
294
295end module mom_stoch_eos