MOM_stochastics.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!> Top-level module for the MOM6 ocean model in coupled mode.
7
8! This is the top level module for the MOM6 ocean model. It contains routines
9! for initialization, update, and writing restart of stochastic physics. This
10! particular version wraps all of the calls for MOM6 in the calls that had
11! been used for MOM4.
12!
13use mom_coms, only : get_pelist
14use mom_debugging, only : hchksum, uvchksum, qchksum
15use mom_diag_mediator, only : register_diag_field, diag_ctrl, time_type, post_data
17use mom_domains, only : pass_var, pass_vector, corner, scalar_pair
18use mom_domains, only : root_pe, num_pes
19use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
21use mom_file_parser, only : get_param, log_version, close_param_file, param_file_type
22use mom_grid, only : ocean_grid_type
26use mom_eos, only : calculate_density, eos_domain
28use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain
29use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain
30
31#include <MOM_memory.h>
32
33implicit none ; private
34
36
37!> This control structure holds parameters for the MOM_stochastics module
38type, public:: stochastic_cs
39 logical :: do_sppt !< If true, stochastically perturb the diabatic
40 logical :: do_skeb !< If true, stochastically perturb the horizontal velocity
41 logical :: skeb_use_gm !< If true, adds GM work to the amplitude of SKEBS
42 logical :: skeb_use_frict !< If true, adds viscous dissipation rate to the amplitude of SKEBS
43 logical :: pert_epbl !< If true, then randomly perturb the KE dissipation and genration terms
44 integer :: id_sppt_wts = -1 !< Diagnostic id for SPPT
45 integer :: id_skeb_wts = -1 !< Diagnostic id for SKEB
46 integer :: id_skebu = -1 !< Diagnostic id for SKEB
47 integer :: id_skebv = -1 !< Diagnostic id for SKEB
48 integer :: id_diss = -1 !< Diagnostic id for SKEB
49 integer :: skeb_npass = -1 !< number of passes of the 9-point smoother for the dissipation estimate
50 integer :: id_psi = -1 !< Diagnostic id for SPPT
51 integer :: id_epbl1_wts = -1 !< Diagnostic id for epbl generation perturbation
52 integer :: id_epbl2_wts = -1 !< Diagnostic id for epbl dissipation perturbation
53 integer :: id_skeb_taperu = -1 !< Diagnostic id for u taper of SKEB velocity increment
54 integer :: id_skeb_taperv = -1 !< Diagnostic id for v taper of SKEB velocity increment
55 real :: skeb_gm_coef !< If skeb_use_gm is true, then skeb_gm_coef * GM_work is added to the
56 !! dissipation rate used to set the amplitude of SKEBS [nondim]
57 real :: skeb_frict_coef !< If skeb_use_frict is true, then skeb_gm_coef * GM_work is added to the
58 !! dissipation rate used to set the amplitude of SKEBS [nondim]
59 real, allocatable :: skeb_diss(:,:,:) !< Dissipation rate used to set amplitude of SKEBS [L2 T-3 ~> m2 s-3]
60 !! Index into this at h points.
61 integer :: answer_date !< The vintage of the order of arithmetic in the stochastics
62 !! calculations. Values below 20250701 recover the answers from
63 !! early in 2025, while higher values use expressions that have been
64 !! refactored for rotational symmetry, including with FMAs enabled.
65
66 ! stochastic patterns
67 real, allocatable :: sppt_wts(:,:) !< Random pattern for ocean SPPT
68 !! tendencies with a number between 0 and 2 [nondim]
69 real, allocatable :: skeb_wts(:,:) !< Random pattern of lengthscales for ocean SKEB in mks units [m]
70 !! Note that SKEB_wts is set via external code in mks units.
71 real, allocatable :: epbl1_wts(:,:) !< Random pattern for K.E. generation [nondim]
72 real, allocatable :: epbl2_wts(:,:) !< Random pattern for K.E. dissipation [nondim]
73 type(time_type), pointer :: time !< Pointer to model time (needed for sponges)
74 type(diag_ctrl), pointer :: diag=>null() !< A structure that is used to regulate the
75
76 ! Taper array to smoothly zero out the SKEBS velocity increment near land
77 real, allocatable :: tapercu(:,:) !< Taper applied to u component of stochastic
78 !! velocity increment range [0,1], [nondim]
79 real, allocatable :: tapercv(:,:) !< Taper applied to v component of stochastic
80 !! velocity increment range [0,1], [nondim]
81
82 ! Weights for smoothing skeb_diss
83 real allocable_, dimension(NIMEM_,NJMEM_) :: isum_area_wts, & !< One over the 3x3 sum of area_wt [L-2 ~> m-2]
84 area_wt !< Masked h cell areas. [L2 ~> m2]
85
86end type stochastic_cs
87
88contains
89
90!> This subroutine initializes the stochastics physics control structure.
91subroutine stochastics_init(dt, grid, GV, US, CS, param_file, diag, Time)
92 real, intent(in) :: dt !< time step [T ~> s]
93 type(ocean_grid_type), intent(in) :: grid !< horizontal grid information
94 type(verticalgrid_type), intent(in) :: gv !< vertical grid structure
95 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
96 type(stochastic_cs), pointer, intent(inout) :: cs !< stochastic control structure
97 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
98 type(diag_ctrl), target, intent(inout) :: diag !< structure to regulate diagnostic output
99 type(time_type), target :: time !< model time
100
101 ! Local variables
102 integer, allocatable :: pelist(:) ! list of pes for this instance of the ocean
103 integer :: mom_comm ! list of pes for this instance of the ocean
104 integer :: num_procs ! number of processors to pass to stochastic physics
105 integer :: iret ! return code from stochastic physics
106 integer :: pe_zero ! root pe
107 integer :: nxt, nxb ! number of x-points including halo
108 integer :: nyt, nyb ! number of y-points including halo
109 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
110 integer :: i, j, k ! loop indices
111 real :: tmp(grid%isdb:grid%iedb,grid%jsdb:grid%jedb) ! Used to construct tapers [nondim]
112 integer :: taper_width ! Width (in cells) of the taper that brings the stochastic velocity
113 ! increments to 0 at the boundary.
114 real :: sum_area_wts ! A rotationally symmetric sum of the surrounding area weights
115 ! that are used to filter skeb_diss [L2 ~> m2]
116
117 ! This include declares and sets the variable "version".
118# include "version_variable.h"
119 character(len=40) :: mdl = "ocean_stochastics_init" ! This module's name.
120
121 call calltree_enter("stochastic_init(), MOM_stochastics.F90")
122 if (associated(cs)) then
123 call mom_error(warning, "MOM_stochastics_init called with an "// &
124 "associated control structure.")
125 return
126 else ; allocate(cs) ; endif
127
128 cs%Time => time
129 cs%diag => diag
130
131 ! Read all relevant parameters and write them to the model log.
132 call log_version(param_file, mdl, version, "")
133
134 ! get number of processors and PE list for stochastic physics initialization
135 call get_param(param_file, mdl, "DO_SPPT", cs%do_sppt, &
136 "If true, then stochastically perturb the thermodynamic "//&
137 "tendencies of T,S, and h. Amplitude and correlations are "//&
138 "controlled by the nam_stoch namelist in the UFS model only.", &
139 default=.false.)
140 call get_param(param_file, mdl, "DO_SKEB", cs%do_skeb, &
141 "If true, then stochastically perturb the currents "//&
142 "using the stochastic kinetic energy backscatter scheme.",&
143 default=.false.)
144 call get_param(param_file, mdl, "SKEB_NPASS", cs%skeb_npass, &
145 "number of passes of a 9-point smoother of the "//&
146 "dissipation estimate.", default=3, do_not_log=.not.cs%do_skeb)
147 call get_param(param_file, mdl, "SKEB_TAPER_WIDTH", taper_width, &
148 "number of cells over which the stochastic velocity increment "//&
149 "is tapered to zero.", default=4, do_not_log=.not.cs%do_skeb)
150 call get_param(param_file, mdl, "SKEB_USE_GM", cs%skeb_use_gm, &
151 "If true, adds GM work rate to the SKEBS amplitude.", &
152 default=.false., do_not_log=.not.cs%do_skeb)
153 if ((.not. cs%do_skeb) .and. (cs%skeb_use_gm)) call mom_error(fatal, "If SKEB_USE_GM is True "//&
154 "then DO_SKEB must also be True.")
155 call get_param(param_file, mdl, "SKEB_GM_COEF", cs%skeb_gm_coef, &
156 "Fraction of GM work that is added to backscatter rate.", &
157 units="nondim", default=0.0, do_not_log=.not.cs%skeb_use_gm)
158 call get_param(param_file, mdl, "SKEB_USE_FRICT", cs%skeb_use_frict, &
159 "If true, adds horizontal friction dissipation rate "//&
160 "to the SKEBS amplitude.", default=.false., do_not_log=.not.cs%do_skeb)
161 if ((.not. cs%do_skeb) .and. (cs%skeb_use_frict)) call mom_error(fatal, "If SKEB_USE_FRICT is "//&
162 "True then DO_SKEB must also be True.")
163 call get_param(param_file, mdl, "SKEB_FRICT_COEF", cs%skeb_frict_coef, &
164 "Fraction of horizontal friction work that is added to backscatter rate.", &
165 units="nondim", default=0.0, do_not_log=.not.cs%skeb_use_frict)
166 call get_param(param_file, mdl, "PERT_EPBL", cs%pert_epbl, &
167 "If true, then stochastically perturb the kinetic energy "//&
168 "production and dissipation terms. Amplitude and correlations are "//&
169 "controlled by the nam_stoch namelist in the UFS model only.", &
170 default=.false.)
171 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
172 "This sets the default value for the various _ANSWER_DATE parameters.", &
173 default=99991231, do_not_log=.true.)
174 call get_param(param_file, mdl, "STOCHASTICS_ANSWER_DATE", cs%answer_date, &
175 "The vintage of the order of arithmetic in the stochastics calculations. "//&
176 "Values below 20250701 recover the answers from early in 2025, while higher "//&
177 "values use expressions that have been refactored for rotational symmetry.", &
178 default=20250101) !### Change to: default=default_answer_date)
179
180 if (cs%do_sppt .OR. cs%pert_epbl .OR. cs%do_skeb) then
181 num_procs = num_pes()
182 allocate(pelist(num_procs))
183 call get_pelist(pelist,commid = mom_comm)
184 pe_zero = root_pe()
185 nxt = grid%ied - grid%isd + 1
186 nyt = grid%jed - grid%jsd + 1
187 nxb = grid%iedB - grid%isdB + 1
188 nyb = grid%jedB - grid%jsdB + 1
189 call init_stochastic_physics_ocn(dt*us%T_to_s, grid%geoLonT, grid%geoLatT, nxt, nyt, gv%ke, &
190 grid%geoLonBu, grid%geoLatBu, nxb, nyb, &
191 cs%pert_epbl, cs%do_sppt, cs%do_skeb, pe_zero, mom_comm, iret)
192 if (iret/=0) then
193 call mom_error(fatal, "call to init_stochastic_physics_ocn failed")
194 return
195 endif
196
197 if ((cs%do_sppt) .or. (cs%do_skeb)) allocate(cs%sppt_wts(grid%isd:grid%ied,grid%jsd:grid%jed))
198 if (cs%do_skeb) allocate(cs%skeb_wts(grid%isdB:grid%iedB,grid%jsdB:grid%jedB))
199 if (cs%do_skeb) allocate(cs%skeb_diss(grid%isd:grid%ied,grid%jsd:grid%jed,gv%ke), source=0.)
200 if ((cs%pert_epbl) .or. (cs%do_skeb)) then
201 allocate(cs%epbl1_wts(grid%isd:grid%ied,grid%jsd:grid%jed))
202 allocate(cs%epbl2_wts(grid%isd:grid%ied,grid%jsd:grid%jed))
203 endif
204 endif
205
206 cs%id_sppt_wts = register_diag_field('ocean_model', 'sppt_pattern', cs%diag%axesT1, time, &
207 'random pattern for sppt', 'None')
208 cs%id_skeb_wts = register_diag_field('ocean_model', 'skeb_pattern', cs%diag%axesB1, time, &
209 'random pattern for skeb', 'm', conversion=1.0) ! SKEB_wts is set in external code in mks units of [m]
210 cs%id_epbl1_wts = register_diag_field('ocean_model', 'epbl1_wts', cs%diag%axesT1, time, &
211 'random pattern for KE generation', 'None')
212 cs%id_epbl2_wts = register_diag_field('ocean_model', 'epbl2_wts', cs%diag%axesT1, time, &
213 'random pattern for KE dissipation', 'None')
214 cs%id_skebu = register_diag_field('ocean_model', 'skebu', cs%diag%axesCuL, time, &
215 'zonal current perts', 'm s-1', conversion=us%L_T_to_m_s)
216 cs%id_skebv = register_diag_field('ocean_model', 'skebv', cs%diag%axesCvL, time, &
217 'zonal current perts', 'm s-1', conversion=us%L_T_to_m_s)
218 cs%id_diss = register_diag_field('ocean_model', 'skeb_amp', cs%diag%axesTL, time, &
219 'SKEB amplitude', 'm s-1', conversion=us%L_T_to_m_s)
220 cs%id_psi = register_diag_field('ocean_model', 'psi', cs%diag%axesBL, time, &
221 'stream function', 'm2 s-1', conversion=us%L_T_to_m_s*us%L_to_m)
222 cs%id_skeb_taperu = register_static_field('ocean_model', 'skeb_taper_u', cs%diag%axesCu1, &
223 'SKEB taper u', 'None', conversion=1.0, interp_method='none')
224 cs%id_skeb_taperv = register_static_field('ocean_model', 'skeb_taper_v', cs%diag%axesCv1, &
225 'SKEB taper v', 'None', conversion=1.0, interp_method='none')
226
227 ! Initialize the "taper" fields. These fields multiply the components of the stochastic
228 ! velocity increment in such a way as to smoothly taper them to zero at land boundaries.
229 if ((cs%do_skeb) .or. (cs%id_skeb_taperu > 0) .or. (cs%id_skeb_taperv > 0)) then
230 allocate(cs%taperCu(grid%IsdB:grid%IedB,grid%jsd:grid%jed))
231 allocate(cs%taperCv(grid%isd:grid%ied,grid%JsdB:grid%JedB))
232 ! Initialize taper from land mask
233 do j=grid%jsd,grid%jed ; do i=grid%isdB,grid%iedB
234 cs%taperCu(i,j) = grid%mask2dCu(i,j)
235 enddo ; enddo
236 do j=grid%jsdB,grid%jedB ; do i=grid%isd,grid%ied
237 cs%taperCv(i,j) = grid%mask2dCv(i,j)
238 enddo ; enddo
239 ! Extend taper land
240 do k=1,(taper_width / 2)
241 do j=grid%jsc-1,grid%jec+1 ; do i=grid%iscB-1,grid%iecB+1
242 tmp(i,j) = minval(cs%taperCu(i-1:i+1,j-1:j+1))
243 enddo ; enddo
244 do j=grid%jsc,grid%jec ; do i=grid%iscB,grid%iecB
245 cs%taperCu(i,j) = minval(tmp(i-1:i+1,j-1:j+1))
246 enddo ; enddo
247 do j=grid%jscB-1,grid%jecB+1 ; do i=grid%isc-1,grid%iec+1
248 tmp(i,j) = minval(cs%taperCv(i-1:i+1,j-1:j+1))
249 enddo ; enddo
250 do j=grid%jscB,grid%jecB ; do i=grid%isc,grid%iec
251 cs%taperCv(i,j) = minval(tmp(i-1:i+1,j-1:j+1))
252 enddo ; enddo
253 ! Update halo
254 call pass_vector(cs%taperCu, cs%taperCv, grid%Domain, scalar_pair)
255 enddo
256 ! Smooth tapers. Each call smooths twice.
257 do k=1,(taper_width - (taper_width/2))
258 call smooth_x9_uv(grid, cs%taperCu, cs%taperCv, zero_land=.true.)
259 call pass_vector(cs%taperCu, cs%taperCv, grid%Domain, scalar_pair)
260 enddo
261 endif
262
263 !call uvchksum("SKEB taper [uv]", CS%taperCu, CS%taperCv, grid%HI)
264
265 if (cs%id_skeb_taperu > 0) call post_data(cs%id_skeb_taperu, cs%taperCu, cs%diag, .true.)
266 if (cs%id_skeb_taperv > 0) call post_data(cs%id_skeb_taperv, cs%taperCv, cs%diag, .true.)
267
268 ! Initialize the smoothing weights
269 if ((cs%do_skeb) .and. cs%skeb_npass >= 1) then
270 alloc_(cs%area_wt(grid%isd:grid%ied,grid%jsd:grid%jed)) ; cs%area_wt(:,:) = 0.0
271 alloc_(cs%Isum_area_wts(grid%isd:grid%ied,grid%jsd:grid%jed)) ; cs%Isum_area_wts(:,:) = 0.0
272 do j=grid%jsc-2,grid%jec+2 ; do i=grid%isc-2,grid%iec+2
273 cs%area_wt(i,j) = grid%mask2dT(i,j)*grid%areaT(i,j)
274 enddo ; enddo
275 do j=grid%jsc-1,grid%jec+1 ; do i=grid%isc-1,grid%iec+1
276 sum_area_wts = cs%area_wt(i,j) + &
277 (((cs%area_wt(i-1,j) + cs%area_wt(i+1,j)) + (cs%area_wt(i,j-1) + cs%area_wt(i,j+1))) + &
278 ((cs%area_wt(i-1,j-1) + cs%area_wt(i+1,j+1)) + (cs%area_wt(i-1,j+1) + cs%area_wt(i+1,j-1))))
279 cs%Isum_area_wts(i,j) = 1.0 / (sum_area_wts + 1.e-16*us%m_to_L**2)
280 enddo ; enddo
281 endif
282
283 if (cs%do_sppt .OR. cs%pert_epbl .OR. cs%do_skeb) &
284 call mom_mesg(' === COMPLETED MOM STOCHASTIC INITIALIZATION =====')
285
286 call calltree_leave("stochastic_init(), MOM_stochastics.F90")
287
288end subroutine stochastics_init
289
290!> Advances the stochastic patterns one time step
291subroutine update_stochastics(CS)
292 type(stochastic_cs), intent(inout) :: cs !< diabatic control structure
293 call calltree_enter("update_stochastics(), MOM_stochastics.F90")
294
295! update stochastic physics patterns before running next time-step
296 call run_stochastic_physics_ocn(cs%sppt_wts, cs%skeb_wts, cs%epbl1_wts, cs%epbl2_wts)
297
298 call calltree_leave("update_stochastics(), MOM_stochastics.F90")
299
300end subroutine update_stochastics
301
302!> Adds a stochastic increment (backscatter) to the input velocity field
303subroutine apply_skeb(grid, GV, US, CS, uc, vc, thickness, tv, dt, Time_end)
304
305 type(ocean_grid_type), intent(in) :: grid !< ocean grid structure
306 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid
307 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
308 type(stochastic_cs), intent(inout) :: cs !< stochastic control structure
309 real, dimension(SZIB_(grid),SZJ_(grid),SZK_(GV)), intent(inout) :: uc !< zonal velocity [L T-1 ~> m s-1]
310 real, dimension(SZI_(grid),SZJB_(grid),SZK_(GV)), intent(inout) :: vc !< meridional velocity [L T-1 ~> m s-1]
311 real, dimension(SZI_(grid),SZJ_(grid),SZK_(GV)), intent(in) :: thickness !< thickness [H ~> m or kg m-2]
312 type(thermo_var_ptrs), intent(in) :: tv !< points to thermodynamic fields
313 real, intent(in) :: dt !< time increment [T ~> s]
314 type(time_type), intent(in) :: time_end !< Time at the end of the interval
315
316 ! local variables
317 real, dimension(SZIB_(grid),SZJB_(grid),SZK_(GV)) :: psi !< Streamfunction for stochastic velocity increments
318 !! [L2 T-1 ~> m2 s-1]
319 real, dimension(SZIB_(grid),SZJ_(grid) ,SZK_(GV)) :: ustar !< Stochastic u velocity increment [L T-1 ~> m s-1]
320 real, dimension(SZI_(grid) ,SZJB_(grid),SZK_(GV)) :: vstar !< Stochastic v velocity increment [L T-1 ~> m s-1]
321 real, dimension(SZI_(grid),SZJ_(grid)) :: diss_tmp !< Temporary array used in smoothing skeb_diss
322 !! [L2 T-3 ~> m2 s-2]
323 real, dimension(3,3) :: local_weights !< 3x3 stencil weights used in smoothing skeb_diss
324 !! [L2 ~> m2]
325
326 real :: shr ! Horizonal shear [T-1 ~> s-1]
327 real :: ten ! Horizonal tension of the flow [T-1 ~> s-1]
328 real :: tot ! The magnitude of the combined shear and tension [T-1 ~> s-1]
329 real :: kh ! A smooothing factor [nondim]
330 real :: sum_wtd_skeb_diss ! The rotationally symmetric sum of the surrounding values of skeb times
331 ! the area weights used to filter skeb_diss [L4 T-3 ~> m4 s-3]
332 integer :: i, j, k, iter
333 integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
334
335 call calltree_enter("apply_skeb(), MOM_stochastics.F90")
336
337 if ((.not. cs%skeb_use_gm) .and. (.not. cs%skeb_use_frict)) then
338 ! fill in halos with zeros
339 do k=1,gv%ke
340 do j=grid%jsd,grid%jed ; do i=grid%isd,grid%ied
341 cs%skeb_diss(i,j,k) = 0.0
342 enddo ; enddo
343 enddo
344
345 ! kh needs to be scaled
346 kh = 1.0 !(120*111)**2
347 if (cs%answer_date < 20250701) then
348 do k=1,gv%ke
349 do j=grid%jsc,grid%jec ; do i=grid%isc,grid%iec
350 ! Shear in [T-1 ~> s-1]
351 shr = (vc(i,j,k)-vc(i-1,j,k)) * grid%mask2dCv(i,j)*grid%mask2dCv(i-1,j)*grid%IdxCv(i,j) + &
352 (uc(i,j,k)-uc(i,j-1,k)) * grid%mask2dCu(i,j)*grid%mask2dCu(i,j-1)*grid%IdyCu(i,j)
353 ! Tension in [T-1 ~> s-1]
354 ten = (vc(i,j,k)-vc(i-1,j,k)) * grid%mask2dCv(i,j)*grid%mask2dCv(i-1,j)*grid%IdyCv(i,j) + &
355 (uc(i,j,k)-uc(i,j-1,k)) * grid%mask2dCu(i,j)*grid%mask2dCu(i,j-1)*grid%IdxCu(i,j)
356
357 tot = sqrt( shr**2 + ten**2 ) * grid%mask2dT(i,j)
358 cs%skeb_diss(i,j,k) = tot**3 * kh * grid%areaT(i,j) !!**2
359 enddo ; enddo
360 enddo
361 else ! This version has parentheses to preserve rotational symmetry when FMAs are enabled.
362 do k=1,gv%ke
363 do j=grid%jsc,grid%jec ; do i=grid%isc,grid%iec
364 ! Shear in [T-1 ~> s-1]
365 shr = ((vc(i,j,k)-vc(i-1,j,k)) * grid%mask2dCv(i,j)*grid%mask2dCv(i-1,j)*grid%IdxCv(i,j)) + &
366 ((uc(i,j,k)-uc(i,j-1,k)) * grid%mask2dCu(i,j)*grid%mask2dCu(i,j-1)*grid%IdyCu(i,j))
367 ! Tension in [T-1 ~> s-1]
368 ten = ((vc(i,j,k)-vc(i-1,j,k)) * grid%mask2dCv(i,j)*grid%mask2dCv(i-1,j)*grid%IdyCv(i,j)) + &
369 ((uc(i,j,k)-uc(i,j-1,k)) * grid%mask2dCu(i,j)*grid%mask2dCu(i,j-1)*grid%IdxCu(i,j))
370
371 tot = sqrt( shr**2 + ten**2 ) * grid%mask2dT(i,j)
372 cs%skeb_diss(i,j,k) = tot**3 * kh * grid%areaT(i,j) !!**2
373 enddo ; enddo
374 enddo
375 endif
376 endif ! Sets CS%skeb_diss in [L2 T-3 ~> m2 s-3] without GM or FrictWork
377
378 ! smooth dissipation skeb_npass times
379 do iter=1,cs%skeb_npass
380 if (mod(iter,2) == 1) call pass_var(cs%skeb_diss, grid%domain)
381 do k=1,gv%ke
382 if (cs%answer_date < 20250701) then
383 ! Do the filter with expressions that do not preserve rotational symmetry.
384 do j=grid%jsc-1,grid%jec+1 ; do i=grid%isc-1,grid%iec+1
385 local_weights(:,:) = cs%area_wt(i-1:i+1,j-1:j+1)
386 diss_tmp(i,j) = sum(local_weights(:,:)*cs%skeb_diss(i-1:i+1,j-1:j+1,k)) / &
387 (sum(local_weights) + 1.e-16*us%m_to_L**2)
388 enddo ; enddo
389 else
390 ! This spatial filter preserves rotational symmeetry (including with FMAs), but is
391 ! mathematically equivalent to the older sum-based form above
392 do j=grid%jsc-1,grid%jec+1 ; do i=grid%isc-1,grid%iec+1
393 sum_wtd_skeb_diss = cs%skeb_diss(i,j,k) * cs%area_wt(i+1,j) + &
394 ((( (cs%skeb_diss(i-1,j,k) * cs%area_wt(i-1,j)) + (cs%skeb_diss(i+1,j,k) * cs%area_wt(i+1,j)) ) + &
395 ( (cs%skeb_diss(i,j-1,k) * cs%area_wt(i,j-1)) + (cs%skeb_diss(i,j+1,k) * cs%area_wt(i,j+1)) )) + &
396 (( (cs%skeb_diss(i-1,j-1,k) * cs%area_wt(i-1,j-1)) + (cs%skeb_diss(i-1,j-1,k) * cs%area_wt(i+1,j+1)) ) + &
397 ( (cs%skeb_diss(i-1,j+1,k) * cs%area_wt(i-1,j+1)) + (cs%skeb_diss(i+1,j-1,k) * cs%area_wt(i+1,j-1)) )))
398 diss_tmp(i,j) = sum_wtd_skeb_diss * cs%Isum_area_wts(i,j)
399 enddo ; enddo
400 endif
401 do j=grid%jsc-1,grid%jec+1 ; do i=grid%isc-1,grid%iec+1
402 cs%skeb_diss(i,j,k) = grid%mask2dT(i,j) * diss_tmp(i,j)
403 enddo ; enddo
404 enddo
405 enddo
406 call pass_var(cs%skeb_diss, grid%domain)
407
408 ! call hchksum(CS%skeb_diss, "SKEB DISS", grid%HI, haloshift=2, unscale=US%L_T_to_m_s**2*US%s_to_T)
409 ! call qchksum(CS%skeb_wts, "SKEB WTS", grid%HI, haloshift=1) ! SKEB_wts comes in from external code in mks units.
410
411 do k=1,gv%ke
412 do j=grid%jscB-1,grid%jecB ; do i=grid%iscB-1,grid%iecB
413 ! psi has units of [L2 T-1 ~> m2 s-1] because skeb_wts is in mks units of [m].
414 psi(i,j,k) = sqrt(0.25 * dt * max((cs%skeb_diss(i ,j ,k) + cs%skeb_diss(i+1,j+1,k)) + &
415 (cs%skeb_diss(i ,j+1,k) + cs%skeb_diss(i+1,j ,k)), 0.) ) &
416 * us%m_to_L*cs%skeb_wts(i,j)
417 enddo ; enddo
418 enddo
419 !call qchksum(psi,"SKEB PSI", grid%HI, haloshift=1, unscale=US%L_T_to_m_s*US%L_to_m)
420 !call pass_var(psi, grid%domain, position=CORNER)
421 do k=1,gv%ke
422 do j=grid%jsc,grid%jec ; do i=grid%iscB,grid%iecB
423 ustar(i,j,k) = - (psi(i,j,k) - psi(i,j-1,k)) * cs%taperCu(i,j) * grid%IdyCu(i,j)
424 uc(i,j,k) = uc(i,j,k) + ustar(i,j,k)
425 enddo ; enddo
426 do j=grid%jscB,grid%jecB ; do i=grid%isc,grid%iec
427 vstar(i,j,k) = (psi(i,j,k) - psi(i-1,j,k)) * cs%taperCv(i,j) * grid%IdxCv(i,j)
428 vc(i,j,k) = vc(i,j,k) + vstar(i,j,k)
429 enddo ; enddo
430 enddo
431
432 !call uvchksum("SKEB increment [uv]", ustar, vstar, grid%HI, unscale=US%L_T_to_m_s)
433
434 call enable_averages(dt, time_end, cs%diag)
435 if (cs%id_diss > 0) then
436 call post_data(cs%id_diss, sqrt(dt * max(cs%skeb_diss(:,:,:), 0.)), cs%diag)
437 endif
438 if (cs%id_skeb_wts > 0) then
439 call post_data(cs%id_skeb_wts, cs%skeb_wts, cs%diag)
440 endif
441 if (cs%id_skebu > 0) then
442 call post_data(cs%id_skebu, ustar(:,:,:), cs%diag)
443 endif
444 if (cs%id_skebv > 0) then
445 call post_data(cs%id_skebv, vstar(:,:,:), cs%diag)
446 endif
447 if (cs%id_psi > 0) then
448 call post_data(cs%id_psi, psi(:,:,:), cs%diag)
449 endif
450 call disable_averaging(cs%diag)
451 cs%skeb_diss(:,:,:) = 0.0 ! Must zero before next time step.
452
453 call calltree_leave("apply_skeb(), MOM_stochastics.F90")
454
455end subroutine apply_skeb
456
457!> Apply a 9-point smoothing filter twice to a pair of velocity components to reduce
458!! horizontal two-grid-point noise.
459!! Note that this subroutine does not conserve angular momentum, so don't use it
460!! in situations where you need conservation. Also note that it assumes that the
461!! input fields have valid values in the first two halo points upon entry.
462subroutine smooth_x9_uv(G, field_u, field_v, zero_land)
463 type(ocean_grid_type), intent(in) :: G !< Ocean grid
464 real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: field_u !< u-point field to be smoothed [arbitrary]
465 real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: field_v !< v-point field to be smoothed [arbitrary]
466 logical, optional, intent(in) :: zero_land !< If present and false, return the average
467 !! of the surrounding ocean points when
468 !! smoothing, otherwise use a value of 0 for
469 !! land points and include them in the averages.
470
471 ! Local variables.
472 real :: fu_prev(SZIB_(G),SZJ_(G)) ! The value of the u-point field at the previous iteration [arbitrary]
473 real :: fv_prev(SZI_(G),SZJB_(G)) ! The value of the v-point field at the previous iteration [arbitrary]
474 real :: Iwts ! The inverse of the sum of the weights [nondim]
475 logical :: zero_land_val ! The value of the zero_land optional argument or .true. if it is absent.
476 integer :: i, j, s, is, ie, js, je, Isq, Ieq, Jsq, Jeq
477
478 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
479 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
480
481 zero_land_val = .true. ; if (present(zero_land)) zero_land_val = zero_land
482
483 do s=1,0,-1
484 fu_prev(:,:) = field_u(:,:)
485 ! apply smoothing on field_u using rotationally symmetric expressions.
486 do j=js-s,je+s ; do i=isq-s,ieq+s ; if (g%mask2dCu(i,j) > 0.0) then
487 iwts = 0.0625
488 if (.not. zero_land_val) &
489 iwts = 1.0 / ( (4.0*g%mask2dCu(i,j) + &
490 ( 2.0*((g%mask2dCu(i-1,j) + g%mask2dCu(i+1,j)) + &
491 (g%mask2dCu(i,j-1) + g%mask2dCu(i,j+1))) + &
492 ((g%mask2dCu(i-1,j-1) + g%mask2dCu(i+1,j+1)) + &
493 (g%mask2dCu(i-1,j+1) + g%mask2dCu(i+1,j-1))) ) ) + 1.0e-16 )
494 field_u(i,j) = iwts * ( 4.0*g%mask2dCu(i,j) * fu_prev(i,j) &
495 + (2.0*((g%mask2dCu(i-1,j) * fu_prev(i-1,j) + g%mask2dCu(i+1,j) * fu_prev(i+1,j)) + &
496 (g%mask2dCu(i,j-1) * fu_prev(i,j-1) + g%mask2dCu(i,j+1) * fu_prev(i,j+1))) &
497 + ((g%mask2dCu(i-1,j-1) * fu_prev(i-1,j-1) + g%mask2dCu(i+1,j+1) * fu_prev(i+1,j+1)) + &
498 (g%mask2dCu(i-1,j+1) * fu_prev(i-1,j+1) + g%mask2dCu(i+1,j-1) * fu_prev(i-1,j-1))) ))
499 endif ; enddo ; enddo
500
501 fv_prev(:,:) = field_v(:,:)
502 ! apply smoothing on field_v using rotationally symmetric expressions.
503 do j=jsq-s,jeq+s ; do i=is-s,ie+s ; if (g%mask2dCv(i,j) > 0.0) then
504 iwts = 0.0625
505 if (.not. zero_land_val) &
506 iwts = 1.0 / ( (4.0*g%mask2dCv(i,j) + &
507 ( 2.0*((g%mask2dCv(i-1,j) + g%mask2dCv(i+1,j)) + &
508 (g%mask2dCv(i,j-1) + g%mask2dCv(i,j+1))) + &
509 ((g%mask2dCv(i-1,j-1) + g%mask2dCv(i+1,j+1)) + &
510 (g%mask2dCv(i-1,j+1) + g%mask2dCv(i+1,j-1))) ) ) + 1.0e-16 )
511 field_v(i,j) = iwts * ( 4.0*g%mask2dCv(i,j) * fv_prev(i,j) &
512 + (2.0*((g%mask2dCv(i-1,j) * fv_prev(i-1,j) + g%mask2dCv(i+1,j) * fv_prev(i+1,j)) + &
513 (g%mask2dCv(i,j-1) * fv_prev(i,j-1) + g%mask2dCv(i,j+1) * fv_prev(i,j+1))) &
514 + ((g%mask2dCv(i-1,j-1) * fv_prev(i-1,j-1) + g%mask2dCv(i+1,j+1) * fv_prev(i+1,j+1)) + &
515 (g%mask2dCv(i-1,j+1) * fv_prev(i-1,j+1) + g%mask2dCv(i+1,j-1) * fv_prev(i-1,j-1))) ))
516 endif ; enddo ; enddo
517 enddo
518
519end subroutine smooth_x9_uv
520!> \namespace mom_stochastics
521!!
522!! This file contains subroutines that implement some stochastic parameterizations in MOM6.
523!! SPPT perturbations of the tendencies of S and T are turned on using <code>DO_SPPT=True</code>.
524!! Stochastic perturbations in ePBL are turned on using <code>PERT_EPBL=True</code>.
525!! Stochastic kinetic energy backscatter (SKEB) via the Stochastic GM+E scheme is turned on using
526!! <code>DO_SKEB=True</code>. For all three schemes the spatial and temporal correlation structure
527!! of the associated random fields is controlled from the <code>nam_stochy</code> namelist used by
528!! the external <code>stochastic_physics</code> package, which is called by subroutines in this
529!! module.
530!!
531!! The SKEB backscatter can be set in a variety of ways. If <code>SKEB_USE_GM=True</code> then
532!! <code>SKEB_GM_COEF</code> times the GM work rate will be added to the backscatter rate. (The
533!! vertical structure for this component of backscatter is the so-called EBT struct.) If
534!! <code>SKEB_USE_FRICT=True</code> then <code>SKEB_FRICT_COEF</code> times the work rate from
535!! lateral viscosity will be added to the backscatter rate. The code uses the total contribution
536!! from Laplacian and biharmonic viscosities as computed within the horizontal viscosity module.
537!! If neither <code>SKEB_USE_GM</code> nor <code>SKEB_USE_FRICT</code> is true, then the code
538!! computes the dissipation rate as if it came from a lateral harmonic viscosity with
539!! coefficient 1 (MKS units). The only thoroughly tested SKEB option at this point is
540!! <code>SKEB_USE_GM</code>.
541!!
542!! The contributions to the backscatter rate are smoothed before use. One smoothing pass uses a
543!! 3x3 moving average with weights proportional to the h-cell areas. The number of smoothing passes
544!! is controlled by <code>SKEB_NPASS</code>.
545!!
546!! A taper is applied to the SKEB velocity increments (equivalently to the SKEB stochastic forcing).
547!! The taper zeros out the increments near land cells. The width of this taper can be controlled using
548!! <code>SKEB_TAPER_WIDTH</code>.
549end module mom_stochastics