File: C:\NOAA\NEMS_11731\src\atmos\phys\sascnv.f

1           SUBROUTINE SASCNV(IM,IX,KM,JCAP,DELT,DELP,PRSLP,PSP,PHIL,QL,
2     !     SUBROUTINE SASCNV(IM,IX,KM,JCAP,DELT,DEL,PRSL,PHIL,QL,
3          &       Q1,T1,U1,V1,CLDWRK,RN,KBOT,KTOP,KUO,SLIMSK,
4          &       DOT,XKT2,ncloud,ud_mf,dd_mf,dt_mf)
5     ! hchuang code change [r1L]
6     !    &       DOT,XKT2,ncloud)
7     !
8     ! 10/14/2008 Ho-Chun Huang The Cloudmass flux fields was added by Jongil
9     !
10     !  for cloud water version
11     !     parameter(ncloud=0)
12     !     SUBROUTINE SASCNV(KM,JCAP,DELT,DEL,SL,SLK,PS,QL,
13     !    &       Q1,T1,U1,V1,CLDWRK,RN,KBOT,KTOP,KUO,SLIMSK,
14     !    &       DOT,xkt2,ncloud)
15     !
16           USE MACHINE , ONLY : kind_phys
17           USE FUNCPHYS , ONLY : fpvs
18           USE PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP
19          &,             RV => con_RV, FV => con_fvirt, T0C => con_T0C
20          &,             CVAP => con_CVAP, CLIQ => con_CLIQ
21          &,             EPS => con_eps, EPSM1 => con_epsm1
22           implicit none
23     !
24     !
25           integer            IM, IX,  KM, JCAP, ncloud,
26          &                   KBOT(IM), KTOP(IM), KUO(IM)
27           real(kind=kind_phys) DELT
28           real(kind=kind_phys) PSP(IM),    DELP(IX,KM), PRSLP(IX,KM)
29           real(kind=kind_phys) PS(IM),     DEL(IX,KM),  PRSL(IX,KM),
30     !     real(kind=kind_phys)             DEL(IX,KM),  PRSL(IX,KM),
31          &                     QL(IX,KM,2),Q1(IX,KM),   T1(IX,KM),
32          &                     U1(IX,KM),  V1(IX,KM),
33          &                     CLDWRK(IM), RN(IM),      SLIMSK(IM),
34          &                     DOT(IX,KM), XKT2(IM),    PHIL(IX,KM)
35     ! hchuang code change [+1L] mass flux output
36          &,                    ud_mf(IM,KM),dd_mf(IM,KM),dt_mf(IM,KM)
37     !
38           integer              I, INDX, jmn, k, knumb, latd, lond, km1
39     !
40           real(kind=kind_phys) adw,     alpha,   alphal,  alphas,
41          &                     aup,     beta,    betal,   betas,
42          &                     c0,      cpoel,   dellat,  delta,
43          &                     desdt,   deta,    detad,   dg,
44          &                     dh,      dhh,     dlnsig,  dp,
45          &                     dq,      dqsdp,   dqsdt,   dt,
46          &                     dt2,     dtmax,   dtmin,   dv1,
47          &                     dv1q,    dv2,     dv2q,    dv1u,
48          &                     dv1v,    dv2u,    dv2v,    dv3u,
49          &                     dv3v,    dv3,     dv3q,    dvq1,
50          &                     dz,      dz1,     e1,      edtmax,
51          &                     edtmaxl, edtmaxs, el2orc,  elocp,
52          &                     es,      etah,
53          &                     evef,    evfact,  evfactl, fact1,
54          &                     fact2,   factor,  fjcap,   fkm,
55          &                     fuv,     g,       gamma,   onemf,
56          &                     onemfu,  pdetrn,  pdpdwn,  pprime,
57          &                     qc,      qlk,     qrch,    qs,
58          &                     rain,    rfact,   shear,   tem1,
59          &                     tem2,    terr,    val,     val1,
60          &                     val2,    w1,      w1l,     w1s,
61          &                     w2,      w2l,     w2s,     w3,
62          &                     w3l,     w3s,     w4,      w4l,
63          &                     w4s,     xdby,    xpw,     xpwd,
64          &                     xqc,     xqrch,   xlambu,  mbdt,
65          &                     tem
66     !
67     !
68           integer              JMIN(IM), KB(IM), KBCON(IM), KBDTR(IM),
69          &                     KT2(IM),  KTCON(IM), LMIN(IM),
70          &                     kbm(IM),  kbmax(IM), kmax(IM)
71     !
72           real(kind=kind_phys) AA1(IM),     ACRT(IM),   ACRTFCT(IM),
73          &                     DELHBAR(IM), DELQ(IM),   DELQ2(IM),
74          &                     DELQBAR(IM), DELQEV(IM), DELTBAR(IM),
75          &                     DELTV(IM),   DTCONV(IM), EDT(IM),
76          &                     EDTO(IM),    EDTX(IM),   FLD(IM),
77          &                     HCDO(IM),    HKBO(IM),   HMAX(IM),
78          &                     HMIN(IM),    HSBAR(IM),  UCDO(IM),
79          &                     UKBO(IM),    VCDO(IM),   VKBO(IM),
80          &                     PBCDIF(IM),  PDOT(IM),   PO(IM,KM),
81          &                                  PWAVO(IM),  PWEVO(IM),
82     !    &                     PSFC(IM),    PWAVO(IM),  PWEVO(IM),
83          &                     QCDO(IM),    QCOND(IM),  QEVAP(IM),
84          &                     QKBO(IM),    RNTOT(IM),  VSHEAR(IM),
85          &                     XAA0(IM),    XHCD(IM),   XHKB(IM),
86          &                     XK(IM),      XLAMB(IM),  XLAMD(IM),
87          &                     XMB(IM),     XMBMAX(IM), XPWAV(IM),
88          &                     XPWEV(IM),   XQCD(IM),   XQKB(IM)
89     cc
90     C  PHYSICAL PARAMETERS
91           PARAMETER(G=grav)
92           PARAMETER(CPOEL=CP/HVAP,ELOCP=HVAP/CP,
93          &          EL2ORC=HVAP*HVAP/(RV*CP))
94           PARAMETER(TERR=0.,C0=.002,DELTA=fv)
95           PARAMETER(FACT1=(CVAP-CLIQ)/RV,FACT2=HVAP/RV-FACT1*T0C)
96     C  LOCAL VARIABLES AND ARRAYS
97           real(kind=kind_phys) PFLD(IM,KM),    TO(IM,KM),     QO(IM,KM),
98          &                     UO(IM,KM),      VO(IM,KM),     QESO(IM,KM)
99     c  cloud water
100           real(kind=kind_phys) QLKO_KTCON(IM), DELLAL(IM),    TVO(IM,KM),
101          &                     DBYO(IM,KM),    ZO(IM,KM),     SUMZ(IM,KM),
102          &                     SUMH(IM,KM),    HEO(IM,KM),    HESO(IM,KM),
103          &                     QRCD(IM,KM),    DELLAH(IM,KM), DELLAQ(IM,KM),
104          &                     DELLAU(IM,KM),  DELLAV(IM,KM), HCKO(IM,KM),
105          &                     UCKO(IM,KM),    VCKO(IM,KM),   QCKO(IM,KM),
106          &                     ETA(IM,KM),     ETAU(IM,KM),   ETAD(IM,KM),
107          &                     QRCDO(IM,KM),   PWO(IM,KM),    PWDO(IM,KM),
108          &                     RHBAR(IM),      TX1(IM)
109     !
110           LOGICAL TOTFLG, CNVFLG(IM), DWNFLG(IM), DWNFLG2(IM), FLG(IM)
111     !
112           real(kind=kind_phys) PCRIT(15), ACRITT(15), ACRIT(15)
113     cmy      SAVE PCRIT, ACRITT
114           DATA PCRIT/850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,
115          &           350.,300.,250.,200.,150./
116           DATA ACRITT/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216,
117          &           .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
118     C  GDAS DERIVED ACRIT
119     C     DATA ACRITT/.203,.515,.521,.566,.625,.665,.659,.688,
120     C    &            .743,.813,.886,.947,1.138,1.377,1.896/
121     cc
122           real(kind=kind_phys) tf, tcr, tcrf
123           parameter (TF=233.16, TCR=263.16, TCRF=1.0/(TCR-TF)) ! From Lord(1978)
124     !
125     !     parameter (tf=258.16, tcr=273.16, tcrf=1.0/(tcr-tf))
126     !
127           real(kind=kind_phys), parameter :: cons_0=0.0
128     c
129     c--------------------------------------------------------------------
130     !
131     !************************************************************************
132     !     convert input Pa terms to Cb terms  -- Moorthi
133           ps   = psp   * 0.001
134           prsl = prslp * 0.001
135           del  = delp  * 0.001
136     !************************************************************************
137     !
138           km1 = km - 1
139     C  INITIALIZE ARRAYS
140     C
141           DO I=1,IM
142             RN(I)=0.
143             KBOT(I)=KM+1
144             KTOP(I)=0
145     !       KUO(I)=0
146             CNVFLG(I) = .TRUE.
147             DTCONV(I) = 3600.
148             CLDWRK(I) = 0.
149             PDOT(I) = 0.
150             KT2(I) = 0
151             QLKO_KTCON(I) = 0.
152             DELLAL(I) = 0.
153           ENDDO
154     ! hchuang code change [+7L]
155           DO K = 1, KM
156             DO I=1,IM
157               ud_mf(I,k) = 0.
158               dd_mf(I,k) = 0.
159               dt_mf(I,k) = 0.
160             ENDDO
161           ENDDO
162     !!
163           DO K = 1, 15
164             ACRIT(K) = ACRITT(K) * (975. - PCRIT(K))
165           ENDDO
166           DT2 = DELT
167           val   =         1200.
168           dtmin = max(dt2, val )
169           val   =         3600.
170           dtmax = max(dt2, val )
171     C  MODEL TUNABLE PARAMETERS ARE ALL HERE
172           MBDT    = 10.
173           EDTMAXl = .3
174           EDTMAXs = .3
175           ALPHAl  = .5
176           ALPHAs  = .5
177           BETAl   = .15
178           betas   = .15
179           BETAl   = .05
180           betas   = .05
181     c     EVEF    = 0.07
182           evfact  = 0.3
183           evfactl = 0.3
184           PDPDWN  = 0.
185           PDETRN  = 200.
186           xlambu  = 1.e-4
187           fjcap   = (float(jcap) / 126.) ** 2
188           val     =           1.
189           fjcap   = max(fjcap,val)
190           fkm     = (float(km) / 28.) ** 2
191           fkm     = max(fkm,val)
192           W1l     = -8.E-3
193           W2l     = -4.E-2
194           W3l     = -5.E-3
195           W4l     = -5.E-4
196           W1s     = -2.E-4
197           W2s     = -2.E-3
198           W3s     = -1.E-3
199           W4s     = -2.E-5
200     CCCCC IF(IM.EQ.384) THEN
201             LATD  = 92
202             lond  = 189
203     CCCCC ELSEIF(IM.EQ.768) THEN
204     CCCCC   LATD = 80
205     CCCCC ELSE
206     CCCCC   LATD = 0
207     CCCCC ENDIF
208     C
209     C  DEFINE TOP LAYER FOR SEARCH OF THE DOWNDRAFT ORIGINATING LAYER
210     C  AND THE MAXIMUM THETAE FOR UPDRAFT
211     C
212           DO I=1,IM
213             KBMAX(I) = KM
214             KBM(I)   = KM
215             KMAX(I)  = KM
216             TX1(I)   = 1.0 / PS(I)
217           ENDDO
218     !
219           DO K = 1, KM
220             DO I=1,IM
221               IF (prSL(I,K)*tx1(I) .GT. 0.45) KBMAX(I) = K + 1
222               IF (prSL(I,K)*tx1(I) .GT. 0.70) KBM(I)   = K + 1
223               IF (prSL(I,K)*tx1(I) .GT. 0.04) KMAX(I)  = K + 1
224             ENDDO
225           ENDDO
226           DO I=1,IM
227             KBMAX(I) = MIN(KBMAX(I),KMAX(I))
228             KBM(I)   = MIN(KBM(I),KMAX(I))
229           ENDDO
230     C
231     C   CONVERT SURFACE PRESSURE TO MB FROM CB
232     C
233     !!
234           DO K = 1, KM
235             DO I=1,IM
236               if (K .le. kmax(i)) then
237                 PFLD(I,k) = PRSL(I,K) * 10.0
238                 PWO(I,k)  = 0.
239                 PWDO(I,k) = 0.
240                 TO(I,k)   = T1(I,k)
241                 QO(I,k)   = Q1(I,k)
242                 UO(I,k)   = U1(I,k)
243                 VO(I,k)   = V1(I,k)
244                 DBYO(I,k) = 0.
245                 SUMZ(I,k) = 0.
246                 SUMH(I,k) = 0.
247               endif
248             ENDDO
249           ENDDO
250     C
251     C  COLUMN VARIABLES
252     C  P IS PRESSURE OF THE LAYER (MB)
253     C  T IS TEMPERATURE AT T-DT (K)..TN
254     C  Q IS MIXING RATIO AT T-DT (KG/KG)..QN
255     C  TO IS TEMPERATURE AT T+DT (K)... THIS IS AFTER ADVECTION AND TURBULAN
256     C  QO IS MIXING RATIO AT T+DT (KG/KG)..Q1
257     C
258           DO K = 1, KM
259             DO I=1,IM
260               if (k .le. kmax(i)) then
261     !jfe        QESO(I,k) = 10. * FPVS(T1(I,k))
262     !
263                 QESO(I,k) = 0.01 * fpvs(T1(I,K))      ! fpvs is in Pa
264     !
265                 QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k) + EPSM1*QESO(I,k))
266                 val1      =             1.E-8
267                 QESO(I,k) = MAX(QESO(I,k), val1)
268                 val2      =           1.e-10
269                 QO(I,k)   = max(QO(I,k), val2 )
270     c           QO(I,k)   = MIN(QO(I,k),QESO(I,k))
271                 TVO(I,k)  = TO(I,k) + DELTA * TO(I,k) * QO(I,k)
272               endif
273             ENDDO
274           ENDDO
275     C
276     C  HYDROSTATIC HEIGHT ASSUME ZERO TERR
277     C
278           DO K = 1, KM
279             DO I=1,IM
280               ZO(I,k) = PHIL(I,k) / G
281             ENDDO
282           ENDDO
283     C  COMPUTE MOIST STATIC ENERGY
284           DO K = 1, KM
285             DO I=1,IM
286               if (K .le. kmax(i)) then
287     !           tem       = G * ZO(I,k) + CP * TO(I,k)
288                 tem       = PHIL(I,k) + CP * TO(I,k)
289                 HEO(I,k)  = tem  + HVAP * QO(I,k)
290                 HESO(I,k) = tem  + HVAP * QESO(I,k)
291     C           HEO(I,k)  = MIN(HEO(I,k),HESO(I,k))
292               endif
293             ENDDO
294           ENDDO
295     C
296     C  DETERMINE LEVEL WITH LARGEST MOIST STATIC ENERGY
297     C  THIS IS THE LEVEL WHERE UPDRAFT STARTS
298     C
299           DO I=1,IM
300             HMAX(I) = HEO(I,1)
301             KB(I) = 1
302           ENDDO
303     !!
304           DO K = 2, KM
305             DO I=1,IM
306               if (k .le. kbm(i)) then
307                 IF(HEO(I,k).GT.HMAX(I).AND.CNVFLG(I)) THEN
308                   KB(I)   = K
309                   HMAX(I) = HEO(I,k)
310                 ENDIF
311               endif
312             ENDDO
313           ENDDO
314     C     DO K = 1, KMAX - 1
315     C         TOL(k) = .5 * (TO(I,k) + TO(I,k+1))
316     C         QOL(k) = .5 * (QO(I,k) + QO(I,k+1))
317     C         QESOL(I,k) = .5 * (QESO(I,k) + QESO(I,k+1))
318     C         HEOL(I,k) = .5 * (HEO(I,k) + HEO(I,k+1))
319     C         HESOL(I,k) = .5 * (HESO(I,k) + HESO(I,k+1))
320     C     ENDDO
321           DO K = 1, KM1
322             DO I=1,IM
323               if (k .le. kmax(i)-1) then
324                 DZ      = .5 * (ZO(I,k+1) - ZO(I,k))
325                 DP      = .5 * (PFLD(I,k+1) - PFLD(I,k))
326     !jfe        ES      = 10. * FPVS(TO(I,k+1))
327     !
328                 ES      = 0.01 * fpvs(TO(I,K+1))      ! fpvs is in Pa
329     !
330                 PPRIME  = PFLD(I,k+1) + EPSM1 * ES
331                 QS      = EPS * ES / PPRIME
332                 DQSDP   = - QS / PPRIME
333                 DESDT   = ES * (FACT1 / TO(I,k+1) + FACT2 / (TO(I,k+1)**2))
334                 DQSDT   = QS * PFLD(I,k+1) * DESDT / (ES * PPRIME)
335                 GAMMA   = EL2ORC * QESO(I,k+1) / (TO(I,k+1)**2)
336                 DT      = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
337                 DQ      = DQSDT * DT + DQSDP * DP
338                 TO(I,k) = TO(I,k+1) + DT
339                 QO(I,k) = QO(I,k+1) + DQ
340                 PO(I,k) = .5 * (PFLD(I,k) + PFLD(I,k+1))
341               endif
342             ENDDO
343           ENDDO
344     !
345           DO K = 1, KM1
346             DO I=1,IM
347               if (k .le. kmax(I)-1) then
348     !jfe        QESO(I,k) = 10. * FPVS(TO(I,k))
349     !
350                 QESO(I,k) = 0.01 * fpvs(TO(I,K))      ! fpvs is in Pa
351     !
352                 QESO(I,k) = EPS * QESO(I,k) / (PO(I,k) + EPSM1*QESO(I,k))
353                 val1      =             1.E-8
354                 QESO(I,k) = MAX(QESO(I,k), val1)
355                 val2      =           1.e-10
356                 QO(I,k)   = max(QO(I,k), val2 )
357     c           QO(I,k)   = MIN(QO(I,k),QESO(I,k))
358                 HEO(I,k)  = .5 * G * (ZO(I,k) + ZO(I,k+1)) +
359          &                  CP * TO(I,k) + HVAP * QO(I,k)
360                 HESO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) +
361          &                  CP * TO(I,k) + HVAP * QESO(I,k)
362                 UO(I,k)   = .5 * (UO(I,k) + UO(I,k+1))
363                 VO(I,k)   = .5 * (VO(I,k) + VO(I,k+1))
364               endif
365             ENDDO
366           ENDDO
367     c     k = kmax
368     c       HEO(I,k) = HEO(I,k)
369     c       hesol(k) = HESO(I,k)
370     c      IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
371     c        PRINT *, '   HEO ='
372     c        PRINT 6001, (HEO(I,K),K=1,KMAX)
373     c        PRINT *, '   HESO ='
374     c        PRINT 6001, (HESO(I,K),K=1,KMAX)
375     c        PRINT *, '   TO ='
376     c        PRINT 6002, (TO(I,K)-273.16,K=1,KMAX)
377     c        PRINT *, '   QO ='
378     c        PRINT 6003, (QO(I,K),K=1,KMAX)
379     c        PRINT *, '   QSO ='
380     c        PRINT 6003, (QESO(I,K),K=1,KMAX)
381     c      ENDIF
382     C
383     C  LOOK FOR CONVECTIVE CLOUD BASE AS THE LEVEL OF FREE CONVECTION
384     C
385           DO I=1,IM
386             IF(CNVFLG(I)) THEN
387               INDX    = KB(I)
388               HKBO(I) = HEO(I,INDX)
389               QKBO(I) = QO(I,INDX)
390               UKBO(I) = UO(I,INDX)
391               VKBO(I) = VO(I,INDX)
392             ENDIF
393             FLG(I)    = CNVFLG(I)
394             KBCON(I)  = KMAX(I)
395           ENDDO
396     !!
397           DO K = 1, KM
398             DO I=1,IM
399               if (k .le. kbmax(i)) then
400                 IF(FLG(I).AND.K.GT.KB(I)) THEN
401                   HSBAR(I)   = HESO(I,k)
402                   IF(HKBO(I).GT.HSBAR(I)) THEN
403                     FLG(I)   = .FALSE.
404                     KBCON(I) = K
405                   ENDIF
406                 ENDIF
407               endif
408             ENDDO
409           ENDDO
410           DO I=1,IM
411             IF(CNVFLG(I)) THEN
412               PBCDIF(I) = -PFLD(I,KBCON(I)) + PFLD(I,KB(I))
413     !         PDOT(I)   = 10.* DOT(I,KBCON(I))
414               PDOT(I)   = 0.01 * DOT(I,KBCON(I))  ! Now DOT is in Pa/s
415               IF(PBCDIF(I).GT.150.)    CNVFLG(I) = .FALSE.
416               IF(KBCON(I).EQ.KMAX(I))  CNVFLG(I) = .FALSE.
417             ENDIF
418           ENDDO
419     !!
420           TOTFLG = .TRUE.
421           DO I=1,IM
422             TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
423           ENDDO
424           IF(TOTFLG) RETURN
425     C  FOUND LFC, CAN DEFINE REST OF VARIABLES
426      6001 FORMAT(2X,-2P10F12.2)
427      6002 FORMAT(2X,10F12.2)
428      6003 FORMAT(2X,3P10F12.2)
429     C
430     C  DETERMINE ENTRAINMENT RATE BETWEEN KB AND KBCON
431     C
432           DO I = 1, IM
433             alpha = alphas
434             if(SLIMSK(I).eq.1.) alpha = alphal
435             IF(CNVFLG(I)) THEN
436               IF(KB(I).EQ.1) THEN
437                 DZ = .5 * (ZO(I,KBCON(I)) + ZO(I,KBCON(I)-1)) - ZO(I,1)
438               ELSE
439                 DZ = .5 * (ZO(I,KBCON(I)) + ZO(I,KBCON(I)-1))
440          &         - .5 * (ZO(I,KB(I)) + ZO(I,KB(I)-1))
441               ENDIF
442               IF(KBCON(I).NE.KB(I)) THEN
443                 XLAMB(I) = - LOG(ALPHA) / DZ
444               ELSE
445                 XLAMB(I) = 0.
446               ENDIF
447             ENDIF
448           ENDDO
449     C  DETERMINE UPDRAFT MASS FLUX
450           DO K = 1, KM
451             DO I = 1, IM
452               if (k .le. kmax(i) .and. CNVFLG(I)) then
453                 ETA(I,k)  = 1.
454                 ETAU(I,k) = 1.
455               ENDIF
456             ENDDO
457           ENDDO
458           DO K = KM1, 2, -1
459             DO I = 1, IM
460               if (k .le. kbmax(i)) then
461                 IF(CNVFLG(I).AND.K.LT.KBCON(I).AND.K.GE.KB(I)) THEN
462                   DZ        = .5 * (ZO(I,k+1) - ZO(I,k-1))
463                   ETA(I,k)  = ETA(I,k+1) * EXP(-XLAMB(I) * DZ)
464                   ETAU(I,k) = ETA(I,k)
465                 ENDIF
466               endif
467             ENDDO
468           ENDDO
469           DO I = 1, IM
470             IF(CNVFLG(I).AND.KB(I).EQ.1.AND.KBCON(I).GT.1) THEN
471               DZ = .5 * (ZO(I,2) - ZO(I,1))
472               ETA(I,1) = ETA(I,2) * EXP(-XLAMB(I) * DZ)
473               ETAU(I,1) = ETA(I,1)
474             ENDIF
475           ENDDO
476     C
477     C  WORK UP UPDRAFT CLOUD PROPERTIES
478     C
479           DO I = 1, IM
480             IF(CNVFLG(I)) THEN
481               INDX         = KB(I)
482               HCKO(I,INDX) = HKBO(I)
483               QCKO(I,INDX) = QKBO(I)
484               UCKO(I,INDX) = UKBO(I)
485               VCKO(I,INDX) = VKBO(I)
486               PWAVO(I)     = 0.
487             ENDIF
488           ENDDO
489     C
490     C  CLOUD PROPERTY BELOW CLOUD BASE IS MODIFIED BY THE ENTRAINMENT PROCES
491     C
492           DO K = 2, KM1
493             DO I = 1, IM
494               if (k .le. kmax(i)-1) then
495                 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KBCON(I)) THEN
496                   FACTOR = ETA(I,k-1) / ETA(I,k)
497                   ONEMF = 1. - FACTOR
498                   HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF *
499          &                    .5 * (HEO(I,k) + HEO(I,k+1))
500                   UCKO(I,k) = FACTOR * UCKO(I,k-1) + ONEMF *
501          &                    .5 * (UO(I,k) + UO(I,k+1))
502                   VCKO(I,k) = FACTOR * VCKO(I,k-1) + ONEMF *
503          &                    .5 * (VO(I,k) + VO(I,k+1))
504                   DBYO(I,k) = HCKO(I,k) - HESO(I,k)
505                 ENDIF
506                 IF(CNVFLG(I).AND.K.GT.KBCON(I)) THEN
507                   HCKO(I,k) = HCKO(I,k-1)
508                   UCKO(I,k) = UCKO(I,k-1)
509                   VCKO(I,k) = VCKO(I,k-1)
510                   DBYO(I,k) = HCKO(I,k) - HESO(I,k)
511                 ENDIF
512               endif
513             ENDDO
514           ENDDO
515     C  DETERMINE CLOUD TOP
516           DO I = 1, IM
517             FLG(I) = CNVFLG(I)
518             KTCON(I) = 1
519           ENDDO
520     C     DO K = 2, KMAX
521     C       KK = KMAX - K + 1
522     C         IF(DBYO(I,kK).GE.0..AND.FLG(I).AND.KK.GT.KBCON(I)) THEN
523     C           KTCON(I) = KK + 1
524     C           FLG(I) = .FALSE.
525     C         ENDIF
526     C     ENDDO
527           DO K = 2, KM
528             DO I = 1, IM
529               if (k .le. kmax(i)) then
530                 IF(DBYO(I,k).LT.0..AND.FLG(I).AND.K.GT.KBCON(I)) THEN
531                   KTCON(I) = K
532                   FLG(I) = .FALSE.
533                 ENDIF
534               endif
535             ENDDO
536           ENDDO
537           DO I = 1, IM
538             IF(CNVFLG(I).AND.(PFLD(I,KBCON(I)) - PFLD(I,KTCON(I))).LT.150.)
539          &  CNVFLG(I) = .FALSE.
540           ENDDO
541           TOTFLG = .TRUE.
542           DO I = 1, IM
543             TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
544           ENDDO
545           IF(TOTFLG) RETURN
546     C
547     C  SEARCH FOR DOWNDRAFT ORIGINATING LEVEL ABOVE THETA-E MINIMUM
548     C
549           DO I = 1, IM
550             HMIN(I) = HEO(I,KBCON(I))
551             LMIN(I) = KBMAX(I)
552             JMIN(I) = KBMAX(I)
553           ENDDO
554           DO I = 1, IM
555             DO K = KBCON(I), KBMAX(I)
556               IF(HEO(I,k).LT.HMIN(I).AND.CNVFLG(I)) THEN
557                 LMIN(I) = K + 1
558                 HMIN(I) = HEO(I,k)
559               ENDIF
560             ENDDO
561           ENDDO
562     C
563     C  Make sure that JMIN(I) is within the cloud
564     C
565           DO I = 1, IM
566             IF(CNVFLG(I)) THEN
567               JMIN(I) = MIN(LMIN(I),KTCON(I)-1)
568               XMBMAX(I) = .1
569               JMIN(I) = MAX(JMIN(I),KBCON(I)+1)
570             ENDIF
571           ENDDO
572     C
573     C  ENTRAINING CLOUD
574     C
575           do k = 2, km1
576             DO I = 1, IM
577               if (k .le. kmax(i)-1) then
578                 if(CNVFLG(I).and.k.gt.JMIN(I).and.k.le.KTCON(I)) THEN
579                   SUMZ(I,k) = SUMZ(I,k-1) + .5 * (ZO(I,k+1) - ZO(I,k-1))
580                   SUMH(I,k) = SUMH(I,k-1) + .5 * (ZO(I,k+1) - ZO(I,k-1))
581          &                  * HEO(I,k)
582                 ENDIF
583               endif
584             enddo
585           enddo
586     !!
587           DO I = 1, IM
588             IF(CNVFLG(I)) THEN
589     c         call random_number(XKT2)
590     c         call srand(fhour)
591     c         XKT2(I) = rand()
592               KT2(I) = nint(XKT2(I)*float(KTCON(I)-JMIN(I))-.5)+JMIN(I)+1
593     !         KT2(I) = nint(sqrt(XKT2(I))*float(KTCON(I)-JMIN(I))-.5) + JMIN(I) + 1
594     c         KT2(I) = nint(ranf() *float(KTCON(I)-JMIN(I))-.5) + JMIN(I) + 1
595               tem1 = (HCKO(I,JMIN(I)) - HESO(I,KT2(I)))
596               tem2 = (SUMZ(I,KT2(I)) * HESO(I,KT2(I)) - SUMH(I,KT2(I)))
597               if (abs(tem2) .gt. 0.000001) THEN
598                 XLAMB(I) = tem1 / tem2
599               else
600                 CNVFLG(I) = .false.
601               ENDIF
602     !         XLAMB(I) = (HCKO(I,JMIN(I)) - HESO(I,KT2(I)))
603     !    &          / (SUMZ(I,KT2(I)) * HESO(I,KT2(I)) - SUMH(I,KT2(I)))
604               XLAMB(I) = max(XLAMB(I),cons_0)
605               XLAMB(I) = min(XLAMB(I),2.3/SUMZ(I,KT2(I)))
606             ENDIF
607           ENDDO
608     !!
609           DO I = 1, IM
610            DWNFLG(I)  = CNVFLG(I)
611            DWNFLG2(I) = CNVFLG(I)
612            IF(CNVFLG(I)) THEN
613             if(KT2(I).ge.KTCON(I)) DWNFLG(I) = .false.
614           if(XLAMB(I).le.1.e-30.or.HCKO(I,JMIN(I))-HESO(I,KT2(I)).le.1.e-30)
615          &  DWNFLG(I) = .false.
616             do k = JMIN(I), KT2(I)
617               if(DWNFLG(I).and.HEO(I,k).gt.HESO(I,KT2(I))) DWNFLG(I)=.false.
618             enddo
619     c       IF(CNVFLG(I).AND.(PFLD(KBCON(I))-PFLD(KTCON(I))).GT.PDETRN)
620     c    &     DWNFLG(I)=.FALSE.
621             IF(CNVFLG(I).AND.(PFLD(I,KBCON(I))-PFLD(I,KTCON(I))).LT.PDPDWN)
622          &     DWNFLG2(I)=.FALSE.
623            ENDIF
624           ENDDO
625     !!
626           DO K = 2, KM1
627             DO I = 1, IM
628               if (k .le. kmax(i)-1) then
629                 IF(DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KT2(I)) THEN
630                   DZ        = .5 * (ZO(I,k+1) - ZO(I,k-1))
631     c             ETA(I,k)  = ETA(I,k-1) * EXP( XLAMB(I) * DZ)
632     c  to simplify matter, we will take the linear approach here
633     c
634                   ETA(I,k)  = ETA(I,k-1) * (1. + XLAMB(I) * dz)
635                   ETAU(I,k) = ETAU(I,k-1) * (1. + (XLAMB(I)+xlambu) * dz)
636                 ENDIF
637               endif
638             ENDDO
639           ENDDO
640     !!
641           DO K = 2, KM1
642             DO I = 1, IM
643               if (k .le. kmax(i)-1) then
644     c           IF(.NOT.DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KT2(I)) THEN
645                 IF(.NOT.DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KTCON(I)) THEN
646                   DZ        = .5 * (ZO(I,k+1) - ZO(I,k-1))
647                   ETAU(I,k) = ETAU(I,k-1) * (1. + xlambu * dz)
648                 ENDIF
649               endif
650             ENDDO
651           ENDDO
652     c      IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
653     c        PRINT *, ' LMIN(I), KT2(I)=', LMIN(I), KT2(I)
654     c        PRINT *, ' KBOT, KTOP, JMIN(I) =', KBCON(I), KTCON(I), JMIN(I)
655     c      ENDIF
656     c     IF(LAT.EQ.LATD.AND.lon.eq.lond) THEN
657     c       print *, ' xlamb =', xlamb
658     c       print *, ' eta =', (eta(k),k=1,KT2(I))
659     c       print *, ' ETAU =', (ETAU(I,k),k=1,KT2(I))
660     c       print *, ' HCKO =', (HCKO(I,k),k=1,KT2(I))
661     c       print *, ' SUMZ =', (SUMZ(I,k),k=1,KT2(I))
662     c       print *, ' SUMH =', (SUMH(I,k),k=1,KT2(I))
663     c     ENDIF
664           DO I = 1, IM
665             if(DWNFLG(I)) THEN
666               KTCON(I) = KT2(I)
667             ENDIF
668           ENDDO
669     C
670     C  CLOUD PROPERTY ABOVE CLOUD Base IS MODIFIED BY THE DETRAINMENT PROCESS
671     C
672           DO K = 2, KM1
673             DO I = 1, IM
674               if (k .le. kmax(i)-1) then
675     cjfe
676                 IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
677     cjfe      IF(K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
678                   FACTOR    = ETA(I,k-1) / ETA(I,k)
679                   ONEMF     = 1. - FACTOR
680                   fuv       = ETAU(I,k-1) / ETAU(I,k)
681                   onemfu    = 1. - fuv
682                   HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF *
683          &                    .5 * (HEO(I,k) + HEO(I,k+1))
684                   UCKO(I,k) = fuv * UCKO(I,k-1) + ONEMFu *
685          &                    .5 * (UO(I,k) + UO(I,k+1))
686                   VCKO(I,k) = fuv * VCKO(I,k-1) + ONEMFu *
687          &                    .5 * (VO(I,k) + VO(I,k+1))
688                   DBYO(I,k) = HCKO(I,k) - HESO(I,k)
689                 ENDIF
690               endif
691             ENDDO
692           ENDDO
693     c      IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
694     c        PRINT *, ' UCKO=', (UCKO(I,k),k=KBCON(I)+1,KTCON(I))
695     c        PRINT *, ' uenv=', (.5*(UO(I,k)+UO(I,k-1)),k=KBCON(I)+1,KTCON(I))
696     c      ENDIF
697           DO I = 1, IM
698             if(CNVFLG(I).and.DWNFLG2(I).and.JMIN(I).le.KBCON(I))
699          &     THEN
700               CNVFLG(I) = .false.
701               DWNFLG(I) = .false.
702               DWNFLG2(I) = .false.
703             ENDIF
704           ENDDO
705     !!
706           TOTFLG = .TRUE.
707           DO I = 1, IM
708             TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
709           ENDDO
710           IF(TOTFLG) RETURN
711     !!
712     C
713     C  COMPUTE CLOUD MOISTURE PROPERTY AND PRECIPITATION
714     C
715           DO I = 1, IM
716               AA1(I) = 0.
717               RHBAR(I) = 0.
718           ENDDO
719           DO K = 1, KM
720             DO I = 1, IM
721               if (k .le. kmax(i)) then
722                 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LT.KTCON(I)) THEN
723                   DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
724                   DZ1 = (ZO(I,k) - ZO(I,k-1))
725                   GAMMA = EL2ORC * QESO(I,k) / (TO(I,k)**2)
726                   QRCH = QESO(I,k)
727          &             + GAMMA * DBYO(I,k) / (HVAP * (1. + GAMMA))
728                   FACTOR = ETA(I,k-1) / ETA(I,k)
729                   ONEMF = 1. - FACTOR
730                   QCKO(I,k) = FACTOR * QCKO(I,k-1) + ONEMF *
731          &                    .5 * (QO(I,k) + QO(I,k+1))
732                   DQ = ETA(I,k) * QCKO(I,k) - ETA(I,k) * QRCH
733                   RHBAR(I) = RHBAR(I) + QO(I,k) / QESO(I,k)
734     C
735     C  BELOW LFC CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
736     C
737                   IF(DQ.GT.0.) THEN
738                     ETAH = .5 * (ETA(I,k) + ETA(I,k-1))
739                     QLK = DQ / (ETA(I,k) + ETAH * C0 * DZ)
740                     AA1(I) = AA1(I) - DZ1 * G * QLK
741                     QC = QLK + QRCH
742                     PWO(I,k) = ETAH * C0 * DZ * QLK
743                     QCKO(I,k) = QC
744                     PWAVO(I) = PWAVO(I) + PWO(I,k)
745                   ENDIF
746                 ENDIF
747               endif
748             ENDDO
749           ENDDO
750           DO I = 1, IM
751             RHBAR(I) = RHBAR(I) / float(KTCON(I) - KB(I) - 1)
752           ENDDO
753     c
754     c  this section is ready for cloud water
755     c
756           if(ncloud.gt.0) THEN
757     c
758     c  compute liquid and vapor separation at cloud top
759     c
760           DO I = 1, IM
761             k = KTCON(I)
762             IF(CNVFLG(I)) THEN
763               GAMMA = EL2ORC * QESO(I,K) / (TO(I,K)**2)
764               QRCH = QESO(I,K)
765          &         + GAMMA * DBYO(I,K) / (HVAP * (1. + GAMMA))
766               DQ = QCKO(I,K-1) - QRCH
767     C
768     C  CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
769     C
770               IF(DQ.GT.0.) THEN
771                 QLKO_KTCON(I) = dq
772                 QCKO(I,K-1) = QRCH
773               ENDIF
774             ENDIF
775           ENDDO
776           ENDIF
777     C
778     C  CALCULATE CLOUD WORK FUNCTION AT T+DT
779     C
780           DO K = 1, KM
781             DO I = 1, IM
782               if (k .le. kmax(i)) then
783                 IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
784                   DZ1 = ZO(I,k) - ZO(I,k-1)
785                   GAMMA = EL2ORC * QESO(I,k-1) / (TO(I,k-1)**2)
786                   RFACT =  1. + DELTA * CP * GAMMA
787          &                 * TO(I,k-1) / HVAP
788                   AA1(I) = AA1(I) +
789          &                 DZ1 * (G / (CP * TO(I,k-1)))
790          &                 * DBYO(I,k-1) / (1. + GAMMA)
791          &                 * RFACT
792                   val = 0.
793                   AA1(I)=AA1(I)+
794          &                 DZ1 * G * DELTA *
795          &                 MAX(val,(QESO(I,k-1) - QO(I,k-1)))
796                 ENDIF
797               endif
798             ENDDO
799           ENDDO
800           DO I = 1, IM
801             IF(CNVFLG(I).AND.AA1(I).LE.0.) DWNFLG(I)  = .FALSE.
802             IF(CNVFLG(I).AND.AA1(I).LE.0.) DWNFLG2(I) = .FALSE.
803             IF(CNVFLG(I).AND.AA1(I).LE.0.) CNVFLG(I)  = .FALSE.
804           ENDDO
805     !!
806           TOTFLG = .TRUE.
807           DO I = 1, IM
808             TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
809           ENDDO
810           IF(TOTFLG) RETURN
811     !!
812     ccccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
813     ccccc   PRINT *, ' AA1(I) BEFORE DWNDRFT =', AA1(I)
814     ccccc ENDIF
815     C
816     C------- DOWNDRAFT CALCULATIONS
817     C
818     C
819     C--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
820     C
821           DO I = 1, IM
822             IF(CNVFLG(I)) THEN
823               VSHEAR(I) = 0.
824             ENDIF
825           ENDDO
826           DO K = 1, KM
827             DO I = 1, IM
828               if (k .le. kmax(i)) then
829                 IF(K.GE.KB(I).AND.K.LE.KTCON(I).AND.CNVFLG(I)) THEN
830                   shear     = sqrt((UO(I,k+1)-UO(I,k)) ** 2
831          &                       + (VO(I,k+1)-VO(I,k)) ** 2)
832                   VSHEAR(I) = VSHEAR(I) + SHEAR
833                 ENDIF
834               endif
835             ENDDO
836           ENDDO
837           DO I = 1, IM
838             EDT(I) = 0.
839             IF(CNVFLG(I)) THEN
840               KNUMB = KTCON(I) - KB(I) + 1
841               KNUMB = MAX(KNUMB,1)
842               VSHEAR(I) = 1.E3 * VSHEAR(I) / (ZO(I,KTCON(I))-ZO(I,KB(I)))
843               E1=1.591-.639*VSHEAR(I)
844          &       +.0953*(VSHEAR(I)**2)-.00496*(VSHEAR(I)**3)
845               EDT(I)=1.-E1
846               val =         .9
847               EDT(I) = MIN(EDT(I),val)
848               val =         .0
849               EDT(I) = MAX(EDT(I),val)
850               EDTO(I)=EDT(I)
851               EDTX(I)=EDT(I)
852             ENDIF
853           ENDDO
854     C  DETERMINE DETRAINMENT RATE BETWEEN 1 AND KBDTR
855           DO I = 1, IM
856             KBDTR(I) = KBCON(I)
857             beta = betas
858             if(SLIMSK(I).eq.1.) beta = betal
859             IF(CNVFLG(I)) THEN
860               KBDTR(I) = KBCON(I)
861               KBDTR(I) = MAX(KBDTR(I),1)
862               XLAMD(I) = 0.
863               IF(KBDTR(I).GT.1) THEN
864                 DZ = .5 * ZO(I,KBDTR(I)) + .5 * ZO(I,KBDTR(I)-1)
865          &         - ZO(I,1)
866                 XLAMD(I) =  LOG(BETA) / DZ
867               ENDIF
868             ENDIF
869           ENDDO
870     C  DETERMINE DOWNDRAFT MASS FLUX
871           DO K = 1, KM
872             DO I = 1, IM
873               IF(k .le. kmax(i)) then
874                 IF(CNVFLG(I)) THEN
875                   ETAD(I,k) = 1.
876                 ENDIF
877                 QRCDO(I,k) = 0.
878               endif
879             ENDDO
880           ENDDO
881           DO K = KM1, 2, -1
882             DO I = 1, IM
883               if (k .le. kbmax(i)) then
884                 IF(CNVFLG(I).AND.K.LT.KBDTR(I)) THEN
885                   DZ        = .5 * (ZO(I,k+1) - ZO(I,k-1))
886                   ETAD(I,k) = ETAD(I,k+1) * EXP(XLAMD(I) * DZ)
887                 ENDIF
888               endif
889             ENDDO
890           ENDDO
891           K = 1
892           DO I = 1, IM
893             IF(CNVFLG(I).AND.KBDTR(I).GT.1) THEN
894               DZ = .5 * (ZO(I,2) - ZO(I,1))
895               ETAD(I,k) = ETAD(I,k+1) * EXP(XLAMD(I) * DZ)
896             ENDIF
897           ENDDO
898     C
899     C--- DOWNDRAFT MOISTURE PROPERTIES
900     C
901           DO I = 1, IM
902             PWEVO(I) = 0.
903             FLG(I) = CNVFLG(I)
904           ENDDO
905           DO I = 1, IM
906             IF(CNVFLG(I)) THEN
907               JMN = JMIN(I)
908               HCDO(I) = HEO(I,JMN)
909               QCDO(I) = QO(I,JMN)
910               QRCDO(I,JMN) = QESO(I,JMN)
911               UCDO(I) = UO(I,JMN)
912               VCDO(I) = VO(I,JMN)
913             ENDIF
914           ENDDO
915           DO K = KM1, 1, -1
916             DO I = 1, IM
917               if (k .le. kmax(i)-1) then
918                 IF(CNVFLG(I).AND.K.LT.JMIN(I)) THEN
919                   DQ = QESO(I,k)
920                   DT = TO(I,k)
921                   GAMMA      = EL2ORC * DQ / DT**2
922                   DH         = HCDO(I) - HESO(I,k)
923                   QRCDO(I,k) = DQ+(1./HVAP)*(GAMMA/(1.+GAMMA))*DH
924                   DETAD      = ETAD(I,k+1) - ETAD(I,k)
925                   PWDO(I,k)  = ETAD(I,k+1) * QCDO(I) -
926          &                     ETAD(I,k) * QRCDO(I,k)
927                   PWDO(I,k)  = PWDO(I,k) - DETAD *
928          &                    .5 * (QRCDO(I,k) + QRCDO(I,k+1))
929                   QCDO(I)    = QRCDO(I,k)
930                   PWEVO(I)   = PWEVO(I) + PWDO(I,k)
931                 ENDIF
932               endif
933             ENDDO
934           ENDDO
935     C     IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG(I)) THEN
936     C       PRINT *, ' PWAVO(I), PWEVO(I) =', PWAVO(I), PWEVO(I)
937     C     ENDIF
938     C
939     C--- FINAL DOWNDRAFT STRENGTH DEPENDENT ON PRECIP
940     C--- EFFICIENCY (EDT), NORMALIZED CONDENSATE (PWAV), AND
941     C--- EVAPORATE (PWEV)
942     C
943           DO I = 1, IM
944             edtmax = edtmaxl
945             if(SLIMSK(I).eq.0.) edtmax = edtmaxs
946             IF(DWNFLG2(I)) THEN
947               IF(PWEVO(I).LT.0.) THEN
948                 EDTO(I) = -EDTO(I) * PWAVO(I) / PWEVO(I)
949                 EDTO(I) = MIN(EDTO(I),EDTMAX)
950               ELSE
951                 EDTO(I) = 0.
952               ENDIF
953             ELSE
954               EDTO(I) = 0.
955             ENDIF
956           ENDDO
957     C
958     C
959     C--- DOWNDRAFT CLOUDWORK FUNCTIONS
960     C
961     C
962           DO K = KM1, 1, -1
963             DO I = 1, IM
964               if (k .le. kmax(i)-1) then
965                 IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
966                   GAMMA = EL2ORC * QESO(I,k+1) / TO(I,k+1)**2
967                   DHH=HCDO(I)
968                   DT=TO(I,k+1)
969                   DG=GAMMA
970                   DH=HESO(I,k+1)
971                   DZ=-1.*(ZO(I,k+1)-ZO(I,k))
972                   AA1(I)=AA1(I)+EDTO(I)*DZ*(G/(CP*DT))*((DHH-DH)/(1.+DG))
973          &               *(1.+DELTA*CP*DG*DT/HVAP)
974                   val=0.
975                   AA1(I)=AA1(I)+EDTO(I)*
976          &        DZ*G*DELTA*MAX(val,(QESO(I,k+1)-QO(I,k+1)))
977                 ENDIF
978               endif
979             ENDDO
980           ENDDO
981     ccccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG2(I)) THEN
982     ccccc   PRINT *, '  AA1(I) AFTER DWNDRFT =', AA1(I)
983     ccccc ENDIF
984           DO I = 1, IM
985             IF(AA1(I).LE.0.) CNVFLG(I)  = .FALSE.
986             IF(AA1(I).LE.0.) DWNFLG(I)  = .FALSE.
987             IF(AA1(I).LE.0.) DWNFLG2(I) = .FALSE.
988           ENDDO
989     !!
990           TOTFLG = .TRUE.
991           DO I = 1, IM
992             TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
993           ENDDO
994           IF(TOTFLG) RETURN
995     !!
996     C
997     C
998     C--- WHAT WOULD THE CHANGE BE, THAT A CLOUD WITH UNIT MASS
999     C--- WILL DO TO THE ENVIRONMENT?
1000     C
1001           DO K = 1, KM
1002             DO I = 1, IM
1003               IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1004                 DELLAH(I,k) = 0.
1005                 DELLAQ(I,k) = 0.
1006                 DELLAU(I,k) = 0.
1007                 DELLAV(I,k) = 0.
1008               ENDIF
1009             ENDDO
1010           ENDDO
1011           DO I = 1, IM
1012             IF(CNVFLG(I)) THEN
1013               DP = 1000. * DEL(I,1)
1014               DELLAH(I,1) = EDTO(I) * ETAD(I,1) * (HCDO(I)
1015          &                - HEO(I,1)) * G / DP
1016               DELLAQ(I,1) = EDTO(I) * ETAD(I,1) * (QCDO(I)
1017          &                - QO(I,1)) * G / DP
1018               DELLAU(I,1) = EDTO(I) * ETAD(I,1) * (UCDO(I)
1019          &                - UO(I,1)) * G / DP
1020               DELLAV(I,1) = EDTO(I) * ETAD(I,1) * (VCDO(I)
1021          &                - VO(I,1)) * G / DP
1022             ENDIF
1023           ENDDO
1024     C
1025     C--- CHANGED DUE TO SUBSIDENCE AND ENTRAINMENT
1026     C
1027           DO K = 2, KM1
1028             DO I = 1, IM
1029               if (k .le. kmax(i)-1) then
1030                 IF(CNVFLG(I).AND.K.LT.KTCON(I)) THEN
1031                   AUP = 1.
1032                   IF(K.LE.KB(I)) AUP = 0.
1033                   ADW = 1.
1034                   IF(K.GT.JMIN(I)) ADW = 0.
1035                   DV1= HEO(I,k)
1036                   DV2 = .5 * (HEO(I,k) + HEO(I,k+1))
1037                   DV3= HEO(I,k-1)
1038                   DV1Q= QO(I,k)
1039                   DV2Q = .5 * (QO(I,k) + QO(I,k+1))
1040                   DV3Q= QO(I,k-1)
1041                   DV1U= UO(I,k)
1042                   DV2U = .5 * (UO(I,k) + UO(I,k+1))
1043                   DV3U= UO(I,k-1)
1044                   DV1V= VO(I,k)
1045                   DV2V = .5 * (VO(I,k) + VO(I,k+1))
1046                   DV3V= VO(I,k-1)
1047                   DP = 1000. * DEL(I,K)
1048                   DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1049                   DETA = ETA(I,k) - ETA(I,k-1)
1050                   DETAD = ETAD(I,k) - ETAD(I,k-1)
1051                   DELLAH(I,k) = DELLAH(I,k) +
1052          &            ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1
1053          &        - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3
1054          &                    - AUP * DETA * DV2
1055          &                    + ADW * EDTO(I) * DETAD * HCDO(I)) * G / DP
1056                   DELLAQ(I,k) = DELLAQ(I,k) +
1057          &            ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1Q
1058          &        - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3Q
1059          &                    - AUP * DETA * DV2Q
1060          &       +ADW*EDTO(I)*DETAD*.5*(QRCDO(I,k)+QRCDO(I,k-1))) * G / DP
1061                   DELLAU(I,k) = DELLAU(I,k) +
1062          &            ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1U
1063          &        - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3U
1064          &                     - AUP * DETA * DV2U
1065          &                    + ADW * EDTO(I) * DETAD * UCDO(I)
1066          &                    ) * G / DP
1067                   DELLAV(I,k) = DELLAV(I,k) +
1068          &            ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1V
1069          &        - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3V
1070          &                     - AUP * DETA * DV2V
1071          &                    + ADW * EDTO(I) * DETAD * VCDO(I)
1072          &                    ) * G / DP
1073                 ENDIF
1074               endif
1075             ENDDO
1076           ENDDO
1077     C
1078     C------- CLOUD TOP
1079     C
1080           DO I = 1, IM
1081             IF(CNVFLG(I)) THEN
1082               INDX = KTCON(I)
1083               DP = 1000. * DEL(I,INDX)
1084               DV1 = HEO(I,INDX-1)
1085               DELLAH(I,INDX) = ETA(I,INDX-1) *
1086          &                     (HCKO(I,INDX-1) - DV1) * G / DP
1087               DVQ1 = QO(I,INDX-1)
1088               DELLAQ(I,INDX) = ETA(I,INDX-1) *
1089          &                     (QCKO(I,INDX-1) - DVQ1) * G / DP
1090               DV1U = UO(I,INDX-1)
1091               DELLAU(I,INDX) = ETA(I,INDX-1) *
1092          &                     (UCKO(I,INDX-1) - DV1U) * G / DP
1093               DV1V = VO(I,INDX-1)
1094               DELLAV(I,INDX) = ETA(I,INDX-1) *
1095          &                     (VCKO(I,INDX-1) - DV1V) * G / DP
1096     c
1097     c  cloud water
1098     c
1099               DELLAL(I) = ETA(I,INDX-1) * QLKO_KTCON(I) * g / dp
1100             ENDIF
1101           ENDDO
1102     C
1103     C------- FINAL CHANGED VARIABLE PER UNIT MASS FLUX
1104     C
1105           DO K = 1, KM
1106             DO I = 1, IM
1107               if (k .le. kmax(i)) then
1108                 IF(CNVFLG(I).and.k.gt.KTCON(I)) THEN
1109                   QO(I,k) = Q1(I,k)
1110                   TO(I,k) = T1(I,k)
1111                   UO(I,k) = U1(I,k)
1112                   VO(I,k) = V1(I,k)
1113                 ENDIF
1114                 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
1115                   QO(I,k) = DELLAQ(I,k) * MBDT + Q1(I,k)
1116                   DELLAT = (DELLAH(I,k) - HVAP * DELLAQ(I,k)) / CP
1117                   TO(I,k) = DELLAT * MBDT + T1(I,k)
1118                   val   =           1.e-10
1119                   QO(I,k) = max(QO(I,k), val  )
1120                 ENDIF
1121               endif
1122             ENDDO
1123           ENDDO
1124     C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1125     C
1126     C--- THE ABOVE CHANGED ENVIRONMENT IS NOW USED TO CALULATE THE
1127     C--- EFFECT THE ARBITRARY CLOUD (WITH UNIT MASS FLUX)
1128     C--- WOULD HAVE ON THE STABILITY,
1129     C--- WHICH THEN IS USED TO CALCULATE THE REAL MASS FLUX,
1130     C--- NECESSARY TO KEEP THIS CHANGE IN BALANCE WITH THE LARGE-SCALE
1131     C--- DESTABILIZATION.
1132     C
1133     C--- ENVIRONMENTAL CONDITIONS AGAIN, FIRST HEIGHTS
1134     C
1135           DO K = 1, KM
1136             DO I = 1, IM
1137               IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1138     !jfe        QESO(I,k) = 10. * FPVS(TO(I,k))
1139     !
1140                 QESO(I,k) = 0.01 * fpvs(TO(I,K))      ! fpvs is in Pa
1141     !
1142                 QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k)+EPSM1*QESO(I,k))
1143                 val       =             1.E-8
1144                 QESO(I,k) = MAX(QESO(I,k), val )
1145                 TVO(I,k)  = TO(I,k) + DELTA * TO(I,k) * QO(I,k)
1146               ENDIF
1147             ENDDO
1148           ENDDO
1149           DO I = 1, IM
1150             IF(CNVFLG(I)) THEN
1151               XAA0(I) = 0.
1152               XPWAV(I) = 0.
1153             ENDIF
1154           ENDDO
1155     C
1156     C  HYDROSTATIC HEIGHT ASSUME ZERO TERR
1157     C
1158     !     DO I = 1, IM
1159     !       IF(CNVFLG(I)) THEN
1160     !         DLNSIG =  LOG(PRSL(I,1)/PS(I))
1161     !         ZO(I,1) = TERR - DLNSIG * RD / G * TVO(I,1)
1162     !       ENDIF
1163     !     ENDDO
1164     !     DO K = 2, KM
1165     !       DO I = 1, IM
1166     !         IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1167     !           DLNSIG =  LOG(PRSL(I,K) / PRSL(I,K-1))
1168     !           ZO(I,k) = ZO(I,k-1) - DLNSIG * RD / G
1169     !    &             * .5 * (TVO(I,k) + TVO(I,k-1))
1170     !         ENDIF
1171     !       ENDDO
1172     !     ENDDO
1173     C
1174     C--- MOIST STATIC ENERGY
1175     C
1176           DO K = 1, KM1
1177             DO I = 1, IM
1178               IF(k .le. kmax(i)-1 .and. CNVFLG(I)) THEN
1179                 DZ = .5 * (ZO(I,k+1) - ZO(I,k))
1180                 DP = .5 * (PFLD(I,k+1) - PFLD(I,k))
1181     cjfe        ES = 10. * FPVS(TO(I,k+1))
1182     !
1183                 ES = 0.01 * fpvs(TO(I,K+1))      ! fpvs is in Pa
1184     !
1185                 PPRIME = PFLD(I,k+1) + EPSM1 * ES
1186                 QS = EPS * ES / PPRIME
1187                 DQSDP = - QS / PPRIME
1188                 DESDT = ES * (FACT1 / TO(I,k+1) + FACT2 / (TO(I,k+1)**2))
1189                 DQSDT = QS * PFLD(I,k+1) * DESDT / (ES * PPRIME)
1190                 GAMMA = EL2ORC * QESO(I,k+1) / (TO(I,k+1)**2)
1191                 DT = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
1192                 DQ = DQSDT * DT + DQSDP * DP
1193                 TO(I,k) = TO(I,k+1) + DT
1194                 QO(I,k) = QO(I,k+1) + DQ
1195                 PO(I,k) = .5 * (PFLD(I,k) + PFLD(I,k+1))
1196               ENDIF
1197             ENDDO
1198           ENDDO
1199           DO K = 1, KM1
1200             DO I = 1, IM
1201               IF(k .le. kmax(i)-1 .and. CNVFLG(I)) THEN
1202     cjfe        QESO(I,k) = 10. * FPVS(TO(I,k))
1203     !
1204                 QESO(I,k) = 0.01 * fpvs(TO(I,K))      ! fpvs is in Pa
1205     !
1206                 QESO(I,k) = EPS * QESO(I,k) / (PO(I,k) + EPSM1 * QESO(I,k))
1207                 val1      =             1.E-8
1208                 QESO(I,k) = MAX(QESO(I,k), val1)
1209                 val2      =           1.e-10
1210                 QO(I,k)   = max(QO(I,k), val2 )
1211     c           QO(I,k)   = MIN(QO(I,k),QESO(I,k))
1212                 HEO(I,k)   = .5 * G * (ZO(I,k) + ZO(I,k+1)) +
1213          &                    CP * TO(I,k) + HVAP * QO(I,k)
1214                 HESO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) +
1215          &                  CP * TO(I,k) + HVAP * QESO(I,k)
1216               ENDIF
1217             ENDDO
1218           ENDDO
1219           DO I = 1, IM
1220             k = kmax(i)
1221             IF(CNVFLG(I)) THEN
1222               HEO(I,k) = G * ZO(I,k) + CP * TO(I,k) + HVAP * QO(I,k)
1223               HESO(I,k) = G * ZO(I,k) + CP * TO(I,k) + HVAP * QESO(I,k)
1224     c         HEO(I,k) = MIN(HEO(I,k),HESO(I,k))
1225             ENDIF
1226           ENDDO
1227           DO I = 1, IM
1228             IF(CNVFLG(I)) THEN
1229               INDX = KB(I)
1230               XHKB(I) = HEO(I,INDX)
1231               XQKB(I) = QO(I,INDX)
1232               HCKO(I,INDX) = XHKB(I)
1233               QCKO(I,INDX) = XQKB(I)
1234             ENDIF
1235           ENDDO
1236     C
1237     C
1238     C**************************** STATIC CONTROL
1239     C
1240     C
1241     C------- MOISTURE AND CLOUD WORK FUNCTIONS
1242     C
1243           DO K = 2, KM1
1244             DO I = 1, IM
1245               if (k .le. kmax(i)-1) then
1246     C           IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KBCON(I)) THEN
1247                 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KTCON(I)) THEN
1248                   FACTOR = ETA(I,k-1) / ETA(I,k)
1249                   ONEMF = 1. - FACTOR
1250                   HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF *
1251          &                    .5 * (HEO(I,k) + HEO(I,k+1))
1252                 ENDIF
1253     C           IF(CNVFLG(I).AND.K.GT.KBCON(I)) THEN
1254     C             HEO(I,k) = HEO(I,k-1)
1255     C           ENDIF
1256               endif
1257             ENDDO
1258           ENDDO
1259           DO K = 2, KM1
1260             DO I = 1, IM
1261               if (k .le. kmax(i)-1) then
1262                 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LT.KTCON(I)) THEN
1263                   DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1264                   GAMMA = EL2ORC * QESO(I,k) / (TO(I,k)**2)
1265                   XDBY = HCKO(I,k) - HESO(I,k)
1266                   val  =          0.
1267                   XDBY = MAX(XDBY,val)
1268                   XQRCH = QESO(I,k)
1269          &              + GAMMA * XDBY / (HVAP * (1. + GAMMA))
1270                   FACTOR = ETA(I,k-1) / ETA(I,k)
1271                   ONEMF = 1. - FACTOR
1272                   QCKO(I,k) = FACTOR * QCKO(I,k-1) + ONEMF *
1273          &                    .5 * (QO(I,k) + QO(I,k+1))
1274                   DQ = ETA(I,k) * QCKO(I,k) - ETA(I,k) * XQRCH
1275                   IF(DQ.GT.0.) THEN
1276                     ETAH = .5 * (ETA(I,k) + ETA(I,k-1))
1277                     QLK = DQ / (ETA(I,k) + ETAH * C0 * DZ)
1278                     XAA0(I) = XAA0(I) - (ZO(I,k) - ZO(I,k-1)) * G * QLK
1279                     XQC = QLK + XQRCH
1280                     XPW = ETAH * C0 * DZ * QLK
1281                     QCKO(I,k) = XQC
1282                     XPWAV(I) = XPWAV(I) + XPW
1283                   ENDIF
1284                 ENDIF
1285     c           IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LT.KTCON(I)) THEN
1286                 IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1287                   DZ1 = ZO(I,k) - ZO(I,k-1)
1288                   GAMMA = EL2ORC * QESO(I,k-1) / (TO(I,k-1)**2)
1289                   RFACT =  1. + DELTA * CP * GAMMA
1290          &                 * TO(I,k-1) / HVAP
1291                   XDBY = HCKO(I,k-1) - HESO(I,k-1)
1292                   XAA0(I) = XAA0(I)
1293          &                + DZ1 * (G / (CP * TO(I,k-1)))
1294          &                * XDBY / (1. + GAMMA)
1295          &                * RFACT
1296                   val=0.
1297                   XAA0(I)=XAA0(I)+
1298          &                 DZ1 * G * DELTA *
1299          &                 MAX(val,(QESO(I,k-1) - QO(I,k-1)))
1300                 ENDIF
1301               endif
1302             ENDDO
1303           ENDDO
1304     ccccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1305     ccccc   PRINT *, ' XAA BEFORE DWNDRFT =', XAA0(I)
1306     ccccc ENDIF
1307     C
1308     C------- DOWNDRAFT CALCULATIONS
1309     C
1310     C
1311     C--- DOWNDRAFT MOISTURE PROPERTIES
1312     C
1313           DO I = 1, IM
1314             XPWEV(I) = 0.
1315           ENDDO
1316           DO I = 1, IM
1317             IF(DWNFLG2(I)) THEN
1318               JMN = JMIN(I)
1319               XHCD(I) = HEO(I,JMN)
1320               XQCD(I) = QO(I,JMN)
1321               QRCD(I,JMN) = QESO(I,JMN)
1322             ENDIF
1323           ENDDO
1324           DO K = KM1, 1, -1
1325             DO I = 1, IM
1326               if (k .le. kmax(i)-1) then
1327                 IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
1328                   DQ = QESO(I,k)
1329                   DT = TO(I,k)
1330                   GAMMA    = EL2ORC * DQ / DT**2
1331                   DH       = XHCD(I) - HESO(I,k)
1332                   QRCD(I,k)=DQ+(1./HVAP)*(GAMMA/(1.+GAMMA))*DH
1333                   DETAD    = ETAD(I,k+1) - ETAD(I,k)
1334                   XPWD     = ETAD(I,k+1) * QRCD(I,k+1) -
1335          &                   ETAD(I,k) * QRCD(I,k)
1336                   XPWD     = XPWD - DETAD *
1337          &                 .5 * (QRCD(I,k) + QRCD(I,k+1))
1338                   XPWEV(I) = XPWEV(I) + XPWD
1339                 ENDIF
1340               endif
1341             ENDDO
1342           ENDDO
1343     C
1344           DO I = 1, IM
1345             edtmax = edtmaxl
1346             if(SLIMSK(I).eq.0.) edtmax = edtmaxs
1347             IF(DWNFLG2(I)) THEN
1348               IF(XPWEV(I).GE.0.) THEN
1349                 EDTX(I) = 0.
1350               ELSE
1351                 EDTX(I) = -EDTX(I) * XPWAV(I) / XPWEV(I)
1352                 EDTX(I) = MIN(EDTX(I),EDTMAX)
1353               ENDIF
1354             ELSE
1355               EDTX(I) = 0.
1356             ENDIF
1357           ENDDO
1358     C
1359     C
1360     C
1361     C--- DOWNDRAFT CLOUDWORK FUNCTIONS
1362     C
1363     C
1364           DO K = KM1, 1, -1
1365             DO I = 1, IM
1366               if (k .le. kmax(i)-1) then
1367                 IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
1368                   GAMMA = EL2ORC * QESO(I,k+1) / TO(I,k+1)**2
1369                   DHH=XHCD(I)
1370                   DT= TO(I,k+1)
1371                   DG= GAMMA
1372                   DH= HESO(I,k+1)
1373                   DZ=-1.*(ZO(I,k+1)-ZO(I,k))
1374                   XAA0(I)=XAA0(I)+EDTX(I)*DZ*(G/(CP*DT))*((DHH-DH)/(1.+DG))
1375          &                *(1.+DELTA*CP*DG*DT/HVAP)
1376                   val=0.
1377                   XAA0(I)=XAA0(I)+EDTX(I)*
1378          &        DZ*G*DELTA*MAX(val,(QESO(I,k+1)-QO(I,k+1)))
1379                 ENDIF
1380               endif
1381             ENDDO
1382           ENDDO
1383     ccccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG2(I)) THEN
1384     ccccc   PRINT *, '  XAA AFTER DWNDRFT =', XAA0(I)
1385     ccccc ENDIF
1386     C
1387     C  CALCULATE CRITICAL CLOUD WORK FUNCTION
1388     C
1389           DO I = 1, IM
1390             ACRT(I) = 0.
1391             IF(CNVFLG(I)) THEN
1392     C       IF(CNVFLG(I).AND.SLIMSK(I).NE.1.) THEN
1393               IF(PFLD(I,KTCON(I)).LT.PCRIT(15))THEN
1394                 ACRT(I)=ACRIT(15)*(975.-PFLD(I,KTCON(I)))
1395          &              /(975.-PCRIT(15))
1396               ELSE IF(PFLD(I,KTCON(I)).GT.PCRIT(1))THEN
1397                 ACRT(I)=ACRIT(1)
1398               ELSE
1399                 K =  int((850. - PFLD(I,KTCON(I)))/50.) + 2
1400                 K = MIN(K,15)
1401                 K = MAX(K,2)
1402                 ACRT(I)=ACRIT(K)+(ACRIT(K-1)-ACRIT(K))*
1403          *           (PFLD(I,KTCON(I))-PCRIT(K))/(PCRIT(K-1)-PCRIT(K))
1404                ENDIF
1405     C        ELSE
1406     C          ACRT(I) = .5 * (PFLD(I,KBCON(I)) - PFLD(I,KTCON(I)))
1407              ENDIF
1408           ENDDO
1409           DO I = 1, IM
1410             ACRTFCT(I) = 1.
1411             IF(CNVFLG(I)) THEN
1412               if(SLIMSK(I).eq.1.) THEN
1413                 w1 = w1l
1414                 w2 = w2l
1415                 w3 = w3l
1416                 w4 = w4l
1417               else
1418                 w1 = w1s
1419                 w2 = w2s
1420                 w3 = w3s
1421                 w4 = w4s
1422               ENDIF
1423     C       IF(CNVFLG(I).AND.SLIMSK(I).EQ.1.) THEN
1424     C         ACRTFCT(I) = PDOT(I) / W3
1425     c
1426     c  modify critical cloud workfunction by cloud base vertical velocity
1427     c
1428               IF(PDOT(I).LE.W4) THEN
1429                 ACRTFCT(I) = (PDOT(I) - W4) / (W3 - W4)
1430               ELSEIF(PDOT(I).GE.-W4) THEN
1431                 ACRTFCT(I) = - (PDOT(I) + W4) / (W4 - W3)
1432               ELSE
1433                 ACRTFCT(I) = 0.
1434               ENDIF
1435               val1    =             -1.
1436               ACRTFCT(I) = MAX(ACRTFCT(I),val1)
1437               val2    =             1.
1438               ACRTFCT(I) = MIN(ACRTFCT(I),val2)
1439               ACRTFCT(I) = 1. - ACRTFCT(I)
1440     c
1441     c  modify ACRTFCT(I) by colume mean rh if RHBAR(I) is greater than 80 percent
1442     c
1443     c         if(RHBAR(I).ge..8) THEN
1444     c           ACRTFCT(I) = ACRTFCT(I) * (.9 - min(RHBAR(I),.9)) * 10.
1445     c         ENDIF
1446     c
1447     c  modify adjustment time scale by cloud base vertical velocity
1448     c
1449               DTCONV(I) = DT2 + max((1800. - DT2),cons_0) *
1450          &                (PDOT(I) - W2) / (W1 - W2)
1451     c         DTCONV(I) = MAX(DTCONV(I), DT2)
1452     c         DTCONV(I) = 1800. * (PDOT(I) - w2) / (w1 - w2)
1453               DTCONV(I) = max(DTCONV(I),dtmin)
1454               DTCONV(I) = min(DTCONV(I),dtmax)
1455      
1456             ENDIF
1457           ENDDO
1458     C
1459     C--- LARGE SCALE FORCING
1460     C
1461           DO I= 1, IM
1462             FLG(I) = CNVFLG(I)
1463             IF(CNVFLG(I)) THEN
1464     C         F = AA1(I) / DTCONV(I)
1465               FLD(I) = (AA1(I) - ACRT(I) * ACRTFCT(I)) / DTCONV(I)
1466               IF(FLD(I).LE.0.) FLG(I) = .FALSE.
1467             ENDIF
1468             CNVFLG(I) = FLG(I)
1469             IF(CNVFLG(I)) THEN
1470     C         XAA0(I) = MAX(XAA0(I),0.)
1471               XK(I) = (XAA0(I) - AA1(I)) / MBDT
1472               IF(XK(I).GE.0.) FLG(I) = .FALSE.
1473             ENDIF
1474     C
1475     C--- KERNEL, CLOUD BASE MASS FLUX
1476     C
1477             CNVFLG(I) = FLG(I)
1478             IF(CNVFLG(I)) THEN
1479               XMB(I) = -FLD(I) / XK(I)
1480               XMB(I) = MIN(XMB(I),XMBMAX(I))
1481             ENDIF
1482           ENDDO
1483     c      IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1484     c        print *, ' RHBAR(I), ACRTFCT(I) =', RHBAR(I), ACRTFCT(I)
1485     c        PRINT *, '  A1, XA =', AA1(I), XAA0(I)
1486     c        PRINT *, ' XMB(I), ACRT =', XMB(I), ACRT
1487     c      ENDIF
1488           TOTFLG = .TRUE.
1489           DO I = 1, IM
1490             TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1491           ENDDO
1492           IF(TOTFLG) RETURN
1493     c
1494     c  restore t0 and QO to t1 and q1 in case convection stops
1495     c
1496           do k = 1, km
1497             DO I = 1, IM
1498               if (k .le. kmax(i)) then
1499                 TO(I,k) = T1(I,k)
1500                 QO(I,k) = Q1(I,k)
1501     !jfe        QESO(I,k) = 10. * FPVS(T1(I,k))
1502     !
1503                 QESO(I,k) = 0.01 * fpvs(T1(I,K))      ! fpvs is in Pa
1504     !
1505                 QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k) + EPSM1*QESO(I,k))
1506                 val     =             1.E-8
1507                 QESO(I,k) = MAX(QESO(I,k), val )
1508               endif
1509             enddo
1510           enddo
1511     C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1512     C
1513     C--- FEEDBACK: SIMPLY THE CHANGES FROM THE CLOUD WITH UNIT MASS FLUX
1514     C---           MULTIPLIED BY  THE MASS FLUX NECESSARY TO KEEP THE
1515     C---           EQUILIBRIUM WITH THE LARGER-SCALE.
1516     C
1517           DO I = 1, IM
1518             DELHBAR(I) = 0.
1519             DELQBAR(I) = 0.
1520             DELTBAR(I) = 0.
1521             QCOND(I) = 0.
1522           ENDDO
1523           DO K = 1, KM
1524             DO I = 1, IM
1525               if (k .le. kmax(i)) then
1526                 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
1527                   AUP = 1.
1528                   IF(K.Le.KB(I)) AUP = 0.
1529                   ADW = 1.
1530                   IF(K.GT.JMIN(I)) ADW = 0.
1531                   DELLAT = (DELLAH(I,k) - HVAP * DELLAQ(I,k)) / CP
1532                   T1(I,k) = T1(I,k) + DELLAT * XMB(I) * DT2
1533                   Q1(I,k) = Q1(I,k) + DELLAQ(I,k) * XMB(I) * DT2
1534                   U1(I,k) = U1(I,k) + DELLAU(I,k) * XMB(I) * DT2
1535                   V1(I,k) = V1(I,k) + DELLAV(I,k) * XMB(I) * DT2
1536                   DP = 1000. * DEL(I,K)
1537                   DELHBAR(I) = DELHBAR(I) + DELLAH(I,k)*XMB(I)*DP/G
1538                   DELQBAR(I) = DELQBAR(I) + DELLAQ(I,k)*XMB(I)*DP/G
1539                   DELTBAR(I) = DELTBAR(I) + DELLAT*XMB(I)*DP/G
1540                 ENDIF
1541               endif
1542             ENDDO
1543           ENDDO
1544           DO K = 1, KM
1545             DO I = 1, IM
1546               if (k .le. kmax(i)) then
1547                 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
1548     !jfe          QESO(I,k) = 10. * FPVS(T1(I,k))
1549     !
1550                   QESO(I,k) = 0.01 * fpvs(T1(I,K))      ! fpvs is in Pa
1551     !
1552                   QESO(I,k) = EPS * QESO(I,k)/(PFLD(I,k) + EPSM1*QESO(I,k))
1553                   val     =             1.E-8
1554                   QESO(I,k) = MAX(QESO(I,k), val )
1555     c
1556     c  cloud water
1557     c
1558                   if(ncloud.gt.0.and.cnvflg(i).and.k.eq.ktcon(i)) then
1559                     tem  = dellal(i) * xmb(i) * dt2
1560                     tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
1561                     if (ql(i,k,2) .gt. -999.0) then
1562                       ql(i,k,1) = ql(i,k,1) + tem * tem1            ! ice
1563                       ql(i,k,2) = ql(i,k,2) + tem *(1.0-tem1)       ! water
1564                     else
1565                       ql(i,k,1) = ql(i,k,1) + tem
1566                     endif
1567                     dp = 1000. * del(i,k)
1568                     dellal(i) = dellal(i) * xmb(i) * dp / g
1569                   endif
1570     !
1571     !             if(ncloud.gt.0.and.CNVFLG(I).and.k.eq.KTCON(I)) THEN
1572     !               QL(I,k) = QL(I,k) + DELLAL(I) * XMB(I) * dt2
1573     !               dp = 1000. * del(i,k)
1574     !               DELLAL(I) = DELLAL(I) * XMB(I) * dp / g
1575     !             ENDIF
1576     !
1577                 ENDIF
1578               endif
1579             ENDDO
1580           ENDDO
1581     c     IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I) ) THEN
1582     c       PRINT *, ' DELHBAR, DELQBAR, DELTBAR ='
1583     c       PRINT *, DELHBAR, HVAP*DELQBAR, CP*DELTBAR
1584     c       PRINT *, '   DELLBAR ='
1585     c       PRINT 6003,  HVAP*DELLbar
1586     c       PRINT *, '   DELLAQ ='
1587     c       PRINT 6003, (HVAP*DELLAQ(I,k)*XMB(I),K=1,KMAX)
1588     c       PRINT *, '   DELLAT ='
1589     c       PRINT 6003, (DELLAH(i,k)*XMB(I)-HVAP*DELLAQ(I,k)*XMB(I),
1590     c    &               K=1,KMAX)
1591     c     ENDIF
1592           DO I = 1, IM
1593             RNTOT(I) = 0.
1594             DELQEV(I) = 0.
1595             DELQ2(I) = 0.
1596             FLG(I) = CNVFLG(I)
1597           ENDDO
1598           DO K = KM, 1, -1
1599             DO I = 1, IM
1600               if (k .le. kmax(i)) then
1601                 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
1602                   AUP = 1.
1603                   IF(K.Le.KB(I)) AUP = 0.
1604                   ADW = 1.
1605                   IF(K.GT.JMIN(I)) ADW = 0.
1606                   rain =  AUP * PWO(I,k) + ADW * EDTO(I) * PWDO(I,k)
1607                   RNTOT(I) = RNTOT(I) + rain * XMB(I) * .001 * dt2
1608                 ENDIF
1609               endif
1610             ENDDO
1611           ENDDO
1612           DO K = KM, 1, -1
1613             DO I = 1, IM
1614               if (k .le. kmax(i)) then
1615                 DELTV(I) = 0.
1616                 DELQ(I) = 0.
1617                 QEVAP(I) = 0.
1618                 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
1619                   AUP = 1.
1620                   IF(K.Le.KB(I)) AUP = 0.
1621                   ADW = 1.
1622                   IF(K.GT.JMIN(I)) ADW = 0.
1623                   rain =  AUP * PWO(I,k) + ADW * EDTO(I) * PWDO(I,k)
1624                   RN(I) = RN(I) + rain * XMB(I) * .001 * dt2
1625                 ENDIF
1626                 IF(FLG(I).AND.K.LE.KTCON(I)) THEN
1627                   evef = EDT(I) * evfact
1628                   if(SLIMSK(I).eq.1.) evef=EDT(I) * evfactl
1629     !             if(SLIMSK(I).eq.1.) evef=.07
1630     c             if(SLIMSK(I).ne.1.) evef = 0.
1631                   QCOND(I) = EVEF * (Q1(I,k) - QESO(I,k))
1632          &                 / (1. + EL2ORC * QESO(I,k) / T1(I,k)**2)
1633                   DP = 1000. * DEL(I,K)
1634                   IF(RN(I).GT.0..AND.QCOND(I).LT.0.) THEN
1635                     QEVAP(I) = -QCOND(I) * (1.-EXP(-.32*SQRT(DT2*RN(I))))
1636                     QEVAP(I) = MIN(QEVAP(I), RN(I)*1000.*G/DP)
1637                     DELQ2(I) = DELQEV(I) + .001 * QEVAP(I) * dp / g
1638                   ENDIF
1639                   if(RN(I).gt.0..and.QCOND(I).LT.0..and.
1640          &           DELQ2(I).gt.RNTOT(I)) THEN
1641                     QEVAP(I) = 1000.* g * (RNTOT(I) - DELQEV(I)) / dp
1642                     FLG(I) = .false.
1643                   ENDIF
1644                   IF(RN(I).GT.0..AND.QEVAP(I).gt.0.) THEN
1645                     Q1(I,k) = Q1(I,k) + QEVAP(I)
1646                     T1(I,k) = T1(I,k) - ELOCP * QEVAP(I)
1647                     RN(I) = RN(I) - .001 * QEVAP(I) * DP / G
1648                     DELTV(I) = - ELOCP*QEVAP(I)/DT2
1649                     DELQ(I) =  + QEVAP(I)/DT2
1650                     DELQEV(I) = DELQEV(I) + .001*dp*QEVAP(I)/g
1651                   ENDIF
1652                   DELLAQ(I,k) = DELLAQ(I,k) + DELQ(I) / XMB(I)
1653                   DELQBAR(I) = DELQBAR(I) + DELQ(I)*DP/G
1654                   DELTBAR(I) = DELTBAR(I) + DELTV(I)*DP/G
1655                 ENDIF
1656               endif
1657             ENDDO
1658           ENDDO
1659     c      IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I) ) THEN
1660     c        PRINT *, '   DELLAH ='
1661     c        PRINT 6003, (DELLAH(k)*XMB(I),K=1,KMAX)
1662     c        PRINT *, '   DELLAQ ='
1663     c        PRINT 6003, (HVAP*DELLAQ(I,k)*XMB(I),K=1,KMAX)
1664     c        PRINT *, ' DELHBAR, DELQBAR, DELTBAR ='
1665     c        PRINT *, DELHBAR, HVAP*DELQBAR, CP*DELTBAR
1666     c        PRINT *, ' PRECIP =', HVAP*RN(I)*1000./DT2
1667     CCCCC   PRINT *, '   DELLBAR ='
1668     CCCCC   PRINT *,  HVAP*DELLbar
1669     c      ENDIF
1670     C
1671     C  PRECIPITATION RATE CONVERTED TO ACTUAL PRECIP
1672     C  IN UNIT OF M INSTEAD OF KG
1673     C
1674           DO I = 1, IM
1675             IF(CNVFLG(I)) THEN
1676     C
1677     C  IN THE EVENT OF UPPER LEVEL RAIN EVAPORATION AND LOWER LEVEL DOWNDRAF
1678     C    MOISTENING, RN CAN BECOME NEGATIVE, IN THIS CASE, WE BACK OUT OF TH
1679     C    HEATING AND THE MOISTENING
1680     C
1681               if(RN(I).lt.0..and..not.FLG(I)) RN(I) = 0.
1682               IF(RN(I).LE.0.) THEN
1683                 RN(I) = 0.
1684               ELSE
1685                 KTOP(I) = KTCON(I)
1686                 KBOT(I) = KBCON(I)
1687                 KUO(I) = 1
1688                 CLDWRK(I) = AA1(I)
1689               ENDIF
1690             ENDIF
1691           ENDDO
1692           DO K = 1, KM
1693             DO I = 1, IM
1694               if (k .le. kmax(i)) then
1695                 IF(CNVFLG(I).AND.RN(I).LE.0.) THEN
1696                   T1(I,k) = TO(I,k)
1697                   Q1(I,k) = QO(I,k)
1698                 ENDIF
1699               endif
1700             ENDDO
1701           ENDDO
1702     ! hchuang code change [+24L]
1703           DO K = 1, KM
1704             DO I = 1, IM
1705               IF(CNVFLG(I).AND.RN(I).gt.0.) THEN
1706                 if(k.ge.kb(i) .and. k.lt.ktop(i)) then
1707                   ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
1708                 endif
1709               endif
1710             ENDDO
1711           ENDDO
1712           DO I = 1, IM
1713             IF(CNVFLG(I).AND.RN(I).gt.0.) THEN
1714                k = ktop(i)-1
1715                dt_mf(i,k) = ud_mf(i,k)
1716             endif
1717           ENDDO
1718           DO K = 1, KM
1719             DO I = 1, IM
1720               IF(CNVFLG(I).AND.RN(I).gt.0.) THEN
1721                 if(k.ge.1 .and. k.le.jmin(i)) then
1722                   dd_mf(i,k) = edto(i) * etad(i,k) * xmb(i) * dt2
1723                 endif
1724               endif
1725             ENDDO
1726           ENDDO
1727     !!
1728           RETURN
1729           END
1730