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