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_debugging, only : hchksum, uvchksum, qchksum
14use mom_diag_mediator, only : register_diag_field, diag_ctrl, time_type, post_data
16use mom_grid, only : ocean_grid_type
18use mom_domains, only : pass_var, pass_vector, corner, scalar_pair
20use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
22use mom_file_parser, only : get_param, log_version, close_param_file, param_file_type
23use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain
24use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain
25use mom_domains, only : root_pe, num_pes
26use mom_coms, only : get_pelist
27use mom_eos, only : calculate_density, eos_domain
29
30#include <MOM_memory.h>
31
32implicit none ; private
33
35
36!> This control structure holds parameters for the MOM_stochastics module
37type, public:: stochastic_cs
38 logical :: do_sppt !< If true, stochastically perturb the diabatic
39 logical :: do_skeb !< If true, stochastically perturb the horizontal velocity
40 logical :: skeb_use_gm !< If true, adds GM work to the amplitude of SKEBS
41 logical :: skeb_use_frict !< If true, adds viscous dissipation rate to the amplitude of SKEBS
42 logical :: pert_epbl !< If true, then randomly perturb the KE dissipation and genration terms
43 integer :: id_sppt_wts = -1 !< Diagnostic id for SPPT
44 integer :: id_skeb_wts = -1 !< Diagnostic id for SKEB
45 integer :: id_skebu = -1 !< Diagnostic id for SKEB
46 integer :: id_skebv = -1 !< Diagnostic id for SKEB
47 integer :: id_diss = -1 !< Diagnostic id for SKEB
48 integer :: skeb_npass = -1 !< number of passes of the 9-point smoother for the dissipation estimate
49 integer :: id_psi = -1 !< Diagnostic id for SPPT
50 integer :: id_epbl1_wts = -1 !< Diagnostic id for epbl generation perturbation
51 integer :: id_epbl2_wts = -1 !< Diagnostic id for epbl dissipation perturbation
52 integer :: id_skeb_taperu = -1 !< Diagnostic id for u taper of SKEB velocity increment
53 integer :: id_skeb_taperv = -1 !< Diagnostic id for v taper of SKEB velocity increment
54 real :: skeb_gm_coef !< If skeb_use_gm is true, then skeb_gm_coef * GM_work is added to the
55 !! dissipation rate used to set the amplitude of SKEBS [nondim]
56 real :: skeb_frict_coef !< If skeb_use_frict is true, then skeb_gm_coef * GM_work is added to the
57 !! dissipation rate used to set the amplitude of SKEBS [nondim]
58 real, allocatable :: skeb_diss(:,:,:) !< Dissipation rate used to set amplitude of SKEBS [L2 T-3 ~> m2 s-3]
59 !! Index into this at h points.
60 ! stochastic patterns
61 real, allocatable :: sppt_wts(:,:) !< Random pattern for ocean SPPT
62 !! tendencies with a number between 0 and 2 [nondim]
63 real, allocatable :: skeb_wts(:,:) !< Random pattern for ocean SKEB [nondim]
64 real, allocatable :: epbl1_wts(:,:) !< Random pattern for K.E. generation [nondim]
65 real, allocatable :: epbl2_wts(:,:) !< Random pattern for K.E. dissipation [nondim]
66 type(time_type), pointer :: time !< Pointer to model time (needed for sponges)
67 type(diag_ctrl), pointer :: diag=>null() !< A structure that is used to regulate the
68
69 ! Taper array to smoothly zero out the SKEBS velocity increment near land
70 real, allocatable :: tapercu(:,:) !< Taper applied to u component of stochastic
71 !! velocity increment range [0,1], [nondim]
72 real, allocatable :: tapercv(:,:) !< Taper applied to v component of stochastic
73 !! velocity increment range [0,1], [nondim]
74
75end type stochastic_cs
76
77contains
78
79!! This subroutine initializes the stochastics physics control structure.
80subroutine stochastics_init(dt, grid, GV, CS, param_file, diag, Time)
81 real, intent(in) :: dt !< time step [T ~> s]
82 type(ocean_grid_type), intent(in) :: grid !< horizontal grid information
83 type(verticalgrid_type), intent(in) :: gv !< vertical grid structure
84 type(stochastic_cs), pointer, intent(inout) :: cs !< stochastic control structure
85 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
86 type(diag_ctrl), target, intent(inout) :: diag !< structure to regulate diagnostic output
87 type(time_type), target :: time !< model time
88
89 ! Local variables
90 integer, allocatable :: pelist(:) ! list of pes for this instance of the ocean
91 integer :: mom_comm ! list of pes for this instance of the ocean
92 integer :: num_procs ! number of processors to pass to stochastic physics
93 integer :: iret ! return code from stochastic physics
94 integer :: pe_zero ! root pe
95 integer :: nxt, nxb ! number of x-points including halo
96 integer :: nyt, nyb ! number of y-points including halo
97 integer :: i, j, k ! loop indices
98 real :: tmp(grid%isdb:grid%iedb,grid%jsdb:grid%jedb) ! Used to construct tapers
99 integer :: taper_width ! Width (in cells) of the taper that brings the stochastic velocity
100 ! increments to 0 at the boundary.
101
102 ! This include declares and sets the variable "version".
103# include "version_variable.h"
104 character(len=40) :: mdl = "ocean_stochastics_init" ! This module's name.
105
106 call calltree_enter("stochastic_init(), MOM_stochastics.F90")
107 if (associated(cs)) then
108 call mom_error(warning, "MOM_stochastics_init called with an "// &
109 "associated control structure.")
110 return
111 else ; allocate(cs) ; endif
112
113 cs%Time => time
114 cs%diag => diag
115
116 ! Read all relevant parameters and write them to the model log.
117 call log_version(param_file, mdl, version, "")
118
119 ! get number of processors and PE list for stochastic physics initialization
120 call get_param(param_file, mdl, "DO_SPPT", cs%do_sppt, &
121 "If true, then stochastically perturb the thermodynamic "//&
122 "tendencies of T,S, and h. Amplitude and correlations are "//&
123 "controlled by the nam_stoch namelist in the UFS model only.", &
124 default=.false.)
125 call get_param(param_file, mdl, "DO_SKEB", cs%do_skeb, &
126 "If true, then stochastically perturb the currents "//&
127 "using the stochastic kinetic energy backscatter scheme.",&
128 default=.false.)
129 call get_param(param_file, mdl, "SKEB_NPASS", cs%skeb_npass, &
130 "number of passes of a 9-point smoother of the "//&
131 "dissipation estimate.", default=3, do_not_log=.not.cs%do_skeb)
132 call get_param(param_file, mdl, "SKEB_TAPER_WIDTH", taper_width, &
133 "number of cells over which the stochastic velocity increment "//&
134 "is tapered to zero.", default=4, do_not_log=.not.cs%do_skeb)
135 call get_param(param_file, mdl, "SKEB_USE_GM", cs%skeb_use_gm, &
136 "If true, adds GM work rate to the SKEBS amplitude.", &
137 default=.false., do_not_log=.not.cs%do_skeb)
138 if ((.not. cs%do_skeb) .and. (cs%skeb_use_gm)) call mom_error(fatal, "If SKEB_USE_GM is True "//&
139 "then DO_SKEB must also be True.")
140 call get_param(param_file, mdl, "SKEB_GM_COEF", cs%skeb_gm_coef, &
141 "Fraction of GM work that is added to backscatter rate.", &
142 units="nondim", default=0.0, do_not_log=.not.cs%skeb_use_gm)
143 call get_param(param_file, mdl, "SKEB_USE_FRICT", cs%skeb_use_frict, &
144 "If true, adds horizontal friction dissipation rate "//&
145 "to the SKEBS amplitude.", default=.false., do_not_log=.not.cs%do_skeb)
146 if ((.not. cs%do_skeb) .and. (cs%skeb_use_frict)) call mom_error(fatal, "If SKEB_USE_FRICT is "//&
147 "True then DO_SKEB must also be True.")
148 call get_param(param_file, mdl, "SKEB_FRICT_COEF", cs%skeb_frict_coef, &
149 "Fraction of horizontal friction work that is added to backscatter rate.", &
150 units="nondim", default=0.0, do_not_log=.not.cs%skeb_use_frict)
151 call get_param(param_file, mdl, "PERT_EPBL", cs%pert_epbl, &
152 "If true, then stochastically perturb the kinetic energy "//&
153 "production and dissipation terms. Amplitude and correlations are "//&
154 "controlled by the nam_stoch namelist in the UFS model only.", &
155 default=.false.)
156
157 if (cs%do_sppt .OR. cs%pert_epbl .OR. cs%do_skeb) then
158 num_procs = num_pes()
159 allocate(pelist(num_procs))
160 call get_pelist(pelist,commid = mom_comm)
161 pe_zero = root_pe()
162 nxt = grid%ied - grid%isd + 1
163 nyt = grid%jed - grid%jsd + 1
164 nxb = grid%iedB - grid%isdB + 1
165 nyb = grid%jedB - grid%jsdB + 1
166 call init_stochastic_physics_ocn(dt, grid%geoLonT, grid%geoLatT, nxt, nyt, gv%ke, &
167 grid%geoLonBu, grid%geoLatBu, nxb, nyb, &
168 cs%pert_epbl, cs%do_sppt, cs%do_skeb, pe_zero, mom_comm, iret)
169 if (iret/=0) then
170 call mom_error(fatal, "call to init_stochastic_physics_ocn failed")
171 return
172 endif
173
174 if (cs%do_sppt) allocate(cs%sppt_wts(grid%isd:grid%ied,grid%jsd:grid%jed))
175 if (cs%do_skeb) allocate(cs%skeb_wts(grid%isdB:grid%iedB,grid%jsdB:grid%jedB))
176 if (cs%do_skeb) allocate(cs%skeb_diss(grid%isd:grid%ied,grid%jsd:grid%jed,gv%ke), source=0.)
177 if (cs%pert_epbl) then
178 allocate(cs%epbl1_wts(grid%isd:grid%ied,grid%jsd:grid%jed))
179 allocate(cs%epbl2_wts(grid%isd:grid%ied,grid%jsd:grid%jed))
180 endif
181 endif
182
183 cs%id_sppt_wts = register_diag_field('ocean_model', 'sppt_pattern', cs%diag%axesT1, time, &
184 'random pattern for sppt', 'None')
185 cs%id_skeb_wts = register_diag_field('ocean_model', 'skeb_pattern', cs%diag%axesB1, time, &
186 'random pattern for skeb', 'None')
187 cs%id_epbl1_wts = register_diag_field('ocean_model', 'epbl1_wts', cs%diag%axesT1, time, &
188 'random pattern for KE generation', 'None')
189 cs%id_epbl2_wts = register_diag_field('ocean_model', 'epbl2_wts', cs%diag%axesT1, time, &
190 'random pattern for KE dissipation', 'None')
191 cs%id_skebu = register_diag_field('ocean_model', 'skebu', cs%diag%axesCuL, time, &
192 'zonal current perts', 'None')
193 cs%id_skebv = register_diag_field('ocean_model', 'skebv', cs%diag%axesCvL, time, &
194 'zonal current perts', 'None')
195 cs%id_diss = register_diag_field('ocean_model', 'skeb_amp', cs%diag%axesTL, time, &
196 'SKEB amplitude', 'm s-1')
197 cs%id_psi = register_diag_field('ocean_model', 'psi', cs%diag%axesBL, time, &
198 'stream function', 'None')
199 cs%id_skeb_taperu = register_static_field('ocean_model', 'skeb_taper_u', cs%diag%axesCu1, &
200 'SKEB taper u', 'None', interp_method='none')
201 cs%id_skeb_taperv = register_static_field('ocean_model', 'skeb_taper_v', cs%diag%axesCv1, &
202 'SKEB taper v', 'None', interp_method='none')
203
204 ! Initialize the "taper" fields. These fields multiply the components of the stochastic
205 ! velocity increment in such a way as to smoothly taper them to zero at land boundaries.
206 if ((cs%do_skeb) .or. (cs%id_skeb_taperu > 0) .or. (cs%id_skeb_taperv > 0)) then
207 allocate(cs%taperCu(grid%IsdB:grid%IedB,grid%jsd:grid%jed))
208 allocate(cs%taperCv(grid%isd:grid%ied,grid%JsdB:grid%JedB))
209 ! Initialize taper from land mask
210 do j=grid%jsd,grid%jed ; do i=grid%isdB,grid%iedB
211 cs%taperCu(i,j) = grid%mask2dCu(i,j)
212 enddo ; enddo
213 do j=grid%jsdB,grid%jedB ; do i=grid%isd,grid%ied
214 cs%taperCv(i,j) = grid%mask2dCv(i,j)
215 enddo ; enddo
216 ! Extend taper land
217 do k=1,(taper_width / 2)
218 do j=grid%jsc-1,grid%jec+1 ; do i=grid%iscB-1,grid%iecB+1
219 tmp(i,j) = minval(cs%taperCu(i-1:i+1,j-1:j+1))
220 enddo ; enddo
221 do j=grid%jsc,grid%jec ; do i=grid%iscB,grid%iecB
222 cs%taperCu(i,j) = minval(tmp(i-1:i+1,j-1:j+1))
223 enddo ; enddo
224 do j=grid%jscB-1,grid%jecB+1 ; do i=grid%isc-1,grid%iec+1
225 tmp(i,j) = minval(cs%taperCv(i-1:i+1,j-1:j+1))
226 enddo ; enddo
227 do j=grid%jscB,grid%jecB ; do i=grid%isc,grid%iec
228 cs%taperCv(i,j) = minval(tmp(i-1:i+1,j-1:j+1))
229 enddo ; enddo
230 ! Update halo
231 call pass_vector(cs%taperCu, cs%taperCv, grid%Domain, scalar_pair)
232 enddo
233 ! Smooth tapers. Each call smooths twice.
234 do k=1,(taper_width - (taper_width/2))
235 call smooth_x9_uv(grid, cs%taperCu, cs%taperCv, zero_land=.true.)
236 call pass_vector(cs%taperCu, cs%taperCv, grid%Domain, scalar_pair)
237 enddo
238 endif
239
240 !call uvchksum("SKEB taper [uv]", CS%taperCu, CS%taperCv, grid%HI)
241
242 if (cs%id_skeb_taperu > 0) call post_data(cs%id_skeb_taperu, cs%taperCu, cs%diag, .true.)
243 if (cs%id_skeb_taperv > 0) call post_data(cs%id_skeb_taperv, cs%taperCv, cs%diag, .true.)
244
245 if (cs%do_sppt .OR. cs%pert_epbl .OR. cs%do_skeb) &
246 call mom_mesg(' === COMPLETED MOM STOCHASTIC INITIALIZATION =====')
247
248 call calltree_leave("stochastic_init(), MOM_stochastics.F90")
249
250end subroutine stochastics_init
251
252!> update_ocean_model uses the forcing in Ice_ocean_boundary to advance the
253!! ocean model's state from the input value of Ocean_state (which must be for
254!! time time_start_update) for a time interval of Ocean_coupling_time_step,
255!! returning the publicly visible ocean surface properties in Ocean_sfc and
256!! storing the new ocean properties in Ocean_state.
257subroutine update_stochastics(CS)
258 type(stochastic_cs), intent(inout) :: cs !< diabatic control structure
259 call calltree_enter("update_stochastics(), MOM_stochastics.F90")
260
261! update stochastic physics patterns before running next time-step
262 call run_stochastic_physics_ocn(cs%sppt_wts,cs%skeb_wts,cs%epbl1_wts,cs%epbl2_wts)
263
264 call calltree_leave("update_stochastics(), MOM_stochastics.F90")
265
266end subroutine update_stochastics
267
268subroutine apply_skeb(grid,GV,CS,uc,vc,thickness,tv,dt,Time_end)
269
270 type(ocean_grid_type), intent(in) :: grid !< ocean grid structure
271 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid
272 type(stochastic_cs), intent(inout) :: cs !< stochastic control structure
273
274 real, dimension(SZIB_(grid),SZJ_(grid),SZK_(GV)), intent(inout) :: uc !< zonal velocity [L T-1 ~> m s-1]
275 real, dimension(SZI_(grid),SZJB_(grid),SZK_(GV)), intent(inout) :: vc !< meridional velocity [L T-1 ~> m s-1]
276 real, dimension(SZI_(grid),SZJ_(grid),SZK_(GV)), intent(in) :: thickness !< thickness [H ~> m or kg m-2]
277 type(thermo_var_ptrs), intent(in) :: tv !< points to thermodynamic fields
278 real, intent(in) :: dt !< time increment [T ~> s]
279 type(time_type), intent(in) :: time_end !< Time at the end of the interval
280! locals
281
282 real, dimension(SZIB_(grid),SZJB_(grid),SZK_(GV)) :: psi !< Streamfunction for stochastic velocity increments
283 !! [L2 T-1 ~> m2 s-1]
284 real, dimension(SZIB_(grid),SZJ_(grid) ,SZK_(GV)) :: ustar !< Stochastic u velocity increment [L T-1 ~> m s-1]
285 real, dimension(SZI_(grid) ,SZJB_(grid),SZK_(GV)) :: vstar !< Stochastic v velocity increment [L T-1 ~> m s-1]
286 real, dimension(SZI_(grid),SZJ_(grid)) :: diss_tmp !< Temporary array used in smoothing skeb_diss
287 !! [L2 T-3 ~> m2 s-2]
288 real, dimension(3,3) :: local_weights !< 3x3 stencil weights used in smoothing skeb_diss
289 !! [L2 ~> m2]
290
291 real :: shr,ten,tot,kh
292 integer :: i,j,k,iter
293 integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
294
295 call calltree_enter("apply_skeb(), MOM_stochastics.F90")
296
297 if ((.not. cs%skeb_use_gm) .and. (.not. cs%skeb_use_frict)) then
298 ! fill in halos with zeros
299 do k=1,gv%ke
300 do j=grid%jsd,grid%jed ; do i=grid%isd,grid%ied
301 cs%skeb_diss(i,j,k) = 0.0
302 enddo ; enddo
303 enddo
304
305 !kh needs to be scaled
306
307 kh=1!(120*111)**2
308 do k=1,gv%ke
309 do j=grid%jsc,grid%jec ; do i=grid%isc,grid%iec
310 ! Shear
311 shr = (vc(i,j,k)-vc(i-1,j,k))*grid%mask2dCv(i,j)*grid%mask2dCv(i-1,j)*grid%IdxCv(i,j)+&
312 (uc(i,j,k)-uc(i,j-1,k))*grid%mask2dCu(i,j)*grid%mask2dCu(i,j-1)*grid%IdyCu(i,j)
313 ! Tension
314 ten = (vc(i,j,k)-vc(i-1,j,k))*grid%mask2dCv(i,j)*grid%mask2dCv(i-1,j)*grid%IdyCv(i,j)+&
315 (uc(i,j,k)-uc(i,j-1,k))*grid%mask2dCu(i,j)*grid%mask2dCu(i,j-1)*grid%IdxCu(i,j)
316
317 tot = sqrt( shr**2 + ten**2 ) * grid%mask2dT(i,j)
318 cs%skeb_diss(i,j,k) = tot**3 * kh * grid%areaT(i,j)!!**2
319 enddo ; enddo
320 enddo
321 endif ! Sets CS%skeb_diss without GM or FrictWork
322
323 ! smooth dissipation skeb_npass times
324 do iter=1,cs%skeb_npass
325 if (mod(iter,2) == 1) call pass_var(cs%skeb_diss, grid%domain)
326 do k=1,gv%ke
327 do j=grid%jsc-1,grid%jec+1 ; do i=grid%isc-1,grid%iec+1
328 ! This does not preserve rotational symmetry
329 local_weights = grid%mask2dT(i-1:i+1,j-1:j+1)*grid%areaT(i-1:i+1,j-1:j+1)
330 diss_tmp(i,j) = sum(local_weights*cs%skeb_diss(i-1:i+1,j-1:j+1,k)) / &
331 (sum(local_weights) + 1.e-16)
332 enddo ; enddo
333 do j=grid%jsc-1,grid%jec+1 ; do i=grid%isc-1,grid%iec+1
334 if (grid%mask2dT(i,j)==0.) cycle
335 cs%skeb_diss(i,j,k) = diss_tmp(i,j)
336 enddo ; enddo
337 enddo
338 enddo
339 call pass_var(cs%skeb_diss, grid%domain)
340
341 ! call hchksum(CS%skeb_diss, "SKEB DISS", grid%HI, haloshift=2)
342 ! call qchksum(CS%skeb_wts, "SKEB WTS", grid%HI, haloshift=1)
343
344 do k=1,gv%ke
345 do j=grid%jscB-1,grid%jecB ; do i=grid%iscB-1,grid%iecB
346 psi(i,j,k) = sqrt(0.25 * dt * max((cs%skeb_diss(i ,j ,k) + cs%skeb_diss(i+1,j+1,k)) + &
347 (cs%skeb_diss(i ,j+1,k) + cs%skeb_diss(i+1,j ,k)), 0.) ) &
348 * cs%skeb_wts(i,j)
349 enddo ; enddo
350 enddo
351 !call qchksum(psi,"SKEB PSI", grid%HI, haloshift=1)
352 !call pass_var(psi, grid%domain, position=CORNER)
353 do k=1,gv%ke
354 do j=grid%jsc,grid%jec ; do i=grid%iscB,grid%iecB
355 ustar(i,j,k) = - (psi(i,j,k) - psi(i,j-1,k)) * cs%taperCu(i,j) * grid%IdyCu(i,j)
356 uc(i,j,k) = uc(i,j,k) + ustar(i,j,k)
357 enddo ; enddo
358 do j=grid%jscB,grid%jecB ; do i=grid%isc,grid%iec
359 vstar(i,j,k) = (psi(i,j,k) - psi(i-1,j,k)) * cs%taperCv(i,j) * grid%IdxCv(i,j)
360 vc(i,j,k) = vc(i,j,k) + vstar(i,j,k)
361 enddo ; enddo
362 enddo
363
364 !call uvchksum("SKEB increment [uv]", ustar, vstar, grid%HI)
365
366 call enable_averages(dt, time_end, cs%diag)
367 if (cs%id_diss > 0) then
368 call post_data(cs%id_diss, sqrt(dt * max(cs%skeb_diss(:,:,:), 0.)), cs%diag)
369 endif
370 if (cs%id_skeb_wts > 0) then
371 call post_data(cs%id_skeb_wts, cs%skeb_wts, cs%diag)
372 endif
373 if (cs%id_skebu > 0) then
374 call post_data(cs%id_skebu, ustar(:,:,:), cs%diag)
375 endif
376 if (cs%id_skebv > 0) then
377 call post_data(cs%id_skebv, vstar(:,:,:), cs%diag)
378 endif
379 if (cs%id_psi > 0) then
380 call post_data(cs%id_psi, psi(:,:,:), cs%diag)
381 endif
382 call disable_averaging(cs%diag)
383 cs%skeb_diss(:,:,:) = 0.0 ! Must zero before next time step.
384
385 call calltree_leave("apply_skeb(), MOM_stochastics.F90")
386
387end subroutine apply_skeb
388
389!> Apply a 9-point smoothing filter twice to a pair of velocity components to reduce
390!! horizontal two-grid-point noise.
391!! Note that this subroutine does not conserve angular momentum, so don't use it
392!! in situations where you need conservation. Also note that it assumes that the
393!! input fields have valid values in the first two halo points upon entry.
394subroutine smooth_x9_uv(G, field_u, field_v, zero_land)
395 type(ocean_grid_type), intent(in) :: G !< Ocean grid
396 real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: field_u !< u-point field to be smoothed[arbitrary]
397 real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: field_v !< v-point field to be smoothed [arbitrary]
398 logical, optional, intent(in) :: zero_land !< If present and false, return the average
399 !! of the surrounding ocean points when
400 !! smoothing, otherwise use a value of 0 for
401 !! land points and include them in the averages.
402
403 ! Local variables.
404 real :: fu_prev(SZIB_(G),SZJ_(G)) ! The value of the u-point field at the previous iteration [arbitrary]
405 real :: fv_prev(SZI_(G),SZJB_(G)) ! The value of the v-point field at the previous iteration [arbitrary]
406 real :: Iwts ! The inverse of the sum of the weights [nondim]
407 logical :: zero_land_val ! The value of the zero_land optional argument or .true. if it is absent.
408 integer :: i, j, s, is, ie, js, je, Isq, Ieq, Jsq, Jeq
409
410 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
411 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
412
413 zero_land_val = .true. ; if (present(zero_land)) zero_land_val = zero_land
414
415 do s=1,0,-1
416 fu_prev(:,:) = field_u(:,:)
417 ! apply smoothing on field_u using rotationally symmetric expressions.
418 do j=js-s,je+s ; do i=isq-s,ieq+s ; if (g%mask2dCu(i,j) > 0.0) then
419 iwts = 0.0625
420 if (.not. zero_land_val) &
421 iwts = 1.0 / ( (4.0*g%mask2dCu(i,j) + &
422 ( 2.0*((g%mask2dCu(i-1,j) + g%mask2dCu(i+1,j)) + &
423 (g%mask2dCu(i,j-1) + g%mask2dCu(i,j+1))) + &
424 ((g%mask2dCu(i-1,j-1) + g%mask2dCu(i+1,j+1)) + &
425 (g%mask2dCu(i-1,j+1) + g%mask2dCu(i+1,j-1))) ) ) + 1.0e-16 )
426 field_u(i,j) = iwts * ( 4.0*g%mask2dCu(i,j) * fu_prev(i,j) &
427 + (2.0*((g%mask2dCu(i-1,j) * fu_prev(i-1,j) + g%mask2dCu(i+1,j) * fu_prev(i+1,j)) + &
428 (g%mask2dCu(i,j-1) * fu_prev(i,j-1) + g%mask2dCu(i,j+1) * fu_prev(i,j+1))) &
429 + ((g%mask2dCu(i-1,j-1) * fu_prev(i-1,j-1) + g%mask2dCu(i+1,j+1) * fu_prev(i+1,j+1)) + &
430 (g%mask2dCu(i-1,j+1) * fu_prev(i-1,j+1) + g%mask2dCu(i+1,j-1) * fu_prev(i-1,j-1))) ))
431 endif ; enddo ; enddo
432
433 fv_prev(:,:) = field_v(:,:)
434 ! apply smoothing on field_v using rotationally symmetric expressions.
435 do j=jsq-s,jeq+s ; do i=is-s,ie+s ; if (g%mask2dCv(i,j) > 0.0) then
436 iwts = 0.0625
437 if (.not. zero_land_val) &
438 iwts = 1.0 / ( (4.0*g%mask2dCv(i,j) + &
439 ( 2.0*((g%mask2dCv(i-1,j) + g%mask2dCv(i+1,j)) + &
440 (g%mask2dCv(i,j-1) + g%mask2dCv(i,j+1))) + &
441 ((g%mask2dCv(i-1,j-1) + g%mask2dCv(i+1,j+1)) + &
442 (g%mask2dCv(i-1,j+1) + g%mask2dCv(i+1,j-1))) ) ) + 1.0e-16 )
443 field_v(i,j) = iwts * ( 4.0*g%mask2dCv(i,j) * fv_prev(i,j) &
444 + (2.0*((g%mask2dCv(i-1,j) * fv_prev(i-1,j) + g%mask2dCv(i+1,j) * fv_prev(i+1,j)) + &
445 (g%mask2dCv(i,j-1) * fv_prev(i,j-1) + g%mask2dCv(i,j+1) * fv_prev(i,j+1))) &
446 + ((g%mask2dCv(i-1,j-1) * fv_prev(i-1,j-1) + g%mask2dCv(i+1,j+1) * fv_prev(i+1,j+1)) + &
447 (g%mask2dCv(i-1,j+1) * fv_prev(i-1,j+1) + g%mask2dCv(i+1,j-1) * fv_prev(i-1,j-1))) ))
448 endif ; enddo ; enddo
449 enddo
450
451end subroutine smooth_x9_uv
452
453end module mom_stochastics
454