MOM6
Data Types
mom_set_diffusivity Module Reference

Detailed Description

Calculate vertical diffusivity from all mixing processes.

Data Types

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

Function/Subroutine Documentation

◆ add_drag_diffusivity()

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

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

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

Definition at line 1105 of file MOM_set_diffusivity.F90.

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

◆ add_lotw_bbl_diffusivity()

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

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

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

Definition at line 1341 of file MOM_set_diffusivity.F90.

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

◆ add_mlrad_diffusivity()

subroutine mom_set_diffusivity::add_mlrad_diffusivity ( real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
type(forcing), intent(in)  fluxes,
integer, intent(in)  j,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(set_diffusivity_cs), pointer  CS,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  Kd_lay,
real, dimension(szi_(g),szk_(g)), intent(in)  TKE_to_Kd,
real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(inout), optional  Kd_int 
)
private

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

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

Definition at line 1502 of file MOM_set_diffusivity.F90.

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

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

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

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

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

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

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

◆ set_bbl_tke()

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

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

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

Definition at line 1642 of file MOM_set_diffusivity.F90.

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

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

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

◆ 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_in_T,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(set_diffusivity_cs), pointer  CS,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(out)  Kd_lay,
real, dimension(szi_(g),szj_(g),szk_(g)+1), intent(out), optional  Kd_int 
)

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

  1. Shear-driven mixing: two options, Jackson et at. and KPP interior;
  2. Background mixing via CVMix (Bryan-Lewis profile) or the scheme described by Harrison & Hallberg, JPO 2008;
  3. Double-diffusion, old method and new method via CVMix;
  4. Tidal mixing: many options available, see MOM_tidal_mixing.F90; In addition, this subroutine has the option to set the interior vertical viscosity associated with processes 1,2 and 4 listed above, which is stored in viscKv_slow. Vertical viscosity due to shear-driven mixing is passed via viscKv_shear
    Parameters
    [in]gThe ocean's grid structure.
    [in]gvThe ocean's vertical grid structure.
    [in]usA dimensional unit scaling type
    [in]uThe zonal velocity [m s-1].
    [in]vThe meridional velocity [m s-1].
    [in]hLayer thicknesses [H ~> m or kg m-2].
    [in]u_hZonal velocity interpolated to h points [m s-1].
    [in]v_hMeridional velocity interpolated to h points [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]dt_in_tTime increment [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 207 of file MOM_set_diffusivity.F90.

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

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

2214  type(set_diffusivity_CS), pointer :: CS !< Control structure for this module
2215 
2216  if (.not.associated(cs)) return
2217 
2218  call bkgnd_mixing_end(cs%bkgnd_mixing_csp)
2219 
2220  if (cs%user_change_diff) &
2221  call user_change_diff_end(cs%user_change_diff_CSp)
2222 
2223  if (cs%use_CVMix_shear) &
2224  call cvmix_shear_end(cs%CVMix_shear_csp)
2225 
2226  if (cs%use_CVMix_ddiff) &
2227  call cvmix_ddiff_end(cs%CVMix_ddiff_csp)
2228 
2229  if (associated(cs)) deallocate(cs)
2230 

◆ set_diffusivity_init()

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

Definition at line 1885 of file MOM_set_diffusivity.F90.

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