MOM6
Data Types
mom_set_diffusivity Module Reference

Detailed Description

Calculate vertical diffusivity from all mixing processes.

Data Types

type  diffusivity_diags
 This structure has memory for used in calculating diagnostics of diffusivity. More...
 
type  set_diffusivity_cs
 This control structure contains parameters for MOM_set_diffusivity. More...
 
integer id_clock_kappashear
 CPU time clocks. More...
 
integer id_clock_cvmix_ddiff
 CPU time clocks. More...
 
subroutine, public set_diffusivity (u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, G, GV, US, CS, Kd_lay, Kd_int)
 Sets the interior vertical diffusion of scalars due to the following processes: More...
 
subroutine find_tke_to_kd (h, tv, dRho_int, N2_lay, j, dt, G, GV, US, CS, TKE_to_Kd, maxTKE, kb)
 Convert turbulent kinetic energy to diffusivity. More...
 
subroutine find_n2 (h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, N2_lay, N2_int, N2_bot)
 Calculate Brunt-Vaisala frequency, N^2. More...
 
subroutine double_diffusion (tv, h, T_f, S_f, j, G, GV, US, CS, Kd_T_dd, Kd_S_dd)
 This subroutine sets the additional diffusivities of temperature and salinity due to double diffusion, using the same functional form as is used in MOM4.1, and taken from an NCAR technical note (REF?) that updates what was in Large et al. (1994). All the coefficients here should probably be made run-time variables rather than hard-coded constants. More...
 
subroutine add_drag_diffusivity (h, u, v, tv, fluxes, visc, j, TKE_to_Kd, maxTKE, kb, G, GV, US, CS, Kd_lay, Kd_int, Kd_BBL)
 This routine adds diffusion sustained by flow energy extracted by bottom drag. More...
 
subroutine add_lotw_bbl_diffusivity (h, u, v, tv, fluxes, visc, j, N2_int, G, GV, US, CS, Kd_lay, Kd_int, Kd_BBL)
 Calculates a BBL diffusivity use a Prandtl number 1 diffusivitiy with a law of the wall turbulent viscosity, up to a BBL height where the energy used for mixing has consumed the mechanical TKE input. More...
 
subroutine add_mlrad_diffusivity (h, fluxes, j, G, GV, US, CS, Kd_lay, TKE_to_Kd, Kd_int)
 This routine adds effects of mixed layer radiation to the layer diffusivities. More...
 
subroutine, public set_bbl_tke (u, v, h, fluxes, visc, G, GV, US, CS)
 This subroutine calculates several properties related to bottom boundary layer turbulence. More...
 
subroutine set_density_ratios (h, tv, kb, G, GV, US, CS, j, ds_dsp1, rho_0)
 
subroutine, public set_diffusivity_init (Time, G, GV, US, param_file, diag, CS, int_tide_CSp, tm_CSp, halo_TS)
 
subroutine, public set_diffusivity_end (CS)
 Clear pointers and dealocate memory. More...
 

Function/Subroutine Documentation

◆ add_drag_diffusivity()

subroutine mom_set_diffusivity::add_drag_diffusivity ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(in)  v,
type(thermo_var_ptrs), intent(in)  tv,
type(forcing), intent(in)  fluxes,
type(vertvisc_type), intent(in)  visc,
integer, intent(in)  j,
real, dimension( g %isd: g %ied, g %ke), intent(in)  TKE_to_Kd,
real, dimension( g %isd: g %ied, g %ke), intent(in)  maxTKE,
integer, dimension( g %isd: g %ied), intent(in)  kb,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(set_diffusivity_cs), pointer  CS,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  Kd_lay,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke+1), intent(inout)  Kd_int,
real, dimension(:,:,:), pointer  Kd_BBL 
)
private

This routine adds diffusion sustained by flow energy extracted by bottom drag.

Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]usA dimensional unit scaling type
[in]uThe zonal velocity [L T-1 ~> m s-1]
[in]vThe meridional velocity [L T-1 ~> m s-1]
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]tvStructure containing pointers to any available thermodynamic fields.
[in]fluxesA structure of thermodynamic surface fluxes
[in]viscStructure containing vertical viscosities, bottom boundary layer properies, and related fields
[in]jj-index of row to work on
[in]tke_to_kdThe conversion rate between the TKE TKE dissipated within a layer and the diapycnal diffusivity witin that layer, usually (~Rho_0 / (G_Earth * dRho_lay)) [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1]
[in]maxtkeThe energy required to for a layer to entrain to its maximum-realizable thickness [m3 T-3 ~> m3 s-3]
[in]kbIndex of lightest layer denser than the buffer layer, or -1 without a bulk mixed layer
csDiffusivity control structure
[in,out]kd_layThe diapycnal diffusvity in layers,
[in,out]kd_intThe diapycnal diffusvity at interfaces,
kd_bblInterface BBL diffusivity [Z2 T-1 ~> m2 s-1].

Definition at line 1100 of file MOM_set_diffusivity.F90.

1100  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1101  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1102  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1103  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1104  intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1]
1105  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1106  intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]
1107  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1108  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1109  type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
1110  !! thermodynamic fields.
1111  type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes
1112  type(vertvisc_type), intent(in) :: visc !< Structure containing vertical viscosities, bottom
1113  !! boundary layer properies, and related fields
1114  integer, intent(in) :: j !< j-index of row to work on
1115  real, dimension(SZI_(G),SZK_(G)), intent(in) :: TKE_to_Kd !< The conversion rate between the TKE
1116  !! TKE dissipated within a layer and the
1117  !! diapycnal diffusivity witin that layer,
1118  !! usually (~Rho_0 / (G_Earth * dRho_lay))
1119  !! [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1]
1120  real, dimension(SZI_(G),SZK_(G)), intent(in) :: maxTKE !< The energy required to for a layer to entrain
1121  !! to its maximum-realizable thickness [m3 T-3 ~> m3 s-3]
1122  integer, dimension(SZI_(G)), intent(in) :: kb !< Index of lightest layer denser than the buffer
1123  !! layer, or -1 without a bulk mixed layer
1124  type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure
1125  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1126  intent(inout) :: Kd_lay !< The diapycnal diffusvity in layers,
1127  !! [Z2 T-1 ~> m2 s-1].
1128  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
1129  intent(inout) :: Kd_int !< The diapycnal diffusvity at interfaces,
1130  !! [Z2 T-1 ~> m2 s-1].
1131  real, dimension(:,:,:), pointer :: Kd_BBL !< Interface BBL diffusivity [Z2 T-1 ~> m2 s-1].
1132 
1133 ! This routine adds diffusion sustained by flow energy extracted by bottom drag.
1134 
1135  real, dimension(SZK_(G)+1) :: &
1136  Rint ! coordinate density of an interface [R ~> kg m-3]
1137  real, dimension(SZI_(G)) :: &
1138  htot, & ! total thickness above or below a layer, or the
1139  ! integrated thickness in the BBL [Z ~> m].
1140  rho_htot, & ! running integral with depth of density [Z R ~> kg m-2]
1141  gh_sum_top, & ! BBL value of g'h that can be supported by
1142  ! the local ustar, times R0_g [R ~> kg m-2]
1143  rho_top, & ! density at top of the BBL [R ~> kg m-3]
1144  tke, & ! turbulent kinetic energy available to drive
1145  ! bottom-boundary layer mixing in a layer [Z3 T-3 ~> m3 s-3]
1146  i2decay ! inverse of twice the TKE decay scale [Z-1 ~> m-1].
1147 
1148  real :: TKE_to_layer ! TKE used to drive mixing in a layer [Z3 T-3 ~> m3 s-3]
1149  real :: TKE_Ray ! TKE from layer Rayleigh drag used to drive mixing in layer [Z3 T-3 ~> m3 s-3]
1150  real :: TKE_here ! TKE that goes into mixing in this layer [Z3 T-3 ~> m3 s-3]
1151  real :: dRl, dRbot ! temporaries holding density differences [R ~> kg m-3]
1152  real :: cdrag_sqrt ! square root of the drag coefficient [nondim]
1153  real :: ustar_h ! value of ustar at a thickness point [Z T-1 ~> m s-1].
1154  real :: absf ! average absolute Coriolis parameter around a thickness point [T-1 ~> s-1]
1155  real :: R0_g ! Rho0 / G_Earth [R T2 Z-1 m-1 ~> kg s2 m-5]
1156  real :: I_rho0 ! 1 / RHO0 [R-1 ~> m3 kg-1]
1157  real :: delta_Kd ! increment to Kd from the bottom boundary layer mixing [Z2 T-1 ~> m2 s-1].
1158  logical :: Rayleigh_drag ! Set to true if Rayleigh drag velocities
1159  ! defined in visc, on the assumption that this
1160  ! extracted energy also drives diapycnal mixing.
1161 
1162  logical :: domore, do_i(SZI_(G))
1163  logical :: do_diag_Kd_BBL
1164 
1165  integer :: i, k, is, ie, nz, i_rem, kb_min
1166  is = g%isc ; ie = g%iec ; nz = g%ke
1167 
1168  do_diag_kd_bbl = associated(kd_bbl)
1169 
1170  if (.not.(cs%bottomdraglaw .and. (cs%BBL_effic>0.0))) return
1171 
1172  cdrag_sqrt = sqrt(cs%cdrag)
1173  tke_ray = 0.0 ; rayleigh_drag = .false.
1174  if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) rayleigh_drag = .true.
1175 
1176  i_rho0 = 1.0 / (gv%Rho0)
1177  r0_g = gv%Rho0 / (us%L_to_Z**2 * gv%g_Earth)
1178 
1179  do k=2,nz ; rint(k) = 0.5*(gv%Rlay(k-1)+gv%Rlay(k)) ; enddo
1180 
1181  kb_min = max(gv%nk_rho_varies+1,2)
1182 
1183  ! The turbulence decay scale is 0.5*ustar/f from K&E & MOM_vertvisc.F90
1184  ! Any turbulence that makes it into the mixed layers is assumed
1185  ! to be relatively small and is discarded.
1186  do i=is,ie
1187  ustar_h = visc%ustar_BBL(i,j)
1188  if (associated(fluxes%ustar_tidal)) &
1189  ustar_h = ustar_h + fluxes%ustar_tidal(i,j)
1190  absf = 0.25 * ((abs(g%CoriolisBu(i-1,j-1)) + abs(g%CoriolisBu(i,j))) + &
1191  (abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j-1))))
1192  if ((ustar_h > 0.0) .and. (absf > 0.5*cs%IMax_decay*ustar_h)) then
1193  i2decay(i) = absf / ustar_h
1194  else
1195  ! The maximum decay scale should be something of order 200 m.
1196  ! If ustar_h = 0, this is land so this value doesn't matter.
1197  i2decay(i) = 0.5*cs%IMax_decay
1198  endif
1199  tke(i) = ((cs%BBL_effic * cdrag_sqrt) * exp(-i2decay(i)*(gv%H_to_Z*h(i,j,nz))) ) * &
1200  visc%TKE_BBL(i,j)
1201 
1202  if (associated(fluxes%TKE_tidal)) &
1203  tke(i) = tke(i) + fluxes%TKE_tidal(i,j) * i_rho0 * &
1204  (cs%BBL_effic * exp(-i2decay(i)*(gv%H_to_Z*h(i,j,nz))))
1205 
1206  ! Distribute the work over a BBL of depth 20^2 ustar^2 / g' following
1207  ! Killworth & Edwards (1999) and Zilitikevich & Mironov (1996).
1208  ! Rho_top is determined by finding the density where
1209  ! integral(bottom, Z) (rho(z') - rho(Z)) dz' = rho_0 400 ustar^2 / g
1210 
1211  gh_sum_top(i) = r0_g * 400.0 * ustar_h**2
1212 
1213  do_i(i) = (g%mask2dT(i,j) > 0.5)
1214  htot(i) = gv%H_to_Z*h(i,j,nz)
1215  rho_htot(i) = gv%Rlay(nz)*(gv%H_to_Z*h(i,j,nz))
1216  rho_top(i) = gv%Rlay(1)
1217  if (cs%bulkmixedlayer .and. do_i(i)) rho_top(i) = gv%Rlay(kb(i)-1)
1218  enddo
1219 
1220  do k=nz-1,2,-1 ; domore = .false.
1221  do i=is,ie ; if (do_i(i)) then
1222  htot(i) = htot(i) + gv%H_to_Z*h(i,j,k)
1223  rho_htot(i) = rho_htot(i) + gv%Rlay(k)*(gv%H_to_Z*h(i,j,k))
1224  if (htot(i)*gv%Rlay(k-1) <= (rho_htot(i) - gh_sum_top(i))) then
1225  ! The top of the mixing is in the interface atop the current layer.
1226  rho_top(i) = (rho_htot(i) - gh_sum_top(i)) / htot(i)
1227  do_i(i) = .false.
1228  elseif (k <= kb(i)) then ; do_i(i) = .false.
1229  else ; domore = .true. ; endif
1230  endif ; enddo
1231  if (.not.domore) exit
1232  enddo ! k-loop
1233 
1234  do i=is,ie ; do_i(i) = (g%mask2dT(i,j) > 0.5) ; enddo
1235  do k=nz-1,kb_min,-1
1236  i_rem = 0
1237  do i=is,ie ; if (do_i(i)) then
1238  if (k<kb(i)) then ; do_i(i) = .false. ; cycle ; endif
1239  i_rem = i_rem + 1 ! Count the i-rows that are still being worked on.
1240  ! Apply vertical decay of the turbulent energy. This energy is
1241  ! simply lost.
1242  tke(i) = tke(i) * exp(-i2decay(i) * (gv%H_to_Z*(h(i,j,k) + h(i,j,k+1))))
1243 
1244 ! if (maxEnt(i,k) <= 0.0) cycle
1245  if (maxtke(i,k) <= 0.0) cycle
1246 
1247  ! This is an analytic integral where diffusity is a quadratic function of
1248  ! rho that goes asymptotically to 0 at Rho_top (vaguely following KPP?).
1249  if (tke(i) > 0.0) then
1250  if (rint(k) <= rho_top(i)) then
1251  tke_to_layer = tke(i)
1252  else
1253  drl = rint(k+1) - rint(k) ; drbot = rint(k+1) - rho_top(i)
1254  tke_to_layer = tke(i) * drl * &
1255  (3.0*drbot*(rint(k) - rho_top(i)) + drl**2) / (drbot**3)
1256  endif
1257  else ; tke_to_layer = 0.0 ; endif
1258 
1259  ! TKE_Ray has been initialized to 0 above.
1260  if (rayleigh_drag) tke_ray = 0.5*cs%BBL_effic * us%L_to_Z**2 * g%IareaT(i,j) * &
1261  ((g%areaCu(i-1,j) * visc%Ray_u(i-1,j,k) * u(i-1,j,k)**2 + &
1262  g%areaCu(i,j) * visc%Ray_u(i,j,k) * u(i,j,k)**2) + &
1263  (g%areaCv(i,j-1) * visc%Ray_v(i,j-1,k) * v(i,j-1,k)**2 + &
1264  g%areaCv(i,j) * visc%Ray_v(i,j,k) * v(i,j,k)**2))
1265 
1266  if (tke_to_layer + tke_ray > 0.0) then
1267  if (cs%BBL_mixing_as_max) then
1268  if (tke_to_layer + tke_ray > maxtke(i,k)) &
1269  tke_to_layer = maxtke(i,k) - tke_ray
1270 
1271  tke(i) = tke(i) - tke_to_layer
1272 
1273  if (kd_lay(i,j,k) < (tke_to_layer + tke_ray) * tke_to_kd(i,k)) then
1274  delta_kd = (tke_to_layer + tke_ray) * tke_to_kd(i,k) - kd_lay(i,j,k)
1275  if ((cs%Kd_max >= 0.0) .and. (delta_kd > cs%Kd_max)) then
1276  delta_kd = cs%Kd_max
1277  kd_lay(i,j,k) = kd_lay(i,j,k) + delta_kd
1278  else
1279  kd_lay(i,j,k) = (tke_to_layer + tke_ray) * tke_to_kd(i,k)
1280  endif
1281  kd_int(i,j,k) = kd_int(i,j,k) + 0.5 * delta_kd
1282  kd_int(i,j,k+1) = kd_int(i,j,k+1) + 0.5 * delta_kd
1283  if (do_diag_kd_bbl) then
1284  kd_bbl(i,j,k) = kd_bbl(i,j,k) + 0.5 * delta_kd
1285  kd_bbl(i,j,k+1) = kd_bbl(i,j,k+1) + 0.5 * delta_kd
1286  endif
1287  endif
1288  else
1289  if (kd_lay(i,j,k) >= maxtke(i,k) * tke_to_kd(i,k)) then
1290  tke_here = 0.0
1291  tke(i) = tke(i) + tke_ray
1292  elseif (kd_lay(i,j,k) + (tke_to_layer + tke_ray) * tke_to_kd(i,k) > &
1293  maxtke(i,k) * tke_to_kd(i,k)) then
1294  tke_here = ((tke_to_layer + tke_ray) + kd_lay(i,j,k) / tke_to_kd(i,k)) - maxtke(i,k)
1295  tke(i) = (tke(i) - tke_here) + tke_ray
1296  else
1297  tke_here = tke_to_layer + tke_ray
1298  tke(i) = tke(i) - tke_to_layer
1299  endif
1300  if (tke(i) < 0.0) tke(i) = 0.0 ! This should be unnecessary?
1301 
1302  if (tke_here > 0.0) then
1303  delta_kd = tke_here * tke_to_kd(i,k)
1304  if (cs%Kd_max >= 0.0) delta_kd = min(delta_kd, cs%Kd_max)
1305  kd_lay(i,j,k) = kd_lay(i,j,k) + delta_kd
1306  kd_int(i,j,k) = kd_int(i,j,k) + 0.5 * delta_kd
1307  kd_int(i,j,k+1) = kd_int(i,j,k+1) + 0.5 * delta_kd
1308  if (do_diag_kd_bbl) then
1309  kd_bbl(i,j,k) = kd_bbl(i,j,k) + 0.5 * delta_kd
1310  kd_bbl(i,j,k+1) = kd_bbl(i,j,k+1) + 0.5 * delta_kd
1311  endif
1312  endif
1313  endif
1314  endif
1315 
1316  ! This may be risky - in the case that there are exactly zero
1317  ! velocities at 4 neighboring points, but nonzero velocities
1318  ! above the iterations would stop too soon. I don't see how this
1319  ! could happen in practice. RWH
1320  if ((tke(i)<= 0.0) .and. (tke_ray == 0.0)) then
1321  do_i(i) = .false. ; i_rem = i_rem - 1
1322  endif
1323 
1324  endif ; enddo
1325  if (i_rem == 0) exit
1326  enddo ! k-loop
1327 

Referenced by set_diffusivity().

Here is the caller graph for this function:

◆ add_lotw_bbl_diffusivity()

subroutine mom_set_diffusivity::add_lotw_bbl_diffusivity ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(in)  v,
type(thermo_var_ptrs), intent(in)  tv,
type(forcing), intent(in)  fluxes,
type(vertvisc_type), intent(in)  visc,
integer, intent(in)  j,
real, dimension( g %isd: g %ied, g %ke+1), intent(in)  N2_int,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(set_diffusivity_cs), pointer  CS,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  Kd_lay,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke+1), intent(inout)  Kd_int,
real, dimension(:,:,:), pointer  Kd_BBL 
)
private

Calculates a BBL diffusivity use a Prandtl number 1 diffusivitiy with a law of the wall turbulent viscosity, up to a BBL height where the energy used for mixing has consumed the mechanical TKE input.

Parameters
[in]gGrid structure
[in]gvVertical grid structure
[in]usA dimensional unit scaling type
[in]uu component of flow [L T-1 ~> m s-1]
[in]vv component of flow [L T-1 ~> m s-1]
[in]hLayer thickness [H ~> m or kg m-2]
[in]tvStructure containing pointers to any available thermodynamic fields.
[in]fluxesSurface fluxes structure
[in]viscStructure containing vertical viscosities, bottom boundary layer properies, and related fields.
[in]jj-index of row to work on
[in]n2_intSquare of Brunt-Vaisala at interfaces [T-2 ~> s-2]
csDiffusivity control structure
[in,out]kd_layLayer net diffusivity [Z2 T-1 ~> m2 s-1]
[in,out]kd_intInterface net diffusivity [Z2 T-1 ~> m2 s-1]
kd_bblInterface BBL diffusivity [Z2 T-1 ~> m2 s-1]

Definition at line 1335 of file MOM_set_diffusivity.F90.

1335  type(ocean_grid_type), intent(in) :: G !< Grid structure
1336  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1337  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1338  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1339  intent(in) :: u !< u component of flow [L T-1 ~> m s-1]
1340  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1341  intent(in) :: v !< v component of flow [L T-1 ~> m s-1]
1342  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1343  intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
1344  type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
1345  !! thermodynamic fields.
1346  type(forcing), intent(in) :: fluxes !< Surface fluxes structure
1347  type(vertvisc_type), intent(in) :: visc !< Structure containing vertical viscosities, bottom
1348  !! boundary layer properies, and related fields.
1349  integer, intent(in) :: j !< j-index of row to work on
1350  real, dimension(SZI_(G),SZK_(G)+1), &
1351  intent(in) :: N2_int !< Square of Brunt-Vaisala at interfaces [T-2 ~> s-2]
1352  type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure
1353  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1354  intent(inout) :: Kd_lay !< Layer net diffusivity [Z2 T-1 ~> m2 s-1]
1355  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
1356  intent(inout) :: Kd_int !< Interface net diffusivity [Z2 T-1 ~> m2 s-1]
1357  real, dimension(:,:,:), pointer :: Kd_BBL !< Interface BBL diffusivity [Z2 T-1 ~> m2 s-1]
1358 
1359  ! Local variables
1360  real :: TKE_column ! net TKE input into the column [Z3 T-3 ~> m3 s-3]
1361  real :: TKE_remaining ! remaining TKE available for mixing in this layer and above [Z3 T-3 ~> m3 s-3]
1362  real :: TKE_consumed ! TKE used for mixing in this layer [Z3 T-3 ~> m3 s-3]
1363  real :: TKE_Kd_wall ! TKE associated with unlimited law of the wall mixing [Z3 T-3 ~> m3 s-3]
1364  real :: cdrag_sqrt ! square root of the drag coefficient [nondim]
1365  real :: ustar ! value of ustar at a thickness point [Z T-1 ~> m s-1].
1366  real :: ustar2 ! square of ustar, for convenience [Z2 T-2 ~> m2 s-2]
1367  real :: absf ! average absolute value of Coriolis parameter around a thickness point [T-1 ~> s-1]
1368  real :: dh, dhm1 ! thickness of layers k and k-1, respecitvely [Z ~> m].
1369  real :: z_bot ! distance to interface k from bottom [Z ~> m].
1370  real :: D_minus_z ! distance to interface k from surface [Z ~> m].
1371  real :: total_thickness ! total thickness of water column [Z ~> m].
1372  real :: Idecay ! inverse of decay scale used for "Joule heating" loss of TKE with height [Z-1 ~> m-1].
1373  real :: Kd_wall ! Law of the wall diffusivity [Z2 T-1 ~> m2 s-1].
1374  real :: Kd_lower ! diffusivity for lower interface [Z2 T-1 ~> m2 s-1]
1375  real :: ustar_D ! u* x D [Z2 T-1 ~> m2 s-1].
1376  real :: I_Rho0 ! 1 / rho0 [R-1 ~> m3 kg-1]
1377  real :: N2_min ! Minimum value of N2 to use in calculation of TKE_Kd_wall [T-2 ~> s-2]
1378  logical :: Rayleigh_drag ! Set to true if there are Rayleigh drag velocities defined in visc, on
1379  ! the assumption that this extracted energy also drives diapycnal mixing.
1380  integer :: i, k, km1
1381  real, parameter :: von_karm = 0.41 ! Von Karman constant (http://en.wikipedia.org/wiki/Von_Karman_constant)
1382  logical :: do_diag_Kd_BBL
1383 
1384  if (.not.(cs%bottomdraglaw .and. (cs%BBL_effic>0.0))) return
1385  do_diag_kd_bbl = associated(kd_bbl)
1386 
1387  n2_min = 0.
1388  if (cs%LOTW_BBL_use_omega) n2_min = cs%omega**2
1389 
1390  ! Determine whether to add Rayleigh drag contribution to TKE
1391  rayleigh_drag = .false.
1392  if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) rayleigh_drag = .true.
1393  i_rho0 = 1.0 / (gv%Rho0)
1394  cdrag_sqrt = sqrt(cs%cdrag)
1395 
1396  do i=g%isc,g%iec ! Developed in single-column mode
1397 
1398  ! Column-wise parameters.
1399  absf = 0.25 * ((abs(g%CoriolisBu(i-1,j-1)) + abs(g%CoriolisBu(i,j))) + &
1400  (abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j-1)))) ! Non-zero on equator!
1401 
1402  ! u* at the bottom [Z T-1 ~> m s-1].
1403  ustar = visc%ustar_BBL(i,j)
1404  ustar2 = ustar**2
1405  ! In add_drag_diffusivity(), fluxes%ustar_tidal is added in. This might be double counting
1406  ! since ustar_BBL should already include all contributions to u*? -AJA
1407  !### Examine this question of whether there is double counting of fluxes%ustar_tidal.
1408  if (associated(fluxes%ustar_tidal)) ustar = ustar + fluxes%ustar_tidal(i,j)
1409 
1410  ! The maximum decay scale should be something of order 200 m. We use the smaller of u*/f and
1411  ! (IMax_decay)^-1 as the decay scale. If ustar = 0, this is land so this value doesn't matter.
1412  idecay = cs%IMax_decay
1413  if ((ustar > 0.0) .and. (absf > cs%IMax_decay * ustar)) idecay = absf / ustar
1414 
1415  ! Energy input at the bottom [Z3 T-3 ~> m3 s-3].
1416  ! (Note that visc%TKE_BBL is in [Z3 T-3 ~> m3 s-3], set in set_BBL_TKE().)
1417  ! I am still unsure about sqrt(cdrag) in this expressions - AJA
1418  tke_column = cdrag_sqrt * visc%TKE_BBL(i,j)
1419  ! Add in tidal dissipation energy at the bottom [R Z3 T-3 ~> m3 s-3].
1420  ! Note that TKE_tidal is in [R Z3 T-3 ~> W m-2].
1421  if (associated(fluxes%TKE_tidal)) &
1422  tke_column = tke_column + fluxes%TKE_tidal(i,j) * i_rho0
1423  tke_column = cs%BBL_effic * tke_column ! Only use a fraction of the mechanical dissipation for mixing.
1424 
1425  tke_remaining = tke_column
1426  total_thickness = ( sum(h(i,j,:)) + gv%H_subroundoff )* gv%H_to_Z ! Total column thickness [Z ~> m].
1427  ustar_d = ustar * total_thickness
1428  z_bot = 0.
1429  kd_lower = 0. ! Diffusivity on bottom boundary.
1430 
1431  ! Work upwards from the bottom, accumulating work used until it exceeds the available TKE input
1432  ! at the bottom.
1433  do k=g%ke,2,-1
1434  dh = gv%H_to_Z * h(i,j,k) ! Thickness of this level [Z ~> m].
1435  km1 = max(k-1, 1)
1436  dhm1 = gv%H_to_Z * h(i,j,km1) ! Thickness of level above [Z ~> m].
1437 
1438  ! Add in additional energy input from bottom-drag against slopes (sides)
1439  if (rayleigh_drag) tke_remaining = tke_remaining + &
1440  0.5*cs%BBL_effic * us%L_to_Z**2 * g%IareaT(i,j) * &
1441  ((g%areaCu(i-1,j) * visc%Ray_u(i-1,j,k) * u(i-1,j,k)**2 + &
1442  g%areaCu(i,j) * visc%Ray_u(i,j,k) * u(i,j,k)**2) + &
1443  (g%areaCv(i,j-1) * visc%Ray_v(i,j-1,k) * v(i,j-1,k)**2 + &
1444  g%areaCv(i,j) * visc%Ray_v(i,j,k) * v(i,j,k)**2))
1445 
1446  ! Exponentially decay TKE across the thickness of the layer.
1447  ! This is energy loss in addition to work done as mixing, apparently to Joule heating.
1448  tke_remaining = exp(-idecay*dh) * tke_remaining
1449 
1450  z_bot = z_bot + h(i,j,k)*gv%H_to_Z ! Distance between upper interface of layer and the bottom [Z ~> m].
1451  d_minus_z = max(total_thickness - z_bot, 0.) ! Thickness above layer [Z ~> m].
1452 
1453  ! Diffusivity using law of the wall, limited by rotation, at height z [Z2 T-1 ~> m2 s-1].
1454  ! This calculation is at the upper interface of the layer
1455  if ( ustar_d + absf * ( z_bot * d_minus_z ) == 0.) then
1456  kd_wall = 0.
1457  else
1458  kd_wall = ((von_karm * ustar2) * (z_bot * d_minus_z)) &
1459  / (ustar_d + absf * (z_bot * d_minus_z))
1460  endif
1461 
1462  ! TKE associated with Kd_wall [Z3 T-3 ~> m3 s-3].
1463  ! This calculation if for the volume spanning the interface.
1464  tke_kd_wall = kd_wall * 0.5 * (dh + dhm1) * max(n2_int(i,k), n2_min)
1465 
1466  ! Now bound Kd such that the associated TKE is no greater than available TKE for mixing.
1467  if (tke_kd_wall > 0.) then
1468  tke_consumed = min(tke_kd_wall, tke_remaining)
1469  kd_wall = (tke_consumed / tke_kd_wall) * kd_wall ! Scale Kd so that only TKE_consumed is used.
1470  else
1471  ! Either N2=0 or dh = 0.
1472  if (tke_remaining > 0.) then
1473  kd_wall = cs%Kd_max
1474  else
1475  kd_wall = 0.
1476  endif
1477  tke_consumed = 0.
1478  endif
1479 
1480  ! Now use up the appropriate about of TKE associated with the diffusivity chosen
1481  tke_remaining = tke_remaining - tke_consumed ! Note this will be non-negative
1482 
1483  ! Add this BBL diffusivity to the model net diffusivity.
1484  kd_int(i,j,k) = kd_int(i,j,k) + kd_wall
1485  kd_lay(i,j,k) = kd_lay(i,j,k) + 0.5 * (kd_wall + kd_lower)
1486  kd_lower = kd_wall ! Store for next level up.
1487  if (do_diag_kd_bbl) kd_bbl(i,j,k) = kd_wall
1488  enddo ! k
1489  enddo ! i
1490 

Referenced by set_diffusivity().

Here is the caller graph for this function:

◆ add_mlrad_diffusivity()

subroutine mom_set_diffusivity::add_mlrad_diffusivity ( real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
type(forcing), intent(in)  fluxes,
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(set_diffusivity_cs), pointer  CS,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  Kd_lay,
real, dimension(szi_(g),szk_(g)), intent(in)  TKE_to_Kd,
real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(inout), optional  Kd_int 
)
private

This routine adds effects of mixed layer radiation to the layer diffusivities.

Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]usA dimensional unit scaling type
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]fluxesSurface fluxes structure
csDiffusivity control structure
[in,out]kd_layThe diapycnal diffusvity in layers [Z2 T-1 ~> m2 s-1].
[in]jThe j-index to work on
[in]tke_to_kdThe conversion rate between the TKE TKE dissipated within a layer and the diapycnal diffusivity witin that layer, usually (~Rho_0 / (G_Earth * dRho_lay)) [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1]
[in,out]kd_intThe diapycnal diffusvity at interfaces

Definition at line 1495 of file MOM_set_diffusivity.F90.

1495  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1496  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1497  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1498  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1499  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1500  type(forcing), intent(in) :: fluxes !< Surface fluxes structure
1501  type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure
1502  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1503  intent(inout) :: Kd_lay !< The diapycnal diffusvity in layers [Z2 T-1 ~> m2 s-1].
1504  integer, intent(in) :: j !< The j-index to work on
1505  real, dimension(SZI_(G),SZK_(G)), intent(in) :: TKE_to_Kd !< The conversion rate between the TKE
1506  !! TKE dissipated within a layer and the
1507  !! diapycnal diffusivity witin that layer,
1508  !! usually (~Rho_0 / (G_Earth * dRho_lay))
1509  !! [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1]
1510  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
1511  optional, intent(inout) :: Kd_int !< The diapycnal diffusvity at interfaces
1512  !! [Z2 T-1 ~> m2 s-1].
1513 
1514 ! This routine adds effects of mixed layer radiation to the layer diffusivities.
1515 
1516  real, dimension(SZI_(G)) :: h_ml ! Mixed layer thickness [Z ~> m].
1517  real, dimension(SZI_(G)) :: TKE_ml_flux ! Mixed layer TKE flux [Z3 T-3 ~> m3 s-3]
1518  real, dimension(SZI_(G)) :: I_decay ! A decay rate [Z-1 ~> m-1].
1519  real, dimension(SZI_(G)) :: Kd_mlr_ml ! Diffusivities associated with mixed layer radiation [Z2 T-1 ~> m2 s-1].
1520 
1521  real :: f_sq ! The square of the local Coriolis parameter or a related variable [T-2 ~> s-2].
1522  real :: h_ml_sq ! The square of the mixed layer thickness [Z2 ~> m2].
1523  real :: ustar_sq ! ustar squared [Z2 T-2 ~> m2 s-2]
1524  real :: Kd_mlr ! A diffusivity associated with mixed layer turbulence radiation [Z2 T-1 ~> m2 s-1].
1525  real :: C1_6 ! 1/6
1526  real :: Omega2 ! rotation rate squared [T-2 ~> s-2].
1527  real :: z1 ! layer thickness times I_decay [nondim]
1528  real :: dzL ! thickness converted to heights [Z ~> m].
1529  real :: I_decay_len2_TKE ! squared inverse decay lengthscale for
1530  ! TKE, as used in the mixed layer code [Z-2 ~> m-2].
1531  real :: h_neglect ! negligibly small thickness [Z ~> m].
1532 
1533  logical :: do_any, do_i(SZI_(G))
1534  integer :: i, k, is, ie, nz, kml
1535  is = g%isc ; ie = g%iec ; nz = g%ke
1536 
1537  omega2 = cs%omega**2
1538  c1_6 = 1.0 / 6.0
1539  kml = gv%nkml
1540  h_neglect = gv%H_subroundoff*gv%H_to_Z
1541 
1542  if (.not.cs%ML_radiation) return
1543 
1544  do i=is,ie ; h_ml(i) = 0.0 ; do_i(i) = (g%mask2dT(i,j) > 0.5) ; enddo
1545  do k=1,kml ; do i=is,ie ; h_ml(i) = h_ml(i) + gv%H_to_Z*h(i,j,k) ; enddo ; enddo
1546 
1547  do i=is,ie ; if (do_i(i)) then
1548  if (cs%ML_omega_frac >= 1.0) then
1549  f_sq = 4.0 * omega2
1550  else
1551  f_sq = 0.25 * ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
1552  (g%CoriolisBu(i,j-1)**2 + g%CoriolisBu(i-1,j)**2))
1553  if (cs%ML_omega_frac > 0.0) &
1554  f_sq = cs%ML_omega_frac * 4.0 * omega2 + (1.0 - cs%ML_omega_frac) * f_sq
1555  endif
1556 
1557  ustar_sq = max(fluxes%ustar(i,j), cs%ustar_min)**2
1558 
1559  tke_ml_flux(i) = (cs%mstar * cs%ML_rad_coeff) * (ustar_sq * (fluxes%ustar(i,j)))
1560  i_decay_len2_tke = cs%TKE_decay**2 * (f_sq / ustar_sq)
1561 
1562  if (cs%ML_rad_TKE_decay) &
1563  tke_ml_flux(i) = tke_ml_flux(i) * exp(-h_ml(i) * sqrt(i_decay_len2_tke))
1564 
1565  ! Calculate the inverse decay scale
1566  h_ml_sq = (cs%ML_rad_efold_coeff * (h_ml(i)+h_neglect))**2
1567  i_decay(i) = sqrt((i_decay_len2_tke * h_ml_sq + 1.0) / h_ml_sq)
1568 
1569  ! Average the dissipation layer kml+1, using
1570  ! a more accurate Taylor series approximations for very thin layers.
1571  z1 = (gv%H_to_Z*h(i,j,kml+1)) * i_decay(i)
1572  if (z1 > 1e-5) then
1573  kd_mlr = tke_ml_flux(i) * tke_to_kd(i,kml+1) * (1.0 - exp(-z1))
1574  else
1575  kd_mlr = tke_ml_flux(i) * tke_to_kd(i,kml+1) * (z1 * (1.0 - z1 * (0.5 - c1_6 * z1)))
1576  endif
1577  kd_mlr_ml(i) = min(kd_mlr, cs%ML_rad_kd_max)
1578  tke_ml_flux(i) = tke_ml_flux(i) * exp(-z1)
1579  endif ; enddo
1580 
1581  do k=1,kml+1 ; do i=is,ie ; if (do_i(i)) then
1582  kd_lay(i,j,k) = kd_lay(i,j,k) + kd_mlr_ml(i)
1583  endif ; enddo ; enddo
1584  if (present(kd_int)) then
1585  do k=2,kml+1 ; do i=is,ie ; if (do_i(i)) then
1586  kd_int(i,j,k) = kd_int(i,j,k) + kd_mlr_ml(i)
1587  endif ; enddo ; enddo
1588  if (kml<=nz-1) then ; do i=is,ie ; if (do_i(i)) then
1589  kd_int(i,j,kml+2) = kd_int(i,j,kml+2) + 0.5 * kd_mlr_ml(i)
1590  endif ; enddo ; endif
1591  endif
1592 
1593  do k=kml+2,nz-1
1594  do_any = .false.
1595  do i=is,ie ; if (do_i(i)) then
1596  dzl = gv%H_to_Z*h(i,j,k) ; z1 = dzl*i_decay(i)
1597  if (cs%ML_Rad_bug) then
1598  ! These expresssions are dimensionally inconsistent. -RWH
1599  ! This is supposed to be the integrated energy deposited in the layer,
1600  ! not the average over the layer as in these expressions.
1601  if (z1 > 1e-5) then
1602  kd_mlr = (tke_ml_flux(i) * tke_to_kd(i,k)) * & ! Units of Z2 T-1
1603  us%m_to_Z * ((1.0 - exp(-z1)) / dzl) ! Units of m-1
1604  else
1605  kd_mlr = (tke_ml_flux(i) * tke_to_kd(i,k)) * & ! Units of Z2 T-1
1606  us%m_to_Z * (i_decay(i) * (1.0 - z1 * (0.5 - c1_6*z1))) ! Units of m-1
1607  endif
1608  else
1609  if (z1 > 1e-5) then
1610  kd_mlr = (tke_ml_flux(i) * tke_to_kd(i,k)) * (1.0 - exp(-z1))
1611  else
1612  kd_mlr = (tke_ml_flux(i) * tke_to_kd(i,k)) * (z1 * (1.0 - z1 * (0.5 - c1_6*z1)))
1613  endif
1614  endif
1615  kd_mlr = min(kd_mlr, cs%ML_rad_kd_max)
1616  kd_lay(i,j,k) = kd_lay(i,j,k) + kd_mlr
1617  if (present(kd_int)) then
1618  kd_int(i,j,k) = kd_int(i,j,k) + 0.5 * kd_mlr
1619  kd_int(i,j,k+1) = kd_int(i,j,k+1) + 0.5 * kd_mlr
1620  endif
1621 
1622  tke_ml_flux(i) = tke_ml_flux(i) * exp(-z1)
1623  if (tke_ml_flux(i) * i_decay(i) < 0.1 * cs%Kd_min * omega2) then
1624  do_i(i) = .false.
1625  else ; do_any = .true. ; endif
1626  endif ; enddo
1627  if (.not.do_any) exit
1628  enddo
1629 

Referenced by set_diffusivity().

Here is the caller graph for this function:

◆ double_diffusion()

subroutine mom_set_diffusivity::double_diffusion ( type(thermo_var_ptrs), intent(in)  tv,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  T_f,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  S_f,
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(set_diffusivity_cs), pointer  CS,
real, dimension( g %isd: g %ied, g %ke+1), intent(out)  Kd_T_dd,
real, dimension( g %isd: g %ied, g %ke+1), intent(out)  Kd_S_dd 
)
private

This subroutine sets the additional diffusivities of temperature and salinity due to double diffusion, using the same functional form as is used in MOM4.1, and taken from an NCAR technical note (REF?) that updates what was in Large et al. (1994). All the coefficients here should probably be made run-time variables rather than hard-coded constants.

Todo:
Find reference for NCAR tech note above.
Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]tvStructure containing pointers to any available thermodynamic fields; absent fields have NULL ptrs.
[in]hLayer thicknesses [H ~> m or kg m-2].
[in]t_flayer temperatures with the values in massless layers
[in]s_fLayer salinities with values in massless
[in]jMeridional index upon which to work.
csModule control structure.
[out]kd_t_ddInterface double diffusion diapycnal
[out]kd_s_ddInterface double diffusion diapycnal

Definition at line 1015 of file MOM_set_diffusivity.F90.

1015  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1016  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1017  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1018  type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
1019  !! thermodynamic fields; absent fields have NULL ptrs.
1020  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1021  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
1022  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1023  intent(in) :: T_f !< layer temperatures with the values in massless layers
1024  !! filled vertically by diffusion [degC].
1025  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1026  intent(in) :: S_f !< Layer salinities with values in massless
1027  !! layers filled vertically by diffusion [ppt].
1028  integer, intent(in) :: j !< Meridional index upon which to work.
1029  type(set_diffusivity_CS), pointer :: CS !< Module control structure.
1030  real, dimension(SZI_(G),SZK_(G)+1), &
1031  intent(out) :: Kd_T_dd !< Interface double diffusion diapycnal
1032  !! diffusivity for temp [Z2 T-1 ~> m2 s-1].
1033  real, dimension(SZI_(G),SZK_(G)+1), &
1034  intent(out) :: Kd_S_dd !< Interface double diffusion diapycnal
1035  !! diffusivity for saln [Z2 T-1 ~> m2 s-1].
1036 
1037  real, dimension(SZI_(G)) :: &
1038  dRho_dT, & ! partial derivatives of density wrt temp [R degC-1 ~> kg m-3 degC-1]
1039  dRho_dS, & ! partial derivatives of density wrt saln [R ppt-1 ~> kg m-3 ppt-1]
1040  pres, & ! pressure at each interface [Pa]
1041  Temp_int, & ! temperature at interfaces [degC]
1042  Salin_int ! Salinity at interfaces [ppt]
1043 
1044  real :: alpha_dT ! density difference between layers due to temp diffs [R ~> kg m-3]
1045  real :: beta_dS ! density difference between layers due to saln diffs [R ~> kg m-3]
1046 
1047  real :: Rrho ! vertical density ratio [nondim]
1048  real :: diff_dd ! factor for double-diffusion [nondim]
1049  real :: Kd_dd ! The dominant double diffusive diffusivity [Z2 T-1 ~> m2 s-1]
1050  real :: prandtl ! flux ratio for diffusive convection regime
1051 
1052  real, parameter :: Rrho0 = 1.9 ! limit for double-diffusive density ratio [nondim]
1053 
1054  integer :: i, k, is, ie, nz
1055  is = g%isc ; ie = g%iec ; nz = g%ke
1056 
1057  if (associated(tv%eqn_of_state)) then
1058  do i=is,ie
1059  pres(i) = 0.0 ; kd_t_dd(i,1) = 0.0 ; kd_s_dd(i,1) = 0.0
1060  kd_t_dd(i,nz+1) = 0.0 ; kd_s_dd(i,nz+1) = 0.0
1061  enddo
1062  do k=2,nz
1063  do i=is,ie
1064  pres(i) = pres(i) + gv%H_to_Pa*h(i,j,k-1)
1065  temp_int(i) = 0.5 * (t_f(i,j,k-1) + t_f(i,j,k))
1066  salin_int(i) = 0.5 * (s_f(i,j,k-1) + s_f(i,j,k))
1067  enddo
1068  call calculate_density_derivs(temp_int, salin_int, pres, &
1069  drho_dt(:), drho_ds(:), is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
1070 
1071  do i=is,ie
1072  alpha_dt = -1.0*drho_dt(i) * (t_f(i,j,k-1) - t_f(i,j,k))
1073  beta_ds = drho_ds(i) * (s_f(i,j,k-1) - s_f(i,j,k))
1074 
1075  if ((alpha_dt > beta_ds) .and. (beta_ds > 0.0)) then ! salt finger case
1076  rrho = min(alpha_dt / beta_ds, rrho0)
1077  diff_dd = 1.0 - ((rrho-1.0)/(rrho0-1.0))
1078  kd_dd = cs%Max_salt_diff_salt_fingers * diff_dd*diff_dd*diff_dd
1079  kd_t_dd(i,k) = 0.7 * kd_dd
1080  kd_s_dd(i,k) = kd_dd
1081  elseif ((alpha_dt < 0.) .and. (beta_ds < 0.) .and. (alpha_dt > beta_ds)) then ! diffusive convection
1082  rrho = alpha_dt / beta_ds
1083  kd_dd = cs%Kv_molecular * 0.909 * exp(4.6 * exp(-0.54 * (1/rrho - 1)))
1084  prandtl = 0.15*rrho
1085  if (rrho > 0.5) prandtl = (1.85-0.85/rrho)*rrho
1086  kd_t_dd(i,k) = kd_dd
1087  kd_s_dd(i,k) = prandtl * kd_dd
1088  else
1089  kd_t_dd(i,k) = 0.0 ; kd_s_dd(i,k) = 0.0
1090  endif
1091  enddo
1092  enddo
1093  endif
1094 

Referenced by set_diffusivity().

Here is the caller graph for this function:

◆ find_n2()

subroutine mom_set_diffusivity::find_n2 ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  T_f,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  S_f,
type(forcing), intent(in)  fluxes,
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(set_diffusivity_cs), pointer  CS,
real, dimension( g %isd: g %ied, g %ke+1), intent(out)  dRho_int,
real, dimension( g %isd: g %ied, g %ke), intent(out)  N2_lay,
real, dimension( g %isd: g %ied, g %ke+1), intent(out)  N2_int,
real, dimension( g %isd: g %ied), intent(out)  N2_bot 
)
private

Calculate Brunt-Vaisala frequency, N^2.

Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]usA dimensional unit scaling type
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]tvStructure containing pointers to any available thermodynamic fields.
[in]t_flayer temperature with the values in massless layers
[in]s_fLayer salinities with values in massless
[in]fluxesA structure of thermodynamic surface fluxes
[in]jj-index of row to work on
csDiffusivity control structure
[out]drho_intChange in locally referenced potential density
[out]n2_intThe squared buoyancy frequency at the interfaces [T-2 ~> s-2].
[out]n2_layThe squared buoyancy frequency of the layers [T-2 ~> s-2].
[out]n2_botThe near-bottom squared buoyancy frequency [T-2 ~> s-2].

Definition at line 840 of file MOM_set_diffusivity.F90.

840  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
841  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
842  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
843  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
844  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
845  type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
846  !! thermodynamic fields.
847  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
848  intent(in) :: T_f !< layer temperature with the values in massless layers
849  !! filled vertically by diffusion [degC].
850  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
851  intent(in) :: S_f !< Layer salinities with values in massless
852  !! layers filled vertically by diffusion [ppt].
853  type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes
854  integer, intent(in) :: j !< j-index of row to work on
855  type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure
856  real, dimension(SZI_(G),SZK_(G)+1), &
857  intent(out) :: dRho_int !< Change in locally referenced potential density
858  !! across each interface [R ~> kg m-3].
859  real, dimension(SZI_(G),SZK_(G)+1), &
860  intent(out) :: N2_int !< The squared buoyancy frequency at the interfaces [T-2 ~> s-2].
861  real, dimension(SZI_(G),SZK_(G)), &
862  intent(out) :: N2_lay !< The squared buoyancy frequency of the layers [T-2 ~> s-2].
863  real, dimension(SZI_(G)), intent(out) :: N2_bot !< The near-bottom squared buoyancy frequency [T-2 ~> s-2].
864  ! Local variables
865  real, dimension(SZI_(G),SZK_(G)+1) :: &
866  dRho_int_unfilt, & ! unfiltered density differences across interfaces [R ~> kg m-3]
867  dRho_dT, & ! partial derivative of density wrt temp [R degC-1 ~> kg m-3 degC-1]
868  dRho_dS ! partial derivative of density wrt saln [R ppt-1 ~> kg m-3 ppt-1]
869 
870  real, dimension(SZI_(G)) :: &
871  pres, & ! pressure at each interface [Pa]
872  Temp_int, & ! temperature at each interface [degC]
873  Salin_int, & ! salinity at each interface [ppt]
874  drho_bot, & ! A density difference [R ~> kg m-3]
875  h_amp, & ! The topographic roughness amplitude [Z ~> m].
876  hb, & ! The thickness of the bottom layer [Z ~> m].
877  z_from_bot ! The hieght above the bottom [Z ~> m].
878 
879  real :: Rml_base ! density of the deepest variable density layer
880  real :: dz_int ! thickness associated with an interface [Z ~> m].
881  real :: G_Rho0 ! gravitation acceleration divided by Bouss reference density
882  ! times some unit conversion factors [Z T-2 R-1 ~> m4 s-2 kg-1].
883  real :: H_neglect ! negligibly small thickness, in the same units as h.
884 
885  logical :: do_i(SZI_(G)), do_any
886  integer :: i, k, is, ie, nz
887 
888  is = g%isc ; ie = g%iec ; nz = g%ke
889  g_rho0 = (us%L_to_Z**2 * gv%g_Earth) / (gv%Rho0)
890  h_neglect = gv%H_subroundoff
891 
892  ! Find the (limited) density jump across each interface.
893  do i=is,ie
894  drho_int(i,1) = 0.0 ; drho_int(i,nz+1) = 0.0
895  drho_int_unfilt(i,1) = 0.0 ; drho_int_unfilt(i,nz+1) = 0.0
896  enddo
897  if (associated(tv%eqn_of_state)) then
898  if (associated(fluxes%p_surf)) then
899  do i=is,ie ; pres(i) = fluxes%p_surf(i,j) ; enddo
900  else
901  do i=is,ie ; pres(i) = 0.0 ; enddo
902  endif
903  do k=2,nz
904  do i=is,ie
905  pres(i) = pres(i) + gv%H_to_Pa*h(i,j,k-1)
906  temp_int(i) = 0.5 * (t_f(i,j,k) + t_f(i,j,k-1))
907  salin_int(i) = 0.5 * (s_f(i,j,k) + s_f(i,j,k-1))
908  enddo
909  call calculate_density_derivs(temp_int, salin_int, pres, &
910  drho_dt(:,k), drho_ds(:,k), is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
911  do i=is,ie
912  drho_int(i,k) = max(drho_dt(i,k)*(t_f(i,j,k) - t_f(i,j,k-1)) + &
913  drho_ds(i,k)*(s_f(i,j,k) - s_f(i,j,k-1)), 0.0)
914  drho_int_unfilt(i,k) = max(drho_dt(i,k)*(tv%T(i,j,k) - tv%T(i,j,k-1)) + &
915  drho_ds(i,k)*(tv%S(i,j,k) - tv%S(i,j,k-1)), 0.0)
916  enddo
917  enddo
918  else
919  do k=2,nz ; do i=is,ie
920  drho_int(i,k) = gv%Rlay(k) - gv%Rlay(k-1)
921  enddo ; enddo
922  endif
923 
924  ! Set the buoyancy frequencies.
925  do k=1,nz ; do i=is,ie
926  n2_lay(i,k) = g_rho0 * 0.5*(drho_int(i,k) + drho_int(i,k+1)) / &
927  (gv%H_to_Z*(h(i,j,k) + h_neglect))
928  enddo ; enddo
929  do i=is,ie ; n2_int(i,1) = 0.0 ; n2_int(i,nz+1) = 0.0 ; enddo
930  do k=2,nz ; do i=is,ie
931  n2_int(i,k) = g_rho0 * drho_int(i,k) / &
932  (0.5*gv%H_to_Z*(h(i,j,k-1) + h(i,j,k) + h_neglect))
933  enddo ; enddo
934 
935  ! Find the bottom boundary layer stratification, and use this in the deepest layers.
936  do i=is,ie
937  hb(i) = 0.0 ; drho_bot(i) = 0.0
938  z_from_bot(i) = 0.5*gv%H_to_Z*h(i,j,nz)
939  do_i(i) = (g%mask2dT(i,j) > 0.5)
940 
941  if ( (cs%tm_csp%Int_tide_dissipation .or. cs%tm_csp%Lee_wave_dissipation) .and. &
942  .not. cs%tm_csp%use_CVMix_tidal ) then
943  h_amp(i) = sqrt(cs%tm_csp%h2(i,j)) ! for computing Nb
944  else
945  h_amp(i) = 0.0
946  endif
947  enddo
948 
949  do k=nz,2,-1
950  do_any = .false.
951  do i=is,ie ; if (do_i(i)) then
952  dz_int = 0.5*gv%H_to_Z*(h(i,j,k) + h(i,j,k-1))
953  z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above
954 
955  hb(i) = hb(i) + dz_int
956  drho_bot(i) = drho_bot(i) + drho_int(i,k)
957 
958  if (z_from_bot(i) > h_amp(i)) then
959  if (k>2) then
960  ! Always include at least one full layer.
961  hb(i) = hb(i) + 0.5*gv%H_to_Z*(h(i,j,k-1) + h(i,j,k-2))
962  drho_bot(i) = drho_bot(i) + drho_int(i,k-1)
963  endif
964  do_i(i) = .false.
965  else
966  do_any = .true.
967  endif
968  endif ; enddo
969  if (.not.do_any) exit
970  enddo
971 
972  do i=is,ie
973  if (hb(i) > 0.0) then
974  n2_bot(i) = (g_rho0 * drho_bot(i)) / hb(i)
975  else ; n2_bot(i) = 0.0 ; endif
976  z_from_bot(i) = 0.5*gv%H_to_Z*h(i,j,nz)
977  do_i(i) = (g%mask2dT(i,j) > 0.5)
978  enddo
979 
980  do k=nz,2,-1
981  do_any = .false.
982  do i=is,ie ; if (do_i(i)) then
983  dz_int = 0.5*gv%H_to_Z*(h(i,j,k) + h(i,j,k-1))
984  z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above
985 
986  n2_int(i,k) = n2_bot(i)
987  if (k>2) n2_lay(i,k-1) = n2_bot(i)
988 
989  if (z_from_bot(i) > h_amp(i)) then
990  if (k>2) n2_int(i,k-1) = n2_bot(i)
991  do_i(i) = .false.
992  else
993  do_any = .true.
994  endif
995  endif ; enddo
996  if (.not.do_any) exit
997  enddo
998 
999  if (associated(tv%eqn_of_state)) then
1000  do k=1,nz+1 ; do i=is,ie
1001  drho_int(i,k) = drho_int_unfilt(i,k)
1002  enddo ; enddo
1003  endif
1004 

Referenced by set_diffusivity().

Here is the caller graph for this function:

◆ find_tke_to_kd()

subroutine mom_set_diffusivity::find_tke_to_kd ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
real, dimension( g %isd: g %ied, g %ke+1), intent(in)  dRho_int,
real, dimension( g %isd: g %ied, g %ke), intent(in)  N2_lay,
integer, intent(in)  j,
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(set_diffusivity_cs), pointer  CS,
real, dimension( g %isd: g %ied, g %ke), intent(out)  TKE_to_Kd,
real, dimension( g %isd: g %ied, g %ke), intent(out)  maxTKE,
integer, dimension( g %isd: g %ied), intent(out)  kb 
)
private

Convert turbulent kinetic energy to diffusivity.

Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]usA dimensional unit scaling type
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]tvStructure containing pointers to any available thermodynamic fields.
[in]drho_intChange in locally referenced potential density across each interface [R ~> kg m-3].
[in]n2_layThe squared buoyancy frequency of the layers [T-2 ~> s-2].
[in]jj-index of row to work on
[in]dtTime increment [T ~> s].
csDiffusivity control structure
[out]tke_to_kdThe conversion rate between the TKE dissipated within a layer and the diapycnal diffusivity witin that layer, usually (~Rho_0 / (G_Earth * dRho_lay)) [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1]
[out]maxtkeThe energy required to for a layer to entrain to its maximum realizable thickness [Z3 T-3 ~> m3 s-3]
[out]kbIndex of lightest layer denser than the buffer layer, or -1 without a bulk mixed layer.

Definition at line 625 of file MOM_set_diffusivity.F90.

625  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
626  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
627  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
628  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
629  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
630  type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
631  !! thermodynamic fields.
632  real, dimension(SZI_(G),SZK_(G)+1), intent(in) :: dRho_int !< Change in locally referenced potential density
633  !! across each interface [R ~> kg m-3].
634  real, dimension(SZI_(G),SZK_(G)), intent(in) :: N2_lay !< The squared buoyancy frequency of the
635  !! layers [T-2 ~> s-2].
636  integer, intent(in) :: j !< j-index of row to work on
637  real, intent(in) :: dt !< Time increment [T ~> s].
638  type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure
639  real, dimension(SZI_(G),SZK_(G)), intent(out) :: TKE_to_Kd !< The conversion rate between the
640  !! TKE dissipated within a layer and the
641  !! diapycnal diffusivity witin that layer,
642  !! usually (~Rho_0 / (G_Earth * dRho_lay))
643  !! [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1]
644  real, dimension(SZI_(G),SZK_(G)), intent(out) :: maxTKE !< The energy required to for a layer to entrain
645  !! to its maximum realizable thickness [Z3 T-3 ~> m3 s-3]
646  integer, dimension(SZI_(G)), intent(out) :: kb !< Index of lightest layer denser than the buffer
647  !! layer, or -1 without a bulk mixed layer.
648  ! Local variables
649  real, dimension(SZI_(G),SZK_(G)) :: &
650  ds_dsp1, & ! coordinate variable (sigma-2) difference across an
651  ! interface divided by the difference across the interface
652  ! below it [nondim]
653  dsp1_ds, & ! inverse coordinate variable (sigma-2) difference
654  ! across an interface times the difference across the
655  ! interface above it [nondim]
656  rho_0, & ! Layer potential densities relative to surface pressure [R ~> kg m-3]
657  maxent ! maxEnt is the maximum value of entrainment from below (with
658  ! compensating entrainment from above to keep the layer
659  ! density from changing) that will not deplete all of the
660  ! layers above or below a layer within a timestep [Z ~> m].
661  real, dimension(SZI_(G)) :: &
662  htot, & ! total thickness above or below a layer, or the
663  ! integrated thickness in the BBL [Z ~> m].
664  mfkb, & ! total thickness in the mixed and buffer layers
665  ! times ds_dsp1 [Z ~> m].
666  p_ref, & ! array of tv%P_Ref pressures
667  rcv_kmb, & ! coordinate density in the lowest buffer layer [R ~> kg m-3]
668  p_0 ! An array of 0 pressures
669 
670  real :: dh_max ! maximum amount of entrainment a layer could
671  ! undergo before entraining all fluid in the layers
672  ! above or below [Z ~> m].
673  real :: dRho_lay ! density change across a layer [R ~> kg m-3]
674  real :: Omega2 ! rotation rate squared [T-2 ~> s-2]
675  real :: G_Rho0 ! gravitation accel divided by Bouss ref density [Z T-2 R-1 ~> m4 s-2 kg-1]
676  real :: G_IRho0 ! Alternate calculation of G_Rho0 for reproducibility [Z T-2 R-1 ~> m4 s-2 kg-1]
677  real :: I_Rho0 ! inverse of Boussinesq reference density [R-1 ~> m3 kg-1]
678  real :: I_dt ! 1/dt [T-1 ~> s-1]
679  real :: H_neglect ! negligibly small thickness [H ~> m or kg m-2]
680  real :: hN2pO2 ! h (N^2 + Omega^2), in [m3 T-2 Z-2 ~> m s-2].
681  logical :: do_i(SZI_(G))
682 
683  integer :: i, k, is, ie, nz, i_rem, kmb, kb_min
684  is = g%isc ; ie = g%iec ; nz = g%ke
685 
686  i_dt = 1.0 / dt
687  omega2 = cs%omega**2
688  h_neglect = gv%H_subroundoff
689  g_rho0 = (us%L_to_Z**2 * gv%g_Earth) / (gv%Rho0)
690  if (cs%answers_2018) then
691  i_rho0 = 1.0 / (gv%Rho0)
692  g_irho0 = (us%L_to_Z**2 * gv%g_Earth) * i_rho0
693  else
694  g_irho0 = g_rho0
695  endif
696 
697  ! Simple but coordinate-independent estimate of Kd/TKE
698  if (cs%simple_TKE_to_Kd) then
699  do k=1,nz ; do i=is,ie
700  hn2po2 = (gv%H_to_Z * h(i,j,k)) * (n2_lay(i,k) + omega2) ! Units of Z T-2.
701  if (hn2po2>0.) then
702  tke_to_kd(i,k) = 1.0 / hn2po2 ! Units of T2 Z-1.
703  else; tke_to_kd(i,k) = 0.; endif
704  ! The maximum TKE conversion we allow is really a statement
705  ! about the upper diffusivity we allow. Kd_max must be set.
706  maxtke(i,k) = hn2po2 * cs%Kd_max ! Units of Z3 T-3.
707  enddo ; enddo
708  kb(is:ie) = -1 ! kb should not be used by any code in non-layered mode -AJA
709  return
710  endif
711 
712  ! Determine kb - the index of the shallowest active interior layer.
713  if (cs%bulkmixedlayer) then
714  kmb = gv%nk_rho_varies
715  do i=is,ie ; p_0(i) = 0.0 ; p_ref(i) = tv%P_Ref ; enddo
716  do k=1,nz
717  call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_0, rho_0(:,k), &
718  is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
719  enddo
720  call calculate_density(tv%T(:,j,kmb), tv%S(:,j,kmb), p_ref, rcv_kmb, &
721  is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
722 
723  kb_min = kmb+1
724  do i=is,ie
725  ! Determine the next denser layer than the buffer layer in the
726  ! coordinate density (sigma-2).
727  do k=kmb+1,nz-1 ; if (rcv_kmb(i) <= gv%Rlay(k)) exit ; enddo
728  kb(i) = k
729 
730  ! Backtrack, in case there are massive layers above that are stable
731  ! in sigma-0.
732  do k=kb(i)-1,kmb+1,-1
733  if (rho_0(i,kmb) > rho_0(i,k)) exit
734  if (h(i,j,k)>2.0*gv%Angstrom_H) kb(i) = k
735  enddo
736  enddo
737 
738  call set_density_ratios(h, tv, kb, g, gv, us, cs, j, ds_dsp1, rho_0)
739  else ! not bulkmixedlayer
740  kb_min = 2 ; kmb = 0
741  do i=is,ie ; kb(i) = 1 ; enddo
742  call set_density_ratios(h, tv, kb, g, gv, us, cs, j, ds_dsp1)
743  endif
744 
745  ! Determine maxEnt - the maximum permitted entrainment from below by each
746  ! interior layer.
747  do k=2,nz-1 ; do i=is,ie
748  dsp1_ds(i,k) = 1.0 / ds_dsp1(i,k)
749  enddo ; enddo
750  do i=is,ie ; dsp1_ds(i,nz) = 0.0 ; enddo
751 
752  if (cs%bulkmixedlayer) then
753  kmb = gv%nk_rho_varies
754  do i=is,ie
755  htot(i) = gv%H_to_Z*h(i,j,kmb)
756  mfkb(i) = 0.0
757  if (kb(i) < nz) &
758  mfkb(i) = ds_dsp1(i,kb(i)) * (gv%H_to_Z*(h(i,j,kmb) - gv%Angstrom_H))
759  enddo
760  do k=1,kmb-1 ; do i=is,ie
761  htot(i) = htot(i) + gv%H_to_Z*h(i,j,k)
762  mfkb(i) = mfkb(i) + ds_dsp1(i,k+1)*(gv%H_to_Z*(h(i,j,k) - gv%Angstrom_H))
763  enddo ; enddo
764  else
765  do i=is,i
766  maxent(i,1) = 0.0 ; htot(i) = gv%H_to_Z*(h(i,j,1) - gv%Angstrom_H)
767  enddo
768  endif
769  do k=kb_min,nz-1 ; do i=is,ie
770  if (k == kb(i)) then
771  maxent(i,kb(i)) = mfkb(i)
772  elseif (k > kb(i)) then
773  if (cs%answers_2018) then
774  maxent(i,k) = (1.0/dsp1_ds(i,k))*(maxent(i,k-1) + htot(i))
775  else
776  maxent(i,k) = ds_dsp1(i,k)*(maxent(i,k-1) + htot(i))
777  endif
778  htot(i) = htot(i) + gv%H_to_Z*(h(i,j,k) - gv%Angstrom_H)
779  endif
780  enddo ; enddo
781 
782  do i=is,ie
783  htot(i) = gv%H_to_Z*(h(i,j,nz) - gv%Angstrom_H) ; maxent(i,nz) = 0.0
784  do_i(i) = (g%mask2dT(i,j) > 0.5)
785  enddo
786  do k=nz-1,kb_min,-1
787  i_rem = 0
788  do i=is,ie ; if (do_i(i)) then
789  if (k<kb(i)) then ; do_i(i) = .false. ; cycle ; endif
790  i_rem = i_rem + 1 ! Count the i-rows that are still being worked on.
791  maxent(i,k) = min(maxent(i,k),dsp1_ds(i,k+1)*maxent(i,k+1) + htot(i))
792  htot(i) = htot(i) + gv%H_to_Z*(h(i,j,k) - gv%Angstrom_H)
793  endif ; enddo
794  if (i_rem == 0) exit
795  enddo ! k-loop
796 
797  ! Now set maxTKE and TKE_to_Kd.
798  do i=is,ie
799  maxtke(i,1) = 0.0 ; tke_to_kd(i,1) = 0.0
800  maxtke(i,nz) = 0.0 ; tke_to_kd(i,nz) = 0.0
801  enddo
802  do k=2,kmb ; do i=is,ie
803  maxtke(i,k) = 0.0
804  tke_to_kd(i,k) = 1.0 / ((n2_lay(i,k) + omega2) * &
805  (gv%H_to_Z*(h(i,j,k) + h_neglect)))
806  enddo ; enddo
807  do k=kmb+1,kb_min-1 ; do i=is,ie
808  ! These are the properties in the deeper mixed and buffer layers, and
809  ! should perhaps be revisited.
810  maxtke(i,k) = 0.0 ; tke_to_kd(i,k) = 0.0
811  enddo ; enddo
812  do k=kb_min,nz-1 ; do i=is,ie
813  if (k<kb(i)) then
814  maxtke(i,k) = 0.0
815  tke_to_kd(i,k) = 0.0
816  else
817  ! maxTKE is found by determining the kappa that gives maxEnt.
818  ! kappa_max = I_dt * dRho_int(i,K+1) * maxEnt(i,k) * &
819  ! (GV%H_to_Z*h(i,j,k) + dh_max) / dRho_lay
820  ! maxTKE(i,k) = (GV%g_Earth*US%L_to_Z**2) * dRho_lay * kappa_max
821  ! dRho_int should already be non-negative, so the max is redundant?
822  dh_max = maxent(i,k) * (1.0 + dsp1_ds(i,k))
823  drho_lay = 0.5 * max(drho_int(i,k) + drho_int(i,k+1), 0.0)
824  maxtke(i,k) = i_dt * (g_irho0 * &
825  (0.5*max(drho_int(i,k+1) + dsp1_ds(i,k)*drho_int(i,k), 0.0))) * &
826  ((gv%H_to_Z*h(i,j,k) + dh_max) * maxent(i,k))
827  ! TKE_to_Kd should be rho_InSitu / G_Earth * (delta rho_InSitu)
828  ! The omega^2 term in TKE_to_Kd is due to a rescaling of the efficiency of turbulent
829  ! mixing by a factor of N^2 / (N^2 + Omega^2), as proposed by Melet et al., 2013?
830  tke_to_kd(i,k) = 1.0 / (g_rho0 * drho_lay + &
831  cs%omega**2 * gv%H_to_Z*(h(i,j,k) + h_neglect))
832  endif
833  enddo ; enddo
834 

References set_density_ratios().

Referenced by set_diffusivity().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ set_bbl_tke()

subroutine, public mom_set_diffusivity::set_bbl_tke ( real, dimension(szib_(g),szj_(g),szk_(g)), intent(in)  u,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  v,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
type(forcing), intent(in)  fluxes,
type(vertvisc_type), intent(in)  visc,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(set_diffusivity_cs), pointer  CS 
)

This subroutine calculates several properties related to bottom boundary layer turbulence.

Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]usA dimensional unit scaling type
[in]uThe zonal velocity [L T-1 ~> m s-1]
[in]vThe meridional velocity [L T-1 ~> m s-1]
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]fluxesA structure of thermodynamic surface fluxes
[in]viscStructure containing vertical viscosities, bottom boundary layer properies, and related fields.
csDiffusivity control structure

Definition at line 1635 of file MOM_set_diffusivity.F90.

1635  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1636  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1637  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1638  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1639  intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1]
1640  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1641  intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]
1642  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1643  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1644  type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes
1645  type(vertvisc_type), intent(in) :: visc !< Structure containing vertical viscosities, bottom
1646  !! boundary layer properies, and related fields.
1647  type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure
1648 
1649  ! This subroutine calculates several properties related to bottom
1650  ! boundary layer turbulence.
1651 
1652  real, dimension(SZI_(G)) :: &
1653  htot ! total thickness above or below a layer, or the
1654  ! integrated thickness in the BBL [Z ~> m].
1655 
1656  real, dimension(SZIB_(G)) :: &
1657  uhtot, & ! running integral of u in the BBL [Z L T-1 ~> m2 s-1]
1658  ustar, & ! bottom boundary layer turbulence speed [Z T-1 ~> m s-1].
1659  u2_bbl ! square of the mean zonal velocity in the BBL [L2 T-2 ~> m2 s-2]
1660 
1661  real :: vhtot(SZI_(G)) ! running integral of v in the BBL [Z L T-1 ~> m2 s-1]
1662 
1663  real, dimension(SZI_(G),SZJB_(G)) :: &
1664  vstar, & ! ustar at at v-points [Z T-1 ~> m s-1].
1665  v2_bbl ! square of average meridional velocity in BBL [L2 T-2 ~> m2 s-2]
1666 
1667  real :: cdrag_sqrt ! square root of the drag coefficient [nondim]
1668  real :: hvel ! thickness at velocity points [Z ~> m].
1669 
1670  logical :: domore, do_i(SZI_(G))
1671  integer :: i, j, k, is, ie, js, je, nz
1672  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1673 
1674  if (.not.associated(cs)) call mom_error(fatal,"set_BBL_TKE: "//&
1675  "Module must be initialized before it is used.")
1676 
1677  if (.not.cs%bottomdraglaw .or. (cs%BBL_effic<=0.0)) then
1678  if (associated(visc%ustar_BBL)) then
1679  do j=js,je ; do i=is,ie ; visc%ustar_BBL(i,j) = 0.0 ; enddo ; enddo
1680  endif
1681  if (associated(visc%TKE_BBL)) then
1682  do j=js,je ; do i=is,ie ; visc%TKE_BBL(i,j) = 0.0 ; enddo ; enddo
1683  endif
1684  return
1685  endif
1686 
1687  cdrag_sqrt = sqrt(cs%cdrag)
1688 
1689 !$OMP parallel default(none) shared(cdrag_sqrt,is,ie,js,je,nz,visc,CS,G,GV,US,vstar,h,v, &
1690 !$OMP v2_bbl,u) &
1691 !$OMP private(do_i,vhtot,htot,domore,hvel,uhtot,ustar,u2_bbl)
1692 !$OMP do
1693  do j=js-1,je
1694  ! Determine ustar and the square magnitude of the velocity in the
1695  ! bottom boundary layer. Together these give the TKE source and
1696  ! vertical decay scale.
1697  do i=is,ie ; if ((g%mask2dCv(i,j) > 0.5) .and. (cdrag_sqrt*visc%bbl_thick_v(i,j) > 0.0)) then
1698  do_i(i) = .true. ; vhtot(i) = 0.0 ; htot(i) = 0.0
1699  vstar(i,j) = visc%Kv_bbl_v(i,j) / (cdrag_sqrt*visc%bbl_thick_v(i,j))
1700  else
1701  do_i(i) = .false. ; vstar(i,j) = 0.0 ; htot(i) = 0.0
1702  endif ; enddo
1703  do k=nz,1,-1
1704  domore = .false.
1705  do i=is,ie ; if (do_i(i)) then
1706  hvel = 0.5*gv%H_to_Z*(h(i,j,k) + h(i,j+1,k))
1707  if ((htot(i) + hvel) >= visc%bbl_thick_v(i,j)) then
1708  vhtot(i) = vhtot(i) + (visc%bbl_thick_v(i,j) - htot(i))*v(i,j,k)
1709  htot(i) = visc%bbl_thick_v(i,j)
1710  do_i(i) = .false.
1711  else
1712  vhtot(i) = vhtot(i) + hvel*v(i,j,k)
1713  htot(i) = htot(i) + hvel
1714  domore = .true.
1715  endif
1716  endif ; enddo
1717  if (.not.domore) exit
1718  enddo
1719  do i=is,ie ; if ((g%mask2dCv(i,j) > 0.5) .and. (htot(i) > 0.0)) then
1720  v2_bbl(i,j) = (vhtot(i)*vhtot(i))/(htot(i)*htot(i))
1721  else
1722  v2_bbl(i,j) = 0.0
1723  endif ; enddo
1724  enddo
1725 !$OMP do
1726  do j=js,je
1727  do i=is-1,ie ; if ((g%mask2dCu(i,j) > 0.5) .and. (cdrag_sqrt*visc%bbl_thick_u(i,j) > 0.0)) then
1728  do_i(i) = .true. ; uhtot(i) = 0.0 ; htot(i) = 0.0
1729  ustar(i) = visc%Kv_bbl_u(i,j) / (cdrag_sqrt*visc%bbl_thick_u(i,j))
1730  else
1731  do_i(i) = .false. ; ustar(i) = 0.0 ; htot(i) = 0.0
1732  endif ; enddo
1733  do k=nz,1,-1 ; domore = .false.
1734  do i=is-1,ie ; if (do_i(i)) then
1735  hvel = 0.5*gv%H_to_Z*(h(i,j,k) + h(i+1,j,k))
1736  if ((htot(i) + hvel) >= visc%bbl_thick_u(i,j)) then
1737  uhtot(i) = uhtot(i) + (visc%bbl_thick_u(i,j) - htot(i))*u(i,j,k)
1738  htot(i) = visc%bbl_thick_u(i,j)
1739  do_i(i) = .false.
1740  else
1741  uhtot(i) = uhtot(i) + hvel*u(i,j,k)
1742  htot(i) = htot(i) + hvel
1743  domore = .true.
1744  endif
1745  endif ; enddo
1746  if (.not.domore) exit
1747  enddo
1748  do i=is-1,ie ; if ((g%mask2dCu(i,j) > 0.5) .and. (htot(i) > 0.0)) then
1749  u2_bbl(i) = (uhtot(i)*uhtot(i))/(htot(i)*htot(i))
1750  else
1751  u2_bbl(i) = 0.0
1752  endif ; enddo
1753 
1754  do i=is,ie
1755  visc%ustar_BBL(i,j) = sqrt(0.5*g%IareaT(i,j) * &
1756  ((g%areaCu(i-1,j)*(ustar(i-1)*ustar(i-1)) + &
1757  g%areaCu(i,j)*(ustar(i)*ustar(i))) + &
1758  (g%areaCv(i,j-1)*(vstar(i,j-1)*vstar(i,j-1)) + &
1759  g%areaCv(i,j)*(vstar(i,j)*vstar(i,j))) ) )
1760  visc%TKE_BBL(i,j) = us%L_to_Z**2 * &
1761  (((g%areaCu(i-1,j)*(ustar(i-1)*u2_bbl(i-1)) + &
1762  g%areaCu(i,j) * (ustar(i)*u2_bbl(i))) + &
1763  (g%areaCv(i,j-1)*(vstar(i,j-1)*v2_bbl(i,j-1)) + &
1764  g%areaCv(i,j) * (vstar(i,j)*v2_bbl(i,j))) )*g%IareaT(i,j))
1765  enddo
1766  enddo
1767 !$OMP end parallel
1768 

◆ set_density_ratios()

subroutine mom_set_diffusivity::set_density_ratios ( real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
integer, dimension(szi_(g)), intent(in)  kb,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(set_diffusivity_cs), pointer  CS,
integer, intent(in)  j,
real, dimension(szi_(g),szk_(g)), intent(out)  ds_dsp1,
real, dimension(szi_(g),szk_(g)), intent(in), optional  rho_0 
)
private
Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]hLayer thicknesses [H ~> m or kg m-2].
[in]tvStructure containing pointers to any available thermodynamic fields; absent fields have NULL ptrs.
[in]kbIndex of lightest layer denser than the buffer layer, or -1 without a bulk mixed layer.
[in]usA dimensional unit scaling type
csControl structure returned by previous call to diabatic_entrain_init.
[in]jMeridional index upon which to work.
[out]ds_dsp1Coordinate variable (sigma-2) difference across an interface divided by the difference across the interface below it [nondim]
[in]rho_0Layer potential densities relative to

Definition at line 1772 of file MOM_set_diffusivity.F90.

1772  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1773  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1774  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1775  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
1776  type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any
1777  !! available thermodynamic fields; absent
1778  !! fields have NULL ptrs.
1779  integer, dimension(SZI_(G)), intent(in) :: kb !< Index of lightest layer denser than the buffer
1780  !! layer, or -1 without a bulk mixed layer.
1781  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1782  type(set_diffusivity_CS), pointer :: CS !< Control structure returned by previous
1783  !! call to diabatic_entrain_init.
1784  integer, intent(in) :: j !< Meridional index upon which to work.
1785  real, dimension(SZI_(G),SZK_(G)), intent(out) :: ds_dsp1 !< Coordinate variable (sigma-2)
1786  !! difference across an interface divided by
1787  !! the difference across the interface below
1788  !! it [nondim]
1789  real, dimension(SZI_(G),SZK_(G)), &
1790  optional, intent(in) :: rho_0 !< Layer potential densities relative to
1791  !! surface press [R ~> kg m-3].
1792 
1793  ! Local variables
1794  real :: g_R0 ! g_R0 is a rescaled version of g/Rho [L2 Z-1 R-1 T-2 ~> m4 kg-1 s-2]
1795  real :: eps, tmp ! nondimensional temproray variables
1796  real :: a(SZK_(G)), a_0(SZK_(G)) ! nondimensional temporary variables
1797  real :: p_ref(SZI_(G)) ! an array of tv%P_Ref pressures
1798  real :: Rcv(SZI_(G),SZK_(G)) ! coordinate density in the mixed and buffer layers [R ~> kg m-3]
1799  real :: I_Drho ! temporary variable [R-1 ~> m3 kg-1]
1800 
1801  integer :: i, k, k3, is, ie, nz, kmb
1802  is = g%isc ; ie = g%iec ; nz = g%ke
1803 
1804  do k=2,nz-1
1805  if (gv%g_prime(k+1) /= 0.0) then
1806  do i=is,ie
1807  ds_dsp1(i,k) = gv%g_prime(k) / gv%g_prime(k+1)
1808  enddo
1809  else
1810  do i=is,ie
1811  ds_dsp1(i,k) = 1.
1812  enddo
1813  endif
1814  enddo
1815 
1816  if (cs%bulkmixedlayer) then
1817  g_r0 = gv%g_Earth / (gv%Rho0)
1818  kmb = gv%nk_rho_varies
1819  eps = 0.1
1820  do i=is,ie ; p_ref(i) = tv%P_Ref ; enddo
1821  do k=1,kmb
1822  call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, rcv(:,k), &
1823  is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
1824  enddo
1825  do i=is,ie
1826  if (kb(i) <= nz-1) then
1827 ! Set up appropriately limited ratios of the reduced gravities of the
1828 ! interfaces above and below the buffer layer and the next denser layer.
1829  k = kb(i)
1830 
1831  i_drho = g_r0 / gv%g_prime(k+1)
1832  ! The indexing convention for a is appropriate for the interfaces.
1833  do k3=1,kmb
1834  a(k3+1) = (gv%Rlay(k) - rcv(i,k3)) * i_drho
1835  enddo
1836  if ((present(rho_0)) .and. (a(kmb+1) < 2.0*eps*ds_dsp1(i,k))) then
1837 ! If the buffer layer nearly matches the density of the layer below in the
1838 ! coordinate variable (sigma-2), use the sigma-0-based density ratio if it is
1839 ! greater (and stable).
1840  if ((rho_0(i,k) > rho_0(i,kmb)) .and. &
1841  (rho_0(i,k+1) > rho_0(i,k))) then
1842  i_drho = 1.0 / (rho_0(i,k+1)-rho_0(i,k))
1843  a_0(kmb+1) = min((rho_0(i,k)-rho_0(i,kmb)) * i_drho, ds_dsp1(i,k))
1844  if (a_0(kmb+1) > a(kmb+1)) then
1845  do k3=2,kmb
1846  a_0(k3) = a_0(kmb+1) + (rho_0(i,kmb)-rho_0(i,k3-1)) * i_drho
1847  enddo
1848  if (a(kmb+1) <= eps*ds_dsp1(i,k)) then
1849  do k3=2,kmb+1 ; a(k3) = a_0(k3) ; enddo
1850  else
1851 ! Alternative... tmp = 0.5*(1.0 - cos(PI*(a(K2+1)/(eps*ds_dsp1(i,k)) - 1.0)) )
1852  tmp = a(kmb+1)/(eps*ds_dsp1(i,k)) - 1.0
1853  do k3=2,kmb+1 ; a(k3) = tmp*a(k3) + (1.0-tmp)*a_0(k3) ; enddo
1854  endif
1855  endif
1856  endif
1857  endif
1858 
1859  ds_dsp1(i,k) = max(a(kmb+1),1e-5)
1860 
1861  do k3=2,kmb
1862 ! ds_dsp1(i,k3) = MAX(a(k3),1e-5)
1863  ! Deliberately treat convective instabilies of the upper mixed
1864  ! and buffer layers with respect to the deepest buffer layer as
1865  ! though they don't exist. They will be eliminated by the upcoming
1866  ! call to the mixedlayer code anyway.
1867  ! The indexing convention is appropriate for the interfaces.
1868  ds_dsp1(i,k3) = max(a(k3),ds_dsp1(i,k))
1869  enddo
1870  endif ! (kb(i) <= nz-1)
1871  enddo ! I-loop.
1872  endif ! bulkmixedlayer
1873 

Referenced by find_tke_to_kd().

Here is the caller graph for this function:

◆ set_diffusivity()

subroutine, public mom_set_diffusivity::set_diffusivity ( real, dimension(szib_(g),szj_(g),szk_(g)), intent(in)  u,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  v,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  u_h,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  v_h,
type(thermo_var_ptrs), intent(inout)  tv,
type(forcing), intent(in)  fluxes,
type(optics_type), pointer  optics,
type(vertvisc_type), intent(inout)  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(set_diffusivity_cs), pointer  CS,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(out)  Kd_lay,
real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(out), optional  Kd_int 
)

Sets the interior vertical diffusion of scalars due to the following processes:

  1. Shear-driven mixing: two options, Jackson et at. and KPP interior;
  2. Background mixing via CVMix (Bryan-Lewis profile) or the scheme described by Harrison & Hallberg, JPO 2008;
  3. Double-diffusion, old method and new method via CVMix;
  4. Tidal mixing: many options available, see MOM_tidal_mixing.F90; In addition, this subroutine has the option to set the interior vertical viscosity associated with processes 1,2 and 4 listed above, which is stored in viscKv_slow. Vertical viscosity due to shear-driven mixing is passed via viscKv_shear
    Parameters
    [in]gThe ocean's grid structure.
    [in]gvThe ocean's vertical grid structure.
    [in]usA dimensional unit scaling type
    [in]uThe zonal velocity [L T-1 ~> m s-1].
    [in]vThe meridional velocity [L T-1 ~> m s-1].
    [in]hLayer thicknesses [H ~> m or kg m-2].
    [in]u_hZonal velocity interpolated to h points [L T-1 ~> m s-1].
    [in]v_hMeridional velocity interpolated to h points [L T-1 ~> m s-1].
    [in,out]tvStructure with pointers to thermodynamic fields. Out is for tvTempxPmE.
    [in]fluxesA structure of thermodynamic surface fluxes
    opticsA structure describing the optical properties of the ocean.
    [in,out]viscStructure containing vertical viscosities, bottom boundary layer properies, and related fields.
    [in]dtTime increment [T ~> s].
    csModule control structure.
    [out]kd_layDiapycnal diffusivity of each layer [Z2 T-1 ~> m2 s-1].
    [out]kd_intDiapycnal diffusivity at each interface [Z2 T-1 ~> m2 s-1].

Definition at line 206 of file MOM_set_diffusivity.F90.

206  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
207  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
208  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
209  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
210  intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
211  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
212  intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
213  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
214  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
215  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
216  intent(in) :: u_h !< Zonal velocity interpolated to h points [L T-1 ~> m s-1].
217  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
218  intent(in) :: v_h !< Meridional velocity interpolated to h points [L T-1 ~> m s-1].
219  type(thermo_var_ptrs), intent(inout) :: tv !< Structure with pointers to thermodynamic
220  !! fields. Out is for tv%TempxPmE.
221  type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes
222  type(optics_type), pointer :: optics !< A structure describing the optical
223  !! properties of the ocean.
224  type(vertvisc_type), intent(inout) :: visc !< Structure containing vertical viscosities, bottom
225  !! boundary layer properies, and related fields.
226  real, intent(in) :: dt !< Time increment [T ~> s].
227  type(set_diffusivity_CS), pointer :: CS !< Module control structure.
228  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
229  intent(out) :: Kd_lay !< Diapycnal diffusivity of each layer [Z2 T-1 ~> m2 s-1].
230  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
231  optional, intent(out) :: Kd_int !< Diapycnal diffusivity at each interface [Z2 T-1 ~> m2 s-1].
232 
233  ! local variables
234  real, dimension(SZI_(G)) :: &
235  N2_bot ! bottom squared buoyancy frequency [T-2 ~> s-2]
236 
237  type(diffusivity_diags) :: dd ! structure w/ arrays of pointers to avail diags
238 
239  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
240  T_f, S_f ! Temperature and salinity [degC] and [ppt] with
241  ! massless layers filled vertically by diffusion.
242  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
243  T_adj, S_adj ! Temperature and salinity [degC] and [ppt]
244  ! after full convective adjustment.
245 
246  real, dimension(SZI_(G),SZK_(G)) :: &
247  N2_lay, & !< squared buoyancy frequency associated with layers [T-2 ~> s-2]
248  maxTKE, & !< energy required to entrain to h_max [Z3 T-3 ~> m3 s-3]
249  TKE_to_Kd !< conversion rate (~1.0 / (G_Earth + dRho_lay)) between
250  !< TKE dissipated within a layer and Kd in that layer
251  !< [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1]
252 
253  real, dimension(SZI_(G),SZK_(G)+1) :: &
254  N2_int, & !< squared buoyancy frequency associated at interfaces [T-2 ~> s-2]
255  dRho_int, & !< locally ref potential density difference across interfaces [R ~> kg m-3]
256  KT_extra, & !< double difusion diffusivity of temperature [Z2 T-1 ~> m2 s-1]
257  KS_extra !< double difusion diffusivity of salinity [Z2 T-1 ~> m2 s-1]
258 
259  real :: dissip ! local variable for dissipation calculations [Z2 R T-3 ~> W m-3]
260  real :: Omega2 ! squared absolute rotation rate [T-2 ~> s-2]
261 
262  logical :: use_EOS ! If true, compute density from T/S using equation of state.
263  integer :: kb(SZI_(G)) ! The index of the lightest layer denser than the
264  ! buffer layer, or -1 without a bulk mixed layer.
265  logical :: showCallTree ! If true, show the call tree.
266 
267  integer :: i, j, k, is, ie, js, je, nz
268  integer :: isd, ied, jsd, jed
269 
270  real :: kappa_dt_fill ! diffusivity times a timestep used to fill massless layers [Z2 ~> m2]
271 
272  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
273  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
274  showcalltree = calltree_showquery()
275  if (showcalltree) call calltree_enter("set_diffusivity(), MOM_set_diffusivity.F90")
276 
277  if (.not.associated(cs)) call mom_error(fatal,"set_diffusivity: "//&
278  "Module must be initialized before it is used.")
279 
280  if (cs%answers_2018) then
281  ! These hard-coded dimensional parameters are being replaced.
282  kappa_dt_fill = us%m_to_Z**2 * 1.e-3 * 7200.
283  else
284  kappa_dt_fill = cs%Kd_smooth * dt
285  endif
286  omega2 = cs%omega * cs%omega
287 
288  use_eos = associated(tv%eqn_of_state)
289 
290  if ((cs%use_CVMix_ddiff .or. cs%double_diffusion) .and. .not. &
291  (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S))) &
292  call mom_error(fatal, "set_diffusivity: both visc%Kd_extra_T and "//&
293  "visc%Kd_extra_S must be associated when USE_CVMIX_DDIFF or DOUBLE_DIFFUSION are true.")
294 
295  ! Set Kd_lay, Kd_int and Kv_slow to constant values.
296  ! If nothing else is specified, this will be the value used.
297  kd_lay(:,:,:) = cs%Kd
298  kd_int(:,:,:) = cs%Kd
299  if (associated(visc%Kv_slow)) visc%Kv_slow(:,:,:) = cs%Kv
300 
301  ! Set up arrays for diagnostics.
302 
303  if (cs%id_N2 > 0) then
304  allocate(dd%N2_3d(isd:ied,jsd:jed,nz+1)) ; dd%N2_3d(:,:,:) = 0.0
305  endif
306  if (cs%id_Kd_user > 0) then
307  allocate(dd%Kd_user(isd:ied,jsd:jed,nz+1)) ; dd%Kd_user(:,:,:) = 0.0
308  endif
309  if (cs%id_Kd_work > 0) then
310  allocate(dd%Kd_work(isd:ied,jsd:jed,nz)) ; dd%Kd_work(:,:,:) = 0.0
311  endif
312  if (cs%id_maxTKE > 0) then
313  allocate(dd%maxTKE(isd:ied,jsd:jed,nz)) ; dd%maxTKE(:,:,:) = 0.0
314  endif
315  if (cs%id_TKE_to_Kd > 0) then
316  allocate(dd%TKE_to_Kd(isd:ied,jsd:jed,nz)) ; dd%TKE_to_Kd(:,:,:) = 0.0
317  endif
318  if (cs%id_KT_extra > 0) then
319  allocate(dd%KT_extra(isd:ied,jsd:jed,nz+1)) ; dd%KT_extra(:,:,:) = 0.0
320  endif
321  if (cs%id_KS_extra > 0) then
322  allocate(dd%KS_extra(isd:ied,jsd:jed,nz+1)) ; dd%KS_extra(:,:,:) = 0.0
323  endif
324  if (cs%id_Kd_BBL > 0) then
325  allocate(dd%Kd_BBL(isd:ied,jsd:jed,nz+1)) ; dd%Kd_BBL(:,:,:) = 0.0
326  endif
327 
328  ! set up arrays for tidal mixing diagnostics
329  call setup_tidal_diagnostics(g,cs%tm_csp)
330 
331  ! Smooth the properties through massless layers.
332  if (use_eos) then
333  if (cs%debug) then
334  call hchksum(tv%T, "before vert_fill_TS tv%T",g%HI)
335  call hchksum(tv%S, "before vert_fill_TS tv%S",g%HI)
336  call hchksum(h, "before vert_fill_TS h",g%HI, scale=gv%H_to_m)
337  endif
338  call vert_fill_ts(h, tv%T, tv%S, kappa_dt_fill, t_f, s_f, g, gv, larger_h_denom=.true.)
339  if (cs%debug) then
340  call hchksum(tv%T, "after vert_fill_TS tv%T",g%HI)
341  call hchksum(tv%S, "after vert_fill_TS tv%S",g%HI)
342  call hchksum(h, "after vert_fill_TS h",g%HI, scale=gv%H_to_m)
343  endif
344  endif
345 
346  if (cs%useKappaShear) then
347  if (cs%debug) then
348  call hchksum_pair("before calc_KS [uv]_h", u_h, v_h, g%HI, scale=us%L_T_to_m_s)
349  endif
350  call cpu_clock_begin(id_clock_kappashear)
351  if (cs%Vertex_shear) then
352  call full_convection(g, gv, h, tv, t_adj, s_adj, fluxes%p_surf, &
353  (gv%Z_to_H**2)*kappa_dt_fill, halo=1)
354 
355  call calc_kappa_shear_vertex(u, v, h, t_adj, s_adj, tv, fluxes%p_surf, visc%Kd_shear, &
356  visc%TKE_turb, visc%Kv_shear_Bu, dt, g, gv, us, cs%kappaShear_CSp)
357  if (associated(visc%Kv_shear)) visc%Kv_shear(:,:,:) = 0.0 ! needed for other parameterizations
358  if (cs%debug) then
359  call hchksum(visc%Kd_shear, "after calc_KS_vert visc%Kd_shear", g%HI, scale=us%Z2_T_to_m2_s)
360  call bchksum(visc%Kv_shear_Bu, "after calc_KS_vert visc%Kv_shear_Bu", g%HI, scale=us%Z2_T_to_m2_s)
361  call bchksum(visc%TKE_turb, "after calc_KS_vert visc%TKE_turb", g%HI, scale=us%Z_to_m**2*us%s_to_T**2)
362  endif
363  else
364  ! Changes: visc%Kd_shear ; Sets: visc%Kv_shear and visc%TKE_turb
365  call calculate_kappa_shear(u_h, v_h, h, tv, fluxes%p_surf, visc%Kd_shear, visc%TKE_turb, &
366  visc%Kv_shear, dt, g, gv, us, cs%kappaShear_CSp)
367  if (cs%debug) then
368  call hchksum(visc%Kd_shear, "after calc_KS visc%Kd_shear", g%HI, scale=us%Z2_T_to_m2_s)
369  call hchksum(visc%Kv_shear, "after calc_KS visc%Kv_shear", g%HI, scale=us%Z2_T_to_m2_s)
370  call hchksum(visc%TKE_turb, "after calc_KS visc%TKE_turb", g%HI, scale=us%Z_to_m**2*us%s_to_T**2)
371  endif
372  endif
373  call cpu_clock_end(id_clock_kappashear)
374  if (showcalltree) call calltree_waypoint("done with calculate_kappa_shear (set_diffusivity)")
375  elseif (cs%use_CVMix_shear) then
376  !NOTE{BGR}: this needs to be cleaned up. It works in 1D case, but has not been tested outside.
377  call calculate_cvmix_shear(u_h, v_h, h, tv, visc%Kd_shear, visc%Kv_shear, g, gv, us, cs%CVMix_shear_CSp)
378  if (cs%debug) then
379  call hchksum(visc%Kd_shear, "after CVMix_shear visc%Kd_shear", g%HI, scale=us%Z2_T_to_m2_s)
380  call hchksum(visc%Kv_shear, "after CVMix_shear visc%Kv_shear", g%HI, scale=us%Z2_T_to_m2_s)
381  endif
382  elseif (associated(visc%Kv_shear)) then
383  visc%Kv_shear(:,:,:) = 0.0 ! needed if calculate_kappa_shear is not enabled
384  endif
385 
386  ! Calculate the diffusivity, Kd_lay, for each layer. This would be
387  ! the appropriate place to add a depth-dependent parameterization or
388  ! another explicit parameterization of Kd.
389 
390  ! set surface diffusivities (CS%bkgnd_mixing_csp%Kd_sfc)
391  call sfc_bkgnd_mixing(g, us, cs%bkgnd_mixing_csp)
392 
393  !$OMP parallel do default(shared) private(dRho_int, N2_lay, N2_int, N2_bot, KT_extra, &
394  !$OMP KS_extra, TKE_to_Kd, maxTKE, dissip, kb)
395  do j=js,je
396 
397  ! Set up variables related to the stratification.
398  call find_n2(h, tv, t_f, s_f, fluxes, j, g, gv, us, cs, drho_int, n2_lay, n2_int, n2_bot)
399 
400  if (associated(dd%N2_3d)) then
401  do k=1,nz+1 ; do i=is,ie ; dd%N2_3d(i,j,k) = n2_int(i,k) ; enddo ; enddo
402  endif
403 
404  ! Add background mixing
405  call calculate_bkgnd_mixing(h, tv, n2_lay, kd_lay, visc%Kv_slow, j, g, gv, us, cs%bkgnd_mixing_csp)
406 
407  ! Double-diffusion (old method)
408  if (cs%double_diffusion) then
409  call double_diffusion(tv, h, t_f, s_f, j, g, gv, us, cs, kt_extra, ks_extra)
410  do k=2,nz ; do i=is,ie
411  if (ks_extra(i,k) > kt_extra(i,k)) then ! salt fingering
412  kd_lay(i,j,k-1) = kd_lay(i,j,k-1) + 0.5 * kt_extra(i,k)
413  kd_lay(i,j,k) = kd_lay(i,j,k) + 0.5 * kt_extra(i,k)
414  visc%Kd_extra_S(i,j,k) = (ks_extra(i,k) - kt_extra(i,k))
415  visc%Kd_extra_T(i,j,k) = 0.0
416  elseif (kt_extra(i,k) > 0.0) then ! double-diffusive convection
417  kd_lay(i,j,k-1) = kd_lay(i,j,k-1) + 0.5 * ks_extra(i,k)
418  kd_lay(i,j,k) = kd_lay(i,j,k) + 0.5 * ks_extra(i,k)
419  visc%Kd_extra_T(i,j,k) = (kt_extra(i,k) - ks_extra(i,k))
420  visc%Kd_extra_S(i,j,k) = 0.0
421  else ! There is no double diffusion at this interface.
422  visc%Kd_extra_T(i,j,k) = 0.0
423  visc%Kd_extra_S(i,j,k) = 0.0
424  endif
425  enddo ; enddo
426  if (associated(dd%KT_extra)) then ; do k=1,nz+1 ; do i=is,ie
427  dd%KT_extra(i,j,k) = kt_extra(i,k)
428  enddo ; enddo ; endif
429 
430  if (associated(dd%KS_extra)) then ; do k=1,nz+1 ; do i=is,ie
431  dd%KS_extra(i,j,k) = ks_extra(i,k)
432  enddo ; enddo ; endif
433  endif
434 
435  ! Apply double diffusion via CVMix
436  ! GMM, we need to pass HBL to compute_ddiff_coeffs, but it is not yet available.
437  if (cs%use_CVMix_ddiff) then
438  call cpu_clock_begin(id_clock_cvmix_ddiff)
439  call compute_ddiff_coeffs(h, tv, g, gv, us, j, visc%Kd_extra_T, visc%Kd_extra_S, cs%CVMix_ddiff_csp)
440  call cpu_clock_end(id_clock_cvmix_ddiff)
441  endif
442 
443  ! Add the input turbulent diffusivity.
444  if (cs%useKappaShear .or. cs%use_CVMix_shear) then
445  if (present(kd_int)) then
446  do k=2,nz ; do i=is,ie
447  kd_int(i,j,k) = visc%Kd_shear(i,j,k) + 0.5 * (kd_lay(i,j,k-1) + kd_lay(i,j,k))
448  enddo ; enddo
449  do i=is,ie
450  kd_int(i,j,1) = visc%Kd_shear(i,j,1) ! This isn't actually used. It could be 0.
451  kd_int(i,j,nz+1) = 0.0
452  enddo
453  endif
454  do k=1,nz ; do i=is,ie
455  kd_lay(i,j,k) = kd_lay(i,j,k) + 0.5 * (visc%Kd_shear(i,j,k) + visc%Kd_shear(i,j,k+1))
456  enddo ; enddo
457  else
458  if (present(kd_int)) then
459  do i=is,ie
460  kd_int(i,j,1) = kd_lay(i,j,1) ; kd_int(i,j,nz+1) = 0.0
461  enddo
462  do k=2,nz ; do i=is,ie
463  kd_int(i,j,k) = 0.5 * (kd_lay(i,j,k-1) + kd_lay(i,j,k))
464  enddo ; enddo
465  endif
466  endif
467 
468  call find_tke_to_kd(h, tv, drho_int, n2_lay, j, dt, g, gv, us, cs, tke_to_kd, &
469  maxtke, kb)
470  if (associated(dd%maxTKE)) then ; do k=1,nz ; do i=is,ie
471  dd%maxTKE(i,j,k) = maxtke(i,k)
472  enddo ; enddo ; endif
473  if (associated(dd%TKE_to_Kd)) then ; do k=1,nz ; do i=is,ie
474  dd%TKE_to_Kd(i,j,k) = tke_to_kd(i,k)
475  enddo ; enddo ; endif
476 
477  ! Add the ML_Rad diffusivity.
478  if (cs%ML_radiation) &
479  call add_mlrad_diffusivity(h, fluxes, j, g, gv, us, cs, kd_lay, tke_to_kd, kd_int)
480 
481  ! Add the Nikurashin and / or tidal bottom-driven mixing
482  call calculate_tidal_mixing(h, n2_bot, j, tke_to_kd, maxtke, g, gv, us, cs%tm_csp, &
483  n2_lay, n2_int, kd_lay, kd_int, cs%Kd_max, visc%Kv_slow)
484 
485  ! This adds the diffusion sustained by the energy extracted from the flow
486  ! by the bottom drag.
487  if (cs%bottomdraglaw .and. (cs%BBL_effic>0.0)) then
488  if (cs%use_LOTW_BBL_diffusivity) then
489  call add_lotw_bbl_diffusivity(h, u, v, tv, fluxes, visc, j, n2_int, g, gv, us, cs, &
490  kd_lay, kd_int, dd%Kd_BBL)
491  else
492  call add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, tke_to_kd, &
493  maxtke, kb, g, gv, us, cs, kd_lay, kd_int, dd%Kd_BBL)
494  endif
495  endif
496 
497  if (cs%limit_dissipation) then
498  ! This calculates the dissipation ONLY from Kd calculated in this routine
499  ! dissip has units of W/m3 (= kg/m3 * m2/s * 1/s2)
500  ! 1) a global constant,
501  ! 2) a dissipation proportional to N (aka Gargett) and
502  ! 3) dissipation corresponding to a (nearly) constant diffusivity.
503  do k=2,nz-1 ; do i=is,ie
504  dissip = max( cs%dissip_min, & ! Const. floor on dissip.
505  cs%dissip_N0 + cs%dissip_N1 * sqrt(n2_lay(i,k)), & ! Floor aka Gargett
506  cs%dissip_N2 * n2_lay(i,k)) ! Floor of Kd_min*rho0/F_Ri
507  kd_lay(i,j,k) = max(kd_lay(i,j,k) , & ! Apply floor to Kd
508  dissip * (cs%FluxRi_max / (gv%Rho0 * (n2_lay(i,k) + omega2))))
509  enddo ; enddo
510 
511  if (present(kd_int)) then ; do k=2,nz ; do i=is,ie
512  dissip = max( cs%dissip_min, & ! Const. floor on dissip.
513  cs%dissip_N0 + cs%dissip_N1 * sqrt(n2_int(i,k)), & ! Floor aka Gargett
514  cs%dissip_N2 * n2_int(i,k)) ! Floor of Kd_min*rho0/F_Ri
515  kd_int(i,j,k) = max(kd_int(i,j,k) , & ! Apply floor to Kd
516  dissip * (cs%FluxRi_max / (gv%Rho0 * (n2_int(i,k) + omega2))))
517  enddo ; enddo ; endif
518  endif
519 
520  if (associated(dd%Kd_work)) then
521  do k=1,nz ; do i=is,ie
522  dd%Kd_Work(i,j,k) = gv%Rho0 * kd_lay(i,j,k) * n2_lay(i,k) * &
523  gv%H_to_Z*h(i,j,k) ! Watt m-2 s = kg s-3
524  enddo ; enddo
525  endif
526  enddo ! j-loop
527 
528  if (cs%debug) then
529  call hchksum(kd_lay ,"Kd_lay", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
530 
531  if (cs%useKappaShear) call hchksum(visc%Kd_shear, "Turbulent Kd", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
532 
533  if (cs%use_CVMix_ddiff) then
534  call hchksum(visc%Kd_extra_T, "MOM_set_diffusivity: Kd_extra_T", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
535  call hchksum(visc%Kd_extra_S, "MOM_set_diffusivity: Kd_extra_S", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
536  endif
537 
538  if (associated(visc%kv_bbl_u) .and. associated(visc%kv_bbl_v)) then
539  call uvchksum("BBL Kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, &
540  g%HI, 0, symmetric=.true., scale=us%Z2_T_to_m2_s)
541  endif
542 
543  if (associated(visc%bbl_thick_u) .and. associated(visc%bbl_thick_v)) then
544  call uvchksum("BBL bbl_thick_[uv]", visc%bbl_thick_u, &
545  visc%bbl_thick_v, g%HI, 0, symmetric=.true., scale=us%Z_to_m)
546  endif
547 
548  if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) then
549  call uvchksum("Ray_[uv]", visc%Ray_u, visc%Ray_v, g%HI, 0, symmetric=.true., scale=us%Z_to_m*us%s_to_T)
550  endif
551 
552  endif
553 
554  if (cs%Kd_add > 0.0) then
555  if (present(kd_int)) then
556  !$OMP parallel do default(shared)
557  do k=1,nz ; do j=js,je ; do i=is,ie
558  kd_int(i,j,k) = kd_int(i,j,k) + cs%Kd_add
559  kd_lay(i,j,k) = kd_lay(i,j,k) + cs%Kd_add
560  enddo ; enddo ; enddo
561  else
562  !$OMP parallel do default(shared)
563  do k=1,nz ; do j=js,je ; do i=is,ie
564  kd_lay(i,j,k) = kd_lay(i,j,k) + cs%Kd_add
565  enddo ; enddo ; enddo
566  endif
567  endif
568 
569  if (cs%user_change_diff) then
570  call user_change_diff(h, tv, g, gv, cs%user_change_diff_CSp, kd_lay, kd_int, &
571  t_f, s_f, dd%Kd_user)
572  endif
573 
574  ! post diagnostics
575 
576  ! background mixing
577  if (cs%bkgnd_mixing_csp%id_kd_bkgnd > 0) &
578  call post_data(cs%bkgnd_mixing_csp%id_kd_bkgnd, cs%bkgnd_mixing_csp%kd_bkgnd, cs%bkgnd_mixing_csp%diag)
579  if (cs%bkgnd_mixing_csp%id_kv_bkgnd > 0) &
580  call post_data(cs%bkgnd_mixing_csp%id_kv_bkgnd, cs%bkgnd_mixing_csp%kv_bkgnd, cs%bkgnd_mixing_csp%diag)
581 
582  ! double diffusive mixing
583  if (cs%CVMix_ddiff_csp%id_KT_extra > 0) &
584  call post_data(cs%CVMix_ddiff_csp%id_KT_extra, visc%Kd_extra_T, cs%CVMix_ddiff_csp%diag)
585  if (cs%CVMix_ddiff_csp%id_KS_extra > 0) &
586  call post_data(cs%CVMix_ddiff_csp%id_KS_extra, visc%Kd_extra_S, cs%CVMix_ddiff_csp%diag)
587  if (cs%CVMix_ddiff_csp%id_R_rho > 0) &
588  call post_data(cs%CVMix_ddiff_csp%id_R_rho, cs%CVMix_ddiff_csp%R_rho, cs%CVMix_ddiff_csp%diag)
589 
590  if (cs%id_Kd_layer > 0) call post_data(cs%id_Kd_layer, kd_lay, cs%diag)
591 
592  ! tidal mixing
593  call post_tidal_diagnostics(g,gv,h,cs%tm_csp)
594 
595  if (cs%tm_csp%Int_tide_dissipation .or. cs%tm_csp%Lee_wave_dissipation .or. &
596  cs%tm_csp%Lowmode_itidal_dissipation) then
597 
598  if (cs%id_N2 > 0) call post_data(cs%id_N2, dd%N2_3d, cs%diag)
599  if (cs%id_Kd_user > 0) call post_data(cs%id_Kd_user, dd%Kd_user, cs%diag)
600  if (cs%id_Kd_Work > 0) call post_data(cs%id_Kd_Work, dd%Kd_Work, cs%diag)
601  if (cs%id_maxTKE > 0) call post_data(cs%id_maxTKE, dd%maxTKE, cs%diag)
602  if (cs%id_TKE_to_Kd > 0) call post_data(cs%id_TKE_to_Kd, dd%TKE_to_Kd, cs%diag)
603  endif
604 
605  if (cs%id_KT_extra > 0) call post_data(cs%id_KT_extra, dd%KT_extra, cs%diag)
606  if (cs%id_KS_extra > 0) call post_data(cs%id_KS_extra, dd%KS_extra, cs%diag)
607  if (cs%id_Kd_BBL > 0) call post_data(cs%id_Kd_BBL, dd%Kd_BBL, cs%diag)
608 
609  if (associated(dd%N2_3d)) deallocate(dd%N2_3d)
610  if (associated(dd%Kd_work)) deallocate(dd%Kd_work)
611  if (associated(dd%Kd_user)) deallocate(dd%Kd_user)
612  if (associated(dd%maxTKE)) deallocate(dd%maxTKE)
613  if (associated(dd%TKE_to_Kd)) deallocate(dd%TKE_to_Kd)
614  if (associated(dd%KT_extra)) deallocate(dd%KT_extra)
615  if (associated(dd%KS_extra)) deallocate(dd%KS_extra)
616  if (associated(dd%Kd_BBL)) deallocate(dd%Kd_BBL)
617 
618  if (showcalltree) call calltree_leave("set_diffusivity()")
619 

References add_drag_diffusivity(), add_lotw_bbl_diffusivity(), add_mlrad_diffusivity(), mom_kappa_shear::calc_kappa_shear_vertex(), mom_error_handler::calltree_enter(), mom_error_handler::calltree_leave(), mom_error_handler::calltree_waypoint(), mom_cvmix_ddiff::compute_ddiff_coeffs(), double_diffusion(), find_n2(), find_tke_to_kd(), mom_full_convection::full_convection(), id_clock_cvmix_ddiff, id_clock_kappashear, mom_tidal_mixing::post_tidal_diagnostics(), mom_tidal_mixing::setup_tidal_diagnostics(), mom_bkgnd_mixing::sfc_bkgnd_mixing(), and mom_isopycnal_slopes::vert_fill_ts().

Here is the call graph for this function:

◆ set_diffusivity_end()

subroutine, public mom_set_diffusivity::set_diffusivity_end ( type(set_diffusivity_cs), pointer  CS)

Clear pointers and dealocate memory.

Parameters
csControl structure for this module

Definition at line 2207 of file MOM_set_diffusivity.F90.

2207  type(set_diffusivity_CS), pointer :: CS !< Control structure for this module
2208 
2209  if (.not.associated(cs)) return
2210 
2211  call bkgnd_mixing_end(cs%bkgnd_mixing_csp)
2212 
2213  if (cs%user_change_diff) &
2214  call user_change_diff_end(cs%user_change_diff_CSp)
2215 
2216  if (cs%use_CVMix_shear) &
2217  call cvmix_shear_end(cs%CVMix_shear_csp)
2218 
2219  if (cs%use_CVMix_ddiff) &
2220  call cvmix_ddiff_end(cs%CVMix_ddiff_csp)
2221 
2222  if (associated(cs)) deallocate(cs)
2223 

References mom_bkgnd_mixing::bkgnd_mixing_end(), mom_cvmix_shear::cvmix_shear_end(), and user_change_diffusivity::user_change_diff_end().

Here is the call graph for this function:

◆ set_diffusivity_init()

subroutine, public mom_set_diffusivity::set_diffusivity_init ( type(time_type), intent(in)  Time,
type(ocean_grid_type), intent(inout)  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(set_diffusivity_cs), pointer  CS,
type(int_tide_cs), pointer  int_tide_CSp,
type(tidal_mixing_cs), pointer  tm_CSp,
integer, intent(out), optional  halo_TS 
)
Parameters
[in]timeThe current model time
[in,out]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]param_fileA structure to parse for run-time parameters.
[in,out]diagA structure used to regulate diagnostic output.
cspointer set to point to the module control structure.
int_tide_csppointer to the internal tides control structure (BDM)
tm_csppointer to tidal mixing control structure
[out]halo_tsThe halo size of tracer points that must be valid for the calculations in set_diffusivity.

Definition at line 1878 of file MOM_set_diffusivity.F90.

1878  type(time_type), intent(in) :: Time !< The current model time
1879  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
1880  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1881  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1882  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
1883  !! parameters.
1884  type(diag_ctrl), target, intent(inout) :: diag !< A structure used to regulate diagnostic output.
1885  type(set_diffusivity_CS), pointer :: CS !< pointer set to point to the module control
1886  !! structure.
1887  type(int_tide_CS), pointer :: int_tide_CSp !< pointer to the internal tides control
1888  !! structure (BDM)
1889  type(tidal_mixing_cs), pointer :: tm_csp !< pointer to tidal mixing control
1890  !! structure
1891  integer, optional, intent(out) :: halo_TS !< The halo size of tracer points that must be
1892  !! valid for the calculations in set_diffusivity.
1893 
1894  ! Local variables
1895  real :: decay_length
1896  logical :: ML_use_omega
1897  logical :: default_2018_answers
1898  ! This include declares and sets the variable "version".
1899 # include "version_variable.h"
1900  character(len=40) :: mdl = "MOM_set_diffusivity" ! This module's name.
1901  real :: omega_frac_dflt
1902  integer :: i, j, is, ie, js, je
1903  integer :: isd, ied, jsd, jed
1904 
1905  if (associated(cs)) then
1906  call mom_error(warning, "diabatic_entrain_init called with an associated "// &
1907  "control structure.")
1908  return
1909  endif
1910  allocate(cs)
1911 
1912  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1913  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1914 
1915  cs%diag => diag
1916  if (associated(int_tide_csp)) cs%int_tide_CSp => int_tide_csp
1917  if (associated(tm_csp)) cs%tm_csp => tm_csp
1918 
1919  ! These default values always need to be set.
1920  cs%BBL_mixing_as_max = .true.
1921  cs%cdrag = 0.003 ; cs%BBL_effic = 0.0
1922  cs%bulkmixedlayer = (gv%nkml > 0)
1923 
1924  ! Read all relevant parameters and write them to the model log.
1925  call log_version(param_file, mdl, version, "")
1926 
1927  call get_param(param_file, mdl, "INPUTDIR", cs%inputdir, default=".")
1928  cs%inputdir = slasher(cs%inputdir)
1929  call get_param(param_file, mdl, "FLUX_RI_MAX", cs%FluxRi_max, &
1930  "The flux Richardson number where the stratification is "//&
1931  "large enough that N2 > omega2. The full expression for "//&
1932  "the Flux Richardson number is usually "//&
1933  "FLUX_RI_MAX*N2/(N2+OMEGA2).", default=0.2)
1934  call get_param(param_file, mdl, "OMEGA", cs%omega, &
1935  "The rotation rate of the earth.", units="s-1", &
1936  default=7.2921e-5, scale=us%T_to_s)
1937 
1938  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
1939  "This sets the default value for the various _2018_ANSWERS parameters.", &
1940  default=.true.)
1941  call get_param(param_file, mdl, "SET_DIFF_2018_ANSWERS", cs%answers_2018, &
1942  "If true, use the order of arithmetic and expressions that recover the "//&
1943  "answers from the end of 2018. Otherwise, use updated and more robust "//&
1944  "forms of the same expressions.", default=default_2018_answers)
1945 
1946  call get_param(param_file, mdl, "ML_RADIATION", cs%ML_radiation, &
1947  "If true, allow a fraction of TKE available from wind "//&
1948  "work to penetrate below the base of the mixed layer "//&
1949  "with a vertical decay scale determined by the minimum "//&
1950  "of: (1) The depth of the mixed layer, (2) an Ekman "//&
1951  "length scale.", default=.false.)
1952  if (cs%ML_radiation) then
1953  ! This give a minimum decay scale that is typically much less than Angstrom.
1954  cs%ustar_min = 2e-4 * cs%omega * (gv%Angstrom_Z + gv%H_subroundoff*gv%H_to_Z)
1955 
1956  call get_param(param_file, mdl, "ML_RAD_EFOLD_COEFF", cs%ML_rad_efold_coeff, &
1957  "A coefficient that is used to scale the penetration "//&
1958  "depth for turbulence below the base of the mixed layer. "//&
1959  "This is only used if ML_RADIATION is true.", units="nondim", &
1960  default=0.2)
1961  call get_param(param_file, mdl, "ML_RAD_BUG", cs%ML_rad_bug, &
1962  "If true use code with a bug that reduces the energy available "//&
1963  "in the transition layer by a factor of the inverse of the energy "//&
1964  "deposition lenthscale (in m).", default=.true.)
1965  call get_param(param_file, mdl, "ML_RAD_KD_MAX", cs%ML_rad_kd_max, &
1966  "The maximum diapycnal diffusivity due to turbulence "//&
1967  "radiated from the base of the mixed layer. "//&
1968  "This is only used if ML_RADIATION is true.", &
1969  units="m2 s-1", default=1.0e-3, &
1970  scale=us%m2_s_to_Z2_T)
1971  call get_param(param_file, mdl, "ML_RAD_COEFF", cs%ML_rad_coeff, &
1972  "The coefficient which scales MSTAR*USTAR^3 to obtain "//&
1973  "the energy available for mixing below the base of the "//&
1974  "mixed layer. This is only used if ML_RADIATION is true.", &
1975  units="nondim", default=0.2)
1976  call get_param(param_file, mdl, "ML_RAD_APPLY_TKE_DECAY", cs%ML_rad_TKE_decay, &
1977  "If true, apply the same exponential decay to ML_rad as "//&
1978  "is applied to the other surface sources of TKE in the "//&
1979  "mixed layer code. This is only used if ML_RADIATION is true.", &
1980  default=.true.)
1981  call get_param(param_file, mdl, "MSTAR", cs%mstar, &
1982  "The ratio of the friction velocity cubed to the TKE "//&
1983  "input to the mixed layer.", "units=nondim", default=1.2)
1984  call get_param(param_file, mdl, "TKE_DECAY", cs%TKE_decay, &
1985  "The ratio of the natural Ekman depth to the TKE decay scale.", &
1986  units="nondim", default=2.5)
1987  call get_param(param_file, mdl, "ML_USE_OMEGA", ml_use_omega, &
1988  "If true, use the absolute rotation rate instead of the "//&
1989  "vertical component of rotation when setting the decay "//&
1990  "scale for turbulence.", default=.false., do_not_log=.true.)
1991  omega_frac_dflt = 0.0
1992  if (ml_use_omega) then
1993  call mom_error(warning, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
1994  omega_frac_dflt = 1.0
1995  endif
1996  call get_param(param_file, mdl, "ML_OMEGA_FRAC", cs%ML_omega_frac, &
1997  "When setting the decay scale for turbulence, use this "//&
1998  "fraction of the absolute rotation rate blended with the "//&
1999  "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
2000  units="nondim", default=omega_frac_dflt)
2001  endif
2002 
2003  call get_param(param_file, mdl, "BOTTOMDRAGLAW", cs%bottomdraglaw, &
2004  "If true, the bottom stress is calculated with a drag "//&
2005  "law of the form c_drag*|u|*u. The velocity magnitude "//&
2006  "may be an assumed value or it may be based on the "//&
2007  "actual velocity in the bottommost HBBL, depending on "//&
2008  "LINEAR_DRAG.", default=.true.)
2009  if (cs%bottomdraglaw) then
2010  call get_param(param_file, mdl, "CDRAG", cs%cdrag, &
2011  "The drag coefficient relating the magnitude of the "//&
2012  "velocity field to the bottom stress. CDRAG is only used "//&
2013  "if BOTTOMDRAGLAW is true.", units="nondim", default=0.003)
2014  call get_param(param_file, mdl, "BBL_EFFIC", cs%BBL_effic, &
2015  "The efficiency with which the energy extracted by "//&
2016  "bottom drag drives BBL diffusion. This is only "//&
2017  "used if BOTTOMDRAGLAW is true.", units="nondim", default=0.20)
2018  call get_param(param_file, mdl, "BBL_MIXING_MAX_DECAY", decay_length, &
2019  "The maximum decay scale for the BBL diffusion, or 0 to allow the mixing "//&
2020  "to penetrate as far as stratification and rotation permit. The default "//&
2021  "for now is 200 m. This is only used if BOTTOMDRAGLAW is true.", &
2022  units="m", default=200.0, scale=us%m_to_Z)
2023 
2024  cs%IMax_decay = 0.0
2025  if (decay_length > 0.0) cs%IMax_decay = 1.0/decay_length
2026  call get_param(param_file, mdl, "BBL_MIXING_AS_MAX", cs%BBL_mixing_as_max, &
2027  "If true, take the maximum of the diffusivity from the "//&
2028  "BBL mixing and the other diffusivities. Otherwise, "//&
2029  "diffusivity from the BBL_mixing is simply added.", &
2030  default=.true.)
2031  call get_param(param_file, mdl, "USE_LOTW_BBL_DIFFUSIVITY", cs%use_LOTW_BBL_diffusivity, &
2032  "If true, uses a simple, imprecise but non-coordinate dependent, model "//&
2033  "of BBL mixing diffusivity based on Law of the Wall. Otherwise, uses "//&
2034  "the original BBL scheme.", default=.false.)
2035  if (cs%use_LOTW_BBL_diffusivity) then
2036  call get_param(param_file, mdl, "LOTW_BBL_USE_OMEGA", cs%LOTW_BBL_use_omega, &
2037  "If true, use the maximum of Omega and N for the TKE to diffusion "//&
2038  "calculation. Otherwise, N is N.", default=.true.)
2039  endif
2040  else
2041  cs%use_LOTW_BBL_diffusivity = .false. ! This parameterization depends on a u* from viscous BBL
2042  endif
2043  cs%id_Kd_BBL = register_diag_field('ocean_model', 'Kd_BBL', diag%axesTi, time, &
2044  'Bottom Boundary Layer Diffusivity', 'm2 s-1', &
2045  conversion=us%Z2_T_to_m2_s)
2046  call get_param(param_file, mdl, "SIMPLE_TKE_TO_KD", cs%simple_TKE_to_Kd, &
2047  "If true, uses a simple estimate of Kd/TKE that will "//&
2048  "work for arbitrary vertical coordinates. If false, "//&
2049  "calculates Kd/TKE and bounds based on exact energetics "//&
2050  "for an isopycnal layer-formulation.", &
2051  default=.false.)
2052 
2053  ! set params releted to the background mixing
2054  call bkgnd_mixing_init(time, g, gv, us, param_file, cs%diag, cs%bkgnd_mixing_csp)
2055 
2056  call get_param(param_file, mdl, "KV", cs%Kv, &
2057  "The background kinematic viscosity in the interior. "//&
2058  "The molecular value, ~1e-6 m2 s-1, may be used.", &
2059  units="m2 s-1", scale=us%m2_s_to_Z2_T, &
2060  fail_if_missing=.true.)
2061 
2062  call get_param(param_file, mdl, "KD", cs%Kd, &
2063  "The background diapycnal diffusivity of density in the "//&
2064  "interior. Zero or the molecular value, ~1e-7 m2 s-1, "//&
2065  "may be used.", units="m2 s-1", scale=us%m2_s_to_Z2_T, &
2066  fail_if_missing=.true.)
2067  call get_param(param_file, mdl, "KD_MIN", cs%Kd_min, &
2068  "The minimum diapycnal diffusivity.", &
2069  units="m2 s-1", default=0.01*cs%Kd*us%Z2_T_to_m2_s, &
2070  scale=us%m2_s_to_Z2_T)
2071  call get_param(param_file, mdl, "KD_MAX", cs%Kd_max, &
2072  "The maximum permitted increment for the diapycnal "//&
2073  "diffusivity from TKE-based parameterizations, or a "//&
2074  "negative value for no limit.", units="m2 s-1", default=-1.0, &
2075  scale=us%m2_s_to_Z2_T)
2076  if (cs%simple_TKE_to_Kd .and. cs%Kd_max<=0.) call mom_error(fatal, &
2077  "set_diffusivity_init: To use SIMPLE_TKE_TO_KD, KD_MAX must be set to >0.")
2078  call get_param(param_file, mdl, "KD_ADD", cs%Kd_add, &
2079  "A uniform diapycnal diffusivity that is added "//&
2080  "everywhere without any filtering or scaling.", &
2081  units="m2 s-1", default=0.0, scale=us%m2_s_to_Z2_T)
2082  if (cs%use_LOTW_BBL_diffusivity .and. cs%Kd_max<=0.) call mom_error(fatal, &
2083  "set_diffusivity_init: KD_MAX must be set (positive) when "// &
2084  "USE_LOTW_BBL_DIFFUSIVITY=True.")
2085  call get_param(param_file, mdl, "KD_SMOOTH", cs%Kd_smooth, &
2086  "A diapycnal diffusivity that is used to interpolate "//&
2087  "more sensible values of T & S into thin layers.", &
2088  default=1.0e-6, scale=us%m_to_Z**2*us%T_to_s)
2089 
2090  call get_param(param_file, mdl, "DEBUG", cs%debug, &
2091  "If true, write out verbose debugging data.", &
2092  default=.false., debuggingparam=.true.)
2093 
2094  call get_param(param_file, mdl, "USER_CHANGE_DIFFUSIVITY", cs%user_change_diff, &
2095  "If true, call user-defined code to change the diffusivity.", &
2096  default=.false.)
2097 
2098  call get_param(param_file, mdl, "DISSIPATION_MIN", cs%dissip_min, &
2099  "The minimum dissipation by which to determine a lower "//&
2100  "bound of Kd (a floor).", units="W m-3", default=0.0, &
2101  scale=us%kg_m3_to_R*us%m2_s_to_Z2_T*(us%T_to_s**2))
2102  call get_param(param_file, mdl, "DISSIPATION_N0", cs%dissip_N0, &
2103  "The intercept when N=0 of the N-dependent expression "//&
2104  "used to set a minimum dissipation by which to determine "//&
2105  "a lower bound of Kd (a floor): A in eps_min = A + B*N.", &
2106  units="W m-3", default=0.0, &
2107  scale=us%kg_m3_to_R*us%m2_s_to_Z2_T*(us%T_to_s**2))
2108  call get_param(param_file, mdl, "DISSIPATION_N1", cs%dissip_N1, &
2109  "The coefficient multiplying N, following Gargett, used to "//&
2110  "set a minimum dissipation by which to determine a lower "//&
2111  "bound of Kd (a floor): B in eps_min = A + B*N", &
2112  units="J m-3", default=0.0, scale=us%kg_m3_to_R*us%m2_s_to_Z2_T*us%T_to_s)
2113  call get_param(param_file, mdl, "DISSIPATION_KD_MIN", cs%dissip_Kd_min, &
2114  "The minimum vertical diffusivity applied as a floor.", &
2115  units="m2 s-1", default=0.0, scale=us%m2_s_to_Z2_T)
2116 
2117  cs%limit_dissipation = (cs%dissip_min>0.) .or. (cs%dissip_N1>0.) .or. &
2118  (cs%dissip_N0>0.) .or. (cs%dissip_Kd_min>0.)
2119  cs%dissip_N2 = 0.0
2120  if (cs%FluxRi_max > 0.0) &
2121  cs%dissip_N2 = cs%dissip_Kd_min * gv%Rho0 / cs%FluxRi_max
2122 
2123  cs%id_Kd_layer = register_diag_field('ocean_model', 'Kd_layer', diag%axesTL, time, &
2124  'Diapycnal diffusivity of layers (as set)', 'm2 s-1', &
2125  conversion=us%Z2_T_to_m2_s)
2126 
2127 
2128  if (cs%tm_csp%Int_tide_dissipation .or. cs%tm_csp%Lee_wave_dissipation .or. &
2129  cs%tm_csp%Lowmode_itidal_dissipation) then
2130 
2131  cs%id_Kd_Work = register_diag_field('ocean_model', 'Kd_Work', diag%axesTL, time, &
2132  'Work done by Diapycnal Mixing', 'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2133  cs%id_maxTKE = register_diag_field('ocean_model', 'maxTKE', diag%axesTL, time, &
2134  'Maximum layer TKE', 'm3 s-3', conversion=(us%Z_to_m**3*us%s_to_T**3))
2135  cs%id_TKE_to_Kd = register_diag_field('ocean_model', 'TKE_to_Kd', diag%axesTL, time, &
2136  'Convert TKE to Kd', 's2 m', &
2137  conversion=us%Z2_T_to_m2_s*(us%m_to_Z**3*us%T_to_s**3))
2138  cs%id_N2 = register_diag_field('ocean_model', 'N2', diag%axesTi, time, &
2139  'Buoyancy frequency squared', 's-2', cmor_field_name='obvfsq', &
2140  cmor_long_name='Square of seawater buoyancy frequency', &
2141  cmor_standard_name='square_of_brunt_vaisala_frequency_in_sea_water', &
2142  conversion=us%s_to_T**2)
2143 
2144  if (cs%user_change_diff) &
2145  cs%id_Kd_user = register_diag_field('ocean_model', 'Kd_user', diag%axesTi, time, &
2146  'User-specified Extra Diffusivity', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
2147  endif
2148 
2149  call get_param(param_file, mdl, "DOUBLE_DIFFUSION", cs%double_diffusion, &
2150  "If true, increase diffusivites for temperature or salt "//&
2151  "based on double-diffusive parameterization from MOM4/KPP.", &
2152  default=.false.)
2153 
2154  if (cs%double_diffusion) then
2155  call get_param(param_file, mdl, "MAX_RRHO_SALT_FINGERS", cs%Max_Rrho_salt_fingers, &
2156  "Maximum density ratio for salt fingering regime.", &
2157  default=2.55, units="nondim")
2158  call get_param(param_file, mdl, "MAX_SALT_DIFF_SALT_FINGERS", cs%Max_salt_diff_salt_fingers, &
2159  "Maximum salt diffusivity for salt fingering regime.", &
2160  default=1.e-4, units="m2 s-1", scale=us%m2_s_to_Z2_T)
2161  call get_param(param_file, mdl, "KV_MOLECULAR", cs%Kv_molecular, &
2162  "Molecular viscosity for calculation of fluxes under "//&
2163  "double-diffusive convection.", default=1.5e-6, units="m2 s-1", &
2164  scale=us%m2_s_to_Z2_T)
2165  ! The default molecular viscosity follows the CCSM4.0 and MOM4p1 defaults.
2166 
2167  cs%id_KT_extra = register_diag_field('ocean_model', 'KT_extra', diag%axesTi, time, &
2168  'Double-diffusive diffusivity for temperature', 'm2 s-1', &
2169  conversion=us%Z2_T_to_m2_s)
2170 
2171  cs%id_KS_extra = register_diag_field('ocean_model', 'KS_extra', diag%axesTi, time, &
2172  'Double-diffusive diffusivity for salinity', 'm2 s-1', &
2173  conversion=us%Z2_T_to_m2_s)
2174  endif ! old double-diffusion
2175 
2176  if (cs%user_change_diff) then
2177  call user_change_diff_init(time, g, gv, us, param_file, diag, cs%user_change_diff_CSp)
2178  endif
2179 
2180  if (cs%tm_csp%Int_tide_dissipation .and. cs%bkgnd_mixing_csp%Bryan_Lewis_diffusivity) &
2181  call mom_error(fatal,"MOM_Set_Diffusivity: "// &
2182  "Bryan-Lewis and internal tidal dissipation are both enabled. Choose one.")
2183 
2184  cs%useKappaShear = kappa_shear_init(time, g, gv, us, param_file, cs%diag, cs%kappaShear_CSp)
2185  if (cs%useKappaShear) cs%Vertex_Shear = kappa_shear_at_vertex(param_file)
2186 
2187  if (cs%useKappaShear) &
2188  id_clock_kappashear = cpu_clock_id('(Ocean kappa_shear)', grain=clock_module)
2189 
2190  ! CVMix shear-driven mixing
2191  cs%use_CVMix_shear = cvmix_shear_init(time, g, gv, us, param_file, cs%diag, cs%CVMix_shear_csp)
2192 
2193  ! CVMix double diffusion mixing
2194  cs%use_CVMix_ddiff = cvmix_ddiff_init(time, g, gv, us, param_file, cs%diag, cs%CVMix_ddiff_csp)
2195  if (cs%use_CVMix_ddiff) &
2196  id_clock_cvmix_ddiff = cpu_clock_id('(Double diffusion via CVMix)', grain=clock_module)
2197 
2198  if (present(halo_ts)) then
2199  halo_ts = 0
2200  if (cs%Vertex_Shear) halo_ts = 1
2201  endif
2202 

References id_clock_cvmix_ddiff, id_clock_kappashear, mom_kappa_shear::kappa_shear_at_vertex(), and mom_diag_mediator::register_diag_field().

Here is the call graph for this function:

Variable Documentation

◆ id_clock_cvmix_ddiff

integer mom_set_diffusivity::id_clock_cvmix_ddiff
private

CPU time clocks.

Definition at line 189 of file MOM_set_diffusivity.F90.

Referenced by set_diffusivity(), and set_diffusivity_init().

◆ id_clock_kappashear

integer mom_set_diffusivity::id_clock_kappashear
private

CPU time clocks.

Definition at line 189 of file MOM_set_diffusivity.F90.

189 integer :: id_clock_kappaShear, id_clock_CVMix_ddiff

Referenced by set_diffusivity(), and set_diffusivity_init().