source: LMDZ.3.3/trunk/libf/phylmd/cv_routines.F @ 416

Last change on this file since 416 was 416, checked in by lmdzadmin, 22 years ago

Inclusion initiale

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