Context Navigation

cv_routines.F @ 965

Last change on this file since 965 was 868, checked in by Laurent Fairhead, 17 years ago
Preparation du remplacement de la physique utilisee pour l'exercice IPCC_AR4 par la version de la physique avec thermique. On garde le repertoire phylmd pour un petit moment pour que les utilisateurs ne soient pas trop perdus ... phy_IPCC_AR4 = phylmd LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 55.0 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE cv_param(nd)
5	implicit none
6
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))+cpvq(i,k)
102	cpx(i,k)=cpd(1.0-q(i,k))+clq(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	! cptp+Lqp+phi=cptnk+Lqnk+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))+clqnk(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+alvalvqg/(rrvticb(i)ticb(i))
321	s=1./s
322	ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+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.112exp(17.67tc/denom)
329	else
330	es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
331	endif
332	qg=epses/(p(i,icb(i))-es(1.-eps))
333	c
334	c Second iteration.
335	c
336	s=cpd+alvalvqg/(rrvticb(i)ticb(i))
337	s=1./s
338	ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+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.112exp(17.67tc/denom)
345	else
346	es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
347	end if
348	qg=epses/(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.+rgepsi)
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 cptp+Lqp+phi=cptnk+Lqnk+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))+clqnk(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+alvalvqg/(rrvt(i,k)t(i,k))
566	s=1./s
567	ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+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.112exp(17.67tc/denom)
574	else
575	es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
576	endif
577	qg=epses/(p(i,k)-es(1.-eps))
578	c
579	c Second iteration.
580	c
581	s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
582	s=1./s
583	ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+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.112exp(17.67tc/denom)
590	else
591	es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
592	endif
593	qg=epses/(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.+rgepsi)
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	& -rrdtvp(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.1alpha*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	& +entp0.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	& /(rrvt(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)-cwatbf2)
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)=gep(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)+gawatment(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=coeffph(ij,i)(qs(ij,i)-qsm)
1246	& /(1.0e4+2.0e3ph(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(b6b6+4.*c6))
1254	evap(ij,i)=sigtafacrevap
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.ginvlv(ij,i)sigdevap(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	real delt
1333	real t(nloc,nd), q(nloc,nd), u(nloc,nd), v(nloc,nd)
1334	real gz(nloc,nd)
1335	real p(nloc,nd), ph(nloc,nd+1), h(nloc,nd)
1336	real hp(nloc,nd), lv(nloc,nd)
1337	real cpn(nloc,nd), ep(nloc,nd), clw(nloc,nd), frac(nloc)
1338	real m(nloc,nd), mp(nloc,nd), qp(nloc,nd)
1339	real up(nloc,nd), vp(nloc,nd)
1340	real wt(nloc,nd), water(nloc,nd), evap(nloc,nd)
1341	real ment(nloc,nd,nd), qent(nloc,nd,nd), elij(nloc,nd,nd)
1342	real uent(nloc,nd,nd), vent(nloc,nd,nd)
1343	real tv(nloc,nd), tvp(nloc,nd)
1344
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	cc precip(i)=precip(i)+wt(i,1)sigdwater(i,1)3600.24000.
1390	cc & /(rowl*g)
1391	cc precip(i)=precip(i)*delt/86400.
1392	precip(i) = wt(i,1)sigdwater(i,1)*86400/g
1393	endif
1394	1190 continue
1395	c
1396	c
1397	c * Calculate downdraft velocity scale and surface temperature and *
1398	c * water vapor fluctuations *
1399	c
1400	do i=1,ncum
1401	wd(i)=betadabs(mp(i,icb(i)))0.01rrdt(i,icb(i))
1402	: /(sigd*p(i,icb(i)))
1403	qprime(i)=0.5*(qp(i,1)-q(i,1))
1404	tprime(i)=lv0*qprime(i)/cpd
1405	enddo
1406	c
1407	c * Calculate tendencies of lowest level potential temperature *
1408	c * and mixing ratio *
1409	c
1410	do 1200 i=1,ncum
1411	work(i)=0.01/(ph(i,1)-ph(i,2))
1412	am(i)=0.0
1413	1200 continue
1414	do 1220 k=2,nl
1415	do 1210 i=1,ncum
1416	if((nk(i).eq.1).and.(k.le.inb(i)).and.(nk(i).eq.1))then
1417	am(i)=am(i)+m(i,k)
1418	endif
1419	1210 continue
1420	1220 continue
1421	do 1240 i=1,ncum
1422	if((gwork(i)am(i)).ge.delti)iflag(i)=1
1423	ft(i,1)=ft(i,1)+gwork(i)am(i)*(t(i,2)-t(i,1)
1424	& +(gz(i,2)-gz(i,1))/cpn(i,1))
1425	ft(i,1)=ft(i,1)-lvcp(i,1)sigdevap(i,1)
1426	ft(i,1)=ft(i,1)+sigdwt(i,2)(cl-cpd)water(i,2)(t(i,2)
1427	& -t(i,1))*work(i)/cpn(i,1)
1428	fq(i,1)=fq(i,1)+gmp(i,2)(qp(i,2)-q(i,1))*
1429	& work(i)+sigd*evap(i,1)
1430	fq(i,1)=fq(i,1)+gam(i)(q(i,2)-q(i,1))*work(i)
1431	fu(i,1)=fu(i,1)+gwork(i)(mp(i,2)*(up(i,2)-u(i,1))
1432	& +am(i)*(u(i,2)-u(i,1)))
1433	fv(i,1)=fv(i,1)+gwork(i)(mp(i,2)*(vp(i,2)-v(i,1))
1434	& +am(i)*(v(i,2)-v(i,1)))
1435	1240 continue
1436	do 1260 j=2,nl
1437	do 1250 i=1,ncum
1438	if(j.le.inb(i))then
1439	fq(i,1)=fq(i,1)
1440	& +gwork(i)ment(i,j,1)*(qent(i,j,1)-q(i,1))
1441	fu(i,1)=fu(i,1)
1442	& +gwork(i)ment(i,j,1)*(uent(i,j,1)-u(i,1))
1443	fv(i,1)=fv(i,1)
1444	& +gwork(i)ment(i,j,1)*(vent(i,j,1)-v(i,1))
1445	endif
1446	1250 continue
1447	1260 continue
1448	c
1449	c * Calculate tendencies of potential temperature and mixing ratio *
1450	c * at levels above the lowest level *
1451	c
1452	c * First find the net saturated updraft and downdraft mass fluxes *
1453	c * through each level *
1454	c
1455	do 1500 i=2,nl+1
1456	c
1457	num1=0
1458	do 1265 ij=1,ncum
1459	if(i.le.inb(ij))num1=num1+1
1460	1265 continue
1461	if(num1.le.0)go to 1500
1462	c
1463	call zilch(amp1,ncum)
1464	call zilch(ad,ncum)
1465	c
1466	do 1280 k=i+1,nl+1
1467	do 1270 ij=1,ncum
1468	if((i.ge.nk(ij)).and.(i.le.inb(ij))
1469	& .and.(k.le.(inb(ij)+1)))then
1470	amp1(ij)=amp1(ij)+m(ij,k)
1471	endif
1472	1270 continue
1473	1280 continue
1474	c
1475	do 1310 k=1,i
1476	do 1300 j=i+1,nl+1
1477	do 1290 ij=1,ncum
1478	if((j.le.(inb(ij)+1)).and.(i.le.inb(ij)))then
1479	amp1(ij)=amp1(ij)+ment(ij,k,j)
1480	endif
1481	1290 continue
1482	1300 continue
1483	1310 continue
1484	do 1340 k=1,i-1
1485	do 1330 j=i,nl+1
1486	do 1320 ij=1,ncum
1487	if((i.le.inb(ij)).and.(j.le.inb(ij)))then
1488	ad(ij)=ad(ij)+ment(ij,j,k)
1489	endif
1490	1320 continue
1491	1330 continue
1492	1340 continue
1493	c
1494	do 1350 ij=1,ncum
1495	if(i.le.inb(ij))then
1496	dpinv=0.01/(ph(ij,i)-ph(ij,i+1))
1497	cpinv=1.0/cpn(ij,i)
1498	c
1499	ft(ij,i)=ft(ij,i)
1500	& +gdpinv(amp1(ij)*(t(ij,i+1)-t(ij,i)
1501	& +(gz(ij,i+1)-gz(ij,i))*cpinv)
1502	& -ad(ij)(t(ij,i)-t(ij,i-1)+(gz(ij,i)-gz(ij,i-1))cpinv))
1503	& -sigdlvcp(ij,i)evap(ij,i)
1504	ft(ij,i)=ft(ij,i)+gdpinvment(ij,i,i)*(hp(ij,i)-h(ij,i)+
1505	& t(ij,i)(cpv-cpd)(q(ij,i)-qent(ij,i,i)))*cpinv
1506	ft(ij,i)=ft(ij,i)+sigdwt(ij,i+1)(cl-cpd)water(ij,i+1)
1507	& (t(ij,i+1)-t(ij,i))dpinvcpinv
1508	fq(ij,i)=fq(ij,i)+gdpinv(amp1(ij)*(q(ij,i+1)-q(ij,i))-
1509	& ad(ij)*(q(ij,i)-q(ij,i-1)))
1510	fu(ij,i)=fu(ij,i)+gdpinv(amp1(ij)*(u(ij,i+1)-u(ij,i))-
1511	& ad(ij)*(u(ij,i)-u(ij,i-1)))
1512	fv(ij,i)=fv(ij,i)+gdpinv(amp1(ij)*(v(ij,i+1)-v(ij,i))-
1513	& ad(ij)*(v(ij,i)-v(ij,i-1)))
1514	endif
1515	1350 continue
1516	do 1370 k=1,i-1
1517	do 1360 ij=1,ncum
1518	if(i.le.inb(ij))then
1519	awat=elij(ij,k,i)-(1.-ep(ij,i))*clw(ij,i)
1520	awat=max(awat,0.0)
1521	fq(ij,i)=fq(ij,i)
1522	& +gdpinvment(ij,k,i)*(qent(ij,k,i)-awat-q(ij,i))
1523	fu(ij,i)=fu(ij,i)
1524	& +gdpinvment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
1525	fv(ij,i)=fv(ij,i)
1526	& +gdpinvment(ij,k,i)*(vent(ij,k,i)-v(ij,i))
1527	c (saturated updrafts resulting from mixing) ! cld
1528	qcond(ij,i)=qcond(ij,i)+(elij(ij,k,i)-awat) ! cld
1529	nqcond(ij,i)=nqcond(ij,i)+1. ! cld
1530	endif
1531	1360 continue
1532	1370 continue
1533	do 1390 k=i,nl+1
1534	do 1380 ij=1,ncum
1535	if((i.le.inb(ij)).and.(k.le.inb(ij)))then
1536	fq(ij,i)=fq(ij,i)
1537	& +gdpinvment(ij,k,i)*(qent(ij,k,i)-q(ij,i))
1538	fu(ij,i)=fu(ij,i)
1539	& +gdpinvment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
1540	fv(ij,i)=fv(ij,i)
1541	& +gdpinvment(ij,k,i)*(vent(ij,k,i)-v(ij,i))
1542	endif
1543	1380 continue
1544	1390 continue
1545	do 1400 ij=1,ncum
1546	if(i.le.inb(ij))then
1547	fq(ij,i)=fq(ij,i)
1548	& +sigdevap(ij,i)+g(mp(ij,i+1)*
1549	& (qp(ij,i+1)-q(ij,i))
1550	& -mp(ij,i)(qp(ij,i)-q(ij,i-1)))dpinv
1551	fu(ij,i)=fu(ij,i)
1552	& +g(mp(ij,i+1)(up(ij,i+1)-u(ij,i))-mp(ij,i)*
1553	& (up(ij,i)-u(ij,i-1)))*dpinv
1554	fv(ij,i)=fv(ij,i)
1555	& +g(mp(ij,i+1)(vp(ij,i+1)-v(ij,i))-mp(ij,i)*
1556	& (vp(ij,i)-v(ij,i-1)))*dpinv
1557	C (saturated downdrafts resulting from mixing) ! cld
1558	do k=i+1,inb(ij) ! cld
1559	qcond(ij,i)=qcond(ij,i)+elij(ij,k,i) ! cld
1560	nqcond(ij,i)=nqcond(ij,i)+1. ! cld
1561	enddo ! cld
1562	C (particular case: no detraining level is found) ! cld
1563	if (nent(ij,i).eq.0) then ! cld
1564	qcond(ij,i)=qcond(ij,i)+(1.-ep(ij,i))*clw(ij,i) ! cld
1565	nqcond(ij,i)=nqcond(ij,i)+1. ! cld
1566	endif ! cld
1567	if (nqcond(ij,i).ne.0.) then ! cld
1568	qcond(ij,i)=qcond(ij,i)/nqcond(ij,i) ! cld
1569	endif ! cld
1570	endif
1571	1400 continue
1572	1500 continue
1573	c
1574	c * Adjust tendencies at top of convection layer to reflect *
1575	c * actual position of the level zero cape *
1576	c
1577	do 503 ij=1,ncum
1578	fqold=fq(ij,inb(ij))
1579	fq(ij,inb(ij))=fq(ij,inb(ij))*(1.-frac(ij))
1580	fq(ij,inb(ij)-1)=fq(ij,inb(ij)-1)
1581	& +frac(ij)fqold((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
1582	1 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*lv(ij,inb(ij))
1583	& /lv(ij,inb(ij)-1)
1584	ftold=ft(ij,inb(ij))
1585	ft(ij,inb(ij))=ft(ij,inb(ij))*(1.-frac(ij))
1586	ft(ij,inb(ij)-1)=ft(ij,inb(ij)-1)
1587	& +frac(ij)ftold((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
1588	1 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*cpn(ij,inb(ij))
1589	& /cpn(ij,inb(ij)-1)
1590	fuold=fu(ij,inb(ij))
1591	fu(ij,inb(ij))=fu(ij,inb(ij))*(1.-frac(ij))
1592	fu(ij,inb(ij)-1)=fu(ij,inb(ij)-1)
1593	& +frac(ij)fuold((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
1594	1 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
1595	fvold=fv(ij,inb(ij))
1596	fv(ij,inb(ij))=fv(ij,inb(ij))*(1.-frac(ij))
1597	fv(ij,inb(ij)-1)=fv(ij,inb(ij)-1)
1598	& +frac(ij)fvold((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
1599	1 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
1600	503 continue
1601	c
1602	c * Very slightly adjust tendencies to force exact *
1603	c * enthalpy, momentum and tracer conservation *
1604	c
1605	do 682 ij=1,ncum
1606	ents(ij)=0.0
1607	uav(ij)=0.0
1608	vav(ij)=0.0
1609	do 681 i=1,inb(ij)
1610	ents(ij)=ents(ij)
1611	& +(cpn(ij,i)ft(ij,i)+lv(ij,i)fq(ij,i))*(ph(ij,i)-ph(ij,i+1))
1612	uav(ij)=uav(ij)+fu(ij,i)*(ph(ij,i)-ph(ij,i+1))
1613	vav(ij)=vav(ij)+fv(ij,i)*(ph(ij,i)-ph(ij,i+1))
1614	681 continue
1615	682 continue
1616	do 683 ij=1,ncum
1617	ents(ij)=ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1618	uav(ij)=uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1619	vav(ij)=vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1620	683 continue
1621	do 642 ij=1,ncum
1622	do 641 i=1,inb(ij)
1623	ft(ij,i)=ft(ij,i)-ents(ij)/cpn(ij,i)
1624	fu(ij,i)=(1.-cu)*(fu(ij,i)-uav(ij))
1625	fv(ij,i)=(1.-cu)*(fv(ij,i)-vav(ij))
1626	641 continue
1627	642 continue
1628	c
1629	do 1810 k=1,nl+1
1630	do 1800 i=1,ncum
1631	if((q(i,k)+delt*fq(i,k)).lt.0.0)iflag(i)=10
1632	1800 continue
1633	1810 continue
1634	c
1635	c
1636	do 1900 i=1,ncum
1637	if(iflag(i).gt.2)then
1638	precip(i)=0.0
1639	cbmf(i)=0.0
1640	endif
1641	1900 continue
1642	do 1920 k=1,nl
1643	do 1910 i=1,ncum
1644	if(iflag(i).gt.2)then
1645	ft(i,k)=0.0
1646	fq(i,k)=0.0
1647	fu(i,k)=0.0
1648	fv(i,k)=0.0
1649	qcondc(i,k)=0.0 ! cld
1650	endif
1651	1910 continue
1652	1920 continue
1653
1654	do k=1,nl+1
1655	do i=1,ncum
1656	Ma(i,k) = 0.
1657	enddo
1658	enddo
1659	do k=nl,1,-1
1660	do i=1,ncum
1661	Ma(i,k) = Ma(i,k+1)+m(i,k)
1662	enddo
1663	enddo
1664
1665	c
1666	c * diagnose the in-cloud mixing ratio * ! cld
1667	c * of condensed water * ! cld
1668	c ! cld
1669	DO ij=1,ncum ! cld
1670	do i=1,nd ! cld
1671	mac(ij,i)=0.0 ! cld
1672	wa(ij,i)=0.0 ! cld
1673	siga(ij,i)=0.0 ! cld
1674	enddo ! cld
1675	do i=nk(ij),inb(ij) ! cld
1676	do k=i+1,inb(ij)+1 ! cld
1677	mac(ij,i)=mac(ij,i)+m(ij,k) ! cld
1678	enddo ! cld
1679	enddo ! cld
1680	do i=icb(ij),inb(ij)-1 ! cld
1681	ax(ij,i)=0. ! cld
1682	do j=icb(ij),i ! cld
1683	ax(ij,i)=ax(ij,i)+rrd*(tvp(ij,j)-tv(ij,j)) ! cld
1684	: *(ph(ij,j)-ph(ij,j+1))/p(ij,j) ! cld
1685	enddo ! cld
1686	if (ax(ij,i).gt.0.0) then ! cld
1687	wa(ij,i)=sqrt(2.*ax(ij,i)) ! cld
1688	endif ! cld
1689	enddo ! cld
1690	do i=1,nl ! cld
1691	if (wa(ij,i).gt.0.0) ! cld
1692	: siga(ij,i)=mac(ij,i)/wa(ij,i) ! cld
1693	: rrdtvp(ij,i)/p(ij,i)/100./delta ! cld
1694	siga(ij,i) = min(siga(ij,i),1.0) ! cld
1695	qcondc(ij,i)=siga(ij,i)clw(ij,i)(1.-ep(ij,i)) ! cld
1696	: + (1.-siga(ij,i))*qcond(ij,i) ! cld
1697	enddo ! cld
1698	ENDDO ! cld
1699
1700	return
1701	end
1702
1703	SUBROUTINE cv_uncompress(nloc,len,ncum,nd,idcum
1704	: ,iflag
1705	: ,precip,cbmf
1706	: ,ft,fq,fu,fv
1707	: ,Ma,qcondc
1708	: ,iflag1
1709	: ,precip1,cbmf1
1710	: ,ft1,fq1,fu1,fv1
1711	: ,Ma1,qcondc1
1712	: )
1713	implicit none
1714
1715	#include "cvparam.h"
1716
1717	c inputs:
1718	integer len, ncum, nd, nloc
1719	integer idcum(nloc)
1720	integer iflag(nloc)
1721	real precip(nloc), cbmf(nloc)
1722	real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
1723	real Ma(nloc,nd)
1724	real qcondc(nloc,nd) !cld
1725
1726	c outputs:
1727	integer iflag1(len)
1728	real precip1(len), cbmf1(len)
1729	real ft1(len,nd), fq1(len,nd), fu1(len,nd), fv1(len,nd)
1730	real Ma1(len,nd)
1731	real qcondc1(len,nd) !cld
1732
1733	c local variables:
1734	integer i,k
1735
1736	do 2000 i=1,ncum
1737	precip1(idcum(i))=precip(i)
1738	cbmf1(idcum(i))=cbmf(i)
1739	iflag1(idcum(i))=iflag(i)
1740	2000 continue
1741
1742	do 2020 k=1,nl
1743	do 2010 i=1,ncum
1744	ft1(idcum(i),k)=ft(i,k)
1745	fq1(idcum(i),k)=fq(i,k)
1746	fu1(idcum(i),k)=fu(i,k)
1747	fv1(idcum(i),k)=fv(i,k)
1748	Ma1(idcum(i),k)=Ma(i,k)
1749	qcondc1(idcum(i),k)=qcondc(i,k)
1750	2010 continue
1751	2020 continue
1752
1753	return
1754	end
1755

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format