source: LMDZ4/trunk/libf/phylmd/convect2.F @ 5166

Last change on this file since 5166 was 1403, checked in by Laurent Fairhead, 14 years ago

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

  • IPSLCM5A configuration (comparison between trunk and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent
  • MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

  • IPSLCM5A configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent
  • MH07 configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files

Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

  • pour une config IPSLCM5A (comparaison tronc/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux (à part sur RN et Pb)
    • fichiers histoire égaux
  • pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux
    • fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

  • pour une config IPSLCM5A:
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc
  • pour une config MH07
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc

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