← back to index

src/parameterizations/stochastic/MOM_stochastics.F90

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!> Top-level module for the MOM6 ocean model in coupled mode.
6module MOM_stochastics
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
15use MOM_diag_mediator, only : register_static_field, enable_averages, disable_averaging
16use MOM_grid, only : ocean_grid_type
17use MOM_variables, only : thermo_var_ptrs
18use MOM_domains, only : pass_var, pass_vector, CORNER, SCALAR_PAIR
19use MOM_verticalGrid, only : verticalGrid_type
20use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe
21use MOM_error_handler, only : callTree_enter, callTree_leave
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
28use stochastic_physics, only : init_stochastic_physics_ocn, run_stochastic_physics_ocn
29
30#include <MOM_memory.h>
31
32implicit none ; private
33
34public stochastics_init, update_stochastics, apply_skeb
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.
801subroutine 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
901 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
982 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
1061 call callTree_enter("stochastic_init(), MOM_stochastics.F90")
1071 if (associated(CS)) then
108 call MOM_error(WARNING, "MOM_stochastics_init called with an "// &
1090 "associated control structure.")
1100 return
1111 else ; allocate(CS) ; endif
112
1131 CS%Time => Time
1141 CS%diag => diag
115
116 ! Read all relevant parameters and write them to the model log.
1171 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.", &
1241 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.",&
1281 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 "//&
1311 "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 "//&
1341 "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.", &
1371 default=.false., do_not_log=.not.CS%do_skeb)
1381 if ((.not. CS%do_skeb) .and. (CS%skeb_use_gm)) call MOM_error(FATAL, "If SKEB_USE_GM is True "//&
1390 "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.", &
1421 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 "//&
1451 "to the SKEBS amplitude.", default=.false., do_not_log=.not.CS%do_skeb)
1461 if ((.not. CS%do_skeb) .and. (CS%skeb_use_frict)) call MOM_error(FATAL, "If SKEB_USE_FRICT is "//&
1470 "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.", &
1501 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.", &
1551 default=.false.)
156
1571 if (CS%do_sppt .OR. CS%pert_epbl .OR. CS%do_skeb) then
1580 num_procs = num_PEs()
1590 allocate(pelist(num_procs))
1600 call Get_PElist(pelist,commID = mom_comm)
1610 pe_zero = root_PE()
1620 nxT = grid%ied - grid%isd + 1
1630 nyT = grid%jed - grid%jsd + 1
1640 nxB = grid%iedB - grid%isdB + 1
1650 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, &
1680 CS%pert_epbl, CS%do_sppt, CS%do_skeb, pe_zero, mom_comm, iret)
1690 if (iret/=0) then
1700 call MOM_error(FATAL, "call to init_stochastic_physics_ocn failed")
1710 return
172 endif
173
1740 if (CS%do_sppt) allocate(CS%sppt_wts(grid%isd:grid%ied,grid%jsd:grid%jed))
1750 if (CS%do_skeb) allocate(CS%skeb_wts(grid%isdB:grid%iedB,grid%jsdB:grid%jedB))
1760 if (CS%do_skeb) allocate(CS%skeb_diss(grid%isd:grid%ied,grid%jsd:grid%jed,GV%ke), source=0.)
1770 if (CS%pert_epbl) then
1780 allocate(CS%epbl1_wts(grid%isd:grid%ied,grid%jsd:grid%jed))
1790 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, &
1841 'random pattern for sppt', 'None')
185 CS%id_skeb_wts = register_diag_field('ocean_model', 'skeb_pattern', CS%diag%axesB1, Time, &
1861 'random pattern for skeb', 'None')
187 CS%id_epbl1_wts = register_diag_field('ocean_model', 'epbl1_wts', CS%diag%axesT1, Time, &
1881 'random pattern for KE generation', 'None')
189 CS%id_epbl2_wts = register_diag_field('ocean_model', 'epbl2_wts', CS%diag%axesT1, Time, &
1901 'random pattern for KE dissipation', 'None')
191 CS%id_skebu = register_diag_field('ocean_model', 'skebu', CS%diag%axesCuL, Time, &
1921 'zonal current perts', 'None')
193 CS%id_skebv = register_diag_field('ocean_model', 'skebv', CS%diag%axesCvL, Time, &
1941 'zonal current perts', 'None')
195 CS%id_diss = register_diag_field('ocean_model', 'skeb_amp', CS%diag%axesTL, Time, &
1961 'SKEB amplitude', 'm s-1')
197 CS%id_psi = register_diag_field('ocean_model', 'psi', CS%diag%axesBL, Time, &
1981 'stream function', 'None')
199 CS%id_skeb_taperu = register_static_field('ocean_model', 'skeb_taper_u', CS%diag%axesCu1, &
2001 'SKEB taper u', 'None', interp_method='none')
201 CS%id_skeb_taperv = register_static_field('ocean_model', 'skeb_taper_v', CS%diag%axesCv1, &
2021 '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.
2061 if ((CS%do_skeb) .or. (CS%id_skeb_taperu > 0) .or. (CS%id_skeb_taperv > 0)) then
2070 allocate(CS%taperCu(grid%IsdB:grid%IedB,grid%jsd:grid%jed))
2080 allocate(CS%taperCv(grid%isd:grid%ied,grid%JsdB:grid%JedB))
209 ! Initialize taper from land mask
2100 do j=grid%jsd,grid%jed ; do I=grid%isdB,grid%iedB
2110 CS%taperCu(I,j) = grid%mask2dCu(I,j)
212 enddo ; enddo
2130 do J=grid%jsdB,grid%jedB ; do i=grid%isd,grid%ied
2140 CS%taperCv(i,J) = grid%mask2dCv(i,J)
215 enddo ; enddo
216 ! Extend taper land
2170 do k=1,(taper_width / 2)
2180 do j=grid%jsc-1,grid%jec+1 ; do I=grid%iscB-1,grid%iecB+1
2190 tmp(I,j) = minval(CS%taperCu(I-1:I+1,j-1:j+1))
220 enddo ; enddo
2210 do j=grid%jsc,grid%jec ; do I=grid%iscB,grid%iecB
2220 CS%taperCu(I,j) = minval(tmp(I-1:I+1,j-1:j+1))
223 enddo ; enddo
2240 do J=grid%jscB-1,grid%jecB+1 ; do i=grid%isc-1,grid%iec+1
2250 tmp(i,J) = minval(CS%taperCv(i-1:i+1,J-1:J+1))
226 enddo ; enddo
2270 do J=grid%jscB,grid%jecB ; do i=grid%isc,grid%iec
2280 CS%taperCv(i,J) = minval(tmp(i-1:i+1,J-1:J+1))
229 enddo ; enddo
230 ! Update halo
2310 call pass_vector(CS%taperCu, CS%taperCv, grid%Domain, SCALAR_PAIR)
232 enddo
233 ! Smooth tapers. Each call smooths twice.
2340 do k=1,(taper_width - (taper_width/2))
2350 call smooth_x9_uv(grid, CS%taperCu, CS%taperCv, zero_land=.true.)
2360 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
2421 if (CS%id_skeb_taperu > 0) call post_data(CS%id_skeb_taperu, CS%taperCu, CS%diag, .true.)
2431 if (CS%id_skeb_taperv > 0) call post_data(CS%id_skeb_taperv, CS%taperCv, CS%diag, .true.)
244
2451 if (CS%do_sppt .OR. CS%pert_epbl .OR. CS%do_skeb) &
2460 call MOM_mesg(' === COMPLETED MOM STOCHASTIC INITIALIZATION =====')
247
2481 call callTree_leave("stochastic_init(), MOM_stochastics.F90")
249
2501end 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.
2570subroutine update_stochastics(CS)
258 type(stochastic_CS), intent(inout) :: CS !< diabatic control structure
2590 call callTree_enter("update_stochastics(), MOM_stochastics.F90")
260
261! update stochastic physics patterns before running next time-step
2620 call run_stochastic_physics_ocn(CS%sppt_wts,CS%skeb_wts,CS%epbl1_wts,CS%epbl2_wts)
263
2640 call callTree_leave("update_stochastics(), MOM_stochastics.F90")
265
2660end subroutine update_stochastics
267
2680subroutine 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
2820 real, dimension(SZIB_(grid),SZJB_(grid),SZK_(GV)) :: psi !< Streamfunction for stochastic velocity increments
283 !! [L2 T-1 ~> m2 s-1]
2840 real, dimension(SZIB_(grid),SZJ_(grid) ,SZK_(GV)) :: ustar !< Stochastic u velocity increment [L T-1 ~> m s-1]
2850 real, dimension(SZI_(grid) ,SZJB_(grid),SZK_(GV)) :: vstar !< Stochastic v velocity increment [L T-1 ~> m s-1]
2860 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
2950 call callTree_enter("apply_skeb(), MOM_stochastics.F90")
296
2970 if ((.not. CS%skeb_use_gm) .and. (.not. CS%skeb_use_frict)) then
298 ! fill in halos with zeros
2990 do k=1,GV%ke
3000 do j=grid%jsd,grid%jed ; do i=grid%isd,grid%ied
3010 CS%skeb_diss(i,j,k) = 0.0
302 enddo ; enddo
303 enddo
304
305 !kh needs to be scaled
306
3070 kh=1!(120*111)**2
3080 do k=1,GV%ke
3090 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)+&
3120 (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)+&
3150 (uc(I,j,k)-uc(I,j-1,k))*grid%mask2dCu(I,j)*grid%mask2dCu(I,j-1)*grid%IdxCu(I,j)
316
3170 tot = sqrt( shr**2 + ten**2 ) * grid%mask2dT(i,j)
3180 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
3240 do iter=1,CS%skeb_npass
3250 if (mod(iter,2) == 1) call pass_var(CS%skeb_diss, grid%domain)
3260 do k=1,GV%ke
3270 do j=grid%jsc-1,grid%jec+1 ; do i=grid%isc-1,grid%iec+1
328 ! This does not preserve rotational symmetry
3290 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)) / &
3310 (sum(local_weights) + 1.E-16)
332 enddo ; enddo
3330 do j=grid%jsc-1,grid%jec+1 ; do i=grid%isc-1,grid%iec+1
3340 if (grid%mask2dT(i,j)==0.) cycle
3350 CS%skeb_diss(i,j,k) = diss_tmp(i,j)
336 enddo ; enddo
337 enddo
338 enddo
3390 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
3440 do k=1,GV%ke
3450 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.) ) &
3480 * 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)
3530 do k=1,GV%ke
3540 do j=grid%jsc,grid%jec ; do I=grid%iscB,grid%iecB
3550 ustar(I,j,k) = - (psi(I,J,k) - psi(I,J-1,k)) * CS%taperCu(I,j) * grid%IdyCu(I,j)
3560 uc(I,j,k) = uc(I,j,k) + ustar(I,j,k)
357 enddo ; enddo
3580 do J=grid%jscB,grid%jecB ; do i=grid%isc,grid%iec
3590 vstar(i,J,k) = (psi(I,J,k) - psi(I-1,J,k)) * CS%taperCv(i,J) * grid%IdxCv(i,J)
3600 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
3660 call enable_averages(dt, Time_end, CS%diag)
3670 if (CS%id_diss > 0) then
3680 call post_data(CS%id_diss, sqrt(dt * max(CS%skeb_diss(:,:,:), 0.)), CS%diag)
369 endif
3700 if (CS%id_skeb_wts > 0) then
3710 call post_data(CS%id_skeb_wts, CS%skeb_wts, CS%diag)
372 endif
3730 if (CS%id_skebu > 0) then
3740 call post_data(CS%id_skebu, ustar(:,:,:), CS%diag)
375 endif
3760 if (CS%id_skebv > 0) then
3770 call post_data(CS%id_skebv, vstar(:,:,:), CS%diag)
378 endif
3790 if (CS%id_psi > 0) then
3800 call post_data(CS%id_psi, psi(:,:,:), CS%diag)
381 endif
3820 call disable_averaging(CS%diag)
3830 CS%skeb_diss(:,:,:) = 0.0 ! Must zero before next time step.
384
3850 call callTree_leave("apply_skeb(), MOM_stochastics.F90")
386
3870end 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.
3940subroutine 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.
4040 real :: fu_prev(SZIB_(G),SZJ_(G)) ! The value of the u-point field at the previous iteration [arbitrary]
4050 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
4100 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
4110 Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
412
4130 zero_land_val = .true. ; if (present(zero_land)) zero_land_val = zero_land
414
4150 do s=1,0,-1
4160 fu_prev(:,:) = field_u(:,:)
417 ! apply smoothing on field_u using rotationally symmetric expressions.
4180 do j=js-s,je+s ; do I=Isq-s,Ieq+s ; if (G%mask2dCu(I,j) > 0.0) then
4190 Iwts = 0.0625
4200 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)) + &
4250 (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)) + &
4300 (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
4330 fv_prev(:,:) = field_v(:,:)
434 ! apply smoothing on field_v using rotationally symmetric expressions.
4350 do J=Jsq-s,Jeq+s ; do i=is-s,ie+s ; if (G%mask2dCv(i,J) > 0.0) then
4360 Iwts = 0.0625
4370 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)) + &
4420 (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)) + &
4470 (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
4510end subroutine smooth_x9_uv
452
4530end module MOM_stochastics
454