source: LMDZ5/trunk/libf/phylmd/cv_routines.F @ 1434

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

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

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

Local machine: bench configuration, 32x24x11, gfortran

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

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

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

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

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

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

En local: 32x24x11, config bench/gfortran

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

Sur brodie, 96x95x39 sur 4 proc:

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

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