MOM6
Data Types | Functions/Subroutines
mom_entrain_diffusive Module Reference

Data Types

type  entrain_diffusive_cs
 

Functions/Subroutines

subroutine, public entrainment_diffusive (u, v, h, tv, fluxes, dt, G, GV, 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, 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)
 
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)
 
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)
 
subroutine, public entrain_diffusive_init (Time, G, GV, param_file, diag, CS)
 
subroutine, public entrain_diffusive_end (CS)
 

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, in m or kg m-2 (abbreviated as H below).
[in]srefReference potential vorticity (in kg m-3?).
[in]ent_blThe average entrainment upward and downward across each interface around the buffer layers, in H.
[in]e_kbThe entrainment by the top interior layer, in H.
[in]ieThe range of i-indices 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. dSkb > 0.
[in,out]ddskb_deThe partial derivative of dSkb with E, in kg m-3 H-1.
[in,out]dslayThe limited potential density difference across the topmost interior layer. 0 < dSkb
[in,out]ddslay_deThe partial derivative of dSlay with E, in kg m-3 H-1.
[in]do_i_inIf present, determines which columns are worked on.

Definition at line 1284 of file MOM_entrain_diffusive.F90.

References mom_error_handler::mom_error().

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

1284  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
1285  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid
1286  !! structure.
1287  real, dimension(SZI_(G),SZK_(G)), intent(in) :: h_bl !< Layer thickness, in m or kg m-2
1288  !! (abbreviated as H below).
1289  real, dimension(SZI_(G),SZK_(G)), intent(in) :: sref !< Reference potential vorticity
1290  !! (in kg m-3?).
1291  real, dimension(SZI_(G),SZK_(G)), intent(in) :: ent_bl !< The average entrainment upward and
1292  !! downward across each interface
1293  !! around the buffer layers, in H.
1294  real, dimension(SZI_(G)), intent(in) :: e_kb !< The entrainment by the top interior
1295  !! layer, in H.
1296  integer, intent(in) :: is, ie !< The range of i-indices to work on.
1297  integer, intent(in) :: kmb !< The number of mixed and buffer
1298  !! layers.
1299  logical, intent(in) :: limit !< If true, limit dSkb and dSlay to
1300  !! avoid negative values.
1301  real, dimension(SZI_(G)), intent(inout) :: dskb !< The limited potential density
1302  !! difference across the interface
1303  !! between the bottommost buffer layer
1304  !! and the topmost interior layer.
1305  !! dSkb > 0.
1306  real, dimension(SZI_(G)), optional, intent(inout) :: ddskb_de !< The partial derivative of dSkb
1307  !! with E, in kg m-3 H-1.
1308  real, dimension(SZI_(G)), optional, intent(inout) :: dslay !< The limited potential density
1309  !! difference across the topmost
1310  !! interior layer. 0 < dSkb
1311  real, dimension(SZI_(G)), optional, intent(inout) :: ddslay_de !< The partial derivative of dSlay
1312  !! with E, in kg m-3 H-1.
1313  real, dimension(SZI_(G)), optional, intent(inout) :: ds_anom_lim
1314  logical, dimension(SZI_(G)), optional, intent(in) :: do_i_in !< If present, determines which
1315  !! columns are worked on.
1316 ! Arguments: h_bl - Layer thickness, in m or kg m-2 (abbreviated as H below).
1317 ! (in) Sref - Reference potential vorticity (in kg m-3?)
1318 ! (in) Ent_bl - The average entrainment upward and downward across
1319 ! each interface around the buffer layers, in H.
1320 ! (in) E_kb - The entrainment by the top interior layer, in H.
1321 ! (in) is, ie - The range of i-indices to work on.
1322 ! (in) kmb - The number of mixed and buffer layers.
1323 ! (in) G - The ocean's grid structure.
1324 ! (in) GV - The ocean's vertical grid structure.
1325 ! (in) limit - If true, limit dSkb and dSlay to avoid negative values.
1326 ! (out) dSkb - The limited potential density difference across the
1327 ! interface between the bottommost buffer layer and the
1328 ! topmost interior layer. dSkb > 0.
1329 ! (out,opt) dSlay - The limited potential density difference across the
1330 ! topmost interior layer. 0 < dSkb
1331 ! (out,opt) ddSkb_dE - The partial derivative of dSkb with E, in kg m-3 H-1.
1332 ! (out,opt) ddSlay_dE - The partial derivative of dSlay with E, in kg m-3 H-1.
1333 ! (in,opt) do_i_in - If present, determines which columns are worked on.
1334 ! Note that dSkb, ddSkb_dE, dSlay, ddSlay_dE, and dS_anom_lim are declared
1335 ! intent(inout) because they should not change where do_i_in is false.
1336 
1337 ! This subroutine determines the reference density difference between the
1338 ! bottommost buffer layer and the first interior after the mixing between mixed
1339 ! and buffer layers and mixing with the layer below. Within the mixed and buffer
1340 ! layers, entrainment from the layer above is increased when it is necessary to
1341 ! keep the layers from developing a negative thickness; otherwise it equals
1342 ! Ent_bl. At each interface, the upward and downward fluxes average out to
1343 ! Ent_bl, unless entrainment by the layer below is larger than twice Ent_bl.
1344 ! The density difference across the first interior layer may also be returned.
1345 ! It could also be limited to avoid negative values or values that greatly
1346 ! exceed the density differences across an interface.
1347 ! Additionally, the partial derivatives of dSkb and dSlay with E_kb could
1348 ! also be returned.
1349  real, dimension(SZI_(G),SZK_(G)) :: &
1350  b1, c1, & ! b1 and c1 are variables used by the tridiagonal solver.
1351  s, ds_de, & ! The coordinate density and its derivative with R.
1352  ea, dea_de, & ! The entrainment from above and its derivative with R.
1353  eb, deb_de ! The entrainment from below and its derivative with R.
1354  real :: deriv_dskb(szi_(g))
1355  real :: d1(szi_(g)) ! d1 = 1.0-c1 is also used by the tridiagonal solver.
1356  real :: src ! A source term for dS_dR.
1357  real :: h1 ! The thickness in excess of the minimum that will remain
1358  ! after exchange with the layer below, in m or kg m-2.
1359  logical, dimension(SZI_(G)) :: do_i
1360  real :: h_neglect ! A thickness that is so small it is usually lost
1361  ! in roundoff and can be neglected, in H.
1362  real :: h_tr ! h_tr is h at tracer points with a tiny thickness
1363  ! added to ensure positive definiteness, in m or kg m-2.
1364  real :: b_denom_1 ! The first term in the denominator of b1 in m or kg m-2.
1365  real :: rat
1366  real :: ds_kbp1, ids_kbp1
1367  real :: deriv_dslay
1368  real :: inv_term ! Nondimensional.
1369  real :: f1, df1_drat ! Nondimensional temporary variables.
1370  real :: z, dz_drat, f2, df2_dz, expz ! Nondimensional temporary variables.
1371  real :: eps_dslay, eps_dskb ! Small nondimensional constants.
1372  integer :: i, k
1373 
1374  if (present(ddslay_de) .and. .not.present(dslay)) call mom_error(fatal, &
1375  "In deterimine_dSkb, ddSLay_dE may only be present if dSlay is.")
1376 
1377  h_neglect = gv%H_subroundoff
1378 
1379  do i=is,ie
1380  ea(i,kmb+1) = e_kb(i) ; dea_de(i,kmb+1) = 1.0
1381  s(i,kmb+1) = sref(i,kmb+1) ; ds_de(i,kmb+1) = 0.0
1382  b1(i,kmb+1) = 0.0
1383  d1(i) = 1.0
1384  do_i(i) = .true.
1385  enddo
1386  if (present(do_i_in)) then
1387  do i=is,ie ; do_i(i) = do_i_in(i) ; enddo
1388  endif
1389  do k=kmb,1,-1 ; do i=is,ie
1390  if (do_i(i)) then
1391  ! The do_i test here is only for efficiency.
1392  ! Determine the entrainment from below for each buffer layer.
1393  if (2.0*ent_bl(i,k+1) > ea(i,k+1)) then
1394  eb(i,k) = 2.0*ent_bl(i,k+1) - ea(i,k+1) ; deb_de(i,k) = -dea_de(i,k+1)
1395  else
1396  eb(i,k) = 0.0 ; deb_de(i,k) = 0.0
1397  endif
1398 
1399  ! Determine the entrainment from above for each buffer layer.
1400  h1 = (h_bl(i,k) - gv%Angstrom) + (eb(i,k) - ea(i,k+1))
1401  if (h1 >= 0.0) then
1402  ea(i,k) = ent_bl(i,k) ; dea_de(i,k) = 0.0
1403  elseif (ent_bl(i,k) + 0.5*h1 >= 0.0) then
1404  ea(i,k) = ent_bl(i,k) - 0.5*h1
1405  dea_de(i,k) = 0.5*(dea_de(i,k+1) - deb_de(i,k))
1406  else
1407  ea(i,k) = -h1
1408  dea_de(i,k) = dea_de(i,k+1) - deb_de(i,k)
1409  endif
1410  else
1411  ea(i,k) = 0.0 ; dea_de(i,k) = 0.0 ; eb(i,k) = 0.0 ; deb_de(i,k) = 0.0
1412  endif
1413 
1414  ! This is the first-pass of a tridiagonal solver for S.
1415  h_tr = h_bl(i,k) + h_neglect
1416  c1(i,k) = ea(i,k+1) * b1(i,k+1)
1417  b_denom_1 = (h_tr + d1(i)*eb(i,k))
1418  b1(i,k) = 1.0 / (b_denom_1 + ea(i,k))
1419  d1(i) = b_denom_1 * b1(i,k)
1420 
1421  s(i,k) = (h_tr*sref(i,k) + eb(i,k)*s(i,k+1)) * b1(i,k)
1422  enddo ; enddo
1423  do k=2,kmb ; do i=is,ie
1424  s(i,k) = s(i,k) + c1(i,k-1)*s(i,k-1)
1425  enddo ; enddo
1426 
1427  if (present(ddskb_de) .or. present(ddslay_de)) then
1428  ! These two tridiagonal solvers cannot be combined because the solutions for
1429  ! S are required as a source for dS_dE.
1430  do k=kmb,2,-1 ; do i=is,ie
1431  if (do_i(i) .and. (dea_de(i,k) - deb_de(i,k) > 0.0)) then
1432  src = (((s(i,k+1) - sref(i,k)) * (h_bl(i,k) + h_neglect) + &
1433  (s(i,k+1) - s(i,k-1)) * ea(i,k)) * deb_de(i,k) - &
1434  ((sref(i,k) - s(i,k-1)) * h_bl(i,k) + &
1435  (s(i,k+1) - s(i,k-1)) * eb(i,k)) * dea_de(i,k)) / &
1436  ((h_bl(i,k) + h_neglect + ea(i,k)) + eb(i,k))
1437  else ; src = 0.0 ; endif
1438  ds_de(i,k) = (src + eb(i,k)*ds_de(i,k+1)) * b1(i,k)
1439  enddo ; enddo
1440  do i=is,ie
1441  if (do_i(i) .and. (deb_de(i,1) < 0.0)) then
1442  src = (((s(i,2) - sref(i,1)) * (h_bl(i,1) + h_neglect)) * deb_de(i,1)) / &
1443  (h_bl(i,1) + h_neglect + eb(i,1))
1444  else ; src = 0.0 ; endif
1445  ds_de(i,1) = (src + eb(i,1)*ds_de(i,2)) * b1(i,1)
1446  enddo
1447  do k=2,kmb ; do i=is,ie
1448  ds_de(i,k) = ds_de(i,k) + c1(i,k-1)*ds_de(i,k-1)
1449  enddo ; enddo
1450  endif
1451 
1452  ! Now, apply any limiting and return the requested variables.
1453 
1454  eps_dskb = 1.0e-6 ! Should be a small, nondimensional, positive number.
1455  if (.not.limit) then
1456  do i=is,ie ; if (do_i(i)) then
1457  dskb(i) = sref(i,kmb+1) - s(i,kmb)
1458  endif ; enddo
1459  if (present(ddskb_de)) then ; do i=is,ie ; if (do_i(i)) then
1460  ddskb_de(i) = -1.0*ds_de(i,kmb)
1461  endif ; enddo ; endif
1462 
1463  if (present(dslay)) then ; do i=is,ie ; if (do_i(i)) then
1464  dslay(i) = 0.5 * (sref(i,kmb+2) - s(i,kmb))
1465  endif ; enddo ; endif
1466  if (present(ddslay_de)) then ; do i=is,ie ; if (do_i(i)) then
1467  ddslay_de(i) = -0.5*ds_de(i,kmb)
1468  endif ; enddo ; endif
1469  else
1470  do i=is,ie ; if (do_i(i)) then
1471  ! Need to ensure that 0 < dSkb <= S_kb - Sbl
1472  if (sref(i,kmb+1) - s(i,kmb) < eps_dskb*(sref(i,kmb+2) - sref(i,kmb+1))) then
1473  dskb(i) = eps_dskb * (sref(i,kmb+2) - sref(i,kmb+1)) ; deriv_dskb(i) = 0.0
1474  else
1475  dskb(i) = sref(i,kmb+1) - s(i,kmb) ; deriv_dskb(i) = -1.0
1476  endif
1477  if (present(ddskb_de)) ddskb_de(i) = deriv_dskb(i)*ds_de(i,kmb)
1478  endif ; enddo
1479 
1480  if (present(dslay)) then
1481  dz_drat = 1000.0 ! The limit of large dz_drat the same as choosing a
1482  ! Heaviside function.
1483  eps_dslay = 1.0e-10 ! Should be ~= GV%Angstrom / sqrt(Kd*dt)
1484  do i=is,ie ; if (do_i(i)) then
1485  ds_kbp1 = sref(i,kmb+2) - sref(i,kmb+1)
1486  ids_kbp1 = 1.0 / (sref(i,kmb+2) - sref(i,kmb+1))
1487  rat = (sref(i,kmb+1) - s(i,kmb)) * ids_kbp1
1488  ! Need to ensure that 0 < dSLay <= 2*dSkb
1489  if (rat < 0.5) then
1490  ! The coefficients here are chosen so that at rat = 0.5, the value (1.5)
1491  ! and first derivative (-0.5) match with the "typical" case (next).
1492  ! The functional form here is arbitrary.
1493  ! f1 provides a reasonable profile that matches the value and derivative
1494  ! of the "typical" case at rat = 0.5, and has a maximum of less than 2.
1495  inv_term = 1.0 / (1.0-rat)
1496  f1 = 2.0 - 0.125*(inv_term**2)
1497  df1_drat = - 0.25*(inv_term**3)
1498 
1499  ! f2 ensures that dSLay goes to 0 rapidly if rat is significantly
1500  ! negative.
1501  z = dz_drat * rat + 4.0 ! The 4 here gives f2(0) = 0.982.
1502  if (z >= 18.0) then ; f2 = 1.0 ; df2_dz = 0.0
1503  elseif (z <= -58.0) then ; f2 = eps_dslay ; df2_dz = 0.0
1504  else
1505  expz = exp(z) ; inv_term = 1.0 / (1.0 + expz)
1506  f2 = (eps_dslay + expz) * inv_term
1507  df2_dz = (1.0 - eps_dslay) * expz * inv_term**2
1508  endif
1509 
1510  dslay(i) = dskb(i) * f1 * f2
1511  deriv_dslay = deriv_dskb(i) * (f1 * f2) - (dskb(i)*ids_kbp1) * &
1512  (df1_drat*f2 + f1 * dz_drat * df2_dz)
1513  elseif (dskb(i) <= 3.0*ds_kbp1) then
1514  ! This is the "typical" case.
1515  dslay(i) = 0.5 * (dskb(i) + ds_kbp1)
1516  deriv_dslay = 0.5 * deriv_dskb(i) ! = -0.5
1517  else
1518  dslay(i) = 2.0*ds_kbp1
1519  deriv_dslay = 0.0
1520  endif
1521  if (present(ddslay_de)) ddslay_de(i) = deriv_dslay*ds_de(i,kmb)
1522  endif ; enddo
1523  endif ! present(dSlay)
1524  endif ! Not limited.
1525 
1526  if (present(ds_anom_lim)) then ; do i=is,ie ; if (do_i(i)) then
1527  ds_anom_lim(i) = max(0.0, eps_dskb * (sref(i,kmb+2) - sref(i,kmb+1)) - &
1528  (sref(i,kmb+1) - s(i,kmb)) )
1529  endif ; enddo ; endif
1530 
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( g %isd: g %ied, g %ke), intent(in)  h_bl,
real, dimension( g %isd: g %ied), intent(in)  dtKd_kb,
real, dimension( g %isd: g %ied, g %ke), intent(in)  Sref,
real, dimension( g %isd: g %ied), intent(in)  I_dSkbp1,
real, dimension( g %isd: g %ied, g %ke), intent(in)  Ent_bl,
real, dimension( g %isd: g %ied), intent(in)  ea_kbp1,
real, dimension( g %isd: g %ied), intent(in)  min_eakb,
real, dimension( g %isd: g %ied), intent(in)  max_eakb,
integer, intent(in)  kmb,
integer, intent(in)  is,
integer, intent(in)  ie,
logical, dimension( g %isd: g %ied), 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(szi_(g)), intent(inout)  Ent,
real, dimension(szi_(g)), intent(out), optional  error,
real, dimension(szi_(g)), intent(in), optional  err_min_eakb0,
real, dimension(szi_(g)), intent(in), optional  err_max_eakb0,
real, dimension(szi_(g)), intent(out), optional  F_kb,
real, dimension(szi_(g)), intent(out), optional  dFdfm_kb 
)
private
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, in units of m or kg m-2 (abbreviated as H below).
[in]srefThe coordinate reference potential density, with the value of the topmost interior layer at layer kmb+1, in units of kg m-3.
[in]ent_blThe average entrainment upward and downward across each interface around the buffer layers, in H.
[in]i_dskbp1The inverse of the difference in reference potential density across the base of the uppermost interior layer, in units of m3 kg-1.
[in]dtkd_kbThe diapycnal diffusivity in the top interior layer times the time step, in H2.
[in]ea_kbp1The entrainment from above by layer kb+1, in H.
[in]min_eakbThe minimum permissible rate of entrainment, in H.
[in]max_eakbThe maximum permissible rate of entrainment, in H.
[in]ieThe range of i-indices 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, in H. The input value is the first guess.
[out]errorThe error (locally defined in this routine) associated with the returned solution.
[in]err_max_eakb0The errors (locally defined) associated with min_eakb and max_eakb when ea_kbp1 = 0, returned from a previous call to this routine.
[out]f_kbThe entrainment from below by the uppermost interior layer corresponding to the returned value of Ent, in H.
[out]dfdfm_kbThe partial derivative of F_kb with ea_kbp1, nondim.

Definition at line 1653 of file MOM_entrain_diffusive.F90.

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

Referenced by entrainment_diffusive().

1653  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
1654  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
1655  real, dimension(SZI_(G),SZK_(G)), intent(in) :: h_bl !< Layer thickness, with the top
1656  !! interior layer at k-index kmb+1, in
1657  !! units of m or kg m-2
1658  !! (abbreviated as H below).
1659  real, dimension(SZI_(G),SZK_(G)), intent(in) :: sref !< The coordinate reference potential
1660  !! density, with the value of the
1661  !! topmost interior layer at layer
1662  !! kmb+1, in units of kg m-3.
1663  real, dimension(SZI_(G),SZK_(G)), intent(in) :: ent_bl !< The average entrainment upward and
1664  !! downward across each interface around
1665  !! the buffer layers, in H.
1666  real, dimension(SZI_(G)), intent(in) :: i_dskbp1 !< The inverse of the difference in
1667  !! reference potential density across
1668  !! the base of the uppermost interior
1669  !! layer, in units of m3 kg-1.
1670  real, dimension(SZI_(G)), intent(in) :: dtkd_kb !< The diapycnal diffusivity in the top
1671  !! interior layer times the time step,
1672  !! in H2.
1673  real, dimension(SZI_(G)), intent(in) :: ea_kbp1 !< The entrainment from above by layer
1674  !! kb+1, in H.
1675  real, dimension(SZI_(G)), intent(in) :: min_eakb !< The minimum permissible rate of
1676  !! entrainment, in H.
1677  real, dimension(SZI_(G)), intent(in) :: max_eakb !< The maximum permissible rate of
1678  !! entrainment, in H.
1679  integer, intent(in) :: kmb
1680  integer, intent(in) :: is, ie !< The range of i-indices to work on.
1681  logical, dimension(SZI_(G)), intent(in) :: do_i !< A logical variable indicating which
1682  !! i-points to work on.
1683  type(entrain_diffusive_cs), pointer :: cs !< This module's control structure.
1684  real, dimension(SZI_(G)), intent(inout) :: ent !< The entrainment rate of the uppermost
1685  !! interior layer, in H. The input value
1686  !! is the first guess.
1687  real, dimension(SZI_(G)), intent(out), optional :: error !< The error (locally defined in this
1688  !! routine) associated with the returned
1689  !! solution.
1690  real, dimension(SZI_(G)), intent(in), optional :: err_min_eakb0, err_max_eakb0 !< The errors
1691  !! (locally defined) associated with
1692  !! min_eakb and max_eakb when ea_kbp1
1693  !! = 0, returned from a previous call
1694  !! to this routine.
1695  real, dimension(SZI_(G)), intent(out), optional :: f_kb !< The entrainment from below by the
1696  !! uppermost interior layer
1697  !! corresponding to the returned
1698  !! value of Ent, in H.
1699  real, dimension(SZI_(G)), intent(out), optional :: dfdfm_kb !< The partial derivative of F_kb with
1700  !! ea_kbp1, nondim.
1701 
1702 ! Arguments: h_bl - Layer thickness, with the top interior layer at k-index
1703 ! kmb+1, in units of m or kg m-2 (abbreviated as H below).
1704 ! (in) dtKd_kb - The diapycnal diffusivity in the top interior layer times
1705 ! the time step, in H2.
1706 ! (in) Sref - The coordinate reference potential density, with the
1707 ! value of the topmost interior layer at layer kmb+1,
1708 ! in units of kg m-3.
1709 ! (in) I_dSkbp1 - The inverse of the difference in reference potential
1710 ! density across the base of the uppermost interior layer,
1711 ! in units of m3 kg-1.
1712 ! (in) Ent_bl - The average entrainment upward and downward across
1713 ! each interface around the buffer layers, in H.
1714 ! (in) ea_kbp1 - The entrainment from above by layer kb+1, in H.
1715 ! (in) min_eakb - The minimum permissible rate of entrainment, in H.
1716 ! (in) max_eakb - The maximum permissible rate of entrainment, in H.
1717 ! (in) is, ie - The range of i-indices to work on.
1718 ! (in) do_i - A logical variable indicating which i-points to work on.
1719 ! (in) G - The ocean's grid structure.
1720 ! (in) GV - The ocean's vertical grid structure.
1721 ! (in) CS - This module's control structure.
1722 ! (in/out) Ent - The entrainment rate of the uppermost interior layer, in H.
1723 ! The input value is the first guess.
1724 ! (out,opt) error - The error (locally defined in this routine) associated with
1725 ! the returned solution.
1726 ! (in,opt) error_min_eakb0, error_max_eakb0 - The errors (locally defined)
1727 ! associated with min_eakb and max_eakb when ea_kbp1 = 0,
1728 ! returned from a previous call to this routine.
1729 ! (out,opt) F_kb - The entrainment from below by the uppermost interior layer
1730 ! corresponding to the returned value of Ent, in H.
1731 ! (out,out) dFdfm_kb - The partial derivative of F_kb with ea_kbp1, nondim.
1732 
1733 ! This subroutine determines the entrainment from above by the top interior
1734 ! layer (labeled kb elsewhere) given an entrainment by the layer below it,
1735 ! constrained to be within the provided bounds.
1736  real, dimension(SZI_(G)) :: &
1737  ds_kb, & ! The coordinate-density difference between the
1738  ! layer kb and deepest buffer layer, limited to
1739  ! ensure that it is positive, in kg m-3.
1740  ds_lay, & ! The coordinate-density difference across layer
1741  ! kb, limited to ensure that it is positive and not
1742  ! too much bigger than dS_kb or dS_kbp1, in kg m-3.
1743  ddskb_de, ddslay_de, & ! The derivatives of dS_kb and dS_Lay with E,
1744  ! in units of kg m-3 H-1.
1745  derror_de, & ! The derivative of err with E, in H.
1746  err, & ! The "error" whose zero is being sought, in H2.
1747  e_min, e_max, & ! The minimum and maximum values of E, in H.
1748  error_mine, error_maxe ! err when E = E_min or E = E_max, in H2.
1749  real :: err_est ! An estimate of what err will be, in H2.
1750  real :: el ! 1 or 0, depending on whether increases in E lead
1751  ! to decreases in the entrainment from below by the
1752  ! deepest buffer layer.
1753  real :: fa, fk, fm, fr ! Temporary variables used to calculate err, in ND, H2, H, H.
1754  real :: tolerance ! The tolerance within which E must be converged, in H.
1755  real :: e_prev ! The previous value of E, in H.
1756  logical, dimension(SZI_(G)) :: false_position ! If true, the false position
1757  ! method might be used for the next iteration.
1758  logical, dimension(SZI_(G)) :: redo_i ! If true, more work is needed on this column.
1759  logical :: do_any
1760  real, parameter :: large_val = 1.0e30
1761  integer :: i, it
1762  integer, parameter :: maxit = 30
1763 
1764  if (.not.cs%bulkmixedlayer) then
1765  call mom_error(fatal, "determine_Ea_kb should not be called "//&
1766  "unless BULKMIXEDLAYER is defined.")
1767  endif
1768  tolerance = gv%m_to_H * cs%Tolerance_Ent
1769 
1770  do i=is,ie ; redo_i(i) = do_i(i) ; enddo
1771 
1772  do i=is,ie ; if (do_i(i)) then
1773  ! The first guess of Ent was the value from the previous iteration.
1774 
1775  ! These were previously calculated and provide good limits and estimates
1776  ! of the errors there. By construction the errors increase with R*ea_kbp1.
1777  e_min(i) = min_eakb(i) ; e_max(i) = max_eakb(i)
1778  error_mine(i) = -large_val ; error_maxe(i) = large_val
1779  false_position(i) = .true. ! Used to alternate between false_position and
1780  ! bisection when Newton's method isn't working.
1781  if (present(err_min_eakb0)) error_mine(i) = err_min_eakb0(i) - e_min(i) * ea_kbp1(i)
1782  if (present(err_max_eakb0)) error_maxe(i) = err_max_eakb0(i) - e_max(i) * ea_kbp1(i)
1783 
1784  if ((error_maxe(i) <= 0.0) .or. (error_mine(i) >= 0.0)) then
1785  ! The root is not bracketed and one of the limiting values should be used.
1786  if (error_maxe(i) <= 0.0) then
1787  ! The errors decrease with E*ea_kbp1, so E_max is the best solution.
1788  ent(i) = e_max(i) ; err(i) = error_maxe(i)
1789  else ! error_minE >= 0 is equivalent to ea_kbp1 = 0.0.
1790  ent(i) = e_min(i) ; err(i) = error_mine(i)
1791  endif
1792  derror_de(i) = 0.0
1793  redo_i(i) = .false.
1794  endif
1795  endif ; enddo ! End of i-loop
1796 
1797  do it = 1,maxit
1798  do_any = .false. ; do i=is,ie ; if (redo_i(i)) do_any = .true. ; enddo
1799  if (.not.do_any) exit
1800  call determine_dskb(h_bl, sref, ent_bl, ent, is, ie, kmb, g, gv, .true., ds_kb, &
1801  ddskb_de, ds_lay, ddslay_de, do_i_in = redo_i)
1802  do i=is,ie ; if (redo_i(i)) then
1803  ! The correct root is bracketed between E_min and E_max.
1804  ! Note the following limits: Ent >= 0 ; fa > 1 ; fk > 0
1805  el = 0.0 ; if (2.0*ent_bl(i,kmb+1) >= ent(i)) el = 1.0
1806  fa = (1.0 + el) + ds_kb(i)*i_dskbp1(i)
1807  fk = dtkd_kb(i) * (ds_lay(i)/ds_kb(i))
1808  fm = (ea_kbp1(i) - h_bl(i,kmb+1)) + el*2.0*ent_bl(i,kmb+1)
1809  if (fm > -gv%Angstrom) fm = fm + gv%Angstrom ! This could be smooth if need be.
1810  err(i) = (fa * ent(i)**2 - fm * ent(i)) - fk
1811  derror_de(i) = ((2.0*fa + (ddskb_de(i)*i_dskbp1(i))*ent(i))*ent(i) - fm) - &
1812  dtkd_kb(i) * (ddslay_de(i)*ds_kb(i) - ddskb_de(i)*ds_lay(i))/(ds_kb(i)**2)
1813 
1814  if (err(i) == 0.0) then
1815  redo_i(i) = .false. ; cycle
1816  elseif (err(i) > 0.0) then
1817  e_max(i) = ent(i) ; error_maxe(i) = err(i)
1818  else
1819  e_min(i) = ent(i) ; error_mine(i) = err(i)
1820  endif
1821 
1822  e_prev = ent(i)
1823  if ((it == 1) .or. (derror_de(i) <= 0.0)) then
1824  ! Assuming that the coefficients of the quadratic equation are correct
1825  ! will usually give a very good first guess. Also, if derror_dE < 0.0,
1826  ! R is on the wrong side of the approximate parabola. In either case,
1827  ! try assuming that the error is approximately a parabola and solve.
1828  fr = sqrt(fm**2 + 4.0*fa*fk)
1829  if (fm >= 0.0) then
1830  ent(i) = (fm + fr) / (2.0 * fa)
1831  else
1832  ent(i) = (2.0 * fk) / (fr - fm)
1833  endif
1834  ! But make sure that the root stays bracketed, bisecting if needed.
1835  if ((ent(i) > e_max(i)) .or. (ent(i) < e_min(i))) &
1836  ent(i) = 0.5*(e_max(i) + e_min(i))
1837  elseif (((e_max(i)-ent(i))*derror_de(i) > -err(i)) .and. &
1838  ((ent(i)-e_min(i))*derror_de(i) > err(i)) ) then
1839  ! Use Newton's method for the next estimate, provided it will
1840  ! remain bracketed between Rmin and Rmax.
1841  ent(i) = ent(i) - err(i) / derror_de(i)
1842  elseif (false_position(i) .and. &
1843  (error_maxe(i) - error_mine(i) < 0.9*large_val)) then
1844  ! Use the false postion method if there are decent error estimates.
1845  ent(i) = e_min(i) + (e_max(i)-e_min(i)) * &
1846  (-error_mine(i)/(error_maxe(i) - error_mine(i)))
1847  false_position(i) = .false.
1848  else ! Bisect as a last resort or if the false position method was used last.
1849  ent(i) = 0.5*(e_max(i) + e_min(i))
1850  false_position(i) = .true.
1851  endif
1852 
1853  if (abs(e_prev - ent(i)) < tolerance) then
1854  err_est = err(i) + (ent(i) - e_prev) * derror_de(i)
1855  if ((it > 1) .or. (err_est*err(i) <= 0.0) .or. &
1856  (abs(err_est) < abs(tolerance*derror_de(i)))) redo_i(i) = .false.
1857  endif
1858 
1859  endif ; enddo ! End of i-loop
1860  enddo ! End of iterations to determine Ent(i).
1861 
1862  ! Update the value of dS_kb for consistency with Ent.
1863  if (present(f_kb) .or. present(dfdfm_kb)) &
1864  call determine_dskb(h_bl, sref, ent_bl, ent, is, ie, kmb, g, gv, .true., &
1865  ds_kb, do_i_in = do_i)
1866 
1867  if (present(f_kb)) then ; do i=is,ie ; if (do_i(i)) then
1868  f_kb(i) = ent(i) * (ds_kb(i) * i_dskbp1(i))
1869  endif ; enddo ; endif
1870  if (present(error)) then ; do i=is,ie ; if (do_i(i)) then
1871  error(i) = err(i)
1872  endif ; enddo ; endif
1873  if (present(dfdfm_kb)) then ; do i=is,ie ; if (do_i(i)) then
1874  ! derror_dE and ddSkb_dE are _not_ recalculated here, since dFdfm_kb is
1875  ! only used in Newton's method, and slightly increasing the accuracy of the
1876  ! estimate is unlikely to speed convergence.
1877  if (derror_de(i) > 0.0) then
1878  dfdfm_kb(i) = ((ds_kb(i) + ent(i) * ddskb_de(i)) * i_dskbp1(i)) * &
1879  (ent(i) / derror_de(i))
1880  else ! Use Adcroft's division by 0 convention.
1881  dfdfm_kb(i) = 0.0
1882  endif
1883  endif ; enddo ; endif
1884 
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)

Definition at line 2268 of file MOM_entrain_diffusive.F90.

Referenced by mom_diabatic_driver::diabatic_driver_end().

2268  type(entrain_diffusive_cs), pointer :: cs
2269 
2270  if (associated(cs)) deallocate(cs)
2271 
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(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(entrain_diffusive_cs), pointer  CS 
)
Parameters
[in]timeThe current model time.
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[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 2206 of file MOM_entrain_diffusive.F90.

References mom_error_handler::mom_error().

2206  type(time_type), intent(in) :: time !< The current model time.
2207  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
2208  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
2209  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
2210  !! parameters.
2211  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
2212  !! output.
2213  type(entrain_diffusive_cs), pointer :: cs !< A pointer that is set to point to the control
2214  !! structure.
2215 ! for this module
2216 ! Arguments: Time - The current model time.
2217 ! (in) G - The ocean's grid structure.
2218 ! (in) GV - The ocean's vertical grid structure.
2219 ! (in) param_file - A structure indicating the open file to parse for
2220 ! model parameter values.
2221 ! (in) diag - A structure that is used to regulate diagnostic output.
2222 ! (in/out) CS - A pointer that is set to point to the control structure
2223 ! for this module
2224  real :: decay_length, dt, kd
2225 ! This include declares and sets the variable "version".
2226 #include "version_variable.h"
2227  character(len=40) :: mod = "MOM_entrain_diffusive" ! This module's name.
2228 
2229  if (associated(cs)) then
2230  call mom_error(warning, "entrain_diffusive_init called with an associated "// &
2231  "control structure.")
2232  return
2233  endif
2234  allocate(cs)
2235 
2236  cs%diag => diag
2237 
2238  cs%bulkmixedlayer = (gv%nkml > 0)
2239 
2240 ! Set default, read and log parameters
2241  call log_version(param_file, mod, version, "")
2242  call get_param(param_file, mod, "CORRECT_DENSITY", cs%correct_density, &
2243  "If true, and USE_EOS is true, the layer densities are \n"//&
2244  "restored toward their target values by the diapycnal \n"//&
2245  "mixing, as described in Hallberg (MWR, 2000).", &
2246  default=.true.)
2247  call get_param(param_file, mod, "MAX_ENT_IT", cs%max_ent_it, &
2248  "The maximum number of iterations that may be used to \n"//&
2249  "calculate the interior diapycnal entrainment.", default=5)
2250 ! In this module, KD is only used to set the default for TOLERANCE_ENT. (m2 s-1)
2251  call get_param(param_file, mod, "KD", kd, fail_if_missing=.true.)
2252  call get_param(param_file, mod, "DT", dt, &
2253  "The (baroclinic) dynamics time step.", units = "s", &
2254  fail_if_missing=.true.)
2255 ! CS%Tolerance_Ent = MAX(100.0*GV%Angstrom,1.0e-4*sqrt(dt*Kd)) !
2256  call get_param(param_file, mod, "TOLERANCE_ENT", cs%Tolerance_Ent, &
2257  "The tolerance with which to solve for entrainment values.", &
2258  units="m", default=max(100.0*gv%Angstrom,1.0e-4*sqrt(dt*kd)))
2259 
2260  cs%id_Kd = register_diag_field('ocean_model', 'Kd_effective', diag%axesTL, time, &
2261  'Diapycnal diffusivity as applied', 'meter2 second-1')
2262  cs%id_diff_work = register_diag_field('ocean_model', 'diff_work', diag%axesTi, time, &
2263  'Work actually done by diapycnal diffusion across each interface', 'W m-2')
2264 
Here is the call graph for this function:

◆ entrainment_diffusive()

subroutine, public mom_entrain_diffusive::entrainment_diffusive ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(in)  v,
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(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]uThe zonal velocity, in m s-1.
[in]vThe meridional velocity, in m s-1.
[in]hLayer thicknesses, in H (usually 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 in 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 111 of file MOM_entrain_diffusive.F90.

References mom_eos::calculate_density_derivs(), 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().

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

Definition at line 1536 of file MOM_entrain_diffusive.F90.

References determine_dskb(), and find_maxf_kb().

Referenced by entrainment_diffusive().

1536  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
1537  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
1538  real, dimension(SZI_(G),SZK_(G)), intent(in) :: h_bl, sref, ent_bl
1539  real, dimension(SZI_(G)), intent(in) :: i_dskbp1, f_kb
1540  integer, intent(in) :: kmb, i
1541  type(entrain_diffusive_cs), pointer :: cs
1542  real, dimension(SZI_(G)), intent(inout) :: ea_kb
1543  real, optional, intent(in) :: tol_in
1544 
1545  ! Given an entrainment from below for layer kb, determine a consistent
1546  ! entrainment from above, such that dSkb * ea_kb = dSkbp1 * F_kb. The input
1547  ! value of ea_kb is both the maximum value that can be obtained and the first
1548  ! guess of the iterations. Also, make sure that ea_kb is an under-estimate
1549  real :: max_ea, min_ea
1550  real :: err, err_min, err_max
1551  real :: derr_dea
1552  real :: val, tolerance, tol1
1553  real :: ea_prev
1554  real :: ds_kbp1
1555  logical :: bisect_next, newton
1556  real, dimension(SZI_(G)) :: ds_kb
1557  real, dimension(SZI_(G)) :: maxf, ent_maxf, zeros
1558  real, dimension(SZI_(G)) :: ddskb_de
1559  integer :: it
1560  integer, parameter :: maxit = 30
1561 
1562  ds_kbp1 = sref(i,kmb+2) - sref(i,kmb+1)
1563  max_ea = ea_kb(i) ; min_ea = 0.0
1564  val = ds_kbp1 * f_kb(i)
1565  err_min = -val
1566 
1567  tolerance = gv%m_to_H * cs%Tolerance_Ent
1568  if (present(tol_in)) tolerance = tol_in
1569  bisect_next = .true.
1570 
1571  call determine_dskb(h_bl, sref, ent_bl, ea_kb, i, i, kmb, g, gv, .true., &
1572  ds_kb, ddskb_de)
1573 
1574  err = ds_kb(i) * ea_kb(i) - val
1575  derr_dea = ds_kb(i) + ddskb_de(i) * ea_kb(i)
1576  ! Return if Newton's method on the first guess would give a tolerably small
1577  ! change in the value of ea_kb.
1578  if ((err <= 0.0) .and. (abs(err) <= tolerance*abs(derr_dea))) return
1579 
1580  if (err == 0.0) then ; return ! The exact solution on the first guess...
1581  elseif (err > 0.0) then ! The root is properly bracketed.
1582  max_ea = ea_kb(i) ; err_max = err
1583  ! Use Newton's method (if it stays bounded) or the false position method
1584  ! to find the next value.
1585  if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i) - min_ea) > err) .and. &
1586  (derr_dea*(max_ea - ea_kb(i)) > -1.0*err)) then
1587  ea_kb(i) = ea_kb(i) - err / derr_dea
1588  else ! Use the bisection for the next guess.
1589  ea_kb(i) = 0.5*(max_ea+min_ea)
1590  endif
1591  else
1592  ! Try to bracket the root first. If unable to bracket the root, return
1593  ! the maximum.
1594  zeros(i) = 0.0
1595  call find_maxf_kb(h_bl, sref, ent_bl, i_dskbp1, zeros, ea_kb, &
1596  kmb, i, i, g, gv, cs, maxf, ent_maxf, f_thresh = f_kb)
1597  err_max = ds_kbp1 * maxf(i) - val
1598  ! If err_max is negative, there is no good solution, so use the maximum
1599  ! value of F in the valid range.
1600  if (err_max <= 0.0) then
1601  ea_kb(i) = ent_maxf(i) ; return
1602  else
1603  max_ea = ent_maxf(i)
1604  ea_kb(i) = 0.5*(max_ea+min_ea) ! Use bisection for the next guess.
1605  endif
1606  endif
1607 
1608  ! Exit if the range between max_ea and min_ea already acceptable.
1609  ! if (abs(max_ea - min_ea) < 0.1*tolerance) return
1610 
1611  do it = 1, maxit
1612  call determine_dskb(h_bl, sref, ent_bl, ea_kb, i, i, kmb, g, gv, .true., &
1613  ds_kb, ddskb_de)
1614 
1615  err = ds_kb(i) * ea_kb(i) - val
1616  derr_dea = ds_kb(i) + ddskb_de(i) * ea_kb(i)
1617 
1618  ea_prev = ea_kb(i)
1619  ! Use Newton's method or the false position method to find the next value.
1620  newton = .false.
1621  if (err > 0.0) then
1622  max_ea = ea_kb(i) ; err_max = err
1623  if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i)-min_ea) > err)) newton = .true.
1624  else
1625  min_ea = ea_kb(i) ; err_min = err
1626  if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i)-max_ea) < err)) newton = .true.
1627  endif
1628 
1629  if (newton) then
1630  ea_kb(i) = ea_kb(i) - err / derr_dea
1631  elseif (bisect_next) then ! Use bisection to reduce the range.
1632  ea_kb(i) = 0.5*(max_ea+min_ea)
1633  bisect_next = .false.
1634  else ! Use the false-position method for the next guess.
1635  ea_kb(i) = min_ea + (max_ea-min_ea) * (err_min/(err_min - err_max))
1636  bisect_next = .true.
1637  endif
1638 
1639  tol1 = tolerance ; if (err > 0.0) tol1 = 0.099*tolerance
1640  if (ds_kb(i) <= ds_kbp1) then
1641  if (abs(ea_kb(i) - ea_prev) <= tol1) return
1642  else
1643  if (ds_kbp1*abs(ea_kb(i) - ea_prev) <= ds_kb(i)*tol1) return
1644  endif
1645  enddo
1646 
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(szi_(g),szk_(g)), intent(in)  dsp1_ds,
real, dimension(szi_(g)), intent(in)  eakb,
real, dimension(szi_(g),szk_(g)), intent(in)  Ent_bl,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  ea,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  eb,
logical, dimension(szi_(g)), 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]hLayer thicknesses, in H (usually m or kg m-2)

Definition at line 979 of file MOM_entrain_diffusive.F90.

Referenced by entrainment_diffusive().

979  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
980  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
981  real, dimension(SZI_(G),SZK_(G)), intent(in) :: f
982  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2)
983  integer, dimension(SZI_(G)), intent(in) :: kb
984  integer, intent(in) :: kmb, j
985  type(entrain_diffusive_cs), intent(in) :: cs
986  real, dimension(SZI_(G),SZK_(G)), intent(in) :: dsp1_ds
987  real, dimension(SZI_(G)), intent(in) :: eakb
988  real, dimension(SZI_(G),SZK_(G)), intent(in) :: ent_bl
989  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: ea, eb
990  logical, dimension(SZI_(G)), optional, intent(in) :: do_i_in
991 ! This subroutine calculates the actual entrainments (ea and eb) and the
992 ! amount of surface forcing that is applied to each layer if there is no bulk
993 ! mixed layer.
994 
995  real :: h1 ! The thickness in excess of the minimum that will remain
996  ! after exchange with the layer below, in m or kg m-2.
997  logical :: do_i(szi_(g))
998  integer :: i, k, is, ie, nz
999 
1000  is = g%isc ; ie = g%iec ; nz = g%ke
1001 
1002  if (present(do_i_in)) then
1003  do i=is,ie ; do_i(i) = do_i_in(i) ; enddo
1004  do i=g%isc,g%iec ; if (do_i(i)) then
1005  is = i ; exit
1006  endif ; enddo
1007  do i=g%iec,g%isc,-1 ; if (do_i(i)) then
1008  ie = i ; exit
1009  endif ; enddo
1010  else
1011  do i=is,ie ; do_i(i) = .true. ; enddo
1012  endif
1013 
1014  do i=is,ie
1015  ea(i,j,nz) = 0.0 ; eb(i,j,nz) = 0.0
1016  enddo
1017  if (cs%bulkmixedlayer) then
1018  do i=is,ie
1019  eb(i,j,kmb) = max(2.0*ent_bl(i,kmb+1) - eakb(i), 0.0)
1020  enddo
1021  do k=nz-1,kmb+1,-1 ; do i=is,ie ; if (do_i(i)) then
1022  if (k > kb(i)) then
1023  ! With a bulk mixed layer, surface buoyancy fluxes are applied
1024  ! elsewhere, so F should always be nonnegative.
1025  ea(i,j,k) = dsp1_ds(i,k)*f(i,k)
1026  eb(i,j,k) = f(i,k)
1027  else if (k == kb(i)) then
1028  ea(i,j,k) = eakb(i)
1029  eb(i,j,k) = f(i,k)
1030  elseif (k == kb(i)-1) then
1031  ea(i,j,k) = ea(i,j,k+1)
1032  eb(i,j,k) = eb(i,j,kmb)
1033  else
1034  ea(i,j,k) = ea(i,j,k+1)
1035  ! Add the entrainment of the thin interior layers to eb going
1036  ! up into the buffer layer.
1037  eb(i,j,k) = eb(i,j,k+1) + max(0.0, h(i,j,k+1) - gv%Angstrom)
1038  endif
1039  endif ; enddo ; enddo
1040  k = kmb
1041  do i=is,ie ; if (do_i(i)) then
1042  ! Adjust the previously calculated entrainment from below by the deepest
1043  ! buffer layer to account for entrainment of thin interior layers .
1044  if (kb(i) > kmb+1) &
1045  eb(i,j,k) = eb(i,j,k+1) + max(0.0, h(i,j,k+1) - gv%Angstrom)
1046 
1047  ! Determine the entrainment from above for each buffer layer.
1048  h1 = (h(i,j,k) - gv%Angstrom) + (eb(i,j,k) - ea(i,j,k+1))
1049  ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)
1050  endif ; enddo
1051  do k=kmb-1,2,-1 ; do i=is,ie ; if (do_i(i)) then
1052  ! Determine the entrainment from below for each buffer layer.
1053  eb(i,j,k) = max(2.0*ent_bl(i,k+1) - ea(i,j,k+1), 0.0)
1054 
1055  ! Determine the entrainment from above for each buffer layer.
1056  h1 = (h(i,j,k) - gv%Angstrom) + (eb(i,j,k) - ea(i,j,k+1))
1057  ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)
1058 ! if (h1 >= 0.0) then ; ea(i,j,k) = Ent_bl(i,K)
1059 ! elseif (Ent_bl(i,K)+0.5*h1 >= 0.0) then ; ea(i,j,k) = Ent_bl(i,K)-0.5*h1
1060 ! else ; ea(i,j,k) = -h1 ; endif
1061  endif ; enddo ; enddo
1062  do i=is,ie ; if (do_i(i)) then
1063  eb(i,j,1) = max(2.0*ent_bl(i,2) - ea(i,j,2), 0.0)
1064  ea(i,j,1) = 0.0
1065  endif ; enddo
1066  else ! not BULKMIXEDLAYER
1067  ! Calculate the entrainment by each layer from above and below.
1068  ! Entrainment is always positive, but F may be negative due to
1069  ! surface buoyancy fluxes.
1070  do i=is,ie
1071  ea(i,j,1) = 0.0 ; eb(i,j,1) = max(f(i,1),0.0)
1072  ea(i,j,2) = dsp1_ds(i,2)*f(i,2) - min(f(i,1),0.0)
1073  enddo
1074 
1075  do k=2,nz-1 ; do i=is,ie
1076  eb(i,j,k) = max(f(i,k),0.0)
1077  ea(i,j,k+1) = dsp1_ds(i,k+1)*f(i,k+1) - (f(i,k)-eb(i,j,k))
1078  if (ea(i,j,k+1) < 0.0) then
1079  eb(i,j,k) = eb(i,j,k) - ea(i,j,k+1)
1080  ea(i,j,k+1) = 0.0
1081  endif
1082  enddo ; enddo
1083  endif ! end BULKMIXEDLAYER
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
Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]h_blLayer thickness, in m or kg m-2
[in]srefReference potential density (in 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, in units of m3 kg-1.
[in]min_ent_inThe minimum value of ent to search, in H.
[in]max_ent_inThe maximum value of ent to search, in H.
[in]ieThe range of i-indices 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, in H.
[out]ent_maxfThe value of ent at that maximum, in H.
[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 1890 of file MOM_entrain_diffusive.F90.

References determine_dskb().

Referenced by entrainment_diffusive(), and f_kb_to_ea_kb().

1890  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
1891  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
1892  real, dimension(SZI_(G),SZK_(G)), &
1893  intent(in) :: h_bl !< Layer thickness, in m or kg m-2
1894  !! (abbreviated as H below).
1895  real, dimension(SZI_(G),SZK_(G)), &
1896  intent(in) :: sref !< Reference potential density (in kg m-3?).
1897  real, dimension(SZI_(G),SZK_(G)), &
1898  intent(in) :: ent_bl !< The average entrainment upward and
1899  !! downward across each interface around
1900  !! the buffer layers, in H.
1901  real, dimension(SZI_(G)), intent(in) :: i_dskbp1 !< The inverse of the difference in
1902  !! reference potential density across the
1903  !! base of the uppermost interior layer,
1904  !! in units of m3 kg-1.
1905  real, dimension(SZI_(G)), intent(in) :: min_ent_in !< The minimum value of ent to search,
1906  !! in H.
1907  real, dimension(SZI_(G)), intent(in) :: max_ent_in !< The maximum value of ent to search,
1908  !! in H.
1909  integer, intent(in) :: kmb
1910  integer, intent(in) :: is, ie !< The range of i-indices to work on.
1911  type(entrain_diffusive_cs), pointer :: cs !< This module's control structure.
1912  real, dimension(SZI_(G)), intent(out) :: maxf !< The maximum value of F
1913  !! = ent*ds_kb*I_dSkbp1 found in the range
1914  !! min_ent < ent < max_ent, in H.
1915  real, dimension(SZI_(G)), intent(out), &
1916  optional :: ent_maxf !< The value of ent at that maximum, in H.
1917  logical, dimension(SZI_(G)), intent(in), &
1918  optional :: do_i_in !< A logical array indicating which columns
1919  !! to work on.
1920  real, dimension(SZI_(G)), intent(out), &
1921  optional :: f_lim_maxent !< If present, do not apply the limit in
1922  !! finding the maximum value, but return the
1923  !! limited value at ent=max_ent_in in this
1924  !! array, in H.
1925  real, dimension(SZI_(G)), intent(in), &
1926  optional :: f_thresh !< If F_thresh is present, return the first
1927  !! value found that has F > F_thresh, or
1928  !! the maximum.
1929 
1930 ! Arguments: h_bl - Layer thickness, in m or kg m-2 (abbreviated as H below).
1931 ! (in) Sref - Reference potential density (in kg m-3?)
1932 ! (in) Ent_bl - The average entrainment upward and downward across
1933 ! each interface around the buffer layers, in H.
1934 ! (in) I_dSkbp1 - The inverse of the difference in reference potential
1935 ! density across the base of the uppermost interior layer,
1936 ! in units of m3 kg-1.
1937 ! (in) min_ent_in - The minimum value of ent to search, in H.
1938 ! (in) max_ent_in - The maximum value of ent to search, in H.
1939 ! (in) is, ie - The range of i-indices to work on.
1940 ! (in) G - The ocean's grid structure.
1941 ! (in) GV - The ocean's vertical grid structure.
1942 ! (in) CS - This module's control structure.
1943 ! (out) maxF - The maximum value of F = ent*ds_kb*I_dSkbp1 found in the
1944 ! range min_ent < ent < max_ent, in H.
1945 ! (out,opt) ent_maxF - The value of ent at that maximum, in H.
1946 ! (in, opt) do_i_in - A logical array indicating which columns to work on.
1947 ! (out,opt) F_lim_maxent - If present, do not apply the limit in finding the
1948 ! maximum value, but return the limited value at
1949 ! ent=max_ent_in in this array, in H.
1950 ! (in, opt) F_thresh - If F_thresh is present, return the first value found
1951 ! that has F > F_thresh, or the maximum.
1952 
1953 ! Maximize F = ent*ds_kb*I_dSkbp1 in the range min_ent < ent < max_ent.
1954 ! ds_kb may itself be limited to positive values in determine_dSkb, which gives
1955 ! the prospect of two local maxima in the range - one at max_ent_in with that
1956 ! minimum value of ds_kb, and the other due to the unlimited (potentially
1957 ! negative) value. It is faster to find the true maximum by first finding the
1958 ! unlimited maximum and comparing it to the limited value at max_ent_in.
1959  real, dimension(SZI_(G)) :: &
1960  ent, &
1961  minent, maxent, ent_best, &
1962  f_max_ent_in, &
1963  f_maxent, f_minent, f, f_best, &
1964  df_dent, df_de_max, df_de_min, df_de_best, &
1965  ds_kb, ds_kb_lim, ddskb_de, ds_anom_lim, &
1966  chg_prev, chg_pre_prev
1967  real :: df_de_mean, maxslope, minslope
1968  real :: tolerance
1969  real :: ratio_select_end
1970  real :: rat, max_chg, min_chg, chg1, chg2, chg
1971  logical, dimension(SZI_(G)) :: do_i, last_it, need_bracket, may_use_best
1972  logical :: doany, ok1, ok2, bisect, new_min_bound
1973  integer :: i, it, is1, ie1
1974  integer, parameter :: maxit = 20
1975 
1976  tolerance = gv%m_to_H * cs%Tolerance_Ent
1977 
1978  if (present(do_i_in)) then
1979  do i=is,ie ; do_i(i) = do_i_in(i) ; enddo
1980  else
1981  do i=is,ie ; do_i(i) = .true. ; enddo
1982  endif
1983 
1984  ! The most likely value is at max_ent.
1985  call determine_dskb(h_bl, sref, ent_bl, max_ent_in, is, ie, kmb, g, gv, .false., &
1986  ds_kb, ddskb_de , ds_anom_lim=ds_anom_lim)
1987  ie1 = is-1 ; doany = .false.
1988  do i=is,ie
1989  ds_kb_lim(i) = ds_kb(i) + ds_anom_lim(i)
1990  f_max_ent_in(i) = max_ent_in(i)*ds_kb_lim(i)*i_dskbp1(i)
1991  maxent(i) = max_ent_in(i) ; minent(i) = min_ent_in(i)
1992  if ((abs(maxent(i) - minent(i)) < tolerance) .or. (.not.do_i(i))) then
1993  f_best(i) = max_ent_in(i)*ds_kb(i)*i_dskbp1(i)
1994  ent_best(i) = max_ent_in(i) ; ent(i) = max_ent_in(i)
1995  do_i(i) = .false.
1996  else
1997  f_maxent(i) = maxent(i) * ds_kb(i) * i_dskbp1(i)
1998  df_de_max(i) = (ds_kb(i) + maxent(i)*ddskb_de(i)) * i_dskbp1(i)
1999  doany = .true. ; last_it(i) = .false. ; need_bracket(i) = .true.
2000  endif
2001  enddo
2002 
2003  if (doany) then
2004  ie1 = is-1 ; do i=is,ie ; if (do_i(i)) ie1 = i ; enddo
2005  do i=ie1,is,-1 ; if (do_i(i)) is1 = i ; enddo
2006  ! Find the value of F and its derivative at min_ent.
2007  call determine_dskb(h_bl, sref, ent_bl, minent, is1, ie1, kmb, g, gv, .false., &
2008  ds_kb, ddskb_de, do_i_in = do_i)
2009  do i=is1,ie1 ; if (do_i(i)) then
2010  f_minent(i) = minent(i) * ds_kb(i) * i_dskbp1(i)
2011  df_de_min(i) = (ds_kb(i) + minent(i)*ddskb_de(i)) * i_dskbp1(i)
2012  endif ; enddo
2013 
2014  ratio_select_end = 0.9
2015  do it=1,maxit
2016  ratio_select_end = 0.5*ratio_select_end
2017  do i=is1,ie1 ; if (do_i(i)) then
2018  if (need_bracket(i)) then
2019  df_de_mean = (f_maxent(i) - f_minent(i)) / (maxent(i) - minent(i))
2020  maxslope = max(df_de_mean, df_de_min(i), df_de_max(i))
2021  minslope = min(df_de_mean, df_de_min(i), df_de_max(i))
2022  if (f_minent(i) >= f_maxent(i)) then
2023  if (df_de_min(i) > 0.0) then ; rat = 0.02 ! A small step should bracket the soln.
2024  elseif (maxslope < ratio_select_end*minslope) then
2025  ! The maximum of F is at minent.
2026  f_best(i) = f_minent(i) ; ent_best(i) = minent(i) ; rat = 0.0
2027  do_i(i) = .false.
2028  else ; rat = 0.382 ; endif ! Use the golden ratio
2029  else
2030  if (df_de_max(i) < 0.0) then ; rat = 0.98 ! A small step should bracket the soln.
2031  elseif (minslope > ratio_select_end*maxslope) then
2032  ! The maximum of F is at maxent.
2033  f_best(i) = f_maxent(i) ; ent_best(i) = maxent(i) ; rat = 1.0
2034  do_i(i) = .false.
2035  else ; rat = 0.618 ; endif ! Use the golden ratio
2036  endif
2037 
2038  if (rat >= 0.0) ent(i) = rat*maxent(i) + (1.0-rat)*minent(i)
2039  if (((maxent(i) - minent(i)) < tolerance) .or. (it==maxit)) &
2040  last_it(i) = .true.
2041  else ! The maximum is bracketed by minent, ent_best, and maxent.
2042  chg1 = 2.0*(maxent(i) - minent(i)) ; chg2 = chg1
2043  if (df_de_best(i) > 0) then
2044  max_chg = maxent(i) - ent_best(i) ; min_chg = 0.0
2045  else
2046  max_chg = 0.0 ; min_chg = minent(i) - ent_best(i) ! < 0
2047  endif
2048  if (max_chg - min_chg < 2.0*tolerance) last_it(i) = .true.
2049  if (df_de_max(i) /= df_de_best(i)) &
2050  chg1 = (maxent(i) - ent_best(i))*df_de_best(i) / &
2051  (df_de_best(i) - df_de_max(i))
2052  if (df_de_min(i) /= df_de_best(i)) &
2053  chg2 = (minent(i) - ent_best(i))*df_de_best(i) / &
2054  (df_de_best(i) - df_de_min(i))
2055  ok1 = ((chg1 < max_chg) .and. (chg1 > min_chg))
2056  ok2 = ((chg2 < max_chg) .and. (chg2 > min_chg))
2057  if (.not.(ok1 .or. ok2)) then ; bisect = .true. ; else
2058  if (ok1 .and. ok2) then ! Take the acceptable smaller change.
2059  chg = chg1 ; if (abs(chg2) < abs(chg1)) chg = chg2
2060  elseif (ok1) then ; chg = chg1
2061  else ; chg = chg2 ; endif
2062  if (abs(chg) > 0.5*abs(chg_pre_prev(i))) then ; bisect = .true.
2063  else ; bisect = .false. ; endif
2064  endif
2065  chg_pre_prev(i) = chg_prev(i)
2066  if (bisect) then
2067  if (df_de_best(i) > 0.0) then
2068  ent(i) = 0.5*(maxent(i) + ent_best(i))
2069  chg_prev(i) = 0.5*(maxent(i) - ent_best(i))
2070  else
2071  ent(i) = 0.5*(minent(i) + ent_best(i))
2072  chg_prev(i) = 0.5*(minent(i) - ent_best(i))
2073  endif
2074  else
2075  if (abs(chg) < tolerance) chg = sign(tolerance,chg)
2076  ent(i) = ent_best(i) + chg
2077  chg_prev(i) = chg
2078  endif
2079  endif
2080  endif ; enddo
2081 
2082  if (mod(it,3) == 0) then ! Re-determine the loop bounds.
2083  ie1 = is-1 ; do i=is1,ie ; if (do_i(i)) ie1 = i ; enddo
2084  do i=ie1,is,-1 ; if (do_i(i)) is1 = i ; enddo
2085  endif
2086 
2087  call determine_dskb(h_bl, sref, ent_bl, ent, is1, ie1, kmb, g, gv, .false., &
2088  ds_kb, ddskb_de, do_i_in = do_i)
2089  do i=is1,ie1 ; if (do_i(i)) then
2090  f(i) = ent(i)*ds_kb(i)*i_dskbp1(i)
2091  df_dent(i) = (ds_kb(i) + ent(i)*ddskb_de(i)) * i_dskbp1(i)
2092  endif ; enddo
2093 
2094  if (present(f_thresh)) then ; do i=is1,ie1 ; if (do_i(i)) then
2095  if (f(i) >= f_thresh(i)) then
2096  f_best(i) = f(i) ; ent_best(i) = ent(i) ; do_i(i) = .false.
2097  endif
2098  endif ; enddo ; endif
2099 
2100  doany = .false.
2101  do i=is1,ie1 ; if (do_i(i)) then
2102  if (.not.last_it(i)) doany = .true.
2103  if (last_it(i)) then
2104  if (need_bracket(i)) then
2105  if ((f(i) > f_maxent(i)) .and. (f(i) > f_minent(i))) then
2106  f_best(i) = f(i) ; ent_best(i) = ent(i)
2107  elseif (f_maxent(i) > f_minent(i)) then
2108  f_best(i) = f_maxent(i) ; ent_best(i) = maxent(i)
2109  else
2110  f_best(i) = f_minent(i) ; ent_best(i) = minent(i)
2111  endif
2112  elseif (f(i) > f_best(i)) then
2113  f_best(i) = f(i) ; ent_best(i) = ent(i)
2114  endif
2115  do_i(i) = .false.
2116  elseif (need_bracket(i)) then
2117  if ((f(i) > f_maxent(i)) .and. (f(i) > f_minent(i))) then
2118  need_bracket(i) = .false. ! The maximum is now bracketed.
2119  chg_prev(i) = (maxent(i) - minent(i))
2120  chg_pre_prev(i) = 2.0*chg_prev(i)
2121  ent_best(i) = ent(i) ; f_best(i) = f(i) ; df_de_best(i) = df_dent(i)
2122  elseif ((f(i) <= f_maxent(i)) .and. (f(i) > f_minent(i))) then
2123  new_min_bound = .true. ! We have a new minimum bound.
2124  elseif ((f(i) <= f_maxent(i)) .and. (f(i) > f_minent(i))) then
2125  new_min_bound = .false. ! We have a new maximum bound.
2126  else ! This case would bracket a minimum. Wierd.
2127  ! Unless the derivative indicates that there is a maximum near the
2128  ! lower bound, try keeping the end with the larger value of F;
2129  ! in a tie keep the minimum as the answer here will be compared
2130  ! with the maximum input value later.
2131  new_min_bound = .true.
2132  if (df_de_min(i) > 0.0 .or. (f_minent(i) >= f_maxent(i))) &
2133  new_min_bound = .false.
2134  endif
2135  if (need_bracket(i)) then ! Still not bracketed.
2136  if (new_min_bound) then
2137  minent(i) = ent(i) ; f_minent(i) = f(i) ; df_de_min(i) = df_dent(i)
2138  else
2139  maxent(i) = ent(i) ; f_maxent(i) = f(i) ; df_de_max(i) = df_dent(i)
2140  endif
2141  endif
2142  else ! The root was previously bracketed.
2143  if (f(i) >= f_best(i)) then ! There is a new maximum.
2144  if (ent(i) > ent_best(i)) then ! Replace minent with ent_prev.
2145  minent(i) = ent_best(i) ; f_minent(i) = f_best(i) ; df_de_min(i) = df_de_best(i)
2146  else ! Replace maxent with ent_best.
2147  maxent(i) = ent_best(i) ; f_maxent(i) = f_best(i) ; df_de_max(i) = df_de_best(i)
2148  endif
2149  ent_best(i) = ent(i) ; f_best(i) = f(i) ; df_de_best(i) = df_dent(i)
2150  else
2151  if (ent(i) < ent_best(i)) then ! Replace the minent with ent.
2152  minent(i) = ent(i) ; f_minent(i) = f(i) ; df_de_min(i) = df_dent(i)
2153  else ! Replace maxent with ent_prev.
2154  maxent(i) = ent(i) ; f_maxent(i) = f(i) ; df_de_max(i) = df_dent(i)
2155  endif
2156  endif
2157  if ((maxent(i) - minent(i)) <= tolerance) do_i(i) = .false. ! Done.
2158  endif ! need_bracket.
2159  endif ; enddo
2160  if (.not.doany) exit
2161  enddo
2162  endif
2163 
2164  if (present(f_lim_maxent)) then
2165  ! Return the unlimited maximum in maxF, and the limited value of F at maxent.
2166  do i=is,ie
2167  maxf(i) = f_best(i)
2168  f_lim_maxent(i) = f_max_ent_in(i)
2169  if (present(ent_maxf)) ent_maxf(i) = ent_best(i)
2170  enddo
2171  else
2172  ! Now compare the two? potential maxima using the limited value of dF_kb.
2173  doany = .false.
2174  do i=is,ie
2175  may_use_best(i) = (ent_best(i) /= max_ent_in(i))
2176  if (may_use_best(i)) doany = .true.
2177  enddo
2178  if (doany) then
2179  ! For efficiency, could save previous value of dS_anom_lim_best?
2180  call determine_dskb(h_bl, sref, ent_bl, ent_best, is, ie, kmb, g, gv, .true., &
2181  ds_kb_lim)
2182  do i=is,ie
2183  f_best(i) = ent_best(i)*ds_kb_lim(i)*i_dskbp1(i)
2184  ! The second test seems necessary because of roundoff differences that
2185  ! can arise during compilation.
2186  if ((f_best(i) > f_max_ent_in(i)) .and. (may_use_best(i))) then
2187  maxf(i) = f_best(i)
2188  if (present(ent_maxf)) ent_maxf(i) = ent_best(i)
2189  else
2190  maxf(i) = f_max_ent_in(i)
2191  if (present(ent_maxf)) ent_maxf(i) = max_ent_in(i)
2192  endif
2193  enddo
2194  else
2195  ! All of the maxima are at the maximum entrainment.
2196  do i=is,ie ; maxf(i) = f_max_ent_in(i) ; enddo
2197  if (present(ent_maxf)) then
2198  do i=is,ie ; ent_maxf(i) = max_ent_in(i) ; enddo
2199  endif
2200  endif
2201  endif
2202 
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(szi_(g),szj_(g),szk_(g)), intent(in)  h,
real, dimension(szi_(g),szk_(g)+1), intent(in)  dtKd_int,
type(thermo_var_ptrs), intent(in)  tv,
integer, dimension(szi_(g)), intent(inout)  kb,
integer, intent(in)  kmb,
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,
integer, intent(in)  j,
real, dimension( g %isd: g %ied, g %ke+1), intent(out)  Ent_bl,
real, dimension( g %isd: g %ied, g %ke), intent(out)  Sref,
real, dimension( g %isd: g %ied, g %ke), 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, in H
[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]do_iA logical variable indicating which i-points to work on.
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 - 1000 for each layer, in kg m-3.
[out]h_blThe thickness of each layer, in H.

Definition at line 1091 of file MOM_entrain_diffusive.F90.

Referenced by entrainment_diffusive().

1091  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
1092  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
1093  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1094  intent(in) :: h !< Layer thicknesses, in H
1095  !! (usually m or kg m-2).
1096  real, dimension(SZI_(G),SZK_(G)+1), &
1097  intent(in) :: dtkd_int !< The diapycnal diffusivity across
1098  !! each interface times the time step,
1099  !! in H2.
1100  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
1101  !! available thermodynamic fields. Absent
1102  !! fields have NULL ptrs.
1103  integer, dimension(SZI_(G)), intent(inout) :: kb !< The index of the lightest layer denser
1104  !! than the buffer layer or 1 if there is
1105  !! no buffer layer.
1106  integer, intent(in) :: kmb
1107  logical, dimension(SZI_(G)), intent(in) :: do_i !< A logical variable indicating which
1108  !! i-points to work on.
1109  type(entrain_diffusive_cs), pointer :: cs !< This module's control structure.
1110  integer, intent(in) :: j !< The meridional index upon which
1111  !! to work.
1112  real, dimension(SZI_(G),SZK_(G)+1), &
1113  intent(out) :: ent_bl !< The average entrainment upward and
1114  !! downward across each interface around
1115  !! the buffer layers, in H.
1116  real, dimension(SZI_(G),SZK_(G)), intent(out) :: sref !< The coordinate potential density -
1117  !! 1000 for each layer, in kg m-3.
1118  real, dimension(SZI_(G),SZK_(G)), intent(out) :: h_bl !< The thickness of each layer, in H.
1119 
1120 ! Arguments: h - Layer thickness, in m or kg m-2 (abbreviated as H below).
1121 ! (in) dtKd_int - The diapycnal diffusivity across each interface times
1122 ! the time step, in H2.
1123 ! (in) tv - A structure containing pointers to any available
1124 ! thermodynamic fields. Absent fields have NULL ptrs.
1125 ! (in) kb - The index of the lightest layer denser than the
1126 ! buffer layer or 1 if there is no buffer layer.
1127 ! (in) do_i - A logical variable indicating which i-points to work on.
1128 ! (in) G - The ocean's grid structure.
1129 ! (in) GV - The ocean's vertical grid structure.
1130 ! (in) CS - This module's control structure.
1131 ! (in) j - The meridional index upon which to work.
1132 ! (out) Ent_bl - The average entrainment upward and downward across
1133 ! each interface around the buffer layers, in H.
1134 ! (out) Sref - The coordinate potential density - 1000 for each layer,
1135 ! in kg m-3.
1136 ! (out) h_bl - The thickness of each layer, in H.
1137 
1138 ! This subroutine sets the average entrainment across each of the interfaces
1139 ! between buffer layers within a timestep. It also causes thin and relatively
1140 ! light interior layers to be entrained by the deepest buffer layer.
1141 ! Also find the initial coordinate potential densities (Sref) of each layer.
1142 
1143 ! Does there need to be limiting when the layers below are all thin?
1144  real, dimension(SZI_(G)) :: &
1145  b1, d1, & ! Variables used by the tridiagonal solver, in H-1 and ND.
1146  rcv, & ! Value of the coordinate variable (potential density)
1147  ! based on the simulated T and S and P_Ref, kg m-3.
1148  pres, & ! Reference pressure (P_Ref) in Pa.
1149  frac_rem, & ! The fraction of the diffusion remaining, ND.
1150  h_interior ! The interior thickness available for entrainment, in H.
1151  real, dimension(SZI_(G), SZK_(G)) :: &
1152  s_est ! An estimate of the coordinate potential density - 1000 after
1153  ! entrainment for each layer, in kg m-3.
1154  real :: max_ent ! The maximum possible entrainment, in H.
1155  real :: dh ! An available thickness, in H.
1156  real :: kd_x_dt ! The diffusion that remains after thin layers are
1157  ! entrained, in H2.
1158  real :: h_neglect ! A thickness that is so small it is usually lost
1159  ! in roundoff and can be neglected, in H.
1160  integer :: i, k, is, ie, nz
1161  is = g%isc ; ie = g%iec ; nz = g%ke
1162 
1163 ! max_ent = 1.0e14*GV%Angstrom ! This is set to avoid roundoff problems.
1164  max_ent = 1.0e4*gv%m_to_H
1165  h_neglect = gv%H_subroundoff
1166 
1167  do i=is,ie ; pres(i) = tv%P_Ref ; enddo
1168  do k=1,kmb
1169  call calculate_density(tv%T(is:ie,j,k), tv%S(is:ie,j,k), pres(is:ie), &
1170  rcv(is:ie), 1, ie-is+1, tv%eqn_of_state)
1171  do i=is,ie
1172  h_bl(i,k) = h(i,j,k) + h_neglect
1173  sref(i,k) = rcv(i) - 1000.0
1174  enddo
1175  enddo
1176 
1177  do i=is,ie
1178  h_interior(i) = 0.0 ; ent_bl(i,1) = 0.0
1179 ! if (kb(i) > nz) Ent_bl(i,Kmb+1) = 0.0
1180  enddo
1181 
1182  do k=2,kmb ; do i=is,ie
1183  if (do_i(i)) then
1184  ent_bl(i,k) = min(2.0 * dtkd_int(i,k) / (h(i,j,k-1) + h(i,j,k) + h_neglect), &
1185  max_ent)
1186  else ; ent_bl(i,k) = 0.0 ; endif
1187  enddo ; enddo
1188 
1189  ! Determine the coordinate density of the bottommost buffer layer if there
1190  ! is no entrainment from the layers below. This is a partial solver, based
1191  ! on the first pass of a tridiagonal solver, as the values in the upper buffer
1192  ! layers are not needed.
1193 
1194  do i=is,ie
1195  b1(i) = 1.0 / (h_bl(i,1) + ent_bl(i,2))
1196  d1(i) = h_bl(i,1) * b1(i) ! = 1.0 - Ent_bl(i,2)*b1(i)
1197  s_est(i,1) = (h_bl(i,1)*sref(i,1)) * b1(i)
1198  enddo
1199  do k=2,kmb-1 ; do i=is,ie
1200  b1(i) = 1.0 / ((h_bl(i,k) + ent_bl(i,k+1)) + d1(i)*ent_bl(i,k))
1201  d1(i) = (h_bl(i,k) + d1(i)*ent_bl(i,k)) * b1(i) ! = 1.0 - Ent_bl(i,K+1)*b1(i)
1202  s_est(i,k) = (h_bl(i,k)*sref(i,k) + ent_bl(i,k)*s_est(i,k-1)) * b1(i)
1203  enddo ; enddo
1204  do i=is,ie
1205  s_est(i,kmb) = (h_bl(i,kmb)*sref(i,kmb) + ent_bl(i,kmb)*s_est(i,kmb-1)) / &
1206  (h_bl(i,kmb) + d1(i)*ent_bl(i,kmb))
1207  frac_rem(i) = 1.0
1208  enddo
1209 
1210  ! Entrain any thin interior layers that are lighter (in the coordinate
1211  ! potential density) than the deepest buffer layer will be, and adjust kb.
1212  do i=is,ie ; kb(i) = nz+1 ; if (do_i(i)) kb(i) = kmb+1 ; enddo
1213 
1214  do k=kmb+1,nz ; do i=is,ie ; if (do_i(i)) then
1215  if ((k == kb(i)) .and. (s_est(i,kmb) > (gv%Rlay(k) - 1000.0))) then
1216  if (4.0*dtkd_int(i,kmb+1)*frac_rem(i) > &
1217  (h_bl(i,kmb) + h(i,j,k)) * (h(i,j,k) - gv%Angstrom)) then
1218  ! Entrain this layer into the buffer layer and move kb down.
1219  dh = max((h(i,j,k) - gv%Angstrom), 0.0)
1220  if (dh > 0.0) then
1221  frac_rem(i) = frac_rem(i) - ((h_bl(i,kmb) + h(i,j,k)) * dh) / &
1222  (4.0*dtkd_int(i,kmb+1))
1223  sref(i,kmb) = (h_bl(i,kmb)*sref(i,kmb) + dh*(gv%Rlay(k)-1000.0)) / &
1224  (h_bl(i,kmb) + dh)
1225  h_bl(i,kmb) = h_bl(i,kmb) + dh
1226  s_est(i,kmb) = (h_bl(i,kmb)*sref(i,kmb) + ent_bl(i,kmb)*s_est(i,kmb-1)) / &
1227  (h_bl(i,kmb) + d1(i)*ent_bl(i,kmb))
1228  endif
1229  kb(i) = kb(i) + 1
1230  endif
1231  endif
1232  endif ; enddo ; enddo
1233 
1234  ! This is where variables are be set up with a different vertical grid
1235  ! in which the (newly?) massless layers are taken out.
1236  do k=nz,kmb+1,-1 ; do i=is,ie
1237  if (k >= kb(i)) h_interior(i) = h_interior(i) + (h(i,j,k)-gv%Angstrom)
1238  if (k==kb(i)) then
1239  h_bl(i,kmb+1) = h(i,j,k) ; sref(i,kmb+1) = gv%Rlay(k) - 1000.0
1240  elseif (k==kb(i)+1) then
1241  h_bl(i,kmb+2) = h(i,j,k) ; sref(i,kmb+2) = gv%Rlay(k) - 1000.0
1242  endif
1243  enddo ; enddo
1244  do i=is,ie ; if (kb(i) >= nz) then
1245  h_bl(i,kmb+1) = h(i,j,nz)
1246  sref(i,kmb+1) = gv%Rlay(nz) - 1000.0
1247  h_bl(i,kmb+2) = gv%Angstrom
1248  sref(i,kmb+2) = sref(i,kmb+1) + (gv%Rlay(nz) - gv%Rlay(nz-1))
1249  endif ; enddo
1250 
1251  ! Perhaps we should revisit the way that the average entrainment between the
1252  ! buffer layer and the interior is calculated so that it is not unduly
1253  ! limited when the layers are less than sqrt(Kd * dt) thick?
1254  do i=is,ie ; if (do_i(i)) then
1255  kd_x_dt = frac_rem(i) * dtkd_int(i,kmb+1)
1256  if ((kd_x_dt <= 0.0) .or. (h_interior(i) <= 0.0)) then
1257  ent_bl(i,kmb+1) = 0.0
1258  else
1259  ! If the combined layers are exceptionally thin, use sqrt(Kd*dt) as the
1260  ! estimate of the thickness in the denominator of the thickness diffusion.
1261  ent_bl(i,kmb+1) = min(0.5*h_interior(i), sqrt(kd_x_dt), &
1262  kd_x_dt / (0.5*(h_bl(i,kmb) + h_bl(i,kmb+1))))
1263  endif
1264  else
1265  ent_bl(i,kmb+1) = 0.0
1266  endif ; enddo
1267 
Here is the caller graph for this function: