MOM6
Data Types
mom_bulk_mixed_layer Module Reference

Detailed Description

Build mixed layer parameterization.

By Robert Hallberg, 1997 - 2005.

This file contains the subroutine (bulkmixedlayer) that implements a Kraus-Turner-like bulk mixed layer, based on the work of various people, as described in the review paper by Niiler and Kraus (1979), with particular attention to the form proposed by Oberhuber (JPO, 1993, 808-829), with an extension to a refied bulk mixed layer as described in Hallberg (Aha Huliko'a, 2003). The physical processes portrayed in this subroutine include convective adjustment and mixed layer entrainment and detrainment. Penetrating shortwave radiation and an exponential decay of TKE fluxes are also supported by this subroutine. Several constants can alternately be set to give a traditional Kraus-Turner mixed layer scheme, although that is not the preferred option. The physical processes and arguments are described in detail below.

Data Types

type  bulkmixedlayer_cs
 The control structure with parameters for the MOM_bulk_mixed_layer module. More...
 
integer id_clock_detrain =0
 CPU clock IDs. More...
 
integer id_clock_mech =0
 CPU clock IDs. More...
 
integer id_clock_conv =0
 CPU clock IDs. More...
 
integer id_clock_adjustment =0
 CPU clock IDs. More...
 
integer id_clock_eos =0
 CPU clock IDs. More...
 
integer id_clock_resort =0
 CPU clock IDs. More...
 
integer id_clock_pass =0
 CPU clock IDs. More...
 
subroutine, public bulkmixedlayer (h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, 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, US, 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, US, 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, US, 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, US, 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, US, 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, US, 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, US, param_file, diag, CS)
 This subroutine initializes the MOM bulk mixed layer module. More...
 
real function ef4 (Ht, En, I_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/((Ht+E)*Ht). More...
 

Function/Subroutine Documentation

◆ bulkmixedlayer()

subroutine, public mom_bulk_mixed_layer::bulkmixedlayer ( real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(inout)  h_3d,
real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)  u_3d,
real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)  v_3d,
type(thermo_var_ptrs), intent(inout)  tv,
type(forcing), intent(inout)  fluxes,
real, intent(in)  dt,
real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(inout)  ea,
real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(inout)  eb,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
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]usA dimensional unit scaling type
[in,out]h_3dLayer thickness [H ~> 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 [T ~> 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 [m-1].
hmlActive mixed layer depth [m].
[in]aggregate_fw_forcingIf true, the net incoming and outgoing surface freshwater fluxes are combined before being applied, instead of being applied separately.
[in]dt_diagThe diagnostic time step, which may be less than dt if there are two callse to mixedlayer [T ~> 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 189 of file MOM_bulk_mixed_layer.F90.

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

References mom_opacity::absorbremainingsw(), convective_adjustment(), mom_diag_mediator::diag_update_remap_grids(), mom_domains::do_group_pass(), mom_opacity::extract_optics_slice(), 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::layered_diabatic().

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(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(bulkmixedlayer_cs), pointer  CS 
)

This subroutine initializes the MOM bulk mixed layer module.

Parameters
[in]timeThe model's clock with the current time.
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]param_fileA structure to parse for run-time parameters.
[in,out]diagA structure that is used to regulate diagnostic output.
csA pointer that is set to point to the control structure for this module.

Definition at line 3382 of file MOM_bulk_mixed_layer.F90.

3382  type(time_type), target, intent(in) :: Time !< The model's clock with the current time.
3383  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
3384  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
3385  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3386  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
3387  !! parameters.
3388  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
3389  !! output.
3390  type(bulkmixedlayer_CS), pointer :: CS !< A pointer that is set to point to the control
3391  !! structure for this module.
3392 ! Arguments: Time - The current model time.
3393 ! (in) G - The ocean's grid structure.
3394 ! (in) GV - The ocean's vertical grid structure.
3395 ! (in) param_file - A structure indicating the open file to parse for
3396 ! model parameter values.
3397 ! (in) diag - A structure that is used to regulate diagnostic output.
3398 ! (in/out) CS - A pointer that is set to point to the control structure
3399 ! for this module
3400 ! This include declares and sets the variable "version".
3401 #include "version_variable.h"
3402  character(len=40) :: mdl = "MOM_mixed_layer" ! This module's name.
3403  real :: BL_detrain_time_dflt ! The default value for BUFFER_LAY_DETRAIN_TIME [s]
3404  real :: omega_frac_dflt, ustar_min_dflt, Hmix_min_m
3405  integer :: isd, ied, jsd, jed
3406  logical :: use_temperature, use_omega
3407  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3408 
3409  if (associated(cs)) then
3410  call mom_error(warning, "mixedlayer_init called with an associated"// &
3411  "associated control structure.")
3412  return
3413  else ; allocate(cs) ; endif
3414 
3415  cs%diag => diag
3416  cs%Time => time
3417 
3418  if (gv%nkml < 1) return
3419 
3420 ! Set default, read and log parameters
3421  call log_version(param_file, mdl, version, "")
3422 
3423  cs%nkml = gv%nkml
3424  call log_param(param_file, mdl, "NKML", cs%nkml, &
3425  "The number of sublayers within the mixed layer if "//&
3426  "BULKMIXEDLAYER is true.", units="nondim", default=2)
3427  cs%nkbl = gv%nk_rho_varies - gv%nkml
3428  call log_param(param_file, mdl, "NKBL", cs%nkbl, &
3429  "The number of variable density buffer layers if "//&
3430  "BULKMIXEDLAYER is true.", units="nondim", default=2)
3431  call get_param(param_file, mdl, "MSTAR", cs%mstar, &
3432  "The ratio of the friction velocity cubed to the TKE "//&
3433  "input to the mixed layer.", "units=nondim", default=1.2)
3434  call get_param(param_file, mdl, "NSTAR", cs%nstar, &
3435  "The portion of the buoyant potential energy imparted by "//&
3436  "surface fluxes that is available to drive entrainment "//&
3437  "at the base of mixed layer when that energy is positive.", &
3438  units="nondim", default=0.15)
3439  call get_param(param_file, mdl, "BULK_RI_ML", cs%bulk_Ri_ML, &
3440  "The efficiency with which mean kinetic energy released "//&
3441  "by mechanically forced entrainment of the mixed layer "//&
3442  "is converted to turbulent kinetic energy.", units="nondim",&
3443  fail_if_missing=.true.)
3444  call get_param(param_file, mdl, "ABSORB_ALL_SW", cs%absorb_all_sw, &
3445  "If true, all shortwave radiation is absorbed by the "//&
3446  "ocean, instead of passing through to the bottom mud.", &
3447  default=.false.)
3448  call get_param(param_file, mdl, "TKE_DECAY", cs%TKE_decay, &
3449  "TKE_DECAY relates the vertical rate of decay of the "//&
3450  "TKE available for mechanical entrainment to the natural "//&
3451  "Ekman depth.", units="nondim", default=2.5)
3452  call get_param(param_file, mdl, "NSTAR2", cs%nstar2, &
3453  "The portion of any potential energy released by "//&
3454  "convective adjustment that is available to drive "//&
3455  "entrainment at the base of mixed layer. By default "//&
3456  "NSTAR2=NSTAR.", units="nondim", default=cs%nstar)
3457  call get_param(param_file, mdl, "BULK_RI_CONVECTIVE", cs%bulk_Ri_convective, &
3458  "The efficiency with which convectively released mean "//&
3459  "kinetic energy is converted to turbulent kinetic "//&
3460  "energy. By default BULK_RI_CONVECTIVE=BULK_RI_ML.", &
3461  units="nondim", default=cs%bulk_Ri_ML)
3462  call get_param(param_file, mdl, "HMIX_MIN", cs%Hmix_min, &
3463  "The minimum mixed layer depth if the mixed layer depth "//&
3464  "is determined dynamically.", units="m", default=0.0, scale=gv%m_to_H, &
3465  unscaled=hmix_min_m)
3466 
3467  call get_param(param_file, mdl, "LIMIT_BUFFER_DETRAIN", cs%limit_det, &
3468  "If true, limit the detrainment from the buffer layers "//&
3469  "to not be too different from the neighbors.", default=.false.)
3470  call get_param(param_file, mdl, "ALLOWED_DETRAIN_TEMP_CHG", cs%Allowed_T_chg, &
3471  "The amount by which temperature is allowed to exceed "//&
3472  "previous values during detrainment.", units="K", default=0.5)
3473  call get_param(param_file, mdl, "ALLOWED_DETRAIN_SALT_CHG", cs%Allowed_S_chg, &
3474  "The amount by which salinity is allowed to exceed "//&
3475  "previous values during detrainment.", units="PSU", default=0.1)
3476  call get_param(param_file, mdl, "ML_DT_DS_WEIGHT", cs%dT_dS_wt, &
3477  "When forced to extrapolate T & S to match the layer "//&
3478  "densities, this factor (in deg C / PSU) is combined "//&
3479  "with the derivatives of density with T & S to determine "//&
3480  "what direction is orthogonal to density contours. It "//&
3481  "should be a typical value of (dR/dS) / (dR/dT) in "//&
3482  "oceanic profiles.", units="degC PSU-1", default=6.0)
3483  call get_param(param_file, mdl, "BUFFER_LAYER_EXTRAP_LIMIT", cs%BL_extrap_lim, &
3484  "A limit on the density range over which extrapolation "//&
3485  "can occur when detraining from the buffer layers, "//&
3486  "relative to the density range within the mixed and "//&
3487  "buffer layers, when the detrainment is going into the "//&
3488  "lightest interior layer, nondimensional, or a negative "//&
3489  "value not to apply this limit.", units="nondim", default = -1.0)
3490  call get_param(param_file, mdl, "BUFFER_LAYER_HMIN_THICK", cs%Hbuffer_min, &
3491  "The minimum buffer layer thickness when the mixed layer is very thick.", &
3492  units="m", default=5.0, scale=gv%m_to_H)
3493  call get_param(param_file, mdl, "BUFFER_LAYER_HMIN_REL", cs%Hbuffer_rel_min, &
3494  "The minimum buffer layer thickness relative to the combined mixed "//&
3495  "land buffer ayer thicknesses when they are thin.", &
3496  units="nondim", default=0.1/cs%nkbl)
3497  bl_detrain_time_dflt = 4.0*3600.0 ; if (cs%nkbl==1) bl_detrain_time_dflt = 86400.0*30.0
3498  call get_param(param_file, mdl, "BUFFER_LAY_DETRAIN_TIME", cs%BL_detrain_time, &
3499  "A timescale that characterizes buffer layer detrainment events.", &
3500  units="s", default=bl_detrain_time_dflt, scale=us%s_to_T)
3501  call get_param(param_file, mdl, "BUFFER_SPLIT_RHO_TOL", cs%BL_split_rho_tol, &
3502  "The fractional tolerance for matching layer target densities when splitting "//&
3503  "layers to deal with massive interior layers that are lighter than one of the "//&
3504  "mixed or buffer layers.", units="nondim", default=0.1)
3505 
3506  call get_param(param_file, mdl, "DEPTH_LIMIT_FLUXES", cs%H_limit_fluxes, &
3507  "The surface fluxes are scaled away when the total ocean "//&
3508  "depth is less than DEPTH_LIMIT_FLUXES.", &
3509  units="m", default=0.1*hmix_min_m, scale=gv%m_to_H)
3510  call get_param(param_file, mdl, "OMEGA", cs%omega, &
3511  "The rotation rate of the earth.", &
3512  default=7.2921e-5, units="s-1", scale=us%T_to_s)
3513  call get_param(param_file, mdl, "ML_USE_OMEGA", use_omega, &
3514  "If true, use the absolute rotation rate instead of the "//&
3515  "vertical component of rotation when setting the decay "//&
3516  "scale for turbulence.", default=.false., do_not_log=.true.)
3517  omega_frac_dflt = 0.0
3518  if (use_omega) then
3519  call mom_error(warning, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
3520  omega_frac_dflt = 1.0
3521  endif
3522  call get_param(param_file, mdl, "ML_OMEGA_FRAC", cs%omega_frac, &
3523  "When setting the decay scale for turbulence, use this "//&
3524  "fraction of the absolute rotation rate blended with the "//&
3525  "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
3526  units="nondim", default=omega_frac_dflt)
3527  call get_param(param_file, mdl, "ML_RESORT", cs%ML_resort, &
3528  "If true, resort the topmost layers by potential density "//&
3529  "before the mixed layer calculations.", default=.false.)
3530  if (cs%ML_resort) &
3531  call get_param(param_file, mdl, "ML_PRESORT_NK_CONV_ADJ", cs%ML_presort_nz_conv_adj, &
3532  "Convectively mix the first ML_PRESORT_NK_CONV_ADJ "//&
3533  "layers before sorting when ML_RESORT is true.", &
3534  units="nondim", default=0, fail_if_missing=.true.) ! Fail added by AJA.
3535  ! This gives a minimum decay scale that is typically much less than Angstrom.
3536  ustar_min_dflt = 2e-4*us%s_to_T*cs%omega*(gv%Angstrom_m + gv%H_to_m*gv%H_subroundoff)
3537  call get_param(param_file, mdl, "BML_USTAR_MIN", cs%ustar_min, &
3538  "The minimum value of ustar that should be used by the "//&
3539  "bulk mixed layer model in setting vertical TKE decay "//&
3540  "scales. This must be greater than 0.", units="m s-1", &
3541  default=ustar_min_dflt, scale=us%m_to_Z*us%T_to_s)
3542  if (cs%ustar_min<=0.0) call mom_error(fatal, "BML_USTAR_MIN must be positive.")
3543 
3544  call get_param(param_file, mdl, "RESOLVE_EKMAN", cs%Resolve_Ekman, &
3545  "If true, the NKML>1 layers in the mixed layer are "//&
3546  "chosen to optimally represent the impact of the Ekman "//&
3547  "transport on the mixed layer TKE budget. Otherwise, "//&
3548  "the sublayers are distributed uniformly through the "//&
3549  "mixed layer.", default=.false.)
3550  call get_param(param_file, mdl, "CORRECT_ABSORPTION_DEPTH", cs%correct_absorption, &
3551  "If true, the average depth at which penetrating shortwave "//&
3552  "radiation is absorbed is adjusted to match the average "//&
3553  "heating depth of an exponential profile by moving some "//&
3554  "of the heating upward in the water column.", default=.false.)
3555  call get_param(param_file, mdl, "DO_RIVERMIX", cs%do_rivermix, &
3556  "If true, apply additional mixing wherever there is "//&
3557  "runoff, so that it is mixed down to RIVERMIX_DEPTH, "//&
3558  "if the ocean is that deep.", default=.false.)
3559  if (cs%do_rivermix) &
3560  call get_param(param_file, mdl, "RIVERMIX_DEPTH", cs%rivermix_depth, &
3561  "The depth to which rivers are mixed if DO_RIVERMIX is "//&
3562  "defined.", units="m", default=0.0, scale=us%m_to_Z)
3563  call get_param(param_file, mdl, "USE_RIVER_HEAT_CONTENT", cs%use_river_heat_content, &
3564  "If true, use the fluxes%runoff_Hflx field to set the "//&
3565  "heat carried by runoff, instead of using SST*CP*liq_runoff.", &
3566  default=.false.)
3567  call get_param(param_file, mdl, "USE_CALVING_HEAT_CONTENT", cs%use_calving_heat_content, &
3568  "If true, use the fluxes%calving_Hflx field to set the "//&
3569  "heat carried by runoff, instead of using SST*CP*froz_runoff.", &
3570  default=.false.)
3571 
3572  call get_param(param_file, mdl, "ALLOW_CLOCKS_IN_OMP_LOOPS", &
3573  cs%allow_clocks_in_omp_loops, &
3574  "If true, clocks can be called from inside loops that can "//&
3575  "be threaded. To run with multiple threads, set to False.", &
3576  default=.true.)
3577 
3578  cs%id_ML_depth = register_diag_field('ocean_model', 'h_ML', diag%axesT1, &
3579  time, 'Surface mixed layer depth', 'm')
3580  cs%id_TKE_wind = register_diag_field('ocean_model', 'TKE_wind', diag%axesT1, &
3581  time, 'Wind-stirring source of mixed layer TKE', &
3582  'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3583  cs%id_TKE_RiBulk = register_diag_field('ocean_model', 'TKE_RiBulk', diag%axesT1, &
3584  time, 'Mean kinetic energy source of mixed layer TKE', &
3585  'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3586  cs%id_TKE_conv = register_diag_field('ocean_model', 'TKE_conv', diag%axesT1, &
3587  time, 'Convective source of mixed layer TKE', 'm3 s-3', &
3588  conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3589  cs%id_TKE_pen_SW = register_diag_field('ocean_model', 'TKE_pen_SW', diag%axesT1, &
3590  time, 'TKE consumed by mixing penetrative shortwave radation through the mixed layer', &
3591  'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3592  cs%id_TKE_mixing = register_diag_field('ocean_model', 'TKE_mixing', diag%axesT1, &
3593  time, 'TKE consumed by mixing that deepens the mixed layer', &
3594  'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3595  cs%id_TKE_mech_decay = register_diag_field('ocean_model', 'TKE_mech_decay', diag%axesT1, &
3596  time, 'Mechanical energy decay sink of mixed layer TKE', &
3597  'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3598  cs%id_TKE_conv_decay = register_diag_field('ocean_model', 'TKE_conv_decay', diag%axesT1, &
3599  time, 'Convective energy decay sink of mixed layer TKE', &
3600  'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3601  cs%id_TKE_conv_s2 = register_diag_field('ocean_model', 'TKE_conv_s2', diag%axesT1, &
3602  time, 'Spurious source of mixed layer TKE from sigma2', &
3603  'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3604  cs%id_PE_detrain = register_diag_field('ocean_model', 'PE_detrain', diag%axesT1, &
3605  time, 'Spurious source of potential energy from mixed layer detrainment', &
3606  'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3607  cs%id_PE_detrain2 = register_diag_field('ocean_model', 'PE_detrain2', diag%axesT1, &
3608  time, 'Spurious source of potential energy from mixed layer only detrainment', &
3609  'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
3610  cs%id_h_mismatch = register_diag_field('ocean_model', 'h_miss_ML', diag%axesT1, &
3611  time, 'Summed absolute mismatch in entrainment terms', 'm', conversion=us%Z_to_m)
3612  cs%id_Hsfc_used = register_diag_field('ocean_model', 'Hs_used', diag%axesT1, &
3613  time, 'Surface region thickness that is used', 'm', conversion=us%Z_to_m)
3614  cs%id_Hsfc_max = register_diag_field('ocean_model', 'Hs_max', diag%axesT1, &
3615  time, 'Maximum surface region thickness', 'm', conversion=us%Z_to_m)
3616  cs%id_Hsfc_min = register_diag_field('ocean_model', 'Hs_min', diag%axesT1, &
3617  time, 'Minimum surface region thickness', 'm', conversion=us%Z_to_m)
3618  !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?
3619  if (cs%limit_det .or. (cs%id_Hsfc_min > 0)) then
3620  call get_param(param_file, mdl, "LIMIT_BUFFER_DET_DH_SFC", cs%lim_det_dH_sfc, &
3621  "The fractional limit in the change between grid points "//&
3622  "of the surface region (mixed & buffer layer) thickness.", &
3623  units="nondim", default=0.5)
3624  call get_param(param_file, mdl, "LIMIT_BUFFER_DET_DH_BATHY", cs%lim_det_dH_bathy, &
3625  "The fraction of the total depth by which the thickness "//&
3626  "of the surface region (mixed & buffer layer) is allowed "//&
3627  "to change between grid points.", units="nondim", default=0.2)
3628  endif
3629 
3630  call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
3631  "If true, temperature and salinity are used as state "//&
3632  "variables.", default=.true.)
3633  cs%nsw = 0
3634  if (use_temperature) then
3635  call get_param(param_file, mdl, "PEN_SW_NBANDS", cs%nsw, default=1)
3636  endif
3637 
3638 
3639  if (max(cs%id_TKE_wind, cs%id_TKE_RiBulk, cs%id_TKE_conv, cs%id_TKE_mixing, &
3640  cs%id_TKE_pen_SW, cs%id_TKE_mech_decay, cs%id_TKE_conv_decay) > 0) then
3641  call safe_alloc_alloc(cs%diag_TKE_wind, isd, ied, jsd, jed)
3642  call safe_alloc_alloc(cs%diag_TKE_RiBulk, isd, ied, jsd, jed)
3643  call safe_alloc_alloc(cs%diag_TKE_conv, isd, ied, jsd, jed)
3644  call safe_alloc_alloc(cs%diag_TKE_pen_SW, isd, ied, jsd, jed)
3645  call safe_alloc_alloc(cs%diag_TKE_mixing, isd, ied, jsd, jed)
3646  call safe_alloc_alloc(cs%diag_TKE_mech_decay, isd, ied, jsd, jed)
3647  call safe_alloc_alloc(cs%diag_TKE_conv_decay, isd, ied, jsd, jed)
3648  call safe_alloc_alloc(cs%diag_TKE_conv_s2, isd, ied, jsd, jed)
3649 
3650  cs%TKE_diagnostics = .true.
3651  endif
3652  if (cs%id_PE_detrain > 0) call safe_alloc_alloc(cs%diag_PE_detrain, isd, ied, jsd, jed)
3653  if (cs%id_PE_detrain2 > 0) call safe_alloc_alloc(cs%diag_PE_detrain2, isd, ied, jsd, jed)
3654  if (cs%id_ML_depth > 0) call safe_alloc_alloc(cs%ML_depth, isd, ied, jsd, jed)
3655 
3656  if (cs%allow_clocks_in_omp_loops) then
3657  id_clock_detrain = cpu_clock_id('(Ocean mixed layer detrain)', grain=clock_routine)
3658  id_clock_mech = cpu_clock_id('(Ocean mixed layer mechanical entrainment)', grain=clock_routine)
3659  id_clock_conv = cpu_clock_id('(Ocean mixed layer convection)', grain=clock_routine)
3660  if (cs%ML_resort) then
3661  id_clock_resort = cpu_clock_id('(Ocean mixed layer resorting)', grain=clock_routine)
3662  else
3663  id_clock_adjustment = cpu_clock_id('(Ocean mixed layer convective adjustment)', grain=clock_routine)
3664  endif
3665  id_clock_eos = cpu_clock_id('(Ocean mixed layer EOS)', grain=clock_routine)
3666  endif
3667 
3668  if (cs%limit_det .or. (cs%id_Hsfc_min > 0)) &
3669  id_clock_pass = cpu_clock_id('(Ocean mixed layer halo updates)', grain=clock_routine)
3670 
3671 
3672 ! if (CS%limit_det) then
3673 ! endif
3674 

References id_clock_adjustment, id_clock_conv, id_clock_detrain, id_clock_eos, id_clock_mech, id_clock_pass, id_clock_resort, and mom_error_handler::mom_error().

Referenced by mom_diabatic_driver::diabatic_driver_init().

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

◆ convective_adjustment()

subroutine mom_bulk_mixed_layer::convective_adjustment ( real, dimension( g %isd: g %ied, gv %ke), intent(inout)  h,
real, dimension( g %isd: g %ied, gv %ke), intent(inout)  u,
real, dimension( g %isd: g %ied, gv %ke), intent(inout)  v,
real, dimension( g %isd: g %ied, gv %ke), intent(inout)  R0,
real, dimension( g %isd: g %ied, gv %ke), intent(inout)  Rcv,
real, dimension( g %isd: g %ied, gv %ke), intent(inout)  T,
real, dimension( g %isd: g %ied, gv %ke), intent(inout)  S,
real, dimension( g %isd: g %ied, gv %ke), intent(in)  eps,
real, dimension( g %isd: g %ied, gv %ke), intent(inout)  d_eb,
real, dimension( g %isd: g %ied, gv %ke), intent(out)  dKE_CA,
real, dimension( g %isd: g %ied, gv %ke), intent(out)  cTKE,
integer, intent(in)  j,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(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 [H ~> m or kg m-2]. The units of h are referred to as H below.
[in,out]uZonal velocities interpolated to h points [L T-1 ~> m s-1].
[in,out]vZonal velocities interpolated to h points [L T-1 ~> m s-1].
[in,out]tLayer temperatures [degC].
[in,out]sLayer salinities [ppt].
[in,out]r0Potential density referenced to surface pressure [R ~> kg m-3].
[in,out]rcvThe coordinate defining potential density [R ~> kg m-3].
[in,out]d_ebThe downward increase across a layer in the entrainment from below [H ~> m or kg m-2]. Positive values go with mass gain by a layer.
[in]epsThe negligibly small amount of water that will be left in each layer [H ~> m or kg m-2].
[out]dke_caThe vertically integrated change in kinetic energy due to convective adjustment [Z L2 T-2 ~> m3 s-2].
[out]ctkeThe buoyant turbulent kinetic energy source due to convective adjustment [Z L2 T-2 ~> m3 s-2].
[in]jThe j-index to work on.
[in]usA dimensional unit scaling type
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 803 of file MOM_bulk_mixed_layer.F90.

803  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
804  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
805  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2].
806  !! The units of h are referred to as H below.
807  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: u !< Zonal velocities interpolated to h
808  !! points [L T-1 ~> m s-1].
809  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: v !< Zonal velocities interpolated to h
810  !! points [L T-1 ~> m s-1].
811  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: T !< Layer temperatures [degC].
812  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: S !< Layer salinities [ppt].
813  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: R0 !< Potential density referenced to
814  !! surface pressure [R ~> kg m-3].
815  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: Rcv !< The coordinate defining potential
816  !! density [R ~> kg m-3].
817  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_eb !< The downward increase across a layer
818  !! in the entrainment from below [H ~> m or kg m-2].
819  !! Positive values go with mass gain by
820  !! a layer.
821  real, dimension(SZI_(G),SZK_(GV)), intent(in) :: eps !< The negligibly small amount of water
822  !! that will be left in each layer [H ~> m or kg m-2].
823  real, dimension(SZI_(G),SZK_(GV)), intent(out) :: dKE_CA !< The vertically integrated change in
824  !! kinetic energy due to convective
825  !! adjustment [Z L2 T-2 ~> m3 s-2].
826  real, dimension(SZI_(G),SZK_(GV)), intent(out) :: cTKE !< The buoyant turbulent kinetic energy
827  !! source due to convective adjustment
828  !! [Z L2 T-2 ~> m3 s-2].
829  integer, intent(in) :: j !< The j-index to work on.
830  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
831  type(bulkmixedlayer_CS), pointer :: CS !< The control structure for this module.
832  integer, optional, intent(in) :: nz_conv !< If present, the number of layers
833  !! over which to do convective adjustment
834  !! (perhaps CS%nkml).
835 
836 ! This subroutine does instantaneous convective entrainment into the buffer
837 ! layers and mixed layers to remove hydrostatic instabilities. Any water that
838 ! is lighter than currently in the mixed- or buffer- layer is entrained.
839 
840  ! Local variables
841  real, dimension(SZI_(G)) :: &
842  htot, & ! The total depth of the layers being considered for
843  ! entrainment [H ~> m or kg m-2].
844  r0_tot, & ! The integrated potential density referenced to the surface
845  ! of the layers which are fully entrained [H R ~> kg m-2 or kg2 m-5].
846  rcv_tot, & ! The integrated coordinate value potential density of the
847  ! layers that are fully entrained [H R ~> kg m-2 or kg2 m-5].
848  ttot, & ! The integrated temperature of layers which are fully
849  ! entrained [degC H ~> degC m or degC kg m-2].
850  stot, & ! The integrated salt of layers which are fully entrained
851  ! [H ppt ~> m ppt or ppt kg m-2].
852  uhtot, & ! The depth integrated zonal and meridional velocities in
853  vhtot, & ! the mixed layer [H L T-1 ~> m2 s-1 or kg m-1 s-1].
854  ke_orig, & ! The total mean kinetic energy in the mixed layer before
855  ! convection, [H L2 T-2 ~> H m2 s-2].
856  h_orig_k1 ! The depth of layer k1 before convective adjustment [H ~> m or kg m-2].
857  real :: h_ent ! The thickness from a layer that is entrained [H ~> m or kg m-2].
858  real :: Ih ! The inverse of a thickness [H-1 ~> m-1 or m2 kg-1].
859  real :: g_H2_2Rho0 ! Half the gravitational acceleration times the square of
860  ! the conversion from H to Z divided by the mean density,
861  ! in [L2 Z T-3 H-2 R-1 ~> m4 s-3 kg-1 or m10 s-3 kg-3].
862  integer :: is, ie, nz, i, k, k1, nzc, nkmb
863 
864  is = g%isc ; ie = g%iec ; nz = gv%ke
865  g_h2_2rho0 = (gv%g_Earth * gv%H_to_Z**2) / (2.0 * gv%Rho0)
866  nzc = nz ; if (present(nz_conv)) nzc = nz_conv
867  nkmb = cs%nkml+cs%nkbl
868 
869 ! Undergo instantaneous entrainment into the buffer layers and mixed layers
870 ! to remove hydrostatic instabilities. Any water that is lighter than currently
871 ! in the layer is entrained.
872  do k1=min(nzc-1,nkmb),1,-1
873  do i=is,ie
874  h_orig_k1(i) = h(i,k1)
875  ke_orig(i) = 0.5*h(i,k1)*(u(i,k1)**2 + v(i,k1)**2)
876  uhtot(i) = h(i,k1)*u(i,k1) ; vhtot(i) = h(i,k1)*v(i,k1)
877  r0_tot(i) = r0(i,k1) * h(i,k1)
878  ctke(i,k1) = 0.0 ; dke_ca(i,k1) = 0.0
879 
880  rcv_tot(i) = rcv(i,k1) * h(i,k1)
881  ttot(i) = t(i,k1) * h(i,k1) ; stot(i) = s(i,k1) * h(i,k1)
882  enddo
883  do k=k1+1,nzc
884  do i=is,ie
885  if ((h(i,k) > eps(i,k)) .and. (r0_tot(i) > h(i,k1)*r0(i,k))) then
886  h_ent = h(i,k)-eps(i,k)
887  ctke(i,k1) = ctke(i,k1) + h_ent * g_h2_2rho0 * &
888  (r0_tot(i) - h(i,k1)*r0(i,k)) * cs%nstar2
889  if (k < nkmb) then
890  ctke(i,k1) = ctke(i,k1) + ctke(i,k)
891  dke_ca(i,k1) = dke_ca(i,k1) + dke_ca(i,k)
892  endif
893  r0_tot(i) = r0_tot(i) + h_ent * r0(i,k)
894  ke_orig(i) = ke_orig(i) + 0.5*h_ent* &
895  (u(i,k)*u(i,k) + v(i,k)*v(i,k))
896  uhtot(i) = uhtot(i) + h_ent*u(i,k)
897  vhtot(i) = vhtot(i) + h_ent*v(i,k)
898 
899  rcv_tot(i) = rcv_tot(i) + h_ent * rcv(i,k)
900  ttot(i) = ttot(i) + h_ent * t(i,k)
901  stot(i) = stot(i) + h_ent * s(i,k)
902  h(i,k1) = h(i,k1) + h_ent ; h(i,k) = eps(i,k)
903 
904  d_eb(i,k) = d_eb(i,k) - h_ent
905  d_eb(i,k1) = d_eb(i,k1) + h_ent
906  endif
907  enddo
908  enddo
909 ! Determine the temperature, salinity, and velocities of the mixed or buffer
910 ! layer in question, if it has entrained.
911  do i=is,ie ; if (h(i,k1) > h_orig_k1(i)) then
912  ih = 1.0 / h(i,k1)
913  r0(i,k1) = r0_tot(i) * ih
914  u(i,k1) = uhtot(i) * ih ; v(i,k1) = vhtot(i) * ih
915  dke_ca(i,k1) = dke_ca(i,k1) + gv%H_to_Z * (cs%bulk_Ri_convective * &
916  (ke_orig(i) - 0.5*h(i,k1)*(u(i,k1)**2 + v(i,k1)**2)))
917  rcv(i,k1) = rcv_tot(i) * ih
918  t(i,k1) = ttot(i) * ih ; s(i,k1) = stot(i) * ih
919  endif ; enddo
920  enddo
921 ! If lower mixed or buffer layers are massless, give them the properties of the
922 ! layer above.
923  do k=2,min(nzc,nkmb) ; do i=is,ie ; if (h(i,k) == 0.0) then
924  r0(i,k) = r0(i,k-1)
925  rcv(i,k) = rcv(i,k-1) ; t(i,k) = t(i,k-1) ; s(i,k) = s(i,k-1)
926  endif ; enddo ; enddo
927 

Referenced by bulkmixedlayer().

Here is the caller graph for this function:

◆ ef4()

real function mom_bulk_mixed_layer::ef4 ( real, intent(in)  Ht,
real, intent(in)  En,
real, intent(in)  I_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/((Ht+E)*Ht).

Parameters
[in]htTotal thickness [H ~> m or kg m-2].
[in]enEntrainment [H ~> m or kg m-2].
[in]i_lThe e-folding scale [H-1 ~> m-1 or m2 kg-1]
[in,out]dr_deThe partial derivative of the result R with E [H-2 ~> m-2 or m4 kg-2].
Returns
The integral [H-1 ~> m-1 or m2 kg-1].

Definition at line 3683 of file MOM_bulk_mixed_layer.F90.

3683  real, intent(in) :: Ht !< Total thickness [H ~> m or kg m-2].
3684  real, intent(in) :: En !< Entrainment [H ~> m or kg m-2].
3685  real, intent(in) :: I_L !< The e-folding scale [H-1 ~> m-1 or m2 kg-1]
3686  real, optional, intent(inout) :: dR_de !< The partial derivative of the result R with E [H-2 ~> m-2 or m4 kg-2].
3687  real :: EF4 !< The integral [H-1 ~> m-1 or m2 kg-1].
3688 
3689  ! Local variables
3690  real :: exp_LHpE ! A nondimensional exponential decay.
3691  real :: I_HpE ! An inverse thickness plus entrainment [H-1 ~> m-1 or m2 kg-1].
3692  real :: Res ! The result of the integral above [H-1 ~> m-1 or m2 kg-1].
3693 
3694  exp_lhpe = exp(-i_l*(en+ht))
3695  i_hpe = 1.0/(ht+en)
3696  res = exp_lhpe * (en*i_hpe/ht - 0.5*i_l*log(ht*i_hpe) + 0.5*i_l*i_l*en)
3697  if (PRESENT(dr_de)) &
3698  dr_de = -i_l*res + exp_lhpe*(i_hpe*i_hpe + 0.5*i_l*i_hpe + 0.5*i_l*i_l)
3699  ef4 = res
3700 

Referenced by mechanical_entrainment().

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(unit_scale_type), intent(in)  US,
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]usA dimensional unit scaling type
[in]htotThe accumulated mixed layer thickness [H ~> m or kg m-2]
[in]h_caThe mixed layer depth after convective adjustment [H ~> m or kg m-2].
[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 [Z L2 T-2 ~> m3 s-2].
[in]dke_fcThe vertically integrated change in kinetic energy due to free convection [Z L2 T-2 ~> 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 [Z m2 T-2 ~> m3 s-2].
[out]idecay_len_tkeThe inverse of the vertical decay scale for TKE [H-1 ~> m-1 or m2 kg-1].
[in]tke_riverThe source of turbulent kinetic energy available for driving mixing at river mouths [Z L2 T-3 ~> m3 s-3].
[out]cmkeCoefficients of HpE and HpE^2 in calculating the denominator of MKE_rate, [H-1 ~> m-1 or m2 kg-1] and [H-2 ~> m-2 or m4 kg-2].
[in]dtThe time step [T ~> s].
[in]idt_diagThe inverse of the accumulated diagnostic time interval [T-1 ~> s-1].
[in]jThe j-index to work on.
[in]ksortThe density-sorted k-indicies.
csThe control structure for this module.

Definition at line 1308 of file MOM_bulk_mixed_layer.F90.

1308  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1309  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1310  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1311  real, dimension(SZI_(G)), intent(in) :: htot !< The accumulated mixed layer thickness
1312  !! [H ~> m or kg m-2]
1313  real, dimension(SZI_(G)), intent(in) :: h_CA !< The mixed layer depth after convective
1314  !! adjustment [H ~> m or kg m-2].
1315  type(forcing), intent(in) :: fluxes !< A structure containing pointers to any
1316  !! possible forcing fields. Unused fields
1317  !! have NULL ptrs.
1318  real, dimension(SZI_(G)), intent(inout) :: Conv_En !< The buoyant turbulent kinetic energy source
1319  !! due to free convection [Z L2 T-2 ~> m3 s-2].
1320  real, dimension(SZI_(G)), intent(in) :: dKE_FC !< The vertically integrated change in
1321  !! kinetic energy due to free convection
1322  !! [Z L2 T-2 ~> m3 s-2].
1323  real, dimension(SZI_(G),SZK_(GV)), &
1324  intent(in) :: cTKE !< The buoyant turbulent kinetic energy
1325  !! source due to convective adjustment
1326  !! [Z L2 T-2 ~> m3 s-2].
1327  real, dimension(SZI_(G),SZK_(GV)), &
1328  intent(in) :: dKE_CA !< The vertically integrated change in
1329  !! kinetic energy due to convective
1330  !! adjustment [Z L2 T-2 ~> m3 s-2].
1331  real, dimension(SZI_(G)), intent(out) :: TKE !< The turbulent kinetic energy available for
1332  !! mixing over a time step [Z m2 T-2 ~> m3 s-2].
1333  real, dimension(SZI_(G)), intent(out) :: Idecay_len_TKE !< The inverse of the vertical decay
1334  !! scale for TKE [H-1 ~> m-1 or m2 kg-1].
1335  real, dimension(SZI_(G)), intent(in) :: TKE_river !< The source of turbulent kinetic energy
1336  !! available for driving mixing at river mouths
1337  !! [Z L2 T-3 ~> m3 s-3].
1338  real, dimension(2,SZI_(G)), intent(out) :: cMKE !< Coefficients of HpE and HpE^2 in
1339  !! calculating the denominator of MKE_rate,
1340  !! [H-1 ~> m-1 or m2 kg-1] and [H-2 ~> m-2 or m4 kg-2].
1341  real, intent(in) :: dt !< The time step [T ~> s].
1342  real, intent(in) :: Idt_diag !< The inverse of the accumulated diagnostic
1343  !! time interval [T-1 ~> s-1].
1344  integer, intent(in) :: j !< The j-index to work on.
1345  integer, dimension(SZI_(G),SZK_(GV)), &
1346  intent(in) :: ksort !< The density-sorted k-indicies.
1347  type(bulkmixedlayer_CS), pointer :: CS !< The control structure for this module.
1348 
1349 ! This subroutine determines the TKE available at the depth of free
1350 ! convection to drive mechanical entrainment.
1351 
1352  ! Local variables
1353  real :: dKE_conv ! The change in mean kinetic energy due to all convection [Z L2 T-2 ~> m3 s-2].
1354  real :: nstar_FC ! The effective efficiency with which the energy released by
1355  ! free convection is converted to TKE, often ~0.2 [nondim].
1356  real :: nstar_CA ! The effective efficiency with which the energy released by
1357  ! convective adjustment is converted to TKE, often ~0.2 [nondim].
1358  real :: TKE_CA ! The potential energy released by convective adjustment if
1359  ! that release is positive [Z L2 T-2 ~> m3 s-2].
1360  real :: MKE_rate_CA ! MKE_rate for convective adjustment [nondim], 0 to 1.
1361  real :: MKE_rate_FC ! MKE_rate for free convection [nondim], 0 to 1.
1362  real :: totEn_Z ! The total potential energy released by convection, [Z3 T-2 ~> m3 s-2].
1363  real :: Ih ! The inverse of a thickness [H-1 ~> m-1 or m2 kg-1].
1364  real :: exp_kh ! The nondimensional decay of TKE across a layer [nondim].
1365  real :: absf ! The absolute value of f averaged to thickness points [T-1 ~> s-1].
1366  real :: U_star ! The friction velocity [Z T-1 ~> m s-1].
1367  real :: absf_Ustar ! The absolute value of f divided by U_star [Z-1 ~> m-1].
1368  real :: wind_TKE_src ! The surface wind source of TKE [Z L2 T-3 ~> m3 s-3].
1369  real :: diag_wt ! The ratio of the current timestep to the diagnostic
1370  ! timestep (which may include 2 calls) [nondim].
1371  integer :: is, ie, nz, i
1372 
1373  is = g%isc ; ie = g%iec ; nz = gv%ke
1374  diag_wt = dt * idt_diag
1375 
1376  if (cs%omega_frac >= 1.0) absf = 2.0*cs%omega
1377  do i=is,ie
1378  u_star = fluxes%ustar(i,j)
1379  if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then
1380  if (fluxes%frac_shelf_h(i,j) > 0.0) &
1381  u_star = (1.0 - fluxes%frac_shelf_h(i,j)) * u_star + &
1382  fluxes%frac_shelf_h(i,j) * fluxes%ustar_shelf(i,j)
1383  endif
1384 
1385  if (u_star < cs%ustar_min) u_star = cs%ustar_min
1386  if (cs%omega_frac < 1.0) then
1387  absf = 0.25*((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
1388  (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i-1,j))))
1389  if (cs%omega_frac > 0.0) &
1390  absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
1391  endif
1392  absf_ustar = absf / u_star
1393  idecay_len_tke(i) = (absf_ustar * cs%TKE_decay) * gv%H_to_Z
1394 
1395 ! The first number in the denominator could be anywhere up to 16. The
1396 ! value of 3 was chosen to minimize the time-step dependence of the amount
1397 ! of shear-driven mixing in 10 days of a 1-degree global model, emphasizing
1398 ! the equatorial areas. Although it is not cast as a parameter, it should
1399 ! be considered an empirical parameter, and it might depend strongly on the
1400 ! number of sublayers in the mixed layer and their locations.
1401 ! The 0.41 is VonKarman's constant. This equation assumes that small & large
1402 ! scales contribute to mixed layer deepening at similar rates, even though
1403 ! small scales are dissipated more rapidly (implying they are less efficient).
1404 ! Ih = 1.0/(16.0*0.41*U_star*dt)
1405  ih = gv%H_to_Z/(3.0*0.41*u_star*dt)
1406  cmke(1,i) = 4.0 * ih ; cmke(2,i) = (absf_ustar*gv%H_to_Z) * ih
1407 
1408  if (idecay_len_tke(i) > 0.0) then
1409  exp_kh = exp(-htot(i)*idecay_len_tke(i))
1410  else
1411  exp_kh = 1.0
1412  endif
1413 
1414 ! Here nstar is a function of the natural Rossby number 0.2/(1+0.2/Ro), based
1415 ! on a curve fit from the data of Wang (GRL, 2003).
1416 ! Note: Ro = 1.0/sqrt(0.5 * dt * (absf*htot(i))**3 / totEn)
1417  if (conv_en(i) < 0.0) conv_en(i) = 0.0
1418  if (ctke(i,1) > 0.0) then ; tke_ca = ctke(i,1) ; else ; tke_ca = 0.0 ; endif
1419  if ((htot(i) >= h_ca(i)) .or. (tke_ca == 0.0)) then
1420  toten_z = us%L_to_Z**2 * (conv_en(i) + tke_ca)
1421 
1422  if (toten_z > 0.0) then
1423  nstar_fc = cs%nstar * toten_z / (toten_z + 0.2 * &
1424  sqrt(0.5 * dt * (absf*(htot(i)*gv%H_to_Z))**3 * toten_z))
1425  else
1426  nstar_fc = cs%nstar
1427  endif
1428  nstar_ca = nstar_fc
1429  else
1430  ! This reconstructs the Buoyancy flux within the topmost htot of water.
1431  if (conv_en(i) > 0.0) then
1432  toten_z = us%L_to_Z**2 * (conv_en(i) + tke_ca * (htot(i) / h_ca(i)) )
1433  nstar_fc = cs%nstar * toten_z / (toten_z + 0.2 * &
1434  sqrt(0.5 * dt * (absf*(htot(i)*gv%H_to_Z))**3 * toten_z))
1435  else
1436  nstar_fc = cs%nstar
1437  endif
1438 
1439  toten_z = us%L_to_Z**2 * (conv_en(i) + tke_ca)
1440  if (tke_ca > 0.0) then
1441  nstar_ca = cs%nstar * toten_z / (toten_z + 0.2 * &
1442  sqrt(0.5 * dt * (absf*(h_ca(i)*gv%H_to_Z))**3 * toten_z))
1443  else
1444  nstar_ca = cs%nstar
1445  endif
1446  endif
1447 
1448  if (dke_fc(i) + dke_ca(i,1) > 0.0) then
1449  if (htot(i) >= h_ca(i)) then
1450  mke_rate_fc = 1.0 / (1.0 + htot(i)*(cmke(1,i) + cmke(2,i)*htot(i)) )
1451  mke_rate_ca = mke_rate_fc
1452  else
1453  mke_rate_fc = 1.0 / (1.0 + htot(i)*(cmke(1,i) + cmke(2,i)*htot(i)) )
1454  mke_rate_ca = 1.0 / (1.0 + h_ca(i)*(cmke(1,i) + cmke(2,i)*h_ca(i)) )
1455  endif
1456  else
1457  ! This branch just saves unnecessary calculations.
1458  mke_rate_fc = 1.0 ; mke_rate_ca = 1.0
1459  endif
1460 
1461  dke_conv = dke_ca(i,1) * mke_rate_ca + dke_fc(i) * mke_rate_fc
1462 ! At this point, it is assumed that cTKE is positive and stored in TKE_CA!
1463 ! Note: Removed factor of 2 in u*^3 terms.
1464  tke(i) = (dt*cs%mstar)*((us%Z_to_L**2*(u_star*u_star*u_star))*exp_kh) + &
1465  (exp_kh * dke_conv + nstar_fc*conv_en(i) + nstar_ca * tke_ca)
1466 
1467  if (cs%do_rivermix) then ! Add additional TKE at river mouths
1468  tke(i) = tke(i) + tke_river(i)*dt*exp_kh
1469  endif
1470 
1471  if (cs%TKE_diagnostics) then
1472  wind_tke_src = cs%mstar*(us%Z_to_L**2*u_star*u_star*u_star) * diag_wt
1473  cs%diag_TKE_wind(i,j) = cs%diag_TKE_wind(i,j) + &
1474  ( wind_tke_src + tke_river(i) * diag_wt )
1475  cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + dke_conv*idt_diag
1476  cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
1477  (exp_kh-1.0)*(wind_tke_src + dke_conv*idt_diag)
1478  cs%diag_TKE_conv(i,j) = cs%diag_TKE_conv(i,j) + &
1479  idt_diag * (nstar_fc*conv_en(i) + nstar_ca*tke_ca)
1480  cs%diag_TKE_conv_decay(i,j) = cs%diag_TKE_conv_decay(i,j) + &
1481  idt_diag * ((cs%nstar-nstar_fc)*conv_en(i) + (cs%nstar-nstar_ca)*tke_ca)
1482  cs%diag_TKE_conv_s2(i,j) = cs%diag_TKE_conv_s2(i,j) + &
1483  idt_diag * (ctke(i,1)-tke_ca)
1484  endif
1485  enddo
1486 

Referenced by bulkmixedlayer().

Here is the caller graph for this function:

◆ mechanical_entrainment()

subroutine mom_bulk_mixed_layer::mechanical_entrainment ( real, dimension( g %isd: g %ied, gv %ke), intent(inout)  h,
real, dimension( g %isd: g %ied, gv %ke), intent(inout)  d_eb,
real, dimension( g %isd: g %ied), intent(inout)  htot,
real, dimension( g %isd: g %ied), intent(inout)  Ttot,
real, dimension( g %isd: g %ied), intent(inout)  Stot,
real, dimension( g %isd: g %ied), intent(inout)  uhtot,
real, dimension( g %isd: g %ied), intent(inout)  vhtot,
real, dimension( g %isd: g %ied), intent(inout)  R0_tot,
real, dimension( g %isd: g %ied), intent(inout)  Rcv_tot,
real, dimension( g %isd: g %ied, gv %ke), intent(in)  u,
real, dimension( g %isd: g %ied, gv %ke), intent(in)  v,
real, dimension( g %isd: g %ied, gv %ke), intent(in)  T,
real, dimension( g %isd: g %ied, gv %ke), intent(in)  S,
real, dimension( g %isd: g %ied, gv %ke), intent(in)  R0,
real, dimension( g %isd: g %ied, gv %ke), intent(in)  Rcv,
real, dimension( g %isd: g %ied, gv %ke), intent(in)  eps,
real, dimension( g %isd: g %ied), intent(in)  dR0_dT,
real, dimension( g %isd: g %ied), intent(in)  dRcv_dT,
real, dimension(2, g %isd: g %ied), intent(in)  cMKE,
real, intent(in)  Idt_diag,
integer, intent(in)  nsw,
real, dimension(max(nsw,1), g %isd: g %ied), intent(inout)  Pen_SW_bnd,
real, dimension(max(nsw,1), g %isd: g %ied, gv %ke), intent(in)  opacity_band,
real, dimension( g %isd: g %ied), intent(inout)  TKE,
real, dimension( g %isd: g %ied), intent(inout)  Idecay_len_TKE,
integer, intent(in)  j,
integer, dimension( g %isd: g %ied, gv %ke), intent(in)  ksort,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
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]usA dimensional unit scaling type
[in,out]hLayer thickness [H ~> m or kg m-2].
[in,out]d_ebThe downward increase across a layer in the
[in,out]htotThe accumlated mixed layer thickness [H ~> m or kg m-2].
[in,out]ttotThe depth integrated mixed layer temperature [degC H ~> degC m or degC kg m-2].
[in,out]stotThe depth integrated mixed layer salinity [ppt H ~> ppt m or ppt kg m-2].
[in,out]uhtotThe depth integrated mixed layer zonal velocity [H L T-1 ~> m2 s-1 or kg m-1 s-1].
[in,out]vhtotThe integrated mixed layer meridional velocity [H L T-1 ~> m2 s-1 or kg m-1 s-1].
[in,out]r0_totThe integrated mixed layer potential density referenced to 0 pressure [H R ~> kg m-2 or kg2 m-5].
[in,out]rcv_totThe integrated mixed layer coordinate variable potential density [H R ~> kg m-2 or kg2 m-5].
[in]uZonal velocities interpolated to h points [L T-1 ~> m s-1].
[in]vZonal velocities interpolated to h points [L T-1 ~> m s-1].
[in]tLayer temperatures [degC].
[in]sLayer salinities [ppt].
[in]r0Potential density referenced to
[in]rcvThe coordinate defining potential
[in]epsThe negligibly small amount of water
[in]dr0_dtThe partial derivative of R0 with respect to temperature [R degC-1 ~> kg m-3 degC-1].
[in]drcv_dtThe partial derivative of Rcv with respect to temperature [R degC-1 ~> kg m-3 degC-1].
[in]cmkeCoefficients of HpE and HpE^2 used in calculating the denominator of MKE_rate; the two elements have differing units of [H-1 ~> m-1 or m2 kg-1] and [H-2 ~> m-2 or m4 kg-2].
[in]idt_diagThe inverse of the accumulated diagnostic time interval [T-1 ~> s-1].
[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 [degC H ~> degC m or degC kg m-2].
[in]opacity_bandThe opacity in each band of penetrating shortwave radiation [H-1 ~> m-1 or m2 kg-1].
[in,out]tkeThe turbulent kinetic energy available for mixing over a time step [Z m2 T-2 ~> m3 s-2].
[in,out]idecay_len_tkeThe vertical TKE decay rate [H-1 ~> m-1 or m2 kg-1].
[in]jThe j-index to work on.
[in]ksortThe density-sorted k-indicies.
csThe control structure for this module.

Definition at line 1495 of file MOM_bulk_mixed_layer.F90.

1495  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1496  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1497  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1498  real, dimension(SZI_(G),SZK_(GV)), &
1499  intent(inout) :: h !< Layer thickness [H ~> m or kg m-2].
1500  real, dimension(SZI_(G),SZK_(GV)), &
1501  intent(inout) :: d_eb !< The downward increase across a layer in the
1502  !! layer in the entrainment from below [H ~> m or kg m-2].
1503  !! Positive values go with mass gain by a layer.
1504  real, dimension(SZI_(G)), intent(inout) :: htot !< The accumlated mixed layer thickness [H ~> m or kg m-2].
1505  real, dimension(SZI_(G)), intent(inout) :: Ttot !< The depth integrated mixed layer temperature
1506  !! [degC H ~> degC m or degC kg m-2].
1507  real, dimension(SZI_(G)), intent(inout) :: Stot !< The depth integrated mixed layer salinity
1508  !! [ppt H ~> ppt m or ppt kg m-2].
1509  real, dimension(SZI_(G)), intent(inout) :: uhtot !< The depth integrated mixed layer zonal
1510  !! velocity [H L T-1 ~> m2 s-1 or kg m-1 s-1].
1511  real, dimension(SZI_(G)), intent(inout) :: vhtot !< The integrated mixed layer meridional
1512  !! velocity [H L T-1 ~> m2 s-1 or kg m-1 s-1].
1513  real, dimension(SZI_(G)), intent(inout) :: R0_tot !< The integrated mixed layer potential density
1514  !! referenced to 0 pressure [H R ~> kg m-2 or kg2 m-5].
1515  real, dimension(SZI_(G)), intent(inout) :: Rcv_tot !< The integrated mixed layer coordinate variable
1516  !! potential density [H R ~> kg m-2 or kg2 m-5].
1517  real, dimension(SZI_(G),SZK_(GV)), &
1518  intent(in) :: u !< Zonal velocities interpolated to h points [L T-1 ~> m s-1].
1519  real, dimension(SZI_(G),SZK_(GV)), &
1520  intent(in) :: v !< Zonal velocities interpolated to h points [L T-1 ~> m s-1].
1521  real, dimension(SZI_(G),SZK_(GV)), &
1522  intent(in) :: T !< Layer temperatures [degC].
1523  real, dimension(SZI_(G),SZK_(GV)), &
1524  intent(in) :: S !< Layer salinities [ppt].
1525  real, dimension(SZI_(G),SZK_(GV)), &
1526  intent(in) :: R0 !< Potential density referenced to
1527  !! surface pressure [R ~> kg m-3].
1528  real, dimension(SZI_(G),SZK_(GV)), &
1529  intent(in) :: Rcv !< The coordinate defining potential
1530  !! density [R ~> kg m-3].
1531  real, dimension(SZI_(G),SZK_(GV)), &
1532  intent(in) :: eps !< The negligibly small amount of water
1533  !! that will be left in each layer [H ~> m or kg m-2].
1534  real, dimension(SZI_(G)), intent(in) :: dR0_dT !< The partial derivative of R0 with respect to
1535  !! temperature [R degC-1 ~> kg m-3 degC-1].
1536  real, dimension(SZI_(G)), intent(in) :: dRcv_dT !< The partial derivative of Rcv with respect to
1537  !! temperature [R degC-1 ~> kg m-3 degC-1].
1538  real, dimension(2,SZI_(G)), intent(in) :: cMKE !< Coefficients of HpE and HpE^2 used in calculating the
1539  !! denominator of MKE_rate; the two elements have differing
1540  !! units of [H-1 ~> m-1 or m2 kg-1] and [H-2 ~> m-2 or m4 kg-2].
1541  real, intent(in) :: Idt_diag !< The inverse of the accumulated diagnostic
1542  !! time interval [T-1 ~> s-1].
1543  integer, intent(in) :: nsw !< The number of bands of penetrating
1544  !! shortwave radiation.
1545  real, dimension(max(nsw,1),SZI_(G)), intent(inout) :: Pen_SW_bnd !< The penetrating shortwave
1546  !! heating at the sea surface in each penetrating
1547  !! band [degC H ~> degC m or degC kg m-2].
1548  real, dimension(max(nsw,1),SZI_(G),SZK_(GV)), intent(in) :: opacity_band !< The opacity in each band of
1549  !! penetrating shortwave radiation [H-1 ~> m-1 or m2 kg-1].
1550  real, dimension(SZI_(G)), intent(inout) :: TKE !< The turbulent kinetic energy
1551  !! available for mixing over a time
1552  !! step [Z m2 T-2 ~> m3 s-2].
1553  real, dimension(SZI_(G)), intent(inout) :: Idecay_len_TKE !< The vertical TKE decay rate [H-1 ~> m-1 or m2 kg-1].
1554  integer, intent(in) :: j !< The j-index to work on.
1555  integer, dimension(SZI_(G),SZK_(GV)), &
1556  intent(in) :: ksort !< The density-sorted k-indicies.
1557  type(bulkmixedlayer_CS), pointer :: CS !< The control structure for this module.
1558 
1559 ! This subroutine calculates mechanically driven entrainment.
1560 
1561  ! Local variables
1562  real :: SW_trans ! The fraction of shortwave radiation that is not
1563  ! absorbed in a layer, nondimensional.
1564  real :: Pen_absorbed ! The amount of penetrative shortwave radiation
1565  ! that is absorbed in a layer [degC H ~> degC m or degC kg m-2].
1566  real :: h_avail ! The thickness in a layer available for entrainment [H ~> m or kg m-2].
1567  real :: h_ent ! The thickness from a layer that is entrained [H ~> m or kg m-2].
1568  real :: h_min, h_max ! Limits on the solution for h_ent [H ~> m or kg m-2].
1569  real :: dh_Newt ! The Newton's method estimate of the change in
1570  ! h_ent between iterations [H ~> m or kg m-2].
1571  real :: MKE_rate ! The fraction of the energy in resolved shears
1572  ! within the mixed layer that will be eliminated
1573  ! within a timestep, nondim, 0 to 1.
1574  real :: HpE ! The current thickness plus entrainment [H ~> m or kg m-2].
1575  real :: g_H_2Rho0 ! Half the gravitational acceleration times the
1576  ! conversion from H to m divided by the mean density,
1577  ! in [L2 T-2 H-1 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2].
1578  real :: TKE_full_ent ! The TKE remaining if a layer is fully entrained
1579  ! [Z L2 T-2 ~> m3 s-2].
1580  real :: dRL ! Work required to mix water from the next layer
1581  ! across the mixed layer [L2 T-2 ~> L2 s-2].
1582  real :: Pen_En_Contrib ! Penetrating SW contributions to the changes in
1583  ! TKE, divided by layer thickness in m [L2 T2 ~> m2 s-2].
1584  real :: Cpen1 ! A temporary variable [L2 T-2 ~> m2 s-2].
1585  real :: dMKE ! A temporary variable related to the release of mean
1586  ! kinetic energy [H Z L2 T-2 ~> m4 s-2 or kg m s-2]
1587  real :: TKE_ent ! The TKE that remains if h_ent were entrained [Z L2 T-2 ~> m3 s-2].
1588  real :: TKE_ent1 ! The TKE that would remain, without considering the
1589  ! release of mean kinetic energy [Z L2 T-2 ~> m3 s-2].
1590  real :: dTKE_dh ! The partial derivative of TKE with h_ent [Z L2 T-2 H-1 ~> m2 s-2 or m5 s-2 kg-1].
1591  real :: Pen_dTKE_dh_Contrib ! The penetrating shortwave contribution to
1592  ! dTKE_dh [L2 T-2 ~> m2 s-2].
1593  real :: EF4_val ! The result of EF4() (see later) [H-1 ~> m-1 or m2 kg-1].
1594  real :: h_neglect ! A thickness that is so small it is usually lost
1595  ! in roundoff and can be neglected [H ~> m or kg m-2].
1596  real :: dEF4_dh ! The partial derivative of EF4 with h [H-2 ~> m-2 or m4 kg-2].
1597  real :: Pen_En1 ! A nondimensional temporary variable.
1598  real :: kh, exp_kh ! Nondimensional temporary variables related to the
1599  real :: f1_kh ! fractional decay of TKE across a layer.
1600  real :: x1, e_x1 ! Nondimensional temporary variables related to
1601  real :: f1_x1, f2_x1 ! the relative decay of TKE and SW radiation across
1602  real :: f3_x1 ! a layer, and exponential-related functions of x1.
1603  real :: E_HxHpE ! Entrainment divided by the product of the new and old
1604  ! thicknesses [H-1 ~> m-1 or m2 kg-1].
1605  real :: Hmix_min ! The minimum mixed layer depth [H ~> m or kg m-2].
1606  real :: opacity
1607  real :: C1_3, C1_6, C1_24 ! 1/3, 1/6, and 1/24.
1608  integer :: is, ie, nz, i, k, ks, itt, n
1609 
1610  c1_3 = 1.0/3.0 ; c1_6 = 1.0/6.0 ; c1_24 = 1.0/24.0
1611  g_h_2rho0 = (gv%g_Earth * gv%H_to_Z) / (2.0 * gv%Rho0)
1612  hmix_min = cs%Hmix_min
1613  h_neglect = gv%H_subroundoff
1614  is = g%isc ; ie = g%iec ; nz = gv%ke
1615 
1616  do ks=1,nz
1617 
1618  do i=is,ie ; if (ksort(i,ks) > 0) then
1619  k = ksort(i,ks)
1620 
1621  h_avail = h(i,k) - eps(i,k)
1622  if ((h_avail > 0.) .and. ((tke(i) > 0.) .or. (htot(i) < hmix_min))) then
1623  drl = g_h_2rho0 * (r0(i,k)*htot(i) - r0_tot(i) )
1624  dmke = (gv%H_to_Z * cs%bulk_Ri_ML) * 0.5 * &
1625  ((uhtot(i)-u(i,k)*htot(i))**2 + (vhtot(i)-v(i,k)*htot(i))**2)
1626 
1627 ! Find the TKE that would remain if the entire layer were entrained.
1628  kh = idecay_len_tke(i)*h_avail ; exp_kh = exp(-kh)
1629  if (kh >= 2.0e-5) then ; f1_kh = (1.0-exp_kh) / kh
1630  else ; f1_kh = (1.0 - kh*(0.5 - c1_6*kh)) ; endif
1631 
1632  pen_en_contrib = 0.0
1633  do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1634  opacity = opacity_band(n,i,k)
1635 ! Two different forms are used here to make sure that only negative
1636 ! values are taken into exponentials to avoid excessively large
1637 ! numbers. They are, of course, mathematically identical.
1638  if (idecay_len_tke(i) > opacity) then
1639  x1 = (idecay_len_tke(i) - opacity) * h_avail
1640  if (x1 >= 2.0e-5) then
1641  e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1642  f3_x1 = ((e_x1-(1.0-x1))/(x1*x1))
1643  else
1644  f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1645  f3_x1 = (0.5 - x1*(c1_6 - c1_24*x1))
1646  endif
1647 
1648  pen_en1 = exp(-opacity*h_avail) * &
1649  ((1.0+opacity*htot(i))*f1_x1 + opacity*h_avail*f3_x1)
1650  else
1651  x1 = (opacity - idecay_len_tke(i)) * h_avail
1652  if (x1 >= 2.0e-5) then
1653  e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1654  f2_x1 = ((1.0-(1.0+x1)*e_x1)/(x1*x1))
1655  else
1656  f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1657  f2_x1 = (0.5 - x1*(c1_3 - 0.125*x1))
1658  endif
1659 
1660  pen_en1 = exp_kh * ((1.0+opacity*htot(i))*f1_x1 + &
1661  opacity*h_avail*f2_x1)
1662  endif
1663  pen_en_contrib = pen_en_contrib + &
1664  (g_h_2rho0*dr0_dt(i)*pen_sw_bnd(n,i)) * (pen_en1 - f1_kh)
1665  endif ; enddo
1666 
1667  hpe = htot(i)+h_avail
1668  mke_rate = 1.0/(1.0 + (cmke(1,i)*hpe + cmke(2,i)*hpe**2))
1669  ef4_val = ef4(htot(i)+h_neglect,h_avail,idecay_len_tke(i))
1670  tke_full_ent = (exp_kh*tke(i) - (h_avail*gv%H_to_Z)*(drl*f1_kh + pen_en_contrib)) + &
1671  mke_rate*dmke*ef4_val
1672  if ((tke_full_ent >= 0.0) .or. (h_avail+htot(i) <= hmix_min)) then
1673  ! The layer will be fully entrained.
1674  h_ent = h_avail
1675 
1676  if (cs%TKE_diagnostics) then
1677  e_hxhpe = h_ent / ((htot(i)+h_neglect)*(htot(i)+h_ent+h_neglect))
1678  cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
1679  idt_diag * ((exp_kh-1.0)* tke(i) + (h_ent*gv%H_to_Z)*drl*(1.0-f1_kh) + &
1680  mke_rate*dmke*(ef4_val-e_hxhpe))
1681  cs%diag_TKE_mixing(i,j) = cs%diag_TKE_mixing(i,j) - &
1682  idt_diag*(gv%H_to_Z*h_ent)*drl
1683  cs%diag_TKE_pen_SW(i,j) = cs%diag_TKE_pen_SW(i,j) - &
1684  idt_diag*(gv%H_to_Z*h_ent)*pen_en_contrib
1685  cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + &
1686  idt_diag*mke_rate*dmke*e_hxhpe
1687  endif
1688 
1689  tke(i) = tke_full_ent
1690  !### The minimum TKE value in this line may be problematically small.
1691  if (tke(i) <= 0.0) tke(i) = 1.0e-150*us%m_to_Z*us%m_s_to_L_T**2
1692  else
1693 ! The layer is only partially entrained. The amount that will be
1694 ! entrained is determined iteratively. No further layers will be
1695 ! entrained.
1696  h_min = 0.0 ; h_max = h_avail
1697  if (tke(i) <= 0.0) then
1698  h_ent = 0.0
1699  else
1700  h_ent = h_avail * tke(i) / (tke(i) - tke_full_ent)
1701 
1702  do itt=1,15
1703  ! Evaluate the TKE that would remain if h_ent were entrained.
1704 
1705  kh = idecay_len_tke(i)*h_ent ; exp_kh = exp(-kh)
1706  if (kh >= 2.0e-5) then
1707  f1_kh = (1.0-exp_kh) / kh
1708  else
1709  f1_kh = (1.0 - kh*(0.5 - c1_6*kh))
1710  endif
1711 
1712 
1713  pen_en_contrib = 0.0 ; pen_dtke_dh_contrib = 0.0
1714  do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1715  ! Two different forms are used here to make sure that only negative
1716  ! values are taken into exponentials to avoid excessively large
1717  ! numbers. They are, of course, mathematically identical.
1718  opacity = opacity_band(n,i,k)
1719  sw_trans = exp(-h_ent*opacity)
1720  if (idecay_len_tke(i) > opacity) then
1721  x1 = (idecay_len_tke(i) - opacity) * h_ent
1722  if (x1 >= 2.0e-5) then
1723  e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1724  f3_x1 = ((e_x1-(1.0-x1))/(x1*x1))
1725  else
1726  f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1727  f3_x1 = (0.5 - x1*(c1_6 - c1_24*x1))
1728  endif
1729  pen_en1 = sw_trans * ((1.0+opacity*htot(i))*f1_x1 + &
1730  opacity*h_ent*f3_x1)
1731  else
1732  x1 = (opacity - idecay_len_tke(i)) * h_ent
1733  if (x1 >= 2.0e-5) then
1734  e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
1735  f2_x1 = ((1.0-(1.0+x1)*e_x1)/(x1*x1))
1736  else
1737  f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
1738  f2_x1 = (0.5 - x1*(c1_3 - 0.125*x1))
1739  endif
1740 
1741  pen_en1 = exp_kh * ((1.0+opacity*htot(i))*f1_x1 + &
1742  opacity*h_ent*f2_x1)
1743  endif
1744  cpen1 = g_h_2rho0*dr0_dt(i)*pen_sw_bnd(n,i)
1745  pen_en_contrib = pen_en_contrib + cpen1*(pen_en1 - f1_kh)
1746  pen_dtke_dh_contrib = pen_dtke_dh_contrib + &
1747  cpen1*((1.0-sw_trans) - opacity*(htot(i) + h_ent)*sw_trans)
1748  endif ; enddo ! (Pen_SW_bnd(n,i) > 0.0)
1749 
1750  tke_ent1 = exp_kh* tke(i) - (h_ent*gv%H_to_Z)*(drl*f1_kh + pen_en_contrib)
1751  ef4_val = ef4(htot(i)+h_neglect,h_ent,idecay_len_tke(i),def4_dh)
1752  hpe = htot(i)+h_ent
1753  mke_rate = 1.0/(1.0 + (cmke(1,i)*hpe + cmke(2,i)*hpe**2))
1754  tke_ent = tke_ent1 + dmke*ef4_val*mke_rate
1755  ! TKE_ent is the TKE that would remain if h_ent were entrained.
1756 
1757  dtke_dh = ((-idecay_len_tke(i)*tke_ent1 - drl*gv%H_to_Z) + &
1758  pen_dtke_dh_contrib*gv%H_to_Z) + dmke * mke_rate* &
1759  (def4_dh - ef4_val*mke_rate*(cmke(1,i)+2.0*cmke(2,i)*hpe))
1760  ! dh_Newt = -TKE_ent / dTKE_dh
1761  ! Bisect if the Newton's method prediction is outside of the bounded range.
1762  if (tke_ent > 0.0) then
1763  if ((h_max-h_ent)*(-dtke_dh) > tke_ent) then
1764  dh_newt = -tke_ent / dtke_dh
1765  else
1766  dh_newt = 0.5*(h_max-h_ent)
1767  endif
1768  h_min = h_ent
1769  else
1770  if ((h_min-h_ent)*(-dtke_dh) < tke_ent) then
1771  dh_newt = -tke_ent / dtke_dh
1772  else
1773  dh_newt = 0.5*(h_min-h_ent)
1774  endif
1775  h_max = h_ent
1776  endif
1777  h_ent = h_ent + dh_newt
1778 
1779  if (abs(dh_newt) < 0.2*gv%Angstrom_H) exit
1780  enddo
1781  endif
1782 
1783  if (h_ent < hmix_min-htot(i)) h_ent = hmix_min - htot(i)
1784 
1785  if (cs%TKE_diagnostics) then
1786  hpe = htot(i)+h_ent
1787  mke_rate = 1.0/(1.0 + cmke(1,i)*hpe + cmke(2,i)*hpe**2)
1788  ef4_val = ef4(htot(i)+h_neglect,h_ent,idecay_len_tke(i))
1789 
1790  e_hxhpe = h_ent / ((htot(i)+h_neglect)*(hpe+h_neglect))
1791  cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
1792  idt_diag * ((exp_kh-1.0)* tke(i) + (h_ent*gv%H_to_Z)*drl*(1.0-f1_kh) + &
1793  dmke*mke_rate*(ef4_val-e_hxhpe))
1794  cs%diag_TKE_mixing(i,j) = cs%diag_TKE_mixing(i,j) - &
1795  idt_diag*(h_ent*gv%H_to_Z)*drl
1796  cs%diag_TKE_pen_SW(i,j) = cs%diag_TKE_pen_SW(i,j) - &
1797  idt_diag*(h_ent*gv%H_to_Z)*pen_en_contrib
1798  cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + &
1799  idt_diag*dmke*mke_rate*e_hxhpe
1800  endif
1801 
1802  tke(i) = 0.0
1803  endif ! TKE_full_ent > 0.0
1804 
1805  pen_absorbed = 0.0
1806  do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1807  sw_trans = exp(-h_ent*opacity_band(n,i,k))
1808  pen_absorbed = pen_absorbed + pen_sw_bnd(n,i) * (1.0 - sw_trans)
1809  pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
1810  endif ; enddo
1811 
1812  htot(i) = htot(i) + h_ent
1813  r0_tot(i) = r0_tot(i) + (h_ent * r0(i,k) + pen_absorbed*dr0_dt(i))
1814  h(i,k) = h(i,k) - h_ent
1815  d_eb(i,k) = d_eb(i,k) - h_ent
1816 
1817  stot(i) = stot(i) + h_ent * s(i,k)
1818  ttot(i) = ttot(i) + (h_ent * t(i,k) + pen_absorbed)
1819  rcv_tot(i) = rcv_tot(i) + (h_ent*rcv(i,k) + pen_absorbed*drcv_dt(i))
1820 
1821  uhtot(i) = uhtot(i) + u(i,k)*h_ent
1822  vhtot(i) = vhtot(i) + v(i,k)*h_ent
1823  endif ! h_avail > 0.0 .AND TKE(i) > 0.0
1824 
1825  endif ; enddo ! i loop
1826  enddo ! k loop
1827 

References ef4().

Referenced by bulkmixedlayer().

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( g %isd: g %ied, gv %ke), intent(inout)  h,
real, dimension( g %isd: g %ied, gv %ke), intent(inout)  d_eb,
real, dimension( g %isd: g %ied), intent(out)  htot,
real, dimension( g %isd: g %ied), intent(out)  Ttot,
real, dimension( g %isd: g %ied), intent(out)  Stot,
real, dimension( g %isd: g %ied), intent(out)  uhtot,
real, dimension( g %isd: g %ied), intent(out)  vhtot,
real, dimension( g %isd: g %ied), intent(out)  R0_tot,
real, dimension( g %isd: g %ied), intent(out)  Rcv_tot,
real, dimension( g %isd: g %ied, gv %ke), intent(in)  u,
real, dimension( g %isd: g %ied, gv %ke), intent(in)  v,
real, dimension( g %isd: g %ied, gv %ke), intent(in)  T,
real, dimension( g %isd: g %ied, gv %ke), intent(in)  S,
real, dimension( g %isd: g %ied, gv %ke), intent(in)  R0,
real, dimension( g %isd: g %ied, gv %ke), intent(in)  Rcv,
real, dimension( g %isd: g %ied, gv %ke), intent(in)  eps,
real, dimension( g %isd: g %ied), intent(in)  dR0_dT,
real, dimension( g %isd: g %ied), intent(in)  dRcv_dT,
real, dimension( g %isd: g %ied), intent(in)  dR0_dS,
real, dimension( g %isd: g %ied), intent(in)  dRcv_dS,
real, dimension( g %isd: g %ied), intent(in)  netMassInOut,
real, dimension( g %isd: g %ied), intent(in)  netMassOut,
real, dimension( g %isd: g %ied), intent(in)  Net_heat,
real, dimension( g %isd: g %ied), intent(in)  Net_salt,
integer, intent(in)  nsw,
real, dimension(max(nsw,1), g %isd: g %ied), intent(inout)  Pen_SW_bnd,
real, dimension(max(nsw,1), g %isd: g %ied, gv %ke), intent(in)  opacity_band,
real, dimension( g %isd: g %ied), intent(out)  Conv_En,
real, dimension( g %isd: g %ied), intent(out)  dKE_FC,
integer, intent(in)  j,
integer, dimension( g %isd: g %ied, gv %ke), intent(in)  ksort,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
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 [H ~> m or kg m-2].
[in,out]d_ebThe downward increase across a layer in the
[out]htotThe accumulated mixed layer thickness [H ~> m or kg m-2].
[out]ttotThe depth integrated mixed layer temperature [degC H ~> degC m or degC kg m-2].
[out]stotThe depth integrated mixed layer salinity [ppt H ~> ppt m or ppt kg m-2].
[out]uhtotThe depth integrated mixed layer zonal velocity [H L T-1 ~> m2 s-1 or kg m-1 s-1].
[out]vhtotThe integrated mixed layer meridional velocity [H L T-1 ~> m2 s-1 or kg m-1 s-1].
[out]r0_totThe integrated mixed layer potential density referenced to 0 pressure [H R ~> kg m-2 or kg2 m-5].
[out]rcv_totThe integrated mixed layer coordinate variable potential density [H R ~> kg m-2 or kg2 m-5].
[in]uZonal velocities interpolated to h points [L T-1 ~> m s-1].
[in]vZonal velocities interpolated to h points [L T-1 ~> m s-1].
[in]tLayer temperatures [degC].
[in]sLayer salinities [ppt].
[in]r0Potential density referenced to
[in]rcvThe coordinate defining potential
[in]epsThe negligibly small amount of water
[in]dr0_dtThe partial derivative of R0 with respect to temperature [R degC-1 ~> kg m-3 degC-1].
[in]drcv_dtThe partial derivative of Rcv with respect to temperature [R degC-1 ~> kg m-3 degC-1].
[in]dr0_dsThe partial derivative of R0 with respect to salinity [R ppt-1 ~> kg m-3 ppt-1].
[in]drcv_dsThe partial derivative of Rcv with respect to salinity [R ppt-1 ~> kg m-3 ppt-1].
[in]netmassinoutThe net mass flux (if non-Boussinesq) or volume flux (if Boussinesq) into the ocean within a time step [H ~> m or kg m-2]. (I.e. P+R-E.)
[in]netmassoutThe mass or volume flux out of the ocean within a time step [H ~> m or kg m-2].
[in]net_heatThe net heating at the surface over a time step [degC H ~> degC m or degC kg m-2]. Any penetrating shortwave radiation is not included in Net_heat.
[in]net_saltThe net surface salt flux into the ocean over a time step [ppt H ~> ppt m or ppt 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 [degC H ~> degC m or degC kg m-2].
[in]opacity_bandThe opacity in each band of penetrating shortwave radiation [H-1 ~> m-1 or m2 kg-1].
[out]conv_enThe buoyant turbulent kinetic energy source due to free convection [Z L2 T-2 ~> m3 s-2].
[out]dke_fcThe vertically integrated change in kinetic energy due to free convection [Z L2 T-2 ~> m3 s-2].
[in]jThe j-index to work on.
[in]ksortThe density-sorted k-indices.
[in]usA dimensional unit scaling type
csThe control structure for this module.
[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 [T ~> s].
[in]aggregate_fw_forcingIf true, the net incoming and outgoing surface freshwater fluxes are combined before being applied, instead of being applied separately.

Definition at line 940 of file MOM_bulk_mixed_layer.F90.

940  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
941  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
942  real, dimension(SZI_(G),SZK_(GV)), &
943  intent(inout) :: h !< Layer thickness [H ~> m or kg m-2].
944  !! The units of h are referred to as H below.
945  real, dimension(SZI_(G),SZK_(GV)), &
946  intent(inout) :: d_eb !< The downward increase across a layer in the
947  !! layer in the entrainment from below [H ~> m or kg m-2].
948  !! Positive values go with mass gain by a layer.
949  real, dimension(SZI_(G)), intent(out) :: htot !< The accumulated mixed layer thickness [H ~> m or kg m-2].
950  real, dimension(SZI_(G)), intent(out) :: Ttot !< The depth integrated mixed layer temperature
951  !! [degC H ~> degC m or degC kg m-2].
952  real, dimension(SZI_(G)), intent(out) :: Stot !< The depth integrated mixed layer salinity
953  !! [ppt H ~> ppt m or ppt kg m-2].
954  real, dimension(SZI_(G)), intent(out) :: uhtot !< The depth integrated mixed layer zonal
955  !! velocity [H L T-1 ~> m2 s-1 or kg m-1 s-1].
956  real, dimension(SZI_(G)), intent(out) :: vhtot !< The integrated mixed layer meridional
957  !! velocity [H L T-1 ~> m2 s-1 or kg m-1 s-1].
958  real, dimension(SZI_(G)), intent(out) :: R0_tot !< The integrated mixed layer potential density referenced
959  !! to 0 pressure [H R ~> kg m-2 or kg2 m-5].
960  real, dimension(SZI_(G)), intent(out) :: Rcv_tot !< The integrated mixed layer coordinate
961  !! variable potential density [H R ~> kg m-2 or kg2 m-5].
962  real, dimension(SZI_(G),SZK_(GV)), &
963  intent(in) :: u !< Zonal velocities interpolated to h points [L T-1 ~> m s-1].
964  real, dimension(SZI_(G),SZK_(GV)), &
965  intent(in) :: v !< Zonal velocities interpolated to h points [L T-1 ~> m s-1].
966  real, dimension(SZI_(G),SZK_(GV)), &
967  intent(in) :: T !< Layer temperatures [degC].
968  real, dimension(SZI_(G),SZK_(GV)), &
969  intent(in) :: S !< Layer salinities [ppt].
970  real, dimension(SZI_(G),SZK_(GV)), &
971  intent(in) :: R0 !< Potential density referenced to
972  !! surface pressure [R ~> kg m-3].
973  real, dimension(SZI_(G),SZK_(GV)), &
974  intent(in) :: Rcv !< The coordinate defining potential
975  !! density [R ~> kg m-3].
976  real, dimension(SZI_(G),SZK_(GV)), &
977  intent(in) :: eps !< The negligibly small amount of water
978  !! that will be left in each layer [H ~> m or kg m-2].
979  real, dimension(SZI_(G)), intent(in) :: dR0_dT !< The partial derivative of R0 with respect to
980  !! temperature [R degC-1 ~> kg m-3 degC-1].
981  real, dimension(SZI_(G)), intent(in) :: dRcv_dT !< The partial derivative of Rcv with respect to
982  !! temperature [R degC-1 ~> kg m-3 degC-1].
983  real, dimension(SZI_(G)), intent(in) :: dR0_dS !< The partial derivative of R0 with respect to
984  !! salinity [R ppt-1 ~> kg m-3 ppt-1].
985  real, dimension(SZI_(G)), intent(in) :: dRcv_dS !< The partial derivative of Rcv with respect to
986  !! salinity [R ppt-1 ~> kg m-3 ppt-1].
987  real, dimension(SZI_(G)), intent(in) :: netMassInOut !< The net mass flux (if non-Boussinesq)
988  !! or volume flux (if Boussinesq) into the ocean
989  !! within a time step [H ~> m or kg m-2]. (I.e. P+R-E.)
990  real, dimension(SZI_(G)), intent(in) :: netMassOut !< The mass or volume flux out of the ocean
991  !! within a time step [H ~> m or kg m-2].
992  real, dimension(SZI_(G)), intent(in) :: Net_heat !< The net heating at the surface over a time
993  !! step [degC H ~> degC m or degC kg m-2]. Any penetrating
994  !! shortwave radiation is not included in Net_heat.
995  real, dimension(SZI_(G)), intent(in) :: Net_salt !< The net surface salt flux into the ocean
996  !! over a time step [ppt H ~> ppt m or ppt kg m-2].
997  integer, intent(in) :: nsw !< The number of bands of penetrating
998  !! shortwave radiation.
999  real, dimension(max(nsw,1),SZI_(G)), intent(inout) :: Pen_SW_bnd !< The penetrating shortwave
1000  !! heating at the sea surface in each penetrating
1001  !! band [degC H ~> degC m or degC kg m-2].
1002  real, dimension(max(nsw,1),SZI_(G),SZK_(GV)), intent(in) :: opacity_band !< The opacity in each band of
1003  !! penetrating shortwave radiation [H-1 ~> m-1 or m2 kg-1].
1004  real, dimension(SZI_(G)), intent(out) :: Conv_En !< The buoyant turbulent kinetic energy source
1005  !! due to free convection [Z L2 T-2 ~> m3 s-2].
1006  real, dimension(SZI_(G)), intent(out) :: dKE_FC !< The vertically integrated change in kinetic
1007  !! energy due to free convection [Z L2 T-2 ~> m3 s-2].
1008  integer, intent(in) :: j !< The j-index to work on.
1009  integer, dimension(SZI_(G),SZK_(GV)), &
1010  intent(in) :: ksort !< The density-sorted k-indices.
1011  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1012  type(bulkmixedlayer_CS), pointer :: CS !< The control structure for this module.
1013  type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers to any
1014  !! available thermodynamic fields. Absent
1015  !! fields have NULL ptrs.
1016  type(forcing), intent(inout) :: fluxes !< A structure containing pointers to any
1017  !! possible forcing fields. Unused fields
1018  !! have NULL ptrs.
1019  real, intent(in) :: dt !< Time increment [T ~> s].
1020  logical, intent(in) :: aggregate_FW_forcing !< If true, the net incoming and
1021  !! outgoing surface freshwater fluxes are
1022  !! combined before being applied, instead of
1023  !! being applied separately.
1024 
1025 ! This subroutine causes the mixed layer to entrain to the depth of free
1026 ! convection. The depth of free convection is the shallowest depth at which the
1027 ! fluid is denser than the average of the fluid above.
1028 
1029  ! Local variables
1030  real, dimension(SZI_(G)) :: &
1031  massOutRem, & ! Evaporation that remains to be supplied [H ~> m or kg m-2].
1032  netMassIn ! mass entering through ocean surface [H ~> m or kg m-2]
1033  real :: SW_trans ! The fraction of shortwave radiation
1034  ! that is not absorbed in a layer [nondim].
1035  real :: Pen_absorbed ! The amount of penetrative shortwave radiation
1036  ! that is absorbed in a layer [degC H ~> degC m or degC kg m-2].
1037  real :: h_avail ! The thickness in a layer available for
1038  ! entrainment [H ~> m or kg m-2].
1039  real :: h_ent ! The thickness from a layer that is entrained [H ~> m or kg m-2].
1040  real :: T_precip ! The temperature of the precipitation [degC].
1041  real :: C1_3, C1_6 ! 1/3 and 1/6.
1042  real :: En_fn, Frac, x1 ! Nondimensional temporary variables.
1043  real :: dr, dr0 ! Temporary variables [R H ~> kg m-2 or kg2 m-5].
1044  real :: dr_ent, dr_comp ! Temporary variables [R H ~> kg m-2 or kg2 m-5].
1045  real :: dr_dh ! The partial derivative of dr_ent with h_ent [R ~> kg m-3].
1046  real :: h_min, h_max ! The minimum, maximum, and previous estimates for
1047  real :: h_prev ! h_ent [H ~> m or kg m-2].
1048  real :: h_evap ! The thickness that is evaporated [H ~> m or kg m-2].
1049  real :: dh_Newt ! The Newton's method estimate of the change in
1050  ! h_ent between iterations [H ~> m or kg m-2].
1051  real :: g_H2_2Rho0 ! Half the gravitational acceleration times the square of
1052  ! the conversion from H to Z divided by the mean density,
1053  ! [L2 Z T-3 H-2 R-1 ~> m4 s-3 kg-1 or m10 s-3 kg-3].
1054  real :: Angstrom ! The minimum layer thickness [H ~> m or kg m-2].
1055  real :: opacity ! The opacity converted to inverse thickness units [H-1 ~> m-1 or m2 kg-1]
1056  real :: sum_Pen_En ! The potential energy change due to penetrating
1057  ! shortwave radiation, integrated over a layer
1058  ! [H R ~> kg m-2 or kg2 m-5].
1059  real :: Idt ! 1.0/dt [T-1 ~> s-1]
1060  real :: netHeatOut ! accumulated heat content of mass leaving ocean
1061  integer :: is, ie, nz, i, k, ks, itt, n
1062  real, dimension(max(nsw,1)) :: &
1063  C2, & ! Temporary variable R H-1 ~> kg m-4 or m-1].
1064  r_SW_top ! Temporary variables [H R ~> kg m-2 or kg2 m-5].
1065 
1066  angstrom = gv%Angstrom_H
1067  c1_3 = 1.0/3.0 ; c1_6 = 1.0/6.0
1068  g_h2_2rho0 = (gv%g_Earth * gv%H_to_Z**2) / (2.0 * gv%Rho0)
1069  idt = 1.0 / dt
1070  is = g%isc ; ie = g%iec ; nz = gv%ke
1071 
1072  do i=is,ie ; if (ksort(i,1) > 0) then
1073  k = ksort(i,1)
1074 
1075  if (aggregate_fw_forcing) then
1076  massoutrem(i) = 0.0
1077  if (netmassinout(i) < 0.0) massoutrem(i) = -netmassinout(i)
1078  netmassin(i) = netmassinout(i) + massoutrem(i)
1079  else
1080  massoutrem(i) = -netmassout(i)
1081  netmassin(i) = netmassinout(i) - netmassout(i)
1082  endif
1083 
1084  ! htot is an Angstrom (taken from layer 1) plus any net precipitation.
1085  h_ent = max(min(angstrom,h(i,k)-eps(i,k)),0.0)
1086  htot(i) = h_ent + netmassin(i)
1087  h(i,k) = h(i,k) - h_ent
1088  d_eb(i,k) = d_eb(i,k) - h_ent
1089 
1090  pen_absorbed = 0.0
1091  do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1092  sw_trans = exp(-htot(i)*opacity_band(n,i,k))
1093  pen_absorbed = pen_absorbed + pen_sw_bnd(n,i) * (1.0-sw_trans)
1094  pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
1095  endif ; enddo
1096 
1097  ! Precipitation is assumed to have the same temperature and velocity
1098  ! as layer 1. Because layer 1 might not be the topmost layer, this
1099  ! involves multiple terms.
1100  t_precip = t(i,1)
1101  ttot(i) = (net_heat(i) + (netmassin(i) * t_precip + h_ent * t(i,k))) + &
1102  pen_absorbed
1103  ! Net_heat contains both heat fluxes and the heat content of mass fluxes.
1104  !! Ttot(i) = netMassIn(i) * T_precip + h_ent * T(i,k)
1105  !! Ttot(i) = Net_heat(i) + Ttot(i)
1106  !! Ttot(i) = Ttot(i) + Pen_absorbed
1107  ! smg:
1108  ! Ttot(i) = (Net_heat(i) + (h_ent * T(i,k))) + Pen_absorbed
1109  stot(i) = h_ent*s(i,k) + net_salt(i)
1110  uhtot(i) = u(i,1)*netmassin(i) + u(i,k)*h_ent
1111  vhtot(i) = v(i,1)*netmassin(i) + v(i,k)*h_ent
1112  r0_tot(i) = (h_ent*r0(i,k) + netmassin(i)*r0(i,1)) + &
1113 ! dR0_dT(i)*netMassIn(i)*(T_precip - T(i,1)) + &
1114  (dr0_dt(i)*(net_heat(i) + pen_absorbed) - &
1115  dr0_ds(i) * (netmassin(i) * s(i,1) - net_salt(i)))
1116  rcv_tot(i) = (h_ent*rcv(i,k) + netmassin(i)*rcv(i,1)) + &
1117 ! dRcv_dT(i)*netMassIn(i)*(T_precip - T(i,1)) + &
1118  (drcv_dt(i)*(net_heat(i) + pen_absorbed) - &
1119  drcv_ds(i) * (netmassin(i) * s(i,1) - net_salt(i)))
1120  conv_en(i) = 0.0 ; dke_fc(i) = 0.0
1121  if (associated(fluxes%heat_content_massin)) &
1122  fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + &
1123  t_precip * netmassin(i) * gv%H_to_RZ * fluxes%C_p * idt
1124  if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
1125  t_precip * netmassin(i) * gv%H_to_RZ
1126  endif ; enddo
1127 
1128  ! Now do netMassOut case in this block.
1129  ! At this point htot contains an Angstrom of fluid from layer 0 plus netMassIn.
1130  do ks=1,nz
1131  do i=is,ie ; if (ksort(i,ks) > 0) then
1132  k = ksort(i,ks)
1133 
1134  if ((htot(i) < angstrom) .and. (h(i,k) > eps(i,k))) then
1135  ! If less than an Angstrom was available from the layers above plus
1136  ! any precipitation, add more fluid from this layer.
1137  h_ent = min(angstrom-htot(i), h(i,k)-eps(i,k))
1138  htot(i) = htot(i) + h_ent
1139  h(i,k) = h(i,k) - h_ent
1140  d_eb(i,k) = d_eb(i,k) - h_ent
1141 
1142  r0_tot(i) = r0_tot(i) + h_ent*r0(i,k)
1143  uhtot(i) = uhtot(i) + h_ent*u(i,k)
1144  vhtot(i) = vhtot(i) + h_ent*v(i,k)
1145 
1146  rcv_tot(i) = rcv_tot(i) + h_ent*rcv(i,k)
1147  ttot(i) = ttot(i) + h_ent*t(i,k)
1148  stot(i) = stot(i) + h_ent*s(i,k)
1149  endif
1150 
1151  ! Water is removed from the topmost layers with any mass.
1152  ! We may lose layers if they are thin enough.
1153  ! The salt that is left behind goes into Stot.
1154  if ((massoutrem(i) > 0.0) .and. (h(i,k) > eps(i,k))) then
1155  if (massoutrem(i) > (h(i,k) - eps(i,k))) then
1156  h_evap = h(i,k) - eps(i,k)
1157  h(i,k) = eps(i,k)
1158  massoutrem(i) = massoutrem(i) - h_evap
1159  else
1160  h_evap = massoutrem(i)
1161  h(i,k) = h(i,k) - h_evap
1162  massoutrem(i) = 0.0
1163  endif
1164 
1165  stot(i) = stot(i) + h_evap*s(i,k)
1166  r0_tot(i) = r0_tot(i) + dr0_ds(i)*h_evap*s(i,k)
1167  rcv_tot(i) = rcv_tot(i) + drcv_ds(i)*h_evap*s(i,k)
1168  d_eb(i,k) = d_eb(i,k) - h_evap
1169 
1170  ! smg: when resolve the A=B code, we will set
1171  ! heat_content_massout = heat_content_massout - T(i,k)*h_evap*GV%H_to_RZ*fluxes%C_p*Idt
1172  ! by uncommenting the lines here.
1173  ! we will also then completely remove TempXpme from the model.
1174  if (associated(fluxes%heat_content_massout)) &
1175  fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) - &
1176  t(i,k)*h_evap*gv%H_to_RZ * fluxes%C_p * idt
1177  if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) - &
1178  t(i,k)*h_evap*gv%H_to_RZ
1179 
1180  endif
1181 
1182  ! The following section calculates how much fluid will be entrained.
1183  h_avail = h(i,k) - eps(i,k)
1184  if (h_avail > 0.0) then
1185  dr = r0_tot(i) - htot(i)*r0(i,k)
1186  h_ent = 0.0
1187 
1188  dr0 = dr
1189  do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1190  dr0 = dr0 - (dr0_dt(i)*pen_sw_bnd(n,i)) * &
1191  opacity_band(n,i,k)*htot(i)
1192  endif ; enddo
1193 
1194  ! Some entrainment will occur from this layer.
1195  if (dr0 > 0.0) then
1196  dr_comp = dr
1197  do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1198  ! Compare the density at the bottom of a layer with the
1199  ! density averaged over the mixed layer and that layer.
1200  opacity = opacity_band(n,i,k)
1201  sw_trans = exp(-h_avail*opacity)
1202  dr_comp = dr_comp + (dr0_dt(i)*pen_sw_bnd(n,i)) * &
1203  ((1.0 - sw_trans) - opacity*(htot(i)+h_avail)*sw_trans)
1204  endif ; enddo
1205  if (dr_comp >= 0.0) then
1206  ! The entire layer is entrained.
1207  h_ent = h_avail
1208  else
1209  ! The layer is partially entrained. Iterate to determine how much
1210  ! entrainment occurs. Solve for the h_ent at which dr_ent = 0.
1211 
1212  ! Instead of assuming that the curve is linear between the two end
1213  ! points, assume that the change is concentrated near small values
1214  ! of entrainment. On average, this saves about 1 iteration.
1215  frac = dr0 / (dr0 - dr_comp)
1216  h_ent = h_avail * frac*frac
1217  h_min = 0.0 ; h_max = h_avail
1218 
1219  do n=1,nsw
1220  r_sw_top(n) = dr0_dt(i) * pen_sw_bnd(n,i)
1221  c2(n) = r_sw_top(n) * opacity_band(n,i,k)**2
1222  enddo
1223  do itt=1,10
1224  dr_ent = dr ; dr_dh = 0.0
1225  do n=1,nsw
1226  opacity = opacity_band(n,i,k)
1227  sw_trans = exp(-h_ent*opacity)
1228  dr_ent = dr_ent + r_sw_top(n) * ((1.0 - sw_trans) - &
1229  opacity*(htot(i)+h_ent)*sw_trans)
1230  dr_dh = dr_dh + c2(n) * (htot(i)+h_ent) * sw_trans
1231  enddo
1232 
1233  if (dr_ent > 0.0) then
1234  h_min = h_ent
1235  else
1236  h_max = h_ent
1237  endif
1238 
1239  dh_newt = -dr_ent / dr_dh
1240  h_prev = h_ent ; h_ent = h_prev+dh_newt
1241  if (h_ent > h_max) then
1242  h_ent = 0.5*(h_prev+h_max)
1243  elseif (h_ent < h_min) then
1244  h_ent = 0.5*(h_prev+h_min)
1245  endif
1246 
1247  if (abs(dh_newt) < 0.2*angstrom) exit
1248  enddo
1249 
1250  endif
1251 
1252  ! Now that the amount of entrainment (h_ent) has been determined,
1253  ! calculate changes in various terms.
1254  sum_pen_en = 0.0 ; pen_absorbed = 0.0
1255  do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
1256  opacity = opacity_band(n,i,k)
1257  sw_trans = exp(-h_ent*opacity)
1258 
1259  x1 = h_ent*opacity
1260  if (x1 < 2.0e-5) then
1261  en_fn = (opacity*htot(i)*(1.0 - 0.5*(x1 - c1_3*x1)) + &
1262  x1*x1*c1_6)
1263  else
1264  en_fn = ((opacity*htot(i) + 2.0) * &
1265  ((1.0-sw_trans) / x1) - 1.0 + sw_trans)
1266  endif
1267  sum_pen_en = sum_pen_en - (dr0_dt(i)*pen_sw_bnd(n,i)) * en_fn
1268 
1269  pen_absorbed = pen_absorbed + pen_sw_bnd(n,i) * (1.0 - sw_trans)
1270  pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
1271  endif ; enddo
1272 
1273  conv_en(i) = conv_en(i) + g_h2_2rho0 * h_ent * &
1274  ( (r0_tot(i) - r0(i,k)*htot(i)) + sum_pen_en )
1275 
1276  r0_tot(i) = r0_tot(i) + (h_ent * r0(i,k) + pen_absorbed*dr0_dt(i))
1277  stot(i) = stot(i) + h_ent * s(i,k)
1278  ttot(i) = ttot(i) + (h_ent * t(i,k) + pen_absorbed)
1279  rcv_tot(i) = rcv_tot(i) + (h_ent * rcv(i,k) + pen_absorbed*drcv_dt(i))
1280  endif ! dr0 > 0.0
1281 
1282  if (h_ent > 0.0) then
1283  if (htot(i) > 0.0) &
1284  dke_fc(i) = dke_fc(i) + cs%bulk_Ri_convective * 0.5 * &
1285  ((gv%H_to_Z*h_ent) / (htot(i)*(h_ent+htot(i)))) * &
1286  ((uhtot(i)-u(i,k)*htot(i))**2 + (vhtot(i)-v(i,k)*htot(i))**2)
1287 
1288  htot(i) = htot(i) + h_ent
1289  h(i,k) = h(i,k) - h_ent
1290  d_eb(i,k) = d_eb(i,k) - h_ent
1291  uhtot(i) = u(i,k)*h_ent ; vhtot(i) = v(i,k)*h_ent
1292  !### I think that the line above should instead be:
1293  ! uhtot(i) = uhtot(i) + h_ent*u(i,k) ; vhtot(i) = vhtot(i) + h_ent*v(i,k)
1294  endif
1295 
1296 
1297  endif ! h_avail>0
1298  endif ; enddo ! i loop
1299  enddo ! k loop
1300 

Referenced by bulkmixedlayer().

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(unit_scale_type), intent(in)  US,
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 [H ~> m or kg m-2]. Layer 0 is the new mixed layer.
[in,out]tPotential temperature [degC].
[in,out]sSalinity [ppt].
[in,out]r0Potential density referenced to surface pressure [R ~> kg m-3].
[in,out]rcvThe coordinate defining potential density [R ~> kg m-3].
[in]rcvtgtThe target value of Rcv for each layer [R ~> kg m-3].
[in]dtTime increment [T ~> s].
[in]dt_diagThe accumulated time interval for diagnostics [T ~> s].
[in,out]d_eaThe upward increase across a layer in the entrainment from above [H ~> m or kg m-2]. Positive d_ea goes with layer thickness increases.
[in,out]d_ebThe downward increase across a layer in the entrainment from below [H ~> m or kg m-2]. Positive values go with mass gain by a layer.
[in]jThe meridional row to work on.
[in]usA dimensional unit scaling type
csThe control structure returned by a previous call to mixedlayer_init.
[in]drcv_dtThe partial derivative of coordinate defining potential density with potential temperature [R degC-1 ~> kg m-3 degC-1].
[in]drcv_dsThe partial derivative of coordinate defining potential density with salinity [R ppt-1 ~> kg m-3 ppt-1].
[in]max_bl_detIf non-negative, the maximum detrainment permitted from the buffer layers [H ~> m or kg m-2].

Definition at line 3097 of file MOM_bulk_mixed_layer.F90.

3097  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
3098  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
3099  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2].
3100  !! Layer 0 is the new mixed layer.
3101  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: T !< Potential temperature [degC].
3102  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: S !< Salinity [ppt].
3103  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: R0 !< Potential density referenced to
3104  !! surface pressure [R ~> kg m-3].
3105  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: Rcv !< The coordinate defining potential
3106  !! density [R ~> kg m-3].
3107  real, dimension(SZK_(GV)), intent(in) :: RcvTgt !< The target value of Rcv for each
3108  !! layer [R ~> kg m-3].
3109  real, intent(in) :: dt !< Time increment [T ~> s].
3110  real, intent(in) :: dt_diag !< The accumulated time interval for
3111  !! diagnostics [T ~> s].
3112  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_ea !< The upward increase across a layer in
3113  !! the entrainment from above
3114  !! [H ~> m or kg m-2]. Positive d_ea
3115  !! goes with layer thickness increases.
3116  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_eb !< The downward increase across a layer
3117  !! in the entrainment from below [H ~> m or kg m-2].
3118  !! Positive values go with mass gain by
3119  !! a layer.
3120  integer, intent(in) :: j !< The meridional row to work on.
3121  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3122  type(bulkmixedlayer_CS), pointer :: CS !< The control structure returned by a
3123  !! previous call to mixedlayer_init.
3124  real, dimension(SZI_(G)), intent(in) :: dRcv_dT !< The partial derivative of
3125  !! coordinate defining potential density
3126  !! with potential temperature
3127  !! [R degC-1 ~> kg m-3 degC-1].
3128  real, dimension(SZI_(G)), intent(in) :: dRcv_dS !< The partial derivative of
3129  !! coordinate defining potential density
3130  !! with salinity [R ppt-1 ~> kg m-3 ppt-1].
3131  real, dimension(SZI_(G)), intent(in) :: max_BL_det !< If non-negative, the maximum
3132  !! detrainment permitted from the buffer
3133  !! layers [H ~> m or kg m-2].
3134 
3135  ! Local variables
3136  real :: Ih ! The inverse of a thickness [H-1 ~> m-1 or m2 kg-1].
3137  real :: h_ent ! The thickness from a layer that is
3138  ! entrained [H ~> m or kg m-2].
3139  real :: max_det_rem(SZI_(G)) ! Remaining permitted detrainment [H ~> m or kg m-2].
3140  real :: detrain(SZI_(G)) ! The thickness of fluid to detrain
3141  ! from the mixed layer [H ~> m or kg m-2].
3142  real :: dT_dR, dS_dR, dRml, dR0_dRcv, dT_dS_wt2
3143  real :: I_denom ! A work variable [ppt2 R-2 ~> ppt2 m6 kg-2].
3144  real :: Sdown, Tdown
3145  real :: dt_Time ! The timestep divided by the detrainment timescale [nondim].
3146  real :: g_H2_2Rho0dt ! Half the gravitational acceleration times the square of the
3147  ! conversion from H to m divided by the mean density times the time
3148  ! step [L2 Z T-3 H-2 R-1 ~> m4 s-3 kg-1 or m10 s-3 kg-3].
3149  real :: g_H2_2dt ! Half the gravitational acceleration times the square of the
3150  ! conversion from H to Z divided by the diagnostic time step
3151  ! [L2 Z H-2 T-3 ~> m s-3 or m7 kg-2 s-3].
3152 
3153  logical :: splittable_BL(SZI_(G)), orthogonal_extrap
3154  real :: x1
3155 
3156  integer :: i, is, ie, k, k1, nkmb, nz
3157  is = g%isc ; ie = g%iec ; nz = gv%ke
3158  nkmb = cs%nkml+cs%nkbl
3159  if (cs%nkbl /= 1) call mom_error(fatal,"MOM_mixed_layer: "// &
3160  "CS%nkbl must be 1 in mixedlayer_detrain_1.")
3161 
3162  dt_time = dt / cs%BL_detrain_time
3163  g_h2_2rho0dt = (gv%g_Earth * gv%H_to_Z**2) / (2.0 * gv%Rho0 * dt_diag)
3164  g_h2_2dt = (gv%g_Earth * gv%H_to_Z**2) / (2.0 * dt_diag)
3165 
3166  ! Move detrained water into the buffer layer.
3167  do k=1,cs%nkml
3168  do i=is,ie ; if (h(i,k) > 0.0) then
3169  ih = 1.0 / (h(i,nkmb) + h(i,k))
3170  if (cs%TKE_diagnostics) &
3171  cs%diag_TKE_conv_s2(i,j) = cs%diag_TKE_conv_s2(i,j) + &
3172  g_h2_2rho0dt * h(i,k) * h(i,nkmb) * (r0(i,nkmb) - r0(i,k))
3173  if (allocated(cs%diag_PE_detrain)) &
3174  cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + &
3175  g_h2_2dt * h(i,k) * h(i,nkmb) * (r0(i,nkmb) - r0(i,k))
3176  if (allocated(cs%diag_PE_detrain2)) &
3177  cs%diag_PE_detrain2(i,j) = cs%diag_PE_detrain2(i,j) + &
3178  g_h2_2dt * h(i,k) * h(i,nkmb) * (r0(i,nkmb) - r0(i,k))
3179 
3180  r0(i,nkmb) = (r0(i,nkmb)*h(i,nkmb) + r0(i,k)*h(i,k)) * ih
3181  rcv(i,nkmb) = (rcv(i,nkmb)*h(i,nkmb) + rcv(i,k)*h(i,k)) * ih
3182  t(i,nkmb) = (t(i,nkmb)*h(i,nkmb) + t(i,k)*h(i,k)) * ih
3183  s(i,nkmb) = (s(i,nkmb)*h(i,nkmb) + s(i,k)*h(i,k)) * ih
3184 
3185  d_ea(i,k) = d_ea(i,k) - h(i,k)
3186  d_ea(i,nkmb) = d_ea(i,nkmb) + h(i,k)
3187  h(i,nkmb) = h(i,nkmb) + h(i,k)
3188  h(i,k) = 0.0
3189  endif ; enddo
3190  enddo
3191 
3192  do i=is,ie
3193  max_det_rem(i) = 10.0 * h(i,nkmb)
3194  if (max_bl_det(i) >= 0.0) max_det_rem(i) = max_bl_det(i)
3195  enddo
3196 
3197 ! If the mixed layer was denser than the densest interior layer,
3198 ! but is now lighter than this layer, leaving a buffer layer that
3199 ! is denser than this layer, there are problems. This should prob-
3200 ! ably be considered a case of an inadequate choice of resolution in
3201 ! density space and should be avoided. To make the model run sens-
3202 ! ibly in this case, it will make the mixed layer denser while making
3203 ! the buffer layer the density of the densest interior layer (pro-
3204 ! vided that the this will not make the mixed layer denser than the
3205 ! interior layer). Otherwise, make the mixed layer the same density
3206 ! as the densest interior layer and lighten the buffer layer with
3207 ! the released buoyancy. With multiple buffer layers, much more
3208 ! graceful options are available.
3209  do i=is,ie ; if (h(i,nkmb) > 0.0) then
3210  if ((r0(i,0)<r0(i,nz)) .and. (r0(i,nz)<r0(i,nkmb))) then
3211  if ((r0(i,nz)-r0(i,0))*h(i,0) > &
3212  (r0(i,nkmb)-r0(i,nz))*h(i,nkmb)) then
3213  detrain(i) = (r0(i,nkmb)-r0(i,nz))*h(i,nkmb) / (r0(i,nkmb)-r0(i,0))
3214  else
3215  detrain(i) = (r0(i,nz)-r0(i,0))*h(i,0) / (r0(i,nkmb)-r0(i,0))
3216  endif
3217 
3218  d_eb(i,cs%nkml) = d_eb(i,cs%nkml) + detrain(i)
3219  d_ea(i,cs%nkml) = d_ea(i,cs%nkml) - detrain(i)
3220  d_eb(i,nkmb) = d_eb(i,nkmb) - detrain(i)
3221  d_ea(i,nkmb) = d_ea(i,nkmb) + detrain(i)
3222 
3223  if (allocated(cs%diag_PE_detrain)) cs%diag_PE_detrain(i,j) = &
3224  cs%diag_PE_detrain(i,j) + g_h2_2dt * detrain(i)* &
3225  (h(i,0) + h(i,nkmb)) * (r0(i,nkmb) - r0(i,0))
3226  x1 = r0(i,0)
3227  r0(i,0) = r0(i,0) - detrain(i)*(r0(i,0)-r0(i,nkmb)) / h(i,0)
3228  r0(i,nkmb) = r0(i,nkmb) - detrain(i)*(r0(i,nkmb)-x1) / h(i,nkmb)
3229  x1 = rcv(i,0)
3230  rcv(i,0) = rcv(i,0) - detrain(i)*(rcv(i,0)-rcv(i,nkmb)) / h(i,0)
3231  rcv(i,nkmb) = rcv(i,nkmb) - detrain(i)*(rcv(i,nkmb)-x1) / h(i,nkmb)
3232  x1 = t(i,0)
3233  t(i,0) = t(i,0) - detrain(i)*(t(i,0)-t(i,nkmb)) / h(i,0)
3234  t(i,nkmb) = t(i,nkmb) - detrain(i)*(t(i,nkmb)-x1) / h(i,nkmb)
3235  x1 = s(i,0)
3236  s(i,0) = s(i,0) - detrain(i)*(s(i,0)-s(i,nkmb)) / h(i,0)
3237  s(i,nkmb) = s(i,nkmb) - detrain(i)*(s(i,nkmb)-x1) / h(i,nkmb)
3238 
3239  endif
3240  endif ; enddo
3241 
3242  ! Move water out of the buffer layer, if convenient.
3243 ! Split the buffer layer if possible, and replace the buffer layer
3244 ! with a small amount of fluid from the mixed layer.
3245 ! This is the exponential-in-time splitting, circa 2005.
3246  do i=is,ie
3247  if (h(i,nkmb) > 0.0) then ; splittable_bl(i) = .true.
3248  else ; splittable_bl(i) = .false. ; endif
3249  enddo
3250 
3251  dt_ds_wt2 = cs%dT_dS_wt**2
3252 
3253  do k=nz-1,nkmb+1,-1 ; do i=is,ie
3254  if (splittable_bl(i)) then
3255  if (rcvtgt(k)<=rcv(i,nkmb)) then
3256 ! Estimate dR/drho, dTheta/dR, and dS/dR, where R is the coordinate variable
3257 ! and rho is in-situ (or surface) potential density.
3258 ! There is no "right" way to do this, so this keeps things reasonable, if
3259 ! slightly arbitrary.
3260  splittable_bl(i) = .false.
3261 
3262  k1 = k+1 ; orthogonal_extrap = .false.
3263  ! Here we try to find a massive layer to use for interpolating the
3264  ! temperature and salinity. If none is available a pseudo-orthogonal
3265  ! extrapolation is used. The 10.0 and 0.9 in the following are
3266  ! arbitrary but probably about right.
3267  if ((h(i,k+1) < 10.0*gv%Angstrom_H) .or. &
3268  ((rcvtgt(k+1)-rcv(i,nkmb)) >= 0.9*(rcv(i,k1) - rcv(i,0)))) then
3269  if (k>=nz-1) then ; orthogonal_extrap = .true.
3270  elseif ((h(i,k+2) <= 10.0*gv%Angstrom_H) .and. &
3271  ((rcvtgt(k+1)-rcv(i,nkmb)) < 0.9*(rcv(i,k+2)-rcv(i,0)))) then
3272  k1 = k+2
3273  else ; orthogonal_extrap = .true. ; endif
3274  endif
3275 
3276  if ((r0(i,0) >= r0(i,k1)) .or. (rcv(i,0) >= rcv(i,nkmb))) cycle
3277  ! In this case there is an inversion of in-situ density relative to
3278  ! the coordinate variable. Do not detrain from the buffer layer.
3279 
3280  if (orthogonal_extrap) then
3281  ! 36 here is a typical oceanic value of (dR/dS) / (dR/dT) - it says
3282  ! that the relative weights of T & S changes is a plausible 6:1.
3283  ! Also, this was coded on Athena's 6th birthday!
3284  i_denom = 1.0 / (drcv_ds(i)**2 + dt_ds_wt2*drcv_dt(i)**2)
3285  dt_dr = dt_ds_wt2*drcv_dt(i) * i_denom
3286  ds_dr = drcv_ds(i) * i_denom
3287  else
3288  dt_dr = (t(i,0) - t(i,k1)) / (rcv(i,0) - rcv(i,k1))
3289  ds_dr = (s(i,0) - s(i,k1)) / (rcv(i,0) - rcv(i,k1))
3290  endif
3291  drml = dt_time * (r0(i,nkmb) - r0(i,0)) * &
3292  (rcv(i,0) - rcv(i,k1)) / (r0(i,0) - r0(i,k1))
3293  ! Once again, there is an apparent density inversion in Rcv.
3294  if (drml < 0.0) cycle
3295  dr0_drcv = (r0(i,0) - r0(i,k1)) / (rcv(i,0) - rcv(i,k1))
3296 
3297  if ((rcv(i,nkmb) - drml < rcvtgt(k)) .and. (max_det_rem(i) > h(i,nkmb))) then
3298  ! In this case, the buffer layer is split into two isopycnal layers.
3299  detrain(i) = h(i,nkmb)*(rcv(i,nkmb) - rcvtgt(k)) / &
3300  (rcvtgt(k+1) - rcvtgt(k))
3301 
3302  if (allocated(cs%diag_PE_detrain)) cs%diag_PE_detrain(i,j) = &
3303  cs%diag_PE_detrain(i,j) - g_h2_2dt * detrain(i) * &
3304  (h(i,nkmb)-detrain(i)) * (rcvtgt(k+1) - rcvtgt(k)) * dr0_drcv
3305 
3306  tdown = detrain(i) * (t(i,nkmb) + dt_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3307  t(i,k) = (h(i,k) * t(i,k) + &
3308  (h(i,nkmb) * t(i,nkmb) - tdown)) / &
3309  (h(i,k) + (h(i,nkmb) - detrain(i)))
3310  t(i,k+1) = (h(i,k+1) * t(i,k+1) + tdown)/ &
3311  (h(i,k+1) + detrain(i))
3312  t(i,nkmb) = t(i,0)
3313  sdown = detrain(i) * (s(i,nkmb) + ds_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3314  s(i,k) = (h(i,k) * s(i,k) + &
3315  (h(i,nkmb) * s(i,nkmb) - sdown)) / &
3316  (h(i,k) + (h(i,nkmb) - detrain(i)))
3317  s(i,k+1) = (h(i,k+1) * s(i,k+1) + sdown)/ &
3318  (h(i,k+1) + detrain(i))
3319  s(i,nkmb) = s(i,0)
3320  rcv(i,nkmb) = rcv(i,0)
3321 
3322  d_ea(i,k+1) = d_ea(i,k+1) + detrain(i)
3323  d_ea(i,k) = d_ea(i,k) + (h(i,nkmb) - detrain(i))
3324  d_ea(i,nkmb) = d_ea(i,nkmb) - h(i,nkmb)
3325 
3326  h(i,k+1) = h(i,k+1) + detrain(i)
3327  h(i,k) = h(i,k) + h(i,nkmb) - detrain(i)
3328  h(i,nkmb) = 0.0
3329  else
3330  ! Here only part of the buffer layer is moved into the interior.
3331  detrain(i) = h(i,nkmb) * drml / (rcvtgt(k+1) - rcv(i,nkmb) + drml)
3332  if (detrain(i) > max_det_rem(i)) detrain(i) = max_det_rem(i)
3333  ih = 1.0 / (h(i,k+1) + detrain(i))
3334 
3335  tdown = (t(i,nkmb) + dt_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3336  t(i,nkmb) = t(i,nkmb) - dt_dr * drml
3337  t(i,k+1) = (h(i,k+1) * t(i,k+1) + detrain(i) * tdown) * ih
3338  sdown = (s(i,nkmb) + ds_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
3339 ! The following two expressions updating S(nkmb) are mathematically identical.
3340 ! S(i,nkmb) = (h(i,nkmb) * S(i,nkmb) - detrain(i) * Sdown) / &
3341 ! (h(i,nkmb) - detrain(i))
3342  s(i,nkmb) = s(i,nkmb) - ds_dr * drml
3343  s(i,k+1) = (h(i,k+1) * s(i,k+1) + detrain(i) * sdown) * ih
3344 
3345  d_ea(i,k+1) = d_ea(i,k+1) + detrain(i)
3346  d_ea(i,nkmb) = d_ea(i,nkmb) - detrain(i)
3347 
3348  h(i,k+1) = h(i,k+1) + detrain(i)
3349  h(i,nkmb) = h(i,nkmb) - detrain(i)
3350 
3351  if (allocated(cs%diag_PE_detrain)) cs%diag_PE_detrain(i,j) = &
3352  cs%diag_PE_detrain(i,j) - g_h2_2dt * detrain(i) * dr0_drcv * &
3353  (h(i,nkmb)-detrain(i)) * (rcvtgt(k+1) - rcv(i,nkmb) + drml)
3354  endif
3355  endif ! RcvTgt(k)<=Rcv(i,nkmb)
3356  endif ! splittable_BL
3357  enddo ; enddo ! i & k loops
3358 
3359 ! The numerical behavior of the buffer layer is dramatically improved
3360 ! if it is always at least a small fraction (say 10%) of the thickness
3361 ! of the mixed layer. As the physical distinction between the mixed
3362 ! and buffer layers is vague anyway, this seems hard to argue against.
3363  do i=is,ie
3364  if (h(i,nkmb) < 0.1*h(i,0)) then
3365  h_ent = 0.1*h(i,0) - h(i,nkmb)
3366  ih = 10.0/h(i,0)
3367  t(i,nkmb) = (h(i,nkmb)*t(i,nkmb) + h_ent*t(i,0)) * ih
3368  s(i,nkmb) = (h(i,nkmb)*s(i,nkmb) + h_ent*s(i,0)) * ih
3369 
3370  d_ea(i,1) = d_ea(i,1) - h_ent
3371  d_ea(i,nkmb) = d_ea(i,nkmb) + h_ent
3372 
3373  h(i,0) = h(i,0) - h_ent
3374  h(i,nkmb) = h(i,nkmb) + h_ent
3375  endif
3376  enddo
3377 

References mom_error_handler::mom_error().

Referenced by bulkmixedlayer().

Here is the call graph for this function:
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(unit_scale_type), intent(in)  US,
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 [H ~> m or kg m-2]. Layer 0 is the new mixed layer.
[in,out]tPotential temperature [degC].
[in,out]sSalinity [ppt].
[in,out]r0Potential density referenced to surface pressure [R ~> kg m-3].
[in,out]rcvThe coordinate defining potential density [R ~> kg m-3].
[in]rcvtgtThe target value of Rcv for each layer [R ~> kg m-3].
[in]dtTime increment [T ~> s].
[in]dt_diagThe diagnostic time step [T ~> s].
[in,out]d_eaThe upward increase across a layer in the entrainment from above [H ~> m or kg m-2]. Positive d_ea goes with layer thickness increases.
[in]jThe meridional row to work on.
[in]usA dimensional unit scaling type
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, [R degC-1 ~> kg m-3 degC-1].
[in]dr0_dsThe partial derivative of cpotential density referenced to the surface with salinity [R ppt-1 ~> kg m-3 ppt-1].
[in]drcv_dtThe partial derivative of coordinate defining potential density with potential temperature, [R degC-1 ~> kg m-3 degC-1].
[in]drcv_dsThe partial derivative of coordinate defining potential density with salinity [R ppt-1 ~> kg m-3 ppt-1].
[in]max_bl_detIf non-negative, the maximum detrainment permitted from the buffer layers [H ~> m or kg m-2].

Definition at line 2206 of file MOM_bulk_mixed_layer.F90.

2206  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2207  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2208  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2].
2209  !! Layer 0 is the new mixed layer.
2210  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: T !< Potential temperature [degC].
2211  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: S !< Salinity [ppt].
2212  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: R0 !< Potential density referenced to
2213  !! surface pressure [R ~> kg m-3].
2214  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: Rcv !< The coordinate defining potential
2215  !! density [R ~> kg m-3].
2216  real, dimension(SZK_(GV)), intent(in) :: RcvTgt !< The target value of Rcv for each
2217  !! layer [R ~> kg m-3].
2218  real, intent(in) :: dt !< Time increment [T ~> s].
2219  real, intent(in) :: dt_diag !< The diagnostic time step [T ~> s].
2220  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_ea !< The upward increase across a layer in
2221  !! the entrainment from above
2222  !! [H ~> m or kg m-2]. Positive d_ea
2223  !! goes with layer thickness increases.
2224  integer, intent(in) :: j !< The meridional row to work on.
2225  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2226  type(bulkmixedlayer_CS), pointer :: CS !< The control structure returned by a
2227  !! previous call to mixedlayer_init.
2228  real, dimension(SZI_(G)), intent(in) :: dR0_dT !< The partial derivative of
2229  !! potential density referenced to the
2230  !! surface with potential temperature,
2231  !! [R degC-1 ~> kg m-3 degC-1].
2232  real, dimension(SZI_(G)), intent(in) :: dR0_dS !< The partial derivative of
2233  !! cpotential density referenced to the
2234  !! surface with salinity
2235  !! [R ppt-1 ~> kg m-3 ppt-1].
2236  real, dimension(SZI_(G)), intent(in) :: dRcv_dT !< The partial derivative of
2237  !! coordinate defining potential density
2238  !! with potential temperature,
2239  !! [R degC-1 ~> kg m-3 degC-1].
2240  real, dimension(SZI_(G)), intent(in) :: dRcv_dS !< The partial derivative of
2241  !! coordinate defining potential density
2242  !! with salinity [R ppt-1 ~> kg m-3 ppt-1].
2243  real, dimension(SZI_(G)), intent(in) :: max_BL_det !< If non-negative, the maximum
2244  !! detrainment permitted from the buffer
2245  !! layers [H ~> m or kg m-2].
2246 
2247 ! This subroutine moves any water left in the former mixed layers into the
2248 ! two buffer layers and may also move buffer layer water into the interior
2249 ! isopycnal layers.
2250 
2251  ! Local variables
2252  real :: h_to_bl ! The total thickness detrained to the buffer
2253  ! layers [H ~> m or kg m-2].
2254  real :: R0_to_bl ! The depth integrated amount of R0 that is detrained to the
2255  ! buffer layer [H R ~> kg m-2 or kg2 m-5]
2256  real :: Rcv_to_bl ! The depth integrated amount of Rcv that is detrained to the
2257  ! buffer layer [H R ~> kg m-2 or kg2 m-5]
2258  real :: T_to_bl ! The depth integrated amount of T that is detrained to the
2259  ! buffer layer [degC H ~> degC m or degC kg m-2]
2260  real :: S_to_bl ! The depth integrated amount of S that is detrained to the
2261  ! buffer layer [ppt H ~> ppt m or ppt kg m-2]
2262  real :: h_min_bl ! The minimum buffer layer thickness [H ~> m or kg m-2].
2263 
2264  real :: h1, h2 ! Scalar variables holding the values of
2265  ! h(i,CS%nkml+1) and h(i,CS%nkml+2) [H ~> m or kg m-2].
2266  real :: h1_avail ! The thickness of the upper buffer layer
2267  ! available to move into the lower buffer
2268  ! layer [H ~> m or kg m-2].
2269  real :: stays ! stays is the thickness of the upper buffer
2270  ! layer that remains there [H ~> m or kg m-2].
2271  real :: stays_min, stays_max ! The minimum and maximum permitted values of
2272  ! stays [H ~> m or kg m-2].
2273 
2274  logical :: mergeable_bl ! If true, it is an option to combine the two
2275  ! buffer layers and create water that matches
2276  ! the target density of an interior layer.
2277  real :: stays_merge ! If the two buffer layers can be combined
2278  ! stays_merge is the thickness of the upper
2279  ! layer that remains [H ~> m or kg m-2].
2280  real :: stays_min_merge ! The minimum allowed value of stays_merge [H ~> m or kg m-2].
2281 
2282  real :: dR0_2dz, dRcv_2dz ! Half the vertical gradients of R0 and Rcv [R H-1 ~> kg m-4 or m-1]
2283 ! real :: dT_2dz, dS_2dz ! Half the vertical gradients of T and S, in degC H-1, and ppt H-1.
2284  real :: scale_slope ! A nondimensional number < 1 used to scale down
2285  ! the slope within the upper buffer layer when
2286  ! water MUST be detrained to the lower layer.
2287 
2288  real :: dPE_extrap ! The potential energy change due to dispersive
2289  ! advection or mixing layers, divided by
2290  ! rho_0*g [H2 ~> m2 or kg2 m-4].
2291  real :: dPE_det, dPE_merge ! The energy required to mix the detrained water
2292  ! into the buffer layer or the merge the two
2293  ! buffer layers [R H2 L2 Z-1 T-2 ~> J m-2 or J kg2 m-8].
2294 
2295  real :: h_from_ml ! The amount of additional water that must be
2296  ! drawn from the mixed layer [H ~> m or kg m-2].
2297  real :: h_det_h2 ! The amount of detrained water and mixed layer
2298  ! water that will go directly into the lower
2299  ! buffer layer [H ~> m or kg m-2].
2300  real :: h_det_to_h2, h_ml_to_h2 ! All of the variables hA_to_hB are the thickness fluxes
2301  real :: h_det_to_h1, h_ml_to_h1 ! from one layer to another [H ~> m or kg m-2],
2302  real :: h1_to_h2, h1_to_k0 ! with h_det the detrained water, h_ml
2303  real :: h2_to_k1, h2_to_k1_rem ! the actively mixed layer, h1 and h2 the upper
2304  ! and lower buffer layers, and k0 and k1 the
2305  ! interior layers that are just lighter and
2306  ! just denser than the lower buffer layer.
2307 
2308  real :: R0_det, T_det, S_det ! Detrained values of R0 [R ~> kg m-3], T [degC], and S [ppt].
2309  real :: Rcv_stays, R0_stays ! Values of Rcv and R0 that stay in a layer.
2310  real :: T_stays, S_stays ! Values of T and S that stay in a layer.
2311  real :: dSpice_det, dSpice_stays! The spiciness difference between an original
2312  ! buffer layer and the water that moves into
2313  ! an interior layer or that stays in that
2314  ! layer [R ~> kg m-3].
2315  real :: dSpice_lim, dSpice_lim2 ! Limits to the spiciness difference between
2316  ! the lower buffer layer and the water that
2317  ! moves into an interior layer [R ~> kg m-3].
2318  real :: dSpice_2dz ! The vertical gradient of spiciness used for
2319  ! advection [R H-1 ~> kg m-4 or m-1].
2320 
2321  real :: dPE_ratio ! Multiplier of dPE_det at which merging is
2322  ! permitted - here (detrainment_per_day/dt)*30
2323  ! days?
2324  real :: num_events ! The number of detrainment events over which
2325  ! to prefer merging the buffer layers.
2326  real :: dPE_time_ratio ! Larger of 1 and the detrainment timescale over dt [nondim].
2327  real :: dT_dS_gauge, dS_dT_gauge ! The relative scales of temperature and
2328  ! salinity changes in defining spiciness, in
2329  ! [degC ppt-1] and [ppt degC-1].
2330  real :: I_denom ! A work variable with units of [ppt2 R-2 ~> ppt2 m6 kg-2].
2331 
2332  real :: g_2 ! 1/2 g_Earth [L2 Z-1 T-2 ~> m s-2].
2333  real :: Rho0xG ! Rho0 times G_Earth [R L2 Z-1 T-2 ~> kg m-2 s-2].
2334  real :: I2Rho0 ! 1 / (2 Rho0) [R-1 ~> m3 kg-1].
2335  real :: Idt_H2 ! The square of the conversion from thickness to Z
2336  ! divided by the time step [Z2 H-2 T-1 ~> s-1 or m6 kg-2 s-1].
2337  logical :: stable_Rcv ! If true, the buffer layers are stable with
2338  ! respect to the coordinate potential density.
2339  real :: h_neglect ! A thickness that is so small it is usually lost
2340  ! in roundoff and can be neglected [H ~> m or kg m-2].
2341 
2342  real :: s1en ! A work variable [H2 L2 kg m-1 T-3 ~> kg m3 s-3 or kg3 m-3 s-3].
2343  real :: s1, s2, bh0 ! Work variables [H ~> m or kg m-2].
2344  real :: s3sq ! A work variable [H2 ~> m2 or kg2 m-4].
2345  real :: I_ya, b1 ! Nondimensional work variables.
2346  real :: Ih, Ihdet, Ih1f, Ih2f ! Assorted inverse thickness work variables,
2347  real :: Ihk0, Ihk1, Ih12 ! all in [H-1 ~> m-1 or m2 kg-1].
2348  real :: dR1, dR2, dR2b, dRk1 ! Assorted density difference work variables,
2349  real :: dR0, dR21, dRcv ! all in [R ~> kg m-3].
2350  real :: dRcv_stays, dRcv_det, dRcv_lim
2351  real :: Angstrom ! The minumum layer thickness [H ~> m or kg m-2].
2352 
2353  real :: h2_to_k1_lim, T_new, S_new, T_max, T_min, S_max, S_min
2354  character(len=200) :: mesg
2355 
2356  integer :: i, k, k0, k1, is, ie, nz, kb1, kb2, nkmb
2357  is = g%isc ; ie = g%iec ; nz = gv%ke
2358  kb1 = cs%nkml+1; kb2 = cs%nkml+2
2359  nkmb = cs%nkml+cs%nkbl
2360  h_neglect = gv%H_subroundoff
2361  g_2 = 0.5 * gv%g_Earth
2362  rho0xg = gv%Rho0 * gv%g_Earth
2363  idt_h2 = gv%H_to_Z**2 / dt_diag
2364  i2rho0 = 0.5 / (gv%Rho0)
2365  angstrom = gv%Angstrom_H
2366 
2367  ! This is hard coding of arbitrary and dimensional numbers.
2368  dt_ds_gauge = cs%dT_dS_wt ; ds_dt_gauge = 1.0 / dt_ds_gauge
2369  num_events = 10.0
2370 
2371  if (cs%nkbl /= 2) call mom_error(fatal, "MOM_mixed_layer"// &
2372  "CS%nkbl must be 2 in mixedlayer_detrain_2.")
2373 
2374  if (dt < cs%BL_detrain_time) then ; dpe_time_ratio = cs%BL_detrain_time / (dt)
2375  else ; dpe_time_ratio = 1.0 ; endif
2376 
2377  do i=is,ie
2378 
2379  ! Determine all of the properties being detrained from the mixed layer.
2380 
2381  ! As coded this has the k and i loop orders switched, but k is CS%nkml is
2382  ! often just 1 or 2, so this seems like it should not be a problem, especially
2383  ! since it means that a number of variables can now be scalars, not arrays.
2384  h_to_bl = 0.0 ; r0_to_bl = 0.0
2385  rcv_to_bl = 0.0 ; t_to_bl = 0.0 ; s_to_bl = 0.0
2386 
2387  do k=1,cs%nkml ; if (h(i,k) > 0.0) then
2388  h_to_bl = h_to_bl + h(i,k)
2389  r0_to_bl = r0_to_bl + r0(i,k)*h(i,k)
2390 
2391  rcv_to_bl = rcv_to_bl + rcv(i,k)*h(i,k)
2392  t_to_bl = t_to_bl + t(i,k)*h(i,k)
2393  s_to_bl = s_to_bl + s(i,k)*h(i,k)
2394 
2395  d_ea(i,k) = d_ea(i,k) - h(i,k)
2396  h(i,k) = 0.0
2397  endif ; enddo
2398  if (h_to_bl > 0.0) then ; r0_det = r0_to_bl / h_to_bl
2399  else ; r0_det = r0(i,0) ; endif
2400 
2401  ! This code does both downward detrainment from both the mixed layer and the
2402  ! buffer layers.
2403  ! Several considerations apply in detraining water into the interior:
2404  ! (1) Water only moves into the interior from the deeper buffer layer,
2405  ! so the deeper buffer layer must have some mass.
2406  ! (2) The upper buffer layer must have some mass so the extrapolation of
2407  ! density is meaningful (i.e. there is not detrainment from the buffer
2408  ! layers when there is strong mixed layer entrainment).
2409  ! (3) The lower buffer layer density extrapolated to its base with a
2410  ! linear fit between the two layers must exceed the density of the
2411  ! next denser interior layer.
2412  ! (4) The average extroplated coordinate density that is moved into the
2413  ! isopycnal interior matches the target value for that layer.
2414  ! (5) The potential energy change is calculated and might be used later
2415  ! to allow the upper buffer layer to mix more into the lower buffer
2416  ! layer.
2417 
2418  ! Determine whether more must be detrained from the mixed layer to keep a
2419  ! minimal amount of mass in the buffer layers. In this case the 5% of the
2420  ! mixed layer thickness is hard-coded, but probably shouldn't be!
2421  h_min_bl = min(cs%Hbuffer_min, cs%Hbuffer_rel_min*h(i,0))
2422 
2423  stable_rcv = .true.
2424  if (((r0(i,kb2)-r0(i,kb1)) * (rcv(i,kb2)-rcv(i,kb1)) <= 0.0)) &
2425  stable_rcv = .false.
2426 
2427  h1 = h(i,kb1) ; h2 = h(i,kb2)
2428 
2429  h2_to_k1_rem = (h1 + h2) + h_to_bl
2430  if ((max_bl_det(i) >= 0.0) .and. (h2_to_k1_rem > max_bl_det(i))) &
2431  h2_to_k1_rem = max_bl_det(i)
2432 
2433 
2434  if ((h2 == 0.0) .and. (h1 > 0.0)) then
2435  ! The lower buffer layer has been eliminated either by convective
2436  ! adjustment or entrainment from the interior, and its current properties
2437  ! are not meaningful, but may later be used to determine the properties of
2438  ! waters moving into the lower buffer layer. So the properties of the
2439  ! lower buffer layer are set to be between those of the upper buffer layer
2440  ! and the next denser interior layer, measured by R0. This probably does
2441  ! not happen very often, so I am not too worried about the inefficiency of
2442  ! the following loop.
2443  do k1=kb2+1,nz ; if (h(i,k1) > 2.0*angstrom) exit ; enddo
2444 
2445  r0(i,kb2) = r0(i,kb1)
2446 
2447  rcv(i,kb2)=rcv(i,kb1) ; t(i,kb2)=t(i,kb1) ; s(i,kb2)=s(i,kb1)
2448 
2449 
2450  if (k1 <= nz) then ; if (r0(i,k1) >= r0(i,kb1)) then
2451  r0(i,kb2) = 0.5*(r0(i,kb1)+r0(i,k1))
2452 
2453  rcv(i,kb2) = 0.5*(rcv(i,kb1)+rcv(i,k1))
2454  t(i,kb2) = 0.5*(t(i,kb1)+t(i,k1))
2455  s(i,kb2) = 0.5*(s(i,kb1)+s(i,k1))
2456  endif ; endif
2457  endif ! (h2 = 0 && h1 > 0)
2458 
2459  dpe_extrap = 0.0 ; dpe_merge = 0.0
2460  mergeable_bl = .false.
2461  if ((h1 > 0.0) .and. (h2 > 0.0) .and. (h_to_bl > 0.0) .and. &
2462  (stable_rcv)) then
2463  ! Check whether it is permissible for the buffer layers to detrain
2464  ! into the interior isopycnal layers.
2465 
2466  ! Determine the layer that has the lightest target density that is
2467  ! denser than the lowermost buffer layer.
2468  do k1=kb2+1,nz ; if (rcvtgt(k1) >= rcv(i,kb2)) exit ; enddo ; k0 = k1-1
2469  dr1 = rcvtgt(k0)-rcv(i,kb1) ; dr2 = rcv(i,kb2)-rcvtgt(k0)
2470 
2471  ! Use an energy-balanced combination of downwind advection into the next
2472  ! denser interior layer and upwind advection from the upper buffer layer
2473  ! into the lower one, each with an energy change that equals that required
2474  ! to mix the detrained water with the upper buffer layer.
2475  h1_avail = h1 - max(0.0,h_min_bl-h_to_bl)
2476  if ((k1<=nz) .and. (h2 > h_min_bl) .and. (h1_avail > 0.0) .and. &
2477  (r0(i,kb1) < r0(i,kb2)) .and. (h_to_bl*r0(i,kb1) > r0_to_bl)) then
2478  drk1 = (rcvtgt(k1) - rcv(i,kb2)) * (r0(i,kb2) - r0(i,kb1)) / &
2479  (rcv(i,kb2) - rcv(i,kb1))
2480  b1 = drk1 / (r0(i,kb2) - r0(i,kb1))
2481  ! b1 = RcvTgt(k1) - Rcv(i,kb2)) / (Rcv(i,kb2) - Rcv(i,kb1))
2482 
2483  ! Apply several limits to the detrainment.
2484  ! Entrain less than the mass in h2, and keep the base of the buffer
2485  ! layers from becoming shallower than any neighbors.
2486  h2_to_k1 = min(h2 - h_min_bl, h2_to_k1_rem)
2487  ! Balance downwind advection of density into the layer below the
2488  ! buffer layers with upwind advection from the layer above.
2489  if (h2_to_k1*(h1_avail + b1*(h1_avail + h2)) > h2*h1_avail) &
2490  h2_to_k1 = (h2*h1_avail) / (h1_avail + b1*(h1_avail + h2))
2491  if (h2_to_k1*(drk1 * h2) > (h_to_bl*r0(i,kb1) - r0_to_bl) * h1) &
2492  h2_to_k1 = (h_to_bl*r0(i,kb1) - r0_to_bl) * h1 / (drk1 * h2)
2493 
2494  if ((k1==kb2+1) .and. (cs%BL_extrap_lim > 0.)) then
2495  ! Simply do not detrain very light water into the lightest isopycnal
2496  ! coordinate layers if the density jump is too large.
2497  drcv_lim = rcv(i,kb2)-rcv(i,0)
2498  do k=1,kb2 ; drcv_lim = max(drcv_lim, rcv(i,kb2)-rcv(i,k)) ; enddo
2499  drcv_lim = cs%BL_extrap_lim*drcv_lim
2500  if ((rcvtgt(k1) - rcv(i,kb2)) >= drcv_lim) then
2501  h2_to_k1 = 0.0
2502  elseif ((rcvtgt(k1) - rcv(i,kb2)) > 0.5*drcv_lim) then
2503  h2_to_k1 = h2_to_k1 * (2.0 - 2.0*((rcvtgt(k1) - rcv(i,kb2)) / drcv_lim))
2504  endif
2505  endif
2506 
2507  drcv = (rcvtgt(k1) - rcv(i,kb2))
2508 
2509  ! Use 2nd order upwind advection of spiciness, limited by the values
2510  ! in deeper thick layers to determine the detrained temperature and
2511  ! salinity.
2512  dspice_det = (ds_dt_gauge*drcv_ds(i)*(t(i,kb2)-t(i,kb1)) - &
2513  dt_ds_gauge*drcv_dt(i)*(s(i,kb2)-s(i,kb1))) * &
2514  (h2 - h2_to_k1) / (h1 + h2)
2515  dspice_lim = 0.0
2516  if (h(i,k1) > 10.0*angstrom) then
2517  dspice_lim = ds_dt_gauge*drcv_ds(i)*(t(i,k1)-t(i,kb2)) - &
2518  dt_ds_gauge*drcv_dt(i)*(s(i,k1)-s(i,kb2))
2519  if (dspice_det*dspice_lim <= 0.0) dspice_lim = 0.0
2520  endif
2521  if (k1<nz) then ; if (h(i,k1+1) > 10.0*angstrom) then
2522  dspice_lim2 = ds_dt_gauge*drcv_ds(i)*(t(i,k1+1)-t(i,kb2)) - &
2523  dt_ds_gauge*drcv_dt(i)*(s(i,k1+1)-s(i,kb2))
2524  if ((dspice_det*dspice_lim2 > 0.0) .and. &
2525  (abs(dspice_lim2) > abs(dspice_lim))) dspice_lim = dspice_lim2
2526  endif ; endif
2527  if (abs(dspice_det) > abs(dspice_lim)) dspice_det = dspice_lim
2528 
2529  i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
2530  t_det = t(i,kb2) + dt_ds_gauge * i_denom * &
2531  (dt_ds_gauge * drcv_dt(i) * drcv + drcv_ds(i) * dspice_det)
2532  s_det = s(i,kb2) + i_denom * &
2533  (drcv_ds(i) * drcv - dt_ds_gauge * drcv_dt(i) * dspice_det)
2534  ! The detrained values of R0 are based on changes in T and S.
2535  r0_det = r0(i,kb2) + (t_det-t(i,kb2)) * dr0_dt(i) + &
2536  (s_det-s(i,kb2)) * dr0_ds(i)
2537 
2538  if (cs%BL_extrap_lim >= 0.) then
2539  ! Only do this detrainment if the new layer's temperature and salinity
2540  ! are not too far outside of the range of previous values.
2541  if (h(i,k1) > 10.0*angstrom) then
2542  t_min = min(t(i,kb1), t(i,kb2), t(i,k1)) - cs%Allowed_T_chg
2543  t_max = max(t(i,kb1), t(i,kb2), t(i,k1)) + cs%Allowed_T_chg
2544  s_min = min(s(i,kb1), s(i,kb2), s(i,k1)) - cs%Allowed_S_chg
2545  s_max = max(s(i,kb1), s(i,kb2), s(i,k1)) + cs%Allowed_S_chg
2546  else
2547  t_min = min(t(i,kb1), t(i,kb2)) - cs%Allowed_T_chg
2548  t_max = max(t(i,kb1), t(i,kb2)) + cs%Allowed_T_chg
2549  s_min = min(s(i,kb1), s(i,kb2)) - cs%Allowed_S_chg
2550  s_max = max(s(i,kb1), s(i,kb2)) + cs%Allowed_S_chg
2551  endif
2552  ihk1 = 1.0 / (h(i,k1) + h2_to_k1)
2553  t_new = (h(i,k1)*t(i,k1) + h2_to_k1*t_det) * ihk1
2554  s_new = (h(i,k1)*s(i,k1) + h2_to_k1*s_det) * ihk1
2555  ! A less restrictive limit might be used here.
2556  if ((t_new < t_min) .or. (t_new > t_max) .or. &
2557  (s_new < s_min) .or. (s_new > s_max)) &
2558  h2_to_k1 = 0.0
2559  endif
2560 
2561  h1_to_h2 = b1*h2*h2_to_k1 / (h2 - (1.0+b1)*h2_to_k1)
2562 
2563  ihk1 = 1.0 / (h(i,k1) + h_neglect + h2_to_k1)
2564  ih2f = 1.0 / ((h(i,kb2) - h2_to_k1) + h1_to_h2)
2565 
2566  rcv(i,kb2) = ((h(i,kb2)*rcv(i,kb2) - h2_to_k1*rcvtgt(k1)) + &
2567  h1_to_h2*rcv(i,kb1))*ih2f
2568  rcv(i,k1) = ((h(i,k1)+h_neglect)*rcv(i,k1) + h2_to_k1*rcvtgt(k1)) * ihk1
2569 
2570  t(i,kb2) = ((h(i,kb2)*t(i,kb2) - h2_to_k1*t_det) + &
2571  h1_to_h2*t(i,kb1)) * ih2f
2572  t(i,k1) = ((h(i,k1)+h_neglect)*t(i,k1) + h2_to_k1*t_det) * ihk1
2573 
2574  s(i,kb2) = ((h(i,kb2)*s(i,kb2) - h2_to_k1*s_det) + &
2575  h1_to_h2*s(i,kb1)) * ih2f
2576  s(i,k1) = ((h(i,k1)+h_neglect)*s(i,k1) + h2_to_k1*s_det) * ihk1
2577 
2578  ! Changes in R0 are based on changes in T and S.
2579  r0(i,kb2) = ((h(i,kb2)*r0(i,kb2) - h2_to_k1*r0_det) + &
2580  h1_to_h2*r0(i,kb1)) * ih2f
2581  r0(i,k1) = ((h(i,k1)+h_neglect)*r0(i,k1) + h2_to_k1*r0_det) * ihk1
2582 
2583  h(i,kb1) = h(i,kb1) - h1_to_h2 ; h1 = h(i,kb1)
2584  h(i,kb2) = (h(i,kb2) - h2_to_k1) + h1_to_h2 ; h2 = h(i,kb2)
2585  h(i,k1) = h(i,k1) + h2_to_k1
2586 
2587  d_ea(i,kb1) = d_ea(i,kb1) - h1_to_h2
2588  d_ea(i,kb2) = (d_ea(i,kb2) - h2_to_k1) + h1_to_h2
2589  d_ea(i,k1) = d_ea(i,k1) + h2_to_k1
2590  h2_to_k1_rem = max(h2_to_k1_rem - h2_to_k1, 0.0)
2591 
2592  ! The lower buffer layer has become lighter - it may be necessary to
2593  ! adjust k1 lighter.
2594  if ((k1>kb2+1) .and. (rcvtgt(k1-1) >= rcv(i,kb2))) then
2595  do k1=k1,kb2+1,-1 ; if (rcvtgt(k1-1) < rcv(i,kb2)) exit ; enddo
2596  endif
2597  endif
2598 
2599  k0 = k1-1
2600  dr1 = rcvtgt(k0)-rcv(i,kb1) ; dr2 = rcv(i,kb2)-rcvtgt(k0)
2601 
2602  if ((k0>kb2) .and. (dr1 > 0.0) .and. (h1 > h_min_bl) .and. &
2603  (h2*dr2 < h1*dr1) .and. (r0(i,kb2) > r0(i,kb1))) then
2604  ! An interior isopycnal layer (k0) is intermediate in density between
2605  ! the two buffer layers, and there can be detrainment. The entire
2606  ! lower buffer layer is combined with a portion of the upper buffer
2607  ! layer to match the target density of layer k0.
2608  stays_merge = 2.0*(h1+h2)*(h1*dr1 - h2*dr2) / &
2609  ((dr1+dr2)*h1 + dr1*(h1+h2) + &
2610  sqrt((dr2*h1-dr1*h2)**2 + 4*(h1+h2)*h2*(dr1+dr2)*dr2))
2611 
2612  stays_min_merge = max(h_min_bl, 2.0*h_min_bl - h_to_bl, &
2613  h1 - (h1+h2)*(r0(i,kb1) - r0_det) / (r0(i,kb2) - r0(i,kb1)))
2614  if ((stays_merge > stays_min_merge) .and. &
2615  (stays_merge + h2_to_k1_rem >= h1 + h2)) then
2616  mergeable_bl = .true.
2617  dpe_merge = g_2*(r0(i,kb2)-r0(i,kb1))*(h1-stays_merge)*(h2-stays_merge)
2618  endif
2619  endif
2620 
2621  if ((k1<=nz).and.(.not.mergeable_bl)) then
2622  ! Check whether linear extrapolation of density (i.e. 2nd order upwind
2623  ! advection) will allow some of the lower buffer layer to detrain into
2624  ! the next denser interior layer (k1).
2625  dr2b = rcvtgt(k1)-rcv(i,kb2) ; dr21 = rcv(i,kb2) - rcv(i,kb1)
2626  if (dr2b*(h1+h2) < h2*dr21) then
2627  ! Some of layer kb2 is denser than k1.
2628  h2_to_k1 = min(h2 - (h1+h2) * dr2b / dr21, h2_to_k1_rem)
2629 
2630  if (h2 > h2_to_k1) then
2631  drcv = (rcvtgt(k1) - rcv(i,kb2))
2632 
2633  ! Use 2nd order upwind advection of spiciness, limited by the values
2634  ! in deeper thick layers to determine the detrained temperature and
2635  ! salinity.
2636  dspice_det = (ds_dt_gauge*drcv_ds(i)*(t(i,kb2)-t(i,kb1)) - &
2637  dt_ds_gauge*drcv_dt(i)*(s(i,kb2)-s(i,kb1))) * &
2638  (h2 - h2_to_k1) / (h1 + h2)
2639  dspice_lim = 0.0
2640  if (h(i,k1) > 10.0*angstrom) then
2641  dspice_lim = ds_dt_gauge*drcv_ds(i)*(t(i,k1)-t(i,kb2)) - &
2642  dt_ds_gauge*drcv_dt(i)*(s(i,k1)-s(i,kb2))
2643  if (dspice_det*dspice_lim <= 0.0) dspice_lim = 0.0
2644  endif
2645  if (k1<nz) then; if (h(i,k1+1) > 10.0*angstrom) then
2646  dspice_lim2 = ds_dt_gauge*drcv_ds(i)*(t(i,k1+1)-t(i,kb2)) - &
2647  dt_ds_gauge*drcv_dt(i)*(s(i,k1+1)-s(i,kb2))
2648  if ((dspice_det*dspice_lim2 > 0.0) .and. &
2649  (abs(dspice_lim2) > abs(dspice_lim))) dspice_lim = dspice_lim2
2650  endif; endif
2651  if (abs(dspice_det) > abs(dspice_lim)) dspice_det = dspice_lim
2652 
2653  i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
2654  t_det = t(i,kb2) + dt_ds_gauge * i_denom * &
2655  (dt_ds_gauge * drcv_dt(i) * drcv + drcv_ds(i) * dspice_det)
2656  s_det = s(i,kb2) + i_denom * &
2657  (drcv_ds(i) * drcv - dt_ds_gauge * drcv_dt(i) * dspice_det)
2658  ! The detrained values of R0 are based on changes in T and S.
2659  r0_det = r0(i,kb2) + (t_det-t(i,kb2)) * dr0_dt(i) + &
2660  (s_det-s(i,kb2)) * dr0_ds(i)
2661 
2662  ! Now that the properties of the detrained water are known,
2663  ! potentially limit the amount of water that is detrained to
2664  ! avoid creating unphysical properties in the remaining water.
2665  ih2f = 1.0 / (h2 - h2_to_k1)
2666 
2667  t_min = min(t(i,kb2), t(i,kb1)) - cs%Allowed_T_chg
2668  t_max = max(t(i,kb2), t(i,kb1)) + cs%Allowed_T_chg
2669  t_new = (h2*t(i,kb2) - h2_to_k1*t_det)*ih2f
2670  if (t_new < t_min) then
2671  h2_to_k1_lim = h2 * (t(i,kb2) - t_min) / (t_det - t_min)
2672 ! write(mesg,'("Low temperature limits det to ", &
2673 ! & 1pe12.5, " from ", 1pe12.5, " at ", 1pg11.4,"E, ",1pg11.4,"N. T=", &
2674 ! & 5(1pe12.5))') &
2675 ! h2_to_k1_lim, h2_to_k1, G%geoLonT(i,j), G%geoLatT(i,j), &
2676 ! T_new, T(i,kb2), T(i,kb1), T_det, T_new-T_min
2677 ! call MOM_error(WARNING, mesg)
2678  h2_to_k1 = h2_to_k1_lim
2679  ih2f = 1.0 / (h2 - h2_to_k1)
2680  elseif (t_new > t_max) then
2681  h2_to_k1_lim = h2 * (t(i,kb2) - t_max) / (t_det - t_max)
2682 ! write(mesg,'("High temperature limits det to ", &
2683 ! & 1pe12.5, " from ", 1pe12.5, " at ", 1pg11.4,"E, ",1pg11.4,"N. T=", &
2684 ! & 5(1pe12.5))') &
2685 ! h2_to_k1_lim, h2_to_k1, G%geoLonT(i,j), G%geoLatT(i,j), &
2686 ! T_new, T(i,kb2), T(i,kb1), T_det, T_new-T_max
2687 ! call MOM_error(WARNING, mesg)
2688  h2_to_k1 = h2_to_k1_lim
2689  ih2f = 1.0 / (h2 - h2_to_k1)
2690  endif
2691  s_min = max(min(s(i,kb2), s(i,kb1)) - cs%Allowed_S_chg, 0.0)
2692  s_max = max(s(i,kb2), s(i,kb1)) + cs%Allowed_S_chg
2693  s_new = (h2*s(i,kb2) - h2_to_k1*s_det)*ih2f
2694  if (s_new < s_min) then
2695  h2_to_k1_lim = h2 * (s(i,kb2) - s_min) / (s_det - s_min)
2696 ! write(mesg,'("Low salinity limits det to ", &
2697 ! & 1pe12.5, " from ", 1pe12.5, " at ", 1pg11.4,"E, ",1pg11.4,"N. S=", &
2698 ! & 5(1pe12.5))') &
2699 ! h2_to_k1_lim, h2_to_k1, G%geoLonT(i,j), G%geoLatT(i,j), &
2700 ! S_new, S(i,kb2), S(i,kb1), S_det, S_new-S_min
2701 ! call MOM_error(WARNING, mesg)
2702  h2_to_k1 = h2_to_k1_lim
2703  ih2f = 1.0 / (h2 - h2_to_k1)
2704  elseif (s_new > s_max) then
2705  h2_to_k1_lim = h2 * (s(i,kb2) - s_max) / (s_det - s_max)
2706 ! write(mesg,'("High salinity limits det to ", &
2707 ! & 1pe12.5, " from ", 1pe12.5, " at ", 1pg11.4,"E, ",1pg11.4,"N. S=", &
2708 ! & 5(1pe12.5))') &
2709 ! h2_to_k1_lim, h2_to_k1, G%geoLonT(i,j), G%geoLatT(i,j), &
2710 ! S_new, S(i,kb2), S(i,kb1), S_det, S_new-S_max
2711 ! call MOM_error(WARNING, mesg)
2712  h2_to_k1 = h2_to_k1_lim
2713  ih2f = 1.0 / (h2 - h2_to_k1)
2714  endif
2715 
2716  ihk1 = 1.0 / (h(i,k1) + h_neglect + h2_to_k1)
2717  rcv(i,k1) = ((h(i,k1)+h_neglect)*rcv(i,k1) + h2_to_k1*rcvtgt(k1)) * ihk1
2718  rcv(i,kb2) = rcv(i,kb2) - h2_to_k1*drcv*ih2f
2719 
2720  t(i,kb2) = (h2*t(i,kb2) - h2_to_k1*t_det)*ih2f
2721  t(i,k1) = ((h(i,k1)+h_neglect)*t(i,k1) + h2_to_k1*t_det) * ihk1
2722 
2723  s(i,kb2) = (h2*s(i,kb2) - h2_to_k1*s_det) * ih2f
2724  s(i,k1) = ((h(i,k1)+h_neglect)*s(i,k1) + h2_to_k1*s_det) * ihk1
2725 
2726  ! Changes in R0 are based on changes in T and S.
2727  r0(i,kb2) = (h2*r0(i,kb2) - h2_to_k1*r0_det) * ih2f
2728  r0(i,k1) = ((h(i,k1)+h_neglect)*r0(i,k1) + h2_to_k1*r0_det) * ihk1
2729  else
2730  ! h2==h2_to_k1 can happen if dR2b = 0 exactly, but this is very
2731  ! unlikely. In this case the entirety of layer kb2 is detrained.
2732  h2_to_k1 = h2 ! These 2 lines are probably unnecessary.
2733  ihk1 = 1.0 / (h(i,k1) + h2)
2734 
2735  rcv(i,k1) = (h(i,k1)*rcv(i,k1) + h2*rcv(i,kb2)) * ihk1
2736  t(i,k1) = (h(i,k1)*t(i,k1) + h2*t(i,kb2)) * ihk1
2737  s(i,k1) = (h(i,k1)*s(i,k1) + h2*s(i,kb2)) * ihk1
2738  r0(i,k1) = (h(i,k1)*r0(i,k1) + h2*r0(i,kb2)) * ihk1
2739  endif
2740 
2741  h(i,k1) = h(i,k1) + h2_to_k1
2742  h(i,kb2) = h(i,kb2) - h2_to_k1 ; h2 = h(i,kb2)
2743  ! dPE_extrap should be positive here.
2744  dpe_extrap = i2rho0*(r0_det-r0(i,kb2))*h2_to_k1*h2
2745 
2746  d_ea(i,kb2) = d_ea(i,kb2) - h2_to_k1
2747  d_ea(i,k1) = d_ea(i,k1) + h2_to_k1
2748  h2_to_k1_rem = max(h2_to_k1_rem - h2_to_k1, 0.0)
2749  endif
2750  endif ! Detrainment by extrapolation.
2751 
2752  endif ! Detrainment to the interior at all.
2753 
2754  ! Does some of the detrained water go into the lower buffer layer?
2755  h_det_h2 = max(h_min_bl-(h1+h2), 0.0)
2756  if (h_det_h2 > 0.0) then
2757  ! Detrained water will go into both upper and lower buffer layers.
2758  ! h(kb2) will be h_min_bl, but h(kb1) may be larger if there was already
2759  ! ample detrainment; all water in layer kb1 moves into layer kb2.
2760 
2761  ! Determine the fluxes between the various layers.
2762  h_det_to_h2 = min(h_to_bl, h_det_h2)
2763  h_ml_to_h2 = h_det_h2 - h_det_to_h2
2764  h_det_to_h1 = h_to_bl - h_det_to_h2
2765  h_ml_to_h1 = max(h_min_bl-h_det_to_h1,0.0)
2766 
2767  ih = 1.0/h_min_bl
2768  ihdet = 0.0 ; if (h_to_bl > 0.0) ihdet = 1.0 / h_to_bl
2769  ih1f = 1.0 / (h_det_to_h1 + h_ml_to_h1)
2770 
2771  r0(i,kb2) = ((h2*r0(i,kb2) + h1*r0(i,kb1)) + &
2772  (h_det_to_h2*r0_to_bl*ihdet + h_ml_to_h2*r0(i,0))) * ih
2773  r0(i,kb1) = (h_det_to_h1*r0_to_bl*ihdet + h_ml_to_h1*r0(i,0)) * ih1f
2774 
2775  rcv(i,kb2) = ((h2*rcv(i,kb2) + h1*rcv(i,kb1)) + &
2776  (h_det_to_h2*rcv_to_bl*ihdet + h_ml_to_h2*rcv(i,0))) * ih
2777  rcv(i,kb1) = (h_det_to_h1*rcv_to_bl*ihdet + h_ml_to_h1*rcv(i,0)) * ih1f
2778 
2779  t(i,kb2) = ((h2*t(i,kb2) + h1*t(i,kb1)) + &
2780  (h_det_to_h2*t_to_bl*ihdet + h_ml_to_h2*t(i,0))) * ih
2781  t(i,kb1) = (h_det_to_h1*t_to_bl*ihdet + h_ml_to_h1*t(i,0)) * ih1f
2782 
2783  s(i,kb2) = ((h2*s(i,kb2) + h1*s(i,kb1)) + &
2784  (h_det_to_h2*s_to_bl*ihdet + h_ml_to_h2*s(i,0))) * ih
2785  s(i,kb1) = (h_det_to_h1*s_to_bl*ihdet + h_ml_to_h1*s(i,0)) * ih1f
2786 
2787  ! Recall that h1 = h(i,kb1) & h2 = h(i,kb2).
2788  d_ea(i,1) = d_ea(i,1) - (h_ml_to_h1 + h_ml_to_h2)
2789  d_ea(i,kb1) = d_ea(i,kb1) + ((h_det_to_h1 + h_ml_to_h1) - h1)
2790  d_ea(i,kb2) = d_ea(i,kb2) + (h_min_bl - h2)
2791 
2792  h(i,kb1) = h_det_to_h1 + h_ml_to_h1 ; h(i,kb2) = h_min_bl
2793  h(i,0) = h(i,0) - (h_ml_to_h1 + h_ml_to_h2)
2794 
2795 
2796  if (allocated(cs%diag_PE_detrain) .or. allocated(cs%diag_PE_detrain2)) then
2797  r0_det = r0_to_bl*ihdet
2798  s1en = g_2 * idt_h2 * ( ((r0(i,kb2)-r0(i,kb1))*h1*h2 + &
2799  h_det_to_h2*( (r0(i,kb1)-r0_det)*h1 + (r0(i,kb2)-r0_det)*h2 ) + &
2800  h_ml_to_h2*( (r0(i,kb2)-r0(i,0))*h2 + (r0(i,kb1)-r0(i,0))*h1 + &
2801  (r0_det-r0(i,0))*h_det_to_h2 ) + &
2802  h_det_to_h1*h_ml_to_h1*(r0_det-r0(i,0))) - 2.0*gv%Rho0*dpe_extrap )
2803 
2804  if (allocated(cs%diag_PE_detrain)) &
2805  cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + s1en
2806 
2807  if (allocated(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
2808  cs%diag_PE_detrain2(i,j) + s1en + idt_h2*rho0xg*dpe_extrap
2809  endif
2810 
2811  elseif ((h_to_bl > 0.0) .or. (h1 < h_min_bl) .or. (h2 < h_min_bl)) then
2812  ! Determine how much of the upper buffer layer will be moved into
2813  ! the lower buffer layer and the properties with which it is moving.
2814  ! This implementation assumes a 2nd-order upwind advection of density
2815  ! from the uppermost buffer layer into the next one down.
2816  h_from_ml = h_min_bl + max(h_min_bl-h2,0.0) - h1 - h_to_bl
2817  if (h_from_ml > 0.0) then
2818  ! Some water needs to be moved from the mixed layer so that the upper
2819  ! (and perhaps lower) buffer layers exceed their minimum thicknesses.
2820  dpe_extrap = dpe_extrap - i2rho0*h_from_ml*(r0_to_bl - r0(i,0)*h_to_bl)
2821  r0_to_bl = r0_to_bl + h_from_ml*r0(i,0)
2822  rcv_to_bl = rcv_to_bl + h_from_ml*rcv(i,0)
2823  t_to_bl = t_to_bl + h_from_ml*t(i,0)
2824  s_to_bl = s_to_bl + h_from_ml*s(i,0)
2825 
2826  h_to_bl = h_to_bl + h_from_ml
2827  h(i,0) = h(i,0) - h_from_ml
2828  d_ea(i,1) = d_ea(i,1) - h_from_ml
2829  endif
2830 
2831  ! The absolute value should be unnecessary and 1e9 is just a large number.
2832  b1 = 1.0e9
2833  if (r0(i,kb2) - r0(i,kb1) > 1.0e-9*abs(r0(i,kb1) - r0_det)) &
2834  b1 = abs(r0(i,kb1) - r0_det) / (r0(i,kb2) - r0(i,kb1))
2835  stays_min = max((1.0-b1)*h1 - b1*h2, 0.0, h_min_bl - h_to_bl)
2836  stays_max = h1 - max(h_min_bl-h2,0.0)
2837 
2838  scale_slope = 1.0
2839  if (stays_max <= stays_min) then
2840  stays = stays_max
2841  mergeable_bl = .false.
2842  if (stays_max < h1) scale_slope = (h1 - stays_min) / (h1 - stays_max)
2843  else
2844  ! There are numerous temporary variables used here that should not be
2845  ! used outside of this "else" branch: s1, s2, s3sq, I_ya, bh0
2846  bh0 = b1*h_to_bl
2847  i_ya = (h1 + h2) / ((h1 + h2) + h_to_bl)
2848  ! s1 is the amount staying that minimizes the PE increase.
2849  s1 = 0.5*(h1 + (h2 - bh0) * i_ya) ; s2 = h1 - s1
2850 
2851  if (s2 < 0.0) then
2852  ! The energy released by detrainment from the lower buffer layer can be
2853  ! used to mix water from the upper buffer layer into the lower one.
2854  s3sq = i_ya*max(bh0*h1-dpe_extrap, 0.0)
2855  else
2856  s3sq = i_ya*(bh0*h1-min(dpe_extrap,0.0))
2857  endif
2858 
2859  if (s3sq == 0.0) then
2860  ! There is a simple, exact solution to the quadratic equation, namely:
2861  stays = h1 ! This will revert to stays_max later.
2862  elseif (s2*s2 <= s3sq) then
2863  ! There is no solution with 0 PE change - use the minimum energy input.
2864  stays = s1
2865  else
2866  ! The following choose the solutions that are continuous with all water
2867  ! staying in the upper buffer layer when there is no detrainment,
2868  ! namely the + root when s2>0 and the - root otherwise. They also
2869  ! carefully avoid differencing large numbers, using s2 = (h1-s).
2870  if (bh0 <= 0.0) then ; stays = h1
2871  elseif (s2 > 0.0) then
2872  ! stays = s + sqrt(s2*s2 - s3sq) ! Note that s2 = h1-s
2873  if (s1 >= stays_max) then ; stays = stays_max
2874  elseif (s1 >= 0.0) then ; stays = s1 + sqrt(s2*s2 - s3sq)
2875  else ; stays = (h1*(s2-s1) - s3sq) / (-s1 + sqrt(s2*s2 - s3sq))
2876  endif
2877  else
2878  ! stays = s - sqrt(s2*s2 - s3sq) ! Note that s2 = h1-s & stays_min >= 0
2879  if (s1 <= stays_min) then ; stays = stays_min
2880  else ; stays = (h1*(s1-s2) + s3sq) / (s1 + sqrt(s2*s2 - s3sq))
2881  endif
2882  endif
2883  endif
2884 
2885  ! Limit the amount that stays so that the motion of water is from the
2886  ! upper buffer layer into the lower, but no more than is in the upper
2887  ! layer, and the water left in the upper layer is no lighter than the
2888  ! detrained water.
2889  if (stays >= stays_max) then ; stays = stays_max
2890  elseif (stays < stays_min) then ; stays = stays_min
2891  endif
2892  endif
2893 
2894  dpe_det = g_2*((r0(i,kb1)*h_to_bl - r0_to_bl)*stays + &
2895  (r0(i,kb2)-r0(i,kb1)) * (h1-stays) * &
2896  (h2 - scale_slope*stays*((h1+h2)+h_to_bl)/(h1+h2)) ) - &
2897  rho0xg*dpe_extrap
2898 
2899  if (dpe_time_ratio*h_to_bl > h_to_bl+h(i,0)) then
2900  dpe_ratio = (h_to_bl+h(i,0)) / h_to_bl
2901  else
2902  dpe_ratio = dpe_time_ratio
2903  endif
2904 
2905  if ((mergeable_bl) .and. (num_events*dpe_ratio*dpe_det > dpe_merge)) then
2906  ! It is energetically preferable to merge the two buffer layers, detrain
2907  ! them into interior layer (k0), move the remaining upper buffer layer
2908  ! water into the lower buffer layer, and detrain undiluted into the
2909  ! upper buffer layer.
2910  h1_to_k0 = (h1-stays_merge)
2911  stays = max(h_min_bl-h_to_bl,0.0)
2912  h1_to_h2 = stays_merge - stays
2913 
2914  ihk0 = 1.0 / ((h1_to_k0 + h2) + h(i,k0))
2915  ih1f = 1.0 / (h_to_bl + stays); ih2f = 1.0 / h1_to_h2
2916  ih12 = 1.0 / (h1 + h2)
2917 
2918  drcv_2dz = (rcv(i,kb1) - rcv(i,kb2)) * ih12
2919  drcv_stays = drcv_2dz*(h1_to_k0 + h1_to_h2)
2920  drcv_det = - drcv_2dz*(stays + h1_to_h2)
2921  rcv(i,k0) = ((h1_to_k0*(rcv(i,kb1) + drcv_det) + &
2922  h2*rcv(i,kb2)) + h(i,k0)*rcv(i,k0)) * ihk0
2923  rcv(i,kb2) = rcv(i,kb1) + drcv_2dz*(h1_to_k0-stays)
2924  rcv(i,kb1) = (rcv_to_bl + stays*(rcv(i,kb1) + drcv_stays)) * ih1f
2925 
2926  ! Use 2nd order upwind advection of spiciness, limited by the value in
2927  ! the water from the mixed layer to determine the temperature and
2928  ! salinity of the water that stays in the buffer layers.
2929  i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
2930  dspice_2dz = (ds_dt_gauge*drcv_ds(i)*(t(i,kb1)-t(i,kb2)) - &
2931  dt_ds_gauge*drcv_dt(i)*(s(i,kb1)-s(i,kb2))) * ih12
2932  dspice_lim = (ds_dt_gauge*dr0_ds(i)*(t_to_bl-t(i,kb1)*h_to_bl) - &
2933  dt_ds_gauge*dr0_dt(i)*(s_to_bl-s(i,kb1)*h_to_bl)) / h_to_bl
2934  if (dspice_lim * dspice_2dz <= 0.0) dspice_2dz = 0.0
2935 
2936  if (stays > 0.0) then
2937  ! Limit the spiciness of the water that stays in the upper buffer layer.
2938  if (abs(dspice_lim) < abs(dspice_2dz*(h1_to_k0 + h1_to_h2))) &
2939  dspice_2dz = dspice_lim/(h1_to_k0 + h1_to_h2)
2940 
2941  dspice_stays = dspice_2dz*(h1_to_k0 + h1_to_h2)
2942  t_stays = t(i,kb1) + dt_ds_gauge * i_denom * &
2943  (dt_ds_gauge * drcv_dt(i) * drcv_stays + drcv_ds(i) * dspice_stays)
2944  s_stays = s(i,kb1) + i_denom * &
2945  (drcv_ds(i) * drcv_stays - dt_ds_gauge * drcv_dt(i) * dspice_stays)
2946  ! The values of R0 are based on changes in T and S.
2947  r0_stays = r0(i,kb1) + (t_stays-t(i,kb1)) * dr0_dt(i) + &
2948  (s_stays-s(i,kb1)) * dr0_ds(i)
2949  else
2950  ! Limit the spiciness of the water that moves into the lower buffer layer.
2951  if (abs(dspice_lim) < abs(dspice_2dz*h1_to_k0)) &
2952  dspice_2dz = dspice_lim/h1_to_k0
2953  ! These will be multiplied by 0 later.
2954  t_stays = 0.0 ; s_stays = 0.0 ; r0_stays = 0.0
2955  endif
2956 
2957  dspice_det = - dspice_2dz*(stays + h1_to_h2)
2958  t_det = t(i,kb1) + dt_ds_gauge * i_denom * &
2959  (dt_ds_gauge * drcv_dt(i) * drcv_det + drcv_ds(i) * dspice_det)
2960  s_det = s(i,kb1) + i_denom * &
2961  (drcv_ds(i) * drcv_det - dt_ds_gauge * drcv_dt(i) * dspice_det)
2962  ! The values of R0 are based on changes in T and S.
2963  r0_det = r0(i,kb1) + (t_det-t(i,kb1)) * dr0_dt(i) + &
2964  (s_det-s(i,kb1)) * dr0_ds(i)
2965 
2966  t(i,k0) = ((h1_to_k0*t_det + h2*t(i,kb2)) + h(i,k0)*t(i,k0)) * ihk0
2967  t(i,kb2) = (h1*t(i,kb1) - stays*t_stays - h1_to_k0*t_det) * ih2f
2968  t(i,kb1) = (t_to_bl + stays*t_stays) * ih1f
2969 
2970  s(i,k0) = ((h1_to_k0*s_det + h2*s(i,kb2)) + h(i,k0)*s(i,k0)) * ihk0
2971  s(i,kb2) = (h1*s(i,kb1) - stays*s_stays - h1_to_k0*s_det) * ih2f
2972  s(i,kb1) = (s_to_bl + stays*s_stays) * ih1f
2973 
2974  r0(i,k0) = ((h1_to_k0*r0_det + h2*r0(i,kb2)) + h(i,k0)*r0(i,k0)) * ihk0
2975  r0(i,kb2) = (h1*r0(i,kb1) - stays*r0_stays - h1_to_k0*r0_det) * ih2f
2976  r0(i,kb1) = (r0_to_bl + stays*r0_stays) * ih1f
2977 
2978 ! ! The following is 2nd-order upwind advection without limiters.
2979 ! dT_2dz = (T(i,kb1) - T(i,kb2)) * Ih12
2980 ! T(i,k0) = (h1_to_k0*(T(i,kb1) - dT_2dz*(stays+h1_to_h2)) + &
2981 ! h2*T(i,kb2) + h(i,k0)*T(i,k0)) * Ihk0
2982 ! T(i,kb2) = T(i,kb1) + dT_2dz*(h1_to_k0-stays)
2983 ! T(i,kb1) = (T_to_bl + stays*(T(i,kb1) + &
2984 ! dT_2dz*(h1_to_k0 + h1_to_h2))) * Ih1f
2985 ! dS_2dz = (S(i,kb1) - S(i,kb2)) * Ih12
2986 ! S(i,k0) = (h1_to_k0*(S(i,kb1) - dS_2dz*(stays+h1_to_h2)) + &
2987 ! h2*S(i,kb2) + h(i,k0)*S(i,k0)) * Ihk0
2988 ! S(i,kb2) = S(i,kb1) + dS_2dz*(h1_to_k0-stays)
2989 ! S(i,kb1) = (S_to_bl + stays*(S(i,kb1) + &
2990 ! dS_2dz*(h1_to_k0 + h1_to_h2))) * Ih1f
2991 ! dR0_2dz = (R0(i,kb1) - R0(i,kb2)) * Ih12
2992 ! R0(i,k0) = (h1_to_k0*(R0(i,kb1) - dR0_2dz*(stays+h1_to_h2)) + &
2993 ! h2*R0(i,kb2) + h(i,k0)*R0(i,k0)) * Ihk0
2994 ! R0(i,kb2) = R0(i,kb1) + dR0_2dz*(h1_to_k0-stays)
2995 ! R0(i,kb1) = (R0_to_bl + stays*(R0(i,kb1) + &
2996 ! dR0_2dz*(h1_to_k0 + h1_to_h2))) * Ih1f
2997 
2998  d_ea(i,kb1) = (d_ea(i,kb1) + h_to_bl) + (stays - h1)
2999  d_ea(i,kb2) = d_ea(i,kb2) + (h1_to_h2 - h2)
3000  d_ea(i,k0) = d_ea(i,k0) + (h1_to_k0 + h2)
3001 
3002  h(i,kb1) = stays + h_to_bl
3003  h(i,kb2) = h1_to_h2
3004  h(i,k0) = h(i,k0) + (h1_to_k0 + h2)
3005  if (allocated(cs%diag_PE_detrain)) &
3006  cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + idt_h2*dpe_merge
3007  if (allocated(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
3008  cs%diag_PE_detrain2(i,j) + idt_h2*(dpe_det + rho0xg*dpe_extrap)
3009  else ! Not mergeable_bl.
3010  ! There is no further detrainment from the buffer layers, and the
3011  ! upper buffer layer water is distributed optimally between the
3012  ! upper and lower buffer layer.
3013  h1_to_h2 = h1 - stays
3014  ih1f = 1.0 / (h_to_bl + stays) ; ih2f = 1.0 / (h2 + h1_to_h2)
3015  ih = 1.0 / (h1 + h2)
3016  dr0_2dz = (r0(i,kb1) - r0(i,kb2)) * ih
3017  r0(i,kb2) = (h2*r0(i,kb2) + h1_to_h2*(r0(i,kb1) - &
3018  scale_slope*dr0_2dz*stays)) * ih2f
3019  r0(i,kb1) = (r0_to_bl + stays*(r0(i,kb1) + &
3020  scale_slope*dr0_2dz*h1_to_h2)) * ih1f
3021 
3022  ! Use 2nd order upwind advection of spiciness, limited by the value
3023  ! in the detrained water to determine the detrained temperature and
3024  ! salinity.
3025  dr0 = scale_slope*dr0_2dz*h1_to_h2
3026  dspice_stays = (ds_dt_gauge*dr0_ds(i)*(t(i,kb1)-t(i,kb2)) - &
3027  dt_ds_gauge*dr0_dt(i)*(s(i,kb1)-s(i,kb2))) * &
3028  scale_slope*h1_to_h2 * ih
3029  if (h_to_bl > 0.0) then
3030  dspice_lim = (ds_dt_gauge*dr0_ds(i)*(t_to_bl-t(i,kb1)*h_to_bl) - &
3031  dt_ds_gauge*dr0_dt(i)*(s_to_bl-s(i,kb1)*h_to_bl)) /&
3032  h_to_bl
3033  else
3034  dspice_lim = ds_dt_gauge*dr0_ds(i)*(t(i,0)-t(i,kb1)) - &
3035  dt_ds_gauge*dr0_dt(i)*(s(i,0)-s(i,kb1))
3036  endif
3037  if (dspice_stays*dspice_lim <= 0.0) then
3038  dspice_stays = 0.0
3039  elseif (abs(dspice_stays) > abs(dspice_lim)) then
3040  dspice_stays = dspice_lim
3041  endif
3042  i_denom = 1.0 / (dr0_ds(i)**2 + (dt_ds_gauge*dr0_dt(i))**2)
3043  t_stays = t(i,kb1) + dt_ds_gauge * i_denom * &
3044  (dt_ds_gauge * dr0_dt(i) * dr0 + dr0_ds(i) * dspice_stays)
3045  s_stays = s(i,kb1) + i_denom * &
3046  (dr0_ds(i) * dr0 - dt_ds_gauge * dr0_dt(i) * dspice_stays)
3047  ! The detrained values of Rcv are based on changes in T and S.
3048  rcv_stays = rcv(i,kb1) + (t_stays-t(i,kb1)) * drcv_dt(i) + &
3049  (s_stays-s(i,kb1)) * drcv_ds(i)
3050 
3051  t(i,kb2) = (h2*t(i,kb2) + h1*t(i,kb1) - t_stays*stays) * ih2f
3052  t(i,kb1) = (t_to_bl + stays*t_stays) * ih1f
3053  s(i,kb2) = (h2*s(i,kb2) + h1*s(i,kb1) - s_stays*stays) * ih2f
3054  s(i,kb1) = (s_to_bl + stays*s_stays) * ih1f
3055  rcv(i,kb2) = (h2*rcv(i,kb2) + h1*rcv(i,kb1) - rcv_stays*stays) * ih2f
3056  rcv(i,kb1) = (rcv_to_bl + stays*rcv_stays) * ih1f
3057 
3058 ! ! The following is 2nd-order upwind advection without limiters.
3059 ! dRcv_2dz = (Rcv(i,kb1) - Rcv(i,kb2)) * Ih
3060 ! dRcv = scale_slope*dRcv_2dz*h1_to_h2
3061 ! Rcv(i,kb2) = (h2*Rcv(i,kb2) + h1_to_h2*(Rcv(i,kb1) - &
3062 ! scale_slope*dRcv_2dz*stays)) * Ih2f
3063 ! Rcv(i,kb1) = (Rcv_to_bl + stays*(Rcv(i,kb1) + dRcv)) * Ih1f
3064 ! dT_2dz = (T(i,kb1) - T(i,kb2)) * Ih
3065 ! T(i,kb2) = (h2*T(i,kb2) + h1_to_h2*(T(i,kb1) - &
3066 ! scale_slope*dT_2dz*stays)) * Ih2f
3067 ! T(i,kb1) = (T_to_bl + stays*(T(i,kb1) + &
3068 ! scale_slope*dT_2dz*h1_to_h2)) * Ih1f
3069 ! dS_2dz = (S(i,kb1) - S(i,kb2)) * Ih
3070 ! S(i,kb2) = (h2*S(i,kb2) + h1_to_h2*(S(i,kb1) - &
3071 ! scale_slope*dS_2dz*stays)) * Ih2f
3072 ! S(i,kb1) = (S_to_bl + stays*(S(i,kb1) + &
3073 ! scale_slope*dS_2dz*h1_to_h2)) * Ih1f
3074 
3075  d_ea(i,kb1) = d_ea(i,kb1) + ((stays - h1) + h_to_bl)
3076  d_ea(i,kb2) = d_ea(i,kb2) + h1_to_h2
3077 
3078  h(i,kb1) = stays + h_to_bl
3079  h(i,kb2) = h(i,kb2) + h1_to_h2
3080 
3081  if (allocated(cs%diag_PE_detrain)) &
3082  cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + idt_h2*dpe_det
3083  if (allocated(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
3084  cs%diag_PE_detrain2(i,j) + idt_h2*(dpe_det + rho0xg*dpe_extrap)
3085  endif
3086  endif ! End of detrainment...
3087 
3088  enddo ! i loop
3089 

References mom_error_handler::mom_error().

Referenced by bulkmixedlayer().

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( g %isd: g %ied,0: gv %ke), intent(inout)  h,
real, dimension( g %isd: g %ied,0: gv %ke), intent(inout)  T,
real, dimension( g %isd: g %ied,0: gv %ke), intent(inout)  S,
real, dimension( g %isd: g %ied,0: gv %ke), intent(inout)  R0,
real, dimension( g %isd: g %ied,0: gv %ke), intent(inout)  Rcv,
real, dimension( gv %ke), intent(in)  RcvTgt,
real, dimension( g %isd: g %ied, gv %ke), intent(inout)  eps,
real, dimension( g %isd: g %ied, gv %ke), intent(inout)  d_ea,
real, dimension( g %isd: g %ied, gv %ke), intent(inout)  d_eb,
integer, dimension( g %isd: g %ied, gv %ke), intent(in)  ksort,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(bulkmixedlayer_cs), pointer  CS,
real, dimension( g %isd: g %ied), intent(in)  dR0_dT,
real, dimension( g %isd: g %ied), intent(in)  dR0_dS,
real, dimension( g %isd: g %ied), intent(in)  dRcv_dT,
real, dimension( g %isd: g %ied), 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 [H ~> m or kg m-2]. Layer 0 is the new mixed layer.
[in,out]tLayer temperatures [degC].
[in,out]sLayer salinities [ppt].
[in,out]r0Potential density referenced to surface pressure [R ~> kg m-3].
[in,out]rcvThe coordinate defining potential density [R ~> kg m-3].
[in]rcvtgtThe target value of Rcv for each layer [R ~> kg m-3].
[in,out]epsThe (small) thickness that must remain in each layer [H ~> m or kg m-2].
[in,out]d_eaThe upward increase across a layer in the entrainment from above [H ~> m or kg m-2]. Positive d_ea goes with layer thickness increases.
[in,out]d_ebThe downward increase across a layer in the entrainment from below [H ~> m or kg m-2]. 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 [R degC-1 ~> kg m-3 degC-1].
[in]dr0_dsThe partial derivative of cpotential density referenced to the surface with salinity, [R ppt-1 ~> kg m-3 ppt-1].
[in]drcv_dtThe partial derivative of coordinate defining potential density with potential temperature [R degC-1 ~> kg m-3 degC-1].
[in]drcv_dsThe partial derivative of coordinate defining potential density with salinity, [R ppt-1 ~> kg m-3 ppt-1].

Definition at line 1885 of file MOM_bulk_mixed_layer.F90.

1885  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1886  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid
1887  !! structure.
1888  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2].
1889  !! Layer 0 is the new mixed layer.
1890  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: T !< Layer temperatures [degC].
1891  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: S !< Layer salinities [ppt].
1892  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: R0 !< Potential density referenced to
1893  !! surface pressure [R ~> kg m-3].
1894  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: Rcv !< The coordinate defining
1895  !! potential density [R ~> kg m-3].
1896  real, dimension(SZK_(GV)), intent(in) :: RcvTgt !< The target value of Rcv for each
1897  !! layer [R ~> kg m-3].
1898  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: eps !< The (small) thickness that must
1899  !! remain in each layer [H ~> m or kg m-2].
1900  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_ea !< The upward increase across a
1901  !! layer in the entrainment from
1902  !! above [H ~> m or kg m-2].
1903  !! Positive d_ea goes with layer
1904  !! thickness increases.
1905  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_eb !< The downward increase across a
1906  !! layer in the entrainment from
1907  !! below [H ~> m or kg m-2]. Positive values go
1908  !! with mass gain by a layer.
1909  integer, dimension(SZI_(G),SZK_(GV)), intent(in) :: ksort !< The density-sorted k-indicies.
1910  type(bulkmixedlayer_CS), pointer :: CS !< The control structure for this
1911  !! module.
1912  real, dimension(SZI_(G)), intent(in) :: dR0_dT !< The partial derivative of
1913  !! potential density referenced
1914  !! to the surface with potential
1915  !! temperature [R degC-1 ~> kg m-3 degC-1].
1916  real, dimension(SZI_(G)), intent(in) :: dR0_dS !< The partial derivative of
1917  !! cpotential density referenced
1918  !! to the surface with salinity,
1919  !! [R ppt-1 ~> kg m-3 ppt-1].
1920  real, dimension(SZI_(G)), intent(in) :: dRcv_dT !< The partial derivative of
1921  !! coordinate defining potential
1922  !! density with potential
1923  !! temperature [R degC-1 ~> kg m-3 degC-1].
1924  real, dimension(SZI_(G)), intent(in) :: dRcv_dS !< The partial derivative of
1925  !! coordinate defining potential
1926  !! density with salinity,
1927  !! [R ppt-1 ~> kg m-3 ppt-1].
1928 
1929 ! If there are no massive light layers above the deepest of the mixed- and
1930 ! buffer layers, do nothing (except perhaps to reshuffle these layers).
1931 ! If there are nkbl or fewer layers above the deepest mixed- or buffer-
1932 ! layers, move them (in sorted order) into the buffer layers, even if they
1933 ! were previously interior layers.
1934 ! If there are interior layers that are intermediate in density (both in-situ
1935 ! and the coordinate density (sigma-2)) between the newly forming mixed layer
1936 ! and a residual buffer- or mixed layer, and the number of massive layers above
1937 ! the deepest massive buffer or mixed layer is greater than nkbl, then split
1938 ! those buffer layers into peices that match the target density of the two
1939 ! nearest interior layers.
1940 ! Otherwise, if there are more than nkbl+1 remaining massive layers
1941 
1942  ! Local variables
1943  real :: h_move, h_tgt_old, I_hnew
1944  real :: dT_dS_wt2, dT_dR, dS_dR, I_denom
1945  real :: Rcv_int
1946  real :: T_up, S_up, R0_up, I_hup, h_to_up
1947  real :: T_dn, S_dn, R0_dn, I_hdn, h_to_dn
1948  real :: wt_dn
1949  real :: dR1, dR2
1950  real :: dPE, hmin, min_dPE, min_hmin
1951  real, dimension(SZK_(GV)) :: &
1952  h_tmp, R0_tmp, T_tmp, S_tmp, Rcv_tmp
1953  integer :: ks_min
1954  logical :: sorted, leave_in_layer
1955  integer :: ks_deep(SZI_(G)), k_count(SZI_(G)), ks2_reverse(SZI_(G), SZK_(GV))
1956  integer :: ks2(SZK_(GV))
1957  integer :: i, k, ks, is, ie, nz, k1, k2, k_tgt, k_src, k_int_top
1958  integer :: nks, nkmb, num_interior, top_interior_ks
1959 
1960  is = g%isc ; ie = g%iec ; nz = gv%ke
1961  nkmb = cs%nkml+cs%nkbl
1962 
1963  dt_ds_wt2 = cs%dT_dS_wt**2
1964 
1965 ! Find out how many massive layers are above the deepest buffer or mixed layer.
1966  do i=is,ie ; ks_deep(i) = -1 ; k_count(i) = 0 ; enddo
1967  do ks=nz,1,-1 ; do i=is,ie ; if (ksort(i,ks) > 0) then
1968  k = ksort(i,ks)
1969 
1970  if (h(i,k) > eps(i,k)) then
1971  if (ks_deep(i) == -1) then
1972  if (k <= nkmb) then
1973  ks_deep(i) = ks ; k_count(i) = k_count(i) + 1
1974  ks2_reverse(i,k_count(i)) = k
1975  endif
1976  else
1977  k_count(i) = k_count(i) + 1
1978  ks2_reverse(i,k_count(i)) = k
1979  endif
1980  endif
1981  endif ; enddo ; enddo
1982 
1983  do i=is,ie ; if (k_count(i) > 1) then
1984  ! This column might need to be reshuffled.
1985  nks = k_count(i)
1986 
1987  ! Put ks2 in the right order and check whether reshuffling is needed.
1988  sorted = .true.
1989  ks2(nks) = ks2_reverse(i,1)
1990  do ks=nks-1,1,-1
1991  ks2(ks) = ks2_reverse(i,1+nks-ks)
1992  if (ks2(ks) > ks2(ks+1)) sorted = .false.
1993  enddo
1994 
1995  ! Go to the next column of no reshuffling is needed.
1996  if (sorted) cycle
1997 
1998  ! Find out how many interior layers are being reshuffled. If none,
1999  ! then this is a simple swapping procedure.
2000  num_interior = 0 ; top_interior_ks = 0
2001  do ks=nks,1,-1 ; if (ks2(ks) > nkmb) then
2002  num_interior = num_interior+1 ; top_interior_ks = ks
2003  endif ; enddo
2004 
2005  if (num_interior >= 1) then
2006  ! Find the lightest interior layer with a target coordinate density
2007  ! greater than the newly forming mixed layer.
2008  do k=nkmb+1,nz ; if (rcv(i,0) < rcvtgt(k)) exit ; enddo
2009  k_int_top = k ; rcv_int = rcvtgt(k)
2010 
2011  i_denom = 1.0 / (drcv_ds(i)**2 + dt_ds_wt2*drcv_dt(i)**2)
2012  dt_dr = dt_ds_wt2*drcv_dt(i) * i_denom
2013  ds_dr = drcv_ds(i) * i_denom
2014 
2015 
2016  ! Examine whether layers can be split out of existence. Stop when there
2017  ! is a layer that cannot be handled this way, or when the topmost formerly
2018  ! interior layer has been dealt with.
2019  do ks = nks,top_interior_ks,-1
2020  k = ks2(ks)
2021  leave_in_layer = .false.
2022  if ((k > nkmb) .and. (rcv(i,k) <= rcvtgt(k))) then
2023  if (rcvtgt(k)-rcv(i,k) < cs%BL_split_rho_tol*(rcvtgt(k) - rcvtgt(k-1))) &
2024  leave_in_layer = .true.
2025  elseif (k > nkmb) then
2026  if (rcv(i,k)-rcvtgt(k) < cs%BL_split_rho_tol*(rcvtgt(k+1) - rcvtgt(k))) &
2027  leave_in_layer = .true.
2028  endif
2029 
2030  if (leave_in_layer) then
2031  ! Just drop this layer from the sorted list.
2032  nks = nks-1
2033  elseif (rcv(i,k) < rcv_int) then
2034  ! There is no interior layer with a target density that is intermediate
2035  ! between this layer and the mixed layer.
2036  exit
2037  else
2038  ! Try splitting the layer between two interior isopycnal layers.
2039  ! Find the target densities that bracket this layer.
2040  do k2=k_int_top+1,nz ; if (rcv(i,k) < rcvtgt(k2)) exit ; enddo
2041  if (k2>nz) exit
2042 
2043  ! This layer is bracketed in density between layers k2-1 and k2.
2044 
2045  dr1 = (rcvtgt(k2-1) - rcv(i,k)) ; dr2 = (rcvtgt(k2) - rcv(i,k))
2046  t_up = t(i,k) + dt_dr * dr1
2047  s_up = s(i,k) + ds_dr * dr1
2048  t_dn = t(i,k) + dt_dr * dr2
2049  s_dn = s(i,k) + ds_dr * dr2
2050 
2051  r0_up = r0(i,k) + (dt_dr*dr0_dt(i) + ds_dr*dr0_ds(i)) * dr1
2052  r0_dn = r0(i,k) + (dt_dr*dr0_dt(i) + ds_dr*dr0_ds(i)) * dr2
2053 
2054  ! Make sure the new properties are acceptable.
2055  if ((r0_up > r0(i,0)) .or. (r0_dn > r0(i,0))) &
2056  ! Avoid creating obviously unstable profiles.
2057  exit
2058 
2059  wt_dn = (rcv(i,k) - rcvtgt(k2-1)) / (rcvtgt(k2) - rcvtgt(k2-1))
2060  h_to_up = (h(i,k)-eps(i,k)) * (1.0 - wt_dn)
2061  h_to_dn = (h(i,k)-eps(i,k)) * wt_dn
2062 
2063  i_hup = 1.0 / (h(i,k2-1) + h_to_up)
2064  i_hdn = 1.0 / (h(i,k2) + h_to_dn)
2065  r0(i,k2-1) = (r0(i,k2)*h(i,k2-1) + r0_up*h_to_up) * i_hup
2066  r0(i,k2) = (r0(i,k2)*h(i,k2) + r0_dn*h_to_dn) * i_hdn
2067 
2068  t(i,k2-1) = (t(i,k2)*h(i,k2-1) + t_up*h_to_up) * i_hup
2069  t(i,k2) = (t(i,k2)*h(i,k2) + t_dn*h_to_dn) * i_hdn
2070  s(i,k2-1) = (s(i,k2)*h(i,k2-1) + s_up*h_to_up) * i_hup
2071  s(i,k2) = (s(i,k2)*h(i,k2) + s_dn*h_to_dn) * i_hdn
2072  rcv(i,k2-1) = (rcv(i,k2)*h(i,k2-1) + rcvtgt(k2-1)*h_to_up) * i_hup
2073  rcv(i,k2) = (rcv(i,k2)*h(i,k2) + rcvtgt(k2)*h_to_dn) * i_hdn
2074 
2075  h(i,k) = eps(i,k)
2076  h(i,k2) = h(i,k2) + h_to_dn
2077  h(i,k2-1) = h(i,k2-1) + h_to_up
2078 
2079  if (k > k2-1) then
2080  d_eb(i,k) = d_eb(i,k) - h_to_up
2081  d_eb(i,k2-1) = d_eb(i,k2-1) + h_to_up
2082  elseif (k < k2-1) then
2083  d_ea(i,k) = d_ea(i,k) - h_to_up
2084  d_ea(i,k2-1) = d_ea(i,k2-1) + h_to_up
2085  endif
2086  if (k > k2) then
2087  d_eb(i,k) = d_eb(i,k) - h_to_dn
2088  d_eb(i,k2) = d_eb(i,k2) + h_to_dn
2089  elseif (k < k2) then
2090  d_ea(i,k) = d_ea(i,k) - h_to_dn
2091  d_ea(i,k2) = d_ea(i,k2) + h_to_dn
2092  endif
2093  nks = nks-1
2094  endif
2095  enddo
2096  endif
2097 
2098  do while (nks > nkmb)
2099  ! Having already tried to move surface layers into the interior, there
2100  ! are still too many layers, and layers must be merged until nks=nkmb.
2101  ! Examine every merger of a layer with its neighbors, and merge the ones
2102  ! that increase the potential energy the least. If there are layers
2103  ! with (apparently?) negative potential energy change, choose the one
2104  ! with the smallest total thickness. Repeat until nkmb layers remain.
2105  ! Choose the smaller value for the remaining index for convenience.
2106 
2107  ks_min = -1 ; min_dpe = 1.0 ; min_hmin = 0.0
2108  do ks=1,nks-1
2109  k1 = ks2(ks) ; k2 = ks2(ks+1)
2110  dpe = max(0.0, (r0(i,k2)-r0(i,k1)) * h(i,k1) * h(i,k2))
2111  hmin = min(h(i,k1)-eps(i,k1), h(i,k2)-eps(i,k2))
2112  if ((ks_min < 0) .or. (dpe < min_dpe) .or. &
2113  ((dpe <= 0.0) .and. (hmin < min_hmin))) then
2114  ks_min = ks ; min_dpe = dpe ; min_hmin = hmin
2115  endif
2116  enddo
2117 
2118  ! Now merge the two layers that do the least damage.
2119  k1 = ks2(ks_min) ; k2 = ks2(ks_min+1)
2120  if (k1 < k2) then ; k_tgt = k1 ; k_src = k2
2121  else ; k_tgt = k2 ; k_src = k1 ; ks2(ks_min) = k2 ; endif
2122 
2123  h_tgt_old = h(i,k_tgt)
2124  h_move = h(i,k_src)-eps(i,k_src)
2125  h(i,k_src) = eps(i,k_src)
2126  h(i,k_tgt) = h(i,k_tgt) + h_move
2127  i_hnew = 1.0 / (h(i,k_tgt))
2128  r0(i,k_tgt) = (r0(i,k_tgt)*h_tgt_old + r0(i,k_src)*h_move) * i_hnew
2129 
2130  t(i,k_tgt) = (t(i,k_tgt)*h_tgt_old + t(i,k_src)*h_move) * i_hnew
2131  s(i,k_tgt) = (s(i,k_tgt)*h_tgt_old + s(i,k_src)*h_move) * i_hnew
2132  rcv(i,k_tgt) = (rcv(i,k_tgt)*h_tgt_old + rcv(i,k_src)*h_move) * i_hnew
2133 
2134  d_eb(i,k_src) = d_eb(i,k_src) - h_move
2135  d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_move
2136 
2137  ! Remove the newly missing layer from the sorted list.
2138  do ks=ks_min+1,nks ; ks2(ks) = ks2(ks+1) ; enddo
2139  nks = nks-1
2140  enddo
2141 
2142  ! Check again whether the layers are sorted, and go on to the next column
2143  ! if they are.
2144  sorted = .true.
2145  do ks=1,nks-1 ; if (ks2(ks) > ks2(ks+1)) sorted = .false. ; enddo
2146  if (sorted) cycle
2147 
2148  if (nks > 1) then
2149  ! At this point, actually swap the properties of the layers, and place
2150  ! the remaining layers in order starting with nkmb.
2151 
2152  ! Save all the properties of the nkmb layers that might be replaced.
2153  do k=1,nkmb
2154  h_tmp(k) = h(i,k) ; r0_tmp(k) = r0(i,k)
2155  t_tmp(k) = t(i,k) ; s_tmp(k) = s(i,k) ; rcv_tmp(k) = rcv(i,k)
2156 
2157  h(i,k) = 0.0
2158  enddo
2159 
2160  do ks=nks,1,-1
2161  k_tgt = nkmb - nks + ks ; k_src = ks2(ks)
2162  if (k_tgt == k_src) then
2163  h(i,k_tgt) = h_tmp(k_tgt) ! This layer doesn't move, so put the water back.
2164  cycle
2165  endif
2166 
2167  ! Note below that eps=0 for k<=nkmb.
2168  if (k_src > nkmb) then
2169  h_move = h(i,k_src)-eps(i,k_src)
2170  h(i,k_src) = eps(i,k_src)
2171  h(i,k_tgt) = h_move
2172  r0(i,k_tgt) = r0(i,k_src)
2173 
2174  t(i,k_tgt) = t(i,k_src) ; s(i,k_tgt) = s(i,k_src)
2175  rcv(i,k_tgt) = rcv(i,k_src)
2176 
2177  d_eb(i,k_src) = d_eb(i,k_src) - h_move
2178  d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_move
2179  else
2180  h(i,k_tgt) = h_tmp(k_src)
2181  r0(i,k_tgt) = r0_tmp(k_src)
2182 
2183  t(i,k_tgt) = t_tmp(k_src) ; s(i,k_tgt) = s_tmp(k_src)
2184  rcv(i,k_tgt) = rcv_tmp(k_src)
2185 
2186  if (k_src > k_tgt) then
2187  d_eb(i,k_src) = d_eb(i,k_src) - h_tmp(k_src)
2188  d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_tmp(k_src)
2189  else
2190  d_ea(i,k_src) = d_ea(i,k_src) - h_tmp(k_src)
2191  d_ea(i,k_tgt) = d_ea(i,k_tgt) + h_tmp(k_src)
2192  endif
2193  endif
2194  enddo
2195  endif
2196 
2197  endif ; enddo
2198 

Referenced by bulkmixedlayer().

Here is the caller graph for this function:

◆ sort_ml()

subroutine mom_bulk_mixed_layer::sort_ml ( real, dimension( g %isd: g %ied, gv %ke), intent(in)  h,
real, dimension( g %isd: g %ied, gv %ke), intent(in)  R0,
real, dimension( g %isd: g %ied, gv %ke), intent(in)  eps,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(bulkmixedlayer_cs), pointer  CS,
integer, dimension( g %isd: g %ied, gv %ke), 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 [H ~> m or kg m-2].
[in]r0The potential density used to sort the layers [R ~> kg m-3].
[in]epsThe (small) thickness that must remain in each layer [H ~> m or kg m-2].
csThe control structure returned by a previous call to mixedlayer_init.
[out]ksortThe k-index to use in the sort.

Definition at line 1833 of file MOM_bulk_mixed_layer.F90.

1833  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1834  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1835  real, dimension(SZI_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2].
1836  real, dimension(SZI_(G),SZK_(GV)), intent(in) :: R0 !< The potential density used to sort
1837  !! the layers [R ~> kg m-3].
1838  real, dimension(SZI_(G),SZK_(GV)), intent(in) :: eps !< The (small) thickness that must
1839  !! remain in each layer [H ~> m or kg m-2].
1840  type(bulkmixedlayer_CS), pointer :: CS !< The control structure returned by a
1841  !! previous call to mixedlayer_init.
1842  integer, dimension(SZI_(G),SZK_(GV)), intent(out) :: ksort !< The k-index to use in the sort.
1843 
1844  ! Local variables
1845  real :: R0sort(SZI_(G),SZK_(GV))
1846  integer :: nsort(SZI_(G))
1847  logical :: done_sorting(SZI_(G))
1848  integer :: i, k, ks, is, ie, nz, nkmb
1849 
1850  is = g%isc ; ie = g%iec ; nz = gv%ke
1851  nkmb = cs%nkml+cs%nkbl
1852 
1853  ! Come up with a sorted index list of layers with increasing R0.
1854  ! Assume that the layers below nkmb are already stably stratified.
1855  ! Only layers that are thicker than eps are in the list. Extra elements
1856  ! have an index of -1.
1857 
1858  ! This is coded using straight insertion, on the assumption that the
1859  ! layers are usually in the right order (or massless) anyway.
1860 
1861  do k=1,nz ; do i=is,ie ; ksort(i,k) = -1 ; enddo ; enddo
1862 
1863  do i=is,ie ; nsort(i) = 0 ; done_sorting(i) = .false. ; enddo
1864  do k=1,nz ; do i=is,ie ; if (h(i,k) > eps(i,k)) then
1865  if (done_sorting(i)) then ; ks = nsort(i) ; else
1866  do ks=nsort(i),1,-1
1867  if (r0(i,k) >= r0sort(i,ks)) exit
1868  r0sort(i,ks+1) = r0sort(i,ks) ; ksort(i,ks+1) = ksort(i,ks)
1869  enddo
1870  if ((k > nkmb) .and. (ks == nsort(i))) done_sorting(i) = .true.
1871  endif
1872 
1873  ksort(i,ks+1) = k
1874  r0sort(i,ks+1) = r0(i,k)
1875  nsort(i) = nsort(i) + 1
1876  endif ; enddo ; enddo
1877 

Referenced by bulkmixedlayer().

Here is the caller graph for this function:

Variable Documentation

◆ id_clock_adjustment

integer mom_bulk_mixed_layer::id_clock_adjustment =0
private

CPU clock IDs.

Definition at line 151 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

CPU clock IDs.

Definition at line 151 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

CPU clock IDs.

Definition at line 151 of file MOM_bulk_mixed_layer.F90.

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

Referenced by bulkmixedlayer(), and bulkmixedlayer_init().

◆ id_clock_eos

integer mom_bulk_mixed_layer::id_clock_eos =0
private

CPU clock IDs.

Definition at line 152 of file MOM_bulk_mixed_layer.F90.

152 integer :: id_clock_EOS=0, id_clock_resort=0, id_clock_pass=0

Referenced by bulkmixedlayer(), and bulkmixedlayer_init().

◆ id_clock_mech

integer mom_bulk_mixed_layer::id_clock_mech =0
private

CPU clock IDs.

Definition at line 151 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

CPU clock IDs.

Definition at line 152 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

CPU clock IDs.

Definition at line 152 of file MOM_bulk_mixed_layer.F90.

Referenced by bulkmixedlayer(), and bulkmixedlayer_init().