source: LMDZ.3.3/branches/rel-LF/libf/phylmd/convect2.F @ 456

Last change on this file since 456 was 268, checked in by lmdzadmin, 23 years ago

Initial version in branch
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 44.9 KB
Line 
1      subroutine convect2(ncum,idcum,len,nd,ndp1,nl,minorig,
2     &                 nk1,icb1,
3     &                 t1,q1,qs1,u1,v1,gz1,tv1,tp1,tvp1,clw1,h1,
4     &                 lv1,cpn1,p1,ph1,ft1,fq1,fu1,fv1,
5     &                 tnk1,qnk1,gznk1,plcl1,
6     &                 precip1,cbmf1,iflag1,
7     &                 delt,cpd,cpv,cl,rv,rd,lv0,g,
8     &                 sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp,
9     &                 alpha,entp,coeffs,coeffr,omtrain,cu,Ma)
10C.............................START PROLOGUE............................
11C
12C SCCS IDENTIFICATION:  @(#)convect2.f  1.2 05/18/00
13C                       22:06:22 /h/cm/library/nogaps4/src/sub/fcst/convect2.f_v
14C
15C CONFIGURATION IDENTIFICATION:  None
16C
17C MODULE NAME:  convect2
18C
19C DESCRIPTION:
20C
21C convect1     The Emanuel Cumulus Convection Scheme - compute tendencies
22C
23C CONTRACT NUMBER AND TITLE:  None
24C
25C REFERENCES: Programmers  K. Emanuel (MIT), Timothy F. Hogan, M. Peng (NRL)
26C
27C CLASSIFICATION:  Unclassified
28C
29C RESTRICTIONS: None
30C
31C COMPILER DEPENDENCIES: FORTRAN 77, FORTRAN 90
32C
33C COMPILE OPTIONS: Fortran 77: -Zu -Wf"-ei -o aggress"
34C                  Fortran 90: -O vector3,scalar3,task1,aggress,overindex  -ei -r 2
35C
36C LIBRARIES OF RESIDENCE: /a/ops/lib/libfcst159.a
37C
38C USAGE: call convect2(ncum,idcum,len,nd,nl,minorig,
39C    &                 nk1,icb1,
40C    &                 t1,q1,qs1,u1,v1,gz1,tv1,tp1,tvp1,clw1,h1,
41C    &                 lv1,cpn1,p1,ph1,ft1,fq1,fu1,fv1,
42C    &                 tnk1,qnk1,gznk1,plcl1,
43C    &                 precip1,cbmf1,iflag1,
44C    &                 delt,cpd,cpv,cl,rv,rd,lv0,g,
45C    &                 sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp,
46C    &                 alpha,entp,coeffs,coeffr,omtrain,cu)
47C
48C PARAMETERS:
49C      Name            Type         Usage            Description
50C   ----------      ----------     -------  ----------------------------
51C
52C      ncum          Integer        Input        number of cumulus points
53C      idcum         Integer        Input        index of cumulus point
54C      len           Integer        Input        first dimension
55C      nd            Integer        Input        total vertical dimension
56C      ndp1          Integer        Input        nd + 1
57C      nl            Integer        Input        vertical dimension for convection
58C      minorig       Integer        Input        First level where convection is allow to begin
59C      nk1           Integer        Input        First level of convection
60C      ncb1          Integer        Input        Level of free convection
61C      t1            Real           Input        temperature
62C      q1            Real           Input        specific hum
63C      qs1           Real           Input        sat specific hum
64C      u1            Real           Input        u-wind
65C      v1            Real           Input        v-wind
66C      gz1           Real           Inout        geop
67C      tv1           Real           Input        virtual temp
68C      tp1           Real           Input
69C      clw1          Real           Inout        cloud liquid water
70C      h1            Real           Inout
71C      lv1           Real           Inout
72C      cpn1          Real           Inout
73C      p1            Real           Input        full level pressure
74C      ph1           Real           Input        half level pressure
75C      ft1           Real           Output       temp tend
76C      fq1           Real           Output       spec hum tend
77C      fu1           Real           Output       u-wind tend
78C      fv1           Real           Output       v-wind tend
79C      precip1       Real           Output       prec
80C      cbmf1         Real           In/Out       cumulus mass flux
81C      iflag1        Integer        Output       iflag on latitude strip
82C      delt          Real           Input        time step
83C      cpd           Integer        Input        See description below
84C      cpv           Integer        Input        See description below
85C      cl            Integer        Input        See description below
86C      rv            Integer        Input        See description below
87C      rd            Integer        Input        See description below
88C      lv0           Integer        Input        See description below
89C      g             Integer        Input        See description below
90C      sigs          Integer        Input        See description below
91C      sigd          Integer        Input        See description below
92C      elcrit        Integer        Input        See description below
93C      tlcrit        Integer        Input        See description below
94C      omtsnow       Integer        Input        See description below
95C      dtmax         Integer        Input        See description below
96C      damp          Integer        Input        See description below
97C      alpha         Integer        Input        See description below
98C      ent           Integer        Input        See description below
99C      coeffs        Integer        Input        See description below
100C      coeffr        Integer        Input        See description below
101C      omtrain       Integer        Input        See description below
102C      cu            Integer        Input        See description below
103C
104C COMMON BLOCKS:
105C      Block      Name     Type    Usage              Notes
106C     --------  --------   ----    ------   ------------------------
107C
108C FILES: None
109C
110C DATA BASES: None
111C
112C NON-FILE INPUT/OUTPUT: None
113C
114C ERROR CONDITIONS: None
115C
116C ADDITIONAL COMMENTS: None
117C
118C.................MAINTENANCE SECTION................................
119C
120C MODULES CALLED:
121C         Name           Description
122C         zilch        Zero out an array
123C        -------     ----------------------
124C LOCAL VARIABLES AND
125C          STRUCTURES:
126C Name     Type    Description
127C -------  ------  -----------
128C See Comments Below
129C
130C i        Integer loop index
131C k        Integer loop index
132c
133C METHOD:
134C
135C See Emanuel, K. and M. Zivkovic-Rothman, 2000: Development and evaluation of a
136C       convective scheme for use in climate models.
137C
138C FILES: None
139C
140C INCLUDE FILES: None
141C
142C MAKEFILE: /a/ops/met/nogaps/src/sub/fcst/fcst159lib.mak
143C
144C..............................END PROLOGUE.............................
145c
146c
147      implicit none
148c
149#include "dimensions.h"
150#include "dimphy.h"
151c
152      integer kmax2,imax2,kmin2,imin2
153      real ftmax2,ftmin2
154      integer kmax,imax,kmin,imin
155      real ftmax,ftmin
156c
157      integer ncum
158      integer idcum(len)
159      integer len
160      integer nd
161      integer ndp1
162      integer nl
163      integer minorig
164      integer nk1(len)
165      integer icb1(len)
166      real t1(len,nd)
167      real q1(len,nd)
168      real qs1(len,nd)
169      real u1(len,nd)
170      real v1(len,nd)
171      real gz1(len,nd)
172      real tv1(len,nd)
173      real tp1(len,nd)
174      real tvp1(len,nd)
175      real clw1(len,nd)
176      real h1(len,nd)
177      real lv1(len,nd)
178      real cpn1(len,nd)
179      real p1(len,nd)
180      real ph1(len,ndp1)
181      real ft1(len,nd)
182      real fq1(len,nd)
183      real fu1(len,nd)
184      real fv1(len,nd)
185      real tnk1(len)
186      real qnk1(len)
187      real gznk1(len)
188      real precip1(len)
189      real cbmf1(len)
190      real plcl1(len)
191      integer iflag1(len)
192      real delt
193      real cpd
194      real cpv
195      real cl
196      real rv
197      real rd
198      real lv0
199      real g
200      real sigs    ! SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE
201      real sigd    ! SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT
202      real elcrit  ! ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm)
203      real tlcrit  ! TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-
204c                     CONVERSION THRESHOLD IS ASSUMED TO BE ZERO
205      real omtsnow ! OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW
206      real dtmax   ! DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION
207c                    A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC.
208      real damp
209      real alpha
210      real entp    ! ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT FORMULATION
211      real coeffs  ! COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION OF SNOW
212      real coeffr  ! COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION OF RAIN
213      real omtrain ! OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN
214      real cu      ! CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM TRANSPORT
215c
216      real Ma(len,nd)
217c
218C
219C   *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***
220C   ***  TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-        ***
221C   ***       CONVERSION THRESHOLD IS ASSUMED TO BE ZERO             ***
222C   ***     (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY            ***
223C   ***               BETWEEN 0 C AND TLCRIT)                        ***
224C   ***   ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT       ***
225C   ***                       FORMULATION                            ***
226C   ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
227C   ***  SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE       ***
228C   ***                        OF CLOUD                              ***
229C   ***        OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN       ***
230C   ***     OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW          ***
231C   ***  COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
232C   ***                          OF RAIN                             ***
233C   ***  COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
234C   ***                          OF SNOW                             ***
235C   ***     CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM      ***
236C   ***                         TRANSPORT                            ***
237C   ***    DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION    ***
238C   ***        A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC      ***
239C   ***    ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF    ***
240C   ***                 APPROACH TO QUASI-EQUILIBRIUM                ***
241C   ***   (THEIR STANDARD VALUES ARE  0.20 AND 0.1, RESPECTIVELY)    ***
242C   ***                   (DAMP MUST BE LESS THAN 1)                 ***
243c
244c Local arrays.
245c
246      real work(ncum)
247      real t(ncum,klev)
248      real q(ncum,klev)
249      real qs(ncum,klev)
250      real u(ncum,klev)
251      real v(ncum,klev)
252      real gz(ncum,klev)
253      real h(ncum,klev)
254      real lv(ncum,klev)
255      real cpn(ncum,klev)
256      real p(ncum,klev)
257      real ph(ncum,klev)
258      real ft(ncum,klev)
259      real fq(ncum,klev)
260      real fu(ncum,klev)
261      real fv(ncum,klev)
262      real precip(ncum)
263      real cbmf(ncum)
264      real plcl(ncum)
265      real tnk(ncum)
266      real qnk(ncum)
267      real gznk(ncum)
268      real tv(ncum,klev)
269      real tp(ncum,klev)
270      real tvp(ncum,klev)
271      real clw(ncum,klev)
272c     real det(ncum,klev)
273      real dph(ncum,klev)
274c     real wd(ncum)
275c     real tprime(ncum)
276c     real qprime(ncum)
277      real ah0(ncum)
278      real ep(ncum,klev)
279      real sigp(ncum,klev)
280      integer nent(ncum,klev)
281      real water(ncum,klev)
282      real evap(ncum,klev)
283      real mp(ncum,klev)
284      real m(ncum,klev)
285      real qti
286      real wt(ncum,klev)
287      real hp(ncum,klev)
288      real lvcp(ncum,klev)
289      real elij(ncum,klev,klev)
290      real ment(ncum,klev,klev)
291      real sij(ncum,klev,klev)
292      real qent(ncum,klev,klev)
293      real uent(ncum,klev,klev)
294      real vent(ncum,klev,klev)
295      real qp(ncum,klev)
296      real up(ncum,klev)
297      real vp(ncum,klev)
298      real cape(ncum)
299      real capem(ncum)
300      real frac(ncum)
301      real dtpbl(ncum)
302      real tvpplcl(ncum)
303      real tvaplcl(ncum)
304      real dtmin(ncum)
305      real w3d(ncum,klev)
306      real am(ncum)
307      real ents(ncum)
308      real uav(ncum)
309      real vav(ncum)
310c
311      integer iflag(ncum)
312      integer nk(ncum)
313      integer icb(ncum)
314      integer inb(ncum)
315      integer inb1(ncum)
316      integer jtt(ncum)
317c
318      integer nn,i,k,n,icbmax,nlp,j
319      integer ij
320      integer nn2,nn3
321      real clmcpv
322      real clmcpd
323      real cpdmcp
324      real cpvmcpd
325      real eps
326      real epsi
327      real epsim1
328      real tg,qg,s,alv,tc,ahg,denom,es,rg,ginv,rowl
329      real delti
330      real tca,elacrit
331      real by,defrac
332c     real byp
333      real byp(ncum)
334      logical lcape(ncum)
335      real dbo
336      real bf2,anum,dei,altem,cwat,stemp
337      real alt,qp1,smid,sjmax,sjmin
338      real delp,delm
339      real awat,coeff,afac,revap,dhdp,fac,qstm,rat
340      real qsm,sigt,b6,c6
341      real dpinv,cpinv
342      real fqold,ftold,fuold,fvold
343      real wdtrain(ncum),xxx
344      real bsum(ncum,klev)
345      real asij(ncum)
346      real smin(ncum)
347      real scrit(ncum)
348c     real amp1,ad
349      real amp1(ncum),ad(ncum)
350      logical lwork(ncum)
351      integer num1,num2
352c
353c     print*,'cpd en entree de convect2 ',cpd
354      nlp=nl+1
355c
356      rowl=1000.0
357      ginv=1.0/g
358      delti=1.0/delt
359c
360c Define some thermodynamic variables.
361c
362      clmcpv=cl-cpv
363      clmcpd=cl-cpd
364      cpdmcp=cpd-cpv
365      cpvmcpd=cpv-cpd
366      eps=rd/rv
367      epsi=1.0/eps
368      epsim1=epsi-1.0
369c
370c Compress the fields.
371c
372      do 110 k=1,nl+1
373       nn=0
374        do 100 i=1,len
375          if(iflag1(i).eq.0)then
376            nn=nn+1
377            t(nn,k)=t1(i,k)
378            q(nn,k)=q1(i,k)
379            qs(nn,k)=qs1(i,k)
380            u(nn,k)=u1(i,k)
381            v(nn,k)=v1(i,k)
382            gz(nn,k)=gz1(i,k)
383            h(nn,k)=h1(i,k)
384            lv(nn,k)=lv1(i,k)
385            cpn(nn,k)=cpn1(i,k)
386            p(nn,k)=p1(i,k)
387            ph(nn,k)=ph1(i,k)
388            tv(nn,k)=tv1(i,k)
389            tp(nn,k)=tp1(i,k)
390            tvp(nn,k)=tvp1(i,k)
391            clw(nn,k)=clw1(i,k)
392          endif
393 100    continue
394c       print*,'100 ncum,nn',ncum,nn
395 110  continue
396      nn=0
397      do 150 i=1,len
398        if(iflag1(i).eq.0)then
399          nn=nn+1
400          cbmf(nn)=cbmf1(i)
401          plcl(nn)=plcl1(i)
402          tnk(nn)=tnk1(i)
403          qnk(nn)=qnk1(i)
404          gznk(nn)=gznk1(i)
405          nk(nn)=nk1(i)
406          icb(nn)=icb1(i)
407          iflag(nn)=iflag1(i)
408        endif
409 150  continue
410c       print*,'150 ncum,nn',ncum,nn
411c
412c Initialize the tendencies, det, wd, tprime, qprime.
413c
414      do 170 k=1,nl
415        do 160 i=1,ncum
416c         det(i,k)=0.0
417          ft(i,k)=0.0
418          fu(i,k)=0.0
419          fv(i,k)=0.0
420          fq(i,k)=0.0
421          dph(i,k)=ph(i,k)-ph(i,k+1)
422          ep(i,k)=0.0
423          sigp(i,k)=sigs
424 160    continue
425 170  continue
426      do 180 i=1,ncum
427c      wd(i)=0.0
428c      tprime(i)=0.0
429c      qprime(i)=0.0
430       precip(i)=0.0
431       ft(i,nl+1)=0.0
432       fu(i,nl+1)=0.0
433       fv(i,nl+1)=0.0
434       fq(i,nl+1)=0.0
435 180  continue
436c
437c Compute icbmax.
438c
439      icbmax=2
440      do 230 i=1,ncum
441        icbmax=max(icbmax,icb(i))
442 230  continue
443c
444c
445!=====================================================================
446! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
447!=====================================================================
448c
449c ---       The procedure is to solve the equation.
450c              cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
451c
452c   ***  Calculate certain parcel quantities, including static energy   ***
453c
454c
455      do 240 i=1,ncum
456        ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i)
457     &         +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i)
458 240  continue
459c
460c
461c    ***  Find lifted parcel quantities above cloud base    ***
462c
463c
464        do 300 k=minorig+1,nl
465          do 290 i=1,ncum
466            if(k.ge.(icb(i)+1))then
467              tg=t(i,k)
468              qg=qs(i,k)
469              alv=lv0-clmcpv*(t(i,k)-273.15)
470c
471c First iteration.
472c
473               s=cpd+alv*alv*qg/(rv*t(i,k)*t(i,k))
474               s=1./s
475               ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
476               tg=tg+s*(ah0(i)-ahg)
477               tg=max(tg,35.0)
478               tc=tg-273.15
479               denom=243.5+tc
480               if(tc.ge.0.0)then
481                es=6.112*exp(17.67*tc/denom)
482               else
483                es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
484               endif
485                qg=eps*es/(p(i,k)-es*(1.-eps))
486c
487c Second iteration.
488c
489               s=cpd+alv*alv*qg/(rv*t(i,k)*t(i,k))
490               s=1./s
491               ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
492               tg=tg+s*(ah0(i)-ahg)
493               tg=max(tg,35.0)
494               tc=tg-273.15
495               denom=243.5+tc
496               if(tc.ge.0.0)then
497                es=6.112*exp(17.67*tc/denom)
498               else
499                es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
500               endif
501                qg=eps*es/(p(i,k)-es*(1.-eps))
502c
503               alv=lv0-clmcpv*(t(i,k)-273.15)
504c      print*,'cpd dans convect2 ',cpd
505c      print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
506c      print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
507               tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)
508     &                  /cpd
509c              if (.not.cpd.gt.1000.) then
510c                  print*,'CPD=',cpd
511c                  stop
512c              endif
513               clw(i,k)=qnk(i)-qg
514               clw(i,k)=max(0.0,clw(i,k))
515               rg=qg/(1.-qnk(i))
516               tvp(i,k)=tp(i,k)*(1.+rg*epsi)
517            endif
518  290     continue
519  300   continue
520c
521!=====================================================================
522! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
523! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
524! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
525!=====================================================================
526c
527      do 320 k=minorig+1,nl
528        do 310 i=1,ncum
529          if(k.ge.(nk(i)+1))then
530            tca=tp(i,k)-273.15
531            if(tca.ge.0.0)then
532              elacrit=elcrit
533            else
534              elacrit=elcrit*(1.0-tca/tlcrit)
535            endif
536            elacrit=max(elacrit,0.0)
537            ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)
538            ep(i,k)=max(ep(i,k),0.0 )
539            ep(i,k)=min(ep(i,k),1.0 )
540            sigp(i,k)=sigs
541          endif
542 310    continue
543 320  continue
544c
545!=====================================================================
546! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
547! --- VIRTUAL TEMPERATURE
548!=====================================================================
549c
550      do 340 k=minorig+1,nl
551        do 330 i=1,ncum
552        if(k.ge.(icb(i)+1))then
553          tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
554c         print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
555c         print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
556        endif
557 330    continue
558 340  continue
559      do 350 i=1,ncum
560       tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
561 350  continue
562c
563c
564c=====================================================================
565c --- NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
566c=====================================================================
567c
568        do 360 i=1,ncum*nlp
569          nent(i,1)=0
570          water(i,1)=0.0
571          evap(i,1)=0.0
572          mp(i,1)=0.0
573          m(i,1)=0.0
574          wt(i,1)=omtsnow
575          hp(i,1)=h(i,1)
576c         if(.not.cpn(i,1).gt.900.) then
577c         print*,'i,lv(i,1),cpn(i,1)'
578c         print*, i,lv(i,1),cpn(i,1)
579c         k=(i-1)/ncum+1
580c         print*,'i,k',mod(i,ncum),k,'  cpn',cpn(mod(i,ncum),k)
581c         stop
582c         endif
583          lvcp(i,1)=lv(i,1)/cpn(i,1)
584 360    continue
585c
586      do 380 i=1,ncum*nlp*nlp
587        elij(i,1,1)=0.0
588        ment(i,1,1)=0.0
589        sij(i,1,1)=0.0
590 380  continue
591c
592      do 400 k=1,nlp
593       do 390 j=1,nlp
594          do 385 i=1,ncum
595            qent(i,k,j)=q(i,j)
596            uent(i,k,j)=u(i,j)
597            vent(i,k,j)=v(i,j)
598 385      continue
599 390    continue
600 400  continue
601c
602      do 420 i=1,ncum
603        qp(i,1)=q(i,1)
604        up(i,1)=u(i,1)
605        vp(i,1)=v(i,1)
606 420  continue
607      do 440 k=2,nlp
608        do 430 i=1,ncum
609          qp(i,k)=q(i,k-1)
610          up(i,k)=u(i,k-1)
611          vp(i,k)=v(i,k-1)
612 430    continue
613 440  continue
614c
615c=====================================================================
616c  --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
617c  --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
618c  --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
619c=====================================================================
620c
621      do 510 i=1,ncum
622        cape(i)=0.0
623        capem(i)=0.0
624        inb(i)=icb(i)+1
625        inb1(i)=inb(i)
626 510  continue
627c
628c Originial Code
629c
630c     do 530 k=minorig+1,nl-1
631c       do 520 i=1,ncum
632c         if(k.ge.(icb(i)+1))then
633c           by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
634c           byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
635c           cape(i)=cape(i)+by
636c           if(by.ge.0.0)inb1(i)=k+1
637c           if(cape(i).gt.0.0)then
638c             inb(i)=k+1
639c             capem(i)=cape(i)
640c           endif
641c         endif
642c520    continue
643c530  continue
644c     do 540 i=1,ncum
645c         byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
646c         cape(i)=capem(i)+byp
647c         defrac=capem(i)-cape(i)
648c         defrac=max(defrac,0.001)
649c         frac(i)=-cape(i)/defrac
650c         frac(i)=min(frac(i),1.0)
651c         frac(i)=max(frac(i),0.0)
652c540   continue
653c
654c K Emanuel fix
655c
656c     call zilch(byp,ncum)
657c     do 530 k=minorig+1,nl-1
658c       do 520 i=1,ncum
659c         if(k.ge.(icb(i)+1))then
660c           by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
661c           cape(i)=cape(i)+by
662c           if(by.ge.0.0)inb1(i)=k+1
663c           if(cape(i).gt.0.0)then
664c             inb(i)=k+1
665c             capem(i)=cape(i)
666c             byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
667c           endif
668c         endif
669c520    continue
670c530  continue
671c     do 540 i=1,ncum
672c         inb(i)=max(inb(i),inb1(i))
673c         cape(i)=capem(i)+byp(i)
674c         defrac=capem(i)-cape(i)
675c         defrac=max(defrac,0.001)
676c         frac(i)=-cape(i)/defrac
677c         frac(i)=min(frac(i),1.0)
678c         frac(i)=max(frac(i),0.0)
679c540   continue
680c
681c J Teixeira fix
682c
683      call zilch(byp,ncum)
684      do 515 i=1,ncum
685        lcape(i)=.true.
686 515  continue
687      do 530 k=minorig+1,nl-1
688        do 520 i=1,ncum
689          if(cape(i).lt.0.0)lcape(i)=.false.
690          if((k.ge.(icb(i)+1)).and.lcape(i))then
691            by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
692            byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
693            cape(i)=cape(i)+by
694            if(by.ge.0.0)inb1(i)=k+1
695            if(cape(i).gt.0.0)then
696              inb(i)=k+1
697              capem(i)=cape(i)
698            endif
699          endif
700 520    continue
701 530  continue
702      do 540 i=1,ncum
703          cape(i)=capem(i)+byp(i)
704          defrac=capem(i)-cape(i)
705          defrac=max(defrac,0.001)
706          frac(i)=-cape(i)/defrac
707          frac(i)=min(frac(i),1.0)
708          frac(i)=max(frac(i),0.0)
709 540  continue
710c
711c=====================================================================
712c ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
713c=====================================================================
714c
715      do 600 k=minorig+1,nl
716        do 590 i=1,ncum
717        if((k.ge.icb(i)).and.(k.le.inb(i)))then
718          hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k)
719        endif
720 590    continue
721 600  continue
722c
723c=====================================================================
724c ---  CALCULATE CLOUD BASE MASS FLUX AND RATES OF MIXING, M(I),
725c --- AT EACH MODEL LEVEL
726c=====================================================================
727c
728c tvpplcl = parcel temperature lifted adiabatically from level
729c           icb-1 to the LCL.
730c tvaplcl = virtual temperature at the LCL.
731c
732      do 610 i=1,ncum
733        dtpbl(i)=0.0
734        tvpplcl(i)=tvp(i,icb(i)-1)
735     &  -rd*tvp(i,icb(i)-1)*(p(i,icb(i)-1)-plcl(i))
736     &  /(cpn(i,icb(i)-1)*p(i,icb(i)-1))
737        tvaplcl(i)=tv(i,icb(i))
738     &  +(tvp(i,icb(i))-tvp(i,icb(i)+1))*(plcl(i)-p(i,icb(i)))
739     &  /(p(i,icb(i))-p(i,icb(i)+1))
740 610  continue
741c
742c-------------------------------------------------------------------
743c --- Interpolate difference between lifted parcel and
744c --- environmental temperatures to lifted condensation level
745c-------------------------------------------------------------------
746c
747c dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
748c
749      do 630 k=minorig,icbmax
750        do 620 i=1,ncum
751        if((k.ge.nk(i)).and.(k.le.(icb(i)-1)))then
752          dtpbl(i)=dtpbl(i)+(tvp(i,k)-tv(i,k))*dph(i,k)
753        endif
754 620    continue
755 630  continue
756      do 640 i=1,ncum
757        dtpbl(i)=dtpbl(i)/(ph(i,nk(i))-ph(i,icb(i)))
758        dtmin(i)=tvpplcl(i)-tvaplcl(i)+dtmax+dtpbl(i)
759 640  continue
760c
761c-------------------------------------------------------------------
762c --- Adjust cloud base mass flux
763c-------------------------------------------------------------------
764c
765      do 650 i=1,ncum
766       work(i)=cbmf(i)
767       cbmf(i)=max(0.0,(1.0-damp)*cbmf(i)+0.1*alpha*dtmin(i))
768       if((work(i).eq.0.0).and.(cbmf(i).eq.0.0))then
769         iflag(i)=3
770       endif
771 650  continue
772c
773c-------------------------------------------------------------------
774c --- Calculate rates of mixing,  m(i)
775c-------------------------------------------------------------------
776c
777      call zilch(work,ncum)
778c
779      do 670 j=minorig+1,nl
780        do 660 i=1,ncum
781          if((j.ge.(icb(i)+1)).and.(j.le.inb(i)))then
782             k=min(j,inb1(i))
783             dbo=abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1))
784     &       +entp*0.04*(ph(i,k)-ph(i,k+1))
785             work(i)=work(i)+dbo
786             m(i,j)=cbmf(i)*dbo
787          endif
788 660    continue
789 670  continue
790      do 690 k=minorig+1,nl
791        do 680 i=1,ncum
792          if((k.ge.(icb(i)+1)).and.(k.le.inb(i)))then
793            m(i,k)=m(i,k)/work(i)
794          endif
795 680    continue
796 690  continue
797c
798c
799c=====================================================================
800c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
801c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
802c --- FRACTION (sij)
803c=====================================================================
804c
805c
806       do 750 i=minorig+1,nl
807         do 710 j=minorig+1,nl
808           do 700 ij=1,ncum
809             if((i.ge.(icb(ij)+1)).and.(j.ge.icb(ij))
810     &         .and.(i.le.inb(ij)).and.(j.le.inb(ij)))then
811               qti=qnk(ij)-ep(ij,i)*clw(ij,i)
812               bf2=1.+lv(ij,j)*lv(ij,j)*qs(ij,j)
813     &         /(rv*t(ij,j)*t(ij,j)*cpd)
814               anum=h(ij,j)-hp(ij,i)+(cpv-cpd)*t(ij,j)*(qti-q(ij,j))
815               denom=h(ij,i)-hp(ij,i)+(cpd-cpv)*(q(ij,i)-qti)*t(ij,j)
816               dei=denom
817               if(abs(dei).lt.0.01)dei=0.01
818               sij(ij,i,j)=anum/dei
819               sij(ij,i,i)=1.0
820               altem=sij(ij,i,j)*q(ij,i)+(1.-sij(ij,i,j))*qti-qs(ij,j)
821               altem=altem/bf2
822               cwat=clw(ij,j)*(1.-ep(ij,j))
823               stemp=sij(ij,i,j)
824               if((stemp.lt.0.0.or.stemp.gt.1.0.or.
825     1           altem.gt.cwat).and.j.gt.i)then
826                 anum=anum-lv(ij,j)*(qti-qs(ij,j)-cwat*bf2)
827                 denom=denom+lv(ij,j)*(q(ij,i)-qti)
828                 if(abs(denom).lt.0.01)denom=0.01
829                 sij(ij,i,j)=anum/denom
830                 altem=sij(ij,i,j)*q(ij,i)+(1.-sij(ij,i,j))*qti-qs(ij,j)
831                 altem=altem-(bf2-1.)*cwat
832               endif
833               if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then
834                 qent(ij,i,j)=sij(ij,i,j)*q(ij,i)
835     &                        +(1.-sij(ij,i,j))*qti
836                 uent(ij,i,j)=sij(ij,i,j)*u(ij,i)
837     &                        +(1.-sij(ij,i,j))*u(ij,nk(ij))
838                 vent(ij,i,j)=sij(ij,i,j)*v(ij,i)
839     &                        +(1.-sij(ij,i,j))*v(ij,nk(ij))
840                 elij(ij,i,j)=altem
841                 elij(ij,i,j)=max(0.0,elij(ij,i,j))
842                 ment(ij,i,j)=m(ij,i)/(1.-sij(ij,i,j))
843                 nent(ij,i)=nent(ij,i)+1
844               endif
845             sij(ij,i,j)=max(0.0,sij(ij,i,j))
846             sij(ij,i,j)=min(1.0,sij(ij,i,j))
847             endif
848  700      continue
849  710    continue
850c
851c   ***   If no air can entrain at level i assume that updraft detrains  ***
852c   ***   at that level and calculate detrained air flux and properties  ***
853c
854           do 740 ij=1,ncum
855             if((i.ge.(icb(ij)+1)).and.(i.le.inb(ij))
856     &       .and.(nent(ij,i).eq.0))then
857               ment(ij,i,i)=m(ij,i)
858               qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
859               uent(ij,i,i)=u(ij,nk(ij))
860               vent(ij,i,i)=v(ij,nk(ij))
861               elij(ij,i,i)=clw(ij,i)
862               sij(ij,i,i)=1.0
863             endif
864 740       continue
865 750   continue
866c
867      do 770 i=1,ncum
868        sij(i,inb(i),inb(i))=1.0
869 770  continue
870c
871c=====================================================================
872c   ---  NORMALIZE ENTRAINED AIR MASS FLUXES
873c   ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
874c=====================================================================
875c
876c
877       call zilch(bsum,ncum*nlp)
878       do 780 ij=1,ncum
879         lwork(ij)=.false.
880 780   continue
881       do 789 i=minorig+1,nl
882c
883         num1=0
884         do 779 ij=1,ncum
885           if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))num1=num1+1
886 779     continue
887         if(num1.le.0)go to 789
888c
889           do 781 ij=1,ncum
890             if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))then
891                lwork(ij)=(nent(ij,i).ne.0)
892                qp1=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
893                anum=h(ij,i)-hp(ij,i)-lv(ij,i)*(qp1-qs(ij,i))
894                denom=h(ij,i)-hp(ij,i)+lv(ij,i)*(q(ij,i)-qp1)
895                if(abs(denom).lt.0.01)denom=0.01
896                scrit(ij)=anum/denom
897                alt=qp1-qs(ij,i)+scrit(ij)*(q(ij,i)-qp1)
898                if(scrit(ij).lt.0.0.or.alt.lt.0.0)scrit(ij)=1.0
899                asij(ij)=0.0
900                smin(ij)=1.0
901             endif
902 781       continue
903         do 783 j=minorig,nl
904c
905         num2=0
906         do 778 ij=1,ncum
907             if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
908     &       .and.(j.ge.icb(ij)).and.(j.le.inb(ij))
909     &       .and.lwork(ij))num2=num2+1
910 778     continue
911         if(num2.le.0)go to 783
912c
913           do 782 ij=1,ncum
914             if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
915     &       .and.(j.ge.icb(ij)).and.(j.le.inb(ij)).and.lwork(ij))then
916                  if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then
917                    if(j.gt.i)then
918                      smid=min(sij(ij,i,j),scrit(ij))
919                      sjmax=smid
920                      sjmin=smid
921                        if(smid.lt.smin(ij)
922     &                  .and.sij(ij,i,j+1).lt.smid)then
923                          smin(ij)=smid
924                          sjmax=min(sij(ij,i,j+1),sij(ij,i,j),scrit(ij))
925                          sjmin=max(sij(ij,i,j-1),sij(ij,i,j))
926                          sjmin=min(sjmin,scrit(ij))
927                        endif
928                    else
929                      sjmax=max(sij(ij,i,j+1),scrit(ij))
930                      smid=max(sij(ij,i,j),scrit(ij))
931                      sjmin=0.0
932                      if(j.gt.1)sjmin=sij(ij,i,j-1)
933                      sjmin=max(sjmin,scrit(ij))
934                    endif
935                    delp=abs(sjmax-smid)
936                    delm=abs(sjmin-smid)
937                    asij(ij)=asij(ij)+(delp+delm)
938     &                           *(ph(ij,j)-ph(ij,j+1))
939                    ment(ij,i,j)=ment(ij,i,j)*(delp+delm)
940     &                           *(ph(ij,j)-ph(ij,j+1))
941                  endif
942              endif
943  782    continue
944  783    continue
945            do 784 ij=1,ncum
946            if((i.ge.icb(ij)+1).and.(i.le.inb(ij)).and.lwork(ij))then
947               asij(ij)=max(1.0e-21,asij(ij))
948               asij(ij)=1.0/asij(ij)
949               bsum(ij,i)=0.0
950            endif
951 784        continue
952            do 786 j=minorig,nl+1
953              do 785 ij=1,ncum
954                if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
955     &          .and.(j.ge.icb(ij)).and.(j.le.inb(ij))
956     &          .and.lwork(ij))then
957                   ment(ij,i,j)=ment(ij,i,j)*asij(ij)
958                   bsum(ij,i)=bsum(ij,i)+ment(ij,i,j)
959                endif
960 785     continue
961 786     continue
962             do 787 ij=1,ncum
963               if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
964     &         .and.(bsum(ij,i).lt.1.0e-18).and.lwork(ij))then
965                 nent(ij,i)=0
966                 ment(ij,i,i)=m(ij,i)
967                 qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
968                 uent(ij,i,i)=u(ij,nk(ij))
969                 vent(ij,i,i)=v(ij,nk(ij))
970                 elij(ij,i,i)=clw(ij,i)
971                 sij(ij,i,i)=1.0
972               endif
973  787        continue
974  789  continue
975c
976c=====================================================================
977c --- PRECIPITATING DOWNDRAFT CALCULATION
978c=====================================================================
979c
980c   ***  Check whether ep(inb)=0, if so, skip precipitating    ***
981c   ***             downdraft calculation                      ***
982c
983c
984c   ***  Integrate liquid water equation to find condensed water   ***
985c   ***                and condensed water flux                    ***
986c
987c
988      do 890 i=1,ncum
989        jtt(i)=2
990        if(ep(i,inb(i)).le.0.0001)iflag(i)=2
991        if(iflag(i).eq.0)then
992          lwork(i)=.true.
993        else
994          lwork(i)=.false.
995        endif
996 890  continue
997c
998c    ***                    Begin downdraft loop                    ***
999c
1000c
1001        call zilch(wdtrain,ncum)
1002        do 899 i=nl+1,1,-1
1003c
1004          num1=0
1005          do 879 ij=1,ncum
1006            if((i.le.inb(ij)).and.lwork(ij))num1=num1+1
1007 879      continue
1008          if(num1.le.0)go to 899
1009c
1010c
1011c    ***        Calculate detrained precipitation             ***
1012c
1013          do 891 ij=1,ncum
1014            if((i.le.inb(ij)).and.(lwork(ij)))then
1015            wdtrain(ij)=g*ep(ij,i)*m(ij,i)*clw(ij,i)
1016            endif
1017 891      continue
1018c
1019          if(i.gt.1)then
1020            do 893 j=1,i-1
1021              do 892 ij=1,ncum
1022                if((i.le.inb(ij)).and.(lwork(ij)))then
1023                  awat=elij(ij,j,i)-(1.-ep(ij,i))*clw(ij,i)
1024                  awat=max(0.0,awat)
1025                  wdtrain(ij)=wdtrain(ij)+g*awat*ment(ij,j,i)
1026                endif
1027 892          continue
1028 893      continue
1029          endif
1030c
1031c    ***    Find rain water and evaporation using provisional   ***
1032c    ***              estimates of qp(i)and qp(i-1)             ***
1033c
1034c
1035c  ***  Value of terminal velocity and coeffecient of evaporation for snow   ***
1036c
1037          do 894 ij=1,ncum
1038            if((i.le.inb(ij)).and.(lwork(ij)))then
1039            coeff=coeffs
1040            wt(ij,i)=omtsnow
1041c
1042c  ***  Value of terminal velocity and coeffecient of evaporation for rain   ***
1043c
1044            if(t(ij,i).gt.273.0)then
1045              coeff=coeffr
1046              wt(ij,i)=omtrain
1047            endif
1048            qsm=0.5*(q(ij,i)+qp(ij,i+1))
1049            afac=coeff*ph(ij,i)*(qs(ij,i)-qsm)
1050     &       /(1.0e4+2.0e3*ph(ij,i)*qs(ij,i))
1051            afac=max(afac,0.0)
1052            sigt=sigp(ij,i)
1053            sigt=max(0.0,sigt)
1054            sigt=min(1.0,sigt)
1055            b6=100.*(ph(ij,i)-ph(ij,i+1))*sigt*afac/wt(ij,i)
1056            c6=(water(ij,i+1)*wt(ij,i+1)+wdtrain(ij)/sigd)/wt(ij,i)
1057            revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
1058            evap(ij,i)=sigt*afac*revap
1059            water(ij,i)=revap*revap
1060c
1061c    ***  Calculate precipitating downdraft mass flux under     ***
1062c    ***              hydrostatic approximation                 ***
1063c
1064            if(i.gt.1)then
1065              dhdp=(h(ij,i)-h(ij,i-1))/(p(ij,i-1)-p(ij,i))
1066              dhdp=max(dhdp,10.0)
1067              mp(ij,i)=100.*ginv*lv(ij,i)*sigd*evap(ij,i)/dhdp
1068              mp(ij,i)=max(mp(ij,i),0.0)
1069c
1070c   ***   Add small amount of inertia to downdraft              ***
1071c
1072              fac=20.0/(ph(ij,i-1)-ph(ij,i))
1073              mp(ij,i)=(fac*mp(ij,i+1)+mp(ij,i))/(1.+fac)
1074c
1075c    ***      Force mp to decrease linearly to zero                 ***
1076c    ***      between about 950 mb and the surface                  ***
1077c
1078              if(p(ij,i).gt.(0.949*p(ij,1)))then
1079                 jtt(ij)=max(jtt(ij),i)
1080                 mp(ij,i)=mp(ij,jtt(ij))*(p(ij,1)-p(ij,i))
1081     &           /(p(ij,1)-p(ij,jtt(ij)))
1082              endif
1083            endif
1084c
1085c    ***       Find mixing ratio of precipitating downdraft     ***
1086c
1087            if(i.ne.inb(ij))then
1088              if(i.eq.1)then
1089                qstm=qs(ij,1)
1090              else
1091                qstm=qs(ij,i-1)
1092              endif
1093              if(mp(ij,i).gt.mp(ij,i+1))then
1094                 rat=mp(ij,i+1)/mp(ij,i)
1095                 qp(ij,i)=qp(ij,i+1)*rat+q(ij,i)*(1.0-rat)+100.*ginv*
1096     &             sigd*(ph(ij,i)-ph(ij,i+1))*(evap(ij,i)/mp(ij,i))
1097                 up(ij,i)=up(ij,i+1)*rat+u(ij,i)*(1.-rat)
1098                 vp(ij,i)=vp(ij,i+1)*rat+v(ij,i)*(1.-rat)
1099               else
1100                 if(mp(ij,i+1).gt.0.0)then
1101                   qp(ij,i)=(gz(ij,i+1)-gz(ij,i)
1102     &               +qp(ij,i+1)*(lv(ij,i+1)+t(ij,i+1)
1103     &               *(cl-cpd))+cpd*(t(ij,i+1)-t(ij,i)))
1104     &               /(lv(ij,i)+t(ij,i)*(cl-cpd))
1105                   up(ij,i)=up(ij,i+1)
1106                   vp(ij,i)=vp(ij,i+1)
1107                 endif
1108              endif
1109              qp(ij,i)=min(qp(ij,i),qstm)
1110              qp(ij,i)=max(qp(ij,i),0.0)
1111            endif
1112            endif
1113 894      continue
1114 899    continue
1115c
1116c   ***  Calculate surface precipitation in mm/day     ***
1117c
1118        do 1190 i=1,ncum
1119          if(iflag(i).le.1)then
1120cc            precip(i)=precip(i)+wt(i,1)*sigd*water(i,1)*3600.*24000.
1121cc     &                /(rowl*g)
1122cc            precip(i)=precip(i)*delt/86400.
1123            precip(i) = wt(i,1)*sigd*water(i,1)*86400/g
1124          endif
1125 1190   continue
1126c
1127c
1128c   ***  Calculate downdraft velocity scale and surface temperature and  ***
1129c   ***                    water vapor fluctuations                      ***
1130c
1131c     wd=beta*abs(mp(icb))*0.01*rd*t(icb)/(sigd*p(icb))
1132c     qprime=0.5*(qp(1)-q(1))
1133c     tprime=lv0*qprime/cpd
1134c
1135c   ***  Calculate tendencies of lowest level potential temperature  ***
1136c   ***                      and mixing ratio                        ***
1137c
1138        do 1200 i=1,ncum
1139          work(i)=0.01/(ph(i,1)-ph(i,2))
1140          am(i)=0.0
1141 1200   continue
1142        do 1220 k=2,nl
1143          do 1210 i=1,ncum
1144            if((nk(i).eq.1).and.(k.le.inb(i)).and.(nk(i).eq.1))then
1145              am(i)=am(i)+m(i,k)
1146            endif
1147 1210     continue
1148 1220   continue
1149        do 1240 i=1,ncum
1150          if((g*work(i)*am(i)).ge.delti)iflag(i)=1
1151          ft(i,1)=ft(i,1)+g*work(i)*am(i)*(t(i,2)-t(i,1)
1152     &    +(gz(i,2)-gz(i,1))/cpn(i,1))
1153          ft(i,1)=ft(i,1)-lvcp(i,1)*sigd*evap(i,1)
1154          ft(i,1)=ft(i,1)+sigd*wt(i,2)*(cl-cpd)*water(i,2)*(t(i,2)
1155     &     -t(i,1))*work(i)/cpn(i,1)
1156          fq(i,1)=fq(i,1)+g*mp(i,2)*(qp(i,2)-q(i,1))*
1157     &    work(i)+sigd*evap(i,1)
1158          fq(i,1)=fq(i,1)+g*am(i)*(q(i,2)-q(i,1))*work(i)
1159          fu(i,1)=fu(i,1)+g*work(i)*(mp(i,2)*(up(i,2)-u(i,1))
1160     &    +am(i)*(u(i,2)-u(i,1)))
1161          fv(i,1)=fv(i,1)+g*work(i)*(mp(i,2)*(vp(i,2)-v(i,1))
1162     &    +am(i)*(v(i,2)-v(i,1)))
1163 1240   continue
1164        do 1260 j=2,nl
1165           do 1250 i=1,ncum
1166             if(j.le.inb(i))then
1167               fq(i,1)=fq(i,1)
1168     &                 +g*work(i)*ment(i,j,1)*(qent(i,j,1)-q(i,1))
1169               fu(i,1)=fu(i,1)
1170     &                 +g*work(i)*ment(i,j,1)*(uent(i,j,1)-u(i,1))
1171               fv(i,1)=fv(i,1)
1172     &                 +g*work(i)*ment(i,j,1)*(vent(i,j,1)-v(i,1))
1173             endif
1174 1250      continue
1175 1260   continue
1176c
1177c   ***  Calculate tendencies of potential temperature and mixing ratio  ***
1178c   ***               at levels above the lowest level                   ***
1179c
1180c   ***  First find the net saturated updraft and downdraft mass fluxes  ***
1181c   ***                      through each level                          ***
1182c
1183        do 1500 i=2,nl+1
1184c
1185          num1=0
1186          do 1265 ij=1,ncum
1187            if(i.le.inb(ij))num1=num1+1
1188 1265     continue
1189          if(num1.le.0)go to 1500
1190c
1191          call zilch(amp1,ncum)
1192          call zilch(ad,ncum)
1193c
1194          do 1280 k=i+1,nl+1
1195            do 1270 ij=1,ncum
1196              if((i.ge.nk(ij)).and.(i.le.inb(ij))
1197     &            .and.(k.le.(inb(ij)+1)))then
1198                amp1(ij)=amp1(ij)+m(ij,k)
1199              endif
1200 1270         continue
1201 1280     continue
1202c
1203          do 1310 k=1,i
1204            do 1300 j=i+1,nl+1
1205               do 1290 ij=1,ncum
1206                 if((j.le.(inb(ij)+1)).and.(i.le.inb(ij)))then
1207                   amp1(ij)=amp1(ij)+ment(ij,k,j)
1208                 endif
1209 1290          continue
1210 1300       continue
1211 1310     continue
1212          do 1340 k=1,i-1
1213            do 1330 j=i,nl+1
1214              do 1320 ij=1,ncum
1215                if((i.le.inb(ij)).and.(j.le.inb(ij)))then
1216                   ad(ij)=ad(ij)+ment(ij,j,k)
1217                endif
1218 1320         continue
1219 1330       continue
1220 1340     continue
1221c
1222          do 1350 ij=1,ncum
1223          if(i.le.inb(ij))then
1224            dpinv=0.01/(ph(ij,i)-ph(ij,i+1))
1225            cpinv=1.0/cpn(ij,i)
1226c
1227            ft(ij,i)=ft(ij,i)
1228     &       +g*dpinv*(amp1(ij)*(t(ij,i+1)-t(ij,i)
1229     &       +(gz(ij,i+1)-gz(ij,i))*cpinv)
1230     &       -ad(ij)*(t(ij,i)-t(ij,i-1)+(gz(ij,i)-gz(ij,i-1))*cpinv))
1231     &       -sigd*lvcp(ij,i)*evap(ij,i)
1232            ft(ij,i)=ft(ij,i)+g*dpinv*ment(ij,i,i)*(hp(ij,i)-h(ij,i)+
1233     &        t(ij,i)*(cpv-cpd)*(q(ij,i)-qent(ij,i,i)))*cpinv
1234            ft(ij,i)=ft(ij,i)+sigd*wt(ij,i+1)*(cl-cpd)*water(ij,i+1)*
1235     &        (t(ij,i+1)-t(ij,i))*dpinv*cpinv
1236            fq(ij,i)=fq(ij,i)+g*dpinv*(amp1(ij)*(q(ij,i+1)-q(ij,i))-
1237     &        ad(ij)*(q(ij,i)-q(ij,i-1)))
1238            fu(ij,i)=fu(ij,i)+g*dpinv*(amp1(ij)*(u(ij,i+1)-u(ij,i))-
1239     &        ad(ij)*(u(ij,i)-u(ij,i-1)))
1240            fv(ij,i)=fv(ij,i)+g*dpinv*(amp1(ij)*(v(ij,i+1)-v(ij,i))-
1241     &        ad(ij)*(v(ij,i)-v(ij,i-1)))
1242         endif
1243 1350    continue
1244         do 1370 k=1,i-1
1245           do 1360 ij=1,ncum
1246             if(i.le.inb(ij))then
1247               awat=elij(ij,k,i)-(1.-ep(ij,i))*clw(ij,i)
1248               awat=max(awat,0.0)
1249               fq(ij,i)=fq(ij,i)
1250     &         +g*dpinv*ment(ij,k,i)*(qent(ij,k,i)-awat-q(ij,i))
1251               fu(ij,i)=fu(ij,i)
1252     &         +g*dpinv*ment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
1253               fv(ij,i)=fv(ij,i)
1254     &         +g*dpinv*ment(ij,k,i)*(vent(ij,k,i)-v(ij,i))
1255             endif
1256 1360      continue
1257 1370    continue
1258         do 1390 k=i,nl+1
1259           do 1380 ij=1,ncum
1260             if((i.le.inb(ij)).and.(k.le.inb(ij)))then
1261               fq(ij,i)=fq(ij,i)
1262     &                  +g*dpinv*ment(ij,k,i)*(qent(ij,k,i)-q(ij,i))
1263               fu(ij,i)=fu(ij,i)
1264     &                  +g*dpinv*ment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
1265               fv(ij,i)=fv(ij,i)
1266     &                  +g*dpinv*ment(ij,k,i)*(vent(ij,k,i)-v(ij,i))
1267             endif
1268 1380      continue
1269 1390    continue
1270          do 1400 ij=1,ncum
1271           if(i.le.inb(ij))then
1272             fq(ij,i)=fq(ij,i)
1273     &                +sigd*evap(ij,i)+g*(mp(ij,i+1)*
1274     &                (qp(ij,i+1)-q(ij,i))
1275     &                -mp(ij,i)*(qp(ij,i)-q(ij,i-1)))*dpinv
1276             fu(ij,i)=fu(ij,i)
1277     &                +g*(mp(ij,i+1)*(up(ij,i+1)-u(ij,i))-mp(ij,i)*
1278     &                (up(ij,i)-u(ij,i-1)))*dpinv
1279             fv(ij,i)=fv(ij,i)
1280     &               +g*(mp(ij,i+1)*(vp(ij,i+1)-v(ij,i))-mp(ij,i)*
1281     &               (vp(ij,i)-v(ij,i-1)))*dpinv
1282           endif
1283 1400     continue
1284 1500   continue
1285c
1286c   *** Adjust tendencies at top of convection layer to reflect  ***
1287c   ***       actual position of the level zero cape             ***
1288c
1289        do 503 ij=1,ncum
1290        fqold=fq(ij,inb(ij))
1291        fq(ij,inb(ij))=fq(ij,inb(ij))*(1.-frac(ij))
1292        fq(ij,inb(ij)-1)=fq(ij,inb(ij)-1)
1293     &   +frac(ij)*fqold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
1294     1   (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*lv(ij,inb(ij))
1295     &   /lv(ij,inb(ij)-1)
1296        ftold=ft(ij,inb(ij))
1297        ft(ij,inb(ij))=ft(ij,inb(ij))*(1.-frac(ij))
1298        ft(ij,inb(ij)-1)=ft(ij,inb(ij)-1)
1299     &   +frac(ij)*ftold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
1300     1   (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*cpn(ij,inb(ij))
1301     &   /cpn(ij,inb(ij)-1)
1302        fuold=fu(ij,inb(ij))
1303        fu(ij,inb(ij))=fu(ij,inb(ij))*(1.-frac(ij))
1304        fu(ij,inb(ij)-1)=fu(ij,inb(ij)-1)
1305     &   +frac(ij)*fuold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
1306     1   (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
1307        fvold=fv(ij,inb(ij))
1308        fv(ij,inb(ij))=fv(ij,inb(ij))*(1.-frac(ij))
1309        fv(ij,inb(ij)-1)=fv(ij,inb(ij)-1)
1310     &  +frac(ij)*fvold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
1311     1   (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
1312 503    continue
1313c
1314c   ***   Very slightly adjust tendencies to force exact   ***
1315c   ***     enthalpy, momentum and tracer conservation     ***
1316c
1317        do 682 ij=1,ncum
1318        ents(ij)=0.0
1319        uav(ij)=0.0
1320        vav(ij)=0.0
1321        do 681 i=1,inb(ij)
1322         ents(ij)=ents(ij)
1323     &  +(cpn(ij,i)*ft(ij,i)+lv(ij,i)*fq(ij,i))*(ph(ij,i)-ph(ij,i+1))   
1324         uav(ij)=uav(ij)+fu(ij,i)*(ph(ij,i)-ph(ij,i+1))
1325         vav(ij)=vav(ij)+fv(ij,i)*(ph(ij,i)-ph(ij,i+1))
1326  681   continue
1327  682   continue
1328        do 683 ij=1,ncum
1329        ents(ij)=ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1330        uav(ij)=uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1331        vav(ij)=vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1332 683    continue
1333        do 642 ij=1,ncum
1334        do 641 i=1,inb(ij)
1335         ft(ij,i)=ft(ij,i)-ents(ij)/cpn(ij,i)
1336         fu(ij,i)=(1.-cu)*(fu(ij,i)-uav(ij))
1337         fv(ij,i)=(1.-cu)*(fv(ij,i)-vav(ij))
1338  641   continue
1339 642    continue
1340c
1341        do 1810 k=1,nl+1
1342          do 1800 i=1,ncum
1343            if((q(i,k)+delt*fq(i,k)).lt.0.0)iflag(i)=10
1344 1800     continue
1345 1810   continue
1346c
1347c
1348        do 1900 i=1,ncum
1349          if(iflag(i).gt.2)then
1350          precip(i)=0.0
1351          cbmf(i)=0.0
1352          endif
1353 1900   continue
1354        do 1920 k=1,nl
1355         do 1910 i=1,ncum
1356           if(iflag(i).gt.2)then
1357             ft(i,k)=0.0
1358             fq(i,k)=0.0
1359             fu(i,k)=0.0
1360             fv(i,k)=0.0
1361           endif
1362 1910    continue
1363 1920   continue
1364        do 2000 i=1,ncum
1365         precip1(idcum(i))=precip(i)
1366         cbmf1(idcum(i))=cbmf(i)
1367         iflag1(idcum(i))=iflag(i)
1368 2000   continue
1369        do 2020 k=1,nl
1370          do 2010 i=1,ncum
1371            ft1(idcum(i),k)=ft(i,k)
1372            fq1(idcum(i),k)=fq(i,k)
1373            fu1(idcum(i),k)=fu(i,k)
1374            fv1(idcum(i),k)=fv(i,k)
1375 2010     continue
1376 2020   continue
1377c
1378      DO k=1,nd
1379        DO i=1,len
1380         Ma(i,k) = 0.
1381        ENDDO
1382      ENDDO
1383      DO k=nl,1,-1
1384        DO i=1,ncum
1385          Ma(i,k) = Ma(i,k+1)+m(i,k)
1386        ENDDO
1387      ENDDO
1388c
1389        return
1390        end
1391
Note: See TracBrowser for help on using the repository browser.