MOM_Zanna_Bolton.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!> Calculates Zanna and Bolton 2020 parameterization
6!! Implemented by Perezhogin P.A. Contact: pperezhogin@gmail.com
8
11use mom_diag_mediator, only : diag_ctrl, time_type
12use mom_file_parser, only : get_param, log_version, param_file_type
14use mom_diag_mediator, only : post_data, register_diag_field
15use mom_domains, only : create_group_pass, do_group_pass, group_pass_type, &
16 start_group_pass, complete_group_pass
17use mom_domains, only : to_north, to_east
18use mom_domains, only : pass_var, corner
19use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
20use mom_cpu_clock, only : clock_module, clock_routine
22
23implicit none ; private
24
25#include <MOM_memory.h>
26
28
29!> Control structure for Zanna-Bolton-2020 parameterization.
30type, public :: zb2020_cs ; private
31 ! Parameters
32 real :: amplitude !< The nondimensional scaling factor in ZB model,
33 !! typically 0.1 - 10 [nondim].
34 integer :: zb_type !< Select how to compute the trace part of ZB model:
35 !! 0 - both deviatoric and trace components are computed
36 !! 1 - only deviatoric component is computed
37 !! 2 - only trace component is computed
38 integer :: zb_cons !< Select a discretization scheme for ZB model
39 !! 0 - non-conservative scheme
40 !! 1 - conservative scheme for deviatoric component
41 integer :: hpf_iter !< Number of sharpening passes for the Velocity Gradient (VG) components
42 !! in ZB model.
43 integer :: stress_iter !< Number of smoothing passes for the Stress tensor components
44 !! in ZB model.
45 real :: klower_r_diss !< Attenuation of
46 !! the ZB parameterization in the regions of
47 !! geostrophically-unbalanced flows (Klower 2018, Juricke2020,2019)
48 !! Subgrid stress is multiplied by 1/(1+(shear/(f*R_diss)))
49 !! R_diss=-1: attenuation is not used; typical value R_diss=1.0 [nondim]
50 integer :: klower_shear !< Type of expression for shear in Klower formula
51 !! 0: sqrt(sh_xx**2 + sh_xy**2)
52 !! 1: sqrt(sh_xx**2 + sh_xy**2 + vort_xy**2)
53 integer :: marching_halo !< The number of filter iterations per a single MPI
54 !! exchange
55
56 real, dimension(:,:,:), allocatable :: &
57 sh_xx, & !< Horizontal tension (du/dx - dv/dy) in h (CENTER)
58 !! points including metric terms [T-1 ~> s-1]
59 sh_xy, & !< Horizontal shearing strain (du/dy + dv/dx) in q (CORNER)
60 !! points including metric terms [T-1 ~> s-1]
61 vort_xy, & !< Vertical vorticity (dv/dx - du/dy) in q (CORNER)
62 !! points including metric terms [T-1 ~> s-1]
63 hq !< Thickness in CORNER points [H ~> m or kg m-2]
64
65 real, dimension(:,:,:), allocatable :: &
66 txx, & !< Subgrid stress xx component in h [L2 T-2 ~> m2 s-2]
67 tyy, & !< Subgrid stress yy component in h [L2 T-2 ~> m2 s-2]
68 txy !< Subgrid stress xy component in q [L2 T-2 ~> m2 s-2]
69
70 real, dimension(:,:), allocatable :: &
71 kappa_h, & !< Scaling coefficient in h points [L2 ~> m2]
72 kappa_q !< Scaling coefficient in q points [L2 ~> m2]
73
74 real, allocatable :: &
75 icoriolis_h(:,:), & !< Inverse Coriolis parameter at h points [T ~> s]
76 c_diss(:,:,:) !< Attenuation parameter at h points
77 !! (Klower 2018, Juricke2019,2020) [nondim]
78
79 real, dimension(:,:), allocatable :: &
80 maskw_h, & !< Mask of land point at h points multiplied by filter weight [nondim]
81 maskw_q !< Same mask but for q points [nondim]
82
83 logical :: use_ann !< If True, momentum fluxes are inferred with ANN
84 integer :: stencil_size !< Default is 3x3
85 type(ann_cs) :: ann_tall !< ANN instance for off-diagonal and diagonal stress
86 character(len=200) :: ann_file_tall !< Path to netcdf file with ANN
87 real :: subroundoff_shear !< Small dimensional constant for save division by zero [T-1 ~> s-1]
88
89 type(diag_ctrl), pointer :: diag => null() !< A type that regulates diagnostics output
90 !>@{ Diagnostic handles
91 integer :: id_zb2020u = -1, id_zb2020v = -1, id_ke_zb2020 = -1
92 integer :: id_txx = -1
93 integer :: id_tyy = -1
94 integer :: id_txy = -1
95 integer :: id_cdiss = -1
96 !>@}
97
98 !>@{ CPU time clock IDs
99 integer :: id_clock_module
100 integer :: id_clock_copy
101 integer :: id_clock_cdiss
102 integer :: id_clock_stress
103 integer :: id_clock_stress_ann
104 integer :: id_clock_divergence
105 integer :: id_clock_mpi
106 integer :: id_clock_filter
107 integer :: id_clock_post
108 integer :: id_clock_source
109 !>@}
110
111 !>@{ MPI group passes
112 type(group_pass_type) :: &
113 pass_tq, pass_th, & !< handles for halo passes of Txy and Txx, Tyy
114 pass_xx, pass_xy !< handles for halo passes of sh_xx and sh_xy, vort_xy
115 integer :: stress_halo = -1, & !< The halo size in filter of the stress tensor
116 hpf_halo = -1 !< The halo size in filter of the velocity gradient
117 !>@}
118
119end type zb2020_cs
120
121contains
122
123!> Read parameters, allocate and precompute arrays,
124!! register diagnosicts used in Zanna_Bolton_2020().
125subroutine zb2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020)
126 type(time_type), intent(in) :: time !< The current model time.
127 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
128 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
129 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
130 type(param_file_type), intent(in) :: param_file !< Parameter file parser structure.
131 type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics structure.
132 type(zb2020_cs), intent(inout) :: cs !< ZB2020 control structure.
133 logical, intent(out) :: use_zb2020 !< If true, turns on ZB scheme.
134
135 real :: subroundoff_cor ! A negligible parameter which avoids division by zero
136 ! but small compared to Coriolis parameter [T-1 ~> s-1]
137
138 integer :: is, ie, js, je, isq, ieq, jsq, jeq
139 integer :: i, j
140
141 ! This include declares and sets the variable "version".
142#include "version_variable.h"
143 character(len=40) :: mdl = "MOM_Zanna_Bolton" ! This module's name.
144
145 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
146 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
147
148 call log_version(param_file, mdl, version, "")
149
150 call get_param(param_file, mdl, "USE_ZB2020", use_zb2020, &
151 "If true, turns on Zanna-Bolton-2020 (ZB) " //&
152 "subgrid momentum parameterization of mesoscale eddies.", default=.false.)
153 if (.not. use_zb2020) return
154
155 call get_param(param_file, mdl, "ZB2020_USE_ANN", cs%use_ann, &
156 "ANN inference of momentum fluxes", default=.false.)
157
158 call get_param(param_file, mdl, "ZB2020_ANN_STENCIL_SIZE", cs%stencil_size, &
159 "ANN stencil size", default=3)
160
161 call get_param(param_file, mdl, "ZB2020_ANN_FILE_TALL", cs%ann_file_Tall, &
162 "ANN parameters for prediction of Txy, Txx and Tyy netcdf input", &
163 default="INPUT/EXP1/Tall.nc")
164
165 call get_param(param_file, mdl, "ZB_SCALING", cs%amplitude, &
166 "The nondimensional scaling factor in ZB model, " //&
167 "typically 0.5-2.5", units="nondim", default=0.5)
168
169 call get_param(param_file, mdl, "ZB_TRACE_MODE", cs%ZB_type, &
170 "Select how to compute the trace part of ZB model:\n" //&
171 "\t 0 - both deviatoric and trace components are computed\n" //&
172 "\t 1 - only deviatoric component is computed\n" //&
173 "\t 2 - only trace component is computed", default=0)
174
175 call get_param(param_file, mdl, "ZB_SCHEME", cs%ZB_cons, &
176 "Select a discretization scheme for ZB model:\n" //&
177 "\t 0 - non-conservative scheme\n" //&
178 "\t 1 - conservative scheme for deviatoric component", default=1)
179
180 call get_param(param_file, mdl, "VG_SHARP_PASS", cs%HPF_iter, &
181 "Number of sharpening passes for the Velocity Gradient (VG) components " //&
182 "in ZB model.", default=0)
183
184 call get_param(param_file, mdl, "STRESS_SMOOTH_PASS", cs%Stress_iter, &
185 "Number of smoothing passes for the Stress tensor components " //&
186 "in ZB model.", default=0)
187
188 call get_param(param_file, mdl, "ZB_KLOWER_R_DISS", cs%Klower_R_diss, &
189 "Attenuation of " //&
190 "the ZB parameterization in the regions of " //&
191 "geostrophically-unbalanced flows (Klower 2018, Juricke2020,2019). " //&
192 "Subgrid stress is multiplied by 1/(1+(shear/(f*R_diss))):\n" //&
193 "\t R_diss=-1. - attenuation is not used\n\t R_diss= 1. - typical value", &
194 units="nondim", default=-1.)
195
196 call get_param(param_file, mdl, "ZB_KLOWER_SHEAR", cs%Klower_shear, &
197 "Type of expression for shear in Klower formula:\n" //&
198 "\t 0: sqrt(sh_xx**2 + sh_xy**2)\n" //&
199 "\t 1: sqrt(sh_xx**2 + sh_xy**2 + vort_xy**2)", &
200 default=1, do_not_log=.not.cs%Klower_R_diss>0)
201
202 call get_param(param_file, mdl, "ZB_MARCHING_HALO", cs%Marching_halo, &
203 "The number of filter iterations per single MPI " //&
204 "exchange", default=4, do_not_log=(cs%Stress_iter==0).and.(cs%HPF_iter==0))
205
206 ! Register fields for output from this module.
207 cs%diag => diag
208
209 cs%id_ZB2020u = register_diag_field('ocean_model', 'ZB2020u', diag%axesCuL, time, &
210 'Zonal Acceleration from Zanna-Bolton 2020', 'm s-2', conversion=us%L_T2_to_m_s2)
211 cs%id_ZB2020v = register_diag_field('ocean_model', 'ZB2020v', diag%axesCvL, time, &
212 'Meridional Acceleration from Zanna-Bolton 2020', 'm s-2', conversion=us%L_T2_to_m_s2)
213 cs%id_KE_ZB2020 = register_diag_field('ocean_model', 'KE_ZB2020', diag%axesTL, time, &
214 'Kinetic Energy Source from Horizontal Viscosity', &
215 'm3 s-3', conversion=gv%H_to_m*(us%L_T_to_m_s**2)*us%s_to_T)
216
217 cs%id_Txx = register_diag_field('ocean_model', 'Txx', diag%axesTL, time, &
218 'Diagonal term (Txx) in the ZB stress tensor', 'm2 s-2', conversion=us%L_T_to_m_s**2)
219
220 cs%id_Tyy = register_diag_field('ocean_model', 'Tyy', diag%axesTL, time, &
221 'Diagonal term (Tyy) in the ZB stress tensor', 'm2 s-2', conversion=us%L_T_to_m_s**2)
222
223 cs%id_Txy = register_diag_field('ocean_model', 'Txy', diag%axesBL, time, &
224 'Off-diagonal term (Txy) in the ZB stress tensor', 'm2 s-2', conversion=us%L_T_to_m_s**2)
225
226 if (cs%Klower_R_diss > 0) then
227 cs%id_cdiss = register_diag_field('ocean_model', 'c_diss', diag%axesTL, time, &
228 'Klower (2018) attenuation coefficient', 'nondim')
229 endif
230
231 ! Clock IDs
232 ! Only module is measured with syncronization. While smaller
233 ! parts are measured without - because these are nested clocks.
234 cs%id_clock_module = cpu_clock_id('(Ocean Zanna-Bolton-2020)', grain=clock_module)
235 cs%id_clock_copy = cpu_clock_id('(ZB2020 copy fields)', grain=clock_routine, sync=.false.)
236 cs%id_clock_cdiss = cpu_clock_id('(ZB2020 compute c_diss)', grain=clock_routine, sync=.false.)
237 cs%id_clock_stress = cpu_clock_id('(ZB2020 compute stress)', grain=clock_routine, sync=.false.)
238 cs%id_clock_stress_ANN = cpu_clock_id('(ZB2020 compute stress ANN)', grain=clock_routine, sync=.false.)
239 cs%id_clock_divergence = cpu_clock_id('(ZB2020 compute divergence)', grain=clock_routine, sync=.false.)
240 cs%id_clock_mpi = cpu_clock_id('(ZB2020 filter MPI exchanges)', grain=clock_routine, sync=.false.)
241 cs%id_clock_filter = cpu_clock_id('(ZB2020 filter no MPI)', grain=clock_routine, sync=.false.)
242 cs%id_clock_post = cpu_clock_id('(ZB2020 post data)', grain=clock_routine, sync=.false.)
243 cs%id_clock_source = cpu_clock_id('(ZB2020 compute energy source)', grain=clock_routine, sync=.false.)
244
245 cs%subroundoff_shear = 1e-30 * us%T_to_s
246 if (cs%use_ann) then
247 call ann_init(cs%ann_Tall, cs%ann_file_Tall)
248 endif
249
250 ! Allocate memory
251 ! We set the stress tensor and velocity gradient tensor to zero
252 ! with full halo because they potentially may be filtered
253 ! with marching halo algorithm
254 allocate(cs%sh_xx(szi_(g),szj_(g),szk_(gv)), source=0.)
255 allocate(cs%sh_xy(szib_(g),szjb_(g),szk_(gv)), source=0.)
256 allocate(cs%vort_xy(szib_(g),szjb_(g),szk_(gv)), source=0.)
257 allocate(cs%hq(szib_(g),szjb_(g),szk_(gv)))
258
259 allocate(cs%Txx(szi_(g),szj_(g),szk_(gv)), source=0.)
260 allocate(cs%Tyy(szi_(g),szj_(g),szk_(gv)), source=0.)
261 allocate(cs%Txy(szib_(g),szjb_(g),szk_(gv)), source=0.)
262 allocate(cs%kappa_h(szi_(g),szj_(g)))
263 allocate(cs%kappa_q(szib_(g),szjb_(g)))
264
265 ! Precomputing the scaling coefficient
266 ! Mask is included to automatically satisfy B.C.
267 do j=js-2,je+2 ; do i=is-2,ie+2
268 cs%kappa_h(i,j) = -cs%amplitude * g%areaT(i,j) * g%mask2dT(i,j)
269 enddo ; enddo
270
271 do j=jsq-2,jeq+2 ; do i=isq-2,ieq+2
272 cs%kappa_q(i,j) = -cs%amplitude * g%areaBu(i,j) * g%mask2dBu(i,j)
273 enddo ; enddo
274
275 if (cs%Klower_R_diss > 0) then
276 allocate(cs%ICoriolis_h(szi_(g),szj_(g)))
277 allocate(cs%c_diss(szi_(g),szj_(g),szk_(gv)))
278
279 subroundoff_cor = 1e-30 * us%T_to_s
280 ! Precomputing 1/(f * R_diss)
281 do j=js-1,je+1 ; do i=is-1,ie+1
282 cs%ICoriolis_h(i,j) = 1. / ((abs(0.25 * ((g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j-1)) &
283 + (g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j-1)))) + subroundoff_cor) &
284 * cs%Klower_R_diss)
285 enddo ; enddo
286 endif
287
288 if (cs%Stress_iter > 0 .or. cs%HPF_iter > 0) then
289 ! Include 1/16. factor to the mask for filter implementation
290 allocate(cs%maskw_h(szi_(g),szj_(g))) ; cs%maskw_h(:,:) = g%mask2dT(:,:) * 0.0625
291 allocate(cs%maskw_q(szib_(g),szjb_(g))) ; cs%maskw_q(:,:) = g%mask2dBu(:,:) * 0.0625
292 endif
293
294 ! Initialize MPI group passes
295 if (cs%Stress_iter > 0) then
296 ! reduce size of halo exchange accordingly to
297 ! Marching halo, number of iterations and the array size
298 ! But let exchange width be at least 1
299 cs%Stress_halo = max(min(cs%Marching_halo, cs%Stress_iter, &
300 g%Domain%nihalo, g%Domain%njhalo), 1)
301
302 call create_group_pass(cs%pass_Tq, cs%Txy, g%Domain, halo=cs%Stress_halo, &
303 position=corner)
304 call create_group_pass(cs%pass_Th, cs%Txx, g%Domain, halo=cs%Stress_halo)
305 call create_group_pass(cs%pass_Th, cs%Tyy, g%Domain, halo=cs%Stress_halo)
306 endif
307
308 if (cs%HPF_iter > 0) then
309 ! The minimum halo size is 2 because it is requirement for the
310 ! outputs of function filter_velocity_gradients
311 cs%HPF_halo = max(min(cs%Marching_halo, cs%HPF_iter, &
312 g%Domain%nihalo, g%Domain%njhalo), 2)
313
314 call create_group_pass(cs%pass_xx, cs%sh_xx, g%Domain, halo=cs%HPF_halo)
315 call create_group_pass(cs%pass_xy, cs%sh_xy, g%Domain, halo=cs%HPF_halo, &
316 position=corner)
317 call create_group_pass(cs%pass_xy, cs%vort_xy, g%Domain, halo=cs%HPF_halo, &
318 position=corner)
319 endif
320
321end subroutine zb2020_init
322
323!> Deallocate any variables allocated in ZB_2020_init
324subroutine zb2020_end(CS)
325 type(zb2020_cs), intent(inout) :: cs !< ZB2020 control structure.
326
327 deallocate(cs%sh_xx)
328 deallocate(cs%sh_xy)
329 deallocate(cs%vort_xy)
330 deallocate(cs%hq)
331
332 deallocate(cs%Txx)
333 deallocate(cs%Tyy)
334 deallocate(cs%Txy)
335 deallocate(cs%kappa_h)
336 deallocate(cs%kappa_q)
337
338 if (cs%Klower_R_diss > 0) then
339 deallocate(cs%ICoriolis_h)
340 deallocate(cs%c_diss)
341 endif
342
343 if (cs%Stress_iter > 0 .or. cs%HPF_iter > 0) then
344 deallocate(cs%maskw_h)
345 deallocate(cs%maskw_q)
346 endif
347
348 if (cs%use_ann) then
349 call ann_end(cs%ann_Tall)
350 endif
351
352end subroutine zb2020_end
353
354!> Save precomputed velocity gradients and thickness
355!! from the horizontal eddy viscosity module
356!! We save as much halo for velocity gradients as possible
357!! In symmetric (preferable) memory model: halo 2 for sh_xx
358!! and halo 1 for sh_xy and vort_xy
359!! We apply zero boundary conditions to velocity gradients
360!! which is required for filtering operations
361subroutine zb2020_copy_gradient_and_thickness(sh_xx, sh_xy, vort_xy, hq, &
362 G, GV, CS, k)
363 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
364 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
365 type(zb2020_cs), intent(inout) :: cs !< ZB2020 control structure.
366
367 real, dimension(SZIB_(G),SZJB_(G)), &
368 intent(in) :: sh_xy !< horizontal shearing strain (du/dy + dv/dx)
369 !! including metric terms [T-1 ~> s-1]
370 real, dimension(SZIB_(G),SZJB_(G)), &
371 intent(in) :: vort_xy !< Vertical vorticity (dv/dx - du/dy)
372 !! including metric terms [T-1 ~> s-1]
373 real, dimension(SZIB_(G),SZJB_(G)), &
374 intent(in) :: hq !< harmonic mean of the harmonic means
375 !! of the u- & v point thicknesses [H ~> m or kg m-2]
376
377 real, dimension(SZI_(G),SZJ_(G)), &
378 intent(in) :: sh_xx !< horizontal tension (du/dx - dv/dy)
379 !! including metric terms [T-1 ~> s-1]
380
381 integer, intent(in) :: k !< The vertical index of the layer to be passed.
382
383 integer :: is, ie, js, je, isq, ieq, jsq, jeq
384 integer :: i, j
385
386 call cpu_clock_begin(cs%id_clock_copy)
387
388 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
389 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
390
391 do j=js-1,jeq ; do i=is-1,ieq
392 cs%hq(i,j,k) = hq(i,j)
393 enddo ; enddo
394
395 ! No physical B.C. is required for
396 ! sh_xx in ZB2020. However, filtering
397 ! may require BC
398 do j=jsq-1,je+2 ; do i=isq-1,ie+2
399 cs%sh_xx(i,j,k) = sh_xx(i,j) * g%mask2dT(i,j)
400 enddo ; enddo
401
402 ! We multiply by mask to remove
403 ! implicit dependence on CS%no_slip
404 ! flag in hor_visc module
405 do j=js-2,jeq+1 ; do i=is-2,ieq+1
406 cs%sh_xy(i,j,k) = sh_xy(i,j) * g%mask2dBu(i,j)
407 enddo ; enddo
408
409 do j=js-2,jeq+1 ; do i=is-2,ieq+1
410 cs%vort_xy(i,j,k) = vort_xy(i,j) * g%mask2dBu(i,j)
411 enddo ; enddo
412
413 call cpu_clock_end(cs%id_clock_copy)
414
416
417!> Baroclinic Zanna-Bolton-2020 parameterization, see
418!! eq. 6 in https://laurezanna.github.io/files/Zanna-Bolton-2020.pdf
419!! We compute the lateral stress tensor according to ZB2020 model
420!! and update the acceleration due to eddy viscosity (diffu, diffv)
421!! as follows:
422!! diffu = diffu + ZB2020u
423!! diffv = diffv + ZB2020v
424subroutine zb2020_lateral_stress(u, v, h, diffu, diffv, G, GV, CS, &
425 dx2h, dy2h, dx2q, dy2q)
426 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
427 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
428 type(zb2020_cs), intent(inout) :: cs !< ZB2020 control structure.
429
430 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
431 intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
432 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
433 intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
434 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
435 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
436
437 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
438 intent(inout) :: diffu !< Zonal acceleration due to eddy viscosity.
439 !! It is updated with ZB closure [L T-2 ~> m s-2]
440 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
441 intent(inout) :: diffv !< Meridional acceleration due to eddy viscosity.
442 !! It is updated with ZB closure [L T-2 ~> m s-2]
443
444 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: dx2h !< dx^2 at h points [L2 ~> m2]
445 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: dy2h !< dy^2 at h points [L2 ~> m2]
446
447 real, dimension(SZIB_(G),SZJB_(G)), intent(in) :: dx2q !< dx^2 at q points [L2 ~> m2]
448 real, dimension(SZIB_(G),SZJB_(G)), intent(in) :: dy2q !< dy^2 at q points [L2 ~> m2]
449
450 call cpu_clock_begin(cs%id_clock_module)
451
452 ! Compute attenuation if specified
453 call compute_c_diss(g, gv, cs)
454
455 ! Sharpen velocity gradients if specified
456 call filter_velocity_gradients(g, gv, cs)
457
458 ! Compute the stress tensor given the
459 ! (optionally sharpened) velocity gradients
460 if (cs%use_ann) then
461 call compute_stress_ann_collocated(g, gv, cs)
462 else
463 call compute_stress(g, gv, cs)
464 endif
465
466 ! Smooth the stress tensor if specified
467 call filter_stress(g, gv, cs)
468
469 ! Update the acceleration due to eddy viscosity (diffu, diffv)
470 ! with the ZB2020 lateral parameterization
471 call compute_stress_divergence(u, v, h, diffu, diffv, &
472 dx2h, dy2h, dx2q, dy2q, &
473 g, gv, cs)
474
475 call cpu_clock_begin(cs%id_clock_post)
476 if (cs%id_Txx>0) call post_data(cs%id_Txx, cs%Txx, cs%diag)
477 if (cs%id_Tyy>0) call post_data(cs%id_Tyy, cs%Tyy, cs%diag)
478 if (cs%id_Txy>0) call post_data(cs%id_Txy, cs%Txy, cs%diag)
479
480 if (cs%id_cdiss>0) call post_data(cs%id_cdiss, cs%c_diss, cs%diag)
481 call cpu_clock_end(cs%id_clock_post)
482
483 call cpu_clock_end(cs%id_clock_module)
484
485end subroutine zb2020_lateral_stress
486
487!> Compute the attenuation parameter similarly
488!! to Klower2018, Juricke2019,2020: c_diss = 1/(1+(shear/(f*R_diss)))
489!! where shear = sqrt(sh_xx**2 + sh_xy**2) or shear = sqrt(sh_xx**2 + sh_xy**2 + vort_xy**2)
490!! In symmetric memory model, components of velocity gradient tensor
491!! should have halo 1 and zero boundary conditions. The result: c_diss having halo 1.
492subroutine compute_c_diss(G, GV, CS)
493 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
494 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
495 type(zb2020_cs), intent(inout) :: CS !< ZB2020 control structure.
496
497 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
498 integer :: i, j, k
499
500 real :: shear ! Shear in Klower2018 formula at h points [T-1 ~> s-1]
501
502 if (.not. cs%Klower_R_diss > 0) &
503 return
504
505 call cpu_clock_begin(cs%id_clock_cdiss)
506
507 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
508 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
509
510 do k=1,nz
511
512 ! sqrt(sh_xx**2 + sh_xy**2)
513 if (cs%Klower_shear == 0) then
514 do j=js-1,je+1 ; do i=is-1,ie+1
515 shear = sqrt(cs%sh_xx(i,j,k)**2 + 0.25 * ( &
516 ((cs%sh_xy(i-1,j-1,k)**2) + (cs%sh_xy(i,j ,k)**2)) &
517 + ((cs%sh_xy(i-1,j ,k)**2) + (cs%sh_xy(i,j-1,k)**2)) &
518 ))
519 cs%c_diss(i,j,k) = 1. / (1. + shear * cs%ICoriolis_h(i,j))
520 enddo ; enddo
521
522 ! sqrt(sh_xx**2 + sh_xy**2 + vort_xy**2)
523 elseif (cs%Klower_shear == 1) then
524 do j=js-1,je+1 ; do i=is-1,ie+1
525 shear = sqrt(cs%sh_xx(i,j,k)**2 + 0.25 * ( &
526 ((cs%sh_xy(i-1,j-1,k)**2 + cs%vort_xy(i-1,j-1,k)**2) &
527 + (cs%sh_xy(i,j,k)**2 + cs%vort_xy(i,j,k)**2)) &
528 + ((cs%sh_xy(i-1,j,k)**2 + cs%vort_xy(i-1,j,k)**2) &
529 + (cs%sh_xy(i,j-1,k)**2 + cs%vort_xy(i,j-1,k)**2)) &
530 ))
531 cs%c_diss(i,j,k) = 1. / (1. + shear * cs%ICoriolis_h(i,j))
532 enddo ; enddo
533 endif
534
535 enddo ! end of k loop
536
537 call cpu_clock_end(cs%id_clock_cdiss)
538
539end subroutine compute_c_diss
540
541!> Compute stress tensor T =
542!! (Txx, Txy;
543!! Txy, Tyy)
544!! Which consists of the deviatoric and trace components, respectively:
545!! T = (-vort_xy * sh_xy, vort_xy * sh_xx;
546!! vort_xy * sh_xx, vort_xy * sh_xy) +
547!! 1/2 * (vort_xy^2 + sh_xy^2 + sh_xx^2, 0;
548!! 0, vort_xy^2 + sh_xy^2 + sh_xx^2)
549!! This stress tensor is multiplied by precomputed kappa=-CS%amplitude * G%area:
550!! T -> T * kappa
551!! The sign of the stress tensor is such that (neglecting h):
552!! (du/dt, dv/dt) = div(T)
553!! In symmetric memory model: sh_xy and vort_xy should have halo 1
554!! and zero B.C.; sh_xx should have halo 2 and zero B.C.
555!! Result: Txx, Tyy, Txy with halo 1 and zero B.C.
556subroutine compute_stress(G, GV, CS)
557 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
558 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
559 type(zb2020_cs), intent(inout) :: CS !< ZB2020 control structure.
560
561 real :: &
562 vort_xy_h, & ! Vorticity interpolated to h point [T-1 ~> s-1]
563 sh_xy_h ! Shearing strain interpolated to h point [T-1 ~> s-1]
564
565 real :: &
566 sh_xx_q ! Horizontal tension interpolated to q point [T-1 ~> s-1]
567
568 ! Local variables
569 real :: sum_sq ! 1/2*(vort_xy^2 + sh_xy^2 + sh_xx^2) in h point [T-2 ~> s-2]
570 real :: vort_sh ! vort_xy*sh_xy in h point [T-2 ~> s-2]
571
572 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
573 integer :: i, j, k
574
575 logical :: sum_sq_flag ! Flag to compute trace
576 logical :: vort_sh_scheme_0, vort_sh_scheme_1 ! Flags to compute diagonal trace-free part
577
578 call cpu_clock_begin(cs%id_clock_stress)
579
580 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
581 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
582
583 sum_sq = 0.
584 vort_sh = 0.
585
586 sum_sq_flag = cs%ZB_type /= 1
587 vort_sh_scheme_0 = cs%ZB_type /= 2 .and. cs%ZB_cons == 0
588 vort_sh_scheme_1 = cs%ZB_type /= 2 .and. cs%ZB_cons == 1
589
590 do k=1,nz
591
592 ! compute Txx, Tyy tensor
593 do j=js-1,je+1 ; do i=is-1,ie+1
594 ! It is assumed that B.C. is applied to sh_xy and vort_xy
595 sh_xy_h = 0.25 * ( (cs%sh_xy(i-1,j-1,k) + cs%sh_xy(i,j,k)) &
596 + (cs%sh_xy(i-1,j,k) + cs%sh_xy(i,j-1,k)) )
597
598 vort_xy_h = 0.25 * ( (cs%vort_xy(i-1,j-1,k) + cs%vort_xy(i,j,k)) &
599 + (cs%vort_xy(i-1,j,k) + cs%vort_xy(i,j-1,k)) )
600
601 if (sum_sq_flag) then
602 sum_sq = 0.5 * &
603 ((vort_xy_h * vort_xy_h &
604 + sh_xy_h * sh_xy_h) &
605 + cs%sh_xx(i,j,k) * cs%sh_xx(i,j,k) &
606 )
607 endif
608
609 if (vort_sh_scheme_0) &
610 vort_sh = vort_xy_h * sh_xy_h
611
612 if (vort_sh_scheme_1) then
613 ! It is assumed that B.C. is applied to sh_xy and vort_xy
614 vort_sh = 0.25 * ( &
615 (((g%areaBu(i-1,j-1) * cs%vort_xy(i-1,j-1,k)) * cs%sh_xy(i-1,j-1,k)) + &
616 ((g%areaBu(i ,j ) * cs%vort_xy(i ,j ,k)) * cs%sh_xy(i ,j ,k))) + &
617 (((g%areaBu(i-1,j ) * cs%vort_xy(i-1,j ,k)) * cs%sh_xy(i-1,j ,k)) + &
618 ((g%areaBu(i ,j-1) * cs%vort_xy(i ,j-1,k)) * cs%sh_xy(i ,j-1,k))) &
619 ) * g%IareaT(i,j)
620 endif
621
622 ! B.C. is already applied in kappa_h
623 cs%Txx(i,j,k) = cs%kappa_h(i,j) * (- vort_sh + sum_sq)
624 cs%Tyy(i,j,k) = cs%kappa_h(i,j) * (+ vort_sh + sum_sq)
625
626 enddo ; enddo
627
628 ! Here we assume that Txy is initialized to zero
629 if (cs%ZB_type /= 2) then
630 do j=jsq-1,jeq+1 ; do i=isq-1,ieq+1
631 sh_xx_q = 0.25 * ( (cs%sh_xx(i+1,j+1,k) + cs%sh_xx(i,j,k)) &
632 + (cs%sh_xx(i+1,j,k) + cs%sh_xx(i,j+1,k)))
633 ! B.C. is already applied in kappa_q
634 cs%Txy(i,j,k) = cs%kappa_q(i,j) * (cs%vort_xy(i,j,k) * sh_xx_q)
635
636 enddo ; enddo
637 endif
638
639 enddo ! end of k loop
640
641 call cpu_clock_end(cs%id_clock_stress)
642
643end subroutine compute_stress
644
645!> Compute stress tensor T =
646!! (Txx, Txy;
647!! Txy, Tyy)
648!! with ANN in non-dimensional form:
649!! T = dx^2 * |grad V|^2 * ANN(grad V / |grad V|)
650!! The sign of the stress tensor is such that:
651!! (du/dt, dv/dt) = 1/h * div(h * T)
652!! Algorithm:
653!! 1) Interpolate input features (sh_xy, sh_xx, vort_xy) to grid centers
654!! 2) Compute norm of velocity gradients on a stencil
655!! 3) Non-dimensionalize input features
656!! 4) Make ANN inference in grid centers
657!! 5) Restore physical dimensionality and interpolate Txy back to corners
658subroutine compute_stress_ann_collocated(G, GV, CS)
659 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
660 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
661 type(zb2020_cs), intent(inout) :: CS !< ZB2020 control structure.
662
663 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
664 integer :: i, j, k, m
665 integer :: ii, jj
666 integer :: nij
667
668 real, allocatable :: x(:,:) ! Vector of non-dimensional input features
669 ! number of horizontal grid points x
670 ! (sh_xy, sh_xx, vort_xy) on a stencil [nondim]
671 real, allocatable :: y(:,:) ! Vector of nondimensional
672 ! output features number of horizontal grid points x
673 ! (Txy,Txx,Tyy) [nondim]
674 real :: yy(3) ! Vector of dimensional
675 ! output features (Txy,Txx,Tyy) [L2 T-2 ~> m2 s-2]
676 real :: tmp ! Temporal value of squared norm [T-2 ~> s-2]
677 integer :: offset ! Half the stencil size. Used for selection
678 integer :: stencil_points ! The number of points after flattening
679
680 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
681 sh_xy_h, & ! sh_xy interpolated to the center [T-1 ~> s-1]
682 vort_xy_h, & ! vort_xy interpolated to the center [T-1 ~> s-1]
683 norm_h ! Norm of input feautres in center points [T-1 ~> s-1]
684
685 real, dimension(SZI_(G),SZJ_(G)) :: &
686 sqr_h, & ! Squared norm of velocity gradients in center points [T-2 ~> s-2]
687 Txy ! Predicted Txy in center points [T-1 ~> s-1]
688
689 call cpu_clock_begin(cs%id_clock_stress_ANN)
690
691 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
692 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
693
694 ! Number of horizontal grid points in ANN inference loop below
695 nij = (ie - is + 5) * (je - js + 5)
696 allocate(x(nij, 3 * cs%stencil_size**2))
697 allocate(y(nij, 3))
698
699 sh_xy_h = 0.
700 vort_xy_h = 0.
701 norm_h = 0.
702
703 call pass_var(cs%sh_xy, g%Domain, clock=cs%id_clock_mpi, position=corner)
704 call pass_var(cs%sh_xx, g%Domain, clock=cs%id_clock_mpi)
705 call pass_var(cs%vort_xy, g%Domain, clock=cs%id_clock_mpi, position=corner)
706
707 offset = (cs%stencil_size-1)/2
708 stencil_points = cs%stencil_size**2
709
710 ! Interpolate input features
711 do k=1,nz
712 do j=js-2,je+2 ; do i=is-2,ie+2
713 ! It is assumed that B.C. is applied to sh_xy and vort_xy
714 sh_xy_h(i,j,k) = 0.25 * ( (cs%sh_xy(i-1,j-1,k) + cs%sh_xy(i,j,k)) &
715 + (cs%sh_xy(i-1,j,k) + cs%sh_xy(i,j-1,k)) )
716
717 vort_xy_h(i,j,k) = 0.25 * ( (cs%vort_xy(i-1,j-1,k) + cs%vort_xy(i,j,k)) &
718 + (cs%vort_xy(i-1,j,k) + cs%vort_xy(i,j-1,k)) )
719
720 sqr_h(i,j) = (((cs%sh_xx(i,j,k)**2) + (sh_xy_h(i,j,k)**2)) + (vort_xy_h(i,j,k)**2)) * g%mask2dT(i,j)
721 enddo ; enddo
722
723 do j=js,je ; do i=is,ie
724 tmp = 0.0
725 do jj=j-offset,j+offset ; do ii=i-offset,i+offset
726 tmp = tmp + sqr_h(ii,jj)
727 enddo ; enddo
728 norm_h(i,j,k) = sqrt(tmp)
729 enddo ; enddo
730 enddo
731
732 call pass_var(sh_xy_h, g%Domain, clock=cs%id_clock_mpi)
733 call pass_var(vort_xy_h, g%Domain, clock=cs%id_clock_mpi)
734 call pass_var(norm_h, g%Domain, clock=cs%id_clock_mpi)
735
736 do k=1,nz
737 m = 0
738 do j=js-2,je+2 ; do i=is-2,ie+2
739 m = m + 1
740 x(m,1:stencil_points) = &
741 reshape(sh_xy_h(i-offset:i+offset, &
742 j-offset:j+offset,k), (/stencil_points/))
743 x(m,stencil_points+1:2*stencil_points) = &
744 reshape(cs%sh_xx(i-offset:i+offset, &
745 j-offset:j+offset,k), (/stencil_points/))
746 x(m,2*stencil_points+1:3*stencil_points) = &
747 reshape(vort_xy_h(i-offset:i+offset, &
748 j-offset:j+offset,k), (/stencil_points/))
749
750 x(m,:) = x(m,:) / (norm_h(i,j,k) + cs%subroundoff_shear)
751 enddo ; enddo
752
753 call ann_apply_array_sio(nij, x, y, cs%ann_Tall)
754
755 m = 0
756 do j=js-2,je+2 ; do i=is-2,ie+2
757 m = m+1
758 yy(:) = y(m, :) * norm_h(i,j,k) * norm_h(i,j,k) * cs%kappa_h(i,j)
759
760 txy(i,j) = yy(1)
761 cs%Txx(i,j,k) = yy(2)
762 cs%Tyy(i,j,k) = yy(3)
763 enddo ; enddo
764
765 do j=jsq-1,jeq+1 ; do i=isq-1,ieq+1
766 cs%Txy(i,j,k) = 0.25 * ( (txy(i+1,j+1) + txy(i,j)) &
767 + (txy(i+1,j) + txy(i,j+1))) * g%mask2dBu(i,j)
768 enddo ; enddo
769
770 enddo ! end of k loop
771
772 call pass_var(cs%Txy, g%Domain, clock=cs%id_clock_mpi, position=corner)
773 call pass_var(cs%Txx, g%Domain, clock=cs%id_clock_mpi)
774 call pass_var(cs%Tyy, g%Domain, clock=cs%id_clock_mpi)
775
776 deallocate(x)
777 deallocate(y)
778
779 call cpu_clock_end(cs%id_clock_stress_ANN)
780
782
783!> Compute the divergence of subgrid stress
784!! weighted with thickness, i.e.
785!! (fx,fy) = 1/h Div(h * [Txx, Txy; Txy, Tyy])
786!! and update the acceleration due to eddy viscosity as
787!! diffu = diffu + dx; diffv = diffv + dy
788!! Optionally, before computing the divergence, we attenuate the stress
789!! according to the Klower formula.
790!! In symmetric memory model: Txx, Tyy, Txy, c_diss should have halo 1
791!! with applied zero B.C.
792subroutine compute_stress_divergence(u, v, h, diffu, diffv, dx2h, dy2h, dx2q, dy2q, G, GV, CS)
793 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
794 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
795 type(zb2020_cs), intent(in) :: CS !< ZB2020 control structure.
796 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
797 intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
798 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
799 intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
800 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
801 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
802 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
803 intent(out) :: diffu !< Zonal acceleration due to convergence of
804 !! along-coordinate stress tensor [L T-2 ~> m s-2]
805 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
806 intent(out) :: diffv !< Meridional acceleration due to convergence
807 !! of along-coordinate stress tensor [L T-2 ~> m s-2]
808 real, dimension(SZI_(G),SZJ_(G)), &
809 intent(in) :: dx2h !< dx^2 at h points [L2 ~> m2]
810 real, dimension(SZI_(G),SZJ_(G)), &
811 intent(in) :: dy2h !< dy^2 at h points [L2 ~> m2]
812 real, dimension(SZIB_(G),SZJB_(G)), &
813 intent(in) :: dx2q !< dx^2 at q points [L2 ~> m2]
814 real, dimension(SZIB_(G),SZJB_(G)), &
815 intent(in) :: dy2q !< dy^2 at q points [L2 ~> m2]
816
817 ! Local variables
818 real, dimension(SZI_(G),SZJ_(G)) :: &
819 Mxx, & ! Subgrid stress Txx multiplied by thickness and dy^2 [H L4 T-2 ~> m5 s-2]
820 Myy ! Subgrid stress Tyy multiplied by thickness and dx^2 [H L4 T-2 ~> m5 s-2]
821
822 real, dimension(SZIB_(G),SZJB_(G)) :: &
823 Mxy ! Subgrid stress Txy multiplied by thickness [H L2 T-2 ~> m3 s-2]
824 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &
825 ZB2020u !< Zonal acceleration due to convergence of
826 !! along-coordinate stress tensor for ZB model
827 !! [L T-2 ~> m s-2]
828 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: &
829 ZB2020v !< Meridional acceleration due to convergence
830 !! of along-coordinate stress tensor for ZB model
831 !! [L T-2 ~> m s-2]
832
833 real :: h_u ! Thickness interpolated to u points [H ~> m or kg m-2].
834 real :: h_v ! Thickness interpolated to v points [H ~> m or kg m-2].
835 real :: fx ! Zonal acceleration [L T-2 ~> m s-2]
836 real :: fy ! Meridional acceleration [L T-2 ~> m s-2]
837
838 real :: h_neglect ! Thickness so small it can be lost in
839 ! roundoff and so neglected [H ~> m or kg m-2]
840
841 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
842 integer :: i, j, k
843 logical :: save_ZB2020u, save_ZB2020v ! Save the acceleration due to ZB2020 model
844
845 call cpu_clock_begin(cs%id_clock_divergence)
846
847 save_zb2020u = (cs%id_ZB2020u > 0) .or. (cs%id_KE_ZB2020 > 0)
848 save_zb2020v = (cs%id_ZB2020v > 0) .or. (cs%id_KE_ZB2020 > 0)
849
850 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
851 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
852
853 h_neglect = gv%H_subroundoff
854
855 do k=1,nz
856 if (cs%Klower_R_diss > 0) then
857 do j=js-1,jeq ; do i=is-1,ieq
858 mxy(i,j) = (cs%Txy(i,j,k) * &
859 (0.25 * ( (cs%c_diss(i,j ,k) + cs%c_diss(i+1,j+1,k)) &
860 + (cs%c_diss(i,j+1,k) + cs%c_diss(i+1,j ,k))) &
861 ) &
862 ) * cs%hq(i,j,k)
863 enddo ; enddo
864 else
865 do j=js-1,jeq ; do i=is-1,ieq
866 mxy(i,j) = cs%Txy(i,j,k) * cs%hq(i,j,k)
867 enddo ; enddo
868 endif
869
870 if (cs%Klower_R_diss > 0) then
871 do j=js-1,je+1 ; do i=is-1,ie+1
872 mxx(i,j) = ((cs%Txx(i,j,k) * cs%c_diss(i,j,k)) * h(i,j,k)) * dy2h(i,j)
873 myy(i,j) = ((cs%Tyy(i,j,k) * cs%c_diss(i,j,k)) * h(i,j,k)) * dx2h(i,j)
874 enddo ; enddo
875 else
876 do j=js-1,je+1 ; do i=is-1,ie+1
877 mxx(i,j) = ((cs%Txx(i,j,k)) * h(i,j,k)) * dy2h(i,j)
878 myy(i,j) = ((cs%Tyy(i,j,k)) * h(i,j,k)) * dx2h(i,j)
879 enddo ; enddo
880 endif
881
882 ! Evaluate du/dt=1/h x.Div(h T) (Line 1495 of MOM_hor_visc.F90)
883 do j=js,je ; do i=isq,ieq
884 h_u = 0.5 * (g%mask2dT(i,j)*h(i,j,k) + g%mask2dT(i+1,j)*h(i+1,j,k)) + h_neglect
885 fx = ((g%IdyCu(i,j)*(mxx(i+1,j) - mxx(i,j)) + &
886 g%IdxCu(i,j)*((dx2q(i,j)*mxy(i,j)) - (dx2q(i,j-1)*mxy(i,j-1)))) * &
887 g%IareaCu(i,j)) / h_u
888 diffu(i,j,k) = diffu(i,j,k) + fx
889 if (save_zb2020u) &
890 zb2020u(i,j,k) = fx
891 enddo ; enddo
892
893 ! Evaluate dv/dt=1/h y.Div(h T) (Line 1517 of MOM_hor_visc.F90)
894 do j=jsq,jeq ; do i=is,ie
895 h_v = 0.5 * (g%mask2dT(i,j)*h(i,j,k) + g%mask2dT(i,j+1)*h(i,j+1,k)) + h_neglect
896 fy = ((g%IdxCv(i,j)*(myy(i,j+1) - myy(i,j)) + &
897 g%IdyCv(i,j)*((dy2q(i,j)*mxy(i,j)) - (dy2q(i-1,j)*mxy(i-1,j)))) * &
898 g%IareaCv(i,j)) / h_v
899 diffv(i,j,k) = diffv(i,j,k) + fy
900 if (save_zb2020v) &
901 zb2020v(i,j,k) = fy
902 enddo ; enddo
903
904 enddo ! end of k loop
905
906 call cpu_clock_end(cs%id_clock_divergence)
907
908 call cpu_clock_begin(cs%id_clock_post)
909 if (cs%id_ZB2020u>0) call post_data(cs%id_ZB2020u, zb2020u, cs%diag)
910 if (cs%id_ZB2020v>0) call post_data(cs%id_ZB2020v, zb2020v, cs%diag)
911 call cpu_clock_end(cs%id_clock_post)
912
913 call compute_energy_source(u, v, h, zb2020u, zb2020v, g, gv, cs)
914
915end subroutine compute_stress_divergence
916
917!> Filtering of the velocity gradients sh_xx, sh_xy, vort_xy.
918!! Here instead of smoothing we do sharpening, i.e.
919!! return (initial - smoothed) fields.
920!! The algorithm: marching halo with non-blocking grouped MPI
921!! exchanges. The input array sh_xx should have halo 2 with
922!! applied zero B.C. The arrays sh_xy and vort_xy should have
923!! halo 1 with applied B.C. The output have the same halo and B.C.
924subroutine filter_velocity_gradients(G, GV, CS)
925 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
926 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
927 type(zb2020_cs), intent(inout) :: CS !< ZB2020 control structure.
928
929 real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: &
930 sh_xx ! Copy of CS%sh_xx [T-1 ~> s-1]
931 real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)) :: &
932 sh_xy, vort_xy ! Copy of CS%sh_xy and CS%vort_xy [T-1 ~> s-1]
933
934 integer :: xx_halo, xy_halo, vort_halo ! currently available halo for gradient components
935 integer :: xx_iter, xy_iter, vort_iter ! remaining number of iterations
936 integer :: niter ! required number of iterations
937
938 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
939 integer :: i, j, k
940
941 niter = cs%HPF_iter
942
943 if (niter == 0) return
944
945 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
946 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
947
948 if (.not. g%symmetric) &
949 call do_group_pass(cs%pass_xx, g%Domain, &
950 clock=cs%id_clock_mpi)
951
952 ! This is just copy of the array
953 call cpu_clock_begin(cs%id_clock_filter)
954 do k=1,nz
955 ! Halo of size 2 is valid
956 do j=js-2,je+2 ; do i=is-2,ie+2
957 sh_xx(i,j,k) = cs%sh_xx(i,j,k)
958 enddo ; enddo
959 ! Only halo of size 1 is valid
960 do j=jsq-1,jeq+1 ; do i=isq-1,ieq+1
961 sh_xy(i,j,k) = cs%sh_xy(i,j,k)
962 vort_xy(i,j,k) = cs%vort_xy(i,j,k)
963 enddo ; enddo
964 enddo
965 call cpu_clock_end(cs%id_clock_filter)
966
967 xx_halo = 2 ; xy_halo = 1 ; vort_halo = 1
968 xx_iter = niter ; xy_iter = niter ; vort_iter = niter
969
970 do while &
971 (xx_iter > 0 .or. xy_iter > 0 .or. & ! filter iterations remain to be done
972 xx_halo < 2 .or. xy_halo < 1) ! there is no halo for VG tensor
973
974 ! ---------- filtering sh_xx ---------
975 if (xx_halo < 2) then
976 call complete_group_pass(cs%pass_xx, g%Domain, clock=cs%id_clock_mpi)
977 xx_halo = cs%HPF_halo
978 endif
979
980 call filter_hq(g, gv, cs, xx_halo, xx_iter, h=cs%sh_xx)
981
982 if (xx_halo < 2) &
983 call start_group_pass(cs%pass_xx, g%Domain, clock=cs%id_clock_mpi)
984
985 ! ------ filtering sh_xy, vort_xy ----
986 if (xy_halo < 1) then
987 call complete_group_pass(cs%pass_xy, g%Domain, clock=cs%id_clock_mpi)
988 xy_halo = cs%HPF_halo ; vort_halo = cs%HPF_halo
989 endif
990
991 call filter_hq(g, gv, cs, xy_halo, xy_iter, q=cs%sh_xy)
992 call filter_hq(g, gv, cs, vort_halo, vort_iter, q=cs%vort_xy)
993
994 if (xy_halo < 1) &
995 call start_group_pass(cs%pass_xy, g%Domain, clock=cs%id_clock_mpi)
996
997 enddo
998
999 ! We implement sharpening by computing residual
1000 ! B.C. are already applied to all fields
1001 call cpu_clock_begin(cs%id_clock_filter)
1002 do k=1,nz
1003 do j=js-2,je+2 ; do i=is-2,ie+2
1004 cs%sh_xx(i,j,k) = sh_xx(i,j,k) - cs%sh_xx(i,j,k)
1005 enddo ; enddo
1006 do j=jsq-1,jeq+1 ; do i=isq-1,ieq+1
1007 cs%sh_xy(i,j,k) = sh_xy(i,j,k) - cs%sh_xy(i,j,k)
1008 cs%vort_xy(i,j,k) = vort_xy(i,j,k) - cs%vort_xy(i,j,k)
1009 enddo ; enddo
1010 enddo
1011 call cpu_clock_end(cs%id_clock_filter)
1012
1013 if (.not. g%symmetric) &
1014 call do_group_pass(cs%pass_xy, g%Domain, &
1015 clock=cs%id_clock_mpi)
1016
1017end subroutine filter_velocity_gradients
1018
1019!> Filtering of the stress tensor Txx, Tyy, Txy.
1020!! The algorithm: marching halo with non-blocking grouped MPI
1021!! exchanges. The input arrays (Txx, Tyy, Txy) must have halo 1
1022!! with zero B.C. applied. The output have the same halo and B.C.
1023subroutine filter_stress(G, GV, CS)
1024 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1025 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1026 type(zb2020_cs), intent(inout) :: CS !< ZB2020 control structure.
1027
1028 integer :: Txx_halo, Tyy_halo, Txy_halo ! currently available halo for stress components
1029 integer :: Txx_iter, Tyy_iter, Txy_iter ! remaining number of iterations
1030 integer :: niter ! required number of iterations
1031
1032 niter = cs%Stress_iter
1033
1034 if (niter == 0) return
1035
1036 txx_halo = 1 ; tyy_halo = 1 ; txy_halo = 1 ; ! these are required halo for Txx, Tyy, Txy
1037 txx_iter = niter ; tyy_iter = niter ; txy_iter = niter
1038
1039 do while &
1040 (txx_iter > 0 .or. txy_iter > 0 .or. & ! filter iterations remain to be done
1041 txx_halo < 1 .or. txy_halo < 1) ! there is no halo for Txx or Txy
1042
1043 ! ---------- filtering Txy -----------
1044 if (txy_halo < 1) then
1045 call complete_group_pass(cs%pass_Tq, g%Domain, clock=cs%id_clock_mpi)
1046 txy_halo = cs%Stress_halo
1047 endif
1048
1049 call filter_hq(g, gv, cs, txy_halo, txy_iter, q=cs%Txy)
1050
1051 if (txy_halo < 1) &
1052 call start_group_pass(cs%pass_Tq, g%Domain, clock=cs%id_clock_mpi)
1053
1054 ! ------- filtering Txx, Tyy ---------
1055 if (txx_halo < 1) then
1056 call complete_group_pass(cs%pass_Th, g%Domain, clock=cs%id_clock_mpi)
1057 txx_halo = cs%Stress_halo ; tyy_halo = cs%Stress_halo
1058 endif
1059
1060 call filter_hq(g, gv, cs, txx_halo, txx_iter, h=cs%Txx)
1061 call filter_hq(g, gv, cs, tyy_halo, tyy_iter, h=cs%Tyy)
1062
1063 if (txx_halo < 1) &
1064 call start_group_pass(cs%pass_Th, g%Domain, clock=cs%id_clock_mpi)
1065
1066 enddo
1067
1068end subroutine filter_stress
1069
1070!> Wrapper for filter_3D function. The border indices for q and h
1071!! arrays are substituted.
1072subroutine filter_hq(G, GV, CS, current_halo, remaining_iterations, q, h)
1073 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1074 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1075 type(zb2020_cs), intent(in) :: CS !< ZB2020 control structure.
1076 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, &
1077 intent(inout) :: h !< Input/output array in h points [arbitrary]
1078 real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)), optional, &
1079 intent(inout) :: q !< Input/output array in q points [arbitrary]
1080 integer, intent(inout) :: current_halo !< Currently available halo points
1081 integer, intent(inout) :: remaining_iterations !< The number of iterations to perform
1082
1083 logical :: direction ! The direction of the first 1D filter
1084
1085 direction = (mod(g%first_direction,2) == 0)
1086
1087 call cpu_clock_begin(cs%id_clock_filter)
1088
1089 if (present(h)) then
1090 call filter_3d(h, cs%maskw_h, &
1091 g%isd, g%ied, g%jsd, g%jed, &
1092 g%isc, g%iec, g%jsc, g%jec, gv%ke, &
1093 current_halo, remaining_iterations, &
1094 direction)
1095 endif
1096
1097 if (present(q)) then
1098 call filter_3d(q, cs%maskw_q, &
1099 g%IsdB, g%IedB, g%JsdB, g%JedB, &
1100 g%IscB, g%IecB, g%JscB, g%JecB, gv%ke, &
1101 current_halo, remaining_iterations, &
1102 direction)
1103 endif
1104
1105 call cpu_clock_end(cs%id_clock_filter)
1106end subroutine filter_hq
1107
1108!> Spatial lateral filter applied to 3D array. The lateral filter is given
1109!! by the convolutional kernel:
1110!! [1 2 1]
1111!! C = |2 4 2| * 1/16
1112!! [1 2 1]
1113!! The fast algorithm decomposes the 2D filter into two 1D filters as follows:
1114!! [1]
1115!! C = |2| * [1 2 1] * 1/16
1116!! [1]
1117!! The input array must have zero B.C. applied. B.C. is applied for output array.
1118!! Note that maskw contains both land mask and 1/16 factor.
1119!! Filter implements marching halo. The available halo is specified and as many
1120!! filter iterations as possible and as needed are performed.
1121subroutine filter_3d(x, maskw, isd, ied, jsd, jed, is, ie, js, je, nz, &
1122 current_halo, remaining_iterations, &
1123 direction)
1124 integer, intent(in) :: isd !< Indices of array size
1125 integer, intent(in) :: ied !< Indices of array size
1126 integer, intent(in) :: jsd !< Indices of array size
1127 integer, intent(in) :: jed !< Indices of array size
1128 integer, intent(in) :: is !< Indices of owned points
1129 integer, intent(in) :: ie !< Indices of owned points
1130 integer, intent(in) :: js !< Indices of owned points
1131 integer, intent(in) :: je !< Indices of owned points
1132 integer, intent(in) :: nz !< Vertical array size
1133 real, dimension(isd:ied,jsd:jed,nz), &
1134 intent(inout) :: x !< Input/output array [arbitrary]
1135 real, dimension(isd:ied,jsd:jed), &
1136 intent(in) :: maskw !< Mask array of land points divided by 16 [nondim]
1137 integer, intent(inout) :: current_halo !< Currently available halo points
1138 integer, intent(inout) :: remaining_iterations !< The number of iterations to perform
1139 logical, intent(in) :: direction !< The direction of the first 1D filter
1140
1141 real, parameter :: weight = 2. ! Filter weight [nondim]
1142 integer :: i, j, k, iter, niter, halo
1143
1144 real :: tmp(isd:ied, jsd:jed) ! Array with temporary results [arbitrary]
1145
1146 ! Do as many iterations as needed and possible
1147 niter = min(current_halo, remaining_iterations)
1148 if (niter == 0) return ! nothing to do
1149
1150 ! Update remaining iterations
1151 remaining_iterations = remaining_iterations - niter
1152 ! Update halo information
1153 current_halo = current_halo - niter
1154
1155 do k=1,nz
1156 halo = niter-1 + &
1157 current_halo ! Save as many halo points as possible
1158 do iter=1,niter
1159
1160 if (direction) then
1161 do j = js-halo, je+halo ; do i = is-halo-1, ie+halo+1
1162 tmp(i,j) = weight * x(i,j,k) + (x(i,j-1,k) + x(i,j+1,k))
1163 enddo ; enddo
1164
1165 do j = js-halo, je+halo ; do i = is-halo, ie+halo
1166 x(i,j,k) = (weight * tmp(i,j) + (tmp(i-1,j) + tmp(i+1,j))) * maskw(i,j)
1167 enddo ; enddo
1168 else
1169 do j = js-halo-1, je+halo+1 ; do i = is-halo, ie+halo
1170 tmp(i,j) = weight * x(i,j,k) + (x(i-1,j,k) + x(i+1,j,k))
1171 enddo ; enddo
1172
1173 do j = js-halo, je+halo ; do i = is-halo, ie+halo
1174 x(i,j,k) = (weight * tmp(i,j) + (tmp(i,j-1) + tmp(i,j+1))) * maskw(i,j)
1175 enddo ; enddo
1176 endif
1177
1178 halo = halo - 1
1179 enddo
1180 enddo
1181
1182end subroutine filter_3d
1183
1184!> Computes the 3D energy source term for the ZB2020 scheme
1185!! similarly to MOM_diagnostics.F90, specifically 1125 line.
1186subroutine compute_energy_source(u, v, h, fx, fy, G, GV, CS)
1187 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1188 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1189 type(zb2020_cs), intent(in) :: CS !< ZB2020 control structure.
1190
1191 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1192 intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
1193 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1194 intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
1195 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1196 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
1197
1198 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1199 intent(in) :: fx !< Zonal acceleration due to convergence of
1200 !! along-coordinate stress tensor [L T-2 ~> m s-2]
1201 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1202 intent(in) :: fy !< Meridional acceleration due to convergence
1203 !! of along-coordinate stress tensor [L T-2 ~> m s-2]
1204
1205 real :: KE_term(SZI_(G),SZJ_(G),SZK_(GV)) ! A term in the kinetic energy budget
1206 ! [H L2 T-3 ~> m3 s-3 or W m-2]
1207 real :: KE_u(SZIB_(G),SZJ_(G)) ! The area integral of a KE term in a layer at u-points
1208 ! [H L4 T-3 ~> m5 s-3 or kg m2 s-3]
1209 real :: KE_v(SZI_(G),SZJB_(G)) ! The area integral of a KE term in a layer at v-points
1210 ! [H L4 T-3 ~> m5 s-3 or kg m2 s-3]
1211
1212 real :: uh ! Transport through zonal faces = u*h*dy,
1213 ! [H L2 T-1 ~> m3 s-1 or kg s-1].
1214 real :: vh ! Transport through meridional faces = v*h*dx,
1215 ! [H L2 T-1 ~> m3 s-1 or kg s-1].
1216
1217 type(group_pass_type) :: pass_KE_uv ! A handle used for group halo passes
1218
1219 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1220 integer :: i, j, k
1221
1222 if (cs%id_KE_ZB2020 > 0) then
1223 call cpu_clock_begin(cs%id_clock_source)
1224 call create_group_pass(pass_ke_uv, ke_u, ke_v, g%Domain, to_north+to_east)
1225
1226 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1227 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1228
1229 ke_term(:,:,:) = 0.
1230 ! Calculate the KE source from Zanna-Bolton2020 [H L2 T-3 ~> m3 s-3].
1231 do k=1,nz
1232 ke_u(:,:) = 0.
1233 ke_v(:,:) = 0.
1234 do j=js,je ; do i=isq,ieq
1235 uh = u(i,j,k) * 0.5 * (g%mask2dT(i,j)*h(i,j,k) + g%mask2dT(i+1,j)*h(i+1,j,k)) * &
1236 g%dyCu(i,j)
1237 ke_u(i,j) = uh * g%dxCu(i,j) * fx(i,j,k)
1238 enddo ; enddo
1239 do j=jsq,jeq ; do i=is,ie
1240 vh = v(i,j,k) * 0.5 * (g%mask2dT(i,j)*h(i,j,k) + g%mask2dT(i,j+1)*h(i,j+1,k)) * &
1241 g%dxCv(i,j)
1242 ke_v(i,j) = vh * g%dyCv(i,j) * fy(i,j,k)
1243 enddo ; enddo
1244 call do_group_pass(pass_ke_uv, g%domain, clock=cs%id_clock_mpi)
1245 do j=js,je ; do i=is,ie
1246 ke_term(i,j,k) = 0.5 * g%IareaT(i,j) &
1247 * ((ke_u(i,j) + ke_u(i-1,j)) + (ke_v(i,j) + ke_v(i,j-1)))
1248 enddo ; enddo
1249 enddo
1250
1251 call cpu_clock_end(cs%id_clock_source)
1252
1253 call cpu_clock_begin(cs%id_clock_post)
1254 call post_data(cs%id_KE_ZB2020, ke_term, cs%diag)
1255 call cpu_clock_end(cs%id_clock_post)
1256 endif
1257
1258end subroutine compute_energy_source
1259
1260end module mom_zanna_bolton