MOM6
mom_entrain_diffusive Module Reference

Detailed Description

Diapycnal mixing and advection in isopycnal mode.

By Robert Hallberg, September 1997 - July 2000

This file contains the subroutines that implement diapycnal mixing and advection in isopycnal layers. The main subroutine, calculate_entrainment, returns the entrainment by each layer across the interfaces above and below it. These are calculated subject to the constraints that no layers can be driven to neg- ative thickness and that the each layer maintains its target density, using the scheme described in Hallberg (MWR 2000). There may or may not be a bulk mixed layer above the isopycnal layers. The solution is iterated until the change in the entrainment between successive iterations is less than some small tolerance.

The dual-stream entrainment scheme of MacDougall and Dewar (JPO 1997) is used for combined diapycnal advection and diffusion, modified as described in Hallberg (MWR 2000) to be solved implicitly in time. Any profile of diffusivities may be used. Diapycnal advection is fundamentally the residual of diapycnal diffusion, so the fully implicit upwind differencing scheme that is used is entirely appropriate. The downward buoyancy flux in each layer is determined from an implicit calculation based on the previously calculated flux of the layer above and an estim- ated flux in the layer below. This flux is subject to the foll- owing conditions: (1) the flux in the top and bottom layers are set by the boundary conditions, and (2) no layer may be driven below an Angstrom thickness. If there is a bulk mixed layer, the mixed and buffer layers are treated as Eulerian layers, whose thicknesses only change due to entrainment by the interior layers.

Data Types

type  entrain_diffusive_cs
 The control structure holding parametes for the MOM_entrain_diffusive module. More...
 

Functions/Subroutines

subroutine, public entrainment_diffusive (h, tv, fluxes, dt, G, GV, US, CS, ea, eb, kb_out, Kd_Lay, Kd_int)
 This subroutine calculates ea and eb, the rates at which a layer entrains from the layers above and below. The entrainment rates are proportional to the buoyancy flux in a layer and inversely proportional to the density differences between layers. The scheme that is used here is described in detail in Hallberg, Mon. Wea. Rev. 2000. More...
 
subroutine f_to_ent (F, h, kb, kmb, j, G, GV, CS, dsp1_ds, eakb, Ent_bl, ea, eb, do_i_in)
 This subroutine calculates the actual entrainments (ea and eb) and the amount of surface forcing that is applied to each layer if there is no bulk mixed layer. More...
 
subroutine set_ent_bl (h, dtKd_int, tv, kb, kmb, do_i, G, GV, US, CS, j, Ent_bl, Sref, h_bl)
 This subroutine sets the average entrainment across each of the interfaces between buffer layers within a timestep. It also causes thin and relatively light interior layers to be entrained by the deepest buffer layer. Also find the initial coordinate potential densities (Sref) of each layer. More...
 
subroutine determine_dskb (h_bl, Sref, Ent_bl, E_kb, is, ie, kmb, G, GV, limit, dSkb, ddSkb_dE, dSlay, ddSlay_dE, dS_anom_lim, do_i_in)
 This subroutine determines the reference density difference between the bottommost buffer layer and the first interior after the mixing between mixed and buffer layers and mixing with the layer below. Within the mixed and buffer layers, entrainment from the layer above is increased when it is necessary to keep the layers from developing a negative thickness; otherwise it equals Ent_bl. At each interface, the upward and downward fluxes average out to Ent_bl, unless entrainment by the layer below is larger than twice Ent_bl. The density difference across the first interior layer may also be returned. It could also be limited to avoid negative values or values that greatly exceed the density differences across an interface. Additionally, the partial derivatives of dSkb and dSlay with E_kb could also be returned. More...
 
subroutine f_kb_to_ea_kb (h_bl, Sref, Ent_bl, I_dSkbp1, F_kb, kmb, i, G, GV, CS, ea_kb, tol_in)
 Given an entrainment from below for layer kb, determine a consistent entrainment from above, such that dSkb * ea_kb = dSkbp1 * F_kb. The input value of ea_kb is both the maximum value that can be obtained and the first guess of the iterations. Ideally ea_kb should be an under-estimate. More...
 
subroutine determine_ea_kb (h_bl, dtKd_kb, Sref, I_dSkbp1, Ent_bl, ea_kbp1, min_eakb, max_eakb, kmb, is, ie, do_i, G, GV, CS, Ent, error, err_min_eakb0, err_max_eakb0, F_kb, dFdfm_kb)
 This subroutine determines the entrainment from above by the top interior layer (labeled kb elsewhere) given an entrainment by the layer below it, constrained to be within the provided bounds. More...
 
subroutine find_maxf_kb (h_bl, Sref, Ent_bl, I_dSkbp1, min_ent_in, max_ent_in, kmb, is, ie, G, GV, CS, maxF, ent_maxF, do_i_in, F_lim_maxent, F_thresh)
 Maximize F = ent*ds_kb*I_dSkbp1 in the range min_ent < ent < max_ent. More...
 
subroutine, public entrain_diffusive_init (Time, G, GV, US, param_file, diag, CS)
 This subroutine initializes the parameters and memory associated with the entrain_diffusive module. More...
 
subroutine, public entrain_diffusive_end (CS)
 This subroutine cleans up and deallocates any memory associated with the entrain_diffusive module. More...
 

Function/Subroutine Documentation

◆ determine_dskb()

subroutine mom_entrain_diffusive::determine_dskb ( real, dimension(szi_(g),szk_(g)), intent(in)  h_bl,
real, dimension(szi_(g),szk_(g)), intent(in)  Sref,
real, dimension(szi_(g),szk_(g)), intent(in)  Ent_bl,
real, dimension(szi_(g)), intent(in)  E_kb,
integer, intent(in)  is,
integer, intent(in)  ie,
integer, intent(in)  kmb,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
logical, intent(in)  limit,
real, dimension(szi_(g)), intent(inout)  dSkb,
real, dimension(szi_(g)), intent(inout), optional  ddSkb_dE,
real, dimension(szi_(g)), intent(inout), optional  dSlay,
real, dimension(szi_(g)), intent(inout), optional  ddSlay_dE,
real, dimension(szi_(g)), intent(inout), optional  dS_anom_lim,
logical, dimension(szi_(g)), intent(in), optional  do_i_in 
)
private

This subroutine determines the reference density difference between the bottommost buffer layer and the first interior after the mixing between mixed and buffer layers and mixing with the layer below. Within the mixed and buffer layers, entrainment from the layer above is increased when it is necessary to keep the layers from developing a negative thickness; otherwise it equals Ent_bl. At each interface, the upward and downward fluxes average out to Ent_bl, unless entrainment by the layer below is larger than twice Ent_bl. The density difference across the first interior layer may also be returned. It could also be limited to avoid negative values or values that greatly exceed the density differences across an interface. Additionally, the partial derivatives of dSkb and dSlay with E_kb could also be returned.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]h_blLayer thickness [H ~> m or kg m-2]
[in]srefReference potential density [R ~> kg m-3]
[in]ent_blThe average entrainment upward and downward across each interface around the buffer layers [H ~> m or kg m-2].
[in]e_kbThe entrainment by the top interior layer [H ~> m or kg m-2].
[in]isThe start of the i-index range to work on.
[in]ieThe end of the i-index range to work on.
[in]kmbThe number of mixed and buffer layers.
[in]limitIf true, limit dSkb and dSlay to avoid negative values.
[in,out]dskbThe limited potential density difference across the interface between the bottommost buffer layer and the topmost interior layer. [R ~> kg m-3] dSkb > 0.
[in,out]ddskb_deThe partial derivative of dSkb with E [R H-1 ~> kg m-4 or m-1].
[in,out]dslayThe limited potential density difference across the topmost interior layer. 0 < dSkb [R ~> kg m-3]
[in,out]ddslay_deThe partial derivative of dSlay with E [R H-1 ~> kg m-4 or m-1].
[in,out]ds_anom_limA limiting value to use for the density anomalies below the buffer layer [R ~> kg m-3].
[in]do_i_inIf present, determines which columns are worked on.

Definition at line 1204 of file MOM_entrain_diffusive.F90.

1204  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1205  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid
1206  !! structure.
1207  real, dimension(SZI_(G),SZK_(G)), intent(in) :: h_bl !< Layer thickness [H ~> m or kg m-2]
1208  real, dimension(SZI_(G),SZK_(G)), intent(in) :: Sref !< Reference potential density [R ~> kg m-3]
1209  real, dimension(SZI_(G),SZK_(G)), intent(in) :: Ent_bl !< The average entrainment upward and
1210  !! downward across each interface
1211  !! around the buffer layers [H ~> m or kg m-2].
1212  real, dimension(SZI_(G)), intent(in) :: E_kb !< The entrainment by the top interior
1213  !! layer [H ~> m or kg m-2].
1214  integer, intent(in) :: is !< The start of the i-index range to work on.
1215  integer, intent(in) :: ie !< The end of the i-index range to work on.
1216  integer, intent(in) :: kmb !< The number of mixed and buffer layers.
1217  logical, intent(in) :: limit !< If true, limit dSkb and dSlay to
1218  !! avoid negative values.
1219  real, dimension(SZI_(G)), intent(inout) :: dSkb !< The limited potential density
1220  !! difference across the interface
1221  !! between the bottommost buffer layer
1222  !! and the topmost interior layer. [R ~> kg m-3]
1223  !! dSkb > 0.
1224  real, dimension(SZI_(G)), optional, intent(inout) :: ddSkb_dE !< The partial derivative of dSkb
1225  !! with E [R H-1 ~> kg m-4 or m-1].
1226  real, dimension(SZI_(G)), optional, intent(inout) :: dSlay !< The limited potential density
1227  !! difference across the topmost
1228  !! interior layer. 0 < dSkb [R ~> kg m-3]
1229  real, dimension(SZI_(G)), optional, intent(inout) :: ddSlay_dE !< The partial derivative of dSlay
1230  !! with E [R H-1 ~> kg m-4 or m-1].
1231  real, dimension(SZI_(G)), optional, intent(inout) :: dS_anom_lim !< A limiting value to use for
1232  !! the density anomalies below the
1233  !! buffer layer [R ~> kg m-3].
1234  logical, dimension(SZI_(G)), optional, intent(in) :: do_i_in !< If present, determines which
1235  !! columns are worked on.
1236 
1237 ! Note that dSkb, ddSkb_dE, dSlay, ddSlay_dE, and dS_anom_lim are declared
1238 ! intent inout because they should not change where do_i_in is false.
1239 
1240 ! This subroutine determines the reference density difference between the
1241 ! bottommost buffer layer and the first interior after the mixing between mixed
1242 ! and buffer layers and mixing with the layer below. Within the mixed and buffer
1243 ! layers, entrainment from the layer above is increased when it is necessary to
1244 ! keep the layers from developing a negative thickness; otherwise it equals
1245 ! Ent_bl. At each interface, the upward and downward fluxes average out to
1246 ! Ent_bl, unless entrainment by the layer below is larger than twice Ent_bl.
1247 ! The density difference across the first interior layer may also be returned.
1248 ! It could also be limited to avoid negative values or values that greatly
1249 ! exceed the density differences across an interface.
1250 ! Additionally, the partial derivatives of dSkb and dSlay with E_kb could
1251 ! also be returned.
1252 
1253  ! Local variables
1254  real, dimension(SZI_(G),SZK_(G)) :: &
1255  b1, c1, & ! b1 and c1 are variables used by the tridiagonal solver.
1256  S, dS_dE, & ! The coordinate density [R ~> kg m-3] and its derivative with E.
1257  ea, dea_dE, & ! The entrainment from above and its derivative with E.
1258  eb, deb_dE ! The entrainment from below and its derivative with E.
1259  real :: deriv_dSkb(SZI_(G))
1260  real :: d1(SZI_(G)) ! d1 = 1.0-c1 is also used by the tridiagonal solver.
1261  real :: src ! A source term for dS_dR.
1262  real :: h1 ! The thickness in excess of the minimum that will remain
1263  ! after exchange with the layer below [H ~> m or kg m-2].
1264  logical, dimension(SZI_(G)) :: do_i
1265  real :: h_neglect ! A thickness that is so small it is usually lost
1266  ! in roundoff and can be neglected [H ~> m or kg m-2].
1267  real :: h_tr ! h_tr is h at tracer points with a tiny thickness
1268  ! added to ensure positive definiteness [H ~> m or kg m-2].
1269  real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2].
1270  real :: rat
1271  real :: dS_kbp1, IdS_kbp1
1272  real :: deriv_dSLay
1273  real :: Inv_term ! [nondim]
1274  real :: f1, df1_drat ! Temporary variables [nondim].
1275  real :: z, dz_drat, f2, df2_dz, expz ! Temporary variables [nondim].
1276  real :: eps_dSLay, eps_dSkb ! Small nondimensional constants.
1277  integer :: i, k
1278 
1279  if (present(ddslay_de) .and. .not.present(dslay)) call mom_error(fatal, &
1280  "In deterimine_dSkb, ddSLay_dE may only be present if dSlay is.")
1281 
1282  h_neglect = gv%H_subroundoff
1283 
1284  do i=is,ie
1285  ea(i,kmb+1) = e_kb(i) ; dea_de(i,kmb+1) = 1.0
1286  s(i,kmb+1) = sref(i,kmb+1) ; ds_de(i,kmb+1) = 0.0
1287  b1(i,kmb+1) = 0.0
1288  d1(i) = 1.0
1289  do_i(i) = .true.
1290  enddo
1291  if (present(do_i_in)) then
1292  do i=is,ie ; do_i(i) = do_i_in(i) ; enddo
1293  endif
1294  do k=kmb,1,-1 ; do i=is,ie
1295  if (do_i(i)) then
1296  ! The do_i test here is only for efficiency.
1297  ! Determine the entrainment from below for each buffer layer.
1298  if (2.0*ent_bl(i,k+1) > ea(i,k+1)) then
1299  eb(i,k) = 2.0*ent_bl(i,k+1) - ea(i,k+1) ; deb_de(i,k) = -dea_de(i,k+1)
1300  else
1301  eb(i,k) = 0.0 ; deb_de(i,k) = 0.0
1302  endif
1303 
1304  ! Determine the entrainment from above for each buffer layer.
1305  h1 = (h_bl(i,k) - gv%Angstrom_H) + (eb(i,k) - ea(i,k+1))
1306  if (h1 >= 0.0) then
1307  ea(i,k) = ent_bl(i,k) ; dea_de(i,k) = 0.0
1308  elseif (ent_bl(i,k) + 0.5*h1 >= 0.0) then
1309  ea(i,k) = ent_bl(i,k) - 0.5*h1
1310  dea_de(i,k) = 0.5*(dea_de(i,k+1) - deb_de(i,k))
1311  else
1312  ea(i,k) = -h1
1313  dea_de(i,k) = dea_de(i,k+1) - deb_de(i,k)
1314  endif
1315  else
1316  ea(i,k) = 0.0 ; dea_de(i,k) = 0.0 ; eb(i,k) = 0.0 ; deb_de(i,k) = 0.0
1317  endif
1318 
1319  ! This is the first-pass of a tridiagonal solver for S.
1320  h_tr = h_bl(i,k) + h_neglect
1321  c1(i,k) = ea(i,k+1) * b1(i,k+1)
1322  b_denom_1 = (h_tr + d1(i)*eb(i,k))
1323  b1(i,k) = 1.0 / (b_denom_1 + ea(i,k))
1324  d1(i) = b_denom_1 * b1(i,k)
1325 
1326  s(i,k) = (h_tr*sref(i,k) + eb(i,k)*s(i,k+1)) * b1(i,k)
1327  enddo ; enddo
1328  do k=2,kmb ; do i=is,ie
1329  s(i,k) = s(i,k) + c1(i,k-1)*s(i,k-1)
1330  enddo ; enddo
1331 
1332  if (present(ddskb_de) .or. present(ddslay_de)) then
1333  ! These two tridiagonal solvers cannot be combined because the solutions for
1334  ! S are required as a source for dS_dE.
1335  do k=kmb,2,-1 ; do i=is,ie
1336  if (do_i(i) .and. (dea_de(i,k) - deb_de(i,k) > 0.0)) then
1337  src = (((s(i,k+1) - sref(i,k)) * (h_bl(i,k) + h_neglect) + &
1338  (s(i,k+1) - s(i,k-1)) * ea(i,k)) * deb_de(i,k) - &
1339  ((sref(i,k) - s(i,k-1)) * h_bl(i,k) + &
1340  (s(i,k+1) - s(i,k-1)) * eb(i,k)) * dea_de(i,k)) / &
1341  ((h_bl(i,k) + h_neglect + ea(i,k)) + eb(i,k))
1342  else ; src = 0.0 ; endif
1343  ds_de(i,k) = (src + eb(i,k)*ds_de(i,k+1)) * b1(i,k)
1344  enddo ; enddo
1345  do i=is,ie
1346  if (do_i(i) .and. (deb_de(i,1) < 0.0)) then
1347  src = (((s(i,2) - sref(i,1)) * (h_bl(i,1) + h_neglect)) * deb_de(i,1)) / &
1348  (h_bl(i,1) + h_neglect + eb(i,1))
1349  else ; src = 0.0 ; endif
1350  ds_de(i,1) = (src + eb(i,1)*ds_de(i,2)) * b1(i,1)
1351  enddo
1352  do k=2,kmb ; do i=is,ie
1353  ds_de(i,k) = ds_de(i,k) + c1(i,k-1)*ds_de(i,k-1)
1354  enddo ; enddo
1355  endif
1356 
1357  ! Now, apply any limiting and return the requested variables.
1358 
1359  eps_dskb = 1.0e-6 ! Should be a small, nondimensional, positive number.
1360  if (.not.limit) then
1361  do i=is,ie ; if (do_i(i)) then
1362  dskb(i) = sref(i,kmb+1) - s(i,kmb)
1363  endif ; enddo
1364  if (present(ddskb_de)) then ; do i=is,ie ; if (do_i(i)) then
1365  ddskb_de(i) = -1.0*ds_de(i,kmb)
1366  endif ; enddo ; endif
1367 
1368  if (present(dslay)) then ; do i=is,ie ; if (do_i(i)) then
1369  dslay(i) = 0.5 * (sref(i,kmb+2) - s(i,kmb))
1370  endif ; enddo ; endif
1371  if (present(ddslay_de)) then ; do i=is,ie ; if (do_i(i)) then
1372  ddslay_de(i) = -0.5*ds_de(i,kmb)
1373  endif ; enddo ; endif
1374  else
1375  do i=is,ie ; if (do_i(i)) then
1376  ! Need to ensure that 0 < dSkb <= S_kb - Sbl
1377  if (sref(i,kmb+1) - s(i,kmb) < eps_dskb*(sref(i,kmb+2) - sref(i,kmb+1))) then
1378  dskb(i) = eps_dskb * (sref(i,kmb+2) - sref(i,kmb+1)) ; deriv_dskb(i) = 0.0
1379  else
1380  dskb(i) = sref(i,kmb+1) - s(i,kmb) ; deriv_dskb(i) = -1.0
1381  endif
1382  if (present(ddskb_de)) ddskb_de(i) = deriv_dskb(i)*ds_de(i,kmb)
1383  endif ; enddo
1384 
1385  if (present(dslay)) then
1386  dz_drat = 1000.0 ! The limit of large dz_drat the same as choosing a
1387  ! Heaviside function.
1388  eps_dslay = 1.0e-10 ! Should be ~= GV%Angstrom_H / sqrt(Kd*dt)
1389  do i=is,ie ; if (do_i(i)) then
1390  ds_kbp1 = sref(i,kmb+2) - sref(i,kmb+1)
1391  ids_kbp1 = 1.0 / (sref(i,kmb+2) - sref(i,kmb+1))
1392  rat = (sref(i,kmb+1) - s(i,kmb)) * ids_kbp1
1393  ! Need to ensure that 0 < dSLay <= 2*dSkb
1394  if (rat < 0.5) then
1395  ! The coefficients here are chosen so that at rat = 0.5, the value (1.5)
1396  ! and first derivative (-0.5) match with the "typical" case (next).
1397  ! The functional form here is arbitrary.
1398  ! f1 provides a reasonable profile that matches the value and derivative
1399  ! of the "typical" case at rat = 0.5, and has a maximum of less than 2.
1400  inv_term = 1.0 / (1.0-rat)
1401  f1 = 2.0 - 0.125*(inv_term**2)
1402  df1_drat = - 0.25*(inv_term**3)
1403 
1404  ! f2 ensures that dSLay goes to 0 rapidly if rat is significantly
1405  ! negative.
1406  z = dz_drat * rat + 4.0 ! The 4 here gives f2(0) = 0.982.
1407  if (z >= 18.0) then ; f2 = 1.0 ; df2_dz = 0.0
1408  elseif (z <= -58.0) then ; f2 = eps_dslay ; df2_dz = 0.0
1409  else
1410  expz = exp(z) ; inv_term = 1.0 / (1.0 + expz)
1411  f2 = (eps_dslay + expz) * inv_term
1412  df2_dz = (1.0 - eps_dslay) * expz * inv_term**2
1413  endif
1414 
1415  dslay(i) = dskb(i) * f1 * f2
1416  deriv_dslay = deriv_dskb(i) * (f1 * f2) - (dskb(i)*ids_kbp1) * &
1417  (df1_drat*f2 + f1 * dz_drat * df2_dz)
1418  elseif (dskb(i) <= 3.0*ds_kbp1) then
1419  ! This is the "typical" case.
1420  dslay(i) = 0.5 * (dskb(i) + ds_kbp1)
1421  deriv_dslay = 0.5 * deriv_dskb(i) ! = -0.5
1422  else
1423  dslay(i) = 2.0*ds_kbp1
1424  deriv_dslay = 0.0
1425  endif
1426  if (present(ddslay_de)) ddslay_de(i) = deriv_dslay*ds_de(i,kmb)
1427  endif ; enddo
1428  endif ! present(dSlay)
1429  endif ! Not limited.
1430 
1431  if (present(ds_anom_lim)) then ; do i=is,ie ; if (do_i(i)) then
1432  ds_anom_lim(i) = max(0.0, eps_dskb * (sref(i,kmb+2) - sref(i,kmb+1)) - &
1433  (sref(i,kmb+1) - s(i,kmb)) )
1434  endif ; enddo ; endif
1435 

References mom_error_handler::mom_error().

Referenced by determine_ea_kb(), entrainment_diffusive(), f_kb_to_ea_kb(), and find_maxf_kb().

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

◆ determine_ea_kb()

subroutine mom_entrain_diffusive::determine_ea_kb ( real, dimension(szi_(g),szk_(g)), intent(in)  h_bl,
real, dimension(szi_(g)), intent(in)  dtKd_kb,
real, dimension(szi_(g),szk_(g)), intent(in)  Sref,
real, dimension(szi_(g)), intent(in)  I_dSkbp1,
real, dimension(szi_(g),szk_(g)), intent(in)  Ent_bl,
real, dimension(szi_(g)), intent(in)  ea_kbp1,
real, dimension(szi_(g)), intent(in)  min_eakb,
real, dimension(szi_(g)), intent(in)  max_eakb,
integer, intent(in)  kmb,
integer, intent(in)  is,
integer, intent(in)  ie,
logical, dimension(szi_(g)), intent(in)  do_i,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(entrain_diffusive_cs), pointer  CS,
real, dimension( g %isd: g %ied), intent(inout)  Ent,
real, dimension( g %isd: g %ied), intent(out), optional  error,
real, dimension( g %isd: g %ied), intent(in), optional  err_min_eakb0,
real, dimension( g %isd: g %ied), intent(in), optional  err_max_eakb0,
real, dimension( g %isd: g %ied), intent(out), optional  F_kb,
real, dimension( g %isd: g %ied), intent(out), optional  dFdfm_kb 
)
private

This subroutine determines the entrainment from above by the top interior layer (labeled kb elsewhere) given an entrainment by the layer below it, constrained to be within the provided bounds.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]h_blLayer thickness, with the top interior layer at k-index kmb+1 [H ~> m or kg m-2].
[in]srefThe coordinate reference potential density, with the value of the topmost interior layer at layer kmb+1 [R ~> kg m-3].
[in]ent_blThe average entrainment upward and downward across each interface around the buffer layers [H ~> m or kg m-2].
[in]i_dskbp1The inverse of the difference in reference potential density across the base of the uppermost interior layer [R-1 ~> m3 kg-1].
[in]dtkd_kbThe diapycnal diffusivity in the top interior layer times the time step [H2 ~> m2 or kg2 m-4].
[in]ea_kbp1The entrainment from above by layer kb+1 [H ~> m or kg m-2].
[in]min_eakbThe minimum permissible rate of entrainment [H ~> m or kg m-2].
[in]max_eakbThe maximum permissible rate of entrainment [H ~> m or kg m-2].
[in]kmbThe number of mixed and buffer layers.
[in]isThe start of the i-index range to work on.
[in]ieThe end of the i-index range to work on.
[in]do_iA logical variable indicating which i-points to work on.
csThis module's control structure.
[in,out]entThe entrainment rate of the uppermost interior layer [H ~> m or kg m-2]. The input value is the first guess.
[out]errorThe error (locally defined in this routine) associated with the returned solution.
[in]err_min_eakb0The errors (locally defined) associated with min_eakb when ea_kbp1 = 0, returned from a previous call to this fn.
[in]err_max_eakb0The errors (locally defined) associated with min_eakb when ea_kbp1 = 0, returned from a previous call to this fn.
[out]f_kbThe entrainment from below by the uppermost interior layer corresponding to the returned value of Ent [H ~> m or kg m-2].
[out]dfdfm_kbThe partial derivative of F_kb with ea_kbp1 [nondim].

Definition at line 1577 of file MOM_entrain_diffusive.F90.

1577  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1578  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1579  real, dimension(SZI_(G),SZK_(G)), intent(in) :: h_bl !< Layer thickness, with the top interior
1580  !! layer at k-index kmb+1 [H ~> m or kg m-2].
1581  real, dimension(SZI_(G),SZK_(G)), intent(in) :: Sref !< The coordinate reference potential
1582  !! density, with the value of the
1583  !! topmost interior layer at layer
1584  !! kmb+1 [R ~> kg m-3].
1585  real, dimension(SZI_(G),SZK_(G)), intent(in) :: Ent_bl !< The average entrainment upward and
1586  !! downward across each interface around
1587  !! the buffer layers [H ~> m or kg m-2].
1588  real, dimension(SZI_(G)), intent(in) :: I_dSkbp1 !< The inverse of the difference in
1589  !! reference potential density across
1590  !! the base of the uppermost interior
1591  !! layer [R-1 ~> m3 kg-1].
1592  real, dimension(SZI_(G)), intent(in) :: dtKd_kb !< The diapycnal diffusivity in the top
1593  !! interior layer times the time step
1594  !! [H2 ~> m2 or kg2 m-4].
1595  real, dimension(SZI_(G)), intent(in) :: ea_kbp1 !< The entrainment from above by layer
1596  !! kb+1 [H ~> m or kg m-2].
1597  real, dimension(SZI_(G)), intent(in) :: min_eakb !< The minimum permissible rate of
1598  !! entrainment [H ~> m or kg m-2].
1599  real, dimension(SZI_(G)), intent(in) :: max_eakb !< The maximum permissible rate of
1600  !! entrainment [H ~> m or kg m-2].
1601  integer, intent(in) :: kmb !< The number of mixed and buffer layers.
1602  integer, intent(in) :: is !< The start of the i-index range to work on.
1603  integer, intent(in) :: ie !< The end of the i-index range to work on.
1604  logical, dimension(SZI_(G)), intent(in) :: do_i !< A logical variable indicating which
1605  !! i-points to work on.
1606  type(entrain_diffusive_CS), pointer :: CS !< This module's control structure.
1607  real, dimension(SZI_(G)), intent(inout) :: Ent !< The entrainment rate of the uppermost
1608  !! interior layer [H ~> m or kg m-2].
1609  !! The input value is the first guess.
1610  real, dimension(SZI_(G)), optional, intent(out) :: error !< The error (locally defined in this
1611  !! routine) associated with the returned
1612  !! solution.
1613  real, dimension(SZI_(G)), optional, intent(in) :: err_min_eakb0 !< The errors (locally defined)
1614  !! associated with min_eakb when ea_kbp1 = 0,
1615  !! returned from a previous call to this fn.
1616  real, dimension(SZI_(G)), optional, intent(in) :: err_max_eakb0 !< The errors (locally defined)
1617  !! associated with min_eakb when ea_kbp1 = 0,
1618  !! returned from a previous call to this fn.
1619  real, dimension(SZI_(G)), optional, intent(out) :: F_kb !< The entrainment from below by the
1620  !! uppermost interior layer
1621  !! corresponding to the returned
1622  !! value of Ent [H ~> m or kg m-2].
1623  real, dimension(SZI_(G)), optional, intent(out) :: dFdfm_kb !< The partial derivative of F_kb with
1624  !! ea_kbp1 [nondim].
1625 
1626 ! This subroutine determines the entrainment from above by the top interior
1627 ! layer (labeled kb elsewhere) given an entrainment by the layer below it,
1628 ! constrained to be within the provided bounds.
1629 
1630  ! Local variables
1631  real, dimension(SZI_(G)) :: &
1632  dS_kb, & ! The coordinate-density difference between the
1633  ! layer kb and deepest buffer layer, limited to
1634  ! ensure that it is positive [R ~> kg m-3].
1635  ds_lay, & ! The coordinate-density difference across layer
1636  ! kb, limited to ensure that it is positive and not
1637  ! too much bigger than dS_kb or dS_kbp1 [R ~> kg m-3].
1638  ddskb_de, ddslay_de, & ! The derivatives of dS_kb and dS_Lay with E
1639  ! [R H-1 ~> kg m-4 or m-1].
1640  derror_de, & ! The derivative of err with E [H ~> m or kg m-2].
1641  err, & ! The "error" whose zero is being sought [H2 ~> m2 or kg2 m-4].
1642  e_min, e_max, & ! The minimum and maximum values of E [H ~> m or kg m-2].
1643  error_mine, error_maxe ! err when E = E_min or E = E_max [H2 ~> m2 or kg2 m-4].
1644  real :: err_est ! An estimate of what err will be [H2 ~> m2 or kg2 m-4].
1645  real :: eL ! 1 or 0, depending on whether increases in E lead
1646  ! to decreases in the entrainment from below by the
1647  ! deepest buffer layer [nondim].
1648  real :: fa ! Temporary variable used to calculate err [nondim].
1649  real :: fk ! Temporary variable used to calculate err [H2 ~> m2 or kg2 m-4].
1650  real :: fm, fr ! Temporary variables used to calculate err [H ~> m or kg m-2].
1651  real :: tolerance ! The tolerance within which E must be converged [H ~> m or kg m-2].
1652  real :: E_prev ! The previous value of E [H ~> m or kg m-2].
1653  logical, dimension(SZI_(G)) :: false_position ! If true, the false position
1654  ! method might be used for the next iteration.
1655  logical, dimension(SZI_(G)) :: redo_i ! If true, more work is needed on this column.
1656  logical :: do_any
1657  real :: large_err ! A large error measure [H2 ~> m2 or kg2 m-4].
1658  integer :: i, it
1659  integer, parameter :: MAXIT = 30
1660 
1661  if (.not.cs%bulkmixedlayer) then
1662  call mom_error(fatal, "determine_Ea_kb should not be called "//&
1663  "unless BULKMIXEDLAYER is defined.")
1664  endif
1665  tolerance = cs%Tolerance_Ent
1666  large_err = gv%m_to_H**2 * 1.0e30
1667 
1668  do i=is,ie ; redo_i(i) = do_i(i) ; enddo
1669 
1670  do i=is,ie ; if (do_i(i)) then
1671  ! The first guess of Ent was the value from the previous iteration.
1672 
1673  ! These were previously calculated and provide good limits and estimates
1674  ! of the errors there. By construction the errors increase with R*ea_kbp1.
1675  e_min(i) = min_eakb(i) ; e_max(i) = max_eakb(i)
1676  error_mine(i) = -large_err ; error_maxe(i) = large_err
1677  false_position(i) = .true. ! Used to alternate between false_position and
1678  ! bisection when Newton's method isn't working.
1679  if (present(err_min_eakb0)) error_mine(i) = err_min_eakb0(i) - e_min(i) * ea_kbp1(i)
1680  if (present(err_max_eakb0)) error_maxe(i) = err_max_eakb0(i) - e_max(i) * ea_kbp1(i)
1681 
1682  if ((error_maxe(i) <= 0.0) .or. (error_mine(i) >= 0.0)) then
1683  ! The root is not bracketed and one of the limiting values should be used.
1684  if (error_maxe(i) <= 0.0) then
1685  ! The errors decrease with E*ea_kbp1, so E_max is the best solution.
1686  ent(i) = e_max(i) ; err(i) = error_maxe(i)
1687  else ! error_minE >= 0 is equivalent to ea_kbp1 = 0.0.
1688  ent(i) = e_min(i) ; err(i) = error_mine(i)
1689  endif
1690  derror_de(i) = 0.0
1691  redo_i(i) = .false.
1692  endif
1693  endif ; enddo ! End of i-loop
1694 
1695  do it = 1,maxit
1696  do_any = .false. ; do i=is,ie ; if (redo_i(i)) do_any = .true. ; enddo
1697  if (.not.do_any) exit
1698  call determine_dskb(h_bl, sref, ent_bl, ent, is, ie, kmb, g, gv, .true., ds_kb, &
1699  ddskb_de, ds_lay, ddslay_de, do_i_in = redo_i)
1700  do i=is,ie ; if (redo_i(i)) then
1701  ! The correct root is bracketed between E_min and E_max.
1702  ! Note the following limits: Ent >= 0 ; fa > 1 ; fk > 0
1703  el = 0.0 ; if (2.0*ent_bl(i,kmb+1) >= ent(i)) el = 1.0
1704  fa = (1.0 + el) + ds_kb(i)*i_dskbp1(i)
1705  fk = dtkd_kb(i) * (ds_lay(i)/ds_kb(i))
1706  fm = (ea_kbp1(i) - h_bl(i,kmb+1)) + el*2.0*ent_bl(i,kmb+1)
1707  if (fm > -gv%Angstrom_H) fm = fm + gv%Angstrom_H ! This could be smooth if need be.
1708  err(i) = (fa * ent(i)**2 - fm * ent(i)) - fk
1709  derror_de(i) = ((2.0*fa + (ddskb_de(i)*i_dskbp1(i))*ent(i))*ent(i) - fm) - &
1710  dtkd_kb(i) * (ddslay_de(i)*ds_kb(i) - ddskb_de(i)*ds_lay(i))/(ds_kb(i)**2)
1711 
1712  if (err(i) == 0.0) then
1713  redo_i(i) = .false. ; cycle
1714  elseif (err(i) > 0.0) then
1715  e_max(i) = ent(i) ; error_maxe(i) = err(i)
1716  else
1717  e_min(i) = ent(i) ; error_mine(i) = err(i)
1718  endif
1719 
1720  e_prev = ent(i)
1721  if ((it == 1) .or. (derror_de(i) <= 0.0)) then
1722  ! Assuming that the coefficients of the quadratic equation are correct
1723  ! will usually give a very good first guess. Also, if derror_dE < 0.0,
1724  ! R is on the wrong side of the approximate parabola. In either case,
1725  ! try assuming that the error is approximately a parabola and solve.
1726  fr = sqrt(fm**2 + 4.0*fa*fk)
1727  if (fm >= 0.0) then
1728  ent(i) = (fm + fr) / (2.0 * fa)
1729  else
1730  ent(i) = (2.0 * fk) / (fr - fm)
1731  endif
1732  ! But make sure that the root stays bracketed, bisecting if needed.
1733  if ((ent(i) > e_max(i)) .or. (ent(i) < e_min(i))) &
1734  ent(i) = 0.5*(e_max(i) + e_min(i))
1735  elseif (((e_max(i)-ent(i))*derror_de(i) > -err(i)) .and. &
1736  ((ent(i)-e_min(i))*derror_de(i) > err(i)) ) then
1737  ! Use Newton's method for the next estimate, provided it will
1738  ! remain bracketed between Rmin and Rmax.
1739  ent(i) = ent(i) - err(i) / derror_de(i)
1740  elseif (false_position(i) .and. &
1741  (error_maxe(i) - error_mine(i) < 0.9*large_err)) then
1742  ! Use the false postion method if there are decent error estimates.
1743  ent(i) = e_min(i) + (e_max(i)-e_min(i)) * &
1744  (-error_mine(i)/(error_maxe(i) - error_mine(i)))
1745  false_position(i) = .false.
1746  else ! Bisect as a last resort or if the false position method was used last.
1747  ent(i) = 0.5*(e_max(i) + e_min(i))
1748  false_position(i) = .true.
1749  endif
1750 
1751  if (abs(e_prev - ent(i)) < tolerance) then
1752  err_est = err(i) + (ent(i) - e_prev) * derror_de(i)
1753  if ((it > 1) .or. (err_est*err(i) <= 0.0) .or. &
1754  (abs(err_est) < abs(tolerance*derror_de(i)))) redo_i(i) = .false.
1755  endif
1756 
1757  endif ; enddo ! End of i-loop
1758  enddo ! End of iterations to determine Ent(i).
1759 
1760  ! Update the value of dS_kb for consistency with Ent.
1761  if (present(f_kb) .or. present(dfdfm_kb)) &
1762  call determine_dskb(h_bl, sref, ent_bl, ent, is, ie, kmb, g, gv, .true., &
1763  ds_kb, do_i_in = do_i)
1764 
1765  if (present(f_kb)) then ; do i=is,ie ; if (do_i(i)) then
1766  f_kb(i) = ent(i) * (ds_kb(i) * i_dskbp1(i))
1767  endif ; enddo ; endif
1768  if (present(error)) then ; do i=is,ie ; if (do_i(i)) then
1769  error(i) = err(i)
1770  endif ; enddo ; endif
1771  if (present(dfdfm_kb)) then ; do i=is,ie ; if (do_i(i)) then
1772  ! derror_dE and ddSkb_dE are _not_ recalculated here, since dFdfm_kb is
1773  ! only used in Newton's method, and slightly increasing the accuracy of the
1774  ! estimate is unlikely to speed convergence.
1775  if (derror_de(i) > 0.0) then
1776  dfdfm_kb(i) = ((ds_kb(i) + ent(i) * ddskb_de(i)) * i_dskbp1(i)) * &
1777  (ent(i) / derror_de(i))
1778  else ! Use Adcroft's division by 0 convention.
1779  dfdfm_kb(i) = 0.0
1780  endif
1781  endif ; enddo ; endif
1782 

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

Referenced by entrainment_diffusive().

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

◆ entrain_diffusive_end()

subroutine, public mom_entrain_diffusive::entrain_diffusive_end ( type(entrain_diffusive_cs), pointer  CS)

This subroutine cleans up and deallocates any memory associated with the entrain_diffusive module.

Parameters
csA pointer to the control structure for this module that will be deallocated.

Definition at line 2152 of file MOM_entrain_diffusive.F90.

2152  type(entrain_diffusive_CS), pointer :: CS !< A pointer to the control structure for this
2153  !! module that will be deallocated.
2154  if (associated(cs)) deallocate(cs)
2155 

Referenced by mom_diabatic_driver::diabatic_driver_end().

Here is the caller graph for this function:

◆ entrain_diffusive_init()

subroutine, public mom_entrain_diffusive::entrain_diffusive_init ( type(time_type), intent(in)  Time,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(entrain_diffusive_cs), pointer  CS 
)

This subroutine initializes the parameters and memory associated with the entrain_diffusive module.

Parameters
[in]timeThe current model time.
[in]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 that is used to regulate diagnostic output.
csA pointer that is set to point to the control structure.

Definition at line 2084 of file MOM_entrain_diffusive.F90.

2084  type(time_type), intent(in) :: Time !< The current model time.
2085  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2086  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2087  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2088  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
2089  !! parameters.
2090  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
2091  !! output.
2092  type(entrain_diffusive_CS), pointer :: CS !< A pointer that is set to point to the control
2093  !! structure.
2094 ! for this module
2095 ! Arguments: Time - The current model time.
2096 ! (in) G - The ocean's grid structure.
2097 ! (in) GV - The ocean's vertical grid structure.
2098 ! (in) param_file - A structure indicating the open file to parse for
2099 ! model parameter values.
2100 ! (in) diag - A structure that is used to regulate diagnostic output.
2101 ! (in/out) CS - A pointer that is set to point to the control structure
2102 ! for this module
2103  real :: decay_length, dt, Kd
2104 ! This include declares and sets the variable "version".
2105 #include "version_variable.h"
2106  character(len=40) :: mdl = "MOM_entrain_diffusive" ! This module's name.
2107 
2108  if (associated(cs)) then
2109  call mom_error(warning, "entrain_diffusive_init called with an associated "// &
2110  "control structure.")
2111  return
2112  endif
2113  allocate(cs)
2114 
2115  cs%diag => diag
2116 
2117  cs%bulkmixedlayer = (gv%nkml > 0)
2118 
2119 ! Set default, read and log parameters
2120  call log_version(param_file, mdl, version, "")
2121  call get_param(param_file, mdl, "CORRECT_DENSITY", cs%correct_density, &
2122  "If true, and USE_EOS is true, the layer densities are "//&
2123  "restored toward their target values by the diapycnal "//&
2124  "mixing, as described in Hallberg (MWR, 2000).", &
2125  default=.true.)
2126  call get_param(param_file, mdl, "MAX_ENT_IT", cs%max_ent_it, &
2127  "The maximum number of iterations that may be used to "//&
2128  "calculate the interior diapycnal entrainment.", default=5)
2129 ! In this module, KD is only used to set the default for TOLERANCE_ENT. [m2 s-1]
2130  call get_param(param_file, mdl, "KD", kd, fail_if_missing=.true.)
2131  call get_param(param_file, mdl, "DT", dt, &
2132  "The (baroclinic) dynamics time step.", units = "s", &
2133  fail_if_missing=.true.)
2134 ! CS%Tolerance_Ent = MAX(100.0*GV%Angstrom_H,1.0e-4*sqrt(dt*Kd)) !
2135  call get_param(param_file, mdl, "TOLERANCE_ENT", cs%Tolerance_Ent, &
2136  "The tolerance with which to solve for entrainment values.", &
2137  units="m", default=max(100.0*gv%Angstrom_m,1.0e-4*sqrt(dt*kd)), scale=gv%m_to_H)
2138 
2139  cs%Rho_sig_off = 1000.0*us%kg_m3_to_R
2140 
2141  cs%id_Kd = register_diag_field('ocean_model', 'Kd_effective', diag%axesTL, time, &
2142  'Diapycnal diffusivity as applied', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
2143  cs%id_diff_work = register_diag_field('ocean_model', 'diff_work', diag%axesTi, time, &
2144  'Work actually done by diapycnal diffusion across each interface', 'W m-2', &
2145  conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2146 

References mom_error_handler::mom_error().

Here is the call graph for this function:

◆ entrainment_diffusive()

subroutine, public mom_entrain_diffusive::entrainment_diffusive ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
type(forcing), intent(in)  fluxes,
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(entrain_diffusive_cs), pointer  CS,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(out)  ea,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(out)  eb,
integer, dimension( g %isd: g %ied, g %jsd: g %jed), intent(inout), optional  kb_out,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in), optional  Kd_Lay,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke+1), intent(in), optional  Kd_int 
)

This subroutine calculates ea and eb, the rates at which a layer entrains from the layers above and below. The entrainment rates are proportional to the buoyancy flux in a layer and inversely proportional to the density differences between layers. The scheme that is used here is described in detail in Hallberg, Mon. Wea. Rev. 2000.

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]tvA structure containing pointers to any available thermodynamic fields. Absent fields have NULL ptrs.
[in]fluxesA structure of surface fluxes that may be used.
[in]dtThe time increment [T ~> s].
csThe control structure returned by a previous call to entrain_diffusive_init.
[out]eaThe amount of fluid entrained from the layer
[out]ebThe amount of fluid entrained from the layer
[in,out]kb_outThe index of the lightest layer denser than
[in]kd_layThe diapycnal diffusivity of layers
[in]kd_intThe diapycnal diffusivity of interfaces

Definition at line 54 of file MOM_entrain_diffusive.F90.

54  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
55  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
56  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
57  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
58  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
59  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available
60  !! thermodynamic fields. Absent fields have NULL
61  !! ptrs.
62  type(forcing), intent(in) :: fluxes !< A structure of surface fluxes that may
63  !! be used.
64  real, intent(in) :: dt !< The time increment [T ~> s].
65  type(entrain_diffusive_CS), pointer :: CS !< The control structure returned by a previous
66  !! call to entrain_diffusive_init.
67  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
68  intent(out) :: ea !< The amount of fluid entrained from the layer
69  !! above within this time step [H ~> m or kg m-2].
70  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
71  intent(out) :: eb !< The amount of fluid entrained from the layer
72  !! below within this time step [H ~> m or kg m-2].
73  integer, dimension(SZI_(G),SZJ_(G)), &
74  optional, intent(inout) :: kb_out !< The index of the lightest layer denser than
75  !! the buffer layer.
76  ! At least one of the two following arguments must be present.
77  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
78  optional, intent(in) :: Kd_Lay !< The diapycnal diffusivity of layers
79  !! [Z2 T-1 ~> m2 s-1].
80  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
81  optional, intent(in) :: Kd_int !< The diapycnal diffusivity of interfaces
82  !! [Z2 T-1 ~> m2 s-1].
83 
84 ! This subroutine calculates ea and eb, the rates at which a layer entrains
85 ! from the layers above and below. The entrainment rates are proportional to
86 ! the buoyancy flux in a layer and inversely proportional to the density
87 ! differences between layers. The scheme that is used here is described in
88 ! detail in Hallberg, Mon. Wea. Rev. 2000.
89 
90  real, dimension(SZI_(G),SZK_(G)) :: &
91  dtKd ! The layer diapycnal diffusivity times the time step [H2 ~> m2 or kg2 m-4].
92  real, dimension(SZI_(G),SZK_(G)+1) :: &
93  dtKd_int ! The diapycnal diffusivity at the interfaces times the time step [H2 ~> m2 or kg2 m-4]
94  real, dimension(SZI_(G),SZK_(G)) :: &
95  F, & ! The density flux through a layer within a time step divided by the
96  ! density difference across the interface below the layer [H ~> m or kg m-2].
97  maxf, & ! maxF is the maximum value of F that will not deplete all of the
98  ! layers above or below a layer within a timestep [H ~> m or kg m-2].
99  minf, & ! minF is the minimum flux that should be expected in the absence of
100  ! interactions between layers [H ~> m or kg m-2].
101  fprev, &! The previous estimate of F [H ~> m or kg m-2].
102  dfdfm, &! The partial derivative of F with respect to changes in F of the
103  ! neighboring layers. [nondim]
104  h_guess ! An estimate of the layer thicknesses after entrainment, but
105  ! before the entrainments are adjusted to drive the layer
106  ! densities toward their target values [H ~> m or kg m-2].
107  real, dimension(SZI_(G),SZK_(G)+1) :: &
108  Ent_bl ! The average entrainment upward and downward across
109  ! each interface around the buffer layers [H ~> m or kg m-2].
110  real, allocatable, dimension(:,:,:) :: &
111  Kd_eff, & ! The effective diffusivity that actually applies to each
112  ! layer after the effects of boundary conditions are
113  ! considered [Z2 T-1 ~> m2 s-1].
114  diff_work ! The work actually done by diffusion across each
115  ! interface [R Z3 T-3 ~> W m-2]. Sum vertically for the total work.
116 
117  real :: hm, fm, fr, fk ! Work variables with units of H, H, H, and H2.
118 
119  real :: b1(SZI_(G)) ! b1 and c1 are variables used by the
120  real :: c1(SZI_(G),SZK_(G)) ! tridiagonal solver.
121 
122  real, dimension(SZI_(G)) :: &
123  htot, & ! The total thickness above or below a layer [H ~> m or kg m-2].
124  Rcv, & ! Value of the coordinate variable (potential density)
125  ! based on the simulated T and S and P_Ref [R ~> kg m-3].
126  pres, & ! Reference pressure (P_Ref) [Pa].
127  eakb, & ! The entrainment from above by the layer below the buffer
128  ! layer (i.e. layer kb) [H ~> m or kg m-2].
129  ea_kbp1, & ! The entrainment from above by layer kb+1 [H ~> m or kg m-2].
130  eb_kmb, & ! The entrainment from below by the deepest buffer layer [H ~> m or kg m-2].
131  ds_kb, & ! The reference potential density difference across the
132  ! interface between the buffer layers and layer kb [R ~> kg m-3].
133  ds_anom_lim, &! The amount by which dS_kb is reduced when limits are
134  ! applied [R ~> kg m-3].
135  i_dskbp1, & ! The inverse of the potential density difference across the
136  ! interface below layer kb [R-1 ~> m3 kg-1].
137  dtkd_kb, & ! The diapycnal diffusivity in layer kb times the time step
138  ! [H2 ~> m2 or kg2 m-4].
139  maxf_correct, & ! An amount by which to correct maxF due to excessive
140  ! surface heat loss [H ~> m or kg m-2].
141  zeros, & ! An array of all zeros. (Usually used with [H ~> m or kg m-2].)
142  max_eakb, & ! The maximum value of eakb that might be realized [H ~> m or kg m-2].
143  min_eakb, & ! The minimum value of eakb that might be realized [H ~> m or kg m-2].
144  err_max_eakb0, & ! The value of error returned by determine_Ea_kb
145  err_min_eakb0, & ! when eakb = min_eakb and max_eakb and ea_kbp1 = 0.
146  err_eakb0, & ! A value of error returned by determine_Ea_kb.
147  f_kb, & ! The value of F in layer kb, or equivalently the entrainment
148  ! from below by layer kb [H ~> m or kg m-2].
149  dfdfm_kb, & ! The partial derivative of F with fm [nondim]. See dFdfm.
150  maxf_kb, & ! The maximum value of F_kb that might be realized [H ~> m or kg m-2].
151  eakb_maxf, & ! The value of eakb that gives F_kb=maxF_kb [H ~> m or kg m-2].
152  f_kb_maxent ! The value of F_kb when eakb = max_eakb [H ~> m or kg m-2].
153  real, dimension(SZI_(G),SZK_(G)) :: &
154  Sref, & ! The reference potential density of the mixed and buffer layers,
155  ! and of the two lightest interior layers (kb and kb+1) copied
156  ! into layers kmb+1 and kmb+2 [R ~> kg m-3].
157  h_bl ! The thicknesses of the mixed and buffer layers, and of the two
158  ! lightest interior layers (kb and kb+1) copied into layers kmb+1
159  ! and kmb+2 [H ~> m or kg m-2].
160 
161  real, dimension(SZI_(G),SZK_(G)) :: &
162  ds_dsp1, & ! The coordinate variable (sigma-2) difference across an
163  ! interface divided by the difference across the interface
164  ! below it. [nondim]
165  dsp1_ds, & ! The inverse coordinate variable (sigma-2) difference
166  ! across an interface times the difference across the
167  ! interface above it. [nondim]
168  i2p2dsp1_ds, & ! 1 / (2 + 2 * ds_k+1 / ds_k). [nondim]
169  grats ! 2*(2 + ds_k+1 / ds_k + ds_k / ds_k+1) =
170  ! 4*ds_Lay*(1/ds_k + 1/ds_k+1). [nondim]
171 
172  real :: dRHo ! The change in locally referenced potential density between
173  ! the layers above and below an interface [R ~> kg m-3].
174  real :: g_2dt ! 0.5 * G_Earth / dt, times unit conversion factors
175  ! [m3 H-2 s-2 T-1 ~> m s-3 or m7 kg-2 s-3].
176  real, dimension(SZI_(G)) :: &
177  pressure, & ! The pressure at an interface [Pa].
178  T_eos, S_eos, & ! The potential temperature and salinity at which to
179  ! evaluate dRho_dT and dRho_dS [degC] and [ppt].
180  drho_dt, drho_ds ! The partial derivatives of potential density with temperature and
181  ! salinity, [R degC-1 ~> kg m-3 degC-1] and [R ppt-1 ~> kg m-3 ppt-1].
182 
183  real :: tolerance ! The tolerance within which E must be converged [H ~> m or kg m-2].
184  real :: Angstrom ! The minimum layer thickness [H ~> m or kg m-2].
185  real :: h_neglect ! A thickness that is so small it is usually lost
186  ! in roundoff and can be neglected [H ~> m or kg m-2].
187  real :: F_cor ! A correction to the amount of F that is used to
188  ! entrain from the layer above [H ~> m or kg m-2].
189  real :: Kd_here ! The effective diapycnal diffusivity [H2 s-1 ~> m2 s-1 or kg2 m-4 s-1].
190  real :: h_avail ! The thickness that is available for entrainment [H ~> m or kg m-2].
191  real :: dS_kb_eff ! The value of dS_kb after limiting is taken into account.
192  real :: Rho_cor ! The depth-integrated potential density anomaly that
193  ! needs to be corrected for [H kg m-3 ~> kg m-2 or kg2 m-5].
194  real :: ea_cor ! The corrective adjustment to eakb [H ~> m or kg m-2].
195  real :: h1 ! The layer thickness after entrainment through the
196  ! interface below is taken into account [H ~> m or kg m-2].
197  real :: Idt ! The inverse of the time step [T-1 ~> s-1].
198 
199  logical :: do_any
200  logical :: do_entrain_eakb ! True if buffer layer is entrained
201  logical :: do_i(SZI_(G)), did_i(SZI_(G)), reiterate, correct_density
202  integer :: it, i, j, k, is, ie, js, je, nz, K2, kmb
203  integer :: kb(SZI_(G)) ! The value of kb in row j.
204  integer :: kb_min ! The minimum value of kb in the current j-row.
205  integer :: kb_min_act ! The minimum active value of kb in the current j-row.
206  integer :: is1, ie1 ! The minimum and maximum active values of i in the current j-row.
207  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
208  angstrom = gv%Angstrom_H
209  h_neglect = gv%H_subroundoff
210 
211  if (.not. associated(cs)) call mom_error(fatal, &
212  "MOM_entrain_diffusive: Module must be initialized before it is used.")
213 
214  if (.not.(present(kd_lay) .or. present(kd_int))) call mom_error(fatal, &
215  "MOM_entrain_diffusive: Either Kd_Lay or Kd_int must be present in call.")
216 
217  if ((.not.cs%bulkmixedlayer .and. .not.associated(fluxes%buoy)) .and. &
218  (associated(fluxes%lprec) .or. associated(fluxes%evap) .or. &
219  associated(fluxes%sens) .or. associated(fluxes%sw))) then
220  if (is_root_pe()) call mom_error(note, "Calculate_Entrainment: &
221  &The code to handle evaporation and precipitation without &
222  &a bulk mixed layer has not been implemented.")
223  if (is_root_pe()) call mom_error(fatal, &
224  "Either define BULKMIXEDLAYER in MOM_input or use fluxes%buoy &
225  &and a linear equation of state to drive the model.")
226  endif
227 
228  tolerance = cs%Tolerance_Ent
229  kmb = gv%nk_rho_varies
230  k2 = max(kmb+1,2) ; kb_min = k2
231  if (.not. cs%bulkmixedlayer) then
232  kb(:) = 1
233  ! These lines fill in values that are arbitrary, but needed because
234  ! they are used to normalize the buoyancy flux in layer nz.
235  do i=is,ie ; ds_dsp1(i,nz) = 2.0 ; dsp1_ds(i,nz) = 0.5 ; enddo
236  else
237  kb(:) = 0
238  do i=is,ie ; ds_dsp1(i,nz) = 0.0 ; dsp1_ds(i,nz) = 0.0 ; enddo
239  endif
240 
241  if (cs%id_diff_work > 0) allocate(diff_work(g%isd:g%ied,g%jsd:g%jed,nz+1))
242  if (cs%id_Kd > 0) allocate(kd_eff(g%isd:g%ied,g%jsd:g%jed,nz))
243 
244  correct_density = (cs%correct_density .and. associated(tv%eqn_of_state))
245  if (correct_density) then
246  pres(:) = tv%P_Ref
247  else
248  pres(:) = 0.0
249  endif
250 
251  !$OMP parallel do default(none) shared(is,ie,js,je,nz,Kd_Lay,G,GV,US,dt,CS,h,tv, &
252  !$OMP kmb,Angstrom,fluxes,K2,h_neglect,tolerance, &
253  !$OMP ea,eb,correct_density,Kd_int,Kd_eff, &
254  !$OMP diff_work,g_2dt, kb_out) &
255  !$OMP firstprivate(kb,ds_dsp1,dsp1_ds,pres,kb_min) &
256  !$OMP private(dtKd,dtKd_int,do_i,Ent_bl,dtKd_kb,h_bl, &
257  !$OMP I2p2dsp1_ds,grats,htot,max_eakb,I_dSkbp1, &
258  !$OMP zeros,maxF_kb,maxF,ea_kbp1,eakb,Sref, &
259  !$OMP maxF_correct,do_any,do_entrain_eakb, &
260  !$OMP err_min_eakb0,err_max_eakb0,eakb_maxF, &
261  !$OMP min_eakb,err_eakb0,F,minF,hm,fk,F_kb_maxent,&
262  !$OMP F_kb,is1,ie1,kb_min_act,dFdfm_kb,b1,dFdfm, &
263  !$OMP Fprev,fm,fr,c1,reiterate,eb_kmb,did_i, &
264  !$OMP h_avail,h_guess,dS_kb,Rcv,F_cor,dS_kb_eff, &
265  !$OMP Rho_cor,ea_cor,h1,Idt,Kd_here,pressure, &
266  !$OMP T_eos,S_eos,dRho_dT,dRho_dS,dRho,dS_anom_lim)
267  do j=js,je
268  do i=is,ie ; kb(i) = 1 ; enddo
269 
270  if (present(kd_lay)) then
271  do k=1,nz ; do i=is,ie
272  dtkd(i,k) = gv%Z_to_H**2 * (dt * kd_lay(i,j,k))
273  enddo ; enddo
274  if (present(kd_int)) then
275  do k=1,nz+1 ; do i=is,ie
276  dtkd_int(i,k) = gv%Z_to_H**2 * (dt * kd_int(i,j,k))
277  enddo ; enddo
278  else
279  do k=2,nz ; do i=is,ie
280  dtkd_int(i,k) = gv%Z_to_H**2 * (0.5 * dt * (kd_lay(i,j,k-1) + kd_lay(i,j,k)))
281  enddo ; enddo
282  endif
283  else ! Kd_int must be present, or there already would have been an error.
284  do k=1,nz ; do i=is,ie
285  dtkd(i,k) = gv%Z_to_H**2 * (0.5 * dt * (kd_int(i,j,k)+kd_int(i,j,k+1)))
286  enddo ; enddo
287  dO k=1,nz+1 ; do i=is,ie
288  dtkd_int(i,k) = gv%Z_to_H**2 * (dt * kd_int(i,j,k))
289  enddo ; enddo
290  endif
291 
292  do i=is,ie ; do_i(i) = (g%mask2dT(i,j) > 0.5) ; enddo
293  do i=is,ie ; ds_dsp1(i,nz) = 0.0 ; enddo
294  do i=is,ie ; dsp1_ds(i,nz) = 0.0 ; enddo
295 
296  do k=2,nz-1 ; do i=is,ie
297  ds_dsp1(i,k) = gv%g_prime(k) / gv%g_prime(k+1)
298  enddo ; enddo
299 
300  if (cs%bulkmixedlayer) then
301  ! This subroutine determines the averaged entrainment across each
302  ! interface and causes thin and relatively light interior layers to be
303  ! entrained by the deepest buffer layer. This also determines kb.
304  call set_ent_bl(h, dtkd_int, tv, kb, kmb, do_i, g, gv, us, cs, j, ent_bl, sref, h_bl)
305 
306  do i=is,ie
307  dtkd_kb(i) = 0.0 ; if (kb(i) < nz) dtkd_kb(i) = dtkd(i,kb(i))
308  enddo
309  else
310  do i=is,ie ; ent_bl(i,kmb+1) = 0.0 ; enddo
311  endif
312 
313  do k=2,nz-1 ; do i=is,ie
314  dsp1_ds(i,k) = 1.0 / ds_dsp1(i,k)
315  i2p2dsp1_ds(i,k) = 0.5/(1.0+dsp1_ds(i,k))
316  grats(i,k) = 2.0*(2.0+(dsp1_ds(i,k)+ds_dsp1(i,k)))
317  enddo ; enddo
318 
319 ! Determine the maximum flux, maxF, for each of the isopycnal layers.
320 ! Also determine when the fluxes start entraining
321 ! from various buffer or mixed layers, where appropriate.
322  if (cs%bulkmixedlayer) then
323  kb_min = nz
324  do i=is,ie
325  htot(i) = h(i,j,1) - angstrom
326  enddo
327  do k=2,kmb ; do i=is,ie
328  htot(i) = htot(i) + (h(i,j,k) - angstrom)
329  enddo ; enddo
330  do i=is,ie
331  max_eakb(i) = max(ent_bl(i,kmb+1) + 0.5*htot(i), htot(i))
332  i_dskbp1(i) = 1.0 / (sref(i,kmb+2) - sref(i,kmb+1))
333  zeros(i) = 0.0
334  enddo
335 
336  ! Find the maximum amount of entrainment from below that the top
337  ! interior layer could exhibit (maxF(i,kb)), approximating that
338  ! entrainment by (eakb*max(dS_kb/dSkbp1,0)). eakb is in the range
339  ! from 0 to max_eakb.
340  call find_maxf_kb(h_bl, sref, ent_bl, i_dskbp1, zeros, max_eakb, kmb, &
341  is, ie, g, gv, cs, maxf_kb, eakb_maxf, do_i, f_kb_maxent)
342  do i=is,ie ; if (kb(i) <= nz) then
343  maxf(i,kb(i)) = max(0.0, maxf_kb(i), f_kb_maxent(i))
344  if ((maxf_kb(i) > f_kb_maxent(i)) .and. (eakb_maxf(i) >= htot(i))) &
345  max_eakb(i) = eakb_maxf(i)
346  endif ; enddo
347 
348  do i=is,ie ; ea_kbp1(i) = 0.0 ; eakb(i) = max_eakb(i) ; enddo
349  call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, &
350  max_eakb, max_eakb, kmb, is, ie, do_i, g, gv, cs, eakb, &
351  error=err_max_eakb0, f_kb=f_kb)
352 
353  ! The maximum value of F(kb) between htot and max_eakb determines
354  ! what maxF(kb+1) should be.
355  do i=is,ie ; min_eakb(i) = min(htot(i), max_eakb(i)) ; enddo
356  call find_maxf_kb(h_bl, sref, ent_bl, i_dskbp1, min_eakb, max_eakb, &
357  kmb, is, ie, g, gv, cs, f_kb_maxent, do_i_in = do_i)
358 
359  do i=is,ie
360  do_entrain_eakb = .false.
361  ! If error_max_eakb0 < 0, then buffer layers are always all entrained
362  if (do_i(i)) then ; if (err_max_eakb0(i) < 0.0) then
363  do_entrain_eakb = .true.
364  endif ; endif
365 
366  if (do_entrain_eakb) then
367  eakb(i) = max_eakb(i) ; min_eakb(i) = max_eakb(i)
368  else
369  eakb(i) = 0.0 ; min_eakb(i) = 0.0
370  endif
371  enddo
372 
373  ! Find the amount of entrainment of the buffer layers that would occur
374  ! if there were no entrainment by the deeper interior layers. Also find
375  ! how much entrainment of the deeper layers would occur.
376  call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, &
377  zeros, max_eakb, kmb, is, ie, do_i, g, gv, cs, min_eakb, &
378  error=err_min_eakb0, f_kb=f_kb, err_max_eakb0=err_max_eakb0)
379  ! Error_min_eakb0 should be ~0 unless error_max_eakb0 < 0.
380  do i=is,ie ; if ((kb(i)<nz) .and. (kb_min>kb(i))) kb_min = kb(i) ; enddo
381  else
382  ! Without a bulk mixed layer, surface fluxes are applied in this
383  ! subroutine. (Otherwise, they are handled in mixedlayer.)
384  ! Initially the maximum flux in layer zero is given by the surface
385  ! buoyancy flux. It will later be limited if the surface flux is
386  ! too large. Here buoy is the surface buoyancy flux.
387  do i=is,ie
388  maxf(i,1) = 0.0
389  htot(i) = h(i,j,1) - angstrom
390  enddo
391  if (associated(fluxes%buoy)) then ; do i=is,ie
392  maxf(i,1) = gv%Z_to_H * (dt*fluxes%buoy(i,j)) / gv%g_prime(2)
393  enddo ; endif
394  endif
395 
396 ! The following code calculates the maximum flux, maxF, for the interior
397 ! layers.
398  do k=kb_min,nz-1 ; do i=is,ie
399  if ((k == kb(i)+1) .and. cs%bulkmixedlayer) then
400  maxf(i,k) = ds_dsp1(i,k)*(f_kb_maxent(i) + htot(i))
401  htot(i) = htot(i) + (h(i,j,k) - angstrom)
402  elseif (k > kb(i)) then
403  maxf(i,k) = ds_dsp1(i,k)*(maxf(i,k-1) + htot(i))
404  htot(i) = htot(i) + (h(i,j,k) - angstrom)
405  endif
406  enddo ; enddo
407  do i=is,ie
408  maxf(i,nz) = 0.0
409  if (.not.cs%bulkmixedlayer) then
410  maxf_correct(i) = max(0.0, -(maxf(i,nz-1) + htot(i)))
411  endif
412  htot(i) = h(i,j,nz) - angstrom
413  enddo
414  if (.not.cs%bulkmixedlayer) then
415  do_any = .false. ; do i=is,ie ; if (maxf_correct(i) > 0.0) do_any = .true. ; enddo
416  if (do_any) then
417  do k=nz-1,1,-1 ; do i=is,ie
418  maxf(i,k) = maxf(i,k) + maxf_correct(i)
419  maxf_correct(i) = maxf_correct(i) * dsp1_ds(i,k)
420  enddo ; enddo
421  endif
422  endif
423  do k=nz-1,kb_min,-1 ; do i=is,ie ; if (do_i(i)) then
424  if (k >= kb(i)) then
425  maxf(i,k) = min(maxf(i,k),dsp1_ds(i,k+1)*maxf(i,k+1) + htot(i))
426  htot(i) = htot(i) + (h(i,j,k) - angstrom)
427  endif
428  if (k == kb(i)) then
429  if ((maxf(i,k) < f_kb(i)) .or. (maxf(i,k) < maxf_kb(i)) &
430  .and. (eakb_maxf(i) <= max_eakb(i))) then
431  ! In this case, too much was being entrained by the topmost interior
432  ! layer, even with the minimum initial estimate. The buffer layer
433  ! will always entrain the maximum amount.
434  f_kb(i) = maxf(i,k)
435  if ((f_kb(i) <= maxf_kb(i)) .and. (eakb_maxf(i) <= max_eakb(i))) then
436  eakb(i) = eakb_maxf(i)
437  else
438  eakb(i) = max_eakb(i)
439  endif
440  call f_kb_to_ea_kb(h_bl, sref, ent_bl, i_dskbp1, f_kb, kmb, i, &
441  g, gv, cs, eakb, angstrom)
442  if ((eakb(i) < max_eakb(i)) .or. (eakb(i) < min_eakb(i))) then
443  call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, zeros, &
444  eakb, eakb, kmb, i, i, do_i, g, gv, cs, eakb, &
445  error=err_eakb0)
446  if (eakb(i) < max_eakb(i)) then
447  max_eakb(i) = eakb(i) ; err_max_eakb0(i) = err_eakb0(i)
448  endif
449  if (eakb(i) < min_eakb(i)) then
450  min_eakb(i) = eakb(i) ; err_min_eakb0(i) = err_eakb0(i)
451  endif
452  endif
453  endif
454  endif
455  endif ; enddo ; enddo
456  if (.not.cs%bulkmixedlayer) then
457  do i=is,ie
458  maxf(i,1) = min(maxf(i,1),dsp1_ds(i,2)*maxf(i,2) + htot(i))
459  enddo
460  endif
461 
462 ! The following code provides an initial estimate of the flux in
463 ! each layer, F. The initial guess for the layer diffusive flux is
464 ! the smaller of a forward discretization or the maximum diffusive
465 ! flux starting from zero thickness in one time step without
466 ! considering adjacent layers.
467  do i=is,ie
468  f(i,1) = maxf(i,1)
469  f(i,nz) = maxf(i,nz) ; minf(i,nz) = 0.0
470  enddo
471  do k=nz-1,k2,-1
472  do i=is,ie
473  if ((k==kb(i)) .and. (do_i(i))) then
474  eakb(i) = min_eakb(i)
475  minf(i,k) = 0.0
476  elseif ((k>kb(i)) .and. (do_i(i))) then
477 ! Here the layer flux is estimated, assuming no entrainment from
478 ! the surrounding layers. The estimate is a forward (steady) flux,
479 ! limited by the maximum flux for a layer starting with zero
480 ! thickness. This is often a good guess and leads to few iterations.
481  hm = h(i,j,k) + h_neglect
482  ! Note: Tried sqrt((0.5*ds_dsp1(i,k))*dtKd(i,k)) for the second limit,
483  ! but it usually doesn't work as well.
484  f(i,k) = min(maxf(i,k), sqrt(ds_dsp1(i,k)*dtkd(i,k)), &
485  0.5*(ds_dsp1(i,k)+1.0) * (dtkd(i,k) / hm))
486 
487 ! Calculate the minimum flux that can be expected if there is no entrainment
488 ! from the neighboring layers. The 0.9 is used to give used to give a 10%
489 ! known error tolerance.
490  fk = dtkd(i,k) * grats(i,k)
491  minf(i,k) = min(maxf(i,k), &
492  0.9*(i2p2dsp1_ds(i,k) * fk / (hm + sqrt(hm*hm + fk))))
493  if (k==kb(i)) minf(i,k) = 0.0 ! BACKWARD COMPATIBILITY - DELETE LATER?
494  else
495  f(i,k) = 0.0
496  minf(i,k) = 0.0
497  endif
498  enddo ! end of i loop
499  enddo ! end of k loop
500 
501  ! This is where the fluxes are actually calculated.
502 
503  is1 = ie+1 ; ie1 = is-1
504  do i=is,ie ; if (do_i(i)) then ; is1 = i ; exit ; endif ; enddo
505  do i=ie,is,-1 ; if (do_i(i)) then ; ie1 = i ; exit ; endif ; enddo
506 
507  if (cs%bulkmixedlayer) then
508  kb_min_act = nz
509  do i=is,ie
510  if (do_i(i) .and. (kb(i) < kb_min_act)) kb_min_act = kb(i)
511  enddo
512  ! Solve for the entrainment rate from above in the topmost interior
513  ! layer, eakb, such that
514  ! eakb*dS_implicit = dt*Kd*dS_layer_implicit / h_implicit.
515  do i=is1,ie1
516  ea_kbp1(i) = 0.0
517  if (do_i(i) .and. (kb(i) < nz)) &
518  ea_kbp1(i) = dsp1_ds(i,kb(i)+1)*f(i,kb(i)+1)
519  enddo
520  call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, min_eakb, &
521  max_eakb, kmb, is1, ie1, do_i, g, gv, cs, eakb, f_kb=f_kb, &
522  err_max_eakb0=err_max_eakb0, err_min_eakb0=err_min_eakb0, &
523  dfdfm_kb=dfdfm_kb)
524  else
525  kb_min_act = kb_min
526  endif
527 
528  do it=0,cs%max_ent_it-1
529  do i=is1,ie1 ; if (do_i(i)) then
530  if (.not.cs%bulkmixedlayer) f(i,1) = min(f(i,1),maxf(i,1))
531  b1(i) = 1.0
532  endif ; enddo ! end of i loop
533 
534  ! F_kb has already been found for this iteration, either above or at
535  ! the end of the code for the previous iteration.
536  do k=kb_min_act,nz-1 ; do i=is1,ie1 ; if (do_i(i) .and. (k>=kb(i))) then
537  ! Calculate the flux in layer k.
538  if (cs%bulkmixedlayer .and. (k==kb(i))) then
539  f(i,k) = f_kb(i)
540  dfdfm(i,k) = dfdfm_kb(i)
541  else ! k > kb(i)
542  fprev(i,k) = f(i,k)
543  fm = (f(i,k-1) - h(i,j,k)) + dsp1_ds(i,k+1)*f(i,k+1)
544  fk = grats(i,k)*dtkd(i,k)
545  fr = sqrt(fm*fm + fk)
546 
547  if (fm>=0) then
548  f(i,k) = min(maxf(i,k), i2p2dsp1_ds(i,k) * (fm+fr))
549  else
550  f(i,k) = min(maxf(i,k), i2p2dsp1_ds(i,k) * (fk / (-1.0*fm+fr)))
551  endif
552 
553  if ((f(i,k) >= maxf(i,k)) .or. (fr == 0.0)) then
554  dfdfm(i,k) = 0.0
555  else
556  dfdfm(i,k) = i2p2dsp1_ds(i,k) * ((fr + fm) / fr)
557  endif
558 
559  if (k > k2) then
560  ! This is part of a tridiagonal solver for the actual flux.
561  c1(i,k) = dfdfm(i,k-1)*(dsp1_ds(i,k)*b1(i))
562  b1(i) = 1.0 / (1.0 - c1(i,k)*dfdfm(i,k))
563  f(i,k) = min(b1(i)*(f(i,k)-fprev(i,k)) + fprev(i,k), maxf(i,k))
564  if (f(i,k) >= maxf(i,k)) dfdfm(i,k) = 0.0
565  endif
566  endif
567  endif ; enddo ; enddo
568 
569  do k=nz-2,kb_min_act,-1 ; do i=is1,ie1
570  if (do_i(i) .and. (k > kb(i))) &
571  f(i,k) = min((f(i,k)+c1(i,k+1)*(f(i,k+1)-fprev(i,k+1))),maxf(i,k))
572  enddo ; enddo
573 
574  if (cs%bulkmixedlayer) then
575  do i=is1,ie1
576  if (do_i(i) .and. (kb(i) < nz)) then
577  ! F will be increased to minF later.
578  ea_kbp1(i) = dsp1_ds(i,kb(i)+1)*max(f(i,kb(i)+1), minf(i,kb(i)+1))
579  else
580  ea_kbp1(i) = 0.0
581  endif
582  enddo
583  call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, min_eakb, &
584  max_eakb, kmb, is1, ie1, do_i, g, gv, cs, eakb, f_kb=f_kb, &
585  err_max_eakb0=err_max_eakb0, err_min_eakb0=err_min_eakb0, &
586  dfdfm_kb=dfdfm_kb)
587  do i=is1,ie1
588  if (do_i(i) .and. (kb(i) < nz)) f(i,kb(i)) = f_kb(i)
589  enddo
590  endif
591 
592 ! Determine whether to do another iteration.
593  if (it < cs%max_ent_it-1) then
594 
595  reiterate = .false.
596  if (cs%bulkmixedlayer) then ; do i=is1,ie1 ; if (do_i(i)) then
597  eb_kmb(i) = max(2.0*ent_bl(i,kmb+1) - eakb(i), 0.0)
598  endif ; enddo ; endif
599  do i=is1,ie1
600  did_i(i) = do_i(i) ; do_i(i) = .false.
601  enddo
602  do k=kb_min_act,nz-1 ; do i=is1,ie1
603  if (did_i(i) .and. (k >= kb(i))) then
604  if (f(i,k) < minf(i,k)) then
605  f(i,k) = minf(i,k)
606  do_i(i) = .true. ; reiterate = .true.
607  elseif (k > kb(i)) then
608  if ((abs(f(i,k) - fprev(i,k)) > tolerance) .or. &
609  ((h(i,j,k) + ((1.0+dsp1_ds(i,k))*f(i,k) - &
610  (f(i,k-1) + dsp1_ds(i,k+1)*f(i,k+1)))) < 0.9*angstrom)) then
611  do_i(i) = .true. ; reiterate = .true.
612  endif
613  else ! (k == kb(i))
614 ! A more complicated test is required for the layer beneath the buffer layer,
615 ! since its flux may be partially used to entrain directly from the mixed layer.
616 ! Negative fluxes should not occur with the bulk mixed layer.
617  if (h(i,j,k) + ((f(i,k) + eakb(i)) - &
618  (eb_kmb(i) + dsp1_ds(i,k+1)*f(i,k+1))) < 0.9*angstrom) then
619  do_i(i) = .true. ; reiterate = .true.
620  endif
621  endif
622  endif
623  enddo ; enddo
624  if (.not.reiterate) exit
625  endif ! end of if (it < CS%max_ent_it-1)
626  enddo ! end of it loop
627 ! This is the end of the section that might be iterated.
628 
629 
630  if (it == (cs%max_ent_it)) then
631  ! Limit the flux so that the layer below is not depleted.
632  ! This should only be applied to the last iteration.
633  do i=is1,ie1 ; if (do_i(i)) then
634  f(i,nz-1) = max(f(i,nz-1), min(minf(i,nz-1), 0.0))
635  if (kb(i) >= nz-1) then ; ea_kbp1(i) = 0.0 ; endif
636  endif ; enddo
637  do k=nz-2,kb_min_act,-1 ; do i=is1,ie1 ; if (do_i(i)) then
638  if (k>kb(i)) then
639  f(i,k) = min(max(minf(i,k),f(i,k)), (dsp1_ds(i,k+1)*f(i,k+1) + &
640  max(((f(i,k+1)-dsp1_ds(i,k+2)*f(i,k+2)) + &
641  (h(i,j,k+1) - angstrom)), 0.5*(h(i,j,k+1)-angstrom))))
642  elseif (k==kb(i)) then
643  ea_kbp1(i) = dsp1_ds(i,k+1)*f(i,k+1)
644  h_avail = dsp1_ds(i,k+1)*f(i,k+1) + max(0.5*(h(i,j,k+1)-angstrom), &
645  ((f(i,k+1)-dsp1_ds(i,k+2)*f(i,k+2)) + (h(i,j,k+1) - angstrom)))
646  if ((f(i,k) > 0.0) .and. (f(i,k) > h_avail)) then
647  f_kb(i) = max(0.0, h_avail)
648  f(i,k) = f_kb(i)
649  if ((f_kb(i) < maxf_kb(i)) .and. (eakb_maxf(i) <= eakb(i))) &
650  eakb(i) = eakb_maxf(i)
651  call f_kb_to_ea_kb(h_bl, sref, ent_bl, i_dskbp1, f_kb, kmb, i, &
652  g, gv, cs, eakb)
653  endif
654  endif
655  endif ; enddo ; enddo
656 
657 
658  if (cs%bulkmixedlayer) then ; do i=is1,ie1
659  if (do_i(i) .and. (kb(i) < nz)) then
660  h_avail = eakb(i) + max(0.5*(h_bl(i,kmb+1)-angstrom), &
661  (f_kb(i)-ea_kbp1(i)) + (h_bl(i,kmb+1)-angstrom))
662  ! Ensure that 0 < eb_kmb < h_avail.
663  ent_bl(i,kmb+1) = min(ent_bl(i,kmb+1),0.5*(eakb(i) + h_avail))
664 
665  eb_kmb(i) = max(2.0*ent_bl(i,kmb+1) - eakb(i), 0.0)
666  endif
667  enddo ; endif
668 
669  ! Limit the flux so that the layer above is not depleted.
670  do k=kb_min_act+1,nz-1 ; do i=is1,ie1 ; if (do_i(i)) then
671  if ((.not.cs%bulkmixedlayer) .or. (k > kb(i)+1)) then
672  f(i,k) = min(f(i,k), ds_dsp1(i,k)*( ((f(i,k-1) + &
673  dsp1_ds(i,k-1)*f(i,k-1)) - f(i,k-2)) + (h(i,j,k-1) - angstrom)))
674  f(i,k) = max(f(i,k),min(minf(i,k),0.0))
675  elseif (k == kb(i)+1) then
676  f(i,k) = min(f(i,k), ds_dsp1(i,k)*( ((f(i,k-1) + eakb(i)) - &
677  eb_kmb(i)) + (h(i,j,k-1) - angstrom)))
678  f(i,k) = max(f(i,k),min(minf(i,k),0.0))
679  endif
680  endif ; enddo ; enddo
681  endif ! (it == (CS%max_ent_it))
682 
683  call f_to_ent(f, h, kb, kmb, j, g, gv, cs, dsp1_ds, eakb, ent_bl, ea, eb)
684 
685  ! Calculate the layer thicknesses after the entrainment to constrain the
686  ! corrective fluxes.
687  if (correct_density) then
688  do i=is,ie
689  h_guess(i,1) = (h(i,j,1) - angstrom) + (eb(i,j,1) - ea(i,j,2))
690  h_guess(i,nz) = (h(i,j,nz) - angstrom) + (ea(i,j,nz) - eb(i,j,nz-1))
691  if (h_guess(i,1) < 0.0) h_guess(i,1) = 0.0
692  if (h_guess(i,nz) < 0.0) h_guess(i,nz) = 0.0
693  enddo
694  do k=2,nz-1 ; do i=is,ie
695  h_guess(i,k) = (h(i,j,k) - angstrom) + ((ea(i,j,k) - eb(i,j,k-1)) + &
696  (eb(i,j,k) - ea(i,j,k+1)))
697  if (h_guess(i,k) < 0.0) h_guess(i,k) = 0.0
698  enddo ; enddo
699  if (cs%bulkmixedlayer) then
700  call determine_dskb(h_bl, sref, ent_bl, eakb, is, ie, kmb, g, gv, &
701  .true., ds_kb, ds_anom_lim=ds_anom_lim)
702  do k=nz-1,kb_min,-1
703  call calculate_density(tv%T(is:ie,j,k), tv%S(is:ie,j,k), pres(is:ie), &
704  rcv(is:ie), 1, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
705  do i=is,ie
706  if ((k>kb(i)) .and. (f(i,k) > 0.0)) then
707  ! Within a time step, a layer may entrain no more than its
708  ! thickness for correction. This limitation should apply
709  ! extremely rarely, but precludes undesirable behavior.
710  ! Note: Corrected a sign/logic error & factor of 2 error, and
711  ! the layers tracked the target density better, mostly due to
712  ! the factor of 2 error.
713  f_cor = h(i,j,k) * min(1.0 , max(-ds_dsp1(i,k), &
714  (gv%Rlay(k) - rcv(i)) / (gv%Rlay(k+1)-gv%Rlay(k))) )
715 
716  ! Ensure that (1) Entrainments are positive, (2) Corrections in
717  ! a layer cannot deplete the layer itself (very generously), and
718  ! (3) a layer can take no more than a quarter the mass of its
719  ! neighbor.
720  if (f_cor > 0.0) then
721  f_cor = min(f_cor, 0.9*f(i,k), ds_dsp1(i,k)*0.5*h_guess(i,k), &
722  0.25*h_guess(i,k+1))
723  else
724  f_cor = -min(-f_cor, 0.9*f(i,k), 0.5*h_guess(i,k), &
725  0.25*ds_dsp1(i,k)*h_guess(i,k-1) )
726  endif
727 
728  ea(i,j,k) = ea(i,j,k) - dsp1_ds(i,k)*f_cor
729  eb(i,j,k) = eb(i,j,k) + f_cor
730  elseif ((k==kb(i)) .and. (f(i,k) > 0.0)) then
731  ! Rho_cor is the density anomaly that needs to be corrected,
732  ! taking into account that the true potential density of the
733  ! deepest buffer layer is not exactly what is returned as dS_kb.
734  ds_kb_eff = 2.0*ds_kb(i) - ds_anom_lim(i) ! Could be negative!!!
735  rho_cor = h(i,j,k) * (gv%Rlay(k)-rcv(i)) + eakb(i)*ds_anom_lim(i)
736 
737  ! Ensure that -.9*eakb < ea_cor < .9*eakb
738  if (abs(rho_cor) < abs(0.9*eakb(i)*ds_kb_eff)) then
739  ea_cor = -rho_cor / ds_kb_eff
740  else
741  ea_cor = sign(0.9*eakb(i),-rho_cor*ds_kb_eff)
742  endif
743 
744  if (ea_cor > 0.0) then
745  ! Ensure that -F_cor < 0.5*h_guess
746  ea_cor = min(ea_cor, 0.5*(max_eakb(i) - eakb(i)), &
747  0.5*h_guess(i,k) / (ds_kb(i) * i_dskbp1(i)))
748  else
749  ! Ensure that -ea_cor < 0.5*h_guess & F_cor < 0.25*h_guess(k+1)
750  ea_cor = -min(-ea_cor, 0.5*h_guess(i,k), &
751  0.25*h_guess(i,k+1) / (ds_kb(i) * i_dskbp1(i)))
752  endif
753 
754  ea(i,j,k) = ea(i,j,k) + ea_cor
755  eb(i,j,k) = eb(i,j,k) - (ds_kb(i) * i_dskbp1(i)) * ea_cor
756  elseif (k < kb(i)) then
757  ! Repetative, unless ea(kb) has been corrected.
758  ea(i,j,k) = ea(i,j,k+1)
759  endif
760  enddo
761  enddo
762  do k=kb_min-1,k2,-1 ; do i=is,ie
763  ea(i,j,k) = ea(i,j,k+1)
764  enddo ; enddo
765 
766  ! Repetative, unless ea(kb) has been corrected.
767  k=kmb
768  do i=is,ie
769  ! Do not adjust eb through the base of the buffer layers, but it
770  ! may be necessary to change entrainment from above.
771  h1 = (h(i,j,k) - angstrom) + (eb(i,j,k) - ea(i,j,k+1))
772  ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)
773  enddo
774  do k=kmb-1,2,-1 ; do i=is,ie
775  ! Determine the entrainment from below for each buffer layer.
776  eb(i,j,k) = max(2.0*ent_bl(i,k+1) - ea(i,j,k+1), 0.0)
777 
778  ! Determine the entrainment from above for each buffer layer.
779  h1 = (h(i,j,k) - angstrom) + (eb(i,j,k) - ea(i,j,k+1))
780  ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)
781  enddo ; enddo
782  do i=is,ie
783  eb(i,j,1) = max(2.0*ent_bl(i,2) - ea(i,j,2), 0.0)
784  enddo
785 
786  else ! not bulkmixedlayer
787  do k=k2,nz-1
788  call calculate_density(tv%T(is:ie,j,k), tv%S(is:ie,j,k), pres(is:ie), &
789  rcv(is:ie), 1, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
790  do i=is,ie ; if (f(i,k) > 0.0) then
791  ! Within a time step, a layer may entrain no more than
792  ! its thickness for correction. This limitation should
793  ! apply extremely rarely, but precludes undesirable
794  ! behavior.
795  f_cor = h(i,j,k) * min(dsp1_ds(i,k) , max(-1.0, &
796  (gv%Rlay(k) - rcv(i)) / (gv%Rlay(k+1)-gv%Rlay(k))) )
797 
798  ! Ensure that (1) Entrainments are positive, (2) Corrections in
799  ! a layer cannot deplete the layer itself (very generously), and
800  ! (3) a layer can take no more than a quarter the mass of its
801  ! neighbor.
802  if (f_cor >= 0.0) then
803  f_cor = min(f_cor, 0.9*f(i,k), 0.5*dsp1_ds(i,k)*h_guess(i,k), &
804  0.25*h_guess(i,k+1))
805  else
806  f_cor = -min(-f_cor, 0.9*f(i,k), 0.5*h_guess(i,k), &
807  0.25*ds_dsp1(i,k)*h_guess(i,k-1) )
808  endif
809 
810  ea(i,j,k) = ea(i,j,k) - dsp1_ds(i,k)*f_cor
811  eb(i,j,k) = eb(i,j,k) + f_cor
812  endif ; enddo
813  enddo
814  endif
815 
816  endif ! correct_density
817 
818  if (cs%id_Kd > 0) then
819  idt = gv%H_to_Z**2 / dt
820  do k=2,nz-1 ; do i=is,ie
821  if (k<kb(i)) then ; kd_here = 0.0 ; else
822  kd_here = f(i,k) * ( h(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + &
823  (eb(i,j,k) - ea(i,j,k+1))) ) / (i2p2dsp1_ds(i,k) * grats(i,k))
824  endif
825 
826  kd_eff(i,j,k) = max(dtkd(i,k), kd_here)*idt
827  enddo ; enddo
828  do i=is,ie
829  kd_eff(i,j,1) = dtkd(i,1)*idt
830  kd_eff(i,j,nz) = dtkd(i,nz)*idt
831  enddo
832  endif
833 
834  if (cs%id_diff_work > 0) then
835  g_2dt = 0.5 * gv%H_to_Z**2*us%L_to_Z**2 * (gv%g_Earth / dt)
836  do i=is,ie ; diff_work(i,j,1) = 0.0 ; diff_work(i,j,nz+1) = 0.0 ; enddo
837  if (associated(tv%eqn_of_state)) then
838  if (associated(fluxes%p_surf)) then
839  do i=is,ie ; pressure(i) = fluxes%p_surf(i,j) ; enddo
840  else
841  do i=is,ie ; pressure(i) = 0.0 ; enddo
842  endif
843  do k=2,nz
844  do i=is,ie ; pressure(i) = pressure(i) + gv%H_to_Pa*h(i,j,k-1) ; enddo
845  do i=is,ie
846  if (k==kb(i)) then
847  t_eos(i) = 0.5*(tv%T(i,j,kmb) + tv%T(i,j,k))
848  s_eos(i) = 0.5*(tv%S(i,j,kmb) + tv%S(i,j,k))
849  else
850  t_eos(i) = 0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))
851  s_eos(i) = 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
852  endif
853  enddo
854  call calculate_density_derivs(t_eos, s_eos, pressure, &
855  drho_dt, drho_ds, is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
856  do i=is,ie
857  if ((k>kmb) .and. (k<kb(i))) then ; diff_work(i,j,k) = 0.0
858  else
859  if (k==kb(i)) then
860  drho = drho_dt(i) * (tv%T(i,j,k)-tv%T(i,j,kmb)) + &
861  drho_ds(i) * (tv%S(i,j,k)-tv%S(i,j,kmb))
862  else
863  drho = drho_dt(i) * (tv%T(i,j,k)-tv%T(i,j,k-1)) + &
864  drho_ds(i) * (tv%S(i,j,k)-tv%S(i,j,k-1))
865  endif
866  diff_work(i,j,k) = g_2dt * drho * &
867  (ea(i,j,k) * (h(i,j,k) + ea(i,j,k)) + &
868  eb(i,j,k-1)*(h(i,j,k-1) + eb(i,j,k-1)))
869  endif
870  enddo
871  enddo
872  else
873  do k=2,nz ; do i=is,ie
874  diff_work(i,j,k) = g_2dt * (gv%Rlay(k)-gv%Rlay(k-1)) * &
875  (ea(i,j,k) * (h(i,j,k) + ea(i,j,k)) + &
876  eb(i,j,k-1)*(h(i,j,k-1) + eb(i,j,k-1)))
877  enddo ; enddo
878  endif
879  endif
880 
881  if (present(kb_out)) then
882  do i=is,ie ; kb_out(i,j) = kb(i) ; enddo
883  endif
884 
885  enddo ! end of j loop
886 
887 ! Offer diagnostic fields for averaging.
888  if (cs%id_Kd > 0) call post_data(cs%id_Kd, kd_eff, cs%diag)
889  if (cs%id_Kd > 0) deallocate(kd_eff)
890  if (cs%id_diff_work > 0) call post_data(cs%id_diff_work, diff_work, cs%diag)
891  if (cs%id_diff_work > 0) deallocate(diff_work)
892 

References determine_dskb(), determine_ea_kb(), f_kb_to_ea_kb(), f_to_ent(), find_maxf_kb(), mom_error_handler::is_root_pe(), mom_error_handler::mom_error(), and set_ent_bl().

Here is the call graph for this function:

◆ f_kb_to_ea_kb()

subroutine mom_entrain_diffusive::f_kb_to_ea_kb ( real, dimension(szi_(g),szk_(g)), intent(in)  h_bl,
real, dimension(szi_(g),szk_(g)), intent(in)  Sref,
real, dimension(szi_(g),szk_(g)), intent(in)  Ent_bl,
real, dimension(szi_(g)), intent(in)  I_dSkbp1,
real, dimension(szi_(g)), intent(in)  F_kb,
integer, intent(in)  kmb,
integer, intent(in)  i,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(entrain_diffusive_cs), pointer  CS,
real, dimension( g %isd: g %ied), intent(inout)  ea_kb,
real, intent(in), optional  tol_in 
)
private

Given an entrainment from below for layer kb, determine a consistent entrainment from above, such that dSkb * ea_kb = dSkbp1 * F_kb. The input value of ea_kb is both the maximum value that can be obtained and the first guess of the iterations. Ideally ea_kb should be an under-estimate.

Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]h_blLayer thickness, with the top interior
[in]srefThe coordinate reference potential density,
[in]ent_blThe average entrainment upward and downward
[in]i_dskbp1The inverse of the difference in reference potential density across the base of the uppermost interior layer [R-1 ~> m3 kg-1].
[in]f_kbThe entrainment from below by the uppermost interior layer [H ~> m or kg m-2]
[in]kmbThe number of mixed and buffer layers.
[in]iThe i-index to work on
csThis module's control structure.
[in,out]ea_kbThe entrainment from above by the layer below the buffer layer (i.e. layer kb) [H ~> m or kg m-2].
[in]tol_inA tolerance for the iterative determination of the entrainment [H ~> m or kg m-2].

Definition at line 1444 of file MOM_entrain_diffusive.F90.

1444  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1445  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1446  real, dimension(SZI_(G),SZK_(G)), &
1447  intent(in) :: h_bl !< Layer thickness, with the top interior
1448  !! layer at k-index kmb+1 [H ~> m or kg m-2].
1449  real, dimension(SZI_(G),SZK_(G)), &
1450  intent(in) :: Sref !< The coordinate reference potential density,
1451  !! with the value of the topmost interior layer
1452  !! at index kmb+1 [R ~> kg m-3].
1453  real, dimension(SZI_(G),SZK_(G)), &
1454  intent(in) :: Ent_bl !< The average entrainment upward and downward
1455  !! across each interface around the buffer layers,
1456  !! [H ~> m or kg m-2].
1457  real, dimension(SZI_(G)), intent(in) :: I_dSkbp1 !< The inverse of the difference in reference
1458  !! potential density across the base of the
1459  !! uppermost interior layer [R-1 ~> m3 kg-1].
1460  real, dimension(SZI_(G)), intent(in) :: F_kb !< The entrainment from below by the
1461  !! uppermost interior layer [H ~> m or kg m-2]
1462  integer, intent(in) :: kmb !< The number of mixed and buffer layers.
1463  integer, intent(in) :: i !< The i-index to work on
1464  type(entrain_diffusive_CS), pointer :: CS !< This module's control structure.
1465  real, dimension(SZI_(G)), intent(inout) :: ea_kb !< The entrainment from above by the layer below
1466  !! the buffer layer (i.e. layer kb) [H ~> m or kg m-2].
1467  real, optional, intent(in) :: tol_in !< A tolerance for the iterative determination
1468  !! of the entrainment [H ~> m or kg m-2].
1469 
1470  real :: max_ea, min_ea
1471  real :: err, err_min, err_max
1472  real :: derr_dea
1473  real :: val, tolerance, tol1
1474  real :: ea_prev
1475  real :: dS_kbp1
1476  logical :: bisect_next, Newton
1477  real, dimension(SZI_(G)) :: dS_kb
1478  real, dimension(SZI_(G)) :: maxF, ent_maxF, zeros
1479  real, dimension(SZI_(G)) :: ddSkb_dE
1480  integer :: it
1481  integer, parameter :: MAXIT = 30
1482 
1483  ds_kbp1 = sref(i,kmb+2) - sref(i,kmb+1)
1484  max_ea = ea_kb(i) ; min_ea = 0.0
1485  val = ds_kbp1 * f_kb(i)
1486  err_min = -val
1487 
1488  tolerance = cs%Tolerance_Ent
1489  if (present(tol_in)) tolerance = tol_in
1490  bisect_next = .true.
1491 
1492  call determine_dskb(h_bl, sref, ent_bl, ea_kb, i, i, kmb, g, gv, .true., &
1493  ds_kb, ddskb_de)
1494 
1495  err = ds_kb(i) * ea_kb(i) - val
1496  derr_dea = ds_kb(i) + ddskb_de(i) * ea_kb(i)
1497  ! Return if Newton's method on the first guess would give a tolerably small
1498  ! change in the value of ea_kb.
1499  if ((err <= 0.0) .and. (abs(err) <= tolerance*abs(derr_dea))) return
1500 
1501  if (err == 0.0) then ; return ! The exact solution on the first guess...
1502  elseif (err > 0.0) then ! The root is properly bracketed.
1503  max_ea = ea_kb(i) ; err_max = err
1504  ! Use Newton's method (if it stays bounded) or the false position method
1505  ! to find the next value.
1506  if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i) - min_ea) > err) .and. &
1507  (derr_dea*(max_ea - ea_kb(i)) > -1.0*err)) then
1508  ea_kb(i) = ea_kb(i) - err / derr_dea
1509  else ! Use the bisection for the next guess.
1510  ea_kb(i) = 0.5*(max_ea+min_ea)
1511  endif
1512  else
1513  ! Try to bracket the root first. If unable to bracket the root, return
1514  ! the maximum.
1515  zeros(i) = 0.0
1516  call find_maxf_kb(h_bl, sref, ent_bl, i_dskbp1, zeros, ea_kb, &
1517  kmb, i, i, g, gv, cs, maxf, ent_maxf, f_thresh = f_kb)
1518  err_max = ds_kbp1 * maxf(i) - val
1519  ! If err_max is negative, there is no good solution, so use the maximum
1520  ! value of F in the valid range.
1521  if (err_max <= 0.0) then
1522  ea_kb(i) = ent_maxf(i) ; return
1523  else
1524  max_ea = ent_maxf(i)
1525  ea_kb(i) = 0.5*(max_ea+min_ea) ! Use bisection for the next guess.
1526  endif
1527  endif
1528 
1529  ! Exit if the range between max_ea and min_ea already acceptable.
1530  ! if (abs(max_ea - min_ea) < 0.1*tolerance) return
1531 
1532  do it = 1, maxit
1533  call determine_dskb(h_bl, sref, ent_bl, ea_kb, i, i, kmb, g, gv, .true., &
1534  ds_kb, ddskb_de)
1535 
1536  err = ds_kb(i) * ea_kb(i) - val
1537  derr_dea = ds_kb(i) + ddskb_de(i) * ea_kb(i)
1538 
1539  ea_prev = ea_kb(i)
1540  ! Use Newton's method or the false position method to find the next value.
1541  newton = .false.
1542  if (err > 0.0) then
1543  max_ea = ea_kb(i) ; err_max = err
1544  if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i)-min_ea) > err)) newton = .true.
1545  else
1546  min_ea = ea_kb(i) ; err_min = err
1547  if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i)-max_ea) < err)) newton = .true.
1548  endif
1549 
1550  if (newton) then
1551  ea_kb(i) = ea_kb(i) - err / derr_dea
1552  elseif (bisect_next) then ! Use bisection to reduce the range.
1553  ea_kb(i) = 0.5*(max_ea+min_ea)
1554  bisect_next = .false.
1555  else ! Use the false-position method for the next guess.
1556  ea_kb(i) = min_ea + (max_ea-min_ea) * (err_min/(err_min - err_max))
1557  bisect_next = .true.
1558  endif
1559 
1560  tol1 = tolerance ; if (err > 0.0) tol1 = 0.099*tolerance
1561  if (ds_kb(i) <= ds_kbp1) then
1562  if (abs(ea_kb(i) - ea_prev) <= tol1) return
1563  else
1564  if (ds_kbp1*abs(ea_kb(i) - ea_prev) <= ds_kb(i)*tol1) return
1565  endif
1566  enddo
1567 

References determine_dskb(), and find_maxf_kb().

Referenced by entrainment_diffusive().

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

◆ f_to_ent()

subroutine mom_entrain_diffusive::f_to_ent ( real, dimension(szi_(g),szk_(g)), intent(in)  F,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
integer, dimension(szi_(g)), intent(in)  kb,
integer, intent(in)  kmb,
integer, intent(in)  j,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(entrain_diffusive_cs), intent(in)  CS,
real, dimension( g %isd: g %ied, g %ke), intent(in)  dsp1_ds,
real, dimension( g %isd: g %ied), intent(in)  eakb,
real, dimension( g %isd: g %ied, g %ke), intent(in)  Ent_bl,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  ea,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  eb,
logical, dimension( g %isd: g %ied), intent(in), optional  do_i_in 
)
private

This subroutine calculates the actual entrainments (ea and eb) and the amount of surface forcing that is applied to each layer if there is no bulk mixed layer.

Parameters
[in]gThe ocean's grid structure
[in]gvThe ocean's vertical grid structure
[in]fThe density flux through a layer within a time step divided by the density difference across the interface below the layer [H ~> m or kg m-2].
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]kbThe index of the lightest layer denser than the deepest buffer layer.
[in]kmbThe number of mixed and buffer layers.
[in]jThe meridional index upon which to work.
[in]csThis module's control structure.
[in]dsp1_dsThe ratio of coordinate variable differences across the interfaces below a layer over the difference across the interface above the layer.
[in]eakbThe entrainment from above by the layer below the buffer layer [H ~> m or kg m-2].
[in]ent_blThe average entrainment upward and downward across each interface around the buffer layers [H ~> m or kg m-2].
[in,out]eaThe amount of fluid entrained from the layer
[in,out]ebThe amount of fluid entrained from the layer
[in]do_i_inIndicates which i-points to work on.

Definition at line 899 of file MOM_entrain_diffusive.F90.

899  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
900  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
901  real, dimension(SZI_(G),SZK_(G)), intent(in) :: F !< The density flux through a layer within
902  !! a time step divided by the density
903  !! difference across the interface below
904  !! the layer [H ~> m or kg m-2].
905  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
906  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
907  integer, dimension(SZI_(G)), intent(in) :: kb !< The index of the lightest layer denser than
908  !! the deepest buffer layer.
909  integer, intent(in) :: kmb !< The number of mixed and buffer layers.
910  integer, intent(in) :: j !< The meridional index upon which to work.
911  type(entrain_diffusive_CS), intent(in) :: CS !< This module's control structure.
912  real, dimension(SZI_(G),SZK_(G)), intent(in) :: dsp1_ds !< The ratio of coordinate variable
913  !! differences across the interfaces below
914  !! a layer over the difference across the
915  !! interface above the layer.
916  real, dimension(SZI_(G)), intent(in) :: eakb !< The entrainment from above by the layer
917  !! below the buffer layer [H ~> m or kg m-2].
918  real, dimension(SZI_(G),SZK_(G)), intent(in) :: Ent_bl !< The average entrainment upward and
919  !! downward across each interface around
920  !! the buffer layers [H ~> m or kg m-2].
921  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
922  intent(inout) :: ea !< The amount of fluid entrained from the layer
923  !! above within this time step [H ~> m or kg m-2].
924  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
925  intent(inout) :: eb !< The amount of fluid entrained from the layer
926  !! below within this time step [H ~> m or kg m-2].
927  logical, dimension(SZI_(G)), &
928  optional, intent(in) :: do_i_in !< Indicates which i-points to work on.
929 ! This subroutine calculates the actual entrainments (ea and eb) and the
930 ! amount of surface forcing that is applied to each layer if there is no bulk
931 ! mixed layer.
932 
933  real :: h1 ! The thickness in excess of the minimum that will remain
934  ! after exchange with the layer below [H ~> m or kg m-2].
935  logical :: do_i(SZI_(G))
936  integer :: i, k, is, ie, nz
937 
938  is = g%isc ; ie = g%iec ; nz = g%ke
939 
940  if (present(do_i_in)) then
941  do i=is,ie ; do_i(i) = do_i_in(i) ; enddo
942  do i=g%isc,g%iec ; if (do_i(i)) then
943  is = i ; exit
944  endif ; enddo
945  do i=g%iec,g%isc,-1 ; if (do_i(i)) then
946  ie = i ; exit
947  endif ; enddo
948  else
949  do i=is,ie ; do_i(i) = .true. ; enddo
950  endif
951 
952  do i=is,ie
953  ea(i,j,nz) = 0.0 ; eb(i,j,nz) = 0.0
954  enddo
955  if (cs%bulkmixedlayer) then
956  do i=is,ie
957  eb(i,j,kmb) = max(2.0*ent_bl(i,kmb+1) - eakb(i), 0.0)
958  enddo
959  do k=nz-1,kmb+1,-1 ; do i=is,ie ; if (do_i(i)) then
960  if (k > kb(i)) then
961  ! With a bulk mixed layer, surface buoyancy fluxes are applied
962  ! elsewhere, so F should always be nonnegative.
963  ea(i,j,k) = dsp1_ds(i,k)*f(i,k)
964  eb(i,j,k) = f(i,k)
965  elseif (k == kb(i)) then
966  ea(i,j,k) = eakb(i)
967  eb(i,j,k) = f(i,k)
968  elseif (k == kb(i)-1) then
969  ea(i,j,k) = ea(i,j,k+1)
970  eb(i,j,k) = eb(i,j,kmb)
971  else
972  ea(i,j,k) = ea(i,j,k+1)
973  ! Add the entrainment of the thin interior layers to eb going
974  ! up into the buffer layer.
975  eb(i,j,k) = eb(i,j,k+1) + max(0.0, h(i,j,k+1) - gv%Angstrom_H)
976  endif
977  endif ; enddo ; enddo
978  k = kmb
979  do i=is,ie ; if (do_i(i)) then
980  ! Adjust the previously calculated entrainment from below by the deepest
981  ! buffer layer to account for entrainment of thin interior layers .
982  if (kb(i) > kmb+1) &
983  eb(i,j,k) = eb(i,j,k+1) + max(0.0, h(i,j,k+1) - gv%Angstrom_H)
984 
985  ! Determine the entrainment from above for each buffer layer.
986  h1 = (h(i,j,k) - gv%Angstrom_H) + (eb(i,j,k) - ea(i,j,k+1))
987  ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)
988  endif ; enddo
989  do k=kmb-1,2,-1 ; do i=is,ie ; if (do_i(i)) then
990  ! Determine the entrainment from below for each buffer layer.
991  eb(i,j,k) = max(2.0*ent_bl(i,k+1) - ea(i,j,k+1), 0.0)
992 
993  ! Determine the entrainment from above for each buffer layer.
994  h1 = (h(i,j,k) - gv%Angstrom_H) + (eb(i,j,k) - ea(i,j,k+1))
995  ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)
996 ! if (h1 >= 0.0) then ; ea(i,j,k) = Ent_bl(i,K)
997 ! elseif (Ent_bl(i,K)+0.5*h1 >= 0.0) then ; ea(i,j,k) = Ent_bl(i,K)-0.5*h1
998 ! else ; ea(i,j,k) = -h1 ; endif
999  endif ; enddo ; enddo
1000  do i=is,ie ; if (do_i(i)) then
1001  eb(i,j,1) = max(2.0*ent_bl(i,2) - ea(i,j,2), 0.0)
1002  ea(i,j,1) = 0.0
1003  endif ; enddo
1004  else ! not BULKMIXEDLAYER
1005  ! Calculate the entrainment by each layer from above and below.
1006  ! Entrainment is always positive, but F may be negative due to
1007  ! surface buoyancy fluxes.
1008  do i=is,ie
1009  ea(i,j,1) = 0.0 ; eb(i,j,1) = max(f(i,1),0.0)
1010  ea(i,j,2) = dsp1_ds(i,2)*f(i,2) - min(f(i,1),0.0)
1011  enddo
1012 
1013  do k=2,nz-1 ; do i=is,ie
1014  eb(i,j,k) = max(f(i,k),0.0)
1015  ea(i,j,k+1) = dsp1_ds(i,k+1)*f(i,k+1) - (f(i,k)-eb(i,j,k))
1016  if (ea(i,j,k+1) < 0.0) then
1017  eb(i,j,k) = eb(i,j,k) - ea(i,j,k+1)
1018  ea(i,j,k+1) = 0.0
1019  endif
1020  enddo ; enddo
1021  endif ! end BULKMIXEDLAYER

Referenced by entrainment_diffusive().

Here is the caller graph for this function:

◆ find_maxf_kb()

subroutine mom_entrain_diffusive::find_maxf_kb ( real, dimension(szi_(g),szk_(g)), intent(in)  h_bl,
real, dimension(szi_(g),szk_(g)), intent(in)  Sref,
real, dimension(szi_(g),szk_(g)), intent(in)  Ent_bl,
real, dimension(szi_(g)), intent(in)  I_dSkbp1,
real, dimension(szi_(g)), intent(in)  min_ent_in,
real, dimension(szi_(g)), intent(in)  max_ent_in,
integer, intent(in)  kmb,
integer, intent(in)  is,
integer, intent(in)  ie,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(entrain_diffusive_cs), pointer  CS,
real, dimension( g %isd: g %ied), intent(out)  maxF,
real, dimension( g %isd: g %ied), intent(out), optional  ent_maxF,
logical, dimension( g %isd: g %ied), intent(in), optional  do_i_in,
real, dimension( g %isd: g %ied), intent(out), optional  F_lim_maxent,
real, dimension( g %isd: g %ied), intent(in), optional  F_thresh 
)
private

Maximize F = ent*ds_kb*I_dSkbp1 in the range min_ent < ent < max_ent.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]h_blLayer thickness [H ~> m or kg m-2]
[in]srefReference potential density [R ~> kg m-3].
[in]ent_blThe average entrainment upward and
[in]i_dskbp1The inverse of the difference in reference potential density across the base of the uppermost interior layer [R-1 ~> m3 kg-1].
[in]min_ent_inThe minimum value of ent to search, [H ~> m or kg m-2].
[in]max_ent_inThe maximum value of ent to search, [H ~> m or kg m-2].
[in]kmbThe number of mixed and buffer layers.
[in]isThe start of the i-index range to work on.
[in]ieThe end of the i-index range to work on.
csThis module's control structure.
[out]maxfThe maximum value of F = ent*ds_kb*I_dSkbp1 found in the range min_ent < ent < max_ent [H ~> m or kg m-2].
[out]ent_maxfThe value of ent at that maximum [H ~> m or kg m-2].
[in]do_i_inA logical array indicating which columns
[out]f_lim_maxentIf present, do not apply the limit in
[in]f_threshIf F_thresh is present, return the first

Definition at line 1789 of file MOM_entrain_diffusive.F90.

1789  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1790  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1791  real, dimension(SZI_(G),SZK_(G)), &
1792  intent(in) :: h_bl !< Layer thickness [H ~> m or kg m-2]
1793  real, dimension(SZI_(G),SZK_(G)), &
1794  intent(in) :: Sref !< Reference potential density [R ~> kg m-3].
1795  real, dimension(SZI_(G),SZK_(G)), &
1796  intent(in) :: Ent_bl !< The average entrainment upward and
1797  !! downward across each interface around
1798  !! the buffer layers [H ~> m or kg m-2].
1799  real, dimension(SZI_(G)), intent(in) :: I_dSkbp1 !< The inverse of the difference in
1800  !! reference potential density across the
1801  !! base of the uppermost interior layer
1802  !! [R-1 ~> m3 kg-1].
1803  real, dimension(SZI_(G)), intent(in) :: min_ent_in !< The minimum value of ent to search,
1804  !! [H ~> m or kg m-2].
1805  real, dimension(SZI_(G)), intent(in) :: max_ent_in !< The maximum value of ent to search,
1806  !! [H ~> m or kg m-2].
1807  integer, intent(in) :: kmb !< The number of mixed and buffer layers.
1808  integer, intent(in) :: is !< The start of the i-index range to work on.
1809  integer, intent(in) :: ie !< The end of the i-index range to work on.
1810  type(entrain_diffusive_CS), pointer :: CS !< This module's control structure.
1811  real, dimension(SZI_(G)), intent(out) :: maxF !< The maximum value of F
1812  !! = ent*ds_kb*I_dSkbp1 found in the range
1813  !! min_ent < ent < max_ent [H ~> m or kg m-2].
1814  real, dimension(SZI_(G)), &
1815  optional, intent(out) :: ent_maxF !< The value of ent at that maximum [H ~> m or kg m-2].
1816  logical, dimension(SZI_(G)), &
1817  optional, intent(in) :: do_i_in !< A logical array indicating which columns
1818  !! to work on.
1819  real, dimension(SZI_(G)), &
1820  optional, intent(out) :: F_lim_maxent !< If present, do not apply the limit in
1821  !! finding the maximum value, but return the
1822  !! limited value at ent=max_ent_in in this
1823  !! array [H ~> m or kg m-2].
1824  real, dimension(SZI_(G)), &
1825  optional, intent(in) :: F_thresh !< If F_thresh is present, return the first
1826  !! value found that has F > F_thresh, or
1827  !! the maximum.
1828 
1829 ! Maximize F = ent*ds_kb*I_dSkbp1 in the range min_ent < ent < max_ent.
1830 ! ds_kb may itself be limited to positive values in determine_dSkb, which gives
1831 ! the prospect of two local maxima in the range - one at max_ent_in with that
1832 ! minimum value of ds_kb, and the other due to the unlimited (potentially
1833 ! negative) value. It is faster to find the true maximum by first finding the
1834 ! unlimited maximum and comparing it to the limited value at max_ent_in.
1835  real, dimension(SZI_(G)) :: &
1836  ent, &
1837  minent, maxent, ent_best, &
1838  F_max_ent_in, &
1839  F_maxent, F_minent, F, F_best, &
1840  dF_dent, dF_dE_max, dF_dE_min, dF_dE_best, &
1841  dS_kb, dS_kb_lim, ddSkb_dE, dS_anom_lim, &
1842  chg_prev, chg_pre_prev
1843  real :: dF_dE_mean, maxslope, minslope
1844  real :: tolerance
1845  real :: ratio_select_end
1846  real :: rat, max_chg, min_chg, chg1, chg2, chg
1847  logical, dimension(SZI_(G)) :: do_i, last_it, need_bracket, may_use_best
1848  logical :: doany, OK1, OK2, bisect, new_min_bound
1849  integer :: i, it, is1, ie1
1850  integer, parameter :: MAXIT = 20
1851 
1852  tolerance = cs%Tolerance_Ent
1853 
1854  if (present(do_i_in)) then
1855  do i=is,ie ; do_i(i) = do_i_in(i) ; enddo
1856  else
1857  do i=is,ie ; do_i(i) = .true. ; enddo
1858  endif
1859 
1860  ! The most likely value is at max_ent.
1861  call determine_dskb(h_bl, sref, ent_bl, max_ent_in, is, ie, kmb, g, gv, .false., &
1862  ds_kb, ddskb_de, ds_anom_lim=ds_anom_lim)
1863  ie1 = is-1 ; doany = .false.
1864  do i=is,ie
1865  ds_kb_lim(i) = ds_kb(i) + ds_anom_lim(i)
1866  f_max_ent_in(i) = max_ent_in(i)*ds_kb_lim(i)*i_dskbp1(i)
1867  maxent(i) = max_ent_in(i) ; minent(i) = min_ent_in(i)
1868  if ((abs(maxent(i) - minent(i)) < tolerance) .or. (.not.do_i(i))) then
1869  f_best(i) = max_ent_in(i)*ds_kb(i)*i_dskbp1(i)
1870  ent_best(i) = max_ent_in(i) ; ent(i) = max_ent_in(i)
1871  do_i(i) = .false.
1872  else
1873  f_maxent(i) = maxent(i) * ds_kb(i) * i_dskbp1(i)
1874  df_de_max(i) = (ds_kb(i) + maxent(i)*ddskb_de(i)) * i_dskbp1(i)
1875  doany = .true. ; last_it(i) = .false. ; need_bracket(i) = .true.
1876  endif
1877  enddo
1878 
1879  if (doany) then
1880  ie1 = is-1 ; do i=is,ie ; if (do_i(i)) ie1 = i ; enddo
1881  do i=ie1,is,-1 ; if (do_i(i)) is1 = i ; enddo
1882  ! Find the value of F and its derivative at min_ent.
1883  call determine_dskb(h_bl, sref, ent_bl, minent, is1, ie1, kmb, g, gv, .false., &
1884  ds_kb, ddskb_de, do_i_in = do_i)
1885  do i=is1,ie1 ; if (do_i(i)) then
1886  f_minent(i) = minent(i) * ds_kb(i) * i_dskbp1(i)
1887  df_de_min(i) = (ds_kb(i) + minent(i)*ddskb_de(i)) * i_dskbp1(i)
1888  endif ; enddo
1889 
1890  ratio_select_end = 0.9
1891  do it=1,maxit
1892  ratio_select_end = 0.5*ratio_select_end
1893  do i=is1,ie1 ; if (do_i(i)) then
1894  if (need_bracket(i)) then
1895  df_de_mean = (f_maxent(i) - f_minent(i)) / (maxent(i) - minent(i))
1896  maxslope = max(df_de_mean, df_de_min(i), df_de_max(i))
1897  minslope = min(df_de_mean, df_de_min(i), df_de_max(i))
1898  if (f_minent(i) >= f_maxent(i)) then
1899  if (df_de_min(i) > 0.0) then ; rat = 0.02 ! A small step should bracket the soln.
1900  elseif (maxslope < ratio_select_end*minslope) then
1901  ! The maximum of F is at minent.
1902  f_best(i) = f_minent(i) ; ent_best(i) = minent(i) ; rat = 0.0
1903  do_i(i) = .false.
1904  else ; rat = 0.382 ; endif ! Use the golden ratio
1905  else
1906  if (df_de_max(i) < 0.0) then ; rat = 0.98 ! A small step should bracket the soln.
1907  elseif (minslope > ratio_select_end*maxslope) then
1908  ! The maximum of F is at maxent.
1909  f_best(i) = f_maxent(i) ; ent_best(i) = maxent(i) ; rat = 1.0
1910  do_i(i) = .false.
1911  else ; rat = 0.618 ; endif ! Use the golden ratio
1912  endif
1913 
1914  if (rat >= 0.0) ent(i) = rat*maxent(i) + (1.0-rat)*minent(i)
1915  if (((maxent(i) - minent(i)) < tolerance) .or. (it==maxit)) &
1916  last_it(i) = .true.
1917  else ! The maximum is bracketed by minent, ent_best, and maxent.
1918  chg1 = 2.0*(maxent(i) - minent(i)) ; chg2 = chg1
1919  if (df_de_best(i) > 0) then
1920  max_chg = maxent(i) - ent_best(i) ; min_chg = 0.0
1921  else
1922  max_chg = 0.0 ; min_chg = minent(i) - ent_best(i) ! < 0
1923  endif
1924  if (max_chg - min_chg < 2.0*tolerance) last_it(i) = .true.
1925  if (df_de_max(i) /= df_de_best(i)) &
1926  chg1 = (maxent(i) - ent_best(i))*df_de_best(i) / &
1927  (df_de_best(i) - df_de_max(i))
1928  if (df_de_min(i) /= df_de_best(i)) &
1929  chg2 = (minent(i) - ent_best(i))*df_de_best(i) / &
1930  (df_de_best(i) - df_de_min(i))
1931  ok1 = ((chg1 < max_chg) .and. (chg1 > min_chg))
1932  ok2 = ((chg2 < max_chg) .and. (chg2 > min_chg))
1933  if (.not.(ok1 .or. ok2)) then ; bisect = .true. ; else
1934  if (ok1 .and. ok2) then ! Take the acceptable smaller change.
1935  chg = chg1 ; if (abs(chg2) < abs(chg1)) chg = chg2
1936  elseif (ok1) then ; chg = chg1
1937  else ; chg = chg2 ; endif
1938  if (abs(chg) > 0.5*abs(chg_pre_prev(i))) then ; bisect = .true.
1939  else ; bisect = .false. ; endif
1940  endif
1941  chg_pre_prev(i) = chg_prev(i)
1942  if (bisect) then
1943  if (df_de_best(i) > 0.0) then
1944  ent(i) = 0.5*(maxent(i) + ent_best(i))
1945  chg_prev(i) = 0.5*(maxent(i) - ent_best(i))
1946  else
1947  ent(i) = 0.5*(minent(i) + ent_best(i))
1948  chg_prev(i) = 0.5*(minent(i) - ent_best(i))
1949  endif
1950  else
1951  if (abs(chg) < tolerance) chg = sign(tolerance,chg)
1952  ent(i) = ent_best(i) + chg
1953  chg_prev(i) = chg
1954  endif
1955  endif
1956  endif ; enddo
1957 
1958  if (mod(it,3) == 0) then ! Re-determine the loop bounds.
1959  ie1 = is-1 ; do i=is1,ie ; if (do_i(i)) ie1 = i ; enddo
1960  do i=ie1,is,-1 ; if (do_i(i)) is1 = i ; enddo
1961  endif
1962 
1963  call determine_dskb(h_bl, sref, ent_bl, ent, is1, ie1, kmb, g, gv, .false., &
1964  ds_kb, ddskb_de, do_i_in = do_i)
1965  do i=is1,ie1 ; if (do_i(i)) then
1966  f(i) = ent(i)*ds_kb(i)*i_dskbp1(i)
1967  df_dent(i) = (ds_kb(i) + ent(i)*ddskb_de(i)) * i_dskbp1(i)
1968  endif ; enddo
1969 
1970  if (present(f_thresh)) then ; do i=is1,ie1 ; if (do_i(i)) then
1971  if (f(i) >= f_thresh(i)) then
1972  f_best(i) = f(i) ; ent_best(i) = ent(i) ; do_i(i) = .false.
1973  endif
1974  endif ; enddo ; endif
1975 
1976  doany = .false.
1977  do i=is1,ie1 ; if (do_i(i)) then
1978  if (.not.last_it(i)) doany = .true.
1979  if (last_it(i)) then
1980  if (need_bracket(i)) then
1981  if ((f(i) > f_maxent(i)) .and. (f(i) > f_minent(i))) then
1982  f_best(i) = f(i) ; ent_best(i) = ent(i)
1983  elseif (f_maxent(i) > f_minent(i)) then
1984  f_best(i) = f_maxent(i) ; ent_best(i) = maxent(i)
1985  else
1986  f_best(i) = f_minent(i) ; ent_best(i) = minent(i)
1987  endif
1988  elseif (f(i) > f_best(i)) then
1989  f_best(i) = f(i) ; ent_best(i) = ent(i)
1990  endif
1991  do_i(i) = .false.
1992  elseif (need_bracket(i)) then
1993  if ((f(i) > f_maxent(i)) .and. (f(i) > f_minent(i))) then
1994  need_bracket(i) = .false. ! The maximum is now bracketed.
1995  chg_prev(i) = (maxent(i) - minent(i))
1996  chg_pre_prev(i) = 2.0*chg_prev(i)
1997  ent_best(i) = ent(i) ; f_best(i) = f(i) ; df_de_best(i) = df_dent(i)
1998  elseif ((f(i) <= f_maxent(i)) .and. (f(i) > f_minent(i))) then
1999  new_min_bound = .true. ! We have a new minimum bound.
2000  elseif ((f(i) <= f_maxent(i)) .and. (f(i) > f_minent(i))) then
2001  new_min_bound = .false. ! We have a new maximum bound.
2002  else ! This case would bracket a minimum. Wierd.
2003  ! Unless the derivative indicates that there is a maximum near the
2004  ! lower bound, try keeping the end with the larger value of F
2005  ! in a tie keep the minimum as the answer here will be compared
2006  ! with the maximum input value later.
2007  new_min_bound = .true.
2008  if (df_de_min(i) > 0.0 .or. (f_minent(i) >= f_maxent(i))) &
2009  new_min_bound = .false.
2010  endif
2011  if (need_bracket(i)) then ! Still not bracketed.
2012  if (new_min_bound) then
2013  minent(i) = ent(i) ; f_minent(i) = f(i) ; df_de_min(i) = df_dent(i)
2014  else
2015  maxent(i) = ent(i) ; f_maxent(i) = f(i) ; df_de_max(i) = df_dent(i)
2016  endif
2017  endif
2018  else ! The root was previously bracketed.
2019  if (f(i) >= f_best(i)) then ! There is a new maximum.
2020  if (ent(i) > ent_best(i)) then ! Replace minent with ent_prev.
2021  minent(i) = ent_best(i) ; f_minent(i) = f_best(i) ; df_de_min(i) = df_de_best(i)
2022  else ! Replace maxent with ent_best.
2023  maxent(i) = ent_best(i) ; f_maxent(i) = f_best(i) ; df_de_max(i) = df_de_best(i)
2024  endif
2025  ent_best(i) = ent(i) ; f_best(i) = f(i) ; df_de_best(i) = df_dent(i)
2026  else
2027  if (ent(i) < ent_best(i)) then ! Replace the minent with ent.
2028  minent(i) = ent(i) ; f_minent(i) = f(i) ; df_de_min(i) = df_dent(i)
2029  else ! Replace maxent with ent_prev.
2030  maxent(i) = ent(i) ; f_maxent(i) = f(i) ; df_de_max(i) = df_dent(i)
2031  endif
2032  endif
2033  if ((maxent(i) - minent(i)) <= tolerance) do_i(i) = .false. ! Done.
2034  endif ! need_bracket.
2035  endif ; enddo
2036  if (.not.doany) exit
2037  enddo
2038  endif
2039 
2040  if (present(f_lim_maxent)) then
2041  ! Return the unlimited maximum in maxF, and the limited value of F at maxent.
2042  do i=is,ie
2043  maxf(i) = f_best(i)
2044  f_lim_maxent(i) = f_max_ent_in(i)
2045  if (present(ent_maxf)) ent_maxf(i) = ent_best(i)
2046  enddo
2047  else
2048  ! Now compare the two? potential maxima using the limited value of dF_kb.
2049  doany = .false.
2050  do i=is,ie
2051  may_use_best(i) = (ent_best(i) /= max_ent_in(i))
2052  if (may_use_best(i)) doany = .true.
2053  enddo
2054  if (doany) then
2055  ! For efficiency, could save previous value of dS_anom_lim_best?
2056  call determine_dskb(h_bl, sref, ent_bl, ent_best, is, ie, kmb, g, gv, .true., &
2057  ds_kb_lim)
2058  do i=is,ie
2059  f_best(i) = ent_best(i)*ds_kb_lim(i)*i_dskbp1(i)
2060  ! The second test seems necessary because of roundoff differences that
2061  ! can arise during compilation.
2062  if ((f_best(i) > f_max_ent_in(i)) .and. (may_use_best(i))) then
2063  maxf(i) = f_best(i)
2064  if (present(ent_maxf)) ent_maxf(i) = ent_best(i)
2065  else
2066  maxf(i) = f_max_ent_in(i)
2067  if (present(ent_maxf)) ent_maxf(i) = max_ent_in(i)
2068  endif
2069  enddo
2070  else
2071  ! All of the maxima are at the maximum entrainment.
2072  do i=is,ie ; maxf(i) = f_max_ent_in(i) ; enddo
2073  if (present(ent_maxf)) then
2074  do i=is,ie ; ent_maxf(i) = max_ent_in(i) ; enddo
2075  endif
2076  endif
2077  endif
2078 

References determine_dskb().

Referenced by entrainment_diffusive(), and f_kb_to_ea_kb().

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

◆ set_ent_bl()

subroutine mom_entrain_diffusive::set_ent_bl ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
real, dimension( g %isd: g %ied, g %ke+1), intent(in)  dtKd_int,
type(thermo_var_ptrs), intent(in)  tv,
integer, dimension( g %isd: g %ied), intent(inout)  kb,
integer, intent(in)  kmb,
logical, dimension( g %isd: g %ied), intent(in)  do_i,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(entrain_diffusive_cs), pointer  CS,
integer, intent(in)  j,
real, dimension(szi_(g),szk_(g)+1), intent(out)  Ent_bl,
real, dimension(szi_(g),szk_(g)), intent(out)  Sref,
real, dimension(szi_(g),szk_(g)), intent(out)  h_bl 
)
private

This subroutine sets the average entrainment across each of the interfaces between buffer layers within a timestep. It also causes thin and relatively light interior layers to be entrained by the deepest buffer layer. Also find the initial coordinate potential densities (Sref) of each layer.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]dtkd_intThe diapycnal diffusivity across
[in]tvA structure containing pointers to any available thermodynamic fields. Absent fields have NULL ptrs.
[in,out]kbThe index of the lightest layer denser than the buffer layer or 1 if there is no buffer layer.
[in]kmbThe number of mixed and buffer layers.
[in]do_iA logical variable indicating which i-points to work on.
[in]usA dimensional unit scaling type
csThis module's control structure.
[in]jThe meridional index upon which to work.
[out]ent_blThe average entrainment upward and
[out]srefThe coordinate potential density minus 1000 for each layer [R ~> kg m-3].
[out]h_blThe thickness of each layer [H ~> m or kg m-2].

Definition at line 1029 of file MOM_entrain_diffusive.F90.

1029  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1030  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1031  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1032  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1033  real, dimension(SZI_(G),SZK_(G)+1), &
1034  intent(in) :: dtKd_int !< The diapycnal diffusivity across
1035  !! each interface times the time step
1036  !! [H2 ~> m2 or kg2 m-4].
1037  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
1038  !! available thermodynamic fields. Absent
1039  !! fields have NULL ptrs.
1040  integer, dimension(SZI_(G)), intent(inout) :: kb !< The index of the lightest layer denser
1041  !! than the buffer layer or 1 if there is
1042  !! no buffer layer.
1043  integer, intent(in) :: kmb !< The number of mixed and buffer layers.
1044  logical, dimension(SZI_(G)), intent(in) :: do_i !< A logical variable indicating which
1045  !! i-points to work on.
1046  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1047  type(entrain_diffusive_CS), pointer :: CS !< This module's control structure.
1048  integer, intent(in) :: j !< The meridional index upon which to work.
1049  real, dimension(SZI_(G),SZK_(G)+1), &
1050  intent(out) :: Ent_bl !< The average entrainment upward and
1051  !! downward across each interface around
1052  !! the buffer layers [H ~> m or kg m-2].
1053  real, dimension(SZI_(G),SZK_(G)), intent(out) :: Sref !< The coordinate potential density minus
1054  !! 1000 for each layer [R ~> kg m-3].
1055  real, dimension(SZI_(G),SZK_(G)), intent(out) :: h_bl !< The thickness of each layer [H ~> m or kg m-2].
1056 
1057 ! This subroutine sets the average entrainment across each of the interfaces
1058 ! between buffer layers within a timestep. It also causes thin and relatively
1059 ! light interior layers to be entrained by the deepest buffer layer.
1060 ! Also find the initial coordinate potential densities (Sref) of each layer.
1061 ! Does there need to be limiting when the layers below are all thin?
1062 
1063  ! Local variables
1064  real, dimension(SZI_(G)) :: &
1065  b1, d1, & ! Variables used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1] and [nondim].
1066  Rcv, & ! Value of the coordinate variable (potential density)
1067  ! based on the simulated T and S and P_Ref [R ~> kg m-3].
1068  pres, & ! Reference pressure (P_Ref) [Pa].
1069  frac_rem, & ! The fraction of the diffusion remaining [nondim].
1070  h_interior ! The interior thickness available for entrainment [H ~> m or kg m-2].
1071  real, dimension(SZI_(G), SZK_(G)) :: &
1072  S_est ! An estimate of the coordinate potential density - 1000 after
1073  ! entrainment for each layer [R ~> kg m-3].
1074  real :: max_ent ! The maximum possible entrainment [H ~> m or kg m-2].
1075  real :: dh ! An available thickness [H ~> m or kg m-2].
1076  real :: Kd_x_dt ! The diffusion that remains after thin layers are
1077  ! entrained [H2 ~> m2 or kg2 m-4].
1078  real :: h_neglect ! A thickness that is so small it is usually lost
1079  ! in roundoff and can be neglected [H ~> m or kg m-2].
1080  integer :: i, k, is, ie, nz
1081  is = g%isc ; ie = g%iec ; nz = g%ke
1082 
1083 ! max_ent = 1.0e14*GV%Angstrom_H ! This is set to avoid roundoff problems.
1084  max_ent = 1.0e4*gv%m_to_H
1085  h_neglect = gv%H_subroundoff
1086 
1087  do i=is,ie ; pres(i) = tv%P_Ref ; enddo
1088  do k=1,kmb
1089  call calculate_density(tv%T(is:ie,j,k), tv%S(is:ie,j,k), pres(is:ie), &
1090  rcv(is:ie), 1, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
1091  do i=is,ie
1092  h_bl(i,k) = h(i,j,k) + h_neglect
1093  sref(i,k) = rcv(i) - cs%Rho_sig_off
1094  enddo
1095  enddo
1096 
1097  do i=is,ie
1098  h_interior(i) = 0.0 ; ent_bl(i,1) = 0.0
1099 ! if (kb(i) > nz) Ent_bl(i,Kmb+1) = 0.0
1100  enddo
1101 
1102  do k=2,kmb ; do i=is,ie
1103  if (do_i(i)) then
1104  ent_bl(i,k) = min(2.0 * dtkd_int(i,k) / (h(i,j,k-1) + h(i,j,k) + h_neglect), &
1105  max_ent)
1106  else ; ent_bl(i,k) = 0.0 ; endif
1107  enddo ; enddo
1108 
1109  ! Determine the coordinate density of the bottommost buffer layer if there
1110  ! is no entrainment from the layers below. This is a partial solver, based
1111  ! on the first pass of a tridiagonal solver, as the values in the upper buffer
1112  ! layers are not needed.
1113 
1114  do i=is,ie
1115  b1(i) = 1.0 / (h_bl(i,1) + ent_bl(i,2))
1116  d1(i) = h_bl(i,1) * b1(i) ! = 1.0 - Ent_bl(i,2)*b1(i)
1117  s_est(i,1) = (h_bl(i,1)*sref(i,1)) * b1(i)
1118  enddo
1119  do k=2,kmb-1 ; do i=is,ie
1120  b1(i) = 1.0 / ((h_bl(i,k) + ent_bl(i,k+1)) + d1(i)*ent_bl(i,k))
1121  d1(i) = (h_bl(i,k) + d1(i)*ent_bl(i,k)) * b1(i) ! = 1.0 - Ent_bl(i,K+1)*b1(i)
1122  s_est(i,k) = (h_bl(i,k)*sref(i,k) + ent_bl(i,k)*s_est(i,k-1)) * b1(i)
1123  enddo ; enddo
1124  do i=is,ie
1125  s_est(i,kmb) = (h_bl(i,kmb)*sref(i,kmb) + ent_bl(i,kmb)*s_est(i,kmb-1)) / &
1126  (h_bl(i,kmb) + d1(i)*ent_bl(i,kmb))
1127  frac_rem(i) = 1.0
1128  enddo
1129 
1130  ! Entrain any thin interior layers that are lighter (in the coordinate
1131  ! potential density) than the deepest buffer layer will be, and adjust kb.
1132  do i=is,ie ; kb(i) = nz+1 ; if (do_i(i)) kb(i) = kmb+1 ; enddo
1133 
1134  do k=kmb+1,nz ; do i=is,ie ; if (do_i(i)) then
1135  if ((k == kb(i)) .and. (s_est(i,kmb) > (gv%Rlay(k) - cs%Rho_sig_off))) then
1136  if (4.0*dtkd_int(i,kmb+1)*frac_rem(i) > &
1137  (h_bl(i,kmb) + h(i,j,k)) * (h(i,j,k) - gv%Angstrom_H)) then
1138  ! Entrain this layer into the buffer layer and move kb down.
1139  dh = max((h(i,j,k) - gv%Angstrom_H), 0.0)
1140  if (dh > 0.0) then
1141  frac_rem(i) = frac_rem(i) - ((h_bl(i,kmb) + h(i,j,k)) * dh) / &
1142  (4.0*dtkd_int(i,kmb+1))
1143  sref(i,kmb) = (h_bl(i,kmb)*sref(i,kmb) + dh*(gv%Rlay(k)-cs%Rho_sig_off)) / &
1144  (h_bl(i,kmb) + dh)
1145  h_bl(i,kmb) = h_bl(i,kmb) + dh
1146  s_est(i,kmb) = (h_bl(i,kmb)*sref(i,kmb) + ent_bl(i,kmb)*s_est(i,kmb-1)) / &
1147  (h_bl(i,kmb) + d1(i)*ent_bl(i,kmb))
1148  endif
1149  kb(i) = kb(i) + 1
1150  endif
1151  endif
1152  endif ; enddo ; enddo
1153 
1154  ! This is where variables are be set up with a different vertical grid
1155  ! in which the (newly?) massless layers are taken out.
1156  do k=nz,kmb+1,-1 ; do i=is,ie
1157  if (k >= kb(i)) h_interior(i) = h_interior(i) + (h(i,j,k)-gv%Angstrom_H)
1158  if (k==kb(i)) then
1159  h_bl(i,kmb+1) = h(i,j,k) ; sref(i,kmb+1) = gv%Rlay(k) - cs%Rho_sig_off
1160  elseif (k==kb(i)+1) then
1161  h_bl(i,kmb+2) = h(i,j,k) ; sref(i,kmb+2) = gv%Rlay(k) - cs%Rho_sig_off
1162  endif
1163  enddo ; enddo
1164  do i=is,ie ; if (kb(i) >= nz) then
1165  h_bl(i,kmb+1) = h(i,j,nz)
1166  sref(i,kmb+1) = gv%Rlay(nz) - cs%Rho_sig_off
1167  h_bl(i,kmb+2) = gv%Angstrom_H
1168  sref(i,kmb+2) = sref(i,kmb+1) + (gv%Rlay(nz) - gv%Rlay(nz-1))
1169  endif ; enddo
1170 
1171  ! Perhaps we should revisit the way that the average entrainment between the
1172  ! buffer layer and the interior is calculated so that it is not unduly
1173  ! limited when the layers are less than sqrt(Kd * dt) thick?
1174  do i=is,ie ; if (do_i(i)) then
1175  kd_x_dt = frac_rem(i) * dtkd_int(i,kmb+1)
1176  if ((kd_x_dt <= 0.0) .or. (h_interior(i) <= 0.0)) then
1177  ent_bl(i,kmb+1) = 0.0
1178  else
1179  ! If the combined layers are exceptionally thin, use sqrt(Kd*dt) as the
1180  ! estimate of the thickness in the denominator of the thickness diffusion.
1181  ent_bl(i,kmb+1) = min(0.5*h_interior(i), sqrt(kd_x_dt), &
1182  kd_x_dt / (0.5*(h_bl(i,kmb) + h_bl(i,kmb+1))))
1183  endif
1184  else
1185  ent_bl(i,kmb+1) = 0.0
1186  endif ; enddo
1187 

Referenced by entrainment_diffusive().

Here is the caller graph for this function: