Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Issue with MSKF when used with aercu_opt = 1 or 2 #1873

Open
sael9740 opened this issue May 24, 2023 · 8 comments
Open

Issue with MSKF when used with aercu_opt = 1 or 2 #1873

sael9740 opened this issue May 24, 2023 · 8 comments
Assignees

Comments

@sael9740
Copy link
Contributor

While working with the MSKF/Morrison AA combo, I noticed that identical runs (same build, same arch, same number of MPI ranks, etc.) were not producing bit for bit results like I typically expect. This seemed to be the case for aercu_opt = 1 and aercu_opt = 2 but not for aercu_opt = 0. When digging deeper I noticed that QSATu is never initialized at the lower levels even though it is used later on in calculations for other variables:

QSATu initialization

updraft:    DO NK=K,KL-1
...
        QSATu(NK) = QSu/(1.+QSu)   ! saturated specific hum
...
            END DO updraft

QSATu used to calculate su, QSATZM and gamhat

            DO KQ=KTS,KTE
...
                   su(1,JK) = oldTU(KQ)*(1.0+0.622*QSATu(KQ)) + (G*zf_wrf(KQ))/CP !TWG updraft temperature calulation
...
                   QSATZM(1,JK) = QSATu(KQ)
...
                  gamhat(1,JK) = QSATu(KQ)*(1.+QSATu(KQ)/eps1u)   &
                 *eps1u*alatent/(R*oldTU(KQ)**2)*alatent/CP

              END DO

I am thinking that we need to add a QSATu(KQ) = ?? in the following loop to ensure it is initialized at the lower levels below NK=K:

WRF/phys/module_cu_mskf.F

Lines 4717 to 4739 in 21c7214

IF (aercu_opt.gt.0) THEN
zf_wrf(0) = 0.0 ! ground
DO KQ=KTS,KTE
zf_wrf(KQ) = zf_wrf(KQ-1)+DZQ(KQ)
Aqnewlq(kq) = 0.0
Aqnewic(kq) = 0.0
rprd(1,kq) = 0.0
wump(1,kq) =0.0
ncmp(1,kq) =0.0
nimp(1,kq) =0.0
sprd(1,kq) =0.0
frz(1,kq) =0.0
jk = kq
muu(1,JK) = 0.0
duu(1,JK) =0.0
EUU(1,JK) =0.0
cmel(1,JK) =0.0
cmei(1,JK) =0.0
oldTU(kq) = t0(kq)
oldQU(kq) = Q0(kq)
End do
Miter = 0
END IF

I confirmed that adding a QSATu(KQ) = 0.0 here produces B4B results between runs but I'm wondering if the NCAR physics folks can help with what exactly this should be initialized to at the lower levels?

Steps to reproduce the behavior:

  1. Use compiler and version 'nvidia' (probably does the same thing with other compilers)
  2. Use namelist options 'cu_physics = 11, mp_physics = 40, aercu_opt = 1 or 2'

@weiwangncar , @dudhia

@weiwangncar
Copy link
Collaborator

@sael9740 I will take a look next week. Thanks for finding this.

@weiwangncar
Copy link
Collaborator

@sael9740 I tend to agree with you on your proposed fix to initialize qsatu in the loop you identified. Did you look at the model output with and without this fix to see if there is much impact?

@sael9740
Copy link
Contributor Author

sael9740 commented May 30, 2023

@weiwangncar

I tend to agree with you on your proposed fix to initialize qsatu in the loop you identified.

Is setting it to zero the correct thing to do? Or should it be initialized some other way here?

Did you look at the model output with and without this fix to see if there is much impact?

Assuming that we initialize qsat to zero it looks like it does have a significant impact in some cases. FYI - I typically do a three-way comparison of:

  1. WRF built and run with -O3 (with bug fix)
  2. WRF built and run with -O0 (with bug fix)
  3. WRF built and run with -O3

This allows me to understand how much of an impact the changes have relative to the differences that occur due to simple floating-point error propagation, which for our purposes with GPU-WRF are extremely important but I'm explaining this here so that you can understand the plot I'm attaching:

jan2000 1dom 1h mp_physics40 cu_physics11 aercu_opt2 RAINC

The top row plots the fields for each run and the bottom row plots the differences between each pair. I'm seeing that the differences due to the code modification (bottom center plot) are notably larger than the differences between the -O3 and -O0 compilation (bottom left plot).

@weiwangncar
Copy link
Collaborator

@sael9740 I contacted the code contributor, and he suggested to initialize qsatu after line 4365 in v4.5 as

      DO K=1,KX
! initialize Qsatu
      QSATu(K) = 0.0
!
!  Saturation vapor pressure (ES) is calculated following Buck (1981)

Can you see if this works for bit-for-bit results? I think it should. Thanks!

@sael9740
Copy link
Contributor Author

@weiwangncar - I just checked and it looks like adding QSATu(K) = 0.0 in the loop at phys/module_cu_mskf.F line 4365 also produces bit-for-bit results and has very similar effects to what we initially proposed:

jan2000 1dom 1h mp_physics40 cu_physics11 aercu_opt2 RAINC

This does produce different results than what we initially proposed for initializing QSATu but that makes sense after realizing all of the code I mentioned earlier is wrapped in the usl loop

WRF/phys/module_cu_mskf.F

Lines 4424 to 5369 in 21c7214

usl: DO
NU = NU+1
IF(NU.GT.NCHECK)THEN
IF(ISHALL.EQ.1)THEN
CHMAX = 0.
NCHM = 0
DO NK = 1,NCHECK
NNN=KCHECK(NK)
IF(CLDHGT(NNN).GT.CHMAX)THEN
NCHM = NNN
NUCHM = NK
CHMAX = CLDHGT(NNN)
ENDIF
ENDDO
NU = NUCHM-1
FBFRC=1.
CYCLE usl
ELSE
RETURN
ENDIF
ENDIF
KMIX = KCHECK(NU)
LOW=KMIX
!...
LC = LOW
!
!...ASSUME THAT IN ORDER TO SUPPORT A DEEP UPDRAFT YOU NEED A LAYER OF
!...UNSTABLE AIR AT LEAST 50 mb DEEP...TO APPROXIMATE THIS, ISOLATE A
!...GROUP OF ADJACENT INDIVIDUAL MODEL LAYERS, WITH THE BASE AT LEVEL
!...LC, SUCH THAT THE COMBINED DEPTH OF THESE LAYERS IS AT LEAST 50 mb..
!
NLAYRS=0
DPTHMX=0.
NK=LC-1
IF ( NK+1 .LT. KTS ) THEN
WRITE(message,*)'WOULD GO OFF BOTTOM: MSKF_PARA I,J,NK',I,J,NK
CALL wrf_message (TRIM(message))
ELSE
DO
NK=NK+1
IF ( NK .GT. KTE ) THEN
WRITE(message,*)'WOULD GO OFF TOP: MSKF_PARA I,J,DPTHMX,DPMIN',I,J,DPTHMX,DPMIN
CALL wrf_message (TRIM(message))
EXIT
ENDIF
DPTHMX=DPTHMX+DP(NK)
NLAYRS=NLAYRS+1
IF(DPTHMX.GT.DPMIN)THEN
EXIT
ENDIF
END DO
ENDIF
IF(DPTHMX.LT.DPMIN)THEN
RETURN
ENDIF
KPBL=LC+NLAYRS-1
!
!...********************************************************
!...for computational simplicity without much loss in accuracy,
!...mix temperature instead of theta for evaluating convective
!...initiation (triggering) potential...
! THMIX=0.
TMIX=0.
QMIX=0.
ZMIX=0.
PMIX=0.
!
!...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
!...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
!...LAYERS...
!
!cdir novector
DO NK=LC,KPBL
TMIX=TMIX+DP(NK)*T0(NK)
QMIX=QMIX+DP(NK)*Q0(NK)
ZMIX=ZMIX+DP(NK)*Z0(NK)
PMIX=PMIX+DP(NK)*P0(NK)
ENDDO
! THMIX=THMIX/DPTHMX
TMIX=TMIX/DPTHMX
QMIX=QMIX/DPTHMX
ZMIX=ZMIX/DPTHMX
PMIX=PMIX/DPTHMX
EMIX=QMIX*PMIX/(0.622+QMIX)
!
!...FIND THE TEMPERATURE OF THE MIXTURE AT ITS LCL...
!
! TLOG=ALOG(EMIX/ALIQ)
! ...calculate dewpoint using lookup table...
!
astrt=1.e-3
ainc=0.075
a1=emix/aliq
tp=(a1-astrt)/ainc
indlu=int(tp)+1
value=(indlu-1)*ainc+astrt
aintrp=(a1-value)/ainc
tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX-TDPT)
TLCL=AMIN1(TLCL,TMIX)
TVLCL=TLCL*(1.+0.608*QMIX)
ZLCL = ZMIX+(TLCL-TMIX)/GDRY
! NK = LC-1
! DO
! NK = NK+1
! KLCL=NK
! IF(ZLCL.LE.Z0(NK) .or. NK.GT.KL)THEN
! EXIT
! ENDIF
! ENDDO
! IF(NK.GT.KL)THEN
! RETURN
! ENDIF
DO NK = LC, KL
KLCL = NK
IF ( ZLCL.LE.Z0(NK) ) EXIT
END DO
IF ( ZLCL.GT.Z0(KL) ) RETURN
K=KLCL-1
! calculate DLP using Z instead of log(P)
DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
!
!...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
!
TENV=T0(K)+(T0(KLCL)-T0(K))*DLP
QENV=Q0(K)+(Q0(KLCL)-Q0(K))*DLP
TVEN=TENV*(1.+0.608*QENV)
!
! ww: this needs to be initialized
DTRH = 0.
! Bechtold 2001 trigger with my Beta parameter
DTLCL = W0AVG1D(KLCL)/Scale_Fac
if(DTLCL.lt.0.0) then
tempKay = -1.0
DTLCL = tempKay * DTLCL
DTLCL = (DTLCL)**0.3333
else
tempKay = 1.0
DTLCL = tempKay * DTLCL
DTLCL = (DTLCL)**0.3333
end if
DTLCL = 6.0 * tempKay * DTLCL
!
! old trigger
! Stick with the old trigger for now... CGM July 2015
!
IF(ZLCL.LT.2.E3)THEN ! Kain (2004) Eq. 2
WKLCL=0.02*ZLCL/2.E3
ELSE
WKLCL=0.02 ! units of m/s
ENDIF
!TWG.ckay c
if(DX.GE.25.E3) then
WKL=(W0AVG1D(K)+(W0AVG1D(KLCL)-W0AVG1D(K))*DLP)*DX/25.E3-WKLCL
else
WKL=(W0AVG1D(K)+(W0AVG1D(KLCL)-W0AVG1D(K))*DLP)-WKLCL
end if
!TWG ckay, Modified WKL
IF(WKL.LT.0.0001)THEN
DTLCL=0.
ELSE
DTLCL=4.64*WKL**0.33 ! Kain (2004) Eq. 1
ENDIF
! IF(ISHALL.EQ.1)IPRNT=.TRUE.
! IPRNT=.TRUE.
! IF(TLCL+DTLCL.GT.TENV)GOTO 45
IF(TLCL+DTLCL.LT.TENV)THEN
!
! Parcel not buoyant, CYCLE back to start of trigger and evaluate next potential
! USL...
!
CYCLE usl
!
ELSE ! Parcel is buoyant, determine updraft
!
!...CONVECTIVE TRIGGERING CRITERIA HAS BEEN SATISFIED...COMPUTE
!...EQUIVALENT POTENTIAL TEMPERATURE
!...(THETEU) AND VERTICAL VELOCITY OF THE RISING PARCEL AT THE LCL...
!
CALL ENVIRTHT(PMIX,TMIX,QMIX,THETEU(K),ALIQ,BLIQ,CLIQ,DLIQ)
!
!...modify calculation of initial parcel vertical velocity...jsk 11/26/97
!
DTTOT = DTLCL+DTRH
IF(DTTOT.GT.1.E-4)THEN
GDT=2.*G*DTTOT*500./TVEN ! Kain (2004) Eq. 3 (sort of)
WLCL=1.+0.5*SQRT(GDT)
WLCL = AMIN1(WLCL,3.)
ELSE
WLCL=1.
ENDIF
PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
WTW=WLCL*WLCL
!
TVLCL=TLCL*(1.+0.608*QMIX)
RHOLCL=PLCL/(R*TVLCL)
!
LCL=KLCL
LET=LCL
!ckay
! new formulation based on the LCL replacing the cloud radius concept
!introduce LCL instead of RAD based on WKL here
RAD = ZLCL
!ckay Dec20
sourceht = Z0(KPBL)
RAD = amax1(sourceht, RAD)
RAD = AMIN1(4000.,RAD) ! max trap
RAD = AMAX1(500.,RAD) ! min trap
!
!*******************************************************************
! *
! COMPUTE UPDRAFT PROPERTIES *
! *
!*******************************************************************
!
!
!...
!...ESTIMATE INITIAL UPDRAFT MASS FLUX (UMF(K))...
!
WU(K)=WLCL
AU0=0.01*DXSQ
UMF(K)=RHOLCL*AU0
VMFLCL=UMF(K)
UPOLD=VMFLCL
UPNEW=UPOLD
!
!...RATIO2 IS THE DEGREE OF GLACIATION IN THE CLOUD (0 TO 1),
!...UER IS THE ENVIR ENTRAINMENT RATE, ABE IS AVAILABLE
!...BUOYANT ENERGY, TRPPT IS THE TOTAL RATE OF PRECIPITATION
!...PRODUCTION...
!
RATIO2(K)=0.
UER(K)=0.
ABE=0.
TRPPT=0.
TU(K)=TLCL
TVU(K)=TVLCL
QU(K)=QMIX
EQFRC(K)=1.
QLIQ(K)=0.
QICE(K)=0.
IF (aercu_opt .GT. 0) THEN
QRAIN(K)=0.
QSNOW(K)=0.
NLIQ(K)=0.
NICE(K)=0.
NRAIN(K)=0.
NSNOW(K)=0.
CCN(K)=0.
EFFCH(K) = 2.5
EFFIH(K) = 4.99
EFFSH(K) = 9.99
END IF
QLQOUT(K)=0.
QICOUT(K)=0.
DETLQ(K)=0.
DETIC(K)=0.
PPTLIQ(K)=0.
PPTICE(K)=0.
IFLAG=0
!
!...TTEMP IS USED DURING CALCULATION OF THE LINEAR GLACIATION
!...PROCESS; IT IS INITIALLY SET TO THE TEMPERATURE AT WHICH
!...FREEZING IS SPECIFIED TO BEGIN. WITHIN THE GLACIATION
!...INTERVAL, IT IS SET EQUAL TO THE UPDRAFT TEMP AT THE
!...PREVIOUS MODEL LEVEL...
!
TTEMP=TTFRZ
!
!...ENTER THE LOOP FOR UPDRAFT CALCULATIONS...CALCULATE UPDRAFT TEMP,
!...MIXING RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND
!...MOISTURE, PRECIPITATION RATES AT EACH MODEL LEVEL...
!
! **1 variables indicate the bottom of a model layer
! **2 variables indicate the top of a model layer
!
EE1=1.
UD1=0.
REI = 0.
DILBE = 0.
!dkay
IF (aercu_opt.gt.0) THEN
zf_wrf(0) = 0.0 ! ground
DO KQ=KTS,KTE
zf_wrf(KQ) = zf_wrf(KQ-1)+DZQ(KQ)
Aqnewlq(kq) = 0.0
Aqnewic(kq) = 0.0
rprd(1,kq) = 0.0
wump(1,kq) =0.0
ncmp(1,kq) =0.0
nimp(1,kq) =0.0
sprd(1,kq) =0.0
frz(1,kq) =0.0
jk = kq
muu(1,JK) = 0.0
duu(1,JK) =0.0
EUU(1,JK) =0.0
cmel(1,JK) =0.0
cmei(1,JK) =0.0
oldTU(kq) = t0(kq)
oldQU(kq) = Q0(kq)
End do
Miter = 0
END IF
updraft: DO NK=K,KL-1
NK1=NK+1
RATIO2(NK1)=RATIO2(NK)
FRC1=0.
TU(NK1)=T0(NK1)
THETEU(NK1)=THETEU(NK)
QU(NK1)=QU(NK)
!dkay
IF (aercu_opt.gt.0) THEN
oldQU(NK) = QU(NK)
oldTU(NK) = TU(NK)
END IF
!dkay
QLIQ(NK1)=QLIQ(NK)
QICE(NK1)=QICE(NK)
call tpmix2(p0(nk1),theteu(nk1),tu(nk1),qu(nk1),qliq(nk1), &
qice(nk1),qnewlq,qnewic,XLV1,XLV0,QSu)
!dkay QSu has been added to the tpmix2
!dkay
! saturation value of Q of updraft for use with gamma hat in DCCMP routine
! IF (aercu_opt.gt.0) THEN
! QSATu(NK) = QSu/(1.+QSu) ! saturated specific hum
! Aqnewlq(NK) = qnewlq
! Aqnewic(NK) = qnewic
! Aqnewlq(NK) = qnewlq + Qliq(nk ) This is to be removed
! Aqnewic(NK) = qnewic + Qice(nk )
!dkaydec26
! if(TU(NK).le.273.) then
! Aqnewlq(NK) = 0.0
! Aqnewic(NK) = qnewlq + qnewic
! else
! Aqnewlq(NK) = qnewlq + qnewic
! Aqnewic(NK) = 0.0
! end if
! END IF
!
!
!...CHECK TO SEE IF UPDRAFT TEMP IS ABOVE THE TEMPERATURE AT WHICH
!/dec26...GLACIATION IS ASSUMED TO INITIATE; IF IT IS, CALCULATE THE
!...FRACTION OF REMAINING LIQUID WATER TO FREEZE...TTFRZ IS THE
!...TEMP AT WHICH FREEZING BEGINS, TBFRZ THE TEMP BELOW WHICH ALL
!...LIQUID WATER IS FROZEN AT EACH LEVEL...
!
IF(TU(NK1).LE.TTFRZ)THEN
IF(TU(NK1).GT.TBFRZ)THEN
IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
FRC1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
ELSE
FRC1=1.
IFLAG=1
ENDIF
TTEMP=TU(NK1)
!print*,'OLD FRC1',FRC1
!TWG Mar 2017
!Refine FRC1 to match super cooled water path estimate from Hu et al. [2010]
IF (aercu_opt.gt.0) THEN
IF(TU(NK1).GT.TBFRZ)THEN
TC_HU10 = TU(NK1)-273.15
SF_HU10 = -1.0*(P1_HU10+(P2_HU10*TC_HU10)+(P3_HU10*(TC_HU10**2))+ &
(P4_HU10*(TC_HU10**3))+(P5_HU10*(TC_HU10**4))+(P6_HU10*(TC_HU10**5)))
FRC1 = 1.0 - (1.0/(1.0 + EXP(SF_HU10)))
ELSE
FRC1=1.
IFLAG=1
ENDIF
END IF
!END TWG
!
! DETERMINE THE EFFECTS OF LIQUID WATER FREEZING WHEN TEMPERATURE
!...IS BELOW TTFRZ...
!
QFRZ = (QLIQ(NK1)+QNEWLQ)*FRC1
QNEWIC=QNEWIC+QNEWLQ*FRC1
QNEWLQ=QNEWLQ-QNEWLQ*FRC1
QICE(NK1) = QICE(NK1)+QLIQ(NK1)*FRC1
QLIQ(NK1) = QLIQ(NK1)-QLIQ(NK1)*FRC1
CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),QFRZ, &
QICE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
ENDIF
TVU(NK1)=TU(NK1)*(1.+0.608*QU(NK1))
IF (aercu_opt.gt.0) THEN
QSATu(NK) = QSu/(1.+QSu) ! saturated specific hum
Aqnewlq(NK) = qnewlq
Aqnewic(NK) = qnewic
END IF
!
! CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
!
IF(NK.EQ.K)THEN
BE=(TVLCL+TVU(NK1))/(TVEN+TV0(NK1))-1.
BOTERM=2.*(Z0(NK1)-ZLCL)*G*BE/1.5
DZZ=Z0(NK1)-ZLCL
ELSE
BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
BOTERM=2.*DZA(NK)*G*BE/1.5
DZZ=DZA(NK)
ENDIF
ENTERM=2.*REI*WTW/UPOLD
!
! ckay
! using corrected RATE_kay for Test simulation #2... CGM July 2015
!
IF(DX.GE.24.999E3) then
RATE_kay = RATE
else
RATE_kay = RATE / Scale_Fac
end if
CALL CONDLOAD(QLIQ(NK1),QICE(NK1),WTW,DZZ,BOTERM,ENTERM, &
RATE_kay,QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1),G)
!
!...IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND,
!...IF CLOUD IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS...
!
IF(WTW.LT.1.E-3)THEN
EXIT
ELSE
WU(NK1)=SQRT(WTW)
ENDIF
!...Calculate value of THETA-E in environment to entrain into updraft...
!
CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
!
!...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
! New formulation for entrainment
!ckay introduce DX dependcy for the TOKIOKA Parameter =0.03
!ckay Kim et al 2011; Kang et al 2009; Lin et al 2013; GCM findings
TOKIOKA = 0.03
TOKIOKA = TOKIOKA * Scale_Fac
REI=VMFLCL*DP(NK1)*TOKIOKA/RAD
!ckay
TVQU(NK1)=TU(NK1)*(1.+0.608*QU(NK1)-QLIQ(NK1)-QICE(NK1))
IF(NK.EQ.K)THEN
DILBE=((TVLCL+TVQU(NK1))/(TVEN+TV0(NK1))-1.)*DZZ
ELSE
DILBE=((TVQU(NK)+TVQU(NK1))/(TV0(NK)+TV0(NK1))-1.)*DZZ
ENDIF
IF(DILBE.GT.0.)ABE=ABE+DILBE*G
!
!...IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, MINIMAL
!...ENTRAINMENT (0.5*REI) IS IMPOSED...
!
IF(TVQU(NK1).LE.TV0(NK1))THEN ! Entrain/Detrain IF BLOCK
EE2=0.5 ! Kain (2004) Eq. 4
UD2=1.
EQFRC(NK1)=0.
ELSE
LET=NK1
TTMP=TVQU(NK1)
!
!...DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT AND ENVIRONMENTAL AIR...
!
F1=0.95
F2=1.-F1
THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
QTMP=F1*Q0(NK1)+F2*QU(NK1)
TMPLIQ=F2*QLIQ(NK1)
TMPICE=F2*QICE(NK1)
call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice, &
qnewlq,qnewic,XLV1,XLV0,QSu)
TU95=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
IF(TU95.GT.TV0(NK1))THEN
EE2=1.
UD2=0.
EQFRC(NK1)=1.0
ELSE
F1=0.10
F2=1.-F1
THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
QTMP=F1*Q0(NK1)+F2*QU(NK1)
TMPLIQ=F2*QLIQ(NK1)
TMPICE=F2*QICE(NK1)
call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice, &
qnewlq,qnewic,XLV1,XLV0,QSu)
TU10=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
TVDIFF = ABS(TU10-TVQU(NK1))
IF(TVDIFF.LT.1.e-3)THEN
EE2=1.
UD2=0.
EQFRC(NK1)=1.0
ELSE
EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1))
EQFRC(NK1)=AMAX1(0.0,EQFRC(NK1))
EQFRC(NK1)=AMIN1(1.0,EQFRC(NK1))
IF(EQFRC(NK1).EQ.1)THEN
EE2=1.
UD2=0.
ELSEIF(EQFRC(NK1).EQ.0.)THEN
EE2=0.
UD2=1.
ELSE
!
!...SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
! FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES...
!
CALL PROF5(EQFRC(NK1),EE2,UD2)
ENDIF
ENDIF
ENDIF
ENDIF ! End of Entrain/Detrain IF BLOCK
!
!
!...NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE FRACTIONAL
! VALUES IN THE LAYER...
!
EE2 = AMAX1(EE2,0.5)
UD2 = 1.5*UD2
UER(NK1)=0.5*REI*(EE1+EE2)
UDR(NK1)=0.5*REI*(UD1+UD2)
!
!...IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
! UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATIONS...
!
IF(UMF(NK)-UDR(NK1).LT.10.)THEN
!
!...IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL UPD MASS
! FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT THE PREVIOUS MODEL LVL..
! First, correct ABE calculation if needed...
!
IF(DILBE.GT.0.)THEN
ABE=ABE-DILBE*G
ENDIF
LET=NK
! WRITE(98,1015)P0(NK1)/100.
EXIT
ELSE
EE1=EE2
UD1=UD2
UPOLD=UMF(NK)-UDR(NK1)
UPNEW=UPOLD+UER(NK1)
UMF(NK1)=UPNEW
DILFRC(NK1) = UPNEW/UPOLD
!
!...DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND
!...ICE IN THE DETRAINING UPDRAFT MASS...
!
DETLQ(NK1)=QLIQ(NK1)*UDR(NK1)
DETIC(NK1)=QICE(NK1)*UDR(NK1)
QDT(NK1)=QU(NK1)
QU(NK1)=(UPOLD*QU(NK1)+UER(NK1)*Q0(NK1))/UPNEW
THETEU(NK1)=(THETEU(NK1)*UPOLD+THETEE(NK1)*UER(NK1))/UPNEW
QLIQ(NK1)=QLIQ(NK1)*UPOLD/UPNEW
QICE(NK1)=QICE(NK1)*UPOLD/UPNEW
!
!...PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF
!...LIQUID PRECIP AT A GIVEN MODEL LVL, PPTICE THE SAME FOR ICE,
!...TRPPT IS THE TOTAL RATE OF PRODUCTION OF PRECIP UP TO THE
!...CURRENT MODEL LEVEL...
!
PPTLIQ(NK1)=QLQOUT(NK1)*UMF(NK)
PPTICE(NK1)=QICOUT(NK1)*UMF(NK)
!
TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
IF(NK1.LE.KPBL)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
ENDIF
!dkay
IF (aercu_opt.gt.0) THEN
eps1u = 0.622
alatent = 2.54E6
KQ = NK
JK = KTE-KQ+1
muu(1,JK) = UMF(KQ)/VMFLCL ! normalized updraft mass flux
duu(1,JK) = UDR(KQ)/DZQ(KQ)/VMFLCL ! fractional detrainment rate in units of per meter
EUU(1,JK) = UER(KQ)/DZQ(KQ)/VMFLCL ! normalized entrainment rate in unts of per meter
cmel(1,JK) = muu(1,JK)*AQNEWLQ(KQ)/DZQ(KQ)
cmei(1,JK) = muu(1,JK)*AQNEWIC(KQ)/DZQ(KQ)
gamhat(1,JK) = QSATu(KQ)*(1.+QSATu(KQ)/eps1u) &
*eps1u*alatent/(R*oldTU(KQ)**2)*alatent/CP
wu_mskf_act(JK) = WU(KQ) ! kf updraft velocity incloud
qc_mskf_act(JK) = AQNEWLQ(KQ)
qi_mskf_act(JK)=AQNEWIC(KQ)
END IF
!end dkay
!
END DO updraft
!dkay
IF (aercu_opt.gt.0) THEN
Zfu(1,KTE+1) = 0.0
CPin = CP
EPSI0(1) = 2.0E-4
DO KQ=KTS,KTE
JK = KTE-KQ+1
zfu(1,JK) = zf_wrf(KQ)
su(1,JK) = oldTU(KQ)*(1.0+0.622*QSATu(KQ)) + (G*zf_wrf(KQ))/CP !TWG updraft temperature calulation
quu(1,JK) = oldQU(KQ)/(1.+oldQU(KQ)) ! specific humidity of updraft
pru(1,JK) = P0(KQ)/100.0 ! in millibars
TEE(1,JK) = T0(KQ) ! ccc
QEE(1,JK) = Q0(KQ)/(1.+Q0(KQ)) ! specific humidity of environment
qee(1,JK) = oldQU(KQ)/(1.+oldQU(KQ)) ! specific humidity of updraft
QSATZM(1,JK) = QSATu(KQ)
!
!psh: Now, using aerosol concs from CESM
!
denSplume = P0(KQ)/(R*oldTU(KQ))
psh_fac = 1.0E-09/denSplume ! convert ug/m3 to kg/kg
aer_mmr(1,JK, 1) = aercu_fct*aerocu(I,KQ,J, 6)*psh_fac
aer_mmr(1,JK, 2) = aercu_fct*aerocu(I,KQ,J, 5)*psh_fac
aer_mmr(1,JK, 3) = aercu_fct*1.44*aerocu(I,KQ,J, 1)*psh_fac
aer_mmr(1,JK, 4) = aercu_fct*1.44*aerocu(I,KQ,J, 2)*psh_fac
aer_mmr(1,JK, 5) = aercu_fct*1.44*aerocu(I,KQ,J, 3)*psh_fac
aer_mmr(1,JK, 6) = aercu_fct*1.44*aerocu(I,KQ,J, 4)*psh_fac
aer_mmr(1,JK, 7) = aercu_fct*1.54*aerocu(I,KQ,J, 9)*psh_fac
aer_mmr(1,JK, 8) = aercu_fct*1.37*aerocu(I,KQ,J, 7)*psh_fac
aer_mmr(1,JK, 9) = aercu_fct*1.25*aerocu(I,KQ,J,10)*psh_fac
aer_mmr(1,JK,10) = aercu_fct*1.37*aerocu(I,KQ,J, 8)*psh_fac
!psh
gamhat(1,JK) = QSATu(KQ)*(1.+QSATu(KQ)/eps1u) &
*eps1u*alatent/(R*oldTU(KQ)**2)*alatent/CP
END DO
JTT(1) = KX-NK+1
JBB(1) = KX-K+1 ! updraft base level =====>>> flipped for CAM5 indexing
if(jtt(1).gt.jbb(1)) then
JTT(1) = JBB(1)
end if
JLCL(1) = JBB(1) - 1
msg1 = 0
il2g = 1
grav = G
Rdry = R
DTzmp = DT
! print *,'jtt,jbb=', JTT(1), JBB(1)
!dkay: call the new DCCMP scheme here
NLEVZM = KTE-KTS+1 ! this is equal to pver in zm_mp
NLEVZMP1 = NLEVZM + 1 ! pverp
if(jtt(1).eq.1) then
print *,' cloud bottom is on ground!'
print*,'I ',I,' J ',J
CALL wrf_error_fatal ('MSKF Cloud Bottom IS ON THE GROUND, diags' )
end if
if(jbb(1).eq.KTE) then
print *,' cloud top went through the roof!'
print *,'JTT, jbb, jlcl=',JTT(1),JBB(1),JLCL(1)
CALL wrf_error_fatal ( 'MSKF CLOUD TOP WENT OVER MODEL TOP, diags' )
end if
if(DCCMP) then
! do kq=KTE,1,-1
! print *,'wrf dz=',dzq(kq),(KTE-KQ+1)
! end do
call mskf_mphy(su,quu,muu,duu,cmel,cmei,zfu, pru,tee,qee,epsi0, &
jbb,jtt,jlcl, msg1,il2g, grav, cpin, rdry,zmqliq,zmqice,zmqrain,zmqsnow,&
rprd,wump, euu, ncmp,nimp,nrmp,nsmp,zmccn,sprd, frz, aer_mmr, dtzmp, &
NLEVZM,NLEVZMP1,gamhat,qsatzm,wu_mskf_act,qc_mskf_act,qi_mskf_act,effc,effi,effs)
end if
Itest = 0
if(Itest.eq.1) then
write(121,*) 'k,nk, kq,jk,su,quuE3,muu,duu,cmel,zfu,pru,tee,&
&qeeE3,zmqliqE4,zmqiceE4,rprd,wump,euu,ncmp,nimp,sprd,frz'
do kq=K,NK
JK = KTE-KQ+1
write (121,2021) k,nk,kq,jk,su(1,jk),quu(1,jk)*1000,muu(1,jk),duu(1,jk),cmel(1,jk)
write (121,2022) zfu(1,jk),pru(1,jk),tee(1,jk),qee(1,jk)*1000,zmqliq(1,jk)*1.e3,zmqice(1,jk)*1.e3
write (121,2022) rprd(1,jk),wump(1,jk),euu(1,jk),ncmp(1,jk),nimp(1,jk),sprd(1,jk)
write (121,2023) frz(1,jk)
2021 format(4I3,6(1x,E13.6))
2022 format(6(1x,e13.6))
2023 format(2(1x,e13.6))
end do
end if ! itest
if(DCCMP) then
do kq=KTS,KTE
QLIQ(KQ) = 0.0
QICE(KQ) = 0.0
QRAIN(KQ) = 0.0
QSNOW(KQ) = 0.0
NLIQ(KQ) = 0.0
NICE(KQ) = 0.0
NRAIN(KQ) = 0.0
NSNOW(KQ) = 0.0
CCN(KQ) = 0.0
EFFCH(KQ) = 2.51
EFFIH(KQ) = 4.99
EFFSH(KQ) = 9.99
PPTLIQ(KQ)=0.0 ! nov23
PPTICE(KQ)=0.0 ! nov23
QLQOUT(KQ)=0.0 ! nov23
QICOUT(KQ)=0.0 ! nov23
DETLQ(KQ)=0.0 ! dec26
DETIC(KQ)=0.0 ! dec26
end do
TRPPT = 0.0
DO KQ=KTS, KTE
JK = KX-KQ+1
! print *,'kf qliq=', QLIQ(KQ)
QLIQ(KQ) = max(0._r8,zmqliq(1,JK))
QICE(KQ) = max(0._r8,zmqice(1,JK))
!TWG 06/14/16
QRAIN(KQ) = max(0._r8,zmqrain(1,JK))
QSNOW(KQ) = max(0._r8,zmqsnow(1,JK))
NLIQ(KQ) = max(0._r8,ncmp(1,JK))
NICE(KQ) = max(0._r8,nimp(1,JK))
NRAIN(KQ) = max(0._r8,nrmp(1,JK))
NSNOW(KQ) = max(0._r8,nsmp(1,JK))
CCN(KQ) = max(0._r8,zmccn(1,JK))
EFFCH(KQ) = max(2.49_r8, min(effc(1,JK), 50._r8))
EFFIH(KQ) = max(4.99_r8, min(effi(1,JK), 125._r8))
EFFSH(KQ) = max(9.99_r8, min(effs(1,JK), 999._r8))
! END TWG
DETLQ(KQ)= QLIQ(KQ)*UDR(KQ)
DETIC(KQ)= QICE(KQ)*UDR(KQ)
! print *,'zm qliq=', QLIQ(KQ)
densPlume = PPTLIQ(KQ)
!nov23
if(rprd(1,JK).lt.0.0) rprd(1,JK) = 0.0
if(sprd(1,JK).lt.0.0) sprd(1,JK) = 0.0
QLQOUT(KQ)=rprd(1,JK)*dzq(KQ)
QICOUT(KQ)=sprd(1,JK)*dzq(KQ)
PPTLIQ(KQ)=QLQOUT(KQ)*VMFLCL ! check this out
PPTICE(KQ)=QICOUT(KQ)*VMFLCL ! ditto
TRPPT=TRPPT+PPTLIQ(KQ)+PPTICE(KQ)
! if(densPlume.gt.0.0) then
! print *,'zm pptliq=', &
! PPTLIQ(KQ),'kf pptliq=',oldPPTLIQ(kq),'KQ=',KQ
! end if
! print *,'zmQliqout=',kq,densPlume,VMFLCL,pptliq(kq)
! if((i.ge.60.and.i.le.65).and.(j.ge.60.and.j.le.65)) then
! print *, 'KF & MP qliq=', Aqnewlq(nk),QLIQ(NK)
! print *, 'mu & du=', muu(1,JK), duu(1,JK)
! print *,'i,j,k=',I,J,KQ
! end if
END DO
end if ! dccmp
!dkay
2999 CONTINUE
END IF
!
!...CHECK CLOUD DEPTH...IF CLOUD IS TALL ENOUGH, ESTIMATE THE EQUILIBRIUM
! TEMPERATURE LEVEL (LET) AND ADJUST MASS FLUX PROFILE AT CLOUD TOP SO
! THAT MASS FLUX DECREASES TO ZERO AS A LINEAR FUNCTION OF PRESSURE BETWEEN
! THE LET AND CLOUD TOP...
!
!...LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL VELOCITY
! FIRST BECOMES NEGATIVE...
!
LTOP=NK
CLDHGT(LC)=Z0(LTOP)-ZLCL
!
!...Instead of using the same minimum cloud height (for deep convection)
!...everywhere, try specifying minimum cloud depth as a function of TLCL...
!
! Kain (2004) Eq. 7
!
IF(TLCL.GT.293.)THEN
CHMIN = 4.E3
ELSEIF(TLCL.LE.293. .and. TLCL.GE.273)THEN
CHMIN = 2.E3 + 100.*(TLCL-273.)
ELSEIF(TLCL.LT.273.)THEN
CHMIN = 2.E3
ENDIF
!ckay
DO NK=K,LTOP
qc_KF(I,NK,J)=QLIQ(NK)
qi_KF(I,NK,J)=QICE(NK)
! TWG 06/14/16
IF (aercu_opt .GT. 0) THEN
qr_KF(I,NK,J)=QRAIN(NK)
qs_KF(I,NK,J)=QSNOW(NK)
nc_KF(I,NK,J)=NLIQ(NK)
ni_KF(I,NK,J)=NICE(NK)
nr_KF(I,NK,J)=NRAIN(NK)
ns_KF(I,NK,J)=NSNOW(NK)
ccn_KF(I,NK,J)=CCN(NK)
EFCS(I,NK,J)=MAX(2.49, MIN(EFFCH(NK), 50.))
EFIS(I,NK,J)=MAX(4.99, MIN(EFFIH(NK), 120.))
EFSS(I,NK,J)=MAX(9.99, MIN(EFFSH(NK), 999.))
END IF
! END TWG
END DO
!ckay: if mean env RH with respect to water/ice is over 99% then dont allow KF
!ckay: added saturation w.r.to ice june 10, 2015
! to avoid double counting
envRHavg = 0.0
DO NK=K-1,LTOP+1
if(T0(NK).LE.273.16) then
envEsat = 6.112*exp(21.87*(T0(NK)-273.16)/(T0(NK)-7.66))
else
envEsat = 6.112*exp(17.67*(T0(NK)-273.16)/(243.5+T0(NK)-273.16))
end if
envEsat = envEsat * 100.0 ! to hPa
envQsat = 0.622*envEsat/(P0(NK)-envEsat)
envRH = Q0(NK)/envQsat
if(NK.GT.K.and.envRH.LT.0.99) then
envRHavg = 0.0
goto 2020
end if
envRHavg = envRHavg + envRH
END DO
!ckay ; get vertically averaged envRHavg
envRHavg = envRHavg / float(LTOP-K+1+2)
2020 continue
!
!...If cloud top height is less than the specified minimum for deep
!...convection, save value to consider this level as source for
!...shallow convection, go back up to check next level...
!
!...Try specifying minimum cloud depth as a function of TLCL...
!
!
!...DO NOT ALLOW ANY CLOUD FROM THIS LAYER IF:
!
!... 1.) if there is no CAPE, or
!... 2.) cloud top is at model level just above LCL, or
!... 3.) cloud top is within updraft source layer, or
!... 4.) cloud-top detrainment layer begins within
!... updraft source layer.
!...ckay 5.) if the environment is supersaturated i.e., RH > 100%
!...ckay For now, with respect to water
!
IF(LTOP.LE.KLCL .or. LTOP.LE.KPBL .or. LET+1.LE.KPBL &
.or. envRHavg.ge.1.01)THEN ! No Convection Allowed
!ckay
CLDHGT(LC)=0.
DO NK=K,LTOP
UMF(NK)=0.
UDR(NK)=0.
UER(NK)=0.
DETLQ(NK)=0.
DETIC(NK)=0.
PPTLIQ(NK)=0.
PPTICE(NK)=0.
!ckay
cldfra_dp_KF(I,NK,J)=0.
cldfra_sh_KF(I,NK,J)=0.
qc_KF(I,NK,J)=0.
qi_KF(I,NK,J)=0.
!TWG 06/14/16
IF (aercu_opt .GT. 0) THEN
qr_KF(I,NK,J)=0.
qs_KF(I,NK,J)=0.
nc_KF(I,NK,J)=0.
ni_KF(I,NK,J)=0.
nr_KF(I,NK,J)=0.
ns_KF(I,NK,J)=0.
ccn_KF(I,NK,J)=0.
EFCS(I,NK,J)=2.51
EFIS(I,NK,J)=5.01
EFSS(I,NK,J)=10.01
END IF
! END TWG
ENDDO
!
ELSEIF(CLDHGT(LC).GT.CHMIN .and. ABE.GT.1)THEN ! Deep Convection allowed
ISHALL=0
!ckay
DO NK=K,LTOP
cldfra_sh_KF(I,NK,J)=0.
ENDDO
EXIT usl
ELSE
!
!...TO DISALLOW SHALLOW CONVECTION, COMMENT OUT NEXT LINE !!!!!!!!
ISHALL = 1
!ckay
DO NK=K,LTOP
cldfra_dp_KF(I,NK,J)=0.
ENDDO
IF(NU.EQ.NUCHM)THEN
EXIT usl ! Shallow Convection from this layer
ELSE
! Remember this layer (by virtue of non-zero CLDHGT) as potential shallow-cloud layer
DO NK=K,LTOP
UMF(NK)=0.
UDR(NK)=0.
UER(NK)=0.
DETLQ(NK)=0.
DETIC(NK)=0.
PPTLIQ(NK)=0.
PPTICE(NK)=0.
!ckay
cldfra_dp_KF(I,NK,J)=0.
cldfra_sh_KF(I,NK,J)=0.
qc_KF(I,NK,J)=0.
qi_KF(I,NK,J)=0.
!TWG 06/14/16
IF (aercu_opt .GT. 0) THEN
qr_KF(I,NK,J)=0.
qs_KF(I,NK,J)=0.
nc_KF(I,NK,J)=0.
ni_KF(I,NK,J)=0.
nr_KF(I,NK,J)=0.
ns_KF(I,NK,J)=0.
ccn_KF(I,NK,J)=0.
EFCS(I,NK,J)=2.51
EFIS(I,NK,J)=5.01
EFSS(I,NK,J)=10.01
END IF
! END TWG
ENDDO
ENDIF
ENDIF
ENDIF ! for trigger
END DO usl

and the starting index K for the updraft loop

WRF/phys/module_cu_mskf.F

Lines 4741 to 5026 in 21c7214

updraft: DO NK=K,KL-1
NK1=NK+1
RATIO2(NK1)=RATIO2(NK)
FRC1=0.
TU(NK1)=T0(NK1)
THETEU(NK1)=THETEU(NK)
QU(NK1)=QU(NK)
!dkay
IF (aercu_opt.gt.0) THEN
oldQU(NK) = QU(NK)
oldTU(NK) = TU(NK)
END IF
!dkay
QLIQ(NK1)=QLIQ(NK)
QICE(NK1)=QICE(NK)
call tpmix2(p0(nk1),theteu(nk1),tu(nk1),qu(nk1),qliq(nk1), &
qice(nk1),qnewlq,qnewic,XLV1,XLV0,QSu)
!dkay QSu has been added to the tpmix2
!dkay
! saturation value of Q of updraft for use with gamma hat in DCCMP routine
! IF (aercu_opt.gt.0) THEN
! QSATu(NK) = QSu/(1.+QSu) ! saturated specific hum
! Aqnewlq(NK) = qnewlq
! Aqnewic(NK) = qnewic
! Aqnewlq(NK) = qnewlq + Qliq(nk ) This is to be removed
! Aqnewic(NK) = qnewic + Qice(nk )
!dkaydec26
! if(TU(NK).le.273.) then
! Aqnewlq(NK) = 0.0
! Aqnewic(NK) = qnewlq + qnewic
! else
! Aqnewlq(NK) = qnewlq + qnewic
! Aqnewic(NK) = 0.0
! end if
! END IF
!
!
!...CHECK TO SEE IF UPDRAFT TEMP IS ABOVE THE TEMPERATURE AT WHICH
!/dec26...GLACIATION IS ASSUMED TO INITIATE; IF IT IS, CALCULATE THE
!...FRACTION OF REMAINING LIQUID WATER TO FREEZE...TTFRZ IS THE
!...TEMP AT WHICH FREEZING BEGINS, TBFRZ THE TEMP BELOW WHICH ALL
!...LIQUID WATER IS FROZEN AT EACH LEVEL...
!
IF(TU(NK1).LE.TTFRZ)THEN
IF(TU(NK1).GT.TBFRZ)THEN
IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
FRC1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
ELSE
FRC1=1.
IFLAG=1
ENDIF
TTEMP=TU(NK1)
!print*,'OLD FRC1',FRC1
!TWG Mar 2017
!Refine FRC1 to match super cooled water path estimate from Hu et al. [2010]
IF (aercu_opt.gt.0) THEN
IF(TU(NK1).GT.TBFRZ)THEN
TC_HU10 = TU(NK1)-273.15
SF_HU10 = -1.0*(P1_HU10+(P2_HU10*TC_HU10)+(P3_HU10*(TC_HU10**2))+ &
(P4_HU10*(TC_HU10**3))+(P5_HU10*(TC_HU10**4))+(P6_HU10*(TC_HU10**5)))
FRC1 = 1.0 - (1.0/(1.0 + EXP(SF_HU10)))
ELSE
FRC1=1.
IFLAG=1
ENDIF
END IF
!END TWG
!
! DETERMINE THE EFFECTS OF LIQUID WATER FREEZING WHEN TEMPERATURE
!...IS BELOW TTFRZ...
!
QFRZ = (QLIQ(NK1)+QNEWLQ)*FRC1
QNEWIC=QNEWIC+QNEWLQ*FRC1
QNEWLQ=QNEWLQ-QNEWLQ*FRC1
QICE(NK1) = QICE(NK1)+QLIQ(NK1)*FRC1
QLIQ(NK1) = QLIQ(NK1)-QLIQ(NK1)*FRC1
CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),QFRZ, &
QICE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
ENDIF
TVU(NK1)=TU(NK1)*(1.+0.608*QU(NK1))
IF (aercu_opt.gt.0) THEN
QSATu(NK) = QSu/(1.+QSu) ! saturated specific hum
Aqnewlq(NK) = qnewlq
Aqnewic(NK) = qnewic
END IF
!
! CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
!
IF(NK.EQ.K)THEN
BE=(TVLCL+TVU(NK1))/(TVEN+TV0(NK1))-1.
BOTERM=2.*(Z0(NK1)-ZLCL)*G*BE/1.5
DZZ=Z0(NK1)-ZLCL
ELSE
BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
BOTERM=2.*DZA(NK)*G*BE/1.5
DZZ=DZA(NK)
ENDIF
ENTERM=2.*REI*WTW/UPOLD
!
! ckay
! using corrected RATE_kay for Test simulation #2... CGM July 2015
!
IF(DX.GE.24.999E3) then
RATE_kay = RATE
else
RATE_kay = RATE / Scale_Fac
end if
CALL CONDLOAD(QLIQ(NK1),QICE(NK1),WTW,DZZ,BOTERM,ENTERM, &
RATE_kay,QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1),G)
!
!...IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND,
!...IF CLOUD IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS...
!
IF(WTW.LT.1.E-3)THEN
EXIT
ELSE
WU(NK1)=SQRT(WTW)
ENDIF
!...Calculate value of THETA-E in environment to entrain into updraft...
!
CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
!
!...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
! New formulation for entrainment
!ckay introduce DX dependcy for the TOKIOKA Parameter =0.03
!ckay Kim et al 2011; Kang et al 2009; Lin et al 2013; GCM findings
TOKIOKA = 0.03
TOKIOKA = TOKIOKA * Scale_Fac
REI=VMFLCL*DP(NK1)*TOKIOKA/RAD
!ckay
TVQU(NK1)=TU(NK1)*(1.+0.608*QU(NK1)-QLIQ(NK1)-QICE(NK1))
IF(NK.EQ.K)THEN
DILBE=((TVLCL+TVQU(NK1))/(TVEN+TV0(NK1))-1.)*DZZ
ELSE
DILBE=((TVQU(NK)+TVQU(NK1))/(TV0(NK)+TV0(NK1))-1.)*DZZ
ENDIF
IF(DILBE.GT.0.)ABE=ABE+DILBE*G
!
!...IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, MINIMAL
!...ENTRAINMENT (0.5*REI) IS IMPOSED...
!
IF(TVQU(NK1).LE.TV0(NK1))THEN ! Entrain/Detrain IF BLOCK
EE2=0.5 ! Kain (2004) Eq. 4
UD2=1.
EQFRC(NK1)=0.
ELSE
LET=NK1
TTMP=TVQU(NK1)
!
!...DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT AND ENVIRONMENTAL AIR...
!
F1=0.95
F2=1.-F1
THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
QTMP=F1*Q0(NK1)+F2*QU(NK1)
TMPLIQ=F2*QLIQ(NK1)
TMPICE=F2*QICE(NK1)
call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice, &
qnewlq,qnewic,XLV1,XLV0,QSu)
TU95=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
IF(TU95.GT.TV0(NK1))THEN
EE2=1.
UD2=0.
EQFRC(NK1)=1.0
ELSE
F1=0.10
F2=1.-F1
THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
QTMP=F1*Q0(NK1)+F2*QU(NK1)
TMPLIQ=F2*QLIQ(NK1)
TMPICE=F2*QICE(NK1)
call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice, &
qnewlq,qnewic,XLV1,XLV0,QSu)
TU10=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
TVDIFF = ABS(TU10-TVQU(NK1))
IF(TVDIFF.LT.1.e-3)THEN
EE2=1.
UD2=0.
EQFRC(NK1)=1.0
ELSE
EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1))
EQFRC(NK1)=AMAX1(0.0,EQFRC(NK1))
EQFRC(NK1)=AMIN1(1.0,EQFRC(NK1))
IF(EQFRC(NK1).EQ.1)THEN
EE2=1.
UD2=0.
ELSEIF(EQFRC(NK1).EQ.0.)THEN
EE2=0.
UD2=1.
ELSE
!
!...SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
! FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES...
!
CALL PROF5(EQFRC(NK1),EE2,UD2)
ENDIF
ENDIF
ENDIF
ENDIF ! End of Entrain/Detrain IF BLOCK
!
!
!...NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE FRACTIONAL
! VALUES IN THE LAYER...
!
EE2 = AMAX1(EE2,0.5)
UD2 = 1.5*UD2
UER(NK1)=0.5*REI*(EE1+EE2)
UDR(NK1)=0.5*REI*(UD1+UD2)
!
!...IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
! UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATIONS...
!
IF(UMF(NK)-UDR(NK1).LT.10.)THEN
!
!...IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL UPD MASS
! FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT THE PREVIOUS MODEL LVL..
! First, correct ABE calculation if needed...
!
IF(DILBE.GT.0.)THEN
ABE=ABE-DILBE*G
ENDIF
LET=NK
! WRITE(98,1015)P0(NK1)/100.
EXIT
ELSE
EE1=EE2
UD1=UD2
UPOLD=UMF(NK)-UDR(NK1)
UPNEW=UPOLD+UER(NK1)
UMF(NK1)=UPNEW
DILFRC(NK1) = UPNEW/UPOLD
!
!...DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND
!...ICE IN THE DETRAINING UPDRAFT MASS...
!
DETLQ(NK1)=QLIQ(NK1)*UDR(NK1)
DETIC(NK1)=QICE(NK1)*UDR(NK1)
QDT(NK1)=QU(NK1)
QU(NK1)=(UPOLD*QU(NK1)+UER(NK1)*Q0(NK1))/UPNEW
THETEU(NK1)=(THETEU(NK1)*UPOLD+THETEE(NK1)*UER(NK1))/UPNEW
QLIQ(NK1)=QLIQ(NK1)*UPOLD/UPNEW
QICE(NK1)=QICE(NK1)*UPOLD/UPNEW
!
!...PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF
!...LIQUID PRECIP AT A GIVEN MODEL LVL, PPTICE THE SAME FOR ICE,
!...TRPPT IS THE TOTAL RATE OF PRODUCTION OF PRECIP UP TO THE
!...CURRENT MODEL LEVEL...
!
PPTLIQ(NK1)=QLQOUT(NK1)*UMF(NK)
PPTICE(NK1)=QICOUT(NK1)*UMF(NK)
!
TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
IF(NK1.LE.KPBL)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
ENDIF
!dkay
IF (aercu_opt.gt.0) THEN
eps1u = 0.622
alatent = 2.54E6
KQ = NK
JK = KTE-KQ+1
muu(1,JK) = UMF(KQ)/VMFLCL ! normalized updraft mass flux
duu(1,JK) = UDR(KQ)/DZQ(KQ)/VMFLCL ! fractional detrainment rate in units of per meter
EUU(1,JK) = UER(KQ)/DZQ(KQ)/VMFLCL ! normalized entrainment rate in unts of per meter
cmel(1,JK) = muu(1,JK)*AQNEWLQ(KQ)/DZQ(KQ)
cmei(1,JK) = muu(1,JK)*AQNEWIC(KQ)/DZQ(KQ)
gamhat(1,JK) = QSATu(KQ)*(1.+QSATu(KQ)/eps1u) &
*eps1u*alatent/(R*oldTU(KQ)**2)*alatent/CP
wu_mskf_act(JK) = WU(KQ) ! kf updraft velocity incloud
qc_mskf_act(JK) = AQNEWLQ(KQ)
qi_mskf_act(JK)=AQNEWIC(KQ)
END IF
!end dkay
!
END DO updraft

can be different for different iterations of the usl loop. Therefore by initializing QSATu prior to the ust loop the assignments in

WRF/phys/module_cu_mskf.F

Lines 5035 to 5067 in 21c7214

DO KQ=KTS,KTE
JK = KTE-KQ+1
zfu(1,JK) = zf_wrf(KQ)
su(1,JK) = oldTU(KQ)*(1.0+0.622*QSATu(KQ)) + (G*zf_wrf(KQ))/CP !TWG updraft temperature calulation
quu(1,JK) = oldQU(KQ)/(1.+oldQU(KQ)) ! specific humidity of updraft
pru(1,JK) = P0(KQ)/100.0 ! in millibars
TEE(1,JK) = T0(KQ) ! ccc
QEE(1,JK) = Q0(KQ)/(1.+Q0(KQ)) ! specific humidity of environment
qee(1,JK) = oldQU(KQ)/(1.+oldQU(KQ)) ! specific humidity of updraft
QSATZM(1,JK) = QSATu(KQ)
!
!psh: Now, using aerosol concs from CESM
!
denSplume = P0(KQ)/(R*oldTU(KQ))
psh_fac = 1.0E-09/denSplume ! convert ug/m3 to kg/kg
aer_mmr(1,JK, 1) = aercu_fct*aerocu(I,KQ,J, 6)*psh_fac
aer_mmr(1,JK, 2) = aercu_fct*aerocu(I,KQ,J, 5)*psh_fac
aer_mmr(1,JK, 3) = aercu_fct*1.44*aerocu(I,KQ,J, 1)*psh_fac
aer_mmr(1,JK, 4) = aercu_fct*1.44*aerocu(I,KQ,J, 2)*psh_fac
aer_mmr(1,JK, 5) = aercu_fct*1.44*aerocu(I,KQ,J, 3)*psh_fac
aer_mmr(1,JK, 6) = aercu_fct*1.44*aerocu(I,KQ,J, 4)*psh_fac
aer_mmr(1,JK, 7) = aercu_fct*1.54*aerocu(I,KQ,J, 9)*psh_fac
aer_mmr(1,JK, 8) = aercu_fct*1.37*aerocu(I,KQ,J, 7)*psh_fac
aer_mmr(1,JK, 9) = aercu_fct*1.25*aerocu(I,KQ,J,10)*psh_fac
aer_mmr(1,JK,10) = aercu_fct*1.37*aerocu(I,KQ,J, 8)*psh_fac
!psh
gamhat(1,JK) = QSATu(KQ)*(1.+QSATu(KQ)/eps1u) &
*eps1u*alatent/(R*oldTU(KQ)**2)*alatent/CP
END DO

may be using some QSATu values that were initialized in a previous iteration of the usl loop.

My question now though is that if this is the contributor's intention for QSATu then why are we treating variables such as Aqnewlq, Aqnewic, etc. differently by initializing them here:

WRF/phys/module_cu_mskf.F

Lines 4717 to 4739 in 21c7214

IF (aercu_opt.gt.0) THEN
zf_wrf(0) = 0.0 ! ground
DO KQ=KTS,KTE
zf_wrf(KQ) = zf_wrf(KQ-1)+DZQ(KQ)
Aqnewlq(kq) = 0.0
Aqnewic(kq) = 0.0
rprd(1,kq) = 0.0
wump(1,kq) =0.0
ncmp(1,kq) =0.0
nimp(1,kq) =0.0
sprd(1,kq) =0.0
frz(1,kq) =0.0
jk = kq
muu(1,JK) = 0.0
duu(1,JK) =0.0
EUU(1,JK) =0.0
cmel(1,JK) =0.0
cmei(1,JK) =0.0
oldTU(kq) = t0(kq)
oldQU(kq) = Q0(kq)
End do
Miter = 0
END IF

@weiwangncar
Copy link
Collaborator

@sael9740 The idea here is that the author may want to do further development of the scheme using qsatu without aerosol, and the use of qsatu may go outside the cloud layer (as you've found out). The variables aqnewlq and aqnewic are strictly in-cloud variables.

@sael9740
Copy link
Contributor Author

sael9740 commented Jun 6, 2023

@weiwangncar - That makes sense. Thanks for the explanation!

Should I create a pull request?

@weiwangncar
Copy link
Collaborator

@sael9740 Yes, please!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants