Context Navigation

cv_routines.F @ 1373

Last change on this file since 1373 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

tous les FLOAT() sont remplacés par des REAL()
tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

all FLOAT() instructions have been replaced by REAL() instructions
all STOP instructions in phylmd have been replaced by calls to abort_gcm
the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
$Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 55.3 KB

Line
1	!
2	! $Id: cv_routines.F 1299 2010-01-20 14:27:21Z musat $
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))+cpvq(i,k)
104	cpx(i,k)=cpd(1.0-q(i,k))+clq(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	! cptp+Lqp+phi=cptnk+Lqnk+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))+clqnk(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+alvalvqg/(rrvticb(i)ticb(i))
323	s=1./s
324	ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+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.112exp(17.67tc/denom)
331	else
332	es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
333	endif
334	qg=epses/(p(i,icb(i))-es(1.-eps))
335	c
336	c Second iteration.
337	c
338	s=cpd+alvalvqg/(rrvticb(i)ticb(i))
339	s=1./s
340	ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+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.112exp(17.67tc/denom)
347	else
348	es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
349	end if
350	qg=epses/(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.+rgepsi)
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 cptp+Lqp+phi=cptnk+Lqnk+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))+clqnk(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+alvalvqg/(rrvt(i,k)t(i,k))
573	s=1./s
574	ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+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.112exp(17.67tc/denom)
581	else
582	es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
583	endif
584	qg=epses/(p(i,k)-es(1.-eps))
585	c
586	c Second iteration.
587	c
588	s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
589	s=1./s
590	ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+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.112exp(17.67tc/denom)
597	else
598	es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
599	endif
600	qg=epses/(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.+rgepsi)
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	& -rrdtvp(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.1alpha*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	& +entp0.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	& /(rrvt(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)-cwatbf2)
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)=gep(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)+gawatment(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=coeffph(ij,i)(qs(ij,i)-qsm)
1253	& /(1.0e4+2.0e3ph(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(b6b6+4.*c6))
1261	evap(ij,i)=sigtafacrevap
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.ginvlv(ij,i)sigdevap(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)sigdwater(i,1)3600.24000.
1397	cc & /(rowl*g)
1398	cc precip(i)=precip(i)*delt/86400.
1399	precip(i) = wt(i,1)sigdwater(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)=betadabs(mp(i,icb(i)))0.01rrdt(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((gwork(i)am(i)).ge.delti)iflag(i)=1
1430	ft(i,1)=ft(i,1)+gwork(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)sigdevap(i,1)
1433	ft(i,1)=ft(i,1)+sigdwt(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)+gmp(i,2)(qp(i,2)-q(i,1))*
1436	& work(i)+sigd*evap(i,1)
1437	fq(i,1)=fq(i,1)+gam(i)(q(i,2)-q(i,1))*work(i)
1438	fu(i,1)=fu(i,1)+gwork(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)+gwork(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	& +gwork(i)ment(i,j,1)*(qent(i,j,1)-q(i,1))
1448	fu(i,1)=fu(i,1)
1449	& +gwork(i)ment(i,j,1)*(uent(i,j,1)-u(i,1))
1450	fv(i,1)=fv(i,1)
1451	& +gwork(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	& +gdpinv(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	& -sigdlvcp(ij,i)evap(ij,i)
1511	ft(ij,i)=ft(ij,i)+gdpinvment(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)+sigdwt(ij,i+1)(cl-cpd)water(ij,i+1)
1514	& (t(ij,i+1)-t(ij,i))dpinvcpinv
1515	fq(ij,i)=fq(ij,i)+gdpinv(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)+gdpinv(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)+gdpinv(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	& +gdpinvment(ij,k,i)*(qent(ij,k,i)-awat-q(ij,i))
1530	fu(ij,i)=fu(ij,i)
1531	& +gdpinvment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
1532	fv(ij,i)=fv(ij,i)
1533	& +gdpinvment(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	& +gdpinvment(ij,k,i)*(qent(ij,k,i)-q(ij,i))
1545	fu(ij,i)=fu(ij,i)
1546	& +gdpinvment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
1547	fv(ij,i)=fv(ij,i)
1548	& +gdpinvment(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	& +sigdevap(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	: rrdtvp(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

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

Download in other formats:

Original Format