source: LMDZ4/trunk/libf/phylmd/cv_routines.F @ 1227

Last change on this file since 1227 was 524, checked in by lmdzadmin, 20 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 55.0 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE cv_param(nd)
5      implicit none
6
7c------------------------------------------------------------
8c Set parameters for convectL
9c (includes microphysical parameters and parameters that
10c  control the rate of approach to quasi-equilibrium)
11c------------------------------------------------------------
12
13C   *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***
14C   ***  TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-        ***
15C   ***       CONVERSION THRESHOLD IS ASSUMED TO BE ZERO             ***
16C   ***     (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY            ***
17C   ***               BETWEEN 0 C AND TLCRIT)                        ***
18C   ***   ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT       ***
19C   ***                       FORMULATION                            ***
20C   ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
21C   ***  SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE       ***
22C   ***                        OF CLOUD                              ***
23C   ***        OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN       ***
24C   ***     OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW          ***
25C   ***  COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
26C   ***                          OF RAIN                             ***
27C   ***  COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
28C   ***                          OF SNOW                             ***
29C   ***     CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM      ***
30C   ***                         TRANSPORT                            ***
31C   ***    DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION    ***
32C   ***        A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC      ***
33C   ***    ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF    ***
34C   ***                 APPROACH TO QUASI-EQUILIBRIUM                ***
35C   ***   (THEIR STANDARD VALUES ARE  0.20 AND 0.1, RESPECTIVELY)    ***
36C   ***                   (DAMP MUST BE LESS THAN 1)                 ***
37
38#include "cvparam.h"
39      integer nd
40
41c noff: integer limit for convection (nd-noff)
42c minorig: First level of convection
43
44      noff = 2
45      minorig = 2
46
47      nl=nd-noff
48      nlp=nl+1
49      nlm=nl-1
50
51      elcrit=0.0011
52      tlcrit=-55.0
53      entp=1.5
54      sigs=0.12
55      sigd=0.05
56      omtrain=50.0
57      omtsnow=5.5
58      coeffr=1.0
59      coeffs=0.8
60      dtmax=0.9
61c
62      cu=0.70
63c
64      betad=10.0
65c
66      damp=0.1
67      alpha=0.2
68c
69      delta=0.01  ! cld
70c
71      return
72      end
73
74      SUBROUTINE cv_prelim(len,nd,ndp1,t,q,p,ph
75     :                    ,lv,cpn,tv,gz,h,hm)
76      implicit none
77
78!=====================================================================
79! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
80!=====================================================================
81
82c inputs:
83      integer len, nd, ndp1
84      real t(len,nd), q(len,nd), p(len,nd), ph(len,ndp1)
85
86c outputs:
87      real lv(len,nd), cpn(len,nd), tv(len,nd)
88      real gz(len,nd), h(len,nd), hm(len,nd)
89
90c local variables:
91      integer k, i
92      real cpx(len,nd)
93
94#include "cvthermo.h"
95#include "cvparam.h"
96
97
98      do 110 k=1,nlp
99        do 100 i=1,len
100          lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
101          cpn(i,k)=cpd*(1.0-q(i,k))+cpv*q(i,k)
102          cpx(i,k)=cpd*(1.0-q(i,k))+cl*q(i,k)
103          tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1)
104 100    continue
105 110  continue
106c
107c gz = phi at the full levels (same as p).
108c
109      do 120 i=1,len
110        gz(i,1)=0.0
111 120  continue
112      do 140 k=2,nlp
113        do 130 i=1,len
114          gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
115     &         *(p(i,k-1)-p(i,k))/ph(i,k)
116 130    continue
117 140  continue
118c
119c h  = phi + cpT (dry static energy).
120c hm = phi + cp(T-Tbase)+Lq
121c
122      do 170 k=1,nlp
123        do 160 i=1,len
124          h(i,k)=gz(i,k)+cpn(i,k)*t(i,k)
125          hm(i,k)=gz(i,k)+cpx(i,k)*(t(i,k)-t(i,1))+lv(i,k)*q(i,k)
126 160    continue
127 170  continue
128
129      return
130      end
131
132      SUBROUTINE cv_feed(len,nd,t,q,qs,p,hm,gz
133     :                  ,nk,icb,icbmax,iflag,tnk,qnk,gznk,plcl)
134      implicit none
135
136C================================================================
137C Purpose: CONVECTIVE FEED
138C================================================================
139
140#include "cvparam.h"
141
142c inputs:
143          integer len, nd
144      real t(len,nd), q(len,nd), qs(len,nd), p(len,nd)
145      real hm(len,nd), gz(len,nd)
146
147c outputs:
148          integer iflag(len), nk(len), icb(len), icbmax
149      real tnk(len), qnk(len), gznk(len), plcl(len)
150
151c local variables:
152      integer i, k
153      integer ihmin(len)
154      real work(len)
155      real pnk(len), qsnk(len), rh(len), chi(len)
156
157!-------------------------------------------------------------------
158! --- Find level of minimum moist static energy
159! --- If level of minimum moist static energy coincides with
160! --- or is lower than minimum allowable parcel origin level,
161! --- set iflag to 6.
162!-------------------------------------------------------------------
163
164      do 180 i=1,len
165       work(i)=1.0e12
166       ihmin(i)=nl
167 180  continue
168      do 200 k=2,nlp
169        do 190 i=1,len
170         if((hm(i,k).lt.work(i)).and.
171     &      (hm(i,k).lt.hm(i,k-1)))then
172           work(i)=hm(i,k)
173           ihmin(i)=k
174         endif
175 190    continue
176 200  continue
177      do 210 i=1,len
178        ihmin(i)=min(ihmin(i),nlm)
179        if(ihmin(i).le.minorig)then
180          iflag(i)=6
181        endif
182 210  continue
183c
184!-------------------------------------------------------------------
185! --- Find that model level below the level of minimum moist static
186! --- energy that has the maximum value of moist static energy
187!-------------------------------------------------------------------
188 
189      do 220 i=1,len
190       work(i)=hm(i,minorig)
191       nk(i)=minorig
192 220  continue
193      do 240 k=minorig+1,nl
194        do 230 i=1,len
195         if((hm(i,k).gt.work(i)).and.(k.le.ihmin(i)))then
196           work(i)=hm(i,k)
197           nk(i)=k
198         endif
199 230     continue
200 240  continue
201!-------------------------------------------------------------------
202! --- Check whether parcel level temperature and specific humidity
203! --- are reasonable
204!-------------------------------------------------------------------
205       do 250 i=1,len
206       if(((t(i,nk(i)).lt.250.0).or.
207     &      (q(i,nk(i)).le.0.0).or.
208     &      (p(i,ihmin(i)).lt.400.0)).and.
209     &      (iflag(i).eq.0))iflag(i)=7
210 250   continue
211!-------------------------------------------------------------------
212! --- Calculate lifted condensation level of air at parcel origin level
213! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
214!-------------------------------------------------------------------
215       do 260 i=1,len
216        tnk(i)=t(i,nk(i))
217        qnk(i)=q(i,nk(i))
218        gznk(i)=gz(i,nk(i))
219        pnk(i)=p(i,nk(i))
220        qsnk(i)=qs(i,nk(i))
221c
222        rh(i)=qnk(i)/qsnk(i)
223        rh(i)=min(1.0,rh(i))
224        chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
225        plcl(i)=pnk(i)*(rh(i)**chi(i))
226        if(((plcl(i).lt.200.0).or.(plcl(i).ge.2000.0))
227     &   .and.(iflag(i).eq.0))iflag(i)=8
228 260   continue
229!-------------------------------------------------------------------
230! --- Calculate first level above lcl (=icb)
231!-------------------------------------------------------------------
232      do 270 i=1,len
233       icb(i)=nlm
234 270  continue
235c
236      do 290 k=minorig,nl
237        do 280 i=1,len
238          if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i)))
239     &    icb(i)=min(icb(i),k)
240 280    continue
241 290  continue
242c
243      do 300 i=1,len
244        if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
245 300  continue
246c
247c Compute icbmax.
248c
249      icbmax=2
250      do 310 i=1,len
251        icbmax=max(icbmax,icb(i))
252 310  continue
253
254      return
255      end
256
257      SUBROUTINE cv_undilute1(len,nd,t,q,qs,gz,p,nk,icb,icbmax
258     :                       ,tp,tvp,clw)
259      implicit none
260
261#include "cvthermo.h"
262#include "cvparam.h"
263
264c inputs:
265      integer len, nd
266      integer nk(len), icb(len), icbmax
267      real t(len,nd), q(len,nd), qs(len,nd), gz(len,nd)
268      real p(len,nd)
269
270c outputs:
271      real tp(len,nd), tvp(len,nd), clw(len,nd)
272
273c local variables:
274      integer i, k
275      real tg, qg, alv, s, ahg, tc, denom, es, rg
276      real ah0(len), cpp(len)
277      real tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)
278
279!-------------------------------------------------------------------
280! --- Calculates the lifted parcel virtual temperature at nk,
281! --- the actual temperature, and the adiabatic
282! --- liquid water content. The procedure is to solve the equation.
283!     cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
284!-------------------------------------------------------------------
285
286      do 320 i=1,len
287        tnk(i)=t(i,nk(i))
288        qnk(i)=q(i,nk(i))
289        gznk(i)=gz(i,nk(i))
290        ticb(i)=t(i,icb(i))
291        gzicb(i)=gz(i,icb(i))
292 320  continue
293c
294c   ***  Calculate certain parcel quantities, including static energy   ***
295c
296      do 330 i=1,len
297        ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i)
298     &         +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i)
299        cpp(i)=cpd*(1.-qnk(i))+qnk(i)*cpv
300 330  continue
301c
302c   ***   Calculate lifted parcel quantities below cloud base   ***
303c
304        do 350 k=minorig,icbmax-1
305          do 340 i=1,len
306           tp(i,k)=tnk(i)-(gz(i,k)-gznk(i))/cpp(i)
307           tvp(i,k)=tp(i,k)*(1.+qnk(i)*epsi)
308  340     continue
309  350   continue
310c
311c    ***  Find lifted parcel quantities above cloud base    ***
312c
313        do 360 i=1,len
314         tg=ticb(i)
315         qg=qs(i,icb(i))
316         alv=lv0-clmcpv*(ticb(i)-t0)
317c
318c First iteration.
319c
320          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
321          s=1./s
322          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
323          tg=tg+s*(ah0(i)-ahg)
324          tg=max(tg,35.0)
325          tc=tg-t0
326          denom=243.5+tc
327          if(tc.ge.0.0)then
328           es=6.112*exp(17.67*tc/denom)
329          else
330           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
331          endif
332          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
333c
334c Second iteration.
335c
336          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
337          s=1./s
338          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
339          tg=tg+s*(ah0(i)-ahg)
340          tg=max(tg,35.0)
341          tc=tg-t0
342          denom=243.5+tc
343          if(tc.ge.0.0)then
344           es=6.112*exp(17.67*tc/denom)
345          else
346           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
347          end if
348          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
349c
350         alv=lv0-clmcpv*(ticb(i)-273.15)
351         tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
352     &   -gz(i,icb(i))-alv*qg)/cpd
353         clw(i,icb(i))=qnk(i)-qg
354         clw(i,icb(i))=max(0.0,clw(i,icb(i)))
355         rg=qg/(1.-qnk(i))
356         tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
357  360   continue
358c
359      do 380 k=minorig,icbmax
360       do 370 i=1,len
361         tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
362 370   continue
363 380  continue
364c
365      return
366      end
367
368      SUBROUTINE cv_trigger(len,nd,icb,cbmf,tv,tvp,iflag)
369      implicit none
370
371!-------------------------------------------------------------------
372! --- Test for instability.
373! --- If there was no convection at last time step and parcel
374! --- is stable at icb, then set iflag to 4.
375!-------------------------------------------------------------------
376 
377#include "cvparam.h"
378
379c inputs:
380       integer len, nd, icb(len)
381       real cbmf(len), tv(len,nd), tvp(len,nd)
382
383c outputs:
384       integer iflag(len) ! also an input
385
386c local variables:
387       integer i
388
389
390      do 390 i=1,len
391        if((cbmf(i).eq.0.0) .and.(iflag(i).eq.0).and.
392     &  (tvp(i,icb(i)).le.(tv(i,icb(i))-dtmax)))iflag(i)=4
393 390  continue
394 
395      return
396      end
397
398      SUBROUTINE cv_compress( len,nloc,ncum,nd
399     :   ,iflag1,nk1,icb1
400     :   ,cbmf1,plcl1,tnk1,qnk1,gznk1
401     :   ,t1,q1,qs1,u1,v1,gz1
402     :   ,h1,lv1,cpn1,p1,ph1,tv1,tp1,tvp1,clw1
403     o   ,iflag,nk,icb
404     o   ,cbmf,plcl,tnk,qnk,gznk
405     o   ,t,q,qs,u,v,gz,h,lv,cpn,p,ph,tv,tp,tvp,clw
406     o   ,dph          )
407      implicit none
408
409#include "cvparam.h"
410
411c inputs:
412      integer len,ncum,nd,nloc
413      integer iflag1(len),nk1(len),icb1(len)
414      real cbmf1(len),plcl1(len),tnk1(len),qnk1(len),gznk1(len)
415      real t1(len,nd),q1(len,nd),qs1(len,nd),u1(len,nd),v1(len,nd)
416      real gz1(len,nd),h1(len,nd),lv1(len,nd),cpn1(len,nd)
417      real p1(len,nd),ph1(len,nd+1),tv1(len,nd),tp1(len,nd)
418      real tvp1(len,nd),clw1(len,nd)
419
420c outputs:
421      integer iflag(nloc),nk(nloc),icb(nloc)
422      real cbmf(nloc),plcl(nloc),tnk(nloc),qnk(nloc),gznk(nloc)
423      real t(nloc,nd),q(nloc,nd),qs(nloc,nd),u(nloc,nd),v(nloc,nd)
424      real gz(nloc,nd),h(nloc,nd),lv(nloc,nd),cpn(nloc,nd)
425      real p(nloc,nd),ph(nloc,nd+1),tv(nloc,nd),tp(nloc,nd)
426      real tvp(nloc,nd),clw(nloc,nd)
427      real dph(nloc,nd)
428
429c local variables:
430      integer i,k,nn
431
432
433      do 110 k=1,nl+1
434       nn=0
435      do 100 i=1,len
436      if(iflag1(i).eq.0)then
437        nn=nn+1
438        t(nn,k)=t1(i,k)
439        q(nn,k)=q1(i,k)
440        qs(nn,k)=qs1(i,k)
441        u(nn,k)=u1(i,k)
442        v(nn,k)=v1(i,k)
443        gz(nn,k)=gz1(i,k)
444        h(nn,k)=h1(i,k)
445        lv(nn,k)=lv1(i,k)
446        cpn(nn,k)=cpn1(i,k)
447        p(nn,k)=p1(i,k)
448        ph(nn,k)=ph1(i,k)
449        tv(nn,k)=tv1(i,k)
450        tp(nn,k)=tp1(i,k)
451        tvp(nn,k)=tvp1(i,k)
452        clw(nn,k)=clw1(i,k)
453      endif
454 100    continue
455 110  continue
456
457      if (nn.ne.ncum) then
458         print*,'strange! nn not equal to ncum: ',nn,ncum
459         stop
460      endif
461
462      nn=0
463      do 150 i=1,len
464      if(iflag1(i).eq.0)then
465      nn=nn+1
466      cbmf(nn)=cbmf1(i)
467      plcl(nn)=plcl1(i)
468      tnk(nn)=tnk1(i)
469      qnk(nn)=qnk1(i)
470      gznk(nn)=gznk1(i)
471      nk(nn)=nk1(i)
472      icb(nn)=icb1(i)
473      iflag(nn)=iflag1(i)
474      endif
475 150  continue
476
477      do 170 k=1,nl
478       do 160 i=1,ncum
479        dph(i,k)=ph(i,k)-ph(i,k+1)
480 160   continue
481 170  continue
482
483      return
484      end
485
486      SUBROUTINE cv_undilute2(nloc,ncum,nd,icb,nk
487     :                       ,tnk,qnk,gznk,t,q,qs,gz
488     :                       ,p,dph,h,tv,lv
489     o                       ,inb,inb1,tp,tvp,clw,hp,ep,sigp,frac)
490      implicit none
491
492C---------------------------------------------------------------------
493C Purpose:
494C     FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
495C     &
496C     COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
497C     FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
498C     &
499C     FIND THE LEVEL OF NEUTRAL BUOYANCY
500C---------------------------------------------------------------------
501
502#include "cvthermo.h"
503#include "cvparam.h"
504
505c inputs:
506      integer ncum, nd, nloc
507      integer icb(nloc), nk(nloc)
508      real t(nloc,nd), q(nloc,nd), qs(nloc,nd), gz(nloc,nd)
509      real p(nloc,nd), dph(nloc,nd)
510      real tnk(nloc), qnk(nloc), gznk(nloc)
511      real lv(nloc,nd), tv(nloc,nd), h(nloc,nd)
512
513c outputs:
514      integer inb(nloc), inb1(nloc)
515      real tp(nloc,nd), tvp(nloc,nd), clw(nloc,nd)
516      real ep(nloc,nd), sigp(nloc,nd), hp(nloc,nd)
517      real frac(nloc)
518
519c local variables:
520      integer i, k
521      real tg,qg,ahg,alv,s,tc,es,denom,rg,tca,elacrit
522      real by, defrac
523      real ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
524      logical lcape(nloc)
525
526!=====================================================================
527! --- SOME INITIALIZATIONS
528!=====================================================================
529
530      do 170 k=1,nl
531      do 160 i=1,ncum
532       ep(i,k)=0.0
533       sigp(i,k)=sigs
534 160  continue
535 170  continue
536
537!=====================================================================
538! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
539!=====================================================================
540c
541c ---       The procedure is to solve the equation.
542c              cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
543c
544c   ***  Calculate certain parcel quantities, including static energy   ***
545c
546c
547      do 240 i=1,ncum
548         ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i)
549     &         +qnk(i)*(lv0-clmcpv*(tnk(i)-t0))+gznk(i)
550 240  continue
551c
552c
553c    ***  Find lifted parcel quantities above cloud base    ***
554c
555c
556        do 300 k=minorig+1,nl
557          do 290 i=1,ncum
558            if(k.ge.(icb(i)+1))then
559              tg=t(i,k)
560              qg=qs(i,k)
561              alv=lv0-clmcpv*(t(i,k)-t0)
562c
563c First iteration.
564c
565               s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
566               s=1./s
567               ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
568               tg=tg+s*(ah0(i)-ahg)
569               tg=max(tg,35.0)
570               tc=tg-t0
571               denom=243.5+tc
572               if(tc.ge.0.0)then
573                        es=6.112*exp(17.67*tc/denom)
574               else
575                        es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
576               endif
577                        qg=eps*es/(p(i,k)-es*(1.-eps))
578c
579c Second iteration.
580c
581               s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
582               s=1./s
583               ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
584               tg=tg+s*(ah0(i)-ahg)
585               tg=max(tg,35.0)
586               tc=tg-t0
587               denom=243.5+tc
588               if(tc.ge.0.0)then
589                        es=6.112*exp(17.67*tc/denom)
590               else
591                        es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
592               endif
593                        qg=eps*es/(p(i,k)-es*(1.-eps))
594c
595               alv=lv0-clmcpv*(t(i,k)-t0)
596c      print*,'cpd dans convect2 ',cpd
597c      print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
598c      print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
599        tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
600c              if (.not.cpd.gt.1000.) then
601c                  print*,'CPD=',cpd
602c                  stop
603c              endif
604               clw(i,k)=qnk(i)-qg
605               clw(i,k)=max(0.0,clw(i,k))
606               rg=qg/(1.-qnk(i))
607               tvp(i,k)=tp(i,k)*(1.+rg*epsi)
608            endif
609  290     continue
610  300   continue
611c
612!=====================================================================
613! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
614! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
615! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
616!=====================================================================
617c
618      do 320 k=minorig+1,nl
619        do 310 i=1,ncum
620          if(k.ge.(nk(i)+1))then
621            tca=tp(i,k)-t0
622            if(tca.ge.0.0)then
623              elacrit=elcrit
624            else
625              elacrit=elcrit*(1.0-tca/tlcrit)
626            endif
627            elacrit=max(elacrit,0.0)
628            ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)
629            ep(i,k)=max(ep(i,k),0.0 )
630            ep(i,k)=min(ep(i,k),1.0 )
631            sigp(i,k)=sigs
632          endif
633 310    continue
634 320  continue
635c
636!=====================================================================
637! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
638! --- VIRTUAL TEMPERATURE
639!=====================================================================
640c
641      do 340 k=minorig+1,nl
642        do 330 i=1,ncum
643        if(k.ge.(icb(i)+1))then
644          tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
645c         print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
646c         print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
647        endif
648 330    continue
649 340  continue
650      do 350 i=1,ncum
651       tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
652 350  continue
653c
654c=====================================================================
655c  --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
656c  --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
657c  --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
658c=====================================================================
659c
660      do 510 i=1,ncum
661        cape(i)=0.0
662        capem(i)=0.0
663        inb(i)=icb(i)+1
664        inb1(i)=inb(i)
665 510  continue
666c
667c Originial Code
668c
669c     do 530 k=minorig+1,nl-1
670c       do 520 i=1,ncum
671c         if(k.ge.(icb(i)+1))then
672c           by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
673c           byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
674c           cape(i)=cape(i)+by
675c           if(by.ge.0.0)inb1(i)=k+1
676c           if(cape(i).gt.0.0)then
677c             inb(i)=k+1
678c             capem(i)=cape(i)
679c           endif
680c         endif
681c520    continue
682c530  continue
683c     do 540 i=1,ncum
684c         byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
685c         cape(i)=capem(i)+byp
686c         defrac=capem(i)-cape(i)
687c         defrac=max(defrac,0.001)
688c         frac(i)=-cape(i)/defrac
689c         frac(i)=min(frac(i),1.0)
690c         frac(i)=max(frac(i),0.0)
691c540   continue
692c
693c K Emanuel fix
694c
695c     call zilch(byp,ncum)
696c     do 530 k=minorig+1,nl-1
697c       do 520 i=1,ncum
698c         if(k.ge.(icb(i)+1))then
699c           by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
700c           cape(i)=cape(i)+by
701c           if(by.ge.0.0)inb1(i)=k+1
702c           if(cape(i).gt.0.0)then
703c             inb(i)=k+1
704c             capem(i)=cape(i)
705c             byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
706c           endif
707c         endif
708c520    continue
709c530  continue
710c     do 540 i=1,ncum
711c         inb(i)=max(inb(i),inb1(i))
712c         cape(i)=capem(i)+byp(i)
713c         defrac=capem(i)-cape(i)
714c         defrac=max(defrac,0.001)
715c         frac(i)=-cape(i)/defrac
716c         frac(i)=min(frac(i),1.0)
717c         frac(i)=max(frac(i),0.0)
718c540   continue
719c
720c J Teixeira fix
721c
722      call zilch(byp,ncum)
723      do 515 i=1,ncum
724        lcape(i)=.true.
725 515  continue
726      do 530 k=minorig+1,nl-1
727        do 520 i=1,ncum
728          if(cape(i).lt.0.0)lcape(i)=.false.
729          if((k.ge.(icb(i)+1)).and.lcape(i))then
730            by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
731            byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
732            cape(i)=cape(i)+by
733            if(by.ge.0.0)inb1(i)=k+1
734            if(cape(i).gt.0.0)then
735              inb(i)=k+1
736              capem(i)=cape(i)
737            endif
738          endif
739 520    continue
740 530  continue
741      do 540 i=1,ncum
742          cape(i)=capem(i)+byp(i)
743          defrac=capem(i)-cape(i)
744          defrac=max(defrac,0.001)
745          frac(i)=-cape(i)/defrac
746          frac(i)=min(frac(i),1.0)
747          frac(i)=max(frac(i),0.0)
748 540  continue
749c
750c=====================================================================
751c ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
752c=====================================================================
753c
754c initialization:
755      do i=1,ncum*nlp
756       hp(i,1)=h(i,1)
757      enddo
758
759      do 600 k=minorig+1,nl
760        do 590 i=1,ncum
761        if((k.ge.icb(i)).and.(k.le.inb(i)))then
762          hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k)
763        endif
764 590    continue
765 600  continue
766c
767        return
768        end
769c
770      SUBROUTINE cv_closure(nloc,ncum,nd,nk,icb
771     :                     ,tv,tvp,p,ph,dph,plcl,cpn
772     :                     ,iflag,cbmf)
773      implicit none
774
775c inputs:
776      integer ncum, nd, nloc
777      integer nk(nloc), icb(nloc)
778      real tv(nloc,nd), tvp(nloc,nd), p(nloc,nd), dph(nloc,nd)
779      real ph(nloc,nd+1) ! caution nd instead ndp1 to be consistent...
780      real plcl(nloc), cpn(nloc,nd)
781
782c outputs:
783      integer iflag(nloc)
784      real cbmf(nloc) ! also an input
785
786c local variables:
787      integer i, k, icbmax
788      real dtpbl(nloc), dtmin(nloc), tvpplcl(nloc), tvaplcl(nloc)
789      real work(nloc)
790
791#include "cvthermo.h"
792#include "cvparam.h"
793
794c-------------------------------------------------------------------
795c Compute icbmax.
796c-------------------------------------------------------------------
797
798      icbmax=2
799      do 230 i=1,ncum
800       icbmax=max(icbmax,icb(i))
801 230  continue
802
803c=====================================================================
804c ---  CALCULATE CLOUD BASE MASS FLUX
805c=====================================================================
806c
807c tvpplcl = parcel temperature lifted adiabatically from level
808c           icb-1 to the LCL.
809c tvaplcl = virtual temperature at the LCL.
810c
811      do 610 i=1,ncum
812        dtpbl(i)=0.0
813        tvpplcl(i)=tvp(i,icb(i)-1)
814     &  -rrd*tvp(i,icb(i)-1)*(p(i,icb(i)-1)-plcl(i))
815     &  /(cpn(i,icb(i)-1)*p(i,icb(i)-1))
816        tvaplcl(i)=tv(i,icb(i))
817     &  +(tvp(i,icb(i))-tvp(i,icb(i)+1))*(plcl(i)-p(i,icb(i)))
818     &  /(p(i,icb(i))-p(i,icb(i)+1))
819 610  continue
820
821c-------------------------------------------------------------------
822c --- Interpolate difference between lifted parcel and
823c --- environmental temperatures to lifted condensation level
824c-------------------------------------------------------------------
825c
826c dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
827c
828      do 630 k=minorig,icbmax
829        do 620 i=1,ncum
830        if((k.ge.nk(i)).and.(k.le.(icb(i)-1)))then
831          dtpbl(i)=dtpbl(i)+(tvp(i,k)-tv(i,k))*dph(i,k)
832        endif
833 620    continue
834 630  continue
835      do 640 i=1,ncum
836        dtpbl(i)=dtpbl(i)/(ph(i,nk(i))-ph(i,icb(i)))
837        dtmin(i)=tvpplcl(i)-tvaplcl(i)+dtmax+dtpbl(i)
838 640  continue
839c
840c-------------------------------------------------------------------
841c --- Adjust cloud base mass flux
842c-------------------------------------------------------------------
843c
844      do 650 i=1,ncum
845       work(i)=cbmf(i)
846       cbmf(i)=max(0.0,(1.0-damp)*cbmf(i)+0.1*alpha*dtmin(i))
847       if((work(i).eq.0.0).and.(cbmf(i).eq.0.0))then
848         iflag(i)=3
849       endif
850 650  continue
851
852       return
853       end
854
855      SUBROUTINE cv_mixing(nloc,ncum,nd,icb,nk,inb,inb1
856     :                    ,ph,t,q,qs,u,v,h,lv,qnk
857     :                    ,hp,tv,tvp,ep,clw,cbmf
858     :                    ,m,ment,qent,uent,vent,nent,sij,elij)
859      implicit none
860
861#include "cvthermo.h"
862#include "cvparam.h"
863
864c inputs:
865      integer ncum, nd, nloc
866      integer icb(nloc), inb(nloc), inb1(nloc), nk(nloc)
867      real cbmf(nloc), qnk(nloc)
868      real ph(nloc,nd+1)
869      real t(nloc,nd), q(nloc,nd), qs(nloc,nd), lv(nloc,nd)
870      real u(nloc,nd), v(nloc,nd), h(nloc,nd), hp(nloc,nd)
871      real tv(nloc,nd), tvp(nloc,nd), ep(nloc,nd), clw(nloc,nd)
872
873c outputs:
874      integer nent(nloc,nd)
875      real m(nloc,nd), ment(nloc,nd,nd), qent(nloc,nd,nd)
876      real uent(nloc,nd,nd), vent(nloc,nd,nd)
877      real sij(nloc,nd,nd), elij(nloc,nd,nd)
878
879c local variables:
880      integer i, j, k, ij
881      integer num1, num2
882      real dbo, qti, bf2, anum, denom, dei, altem, cwat, stemp
883      real alt, qp1, smid, sjmin, sjmax, delp, delm
884      real work(nloc), asij(nloc), smin(nloc), scrit(nloc)
885      real bsum(nloc,nd)
886      logical lwork(nloc)
887
888c=====================================================================
889c --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
890c=====================================================================
891c
892        do 360 i=1,ncum*nlp
893          nent(i,1)=0
894          m(i,1)=0.0
895 360    continue
896c
897      do 400 k=1,nlp
898       do 390 j=1,nlp
899          do 385 i=1,ncum
900            qent(i,k,j)=q(i,j)
901            uent(i,k,j)=u(i,j)
902            vent(i,k,j)=v(i,j)
903            elij(i,k,j)=0.0
904            ment(i,k,j)=0.0
905            sij(i,k,j)=0.0
906 385      continue
907 390    continue
908 400  continue
909c
910c-------------------------------------------------------------------
911c --- Calculate rates of mixing,  m(i)
912c-------------------------------------------------------------------
913c
914      call zilch(work,ncum)
915c
916      do 670 j=minorig+1,nl
917        do 660 i=1,ncum
918          if((j.ge.(icb(i)+1)).and.(j.le.inb(i)))then
919             k=min(j,inb1(i))
920             dbo=abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1))
921     &       +entp*0.04*(ph(i,k)-ph(i,k+1))
922             work(i)=work(i)+dbo
923             m(i,j)=cbmf(i)*dbo
924          endif
925 660    continue
926 670  continue
927      do 690 k=minorig+1,nl
928        do 680 i=1,ncum
929          if((k.ge.(icb(i)+1)).and.(k.le.inb(i)))then
930            m(i,k)=m(i,k)/work(i)
931          endif
932 680    continue
933 690  continue
934c
935c
936c=====================================================================
937c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
938c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
939c --- FRACTION (sij)
940c=====================================================================
941c
942c
943       do 750 i=minorig+1,nl
944         do 710 j=minorig+1,nl
945           do 700 ij=1,ncum
946             if((i.ge.(icb(ij)+1)).and.(j.ge.icb(ij))
947     &         .and.(i.le.inb(ij)).and.(j.le.inb(ij)))then
948               qti=qnk(ij)-ep(ij,i)*clw(ij,i)
949               bf2=1.+lv(ij,j)*lv(ij,j)*qs(ij,j)
950     &         /(rrv*t(ij,j)*t(ij,j)*cpd)
951               anum=h(ij,j)-hp(ij,i)+(cpv-cpd)*t(ij,j)*(qti-q(ij,j))
952               denom=h(ij,i)-hp(ij,i)+(cpd-cpv)*(q(ij,i)-qti)*t(ij,j)
953               dei=denom
954               if(abs(dei).lt.0.01)dei=0.01
955               sij(ij,i,j)=anum/dei
956               sij(ij,i,i)=1.0
957               altem=sij(ij,i,j)*q(ij,i)+(1.-sij(ij,i,j))*qti-qs(ij,j)
958               altem=altem/bf2
959               cwat=clw(ij,j)*(1.-ep(ij,j))
960               stemp=sij(ij,i,j)
961               if((stemp.lt.0.0.or.stemp.gt.1.0.or.
962     1           altem.gt.cwat).and.j.gt.i)then
963                 anum=anum-lv(ij,j)*(qti-qs(ij,j)-cwat*bf2)
964                 denom=denom+lv(ij,j)*(q(ij,i)-qti)
965                 if(abs(denom).lt.0.01)denom=0.01
966                 sij(ij,i,j)=anum/denom
967                 altem=sij(ij,i,j)*q(ij,i)+(1.-sij(ij,i,j))*qti-qs(ij,j)
968                 altem=altem-(bf2-1.)*cwat
969               endif
970               if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then
971                 qent(ij,i,j)=sij(ij,i,j)*q(ij,i)
972     &                        +(1.-sij(ij,i,j))*qti
973                 uent(ij,i,j)=sij(ij,i,j)*u(ij,i)
974     &                        +(1.-sij(ij,i,j))*u(ij,nk(ij))
975                 vent(ij,i,j)=sij(ij,i,j)*v(ij,i)
976     &                        +(1.-sij(ij,i,j))*v(ij,nk(ij))
977                 elij(ij,i,j)=altem
978                 elij(ij,i,j)=max(0.0,elij(ij,i,j))
979                 ment(ij,i,j)=m(ij,i)/(1.-sij(ij,i,j))
980                 nent(ij,i)=nent(ij,i)+1
981               endif
982             sij(ij,i,j)=max(0.0,sij(ij,i,j))
983             sij(ij,i,j)=min(1.0,sij(ij,i,j))
984             endif
985  700      continue
986  710    continue
987c
988c   ***   If no air can entrain at level i assume that updraft detrains  ***
989c   ***   at that level and calculate detrained air flux and properties  ***
990c
991           do 740 ij=1,ncum
992             if((i.ge.(icb(ij)+1)).and.(i.le.inb(ij))
993     &       .and.(nent(ij,i).eq.0))then
994               ment(ij,i,i)=m(ij,i)
995               qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
996               uent(ij,i,i)=u(ij,nk(ij))
997               vent(ij,i,i)=v(ij,nk(ij))
998               elij(ij,i,i)=clw(ij,i)
999               sij(ij,i,i)=1.0
1000             endif
1001 740       continue
1002 750   continue
1003c
1004      do 770 i=1,ncum
1005        sij(i,inb(i),inb(i))=1.0
1006 770  continue
1007c
1008c=====================================================================
1009c   ---  NORMALIZE ENTRAINED AIR MASS FLUXES
1010c   ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
1011c=====================================================================
1012c
1013       call zilch(bsum,ncum*nlp)
1014       do 780 ij=1,ncum
1015         lwork(ij)=.false.
1016 780   continue
1017       do 789 i=minorig+1,nl
1018c
1019         num1=0
1020         do 779 ij=1,ncum
1021           if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))num1=num1+1
1022 779     continue
1023         if(num1.le.0)go to 789
1024c
1025           do 781 ij=1,ncum
1026             if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))then
1027                lwork(ij)=(nent(ij,i).ne.0)
1028                qp1=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
1029                anum=h(ij,i)-hp(ij,i)-lv(ij,i)*(qp1-qs(ij,i))
1030                denom=h(ij,i)-hp(ij,i)+lv(ij,i)*(q(ij,i)-qp1)
1031                if(abs(denom).lt.0.01)denom=0.01
1032                scrit(ij)=anum/denom
1033                alt=qp1-qs(ij,i)+scrit(ij)*(q(ij,i)-qp1)
1034                if(scrit(ij).lt.0.0.or.alt.lt.0.0)scrit(ij)=1.0
1035                asij(ij)=0.0
1036                smin(ij)=1.0
1037             endif
1038 781       continue
1039         do 783 j=minorig,nl
1040c
1041         num2=0
1042         do 778 ij=1,ncum
1043             if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
1044     &       .and.(j.ge.icb(ij)).and.(j.le.inb(ij))
1045     &       .and.lwork(ij))num2=num2+1
1046 778     continue
1047         if(num2.le.0)go to 783
1048c
1049           do 782 ij=1,ncum
1050             if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
1051     &       .and.(j.ge.icb(ij)).and.(j.le.inb(ij)).and.lwork(ij))then
1052                  if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then
1053                    if(j.gt.i)then
1054                      smid=min(sij(ij,i,j),scrit(ij))
1055                      sjmax=smid
1056                      sjmin=smid
1057                        if(smid.lt.smin(ij)
1058     &                  .and.sij(ij,i,j+1).lt.smid)then
1059                          smin(ij)=smid
1060                          sjmax=min(sij(ij,i,j+1),sij(ij,i,j),scrit(ij))
1061                          sjmin=max(sij(ij,i,j-1),sij(ij,i,j))
1062                          sjmin=min(sjmin,scrit(ij))
1063                        endif
1064                    else
1065                      sjmax=max(sij(ij,i,j+1),scrit(ij))
1066                      smid=max(sij(ij,i,j),scrit(ij))
1067                      sjmin=0.0
1068                      if(j.gt.1)sjmin=sij(ij,i,j-1)
1069                      sjmin=max(sjmin,scrit(ij))
1070                    endif
1071                    delp=abs(sjmax-smid)
1072                    delm=abs(sjmin-smid)
1073                    asij(ij)=asij(ij)+(delp+delm)
1074     &                           *(ph(ij,j)-ph(ij,j+1))
1075                    ment(ij,i,j)=ment(ij,i,j)*(delp+delm)
1076     &                           *(ph(ij,j)-ph(ij,j+1))
1077                  endif
1078              endif
1079  782    continue
1080  783    continue
1081            do 784 ij=1,ncum
1082            if((i.ge.icb(ij)+1).and.(i.le.inb(ij)).and.lwork(ij))then
1083               asij(ij)=max(1.0e-21,asij(ij))
1084               asij(ij)=1.0/asij(ij)
1085               bsum(ij,i)=0.0
1086            endif
1087 784        continue
1088            do 786 j=minorig,nl+1
1089              do 785 ij=1,ncum
1090                if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
1091     &          .and.(j.ge.icb(ij)).and.(j.le.inb(ij))
1092     &          .and.lwork(ij))then
1093                   ment(ij,i,j)=ment(ij,i,j)*asij(ij)
1094                   bsum(ij,i)=bsum(ij,i)+ment(ij,i,j)
1095                endif
1096 785     continue
1097 786     continue
1098             do 787 ij=1,ncum
1099               if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
1100     &         .and.(bsum(ij,i).lt.1.0e-18).and.lwork(ij))then
1101                 nent(ij,i)=0
1102                 ment(ij,i,i)=m(ij,i)
1103                 qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
1104                 uent(ij,i,i)=u(ij,nk(ij))
1105                 vent(ij,i,i)=v(ij,nk(ij))
1106                 elij(ij,i,i)=clw(ij,i)
1107                 sij(ij,i,i)=1.0
1108               endif
1109  787        continue
1110  789  continue
1111c
1112       return
1113       end
1114
1115      SUBROUTINE cv_unsat(nloc,ncum,nd,inb,t,q,qs,gz,u,v,p,ph
1116     :                  ,h,lv,ep,sigp,clw,m,ment,elij
1117     :                  ,iflag,mp,qp,up,vp,wt,water,evap)
1118      implicit none
1119
1120
1121#include "cvthermo.h"
1122#include "cvparam.h"
1123
1124c inputs:
1125      integer ncum, nd, nloc
1126      integer inb(nloc)
1127      real t(nloc,nd), q(nloc,nd), qs(nloc,nd)
1128      real gz(nloc,nd), u(nloc,nd), v(nloc,nd)
1129      real p(nloc,nd), ph(nloc,nd+1), h(nloc,nd)
1130      real lv(nloc,nd), ep(nloc,nd), sigp(nloc,nd), clw(nloc,nd)
1131      real m(nloc,nd), ment(nloc,nd,nd), elij(nloc,nd,nd)
1132
1133c outputs:
1134      integer iflag(nloc) ! also an input
1135      real mp(nloc,nd), qp(nloc,nd), up(nloc,nd), vp(nloc,nd)
1136      real water(nloc,nd), evap(nloc,nd), wt(nloc,nd)
1137
1138c local variables:
1139      integer i,j,k,ij,num1
1140      integer jtt(nloc)
1141      real awat, coeff, qsm, afac, sigt, b6, c6, revap
1142      real dhdp, fac, qstm, rat
1143      real wdtrain(nloc)
1144      logical lwork(nloc)
1145
1146c=====================================================================
1147c --- PRECIPITATING DOWNDRAFT CALCULATION
1148c=====================================================================
1149c
1150c Initializations:
1151c
1152         do i = 1, ncum
1153         do k = 1, nl+1
1154          wt(i,k) = omtsnow
1155          mp(i,k) = 0.0
1156          evap(i,k) = 0.0
1157          water(i,k) = 0.0
1158         enddo
1159         enddo
1160
1161         do 420 i=1,ncum
1162          qp(i,1)=q(i,1)
1163          up(i,1)=u(i,1)
1164          vp(i,1)=v(i,1)
1165 420     continue
1166
1167         do 440 k=2,nl+1
1168         do 430 i=1,ncum
1169          qp(i,k)=q(i,k-1)
1170          up(i,k)=u(i,k-1)
1171          vp(i,k)=v(i,k-1)
1172 430     continue
1173 440     continue
1174
1175
1176c   ***  Check whether ep(inb)=0, if so, skip precipitating    ***
1177c   ***             downdraft calculation                      ***
1178c
1179c
1180c   ***  Integrate liquid water equation to find condensed water   ***
1181c   ***                and condensed water flux                    ***
1182c
1183c
1184      do 890 i=1,ncum
1185        jtt(i)=2
1186        if(ep(i,inb(i)).le.0.0001)iflag(i)=2
1187        if(iflag(i).eq.0)then
1188          lwork(i)=.true.
1189        else
1190          lwork(i)=.false.
1191        endif
1192 890  continue
1193c
1194c    ***                    Begin downdraft loop                    ***
1195c
1196c
1197        call zilch(wdtrain,ncum)
1198        do 899 i=nl+1,1,-1
1199c
1200          num1=0
1201          do 879 ij=1,ncum
1202            if((i.le.inb(ij)).and.lwork(ij))num1=num1+1
1203 879      continue
1204          if(num1.le.0)go to 899
1205c
1206c
1207c    ***        Calculate detrained precipitation             ***
1208c
1209          do 891 ij=1,ncum
1210            if((i.le.inb(ij)).and.(lwork(ij)))then
1211            wdtrain(ij)=g*ep(ij,i)*m(ij,i)*clw(ij,i)
1212            endif
1213 891      continue
1214c
1215          if(i.gt.1)then
1216            do 893 j=1,i-1
1217              do 892 ij=1,ncum
1218                if((i.le.inb(ij)).and.(lwork(ij)))then
1219                  awat=elij(ij,j,i)-(1.-ep(ij,i))*clw(ij,i)
1220                  awat=max(0.0,awat)
1221                  wdtrain(ij)=wdtrain(ij)+g*awat*ment(ij,j,i)
1222                endif
1223 892          continue
1224 893      continue
1225          endif
1226c
1227c    ***    Find rain water and evaporation using provisional   ***
1228c    ***              estimates of qp(i)and qp(i-1)             ***
1229c
1230c
1231c  ***  Value of terminal velocity and coeffecient of evaporation for snow   ***
1232c
1233          do 894 ij=1,ncum
1234            if((i.le.inb(ij)).and.(lwork(ij)))then
1235            coeff=coeffs
1236            wt(ij,i)=omtsnow
1237c
1238c  ***  Value of terminal velocity and coeffecient of evaporation for rain   ***
1239c
1240            if(t(ij,i).gt.273.0)then
1241              coeff=coeffr
1242              wt(ij,i)=omtrain
1243            endif
1244            qsm=0.5*(q(ij,i)+qp(ij,i+1))
1245            afac=coeff*ph(ij,i)*(qs(ij,i)-qsm)
1246     &       /(1.0e4+2.0e3*ph(ij,i)*qs(ij,i))
1247            afac=max(afac,0.0)
1248            sigt=sigp(ij,i)
1249            sigt=max(0.0,sigt)
1250            sigt=min(1.0,sigt)
1251            b6=100.*(ph(ij,i)-ph(ij,i+1))*sigt*afac/wt(ij,i)
1252            c6=(water(ij,i+1)*wt(ij,i+1)+wdtrain(ij)/sigd)/wt(ij,i)
1253            revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
1254            evap(ij,i)=sigt*afac*revap
1255            water(ij,i)=revap*revap
1256c
1257c    ***  Calculate precipitating downdraft mass flux under     ***
1258c    ***              hydrostatic approximation                 ***
1259c
1260            if(i.gt.1)then
1261              dhdp=(h(ij,i)-h(ij,i-1))/(p(ij,i-1)-p(ij,i))
1262              dhdp=max(dhdp,10.0)
1263              mp(ij,i)=100.*ginv*lv(ij,i)*sigd*evap(ij,i)/dhdp
1264              mp(ij,i)=max(mp(ij,i),0.0)
1265c
1266c   ***   Add small amount of inertia to downdraft              ***
1267c
1268              fac=20.0/(ph(ij,i-1)-ph(ij,i))
1269              mp(ij,i)=(fac*mp(ij,i+1)+mp(ij,i))/(1.+fac)
1270c
1271c    ***      Force mp to decrease linearly to zero                 ***
1272c    ***      between about 950 mb and the surface                  ***
1273c
1274              if(p(ij,i).gt.(0.949*p(ij,1)))then
1275                 jtt(ij)=max(jtt(ij),i)
1276                 mp(ij,i)=mp(ij,jtt(ij))*(p(ij,1)-p(ij,i))
1277     &           /(p(ij,1)-p(ij,jtt(ij)))
1278              endif
1279            endif
1280c
1281c    ***       Find mixing ratio of precipitating downdraft     ***
1282c
1283            if(i.ne.inb(ij))then
1284              if(i.eq.1)then
1285                qstm=qs(ij,1)
1286              else
1287                qstm=qs(ij,i-1)
1288              endif
1289              if(mp(ij,i).gt.mp(ij,i+1))then
1290                 rat=mp(ij,i+1)/mp(ij,i)
1291                 qp(ij,i)=qp(ij,i+1)*rat+q(ij,i)*(1.0-rat)+100.*ginv*
1292     &             sigd*(ph(ij,i)-ph(ij,i+1))*(evap(ij,i)/mp(ij,i))
1293                 up(ij,i)=up(ij,i+1)*rat+u(ij,i)*(1.-rat)
1294                 vp(ij,i)=vp(ij,i+1)*rat+v(ij,i)*(1.-rat)
1295               else
1296                 if(mp(ij,i+1).gt.0.0)then
1297                   qp(ij,i)=(gz(ij,i+1)-gz(ij,i)
1298     &               +qp(ij,i+1)*(lv(ij,i+1)+t(ij,i+1)
1299     &               *(cl-cpd))+cpd*(t(ij,i+1)-t(ij,i)))
1300     &               /(lv(ij,i)+t(ij,i)*(cl-cpd))
1301                   up(ij,i)=up(ij,i+1)
1302                   vp(ij,i)=vp(ij,i+1)
1303                 endif
1304              endif
1305              qp(ij,i)=min(qp(ij,i),qstm)
1306              qp(ij,i)=max(qp(ij,i),0.0)
1307            endif
1308            endif
1309 894      continue
1310 899    continue
1311c
1312        return
1313        end
1314
1315      SUBROUTINE cv_yield(nloc,ncum,nd,nk,icb,inb,delt
1316     :             ,t,q,u,v,gz,p,ph,h,hp,lv,cpn
1317     :             ,ep,clw,frac,m,mp,qp,up,vp
1318     :             ,wt,water,evap
1319     :             ,ment,qent,uent,vent,nent,elij
1320     :             ,tv,tvp
1321     o             ,iflag,wd,qprime,tprime
1322     o             ,precip,cbmf,ft,fq,fu,fv,Ma,qcondc)
1323      implicit none
1324
1325#include "cvthermo.h"
1326#include "cvparam.h"
1327
1328c inputs
1329      integer ncum, nd, nloc
1330      integer nk(nloc), icb(nloc), inb(nloc)
1331      integer nent(nloc,nd)
1332      real delt
1333      real t(nloc,nd), q(nloc,nd), u(nloc,nd), v(nloc,nd)
1334      real gz(nloc,nd)
1335      real p(nloc,nd), ph(nloc,nd+1), h(nloc,nd)
1336      real hp(nloc,nd), lv(nloc,nd)
1337      real cpn(nloc,nd), ep(nloc,nd), clw(nloc,nd), frac(nloc)
1338      real m(nloc,nd), mp(nloc,nd), qp(nloc,nd)
1339      real up(nloc,nd), vp(nloc,nd)
1340      real wt(nloc,nd), water(nloc,nd), evap(nloc,nd)
1341      real ment(nloc,nd,nd), qent(nloc,nd,nd), elij(nloc,nd,nd)
1342      real uent(nloc,nd,nd), vent(nloc,nd,nd)
1343      real tv(nloc,nd), tvp(nloc,nd)
1344
1345c outputs
1346      integer iflag(nloc)  ! also an input
1347      real cbmf(nloc)      ! also an input
1348      real wd(nloc), tprime(nloc), qprime(nloc)
1349      real precip(nloc)
1350      real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
1351      real Ma(nloc,nd)
1352      real qcondc(nloc,nd)
1353
1354c local variables
1355      integer i,j,ij,k,num1
1356      real dpinv,cpinv,awat,fqold,ftold,fuold,fvold,delti
1357      real work(nloc), am(nloc),amp1(nloc),ad(nloc)
1358      real ents(nloc), uav(nloc),vav(nloc),lvcp(nloc,nd)
1359      real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd) ! cld
1360      real siga(nloc,nd), ax(nloc,nd), mac(nloc,nd)     ! cld
1361
1362 
1363c -- initializations:
1364
1365      delti = 1.0/delt
1366
1367      do 160 i=1,ncum
1368      precip(i)=0.0
1369      wd(i)=0.0
1370      tprime(i)=0.0
1371      qprime(i)=0.0
1372       do 170 k=1,nl+1
1373        ft(i,k)=0.0
1374        fu(i,k)=0.0
1375        fv(i,k)=0.0
1376        fq(i,k)=0.0
1377        lvcp(i,k)=lv(i,k)/cpn(i,k)
1378        qcondc(i,k)=0.0              ! cld
1379        qcond(i,k)=0.0               ! cld
1380        nqcond(i,k)=0.0              ! cld
1381 170   continue
1382 160  continue
1383
1384c
1385c   ***  Calculate surface precipitation in mm/day     ***
1386c
1387        do 1190 i=1,ncum
1388          if(iflag(i).le.1)then
1389cc            precip(i)=precip(i)+wt(i,1)*sigd*water(i,1)*3600.*24000.
1390cc     &                /(rowl*g)
1391cc            precip(i)=precip(i)*delt/86400.
1392            precip(i) = wt(i,1)*sigd*water(i,1)*86400/g
1393          endif
1394 1190   continue
1395c
1396c
1397c   ***  Calculate downdraft velocity scale and surface temperature and  ***
1398c   ***                    water vapor fluctuations                      ***
1399c
1400      do i=1,ncum
1401       wd(i)=betad*abs(mp(i,icb(i)))*0.01*rrd*t(i,icb(i))
1402     :           /(sigd*p(i,icb(i)))
1403       qprime(i)=0.5*(qp(i,1)-q(i,1))
1404       tprime(i)=lv0*qprime(i)/cpd
1405      enddo
1406c
1407c   ***  Calculate tendencies of lowest level potential temperature  ***
1408c   ***                      and mixing ratio                        ***
1409c
1410        do 1200 i=1,ncum
1411          work(i)=0.01/(ph(i,1)-ph(i,2))
1412          am(i)=0.0
1413 1200   continue
1414        do 1220 k=2,nl
1415          do 1210 i=1,ncum
1416            if((nk(i).eq.1).and.(k.le.inb(i)).and.(nk(i).eq.1))then
1417              am(i)=am(i)+m(i,k)
1418            endif
1419 1210     continue
1420 1220   continue
1421        do 1240 i=1,ncum
1422          if((g*work(i)*am(i)).ge.delti)iflag(i)=1
1423          ft(i,1)=ft(i,1)+g*work(i)*am(i)*(t(i,2)-t(i,1)
1424     &    +(gz(i,2)-gz(i,1))/cpn(i,1))
1425          ft(i,1)=ft(i,1)-lvcp(i,1)*sigd*evap(i,1)
1426          ft(i,1)=ft(i,1)+sigd*wt(i,2)*(cl-cpd)*water(i,2)*(t(i,2)
1427     &     -t(i,1))*work(i)/cpn(i,1)
1428          fq(i,1)=fq(i,1)+g*mp(i,2)*(qp(i,2)-q(i,1))*
1429     &    work(i)+sigd*evap(i,1)
1430          fq(i,1)=fq(i,1)+g*am(i)*(q(i,2)-q(i,1))*work(i)
1431          fu(i,1)=fu(i,1)+g*work(i)*(mp(i,2)*(up(i,2)-u(i,1))
1432     &    +am(i)*(u(i,2)-u(i,1)))
1433          fv(i,1)=fv(i,1)+g*work(i)*(mp(i,2)*(vp(i,2)-v(i,1))
1434     &    +am(i)*(v(i,2)-v(i,1)))
1435 1240   continue
1436        do 1260 j=2,nl
1437           do 1250 i=1,ncum
1438             if(j.le.inb(i))then
1439               fq(i,1)=fq(i,1)
1440     &                 +g*work(i)*ment(i,j,1)*(qent(i,j,1)-q(i,1))
1441               fu(i,1)=fu(i,1)
1442     &                 +g*work(i)*ment(i,j,1)*(uent(i,j,1)-u(i,1))
1443               fv(i,1)=fv(i,1)
1444     &                 +g*work(i)*ment(i,j,1)*(vent(i,j,1)-v(i,1))
1445             endif
1446 1250      continue
1447 1260   continue
1448c
1449c   ***  Calculate tendencies of potential temperature and mixing ratio  ***
1450c   ***               at levels above the lowest level                   ***
1451c
1452c   ***  First find the net saturated updraft and downdraft mass fluxes  ***
1453c   ***                      through each level                          ***
1454c
1455        do 1500 i=2,nl+1
1456c
1457          num1=0
1458          do 1265 ij=1,ncum
1459            if(i.le.inb(ij))num1=num1+1
1460 1265     continue
1461          if(num1.le.0)go to 1500
1462c
1463          call zilch(amp1,ncum)
1464          call zilch(ad,ncum)
1465c
1466          do 1280 k=i+1,nl+1
1467            do 1270 ij=1,ncum
1468              if((i.ge.nk(ij)).and.(i.le.inb(ij))
1469     &            .and.(k.le.(inb(ij)+1)))then
1470                amp1(ij)=amp1(ij)+m(ij,k)
1471              endif
1472 1270         continue
1473 1280     continue
1474c
1475          do 1310 k=1,i
1476            do 1300 j=i+1,nl+1
1477               do 1290 ij=1,ncum
1478                 if((j.le.(inb(ij)+1)).and.(i.le.inb(ij)))then
1479                   amp1(ij)=amp1(ij)+ment(ij,k,j)
1480                 endif
1481 1290          continue
1482 1300       continue
1483 1310     continue
1484          do 1340 k=1,i-1
1485            do 1330 j=i,nl+1
1486              do 1320 ij=1,ncum
1487                if((i.le.inb(ij)).and.(j.le.inb(ij)))then
1488                   ad(ij)=ad(ij)+ment(ij,j,k)
1489                endif
1490 1320         continue
1491 1330       continue
1492 1340     continue
1493c
1494          do 1350 ij=1,ncum
1495          if(i.le.inb(ij))then
1496            dpinv=0.01/(ph(ij,i)-ph(ij,i+1))
1497            cpinv=1.0/cpn(ij,i)
1498c
1499            ft(ij,i)=ft(ij,i)
1500     &       +g*dpinv*(amp1(ij)*(t(ij,i+1)-t(ij,i)
1501     &       +(gz(ij,i+1)-gz(ij,i))*cpinv)
1502     &       -ad(ij)*(t(ij,i)-t(ij,i-1)+(gz(ij,i)-gz(ij,i-1))*cpinv))
1503     &       -sigd*lvcp(ij,i)*evap(ij,i)
1504            ft(ij,i)=ft(ij,i)+g*dpinv*ment(ij,i,i)*(hp(ij,i)-h(ij,i)+
1505     &        t(ij,i)*(cpv-cpd)*(q(ij,i)-qent(ij,i,i)))*cpinv
1506            ft(ij,i)=ft(ij,i)+sigd*wt(ij,i+1)*(cl-cpd)*water(ij,i+1)*
1507     &        (t(ij,i+1)-t(ij,i))*dpinv*cpinv
1508            fq(ij,i)=fq(ij,i)+g*dpinv*(amp1(ij)*(q(ij,i+1)-q(ij,i))-
1509     &        ad(ij)*(q(ij,i)-q(ij,i-1)))
1510            fu(ij,i)=fu(ij,i)+g*dpinv*(amp1(ij)*(u(ij,i+1)-u(ij,i))-
1511     &        ad(ij)*(u(ij,i)-u(ij,i-1)))
1512            fv(ij,i)=fv(ij,i)+g*dpinv*(amp1(ij)*(v(ij,i+1)-v(ij,i))-
1513     &        ad(ij)*(v(ij,i)-v(ij,i-1)))
1514         endif
1515 1350    continue
1516         do 1370 k=1,i-1
1517           do 1360 ij=1,ncum
1518             if(i.le.inb(ij))then
1519               awat=elij(ij,k,i)-(1.-ep(ij,i))*clw(ij,i)
1520               awat=max(awat,0.0)
1521               fq(ij,i)=fq(ij,i)
1522     &         +g*dpinv*ment(ij,k,i)*(qent(ij,k,i)-awat-q(ij,i))
1523               fu(ij,i)=fu(ij,i)
1524     &         +g*dpinv*ment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
1525               fv(ij,i)=fv(ij,i)
1526     &         +g*dpinv*ment(ij,k,i)*(vent(ij,k,i)-v(ij,i))
1527c (saturated updrafts resulting from mixing)               ! cld
1528               qcond(ij,i)=qcond(ij,i)+(elij(ij,k,i)-awat) ! cld
1529               nqcond(ij,i)=nqcond(ij,i)+1.                ! cld
1530             endif
1531 1360      continue
1532 1370    continue
1533         do 1390 k=i,nl+1
1534           do 1380 ij=1,ncum
1535             if((i.le.inb(ij)).and.(k.le.inb(ij)))then
1536               fq(ij,i)=fq(ij,i)
1537     &                  +g*dpinv*ment(ij,k,i)*(qent(ij,k,i)-q(ij,i))
1538               fu(ij,i)=fu(ij,i)
1539     &                  +g*dpinv*ment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
1540               fv(ij,i)=fv(ij,i)
1541     &                  +g*dpinv*ment(ij,k,i)*(vent(ij,k,i)-v(ij,i))
1542             endif
1543 1380      continue
1544 1390    continue
1545          do 1400 ij=1,ncum
1546           if(i.le.inb(ij))then
1547             fq(ij,i)=fq(ij,i)
1548     &                +sigd*evap(ij,i)+g*(mp(ij,i+1)*
1549     &                (qp(ij,i+1)-q(ij,i))
1550     &                -mp(ij,i)*(qp(ij,i)-q(ij,i-1)))*dpinv
1551             fu(ij,i)=fu(ij,i)
1552     &                +g*(mp(ij,i+1)*(up(ij,i+1)-u(ij,i))-mp(ij,i)*
1553     &                (up(ij,i)-u(ij,i-1)))*dpinv
1554             fv(ij,i)=fv(ij,i)
1555     &               +g*(mp(ij,i+1)*(vp(ij,i+1)-v(ij,i))-mp(ij,i)*
1556     &               (vp(ij,i)-v(ij,i-1)))*dpinv
1557C (saturated downdrafts resulting from mixing)               ! cld
1558            do k=i+1,inb(ij)                                 ! cld
1559             qcond(ij,i)=qcond(ij,i)+elij(ij,k,i)            ! cld
1560             nqcond(ij,i)=nqcond(ij,i)+1.                    ! cld
1561            enddo                                            ! cld
1562C (particular case: no detraining level is found)            ! cld
1563            if (nent(ij,i).eq.0) then                        ! cld
1564             qcond(ij,i)=qcond(ij,i)+(1.-ep(ij,i))*clw(ij,i) ! cld
1565             nqcond(ij,i)=nqcond(ij,i)+1.                    ! cld
1566            endif                                            ! cld
1567            if (nqcond(ij,i).ne.0.) then                     ! cld
1568             qcond(ij,i)=qcond(ij,i)/nqcond(ij,i)            ! cld
1569            endif                                            ! cld
1570           endif
1571 1400     continue
1572 1500   continue
1573c
1574c   *** Adjust tendencies at top of convection layer to reflect  ***
1575c   ***       actual position of the level zero cape             ***
1576c
1577        do 503 ij=1,ncum
1578        fqold=fq(ij,inb(ij))
1579        fq(ij,inb(ij))=fq(ij,inb(ij))*(1.-frac(ij))
1580        fq(ij,inb(ij)-1)=fq(ij,inb(ij)-1)
1581     &   +frac(ij)*fqold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
1582     1   (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*lv(ij,inb(ij))
1583     &   /lv(ij,inb(ij)-1)
1584        ftold=ft(ij,inb(ij))
1585        ft(ij,inb(ij))=ft(ij,inb(ij))*(1.-frac(ij))
1586        ft(ij,inb(ij)-1)=ft(ij,inb(ij)-1)
1587     &   +frac(ij)*ftold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
1588     1   (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*cpn(ij,inb(ij))
1589     &   /cpn(ij,inb(ij)-1)
1590        fuold=fu(ij,inb(ij))
1591        fu(ij,inb(ij))=fu(ij,inb(ij))*(1.-frac(ij))
1592        fu(ij,inb(ij)-1)=fu(ij,inb(ij)-1)
1593     &   +frac(ij)*fuold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
1594     1   (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
1595        fvold=fv(ij,inb(ij))
1596        fv(ij,inb(ij))=fv(ij,inb(ij))*(1.-frac(ij))
1597        fv(ij,inb(ij)-1)=fv(ij,inb(ij)-1)
1598     &  +frac(ij)*fvold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
1599     1   (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
1600 503    continue
1601c
1602c   ***   Very slightly adjust tendencies to force exact   ***
1603c   ***     enthalpy, momentum and tracer conservation     ***
1604c
1605        do 682 ij=1,ncum
1606        ents(ij)=0.0
1607        uav(ij)=0.0
1608        vav(ij)=0.0
1609        do 681 i=1,inb(ij)
1610         ents(ij)=ents(ij)
1611     &  +(cpn(ij,i)*ft(ij,i)+lv(ij,i)*fq(ij,i))*(ph(ij,i)-ph(ij,i+1))   
1612         uav(ij)=uav(ij)+fu(ij,i)*(ph(ij,i)-ph(ij,i+1))
1613         vav(ij)=vav(ij)+fv(ij,i)*(ph(ij,i)-ph(ij,i+1))
1614  681   continue
1615  682   continue
1616        do 683 ij=1,ncum
1617        ents(ij)=ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1618        uav(ij)=uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1619        vav(ij)=vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1620 683    continue
1621        do 642 ij=1,ncum
1622        do 641 i=1,inb(ij)
1623         ft(ij,i)=ft(ij,i)-ents(ij)/cpn(ij,i)
1624         fu(ij,i)=(1.-cu)*(fu(ij,i)-uav(ij))
1625         fv(ij,i)=(1.-cu)*(fv(ij,i)-vav(ij))
1626  641   continue
1627 642    continue
1628c
1629        do 1810 k=1,nl+1
1630          do 1800 i=1,ncum
1631            if((q(i,k)+delt*fq(i,k)).lt.0.0)iflag(i)=10
1632 1800     continue
1633 1810   continue
1634c
1635c
1636        do 1900 i=1,ncum
1637          if(iflag(i).gt.2)then
1638          precip(i)=0.0
1639          cbmf(i)=0.0
1640          endif
1641 1900   continue
1642        do 1920 k=1,nl
1643         do 1910 i=1,ncum
1644           if(iflag(i).gt.2)then
1645             ft(i,k)=0.0
1646             fq(i,k)=0.0
1647             fu(i,k)=0.0
1648             fv(i,k)=0.0
1649             qcondc(i,k)=0.0                               ! cld
1650           endif
1651 1910    continue
1652 1920   continue
1653
1654        do k=1,nl+1
1655        do i=1,ncum
1656          Ma(i,k) = 0.
1657        enddo
1658        enddo
1659        do k=nl,1,-1
1660        do i=1,ncum
1661          Ma(i,k) = Ma(i,k+1)+m(i,k)
1662        enddo
1663        enddo
1664
1665c
1666c   *** diagnose the in-cloud mixing ratio   ***            ! cld
1667c   ***           of condensed water         ***            ! cld
1668c                                                           ! cld
1669      DO ij=1,ncum                                          ! cld   
1670       do i=1,nd                                            ! cld
1671        mac(ij,i)=0.0                                       ! cld   
1672        wa(ij,i)=0.0                                        ! cld
1673        siga(ij,i)=0.0                                      ! cld
1674       enddo                                                ! cld
1675       do i=nk(ij),inb(ij)                                  ! cld
1676       do k=i+1,inb(ij)+1                                   ! cld
1677        mac(ij,i)=mac(ij,i)+m(ij,k)                         ! cld
1678       enddo                                                ! cld
1679       enddo                                                ! cld
1680       do i=icb(ij),inb(ij)-1                               ! cld
1681        ax(ij,i)=0.                                         ! cld
1682        do j=icb(ij),i                                      ! cld
1683         ax(ij,i)=ax(ij,i)+rrd*(tvp(ij,j)-tv(ij,j))         ! cld   
1684     :       *(ph(ij,j)-ph(ij,j+1))/p(ij,j)                 ! cld   
1685        enddo                                               ! cld
1686        if (ax(ij,i).gt.0.0) then                           ! cld   
1687         wa(ij,i)=sqrt(2.*ax(ij,i))                         ! cld
1688        endif                                               ! cld
1689       enddo                                                ! cld
1690       do i=1,nl                                            ! cld
1691        if (wa(ij,i).gt.0.0)                                ! cld
1692     :    siga(ij,i)=mac(ij,i)/wa(ij,i)                     ! cld   
1693     :        *rrd*tvp(ij,i)/p(ij,i)/100./delta             ! cld   
1694        siga(ij,i) = min(siga(ij,i),1.0)                    ! cld
1695        qcondc(ij,i)=siga(ij,i)*clw(ij,i)*(1.-ep(ij,i))     ! cld   
1696     :          + (1.-siga(ij,i))*qcond(ij,i)               ! cld   
1697       enddo                                                ! cld
1698      ENDDO                                                 ! cld   
1699
1700        return
1701        end
1702
1703      SUBROUTINE cv_uncompress(nloc,len,ncum,nd,idcum
1704     :         ,iflag
1705     :         ,precip,cbmf
1706     :         ,ft,fq,fu,fv
1707     :         ,Ma,qcondc           
1708     :         ,iflag1
1709     :         ,precip1,cbmf1
1710     :         ,ft1,fq1,fu1,fv1
1711     :         ,Ma1,qcondc1           
1712     :                               )
1713      implicit none
1714
1715#include "cvparam.h"
1716
1717c inputs:
1718      integer len, ncum, nd, nloc
1719      integer idcum(nloc)
1720      integer iflag(nloc)
1721      real precip(nloc), cbmf(nloc)
1722      real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
1723      real Ma(nloc,nd)
1724      real qcondc(nloc,nd) !cld
1725
1726c outputs:
1727      integer iflag1(len)
1728      real precip1(len), cbmf1(len)
1729      real ft1(len,nd), fq1(len,nd), fu1(len,nd), fv1(len,nd)
1730      real Ma1(len,nd)
1731      real qcondc1(len,nd) !cld
1732
1733c local variables:
1734      integer i,k
1735
1736        do 2000 i=1,ncum
1737         precip1(idcum(i))=precip(i)
1738         cbmf1(idcum(i))=cbmf(i)
1739         iflag1(idcum(i))=iflag(i)
1740 2000   continue
1741
1742        do 2020 k=1,nl
1743          do 2010 i=1,ncum
1744            ft1(idcum(i),k)=ft(i,k)
1745            fq1(idcum(i),k)=fq(i,k)
1746            fu1(idcum(i),k)=fu(i,k)
1747            fv1(idcum(i),k)=fv(i,k)
1748            Ma1(idcum(i),k)=Ma(i,k)
1749            qcondc1(idcum(i),k)=qcondc(i,k)
1750 2010     continue
1751 2020   continue
1752
1753        return
1754        end
1755
Note: See TracBrowser for help on using the repository browser.