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