MOM6
Data Types | Functions/Subroutines
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...
 

Functions/Subroutines

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_BBL, Kd_lay, Kd_int)
 Calculates a BBL diffusivity use a Prandtl number 1 diffusivity 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, TKE_to_Kd, Kd_lay, 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, OBC)
 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, halo_TS)
 
subroutine, public set_diffusivity_end (CS)
 Clear pointers and dealocate memory. More...
 

Variables

integer id_clock_kappashear
 CPU time clocks.
 
integer id_clock_cvmix_ddiff
 CPU time clocks.
 

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 %ke), intent(inout)  Kd_lay,
real, dimension( g %isd: g %ied, g %ke+1), intent(inout), optional  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 [Z3 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 diffusivity in layers, [Z2 T-1 ~> m2 s-1].
[in,out]kd_intThe diapycnal diffusivity at interfaces,
kd_bblInterface BBL diffusivity [Z2 T-1 ~> m2 s-1].

Definition at line 1165 of file MOM_set_diffusivity.F90.

1165  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
1166  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
1167  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1168  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1169  intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1]
1170  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1171  intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]
1172  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1173  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1174  type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
1175  !! thermodynamic fields.
1176  type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes
1177  type(vertvisc_type), intent(in) :: visc !< Structure containing vertical viscosities, bottom
1178  !! boundary layer properies, and related fields
1179  integer, intent(in) :: j !< j-index of row to work on
1180  real, dimension(SZI_(G),SZK_(G)), intent(in) :: tke_to_kd !< The conversion rate between the TKE
1181  !! TKE dissipated within a layer and the
1182  !! diapycnal diffusivity witin that layer,
1183  !! usually (~Rho_0 / (G_Earth * dRho_lay))
1184  !! [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1]
1185  real, dimension(SZI_(G),SZK_(G)), intent(in) :: maxtke !< The energy required to for a layer to entrain
1186  !! to its maximum-realizable thickness [Z3 T-3 ~> m3 s-3]
1187  integer, dimension(SZI_(G)), intent(in) :: kb !< Index of lightest layer denser than the buffer
1188  !! layer, or -1 without a bulk mixed layer
1189  type(set_diffusivity_cs), pointer :: cs !< Diffusivity control structure
1190  real, dimension(SZI_(G),SZK_(G)), intent(inout) :: kd_lay !< The diapycnal diffusivity in layers,
1191  !! [Z2 T-1 ~> m2 s-1].
1192  real, dimension(SZI_(G),SZK_(G)+1), &
1193  optional, intent(inout) :: kd_int !< The diapycnal diffusivity at interfaces,
1194  !! [Z2 T-1 ~> m2 s-1].
1195  real, dimension(:,:,:), pointer :: kd_bbl !< Interface BBL diffusivity [Z2 T-1 ~> m2 s-1].
1196 
1197 ! This routine adds diffusion sustained by flow energy extracted by bottom drag.
1198 
1199  real, dimension(SZK_(G)+1) :: &
1200  rint ! coordinate density of an interface [R ~> kg m-3]
1201  real, dimension(SZI_(G)) :: &
1202  htot, & ! total thickness above or below a layer, or the
1203  ! integrated thickness in the BBL [Z ~> m].
1204  rho_htot, & ! running integral with depth of density [Z R ~> kg m-2]
1205  gh_sum_top, & ! BBL value of g'h that can be supported by
1206  ! the local ustar, times R0_g [R ~> kg m-2]
1207  rho_top, & ! density at top of the BBL [R ~> kg m-3]
1208  tke, & ! turbulent kinetic energy available to drive
1209  ! bottom-boundary layer mixing in a layer [Z3 T-3 ~> m3 s-3]
1210  i2decay ! inverse of twice the TKE decay scale [Z-1 ~> m-1].
1211 
1212  real :: tke_to_layer ! TKE used to drive mixing in a layer [Z3 T-3 ~> m3 s-3]
1213  real :: tke_ray ! TKE from layer Rayleigh drag used to drive mixing in layer [Z3 T-3 ~> m3 s-3]
1214  real :: tke_here ! TKE that goes into mixing in this layer [Z3 T-3 ~> m3 s-3]
1215  real :: drl, drbot ! temporaries holding density differences [R ~> kg m-3]
1216  real :: cdrag_sqrt ! square root of the drag coefficient [nondim]
1217  real :: ustar_h ! value of ustar at a thickness point [Z T-1 ~> m s-1].
1218  real :: absf ! average absolute Coriolis parameter around a thickness point [T-1 ~> s-1]
1219  real :: r0_g ! Rho0 / G_Earth [R T2 Z-1 m-1 ~> kg s2 m-5]
1220  real :: i_rho0 ! 1 / RHO0 [R-1 ~> m3 kg-1]
1221  real :: delta_kd ! increment to Kd from the bottom boundary layer mixing [Z2 T-1 ~> m2 s-1].
1222  logical :: rayleigh_drag ! Set to true if Rayleigh drag velocities
1223  ! defined in visc, on the assumption that this
1224  ! extracted energy also drives diapycnal mixing.
1225 
1226  logical :: domore, do_i(szi_(g))
1227  logical :: do_diag_kd_bbl
1228 
1229  integer :: i, k, is, ie, nz, i_rem, kb_min
1230  is = g%isc ; ie = g%iec ; nz = g%ke
1231 
1232  do_diag_kd_bbl = associated(kd_bbl)
1233 
1234  if (.not.(cs%bottomdraglaw .and. (cs%BBL_effic>0.0))) return
1235 
1236  cdrag_sqrt = sqrt(cs%cdrag)
1237  tke_ray = 0.0 ; rayleigh_drag = .false.
1238  if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) rayleigh_drag = .true.
1239 
1240  i_rho0 = 1.0 / (gv%Rho0)
1241  r0_g = gv%Rho0 / (us%L_to_Z**2 * gv%g_Earth)
1242 
1243  do k=2,nz ; rint(k) = 0.5*(gv%Rlay(k-1)+gv%Rlay(k)) ; enddo
1244 
1245  kb_min = max(gv%nk_rho_varies+1,2)
1246 
1247  ! The turbulence decay scale is 0.5*ustar/f from K&E & MOM_vertvisc.F90
1248  ! Any turbulence that makes it into the mixed layers is assumed
1249  ! to be relatively small and is discarded.
1250  do i=is,ie
1251  ustar_h = visc%ustar_BBL(i,j)
1252  if (associated(fluxes%ustar_tidal)) &
1253  ustar_h = ustar_h + fluxes%ustar_tidal(i,j)
1254  absf = 0.25 * ((abs(g%CoriolisBu(i-1,j-1)) + abs(g%CoriolisBu(i,j))) + &
1255  (abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j-1))))
1256  if ((ustar_h > 0.0) .and. (absf > 0.5*cs%IMax_decay*ustar_h)) then
1257  i2decay(i) = absf / ustar_h
1258  else
1259  ! The maximum decay scale should be something of order 200 m.
1260  ! If ustar_h = 0, this is land so this value doesn't matter.
1261  i2decay(i) = 0.5*cs%IMax_decay
1262  endif
1263  tke(i) = ((cs%BBL_effic * cdrag_sqrt) * exp(-i2decay(i)*(gv%H_to_Z*h(i,j,nz))) ) * &
1264  visc%TKE_BBL(i,j)
1265 
1266  if (associated(fluxes%TKE_tidal)) &
1267  tke(i) = tke(i) + fluxes%TKE_tidal(i,j) * i_rho0 * &
1268  (cs%BBL_effic * exp(-i2decay(i)*(gv%H_to_Z*h(i,j,nz))))
1269 
1270  ! Distribute the work over a BBL of depth 20^2 ustar^2 / g' following
1271  ! Killworth & Edwards (1999) and Zilitikevich & Mironov (1996).
1272  ! Rho_top is determined by finding the density where
1273  ! integral(bottom, Z) (rho(z') - rho(Z)) dz' = rho_0 400 ustar^2 / g
1274 
1275  gh_sum_top(i) = r0_g * 400.0 * ustar_h**2
1276 
1277  do_i(i) = (g%mask2dT(i,j) > 0.5)
1278  htot(i) = gv%H_to_Z*h(i,j,nz)
1279  rho_htot(i) = gv%Rlay(nz)*(gv%H_to_Z*h(i,j,nz))
1280  rho_top(i) = gv%Rlay(1)
1281  if (cs%bulkmixedlayer .and. do_i(i)) rho_top(i) = gv%Rlay(kb(i)-1)
1282  enddo
1283 
1284  do k=nz-1,2,-1 ; domore = .false.
1285  do i=is,ie ; if (do_i(i)) then
1286  htot(i) = htot(i) + gv%H_to_Z*h(i,j,k)
1287  rho_htot(i) = rho_htot(i) + gv%Rlay(k)*(gv%H_to_Z*h(i,j,k))
1288  if (htot(i)*gv%Rlay(k-1) <= (rho_htot(i) - gh_sum_top(i))) then
1289  ! The top of the mixing is in the interface atop the current layer.
1290  rho_top(i) = (rho_htot(i) - gh_sum_top(i)) / htot(i)
1291  do_i(i) = .false.
1292  elseif (k <= kb(i)) then ; do_i(i) = .false.
1293  else ; domore = .true. ; endif
1294  endif ; enddo
1295  if (.not.domore) exit
1296  enddo ! k-loop
1297 
1298  do i=is,ie ; do_i(i) = (g%mask2dT(i,j) > 0.5) ; enddo
1299  do k=nz-1,kb_min,-1
1300  i_rem = 0
1301  do i=is,ie ; if (do_i(i)) then
1302  if (k<kb(i)) then ; do_i(i) = .false. ; cycle ; endif
1303  i_rem = i_rem + 1 ! Count the i-rows that are still being worked on.
1304  ! Apply vertical decay of the turbulent energy. This energy is
1305  ! simply lost.
1306  tke(i) = tke(i) * exp(-i2decay(i) * (gv%H_to_Z*(h(i,j,k) + h(i,j,k+1))))
1307 
1308 ! if (maxEnt(i,k) <= 0.0) cycle
1309  if (maxtke(i,k) <= 0.0) cycle
1310 
1311  ! This is an analytic integral where diffusity is a quadratic function of
1312  ! rho that goes asymptotically to 0 at Rho_top (vaguely following KPP?).
1313  if (tke(i) > 0.0) then
1314  if (rint(k) <= rho_top(i)) then
1315  tke_to_layer = tke(i)
1316  else
1317  drl = rint(k+1) - rint(k) ; drbot = rint(k+1) - rho_top(i)
1318  tke_to_layer = tke(i) * drl * &
1319  (3.0*drbot*(rint(k) - rho_top(i)) + drl**2) / (drbot**3)
1320  endif
1321  else ; tke_to_layer = 0.0 ; endif
1322 
1323  ! TKE_Ray has been initialized to 0 above.
1324  if (rayleigh_drag) tke_ray = 0.5*cs%BBL_effic * us%L_to_Z**2 * g%IareaT(i,j) * &
1325  ((g%areaCu(i-1,j) * visc%Ray_u(i-1,j,k) * u(i-1,j,k)**2 + &
1326  g%areaCu(i,j) * visc%Ray_u(i,j,k) * u(i,j,k)**2) + &
1327  (g%areaCv(i,j-1) * visc%Ray_v(i,j-1,k) * v(i,j-1,k)**2 + &
1328  g%areaCv(i,j) * visc%Ray_v(i,j,k) * v(i,j,k)**2))
1329 
1330  if (tke_to_layer + tke_ray > 0.0) then
1331  if (cs%BBL_mixing_as_max) then
1332  if (tke_to_layer + tke_ray > maxtke(i,k)) &
1333  tke_to_layer = maxtke(i,k) - tke_ray
1334 
1335  tke(i) = tke(i) - tke_to_layer
1336 
1337  if (kd_lay(i,k) < (tke_to_layer + tke_ray) * tke_to_kd(i,k)) then
1338  delta_kd = (tke_to_layer + tke_ray) * tke_to_kd(i,k) - kd_lay(i,k)
1339  if ((cs%Kd_max >= 0.0) .and. (delta_kd > cs%Kd_max)) then
1340  delta_kd = cs%Kd_max
1341  kd_lay(i,k) = kd_lay(i,k) + delta_kd
1342  else
1343  kd_lay(i,k) = (tke_to_layer + tke_ray) * tke_to_kd(i,k)
1344  endif
1345  if (present(kd_int)) then
1346  kd_int(i,k) = kd_int(i,k) + 0.5 * delta_kd
1347  kd_int(i,k+1) = kd_int(i,k+1) + 0.5 * delta_kd
1348  endif
1349  if (do_diag_kd_bbl) then
1350  kd_bbl(i,j,k) = kd_bbl(i,j,k) + 0.5 * delta_kd
1351  kd_bbl(i,j,k+1) = kd_bbl(i,j,k+1) + 0.5 * delta_kd
1352  endif
1353  endif
1354  else
1355  if (kd_lay(i,k) >= maxtke(i,k) * tke_to_kd(i,k)) then
1356  tke_here = 0.0
1357  tke(i) = tke(i) + tke_ray
1358  elseif (kd_lay(i,k) + (tke_to_layer + tke_ray) * tke_to_kd(i,k) > &
1359  maxtke(i,k) * tke_to_kd(i,k)) then
1360  tke_here = ((tke_to_layer + tke_ray) + kd_lay(i,k) / tke_to_kd(i,k)) - maxtke(i,k)
1361  tke(i) = (tke(i) - tke_here) + tke_ray
1362  else
1363  tke_here = tke_to_layer + tke_ray
1364  tke(i) = tke(i) - tke_to_layer
1365  endif
1366  if (tke(i) < 0.0) tke(i) = 0.0 ! This should be unnecessary?
1367 
1368  if (tke_here > 0.0) then
1369  delta_kd = tke_here * tke_to_kd(i,k)
1370  if (cs%Kd_max >= 0.0) delta_kd = min(delta_kd, cs%Kd_max)
1371  kd_lay(i,k) = kd_lay(i,k) + delta_kd
1372  if (present(kd_int)) then
1373  kd_int(i,k) = kd_int(i,k) + 0.5 * delta_kd
1374  kd_int(i,k+1) = kd_int(i,k+1) + 0.5 * delta_kd
1375  endif
1376  if (do_diag_kd_bbl) then
1377  kd_bbl(i,j,k) = kd_bbl(i,j,k) + 0.5 * delta_kd
1378  kd_bbl(i,j,k+1) = kd_bbl(i,j,k+1) + 0.5 * delta_kd
1379  endif
1380  endif
1381  endif
1382  endif
1383 
1384  ! This may be risky - in the case that there are exactly zero
1385  ! velocities at 4 neighboring points, but nonzero velocities
1386  ! above the iterations would stop too soon. I don't see how this
1387  ! could happen in practice. RWH
1388  if ((tke(i)<= 0.0) .and. (tke_ray == 0.0)) then
1389  do_i(i) = .false. ; i_rem = i_rem - 1
1390  endif
1391 
1392  endif ; enddo
1393  if (i_rem == 0) exit
1394  enddo ! k-loop
1395 

◆ 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(:,:,:), pointer  Kd_BBL,
real, dimension( g %isd: g %ied, g %ke), intent(inout), optional  Kd_lay,
real, dimension( g %isd: g %ied, g %ke+1), intent(inout), optional  Kd_int 
)
private

Calculates a BBL diffusivity use a Prandtl number 1 diffusivity 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
kd_bblInterface BBL diffusivity [Z2 T-1 ~> m2 s-1]
[in,out]kd_layLayer net diffusivity [Z2 T-1 ~> m2 s-1]
[in,out]kd_intInterface net diffusivity [Z2 T-1 ~> m2 s-1]

Definition at line 1403 of file MOM_set_diffusivity.F90.

1403  type(ocean_grid_type), intent(in) :: g !< Grid structure
1404  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
1405  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1406  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1407  intent(in) :: u !< u component of flow [L T-1 ~> m s-1]
1408  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1409  intent(in) :: v !< v component of flow [L T-1 ~> m s-1]
1410  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1411  intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
1412  type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
1413  !! thermodynamic fields.
1414  type(forcing), intent(in) :: fluxes !< Surface fluxes structure
1415  type(vertvisc_type), intent(in) :: visc !< Structure containing vertical viscosities, bottom
1416  !! boundary layer properies, and related fields.
1417  integer, intent(in) :: j !< j-index of row to work on
1418  real, dimension(SZI_(G),SZK_(G)+1), &
1419  intent(in) :: n2_int !< Square of Brunt-Vaisala at interfaces [T-2 ~> s-2]
1420  type(set_diffusivity_cs), pointer :: cs !< Diffusivity control structure
1421  real, dimension(:,:,:), pointer :: kd_bbl !< Interface BBL diffusivity [Z2 T-1 ~> m2 s-1]
1422  real, dimension(SZI_(G),SZK_(G)), &
1423  optional, intent(inout) :: kd_lay !< Layer net diffusivity [Z2 T-1 ~> m2 s-1]
1424  real, dimension(SZI_(G),SZK_(G)+1), &
1425  optional, intent(inout) :: kd_int !< Interface net diffusivity [Z2 T-1 ~> m2 s-1]
1426 
1427  ! Local variables
1428  real :: tke_column ! net TKE input into the column [Z3 T-3 ~> m3 s-3]
1429  real :: tke_remaining ! remaining TKE available for mixing in this layer and above [Z3 T-3 ~> m3 s-3]
1430  real :: tke_consumed ! TKE used for mixing in this layer [Z3 T-3 ~> m3 s-3]
1431  real :: tke_kd_wall ! TKE associated with unlimited law of the wall mixing [Z3 T-3 ~> m3 s-3]
1432  real :: cdrag_sqrt ! square root of the drag coefficient [nondim]
1433  real :: ustar ! value of ustar at a thickness point [Z T-1 ~> m s-1].
1434  real :: ustar2 ! square of ustar, for convenience [Z2 T-2 ~> m2 s-2]
1435  real :: absf ! average absolute value of Coriolis parameter around a thickness point [T-1 ~> s-1]
1436  real :: dh, dhm1 ! thickness of layers k and k-1, respecitvely [Z ~> m].
1437  real :: z_bot ! distance to interface k from bottom [Z ~> m].
1438  real :: d_minus_z ! distance to interface k from surface [Z ~> m].
1439  real :: total_thickness ! total thickness of water column [Z ~> m].
1440  real :: idecay ! inverse of decay scale used for "Joule heating" loss of TKE with height [Z-1 ~> m-1].
1441  real :: kd_wall ! Law of the wall diffusivity [Z2 T-1 ~> m2 s-1].
1442  real :: kd_lower ! diffusivity for lower interface [Z2 T-1 ~> m2 s-1]
1443  real :: ustar_d ! u* x D [Z2 T-1 ~> m2 s-1].
1444  real :: i_rho0 ! 1 / rho0 [R-1 ~> m3 kg-1]
1445  real :: n2_min ! Minimum value of N2 to use in calculation of TKE_Kd_wall [T-2 ~> s-2]
1446  logical :: rayleigh_drag ! Set to true if there are Rayleigh drag velocities defined in visc, on
1447  ! the assumption that this extracted energy also drives diapycnal mixing.
1448  integer :: i, k, km1
1449  real, parameter :: von_karm = 0.41 ! Von Karman constant (http://en.wikipedia.org/wiki/Von_Karman_constant)
1450  logical :: do_diag_kd_bbl
1451 
1452  if (.not.(cs%bottomdraglaw .and. (cs%BBL_effic>0.0))) return
1453  do_diag_kd_bbl = associated(kd_bbl)
1454 
1455  n2_min = 0.
1456  if (cs%LOTW_BBL_use_omega) n2_min = cs%omega**2
1457 
1458  ! Determine whether to add Rayleigh drag contribution to TKE
1459  rayleigh_drag = .false.
1460  if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) rayleigh_drag = .true.
1461  i_rho0 = 1.0 / (gv%Rho0)
1462  cdrag_sqrt = sqrt(cs%cdrag)
1463 
1464  do i=g%isc,g%iec ! Developed in single-column mode
1465 
1466  ! Column-wise parameters.
1467  absf = 0.25 * ((abs(g%CoriolisBu(i-1,j-1)) + abs(g%CoriolisBu(i,j))) + &
1468  (abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j-1)))) ! Non-zero on equator!
1469 
1470  ! u* at the bottom [Z T-1 ~> m s-1].
1471  ustar = visc%ustar_BBL(i,j)
1472  ustar2 = ustar**2
1473  ! In add_drag_diffusivity(), fluxes%ustar_tidal is added in. This might be double counting
1474  ! since ustar_BBL should already include all contributions to u*? -AJA
1475  !### Examine the question of whether there is double counting of fluxes%ustar_tidal.
1476  if (associated(fluxes%ustar_tidal)) ustar = ustar + fluxes%ustar_tidal(i,j)
1477 
1478  ! The maximum decay scale should be something of order 200 m. We use the smaller of u*/f and
1479  ! (IMax_decay)^-1 as the decay scale. If ustar = 0, this is land so this value doesn't matter.
1480  idecay = cs%IMax_decay
1481  if ((ustar > 0.0) .and. (absf > cs%IMax_decay * ustar)) idecay = absf / ustar
1482 
1483  ! Energy input at the bottom [Z3 T-3 ~> m3 s-3].
1484  ! (Note that visc%TKE_BBL is in [Z3 T-3 ~> m3 s-3], set in set_BBL_TKE().)
1485  ! I am still unsure about sqrt(cdrag) in this expressions - AJA
1486  tke_column = cdrag_sqrt * visc%TKE_BBL(i,j)
1487  ! Add in tidal dissipation energy at the bottom [R Z3 T-3 ~> m3 s-3].
1488  ! Note that TKE_tidal is in [R Z3 T-3 ~> W m-2].
1489  if (associated(fluxes%TKE_tidal)) &
1490  tke_column = tke_column + fluxes%TKE_tidal(i,j) * i_rho0
1491  tke_column = cs%BBL_effic * tke_column ! Only use a fraction of the mechanical dissipation for mixing.
1492 
1493  tke_remaining = tke_column
1494  total_thickness = ( sum(h(i,j,:)) + gv%H_subroundoff )* gv%H_to_Z ! Total column thickness [Z ~> m].
1495  ustar_d = ustar * total_thickness
1496  z_bot = 0.
1497  kd_lower = 0. ! Diffusivity on bottom boundary.
1498 
1499  ! Work upwards from the bottom, accumulating work used until it exceeds the available TKE input
1500  ! at the bottom.
1501  do k=g%ke,2,-1
1502  dh = gv%H_to_Z * h(i,j,k) ! Thickness of this level [Z ~> m].
1503  km1 = max(k-1, 1)
1504  dhm1 = gv%H_to_Z * h(i,j,km1) ! Thickness of level above [Z ~> m].
1505 
1506  ! Add in additional energy input from bottom-drag against slopes (sides)
1507  if (rayleigh_drag) tke_remaining = tke_remaining + &
1508  0.5*cs%BBL_effic * us%L_to_Z**2 * g%IareaT(i,j) * &
1509  ((g%areaCu(i-1,j) * visc%Ray_u(i-1,j,k) * u(i-1,j,k)**2 + &
1510  g%areaCu(i,j) * visc%Ray_u(i,j,k) * u(i,j,k)**2) + &
1511  (g%areaCv(i,j-1) * visc%Ray_v(i,j-1,k) * v(i,j-1,k)**2 + &
1512  g%areaCv(i,j) * visc%Ray_v(i,j,k) * v(i,j,k)**2))
1513 
1514  ! Exponentially decay TKE across the thickness of the layer.
1515  ! This is energy loss in addition to work done as mixing, apparently to Joule heating.
1516  tke_remaining = exp(-idecay*dh) * tke_remaining
1517 
1518  z_bot = z_bot + h(i,j,k)*gv%H_to_Z ! Distance between upper interface of layer and the bottom [Z ~> m].
1519  d_minus_z = max(total_thickness - z_bot, 0.) ! Thickness above layer [Z ~> m].
1520 
1521  ! Diffusivity using law of the wall, limited by rotation, at height z [Z2 T-1 ~> m2 s-1].
1522  ! This calculation is at the upper interface of the layer
1523  if ( ustar_d + absf * ( z_bot * d_minus_z ) == 0.) then
1524  kd_wall = 0.
1525  else
1526  kd_wall = ((von_karm * ustar2) * (z_bot * d_minus_z)) &
1527  / (ustar_d + absf * (z_bot * d_minus_z))
1528  endif
1529 
1530  ! TKE associated with Kd_wall [Z3 T-3 ~> m3 s-3].
1531  ! This calculation if for the volume spanning the interface.
1532  tke_kd_wall = kd_wall * 0.5 * (dh + dhm1) * max(n2_int(i,k), n2_min)
1533 
1534  ! Now bound Kd such that the associated TKE is no greater than available TKE for mixing.
1535  if (tke_kd_wall > 0.) then
1536  tke_consumed = min(tke_kd_wall, tke_remaining)
1537  kd_wall = (tke_consumed / tke_kd_wall) * kd_wall ! Scale Kd so that only TKE_consumed is used.
1538  else
1539  ! Either N2=0 or dh = 0.
1540  if (tke_remaining > 0.) then
1541  kd_wall = cs%Kd_max
1542  else
1543  kd_wall = 0.
1544  endif
1545  tke_consumed = 0.
1546  endif
1547 
1548  ! Now use up the appropriate about of TKE associated with the diffusivity chosen
1549  tke_remaining = tke_remaining - tke_consumed ! Note this will be non-negative
1550 
1551  ! Add this BBL diffusivity to the model net diffusivity.
1552  if (present(kd_int)) kd_int(i,k) = kd_int(i,k) + kd_wall
1553  if (present(kd_lay)) kd_lay(i,k) = kd_lay(i,k) + 0.5 * (kd_wall + kd_lower)
1554  kd_lower = kd_wall ! Store for next layer up.
1555  if (do_diag_kd_bbl) kd_bbl(i,j,k) = kd_wall
1556  enddo ! k
1557  enddo ! i
1558 

◆ 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),szk_(g)), intent(in)  TKE_to_Kd,
real, dimension(szi_(g),szk_(g)), intent(inout), optional  Kd_lay,
real, dimension(szi_(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
[in]jThe j-index to work on
csDiffusivity control structure
[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_layThe diapycnal diffusivity in layers [Z2 T-1 ~> m2 s-1].
[in,out]kd_intThe diapycnal diffusivity at interfaces

Definition at line 1563 of file MOM_set_diffusivity.F90.

1563  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
1564  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
1565  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1566  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1567  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1568  type(forcing), intent(in) :: fluxes !< Surface fluxes structure
1569  integer, intent(in) :: j !< The j-index to work on
1570  type(set_diffusivity_cs), pointer :: cs !< Diffusivity control structure
1571  real, dimension(SZI_(G),SZK_(G)), intent(in) :: tke_to_kd !< The conversion rate between the TKE
1572  !! TKE dissipated within a layer and the
1573  !! diapycnal diffusivity witin that layer,
1574  !! usually (~Rho_0 / (G_Earth * dRho_lay))
1575  !! [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1]
1576  real, dimension(SZI_(G),SZK_(G)), &
1577  optional, intent(inout) :: kd_lay !< The diapycnal diffusivity in layers [Z2 T-1 ~> m2 s-1].
1578  real, dimension(SZI_(G),SZK_(G)+1), &
1579  optional, intent(inout) :: kd_int !< The diapycnal diffusivity at interfaces
1580  !! [Z2 T-1 ~> m2 s-1].
1581 
1582 ! This routine adds effects of mixed layer radiation to the layer diffusivities.
1583 
1584  real, dimension(SZI_(G)) :: h_ml ! Mixed layer thickness [Z ~> m].
1585  real, dimension(SZI_(G)) :: tke_ml_flux ! Mixed layer TKE flux [Z3 T-3 ~> m3 s-3]
1586  real, dimension(SZI_(G)) :: i_decay ! A decay rate [Z-1 ~> m-1].
1587  real, dimension(SZI_(G)) :: kd_mlr_ml ! Diffusivities associated with mixed layer radiation [Z2 T-1 ~> m2 s-1].
1588 
1589  real :: f_sq ! The square of the local Coriolis parameter or a related variable [T-2 ~> s-2].
1590  real :: h_ml_sq ! The square of the mixed layer thickness [Z2 ~> m2].
1591  real :: ustar_sq ! ustar squared [Z2 T-2 ~> m2 s-2]
1592  real :: kd_mlr ! A diffusivity associated with mixed layer turbulence radiation [Z2 T-1 ~> m2 s-1].
1593  real :: c1_6 ! 1/6
1594  real :: omega2 ! rotation rate squared [T-2 ~> s-2].
1595  real :: z1 ! layer thickness times I_decay [nondim]
1596  real :: dzl ! thickness converted to heights [Z ~> m].
1597  real :: i_decay_len2_tke ! squared inverse decay lengthscale for
1598  ! TKE, as used in the mixed layer code [Z-2 ~> m-2].
1599  real :: h_neglect ! negligibly small thickness [Z ~> m].
1600 
1601  logical :: do_any, do_i(szi_(g))
1602  integer :: i, k, is, ie, nz, kml
1603  is = g%isc ; ie = g%iec ; nz = g%ke
1604 
1605  omega2 = cs%omega**2
1606  c1_6 = 1.0 / 6.0
1607  kml = gv%nkml
1608  h_neglect = gv%H_subroundoff*gv%H_to_Z
1609 
1610  if (.not.cs%ML_radiation) return
1611 
1612  do i=is,ie ; h_ml(i) = 0.0 ; do_i(i) = (g%mask2dT(i,j) > 0.5) ; enddo
1613  do k=1,kml ; do i=is,ie ; h_ml(i) = h_ml(i) + gv%H_to_Z*h(i,j,k) ; enddo ; enddo
1614 
1615  do i=is,ie ; if (do_i(i)) then
1616  if (cs%ML_omega_frac >= 1.0) then
1617  f_sq = 4.0 * omega2
1618  else
1619  f_sq = 0.25 * ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
1620  (g%CoriolisBu(i,j-1)**2 + g%CoriolisBu(i-1,j)**2))
1621  if (cs%ML_omega_frac > 0.0) &
1622  f_sq = cs%ML_omega_frac * 4.0 * omega2 + (1.0 - cs%ML_omega_frac) * f_sq
1623  endif
1624 
1625  ustar_sq = max(fluxes%ustar(i,j), cs%ustar_min)**2
1626 
1627  tke_ml_flux(i) = (cs%mstar * cs%ML_rad_coeff) * (ustar_sq * (fluxes%ustar(i,j)))
1628  i_decay_len2_tke = cs%TKE_decay**2 * (f_sq / ustar_sq)
1629 
1630  if (cs%ML_rad_TKE_decay) &
1631  tke_ml_flux(i) = tke_ml_flux(i) * exp(-h_ml(i) * sqrt(i_decay_len2_tke))
1632 
1633  ! Calculate the inverse decay scale
1634  h_ml_sq = (cs%ML_rad_efold_coeff * (h_ml(i)+h_neglect))**2
1635  i_decay(i) = sqrt((i_decay_len2_tke * h_ml_sq + 1.0) / h_ml_sq)
1636 
1637  ! Average the dissipation layer kml+1, using
1638  ! a more accurate Taylor series approximations for very thin layers.
1639  z1 = (gv%H_to_Z*h(i,j,kml+1)) * i_decay(i)
1640  if (z1 > 1e-5) then
1641  kd_mlr = tke_ml_flux(i) * tke_to_kd(i,kml+1) * (1.0 - exp(-z1))
1642  else
1643  kd_mlr = tke_ml_flux(i) * tke_to_kd(i,kml+1) * (z1 * (1.0 - z1 * (0.5 - c1_6 * z1)))
1644  endif
1645  kd_mlr_ml(i) = min(kd_mlr, cs%ML_rad_kd_max)
1646  tke_ml_flux(i) = tke_ml_flux(i) * exp(-z1)
1647  endif ; enddo
1648 
1649  if (present(kd_lay)) then
1650  do k=1,kml+1 ; do i=is,ie ; if (do_i(i)) then
1651  kd_lay(i,k) = kd_lay(i,k) + kd_mlr_ml(i)
1652  endif ; enddo ; enddo
1653  endif
1654  if (present(kd_int)) then
1655  do k=2,kml+1 ; do i=is,ie ; if (do_i(i)) then
1656  kd_int(i,k) = kd_int(i,k) + kd_mlr_ml(i)
1657  endif ; enddo ; enddo
1658  if (kml<=nz-1) then ; do i=is,ie ; if (do_i(i)) then
1659  kd_int(i,kml+2) = kd_int(i,kml+2) + 0.5 * kd_mlr_ml(i)
1660  endif ; enddo ; endif
1661  endif
1662 
1663  do k=kml+2,nz-1
1664  do_any = .false.
1665  do i=is,ie ; if (do_i(i)) then
1666  dzl = gv%H_to_Z*h(i,j,k) ; z1 = dzl*i_decay(i)
1667  if (cs%ML_Rad_bug) then
1668  ! These expresssions are dimensionally inconsistent. -RWH
1669  ! This is supposed to be the integrated energy deposited in the layer,
1670  ! not the average over the layer as in these expressions.
1671  if (z1 > 1e-5) then
1672  kd_mlr = (tke_ml_flux(i) * tke_to_kd(i,k)) * & ! Units of Z2 T-1
1673  us%m_to_Z * ((1.0 - exp(-z1)) / dzl) ! Units of m-1
1674  else
1675  kd_mlr = (tke_ml_flux(i) * tke_to_kd(i,k)) * & ! Units of Z2 T-1
1676  us%m_to_Z * (i_decay(i) * (1.0 - z1 * (0.5 - c1_6*z1))) ! Units of m-1
1677  endif
1678  else
1679  if (z1 > 1e-5) then
1680  kd_mlr = (tke_ml_flux(i) * tke_to_kd(i,k)) * (1.0 - exp(-z1))
1681  else
1682  kd_mlr = (tke_ml_flux(i) * tke_to_kd(i,k)) * (z1 * (1.0 - z1 * (0.5 - c1_6*z1)))
1683  endif
1684  endif
1685  kd_mlr = min(kd_mlr, cs%ML_rad_kd_max)
1686  if (present(kd_lay)) then
1687  kd_lay(i,k) = kd_lay(i,k) + kd_mlr
1688  endif
1689  if (present(kd_int)) then
1690  kd_int(i,k) = kd_int(i,k) + 0.5 * kd_mlr
1691  kd_int(i,k+1) = kd_int(i,k+1) + 0.5 * kd_mlr
1692  endif
1693 
1694  tke_ml_flux(i) = tke_ml_flux(i) * exp(-z1)
1695  if (tke_ml_flux(i) * i_decay(i) < 0.1 * cs%Kd_min * omega2) then
1696  do_i(i) = .false.
1697  else ; do_any = .true. ; endif
1698  endif ; enddo
1699  if (.not.do_any) exit
1700  enddo
1701 

◆ 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 1077 of file MOM_set_diffusivity.F90.

1077  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
1078  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
1079  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1080  type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
1081  !! thermodynamic fields; absent fields have NULL ptrs.
1082  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1083  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
1084  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1085  intent(in) :: t_f !< layer temperatures with the values in massless layers
1086  !! filled vertically by diffusion [degC].
1087  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1088  intent(in) :: s_f !< Layer salinities with values in massless
1089  !! layers filled vertically by diffusion [ppt].
1090  integer, intent(in) :: j !< Meridional index upon which to work.
1091  type(set_diffusivity_cs), pointer :: cs !< Module control structure.
1092  real, dimension(SZI_(G),SZK_(G)+1), &
1093  intent(out) :: kd_t_dd !< Interface double diffusion diapycnal
1094  !! diffusivity for temp [Z2 T-1 ~> m2 s-1].
1095  real, dimension(SZI_(G),SZK_(G)+1), &
1096  intent(out) :: kd_s_dd !< Interface double diffusion diapycnal
1097  !! diffusivity for saln [Z2 T-1 ~> m2 s-1].
1098 
1099  real, dimension(SZI_(G)) :: &
1100  drho_dt, & ! partial derivatives of density wrt temp [R degC-1 ~> kg m-3 degC-1]
1101  drho_ds, & ! partial derivatives of density wrt saln [R ppt-1 ~> kg m-3 ppt-1]
1102  pres, & ! pressure at each interface [R L2 T-2 ~> Pa]
1103  temp_int, & ! temperature at interfaces [degC]
1104  salin_int ! Salinity at interfaces [ppt]
1105 
1106  real :: alpha_dt ! density difference between layers due to temp diffs [R ~> kg m-3]
1107  real :: beta_ds ! density difference between layers due to saln diffs [R ~> kg m-3]
1108 
1109  real :: rrho ! vertical density ratio [nondim]
1110  real :: diff_dd ! factor for double-diffusion [nondim]
1111  real :: kd_dd ! The dominant double diffusive diffusivity [Z2 T-1 ~> m2 s-1]
1112  real :: prandtl ! flux ratio for diffusive convection regime
1113 
1114  real, parameter :: rrho0 = 1.9 ! limit for double-diffusive density ratio [nondim]
1115 
1116  integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
1117  integer :: i, k, is, ie, nz
1118  is = g%isc ; ie = g%iec ; nz = g%ke
1119 
1120  if (associated(tv%eqn_of_state)) then
1121  do i=is,ie
1122  pres(i) = 0.0 ; kd_t_dd(i,1) = 0.0 ; kd_s_dd(i,1) = 0.0
1123  kd_t_dd(i,nz+1) = 0.0 ; kd_s_dd(i,nz+1) = 0.0
1124  enddo
1125  if (associated(tv%p_surf)) then ; do i=is,ie ; pres(i) = tv%p_surf(i,j) ; enddo ; endif
1126  eosdom(:) = eos_domain(g%HI)
1127  do k=2,nz
1128  do i=is,ie
1129  pres(i) = pres(i) + (gv%g_Earth*gv%H_to_RZ)*h(i,j,k-1)
1130  temp_int(i) = 0.5 * (t_f(i,j,k-1) + t_f(i,j,k))
1131  salin_int(i) = 0.5 * (s_f(i,j,k-1) + s_f(i,j,k))
1132  enddo
1133  call calculate_density_derivs(temp_int, salin_int, pres, drho_dt, drho_ds, &
1134  tv%eqn_of_state, eosdom)
1135 
1136  do i=is,ie
1137  alpha_dt = -1.0*drho_dt(i) * (t_f(i,j,k-1) - t_f(i,j,k))
1138  beta_ds = drho_ds(i) * (s_f(i,j,k-1) - s_f(i,j,k))
1139 
1140  if ((alpha_dt > beta_ds) .and. (beta_ds > 0.0)) then ! salt finger case
1141  rrho = min(alpha_dt / beta_ds, rrho0)
1142  diff_dd = 1.0 - ((rrho-1.0)/(rrho0-1.0))
1143  kd_dd = cs%Max_salt_diff_salt_fingers * diff_dd*diff_dd*diff_dd
1144  kd_t_dd(i,k) = 0.7 * kd_dd
1145  kd_s_dd(i,k) = kd_dd
1146  elseif ((alpha_dt < 0.) .and. (beta_ds < 0.) .and. (alpha_dt > beta_ds)) then ! diffusive convection
1147  rrho = alpha_dt / beta_ds
1148  kd_dd = cs%Kv_molecular * 0.909 * exp(4.6 * exp(-0.54 * (1/rrho - 1)))
1149  prandtl = 0.15*rrho
1150  if (rrho > 0.5) prandtl = (1.85-0.85/rrho)*rrho
1151  kd_t_dd(i,k) = kd_dd
1152  kd_s_dd(i,k) = prandtl * kd_dd
1153  else
1154  kd_t_dd(i,k) = 0.0 ; kd_s_dd(i,k) = 0.0
1155  endif
1156  enddo
1157  enddo
1158  endif
1159 

◆ 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 906 of file MOM_set_diffusivity.F90.

906  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
907  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
908  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
909  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
910  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
911  type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
912  !! thermodynamic fields.
913  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
914  intent(in) :: t_f !< layer temperature with the values in massless layers
915  !! filled vertically by diffusion [degC].
916  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
917  intent(in) :: s_f !< Layer salinities with values in massless
918  !! layers filled vertically by diffusion [ppt].
919  type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes
920  integer, intent(in) :: j !< j-index of row to work on
921  type(set_diffusivity_cs), pointer :: cs !< Diffusivity control structure
922  real, dimension(SZI_(G),SZK_(G)+1), &
923  intent(out) :: drho_int !< Change in locally referenced potential density
924  !! across each interface [R ~> kg m-3].
925  real, dimension(SZI_(G),SZK_(G)+1), &
926  intent(out) :: n2_int !< The squared buoyancy frequency at the interfaces [T-2 ~> s-2].
927  real, dimension(SZI_(G),SZK_(G)), &
928  intent(out) :: n2_lay !< The squared buoyancy frequency of the layers [T-2 ~> s-2].
929  real, dimension(SZI_(G)), intent(out) :: n2_bot !< The near-bottom squared buoyancy frequency [T-2 ~> s-2].
930  ! Local variables
931  real, dimension(SZI_(G),SZK_(G)+1) :: &
932  drho_int_unfilt, & ! unfiltered density differences across interfaces [R ~> kg m-3]
933  drho_dt, & ! partial derivative of density wrt temp [R degC-1 ~> kg m-3 degC-1]
934  drho_ds ! partial derivative of density wrt saln [R ppt-1 ~> kg m-3 ppt-1]
935 
936  real, dimension(SZI_(G)) :: &
937  pres, & ! pressure at each interface [R L2 T-2 ~> Pa]
938  temp_int, & ! temperature at each interface [degC]
939  salin_int, & ! salinity at each interface [ppt]
940  drho_bot, & ! A density difference [R ~> kg m-3]
941  h_amp, & ! The topographic roughness amplitude [Z ~> m].
942  hb, & ! The thickness of the bottom layer [Z ~> m].
943  z_from_bot ! The hieght above the bottom [Z ~> m].
944 
945  real :: rml_base ! density of the deepest variable density layer
946  real :: dz_int ! thickness associated with an interface [Z ~> m].
947  real :: g_rho0 ! gravitation acceleration divided by Bouss reference density
948  ! times some unit conversion factors [Z T-2 R-1 ~> m4 s-2 kg-1].
949  real :: h_neglect ! negligibly small thickness, in the same units as h.
950 
951  logical :: do_i(szi_(g)), do_any
952  integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
953  integer :: i, k, is, ie, nz
954 
955  is = g%isc ; ie = g%iec ; nz = g%ke
956  g_rho0 = (us%L_to_Z**2 * gv%g_Earth) / (gv%Rho0)
957  h_neglect = gv%H_subroundoff
958 
959  ! Find the (limited) density jump across each interface.
960  do i=is,ie
961  drho_int(i,1) = 0.0 ; drho_int(i,nz+1) = 0.0
962  drho_int_unfilt(i,1) = 0.0 ; drho_int_unfilt(i,nz+1) = 0.0
963  enddo
964  if (associated(tv%eqn_of_state)) then
965  if (associated(fluxes%p_surf)) then
966  do i=is,ie ; pres(i) = fluxes%p_surf(i,j) ; enddo
967  else
968  do i=is,ie ; pres(i) = 0.0 ; enddo
969  endif
970  eosdom(:) = eos_domain(g%HI)
971  do k=2,nz
972  do i=is,ie
973  pres(i) = pres(i) + (gv%g_Earth*gv%H_to_RZ)*h(i,j,k-1)
974  temp_int(i) = 0.5 * (t_f(i,j,k) + t_f(i,j,k-1))
975  salin_int(i) = 0.5 * (s_f(i,j,k) + s_f(i,j,k-1))
976  enddo
977  call calculate_density_derivs(temp_int, salin_int, pres, drho_dt(:,k), drho_ds(:,k), &
978  tv%eqn_of_state, eosdom)
979  do i=is,ie
980  drho_int(i,k) = max(drho_dt(i,k)*(t_f(i,j,k) - t_f(i,j,k-1)) + &
981  drho_ds(i,k)*(s_f(i,j,k) - s_f(i,j,k-1)), 0.0)
982  drho_int_unfilt(i,k) = max(drho_dt(i,k)*(tv%T(i,j,k) - tv%T(i,j,k-1)) + &
983  drho_ds(i,k)*(tv%S(i,j,k) - tv%S(i,j,k-1)), 0.0)
984  enddo
985  enddo
986  else
987  do k=2,nz ; do i=is,ie
988  drho_int(i,k) = gv%Rlay(k) - gv%Rlay(k-1)
989  enddo ; enddo
990  endif
991 
992  ! Set the buoyancy frequencies.
993  do k=1,nz ; do i=is,ie
994  n2_lay(i,k) = g_rho0 * 0.5*(drho_int(i,k) + drho_int(i,k+1)) / &
995  (gv%H_to_Z*(h(i,j,k) + h_neglect))
996  enddo ; enddo
997  do i=is,ie ; n2_int(i,1) = 0.0 ; n2_int(i,nz+1) = 0.0 ; enddo
998  do k=2,nz ; do i=is,ie
999  n2_int(i,k) = g_rho0 * drho_int(i,k) / &
1000  (0.5*gv%H_to_Z*(h(i,j,k-1) + h(i,j,k) + h_neglect))
1001  enddo ; enddo
1002 
1003  ! Find the bottom boundary layer stratification, and use this in the deepest layers.
1004  do i=is,ie
1005  hb(i) = 0.0 ; drho_bot(i) = 0.0 ; h_amp(i) = 0.0
1006  z_from_bot(i) = 0.5*gv%H_to_Z*h(i,j,nz)
1007  do_i(i) = (g%mask2dT(i,j) > 0.5)
1008  enddo
1009  if (cs%use_tidal_mixing) call tidal_mixing_h_amp(h_amp, g, j, cs%tidal_mixing_CSp)
1010 
1011  do k=nz,2,-1
1012  do_any = .false.
1013  do i=is,ie ; if (do_i(i)) then
1014  dz_int = 0.5*gv%H_to_Z*(h(i,j,k) + h(i,j,k-1))
1015  z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above
1016 
1017  hb(i) = hb(i) + dz_int
1018  drho_bot(i) = drho_bot(i) + drho_int(i,k)
1019 
1020  if (z_from_bot(i) > h_amp(i)) then
1021  if (k>2) then
1022  ! Always include at least one full layer.
1023  hb(i) = hb(i) + 0.5*gv%H_to_Z*(h(i,j,k-1) + h(i,j,k-2))
1024  drho_bot(i) = drho_bot(i) + drho_int(i,k-1)
1025  endif
1026  do_i(i) = .false.
1027  else
1028  do_any = .true.
1029  endif
1030  endif ; enddo
1031  if (.not.do_any) exit
1032  enddo
1033 
1034  do i=is,ie
1035  if (hb(i) > 0.0) then
1036  n2_bot(i) = (g_rho0 * drho_bot(i)) / hb(i)
1037  else ; n2_bot(i) = 0.0 ; endif
1038  z_from_bot(i) = 0.5*gv%H_to_Z*h(i,j,nz)
1039  do_i(i) = (g%mask2dT(i,j) > 0.5)
1040  enddo
1041 
1042  do k=nz,2,-1
1043  do_any = .false.
1044  do i=is,ie ; if (do_i(i)) then
1045  dz_int = 0.5*gv%H_to_Z*(h(i,j,k) + h(i,j,k-1))
1046  z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above
1047 
1048  n2_int(i,k) = n2_bot(i)
1049  if (k>2) n2_lay(i,k-1) = n2_bot(i)
1050 
1051  if (z_from_bot(i) > h_amp(i)) then
1052  if (k>2) n2_int(i,k-1) = n2_bot(i)
1053  do_i(i) = .false.
1054  else
1055  do_any = .true.
1056  endif
1057  endif ; enddo
1058  if (.not.do_any) exit
1059  enddo
1060 
1061  if (associated(tv%eqn_of_state)) then
1062  do k=1,nz+1 ; do i=is,ie
1063  drho_int(i,k) = drho_int_unfilt(i,k)
1064  enddo ; enddo
1065  endif
1066 

◆ 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 692 of file MOM_set_diffusivity.F90.

692  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
693  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
694  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
695  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
696  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
697  type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
698  !! thermodynamic fields.
699  real, dimension(SZI_(G),SZK_(G)+1), intent(in) :: drho_int !< Change in locally referenced potential density
700  !! across each interface [R ~> kg m-3].
701  real, dimension(SZI_(G),SZK_(G)), intent(in) :: n2_lay !< The squared buoyancy frequency of the
702  !! layers [T-2 ~> s-2].
703  integer, intent(in) :: j !< j-index of row to work on
704  real, intent(in) :: dt !< Time increment [T ~> s].
705  type(set_diffusivity_cs), pointer :: cs !< Diffusivity control structure
706  real, dimension(SZI_(G),SZK_(G)), intent(out) :: tke_to_kd !< The conversion rate between the
707  !! TKE dissipated within a layer and the
708  !! diapycnal diffusivity witin that layer,
709  !! usually (~Rho_0 / (G_Earth * dRho_lay))
710  !! [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1]
711  real, dimension(SZI_(G),SZK_(G)), intent(out) :: maxtke !< The energy required to for a layer to entrain
712  !! to its maximum realizable thickness [Z3 T-3 ~> m3 s-3]
713  integer, dimension(SZI_(G)), intent(out) :: kb !< Index of lightest layer denser than the buffer
714  !! layer, or -1 without a bulk mixed layer.
715  ! Local variables
716  real, dimension(SZI_(G),SZK_(G)) :: &
717  ds_dsp1, & ! coordinate variable (sigma-2) difference across an
718  ! interface divided by the difference across the interface
719  ! below it [nondim]
720  dsp1_ds, & ! inverse coordinate variable (sigma-2) difference
721  ! across an interface times the difference across the
722  ! interface above it [nondim]
723  rho_0, & ! Layer potential densities relative to surface pressure [R ~> kg m-3]
724  maxent ! maxEnt is the maximum value of entrainment from below (with
725  ! compensating entrainment from above to keep the layer
726  ! density from changing) that will not deplete all of the
727  ! layers above or below a layer within a timestep [Z ~> m].
728  real, dimension(SZI_(G)) :: &
729  htot, & ! total thickness above or below a layer, or the
730  ! integrated thickness in the BBL [Z ~> m].
731  mfkb, & ! total thickness in the mixed and buffer layers times ds_dsp1 [Z ~> m].
732  p_ref, & ! array of tv%P_Ref pressures [R L2 T-2 ~> Pa]
733  rcv_kmb, & ! coordinate density in the lowest buffer layer [R ~> kg m-3]
734  p_0 ! An array of 0 pressures [R L2 T-2 ~> Pa]
735 
736  real :: dh_max ! maximum amount of entrainment a layer could
737  ! undergo before entraining all fluid in the layers
738  ! above or below [Z ~> m].
739  real :: drho_lay ! density change across a layer [R ~> kg m-3]
740  real :: omega2 ! rotation rate squared [T-2 ~> s-2]
741  real :: g_rho0 ! gravitation accel divided by Bouss ref density [Z T-2 R-1 ~> m4 s-2 kg-1]
742  real :: g_irho0 ! Alternate calculation of G_Rho0 for reproducibility [Z T-2 R-1 ~> m4 s-2 kg-1]
743  real :: i_rho0 ! inverse of Boussinesq reference density [R-1 ~> m3 kg-1]
744  real :: i_dt ! 1/dt [T-1 ~> s-1]
745  real :: h_neglect ! negligibly small thickness [H ~> m or kg m-2]
746  real :: hn2po2 ! h (N^2 + Omega^2), in [Z T-2 ~> m s-2].
747  logical :: do_i(szi_(g))
748 
749  integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
750  integer :: i, k, is, ie, nz, i_rem, kmb, kb_min
751  is = g%isc ; ie = g%iec ; nz = g%ke
752 
753  i_dt = 1.0 / dt
754  omega2 = cs%omega**2
755  h_neglect = gv%H_subroundoff
756  g_rho0 = (us%L_to_Z**2 * gv%g_Earth) / (gv%Rho0)
757  if (cs%answers_2018) then
758  i_rho0 = 1.0 / (gv%Rho0)
759  g_irho0 = (us%L_to_Z**2 * gv%g_Earth) * i_rho0
760  else
761  g_irho0 = g_rho0
762  endif
763 
764  ! Simple but coordinate-independent estimate of Kd/TKE
765  if (cs%simple_TKE_to_Kd) then
766  do k=1,nz ; do i=is,ie
767  hn2po2 = (gv%H_to_Z * h(i,j,k)) * (n2_lay(i,k) + omega2) ! Units of Z T-2.
768  if (hn2po2>0.) then
769  tke_to_kd(i,k) = 1.0 / hn2po2 ! Units of T2 Z-1.
770  else; tke_to_kd(i,k) = 0.; endif
771  ! The maximum TKE conversion we allow is really a statement
772  ! about the upper diffusivity we allow. Kd_max must be set.
773  maxtke(i,k) = hn2po2 * cs%Kd_max ! Units of Z3 T-3.
774  enddo ; enddo
775  kb(is:ie) = -1 ! kb should not be used by any code in non-layered mode -AJA
776  return
777  endif
778 
779  ! Determine kb - the index of the shallowest active interior layer.
780  if (cs%bulkmixedlayer) then
781  kmb = gv%nk_rho_varies
782  do i=is,ie ; p_0(i) = 0.0 ; p_ref(i) = tv%P_Ref ; enddo
783  eosdom(:) = eos_domain(g%HI)
784  do k=1,nz
785  call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_0, rho_0(:,k), tv%eqn_of_state, eosdom)
786  enddo
787  call calculate_density(tv%T(:,j,kmb), tv%S(:,j,kmb), p_ref, rcv_kmb, tv%eqn_of_state, eosdom)
788 
789  kb_min = kmb+1
790  do i=is,ie
791  ! Determine the next denser layer than the buffer layer in the
792  ! coordinate density (sigma-2).
793  do k=kmb+1,nz-1 ; if (rcv_kmb(i) <= gv%Rlay(k)) exit ; enddo
794  kb(i) = k
795 
796  ! Backtrack, in case there are massive layers above that are stable
797  ! in sigma-0.
798  do k=kb(i)-1,kmb+1,-1
799  if (rho_0(i,kmb) > rho_0(i,k)) exit
800  if (h(i,j,k)>2.0*gv%Angstrom_H) kb(i) = k
801  enddo
802  enddo
803 
804  call set_density_ratios(h, tv, kb, g, gv, us, cs, j, ds_dsp1, rho_0)
805  else ! not bulkmixedlayer
806  kb_min = 2 ; kmb = 0
807  do i=is,ie ; kb(i) = 1 ; enddo
808  call set_density_ratios(h, tv, kb, g, gv, us, cs, j, ds_dsp1)
809  endif
810 
811  ! Determine maxEnt - the maximum permitted entrainment from below by each
812  ! interior layer.
813  do k=2,nz-1 ; do i=is,ie
814  dsp1_ds(i,k) = 1.0 / ds_dsp1(i,k)
815  enddo ; enddo
816  do i=is,ie ; dsp1_ds(i,nz) = 0.0 ; enddo
817 
818  if (cs%bulkmixedlayer) then
819  kmb = gv%nk_rho_varies
820  do i=is,ie
821  htot(i) = gv%H_to_Z*h(i,j,kmb)
822  mfkb(i) = 0.0
823  if (kb(i) < nz) &
824  mfkb(i) = ds_dsp1(i,kb(i)) * (gv%H_to_Z*(h(i,j,kmb) - gv%Angstrom_H))
825  enddo
826  do k=1,kmb-1 ; do i=is,ie
827  htot(i) = htot(i) + gv%H_to_Z*h(i,j,k)
828  mfkb(i) = mfkb(i) + ds_dsp1(i,k+1)*(gv%H_to_Z*(h(i,j,k) - gv%Angstrom_H))
829  enddo ; enddo
830  else
831  do i=is,i
832  maxent(i,1) = 0.0 ; htot(i) = gv%H_to_Z*(h(i,j,1) - gv%Angstrom_H)
833  enddo
834  endif
835  do k=kb_min,nz-1 ; do i=is,ie
836  if (k == kb(i)) then
837  maxent(i,kb(i)) = mfkb(i)
838  elseif (k > kb(i)) then
839  if (cs%answers_2018) then
840  maxent(i,k) = (1.0/dsp1_ds(i,k))*(maxent(i,k-1) + htot(i))
841  else
842  maxent(i,k) = ds_dsp1(i,k)*(maxent(i,k-1) + htot(i))
843  endif
844  htot(i) = htot(i) + gv%H_to_Z*(h(i,j,k) - gv%Angstrom_H)
845  endif
846  enddo ; enddo
847 
848  do i=is,ie
849  htot(i) = gv%H_to_Z*(h(i,j,nz) - gv%Angstrom_H) ; maxent(i,nz) = 0.0
850  do_i(i) = (g%mask2dT(i,j) > 0.5)
851  enddo
852  do k=nz-1,kb_min,-1
853  i_rem = 0
854  do i=is,ie ; if (do_i(i)) then
855  if (k<kb(i)) then ; do_i(i) = .false. ; cycle ; endif
856  i_rem = i_rem + 1 ! Count the i-rows that are still being worked on.
857  maxent(i,k) = min(maxent(i,k),dsp1_ds(i,k+1)*maxent(i,k+1) + htot(i))
858  htot(i) = htot(i) + gv%H_to_Z*(h(i,j,k) - gv%Angstrom_H)
859  endif ; enddo
860  if (i_rem == 0) exit
861  enddo ! k-loop
862 
863  ! Now set maxTKE and TKE_to_Kd.
864  do i=is,ie
865  maxtke(i,1) = 0.0 ; tke_to_kd(i,1) = 0.0
866  maxtke(i,nz) = 0.0 ; tke_to_kd(i,nz) = 0.0
867  enddo
868  do k=2,kmb ; do i=is,ie
869  maxtke(i,k) = 0.0
870  tke_to_kd(i,k) = 1.0 / ((n2_lay(i,k) + omega2) * &
871  (gv%H_to_Z*(h(i,j,k) + h_neglect)))
872  enddo ; enddo
873  do k=kmb+1,kb_min-1 ; do i=is,ie
874  ! These are the properties in the deeper mixed and buffer layers, and
875  ! should perhaps be revisited.
876  maxtke(i,k) = 0.0 ; tke_to_kd(i,k) = 0.0
877  enddo ; enddo
878  do k=kb_min,nz-1 ; do i=is,ie
879  if (k<kb(i)) then
880  maxtke(i,k) = 0.0
881  tke_to_kd(i,k) = 0.0
882  else
883  ! maxTKE is found by determining the kappa that gives maxEnt.
884  ! kappa_max = I_dt * dRho_int(i,K+1) * maxEnt(i,k) * &
885  ! (GV%H_to_Z*h(i,j,k) + dh_max) / dRho_lay
886  ! maxTKE(i,k) = (GV%g_Earth*US%L_to_Z**2) * dRho_lay * kappa_max
887  ! dRho_int should already be non-negative, so the max is redundant?
888  dh_max = maxent(i,k) * (1.0 + dsp1_ds(i,k))
889  drho_lay = 0.5 * max(drho_int(i,k) + drho_int(i,k+1), 0.0)
890  maxtke(i,k) = i_dt * (g_irho0 * &
891  (0.5*max(drho_int(i,k+1) + dsp1_ds(i,k)*drho_int(i,k), 0.0))) * &
892  ((gv%H_to_Z*h(i,j,k) + dh_max) * maxent(i,k))
893  ! TKE_to_Kd should be rho_InSitu / G_Earth * (delta rho_InSitu)
894  ! The omega^2 term in TKE_to_Kd is due to a rescaling of the efficiency of turbulent
895  ! mixing by a factor of N^2 / (N^2 + Omega^2), as proposed by Melet et al., 2013?
896  tke_to_kd(i,k) = 1.0 / (g_rho0 * drho_lay + &
897  cs%omega**2 * gv%H_to_Z*(h(i,j,k) + h_neglect))
898  endif
899  enddo ; enddo
900 

◆ 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,
type(ocean_obc_type), optional, pointer  OBC 
)

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
obcOpen boundaries control structure.

Definition at line 1707 of file MOM_set_diffusivity.F90.

1707  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
1708  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
1709  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1710  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1711  intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1]
1712  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1713  intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]
1714  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1715  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1716  type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes
1717  type(vertvisc_type), intent(in) :: visc !< Structure containing vertical viscosities, bottom
1718  !! boundary layer properies, and related fields.
1719  type(set_diffusivity_cs), pointer :: cs !< Diffusivity control structure
1720  type(ocean_obc_type), optional, pointer :: obc !< Open boundaries control structure.
1721 
1722  ! This subroutine calculates several properties related to bottom
1723  ! boundary layer turbulence.
1724 
1725  real, dimension(SZI_(G)) :: &
1726  htot ! total thickness above or below a layer, or the
1727  ! integrated thickness in the BBL [Z ~> m].
1728 
1729  real, dimension(SZIB_(G)) :: &
1730  uhtot, & ! running integral of u in the BBL [Z L T-1 ~> m2 s-1]
1731  ustar, & ! bottom boundary layer turbulence speed [Z T-1 ~> m s-1].
1732  u2_bbl ! square of the mean zonal velocity in the BBL [L2 T-2 ~> m2 s-2]
1733 
1734  real :: vhtot(szi_(g)) ! running integral of v in the BBL [Z L T-1 ~> m2 s-1]
1735 
1736  real, dimension(SZI_(G),SZJB_(G)) :: &
1737  vstar, & ! ustar at at v-points [Z T-1 ~> m s-1].
1738  v2_bbl ! square of average meridional velocity in BBL [L2 T-2 ~> m2 s-2]
1739 
1740  real :: cdrag_sqrt ! square root of the drag coefficient [nondim]
1741  real :: hvel ! thickness at velocity points [Z ~> m].
1742 
1743  logical :: domore, do_i(szi_(g))
1744  integer :: i, j, k, is, ie, js, je, nz
1745  integer :: l_seg
1746  logical :: local_open_u_bc, local_open_v_bc
1747  logical :: has_obc
1748 
1749  local_open_u_bc = .false.
1750  local_open_v_bc = .false.
1751  if (present(obc)) then ; if (associated(obc)) then
1752  local_open_u_bc = obc%open_u_BCs_exist_globally
1753  local_open_v_bc = obc%open_v_BCs_exist_globally
1754  endif ; endif
1755 
1756  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1757 
1758  if (.not.associated(cs)) call mom_error(fatal,"set_BBL_TKE: "//&
1759  "Module must be initialized before it is used.")
1760 
1761  if (.not.cs%bottomdraglaw .or. (cs%BBL_effic<=0.0)) then
1762  if (associated(visc%ustar_BBL)) then
1763  do j=js,je ; do i=is,ie ; visc%ustar_BBL(i,j) = 0.0 ; enddo ; enddo
1764  endif
1765  if (associated(visc%TKE_BBL)) then
1766  do j=js,je ; do i=is,ie ; visc%TKE_BBL(i,j) = 0.0 ; enddo ; enddo
1767  endif
1768  return
1769  endif
1770 
1771  cdrag_sqrt = sqrt(cs%cdrag)
1772 
1773  !$OMP parallel default(shared) private(do_i,vhtot,htot,domore,hvel,uhtot,ustar,u2_bbl)
1774  !$OMP do
1775  do j=js-1,je
1776  ! Determine ustar and the square magnitude of the velocity in the
1777  ! bottom boundary layer. Together these give the TKE source and
1778  ! vertical decay scale.
1779  do i=is,ie ; if ((g%mask2dCv(i,j) > 0.5) .and. (cdrag_sqrt*visc%bbl_thick_v(i,j) > 0.0)) then
1780  do_i(i) = .true. ; vhtot(i) = 0.0 ; htot(i) = 0.0
1781  vstar(i,j) = visc%Kv_bbl_v(i,j) / (cdrag_sqrt*visc%bbl_thick_v(i,j))
1782  else
1783  do_i(i) = .false. ; vstar(i,j) = 0.0 ; htot(i) = 0.0
1784  endif ; enddo
1785  do k=nz,1,-1
1786  domore = .false.
1787  do i=is,ie ; if (do_i(i)) then
1788  ! Determine if grid point is an OBC
1789  has_obc = .false.
1790  if (local_open_v_bc) then
1791  l_seg = obc%segnum_v(i,j)
1792  if (l_seg /= obc_none) then
1793  has_obc = obc%segment(l_seg)%open
1794  endif
1795  endif
1796 
1797  ! Compute h based on OBC state
1798  if (has_obc) then
1799  if (obc%segment(l_seg)%direction == obc_direction_n) then
1800  hvel = gv%H_to_Z*h(i,j,k)
1801  else
1802  hvel = gv%H_to_Z*h(i,j+1,k)
1803  endif
1804  else
1805  hvel = 0.5*gv%H_to_Z*(h(i,j,k) + h(i,j+1,k))
1806  endif
1807 
1808  if ((htot(i) + hvel) >= visc%bbl_thick_v(i,j)) then
1809  vhtot(i) = vhtot(i) + (visc%bbl_thick_v(i,j) - htot(i))*v(i,j,k)
1810  htot(i) = visc%bbl_thick_v(i,j)
1811  do_i(i) = .false.
1812  else
1813  vhtot(i) = vhtot(i) + hvel*v(i,j,k)
1814  htot(i) = htot(i) + hvel
1815  domore = .true.
1816  endif
1817  endif ; enddo
1818  if (.not.domore) exit
1819  enddo
1820  do i=is,ie ; if ((g%mask2dCv(i,j) > 0.5) .and. (htot(i) > 0.0)) then
1821  v2_bbl(i,j) = (vhtot(i)*vhtot(i))/(htot(i)*htot(i))
1822  else
1823  v2_bbl(i,j) = 0.0
1824  endif ; enddo
1825  enddo
1826  !$OMP do
1827  do j=js,je
1828  do i=is-1,ie ; if ((g%mask2dCu(i,j) > 0.5) .and. (cdrag_sqrt*visc%bbl_thick_u(i,j) > 0.0)) then
1829  do_i(i) = .true. ; uhtot(i) = 0.0 ; htot(i) = 0.0
1830  ustar(i) = visc%Kv_bbl_u(i,j) / (cdrag_sqrt*visc%bbl_thick_u(i,j))
1831  else
1832  do_i(i) = .false. ; ustar(i) = 0.0 ; htot(i) = 0.0
1833  endif ; enddo
1834  do k=nz,1,-1 ; domore = .false.
1835  do i=is-1,ie ; if (do_i(i)) then
1836  ! Determine if grid point is an OBC
1837  has_obc = .false.
1838  if (local_open_u_bc) then
1839  l_seg = obc%segnum_u(i,j)
1840  if (l_seg /= obc_none) then
1841  has_obc = obc%segment(l_seg)%open
1842  endif
1843  endif
1844 
1845  ! Compute h based on OBC state
1846  if (has_obc) then
1847  if (obc%segment(l_seg)%direction == obc_direction_e) then
1848  hvel = gv%H_to_Z*h(i,j,k)
1849  else ! OBC_DIRECTION_W
1850  hvel = gv%H_to_Z*h(i+1,j,k)
1851  endif
1852  else
1853  hvel = 0.5*gv%H_to_Z*(h(i,j,k) + h(i+1,j,k))
1854  endif
1855 
1856  if ((htot(i) + hvel) >= visc%bbl_thick_u(i,j)) then
1857  uhtot(i) = uhtot(i) + (visc%bbl_thick_u(i,j) - htot(i))*u(i,j,k)
1858  htot(i) = visc%bbl_thick_u(i,j)
1859  do_i(i) = .false.
1860  else
1861  uhtot(i) = uhtot(i) + hvel*u(i,j,k)
1862  htot(i) = htot(i) + hvel
1863  domore = .true.
1864  endif
1865  endif ; enddo
1866  if (.not.domore) exit
1867  enddo
1868  do i=is-1,ie ; if ((g%mask2dCu(i,j) > 0.5) .and. (htot(i) > 0.0)) then
1869  u2_bbl(i) = (uhtot(i)*uhtot(i))/(htot(i)*htot(i))
1870  else
1871  u2_bbl(i) = 0.0
1872  endif ; enddo
1873 
1874  do i=is,ie
1875  visc%ustar_BBL(i,j) = sqrt(0.5*g%IareaT(i,j) * &
1876  ((g%areaCu(i-1,j)*(ustar(i-1)*ustar(i-1)) + &
1877  g%areaCu(i,j)*(ustar(i)*ustar(i))) + &
1878  (g%areaCv(i,j-1)*(vstar(i,j-1)*vstar(i,j-1)) + &
1879  g%areaCv(i,j)*(vstar(i,j)*vstar(i,j))) ) )
1880  visc%TKE_BBL(i,j) = us%L_to_Z**2 * &
1881  (((g%areaCu(i-1,j)*(ustar(i-1)*u2_bbl(i-1)) + &
1882  g%areaCu(i,j) * (ustar(i)*u2_bbl(i))) + &
1883  (g%areaCv(i,j-1)*(vstar(i,j-1)*v2_bbl(i,j-1)) + &
1884  g%areaCv(i,j) * (vstar(i,j)*v2_bbl(i,j))) )*g%IareaT(i,j))
1885  enddo
1886  enddo
1887  !$OMP end parallel
1888 

◆ 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 1892 of file MOM_set_diffusivity.F90.

1892  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
1893  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
1894  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1895  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
1896  type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any
1897  !! available thermodynamic fields; absent
1898  !! fields have NULL ptrs.
1899  integer, dimension(SZI_(G)), intent(in) :: kb !< Index of lightest layer denser than the buffer
1900  !! layer, or -1 without a bulk mixed layer.
1901  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1902  type(set_diffusivity_cs), pointer :: cs !< Control structure returned by previous
1903  !! call to diabatic_entrain_init.
1904  integer, intent(in) :: j !< Meridional index upon which to work.
1905  real, dimension(SZI_(G),SZK_(G)), intent(out) :: ds_dsp1 !< Coordinate variable (sigma-2)
1906  !! difference across an interface divided by
1907  !! the difference across the interface below
1908  !! it [nondim]
1909  real, dimension(SZI_(G),SZK_(G)), &
1910  optional, intent(in) :: rho_0 !< Layer potential densities relative to
1911  !! surface press [R ~> kg m-3].
1912 
1913  ! Local variables
1914  real :: g_r0 ! g_R0 is a rescaled version of g/Rho [L2 Z-1 R-1 T-2 ~> m4 kg-1 s-2]
1915  real :: eps, tmp ! nondimensional temproray variables
1916  real :: a(szk_(g)), a_0(szk_(g)) ! nondimensional temporary variables
1917  real :: p_ref(szi_(g)) ! an array of tv%P_Ref pressures [R L2 T-2 ~> Pa]
1918  real :: rcv(szi_(g),szk_(g)) ! coordinate density in the mixed and buffer layers [R ~> kg m-3]
1919  real :: i_drho ! temporary variable [R-1 ~> m3 kg-1]
1920 
1921  integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
1922  integer :: i, k, k3, is, ie, nz, kmb
1923  is = g%isc ; ie = g%iec ; nz = g%ke
1924 
1925  do k=2,nz-1
1926  if (gv%g_prime(k+1) /= 0.0) then
1927  do i=is,ie
1928  ds_dsp1(i,k) = gv%g_prime(k) / gv%g_prime(k+1)
1929  enddo
1930  else
1931  do i=is,ie
1932  ds_dsp1(i,k) = 1.
1933  enddo
1934  endif
1935  enddo
1936 
1937  if (cs%bulkmixedlayer) then
1938  g_r0 = gv%g_Earth / (gv%Rho0)
1939  kmb = gv%nk_rho_varies
1940  eps = 0.1
1941  do i=is,ie ; p_ref(i) = tv%P_Ref ; enddo
1942  eosdom(:) = eos_domain(g%HI)
1943  do k=1,kmb
1944  call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, rcv(:,k), tv%eqn_of_state, eosdom)
1945  enddo
1946  do i=is,ie
1947  if (kb(i) <= nz-1) then
1948 ! Set up appropriately limited ratios of the reduced gravities of the
1949 ! interfaces above and below the buffer layer and the next denser layer.
1950  k = kb(i)
1951 
1952  i_drho = g_r0 / gv%g_prime(k+1)
1953  ! The indexing convention for a is appropriate for the interfaces.
1954  do k3=1,kmb
1955  a(k3+1) = (gv%Rlay(k) - rcv(i,k3)) * i_drho
1956  enddo
1957  if ((present(rho_0)) .and. (a(kmb+1) < 2.0*eps*ds_dsp1(i,k))) then
1958 ! If the buffer layer nearly matches the density of the layer below in the
1959 ! coordinate variable (sigma-2), use the sigma-0-based density ratio if it is
1960 ! greater (and stable).
1961  if ((rho_0(i,k) > rho_0(i,kmb)) .and. &
1962  (rho_0(i,k+1) > rho_0(i,k))) then
1963  i_drho = 1.0 / (rho_0(i,k+1)-rho_0(i,k))
1964  a_0(kmb+1) = min((rho_0(i,k)-rho_0(i,kmb)) * i_drho, ds_dsp1(i,k))
1965  if (a_0(kmb+1) > a(kmb+1)) then
1966  do k3=2,kmb
1967  a_0(k3) = a_0(kmb+1) + (rho_0(i,kmb)-rho_0(i,k3-1)) * i_drho
1968  enddo
1969  if (a(kmb+1) <= eps*ds_dsp1(i,k)) then
1970  do k3=2,kmb+1 ; a(k3) = a_0(k3) ; enddo
1971  else
1972 ! Alternative... tmp = 0.5*(1.0 - cos(PI*(a(K2+1)/(eps*ds_dsp1(i,k)) - 1.0)) )
1973  tmp = a(kmb+1)/(eps*ds_dsp1(i,k)) - 1.0
1974  do k3=2,kmb+1 ; a(k3) = tmp*a(k3) + (1.0-tmp)*a_0(k3) ; enddo
1975  endif
1976  endif
1977  endif
1978  endif
1979 
1980  ds_dsp1(i,k) = max(a(kmb+1),1e-5)
1981 
1982  do k3=2,kmb
1983 ! ds_dsp1(i,k3) = MAX(a(k3),1e-5)
1984  ! Deliberately treat convective instabilies of the upper mixed
1985  ! and buffer layers with respect to the deepest buffer layer as
1986  ! though they don't exist. They will be eliminated by the upcoming
1987  ! call to the mixedlayer code anyway.
1988  ! The indexing convention is appropriate for the interfaces.
1989  ds_dsp1(i,k3) = max(a(k3),ds_dsp1(i,k))
1990  enddo
1991  endif ! (kb(i) <= nz-1)
1992  enddo ! I-loop.
1993  endif ! bulkmixedlayer
1994 

◆ 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), optional  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 213 of file MOM_set_diffusivity.F90.

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

◆ 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 2322 of file MOM_set_diffusivity.F90.

2322  type(set_diffusivity_cs), pointer :: cs !< Control structure for this module
2323 
2324  if (.not.associated(cs)) return
2325 
2326  call bkgnd_mixing_end(cs%bkgnd_mixing_csp)
2327 
2328  if (cs%use_tidal_mixing) call tidal_mixing_end(cs%tidal_mixing_CSp)
2329 
2330  if (cs%user_change_diff) call user_change_diff_end(cs%user_change_diff_CSp)
2331 
2332  if (cs%use_CVMix_shear) call cvmix_shear_end(cs%CVMix_shear_csp)
2333 
2334  if (cs%use_CVMix_ddiff) call cvmix_ddiff_end(cs%CVMix_ddiff_csp)
2335 
2336  if (associated(cs)) deallocate(cs)
2337 

◆ 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,
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_cspA pointer to the internal tides control structure
[out]halo_tsThe halo size of tracer points that must be valid for the calculations in set_diffusivity.

Definition at line 1998 of file MOM_set_diffusivity.F90.

1998  type(time_type), intent(in) :: time !< The current model time
1999  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
2000  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
2001  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
2002  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
2003  !! parameters.
2004  type(diag_ctrl), target, intent(inout) :: diag !< A structure used to regulate diagnostic output.
2005  type(set_diffusivity_cs), pointer :: cs !< pointer set to point to the module control
2006  !! structure.
2007  type(int_tide_cs), pointer :: int_tide_csp !< A pointer to the internal tides control
2008  !! structure
2009  integer, optional, intent(out) :: halo_ts !< The halo size of tracer points that must be
2010  !! valid for the calculations in set_diffusivity.
2011 
2012  ! Local variables
2013  real :: decay_length
2014  logical :: ml_use_omega
2015  logical :: default_2018_answers
2016  ! This include declares and sets the variable "version".
2017 # include "version_variable.h"
2018  character(len=40) :: mdl = "MOM_set_diffusivity" ! This module's name.
2019  real :: omega_frac_dflt
2020  logical :: bryan_lewis_diffusivity ! If true, the background diapycnal diffusivity uses
2021  ! the Bryan-Lewis (1979) style tanh profile.
2022  integer :: i, j, is, ie, js, je
2023  integer :: isd, ied, jsd, jed
2024 
2025  if (associated(cs)) then
2026  call mom_error(warning, "diabatic_entrain_init called with an associated "// &
2027  "control structure.")
2028  return
2029  endif
2030  allocate(cs)
2031 
2032  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
2033  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2034 
2035  cs%diag => diag
2036  if (associated(int_tide_csp)) cs%int_tide_CSp => int_tide_csp
2037 
2038  ! These default values always need to be set.
2039  cs%BBL_mixing_as_max = .true.
2040  cs%cdrag = 0.003 ; cs%BBL_effic = 0.0
2041  cs%bulkmixedlayer = (gv%nkml > 0)
2042 
2043  ! Read all relevant parameters and write them to the model log.
2044  call log_version(param_file, mdl, version, "")
2045 
2046  call get_param(param_file, mdl, "INPUTDIR", cs%inputdir, default=".")
2047  cs%inputdir = slasher(cs%inputdir)
2048  call get_param(param_file, mdl, "FLUX_RI_MAX", cs%FluxRi_max, &
2049  "The flux Richardson number where the stratification is "//&
2050  "large enough that N2 > omega2. The full expression for "//&
2051  "the Flux Richardson number is usually "//&
2052  "FLUX_RI_MAX*N2/(N2+OMEGA2).", units="nondim", default=0.2)
2053  call get_param(param_file, mdl, "OMEGA", cs%omega, &
2054  "The rotation rate of the earth.", units="s-1", default=7.2921e-5, scale=us%T_to_s)
2055 
2056  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
2057  "This sets the default value for the various _2018_ANSWERS parameters.", &
2058  default=.false.)
2059  call get_param(param_file, mdl, "SET_DIFF_2018_ANSWERS", cs%answers_2018, &
2060  "If true, use the order of arithmetic and expressions that recover the "//&
2061  "answers from the end of 2018. Otherwise, use updated and more robust "//&
2062  "forms of the same expressions.", default=default_2018_answers)
2063 
2064  ! CS%use_tidal_mixing is set to True if an internal tidal dissipation scheme is to be used.
2065  cs%use_tidal_mixing = tidal_mixing_init(time, g, gv, us, param_file, diag, &
2066  cs%tidal_mixing_CSp)
2067 
2068  call get_param(param_file, mdl, "ML_RADIATION", cs%ML_radiation, &
2069  "If true, allow a fraction of TKE available from wind "//&
2070  "work to penetrate below the base of the mixed layer "//&
2071  "with a vertical decay scale determined by the minimum "//&
2072  "of: (1) The depth of the mixed layer, (2) an Ekman "//&
2073  "length scale.", default=.false.)
2074  if (cs%ML_radiation) then
2075  ! This give a minimum decay scale that is typically much less than Angstrom.
2076  cs%ustar_min = 2e-4 * cs%omega * (gv%Angstrom_Z + gv%H_subroundoff*gv%H_to_Z)
2077 
2078  call get_param(param_file, mdl, "ML_RAD_EFOLD_COEFF", cs%ML_rad_efold_coeff, &
2079  "A coefficient that is used to scale the penetration "//&
2080  "depth for turbulence below the base of the mixed layer. "//&
2081  "This is only used if ML_RADIATION is true.", units="nondim", default=0.2)
2082  call get_param(param_file, mdl, "ML_RAD_BUG", cs%ML_rad_bug, &
2083  "If true use code with a bug that reduces the energy available "//&
2084  "in the transition layer by a factor of the inverse of the energy "//&
2085  "deposition lenthscale (in m).", default=.false.)
2086  call get_param(param_file, mdl, "ML_RAD_KD_MAX", cs%ML_rad_kd_max, &
2087  "The maximum diapycnal diffusivity due to turbulence "//&
2088  "radiated from the base of the mixed layer. "//&
2089  "This is only used if ML_RADIATION is true.", &
2090  units="m2 s-1", default=1.0e-3, scale=us%m2_s_to_Z2_T)
2091  call get_param(param_file, mdl, "ML_RAD_COEFF", cs%ML_rad_coeff, &
2092  "The coefficient which scales MSTAR*USTAR^3 to obtain "//&
2093  "the energy available for mixing below the base of the "//&
2094  "mixed layer. This is only used if ML_RADIATION is true.", &
2095  units="nondim", default=0.2)
2096  call get_param(param_file, mdl, "ML_RAD_APPLY_TKE_DECAY", cs%ML_rad_TKE_decay, &
2097  "If true, apply the same exponential decay to ML_rad as "//&
2098  "is applied to the other surface sources of TKE in the "//&
2099  "mixed layer code. This is only used if ML_RADIATION is true.", default=.true.)
2100  call get_param(param_file, mdl, "MSTAR", cs%mstar, &
2101  "The ratio of the friction velocity cubed to the TKE "//&
2102  "input to the mixed layer.", units="nondim", default=1.2)
2103  call get_param(param_file, mdl, "TKE_DECAY", cs%TKE_decay, &
2104  "The ratio of the natural Ekman depth to the TKE decay scale.", &
2105  units="nondim", default=2.5)
2106  call get_param(param_file, mdl, "ML_USE_OMEGA", ml_use_omega, &
2107  "If true, use the absolute rotation rate instead of the "//&
2108  "vertical component of rotation when setting the decay "//&
2109  "scale for turbulence.", default=.false., do_not_log=.true.)
2110  omega_frac_dflt = 0.0
2111  if (ml_use_omega) then
2112  call mom_error(warning, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
2113  omega_frac_dflt = 1.0
2114  endif
2115  call get_param(param_file, mdl, "ML_OMEGA_FRAC", cs%ML_omega_frac, &
2116  "When setting the decay scale for turbulence, use this "//&
2117  "fraction of the absolute rotation rate blended with the "//&
2118  "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
2119  units="nondim", default=omega_frac_dflt)
2120  endif
2121 
2122  call get_param(param_file, mdl, "BOTTOMDRAGLAW", cs%bottomdraglaw, &
2123  "If true, the bottom stress is calculated with a drag "//&
2124  "law of the form c_drag*|u|*u. The velocity magnitude "//&
2125  "may be an assumed value or it may be based on the actual "//&
2126  "velocity in the bottommost HBBL, depending on LINEAR_DRAG.", default=.true.)
2127  if (cs%bottomdraglaw) then
2128  call get_param(param_file, mdl, "CDRAG", cs%cdrag, &
2129  "The drag coefficient relating the magnitude of the "//&
2130  "velocity field to the bottom stress. CDRAG is only used "//&
2131  "if BOTTOMDRAGLAW is true.", units="nondim", default=0.003)
2132  call get_param(param_file, mdl, "BBL_EFFIC", cs%BBL_effic, &
2133  "The efficiency with which the energy extracted by "//&
2134  "bottom drag drives BBL diffusion. This is only "//&
2135  "used if BOTTOMDRAGLAW is true.", units="nondim", default=0.20)
2136  call get_param(param_file, mdl, "BBL_MIXING_MAX_DECAY", decay_length, &
2137  "The maximum decay scale for the BBL diffusion, or 0 to allow the mixing "//&
2138  "to penetrate as far as stratification and rotation permit. The default "//&
2139  "for now is 200 m. This is only used if BOTTOMDRAGLAW is true.", &
2140  units="m", default=200.0, scale=us%m_to_Z)
2141 
2142  cs%IMax_decay = 0.0
2143  if (decay_length > 0.0) cs%IMax_decay = 1.0/decay_length
2144  call get_param(param_file, mdl, "BBL_MIXING_AS_MAX", cs%BBL_mixing_as_max, &
2145  "If true, take the maximum of the diffusivity from the "//&
2146  "BBL mixing and the other diffusivities. Otherwise, "//&
2147  "diffusivity from the BBL_mixing is simply added.", &
2148  default=.true.)
2149  call get_param(param_file, mdl, "USE_LOTW_BBL_DIFFUSIVITY", cs%use_LOTW_BBL_diffusivity, &
2150  "If true, uses a simple, imprecise but non-coordinate dependent, model "//&
2151  "of BBL mixing diffusivity based on Law of the Wall. Otherwise, uses "//&
2152  "the original BBL scheme.", default=.false.)
2153  if (cs%use_LOTW_BBL_diffusivity) then
2154  call get_param(param_file, mdl, "LOTW_BBL_USE_OMEGA", cs%LOTW_BBL_use_omega, &
2155  "If true, use the maximum of Omega and N for the TKE to diffusion "//&
2156  "calculation. Otherwise, N is N.", default=.true.)
2157  endif
2158  else
2159  cs%use_LOTW_BBL_diffusivity = .false. ! This parameterization depends on a u* from viscous BBL
2160  endif
2161  cs%id_Kd_BBL = register_diag_field('ocean_model', 'Kd_BBL', diag%axesTi, time, &
2162  'Bottom Boundary Layer Diffusivity', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
2163  call get_param(param_file, mdl, "SIMPLE_TKE_TO_KD", cs%simple_TKE_to_Kd, &
2164  "If true, uses a simple estimate of Kd/TKE that will "//&
2165  "work for arbitrary vertical coordinates. If false, "//&
2166  "calculates Kd/TKE and bounds based on exact energetics "//&
2167  "for an isopycnal layer-formulation.", default=.false.)
2168 
2169  ! set params related to the background mixing
2170  call bkgnd_mixing_init(time, g, gv, us, param_file, cs%diag, cs%bkgnd_mixing_csp)
2171 
2172  call get_param(param_file, mdl, "KV", cs%Kv, &
2173  "The background kinematic viscosity in the interior. "//&
2174  "The molecular value, ~1e-6 m2 s-1, may be used.", &
2175  units="m2 s-1", scale=us%m2_s_to_Z2_T, fail_if_missing=.true.)
2176 
2177  call get_param(param_file, mdl, "KD", cs%Kd, &
2178  "The background diapycnal diffusivity of density in the "//&
2179  "interior. Zero or the molecular value, ~1e-7 m2 s-1, "//&
2180  "may be used.", default=0.0, units="m2 s-1", scale=us%m2_s_to_Z2_T)
2181  call get_param(param_file, mdl, "KD_MIN", cs%Kd_min, &
2182  "The minimum diapycnal diffusivity.", &
2183  units="m2 s-1", default=0.01*cs%Kd*us%Z2_T_to_m2_s, scale=us%m2_s_to_Z2_T)
2184  call get_param(param_file, mdl, "KD_MAX", cs%Kd_max, &
2185  "The maximum permitted increment for the diapycnal "//&
2186  "diffusivity from TKE-based parameterizations, or a negative "//&
2187  "value for no limit.", units="m2 s-1", default=-1.0, scale=us%m2_s_to_Z2_T)
2188  if (cs%simple_TKE_to_Kd .and. cs%Kd_max<=0.) call mom_error(fatal, &
2189  "set_diffusivity_init: To use SIMPLE_TKE_TO_KD, KD_MAX must be set to >0.")
2190  call get_param(param_file, mdl, "KD_ADD", cs%Kd_add, &
2191  "A uniform diapycnal diffusivity that is added "//&
2192  "everywhere without any filtering or scaling.", &
2193  units="m2 s-1", default=0.0, scale=us%m2_s_to_Z2_T)
2194  if (cs%use_LOTW_BBL_diffusivity .and. cs%Kd_max<=0.) call mom_error(fatal, &
2195  "set_diffusivity_init: KD_MAX must be set (positive) when "// &
2196  "USE_LOTW_BBL_DIFFUSIVITY=True.")
2197  call get_param(param_file, mdl, "KD_SMOOTH", cs%Kd_smooth, &
2198  "A diapycnal diffusivity that is used to interpolate "//&
2199  "more sensible values of T & S into thin layers.", &
2200  units="m2 s-1", default=1.0e-6, scale=us%m2_s_to_Z2_T)
2201 
2202  call get_param(param_file, mdl, "DEBUG", cs%debug, &
2203  "If true, write out verbose debugging data.", &
2204  default=.false., debuggingparam=.true.)
2205 
2206  call get_param(param_file, mdl, "USER_CHANGE_DIFFUSIVITY", cs%user_change_diff, &
2207  "If true, call user-defined code to change the diffusivity.", default=.false.)
2208 
2209  call get_param(param_file, mdl, "DISSIPATION_MIN", cs%dissip_min, &
2210  "The minimum dissipation by which to determine a lower "//&
2211  "bound of Kd (a floor).", &
2212  units="W m-3", default=0.0, scale=us%W_m2_to_RZ3_T3*us%Z_to_m)
2213  call get_param(param_file, mdl, "DISSIPATION_N0", cs%dissip_N0, &
2214  "The intercept when N=0 of the N-dependent expression "//&
2215  "used to set a minimum dissipation by which to determine "//&
2216  "a lower bound of Kd (a floor): A in eps_min = A + B*N.", &
2217  units="W m-3", default=0.0, scale=us%W_m2_to_RZ3_T3*us%Z_to_m)
2218  call get_param(param_file, mdl, "DISSIPATION_N1", cs%dissip_N1, &
2219  "The coefficient multiplying N, following Gargett, used to "//&
2220  "set a minimum dissipation by which to determine a lower "//&
2221  "bound of Kd (a floor): B in eps_min = A + B*N", &
2222  units="J m-3", default=0.0, scale=us%W_m2_to_RZ3_T3*us%Z_to_m*us%s_to_T)
2223  call get_param(param_file, mdl, "DISSIPATION_KD_MIN", cs%dissip_Kd_min, &
2224  "The minimum vertical diffusivity applied as a floor.", &
2225  units="m2 s-1", default=0.0, scale=us%m2_s_to_Z2_T)
2226 
2227  cs%limit_dissipation = (cs%dissip_min>0.) .or. (cs%dissip_N1>0.) .or. &
2228  (cs%dissip_N0>0.) .or. (cs%dissip_Kd_min>0.)
2229  cs%dissip_N2 = 0.0
2230  if (cs%FluxRi_max > 0.0) &
2231  cs%dissip_N2 = cs%dissip_Kd_min * gv%Rho0 / cs%FluxRi_max
2232 
2233  cs%id_Kd_bkgnd = register_diag_field('ocean_model', 'Kd_bkgnd', diag%axesTi, time, &
2234  'Background diffusivity added by MOM_bkgnd_mixing module', 'm2/s', conversion=us%Z2_T_to_m2_s)
2235  cs%id_Kv_bkgnd = register_diag_field('ocean_model', 'Kv_bkgnd', diag%axesTi, time, &
2236  'Background viscosity added by MOM_bkgnd_mixing module', 'm2/s', conversion=us%Z2_T_to_m2_s)
2237 
2238  cs%id_Kd_layer = register_diag_field('ocean_model', 'Kd_layer', diag%axesTL, time, &
2239  'Diapycnal diffusivity of layers (as set)', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
2240 
2241  if (cs%use_tidal_mixing) then
2242  cs%id_Kd_Work = register_diag_field('ocean_model', 'Kd_Work', diag%axesTL, time, &
2243  'Work done by Diapycnal Mixing', 'W m-2', conversion=us%RZ3_T3_to_W_m2)
2244  cs%id_maxTKE = register_diag_field('ocean_model', 'maxTKE', diag%axesTL, time, &
2245  'Maximum layer TKE', 'm3 s-3', conversion=(us%Z_to_m**3*us%s_to_T**3))
2246  cs%id_TKE_to_Kd = register_diag_field('ocean_model', 'TKE_to_Kd', diag%axesTL, time, &
2247  'Convert TKE to Kd', 's2 m', conversion=us%Z2_T_to_m2_s*(us%m_to_Z**3*us%T_to_s**3))
2248  cs%id_N2 = register_diag_field('ocean_model', 'N2', diag%axesTi, time, &
2249  'Buoyancy frequency squared', 's-2', conversion=us%s_to_T**2, cmor_field_name='obvfsq', &
2250  cmor_long_name='Square of seawater buoyancy frequency', &
2251  cmor_standard_name='square_of_brunt_vaisala_frequency_in_sea_water')
2252  endif
2253 
2254  if (cs%user_change_diff) &
2255  cs%id_Kd_user = register_diag_field('ocean_model', 'Kd_user', diag%axesTi, time, &
2256  'User-specified Extra Diffusivity', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
2257 
2258  call get_param(param_file, mdl, "DOUBLE_DIFFUSION", cs%double_diffusion, &
2259  "If true, increase diffusivites for temperature or salt "//&
2260  "based on double-diffusive parameterization from MOM4/KPP.", &
2261  default=.false.)
2262 
2263  if (cs%double_diffusion) then
2264  call get_param(param_file, mdl, "MAX_RRHO_SALT_FINGERS", cs%Max_Rrho_salt_fingers, &
2265  "Maximum density ratio for salt fingering regime.", &
2266  default=2.55, units="nondim")
2267  call get_param(param_file, mdl, "MAX_SALT_DIFF_SALT_FINGERS", cs%Max_salt_diff_salt_fingers, &
2268  "Maximum salt diffusivity for salt fingering regime.", &
2269  default=1.e-4, units="m2 s-1", scale=us%m2_s_to_Z2_T)
2270  call get_param(param_file, mdl, "KV_MOLECULAR", cs%Kv_molecular, &
2271  "Molecular viscosity for calculation of fluxes under double-diffusive "//&
2272  "convection.", default=1.5e-6, units="m2 s-1", scale=us%m2_s_to_Z2_T)
2273  ! The default molecular viscosity follows the CCSM4.0 and MOM4p1 defaults.
2274  endif ! old double-diffusion
2275 
2276  if (cs%user_change_diff) then
2277  call user_change_diff_init(time, g, gv, us, param_file, diag, cs%user_change_diff_CSp)
2278  endif
2279 
2280  call get_param(param_file, mdl, "BRYAN_LEWIS_DIFFUSIVITY", bryan_lewis_diffusivity, &
2281  "If true, use a Bryan & Lewis (JGR 1979) like tanh "//&
2282  "profile of background diapycnal diffusivity with depth. "//&
2283  "This is done via CVMix.", default=.false., do_not_log=.true.)
2284  if (cs%use_tidal_mixing .and. bryan_lewis_diffusivity) &
2285  call mom_error(fatal,"MOM_Set_Diffusivity: "// &
2286  "Bryan-Lewis and internal tidal dissipation are both enabled. Choose one.")
2287 
2288  cs%useKappaShear = kappa_shear_init(time, g, gv, us, param_file, cs%diag, cs%kappaShear_CSp)
2289  cs%Vertex_Shear = kappa_shear_at_vertex(param_file)
2290 
2291  if (cs%useKappaShear) &
2292  id_clock_kappashear = cpu_clock_id('(Ocean kappa_shear)', grain=clock_module)
2293 
2294  ! CVMix shear-driven mixing
2295  cs%use_CVMix_shear = cvmix_shear_init(time, g, gv, us, param_file, cs%diag, cs%CVMix_shear_csp)
2296 
2297  ! CVMix double diffusion mixing
2298  cs%use_CVMix_ddiff = cvmix_ddiff_init(time, g, gv, us, param_file, cs%diag, cs%CVMix_ddiff_csp)
2299  if (cs%use_CVMix_ddiff) &
2300  id_clock_cvmix_ddiff = cpu_clock_id('(Double diffusion via CVMix)', grain=clock_module)
2301 
2302  if (cs%double_diffusion .or. cs%use_CVMix_ddiff) then
2303  cs%id_KT_extra = register_diag_field('ocean_model', 'KT_extra', diag%axesTi, time, &
2304  'Double-diffusive diffusivity for temperature', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
2305  cs%id_KS_extra = register_diag_field('ocean_model', 'KS_extra', diag%axesTi, time, &
2306  'Double-diffusive diffusivity for salinity', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
2307  endif
2308  if (cs%use_CVMix_ddiff) then
2309  cs%id_R_rho = register_diag_field('ocean_model', 'R_rho', diag%axesTi, time, &
2310  'Double-diffusion density ratio', 'nondim')
2311  endif
2312 
2313  if (present(halo_ts)) then
2314  halo_ts = 0
2315  if (cs%Vertex_Shear) halo_ts = 1
2316  endif
2317