MOM6
MOM_vert_friction.F90
1 !> Implements vertical viscosity (vertvisc)
2 module mom_vert_friction
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 use mom_domains, only : pass_var, to_all, omit_corners
6 use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
7 use mom_diag_mediator, only : diag_ctrl
8 use mom_debugging, only : uvchksum, hchksum
9 use mom_error_handler, only : mom_error, fatal, warning, note
10 use mom_file_parser, only : get_param, log_version, param_file_type
11 use mom_forcing_type, only : mech_forcing
12 use mom_get_input, only : directories
13 use mom_grid, only : ocean_grid_type
14 use mom_open_boundary, only : ocean_obc_type, obc_simple, obc_none, obc_direction_e
15 use mom_open_boundary, only : obc_direction_w, obc_direction_n, obc_direction_s
16 use mom_pointaccel, only : write_u_accel, write_v_accel, pointaccel_init
17 use mom_pointaccel, only : pointaccel_cs
18 use mom_time_manager, only : time_type, time_type_to_real, operator(-)
19 use mom_unit_scaling, only : unit_scale_type
20 use mom_variables, only : thermo_var_ptrs, vertvisc_type
21 use mom_variables, only : cont_diag_ptrs, accel_diag_ptrs
22 use mom_variables, only : ocean_internal_state
23 use mom_verticalgrid, only : verticalgrid_type
24 use mom_wave_interface, only : wave_parameters_cs
25 implicit none ; private
26 
27 #include <MOM_memory.h>
28 
29 public vertvisc, vertvisc_remnant, vertvisc_coef
30 public vertvisc_limit_vel, vertvisc_init, vertvisc_end
31 public updatecfltruncationvalue
32 
33 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
34 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
35 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
36 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
37 
38 !> The control structure with parameters and memory for the MOM_vert_friction module
39 type, public :: vertvisc_cs ; private
40  real :: hmix !< The mixed layer thickness in thickness units [H ~> m or kg m-2].
41  real :: hmix_stress !< The mixed layer thickness over which the wind
42  !! stress is applied with direct_stress [H ~> m or kg m-2].
43  real :: kvml !< The mixed layer vertical viscosity [Z2 T-1 ~> m2 s-1].
44  real :: kv !< The interior vertical viscosity [Z2 T-1 ~> m2 s-1].
45  real :: hbbl !< The static bottom boundary layer thickness [H ~> m or kg m-2].
46  real :: kvbbl !< The vertical viscosity in the bottom boundary
47  !! layer [Z2 T-1 ~> m2 s-1].
48 
49  real :: maxvel !< Velocity components greater than maxvel are truncated [L T-1 ~> m s-1].
50  real :: vel_underflow !< Velocity components smaller than vel_underflow
51  !! are set to 0 [L T-1 ~> m s-1].
52  logical :: cfl_based_trunc !< If true, base truncations on CFL numbers, not
53  !! absolute velocities.
54  real :: cfl_trunc !< Velocity components will be truncated when they
55  !! are large enough that the corresponding CFL number
56  !! exceeds this value, nondim.
57  real :: cfl_report !< The value of the CFL number that will cause the
58  !! accelerations to be reported, nondim. CFL_report
59  !! will often equal CFL_trunc.
60  real :: truncramptime !< The time-scale over which to ramp up the value of
61  !! CFL_trunc from CFL_truncS to CFL_truncE
62  real :: cfl_truncs !< The start value of CFL_trunc
63  real :: cfl_trunce !< The end/target value of CFL_trunc
64  logical :: cflrampingisactivated = .false. !< True if the ramping has been initialized
65  type(time_type) :: rampstarttime !< The time at which the ramping of CFL_trunc starts
66 
67  real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NK_INTERFACE_) :: &
68  a_u !< The u-drag coefficient across an interface [Z T-1 ~> m s-1].
69  real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
70  h_u !< The effective layer thickness at u-points [H ~> m or kg m-2].
71  real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NK_INTERFACE_) :: &
72  a_v !< The v-drag coefficient across an interface [Z T-1 ~> m s-1].
73  real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
74  h_v !< The effective layer thickness at v-points [H ~> m or kg m-2].
75  real, pointer, dimension(:,:) :: a1_shelf_u => null() !< The u-momentum coupling coefficient under
76  !! ice shelves [Z T-1 ~> m s-1]. Retained to determine stress under shelves.
77  real, pointer, dimension(:,:) :: a1_shelf_v => null() !< The v-momentum coupling coefficient under
78  !! ice shelves [Z T-1 ~> m s-1]. Retained to determine stress under shelves.
79 
80  logical :: split !< If true, use the split time stepping scheme.
81  logical :: bottomdraglaw !< If true, the bottom stress is calculated with a
82  !! drag law c_drag*|u|*u. The velocity magnitude
83  !! may be an assumed value or it may be based on the
84  !! actual velocity in the bottommost HBBL, depending
85  !! on whether linear_drag is true.
86  logical :: channel_drag !< If true, the drag is exerted directly on each
87  !! layer according to what fraction of the bottom
88  !! they overlie.
89  logical :: harmonic_visc !< If true, the harmonic mean thicknesses are used
90  !! to calculate the viscous coupling between layers
91  !! except near the bottom. Otherwise the arithmetic
92  !! mean thickness is used except near the bottom.
93  real :: harm_bl_val !< A scale to determine when water is in the boundary
94  !! layers based solely on harmonic mean thicknesses
95  !! for the purpose of determining the extent to which
96  !! the thicknesses used in the viscosities are upwinded.
97  logical :: direct_stress !< If true, the wind stress is distributed over the
98  !! topmost Hmix_stress of fluid and KVML may be very small.
99  logical :: dynamic_viscous_ml !< If true, use the results from a dynamic
100  !! calculation, perhaps based on a bulk Richardson
101  !! number criterion, to determine the mixed layer
102  !! thickness for viscosity.
103  logical :: answers_2018 !< If true, use the order of arithmetic and expressions that recover the
104  !! answers from the end of 2018. Otherwise, use expressions that do not
105  !! use an arbitary and hard-coded maximum viscous coupling coefficient
106  !! between layers.
107  logical :: debug !< If true, write verbose checksums for debugging purposes.
108  integer :: nkml !< The number of layers in the mixed layer.
109  integer, pointer :: ntrunc !< The number of times the velocity has been
110  !! truncated since the last call to write_energy.
111  character(len=200) :: u_trunc_file !< The complete path to a file in which a column of
112  !! u-accelerations are written if velocity truncations occur.
113  character(len=200) :: v_trunc_file !< The complete path to a file in which a column of
114  !! v-accelerations are written if velocity truncations occur.
115  logical :: stokesmixing !< If true, do Stokes drift mixing via the Lagrangian current
116  !! (Eulerian plus Stokes drift). False by default and set
117  !! via STOKES_MIXING_COMBINED.
118 
119  type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
120  !! timing of diagnostic output.
121 
122  !>@{ Diagnostic identifiers
123  integer :: id_du_dt_visc = -1, id_dv_dt_visc = -1, id_au_vv = -1, id_av_vv = -1
124  integer :: id_h_u = -1, id_h_v = -1, id_hml_u = -1 , id_hml_v = -1
125  integer :: id_taux_bot = -1, id_tauy_bot = -1
126  integer :: id_kv_slow = -1, id_kv_u = -1, id_kv_v = -1
127  ! integer :: id_hf_du_dt_visc = -1, id_hf_dv_dt_visc = -1
128  integer :: id_hf_du_dt_visc_2d = -1, id_hf_dv_dt_visc_2d = -1
129  !>@}
130 
131  type(pointaccel_cs), pointer :: pointaccel_csp => null() !< A pointer to the control structure
132  !! for recording accelerations leading to velocity truncations
133 
134  ! real, pointer :: hf_du_dt_visc(:,:,:) => NULL() ! Zonal friction accel. x fract. thickness [L T-2 ~> m s-2].
135  ! real, pointer :: hf_dv_dt_visc(:,:,:) => NULL() ! Merdional friction accel. x fract. thickness [L T-2 ~> m s-2].
136  ! 3D diagnostics hf_du(dv)_dt_visc are commented because there is no clarity on proper remapping grid option.
137  ! The code is retained for degugging purposes in the future.
138 
139 end type vertvisc_cs
140 
141 contains
142 
143 !> Perform a fully implicit vertical diffusion
144 !! of momentum. Stress top and bottom boundary conditions are used.
145 !!
146 !! This is solving the tridiagonal system
147 !! \f[ \left(h_k + a_{k + 1/2} + a_{k - 1/2} + r_k\right) u_k^{n+1}
148 !! = h_k u_k^n + a_{k + 1/2} u_{k+1}^{n+1} + a_{k - 1/2} u_{k-1}^{n+1} \f]
149 !! where \f$a_{k + 1/2} = \Delta t \nu_{k + 1/2} / h_{k + 1/2}\f$
150 !! is the <em>interfacial coupling thickness per time step</em>,
151 !! encompassing background viscosity as well as contributions from
152 !! enhanced mixed and bottom layer viscosities.
153 !! $r_k$ is a Rayleight drag term due to channel drag.
154 !! There is an additional stress term on the right-hand side
155 !! if DIRECT_STRESS is true, applied to the surface layer.
156 
157 subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
158  taux_bot, tauy_bot, Waves)
159  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
160  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
161  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
162  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
163  intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]
164  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
165  intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]
166  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
167  intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
168  type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
169  type(vertvisc_type), intent(inout) :: visc !< Viscosities and bottom drag
170  real, intent(in) :: dt !< Time increment [T ~> s]
171  type(ocean_obc_type), pointer :: obc !< Open boundary condition structure
172  type(accel_diag_ptrs), intent(inout) :: adp !< Accelerations in the momentum
173  !! equations for diagnostics
174  type(cont_diag_ptrs), intent(inout) :: cdp !< Continuity equation terms
175  type(vertvisc_cs), pointer :: cs !< Vertical viscosity control structure
176  real, dimension(SZIB_(G),SZJ_(G)), &
177  optional, intent(out) :: taux_bot !< Zonal bottom stress from ocean to
178  !! rock [R L Z T-2 ~> Pa]
179  real, dimension(SZI_(G),SZJB_(G)), &
180  optional, intent(out) :: tauy_bot !< Meridional bottom stress from ocean to
181  !! rock [R L Z T-2 ~> Pa]
182  type(wave_parameters_cs), &
183  optional, pointer :: waves !< Container for wave/Stokes information
184 
185  ! Fields from forces used in this subroutine:
186  ! taux: Zonal wind stress [R L Z T-2 ~> Pa].
187  ! tauy: Meridional wind stress [R L Z T-2 ~> Pa].
188 
189  ! Local variables
190 
191  real :: b1(szib_(g)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
192  real :: c1(szib_(g),szk_(g)) ! A variable used by the tridiagonal solver [nondim].
193  real :: d1(szib_(g)) ! d1=1-c1 is used by the tridiagonal solver [nondim].
194  real :: ray(szib_(g),szk_(g)) ! Ray is the Rayleigh-drag velocity [Z T-1 ~> m s-1].
195  real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2].
196 
197  real :: hmix ! The mixed layer thickness over which stress
198  ! is applied with direct_stress [H ~> m or kg m-2].
199  real :: i_hmix ! The inverse of Hmix [H-1 ~> m-1 or m2 kg-1].
200  real :: idt ! The inverse of the time step [T-1 ~> s-1].
201  real :: dt_rho0 ! The time step divided by the mean density [T H Z-1 R-1 ~> s m3 kg-1 or s].
202  real :: dt_z_to_h ! The time step times the conversion from Z to the
203  ! units of thickness - [T H Z-1 ~> s or s kg m-3].
204  real :: h_neglect ! A thickness that is so small it is usually lost
205  ! in roundoff and can be neglected [H ~> m or kg m-2].
206 
207  real :: stress ! The surface stress times the time step, divided
208  ! by the density [H L T-1 ~> m2 s-1 or kg m-1 s-1].
209  real :: zds, hfr, h_a ! Temporary variables used with direct_stress.
210  real :: surface_stress(szib_(g))! The same as stress, unless the wind stress
211  ! stress is applied as a body force [H L T-1 ~> m2 s-1 or kg m-1 s-1].
212 
213  real, allocatable, dimension(:,:) :: hf_du_dt_visc_2d ! Depth sum of hf_du_dt_visc [L T-2 ~> m s-2]
214  real, allocatable, dimension(:,:) :: hf_dv_dt_visc_2d ! Depth sum of hf_dv_dt_visc [L T-2 ~> m s-2]
215 
216  logical :: do_i(szib_(g))
217  logical :: dostokesmixing
218 
219  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, n
220  is = g%isc ; ie = g%iec; js = g%jsc; je = g%jec
221  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
222 
223  if (.not.associated(cs)) call mom_error(fatal,"MOM_vert_friction(visc): "// &
224  "Module must be initialized before it is used.")
225 
226  if (cs%direct_stress) then
227  hmix = cs%Hmix_stress
228  i_hmix = 1.0 / hmix
229  endif
230  dt_rho0 = dt / gv%H_to_RZ
231  dt_z_to_h = dt*gv%Z_to_H
232  h_neglect = gv%H_subroundoff
233  idt = 1.0 / dt
234 
235  !Check if Stokes mixing allowed if requested (present and associated)
236  dostokesmixing=.false.
237  if (cs%StokesMixing) then
238  if (present(waves)) dostokesmixing = associated(waves)
239  if (.not. dostokesmixing) &
240  call mom_error(fatal,"Stokes Mixing called without allocated"//&
241  "Waves Control Structure")
242  endif
243 
244  do k=1,nz ; do i=isq,ieq ; ray(i,k) = 0.0 ; enddo ; enddo
245 
246  ! Update the zonal velocity component using a modification of a standard
247  ! tridagonal solver.
248 
249  !$OMP parallel do default(shared) firstprivate(Ray) &
250  !$OMP private(do_i,surface_stress,zDS,stress,h_a,hfr, &
251  !$OMP b_denom_1,b1,d1,c1)
252  do j=g%jsc,g%jec
253  do i=isq,ieq ; do_i(i) = (g%mask2dCu(i,j) > 0) ; enddo
254 
255  ! When mixing down Eulerian current + Stokes drift add before calling solver
256  if (dostokesmixing) then ; do k=1,nz ; do i=isq,ieq
257  if (do_i(i)) u(i,j,k) = u(i,j,k) + us%m_s_to_L_T*waves%Us_x(i,j,k)
258  enddo ; enddo ; endif
259 
260  if (associated(adp%du_dt_visc)) then ; do k=1,nz ; do i=isq,ieq
261  adp%du_dt_visc(i,j,k) = u(i,j,k)
262  enddo ; enddo ; endif
263 
264 ! One option is to have the wind stress applied as a body force
265 ! over the topmost Hmix fluid. If DIRECT_STRESS is not defined,
266 ! the wind stress is applied as a stress boundary condition.
267  if (cs%direct_stress) then
268  do i=isq,ieq ; if (do_i(i)) then
269  surface_stress(i) = 0.0
270  zds = 0.0
271  stress = dt_rho0 * forces%taux(i,j)
272  do k=1,nz
273  h_a = 0.5 * (h(i,j,k) + h(i+1,j,k)) + h_neglect
274  hfr = 1.0 ; if ((zds+h_a) > hmix) hfr = (hmix - zds) / h_a
275  u(i,j,k) = u(i,j,k) + i_hmix * hfr * stress
276  zds = zds + h_a ; if (zds >= hmix) exit
277  enddo
278  endif ; enddo ! end of i loop
279  else ; do i=isq,ieq
280  surface_stress(i) = dt_rho0 * (g%mask2dCu(i,j)*forces%taux(i,j))
281  enddo ; endif ! direct_stress
282 
283  if (cs%Channel_drag) then ; do k=1,nz ; do i=isq,ieq
284  ray(i,k) = visc%Ray_u(i,j,k)
285  enddo ; enddo ; endif
286 
287  ! perform forward elimination on the tridiagonal system
288  !
289  ! denote the diagonal of the system as b_k, the subdiagonal as a_k
290  ! and the superdiagonal as c_k. The right-hand side terms are d_k.
291  !
292  ! ignoring the rayleigh drag contribution,
293  ! we have a_k = -dt_Z_to_H * a_u(k)
294  ! b_k = h_u(k) + dt_Z_to_H * (a_u(k) + a_u(k+1))
295  ! c_k = -dt_Z_to_H * a_u(k+1)
296  !
297  ! for forward elimination, we want to:
298  ! calculate c'_k = - c_k / (b_k + a_k c'_(k-1))
299  ! and d'_k = (d_k - a_k d'_(k-1)) / (b_k + a_k c'_(k-1))
300  ! where c'_1 = c_1/b_1 and d'_1 = d_1/b_1
301  ! (see Thomas' tridiagonal matrix algorithm)
302  !
303  ! b1 is the denominator term 1 / (b_k + a_k c'_(k-1))
304  ! b_denom_1 is (b_k + a_k + c_k) - a_k(1 - c'_(k-1))
305  ! = (b_k + c_k + c'_(k-1))
306  ! this is done so that d1 = b1 * b_denom_1 = 1 - c'_(k-1)
307  ! c1(k) is -c'_(k - 1)
308  ! and the right-hand-side is destructively updated to be d'_k
309  !
310  do i=isq,ieq ; if (do_i(i)) then
311  b_denom_1 = cs%h_u(i,j,1) + dt_z_to_h * (ray(i,1) + cs%a_u(i,j,1))
312  b1(i) = 1.0 / (b_denom_1 + dt_z_to_h*cs%a_u(i,j,2))
313  d1(i) = b_denom_1 * b1(i)
314  u(i,j,1) = b1(i) * (cs%h_u(i,j,1) * u(i,j,1) + surface_stress(i))
315  endif ; enddo
316  do k=2,nz ; do i=isq,ieq ; if (do_i(i)) then
317  c1(i,k) = dt_z_to_h * cs%a_u(i,j,k) * b1(i)
318  b_denom_1 = cs%h_u(i,j,k) + dt_z_to_h * (ray(i,k) + cs%a_u(i,j,k)*d1(i))
319  b1(i) = 1.0 / (b_denom_1 + dt_z_to_h * cs%a_u(i,j,k+1))
320  d1(i) = b_denom_1 * b1(i)
321  u(i,j,k) = (cs%h_u(i,j,k) * u(i,j,k) + &
322  dt_z_to_h * cs%a_u(i,j,k) * u(i,j,k-1)) * b1(i)
323  endif ; enddo ; enddo
324 
325  ! back substitute to solve for the new velocities
326  ! u_k = d'_k - c'_k x_(k+1)
327  do k=nz-1,1,-1 ; do i=isq,ieq ; if (do_i(i)) then
328  u(i,j,k) = u(i,j,k) + c1(i,k+1) * u(i,j,k+1)
329  endif ; enddo ; enddo ! i and k loops
330 
331  if (associated(adp%du_dt_visc)) then ; do k=1,nz ; do i=isq,ieq
332  adp%du_dt_visc(i,j,k) = (u(i,j,k) - adp%du_dt_visc(i,j,k))*idt
333  enddo ; enddo ; endif
334 
335  if (associated(visc%taux_shelf)) then ; do i=isq,ieq
336  visc%taux_shelf(i,j) = -gv%Rho0*cs%a1_shelf_u(i,j)*u(i,j,1) ! - u_shelf?
337  enddo ; endif
338 
339  if (PRESENT(taux_bot)) then
340  do i=isq,ieq
341  taux_bot(i,j) = gv%Rho0 * (u(i,j,nz)*cs%a_u(i,j,nz+1))
342  enddo
343  if (cs%Channel_drag) then ; do k=1,nz ; do i=isq,ieq
344  taux_bot(i,j) = taux_bot(i,j) + gv%Rho0 * (ray(i,k)*u(i,j,k))
345  enddo ; enddo ; endif
346  endif
347 
348  ! When mixing down Eulerian current + Stokes drift subtract after calling solver
349  if (dostokesmixing) then ; do k=1,nz ; do i=isq,ieq
350  if (do_i(i)) u(i,j,k) = u(i,j,k) - us%m_s_to_L_T*waves%Us_x(i,j,k)
351  enddo ; enddo ; endif
352 
353  enddo ! end u-component j loop
354 
355  ! Now work on the meridional velocity component.
356 
357  !$OMP parallel do default(shared) firstprivate(Ray) &
358  !$OMP private(do_i,surface_stress,zDS,stress,h_a,hfr, &
359  !$OMP b_denom_1,b1,d1,c1)
360  do j=jsq,jeq
361  do i=is,ie ; do_i(i) = (g%mask2dCv(i,j) > 0) ; enddo
362 
363  ! When mixing down Eulerian current + Stokes drift add before calling solver
364  if (dostokesmixing) then ; do k=1,nz ; do i=is,ie
365  if (do_i(i)) v(i,j,k) = v(i,j,k) + us%m_s_to_L_T*waves%Us_y(i,j,k)
366  enddo ; enddo ; endif
367 
368  if (associated(adp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie
369  adp%dv_dt_visc(i,j,k) = v(i,j,k)
370  enddo ; enddo ; endif
371 
372 ! One option is to have the wind stress applied as a body force
373 ! over the topmost Hmix fluid. If DIRECT_STRESS is not defined,
374 ! the wind stress is applied as a stress boundary condition.
375  if (cs%direct_stress) then
376  do i=is,ie ; if (do_i(i)) then
377  surface_stress(i) = 0.0
378  zds = 0.0
379  stress = dt_rho0 * forces%tauy(i,j)
380  do k=1,nz
381  h_a = 0.5 * (h(i,j,k) + h(i,j+1,k)) + h_neglect
382  hfr = 1.0 ; if ((zds+h_a) > hmix) hfr = (hmix - zds) / h_a
383  v(i,j,k) = v(i,j,k) + i_hmix * hfr * stress
384  zds = zds + h_a ; if (zds >= hmix) exit
385  enddo
386  endif ; enddo ! end of i loop
387  else ; do i=is,ie
388  surface_stress(i) = dt_rho0 * (g%mask2dCv(i,j)*forces%tauy(i,j))
389  enddo ; endif ! direct_stress
390 
391  if (cs%Channel_drag) then ; do k=1,nz ; do i=is,ie
392  ray(i,k) = visc%Ray_v(i,j,k)
393  enddo ; enddo ; endif
394 
395  do i=is,ie ; if (do_i(i)) then
396  b_denom_1 = cs%h_v(i,j,1) + dt_z_to_h * (ray(i,1) + cs%a_v(i,j,1))
397  b1(i) = 1.0 / (b_denom_1 + dt_z_to_h*cs%a_v(i,j,2))
398  d1(i) = b_denom_1 * b1(i)
399  v(i,j,1) = b1(i) * (cs%h_v(i,j,1) * v(i,j,1) + surface_stress(i))
400  endif ; enddo
401  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
402  c1(i,k) = dt_z_to_h * cs%a_v(i,j,k) * b1(i)
403  b_denom_1 = cs%h_v(i,j,k) + dt_z_to_h * (ray(i,k) + cs%a_v(i,j,k)*d1(i))
404  b1(i) = 1.0 / (b_denom_1 + dt_z_to_h * cs%a_v(i,j,k+1))
405  d1(i) = b_denom_1 * b1(i)
406  v(i,j,k) = (cs%h_v(i,j,k) * v(i,j,k) + dt_z_to_h * cs%a_v(i,j,k) * v(i,j,k-1)) * b1(i)
407  endif ; enddo ; enddo
408  do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then
409  v(i,j,k) = v(i,j,k) + c1(i,k+1) * v(i,j,k+1)
410  endif ; enddo ; enddo ! i and k loops
411 
412  if (associated(adp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie
413  adp%dv_dt_visc(i,j,k) = (v(i,j,k) - adp%dv_dt_visc(i,j,k))*idt
414  enddo ; enddo ; endif
415 
416  if (associated(visc%tauy_shelf)) then ; do i=is,ie
417  visc%tauy_shelf(i,j) = -gv%Rho0*cs%a1_shelf_v(i,j)*v(i,j,1) ! - v_shelf?
418  enddo ; endif
419 
420  if (present(tauy_bot)) then
421  do i=is,ie
422  tauy_bot(i,j) = gv%Rho0 * (v(i,j,nz)*cs%a_v(i,j,nz+1))
423  enddo
424  if (cs%Channel_drag) then ; do k=1,nz ; do i=is,ie
425  tauy_bot(i,j) = tauy_bot(i,j) + gv%Rho0 * (ray(i,k)*v(i,j,k))
426  enddo ; enddo ; endif
427  endif
428 
429  ! When mixing down Eulerian current + Stokes drift subtract after calling solver
430  if (dostokesmixing) then ; do k=1,nz ; do i=is,ie
431  if (do_i(i)) v(i,j,k) = v(i,j,k) - us%m_s_to_L_T*waves%Us_y(i,j,k)
432  enddo ; enddo ; endif
433 
434  enddo ! end of v-component J loop
435 
436  call vertvisc_limit_vel(u, v, h, adp, cdp, forces, visc, dt, g, gv, us, cs)
437 
438  ! Here the velocities associated with open boundary conditions are applied.
439  if (associated(obc)) then
440  do n=1,obc%number_of_segments
441  if (obc%segment(n)%specified) then
442  if (obc%segment(n)%is_N_or_S) then
443  j = obc%segment(n)%HI%JsdB
444  do k=1,nz ; do i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied
445  v(i,j,k) = obc%segment(n)%normal_vel(i,j,k)
446  enddo ; enddo
447  elseif (obc%segment(n)%is_E_or_W) then
448  i = obc%segment(n)%HI%IsdB
449  do k=1,nz ; do j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed
450  u(i,j,k) = obc%segment(n)%normal_vel(i,j,k)
451  enddo ; enddo
452  endif
453  endif
454  enddo
455  endif
456 
457 ! Offer diagnostic fields for averaging.
458  if (cs%id_du_dt_visc > 0) &
459  call post_data(cs%id_du_dt_visc, adp%du_dt_visc, cs%diag)
460  if (cs%id_dv_dt_visc > 0) &
461  call post_data(cs%id_dv_dt_visc, adp%dv_dt_visc, cs%diag)
462  if (present(taux_bot) .and. (cs%id_taux_bot > 0)) &
463  call post_data(cs%id_taux_bot, taux_bot, cs%diag)
464  if (present(tauy_bot) .and. (cs%id_tauy_bot > 0)) &
465  call post_data(cs%id_tauy_bot, tauy_bot, cs%diag)
466 
467  ! Diagnostics for terms multiplied by fractional thicknesses
468 
469  ! 3D diagnostics hf_du(dv)_dt_visc are commented because there is no clarity on proper remapping grid option.
470  ! The code is retained for degugging purposes in the future.
471  !if (CS%id_hf_du_dt_visc > 0) then
472  ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq
473  ! CS%hf_du_dt_visc(I,j,k) = ADp%du_dt_visc(I,j,k) * ADp%diag_hfrac_u(I,j,k)
474  ! enddo ; enddo ; enddo
475  ! call post_data(CS%id_hf_du_dt_visc, CS%hf_du_dt_visc, CS%diag)
476  !endif
477  !if (CS%id_hf_dv_dt_visc > 0) then
478  ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
479  ! CS%hf_dv_dt_visc(i,J,k) = ADp%dv_dt_visc(i,J,k) * ADp%diag_hfrac_v(i,J,k)
480  ! enddo ; enddo ; enddo
481  ! call post_data(CS%id_hf_dv_dt_visc, CS%hf_dv_dt_visc, CS%diag)
482  !endif
483  if (cs%id_hf_du_dt_visc_2d > 0) then
484  allocate(hf_du_dt_visc_2d(g%IsdB:g%IedB,g%jsd:g%jed))
485  hf_du_dt_visc_2d(:,:) = 0.0
486  do k=1,nz ; do j=js,je ; do i=isq,ieq
487  hf_du_dt_visc_2d(i,j) = hf_du_dt_visc_2d(i,j) + adp%du_dt_visc(i,j,k) * adp%diag_hfrac_u(i,j,k)
488  enddo ; enddo ; enddo
489  call post_data(cs%id_hf_du_dt_visc_2d, hf_du_dt_visc_2d, cs%diag)
490  deallocate(hf_du_dt_visc_2d)
491  endif
492  if (cs%id_hf_dv_dt_visc_2d > 0) then
493  allocate(hf_dv_dt_visc_2d(g%isd:g%ied,g%JsdB:g%JedB))
494  hf_dv_dt_visc_2d(:,:) = 0.0
495  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
496  hf_dv_dt_visc_2d(i,j) = hf_dv_dt_visc_2d(i,j) + adp%dv_dt_visc(i,j,k) * adp%diag_hfrac_v(i,j,k)
497  enddo ; enddo ; enddo
498  call post_data(cs%id_hf_dv_dt_visc_2d, hf_dv_dt_visc_2d, cs%diag)
499  deallocate(hf_dv_dt_visc_2d)
500  endif
501 
502 end subroutine vertvisc
503 
504 !> Calculate the fraction of momentum originally in a layer that remains
505 !! after a time-step of viscosity, and the fraction of a time-step's
506 !! worth of barotropic acceleration that a layer experiences after
507 !! viscosity is applied.
508 subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS)
509  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
510  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
511  type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag
512  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
513  intent(inout) :: visc_rem_u !< Fraction of a time-step's worth of a
514  !! barotopic acceleration that a layer experiences after
515  !! viscosity is applied in the zonal direction [nondim]
516  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
517  intent(inout) :: visc_rem_v !< Fraction of a time-step's worth of a
518  !! barotopic acceleration that a layer experiences after
519  !! viscosity is applied in the meridional direction [nondim]
520  real, intent(in) :: dt !< Time increment [T ~> s]
521  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
522  type(vertvisc_cs), pointer :: cs !< Vertical viscosity control structure
523 
524  ! Local variables
525 
526  real :: b1(szib_(g)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
527  real :: c1(szib_(g),szk_(g)) ! A variable used by the tridiagonal solver [nondim].
528  real :: d1(szib_(g)) ! d1=1-c1 is used by the tridiagonal solver [nondim].
529  real :: ray(szib_(g),szk_(g)) ! Ray is the Rayleigh-drag velocity [Z T-1 ~> m s-1].
530  real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2].
531  real :: dt_z_to_h ! The time step times the conversion from Z to the
532  ! units of thickness [T H Z-1 ~> s or s kg m-3].
533  logical :: do_i(szib_(g))
534 
535  integer :: i, j, k, is, ie, isq, ieq, jsq, jeq, nz
536  is = g%isc ; ie = g%iec
537  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
538 
539  if (.not.associated(cs)) call mom_error(fatal,"MOM_vert_friction(visc): "// &
540  "Module must be initialized before it is used.")
541 
542  dt_z_to_h = dt*gv%Z_to_H
543 
544  do k=1,nz ; do i=isq,ieq ; ray(i,k) = 0.0 ; enddo ; enddo
545 
546  ! Find the zonal viscous using a modification of a standard tridagonal solver.
547 !$OMP parallel do default(none) shared(G,Isq,Ieq,CS,nz,visc,dt_Z_to_H,visc_rem_u) &
548 !$OMP firstprivate(Ray) &
549 !$OMP private(do_i,b_denom_1,b1,d1,c1)
550  do j=g%jsc,g%jec
551  do i=isq,ieq ; do_i(i) = (g%mask2dCu(i,j) > 0) ; enddo
552 
553  if (cs%Channel_drag) then ; do k=1,nz ; do i=isq,ieq
554  ray(i,k) = visc%Ray_u(i,j,k)
555  enddo ; enddo ; endif
556 
557  do i=isq,ieq ; if (do_i(i)) then
558  b_denom_1 = cs%h_u(i,j,1) + dt_z_to_h * (ray(i,1) + cs%a_u(i,j,1))
559  b1(i) = 1.0 / (b_denom_1 + dt_z_to_h*cs%a_u(i,j,2))
560  d1(i) = b_denom_1 * b1(i)
561  visc_rem_u(i,j,1) = b1(i) * cs%h_u(i,j,1)
562  endif ; enddo
563  do k=2,nz ; do i=isq,ieq ; if (do_i(i)) then
564  c1(i,k) = dt_z_to_h * cs%a_u(i,j,k)*b1(i)
565  b_denom_1 = cs%h_u(i,j,k) + dt_z_to_h * (ray(i,k) + cs%a_u(i,j,k)*d1(i))
566  b1(i) = 1.0 / (b_denom_1 + dt_z_to_h * cs%a_u(i,j,k+1))
567  d1(i) = b_denom_1 * b1(i)
568  visc_rem_u(i,j,k) = (cs%h_u(i,j,k) + dt_z_to_h * cs%a_u(i,j,k) * visc_rem_u(i,j,k-1)) * b1(i)
569  endif ; enddo ; enddo
570  do k=nz-1,1,-1 ; do i=isq,ieq ; if (do_i(i)) then
571  visc_rem_u(i,j,k) = visc_rem_u(i,j,k) + c1(i,k+1)*visc_rem_u(i,j,k+1)
572 
573  endif ; enddo ; enddo ! i and k loops
574 
575  enddo ! end u-component j loop
576 
577  ! Now find the meridional viscous using a modification.
578 !$OMP parallel do default(none) shared(Jsq,Jeq,is,ie,G,CS,visc,dt_Z_to_H,visc_rem_v,nz) &
579 !$OMP firstprivate(Ray) &
580 !$OMP private(do_i,b_denom_1,b1,d1,c1)
581  do j=jsq,jeq
582  do i=is,ie ; do_i(i) = (g%mask2dCv(i,j) > 0) ; enddo
583 
584  if (cs%Channel_drag) then ; do k=1,nz ; do i=is,ie
585  ray(i,k) = visc%Ray_v(i,j,k)
586  enddo ; enddo ; endif
587 
588  do i=is,ie ; if (do_i(i)) then
589  b_denom_1 = cs%h_v(i,j,1) + dt_z_to_h * (ray(i,1) + cs%a_v(i,j,1))
590  b1(i) = 1.0 / (b_denom_1 + dt_z_to_h*cs%a_v(i,j,2))
591  d1(i) = b_denom_1 * b1(i)
592  visc_rem_v(i,j,1) = b1(i) * cs%h_v(i,j,1)
593  endif ; enddo
594  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
595  c1(i,k) = dt_z_to_h * cs%a_v(i,j,k)*b1(i)
596  b_denom_1 = cs%h_v(i,j,k) + dt_z_to_h * (ray(i,k) + cs%a_v(i,j,k)*d1(i))
597  b1(i) = 1.0 / (b_denom_1 + dt_z_to_h * cs%a_v(i,j,k+1))
598  d1(i) = b_denom_1 * b1(i)
599  visc_rem_v(i,j,k) = (cs%h_v(i,j,k) + dt_z_to_h * cs%a_v(i,j,k) * visc_rem_v(i,j,k-1)) * b1(i)
600  endif ; enddo ; enddo
601  do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then
602  visc_rem_v(i,j,k) = visc_rem_v(i,j,k) + c1(i,k+1)*visc_rem_v(i,j,k+1)
603  endif ; enddo ; enddo ! i and k loops
604  enddo ! end of v-component J loop
605 
606  if (cs%debug) then
607  call uvchksum("visc_rem_[uv]", visc_rem_u, visc_rem_v, g%HI, haloshift=0, &
608  scalar_pair=.true.)
609  endif
610 
611 end subroutine vertvisc_remnant
612 
613 
614 !> Calculate the coupling coefficients (CS%a_u and CS%a_v)
615 !! and effective layer thicknesses (CS%h_u and CS%h_v) for later use in the
616 !! applying the implicit vertical viscosity via vertvisc().
617 subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
618  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
619  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
620  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
621  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
622  intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]
623  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
624  intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]
625  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
626  intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
627  type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
628  type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag
629  real, intent(in) :: dt !< Time increment [T ~> s]
630  type(vertvisc_cs), pointer :: cs !< Vertical viscosity control structure
631  type(ocean_obc_type), pointer :: obc !< Open boundary condition structure
632 
633  ! Field from forces used in this subroutine:
634  ! ustar: the friction velocity [Z T-1 ~> m s-1], used here as the mixing
635  ! velocity in the mixed layer if NKML > 1 in a bulk mixed layer.
636 
637  ! Local variables
638 
639  real, dimension(SZIB_(G),SZK_(G)) :: &
640  h_harm, & ! Harmonic mean of the thicknesses around a velocity grid point,
641  ! given by 2*(h+ * h-)/(h+ + h-) [H ~> m or kg m-2].
642  h_arith, & ! The arithmetic mean thickness [H ~> m or kg m-2].
643  h_delta, & ! The lateral difference of thickness [H ~> m or kg m-2].
644  hvel, & ! hvel is the thickness used at a velocity grid point [H ~> m or kg m-2].
645  hvel_shelf ! The equivalent of hvel under shelves [H ~> m or kg m-2].
646  real, dimension(SZIB_(G),SZK_(G)+1) :: &
647  a_cpl, & ! The drag coefficients across interfaces [Z T-1 ~> m s-1]. a_cpl times
648  ! the velocity difference gives the stress across an interface.
649  a_shelf, & ! The drag coefficients across interfaces in water columns under
650  ! ice shelves [Z T-1 ~> m s-1].
651  z_i ! An estimate of each interface's height above the bottom,
652  ! normalized by the bottom boundary layer thickness, nondim.
653  real, dimension(SZIB_(G)) :: &
654  kv_bbl, & ! The bottom boundary layer viscosity [Z2 T-1 ~> m2 s-1].
655  bbl_thick, & ! The bottom boundary layer thickness [H ~> m or kg m-2].
656  i_hbbl, & ! The inverse of the bottom boundary layer thickness [H-1 ~> m-1 or m2 kg-1].
657  i_htbl, & ! The inverse of the top boundary layer thickness [H-1 ~> m-1 or m2 kg-1].
658  zcol1, & ! The height of the interfaces to the north and south of a
659  zcol2, & ! v-point [H ~> m or kg m-2].
660  ztop_min, & ! The deeper of the two adjacent surface heights [H ~> m or kg m-2].
661  dmin, & ! The shallower of the two adjacent bottom depths converted to
662  ! thickness units [H ~> m or kg m-2].
663  zh, & ! An estimate of the interface's distance from the bottom
664  ! based on harmonic mean thicknesses [H ~> m or kg m-2].
665  h_ml ! The mixed layer depth [H ~> m or kg m-2].
666  real, allocatable, dimension(:,:) :: hml_u ! Diagnostic of the mixed layer depth at u points [H ~> m or kg m-2].
667  real, allocatable, dimension(:,:) :: hml_v ! Diagnostic of the mixed layer depth at v points [H ~> m or kg m-2].
668  real, allocatable, dimension(:,:,:) :: kv_u !< Total vertical viscosity at u-points [Z2 T-1 ~> m2 s-1].
669  real, allocatable, dimension(:,:,:) :: kv_v !< Total vertical viscosity at v-points [Z2 T-1 ~> m2 s-1].
670  real :: zcol(szi_(g)) ! The height of an interface at h-points [H ~> m or kg m-2].
671  real :: botfn ! A function which goes from 1 at the bottom to 0 much more
672  ! than Hbbl into the interior.
673  real :: topfn ! A function which goes from 1 at the top to 0 much more
674  ! than Htbl into the interior.
675  real :: z2 ! The distance from the bottom, normalized by Hbbl, nondim.
676  real :: z2_wt ! A nondimensional (0-1) weight used when calculating z2.
677  real :: z_clear ! The clearance of an interface above the surrounding topography [H ~> m or kg m-2].
678  real :: a_cpl_max ! The maximum drag doefficient across interfaces, set so that it will be
679  ! representable as a 32-bit float in MKS units [Z T-1 ~> m s-1]
680  real :: h_neglect ! A thickness that is so small it is usually lost
681  ! in roundoff and can be neglected [H ~> m or kg m-2].
682 
683  real :: i_valbl ! The inverse of a scaling factor determining when water is
684  ! still within the boundary layer, as determined by the sum
685  ! of the harmonic mean thicknesses.
686  logical, dimension(SZIB_(G)) :: do_i, do_i_shelf
687  logical :: do_any_shelf
688  integer, dimension(SZIB_(G)) :: &
689  zi_dir ! A trinary logical array indicating which thicknesses to use for
690  ! finding z_clear.
691  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
692  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
693  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
694 
695  if (.not.associated(cs)) call mom_error(fatal,"MOM_vert_friction(coef): "// &
696  "Module must be initialized before it is used.")
697 
698  h_neglect = gv%H_subroundoff
699  a_cpl_max = 1.0e37 * us%m_to_Z * us%T_to_s
700  i_hbbl(:) = 1.0 / (cs%Hbbl + h_neglect)
701  i_valbl = 0.0 ; if (cs%harm_BL_val > 0.0) i_valbl = 1.0 / cs%harm_BL_val
702 
703  if (cs%id_Kv_u > 0) then
704  allocate(kv_u(g%IsdB:g%IedB,g%jsd:g%jed,g%ke)) ; kv_u(:,:,:) = 0.0
705  endif
706 
707  if (cs%id_Kv_v > 0) then
708  allocate(kv_v(g%isd:g%ied,g%JsdB:g%JedB,g%ke)) ; kv_v(:,:,:) = 0.0
709  endif
710 
711  if (cs%debug .or. (cs%id_hML_u > 0)) then
712  allocate(hml_u(g%IsdB:g%IedB,g%jsd:g%jed)) ; hml_u(:,:) = 0.0
713  endif
714  if (cs%debug .or. (cs%id_hML_v > 0)) then
715  allocate(hml_v(g%isd:g%ied,g%JsdB:g%JedB)) ; hml_v(:,:) = 0.0
716  endif
717 
718  if ((associated(visc%taux_shelf) .or. associated(forces%frac_shelf_u)) .and. &
719  .not.associated(cs%a1_shelf_u)) then
720  allocate(cs%a1_shelf_u(g%IsdB:g%IedB,g%jsd:g%jed)) ; cs%a1_shelf_u(:,:)=0.0
721  endif
722  if ((associated(visc%tauy_shelf) .or. associated(forces%frac_shelf_v)) .and. &
723  .not.associated(cs%a1_shelf_v)) then
724  allocate(cs%a1_shelf_v(g%isd:g%ied,g%JsdB:g%JedB)) ; cs%a1_shelf_v(:,:)=0.0
725  endif
726 
727  !$OMP parallel do default(private) shared(G,GV,US,CS,visc,Isq,Ieq,nz,u,h,forces,hML_u, &
728  !$OMP OBC,h_neglect,dt,I_valBL,Kv_u,a_cpl_max) &
729  !$OMP firstprivate(i_hbbl)
730  do j=g%Jsc,g%Jec
731  do i=isq,ieq ; do_i(i) = (g%mask2dCu(i,j) > 0) ; enddo
732 
733  if (cs%bottomdraglaw) then ; do i=isq,ieq
734  kv_bbl(i) = visc%Kv_bbl_u(i,j)
735  bbl_thick(i) = visc%bbl_thick_u(i,j) * gv%Z_to_H + h_neglect
736  if (do_i(i)) i_hbbl(i) = 1.0 / bbl_thick(i)
737  enddo ; endif
738 
739  do k=1,nz ; do i=isq,ieq ; if (do_i(i)) then
740  h_harm(i,k) = 2.0*h(i,j,k)*h(i+1,j,k) / (h(i,j,k)+h(i+1,j,k)+h_neglect)
741  h_arith(i,k) = 0.5*(h(i+1,j,k)+h(i,j,k))
742  h_delta(i,k) = h(i+1,j,k) - h(i,j,k)
743  endif ; enddo ; enddo
744  do i=isq,ieq
745  dmin(i) = min(g%bathyT(i,j), g%bathyT(i+1,j)) * gv%Z_to_H
746  zi_dir(i) = 0
747  enddo
748 
749  ! Project thickness outward across OBCs using a zero-gradient condition.
750  if (associated(obc)) then ; if (obc%number_of_segments > 0) then
751  do i=isq,ieq ; if (do_i(i) .and. (obc%segnum_u(i,j) /= obc_none)) then
752  if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
753  do k=1,nz ; h_harm(i,k) = h(i,j,k) ; h_arith(i,k) = h(i,j,k) ; h_delta(i,k) = 0. ; enddo
754  dmin(i) = g%bathyT(i,j) * gv%Z_to_H
755  zi_dir(i) = -1
756  elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
757  do k=1,nz ; h_harm(i,k) = h(i+1,j,k) ; h_arith(i,k) = h(i+1,j,k) ; h_delta(i,k) = 0. ; enddo
758  dmin(i) = g%bathyT(i+1,j) * gv%Z_to_H
759  zi_dir(i) = 1
760  endif
761  endif ; enddo
762  endif ; endif
763 
764 ! The following block calculates the thicknesses at velocity
765 ! grid points for the vertical viscosity (hvel). Near the
766 ! bottom an upwind biased thickness is used to control the effect
767 ! of spurious Montgomery potential gradients at the bottom where
768 ! nearly massless layers layers ride over the topography.
769  if (cs%harmonic_visc) then
770  do i=isq,ieq ; z_i(i,nz+1) = 0.0 ; enddo
771  do k=nz,1,-1 ; do i=isq,ieq ; if (do_i(i)) then
772  hvel(i,k) = h_harm(i,k)
773  if (u(i,j,k) * h_delta(i,k) < 0) then
774  z2 = z_i(i,k+1) ; botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
775  hvel(i,k) = (1.0-botfn)*h_harm(i,k) + botfn*h_arith(i,k)
776  endif
777  z_i(i,k) = z_i(i,k+1) + h_harm(i,k)*i_hbbl(i)
778  endif ; enddo ; enddo ! i & k loops
779  else ! Not harmonic_visc
780  do i=isq,ieq ; zh(i) = 0.0 ; z_i(i,nz+1) = 0.0 ; enddo
781  do i=isq,ieq+1 ; zcol(i) = -g%bathyT(i,j) * gv%Z_to_H ; enddo
782  do k=nz,1,-1
783  do i=isq,ieq+1 ; zcol(i) = zcol(i) + h(i,j,k) ; enddo
784  do i=isq,ieq ; if (do_i(i)) then
785  zh(i) = zh(i) + h_harm(i,k)
786 
787  z_clear = max(zcol(i),zcol(i+1)) + dmin(i)
788  if (zi_dir(i) < 0) z_clear = zcol(i) + dmin(i)
789  if (zi_dir(i) > 0) z_clear = zcol(i+1) + dmin(i)
790 
791  z_i(i,k) = max(zh(i), z_clear) * i_hbbl(i)
792 
793  hvel(i,k) = h_arith(i,k)
794  if (u(i,j,k) * h_delta(i,k) > 0) then
795  if (zh(i) * i_hbbl(i) < cs%harm_BL_val) then
796  hvel(i,k) = h_harm(i,k)
797  else
798  z2_wt = 1.0 ; if (zh(i) * i_hbbl(i) < 2.0*cs%harm_BL_val) &
799  z2_wt = max(0.0, min(1.0, zh(i) * i_hbbl(i) * i_valbl - 1.0))
800  z2 = z2_wt * (max(zh(i), z_clear) * i_hbbl(i))
801  botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
802  hvel(i,k) = (1.0-botfn)*h_arith(i,k) + botfn*h_harm(i,k)
803  endif
804  endif
805 
806  endif ; enddo ! i loop
807  enddo ! k loop
808  endif
809 
810  call find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_ml, &
811  dt, j, g, gv, us, cs, visc, forces, work_on_u=.true., obc=obc)
812  if (allocated(hml_u)) then
813  do i=isq,ieq ; if (do_i(i)) then ; hml_u(i,j) = h_ml(i) ; endif ; enddo
814  endif
815 
816  do_any_shelf = .false.
817  if (associated(forces%frac_shelf_u)) then
818  do i=isq,ieq
819  cs%a1_shelf_u(i,j) = 0.0
820  do_i_shelf(i) = (do_i(i) .and. forces%frac_shelf_u(i,j) > 0.0)
821  if (do_i_shelf(i)) do_any_shelf = .true.
822  enddo
823  if (do_any_shelf) then
824  if (cs%harmonic_visc) then
825  call find_coupling_coef(a_shelf, hvel, do_i_shelf, h_harm, bbl_thick, &
826  kv_bbl, z_i, h_ml, dt, j, g, gv, us, cs, &
827  visc, forces, work_on_u=.true., obc=obc, shelf=.true.)
828  else ! Find upwind-biased thickness near the surface.
829  ! Perhaps this needs to be done more carefully, via find_eta.
830  do i=isq,ieq ; if (do_i_shelf(i)) then
831  zh(i) = 0.0 ; ztop_min(i) = min(zcol(i), zcol(i+1))
832  i_htbl(i) = 1.0 / (visc%tbl_thick_shelf_u(i,j)*gv%Z_to_H + h_neglect)
833  endif ; enddo
834  do k=1,nz
835  do i=isq,ieq+1 ; zcol(i) = zcol(i) - h(i,j,k) ; enddo
836  do i=isq,ieq ; if (do_i_shelf(i)) then
837  zh(i) = zh(i) + h_harm(i,k)
838 
839  hvel_shelf(i,k) = hvel(i,k)
840  if (u(i,j,k) * h_delta(i,k) > 0) then
841  if (zh(i) * i_htbl(i) < cs%harm_BL_val) then
842  hvel_shelf(i,k) = min(hvel(i,k), h_harm(i,k))
843  else
844  z2_wt = 1.0 ; if (zh(i) * i_htbl(i) < 2.0*cs%harm_BL_val) &
845  z2_wt = max(0.0, min(1.0, zh(i) * i_htbl(i) * i_valbl - 1.0))
846  z2 = z2_wt * (max(zh(i), ztop_min(i) - min(zcol(i),zcol(i+1))) * i_htbl(i))
847  topfn = 1.0 / (1.0 + 0.09*z2**6)
848  hvel_shelf(i,k) = min(hvel(i,k), (1.0-topfn)*h_arith(i,k) + topfn*h_harm(i,k))
849  endif
850  endif
851  endif ; enddo
852  enddo
853  call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, h_harm, &
854  bbl_thick, kv_bbl, z_i, h_ml, dt, j, g, gv, us, cs, &
855  visc, forces, work_on_u=.true., obc=obc, shelf=.true.)
856  endif
857  do i=isq,ieq ; if (do_i_shelf(i)) cs%a1_shelf_u(i,j) = a_shelf(i,1) ; enddo
858  endif
859  endif
860 
861  if (do_any_shelf) then
862  do k=1,nz+1 ; do i=isq,ieq ; if (do_i_shelf(i)) then
863  cs%a_u(i,j,k) = min(a_cpl_max, forces%frac_shelf_u(i,j) * a_shelf(i,k) + &
864  (1.0-forces%frac_shelf_u(i,j)) * a_cpl(i,k))
865 ! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH
866 ! CS%a_u(I,j,K) = min(a_cpl_max, forces%frac_shelf_u(I,j) * max(a_shelf(I,K), a_cpl(I,K)) + &
867 ! (1.0-forces%frac_shelf_u(I,j)) * a_cpl(I,K))
868  elseif (do_i(i)) then
869  cs%a_u(i,j,k) = min(a_cpl_max, a_cpl(i,k))
870  endif ; enddo ; enddo
871  do k=1,nz ; do i=isq,ieq ; if (do_i_shelf(i)) then
872  ! Should we instead take the inverse of the average of the inverses?
873  cs%h_u(i,j,k) = forces%frac_shelf_u(i,j) * hvel_shelf(i,k) + &
874  (1.0-forces%frac_shelf_u(i,j)) * hvel(i,k) + h_neglect
875  elseif (do_i(i)) then
876  cs%h_u(i,j,k) = hvel(i,k) + h_neglect
877  endif ; enddo ; enddo
878  else
879  do k=1,nz+1 ; do i=isq,ieq ; if (do_i(i)) cs%a_u(i,j,k) = min(a_cpl_max, a_cpl(i,k)) ; enddo ; enddo
880  do k=1,nz ; do i=isq,ieq ; if (do_i(i)) cs%h_u(i,j,k) = hvel(i,k) + h_neglect ; enddo ; enddo
881  endif
882 
883  ! Diagnose total Kv at u-points
884  if (cs%id_Kv_u > 0) then
885  do k=1,nz ; do i=isq,ieq
886  if (do_i(i)) kv_u(i,j,k) = 0.5 * gv%H_to_Z*(cs%a_u(i,j,k)+cs%a_u(i,j,k+1)) * cs%h_u(i,j,k)
887  enddo ; enddo
888  endif
889 
890  enddo
891 
892 
893  ! Now work on v-points.
894  !$OMP parallel do default(private) shared(G,GV,CS,US,visc,is,ie,Jsq,Jeq,nz,v,h,forces,hML_v, &
895  !$OMP OBC,h_neglect,dt,I_valBL,Kv_v,a_cpl_max) &
896  !$OMP firstprivate(i_hbbl)
897  do j=jsq,jeq
898  do i=is,ie ; do_i(i) = (g%mask2dCv(i,j) > 0) ; enddo
899 
900  if (cs%bottomdraglaw) then ; do i=is,ie
901  kv_bbl(i) = visc%Kv_bbl_v(i,j)
902  bbl_thick(i) = visc%bbl_thick_v(i,j) * gv%Z_to_H + h_neglect
903  if (do_i(i)) i_hbbl(i) = 1.0 / bbl_thick(i)
904  enddo ; endif
905 
906  do k=1,nz ; do i=is,ie ; if (do_i(i)) then
907  h_harm(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / (h(i,j,k)+h(i,j+1,k)+h_neglect)
908  h_arith(i,k) = 0.5*(h(i,j+1,k)+h(i,j,k))
909  h_delta(i,k) = h(i,j+1,k) - h(i,j,k)
910  endif ; enddo ; enddo
911  do i=is,ie
912  dmin(i) = min(g%bathyT(i,j), g%bathyT(i,j+1)) * gv%Z_to_H
913  zi_dir(i) = 0
914  enddo
915 
916  ! Project thickness outward across OBCs using a zero-gradient condition.
917  if (associated(obc)) then ; if (obc%number_of_segments > 0) then
918  do i=is,ie ; if (do_i(i) .and. (obc%segnum_v(i,j) /= obc_none)) then
919  if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
920  do k=1,nz ; h_harm(i,k) = h(i,j,k) ; h_arith(i,k) = h(i,j,k) ; h_delta(i,k) = 0. ; enddo
921  dmin(i) = g%bathyT(i,j) * gv%Z_to_H
922  zi_dir(i) = -1
923  elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
924  do k=1,nz ; h_harm(i,k) = h(i,j+1,k) ; h_arith(i,k) = h(i,j+1,k) ; h_delta(i,k) = 0. ; enddo
925  dmin(i) = g%bathyT(i,j+1) * gv%Z_to_H
926  zi_dir(i) = 1
927  endif
928  endif ; enddo
929  endif ; endif
930 
931 ! The following block calculates the thicknesses at velocity
932 ! grid points for the vertical viscosity (hvel). Near the
933 ! bottom an upwind biased thickness is used to control the effect
934 ! of spurious Montgomery potential gradients at the bottom where
935 ! nearly massless layers layers ride over the topography.
936  if (cs%harmonic_visc) then
937  do i=is,ie ; z_i(i,nz+1) = 0.0 ; enddo
938 
939  do k=nz,1,-1 ; do i=is,ie ; if (do_i(i)) then
940  hvel(i,k) = h_harm(i,k)
941  if (v(i,j,k) * h_delta(i,k) < 0) then
942  z2 = z_i(i,k+1) ; botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
943  hvel(i,k) = (1.0-botfn)*h_harm(i,k) + botfn*h_arith(i,k)
944  endif
945  z_i(i,k) = z_i(i,k+1) + h_harm(i,k)*i_hbbl(i)
946  endif ; enddo ; enddo ! i & k loops
947  else ! Not harmonic_visc
948  do i=is,ie
949  zh(i) = 0.0 ; z_i(i,nz+1) = 0.0
950  zcol1(i) = -g%bathyT(i,j) * gv%Z_to_H
951  zcol2(i) = -g%bathyT(i,j+1) * gv%Z_to_H
952  enddo
953  do k=nz,1,-1 ; do i=is,ie ; if (do_i(i)) then
954  zh(i) = zh(i) + h_harm(i,k)
955  zcol1(i) = zcol1(i) + h(i,j,k) ; zcol2(i) = zcol2(i) + h(i,j+1,k)
956 
957  z_clear = max(zcol1(i),zcol2(i)) + dmin(i)
958  if (zi_dir(i) < 0) z_clear = zcol1(i) + dmin(i)
959  if (zi_dir(i) > 0) z_clear = zcol2(i) + dmin(i)
960 
961  z_i(i,k) = max(zh(i), z_clear) * i_hbbl(i)
962 
963  hvel(i,k) = h_arith(i,k)
964  if (v(i,j,k) * h_delta(i,k) > 0) then
965  if (zh(i) * i_hbbl(i) < cs%harm_BL_val) then
966  hvel(i,k) = h_harm(i,k)
967  else
968  z2_wt = 1.0 ; if (zh(i) * i_hbbl(i) < 2.0*cs%harm_BL_val) &
969  z2_wt = max(0.0, min(1.0, zh(i) * i_hbbl(i) * i_valbl - 1.0))
970  z2 = z2_wt * (max(zh(i), max(zcol1(i),zcol2(i)) + dmin(i)) * i_hbbl(i))
971  botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
972  hvel(i,k) = (1.0-botfn)*h_arith(i,k) + botfn*h_harm(i,k)
973  endif
974  endif
975 
976  endif ; enddo ; enddo ! i & k loops
977  endif
978 
979  call find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_ml, &
980  dt, j, g, gv, us, cs, visc, forces, work_on_u=.false., obc=obc)
981  if ( allocated(hml_v)) then
982  do i=is,ie ; if (do_i(i)) then ; hml_v(i,j) = h_ml(i) ; endif ; enddo
983  endif
984  do_any_shelf = .false.
985  if (associated(forces%frac_shelf_v)) then
986  do i=is,ie
987  cs%a1_shelf_v(i,j) = 0.0
988  do_i_shelf(i) = (do_i(i) .and. forces%frac_shelf_v(i,j) > 0.0)
989  if (do_i_shelf(i)) do_any_shelf = .true.
990  enddo
991  if (do_any_shelf) then
992  if (cs%harmonic_visc) then
993  call find_coupling_coef(a_shelf, hvel, do_i_shelf, h_harm, bbl_thick, &
994  kv_bbl, z_i, h_ml, dt, j, g, gv, us, cs, visc, &
995  forces, work_on_u=.false., obc=obc, shelf=.true.)
996  else ! Find upwind-biased thickness near the surface.
997  ! Perhaps this needs to be done more carefully, via find_eta.
998  do i=is,ie ; if (do_i_shelf(i)) then
999  zh(i) = 0.0 ; ztop_min(i) = min(zcol1(i), zcol2(i))
1000  i_htbl(i) = 1.0 / (visc%tbl_thick_shelf_v(i,j)*gv%Z_to_H + h_neglect)
1001  endif ; enddo
1002  do k=1,nz
1003  do i=is,ie ; if (do_i_shelf(i)) then
1004  zcol1(i) = zcol1(i) - h(i,j,k) ; zcol2(i) = zcol2(i) - h(i,j+1,k)
1005  zh(i) = zh(i) + h_harm(i,k)
1006 
1007  hvel_shelf(i,k) = hvel(i,k)
1008  if (v(i,j,k) * h_delta(i,k) > 0) then
1009  if (zh(i) * i_htbl(i) < cs%harm_BL_val) then
1010  hvel_shelf(i,k) = min(hvel(i,k), h_harm(i,k))
1011  else
1012  z2_wt = 1.0 ; if (zh(i) * i_htbl(i) < 2.0*cs%harm_BL_val) &
1013  z2_wt = max(0.0, min(1.0, zh(i) * i_htbl(i) * i_valbl - 1.0))
1014  z2 = z2_wt * (max(zh(i), ztop_min(i) - min(zcol1(i),zcol2(i))) * i_htbl(i))
1015  topfn = 1.0 / (1.0 + 0.09*z2**6)
1016  hvel_shelf(i,k) = min(hvel(i,k), (1.0-topfn)*h_arith(i,k) + topfn*h_harm(i,k))
1017  endif
1018  endif
1019  endif ; enddo
1020  enddo
1021  call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, h_harm, &
1022  bbl_thick, kv_bbl, z_i, h_ml, dt, j, g, gv, us, cs, &
1023  visc, forces, work_on_u=.false., obc=obc, shelf=.true.)
1024  endif
1025  do i=is,ie ; if (do_i_shelf(i)) cs%a1_shelf_v(i,j) = a_shelf(i,1) ; enddo
1026  endif
1027  endif
1028 
1029  if (do_any_shelf) then
1030  do k=1,nz+1 ; do i=is,ie ; if (do_i_shelf(i)) then
1031  cs%a_v(i,j,k) = min(a_cpl_max, forces%frac_shelf_v(i,j) * a_shelf(i,k) + &
1032  (1.0-forces%frac_shelf_v(i,j)) * a_cpl(i,k))
1033 ! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH
1034 ! CS%a_v(i,J,K) = min(a_cpl_max, forces%frac_shelf_v(i,J) * max(a_shelf(i,K), a_cpl(i,K)) + &
1035  ! (1.0-forces%frac_shelf_v(i,J)) * a_cpl(i,K))
1036  elseif (do_i(i)) then
1037  cs%a_v(i,j,k) = min(a_cpl_max, a_cpl(i,k))
1038  endif ; enddo ; enddo
1039  do k=1,nz ; do i=is,ie ; if (do_i_shelf(i)) then
1040  ! Should we instead take the inverse of the average of the inverses?
1041  cs%h_v(i,j,k) = forces%frac_shelf_v(i,j) * hvel_shelf(i,k) + &
1042  (1.0-forces%frac_shelf_v(i,j)) * hvel(i,k) + h_neglect
1043  elseif (do_i(i)) then
1044  cs%h_v(i,j,k) = hvel(i,k) + h_neglect
1045  endif ; enddo ; enddo
1046  else
1047  do k=1,nz+1 ; do i=is,ie ; if (do_i(i)) cs%a_v(i,j,k) = min(a_cpl_max, a_cpl(i,k)) ; enddo ; enddo
1048  do k=1,nz ; do i=is,ie ; if (do_i(i)) cs%h_v(i,j,k) = hvel(i,k) + h_neglect ; enddo ; enddo
1049  endif
1050 
1051  ! Diagnose total Kv at v-points
1052  if (cs%id_Kv_v > 0) then
1053  do k=1,nz ; do i=is,ie
1054  if (do_i(i)) kv_v(i,j,k) = 0.5 * gv%H_to_Z*(cs%a_v(i,j,k)+cs%a_v(i,j,k+1)) * cs%h_v(i,j,k)
1055  enddo ; enddo
1056  endif
1057 
1058  enddo ! end of v-point j loop
1059 
1060  if (cs%debug) then
1061  call uvchksum("vertvisc_coef h_[uv]", cs%h_u, cs%h_v, g%HI, haloshift=0, &
1062  scale=gv%H_to_m, scalar_pair=.true.)
1063  call uvchksum("vertvisc_coef a_[uv]", cs%a_u, cs%a_v, g%HI, haloshift=0, &
1064  scale=us%Z_to_m*us%s_to_T, scalar_pair=.true.)
1065  if (allocated(hml_u) .and. allocated(hml_v)) &
1066  call uvchksum("vertvisc_coef hML_[uv]", hml_u, hml_v, g%HI, &
1067  haloshift=0, scale=gv%H_to_m, scalar_pair=.true.)
1068  endif
1069 
1070 ! Offer diagnostic fields for averaging.
1071  if (associated(visc%Kv_slow) .and. (cs%id_Kv_slow > 0)) &
1072  call post_data(cs%id_Kv_slow, visc%Kv_slow, cs%diag)
1073  if (cs%id_Kv_u > 0) call post_data(cs%id_Kv_u, kv_u, cs%diag)
1074  if (cs%id_Kv_v > 0) call post_data(cs%id_Kv_v, kv_v, cs%diag)
1075  if (cs%id_au_vv > 0) call post_data(cs%id_au_vv, cs%a_u, cs%diag)
1076  if (cs%id_av_vv > 0) call post_data(cs%id_av_vv, cs%a_v, cs%diag)
1077  if (cs%id_h_u > 0) call post_data(cs%id_h_u, cs%h_u, cs%diag)
1078  if (cs%id_h_v > 0) call post_data(cs%id_h_v, cs%h_v, cs%diag)
1079  if (cs%id_hML_u > 0) call post_data(cs%id_hML_u, hml_u, cs%diag)
1080  if (cs%id_hML_v > 0) call post_data(cs%id_hML_v, hml_v, cs%diag)
1081 
1082  if (allocated(hml_u)) deallocate(hml_u)
1083  if (allocated(hml_v)) deallocate(hml_v)
1084 
1085 end subroutine vertvisc_coef
1086 
1087 !> Calculate the 'coupling coefficient' (a_cpl) at the interfaces.
1088 !! If BOTTOMDRAGLAW is defined, the minimum of Hbbl and half the adjacent
1089 !! layer thicknesses are used to calculate a_cpl near the bottom.
1090 subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_ml, &
1091  dt, j, G, GV, US, CS, visc, forces, work_on_u, OBC, shelf)
1092  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1093  type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure
1094  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1095  real, dimension(SZIB_(G),SZK_(GV)+1), &
1096  intent(out) :: a_cpl !< Coupling coefficient across interfaces [Z T-1 ~> m s-1].
1097  real, dimension(SZIB_(G),SZK_(GV)), &
1098  intent(in) :: hvel !< Thickness at velocity points [H ~> m or kg m-2]
1099  logical, dimension(SZIB_(G)), &
1100  intent(in) :: do_i !< If true, determine coupling coefficient for a column
1101  real, dimension(SZIB_(G),SZK_(GV)), &
1102  intent(in) :: h_harm !< Harmonic mean of thicknesses around a velocity
1103  !! grid point [H ~> m or kg m-2]
1104  real, dimension(SZIB_(G)), intent(in) :: bbl_thick !< Bottom boundary layer thickness [H ~> m or kg m-2]
1105  real, dimension(SZIB_(G)), intent(in) :: kv_bbl !< Bottom boundary layer viscosity [Z2 T-1 ~> m2 s-1].
1106  real, dimension(SZIB_(G),SZK_(GV)+1), &
1107  intent(in) :: z_i !< Estimate of interface heights above the bottom,
1108  !! normalized by the bottom boundary layer thickness
1109  real, dimension(SZIB_(G)), intent(out) :: h_ml !< Mixed layer depth [H ~> m or kg m-2]
1110  integer, intent(in) :: j !< j-index to find coupling coefficient for
1111  real, intent(in) :: dt !< Time increment [T ~> s]
1112  type(vertvisc_cs), pointer :: CS !< Vertical viscosity control structure
1113  type(vertvisc_type), intent(in) :: visc !< Structure containing viscosities and bottom drag
1114  type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
1115  logical, intent(in) :: work_on_u !< If true, u-points are being calculated,
1116  !! otherwise they are v-points
1117  type(ocean_obc_type), pointer :: OBC !< Open boundary condition structure
1118  logical, optional, intent(in) :: shelf !< If present and true, use a surface boundary
1119  !! condition appropriate for an ice shelf.
1120 
1121  ! Local variables
1122 
1123  real, dimension(SZIB_(G)) :: &
1124  u_star, & ! ustar at a velocity point [Z T-1 ~> m s-1].
1125  absf, & ! The average of the neighboring absolute values of f [T-1 ~> s-1].
1126 ! h_ml, & ! The mixed layer depth [H ~> m or kg m-2].
1127  nk_visc, & ! The (real) interface index of the base of mixed layer.
1128  z_t, & ! The distance from the top, sometimes normalized
1129  ! by Hmix, [H ~> m or kg m-2] or [nondim].
1130  kv_tbl, & ! The viscosity in a top boundary layer under ice [Z2 T-1 ~> m2 s-1].
1131  tbl_thick
1132  real, dimension(SZIB_(G),SZK_(GV)+1) :: &
1133  Kv_tot, & ! The total viscosity at an interface [Z2 T-1 ~> m2 s-1].
1134  Kv_add ! A viscosity to add [Z2 T-1 ~> m2 s-1].
1135  real :: h_shear ! The distance over which shears occur [H ~> m or kg m-2].
1136  real :: r ! A thickness to compare with Hbbl [H ~> m or kg m-2].
1137  real :: visc_ml ! The mixed layer viscosity [Z2 T-1 ~> m2 s-1].
1138  real :: I_Hmix ! The inverse of the mixed layer thickness [H-1 ~> m-1 or m2 kg-1].
1139  real :: a_ml ! The layer coupling coefficient across an interface in
1140  ! the mixed layer [Z T-1 ~> m s-1].
1141  real :: I_amax ! The inverse of the maximum coupling coefficient [T Z-1 ~> s m-1].
1142  real :: temp1 ! A temporary variable [H Z ~> m2 or kg m-1]
1143  real :: h_neglect ! A thickness that is so small it is usually lost
1144  ! in roundoff and can be neglected [H ~> m or kg m-2].
1145  real :: z2 ! A copy of z_i [nondim]
1146  real :: botfn ! A function that is 1 at the bottom and small far from it [nondim]
1147  real :: topfn ! A function that is 1 at the top and small far from it [nondim]
1148  real :: kv_top ! A viscosity associated with the top boundary layer [Z2 T-1 ~> m2 s-1]
1149  logical :: do_shelf, do_OBCs
1150  integer :: i, k, is, ie, max_nk
1151  integer :: nz
1152 
1153  a_cpl(:,:) = 0.0
1154  kv_tot(:,:) = 0.0
1155 
1156  if (work_on_u) then ; is = g%IscB ; ie = g%IecB
1157  else ; is = g%isc ; ie = g%iec ; endif
1158  nz = g%ke
1159  h_neglect = gv%H_subroundoff
1160 
1161  if (cs%answers_2018) then
1162  ! The maximum coupling coefficent was originally introduced to avoid
1163  ! truncation error problems in the tridiagonal solver. Effectively, the 1e-10
1164  ! sets the maximum coupling coefficient increment to 1e10 m per timestep.
1165  i_amax = (1.0e-10*us%Z_to_m) * dt
1166  else
1167  i_amax = 0.0
1168  endif
1169 
1170  do_shelf = .false. ; if (present(shelf)) do_shelf = shelf
1171  do_obcs = .false.
1172  if (associated(obc)) then ; do_obcs = (obc%number_of_segments > 0) ; endif
1173  h_ml(:) = 0.0
1174 
1175 ! The following loop calculates the vertical average velocity and
1176 ! surface mixed layer contributions to the vertical viscosity.
1177  do i=is,ie ; kv_tot(i,1) = 0.0 ; enddo
1178  if ((gv%nkml>0) .or. do_shelf) then ; do k=2,nz ; do i=is,ie
1179  if (do_i(i)) kv_tot(i,k) = cs%Kv
1180  enddo ; enddo ; else
1181  i_hmix = 1.0 / (cs%Hmix + h_neglect)
1182  do i=is,ie ; z_t(i) = h_neglect*i_hmix ; enddo
1183  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1184  z_t(i) = z_t(i) + h_harm(i,k-1)*i_hmix
1185  kv_tot(i,k) = cs%Kv + cs%Kvml / ((z_t(i)*z_t(i)) * &
1186  (1.0 + 0.09*z_t(i)*z_t(i)*z_t(i)*z_t(i)*z_t(i)*z_t(i)))
1187  endif ; enddo ; enddo
1188  endif
1189 
1190  do i=is,ie ; if (do_i(i)) then
1191  if (cs%bottomdraglaw) then
1192  r = hvel(i,nz)*0.5
1193  if (r < bbl_thick(i)) then
1194  a_cpl(i,nz+1) = kv_bbl(i) / (i_amax*kv_bbl(i) + (r+h_neglect)*gv%H_to_Z)
1195  else
1196  a_cpl(i,nz+1) = kv_bbl(i) / (i_amax*kv_bbl(i) + (bbl_thick(i)+h_neglect)*gv%H_to_Z)
1197  endif
1198  else
1199  a_cpl(i,nz+1) = cs%Kvbbl / ((0.5*hvel(i,nz)+h_neglect)*gv%H_to_Z + i_amax*cs%Kvbbl)
1200  endif
1201  endif ; enddo
1202 
1203  if (associated(visc%Kv_shear)) then
1204  ! The factor of 2 that used to be required in the viscosities is no longer needed.
1205  if (work_on_u) then
1206  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1207  kv_add(i,k) = 0.5*(visc%Kv_shear(i,j,k) + visc%Kv_shear(i+1,j,k))
1208  endif ; enddo ; enddo
1209  if (do_obcs) then
1210  do i=is,ie ; if (do_i(i) .and. (obc%segnum_u(i,j) /= obc_none)) then
1211  if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
1212  do k=2,nz ; kv_add(i,k) = visc%Kv_shear(i,j,k) ; enddo
1213  elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
1214  do k=2,nz ; kv_add(i,k) = visc%Kv_shear(i+1,j,k) ; enddo
1215  endif
1216  endif ; enddo
1217  endif
1218  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1219  kv_tot(i,k) = kv_tot(i,k) + kv_add(i,k)
1220  endif ; enddo ; enddo
1221  else
1222  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1223  kv_add(i,k) = 0.5*(visc%Kv_shear(i,j,k) + visc%Kv_shear(i,j+1,k))
1224  endif ; enddo ; enddo
1225  if (do_obcs) then
1226  do i=is,ie ; if (do_i(i) .and. (obc%segnum_v(i,j) /= obc_none)) then
1227  if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
1228  do k=2,nz ; kv_add(i,k) = visc%Kv_shear(i,j,k) ; enddo
1229  elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
1230  do k=2,nz ; kv_add(i,k) = visc%Kv_shear(i,j+1,k) ; enddo
1231  endif
1232  endif ; enddo
1233  endif
1234  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1235  kv_tot(i,k) = kv_tot(i,k) + kv_add(i,k)
1236  endif ; enddo ; enddo
1237  endif
1238  endif
1239 
1240  if (associated(visc%Kv_shear_Bu)) then
1241  if (work_on_u) then
1242  do k=2,nz ; do i=is,ie ; If (do_i(i)) then
1243  kv_tot(i,k) = kv_tot(i,k) + (0.5)*(visc%Kv_shear_Bu(i,j-1,k) + visc%Kv_shear_Bu(i,j,k))
1244  endif ; enddo ; enddo
1245  else
1246  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1247  kv_tot(i,k) = kv_tot(i,k) + (0.5)*(visc%Kv_shear_Bu(i-1,j,k) + visc%Kv_shear_Bu(i,j,k))
1248  endif ; enddo ; enddo
1249  endif
1250  endif
1251 
1252  do k=nz,2,-1 ; do i=is,ie ; if (do_i(i)) then
1253  ! botfn determines when a point is within the influence of the bottom
1254  ! boundary layer, going from 1 at the bottom to 0 in the interior.
1255  z2 = z_i(i,k)
1256  botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
1257 
1258  if (cs%bottomdraglaw) then
1259  kv_tot(i,k) = kv_tot(i,k) + (kv_bbl(i) - cs%Kv)*botfn
1260  r = 0.5*(hvel(i,k) + hvel(i,k-1))
1261  if (r > bbl_thick(i)) then
1262  h_shear = ((1.0 - botfn) * r + botfn*bbl_thick(i)) + h_neglect
1263  else
1264  h_shear = r + h_neglect
1265  endif
1266  else
1267  kv_tot(i,k) = kv_tot(i,k) + (cs%Kvbbl-cs%Kv)*botfn
1268  h_shear = 0.5*(hvel(i,k) + hvel(i,k-1) + h_neglect)
1269  endif
1270 
1271  ! Calculate the coupling coefficients from the viscosities.
1272  a_cpl(i,k) = kv_tot(i,k) / (h_shear*gv%H_to_Z + i_amax*kv_tot(i,k))
1273  endif ; enddo ; enddo ! i & k loops
1274 
1275  if (do_shelf) then
1276  ! Set the coefficients to include the no-slip surface stress.
1277  do i=is,ie ; if (do_i(i)) then
1278  if (work_on_u) then
1279  kv_tbl(i) = visc%Kv_tbl_shelf_u(i,j)
1280  tbl_thick(i) = visc%tbl_thick_shelf_u(i,j) * gv%Z_to_H + h_neglect
1281  else
1282  kv_tbl(i) = visc%Kv_tbl_shelf_v(i,j)
1283  tbl_thick(i) = visc%tbl_thick_shelf_v(i,j) * gv%Z_to_H + h_neglect
1284  endif
1285  z_t(i) = 0.0
1286 
1287  ! If a_cpl(i,1) were not already 0, it would be added here.
1288  if (0.5*hvel(i,1) > tbl_thick(i)) then
1289  a_cpl(i,1) = kv_tbl(i) / (tbl_thick(i)*gv%H_to_Z + i_amax*kv_tbl(i))
1290  else
1291  a_cpl(i,1) = kv_tbl(i) / ((0.5*hvel(i,1)+h_neglect)*gv%H_to_Z + i_amax*kv_tbl(i))
1292  endif
1293  endif ; enddo
1294 
1295  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1296  z_t(i) = z_t(i) + hvel(i,k-1) / tbl_thick(i)
1297  topfn = 1.0 / (1.0 + 0.09 * z_t(i)**6)
1298 
1299  r = 0.5*(hvel(i,k)+hvel(i,k-1))
1300  if (r > tbl_thick(i)) then
1301  h_shear = ((1.0 - topfn) * r + topfn*tbl_thick(i)) + h_neglect
1302  else
1303  h_shear = r + h_neglect
1304  endif
1305 
1306  kv_top = topfn * kv_tbl(i)
1307  a_cpl(i,k) = a_cpl(i,k) + kv_top / (h_shear*gv%H_to_Z + i_amax*kv_top)
1308  endif ; enddo ; enddo
1309  elseif (cs%dynamic_viscous_ML .or. (gv%nkml>0)) then
1310  max_nk = 0
1311  do i=is,ie ; if (do_i(i)) then
1312  if (gv%nkml>0) nk_visc(i) = real(gv%nkml+1)
1313  if (work_on_u) then
1314  u_star(i) = 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j))
1315  absf(i) = 0.5*(abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i,j)))
1316  if (cs%dynamic_viscous_ML) nk_visc(i) = visc%nkml_visc_u(i,j) + 1
1317  else
1318  u_star(i) = 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1))
1319  absf(i) = 0.5*(abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j)))
1320  if (cs%dynamic_viscous_ML) nk_visc(i) = visc%nkml_visc_v(i,j) + 1
1321  endif
1322  h_ml(i) = h_neglect ; z_t(i) = 0.0
1323  max_nk = max(max_nk,ceiling(nk_visc(i) - 1.0))
1324  endif ; enddo
1325 
1326  if (do_obcs) then ; if (work_on_u) then
1327  do i=is,ie ; if (do_i(i) .and. (obc%segnum_u(i,j) /= obc_none)) then
1328  if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) &
1329  u_star(i) = forces%ustar(i,j)
1330  if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) &
1331  u_star(i) = forces%ustar(i+1,j)
1332  endif ; enddo
1333  else
1334  do i=is,ie ; if (do_i(i) .and. (obc%segnum_v(i,j) /= obc_none)) then
1335  if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) &
1336  u_star(i) = forces%ustar(i,j)
1337  if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) &
1338  u_star(i) = forces%ustar(i,j+1)
1339  endif ; enddo
1340  endif ; endif
1341 
1342  do k=1,max_nk ; do i=is,ie ; if (do_i(i)) then
1343  if (k+1 <= nk_visc(i)) then ! This layer is all in the ML.
1344  h_ml(i) = h_ml(i) + hvel(i,k)
1345  elseif (k < nk_visc(i)) then ! Part of this layer is in the ML.
1346  h_ml(i) = h_ml(i) + (nk_visc(i) - k) * hvel(i,k)
1347  endif
1348  endif ; enddo ; enddo
1349 
1350  do k=2,max_nk ; do i=is,ie ; if (do_i(i)) then ; if (k < nk_visc(i)) then
1351  ! Set the viscosity at the interfaces.
1352  z_t(i) = z_t(i) + hvel(i,k-1)
1353  temp1 = (z_t(i)*h_ml(i) - z_t(i)*z_t(i))*gv%H_to_Z
1354  ! This viscosity is set to go to 0 at the mixed layer top and bottom (in a log-layer)
1355  ! and be further limited by rotation to give the natural Ekman length.
1356  visc_ml = u_star(i) * 0.41 * (temp1*u_star(i)) / (absf(i)*temp1 + (h_ml(i)+h_neglect)*u_star(i))
1357  a_ml = visc_ml / (0.25*(hvel(i,k)+hvel(i,k-1) + h_neglect) * gv%H_to_Z + 0.5*i_amax*visc_ml)
1358  ! Choose the largest estimate of a.
1359  if (a_ml > a_cpl(i,k)) a_cpl(i,k) = a_ml
1360  endif ; endif ; enddo ; enddo
1361  endif
1362 
1363 end subroutine find_coupling_coef
1364 
1365 !> Velocity components which exceed a threshold for physically reasonable values
1366 !! are truncated. Optionally, any column with excessive velocities may be sent
1367 !! to a diagnostic reporting subroutine.
1368 subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS)
1369  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
1370  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
1371  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1372  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1373  intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]
1374  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1375  intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]
1376  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1377  intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
1378  type(accel_diag_ptrs), intent(in) :: adp !< Acceleration diagnostic pointers
1379  type(cont_diag_ptrs), intent(in) :: cdp !< Continuity diagnostic pointers
1380  type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
1381  type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag
1382  real, intent(in) :: dt !< Time increment [T ~> s]
1383  type(vertvisc_cs), pointer :: cs !< Vertical viscosity control structure
1384 
1385  ! Local variables
1386 
1387  real :: maxvel ! Velocities components greater than maxvel
1388  real :: truncvel ! are truncated to truncvel, both [L T-1 ~> m s-1].
1389  real :: cfl ! The local CFL number.
1390  real :: h_report ! A thickness below which not to report truncations.
1391  real :: dt_rho0 ! The timestep divided by the Boussinesq density [m2 T2 s-1 L-1 Z-1 R-1 ~> s m3 kg-1].
1392  real :: vel_report(szib_(g),szjb_(g)) ! The velocity to report [L T-1 ~> m s-1]
1393  real :: u_old(szib_(g),szj_(g),szk_(g)) ! The previous u-velocity [L T-1 ~> m s-1]
1394  real :: v_old(szi_(g),szjb_(g),szk_(g)) ! The previous v-velocity [L T-1 ~> m s-1]
1395  logical :: trunc_any, dowrite(szib_(g),szjb_(g))
1396  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
1397  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1398  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1399 
1400  maxvel = cs%maxvel
1401  truncvel = 0.9*maxvel
1402  h_report = 6.0 * gv%Angstrom_H
1403  dt_rho0 = (us%L_T_to_m_s*us%Z_to_m) * dt / gv%Rho0
1404 
1405  if (len_trim(cs%u_trunc_file) > 0) then
1406  !$OMP parallel do default(shared) private(trunc_any,CFL)
1407  do j=js,je
1408  trunc_any = .false.
1409  do i=isq,ieq ; dowrite(i,j) = .false. ; enddo
1410  if (cs%CFL_based_trunc) then
1411  do i=isq,ieq ; vel_report(i,j) = 3.0e8*us%m_s_to_L_T ; enddo ! Speed of light default.
1412  do k=1,nz ; do i=isq,ieq
1413  if (abs(u(i,j,k)) < cs%vel_underflow) u(i,j,k) = 0.0
1414  if (u(i,j,k) < 0.0) then
1415  cfl = (-u(i,j,k) * dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
1416  else
1417  cfl = (u(i,j,k) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
1418  endif
1419  if (cfl > cs%CFL_trunc) trunc_any = .true.
1420  if (cfl > cs%CFL_report) then
1421  dowrite(i,j) = .true.
1422  vel_report(i,j) = min(vel_report(i,j), abs(u(i,j,k)))
1423  endif
1424  enddo ; enddo
1425  else
1426  do i=isq,ieq; vel_report(i,j) = maxvel; enddo
1427  do k=1,nz ; do i=isq,ieq
1428  if (abs(u(i,j,k)) < cs%vel_underflow) then ; u(i,j,k) = 0.0
1429  elseif (abs(u(i,j,k)) > maxvel) then
1430  dowrite(i,j) = .true. ; trunc_any = .true.
1431  endif
1432  enddo ; enddo
1433  endif
1434 
1435  do i=isq,ieq ; if (dowrite(i,j)) then
1436  u_old(i,j,:) = u(i,j,:)
1437  endif ; enddo
1438 
1439  if (trunc_any) then ; if (cs%CFL_based_trunc) then
1440  do k=1,nz ; do i=isq,ieq
1441  if ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i+1,j) < -cs%CFL_trunc) then
1442  u(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i+1,j) / (dt * g%dy_Cu(i,j)))
1443  if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1444  elseif ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then
1445  u(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dy_Cu(i,j)))
1446  if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1447  endif
1448  enddo ; enddo
1449  else
1450  do k=1,nz ; do i=isq,ieq ; if (abs(u(i,j,k)) > maxvel) then
1451  u(i,j,k) = sign(truncvel,u(i,j,k))
1452  if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1453  endif ; enddo ; enddo
1454  endif ; endif
1455  enddo ! j-loop
1456  else ! Do not report accelerations leading to large velocities.
1457  if (cs%CFL_based_trunc) then
1458 !$OMP parallel do default(none) shared(nz,js,je,Isq,Ieq,u,dt,G,CS,h,H_report)
1459  do k=1,nz ; do j=js,je ; do i=isq,ieq
1460  if (abs(u(i,j,k)) < cs%vel_underflow) then ; u(i,j,k) = 0.0
1461  elseif ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i+1,j) < -cs%CFL_trunc) then
1462  u(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i+1,j) / (dt * g%dy_Cu(i,j)))
1463  if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1464  elseif ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then
1465  u(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dy_Cu(i,j)))
1466  if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1467  endif
1468  enddo ; enddo ; enddo
1469  else
1470 !$OMP parallel do default(none) shared(nz,js,je,Isq,Ieq,u,G,CS,truncvel,maxvel,h,H_report)
1471  do k=1,nz ; do j=js,je ; do i=isq,ieq
1472  if (abs(u(i,j,k)) < cs%vel_underflow) then ; u(i,j,k) = 0.0
1473  elseif (abs(u(i,j,k)) > maxvel) then
1474  u(i,j,k) = sign(truncvel, u(i,j,k))
1475  if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1476  endif
1477  enddo ; enddo ; enddo
1478  endif
1479  endif
1480 
1481  if (len_trim(cs%u_trunc_file) > 0) then
1482  do j=js,je; do i=isq,ieq ; if (dowrite(i,j)) then
1483 ! Here the diagnostic reporting subroutines are called if
1484 ! unphysically large values were found.
1485  call write_u_accel(i, j, u_old, h, adp, cdp, dt, g, gv, us, cs%PointAccel_CSp, &
1486  vel_report(i,j), forces%taux(i,j)*dt_rho0, a=cs%a_u, hv=cs%h_u)
1487  endif ; enddo ; enddo
1488  endif
1489 
1490  if (len_trim(cs%v_trunc_file) > 0) then
1491  !$OMP parallel do default(shared) private(trunc_any,CFL)
1492  do j=jsq,jeq
1493  trunc_any = .false.
1494  do i=is,ie ; dowrite(i,j) = .false. ; enddo
1495  if (cs%CFL_based_trunc) then
1496  do i=is,ie ; vel_report(i,j) = 3.0e8*us%m_s_to_L_T ; enddo ! Speed of light default.
1497  do k=1,nz ; do i=is,ie
1498  if (abs(v(i,j,k)) < cs%vel_underflow) v(i,j,k) = 0.0
1499  if (v(i,j,k) < 0.0) then
1500  cfl = (-v(i,j,k) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
1501  else
1502  cfl = (v(i,j,k) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
1503  endif
1504  if (cfl > cs%CFL_trunc) trunc_any = .true.
1505  if (cfl > cs%CFL_report) then
1506  dowrite(i,j) = .true.
1507  vel_report(i,j) = min(vel_report(i,j), abs(v(i,j,k)))
1508  endif
1509  enddo ; enddo
1510  else
1511  do i=is,ie ; vel_report(i,j) = maxvel ; enddo
1512  do k=1,nz ; do i=is,ie
1513  if (abs(v(i,j,k)) < cs%vel_underflow) then ; v(i,j,k) = 0.0
1514  elseif (abs(v(i,j,k)) > maxvel) then
1515  dowrite(i,j) = .true. ; trunc_any = .true.
1516  endif
1517  enddo ; enddo
1518  endif
1519 
1520  do i=is,ie ; if (dowrite(i,j)) then
1521  v_old(i,j,:) = v(i,j,:)
1522  endif ; enddo
1523 
1524  if (trunc_any) then ; if (cs%CFL_based_trunc) then
1525  do k=1,nz; do i=is,ie
1526  if ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j+1) < -cs%CFL_trunc) then
1527  v(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i,j+1) / (dt * g%dx_Cv(i,j)))
1528  if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1529  elseif ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then
1530  v(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dx_Cv(i,j)))
1531  if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1532  endif
1533  enddo ; enddo
1534  else
1535  do k=1,nz ; do i=is,ie ; if (abs(v(i,j,k)) > maxvel) then
1536  v(i,j,k) = sign(truncvel,v(i,j,k))
1537  if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1538  endif ; enddo ; enddo
1539  endif ; endif
1540  enddo ! J-loop
1541  else ! Do not report accelerations leading to large velocities.
1542  if (cs%CFL_based_trunc) then
1543  !$OMP parallel do default(shared)
1544  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
1545  if (abs(v(i,j,k)) < cs%vel_underflow) then ; v(i,j,k) = 0.0
1546  elseif ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j+1) < -cs%CFL_trunc) then
1547  v(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i,j+1) / (dt * g%dx_Cv(i,j)))
1548  if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1549  elseif ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then
1550  v(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dx_Cv(i,j)))
1551  if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1552  endif
1553  enddo ; enddo ; enddo
1554  else
1555  !$OMP parallel do default(shared)
1556  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
1557  if (abs(v(i,j,k)) < cs%vel_underflow) then ; v(i,j,k) = 0.0
1558  elseif (abs(v(i,j,k)) > maxvel) then
1559  v(i,j,k) = sign(truncvel, v(i,j,k))
1560  if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1561  endif
1562  enddo ; enddo ; enddo
1563  endif
1564  endif
1565 
1566  if (len_trim(cs%v_trunc_file) > 0) then
1567  do j=jsq,jeq; do i=is,ie ; if (dowrite(i,j)) then
1568 ! Here the diagnostic reporting subroutines are called if
1569 ! unphysically large values were found.
1570  call write_v_accel(i, j, v_old, h, adp, cdp, dt, g, gv, us, cs%PointAccel_CSp, &
1571  vel_report(i,j), forces%tauy(i,j)*dt_rho0, a=cs%a_v, hv=cs%h_v)
1572  endif ; enddo ; enddo
1573  endif
1574 
1575 end subroutine vertvisc_limit_vel
1576 
1577 !> Initialize the vertical friction module
1578 subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
1579  ntrunc, CS)
1580  type(ocean_internal_state), &
1581  target, intent(in) :: mis !< The "MOM Internal State", a set of pointers
1582  !! to the fields and accelerations that make
1583  !! up the ocean's physical state
1584  type(time_type), target, intent(in) :: time !< Current model time
1585  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
1586  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
1587  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1588  type(param_file_type), intent(in) :: param_file !< File to parse for parameters
1589  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostic control structure
1590  type(accel_diag_ptrs), intent(inout) :: adp !< Acceleration diagnostic pointers
1591  type(directories), intent(in) :: dirs !< Relevant directory paths
1592  integer, target, intent(inout) :: ntrunc !< Number of velocity truncations
1593  type(vertvisc_cs), pointer :: cs !< Vertical viscosity control structure
1594 
1595  ! Local variables
1596 
1597  real :: hmix_str_dflt
1598  real :: kv_dflt ! A default viscosity [m2 s-1].
1599  real :: hmix_m ! A boundary layer thickness [m].
1600  logical :: default_2018_answers
1601  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz
1602 ! This include declares and sets the variable "version".
1603 #include "version_variable.h"
1604  character(len=40) :: mdl = "MOM_vert_friction" ! This module's name.
1605  character(len=40) :: thickness_units
1606 
1607  if (associated(cs)) then
1608  call mom_error(warning, "vertvisc_init called with an associated "// &
1609  "control structure.")
1610  return
1611  endif
1612  allocate(cs)
1613 
1614  if (gv%Boussinesq) then; thickness_units = "m"
1615  else; thickness_units = "kg m-2"; endif
1616 
1617  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
1618  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1619 
1620  cs%diag => diag ; cs%ntrunc => ntrunc ; ntrunc = 0
1621 
1622 ! Default, read and log parameters
1623  call log_version(param_file, mdl, version, "", log_to_all=.true., debugging=.true.)
1624  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
1625  "This sets the default value for the various _2018_ANSWERS parameters.", &
1626  default=.false.)
1627  call get_param(param_file, mdl, "VERT_FRICTION_2018_ANSWERS", cs%answers_2018, &
1628  "If true, use the order of arithmetic and expressions that recover the answers "//&
1629  "from the end of 2018. Otherwise, use expressions that do not use an arbitrary "//&
1630  "hard-coded maximum viscous coupling coefficient between layers.", &
1631  default=default_2018_answers)
1632  call get_param(param_file, mdl, "BOTTOMDRAGLAW", cs%bottomdraglaw, &
1633  "If true, the bottom stress is calculated with a drag "//&
1634  "law of the form c_drag*|u|*u. The velocity magnitude "//&
1635  "may be an assumed value or it may be based on the "//&
1636  "actual velocity in the bottommost HBBL, depending on "//&
1637  "LINEAR_DRAG.", default=.true.)
1638  call get_param(param_file, mdl, "CHANNEL_DRAG", cs%Channel_drag, &
1639  "If true, the bottom drag is exerted directly on each "//&
1640  "layer proportional to the fraction of the bottom it "//&
1641  "overlies.", default=.false.)
1642  call get_param(param_file, mdl, "DIRECT_STRESS", cs%direct_stress, &
1643  "If true, the wind stress is distributed over the "//&
1644  "topmost HMIX_STRESS of fluid (like in HYCOM), and KVML "//&
1645  "may be set to a very small value.", default=.false.)
1646  call get_param(param_file, mdl, "DYNAMIC_VISCOUS_ML", cs%dynamic_viscous_ML, &
1647  "If true, use a bulk Richardson number criterion to "//&
1648  "determine the mixed layer thickness for viscosity.", &
1649  default=.false.)
1650  call get_param(param_file, mdl, "U_TRUNC_FILE", cs%u_trunc_file, &
1651  "The absolute path to a file into which the accelerations "//&
1652  "leading to zonal velocity truncations are written. "//&
1653  "Undefine this for efficiency if this diagnostic is not "//&
1654  "needed.", default=" ", debuggingparam=.true.)
1655  call get_param(param_file, mdl, "V_TRUNC_FILE", cs%v_trunc_file, &
1656  "The absolute path to a file into which the accelerations "//&
1657  "leading to meridional velocity truncations are written. "//&
1658  "Undefine this for efficiency if this diagnostic is not "//&
1659  "needed.", default=" ", debuggingparam=.true.)
1660  call get_param(param_file, mdl, "HARMONIC_VISC", cs%harmonic_visc, &
1661  "If true, use the harmonic mean thicknesses for "//&
1662  "calculating the vertical viscosity.", default=.false.)
1663  call get_param(param_file, mdl, "HARMONIC_BL_SCALE", cs%harm_BL_val, &
1664  "A scale to determine when water is in the boundary "//&
1665  "layers based solely on harmonic mean thicknesses for "//&
1666  "the purpose of determining the extent to which the "//&
1667  "thicknesses used in the viscosities are upwinded.", &
1668  default=0.0, units="nondim")
1669  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
1670 
1671  if (gv%nkml < 1) &
1672  call get_param(param_file, mdl, "HMIX_FIXED", cs%Hmix, &
1673  "The prescribed depth over which the near-surface "//&
1674  "viscosity and diffusivity are elevated when the bulk "//&
1675  "mixed layer is not used.", units="m", scale=gv%m_to_H, &
1676  unscaled=hmix_m, fail_if_missing=.true.)
1677  if (cs%direct_stress) then
1678  if (gv%nkml < 1) then
1679  call get_param(param_file, mdl, "HMIX_STRESS", cs%Hmix_stress, &
1680  "The depth over which the wind stress is applied if "//&
1681  "DIRECT_STRESS is true.", units="m", default=hmix_m, scale=gv%m_to_H)
1682  else
1683  call get_param(param_file, mdl, "HMIX_STRESS", cs%Hmix_stress, &
1684  "The depth over which the wind stress is applied if "//&
1685  "DIRECT_STRESS is true.", units="m", fail_if_missing=.true., scale=gv%m_to_H)
1686  endif
1687  if (cs%Hmix_stress <= 0.0) call mom_error(fatal, "vertvisc_init: " // &
1688  "HMIX_STRESS must be set to a positive value if DIRECT_STRESS is true.")
1689  endif
1690  call get_param(param_file, mdl, "KV", cs%Kv, &
1691  "The background kinematic viscosity in the interior. "//&
1692  "The molecular value, ~1e-6 m2 s-1, may be used.", &
1693  units="m2 s-1", fail_if_missing=.true., scale=us%m2_s_to_Z2_T, unscaled=kv_dflt)
1694 
1695  if (gv%nkml < 1) call get_param(param_file, mdl, "KVML", cs%Kvml, &
1696  "The kinematic viscosity in the mixed layer. A typical "//&
1697  "value is ~1e-2 m2 s-1. KVML is not used if "//&
1698  "BULKMIXEDLAYER is true. The default is set by KV.", &
1699  units="m2 s-1", default=kv_dflt, scale=us%m2_s_to_Z2_T)
1700  if (.not.cs%bottomdraglaw) call get_param(param_file, mdl, "KVBBL", cs%Kvbbl, &
1701  "The kinematic viscosity in the benthic boundary layer. "//&
1702  "A typical value is ~1e-2 m2 s-1. KVBBL is not used if "//&
1703  "BOTTOMDRAGLAW is true. The default is set by KV.", &
1704  units="m2 s-1", default=kv_dflt, scale=us%m2_s_to_Z2_T)
1705  call get_param(param_file, mdl, "HBBL", cs%Hbbl, &
1706  "The thickness of a bottom boundary layer with a "//&
1707  "viscosity of KVBBL if BOTTOMDRAGLAW is not defined, or "//&
1708  "the thickness over which near-bottom velocities are "//&
1709  "averaged for the drag law if BOTTOMDRAGLAW is defined "//&
1710  "but LINEAR_DRAG is not.", units="m", fail_if_missing=.true., scale=gv%m_to_H)
1711  call get_param(param_file, mdl, "MAXVEL", cs%maxvel, &
1712  "The maximum velocity allowed before the velocity "//&
1713  "components are truncated.", units="m s-1", default=3.0e8, scale=us%m_s_to_L_T)
1714  call get_param(param_file, mdl, "CFL_BASED_TRUNCATIONS", cs%CFL_based_trunc, &
1715  "If true, base truncations on the CFL number, and not an "//&
1716  "absolute speed.", default=.true.)
1717  call get_param(param_file, mdl, "CFL_TRUNCATE", cs%CFL_trunc, &
1718  "The value of the CFL number that will cause velocity "//&
1719  "components to be truncated; instability can occur past 0.5.", &
1720  units="nondim", default=0.5)
1721  call get_param(param_file, mdl, "CFL_REPORT", cs%CFL_report, &
1722  "The value of the CFL number that causes accelerations "//&
1723  "to be reported; the default is CFL_TRUNCATE.", &
1724  units="nondim", default=cs%CFL_trunc)
1725  call get_param(param_file, mdl, "CFL_TRUNCATE_RAMP_TIME", cs%truncRampTime, &
1726  "The time over which the CFL truncation value is ramped "//&
1727  "up at the beginning of the run.", &
1728  units="s", default=0.)
1729  cs%CFL_truncE = cs%CFL_trunc
1730  call get_param(param_file, mdl, "CFL_TRUNCATE_START", cs%CFL_truncS, &
1731  "The start value of the truncation CFL number used when "//&
1732  "ramping up CFL_TRUNC.", &
1733  units="nondim", default=0.)
1734  call get_param(param_file, mdl, "STOKES_MIXING_COMBINED", cs%StokesMixing, &
1735  "Flag to use Stokes drift Mixing via the Lagrangian "//&
1736  " current (Eulerian plus Stokes drift). "//&
1737  " Still needs work and testing, so not recommended for use.",&
1738  default=.false.)
1739  !BGR 04/04/2018{
1740  ! StokesMixing is required for MOM6 for some Langmuir mixing parameterization.
1741  ! The code used here has not been developed for vanishing layers or in
1742  ! conjunction with any bottom friction. Therefore, the following line is
1743  ! added so this functionality cannot be used without user intervention in
1744  ! the code. This will prevent general use of this functionality until proper
1745  ! care is given to the previously mentioned issues. Comment out the following
1746  ! MOM_error to use, but do so at your own risk and with these points in mind.
1747  !}
1748  if (cs%StokesMixing) then
1749  call mom_error(fatal, "Stokes mixing requires user intervention in the code.\n"//&
1750  " Model now exiting. See MOM_vert_friction.F90 for \n"//&
1751  " details (search 'BGR 04/04/2018' to locate comment).")
1752  endif
1753  call get_param(param_file, mdl, "VEL_UNDERFLOW", cs%vel_underflow, &
1754  "A negligibly small velocity magnitude below which velocity "//&
1755  "components are set to 0. A reasonable value might be "//&
1756  "1e-30 m/s, which is less than an Angstrom divided by "//&
1757  "the age of the universe.", units="m s-1", default=0.0, scale=us%m_s_to_L_T)
1758 
1759  alloc_(cs%a_u(isdb:iedb,jsd:jed,nz+1)) ; cs%a_u(:,:,:) = 0.0
1760  alloc_(cs%h_u(isdb:iedb,jsd:jed,nz)) ; cs%h_u(:,:,:) = 0.0
1761  alloc_(cs%a_v(isd:ied,jsdb:jedb,nz+1)) ; cs%a_v(:,:,:) = 0.0
1762  alloc_(cs%h_v(isd:ied,jsdb:jedb,nz)) ; cs%h_v(:,:,:) = 0.0
1763 
1764  cs%id_Kv_slow = register_diag_field('ocean_model', 'Kv_slow', diag%axesTi, time, &
1765  'Slow varying vertical viscosity', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
1766 
1767  cs%id_Kv_u = register_diag_field('ocean_model', 'Kv_u', diag%axesCuL, time, &
1768  'Total vertical viscosity at u-points', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
1769 
1770  cs%id_Kv_v = register_diag_field('ocean_model', 'Kv_v', diag%axesCvL, time, &
1771  'Total vertical viscosity at v-points', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
1772 
1773  cs%id_au_vv = register_diag_field('ocean_model', 'au_visc', diag%axesCui, time, &
1774  'Zonal Viscous Vertical Coupling Coefficient', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
1775 
1776  cs%id_av_vv = register_diag_field('ocean_model', 'av_visc', diag%axesCvi, time, &
1777  'Meridional Viscous Vertical Coupling Coefficient', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
1778 
1779  cs%id_h_u = register_diag_field('ocean_model', 'Hu_visc', diag%axesCuL, time, &
1780  'Thickness at Zonal Velocity Points for Viscosity', thickness_units, &
1781  conversion=gv%H_to_m)
1782 
1783  cs%id_h_v = register_diag_field('ocean_model', 'Hv_visc', diag%axesCvL, time, &
1784  'Thickness at Meridional Velocity Points for Viscosity', thickness_units, &
1785  conversion=gv%H_to_m)
1786 
1787  cs%id_hML_u = register_diag_field('ocean_model', 'HMLu_visc', diag%axesCu1, time, &
1788  'Mixed Layer Thickness at Zonal Velocity Points for Viscosity', thickness_units, &
1789  conversion=gv%H_to_m)
1790 
1791  cs%id_hML_v = register_diag_field('ocean_model', 'HMLv_visc', diag%axesCv1, time, &
1792  'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', thickness_units, &
1793  conversion=gv%H_to_m)
1794 
1795  cs%id_du_dt_visc = register_diag_field('ocean_model', 'du_dt_visc', diag%axesCuL, &
1796  time, 'Zonal Acceleration from Vertical Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)
1797  if (cs%id_du_dt_visc > 0) call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)
1798  cs%id_dv_dt_visc = register_diag_field('ocean_model', 'dv_dt_visc', diag%axesCvL, &
1799  time, 'Meridional Acceleration from Vertical Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)
1800  if (cs%id_dv_dt_visc > 0) call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)
1801 
1802  cs%id_taux_bot = register_diag_field('ocean_model', 'taux_bot', diag%axesCu1, &
1803  time, 'Zonal Bottom Stress from Ocean to Earth', 'Pa', &
1804  conversion=us%RZ_to_kg_m2*us%L_T2_to_m_s2)
1805  cs%id_tauy_bot = register_diag_field('ocean_model', 'tauy_bot', diag%axesCv1, &
1806  time, 'Meridional Bottom Stress from Ocean to Earth', 'Pa', &
1807  conversion=us%RZ_to_kg_m2*us%L_T2_to_m_s2)
1808 
1809  !CS%id_hf_du_dt_visc = register_diag_field('ocean_model', 'hf_du_dt_visc', diag%axesCuL, Time, &
1810  ! 'Fractional Thickness-weighted Zonal Acceleration from Vertical Viscosity', 'm s-2', &
1811  ! v_extensive=.true., conversion=US%L_T2_to_m_s2)
1812  !if (CS%id_hf_du_dt_visc > 0) then
1813  ! call safe_alloc_ptr(CS%hf_du_dt_visc,IsdB,IedB,jsd,jed,nz)
1814  ! call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
1815  ! call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1816  !endif
1817 
1818  !CS%id_hf_dv_dt_visc = register_diag_field('ocean_model', 'hf_dv_dt_visc', diag%axesCvL, Time, &
1819  ! 'Fractional Thickness-weighted Meridional Acceleration from Vertical Viscosity', 'm s-2', &
1820  ! v_extensive=.true., conversion=US%L_T2_to_m_s2)
1821  !if (CS%id_hf_dv_dt_visc > 0) then
1822  ! call safe_alloc_ptr(CS%hf_dv_dt_visc,isd,ied,JsdB,JedB,nz)
1823  ! call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
1824  ! call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
1825  !endif
1826 
1827  cs%id_hf_du_dt_visc_2d = register_diag_field('ocean_model', 'hf_du_dt_visc_2d', diag%axesCu1, time, &
1828  'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Vertical Viscosity', 'm s-2', &
1829  conversion=us%L_T2_to_m_s2)
1830  if (cs%id_hf_du_dt_visc_2d > 0) then
1831  call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)
1832  call safe_alloc_ptr(adp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
1833  endif
1834 
1835  cs%id_hf_dv_dt_visc_2d = register_diag_field('ocean_model', 'hf_dv_dt_visc_2d', diag%axesCv1, time, &
1836  'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Vertical Viscosity', 'm s-2', &
1837  conversion=us%L_T2_to_m_s2)
1838  if (cs%id_hf_dv_dt_visc_2d > 0) then
1839  call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)
1840  call safe_alloc_ptr(adp%diag_hfrac_v,isd,ied,jsd,jedb,nz)
1841  endif
1842 
1843  if ((len_trim(cs%u_trunc_file) > 0) .or. (len_trim(cs%v_trunc_file) > 0)) &
1844  call pointaccel_init(mis, time, g, param_file, diag, dirs, cs%PointAccel_CSp)
1845 
1846 end subroutine vertvisc_init
1847 
1848 !> Update the CFL truncation value as a function of time.
1849 !! If called with the optional argument activate=.true., record the
1850 !! value of Time as the beginning of the ramp period.
1851 subroutine updatecfltruncationvalue(Time, CS, activate)
1852  type(time_type), target, intent(in) :: time !< Current model time
1853  type(vertvisc_cs), pointer :: cs !< Vertical viscosity control structure
1854  logical, optional, intent(in) :: activate !< Specifiy whether to record the value of
1855  !! Time as the beginning of the ramp period
1856 
1857  ! Local variables
1858  real :: deltatime, wghta
1859  character(len=12) :: msg
1860 
1861  if (cs%truncRampTime==0.) return ! This indicates to ramping is turned off
1862 
1863  ! We use the optional argument to indicate this Time should be recorded as the
1864  ! beginning of the ramp-up period.
1865  if (present(activate)) then
1866  if (activate) then
1867  cs%rampStartTime = time ! Record the current time
1868  cs%CFLrampingIsActivated = .true.
1869  endif
1870  endif
1871  if (.not.cs%CFLrampingIsActivated) return
1872  deltatime = max( 0., time_type_to_real( time - cs%rampStartTime ) )
1873  if (deltatime >= cs%truncRampTime) then
1874  cs%CFL_trunc = cs%CFL_truncE
1875  cs%truncRampTime = 0. ! This turns off ramping after this call
1876  else
1877  wghta = min( 1., deltatime / cs%truncRampTime ) ! Linear profile in time
1878  !wghtA = wghtA*wghtA ! Convert linear profile to parabolic profile in time
1879  !wghtA = wghtA*wghtA*(3. - 2.*wghtA) ! Convert linear profile to cosine profile
1880  wghta = 1. - ( (1. - wghta)**2 ) ! Convert linear profiel to nverted parabolic profile
1881  cs%CFL_trunc = cs%CFL_truncS + wghta * ( cs%CFL_truncE - cs%CFL_truncS )
1882  endif
1883  write(msg(1:12),'(es12.3)') cs%CFL_trunc
1884  call mom_error(note, "MOM_vert_friction: updateCFLtruncationValue set CFL"// &
1885  " limit to "//trim(msg))
1886 end subroutine updatecfltruncationvalue
1887 
1888 !> Clean up and deallocate the vertical friction module
1889 subroutine vertvisc_end(CS)
1890  type(vertvisc_cs), pointer :: cs !< Vertical viscosity control structure that
1891  !! will be deallocated in this subroutine.
1892 
1893  dealloc_(cs%a_u) ; dealloc_(cs%h_u)
1894  dealloc_(cs%a_v) ; dealloc_(cs%h_v)
1895  if (associated(cs%a1_shelf_u)) deallocate(cs%a1_shelf_u)
1896  if (associated(cs%a1_shelf_v)) deallocate(cs%a1_shelf_v)
1897  deallocate(cs)
1898 end subroutine vertvisc_end
1899 
1900 !> \namespace mom_vert_friction
1901 !! \author Robert Hallberg
1902 !! \date April 1994 - October 2006
1903 !!
1904 !! The vertical diffusion of momentum is fully implicit. This is
1905 !! necessary to allow for vanishingly small layers. The coupling
1906 !! is based on the distance between the centers of adjacent layers,
1907 !! except where a layer is close to the bottom compared with a
1908 !! bottom boundary layer thickness when a bottom drag law is used.
1909 !! A stress top b.c. and a no slip bottom b.c. are used. There
1910 !! is no limit on the time step for vertvisc.
1911 !!
1912 !! Near the bottom, the horizontal thickness interpolation scheme
1913 !! changes to an upwind biased estimate to control the effect of
1914 !! spurious Montgomery potential gradients at the bottom where
1915 !! nearly massless layers layers ride over the topography. Within a
1916 !! few boundary layer depths of the bottom, the harmonic mean
1917 !! thickness (i.e. (2 h+ h-) / (h+ + h-) ) is used if the velocity
1918 !! is from the thinner side and the arithmetic mean thickness
1919 !! (i.e. (h+ + h-)/2) is used if the velocity is from the thicker
1920 !! side. Both of these thickness estimates are second order
1921 !! accurate. Above this the arithmetic mean thickness is used.
1922 !!
1923 !! In addition, vertvisc truncates any velocity component that
1924 !! exceeds maxvel to truncvel. This basically keeps instabilities
1925 !! spatially localized. The number of times the velocity is
1926 !! truncated is reported each time the energies are saved, and if
1927 !! exceeds CS%Maxtrunc the model will stop itself and change the time
1928 !! to a large value. This has proven very useful in (1) diagnosing
1929 !! model failures and (2) letting the model settle down to a
1930 !! meaningful integration from a poorly specified initial condition.
1931 !!
1932 !! The same code is used for the two velocity components, by
1933 !! indirectly referencing the velocities and defining a handful of
1934 !! direction-specific defined variables.
1935 !!
1936 !! Macros written all in capital letters are defined in MOM_memory.h.
1937 !!
1938 !! A small fragment of the grid is shown below:
1939 !! \verbatim
1940 !! j+1 x ^ x ^ x At x: q
1941 !! j+1 > o > o > At ^: v, frhatv, tauy
1942 !! j x ^ x ^ x At >: u, frhatu, taux
1943 !! j > o > o > At o: h
1944 !! j-1 x ^ x ^ x
1945 !! i-1 i i+1 At x & ^:
1946 !! i i+1 At > & o:
1947 !! \endverbatim
1948 !!
1949 !! The boundaries always run through q grid points (x).
1950 end module mom_vert_friction
mom_time_manager
Wraps the FMS time manager functions.
Definition: MOM_time_manager.F90:2
mom_forcing_type::mech_forcing
Structure that contains pointers to the mechanical forcing at the surface used to drive the liquid oc...
Definition: MOM_forcing_type.F90:204
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:80
mom_get_input::directories
Container for paths and parameter file names.
Definition: MOM_get_input.F90:20
mom_domains::pass_var
Do a halo update on an array.
Definition: MOM_domains.F90:54
mom_pointaccel::pointaccel_cs
The control structure for the MOM_PointAccel module.
Definition: MOM_PointAccel.F90:32
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_vert_friction::vertvisc_cs
The control structure with parameters and memory for the MOM_vert_friction module.
Definition: MOM_vert_friction.F90:39
mom_get_input
Reads the only Fortran name list needed to boot-strap the model.
Definition: MOM_get_input.F90:6
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_forcing_type
This module implements boundary forcing for MOM6.
Definition: MOM_forcing_type.F90:2
mom_wave_interface
Interface for surface waves.
Definition: MOM_wave_interface.F90:2
mom_wave_interface::wave_parameters_cs
Container for all surface wave related parameters.
Definition: MOM_wave_interface.F90:47
mom_variables::vertvisc_type
Vertical viscosities, drag coefficients, and related fields.
Definition: MOM_variables.F90:208
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_open_boundary
Controls where open boundary conditions are applied.
Definition: MOM_open_boundary.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_variables::cont_diag_ptrs
Pointers to arrays with transports, which can later be used for derived diagnostics,...
Definition: MOM_variables.F90:193
mom_variables::accel_diag_ptrs
Pointers to arrays with accelerations, which can later be used for derived diagnostics,...
Definition: MOM_variables.F90:161
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_open_boundary::ocean_obc_type
Open-boundary data.
Definition: MOM_open_boundary.F90:218
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_debugging
Provides checksumming functions for debugging.
Definition: MOM_debugging.F90:7
mom_vert_friction
Implements vertical viscosity (vertvisc)
Definition: MOM_vert_friction.F90:2
mom_variables::ocean_internal_state
Pointers to all of the prognostic variables allocated in MOM_variables.F90 and MOM....
Definition: MOM_variables.F90:132
mom_pointaccel
Debug accelerations at a given point.
Definition: MOM_PointAccel.F90:8
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:240