source: LMDZ4/branches/IPSL-CM4_IPCC_patches/libf/phylmd/convect2.F @ 5451

Last change on this file since 5451 was 524, checked in by lmdzadmin, 21 years ago

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