MOM6
Data Types
mom_energetic_pbl Module Reference

Detailed Description

Energetically consistent planetary boundary layer parameterization.

Data Types

type  energetic_pbl_cs
 This control structure holds parameters for the MOM_energetic_PBL module. More...
 
type  epbl_column_diags
 A type for conveniently passing around ePBL diagnostics for a column. More...
 
integer, parameter use_fixed_mstar = 0
 Enumeration values for mstar_Scheme. More...
 
integer, parameter mstar_from_ekman = 2
 The value of mstar_scheme to base mstar on the ratio of the Ekman layer depth to the Obukhov depth. More...
 
integer, parameter mstar_from_rh18 = 3
 The value of mstar_scheme to base mstar of of RH18. More...
 
integer, parameter no_langmuir = 0
 The value of LT_ENHANCE_FORM not use Langmuir turbolence. More...
 
integer, parameter langmuir_rescale = 2
 The value of LT_ENHANCE_FORM to use a multiplicative rescaling of mstar to account for Langmuir turbulence. More...
 
integer, parameter langmuir_add = 3
 The value of LT_ENHANCE_FORM to add a contribution to mstar from Langmuir turblence to other contributions. More...
 
integer, parameter wt_from_croot_tke = 0
 Use a constant times the cube root of remaining TKE to calculate the turbulent velocity. More...
 
integer, parameter wt_from_rh18 = 1
 Use a scheme based on a combination of w* and v* as documented in Reichl & Hallberg (2018) to calculate the turbulent velocity. More...
 
character *(20), parameter constant_string = "CONSTANT"
 Enumeration values for mstar_Scheme. More...
 
character *(20), parameter om4_string = "OM4"
 Enumeration values for mstar_Scheme. More...
 
character *(20), parameter rh18_string = "REICHL_H18"
 Enumeration values for mstar_Scheme. More...
 
character *(20), parameter root_tke_string = "CUBE_ROOT_TKE"
 Enumeration values for mstar_Scheme. More...
 
character *(20), parameter none_string = "NONE"
 Enumeration values for mstar_Scheme. More...
 
character *(20), parameter rescaled_string = "RESCALE"
 Enumeration values for mstar_Scheme. More...
 
character *(20), parameter additive_string = "ADDITIVE"
 Enumeration values for mstar_Scheme. More...
 
subroutine, public energetic_pbl (h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS, dSV_dT, dSV_dS, TKE_forced, buoy_flux, dt_diag, last_call, dT_expected, dS_expected, Waves)
 This subroutine determines the diffusivities from the integrated energetics mixed layer model. It assumes that heating, cooling and freshwater fluxes have already been applied. All calculations are done implicitly, and there is no stability limit on the time step. More...
 
subroutine epbl_column (h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, absf, u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, US, CS, eCD, dt_diag, Waves, G, i, j)
 This subroutine determines the diffusivities from the integrated energetics mixed layer model for a single column of water. More...
 
subroutine find_pe_chg (Kddt_h0, dKddt_h, hp_a, hp_b, Th_a, Sh_a, Th_b, Sh_b, dT_to_dPE_a, dS_to_dPE_a, dT_to_dPE_b, dS_to_dPE_b, pres_Z, dT_to_dColHt_a, dS_to_dColHt_a, dT_to_dColHt_b, dS_to_dColHt_b, PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0, PE_ColHt_cor)
 This subroutine calculates the change in potential energy and or derivatives for several changes in an interfaces's diapycnal diffusivity times a timestep. More...
 
subroutine find_pe_chg_orig (Kddt_h, h_k, b_den_1, dTe_term, dSe_term, dT_km1_t2, dS_km1_t2, dT_to_dPE_k, dS_to_dPE_k, dT_to_dPEa, dS_to_dPEa, pres_Z, dT_to_dColHt_k, dS_to_dColHt_k, dT_to_dColHta, dS_to_dColHta, PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0)
 This subroutine calculates the change in potential energy and or derivatives for several changes in an interfaces's diapycnal diffusivity times a timestep using the original form used in the first version of ePBL. More...
 
subroutine find_mstar (CS, US, Buoyancy_Flux, UStar, UStar_Mean, BLD, Abs_Coriolis, MStar, Langmuir_Number, MStar_LT, Convect_Langmuir_Number)
 This subroutine finds the Mstar value for ePBL. More...
 
subroutine mstar_langmuir (CS, US, Abs_Coriolis, Buoyancy_Flux, UStar, BLD, Langmuir_Number, Mstar, MStar_LT, Convect_Langmuir_Number)
 This subroutine modifies the Mstar value if the Langmuir number is present. More...
 
subroutine, public energetic_pbl_get_mld (CS, MLD, G, US, m_to_MLD_units)
 Copies the ePBL active mixed layer depth into MLD. More...
 
subroutine, public energetic_pbl_init (Time, G, GV, US, param_file, diag, CS)
 This subroutine initializes the energetic_PBL module. More...
 
subroutine, public energetic_pbl_end (CS)
 Clean up and deallocate memory associated with the energetic_PBL module. More...
 

Function/Subroutine Documentation

◆ energetic_pbl()

subroutine, public mom_energetic_pbl::energetic_pbl ( 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+1), intent(out)  Kd_int,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(energetic_pbl_cs), pointer  CS,
real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)  dSV_dT,
real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)  dSV_dS,
real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(in)  TKE_forced,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  buoy_flux,
real, intent(in), optional  dt_diag,
logical, intent(in), optional  last_call,
real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(out), optional  dT_expected,
real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(out), optional  dS_expected,
type(wave_parameters_cs), optional, pointer  Waves 
)

This subroutine determines the diffusivities from the integrated energetics mixed layer model. It assumes that heating, cooling and freshwater fluxes have already been applied. All calculations are done implicitly, and there is no stability limit on the time step.

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 thicknesses [H ~> m or kg m-2].
[in]u_3dZonal velocities interpolated to h points
[in]v_3dZonal velocities interpolated to h points
[in]dsv_dtThe partial derivative of in-situ specific
[in]dsv_dsThe partial derivative of in-situ specific
[in]tke_forcedThe forcing requirements to homogenize the
[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].
[out]kd_intThe diagnosed diffusivities at interfaces
csThe control structure returned by a previous call to mixedlayer_init.
[in]buoy_fluxThe surface buoyancy flux [Z2 T-3 ~> m2 s-3].
[in]dt_diagThe diagnostic time step, which may be less than dt if there are two calls 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.
[out]dt_expectedThe values of temperature change that
[out]ds_expectedThe values of salinity change that
wavesWave CS

Definition at line 243 of file MOM_energetic_PBL.F90.

243  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
244  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
245  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
246  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
247  intent(inout) :: h_3d !< Layer thicknesses [H ~> m or kg m-2].
248  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
249  intent(in) :: u_3d !< Zonal velocities interpolated to h points
250  !! [L T-1 ~> m s-1].
251  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
252  intent(in) :: v_3d !< Zonal velocities interpolated to h points
253  !! [L T-1 ~> m s-1].
254  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
255  intent(in) :: dSV_dT !< The partial derivative of in-situ specific
256  !! volume with potential temperature
257  !! [R-1 degC-1 ~> m3 kg-1 degC-1].
258  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
259  intent(in) :: dSV_dS !< The partial derivative of in-situ specific
260  !! volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
261  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
262  intent(in) :: TKE_forced !< The forcing requirements to homogenize the
263  !! forcing that has been applied to each layer
264  !! [R Z3 T-2 ~> J m-2].
265  type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers to any
266  !! available thermodynamic fields. Absent fields
267  !! have NULL ptrs.
268  type(forcing), intent(inout) :: fluxes !< A structure containing pointers to any
269  !! possible forcing fields. Unused fields have
270  !! NULL ptrs.
271  real, intent(in) :: dt !< Time increment [T ~> s].
272  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
273  intent(out) :: Kd_int !< The diagnosed diffusivities at interfaces
274  !! [Z2 s-1 ~> m2 s-1].
275  type(energetic_PBL_CS), pointer :: CS !< The control structure returned by a previous
276  !! call to mixedlayer_init.
277  real, dimension(SZI_(G),SZJ_(G)), &
278  intent(in) :: buoy_flux !< The surface buoyancy flux [Z2 T-3 ~> m2 s-3].
279  real, optional, intent(in) :: dt_diag !< The diagnostic time step, which may be less
280  !! than dt if there are two calls to mixedlayer [T ~> s].
281  logical, optional, intent(in) :: last_call !< If true, this is the last call to
282  !! mixedlayer in the current time step, so
283  !! diagnostics will be written. The default
284  !! is .true.
285  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
286  optional, intent(out) :: dT_expected !< The values of temperature change that
287  !! should be expected when the returned
288  !! diffusivities are applied [degC].
289  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
290  optional, intent(out) :: dS_expected !< The values of salinity change that
291  !! should be expected when the returned
292  !! diffusivities are applied [ppt].
293  type(wave_parameters_CS), &
294  optional, pointer :: Waves !< Wave CS
295 
296 ! This subroutine determines the diffusivities from the integrated energetics
297 ! mixed layer model. It assumes that heating, cooling and freshwater fluxes
298 ! have already been applied. All calculations are done implicitly, and there
299 ! is no stability limit on the time step.
300 !
301 ! For each interior interface, first discard the TKE to account for mixing
302 ! of shortwave radiation through the next denser cell. Next drive mixing based
303 ! on the local? values of ustar + wstar, subject to available energy. This
304 ! step sets the value of Kd(K). Any remaining energy is then subject to decay
305 ! before being handed off to the next interface. mech_TKE and conv_PErel are treated
306 ! separately for the purposes of decay, but are used proportionately to drive
307 ! mixing.
308 !
309 ! The key parameters for the mixed layer are found in the control structure.
310 ! To use the classic constant mstar mixied layers choose MSTAR_SCHEME=CONSTANT.
311 ! The key parameters then include mstar, nstar, TKE_decay, and conv_decay.
312 ! For the Oberhuber (1993) mixed layer,the values of these are:
313 ! mstar = 1.25, nstar = 1, TKE_decay = 2.5, conv_decay = 0.5
314 ! TKE_decay is 1/kappa in eq. 28 of Oberhuber (1993), while conv_decay is 1/mu.
315 ! For a traditional Kraus-Turner mixed layer, the values are:
316 ! mstar = 1.25, nstar = 0.4, TKE_decay = 0.0, conv_decay = 0.0
317 
318  ! Local variables
319  real, dimension(SZI_(G),SZK_(GV)) :: &
320  h_2d, & ! A 2-d slice of the layer thickness [H ~> m or kg m-2].
321  T_2d, & ! A 2-d slice of the layer temperatures [degC].
322  S_2d, & ! A 2-d slice of the layer salinities [ppt].
323  TKE_forced_2d, & ! A 2-d slice of TKE_forced [R Z3 T-2 ~> J m-2].
324  dSV_dT_2d, & ! A 2-d slice of dSV_dT [R-1 degC-1 ~> m3 kg-1 degC-1].
325  dSV_dS_2d, & ! A 2-d slice of dSV_dS [R-1 ppt-1 ~> m3 kg-1 ppt-1].
326  u_2d, & ! A 2-d slice of the zonal velocity [L T-1 ~> m s-1].
327  v_2d ! A 2-d slice of the meridional velocity [L T-1 ~> m s-1].
328  real, dimension(SZI_(G),SZK_(GV)+1) :: &
329  Kd_2d ! A 2-d version of the diapycnal diffusivity [Z2 T-1 ~> m2 s-1].
330  real, dimension(SZK_(GV)) :: &
331  h, & ! The layer thickness [H ~> m or kg m-2].
332  T0, & ! The initial layer temperatures [degC].
333  S0, & ! The initial layer salinities [ppt].
334  dSV_dT_1d, & ! The partial derivatives of specific volume with temperature [R-1 degC-1 ~> m3 kg-1 degC-1].
335  dSV_dS_1d, & ! The partial derivatives of specific volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
336  TKE_forcing, & ! Forcing of the TKE in the layer coming from TKE_forced [R Z3 T-2 ~> J m-2].
337  u, & ! The zonal velocity [L T-1 ~> m s-1].
338  v ! The meridional velocity [L T-1 ~> m s-1].
339  real, dimension(SZK_(GV)+1) :: &
340  Kd, & ! The diapycnal diffusivity [Z2 T-1 ~> m2 s-1].
341  mixvel, & ! A turbulent mixing veloxity [Z T-1 ~> m s-1].
342  mixlen ! A turbulent mixing length [Z ~> m].
343  real :: h_neglect ! A thickness that is so small it is usually lost
344  ! in roundoff and can be neglected [H ~> m or kg m-2].
345 
346  real :: absf ! The absolute value of f [T-1 ~> s-1].
347  real :: U_star ! The surface friction velocity [Z T-1 ~> m s-1].
348  real :: U_Star_Mean ! The surface friction without gustiness [Z T-1 ~> m s-1].
349  real :: B_Flux ! The surface buoyancy flux [Z2 T-3 ~> m2 s-3]
350  real :: MLD_io ! The mixed layer depth found by ePBL_column [Z ~> m].
351 
352 ! The following are only used for diagnostics.
353  real :: dt__diag ! A copy of dt_diag (if present) or dt [T ~> s].
354  logical :: write_diags ! If true, write out diagnostics with this step.
355  logical :: reset_diags ! If true, zero out the accumulated diagnostics.
356 
357  logical :: debug=.false. ! Change this hard-coded value for debugging.
358  type(ePBL_column_diags) :: eCD ! A container for passing around diagnostics.
359 
360  integer :: i, j, k, is, ie, js, je, nz
361 
362  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
363 
364  if (.not. associated(cs)) call mom_error(fatal, "energetic_PBL: "//&
365  "Module must be initialized before it is used.")
366 
367  if (.not. associated(tv%eqn_of_state)) call mom_error(fatal, &
368  "energetic_PBL: Temperature, salinity and an equation of state "//&
369  "must now be used.")
370  if (.NOT. associated(fluxes%ustar)) call mom_error(fatal, &
371  "energetic_PBL: No surface TKE fluxes (ustar) defined in mixedlayer!")
372  debug = .false. ; if (present(dt_expected) .or. present(ds_expected)) debug = .true.
373 
374  if (debug) allocate(ecd%dT_expect(nz), ecd%dS_expect(nz))
375 
376  h_neglect = gv%H_subroundoff
377 
378  dt__diag = dt ; if (present(dt_diag)) dt__diag = dt_diag
379  write_diags = .true. ; if (present(last_call)) write_diags = last_call
380 
381 
382  ! Determine whether to zero out diagnostics before accumulation.
383  reset_diags = .true.
384  if (present(dt_diag) .and. write_diags .and. (dt__diag > dt)) &
385  reset_diags = .false. ! This is the second call to mixedlayer.
386 
387  if (reset_diags) then
388  if (cs%TKE_diagnostics) then
389 !!OMP parallel do default(none) shared(is,ie,js,je,CS)
390  do j=js,je ; do i=is,ie
391  cs%diag_TKE_wind(i,j) = 0.0 ; cs%diag_TKE_MKE(i,j) = 0.0
392  cs%diag_TKE_conv(i,j) = 0.0 ; cs%diag_TKE_forcing(i,j) = 0.0
393  cs%diag_TKE_mixing(i,j) = 0.0 ; cs%diag_TKE_mech_decay(i,j) = 0.0
394  cs%diag_TKE_conv_decay(i,j) = 0.0 !; CS%diag_TKE_unbalanced(i,j) = 0.0
395  enddo ; enddo
396  endif
397  endif
398  ! if (CS%id_Mixing_Length>0) CS%Mixing_Length(:,:,:) = 0.0
399  ! if (CS%id_Velocity_Scale>0) CS%Velocity_Scale(:,:,:) = 0.0
400 
401 !!OMP parallel do default(private) shared(js,je,nz,is,ie,h_3d,u_3d,v_3d,tv,dt, &
402 !!OMP CS,G,GV,US,fluxes,debug, &
403 !!OMP TKE_forced,dSV_dT,dSV_dS,Kd_int)
404  do j=js,je
405  ! Copy the thicknesses and other fields to 2-d arrays.
406  do k=1,nz ; do i=is,ie
407  h_2d(i,k) = h_3d(i,j,k) ; u_2d(i,k) = u_3d(i,j,k) ; v_2d(i,k) = v_3d(i,j,k)
408  t_2d(i,k) = tv%T(i,j,k) ; s_2d(i,k) = tv%S(i,j,k)
409  tke_forced_2d(i,k) = tke_forced(i,j,k)
410  dsv_dt_2d(i,k) = dsv_dt(i,j,k) ; dsv_ds_2d(i,k) = dsv_ds(i,j,k)
411  enddo ; enddo
412 
413  ! Determine the initial mech_TKE and conv_PErel, including the energy required
414  ! to mix surface heating through the topmost cell, the energy released by mixing
415  ! surface cooling & brine rejection down through the topmost cell, and
416  ! homogenizing the shortwave heating within that cell. This sets the energy
417  ! and ustar and wstar available to drive mixing at the first interior
418  ! interface.
419  do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then
420 
421  ! Copy the thicknesses and other fields to 1-d arrays.
422  do k=1,nz
423  h(k) = h_2d(i,k) + gv%H_subroundoff ; u(k) = u_2d(i,k) ; v(k) = v_2d(i,k)
424  t0(k) = t_2d(i,k) ; s0(k) = s_2d(i,k) ; tke_forcing(k) = tke_forced_2d(i,k)
425  dsv_dt_1d(k) = dsv_dt_2d(i,k) ; dsv_ds_1d(k) = dsv_ds_2d(i,k)
426  enddo
427  do k=1,nz+1 ; kd(k) = 0.0 ; enddo
428 
429  ! Make local copies of surface forcing and process them.
430  u_star = fluxes%ustar(i,j)
431  u_star_mean = fluxes%ustar_gustless(i,j)
432  b_flux = buoy_flux(i,j)
433  if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then
434  if (fluxes%frac_shelf_h(i,j) > 0.0) &
435  u_star = (1.0 - fluxes%frac_shelf_h(i,j)) * u_star + &
436  fluxes%frac_shelf_h(i,j) * fluxes%ustar_shelf(i,j)
437  endif
438  if (u_star < cs%ustar_min) u_star = cs%ustar_min
439  if (cs%omega_frac >= 1.0) then
440  absf = 2.0*cs%omega
441  else
442  absf = 0.25*((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
443  (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i-1,j))))
444  if (cs%omega_frac > 0.0) &
445  absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
446  endif
447 
448  ! Perhaps provide a first guess for MLD based on a stored previous value.
449  mld_io = -1.0
450  if (cs%MLD_iteration_guess .and. (cs%ML_Depth(i,j) > 0.0)) mld_io = cs%ML_Depth(i,j)
451 
452  call epbl_column(h, u, v, t0, s0, dsv_dt_1d, dsv_ds_1d, tke_forcing, b_flux, absf, &
453  u_star, u_star_mean, dt, mld_io, kd, mixvel, mixlen, gv, &
454  us, cs, ecd, dt_diag=dt_diag, waves=waves, g=g, i=i, j=j)
455 
456 
457  ! Copy the diffusivities to a 2-d array.
458  do k=1,nz+1
459  kd_2d(i,k) = kd(k)
460  enddo
461  cs%ML_depth(i,j) = mld_io
462 
463  if (present(dt_expected)) then
464  do k=1,nz ; dt_expected(i,j,k) = ecd%dT_expect(k) ; enddo
465  endif
466  if (present(ds_expected)) then
467  do k=1,nz ; ds_expected(i,j,k) = ecd%dS_expect(k) ; enddo
468  endif
469 
470  if (cs%TKE_diagnostics) then
471  cs%diag_TKE_MKE(i,j) = cs%diag_TKE_MKE(i,j) + ecd%dTKE_MKE
472  cs%diag_TKE_conv(i,j) = cs%diag_TKE_conv(i,j) + ecd%dTKE_conv
473  cs%diag_TKE_forcing(i,j) = cs%diag_TKE_forcing(i,j) + ecd%dTKE_forcing
474  cs%diag_TKE_wind(i,j) = cs%diag_TKE_wind(i,j) + ecd%dTKE_wind
475  cs%diag_TKE_mixing(i,j) = cs%diag_TKE_mixing(i,j) + ecd%dTKE_mixing
476  cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + ecd%dTKE_mech_decay
477  cs%diag_TKE_conv_decay(i,j) = cs%diag_TKE_conv_decay(i,j) + ecd%dTKE_conv_decay
478  ! CS%diag_TKE_unbalanced(i,j) = CS%diag_TKE_unbalanced(i,j) + eCD%dTKE_unbalanced
479  endif
480  ! Write to 3-D for outputing Mixing length and velocity scale.
481  if (cs%id_Mixing_Length>0) then ; do k=1,nz
482  cs%Mixing_Length(i,j,k) = mixlen(k)
483  enddo ; endif
484  if (cs%id_Velocity_Scale>0) then ; do k=1,nz
485  cs%Velocity_Scale(i,j,k) = mixvel(k)
486  enddo ; endif
487  if (allocated(cs%mstar_mix)) cs%mstar_mix(i,j) = ecd%mstar
488  if (allocated(cs%mstar_lt)) cs%mstar_lt(i,j) = ecd%mstar_LT
489  if (allocated(cs%La)) cs%La(i,j) = ecd%LA
490  if (allocated(cs%La_mod)) cs%La_mod(i,j) = ecd%LAmod
491  else ! End of the ocean-point part of the i-loop
492  ! For masked points, Kd_int must still be set (to 0) because it has intent out.
493  do k=1,nz+1 ; kd_2d(i,k) = 0. ; enddo
494  cs%ML_depth(i,j) = 0.0
495 
496  if (present(dt_expected)) then
497  do k=1,nz ; dt_expected(i,j,k) = 0.0 ; enddo
498  endif
499  if (present(ds_expected)) then
500  do k=1,nz ; ds_expected(i,j,k) = 0.0 ; enddo
501  endif
502  endif ; enddo ! Close of i-loop - Note unusual loop order!
503 
504  do k=1,nz+1 ; do i=is,ie ; kd_int(i,j,k) = kd_2d(i,k) ; enddo ; enddo
505 
506  enddo ! j-loop
507 
508  if (write_diags) then
509  if (cs%id_ML_depth > 0) call post_data(cs%id_ML_depth, cs%ML_depth, cs%diag)
510  if (cs%id_TKE_wind > 0) call post_data(cs%id_TKE_wind, cs%diag_TKE_wind, cs%diag)
511  if (cs%id_TKE_MKE > 0) call post_data(cs%id_TKE_MKE, cs%diag_TKE_MKE, cs%diag)
512  if (cs%id_TKE_conv > 0) call post_data(cs%id_TKE_conv, cs%diag_TKE_conv, cs%diag)
513  if (cs%id_TKE_forcing > 0) call post_data(cs%id_TKE_forcing, cs%diag_TKE_forcing, cs%diag)
514  if (cs%id_TKE_mixing > 0) call post_data(cs%id_TKE_mixing, cs%diag_TKE_mixing, cs%diag)
515  if (cs%id_TKE_mech_decay > 0) &
516  call post_data(cs%id_TKE_mech_decay, cs%diag_TKE_mech_decay, cs%diag)
517  if (cs%id_TKE_conv_decay > 0) &
518  call post_data(cs%id_TKE_conv_decay, cs%diag_TKE_conv_decay, cs%diag)
519  if (cs%id_Mixing_Length > 0) call post_data(cs%id_Mixing_Length, cs%Mixing_Length, cs%diag)
520  if (cs%id_Velocity_Scale >0) call post_data(cs%id_Velocity_Scale, cs%Velocity_Scale, cs%diag)
521  if (cs%id_MSTAR_MIX > 0) call post_data(cs%id_MSTAR_MIX, cs%MSTAR_MIX, cs%diag)
522  if (cs%id_LA > 0) call post_data(cs%id_LA, cs%LA, cs%diag)
523  if (cs%id_LA_MOD > 0) call post_data(cs%id_LA_MOD, cs%LA_MOD, cs%diag)
524  if (cs%id_MSTAR_LT > 0) call post_data(cs%id_MSTAR_LT, cs%MSTAR_LT, cs%diag)
525  endif
526 
527  if (debug) deallocate(ecd%dT_expect, ecd%dS_expect)
528 

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

Here is the call graph for this function:

◆ energetic_pbl_end()

subroutine, public mom_energetic_pbl::energetic_pbl_end ( type(energetic_pbl_cs), pointer  CS)

Clean up and deallocate memory associated with the energetic_PBL module.

Parameters
csEnergetic_PBL control structure that will be deallocated in this subroutine.

Definition at line 2379 of file MOM_energetic_PBL.F90.

2379  type(energetic_PBL_CS), pointer :: CS !< Energetic_PBL control structure that
2380  !! will be deallocated in this subroutine.
2381 
2382  if (.not.associated(cs)) return
2383 
2384  if (allocated(cs%ML_depth)) deallocate(cs%ML_depth)
2385  if (allocated(cs%LA)) deallocate(cs%LA)
2386  if (allocated(cs%LA_MOD)) deallocate(cs%LA_MOD)
2387  if (allocated(cs%MSTAR_MIX)) deallocate(cs%MSTAR_MIX)
2388  if (allocated(cs%MSTAR_LT)) deallocate(cs%MSTAR_LT)
2389  if (allocated(cs%diag_TKE_wind)) deallocate(cs%diag_TKE_wind)
2390  if (allocated(cs%diag_TKE_MKE)) deallocate(cs%diag_TKE_MKE)
2391  if (allocated(cs%diag_TKE_conv)) deallocate(cs%diag_TKE_conv)
2392  if (allocated(cs%diag_TKE_forcing)) deallocate(cs%diag_TKE_forcing)
2393  if (allocated(cs%diag_TKE_mixing)) deallocate(cs%diag_TKE_mixing)
2394  if (allocated(cs%diag_TKE_mech_decay)) deallocate(cs%diag_TKE_mech_decay)
2395  if (allocated(cs%diag_TKE_conv_decay)) deallocate(cs%diag_TKE_conv_decay)
2396  if (allocated(cs%Mixing_Length)) deallocate(cs%Mixing_Length)
2397  if (allocated(cs%Velocity_Scale)) deallocate(cs%Velocity_Scale)
2398 
2399  deallocate(cs)
2400 

◆ energetic_pbl_get_mld()

subroutine, public mom_energetic_pbl::energetic_pbl_get_mld ( type(energetic_pbl_cs), pointer  CS,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(out)  MLD,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
real, intent(in), optional  m_to_MLD_units 
)

Copies the ePBL active mixed layer depth into MLD.

Parameters
csControl structure for ePBL
[in]gGrid structure
[in]usA dimensional unit scaling type
[out]mldDepth of ePBL active mixing layer [m or other units]
[in]m_to_mld_unitsA conversion factor to the desired units for MLD

Definition at line 1920 of file MOM_energetic_PBL.F90.

1920  type(energetic_PBL_CS), pointer :: CS !< Control structure for ePBL
1921  type(ocean_grid_type), intent(in) :: G !< Grid structure
1922  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1923  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: MLD !< Depth of ePBL active mixing layer [m or other units]
1924  real, optional, intent(in) :: m_to_MLD_units !< A conversion factor to the
1925  !! desired units for MLD
1926  ! Local variables
1927  real :: scale ! A dimensional rescaling factor
1928  integer :: i,j
1929 
1930  scale = us%Z_to_m ; if (present(m_to_mld_units)) scale = scale * m_to_mld_units
1931 
1932  do j=g%jsc,g%jec ; do i=g%isc,g%iec
1933  mld(i,j) = scale*cs%ML_Depth(i,j)
1934  enddo ; enddo
1935 

Referenced by mom_diabatic_driver::diabatic_ale(), mom_diabatic_driver::diabatic_ale_legacy(), mom_lateral_boundary_diffusion::lateral_boundary_diffusion(), and mom_neutral_diffusion::neutral_diffusion_calc_coeffs().

Here is the caller graph for this function:

◆ energetic_pbl_init()

subroutine, public mom_energetic_pbl::energetic_pbl_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(energetic_pbl_cs), pointer  CS 
)

This subroutine initializes the energetic_PBL module.

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

Definition at line 1941 of file MOM_energetic_PBL.F90.

1941  type(time_type), target, intent(in) :: Time !< The current model time
1942  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1943  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1944  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1945  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1946  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic output
1947  type(energetic_PBL_CS), pointer :: CS !< A pointer that is set to point to the control
1948  !! structure for this module
1949  ! Local variables
1950  ! This include declares and sets the variable "version".
1951 # include "version_variable.h"
1952  character(len=40) :: mdl = "MOM_energetic_PBL" ! This module's name.
1953  character(len=20) :: tmpstr
1954  real :: omega_frac_dflt
1955  real :: R_Z3_T3_to_kg_s3 ! A conversion factor for work diagnostics [kg T3 R-1 Z-3 s-3 ~> nondim]
1956  integer :: isd, ied, jsd, jed
1957  integer :: mstar_mode, LT_enhance, wT_mode
1958  logical :: default_2018_answers
1959  logical :: use_temperature, use_omega
1960  logical :: use_la_windsea
1961  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1962 
1963  if (associated(cs)) then
1964  call mom_error(warning, "mixedlayer_init called with an associated"//&
1965  "associated control structure.")
1966  return
1967  else ; allocate(cs) ; endif
1968 
1969  cs%diag => diag
1970  cs%Time => time
1971 
1972 ! Set default, read and log parameters
1973  call log_version(param_file, mdl, version, "")
1974 
1975 
1976 !/1. General ePBL settings
1977  call get_param(param_file, mdl, "OMEGA", cs%omega, &
1978  "The rotation rate of the earth.", units="s-1", &
1979  default=7.2921e-5, scale=us%T_to_S)
1980  call get_param(param_file, mdl, "ML_USE_OMEGA", use_omega, &
1981  "If true, use the absolute rotation rate instead of the "//&
1982  "vertical component of rotation when setting the decay "//&
1983  "scale for turbulence.", default=.false., do_not_log=.true.)
1984  omega_frac_dflt = 0.0
1985  if (use_omega) then
1986  call mom_error(warning, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
1987  omega_frac_dflt = 1.0
1988  endif
1989  call get_param(param_file, mdl, "ML_OMEGA_FRAC", cs%omega_frac, &
1990  "When setting the decay scale for turbulence, use this "//&
1991  "fraction of the absolute rotation rate blended with the "//&
1992  "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
1993  units="nondim", default=omega_frac_dflt)
1994  call get_param(param_file, mdl, "EKMAN_SCALE_COEF", cs%Ekman_scale_coef, &
1995  "A nondimensional scaling factor controlling the inhibition "//&
1996  "of the diffusive length scale by rotation. Making this larger "//&
1997  "decreases the PBL diffusivity.", units="nondim", default=1.0)
1998  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
1999  "This sets the default value for the various _2018_ANSWERS parameters.", &
2000  default=.true.)
2001  call get_param(param_file, mdl, "EPBL_2018_ANSWERS", cs%answers_2018, &
2002  "If true, use the order of arithmetic and expressions that recover the "//&
2003  "answers from the end of 2018. Otherwise, use updated and more robust "//&
2004  "forms of the same expressions.", default=default_2018_answers)
2005 
2006 
2007  call get_param(param_file, mdl, "EPBL_ORIGINAL_PE_CALC", cs%orig_PE_calc, &
2008  "If true, the ePBL code uses the original form of the "//&
2009  "potential energy change code. Otherwise, the newer "//&
2010  "version that can work with successive increments to the "//&
2011  "diffusivity in upward or downward passes is used.", default=.true.)
2012 
2013  call get_param(param_file, mdl, "MKE_TO_TKE_EFFIC", cs%MKE_to_TKE_effic, &
2014  "The efficiency with which mean kinetic energy released "//&
2015  "by mechanically forced entrainment of the mixed layer "//&
2016  "is converted to turbulent kinetic energy.", units="nondim", &
2017  default=0.0)
2018  call get_param(param_file, mdl, "TKE_DECAY", cs%TKE_decay, &
2019  "TKE_DECAY relates the vertical rate of decay of the "//&
2020  "TKE available for mechanical entrainment to the natural "//&
2021  "Ekman depth.", units="nondim", default=2.5)
2022 
2023 
2024 !/2. Options related to setting MSTAR
2025  call get_param(param_file, mdl, "EPBL_MSTAR_SCHEME", tmpstr, &
2026  "EPBL_MSTAR_SCHEME selects the method for setting mstar. Valid values are: \n"//&
2027  "\t CONSTANT - Use a fixed mstar given by MSTAR \n"//&
2028  "\t OM4 - Use L_Ekman/L_Obukhov in the sabilizing limit, as in OM4 \n"//&
2029  "\t REICHL_H18 - Use the scheme documented in Reichl & Hallberg, 2018.", &
2030  default=constant_string, do_not_log=.true.)
2031  call get_param(param_file, mdl, "MSTAR_MODE", mstar_mode, default=-1)
2032  if (mstar_mode == 0) then
2033  tmpstr = constant_string
2034  call mom_error(warning, "Use EPBL_MSTAR_SCHEME = CONSTANT instead of the archaic MSTAR_MODE = 0.")
2035  elseif (mstar_mode == 1) then
2036  call mom_error(fatal, "You are using a legacy mstar mode in ePBL that has been phased out. "//&
2037  "If you need to use this setting please report this error. Also use "//&
2038  "EPBL_MSTAR_SCHEME to specify the scheme for mstar.")
2039  elseif (mstar_mode == 2) then
2040  tmpstr = om4_string
2041  call mom_error(warning, "Use EPBL_MSTAR_SCHEME = OM4 instead of the archaic MSTAR_MODE = 2.")
2042  elseif (mstar_mode == 3) then
2043  tmpstr = rh18_string
2044  call mom_error(warning, "Use EPBL_MSTAR_SCHEME = REICHL_H18 instead of the archaic MSTAR_MODE = 3.")
2045  elseif (mstar_mode > 3) then
2046  call mom_error(fatal, "An unrecognized value of the obsolete parameter MSTAR_MODE was specified.")
2047  endif
2048  call log_param(param_file, mdl, "EPBL_MSTAR_SCHEME", tmpstr, &
2049  "EPBL_MSTAR_SCHEME selects the method for setting mstar. Valid values are: \n"//&
2050  "\t CONSTANT - Use a fixed mstar given by MSTAR \n"//&
2051  "\t OM4 - Use L_Ekman/L_Obukhov in the sabilizing limit, as in OM4 \n"//&
2052  "\t REICHL_H18 - Use the scheme documented in Reichl & Hallberg, 2018.", &
2053  default=constant_string)
2054  tmpstr = uppercase(tmpstr)
2055  select case (tmpstr)
2056  case (constant_string)
2057  cs%mstar_Scheme = use_fixed_mstar
2058  case (om4_string)
2059  cs%mstar_Scheme = mstar_from_ekman
2060  case (rh18_string)
2061  cs%mstar_Scheme = mstar_from_rh18
2062  case default
2063  call mom_mesg('energetic_PBL_init: EPBL_MSTAR_SCHEME ="'//trim(tmpstr)//'"', 0)
2064  call mom_error(fatal, "energetic_PBL_init: Unrecognized setting "// &
2065  "EPBL_MSTAR_SCHEME = "//trim(tmpstr)//" found in input file.")
2066  end select
2067 
2068  call get_param(param_file, mdl, "MSTAR", cs%fixed_mstar, &
2069  "The ratio of the friction velocity cubed to the TKE input to the "//&
2070  "mixed layer. This option is used if EPBL_MSTAR_SCHEME = CONSTANT.", &
2071  units="nondim", default=1.2, do_not_log=(cs%mstar_scheme/=use_fixed_mstar))
2072  call get_param(param_file, mdl, "MSTAR_CAP", cs%mstar_cap, &
2073  "If this value is positive, it sets the maximum value of mstar "//&
2074  "allowed in ePBL. (This is not used if EPBL_MSTAR_SCHEME = CONSTANT).", &
2075  units="nondim", default=-1.0, do_not_log=(cs%mstar_scheme==use_fixed_mstar))
2076  ! mstar_scheme==MStar_from_Ekman options
2077  call get_param(param_file, mdl, "MSTAR2_COEF1", cs%MSTAR_COEF, &
2078  "Coefficient in computing mstar when rotation and stabilizing "//&
2079  "effects are both important (used if EPBL_MSTAR_SCHEME = OM4).", &
2080  units="nondim", default=0.3, do_not_log=(cs%mstar_scheme/=mstar_from_ekman))
2081  call get_param(param_file, mdl, "MSTAR2_COEF2", cs%C_EK, &
2082  "Coefficient in computing mstar when only rotation limits "// &
2083  "the total mixing (used if EPBL_MSTAR_SCHEME = OM4)", &
2084  units="nondim", default=0.085, do_not_log=(cs%mstar_scheme/=mstar_from_ekman))
2085  ! mstar_scheme==MStar_from_RH18 options
2086  call get_param(param_file, mdl, "RH18_MSTAR_CN1", cs%RH18_mstar_cn1,&
2087  "MSTAR_N coefficient 1 (outter-most coefficient for fit). "//&
2088  "The value of 0.275 is given in RH18. Increasing this "//&
2089  "coefficient increases MSTAR for all values of Hf/ust, but more "//&
2090  "effectively at low values (weakly developed OSBLs).", &
2091  units="nondim", default=0.275, do_not_log=(cs%mstar_scheme/=mstar_from_rh18))
2092  call get_param(param_file, mdl, "RH18_MSTAR_CN2", cs%RH18_mstar_cn2,&
2093  "MSTAR_N coefficient 2 (coefficient outside of exponential decay). "//&
2094  "The value of 8.0 is given in RH18. Increasing this coefficient "//&
2095  "increases MSTAR for all values of HF/ust, with a much more even "//&
2096  "effect across a wide range of Hf/ust than CN1.", &
2097  units="nondim", default=8.0, do_not_log=(cs%mstar_scheme/=mstar_from_rh18))
2098  call get_param(param_file, mdl, "RH18_MSTAR_CN3", cs%RH18_mstar_CN3,&
2099  "MSTAR_N coefficient 3 (exponential decay coefficient). "//&
2100  "The value of -5.0 is given in RH18. Increasing this increases how "//&
2101  "quickly the value of MSTAR decreases as Hf/ust increases.", &
2102  units="nondim", default=-5.0, do_not_log=(cs%mstar_scheme/=mstar_from_rh18))
2103  call get_param(param_file, mdl, "RH18_MSTAR_CS1", cs%RH18_mstar_cs1,&
2104  "MSTAR_S coefficient for RH18 in stabilizing limit. "//&
2105  "The value of 0.2 is given in RH18 and increasing it increases "//&
2106  "MSTAR in the presence of a stabilizing surface buoyancy flux.", &
2107  units="nondim", default=0.2, do_not_log=(cs%mstar_scheme/=mstar_from_rh18))
2108  call get_param(param_file, mdl, "RH18_MSTAR_CS2", cs%RH18_mstar_cs2,&
2109  "MSTAR_S exponent for RH18 in stabilizing limit. "//&
2110  "The value of 0.4 is given in RH18 and increasing it increases MSTAR "//&
2111  "exponentially in the presence of a stabilizing surface buoyancy flux.", &
2112  units="nondim", default=0.4, do_not_log=(cs%mstar_scheme/=mstar_from_rh18))
2113 
2114 
2115 !/ Convective turbulence related options
2116  call get_param(param_file, mdl, "NSTAR", cs%nstar, &
2117  "The portion of the buoyant potential energy imparted by "//&
2118  "surface fluxes that is available to drive entrainment "//&
2119  "at the base of mixed layer when that energy is positive.", &
2120  units="nondim", default=0.2)
2121  call get_param(param_file, mdl, "MSTAR_CONV_ADJ", cs%mstar_convect_coef, &
2122  "Coefficient used for reducing mstar during convection "//&
2123  "due to reduction of stable density gradient.", &
2124  units="nondim", default=0.0)
2125 
2126 !/ Mixing Length Options
2127  !### THIS DEFAULT SHOULD BECOME TRUE.
2128  call get_param(param_file, mdl, "USE_MLD_ITERATION", cs%Use_MLD_iteration, &
2129  "A logical that specifies whether or not to use the "//&
2130  "distance to the bottom of the actively turbulent boundary "//&
2131  "layer to help set the EPBL length scale.", default=.false.)
2132  call get_param(param_file, mdl, "EPBL_TRANSITION_SCALE", cs%transLay_scale, &
2133  "A scale for the mixing length in the transition layer "//&
2134  "at the edge of the boundary layer as a fraction of the "//&
2135  "boundary layer thickness.", units="nondim", default=0.1)
2136  if ( cs%Use_MLD_iteration .and. abs(cs%transLay_scale-0.5) >= 0.5) then
2137  call mom_error(fatal, "If flag USE_MLD_ITERATION is true, then "//&
2138  "EPBL_TRANSITION should be greater than 0 and less than 1.")
2139  endif
2140 
2141  call get_param(param_file, mdl, "MLD_ITERATION_GUESS", cs%MLD_ITERATION_GUESS, &
2142  "A logical that specifies whether or not to use the "//&
2143  "previous timestep MLD as a first guess in the MLD iteration. "//&
2144  "The default is false to facilitate reproducibility.", default=.false.)
2145  call get_param(param_file, mdl, "EPBL_MLD_TOLERANCE", cs%MLD_tol, &
2146  "The tolerance for the iteratively determined mixed "//&
2147  "layer depth. This is only used with USE_MLD_ITERATION.", &
2148  units="meter", default=1.0, scale=us%m_to_Z)
2149  call get_param(param_file, mdl, "EPBL_MLD_MAX_ITS", cs%max_MLD_its, &
2150  "The maximum number of iterations that can be used to find a self-consistent "//&
2151  "mixed layer depth. For now, due to the use of bisection, the maximum number "//&
2152  "iteractions needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.", &
2153  default=20, do_not_log=.not.cs%Use_MLD_iteration)
2154  if (.not.cs%Use_MLD_iteration) cs%Max_MLD_Its = 1
2155  call get_param(param_file, mdl, "EPBL_MIN_MIX_LEN", cs%min_mix_len, &
2156  "The minimum mixing length scale that will be used "//&
2157  "by ePBL. The default (0) does not set a minimum.", &
2158  units="meter", default=0.0, scale=us%m_to_Z)
2159 
2160  call get_param(param_file, mdl, "MIX_LEN_EXPONENT", cs%MixLenExponent, &
2161  "The exponent applied to the ratio of the distance to the MLD "//&
2162  "and the MLD depth which determines the shape of the mixing length. "//&
2163  "This is only used if USE_MLD_ITERATION is True.", &
2164  units="nondim", default=2.0)
2165 
2166 !/ Turbulent velocity scale in mixing coefficient
2167  call get_param(param_file, mdl, "EPBL_VEL_SCALE_SCHEME", tmpstr, &
2168  "Selects the method for translating TKE into turbulent velocities. "//&
2169  "Valid values are: \n"//&
2170  "\t CUBE_ROOT_TKE - A constant times the cube root of remaining TKE. \n"//&
2171  "\t REICHL_H18 - Use the scheme based on a combination of w* and v* as \n"//&
2172  "\t documented in Reichl & Hallberg, 2018.", &
2173  default=root_tke_string, do_not_log=.true.)
2174  call get_param(param_file, mdl, "EPBL_VEL_SCALE_MODE", wt_mode, default=-1)
2175  if (wt_mode == 0) then
2176  tmpstr = root_tke_string
2177  call mom_error(warning, "Use EPBL_VEL_SCALE_SCHEME = CUBE_ROOT_TKE instead of the archaic EPBL_VEL_SCALE_MODE = 0.")
2178  elseif (wt_mode == 1) then
2179  tmpstr = rh18_string
2180  call mom_error(warning, "Use EPBL_VEL_SCALE_SCHEME = REICHL_H18 instead of the archaic EPBL_VEL_SCALE_MODE = 1.")
2181  elseif (wt_mode >= 2) then
2182  call mom_error(fatal, "An unrecognized value of the obsolete parameter EPBL_VEL_SCALE_MODE was specified.")
2183  endif
2184  call log_param(param_file, mdl, "EPBL_VEL_SCALE_SCHEME", tmpstr, &
2185  "Selects the method for translating TKE into turbulent velocities. "//&
2186  "Valid values are: \n"//&
2187  "\t CUBE_ROOT_TKE - A constant times the cube root of remaining TKE. \n"//&
2188  "\t REICHL_H18 - Use the scheme based on a combination of w* and v* as \n"//&
2189  "\t documented in Reichl & Hallberg, 2018.", &
2190  default=root_tke_string)
2191  tmpstr = uppercase(tmpstr)
2192  select case (tmpstr)
2193  case (root_tke_string)
2194  cs%wT_scheme = wt_from_croot_tke
2195  case (rh18_string)
2196  cs%wT_scheme = wt_from_rh18
2197  case default
2198  call mom_mesg('energetic_PBL_init: EPBL_VEL_SCALE_SCHEME ="'//trim(tmpstr)//'"', 0)
2199  call mom_error(fatal, "energetic_PBL_init: Unrecognized setting "// &
2200  "EPBL_VEL_SCALE_SCHEME = "//trim(tmpstr)//" found in input file.")
2201  end select
2202 
2203  call get_param(param_file, mdl, "WSTAR_USTAR_COEF", cs%wstar_ustar_coef, &
2204  "A ratio relating the efficiency with which convectively "//&
2205  "released energy is converted to a turbulent velocity, "//&
2206  "relative to mechanically forced TKE. Making this larger "//&
2207  "increases the BL diffusivity", units="nondim", default=1.0)
2208  call get_param(param_file, mdl, "EPBL_VEL_SCALE_FACTOR", cs%vstar_scale_fac, &
2209  "An overall nondimensional scaling factor for wT. "//&
2210  "Making this larger increases the PBL diffusivity.", &
2211  units="nondim", default=1.0)
2212  call get_param(param_file, mdl, "VSTAR_SURF_FAC", cs%vstar_surf_fac,&
2213  "The proportionality times ustar to set vstar at the surface.", &
2214  units="nondim", default=1.2)
2215 
2216  !/ Options related to Langmuir turbulence
2217  call get_param(param_file, mdl, "USE_LA_LI2016", use_la_windsea, &
2218  "A logical to use the Li et al. 2016 (submitted) formula to "//&
2219  "determine the Langmuir number.", units="nondim", default=.false.)
2220  ! Note this can be activated in other ways, but this preserves the old method.
2221  if (use_la_windsea) then
2222  cs%USE_LT = .true.
2223  else
2224  call get_param(param_file, mdl, "EPBL_LT", cs%USE_LT, &
2225  "A logical to use a LT parameterization.", &
2226  units="nondim", default=.false.)
2227  endif
2228  if (cs%USE_LT) then
2229  call get_param(param_file, mdl, "EPBL_LANGMUIR_SCHEME", tmpstr, &
2230  "EPBL_LANGMUIR_SCHEME selects the method for including Langmuir turbulence. "//&
2231  "Valid values are: \n"//&
2232  "\t NONE - Do not do any extra mixing due to Langmuir turbulence \n"//&
2233  "\t RESCALE - Use a multiplicative rescaling of mstar to account for Langmuir turbulence \n"//&
2234  "\t ADDITIVE - Add a Langmuir turblence contribution to mstar to other contributions", &
2235  default=none_string, do_not_log=.true.)
2236  call get_param(param_file, mdl, "LT_ENHANCE", lt_enhance, default=-1)
2237  if (lt_enhance == 0) then
2238  tmpstr = none_string
2239  call mom_error(warning, "Use EPBL_LANGMUIR_SCHEME = NONE instead of the archaic LT_ENHANCE = 0.")
2240  elseif (lt_enhance == 1) then
2241  call mom_error(fatal, "You are using a legacy LT_ENHANCE mode in ePBL that has been phased out. "//&
2242  "If you need to use this setting please report this error. Also use "//&
2243  "EPBL_LANGMUIR_SCHEME to specify the scheme for mstar.")
2244  elseif (lt_enhance == 2) then
2245  tmpstr = rescaled_string
2246  call mom_error(warning, "Use EPBL_LANGMUIR_SCHEME = RESCALE instead of the archaic LT_ENHANCE = 2.")
2247  elseif (lt_enhance == 3) then
2248  tmpstr = additive_string
2249  call mom_error(warning, "Use EPBL_LANGMUIR_SCHEME = ADDITIVE instead of the archaic LT_ENHANCE = 3.")
2250  elseif (lt_enhance > 3) then
2251  call mom_error(fatal, "An unrecognized value of the obsolete parameter LT_ENHANCE was specified.")
2252  endif
2253  call log_param(param_file, mdl, "EPBL_LANGMUIR_SCHEME", tmpstr, &
2254  "EPBL_LANGMUIR_SCHEME selects the method for including Langmuir turbulence. "//&
2255  "Valid values are: \n"//&
2256  "\t NONE - Do not do any extra mixing due to Langmuir turbulence \n"//&
2257  "\t RESCALE - Use a multiplicative rescaling of mstar to account for Langmuir turbulence \n"//&
2258  "\t ADDITIVE - Add a Langmuir turblence contribution to mstar to other contributions", &
2259  default=none_string)
2260  tmpstr = uppercase(tmpstr)
2261  select case (tmpstr)
2262  case (none_string)
2263  cs%LT_enhance_form = no_langmuir
2264  case (rescaled_string)
2265  cs%LT_enhance_form = langmuir_rescale
2266  case (additive_string)
2267  cs%LT_enhance_form = langmuir_add
2268  case default
2269  call mom_mesg('energetic_PBL_init: EPBL_LANGMUIR_SCHEME ="'//trim(tmpstr)//'"', 0)
2270  call mom_error(fatal, "energetic_PBL_init: Unrecognized setting "// &
2271  "EPBL_LANGMUIR_SCHEME = "//trim(tmpstr)//" found in input file.")
2272  end select
2273 
2274  call get_param(param_file, mdl, "LT_ENHANCE_COEF", cs%LT_ENHANCE_COEF, &
2275  "Coefficient for Langmuir enhancement of mstar", &
2276  units="nondim", default=0.447, do_not_log=(cs%LT_enhance_form==no_langmuir))
2277  call get_param(param_file, mdl, "LT_ENHANCE_EXP", cs%LT_ENHANCE_EXP, &
2278  "Exponent for Langmuir enhancementt of mstar", &
2279  units="nondim", default=-1.33, do_not_log=(cs%LT_enhance_form==no_langmuir))
2280  call get_param(param_file, mdl, "LT_MOD_LAC1", cs%LaC_MLDoEK, &
2281  "Coefficient for modification of Langmuir number due to "//&
2282  "MLD approaching Ekman depth.", &
2283  units="nondim", default=-0.87, do_not_log=(cs%LT_enhance_form==no_langmuir))
2284  call get_param(param_file, mdl, "LT_MOD_LAC2", cs%LaC_MLDoOB_stab, &
2285  "Coefficient for modification of Langmuir number due to "//&
2286  "MLD approaching stable Obukhov depth.", &
2287  units="nondim", default=0.0, do_not_log=(cs%LT_enhance_form==no_langmuir))
2288  call get_param(param_file, mdl, "LT_MOD_LAC3", cs%LaC_MLDoOB_un, &
2289  "Coefficient for modification of Langmuir number due to "//&
2290  "MLD approaching unstable Obukhov depth.", &
2291  units="nondim", default=0.0, do_not_log=(cs%LT_enhance_form==no_langmuir))
2292  call get_param(param_file, mdl, "LT_MOD_LAC4", cs%Lac_EKoOB_stab, &
2293  "Coefficient for modification of Langmuir number due to "//&
2294  "ratio of Ekman to stable Obukhov depth.", &
2295  units="nondim", default=0.95, do_not_log=(cs%LT_enhance_form==no_langmuir))
2296  call get_param(param_file, mdl, "LT_MOD_LAC5", cs%Lac_EKoOB_un, &
2297  "Coefficient for modification of Langmuir number due to "//&
2298  "ratio of Ekman to unstable Obukhov depth.", &
2299  units="nondim", default=0.95, do_not_log=(cs%LT_enhance_form==no_langmuir))
2300  endif
2301 
2302 
2303 !/ Logging parameters
2304  ! This gives a minimum decay scale that is typically much less than Angstrom.
2305  cs%ustar_min = 2e-4*cs%omega*(gv%Angstrom_Z + gv%H_to_Z*gv%H_subroundoff)
2306  call log_param(param_file, mdl, "EPBL_USTAR_MIN", cs%ustar_min*us%Z_to_m*us%s_to_T, &
2307  "The (tiny) minimum friction velocity used within the "//&
2308  "ePBL code, derived from OMEGA and ANGSTROM.", units="m s-1")
2309 
2310 
2311 !/ Checking output flags
2312  r_z3_t3_to_kg_s3 = us%R_to_kg_m3 * us%Z_to_m**3 * us%s_to_T**3
2313  cs%id_ML_depth = register_diag_field('ocean_model', 'ePBL_h_ML', diag%axesT1, &
2314  time, 'Surface boundary layer depth', 'm', conversion=us%Z_to_m, &
2315  cmor_long_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme')
2316  cs%id_TKE_wind = register_diag_field('ocean_model', 'ePBL_TKE_wind', diag%axesT1, &
2317  time, 'Wind-stirring source of mixed layer TKE', 'm3 s-3', conversion=r_z3_t3_to_kg_s3)
2318  cs%id_TKE_MKE = register_diag_field('ocean_model', 'ePBL_TKE_MKE', diag%axesT1, &
2319  time, 'Mean kinetic energy source of mixed layer TKE', 'm3 s-3', conversion=r_z3_t3_to_kg_s3)
2320  cs%id_TKE_conv = register_diag_field('ocean_model', 'ePBL_TKE_conv', diag%axesT1, &
2321  time, 'Convective source of mixed layer TKE', 'm3 s-3', conversion=r_z3_t3_to_kg_s3)
2322  cs%id_TKE_forcing = register_diag_field('ocean_model', 'ePBL_TKE_forcing', diag%axesT1, &
2323  time, 'TKE consumed by mixing surface forcing or penetrative shortwave radation '//&
2324  'through model layers', 'm3 s-3', conversion=r_z3_t3_to_kg_s3)
2325  cs%id_TKE_mixing = register_diag_field('ocean_model', 'ePBL_TKE_mixing', diag%axesT1, &
2326  time, 'TKE consumed by mixing that deepens the mixed layer', 'm3 s-3', conversion=r_z3_t3_to_kg_s3)
2327  cs%id_TKE_mech_decay = register_diag_field('ocean_model', 'ePBL_TKE_mech_decay', diag%axesT1, &
2328  time, 'Mechanical energy decay sink of mixed layer TKE', 'm3 s-3', conversion=r_z3_t3_to_kg_s3)
2329  cs%id_TKE_conv_decay = register_diag_field('ocean_model', 'ePBL_TKE_conv_decay', diag%axesT1, &
2330  time, 'Convective energy decay sink of mixed layer TKE', 'm3 s-3', conversion=r_z3_t3_to_kg_s3)
2331  cs%id_Mixing_Length = register_diag_field('ocean_model', 'Mixing_Length', diag%axesTi, &
2332  time, 'Mixing Length that is used', 'm', conversion=us%Z_to_m)
2333  cs%id_Velocity_Scale = register_diag_field('ocean_model', 'Velocity_Scale', diag%axesTi, &
2334  time, 'Velocity Scale that is used.', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
2335  cs%id_MSTAR_mix = register_diag_field('ocean_model', 'MSTAR', diag%axesT1, &
2336  time, 'Total mstar that is used.', 'nondim')
2337 
2338  if (cs%use_LT) then
2339  cs%id_LA = register_diag_field('ocean_model', 'LA', diag%axesT1, &
2340  time, 'Langmuir number.', 'nondim')
2341  cs%id_LA_mod = register_diag_field('ocean_model', 'LA_MOD', diag%axesT1, &
2342  time, 'Modified Langmuir number.', 'nondim')
2343  cs%id_MSTAR_LT = register_diag_field('ocean_model', 'MSTAR_LT', diag%axesT1, &
2344  time, 'Increase in mstar due to Langmuir Turbulence.', 'nondim')
2345  endif
2346 
2347  call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
2348  "If true, temperature and salinity are used as state "//&
2349  "variables.", default=.true.)
2350 
2351  if (max(cs%id_TKE_wind, cs%id_TKE_MKE, cs%id_TKE_conv, &
2352  cs%id_TKE_mixing, cs%id_TKE_mech_decay, cs%id_TKE_forcing, &
2353  cs%id_TKE_conv_decay) > 0) then
2354  call safe_alloc_alloc(cs%diag_TKE_wind, isd, ied, jsd, jed)
2355  call safe_alloc_alloc(cs%diag_TKE_MKE, isd, ied, jsd, jed)
2356  call safe_alloc_alloc(cs%diag_TKE_conv, isd, ied, jsd, jed)
2357  call safe_alloc_alloc(cs%diag_TKE_forcing, isd, ied, jsd, jed)
2358  call safe_alloc_alloc(cs%diag_TKE_mixing, isd, ied, jsd, jed)
2359  call safe_alloc_alloc(cs%diag_TKE_mech_decay, isd, ied, jsd, jed)
2360  call safe_alloc_alloc(cs%diag_TKE_conv_decay, isd, ied, jsd, jed)
2361 
2362  cs%TKE_diagnostics = .true.
2363  endif
2364  if (cs%id_Velocity_Scale>0) call safe_alloc_alloc(cs%Velocity_Scale, isd, ied, jsd, jed, gv%ke+1)
2365  if (cs%id_Mixing_Length>0) call safe_alloc_alloc(cs%Mixing_Length, isd, ied, jsd, jed, gv%ke+1)
2366 
2367  call safe_alloc_alloc(cs%ML_depth, isd, ied, jsd, jed)
2368  if (max(cs%id_mstar_mix, cs%id_LA, cs%id_LA_mod, cs%id_MSTAR_LT ) >0) then
2369  call safe_alloc_alloc(cs%Mstar_mix, isd, ied, jsd, jed)
2370  call safe_alloc_alloc(cs%LA, isd, ied, jsd, jed)
2371  call safe_alloc_alloc(cs%LA_MOD, isd, ied, jsd, jed)
2372  call safe_alloc_alloc(cs%MSTAR_LT, isd, ied, jsd, jed)
2373  endif
2374 

References additive_string, constant_string, langmuir_add, langmuir_rescale, mom_error_handler::mom_error(), mom_error_handler::mom_mesg(), mstar_from_ekman, mstar_from_rh18, no_langmuir, none_string, om4_string, rescaled_string, rh18_string, root_tke_string, mom_string_functions::uppercase(), use_fixed_mstar, wt_from_croot_tke, and wt_from_rh18.

Here is the call graph for this function:

◆ epbl_column()

subroutine mom_energetic_pbl::epbl_column ( real, dimension(szk_(gv)), intent(in)  h,
real, dimension(szk_(gv)), intent(in)  u,
real, dimension(szk_(gv)), intent(in)  v,
real, dimension(szk_(gv)), intent(in)  T0,
real, dimension(szk_(gv)), intent(in)  S0,
real, dimension(szk_(gv)), intent(in)  dSV_dT,
real, dimension(szk_(gv)), intent(in)  dSV_dS,
real, dimension(szk_(gv)), intent(in)  TKE_forcing,
real, intent(in)  B_flux,
real, intent(in)  absf,
real, intent(in)  u_star,
real, intent(in)  u_star_mean,
real, intent(in)  dt,
real, intent(inout)  MLD_io,
real, dimension(szk_(gv)+1), intent(out)  Kd,
real, dimension(szk_(gv)+1), intent(out)  mixvel,
real, dimension(szk_(gv)+1), intent(out)  mixlen,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(energetic_pbl_cs), pointer  CS,
type(epbl_column_diags), intent(inout)  eCD,
real, intent(in), optional  dt_diag,
type(wave_parameters_cs), optional, pointer  Waves,
type(ocean_grid_type), intent(inout), optional  G,
integer, intent(in), optional  i,
integer, intent(in), optional  j 
)
private

This subroutine determines the diffusivities from the integrated energetics mixed layer model for a single column of water.

Parameters
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]hLayer thicknesses [H ~> m or kg m-2].
[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]t0The initial layer temperatures [degC].
[in]s0The initial layer salinities [ppt].
[in]dsv_dtThe partial derivative of in-situ specific volume with potential temperature [R-1 degC-1 ~> m3 kg-1 degC-1].
[in]dsv_dsThe partial derivative of in-situ specific volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
[in]tke_forcingThe forcing requirements to homogenize the forcing that has been applied to each layer [R Z3 T-2 ~> J m-2].
[in]b_fluxThe surface buoyancy flux [Z2 T-3 ~> m2 s-3]
[in]absfThe absolute value of the Coriolis parameter [T-1 ~> s-1].
[in]u_starThe surface friction velocity [Z T-1 ~> m s-1].
[in]u_star_meanThe surface friction velocity without any contribution from unresolved gustiness [Z T-1 ~> m s-1].
[in,out]mld_ioA first guess at the mixed layer depth on input, and the calculated mixed layer depth on output [Z ~> m].
[in]dtTime increment [T ~> s].
[out]kdThe diagnosed diffusivities at interfaces
[out]mixvelThe mixing velocity scale used in Kd
[out]mixlenThe mixing length scale used in Kd [Z ~> m].
csThe control structure returned by a previous call to mixedlayer_init.
[in,out]ecdA container for passing around diagnostics.
[in]dt_diagThe diagnostic time step, which may be less than dt if there are two calls to mixedlayer [T ~> s].
wavesWave CS for Langmuir turbulence
[in,out]gThe ocean's grid structure.
[in]iThe i-index to work on (used for Waves)
[in]jThe i-index to work on (used for Waves)

Definition at line 538 of file MOM_energetic_PBL.F90.

538  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
539  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
540  real, dimension(SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
541  real, dimension(SZK_(GV)), intent(in) :: u !< Zonal velocities interpolated to h points
542  !! [L T-1 ~> m s-1].
543  real, dimension(SZK_(GV)), intent(in) :: v !< Zonal velocities interpolated to h points
544  !! [L T-1 ~> m s-1].
545  real, dimension(SZK_(GV)), intent(in) :: T0 !< The initial layer temperatures [degC].
546  real, dimension(SZK_(GV)), intent(in) :: S0 !< The initial layer salinities [ppt].
547 
548  real, dimension(SZK_(GV)), intent(in) :: dSV_dT !< The partial derivative of in-situ specific
549  !! volume with potential temperature
550  !! [R-1 degC-1 ~> m3 kg-1 degC-1].
551  real, dimension(SZK_(GV)), intent(in) :: dSV_dS !< The partial derivative of in-situ specific
552  !! volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
553  real, dimension(SZK_(GV)), intent(in) :: TKE_forcing !< The forcing requirements to homogenize the
554  !! forcing that has been applied to each layer
555  !! [R Z3 T-2 ~> J m-2].
556  real, intent(in) :: B_flux !< The surface buoyancy flux [Z2 T-3 ~> m2 s-3]
557  real, intent(in) :: absf !< The absolute value of the Coriolis parameter [T-1 ~> s-1].
558  real, intent(in) :: u_star !< The surface friction velocity [Z T-1 ~> m s-1].
559  real, intent(in) :: u_star_mean !< The surface friction velocity without any
560  !! contribution from unresolved gustiness [Z T-1 ~> m s-1].
561  real, intent(inout) :: MLD_io !< A first guess at the mixed layer depth on input, and
562  !! the calculated mixed layer depth on output [Z ~> m].
563  real, intent(in) :: dt !< Time increment [T ~> s].
564  real, dimension(SZK_(GV)+1), &
565  intent(out) :: Kd !< The diagnosed diffusivities at interfaces
566  !! [Z2 T-1 ~> m2 s-1].
567  real, dimension(SZK_(GV)+1), &
568  intent(out) :: mixvel !< The mixing velocity scale used in Kd
569  !! [Z T-1 ~> m s-1].
570  real, dimension(SZK_(GV)+1), &
571  intent(out) :: mixlen !< The mixing length scale used in Kd [Z ~> m].
572  type(energetic_PBL_CS), pointer :: CS !< The control structure returned by a previous
573  !! call to mixedlayer_init.
574  type(ePBL_column_diags), intent(inout) :: eCD !< A container for passing around diagnostics.
575  real, optional, intent(in) :: dt_diag !< The diagnostic time step, which may be less
576  !! than dt if there are two calls to mixedlayer [T ~> s].
577  type(wave_parameters_CS), &
578  optional, pointer :: Waves !< Wave CS for Langmuir turbulence
579  type(ocean_grid_type), &
580  optional, intent(inout) :: G !< The ocean's grid structure.
581  integer, optional, intent(in) :: i !< The i-index to work on (used for Waves)
582  integer, optional, intent(in) :: j !< The i-index to work on (used for Waves)
583 
584 ! This subroutine determines the diffusivities in a single column from the integrated energetics
585 ! planetary boundary layer (ePBL) model. It assumes that heating, cooling and freshwater fluxes
586 ! have already been applied. All calculations are done implicitly, and there
587 ! is no stability limit on the time step.
588 !
589 ! For each interior interface, first discard the TKE to account for mixing
590 ! of shortwave radiation through the next denser cell. Next drive mixing based
591 ! on the local? values of ustar + wstar, subject to available energy. This
592 ! step sets the value of Kd(K). Any remaining energy is then subject to decay
593 ! before being handed off to the next interface. mech_TKE and conv_PErel are treated
594 ! separately for the purposes of decay, but are used proportionately to drive
595 ! mixing.
596 
597  ! Local variables
598  real, dimension(SZK_(GV)+1) :: &
599  pres_Z, & ! Interface pressures with a rescaling factor to convert interface height
600  ! movements into changes in column potential energy [R Z2 T-2 ~> kg m-1 s-2].
601  hb_hs ! The distance from the bottom over the thickness of the
602  ! water column [nondim].
603  real :: mech_TKE ! The mechanically generated turbulent kinetic energy
604  ! available for mixing over a time step [R Z3 T-2 ~> J m-2].
605  real :: conv_PErel ! The potential energy that has been convectively released
606  ! during this timestep [R Z3 T-2 ~> J m-2]. A portion nstar_FC
607  ! of conv_PErel is available to drive mixing.
608  real :: htot ! The total depth of the layers above an interface [H ~> m or kg m-2].
609  real :: uhtot ! The depth integrated zonal and meridional velocities in the
610  real :: vhtot ! layers above [H L T-1 ~> m2 s-1 or kg m-1 s-1].
611  real :: Idecay_len_TKE ! The inverse of a turbulence decay length scale [H-1 ~> m-1 or m2 kg-1].
612  real :: h_sum ! The total thickness of the water column [H ~> m or kg m-2].
613 
614  real, dimension(SZK_(GV)) :: &
615  dT_to_dColHt, & ! Partial derivative of the total column height with the temperature changes
616  ! within a layer [Z degC-1 ~> m degC-1].
617  ds_to_dcolht, & ! Partial derivative of the total column height with the salinity changes
618  ! within a layer [Z ppt-1 ~> m ppt-1].
619  dt_to_dpe, & ! Partial derivatives of column potential energy with the temperature
620  ! changes within a layer, in [R Z3 T-2 degC-1 ~> J m-2 degC-1].
621  ds_to_dpe, & ! Partial derivatives of column potential energy with the salinity changes
622  ! within a layer, in [R Z3 T-2 ppt-1 ~> J m-2 ppt-1].
623  dt_to_dcolht_a, & ! Partial derivative of the total column height with the temperature changes
624  ! within a layer, including the implicit effects of mixing with layers higher
625  ! in the water column [Z degC-1 ~> m degC-1].
626  ds_to_dcolht_a, & ! Partial derivative of the total column height with the salinity changes
627  ! within a layer, including the implicit effects of mixing with layers higher
628  ! in the water column [Z ppt-1 ~> m ppt-1].
629  dt_to_dpe_a, & ! Partial derivatives of column potential energy with the temperature changes
630  ! within a layer, including the implicit effects of mixing with layers higher
631  ! in the water column [R Z3 T-2 degC-1 ~> J m-2 degC-1].
632  ds_to_dpe_a, & ! Partial derivative of column potential energy with the salinity changes
633  ! within a layer, including the implicit effects of mixing with layers higher
634  ! in the water column [R Z3 T-2 ppt-1 ~> J m-2 ppt-1].
635  c1, & ! c1 is used by the tridiagonal solver [nondim].
636  te, & ! Estimated final values of T in the column [degC].
637  se, & ! Estimated final values of S in the column [ppt].
638  dte, & ! Running (1-way) estimates of temperature change [degC].
639  dse, & ! Running (1-way) estimates of salinity change [ppt].
640  th_a, & ! An effective temperature times a thickness in the layer above, including implicit
641  ! mixing effects with other yet higher layers [degC H ~> degC m or degC kg m-2].
642  sh_a, & ! An effective salinity times a thickness in the layer above, including implicit
643  ! mixing effects with other yet higher layers [ppt H ~> ppt m or ppt kg m-2].
644  th_b, & ! An effective temperature times a thickness in the layer below, including implicit
645  ! mixing effects with other yet lower layers [degC H ~> degC m or degC kg m-2].
646  sh_b ! An effective salinity times a thickness in the layer below, including implicit
647  ! mixing effects with other yet lower layers [ppt H ~> ppt m or ppt kg m-2].
648  real, dimension(SZK_(GV)+1) :: &
649  MixLen_shape, & ! A nondimensional shape factor for the mixing length that
650  ! gives it an appropriate assymptotic value at the bottom of
651  ! the boundary layer.
652  kddt_h ! The diapycnal diffusivity times a timestep divided by the
653  ! average thicknesses around a layer [H ~> m or kg m-2].
654  real :: b1 ! b1 is inverse of the pivot used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
655  real :: hp_a ! An effective pivot thickness of the layer including the effects
656  ! of coupling with layers above [H ~> m or kg m-2]. This is the first term
657  ! in the denominator of b1 in a downward-oriented tridiagonal solver.
658  real :: h_neglect ! A thickness that is so small it is usually lost
659  ! in roundoff and can be neglected [H ~> m or kg m-2].
660  real :: dMass ! The mass per unit area within a layer [Z R ~> kg m-2].
661  real :: dPres ! The hydrostatic pressure change across a layer [R Z2 T-2 ~> kg m-1 s-2 = Pa = J m-3].
662  real :: dMKE_max ! The maximum amount of mean kinetic energy that could be
663  ! converted to turbulent kinetic energy if the velocity in
664  ! the layer below an interface were homogenized with all of
665  ! the water above the interface [R Z3 T-2 ~> J m-2].
666  real :: MKE2_Hharm ! Twice the inverse of the harmonic mean of the thickness
667  ! of a layer and the thickness of the water above, used in
668  ! the MKE conversion equation [H-1 ~> m-1 or m2 kg-1].
669 
670  real :: dt_h ! The timestep divided by the averages of the thicknesses around
671  ! a layer, times a thickness conversion factor [H T m-2 ~> s m-1 or kg s m-4].
672  real :: h_bot ! The distance from the bottom [H ~> m or kg m-2].
673  real :: h_rsum ! The running sum of h from the top [Z ~> m].
674  real :: I_hs ! The inverse of h_sum [H-1 ~> m-1 or m2 kg-1].
675  real :: I_MLD ! The inverse of the current value of MLD [Z-1 ~> m-1].
676  real :: h_tt ! The distance from the surface or up to the next interface
677  ! that did not exhibit turbulent mixing from this scheme plus
678  ! a surface mixing roughness length given by h_tt_min [H ~> m or kg m-2].
679  real :: h_tt_min ! A surface roughness length [H ~> m or kg m-2].
680 
681  real :: C1_3 ! = 1/3.
682  real :: I_dtrho ! 1.0 / (dt * Rho0) times conversion factors in [m3 Z-3 R-1 T2 s-3 ~> m3 kg-1 s-1].
683  ! This is used convert TKE back into ustar^3.
684  real :: vstar ! An in-situ turbulent velocity [Z T-1 ~> m s-1].
685  real :: mstar_total ! The value of mstar used in ePBL [nondim]
686  real :: mstar_LT ! An addition to mstar due to Langmuir turbulence [nondim] (output for diagnostic)
687  real :: MLD_output ! The mixed layer depth output from this routine [Z ~> m].
688  real :: LA ! The value of the Langmuir number [nondim]
689  real :: LAmod ! The modified Langmuir number by convection [nondim]
690  real :: hbs_here ! The local minimum of hb_hs and MixLen_shape, times a
691  ! conversion factor from H to Z [Z H-1 ~> 1 or m3 kg-1].
692  real :: nstar_FC ! The fraction of conv_PErel that can be converted to mixing [nondim].
693  real :: TKE_reduc ! The fraction by which TKE and other energy fields are
694  ! reduced to support mixing [nondim]. between 0 and 1.
695  real :: tot_TKE ! The total TKE available to support mixing at interface K [R Z3 T-2 ~> J m-2].
696  real :: TKE_here ! The total TKE at this point in the algorithm [R Z3 T-2 ~> J m-2].
697  real :: dT_km1_t2 ! A diffusivity-independent term related to the temperature
698  ! change in the layer above the interface [degC].
699  real :: dS_km1_t2 ! A diffusivity-independent term related to the salinity
700  ! change in the layer above the interface [ppt].
701  real :: dTe_term ! A diffusivity-independent term related to the temperature
702  ! change in the layer below the interface [degC H ~> degC m or degC kg m-2].
703  real :: dSe_term ! A diffusivity-independent term related to the salinity
704  ! change in the layer above the interface [ppt H ~> ppt m or ppt kg m-2].
705  real :: dTe_t2 ! A part of dTe_term [degC H ~> degC m or degC kg m-2].
706  real :: dSe_t2 ! A part of dSe_term [ppt H ~> ppt m or ppt kg m-2].
707  real :: dPE_conv ! The convective change in column potential energy [R Z3 T-2 ~> J m-2].
708  real :: MKE_src ! The mean kinetic energy source of TKE due to Kddt_h(K) [R Z3 T-2 ~> J m-2].
709  real :: dMKE_src_dK ! The partial derivative of MKE_src with Kddt_h(K) [R Z3 T-2 H-1 ~> J m-3 or J kg-1].
710  real :: Kd_guess0 ! A first guess of the diapycnal diffusivity [Z2 T-1 ~> m2 s-1].
711  real :: PE_chg_g0 ! The potential energy change when Kd is Kd_guess0 [R Z3 T-2 ~> J m-2]
712  !### The following might be unused.
713  real :: dPEa_dKd_g0 ! The derivative of the change in the potential energy of the column above an interface
714  ! with the diffusivity when the Kd is Kd_guess0 [R Z T-1 ~> J s m-4]
715  real :: Kddt_h_g0 ! The first guess diapycnal diffusivity times a timestep divided
716  ! by the average thicknesses around a layer [H ~> m or kg m-2].
717  real :: PE_chg_max ! The maximum PE change for very large values of Kddt_h(K) [R Z3 T-2 ~> J m-2].
718  real :: dPEc_dKd_Kd0 ! The partial derivative of PE change with Kddt_h(K)
719  ! for very small values of Kddt_h(K) [R Z3 T-2 H-1 ~> J m-3 or J kg-1].
720  real :: PE_chg ! The change in potential energy due to mixing at an
721  ! interface [R Z3 T-2 ~> J m-2], positive for the column increasing
722  ! in potential energy (i.e., consuming TKE).
723  real :: TKE_left ! The amount of turbulent kinetic energy left for the most
724  ! recent guess at Kddt_h(K) [R Z3 T-2 ~> J m-2].
725  real :: dPEc_dKd ! The partial derivative of PE_chg with Kddt_h(K) [R Z3 T-2 H-1 ~> J m-3 or J kg-1].
726  real :: TKE_left_min, TKE_left_max ! Maximum and minimum values of TKE_left [R Z3 T-2 ~> J m-2].
727  real :: Kddt_h_max, Kddt_h_min ! Maximum and minimum values of Kddt_h(K) [H ~> m or kg m-2].
728  real :: Kddt_h_guess ! A guess at the value of Kddt_h(K) [H ~> m or kg m-2].
729  real :: Kddt_h_next ! The next guess at the value of Kddt_h(K) [H ~> m or kg m-2].
730  real :: dKddt_h ! The change between guesses at Kddt_h(K) [H ~> m or kg m-2].
731  real :: dKddt_h_Newt ! The change between guesses at Kddt_h(K) with Newton's method [H ~> m or kg m-2].
732  real :: Kddt_h_newt ! The Newton's method next guess for Kddt_h(K) [H ~> m or kg m-2].
733  real :: exp_kh ! The nondimensional decay of TKE across a layer [nondim].
734  real :: vstar_unit_scale ! A unit converion factor for turbulent velocities [Z T-1 s m-1 ~> 1]
735  logical :: use_Newt ! Use Newton's method for the next guess at Kddt_h(K).
736  logical :: convectively_stable ! If true the water column is convectively stable at this interface.
737  logical :: sfc_connected ! If true the ocean is actively turbulent from the present
738  ! interface all the way up to the surface.
739  logical :: sfc_disconnect ! If true, any turbulence has become disconnected
740  ! from the surface.
741 
742 ! The following are only used for diagnostics.
743  real :: dt__diag ! A copy of dt_diag (if present) or dt [T ~> s].
744  real :: I_dtdiag ! = 1.0 / dt__diag [T-1 ~> s-1].
745 
746  !----------------------------------------------------------------------
747  !/BGR added Aug24,2016 for adding iteration to get boundary layer depth
748  ! - needed to compute new mixing length.
749  real :: MLD_guess, MLD_found ! Mixing Layer depth guessed/found for iteration [Z ~> m].
750  real :: min_MLD ! Iteration bounds [Z ~> m], which are adjusted at each step
751  real :: max_MLD ! - These are initialized based on surface/bottom
752  ! 1. The iteration guesses a value (possibly from prev step or neighbor).
753  ! 2. The iteration checks if value is converged, too shallow, or too deep.
754  ! 3. Based on result adjusts the Max/Min and searches through the water column.
755  ! - If using an accurate guess the iteration is very quick (e.g. if MLD doesn't
756  ! change over timestep). Otherwise it takes 5-10 passes, but has a high
757  ! convergence rate. Other iteration may be tried, but this method seems to
758  ! fail very rarely and the added cost is likely not significant.
759  ! Additionally, when it fails to converge it does so in a reasonable
760  ! manner giving a usable guess. When it does fail, it is due to convection
761  ! within the boundary layer. Likely, a new method e.g. surface_disconnect,
762  ! can improve this.
763  logical :: FIRST_OBL ! Flag for computing "found" Mixing layer depth
764  logical :: OBL_converged ! Flag for convergence of MLD
765  integer :: OBL_it ! Iteration counter
766 
767  real :: Surface_Scale ! Surface decay scale for vstar
768 
769  logical :: debug=.false. ! Change this hard-coded value for debugging.
770 
771  ! The following arrays are used only for debugging purposes.
772  real :: dPE_debug, mixing_debug, taux2, tauy2
773  real, dimension(20) :: TKE_left_itt, PE_chg_itt, Kddt_h_itt, dPEa_dKd_itt, MKE_src_itt
774  real, dimension(SZK_(GV)) :: mech_TKE_k, conv_PErel_k, nstar_k
775  integer, dimension(SZK_(GV)) :: num_itts
776 
777  integer :: k, nz, itt, max_itt
778 
779  nz = gv%ke
780 
781  if (.not. associated(cs)) call mom_error(fatal, "energetic_PBL: "//&
782  "Module must be initialized before it is used.")
783 
784  debug = .false. ; if (allocated(ecd%dT_expect) .or. allocated(ecd%dS_expect)) debug = .true.
785 
786  h_neglect = gv%H_subroundoff
787 
788  c1_3 = 1.0 / 3.0
789  dt__diag = dt ; if (present(dt_diag)) dt__diag = dt_diag
790  i_dtdiag = 1.0 / dt__diag
791  max_itt = 20
792 
793  h_tt_min = 0.0
794  i_dtrho = 0.0 ; if (dt*gv%Rho0 > 0.0) i_dtrho = (us%Z_to_m**3*us%s_to_T**3) / (dt*gv%Rho0)
795  vstar_unit_scale = us%m_to_Z * us%T_to_s
796 
797  mld_guess = mld_io
798 
799 ! Determine the initial mech_TKE and conv_PErel, including the energy required
800 ! to mix surface heating through the topmost cell, the energy released by mixing
801 ! surface cooling & brine rejection down through the topmost cell, and
802 ! homogenizing the shortwave heating within that cell. This sets the energy
803 ! and ustar and wstar available to drive mixing at the first interior
804 ! interface.
805 
806  do k=1,nz+1 ; kd(k) = 0.0 ; enddo
807 
808  pres_z(1) = 0.0
809  do k=1,nz
810  dmass = gv%H_to_RZ * h(k)
811  dpres = us%L_to_Z**2 * gv%g_Earth * dmass ! Equivalent to GV%H_to_Pa * h(k) with rescaling
812  dt_to_dpe(k) = (dmass * (pres_z(k) + 0.5*dpres)) * dsv_dt(k)
813  ds_to_dpe(k) = (dmass * (pres_z(k) + 0.5*dpres)) * dsv_ds(k)
814  dt_to_dcolht(k) = dmass * dsv_dt(k)
815  ds_to_dcolht(k) = dmass * dsv_ds(k)
816 
817  pres_z(k+1) = pres_z(k) + dpres
818  enddo
819 
820  ! Determine the total thickness (h_sum) and the fractional distance from the bottom (hb_hs).
821  h_sum = h_neglect ; do k=1,nz ; h_sum = h_sum + h(k) ; enddo
822  i_hs = 0.0 ; if (h_sum > 0.0) i_hs = 1.0 / h_sum
823  h_bot = 0.0
824  hb_hs(nz+1) = 0.0
825  do k=nz,1,-1
826  h_bot = h_bot + h(k)
827  hb_hs(k) = h_bot * i_hs
828  enddo
829 
830  mld_output = h(1)*gv%H_to_Z
831 
832  !/The following lines are for the iteration over MLD
833  ! max_MLD will initialized as ocean bottom depth
834  max_mld = 0.0 ; do k=1,nz ; max_mld = max_mld + h(k)*gv%H_to_Z ; enddo
835  !min_MLD will initialize as 0.
836  min_mld = 0.0
837 
838  ! If no first guess is provided for MLD, try the middle of the water column
839  if (mld_guess <= min_mld) mld_guess = 0.5 * (min_mld + max_mld)
840 
841  ! Iterate to determine a converged EPBL depth.
842  obl_converged = .false.
843  do obl_it=1,cs%Max_MLD_Its
844 
845  if (.not. obl_converged) then
846  ! If not using MLD_Iteration flag loop to only execute once.
847  if (.not.cs%Use_MLD_iteration) obl_converged = .true.
848 
849  if (debug) then ; mech_tke_k(:) = 0.0 ; conv_perel_k(:) = 0.0 ; endif
850 
851  ! Reset ML_depth
852  mld_output = h(1)*gv%H_to_Z
853  sfc_connected = .true.
854 
855  !/ Here we get MStar, which is the ratio of convective TKE driven mixing to UStar**3
856  if (cs%Use_LT) then
857  call get_langmuir_number(la, g, gv, us, abs(mld_guess), u_star_mean, i, j, &
858  h=h, u_h=u, v_h=v, waves=waves)
859  call find_mstar(cs, us, b_flux, u_star, u_star_mean, mld_guess, absf, &
860  mstar_total, langmuir_number=la, convect_langmuir_number=lamod,&
861  mstar_lt=mstar_lt)
862  else
863  call find_mstar(cs, us, b_flux, u_star, u_star_mean, mld_guess, absf, mstar_total)
864  endif
865 
866  !/ Apply MStar to get mech_TKE
867  if ((cs%answers_2018) .and. (cs%mstar_scheme==use_fixed_mstar)) then
868  mech_tke = (dt*mstar_total*gv%Rho0) * u_star**3
869  else
870  mech_tke = mstar_total * (dt*gv%Rho0* u_star**3)
871  endif
872 
873  if (cs%TKE_diagnostics) then
874  ecd%dTKE_conv = 0.0 ; ecd%dTKE_mixing = 0.0
875  ecd%dTKE_MKE = 0.0 ; ecd%dTKE_mech_decay = 0.0 ; ecd%dTKE_conv_decay = 0.0
876 
877  ecd%dTKE_wind = mech_tke * i_dtdiag
878  if (tke_forcing(1) <= 0.0) then
879  ecd%dTKE_forcing = max(-mech_tke, tke_forcing(1)) * i_dtdiag
880  ! eCD%dTKE_unbalanced = min(0.0, TKE_forcing(1) + mech_TKE) * I_dtdiag
881  else
882  ecd%dTKE_forcing = cs%nstar*tke_forcing(1) * i_dtdiag
883  ! eCD%dTKE_unbalanced = 0.0
884  endif
885  endif
886 
887  if (tke_forcing(1) <= 0.0) then
888  mech_tke = mech_tke + tke_forcing(1)
889  if (mech_tke < 0.0) mech_tke = 0.0
890  conv_perel = 0.0
891  else
892  conv_perel = tke_forcing(1)
893  endif
894 
895 
896  ! Store in 1D arrays for output.
897  do k=1,nz+1 ; mixvel(k) = 0.0 ; mixlen(k) = 0.0 ; enddo
898 
899  ! Determine the mixing shape function MixLen_shape.
900  if ((.not.cs%Use_MLD_iteration) .or. &
901  (cs%transLay_scale >= 1.0) .or. (cs%transLay_scale < 0.0) ) then
902  do k=1,nz+1
903  mixlen_shape(k) = 1.0
904  enddo
905  elseif (mld_guess <= 0.0) then
906  if (cs%transLay_scale > 0.0) then ; do k=1,nz+1
907  mixlen_shape(k) = cs%transLay_scale
908  enddo ; else ; do k=1,nz+1
909  mixlen_shape(k) = 1.0
910  enddo ; endif
911  else
912  ! Reduce the mixing length based on MLD, with a quadratic
913  ! expression that follows KPP.
914  i_mld = 1.0 / mld_guess
915  h_rsum = 0.0
916  mixlen_shape(1) = 1.0
917  do k=2,nz+1
918  h_rsum = h_rsum + h(k-1)*gv%H_to_Z
919  if (cs%MixLenExponent==2.0) then
920  mixlen_shape(k) = cs%transLay_scale + (1.0 - cs%transLay_scale) * &
921  (max(0.0, (mld_guess - h_rsum)*i_mld) )**2 ! CS%MixLenExponent
922  else
923  mixlen_shape(k) = cs%transLay_scale + (1.0 - cs%transLay_scale) * &
924  (max(0.0, (mld_guess - h_rsum)*i_mld) )**cs%MixLenExponent
925  endif
926  enddo
927  endif
928 
929  kd(1) = 0.0 ; kddt_h(1) = 0.0
930  hp_a = h(1)
931  dt_to_dpe_a(1) = dt_to_dpe(1) ; dt_to_dcolht_a(1) = dt_to_dcolht(1)
932  ds_to_dpe_a(1) = ds_to_dpe(1) ; ds_to_dcolht_a(1) = ds_to_dcolht(1)
933 
934  htot = h(1) ; uhtot = u(1)*h(1) ; vhtot = v(1)*h(1)
935 
936  if (debug) then
937  mech_tke_k(1) = mech_tke ; conv_perel_k(1) = conv_perel
938  nstar_k(:) = 0.0 ; nstar_k(1) = cs%nstar ; num_itts(:) = -1
939  endif
940 
941  do k=2,nz
942  ! Apply dissipation to the TKE, here applied as an exponential decay
943  ! due to 3-d turbulent energy being lost to inefficient rotational modes.
944 
945  ! There should be several different "flavors" of TKE that decay at
946  ! different rates. The following form is often used for mechanical
947  ! stirring from the surface, perhaps due to breaking surface gravity
948  ! waves and wind-driven turbulence.
949  idecay_len_tke = (cs%TKE_decay * absf / u_star) * gv%H_to_Z
950  exp_kh = 1.0
951  if (idecay_len_tke > 0.0) exp_kh = exp(-h(k-1)*idecay_len_tke)
952  if (cs%TKE_diagnostics) &
953  ecd%dTKE_mech_decay = ecd%dTKE_mech_decay + (exp_kh-1.0) * mech_tke * i_dtdiag
954  mech_tke = mech_tke * exp_kh
955 
956  ! Accumulate any convectively released potential energy to contribute
957  ! to wstar and to drive penetrating convection.
958  if (tke_forcing(k) > 0.0) then
959  conv_perel = conv_perel + tke_forcing(k)
960  if (cs%TKE_diagnostics) &
961  ecd%dTKE_forcing = ecd%dTKE_forcing + cs%nstar*tke_forcing(k) * i_dtdiag
962  endif
963 
964  if (debug) then
965  mech_tke_k(k) = mech_tke ; conv_perel_k(k) = conv_perel
966  endif
967 
968  ! Determine the total energy
969  nstar_fc = cs%nstar
970  if (cs%nstar * conv_perel > 0.0) then
971  ! Here nstar is a function of the natural Rossby number 0.2/(1+0.2/Ro), based
972  ! on a curve fit from the data of Wang (GRL, 2003).
973  ! Note: Ro = 1.0 / sqrt(0.5 * dt * Rho0 * (absf*htot)**3 / conv_PErel)
974  nstar_fc = cs%nstar * conv_perel / (conv_perel + 0.2 * &
975  sqrt(0.5 * dt * gv%Rho0 * (absf*(htot*gv%H_to_Z))**3 * conv_perel))
976  endif
977 
978  if (debug) nstar_k(k) = nstar_fc
979 
980  tot_tke = mech_tke + nstar_fc * conv_perel
981 
982  ! For each interior interface, first discard the TKE to account for
983  ! mixing of shortwave radiation through the next denser cell.
984  if (tke_forcing(k) < 0.0) then
985  if (tke_forcing(k) + tot_tke < 0.0) then
986  ! The shortwave requirements deplete all the energy in this layer.
987  if (cs%TKE_diagnostics) then
988  ecd%dTKE_mixing = ecd%dTKE_mixing + tot_tke * i_dtdiag
989  ecd%dTKE_forcing = ecd%dTKE_forcing - tot_tke * i_dtdiag
990  ! eCD%dTKE_unbalanced = eCD%dTKE_unbalanced + (TKE_forcing(k) + tot_TKE) * I_dtdiag
991  ecd%dTKE_conv_decay = ecd%dTKE_conv_decay + (cs%nstar-nstar_fc) * conv_perel * i_dtdiag
992  endif
993  tot_tke = 0.0 ; mech_tke = 0.0 ; conv_perel = 0.0
994  else
995  ! Reduce the mechanical and convective TKE proportionately.
996  tke_reduc = (tot_tke + tke_forcing(k)) / tot_tke
997  if (cs%TKE_diagnostics) then
998  ecd%dTKE_mixing = ecd%dTKE_mixing - tke_forcing(k) * i_dtdiag
999  ecd%dTKE_forcing = ecd%dTKE_forcing + tke_forcing(k) * i_dtdiag
1000  ecd%dTKE_conv_decay = ecd%dTKE_conv_decay + &
1001  (1.0-tke_reduc)*(cs%nstar-nstar_fc) * conv_perel * i_dtdiag
1002  endif
1003  tot_tke = tke_reduc*tot_tke ! = tot_TKE + TKE_forcing(k)
1004  mech_tke = tke_reduc*mech_tke
1005  conv_perel = tke_reduc*conv_perel
1006  endif
1007  endif
1008 
1009  ! Precalculate some temporary expressions that are independent of Kddt_h(K).
1010  if (cs%orig_PE_calc) then
1011  if (k==2) then
1012  dte_t2 = 0.0 ; dse_t2 = 0.0
1013  else
1014  dte_t2 = kddt_h(k-1) * ((t0(k-2) - t0(k-1)) + dte(k-2))
1015  dse_t2 = kddt_h(k-1) * ((s0(k-2) - s0(k-1)) + dse(k-2))
1016  endif
1017  endif
1018  dt_h = (gv%Z_to_H**2*dt) / max(0.5*(h(k-1)+h(k)), 1e-15*h_sum)
1019 
1020  ! This tests whether the layers above and below this interface are in
1021  ! a convetively stable configuration, without considering any effects of
1022  ! mixing at higher interfaces. It is an approximation to the more
1023  ! complete test dPEc_dKd_Kd0 >= 0.0, that would include the effects of
1024  ! mixing across interface K-1. The dT_to_dColHt here are effectively
1025  ! mass-weigted estimates of dSV_dT.
1026  convectively_stable = ( 0.0 <= &
1027  ( (dt_to_dcolht(k) + dt_to_dcolht(k-1) ) * (t0(k-1)-t0(k)) + &
1028  (ds_to_dcolht(k) + ds_to_dcolht(k-1) ) * (s0(k-1)-s0(k)) ) )
1029 
1030  if ((mech_tke + conv_perel) <= 0.0 .and. convectively_stable) then
1031  ! Energy is already exhausted, so set Kd = 0 and cycle or exit?
1032  tot_tke = 0.0 ; mech_tke = 0.0 ; conv_perel = 0.0
1033  kd(k) = 0.0 ; kddt_h(k) = 0.0
1034  sfc_disconnect = .true.
1035  ! if (.not.debug) exit
1036 
1037  ! The estimated properties for layer k-1 can be calculated, using
1038  ! greatly simplified expressions when Kddt_h = 0. This enables the
1039  ! tridiagonal solver for the whole column to be completed for debugging
1040  ! purposes, and also allows for something akin to convective adjustment
1041  ! in unstable interior regions?
1042  b1 = 1.0 / hp_a
1043  c1(k) = 0.0
1044  if (cs%orig_PE_calc) then
1045  dte(k-1) = b1 * ( dte_t2 )
1046  dse(k-1) = b1 * ( dse_t2 )
1047  endif
1048 
1049  hp_a = h(k)
1050  dt_to_dpe_a(k) = dt_to_dpe(k)
1051  ds_to_dpe_a(k) = ds_to_dpe(k)
1052  dt_to_dcolht_a(k) = dt_to_dcolht(k)
1053  ds_to_dcolht_a(k) = ds_to_dcolht(k)
1054 
1055  else ! tot_TKE > 0.0 or this is a potentially convectively unstable profile.
1056  sfc_disconnect = .false.
1057 
1058  ! Precalculate some more temporary expressions that are independent of
1059  ! Kddt_h(K).
1060  if (cs%orig_PE_calc) then
1061  if (k==2) then
1062  dt_km1_t2 = (t0(k)-t0(k-1))
1063  ds_km1_t2 = (s0(k)-s0(k-1))
1064  else
1065  dt_km1_t2 = (t0(k)-t0(k-1)) - &
1066  (kddt_h(k-1) / hp_a) * ((t0(k-2) - t0(k-1)) + dte(k-2))
1067  ds_km1_t2 = (s0(k)-s0(k-1)) - &
1068  (kddt_h(k-1) / hp_a) * ((s0(k-2) - s0(k-1)) + dse(k-2))
1069  endif
1070  dte_term = dte_t2 + hp_a * (t0(k-1)-t0(k))
1071  dse_term = dse_t2 + hp_a * (s0(k-1)-s0(k))
1072  else
1073  if (k<=2) then
1074  th_a(k-1) = h(k-1) * t0(k-1) ; sh_a(k-1) = h(k-1) * s0(k-1)
1075  else
1076  th_a(k-1) = h(k-1) * t0(k-1) + kddt_h(k-1) * te(k-2)
1077  sh_a(k-1) = h(k-1) * s0(k-1) + kddt_h(k-1) * se(k-2)
1078  endif
1079  th_b(k) = h(k) * t0(k) ; sh_b(k) = h(k) * s0(k)
1080  endif
1081 
1082  ! Using Pr=1 and the diffusivity at the bottom interface (once it is
1083  ! known), determine how much resolved mean kinetic energy (MKE) will be
1084  ! extracted within a timestep and add a fraction CS%MKE_to_TKE_effic of
1085  ! this to the mTKE budget available for mixing in the next layer.
1086 
1087  if ((cs%MKE_to_TKE_effic > 0.0) .and. (htot*h(k) > 0.0)) then
1088  ! This is the energy that would be available from homogenizing the
1089  ! velocities between layer k and the layers above.
1090  dmke_max = (us%L_to_Z**2*gv%H_to_RZ * cs%MKE_to_TKE_effic) * 0.5 * &
1091  (h(k) / ((htot + h(k))*htot)) * &
1092  ((uhtot-u(k)*htot)**2 + (vhtot-v(k)*htot)**2)
1093  ! A fraction (1-exp(Kddt_h*MKE2_Hharm)) of this energy would be
1094  ! extracted by mixing with a finite viscosity.
1095  mke2_hharm = (htot + h(k) + 2.0*h_neglect) / &
1096  ((htot+h_neglect) * (h(k)+h_neglect))
1097  else
1098  dmke_max = 0.0
1099  mke2_hharm = 0.0
1100  endif
1101 
1102  ! At this point, Kddt_h(K) will be unknown because its value may depend
1103  ! on how much energy is available. mech_TKE might be negative due to
1104  ! contributions from TKE_forced.
1105  h_tt = htot + h_tt_min
1106  tke_here = mech_tke + cs%wstar_ustar_coef*conv_perel
1107  if (tke_here > 0.0) then
1108  if (cs%wT_scheme==wt_from_croot_tke) then
1109  vstar = cs%vstar_scale_fac * vstar_unit_scale * (i_dtrho*tke_here)**c1_3
1110  elseif (cs%wT_scheme==wt_from_rh18) then
1111  surface_scale = max(0.05, 1.0 - htot/mld_guess)
1112  vstar = cs%vstar_scale_fac * surface_scale * (cs%vstar_surf_fac*u_star + &
1113  vstar_unit_scale * (cs%wstar_ustar_coef*conv_perel*i_dtrho)**c1_3)
1114  endif
1115  hbs_here = gv%H_to_Z * min(hb_hs(k), mixlen_shape(k))
1116  mixlen(k) = max(cs%min_mix_len, ((h_tt*hbs_here)*vstar) / &
1117  ((cs%Ekman_scale_coef * absf) * (h_tt*hbs_here) + vstar))
1118  !Note setting Kd_guess0 to vstar * CS%vonKar * mixlen(K) here will
1119  ! change the answers. Therefore, skipping that.
1120  if (.not.cs%Use_MLD_iteration) then
1121  kd_guess0 = vstar * cs%vonKar * ((h_tt*hbs_here)*vstar) / &
1122  ((cs%Ekman_scale_coef * absf) * (h_tt*hbs_here) + vstar)
1123  else
1124  kd_guess0 = vstar * cs%vonKar * mixlen(k)
1125  endif
1126  else
1127  vstar = 0.0 ; kd_guess0 = 0.0
1128  endif
1129  mixvel(k) = vstar ! Track vstar
1130  kddt_h_g0 = kd_guess0 * dt_h
1131 
1132  if (cs%orig_PE_calc) then
1133  call find_pe_chg_orig(kddt_h_g0, h(k), hp_a, dte_term, dse_term, &
1134  dt_km1_t2, ds_km1_t2, dt_to_dpe(k), ds_to_dpe(k), &
1135  dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), &
1136  pres_z(k), dt_to_dcolht(k), ds_to_dcolht(k), &
1137  dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1138  pe_chg=pe_chg_g0, dpec_dkd=dpea_dkd_g0, dpe_max=pe_chg_max, &
1139  dpec_dkd_0=dpec_dkd_kd0 )
1140  else
1141  call find_pe_chg(0.0, kddt_h_g0, hp_a, h(k), &
1142  th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
1143  dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe(k), ds_to_dpe(k), &
1144  pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1145  dt_to_dcolht(k), ds_to_dcolht(k), &
1146  pe_chg=pe_chg_g0, dpec_dkd=dpea_dkd_g0, dpe_max=pe_chg_max, &
1147  dpec_dkd_0=dpec_dkd_kd0 )
1148  endif
1149 
1150  mke_src = dmke_max*(1.0 - exp(-kddt_h_g0 * mke2_hharm))
1151 
1152  ! This block checks out different cases to determine Kd at the present interface.
1153  if ((pe_chg_g0 < 0.0) .or. ((vstar == 0.0) .and. (dpec_dkd_kd0 < 0.0))) then
1154  ! This column is convectively unstable.
1155  if (pe_chg_max <= 0.0) then
1156  ! Does MKE_src need to be included in the calculation of vstar here?
1157  tke_here = mech_tke + cs%wstar_ustar_coef*(conv_perel-pe_chg_max)
1158  if (tke_here > 0.0) then
1159  if (cs%wT_scheme==wt_from_croot_tke) then
1160  vstar = cs%vstar_scale_fac * vstar_unit_scale * (i_dtrho*tke_here)**c1_3
1161  elseif (cs%wT_scheme==wt_from_rh18) then
1162  surface_scale = max(0.05, 1. - htot/mld_guess)
1163  vstar = cs%vstar_scale_fac * surface_scale * (cs%vstar_surf_fac*u_star + &
1164  vstar_unit_scale * (cs%wstar_ustar_coef*conv_perel*i_dtrho)**c1_3)
1165  endif
1166  hbs_here = gv%H_to_Z * min(hb_hs(k), mixlen_shape(k))
1167  mixlen(k) = max(cs%min_mix_len, ((h_tt*hbs_here)*vstar) / &
1168  ((cs%Ekman_scale_coef * absf) * (h_tt*hbs_here) + vstar))
1169  if (.not.cs%Use_MLD_iteration) then
1170  ! Note again (as prev) that using mixlen here
1171  ! instead of redoing the computation will change answers...
1172  kd(k) = vstar * cs%vonKar * ((h_tt*hbs_here)*vstar) / &
1173  ((cs%Ekman_scale_coef * absf) * (h_tt*hbs_here) + vstar)
1174  else
1175  kd(k) = vstar * cs%vonKar * mixlen(k)
1176  endif
1177  else
1178  vstar = 0.0 ; kd(k) = 0.0
1179  endif
1180  mixvel(k) = vstar
1181 
1182  if (cs%orig_PE_calc) then
1183  call find_pe_chg_orig(kd(k)*dt_h, h(k), hp_a, dte_term, dse_term, &
1184  dt_km1_t2, ds_km1_t2, dt_to_dpe(k), ds_to_dpe(k), &
1185  dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), &
1186  pres_z(k), dt_to_dcolht(k), ds_to_dcolht(k), &
1187  dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1188  pe_chg=dpe_conv)
1189  else
1190  call find_pe_chg(0.0, kd(k)*dt_h, hp_a, h(k), &
1191  th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
1192  dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe(k), ds_to_dpe(k), &
1193  pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1194  dt_to_dcolht(k), ds_to_dcolht(k), &
1195  pe_chg=dpe_conv)
1196  endif
1197  ! Should this be iterated to convergence for Kd?
1198  if (dpe_conv > 0.0) then
1199  kd(k) = kd_guess0 ; dpe_conv = pe_chg_g0
1200  else
1201  mke_src = dmke_max*(1.0 - exp(-(kd(k)*dt_h) * mke2_hharm))
1202  endif
1203  else
1204  ! The energy change does not vary monotonically with Kddt_h. Find the maximum?
1205  kd(k) = kd_guess0 ; dpe_conv = pe_chg_g0
1206  endif
1207 
1208  conv_perel = conv_perel - dpe_conv
1209  mech_tke = mech_tke + mke_src
1210  if (cs%TKE_diagnostics) then
1211  ecd%dTKE_conv = ecd%dTKE_conv - cs%nstar*dpe_conv * i_dtdiag
1212  ecd%dTKE_MKE = ecd%dTKE_MKE + mke_src * i_dtdiag
1213  endif
1214  if (sfc_connected) then
1215  mld_output = mld_output + gv%H_to_Z * h(k)
1216  endif
1217 
1218  kddt_h(k) = kd(k) * dt_h
1219  elseif (tot_tke + (mke_src - pe_chg_g0) >= 0.0) then
1220  ! This column is convctively stable and there is energy to support the suggested
1221  ! mixing. Keep that estimate.
1222  kd(k) = kd_guess0
1223  kddt_h(k) = kddt_h_g0
1224 
1225  ! Reduce the mechanical and convective TKE proportionately.
1226  tot_tke = tot_tke + mke_src
1227  tke_reduc = 0.0 ! tot_TKE could be 0 if Convectively_stable is false.
1228  if (tot_tke > 0.0) tke_reduc = (tot_tke - pe_chg_g0) / tot_tke
1229  if (cs%TKE_diagnostics) then
1230  ecd%dTKE_mixing = ecd%dTKE_mixing - pe_chg_g0 * i_dtdiag
1231  ecd%dTKE_MKE = ecd%dTKE_MKE + mke_src * i_dtdiag
1232  ecd%dTKE_conv_decay = ecd%dTKE_conv_decay + &
1233  (1.0-tke_reduc)*(cs%nstar-nstar_fc) * conv_perel * i_dtdiag
1234  endif
1235  tot_tke = tke_reduc*tot_tke
1236  mech_tke = tke_reduc*(mech_tke + mke_src)
1237  conv_perel = tke_reduc*conv_perel
1238  if (sfc_connected) then
1239  mld_output = mld_output + gv%H_to_Z * h(k)
1240  endif
1241 
1242  elseif (tot_tke == 0.0) then
1243  ! This can arise if nstar_FC = 0, but it is not common.
1244  kd(k) = 0.0 ; kddt_h(k) = 0.0
1245  tot_tke = 0.0 ; conv_perel = 0.0 ; mech_tke = 0.0
1246  sfc_disconnect = .true.
1247  else
1248  ! There is not enough energy to support the mixing, so reduce the
1249  ! diffusivity to what can be supported.
1250  kddt_h_max = kddt_h_g0 ; kddt_h_min = 0.0
1251  tke_left_max = tot_tke + (mke_src - pe_chg_g0)
1252  tke_left_min = tot_tke
1253 
1254  ! As a starting guess, take the minimum of a false position estimate
1255  ! and a Newton's method estimate starting from Kddt_h = 0.0.
1256  kddt_h_guess = tot_tke * kddt_h_max / max( pe_chg_g0 - mke_src, &
1257  kddt_h_max * (dpec_dkd_kd0 - dmke_max * mke2_hharm) )
1258  ! The above expression is mathematically the same as the following
1259  ! except it is not susceptible to division by zero when
1260  ! dPEc_dKd_Kd0 = dMKE_max = 0 .
1261  ! Kddt_h_guess = tot_TKE * min( Kddt_h_max / (PE_chg_g0 - MKE_src), &
1262  ! 1.0 / (dPEc_dKd_Kd0 - dMKE_max * MKE2_Hharm) )
1263  if (debug) then
1264  tke_left_itt(:) = 0.0 ; dpea_dkd_itt(:) = 0.0 ; pe_chg_itt(:) = 0.0
1265  mke_src_itt(:) = 0.0 ; kddt_h_itt(:) = 0.0
1266  endif
1267  do itt=1,max_itt
1268  if (cs%orig_PE_calc) then
1269  call find_pe_chg_orig(kddt_h_guess, h(k), hp_a, dte_term, dse_term, &
1270  dt_km1_t2, ds_km1_t2, dt_to_dpe(k), ds_to_dpe(k), &
1271  dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), &
1272  pres_z(k), dt_to_dcolht(k), ds_to_dcolht(k), &
1273  dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1274  pe_chg=pe_chg, dpec_dkd=dpec_dkd )
1275  else
1276  call find_pe_chg(0.0, kddt_h_guess, hp_a, h(k), &
1277  th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
1278  dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe(k), ds_to_dpe(k), &
1279  pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1280  dt_to_dcolht(k), ds_to_dcolht(k), &
1281  pe_chg=dpe_conv)
1282  endif
1283  mke_src = dmke_max * (1.0 - exp(-mke2_hharm * kddt_h_guess))
1284  dmke_src_dk = dmke_max * mke2_hharm * exp(-mke2_hharm * kddt_h_guess)
1285 
1286  tke_left = tot_tke + (mke_src - pe_chg)
1287  if (debug .and. itt<=20) then
1288  kddt_h_itt(itt) = kddt_h_guess ; mke_src_itt(itt) = mke_src
1289  pe_chg_itt(itt) = pe_chg ; dpea_dkd_itt(itt) = dpec_dkd
1290  tke_left_itt(itt) = tke_left
1291  endif
1292  ! Store the new bounding values, bearing in mind that min and max
1293  ! here refer to Kddt_h and dTKE_left/dKddt_h < 0:
1294  if (tke_left >= 0.0) then
1295  kddt_h_min = kddt_h_guess ; tke_left_min = tke_left
1296  else
1297  kddt_h_max = kddt_h_guess ; tke_left_max = tke_left
1298  endif
1299 
1300  ! Try to use Newton's method, but if it would go outside the bracketed
1301  ! values use the false-position method instead.
1302  use_newt = .true.
1303  if (dpec_dkd - dmke_src_dk <= 0.0) then
1304  use_newt = .false.
1305  else
1306  dkddt_h_newt = tke_left / (dpec_dkd - dmke_src_dk)
1307  kddt_h_newt = kddt_h_guess + dkddt_h_newt
1308  if ((kddt_h_newt > kddt_h_max) .or. (kddt_h_newt < kddt_h_min)) &
1309  use_newt = .false.
1310  endif
1311 
1312  if (use_newt) then
1313  kddt_h_next = kddt_h_guess + dkddt_h_newt
1314  dkddt_h = dkddt_h_newt
1315  else
1316  kddt_h_next = (tke_left_max * kddt_h_min - kddt_h_max * tke_left_min) / &
1317  (tke_left_max - tke_left_min)
1318  dkddt_h = kddt_h_next - kddt_h_guess
1319  endif
1320 
1321  if ((abs(dkddt_h) < 1e-9*kddt_h_guess) .or. (itt==max_itt)) then
1322  ! Use the old value so that the energy calculation does not need to be repeated.
1323  if (debug) num_itts(k) = itt
1324  exit
1325  else
1326  kddt_h_guess = kddt_h_next
1327  endif
1328  enddo ! Inner iteration loop on itt.
1329  kd(k) = kddt_h_guess / dt_h ; kddt_h(k) = kd(k) * dt_h
1330 
1331  ! All TKE should have been consumed.
1332  if (cs%TKE_diagnostics) then
1333  ecd%dTKE_mixing = ecd%dTKE_mixing - (tot_tke + mke_src) * i_dtdiag
1334  ecd%dTKE_MKE = ecd%dTKE_MKE + mke_src * i_dtdiag
1335  ecd%dTKE_conv_decay = ecd%dTKE_conv_decay + &
1336  (cs%nstar-nstar_fc) * conv_perel * i_dtdiag
1337  endif
1338 
1339  if (sfc_connected) mld_output = mld_output + &
1340  (pe_chg / (pe_chg_g0)) * gv%H_to_Z * h(k)
1341 
1342  tot_tke = 0.0 ; mech_tke = 0.0 ; conv_perel = 0.0
1343  sfc_disconnect = .true.
1344  endif ! End of convective or forced mixing cases to determine Kd.
1345 
1346  kddt_h(k) = kd(k) * dt_h
1347  ! At this point, the final value of Kddt_h(K) is known, so the
1348  ! estimated properties for layer k-1 can be calculated.
1349  b1 = 1.0 / (hp_a + kddt_h(k))
1350  c1(k) = kddt_h(k) * b1
1351  if (cs%orig_PE_calc) then
1352  dte(k-1) = b1 * ( kddt_h(k)*(t0(k)-t0(k-1)) + dte_t2 )
1353  dse(k-1) = b1 * ( kddt_h(k)*(s0(k)-s0(k-1)) + dse_t2 )
1354  endif
1355 
1356  hp_a = h(k) + (hp_a * b1) * kddt_h(k)
1357  dt_to_dpe_a(k) = dt_to_dpe(k) + c1(k)*dt_to_dpe_a(k-1)
1358  ds_to_dpe_a(k) = ds_to_dpe(k) + c1(k)*ds_to_dpe_a(k-1)
1359  dt_to_dcolht_a(k) = dt_to_dcolht(k) + c1(k)*dt_to_dcolht_a(k-1)
1360  ds_to_dcolht_a(k) = ds_to_dcolht(k) + c1(k)*ds_to_dcolht_a(k-1)
1361 
1362  endif ! tot_TKT > 0.0 branch. Kddt_h(K) has been set.
1363 
1364  ! Store integrated velocities and thicknesses for MKE conversion calculations.
1365  if (sfc_disconnect) then
1366  ! There is no turbulence at this interface, so zero out the running sums.
1367  uhtot = u(k)*h(k)
1368  vhtot = v(k)*h(k)
1369  htot = h(k)
1370  sfc_connected = .false.
1371  else
1372  uhtot = uhtot + u(k)*h(k)
1373  vhtot = vhtot + v(k)*h(k)
1374  htot = htot + h(k)
1375  endif
1376 
1377  if (debug) then
1378  if (k==2) then
1379  te(1) = b1*(h(1)*t0(1))
1380  se(1) = b1*(h(1)*s0(1))
1381  else
1382  te(k-1) = b1 * (h(k-1) * t0(k-1) + kddt_h(k-1) * te(k-2))
1383  se(k-1) = b1 * (h(k-1) * s0(k-1) + kddt_h(k-1) * se(k-2))
1384  endif
1385  endif
1386  enddo
1387  kd(nz+1) = 0.0
1388 
1389  if (debug) then
1390  ! Complete the tridiagonal solve for Te.
1391  b1 = 1.0 / hp_a
1392  te(nz) = b1 * (h(nz) * t0(nz) + kddt_h(nz) * te(nz-1))
1393  se(nz) = b1 * (h(nz) * s0(nz) + kddt_h(nz) * se(nz-1))
1394  ecd%dT_expect(nz) = te(nz) - t0(nz) ; ecd%dS_expect(nz) = se(nz) - s0(nz)
1395  do k=nz-1,1,-1
1396  te(k) = te(k) + c1(k+1)*te(k+1)
1397  se(k) = se(k) + c1(k+1)*se(k+1)
1398  ecd%dT_expect(k) = te(k) - t0(k) ; ecd%dS_expect(k) = se(k) - s0(k)
1399  enddo
1400 
1401  dpe_debug = 0.0
1402  do k=1,nz
1403  dpe_debug = dpe_debug + (dt_to_dpe(k) * (te(k) - t0(k)) + &
1404  ds_to_dpe(k) * (se(k) - s0(k)))
1405  enddo
1406  mixing_debug = dpe_debug * i_dtdiag
1407  endif
1408  k = nz ! This is here to allow a breakpoint to be set.
1409  !/BGR
1410  ! The following lines are used for the iteration
1411  ! note the iteration has been altered to use the value predicted by
1412  ! the TKE threshold (ML_DEPTH). This is because the MSTAR
1413  ! is now dependent on the ML, and therefore the ML needs to be estimated
1414  ! more precisely than the grid spacing.
1415 
1416  !New method uses ML_DEPTH as computed in ePBL routine
1417  mld_found = mld_output
1418  if (mld_found - cs%MLD_tol > mld_guess) then
1419  min_mld = mld_guess
1420  elseif (abs(mld_guess - mld_found) < cs%MLD_tol) then
1421  obl_converged = .true. ! Break convergence loop
1422  else
1423  max_mld = mld_guess ! We know this guess was too deep
1424  endif
1425 
1426  ! For next pass, guess average of minimum and maximum values.
1427  !### We should try using the false position method instead of simple bisection.
1428  mld_guess = 0.5*(min_mld + max_mld)
1429  endif
1430  enddo ! Iteration loop for converged boundary layer thickness.
1431  if (cs%Use_LT) then
1432  ecd%LA = la ; ecd%LAmod = lamod ; ecd%mstar = mstar_total ; ecd%mstar_LT = mstar_lt
1433  else
1434  ecd%LA = 0.0 ; ecd%LAmod = 0.0 ; ecd%mstar = mstar_total ; ecd%mstar_LT = 0.0
1435  endif
1436 
1437  mld_io = mld_output
1438 

References find_mstar(), find_pe_chg(), find_pe_chg_orig(), mom_wave_interface::get_langmuir_number(), mom_error_handler::mom_error(), use_fixed_mstar, wt_from_croot_tke, and wt_from_rh18.

Referenced by energetic_pbl().

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

◆ find_mstar()

subroutine mom_energetic_pbl::find_mstar ( type(energetic_pbl_cs), pointer  CS,
type(unit_scale_type), intent(in)  US,
real, intent(in)  Buoyancy_Flux,
real, intent(in)  UStar,
real, intent(in)  UStar_Mean,
real, intent(in)  BLD,
real, intent(in)  Abs_Coriolis,
real, intent(out)  MStar,
real, intent(in), optional  Langmuir_Number,
real, intent(out), optional  MStar_LT,
real, intent(out), optional  Convect_Langmuir_Number 
)
private

This subroutine finds the Mstar value for ePBL.

Parameters
csEnergetic_PBL control structure.
[in]usA dimensional unit scaling type
[in]ustarustar w/ gustiness [Z T-1 ~> m s-1]
[in]ustar_meanustar w/o gustiness [Z T-1 ~> m s-1]
[in]abs_coriolisabolute value of the Coriolis parameter [T-1 ~> s-1]
[in]buoyancy_fluxBuoyancy flux [Z2 T-3 ~> m2 s-3]
[in]bldboundary layer depth [Z ~> m]
[out]mstarOuput mstar (Mixing/ustar**3) [nondim]
[in]langmuir_numberLangmuir number [nondim]
[out]mstar_ltMstar increase due to Langmuir turbulence [nondim]
[out]convect_langmuir_numberLangmuir number including buoyancy flux [nondim]

Definition at line 1748 of file MOM_energetic_PBL.F90.

1748  type(energetic_PBL_CS), pointer :: CS !< Energetic_PBL control structure.
1749  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1750  real, intent(in) :: UStar !< ustar w/ gustiness [Z T-1 ~> m s-1]
1751  real, intent(in) :: UStar_Mean !< ustar w/o gustiness [Z T-1 ~> m s-1]
1752  real, intent(in) :: Abs_Coriolis !< abolute value of the Coriolis parameter [T-1 ~> s-1]
1753  real, intent(in) :: Buoyancy_Flux !< Buoyancy flux [Z2 T-3 ~> m2 s-3]
1754  real, intent(in) :: BLD !< boundary layer depth [Z ~> m]
1755  real, intent(out) :: Mstar !< Ouput mstar (Mixing/ustar**3) [nondim]
1756  real, optional, intent(in) :: Langmuir_Number !< Langmuir number [nondim]
1757  real, optional, intent(out) :: MStar_LT !< Mstar increase due to Langmuir turbulence [nondim]
1758  real, optional, intent(out) :: Convect_Langmuir_number !< Langmuir number including buoyancy flux [nondim]
1759 
1760  !/ Variables used in computing mstar
1761  real :: MSN_term ! Temporary terms [nondim]
1762  real :: MSCR_term1, MSCR_term2 ! Temporary terms [Z3 T-3 ~> m3 s-3]
1763  real :: MStar_Conv_Red ! Adjustment made to mstar due to convection reducing mechanical mixing [nondim]
1764  real :: MStar_S, MStar_N ! Mstar in (S)tabilizing/(N)ot-stabilizing buoyancy flux [nondim]
1765 
1766  !/ Integer options for how to find mstar
1767 
1768  !/
1769 
1770  if (cs%mstar_scheme == use_fixed_mstar) then
1771  mstar = cs%Fixed_MStar
1772  !/ 1. Get mstar
1773  elseif (cs%mstar_scheme == mstar_from_ekman) then
1774 
1775  if (cs%answers_2018) then
1776  ! The limit for the balance of rotation and stabilizing is f(L_Ekman,L_Obukhov)
1777  mstar_s = cs%MStar_coef*sqrt(max(0.0,buoyancy_flux) / ustar**2 / &
1778  (abs_coriolis + 1.e-10*us%T_to_s) )
1779  ! The limit for rotation (Ekman length) limited mixing
1780  mstar_n = cs%C_Ek * log( max( 1., ustar / (abs_coriolis + 1.e-10*us%T_to_s) / bld ) )
1781  else
1782  ! The limit for the balance of rotation and stabilizing is f(L_Ekman,L_Obukhov)
1783  mstar_s = cs%MSTAR_COEF*sqrt(max(0.0, buoyancy_flux) / (ustar**2 * max(abs_coriolis, 1.e-20*us%T_to_s)))
1784  ! The limit for rotation (Ekman length) limited mixing
1785  mstar_n = 0.0
1786  if (ustar > abs_coriolis * bld) mstar_n = cs%C_EK * log(ustar / (abs_coriolis * bld))
1787  endif
1788 
1789  ! Here 1.25 is about .5/von Karman, which gives the Obukhov limit.
1790  mstar = max(mstar_s, min(1.25, mstar_n))
1791  if (cs%MStar_Cap > 0.0) mstar = min( cs%MStar_Cap,mstar )
1792  elseif ( cs%mstar_scheme == mstar_from_rh18 ) then
1793  if (cs%answers_2018) then
1794  mstar_n = cs%RH18_MStar_cn1 * ( 1.0 - 1.0 / ( 1. + cs%RH18_MStar_cn2 * &
1795  exp( cs%RH18_mstar_CN3 * bld * abs_coriolis / ustar) ) )
1796  else
1797  msn_term = cs%RH18_MStar_cn2 * exp( cs%RH18_mstar_CN3 * bld * abs_coriolis / ustar)
1798  mstar_n = (cs%RH18_MStar_cn1 * msn_term) / ( 1. + msn_term)
1799  endif
1800  mstar_s = cs%RH18_MStar_CS1 * ( max(0.0, buoyancy_flux)**2 * bld / &
1801  ( ustar**5 * max(abs_coriolis,1.e-20*us%T_to_s) ) )**cs%RH18_mstar_cs2
1802  mstar = mstar_n + mstar_s
1803  endif
1804 
1805  !/ 2. Adjust mstar to account for convective turbulence
1806  if (cs%answers_2018) then
1807  mstar_conv_red = 1. - cs%MStar_Convect_coef * (-min(0.0,buoyancy_flux) + 1.e-10*us%T_to_s**3*us%m_to_Z**2) / &
1808  ( (-min(0.0,buoyancy_flux) + 1.e-10*us%T_to_s**3*us%m_to_Z**2) + &
1809  2.0 *mstar * ustar**3 / bld )
1810  else
1811  mscr_term1 = -bld * min(0.0, buoyancy_flux)
1812  mscr_term2 = 2.0*mstar * ustar**3
1813  if ( abs(mscr_term2) > 0.0) then
1814  mstar_conv_red = ((1.-cs%mstar_convect_coef) * mscr_term1 + mscr_term2) / (mscr_term1 + mscr_term2)
1815  else
1816  mstar_conv_red = 1.-cs%mstar_convect_coef
1817  endif
1818  endif
1819 
1820  !/3. Combine various mstar terms to get final value
1821  mstar = mstar * mstar_conv_red
1822 
1823  if (present(langmuir_number)) then
1824  !### In this call, ustar was previously ustar_mean. Is this change deliberate?
1825  call mstar_langmuir(cs, us, abs_coriolis, buoyancy_flux, ustar, bld, langmuir_number, mstar, &
1826  mstar_lt, convect_langmuir_number)
1827  endif
1828 

References mstar_from_ekman, mstar_from_rh18, mstar_langmuir(), and use_fixed_mstar.

Referenced by epbl_column().

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

◆ find_pe_chg()

subroutine mom_energetic_pbl::find_pe_chg ( real, intent(in)  Kddt_h0,
real, intent(in)  dKddt_h,
real, intent(in)  hp_a,
real, intent(in)  hp_b,
real, intent(in)  Th_a,
real, intent(in)  Sh_a,
real, intent(in)  Th_b,
real, intent(in)  Sh_b,
real, intent(in)  dT_to_dPE_a,
real, intent(in)  dS_to_dPE_a,
real, intent(in)  dT_to_dPE_b,
real, intent(in)  dS_to_dPE_b,
real, intent(in)  pres_Z,
real, intent(in)  dT_to_dColHt_a,
real, intent(in)  dS_to_dColHt_a,
real, intent(in)  dT_to_dColHt_b,
real, intent(in)  dS_to_dColHt_b,
real, intent(out), optional  PE_chg,
real, intent(out), optional  dPEc_dKd,
real, intent(out), optional  dPE_max,
real, intent(out), optional  dPEc_dKd_0,
real, intent(out), optional  PE_ColHt_cor 
)
private

This subroutine calculates the change in potential energy and or derivatives for several changes in an interfaces's diapycnal diffusivity times a timestep.

Parameters
[in]kddt_h0The previously used diffusivity at an interface times the time step and divided by the average of the thicknesses around the interface [H ~> m or kg m-2].
[in]dkddt_hThe trial change in the diffusivity at an interface times the time step and divided by the average of the thicknesses around the interface [H ~> m or kg m-2].
[in]hp_aThe effective pivot thickness of the layer above the interface, given by h_k plus a term that is a fraction (determined from the tridiagonal solver) of Kddt_h for the interface above [H ~> m or kg m-2].
[in]hp_bThe effective pivot thickness of the layer below the interface, given by h_k plus a term that is a fraction (determined from the tridiagonal solver) of Kddt_h for the interface above [H ~> m or kg m-2].
[in]th_aAn effective temperature times a thickness in the layer above, including implicit mixing effects with other yet higher layers [degC H ~> degC m or degC kg m-2].
[in]sh_aAn effective salinity times a thickness in the layer above, including implicit mixing effects with other yet higher layers [degC H ~> degC m or degC kg m-2].
[in]th_bAn effective temperature times a thickness in the layer below, including implicit mixfing effects with other yet lower layers [degC H ~> degC m or degC kg m-2].
[in]sh_bAn effective salinity times a thickness in the layer below, including implicit mixing effects with other yet lower layers [degC H ~> degC m or degC kg m-2].
[in]dt_to_dpe_aA factor (pres_lay*mass_lay*dSpec_vol/dT) relating a layer's temperature change to the change in column potential energy, including all implicit diffusive changes in the temperatures of all the layers above [R Z3 T-2 degC-1 ~> J m-2 degC-1].
[in]ds_to_dpe_aA factor (pres_lay*mass_lay*dSpec_vol/dS) relating a layer's salinity change to the change in column potential energy, including all implicit diffusive changes in the salinities of all the layers above [R Z3 T-2 ppt-1 ~> J m-2 ppt-1].
[in]dt_to_dpe_bA factor (pres_lay*mass_lay*dSpec_vol/dT) relating a layer's temperature change to the change in column potential energy, including all implicit diffusive changes in the temperatures of all the layers below [R Z3 T-2 degC-1 ~> J m-2 degC-1].
[in]ds_to_dpe_bA factor (pres_lay*mass_lay*dSpec_vol/dS) relating a layer's salinity change to the change in column potential energy, including all implicit diffusive changes in the salinities of all the layers below [R Z3 T-2 ppt-1 ~> J m-2 ppt-1].
[in]pres_zThe rescaled hydrostatic interface pressure, which relates the changes in column thickness to the energy that is radiated as gravity waves and unavailable to drive mixing [R Z2 T-2 ~> J m-3].
[in]dt_to_dcolht_aA factor (mass_lay*dSColHtc_vol/dT) relating a layer's temperature change to the change in column height, including all implicit diffusive changes in the temperatures of all the layers above [Z degC-1 ~> m degC-1].
[in]ds_to_dcolht_aA factor (mass_lay*dSColHtc_vol/dS) relating a layer's salinity change to the change in column height, including all implicit diffusive changes in the salinities of all the layers above [Z ppt-1 ~> m ppt-1].
[in]dt_to_dcolht_bA factor (mass_lay*dSColHtc_vol/dT) relating a layer's temperature change to the change in column height, including all implicit diffusive changes in the temperatures of all the layers below [Z degC-1 ~> m degC-1].
[in]ds_to_dcolht_bA factor (mass_lay*dSColHtc_vol/dS) relating a layer's salinity change to the change in column height, including all implicit diffusive changes in the salinities of all the layers below [Z ppt-1 ~> m ppt-1].
[out]pe_chgThe change in column potential energy from applying Kddt_h at the present interface [R Z3 T-2 ~> J m-2].
[out]dpec_dkdThe partial derivative of PE_chg with Kddt_h [R Z3 T-2 H-1 ~> J m-3 or J kg-1].
[out]dpe_maxThe maximum change in column potential energy that could be realizedd by applying a huge value of Kddt_h at the present interface [R Z3 T-2 ~> J m-2].
[out]dpec_dkd_0The partial derivative of PE_chg with Kddt_h in the limit where Kddt_h = 0 [R Z3 T-2 H-1 ~> J m-3 or J kg-1].
[out]pe_colht_corThe correction to PE_chg that is made due to a net change in the column height [R Z3 T-2 ~> J m-2].

Definition at line 1447 of file MOM_energetic_PBL.F90.

1447  real, intent(in) :: Kddt_h0 !< The previously used diffusivity at an interface times
1448  !! the time step and divided by the average of the
1449  !! thicknesses around the interface [H ~> m or kg m-2].
1450  real, intent(in) :: dKddt_h !< The trial change in the diffusivity at an interface times
1451  !! the time step and divided by the average of the
1452  !! thicknesses around the interface [H ~> m or kg m-2].
1453  real, intent(in) :: hp_a !< The effective pivot thickness of the layer above the
1454  !! interface, given by h_k plus a term that
1455  !! is a fraction (determined from the tridiagonal solver) of
1456  !! Kddt_h for the interface above [H ~> m or kg m-2].
1457  real, intent(in) :: hp_b !< The effective pivot thickness of the layer below the
1458  !! interface, given by h_k plus a term that
1459  !! is a fraction (determined from the tridiagonal solver) of
1460  !! Kddt_h for the interface above [H ~> m or kg m-2].
1461  real, intent(in) :: Th_a !< An effective temperature times a thickness in the layer
1462  !! above, including implicit mixing effects with other
1463  !! yet higher layers [degC H ~> degC m or degC kg m-2].
1464  real, intent(in) :: Sh_a !< An effective salinity times a thickness in the layer
1465  !! above, including implicit mixing effects with other
1466  !! yet higher layers [degC H ~> degC m or degC kg m-2].
1467  real, intent(in) :: Th_b !< An effective temperature times a thickness in the layer
1468  !! below, including implicit mixfing effects with other
1469  !! yet lower layers [degC H ~> degC m or degC kg m-2].
1470  real, intent(in) :: Sh_b !< An effective salinity times a thickness in the layer
1471  !! below, including implicit mixing effects with other
1472  !! yet lower layers [degC H ~> degC m or degC kg m-2].
1473  real, intent(in) :: dT_to_dPE_a !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating
1474  !! a layer's temperature change to the change in column potential
1475  !! energy, including all implicit diffusive changes in the
1476  !! temperatures of all the layers above [R Z3 T-2 degC-1 ~> J m-2 degC-1].
1477  real, intent(in) :: dS_to_dPE_a !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating
1478  !! a layer's salinity change to the change in column potential
1479  !! energy, including all implicit diffusive changes in the
1480  !! salinities of all the layers above [R Z3 T-2 ppt-1 ~> J m-2 ppt-1].
1481  real, intent(in) :: dT_to_dPE_b !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating
1482  !! a layer's temperature change to the change in column potential
1483  !! energy, including all implicit diffusive changes in the
1484  !! temperatures of all the layers below [R Z3 T-2 degC-1 ~> J m-2 degC-1].
1485  real, intent(in) :: dS_to_dPE_b !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating
1486  !! a layer's salinity change to the change in column potential
1487  !! energy, including all implicit diffusive changes in the
1488  !! salinities of all the layers below [R Z3 T-2 ppt-1 ~> J m-2 ppt-1].
1489  real, intent(in) :: pres_Z !< The rescaled hydrostatic interface pressure, which relates
1490  !! the changes in column thickness to the energy that is radiated
1491  !! as gravity waves and unavailable to drive mixing [R Z2 T-2 ~> J m-3].
1492  real, intent(in) :: dT_to_dColHt_a !< A factor (mass_lay*dSColHtc_vol/dT) relating
1493  !! a layer's temperature change to the change in column
1494  !! height, including all implicit diffusive changes
1495  !! in the temperatures of all the layers above [Z degC-1 ~> m degC-1].
1496  real, intent(in) :: dS_to_dColHt_a !< A factor (mass_lay*dSColHtc_vol/dS) relating
1497  !! a layer's salinity change to the change in column
1498  !! height, including all implicit diffusive changes
1499  !! in the salinities of all the layers above [Z ppt-1 ~> m ppt-1].
1500  real, intent(in) :: dT_to_dColHt_b !< A factor (mass_lay*dSColHtc_vol/dT) relating
1501  !! a layer's temperature change to the change in column
1502  !! height, including all implicit diffusive changes
1503  !! in the temperatures of all the layers below [Z degC-1 ~> m degC-1].
1504  real, intent(in) :: dS_to_dColHt_b !< A factor (mass_lay*dSColHtc_vol/dS) relating
1505  !! a layer's salinity change to the change in column
1506  !! height, including all implicit diffusive changes
1507  !! in the salinities of all the layers below [Z ppt-1 ~> m ppt-1].
1508 
1509  real, optional, intent(out) :: PE_chg !< The change in column potential energy from applying
1510  !! Kddt_h at the present interface [R Z3 T-2 ~> J m-2].
1511  real, optional, intent(out) :: dPEc_dKd !< The partial derivative of PE_chg with Kddt_h
1512  !! [R Z3 T-2 H-1 ~> J m-3 or J kg-1].
1513  real, optional, intent(out) :: dPE_max !< The maximum change in column potential energy that could
1514  !! be realizedd by applying a huge value of Kddt_h at the
1515  !! present interface [R Z3 T-2 ~> J m-2].
1516  real, optional, intent(out) :: dPEc_dKd_0 !< The partial derivative of PE_chg with Kddt_h in the
1517  !! limit where Kddt_h = 0 [R Z3 T-2 H-1 ~> J m-3 or J kg-1].
1518  real, optional, intent(out) :: PE_ColHt_cor !< The correction to PE_chg that is made due to a net
1519  !! change in the column height [R Z3 T-2 ~> J m-2].
1520 
1521  real :: hps ! The sum of the two effective pivot thicknesses [H ~> m or kg m-2].
1522  real :: bdt1 ! A product of the two pivot thicknesses plus a diffusive term [H2 ~> m2 or kg2 m-4].
1523  real :: dT_c ! The core term in the expressions for the temperature changes [degC H2 ~> degC m2 or degC kg2 m-4].
1524  real :: dS_c ! The core term in the expressions for the salinity changes [ppt H2 ~> ppt m2 or ppt kg2 m-4].
1525  real :: PEc_core ! The diffusivity-independent core term in the expressions
1526  ! for the potential energy changes [R Z2 T-2 ~> J m-3].
1527  real :: ColHt_core ! The diffusivity-independent core term in the expressions
1528  ! for the column height changes [H Z ~> m2 or kg m-1].
1529  real :: ColHt_chg ! The change in the column height [H ~> m or kg m-2].
1530  real :: y1_3 ! A local temporary term in [H-3 ~> m-3 or m6 kg-3].
1531  real :: y1_4 ! A local temporary term in [H-4 ~> m-4 or m8 kg-4].
1532 
1533  ! The expression for the change in potential energy used here is derived
1534  ! from the expression for the final estimates of the changes in temperature
1535  ! and salinities, and then extensively manipulated to get it into its most
1536  ! succint form. The derivation is not necessarily obvious, but it demonstrably
1537  ! works by comparison with separate calculations of the energy changes after
1538  ! the tridiagonal solver for the final changes in temperature and salinity are
1539  ! applied.
1540 
1541  hps = hp_a + hp_b
1542  bdt1 = hp_a * hp_b + kddt_h0 * hps
1543  dt_c = hp_a * th_b - hp_b * th_a
1544  ds_c = hp_a * sh_b - hp_b * sh_a
1545  pec_core = hp_b * (dt_to_dpe_a * dt_c + ds_to_dpe_a * ds_c) - &
1546  hp_a * (dt_to_dpe_b * dt_c + ds_to_dpe_b * ds_c)
1547  colht_core = hp_b * (dt_to_dcolht_a * dt_c + ds_to_dcolht_a * ds_c) - &
1548  hp_a * (dt_to_dcolht_b * dt_c + ds_to_dcolht_b * ds_c)
1549 
1550  if (present(pe_chg)) then
1551  ! Find the change in column potential energy due to the change in the
1552  ! diffusivity at this interface by dKddt_h.
1553  y1_3 = dkddt_h / (bdt1 * (bdt1 + dkddt_h * hps))
1554  pe_chg = pec_core * y1_3
1555  colht_chg = colht_core * y1_3
1556  if (colht_chg < 0.0) pe_chg = pe_chg - pres_z * colht_chg
1557  if (present(pe_colht_cor)) pe_colht_cor = -pres_z * min(colht_chg, 0.0)
1558  elseif (present(pe_colht_cor)) then
1559  y1_3 = dkddt_h / (bdt1 * (bdt1 + dkddt_h * hps))
1560  pe_colht_cor = -pres_z * min(colht_core * y1_3, 0.0)
1561  endif
1562 
1563  if (present(dpec_dkd)) then
1564  ! Find the derivative of the potential energy change with dKddt_h.
1565  y1_4 = 1.0 / (bdt1 + dkddt_h * hps)**2
1566  dpec_dkd = pec_core * y1_4
1567  colht_chg = colht_core * y1_4
1568  if (colht_chg < 0.0) dpec_dkd = dpec_dkd - pres_z * colht_chg
1569  endif
1570 
1571  if (present(dpe_max)) then
1572  ! This expression is the limit of PE_chg for infinite dKddt_h.
1573  y1_3 = 1.0 / (bdt1 * hps)
1574  dpe_max = pec_core * y1_3
1575  colht_chg = colht_core * y1_3
1576  if (colht_chg < 0.0) dpe_max = dpe_max - pres_z * colht_chg
1577  endif
1578 
1579  if (present(dpec_dkd_0)) then
1580  ! This expression is the limit of dPEc_dKd for dKddt_h = 0.
1581  y1_4 = 1.0 / bdt1**2
1582  dpec_dkd_0 = pec_core * y1_4
1583  colht_chg = colht_core * y1_4
1584  if (colht_chg < 0.0) dpec_dkd_0 = dpec_dkd_0 - pres_z * colht_chg
1585  endif
1586 

Referenced by epbl_column().

Here is the caller graph for this function:

◆ find_pe_chg_orig()

subroutine mom_energetic_pbl::find_pe_chg_orig ( real, intent(in)  Kddt_h,
real, intent(in)  h_k,
real, intent(in)  b_den_1,
real, intent(in)  dTe_term,
real, intent(in)  dSe_term,
real, intent(in)  dT_km1_t2,
real, intent(in)  dS_km1_t2,
real, intent(in)  dT_to_dPE_k,
real, intent(in)  dS_to_dPE_k,
real, intent(in)  dT_to_dPEa,
real, intent(in)  dS_to_dPEa,
real, intent(in)  pres_Z,
real, intent(in)  dT_to_dColHt_k,
real, intent(in)  dS_to_dColHt_k,
real, intent(in)  dT_to_dColHta,
real, intent(in)  dS_to_dColHta,
real, intent(out), optional  PE_chg,
real, intent(out), optional  dPEc_dKd,
real, intent(out), optional  dPE_max,
real, intent(out), optional  dPEc_dKd_0 
)
private

This subroutine calculates the change in potential energy and or derivatives for several changes in an interfaces's diapycnal diffusivity times a timestep using the original form used in the first version of ePBL.

Parameters
[in]kddt_hThe diffusivity at an interface times the time step and divided by the average of the thicknesses around the interface [H ~> m or kg m-2].
[in]h_kThe thickness of the layer below the interface [H ~> m or kg m-2].
[in]b_den_1The first term in the denominator of the pivot for the tridiagonal solver, given by h_k plus a term that is a fraction (determined from the tridiagonal solver) of Kddt_h for the interface above [H ~> m or kg m-2].
[in]dte_termA diffusivity-independent term related to the temperature change in the layer below the interface [degC H ~> degC m or degC kg m-2].
[in]dse_termA diffusivity-independent term related to the salinity change in the layer below the interface [ppt H ~> ppt m or ppt kg m-2].
[in]dt_km1_t2A diffusivity-independent term related to the temperature change in the layer above the interface [degC].
[in]ds_km1_t2A diffusivity-independent term related to the salinity change in the layer above the interface [ppt].
[in]pres_zThe rescaled hydrostatic interface pressure, which relates the changes in column thickness to the energy that is radiated as gravity waves and unavailable to drive mixing [R Z2 T-2 ~> J m-3].
[in]dt_to_dpe_kA factor (pres_lay*mass_lay*dSpec_vol/dT) relating a layer's temperature change to the change in column potential energy, including all implicit diffusive changes in the temperatures of all the layers below [R Z3 T-2 degC-1 ~> J m-2 degC-1].
[in]ds_to_dpe_kA factor (pres_lay*mass_lay*dSpec_vol/dS) relating a layer's salinity change to the change in column potential energy, including all implicit diffusive changes in the in the salinities of all the layers below [R Z3 T-2 ppt-1 ~> J m-2 ppt-1].
[in]dt_to_dpeaA factor (pres_lay*mass_lay*dSpec_vol/dT) relating a layer's temperature change to the change in column potential energy, including all implicit diffusive changes in the temperatures of all the layers above [R Z3 T-2 degC-1 ~> J m-2 degC-1].
[in]ds_to_dpeaA factor (pres_lay*mass_lay*dSpec_vol/dS) relating a layer's salinity change to the change in column potential energy, including all implicit diffusive changes in the salinities of all the layers above [R Z3 T-2 ppt-1 ~> J m-2 ppt-1].
[in]dt_to_dcolht_kA factor (mass_lay*dSColHtc_vol/dT) relating a layer's temperature change to the change in column height, including all implicit diffusive changes in the temperatures of all the layers below [Z degC-1 ~> m degC-1].
[in]ds_to_dcolht_kA factor (mass_lay*dSColHtc_vol/dS) relating a layer's salinity change to the change in column height, including all implicit diffusive changes in the salinities of all the layers below [Z ppt-1 ~> m ppt-1].
[in]dt_to_dcolhtaA factor (mass_lay*dSColHtc_vol/dT) relating a layer's temperature change to the change in column height, including all implicit diffusive changes in the temperatures of all the layers above [Z degC-1 ~> m degC-1].
[in]ds_to_dcolhtaA factor (mass_lay*dSColHtc_vol/dS) relating a layer's salinity change to the change in column height, including all implicit diffusive changes in the salinities of all the layers above [Z ppt-1 ~> m ppt-1].
[out]pe_chgThe change in column potential energy from applying Kddt_h at the present interface [R Z3 T-2 ~> J m-2].
[out]dpec_dkdThe partial derivative of PE_chg with Kddt_h [R Z3 T-2 H-1 ~> J m-3 or J kg-1].
[out]dpe_maxThe maximum change in column potential energy that could be realizedd by applying a huge value of Kddt_h at the present interface [R Z3 T-2 ~> J m-2].
[out]dpec_dkd_0The partial derivative of PE_chg with Kddt_h in the limit where Kddt_h = 0 [R Z3 T-2 H-1 ~> J m-3 or J kg-1].

Definition at line 1597 of file MOM_energetic_PBL.F90.

1597  real, intent(in) :: Kddt_h !< The diffusivity at an interface times the time step and
1598  !! divided by the average of the thicknesses around the
1599  !! interface [H ~> m or kg m-2].
1600  real, intent(in) :: h_k !< The thickness of the layer below the interface [H ~> m or kg m-2].
1601  real, intent(in) :: b_den_1 !< The first term in the denominator of the pivot
1602  !! for the tridiagonal solver, given by h_k plus a term that
1603  !! is a fraction (determined from the tridiagonal solver) of
1604  !! Kddt_h for the interface above [H ~> m or kg m-2].
1605  real, intent(in) :: dTe_term !< A diffusivity-independent term related to the temperature change
1606  !! in the layer below the interface [degC H ~> degC m or degC kg m-2].
1607  real, intent(in) :: dSe_term !< A diffusivity-independent term related to the salinity change
1608  !! in the layer below the interface [ppt H ~> ppt m or ppt kg m-2].
1609  real, intent(in) :: dT_km1_t2 !< A diffusivity-independent term related to the
1610  !! temperature change in the layer above the interface [degC].
1611  real, intent(in) :: dS_km1_t2 !< A diffusivity-independent term related to the
1612  !! salinity change in the layer above the interface [ppt].
1613  real, intent(in) :: pres_Z !< The rescaled hydrostatic interface pressure, which relates
1614  !! the changes in column thickness to the energy that is radiated
1615  !! as gravity waves and unavailable to drive mixing [R Z2 T-2 ~> J m-3].
1616  real, intent(in) :: dT_to_dPE_k !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating
1617  !! a layer's temperature change to the change in column potential
1618  !! energy, including all implicit diffusive changes in the
1619  !! temperatures of all the layers below [R Z3 T-2 degC-1 ~> J m-2 degC-1].
1620  real, intent(in) :: dS_to_dPE_k !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating
1621  !! a layer's salinity change to the change in column potential
1622  !! energy, including all implicit diffusive changes in the
1623  !! in the salinities of all the layers below [R Z3 T-2 ppt-1 ~> J m-2 ppt-1].
1624  real, intent(in) :: dT_to_dPEa !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating
1625  !! a layer's temperature change to the change in column potential
1626  !! energy, including all implicit diffusive changes in the
1627  !! temperatures of all the layers above [R Z3 T-2 degC-1 ~> J m-2 degC-1].
1628  real, intent(in) :: dS_to_dPEa !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating
1629  !! a layer's salinity change to the change in column potential
1630  !! energy, including all implicit diffusive changes in the
1631  !! salinities of all the layers above [R Z3 T-2 ppt-1 ~> J m-2 ppt-1].
1632  real, intent(in) :: dT_to_dColHt_k !< A factor (mass_lay*dSColHtc_vol/dT) relating
1633  !! a layer's temperature change to the change in column
1634  !! height, including all implicit diffusive changes in the
1635  !! temperatures of all the layers below [Z degC-1 ~> m degC-1].
1636  real, intent(in) :: dS_to_dColHt_k !< A factor (mass_lay*dSColHtc_vol/dS) relating
1637  !! a layer's salinity change to the change in column
1638  !! height, including all implicit diffusive changes
1639  !! in the salinities of all the layers below [Z ppt-1 ~> m ppt-1].
1640  real, intent(in) :: dT_to_dColHta !< A factor (mass_lay*dSColHtc_vol/dT) relating
1641  !! a layer's temperature change to the change in column
1642  !! height, including all implicit diffusive changes
1643  !! in the temperatures of all the layers above [Z degC-1 ~> m degC-1].
1644  real, intent(in) :: dS_to_dColHta !< A factor (mass_lay*dSColHtc_vol/dS) relating
1645  !! a layer's salinity change to the change in column
1646  !! height, including all implicit diffusive changes
1647  !! in the salinities of all the layers above [Z ppt-1 ~> m ppt-1].
1648 
1649  real, optional, intent(out) :: PE_chg !< The change in column potential energy from applying
1650  !! Kddt_h at the present interface [R Z3 T-2 ~> J m-2].
1651  real, optional, intent(out) :: dPEc_dKd !< The partial derivative of PE_chg with Kddt_h
1652  !! [R Z3 T-2 H-1 ~> J m-3 or J kg-1].
1653  real, optional, intent(out) :: dPE_max !< The maximum change in column potential energy that could
1654  !! be realizedd by applying a huge value of Kddt_h at the
1655  !! present interface [R Z3 T-2 ~> J m-2].
1656  real, optional, intent(out) :: dPEc_dKd_0 !< The partial derivative of PE_chg with Kddt_h in the
1657  !! limit where Kddt_h = 0 [R Z3 T-2 H-1 ~> J m-3 or J kg-1].
1658 
1659 ! This subroutine determines the total potential energy change due to mixing
1660 ! at an interface, including all of the implicit effects of the prescribed
1661 ! mixing at interfaces above. Everything here is derived by careful manipulation
1662 ! of the robust tridiagonal solvers used for tracers by MOM6. The results are
1663 ! positive for mixing in a stably stratified environment.
1664 ! The comments describing these arguments are for a downward mixing pass, but
1665 ! this routine can also be used for an upward pass with the sense of direction
1666 ! reversed.
1667 
1668  real :: b1 ! b1 is used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
1669  real :: b1Kd ! Temporary array [nondim]
1670  real :: ColHt_chg ! The change in column thickness [Z ~> m].
1671  real :: dColHt_max ! The change in column thickness for infinite diffusivity [Z ~> m].
1672  real :: dColHt_dKd ! The partial derivative of column thickness with diffusivity [s Z-1 ~> s m-1].
1673  real :: dT_k, dT_km1 ! Temporary arrays [degC].
1674  real :: dS_k, dS_km1 ! Temporary arrays [ppt].
1675  real :: I_Kr_denom ! Temporary array [H-2 ~> m-2 or m4 kg-2]
1676  real :: dKr_dKd ! Nondimensional temporary array.
1677  real :: ddT_k_dKd, ddT_km1_dKd ! Temporary arrays [degC H-1 ~> m-1 or m2 kg-1].
1678  real :: ddS_k_dKd, ddS_km1_dKd ! Temporary arrays [ppt H-1 ~> ppt m-1 or ppt m2 kg-1].
1679 
1680  b1 = 1.0 / (b_den_1 + kddt_h)
1681  b1kd = kddt_h*b1
1682 
1683  ! Start with the temperature change in layer k-1 due to the diffusivity at
1684  ! interface K without considering the effects of changes in layer k.
1685 
1686  ! Calculate the change in PE due to the diffusion at interface K
1687  ! if Kddt_h(K+1) = 0.
1688  i_kr_denom = 1.0 / (h_k*b_den_1 + (b_den_1 + h_k)*kddt_h)
1689 
1690  dt_k = (kddt_h*i_kr_denom) * dte_term
1691  ds_k = (kddt_h*i_kr_denom) * dse_term
1692 
1693  if (present(pe_chg)) then
1694  ! Find the change in energy due to diffusion with strength Kddt_h at this interface.
1695  ! Increment the temperature changes in layer k-1 due the changes in layer k.
1696  dt_km1 = b1kd * ( dt_k + dt_km1_t2 )
1697  ds_km1 = b1kd * ( ds_k + ds_km1_t2 )
1698  pe_chg = (dt_to_dpe_k * dt_k + dt_to_dpea * dt_km1) + &
1699  (ds_to_dpe_k * ds_k + ds_to_dpea * ds_km1)
1700  colht_chg = (dt_to_dcolht_k * dt_k + dt_to_dcolhta * dt_km1) + &
1701  (ds_to_dcolht_k * ds_k + ds_to_dcolhta * ds_km1)
1702  if (colht_chg < 0.0) pe_chg = pe_chg - pres_z * colht_chg
1703  endif
1704 
1705  if (present(dpec_dkd)) then
1706  ! Find the derivatives of the temperature and salinity changes with Kddt_h.
1707  dkr_dkd = (h_k*b_den_1) * i_kr_denom**2
1708 
1709  ddt_k_dkd = dkr_dkd * dte_term
1710  dds_k_dkd = dkr_dkd * dse_term
1711  ddt_km1_dkd = (b1**2 * b_den_1) * ( dt_k + dt_km1_t2 ) + b1kd * ddt_k_dkd
1712  dds_km1_dkd = (b1**2 * b_den_1) * ( ds_k + ds_km1_t2 ) + b1kd * dds_k_dkd
1713 
1714  ! Calculate the partial derivative of Pe_chg with Kddt_h.
1715  dpec_dkd = (dt_to_dpe_k * ddt_k_dkd + dt_to_dpea * ddt_km1_dkd) + &
1716  (ds_to_dpe_k * dds_k_dkd + ds_to_dpea * dds_km1_dkd)
1717  dcolht_dkd = (dt_to_dcolht_k * ddt_k_dkd + dt_to_dcolhta * ddt_km1_dkd) + &
1718  (ds_to_dcolht_k * dds_k_dkd + ds_to_dcolhta * dds_km1_dkd)
1719  if (dcolht_dkd < 0.0) dpec_dkd = dpec_dkd - pres_z * dcolht_dkd
1720  endif
1721 
1722  if (present(dpe_max)) then
1723  ! This expression is the limit of PE_chg for infinite Kddt_h.
1724  dpe_max = (dt_to_dpea * dt_km1_t2 + ds_to_dpea * ds_km1_t2) + &
1725  ((dt_to_dpe_k + dt_to_dpea) * dte_term + &
1726  (ds_to_dpe_k + ds_to_dpea) * dse_term) / (b_den_1 + h_k)
1727  dcolht_max = (dt_to_dcolhta * dt_km1_t2 + ds_to_dcolhta * ds_km1_t2) + &
1728  ((dt_to_dcolht_k + dt_to_dcolhta) * dte_term + &
1729  (ds_to_dcolht_k + ds_to_dcolhta) * dse_term) / (b_den_1 + h_k)
1730  if (dcolht_max < 0.0) dpe_max = dpe_max - pres_z*dcolht_max
1731  endif
1732 
1733  if (present(dpec_dkd_0)) then
1734  ! This expression is the limit of dPEc_dKd for Kddt_h = 0.
1735  dpec_dkd_0 = (dt_to_dpea * dt_km1_t2 + ds_to_dpea * ds_km1_t2) / (b_den_1) + &
1736  (dt_to_dpe_k * dte_term + ds_to_dpe_k * dse_term) / (h_k*b_den_1)
1737  dcolht_dkd = (dt_to_dcolhta * dt_km1_t2 + ds_to_dcolhta * ds_km1_t2) / (b_den_1) + &
1738  (dt_to_dcolht_k * dte_term + ds_to_dcolht_k * dse_term) / (h_k*b_den_1)
1739  if (dcolht_dkd < 0.0) dpec_dkd_0 = dpec_dkd_0 - pres_z*dcolht_dkd
1740  endif
1741 

Referenced by epbl_column().

Here is the caller graph for this function:

◆ mstar_langmuir()

subroutine mom_energetic_pbl::mstar_langmuir ( type(energetic_pbl_cs), pointer  CS,
type(unit_scale_type), intent(in)  US,
real, intent(in)  Abs_Coriolis,
real, intent(in)  Buoyancy_Flux,
real, intent(in)  UStar,
real, intent(in)  BLD,
real, intent(in)  Langmuir_Number,
real, intent(inout)  Mstar,
real, intent(out)  MStar_LT,
real, intent(out)  Convect_Langmuir_Number 
)
private

This subroutine modifies the Mstar value if the Langmuir number is present.

Parameters
csEnergetic_PBL control structure.
[in]usA dimensional unit scaling type
[in]abs_coriolisAbsolute value of the Coriolis parameter [T-1 ~> s-1]
[in]buoyancy_fluxBuoyancy flux [Z2 T-3 ~> m2 s-3]
[in]ustarSurface friction velocity with? gustiness [Z T-1 ~> m s-1]
[in]bldboundary layer depth [Z ~> m]
[in,out]mstarInput/output mstar (Mixing/ustar**3) [nondim]
[in]langmuir_numberLangmuir number [nondim]
[out]mstar_ltMstar increase due to Langmuir turbulence [nondim]
[out]convect_langmuir_numberLangmuir number including buoyancy flux [nondim]

Definition at line 1834 of file MOM_energetic_PBL.F90.

1834  type(energetic_PBL_CS), pointer :: CS !< Energetic_PBL control structure.
1835  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1836  real, intent(in) :: Abs_Coriolis !< Absolute value of the Coriolis parameter [T-1 ~> s-1]
1837  real, intent(in) :: Buoyancy_Flux !< Buoyancy flux [Z2 T-3 ~> m2 s-3]
1838  real, intent(in) :: UStar !< Surface friction velocity with? gustiness [Z T-1 ~> m s-1]
1839  real, intent(in) :: BLD !< boundary layer depth [Z ~> m]
1840  real, intent(inout) :: Mstar !< Input/output mstar (Mixing/ustar**3) [nondim]
1841  real, intent(in) :: Langmuir_Number !< Langmuir number [nondim]
1842  real, intent(out) :: MStar_LT !< Mstar increase due to Langmuir turbulence [nondim]
1843  real, intent(out) :: Convect_Langmuir_number !< Langmuir number including buoyancy flux [nondim]
1844 
1845  !/
1846  real, parameter :: Max_ratio = 1.0e16 ! The maximum value of a nondimensional ratio.
1847  real :: enhance_mstar ! A multiplicative scaling of mstar due to Langmuir turbulence.
1848  real :: mstar_LT_add ! A value that is added to mstar due to Langmuir turbulence.
1849  real :: iL_Ekman ! Inverse of Ekman length scale [Z-1 ~> m-1].
1850  real :: iL_Obukhov ! Inverse of Obukhov length scale [Z-1 ~> m-1].
1851  real :: I_ustar ! The Adcroft reciprocal of ustar [T Z-1 ~> s m-1]
1852  real :: I_f ! The Adcroft reciprocal of the Coriolis parameter [T ~> s]
1853  real :: MLD_Ekman ! The ratio of the mixed layer depth to the Ekman layer depth [nondim].
1854  real :: Ekman_Obukhov ! The Ekman layer thickness divided by the Obukhov depth [nondim].
1855  real :: MLD_Obukhov ! The mixed layer depth divided by the Obukhov depth [nondim].
1856  real :: MLD_Obukhov_stab ! Ratios of length scales where MLD is boundary layer depth [nondim].
1857  real :: Ekman_Obukhov_stab ! >
1858  real :: MLD_Obukhov_un ! Ratios of length scales where MLD is boundary layer depth
1859  real :: Ekman_Obukhov_un ! >
1860 
1861  ! Set default values for no Langmuir effects.
1862  enhance_mstar = 1.0 ; mstar_lt_add = 0.0
1863 
1864  if (cs%LT_Enhance_Form /= no_langmuir) then
1865  ! a. Get parameters for modified LA
1866  if (cs%answers_2018) then
1867  il_ekman = abs_coriolis / ustar
1868  il_obukhov = buoyancy_flux*cs%vonkar / ustar**3
1869  ekman_obukhov_stab = abs(max(0., il_obukhov / (il_ekman + 1.e-10*us%Z_to_m)))
1870  ekman_obukhov_un = abs(min(0., il_obukhov / (il_ekman + 1.e-10*us%Z_to_m)))
1871  mld_obukhov_stab = abs(max(0., bld*il_obukhov))
1872  mld_obukhov_un = abs(min(0., bld*il_obukhov))
1873  mld_ekman = abs( bld*il_ekman )
1874  else
1875  ekman_obukhov = max_ratio ; mld_obukhov = max_ratio ; mld_ekman = max_ratio
1876  i_f = 0.0 ; if (abs(abs_coriolis) > 0.0) i_f = 1.0 / abs_coriolis
1877  i_ustar = 0.0 ; if (abs(ustar) > 0.0) i_ustar = 1.0 / ustar
1878  if (abs(buoyancy_flux*cs%vonkar) < max_ratio*(abs_coriolis * ustar**2)) &
1879  ekman_obukhov = abs(buoyancy_flux*cs%vonkar) * (i_f * i_ustar**2)
1880  if (abs(bld*buoyancy_flux*cs%vonkar) < max_ratio*ustar**3) &
1881  mld_obukhov = abs(bld*buoyancy_flux*cs%vonkar) * i_ustar**3
1882  if (bld*abs_coriolis < max_ratio*ustar) &
1883  mld_ekman = bld*abs_coriolis * i_ustar
1884 
1885  if (buoyancy_flux > 0.0) then
1886  ekman_obukhov_stab = ekman_obukhov ; ekman_obukhov_un = 0.0
1887  mld_obukhov_stab = mld_obukhov ; mld_obukhov_un = 0.0
1888  else
1889  ekman_obukhov_un = ekman_obukhov ; ekman_obukhov_stab = 0.0
1890  mld_obukhov_un = mld_obukhov ; mld_obukhov_stab = 0.0
1891  endif
1892  endif
1893 
1894  ! b. Adjust LA based on various parameters.
1895  ! Assumes linear factors based on length scale ratios to adjust LA
1896  ! Note when these coefficients are set to 0 recovers simple LA.
1897  convect_langmuir_number = langmuir_number * &
1898  ( (1.0 + max(-0.5, cs%LaC_MLDoEK * mld_ekman)) + &
1899  ((cs%LaC_EKoOB_stab * ekman_obukhov_stab + cs%LaC_EKoOB_un * ekman_obukhov_un) + &
1900  (cs%LaC_MLDoOB_stab * mld_obukhov_stab + cs%LaC_MLDoOB_un * mld_obukhov_un)) )
1901 
1902  if (cs%LT_Enhance_Form == langmuir_rescale) then
1903  ! Enhancement is multiplied (added mst_lt set to 0)
1904  enhance_mstar = min(cs%Max_Enhance_M, &
1905  (1. + cs%LT_ENHANCE_COEF * convect_langmuir_number**cs%LT_ENHANCE_EXP) )
1906  elseif (cs%LT_ENHANCE_Form == langmuir_add) then
1907  ! or Enhancement is additive (multiplied enhance_m set to 1)
1908  mstar_lt_add = cs%LT_ENHANCE_COEF * convect_langmuir_number**cs%LT_ENHANCE_EXP
1909  endif
1910  endif
1911 
1912  mstar_lt = (enhance_mstar - 1.0)*mstar + mstar_lt_add ! Diagnose the full increase in mstar.
1913  mstar = mstar*enhance_mstar + mstar_lt_add
1914 

References langmuir_add, langmuir_rescale, and no_langmuir.

Referenced by find_mstar().

Here is the caller graph for this function:

Variable Documentation

◆ additive_string

character*(20), parameter mom_energetic_pbl::additive_string = "ADDITIVE"
private

Enumeration values for mstar_Scheme.

The value of mstar_scheme to use a constant mstar

Definition at line 217 of file MOM_energetic_PBL.F90.

217 character*(20), parameter :: ADDITIVE_STRING = "ADDITIVE"

Referenced by energetic_pbl_init().

◆ constant_string

character*(20), parameter mom_energetic_pbl::constant_string = "CONSTANT"
private

Enumeration values for mstar_Scheme.

The value of mstar_scheme to use a constant mstar

Definition at line 211 of file MOM_energetic_PBL.F90.

211 character*(20), parameter :: CONSTANT_STRING = "CONSTANT"

Referenced by energetic_pbl_init().

◆ langmuir_add

integer, parameter mom_energetic_pbl::langmuir_add = 3
private

The value of LT_ENHANCE_FORM to add a contribution to mstar from Langmuir turblence to other contributions.

Definition at line 204 of file MOM_energetic_PBL.F90.

204 integer, parameter :: Langmuir_add = 3 !< The value of LT_ENHANCE_FORM to add a contribution to

Referenced by energetic_pbl_init(), and mstar_langmuir().

◆ langmuir_rescale

integer, parameter mom_energetic_pbl::langmuir_rescale = 2
private

The value of LT_ENHANCE_FORM to use a multiplicative rescaling of mstar to account for Langmuir turbulence.

Definition at line 202 of file MOM_energetic_PBL.F90.

202 integer, parameter :: Langmuir_rescale = 2 !< The value of LT_ENHANCE_FORM to use a multiplicative

Referenced by energetic_pbl_init(), and mstar_langmuir().

◆ mstar_from_ekman

integer, parameter mom_energetic_pbl::mstar_from_ekman = 2
private

The value of mstar_scheme to base mstar on the ratio of the Ekman layer depth to the Obukhov depth.

Definition at line 198 of file MOM_energetic_PBL.F90.

198 integer, parameter :: MStar_from_Ekman = 2 !< The value of mstar_scheme to base mstar on the ratio

Referenced by energetic_pbl_init(), and find_mstar().

◆ mstar_from_rh18

integer, parameter mom_energetic_pbl::mstar_from_rh18 = 3
private

The value of mstar_scheme to base mstar of of RH18.

Definition at line 200 of file MOM_energetic_PBL.F90.

200 integer, parameter :: MStar_from_RH18 = 3 !< The value of mstar_scheme to base mstar of of RH18

Referenced by energetic_pbl_init(), and find_mstar().

◆ no_langmuir

integer, parameter mom_energetic_pbl::no_langmuir = 0
private

The value of LT_ENHANCE_FORM not use Langmuir turbolence.

Definition at line 201 of file MOM_energetic_PBL.F90.

201 integer, parameter :: No_Langmuir = 0 !< The value of LT_ENHANCE_FORM not use Langmuir turbolence.

Referenced by energetic_pbl_init(), and mstar_langmuir().

◆ none_string

character*(20), parameter mom_energetic_pbl::none_string = "NONE"
private

Enumeration values for mstar_Scheme.

The value of mstar_scheme to use a constant mstar

Definition at line 215 of file MOM_energetic_PBL.F90.

215 character*(20), parameter :: NONE_STRING = "NONE"

Referenced by energetic_pbl_init().

◆ om4_string

character*(20), parameter mom_energetic_pbl::om4_string = "OM4"
private

Enumeration values for mstar_Scheme.

The value of mstar_scheme to use a constant mstar

Definition at line 212 of file MOM_energetic_PBL.F90.

212 character*(20), parameter :: OM4_STRING = "OM4"

Referenced by energetic_pbl_init().

◆ rescaled_string

character*(20), parameter mom_energetic_pbl::rescaled_string = "RESCALE"
private

Enumeration values for mstar_Scheme.

The value of mstar_scheme to use a constant mstar

Definition at line 216 of file MOM_energetic_PBL.F90.

216 character*(20), parameter :: RESCALED_STRING = "RESCALE"

Referenced by energetic_pbl_init().

◆ rh18_string

character*(20), parameter mom_energetic_pbl::rh18_string = "REICHL_H18"
private

Enumeration values for mstar_Scheme.

The value of mstar_scheme to use a constant mstar

Definition at line 213 of file MOM_energetic_PBL.F90.

213 character*(20), parameter :: RH18_STRING = "REICHL_H18"

Referenced by energetic_pbl_init().

◆ root_tke_string

character*(20), parameter mom_energetic_pbl::root_tke_string = "CUBE_ROOT_TKE"
private

Enumeration values for mstar_Scheme.

The value of mstar_scheme to use a constant mstar

Definition at line 214 of file MOM_energetic_PBL.F90.

214 character*(20), parameter :: ROOT_TKE_STRING = "CUBE_ROOT_TKE"

Referenced by energetic_pbl_init().

◆ use_fixed_mstar

integer, parameter mom_energetic_pbl::use_fixed_mstar = 0

Enumeration values for mstar_Scheme.

The value of mstar_scheme to use a constant mstar

Definition at line 197 of file MOM_energetic_PBL.F90.

197 integer, parameter :: Use_Fixed_MStar = 0 !< The value of mstar_scheme to use a constant mstar

Referenced by energetic_pbl_init(), epbl_column(), and find_mstar().

◆ wt_from_croot_tke

integer, parameter mom_energetic_pbl::wt_from_croot_tke = 0
private

Use a constant times the cube root of remaining TKE to calculate the turbulent velocity.

Definition at line 206 of file MOM_energetic_PBL.F90.

206 integer, parameter :: wT_from_cRoot_TKE = 0 !< Use a constant times the cube root of remaining TKE

Referenced by energetic_pbl_init(), and epbl_column().

◆ wt_from_rh18

integer, parameter mom_energetic_pbl::wt_from_rh18 = 1
private

Use a scheme based on a combination of w* and v* as documented in Reichl & Hallberg (2018) to calculate the turbulent velocity.

Definition at line 208 of file MOM_energetic_PBL.F90.

208 integer, parameter :: wT_from_RH18 = 1 !< Use a scheme based on a combination of w* and v* as

Referenced by energetic_pbl_init(), and epbl_column().