Skip to content

Commit d20b989

Browse files
mnlevy1981mnlevy1981
authored
Merge pull request #106 from gustavo-marques/update_stokes_most
Update stokes most
2 parents aaff195 + 971b366 commit d20b989

File tree

3 files changed

+315
-73
lines changed

3 files changed

+315
-73
lines changed

‎src/drivers/cvmix_kpp_drv.F90‎

Lines changed: 16 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -62,7 +62,7 @@ Subroutine cvmix_kpp_driver()
6262
real(cvmix_r8) :: hmix1, hmix5, hmix7, ri_crit, layer_thick1, &
6363
layer_thick4, layer_thick5, layer_thick7, OBL_depth4, &
6464
OBL_depth5, OBL_depth7, N, Nsqr
65-
real(cvmix_r8) :: StokesXI
65+
real(cvmix_r8) :: StokesXI, BEdE_ER, PU_TKE, PS_TKE, PB_TKE
6666
real(cvmix_r8) :: kOBL_depth, Bslope, Vslope
6767
real(cvmix_r8) :: sigma6, OBL_depth6, surf_buoy_force6, surf_fric_vel6, &
6868
vonkarman6, tao, rho0, grav, alpha, Qnonpen, Cp0, &
@@ -690,11 +690,16 @@ Subroutine cvmix_kpp_driver()
690690
vSbar(nlev7), &
691691
source=cvmix_one)
692692
StokesXI = cvmix_one
693+
BEdE_ER = cvmix_one
694+
PU_TKE = cvmix_one
695+
PS_TKE = cvmix_one
696+
PB_TKE = cvmix_one
693697
call cvmix_kpp_compute_StokesXi(zi=zt, &
694698
zk=zw_iface, &
695699
kSL=(nlev7+1)/2, &
696-
SLDepth=cvmix_zero, &
700+
SLDepth=cvmix_one, &
697701
surf_buoy_force=cvmix_one, &
702+
surfBuoy_NS=cvmix_one, &
698703
surf_fric_vel=cvmix_one, &
699704
omega_w2x=cvmix_one, &
700705
uE=ones, &
@@ -703,11 +708,15 @@ Subroutine cvmix_kpp_driver()
703708
vS=vS, &
704709
uSbar=uSbar, &
705710
vSbar=vSbar, &
706-
uS_SLD=cvmix_one, &
707-
vS_SLD=cvmix_one, &
708-
uSbar_SLD=cvmix_one, &
709-
vSbar_SLD=cvmix_one, &
710-
StokesXI=StokesXI)
711+
uS_SL=cvmix_one, &
712+
vS_SL=cvmix_one, &
713+
uSb_SL=cvmix_one, &
714+
vSb_SL=cvmix_one, &
715+
StokesXI=StokesXI, &
716+
BEdE_ER=BEdE_ER, &
717+
PU_TKE=PU_TKE, &
718+
PS_TKE=PS_TKE, &
719+
PB_TKE=PB_TKE)
711720
call cvmix_kpp_compute_OBL_depth(Ri_bulk, zw_iface, OBL_depth7, &
712721
kOBL_depth, zt, surf_fric=cvmix_one, &
713722
surf_buoy=ones, Xi=ones)

‎src/shared/cvmix_kpp.F90‎

Lines changed: 203 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -79,7 +79,7 @@ module cvmix_kpp
7979
public :: cvmix_kpp_compute_bulk_Richardson
8080
public :: cvmix_kpp_compute_turbulent_scales
8181
public :: cvmix_kpp_compute_unresolved_shear
82-
public :: cvmix_kpp_composite_Gshape !STOKES_MOST
82+
public :: cvmix_kpp_composite_Gshape !STOKES_MOST
8383
public :: cvmix_kpp_compute_StokesXi !STOKES_MOST
8484
! These are public for testing, may end up private later
8585
public :: cvmix_kpp_compute_shape_function_coeffs
@@ -88,7 +88,7 @@ module cvmix_kpp
8888
public :: cvmix_kpp_compute_nu_at_OBL_depth_LMD94
8989
public :: cvmix_kpp_EFactor_model
9090
public :: cvmix_kpp_ustokes_SL_model
91-
91+
public :: cvmix_kpp_compute_ER_depth
9292

9393
interface cvmix_coeffs_kpp
9494
module procedure cvmix_coeffs_kpp_low
@@ -3354,7 +3354,7 @@ subroutine cvmix_kpp_composite_Gshape(sigma , Gat1, Gsig, dGdsig)
33543354

33553355
G_1 = MAX( cvmix_zero , MIN( Gat1 , G_m ) )
33563356

3357-
if (sigma .lT. sig_m) then
3357+
if (sigma .lt. sig_m) then
33583358
sig = MAX( sigma , cvmix_zero)
33593359
Gsig = sig + sig * sig * (a2Gsig + a3Gsig * sig)
33603360
dGdsig = cvmix_one + sig * (2.0_cvmix_r8 * a2Gsig + 3.0_cvmix_r8 * a3Gsig * sig)
@@ -3374,78 +3374,73 @@ end subroutine cvmix_kpp_composite_Gshape
33743374
! !IROUTINE: cvmix_coeffs_bkgnd_wrap
33753375
! !INTERFACE:
33763376

3377-
subroutine cvmix_kpp_compute_StokesXi (zi, zk, kSL, SLDepth, &
3378-
surf_buoy_force, surf_fric_vel, omega_w2x, uE, vE, uS, vS, &
3379-
uSbar, vSbar, uS_SLD, vS_SLD, uSbar_SLD, vSbar_SLD, &
3380-
StokesXI, CVmix_kpp_params_user)
3381-
33823377
! !DESCRIPTION:
3383-
! Compute the Stokes similarity parameter, StokesXI, and Entrainment Rule, BEdE\_ER, from
3384-
! surface layer integrated TKE production terms as parameterized in
3385-
! Large et al., 2020 (doi:10.1175/JPO-D-20-0308.1)
3378+
! Compute the Stokes similarity parameter StokesXI, and Entrainment Rule BEdE,
3379+
! and SL integrated tke production terms as
3380+
! parameterized in Large et al., 2021 (doi:10.1175/JPO-D-20-0308.1)
33863381
!\\
33873382
!\\
33883383
3384+
subroutine cvmix_kpp_compute_StokesXi (zi, zk, kSL, SLDepth, surf_buoy_force, &
3385+
surfBuoy_NS,surf_fric_vel, omega_w2x, uE, vE, uS, vS, uSbar, vSbar, &
3386+
uS_SL, vS_SL, uSb_SL, vSb_SL, &
3387+
StokesXI,BEdE_ER,PU_TKE,PS_TKE,PB_TKE,CVmix_kpp_params_user)
3388+
33893389
! !INPUT PARAMETERS:
3390-
real(cvmix_r8), dimension(:), intent(in) :: zi, zk !< Cell interface and center heights < 0
3390+
real(cvmix_r8), dimension(:), intent(in) :: zi, zk !< Cell interface and center heights <= 0
33913391
integer, intent(in) :: kSL !< cell index of Surface Layer Depth
3392-
real(cvmix_r8), intent(in) :: SLDepth !< Surface Layer Depth Integration limit
3393-
real(cvmix_r8), intent(in) :: surf_buoy_force !< Surface buoyancy flux forcing
3392+
real(cvmix_r8), intent(in) :: SLDepth !< Surface Layer Depth > 0
3393+
real(cvmix_r8), intent(in) :: surf_buoy_force,surfBuoy_NS !< Surface buoyancy flux forcing, non-solar
33943394
real(cvmix_r8), intent(in) :: surf_fric_vel, omega_w2x !< Surface wind forcing from x-axis
3395-
real(cvmix_r8), dimension(:), intent(in) :: uE, vE !< Eulerian velocity at centers
3396-
real(cvmix_r8), intent(in) :: uS_SLD, vS_SLD !< Stokes drift at SLDepth
3397-
real(cvmix_r8), intent(in) :: uSbar_SLD, vSbar_SLD !< Average Stokes drift cell kSL to SLDepth
3398-
3395+
real(cvmix_r8), dimension(:), intent(in) :: uE, vE !< Eulerian velocity at cell centers
3396+
real(cvmix_r8), dimension(:), intent(in) :: uS, vS !< Stokes drift at interfaces
3397+
real(cvmix_r8), dimension(:), intent(in) :: uSbar, vSbar !< Cell average Stokes drift
3398+
real(cvmix_r8), intent(in) :: uS_SL, vS_SL !< Stokes drift at SLDepth
3399+
real(cvmix_r8), intent(in) :: uSb_SL, vSb_SL !< Average Stokes drift cell top to SLDepth
33993400
type(cvmix_kpp_params_type), intent(in), optional, target :: CVmix_kpp_params_user
34003401

3401-
! !INPUT/OUTPUT PARAMETERS:
3402-
real(cvmix_r8), dimension(:), intent(inout) :: uS, vS !< Stokes drift at interfaces
3403-
real(cvmix_r8), dimension(:), intent(inout) :: uSbar, vSbar !< Cell average Stokes drift
3402+
! !OUTPUT PARAMETERS:
34043403
real(cvmix_r8), intent(inout) :: StokesXI !< Stokes similarity parameter
3404+
real(cvmix_r8), intent(inout) :: BEdE_ER ! Entrainment rule product [m3 s-3]
3405+
real(cvmix_r8), intent(inout) :: PU_TKE,PS_TKE,PB_TKE ! Shear, Stkes, Buoyancy SL TKE Production [m3 s-3]
34053406

34063407
!EOP
34073408
!BOC
34083409

3410+
! !LOCAL PARAMETERS:
34093411
type(cvmix_kpp_params_type), pointer :: CVmix_kpp_params_in
3410-
real(cvmix_r8) :: PU, PS , PB ! surface layer TKE production terms
3411-
real(cvmix_r8) :: uS_TMP, vS_TMP, uSbar_TMP, vSbar_TMP ! Temporary Store
3412+
real(cvmix_r8), parameter :: CempCGm = 3.5_cvmix_r8 ! Coeff. relating cross-shear mtm flux to sfc stress [nondim]
3413+
real(cvmix_r8), parameter :: CempCGs = 4.7_cvmix_r8 ! Coeff. relating non-local scalar flux to sfc flux [nondim]
3414+
real(cvmix_r8), parameter :: Cu = 0.023_cvmix_r8 ! TKE shear production weight [nondim]
3415+
real(cvmix_r8), parameter :: Cs = 0.038_cvmix_r8 ! TKE Stokes production weight [nondim]
3416+
real(cvmix_r8), parameter :: Cb = 0.96_cvmix_r8 ! TKE buoyancy production weight [nondim]
3417+
real(cvmix_r8) :: PBfact ! Ratio of TKE surface layer production to w*^3 [nondim]
3418+
real(cvmix_r8) :: PU, PS , PB ! Surface layer TKE production terms increments [m3 s-3]
34123419
real(cvmix_r8) :: ustar, delH, delU, delV, omega_E2x, cosOmega, sinOmega
34133420
real(cvmix_r8) :: BLDepth, TauMAG, TauCG, TauDG, taux0, tauy0, Stk0 , Pinc
3414-
real(cvmix_r8) :: PBfact , CempCGm ! Empirical constants
3415-
real(cvmix_r8) :: dtop, tauEtop, tauxtop, tauytop ! Cell top values
3416-
real(cvmix_r8) :: dbot, tauEbot, tauxbot, tauybot, sigbot, Gbot ! Cell bottom values
3417-
integer :: ktmp ! vertical loop
3421+
real(cvmix_r8) :: dtop, tauEtop, tauxtop, tauytop ! Cell top values
3422+
real(cvmix_r8) :: dbot, tauEbot, tauxbot, tauybot, sigbot, Gbot ! Cell bottom values
3423+
integer :: ktmp ! vertical loop
34183424

34193425
CVmix_kpp_params_in => CVmix_kpp_params_saved
34203426
if (present(CVmix_kpp_params_user)) then
34213427
CVmix_kpp_params_in => CVmix_kpp_params_user
34223428
end if
34233429

3424-
if ( CVmix_kpp_params_in%lStokesMOST ) then
3425-
3426-
! Move bottom of cell kSL up to Surface Layer Extent = SLDepth
3427-
uS_TMP = uS(kSL+1)
3428-
vS_TMP = vS(kSL+1)
3429-
uSbar_TMP = uSbar(kSL)
3430-
vSbar_TMP = vSbar(kSL)
3431-
uS(kSL+1) = uS_SLD
3432-
vS(kSL+1) = vS_SLD
3433-
uSbar(kSL)= uSbar_SLD
3434-
vSbar(kSL)= vSbar_SLD
3435-
3436-
CempCGm= 3.5_cvmix_r8
3437-
34383430
ustar = MAX( surf_fric_vel , 1.e-4_cvmix_r8 ) ! > 0
34393431
taux0 = ustar**2 * cos(omega_w2x)
34403432
tauy0 = ustar**2 * sin(omega_w2x)
34413433
Stk0 = sqrt( uS(1)**2 + vS(1)**2 )
3442-
BLDepth = SLDepth / CVmix_kpp_params_in%surf_layer_ext
3434+
BLdepth = SLDepth / CVmix_kpp_params_in%surf_layer_ext
34433435

34443436
! Parameterized Buoyancy production of TKE
3445-
PBfact = 0.110_cvmix_r8
3446-
PB = PBfact * MAX( -surf_buoy_force * BLdepth , cvmix_zero )
3437+
PBfact = 0.0893759_cvmix_r8
3438+
PB = MAX( -surfBuoy_NS * BLdepth * PBfact , cvmix_zero )
3439+
PBfact = 0.00215_cvmix_r8 * CempCGs
3440+
PB = PB + PBfact * BLdepth * ( abs(surf_buoy_force) - surf_buoy_force )
34473441

3448-
! Compute Both Shear Production Terms down from Surface = initial top values
3442+
cosOmega = 0.0
3443+
sinOmega = 0.0
34493444
PU = 0.0
34503445
PS = 0.0
34513446
dtop = 0.0
@@ -3455,18 +3450,19 @@ subroutine cvmix_kpp_compute_StokesXi (zi, zk, kSL, SLDepth, &
34553450
tauxtop = taux0
34563451
tauytop = tauy0
34573452

3458-
do ktmp = 1, kSL
3459-
! SLdepth can be between cell interfaces kSL and kSL+1
3460-
delH = min( max(cvmix_zero, SLdepth - dtop), (zi(ktmp) - zi(ktmp+1) ) )
3461-
dbot = MIN( dtop + delH , SLdepth)
3462-
sigbot = dbot / BLdepth
3463-
Gbot = cvmix_kpp_composite_shape(sigbot)
3464-
TauMAG = ustar * ustar * Gbot / sigbot
3453+
! Integrate Both Shear Production Terms from Surface to top of surface layer cell
3454+
do ktmp = 1, kSL-1
34653455
delU = uE(ktmp) - uE(ktmp+1)
34663456
delV = vE(ktmp) - vE(ktmp+1)
34673457
Omega_E2x= atan2( delV , delU )
34683458
cosOmega = cos(Omega_E2x)
34693459
sinOmega = sin(Omega_E2x)
3460+
3461+
delH = zi(ktmp) - zi(ktmp+1)
3462+
dbot = dtop + delH
3463+
sigbot = dbot / BLdepth
3464+
Gbot = cvmix_kpp_composite_shape(sigbot)
3465+
TauMAG = ustar * ustar * Gbot / sigbot
34703466
tauCG = CempCGm * Gbot * (taux0 * cosOmega - tauy0 * sinOmega)
34713467
! tauDG = sqrt( TauMAG**2 - tauCG**2 ) ! G
34723468
tauDG = TauMAG ! E
@@ -3490,22 +3486,166 @@ subroutine cvmix_kpp_compute_StokesXi (zi, zk, kSL, SLDepth, &
34903486
tauEtop = tauEbot
34913487
enddo
34923488

3493-
! Compute Stokes similarity parameter
3494-
StokesXI = PS / MAX( PU + PS + PB , 1.e-12_cvmix_r8 )
3489+
! Integrate from top of surface layer cell to Surface layer Depth
3490+
delH = SLDepth + zi(ktmp)
3491+
sigbot = CVmix_kpp_params_in%surf_layer_ext
3492+
Gbot = cvmix_kpp_composite_shape(sigbot)
3493+
TauMAG = ustar * ustar * Gbot / sigbot
3494+
tauCG = CempCGm * Gbot * (taux0 * cosOmega - tauy0 * sinOmega)
3495+
! tauDG = sqrt( TauMAG**2 - tauCG**2 ) ! G
3496+
tauDG = TauMAG ! E
3497+
tauxbot = tauDG * cosOmega - tauCG * sinOmega
3498+
tauybot = tauDG * sinOmega + tauCG * cosOmega
3499+
tauEbot = (tauxbot * delU + tauybot * delV) / (zk(ktmp) - zk(ktmp+1) )
3500+
! Increment Eulerian Shear Production
3501+
Pinc = 0.5_cvmix_r8 * (tauEbot + tauEtop) * delH
3502+
PU = PU + MAX( Pinc , cvmix_zero )
3503+
3504+
! Increment Stokes Shear Production
3505+
Pinc = tauxtop*uS(ktmp) - tauxbot*uS_SL + tauytop*vS(ktmp) - tauybot*vS_SL
3506+
Pinc = Pinc - (tauxtop-tauxbot) * uSb_SL - (tauytop-tauybot) * vSb_SL
3507+
PS = PS + MAX( Pinc , cvmix_zero )
3508+
3509+
! Compute Stokes similarity parameter, Entrainment Rule, and TKE prodction Rates
3510+
if ( (PU + PB + PS) > cvmix_zero ) then
3511+
StokesXI = PS / (PU + PS + PB)
3512+
else
3513+
StokesXI = cvmix_zero
3514+
endif
3515+
BEdE_ER = MAX( ( Cu*PU + Cs*PS + Cb*PB ) , cvmix_zero )
3516+
PU_TKE = PU
3517+
PS_TKE = PS
3518+
PB_TKE = PB
3519+
3520+
!EOC
3521+
3522+
end subroutine cvmix_kpp_compute_StokesXi
3523+
3524+
!BOP
3525+
3526+
! !DESCRIPTION:
3527+
! Use Entrainment Rule, BEdE_ER, To Find Entrainment Flux and Depth
3528+
! for each assumed OBL_depth = cell centers,
3529+
! until the boundary layer depth, ERdepth > 0 kER_depth are determined, OR
3530+
! if there is no viable solution ERdepth = -1 , kER_depth=-1
3531+
3532+
subroutine cvmix_kpp_compute_ER_depth( z_inter,Nsq,OBL_depth, &
3533+
uStar,Bsfc_ns,surfBuoy,StokesXI,BEdE_ER,ERdepth, &
3534+
CVMix_kpp_params_user)
3535+
3536+
! !INPUT PARAMETERS:
3537+
real(cvmix_r8), dimension(:), intent(in) :: &
3538+
z_inter, & ! Interface heights <= 0 [m]
3539+
Nsq ! Column of Buoyancy Gradients at interfaces
3540+
real(cvmix_r8), dimension(:), intent(in) :: &
3541+
OBL_depth, & ! Array of assumed OBL depths >0 at cell centers [m]
3542+
surfBuoy, & ! surface Buoyancy flux surface to OBL_depth
3543+
StokesXI, & ! Stokes similarity parameter given OBL_depth
3544+
BEdE_ER ! Parameterized Entrainment Rule given OBL_depth
3545+
real(cvmix_r8), intent(in) :: uStar ! surface friction velocity [m s-1]
3546+
real(cvmix_r8), intent(in) :: Bsfc_ns ! non-solar surface buoyancy flux boundary condition
3547+
3548+
type(cvmix_kpp_params_type), optional, target, intent(in) :: CVmix_kpp_params_user
34953549

3550+
! !OUTPUT PARAMETERS:
3551+
real(cvmix_r8), intent(out) :: ERdepth ! ER Boundary Layer Depth [m]
34963552

3497-
! Restore bottom of cell kSL at zi(kSL+1) with stored Stokes Drift ; ditto average over cell kSL
3498-
uS(kSL+1) = uS_TMP
3499-
vS(kSL+1) = vS_TMP
3500-
uSbar(kSL) = uSbar_TMP
3501-
vSbar(kSL) = vSbar_TMP
3553+
!EOP
3554+
!BOC
35023555

3503-
else ! not lStokesMOST
3504-
StokesXI = cvmix_zero
3505-
end if
3556+
! Local variables
3557+
integer :: iz, nlev , kbl , kinv
3558+
real(cvmix_r8), dimension(size(OBL_depth)+1) :: zdepth, BEdE ! surface then cell-centers<0
3559+
real(cvmix_r8), dimension(size(z_inter)+1) :: sigma, Bflux ! interface values
3560+
real(cvmix_r8) :: ws_i, Cemp_Rs, Gsig_i, Fdelrho, Bnonlocal, sigE, maxNsq, BEnt
3561+
real(kind=cvmix_r8), dimension(4) :: coeffs
3562+
type(cvmix_kpp_params_type), pointer :: CVmix_kpp_params_in
3563+
3564+
CVmix_kpp_params_in => CVmix_kpp_params_saved
3565+
if (present(CVmix_kpp_params_user)) then
3566+
CVmix_kpp_params_in => CVmix_kpp_params_user
3567+
end if
3568+
3569+
nlev = size(OBL_depth)
3570+
Cemp_Rs = 4.7_cvmix_r8
3571+
Fdelrho = cvmix_one
3572+
maxNsq = 0.0
3573+
do kbl = 2, nlev+1
3574+
if ( Nsq(kbl) > maxNsq ) then
3575+
kinv = kbl
3576+
maxNsq = Nsq(kbl)
3577+
endif
3578+
enddo
3579+
3580+
! Set default output values (no solution)
3581+
ERdepth = -cvmix_one
3582+
! Set surface values
3583+
zdepth(1) = cvmix_zero
3584+
BEdE(1) = cvmix_zero
3585+
sigma(:) = cvmix_zero
3586+
Bflux(1) = Bsfc_ns ! non-solar surface buoyancy boundary condition for all kbl
3587+
! Set OBL_depth(1)=top cell center values
3588+
zdepth(2) = -OBL_depth(1)
3589+
BEdE(2) = cvmix_zero
3590+
3591+
do kbl = 2, nlev
3592+
zdepth(kbl+1)= -OBL_depth(kbl)
3593+
BEdE(kbl+1) = cvmix_zero
3594+
3595+
sigma(kbl+1) = cvmix_one
3596+
Bflux(kbl+1) = cvmix_zero
3597+
sigma(kbl+2) = -z_inter(kbl+1)/OBL_depth(kbl) ! > 1.0
3598+
Bflux(kbl+2) = cvmix_zero
3599+
Bnonlocal = 0.5_cvmix_r8 * Cemp_Rs * (abs(surfBuoy(kbl)) - surfBuoy(kbl))
3600+
3601+
do iz = kbl,1,-1
3602+
if (iz .gt. 1) then
3603+
sigma(iz) = -z_inter(iz)/OBL_depth(kbl) ! < 1.0
3604+
call cvmix_kpp_compute_turbulent_scales(sigma(iz), OBL_depth(kbl), surfBuoy(kbl), uStar, StokesXI(kbl), & !0d
3605+
w_s=ws_i , CVmix_kpp_params_user=CVmix_kpp_params_user)
3606+
Gsig_i = cvmix_kpp_composite_shape(sigma(iz))
3607+
Bflux(iz) = Gsig_i * (OBL_depth(kbl) * ws_i * Nsq(iz) - Bnonlocal)
3608+
endif
3609+
3610+
! find the peak
3611+
if ( (Bflux(iz+1) .gt. Bflux(iz+2)) .and. (Bflux(iz+1) .ge. Bflux(iz)) .and. &
3612+
(Bflux(iz+1) .gt. cvmix_zero) ) then
3613+
! Find sigE (the root of the derivative of the quadratic polynomial
3614+
! interpolating (sigma(i), Bflux(i)) for i in [iz, iz+1, iz+2])
3615+
! Also find BEnt (value of quadratic at sigE)
3616+
call cvmix_math_poly_interp(coeffs, sigma(iz:iz+2), Bflux(iz:iz+2))
3617+
! coeffs(3) = 0 => sigma(iz), sigma(iz+1), and sigma(iz+2) are not unique values
3618+
! so there interpolation returned a linear equation. In this case we select
3619+
! (sigma(iz+1), Bflux(iz+1)) as the maximum.
3620+
if (coeffs(3) .eq. cvmix_zero) then
3621+
sigE = sigma(iz+1)
3622+
Bent = Bflux(iz+1)
3623+
else
3624+
sigE = -0.5_cvmix_r8 * (coeffs(2) / coeffs(3))
3625+
Bent = cvmix_math_evaluate_cubic(coeffs, sigE)
3626+
endif
3627+
Fdelrho = cvmix_one
3628+
BEdE(kbl+1) = Fdelrho*BEnt*sigE*OBL_depth(kbl)
3629+
exit ! No exit leaves BEdE(kbl+1) = cvmix_zero
3630+
endif
3631+
enddo
3632+
3633+
if ( (BEdE(kbl+1) > BEdE_ER(kbl)) .and. (Bsfc_ns < cvmix_zero) ) then
3634+
call cvmix_math_poly_interp(coeffs, CVmix_kpp_params_in%interp_type, &
3635+
zdepth(kbl:kbl+1) , BEdE(kbl:kbl+1), &
3636+
zdepth(kbl-1) , BEdE(kbl-1) ) ! surface values if kbl=2;
3637+
coeffs(1) = coeffs(1) - BEdE_ER(kbl)
3638+
ERdepth = -cvmix_math_cubic_root_find(coeffs, 0.5_cvmix_r8 * &
3639+
(zdepth(kbl)+zdepth(kbl+1)))
3640+
3641+
exit
3642+
endif
3643+
3644+
enddo
35063645

35073646
!EOC
35083647

3509-
end subroutine cvmix_kpp_compute_StokesXi
3648+
end subroutine cvmix_kpp_compute_ER_depth
3649+
35103650

35113651
end module cvmix_kpp

0 commit comments

Comments
 (0)