Context Navigation

cv_routines.F @ 5425

Last change on this file since 5425 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

Rev	Line
[524]	1	!
[1299]	2	! $Id: cv_routines.F 1299 2010-01-20 14:27:21Z jyg $
[524]	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
[1299]	40	CHARACTER (LEN=20) :: modname='cv_routines'
	41	CHARACTER (LEN=80) :: abort_message
[524]	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
[1299]	433	CHARACTER (LEN=20) :: modname='cv_compress'
	434	CHARACTER (LEN=80) :: abort_message
[524]	435
[1299]	436	include 'iniprint.h'
[524]	437
[1299]	438
[524]	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
[1299]	464	write(lunout,*)'strange! nn not equal to ncum: ',nn,ncum
	465	abort_message = ''
	466	CALL abort_gcm (modname,abort_message,1)
[524]	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