Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
94 changes: 48 additions & 46 deletions physics/PBL/SATMEDMF/satmedmfvdifq.F
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ module satmedmfvdifq
!!
!! Incorporate the TTE-EDMF; if (tte_edmf=.true.),
!! TKE-EDMF scheme becomes TTE-EDMF scheme and the variable 'te'
!! is read as TTE; if (tte_edmf=.false.), the variable 'te' is
!! is read as TTE; if (tte_edmf=.false.), the variable 'te' is
!! read as TKE, 5/22/2025
!!
!!
Expand Down Expand Up @@ -85,25 +85,25 @@ end subroutine satmedmfvdifq_init
!! (mfpbltq.f), is represented using a mass flux approach (Siebesma et al.(2007) \cite Siebesma_2007 ).
!! -# A mass-flux approach is also used to represent the stratocumulus-top-induced turbulence
!! (mfscuq.f).
!! \section detail_satmedmfvidfq GFS satmedmfvdifq Detailed Algorithm
!! \section detail_satmedmfvidfq GFS satmedmfvdifq Detailed Algorithm
subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
& ntiw,ntke,grav,pi,rd,cp,rv,hvap,hfus,fv,eps,epsm1, &
!The following three variables are for SA-3D-TKE
& def_1,def_2,def_3,sa3dtke,dku3d_h,dku3d_e, &
& dv,du,tdt,rtg,u1,v1,t1,q1,usfco,vsfco,icplocn2atm, &
& dv,du,tdt,rtg,u1,v1,t1,q1,usfco,vsfco,use_oceanuv, &
& swh,hlw,xmu,garea,zvfun,sigmaf, &
& psk,rbsoil,zorl,u10m,v10m,fm,fh, &
& tsea,heat,evap,stress,spd1,kpbl, &
& prsi,del,prsl,prslk,phii,phil,delt,tte_edmf, &
& dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl,dkt,dku,tkeh, &
& kinver,xkzm_m,xkzm_h,xkzm_s,dspfac,bl_upfr,bl_dnfr, &
& rlmx,elmx,sfc_rlm,tc_pbl,use_lpt, &
!IVAI: canopy inputs from AQM
!IVAI: canopy inputs from AQM
& do_canopy, cplaqm, claie, cfch, cfrt, cclu, cpopu, &
!IVAI
& ntqv,dtend,dtidx,index_of_temperature,index_of_x_wind, &
& index_of_y_wind,index_of_process_pbl,gen_tend,ldiag3d, &
& errmsg,errflg)
& errmsg,errflg)
!IVAI: aux arrays
! & naux2d,naux3d,aux2d,aux3d)

Expand Down Expand Up @@ -158,7 +158,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
& dtend
integer, intent(in) :: dtidx(:,:), index_of_temperature, &
& index_of_x_wind, index_of_y_wind, index_of_process_pbl
integer, intent(in) :: icplocn2atm
logical, intent(in) :: use_oceanuv
real(kind=kind_phys), intent(out) :: &
& dusfc(:), dvsfc(:), &
& dtsfc(:), dqsfc(:), &
Expand Down Expand Up @@ -192,7 +192,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
integer lcld(im),kcld(im),krad(im),mrad(im)
integer kx1(im), kb1(im), kpblx(im)
!
real(kind=kind_phys) te(im,km), tei(im,km-1), tke(im,km),
real(kind=kind_phys) te(im,km), tei(im,km-1), tke(im,km),
& tteh(im,km), tesq(im,km-1),e2(im,0:km)
!
real(kind=kind_phys) theta(im,km),thvx(im,km), thlvx(im,km),
Expand Down Expand Up @@ -227,7 +227,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
!
real(kind=kind_phys) elm(im,km), ele(im,km),
& ckz(im,km), chz(im,km),
& diss(im,km-1),prod(im,km-1),
& diss(im,km-1),prod(im,km-1),
& bf(im,km-1), shr2(im,km-1), wush(im,km),
& xlamue(im,km-1), xlamde(im,km-1),
& gotvx(im,km), rlam(im,km-1)
Expand Down Expand Up @@ -317,7 +317,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
real(kind=kind_phys) h1

real(kind=kind_phys) bfac, mffac

real(kind=kind_phys) qice(im,km),qliq(im,km)

!PCC_CANOPY------------------------------------
Expand Down Expand Up @@ -438,7 +438,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
km1 = km - 1
kmpbl = km / 2
kmscu = km / 2
!> - Compute physical height of the layer centers and interfaces from
!> - Compute physical height of the layer centers and interfaces from
!! the geopotential height (\p zi and \p zl)
do k=1,km
do i=1,im
Expand Down Expand Up @@ -466,7 +466,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
do i=1,im
gdx(i) = sqrt(garea(i))
enddo
!> - Initialize tke value at vertical layer centers and interfaces
!> - Initialize tke value at vertical layer centers and interfaces
!! from tracer (\p tke and \p tkeh)
do k=1,km
do i=1,im
Expand Down Expand Up @@ -583,11 +583,11 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
kcld(i) = km1
endif
enddo
!
!
!> - Compute a function for green vegetation fraction and surface roughness.
!! Entrainment rate in updraft is a function of vegetation fraction and surface
!! roughness length
!
!
do i = 1,im
tem = (sigmaf(i) - vegflo) / (vegfup - vegflo)
tem = min(max(tem, 0.), 1.)
Expand All @@ -598,7 +598,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
enddo
!
!> - Compute \f$\theta\f$(theta), and \f$q_l\f$(qlx), \f$\theta_e\f$(thetae),
!! \f$\theta_v\f$(thvx),\f$\theta_{l,v}\f$ (thlvx) including ice water
!! \f$\theta_v\f$(thvx),\f$\theta_{l,v}\f$ (thlvx) including ice water
do k=1,km
do i=1,im
pix(i,k) = psk(i) / prslk(i,k)
Expand Down Expand Up @@ -637,7 +637,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
enddo
enddo
!
!> - Compute an empirical cloud fraction based on
!> - Compute an empirical cloud fraction based on
!! Xu and Randall (1996) \cite xu_and_randall_1996
do k = 1, km
do i = 1, im
Expand Down Expand Up @@ -688,7 +688,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!> - Initialize diffusion coefficients to 0 and calculate the total
!> - Initialize diffusion coefficients to 0 and calculate the total
!! radiative heating rate (dku, dkt, radx)
do k=1,km
do i=1,im
Expand All @@ -705,7 +705,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
enddo
!> - Compute stable/unstable PBL flag (pblflg) based on the total
!! surface energy flux (\e false if the total surface energy flux
!! is into the surface)
!! is into the surface)
do i = 1,im
sflux(i) = heat(i) + evap(i)*fv*theta(i,1)
if(.not.sfcflg(i) .or. sflux(i) <= 0.) pblflg(i)=.false.
Expand All @@ -716,12 +716,12 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
!! - Compute critical bulk Richardson number (\f$Rb_{cr}\f$) (crb)
!! - For the unstable PBL, crb is a constant (0.25)
!! - For the stable boundary layer (SBL), \f$Rb_{cr}\f$ varies
!! with the surface Rossby number, \f$R_{0}\f$, as given by
!! with the surface Rossby number, \f$R_{0}\f$, as given by
!! Vickers and Mahrt (2004) \cite Vickers_2004
!! \f[
!! Rb_{cr}=0.16(10^{-7}R_{0})^{-0.18}
!! \f]
!! \f[
!! \f[
!! R_{0}=\frac{U_{10}}{f_{0}z_{0}}
!! \f]
!! where \f$U_{10}\f$ is the wind speed at 10m above the ground surface,
Expand Down Expand Up @@ -753,7 +753,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
ustar(i) = sqrt(stress(i))
enddo
!
!> - Compute buoyancy \f$\frac{\partial \theta_v}{\partial z}\f$ (bf)
!> - Compute buoyancy \f$\frac{\partial \theta_v}{\partial z}\f$ (bf)
!! and the wind shear squared (shr2)
!
do k = 1, km1
Expand Down Expand Up @@ -980,7 +980,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
!
!> - The \f$z/L\f$ (zol) is used as the stability criterion for the PBL.Currently,
!! strong unstable (convective) PBL for \f$z/L < -0.02\f$ and weakly and moderately
!! unstable PBL for \f$0>z/L>-0.02\f$
!! unstable PBL for \f$0>z/L>-0.02\f$
!> - Compute the velocity scale \f$w_s\f$ (wscale) (eqn 22 of Han et al. 2019). It
!! is represented by the value scaled at the top of the surface layer:
!! \f[
Expand Down Expand Up @@ -1172,11 +1172,11 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
do i = 1, im
if(scuflg(i) .and. kcld(i)==km1) scuflg(i)=.false.
enddo
!> - Starting at the PBL top and going downward, if the level is less
!> - Starting at the PBL top and going downward, if the level is less
!! than the cloud top, find the level of the minimum radiative heating
!! rate wihin the cloud. If the level of the minimum is the lowest model
!! level or the minimum radiative heating rate is positive, then set
!! scuflg to F.
!! scuflg to F.
do i = 1, im
flg(i)=scuflg(i)
enddo
Expand All @@ -1202,7 +1202,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> ## Compute components for mass flux mixing by large thermals
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> - If the PBL is convective, the updraft properties are initialized
!> - If the PBL is convective, the updraft properties are initialized
!! to be the same as the state variables.
do k = 1, km
do i = 1, im
Expand Down Expand Up @@ -1236,7 +1236,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
& pcnvflg,zl,zm,q1,t1,u1,v1,plyr,pix,thlx,thvx,
& gdx,hpbl,kpbl,vpert,buou,wush,tkemean,vez0fun,xmf,
& tcko,qcko,ucko,vcko,xlamue,bl_upfr)
!> - Call mfscuq(), which is a new mass-flux parameterization for
!> - Call mfscuq(), which is a new mass-flux parameterization for
!! stratocumulus-top-induced turbulence mixing. For details of the mfscuq subroutine, step into its documentation ::mfscuq
call mfscuq(im,im,km,kmscu,ntcw,ntrac1,dt2,
& scuflg,zl,zm,q1,t1,u1,v1,plyr,pix,
Expand Down Expand Up @@ -1389,7 +1389,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
!! l_2=min(l_{up},l_{down})
!!\f]
!! and dissipation length scale \f$l_d\f$ is given by:
!!\f[
!!\f[
!! l_d=(l_{up}l_{down})^{1/2}
!!\f]
!! where \f$l_{up}\f$ and \f$l_{down}\f$ are the distances that a parcel
Expand All @@ -1413,7 +1413,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
!
enddo
enddo
!> - Compute the surface layer length scale (\f$l_1\f$) following
!> - Compute the surface layer length scale (\f$l_1\f$) following
!! Nakanishi (2001) \cite Nakanish_2001 (eqn 9 of Han et al.(2019) \cite Han_2019)
do k = 1, km1
do i = 1, im
Expand Down Expand Up @@ -1581,7 +1581,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
enddo
enddo
!
! eddy diffusivities at model interface (zm level) in LES scale
! eddy diffusivities at model interface (zm level) in LES scale
!
do k = 1, km1
do i = 1, im
Expand Down Expand Up @@ -1614,7 +1614,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
do k = 2, kmpbl
do i = 1, im
if(tke(i,k) > tkemax(i)) then
tkemax(i) = tke(i,k)
tkemax(i) = tke(i,k)
ktkemax(i) = k
endif
enddo
Expand Down Expand Up @@ -1673,7 +1673,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
enddo
enddo
!
! eddy diffusivities at model layer (zl level) in LES scale
! eddy diffusivities at model layer (zl level) in LES scale
!
do k = 1, km1
do i = 1, im
Expand Down Expand Up @@ -1717,7 +1717,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &

!PCC_CANOPY------------------------------------
kount=0 !IVAI
if (do_canopy .and. cplaqm) then
if (do_canopy .and. cplaqm) then

!IVAI
! print*, 'SATMEDMFVDIFQ_RUN: CLAIE = ', claie(:)
Expand Down Expand Up @@ -1774,7 +1774,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
! & .OR. MAX(0.0, 1.0 - XCANOPYFRT) .GT. 0.5
! & .OR. POPU .GT. 10000.0
! & .OR. EXP(-0.5*XCANOPYLAI*XCANOPYCLU).GT. 0.45
! & .AND. FCH .LT. 18.0 ) THEN
! & .AND. FCH .LT. 18.0 ) THEN

dkt(i,k)= dkt(i,k)
dkq(i,k)= dkq(i,k)
Expand Down Expand Up @@ -1811,9 +1811,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
IF ( ZCAN/FCH .GT. 1.25 ) THEN !SIGMACAN = Eulerian vertical velocity variance
SIGMACAN = 1.25*ustar(i)
END IF
IF ( ZCAN/FCH .GE. 0.175
IF ( ZCAN/FCH .GE. 0.175
& .AND. ZCAN/FCH .LE. 1.25 ) THEN
SIGMACAN = ustar(i) * ( 0.75 +
SIGMACAN = ustar(i) * ( 0.75 +
& (0.5 * COS((PI/1.06818) *
& (1.25 - (ZCAN/FCH)))) )
END IF
Expand All @@ -1825,9 +1825,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
IF ( ZCAN/FCH .GT. 1.25 ) THEN
SIGMACAN = 1.0*ustar(i)
END IF
IF ( ZCAN/FCH .GE. 0.175
IF ( ZCAN/FCH .GE. 0.175
& .AND. ZCAN/FCH .LE. 1.25 ) THEN
SIGMACAN = ustar(i) * ( 0.625 +
SIGMACAN = ustar(i) * ( 0.625 +
& (0.375* COS((PI/1.06818) *
& (1.25 - (ZCAN/FCH)))) )
END IF
Expand All @@ -1839,12 +1839,12 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
IF ( ZCAN/FCH .GT. 1.25 ) THEN
SIGMACAN = 0.25*(4.375 - (3.75*HOL))*ustar(i)
END IF
IF ( ZCAN/FCH .GE. 0.175
IF ( ZCAN/FCH .GE. 0.175
& .AND. ZCAN/FCH .LE. 1.25 ) THEN
RRCAN=4.375-(3.75*HOL)
AACAN=(0.125*RRCAN) + 0.125
BBCAN=(0.125*RRCAN) - 0.125
SIGMACAN = ustar(i) * ( AACAN +
SIGMACAN = ustar(i) * ( AACAN +
& (BBCAN * COS((PI/1.06818) *
& (1.25 - (ZCAN/FCH)))) )
END IF
Expand Down Expand Up @@ -2051,7 +2051,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
enddo
!
!----------------------------------------------------------------------
!> - First predict te due to te production & dissipation(diss)
!> - First predict te due to te production & dissipation(diss)
!
if(sa3dtke) then
!The following is for SA-3D-TKE
Expand Down Expand Up @@ -3051,16 +3051,18 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
! dvsfc(i) = dvsfc(i)+conw*del(i,k)*vtend
enddo
enddo
do i = 1,im
if(icplocn2atm == 0) then
dusfc(i) = -1.*rho_a(i)*stress(i)*u1(i,1)/spd1(i)
dvsfc(i) = -1.*rho_a(i)*stress(i)*v1(i,1)/spd1(i)
else if (icplocn2atm ==1) then
if (use_oceanuv) then
do i = 1,im
spd1_m=sqrt( (u1(i,1)-usfco(i))**2+(v1(i,1)-vsfco(i))**2 )
dusfc(i) = -1.*rho_a(i)*stress(i)*(u1(i,1)-usfco(i))/spd1_m
dvsfc(i) = -1.*rho_a(i)*stress(i)*(v1(i,1)-vsfco(i))/spd1_m
endif
enddo
enddo
else
do i = 1,im
dusfc(i) = -1.*rho_a(i)*stress(i)*u1(i,1)/spd1(i)
dvsfc(i) = -1.*rho_a(i)*stress(i)*v1(i,1)/spd1(i)
enddo
endif
!
if(ldiag3d .and. .not. gen_tend) then
idtend = dtidx(index_of_x_wind,index_of_process_pbl)
Expand Down
8 changes: 4 additions & 4 deletions physics/PBL/SATMEDMF/satmedmfvdifq.meta
Original file line number Diff line number Diff line change
Expand Up @@ -241,12 +241,12 @@
type = real
kind = kind_phys
intent = in
[icplocn2atm]
standard_name = control_for_air_sea_flux_computation_over_water
[use_oceanuv]
standard_name = do_air_sea_flux_computation_over_water
long_name = air-sea flux option
units = 1
units = flag
dimensions = ()
type = integer
type = logical
intent = in
[t1]
standard_name = air_temperature
Expand Down
Loading