MOM6
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]gOcean grid structure
[in]gvOcean vertical grid structure
[in]usA dimensional unit scaling type
[out]a_cplCoupling coefficient across interfaces [Z T-1 ~> m s-1].
[in]hvelThickness at velocity points [H ~> m or kg m-2]
[in]do_iIf true, determine coupling coefficient for a column
[in]h_harmHarmonic mean of thicknesses around a velocity
[in]bbl_thickBottom boundary layer thickness [H ~> m or kg m-2]
[in]kv_bblBottom boundary layer viscosity [Z2 T-1 ~> m2 s-1].
[in]z_iEstimate of interface heights above the bottom,
[out]h_mlMixed layer depth [H ~> m or kg m-2]
[in]jj-index to find coupling coefficient for
[in]dtTime increment [s]
csVertical viscosity control structure
[in]viscStructure containing viscosities and bottom drag
[in]forcesA structure with the driving mechanical forces
[in]work_on_uIf true, u-points are being calculated, otherwise they are v-points
obcOpen boundary condition structure
[in]shelfIf present and true, use a surface boundary condition appropriate for an ice shelf.

Definition at line 1037 of file MOM_vert_friction.F90.

1037  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1038  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1039  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1040  real, dimension(SZIB_(G),SZK_(GV)+1), &
1041  intent(out) :: a_cpl !< Coupling coefficient across interfaces [Z T-1 ~> m s-1].
1042  real, dimension(SZIB_(G),SZK_(GV)), &
1043  intent(in) :: hvel !< Thickness at velocity points [H ~> m or kg m-2]
1044  logical, dimension(SZIB_(G)), &
1045  intent(in) :: do_i !< If true, determine coupling coefficient for a column
1046  real, dimension(SZIB_(G),SZK_(GV)), &
1047  intent(in) :: h_harm !< Harmonic mean of thicknesses around a velocity
1048  !! grid point [H ~> m or kg m-2]
1049  real, dimension(SZIB_(G)), intent(in) :: bbl_thick !< Bottom boundary layer thickness [H ~> m or kg m-2]
1050  real, dimension(SZIB_(G)), intent(in) :: kv_bbl !< Bottom boundary layer viscosity [Z2 T-1 ~> m2 s-1].
1051  real, dimension(SZIB_(G),SZK_(GV)+1), &
1052  intent(in) :: z_i !< Estimate of interface heights above the bottom,
1053  !! normalized by the bottom boundary layer thickness
1054  real, dimension(SZIB_(G)), intent(out) :: h_ml !< Mixed layer depth [H ~> m or kg m-2]
1055  integer, intent(in) :: j !< j-index to find coupling coefficient for
1056  real, intent(in) :: dt !< Time increment [s]
1057  type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
1058  type(vertvisc_type), intent(in) :: visc !< Structure containing viscosities and bottom drag
1059  type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
1060  logical, intent(in) :: work_on_u !< If true, u-points are being calculated,
1061  !! otherwise they are v-points
1062  type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure
1063  logical, optional, intent(in) :: shelf !< If present and true, use a surface boundary
1064  !! condition appropriate for an ice shelf.
1065 
1066  ! Local variables
1067 
1068  real, dimension(SZIB_(G)) :: &
1069  u_star, & ! ustar at a velocity point [Z T-1 ~> m s-1].
1070  absf, & ! The average of the neighboring absolute values of f [T-1 ~> s-1].
1071 ! h_ml, & ! The mixed layer depth [H ~> m or kg m-2].
1072  nk_visc, & ! The (real) interface index of the base of mixed layer.
1073  z_t, & ! The distance from the top, sometimes normalized
1074  ! by Hmix, [H ~> m or kg m-2] or [nondim].
1075  kv_tbl, & ! The viscosity in a top boundary layer under ice [Z2 T-1 ~> m2 s-1].
1076  tbl_thick
1077  real, dimension(SZIB_(G),SZK_(GV)+1) :: &
1078  Kv_tot, & ! The total viscosity at an interface [Z2 T-1 ~> m2 s-1].
1079  Kv_add ! A viscosity to add [Z2 T-1 ~> m2 s-1].
1080  real :: h_shear ! The distance over which shears occur [H ~> m or kg m-2].
1081  real :: r ! A thickness to compare with Hbbl [H ~> m or kg m-2].
1082  real :: visc_ml ! The mixed layer viscosity [Z2 T-1 ~> m2 s-1].
1083  real :: I_Hmix ! The inverse of the mixed layer thickness [H-1 ~> m-1 or m2 kg-1].
1084  real :: a_ml ! The layer coupling coefficient across an interface in
1085  ! the mixed layer [Z T-1 ~> m s-1].
1086  real :: I_amax ! The inverse of the maximum coupling coefficient [T s-1 Z-1 ~> m-1].???
1087  real :: temp1 ! A temporary variable [H Z ~> m2 or kg m-1]
1088  real :: h_neglect ! A thickness that is so small it is usually lost
1089  ! in roundoff and can be neglected [H ~> m or kg m-2].
1090  real :: z2 ! A copy of z_i [nondim]
1091  real :: topfn ! A function that is 1 at the top and small far from it [nondim]
1092  real :: a_top ! A viscosity associated with the top boundary layer [Z2 T-1 ~> m2 s-1]
1093  logical :: do_shelf, do_OBCs
1094  integer :: i, k, is, ie, max_nk
1095  integer :: nz
1096  real :: botfn
1097 
1098  a_cpl(:,:) = 0.0
1099  kv_tot(:,:) = 0.0
1100 
1101  if (work_on_u) then ; is = g%IscB ; ie = g%IecB
1102  else ; is = g%isc ; ie = g%iec ; endif
1103  nz = g%ke
1104  h_neglect = gv%H_subroundoff
1105 
1106  ! The maximum coupling coefficent was originally introduced to avoid
1107  ! truncation error problems in the tridiagonal solver. Effectively, the 1e-10
1108  ! sets the maximum coupling coefficient increment to 1e10 m per timestep.
1109  i_amax = (1.0e-10*us%Z_to_m) * dt*us%s_to_T
1110 
1111  do_shelf = .false. ; if (present(shelf)) do_shelf = shelf
1112  do_obcs = .false.
1113  if (associated(obc)) then ; do_obcs = (obc%number_of_segments > 0) ; endif
1114  h_ml(:) = 0.0
1115 
1116 ! The following loop calculates the vertical average velocity and
1117 ! surface mixed layer contributions to the vertical viscosity.
1118  do i=is,ie ; kv_tot(i,1) = 0.0 ; enddo
1119  if ((gv%nkml>0) .or. do_shelf) then ; do k=2,nz ; do i=is,ie
1120  if (do_i(i)) kv_tot(i,k) = cs%Kv
1121  enddo ; enddo ; else
1122  i_hmix = 1.0 / (cs%Hmix + h_neglect)
1123  do i=is,ie ; z_t(i) = h_neglect*i_hmix ; enddo
1124  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1125  z_t(i) = z_t(i) + h_harm(i,k-1)*i_hmix
1126  kv_tot(i,k) = cs%Kv + cs%Kvml / ((z_t(i)*z_t(i)) * &
1127  (1.0 + 0.09*z_t(i)*z_t(i)*z_t(i)*z_t(i)*z_t(i)*z_t(i)))
1128  endif ; enddo ; enddo
1129  endif
1130 
1131  do i=is,ie ; if (do_i(i)) then
1132  if (cs%bottomdraglaw) then
1133  r = hvel(i,nz)*0.5
1134  if (r < bbl_thick(i)) then
1135  a_cpl(i,nz+1) = kv_bbl(i) / (i_amax*kv_bbl(i) + r*gv%H_to_Z)
1136  else
1137  a_cpl(i,nz+1) = kv_bbl(i) / (i_amax*kv_bbl(i) + bbl_thick(i)*gv%H_to_Z)
1138  endif
1139  else
1140  a_cpl(i,nz+1) = cs%Kvbbl / (0.5*hvel(i,nz)*gv%H_to_Z + i_amax*cs%Kvbbl)
1141  endif
1142  endif ; enddo
1143 
1144  if (associated(visc%Kv_shear)) then
1145  ! The factor of 2 that used to be required in the viscosities is no longer needed.
1146  if (work_on_u) then
1147  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1148  kv_add(i,k) = 0.5*(visc%Kv_shear(i,j,k) + visc%Kv_shear(i+1,j,k))
1149  endif ; enddo ; enddo
1150  if (do_obcs) then
1151  do i=is,ie ; if (do_i(i) .and. (obc%segnum_u(i,j) /= obc_none)) then
1152  if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
1153  do k=2,nz ; kv_add(i,k) = visc%Kv_shear(i,j,k) ; enddo
1154  elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
1155  do k=2,nz ; kv_add(i,k) = visc%Kv_shear(i+1,j,k) ; enddo
1156  endif
1157  endif ; enddo
1158  endif
1159  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1160  kv_tot(i,k) = kv_tot(i,k) + kv_add(i,k)
1161  endif ; enddo ; enddo
1162  else
1163  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1164  kv_add(i,k) = 0.5*(visc%Kv_shear(i,j,k) + visc%Kv_shear(i,j+1,k))
1165  endif ; enddo ; enddo
1166  if (do_obcs) then
1167  do i=is,ie ; if (do_i(i) .and. (obc%segnum_v(i,j) /= obc_none)) then
1168  if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
1169  do k=2,nz ; kv_add(i,k) = visc%Kv_shear(i,j,k) ; enddo
1170  elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
1171  do k=2,nz ; kv_add(i,k) = visc%Kv_shear(i,j+1,k) ; enddo
1172  endif
1173  endif ; enddo
1174  endif
1175  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1176  kv_tot(i,k) = kv_tot(i,k) + kv_add(i,k)
1177  endif ; enddo ; enddo
1178  endif
1179  endif
1180 
1181  if (associated(visc%Kv_shear_Bu)) then
1182  if (work_on_u) then
1183  do k=2,nz ; do i=is,ie ; If (do_i(i)) then
1184  kv_tot(i,k) = kv_tot(i,k) + (0.5)*(visc%Kv_shear_Bu(i,j-1,k) + visc%Kv_shear_Bu(i,j,k))
1185  endif ; enddo ; enddo
1186  else
1187  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1188  kv_tot(i,k) = kv_tot(i,k) + (0.5)*(visc%Kv_shear_Bu(i-1,j,k) + visc%Kv_shear_Bu(i,j,k))
1189  endif ; enddo ; enddo
1190  endif
1191  endif
1192 
1193  ! add "slow" varying vertical viscosity (e.g., from background, tidal etc)
1194  if (associated(visc%Kv_slow) .and. (visc%add_Kv_slow)) then
1195  ! GMM/ A factor of 2 is also needed here, see comment above from BGR.
1196  if (work_on_u) then
1197  !### Incrementing Kv_add here will cause visc%Kv_shear to be double counted. - RWH
1198  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1199  kv_add(i,k) = kv_add(i,k) + 0.5 * (visc%Kv_slow(i,j,k) + visc%Kv_slow(i+1,j,k))
1200  ! Should be : Kv_add(I,K) = 0.5 * (visc%Kv_slow(i,j,k) + visc%Kv_slow(i+1,j,k))
1201  endif ; enddo ; enddo
1202  !### I am pretty sure that this code is double counting viscosity at OBC points! - RWH
1203  if (do_obcs) then
1204  do i=is,ie ; if (do_i(i) .and. (obc%segnum_u(i,j) /= obc_none)) then
1205  if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
1206  do k=2,nz ; kv_add(i,k) = kv_add(i,k) + visc%Kv_slow(i,j,k) ; enddo
1207  ! Should be : do K=2,nz ; Kv_add(I,K) = visc%Kv_slow(i,j,k) ; enddo
1208  elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
1209  do k=2,nz ; kv_add(i,k) = kv_add(i,k) + visc%Kv_slow(i+1,j,k) ; enddo
1210  ! Should be : do K=2,nz ; Kv_add(I,K) = visc%Kv_slow(i+1,j,k) ; enddo
1211  endif
1212  endif ; enddo
1213  endif
1214  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1215  kv_tot(i,k) = kv_tot(i,k) + kv_add(i,k)
1216  endif ; enddo ; enddo
1217  else
1218  !### Incrementing Kv_add here will cause visc%Kv_shear to be double counted. - RWH
1219  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1220  kv_add(i,k) = kv_add(i,k) + 0.5*(visc%Kv_slow(i,j,k) + visc%Kv_slow(i,j+1,k))
1221  endif ; enddo ; enddo
1222  !### I am pretty sure that this code is double counting viscosity at OBC points! - RWH
1223  if (do_obcs) then
1224  do i=is,ie ; if (do_i(i) .and. (obc%segnum_v(i,j) /= obc_none)) then
1225  if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
1226  do k=2,nz ; kv_add(i,k) = kv_add(i,k) + visc%Kv_slow(i,j,k) ; enddo
1227  elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
1228  do k=2,nz ; kv_add(i,k) = kv_add(i,k) + visc%Kv_slow(i,j+1,k) ; enddo
1229  endif
1230  endif ; enddo
1231  endif
1232  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1233  kv_tot(i,k) = kv_tot(i,k) + kv_add(i,k)
1234  endif ; enddo ; enddo
1235  endif
1236  endif
1237 
1238  do k=nz,2,-1 ; do i=is,ie ; if (do_i(i)) then
1239  ! botfn determines when a point is within the influence of the bottom
1240  ! boundary layer, going from 1 at the bottom to 0 in the interior.
1241  z2 = z_i(i,k)
1242  botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
1243 
1244  if (cs%bottomdraglaw) then
1245  kv_tot(i,k) = kv_tot(i,k) + (kv_bbl(i) - cs%Kv)*botfn
1246  r = 0.5*(hvel(i,k) + hvel(i,k-1))
1247  if (r > bbl_thick(i)) then
1248  h_shear = ((1.0 - botfn) * r + botfn*bbl_thick(i))
1249  else
1250  h_shear = r
1251  endif
1252  else
1253  kv_tot(i,k) = kv_tot(i,k) + (cs%Kvbbl-cs%Kv)*botfn
1254  h_shear = 0.5*(hvel(i,k) + hvel(i,k-1) + h_neglect)
1255  endif
1256 
1257  ! Calculate the coupling coefficients from the viscosities.
1258  a_cpl(i,k) = kv_tot(i,k) / (h_shear*gv%H_to_Z + i_amax*kv_tot(i,k))
1259  endif ; enddo ; enddo ! i & k loops
1260 
1261  if (do_shelf) then
1262  ! Set the coefficients to include the no-slip surface stress.
1263  do i=is,ie ; if (do_i(i)) then
1264  if (work_on_u) then
1265  kv_tbl(i) = visc%Kv_tbl_shelf_u(i,j)
1266  tbl_thick(i) = visc%tbl_thick_shelf_u(i,j) * gv%Z_to_H
1267  else
1268  kv_tbl(i) = visc%Kv_tbl_shelf_v(i,j)
1269  tbl_thick(i) = visc%tbl_thick_shelf_v(i,j) * gv%Z_to_H
1270  endif
1271  z_t(i) = 0.0
1272 
1273  ! If a_cpl(i,1) were not already 0, it would be added here.
1274  if (0.5*hvel(i,1) > tbl_thick(i)) then
1275  a_cpl(i,1) = kv_tbl(i) / (tbl_thick(i)*gv%H_to_Z + i_amax*kv_tbl(i))
1276  else
1277  a_cpl(i,1) = kv_tbl(i) / (0.5*hvel(i,1)*gv%H_to_Z + i_amax*kv_tbl(i))
1278  endif
1279  endif ; enddo
1280 
1281  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
1282  z_t(i) = z_t(i) + hvel(i,k-1) / tbl_thick(i)
1283  topfn = 1.0 / (1.0 + 0.09 * z_t(i)**6)
1284 
1285  r = 0.5*(hvel(i,k)+hvel(i,k-1))
1286  if (r > tbl_thick(i)) then
1287  h_shear = ((1.0 - topfn) * r + topfn*tbl_thick(i))
1288  else
1289  h_shear = r
1290  endif
1291 
1292  a_top = topfn * kv_tbl(i)
1293  a_cpl(i,k) = a_cpl(i,k) + a_top / (h_shear*gv%H_to_Z + i_amax*a_top)
1294  endif ; enddo ; enddo
1295  elseif (cs%dynamic_viscous_ML .or. (gv%nkml>0)) then
1296  max_nk = 0
1297  do i=is,ie ; if (do_i(i)) then
1298  if (gv%nkml>0) nk_visc(i) = real(gv%nkml+1)
1299  if (work_on_u) then
1300  u_star(i) = 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j))
1301  absf(i) = 0.5*(abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i,j)))
1302  if (cs%dynamic_viscous_ML) nk_visc(i) = visc%nkml_visc_u(i,j) + 1
1303  else
1304  u_star(i) = 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1))
1305  absf(i) = 0.5*(abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j)))
1306  if (cs%dynamic_viscous_ML) nk_visc(i) = visc%nkml_visc_v(i,j) + 1
1307  endif
1308  h_ml(i) = h_neglect ; z_t(i) = 0.0
1309  max_nk = max(max_nk,ceiling(nk_visc(i) - 1.0))
1310  endif ; enddo
1311 
1312  if (do_obcs) then ; if (work_on_u) then
1313  do i=is,ie ; if (do_i(i) .and. (obc%segnum_u(i,j) /= obc_none)) then
1314  if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) &
1315  u_star(i) = forces%ustar(i,j)
1316  if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) &
1317  u_star(i) = forces%ustar(i+1,j)
1318  endif ; enddo
1319  else
1320  do i=is,ie ; if (do_i(i) .and. (obc%segnum_v(i,j) /= obc_none)) then
1321  if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) &
1322  u_star(i) = forces%ustar(i,j)
1323  if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) &
1324  u_star(i) = forces%ustar(i,j+1)
1325  endif ; enddo
1326  endif ; endif
1327 
1328  do k=1,max_nk ; do i=is,ie ; if (do_i(i)) then
1329  if (k+1 <= nk_visc(i)) then ! This layer is all in the ML.
1330  h_ml(i) = h_ml(i) + hvel(i,k)
1331  elseif (k < nk_visc(i)) then ! Part of this layer is in the ML.
1332  h_ml(i) = h_ml(i) + (nk_visc(i) - k) * hvel(i,k)
1333  endif
1334  endif ; enddo ; enddo
1335 
1336  do k=2,max_nk ; do i=is,ie ; if (do_i(i)) then ; if (k < nk_visc(i)) then
1337  ! Set the viscosity at the interfaces.
1338  z_t(i) = z_t(i) + hvel(i,k-1)
1339  temp1 = (z_t(i)*h_ml(i) - z_t(i)*z_t(i))*gv%H_to_Z
1340  ! This viscosity is set to go to 0 at the mixed layer top and bottom (in a log-layer)
1341  ! and be further limited by rotation to give the natural Ekman length.
1342  visc_ml = u_star(i) * 0.41 * (temp1*u_star(i)) / (absf(i)*temp1 + h_ml(i)*u_star(i))
1343  a_ml = visc_ml / (0.25*(hvel(i,k)+hvel(i,k-1) + h_neglect) * gv%H_to_Z + 0.5*i_amax*visc_ml)
1344  ! Choose the largest estimate of a.
1345  if (a_ml > a_cpl(i,k)) a_cpl(i,k) = a_ml
1346  endif ; endif ; enddo ; enddo
1347  endif
1348 

◆ 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]timeCurrent model time
csVertical viscosity control structure
[in]activateSpecifiy whether to record the value of Time as the beginning of the ramp period

Definition at line 1791 of file MOM_vert_friction.F90.

1791  type(time_type), target, intent(in) :: Time !< Current model time
1792  type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
1793  logical, optional, intent(in) :: activate !< Specifiy whether to record the value of
1794  !! Time as the beginning of the ramp period
1795 
1796  ! Local variables
1797  real :: deltaTime, wghtA
1798  character(len=12) :: msg
1799 
1800  if (cs%truncRampTime==0.) return ! This indicates to ramping is turned off
1801 
1802  ! We use the optional argument to indicate this Time should be recorded as the
1803  ! beginning of the ramp-up period.
1804  if (present(activate)) then
1805  if (activate) then
1806  cs%rampStartTime = time ! Record the current time
1807  cs%CFLrampingIsActivated = .true.
1808  endif
1809  endif
1810  if (.not.cs%CFLrampingIsActivated) return
1811  deltatime = max( 0., time_type_to_real( time - cs%rampStartTime ) )
1812  if (deltatime >= cs%truncRampTime) then
1813  cs%CFL_trunc = cs%CFL_truncE
1814  cs%truncRampTime = 0. ! This turns off ramping after this call
1815  else
1816  wghta = min( 1., deltatime / cs%truncRampTime ) ! Linear profile in time
1817  !wghtA = wghtA*wghtA ! Convert linear profile to parabolic profile in time
1818  !wghtA = wghtA*wghtA*(3. - 2.*wghtA) ! Convert linear profile to cosine profile
1819  wghta = 1. - ( (1. - wghta)**2 ) ! Convert linear profiel to nverted parabolic profile
1820  cs%CFL_trunc = cs%CFL_truncS + wghta * ( cs%CFL_truncE - cs%CFL_truncS )
1821  endif
1822  write(msg(1:12),'(es12.3)') cs%CFL_trunc
1823  call mom_error(note, "MOM_vert_friction: updateCFLtruncationValue set CFL"// &
1824  " limit to "//trim(msg))

◆ 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]gOcean grid structure
[in]gvOcean vertical grid structure
[in]usA dimensional unit scaling type
[in,out]uZonal velocity [m s-1]
[in,out]vMeridional velocity [m s-1]
[in]hLayer thickness [H ~> m or kg m-2]
[in]forcesA structure with the driving mechanical forces
[in,out]viscViscosities and bottom drag
[in]dtTime increment [s]
obcOpen boundary condition structure
[in,out]adpAccelerations in the momentum equations for diagnostics
[in,out]cdpContinuity equation terms
csVertical viscosity control structure
[out]taux_botZonal bottom stress from ocean to rock [kg Z s-2 m-2 ~> Pa]
[out]tauy_botMeridional bottom stress from ocean to rock [kg Z s-2 m-2 ~> Pa]
wavesContainer for wave/Stokes information

Definition at line 147 of file MOM_vert_friction.F90.

147  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
148  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
149  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
150  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
151  intent(inout) :: u !< Zonal velocity [m s-1]
152  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
153  intent(inout) :: v !< Meridional velocity [m s-1]
154  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
155  intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
156  type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
157  type(vertvisc_type), intent(inout) :: visc !< Viscosities and bottom drag
158  real, intent(in) :: dt !< Time increment [s]
159  type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure
160  type(accel_diag_ptrs), intent(inout) :: ADp !< Accelerations in the momentum
161  !! equations for diagnostics
162  type(cont_diag_ptrs), intent(inout) :: CDp !< Continuity equation terms
163  type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
164  real, dimension(SZIB_(G),SZJ_(G)), &
165  optional, intent(out) :: taux_bot !< Zonal bottom stress from ocean to rock [kg Z s-2 m-2 ~> Pa]
166  real, dimension(SZI_(G),SZJB_(G)), &
167  optional, intent(out) :: tauy_bot !< Meridional bottom stress from ocean to rock [kg Z s-2 m-2 ~> Pa]
168  type(wave_parameters_CS), &
169  optional, pointer :: Waves !< Container for wave/Stokes information
170 
171  ! Fields from forces used in this subroutine:
172  ! taux: Zonal wind stress [Pa].
173  ! tauy: Meridional wind stress [Pa].
174 
175  ! Local variables
176 
177  real :: b1(SZIB_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
178  real :: c1(SZIB_(G),SZK_(G)) ! A variable used by the tridiagonal solver [nondim].
179  real :: d1(SZIB_(G)) ! d1=1-c1 is used by the tridiagonal solver [nondim].
180  real :: Ray(SZIB_(G),SZK_(G)) ! Ray is the Rayleigh-drag velocity [Z T-1 ~> m s-1].
181  real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2].
182 
183  real :: Hmix ! The mixed layer thickness over which stress
184  ! is applied with direct_stress [H ~> m or kg m-2].
185  real :: I_Hmix ! The inverse of Hmix [H-1 ~> m-1 or m2 kg-1].
186  real :: Idt ! The inverse of the time step [s-1].
187  real :: dt_Rho0 ! The time step divided by the mean density [s m3 kg-1].
188  real :: Rho0 ! A density used to convert drag laws into stress in Pa [kg m-3].
189  real :: dt_Z_to_H ! The time step times the conversion from Z to the
190  ! units of thickness - [T H Z-1 ~> s or s kg m-3].
191  real :: h_neglect ! A thickness that is so small it is usually lost
192  ! in roundoff and can be neglected [H ~> m or kg m-2].
193 
194  real :: stress ! The surface stress times the time step, divided
195  ! by the density [m2 s-1].
196  real :: zDS, hfr, h_a ! Temporary variables used with direct_stress.
197  real :: surface_stress(SZIB_(G))! The same as stress, unless the wind stress
198  ! stress is applied as a body force [m2 s-1].
199 
200  logical :: do_i(SZIB_(G))
201  logical :: DoStokesMixing
202 
203  integer :: i, j, k, is, ie, Isq, Ieq, Jsq, Jeq, nz, n
204  is = g%isc ; ie = g%iec
205  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
206 
207  if (.not.associated(cs)) call mom_error(fatal,"MOM_vert_friction(visc): "// &
208  "Module must be initialized before it is used.")
209 
210  if (cs%direct_stress) then
211  hmix = cs%Hmix_stress
212  i_hmix = 1.0 / hmix
213  endif
214  dt_rho0 = dt/gv%H_to_kg_m2
215  dt_z_to_h = us%s_to_T*dt*gv%Z_to_H
216  rho0 = gv%Rho0
217  h_neglect = gv%H_subroundoff
218  idt = 1.0 / dt
219 
220  !Check if Stokes mixing allowed if requested (present and associated)
221  dostokesmixing=.false.
222  if (cs%StokesMixing) then
223  if (present(waves)) dostokesmixing = associated(waves)
224  if (.not. dostokesmixing) &
225  call mom_error(fatal,"Stokes Mixing called without allocated"//&
226  "Waves Control Structure")
227  endif
228 
229  do k=1,nz ; do i=isq,ieq ; ray(i,k) = 0.0 ; enddo ; enddo
230 
231  ! Update the zonal velocity component using a modification of a standard
232  ! tridagonal solver.
233 
234  !$OMP parallel do default(shared) firstprivate(Ray) &
235  !$OMP private(do_i,surface_stress,zDS,stress,h_a,hfr, &
236  !$OMP b_denom_1,b1,d1,c1)
237  do j=g%jsc,g%jec
238  do i=isq,ieq ; do_i(i) = (g%mask2dCu(i,j) > 0) ; enddo
239 
240  ! When mixing down Eulerian current + Stokes drift add before calling solver
241  if (dostokesmixing) then ; do k=1,nz ; do i=isq,ieq
242  if (do_i(i)) u(i,j,k) = u(i,j,k) + waves%Us_x(i,j,k)
243  enddo ; enddo ; endif
244 
245  if (associated(adp%du_dt_visc)) then ; do k=1,nz ; do i=isq,ieq
246  adp%du_dt_visc(i,j,k) = u(i,j,k)
247  enddo ; enddo ; endif
248 
249 ! One option is to have the wind stress applied as a body force
250 ! over the topmost Hmix fluid. If DIRECT_STRESS is not defined,
251 ! the wind stress is applied as a stress boundary condition.
252  if (cs%direct_stress) then
253  do i=isq,ieq ; if (do_i(i)) then
254  surface_stress(i) = 0.0
255  zds = 0.0
256  stress = dt_rho0 * forces%taux(i,j)
257  do k=1,nz
258  h_a = 0.5 * (h(i,j,k) + h(i+1,j,k)) + h_neglect
259  hfr = 1.0 ; if ((zds+h_a) > hmix) hfr = (hmix - zds) / h_a
260  u(i,j,k) = u(i,j,k) + i_hmix * hfr * stress
261  zds = zds + h_a ; if (zds >= hmix) exit
262  enddo
263  endif ; enddo ! end of i loop
264  else ; do i=isq,ieq
265  surface_stress(i) = dt_rho0 * (g%mask2dCu(i,j)*forces%taux(i,j))
266  enddo ; endif ! direct_stress
267 
268  if (cs%Channel_drag) then ; do k=1,nz ; do i=isq,ieq
269  ray(i,k) = visc%Ray_u(i,j,k)
270  enddo ; enddo ; endif
271 
272  ! perform forward elimination on the tridiagonal system
273  !
274  ! denote the diagonal of the system as b_k, the subdiagonal as a_k
275  ! and the superdiagonal as c_k. The right-hand side terms are d_k.
276  !
277  ! ignoring the rayleigh drag contribution,
278  ! we have a_k = -dt_Z_to_H * a_u(k)
279  ! b_k = h_u(k) + dt_Z_to_H * (a_u(k) + a_u(k+1))
280  ! c_k = -dt_Z_to_H * a_u(k+1)
281  !
282  ! for forward elimination, we want to:
283  ! calculate c'_k = - c_k / (b_k + a_k c'_(k-1))
284  ! and d'_k = (d_k - a_k d'_(k-1)) / (b_k + a_k c'_(k-1))
285  ! where c'_1 = c_1/b_1 and d'_1 = d_1/b_1
286  ! (see Thomas' tridiagonal matrix algorithm)
287  !
288  ! b1 is the denominator term 1 / (b_k + a_k c'_(k-1))
289  ! b_denom_1 is (b_k + a_k + c_k) - a_k(1 - c'_(k-1))
290  ! = (b_k + c_k + c'_(k-1))
291  ! this is done so that d1 = b1 * b_denom_1 = 1 - c'_(k-1)
292  ! c1(k) is -c'_(k - 1)
293  ! and the right-hand-side is destructively updated to be d'_k
294  !
295  do i=isq,ieq ; if (do_i(i)) then
296  b_denom_1 = cs%h_u(i,j,1) + dt_z_to_h * (ray(i,1) + cs%a_u(i,j,1))
297  b1(i) = 1.0 / (b_denom_1 + dt_z_to_h*cs%a_u(i,j,2))
298  d1(i) = b_denom_1 * b1(i)
299  u(i,j,1) = b1(i) * (cs%h_u(i,j,1) * u(i,j,1) + surface_stress(i))
300  endif ; enddo
301  do k=2,nz ; do i=isq,ieq ; if (do_i(i)) then
302  c1(i,k) = dt_z_to_h * cs%a_u(i,j,k) * b1(i)
303  b_denom_1 = cs%h_u(i,j,k) + dt_z_to_h * (ray(i,k) + cs%a_u(i,j,k)*d1(i))
304  b1(i) = 1.0 / (b_denom_1 + dt_z_to_h * cs%a_u(i,j,k+1))
305  d1(i) = b_denom_1 * b1(i)
306  u(i,j,k) = (cs%h_u(i,j,k) * u(i,j,k) + &
307  dt_z_to_h * cs%a_u(i,j,k) * u(i,j,k-1)) * b1(i)
308  endif ; enddo ; enddo
309 
310  ! back substitute to solve for the new velocities
311  ! u_k = d'_k - c'_k x_(k+1)
312  do k=nz-1,1,-1 ; do i=isq,ieq ; if (do_i(i)) then
313  u(i,j,k) = u(i,j,k) + c1(i,k+1) * u(i,j,k+1)
314  endif ; enddo ; enddo ! i and k loops
315 
316  if (associated(adp%du_dt_visc)) then ; do k=1,nz ; do i=isq,ieq
317  adp%du_dt_visc(i,j,k) = (u(i,j,k) - adp%du_dt_visc(i,j,k))*idt
318  enddo ; enddo ; endif
319 
320  if (associated(visc%taux_shelf)) then ; do i=isq,ieq
321  visc%taux_shelf(i,j) = -rho0*us%s_to_T*cs%a1_shelf_u(i,j)*u(i,j,1) ! - u_shelf?
322  enddo ; endif
323 
324  if (PRESENT(taux_bot)) then
325  do i=isq,ieq
326  taux_bot(i,j) = rho0 * (u(i,j,nz)*us%s_to_T*cs%a_u(i,j,nz+1))
327  enddo
328  if (cs%Channel_drag) then ; do k=1,nz ; do i=isq,ieq
329  taux_bot(i,j) = taux_bot(i,j) + rho0 * (us%s_to_T*ray(i,k)*u(i,j,k))
330  enddo ; enddo ; endif
331  endif
332 
333  ! When mixing down Eulerian current + Stokes drift subtract after calling solver
334  if (dostokesmixing) then ; do k=1,nz ; do i=isq,ieq
335  if (do_i(i)) u(i,j,k) = u(i,j,k) - waves%Us_x(i,j,k)
336  enddo ; enddo ; endif
337 
338  enddo ! end u-component j loop
339 
340  ! Now work on the meridional velocity component.
341 
342  !$OMP parallel do default(shared) firstprivate(Ray) &
343  !$OMP private(do_i,surface_stress,zDS,stress,h_a,hfr, &
344  !$OMP b_denom_1,b1,d1,c1)
345  do j=jsq,jeq
346  do i=is,ie ; do_i(i) = (g%mask2dCv(i,j) > 0) ; enddo
347 
348  ! When mixing down Eulerian current + Stokes drift add before calling solver
349  if (dostokesmixing) then ; do k=1,nz ; do i=is,ie
350  if (do_i(i)) v(i,j,k) = v(i,j,k) + waves%Us_y(i,j,k)
351  enddo ; enddo ; endif
352 
353  if (associated(adp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie
354  adp%dv_dt_visc(i,j,k) = v(i,j,k)
355  enddo ; enddo ; endif
356 
357 ! One option is to have the wind stress applied as a body force
358 ! over the topmost Hmix fluid. If DIRECT_STRESS is not defined,
359 ! the wind stress is applied as a stress boundary condition.
360  if (cs%direct_stress) then
361  do i=is,ie ; if (do_i(i)) then
362  surface_stress(i) = 0.0
363  zds = 0.0
364  stress = dt_rho0 * forces%tauy(i,j)
365  do k=1,nz
366  h_a = 0.5 * (h(i,j,k) + h(i,j+1,k))
367  hfr = 1.0 ; if ((zds+h_a) > hmix) hfr = (hmix - zds) / h_a
368  v(i,j,k) = v(i,j,k) + i_hmix * hfr * stress
369  zds = zds + h_a ; if (zds >= hmix) exit
370  enddo
371  endif ; enddo ! end of i loop
372  else ; do i=is,ie
373  surface_stress(i) = dt_rho0 * (g%mask2dCv(i,j)*forces%tauy(i,j))
374  enddo ; endif ! direct_stress
375 
376  if (cs%Channel_drag) then ; do k=1,nz ; do i=is,ie
377  ray(i,k) = visc%Ray_v(i,j,k)
378  enddo ; enddo ; endif
379 
380  do i=is,ie ; if (do_i(i)) then
381  b_denom_1 = cs%h_v(i,j,1) + dt_z_to_h * (ray(i,1) + cs%a_v(i,j,1))
382  b1(i) = 1.0 / (b_denom_1 + dt_z_to_h*cs%a_v(i,j,2))
383  d1(i) = b_denom_1 * b1(i)
384  v(i,j,1) = b1(i) * (cs%h_v(i,j,1) * v(i,j,1) + surface_stress(i))
385  endif ; enddo
386  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
387  c1(i,k) = dt_z_to_h * cs%a_v(i,j,k) * b1(i)
388  b_denom_1 = cs%h_v(i,j,k) + dt_z_to_h * (ray(i,k) + cs%a_v(i,j,k)*d1(i))
389  b1(i) = 1.0 / (b_denom_1 + dt_z_to_h * cs%a_v(i,j,k+1))
390  d1(i) = b_denom_1 * b1(i)
391  v(i,j,k) = (cs%h_v(i,j,k) * v(i,j,k) + dt_z_to_h * cs%a_v(i,j,k) * v(i,j,k-1)) * b1(i)
392  endif ; enddo ; enddo
393  do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then
394  v(i,j,k) = v(i,j,k) + c1(i,k+1) * v(i,j,k+1)
395  endif ; enddo ; enddo ! i and k loops
396 
397  if (associated(adp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie
398  adp%dv_dt_visc(i,j,k) = (v(i,j,k) - adp%dv_dt_visc(i,j,k))*idt
399  enddo ; enddo ; endif
400 
401  if (associated(visc%tauy_shelf)) then ; do i=is,ie
402  visc%tauy_shelf(i,j) = -rho0*us%s_to_T*cs%a1_shelf_v(i,j)*v(i,j,1) ! - v_shelf?
403  enddo ; endif
404 
405  if (present(tauy_bot)) then
406  do i=is,ie
407  tauy_bot(i,j) = rho0 * (v(i,j,nz)*us%s_to_T*cs%a_v(i,j,nz+1))
408  enddo
409  if (cs%Channel_drag) then ; do k=1,nz ; do i=is,ie
410  tauy_bot(i,j) = tauy_bot(i,j) + rho0 * (us%s_to_T*ray(i,k)*v(i,j,k))
411  enddo ; enddo ; endif
412  endif
413 
414  ! When mixing down Eulerian current + Stokes drift subtract after calling solver
415  if (dostokesmixing) then ; do k=1,nz ; do i=is,ie
416  if (do_i(i)) v(i,j,k) = v(i,j,k) - waves%Us_y(i,j,k)
417  enddo ; enddo ; endif
418 
419  enddo ! end of v-component J loop
420 
421  call vertvisc_limit_vel(u, v, h, adp, cdp, forces, visc, dt, g, gv, us, cs)
422 
423  ! Here the velocities associated with open boundary conditions are applied.
424  if (associated(obc)) then
425  do n=1,obc%number_of_segments
426  if (obc%segment(n)%specified) then
427  if (obc%segment(n)%is_N_or_S) then
428  j = obc%segment(n)%HI%JsdB
429  do k=1,nz ; do i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied
430  v(i,j,k) = obc%segment(n)%normal_vel(i,j,k)
431  enddo ; enddo
432  elseif (obc%segment(n)%is_E_or_W) then
433  i = obc%segment(n)%HI%IsdB
434  do k=1,nz ; do j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed
435  u(i,j,k) = obc%segment(n)%normal_vel(i,j,k)
436  enddo ; enddo
437  endif
438  endif
439  enddo
440  endif
441 
442 ! Offer diagnostic fields for averaging.
443  if (cs%id_du_dt_visc > 0) &
444  call post_data(cs%id_du_dt_visc, adp%du_dt_visc, cs%diag)
445  if (cs%id_dv_dt_visc > 0) &
446  call post_data(cs%id_dv_dt_visc, adp%dv_dt_visc, cs%diag)
447  if (present(taux_bot) .and. (cs%id_taux_bot > 0)) &
448  call post_data(cs%id_taux_bot, taux_bot, cs%diag)
449  if (present(tauy_bot) .and. (cs%id_tauy_bot > 0)) &
450  call post_data(cs%id_tauy_bot, tauy_bot, cs%diag)
451 

◆ 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]gOcean grid structure
[in]gvOcean vertical grid structure
[in]usA dimensional unit scaling type
[in]uZonal velocity [m s-1]
[in]vMeridional velocity [m s-1]
[in]hLayer thickness [H ~> m or kg m-2]
[in]forcesA structure with the driving mechanical forces
[in]viscViscosities and bottom drag
[in]dtTime increment [s]
csVertical viscosity control structure
obcOpen boundary condition structure

Definition at line 567 of file MOM_vert_friction.F90.

567  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
568  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
569  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
570  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
571  intent(in) :: u !< Zonal velocity [m s-1]
572  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
573  intent(in) :: v !< Meridional velocity [m s-1]
574  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
575  intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
576  type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
577  type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag
578  real, intent(in) :: dt !< Time increment [s]
579  type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
580  type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure
581 
582  ! Field from forces used in this subroutine:
583  ! ustar: the friction velocity [m s-1], used here as the mixing
584  ! velocity in the mixed layer if NKML > 1 in a bulk mixed layer.
585 
586  ! Local variables
587 
588  real, dimension(SZIB_(G),SZK_(G)) :: &
589  h_harm, & ! Harmonic mean of the thicknesses around a velocity grid point,
590  ! given by 2*(h+ * h-)/(h+ + h-) [H ~> m or kg m-2].
591  h_arith, & ! The arithmetic mean thickness [H ~> m or kg m-2].
592  h_delta, & ! The lateral difference of thickness [H ~> m or kg m-2].
593  hvel, & ! hvel is the thickness used at a velocity grid point [H ~> m or kg m-2].
594  hvel_shelf ! The equivalent of hvel under shelves [H ~> m or kg m-2].
595  real, dimension(SZIB_(G),SZK_(G)+1) :: &
596  a_cpl, & ! The drag coefficients across interfaces [Z T-1 ~> m s-1]. a_cpl times
597  ! the velocity difference gives the stress across an interface.
598  a_shelf, & ! The drag coefficients across interfaces in water columns under
599  ! ice shelves [Z T-1 ~> m s-1].
600  z_i ! An estimate of each interface's height above the bottom,
601  ! normalized by the bottom boundary layer thickness, nondim.
602  real, dimension(SZIB_(G)) :: &
603  kv_bbl, & ! The bottom boundary layer viscosity [Z2 T-1 ~> m2 s-1].
604  bbl_thick, & ! The bottom boundary layer thickness [H ~> m or kg m-2].
605  I_Hbbl, & ! The inverse of the bottom boundary layer thickness [H-1 ~> m-1 or m2 kg-1].
606  I_Htbl, & ! The inverse of the top boundary layer thickness [H-1 ~> m-1 or m2 kg-1].
607  zcol1, & ! The height of the interfaces to the north and south of a
608  zcol2, & ! v-point [H ~> m or kg m-2].
609  Ztop_min, & ! The deeper of the two adjacent surface heights [H ~> m or kg m-2].
610  Dmin, & ! The shallower of the two adjacent bottom depths converted to
611  ! thickness units [H ~> m or kg m-2].
612  zh, & ! An estimate of the interface's distance from the bottom
613  ! based on harmonic mean thicknesses [H ~> m or kg m-2].
614  h_ml ! The mixed layer depth [H ~> m or kg m-2].
615  real, allocatable, dimension(:,:) :: hML_u ! Diagnostic of the mixed layer depth at u points [H ~> m or kg m-2].
616  real, allocatable, dimension(:,:) :: hML_v ! Diagnostic of the mixed layer depth at v points [H ~> m or kg m-2].
617  real, allocatable, dimension(:,:,:) :: Kv_u !< Total vertical viscosity at u-points [Z2 T-1 ~> m2 s-1].
618  real, allocatable, dimension(:,:,:) :: Kv_v !< Total vertical viscosity at v-points [Z2 T-1 ~> m2 s-1].
619  real :: zcol(SZI_(G)) ! The height of an interface at h-points [H ~> m or kg m-2].
620  real :: botfn ! A function which goes from 1 at the bottom to 0 much more
621  ! than Hbbl into the interior.
622  real :: topfn ! A function which goes from 1 at the top to 0 much more
623  ! than Htbl into the interior.
624  real :: z2 ! The distance from the bottom, normalized by Hbbl, nondim.
625  real :: z2_wt ! A nondimensional (0-1) weight used when calculating z2.
626  real :: z_clear ! The clearance of an interface above the surrounding topography [H ~> m or kg m-2].
627  real :: h_neglect ! A thickness that is so small it is usually lost
628  ! in roundoff and can be neglected [H ~> m or kg m-2].
629 
630  real :: I_valBL ! The inverse of a scaling factor determining when water is
631  ! still within the boundary layer, as determined by the sum
632  ! of the harmonic mean thicknesses.
633  logical, dimension(SZIB_(G)) :: do_i, do_i_shelf
634  logical :: do_any_shelf
635  integer, dimension(SZIB_(G)) :: &
636  zi_dir ! A trinary logical array indicating which thicknesses to use for
637  ! finding z_clear.
638  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
639  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
640  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
641 
642  if (.not.associated(cs)) call mom_error(fatal,"MOM_vert_friction(coef): "// &
643  "Module must be initialized before it is used.")
644 
645  h_neglect = gv%H_subroundoff
646  i_hbbl(:) = 1.0 / (cs%Hbbl + h_neglect)
647  i_valbl = 0.0 ; if (cs%harm_BL_val > 0.0) i_valbl = 1.0 / cs%harm_BL_val
648 
649  if (cs%id_Kv_u > 0) then
650  allocate(kv_u(g%IsdB:g%IedB,g%jsd:g%jed,g%ke)) ; kv_u(:,:,:) = 0.0
651  endif
652 
653  if (cs%id_Kv_v > 0) then
654  allocate(kv_v(g%isd:g%ied,g%JsdB:g%JedB,g%ke)) ; kv_v(:,:,:) = 0.0
655  endif
656 
657  if (cs%debug .or. (cs%id_hML_u > 0)) then
658  allocate(hml_u(g%IsdB:g%IedB,g%jsd:g%jed)) ; hml_u(:,:) = 0.0
659  endif
660  if (cs%debug .or. (cs%id_hML_v > 0)) then
661  allocate(hml_v(g%isd:g%ied,g%JsdB:g%JedB)) ; hml_v(:,:) = 0.0
662  endif
663 
664  if ((associated(visc%taux_shelf) .or. associated(forces%frac_shelf_u)) .and. &
665  .not.associated(cs%a1_shelf_u)) then
666  allocate(cs%a1_shelf_u(g%IsdB:g%IedB,g%jsd:g%jed)) ; cs%a1_shelf_u(:,:)=0.0
667  endif
668  if ((associated(visc%tauy_shelf) .or. associated(forces%frac_shelf_v)) .and. &
669  .not.associated(cs%a1_shelf_v)) then
670  allocate(cs%a1_shelf_v(g%isd:g%ied,g%JsdB:g%JedB)) ; cs%a1_shelf_v(:,:)=0.0
671  endif
672 
673  !$OMP parallel do default(private) shared(G,GV,US,CS,visc,Isq,Ieq,nz,u,h,forces,hML_u, &
674  !$OMP OBC,h_neglect,dt,I_valBL,Kv_u) &
675  !$OMP firstprivate(i_hbbl)
676  do j=g%Jsc,g%Jec
677  do i=isq,ieq ; do_i(i) = (g%mask2dCu(i,j) > 0) ; enddo
678 
679  if (cs%bottomdraglaw) then ; do i=isq,ieq
680  kv_bbl(i) = visc%Kv_bbl_u(i,j)
681  bbl_thick(i) = visc%bbl_thick_u(i,j) * gv%Z_to_H
682  if (do_i(i)) i_hbbl(i) = 1.0 / (bbl_thick(i) + h_neglect)
683  enddo ; endif
684 
685  do k=1,nz ; do i=isq,ieq ; if (do_i(i)) then
686  h_harm(i,k) = 2.0*h(i,j,k)*h(i+1,j,k) / (h(i,j,k)+h(i+1,j,k)+h_neglect)
687  h_arith(i,k) = 0.5*(h(i+1,j,k)+h(i,j,k))
688  h_delta(i,k) = h(i+1,j,k) - h(i,j,k)
689  endif ; enddo ; enddo
690  do i=isq,ieq
691  dmin(i) = min(g%bathyT(i,j), g%bathyT(i+1,j)) * gv%Z_to_H
692  zi_dir(i) = 0
693  enddo
694 
695  ! Project thickness outward across OBCs using a zero-gradient condition.
696  if (associated(obc)) then ; if (obc%number_of_segments > 0) then
697  do i=isq,ieq ; if (do_i(i) .and. (obc%segnum_u(i,j) /= obc_none)) then
698  if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
699  do k=1,nz ; h_harm(i,k) = h(i,j,k) ; h_arith(i,k) = h(i,j,k) ; h_delta(i,k) = 0. ; enddo
700  dmin(i) = g%bathyT(i,j) * gv%Z_to_H
701  zi_dir(i) = -1
702  elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
703  do k=1,nz ; h_harm(i,k) = h(i+1,j,k) ; h_arith(i,k) = h(i+1,j,k) ; h_delta(i,k) = 0. ; enddo
704  dmin(i) = g%bathyT(i+1,j) * gv%Z_to_H
705  zi_dir(i) = 1
706  endif
707  endif ; enddo
708  endif ; endif
709 
710 ! The following block calculates the thicknesses at velocity
711 ! grid points for the vertical viscosity (hvel). Near the
712 ! bottom an upwind biased thickness is used to control the effect
713 ! of spurious Montgomery potential gradients at the bottom where
714 ! nearly massless layers layers ride over the topography.
715  if (cs%harmonic_visc) then
716  do i=isq,ieq ; z_i(i,nz+1) = 0.0 ; enddo
717  do k=nz,1,-1 ; do i=isq,ieq ; if (do_i(i)) then
718  hvel(i,k) = h_harm(i,k)
719  if (u(i,j,k) * h_delta(i,k) < 0) then
720  z2 = z_i(i,k+1) ; botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
721  hvel(i,k) = (1.0-botfn)*h_harm(i,k) + botfn*h_arith(i,k)
722  endif
723  z_i(i,k) = z_i(i,k+1) + h_harm(i,k)*i_hbbl(i)
724  endif ; enddo ; enddo ! i & k loops
725  else ! Not harmonic_visc
726  do i=isq,ieq ; zh(i) = 0.0 ; z_i(i,nz+1) = 0.0 ; enddo
727  do i=isq,ieq+1 ; zcol(i) = -g%bathyT(i,j) * gv%Z_to_H ; enddo
728  do k=nz,1,-1
729  do i=isq,ieq+1 ; zcol(i) = zcol(i) + h(i,j,k) ; enddo
730  do i=isq,ieq ; if (do_i(i)) then
731  zh(i) = zh(i) + h_harm(i,k)
732 
733  z_clear = max(zcol(i),zcol(i+1)) + dmin(i)
734  if (zi_dir(i) < 0) z_clear = zcol(i) + dmin(i)
735  if (zi_dir(i) > 0) z_clear = zcol(i+1) + dmin(i)
736 
737  z_i(i,k) = max(zh(i), z_clear) * i_hbbl(i)
738 
739  hvel(i,k) = h_arith(i,k)
740  if (u(i,j,k) * h_delta(i,k) > 0) then
741  if (zh(i) * i_hbbl(i) < cs%harm_BL_val) then
742  hvel(i,k) = h_harm(i,k)
743  else
744  z2_wt = 1.0 ; if (zh(i) * i_hbbl(i) < 2.0*cs%harm_BL_val) &
745  z2_wt = max(0.0, min(1.0, zh(i) * i_hbbl(i) * i_valbl - 1.0))
746  z2 = z2_wt * (max(zh(i), z_clear) * i_hbbl(i))
747  botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
748  hvel(i,k) = (1.0-botfn)*h_arith(i,k) + botfn*h_harm(i,k)
749  endif
750  endif
751 
752  endif ; enddo ! i loop
753  enddo ! k loop
754  endif
755 
756  call find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_ml, &
757  dt, j, g, gv, us, cs, visc, forces, work_on_u=.true., obc=obc)
758  if (allocated(hml_u)) then
759  do i=isq,ieq ; if (do_i(i)) then ; hml_u(i,j) = h_ml(i) ; endif ; enddo
760  endif
761 
762  do_any_shelf = .false.
763  if (associated(forces%frac_shelf_u)) then
764  do i=isq,ieq
765  cs%a1_shelf_u(i,j) = 0.0
766  do_i_shelf(i) = (do_i(i) .and. forces%frac_shelf_u(i,j) > 0.0)
767  if (do_i_shelf(i)) do_any_shelf = .true.
768  enddo
769  if (do_any_shelf) then
770  if (cs%harmonic_visc) then
771  call find_coupling_coef(a_shelf, hvel, do_i_shelf, h_harm, bbl_thick, &
772  kv_bbl, z_i, h_ml, dt, j, g, gv, us, cs, &
773  visc, forces, work_on_u=.true., obc=obc, shelf=.true.)
774  else ! Find upwind-biased thickness near the surface.
775  ! Perhaps this needs to be done more carefully, via find_eta.
776  do i=isq,ieq ; if (do_i_shelf(i)) then
777  zh(i) = 0.0 ; ztop_min(i) = min(zcol(i), zcol(i+1))
778  i_htbl(i) = 1.0 / (visc%tbl_thick_shelf_u(i,j)*gv%Z_to_H + h_neglect)
779  endif ; enddo
780  do k=1,nz
781  do i=isq,ieq+1 ; zcol(i) = zcol(i) - h(i,j,k) ; enddo
782  do i=isq,ieq ; if (do_i_shelf(i)) then
783  zh(i) = zh(i) + h_harm(i,k)
784 
785  hvel_shelf(i,k) = hvel(i,k)
786  if (u(i,j,k) * h_delta(i,k) > 0) then
787  if (zh(i) * i_htbl(i) < cs%harm_BL_val) then
788  hvel_shelf(i,k) = min(hvel(i,k), h_harm(i,k))
789  else
790  z2_wt = 1.0 ; if (zh(i) * i_htbl(i) < 2.0*cs%harm_BL_val) &
791  z2_wt = max(0.0, min(1.0, zh(i) * i_htbl(i) * i_valbl - 1.0))
792  z2 = z2_wt * (max(zh(i), ztop_min(i) - min(zcol(i),zcol(i+1))) * i_htbl(i))
793  topfn = 1.0 / (1.0 + 0.09*z2**6)
794  hvel_shelf(i,k) = min(hvel(i,k), (1.0-topfn)*h_arith(i,k) + topfn*h_harm(i,k))
795  endif
796  endif
797  endif ; enddo
798  enddo
799  call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, h_harm, &
800  bbl_thick, kv_bbl, z_i, h_ml, dt, j, g, gv, us, cs, &
801  visc, forces, work_on_u=.true., obc=obc, shelf=.true.)
802  endif
803  do i=isq,ieq ; if (do_i_shelf(i)) cs%a1_shelf_u(i,j) = a_shelf(i,1) ; enddo
804  endif
805  endif
806 
807  if (do_any_shelf) then
808  do k=1,nz+1 ; do i=isq,ieq ; if (do_i_shelf(i)) then
809  cs%a_u(i,j,k) = forces%frac_shelf_u(i,j) * a_shelf(i,k) + &
810  (1.0-forces%frac_shelf_u(i,j)) * a_cpl(i,k)
811 ! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH
812 ! CS%a_u(I,j,K) = forces%frac_shelf_u(I,j) * max(a_shelf(I,K), a_cpl(I,K)) + &
813 ! (1.0-forces%frac_shelf_u(I,j)) * a_cpl(I,K)
814  elseif (do_i(i)) then
815  cs%a_u(i,j,k) = a_cpl(i,k)
816  endif ; enddo ; enddo
817  do k=1,nz ; do i=isq,ieq ; if (do_i_shelf(i)) then
818  ! Should we instead take the inverse of the average of the inverses?
819  cs%h_u(i,j,k) = forces%frac_shelf_u(i,j) * hvel_shelf(i,k) + &
820  (1.0-forces%frac_shelf_u(i,j)) * hvel(i,k)
821  elseif (do_i(i)) then
822  cs%h_u(i,j,k) = hvel(i,k)
823  endif ; enddo ; enddo
824  else
825  do k=1,nz+1 ; do i=isq,ieq ; if (do_i(i)) cs%a_u(i,j,k) = a_cpl(i,k) ; enddo ; enddo
826  do k=1,nz ; do i=isq,ieq ; if (do_i(i)) cs%h_u(i,j,k) = hvel(i,k) ; enddo ; enddo
827  endif
828 
829  ! Diagnose total Kv at u-points
830  if (cs%id_Kv_u > 0) then
831  do k=1,nz ; do i=isq,ieq
832  if (do_i(i)) kv_u(i,j,k) = 0.5 * gv%H_to_Z*(cs%a_u(i,j,k)+cs%a_u(i,j,k+1)) * cs%h_u(i,j,k)
833  enddo ; enddo
834  endif
835 
836  enddo
837 
838 
839  ! Now work on v-points.
840  !$OMP parallel do default(private) shared(G,GV,CS,US,visc,is,ie,Jsq,Jeq,nz,v,h,forces,hML_v, &
841  !$OMP OBC,h_neglect,dt,I_valBL,Kv_v) &
842  !$OMP firstprivate(i_hbbl)
843  do j=jsq,jeq
844  do i=is,ie ; do_i(i) = (g%mask2dCv(i,j) > 0) ; enddo
845 
846  if (cs%bottomdraglaw) then ; do i=is,ie
847  kv_bbl(i) = visc%Kv_bbl_v(i,j)
848  bbl_thick(i) = visc%bbl_thick_v(i,j) * gv%Z_to_H
849  if (do_i(i)) i_hbbl(i) = 1.0 / bbl_thick(i)
850  enddo ; endif
851 
852  do k=1,nz ; do i=is,ie ; if (do_i(i)) then
853  h_harm(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / (h(i,j,k)+h(i,j+1,k)+h_neglect)
854  h_arith(i,k) = 0.5*(h(i,j+1,k)+h(i,j,k))
855  h_delta(i,k) = h(i,j+1,k) - h(i,j,k)
856  endif ; enddo ; enddo
857  do i=is,ie
858  dmin(i) = min(g%bathyT(i,j), g%bathyT(i,j+1)) * gv%Z_to_H
859  zi_dir(i) = 0
860  enddo
861 
862  ! Project thickness outward across OBCs using a zero-gradient condition.
863  if (associated(obc)) then ; if (obc%number_of_segments > 0) then
864  do i=is,ie ; if (do_i(i) .and. (obc%segnum_v(i,j) /= obc_none)) then
865  if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
866  do k=1,nz ; h_harm(i,k) = h(i,j,k) ; h_arith(i,k) = h(i,j,k) ; h_delta(i,k) = 0. ; enddo
867  dmin(i) = g%bathyT(i,j) * gv%Z_to_H
868  zi_dir(i) = -1
869  elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
870  do k=1,nz ; h_harm(i,k) = h(i,j+1,k) ; h_arith(i,k) = h(i,j+1,k) ; h_delta(i,k) = 0. ; enddo
871  dmin(i) = g%bathyT(i,j+1) * gv%Z_to_H
872  zi_dir(i) = 1
873  endif
874  endif ; enddo
875  endif ; endif
876 
877 ! The following block calculates the thicknesses at velocity
878 ! grid points for the vertical viscosity (hvel). Near the
879 ! bottom an upwind biased thickness is used to control the effect
880 ! of spurious Montgomery potential gradients at the bottom where
881 ! nearly massless layers layers ride over the topography.
882  if (cs%harmonic_visc) then
883  do i=is,ie ; z_i(i,nz+1) = 0.0 ; enddo
884 
885  do k=nz,1,-1 ; do i=is,ie ; if (do_i(i)) then
886  hvel(i,k) = h_harm(i,k)
887  if (v(i,j,k) * h_delta(i,k) < 0) then
888  z2 = z_i(i,k+1) ; botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
889  hvel(i,k) = (1.0-botfn)*h_harm(i,k) + botfn*h_arith(i,k)
890  endif
891  z_i(i,k) = z_i(i,k+1) + h_harm(i,k)*i_hbbl(i)
892  endif ; enddo ; enddo ! i & k loops
893  else ! Not harmonic_visc
894  do i=is,ie
895  zh(i) = 0.0 ; z_i(i,nz+1) = 0.0
896  zcol1(i) = -g%bathyT(i,j) * gv%Z_to_H
897  zcol2(i) = -g%bathyT(i,j+1) * gv%Z_to_H
898  enddo
899  do k=nz,1,-1 ; do i=is,ie ; if (do_i(i)) then
900  zh(i) = zh(i) + h_harm(i,k)
901  zcol1(i) = zcol1(i) + h(i,j,k) ; zcol2(i) = zcol2(i) + h(i,j+1,k)
902 
903  z_clear = max(zcol1(i),zcol2(i)) + dmin(i)
904  if (zi_dir(i) < 0) z_clear = zcol1(i) + dmin(i)
905  if (zi_dir(i) > 0) z_clear = zcol2(i) + dmin(i)
906 
907  z_i(i,k) = max(zh(i), z_clear) * i_hbbl(i)
908 
909  hvel(i,k) = h_arith(i,k)
910  if (v(i,j,k) * h_delta(i,k) > 0) then
911  if (zh(i) * i_hbbl(i) < cs%harm_BL_val) then
912  hvel(i,k) = h_harm(i,k)
913  else
914  z2_wt = 1.0 ; if (zh(i) * i_hbbl(i) < 2.0*cs%harm_BL_val) &
915  z2_wt = max(0.0, min(1.0, zh(i) * i_hbbl(i) * i_valbl - 1.0))
916  z2 = z2_wt * (max(zh(i), max(zcol1(i),zcol2(i)) + dmin(i)) * i_hbbl(i))
917  botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2)
918  hvel(i,k) = (1.0-botfn)*h_arith(i,k) + botfn*h_harm(i,k)
919  endif
920  endif
921 
922  endif ; enddo ; enddo ! i & k loops
923  endif
924 
925  call find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_ml, &
926  dt, j, g, gv, us, cs, visc, forces, work_on_u=.false., obc=obc)
927  if ( allocated(hml_v)) then
928  do i=is,ie ; if (do_i(i)) then ; hml_v(i,j) = h_ml(i) ; endif ; enddo
929  endif
930  do_any_shelf = .false.
931  if (associated(forces%frac_shelf_v)) then
932  do i=is,ie
933  cs%a1_shelf_v(i,j) = 0.0
934  do_i_shelf(i) = (do_i(i) .and. forces%frac_shelf_v(i,j) > 0.0)
935  if (do_i_shelf(i)) do_any_shelf = .true.
936  enddo
937  if (do_any_shelf) then
938  if (cs%harmonic_visc) then
939  call find_coupling_coef(a_shelf, hvel, do_i_shelf, h_harm, bbl_thick, &
940  kv_bbl, z_i, h_ml, dt, j, g, gv, us, cs, visc, &
941  forces, work_on_u=.false., obc=obc, shelf=.true.)
942  else ! Find upwind-biased thickness near the surface.
943  ! Perhaps this needs to be done more carefully, via find_eta.
944  do i=is,ie ; if (do_i_shelf(i)) then
945  zh(i) = 0.0 ; ztop_min(i) = min(zcol1(i), zcol2(i))
946  i_htbl(i) = 1.0 / (visc%tbl_thick_shelf_v(i,j)*gv%Z_to_H + h_neglect)
947  endif ; enddo
948  do k=1,nz
949  do i=is,ie ; if (do_i_shelf(i)) then
950  zcol1(i) = zcol1(i) - h(i,j,k) ; zcol2(i) = zcol2(i) - h(i,j+1,k)
951  zh(i) = zh(i) + h_harm(i,k)
952 
953  hvel_shelf(i,k) = hvel(i,k)
954  if (v(i,j,k) * h_delta(i,k) > 0) then
955  if (zh(i) * i_htbl(i) < cs%harm_BL_val) then
956  hvel_shelf(i,k) = min(hvel(i,k), h_harm(i,k))
957  else
958  z2_wt = 1.0 ; if (zh(i) * i_htbl(i) < 2.0*cs%harm_BL_val) &
959  z2_wt = max(0.0, min(1.0, zh(i) * i_htbl(i) * i_valbl - 1.0))
960  z2 = z2_wt * (max(zh(i), ztop_min(i) - min(zcol1(i),zcol2(i))) * i_htbl(i))
961  topfn = 1.0 / (1.0 + 0.09*z2**6)
962  hvel_shelf(i,k) = min(hvel(i,k), (1.0-topfn)*h_arith(i,k) + topfn*h_harm(i,k))
963  endif
964  endif
965  endif ; enddo
966  enddo
967  call find_coupling_coef(a_shelf, hvel_shelf, do_i_shelf, h_harm, &
968  bbl_thick, kv_bbl, z_i, h_ml, dt, j, g, gv, us, cs, &
969  visc, forces, work_on_u=.false., obc=obc, shelf=.true.)
970  endif
971  do i=is,ie ; if (do_i_shelf(i)) cs%a1_shelf_v(i,j) = a_shelf(i,1) ; enddo
972  endif
973  endif
974 
975  if (do_any_shelf) then
976  do k=1,nz+1 ; do i=is,ie ; if (do_i_shelf(i)) then
977  cs%a_v(i,j,k) = forces%frac_shelf_v(i,j) * a_shelf(i,k) + &
978  (1.0-forces%frac_shelf_v(i,j)) * a_cpl(i,k)
979 ! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH
980 ! CS%a_v(i,J,K) = forces%frac_shelf_v(i,J) * max(a_shelf(i,K), a_cpl(i,K)) + &
981 ! (1.0-forces%frac_shelf_v(i,J)) * a_cpl(i,K)
982  elseif (do_i(i)) then
983  cs%a_v(i,j,k) = a_cpl(i,k)
984  endif ; enddo ; enddo
985  do k=1,nz ; do i=is,ie ; if (do_i_shelf(i)) then
986  ! Should we instead take the inverse of the average of the inverses?
987  cs%h_v(i,j,k) = forces%frac_shelf_v(i,j) * hvel_shelf(i,k) + &
988  (1.0-forces%frac_shelf_v(i,j)) * hvel(i,k)
989  elseif (do_i(i)) then
990  cs%h_v(i,j,k) = hvel(i,k)
991  endif ; enddo ; enddo
992  else
993  do k=1,nz+1 ; do i=is,ie ; if (do_i(i)) cs%a_v(i,j,k) = a_cpl(i,k) ; enddo ; enddo
994  do k=1,nz ; do i=is,ie ; if (do_i(i)) cs%h_v(i,j,k) = hvel(i,k) ; enddo ; enddo
995  endif
996 
997  ! Diagnose total Kv at v-points
998  if (cs%id_Kv_v > 0) then
999  do k=1,nz ; do i=is,ie
1000  if (do_i(i)) kv_v(i,j,k) = 0.5 * gv%H_to_Z*(cs%a_v(i,j,k)+cs%a_v(i,j,k+1)) * cs%h_v(i,j,k)
1001  enddo ; enddo
1002  endif
1003 
1004  enddo ! end of v-point j loop
1005 
1006  if (cs%debug) then
1007  call uvchksum("vertvisc_coef h_[uv]", cs%h_u, &
1008  cs%h_v, g%HI,haloshift=0, scale=gv%H_to_m*us%s_to_T)
1009  call uvchksum("vertvisc_coef a_[uv]", cs%a_u, &
1010  cs%a_v, g%HI, haloshift=0, scale=us%Z_to_m*us%s_to_T)
1011  if (allocated(hml_u) .and. allocated(hml_v)) &
1012  call uvchksum("vertvisc_coef hML_[uv]", hml_u, hml_v, &
1013  g%HI, haloshift=0, scale=gv%H_to_m)
1014  endif
1015 
1016 ! Offer diagnostic fields for averaging.
1017  if (cs%id_Kv_slow > 0) call post_data(cs%id_Kv_slow, visc%Kv_slow, cs%diag)
1018  if (cs%id_Kv_u > 0) call post_data(cs%id_Kv_u, kv_u, cs%diag)
1019  if (cs%id_Kv_v > 0) call post_data(cs%id_Kv_v, kv_v, cs%diag)
1020  if (cs%id_au_vv > 0) call post_data(cs%id_au_vv, cs%a_u, cs%diag)
1021  if (cs%id_av_vv > 0) call post_data(cs%id_av_vv, cs%a_v, cs%diag)
1022  if (cs%id_h_u > 0) call post_data(cs%id_h_u, cs%h_u, cs%diag)
1023  if (cs%id_h_v > 0) call post_data(cs%id_h_v, cs%h_v, cs%diag)
1024  if (cs%id_hML_u > 0) call post_data(cs%id_hML_u, hml_u, cs%diag)
1025  if (cs%id_hML_v > 0) call post_data(cs%id_hML_v, hml_v, cs%diag)
1026 
1027  if (allocated(hml_u)) deallocate(hml_u)
1028  if (allocated(hml_v)) deallocate(hml_v)
1029 

◆ vertvisc_end()

subroutine, public mom_vert_friction::vertvisc_end ( type(vertvisc_cs), pointer  CS)

Clean up and deallocate the vertical friction module.

Parameters
csVertical viscosity control structure that will be deallocated in this subroutine.

Definition at line 1829 of file MOM_vert_friction.F90.

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

◆ 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]misThe "MOM Internal State", a set of pointers
[in]timeCurrent model time
[in]gOcean grid structure
[in]gvOcean vertical grid structure
[in]usA dimensional unit scaling type
[in]param_fileFile to parse for parameters
[in,out]diagDiagnostic control structure
[in,out]adpAcceleration diagnostic pointers
[in]dirsRelevant directory paths
[in,out]ntruncNumber of velocity truncations
csVertical viscosity control structure

Definition at line 1566 of file MOM_vert_friction.F90.

1566  type(ocean_internal_state), &
1567  target, intent(in) :: MIS !< The "MOM Internal State", a set of pointers
1568  !! to the fields and accelerations that make
1569  !! up the ocean's physical state
1570  type(time_type), target, intent(in) :: Time !< Current model time
1571  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1572  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1573  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1574  type(param_file_type), intent(in) :: param_file !< File to parse for parameters
1575  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostic control structure
1576  type(accel_diag_ptrs), intent(inout) :: ADp !< Acceleration diagnostic pointers
1577  type(directories), intent(in) :: dirs !< Relevant directory paths
1578  integer, target, intent(inout) :: ntrunc !< Number of velocity truncations
1579  type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
1580 
1581  ! Local variables
1582 
1583  real :: hmix_str_dflt
1584  real :: Kv_dflt ! A default viscosity [m2 s-1].
1585  real :: Hmix_m ! A boundary layer thickness [m].
1586  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
1587 ! This include declares and sets the variable "version".
1588 #include "version_variable.h"
1589  character(len=40) :: mdl = "MOM_vert_friction" ! This module's name.
1590  character(len=40) :: thickness_units
1591 
1592  if (associated(cs)) then
1593  call mom_error(warning, "vertvisc_init called with an associated "// &
1594  "control structure.")
1595  return
1596  endif
1597  allocate(cs)
1598 
1599  if (gv%Boussinesq) then; thickness_units = "m"
1600  else; thickness_units = "kg m-2"; endif
1601 
1602  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
1603  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1604 
1605  cs%diag => diag ; cs%ntrunc => ntrunc ; ntrunc = 0
1606 
1607 ! Default, read and log parameters
1608  call log_version(param_file, mdl, version, "")
1609  call get_param(param_file, mdl, "BOTTOMDRAGLAW", cs%bottomdraglaw, &
1610  "If true, the bottom stress is calculated with a drag "//&
1611  "law of the form c_drag*|u|*u. The velocity magnitude "//&
1612  "may be an assumed value or it may be based on the "//&
1613  "actual velocity in the bottommost HBBL, depending on "//&
1614  "LINEAR_DRAG.", default=.true.)
1615  call get_param(param_file, mdl, "CHANNEL_DRAG", cs%Channel_drag, &
1616  "If true, the bottom drag is exerted directly on each "//&
1617  "layer proportional to the fraction of the bottom it "//&
1618  "overlies.", default=.false.)
1619  call get_param(param_file, mdl, "DIRECT_STRESS", cs%direct_stress, &
1620  "If true, the wind stress is distributed over the "//&
1621  "topmost HMIX_STRESS of fluid (like in HYCOM), and KVML "//&
1622  "may be set to a very small value.", default=.false.)
1623  call get_param(param_file, mdl, "DYNAMIC_VISCOUS_ML", cs%dynamic_viscous_ML, &
1624  "If true, use a bulk Richardson number criterion to "//&
1625  "determine the mixed layer thickness for viscosity.", &
1626  default=.false.)
1627  call get_param(param_file, mdl, "U_TRUNC_FILE", cs%u_trunc_file, &
1628  "The absolute path to a file into which the accelerations "//&
1629  "leading to zonal velocity truncations are written. "//&
1630  "Undefine this for efficiency if this diagnostic is not "//&
1631  "needed.", default=" ", debuggingparam=.true.)
1632  call get_param(param_file, mdl, "V_TRUNC_FILE", cs%v_trunc_file, &
1633  "The absolute path to a file into which the accelerations "//&
1634  "leading to meridional velocity truncations are written. "//&
1635  "Undefine this for efficiency if this diagnostic is not "//&
1636  "needed.", default=" ", debuggingparam=.true.)
1637  call get_param(param_file, mdl, "HARMONIC_VISC", cs%harmonic_visc, &
1638  "If true, use the harmonic mean thicknesses for "//&
1639  "calculating the vertical viscosity.", default=.false.)
1640  call get_param(param_file, mdl, "HARMONIC_BL_SCALE", cs%harm_BL_val, &
1641  "A scale to determine when water is in the boundary "//&
1642  "layers based solely on harmonic mean thicknesses for "//&
1643  "the purpose of determining the extent to which the "//&
1644  "thicknesses used in the viscosities are upwinded.", &
1645  default=0.0, units="nondim")
1646  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
1647 
1648  if (gv%nkml < 1) &
1649  call get_param(param_file, mdl, "HMIX_FIXED", cs%Hmix, &
1650  "The prescribed depth over which the near-surface "//&
1651  "viscosity and diffusivity are elevated when the bulk "//&
1652  "mixed layer is not used.", units="m", scale=gv%m_to_H, &
1653  unscaled=hmix_m, fail_if_missing=.true.)
1654  if (cs%direct_stress) then
1655  if (gv%nkml < 1) then
1656  call get_param(param_file, mdl, "HMIX_STRESS", cs%Hmix_stress, &
1657  "The depth over which the wind stress is applied if "//&
1658  "DIRECT_STRESS is true.", units="m", default=hmix_m, scale=gv%m_to_H)
1659  else
1660  call get_param(param_file, mdl, "HMIX_STRESS", cs%Hmix_stress, &
1661  "The depth over which the wind stress is applied if "//&
1662  "DIRECT_STRESS is true.", units="m", fail_if_missing=.true., scale=gv%m_to_H)
1663  endif
1664  if (cs%Hmix_stress <= 0.0) call mom_error(fatal, "vertvisc_init: " // &
1665  "HMIX_STRESS must be set to a positive value if DIRECT_STRESS is true.")
1666  endif
1667  call get_param(param_file, mdl, "KV", cs%Kv, &
1668  "The background kinematic viscosity in the interior. "//&
1669  "The molecular value, ~1e-6 m2 s-1, may be used.", &
1670  units="m2 s-1", fail_if_missing=.true., scale=us%m2_s_to_Z2_T, unscaled=kv_dflt)
1671 
1672  if (gv%nkml < 1) call get_param(param_file, mdl, "KVML", cs%Kvml, &
1673  "The kinematic viscosity in the mixed layer. A typical "//&
1674  "value is ~1e-2 m2 s-1. KVML is not used if "//&
1675  "BULKMIXEDLAYER is true. The default is set by KV.", &
1676  units="m2 s-1", default=kv_dflt, scale=us%m2_s_to_Z2_T)
1677  if (.not.cs%bottomdraglaw) call get_param(param_file, mdl, "KVBBL", cs%Kvbbl, &
1678  "The kinematic viscosity in the benthic boundary layer. "//&
1679  "A typical value is ~1e-2 m2 s-1. KVBBL is not used if "//&
1680  "BOTTOMDRAGLAW is true. The default is set by KV.", &
1681  units="m2 s-1", default=kv_dflt, scale=us%m2_s_to_Z2_T)
1682  call get_param(param_file, mdl, "HBBL", cs%Hbbl, &
1683  "The thickness of a bottom boundary layer with a "//&
1684  "viscosity of KVBBL if BOTTOMDRAGLAW is not defined, or "//&
1685  "the thickness over which near-bottom velocities are "//&
1686  "averaged for the drag law if BOTTOMDRAGLAW is defined "//&
1687  "but LINEAR_DRAG is not.", units="m", fail_if_missing=.true., scale=gv%m_to_H)
1688  call get_param(param_file, mdl, "MAXVEL", cs%maxvel, &
1689  "The maximum velocity allowed before the velocity "//&
1690  "components are truncated.", units="m s-1", default=3.0e8)
1691  call get_param(param_file, mdl, "CFL_BASED_TRUNCATIONS", cs%CFL_based_trunc, &
1692  "If true, base truncations on the CFL number, and not an "//&
1693  "absolute speed.", default=.true.)
1694  call get_param(param_file, mdl, "CFL_TRUNCATE", cs%CFL_trunc, &
1695  "The value of the CFL number that will cause velocity "//&
1696  "components to be truncated; instability can occur past 0.5.", &
1697  units="nondim", default=0.5)
1698  call get_param(param_file, mdl, "CFL_REPORT", cs%CFL_report, &
1699  "The value of the CFL number that causes accelerations "//&
1700  "to be reported; the default is CFL_TRUNCATE.", &
1701  units="nondim", default=cs%CFL_trunc)
1702  call get_param(param_file, mdl, "CFL_TRUNCATE_RAMP_TIME", cs%truncRampTime, &
1703  "The time over which the CFL truncation value is ramped "//&
1704  "up at the beginning of the run.", &
1705  units="s", default=0.)
1706  cs%CFL_truncE = cs%CFL_trunc
1707  call get_param(param_file, mdl, "CFL_TRUNCATE_START", cs%CFL_truncS, &
1708  "The start value of the truncation CFL number used when "//&
1709  "ramping up CFL_TRUNC.", &
1710  units="nondim", default=0.)
1711  call get_param(param_file, mdl, "STOKES_MIXING_COMBINED", cs%StokesMixing, &
1712  "Flag to use Stokes drift Mixing via the Lagrangian "//&
1713  " current (Eulerian plus Stokes drift). "//&
1714  " Still needs work and testing, so not recommended for use.",&
1715  default=.false.)
1716  !BGR 04/04/2018{
1717  ! StokesMixing is required for MOM6 for some Langmuir mixing parameterization.
1718  ! The code used here has not been developed for vanishing layers or in
1719  ! conjunction with any bottom friction. Therefore, the following line is
1720  ! added so this functionality cannot be used without user intervention in
1721  ! the code. This will prevent general use of this functionality until proper
1722  ! care is given to the previously mentioned issues. Comment out the following
1723  ! MOM_error to use, but do so at your own risk and with these points in mind.
1724  !}
1725  if (cs%StokesMixing) then
1726  call mom_error(fatal, "Stokes mixing requires user intervention in the code.\n"//&
1727  " Model now exiting. See MOM_vert_friction.F90 for \n"//&
1728  " details (search 'BGR 04/04/2018' to locate comment).")
1729  endif
1730  call get_param(param_file, mdl, "VEL_UNDERFLOW", cs%vel_underflow, &
1731  "A negligibly small velocity magnitude below which velocity "//&
1732  "components are set to 0. A reasonable value might be "//&
1733  "1e-30 m/s, which is less than an Angstrom divided by "//&
1734  "the age of the universe.", units="m s-1", default=0.0)
1735 
1736  alloc_(cs%a_u(isdb:iedb,jsd:jed,nz+1)) ; cs%a_u(:,:,:) = 0.0
1737  alloc_(cs%h_u(isdb:iedb,jsd:jed,nz)) ; cs%h_u(:,:,:) = 0.0
1738  alloc_(cs%a_v(isd:ied,jsdb:jedb,nz+1)) ; cs%a_v(:,:,:) = 0.0
1739  alloc_(cs%h_v(isd:ied,jsdb:jedb,nz)) ; cs%h_v(:,:,:) = 0.0
1740 
1741  cs%id_Kv_slow = register_diag_field('ocean_model', 'Kv_slow', diag%axesTi, time, &
1742  'Slow varying vertical viscosity', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
1743 
1744  cs%id_Kv_u = register_diag_field('ocean_model', 'Kv_u', diag%axesCuL, time, &
1745  'Total vertical viscosity at u-points', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
1746 
1747  cs%id_Kv_v = register_diag_field('ocean_model', 'Kv_v', diag%axesCvL, time, &
1748  'Total vertical viscosity at v-points', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
1749 
1750  cs%id_au_vv = register_diag_field('ocean_model', 'au_visc', diag%axesCui, time, &
1751  'Zonal Viscous Vertical Coupling Coefficient', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
1752 
1753  cs%id_av_vv = register_diag_field('ocean_model', 'av_visc', diag%axesCvi, time, &
1754  'Meridional Viscous Vertical Coupling Coefficient', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
1755 
1756  cs%id_h_u = register_diag_field('ocean_model', 'Hu_visc', diag%axesCuL, time, &
1757  'Thickness at Zonal Velocity Points for Viscosity', thickness_units)
1758 
1759  cs%id_h_v = register_diag_field('ocean_model', 'Hv_visc', diag%axesCvL, time, &
1760  'Thickness at Meridional Velocity Points for Viscosity', thickness_units)
1761 
1762  cs%id_hML_u = register_diag_field('ocean_model', 'HMLu_visc', diag%axesCu1, time, &
1763  'Mixed Layer Thickness at Zonal Velocity Points for Viscosity', thickness_units)
1764 
1765  cs%id_hML_v = register_diag_field('ocean_model', 'HMLv_visc', diag%axesCv1, time, &
1766  'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', thickness_units)
1767 
1768  cs%id_du_dt_visc = register_diag_field('ocean_model', 'du_dt_visc', diag%axesCuL, &
1769  time, 'Zonal Acceleration from Vertical Viscosity', 'm s-2')
1770  if (cs%id_du_dt_visc > 0) call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)
1771  cs%id_dv_dt_visc = register_diag_field('ocean_model', 'dv_dt_visc', diag%axesCvL, &
1772  time, 'Meridional Acceleration from Vertical Viscosity', 'm s-2')
1773  if (cs%id_dv_dt_visc > 0) call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)
1774 
1775  cs%id_taux_bot = register_diag_field('ocean_model', 'taux_bot', diag%axesCu1, &
1776  time, 'Zonal Bottom Stress from Ocean to Earth', 'Pa', &
1777  conversion=us%Z_to_m)
1778  cs%id_tauy_bot = register_diag_field('ocean_model', 'tauy_bot', diag%axesCv1, &
1779  time, 'Meridional Bottom Stress from Ocean to Earth', 'Pa', &
1780  conversion=us%Z_to_m)
1781 
1782  if ((len_trim(cs%u_trunc_file) > 0) .or. (len_trim(cs%v_trunc_file) > 0)) &
1783  call pointaccel_init(mis, time, g, param_file, diag, dirs, cs%PointAccel_CSp)
1784 

◆ 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]gOcean grid structure
[in]gvOcean vertical grid structure
[in]usA dimensional unit scaling type
[in,out]uZonal velocity [m s-1]
[in,out]vMeridional velocity [m s-1]
[in]hLayer thickness [H ~> m or kg m-2]
[in]adpAcceleration diagnostic pointers
[in]cdpContinuity diagnostic pointers
[in]forcesA structure with the driving mechanical forces
[in]viscViscosities and bottom drag
[in]dtTime increment [s]
csVertical viscosity control structure

Definition at line 1355 of file MOM_vert_friction.F90.

1355  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1356  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1357  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1358  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1359  intent(inout) :: u !< Zonal velocity [m s-1]
1360  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1361  intent(inout) :: v !< Meridional velocity [m s-1]
1362  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1363  intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
1364  type(accel_diag_ptrs), intent(in) :: ADp !< Acceleration diagnostic pointers
1365  type(cont_diag_ptrs), intent(in) :: CDp !< Continuity diagnostic pointers
1366  type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
1367  type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag
1368  real, intent(in) :: dt !< Time increment [s]
1369  type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
1370 
1371  ! Local variables
1372 
1373  real :: maxvel ! Velocities components greater than maxvel
1374  real :: truncvel ! are truncated to truncvel, both [m s-1].
1375  real :: CFL ! The local CFL number.
1376  real :: H_report ! A thickness below which not to report truncations.
1377  real :: dt_Rho0 ! The timestep divided by the Boussinesq density [s m3 kg-1].
1378  real :: vel_report(SZIB_(G),SZJB_(G))
1379  real :: u_old(SZIB_(G),SZJ_(G),SZK_(G))
1380  real :: v_old(SZI_(G),SZJB_(G),SZK_(G))
1381  logical :: trunc_any, dowrite(SZIB_(G),SZJB_(G))
1382  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1383  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1384  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1385 
1386  maxvel = cs%maxvel
1387  truncvel = 0.9*maxvel
1388  h_report = 6.0 * gv%Angstrom_H
1389  dt_rho0 = dt / gv%Rho0
1390 
1391  if (len_trim(cs%u_trunc_file) > 0) then
1392  !$OMP parallel do default(shared) private(trunc_any,CFL)
1393  do j=js,je
1394  trunc_any = .false.
1395  do i=isq,ieq ; dowrite(i,j) = .false. ; enddo
1396  if (cs%CFL_based_trunc) then
1397  do i=isq,ieq ; vel_report(i,j) = 3.0e8 ; enddo ! Speed of light default.
1398  do k=1,nz ; do i=isq,ieq
1399  if (abs(u(i,j,k)) < cs%vel_underflow) u(i,j,k) = 0.0
1400  if (u(i,j,k) < 0.0) then
1401  cfl = (-u(i,j,k) * dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
1402  else
1403  cfl = (u(i,j,k) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
1404  endif
1405  if (cfl > cs%CFL_trunc) trunc_any = .true.
1406  if (cfl > cs%CFL_report) then
1407  dowrite(i,j) = .true.
1408  vel_report(i,j) = min(vel_report(i,j), abs(u(i,j,k)))
1409  endif
1410  enddo ; enddo
1411  else
1412  do i=isq,ieq; vel_report(i,j) = maxvel; enddo
1413  do k=1,nz ; do i=isq,ieq
1414  if (abs(u(i,j,k)) < cs%vel_underflow) then ; u(i,j,k) = 0.0
1415  elseif (abs(u(i,j,k)) > maxvel) then
1416  dowrite(i,j) = .true. ; trunc_any = .true.
1417  endif
1418  enddo ; enddo
1419  endif
1420 
1421  do i=isq,ieq ; if (dowrite(i,j)) then
1422  u_old(i,j,:) = u(i,j,:)
1423  endif ; enddo
1424 
1425  if (trunc_any) then ; if (cs%CFL_based_trunc) then
1426  do k=1,nz ; do i=isq,ieq
1427  if ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i+1,j) < -cs%CFL_trunc) then
1428  u(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i+1,j) / (dt * g%dy_Cu(i,j)))
1429  if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1430  elseif ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then
1431  u(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dy_Cu(i,j)))
1432  if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1433  endif
1434  enddo ; enddo
1435  else
1436  do k=1,nz ; do i=isq,ieq ; if (abs(u(i,j,k)) > maxvel) then
1437  u(i,j,k) = sign(truncvel,u(i,j,k))
1438  if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1439  endif ; enddo ; enddo
1440  endif ; endif
1441  enddo ! j-loop
1442  else ! Do not report accelerations leading to large velocities.
1443  if (cs%CFL_based_trunc) then
1444 !$OMP parallel do default(none) shared(nz,js,je,Isq,Ieq,u,dt,G,CS,h,H_report)
1445  do k=1,nz ; do j=js,je ; do i=isq,ieq
1446  if (abs(u(i,j,k)) < cs%vel_underflow) then ; u(i,j,k) = 0.0
1447  elseif ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i+1,j) < -cs%CFL_trunc) then
1448  u(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i+1,j) / (dt * g%dy_Cu(i,j)))
1449  if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1450  elseif ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then
1451  u(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dy_Cu(i,j)))
1452  if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1453  endif
1454  enddo ; enddo ; enddo
1455  else
1456 !$OMP parallel do default(none) shared(nz,js,je,Isq,Ieq,u,G,CS,truncvel,maxvel,h,H_report)
1457  do k=1,nz ; do j=js,je ; do i=isq,ieq
1458  if (abs(u(i,j,k)) < cs%vel_underflow) then ; u(i,j,k) = 0.0
1459  elseif (abs(u(i,j,k)) > maxvel) then
1460  u(i,j,k) = sign(truncvel,u(i,j,k))
1461  if (h(i,j,k) + h(i+1,j,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1462  endif
1463  enddo ; enddo ; enddo
1464  endif
1465  endif
1466 
1467  if (len_trim(cs%u_trunc_file) > 0) then
1468  do j=js,je; do i=isq,ieq ; if (dowrite(i,j)) then
1469 ! Here the diagnostic reporting subroutines are called if
1470 ! unphysically large values were found.
1471  call write_u_accel(i, j, u_old, h, adp, cdp, dt, g, gv, us, cs%PointAccel_CSp, &
1472  vel_report(i,j), forces%taux(i,j)*dt_rho0, a=cs%a_u, hv=cs%h_u)
1473  endif ; enddo ; enddo
1474  endif
1475 
1476  if (len_trim(cs%v_trunc_file) > 0) then
1477  !$OMP parallel do default(shared) private(trunc_any,CFL)
1478  do j=jsq,jeq
1479  trunc_any = .false.
1480  do i=is,ie ; dowrite(i,j) = .false. ; enddo
1481  if (cs%CFL_based_trunc) then
1482  do i=is,ie ; vel_report(i,j) = 3.0e8 ; enddo ! Speed of light default.
1483  do k=1,nz ; do i=is,ie
1484  if (abs(v(i,j,k)) < cs%vel_underflow) v(i,j,k) = 0.0
1485  if (v(i,j,k) < 0.0) then
1486  cfl = (-v(i,j,k) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
1487  else
1488  cfl = (v(i,j,k) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
1489  endif
1490  if (cfl > cs%CFL_trunc) trunc_any = .true.
1491  if (cfl > cs%CFL_report) then
1492  dowrite(i,j) = .true.
1493  vel_report(i,j) = min(vel_report(i,j), abs(v(i,j,k)))
1494  endif
1495  enddo ; enddo
1496  else
1497  do i=is,ie ; vel_report(i,j) = maxvel ; enddo
1498  do k=1,nz ; do i=is,ie
1499  if (abs(v(i,j,k)) < cs%vel_underflow) then ; v(i,j,k) = 0.0
1500  elseif (abs(v(i,j,k)) > maxvel) then
1501  dowrite(i,j) = .true. ; trunc_any = .true.
1502  endif
1503  enddo ; enddo
1504  endif
1505 
1506  do i=is,ie ; if (dowrite(i,j)) then
1507  v_old(i,j,:) = v(i,j,:)
1508  endif ; enddo
1509 
1510  if (trunc_any) then ; if (cs%CFL_based_trunc) then
1511  do k=1,nz; do i=is,ie
1512  if ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j+1) < -cs%CFL_trunc) then
1513  v(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i,j+1) / (dt * g%dx_Cv(i,j)))
1514  if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1515  elseif ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then
1516  v(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dx_Cv(i,j)))
1517  if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1518  endif
1519  enddo ; enddo
1520  else
1521  do k=1,nz ; do i=is,ie ; if (abs(v(i,j,k)) > maxvel) then
1522  v(i,j,k) = sign(truncvel,v(i,j,k))
1523  if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1524  endif ; enddo ; enddo
1525  endif ; endif
1526  enddo ! J-loop
1527  else ! Do not report accelerations leading to large velocities.
1528  if (cs%CFL_based_trunc) then
1529  !$OMP parallel do default(shared)
1530  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
1531  if (abs(v(i,j,k)) < cs%vel_underflow) then ; v(i,j,k) = 0.0
1532  elseif ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j+1) < -cs%CFL_trunc) then
1533  v(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i,j+1) / (dt * g%dx_Cv(i,j)))
1534  if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1535  elseif ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then
1536  v(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dx_Cv(i,j)))
1537  if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1538  endif
1539  enddo ; enddo ; enddo
1540  else
1541  !$OMP parallel do default(shared)
1542  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
1543  if (abs(v(i,j,k)) < cs%vel_underflow) then ; v(i,j,k) = 0.0
1544  elseif (abs(v(i,j,k)) > maxvel) then
1545  v(i,j,k) = sign(truncvel,v(i,j,k))
1546  if (h(i,j,k) + h(i,j+1,k) > h_report) cs%ntrunc = cs%ntrunc + 1
1547  endif
1548  enddo ; enddo ; enddo
1549  endif
1550  endif
1551 
1552  if (len_trim(cs%v_trunc_file) > 0) then
1553  do j=jsq,jeq; do i=is,ie ; if (dowrite(i,j)) then
1554 ! Here the diagnostic reporting subroutines are called if
1555 ! unphysically large values were found.
1556  call write_v_accel(i, j, v_old, h, adp, cdp, dt, g, gv, us, cs%PointAccel_CSp, &
1557  vel_report(i,j), forces%tauy(i,j)*dt_rho0, a=cs%a_v, hv=cs%h_v)
1558  endif ; enddo ; enddo
1559  endif
1560 

◆ 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]gOcean grid structure
[in]gvOcean vertical grid structure
[in]viscViscosities and bottom drag
[in,out]visc_rem_uFraction of a time-step's worth of a
[in,out]visc_rem_vFraction of a time-step's worth of a
[in]dtTime increment [s]
[in]usA dimensional unit scaling type
csVertical viscosity control structure

Definition at line 459 of file MOM_vert_friction.F90.

459  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
460  type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
461  type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag
462  real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
463  intent(inout) :: visc_rem_u !< Fraction of a time-step's worth of a
464  !! barotopic acceleration that a layer experiences after
465  !! viscosity is applied in the zonal direction [nondim]
466  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
467  intent(inout) :: visc_rem_v !< Fraction of a time-step's worth of a
468  !! barotopic acceleration that a layer experiences after
469  !! viscosity is applied in the meridional direction [nondim]
470  real, intent(in) :: dt !< Time increment [s]
471  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
472  type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
473 
474  ! Local variables
475 
476  real :: b1(SZIB_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
477  real :: c1(SZIB_(G),SZK_(G)) ! A variable used by the tridiagonal solver [nondim].
478  real :: d1(SZIB_(G)) ! d1=1-c1 is used by the tridiagonal solver [nondim].
479  real :: Ray(SZIB_(G),SZK_(G)) ! Ray is the Rayleigh-drag velocity [Z T-1 ~> m s-1].
480  real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2].
481  real :: dt_Z_to_H ! The time step times the conversion from Z to the
482  ! units of thickness [T H Z-1 ~> s or s kg m-3].
483  logical :: do_i(SZIB_(G))
484 
485  integer :: i, j, k, is, ie, Isq, Ieq, Jsq, Jeq, nz
486  is = g%isc ; ie = g%iec
487  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
488 
489  if (.not.associated(cs)) call mom_error(fatal,"MOM_vert_friction(visc): "// &
490  "Module must be initialized before it is used.")
491 
492  dt_z_to_h = us%s_to_T*dt*gv%Z_to_H
493 
494  do k=1,nz ; do i=isq,ieq ; ray(i,k) = 0.0 ; enddo ; enddo
495 
496  ! Find the zonal viscous using a modification of a standard tridagonal solver.
497 !$OMP parallel do default(none) shared(G,Isq,Ieq,CS,nz,visc,dt_Z_to_H,visc_rem_u) &
498 !$OMP firstprivate(Ray) &
499 !$OMP private(do_i,b_denom_1,b1,d1,c1)
500  do j=g%jsc,g%jec
501  do i=isq,ieq ; do_i(i) = (g%mask2dCu(i,j) > 0) ; enddo
502 
503  if (cs%Channel_drag) then ; do k=1,nz ; do i=isq,ieq
504  ray(i,k) = visc%Ray_u(i,j,k)
505  enddo ; enddo ; endif
506 
507  do i=isq,ieq ; if (do_i(i)) then
508  b_denom_1 = cs%h_u(i,j,1) + dt_z_to_h * (ray(i,1) + cs%a_u(i,j,1))
509  b1(i) = 1.0 / (b_denom_1 + dt_z_to_h*cs%a_u(i,j,2))
510  d1(i) = b_denom_1 * b1(i)
511  visc_rem_u(i,j,1) = b1(i) * cs%h_u(i,j,1)
512  endif ; enddo
513  do k=2,nz ; do i=isq,ieq ; if (do_i(i)) then
514  c1(i,k) = dt_z_to_h * cs%a_u(i,j,k)*b1(i)
515  b_denom_1 = cs%h_u(i,j,k) + dt_z_to_h * (ray(i,k) + cs%a_u(i,j,k)*d1(i))
516  b1(i) = 1.0 / (b_denom_1 + dt_z_to_h * cs%a_u(i,j,k+1))
517  d1(i) = b_denom_1 * b1(i)
518  visc_rem_u(i,j,k) = (cs%h_u(i,j,k) + dt_z_to_h * cs%a_u(i,j,k) * visc_rem_u(i,j,k-1)) * b1(i)
519  endif ; enddo ; enddo
520  do k=nz-1,1,-1 ; do i=isq,ieq ; if (do_i(i)) then
521  visc_rem_u(i,j,k) = visc_rem_u(i,j,k) + c1(i,k+1)*visc_rem_u(i,j,k+1)
522 
523  endif ; enddo ; enddo ! i and k loops
524 
525  enddo ! end u-component j loop
526 
527  ! Now find the meridional viscous using a modification.
528 !$OMP parallel do default(none) shared(Jsq,Jeq,is,ie,G,CS,visc,dt_Z_to_H,visc_rem_v,nz) &
529 !$OMP firstprivate(Ray) &
530 !$OMP private(do_i,b_denom_1,b1,d1,c1)
531  do j=jsq,jeq
532  do i=is,ie ; do_i(i) = (g%mask2dCv(i,j) > 0) ; enddo
533 
534  if (cs%Channel_drag) then ; do k=1,nz ; do i=is,ie
535  ray(i,k) = visc%Ray_v(i,j,k)
536  enddo ; enddo ; endif
537 
538  do i=is,ie ; if (do_i(i)) then
539  b_denom_1 = cs%h_v(i,j,1) + dt_z_to_h * (ray(i,1) + cs%a_v(i,j,1))
540  b1(i) = 1.0 / (b_denom_1 + dt_z_to_h*cs%a_v(i,j,2))
541  d1(i) = b_denom_1 * b1(i)
542  visc_rem_v(i,j,1) = b1(i) * cs%h_v(i,j,1)
543  endif ; enddo
544  do k=2,nz ; do i=is,ie ; if (do_i(i)) then
545  c1(i,k) = dt_z_to_h * cs%a_v(i,j,k)*b1(i)
546  b_denom_1 = cs%h_v(i,j,k) + dt_z_to_h * (ray(i,k) + cs%a_v(i,j,k)*d1(i))
547  b1(i) = 1.0 / (b_denom_1 + dt_z_to_h * cs%a_v(i,j,k+1))
548  d1(i) = b_denom_1 * b1(i)
549  visc_rem_v(i,j,k) = (cs%h_v(i,j,k) + dt_z_to_h * cs%a_v(i,j,k) * visc_rem_v(i,j,k-1)) * b1(i)
550  endif ; enddo ; enddo
551  do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then
552  visc_rem_v(i,j,k) = visc_rem_v(i,j,k) + c1(i,k+1)*visc_rem_v(i,j,k+1)
553  endif ; enddo ; enddo ! i and k loops
554  enddo ! end of v-component J loop
555 
556  if (cs%debug) then
557  call uvchksum("visc_rem_[uv]", visc_rem_u, visc_rem_v, g%HI,haloshift=0)
558  endif
559