mom_vert_friction Module Reference

Detailed Description

Implements vertical viscosity (vertvisc)

Author

Robert Hallberg

Date

April 1994 - October 2006

The vertical diffusion of momentum is fully implicit. This is necessary to allow for vanishingly small layers. The coupling is based on the distance between the centers of adjacent layers, except where a layer is close to the bottom compared with a bottom boundary layer thickness when a bottom drag law is used. A stress top b.c. and a no slip bottom b.c. are used. There is no limit on the time step for vertvisc.

Near the bottom, the horizontal thickness interpolation scheme changes to an upwind biased estimate to control the effect of spurious Montgomery potential gradients at the bottom where nearly massless layers layers ride over the topography. Within a few boundary layer depths of the bottom, the harmonic mean thickness (i.e. (2 h+ h-) / (h+ + h-) ) is used if the velocity is from the thinner side and the arithmetic mean thickness (i.e. (h+ + h-)/2) is used if the velocity is from the thicker side. Both of these thickness estimates are second order accurate. Above this the arithmetic mean thickness is used.

In addition, vertvisc truncates any velocity component that exceeds maxvel to truncvel. This basically keeps instabilities spatially localized. The number of times the velocity is truncated is reported each time the energies are saved, and if exceeds CSMaxtrunc the model will stop itself and change the time to a large value. This has proven very useful in (1) diagnosing model failures and (2) letting the model settle down to a meaningful integration from a poorly specified initial condition.

The same code is used for the two velocity components, by indirectly referencing the velocities and defining a handful of direction-specific defined variables.

Macros written all in capital letters are defined in MOM_memory.h.

A small fragment of the grid is shown below:

    j+1  x ^ x ^ x   At x:  q
    j+1  > o > o >   At ^:  v, frhatv, tauy
    j    x ^ x ^ x   At >:  u, frhatu, taux
    j    > o > o >   At o:  h
    j-1  x ^ x ^ x
        i-1  i  i+1  At x & ^:
           i  i+1    At > & o:

The boundaries always run through q grid points (x).

Data Types
type	vertvisc_cs
	The control structure with parameters and memory for the MOM_vert_friction module. More...

Functions/Subroutines
subroutine, public	vertvisc (u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, taux_bot, tauy_bot, Waves)
	Perform a fully implicit vertical diffusion of momentum. Stress top and bottom boundary conditions are used. More...
subroutine, public	vertvisc_remnant (visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS)
	Calculate the fraction of momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of barotropic acceleration that a layer experiences after viscosity is applied. More...
subroutine, public	vertvisc_coef (u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
	Calculate the coupling coefficients (CSa_u and CSa_v) and effective layer thicknesses (CSh_u and CSh_v) for later use in the applying the implicit vertical viscosity via vertvisc(). More...
subroutine	find_coupling_coef (a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_ml, dt, j, G, GV, US, CS, visc, forces, work_on_u, OBC, shelf)
	Calculate the 'coupling coefficient' (a_cpl) at the interfaces. If BOTTOMDRAGLAW is defined, the minimum of Hbbl and half the adjacent layer thicknesses are used to calculate a_cpl near the bottom. More...
subroutine, public	vertvisc_limit_vel (u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS)
	Velocity components which exceed a threshold for physically reasonable values are truncated. Optionally, any column with excessive velocities may be sent to a diagnostic reporting subroutine. More...
subroutine, public	vertvisc_init (MIS, Time, G, GV, US, param_file, diag, ADp, dirs, ntrunc, CS)
	Initialize the vertical friction module. More...
subroutine, public	updatecfltruncationvalue (Time, CS, activate)
	Update the CFL truncation value as a function of time. If called with the optional argument activate=.true., record the value of Time as the beginning of the ramp period. More...
subroutine, public	vertvisc_end (CS)
	Clean up and deallocate the vertical friction module. More...

Function/Subroutine Documentation

find_coupling_coef()

subroutine mom_vert_friction::find_coupling_coef ( real, dimension(szib_(g),szk_(gv)+1), intent(out)  a_cpl, real, dimension(szib_(g),szk_(gv)), intent(in)  hvel, logical, dimension(szib_(g)), intent(in)  do_i, real, dimension(szib_(g),szk_(gv)), intent(in)  h_harm, real, dimension(szib_(g)), intent(in)  bbl_thick, real, dimension(szib_(g)), intent(in)  kv_bbl, real, dimension(szib_(g),szk_(gv)+1), intent(in)  z_i, real, dimension(szib_(g)), intent(out)  h_ml, real, intent(in)  dt, integer, intent(in)  j, type(ocean_grid_type), intent(in)  G, type(verticalgrid_type), intent(in)  GV, type(unit_scale_type), intent(in)  US, type(vertvisc_cs), pointer  CS, type(vertvisc_type), intent(in)  visc, type(mech_forcing), intent(in)  forces, logical, intent(in)  work_on_u, type(ocean_obc_type), pointer  OBC, logical, intent(in), optional  shelf  )

private

Calculate the 'coupling coefficient' (a_cpl) at the interfaces. If BOTTOMDRAGLAW is defined, the minimum of Hbbl and half the adjacent layer thicknesses are used to calculate a_cpl near the bottom.

Parameters

[in]	g	Ocean grid structure
[in]	gv	Ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[out]	a_cpl	Coupling coefficient across interfaces [Z T-1 ~> m s-1].
[in]	hvel	Thickness at velocity points [H ~> m or kg m-2]
[in]	do_i	If true, determine coupling coefficient for a column
[in]	h_harm	Harmonic mean of thicknesses around a velocity
[in]	bbl_thick	Bottom boundary layer thickness [H ~> m or kg m-2]
[in]	kv_bbl	Bottom boundary layer viscosity [Z2 T-1 ~> m2 s-1].
[in]	z_i	Estimate of interface heights above the bottom,
[out]	h_ml	Mixed layer depth [H ~> m or kg m-2]
[in]	j	j-index to find coupling coefficient for
[in]	dt	Time increment [T ~> s]
	cs	Vertical viscosity control structure
[in]	visc	Structure containing viscosities and bottom drag
[in]	forces	A structure with the driving mechanical forces
[in]	work_on_u	If true, u-points are being calculated, otherwise they are v-points
	obc	Open boundary condition structure
[in]	shelf	If present and true, use a surface boundary condition appropriate for an ice shelf.

Definition at line 1039 of file MOM_vert_friction.F90.

  type(ocean_grid_type),     intent(in)  :: G  !< Ocean grid structure
  type(verticalGrid_type),   intent(in)  :: GV !< Ocean vertical grid structure
  type(unit_scale_type),     intent(in)  :: US !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZK_(GV)+1), &
                             intent(out) :: a_cpl !< Coupling coefficient across interfaces [Z T-1 ~> m s-1].
  real, dimension(SZIB_(G),SZK_(GV)), &
                             intent(in)  :: hvel !< Thickness at velocity points [H ~> m or kg m-2]
  logical, dimension(SZIB_(G)), &
                             intent(in)  :: do_i !< If true, determine coupling coefficient for a column
  real, dimension(SZIB_(G),SZK_(GV)), &
                             intent(in)  :: h_harm !< Harmonic mean of thicknesses around a velocity
                                                   !! grid point [H ~> m or kg m-2]
  real, dimension(SZIB_(G)), intent(in)  :: bbl_thick !< Bottom boundary layer thickness [H ~> m or kg m-2]
  real, dimension(SZIB_(G)), intent(in)  :: kv_bbl !< Bottom boundary layer viscosity [Z2 T-1 ~> m2 s-1].
  real, dimension(SZIB_(G),SZK_(GV)+1), &
                             intent(in)  :: z_i  !< Estimate of interface heights above the bottom,
                                                 !! normalized by the bottom boundary layer thickness
  real, dimension(SZIB_(G)), intent(out) :: h_ml !< Mixed layer depth [H ~> m or kg m-2]
  integer,                   intent(in)  :: j    !< j-index to find coupling coefficient for
  real,                      intent(in)  :: dt   !< Time increment [T ~> s]
  type(vertvisc_CS),         pointer     :: CS   !< Vertical viscosity control structure
  type(vertvisc_type),       intent(in)  :: visc !< Structure containing viscosities and bottom drag
  type(mech_forcing),        intent(in)  :: forces !< A structure with the driving mechanical forces
  logical,                   intent(in)  :: work_on_u !< If true, u-points are being calculated,
                                                  !! otherwise they are v-points
  type(ocean_OBC_type),      pointer     :: OBC   !< Open boundary condition structure
  logical,         optional, intent(in)  :: shelf !< If present and true, use a surface boundary
                                                  !! condition appropriate for an ice shelf.
 
  ! Local variables
 
  real, dimension(SZIB_(G)) :: &
    u_star, &   ! ustar at a velocity point [Z T-1 ~> m s-1].
    absf, &     ! The average of the neighboring absolute values of f [T-1 ~> s-1].
!      h_ml, &  ! The mixed layer depth [H ~> m or kg m-2].
    nk_visc, &  ! The (real) interface index of the base of mixed layer.
    z_t, &      ! The distance from the top, sometimes normalized
                ! by Hmix, [H ~> m or kg m-2] or [nondim].
    kv_tbl, &   ! The viscosity in a top boundary layer under ice [Z2 T-1 ~> m2 s-1].
    tbl_thick
  real, dimension(SZIB_(G),SZK_(GV)+1) :: &
    Kv_tot, &   ! The total viscosity at an interface [Z2 T-1 ~> m2 s-1].
    Kv_add      ! A viscosity to add [Z2 T-1 ~> m2 s-1].
  real :: h_shear ! The distance over which shears occur [H ~> m or kg m-2].
  real :: r       ! A thickness to compare with Hbbl [H ~> m or kg m-2].
  real :: visc_ml ! The mixed layer viscosity [Z2 T-1 ~> m2 s-1].
  real :: I_Hmix  ! The inverse of the mixed layer thickness [H-1 ~> m-1 or m2 kg-1].
  real :: a_ml    ! The layer coupling coefficient across an interface in
                  ! the mixed layer [Z T-1 ~> m s-1].
  real :: I_amax  ! The inverse of the maximum coupling coefficient [T Z-1 ~> s m-1].
  real :: temp1   ! A temporary variable [H Z ~> m2 or kg m-1]
  real :: h_neglect   ! A thickness that is so small it is usually lost
                      ! in roundoff and can be neglected [H ~> m or kg m-2].
  real :: z2      ! A copy of z_i [nondim]
  real :: topfn   ! A function that is 1 at the top and small far from it [nondim]
  real :: kv_top  ! A viscosity associated with the top boundary layer [Z2 T-1 ~> m2 s-1]
  logical :: do_shelf, do_OBCs
  integer :: i, k, is, ie, max_nk
  integer :: nz
  real    :: botfn
 
  a_cpl(:,:) = 0.0
  kv_tot(:,:) = 0.0
 
  if (work_on_u) then ; is = g%IscB ; ie = g%IecB
  else ; is = g%isc ; ie = g%iec ; endif
  nz = g%ke
  h_neglect = gv%H_subroundoff
 
  !   The maximum coupling coefficent was originally introduced to avoid
  ! truncation error problems in the tridiagonal solver. Effectively, the 1e-10
  ! sets the maximum coupling coefficient increment to 1e10 m per timestep.
  if (cs%answers_2018) then
    i_amax = (1.0e-10*us%Z_to_m) * dt
  else
    i_amax = 0.0
  endif
 
  do_shelf = .false. ; if (present(shelf)) do_shelf = shelf
  do_obcs = .false.
  if (associated(obc)) then ; do_obcs = (obc%number_of_segments > 0) ; endif
  h_ml(:) = 0.0
 
!    The following loop calculates the vertical average velocity and
!  surface mixed layer contributions to the vertical viscosity.
  do i=is,ie ; kv_tot(i,1) = 0.0 ; enddo
  if ((gv%nkml>0) .or. do_shelf) then ; do k=2,nz ; do i=is,ie
    if (do_i(i)) kv_tot(i,k) = cs%Kv
  enddo ; enddo ; else
    i_hmix = 1.0 / (cs%Hmix + h_neglect)
    do i=is,ie ; z_t(i) = h_neglect*i_hmix ; enddo
    do k=2,nz ; do i=is,ie ; if (do_i(i)) then
      z_t(i) = z_t(i) + h_harm(i,k-1)*i_hmix
      kv_tot(i,k) = cs%Kv + cs%Kvml / ((z_t(i)*z_t(i)) *  &
               (1.0 + 0.09*z_t(i)*z_t(i)*z_t(i)*z_t(i)*z_t(i)*z_t(i)))
    endif ; enddo ; enddo
  endif
 
  do i=is,ie ; if (do_i(i)) then
    if (cs%bottomdraglaw) then
      r = hvel(i,nz)*0.5
      if (r < bbl_thick(i)) then
        a_cpl(i,nz+1) = kv_bbl(i) / (i_amax*kv_bbl(i) + (r+h_neglect)*gv%H_to_Z)
      else
        a_cpl(i,nz+1) = kv_bbl(i) / (i_amax*kv_bbl(i) + (bbl_thick(i)+h_neglect)*gv%H_to_Z)
      endif
    else
      a_cpl(i,nz+1) = cs%Kvbbl / ((0.5*hvel(i,nz)+h_neglect)*gv%H_to_Z + i_amax*cs%Kvbbl)
    endif
  endif ; enddo
 
  if (associated(visc%Kv_shear)) then
    ! The factor of 2 that used to be required in the viscosities is no longer needed.
    if (work_on_u) then
      do k=2,nz ; do i=is,ie ; if (do_i(i)) then
        kv_add(i,k) = 0.5*(visc%Kv_shear(i,j,k) + visc%Kv_shear(i+1,j,k))
      endif ; enddo ; enddo
      if (do_obcs) then
        do i=is,ie ; if (do_i(i) .and. (obc%segnum_u(i,j) /= obc_none)) then
          if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
            do k=2,nz ; kv_add(i,k) = visc%Kv_shear(i,j,k) ; enddo
          elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
            do k=2,nz ; kv_add(i,k) = visc%Kv_shear(i+1,j,k) ; enddo
          endif
        endif ; enddo
      endif
      do k=2,nz ; do i=is,ie ; if (do_i(i)) then
        kv_tot(i,k) = kv_tot(i,k) + kv_add(i,k)
      endif ; enddo ; enddo
    else
      do k=2,nz ; do i=is,ie ; if (do_i(i)) then
        kv_add(i,k) = 0.5*(visc%Kv_shear(i,j,k) + visc%Kv_shear(i,j+1,k))
      endif ; enddo ; enddo
      if (do_obcs) then
        do i=is,ie ; if (do_i(i) .and. (obc%segnum_v(i,j) /= obc_none)) then
          if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
            do k=2,nz ; kv_add(i,k) = visc%Kv_shear(i,j,k) ; enddo
          elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
            do k=2,nz ; kv_add(i,k) = visc%Kv_shear(i,j+1,k) ; enddo
          endif
        endif ; enddo
      endif
      do k=2,nz ; do i=is,ie ; if (do_i(i)) then
        kv_tot(i,k) = kv_tot(i,k) + kv_add(i,k)
      endif ; enddo ; enddo
    endif
  endif
 
  if (associated(visc%Kv_shear_Bu)) then
    if (work_on_u) then
      do k=2,nz ; do i=is,ie ; If (do_i(i)) then
        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))
      endif ; enddo ; enddo
    else
      do k=2,nz ; do i=is,ie ; if (do_i(i)) then
        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))
      endif ; enddo ; enddo
    endif
  endif
 
  do k=nz,2,-1 ; do i=is,ie ; if (do_i(i)) then
    !    botfn determines when a point is within the influence of the bottom
    !  boundary layer, going from 1 at the bottom to 0 in the interior.
    z2 = z_i(i,k)
    botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
 
    if (cs%bottomdraglaw) then
      kv_tot(i,k) = kv_tot(i,k) + (kv_bbl(i) - cs%Kv)*botfn
      r = 0.5*(hvel(i,k) + hvel(i,k-1))
      if (r > bbl_thick(i)) then
        h_shear = ((1.0 - botfn) * r + botfn*bbl_thick(i)) + h_neglect
      else
        h_shear = r + h_neglect
      endif
    else
      kv_tot(i,k) = kv_tot(i,k) + (cs%Kvbbl-cs%Kv)*botfn
      h_shear = 0.5*(hvel(i,k) + hvel(i,k-1) + h_neglect)
    endif
 
    ! Calculate the coupling coefficients from the viscosities.
    a_cpl(i,k) = kv_tot(i,k) / (h_shear*gv%H_to_Z + i_amax*kv_tot(i,k))
  endif ; enddo ; enddo ! i & k loops
 
  if (do_shelf) then
    ! Set the coefficients to include the no-slip surface stress.
    do i=is,ie ; if (do_i(i)) then
      if (work_on_u) then
        kv_tbl(i) = visc%Kv_tbl_shelf_u(i,j)
        tbl_thick(i) = visc%tbl_thick_shelf_u(i,j) * gv%Z_to_H + h_neglect
      else
        kv_tbl(i) = visc%Kv_tbl_shelf_v(i,j)
        tbl_thick(i) = visc%tbl_thick_shelf_v(i,j) * gv%Z_to_H + h_neglect
      endif
      z_t(i) = 0.0
 
      ! If a_cpl(i,1) were not already 0, it would be added here.
      if (0.5*hvel(i,1) > tbl_thick(i)) then
        a_cpl(i,1) = kv_tbl(i) / (tbl_thick(i)*gv%H_to_Z + i_amax*kv_tbl(i))
      else
        a_cpl(i,1) = kv_tbl(i) / ((0.5*hvel(i,1)+h_neglect)*gv%H_to_Z + i_amax*kv_tbl(i))
      endif
    endif ; enddo
 
    do k=2,nz ; do i=is,ie ;  if (do_i(i)) then
      z_t(i) = z_t(i) + hvel(i,k-1) / tbl_thick(i)
      topfn = 1.0 / (1.0 + 0.09 * z_t(i)**6)
 
      r = 0.5*(hvel(i,k)+hvel(i,k-1))
      if (r > tbl_thick(i)) then
        h_shear = ((1.0 - topfn) * r + topfn*tbl_thick(i)) + h_neglect
      else
        h_shear = r + h_neglect
      endif
 
      kv_top = topfn * kv_tbl(i)
      a_cpl(i,k) = a_cpl(i,k) + kv_top / (h_shear*gv%H_to_Z + i_amax*kv_top)
    endif ; enddo ; enddo
  elseif (cs%dynamic_viscous_ML .or. (gv%nkml>0)) then
    max_nk = 0
    do i=is,ie ; if (do_i(i)) then
      if (gv%nkml>0) nk_visc(i) = real(gv%nkml+1)
      if (work_on_u) then
        u_star(i) = 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j))
        absf(i) = 0.5*(abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i,j)))
        if (cs%dynamic_viscous_ML) nk_visc(i) = visc%nkml_visc_u(i,j) + 1
      else
        u_star(i) = 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1))
        absf(i) = 0.5*(abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j)))
        if (cs%dynamic_viscous_ML) nk_visc(i) = visc%nkml_visc_v(i,j) + 1
      endif
      h_ml(i) = h_neglect ; z_t(i) = 0.0
      max_nk = max(max_nk,ceiling(nk_visc(i) - 1.0))
    endif ; enddo
 
    if (do_obcs) then ; if (work_on_u) then
      do i=is,ie ; if (do_i(i) .and. (obc%segnum_u(i,j) /= obc_none)) then
        if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) &
          u_star(i) = forces%ustar(i,j)
        if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) &
          u_star(i) = forces%ustar(i+1,j)
      endif ; enddo
    else
      do i=is,ie ; if (do_i(i) .and. (obc%segnum_v(i,j) /= obc_none)) then
        if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) &
          u_star(i) = forces%ustar(i,j)
        if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) &
          u_star(i) = forces%ustar(i,j+1)
      endif ; enddo
    endif ; endif
 
    do k=1,max_nk ; do i=is,ie ; if (do_i(i)) then
      if (k+1 <= nk_visc(i)) then ! This layer is all in the ML.
        h_ml(i) = h_ml(i) + hvel(i,k)
      elseif (k < nk_visc(i)) then ! Part of this layer is in the ML.
        h_ml(i) = h_ml(i) + (nk_visc(i) - k) * hvel(i,k)
      endif
    endif ; enddo ; enddo
 
    do k=2,max_nk ; do i=is,ie ; if (do_i(i)) then ; if (k < nk_visc(i)) then
      ! Set the viscosity at the interfaces.
      z_t(i) = z_t(i) + hvel(i,k-1)
      temp1 = (z_t(i)*h_ml(i) - z_t(i)*z_t(i))*gv%H_to_Z
      !   This viscosity is set to go to 0 at the mixed layer top and bottom (in a log-layer)
      ! and be further limited by rotation to give the natural Ekman length.
      visc_ml = u_star(i) * 0.41 * (temp1*u_star(i)) / (absf(i)*temp1 + (h_ml(i)+h_neglect)*u_star(i))
      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)
      ! Choose the largest estimate of a.
      if (a_ml > a_cpl(i,k)) a_cpl(i,k) = a_ml
    endif ; endif ; enddo ; enddo
  endif
 

References mom_open_boundary::obc_direction_n, mom_open_boundary::obc_direction_s, and mom_open_boundary::obc_direction_w.

Referenced by vertvisc_coef().

Here is the caller graph for this function:

updatecfltruncationvalue()

subroutine, public mom_vert_friction::updatecfltruncationvalue	(	type(time_type), intent(in), target	Time,
		type(vertvisc_cs), pointer	CS,
		logical, intent(in), optional	activate
	)

Update the CFL truncation value as a function of time. If called with the optional argument activate=.true., record the value of Time as the beginning of the ramp period.

Parameters

[in]	time	Current model time
	cs	Vertical viscosity control structure
[in]	activate	Specifiy whether to record the value of Time as the beginning of the ramp period

Definition at line 1765 of file MOM_vert_friction.F90.

  type(time_type), target, intent(in)    :: Time     !< Current model time
  type(vertvisc_CS),       pointer       :: CS       !< Vertical viscosity control structure
  logical, optional,       intent(in)    :: activate !< Specifiy whether to record the value of
                                                     !! Time as the beginning of the ramp period
 
  ! Local variables
  real :: deltaTime, wghtA
  character(len=12) :: msg
 
  if (cs%truncRampTime==0.) return ! This indicates to ramping is turned off
 
  ! We use the optional argument to indicate this Time should be recorded as the
  ! beginning of the ramp-up period.
  if (present(activate)) then
    if (activate) then
      cs%rampStartTime = time ! Record the current time
      cs%CFLrampingIsActivated = .true.
    endif
  endif
  if (.not.cs%CFLrampingIsActivated) return
  deltatime = max( 0., time_type_to_real( time - cs%rampStartTime ) )
  if (deltatime >= cs%truncRampTime) then
    cs%CFL_trunc = cs%CFL_truncE
    cs%truncRampTime = 0. ! This turns off ramping after this call
  else
    wghta = min( 1., deltatime / cs%truncRampTime ) ! Linear profile in time
    !wghtA = wghtA*wghtA ! Convert linear profile to parabolic profile in time
    !wghtA = wghtA*wghtA*(3. - 2.*wghtA) ! Convert linear profile to cosine profile
    wghta = 1. - ( (1. - wghta)**2 ) ! Convert linear profiel to nverted parabolic profile
    cs%CFL_trunc = cs%CFL_truncS + wghta * ( cs%CFL_truncE - cs%CFL_truncS )
  endif
  write(msg(1:12),'(es12.3)') cs%CFL_trunc
  call mom_error(note, "MOM_vert_friction: updateCFLtruncationValue set CFL"// &
                       " limit to "//trim(msg))

References mom_error_handler::mom_error().

Referenced by mom_dynamics_split_rk2::initialize_dyn_split_rk2(), and mom_dynamics_split_rk2::step_mom_dyn_split_rk2().

Here is the call graph for this function:

Here is the caller graph for this function:

vertvisc()

subroutine, public mom_vert_friction::vertvisc	(	real, dimension( g %isdb: g %iedb, g %jsd: g %jed, gv %ke), intent(inout)	u,
		real, dimension( g %isd: g %ied, g %jsdb: g %jedb, gv %ke), intent(inout)	v,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)	h,
		type(mech_forcing), intent(in)	forces,
		type(vertvisc_type), intent(inout)	visc,
		real, intent(in)	dt,
		type(ocean_obc_type), pointer	OBC,
		type(accel_diag_ptrs), intent(inout)	ADp,
		type(cont_diag_ptrs), intent(inout)	CDp,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(vertvisc_cs), pointer	CS,
		real, dimension( g %isdb: g %iedb, g %jsd: g %jed), intent(out), optional	taux_bot,
		real, dimension( g %isd: g %ied, g %jsdb: g %jedb), intent(out), optional	tauy_bot,
		type(wave_parameters_cs), optional, pointer	Waves
	)

Perform a fully implicit vertical diffusion of momentum. Stress top and bottom boundary conditions are used.

This is solving the tridiagonal system

\[ \left(h_k + a_{k + 1/2} + a_{k - 1/2} + r_k\right) u_k^{n+1} = h_k u_k^n + a_{k + 1/2} u_{k+1}^{n+1} + a_{k - 1/2} u_{k-1}^{n+1} \]

where \(a_{k + 1/2} = \Delta t \nu_{k + 1/2} / h_{k + 1/2}\) is the interfacial coupling thickness per time step, encompassing background viscosity as well as contributions from enhanced mixed and bottom layer viscosities. $r_k$ is a Rayleight drag term due to channel drag. There is an additional stress term on the right-hand side if DIRECT_STRESS is true, applied to the surface layer.

Parameters

[in]	g	Ocean grid structure
[in]	gv	Ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in,out]	u	Zonal velocity [L T-1 ~> m s-1]
[in,out]	v	Meridional velocity [L T-1 ~> m s-1]
[in]	h	Layer thickness [H ~> m or kg m-2]
[in]	forces	A structure with the driving mechanical forces
[in,out]	visc	Viscosities and bottom drag
[in]	dt	Time increment [T ~> s]
	obc	Open boundary condition structure
[in,out]	adp	Accelerations in the momentum equations for diagnostics
[in,out]	cdp	Continuity equation terms
	cs	Vertical viscosity control structure
[out]	taux_bot	Zonal bottom stress from ocean to
[out]	tauy_bot	Meridional bottom stress from ocean to
	waves	Container for wave/Stokes information

Definition at line 151 of file MOM_vert_friction.F90.

  type(ocean_grid_type),   intent(in)    :: G      !< Ocean grid structure
  type(verticalGrid_type), intent(in)    :: GV     !< Ocean vertical grid structure
  type(unit_scale_type),   intent(in)    :: US     !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: u      !< Zonal velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
                           intent(inout) :: v      !< Meridional velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: h      !< Layer thickness [H ~> m or kg m-2]
  type(mech_forcing),    intent(in)      :: forces !< A structure with the driving mechanical forces
  type(vertvisc_type),   intent(inout)   :: visc   !< Viscosities and bottom drag
  real,                  intent(in)      :: dt     !< Time increment [T ~> s]
  type(ocean_OBC_type),  pointer         :: OBC    !< Open boundary condition structure
  type(accel_diag_ptrs), intent(inout)   :: ADp    !< Accelerations in the momentum
                                                   !! equations for diagnostics
  type(cont_diag_ptrs),  intent(inout)   :: CDp    !< Continuity equation terms
  type(vertvisc_CS),     pointer         :: CS     !< Vertical viscosity control structure
  real, dimension(SZIB_(G),SZJ_(G)), &
                   optional, intent(out) :: taux_bot !< Zonal bottom stress from ocean to
                                                     !! rock [R L Z T-2 ~> Pa]
  real, dimension(SZI_(G),SZJB_(G)), &
                   optional, intent(out) :: tauy_bot !< Meridional bottom stress from ocean to
                                                     !! rock [R L Z T-2 ~> Pa]
  type(wave_parameters_CS), &
                   optional, pointer     :: Waves !< Container for wave/Stokes information
 
  ! Fields from forces used in this subroutine:
  !   taux: Zonal wind stress [Pa].
  !   tauy: Meridional wind stress [Pa].
 
  ! Local variables
 
  real :: b1(SZIB_(G))          ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
  real :: c1(SZIB_(G),SZK_(G))  ! A variable used by the tridiagonal solver [nondim].
  real :: d1(SZIB_(G))          ! d1=1-c1 is used by the tridiagonal solver [nondim].
  real :: Ray(SZIB_(G),SZK_(G)) ! Ray is the Rayleigh-drag velocity [Z T-1 ~> m s-1].
  real :: b_denom_1             ! The first term in the denominator of b1 [H ~> m or kg m-2].
 
  real :: Hmix             ! The mixed layer thickness over which stress
                           ! is applied with direct_stress [H ~> m or kg m-2].
  real :: I_Hmix           ! The inverse of Hmix [H-1 ~> m-1 or m2 kg-1].
  real :: Idt              ! The inverse of the time step [T-1 ~> s-1].
  real :: dt_Rho0          ! The time step divided by the mean density [T H Z-1 R-1 ~> s m3 kg-1 or s].
  real :: dt_Z_to_H        ! The time step times the conversion from Z to the
                           ! units of thickness - [T H Z-1 ~> s or s kg m-3].
  real :: h_neglect        ! A thickness that is so small it is usually lost
                           ! in roundoff and can be neglected [H ~> m or kg m-2].
 
  real :: stress           !   The surface stress times the time step, divided
                           ! by the density [H L T-1 ~> m2 s-1 or kg m-1 s-1].
  real :: zDS, hfr, h_a    ! Temporary variables used with direct_stress.
  real :: surface_stress(SZIB_(G))! The same as stress, unless the wind stress
                           ! stress is applied as a body force [H L T-1 ~> m2 s-1 or kg m-1 s-1].
 
  logical :: do_i(SZIB_(G))
  logical :: DoStokesMixing
 
  integer :: i, j, k, is, ie, Isq, Ieq, Jsq, Jeq, nz, n
  is = g%isc ; ie = g%iec
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
 
  if (.not.associated(cs)) call mom_error(fatal,"MOM_vert_friction(visc): "// &
         "Module must be initialized before it is used.")
 
  if (cs%direct_stress) then
    hmix = cs%Hmix_stress
    i_hmix = 1.0 / hmix
  endif
  dt_rho0 = dt / gv%H_to_RZ
  dt_z_to_h = dt*gv%Z_to_H
  h_neglect = gv%H_subroundoff
  idt = 1.0 / dt
 
  !Check if Stokes mixing allowed if requested (present and associated)
  dostokesmixing=.false.
  if (cs%StokesMixing) then
    if (present(waves)) dostokesmixing = associated(waves)
    if (.not. dostokesmixing) &
      call mom_error(fatal,"Stokes Mixing called without allocated"//&
                     "Waves Control Structure")
  endif
 
  do k=1,nz ; do i=isq,ieq ; ray(i,k) = 0.0 ; enddo ; enddo
 
  !   Update the zonal velocity component using a modification of a standard
  ! tridagonal solver.
 
  !$OMP parallel do default(shared) firstprivate(Ray) &
  !$OMP                 private(do_i,surface_stress,zDS,stress,h_a,hfr, &
  !$OMP                         b_denom_1,b1,d1,c1)
  do j=g%jsc,g%jec
    do i=isq,ieq ; do_i(i) = (g%mask2dCu(i,j) > 0) ; enddo
 
    ! When mixing down Eulerian current + Stokes drift add before calling solver
    if (dostokesmixing) then ; do k=1,nz ; do i=isq,ieq
      if (do_i(i)) u(i,j,k) = u(i,j,k) + us%m_s_to_L_T*waves%Us_x(i,j,k)
    enddo ; enddo ; endif
 
    if (associated(adp%du_dt_visc)) then ; do k=1,nz ; do i=isq,ieq
      adp%du_dt_visc(i,j,k) = u(i,j,k)
    enddo ; enddo ; endif
 
!   One option is to have the wind stress applied as a body force
! over the topmost Hmix fluid.  If DIRECT_STRESS is not defined,
! the wind stress is applied as a stress boundary condition.
    if (cs%direct_stress) then
      do i=isq,ieq ; if (do_i(i)) then
        surface_stress(i) = 0.0
        zds = 0.0
        stress = dt_rho0 * forces%taux(i,j)
        do k=1,nz
          h_a = 0.5 * (h(i,j,k) + h(i+1,j,k)) + h_neglect
          hfr = 1.0 ; if ((zds+h_a) > hmix) hfr = (hmix - zds) / h_a
          u(i,j,k) = u(i,j,k) + i_hmix * hfr * stress
          zds = zds + h_a ; if (zds >= hmix) exit
        enddo
      endif ; enddo ! end of i loop
    else ; do i=isq,ieq
      surface_stress(i) = dt_rho0 * (g%mask2dCu(i,j)*forces%taux(i,j))
    enddo ; endif ! direct_stress
 
    if (cs%Channel_drag) then ; do k=1,nz ; do i=isq,ieq
      ray(i,k) = visc%Ray_u(i,j,k)
    enddo ; enddo ; endif
 
    ! perform forward elimination on the tridiagonal system
    !
    ! denote the diagonal of the system as b_k, the subdiagonal as a_k
    ! and the superdiagonal as c_k. The right-hand side terms are d_k.
    !
    ! ignoring the rayleigh drag contribution,
    ! we have a_k = -dt_Z_to_H * a_u(k)
    !         b_k = h_u(k) + dt_Z_to_H * (a_u(k) + a_u(k+1))
    !         c_k = -dt_Z_to_H * a_u(k+1)
    !
    ! for forward elimination, we want to:
    ! calculate c'_k = - c_k                / (b_k + a_k c'_(k-1))
    ! and       d'_k = (d_k - a_k d'_(k-1)) / (b_k + a_k c'_(k-1))
    ! where c'_1 = c_1/b_1 and d'_1 = d_1/b_1
    ! (see Thomas' tridiagonal matrix algorithm)
    !
    ! b1 is the denominator term 1 / (b_k + a_k c'_(k-1))
    ! b_denom_1 is (b_k + a_k + c_k) - a_k(1 - c'_(k-1))
    !            = (b_k + c_k + c'_(k-1))
    ! this is done so that d1 = b1 * b_denom_1 = 1 - c'_(k-1)
    ! c1(k) is -c'_(k - 1)
    ! and the right-hand-side is destructively updated to be d'_k
    !
    do i=isq,ieq ; if (do_i(i)) then
      b_denom_1 = cs%h_u(i,j,1) + dt_z_to_h * (ray(i,1) + cs%a_u(i,j,1))
      b1(i) = 1.0 / (b_denom_1 + dt_z_to_h*cs%a_u(i,j,2))
      d1(i) = b_denom_1 * b1(i)
      u(i,j,1) = b1(i) * (cs%h_u(i,j,1) * u(i,j,1) + surface_stress(i))
    endif ; enddo
    do k=2,nz ; do i=isq,ieq ; if (do_i(i)) then
      c1(i,k) = dt_z_to_h * cs%a_u(i,j,k) * b1(i)
      b_denom_1 = cs%h_u(i,j,k) + dt_z_to_h * (ray(i,k) + cs%a_u(i,j,k)*d1(i))
      b1(i) = 1.0 / (b_denom_1 + dt_z_to_h * cs%a_u(i,j,k+1))
      d1(i) = b_denom_1 * b1(i)
      u(i,j,k) = (cs%h_u(i,j,k) * u(i,j,k) + &
                  dt_z_to_h * cs%a_u(i,j,k) * u(i,j,k-1)) * b1(i)
    endif ; enddo ; enddo
 
    ! back substitute to solve for the new velocities
    ! u_k = d'_k - c'_k x_(k+1)
    do k=nz-1,1,-1 ; do i=isq,ieq ; if (do_i(i)) then
      u(i,j,k) = u(i,j,k) + c1(i,k+1) * u(i,j,k+1)
    endif ; enddo ; enddo ! i and k loops
 
    if (associated(adp%du_dt_visc)) then ; do k=1,nz ; do i=isq,ieq
      adp%du_dt_visc(i,j,k) = (u(i,j,k) - adp%du_dt_visc(i,j,k))*idt
    enddo ; enddo ; endif
 
    if (associated(visc%taux_shelf)) then ; do i=isq,ieq
      visc%taux_shelf(i,j) = -gv%Rho0*cs%a1_shelf_u(i,j)*u(i,j,1) ! - u_shelf?
    enddo ; endif
 
    if (PRESENT(taux_bot)) then
      do i=isq,ieq
        taux_bot(i,j) = gv%Rho0 * (u(i,j,nz)*cs%a_u(i,j,nz+1))
      enddo
      if (cs%Channel_drag) then ; do k=1,nz ; do i=isq,ieq
        taux_bot(i,j) = taux_bot(i,j) + gv%Rho0 * (ray(i,k)*u(i,j,k))
      enddo ; enddo ; endif
    endif
 
    ! When mixing down Eulerian current + Stokes drift subtract after calling solver
    if (dostokesmixing) then ; do k=1,nz ; do i=isq,ieq
      if (do_i(i)) u(i,j,k) = u(i,j,k) - us%m_s_to_L_T*waves%Us_x(i,j,k)
    enddo ; enddo ; endif
 
  enddo ! end u-component j loop
 
  ! Now work on the meridional velocity component.
 
  !$OMP parallel do default(shared) firstprivate(Ray) &
  !$OMP               private(do_i,surface_stress,zDS,stress,h_a,hfr, &
  !$OMP                       b_denom_1,b1,d1,c1)
  do j=jsq,jeq
    do i=is,ie ; do_i(i) = (g%mask2dCv(i,j) > 0) ; enddo
 
    ! When mixing down Eulerian current + Stokes drift add before calling solver
    if (dostokesmixing) then ; do k=1,nz ; do i=is,ie
      if (do_i(i)) v(i,j,k) = v(i,j,k) + us%m_s_to_L_T*waves%Us_y(i,j,k)
    enddo ; enddo ; endif
 
    if (associated(adp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie
      adp%dv_dt_visc(i,j,k) = v(i,j,k)
    enddo ; enddo ; endif
 
!   One option is to have the wind stress applied as a body force
! over the topmost Hmix fluid.  If DIRECT_STRESS is not defined,
! the wind stress is applied as a stress boundary condition.
    if (cs%direct_stress) then
      do i=is,ie ; if (do_i(i)) then
        surface_stress(i) = 0.0
        zds = 0.0
        stress = dt_rho0 * forces%tauy(i,j)
        do k=1,nz
          h_a = 0.5 * (h(i,j,k) + h(i,j+1,k)) + h_neglect
          hfr = 1.0 ; if ((zds+h_a) > hmix) hfr = (hmix - zds) / h_a
          v(i,j,k) = v(i,j,k) + i_hmix * hfr * stress
          zds = zds + h_a ; if (zds >= hmix) exit
        enddo
      endif ; enddo ! end of i loop
    else ; do i=is,ie
      surface_stress(i) = dt_rho0 * (g%mask2dCv(i,j)*forces%tauy(i,j))
    enddo ; endif ! direct_stress
 
    if (cs%Channel_drag) then ; do k=1,nz ; do i=is,ie
      ray(i,k) = visc%Ray_v(i,j,k)
    enddo ; enddo ; endif
 
    do i=is,ie ; if (do_i(i)) then
      b_denom_1 = cs%h_v(i,j,1) + dt_z_to_h * (ray(i,1) + cs%a_v(i,j,1))
      b1(i) = 1.0 / (b_denom_1 + dt_z_to_h*cs%a_v(i,j,2))
      d1(i) = b_denom_1 * b1(i)
      v(i,j,1) = b1(i) * (cs%h_v(i,j,1) * v(i,j,1) + surface_stress(i))
    endif ; enddo
    do k=2,nz ; do i=is,ie ; if (do_i(i)) then
      c1(i,k) = dt_z_to_h * cs%a_v(i,j,k) * b1(i)
      b_denom_1 = cs%h_v(i,j,k) + dt_z_to_h * (ray(i,k) + cs%a_v(i,j,k)*d1(i))
      b1(i) = 1.0 / (b_denom_1 + dt_z_to_h * cs%a_v(i,j,k+1))
      d1(i) = b_denom_1 * b1(i)
      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)
    endif ; enddo ; enddo
    do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then
      v(i,j,k) = v(i,j,k) + c1(i,k+1) * v(i,j,k+1)
    endif ; enddo ; enddo ! i and k loops
 
    if (associated(adp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie
      adp%dv_dt_visc(i,j,k) = (v(i,j,k) - adp%dv_dt_visc(i,j,k))*idt
    enddo ; enddo ; endif
 
    if (associated(visc%tauy_shelf)) then ; do i=is,ie
      visc%tauy_shelf(i,j) = -gv%Rho0*cs%a1_shelf_v(i,j)*v(i,j,1) ! - v_shelf?
    enddo ; endif
 
    if (present(tauy_bot)) then
      do i=is,ie
        tauy_bot(i,j) = gv%Rho0 * (v(i,j,nz)*cs%a_v(i,j,nz+1))
      enddo
      if (cs%Channel_drag) then ; do k=1,nz ; do i=is,ie
        tauy_bot(i,j) = tauy_bot(i,j) + gv%Rho0 * (ray(i,k)*v(i,j,k))
      enddo ; enddo ; endif
    endif
 
    ! When mixing down Eulerian current + Stokes drift subtract after calling solver
    if (dostokesmixing) then ; do k=1,nz ; do i=is,ie
      if (do_i(i)) v(i,j,k) = v(i,j,k) - us%m_s_to_L_T*waves%Us_y(i,j,k)
    enddo ; enddo ; endif
 
  enddo ! end of v-component J loop
 
  call vertvisc_limit_vel(u, v, h, adp, cdp, forces, visc, dt, g, gv, us, cs)
 
  ! Here the velocities associated with open boundary conditions are applied.
  if (associated(obc)) then
    do n=1,obc%number_of_segments
      if (obc%segment(n)%specified) then
        if (obc%segment(n)%is_N_or_S) then
          j = obc%segment(n)%HI%JsdB
          do k=1,nz ; do i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied
            v(i,j,k) = obc%segment(n)%normal_vel(i,j,k)
          enddo ; enddo
        elseif (obc%segment(n)%is_E_or_W) then
          i = obc%segment(n)%HI%IsdB
          do k=1,nz ; do j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed
            u(i,j,k) = obc%segment(n)%normal_vel(i,j,k)
          enddo ; enddo
        endif
      endif
    enddo
  endif
 
! Offer diagnostic fields for averaging.
  if (cs%id_du_dt_visc > 0) &
    call post_data(cs%id_du_dt_visc, adp%du_dt_visc, cs%diag)
  if (cs%id_dv_dt_visc > 0) &
    call post_data(cs%id_dv_dt_visc, adp%dv_dt_visc, cs%diag)
  if (present(taux_bot) .and. (cs%id_taux_bot > 0)) &
    call post_data(cs%id_taux_bot, taux_bot, cs%diag)
  if (present(tauy_bot) .and. (cs%id_tauy_bot > 0)) &
    call post_data(cs%id_tauy_bot, tauy_bot, cs%diag)
 

References mom_error_handler::mom_error(), and vertvisc_limit_vel().

Referenced by mom_dynamics_unsplit::step_mom_dyn_unsplit(), and mom_dynamics_unsplit_rk2::step_mom_dyn_unsplit_rk2().

Here is the call graph for this function:

Here is the caller graph for this function:

vertvisc_coef()

subroutine, public mom_vert_friction::vertvisc_coef	(	real, dimension(szib_(g),szj_(g),szk_(gv)), intent(in)	u,
		real, dimension(szi_(g),szjb_(g),szk_(gv)), intent(in)	v,
		real, dimension(szi_(g),szj_(g),szk_(gv)), intent(in)	h,
		type(mech_forcing), intent(in)	forces,
		type(vertvisc_type), intent(in)	visc,
		real, intent(in)	dt,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(vertvisc_cs), pointer	CS,
		type(ocean_obc_type), pointer	OBC
	)

Calculate the coupling coefficients (CSa_u and CSa_v) and effective layer thicknesses (CSh_u and CSh_v) for later use in the applying the implicit vertical viscosity via vertvisc().

Parameters

[in]	g	Ocean grid structure
[in]	gv	Ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	u	Zonal velocity [L T-1 ~> m s-1]
[in]	v	Meridional velocity [L T-1 ~> m s-1]
[in]	h	Layer thickness [H ~> m or kg m-2]
[in]	forces	A structure with the driving mechanical forces
[in]	visc	Viscosities and bottom drag
[in]	dt	Time increment [T ~> s]
	cs	Vertical viscosity control structure
	obc	Open boundary condition structure

Definition at line 571 of file MOM_vert_friction.F90.

  type(ocean_grid_type),   intent(in)    :: G      !< Ocean grid structure
  type(verticalGrid_type), intent(in)    :: GV     !< Ocean vertical grid structure
  type(unit_scale_type),   intent(in)    :: US     !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: u      !< Zonal velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
                           intent(in)    :: v      !< Meridional velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: h      !< Layer thickness [H ~> m or kg m-2]
  type(mech_forcing),      intent(in)    :: forces !< A structure with the driving mechanical forces
  type(vertvisc_type),     intent(in)    :: visc   !< Viscosities and bottom drag
  real,                    intent(in)    :: dt     !< Time increment [T ~> s]
  type(vertvisc_CS),       pointer       :: CS     !< Vertical viscosity control structure
  type(ocean_OBC_type),    pointer       :: OBC    !< Open boundary condition structure
 
  ! Field from forces used in this subroutine:
  !   ustar: the friction velocity [Z T-1 ~> m s-1], used here as the mixing
  !     velocity in the mixed layer if NKML > 1 in a bulk mixed layer.
 
  ! Local variables
 
  real, dimension(SZIB_(G),SZK_(G)) :: &
    h_harm, &   ! Harmonic mean of the thicknesses around a velocity grid point,
                ! given by 2*(h+ * h-)/(h+ + h-) [H ~> m or kg m-2].
    h_arith, &  ! The arithmetic mean thickness [H ~> m or kg m-2].
    h_delta, &  ! The lateral difference of thickness [H ~> m or kg m-2].
    hvel, &     ! hvel is the thickness used at a velocity grid point [H ~> m or kg m-2].
    hvel_shelf  ! The equivalent of hvel under shelves [H ~> m or kg m-2].
  real, dimension(SZIB_(G),SZK_(G)+1) :: &
    a_cpl, &    ! The drag coefficients across interfaces [Z T-1 ~> m s-1].  a_cpl times
                ! the velocity difference gives the stress across an interface.
    a_shelf, &  ! The drag coefficients across interfaces in water columns under
                ! ice shelves [Z T-1 ~> m s-1].
    z_i         ! An estimate of each interface's height above the bottom,
                ! normalized by the bottom boundary layer thickness, nondim.
  real, dimension(SZIB_(G)) :: &
    kv_bbl, &     ! The bottom boundary layer viscosity [Z2 T-1 ~> m2 s-1].
    bbl_thick, &  ! The bottom boundary layer thickness [H ~> m or kg m-2].
    I_Hbbl, &     ! The inverse of the bottom boundary layer thickness [H-1 ~> m-1 or m2 kg-1].
    I_Htbl, &     ! The inverse of the top boundary layer thickness [H-1 ~> m-1 or m2 kg-1].
    zcol1, &      ! The height of the interfaces to the north and south of a
    zcol2, &      ! v-point [H ~> m or kg m-2].
    Ztop_min, &   ! The deeper of the two adjacent surface heights [H ~> m or kg m-2].
    Dmin, &       ! The shallower of the two adjacent bottom depths converted to
                  ! thickness units [H ~> m or kg m-2].
    zh, &         ! An estimate of the interface's distance from the bottom
                  ! based on harmonic mean thicknesses [H ~> m or kg m-2].
    h_ml          ! The mixed layer depth [H ~> m or kg m-2].
  real, allocatable, dimension(:,:) :: hML_u ! Diagnostic of the mixed layer depth at u points [H ~> m or kg m-2].
  real, allocatable, dimension(:,:) :: hML_v ! Diagnostic of the mixed layer depth at v points [H ~> m or kg m-2].
  real, allocatable, dimension(:,:,:) :: Kv_u !< Total vertical viscosity at u-points [Z2 T-1 ~> m2 s-1].
  real, allocatable, dimension(:,:,:) :: Kv_v !< Total vertical viscosity at v-points [Z2 T-1 ~> m2 s-1].
  real :: zcol(SZI_(G)) ! The height of an interface at h-points [H ~> m or kg m-2].
  real :: botfn   ! A function which goes from 1 at the bottom to 0 much more
                  ! than Hbbl into the interior.
  real :: topfn   ! A function which goes from 1 at the top to 0 much more
                  ! than Htbl into the interior.
  real :: z2      ! The distance from the bottom, normalized by Hbbl, nondim.
  real :: z2_wt   ! A nondimensional (0-1) weight used when calculating z2.
  real :: z_clear ! The clearance of an interface above the surrounding topography [H ~> m or kg m-2].
  real :: h_neglect  ! A thickness that is so small it is usually lost
                     ! in roundoff and can be neglected [H ~> m or kg m-2].
 
  real :: I_valBL ! The inverse of a scaling factor determining when water is
                  ! still within the boundary layer, as determined by the sum
                  ! of the harmonic mean thicknesses.
  logical, dimension(SZIB_(G)) :: do_i, do_i_shelf
  logical :: do_any_shelf
  integer, dimension(SZIB_(G)) :: &
    zi_dir   !  A trinary logical array indicating which thicknesses to use for
             !  finding z_clear.
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
 
  if (.not.associated(cs)) call mom_error(fatal,"MOM_vert_friction(coef): "// &
         "Module must be initialized before it is used.")
 
  h_neglect = gv%H_subroundoff
  i_hbbl(:) = 1.0 / (cs%Hbbl + h_neglect)
  i_valbl = 0.0 ; if (cs%harm_BL_val > 0.0) i_valbl = 1.0 / cs%harm_BL_val
 
  if (cs%id_Kv_u > 0) then
    allocate(kv_u(g%IsdB:g%IedB,g%jsd:g%jed,g%ke)) ; kv_u(:,:,:) = 0.0
  endif
 
  if (cs%id_Kv_v > 0) then
    allocate(kv_v(g%isd:g%ied,g%JsdB:g%JedB,g%ke)) ; kv_v(:,:,:) = 0.0
  endif
 
  if (cs%debug .or. (cs%id_hML_u > 0)) then
    allocate(hml_u(g%IsdB:g%IedB,g%jsd:g%jed)) ; hml_u(:,:) = 0.0
  endif
  if (cs%debug .or. (cs%id_hML_v > 0)) then
    allocate(hml_v(g%isd:g%ied,g%JsdB:g%JedB)) ; hml_v(:,:) = 0.0
  endif
 
  if ((associated(visc%taux_shelf) .or. associated(forces%frac_shelf_u)) .and. &
      .not.associated(cs%a1_shelf_u)) then
    allocate(cs%a1_shelf_u(g%IsdB:g%IedB,g%jsd:g%jed)) ; cs%a1_shelf_u(:,:)=0.0
  endif
  if ((associated(visc%tauy_shelf) .or. associated(forces%frac_shelf_v)) .and. &
      .not.associated(cs%a1_shelf_v)) then
    allocate(cs%a1_shelf_v(g%isd:g%ied,g%JsdB:g%JedB)) ; cs%a1_shelf_v(:,:)=0.0
  endif
 
  !$OMP parallel do default(private) shared(G,GV,US,CS,visc,Isq,Ieq,nz,u,h,forces,hML_u, &
  !$OMP                                     OBC,h_neglect,dt,I_valBL,Kv_u) &
  !$OMP                     firstprivate(i_hbbl)
  do j=g%Jsc,g%Jec
    do i=isq,ieq ; do_i(i) = (g%mask2dCu(i,j) > 0) ; enddo
 
    if (cs%bottomdraglaw) then ; do i=isq,ieq
      kv_bbl(i) = visc%Kv_bbl_u(i,j)
      bbl_thick(i) = visc%bbl_thick_u(i,j) * gv%Z_to_H + h_neglect
      if (do_i(i)) i_hbbl(i) = 1.0 / bbl_thick(i)
    enddo ; endif
 
    do k=1,nz ; do i=isq,ieq ; if (do_i(i)) then
      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)
      h_arith(i,k) = 0.5*(h(i+1,j,k)+h(i,j,k))
      h_delta(i,k) = h(i+1,j,k) - h(i,j,k)
    endif ; enddo ; enddo
    do i=isq,ieq
      dmin(i) = min(g%bathyT(i,j), g%bathyT(i+1,j)) * gv%Z_to_H
      zi_dir(i) = 0
    enddo
 
    ! Project thickness outward across OBCs using a zero-gradient condition.
    if (associated(obc)) then ; if (obc%number_of_segments > 0) then
      do i=isq,ieq ; if (do_i(i) .and. (obc%segnum_u(i,j) /= obc_none)) then
        if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
          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
          dmin(i) = g%bathyT(i,j) * gv%Z_to_H
          zi_dir(i) = -1
        elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
          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
          dmin(i) = g%bathyT(i+1,j) * gv%Z_to_H
          zi_dir(i) = 1
        endif
      endif ; enddo
    endif ; endif
 
!    The following block calculates the thicknesses at velocity
!  grid points for the vertical viscosity (hvel).  Near the
!  bottom an upwind biased thickness is used to control the effect
!  of spurious Montgomery potential gradients at the bottom where
!  nearly massless layers layers ride over the topography.
    if (cs%harmonic_visc) then
      do i=isq,ieq ; z_i(i,nz+1) = 0.0 ; enddo
      do k=nz,1,-1 ; do i=isq,ieq ; if (do_i(i)) then
        hvel(i,k) = h_harm(i,k)
        if (u(i,j,k) * h_delta(i,k) < 0) then
          z2 = z_i(i,k+1) ; botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
          hvel(i,k) = (1.0-botfn)*h_harm(i,k) + botfn*h_arith(i,k)
        endif
        z_i(i,k) =  z_i(i,k+1) + h_harm(i,k)*i_hbbl(i)
      endif ; enddo ; enddo ! i & k loops
    else ! Not harmonic_visc
      do i=isq,ieq ; zh(i) = 0.0 ; z_i(i,nz+1) = 0.0 ; enddo
      do i=isq,ieq+1 ; zcol(i) = -g%bathyT(i,j) * gv%Z_to_H ; enddo
      do k=nz,1,-1
        do i=isq,ieq+1 ; zcol(i) = zcol(i) + h(i,j,k) ; enddo
        do i=isq,ieq ; if (do_i(i)) then
          zh(i) = zh(i) + h_harm(i,k)
 
          z_clear = max(zcol(i),zcol(i+1)) + dmin(i)
          if (zi_dir(i) < 0) z_clear = zcol(i) + dmin(i)
          if (zi_dir(i) > 0) z_clear = zcol(i+1) + dmin(i)
 
          z_i(i,k) = max(zh(i), z_clear) * i_hbbl(i)
 
          hvel(i,k) = h_arith(i,k)
          if (u(i,j,k) * h_delta(i,k) > 0) then
            if (zh(i) * i_hbbl(i) < cs%harm_BL_val) then
              hvel(i,k) = h_harm(i,k)
            else
              z2_wt = 1.0  ; if (zh(i) * i_hbbl(i) < 2.0*cs%harm_BL_val) &
                z2_wt = max(0.0, min(1.0, zh(i) * i_hbbl(i) * i_valbl - 1.0))
              z2 = z2_wt * (max(zh(i), z_clear) * i_hbbl(i))
              botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
              hvel(i,k) = (1.0-botfn)*h_arith(i,k) + botfn*h_harm(i,k)
            endif
          endif
 
        endif ; enddo ! i  loop
      enddo ! k loop
    endif
 
    call find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_ml, &
                            dt, j, g, gv, us, cs, visc, forces, work_on_u=.true., obc=obc)
    if (allocated(hml_u)) then
      do i=isq,ieq ; if (do_i(i)) then ; hml_u(i,j) = h_ml(i) ; endif ; enddo
    endif
 
    do_any_shelf = .false.
    if (associated(forces%frac_shelf_u)) then
      do i=isq,ieq
        cs%a1_shelf_u(i,j) = 0.0
        do_i_shelf(i) = (do_i(i) .and. forces%frac_shelf_u(i,j) > 0.0)
        if (do_i_shelf(i)) do_any_shelf = .true.
      enddo
      if (do_any_shelf) then
        if (cs%harmonic_visc) then
          call find_coupling_coef(a_shelf, hvel, do_i_shelf, h_harm, bbl_thick, &
                                  kv_bbl, z_i, h_ml, dt, j, g, gv, us, cs, &
                                  visc, forces, work_on_u=.true., obc=obc, shelf=.true.)
        else  ! Find upwind-biased thickness near the surface.
          ! Perhaps this needs to be done more carefully, via find_eta.
          do i=isq,ieq ; if (do_i_shelf(i)) then
            zh(i) = 0.0 ; ztop_min(i) = min(zcol(i), zcol(i+1))
            i_htbl(i) = 1.0 / (visc%tbl_thick_shelf_u(i,j)*gv%Z_to_H + h_neglect)
          endif ; enddo
          do k=1,nz
            do i=isq,ieq+1 ; zcol(i) = zcol(i) - h(i,j,k) ; enddo
            do i=isq,ieq ; if (do_i_shelf(i)) then
              zh(i) = zh(i) + h_harm(i,k)
 
              hvel_shelf(i,k) = hvel(i,k)
              if (u(i,j,k) * h_delta(i,k) > 0) then
                if (zh(i) * i_htbl(i) < cs%harm_BL_val) then
                  hvel_shelf(i,k) = min(hvel(i,k), h_harm(i,k))
                else
                  z2_wt = 1.0  ; if (zh(i) * i_htbl(i) < 2.0*cs%harm_BL_val) &
                    z2_wt = max(0.0, min(1.0, zh(i) * i_htbl(i) * i_valbl - 1.0))
                  z2 = z2_wt * (max(zh(i), ztop_min(i) - min(zcol(i),zcol(i+1))) * i_htbl(i))
                  topfn = 1.0 / (1.0 + 0.09*z2**6)
                  hvel_shelf(i,k) = min(hvel(i,k), (1.0-topfn)*h_arith(i,k) + topfn*h_harm(i,k))
                endif
              endif
            endif ; enddo
          enddo
          call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, h_harm, &
                                  bbl_thick, kv_bbl, z_i, h_ml, dt, j, g, gv, us, cs, &
                                  visc, forces, work_on_u=.true., obc=obc, shelf=.true.)
        endif
        do i=isq,ieq ; if (do_i_shelf(i)) cs%a1_shelf_u(i,j) = a_shelf(i,1) ; enddo
      endif
    endif
 
    if (do_any_shelf) then
      do k=1,nz+1 ; do i=isq,ieq ; if (do_i_shelf(i)) then
        cs%a_u(i,j,k) = forces%frac_shelf_u(i,j)  * a_shelf(i,k) + &
                   (1.0-forces%frac_shelf_u(i,j)) * a_cpl(i,k)
! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH
!        CS%a_u(I,j,K) = forces%frac_shelf_u(I,j)  * max(a_shelf(I,K), a_cpl(I,K)) + &
!                   (1.0-forces%frac_shelf_u(I,j)) * a_cpl(I,K)
      elseif (do_i(i)) then
        cs%a_u(i,j,k) = a_cpl(i,k)
      endif ; enddo ; enddo
      do k=1,nz ; do i=isq,ieq ; if (do_i_shelf(i)) then
        ! Should we instead take the inverse of the average of the inverses?
        cs%h_u(i,j,k) = forces%frac_shelf_u(i,j)  * hvel_shelf(i,k) + &
                   (1.0-forces%frac_shelf_u(i,j)) * hvel(i,k)
      elseif (do_i(i)) then
        cs%h_u(i,j,k) = hvel(i,k)
      endif ; enddo ; enddo
    else
      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
      do k=1,nz ; do i=isq,ieq ; if (do_i(i)) cs%h_u(i,j,k) = hvel(i,k) ; enddo ; enddo
    endif
 
    ! Diagnose total Kv at u-points
    if (cs%id_Kv_u > 0) then
      do k=1,nz ; do i=isq,ieq
        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)
      enddo ; enddo
    endif
 
  enddo
 
 
  ! Now work on v-points.
  !$OMP parallel do default(private) shared(G,GV,CS,US,visc,is,ie,Jsq,Jeq,nz,v,h,forces,hML_v, &
  !$OMP                                  OBC,h_neglect,dt,I_valBL,Kv_v) &
  !$OMP                     firstprivate(i_hbbl)
  do j=jsq,jeq
    do i=is,ie ; do_i(i) = (g%mask2dCv(i,j) > 0) ; enddo
 
    if (cs%bottomdraglaw) then ; do i=is,ie
      kv_bbl(i) = visc%Kv_bbl_v(i,j)
      bbl_thick(i) = visc%bbl_thick_v(i,j) * gv%Z_to_H + h_neglect
      if (do_i(i)) i_hbbl(i) = 1.0 / bbl_thick(i)
    enddo ; endif
 
    do k=1,nz ; do i=is,ie ; if (do_i(i)) then
      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)
      h_arith(i,k) = 0.5*(h(i,j+1,k)+h(i,j,k))
      h_delta(i,k) = h(i,j+1,k) - h(i,j,k)
    endif ; enddo ; enddo
    do i=is,ie
      dmin(i) = min(g%bathyT(i,j), g%bathyT(i,j+1)) * gv%Z_to_H
      zi_dir(i) = 0
    enddo
 
    ! Project thickness outward across OBCs using a zero-gradient condition.
    if (associated(obc)) then ; if (obc%number_of_segments > 0) then
      do i=is,ie ; if (do_i(i) .and. (obc%segnum_v(i,j) /= obc_none)) then
        if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
          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
          dmin(i) = g%bathyT(i,j) * gv%Z_to_H
          zi_dir(i) = -1
        elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
          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
          dmin(i) = g%bathyT(i,j+1) * gv%Z_to_H
          zi_dir(i) = 1
        endif
      endif ; enddo
    endif ; endif
 
!    The following block calculates the thicknesses at velocity
!  grid points for the vertical viscosity (hvel).  Near the
!  bottom an upwind biased thickness is used to control the effect
!  of spurious Montgomery potential gradients at the bottom where
!  nearly massless layers layers ride over the topography.
    if (cs%harmonic_visc) then
      do i=is,ie ; z_i(i,nz+1) = 0.0 ; enddo
 
      do k=nz,1,-1 ; do i=is,ie ; if (do_i(i)) then
        hvel(i,k) = h_harm(i,k)
        if (v(i,j,k) * h_delta(i,k) < 0) then
          z2 = z_i(i,k+1) ; botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
          hvel(i,k) = (1.0-botfn)*h_harm(i,k) + botfn*h_arith(i,k)
        endif
        z_i(i,k) = z_i(i,k+1)  + h_harm(i,k)*i_hbbl(i)
      endif ; enddo ; enddo ! i & k loops
    else ! Not harmonic_visc
      do i=is,ie
        zh(i) = 0.0 ; z_i(i,nz+1) = 0.0
        zcol1(i) = -g%bathyT(i,j) * gv%Z_to_H
        zcol2(i) = -g%bathyT(i,j+1) * gv%Z_to_H
      enddo
      do k=nz,1,-1 ; do i=is,ie ; if (do_i(i)) then
        zh(i) = zh(i) + h_harm(i,k)
        zcol1(i) = zcol1(i) + h(i,j,k) ; zcol2(i) = zcol2(i) + h(i,j+1,k)
 
        z_clear = max(zcol1(i),zcol2(i)) + dmin(i)
        if (zi_dir(i) < 0) z_clear = zcol1(i) + dmin(i)
        if (zi_dir(i) > 0) z_clear = zcol2(i) + dmin(i)
 
        z_i(i,k) = max(zh(i), z_clear) * i_hbbl(i)
 
        hvel(i,k) = h_arith(i,k)
        if (v(i,j,k) * h_delta(i,k) > 0) then
          if (zh(i) * i_hbbl(i) < cs%harm_BL_val) then
            hvel(i,k) = h_harm(i,k)
          else
            z2_wt = 1.0  ; if (zh(i) * i_hbbl(i) < 2.0*cs%harm_BL_val) &
              z2_wt = max(0.0, min(1.0, zh(i) * i_hbbl(i) * i_valbl - 1.0))
            z2 = z2_wt * (max(zh(i), max(zcol1(i),zcol2(i)) + dmin(i)) * i_hbbl(i))
            botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
            hvel(i,k) = (1.0-botfn)*h_arith(i,k) + botfn*h_harm(i,k)
          endif
       endif
 
      endif ; enddo ; enddo ! i & k loops
    endif
 
    call find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_ml, &
                            dt, j, g, gv, us, cs, visc, forces, work_on_u=.false., obc=obc)
    if ( allocated(hml_v)) then
       do i=is,ie ; if (do_i(i)) then ; hml_v(i,j) = h_ml(i) ; endif ; enddo
    endif
    do_any_shelf = .false.
    if (associated(forces%frac_shelf_v)) then
      do i=is,ie
        cs%a1_shelf_v(i,j) = 0.0
        do_i_shelf(i) = (do_i(i) .and. forces%frac_shelf_v(i,j) > 0.0)
        if (do_i_shelf(i)) do_any_shelf = .true.
      enddo
      if (do_any_shelf) then
        if (cs%harmonic_visc) then
          call find_coupling_coef(a_shelf, hvel, do_i_shelf, h_harm, bbl_thick, &
                                  kv_bbl, z_i, h_ml, dt, j, g, gv, us, cs, visc, &
                                  forces, work_on_u=.false., obc=obc, shelf=.true.)
        else  ! Find upwind-biased thickness near the surface.
          ! Perhaps this needs to be done more carefully, via find_eta.
          do i=is,ie ; if (do_i_shelf(i)) then
            zh(i) = 0.0 ; ztop_min(i) = min(zcol1(i), zcol2(i))
            i_htbl(i) = 1.0 / (visc%tbl_thick_shelf_v(i,j)*gv%Z_to_H + h_neglect)
          endif ; enddo
          do k=1,nz
            do i=is,ie ; if (do_i_shelf(i)) then
              zcol1(i) = zcol1(i) - h(i,j,k) ; zcol2(i) = zcol2(i) - h(i,j+1,k)
              zh(i) = zh(i) + h_harm(i,k)
 
              hvel_shelf(i,k) = hvel(i,k)
              if (v(i,j,k) * h_delta(i,k) > 0) then
                if (zh(i) * i_htbl(i) < cs%harm_BL_val) then
                  hvel_shelf(i,k) = min(hvel(i,k), h_harm(i,k))
                else
                  z2_wt = 1.0  ; if (zh(i) * i_htbl(i) < 2.0*cs%harm_BL_val) &
                    z2_wt = max(0.0, min(1.0, zh(i) * i_htbl(i) * i_valbl - 1.0))
                  z2 = z2_wt * (max(zh(i), ztop_min(i) - min(zcol1(i),zcol2(i))) * i_htbl(i))
                  topfn = 1.0 / (1.0 + 0.09*z2**6)
                  hvel_shelf(i,k) = min(hvel(i,k), (1.0-topfn)*h_arith(i,k) + topfn*h_harm(i,k))
                endif
             endif
            endif ; enddo
          enddo
          call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, h_harm, &
                                  bbl_thick, kv_bbl, z_i, h_ml, dt, j, g, gv, us, cs, &
                                  visc, forces, work_on_u=.false., obc=obc, shelf=.true.)
        endif
        do i=is,ie ; if (do_i_shelf(i)) cs%a1_shelf_v(i,j) = a_shelf(i,1) ; enddo
      endif
    endif
 
    if (do_any_shelf) then
      do k=1,nz+1 ; do i=is,ie ; if (do_i_shelf(i)) then
        cs%a_v(i,j,k) = forces%frac_shelf_v(i,j)  * a_shelf(i,k) + &
                   (1.0-forces%frac_shelf_v(i,j)) * a_cpl(i,k)
! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH
!        CS%a_v(i,J,K) = forces%frac_shelf_v(i,J)  * max(a_shelf(i,K), a_cpl(i,K)) + &
!                   (1.0-forces%frac_shelf_v(i,J)) * a_cpl(i,K)
      elseif (do_i(i)) then
        cs%a_v(i,j,k) = a_cpl(i,k)
      endif ; enddo ; enddo
      do k=1,nz ; do i=is,ie ; if (do_i_shelf(i)) then
        ! Should we instead take the inverse of the average of the inverses?
        cs%h_v(i,j,k) = forces%frac_shelf_v(i,j)  * hvel_shelf(i,k) + &
                   (1.0-forces%frac_shelf_v(i,j)) * hvel(i,k)
      elseif (do_i(i)) then
        cs%h_v(i,j,k) = hvel(i,k)
      endif ; enddo ; enddo
    else
      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
      do k=1,nz ; do i=is,ie ; if (do_i(i)) cs%h_v(i,j,k) = hvel(i,k) ; enddo ; enddo
    endif
 
    ! Diagnose total Kv at v-points
    if (cs%id_Kv_v > 0) then
      do k=1,nz ; do i=is,ie
        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)
      enddo ; enddo
    endif
 
  enddo ! end of v-point j loop
 
  if (cs%debug) then
    call uvchksum("vertvisc_coef h_[uv]", cs%h_u, cs%h_v, g%HI, haloshift=0, scale=gv%H_to_m)
    call uvchksum("vertvisc_coef a_[uv]", cs%a_u, cs%a_v, g%HI, haloshift=0, scale=us%Z_to_m*us%s_to_T)
    if (allocated(hml_u) .and. allocated(hml_v)) &
      call uvchksum("vertvisc_coef hML_[uv]", hml_u, hml_v, g%HI, haloshift=0, scale=gv%H_to_m)
  endif
 
! Offer diagnostic fields for averaging.
  if (associated(visc%Kv_slow) .and. (cs%id_Kv_slow > 0)) &
      call post_data(cs%id_Kv_slow, visc%Kv_slow, cs%diag)
  if (cs%id_Kv_u > 0) call post_data(cs%id_Kv_u, kv_u, cs%diag)
  if (cs%id_Kv_v > 0) call post_data(cs%id_Kv_v, kv_v, cs%diag)
  if (cs%id_au_vv > 0) call post_data(cs%id_au_vv, cs%a_u, cs%diag)
  if (cs%id_av_vv > 0) call post_data(cs%id_av_vv, cs%a_v, cs%diag)
  if (cs%id_h_u > 0) call post_data(cs%id_h_u, cs%h_u, cs%diag)
  if (cs%id_h_v > 0) call post_data(cs%id_h_v, cs%h_v, cs%diag)
  if (cs%id_hML_u > 0) call post_data(cs%id_hML_u, hml_u, cs%diag)
  if (cs%id_hML_v > 0) call post_data(cs%id_hML_v, hml_v, cs%diag)
 
  if (allocated(hml_u)) deallocate(hml_u)
  if (allocated(hml_v)) deallocate(hml_v)
 

References find_coupling_coef(), mom_error_handler::mom_error(), mom_open_boundary::obc_direction_n, mom_open_boundary::obc_direction_s, and mom_open_boundary::obc_direction_w.

Referenced by mom_dynamics_unsplit::step_mom_dyn_unsplit(), and mom_dynamics_unsplit_rk2::step_mom_dyn_unsplit_rk2().

Here is the call graph for this function:

Here is the caller graph for this function:

vertvisc_end()

subroutine, public mom_vert_friction::vertvisc_end

(

type(vertvisc_cs), pointer

CS

)

Clean up and deallocate the vertical friction module.

Parameters

cs	Vertical viscosity control structure that will be deallocated in this subroutine.

Definition at line 1803 of file MOM_vert_friction.F90.

  type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure that
                                   !! will be deallocated in this subroutine.
 
  dealloc_(cs%a_u) ; dealloc_(cs%h_u)
  dealloc_(cs%a_v) ; dealloc_(cs%h_v)
  if (associated(cs%a1_shelf_u)) deallocate(cs%a1_shelf_u)
  if (associated(cs%a1_shelf_v)) deallocate(cs%a1_shelf_v)
  deallocate(cs)

vertvisc_init()

subroutine, public mom_vert_friction::vertvisc_init	(	type(ocean_internal_state), intent(in), target	MIS,
		type(time_type), intent(in), target	Time,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(param_file_type), intent(in)	param_file,
		type(diag_ctrl), intent(inout), target	diag,
		type(accel_diag_ptrs), intent(inout)	ADp,
		type(directories), intent(in)	dirs,
		integer, intent(inout), target	ntrunc,
		type(vertvisc_cs), pointer	CS
	)

Initialize the vertical friction module.

Parameters

[in]	mis	The "MOM Internal State", a set of pointers
[in]	time	Current model time
[in]	g	Ocean grid structure
[in]	gv	Ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	param_file	File to parse for parameters
[in,out]	diag	Diagnostic control structure
[in,out]	adp	Acceleration diagnostic pointers
[in]	dirs	Relevant directory paths
[in,out]	ntrunc	Number of velocity truncations
	cs	Vertical viscosity control structure

Definition at line 1527 of file MOM_vert_friction.F90.

  type(ocean_internal_state), &
                   target, intent(in)    :: MIS    !< The "MOM Internal State", a set of pointers
                                                   !! to the fields and accelerations that make
                                                   !! up the ocean's physical state
  type(time_type), target, intent(in)    :: Time   !< Current model time
  type(ocean_grid_type),   intent(in)    :: G      !< Ocean grid structure
  type(verticalGrid_type), intent(in)    :: GV     !< Ocean vertical grid structure
  type(unit_scale_type),   intent(in)    :: US     !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file !< File to parse for parameters
  type(diag_ctrl), target, intent(inout) :: diag   !< Diagnostic control structure
  type(accel_diag_ptrs),   intent(inout) :: ADp    !< Acceleration diagnostic pointers
  type(directories),       intent(in)    :: dirs   !< Relevant directory paths
  integer, target,         intent(inout) :: ntrunc !< Number of velocity truncations
  type(vertvisc_CS),       pointer       :: CS     !< Vertical viscosity control structure
 
  ! Local variables
 
  real :: hmix_str_dflt
  real :: Kv_dflt ! A default viscosity [m2 s-1].
  real :: Hmix_m  ! A boundary layer thickness [m].
  logical :: default_2018_answers
  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_vert_friction" ! This module's name.
  character(len=40)  :: thickness_units
 
  if (associated(cs)) then
    call mom_error(warning, "vertvisc_init called with an associated "// &
                            "control structure.")
    return
  endif
  allocate(cs)
 
  if (gv%Boussinesq) then; thickness_units = "m"
  else; thickness_units = "kg m-2"; endif
 
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
 
  cs%diag => diag ; cs%ntrunc => ntrunc ; ntrunc = 0
 
! Default, read and log parameters
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
                 "This sets the default value for the various _2018_ANSWERS parameters.", &
                 default=.true.)
  call get_param(param_file, mdl, "VERT_FRICTION_2018_ANSWERS", cs%answers_2018, &
                 "If true, use the order of arithmetic and expressions that recover the answers "//&
                 "from the end of 2018.  Otherwise, use expressions that do not use an arbitary "//&
                 "and hard-coded maximum viscous coupling coefficient between layers.", &
                 default=default_2018_answers)
  call get_param(param_file, mdl, "BOTTOMDRAGLAW", cs%bottomdraglaw, &
                 "If true, the bottom stress is calculated with a drag "//&
                 "law of the form c_drag*|u|*u. The velocity magnitude "//&
                 "may be an assumed value or it may be based on the "//&
                 "actual velocity in the bottommost HBBL, depending on "//&
                 "LINEAR_DRAG.", default=.true.)
  call get_param(param_file, mdl, "CHANNEL_DRAG", cs%Channel_drag, &
                 "If true, the bottom drag is exerted directly on each "//&
                 "layer proportional to the fraction of the bottom it "//&
                 "overlies.", default=.false.)
  call get_param(param_file, mdl, "DIRECT_STRESS", cs%direct_stress, &
                 "If true, the wind stress is distributed over the "//&
                 "topmost HMIX_STRESS of fluid (like in HYCOM), and KVML "//&
                 "may be set to a very small value.", default=.false.)
  call get_param(param_file, mdl, "DYNAMIC_VISCOUS_ML", cs%dynamic_viscous_ML, &
                 "If true, use a bulk Richardson number criterion to "//&
                 "determine the mixed layer thickness for viscosity.", &
                 default=.false.)
  call get_param(param_file, mdl, "U_TRUNC_FILE", cs%u_trunc_file, &
                 "The absolute path to a file into which the accelerations "//&
                 "leading to zonal velocity truncations are written. "//&
                 "Undefine this for efficiency if this diagnostic is not "//&
                 "needed.", default=" ", debuggingparam=.true.)
  call get_param(param_file, mdl, "V_TRUNC_FILE", cs%v_trunc_file, &
                 "The absolute path to a file into which the accelerations "//&
                 "leading to meridional velocity truncations are written. "//&
                 "Undefine this for efficiency if this diagnostic is not "//&
                 "needed.", default=" ", debuggingparam=.true.)
  call get_param(param_file, mdl, "HARMONIC_VISC", cs%harmonic_visc, &
                 "If true, use the harmonic mean thicknesses for "//&
                 "calculating the vertical viscosity.", default=.false.)
  call get_param(param_file, mdl, "HARMONIC_BL_SCALE", cs%harm_BL_val, &
                 "A scale to determine when water is in the boundary "//&
                 "layers based solely on harmonic mean thicknesses for "//&
                 "the purpose of determining the extent to which the "//&
                 "thicknesses used in the viscosities are upwinded.", &
                 default=0.0, units="nondim")
  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
 
  if (gv%nkml < 1) &
    call get_param(param_file, mdl, "HMIX_FIXED", cs%Hmix, &
                 "The prescribed depth over which the near-surface "//&
                 "viscosity and diffusivity are elevated when the bulk "//&
                 "mixed layer is not used.", units="m", scale=gv%m_to_H, &
                 unscaled=hmix_m, fail_if_missing=.true.)
  if (cs%direct_stress) then
    if (gv%nkml < 1) then
      call get_param(param_file, mdl, "HMIX_STRESS", cs%Hmix_stress, &
                 "The depth over which the wind stress is applied if "//&
                 "DIRECT_STRESS is true.", units="m", default=hmix_m, scale=gv%m_to_H)
    else
      call get_param(param_file, mdl, "HMIX_STRESS", cs%Hmix_stress, &
                 "The depth over which the wind stress is applied if "//&
                 "DIRECT_STRESS is true.", units="m", fail_if_missing=.true., scale=gv%m_to_H)
    endif
    if (cs%Hmix_stress <= 0.0) call mom_error(fatal, "vertvisc_init: " // &
       "HMIX_STRESS must be set to a positive value if DIRECT_STRESS is true.")
  endif
  call get_param(param_file, mdl, "KV", cs%Kv, &
                 "The background kinematic viscosity in the interior. "//&
                 "The molecular value, ~1e-6 m2 s-1, may be used.", &
                 units="m2 s-1", fail_if_missing=.true., scale=us%m2_s_to_Z2_T, unscaled=kv_dflt)
 
  if (gv%nkml < 1) call get_param(param_file, mdl, "KVML", cs%Kvml, &
                 "The kinematic viscosity in the mixed layer.  A typical "//&
                 "value is ~1e-2 m2 s-1. KVML is not used if "//&
                 "BULKMIXEDLAYER is true.  The default is set by KV.", &
                 units="m2 s-1", default=kv_dflt, scale=us%m2_s_to_Z2_T)
  if (.not.cs%bottomdraglaw) call get_param(param_file, mdl, "KVBBL", cs%Kvbbl, &
                 "The kinematic viscosity in the benthic boundary layer. "//&
                 "A typical value is ~1e-2 m2 s-1. KVBBL is not used if "//&
                 "BOTTOMDRAGLAW is true.  The default is set by KV.", &
                 units="m2 s-1", default=kv_dflt, scale=us%m2_s_to_Z2_T)
  call get_param(param_file, mdl, "HBBL", cs%Hbbl, &
                 "The thickness of a bottom boundary layer with a "//&
                 "viscosity of KVBBL if BOTTOMDRAGLAW is not defined, or "//&
                 "the thickness over which near-bottom velocities are "//&
                 "averaged for the drag law if BOTTOMDRAGLAW is defined "//&
                 "but LINEAR_DRAG is not.", units="m", fail_if_missing=.true., scale=gv%m_to_H)
  call get_param(param_file, mdl, "MAXVEL", cs%maxvel, &
                 "The maximum velocity allowed before the velocity "//&
                 "components are truncated.", units="m s-1", default=3.0e8, scale=us%m_s_to_L_T)
  call get_param(param_file, mdl, "CFL_BASED_TRUNCATIONS", cs%CFL_based_trunc, &
                 "If true, base truncations on the CFL number, and not an "//&
                 "absolute speed.", default=.true.)
  call get_param(param_file, mdl, "CFL_TRUNCATE", cs%CFL_trunc, &
                 "The value of the CFL number that will cause velocity "//&
                 "components to be truncated; instability can occur past 0.5.", &
                 units="nondim", default=0.5)
  call get_param(param_file, mdl, "CFL_REPORT", cs%CFL_report, &
                 "The value of the CFL number that causes accelerations "//&
                 "to be reported; the default is CFL_TRUNCATE.", &
                 units="nondim", default=cs%CFL_trunc)
  call get_param(param_file, mdl, "CFL_TRUNCATE_RAMP_TIME", cs%truncRampTime, &
                 "The time over which the CFL truncation value is ramped "//&
                 "up at the beginning of the run.", &
                 units="s", default=0.)
  cs%CFL_truncE = cs%CFL_trunc
  call get_param(param_file, mdl, "CFL_TRUNCATE_START", cs%CFL_truncS, &
                 "The start value of the truncation CFL number used when "//&
                 "ramping up CFL_TRUNC.", &
                 units="nondim", default=0.)
  call get_param(param_file, mdl, "STOKES_MIXING_COMBINED", cs%StokesMixing, &
                 "Flag to use Stokes drift Mixing via the Lagrangian "//&
                 " current (Eulerian plus Stokes drift). "//&
                 " Still needs work and testing, so not recommended for use.",&
                 default=.false.)
  !BGR 04/04/2018{
  ! StokesMixing is required for MOM6 for some Langmuir mixing parameterization.
  !   The code used here has not been developed for vanishing layers or in
  !   conjunction with any bottom friction.  Therefore, the following line is
  !   added so this functionality cannot be used without user intervention in
  !   the code.  This will prevent general use of this functionality until proper
  !   care is given to the previously mentioned issues.  Comment out the following
  !   MOM_error to use, but do so at your own risk and with these points in mind.
  !}
  if (cs%StokesMixing) then
    call mom_error(fatal, "Stokes mixing requires user intervention in the code.\n"//&
                          "  Model now exiting.  See MOM_vert_friction.F90 for \n"//&
                          "  details (search 'BGR 04/04/2018' to locate comment).")
  endif
  call get_param(param_file, mdl, "VEL_UNDERFLOW", cs%vel_underflow, &
                 "A negligibly small velocity magnitude below which velocity "//&
                 "components are set to 0.  A reasonable value might be "//&
                 "1e-30 m/s, which is less than an Angstrom divided by "//&
                 "the age of the universe.", units="m s-1", default=0.0, scale=us%m_s_to_L_T)
 
  alloc_(cs%a_u(isdb:iedb,jsd:jed,nz+1)) ; cs%a_u(:,:,:) = 0.0
  alloc_(cs%h_u(isdb:iedb,jsd:jed,nz))   ; cs%h_u(:,:,:) = 0.0
  alloc_(cs%a_v(isd:ied,jsdb:jedb,nz+1)) ; cs%a_v(:,:,:) = 0.0
  alloc_(cs%h_v(isd:ied,jsdb:jedb,nz))   ; cs%h_v(:,:,:) = 0.0
 
  cs%id_Kv_slow = register_diag_field('ocean_model', 'Kv_slow', diag%axesTi, time, &
     'Slow varying vertical viscosity', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
 
  cs%id_Kv_u = register_diag_field('ocean_model', 'Kv_u', diag%axesCuL, time, &
     'Total vertical viscosity at u-points', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
 
  cs%id_Kv_v = register_diag_field('ocean_model', 'Kv_v', diag%axesCvL, time, &
     'Total vertical viscosity at v-points', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
 
  cs%id_au_vv = register_diag_field('ocean_model', 'au_visc', diag%axesCui, time, &
     'Zonal Viscous Vertical Coupling Coefficient', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
 
  cs%id_av_vv = register_diag_field('ocean_model', 'av_visc', diag%axesCvi, time, &
     'Meridional Viscous Vertical Coupling Coefficient', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
 
  cs%id_h_u = register_diag_field('ocean_model', 'Hu_visc', diag%axesCuL, time, &
     'Thickness at Zonal Velocity Points for Viscosity', thickness_units, &
     conversion=gv%H_to_m)
 
  cs%id_h_v = register_diag_field('ocean_model', 'Hv_visc', diag%axesCvL, time, &
     'Thickness at Meridional Velocity Points for Viscosity', thickness_units, &
     conversion=gv%H_to_m)
 
  cs%id_hML_u = register_diag_field('ocean_model', 'HMLu_visc', diag%axesCu1, time, &
     'Mixed Layer Thickness at Zonal Velocity Points for Viscosity', thickness_units, &
     conversion=gv%H_to_m)
 
  cs%id_hML_v = register_diag_field('ocean_model', 'HMLv_visc', diag%axesCv1, time, &
     'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', thickness_units, &
     conversion=gv%H_to_m)
 
  cs%id_du_dt_visc = register_diag_field('ocean_model', 'du_dt_visc', diag%axesCuL, &
     time, 'Zonal Acceleration from Vertical Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)
  if (cs%id_du_dt_visc > 0) call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)
  cs%id_dv_dt_visc = register_diag_field('ocean_model', 'dv_dt_visc', diag%axesCvL, &
     time, 'Meridional Acceleration from Vertical Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)
  if (cs%id_dv_dt_visc > 0) call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)
 
  cs%id_taux_bot = register_diag_field('ocean_model', 'taux_bot', diag%axesCu1, &
     time, 'Zonal Bottom Stress from Ocean to Earth', 'Pa', &
     conversion=us%R_to_kg_m3*us%L_T2_to_m_s2*us%Z_to_m)
  cs%id_tauy_bot = register_diag_field('ocean_model', 'tauy_bot', diag%axesCv1, &
     time, 'Meridional Bottom Stress from Ocean to Earth', 'Pa', &
     conversion=us%R_to_kg_m3*us%L_T2_to_m_s2*us%Z_to_m)
 
  if ((len_trim(cs%u_trunc_file) > 0) .or. (len_trim(cs%v_trunc_file) > 0)) &
    call pointaccel_init(mis, time, g, param_file, diag, dirs, cs%PointAccel_CSp)
 

References mom_error_handler::mom_error().

Referenced by mom_dynamics_unsplit::initialize_dyn_unsplit(), and mom_dynamics_unsplit_rk2::initialize_dyn_unsplit_rk2().

Here is the call graph for this function:

Here is the caller graph for this function:

vertvisc_limit_vel()

subroutine, public mom_vert_friction::vertvisc_limit_vel	(	real, dimension(szib_(g),szj_(g),szk_(gv)), intent(inout)	u,
		real, dimension(szi_(g),szjb_(g),szk_(gv)), intent(inout)	v,
		real, dimension(szi_(g),szj_(g),szk_(gv)), intent(in)	h,
		type(accel_diag_ptrs), intent(in)	ADp,
		type(cont_diag_ptrs), intent(in)	CDp,
		type(mech_forcing), intent(in)	forces,
		type(vertvisc_type), intent(in)	visc,
		real, intent(in)	dt,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(vertvisc_cs), pointer	CS
	)

Velocity components which exceed a threshold for physically reasonable values are truncated. Optionally, any column with excessive velocities may be sent to a diagnostic reporting subroutine.

Parameters

[in]	g	Ocean grid structure
[in]	gv	Ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in,out]	u	Zonal velocity [L T-1 ~> m s-1]
[in,out]	v	Meridional velocity [L T-1 ~> m s-1]
[in]	h	Layer thickness [H ~> m or kg m-2]
[in]	adp	Acceleration diagnostic pointers
[in]	cdp	Continuity diagnostic pointers
[in]	forces	A structure with the driving mechanical forces
[in]	visc	Viscosities and bottom drag
[in]	dt	Time increment [T ~> s]
	cs	Vertical viscosity control structure

Definition at line 1316 of file MOM_vert_friction.F90.

  type(ocean_grid_type),   intent(in)    :: G      !< Ocean grid structure
  type(verticalGrid_type), intent(in)    :: GV     !< Ocean vertical grid structure
  type(unit_scale_type),   intent(in)    :: US     !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: u      !< Zonal velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
                           intent(inout) :: v      !< Meridional velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: h      !< Layer thickness [H ~> m or kg m-2]
  type(accel_diag_ptrs),   intent(in)    :: ADp    !< Acceleration diagnostic pointers
  type(cont_diag_ptrs),    intent(in)    :: CDp    !< Continuity diagnostic pointers
  type(mech_forcing),      intent(in)    :: forces !< A structure with the driving mechanical forces
  type(vertvisc_type),     intent(in)    :: visc   !< Viscosities and bottom drag
  real,                    intent(in)    :: dt     !< Time increment [T ~> s]
  type(vertvisc_CS),       pointer       :: CS     !< Vertical viscosity control structure
 
  ! Local variables
 
  real :: maxvel           ! Velocities components greater than maxvel
  real :: truncvel         ! are truncated to truncvel, both [L T-1 ~> m s-1].
  real :: CFL              ! The local CFL number.
  real :: H_report         ! A thickness below which not to report truncations.
  real :: dt_Rho0          ! The timestep divided by the Boussinesq density [m2 T2 s-1 L-1 Z-1 R-1 ~> s m3 kg-1].
  real :: vel_report(SZIB_(G),SZJB_(G))   ! The velocity to report [L T-1 ~> m s-1]
  real :: u_old(SZIB_(G),SZJ_(G),SZK_(G)) ! The previous u-velocity [L T-1 ~> m s-1]
  real :: v_old(SZI_(G),SZJB_(G),SZK_(G)) ! The previous v-velocity [L T-1 ~> m s-1]
  logical :: trunc_any, dowrite(SZIB_(G),SZJB_(G))
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
 
  maxvel = cs%maxvel
  truncvel = 0.9*maxvel
  h_report = 6.0 * gv%Angstrom_H
  dt_rho0 = (us%L_T_to_m_s*us%Z_to_m) * dt / gv%Rho0
 
  if (len_trim(cs%u_trunc_file) > 0) then
    !$OMP parallel do default(shared) private(trunc_any,CFL)
    do j=js,je
      trunc_any = .false.
      do i=isq,ieq ; dowrite(i,j) = .false. ; enddo
      if (cs%CFL_based_trunc) then
        do i=isq,ieq ; vel_report(i,j) = 3.0e8*us%m_s_to_L_T ; enddo ! Speed of light default.
        do k=1,nz ; do i=isq,ieq
          if (abs(u(i,j,k)) < cs%vel_underflow) u(i,j,k) = 0.0
          if (u(i,j,k) < 0.0) then
            cfl = (-u(i,j,k) * dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
          else
            cfl = (u(i,j,k) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
          endif
          if (cfl > cs%CFL_trunc) trunc_any = .true.
          if (cfl > cs%CFL_report) then
            dowrite(i,j) = .true.
            vel_report(i,j) = min(vel_report(i,j), abs(u(i,j,k)))
          endif
        enddo ; enddo
      else
        do i=isq,ieq; vel_report(i,j) = maxvel; enddo
        do k=1,nz ; do i=isq,ieq
          if (abs(u(i,j,k)) < cs%vel_underflow) then ; u(i,j,k) = 0.0
          elseif (abs(u(i,j,k)) > maxvel) then
            dowrite(i,j) = .true. ; trunc_any = .true.
          endif
        enddo ; enddo
      endif
 
      do i=isq,ieq ; if (dowrite(i,j)) then
        u_old(i,j,:) = u(i,j,:)
      endif ; enddo
 
      if (trunc_any) then ; if (cs%CFL_based_trunc) then
        do k=1,nz ; do i=isq,ieq
          if ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i+1,j) < -cs%CFL_trunc) then
            u(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i+1,j) / (dt * g%dy_Cu(i,j)))
            if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
          elseif ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then
            u(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dy_Cu(i,j)))
            if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
          endif
        enddo ; enddo
      else
        do k=1,nz ; do i=isq,ieq ; if (abs(u(i,j,k)) > maxvel) then
          u(i,j,k) = sign(truncvel,u(i,j,k))
          if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
        endif ; enddo ; enddo
      endif ; endif
    enddo ! j-loop
  else  ! Do not report accelerations leading to large velocities.
    if (cs%CFL_based_trunc) then
!$OMP parallel do default(none) shared(nz,js,je,Isq,Ieq,u,dt,G,CS,h,H_report)
      do k=1,nz ; do j=js,je ; do i=isq,ieq
        if (abs(u(i,j,k)) < cs%vel_underflow) then ; u(i,j,k) = 0.0
        elseif ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i+1,j) < -cs%CFL_trunc) then
          u(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i+1,j) / (dt * g%dy_Cu(i,j)))
          if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
        elseif ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then
          u(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dy_Cu(i,j)))
          if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
        endif
      enddo ; enddo ; enddo
    else
!$OMP parallel do default(none) shared(nz,js,je,Isq,Ieq,u,G,CS,truncvel,maxvel,h,H_report)
      do k=1,nz ; do j=js,je ; do i=isq,ieq
        if (abs(u(i,j,k)) < cs%vel_underflow) then ; u(i,j,k) = 0.0
        elseif (abs(u(i,j,k)) > maxvel) then
          u(i,j,k) = sign(truncvel, u(i,j,k))
          if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
        endif
      enddo ; enddo ; enddo
    endif
  endif
 
  if (len_trim(cs%u_trunc_file) > 0) then
    do j=js,je; do i=isq,ieq ; if (dowrite(i,j)) then
!   Here the diagnostic reporting subroutines are called if
! unphysically large values were found.
      call write_u_accel(i, j, u_old, h, adp, cdp, dt, g, gv, us, cs%PointAccel_CSp, &
               vel_report(i,j), forces%taux(i,j)*dt_rho0, a=cs%a_u, hv=cs%h_u)
    endif ; enddo ; enddo
  endif
 
  if (len_trim(cs%v_trunc_file) > 0) then
    !$OMP parallel do default(shared) private(trunc_any,CFL)
    do j=jsq,jeq
      trunc_any = .false.
      do i=is,ie ; dowrite(i,j) = .false. ; enddo
      if (cs%CFL_based_trunc) then
        do i=is,ie ; vel_report(i,j) = 3.0e8*us%m_s_to_L_T ; enddo ! Speed of light default.
        do k=1,nz ; do i=is,ie
          if (abs(v(i,j,k)) < cs%vel_underflow) v(i,j,k) = 0.0
          if (v(i,j,k) < 0.0) then
            cfl = (-v(i,j,k) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
          else
            cfl = (v(i,j,k) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
          endif
          if (cfl > cs%CFL_trunc) trunc_any = .true.
          if (cfl > cs%CFL_report) then
            dowrite(i,j) = .true.
            vel_report(i,j) = min(vel_report(i,j), abs(v(i,j,k)))
          endif
        enddo ; enddo
      else
        do i=is,ie ; vel_report(i,j) = maxvel ; enddo
        do k=1,nz ; do i=is,ie
          if (abs(v(i,j,k)) < cs%vel_underflow) then ; v(i,j,k) = 0.0
          elseif (abs(v(i,j,k)) > maxvel) then
            dowrite(i,j) = .true. ; trunc_any = .true.
          endif
        enddo ; enddo
      endif
 
      do i=is,ie ; if (dowrite(i,j)) then
        v_old(i,j,:) = v(i,j,:)
      endif ; enddo
 
      if (trunc_any) then ; if (cs%CFL_based_trunc) then
        do k=1,nz; do i=is,ie
          if ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j+1) < -cs%CFL_trunc) then
            v(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i,j+1) / (dt * g%dx_Cv(i,j)))
            if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
          elseif ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then
            v(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dx_Cv(i,j)))
            if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
          endif
        enddo ; enddo
      else
        do k=1,nz ; do i=is,ie ; if (abs(v(i,j,k)) > maxvel) then
          v(i,j,k) = sign(truncvel,v(i,j,k))
          if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
        endif ; enddo ; enddo
      endif ; endif
    enddo ! J-loop
  else  ! Do not report accelerations leading to large velocities.
    if (cs%CFL_based_trunc) then
      !$OMP parallel do default(shared)
      do k=1,nz ; do j=jsq,jeq ; do i=is,ie
        if (abs(v(i,j,k)) < cs%vel_underflow) then ; v(i,j,k) = 0.0
        elseif ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j+1) < -cs%CFL_trunc) then
          v(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i,j+1) / (dt * g%dx_Cv(i,j)))
          if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
        elseif ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then
          v(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dx_Cv(i,j)))
          if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
        endif
      enddo ; enddo ; enddo
    else
      !$OMP parallel do default(shared)
      do k=1,nz ; do j=jsq,jeq ; do i=is,ie
        if (abs(v(i,j,k)) < cs%vel_underflow) then ; v(i,j,k) = 0.0
        elseif (abs(v(i,j,k)) > maxvel) then
          v(i,j,k) = sign(truncvel, v(i,j,k))
          if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
        endif
      enddo ; enddo ; enddo
    endif
  endif
 
  if (len_trim(cs%v_trunc_file) > 0) then
    do j=jsq,jeq; do i=is,ie ; if (dowrite(i,j)) then
!   Here the diagnostic reporting subroutines are called if
! unphysically large values were found.
      call write_v_accel(i, j, v_old, h, adp, cdp, dt, g, gv, us, cs%PointAccel_CSp, &
               vel_report(i,j), forces%tauy(i,j)*dt_rho0, a=cs%a_v, hv=cs%h_v)
    endif ; enddo ; enddo
  endif
 

Referenced by vertvisc().

Here is the caller graph for this function:

vertvisc_remnant()

subroutine, public mom_vert_friction::vertvisc_remnant	(	type(vertvisc_type), intent(in)	visc,
		real, dimension(szib_(g),szj_(g),szk_(gv)), intent(inout)	visc_rem_u,
		real, dimension( g %isd: g %ied, g %jsdb: g %jedb, gv %ke), intent(inout)	visc_rem_v,
		real, intent(in)	dt,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(vertvisc_cs), pointer	CS
	)

Calculate the fraction of momentum originally in a layer that remains after a time-step of viscosity, and the fraction of a time-step's worth of barotropic acceleration that a layer experiences after viscosity is applied.

Parameters

[in]	g	Ocean grid structure
[in]	gv	Ocean vertical grid structure
[in]	visc	Viscosities and bottom drag
[in,out]	visc_rem_u	Fraction of a time-step's worth of a
[in,out]	visc_rem_v	Fraction of a time-step's worth of a
[in]	dt	Time increment [T ~> s]
[in]	us	A dimensional unit scaling type
	cs	Vertical viscosity control structure

Definition at line 463 of file MOM_vert_friction.F90.

  type(ocean_grid_type), intent(in)   :: G    !< Ocean grid structure
  type(verticalGrid_type), intent(in) :: GV   !< Ocean vertical grid structure
  type(vertvisc_type),   intent(in)   :: visc !< Viscosities and bottom drag
  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
                         intent(inout) :: visc_rem_u !< Fraction of a time-step's worth of a
                                              !! barotopic acceleration that a layer experiences after
                                              !! viscosity is applied in the zonal direction [nondim]
  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
                         intent(inout) :: visc_rem_v !< Fraction of a time-step's worth of a
                                              !! barotopic acceleration that a layer experiences after
                                              !! viscosity is applied in the meridional direction [nondim]
  real,                  intent(in)    :: dt  !< Time increment [T ~> s]
  type(unit_scale_type), intent(in)    :: US  !< A dimensional unit scaling type
  type(vertvisc_CS),     pointer       :: CS  !< Vertical viscosity control structure
 
  ! Local variables
 
  real :: b1(SZIB_(G))          ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
  real :: c1(SZIB_(G),SZK_(G))  ! A variable used by the tridiagonal solver [nondim].
  real :: d1(SZIB_(G))          ! d1=1-c1 is used by the tridiagonal solver [nondim].
  real :: Ray(SZIB_(G),SZK_(G)) ! Ray is the Rayleigh-drag velocity [Z T-1 ~> m s-1].
  real :: b_denom_1   ! The first term in the denominator of b1 [H ~> m or kg m-2].
  real :: dt_Z_to_H        ! The time step times the conversion from Z to the
                           ! units of thickness [T H Z-1 ~> s or s kg m-3].
  logical :: do_i(SZIB_(G))
 
  integer :: i, j, k, is, ie, Isq, Ieq, Jsq, Jeq, nz
  is = g%isc ; ie = g%iec
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
 
  if (.not.associated(cs)) call mom_error(fatal,"MOM_vert_friction(visc): "// &
         "Module must be initialized before it is used.")
 
  dt_z_to_h = dt*gv%Z_to_H
 
  do k=1,nz ; do i=isq,ieq ; ray(i,k) = 0.0 ; enddo ; enddo
 
  ! Find the zonal viscous using a modification of a standard tridagonal solver.
!$OMP parallel do default(none) shared(G,Isq,Ieq,CS,nz,visc,dt_Z_to_H,visc_rem_u) &
!$OMP                     firstprivate(Ray)                                       &
!$OMP                          private(do_i,b_denom_1,b1,d1,c1)
  do j=g%jsc,g%jec
    do i=isq,ieq ; do_i(i) = (g%mask2dCu(i,j) > 0) ; enddo
 
    if (cs%Channel_drag) then ; do k=1,nz ; do i=isq,ieq
      ray(i,k) = visc%Ray_u(i,j,k)
    enddo ; enddo ; endif
 
    do i=isq,ieq ; if (do_i(i)) then
      b_denom_1 = cs%h_u(i,j,1) + dt_z_to_h * (ray(i,1) + cs%a_u(i,j,1))
      b1(i) = 1.0 / (b_denom_1 + dt_z_to_h*cs%a_u(i,j,2))
      d1(i) = b_denom_1 * b1(i)
      visc_rem_u(i,j,1) = b1(i) * cs%h_u(i,j,1)
    endif ; enddo
    do k=2,nz ; do i=isq,ieq ; if (do_i(i)) then
      c1(i,k) = dt_z_to_h * cs%a_u(i,j,k)*b1(i)
      b_denom_1 = cs%h_u(i,j,k) + dt_z_to_h * (ray(i,k) + cs%a_u(i,j,k)*d1(i))
      b1(i) = 1.0 / (b_denom_1 + dt_z_to_h * cs%a_u(i,j,k+1))
      d1(i) = b_denom_1 * b1(i)
      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)
    endif ; enddo ; enddo
    do k=nz-1,1,-1 ; do i=isq,ieq ; if (do_i(i)) then
      visc_rem_u(i,j,k) = visc_rem_u(i,j,k) + c1(i,k+1)*visc_rem_u(i,j,k+1)
 
    endif ; enddo ; enddo ! i and k loops
 
  enddo ! end u-component j loop
 
  ! Now find the meridional viscous using a modification.
!$OMP parallel do default(none) shared(Jsq,Jeq,is,ie,G,CS,visc,dt_Z_to_H,visc_rem_v,nz) &
!$OMP                     firstprivate(Ray)                                             &
!$OMP                          private(do_i,b_denom_1,b1,d1,c1)
  do j=jsq,jeq
    do i=is,ie ; do_i(i) = (g%mask2dCv(i,j) > 0) ; enddo
 
    if (cs%Channel_drag) then ; do k=1,nz ; do i=is,ie
      ray(i,k) = visc%Ray_v(i,j,k)
    enddo ; enddo ; endif
 
    do i=is,ie ; if (do_i(i)) then
      b_denom_1 = cs%h_v(i,j,1) + dt_z_to_h * (ray(i,1) + cs%a_v(i,j,1))
      b1(i) = 1.0 / (b_denom_1 + dt_z_to_h*cs%a_v(i,j,2))
      d1(i) = b_denom_1 * b1(i)
      visc_rem_v(i,j,1) = b1(i) * cs%h_v(i,j,1)
    endif ; enddo
    do k=2,nz ; do i=is,ie ; if (do_i(i)) then
      c1(i,k) = dt_z_to_h * cs%a_v(i,j,k)*b1(i)
      b_denom_1 = cs%h_v(i,j,k) + dt_z_to_h * (ray(i,k) + cs%a_v(i,j,k)*d1(i))
      b1(i) = 1.0 / (b_denom_1 + dt_z_to_h * cs%a_v(i,j,k+1))
      d1(i) = b_denom_1 * b1(i)
      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)
    endif ; enddo ; enddo
    do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then
      visc_rem_v(i,j,k) = visc_rem_v(i,j,k) + c1(i,k+1)*visc_rem_v(i,j,k+1)
    endif ; enddo ; enddo ! i and k loops
  enddo ! end of v-component J loop
 
  if (cs%debug) then
    call uvchksum("visc_rem_[uv]", visc_rem_u, visc_rem_v, g%HI,haloshift=0)
  endif
 

References mom_error_handler::mom_error().

Here is the call graph for this function: