mom_energetic_pbl Module Reference

Data Types
type	energetic_pbl_cs

Functions/Subroutines
subroutine, public	energetic_pbl (h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, dSV_dT, dSV_dS, TKE_forced, Buoy_Flux, dt_diag, last_call, dT_expected, dS_expected)
	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	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, 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, 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, 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, public	energetic_pbl_get_mld (CS, MLD, G)
	Copies the ePBL active mixed layer depth into MLD. More...
subroutine	ust_2_u10_coare3p5 (USTair, U10, GV)
	Computes wind speed from ustar_air based on COARE 3.5 Cd relationship. More...
subroutine	get_la_windsea (ustar, hbl, GV, LA)
subroutine, public	energetic_pbl_init (Time, G, GV, param_file, diag, CS)
subroutine, public	energetic_pbl_end (CS)

Variables
integer	num_msg = 0
integer	max_msg = 2

Function/Subroutine Documentation

energetic_pbl()

subroutine, public mom_energetic_pbl::energetic_pbl	(	real, dimension(szi_(g),szj_(g),szk_(gv)), intent(inout)	h_3d,
		real, dimension(szi_(g),szj_(g),szk_(gv)), intent(in)	u_3d,
		real, dimension(szi_(g),szj_(g),szk_(gv)), intent(in)	v_3d,
		type(thermo_var_ptrs), intent(inout)	tv,
		type(forcing), intent(inout)	fluxes,
		real, intent(in)	dt,
		real, dimension(szi_(g),szj_(g),szk_(gv)+1), intent(out)	Kd_int,
		type(ocean_grid_type), intent(inout)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(energetic_pbl_cs), pointer	CS,
		real, dimension(szi_(g),szj_(g),szk_(gv)), intent(in)	dSV_dT,
		real, dimension(szi_(g),szj_(g),szk_(gv)), intent(in)	dSV_dS,
		real, dimension(szi_(g),szj_(g),szk_(gv)), intent(in)	TKE_forced,
		real, dimension(szi_(g),szj_(g)), intent(in)	Buoy_Flux,
		real, intent(in), optional	dt_diag,
		logical, intent(in), optional	last_call,
		real, dimension(szi_(g),szj_(g),szk_(gv)), intent(out), optional	dT_expected,
		real, dimension(szi_(g),szj_(g),szk_(gv)), intent(out), optional	dS_expected
	)

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]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in,out]	h_3d	Layer thickness, in m or kg m-2.
[in]	u_3d	Zonal velocities interpolated to h points,
[in]	v_3d	Zonal velocities interpolated to h points,
[in]	dsv_dt	The partial derivative of in-situ specific
[in]	dsv_ds	The partial derivative of in-situ specific
[in]	tke_forced	The forcing requirements to homogenize the
[in,out]	tv	A structure containing pointers to any available thermodynamic fields. Absent fields have NULL ptrs.
[in,out]	fluxes	A structure containing pointers to any possible forcing fields. Unused fields have NULL ptrs.
[in]	dt	Time increment, in s.
[out]	kd_int	The diagnosed diffusivities at interfaces,
	cs	The control structure returned by a previous call to mixedlayer_init.
[in]	buoy_flux	The surface buoyancy flux. in m2/s3.
[in]	dt_diag	The diagnostic time step, which may be less than dt if there are two callse to mixedlayer, in s.
[in]	last_call	If true, this is the last call to mixedlayer in the current time step, so diagnostics will be written. The default is .true.

Definition at line 237 of file MOM_energetic_PBL.F90.

References find_pe_chg(), find_pe_chg_orig(), get_la_windsea(), and mom_error_handler::mom_error().

  type(ocean_grid_type),   intent(inout) :: g      !< The ocean's grid structure.
  type(verticalgrid_type), intent(in)    :: gv     !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: h_3d   !< Layer thickness, in m or kg m-2.
                                                   !! (Intent in/out) The units of h are referred
                                                   !! to as H below.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: u_3d   !< Zonal velocities interpolated to h points,
                                                   !! m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: v_3d   !< Zonal velocities interpolated to h points,
                                                   !! m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: dsv_dt !< The partial derivative of in-situ specific
                                                   !! volume with potential temperature,
                                                   !! in m3 kg-1 K-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: dsv_ds !< The partial derivative of in-situ specific
                                                   !! volume with salinity, in m3 kg-1 ppt-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: tke_forced !< The forcing requirements to homogenize the
                                                   !! forcing that has been applied to each layer
                                                   !! through each layer, in J m-2.
  type(thermo_var_ptrs),   intent(inout) :: tv     !< A structure containing pointers to any
                                                   !! available thermodynamic fields. Absent fields
                                                   !! have NULL ptrs.
  type(forcing),           intent(inout) :: fluxes !< A structure containing pointers to any
                                                   !! possible forcing fields. Unused fields have
                                                   !! NULL ptrs.
  real,                    intent(in)    :: dt     !< Time increment, in s.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
                           intent(out)   :: kd_int !< The diagnosed diffusivities at interfaces,
                                                   !! in m2 s-1.
  type(energetic_pbl_cs),  pointer       :: cs     !< The control structure returned by a previous
                                                   !! call to mixedlayer_init.
  real, dimension(SZI_(G),SZJ_(G)), &
                           intent(in)    :: buoy_flux !< The surface buoyancy flux. in m2/s3.
  real,          optional, intent(in)    :: dt_diag   !< The diagnostic time step, which may be less
                                                   !! than dt if there are two callse to
                                                   !! mixedlayer, in s.
  logical,       optional, intent(in)    :: last_call !< If true, this is the last call to
                                                   !! mixedlayer in the current time step, so
                                                   !! diagnostics will be written. The default
                                                   !! is .true.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                 optional, intent(out)   :: dt_expected, ds_expected

!    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.
!
!    For each interior interface, first discard the TKE to account for mixing
! of shortwave radiation through the next denser cell.  Next drive mixing based
! on the local? values of ustar + wstar, subject to available energy.  This
! step sets the value of Kd(K).  Any remaining energy is then subject to decay
! before being handed off to the next interface.  mech_TKE and conv_PErel are treated
! separately for the purposes of decay, but are used proportionately to drive
! mixing.
!
!   The key parameters for the mixed layer are found in the control structure.
! These include mstar, nstar, TKE_decay, and conv_decay.  For the Oberhuber (1993) mixed layer,
! the values of these are:
!      mstar = 1.25,  nstar = 1, TKE_decay = 2.5, conv_decay = 0.5
!  TKE_decay is 1/kappa in eq. 28 of Oberhuber (1993), while
!  conv_decay is 1/mu.
!    For a traditional Kraus-Turner mixed layer, the values are:
!      mstar = 1.25, nstar = 0.4, TKE_decay = 0.0, conv_decay = 0.0

! Arguments: h_3d - Layer thickness, in m or kg m-2. (Intent in/out)
!                   The units of h are referred to as H below.
!  (in)      u_3d - Zonal velocities interpolated to h points, m s-1.
!  (in)      v_3d - Zonal velocities interpolated to h points, m s-1.
!  (in/out)  tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (in)      fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      dt - Time increment, in s.
!  (out)     Kd_int - The diagnosed diffusivities at interfaces, in m2 s-1.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 mixedlayer_init.
!  (in)      dSV_dT - The partial derivative of in-situ specific volume with
!                     potential temperature, in m3 kg-1 K-1.
!  (in)      dSV_dS - The partial derivative of in-situ specific volume with
!                     salinity, in m3 kg-1 ppt-1.
!  (in)      TKE_forced - The forcing requirements to homogenize the forcing
!                 that has been applied to each layer through each layer, in J m-2.
!  (in)      Buoy_Flux - The surface buoyancy flux. in m2/s3.
!  (in,opt)  dt_diag - The diagnostic time step, which may be less than dt
!                      if there are two callse to mixedlayer, in s.
!  (in,opt)  last_call - if true, this is the last call to mixedlayer in the
!                        current time step, so diagnostics will be written.
!                        The default is .true.

  real, dimension(SZI_(G),SZK_(GV)) :: &
    h, &            !   The layer thickness, in H (usually m or kg m-2).
    t, &            !   The layer temperatures, in deg C.
    s, &            !   The layer salinities, in psu.
    u, &            !   The zonal velocity, in m s-1.
    v               !   The meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZK_(GV)+1) :: &
    kd, &           ! The diapycnal diffusivity, in m2 s-1.
    pres, &         ! Interface pressures in Pa.
    hb_hs           ! The distance from the bottom over the thickness of the
                    ! water column, nondim.
  real, dimension(SZI_(G)) :: &
    mech_tke, &     !   The mechanically generated turbulent kinetic energy
                    ! available for mixing over a time step, in J m-2 = kg s-2.
    conv_perel, &   ! The potential energy that has been convectively released
                    ! during this timestep, in J m-2 = kg s-2. A portion nstar_FC
                    ! of conv_PErel is available to drive mixing.
    htot, &         !   The total depth of the layers above an interface, in H.
    uhtot, &        !   The depth integrated zonal and meridional velocities in the
    vhtot, &        ! layers above, in H m s-1.
    mech_tke_top, & !    The value of mech_TKE at the top of the column, in J m-2.
    conv_perel_top, & !  The value of conv_PErel at the top of the column, in J m-2.

    idecay_len_tke, &  ! The inverse of a turbulence decay length scale, in H-1.
    h_sum, &        ! The total thickness of the water column, in H.
    absf            ! The absolute value of f, in s-1.


  real, dimension(SZI_(G),SZK_(GV)) :: &
    dt_to_dcolht, & ! Partial derivatives of the total column height with the temperature
    ds_to_dcolht, & ! and salinity changes within a layer, in m K-1 and m ppt-1.
    dt_to_dpe, &    ! Partial derivatives of column potential energy with the temperature
    ds_to_dpe, &    ! and salinity changes within a layer, in J m-2 K-1 and J m-2 ppt-1.
    dt_to_dcolht_a, & ! Partial derivatives of the total column height with the temperature
    ds_to_dcolht_a, & ! and salinity changes within a layer, including the implicit effects
                    ! of mixing with layers higher in the water colun, in m K-1 and m ppt-1.
    dt_to_dpe_a, &  ! Partial derivatives of column potential energy with the temperature
    ds_to_dpe_a     ! and salinity changes within a layer, including the implicit effects
                    ! of mixing with layers higher in the water column, in
                    ! units of J m-2 K-1 and J m-2 ppt-1.
  real, dimension(SZK_(GV)) :: &
    t0, s0, &       ! Initial values of T and S in the column, in K and ppt.
    te, se, &       ! Estimated final values of T and S in the column, in K and ppt.
    c1, &           ! c1 is used by the tridiagonal solver, ND.
    dte, dse        ! Running (1-way) estimates of temperature and salinity change.
  real, dimension(SZK_(GV)) :: &
    th_a, &         ! An effective temperature times a thickness in the layer above,
                    ! including implicit mixing effects with other yet higher layers, in K H.
    sh_a, &         ! An effective salinity times a thickness in the layer above,
                    ! including implicit mixing effects with other yet higher layers, in K H.
    th_b, &         ! An effective temperature times a thickness in the layer below,
                    ! including implicit mixing effects with other yet lower layers, in K H.
    sh_b            ! An effective salinity times a thickness in the layer below,
                    ! including implicit mixing effects with other yet lower layers, in K H.
  real, dimension(SZI_(G)) :: &
    hp_a            ! An effective pivot thickness of the layer including the effects
                    ! of coupling with layers above, in H.  This is the first term
                    ! in the denominator of b1 in a downward-oriented tridiagonal solver.
  real, dimension(SZK_(GV)+1) :: &
    mixlen_shape, & ! A nondimensional shape factor for the mixing length that
                    ! gives it an appropriate assymptotic value at the bottom of
                    ! the boundary layer.
    kddt_h          ! The diapycnal diffusivity times a timestep divided by the
                    ! average thicknesses around a layer, in H (m or kg m-2).
  real :: b1        ! b1 is inverse of the pivot used by the tridiagonal solver, in H-1.
  real :: h_neglect ! A thickness that is so small it is usually lost
                    ! in roundoff and can be neglected, in H.
  real :: dmass     ! The mass per unit area within a layer, in kg m-2.
  real :: dpres     ! The hydrostatic pressure change across a layer, in Pa.
  real :: dmke_max  ! The maximum amount of mean kinetic energy that could be
                    ! converted to turbulent kinetic energy if the velocity in
                    ! the layer below an interface were homogenized with all of
                    ! the water above the interface, in J m-2 = kg s-2.
  real :: mke2_hharm ! Twice the inverse of the harmonic mean of the thickness
                    ! of a layer and the thickness of the water above, used in
                    ! the MKE conversion equation, in H-1.

  real :: dt_h      ! The timestep divided by the averages of the thicknesses around
                    ! a layer, times a thickness conversion factor, in H s m-2.
  real :: h_bot     ! The distance from the bottom, in H.
  real :: h_rsum    ! The running sum of h from the top, in H.
  real :: i_hs      ! The inverse of h_sum, in H-1.
  real :: i_mld     ! The inverse of the current value of MLD, in H-1.
  real :: h_tt      ! The distance from the surface or up to the next interface
                    ! that did not exhibit turbulent mixing from this scheme plus
                    ! a surface mixing roughness length given by h_tt_min, in H.
  real :: h_tt_min  ! A surface roughness length, in H.

  real :: c1_3      ! = 1/3.
  real :: vonkar    ! The vonKarman constant.
  real :: i_dtrho   ! 1.0 / (dt * Rho0) in m3 kg-1 s-1.  This is
                    ! used convert TKE back into ustar^3.
  real :: u_star    ! The surface friction velocity, in m s-1.
  real :: u_star_mean ! The surface friction without gustiness in m s-1.
  real :: vstar     ! An in-situ turbulent velocity, in m s-1.
  real :: enhance_m ! An enhancement factor for vstar, based here on Langmuir impact.
  real :: la        ! The Langmuir number (non-dim)
  real :: lamod     ! A modified Langmuir number accounting for other parameters.
  real :: hbs_here  ! The local minimum of hb_hs and MixLen_shape, times a
                    ! conversion factor from H to M, in m H-1.
  real :: nstar_fc  ! The fraction of conv_PErel that can be converted to mixing, nondim.
  real :: tke_reduc ! The fraction by which TKE and other energy fields are
                    ! reduced to support mixing, nondim. between 0 and 1.
  real :: tot_tke   ! The total TKE available to support mixing at interface K, in J m-2.
  real :: tke_here  ! The total TKE at this point in the algorithm, in J m-2.
  real :: dt_km1_t2 ! A diffusivity-independent term related to the temperature
                    ! change in the layer above the interface, in K.
  real :: ds_km1_t2 ! A diffusivity-independent term related to the salinity
                    ! change in the layer above the interface, in ppt.
  real :: dte_term  ! A diffusivity-independent term related to the temperature
                    ! change in the layer below the interface, in K H.
  real :: dse_term  ! A diffusivity-independent term related to the salinity
                    ! change in the layer above the interface, in ppt H.
  real :: dte_t2    ! A part of dTe_term, in K H.
  real :: dse_t2    ! A part of dSe_term, in ppt H.
  real :: dpe_conv  ! The convective change in column potential energy, in J m-2.
  real :: mke_src   ! The mean kinetic energy source of TKE due to Kddt_h(K), in J m-2.
  real :: dmke_src_dk  ! The partial derivative of MKE_src with Kddt_h(K), in J m-2 H-1.
  real :: kd_guess0, pe_chg_g0, dpea_dkd_g0, kddt_h_g0
  real :: pe_chg_max   ! The maximum PE change for very large values of Kddt_h(K).
  real :: dpec_dkd_kd0 ! The partial derivative of PE change with Kddt_h(K)
                       ! for very small values of Kddt_h(K), in J m-2 H-1.
  real :: pe_chg    ! The change in potential energy due to mixing at an
                    ! interface, in J m-2, positive for the column increasing
                    ! in potential energy (i.e., consuming TKE).
  real :: tke_left  ! The amount of turbulent kinetic energy left for the most
                    ! recent guess at Kddt_h(K), in J m-2.
  real :: dpec_dkd  ! The partial derivative of PE_chg with Kddt_h(K), in J m-2 H-1.
  real :: tke_left_min, tke_left_max, kddt_h_max, kddt_h_min
  real :: kddt_h_guess ! A guess at the value of Kddt_h(K), in H.
  real :: kddt_h_next  ! The next guess at the value of Kddt_h(K), in H.
  real :: dkddt_h      ! The change between guesses at Kddt_h(K), in H.
  real :: dkddt_h_newt ! The change between guesses at Kddt_h(K) with Newton's method, in H.
  real :: kddt_h_newt  ! The Newton's method next guess for Kddt_h(K), in H.
  real :: exp_kh    ! The nondimensional decay of TKE across a layer, ND.
  logical :: use_newt  ! Use Newton's method for the next guess at Kddt_h(K).
  logical :: convectively_stable
  logical, dimension(SZI_(G)) :: &
    sfc_connected   ! If true the ocean is actively turbulent from the present
                    ! interface all the way up to the surface.
  logical :: sfc_disconnect ! If true, any turbulence has become disconnected
                    ! from the surface.

! The following is only used as a diagnostic.
  real :: dt__diag  ! A copy of dt_diag (if present) or dt, in s.
  real :: idtdr0    !  = 1.0 / (dt__diag * Rho0), in m3 kg-1 s-1.
  real, dimension(SZI_(G),SZJ_(G)) :: &
    hsfc_used   ! The thickness of the surface region after buffer layer
  logical :: write_diags  ! If true, write out diagnostics with this step.
  logical :: reset_diags  ! If true, zero out the accumulated diagnostics.
                ! detrainment, in units of m.
  ! Local column copies of energy change diagnostics, all in J m-2.
  real :: dtke_conv, dtke_forcing, dtke_mixing
  real :: dtke_mke, dtke_mech_decay, dtke_conv_decay
  !----------------------------------------------------------------------
  !/BGR added Aug24,2016 for adding iteration to get boundary layer depth
  !    - needed to compute new mixing length.
  real :: mld_guess, mld_found ! Mixing Layer depth guessed/found for iteration, in m.
  real :: max_mld, min_mld ! Iteration bounds, in m, which are adjusted at each step
                           !  - These are initialized based on surface/bottom
                           !  1. The iteration guesses a value (possibly from
                           !     prev step or neighbor).
                           !  2. The iteration checks if value is converged,
                           !     too shallow, or too deep.
                           !  3. Based on result adjusts the Max/Min
                           !     and searches through the water column.
                           !  - If using an accurate guess the iteration
                           !    is very quick (e.g. if MLD doesn't change
                           !    over timestep).  Otherwise it takes 5-10
                           !    passes, but has a high convergence rate.
                           !    Other iteration may be tried, but this
                           !    method seems to rarely fail and the added
                           !    cost is likely not significant.  Additionally,
                           !    when it fails it does so in a reasonable
                           !    manner giving a usable guess. When it
                           !    does fail, it is due to convection within
                           !    the boundary.  Likely, a new method e.g.
                           !    surface_disconnect, can improve this.
  logical :: first_obl     ! Flag for computing "found" Mixing layer depth
  logical :: obl_converged ! Flag for convergence of MLD
  integer :: obl_it        ! Iteration counter
!### These need to be made into run-time parameters.
  integer :: max_obl_it=20 ! Set maximum number of iterations.  Probably
                           !  best as an input parameter, but then may want
                           !  to use allocatable arrays if storing
                           !  guess/found (as diagnostic); skipping for now.
                           !  In reality, the maximum number of guesses
                           !  needed is set by:
                           !    DEPTH/2^M < DZ
                           !    where M is the number of guesses
                           !    e.g. M=12 for DEPTH=4000m and DZ=1m
  real, dimension(SZK_(GV)+1) :: vstar_used, &      ! 1D arrays used to store
                               mixing_length_used   ! Vstar and Mixing_Length
  !/BGR - remaining variables are related to tracking iteration statistics.
  logical :: obl_it_stats=.false. ! Flag for computing OBL iteration statistics
  REAL :: itguess(20), itresult(20),itmax(20),itmin(20) ! Flag for storing guess/result
                                                        ! should have dim=MAX_OBL_IT
  integer, save :: maxit=0   ! Stores maximum number of iterations
  integer, save :: minit=1e8 ! Stores minimum number of iterations
  integer, save :: sumit=0   ! Stores total iterations (summed over all)
  integer, save :: numit=0   ! Stores number of times iterated
                             !e.g. Average iterations = SUMIT/NUMIT
  integer, save :: converged!
  integer, save :: notconverged!
  !-End BGR iteration parameters-----------------------------------------
  real :: n2_dissipation
  real :: bf_stable ! Buoyancy flux, capped at 0 (negative only)
  real :: bf_unstable ! Buoyancy flux, floored at 0 (positive only)
  real :: stab_scale ! Composite of Stabilizing length scales:
                     !  Ekman scale and Monin-Obukhov scale.
  real :: il_ekman ! Inverse of Ekman length scale
  real :: il_obukhov ! Inverse of Obukhov length scale
  real :: mld_o_ekman          ! >
  real :: mld_o_obukhov_stab   ! Ratios of length scales where MLD is boundary layer depth
  real :: ekman_o_obukhov_stab ! >
  real :: mld_o_obukhov_un   ! Ratios of length scales where MLD is boundary layer depth
  real :: ekman_o_obukhov_un ! >

  real :: c_mo = 1. ! Constant in STAB_SCALE for Monin-Obukhov
  real :: c_ek = 2. ! Constant in STAB_SCALE for Ekman length
  real :: mld_over_stab ! Mixing layer depth divided by STAB_SCALE
  real :: mstar_mix! The value of mstar (Proportionality of TKE to drive mixing to ustar
                    ! cubed) which is computed as a function of latitude, boundary layer depth,
                    ! and the Monin-Obukhov depth.
  real :: mstar_lt ! The added mstar contribution due to Langmuir turbulence
  real :: mstar_conv_adj ! Adjustment made to mstar due to convection reducing mechanical mixing.
  real :: mstar_stab, mstar_rot ! Mstar in each limit, max is used.
  logical :: debug=.false.  ! Change this hard-coded value for debugging.
!  The following arrays are used only for debugging purposes.
  real :: dpe_debug, mixing_debug, taux2, tauy2
  real, dimension(20) :: tke_left_itt, pe_chg_itt, kddt_h_itt, dpea_dkd_itt, mke_src_itt
  real, dimension(SZI_(G),SZK_(GV)) :: &
    mech_tke_k, conv_perel_k
  real, dimension(SZK_(GV)) :: nstar_k
  integer, dimension(SZK_(GV)) :: num_itts

  integer :: i, j, k, is, ie, js, je, nz, itt, max_itt

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

  if (.not. associated(cs)) call mom_error(fatal, "energetic_PBL: "//&
         "Module must be initialized before it is used.")

  if (.not. ASSOCIATED(tv%eqn_of_state)) call mom_error(fatal, &
      "energetic_PBL: Temperature, salinity and an equation of state "//&
      "must now be used.")
  if (.NOT. ASSOCIATED(fluxes%ustar)) call mom_error(fatal, &
      "energetic_PBL: No surface TKE fluxes (ustar) defined in mixedlayer!")
  if (present(dt_expected) .or. present(ds_expected)) debug = .true.

  h_neglect = gv%H_subroundoff

  if(.not.cs%Use_MLD_Iteration) max_obl_it=1
  c1_3 = 1.0 / 3.0
  dt__diag = dt ; if (present(dt_diag)) dt__diag = dt_diag
  idtdr0 = 1.0 / (dt__diag * gv%Rho0)
  write_diags = .true. ; if (present(last_call)) write_diags = last_call
  max_itt = 20

  h_tt_min = 0.0
  vonkar = 0.41
  mstar_mix=cs%MSTAR!Initialize to mstar
  i_dtrho = 0.0 ; if (dt*gv%Rho0 > 0.0) i_dtrho = 1.0 / (dt*gv%Rho0)

  ! Determine whether to zero out diagnostics before accumulation.
  reset_diags = .true.
  if (present(dt_diag) .and. write_diags .and. (dt__diag > dt)) &
    reset_diags = .false.  ! This is the second call to mixedlayer.

  if (reset_diags) then
!!OMP parallel default(none) shared(is,ie,js,je,CS)
    if (cs%TKE_diagnostics) then
!!OMP do
      do j=js,je ; do i=is,ie
        cs%diag_TKE_wind(i,j) = 0.0 ; cs%diag_TKE_MKE(i,j) = 0.0
        cs%diag_TKE_conv(i,j) = 0.0 ; cs%diag_TKE_forcing(i,j) = 0.0
        cs%diag_TKE_mixing(i,j) = 0.0 ; cs%diag_TKE_mech_decay(i,j) = 0.0
        cs%diag_TKE_conv_decay(i,j) = 0.0 !; CS%diag_TKE_unbalanced_forcing(i,j) = 0.0
      enddo ; enddo
    endif
    if (cs%Mixing_Diagnostics) then
      cs%Mixing_Length(:,:,:) = 0.0
      cs%Velocity_Scale(:,:,:) = 0.0
    endif
!!OMP end parallel
  endif


!!OMP parallel do default(none) shared(js,je,nz,is,ie,h_3d,u_3d,v_3d,tv,dt,      &
!!OMP                                  CS,G,GV,fluxes,IdtdR0,                    &
!!OMP                                  TKE_forced,debug,H_neglect,dSV_dT,        &
!!OMP                                  dSV_dS,I_dtrho,C1_3,h_tt_min,vonKar,     &
!!OMP                                  max_itt,Kd_int)                           &
!!OMP                           private(i,j,k,h,u,v,T,S,Kd,mech_TKE_k,conv_PErel_k,    &
!!OMP                                   U_Star,absf,mech_TKE,conv_PErel,nstar_k, &
!!OMP                                   h_sum,I_hs,h_bot,hb_hs,T0,S0,num_itts,   &
!!OMP                                   pres,dMass,dPres,dT_to_dPE,dS_to_dPE,    &
!!OMP                                   dT_to_dColHt,dS_to_dColHt,Kddt_h,hp_a,   &
!!OMP                                   Th_a,Sh_a,Th_b,Sh_b,dT_to_dPE_a,htot,    &
!!OMP                                   dT_to_dColHt_a,dS_to_dColHt_a,uhtot,vhtot, &
!!OMP                                   Idecay_len_TKE,exp_kh,nstar_FC,tot_TKE,  &
!!OMP                                   TKE_reduc,dTe_t2,dSe_t2,dTe,dSe,dt_h,    &
!!OMP                                   Convectively_stable,sfc_disconnect,b1,   &
!!OMP                                   c1,dT_km1_t2,dS_km1_t2,dTe_term,         &
!!OMP                                   dSe_term,MKE2_Hharm,vstar,h_tt,h_rsum,   &
!!OMP                                   Kd_guess0,Kddt_h_g0,dPEc_dKd_Kd0,        &
!!OMP                                   PE_chg_max,dPEa_dKd_g0,PE_chg_g0,        &
!!OMP                                   MKE_src,dPE_conv,Kddt_h_max,Kddt_h_min,  &
!!OMP                                   dTKE_conv, dTKE_forcing, dTKE_mixing,    &
!!OMP                                   dTKE_MKE,dTKE_mech_decay,dTKE_conv_decay,&
!!OMP                                   TKE_left_max,TKE_left_min,Kddt_h_guess,  &
!!OMP                                   TKE_left_itt,dPEa_dKd_itt,PE_chg_itt,    &
!!OMP                                   MKE_src_itt,Kddt_h_itt,dPEc_dKd,PE_chg,  &
!!OMP                                   dMKE_src_dK,TKE_left,use_Newt,           &
!!OMP                                   dKddt_h_Newt,Kddt_h_Newt,Kddt_h_next,    &
!!OMP                                   dKddt_h,Te,Se,Hsfc_used,dS_to_dPE_a,     &
!!OMP                                   dMKE_max,sfc_connected,TKE_here)
  do j=js,je
    ! Copy the thicknesses and other fields to 2-d arrays.
    do k=1,nz ; do i=is,ie
      h(i,k) = h_3d(i,j,k) + h_neglect ; u(i,k) = u_3d(i,j,k) ; v(i,k) = v_3d(i,j,k)
      t(i,k) = tv%T(i,j,k) ; s(i,k) = tv%S(i,j,k)
      kd(i,k) = 0.0
    enddo ; enddo
    do i=is,ie
       cs%ML_depth(i,j) = h(i,1)*gv%H_to_m
       !CS%ML_depth2(i,j) = h(i,1)*GV%H_to_m
       sfc_connected(i) = .true.
    enddo

    if (debug) then
      mech_tke_k(:,:) = 0.0 ; conv_perel_k(:,:) = 0.0
    endif

    !   Determine the initial mech_TKE and conv_PErel, including the energy required
    ! to mix surface heating through the topmost cell, the energy released by mixing
    ! surface cooling & brine rejection down through the topmost cell, and
    ! homogenizing the shortwave heating within that cell.  This sets the energy
    ! and ustar and wstar available to drive mixing at the first interior
    ! interface.
    do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then


      u_star = fluxes%ustar(i,j)
      taux2 = 0.
      if ((g%mask2dCu(i-1,j) + g%mask2dCu(i,j)) > 0) &
        taux2 = (g%mask2dCu(i-1,j)*fluxes%taux(i-1,j)**2 + &
                 g%mask2dCu(i,j)*fluxes%taux(i,j)**2) / (g%mask2dCu(i-1,j) + g%mask2dCu(i,j))
      tauy2 = 0.0
      if ((g%mask2dCv(i,j-1) + g%mask2dCv(i,j)) > 0) &
        tauy2 = (g%mask2dCv(i,j-1)*fluxes%tauy(i,j-1)**2 + &
                 g%mask2dCv(i,j)*fluxes%tauy(i,j)**2) / (g%mask2dCv(i,j-1) + g%mask2dCv(i,j))
      u_star_mean = sqrt(sqrt(taux2 + tauy2)/gv%rho0)
      if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then
        if (fluxes%frac_shelf_h(i,j) > 0.0) &
          u_star = (1.0 - fluxes%frac_shelf_h(i,j)) * u_star + &
                    fluxes%frac_shelf_h(i,j) * fluxes%ustar_shelf(i,j)
      endif
      if (u_star < cs%ustar_min) u_star = cs%ustar_min
      ! Computing Bf w/ limiters.
      bf_stable = max(0.0,buoy_flux(i,j)) ! Positive for stable
      bf_unstable = min(0.0,buoy_flux(i,j)) ! Negative for unstable
      if (cs%omega_frac >= 1.0) then ; absf(i) = 2.0*cs%omega
      else
        absf(i) = 0.25*((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
                     (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i-1,j))))
        if (cs%omega_frac > 0.0) &
          absf(i) = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf(i)**2)
      endif
      ! Computing stability scale which correlates with TKE for mixing, where
      ! TKE for mixing = TKE production minus TKE dissipation
      stab_scale = u_star**2 / ( vonkar * ( c_mo * bf_stable/u_star -  c_ek * u_star * absf(i)))
      ! Inverse of Ekman and Obukhov
      il_ekman   = absf(i)/u_star
      il_obukhov = buoy_flux(i,j)*vonkar/u_star**3

      if (cs%Mstar_Mode.eq.cs%CONST_MSTAR) then
        mech_tke(i) = (dt*cs%mstar*gv%Rho0)*((u_star**3))
        conv_perel(i) = 0.0

        if (cs%TKE_diagnostics) then
          cs%diag_TKE_wind(i,j) = cs%diag_TKE_wind(i,j) + mech_tke(i) * idtdr0
          if (tke_forced(i,j,1) <= 0.0) then
            cs%diag_TKE_forcing(i,j) = cs%diag_TKE_forcing(i,j) + &
                 max(-mech_tke(i), tke_forced(i,j,1)) * idtdr0
            ! CS%diag_TKE_unbalanced_forcing(i,j) = CS%diag_TKE_unbalanced_forcing(i,j) + &
            !     min(0.0, TKE_forced(i,j,1) + mech_TKE(i)) * IdtdR0
          else
            cs%diag_TKE_forcing(i,j) = cs%diag_TKE_forcing(i,j) + cs%nstar*tke_forced(i,j,1) * idtdr0
          endif
        endif

        if (tke_forced(i,j,1) <= 0.0) then
          mech_tke(i) = mech_tke(i) + tke_forced(i,j,1)
          if (mech_tke(i) < 0.0) mech_tke(i) = 0.0
        else
          conv_perel(i) = conv_perel(i) + tke_forced(i,j,1)
        endif

      endif

!    endif ; enddo

!    do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then

      h_sum(i) = h_neglect ; do k=1,nz ; h_sum(i) = h_sum(i) + h(i,k) ; enddo
      i_hs = 0.0 ; if (h_sum(i) > 0.0) i_hs = 1.0 / h_sum(i)

      h_bot = 0.0 ; hb_hs(i,nz+1) = 0.0
      do k=nz,1,-1
        h_bot = h_bot + h(i,k)
        hb_hs(i,k) = h_bot * i_hs
      enddo

      pres(i,1) = 0.0
      do k=1,nz
        dmass = gv%H_to_kg_m2 * h(i,k)
        dpres = gv%g_Earth * dmass
        dt_to_dpe(i,k) = (dmass * (pres(i,k) + 0.5*dpres)) * dsv_dt(i,j,k)
        ds_to_dpe(i,k) = (dmass * (pres(i,k) + 0.5*dpres)) * dsv_ds(i,j,k)
        dt_to_dcolht(i,k) = dmass * dsv_dt(i,j,k)
        ds_to_dcolht(i,k) = dmass * dsv_ds(i,j,k)

        pres(i,k+1) = pres(i,k) + dpres
      enddo

!    endif ; enddo

    ! Note the outer i-loop and inner k-loop loop order!!!
!    do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
      do k=1,nz ; t0(k) = t(i,k) ; s0(k) = s(i,k) ; enddo

      ! Store the initial mechanical TKE and convectively released PE to
      ! enable multiple iterations.
      mech_tke_top(i) = mech_tke(i) ; conv_perel_top(i) = conv_perel(i)

      !/The following lines are for the iteration over MLD
      !{
      ! max_MLD will initialized as ocean bottom depth
      max_mld = 0.0 ; do k=1,nz ; max_mld = max_mld + h(i,k)*gv%H_to_m ; enddo
      min_mld = 0.0 !min_MLD will initialize as 0.
      !/BGR: May add user-input bounds for max/min MLD

      !/BGR: Add MLD_guess based on stored previous value.
      !      note that this is different from ML_Depth already
      !      computed by EPBL, need to figure out why.
      if (cs%MLD_iteration_guess .and. cs%ML_Depth2(i,j) > 1.) then
        !If prev value is present use for guess.
        mld_guess=cs%ML_Depth2(i,j)
      else
        !Otherwise guess middle of water column (or stab_scale if smaller).
        mld_guess = 0.5 * (min_mld+max_mld)
      endif

      ! Iterate up to MAX_OBL_IT times to determine a converged EPBL depth.
      obl_converged = .false.
      ! Initialize ENHANCE_M to 1 and mstar_lt to 0
      enhance_m=1.e0
      mstar_lt = 0.0
      do obl_it=1,max_obl_it ; if (.not. obl_converged) then

        ! Reset ML_depth
        cs%ML_depth(i,j) = h(i,1)*gv%H_to_m
        !CS%ML_depth2(i,j) = h(i,1)*GV%H_to_m

        sfc_connected(i) = .true.

        if (cs%Mstar_Mode.gt.0) then
        ! Note the value of mech_TKE(i) now must be iterated over, so it is moved here
        ! First solve for the TKE to PE length scale
          if (cs%MSTAR_MODE.eq.cs%MLD_o_OBUKHOV) then
            mld_over_stab = mld_guess / stab_scale - cs%MSTAR_XINT
            if ((mld_over_stab) .le. 0.0) then
              !Asymptote to 0 as MLD_over_Stab -> -infinity (always)
              mstar_mix = (cs%MSTAR_B*(mld_over_stab)+cs%MSTAR_A)**(cs%MSTAR_N)
            else
              if (cs%MSTAR_CAP>=0.) then
                if (cs%MSTAR_FLATCAP .OR. (mld_over_stab .le.cs%MSTAR_XINT_UP)) then
                !If using flat cap (or if using asymptotic cap
                !   but within linear regime we can make use of same code)
                  mstar_mix = min(cs%MSTAR_CAP, &
                     cs%MSTAR_SLOPE*(mld_over_stab)+cs%MSTAR_AT_XINT)
                else
                  !Asymptote to MSTAR_CAP as MLD_over_Stab -> infinity
                  mstar_mix = cs%MSTAR_CAP - &
                       (cs%MSTAR_B2*(mld_over_stab-cs%MSTAR_XINT_UP)&
                       +cs%MSTAR_A2)**(cs%MSTAR_N)
                endif
              else
                !No cap if negative cap value given.
                mstar_mix = cs%MSTAR_SLOPE*(mld_over_stab)+cs%MSTAR_AT_XINT
              endif
            endif
          elseif (cs%MSTAR_MODE.eq.cs%EKMAN_o_OBUKHOV) then
            mstar_stab = cs%MSTAR_COEF*sqrt(bf_stable/u_star**2/(absf(i)+1.e-10))
            mstar_rot =  cs%C_EK*log(max(1.,u_star/(absf(i)+1.e-10)/mld_guess))
            if ( cs%MSTAR_CAP.le.0.0) then !No cap.
              mstar_mix = max(mstar_stab,& ! 1st term if balance of rotation and stabilizing
                                    ! the balance is f(L_Ekman,L_Obukhov)
                              min(& ! 2nd term for forced stratification limited
                                 1.25,& !.5/von Karman (Obukhov limit)
                                 ! 3rd term for rotation (Ekman length) limited
                                 mstar_rot))
            else
              mstar_mix = min( & ! Sets a cap. The cap should be large and just
                               !  meant to be a safety net.
                           cs%MSTAR_CAP, &
                           max(mstar_stab,& ! 1st term if balance of rotation and stabilizing
                                 ! the balance is f(L_Ekman,L_Obukhov)
                              min(& ! 2nd term for forced stratification limited
                                 1.25,& !.5/von Karman (Obukhov limit)
                                 ! 3rd term for rotation (Ekman length) limited
                                 mstar_rot)))
            endif!cap for mstar_mode==2
          endif!mstar_mode==1 or ==2
          ! Adjustment for unstable buoyancy flux.
          !  Convection reduces mechanical mixing because there
          !  is less density gradient to mix. (Statically unstable near surface)
          mstar_conv_adj = 1. - cs%CNV_MST_FAC * (-bf_unstable+1.e-10) / &
                                 ( (-bf_unstable+1.e-10)+                 &
                                   2. *mstar_mix *u_star**3 / mld_guess )
          if (cs%Use_LA_windsea) then
            ! 1. Get LA
            call get_la_windsea( u_star_mean, mld_guess*cs%LaDepthRatio, gv, la)
            ! 2. Get parameters for modified LA
            mld_o_ekman = abs(mld_guess*il_ekman)
            mld_o_obukhov_stab = abs(max(0.,mld_guess*il_obukhov))
            mld_o_obukhov_un = abs(min(0.,mld_guess*il_obukhov))
            ekman_o_obukhov_stab = abs(max(0.,il_obukhov/(il_ekman+1.e-10)))
            ekman_o_obukhov_un = abs(min(0.,il_obukhov/(il_ekman+1.e-10)))
            ! 3. Adjust LA based on various parameters.
            !    Assumes linear factors based on length scale ratios to adjust LA
            !    Note when these coefficients are set to 0 recovers simple LA.
            lamod = la * (1.0 + max(-0.5,cs%LaC_MLDoEK * mld_o_ekman) +   &
                 cs%LaC_EKoOB_stab * ekman_o_obukhov_stab + &
                 cs%LaC_EKoOB_un * ekman_o_obukhov_un + &
                 cs%LaC_MLDoOB_stab * mld_o_obukhov_stab  + &
                 cs%LaC_MLDoOB_un * mld_o_obukhov_un )
            if (cs%LT_Enhance_Form==1) then
              !Original w'/ust scaling w/ Van Roekel's scaling
              enhance_m = (1+(1.4*la)**(-2)+(5.4*la)**(-4))**(1.5)
            elseif (cs%LT_Enhance_Form==2) then
              ! Enhancement is multiplied (added mst_lt set to 0)
              enhance_m = min(cs%Max_Enhance_M,(1.+cs%LT_ENHANCE_COEF*lamod**cs%LT_ENHANCE_EXP))
              mstar_lt = 0.0
            elseif (cs%LT_ENHANCE_Form == 3) then
              ! or Enhancement is additive (multiplied enhance_m set to 1)
              mstar_lt = cs%LT_ENHANCE_COEF * lamod**cs%LT_ENHANCE_EXP
              enhance_m = 1.0
            endif
          endif
          !Reset mech_tke and conv_perel values (based on new mstar)
          mech_tke(i) = ( mstar_mix * mstar_conv_adj * enhance_m + mstar_lt) * (dt*gv%Rho0*u_star**3)
          conv_perel(i) = 0.0
          if (cs%TKE_diagnostics) then
            cs%diag_TKE_wind(i,j) = cs%diag_TKE_wind(i,j) + mech_tke(i) * idtdr0
            if (tke_forced(i,j,1) <= 0.0) then
              cs%diag_TKE_forcing(i,j) = cs%diag_TKE_forcing(i,j) + &
                   max(-mech_tke(i), tke_forced(i,j,1)) * idtdr0
              ! CS%diag_TKE_unbalanced_forcing(i,j) = CS%diag_TKE_unbalanced_forcing(i,j) + &
              !     min(0.0, TKE_forced(i,j,1) + mech_TKE(i)) * IdtdR0
            else
              cs%diag_TKE_forcing(i,j) = cs%diag_TKE_forcing(i,j) + cs%nstar*tke_forced(i,j,1) * idtdr0
            endif
          endif

          if (tke_forced(i,j,1) <= 0.0) then
            mech_tke(i) = mech_tke(i) + tke_forced(i,j,1)
            if (mech_tke(i) < 0.0) mech_tke(i) = 0.0
          else
            conv_perel(i) = conv_perel(i) + tke_forced(i,j,1)
          endif
        else
          mech_tke(i) = mech_tke_top(i)*enhance_m ; conv_perel(i) = conv_perel_top(i)
        endif

        if (cs%TKE_diagnostics) then
          dtke_conv = 0.0 ; dtke_forcing = 0.0 ; dtke_mixing = 0.0
          dtke_mke = 0.0 ; dtke_mech_decay = 0.0 ; dtke_conv_decay = 0.0
        endif

        ! Store in 1D arrays cleared out each iteration.  Only write in
        !  3D arrays after convergence.
        do k=1,nz
          vstar_used(k) = 0.0 ; mixing_length_used(k) = 0.0
        enddo
        if (.not.cs%Use_MLD_Iteration) obl_converged = .true.

        if ((.not.cs%Use_MLD_Iteration) .or. &
            (cs%transLay_scale >= 1.0) .or. (cs%transLay_scale < 0.0) ) then
          do k=1,nz+1 ; mixlen_shape(k) = 1.0 ; enddo
        elseif (mld_guess <= 0.0) then
          if (cs%transLay_scale > 0.0) then
            do k=1,nz+1 ; mixlen_shape(k) = cs%transLay_scale ; enddo
          else
            do k=1,nz+1 ; mixlen_shape(k) = 1.0 ; enddo
          endif
        else
          ! Reduce the mixing length based on MLD, with a quadratic
          ! expression that follows KPP.
          i_mld = 1.0 / mld_guess ; h_rsum = 0.0
          mixlen_shape(1) = 1.0
          do k=2,nz+1
            h_rsum = h_rsum + h(i,k-1)
            if (cs%MixLenExponent==2.0)then
              mixlen_shape(k) = cs%transLay_scale + (1.0 - cs%transLay_scale) * &
                   (max(0.0, (mld_guess - h_rsum)*i_mld) )**2!CS%MixLenExponent
            else
              mixlen_shape(k) = cs%transLay_scale + (1.0 - cs%transLay_scale) * &
                   (max(0.0, (mld_guess - h_rsum)*i_mld) )**cs%MixLenExponent
            endif
          enddo
        endif


        kd(i,1) = 0.0 ; kddt_h(1) = 0.0
        hp_a(i) = h(i,1)
        dt_to_dpe_a(i,1) = dt_to_dpe(i,1) ; dt_to_dcolht_a(i,1) = dt_to_dcolht(i,1)
        ds_to_dpe_a(i,1) = ds_to_dpe(i,1) ; ds_to_dcolht_a(i,1) = ds_to_dcolht(i,1)

        htot(i) = h(i,1) ; uhtot(i) = u(i,1)*h(i,1) ; vhtot(i) = v(i,1)*h(i,1)

        if (debug) then
          mech_tke_k(i,1) = mech_tke(i) ; conv_perel_k(i,1) = conv_perel(i)
          nstar_k(:) = 0.0 ; nstar_k(1) = cs%nstar ; num_itts(:) = -1
        endif

        do k=2,nz
          ! Apply dissipation to the TKE, here applied as an exponential decay
          ! due to 3-d turbulent energy being lost to inefficient rotational modes.

          !   There should be several different "flavors" of TKE that decay at
          ! different rates.  The following form is often used for mechanical
          ! stirring from the surface, perhaps due to breaking surface gravity
          ! waves and wind-driven turbulence.
          idecay_len_tke(i) = (cs%TKE_decay * absf(i) / u_star) * gv%H_to_m
          exp_kh = 1.0
          if (idecay_len_tke(i) > 0.0) exp_kh = exp(-h(i,k-1)*idecay_len_tke(i))
          if (cs%TKE_diagnostics) &
            dtke_mech_decay = dtke_mech_decay + (exp_kh-1.0) * mech_tke(i) * idtdr0
          mech_tke(i) = mech_tke(i) * exp_kh

          !   Accumulate any convectively released potential energy to contribute
          ! to wstar and to drive penetrating convection.
          if (tke_forced(i,j,k) > 0.0) then
            conv_perel(i) = conv_perel(i) + tke_forced(i,j,k)
            if (cs%TKE_diagnostics) &
              dtke_forcing = dtke_forcing + cs%nstar*tke_forced(i,j,k) * idtdr0
          endif

          if (debug) then
            mech_tke_k(i,k) = mech_tke(i) ; conv_perel_k(i,k) = conv_perel(i)
          endif

          !  Determine the total energy
          nstar_fc = cs%nstar
          if (cs%nstar * conv_perel(i) > 0.0) then
            ! Here nstar is a function of the natural Rossby number 0.2/(1+0.2/Ro), based
            ! on a curve fit from the data of Wang (GRL, 2003).
            ! Note:         Ro = 1.0 / sqrt(0.5 * dt * Rho0 * (absf*htot(i))**3 / conv_PErel(i))
            nstar_fc = cs%nstar * conv_perel(i) / (conv_perel(i) + 0.2 * &
                            sqrt(0.5 * dt * gv%Rho0 * (absf(i)*(htot(i)*gv%H_to_m))**3 * conv_perel(i)))
          endif
          if (debug) nstar_k(k) = nstar_fc

          tot_tke = mech_tke(i) + nstar_fc * conv_perel(i)

          !   For each interior interface, first discard the TKE to account for
          ! mixing of shortwave radiation through the next denser cell.
          if (tke_forced(i,j,k) < 0.0) then
            if (tke_forced(i,j,k) + tot_tke < 0.0) then
              ! The shortwave requirements deplete all the energy in this layer.
              if (cs%TKE_diagnostics) then
                dtke_mixing = dtke_mixing + tot_tke * idtdr0
                dtke_forcing = dtke_forcing - tot_tke * idtdr0
                ! dTKE_unbalanced_forcing = dTKE_unbalanced_forcing + &
                !     (TKE_forced(i,j,k) + tot_TKE) * IdtdR0
                dtke_conv_decay = dtke_conv_decay + &
                        (cs%nstar-nstar_fc) * conv_perel(i) * idtdr0
              endif
              tot_tke = 0.0 ; mech_tke(i) = 0.0 ; conv_perel(i) = 0.0
            else
              ! Reduce the mechanical and convective TKE proportionately.
              tke_reduc = (tot_tke + tke_forced(i,j,k)) / tot_tke
              if (cs%TKE_diagnostics) then
                dtke_mixing = dtke_mixing - tke_forced(i,j,k) * idtdr0
                dtke_forcing = dtke_forcing + tke_forced(i,j,k) * idtdr0
                dtke_conv_decay = dtke_conv_decay + &
                    (1.0-tke_reduc)*(cs%nstar-nstar_fc) * conv_perel(i) * idtdr0
              endif
              tot_tke = tke_reduc*tot_tke   ! = tot_TKE + TKE_forced(i,j,k)
              mech_tke(i) = tke_reduc*mech_tke(i)
              conv_perel(i) = tke_reduc*conv_perel(i)
            endif
          endif

          ! Precalculate some temporary expressions that are independent of Kddt_h(K).
          if (cs%orig_PE_calc) then
            if (k==2) then
              dte_t2 = 0.0 ; dse_t2 = 0.0
            else
              dte_t2 = kddt_h(k-1) * ((t0(k-2) - t0(k-1)) + dte(k-2))
              dse_t2 = kddt_h(k-1) * ((s0(k-2) - s0(k-1)) + dse(k-2))
            endif
          endif
          dt_h = (gv%m_to_H**2*dt) / max(0.5*(h(i,k-1)+h(i,k)), 1e-15*h_sum(i))

          !   This tests whether the layers above and below this interface are in
          ! a convetively stable configuration, without considering any effects of
          ! mixing at higher interfaces.  It is an approximation to the more
          ! complete test dPEc_dKd_Kd0 >= 0.0, that would include the effects of
          ! mixing across interface K-1.  The dT_to_dColHt here are effectively
          ! mass-weigted estimates of dSV_dT.
          convectively_stable = ( 0.0 <= &
            ( (dt_to_dcolht(i,k) + dt_to_dcolht(i,k-1) ) * (t0(k-1)-t0(k)) + &
              (ds_to_dcolht(i,k) + ds_to_dcolht(i,k-1) ) * (s0(k-1)-s0(k)) ) )

          if ((mech_tke(i) + conv_perel(i)) <= 0.0 .and. convectively_stable) then
            ! Energy is already exhausted, so set Kd = 0 and cycle or exit?
            tot_tke = 0.0 ; mech_tke(i) = 0.0 ; conv_perel(i) = 0.0
            kd(i,k) = 0.0 ; kddt_h(k) = 0.0
            sfc_disconnect = .true.
            ! if (.not.debug) exit

           !   The estimated properties for layer k-1 can be calculated, using
           ! greatly simplified expressions when Kddt_h = 0.  This enables the
           ! tridiagonal solver for the whole column to be completed for debugging
           ! purposes, and also allows for something akin to convective adjustment
           ! in unstable interior regions?
            b1 = 1.0 / hp_a(i)
            c1(k) = 0.0
            if (cs%orig_PE_calc) then
              dte(k-1) = b1 * ( dte_t2 )
              dse(k-1) = b1 * ( dse_t2 )
            endif

            hp_a(i) = h(i,k)
            dt_to_dpe_a(i,k) = dt_to_dpe(i,k)
            ds_to_dpe_a(i,k) = ds_to_dpe(i,k)
            dt_to_dcolht_a(i,k) = dt_to_dcolht(i,k)
            ds_to_dcolht_a(i,k) = ds_to_dcolht(i,k)

          else ! tot_TKE > 0.0 or this is a potentially convectively unstable profile.
            sfc_disconnect = .false.

            ! Precalculate some more temporary expressions that are independent of
            ! Kddt_h(K).
            if (cs%orig_PE_calc) then
              if (k==2) then
                dt_km1_t2 = (t0(k)-t0(k-1))
                ds_km1_t2 = (s0(k)-s0(k-1))
              else
                dt_km1_t2 = (t0(k)-t0(k-1)) - &
                      (kddt_h(k-1) / hp_a(i)) * ((t0(k-2) - t0(k-1)) + dte(k-2))
                ds_km1_t2 = (s0(k)-s0(k-1)) - &
                      (kddt_h(k-1) / hp_a(i)) * ((s0(k-2) - s0(k-1)) + dse(k-2))
              endif
              dte_term = dte_t2 + hp_a(i) * (t0(k-1)-t0(k))
              dse_term = dse_t2 + hp_a(i) * (s0(k-1)-s0(k))
            else
              if (k<=2) then
                th_a(k-1) = h(i,k-1) * t0(k-1) ; sh_a(k-1) = h(i,k-1) * s0(k-1)
              else
                th_a(k-1) = h(i,k-1) * t0(k-1) + kddt_h(k-1) * te(k-2)
                sh_a(k-1) = h(i,k-1) * s0(k-1) + kddt_h(k-1) * se(k-2)
              endif
              th_b(k) = h(i,k) * t0(k) ; sh_b(k) = h(i,k) * s0(k)
            endif

            !   Using Pr=1 and the diffusivity at the bottom interface (once it is
            ! known), determine how much resolved mean kinetic energy (MKE) will be
            ! extracted within a timestep and add a fraction CS%MKE_to_TKE_effic of
            ! this to the mTKE budget available for mixing in the next layer.

            if ((cs%MKE_to_TKE_effic > 0.0) .and. (htot(i)*h(i,k) > 0.0)) then
              ! This is the energy that would be available from homogenizing the
              ! velocities between layer k and the layers above.
              dmke_max = (gv%H_to_kg_m2 * cs%MKE_to_TKE_effic) * 0.5 * &
                  (h(i,k) / ((htot(i) + h(i,k))*htot(i))) * &
                  ((uhtot(i)-u(i,k)*htot(i))**2 + (vhtot(i)-v(i,k)*htot(i))**2)
              ! A fraction (1-exp(Kddt_h*MKE2_Hharm)) of this energy would be
              ! extracted by mixing with a finite viscosity.
              mke2_hharm = (htot(i) + h(i,k) + 2.0*h_neglect) / &
                           ((htot(i)+h_neglect) * (h(i,k)+h_neglect))
            else
              dmke_max = 0.0 ; mke2_hharm = 0.0
            endif

            ! At this point, Kddt_h(K) will be unknown because its value may depend
            ! on how much energy is available.  mech_TKE might be negative due to
            ! contributions from TKE_forced.
            h_tt = htot(i) + h_tt_min
            tke_here = mech_tke(i) + cs%wstar_ustar_coef*conv_perel(i)
            if (tke_here > 0.0) then
              vstar = cs%vstar_scale_fac * (i_dtrho*tke_here)**c1_3
              hbs_here = gv%H_to_m * min(hb_hs(i,k), mixlen_shape(k))
              mixing_length_used(k) = max(cs%min_mix_len,((h_tt*hbs_here)*vstar) / &
                  ((cs%Ekman_scale_coef * absf(i)) * (h_tt*hbs_here) + vstar))
              !Note setting Kd_guess0 to Mixing_Length_Used(K) here will
              ! change the answers.  Therefore, skipping that.
              if (.not.cs%Use_MLD_Iteration) then
                 kd_guess0 = vstar * vonkar *  ((h_tt*hbs_here)*vstar) / &
                  ((cs%Ekman_scale_coef * absf(i)) * (h_tt*hbs_here) + vstar)
              else
                 kd_guess0 = vstar * vonkar * mixing_length_used(k)
              endif
            else
              vstar = 0.0 ; kd_guess0 = 0.0
            endif
            vstar_used(k) = vstar ! Track vstar
            kddt_h_g0 = kd_guess0*dt_h

            if (cs%orig_PE_calc) then
              call find_pe_chg_orig(kddt_h_g0, h(i,k), hp_a(i), dte_term, dse_term, &
                       dt_km1_t2, ds_km1_t2, dt_to_dpe(i,k), ds_to_dpe(i,k), &
                       dt_to_dpe_a(i,k-1), ds_to_dpe_a(i,k-1), &
                       pres(i,k), dt_to_dcolht(i,k), ds_to_dcolht(i,k), &
                       dt_to_dcolht_a(i,k-1), ds_to_dcolht_a(i,k-1), &
                       pe_chg=pe_chg_g0, dpec_dkd=dpea_dkd_g0, dpe_max=pe_chg_max, &
                       dpec_dkd_0=dpec_dkd_kd0 )
            else
              call find_pe_chg(0.0, kddt_h_g0, hp_a(i), h(i,k), &
                         th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
                         dt_to_dpe_a(i,k-1), ds_to_dpe_a(i,k-1), dt_to_dpe(i,k), ds_to_dpe(i,k), &
                         pres(i,k), dt_to_dcolht_a(i,k-1), ds_to_dcolht_a(i,k-1), &
                         dt_to_dcolht(i,k), ds_to_dcolht(i,k), &
                         pe_chg=pe_chg_g0, dpec_dkd=dpea_dkd_g0, dpe_max=pe_chg_max, &
                         dpec_dkd_0=dpec_dkd_kd0 )
            endif

            mke_src = dmke_max*(1.0 - exp(-kddt_h_g0 * mke2_hharm))

            if (pe_chg_g0 .gt. 0.0) then
              !Negative buoyancy (increases PE)
              n2_dissipation = 1.+cs%N2_DISSIPATION_SCALE_NEG
            else
              !Positive buoyancy (decreases PE)
              n2_dissipation = 1.+cs%N2_DISSIPATION_SCALE_POS
            endif

            if ((pe_chg_g0 < 0.0) .or. ((vstar == 0.0) .and. (dpec_dkd_kd0 < 0.0))) then
              ! This column is convectively unstable.
              if (pe_chg_max <= 0.0) then
                ! Does MKE_src need to be included in the calculation of vstar here?
                tke_here = mech_tke(i) + cs%wstar_ustar_coef*(conv_perel(i)-pe_chg_max)
                if (tke_here > 0.0) then
                  vstar = cs%vstar_scale_fac * (i_dtrho*tke_here)**c1_3
                  hbs_here = gv%H_to_m * min(hb_hs(i,k), mixlen_shape(k))
                  mixing_length_used(k) = max(cs%min_mix_len,((h_tt*hbs_here)*vstar) / &
                      ((cs%Ekman_scale_coef * absf(i)) * (h_tt*hbs_here) + vstar))
                  if (.not.cs%Use_MLD_Iteration) then
                  ! Note again (as prev) that using Mixing_Length_Used here
                  !  instead of redoing the computation will change answers...
                     kd(i,k) = vstar * vonkar *  ((h_tt*hbs_here)*vstar) / &
                          ((cs%Ekman_scale_coef * absf(i)) * (h_tt*hbs_here) + vstar)
                  else
                     kd(i,k) = vstar * vonkar * mixing_length_used(k)
                  endif
                else
                  vstar = 0.0 ; kd(i,k) = 0.0
                endif
                vstar_used(k) = vstar

                if (cs%orig_PE_calc) then
                  call find_pe_chg_orig(kd(i,k)*dt_h, h(i,k), hp_a(i), dte_term, dse_term, &
                           dt_km1_t2, ds_km1_t2, dt_to_dpe(i,k), ds_to_dpe(i,k), &
                           dt_to_dpe_a(i,k-1), ds_to_dpe_a(i,k-1), &
                           pres(i,k), dt_to_dcolht(i,k), ds_to_dcolht(i,k), &
                           dt_to_dcolht_a(i,k-1), ds_to_dcolht_a(i,k-1), &
                           pe_chg=dpe_conv)
                else
                  call find_pe_chg(0.0, kd(i,k)*dt_h, hp_a(i), h(i,k), &
                           th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
                           dt_to_dpe_a(i,k-1), ds_to_dpe_a(i,k-1), dt_to_dpe(i,k), ds_to_dpe(i,k), &
                           pres(i,k), dt_to_dcolht_a(i,k-1), ds_to_dcolht_a(i,k-1), &
                           dt_to_dcolht(i,k), ds_to_dcolht(i,k), &
                           pe_chg=dpe_conv)
                endif
                ! Should this be iterated to convergence for Kd?
                if (dpe_conv > 0.0) then
                  kd(i,k) = kd_guess0 ; dpe_conv = pe_chg_g0
                else
                  mke_src = dmke_max*(1.0 - exp(-(kd(i,k)*dt_h) * mke2_hharm))
                endif
              else
                ! The energy change does not vary monotonically with Kddt_h.  Find the maximum?
                kd(i,k) = kd_guess0 ; dpe_conv = pe_chg_g0
              endif
              conv_perel(i) = conv_perel(i) - dpe_conv
              mech_tke(i) = mech_tke(i) + mke_src
              if (cs%TKE_diagnostics) then
                dtke_conv = dtke_conv - cs%nstar*dpe_conv * idtdr0
                dtke_mke = dtke_mke + mke_src * idtdr0
              endif
              if (sfc_connected(i)) then
                cs%ML_depth(i,j) = cs%ML_depth(i,j) + gv%H_to_m * h(i,k)
                !CS%ML_depth2(i,j) = CS%ML_depth2(i,J) + GV%H_to_m * h(i,k)
              endif

              kddt_h(k) = kd(i,k)*dt_h
            elseif (tot_tke + (mke_src - n2_dissipation*pe_chg_g0) >= 0.0) then
              ! There is energy to support the suggested mixing.  Keep that estimate.
              kd(i,k) = kd_guess0
              kddt_h(k) = kddt_h_g0

              ! Reduce the mechanical and convective TKE proportionately.
              tot_tke = tot_tke + mke_src
              tke_reduc = 0.0   ! tot_TKE could be 0 if Convectively_stable is false.
              if (tot_tke > 0.0) tke_reduc = (tot_tke - n2_dissipation*pe_chg_g0) &
                                             / tot_tke
              if (cs%TKE_diagnostics) then
                dtke_mixing = dtke_mixing - pe_chg_g0 * idtdr0
                dtke_mke = dtke_mke + mke_src * idtdr0
                dtke_conv_decay = dtke_conv_decay + &
                    (1.0-tke_reduc)*(cs%nstar-nstar_fc) * conv_perel(i) * idtdr0
              endif
              tot_tke = tke_reduc*tot_tke
              mech_tke(i) = tke_reduc*(mech_tke(i) + mke_src)
              conv_perel(i) = tke_reduc*conv_perel(i)
              if (sfc_connected(i)) then
                cs%ML_depth(i,j) = cs%ML_depth(i,j) + gv%H_to_m * h(i,k)
                !CS%ML_depth2(i,J) = CS%ML_depth2(i,J) + GV%H_to_m * h(i,k)
              endif
            elseif (tot_tke == 0.0) then
              ! This can arise if nstar_FC = 0.
              kd(i,k) = 0.0 ; kddt_h(k) = 0.0
              tot_tke = 0.0 ; conv_perel(i) = 0.0 ; mech_tke(i) = 0.0
              sfc_disconnect = .true.
            else
              ! There is not enough energy to support the mixing, so reduce the
              ! diffusivity to what can be supported.
              kddt_h_max = kddt_h_g0 ; kddt_h_min = 0.0
              tke_left_max = tot_tke + (mke_src - n2_dissipation*pe_chg_g0) ;
              tke_left_min = tot_tke

              ! As a starting guess, take the minimum of a false position estimate
              ! and a Newton's method estimate starting from Kddt_h = 0.0.
              kddt_h_guess = tot_tke * kddt_h_max / max( n2_dissipation*pe_chg_g0 &
                                 - mke_src, kddt_h_max * (dpec_dkd_kd0 - dmke_max *        &
                                 mke2_hharm) )
              ! The above expression is mathematically the same as the following
              ! except it is not susceptible to division by zero when
              !   dPEc_dKd_Kd0 = dMKE_max = 0 .
              !  Kddt_h_guess = tot_TKE * min( Kddt_h_max / (PE_chg_g0 - MKE_src), &
              !                      1.0 / (dPEc_dKd_Kd0 - dMKE_max * MKE2_Hharm) )
              if (debug) then
                tke_left_itt(:) = 0.0 ; dpea_dkd_itt(:) = 0.0 ; pe_chg_itt(:) = 0.0
                mke_src_itt(:) = 0.0 ; kddt_h_itt(:) = 0.0
              endif
              do itt=1,max_itt
                if (cs%orig_PE_calc) then
                  call find_pe_chg_orig(kddt_h_guess, h(i,k), hp_a(i), dte_term, dse_term, &
                           dt_km1_t2, ds_km1_t2, dt_to_dpe(i,k), ds_to_dpe(i,k), &
                           dt_to_dpe_a(i,k-1), ds_to_dpe_a(i,k-1), &
                           pres(i,k), dt_to_dcolht(i,k), ds_to_dcolht(i,k), &
                           dt_to_dcolht_a(i,k-1), ds_to_dcolht_a(i,k-1), &
                           pe_chg=pe_chg, dpec_dkd=dpec_dkd )
                else
                  call find_pe_chg(0.0, kddt_h_guess, hp_a(i), h(i,k), &
                           th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
                           dt_to_dpe_a(i,k-1), ds_to_dpe_a(i,k-1), dt_to_dpe(i,k), ds_to_dpe(i,k), &
                           pres(i,k), dt_to_dcolht_a(i,k-1), ds_to_dcolht_a(i,k-1), &
                           dt_to_dcolht(i,k), ds_to_dcolht(i,k), &
                           pe_chg=dpe_conv)
                endif
                mke_src = dmke_max * (1.0 - exp(-mke2_hharm * kddt_h_guess))
                dmke_src_dk = dmke_max * mke2_hharm * exp(-mke2_hharm * kddt_h_guess)

                tke_left = tot_tke + (mke_src - n2_dissipation*pe_chg)
                if (debug) then
                  kddt_h_itt(itt) = kddt_h_guess ; mke_src_itt(itt) = mke_src
                  pe_chg_itt(itt) = n2_dissipation*pe_chg
                  tke_left_itt(itt) = tke_left
                  dpea_dkd_itt(itt) = dpec_dkd
                endif
                ! Store the new bounding values, bearing in mind that min and max
                ! here refer to Kddt_h and dTKE_left/dKddt_h < 0:
                if (tke_left >= 0.0) then
                  kddt_h_min = kddt_h_guess ; tke_left_min = tke_left
                else
                  kddt_h_max = kddt_h_guess ; tke_left_max = tke_left
                endif

                ! Try to use Newton's method, but if it would go outside the bracketed
                ! values use the false-position method instead.
                use_newt = .true.
                if (dpec_dkd*n2_dissipation - dmke_src_dk <= 0.0) then
                  use_newt = .false.
                else
                  dkddt_h_newt = tke_left / (dpec_dkd*n2_dissipation - dmke_src_dk)
                  kddt_h_newt = kddt_h_guess + dkddt_h_newt
                  if ((kddt_h_newt > kddt_h_max) .or. (kddt_h_newt < kddt_h_min)) &
                    use_newt = .false.
                endif

                if (use_newt) then
                  kddt_h_next = kddt_h_guess + dkddt_h_newt
                  dkddt_h = dkddt_h_newt
                else
                  kddt_h_next = (tke_left_max * kddt_h_min - kddt_h_max * tke_left_min) / &
                                (tke_left_max - tke_left_min)
                  dkddt_h = kddt_h_next - kddt_h_guess
                endif

                if ((abs(dkddt_h) < 1e-9*kddt_h_guess) .or. (itt==max_itt)) then
                  ! Use the old value so that the energy calculation does not need to be repeated.
                  if (debug) num_itts(k) = itt
                  exit
                else
                  kddt_h_guess = kddt_h_next
                endif
              enddo
              kd(i,k) = kddt_h_guess / dt_h ; kddt_h(k) = kd(i,k)*dt_h

              ! All TKE should have been consumed.
              if (cs%TKE_diagnostics) then
                dtke_mixing = dtke_mixing - (tot_tke + mke_src) * idtdr0
                dtke_mke = dtke_mke + mke_src * idtdr0
                dtke_conv_decay = dtke_conv_decay + &
                    (cs%nstar-nstar_fc) * conv_perel(i) * idtdr0
              endif

              if (sfc_connected(i)) cs%ML_depth(i,j) = cs%ML_depth(i,j) + &
                   (pe_chg / pe_chg_g0) * gv%H_to_m * h(i,k)
              tot_tke = 0.0 ; mech_tke(i) = 0.0 ; conv_perel(i) = 0.0
              sfc_disconnect = .true.
            endif

            kddt_h(k) = kd(i,k)*dt_h
           !   At this point, the final value of Kddt_h(K) is known, so the
           ! estimated properties for layer k-1 can be calculated.
            b1 = 1.0 / (hp_a(i) + kddt_h(k))
            c1(k) = kddt_h(k) * b1
            if (cs%orig_PE_calc) then
              dte(k-1) = b1 * ( kddt_h(k)*(t0(k)-t0(k-1)) + dte_t2 )
              dse(k-1) = b1 * ( kddt_h(k)*(s0(k)-s0(k-1)) + dse_t2 )
            endif

            hp_a(i) = h(i,k) + (hp_a(i) * b1) * kddt_h(k)
            dt_to_dpe_a(i,k) = dt_to_dpe(i,k) + c1(k)*dt_to_dpe_a(i,k-1)
            ds_to_dpe_a(i,k) = ds_to_dpe(i,k) + c1(k)*ds_to_dpe_a(i,k-1)
            dt_to_dcolht_a(i,k) = dt_to_dcolht(i,k) + c1(k)*dt_to_dcolht_a(i,k-1)
            ds_to_dcolht_a(i,k) = ds_to_dcolht(i,k) + c1(k)*ds_to_dcolht_a(i,k-1)

          endif  ! tot_TKT > 0.0 branch.  Kddt_h(K) has been set.

          ! Store integrated velocities and thicknesses for MKE conversion calculations.
          if (sfc_disconnect) then
            ! There is no turbulence at this interface, so zero out the running sums.
            uhtot(i) = u(i,k)*h(i,k)
            vhtot(i) = v(i,k)*h(i,k)
            htot(i)  = h(i,k)
            sfc_connected(i) = .false.
          else
            uhtot(i) = uhtot(i) + u(i,k)*h(i,k)
            vhtot(i) = vhtot(i) + v(i,k)*h(i,k)
            htot(i)  = htot(i) + h(i,k)
          endif

          if (debug) then
            if (k==2) then
              te(1) = b1*(h(i,1)*t0(1))
              se(1) = b1*(h(i,1)*s0(1))
            else
              te(k-1) = b1 * (h(i,k-1) * t0(k-1) + kddt_h(k-1) * te(k-2))
              se(k-1) = b1 * (h(i,k-1) * s0(k-1) + kddt_h(k-1) * se(k-2))
            endif
          endif
        enddo
        kd(i,nz+1) = 0.0

        if (debug) then
          ! Complete the tridiagonal solve for Te.
          b1 = 1.0 / hp_a(i)
          te(nz) = b1 * (h(i,nz) * t0(nz) + kddt_h(nz) * te(nz-1))
          se(nz) = b1 * (h(i,nz) * s0(nz) + kddt_h(nz) * se(nz-1))
          do k=nz-1,1,-1
            te(k) = te(k) + c1(k+1)*te(k+1)
            se(k) = se(k) + c1(k+1)*se(k+1)
          enddo
        endif
        if (present(dt_expected)) then
          do k=1,nz ; dt_expected(i,j,k) = te(k) - t0(k) ; enddo
        endif
        if (present(ds_expected)) then
          do k=1,nz ; ds_expected(i,j,k) = se(k) - s0(k) ; enddo
        endif
        if (debug) then
          dpe_debug = 0.0
          do k=1,nz
            dpe_debug = dpe_debug + (dt_to_dpe(i,k) * (te(k) - t0(k)) + &
                                     ds_to_dpe(i,k) * (se(k) - s0(k)))
          enddo
          mixing_debug = dpe_debug * idtdr0
        endif
        k = nz ! This is here to allow a breakpoint to be set.
        !/BGR
        ! The following lines are used for the iteration
        ! note the iteration has been altered to use the value predicted by
        ! the TKE threshold (ML_DEPTH).  This is because the MSTAR
        ! is now dependent on the ML, and therefore the ML needs to be estimated
        ! more precisely than the grid spacing.
        !/
        itmax(obl_it) = max_mld       ! Track max    }
        itmin(obl_it) = min_mld       ! Track min    } For debug purpose
        itguess(obl_it) = mld_guess   ! Track guess  }
        !/
        mld_found=0.0 ; first_obl=.true.
        if (cs%Orig_MLD_iteration) then
          !This is how the iteration was original conducted
          do k=2,nz
            if (first_obl) then !Breaks when OBL found
              if (vstar_used(k) > 1.e-10 .and. k < nz) then
                mld_found = mld_found+h(i,k-1)*gv%H_to_m
              else
                first_obl = .false.
                if (mld_found-cs%MLD_tol > mld_guess) then
                  min_mld = mld_guess
                elseif ((mld_guess-mld_found) < max(cs%MLD_tol,h(i,k-1)*gv%H_to_m)) then
                  obl_converged = .true.!Break convergence loop
                  if (obl_it_stats) then !Compute iteration statistics
                    maxit = max(maxit,obl_it)
                    minit = min(minit,obl_it)
                    sumit = sumit+obl_it
                    numit = numit+1
                    print*,maxit,minit,sumit/numit
                  endif
                  cs%ML_Depth2(i,j) = mld_guess
                else
                  max_mld = mld_guess !We know this guess was too deep
                endif
              endif
            endif
          enddo
        else
          !New method uses ML_DEPTH as computed in ePBL routine
          mld_found=cs%ML_DEPTH(i,j)
          if (mld_found-cs%MLD_tol > mld_guess) then
            min_mld = mld_guess
          elseif (abs(mld_guess-mld_found) < (cs%MLD_tol)) then
            obl_converged = .true.!Break convergence loop
            if (obl_it_stats) then !Compute iteration statistics
              maxit = max(maxit,obl_it)
              minit = min(minit,obl_it)
              sumit = sumit+obl_it
              numit = numit+1
              print*,maxit,minit,sumit/numit
            endif
            cs%ML_Depth2(i,j) = mld_guess
          else
            max_mld = mld_guess !We know this guess was too deep
          endif
        endif
        ! For next pass, guess average of minimum and maximum values.
        mld_guess = min_mld*0.5 + max_mld*0.5
        itresult(obl_it) = mld_found
      endif ; enddo ! Iteration loop for converged boundary layer thickness.
      if (.not.obl_converged) then
        !/Temp output, warn that EPBL didn't converge
        !/Print guess/found for every iteration step
        !print*,'EPBL MLD DID NOT CONVERGE'
        notconverged=notconverged+1
        !do obl_it=1,max_obl_it
        !   print*,ITmin(obl_it),ITmax(obl_it)
        !   print*,obl_it,ITguess(obl_it),ITresult(obl_it)
        !enddo
        !Activate to print out some output when not converged
        !{
        !print*,'Min/Max: ',ITmin(50),ITmax(50)
        !print*,'Guess/result: ',ITguess(50),ITresult(50)
        !print*,'Stats on CPU: ',CONVERGED,NOTCONVERGED,&
        !     real(NOTCONVERGED)/real(CONVERGED)
        !}
        !stop !Kill if not converged during testing.
      else
        converged=converged+1
      endif

      if (cs%TKE_diagnostics) then
        cs%diag_TKE_MKE(i,j) = cs%diag_TKE_MKE(i,j) + dtke_mke
        cs%diag_TKE_conv(i,j) = cs%diag_TKE_conv(i,j) + dtke_conv
        cs%diag_TKE_forcing(i,j) = cs%diag_TKE_forcing(i,j) + dtke_forcing
        cs%diag_TKE_mixing(i,j) = cs%diag_TKE_mixing(i,j) + dtke_mixing
        cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + dtke_mech_decay
        cs%diag_TKE_conv_decay(i,j) = cs%diag_TKE_conv_decay(i,j) + dtke_conv_decay
       ! CS%diag_TKE_unbalanced_forcing(i,j) = CS%diag_TKE_unbalanced_forcing(i,j) + dTKE_unbalanced
      endif
      if (cs%Mixing_Diagnostics) then
        !Write to 3-D for outputing Mixing length and
        !  velocity scale.
        do k=1,nz
          cs%Mixing_Length(i,j,k) = mixing_length_used(k)
          cs%Velocity_Scale(i,j,k) = vstar_used(k)
        enddo
      endif
      if (allocated(cs%Enhance_M)) cs%Enhance_M(i,j) = enhance_m
      if (allocated(cs%mstar_mix)) cs%mstar_mix(i,j) = mstar_mix
      if (allocated(cs%MLD_Obukhov)) cs%MLD_Obukhov(i,j) = (mld_guess*il_obukhov)
      if (allocated(cs%MLD_Ekman)) cs%MLD_Ekman(i,j) = (mld_guess*il_ekman)
      if (allocated(cs%Ekman_Obukhov)) cs%Ekman_Obukhov(i,j) = (il_obukhov/(il_ekman+1.e-10))
      if (allocated(cs%La)) cs%La(i,j) = la
      if (allocated(cs%La_mod)) cs%La_mod(i,j) = lamod
    else
      ! For masked points, Kd_int must still be set (to 0) because it has intent(out).
      do k=1,nz+1
        kd(i,k) = 0.
      enddo
      if (present(dt_expected)) then
        do k=1,nz ; dt_expected(i,j,k) = 0.0 ; enddo
      endif
      if (present(ds_expected)) then
        do k=1,nz ; ds_expected(i,j,k) = 0.0 ; enddo
      endif
    endif ; enddo ; ! Close of i-loop - Note unusual loop order!

    if (cs%id_Hsfc_used > 0) then
      do i=is,ie ; hsfc_used(i,j) = h(i,1)*gv%H_to_m ; enddo
      do k=2,nz ; do i=is,ie
        if (kd(i,k) > 0.0) hsfc_used(i,j) = hsfc_used(i,j) + h(i,k)*gv%H_to_m
      enddo ; enddo
    endif

    do k=1,nz+1 ; do i=is,ie
      kd_int(i,j,k) = kd(i,k)
    enddo ; enddo

  enddo ! j-loop

  if (write_diags) then
    if (cs%id_ML_depth > 0) &
      call post_data(cs%id_ML_depth, cs%ML_depth, cs%diag)
    if (cs%id_TKE_wind > 0) &
      call post_data(cs%id_TKE_wind, cs%diag_TKE_wind, cs%diag)
    if (cs%id_TKE_MKE > 0) &
      call post_data(cs%id_TKE_MKE, cs%diag_TKE_MKE, cs%diag)
    if (cs%id_TKE_conv > 0) &
      call post_data(cs%id_TKE_conv, cs%diag_TKE_conv, cs%diag)
    if (cs%id_TKE_forcing > 0) &
      call post_data(cs%id_TKE_forcing, cs%diag_TKE_forcing, cs%diag)
    if (cs%id_TKE_mixing > 0) &
      call post_data(cs%id_TKE_mixing, cs%diag_TKE_mixing, cs%diag)
    if (cs%id_TKE_mech_decay > 0) &
      call post_data(cs%id_TKE_mech_decay, cs%diag_TKE_mech_decay, cs%diag)
    if (cs%id_TKE_conv_decay > 0) &
      call post_data(cs%id_TKE_conv_decay, cs%diag_TKE_conv_decay, cs%diag)
    if (cs%id_Hsfc_used > 0) &
      call post_data(cs%id_Hsfc_used, hsfc_used, cs%diag)
    if (cs%id_Mixing_Length > 0) &
      call post_data(cs%id_Mixing_Length, cs%Mixing_Length, cs%diag)
    if (cs%id_Velocity_Scale >0) &
      call post_data(cs%id_Velocity_Scale, cs%Velocity_Scale, cs%diag)
    if (cs%id_OSBL >0) &
      call post_data(cs%id_OSBL, cs%ML_Depth2, cs%diag)
    if (cs%id_LT_Enhancement >0) &
      call post_data(cs%id_LT_Enhancement, cs%Enhance_M, cs%diag)
    if (cs%id_MSTAR_MIX >0) &
      call post_data(cs%id_MSTAR_MIX, cs%MSTAR_MIX, cs%diag)
    if (cs%id_MLD_OBUKHOV >0) &
      call post_data(cs%id_MLD_Obukhov, cs%MLD_OBUKHOV, cs%diag)
    if (cs%id_MLD_EKMAN >0) &
      call post_data(cs%id_MLD_Ekman, cs%MLD_EKMAN, cs%diag)
    if (cs%id_Ekman_Obukhov >0) &
      call post_data(cs%id_Ekman_Obukhov, cs%Ekman_Obukhov, cs%diag)
    if (cs%id_LA >0) &
      call post_data(cs%id_LA, cs%LA, cs%diag)
    if (cs%id_LA_MOD >0) &
      call post_data(cs%id_LA_MOD, cs%LA_MOD, cs%diag)

  endif

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

)

Definition at line 2371 of file MOM_energetic_PBL.F90.

  type(energetic_pbl_cs), pointer :: cs

  if (.not.associated(cs)) return

  if (allocated(cs%ML_depth))            deallocate(cs%ML_depth)
  if (allocated(cs%ML_depth2))           deallocate(cs%ML_depth2)
  if (allocated(cs%Enhance_M))           deallocate(cs%Enhance_M)
  if (allocated(cs%MLD_EKMAN))           deallocate(cs%MLD_EKMAN)
  if (allocated(cs%MLD_OBUKHOV))         deallocate(cs%MLD_OBUKHOV)
  if (allocated(cs%EKMAN_OBUKHOV))       deallocate(cs%EKMAN_OBUKHOV)
  if (allocated(cs%LA))                  deallocate(cs%LA)
  if (allocated(cs%LA_mod))              deallocate(cs%LA_mod)
  if (allocated(cs%MSTAR_MIX))           deallocate(cs%MSTAR_MIX)
  if (allocated(cs%diag_TKE_wind))       deallocate(cs%diag_TKE_wind)
  if (allocated(cs%diag_TKE_MKE))        deallocate(cs%diag_TKE_MKE)
  if (allocated(cs%diag_TKE_conv))       deallocate(cs%diag_TKE_conv)
  if (allocated(cs%diag_TKE_forcing))    deallocate(cs%diag_TKE_forcing)
  if (allocated(cs%diag_TKE_mixing))     deallocate(cs%diag_TKE_mixing)
  if (allocated(cs%diag_TKE_mech_decay)) deallocate(cs%diag_TKE_mech_decay)
  if (allocated(cs%diag_TKE_conv_decay)) deallocate(cs%diag_TKE_conv_decay)
  if (allocated(cs%Mixing_Length))       deallocate(cs%Mixing_Length)
  if (allocated(cs%Velocity_Scale))      deallocate(cs%Velocity_Scale)

  deallocate(cs)

energetic_pbl_get_mld()

subroutine, public mom_energetic_pbl::energetic_pbl_get_mld	(	type(energetic_pbl_cs), pointer	CS,
		real, dimension(szi_(g),szj_(g)), intent(out)	MLD,
		type(ocean_grid_type), intent(in)	G
	)

Copies the ePBL active mixed layer depth into MLD.

Parameters

	cs	Control structure for ePBL
[in]	g	Grid structure
[out]	mld	Depth of ePBL active mixing layer

Definition at line 1908 of file MOM_energetic_PBL.F90.

Referenced by mom_diabatic_driver::diabatic().

  type(energetic_pbl_cs),           pointer     :: cs  !< Control structure for ePBL
  type(ocean_grid_type),            intent(in)  :: g   !< Grid structure
  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: mld !< Depth of ePBL active mixing layer
  ! Local variables
  integer :: i,j
  do j = g%jsc, g%jec ; do i = g%isc, g%iec
    mld(i,j) = cs%ML_depth(i,j)
  enddo ; enddo

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

Parameters

[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[in]	param_file	A structure to parse for run-time parameters

Definition at line 2056 of file MOM_energetic_PBL.F90.

  type(time_type), target, intent(in)    :: time
  type(ocean_grid_type),   intent(in)    :: g    !< The ocean's grid structure
  type(verticalgrid_type), intent(in)    :: gv   !< The ocean's vertical grid structure
  type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time parameters
  type(diag_ctrl), target, intent(inout) :: diag
  type(energetic_pbl_cs), pointer        :: cs
! Arguments: Time - The current model time.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      diag - A structure that is used to regulate diagnostic output.
!  (in/out)  CS - A pointer that is set to point to the control structure
!                  for this module
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_energetic_PBL"  ! This module's name.
  real :: omega_frac_dflt
  integer :: isd, ied, jsd, jed
  logical :: use_temperature, use_omega
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

  if (associated(cs)) then
    call mom_error(warning, "mixedlayer_init called with an associated"// &
                            "associated control structure.")
    return
  else ; allocate(cs) ; endif

  cs%diag => diag
  cs%Time => time

! Set default, read and log parameters
  call log_version(param_file, mdl, version, "")

  call get_param(param_file, mdl, "MSTAR_MODE", cs%mstar_mode, &
                 "An integer switch for how to compute MSTAR. \n"//&
                 "    0 for constant MSTAR\n"//&
                 "    1 for MSTAR w/ MLD in stabilizing limit\n"//&
                 "    2 for MSTAR w/ L_E/L_O in stabilizing limit.",&
                 "units=nondim",default=0)
  call get_param(param_file, mdl, "MSTAR", cs%mstar, &
                 "The ratio of the friction velocity cubed to the TKE \n"//&
                 "input to the mixed layer.", "units=nondim", default=1.2)
  call get_param(param_file, mdl, "MIX_LEN_EXPONENT", cs%MixLenExponent, &
                 "The exponent applied to the ratio of the distance to the MLD \n"//&
                 "and the MLD depth which determines the shape of the mixing length.",&
                 "units=nondim", default=2.0)
  call get_param(param_file, mdl, "MSTAR_CAP", cs%mstar_cap, &
                 "Maximum value of mstar allowed in model if non-negative\n"//&
                 "(used if MSTAR_MODE>0).",&
                 "units=nondim", default=-1.0)
  call get_param(param_file, mdl, "MSTAR_CONV_ADJ", cs%cnv_mst_fac, &
                 "Factor used for reducing mstar during convection \n"//&
                 " due to reduction of stable density gradient.",&
                 "units=nondim", default=0.0)
  call get_param(param_file, mdl, "MSTAR_SLOPE", cs%mstar_slope, &
                 "The slope of the linear relationship between mstar \n"//&
                 "and the length scale ratio (used if MSTAR_MODE=1).",&
                 "units=nondim", default=0.85)
  call get_param(param_file, mdl, "MSTAR_XINT", cs%mstar_xint, &
                 "The value of the length scale ratio where the mstar \n"//&
                 "is linear above (used if MSTAR_MODE=1).",&
                 "units=nondim", default=-0.3)
  call get_param(param_file, mdl, "MSTAR_AT_XINT", cs%mstar_at_xint, &
                 "The value of mstar at MSTAR_XINT \n"//&
                 "(used if MSTAR_MODE=1).",&
                 "units=nondim", default=0.095)
  call get_param(param_file, mdl, "MSTAR_FLATCAP", cs%MSTAR_FLATCAP, &
                 "Set false to use asymptotic cap, defaults to true.\n"//&
                 "(used only if MSTAR_MODE=1)"&
                 ,"units=nondim",default=.true.)
  call get_param(param_file, mdl, "MSTAR2_COEF1", cs%MSTAR_COEF, &
                 "Coefficient in computing mstar when rotation and \n"//&
                 " stabilizing effects are both important (used if MSTAR_MODE=2)"&
                  ,"units=nondim",default=0.3)
  call get_param(param_file, mdl, "MSTAR2_COEF2", cs%C_EK, &
                 "Coefficient in computing mstar when only rotation limits \n"//&
                 " the total mixing. (used only if MSTAR_MODE=2)"&
                  ,"units=nondim",default=0.085)
  call get_param(param_file, mdl, "NSTAR", cs%nstar, &
                 "The portion of the buoyant potential energy imparted by \n"//&
                 "surface fluxes that is available to drive entrainment \n"//&
                 "at the base of mixed layer when that energy is positive.", &
                 units="nondim", default=0.2)
  call get_param(param_file, mdl, "MKE_TO_TKE_EFFIC", cs%MKE_to_TKE_effic, &
                 "The efficiency with which mean kinetic energy released \n"//&
                 "by mechanically forced entrainment of the mixed layer \n"//&
                 "is converted to turbulent kinetic energy.", units="nondim", &
                 default=0.0)
  call get_param(param_file, mdl, "TKE_DECAY", cs%TKE_decay, &
                 "TKE_DECAY relates the vertical rate of decay of the \n"//&
                 "TKE available for mechanical entrainment to the natural \n"//&
                 "Ekman depth.", units="nondim", default=2.5)
!  call get_param(param_file, mdl, "HMIX_MIN", CS%Hmix_min, &
!                 "The minimum mixed layer depth if the mixed layer depth \n"//&
!                 "is determined dynamically.", units="m", default=0.0)

  call get_param(param_file, mdl, "OMEGA",cs%omega,              &
                 "The rotation rate of the earth.", units="s-1", &
                 default=7.2921e-5)
  call get_param(param_file, mdl, "ML_USE_OMEGA", use_omega,                  &
                 "If true, use the absolute rotation rate instead of the \n"//&
                 "vertical component of rotation when setting the decay \n"// &
                 "scale for turbulence.", default=.false., do_not_log=.true.)
  omega_frac_dflt = 0.0
  if (use_omega) then
    call mom_error(warning, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
    omega_frac_dflt = 1.0
  endif
  call get_param(param_file, mdl, "ML_OMEGA_FRAC", cs%omega_frac,              &
                 "When setting the decay scale for turbulence, use this \n"//  &
                 "fraction of the absolute rotation rate blended with the \n"//&
                 "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).",      &
                 units="nondim", default=omega_frac_dflt)
  call get_param(param_file, mdl, "WSTAR_USTAR_COEF", cs%wstar_ustar_coef,     &
                 "A ratio relating the efficiency with which convectively \n"//&
                 "released energy is converted to a turbulent velocity, \n"//  &
                 "relative to mechanically forced TKE. Making this larger \n"//&
                 "increases the BL diffusivity", units="nondim", default=1.0)
  call get_param(param_file, mdl, "VSTAR_SCALE_FACTOR", cs%vstar_scale_fac, &
                 "An overall nondimensional scaling factor for v*. \n"//    &
                 "Making this larger decreases the PBL diffusivity.",       &
                 units="nondim", default=1.0)
  call get_param(param_file, mdl, "EKMAN_SCALE_COEF", cs%Ekman_scale_coef,           &
                 "A nondimensional scaling factor controlling the inhibition \n"//   &
                 "of the diffusive length scale by rotation. Making this larger \n"//&
                 "decreases the PBL diffusivity.", units="nondim", default=1.0)
  call get_param(param_file, mdl, "USE_MLD_ITERATION", cs%USE_MLD_ITERATION,    &
                 "A logical that specifies whether or not to use the \n"//      &
                 "distance to the bottom of the actively turblent boundary \n"//&
                 "layer to help set the EPBL length scale.", default=.false.)
  call get_param(param_file, mdl, "ORIG_MLD_ITERATION", cs%ORIG_MLD_ITERATION,  &
                 "A logical that specifies whether or not to use the \n"//      &
                 "old method for determining MLD depth in iteration, which \n"//&
                 "is limited to resolution.", default=.true.)
  call get_param(param_file, mdl, "MLD_ITERATION_GUESS", cs%MLD_ITERATION_GUESS,       &
                 "A logical that specifies whether or not to use the \n"//             &
                 "previous timestep MLD as a first guess in the MLD iteration.\n"//    &
                 "The default is false to facilitate reproducibility.", default=.false.)
  call get_param(param_file, mdl, "EPBL_MLD_TOLERANCE", cs%MLD_tol,         &
                 "The tolerance for the iteratively determined mixed \n"//  &
                 "layer depth.  This is only used with USE_MLD_ITERATION.", &
                 units="meter", default=1.0)
  call get_param(param_file, mdl, "EPBL_MIN_MIX_LEN", cs%min_mix_len,    &
                 "The minimum mixing length scale that will be used \n"//&
                 "by ePBL.  The default (0) does not set a minimum.",    &
                 units="meter", default=0.0)
  call get_param(param_file, mdl, "EPBL_ORIGINAL_PE_CALC", cs%orig_PE_calc,         &
                 "If true, the ePBL code uses the original form of the \n"//        &
                 "potential energy change code.  Otherwise, the newer \n"//         &
                 "version that can work with successive increments to the \n"//     &
                 "diffusivity in upward or downward passes is used.", default=.true.)
  call get_param(param_file, mdl, "EPBL_TRANSITION_SCALE", cs%transLay_scale, &
                 "A scale for the mixing length in the transition layer \n"// &
                 "at the edge of the boundary layer as a fraction of the \n"//&
                 "boundary layer thickness.  The default is 0.1.", &
                 units="nondim", default=0.1)
  if ( cs%USE_MLD_ITERATION .and. abs(cs%transLay_scale-0.5).ge.0.5) then
    call mom_error(fatal, "If flag USE_MLD_ITERATION is true, then "//       &
                 "EPBL_TRANSITION should be greater than 0 and less than 1.")
  endif
  call get_param(param_file, mdl, "N2_DISSIPATION_POS", cs%N2_Dissipation_Scale_Pos, &
                 "A scale for the dissipation of TKE due to stratification \n"//     &
                 "in the boundary layer, applied when local stratification \n"//     &
                 "is positive.  The default is 0, but should probably be ~0.4.",     &
                 units="nondim", default=0.0)
  call get_param(param_file, mdl, "N2_DISSIPATION_NEG", cs%N2_Dissipation_Scale_Neg,&
                 "A scale for the dissipation of TKE due to stratification \n"//    &
                 "in the boundary layer, applied when local stratification \n"//    &
                 "is negative.  The default is 0, but should probably be ~1.",      &
                 units="nondim", default=0.0)
   call get_param(param_file, mdl, "USE_LA_LI2016", cs%USE_LA_Windsea,            &
                 "A logical to use the Li et al. 2016 (submitted) formula to \n"//&
                 " determine the Langmuir number.",                               &
                 units="nondim", default=.false.)
   call get_param(param_file, mdl, "LA_DEPTH_RATIO", cs%LaDepthRatio,                &
                 "The depth (normalized by BLD) to average Stokes drift over in \n"//&
                 " Lanmguir number calculation, where La = sqrt(ust/Stokes).",       &
                 units="nondim",default=0.04)
   call get_param(param_file, mdl, "LT_ENHANCE", cs%LT_ENHANCE_FORM,        &
                 "Integer for Langmuir number mode. \n"//                   &
                 " *Requires USE_LA_LI2016 to be set to True. \n"//         &
                 "Options: 0 - No Langmuir \n"//                            &
                 "         1 - Van Roekel et al. 2014/Li et al., 2016  \n"//&
                 "         2 - Multiplied w/ adjusted La. \n"//             &
                 "         3 - Added w/ adjusted La.",                      &
                 units="nondim", default=0)
   call get_param(param_file, mdl, "LT_ENHANCE_COEF", cs%LT_ENHANCE_COEF, &
                 "Coefficient for Langmuir enhancement if LT_ENHANCE > 1",&
                 units="nondim", default=0.447)
   call get_param(param_file, mdl, "LT_ENHANCE_EXP", cs%LT_ENHANCE_EXP, &
                 "Exponent for Langmuir enhancement if LT_ENHANCE > 1", &
                 units="nondim", default=-1.33)
   call get_param(param_file, mdl, "LT_MOD_LAC1", cs%LaC_MLDoEK,             &
                 "Coefficient for modification of Langmuir number due to\n"//&
                 " MLD approaching Ekman depth if LT_ENHANCE=2.",            &
                 units="nondim", default=-0.87)
   call get_param(param_file, mdl, "LT_MOD_LAC2", cs%LaC_MLDoOB_stab,        &
                 "Coefficient for modification of Langmuir number due to\n"//&
                 " MLD approaching stable Obukhov depth if LT_ENHANCE=2.",   &
                  units="nondim", default=0.0)
   call get_param(param_file, mdl, "LT_MOD_LAC3", cs%LaC_MLDoOB_un,          &
                 "Coefficient for modification of Langmuir number due to\n"//&
                 " MLD approaching unstable Obukhov depth if LT_ENHANCE=2.", &
                  units="nondim", default=0.0)
   call get_param(param_file, mdl, "LT_MOD_LAC4", cs%Lac_EKoOB_stab,         &
                 "Coefficient for modification of Langmuir number due to\n"//&
                 " ratio of Ekman to stable Obukhov depth if LT_ENHANCE=2.", &
                  units="nondim", default=0.95)
   call get_param(param_file, mdl, "LT_MOD_LAC5", cs%Lac_EKoOB_un,            &
                 "Coefficient for modification of Langmuir number due to\n"// &
                 " ratio of Ekman to unstable Obukhov depth if LT_ENHANCE=2.",&
                  units="nondim", default=0.95)
   if (cs%LT_ENHANCE_FORM.gt.0 .and. (.not.cs%USE_LA_Windsea)) then
      call mom_error(fatal, "If flag USE_LA_LI2016 is false, then "//&
                 " LT_ENHANCE must be 0.")
  endif
  ! This gives a minimum decay scale that is typically much less than Angstrom.
  cs%ustar_min = 2e-4*cs%omega*(gv%Angstrom_z + gv%H_to_m*gv%H_subroundoff)
  call log_param(param_file, mdl, "EPBL_USTAR_MIN", cs%ustar_min, &
                 "The (tiny) minimum friction velocity used within the \n"//&
                 "ePBL code, derived from OMEGA and ANGSTROM.", units="meter second-1")

  cs%id_ML_depth = register_diag_field('ocean_model', 'ePBL_h_ML', diag%axesT1, &
      time, 'Surface mixed layer depth', 'meter')
  cs%id_TKE_wind = register_diag_field('ocean_model', 'ePBL_TKE_wind', diag%axesT1, &
      time, 'Wind-stirring source of mixed layer TKE', 'meter3 second-3')
  cs%id_TKE_MKE = register_diag_field('ocean_model', 'ePBL_TKE_MKE', diag%axesT1, &
      time, 'Mean kinetic energy source of mixed layer TKE', 'meter3 second-3')
  cs%id_TKE_conv = register_diag_field('ocean_model', 'ePBL_TKE_conv', diag%axesT1, &
      time, 'Convective source of mixed layer TKE', 'meter3 second-3')
  cs%id_TKE_forcing = register_diag_field('ocean_model', 'ePBL_TKE_forcing', diag%axesT1, &
      time, 'TKE consumed by mixing surface forcing or penetrative shortwave radation'//&
            ' through model layers', 'meter3 second-3')
  cs%id_TKE_mixing = register_diag_field('ocean_model', 'ePBL_TKE_mixing', diag%axesT1, &
      time, 'TKE consumed by mixing that deepens the mixed layer', 'meter3 second-3')
  cs%id_TKE_mech_decay = register_diag_field('ocean_model', 'ePBL_TKE_mech_decay', diag%axesT1, &
      time, 'Mechanical energy decay sink of mixed layer TKE', 'meter3 second-3')
  cs%id_TKE_conv_decay = register_diag_field('ocean_model', 'ePBL_TKE_conv_decay', diag%axesT1, &
      time, 'Convective energy decay sink of mixed layer TKE', 'meter3 second-3')
  cs%id_Hsfc_used = register_diag_field('ocean_model', 'ePBL_Hs_used', diag%axesT1, &
      time, 'Surface region thickness that is used', 'meter')
  cs%id_Mixing_Length = register_diag_field('ocean_model', 'Mixing_Length', diag%axesTi, &
      time, 'Mixing Length that is used', 'meter')
  cs%id_Velocity_Scale = register_diag_field('ocean_model', 'Velocity_Scale', diag%axesTi, &
      time, 'Velocity Scale that is used.', 'meter second-1')
  cs%id_LT_enhancement = register_diag_field('ocean_model', 'LT_Enhancement', diag%axesT1, &
      time, 'LT enhancement that is used.', 'non-dim')
  cs%id_MSTAR_mix = register_diag_field('ocean_model', 'MSTAR', diag%axesT1, &
      time, 'MSTAR that is used.', 'non-dim')
  cs%id_OSBL = register_diag_field('ocean_model', 'ePBL_OSBL', diag%axesT1, &
      time, 'Boundary layer depth from the iteration.', 'meter')
  cs%id_mld_ekman = register_diag_field('ocean_model', 'MLD_EKMAN', diag%axesT1, &
      time, 'Boundary layer depth over Ekman length.', 'meter')
  cs%id_mld_obukhov = register_diag_field('ocean_model', 'MLD_OBUKHOV', diag%axesT1, &
      time, 'Boundary layer depth over Obukhov length.', 'meter')
  cs%id_ekman_obukhov = register_diag_field('ocean_model', 'EKMAN_OBUKHOV', diag%axesT1, &
      time, 'Ekman length over Obukhov length.', 'meter')
  cs%id_LA = register_diag_field('ocean_model', 'LA', diag%axesT1, &
      time, 'Langmuir number.', 'non-dim')
  cs%id_LA_mod = register_diag_field('ocean_model', 'LA_MOD', diag%axesT1, &
      time, 'Modified Langmuir number.', 'non-dim')


  call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
                 "If true, temperature and salinity are used as state \n"//&
                 "variables.", default=.true.)

  if (max(cs%id_TKE_wind, cs%id_TKE_MKE, cs%id_TKE_conv, &
          cs%id_TKE_mixing, cs%id_TKE_mech_decay, cs%id_TKE_forcing, &
          cs%id_TKE_conv_decay) > 0) then
    call safe_alloc_alloc(cs%diag_TKE_wind, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%diag_TKE_MKE, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%diag_TKE_conv, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%diag_TKE_forcing, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%diag_TKE_mixing, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%diag_TKE_mech_decay, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%diag_TKE_conv_decay, isd, ied, jsd, jed)

    cs%TKE_diagnostics = .true.
  endif
  if ((cs%id_Mixing_Length>0) .or. (cs%id_Velocity_Scale>0)) then
    call safe_alloc_alloc(cs%Velocity_Scale,isd,ied,jsd,jed,gv%ke+1)
    call safe_alloc_alloc(cs%Mixing_Length,isd,ied,jsd,jed,gv%ke+1)
    cs%Velocity_Scale(:,:,:) = 0.0
    cs%Mixing_Length(:,:,:) = 0.0
    cs%mixing_diagnostics = .true.
  endif
  call safe_alloc_alloc(cs%ML_depth, isd, ied, jsd, jed)
  call safe_alloc_alloc(cs%ML_depth2, isd, ied, jsd, jed)
  if (max(cs%id_LT_Enhancement, cs%id_mstar_mix,cs%id_mld_ekman, &
       cs%id_ekman_obukhov, cs%id_mld_obukhov, cs%id_LA, cs%id_LA_mod)>0) then
    call safe_alloc_alloc(cs%Mstar_mix, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%Enhance_M, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%MLD_EKMAN, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%MLD_OBUKHOV, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%EKMAN_OBUKHOV, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%LA, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%LA_MOD, isd, ied, jsd, jed)
  endif

  !Fitting coefficients to asymptote twoard 0 as MLD -> Ekman depth
  cs%MSTAR_A = cs%MSTAR_AT_XINT**(1./cs%MSTAR_N)
  cs%MSTAR_B = cs%MSTAR_SLOPE / (cs%MSTAR_N*cs%MSTAR_A**(cs%MSTAR_N-1.))
  !Fitting coefficients to asymptote toward MSTAR_CAP
  !*Fixed to begin asymptote at MSTAR_CAP-0.5 toward MSTAR_CAP
  cs%MSTAR_A2 = 0.5**(1./cs%MSTAR_N)
  cs%MSTAR_B2 = -cs%MSTAR_SLOPE / (cs%MSTAR_N*cs%MSTAR_A2**(cs%MSTAR_N-1))
  !Compute value of X (referenced to MSTAR_XINT) where transition
  ! to asymptotic regime based on value of X where MSTAR=MSTAR_CAP-0.5
  cs%MSTAR_XINT_UP = (cs%MSTAR_CAP-0.5-cs%MSTAR_AT_XINT)/cs%MSTAR_SLOPE

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, 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  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_h0	The previously used diffusivity at an interface times the time step and divided by the average of the thicknesses around the interface, in units of H (m or kg-2).
[in]	dkddt_h	The trial change in the diffusivity at an interface times the time step and divided by the average of the thicknesses around the interface, in units of H (m or kg-2).
[in]	hp_a	The 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, in H.
[in]	hp_b	The 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, in H.
[in]	th_a	An effective temperature times a thickness in the layer above, including implicit mixing effects with other yet higher layers, in K H.
[in]	sh_a	An effective salinity times a thickness in the layer above, including implicit mixing effects with other yet higher layers, in K H.
[in]	th_b	An effective temperature times a thickness in the layer below, including implicit mixing effects with other yet lower layers, in K H.
[in]	sh_b	An effective salinity times a thickness in the layer below, including implicit mixing effects with other yet lower layers, in K H.
[in]	dt_to_dpe_a	A 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, in J m-2 K-1.
[in]	ds_to_dpe_a	A 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, in J m-2 ppt-1.
[in]	dt_to_dpe_b	A 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, in J m-2 K-1.
[in]	ds_to_dpe_b	A 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, in J m-2 ppt-1.
[in]	pres	The hydrostatic interface pressure, which is used to relate the changes in column thickness to the energy that is radiated as gravity waves and unavailable to drive mixing, in Pa.
[in]	dt_to_dcolht_a	A 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, in m K-1.
[in]	ds_to_dcolht_a	A 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, in m ppt-1.
[in]	dt_to_dcolht_b	A 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, in m K-1.
[in]	ds_to_dcolht_b	A 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, in m ppt-1.
[out]	pe_chg	The change in column potential energy from applying Kddt_h at the present interface, in J m-2.
[out]	dpec_dkd	The partial derivative of PE_chg with Kddt_h, in units of J m-2 H-1.
[out]	dpe_max	The maximum change in column potential energy that could be realizedd by applying a huge value of Kddt_h at the present interface, in J m-2.
[out]	dpec_dkd_0	The partial derivative of PE_chg with Kddt_h in the limit where Kddt_h = 0, in J m-2 H-1.
[out]	colht_cor	The correction to PE_chg that is made due to a net change in the column height, in J m-2.

Definition at line 1611 of file MOM_energetic_PBL.F90.

Referenced by energetic_pbl().

  real, intent(in)  :: kddt_h0  !< The previously used diffusivity at an interface times
                                !! the time step and  divided by the average of the
                                !! thicknesses around the interface, in units of H (m or kg-2).
  real, intent(in)  :: dkddt_h  !< The trial change in the diffusivity at an interface times
                                !! the time step and  divided by the average of the
                                !! thicknesses around the interface, in units of H (m or kg-2).
  real, intent(in)  :: hp_a     !< The 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, in H.
  real, intent(in)  :: hp_b     !< The 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, in H.
  real, intent(in)  :: th_a     !< An effective temperature times a thickness in the layer
                                !! above, including implicit mixing effects with other
                                !! yet higher layers, in K H.
  real, intent(in)  :: sh_a     !< An effective salinity times a thickness in the layer
                                !! above, including implicit mixing effects with other
                                !! yet higher layers, in K H.
  real, intent(in)  :: th_b     !< An effective temperature times a thickness in the layer
                                !! below, including implicit mixing effects with other
                                !! yet lower layers, in K H.
  real, intent(in)  :: sh_b     !< An effective salinity times a thickness in the layer
                                !! below, including implicit mixing effects with other
                                !! yet lower layers, in K H.
  real, intent(in)  :: dt_to_dpe_a !< A 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, in J m-2 K-1.
  real, intent(in)  :: ds_to_dpe_a !< A 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, in J m-2 ppt-1.
  real, intent(in)  :: dt_to_dpe_b !< A 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, in J m-2 K-1.
  real, intent(in)  :: ds_to_dpe_b !< A 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, in J m-2 ppt-1.
  real, intent(in)  :: pres     !< The hydrostatic interface pressure, which is used to relate
                                !! the changes in column thickness to the energy that is radiated
                                !! as gravity waves and unavailable to drive mixing, in Pa.
  real, intent(in)  :: dt_to_dcolht_a !< A 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, in m K-1.
  real, intent(in)  :: ds_to_dcolht_a !< A 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, in m ppt-1.
  real, intent(in)  :: dt_to_dcolht_b !< A 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, in m K-1.
  real, intent(in)  :: ds_to_dcolht_b !< A 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, in m ppt-1.

  real, optional, intent(out) :: pe_chg   !< The change in column potential energy from applying
                                          !! Kddt_h at the present interface, in J m-2.
  real, optional, intent(out) :: dpec_dkd !< The partial derivative of PE_chg with Kddt_h,
                                          !! in units of J m-2 H-1.
  real, optional, intent(out) :: dpe_max  !< The maximum change in column potential energy that could
                                          !! be realizedd by applying a huge value of Kddt_h at the
                                          !! present interface, in J m-2.
  real, optional, intent(out) :: dpec_dkd_0 !< The partial derivative of PE_chg with Kddt_h in the
                                            !! limit where Kddt_h = 0, in J m-2 H-1.
  real, optional, intent(out) :: colht_cor  !< The correction to PE_chg that is made due to a net
                                            !! change in the column height, in J m-2.

  real :: hps ! The sum of the two effective pivot thicknesses, in H.
  real :: bdt1 ! A product of the two pivot thicknesses plus a diffusive term, in H2.
  real :: dt_c ! The core term in the expressions for the temperature changes, in K H2.
  real :: ds_c ! The core term in the expressions for the salinity changes, in psu H2.
  real :: pec_core ! The diffusivity-independent core term in the expressions
                   ! for the potential energy changes, J m-3.
  real :: colht_core ! The diffusivity-independent core term in the expressions
                     ! for the column height changes, J m-3.
  real :: colht_chg  ! The change in the column height, in m.
  real :: y1   ! A local temporary term, in units of H-3 or H-4 in various contexts.

  !   The expression for the change in potential energy used here is derived
  ! from the expression for the final estimates of the changes in temperature
  ! and salinities, and then extensively manipulated to get it into its most
  ! succint form. The derivation is not necessarily obvious, but it demonstrably
  ! works by comparison with separate calculations of the energy changes after
  ! the tridiagonal solver for the final changes in temperature and salinity are
  ! applied.

  hps = hp_a + hp_b
  bdt1 = hp_a * hp_b + kddt_h0 * hps
  dt_c = hp_a * th_b - hp_b * th_a
  ds_c = hp_a * sh_b - hp_b * sh_a
  pec_core = hp_b * (dt_to_dpe_a * dt_c + ds_to_dpe_a * ds_c) - &
             hp_a * (dt_to_dpe_b * dt_c + ds_to_dpe_b * ds_c)
  colht_core = hp_b * (dt_to_dcolht_a * dt_c + ds_to_dcolht_a * ds_c) - &
               hp_a * (dt_to_dcolht_b * dt_c + ds_to_dcolht_b * ds_c)

  if (present(pe_chg)) then
    ! Find the change in column potential energy due to the change in the
    ! diffusivity at this interface by dKddt_h.
    y1 = dkddt_h / (bdt1 * (bdt1 + dkddt_h * hps))
    pe_chg = pec_core * y1
    colht_chg = colht_core * y1
    if (colht_chg < 0.0) pe_chg = pe_chg - pres * colht_chg
    if (present(colht_cor)) colht_cor = -pres * min(colht_chg, 0.0)
  else if (present(colht_cor)) then
    y1 = dkddt_h / (bdt1 * (bdt1 + dkddt_h * hps))
    colht_cor = -pres * min(colht_core * y1, 0.0)
  endif

  if (present(dpec_dkd)) then
    ! Find the derivative of the potential energy change with dKddt_h.
    y1 = 1.0 / (bdt1 + dkddt_h * hps)**2
    dpec_dkd = pec_core * y1
    colht_chg = colht_core * y1
    if (colht_chg < 0.0) dpec_dkd = dpec_dkd - pres * colht_chg
  endif

  if (present(dpe_max)) then
    ! This expression is the limit of PE_chg for infinite dKddt_h.
    y1 = 1.0 / (bdt1 * hps)
    dpe_max = pec_core * y1
    colht_chg = colht_core * y1
    if (colht_chg < 0.0) dpe_max = dpe_max - pres * colht_chg
  endif

  if (present(dpec_dkd_0)) then
    ! This expression is the limit of dPEc_dKd for dKddt_h = 0.
    y1 = 1.0 / bdt1**2
    dpec_dkd_0 = pec_core * y1
    colht_chg = colht_core * y1
    if (colht_chg < 0.0) dpec_dkd_0 = dpec_dkd_0 - pres * colht_chg
  endif

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, 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_h	The diffusivity at an interface times the time step and divided by the average of the thicknesses around the interface, in units of H (m or kg-2).
[in]	h_k	The thickness of the layer below the interface, in H.
[in]	b_den_1	The 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, in H.
[in]	dte_term	A diffusivity-independent term related to the temperature change in the layer below the interface, in K H.
[in]	dse_term	A diffusivity-independent term related to the salinity change in the layer below the interface, in ppt H.
[in]	dt_km1_t2	A diffusivity-independent term related to the temperature change in the layer above the interface, in K.
[in]	ds_km1_t2	A diffusivity-independent term related to the salinity change in the layer above the interface, in ppt.
[in]	pres	The hydrostatic interface pressure, which is used to relate the changes in column thickness to the energy that is radiated as gravity waves and unavailable to drive mixing, in Pa.
[in]	dt_to_dpe_k	A 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, in J m-2 K-1.
[in]	ds_to_dpe_k	A 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, in J m-2 ppt-1.
[in]	dt_to_dpea	A 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, in J m-2 K-1.
[in]	ds_to_dpea	A 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, in J m-2 ppt-1.
[in]	dt_to_dcolht_k	A 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, in m K-1.
[in]	ds_to_dcolht_k	A 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, in m ppt-1.
[in]	dt_to_dcolhta	A 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, in m K-1.
[in]	ds_to_dcolhta	A 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, in m ppt-1.
[out]	pe_chg	The change in column potential energy from applying Kddt_h at the present interface, in J m-2.
[out]	dpec_dkd	The partial derivative of PE_chg with Kddt_h, in units of J m-2 H-1.
[out]	dpe_max	The maximum change in column potential energy that could be realizedd by applying a huge value of Kddt_h at the present interface, in J m-2.
[out]	dpec_dkd_0	The partial derivative of PE_chg with Kddt_h in the limit where Kddt_h = 0, in J m-2 H-1.

Definition at line 1760 of file MOM_energetic_PBL.F90.

Referenced by energetic_pbl().

  real, intent(in)  :: kddt_h   !< The diffusivity at an interface times the time step and
                                !! divided by the average of the thicknesses around the
                                !! interface, in units of H (m or kg-2).
  real, intent(in)  :: h_k      !< The thickness of the layer below the interface, in H.
  real, intent(in)  :: b_den_1  !< The 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, in H.
  real, intent(in)  :: dte_term !< A diffusivity-independent term related to the
                                !! temperature change in the layer below the interface, in K H.
  real, intent(in)  :: dse_term !< A diffusivity-independent term related to the
                                !! salinity change in the layer below the interface, in ppt H.
  real, intent(in)  :: dt_km1_t2 !< A diffusivity-independent term related to the
                                 !! temperature change in the layer above the interface, in K.
  real, intent(in)  :: ds_km1_t2 !< A diffusivity-independent term related to the
                                 !! salinity change in the layer above the interface, in ppt.
  real, intent(in)  :: pres      !< The hydrostatic interface pressure, which is used to relate
                                 !! the changes in column thickness to the energy that is radiated
                                 !! as gravity waves and unavailable to drive mixing, in Pa.
  real, intent(in)  :: dt_to_dpe_k !< A 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, in J m-2 K-1.
  real, intent(in)  :: ds_to_dpe_k !< A 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, in J m-2 ppt-1.
  real, intent(in)  :: dt_to_dpea !< A 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, in J m-2 K-1.
  real, intent(in)  :: ds_to_dpea !< A 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, in J m-2 ppt-1.
  real, intent(in)  :: dt_to_dcolht_k !< A 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, in m K-1.
  real, intent(in)  :: ds_to_dcolht_k !< A 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, in m ppt-1.
  real, intent(in)  :: dt_to_dcolhta !< A 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, in m K-1.
  real, intent(in)  :: ds_to_dcolhta !< A 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, in m ppt-1.

  real, optional, intent(out) :: pe_chg   !< The change in column potential energy from applying
                                          !! Kddt_h at the present interface, in J m-2.
  real, optional, intent(out) :: dpec_dkd !< The partial derivative of PE_chg with Kddt_h,
                                          !! in units of J m-2 H-1.
  real, optional, intent(out) :: dpe_max  !< The maximum change in column potential energy that could
                                          !! be realizedd by applying a huge value of Kddt_h at the
                                          !! present interface, in J m-2.
  real, optional, intent(out) :: dpec_dkd_0 !< The partial derivative of PE_chg with Kddt_h in the
                                            !! limit where Kddt_h = 0, in J m-2 H-1.

!   This subroutine determines the total potential energy change due to mixing
! at an interface, including all of the implicit effects of the prescribed
! mixing at interfaces above.  Everything here is derived by careful manipulation
! of the robust tridiagonal solvers used for tracers by MOM6.  The results are
! positive for mixing in a stably stratified environment.
!   The comments describing these arguments are for a downward mixing pass, but
! this routine can also be used for an upward pass with the sense of direction
! reversed.

  real :: b1            ! b1 is used by the tridiagonal solver, in H-1.
  real :: b1kd          ! Temporary array (nondim.)
  real :: colht_chg     ! The change in column thickness in m.
  real :: dcolht_max    ! The change in column thickess for infinite diffusivity, in m.
  real :: dcolht_dkd    ! The partial derivative of column thickess with diffusivity, in s m-1.
  real :: dt_k, dt_km1  ! Temporary arrays in K.
  real :: ds_k, ds_km1  ! Temporary arrays in ppt.
  real :: i_kr_denom, dkr_dkd   ! Temporary arrays in H-2 and nondim.
  real :: ddt_k_dkd, ddt_km1_dkd ! Temporary arrays in K H-1.
  real :: dds_k_dkd, dds_km1_dkd ! Temporary arrays in ppt H-1.

  b1 = 1.0 / (b_den_1 + kddt_h)
  b1kd = kddt_h*b1

  ! Start with the temperature change in layer k-1 due to the diffusivity at
  ! interface K without considering the effects of changes in layer k.

  ! Calculate the change in PE due to the diffusion at interface K
  ! if Kddt_h(K+1) = 0.
  i_kr_denom = 1.0 / (h_k*b_den_1 + (b_den_1 + h_k)*kddt_h)

  dt_k = (kddt_h*i_kr_denom) * dte_term
  ds_k = (kddt_h*i_kr_denom) * dse_term

  if (present(pe_chg)) then
    ! Find the change in energy due to diffusion with strength Kddt_h at this interface.
    ! Increment the temperature changes in layer k-1 due the changes in layer k.
    dt_km1 = b1kd * ( dt_k + dt_km1_t2 )
    ds_km1 = b1kd * ( ds_k + ds_km1_t2 )
    pe_chg = (dt_to_dpe_k * dt_k + dt_to_dpea * dt_km1) + &
             (ds_to_dpe_k * ds_k + ds_to_dpea * ds_km1)
    colht_chg = (dt_to_dcolht_k * dt_k + dt_to_dcolhta * dt_km1) + &
                (ds_to_dcolht_k * ds_k + ds_to_dcolhta * ds_km1)
    if (colht_chg < 0.0) pe_chg = pe_chg - pres * colht_chg
  endif

  if (present(dpec_dkd)) then
    ! Find the derivatives of the temperature and salinity changes with Kddt_h.
    dkr_dkd = (h_k*b_den_1) * i_kr_denom**2

    ddt_k_dkd = dkr_dkd * dte_term
    dds_k_dkd = dkr_dkd * dse_term
    ddt_km1_dkd = (b1**2 * b_den_1) * ( dt_k + dt_km1_t2 ) + b1kd * ddt_k_dkd
    dds_km1_dkd = (b1**2 * b_den_1) * ( ds_k + ds_km1_t2 ) + b1kd * dds_k_dkd

    ! Calculate the partial derivative of Pe_chg with Kddt_h.
    dpec_dkd = (dt_to_dpe_k * ddt_k_dkd + dt_to_dpea * ddt_km1_dkd) + &
               (ds_to_dpe_k * dds_k_dkd + ds_to_dpea * dds_km1_dkd)
    dcolht_dkd = (dt_to_dcolht_k * ddt_k_dkd + dt_to_dcolhta * ddt_km1_dkd) + &
                 (ds_to_dcolht_k * dds_k_dkd + ds_to_dcolhta * dds_km1_dkd)
    if (dcolht_dkd < 0.0) dpec_dkd = dpec_dkd - pres * dcolht_dkd
  endif

  if (present(dpe_max)) then
    ! This expression is the limit of PE_chg for infinite Kddt_h.
    dpe_max = (dt_to_dpea * dt_km1_t2 + ds_to_dpea * ds_km1_t2) + &
              ((dt_to_dpe_k + dt_to_dpea) * dte_term + &
               (ds_to_dpe_k + ds_to_dpea) * dse_term) / (b_den_1 + h_k)
    dcolht_max = (dt_to_dcolhta * dt_km1_t2 + ds_to_dcolhta * ds_km1_t2) + &
              ((dt_to_dcolht_k + dt_to_dcolhta) * dte_term + &
               (ds_to_dcolht_k + ds_to_dcolhta) * dse_term) / (b_den_1 + h_k)
    if (dcolht_max < 0.0) dpe_max = dpe_max - pres*dcolht_max
  endif

  if (present(dpec_dkd_0)) then
    ! This expression is the limit of dPEc_dKd for Kddt_h = 0.
    dpec_dkd_0 = (dt_to_dpea * dt_km1_t2 + ds_to_dpea * ds_km1_t2) / (b_den_1) + &
                 (dt_to_dpe_k * dte_term + ds_to_dpe_k * dse_term) / (h_k*b_den_1)
    dcolht_dkd = (dt_to_dcolhta * dt_km1_t2 + ds_to_dcolhta * ds_km1_t2) / (b_den_1) + &
                 (dt_to_dcolht_k * dte_term + ds_to_dcolht_k * dse_term) / (h_k*b_den_1)
    if (dcolht_dkd < 0.0) dpec_dkd_0 = dpec_dkd_0 - pres*dcolht_dkd
  endif

Here is the caller graph for this function:

get_la_windsea()

subroutine mom_energetic_pbl::get_la_windsea ( real, intent(in)  ustar, real, intent(in)  hbl, type(verticalgrid_type), intent(in)  GV, real, intent(out)  LA  )

private

Parameters

[in]

gv

The ocean's vertical grid structure

Definition at line 1962 of file MOM_energetic_PBL.F90.

References ust_2_u10_coare3p5().

Referenced by energetic_pbl().

! Original description:
! This function returns the enhancement factor, given the 10-meter
! wind (m/s), friction velocity (m/s) and the boundary layer depth (m).
! Update (Jan/25):
! Converted from function to subroutine, now returns Langmuir number.
! Computs 10m wind internally, so only ustar and hbl need passed to
! subroutine.
!
! Qing Li, 160606
! BGR port from CVMix to MOM6 Jan/25/2017
! BGR change output to LA from Efactor
! BGR remove u10 input

! Input
  real, intent(in) :: &
       ! water-side surface friction velocity (m/s)
       ustar, &
       ! boundary layer depth (m)
       hbl
  type(verticalgrid_type), intent(in) :: gv   !< The ocean's vertical grid structure
  real, intent(out) :: la
! Local variables
  ! parameters
  real, parameter :: &
       ! ratio of U19.5 to U10 (Holthuijsen, 2007)
       u19p5_to_u10 = 1.075, &
       ! ratio of mean frequency to peak frequency for
       ! Pierson-Moskowitz spectrum (Webb, 2011)
       fm_to_fp = 1.296, &
       ! ratio of surface Stokes drift to U10
       us_to_u10 = 0.0162, &
       ! loss ratio of Stokes transport
       r_loss = 0.667
  real :: us, hm0, fm, fp, vstokes, kphil, kstar
  real :: z0, z0i, r1, r2, r3, r4, tmp, us_sl, lasl_sqr_i
  real :: pi, u10
  pi = 4.0*atan(1.0)
  ! Computing u10 based on u_star and COARE 3.5 relationships
  call ust_2_u10_coare3p5(ustar*sqrt(gv%Rho0/1.225),u10,gv)
  if (u10 .gt. 0.0 .and. ustar .gt. 0.0) then
    ! surface Stokes drift
    us = us_to_u10*u10
    !
    ! significant wave height from Pierson-Moskowitz
    ! spectrum (Bouws, 1998)
    hm0 = 0.0246 *u10**2
    !
    ! peak frequency (PM, Bouws, 1998)
    tmp = 2.0 * pi * u19p5_to_u10 * u10
    fp = 0.877 * gv%g_Earth / tmp
    !
    ! mean frequency
    fm = fm_to_fp * fp
    !
    ! total Stokes transport (a factor r_loss is applied to account
    !  for the effect of directional spreading, multidirectional waves
    !  and the use of PM peak frequency and PM significant wave height
    !  on estimating the Stokes transport)
    vstokes = 0.125 * pi * r_loss * fm * hm0**2
    !
    ! the general peak wavenumber for Phillips' spectrum
    ! (Breivik et al., 2016) with correction of directional spreading
    kphil = 0.176 * us / vstokes
    !
    ! surface layer averaged Stokes dirft with Stokes drift profile
    ! estimated from Phillips' spectrum (Breivik et al., 2016)
    ! the directional spreading effect from Webb and Fox-Kemper, 2015
    ! is also included
    kstar = kphil * 2.56
    ! surface layer
    !z0 = 0.2 * abs(hbl)
    !BGR hbl now adjusted by averaging ratio before function call.
    z0 = abs(hbl)
    z0i = 1.0 / z0
    ! term 1 to 4
    r1 = ( 0.151 / kphil * z0i -0.84 ) &
         * ( 1.0 - exp(-2.0 * kphil * z0) )
    r2 = -( 0.84 + 0.0591 / kphil * z0i ) &
         *sqrt( 2.0 * pi * kphil * z0 ) &
         *erfc( sqrt( 2.0 * kphil * z0 ) )
    r3 = ( 0.0632 / kstar * z0i + 0.125 ) &
         * (1.0 - exp(-2.0 * kstar * z0) )
    r4 = ( 0.125 + 0.0946 / kstar * z0i ) &
         *sqrt( 2.0 * pi *kstar * z0) &
         *erfc( sqrt( 2.0 * kstar * z0 ) )
    us_sl = us * (0.715 + r1 + r2 + r3 + r4)
    !
    la = sqrt(ustar / us_sl)
  else
    la=1.e8
  endif

Here is the call graph for this function:

Here is the caller graph for this function:

ust_2_u10_coare3p5()

subroutine mom_energetic_pbl::ust_2_u10_coare3p5 ( real, intent(in)  USTair, real, intent(out)  U10, type(verticalgrid_type), intent(in)  GV  )

private

Computes wind speed from ustar_air based on COARE 3.5 Cd relationship.

Parameters

[in]

gv

The ocean's vertical grid structure

Definition at line 1920 of file MOM_energetic_PBL.F90.

Referenced by get_la_windsea().

  real, intent(in)  :: ustair
  type(verticalgrid_type), intent(in) :: gv   !< The ocean's vertical grid structure
  real, intent(out) :: u10
  real, parameter :: vonkar = 0.4
  real, parameter :: nu=1e-6
  real :: z0sm, z0, z0rough, u10a, alpha, cd
  integer :: ct

  ! Uses empirical formula for z0 to convert ustar_air to u10 based on the
  !  COARE 3.5 paper (Edson et al., 2013)
  !alpha=m*U10+b
  !Note in Edson et al. 2013, eq. 13 m is given as 0.017.  However,
  ! m=0.0017 reproduces the curve in their figure 6.

  z0sm = 0.11 * nu / ustair; !Compute z0smooth from ustar guess
  u10 = ustair/sqrt(0.001);  !Guess for u10
  u10a = 1000;

  ct=0
  do while (abs(u10a/u10-1.)>0.001)
    ct=ct+1
    u10a = u10
    alpha = min(0.028,0.0017 * u10 - 0.005)
    z0rough = alpha * ustair**2/gv%g_Earth ! Compute z0rough from ustar guess
    z0=z0sm+z0rough
    cd = ( vonkar / log(10/z0) )**2 ! Compute CD from derived roughness
    u10 = ustair/sqrt(cd);!Compute new u10 from derived CD, while loop
                       ! ends and checks for convergence...CT counter
                       ! makes sure loop doesn't run away if function
                       ! doesn't converge.  This code was produced offline
                       ! and converged rapidly (e.g. 2 cycles)
                       ! for ustar=0.0001:0.0001:10.
    if (ct>20) then
      u10 = ustair/sqrt(0.0015) ! I don't expect to get here, but just
                              !  in case it will output a reasonable value.
      exit
    endif
  enddo
  return

Here is the caller graph for this function:

Variable Documentation

max_msg

integer mom_energetic_pbl::max_msg = 2

private

Definition at line 226 of file MOM_energetic_PBL.F90.

num_msg

integer mom_energetic_pbl::num_msg = 0

Definition at line 226 of file MOM_energetic_PBL.F90.

integer :: num_msg = 0, max_msg = 2