SUBROUTINE MSTCNV(IM,IX,KM,DT,T1,Q1,PRSL,DEL,PRSLK,RAIN &, lprnt,ipr) !! USE MACHINE , ONLY : kind_phys USE FUNCPHYS , ONLY : fpvs, ftdp, fthe, stma, fpkap, 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)) ! ! logical lprnt integer im,ix,km,ipr real(kind=kind_phys) dt,rain(im) ! real(kind=kind_phys) PRSL(IX,KM), PRSLK(IX,KM), DEL(IX,KM) real(kind=kind_phys) T1(IX,KM), Q1(IX,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) 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, pmax integer KMLEV(im,KM), k, kmax, kk(im), ks(im), ke(im), i ! LOGICAL FLG(im), TOPFLG(im), TOTFLG cc cc-------------------------------------------------------------------- cc real(kind=kind_phys) cons_0 !constant real(kind=kind_phys) cons_1pdm8 !constant cc cons_0 = 0.d0 !constant cons_1pdm8 = 1.d-8 !constant cc cc-------------------------------------------------------------------- cc KMAX = 0 DO K = 1, KM pmax = 0.0 do i=1,im pmax = max(prsl(i,k), pmax) enddo IF(pmax .GT. 60.0) KMAX = K + 1 ENDDO if (lprnt) print *,' kmax=',kmax ! cselaDG3 LATD = 57 cselaDG3 LOND = 1 cselaDG3 LATD = 48 cselaDG3 LOND = 176 C C SURFACE PRESSURE UNIT IS CB C 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 TOTFLG = .FALSE. ! 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 ! C C COLUMN VARIABLES C P IS PRESSURE OF THE LAYER (CB) C TO IS TEMPERATURE AT T+DT (K)... THIS IS AFTER ADVECTION AND TURBULAN C QO IS MIXING RATIO AT T+DT (KG/KG)..Q1 C DO K = 1, KMAX do i=1,im P(i,k) = 1000.0 * PRSL(I,K) TO(i,k) = T1(i,k) QO(i,k) = Q1(i,k) enddo ENDDO C C MODEL CONSISTENT SATURATION MIXING RATIO C DO K = 1, KMAX do i=1,im 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) !constant IF(QO(i,k) .GT. QS(i,k)) FLG(i) = .TRUE. 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 C C COMPUTE THETA-E C DO K = 1, KMAX do i=1,im IF(FLG(i)) 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) C C MODIFICATION OF THETA-E FOR SUPER-SATURATION C THE(i,k)= THE(i,k) * (1. + HVAP*MAX(DQ(i,k),cons_0) !constant & / (CP*TO(i,k))) if (lprnt .and. i .eq. ipr) print *,' k=',k,' the=',the(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 C C STARTING POINT OF ADJUSTMENT C k = 1 do i=1,im KK(i) = 0 DQINT(i) = 0. THEINT(i) = 0. THEBAR(i) = 0. PINT(i) = 0. C C FOR CONDITIONALLY UNSTABLE AND SUPERSATURATED LAYERS, C OBTAIN INTEGRATED THETA AND Q-QS C C KMLEV KEEPS TRACK OF THE NUMBER OF LAYERS THAT SATISFIES C THE CONDITION FOR ADJUSTMENT C 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 DO K = 2, KMAX - 1 do i=1,im 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 (PINT(i) .GT. 0.)THEBAR(i) = THEINT(i) / PINT(i) C C IF THE LAYER BELOW SATISFIES THE CONDITION AND THE PRESENT C LAYER IS COLDER THAN THE ADJSUTED THETA-E, C THE LAYER IS INCLUDED IF THE INTEGRATED MOISTURE EXCESS C CAN BE MAINTAINED C 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 C C RESET THE INTEGRAL IF THE LAYER IS NOT IN THE CLOUD C 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 C C When within A CLOUD LAYER, COMPUTE THE MOIST-ADIABATIC C (TO AND QO) USING THE AVERAGED THETA-E AND THE RESULTANT RAIN C 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 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) RAIN(i) = RAIN(i) + RAINLVL(i,k) DQ(i,k) = (QO(i,k) - QS(i,k)) / (1. + & EL2ORC*QS(i,k)/TO(i,k)**2) DPOVG = DEL(i,K) / grav IF (RAIN(i) .GT. 0. .AND. RAINLVL(i,k) .LE. 0.) THEN QEVAP =-DQ(i,k)*(1.-EXP(-0.32*SQRT(DT*RAIN(i)))) 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 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) !constant enddo ! RETURN END