MOM6
Data Types | Functions/Subroutines | Variables
mom_bulk_mixed_layer Module Reference

Data Types

type  bulkmixedlayer_cs
 

Functions/Subroutines

subroutine, public bulkmixedlayer (h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, CS, optics, Hml, aggregate_FW_forcing, dt_diag, last_call)
 This subroutine partially steps the bulk mixed layer model. The following processes are executed, in the order listed. More...
 
subroutine convective_adjustment (h, u, v, R0, Rcv, T, S, eps, d_eb, dKE_CA, cTKE, j, G, GV, CS, nz_conv)
 This subroutine does instantaneous convective entrainment into the buffer layers and mixed layers to remove hydrostatic instabilities. Any water that is lighter than currently in the mixed- or buffer- layer is entrained. More...
 
subroutine mixedlayer_convection (h, d_eb, htot, Ttot, Stot, uhtot, vhtot, R0_tot, Rcv_tot, u, v, T, S, R0, Rcv, eps, dR0_dT, dRcv_dT, dR0_dS, dRcv_dS, netMassInOut, netMassOut, Net_heat, Net_salt, nsw, Pen_SW_bnd, opacity_band, Conv_en, dKE_FC, j, ksort, G, GV, CS, tv, fluxes, dt, aggregate_FW_forcing)
 This subroutine causes the mixed layer to entrain to the depth of free convection. The depth of free convection is the shallowest depth at which the fluid is denser than the average of the fluid above. More...
 
subroutine find_starting_tke (htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, TKE, TKE_river, Idecay_len_TKE, cMKE, dt, Idt_diag, j, ksort, G, GV, CS)
 This subroutine determines the TKE available at the depth of free convection to drive mechanical entrainment. More...
 
subroutine mechanical_entrainment (h, d_eb, htot, Ttot, Stot, uhtot, vhtot, R0_tot, Rcv_tot, u, v, T, S, R0, Rcv, eps, dR0_dT, dRcv_dT, cMKE, Idt_diag, nsw, Pen_SW_bnd, opacity_band, TKE, Idecay_len_TKE, j, ksort, G, GV, CS)
 This subroutine calculates mechanically driven entrainment. More...
 
subroutine sort_ml (h, R0, eps, G, GV, CS, ksort)
 This subroutine generates an array of indices that are sorted by layer density. More...
 
subroutine resort_ml (h, T, S, R0, Rcv, RcvTgt, eps, d_ea, d_eb, ksort, G, GV, CS, dR0_dT, dR0_dS, dRcv_dT, dRcv_dS)
 This subroutine actually moves properties between layers to achieve a resorted state, with all of the resorted water either moved into the correct interior layers or in the top nkmb layers. More...
 
subroutine mixedlayer_detrain_2 (h, T, S, R0, Rcv, RcvTgt, dt, dt_diag, d_ea, j, G, GV, CS, dR0_dT, dR0_dS, dRcv_dT, dRcv_dS, max_BL_det)
 This subroutine moves any water left in the former mixed layers into the two buffer layers and may also move buffer layer water into the interior isopycnal layers. More...
 
subroutine mixedlayer_detrain_1 (h, T, S, R0, Rcv, RcvTgt, dt, dt_diag, d_ea, d_eb, j, G, GV, CS, dRcv_dT, dRcv_dS, max_BL_det)
 This subroutine moves any water left in the former mixed layers into the single buffer layers and may also move buffer layer water into the interior isopycnal layers. More...
 
subroutine, public bulkmixedlayer_init (Time, G, GV, param_file, diag, CS)
 
real function ef4 (H, E, L, dR_de)
 This subroutine returns an approximation to the integral R = exp(-L*(H+E)) integral(LH to L(H+E)) L/(1-(1+x)exp(-x)) dx. The approximation to the integrand is good to within -2% at x~.3 and +25% at x~3.5, but the exponential deemphasizes the importance of large x. When L=0, EF4 returns E/((H+E)*H). More...
 

Variables

integer id_clock_detrain =0
 
integer id_clock_mech =0
 
integer id_clock_conv =0
 
integer id_clock_adjustment =0
 
integer id_clock_eos =0
 
integer id_clock_resort =0
 
integer id_clock_pass =0
 
integer num_msg = 0
 
integer max_msg = 2
 

Function/Subroutine Documentation

◆ bulkmixedlayer()

subroutine, public mom_bulk_mixed_layer::bulkmixedlayer ( real, dimension(szi_(g),szj_(g),szk_(gv)), intent(inout)  h_3d,
real, dimension(szi_(g),szj_(g),szk_(gv)), intent(in)  u_3d,
real, dimension(szi_(g),szj_(g),szk_(gv)), intent(in)  v_3d,
type(thermo_var_ptrs), intent(inout)  tv,
type(forcing), intent(inout)  fluxes,
real, intent(in)  dt,
real, dimension(szi_(g),szj_(g),szk_(gv)), intent(inout)  ea,
real, dimension(szi_(g),szj_(g),szk_(gv)), intent(inout)  eb,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(bulkmixedlayer_cs), pointer  CS,
type(optics_type), pointer  optics,
real, dimension(:,:), pointer  Hml,
logical, intent(in)  aggregate_FW_forcing,
real, intent(in), optional  dt_diag,
logical, intent(in), optional  last_call 
)

This subroutine partially steps the bulk mixed layer model. The following processes are executed, in the order listed.

  1. Undergo convective adjustment into mixed layer.
  2. Apply surface heating and cooling.
  3. Starting from the top, entrain whatever fluid the TKE budget permits. Penetrating shortwave radiation is also applied at this point.
  4. If there is any unentrained fluid that was formerly in the mixed layer, detrain this fluid into the buffer layer. This is equivalent to the mixed layer detraining to the Monin- Obukhov depth.
  5. Divide the fluid in the mixed layer evenly into CSnkml pieces.
  6. Split the buffer layer if appropriate. Layers 1 to nkml are the mixed layer, nkml+1 to nkml+nkbl are the buffer layers. The results of this subroutine are mathematically identical if there are multiple pieces of the mixed layer with the same density or if there is just a single layer. There is no stability limit on the time step.

The key parameters for the mixed layer are found in the control structure. These include mstar, nstar, nstar2, pen_SW_frac, pen_SW_scale, and TKE_decay. For the Oberhuber (1993) mixed layer, the values of these are: pen_SW_frac = 0.42, pen_SW_scale = 15.0 m, mstar = 1.25, nstar = 1, TKE_decay = 2.5, conv_decay = 0.5 TKE_decay is 1/kappa in eq. 28 of Oberhuber (1993), while conv_decay is 1/mu. Conv_decay has been eliminated in favor of the well-calibrated form for the efficiency of penetrating convection from Wang (2003). For a traditional Kraus-Turner mixed layer, the values are: pen_SW_frac = 0.0, pen_SW_scale = 0.0 m, mstar = 1.25, nstar = 0.4, TKE_decay = 0.0, conv_decay = 0.0

Parameters
[in,out]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in,out]h_3dLayer thickness, in m or kg m-2.
[in]u_3dZonal velocities interpolated to h points,
[in]v_3dZonal velocities interpolated to h points,
[in,out]tvA structure containing pointers to any available thermodynamic fields. Absent fields have NULL ptrs.
[in,out]fluxesA structure containing pointers to any possible forcing fields. Unused fields have NULL ptrs.
[in]dtTime increment, in s.
[in,out]eaThe amount of fluid moved downward into a
[in,out]ebThe amount of fluid moved upward into a
csThe control structure returned by a previous call to mixedlayer_init.
opticsThe structure containing the inverse of the vertical absorption decay scale for penetrating shortwave radiation, in m-1.
hmlactive mixed layer depth
[in]dt_diagThe diagnostic time step, which may be less than dt if there are two callse to mixedlayer, in s.
[in]last_callif true, this is the last call to mixedlayer in the current time step, so diagnostics will be written. The default is .true.

Definition at line 227 of file MOM_bulk_mixed_layer.F90.

References mom_shortwave_abs::absorbremainingsw(), mom_eos::calculate_density_derivs(), convective_adjustment(), mom_diag_mediator::diag_update_remap_grids(), mom_domains::do_group_pass(), mom_forcing_type::extractfluxes1d(), find_starting_tke(), id_clock_adjustment, id_clock_conv, id_clock_detrain, id_clock_eos, id_clock_mech, id_clock_pass, id_clock_resort, mechanical_entrainment(), mixedlayer_convection(), mixedlayer_detrain_1(), mixedlayer_detrain_2(), mom_error_handler::mom_error(), resort_ml(), and sort_ml().

Referenced by mom_diabatic_driver::diabatic().

227  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
228  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
229  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
230  intent(inout) :: h_3d !< Layer thickness, in m or kg m-2.
231  !! (Intent in/out) The units of h are
232  !! referred to as H below.
233  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
234  intent(in) :: u_3d !< Zonal velocities interpolated to h points,
235  !! m s-1.
236  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
237  intent(in) :: v_3d !< Zonal velocities interpolated to h points,
238  !! m s-1.
239  type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers to any
240  !! available thermodynamic fields. Absent
241  !! fields have NULL ptrs.
242  type(forcing), intent(inout) :: fluxes !< A structure containing pointers to any
243  !! possible forcing fields. Unused fields
244  !! have NULL ptrs.
245  real, intent(in) :: dt !< Time increment, in s.
246  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
247  intent(inout) :: ea !< The amount of fluid moved downward into a
248  !! layer; this should be increased due to
249  !! mixed layer detrainment, in the same units
250  !! as h - usually m or kg m-2 (i.e., H).
251  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
252  intent(inout) :: eb !< The amount of fluid moved upward into a
253  !! layer; this should be increased due to
254  !! mixed layer entrainment, in the same units
255  !! as h - usually m or kg m-2 (i.e., H).
256  type(bulkmixedlayer_cs), pointer :: cs !< The control structure returned by a
257  !! previous call to mixedlayer_init.
258  type(optics_type), pointer :: optics !< The structure containing the inverse of the
259  !! vertical absorption decay scale for
260  !! penetrating shortwave radiation, in m-1.
261  real, dimension(:,:), pointer :: hml !< active mixed layer depth
262  logical, intent(in) :: aggregate_fw_forcing
263  real, optional, intent(in) :: dt_diag !< The diagnostic time step,
264  !! which may be less than dt if there are
265  !! two callse to mixedlayer, in s.
266  logical, optional, intent(in) :: last_call !< if true, this is the last call
267  !! to mixedlayer in the current time step, so
268  !! diagnostics will be written. The default is
269  !! .true.
270 
271 ! This subroutine partially steps the bulk mixed layer model.
272 ! The following processes are executed, in the order listed.
273 ! 1. Undergo convective adjustment into mixed layer.
274 ! 2. Apply surface heating and cooling.
275 ! 3. Starting from the top, entrain whatever fluid the TKE budget
276 ! permits. Penetrating shortwave radiation is also applied at
277 ! this point.
278 ! 4. If there is any unentrained fluid that was formerly in the
279 ! mixed layer, detrain this fluid into the buffer layer. This
280 ! is equivalent to the mixed layer detraining to the Monin-
281 ! Obukhov depth.
282 ! 5. Divide the fluid in the mixed layer evenly into CS%nkml pieces.
283 ! 6. Split the buffer layer if appropriate.
284 ! Layers 1 to nkml are the mixed layer, nkml+1 to nkml+nkbl are the
285 ! buffer layers. The results of this subroutine are mathematically
286 ! identical if there are multiple pieces of the mixed layer with
287 ! the same density or if there is just a single layer. There is no
288 ! stability limit on the time step.
289 !
290 ! The key parameters for the mixed layer are found in the control structure.
291 ! These include mstar, nstar, nstar2, pen_SW_frac, pen_SW_scale, and TKE_decay.
292 ! For the Oberhuber (1993) mixed layer, the values of these are:
293 ! pen_SW_frac = 0.42, pen_SW_scale = 15.0 m, mstar = 1.25,
294 ! nstar = 1, TKE_decay = 2.5, conv_decay = 0.5
295 ! TKE_decay is 1/kappa in eq. 28 of Oberhuber (1993), while conv_decay is 1/mu.
296 ! Conv_decay has been eliminated in favor of the well-calibrated form for the
297 ! efficiency of penetrating convection from Wang (2003).
298 ! For a traditional Kraus-Turner mixed layer, the values are:
299 ! pen_SW_frac = 0.0, pen_SW_scale = 0.0 m, mstar = 1.25,
300 ! nstar = 0.4, TKE_decay = 0.0, conv_decay = 0.0
301 
302 ! Arguments: h_3d - Layer thickness, in m or kg m-2. (Intent in/out)
303 ! The units of h are referred to as H below.
304 ! (in) u_3d - Zonal velocities interpolated to h points, m s-1.
305 ! (in) v_3d - Zonal velocities interpolated to h points, m s-1.
306 ! (in/out) tv - A structure containing pointers to any available
307 ! thermodynamic fields. Absent fields have NULL ptrs.
308 ! (in) fluxes - A structure containing pointers to any possible
309 ! forcing fields. Unused fields have NULL ptrs.
310 ! (in) dt - Time increment, in s.
311 ! (in/out) ea - The amount of fluid moved downward into a layer; this should
312 ! be increased due to mixed layer detrainment, in the same units
313 ! as h - usually m or kg m-2 (i.e., H).
314 ! (in/out) eb - The amount of fluid moved upward into a layer; this should
315 ! be increased due to mixed layer entrainment, in the same units
316 ! as h - usually m or kg m-2 (i.e., H).
317 ! (in) G - The ocean's grid structure.
318 ! (in) GV - The ocean's vertical grid structure.
319 ! (in) CS - The control structure returned by a previous call to
320 ! mixedlayer_init.
321 ! (in) optics - The structure containing the inverse of the vertical
322 ! absorption decay scale for penetrating shortwave
323 ! radiation, in m-1.
324 ! (in,opt) dt_diag - The diagnostic time step, which may be less than dt
325 ! if there are two callse to mixedlayer, in s.
326 ! (in,opt) last_call - if true, this is the last call to mixedlayer in the
327 ! current time step, so diagnostics will be written.
328 ! The default is .true.
329 
330  real, dimension(SZI_(G),SZK_(GV)) :: &
331  eaml, & ! The amount of fluid moved downward into a layer due to mixed
332  ! mixed layer detrainment, in m. (I.e. entrainment from above.)
333  ebml ! The amount of fluid moved upward into a layer due to mixed
334  ! mixed layer detrainment, in m. (I.e. entrainment from below.)
335 
336  ! If there is resorting, the vertical coordinate for these variables is the
337  ! new, sorted index space. Here layer 0 is an initially massless layer that
338  ! will be used to hold the new mixed layer properties.
339  real, dimension(SZI_(G),SZK0_(GV)) :: &
340  h, & ! The layer thickness, in m or kg m-2.
341  t, & ! The layer temperatures, in deg C.
342  s, & ! The layer salinities, in psu.
343  r0, & ! The potential density referenced to the surface, in kg m-3.
344  rcv ! The coordinate variable potential density, in kg m-3.
345  real, dimension(SZI_(G),SZK_(GV)) :: &
346  u, & ! The zonal velocity, in m s-1.
347  v, & ! The meridional velocity, in m s-1.
348  h_orig, & ! The original thickness in m or kg m-2.
349  d_eb, & ! The downward increase across a layer in the entrainment from
350  ! below, in H. The sign convention is that positive values of
351  ! d_eb correspond to a gain in mass by a layer by upward motion.
352  d_ea, & ! The upward increase across a layer in the entrainment from
353  ! above, in H. The sign convention is that positive values of
354  ! d_ea mean a net gain in mass by a layer from downward motion.
355  eps ! The (small) thickness that must remain in a layer, in H.
356  integer, dimension(SZI_(G),SZK_(GV)) :: &
357  ksort ! The sorted k-index that each original layer goes to.
358  real, dimension(SZI_(G),SZJ_(G)) :: &
359  h_miss ! The summed absolute mismatch, in m.
360  real, dimension(SZI_(G)) :: &
361  tke, & ! The turbulent kinetic energy available for mixing over a
362  ! time step, in m3 s-2.
363  conv_en, & ! The turbulent kinetic energy source due to mixing down to
364  ! the depth of free convection, in m3 s-2.
365  htot, & ! The total depth of the layers being considered for
366  ! entrainment, in H.
367  r0_tot, & ! The integrated potential density referenced to the surface
368  ! of the layers which are fully entrained, in H kg m-3.
369  rcv_tot, & ! The integrated coordinate value potential density of the
370  ! layers that are fully entrained, in H kg m-3.
371  ttot, & ! The integrated temperature of layers which are fully
372  ! entrained, in H K.
373  stot, & ! The integrated salt of layers which are fully entrained,
374  ! in H PSU.
375  uhtot, & ! The depth integrated zonal and meridional velocities in the
376  vhtot, & ! mixed layer, in H m s-1.
377 
378  netmassinout, & ! The net mass flux (if non-Boussinsq) or volume flux (if
379  ! Boussinesq - i.e. the fresh water flux (P+R-E)) into the
380  ! ocean over a time step, in H.
381  netmassout, & ! The mass flux (if non-Boussinsq) or volume flux (if
382  ! Boussinesq) over a time step from evaporating fresh water (H)
383  net_heat, & ! The net heating at the surface over a time step in K H. Any
384  ! penetrating shortwave radiation is not included in Net_heat.
385  net_salt, & ! The surface salt flux into the ocean over a time step, psu H.
386  idecay_len_tke, & ! The inverse of a turbulence decay length scale, in H-1.
387  p_ref, & ! Reference pressure for the potential density governing mixed
388  ! layer dynamics, almost always 0 (or 1e5) Pa.
389  p_ref_cv, & ! Reference pressure for the potential density which defines
390  ! the coordinate variable, set to P_Ref, in Pa.
391  dr0_dt, & ! Partial derivative of the mixed layer potential density with
392  ! temperature, in units of kg m-3 K-1.
393  drcv_dt, & ! Partial derivative of the coordinate variable potential
394  ! density in the mixed layer with temperature, in kg m-3 K-1.
395  dr0_ds, & ! Partial derivative of the mixed layer potential density with
396  ! salinity, in units of kg m-3 psu-1.
397  drcv_ds, & ! Partial derivative of the coordinate variable potential
398  ! density in the mixed layer with salinity, in kg m-3 psu-1.
399  tke_river ! The turbulent kinetic energy available for mixing at rivermouths over a
400  ! time step, in m3 s-2.
401 
402  real, dimension(max(CS%nsw,1),SZI_(G)) :: &
403  pen_sw_bnd ! The penetrating fraction of the shortwave heating integrated
404  ! over a time step in each band, in K H.
405  real, dimension(max(CS%nsw,1),SZI_(G),SZK_(GV)) :: &
406  opacity_band ! The opacity in each band, in H-1. The indicies are band, i, k.
407 
408  real :: cmke(2,szi_(g)) ! Coefficients of HpE and HpE^2 in calculating the
409  ! denominator of MKE_rate, in m-1 and m-2.
410  real :: irho0 ! 1.0 / rho_0
411  real :: inkml, inkmlm1! 1.0 / REAL(nkml) and 1.0 / REAL(nkml-1)
412  real :: ih ! The inverse of a thickness, in H-1.
413  real :: idt ! The inverse of the timestep in s-1.
414  real :: idt_diag ! The inverse of the timestep used for diagnostics in s-1.
415  real :: rmixconst
416 
417  real, dimension(SZI_(G)) :: &
418  dke_fc, & ! The change in mean kinetic energy due to free convection,
419  ! in m3 s-2.
420  h_ca ! The depth to which convective adjustment has gone in H.
421  real, dimension(SZI_(G),SZK_(GV)) :: &
422  dke_ca, & ! The change in mean kinetic energy due to convective
423  ! adjustment, in m3 s-2.
424  ctke ! The turbulent kinetic energy source due to convective
425  ! adjustment, m3 s-2.
426  real, dimension(SZI_(G),SZJ_(G)) :: &
427  hsfc_max, & ! The thickness of the surface region (mixed and buffer layers)
428  ! after entrainment but before any buffer layer detrainment, m.
429  hsfc_used, & ! The thickness of the surface region after buffer layer
430  ! detrainment, in units of m.
431  hsfc_min, & ! The minimum thickness of the surface region based on the
432  ! new mixed layer depth and the previous thickness of the
433  ! neighboring water columns, in m.
434  h_sum, & ! The total thickness of the water column, in H.
435  hmbl_prev ! The previous thickness of the mixed and buffer layers, in H.
436  real, dimension(SZI_(G)) :: &
437  hsfc, & ! The thickness of the surface region (mixed and buffer
438  ! layers before detrainment in to the interior, in H.
439  max_bl_det ! If non-negative, the maximum amount of entrainment from
440  ! the buffer layers that will be allowed this time step, in H.
441  real :: dhsfc, dhd ! Local copies of nondimensional parameters.
442  real :: h_nbr ! A minimum thickness based on neighboring thicknesses, in H.
443 
444  real :: absf_x_h ! The absolute value of f times the mixed layer thickness,
445  ! in units of m s-1.
446  real :: ku_star ! Ustar times the Von Karmen constant, in m s-1.
447  real :: dt__diag ! A copy of dt_diag (if present) or dt, in s.
448  logical :: write_diags ! If true, write out diagnostics with this step.
449  logical :: reset_diags ! If true, zero out the accumulated diagnostics.
450  integer :: i, j, k, is, ie, js, je, nz, nkmb, n
451  integer :: nsw ! The number of bands of penetrating shortwave radiation.
452 
453  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
454 
455  if (.not. associated(cs)) call mom_error(fatal, "MOM_mixed_layer: "//&
456  "Module must be initialized before it is used.")
457  if (gv%nkml < 1) return
458 
459  if (.not. ASSOCIATED(tv%eqn_of_state)) call mom_error(fatal, &
460  "MOM_mixed_layer: Temperature, salinity and an equation of state "//&
461  "must now be used.")
462  if (.NOT. ASSOCIATED(fluxes%ustar)) call mom_error(fatal, &
463  "MOM_mixed_layer: No surface TKE fluxes (ustar) defined in mixedlayer!")
464 
465  nkmb = cs%nkml+cs%nkbl
466  inkml = 1.0 / REAL(cs%nkml)
467  if (cs%nkml > 1) inkmlm1 = 1.0 / REAL(cs%nkml-1)
468 
469  irho0 = 1.0 / gv%Rho0
470  dt__diag = dt ; if (present(dt_diag)) dt__diag = dt_diag
471  idt = 1.0/dt
472  idt_diag = 1.0 / dt__diag
473  write_diags = .true. ; if (present(last_call)) write_diags = last_call
474 
475  p_ref(:) = 0.0 ; p_ref_cv(:) = tv%P_Ref
476 
477 
478  nsw = cs%nsw
479 
480  if (cs%limit_det .or. (cs%id_Hsfc_min > 0)) then
481 !$OMP parallel default(none) shared(is,ie,js,je,nkmb,h_sum,hmbl_prev,h_3d,nz)
482 !$OMP do
483  do j=js-1,je+1 ; do i=is-1,ie+1
484  h_sum(i,j) = 0.0 ; hmbl_prev(i,j) = 0.0
485  enddo ; enddo
486 !$OMP do
487  do j=js-1,je+1
488  do k=1,nkmb ; do i=is-1,ie+1
489  h_sum(i,j) = h_sum(i,j) + h_3d(i,j,k)
490  hmbl_prev(i,j) = hmbl_prev(i,j) + h_3d(i,j,k)
491  enddo ; enddo
492  do k=nkmb+1,nz ; do i=is-1,ie+1
493  h_sum(i,j) = h_sum(i,j) + h_3d(i,j,k)
494  enddo ; enddo
495  enddo
496 !$OMP end parallel
497 
498  call cpu_clock_begin(id_clock_pass)
499  call create_group_pass(cs%pass_h_sum_hmbl_prev, h_sum,g%Domain)
500  call create_group_pass(cs%pass_h_sum_hmbl_prev, hmbl_prev,g%Domain)
501  call do_group_pass(cs%pass_h_sum_hmbl_prev, g%Domain)
502  call cpu_clock_end(id_clock_pass)
503  endif
504 
505  ! Determine whether to zero out diagnostics before accumulation.
506  reset_diags = .true.
507  if (present(dt_diag) .and. write_diags .and. (dt__diag > dt)) &
508  reset_diags = .false. ! This is the second call to mixedlayer.
509 
510  if (reset_diags) then
511 !$OMP parallel default(none) shared(is,ie,js,je,CS)
512  if (cs%TKE_diagnostics) then
513 !$OMP do
514  do j=js,je ; do i=is,ie
515  cs%diag_TKE_wind(i,j) = 0.0 ; cs%diag_TKE_RiBulk(i,j) = 0.0
516  cs%diag_TKE_conv(i,j) = 0.0 ; cs%diag_TKE_pen_SW(i,j) = 0.0
517  cs%diag_TKE_mixing(i,j) = 0.0 ; cs%diag_TKE_mech_decay(i,j) = 0.0
518  cs%diag_TKE_conv_decay(i,j) = 0.0 ; cs%diag_TKE_conv_s2(i,j) = 0.0
519  enddo ; enddo
520  endif
521  if (ALLOCATED(cs%diag_PE_detrain)) then
522 !$OMP do
523  do j=js,je ; do i=is,ie
524  cs%diag_PE_detrain(i,j) = 0.0
525  enddo ; enddo
526  endif
527  if (ALLOCATED(cs%diag_PE_detrain2)) then
528 !$OMP do
529  do j=js,je ; do i=is,ie
530  cs%diag_PE_detrain2(i,j) = 0.0
531  enddo ; enddo
532  endif
533 !$OMP end parallel
534  endif
535 
536  if (cs%ML_resort) then
537  do i=is,ie ; h_ca(i) = 0.0 ; enddo
538  do k=1,nz ; do i=is,ie ; dke_ca(i,k) = 0.0 ; ctke(i,k) = 0.0 ; enddo ; enddo
539  endif
540  max_bl_det(:) = -1
541 
542 !$OMP parallel default(none) shared(is,ie,js,je,nz,h_3d,u_3d,v_3d,nkmb,G,GV,nsw,optics, &
543 !$OMP CS,tv,fluxes,Irho0,dt,Idt_diag,Ih,write_diags, &
544 !$OMP hmbl_prev,h_sum,Hsfc_min,Hsfc_max,dt__diag, &
545 !$OMP Hsfc_used,Inkmlm1,Inkml,ea,eb,h_miss,Hml, &
546 !$OMP id_clock_EOS,id_clock_resort,id_clock_adjustment, &
547 !$OMP id_clock_conv,id_clock_mech,id_clock_detrain,aggregate_FW_forcing ) &
548 !$OMP firstprivate(dKE_CA,cTKE,h_CA,max_BL_det,p_ref,p_ref_cv) &
549 !$OMP private(h,h_orig,u,v,eps,T,S,opacity_band,d_ea,d_eb, &
550 !$OMP dR0_dT,dR0_dS,dRcv_dT,dRcv_dS,R0,Rcv,ksort, &
551 !$OMP RmixConst,TKE_river,netMassInOut, NetMassOut, &
552 !$OMP Net_heat, Net_salt, htot,TKE,Pen_SW_bnd,Ttot,Stot, uhtot,&
553 !$OMP vhtot, R0_tot, Rcv_tot,Conv_en,dKE_FC,Idecay_len_TKE, &
554 !$OMP cMKE,Hsfc,dHsfc,dHD,H_nbr,kU_Star,absf_x_H, &
555 !$OMP ebml,eaml)
556 !$OMP do
557  do j=js,je
558  ! Copy the thicknesses and other fields to 2-d arrays.
559  do k=1,nz ; do i=is,ie
560  h(i,k) = h_3d(i,j,k) ; u(i,k) = u_3d(i,j,k) ; v(i,k) = v_3d(i,j,k)
561  h_orig(i,k) = h_3d(i,j,k)
562  eps(i,k) = 0.0 ; if (k > nkmb) eps(i,k) = gv%Angstrom
563  t(i,k) = tv%T(i,j,k) ; s(i,k) = tv%S(i,j,k)
564  do n=1,nsw
565  opacity_band(n,i,k) = gv%H_to_m*optics%opacity_band(n,i,j,k)
566  enddo
567  enddo ; enddo
568 
569  do k=1,nz ; do i=is,ie
570  d_ea(i,k) = 0.0 ; d_eb(i,k) = 0.0
571  enddo ; enddo
572 
573  if(id_clock_eos>0) call cpu_clock_begin(id_clock_eos)
574  ! Calculate an estimate of the mid-mixed layer pressure (in Pa)
575  do i=is,ie ; p_ref(i) = 0.0 ; enddo
576  do k=1,cs%nkml ; do i=is,ie
577  p_ref(i) = p_ref(i) + 0.5*gv%H_to_Pa*h(i,k)
578  enddo ; enddo
579  call calculate_density_derivs(t(:,1), s(:,1), p_ref, dr0_dt, dr0_ds, &
580  is, ie-is+1, tv%eqn_of_state)
581  call calculate_density_derivs(t(:,1), s(:,1), p_ref_cv, drcv_dt, drcv_ds, &
582  is, ie-is+1, tv%eqn_of_state)
583  do k=1,nz
584  call calculate_density(t(:,k), s(:,k), p_ref, r0(:,k), is, ie-is+1, &
585  tv%eqn_of_state)
586  call calculate_density(t(:,k), s(:,k), p_ref_cv, rcv(:,k), is, &
587  ie-is+1, tv%eqn_of_state)
588  enddo
589  if(id_clock_eos>0) call cpu_clock_end(id_clock_eos)
590 
591  if (cs%ML_resort) then
592  if(id_clock_resort>0) call cpu_clock_begin(id_clock_resort)
593  if (cs%ML_presort_nz_conv_adj > 0) &
594  call convective_adjustment(h(:,1:), u, v, r0(:,1:), rcv(:,1:), t(:,1:), &
595  s(:,1:), eps, d_eb, dke_ca, ctke, j, g, gv, cs, &
596  cs%ML_presort_nz_conv_adj)
597 
598  call sort_ml(h(:,1:), r0(:,1:), eps, g, gv, cs, ksort)
599  if(id_clock_resort>0) call cpu_clock_end(id_clock_resort)
600  else
601  do k=1,nz ; do i=is,ie ; ksort(i,k) = k ; enddo ; enddo
602 
603  if(id_clock_adjustment>0) call cpu_clock_begin(id_clock_adjustment)
604  ! Undergo instantaneous entrainment into the buffer layers and mixed layers
605  ! to remove hydrostatic instabilities. Any water that is lighter than
606  ! currently in the mixed or buffer layer is entrained.
607  call convective_adjustment(h(:,1:), u, v, r0(:,1:), rcv(:,1:), t(:,1:), &
608  s(:,1:), eps, d_eb, dke_ca, ctke, j, g, gv, cs)
609  do i=is,ie ; h_ca(i) = h(i,1) ; enddo
610 
611  if(id_clock_adjustment>0) call cpu_clock_end(id_clock_adjustment)
612  endif
613 
614  if (associated(fluxes%lrunoff) .and. cs%do_rivermix) then
615 
616  ! Here we add an additional source of TKE to the mixed layer where river
617  ! is present to simulate unresolved estuaries. The TKE input is diagnosed
618  ! as follows:
619  ! TKE_river[m3 s-3] = 0.5*rivermix_depth*g*Irho0*drho_ds*
620  ! River*(Samb - Sriver) = CS%mstar*U_star^3
621  ! where River is in units of m s-1.
622  ! Samb = Ambient salinity at the mouth of the estuary
623  ! rivermix_depth = The prescribed depth over which to mix river inflow
624  ! drho_ds = The gradient of density wrt salt at the ambient surface salinity.
625  ! Sriver = 0 (i.e. rivers are assumed to be pure freshwater)
626  rmixconst = 0.5*cs%rivermix_depth*gv%g_Earth*irho0**2
627  do i=is,ie
628  tke_river(i) = max(0.0, rmixconst*dr0_ds(i)* &
629  (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) * s(i,1))
630  enddo
631  else
632  do i=is,ie ; tke_river(i) = 0.0 ; enddo
633  endif
634 
635 
636  if(id_clock_conv>0) call cpu_clock_begin(id_clock_conv)
637 
638  ! The surface forcing is contained in the fluxes type.
639  ! We aggregate the thermodynamic forcing for a time step into the following:
640  ! netMassInOut = water (H units) added/removed via surface fluxes
641  ! netMassOut = water (H units) removed via evaporating surface fluxes
642  ! net_heat = heat (degC * H) via surface fluxes
643  ! net_salt = salt ( g(salt)/m2 for non-Bouss and ppt*m/s for Bouss ) via surface fluxes
644  ! Pen_SW_bnd = components to penetrative shortwave radiation
645  call extractfluxes1d(g, gv, fluxes, optics, nsw, j, dt, &
646  cs%H_limit_fluxes, cs%use_river_heat_content, cs%use_calving_heat_content, &
647  h(:,1:), t(:,1:), netmassinout, netmassout, net_heat, net_salt, pen_sw_bnd,&
648  tv, aggregate_fw_forcing)
649 
650  ! This subroutine causes the mixed layer to entrain to depth of free convection.
651  call mixedlayer_convection(h(:,1:), d_eb, htot, ttot, stot, uhtot, vhtot, &
652  r0_tot, rcv_tot, u, v, t(:,1:), s(:,1:), &
653  r0(:,1:), rcv(:,1:), eps, &
654  dr0_dt, drcv_dt, dr0_ds, drcv_ds, &
655  netmassinout, netmassout, net_heat, net_salt, &
656  nsw, pen_sw_bnd, opacity_band, conv_en, &
657  dke_fc, j, ksort, g, gv, cs, tv, fluxes, dt, &
658  aggregate_fw_forcing)
659 
660  if(id_clock_conv>0) call cpu_clock_end(id_clock_conv)
661 
662  ! Now the mixed layer undergoes mechanically forced entrainment.
663  ! The mixed layer may entrain down to the Monin-Obukhov depth if the
664  ! surface is becoming lighter, and is effectively detraining.
665 
666  ! First the TKE at the depth of free convection that is available
667  ! to drive mixing is calculated.
668  if(id_clock_mech>0) call cpu_clock_begin(id_clock_mech)
669 
670  call find_starting_tke(htot, h_ca, fluxes, conv_en, ctke, dke_fc, dke_ca, &
671  tke, tke_river, idecay_len_tke, cmke, dt, idt_diag, &
672  j, ksort, g, gv, cs)
673 
674  ! Here the mechanically driven entrainment occurs.
675  call mechanical_entrainment(h(:,1:), d_eb, htot, ttot, stot, uhtot, vhtot, &
676  r0_tot, rcv_tot, u, v, t(:,1:), s(:,1:), r0(:,1:), rcv(:,1:), eps, dr0_dt, drcv_dt, &
677  cmke, idt_diag, nsw, pen_sw_bnd, opacity_band, tke, &
678  idecay_len_tke, j, ksort, g, gv, cs)
679 
680  call absorbremainingsw(g, gv, h(:,1:), opacity_band, nsw, j, dt, cs%H_limit_fluxes, &
681  cs%correct_absorption, cs%absorb_all_SW, &
682  t(:,1:), pen_sw_bnd, eps, ksort, htot, ttot)
683 
684  if (cs%TKE_diagnostics) then ; do i=is,ie
685  cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) - idt_diag*tke(i)
686  enddo ; endif
687  if(id_clock_mech>0) call cpu_clock_end(id_clock_mech)
688 
689  ! Calculate the homogeneous mixed layer properties and store them in layer 0.
690  do i=is,ie ; if (htot(i) > 0.0) then
691  ih = 1.0 / htot(i)
692  r0(i,0) = r0_tot(i) * ih ; rcv(i,0) = rcv_tot(i) * ih
693  t(i,0) = ttot(i) * ih ; s(i,0) = stot(i) * ih
694  h(i,0) = htot(i)
695  else ! This may not ever be needed?
696  t(i,0) = t(i,1) ; s(i,0) = s(i,1) ; r0(i,0) = r0(i,1) ; rcv(i,0) = rcv(i,1)
697  h(i,0) = htot(i)
698  endif ; enddo
699  if (write_diags .and. ALLOCATED(cs%ML_depth)) then ; do i=is,ie
700  cs%ML_depth(i,j) = h(i,0) * gv%H_to_m
701  enddo ; endif
702  if (ASSOCIATED(hml)) then ; do i=is,ie
703  hml(i,j) = g%mask2dT(i,j) * (h(i,0) * gv%H_to_m)
704  enddo ; endif
705 
706 ! At this point, return water to the original layers, but constrained to
707 ! still be sorted. After this point, all the water that is in massive
708 ! interior layers will be denser than water remaining in the mixed- and
709 ! buffer-layers. To achieve this, some of these variable density layers
710 ! might be split between two isopycnal layers that are denser than new
711 ! mixed layer or any remaining water from the old mixed- or buffer-layers.
712 ! Alternately, if there are fewer than nkbl of the old buffer or mixed layers
713 ! with any mass, relatively light interior layers might be transferred to
714 ! these unused layers (but not currently in the code).
715 
716  if (cs%ML_resort) then
717  if(id_clock_resort>0) call cpu_clock_begin(id_clock_resort)
718  call resort_ml(h(:,0:), t(:,0:), s(:,0:), r0(:,0:), rcv(:,0:), gv%Rlay, eps, &
719  d_ea, d_eb, ksort, g, gv, cs, dr0_dt, dr0_ds, drcv_dt, drcv_ds)
720  if(id_clock_resort>0) call cpu_clock_end(id_clock_resort)
721  endif
722 
723  if (cs%limit_det .or. (cs%id_Hsfc_max > 0) .or. (cs%id_Hsfc_min > 0)) then
724  do i=is,ie ; hsfc(i) = h(i,0) ; enddo
725  do k=1,nkmb ; do i=is,ie ; hsfc(i) = hsfc(i) + h(i,k) ; enddo ; enddo
726 
727  if (cs%limit_det .or. (cs%id_Hsfc_min > 0)) then
728  dhsfc = cs%lim_det_dH_sfc ; dhd = cs%lim_det_dH_bathy
729  do i=is,ie
730  h_nbr = min(dhsfc*max(hmbl_prev(i-1,j), hmbl_prev(i+1,j), &
731  hmbl_prev(i,j-1), hmbl_prev(i,j+1)), &
732  max(hmbl_prev(i-1,j) - dhd*min(h_sum(i,j),h_sum(i-1,j)), &
733  hmbl_prev(i+1,j) - dhd*min(h_sum(i,j),h_sum(i+1,j)), &
734  hmbl_prev(i,j-1) - dhd*min(h_sum(i,j),h_sum(i,j-1)), &
735  hmbl_prev(i,j+1) - dhd*min(h_sum(i,j),h_sum(i,j+1))) )
736 
737  hsfc_min(i,j) = gv%H_to_m*max(h(i,0), min(hsfc(i), h_nbr))
738 
739  if (cs%limit_det) max_bl_det(i) = max(0.0, hsfc(i)-h_nbr)
740  enddo
741  endif
742 
743  if (cs%id_Hsfc_max > 0) then ; do i=is,ie
744  hsfc_max(i,j) = hsfc(i)*gv%H_to_m
745  enddo ; endif
746  endif
747 
748 ! Move water left in the former mixed layer into the buffer layer and
749 ! from the buffer layer into the interior. These steps might best be
750 ! treated in conjuction.
751  if(id_clock_detrain>0) call cpu_clock_begin(id_clock_detrain)
752  if (cs%nkbl == 1) then
753  call mixedlayer_detrain_1(h(:,0:), t(:,0:), s(:,0:), r0(:,0:), rcv(:,0:), &
754  gv%Rlay, dt, dt__diag, d_ea, d_eb, j, g, gv, cs, &
755  drcv_dt, drcv_ds, max_bl_det)
756  elseif (cs%nkbl == 2) then
757  call mixedlayer_detrain_2(h(:,0:), t(:,0:), s(:,0:), r0(:,0:), rcv(:,0:), &
758  gv%Rlay, dt, dt__diag, d_ea, j, g, gv, cs, &
759  dr0_dt, dr0_ds, drcv_dt, drcv_ds, max_bl_det)
760  else ! CS%nkbl not = 1 or 2
761  ! This code only works with 1 or 2 buffer layers.
762  call mom_error(fatal, "MOM_mixed_layer: CS%nkbl must be 1 or 2 for now.")
763  endif
764  if(id_clock_detrain>0) call cpu_clock_end(id_clock_detrain)
765 
766 
767  if (cs%id_Hsfc_used > 0) then
768  do i=is,ie ; hsfc_used(i,j) = h(i,0)*gv%H_to_m ; enddo
769  do k=cs%nkml+1,nkmb ; do i=is,ie
770  hsfc_used(i,j) = hsfc_used(i,j) + h(i,k)*gv%H_to_m
771  enddo ; enddo
772  endif
773 
774 ! Now set the properties of the layers in the mixed layer in the original
775 ! 3-d variables.
776  if (cs%Resolve_Ekman .and. (cs%nkml>1)) then
777  ! The thickness of the topmost piece of the mixed layer is given by
778  ! h_1 = H / (3 + sqrt(|f|*H^2/2*nu_max)), which asymptotes to the Ekman
779  ! layer depth and 1/3 of the mixed layer depth. This curve has been
780  ! determined to maximize the impact of the Ekman transport in the mixed
781  ! layer TKE budget with nkml=2. With nkml=3, this should also be used,
782  ! as the third piece will then optimally describe mixed layer
783  ! restratification. For nkml>=4 the whole strategy should be revisited.
784  do i=is,ie
785  ku_star = 0.41*fluxes%ustar(i,j) ! Maybe could be replaced with u*+w*?
786  if (associated(fluxes%ustar_shelf) .and. &
787  associated(fluxes%frac_shelf_h)) then
788  if (fluxes%frac_shelf_h(i,j) > 0.0) &
789  ku_star = (1.0 - fluxes%frac_shelf_h(i,j)) * ku_star + &
790  fluxes%frac_shelf_h(i,j) * (0.41*fluxes%ustar_shelf(i,j))
791  endif
792  absf_x_h = 0.25 * h(i,0) * &
793  ((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
794  (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i-1,j))))
795  ! If the mixed layer vertical viscosity specification is changed in
796  ! MOM_vert_friction.F90, this line will have to be modified accordingly.
797  h_3d(i,j,1) = h(i,0) / (3.0 + sqrt(absf_x_h*(absf_x_h + 2.0*ku_star) / &
798  (ku_star**2)) )
799  do k=2,cs%nkml
800  ! The other layers are evenly distributed through the mixed layer.
801  h_3d(i,j,k) = (h(i,0)-h_3d(i,j,1)) * inkmlm1
802  d_ea(i,k) = d_ea(i,k) + h_3d(i,j,k)
803  d_ea(i,1) = d_ea(i,1) - h_3d(i,j,k)
804  enddo
805  enddo
806  else
807  do i=is,ie
808  h_3d(i,j,1) = h(i,0) * inkml
809  enddo
810  do k=2,cs%nkml ; do i=is,ie
811  h_3d(i,j,k) = h(i,0) * inkml
812  d_ea(i,k) = d_ea(i,k) + h_3d(i,j,k)
813  d_ea(i,1) = d_ea(i,1) - h_3d(i,j,k)
814  enddo ; enddo
815  endif
816  do i=is,ie ; h(i,0) = 0.0 ; enddo
817  do k=1,cs%nkml ; do i=is,ie
818  tv%T(i,j,k) = t(i,0) ; tv%S(i,j,k) = s(i,0)
819  enddo ; enddo
820 
821  ! These sum needs to be done in the original layer space.
822 
823  ! The treatment of layer 1 is atypical because evaporation shows up as
824  ! negative ea(i,1), and because all precipitation goes straight into layer 1.
825  ! The code is ordered so that any roundoff errors in ea are lost the surface.
826 ! do i=is,ie ; eaml(i,1) = 0.0 ; enddo
827 ! do k=2,nz ; do i=is,ie ; eaml(i,k) = eaml(i,k-1) - d_ea(i,k-1) ; enddo ; enddo
828 ! do i=is,ie ; eaml(i,1) = netMassInOut(i) ; enddo
829 
830 
831  do i=is,ie
832 ! eaml(i,nz) is derived from h(i,nz) - h_orig(i,nz) = eaml(i,nz) - ebml(i,nz-1)
833  ebml(i,nz) = 0.0
834  eaml(i,nz) = (h(i,nz) - h_orig(i,nz)) - d_eb(i,nz)
835  enddo
836  do k=nz-1,1,-1 ; do i=is,ie
837  ebml(i,k) = ebml(i,k+1) - d_eb(i,k+1)
838  eaml(i,k) = eaml(i,k+1) + d_ea(i,k)
839  enddo ; enddo
840  do i=is,ie ; eaml(i,1) = netmassinout(i) ; enddo
841 
842  ! Copy the interior thicknesses and other fields back to the 3-d arrays.
843  do k=cs%nkml+1,nz ; do i=is,ie
844  h_3d(i,j,k) = h(i,k); tv%T(i,j,k) = t(i,k) ; tv%S(i,j,k) = s(i,k)
845  enddo ; enddo
846 
847  do k=1,nz ; do i=is,ie
848  ea(i,j,k) = ea(i,j,k) + eaml(i,k)
849  eb(i,j,k) = eb(i,j,k) + ebml(i,k)
850  enddo ; enddo
851 
852  if (cs%id_h_mismatch > 0) then
853  do i=is,ie
854  h_miss(i,j) = abs(h_3d(i,j,1) - (h_orig(i,1) + &
855  (eaml(i,1) + (ebml(i,1) - eaml(i,1+1)))))
856  enddo
857  do k=2,nz-1 ; do i=is,ie
858  h_miss(i,j) = h_miss(i,j) + abs(h_3d(i,j,k) - (h_orig(i,k) + &
859  ((eaml(i,k) - ebml(i,k-1)) + (ebml(i,k) - eaml(i,k+1)))))
860  enddo ; enddo
861  do i=is,ie
862  h_miss(i,j) = h_miss(i,j) + abs(h_3d(i,j,nz) - (h_orig(i,nz) + &
863  ((eaml(i,nz) - ebml(i,nz-1)) + ebml(i,nz))))
864  h_miss(i,j) = gv%H_to_m * h_miss(i,j)
865  enddo
866  endif
867 
868  enddo ! j loop
869 
870  ! Whenever thickness changes let the diag manager know, target grids
871  ! for vertical remapping may need to be regenerated.
872  ! This needs to happen after the H update and before the next post_data.
873  call diag_update_remap_grids(cs%diag)
874 
875 !$OMP end parallel
876 
877  if (write_diags) then
878  if (cs%id_ML_depth > 0) &
879  call post_data(cs%id_ML_depth, cs%ML_depth, cs%diag)
880  if (cs%id_TKE_wind > 0) &
881  call post_data(cs%id_TKE_wind, cs%diag_TKE_wind, cs%diag)
882  if (cs%id_TKE_RiBulk > 0) &
883  call post_data(cs%id_TKE_RiBulk, cs%diag_TKE_RiBulk, cs%diag)
884  if (cs%id_TKE_conv > 0) &
885  call post_data(cs%id_TKE_conv, cs%diag_TKE_conv, cs%diag)
886  if (cs%id_TKE_pen_SW > 0) &
887  call post_data(cs%id_TKE_pen_SW, cs%diag_TKE_pen_SW, cs%diag)
888  if (cs%id_TKE_mixing > 0) &
889  call post_data(cs%id_TKE_mixing, cs%diag_TKE_mixing, cs%diag)
890  if (cs%id_TKE_mech_decay > 0) &
891  call post_data(cs%id_TKE_mech_decay, cs%diag_TKE_mech_decay, cs%diag)
892  if (cs%id_TKE_conv_decay > 0) &
893  call post_data(cs%id_TKE_conv_decay, cs%diag_TKE_conv_decay, cs%diag)
894  if (cs%id_TKE_conv_s2 > 0) &
895  call post_data(cs%id_TKE_conv_s2, cs%diag_TKE_conv_s2, cs%diag)
896  if (cs%id_PE_detrain > 0) &
897  call post_data(cs%id_PE_detrain, cs%diag_PE_detrain, cs%diag)
898  if (cs%id_PE_detrain2 > 0) &
899  call post_data(cs%id_PE_detrain2, cs%diag_PE_detrain2, cs%diag)
900  if (cs%id_h_mismatch > 0) &
901  call post_data(cs%id_h_mismatch, h_miss, cs%diag)
902  if (cs%id_Hsfc_used > 0) &
903  call post_data(cs%id_Hsfc_used, hsfc_used, cs%diag)
904  if (cs%id_Hsfc_max > 0) &
905  call post_data(cs%id_Hsfc_max, hsfc_max, cs%diag)
906  if (cs%id_Hsfc_min > 0) &
907  call post_data(cs%id_Hsfc_min, hsfc_min, cs%diag)
908  endif
909 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ bulkmixedlayer_init()

subroutine, public mom_bulk_mixed_layer::bulkmixedlayer_init ( type(time_type), intent(in), target  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(bulkmixedlayer_cs), pointer  CS 
)
Parameters
[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 for this module.

Definition at line 3647 of file MOM_bulk_mixed_layer.F90.

References id_clock_adjustment, id_clock_conv, id_clock_detrain, id_clock_eos, id_clock_mech, id_clock_pass, and id_clock_resort.

Referenced by mom_diabatic_driver::diabatic_driver_init().

3647  type(time_type), target, intent(in) :: time
3648  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
3649  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
3650  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
3651  !! parameters.
3652  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
3653  !! output.
3654  type(bulkmixedlayer_cs), pointer :: cs !< A pointer that is set to point to the control
3655  !! structure for this module.
3656 ! Arguments: Time - The current model time.
3657 ! (in) G - The ocean's grid structure.
3658 ! (in) GV - The ocean's vertical grid structure.
3659 ! (in) param_file - A structure indicating the open file to parse for
3660 ! model parameter values.
3661 ! (in) diag - A structure that is used to regulate diagnostic output.
3662 ! (in/out) CS - A pointer that is set to point to the control structure
3663 ! for this module
3664 ! This include declares and sets the variable "version".
3665 #include "version_variable.h"
3666  character(len=40) :: mdl = "MOM_mixed_layer" ! This module's name.
3667  real :: omega_frac_dflt, ustar_min_dflt
3668  integer :: isd, ied, jsd, jed
3669  logical :: use_temperature, use_omega
3670  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3671 
3672  if (associated(cs)) then
3673  call mom_error(warning, "mixedlayer_init called with an associated"// &
3674  "associated control structure.")
3675  return
3676  else ; allocate(cs) ; endif
3677 
3678  cs%diag => diag
3679  cs%Time => time
3680 
3681  if (gv%nkml < 1) return
3682 
3683 ! Set default, read and log parameters
3684  call log_version(param_file, mdl, version, "")
3685 
3686  cs%nkml = gv%nkml
3687  call log_param(param_file, mdl, "NKML", cs%nkml, &
3688  "The number of sublayers within the mixed layer if \n"//&
3689  "BULKMIXEDLAYER is true.", units="nondim", default=2)
3690  cs%nkbl = gv%nk_rho_varies - gv%nkml
3691  call log_param(param_file, mdl, "NKBL", cs%nkbl, &
3692  "The number of variable density buffer layers if \n"//&
3693  "BULKMIXEDLAYER is true.", units="nondim", default=2)
3694  call get_param(param_file, mdl, "MSTAR", cs%mstar, &
3695  "The ratio of the friction velocity cubed to the TKE \n"//&
3696  "input to the mixed layer.", "units=nondim", default=1.2)
3697  call get_param(param_file, mdl, "NSTAR", cs%nstar, &
3698  "The portion of the buoyant potential energy imparted by \n"//&
3699  "surface fluxes that is available to drive entrainment \n"//&
3700  "at the base of mixed layer when that energy is positive.", &
3701  units="nondim", default=0.15)
3702  call get_param(param_file, mdl, "BULK_RI_ML", cs%bulk_Ri_ML, &
3703  "The efficiency with which mean kinetic energy released \n"//&
3704  "by mechanically forced entrainment of the mixed layer \n"//&
3705  "is converted to turbulent kinetic energy.", units="nondim",&
3706  fail_if_missing=.true.)
3707  call get_param(param_file, mdl, "ABSORB_ALL_SW", cs%absorb_all_sw, &
3708  "If true, all shortwave radiation is absorbed by the \n"//&
3709  "ocean, instead of passing through to the bottom mud.", &
3710  default=.false.)
3711  call get_param(param_file, mdl, "TKE_DECAY", cs%TKE_decay, &
3712  "TKE_DECAY relates the vertical rate of decay of the \n"//&
3713  "TKE available for mechanical entrainment to the natural \n"//&
3714  "Ekman depth.", units="nondim", default=2.5)
3715  call get_param(param_file, mdl, "NSTAR2", cs%nstar2, &
3716  "The portion of any potential energy released by \n"//&
3717  "convective adjustment that is available to drive \n"//&
3718  "entrainment at the base of mixed layer. By default \n"//&
3719  "NSTAR2=NSTAR.", units="nondim", default=cs%nstar)
3720  call get_param(param_file, mdl, "BULK_RI_CONVECTIVE", cs%bulk_Ri_convective, &
3721  "The efficiency with which convectively released mean \n"//&
3722  "kinetic energy is converted to turbulent kinetic \n"//&
3723  "energy. By default BULK_RI_CONVECTIVE=BULK_RI_ML.", &
3724  units="nondim", default=cs%bulk_Ri_ML)
3725  call get_param(param_file, mdl, "HMIX_MIN", cs%Hmix_min, &
3726  "The minimum mixed layer depth if the mixed layer depth \n"//&
3727  "is determined dynamically.", units="m", default=0.0)
3728 
3729  call get_param(param_file, mdl, "LIMIT_BUFFER_DETRAIN", cs%limit_det, &
3730  "If true, limit the detrainment from the buffer layers \n"//&
3731  "to not be too different from the neighbors.", default=.false.)
3732  call get_param(param_file, mdl, "ALLOWED_DETRAIN_TEMP_CHG", cs%Allowed_T_chg, &
3733  "The amount by which temperature is allowed to exceed \n"//&
3734  "previous values during detrainment.", units="K", default=0.5)
3735  call get_param(param_file, mdl, "ALLOWED_DETRAIN_SALT_CHG", cs%Allowed_S_chg, &
3736  "The amount by which salinity is allowed to exceed \n"//&
3737  "previous values during detrainment.", units="PSU", default=0.1)
3738  call get_param(param_file, mdl, "ML_DT_DS_WEIGHT", cs%dT_dS_wt, &
3739  "When forced to extrapolate T & S to match the layer \n"//&
3740  "densities, this factor (in deg C / PSU) is combined \n"//&
3741  "with the derivatives of density with T & S to determine \n"//&
3742  "what direction is orthogonal to density contours. It \n"//&
3743  "should be a typical value of (dR/dS) / (dR/dT) in \n"//&
3744  "oceanic profiles.", units="degC PSU-1", default=6.0)
3745  call get_param(param_file, mdl, "BUFFER_LAYER_EXTRAP_LIMIT", cs%BL_extrap_lim, &
3746  "A limit on the density range over which extrapolation \n"//&
3747  "can occur when detraining from the buffer layers, \n"//&
3748  "relative to the density range within the mixed and \n"//&
3749  "buffer layers, when the detrainment is going into the \n"//&
3750  "lightest interior layer, nondimensional, or a negative \n"//&
3751  "value not to apply this limit.", units="nondim", default = -1.0)
3752  call get_param(param_file, mdl, "DEPTH_LIMIT_FLUXES", cs%H_limit_fluxes, &
3753  "The surface fluxes are scaled away when the total ocean \n"//&
3754  "depth is less than DEPTH_LIMIT_FLUXES.", &
3755  units="m", default=0.1*cs%Hmix_min)
3756  call get_param(param_file, mdl, "OMEGA",cs%omega, &
3757  "The rotation rate of the earth.", units="s-1", &
3758  default=7.2921e-5)
3759  call get_param(param_file, mdl, "ML_USE_OMEGA", use_omega, &
3760  "If true, use the absolute rotation rate instead of the \n"//&
3761  "vertical component of rotation when setting the decay \n"//&
3762  "scale for turbulence.", default=.false., do_not_log=.true.)
3763  omega_frac_dflt = 0.0
3764  if (use_omega) then
3765  call mom_error(warning, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
3766  omega_frac_dflt = 1.0
3767  endif
3768  call get_param(param_file, mdl, "ML_OMEGA_FRAC", cs%omega_frac, &
3769  "When setting the decay scale for turbulence, use this \n"//&
3770  "fraction of the absolute rotation rate blended with the \n"//&
3771  "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
3772  units="nondim", default=omega_frac_dflt)
3773  call get_param(param_file, mdl, "ML_RESORT", cs%ML_resort, &
3774  "If true, resort the topmost layers by potential density \n"//&
3775  "before the mixed layer calculations.", default=.false.)
3776  if (cs%ML_resort) &
3777  call get_param(param_file, mdl, "ML_PRESORT_NK_CONV_ADJ", cs%ML_presort_nz_conv_adj, &
3778  "Convectively mix the first ML_PRESORT_NK_CONV_ADJ \n"//&
3779  "layers before sorting when ML_RESORT is true.", &
3780  units="nondim", default=0, fail_if_missing=.true.) ! Fail added by AJA.
3781  ! This gives a minimum decay scale that is typically much less than Angstrom.
3782  ustar_min_dflt = 2e-4*cs%omega*(gv%Angstrom_z + gv%H_to_m*gv%H_subroundoff)
3783  call get_param(param_file, mdl, "BML_USTAR_MIN", cs%ustar_min, &
3784  "The minimum value of ustar that should be used by the \n"//&
3785  "bulk mixed layer model in setting vertical TKE decay \n"//&
3786  "scales. This must be greater than 0.", units="m s-1", &
3787  default=ustar_min_dflt)
3788  if (cs%ustar_min<=0.0) call mom_error(fatal, "BML_USTAR_MIN must be positive.")
3789 
3790  call get_param(param_file, mdl, "RESOLVE_EKMAN", cs%Resolve_Ekman, &
3791  "If true, the NKML>1 layers in the mixed layer are \n"//&
3792  "chosen to optimally represent the impact of the Ekman \n"//&
3793  "transport on the mixed layer TKE budget. Otherwise, \n"//&
3794  "the sublayers are distributed uniformly through the \n"//&
3795  "mixed layer.", default=.false.)
3796  call get_param(param_file, mdl, "CORRECT_ABSORPTION_DEPTH", cs%correct_absorption, &
3797  "If true, the average depth at which penetrating shortwave \n"//&
3798  "radiation is absorbed is adjusted to match the average \n"//&
3799  "heating depth of an exponential profile by moving some \n"//&
3800  "of the heating upward in the water column.", default=.false.)
3801  call get_param(param_file, mdl, "DO_RIVERMIX", cs%do_rivermix, &
3802  "If true, apply additional mixing whereever there is \n"//&
3803  "runoff, so that it is mixed down to RIVERMIX_DEPTH, \n"//&
3804  "if the ocean is that deep.", default=.false.)
3805  if (cs%do_rivermix) &
3806  call get_param(param_file, mdl, "RIVERMIX_DEPTH", cs%rivermix_depth, &
3807  "The depth to which rivers are mixed if DO_RIVERMIX is \n"//&
3808  "defined.", units="m", default=0.0)
3809  call get_param(param_file, mdl, "USE_RIVER_HEAT_CONTENT", cs%use_river_heat_content, &
3810  "If true, use the fluxes%runoff_Hflx field to set the \n"//&
3811  "heat carried by runoff, instead of using SST*CP*liq_runoff.", &
3812  default=.false.)
3813  call get_param(param_file, mdl, "USE_CALVING_HEAT_CONTENT", cs%use_calving_heat_content, &
3814  "If true, use the fluxes%calving_Hflx field to set the \n"//&
3815  "heat carried by runoff, instead of using SST*CP*froz_runoff.", &
3816  default=.false.)
3817 
3818  call get_param(param_file, mdl, "ALLOW_CLOCKS_IN_OMP_LOOPS", &
3819  cs%allow_clocks_in_omp_loops, &
3820  "If true, clocks can be called from inside loops that can \n"//&
3821  "be threaded. To run with multiple threads, set to False.", &
3822  default=.true.)
3823 
3824  cs%id_ML_depth = register_diag_field('ocean_model', 'h_ML', diag%axesT1, &
3825  time, 'Surface mixed layer depth', 'meter')
3826  cs%id_TKE_wind = register_diag_field('ocean_model', 'TKE_wind', diag%axesT1, &
3827  time, 'Wind-stirring source of mixed layer TKE', 'meter3 second-3')
3828  cs%id_TKE_RiBulk = register_diag_field('ocean_model', 'TKE_RiBulk', diag%axesT1, &
3829  time, 'Mean kinetic energy source of mixed layer TKE', 'meter3 second-3')
3830  cs%id_TKE_conv = register_diag_field('ocean_model', 'TKE_conv', diag%axesT1, &
3831  time, 'Convective source of mixed layer TKE', 'meter3 second-3')
3832  cs%id_TKE_pen_SW = register_diag_field('ocean_model', 'TKE_pen_SW', diag%axesT1, &
3833  time, 'TKE consumed by mixing penetrative shortwave radation through the mixed layer', 'meter3 second-3')
3834  cs%id_TKE_mixing = register_diag_field('ocean_model', 'TKE_mixing', diag%axesT1, &
3835  time, 'TKE consumed by mixing that deepens the mixed layer', 'meter3 second-3')
3836  cs%id_TKE_mech_decay = register_diag_field('ocean_model', 'TKE_mech_decay', diag%axesT1, &
3837  time, 'Mechanical energy decay sink of mixed layer TKE', 'meter3 second-3')
3838  cs%id_TKE_conv_decay = register_diag_field('ocean_model', 'TKE_conv_decay', diag%axesT1, &
3839  time, 'Convective energy decay sink of mixed layer TKE', 'meter3 second-3')
3840  cs%id_TKE_conv_s2 = register_diag_field('ocean_model', 'TKE_conv_s2', diag%axesT1, &
3841  time, 'Spurious source of mixed layer TKE from sigma2', 'meter3 second-3')
3842  cs%id_PE_detrain = register_diag_field('ocean_model', 'PE_detrain', diag%axesT1, &
3843  time, 'Spurious source of potential energy from mixed layer detrainment', 'Watt meter-2')
3844  cs%id_PE_detrain2 = register_diag_field('ocean_model', 'PE_detrain2', diag%axesT1, &
3845  time, 'Spurious source of potential energy from mixed layer only detrainment', 'Watt meter-2')
3846  cs%id_h_mismatch = register_diag_field('ocean_model', 'h_miss_ML', diag%axesT1, &
3847  time, 'Summed absolute mismatch in entrainment terms', 'meter')
3848  cs%id_Hsfc_used = register_diag_field('ocean_model', 'Hs_used', diag%axesT1, &
3849  time, 'Surface region thickness that is used', 'meter')
3850  cs%id_Hsfc_max = register_diag_field('ocean_model', 'Hs_max', diag%axesT1, &
3851  time, 'Maximum surface region thickness', 'meter')
3852  cs%id_Hsfc_min = register_diag_field('ocean_model', 'Hs_min', diag%axesT1, &
3853  time, 'Minimum surface region thickness', 'meter')
3854  !CS%lim_det_dH_sfc = 0.5 ; CS%lim_det_dH_bathy = 0.2 ! Technically these should not get used if limit_det is false?
3855  if (cs%limit_det .or. (cs%id_Hsfc_min > 0)) then
3856  call get_param(param_file, mdl, "LIMIT_BUFFER_DET_DH_SFC", cs%lim_det_dH_sfc, &
3857  "The fractional limit in the change between grid points \n"//&
3858  "of the surface region (mixed & buffer layer) thickness.", &
3859  units="nondim", default=0.5)
3860  call get_param(param_file, mdl, "LIMIT_BUFFER_DET_DH_BATHY", cs%lim_det_dH_bathy, &
3861  "The fraction of the total depth by which the thickness \n"//&
3862  "of the surface region (mixed & buffer layer) is allowed \n"//&
3863  "to change between grid points.", units="nondim", default=0.2)
3864  endif
3865 
3866  call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
3867  "If true, temperature and salinity are used as state \n"//&
3868  "variables.", default=.true.)
3869  cs%nsw = 0
3870  if (use_temperature) then
3871  call get_param(param_file, mdl, "PEN_SW_NBANDS", cs%nsw, default=1)
3872  endif
3873 
3874 
3875  if (max(cs%id_TKE_wind, cs%id_TKE_RiBulk, cs%id_TKE_conv, &
3876  cs%id_TKE_mixing, cs%id_TKE_pen_SW, cs%id_TKE_mech_decay, &
3877  cs%id_TKE_conv_decay) > 0) then
3878  call safe_alloc_alloc(cs%diag_TKE_wind, isd, ied, jsd, jed)
3879  call safe_alloc_alloc(cs%diag_TKE_RiBulk, isd, ied, jsd, jed)
3880  call safe_alloc_alloc(cs%diag_TKE_conv, isd, ied, jsd, jed)
3881  call safe_alloc_alloc(cs%diag_TKE_pen_SW, isd, ied, jsd, jed)
3882  call safe_alloc_alloc(cs%diag_TKE_mixing, isd, ied, jsd, jed)
3883  call safe_alloc_alloc(cs%diag_TKE_mech_decay, isd, ied, jsd, jed)
3884  call safe_alloc_alloc(cs%diag_TKE_conv_decay, isd, ied, jsd, jed)
3885  call safe_alloc_alloc(cs%diag_TKE_conv_s2, isd, ied, jsd, jed)
3886 
3887  cs%TKE_diagnostics = .true.
3888  endif
3889  if (cs%id_PE_detrain > 0) call safe_alloc_alloc(cs%diag_PE_detrain, isd, ied, jsd, jed)
3890  if (cs%id_PE_detrain2 > 0) call safe_alloc_alloc(cs%diag_PE_detrain2, isd, ied, jsd, jed)
3891  if (cs%id_ML_depth > 0) call safe_alloc_alloc(cs%ML_depth, isd, ied, jsd, jed)
3892 
3893  if(cs%allow_clocks_in_omp_loops) then
3894  id_clock_detrain = cpu_clock_id('(Ocean mixed layer detrain)', grain=clock_routine)
3895  id_clock_mech = cpu_clock_id('(Ocean mixed layer mechanical entrainment)', grain=clock_routine)
3896  id_clock_conv = cpu_clock_id('(Ocean mixed layer convection)', grain=clock_routine)
3897  if (cs%ML_resort) then
3898  id_clock_resort = cpu_clock_id('(Ocean mixed layer resorting)', grain=clock_routine)
3899  else
3900  id_clock_adjustment = cpu_clock_id('(Ocean mixed layer convective adjustment)', grain=clock_routine)
3901  endif
3902  id_clock_eos = cpu_clock_id('(Ocean mixed layer EOS)', grain=clock_routine)
3903  endif
3904 
3905  if (cs%limit_det .or. (cs%id_Hsfc_min > 0)) &
3906  id_clock_pass = cpu_clock_id('(Ocean mixed layer halo updates)', grain=clock_routine)
3907 
3908 
3909 ! if (CS%limit_det) then
3910 ! endif
3911 
Here is the caller graph for this function:

◆ convective_adjustment()

subroutine mom_bulk_mixed_layer::convective_adjustment ( real, dimension(szi_(g),szk_(gv)), intent(inout)  h,
real, dimension(szi_(g),szk_(gv)), intent(inout)  u,
real, dimension(szi_(g),szk_(gv)), intent(inout)  v,
real, dimension(szi_(g),szk_(gv)), intent(inout)  R0,
real, dimension(szi_(g),szk_(gv)), intent(inout)  Rcv,
real, dimension(szi_(g),szk_(gv)), intent(inout)  T,
real, dimension(szi_(g),szk_(gv)), intent(inout)  S,
real, dimension(szi_(g),szk_(gv)), intent(in)  eps,
real, dimension(szi_(g),szk_(gv)), intent(inout)  d_eb,
real, dimension(szi_(g),szk_(gv)), intent(out)  dKE_CA,
real, dimension(szi_(g),szk_(gv)), intent(out)  cTKE,
integer, intent(in)  j,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(bulkmixedlayer_cs), pointer  CS,
integer, intent(in), optional  nz_conv 
)
private

This subroutine does instantaneous convective entrainment into the buffer layers and mixed layers to remove hydrostatic instabilities. Any water that is lighter than currently in the mixed- or buffer- layer is entrained.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in,out]hLayer thickness, in m or kg m-2. (Intent in/out) The units of h are referred to as H below.
[in,out]uZonal velocities interpolated to h points, m s-1.
[in,out]vZonal velocities interpolated to h points, m s-1.
[in,out]tLayer temperatures, in deg C.
[in,out]r0Potential density referenced to surface pressure, in kg m-3.
[in,out]rcvThe coordinate defining potential density, in kg m-3.
[in,out]d_ebThe downward increase across a layer in the entrainment from below, in H. Positive values go with mass gain by a layer.
[out]dke_caThe vertically integrated change in kinetic energy due to convective adjustment, in m3 s-2.
[out]ctkeThe buoyant turbulent kinetic energy source due to convective adjustment, in m3 s-2.
[in]jThe j-index to work on.
csThe control structure for this module.
[in]nz_convIf present, the number of layers over which to do convective adjustment (perhaps CSnkml).

Definition at line 917 of file MOM_bulk_mixed_layer.F90.

Referenced by bulkmixedlayer().

917  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
918  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
919  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness, in m or kg m-2.
920  !! (Intent in/out) The units of h are
921  !! referred to as H below.
922  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: u !< Zonal velocities interpolated to h
923  !! points, m s-1.
924  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: v !< Zonal velocities interpolated to h
925  !! points, m s-1.
926  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: t !< Layer temperatures, in deg C.
927  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: s
928  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: r0 !< Potential density referenced to
929  !! surface pressure, in kg m-3.
930  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: rcv !< The coordinate defining potential
931  !! density, in kg m-3.
932  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_eb !< The downward increase across a layer
933  !! in the entrainment from below, in H.
934  !! Positive values go with mass gain by
935  !! a layer.
936  real, dimension(SZI_(G),SZK_(GV)), intent(in) :: eps
937  real, dimension(SZI_(G),SZK_(GV)), intent(out) :: dke_ca !< The vertically integrated change in
938  !! kinetic energy due to convective
939  !! adjustment, in m3 s-2.
940  real, dimension(SZI_(G),SZK_(GV)), intent(out) :: ctke !< The buoyant turbulent kinetic energy
941  !! source due to convective adjustment,
942  !! in m3 s-2.
943  integer, intent(in) :: j !< The j-index to work on.
944  type(bulkmixedlayer_cs), pointer :: cs !< The control structure for this module.
945  integer, optional, intent(in) :: nz_conv !< If present, the number of layers
946  !! over which to do convective adjustment
947  !! (perhaps CS%nkml).
948 
949 ! This subroutine does instantaneous convective entrainment into the buffer
950 ! layers and mixed layers to remove hydrostatic instabilities. Any water that
951 ! is lighter than currently in the mixed- or buffer- layer is entrained.
952 
953 ! Arguments: h - Layer thickness, in m or kg m-2. (Intent in/out) The units
954 ! of h are referred to as H below.
955 ! (in/out) u - Zonal velocities interpolated to h points, m s-1.
956 ! (in/out) v - Zonal velocities interpolated to h points, m s-1.
957 ! (in/out) R0 - Potential density referenced to surface pressure, in kg m-3.
958 ! (in/out) Rcv - The coordinate defining potential density, in kg m-3.
959 ! (in/out) T - Layer temperatures, in deg C.
960 ! (in/out) S - Layer salinities, in psu.
961 ! (in/out) d_eb - The downward increase across a layer in the entrainment from
962 ! below, in H. Positive values go with mass gain by a layer.
963 ! (out) dKE_CA - The vertically integrated change in kinetic energy due
964 ! to convective adjustment, in m3 s-2.
965 ! (out) cTKE - The buoyant turbulent kinetic energy source due to
966 ! convective adjustment, in m3 s-2.
967 ! (in) j - The j-index to work on.
968 ! (in) ksort - The density-sorted k-indicies.
969 ! (in) G - The ocean's grid structure.
970 ! (in) GV - The ocean's vertical grid structure.
971 ! (in) CS - The control structure for this module.
972 ! (in,opt) nz_conv - If present, the number of layers over which to do
973 ! convective adjustment (perhaps CS%nkml).
974  real, dimension(SZI_(G)) :: &
975  htot, & ! The total depth of the layers being considered for
976  ! entrainment, in H.
977  r0_tot, & ! The integrated potential density referenced to the surface
978  ! of the layers which are fully entrained, in H kg m-3.
979  rcv_tot, & ! The integrated coordinate value potential density of the
980  ! layers that are fully entrained, in H kg m-3.
981  ttot, & ! The integrated temperature of layers which are fully
982  ! entrained, in H K.
983  stot, & ! The integrated salt of layers which are fully entrained,
984  ! in H PSU.
985  uhtot, & ! The depth integrated zonal and meridional velocities in
986  vhtot, & ! the mixed layer, in H m s-1.
987  ke_orig, & ! The total mean kinetic energy in the mixed layer before
988  ! convection, H m2 s-2.
989  h_orig_k1 ! The depth of layer k1 before convective adjustment, in H.
990  real :: h_ent ! The thickness from a layer that is entrained, in H.
991  real :: ih ! The inverse of a thickness, in H-1.
992  real :: g_h2_2rho0 ! Half the gravitational acceleration times the
993  ! square of the conversion from H to m divided
994  ! by the mean density, in m6 s-2 H-2 kg-1.
995  integer :: is, ie, nz, i, k, k1, nzc, nkmb
996 
997  is = g%isc ; ie = g%iec ; nz = gv%ke
998  g_h2_2rho0 = (gv%g_Earth * gv%H_to_m**2) / (2.0 * gv%Rho0)
999  nzc = nz ; if (present(nz_conv)) nzc = nz_conv
1000  nkmb = cs%nkml+cs%nkbl
1001 
1002 ! Undergo instantaneous entrainment into the buffer layers and mixed layers
1003 ! to remove hydrostatic instabilities. Any water that is lighter than currently
1004 ! in the layer is entrained.
1005  do k1=min(nzc-1,nkmb),1,-1
1006  do i=is,ie
1007  h_orig_k1(i) = h(i,k1)
1008  ke_orig(i) = 0.5*h(i,k1)*(u(i,k1)**2 + v(i,k1)**2)
1009  uhtot(i) = h(i,k1)*u(i,k1) ; vhtot(i) = h(i,k1)*v(i,k1)
1010  r0_tot(i) = r0(i,k1) * h(i,k1)
1011  ctke(i,k1) = 0.0 ; dke_ca(i,k1) = 0.0
1012 
1013  rcv_tot(i) = rcv(i,k1) * h(i,k1)
1014  ttot(i) = t(i,k1) * h(i,k1) ; stot(i) = s(i,k1) * h(i,k1)
1015  enddo
1016  do k=k1+1,nzc
1017  do i=is,ie
1018  if ((h(i,k) > eps(i,k)) .and. (r0_tot(i) > h(i,k1)*r0(i,k))) then
1019  h_ent = h(i,k)-eps(i,k)
1020  ctke(i,k1) = ctke(i,k1) + (h_ent * g_h2_2rho0 * &
1021  (r0_tot(i) - h(i,k1)*r0(i,k)) * cs%nstar2)
1022  if (k < nkmb) then
1023  ctke(i,k1) = ctke(i,k1) + ctke(i,k)
1024  dke_ca(i,k1) = dke_ca(i,k1) + dke_ca(i,k)
1025  endif
1026  r0_tot(i) = r0_tot(i) + h_ent * r0(i,k)
1027  ke_orig(i) = ke_orig(i) + 0.5*h_ent* &
1028  (u(i,k)*u(i,k) + v(i,k)*v(i,k))
1029  uhtot(i) = uhtot(i) + h_ent*u(i,k)
1030  vhtot(i) = vhtot(i) + h_ent*v(i,k)
1031 
1032  rcv_tot(i) = rcv_tot(i) + h_ent * rcv(i,k)
1033  ttot(i) = ttot(i) + h_ent * t(i,k)
1034  stot(i) = stot(i) + h_ent * s(i,k)
1035  h(i,k1) = h(i,k1) + h_ent ; h(i,k) = eps(i,k)
1036 
1037  d_eb(i,k) = d_eb(i,k) - h_ent
1038  d_eb(i,k1) = d_eb(i,k1) + h_ent
1039  endif
1040  enddo
1041  enddo
1042 ! Determine the temperature, salinity, and velocities of the mixed or buffer
1043 ! layer in question, if it has entrained.
1044  do i=is,ie ; if (h(i,k1) > h_orig_k1(i)) then
1045  ih = 1.0 / h(i,k1)
1046  r0(i,k1) = r0_tot(i) * ih
1047  u(i,k1) = uhtot(i) * ih ; v(i,k1) = vhtot(i) * ih
1048  dke_ca(i,k1) = dke_ca(i,k1) + gv%H_to_m * (cs%bulk_Ri_convective * &
1049  (ke_orig(i) - 0.5*h(i,k1)*(u(i,k1)**2 + v(i,k1)**2)))
1050  rcv(i,k1) = rcv_tot(i) * ih
1051  t(i,k1) = ttot(i) * ih ; s(i,k1) = stot(i) * ih
1052  endif ; enddo
1053  enddo
1054 ! If lower mixed or buffer layers are massless, give them the properties of the
1055 ! layer above.
1056  do k=2,min(nzc,nkmb) ; do i=is,ie ; if (h(i,k) == 0.0) then
1057  r0(i,k) = r0(i,k-1)
1058  rcv(i,k) = rcv(i,k-1) ; t(i,k) = t(i,k-1) ; s(i,k) = s(i,k-1)
1059  endif ; enddo ; enddo
1060 
Here is the caller graph for this function:

◆ ef4()

real function mom_bulk_mixed_layer::ef4 ( real, intent(in)  H,
real, intent(in)  E,
real, intent(in)  L,
real, intent(inout), optional  dR_de 
)
private

This subroutine returns an approximation to the integral R = exp(-L*(H+E)) integral(LH to L(H+E)) L/(1-(1+x)exp(-x)) dx. The approximation to the integrand is good to within -2% at x~.3 and +25% at x~3.5, but the exponential deemphasizes the importance of large x. When L=0, EF4 returns E/((H+E)*H).

Parameters
[in]hTotal thickness, in m or kg m-2. (Intent in) The units of h are referred to as H below.
[in]eEntrainment, in units of H.
[in]lThe e-folding scale in H-1.
[in,out]dr_deThe partial derivative of the result R with E, in H-2.

Definition at line 3920 of file MOM_bulk_mixed_layer.F90.

Referenced by mechanical_entrainment().

3920 real, intent(in) :: h !< Total thickness, in m or kg m-2. (Intent in) The units of h
3921  !! are referred to as H below.
3922 real, intent(in) :: e !< Entrainment, in units of H.
3923 real, intent(in) :: l !< The e-folding scale in H-1.
3924 real, intent(inout), optional :: dr_de !< The partial derivative of the result R with E, in H-2.
3925 real :: ef4
3926 ! This subroutine returns an approximation to the integral
3927 ! R = exp(-L*(H+E)) integral(LH to L(H+E)) L/(1-(1+x)exp(-x)) dx.
3928 ! The approximation to the integrand is good to within -2% at x~.3
3929 ! and +25% at x~3.5, but the exponential deemphasizes the importance of
3930 ! large x. When L=0, EF4 returns E/((H+E)*H).
3931 !
3932 ! Arguments: h - Total thickness, in m or kg m-2. (Intent in) The units
3933 ! of h are referred to as H below.
3934 ! (in) E - Entrainment, in units of H.
3935 ! (in) L - The e-folding scale in H-1.
3936 ! (out) dR_de - the partial derivative of the result R with E, in H-2.
3937 ! (return value) R - The integral, in units of H-1.
3938  real :: exp_lhpe ! A nondimensional exponential decay.
3939  real :: i_hpe ! An inverse thickness plus entrainment, in H-1.
3940  real :: r ! The result of the integral above, in H-1.
3941 
3942  exp_lhpe = exp(-l*(e+h))
3943  i_hpe = 1.0/(h+e)
3944  r = exp_lhpe * (e*i_hpe/h - 0.5*l*log(h*i_hpe) + 0.5*l*l*e)
3945  if (PRESENT(dr_de)) &
3946  dr_de = -l*r + exp_lhpe*(i_hpe*i_hpe + 0.5*l*i_hpe + 0.5*l*l)
3947  ef4 = r
3948 
Here is the caller graph for this function:

◆ find_starting_tke()

subroutine mom_bulk_mixed_layer::find_starting_tke ( real, dimension(szi_(g)), intent(in)  htot,
real, dimension(szi_(g)), intent(in)  h_CA,
type(forcing), intent(in)  fluxes,
real, dimension(szi_(g)), intent(inout)  Conv_En,
real, dimension(szi_(g),szk_(gv)), intent(in)  cTKE,
real, dimension(szi_(g)), intent(in)  dKE_FC,
real, dimension(szi_(g),szk_(gv)), intent(in)  dKE_CA,
real, dimension(szi_(g)), intent(out)  TKE,
real, dimension(szi_(g)), intent(in)  TKE_river,
real, dimension(szi_(g)), intent(out)  Idecay_len_TKE,
real, dimension(2,szi_(g)), intent(out)  cMKE,
real, intent(in)  dt,
real, intent(in)  Idt_diag,
integer, intent(in)  j,
integer, dimension(szi_(g),szk_(gv)), intent(in)  ksort,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(bulkmixedlayer_cs), pointer  CS 
)
private

This subroutine determines the TKE available at the depth of free convection to drive mechanical entrainment.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]htotThe accumlated mixed layer thickness, in m or kg m-2. (Intent in).
[in]h_caThe mixed layer depth after convective adjustment, in H.
[in]fluxesA structure containing pointers to any possible forcing fields. Unused fields have NULL ptrs.
[in,out]conv_enThe buoyant turbulent kinetic energy source due to free convection, in m3 s-2.
[in]dke_fcThe vertically integrated change in kinetic energy due to free convection, in m3 s-2.
[in]ctkeThe buoyant turbulent kinetic energy
[in]dke_caThe vertically integrated change in
[out]tkeThe turbulent kinetic energy available for mixing over a time step, in m3 s-2.
[out]idecay_len_tkeThe inverse of the vertical decay scale for TKE, in H-1.
[out]cmkeCoefficients of HpE and HpE^2 in calculating the denominator of MKE_rate, in H-1 and H-2.
[in]dtThe time step in s.
[in]jThe j-index to work on.
[in]ksortThe density-sorted k-indicies.
csThe control structure for this module.

Definition at line 1431 of file MOM_bulk_mixed_layer.F90.

Referenced by bulkmixedlayer().

1431  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
1432  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
1433  real, dimension(SZI_(G)), intent(in) :: htot !< The accumlated mixed layer thickness, in m
1434  !! or kg m-2. (Intent in).
1435  real, dimension(SZI_(G)), intent(in) :: h_ca !< The mixed layer depth after convective
1436  !! adjustment, in H.
1437  type(forcing), intent(in) :: fluxes !< A structure containing pointers to any
1438  !! possible forcing fields. Unused fields
1439  !! have NULL ptrs.
1440  real, dimension(SZI_(G)), intent(inout) :: conv_en !< The buoyant turbulent kinetic energy
1441  !! source due to free convection,
1442  !! in m3 s-2.
1443  real, dimension(SZI_(G)), intent(in) :: dke_fc !< The vertically integrated change in
1444  !! kinetic energy due to free convection,
1445  !! in m3 s-2.
1446  real, dimension(SZI_(G),SZK_(GV)), &
1447  intent(in) :: ctke !< The buoyant turbulent kinetic energy
1448  !! source due to convective adjustment,
1449  !! in m3 s-2.
1450  real, dimension(SZI_(G),SZK_(GV)), &
1451  intent(in) :: dke_ca !< The vertically integrated change in
1452  !! kinetic energy due to convective
1453  !! adjustment, in m3 s-2.
1454  real, dimension(SZI_(G)), intent(out) :: tke !< The turbulent kinetic energy available for
1455  !! mixing over a time step, in m3 s-2.
1456  real, dimension(SZI_(G)), intent(out) :: idecay_len_tke !< The inverse of the vertical decay
1457  !! scale for TKE, in H-1.
1458  real, dimension(SZI_(G)), intent(in) :: tke_river
1459  real, dimension(2,SZI_(G)), intent(out) :: cmke !< Coefficients of HpE and HpE^2 in
1460  !! calculating the denominator of MKE_rate,
1461  !! in H-1 and H-2.
1462  real, intent(in) :: dt !< The time step in s.
1463  real, intent(in) :: idt_diag
1464  integer, intent(in) :: j !< The j-index to work on.
1465  integer, dimension(SZI_(G),SZK_(GV)), &
1466  intent(in) :: ksort !< The density-sorted k-indicies.
1467  type(bulkmixedlayer_cs), pointer :: cs !< The control structure for this module.
1468 
1469 ! This subroutine determines the TKE available at the depth of free
1470 ! convection to drive mechanical entrainment.
1471 
1472 ! Arguments: htot - The accumlated mixed layer thickness, in m or kg m-2. (Intent in)
1473 ! The units of htot are referred to as H below.
1474 ! (in) h_CA - The mixed layer depth after convective adjustment, in H.
1475 ! (in) fluxes - A structure containing pointers to any possible
1476 ! forcing fields. Unused fields have NULL ptrs.
1477 ! (in) Conv_en - The buoyant turbulent kinetic energy source due to
1478 ! free convection, in m3 s-2.
1479 ! (in) cTKE - The buoyant turbulent kinetic energy source due to
1480 ! convective adjustment, in m3 s-2.
1481 ! (in) dKE_FC - The vertically integrated change in kinetic energy due
1482 ! to free convection, in m3 s-2.
1483 ! (in) dKE_CA - The vertically integrated change in kinetic energy due
1484 ! to convective adjustment, in m3 s-2.
1485 ! (out) TKE - The turbulent kinetic energy available for mixing over a
1486 ! time step, in m3 s-2.
1487 ! (out) Idecay_len_TKE - The inverse of the vertical decay scale for
1488 ! TKE, in H-1.
1489 ! (out) cMKE - Coefficients of HpE and HpE^2 in calculating the
1490 ! denominator of MKE_rate, in H-1 and H-2.
1491 ! (in) dt - The time step in s.
1492 ! (in) j - The j-index to work on.
1493 ! (in) ksort - The density-sorted k-indicies.
1494 ! (in) G - The ocean's grid structure.
1495 ! (in) GV - The ocean's vertical grid structure.
1496 ! (in) CS - The control structure for this module.
1497  real :: dke_conv ! The change in mean kinetic energy due
1498  ! to all convection, in m3 s-2.
1499  real :: nstar_fc ! The effective efficiency with which the energy released by
1500  ! free convection is converted to TKE, often ~0.2, ND.
1501  real :: nstar_ca ! The effective efficiency with which the energy released by
1502  ! convective adjustment is converted to TKE, often ~0.2, ND.
1503  real :: tke_ca ! The potential energy released by convective adjustment if
1504  ! that release is positive, in m3 s2.
1505  real :: mke_rate_ca ! MKE_rate for convective adjustment, ND, 0 to 1.
1506  real :: mke_rate_fc ! MKE_rate for free convection, ND, 0 to 1.
1507  real :: toten ! The total potential energy released by convection, m3 s-2.
1508  real :: ih ! The inverse of a thickness, in H-1.
1509  real :: exp_kh ! The nondimensional decay of TKE across a layer, ND.
1510  real :: absf ! The absolute value of f averaged to thickness points, s-1.
1511  real :: u_star ! The friction velocity in m s-1.
1512  real :: absf_ustar ! The absolute value of f divided by U_star, in m-1.
1513  real :: wind_tke_src ! The surface wind source of TKE, in m3 s-3.
1514  real :: diag_wt ! The ratio of the current timestep to the diagnostic
1515  ! timestep (which may include 2 calls), ND.
1516  integer :: is, ie, nz, i
1517 
1518  is = g%isc ; ie = g%iec ; nz = gv%ke
1519  diag_wt = dt * idt_diag
1520 
1521  if (cs%omega_frac >= 1.0) absf = 2.0*cs%omega
1522  do i=is,ie
1523  u_star = fluxes%ustar(i,j)
1524  if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then
1525  if (fluxes%frac_shelf_h(i,j) > 0.0) &
1526  u_star = (1.0 - fluxes%frac_shelf_h(i,j)) * u_star + &
1527  fluxes%frac_shelf_h(i,j) * fluxes%ustar_shelf(i,j)
1528  endif
1529 
1530  if (u_star < cs%ustar_min) u_star = cs%ustar_min
1531  if (cs%omega_frac < 1.0) then
1532  absf = 0.25*((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
1533  (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i-1,j))))
1534  if (cs%omega_frac > 0.0) &
1535  absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
1536  endif
1537  absf_ustar = absf / u_star
1538  idecay_len_tke(i) = (absf_ustar * cs%TKE_decay) * gv%H_to_m
1539 
1540 ! The first number in the denominator could be anywhere up to 16. The
1541 ! value of 3 was chosen to minimize the time-step dependence of the amount
1542 ! of shear-driven mixing in 10 days of a 1-degree global model, emphasizing
1543 ! the equatorial areas. Although it is not cast as a parameter, it should
1544 ! be considered an empirical parameter, and it might depend strongly on the
1545 ! number of sublayers in the mixed layer and their locations.
1546 ! The 0.41 is VonKarman's constant. This equation assumes that small & large
1547 ! scales contribute to mixed layer deepening at similar rates, even though
1548 ! small scales are dissipated more rapidly (implying they are less efficient).
1549 ! Ih = 1.0/(16.0*0.41*U_star*dt)
1550  ih = gv%H_to_m/(3.0*0.41*u_star*dt)
1551  cmke(1,i) = 4.0 * ih ; cmke(2,i) = (absf_ustar*gv%H_to_m) * ih
1552 
1553  if (idecay_len_tke(i) > 0.0) then
1554  exp_kh = exp(-htot(i)*idecay_len_tke(i))
1555  else
1556  exp_kh = 1.0
1557  endif
1558 
1559 ! Here nstar is a function of the natural Rossby number 0.2/(1+0.2/Ro), based
1560 ! on a curve fit from the data of Wang (GRL, 2003).
1561 ! Note: Ro = 1.0/sqrt(0.5 * dt * (absf*htot(i))**3 / totEn)
1562  if (conv_en(i) < 0.0) conv_en(i) = 0.0
1563  if (ctke(i,1) > 0.0) then ; tke_ca = ctke(i,1) ; else ; tke_ca = 0.0 ; endif
1564  if ((htot(i) >= h_ca(i)) .or. (tke_ca == 0.0)) then
1565  toten = conv_en(i) + tke_ca
1566 
1567  if (toten > 0.0) then
1568  nstar_fc = cs%nstar * toten / (toten + 0.2 * &
1569  sqrt(0.5 * dt * (absf*(htot(i)*gv%H_to_m))**3 * toten))
1570  else
1571  nstar_fc = cs%nstar
1572  endif
1573  nstar_ca = nstar_fc
1574  else
1575  ! This reconstructs the Buoyancy flux within the topmost htot of water.
1576  if (conv_en(i) > 0.0) then
1577  toten = conv_en(i) + tke_ca * (htot(i) / h_ca(i))
1578  nstar_fc = cs%nstar * toten / (toten + 0.2 * &
1579  sqrt(0.5 * dt * (absf*(htot(i)*gv%H_to_m))**3 * toten))
1580  else
1581  nstar_fc = cs%nstar
1582  endif
1583 
1584  toten = conv_en(i) + tke_ca
1585  if (tke_ca > 0.0) then
1586  nstar_ca = cs%nstar * toten / (toten + 0.2 * &
1587  sqrt(0.5 * dt * (absf*(h_ca(i)*gv%H_to_m))**3 * toten))
1588  else
1589  nstar_ca = cs%nstar
1590  endif
1591  endif
1592 
1593  if (dke_fc(i) + dke_ca(i,1) > 0.0) then
1594  if (htot(i) >= h_ca(i)) then
1595  mke_rate_fc = 1.0 / (1.0 + htot(i)*(cmke(1,i) + cmke(2,i)*htot(i)) )
1596  mke_rate_ca = mke_rate_fc
1597  else
1598  mke_rate_fc = 1.0 / (1.0 + htot(i)*(cmke(1,i) + cmke(2,i)*htot(i)) )
1599  mke_rate_ca = 1.0 / (1.0 + h_ca(i)*(cmke(1,i) + cmke(2,i)*h_ca(i)) )
1600  endif
1601  else
1602  ! This branch just saves unnecessary calculations.
1603  mke_rate_fc = 1.0 ; mke_rate_ca = 1.0
1604  endif
1605 
1606  dke_conv = dke_ca(i,1) * mke_rate_ca + dke_fc(i) * mke_rate_fc
1607 ! At this point, it is assumed that cTKE is positive and stored in TKE_CA!
1608 ! Note: Removed factor of 2 in u*^3 terms.
1609  tke(i) = (dt*cs%mstar)*((u_star*u_star*u_star)*exp_kh) + &
1610  (exp_kh * dke_conv + nstar_fc*conv_en(i) + nstar_ca*tke_ca)
1611 ! Add additional TKE at river mouths
1612 
1613  if (cs%do_rivermix) then
1614  tke(i) = tke(i) + tke_river(i)*dt*exp_kh
1615  endif
1616 
1617  if (cs%TKE_diagnostics) then
1618  wind_tke_src = cs%mstar*(u_star*u_star*u_star) * diag_wt
1619  cs%diag_TKE_wind(i,j) = cs%diag_TKE_wind(i,j) + &
1620  wind_tke_src + tke_river(i) * diag_wt
1621  cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + dke_conv*idt_diag
1622  cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
1623  (exp_kh-1.0)*(wind_tke_src + dke_conv*idt_diag)
1624  cs%diag_TKE_conv(i,j) = cs%diag_TKE_conv(i,j) + &
1625  idt_diag*(nstar_fc*conv_en(i) + nstar_ca*tke_ca)
1626  cs%diag_TKE_conv_decay(i,j) = cs%diag_TKE_conv_decay(i,j) + &
1627  idt_diag*((cs%nstar-nstar_fc)*conv_en(i) + (cs%nstar-nstar_ca)*tke_ca)
1628  cs%diag_TKE_conv_s2(i,j) = cs%diag_TKE_conv_s2(i,j) + &
1629  idt_diag*(ctke(i,1)-tke_ca)
1630  endif
1631  enddo
1632 
Here is the caller graph for this function:

◆ mechanical_entrainment()

subroutine mom_bulk_mixed_layer::mechanical_entrainment ( real, dimension(szi_(g),szk_(gv)), intent(inout)  h,
real, dimension(szi_(g),szk_(gv)), intent(inout)  d_eb,
real, dimension(szi_(g)), intent(inout)  htot,
real, dimension(szi_(g)), intent(inout)  Ttot,
real, dimension(szi_(g)), intent(inout)  Stot,
real, dimension(szi_(g)), intent(inout)  uhtot,
real, dimension(szi_(g)), intent(inout)  vhtot,
real, dimension(szi_(g)), intent(inout)  R0_tot,
real, dimension(szi_(g)), intent(inout)  Rcv_tot,
real, dimension(szi_(g),szk_(gv)), intent(in)  u,
real, dimension(szi_(g),szk_(gv)), intent(in)  v,
real, dimension(szi_(g),szk_(gv)), intent(in)  T,
real, dimension(szi_(g),szk_(gv)), intent(in)  S,
real, dimension(szi_(g),szk_(gv)), intent(in)  R0,
real, dimension(szi_(g),szk_(gv)), intent(in)  Rcv,
real, dimension(szi_(g),szk_(gv)), intent(in)  eps,
real, dimension(szi_(g)), intent(in)  dR0_dT,
real, dimension(szi_(g)), intent(in)  dRcv_dT,
real, dimension(2,szi_(g)), intent(in)  cMKE,
real, intent(in)  Idt_diag,
integer, intent(in)  nsw,
real, dimension(:,:), intent(inout)  Pen_SW_bnd,
real, dimension(:,:,:), intent(in)  opacity_band,
real, dimension(szi_(g)), intent(inout)  TKE,
real, dimension(szi_(g)), intent(inout)  Idecay_len_TKE,
integer, intent(in)  j,
integer, dimension(szi_(g),szk_(gv)), intent(in)  ksort,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(bulkmixedlayer_cs), pointer  CS 
)
private

This subroutine calculates mechanically driven entrainment.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in,out]hLayer thickness, in m or kg m-2. (Intent in/out) The units of h are referred to as H below.
[in,out]d_ebThe downward increase across a layer in the entrainment from below, in H. Positive values go with mass gain by a layer.
[in,out]htotThe accumlated mixed layer thickness, in H.
[in,out]ttotThe depth integrated mixed layer temperature, in deg C H.
[in,out]stotThe depth integrated mixed layer salinity, in psu H.
[in,out]uhtotThe depth integrated mixed layer zonal velocity, H m s-1.
[in,out]vhtotThe integrated mixed layer meridional velocity, H m s-1.
[in,out]r0_totThe integrated mixed layer potential density referenced to 0 pressure, in H kg m-3.
[in,out]rcv_totThe integrated mixed layer coordinate variable potential density, in H kg m-3.
[in]nswThe number of bands of penetrating shortwave radiation.
[in,out]pen_sw_bndThe penetrating shortwave heating at the sea surface in each penetrating band, in K H, size nsw x SZI_(G).
[in,out]tkeThe turbulent kinetic energy available for mixing over a time step, in m3 s-2.
[in]jThe j-index to work on.
[in]ksortThe density-sorted k-indicies.
csThe control structure for this module.

Definition at line 1641 of file MOM_bulk_mixed_layer.F90.

References ef4().

Referenced by bulkmixedlayer().

1641  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
1642  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid
1643  !! structure.
1644  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness, in m or kg m-2.
1645  !! (Intent in/out) The units of h are
1646  !! referred to as H below.
1647  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_eb !< The downward increase across a
1648  !! layer in the entrainment from
1649  !! below, in H. Positive values go
1650  !! with mass gain by a layer.
1651  real, dimension(SZI_(G)), intent(inout) :: htot !< The accumlated mixed layer
1652  !! thickness, in H.
1653  real, dimension(SZI_(G)), intent(inout) :: ttot !< The depth integrated mixed layer
1654  !! temperature, in deg C H.
1655  real, dimension(SZI_(G)), intent(inout) :: stot !< The depth integrated mixed layer
1656  !! salinity, in psu H.
1657  real, dimension(SZI_(G)), intent(inout) :: uhtot !< The depth integrated mixed layer
1658  !! zonal velocity, H m s-1.
1659  real, dimension(SZI_(G)), intent(inout) :: vhtot !< The integrated mixed layer
1660  !! meridional velocity, H m s-1.
1661  real, dimension(SZI_(G)), intent(inout) :: r0_tot !< The integrated mixed layer
1662  !! potential density referenced to 0
1663  !! pressure, in H kg m-3.
1664  real, dimension(SZI_(G)), intent(inout) :: rcv_tot !< The integrated mixed layer
1665  !! coordinate variable potential
1666  !! density, in H kg m-3.
1667  real, dimension(SZI_(G),SZK_(GV)), intent(in) :: u, v, t, s, r0, rcv, eps
1668  real, dimension(SZI_(G)), intent(in) :: dr0_dt, drcv_dt
1669  real, dimension(2,SZI_(G)), intent(in) :: cmke
1670  real, intent(in) :: idt_diag
1671  integer, intent(in) :: nsw !< The number of bands of penetrating
1672  !! shortwave radiation.
1673  real, dimension(:,:), intent(inout) :: pen_sw_bnd !< The penetrating shortwave
1674  !! heating at the sea surface in each
1675  !! penetrating band, in K H,
1676  !! size nsw x SZI_(G).
1677  real, dimension(:,:,:), intent(in) :: opacity_band
1678  real, dimension(SZI_(G)), intent(inout) :: tke !< The turbulent kinetic energy
1679  !! available for mixing over a time
1680  !! step, in m3 s-2.
1681  real, dimension(SZI_(G)), intent(inout) :: idecay_len_tke
1682  integer, intent(in) :: j !< The j-index to work on.
1683  integer, dimension(SZI_(G),SZK_(GV)), &
1684  intent(in) :: ksort !< The density-sorted k-indicies.
1685  type(bulkmixedlayer_cs), pointer :: cs !< The control structure for this
1686  !! module.
1687 
1688 ! This subroutine calculates mechanically driven entrainment.
1689 ! Arguments: h - Layer thickness, in m or kg m-2. (Intent in/out) The units
1690 ! of h are referred to as H below.
1691 ! (in/out) d_eb - The downward increase across a layer in the entrainment from
1692 ! below, in H. Positive values go with mass gain by a layer.
1693 ! (in/out) htot - The accumlated mixed layer thickness, in H.
1694 ! (in/out) Ttot - The depth integrated mixed layer temperature, in deg C H.
1695 ! (in/out) Stot - The depth integrated mixed layer salinity, in psu H.
1696 ! (in/out) uhtot - The depth integrated mixed layer zonal velocity, H m s-1.
1697 ! (in/out) vhtot - The integrated mixed layer meridional velocity, H m s-1.
1698 ! (in/out) R0_tot - The integrated mixed layer potential density referenced
1699 ! to 0 pressure, in H kg m-3.
1700 ! (in/out) Rcv_tot - The integrated mixed layer coordinate variable
1701 ! potential density, in H kg m-3.
1702 ! (in) nsw - The number of bands of penetrating shortwave radiation.
1703 ! (in/out) Pen_SW_bnd - The penetrating shortwave heating at the sea surface
1704 ! in each penetrating band, in K H, size nsw x SZI_(G).
1705 ! (in/out) TKE - The turbulent kinetic energy available for mixing over a
1706 ! time step, in m3 s-2.
1707 ! (in) j - The j-index to work on.
1708 ! (in) ksort - The density-sorted k-indicies.
1709 ! (in) G - The ocean's grid structure.
1710 ! (in) GV - The ocean's vertical grid structure.
1711 ! (in) CS - The control structure for this module.
1712 
1713  real :: sw_trans ! The fraction of shortwave radiation that is not
1714  ! absorbed in a layer, nondimensional.
1715  real :: pen_absorbed ! The amount of penetrative shortwave radiation
1716  ! that is absorbed in a layer, in units of K m.
1717  real :: h_avail ! The thickness in a layer available for entrainment in H.
1718  real :: h_ent ! The thickness from a layer that is entrained, in H.
1719  real :: h_min, h_max ! Limits on the solution for h_ent, in H.
1720  real :: dh_newt ! The Newton's method estimate of the change in
1721  ! h_ent between iterations, in H.
1722  real :: mke_rate ! The fraction of the energy in resolved shears
1723  ! within the mixed layer that will be eliminated
1724  ! within a timestep, nondim, 0 to 1.
1725  real :: hpe ! The current thickness plus entrainment, H.
1726  real :: g_h_2rho0 ! Half the gravitational acceleration times the
1727  ! conversion from H to m divided by the mean density,
1728  ! in m5 s-2 H-1 kg-1.
1729  real :: tke_full_ent ! The TKE remaining if a layer is fully entrained,
1730  ! in units of m3 s-2.
1731  real :: drl ! Work required to mix water from the next layer
1732  ! across the mixed layer, in m2 s-2.
1733  real :: pen_en_contrib ! Penetrating SW contributions to the changes in
1734  ! TKE, divided by layer thickness in m, in m2 s-2.
1735  real :: c1 ! A temporary variable in units of m2 s-2.
1736  real :: dmke ! A temporary variable related to the release of mean
1737  ! kinetic energy, with units of H m3 s-2.
1738  real :: tke_ent ! The TKE that remains if h_ent were entrained, in m3 s-2.
1739  real :: tke_ent1 ! The TKE that would remain, without considering the
1740  ! release of mean kinetic energy, in m3 s2.
1741  real :: dtke_dh ! The partial derivative of TKE with h_ent, in m3 s-2 H-1.
1742  real :: pen_dtke_dh_contrib ! The penetrating shortwave contribution to
1743  ! dTKE_dh, in m2 s-2.
1744  real :: ef4_val ! The result of EF4() (see later), in H-1.
1745  real :: h_neglect ! A thickness that is so small it is usually lost
1746  ! in roundoff and can be neglected, in H.
1747  real :: def4_dh ! The partial derivative of EF4 with h, in H-2.
1748  real :: pen_en1 ! A nondimensional temporary variable.
1749  real :: kh, exp_kh ! Nondimensional temporary variables related to the.
1750  real :: f1_kh ! fractional decay of TKE across a layer.
1751  real :: x1, e_x1 ! Nondimensional temporary variables related to
1752  real :: f1_x1, f2_x1 ! the relative decay of TKE and SW radiation across
1753  real :: f3_x1 ! a layer, and exponential-related functions of x1.
1754  real :: e_hxhpe ! Entrainment divided by the product of the new and old
1755  ! thicknesses, in H-1.
1756  real :: hmix_min ! The minimum mixed layer depth in H.
1757  real :: h_to_m ! Local copies of unit conversion factors.
1758  real :: opacity
1759  real :: c1_3, c1_6, c1_24 ! 1/3, 1/6, and 1/24.
1760  integer :: is, ie, nz, i, k, ks, itt, n
1761 
1762  c1_3 = 1.0/3.0 ; c1_6 = 1.0/6.0 ; c1_24 = 1.0/24.0
1763  h_to_m = gv%H_to_m
1764  g_h_2rho0 = (gv%g_Earth * h_to_m) / (2.0 * gv%Rho0)
1765  hmix_min = cs%Hmix_min * gv%m_to_H
1766  h_neglect = gv%H_subroundoff
1767  is = g%isc ; ie = g%iec ; nz = gv%ke
1768 
1769  do ks=1,nz
1770 
1771  do i=is,ie ; if (ksort(i,ks) > 0) then
1772  k = ksort(i,ks)
1773 
1774  h_avail = h(i,k) - eps(i,k)
1775  if ((h_avail > 0.) .and. ((tke(i) > 0.) .or. (htot(i) < hmix_min))) then
1776  drl = g_h_2rho0 * (r0(i,k)*htot(i) - r0_tot(i) )
1777  dmke = (h_to_m * cs%bulk_Ri_ML) * 0.5 * &
1778  ((uhtot(i)-u(i,k)*htot(i))**2 + (vhtot(i)-v(i,k)*htot(i))**2)
1779 
1780 ! Find the TKE that would remain if the entire layer were entrained.
1781  kh = idecay_len_tke(i)*h_avail ; exp_kh = exp(-kh)
1782  if (kh >= 2.0e-5) then ; f1_kh = (1.0-exp_kh) / kh
1783  else ; f1_kh = (1.0 - kh*(0.5 - c1_6*kh)) ; endif
1784 
1785  pen_en_contrib = 0.0
1786  do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1787  opacity = opacity_band(n,i,k)
1788 ! Two different forms are used here to make sure that only negative
1789 ! values are taken into exponentials to avoid excessively large
1790 ! numbers. They are, of course, mathematically identical.
1791  if (idecay_len_tke(i) > opacity) then
1792  x1 = (idecay_len_tke(i) - opacity) * h_avail
1793  if (x1 >= 2.0e-5) then
1794  e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1795  f3_x1 = ((e_x1-(1.0-x1))/(x1*x1))
1796  else
1797  f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1798  f3_x1 = (0.5 - x1*(c1_6 - c1_24*x1))
1799  endif
1800 
1801  pen_en1 = exp(-opacity*h_avail) * &
1802  ((1.0+opacity*htot(i))*f1_x1 + opacity*h_avail*f3_x1)
1803  else
1804  x1 = (opacity - idecay_len_tke(i)) * h_avail
1805  if (x1 >= 2.0e-5) then
1806  e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1807  f2_x1 = ((1.0-(1.0+x1)*e_x1)/(x1*x1))
1808  else
1809  f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1810  f2_x1 = (0.5 - x1*(c1_3 - 0.125*x1))
1811  endif
1812 
1813  pen_en1 = exp_kh * ((1.0+opacity*htot(i))*f1_x1 + &
1814  opacity*h_avail*f2_x1)
1815  endif
1816  pen_en_contrib = pen_en_contrib + &
1817  (g_h_2rho0*dr0_dt(i)*pen_sw_bnd(n,i)) * (pen_en1 - f1_kh)
1818  endif ; enddo
1819 
1820  hpe = htot(i)+h_avail
1821  mke_rate = 1.0/(1.0 + (cmke(1,i)*hpe + cmke(2,i)*hpe**2))
1822  ef4_val = ef4(htot(i)+h_neglect,h_avail,idecay_len_tke(i))
1823  tke_full_ent = (exp_kh*tke(i) - (h_avail*h_to_m)*(drl*f1_kh + pen_en_contrib)) + &
1824  mke_rate*dmke*ef4_val
1825  if ((tke_full_ent >= 0.0) .or. (h_avail+htot(i) <= hmix_min)) then
1826  ! The layer will be fully entrained.
1827  h_ent = h_avail
1828 
1829  if (cs%TKE_diagnostics) then
1830  e_hxhpe = h_ent / ((htot(i)+h_neglect)*(htot(i)+h_ent+h_neglect))
1831  cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
1832  idt_diag * ((exp_kh-1.0)*tke(i) + &
1833  (h_ent*h_to_m)*drl*(1.0-f1_kh) + &
1834  mke_rate*dmke*(ef4_val-e_hxhpe))
1835  cs%diag_TKE_mixing(i,j) = cs%diag_TKE_mixing(i,j) - &
1836  idt_diag*(h_to_m*h_ent)*drl
1837  cs%diag_TKE_pen_SW(i,j) = cs%diag_TKE_pen_SW(i,j) - &
1838  idt_diag*(h_to_m*h_ent)*pen_en_contrib
1839  cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + &
1840  idt_diag*mke_rate*dmke*e_hxhpe
1841  endif
1842 
1843  tke(i) = tke_full_ent
1844  if (tke(i) <= 0.0) tke(i) = 1e-300
1845  else
1846 ! The layer is only partially entrained. The amount that will be
1847 ! entrained is determined iteratively. No further layers will be
1848 ! entrained.
1849  h_min = 0.0 ; h_max = h_avail
1850  if (tke(i) <= 0.0) then
1851  h_ent = 0.0
1852  else
1853  h_ent = h_avail * tke(i) / (tke(i) - tke_full_ent)
1854 
1855  do itt=1,15
1856  ! Evaluate the TKE that would remain if h_ent were entrained.
1857 
1858  kh = idecay_len_tke(i)*h_ent ; exp_kh = exp(-kh)
1859  if (kh >= 2.0e-5) then
1860  f1_kh = (1.0-exp_kh) / kh
1861  else
1862  f1_kh = (1.0 - kh*(0.5 - c1_6*kh))
1863  endif
1864 
1865 
1866  pen_en_contrib = 0.0 ; pen_dtke_dh_contrib = 0.0
1867  do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1868  ! Two different forms are used here to make sure that only negative
1869  ! values are taken into exponentials to avoid excessively large
1870  ! numbers. They are, of course, mathematically identical.
1871  opacity = opacity_band(n,i,k)
1872  sw_trans = exp(-h_ent*opacity)
1873  if (idecay_len_tke(i) > opacity) then
1874  x1 = (idecay_len_tke(i) - opacity) * h_ent
1875  if (x1 >= 2.0e-5) then
1876  e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1877  f3_x1 = ((e_x1-(1.0-x1))/(x1*x1))
1878  else
1879  f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1880  f3_x1 = (0.5 - x1*(c1_6 - c1_24*x1))
1881  endif
1882  pen_en1 = sw_trans * ((1.0+opacity*htot(i))*f1_x1 + &
1883  opacity*h_ent*f3_x1)
1884  else
1885  x1 = (opacity - idecay_len_tke(i)) * h_ent
1886  if (x1 >= 2.0e-5) then
1887  e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1888  f2_x1 = ((1.0-(1.0+x1)*e_x1)/(x1*x1))
1889  else
1890  f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1891  f2_x1 = (0.5 - x1*(c1_3 - 0.125*x1))
1892  endif
1893 
1894  pen_en1 = exp_kh * ((1.0+opacity*htot(i))*f1_x1 + &
1895  opacity*h_ent*f2_x1)
1896  endif
1897  c1 = g_h_2rho0*dr0_dt(i)*pen_sw_bnd(n,i)
1898  pen_en_contrib = pen_en_contrib + c1*(pen_en1 - f1_kh)
1899  pen_dtke_dh_contrib = pen_dtke_dh_contrib + &
1900  c1*((1.0-sw_trans) - opacity*(htot(i) + h_ent)*sw_trans)
1901  endif ; enddo ! (Pen_SW_bnd(n,i) > 0.0)
1902 
1903  tke_ent1 = exp_kh*tke(i) - (h_ent*h_to_m)*(drl*f1_kh + pen_en_contrib)
1904  ef4_val = ef4(htot(i)+h_neglect,h_ent,idecay_len_tke(i),def4_dh)
1905  hpe = htot(i)+h_ent
1906  mke_rate = 1.0/(1.0 + (cmke(1,i)*hpe + cmke(2,i)*hpe**2))
1907  tke_ent = tke_ent1 + dmke*ef4_val*mke_rate
1908  ! TKE_ent is the TKE that would remain if h_ent were entrained.
1909 
1910  dtke_dh = ((-idecay_len_tke(i)*tke_ent1 - drl*h_to_m) + &
1911  pen_dtke_dh_contrib*h_to_m) + dmke * mke_rate* &
1912  (def4_dh - ef4_val*mke_rate*(cmke(1,i)+2.0*cmke(2,i)*hpe))
1913  ! dh_Newt = -TKE_ent / dTKE_dh
1914  ! Bisect if the Newton's method prediction is outside of the bounded range.
1915  if (tke_ent > 0.0) then
1916  if ((h_max-h_ent)*(-dtke_dh) > tke_ent) then
1917  dh_newt = -tke_ent / dtke_dh
1918  else
1919  dh_newt = 0.5*(h_max-h_ent)
1920  endif
1921  h_min = h_ent
1922  else
1923  if ((h_min-h_ent)*(-dtke_dh) < tke_ent) then
1924  dh_newt = -tke_ent / dtke_dh
1925  else
1926  dh_newt = 0.5*(h_min-h_ent)
1927  endif
1928  h_max = h_ent
1929  endif
1930  h_ent = h_ent + dh_newt
1931 
1932  if (abs(dh_newt) < 0.2*gv%Angstrom) exit
1933  enddo
1934  endif
1935 
1936  if (h_ent < hmix_min-htot(i)) h_ent = hmix_min - htot(i)
1937 
1938  if (cs%TKE_diagnostics) then
1939  hpe = htot(i)+h_ent
1940  mke_rate = 1.0/(1.0 + cmke(1,i)*hpe + cmke(2,i)*hpe**2)
1941  ef4_val = ef4(htot(i)+h_neglect,h_ent,idecay_len_tke(i))
1942 
1943  e_hxhpe = h_ent / ((htot(i)+h_neglect)*(hpe+h_neglect))
1944  cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
1945  idt_diag * ((exp_kh-1.0)*tke(i) + &
1946  (h_ent*h_to_m)*drl*(1.0-f1_kh) + &
1947  dmke*mke_rate*(ef4_val-e_hxhpe))
1948  cs%diag_TKE_mixing(i,j) = cs%diag_TKE_mixing(i,j) - &
1949  idt_diag*(h_ent*h_to_m)*drl
1950  cs%diag_TKE_pen_SW(i,j) = cs%diag_TKE_pen_SW(i,j) - &
1951  idt_diag*(h_ent*h_to_m)*pen_en_contrib
1952  cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + &
1953  idt_diag*dmke*mke_rate*e_hxhpe
1954  endif
1955 
1956  tke(i) = 0.0
1957  endif ! TKE_full_ent > 0.0
1958 
1959  pen_absorbed = 0.0
1960  do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1961  sw_trans = exp(-h_ent*opacity_band(n,i,k))
1962  pen_absorbed = pen_absorbed + pen_sw_bnd(n,i) * (1.0 - sw_trans)
1963  pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
1964  endif ; enddo
1965 
1966  htot(i) = htot(i) + h_ent
1967  r0_tot(i) = r0_tot(i) + (h_ent * r0(i,k) + pen_absorbed*dr0_dt(i))
1968  h(i,k) = h(i,k) - h_ent
1969  d_eb(i,k) = d_eb(i,k) - h_ent
1970 
1971  stot(i) = stot(i) + h_ent * s(i,k)
1972  ttot(i) = ttot(i) + (h_ent * t(i,k) + pen_absorbed)
1973  rcv_tot(i) = rcv_tot(i) + (h_ent*rcv(i,k) + pen_absorbed*drcv_dt(i))
1974 
1975  uhtot(i) = uhtot(i) + u(i,k)*h_ent
1976  vhtot(i) = vhtot(i) + v(i,k)*h_ent
1977  endif ! h_avail > 0.0 .AND TKE(i) > 0.0
1978 
1979  endif ; enddo ! i loop
1980  enddo ! k loop
1981 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ mixedlayer_convection()

subroutine mom_bulk_mixed_layer::mixedlayer_convection ( real, dimension(szi_(g),szk_(gv)), intent(inout)  h,
real, dimension(szi_(g),szk_(gv)), intent(inout)  d_eb,
real, dimension(szi_(g)), intent(out)  htot,
real, dimension(szi_(g)), intent(out)  Ttot,
real, dimension(szi_(g)), intent(out)  Stot,
real, dimension(szi_(g)), intent(out)  uhtot,
real, dimension(szi_(g)), intent(out)  vhtot,
real, dimension(szi_(g)), intent(out)  R0_tot,
real, dimension(szi_(g)), intent(out)  Rcv_tot,
real, dimension(szi_(g),szk_(gv)), intent(in)  u,
real, dimension(szi_(g),szk_(gv)), intent(in)  v,
real, dimension(szi_(g),szk_(gv)), intent(in)  T,
real, dimension(szi_(g),szk_(gv)), intent(in)  S,
real, dimension(szi_(g),szk_(gv)), intent(in)  R0,
real, dimension(szi_(g),szk_(gv)), intent(in)  Rcv,
real, dimension(szi_(g),szk_(gv)), intent(in)  eps,
real, dimension(szi_(g)), intent(in)  dR0_dT,
real, dimension(szi_(g)), intent(in)  dRcv_dT,
real, dimension(szi_(g)), intent(in)  dR0_dS,
real, dimension(szi_(g)), intent(in)  dRcv_dS,
real, dimension(szi_(g)), intent(in)  netMassInOut,
real, dimension(szi_(g)), intent(in)  netMassOut,
real, dimension(szi_(g)), intent(in)  Net_heat,
real, dimension(szi_(g)), intent(in)  Net_salt,
integer, intent(in)  nsw,
real, dimension(:,:), intent(inout)  Pen_SW_bnd,
real, dimension(:,:,:), intent(in)  opacity_band,
real, dimension(szi_(g)), intent(out)  Conv_en,
real, dimension(szi_(g)), intent(out)  dKE_FC,
integer, intent(in)  j,
integer, dimension(szi_(g),szk_(gv)), intent(in)  ksort,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(bulkmixedlayer_cs), pointer  CS,
type(thermo_var_ptrs), intent(inout)  tv,
type(forcing), intent(inout)  fluxes,
real, intent(in)  dt,
logical, intent(in)  aggregate_FW_forcing 
)
private

This subroutine causes the mixed layer to entrain to the depth of free convection. The depth of free convection is the shallowest depth at which the fluid is denser than the average of the fluid above.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in,out]hLayer thickness, in m or kg m-2. (Intent in/out) The units of h are referred to as H below.
[in,out]d_ebThe downward increase across a layer in the entrainment from below , in H. Positive values go with mass gain by a layer.
[out]htotThe accumulated mixed layer thickness, in H.
[out]ttotThe depth integrated mixed layer temperature, in deg C H.
[out]stotThe depth integrated mixed layer salinity, in psu H.
[out]uhtotThe depth integrated mixed layer zonal velocity, H m s-1.
[out]vhtotThe integrated mixed layer meridional velocity, H m s-1.
[out]r0_totThe integrated mixed layer potential density referenced to 0 pressure, in kg m-2.
[out]rcv_totThe integrated mixed layer coordinate variable potential density, in kg m-2.
[in]nswThe number of bands of penetrating shortwave radiation.
[in,out]pen_sw_bndThe penetrating shortwave heating at the sea surface in each penetrating band, in K H, size nsw x SZI_(G).
[out]conv_enThe buoyant turbulent kinetic energy source due to free convection, in m3 s-2.
[out]dke_fcThe vertically integrated change in kinetic energy due to free convection, in m3 s-2.
[in]jThe j-index to work on.
[in]ksortThe density-sorted k-indices.
csThe control structure for this module.

Definition at line 1073 of file MOM_bulk_mixed_layer.F90.

Referenced by bulkmixedlayer().

1073  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
1074  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid
1075  !! structure.
1076  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness, in m or kg m-2.
1077  !! (Intent in/out) The units of h are
1078  !! referred to as H below.
1079  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_eb !< The downward increase across a
1080  !! layer in the entrainment from below
1081  !! , in H. Positive values go with
1082  !! mass gain by a layer.
1083  real, dimension(SZI_(G)), intent(out) :: htot !< The accumulated mixed layer
1084  !! thickness, in H.
1085  real, dimension(SZI_(G)), intent(out) :: ttot !< The depth integrated mixed layer
1086  !! temperature, in deg C H.
1087  real, dimension(SZI_(G)), intent(out) :: stot !< The depth integrated mixed layer
1088  !! salinity, in psu H.
1089  real, dimension(SZI_(G)), intent(out) :: uhtot !< The depth integrated mixed layer
1090  !! zonal velocity, H m s-1.
1091  real, dimension(SZI_(G)), intent(out) :: vhtot !< The integrated mixed layer
1092  !! meridional velocity, H m s-1.
1093  real, dimension(SZI_(G)), intent(out) :: r0_tot !< The integrated mixed layer
1094  !! potential density referenced to 0
1095  !! pressure, in kg m-2.
1096  real, dimension(SZI_(G)), intent(out) :: rcv_tot !< The integrated mixed layer
1097  !! coordinate variable potential
1098  !! density, in kg m-2.
1099  real, dimension(SZI_(G),SZK_(GV)), intent(in) :: u, v, t, s, r0, rcv, eps
1100  real, dimension(SZI_(G)), intent(in) :: dr0_dt, drcv_dt, dr0_ds, drcv_ds
1101  real, dimension(SZI_(G)), intent(in) :: netmassinout, netmassout
1102  real, dimension(SZI_(G)), intent(in) :: net_heat, net_salt
1103  integer, intent(in) :: nsw !< The number of bands of penetrating
1104  !! shortwave radiation.
1105  real, dimension(:,:), intent(inout) :: pen_sw_bnd !< The penetrating shortwave
1106  !! heating at the sea surface in each
1107  !! penetrating band, in K H,
1108  !! size nsw x SZI_(G).
1109  real, dimension(:,:,:), intent(in) :: opacity_band
1110  real, dimension(SZI_(G)), intent(out) :: conv_en !< The buoyant turbulent kinetic
1111  !! energy source due to free
1112  !! convection, in m3 s-2.
1113  real, dimension(SZI_(G)), intent(out) :: dke_fc !< The vertically integrated change
1114  !! in kinetic energy due to free
1115  !! convection, in m3 s-2.
1116  integer, intent(in) :: j !< The j-index to work on.
1117  integer, dimension(SZI_(G),SZK_(GV)), &
1118  intent(in) :: ksort !< The density-sorted k-indices.
1119  type(bulkmixedlayer_cs), pointer :: cs !< The control structure for this
1120  !! module.
1121  type(thermo_var_ptrs), intent(inout) :: tv
1122  type(forcing), intent(inout) :: fluxes
1123  real, intent(in) :: dt
1124  logical, intent(in) :: aggregate_fw_forcing
1125 
1126 ! This subroutine causes the mixed layer to entrain to the depth of free
1127 ! convection. The depth of free convection is the shallowest depth at which the
1128 ! fluid is denser than the average of the fluid above.
1129 ! Arguments: h - Layer thickness, in m or kg m-2. (Intent in/out) The units
1130 ! of h are referred to as H below.
1131 ! (in/out) d_eb - The downward increase across a layer in the entrainment from
1132 ! below, in H. Positive values go with mass gain by a layer.
1133 ! (out) htot - The accumulated mixed layer thickness, in H.
1134 ! (out) Ttot - The depth integrated mixed layer temperature, in deg C H.
1135 ! (out) Stot - The depth integrated mixed layer salinity, in psu H.
1136 ! (out) uhtot - The depth integrated mixed layer zonal velocity, H m s-1.
1137 ! (out) vhtot - The integrated mixed layer meridional velocity, H m s-1.
1138 ! (out) R0_tot - The integrated mixed layer potential density referenced
1139 ! to 0 pressure, in kg m-2.
1140 ! (out) Rcv_tot - The integrated mixed layer coordinate variable
1141 ! potential density, in kg m-2.
1142 ! (in) nsw - The number of bands of penetrating shortwave radiation.
1143 ! (out) Pen_SW_bnd - The penetrating shortwave heating at the sea surface
1144 ! in each penetrating band, in K H, size nsw x SZI_(G).
1145 ! (out) Conv_en - The buoyant turbulent kinetic energy source due to
1146 ! free convection, in m3 s-2.
1147 ! (out) dKE_FC - The vertically integrated change in kinetic energy due
1148 ! to free convection, in m3 s-2.
1149 ! (in) j - The j-index to work on.
1150 ! (in) ksort - The density-sorted k-indices.
1151 ! (in) G - The ocean's grid structure.
1152 ! (in) GV - The ocean's vertical grid structure.
1153 ! (in) CS - The control structure for this module.
1154 
1155  real, dimension(SZI_(G)) :: &
1156  massoutrem, & ! Evaporation that remains to be supplied, in H.
1157  netmassin ! mass entering through ocean surface (H)
1158  real :: sw_trans ! The fraction of shortwave radiation
1159  ! that is not absorbed in a layer, ND.
1160  real :: pen_absorbed ! The amount of penetrative shortwave radiation
1161  ! that is absorbed in a layer, in units of K H.
1162  real :: h_avail ! The thickness in a layer available for
1163  ! entrainment, in H.
1164  real :: h_ent ! The thickness from a layer that is entrained, in H.
1165  real :: t_precip ! The temperature of the precipitation, in deg C.
1166  real :: c1_3, c1_6 ! 1/3 and 1/6.
1167  real :: en_fn, frac, x1 ! Nondimensional temporary variables.
1168  real :: dr, dr0 ! Temporary variables with units of kg m-3 H.
1169  real :: dr_ent, dr_comp ! Temporary variables with units of kg m-3 H.
1170  real :: dr_dh ! The partial derivative of dr_ent with h_ent, in kg m-3.
1171  real :: h_min, h_max ! The minimum, maximum, and previous estimates for
1172  real :: h_prev ! h_ent, in H.
1173  real :: h_evap ! The thickness that is evaporated, in H.
1174  real :: dh_newt ! The Newton's method estimate of the change in
1175  ! h_ent between iterations, in H.
1176  real :: g_h2_2rho0 ! Half the gravitational acceleration times the
1177  ! square of the conversion from H to m divided
1178  ! by the mean density, in m6 s-2 H-2 kg-1.
1179  real :: angstrom ! The minimum layer thickness, in H.
1180  real :: opacity ! The opacity converted to units of H-1.
1181  real :: sum_pen_en ! The potential energy change due to penetrating
1182  ! shortwave radiation, integrated over a layer, in
1183  ! H kg m-3.
1184  real :: idt ! 1.0/dt
1185  real :: netheatout ! accumulated heat content of mass leaving ocean
1186  integer :: is, ie, nz, i, k, ks, itt, n
1187  real, dimension(max(nsw,1)) :: &
1188  c2, & ! Temporary variable with units of kg m-3 H-1.
1189  r_sw_top ! Temporary variables with units of H kg m-3.
1190 
1191  angstrom = gv%Angstrom
1192  c1_3 = 1.0/3.0 ; c1_6 = 1.0/6.0
1193  g_h2_2rho0 = (gv%g_Earth * gv%H_to_m**2) / (2.0 * gv%Rho0)
1194  idt = 1.0/dt
1195  is = g%isc ; ie = g%iec ; nz = gv%ke
1196 
1197  do i=is,ie ; if (ksort(i,1) > 0) then
1198  k = ksort(i,1)
1199 
1200  if (aggregate_fw_forcing) then
1201  massoutrem(i) = 0.0
1202  if (netmassinout(i) < 0.0) massoutrem(i) = -netmassinout(i)
1203  netmassin(i) = netmassinout(i) + massoutrem(i)
1204  else
1205  massoutrem(i) = -netmassout(i)
1206  netmassin(i) = netmassinout(i) - netmassout(i)
1207  endif
1208 
1209  ! htot is an Angstrom (taken from layer 1) plus any net precipitation.
1210  h_ent = max(min(angstrom,h(i,k)-eps(i,k)),0.0)
1211  htot(i) = h_ent + netmassin(i)
1212  h(i,k) = h(i,k) - h_ent
1213  d_eb(i,k) = d_eb(i,k) - h_ent
1214 
1215  pen_absorbed = 0.0
1216  do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1217  sw_trans = exp(-htot(i)*opacity_band(n,i,k))
1218  pen_absorbed = pen_absorbed + pen_sw_bnd(n,i) * (1.0-sw_trans)
1219  pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
1220  endif ; enddo
1221 
1222  ! Precipitation is assumed to have the same temperature and velocity
1223  ! as layer 1. Because layer 1 might not be the topmost layer, this
1224  ! involves multiple terms.
1225  t_precip = t(i,1)
1226  ttot(i) = (net_heat(i) + (netmassin(i) * t_precip + h_ent * t(i,k))) + &
1227  pen_absorbed
1228  ! Net_heat contains both heat fluxes and the heat content of mass fluxes.
1229  !! Ttot(i) = netMassIn(i) * T_precip + h_ent * T(i,k)
1230  !! Ttot(i) = Net_heat(i) + Ttot(i)
1231  !! Ttot(i) = Ttot(i) + Pen_absorbed
1232  ! smg:
1233  ! Ttot(i) = (Net_heat(i) + (h_ent * T(i,k))) + Pen_absorbed
1234  stot(i) = h_ent*s(i,k) + net_salt(i)
1235  uhtot(i) = u(i,1)*netmassin(i) + u(i,k)*h_ent
1236  vhtot(i) = v(i,1)*netmassin(i) + v(i,k)*h_ent
1237  r0_tot(i) = (h_ent*r0(i,k) + netmassin(i)*r0(i,1)) + &
1238 ! dR0_dT(i)*netMassIn(i)*(T_precip - T(i,1)) + &
1239  (dr0_dt(i)*(net_heat(i) + pen_absorbed) - &
1240  dr0_ds(i) * (netmassin(i) * s(i,1) - net_salt(i)))
1241  rcv_tot(i) = (h_ent*rcv(i,k) + netmassin(i)*rcv(i,1)) + &
1242 ! dRcv_dT(i)*netMassIn(i)*(T_precip - T(i,1)) + &
1243  (drcv_dt(i)*(net_heat(i) + pen_absorbed) - &
1244  drcv_ds(i) * (netmassin(i) * s(i,1) - net_salt(i)))
1245  conv_en(i) = 0.0 ; dke_fc(i) = 0.0
1246  if(ASSOCIATED(fluxes%heat_content_massin)) &
1247  fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) &
1248  + t_precip * netmassin(i) * gv%H_to_kg_m2 * fluxes%C_p * idt
1249  if (ASSOCIATED(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
1250  t_precip * netmassin(i) * gv%H_to_kg_m2
1251  endif ; enddo
1252 
1253  ! Now do netMassOut case in this block.
1254  ! At this point htot contains an Angstrom of fluid from layer 0 plus netMassIn.
1255  do ks=1,nz
1256  do i=is,ie ; if (ksort(i,ks) > 0) then
1257  k = ksort(i,ks)
1258 
1259  if ((htot(i) < angstrom) .and. (h(i,k) > eps(i,k))) then
1260  ! If less than an Angstrom was available from the layers above plus
1261  ! any precipitation, add more fluid from this layer.
1262  h_ent = min(angstrom-htot(i), h(i,k)-eps(i,k))
1263  htot(i) = htot(i) + h_ent
1264  h(i,k) = h(i,k) - h_ent
1265  d_eb(i,k) = d_eb(i,k) - h_ent
1266 
1267  r0_tot(i) = r0_tot(i) + h_ent*r0(i,k)
1268  uhtot(i) = uhtot(i) + h_ent*u(i,k)
1269  vhtot(i) = vhtot(i) + h_ent*v(i,k)
1270 
1271  rcv_tot(i) = rcv_tot(i) + h_ent*rcv(i,k)
1272  ttot(i) = ttot(i) + h_ent*t(i,k)
1273  stot(i) = stot(i) + h_ent*s(i,k)
1274  endif
1275 
1276  ! Water is removed from the topmost layers with any mass.
1277  ! We may lose layers if they are thin enough.
1278  ! The salt that is left behind goes into Stot.
1279  if ((massoutrem(i) > 0.0) .and. (h(i,k) > eps(i,k))) then
1280  if (massoutrem(i) > (h(i,k) - eps(i,k))) then
1281  h_evap = h(i,k) - eps(i,k)
1282  h(i,k) = eps(i,k)
1283  massoutrem(i) = massoutrem(i) - h_evap
1284  else
1285  h_evap = massoutrem(i)
1286  h(i,k) = h(i,k) - h_evap
1287  massoutrem(i) = 0.0
1288  endif
1289 
1290  stot(i) = stot(i) + h_evap*s(i,k)
1291  r0_tot(i) = r0_tot(i) + dr0_ds(i)*h_evap*s(i,k)
1292  rcv_tot(i) = rcv_tot(i) + drcv_ds(i)*h_evap*s(i,k)
1293  d_eb(i,k) = d_eb(i,k) - h_evap
1294 
1295  ! smg: when resolve the A=B code, we will set
1296  ! heat_content_massout = heat_content_massout - T(i,k)*h_evap*GV%H_to_kg_m2*fluxes%C_p*Idt
1297  ! by uncommenting the lines here.
1298  ! we will also then completely remove TempXpme from the model.
1299  if(ASSOCIATED(fluxes%heat_content_massout)) &
1300  fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) &
1301  - t(i,k)*h_evap*gv%H_to_kg_m2 * fluxes%C_p * idt
1302  if (ASSOCIATED(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) - &
1303  t(i,k)*h_evap*gv%H_to_kg_m2
1304 
1305  endif
1306 
1307  ! The following section calculates how much fluid will be entrained.
1308  h_avail = h(i,k) - eps(i,k)
1309  if (h_avail > 0.0) then
1310  dr = r0_tot(i) - htot(i)*r0(i,k)
1311  h_ent = 0.0
1312 
1313  dr0 = dr
1314  do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1315  dr0 = dr0 - (dr0_dt(i)*pen_sw_bnd(n,i)) * &
1316  opacity_band(n,i,k)*htot(i)
1317  endif ; enddo
1318 
1319  ! Some entrainment will occur from this layer.
1320  if (dr0 > 0.0) then
1321  dr_comp = dr
1322  do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1323  ! Compare the density at the bottom of a layer with the
1324  ! density averaged over the mixed layer and that layer.
1325  opacity = opacity_band(n,i,k)
1326  sw_trans = exp(-h_avail*opacity)
1327  dr_comp = dr_comp + (dr0_dt(i)*pen_sw_bnd(n,i)) * &
1328  ((1.0 - sw_trans) - opacity*(htot(i)+h_avail)*sw_trans)
1329  endif ; enddo
1330  if (dr_comp >= 0.0) then
1331  ! The entire layer is entrained.
1332  h_ent = h_avail
1333  else
1334  ! The layer is partially entrained. Iterate to determine how much
1335  ! entrainment occurs. Solve for the h_ent at which dr_ent = 0.
1336 
1337  ! Instead of assuming that the curve is linear between the two end
1338  ! points, assume that the change is concentrated near small values
1339  ! of entrainment. On average, this saves about 1 iteration.
1340  frac = dr0 / (dr0 - dr_comp)
1341  h_ent = h_avail * frac*frac
1342  h_min = 0.0 ; h_max = h_avail
1343 
1344  do n=1,nsw
1345  r_sw_top(n) = dr0_dt(i) * pen_sw_bnd(n,i)
1346  c2(n) = r_sw_top(n) * opacity_band(n,i,k)**2
1347  enddo
1348  do itt=1,10
1349  dr_ent = dr ; dr_dh = 0.0
1350  do n=1,nsw
1351  opacity = opacity_band(n,i,k)
1352  sw_trans = exp(-h_ent*opacity)
1353  dr_ent = dr_ent + r_sw_top(n) * ((1.0 - sw_trans) - &
1354  opacity*(htot(i)+h_ent)*sw_trans)
1355  dr_dh = dr_dh + c2(n) * (htot(i)+h_ent) * sw_trans
1356  enddo
1357 
1358  if (dr_ent > 0.0) then
1359  h_min = h_ent
1360  else
1361  h_max = h_ent
1362  endif
1363 
1364  dh_newt = -dr_ent / dr_dh
1365  h_prev = h_ent ; h_ent = h_prev+dh_newt
1366  if (h_ent > h_max) then
1367  h_ent = 0.5*(h_prev+h_max)
1368  else if (h_ent < h_min) then
1369  h_ent = 0.5*(h_prev+h_min)
1370  endif
1371 
1372  if (abs(dh_newt) < 0.2*angstrom) exit
1373  enddo
1374 
1375  endif
1376 
1377  ! Now that the amount of entrainment (h_ent) has been determined,
1378  ! calculate changes in various terms.
1379  sum_pen_en = 0.0 ; pen_absorbed = 0.0
1380  do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1381  opacity = opacity_band(n,i,k)
1382  sw_trans = exp(-h_ent*opacity)
1383 
1384  x1 = h_ent*opacity
1385  if (x1 < 2.0e-5) then
1386  en_fn = (opacity*htot(i)*(1.0 - 0.5*(x1 - c1_3*x1)) + &
1387  x1*x1*c1_6)
1388  else
1389  en_fn = ((opacity*htot(i) + 2.0) * &
1390  ((1.0-sw_trans) / x1) - 1.0 + sw_trans)
1391  endif
1392  sum_pen_en = sum_pen_en - (dr0_dt(i)*pen_sw_bnd(n,i)) * en_fn
1393 
1394  pen_absorbed = pen_absorbed + pen_sw_bnd(n,i) * (1.0 - sw_trans)
1395  pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
1396  endif ; enddo
1397 
1398  conv_en(i) = conv_en(i) + g_h2_2rho0 * h_ent * &
1399  ( (r0_tot(i) - r0(i,k)*htot(i)) + sum_pen_en )
1400 
1401  r0_tot(i) = r0_tot(i) + (h_ent * r0(i,k) + pen_absorbed*dr0_dt(i))
1402  stot(i) = stot(i) + h_ent * s(i,k)
1403  ttot(i) = ttot(i) + (h_ent * t(i,k) + pen_absorbed)
1404  rcv_tot(i) = rcv_tot(i) + (h_ent * rcv(i,k) + pen_absorbed*drcv_dt(i))
1405  endif ! dr0 > 0.0
1406 
1407  if (h_ent > 0.0) then
1408  if (htot(i) > 0.0) &
1409  dke_fc(i) = dke_fc(i) + cs%bulk_Ri_convective * 0.5 * &
1410  ((gv%H_to_m*h_ent) / (htot(i)*(h_ent+htot(i)))) * &
1411  ((uhtot(i)-u(i,k)*htot(i))**2 + (vhtot(i)-v(i,k)*htot(i))**2)
1412 
1413  htot(i) = htot(i) + h_ent
1414  h(i,k) = h(i,k) - h_ent
1415  d_eb(i,k) = d_eb(i,k) - h_ent
1416  uhtot(i) = u(i,k)*h_ent ; vhtot(i) = v(i,k)*h_ent
1417  endif
1418 
1419 
1420  endif ! h_avail>0
1421  endif ; enddo ! i loop
1422  enddo ! k loop
1423 
Here is the caller graph for this function:

◆ mixedlayer_detrain_1()

subroutine mom_bulk_mixed_layer::mixedlayer_detrain_1 ( real, dimension(szi_(g),szk0_(gv)), intent(inout)  h,
real, dimension(szi_(g),szk0_(gv)), intent(inout)  T,
real, dimension(szi_(g),szk0_(gv)), intent(inout)  S,
real, dimension(szi_(g),szk0_(gv)), intent(inout)  R0,
real, dimension(szi_(g),szk0_(gv)), intent(inout)  Rcv,
real, dimension(szk_(gv)), intent(in)  RcvTgt,
real, intent(in)  dt,
real, intent(in)  dt_diag,
real, dimension(szi_(g),szk_(gv)), intent(inout)  d_ea,
real, dimension(szi_(g),szk_(gv)), intent(inout)  d_eb,
integer, intent(in)  j,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(bulkmixedlayer_cs), pointer  CS,
real, dimension(szi_(g)), intent(in)  dRcv_dT,
real, dimension(szi_(g)), intent(in)  dRcv_dS,
real, dimension(szi_(g)), intent(in)  max_BL_det 
)
private

This subroutine moves any water left in the former mixed layers into the single buffer layers and may also move buffer layer water into the interior isopycnal layers.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in,out]hLayer thickness, in m or kg m-2. (Intent in/out) The units of h are referred to as H below. Layer 0 is the new mixed layer.
[in,out]tPotential temperature, in C.
[in,out]sSalinity, in psu.
[in,out]r0Potential density referenced to surface pressure, in kg m-3.
[in,out]rcvThe coordinate defining potential density, in kg m-3.
[in]rcvtgtThe target value of Rcv for each layer, in kg m-3.
[in]dtTime increment, in s.
[in,out]d_eaThe upward increase across a layer in the entrainment from above, in m or kg m-2 (H). Positive d_ea goes with layer thickness increases.
[in,out]d_ebThe downward increase across a layer in the entrainment from below, in H. Positive values go with mass gain by a layer.
[in]jThe meridional row to work on.
csThe control structure returned by a previous call to mixedlayer_init.
[in]drcv_dtThe partial derivative of coordinate defining potential density with potential temperature, in kg m-3 K-1.
[in]drcv_dsThe partial derivative of coordinate defining potential density with salinity, in kg m-3 psu-1.
[in]max_bl_detIf non-negative, the maximum detrainment permitted from the buffer layers, in H.

Definition at line 3332 of file MOM_bulk_mixed_layer.F90.

Referenced by bulkmixedlayer().

3332  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
3333  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
3334  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: h !< Layer thickness, in m or kg m-2.
3335  !! (Intent in/out) The units of h are
3336  !! referred to as H below. Layer 0 is
3337  !! the new mixed layer.
3338  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: t !< Potential temperature, in C.
3339  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: s !< Salinity, in psu.
3340  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: r0 !< Potential density referenced to
3341  !! surface pressure, in kg m-3.
3342  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: rcv !< The coordinate defining potential
3343  !! density, in kg m-3.
3344  real, dimension(SZK_(GV)), intent(in) :: rcvtgt !< The target value of Rcv for each
3345  !! layer, in kg m-3.
3346  real, intent(in) :: dt !< Time increment, in s.
3347  real, intent(in) :: dt_diag
3348  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_ea !< The upward increase across a layer in
3349  !! the entrainment from above, in m or
3350  !! kg m-2 (H). Positive d_ea goes with
3351  !! layer thickness increases.
3352  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_eb !< The downward increase across a layer
3353  !! in the entrainment from below, in H.
3354  !! Positive values go with mass gain by
3355  !! a layer.
3356  integer, intent(in) :: j !< The meridional row to work on.
3357  type(bulkmixedlayer_cs), pointer :: cs !< The control structure returned by a
3358  !! previous call to mixedlayer_init.
3359  real, dimension(SZI_(G)), intent(in) :: drcv_dt !< The partial derivative of
3360  !! coordinate defining potential density
3361  !! with potential temperature,
3362  !! in kg m-3 K-1.
3363  real, dimension(SZI_(G)), intent(in) :: drcv_ds !< The partial derivative of
3364  !! coordinate defining potential density
3365  !! with salinity, in kg m-3 psu-1.
3366  real, dimension(SZI_(G)), intent(in) :: max_bl_det !< If non-negative, the maximum
3367  !! detrainment permitted from the buffer
3368  !! layers, in H.
3369 
3370 ! This subroutine moves any water left in the former mixed layers into the
3371 ! single buffer layers and may also move buffer layer water into the interior
3372 ! isopycnal layers.
3373 
3374 ! Arguments: h - Layer thickness, in m or kg m-2. (Intent in/out) The units of
3375 ! h are referred to as H below. Layer 0 is the new mixed layer.
3376 ! (in/out) T - Potential temperature, in C.
3377 ! (in/out) S - Salinity, in psu.
3378 ! (in/out) R0 - Potential density referenced to surface pressure, in kg m-3.
3379 ! (in/out) Rcv - The coordinate defining potential density, in kg m-3.
3380 ! (in) RcvTgt - The target value of Rcv for each layer, in kg m-3.
3381 ! (in) dt - Time increment, in s.
3382 ! (in/out) d_ea - The upward increase across a layer in the entrainment from
3383 ! above, in m or kg m-2 (H). Positive d_ea goes with layer
3384 ! thickness increases.
3385 ! (in/out) d_eb - The downward increase across a layer in the entrainment from
3386 ! below, in H. Positive values go with mass gain by a layer.
3387 ! (in) j - The meridional row to work on.
3388 ! (in) G - The ocean's grid structure.
3389 ! (in) GV - The ocean's vertical grid structure.
3390 ! (in) CS - The control structure returned by a previous call to
3391 ! mixedlayer_init.
3392 ! (in) max_BL_det - If non-negative, the maximum detrainment permitted
3393 ! from the buffer layers, in H.
3394 ! (in/out) dRcv_dT - The partial derivative of coordinate defining potential
3395 ! density with potential temperature, in kg m-3 K-1.
3396 ! (in/out) dRcv_dS - The partial derivative of coordinate defining potential
3397 ! density with salinity, in kg m-3 psu-1.
3398  real :: ih ! The inverse of a thickness, in H-1.
3399  real :: h_ent ! The thickness from a layer that is
3400  ! entrained, in H.
3401  real :: max_det_rem(szi_(g)) ! Remaining permitted detrainment, in H.
3402  real :: detrain(szi_(g)) ! The thickness of fluid to detrain
3403  ! from the mixed layer, in H.
3404  real :: idt ! The inverse of the timestep in s-1.
3405  real :: dt_dr, ds_dr, drml, dr0_drcv, dt_ds_wt2
3406  real :: i_denom ! A work variable with units of psu2 m6 kg-2.
3407  real :: sdown, tdown
3408  real :: dt_time, timescale = 86400.0*30.0! *365.0/12.0
3409  real :: g_h2_2rho0dt ! Half the gravitational acceleration times the
3410  ! square of the conversion from H to m divided
3411  ! by the mean density times the time step, in m6 s-3 H-2 kg-1.
3412  real :: g_h2_2dt ! Half the gravitational acceleration times the
3413  ! square of the conversion from H to m divided
3414  ! by the diagnostic time step, in m3 H-2 s-3.
3415 
3416  logical :: splittable_bl(szi_(g)), orthogonal_extrap
3417  real :: x1
3418 
3419  integer :: i, is, ie, k, k1, nkmb, nz
3420  is = g%isc ; ie = g%iec ; nz = gv%ke
3421  nkmb = cs%nkml+cs%nkbl
3422  if (cs%nkbl /= 1) call mom_error(fatal,"MOM_mixed_layer: "// &
3423  "CS%nkbl must be 1 in mixedlayer_detrain_1.")
3424  idt = 1.0/dt
3425  dt_time = dt/timescale
3426  g_h2_2rho0dt = (gv%g_Earth * gv%H_to_m**2) / (2.0 * gv%Rho0 * dt_diag)
3427  g_h2_2dt = (gv%g_Earth * gv%H_to_m**2) / (2.0 * dt_diag)
3428 
3429  ! Move detrained water into the buffer layer.
3430  do k=1,cs%nkml
3431  do i=is,ie ; if (h(i,k) > 0.0) then
3432  ih = 1.0 / (h(i,nkmb) + h(i,k))
3433  if (cs%TKE_diagnostics) &
3434  cs%diag_TKE_conv_s2(i,j) = cs%diag_TKE_conv_s2(i,j) + &
3435  g_h2_2rho0dt * h(i,k) * h(i,nkmb) * &
3436  (r0(i,nkmb) - r0(i,k))
3437  if (ALLOCATED(cs%diag_PE_detrain)) &
3438  cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + &
3439  g_h2_2dt * h(i,k) * h(i,nkmb) * (r0(i,nkmb) - r0(i,k))
3440  if (ALLOCATED(cs%diag_PE_detrain2)) &
3441  cs%diag_PE_detrain2(i,j) = cs%diag_PE_detrain2(i,j) + &
3442  g_h2_2dt * h(i,k) * h(i,nkmb) * (r0(i,nkmb) - r0(i,k))
3443 
3444  r0(i,nkmb) = (r0(i,nkmb)*h(i,nkmb) + r0(i,k)*h(i,k)) * ih
3445  rcv(i,nkmb) = (rcv(i,nkmb)*h(i,nkmb) + rcv(i,k)*h(i,k)) * ih
3446  t(i,nkmb) = (t(i,nkmb)*h(i,nkmb) + t(i,k)*h(i,k)) * ih
3447  s(i,nkmb) = (s(i,nkmb)*h(i,nkmb) + s(i,k)*h(i,k)) * ih
3448 
3449  d_ea(i,k) = d_ea(i,k) - h(i,k)
3450  d_ea(i,nkmb) = d_ea(i,nkmb) + h(i,k)
3451  h(i,nkmb) = h(i,nkmb) + h(i,k)
3452  h(i,k) = 0.0
3453  endif ; enddo
3454  enddo
3455 
3456  do i=is,ie
3457  max_det_rem(i) = 10.0 * h(i,nkmb)
3458  if (max_bl_det(i) >= 0.0) max_det_rem(i) = max_bl_det(i)
3459  enddo
3460 
3461 ! If the mixed layer was denser than the densest interior layer,
3462 ! but is now lighter than this layer, leaving a buffer layer that
3463 ! is denser than this layer, there are problems. This should prob-
3464 ! ably be considered a case of an inadequate choice of resolution in
3465 ! density space and should be avoided. To make the model run sens-
3466 ! ibly in this case, it will make the mixed layer denser while making
3467 ! the buffer layer the density of the densest interior layer (pro-
3468 ! vided that the this will not make the mixed layer denser than the
3469 ! interior layer). Otherwise, make the mixed layer the same density
3470 ! as the densest interior layer and lighten the buffer layer with
3471 ! the released buoyancy. With multiple buffer layers, much more
3472 ! graceful options are available.
3473  do i=is,ie ; if (h(i,nkmb) > 0.0) then
3474  if ((r0(i,0)<r0(i,nz)) .and. (r0(i,nz)<r0(i,nkmb))) then
3475  if ((r0(i,nz)-r0(i,0))*h(i,0) > &
3476  (r0(i,nkmb)-r0(i,nz))*h(i,nkmb)) then
3477  detrain(i) = (r0(i,nkmb)-r0(i,nz))*h(i,nkmb) / (r0(i,nkmb)-r0(i,0))
3478  else
3479  detrain(i) = (r0(i,nz)-r0(i,0))*h(i,0) / (r0(i,nkmb)-r0(i,0))
3480  endif
3481 
3482  d_eb(i,cs%nkml) = d_eb(i,cs%nkml) + detrain(i)
3483  d_ea(i,cs%nkml) = d_ea(i,cs%nkml) - detrain(i)
3484  d_eb(i,nkmb) = d_eb(i,nkmb) - detrain(i)
3485  d_ea(i,nkmb) = d_ea(i,nkmb) + detrain(i)
3486 
3487  if (ALLOCATED(cs%diag_PE_detrain)) cs%diag_PE_detrain(i,j) = &
3488  cs%diag_PE_detrain(i,j) + g_h2_2dt * detrain(i)* &
3489  (h(i,0) + h(i,nkmb)) * (r0(i,nkmb) - r0(i,0))
3490  x1 = r0(i,0)
3491  r0(i,0) = r0(i,0) - detrain(i)*(r0(i,0)-r0(i,nkmb)) / h(i,0)
3492  r0(i,nkmb) = r0(i,nkmb) - detrain(i)*(r0(i,nkmb)-x1) / h(i,nkmb)
3493  x1 = rcv(i,0)
3494  rcv(i,0) = rcv(i,0) - detrain(i)*(rcv(i,0)-rcv(i,nkmb)) / h(i,0)
3495  rcv(i,nkmb) = rcv(i,nkmb) - detrain(i)*(rcv(i,nkmb)-x1) / h(i,nkmb)
3496  x1 = t(i,0)
3497  t(i,0) = t(i,0) - detrain(i)*(t(i,0)-t(i,nkmb)) / h(i,0)
3498  t(i,nkmb) = t(i,nkmb) - detrain(i)*(t(i,nkmb)-x1) / h(i,nkmb)
3499  x1 = s(i,0)
3500  s(i,0) = s(i,0) - detrain(i)*(s(i,0)-s(i,nkmb)) / h(i,0)
3501  s(i,nkmb) = s(i,nkmb) - detrain(i)*(s(i,nkmb)-x1) / h(i,nkmb)
3502 
3503  endif
3504  endif ; enddo
3505 
3506  ! Move water out of the buffer layer, if convenient.
3507 ! Split the buffer layer if possible, and replace the buffer layer
3508 ! with a small amount of fluid from the mixed layer.
3509 ! This is the exponential-in-time splitting, circa 2005.
3510  do i=is,ie
3511  if (h(i,nkmb) > 0.0) then ; splittable_bl(i) = .true.
3512  else ; splittable_bl(i) = .false. ; endif
3513  enddo
3514 
3515  dt_ds_wt2 = cs%dT_dS_wt**2
3516 
3517 ! dt_Time = dt/Timescale
3518  do k=nz-1,nkmb+1,-1 ; do i=is,ie
3519  if (splittable_bl(i)) then
3520  if (rcvtgt(k)<=rcv(i,nkmb)) then
3521 ! Estimate dR/drho, dTheta/dR, and dS/dR, where R is the coordinate variable
3522 ! and rho is in-situ (or surface) potential density.
3523 ! There is no "right" way to do this, so this keeps things reasonable, if
3524 ! slightly arbitrary.
3525  splittable_bl(i) = .false.
3526 
3527  k1 = k+1 ; orthogonal_extrap = .false.
3528  ! Here we try to find a massive layer to use for interpolating the
3529  ! temperature and salinity. If none is available a pseudo-orthogonal
3530  ! extrapolation is used. The 10.0 and 0.9 in the following are
3531  ! arbitrary but probably about right.
3532  if ((h(i,k+1) < 10.0*gv%Angstrom) .or. &
3533  ((rcvtgt(k+1)-rcv(i,nkmb)) >= 0.9*(rcv(i,k1) - rcv(i,0)))) then
3534  if (k>=nz-1) then ; orthogonal_extrap = .true.
3535  elseif ((h(i,k+2) <= 10.0*gv%Angstrom) .and. &
3536  ((rcvtgt(k+1)-rcv(i,nkmb)) < 0.9*(rcv(i,k+2)-rcv(i,0)))) then
3537  k1 = k+2
3538  else ; orthogonal_extrap = .true. ; endif
3539  endif
3540 
3541  if ((r0(i,0) >= r0(i,k1)) .or. (rcv(i,0) >= rcv(i,nkmb))) cycle
3542  ! In this case there is an inversion of in-situ density relative to
3543  ! the coordinate variable. Do not detrain from the buffer layer.
3544 
3545  if (orthogonal_extrap) then
3546  ! 36 here is a typical oceanic value of (dR/dS) / (dR/dT) - it says
3547  ! that the relative weights of T & S changes is a plausible 6:1.
3548  ! Also, this was coded on Athena's 6th birthday!
3549  i_denom = 1.0 / (drcv_ds(i)**2 + dt_ds_wt2*drcv_dt(i)**2)
3550  dt_dr = dt_ds_wt2*drcv_dt(i) * i_denom
3551  ds_dr = drcv_ds(i) * i_denom
3552  else
3553  dt_dr = (t(i,0) - t(i,k1)) / (rcv(i,0) - rcv(i,k1))
3554  ds_dr = (s(i,0) - s(i,k1)) / (rcv(i,0) - rcv(i,k1))
3555  endif
3556  drml = dt_time * (r0(i,nkmb) - r0(i,0)) * &
3557  (rcv(i,0) - rcv(i,k1)) / (r0(i,0) - r0(i,k1))
3558  ! Once again, there is an apparent density inversion in Rcv.
3559  if (drml < 0.0) cycle
3560  dr0_drcv = (r0(i,0) - r0(i,k1)) / (rcv(i,0) - rcv(i,k1))
3561 
3562  if ((rcv(i,nkmb) - drml < rcvtgt(k)) .and. (max_det_rem(i) > h(i,nkmb))) then
3563  ! In this case, the buffer layer is split into two isopycnal layers.
3564  detrain(i) = h(i,nkmb)*(rcv(i,nkmb) - rcvtgt(k)) / &
3565  (rcvtgt(k+1) - rcvtgt(k))
3566 
3567  if (ALLOCATED(cs%diag_PE_detrain)) cs%diag_PE_detrain(i,j) = &
3568  cs%diag_PE_detrain(i,j) - g_h2_2dt * detrain(i) * &
3569  (h(i,nkmb)-detrain(i)) * (rcvtgt(k+1) - rcvtgt(k)) * dr0_drcv
3570 
3571  tdown = detrain(i) * (t(i,nkmb) + dt_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3572  t(i,k) = (h(i,k) * t(i,k) + &
3573  (h(i,nkmb) * t(i,nkmb) - tdown)) / &
3574  (h(i,k) + (h(i,nkmb) - detrain(i)))
3575  t(i,k+1) = (h(i,k+1) * t(i,k+1) + tdown)/ &
3576  (h(i,k+1) + detrain(i))
3577  t(i,nkmb) = t(i,0)
3578  sdown = detrain(i) * (s(i,nkmb) + ds_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3579  s(i,k) = (h(i,k) * s(i,k) + &
3580  (h(i,nkmb) * s(i,nkmb) - sdown)) / &
3581  (h(i,k) + (h(i,nkmb) - detrain(i)))
3582  s(i,k+1) = (h(i,k+1) * s(i,k+1) + sdown)/ &
3583  (h(i,k+1) + detrain(i))
3584  s(i,nkmb) = s(i,0)
3585  rcv(i,nkmb) = rcv(i,0)
3586 
3587  d_ea(i,k+1) = d_ea(i,k+1) + detrain(i)
3588  d_ea(i,k) = d_ea(i,k) + (h(i,nkmb) - detrain(i))
3589  d_ea(i,nkmb) = d_ea(i,nkmb) - h(i,nkmb)
3590 
3591  h(i,k+1) = h(i,k+1) + detrain(i)
3592  h(i,k) = h(i,k) + h(i,nkmb) - detrain(i)
3593  h(i,nkmb) = 0.0
3594  else
3595  ! Here only part of the buffer layer is moved into the interior.
3596  detrain(i) = h(i,nkmb) * drml / (rcvtgt(k+1) - rcv(i,nkmb) + drml)
3597  if (detrain(i) > max_det_rem(i)) detrain(i) = max_det_rem(i)
3598  ih = 1.0 / (h(i,k+1) + detrain(i))
3599 
3600  tdown = (t(i,nkmb) + dt_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3601  t(i,nkmb) = t(i,nkmb) - dt_dr * drml
3602  t(i,k+1) = (h(i,k+1) * t(i,k+1) + detrain(i) * tdown) * ih
3603  sdown = (s(i,nkmb) + ds_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3604 ! The following two expressions updating S(nkmb) are mathematically identical.
3605 ! S(i,nkmb) = (h(i,nkmb) * S(i,nkmb) - detrain(i) * Sdown) / &
3606 ! (h(i,nkmb) - detrain(i))
3607  s(i,nkmb) = s(i,nkmb) - ds_dr * drml
3608  s(i,k+1) = (h(i,k+1) * s(i,k+1) + detrain(i) * sdown) * ih
3609 
3610  d_ea(i,k+1) = d_ea(i,k+1) + detrain(i)
3611  d_ea(i,nkmb) = d_ea(i,nkmb) - detrain(i)
3612 
3613  h(i,k+1) = h(i,k+1) + detrain(i)
3614  h(i,nkmb) = h(i,nkmb) - detrain(i)
3615 
3616  if (ALLOCATED(cs%diag_PE_detrain)) cs%diag_PE_detrain(i,j) = &
3617  cs%diag_PE_detrain(i,j) - g_h2_2dt * detrain(i) * dr0_drcv * &
3618  (h(i,nkmb)-detrain(i)) * (rcvtgt(k+1) - rcv(i,nkmb) + drml)
3619  endif
3620  endif ! RcvTgt(k)<=Rcv(i,nkmb)
3621  endif ! splittable_BL
3622  enddo ; enddo ! i & k loops
3623 
3624 ! The numerical behavior of the buffer layer is dramatically improved
3625 ! if it is always at least a small fraction (say 10%) of the thickness
3626 ! of the mixed layer. As the physical distinction between the mixed
3627 ! and buffer layers is vague anyway, this seems hard to argue against.
3628  do i=is,ie
3629  if (h(i,nkmb) < 0.1*h(i,0)) then
3630  h_ent = 0.1*h(i,0) - h(i,nkmb)
3631  ih = 10.0/h(i,0)
3632  t(i,nkmb) = (h(i,nkmb)*t(i,nkmb) + h_ent*t(i,0)) * ih
3633  s(i,nkmb) = (h(i,nkmb)*s(i,nkmb) + h_ent*s(i,0)) * ih
3634 
3635  d_ea(i,1) = d_ea(i,1) - h_ent
3636  d_ea(i,nkmb) = d_ea(i,nkmb) + h_ent
3637 
3638  h(i,0) = h(i,0) - h_ent
3639  h(i,nkmb) = h(i,nkmb) + h_ent
3640  endif
3641  enddo
3642 
Here is the caller graph for this function:

◆ mixedlayer_detrain_2()

subroutine mom_bulk_mixed_layer::mixedlayer_detrain_2 ( real, dimension(szi_(g),szk0_(gv)), intent(inout)  h,
real, dimension(szi_(g),szk0_(gv)), intent(inout)  T,
real, dimension(szi_(g),szk0_(gv)), intent(inout)  S,
real, dimension(szi_(g),szk0_(gv)), intent(inout)  R0,
real, dimension(szi_(g),szk0_(gv)), intent(inout)  Rcv,
real, dimension(szk_(gv)), intent(in)  RcvTgt,
real, intent(in)  dt,
real, intent(in)  dt_diag,
real, dimension(szi_(g),szk_(gv)), intent(inout)  d_ea,
integer, intent(in)  j,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(bulkmixedlayer_cs), pointer  CS,
real, dimension(szi_(g)), intent(in)  dR0_dT,
real, dimension(szi_(g)), intent(in)  dR0_dS,
real, dimension(szi_(g)), intent(in)  dRcv_dT,
real, dimension(szi_(g)), intent(in)  dRcv_dS,
real, dimension(szi_(g)), intent(in)  max_BL_det 
)
private

This subroutine moves any water left in the former mixed layers into the two buffer layers and may also move buffer layer water into the interior isopycnal layers.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in,out]hLayer thickness, in m or kg m-2. (Intent in/out) The units of h are referred to as H below. Layer 0 is the new mixed layer.
[in,out]tPotential temperature, in C.
[in,out]sSalinity, in psu.
[in,out]r0Potential density referenced to surface pressure, in kg m-3.
[in,out]rcvThe coordinate defining potential density, in kg m-3.
[in]rcvtgtThe target value of Rcv for each layer, in kg m-3.
[in]dtTime increment, in s.
[in]dt_diagThe diagnostic time step, in s.
[in,out]d_eaThe upward increase across a layer in the entrainment from above, in m or kg m-2 (H). Positive d_ea goes with layer thickness increases.
[in]jThe meridional row to work on.
csThe control structure returned by a previous call to mixedlayer_init.
[in]dr0_dtThe partial derivative of potential density referenced to the surface with potential temperature, in kg m-3 K-1.
[in]dr0_dsThe partial derivative of cpotential density referenced to the surface with salinity, in kg m-3 psu-1.
[in]drcv_dtThe partial derivative of coordinate defining potential density with potential temperature, in kg m-3 K-1.
[in]drcv_dsThe partial derivative of coordinate defining potential density with salinity, in kg m-3 psu-1.
[in]max_bl_detIf non-negative, the maximum detrainment permitted from the buffer layers, in H.

Definition at line 2408 of file MOM_bulk_mixed_layer.F90.

References mom_error_handler::mom_error().

Referenced by bulkmixedlayer().

2408  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
2409  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
2410  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: h !< Layer thickness, in m or kg m-2.
2411  !! (Intent in/out) The units of h are
2412  !! referred to as H below.
2413  !! Layer 0 is the new mixed layer.
2414  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: t !< Potential temperature, in C.
2415  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: s !< Salinity, in psu.
2416  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: r0 !< Potential density referenced to
2417  !! surface pressure, in kg m-3.
2418  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: rcv !< The coordinate defining potential
2419  !! density, in kg m-3.
2420  real, dimension(SZK_(GV)), intent(in) :: rcvtgt !< The target value of Rcv for each
2421  !! layer, in kg m-3.
2422  real, intent(in) :: dt !< Time increment, in s.
2423  real, intent(in) :: dt_diag !< The diagnostic time step, in s.
2424  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_ea !< The upward increase across a layer in
2425  !! the entrainment from above, in m or
2426  !! kg m-2 (H). Positive d_ea goes with
2427  !! layer thickness increases.
2428  integer, intent(in) :: j !< The meridional row to work on.
2429  type(bulkmixedlayer_cs), pointer :: cs !< The control structure returned by a
2430  !! previous call to mixedlayer_init.
2431  real, dimension(SZI_(G)), intent(in) :: dr0_dt !< The partial derivative of
2432  !! potential density referenced to the
2433  !! surface with potential temperature,
2434  !! in kg m-3 K-1.
2435  real, dimension(SZI_(G)), intent(in) :: dr0_ds !< The partial derivative of
2436  !! cpotential density referenced to the
2437  !! surface with salinity,
2438  !! in kg m-3 psu-1.
2439  real, dimension(SZI_(G)), intent(in) :: drcv_dt !< The partial derivative of
2440  !! coordinate defining potential density
2441  !! with potential temperature,
2442  !! in kg m-3 K-1.
2443  real, dimension(SZI_(G)), intent(in) :: drcv_ds !< The partial derivative of
2444  !! coordinate defining potential density
2445  !! with salinity, in kg m-3 psu-1.
2446  real, dimension(SZI_(G)), intent(in) :: max_bl_det !< If non-negative, the maximum
2447  !! detrainment permitted from the buffer
2448  !! layers, in H.
2449 
2450 ! This subroutine moves any water left in the former mixed layers into the
2451 ! two buffer layers and may also move buffer layer water into the interior
2452 ! isopycnal layers.
2453 
2454 ! Arguments: h - Layer thickness, in m or kg m-2. (Intent in/out) The units of
2455 ! h are referred to as H below. Layer 0 is the new mixed layer.
2456 ! (in/out) T - Potential temperature, in C.
2457 ! (in/out) S - Salinity, in psu.
2458 ! (in/out) R0 - Potential density referenced to surface pressure, in kg m-3.
2459 ! (in/out) Rcv - The coordinate defining potential density, in kg m-3.
2460 ! (in) RcvTgt - The target value of Rcv for each layer, in kg m-3.
2461 ! (in) dt - Time increment, in s.
2462 ! (in) dt_diag - The diagnostic time step, in s.
2463 ! (in/out) d_ea - The upward increase across a layer in the entrainment from
2464 ! above, in m or kg m-2 (H). Positive d_ea goes with layer
2465 ! thickness increases.
2466 ! (in) j - The meridional row to work on.
2467 ! (in) G - The ocean's grid structure.
2468 ! (in) GV - The ocean's vertical grid structure.
2469 ! (in) CS - The control structure returned by a previous call to
2470 ! mixedlayer_init.
2471 ! (in) max_BL_det - If non-negative, the maximum detrainment permitted
2472 ! from the buffer layers, in H.
2473 ! (in) dR0_dT - The partial derivative of potential density referenced
2474 ! to the surface with potential temperature, in kg m-3 K-1.
2475 ! (in) dR0_dS - The partial derivative of cpotential density referenced
2476 ! to the surface with salinity, in kg m-3 psu-1.
2477 ! (in) dRcv_dT - The partial derivative of coordinate defining potential
2478 ! density with potential temperature, in kg m-3 K-1.
2479 ! (in) dRcv_dS - The partial derivative of coordinate defining potential
2480 ! density with salinity, in kg m-3 psu-1.
2481  real :: h_to_bl ! The total thickness detrained to the buffer
2482  ! layers, in H (the units of h).
2483  real :: r0_to_bl, rcv_to_bl ! The depth integrated amount of R0, Rcv, T
2484  real :: t_to_bl, s_to_bl ! and S that is detrained to the buffer layer,
2485  ! in H kg m-3, H kg m-3, K H, and psu H.
2486 
2487  real :: h_min_bl ! The minimum buffer layer thickness, in H.
2488  real :: h_min_bl_thick ! The minimum buffer layer thickness when the
2489  ! mixed layer is very large, in H.
2490  real :: h_min_bl_frac_ml = 0.05 ! The minimum buffer layer thickness relative
2491  ! to the total mixed layer thickness for thin
2492  ! mixed layers, nondim., maybe 0.1/CS%nkbl.
2493 
2494  real :: h1, h2 ! Scalar variables holding the values of
2495  ! h(i,CS%nkml+1) and h(i,CS%nkml+2), in H.
2496  real :: h1_avail ! The thickess of the upper buffer layer
2497  ! available to move into the lower buffer
2498  ! layer, in H.
2499  real :: stays ! stays is the thickness of the upper buffer
2500  ! layer that remains there, in units of H.
2501  real :: stays_min, stays_max ! The minimum and maximum permitted values of
2502  ! stays, in units of H.
2503 
2504  logical :: mergeable_bl ! If true, it is an option to combine the two
2505  ! buffer layers and create water that matches
2506  ! the target density of an interior layer.
2507  real :: stays_merge ! If the two buffer layers can be combined
2508  ! stays_merge is the thickness of the upper
2509  ! layer that remains, in units of H.
2510  real :: stays_min_merge ! The minimum allowed value of stays_merge in H.
2511 
2512  real :: dr0_2dz, drcv_2dz ! Half the vertical gradients of R0, Rcv, T, and
2513 ! real :: dT_2dz, dS_2dz ! S, in kg m-4, kg m-4, K m-1, and psu m-1.
2514  real :: scale_slope ! A nondimensional number < 1 used to scale down
2515  ! the slope within the upper buffer layer when
2516  ! water MUST be detrained to the lower layer.
2517 
2518  real :: dpe_extrap ! The potential energy change due to dispersive
2519  ! advection or mixing layers, divided by
2520  ! rho_0*g, in units of H2.
2521  real :: dpe_det, dpe_merge ! The energy required to mix the detrained water
2522  ! into the buffer layer or the merge the two
2523  ! buffer layers, both in units of J H2 m-4.
2524 
2525  real :: h_from_ml ! The amount of additional water that must be
2526  ! drawn from the mixed layer, in H.
2527  real :: h_det_h2 ! The amount of detrained water and mixed layer
2528  ! water that will go directly into the lower
2529  ! buffer layer, in H.
2530  real :: h_det_to_h2, h_ml_to_h2 ! All of the variables hA_to_hB are the
2531  real :: h_det_to_h1, h_ml_to_h1 ! thickess fluxes from one layer to another,
2532  real :: h1_to_h2, h1_to_k0 ! in H, with h_det the detrained water, h_ml
2533  real :: h2_to_k1, h2_to_k1_rem ! the actively mixed layer, h1 and h2 the upper
2534  ! and lower buffer layers, and k0 and k1 the
2535  ! interior layers that are just lighter and
2536  ! just denser than the lower buffer layer.
2537 
2538  real :: r0_det, t_det, s_det ! Detrained values of R0, T, and S.
2539  real :: rcv_stays, r0_stays ! Values of Rcv and R0 that stay in a layer.
2540  real :: t_stays, s_stays ! Values of T and S that stay in a layer.
2541  real :: dspice_det, dspice_stays! The spiciness difference between an original
2542  ! buffer layer and the water that moves into
2543  ! an interior layer or that stays in that
2544  ! layer, in kg m-3.
2545  real :: dspice_lim, dspice_lim2 ! Limits to the spiciness difference between
2546  ! the lower buffer layer and the water that
2547  ! moves into an interior layer, in kg m-3.
2548  real :: dspice_2dz ! The vertical gradient of spiciness used for
2549  ! advection, in kg m-4.
2550 
2551  real :: dpe_ratio ! Multiplier of dPE_det at which merging is
2552  ! permitted - here (detrainment_per_day/dt)*30
2553  ! days?
2554  real :: num_events ! The number of detrainment events over which
2555  ! to prefer merging the buffer layers.
2556  real :: detrainment_timescale ! The typical timescale for a detrainment
2557  ! event, in s.
2558  real :: dpe_time_ratio ! Larger of 1 and the detrainment_timescale
2559  ! over dt, nondimensional.
2560  real :: dt_ds_gauge, ds_dt_gauge ! The relative scales of temperature and
2561  ! salinity changes in defining spiciness, in
2562  ! K psu-1 and psu K-1.
2563  real :: i_denom ! A work variable with units of psu2 m6 kg-2.
2564 
2565  real :: g_2 ! 1/2 G_Earth, in m s-2.
2566  real :: rho0xg ! Rho0 times G_Earth, in kg m-2 s-2.
2567  real :: i2rho0 ! 1 / (2 Rho0), in m3 kg-1.
2568  real :: idt_h2 ! The square of the conversion from thickness
2569  ! to m divided by the time step in m2 H-2 s-1.
2570  logical :: stable_rcv ! If true, the buffer layers are stable with
2571  ! respect to the coordinate potential density.
2572  real :: h_neglect ! A thickness that is so small it is usually lost
2573  ! in roundoff and can be neglected, in H.
2574 
2575  real :: s1en ! A work variable with units of H2 kg m s-3.
2576  real :: s1, s2, bh0 ! Work variables with units of H.
2577  real :: s3sq ! A work variable with units of H2.
2578  real :: i_ya, b1 ! Nondimensional work variables.
2579  real :: ih, ihdet, ih1f, ih2f ! Assorted inverse thickness work variables,
2580  real :: ihk0, ihk1, ih12 ! all with units of H-1.
2581  real :: dr1, dr2, dr2b, drk1 ! Assorted density difference work variables,
2582  real :: dr0, dr21, drcv ! all with units of kg m-3.
2583  real :: drcv_stays, drcv_det, drcv_lim
2584  real :: angstrom ! The minumum layer thickness, in H.
2585 
2586  real :: h2_to_k1_lim, t_new, s_new, t_max, t_min, s_max, s_min
2587  character(len=200) :: mesg
2588 
2589  integer :: i, k, k0, k1, is, ie, nz, kb1, kb2, nkmb
2590  is = g%isc ; ie = g%iec ; nz = gv%ke
2591  kb1 = cs%nkml+1; kb2 = cs%nkml+2
2592  nkmb = cs%nkml+cs%nkbl
2593  h_neglect = gv%H_subroundoff
2594  g_2 = 0.5*gv%g_Earth
2595  rho0xg = gv%Rho0 * gv%g_Earth
2596  idt_h2 = gv%H_to_m**2 / dt_diag
2597  i2rho0 = 0.5 / gv%Rho0
2598  angstrom = gv%Angstrom
2599 
2600  ! This is hard coding of arbitrary and dimensional numbers.
2601  h_min_bl_thick = 5.0 * gv%m_to_H
2602  dt_ds_gauge = cs%dT_dS_wt ; ds_dt_gauge = 1.0 /dt_ds_gauge
2603  num_events = 10.0
2604  detrainment_timescale = 4.0*3600.0
2605 
2606  if (cs%nkbl /= 2) call mom_error(fatal, "MOM_mixed_layer"// &
2607  "CS%nkbl must be 2 in mixedlayer_detrain_2.")
2608 
2609  if (dt < detrainment_timescale) then ; dpe_time_ratio = detrainment_timescale/dt
2610  else ; dpe_time_ratio = 1.0 ; endif
2611 
2612  do i=is,ie
2613 
2614  ! Determine all of the properties being detrained from the mixed layer.
2615 
2616  ! As coded this has the k and i loop orders switched, but k is CS%nkml is
2617  ! often just 1 or 2, so this seems like it should not be a problem, especially
2618  ! since it means that a number of variables can now be scalars, not arrays.
2619  h_to_bl = 0.0 ; r0_to_bl = 0.0
2620  rcv_to_bl = 0.0 ; t_to_bl = 0.0 ; s_to_bl = 0.0
2621 
2622  do k=1,cs%nkml ; if (h(i,k) > 0.0) then
2623  h_to_bl = h_to_bl + h(i,k)
2624  r0_to_bl = r0_to_bl + r0(i,k)*h(i,k)
2625 
2626  rcv_to_bl = rcv_to_bl + rcv(i,k)*h(i,k)
2627  t_to_bl = t_to_bl + t(i,k)*h(i,k)
2628  s_to_bl = s_to_bl + s(i,k)*h(i,k)
2629 
2630  d_ea(i,k) = d_ea(i,k) - h(i,k)
2631  h(i,k) = 0.0
2632  endif ; enddo
2633  if (h_to_bl > 0.0) then ; r0_det = r0_to_bl / h_to_bl
2634  else ; r0_det = r0(i,0) ; endif
2635 
2636  ! This code does both downward detrainment from both the mixed layer and the
2637  ! buffer layers.
2638  ! Several considerations apply in detraining water into the interior:
2639  ! (1) Water only moves into the interior from the deeper buffer layer,
2640  ! so the deeper buffer layer must have some mass.
2641  ! (2) The upper buffer layer must have some mass so the extrapolation of
2642  ! density is meaningful (i.e. there is not detrainment from the buffer
2643  ! layers when there is strong mixed layer entrainment).
2644  ! (3) The lower buffer layer density extrapolated to its base with a
2645  ! linear fit between the two layers must exceed the density of the
2646  ! next denser interior layer.
2647  ! (4) The average extroplated coordinate density that is moved into the
2648  ! isopycnal interior matches the target value for that layer.
2649  ! (5) The potential energy change is calculated and might be used later
2650  ! to allow the upper buffer layer to mix more into the lower buffer
2651  ! layer.
2652 
2653  ! Determine whether more must be detrained from the mixed layer to keep a
2654  ! minimal amount of mass in the buffer layers. In this case the 5% of the
2655  ! mixed layer thickness is hard-coded, but probably shouldn't be!
2656  h_min_bl = min(h_min_bl_thick,h_min_bl_frac_ml*h(i,0))
2657 
2658  stable_rcv = .true.
2659  if (((r0(i,kb2)-r0(i,kb1)) * (rcv(i,kb2)-rcv(i,kb1)) <= 0.0)) &
2660  stable_rcv = .false.
2661 
2662  h1 = h(i,kb1) ; h2 = h(i,kb2)
2663 
2664  h2_to_k1_rem = (h1 + h2) + h_to_bl
2665  if ((max_bl_det(i) >= 0.0) .and. (h2_to_k1_rem > max_bl_det(i))) &
2666  h2_to_k1_rem = max_bl_det(i)
2667 
2668 
2669  if ((h2 == 0.0) .and. (h1 > 0.0)) then
2670  ! The lower buffer layer has been eliminated either by convective
2671  ! adjustment or entrainment from the interior, and its current properties
2672  ! are not meaningful, but may later be used to determine the properties of
2673  ! waters moving into the lower buffer layer. So the properties of the
2674  ! lower buffer layer are set to be between those of the upper buffer layer
2675  ! and the next denser interior layer, measured by R0. This probably does
2676  ! not happen very often, so I am not too worried about the inefficiency of
2677  ! the following loop.
2678  do k1=kb2+1,nz ; if (h(i,k1) > 2.0*angstrom) exit ; enddo
2679 
2680  r0(i,kb2) = r0(i,kb1)
2681 
2682  rcv(i,kb2)=rcv(i,kb1) ; t(i,kb2)=t(i,kb1) ; s(i,kb2)=s(i,kb1)
2683 
2684 
2685  if (k1 <= nz) then ; if (r0(i,k1) >= r0(i,kb1)) then
2686  r0(i,kb2) = 0.5*(r0(i,kb1)+r0(i,k1))
2687 
2688  rcv(i,kb2) = 0.5*(rcv(i,kb1)+rcv(i,k1))
2689  t(i,kb2) = 0.5*(t(i,kb1)+t(i,k1))
2690  s(i,kb2) = 0.5*(s(i,kb1)+s(i,k1))
2691  endif ; endif
2692  endif ! (h2 = 0 && h1 > 0)
2693 
2694  dpe_extrap = 0.0 ; dpe_merge = 0.0
2695  mergeable_bl = .false.
2696  if ((h1 > 0.0) .and. (h2 > 0.0) .and. (h_to_bl > 0.0) .and. &
2697  (stable_rcv)) then
2698  ! Check whether it is permissible for the buffer layers to detrain
2699  ! into the interior isopycnal layers.
2700 
2701  ! Determine the layer that has the lightest target density that is
2702  ! denser than the lowermost buffer layer.
2703  do k1=kb2+1,nz ; if (rcvtgt(k1) >= rcv(i,kb2)) exit ; enddo ; k0 = k1-1
2704  dr1 = rcvtgt(k0)-rcv(i,kb1) ; dr2 = rcv(i,kb2)-rcvtgt(k0)
2705 
2706  ! Use an energy-balanced combination of downwind advection into the next
2707  ! denser interior layer and upwind advection from the upper buffer layer
2708  ! into the lower one, each with an energy change that equals that required
2709  ! to mix the detrained water with the upper buffer layer.
2710  h1_avail = h1 - max(0.0,h_min_bl-h_to_bl)
2711  if ((k1<=nz) .and. (h2 > h_min_bl) .and. (h1_avail > 0.0) .and. &
2712  (r0(i,kb1) < r0(i,kb2)) .and. (h_to_bl*r0(i,kb1) > r0_to_bl)) then
2713  drk1 = (rcvtgt(k1) - rcv(i,kb2)) * (r0(i,kb2) - r0(i,kb1)) / &
2714  (rcv(i,kb2) - rcv(i,kb1))
2715  b1 = drk1 / (r0(i,kb2) - r0(i,kb1))
2716  ! b1 = RcvTgt(k1) - Rcv(i,kb2)) / (Rcv(i,kb2) - Rcv(i,kb1))
2717 
2718  ! Apply several limits to the detrainment.
2719  ! Entrain less than the mass in h2, and keep the base of the buffer
2720  ! layers from becoming shallower than any neighbors.
2721  h2_to_k1 = min(h2 - h_min_bl, h2_to_k1_rem)
2722  ! Balance downwind advection of density into the layer below the
2723  ! buffer layers with upwind advection from the layer above.
2724  if (h2_to_k1*(h1_avail + b1*(h1_avail + h2)) > h2*h1_avail) &
2725  h2_to_k1 = (h2*h1_avail) / (h1_avail + b1*(h1_avail + h2))
2726  if (h2_to_k1*(drk1 * h2) > (h_to_bl*r0(i,kb1) - r0_to_bl) * h1) &
2727  h2_to_k1 = (h_to_bl*r0(i,kb1) - r0_to_bl) * h1 / (drk1 * h2)
2728 
2729  if ((k1==kb2+1) .and. (cs%BL_extrap_lim > 0.)) then
2730  ! Simply do not detrain very light water into the lightest isopycnal
2731  ! coordinate layers if the density jump is too large.
2732  drcv_lim = rcv(i,kb2)-rcv(i,0)
2733  do k=1,kb2 ; drcv_lim = max(drcv_lim, rcv(i,kb2)-rcv(i,k)) ; enddo
2734  drcv_lim = cs%BL_extrap_lim*drcv_lim
2735  if ((rcvtgt(k1) - rcv(i,kb2)) >= drcv_lim) then
2736  h2_to_k1 = 0.0
2737  elseif ((rcvtgt(k1) - rcv(i,kb2)) > 0.5*drcv_lim) then
2738  h2_to_k1 = h2_to_k1 * (2.0 - 2.0*((rcvtgt(k1) - rcv(i,kb2)) / drcv_lim))
2739  endif
2740  endif
2741 
2742  drcv = (rcvtgt(k1) - rcv(i,kb2))
2743 
2744  ! Use 2nd order upwind advection of spiciness, limited by the values
2745  ! in deeper thick layers to determine the detrained temperature and
2746  ! salinity.
2747  dspice_det = (ds_dt_gauge*drcv_ds(i)*(t(i,kb2)-t(i,kb1)) - &
2748  dt_ds_gauge*drcv_dt(i)*(s(i,kb2)-s(i,kb1))) * &
2749  (h2 - h2_to_k1) / (h1 + h2)
2750  dspice_lim = 0.0
2751  if (h(i,k1) > 10.0*angstrom) then
2752  dspice_lim = ds_dt_gauge*drcv_ds(i)*(t(i,k1)-t(i,kb2)) - &
2753  dt_ds_gauge*drcv_dt(i)*(s(i,k1)-s(i,kb2))
2754  if (dspice_det*dspice_lim <= 0.0) dspice_lim = 0.0
2755  endif
2756  if (k1<nz) then ; if (h(i,k1+1) > 10.0*angstrom) then
2757  dspice_lim2 = ds_dt_gauge*drcv_ds(i)*(t(i,k1+1)-t(i,kb2)) - &
2758  dt_ds_gauge*drcv_dt(i)*(s(i,k1+1)-s(i,kb2))
2759  if ((dspice_det*dspice_lim2 > 0.0) .and. &
2760  (abs(dspice_lim2) > abs(dspice_lim))) dspice_lim = dspice_lim2
2761  endif ; endif
2762  if (abs(dspice_det) > abs(dspice_lim)) dspice_det = dspice_lim
2763 
2764  i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
2765  t_det = t(i,kb2) + dt_ds_gauge * i_denom * &
2766  (dt_ds_gauge * drcv_dt(i) * drcv + drcv_ds(i) * dspice_det)
2767  s_det = s(i,kb2) + i_denom * &
2768  (drcv_ds(i) * drcv - dt_ds_gauge * drcv_dt(i) * dspice_det)
2769  ! The detrained values of R0 are based on changes in T and S.
2770  r0_det = r0(i,kb2) + (t_det-t(i,kb2)) * dr0_dt(i) + &
2771  (s_det-s(i,kb2)) * dr0_ds(i)
2772 
2773  if (cs%BL_extrap_lim >= 0.) then
2774  ! Only do this detrainment if the new layer's temperature and salinity
2775  ! are not too far outside of the range of previous values.
2776  if (h(i,k1) > 10.0*angstrom) then
2777  t_min = min(t(i,kb1), t(i,kb2), t(i,k1)) - cs%Allowed_T_chg
2778  t_max = max(t(i,kb1), t(i,kb2), t(i,k1)) + cs%Allowed_T_chg
2779  s_min = min(s(i,kb1), s(i,kb2), s(i,k1)) - cs%Allowed_S_chg
2780  s_max = max(s(i,kb1), s(i,kb2), s(i,k1)) + cs%Allowed_S_chg
2781  else
2782  t_min = min(t(i,kb1), t(i,kb2)) - cs%Allowed_T_chg
2783  t_max = max(t(i,kb1), t(i,kb2)) + cs%Allowed_T_chg
2784  s_min = min(s(i,kb1), s(i,kb2)) - cs%Allowed_S_chg
2785  s_max = max(s(i,kb1), s(i,kb2)) + cs%Allowed_S_chg
2786  endif
2787  ihk1 = 1.0 / (h(i,k1) + h2_to_k1)
2788  t_new = (h(i,k1)*t(i,k1) + h2_to_k1*t_det) * ihk1
2789  s_new = (h(i,k1)*s(i,k1) + h2_to_k1*s_det) * ihk1
2790  ! A less restrictive limit might be used here.
2791  if ((t_new < t_min) .or. (t_new > t_max) .or. &
2792  (s_new < s_min) .or. (s_new > s_max)) &
2793  h2_to_k1 = 0.0
2794  endif
2795 
2796  h1_to_h2 = b1*h2*h2_to_k1 / (h2 - (1.0+b1)*h2_to_k1)
2797 
2798  ihk1 = 1.0 / (h(i,k1) + h_neglect + h2_to_k1)
2799  ih2f = 1.0 / ((h(i,kb2) - h2_to_k1) + h1_to_h2)
2800 
2801  rcv(i,kb2) = ((h(i,kb2)*rcv(i,kb2) - h2_to_k1*rcvtgt(k1)) + &
2802  h1_to_h2*rcv(i,kb1))*ih2f
2803  rcv(i,k1) = ((h(i,k1)+h_neglect)*rcv(i,k1) + h2_to_k1*rcvtgt(k1)) * ihk1
2804 
2805  t(i,kb2) = ((h(i,kb2)*t(i,kb2) - h2_to_k1*t_det) + &
2806  h1_to_h2*t(i,kb1)) * ih2f
2807  t(i,k1) = ((h(i,k1)+h_neglect)*t(i,k1) + h2_to_k1*t_det) * ihk1
2808 
2809  s(i,kb2) = ((h(i,kb2)*s(i,kb2) - h2_to_k1*s_det) + &
2810  h1_to_h2*s(i,kb1)) * ih2f
2811  s(i,k1) = ((h(i,k1)+h_neglect)*s(i,k1) + h2_to_k1*s_det) * ihk1
2812 
2813  ! Changes in R0 are based on changes in T and S.
2814  r0(i,kb2) = ((h(i,kb2)*r0(i,kb2) - h2_to_k1*r0_det) + &
2815  h1_to_h2*r0(i,kb1)) * ih2f
2816  r0(i,k1) = ((h(i,k1)+h_neglect)*r0(i,k1) + h2_to_k1*r0_det) * ihk1
2817 
2818  h(i,kb1) = h(i,kb1) - h1_to_h2 ; h1 = h(i,kb1)
2819  h(i,kb2) = (h(i,kb2) - h2_to_k1) + h1_to_h2 ; h2 = h(i,kb2)
2820  h(i,k1) = h(i,k1) + h2_to_k1
2821 
2822  d_ea(i,kb1) = d_ea(i,kb1) - h1_to_h2
2823  d_ea(i,kb2) = (d_ea(i,kb2) - h2_to_k1) + h1_to_h2
2824  d_ea(i,k1) = d_ea(i,k1) + h2_to_k1
2825  h2_to_k1_rem = max(h2_to_k1_rem - h2_to_k1, 0.0)
2826 
2827  ! The lower buffer layer has become lighter - it may be necessary to
2828  ! adjust k1 lighter.
2829  if ((k1>kb2+1) .and. (rcvtgt(k1-1) >= rcv(i,kb2))) then
2830  do k1=k1,kb2+1,-1 ; if (rcvtgt(k1-1) < rcv(i,kb2)) exit ; enddo
2831  endif
2832  endif
2833 
2834  k0 = k1-1
2835  dr1 = rcvtgt(k0)-rcv(i,kb1) ; dr2 = rcv(i,kb2)-rcvtgt(k0)
2836 
2837  if ((k0>kb2) .and. (dr1 > 0.0) .and. (h1 > h_min_bl) .and. &
2838  (h2*dr2 < h1*dr1) .and. (r0(i,kb2) > r0(i,kb1))) then
2839  ! An interior isopycnal layer (k0) is intermediate in density between
2840  ! the two buffer layers, and there can be detrainment. The entire
2841  ! lower buffer layer is combined with a portion of the upper buffer
2842  ! layer to match the target density of layer k0.
2843  stays_merge = 2.0*(h1+h2)*(h1*dr1 - h2*dr2) / &
2844  ((dr1+dr2)*h1 + dr1*(h1+h2) + &
2845  sqrt((dr2*h1-dr1*h2)**2 + 4*(h1+h2)*h2*(dr1+dr2)*dr2))
2846 
2847  stays_min_merge = max(h_min_bl, 2.0*h_min_bl - h_to_bl, &
2848  h1 - (h1+h2)*(r0(i,kb1) - r0_det) / (r0(i,kb2) - r0(i,kb1)))
2849  if ((stays_merge > stays_min_merge) .and. &
2850  (stays_merge + h2_to_k1_rem >= h1 + h2)) then
2851  mergeable_bl = .true.
2852  dpe_merge = g_2*(r0(i,kb2)-r0(i,kb1))*(h1-stays_merge)*(h2-stays_merge)
2853  endif
2854  endif
2855 
2856  if ((k1<=nz).and.(.not.mergeable_bl)) then
2857  ! Check whether linear extrapolation of density (i.e. 2nd order upwind
2858  ! advection) will allow some of the lower buffer layer to detrain into
2859  ! the next denser interior layer (k1).
2860  dr2b = rcvtgt(k1)-rcv(i,kb2) ; dr21 = rcv(i,kb2) - rcv(i,kb1)
2861  if (dr2b*(h1+h2) < h2*dr21) then
2862  ! Some of layer kb2 is denser than k1.
2863  h2_to_k1 = min(h2 - (h1+h2) * dr2b / dr21, h2_to_k1_rem)
2864 
2865  if (h2 > h2_to_k1) then
2866  drcv = (rcvtgt(k1) - rcv(i,kb2))
2867 
2868  ! Use 2nd order upwind advection of spiciness, limited by the values
2869  ! in deeper thick layers to determine the detrained temperature and
2870  ! salinity.
2871  dspice_det = (ds_dt_gauge*drcv_ds(i)*(t(i,kb2)-t(i,kb1)) - &
2872  dt_ds_gauge*drcv_dt(i)*(s(i,kb2)-s(i,kb1))) * &
2873  (h2 - h2_to_k1) / (h1 + h2)
2874  dspice_lim = 0.0
2875  if (h(i,k1) > 10.0*angstrom) then
2876  dspice_lim = ds_dt_gauge*drcv_ds(i)*(t(i,k1)-t(i,kb2)) - &
2877  dt_ds_gauge*drcv_dt(i)*(s(i,k1)-s(i,kb2))
2878  if (dspice_det*dspice_lim <= 0.0) dspice_lim = 0.0
2879  endif
2880  if (k1<nz) then; if (h(i,k1+1) > 10.0*angstrom) then
2881  dspice_lim2 = ds_dt_gauge*drcv_ds(i)*(t(i,k1+1)-t(i,kb2)) - &
2882  dt_ds_gauge*drcv_dt(i)*(s(i,k1+1)-s(i,kb2))
2883  if ((dspice_det*dspice_lim2 > 0.0) .and. &
2884  (abs(dspice_lim2) > abs(dspice_lim))) dspice_lim = dspice_lim2
2885  endif; endif
2886  if (abs(dspice_det) > abs(dspice_lim)) dspice_det = dspice_lim
2887 
2888  i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
2889  t_det = t(i,kb2) + dt_ds_gauge * i_denom * &
2890  (dt_ds_gauge * drcv_dt(i) * drcv + drcv_ds(i) * dspice_det)
2891  s_det = s(i,kb2) + i_denom * &
2892  (drcv_ds(i) * drcv - dt_ds_gauge * drcv_dt(i) * dspice_det)
2893  ! The detrained values of R0 are based on changes in T and S.
2894  r0_det = r0(i,kb2) + (t_det-t(i,kb2)) * dr0_dt(i) + &
2895  (s_det-s(i,kb2)) * dr0_ds(i)
2896 
2897  ! Now that the properties of the detrained water are known,
2898  ! potentially limit the amount of water that is detrained to
2899  ! avoid creating unphysical properties in the remaining water.
2900  ih2f = 1.0 / (h2 - h2_to_k1)
2901 
2902  t_min = min(t(i,kb2), t(i,kb1)) - cs%Allowed_T_chg
2903  t_max = max(t(i,kb2), t(i,kb1)) + cs%Allowed_T_chg
2904  t_new = (h2*t(i,kb2) - h2_to_k1*t_det)*ih2f
2905  if (t_new < t_min) then
2906  h2_to_k1_lim = h2 * (t(i,kb2) - t_min) / (t_det - t_min)
2907 ! write(mesg,'("Low temperature limits det to ", &
2908 ! & 1pe12.5, " from ", 1pe12.5, " at ", 1pg11.4,"E, ",1pg11.4,"N. T=", &
2909 ! & 5(1pe12.5))') &
2910 ! h2_to_k1_lim, h2_to_k1, G%geoLonT(i,j), G%geoLatT(i,j), &
2911 ! T_new, T(i,kb2), T(i,kb1), T_det, T_new-T_min
2912 ! call MOM_error(WARNING, mesg)
2913  h2_to_k1 = h2_to_k1_lim
2914  ih2f = 1.0 / (h2 - h2_to_k1)
2915  elseif (t_new > t_max) then
2916  h2_to_k1_lim = h2 * (t(i,kb2) - t_max) / (t_det - t_max)
2917 ! write(mesg,'("High temperature limits det to ", &
2918 ! & 1pe12.5, " from ", 1pe12.5, " at ", 1pg11.4,"E, ",1pg11.4,"N. T=", &
2919 ! & 5(1pe12.5))') &
2920 ! h2_to_k1_lim, h2_to_k1, G%geoLonT(i,j), G%geoLatT(i,j), &
2921 ! T_new, T(i,kb2), T(i,kb1), T_det, T_new-T_max
2922 ! call MOM_error(WARNING, mesg)
2923  h2_to_k1 = h2_to_k1_lim
2924  ih2f = 1.0 / (h2 - h2_to_k1)
2925  endif
2926  s_min = max(min(s(i,kb2), s(i,kb1)) - cs%Allowed_S_chg, 0.0)
2927  s_max = max(s(i,kb2), s(i,kb1)) + cs%Allowed_S_chg
2928  s_new = (h2*s(i,kb2) - h2_to_k1*s_det)*ih2f
2929  if (s_new < s_min) then
2930  h2_to_k1_lim = h2 * (s(i,kb2) - s_min) / (s_det - s_min)
2931 ! write(mesg,'("Low salinity limits det to ", &
2932 ! & 1pe12.5, " from ", 1pe12.5, " at ", 1pg11.4,"E, ",1pg11.4,"N. S=", &
2933 ! & 5(1pe12.5))') &
2934 ! h2_to_k1_lim, h2_to_k1, G%geoLonT(i,j), G%geoLatT(i,j), &
2935 ! S_new, S(i,kb2), S(i,kb1), S_det, S_new-S_min
2936 ! call MOM_error(WARNING, mesg)
2937  h2_to_k1 = h2_to_k1_lim
2938  ih2f = 1.0 / (h2 - h2_to_k1)
2939  elseif (s_new > s_max) then
2940  h2_to_k1_lim = h2 * (s(i,kb2) - s_max) / (s_det - s_max)
2941 ! write(mesg,'("High salinity limits det to ", &
2942 ! & 1pe12.5, " from ", 1pe12.5, " at ", 1pg11.4,"E, ",1pg11.4,"N. S=", &
2943 ! & 5(1pe12.5))') &
2944 ! h2_to_k1_lim, h2_to_k1, G%geoLonT(i,j), G%geoLatT(i,j), &
2945 ! S_new, S(i,kb2), S(i,kb1), S_det, S_new-S_max
2946 ! call MOM_error(WARNING, mesg)
2947  h2_to_k1 = h2_to_k1_lim
2948  ih2f = 1.0 / (h2 - h2_to_k1)
2949  endif
2950 
2951  ihk1 = 1.0 / (h(i,k1) + h_neglect + h2_to_k1)
2952  rcv(i,k1) = ((h(i,k1)+h_neglect)*rcv(i,k1) + h2_to_k1*rcvtgt(k1)) * ihk1
2953  rcv(i,kb2) = rcv(i,kb2) - h2_to_k1*drcv*ih2f
2954 
2955  t(i,kb2) = (h2*t(i,kb2) - h2_to_k1*t_det)*ih2f
2956  t(i,k1) = ((h(i,k1)+h_neglect)*t(i,k1) + h2_to_k1*t_det) * ihk1
2957 
2958  s(i,kb2) = (h2*s(i,kb2) - h2_to_k1*s_det) * ih2f
2959  s(i,k1) = ((h(i,k1)+h_neglect)*s(i,k1) + h2_to_k1*s_det) * ihk1
2960 
2961  ! Changes in R0 are based on changes in T and S.
2962  r0(i,kb2) = (h2*r0(i,kb2) - h2_to_k1*r0_det) * ih2f
2963  r0(i,k1) = ((h(i,k1)+h_neglect)*r0(i,k1) + h2_to_k1*r0_det) * ihk1
2964  else
2965  ! h2==h2_to_k1 can happen if dR2b = 0 exactly, but this is very
2966  ! unlikely. In this case the entirety of layer kb2 is detrained.
2967  h2_to_k1 = h2 ! These 2 lines are probably unnecessary.
2968  ihk1 = 1.0 / (h(i,k1) + h2)
2969 
2970  rcv(i,k1) = (h(i,k1)*rcv(i,k1) + h2*rcv(i,kb2)) * ihk1
2971  t(i,k1) = (h(i,k1)*t(i,k1) + h2*t(i,kb2)) * ihk1
2972  s(i,k1) = (h(i,k1)*s(i,k1) + h2*s(i,kb2)) * ihk1
2973  r0(i,k1) = (h(i,k1)*r0(i,k1) + h2*r0(i,kb2)) * ihk1
2974  endif
2975 
2976  h(i,k1) = h(i,k1) + h2_to_k1
2977  h(i,kb2) = h(i,kb2) - h2_to_k1 ; h2 = h(i,kb2)
2978  ! dPE_extrap should be positive here.
2979  dpe_extrap = i2rho0*(r0_det-r0(i,kb2))*h2_to_k1*h2
2980 
2981  d_ea(i,kb2) = d_ea(i,kb2) - h2_to_k1
2982  d_ea(i,k1) = d_ea(i,k1) + h2_to_k1
2983  h2_to_k1_rem = max(h2_to_k1_rem - h2_to_k1, 0.0)
2984  endif
2985  endif ! Detrainment by extrapolation.
2986 
2987  endif ! Detrainment to the interior at all.
2988 
2989  ! Does some of the detrained water go into the lower buffer layer?
2990  h_det_h2 = max(h_min_bl-(h1+h2), 0.0)
2991  if (h_det_h2 > 0.0) then
2992  ! Detrained water will go into both upper and lower buffer layers.
2993  ! h(kb2) will be h_min_bl, but h(kb1) may be larger if there was already
2994  ! ample detrainment; all water in layer kb1 moves into layer kb2.
2995 
2996  ! Determine the fluxes between the various layers.
2997  h_det_to_h2 = min(h_to_bl, h_det_h2)
2998  h_ml_to_h2 = h_det_h2 - h_det_to_h2
2999  h_det_to_h1 = h_to_bl - h_det_to_h2
3000  h_ml_to_h1 = max(h_min_bl-h_det_to_h1,0.0)
3001 
3002  ih = 1.0/h_min_bl;
3003  ihdet = 0.0 ; if (h_to_bl > 0.0) ihdet = 1.0 / h_to_bl
3004  ih1f = 1.0 / (h_det_to_h1 + h_ml_to_h1)
3005 
3006  r0(i,kb2) = ((h2*r0(i,kb2) + h1*r0(i,kb1)) + &
3007  (h_det_to_h2*r0_to_bl*ihdet + h_ml_to_h2*r0(i,0))) * ih
3008  r0(i,kb1) = (h_det_to_h1*r0_to_bl*ihdet + h_ml_to_h1*r0(i,0)) * ih1f
3009 
3010  rcv(i,kb2) = ((h2*rcv(i,kb2) + h1*rcv(i,kb1)) + &
3011  (h_det_to_h2*rcv_to_bl*ihdet + h_ml_to_h2*rcv(i,0))) * ih
3012  rcv(i,kb1) = (h_det_to_h1*rcv_to_bl*ihdet + h_ml_to_h1*rcv(i,0)) * ih1f
3013 
3014  t(i,kb2) = ((h2*t(i,kb2) + h1*t(i,kb1)) + &
3015  (h_det_to_h2*t_to_bl*ihdet + h_ml_to_h2*t(i,0))) * ih
3016  t(i,kb1) = (h_det_to_h1*t_to_bl*ihdet + h_ml_to_h1*t(i,0)) * ih1f
3017 
3018  s(i,kb2) = ((h2*s(i,kb2) + h1*s(i,kb1)) + &
3019  (h_det_to_h2*s_to_bl*ihdet + h_ml_to_h2*s(i,0))) * ih
3020  s(i,kb1) = (h_det_to_h1*s_to_bl*ihdet + h_ml_to_h1*s(i,0)) * ih1f
3021 
3022  ! Recall that h1 = h(i,kb1) & h2 = h(i,kb2).
3023  d_ea(i,1) = d_ea(i,1) - (h_ml_to_h1 + h_ml_to_h2)
3024  d_ea(i,kb1) = d_ea(i,kb1) + ((h_det_to_h1 + h_ml_to_h1) - h1)
3025  d_ea(i,kb2) = d_ea(i,kb2) + (h_min_bl - h2)
3026 
3027  h(i,kb1) = h_det_to_h1 + h_ml_to_h1 ; h(i,kb2) = h_min_bl
3028  h(i,0) = h(i,0) - (h_ml_to_h1 + h_ml_to_h2)
3029 
3030 
3031  if (ALLOCATED(cs%diag_PE_detrain) .or. ALLOCATED(cs%diag_PE_detrain2)) then
3032  r0_det = r0_to_bl*ihdet
3033  s1en = g_2 * idt_h2 * ( ((r0(i,kb2)-r0(i,kb1))*h1*h2 + &
3034  h_det_to_h2*( (r0(i,kb1)-r0_det)*h1 + (r0(i,kb2)-r0_det)*h2 ) + &
3035  h_ml_to_h2*( (r0(i,kb2)-r0(i,0))*h2 + (r0(i,kb1)-r0(i,0))*h1 + &
3036  (r0_det-r0(i,0))*h_det_to_h2 ) + &
3037  h_det_to_h1*h_ml_to_h1*(r0_det-r0(i,0))) - 2.0*gv%Rho0*dpe_extrap )
3038 
3039  if (ALLOCATED(cs%diag_PE_detrain)) &
3040  cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + s1en
3041 
3042  if (ALLOCATED(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
3043  cs%diag_PE_detrain2(i,j) + s1en + idt_h2*rho0xg*dpe_extrap
3044  endif
3045 
3046  elseif ((h_to_bl > 0.0) .or. (h1 < h_min_bl) .or. (h2 < h_min_bl)) then
3047  ! Determine how much of the upper buffer layer will be moved into
3048  ! the lower buffer layer and the properties with which it is moving.
3049  ! This implementation assumes a 2nd-order upwind advection of density
3050  ! from the uppermost buffer layer into the next one down.
3051  h_from_ml = h_min_bl + max(h_min_bl-h2,0.0) - h1 - h_to_bl
3052  if (h_from_ml > 0.0) then
3053  ! Some water needs to be moved from the mixed layer so that the upper
3054  ! (and perhaps lower) buffer layers exceed their minimum thicknesses.
3055  dpe_extrap = dpe_extrap - i2rho0*h_from_ml*(r0_to_bl - r0(i,0)*h_to_bl)
3056  r0_to_bl = r0_to_bl + h_from_ml*r0(i,0)
3057  rcv_to_bl = rcv_to_bl + h_from_ml*rcv(i,0)
3058  t_to_bl = t_to_bl + h_from_ml*t(i,0)
3059  s_to_bl = s_to_bl + h_from_ml*s(i,0)
3060 
3061  h_to_bl = h_to_bl + h_from_ml
3062  h(i,0) = h(i,0) - h_from_ml
3063  d_ea(i,1) = d_ea(i,1) - h_from_ml
3064  endif
3065 
3066  ! The absolute value should be unnecessary and 1e9 is just a large number.
3067  b1 = 1.0e9
3068  if (r0(i,kb2) - r0(i,kb1) > 1.0e-9*abs(r0(i,kb1) - r0_det)) &
3069  b1 = abs(r0(i,kb1) - r0_det) / (r0(i,kb2) - r0(i,kb1))
3070  stays_min = max((1.0-b1)*h1 - b1*h2, 0.0, h_min_bl - h_to_bl)
3071  stays_max = h1 - max(h_min_bl-h2,0.0)
3072 
3073  scale_slope = 1.0
3074  if (stays_max <= stays_min) then
3075  stays = stays_max
3076  mergeable_bl = .false.
3077  if (stays_max < h1) scale_slope = (h1 - stays_min) / (h1 - stays_max)
3078  else
3079  ! There are numerous temporary variables used here that should not be
3080  ! used outside of this "else" branch: s1, s2, s3sq, I_ya, bh0
3081  bh0 = b1*h_to_bl
3082  i_ya = (h1 + h2) / ((h1 + h2) + h_to_bl)
3083  ! s1 is the amount staying that minimizes the PE increase.
3084  s1 = 0.5*(h1 + (h2 - bh0) * i_ya) ; s2 = h1 - s1
3085 
3086  if (s2 < 0.0) then
3087  ! The energy released by detrainment from the lower buffer layer can be
3088  ! used to mix water from the upper buffer layer into the lower one.
3089  s3sq = i_ya*max(bh0*h1-dpe_extrap, 0.0)
3090  else
3091  s3sq = i_ya*(bh0*h1-min(dpe_extrap,0.0))
3092  endif
3093 
3094  if (s3sq == 0.0) then
3095  ! There is a simple, exact solution to the quadratic equation, namely:
3096  stays = h1 ! This will revert to stays_max later.
3097  elseif (s2*s2 <= s3sq) then
3098  ! There is no solution with 0 PE change - use the minimum energy input.
3099  stays = s1
3100  else
3101  ! The following choose the solutions that are continuous with all water
3102  ! staying in the upper buffer layer when there is no detrainment,
3103  ! namely the + root when s2>0 and the - root otherwise. They also
3104  ! carefully avoid differencing large numbers, using s2 = (h1-s).
3105  if (bh0 <= 0.0) then ; stays = h1
3106  elseif (s2 > 0.0) then
3107  ! stays = s + sqrt(s2*s2 - s3sq) ! Note that s2 = h1-s
3108  if (s1 >= stays_max) then ; stays = stays_max
3109  elseif (s1 >= 0.0) then ; stays = s1 + sqrt(s2*s2 - s3sq)
3110  else ; stays = (h1*(s2-s1) - s3sq) / (-s1 + sqrt(s2*s2 - s3sq))
3111  endif
3112  else
3113  ! stays = s - sqrt(s2*s2 - s3sq) ! Note that s2 = h1-s & stays_min >= 0
3114  if (s1 <= stays_min) then ; stays = stays_min
3115  else ; stays = (h1*(s1-s2) + s3sq) / (s1 + sqrt(s2*s2 - s3sq))
3116  endif
3117  endif
3118  endif
3119 
3120  ! Limit the amount that stays so that the motion of water is from the
3121  ! upper buffer layer into the lower, but no more than is in the upper
3122  ! layer, and the water left in the upper layer is no lighter than the
3123  ! detrained water.
3124  if (stays >= stays_max) then ; stays = stays_max
3125  elseif (stays < stays_min) then ; stays = stays_min
3126  endif
3127  endif
3128 
3129  dpe_det = g_2*((r0(i,kb1)*h_to_bl - r0_to_bl)*stays + &
3130  (r0(i,kb2)-r0(i,kb1)) * (h1-stays) * &
3131  (h2 - scale_slope*stays*((h1+h2)+h_to_bl)/(h1+h2)) ) - &
3132  rho0xg*dpe_extrap
3133 
3134  if (dpe_time_ratio*h_to_bl > h_to_bl+h(i,0)) then
3135  dpe_ratio = (h_to_bl+h(i,0)) / h_to_bl
3136  else
3137  dpe_ratio = dpe_time_ratio
3138  endif
3139 
3140  if ((mergeable_bl) .and. (num_events*dpe_ratio*dpe_det > dpe_merge)) then
3141  ! It is energetically preferable to merge the two buffer layers, detrain
3142  ! them into interior layer (k0), move the remaining upper buffer layer
3143  ! water into the lower buffer layer, and detrain undiluted into the
3144  ! upper buffer layer.
3145  h1_to_k0 = (h1-stays_merge)
3146  stays = max(h_min_bl-h_to_bl,0.0)
3147  h1_to_h2 = stays_merge - stays
3148 
3149  ihk0 = 1.0 / ((h1_to_k0 + h2) + h(i,k0))
3150  ih1f = 1.0 / (h_to_bl + stays); ih2f = 1.0 / h1_to_h2
3151  ih12 = 1.0 / (h1 + h2)
3152 
3153  drcv_2dz = (rcv(i,kb1) - rcv(i,kb2)) * ih12
3154  drcv_stays = drcv_2dz*(h1_to_k0 + h1_to_h2)
3155  drcv_det = - drcv_2dz*(stays + h1_to_h2)
3156  rcv(i,k0) = ((h1_to_k0*(rcv(i,kb1) + drcv_det) + &
3157  h2*rcv(i,kb2)) + h(i,k0)*rcv(i,k0)) * ihk0
3158  rcv(i,kb2) = rcv(i,kb1) + drcv_2dz*(h1_to_k0-stays)
3159  rcv(i,kb1) = (rcv_to_bl + stays*(rcv(i,kb1) + drcv_stays)) * ih1f
3160 
3161  ! Use 2nd order upwind advection of spiciness, limited by the value in
3162  ! the water from the mixed layer to determine the temperature and
3163  ! salinity of the water that stays in the buffer layers.
3164  i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
3165  dspice_2dz = (ds_dt_gauge*drcv_ds(i)*(t(i,kb1)-t(i,kb2)) - &
3166  dt_ds_gauge*drcv_dt(i)*(s(i,kb1)-s(i,kb2))) * ih12
3167  dspice_lim = (ds_dt_gauge*dr0_ds(i)*(t_to_bl-t(i,kb1)*h_to_bl) - &
3168  dt_ds_gauge*dr0_dt(i)*(s_to_bl-s(i,kb1)*h_to_bl)) / h_to_bl
3169  if (dspice_lim * dspice_2dz <= 0.0) dspice_2dz = 0.0
3170 
3171  if (stays > 0.0) then
3172  ! Limit the spiciness of the water that stays in the upper buffer layer.
3173  if (abs(dspice_lim) < abs(dspice_2dz*(h1_to_k0 + h1_to_h2))) &
3174  dspice_2dz = dspice_lim/(h1_to_k0 + h1_to_h2)
3175 
3176  dspice_stays = dspice_2dz*(h1_to_k0 + h1_to_h2)
3177  t_stays = t(i,kb1) + dt_ds_gauge * i_denom * &
3178  (dt_ds_gauge * drcv_dt(i) * drcv_stays + drcv_ds(i) * dspice_stays)
3179  s_stays = s(i,kb1) + i_denom * &
3180  (drcv_ds(i) * drcv_stays - dt_ds_gauge * drcv_dt(i) * dspice_stays)
3181  ! The values of R0 are based on changes in T and S.
3182  r0_stays = r0(i,kb1) + (t_stays-t(i,kb1)) * dr0_dt(i) + &
3183  (s_stays-s(i,kb1)) * dr0_ds(i)
3184  else
3185  ! Limit the spiciness of the water that moves into the lower buffer layer.
3186  if (abs(dspice_lim) < abs(dspice_2dz*h1_to_k0)) &
3187  dspice_2dz = dspice_lim/h1_to_k0
3188  ! These will be multiplied by 0 later.
3189  t_stays = 0.0 ; s_stays = 0.0 ; r0_stays = 0.0
3190  endif
3191 
3192  dspice_det = - dspice_2dz*(stays + h1_to_h2)
3193  t_det = t(i,kb1) + dt_ds_gauge * i_denom * &
3194  (dt_ds_gauge * drcv_dt(i) * drcv_det + drcv_ds(i) * dspice_det)
3195  s_det = s(i,kb1) + i_denom * &
3196  (drcv_ds(i) * drcv_det - dt_ds_gauge * drcv_dt(i) * dspice_det)
3197  ! The values of R0 are based on changes in T and S.
3198  r0_det = r0(i,kb1) + (t_det-t(i,kb1)) * dr0_dt(i) + &
3199  (s_det-s(i,kb1)) * dr0_ds(i)
3200 
3201  t(i,k0) = ((h1_to_k0*t_det + h2*t(i,kb2)) + h(i,k0)*t(i,k0)) * ihk0
3202  t(i,kb2) = (h1*t(i,kb1) - stays*t_stays - h1_to_k0*t_det) * ih2f
3203  t(i,kb1) = (t_to_bl + stays*t_stays) * ih1f
3204 
3205  s(i,k0) = ((h1_to_k0*s_det + h2*s(i,kb2)) + h(i,k0)*s(i,k0)) * ihk0
3206  s(i,kb2) = (h1*s(i,kb1) - stays*s_stays - h1_to_k0*s_det) * ih2f
3207  s(i,kb1) = (s_to_bl + stays*s_stays) * ih1f
3208 
3209  r0(i,k0) = ((h1_to_k0*r0_det + h2*r0(i,kb2)) + h(i,k0)*r0(i,k0)) * ihk0
3210  r0(i,kb2) = (h1*r0(i,kb1) - stays*r0_stays - h1_to_k0*r0_det) * ih2f
3211  r0(i,kb1) = (r0_to_bl + stays*r0_stays) * ih1f
3212 
3213 ! ! The following is 2nd-order upwind advection without limiters.
3214 ! dT_2dz = (T(i,kb1) - T(i,kb2)) * Ih12
3215 ! T(i,k0) = (h1_to_k0*(T(i,kb1) - dT_2dz*(stays+h1_to_h2)) + &
3216 ! h2*T(i,kb2) + h(i,k0)*T(i,k0)) * Ihk0
3217 ! T(i,kb2) = T(i,kb1) + dT_2dz*(h1_to_k0-stays)
3218 ! T(i,kb1) = (T_to_bl + stays*(T(i,kb1) + &
3219 ! dT_2dz*(h1_to_k0 + h1_to_h2))) * Ih1f
3220 ! dS_2dz = (S(i,kb1) - S(i,kb2)) * Ih12
3221 ! S(i,k0) = (h1_to_k0*(S(i,kb1) - dS_2dz*(stays+h1_to_h2)) + &
3222 ! h2*S(i,kb2) + h(i,k0)*S(i,k0)) * Ihk0
3223 ! S(i,kb2) = S(i,kb1) + dS_2dz*(h1_to_k0-stays)
3224 ! S(i,kb1) = (S_to_bl + stays*(S(i,kb1) + &
3225 ! dS_2dz*(h1_to_k0 + h1_to_h2))) * Ih1f
3226 ! dR0_2dz = (R0(i,kb1) - R0(i,kb2)) * Ih12
3227 ! R0(i,k0) = (h1_to_k0*(R0(i,kb1) - dR0_2dz*(stays+h1_to_h2)) + &
3228 ! h2*R0(i,kb2) + h(i,k0)*R0(i,k0)) * Ihk0
3229 ! R0(i,kb2) = R0(i,kb1) + dR0_2dz*(h1_to_k0-stays)
3230 ! R0(i,kb1) = (R0_to_bl + stays*(R0(i,kb1) + &
3231 ! dR0_2dz*(h1_to_k0 + h1_to_h2))) * Ih1f
3232 
3233  d_ea(i,kb1) = (d_ea(i,kb1) + h_to_bl) + (stays - h1)
3234  d_ea(i,kb2) = d_ea(i,kb2) + (h1_to_h2 - h2)
3235  d_ea(i,k0) = d_ea(i,k0) + (h1_to_k0 + h2)
3236 
3237  h(i,kb1) = stays + h_to_bl
3238  h(i,kb2) = h1_to_h2
3239  h(i,k0) = h(i,k0) + (h1_to_k0 + h2)
3240  if (ALLOCATED(cs%diag_PE_detrain)) &
3241  cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + idt_h2*dpe_merge
3242  if (ALLOCATED(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
3243  cs%diag_PE_detrain2(i,j) + idt_h2*(dpe_det+rho0xg*dpe_extrap)
3244  else ! Not mergeable_bl.
3245  ! There is no further detrainment from the buffer layers, and the
3246  ! upper buffer layer water is distributed optimally between the
3247  ! upper and lower buffer layer.
3248  h1_to_h2 = h1 - stays
3249  ih1f = 1.0 / (h_to_bl + stays) ; ih2f = 1.0 / (h2 + h1_to_h2)
3250  ih = 1.0 / (h1 + h2)
3251  dr0_2dz = (r0(i,kb1) - r0(i,kb2)) * ih
3252  r0(i,kb2) = (h2*r0(i,kb2) + h1_to_h2*(r0(i,kb1) - &
3253  scale_slope*dr0_2dz*stays)) * ih2f
3254  r0(i,kb1) = (r0_to_bl + stays*(r0(i,kb1) + &
3255  scale_slope*dr0_2dz*h1_to_h2)) * ih1f
3256 
3257  ! Use 2nd order upwind advection of spiciness, limited by the value
3258  ! in the detrained water to determine the detrained temperature and
3259  ! salinity.
3260  dr0 = scale_slope*dr0_2dz*h1_to_h2
3261  dspice_stays = (ds_dt_gauge*dr0_ds(i)*(t(i,kb1)-t(i,kb2)) - &
3262  dt_ds_gauge*dr0_dt(i)*(s(i,kb1)-s(i,kb2))) * &
3263  scale_slope*h1_to_h2 * ih
3264  if (h_to_bl > 0.0) then
3265  dspice_lim = (ds_dt_gauge*dr0_ds(i)*(t_to_bl-t(i,kb1)*h_to_bl) - &
3266  dt_ds_gauge*dr0_dt(i)*(s_to_bl-s(i,kb1)*h_to_bl)) /&
3267  h_to_bl
3268  else
3269  dspice_lim = ds_dt_gauge*dr0_ds(i)*(t(i,0)-t(i,kb1)) - &
3270  dt_ds_gauge*dr0_dt(i)*(s(i,0)-s(i,kb1))
3271  endif
3272  if (dspice_stays*dspice_lim <= 0.0) then
3273  dspice_stays = 0.0
3274  elseif (abs(dspice_stays) > abs(dspice_lim)) then
3275  dspice_stays = dspice_lim
3276  endif
3277  i_denom = 1.0 / (dr0_ds(i)**2 + (dt_ds_gauge*dr0_dt(i))**2)
3278  t_stays = t(i,kb1) + dt_ds_gauge * i_denom * &
3279  (dt_ds_gauge * dr0_dt(i) * dr0 + dr0_ds(i) * dspice_stays)
3280  s_stays = s(i,kb1) + i_denom * &
3281  (dr0_ds(i) * dr0 - dt_ds_gauge * dr0_dt(i) * dspice_stays)
3282  ! The detrained values of Rcv are based on changes in T and S.
3283  rcv_stays = rcv(i,kb1) + (t_stays-t(i,kb1)) * drcv_dt(i) + &
3284  (s_stays-s(i,kb1)) * drcv_ds(i)
3285 
3286  t(i,kb2) = (h2*t(i,kb2) + h1*t(i,kb1) - t_stays*stays) * ih2f
3287  t(i,kb1) = (t_to_bl + stays*t_stays) * ih1f
3288  s(i,kb2) = (h2*s(i,kb2) + h1*s(i,kb1) - s_stays*stays) * ih2f
3289  s(i,kb1) = (s_to_bl + stays*s_stays) * ih1f
3290  rcv(i,kb2) = (h2*rcv(i,kb2) + h1*rcv(i,kb1) - rcv_stays*stays) * ih2f
3291  rcv(i,kb1) = (rcv_to_bl + stays*rcv_stays) * ih1f
3292 
3293 ! ! The following is 2nd-order upwind advection without limiters.
3294 ! dRcv_2dz = (Rcv(i,kb1) - Rcv(i,kb2)) * Ih
3295 ! dRcv = scale_slope*dRcv_2dz*h1_to_h2
3296 ! Rcv(i,kb2) = (h2*Rcv(i,kb2) + h1_to_h2*(Rcv(i,kb1) - &
3297 ! scale_slope*dRcv_2dz*stays)) * Ih2f
3298 ! Rcv(i,kb1) = (Rcv_to_bl + stays*(Rcv(i,kb1) + dRcv)) * Ih1f
3299 ! dT_2dz = (T(i,kb1) - T(i,kb2)) * Ih
3300 ! T(i,kb2) = (h2*T(i,kb2) + h1_to_h2*(T(i,kb1) - &
3301 ! scale_slope*dT_2dz*stays)) * Ih2f
3302 ! T(i,kb1) = (T_to_bl + stays*(T(i,kb1) + &
3303 ! scale_slope*dT_2dz*h1_to_h2)) * Ih1f
3304 ! dS_2dz = (S(i,kb1) - S(i,kb2)) * Ih
3305 ! S(i,kb2) = (h2*S(i,kb2) + h1_to_h2*(S(i,kb1) - &
3306 ! scale_slope*dS_2dz*stays)) * Ih2f
3307 ! S(i,kb1) = (S_to_bl + stays*(S(i,kb1) + &
3308 ! scale_slope*dS_2dz*h1_to_h2)) * Ih1f
3309 
3310  d_ea(i,kb1) = d_ea(i,kb1) + ((stays - h1) + h_to_bl)
3311  d_ea(i,kb2) = d_ea(i,kb2) + h1_to_h2
3312 
3313  h(i,kb1) = stays + h_to_bl
3314  h(i,kb2) = h(i,kb2) + h1_to_h2
3315 
3316  if (ALLOCATED(cs%diag_PE_detrain)) &
3317  cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + idt_h2*dpe_det
3318  if (ALLOCATED(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
3319  cs%diag_PE_detrain2(i,j) + idt_h2*(dpe_det+rho0xg*dpe_extrap)
3320  endif
3321  endif ! End of detrainment...
3322 
3323  enddo ! i loop
3324 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ resort_ml()

subroutine mom_bulk_mixed_layer::resort_ml ( real, dimension(szi_(g),szk0_(gv)), intent(inout)  h,
real, dimension(szi_(g),szk0_(gv)), intent(inout)  T,
real, dimension(szi_(g),szk0_(gv)), intent(inout)  S,
real, dimension(szi_(g),szk0_(gv)), intent(inout)  R0,
real, dimension(szi_(g),szk0_(gv)), intent(inout)  Rcv,
real, dimension(szk_(gv)), intent(in)  RcvTgt,
real, dimension(szi_(g),szk_(gv)), intent(inout)  eps,
real, dimension(szi_(g),szk_(gv)), intent(inout)  d_ea,
real, dimension(szi_(g),szk_(gv)), intent(inout)  d_eb,
integer, dimension(szi_(g),szk_(gv)), intent(in)  ksort,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(bulkmixedlayer_cs), pointer  CS,
real, dimension(szi_(g)), intent(in)  dR0_dT,
real, dimension(szi_(g)), intent(in)  dR0_dS,
real, dimension(szi_(g)), intent(in)  dRcv_dT,
real, dimension(szi_(g)), intent(in)  dRcv_dS 
)
private

This subroutine actually moves properties between layers to achieve a resorted state, with all of the resorted water either moved into the correct interior layers or in the top nkmb layers.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in,out]hLayer thickness, in m or kg m-2. (Intent in/out) The units of h are referred to as H below. Layer 0 is the new mixed layer.
[in,out]tLayer temperatures, in deg C.
[in,out]sLayer salinities, in psu.
[in,out]r0Potential density referenced to surface pressure, in kg m-3.
[in,out]rcvThe coordinate defining potential density, in kg m-3.
[in]rcvtgtThe target value of Rcv for each layer, in kg m-3.
[in,out]epsThe (small) thickness that must remain in each layer, in H.
[in,out]d_eaThe upward increase across a layer in the entrainment from above, in m or kg m-2 (H). Positive d_ea goes with layer thickness increases.
[in,out]d_ebThe downward increase across a layer in the entrainment from below, in H. Positive values go with mass gain by a layer.
[in]ksortThe density-sorted k-indicies.
csThe control structure for this module.
[in]dr0_dtThe partial derivative of potential density referenced to the surface with potential temperature, in kg m-3 K-1.
[in]dr0_dsThe partial derivative of cpotential density referenced to the surface with salinity, in kg m-3 psu-1.
[in]drcv_dtThe partial derivative of coordinate defining potential density with potential temperature, in kg m-3 K-1.
[in]drcv_dsThe partial derivative of coordinate defining potential density with salinity, in kg m-3 psu-1.

Definition at line 2055 of file MOM_bulk_mixed_layer.F90.

Referenced by bulkmixedlayer().

2055  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
2056  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid
2057  !! structure.
2058  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: h !< Layer thickness, in m or kg m-2.
2059  !! (Intent in/out) The units of h
2060  !! are referred to as H below.
2061  !! Layer 0 is the new mixed layer.
2062  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: t !< Layer temperatures, in deg C.
2063  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: s !< Layer salinities, in psu.
2064  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: r0 !< Potential density referenced to
2065  !! surface pressure, in kg m-3.
2066  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: rcv !< The coordinate defining
2067  !! potential density, in kg m-3.
2068  real, dimension(SZK_(GV)), intent(in) :: rcvtgt !< The target value of Rcv for each
2069  !! layer, in kg m-3.
2070  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: eps !< The (small) thickness that must
2071  !! remain in each layer, in H.
2072  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_ea !< The upward increase across a
2073  !! layer in the entrainment from
2074  !! above, in m or kg m-2 (H).
2075  !! Positive d_ea goes with layer
2076  !! thickness increases.
2077  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_eb !< The downward increase across a
2078  !! layer in the entrainment from
2079  !! below, in H. Positive values go
2080  !! with mass gain by a layer.
2081  integer, dimension(SZI_(G),SZK_(GV)), intent(in) :: ksort !< The density-sorted k-indicies.
2082  type(bulkmixedlayer_cs), pointer :: cs !< The control structure for this
2083  !! module.
2084  real, dimension(SZI_(G)), intent(in) :: dr0_dt !< The partial derivative of
2085  !! potential density referenced
2086  !! to the surface with potential
2087  !! temperature, in kg m-3 K-1.
2088  real, dimension(SZI_(G)), intent(in) :: dr0_ds !< The partial derivative of
2089  !! cpotential density referenced
2090  !! to the surface with salinity,
2091  !! in kg m-3 psu-1.
2092  real, dimension(SZI_(G)), intent(in) :: drcv_dt !< The partial derivative of
2093  !! coordinate defining potential
2094  !! density with potential
2095  !! temperature, in kg m-3 K-1.
2096  real, dimension(SZI_(G)), intent(in) :: drcv_ds !< The partial derivative of
2097  !! coordinate defining potential
2098  !! density with salinity,
2099  !! in kg m-3 psu-1.
2100 
2101 ! This subroutine actually moves properties between layers to achieve a
2102 ! resorted state, with all of the resorted water either moved into the correct
2103 ! interior layers or in the top nkmb layers.
2104 !
2105 ! Arguments: h - Layer thickness, in m or kg m-2. (Intent in/out) The units of
2106 ! h are referred to as H below. Layer 0 is the new mixed layer.
2107 ! (in/out) T - Layer temperatures, in deg C.
2108 ! (in/out) S - Layer salinities, in psu.
2109 ! (in/out) R0 - Potential density referenced to surface pressure, in kg m-3.
2110 ! (in/out) Rcv - The coordinate defining potential density, in kg m-3.
2111 ! (in) RcvTgt - The target value of Rcv for each layer, in kg m-3.
2112 ! (in) eps - The (small) thickness that must remain in each layer, in H.
2113 ! (in/out) d_ea - The upward increase across a layer in the entrainment from
2114 ! above, in m or kg m-2 (H). Positive d_ea goes with layer
2115 ! thickness increases.
2116 ! (in/out) d_eb - The downward increase across a layer in the entrainment from
2117 ! below, in H. Positive values go with mass gain by a layer.
2118 ! (in) ksort - The density-sorted k-indicies.
2119 ! (in) G - The ocean's grid structure.
2120 ! (in) GV - The ocean's vertical grid structure.
2121 ! (in) CS - The control structure for this module.
2122 ! (in/out) dR0_dT - The partial derivative of potential density referenced
2123 ! to the surface with potential temperature, in kg m-3 K-1.
2124 ! (in/out) dR0_dS - The partial derivative of cpotential density referenced
2125 ! to the surface with salinity, in kg m-3 psu-1.
2126 ! (in/out) dRcv_dT - The partial derivative of coordinate defining potential
2127 ! density with potential temperature, in kg m-3 K-1.
2128 ! (in/out) dRcv_dS - The partial derivative of coordinate defining potential
2129 ! density with salinity, in kg m-3 psu-1.
2130 
2131 ! If there are no massive light layers above the deepest of the mixed- and
2132 ! buffer layers, do nothing (except perhaps to reshuffle these layers).
2133 ! If there are nkbl or fewer layers above the deepest mixed- or buffer-
2134 ! layers, move them (in sorted order) into the buffer layers, even if they
2135 ! were previously interior layers.
2136 ! If there are interior layers that are intermediate in density (both in-situ
2137 ! and the coordinate density (sigma-2)) between the newly forming mixed layer
2138 ! and a residual buffer- or mixed layer, and the number of massive layers above
2139 ! the deepest massive buffer or mixed layer is greater than nkbl, then split
2140 ! those buffer layers into peices that match the target density of the two
2141 ! nearest interior layers.
2142 ! Otherwise, if there are more than nkbl+1 remaining massive layers
2143  real :: h_move, h_tgt_old, i_hnew
2144  real :: dt_ds_wt2, dt_dr, ds_dr, i_denom
2145  real :: rcv_int
2146  real :: target_match_tol
2147  real :: t_up, s_up, r0_up, i_hup, h_to_up
2148  real :: t_dn, s_dn, r0_dn, i_hdn, h_to_dn
2149  real :: wt_dn
2150  real :: dr1, dr2
2151  real :: dpe, hmin, min_dpe, min_hmin
2152  real, dimension(SZK_(GV)) :: &
2153  h_tmp, r0_tmp, t_tmp, s_tmp, rcv_tmp
2154  integer :: ks_min
2155  logical :: sorted, leave_in_layer
2156  integer :: ks_deep(szi_(g)), k_count(szi_(g)), ks2_reverse(szi_(g), szk_(gv))
2157  integer :: ks2(szk_(gv))
2158  integer :: i, k, ks, is, ie, nz, k1, k2, k_tgt, k_src, k_int_top
2159  integer :: nks, nkmb, num_interior, top_interior_ks
2160 
2161  is = g%isc ; ie = g%iec ; nz = gv%ke
2162  nkmb = cs%nkml+cs%nkbl
2163  target_match_tol = 0.1 ! ### MAKE THIS A PARAMETER.
2164 
2165  dt_ds_wt2 = cs%dT_dS_wt**2
2166 
2167 ! Find out how many massive layers are above the deepest buffer or mixed layer.
2168  do i=is,ie ; ks_deep(i) = -1 ; k_count(i) = 0 ; enddo
2169  do ks=nz,1,-1 ; do i=is,ie ; if (ksort(i,ks) > 0) then
2170  k = ksort(i,ks)
2171 
2172  if (h(i,k) > eps(i,k)) then
2173  if (ks_deep(i) == -1) then
2174  if (k <= nkmb) then
2175  ks_deep(i) = ks ; k_count(i) = k_count(i) + 1
2176  ks2_reverse(i,k_count(i)) = k
2177  endif
2178  else
2179  k_count(i) = k_count(i) + 1
2180  ks2_reverse(i,k_count(i)) = k
2181  endif
2182  endif
2183  endif ; enddo ; enddo
2184 
2185  do i=is,ie ; if (k_count(i) > 1) then
2186  ! This column might need to be reshuffled.
2187  nks = k_count(i)
2188 
2189  ! Put ks2 in the right order and check whether reshuffling is needed.
2190  sorted = .true.
2191  ks2(nks) = ks2_reverse(i,1)
2192  do ks=nks-1,1,-1
2193  ks2(ks) = ks2_reverse(i,1+nks-ks)
2194  if (ks2(ks) > ks2(ks+1)) sorted = .false.
2195  enddo
2196 
2197  ! Go to the next column of no reshuffling is needed.
2198  if (sorted) cycle
2199 
2200  ! Find out how many interior layers are being reshuffled. If none,
2201  ! then this is a simple swapping procedure.
2202  num_interior = 0 ; top_interior_ks = 0
2203  do ks=nks,1,-1 ; if (ks2(ks) > nkmb) then
2204  num_interior = num_interior+1 ; top_interior_ks = ks
2205  endif ; enddo
2206 
2207  if (num_interior >= 1) then
2208  ! Find the lightest interior layer with a target coordinate density
2209  ! greater than the newly forming mixed layer.
2210  do k=nkmb+1,nz ; if (rcv(i,0) < rcvtgt(k)) exit ; enddo
2211  k_int_top = k ; rcv_int = rcvtgt(k)
2212 
2213  i_denom = 1.0 / (drcv_ds(i)**2 + dt_ds_wt2*drcv_dt(i)**2)
2214  dt_dr = dt_ds_wt2*drcv_dt(i) * i_denom
2215  ds_dr = drcv_ds(i) * i_denom
2216 
2217 
2218  ! Examine whether layers can be split out of existence. Stop when there
2219  ! is a layer that cannot be handled this way, or when the topmost formerly
2220  ! interior layer has been dealt with.
2221  do ks = nks,top_interior_ks,-1
2222  k = ks2(ks)
2223  leave_in_layer = .false.
2224  if ((k > nkmb) .and. (rcv(i,k) <= rcvtgt(k))) then
2225  if (rcvtgt(k)-rcv(i,k) < target_match_tol*(rcvtgt(k) - rcvtgt(k-1))) &
2226  leave_in_layer = .true.
2227  elseif (k > nkmb) then
2228  if (rcv(i,k)-rcvtgt(k) < target_match_tol*(rcvtgt(k+1) - rcvtgt(k))) &
2229  leave_in_layer = .true.
2230  endif
2231 
2232  if (leave_in_layer) then
2233  ! Just drop this layer from the sorted list.
2234  nks = nks-1
2235  elseif (rcv(i,k) < rcv_int) then
2236  ! There is no interior layer with a target density that is intermediate
2237  ! between this layer and the mixed layer.
2238  exit
2239  else
2240  ! Try splitting the layer between two interior isopycnal layers.
2241  ! Find the target densities that bracket this layer.
2242  do k2=k_int_top+1,nz ; if (rcv(i,k) < rcvtgt(k2)) exit ; enddo
2243  if (k2>nz) exit
2244 
2245  ! This layer is bracketed in density between layers k2-1 and k2.
2246 
2247  dr1 = (rcvtgt(k2-1) - rcv(i,k)) ; dr2 = (rcvtgt(k2) - rcv(i,k))
2248  t_up = t(i,k) + dt_dr * dr1
2249  s_up = s(i,k) + ds_dr * dr1
2250  t_dn = t(i,k) + dt_dr * dr2
2251  s_dn = s(i,k) + ds_dr * dr2
2252 
2253  r0_up = r0(i,k) + (dt_dr*dr0_dt(i) + ds_dr*dr0_ds(i)) * dr1
2254  r0_dn = r0(i,k) + (dt_dr*dr0_dt(i) + ds_dr*dr0_ds(i)) * dr2
2255 
2256  ! Make sure the new properties are acceptable.
2257  if ((r0_up > r0(i,0)) .or. (r0_dn > r0(i,0))) &
2258  ! Avoid creating obviously unstable profiles.
2259  exit
2260 
2261  wt_dn = (rcv(i,k) - rcvtgt(k2-1)) / (rcvtgt(k2) - rcvtgt(k2-1))
2262  h_to_up = (h(i,k)-eps(i,k)) * (1.0 - wt_dn)
2263  h_to_dn = (h(i,k)-eps(i,k)) * wt_dn
2264 
2265  i_hup = 1.0 / (h(i,k2-1) + h_to_up)
2266  i_hdn = 1.0 / (h(i,k2) + h_to_dn)
2267  r0(i,k2-1) = (r0(i,k2)*h(i,k2-1) + r0_up*h_to_up) * i_hup
2268  r0(i,k2) = (r0(i,k2)*h(i,k2) + r0_dn*h_to_dn) * i_hdn
2269 
2270  t(i,k2-1) = (t(i,k2)*h(i,k2-1) + t_up*h_to_up) * i_hup
2271  t(i,k2) = (t(i,k2)*h(i,k2) + t_dn*h_to_dn) * i_hdn
2272  s(i,k2-1) = (s(i,k2)*h(i,k2-1) + s_up*h_to_up) * i_hup
2273  s(i,k2) = (s(i,k2)*h(i,k2) + s_dn*h_to_dn) * i_hdn
2274  rcv(i,k2-1) = (rcv(i,k2)*h(i,k2-1) + rcvtgt(k2-1)*h_to_up) * i_hup
2275  rcv(i,k2) = (rcv(i,k2)*h(i,k2) + rcvtgt(k2)*h_to_dn) * i_hdn
2276 
2277  h(i,k) = eps(i,k)
2278  h(i,k2) = h(i,k2) + h_to_dn
2279  h(i,k2-1) = h(i,k2-1) + h_to_up
2280 
2281  if (k > k2-1) then
2282  d_eb(i,k) = d_eb(i,k) - h_to_up
2283  d_eb(i,k2-1) = d_eb(i,k2-1) + h_to_up
2284  elseif (k < k2-1) then
2285  d_ea(i,k) = d_ea(i,k) - h_to_up
2286  d_ea(i,k2-1) = d_ea(i,k2-1) + h_to_up
2287  endif
2288  if (k > k2) then
2289  d_eb(i,k) = d_eb(i,k) - h_to_dn
2290  d_eb(i,k2) = d_eb(i,k2) + h_to_dn
2291  elseif (k < k2) then
2292  d_ea(i,k) = d_ea(i,k) - h_to_dn
2293  d_ea(i,k2) = d_ea(i,k2) + h_to_dn
2294  endif
2295  nks = nks-1
2296  endif
2297  enddo
2298  endif
2299 
2300  do while (nks > nkmb)
2301  ! Having already tried to move surface layers into the interior, there
2302  ! are still too many layers, and layers must be merged until nks=nkmb.
2303  ! Examine every merger of a layer with its neighbors, and merge the ones
2304  ! that increase the potential energy the least. If there are layers
2305  ! with (apparently?) negative potential energy change, choose the one
2306  ! with the smallest total thickness. Repeat until nkmb layers remain.
2307  ! Choose the smaller value for the remaining index for convenience.
2308 
2309  ks_min = -1 ; min_dpe = 1.0 ; min_hmin = 0.0
2310  do ks=1,nks-1
2311  k1 = ks2(ks) ; k2 = ks2(ks+1)
2312  dpe = max(0.0, (r0(i,k2)-r0(i,k1)) * h(i,k1) * h(i,k2))
2313  hmin = min(h(i,k1)-eps(i,k1), h(i,k2)-eps(i,k2))
2314  if ((ks_min < 0) .or. (dpe < min_dpe) .or. &
2315  ((dpe <= 0.0) .and. (hmin < min_hmin))) then
2316  ks_min = ks ; min_dpe = dpe ; min_hmin = hmin
2317  endif
2318  enddo
2319 
2320  ! Now merge the two layers that do the least damage.
2321  k1 = ks2(ks_min) ; k2 = ks2(ks_min+1)
2322  if (k1 < k2) then ; k_tgt = k1 ; k_src = k2
2323  else ; k_tgt = k2 ; k_src = k1 ; ks2(ks_min) = k2 ; endif
2324 
2325  h_tgt_old = h(i,k_tgt)
2326  h_move = h(i,k_src)-eps(i,k_src)
2327  h(i,k_src) = eps(i,k_src)
2328  h(i,k_tgt) = h(i,k_tgt) + h_move
2329  i_hnew = 1.0 / (h(i,k_tgt))
2330  r0(i,k_tgt) = (r0(i,k_tgt)*h_tgt_old + r0(i,k_src)*h_move) * i_hnew
2331 
2332  t(i,k_tgt) = (t(i,k_tgt)*h_tgt_old + t(i,k_src)*h_move) * i_hnew
2333  s(i,k_tgt) = (s(i,k_tgt)*h_tgt_old + s(i,k_src)*h_move) * i_hnew
2334  rcv(i,k_tgt) = (rcv(i,k_tgt)*h_tgt_old + rcv(i,k_src)*h_move) * i_hnew
2335 
2336  d_eb(i,k_src) = d_eb(i,k_src) - h_move
2337  d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_move
2338 
2339  ! Remove the newly missing layer from the sorted list.
2340  do ks=ks_min+1,nks ; ks2(ks) = ks2(ks+1) ; enddo
2341  nks = nks-1
2342  enddo
2343 
2344  ! Check again whether the layers are sorted, and go on to the next column
2345  ! if they are.
2346  sorted = .true.
2347  do ks=1,nks-1 ; if (ks2(ks) > ks2(ks+1)) sorted = .false. ; enddo
2348  if (sorted) cycle
2349 
2350  if (nks > 1) then
2351  ! At this point, actually swap the properties of the layers, and place
2352  ! the remaining layers in order starting with nkmb.
2353 
2354  ! Save all the properties of the nkmb layers that might be replaced.
2355  do k=1,nkmb
2356  h_tmp(k) = h(i,k) ; r0_tmp(k) = r0(i,k)
2357  t_tmp(k) = t(i,k) ; s_tmp(k) = s(i,k) ; rcv_tmp(k) = rcv(i,k)
2358 
2359  h(i,k) = 0.0
2360  enddo
2361 
2362  do ks=nks,1,-1
2363  k_tgt = nkmb - nks + ks ; k_src = ks2(ks)
2364  if (k_tgt == k_src) then
2365  h(i,k_tgt) = h_tmp(k_tgt) ! This layer doesn't move, so put the water back.
2366  cycle
2367  endif
2368 
2369  ! Note below that eps=0 for k<=nkmb.
2370  if (k_src > nkmb) then
2371  h_move = h(i,k_src)-eps(i,k_src)
2372  h(i,k_src) = eps(i,k_src)
2373  h(i,k_tgt) = h_move
2374  r0(i,k_tgt) = r0(i,k_src)
2375 
2376  t(i,k_tgt) = t(i,k_src) ; s(i,k_tgt) = s(i,k_src)
2377  rcv(i,k_tgt) = rcv(i,k_src)
2378 
2379  d_eb(i,k_src) = d_eb(i,k_src) - h_move
2380  d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_move
2381  else
2382  h(i,k_tgt) = h_tmp(k_src)
2383  r0(i,k_tgt) = r0_tmp(k_src)
2384 
2385  t(i,k_tgt) = t_tmp(k_src) ; s(i,k_tgt) = s_tmp(k_src)
2386  rcv(i,k_tgt) = rcv_tmp(k_src)
2387 
2388  if (k_src > k_tgt) then
2389  d_eb(i,k_src) = d_eb(i,k_src) - h_tmp(k_src)
2390  d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_tmp(k_src)
2391  else
2392  d_ea(i,k_src) = d_ea(i,k_src) - h_tmp(k_src)
2393  d_ea(i,k_tgt) = d_ea(i,k_tgt) + h_tmp(k_src)
2394  endif
2395  endif
2396  enddo
2397  endif
2398 
2399  endif ; enddo
2400 
Here is the caller graph for this function:

◆ sort_ml()

subroutine mom_bulk_mixed_layer::sort_ml ( real, dimension(szi_(g),szk_(gv)), intent(in)  h,
real, dimension(szi_(g),szk_(gv)), intent(in)  R0,
real, dimension(szi_(g),szk_(gv)), intent(in)  eps,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(bulkmixedlayer_cs), pointer  CS,
integer, dimension(szi_(g),szk_(gv)), intent(out)  ksort 
)
private

This subroutine generates an array of indices that are sorted by layer density.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]hLayer thickness, in m or kg m-2. (Intent in/out) The units of h are referred to as H below.
[in]r0The potential density used to sort the layers, in kg m-3.
[in]epsThe (small) thickness that must remain in each layer, in H.
csThe control structure returned by a previous call to mixedlayer_init.
[out]ksortThe k-index to use in the sort.

Definition at line 1987 of file MOM_bulk_mixed_layer.F90.

Referenced by bulkmixedlayer().

1987  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
1988  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
1989  real, dimension(SZI_(G),SZK_(GV)), intent(in) :: h !< Layer thickness, in m or kg m-2.
1990  !! (Intent in/out) The units of h are
1991  !! referred to as H below.
1992  real, dimension(SZI_(G),SZK_(GV)), intent(in) :: r0 !< The potential density used to sort
1993  !! the layers, in kg m-3.
1994  real, dimension(SZI_(G),SZK_(GV)), intent(in) :: eps !< The (small) thickness that must
1995  !! remain in each layer, in H.
1996  type(bulkmixedlayer_cs), pointer :: cs !< The control structure returned by a
1997  !! previous call to mixedlayer_init.
1998  integer, dimension(SZI_(G),SZK_(GV)), intent(out) :: ksort !< The k-index to use in the sort.
1999 
2000 ! This subroutine generates an array of indices that are sorted by layer
2001 ! density.
2002 
2003 ! Arguments: h - Layer thickness, in m or kg m-2. (Intent in/out) The units
2004 ! of h are referred to as H below.
2005 ! (in) R0 - The potential density used to sort the layers, in kg m-3.
2006 ! (in) eps - The (small) thickness that must remain in each layer, in H.
2007 ! (in) tv - A structure containing pointers to any available
2008 ! thermodynamic fields. Absent fields have NULL ptrs.
2009 ! (in) j - The meridional row to work on.
2010 ! (in) G - The ocean's grid structure.
2011 ! (in) GV - The ocean's vertical grid structure.
2012 ! (in) CS - The control structure returned by a previous call to
2013 ! mixedlayer_init.
2014 ! (out) ksort - The k-index to use in the sort.
2015  real :: r0sort(szi_(g),szk_(gv))
2016  integer :: nsort(szi_(g))
2017  logical :: done_sorting(szi_(g))
2018  integer :: i, k, ks, is, ie, nz, nkmb
2019 
2020  is = g%isc ; ie = g%iec ; nz = gv%ke
2021  nkmb = cs%nkml+cs%nkbl
2022 
2023  ! Come up with a sorted index list of layers with increasing R0.
2024  ! Assume that the layers below nkmb are already stably stratified.
2025  ! Only layers that are thicker than eps are in the list. Extra elements
2026  ! have an index of -1.
2027 
2028  ! This is coded using straight insertion, on the assumption that the
2029  ! layers are usually in the right order (or massless) anyway.
2030 
2031  do k=1,nz ; do i=is,ie ; ksort(i,k) = -1 ; enddo ; enddo
2032 
2033  do i=is,ie ; nsort(i) = 0 ; done_sorting(i) = .false. ; enddo
2034  do k=1,nz ; do i=is,ie ; if (h(i,k) > eps(i,k)) then
2035  if (done_sorting(i)) then ; ks = nsort(i) ; else
2036  do ks=nsort(i),1,-1
2037  if (r0(i,k) >= r0sort(i,ks)) exit
2038  r0sort(i,ks+1) = r0sort(i,ks) ; ksort(i,ks+1) = ksort(i,ks)
2039  enddo
2040  if ((k > nkmb) .and. (ks == nsort(i))) done_sorting(i) = .true.
2041  endif
2042 
2043  ksort(i,ks+1) = k
2044  r0sort(i,ks+1) = r0(i,k)
2045  nsort(i) = nsort(i) + 1
2046  endif ; enddo ; enddo
2047 
Here is the caller graph for this function:

Variable Documentation

◆ id_clock_adjustment

integer mom_bulk_mixed_layer::id_clock_adjustment =0
private

Definition at line 188 of file MOM_bulk_mixed_layer.F90.

Referenced by bulkmixedlayer(), and bulkmixedlayer_init().

◆ id_clock_conv

integer mom_bulk_mixed_layer::id_clock_conv =0
private

Definition at line 188 of file MOM_bulk_mixed_layer.F90.

Referenced by bulkmixedlayer(), and bulkmixedlayer_init().

◆ id_clock_detrain

integer mom_bulk_mixed_layer::id_clock_detrain =0

Definition at line 188 of file MOM_bulk_mixed_layer.F90.

Referenced by bulkmixedlayer(), and bulkmixedlayer_init().

188 integer :: id_clock_detrain=0, id_clock_mech=0, id_clock_conv=0, id_clock_adjustment=0

◆ id_clock_eos

integer mom_bulk_mixed_layer::id_clock_eos =0
private

Definition at line 189 of file MOM_bulk_mixed_layer.F90.

Referenced by bulkmixedlayer(), and bulkmixedlayer_init().

189 integer :: id_clock_eos=0, id_clock_resort=0, id_clock_pass=0

◆ id_clock_mech

integer mom_bulk_mixed_layer::id_clock_mech =0
private

Definition at line 188 of file MOM_bulk_mixed_layer.F90.

Referenced by bulkmixedlayer(), and bulkmixedlayer_init().

◆ id_clock_pass

integer mom_bulk_mixed_layer::id_clock_pass =0
private

Definition at line 189 of file MOM_bulk_mixed_layer.F90.

Referenced by bulkmixedlayer(), and bulkmixedlayer_init().

◆ id_clock_resort

integer mom_bulk_mixed_layer::id_clock_resort =0
private

Definition at line 189 of file MOM_bulk_mixed_layer.F90.

Referenced by bulkmixedlayer(), and bulkmixedlayer_init().

◆ max_msg

integer mom_bulk_mixed_layer::max_msg = 2
private

Definition at line 191 of file MOM_bulk_mixed_layer.F90.

◆ num_msg

integer mom_bulk_mixed_layer::num_msg = 0
private

Definition at line 191 of file MOM_bulk_mixed_layer.F90.

191 integer :: num_msg = 0, max_msg = 2