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