My Project
 All Classes Files Functions Variables Macros
convect2.F
Go to the documentation of this file.
1 !
2 ! $Id: convect2.F 1403 2010-07-01 09:02:53Z fairhead $
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)
13 C.............................START PROLOGUE............................
14 C
15 C SCCS IDENTIFICATION: @(#)convect2.f 1.2 05/18/00
16 C 22:06:22 /h/cm/library/nogaps4/src/sub/fcst/convect2.f_v
17 C
18 C CONFIGURATION IDENTIFICATION: None
19 C
20 C MODULE NAME: convect2
21 C
22 C DESCRIPTION:
23 C
24 C convect1 The Emanuel Cumulus Convection Scheme - compute tendencies
25 C
26 C CONTRACT NUMBER AND TITLE: None
27 C
28 C REFERENCES: Programmers K. Emanuel (MIT), Timothy F. Hogan, M. Peng (NRL)
29 C
30 C CLASSIFICATION: Unclassified
31 C
32 C RESTRICTIONS: None
33 C
34 C COMPILER DEPENDENCIES: FORTRAN 77, FORTRAN 90
35 C
36 C COMPILE OPTIONS: Fortran 77: -Zu -Wf"-ei -o aggress"
37 C Fortran 90: -O vector3,scalar3,task1,aggress,overindex -ei -r 2
38 C
39 C LIBRARIES OF RESIDENCE: /a/ops/lib/libfcst159.a
40 C
41 C USAGE: call convect2(ncum,idcum,len,nd,nl,minorig,
42 C & nk1,icb1,
43 C & t1,q1,qs1,u1,v1,gz1,tv1,tp1,tvp1,clw1,h1,
44 C & lv1,cpn1,p1,ph1,ft1,fq1,fu1,fv1,
45 C & tnk1,qnk1,gznk1,plcl1,
46 C & precip1,cbmf1,iflag1,
47 C & delt,cpd,cpv,cl,rv,rd,lv0,g,
48 C & sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp,
49 C & alpha,entp,coeffs,coeffr,omtrain,cu)
50 C
51 C PARAMETERS:
52 C Name Type Usage Description
53 C ---------- ---------- ------- ----------------------------
54 C
55 C ncum Integer Input number of cumulus points
56 C idcum Integer Input index of cumulus point
57 C len Integer Input first dimension
58 C nd Integer Input total vertical dimension
59 C ndp1 Integer Input nd + 1
60 C nl Integer Input vertical dimension for convection
61 C minorig Integer Input First level where convection is allow to begin
62 C nk1 Integer Input First level of convection
63 C ncb1 Integer Input Level of free convection
64 C t1 Real Input temperature
65 C q1 Real Input specific hum
66 C qs1 Real Input sat specific hum
67 C u1 Real Input u-wind
68 C v1 Real Input v-wind
69 C gz1 Real Inout geop
70 C tv1 Real Input virtual temp
71 C tp1 Real Input
72 C clw1 Real Inout cloud liquid water
73 C h1 Real Inout
74 C lv1 Real Inout
75 C cpn1 Real Inout
76 C p1 Real Input full level pressure
77 C ph1 Real Input half level pressure
78 C ft1 Real Output temp tend
79 C fq1 Real Output spec hum tend
80 C fu1 Real Output u-wind tend
81 C fv1 Real Output v-wind tend
82 C precip1 Real Output prec
83 C cbmf1 Real In/Out cumulus mass flux
84 C iflag1 Integer Output iflag on latitude strip
85 C delt Real Input time step
86 C cpd Integer Input See description below
87 C cpv Integer Input See description below
88 C cl Integer Input See description below
89 C rv Integer Input See description below
90 C rd Integer Input See description below
91 C lv0 Integer Input See description below
92 C g Integer Input See description below
93 C sigs Integer Input See description below
94 C sigd Integer Input See description below
95 C elcrit Integer Input See description below
96 C tlcrit Integer Input See description below
97 C omtsnow Integer Input See description below
98 C dtmax Integer Input See description below
99 C damp Integer Input See description below
100 C alpha Integer Input See description below
101 C ent Integer Input See description below
102 C coeffs Integer Input See description below
103 C coeffr Integer Input See description below
104 C omtrain Integer Input See description below
105 C cu Integer Input See description below
106 C
107 C COMMON BLOCKS:
108 C Block Name Type Usage Notes
109 C -------- -------- ---- ------ ------------------------
110 C
111 C FILES: None
112 C
113 C DATA BASES: None
114 C
115 C NON-FILE INPUT/OUTPUT: None
116 C
117 C ERROR CONDITIONS: None
118 C
119 C ADDITIONAL COMMENTS: None
120 C
121 C.................MAINTENANCE SECTION................................
122 C
123 C MODULES CALLED:
124 C Name Description
125 C zilch Zero out an array
126 C ------- ----------------------
127 C LOCAL VARIABLES AND
128 C STRUCTURES:
129 C Name Type Description
130 C ------- ------ -----------
131 C See Comments Below
132 C
133 C i Integer loop index
134 C k Integer loop index
135 c
136 C METHOD:
137 C
138 C See Emanuel, K. and M. Zivkovic-Rothman, 2000: Development and evaluation of a
139 C convective scheme for use in climate models.
140 C
141 C FILES: None
142 C
143 C INCLUDE FILES: None
144 C
145 C MAKEFILE: /a/ops/met/nogaps/src/sub/fcst/fcst159lib.mak
146 C
147 C..............................END PROLOGUE.............................
148 c
149 c
150  USE dimphy
151  implicit none
152 c
153 cym#include "dimensions.h"
154 cym#include "dimphy.h"
155 c
156  integer kmax2,imax2,kmin2,imin2
157  real ftmax2,ftmin2
158  integer kmax,imax,kmin,imin
159  real ftmax,ftmin
160 c
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-
208 c 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
211 c 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
219 c
220  real ma(len,nd)
221 c
222 C
223 C *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***
224 C *** TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO- ***
225 C *** CONVERSION THRESHOLD IS ASSUMED TO BE ZERO ***
226 C *** (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY ***
227 C *** BETWEEN 0 C AND TLCRIT) ***
228 C *** ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT ***
229 C *** FORMULATION ***
230 C *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT ***
231 C *** SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE ***
232 C *** OF CLOUD ***
233 C *** OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN ***
234 C *** OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW ***
235 C *** COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION ***
236 C *** OF RAIN ***
237 C *** COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION ***
238 C *** OF SNOW ***
239 C *** CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM ***
240 C *** TRANSPORT ***
241 C *** DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION ***
242 C *** A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC ***
243 C *** ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF ***
244 C *** APPROACH TO QUASI-EQUILIBRIUM ***
245 C *** (THEIR STANDARD VALUES ARE 0.20 AND 0.1, RESPECTIVELY) ***
246 C *** (DAMP MUST BE LESS THAN 1) ***
247 c
248 c Local arrays.
249 c
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)
276 c real det(ncum,klev)
277  real dph(ncum,klev)
278 c real wd(ncum)
279 c real tprime(ncum)
280 c 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)
314 c
315  integer iflag(ncum)
316  integer nk(ncum)
317  integer icb(ncum)
318  integer inb(ncum)
319  integer inb1(ncum)
320  integer jtt(ncum)
321 c
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
336 c 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)
352 c real amp1,ad
353  real amp1(ncum),ad(ncum)
354  logical lwork(ncum)
355  integer num1,num2
356 c
357 c print*,'cpd en entree de convect2 ',cpd
358  nlp=nl+1
359 c
360  rowl=1000.0
361  ginv=1.0/g
362  delti=1.0/delt
363 c
364 c Define some thermodynamic variables.
365 c
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
373 c
374 c Compress the fields.
375 c
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
398 c 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
414 c print*,'150 ncum,nn',ncum,nn
415 c
416 c Initialize the tendencies, det, wd, tprime, qprime.
417 c
418  do 170 k=1,nl
419  do 160 i=1,ncum
420 c 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
431 c wd(i)=0.0
432 c tprime(i)=0.0
433 c 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
440 c
441 c Compute icbmax.
442 c
443  icbmax=2
444  do 230 i=1,ncum
445  icbmax=max(icbmax,icb(i))
446  230 continue
447 c
448 c
449 !=====================================================================
450 ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
451 !=====================================================================
452 c
453 c --- The procedure is to solve the equation.
454 c cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
455 c
456 c *** Calculate certain parcel quantities, including static energy ***
457 c
458 c
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
463 c
464 c
465 c *** Find lifted parcel quantities above cloud base ***
466 c
467 c
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)
474 c
475 c First iteration.
476 c
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))
490 c
491 c Second iteration.
492 c
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))
506 c
507  alv=lv0-clmcpv*(t(i,k)-273.15)
508 c print*,'cpd dans convect2 ',cpd
509 c print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
510 c 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
513 c if (.not.cpd.gt.1000.) then
514 c print*,'CPD=',cpd
515 c stop
516 c 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
524 c
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 !=====================================================================
530 c
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
548 c
549 !=====================================================================
550 ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
551 ! --- VIRTUAL TEMPERATURE
552 !=====================================================================
553 c
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))
558 c print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
559 c 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
566 c
567 c
568 c=====================================================================
569 c --- NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
570 c=====================================================================
571 c
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)
580 c if(.not.cpn(i,1).gt.900.) then
581 c print*,'i,lv(i,1),cpn(i,1)'
582 c print*, i,lv(i,1),cpn(i,1)
583 c k=(i-1)/ncum+1
584 c print*,'i,k',mod(i,ncum),k,' cpn',cpn(mod(i,ncum),k)
585 c stop
586 c endif
587  lvcp(i,1)=lv(i,1)/cpn(i,1)
588  360 continue
589 c
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
595 c
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
605 c
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
618 c
619 c=====================================================================
620 c --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
621 c --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
622 c --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
623 c=====================================================================
624 c
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
631 c
632 c Originial Code
633 c
634 c do 530 k=minorig+1,nl-1
635 c do 520 i=1,ncum
636 c if(k.ge.(icb(i)+1))then
637 c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
638 c byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
639 c cape(i)=cape(i)+by
640 c if(by.ge.0.0)inb1(i)=k+1
641 c if(cape(i).gt.0.0)then
642 c inb(i)=k+1
643 c capem(i)=cape(i)
644 c endif
645 c endif
646 c520 continue
647 c530 continue
648 c do 540 i=1,ncum
649 c byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
650 c cape(i)=capem(i)+byp
651 c defrac=capem(i)-cape(i)
652 c defrac=max(defrac,0.001)
653 c frac(i)=-cape(i)/defrac
654 c frac(i)=min(frac(i),1.0)
655 c frac(i)=max(frac(i),0.0)
656 c540 continue
657 c
658 c K Emanuel fix
659 c
660 c call zilch(byp,ncum)
661 c do 530 k=minorig+1,nl-1
662 c do 520 i=1,ncum
663 c if(k.ge.(icb(i)+1))then
664 c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
665 c cape(i)=cape(i)+by
666 c if(by.ge.0.0)inb1(i)=k+1
667 c if(cape(i).gt.0.0)then
668 c inb(i)=k+1
669 c capem(i)=cape(i)
670 c byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
671 c endif
672 c endif
673 c520 continue
674 c530 continue
675 c do 540 i=1,ncum
676 c inb(i)=max(inb(i),inb1(i))
677 c cape(i)=capem(i)+byp(i)
678 c defrac=capem(i)-cape(i)
679 c defrac=max(defrac,0.001)
680 c frac(i)=-cape(i)/defrac
681 c frac(i)=min(frac(i),1.0)
682 c frac(i)=max(frac(i),0.0)
683 c540 continue
684 c
685 c J Teixeira fix
686 c
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
714 c
715 c=====================================================================
716 c --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
717 c=====================================================================
718 c
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
726 c
727 c=====================================================================
728 c --- CALCULATE CLOUD BASE MASS FLUX AND RATES OF MIXING, M(I),
729 c --- AT EACH MODEL LEVEL
730 c=====================================================================
731 c
732 c tvpplcl = parcel temperature lifted adiabatically from level
733 c icb-1 to the LCL.
734 c tvaplcl = virtual temperature at the LCL.
735 c
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
745 c
746 c-------------------------------------------------------------------
747 c --- Interpolate difference between lifted parcel and
748 c --- environmental temperatures to lifted condensation level
749 c-------------------------------------------------------------------
750 c
751 c dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
752 c
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
764 c
765 c-------------------------------------------------------------------
766 c --- Adjust cloud base mass flux
767 c-------------------------------------------------------------------
768 c
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
776 c
777 c-------------------------------------------------------------------
778 c --- Calculate rates of mixing, m(i)
779 c-------------------------------------------------------------------
780 c
781  call zilch(work,ncum)
782 c
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
801 c
802 c
803 c=====================================================================
804 c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
805 c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
806 c --- FRACTION (sij)
807 c=====================================================================
808 c
809 c
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
854 c
855 c *** If no air can entrain at level i assume that updraft detrains ***
856 c *** at that level and calculate detrained air flux and properties ***
857 c
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
870 c
871  do 770 i=1,ncum
872  sij(i,inb(i),inb(i))=1.0
873  770 continue
874 c
875 c=====================================================================
876 c --- NORMALIZE ENTRAINED AIR MASS FLUXES
877 c --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
878 c=====================================================================
879 c
880 c
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
886 c
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
892 c
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
908 c
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
916 c
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
979 c
980 c=====================================================================
981 c --- PRECIPITATING DOWNDRAFT CALCULATION
982 c=====================================================================
983 c
984 c *** Check whether ep(inb)=0, if so, skip precipitating ***
985 c *** downdraft calculation ***
986 c
987 c
988 c *** Integrate liquid water equation to find condensed water ***
989 c *** and condensed water flux ***
990 c
991 c
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
1001 c
1002 c *** Begin downdraft loop ***
1003 c
1004 c
1005  call zilch(wdtrain,ncum)
1006  do 899 i=nl+1,1,-1
1007 c
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
1013 c
1014 c
1015 c *** Calculate detrained precipitation ***
1016 c
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
1022 c
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
1034 c
1035 c *** Find rain water and evaporation using provisional ***
1036 c *** estimates of qp(i)and qp(i-1) ***
1037 c
1038 c
1039 c *** Value of terminal velocity and coeffecient of evaporation for snow ***
1040 c
1041  do 894 ij=1,ncum
1042  if((i.le.inb(ij)).and.(lwork(ij)))then
1043  coeff=coeffs
1044  wt(ij,i)=omtsnow
1045 c
1046 c *** Value of terminal velocity and coeffecient of evaporation for rain ***
1047 c
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
1064 c
1065 c *** Calculate precipitating downdraft mass flux under ***
1066 c *** hydrostatic approximation ***
1067 c
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)
1073 c
1074 c *** Add small amount of inertia to downdraft ***
1075 c
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)
1078 c
1079 c *** Force mp to decrease linearly to zero ***
1080 c *** between about 950 mb and the surface ***
1081 c
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
1088 c
1089 c *** Find mixing ratio of precipitating downdraft ***
1090 c
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
1119 c
1120 c *** Calculate surface precipitation in mm/day ***
1121 c
1122  do 1190 i=1,ncum
1123  if(iflag(i).le.1)then
1124 cc precip(i)=precip(i)+wt(i,1)*sigd*water(i,1)*3600.*24000.
1125 cc & /(rowl*g)
1126 cc precip(i)=precip(i)*delt/86400.
1127  precip(i) = wt(i,1)*sigd*water(i,1)*86400/g
1128  endif
1129  1190 continue
1130 c
1131 c
1132 c *** Calculate downdraft velocity scale and surface temperature and ***
1133 c *** water vapor fluctuations ***
1134 c
1135 c wd=beta*abs(mp(icb))*0.01*rd*t(icb)/(sigd*p(icb))
1136 c qprime=0.5*(qp(1)-q(1))
1137 c tprime=lv0*qprime/cpd
1138 c
1139 c *** Calculate tendencies of lowest level potential temperature ***
1140 c *** and mixing ratio ***
1141 c
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
1180 c
1181 c *** Calculate tendencies of potential temperature and mixing ratio ***
1182 c *** at levels above the lowest level ***
1183 c
1184 c *** First find the net saturated updraft and downdraft mass fluxes ***
1185 c *** through each level ***
1186 c
1187  do 1500 i=2,nl+1
1188 c
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
1194 c
1195  call zilch(amp1,ncum)
1196  call zilch(ad,ncum)
1197 c
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
1206 c
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
1225 c
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)
1230 c
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
1289 c
1290 c *** Adjust tendencies at top of convection layer to reflect ***
1291 c *** actual position of the level zero cape ***
1292 c
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
1317 c
1318 c *** Very slightly adjust tendencies to force exact ***
1319 c *** enthalpy, momentum and tracer conservation ***
1320 c
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
1344 c
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
1350 c
1351 c
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
1381 c
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
1392 c
1393  return
1394  end
1395