SUBROUTINE MSTCNV(IM,IX,KM,DT,T1,Q1,PRSL,DELPA,PRSLK,RAIN,CLW
&, rhc,lprnt,ipr)
!!
USE MACHINE , ONLY : kind_phys
USE FUNCPHYS , ONLY : fpvs, ftdp, fthe, stma, ftlcl
USE PHYSCONS, HVAP => con_HVAP, CP => con_CP, RV => con_RV
&, EPS => con_eps, EPSM1 => con_epsm1, grav => con_g
implicit none
!!
! PHYSICAL PARAMETERS
real(kind=kind_phys) elocp, el2orc
PARAMETER(ELOCP=HVAP/CP, EL2ORC=HVAP*HVAP/(RV*CP))
! real(kind=kind_phys), parameter :: pres_min=40000.0 ! in Pascals
real(kind=kind_phys), parameter :: pres_min=6000.0 ! in Pascals
!
!
logical lprnt
integer im,ix,km,ipr
real(kind=kind_phys) dt,rain(im)
!
real(kind=kind_phys) PRSL(IX,KM), PRSLK(IX,KM), DELPA(IX,KM)
real(kind=kind_phys) T1(IX,KM), Q1(IX,KM), CLW(IX,KM)
&, rhc(im,km)
!
! LOCAL VARIABLES
!
real(kind=kind_phys) P(IM,KM), TO(IM,KM), QO(IM,KM), QS(IM,KM)
&, THE(IM,KM), DQ(IM,KM), RAINLVL(IM,KM)
&, ES(IM,KM), DEL(IM,KM)
real(kind=kind_phys) pint(im), delqbar(im), deltbar(im), dqint(im)
&, ei(im), thebar(im), theint(im)
&, qevap, dpovg, rnevap, slklcl, tdpd
&, thelcl, tlcl, pmin
integer KMLEV(im,KM), k, kmax, kk(im), ks(im), ke(im), i
!
LOGICAL FLG(im), TOPFLG(im), TOTFLG
!
real(kind=kind_phys), parameter :: cons_0=0.d0,
& cons_1pdm8=1.d-8
!
KMAX = 0
DO K = 1, KM
pmin = 0.0
do i=1,im
pmin = min(prsl(i,k), pmin)
enddo
! IF(pmin .GT. 6000.0) KMAX = K + 1 ! input pressure is in Pa
IF(pmin .GE. pres_min) KMAX = K + 1 ! input pressure is in Pa
ENDDO
! if (lprnt) print *,' kmax=',kmax
!
! SURFACE PRESSURE UNIT IS Pa
!
do i=1,im
RAIN(i) = 0.
DELTBAR(i) = 0.
DELQBAR(i) = 0.
FLG(i) = .FALSE.
ks(i) = 0
! ke(i) = kmax + 1
TOPFLG(i) = .FALSE.
enddo
do k=1,kmax
do i=1,im
! if (lprnt .and. i == ipr) print *,' p=',p(i,k)
! &,' ke=',ke(i)
if (p(i,k) >= pres_min) ke(i) = k + 1
enddo
enddo
TOTFLG = .FALSE.
! if (lprnt) print *,' ke0=',ke(ipr),' ks0=',ks(ipr)
! &,' i=',ipr,'p=',p(ipr,kmax)
!
cselaDG3 IF(LAT.EQ.LATD) THEN
cselaDG3 PRINT *, ' T AND Q BEFORE ADJUSTMENT'
cselaDG3 PRINT 6000, (T1(k)-273.16,K=1,KMAX)
cselaDG3 PRINT 6000, (Q1(k)*1.E3,K=1,KMAX)
cselaDG3 PRINT *, ' PS =', PS
cselaDG3 ENDIF
!
!
! COLUMN VARIABLES
! P IS PRESSURE OF THE LAYER (Pa)
! TO IS TEMPERATURE AT T+DT (K)... THIS IS AFTER ADVECTION AND TURBULAN
! QO IS MIXING RATIO AT T+DT (KG/KG)..Q1
!
DO K = 1, KMAX
do i=1,im
P(i,k) = PRSL(I,K)
TO(i,k) = T1(i,k)
QO(i,k) = Q1(i,k)
DEL(i,k) = DELPA(i,k) * 0.001 ! Convert from Pa to kPa (aka Cb)
enddo
ENDDO
! if (lprnt) then
! print *,' toin=',to(ipr,:)
! print *,' pin=',p(ipr,:)
! print *,' prslk=',prslk(ipr,:)
! endif
!
! MODEL CONSISTENT SATURATION MIXING RATIO
!
DO K = 1, KMAX
do i=1,im
if (p(i,k) >= pres_min) then
ES(i,k) = min(P(i,k), FPVS(T1(i,k)))
QS(i,k) = EPS * ES(i,k) / (P(i,k) + EPSM1 * ES(i,k))
QS(i,k) = MAX(QS(i,k),cons_1pdm8)
IF(QO(i,k) .GT. QS(i,k)) FLG(i) = .TRUE.
endif
enddo
ENDDO
do i=1,im
IF (FLG(i)) TOTFLG = .TRUE.
enddo
IF (.NOT. TOTFLG) RETURN
!
DO K = 1, KMAX
do i=1,im
DQ(i,k) = 0.
THE(i,k) = TO(i,k)
enddo
ENDDO
! if (lprnt) print *,' qs=',(qs(ipr,k),k=1,kmax)
!
! COMPUTE THETA-E
!
DO K = 1, KMAX
do i=1,im
IF(FLG(i) .and. p(i,k) >= pres_min) THEN
THE(i,k) = FTHE(TO(i,k),prslk(i,k))
IF (THE(i,k) .EQ. 0.) THEN
THE(i,k) = TO(i,k) / prslk(i,k)
ENDIF
C THE(i,k) = TO(i,k) * ((P(i,k) - ES(i,k))*0.01) ** (-ROCP)
C & * EXP(ELOCP * QS(i,k) / TO(i,k))
DQ(i,k) = QO(i,k)- QS(i,k)
!
! MODIFICATION OF THETA-E FOR SUPER-SATURATION
!
THE(i,k)= THE(i,k) * (1. + HVAP*MAX(DQ(i,k),cons_0)
& / (CP*TO(i,k)))
! if (lprnt .and. i .eq. ipr) print *,' k=',k,' the=',the(i,k)
! &,' dq=',dq(i,k)
ENDIF
enddo
ENDDO
!
cselaDG3 IF(LAT.EQ.LATD.AND.FLG(LOND)) THEN
cselaDG3 PRINT *, ' THETA-E, QS AND DQ BEFORE ADJUSTMENT'
cselaDG3 PRINT 6000, (THE(k)-273.16,K=1,KMAX)
cselaDG3 PRINT 6000, (QS(k)*1.E3,K=1,KMAX)
cselaDG3 PRINT 6000, (DQ(k)*1.E3,K=1,KMAX)
cselaDG3 ENDIF
!
DO K = 1, KMAX
do i=1,im
KMLEV(i,k) = 0
RAINLVL(i,k) = 0.
enddo
ENDDO
!
! STARTING POINT OF ADJUSTMENT
!
k = 1
do i=1,im
KK(i) = 0
DQINT(i) = 0.
THEINT(i) = 0.
THEBAR(i) = 0.
PINT(i) = 0.
!
! FOR CONDITIONALLY UNSTABLE AND SUPERSATURATED LAYERS,
! OBTAIN INTEGRATED THETA AND Q-QS
!
! KMLEV KEEPS TRACK OF THE NUMBER OF LAYERS THAT SATISFIES
! THE CONDITION FOR ADJUSTMENT
!
IF(DQ(i,k).GT.0..AND.THE(i,k).GE.THE(i,K+1).AND.FLG(i)) THEN
DQINT(i) = DQINT(i) + DQ(i,k) * DEL(i,K)
THEINT(i) = THEINT(i) + THE(i,k) * DEL(i,K)
PINT(i) = PINT(i) + DEL(i,K)
KK(i) = KK(i) + 1
KMLEV(i,k) = KK(i)
ENDIF
enddo
! if (lprnt) print *,' kmlev=',kmlev(ipr,k),' k=',k
DO K = 2, KMAX - 1
do i=1,im
if(p(i,k) >= pres_min) then
IF(DQ(i,k).GT.0..AND.THE(i,k).GE.THE(i,K+1).AND.FLG(i)) THEN
DQINT(i) = DQINT(i) + DQ(i,k) * DEL(i,K)
THEINT(i) = THEINT(i) + THE(i,k) * DEL(i,K)
PINT(i) = PINT(i) + DEL(i,K)
KK(i) = KK(i) + 1
KMLEV(i,k) = KK(i)
ENDIF
! if (lprnt) print *,' kmlev=',kmlev(ipr,k),' k=',k
!
IF (PINT(i) .GT. 0.)THEBAR(i) = THEINT(i) / PINT(i)
!
! IF THE LAYER BELOW SATISFIES THE CONDITION AND THE PRESENT
! LAYER IS COLDER THAN THE ADJSUTED THETA-E,
! THE LAYER IS INCLUDED IF THE INTEGRATED MOISTURE EXCESS
! CAN BE MAINTAINED
!
IF (KMLEV(i,k) .EQ.0 .AND. KMLEV(i,K-1) .GT. 0 .AND.
& THEBAR(i) .GE. THE(i,k) .AND. .NOT. TOPFLG(i)) THEN
DQINT(i) = DQINT(i) + DQ(i,k) * DEL(i,K)
! ENDIF
! IF (KMLEV(i,k) .EQ. 0 .AND. KMLEV(i,K-1) .GT. 0 .AND.
! & THEBAR(i) .GE. THE(i,k) .AND. DQINT(i) .GT. 0.
! & .AND. .NOT. TOPFLG(i)) THEN
if (dqint(i) .gt. 0) then
KK(i) = KK(i) + 1
KMLEV(i,k) = KK(i)
TOPFLG(i) = .TRUE.
EI(i) = P(i,k) * QO(i,k) / (EPS - EPSM1 * QO(i,k))
EI(i) = MIN(MAX(EI(i),cons_1pdm8),ES(i,k)) !constant
TDPD = MAX(TO(i,k)-FTDP(EI(i)),cons_0) !constant
TLCL = FTLCL(TO(i,k), TDPD)
SLKLCL = prSLK(i,K) * TLCL / TO(i,k)
THELCL = FTHE(TLCL,SLKLCL)
IF(THELCL.NE.0.) THEN
THE(i,k) = THELCL
C THE(i,k) = TO(i,k) * ((P(i,k) - EI(i))*.01) ** (-ROCP)
C & * EXP(ELOCP * MAX(QO(i,k),1.E-8) / TO(i,k))
ENDIF
THEINT(i) = THEINT(i) + THE(i,k) * DEL(i,K)
PINT(i) = PINT(i) + DEL(i,K)
endif
ENDIF
endif
!
! RESET THE INTEGRAL IF THE LAYER IS NOT IN THE CLOUD
!
! if (lprnt) print *,' kmlev3=',kmlev(ipr,k),' k=',k
IF (KMLEV(i,k) .EQ. 0 .AND. KMLEV(i,K-1) .GT. 0) THEN
THEBAR(i) = THEINT(i) / PINT(i)
DQINT(i) = 0.
THEINT(i) = 0.
PINT(i) = 0.
KK(i) = 0
ks(i) = k - 1
ke(i) = ks(i) - kmlev(i,k-1) + 1
flg(i) = .false.
ENDIF
enddo
enddo
!
! When within A CLOUD LAYER, COMPUTE THE MOIST-ADIABATIC
! (TO AND QO) USING THE AVERAGED THETA-E AND THE RESULTANT RAIN
!
! if (lprnt) print *,' ke=',ke(ipr),' ks=',ks(ipr)
do k = 1, kmax
do i=1,im
if (k .ge. ke(i) .and. k .le. ks(i)) then
CALL STMA(THEBAR(i),PRSLK(i,k),TO(i,k),QO(i,k))
THE(i,k) = THEBAR(i)
QS(i,k) = QO(i,k)
!
DPOVG = DEL(i,K) / grav
RAINLVL(i,k) = (Q1(i,k) - QO(i,k)) * DPOVG
clw(i,k) = clw(i,k) + Q1(i,k) - QO(i,k)
DELTBAR(i) = DELTBAR(i) + (TO(i,k) - T1(i,k)) * DPOVG
& / PRSLK(i,k)
DELQBAR(i) = DELQBAR(i) - RAINLVL(i,K)
! if (lprnt) print *,' k=',k,' to=',to(i,k),' qo=',qo(i,k),
! & ' rainlvl=',rainlvl(i,k)
ENDIF
enddo
ENDDO
!
! EVAPORATION OF FALLING RAIN
!
! DO K = KMAX, 1, -1
! do i=1,im
! T1(i,k) = TO(i,k)
! Q1(i,k) = QO(i,k)
! DPOVG = DEL(i,K) / grav
! IF (RAIN(i) .GT. 0. .AND. RAINLVL(i,k) .LE. 0.) THEN
! DQ(i,k) = (QO(i,k) - QS(i,k)*rhc(i,k))
! & / (1. + EL2ORC*QS(i,k)/(TO(i,k)*TO(i,k)))
! QEVAP =-DQ(i,k)*(1.-EXP(-0.32*SQRT(DT*RAIN(i))))
!! & -clw(i,k)
! RNEVAP = MIN(QEVAP*DPOVG,RAIN(i))
! Q1(i,k) = Q1(i,k)+RNEVAP/DPOVG
! T1(i,k) = T1(i,k)-RNEVAP/DPOVG*ELOCP
! RAIN(i) = RAIN(i) - RNEVAP
! DELTBAR(i) = DELTBAR(i) - RNEVAP * ELOCP
! DELQBAR(i) = DELQBAR(i) + RNEVAP
! ELSE
! RAIN(i) = RAIN(i) + RAINLVL(i,k)
! ENDIF
! enddo
! ENDDO
!
cselaDG3 IF(LAT.EQ.LATD.AND.FLG(LOND)) THEN
cselaDG3 PRINT *, ' THETA-E AFTER ADJUSTMENT'
cselaDG3 PRINT 6000, (THE(k)-273.16,K=1,KMAX)
cselaDG3 PRINT *, ' T AND Q AFTER ADJUSTMENT'
cselaDG3 PRINT 6000, (T1(k)-273.16,K=1,KMAX)
cselaDG3 PRINT 6000, (Q1(k)*1.E3,K=1,KMAX)
cselaDG3 PRINT *, ' DELTBAR, DELQBAR =', DELTBAR*CP,DELQBAR*HVAP
cselaDG3 PRINT *, ' RAIN =', HVAP*RAIN
cselaDG3 ENDIF
!6000 FORMAT(2X,0P,11(F6.2,1H,))
!6100 FORMAT(2X,3P11F7.2)
!
! do i=1,im
! RAIN(i) = MAX(RAIN(i),cons_0)
! enddo
! if (lprnt) print *,' rain_in_mst=',rain(ipr)
!
RETURN
END