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, 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 integer :: i,j
69
70 mom_stoch_eos_init = .false.
71
72 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 "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 "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 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 do_not_log=(cs%stanley_coeff<0.0))
89
90 ! Don't run anything if STANLEY_COEFF < 0
91 if (cs%stanley_coeff >= 0.0) then
92 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 "stoch_EOS_register_restarts() has not been called before MOM_stoch_eos_init().")
95
96 allocate(cs%phi(g%isd:g%ied,g%jsd:g%jed), source=0.0)
97 allocate(cs%l2_inv(g%isd:g%ied,g%jsd:g%jed), source=0.0)
98 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 "Specfied seed for random number sequence ", default=0)
101 call random_2d_constructor(cs%rn_CS, g%HI, time, cs%seed)
102 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 do j=g%jsc,g%jec
105 do i=g%isc,g%iec
106 cs%l2_inv(i,j) = 1.0 / ( (g%dxT(i,j)**2) + (g%dyT(i,j)**2) )
107 enddo
108 enddo
109
110 if (.not.query_initialized(cs%pattern, "stoch_eos_pattern", restart_cs) .or. &
111 is_new_run(restart_cs)) then
112 do j=g%jsc,g%jec ; do i=g%isc,g%iec
113 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 'Parameterized SGS Temperature Variance ', 'None')
120 if (cs%use_stoch_eos) then
121 cs%id_stoch_eos = register_diag_field('ocean_model', 'stoch_eos', diag%axesT1, time, &
122 'random pattern for EOS', 'None')
123 cs%id_stoch_phi = register_diag_field('ocean_model', 'stoch_phi', diag%axesT1, time, &
124 '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 mom_stoch_eos_init = cs%use_stoch_eos .or. (cs%stanley_coeff >= 0.0)
130
131end function mom_stoch_eos_init
132
133!> Register fields related to the stoch_EOS module for resarts
134subroutine 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 "and SGS T variance.", units="nondim", default=-1.0, do_not_log=.true.)
143
144 if (cs%stanley_coeff >= 0.0) then
145 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 "Random pattern for stoch EOS", "nondim")
148 endif
149
150end subroutine stoch_eos_register_restarts
151
152!> Generates a pattern in space and time for the ocean stochastic equation of state
153subroutine 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 if (.not.cs%use_stoch_eos) return
170
171 call random_2d_constructor(cs%rn_CS, g%HI, time, cs%seed)
172 call random_2d_norm(cs%rn_CS, g%HI, cs%rgauss)
173
174 ! advance AR(1)
175 do j=g%jsc,g%jec
176 do i=g%isc,g%iec
177 ubar = 0.5*(u(i,j,1)*g%mask2dCu(i,j)+u(i-1,j,1)*g%mask2dCu(i-1,j))
178 vbar = 0.5*(v(i,j,1)*g%mask2dCv(i,j)+v(i,j-1,1)*g%mask2dCv(i,j-1))
179 phi = exp(-delt*cs%tfac * sqrt(((ubar**2) + (vbar**2))*cs%l2_inv(i,j)))
180 cs%pattern(i,j) = phi*cs%pattern(i,j) + cs%amplitude*sqrt(1-phi**2)*cs%rgauss(i,j)
181 cs%phi(i,j) = phi
182 enddo
183 enddo
184
185end subroutine mom_stoch_eos_run
186
187!> Write out any diagnostics related to this module.
188subroutine 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 if (cs%id_stoch_eos > 0) call post_data(cs%id_stoch_eos, cs%pattern, diag)
194 if (cs%id_stoch_phi > 0) call post_data(cs%id_stoch_phi, cs%phi, diag)
195 if (cs%id_tvar_sgs > 0) call post_data(cs%id_tvar_sgs, tv%varT, diag)
196
197end subroutine post_stoch_eos_diags
198
199!> Computes a parameterization of the SGS temperature variance
200subroutine 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 t, & !> The temperature (or density) [C ~> degC], with the values in
213 !! in massless layers filled vertically by diffusion.
214 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 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 if (.not. associated(tv%varT)) allocate(tv%varT(g%isd:g%ied, g%jsd:g%jed, gv%ke), source=0.0)
227 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.)
228
229 do k=1,g%ke
230 do j=g%jsc,g%jec
231 do i=g%isc,g%iec
232 hl(1) = h(i,j,k) * g%mask2dT(i,j)
233 hl(2) = h(i-1,j,k) * g%mask2dCu(i-1,j)
234 hl(3) = h(i+1,j,k) * g%mask2dCu(i,j)
235 hl(4) = h(i,j-1,k) * g%mask2dCv(i,j-1)
236 hl(5) = h(i,j+1,k) * g%mask2dCv(i,j)
237
238 ! SGS variance in i-direction [C2 ~> degC2]
239 dtdi2 = ( ( g%mask2dCu(i ,j) * (g%IdxCu(i ,j) * ( t(i+1,j,k) - t(i,j,k) )) &
240 + g%mask2dCu(i-1,j) * (g%IdxCu(i-1,j) * ( t(i,j,k) - t(i-1,j,k) )) &
241 ) * g%dxT(i,j) * 0.5 )**2
242 ! SGS variance in j-direction [C2 ~> degC2]
243 dtdj2 = ( ( g%mask2dCv(i,j ) * (g%IdyCv(i,j ) * ( t(i,j+1,k) - t(i,j,k) )) &
244 + g%mask2dCv(i,j-1) * (g%IdyCv(i,j-1) * ( t(i,j,k) - t(i,j-1,k) )) &
245 ) * g%dyT(i,j) * 0.5 )**2
246 tv%varT(i,j,k) = cs%stanley_coeff * ( dtdi2 + dtdj2 )
247 ! Turn off scheme near land
248 tv%varT(i,j,k) = tv%varT(i,j,k) * (minval(hl) / (maxval(hl) + gv%H_subroundoff))
249 enddo
250 enddo
251 enddo
252 ! if stochastic, perturb
253 if (cs%use_stoch_eos) then
254 do k=1,g%ke
255 do j=g%jsc,g%jec
256 do i=g%isc,g%iec
257 tv%varT(i,j,k) = exp(cs%stanley_a * cs%pattern(i,j)) * tv%varT(i,j,k)
258 enddo
259 enddo
260 enddo
261 endif
262end subroutine mom_calc_vart
263
264end module mom_stoch_eos