My Project
 All Classes Files Functions Variables Macros
cv_routines.F
Go to the documentation of this file.
1 !
2 ! $Id: cv_routines.F 1403 2010-07-01 09:02:53Z fairhead $
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  CHARACTER (LEN=20) :: modname='cv_routines'
41  CHARACTER (LEN=80) :: abort_message
42 
43 c noff: integer limit for convection (nd-noff)
44 c 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
63 c
64  cu=0.70
65 c
66  betad=10.0
67 c
68  damp=0.1
69  alpha=0.2
70 c
71  delta=0.01 ! cld
72 c
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 
84 c inputs:
85  integer len, nd, ndp1
86  real t(len,nd), q(len,nd), p(len,nd), ph(len,ndp1)
87 
88 c outputs:
89  real lv(len,nd), cpn(len,nd), tv(len,nd)
90  real gz(len,nd), h(len,nd), hm(len,nd)
91 
92 c 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
108 c
109 c gz = phi at the full levels (same as p).
110 c
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
120 c
121 c h = phi + cpT (dry static energy).
122 c hm = phi + cp(T-Tbase)+Lq
123 c
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 
138 C================================================================
139 C Purpose: CONVECTIVE FEED
140 C================================================================
141 
142 #include "cvparam.h"
143 
144 c 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 
149 c outputs:
150  integer iflag(len), nk(len), icb(len), icbmax
151  real tnk(len), qnk(len), gznk(len), plcl(len)
152 
153 c 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
185 c
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))
223 c
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
237 c
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
244 c
245  do 300 i=1,len
246  if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
247  300 continue
248 c
249 c Compute icbmax.
250 c
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 
266 c 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 
272 c outputs:
273  real tp(len,nd), tvp(len,nd), clw(len,nd)
274 
275 c 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
295 c
296 c *** Calculate certain parcel quantities, including static energy ***
297 c
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
303 c
304 c *** Calculate lifted parcel quantities below cloud base ***
305 c
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
312 c
313 c *** Find lifted parcel quantities above cloud base ***
314 c
315  do 360 i=1,len
316  tg=ticb(i)
317  qg=qs(i,icb(i))
318  alv=lv0-clmcpv*(ticb(i)-t0)
319 c
320 c First iteration.
321 c
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))
335 c
336 c Second iteration.
337 c
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))
351 c
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
360 c
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
366 c
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 
381 c inputs:
382  integer len, nd, icb(len)
383  real cbmf(len), tv(len,nd), tvp(len,nd)
384 
385 c outputs:
386  integer iflag(len) ! also an input
387 
388 c 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 
413 c 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 
422 c 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 
431 c 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 
499 C---------------------------------------------------------------------
500 C Purpose:
501 C FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
502 C &
503 C COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
504 C FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
505 C &
506 C FIND THE LEVEL OF NEUTRAL BUOYANCY
507 C---------------------------------------------------------------------
508 
509 #include "cvthermo.h"
510 #include "cvparam.h"
511 
512 c 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 
520 c 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 
526 c 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 !=====================================================================
547 c
548 c --- The procedure is to solve the equation.
549 c cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
550 c
551 c *** Calculate certain parcel quantities, including static energy ***
552 c
553 c
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
558 c
559 c
560 c *** Find lifted parcel quantities above cloud base ***
561 c
562 c
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)
569 c
570 c First iteration.
571 c
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))
585 c
586 c Second iteration.
587 c
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))
601 c
602  alv=lv0-clmcpv*(t(i,k)-t0)
603 c print*,'cpd dans convect2 ',cpd
604 c print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
605 c 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
607 c if (.not.cpd.gt.1000.) then
608 c print*,'CPD=',cpd
609 c stop
610 c 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
618 c
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 !=====================================================================
624 c
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
642 c
643 !=====================================================================
644 ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
645 ! --- VIRTUAL TEMPERATURE
646 !=====================================================================
647 c
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))
652 c print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
653 c 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
660 c
661 c=====================================================================
662 c --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
663 c --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
664 c --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
665 c=====================================================================
666 c
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
673 c
674 c Originial Code
675 c
676 c do 530 k=minorig+1,nl-1
677 c do 520 i=1,ncum
678 c if(k.ge.(icb(i)+1))then
679 c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
680 c byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
681 c cape(i)=cape(i)+by
682 c if(by.ge.0.0)inb1(i)=k+1
683 c if(cape(i).gt.0.0)then
684 c inb(i)=k+1
685 c capem(i)=cape(i)
686 c endif
687 c endif
688 c520 continue
689 c530 continue
690 c do 540 i=1,ncum
691 c byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
692 c cape(i)=capem(i)+byp
693 c defrac=capem(i)-cape(i)
694 c defrac=max(defrac,0.001)
695 c frac(i)=-cape(i)/defrac
696 c frac(i)=min(frac(i),1.0)
697 c frac(i)=max(frac(i),0.0)
698 c540 continue
699 c
700 c K Emanuel fix
701 c
702 c call zilch(byp,ncum)
703 c do 530 k=minorig+1,nl-1
704 c do 520 i=1,ncum
705 c if(k.ge.(icb(i)+1))then
706 c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
707 c cape(i)=cape(i)+by
708 c if(by.ge.0.0)inb1(i)=k+1
709 c if(cape(i).gt.0.0)then
710 c inb(i)=k+1
711 c capem(i)=cape(i)
712 c byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
713 c endif
714 c endif
715 c520 continue
716 c530 continue
717 c do 540 i=1,ncum
718 c inb(i)=max(inb(i),inb1(i))
719 c cape(i)=capem(i)+byp(i)
720 c defrac=capem(i)-cape(i)
721 c defrac=max(defrac,0.001)
722 c frac(i)=-cape(i)/defrac
723 c frac(i)=min(frac(i),1.0)
724 c frac(i)=max(frac(i),0.0)
725 c540 continue
726 c
727 c J Teixeira fix
728 c
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
756 c
757 c=====================================================================
758 c --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
759 c=====================================================================
760 c
761 c 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
773 c
774  return
775  end
776 c
777  SUBROUTINE cv_closure(nloc,ncum,nd,nk,icb
778  : ,tv,tvp,p,ph,dph,plcl,cpn
779  : ,iflag,cbmf)
780  implicit none
781 
782 c 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 
789 c outputs:
790  integer iflag(nloc)
791  real cbmf(nloc) ! also an input
792 
793 c 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 
801 c-------------------------------------------------------------------
802 c Compute icbmax.
803 c-------------------------------------------------------------------
804 
805  icbmax=2
806  do 230 i=1,ncum
807  icbmax=max(icbmax,icb(i))
808  230 continue
809 
810 c=====================================================================
811 c --- CALCULATE CLOUD BASE MASS FLUX
812 c=====================================================================
813 c
814 c tvpplcl = parcel temperature lifted adiabatically from level
815 c icb-1 to the LCL.
816 c tvaplcl = virtual temperature at the LCL.
817 c
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 
828 c-------------------------------------------------------------------
829 c --- Interpolate difference between lifted parcel and
830 c --- environmental temperatures to lifted condensation level
831 c-------------------------------------------------------------------
832 c
833 c dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
834 c
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
846 c
847 c-------------------------------------------------------------------
848 c --- Adjust cloud base mass flux
849 c-------------------------------------------------------------------
850 c
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 
871 c 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 
880 c 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 
886 c 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 
895 c=====================================================================
896 c --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
897 c=====================================================================
898 c
899  do 360 i=1,ncum*nlp
900  nent(i,1)=0
901  m(i,1)=0.0
902  360 continue
903 c
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
916 c
917 c-------------------------------------------------------------------
918 c --- Calculate rates of mixing, m(i)
919 c-------------------------------------------------------------------
920 c
921  call zilch(work,ncum)
922 c
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
941 c
942 c
943 c=====================================================================
944 c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
945 c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
946 c --- FRACTION (sij)
947 c=====================================================================
948 c
949 c
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
994 c
995 c *** If no air can entrain at level i assume that updraft detrains ***
996 c *** at that level and calculate detrained air flux and properties ***
997 c
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
1010 c
1011  do 770 i=1,ncum
1012  sij(i,inb(i),inb(i))=1.0
1013  770 continue
1014 c
1015 c=====================================================================
1016 c --- NORMALIZE ENTRAINED AIR MASS FLUXES
1017 c --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
1018 c=====================================================================
1019 c
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
1025 c
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
1031 c
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
1047 c
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
1055 c
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
1118 c
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 
1131 c 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 
1140 c 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 
1145 c 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 
1153 c=====================================================================
1154 c --- PRECIPITATING DOWNDRAFT CALCULATION
1155 c=====================================================================
1156 c
1157 c Initializations:
1158 c
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 
1183 c *** Check whether ep(inb)=0, if so, skip precipitating ***
1184 c *** downdraft calculation ***
1185 c
1186 c
1187 c *** Integrate liquid water equation to find condensed water ***
1188 c *** and condensed water flux ***
1189 c
1190 c
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
1200 c
1201 c *** Begin downdraft loop ***
1202 c
1203 c
1204  call zilch(wdtrain,ncum)
1205  do 899 i=nl+1,1,-1
1206 c
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
1212 c
1213 c
1214 c *** Calculate detrained precipitation ***
1215 c
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
1221 c
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
1233 c
1234 c *** Find rain water and evaporation using provisional ***
1235 c *** estimates of qp(i)and qp(i-1) ***
1236 c
1237 c
1238 c *** Value of terminal velocity and coeffecient of evaporation for snow ***
1239 c
1240  do 894 ij=1,ncum
1241  if((i.le.inb(ij)).and.(lwork(ij)))then
1242  coeff=coeffs
1243  wt(ij,i)=omtsnow
1244 c
1245 c *** Value of terminal velocity and coeffecient of evaporation for rain ***
1246 c
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
1263 c
1264 c *** Calculate precipitating downdraft mass flux under ***
1265 c *** hydrostatic approximation ***
1266 c
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)
1272 c
1273 c *** Add small amount of inertia to downdraft ***
1274 c
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)
1277 c
1278 c *** Force mp to decrease linearly to zero ***
1279 c *** between about 950 mb and the surface ***
1280 c
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
1287 c
1288 c *** Find mixing ratio of precipitating downdraft ***
1289 c
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
1318 c
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 
1335 c 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 
1352 c 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 
1361 c 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 
1370 c -- 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 
1391 c
1392 c *** Calculate surface precipitation in mm/day ***
1393 c
1394  do 1190 i=1,ncum
1395  if(iflag(i).le.1)then
1396 cc precip(i)=precip(i)+wt(i,1)*sigd*water(i,1)*3600.*24000.
1397 cc & /(rowl*g)
1398 cc precip(i)=precip(i)*delt/86400.
1399  precip(i) = wt(i,1)*sigd*water(i,1)*86400/g
1400  endif
1401  1190 continue
1402 c
1403 c
1404 c *** Calculate downdraft velocity scale and surface temperature and ***
1405 c *** water vapor fluctuations ***
1406 c
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
1413 c
1414 c *** Calculate tendencies of lowest level potential temperature ***
1415 c *** and mixing ratio ***
1416 c
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
1455 c
1456 c *** Calculate tendencies of potential temperature and mixing ratio ***
1457 c *** at levels above the lowest level ***
1458 c
1459 c *** First find the net saturated updraft and downdraft mass fluxes ***
1460 c *** through each level ***
1461 c
1462  do 1500 i=2,nl+1
1463 c
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
1469 c
1470  call zilch(amp1,ncum)
1471  call zilch(ad,ncum)
1472 c
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
1481 c
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
1500 c
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)
1505 c
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))
1534 c (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
1564 C (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
1569 C (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
1580 c
1581 c *** Adjust tendencies at top of convection layer to reflect ***
1582 c *** actual position of the level zero cape ***
1583 c
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
1608 c
1609 c *** Very slightly adjust tendencies to force exact ***
1610 c *** enthalpy, momentum and tracer conservation ***
1611 c
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
1635 c
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
1641 c
1642 c
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 
1672 c
1673 c *** diagnose the in-cloud mixing ratio *** ! cld
1674 c *** of condensed water *** ! cld
1675 c ! 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 
1724 c 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 
1733 c 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 
1740 c 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