/[lmdze]/trunk/phylmd/CV_routines/cv_param.f
ViewVC logotype

Annotation of /trunk/phylmd/CV_routines/cv_param.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 12 - (hide annotations)
Mon Jul 21 16:05:07 2008 UTC (15 years, 10 months ago) by guez
Original Path: trunk/libf/phylmd/cv_routines.f
File size: 56400 byte(s) 
-- Minor modification of input/output:

Created procedure "read_logic". Variables of module "logic" are read
by "read_logic" instead of "conf_gcm". Variable "offline" of module
"conf_gcm" is read from namelist instead of "*.def".

Deleted arguments "dtime", "co2_ppm_etat0", "solaire_etat0",
"tabcntr0" and local variables "radpas", "tab_cntrl" of
"phyetat0". "phyetat0" does not read "controle" in "startphy.nc" any
longer. "phyetat0" now reads global attribute "itau_phy" from
"startphy.nc". "phyredem" does not create variable "controle" in
"startphy.nc" any longer. "phyredem" now writes global attribute
"itau_phy" of "startphy.nc". Deleted argument "tabcntr0" of
"printflag". Removed diagnostic messages written by "printflag" for
comparison of the variable "controle" of "startphy.nc" and the
variables read from "*.def" or namelist input.

-- Removing unwanted functionality:

Removed variable "lunout" from module "iniprint", replaced everywhere
by standard output.

Removed case "ocean == 'couple'" in "clmain", "interfsurf_hq" and
"physiq". Removed procedure "interfoce_cpl".

-- Should not change anything at run time:

Automated creation of graphs in documentation. More documentation on
input files.

Converted Fortran files to free format: "phyredem.f90", "printflag.f90".

Split module "clesphy" into "clesphys" and "clesphys2".

Removed variables "conser", "leapf", "forward", "apphys", "apdiss" and
"statcl" from module "logic". Added arguments "conser" to "advect",
"leapf" to "integrd". Added local variables "forward", "leapf",
"apphys", "conser", "apdiss" in "leapfrog".

Added intent attributes.

Deleted arguments "dtime" of "phyredem", "pdtime" of "flxdtdq", "sh"
of "phytrac", "dt" of "yamada".

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