Context Navigation

cv30_routines.F @ 3969

Last change on this file since 3969 was 1403, checked in by Laurent Fairhead, 14 years ago

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

IPSLCM5A configuration (comparison between trunk and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

IPSLCM5A configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

MH07 configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files

Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

pour une config IPSLCM5A (comparaison tronc/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux (à part sur RN et Pb)
- fichiers histoire égaux

pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux
- fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

pour une config IPSLCM5A:
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

pour une config MH07
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

Changement sur makegcm et create_make-gcm pour pouvoir prendre en compte des
programmes principaux en *F90

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

File size: 93.8 KB

Rev	Line
[879]	1	!
[1403]	2	! $Id: cv30_routines.F 1403 2010-07-01 09:02:53Z fairhead $
[879]	3	!
	4	c
	5	c
	6	SUBROUTINE cv30_param(nd,delt)
	7	implicit none
	8
	9	c------------------------------------------------------------
	10	c Set parameters for convectL for iflag_con = 3
	11	c------------------------------------------------------------
	12
	13	C
	14	C * PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE *
	15	C * PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO *
	16	C * PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. *
	17	C * EFFICIENCY IS ASSUMED TO BE UNITY *
	18	C * SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *
	19	C * SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *
	20	C * OF CLOUD *
	21	C
	22	C [TAU: CHARACTERISTIC TIMESCALE USED TO COMPUTE ALPHA & BETA]
	23	C * ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF *
	24	C * APPROACH TO QUASI-EQUILIBRIUM *
	25	C * (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) *
	26	C * (BETA MUST BE LESS THAN OR EQUAL TO 1) *
	27	C
	28	C * DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE *
	29	C * APPROACH TO QUASI-EQUILIBRIUM *
	30	C * IT MUST BE LESS THAN 0 *
	31
	32	#include "cv30param.h"
	33	#include "conema3.h"
	34
	35	integer nd
	36	real delt ! timestep (seconds)
	37
	38	c noff: integer limit for convection (nd-noff)
	39	c minorig: First level of convection
	40
	41	c -- limit levels for convection:
	42
	43	noff = 1
	44	minorig = 1
	45	nl=nd-noff
	46	nlp=nl+1
	47	nlm=nl-1
	48
	49	c -- "microphysical" parameters:
	50
	51	sigd = 0.01
	52	spfac = 0.15
	53	pbcrit = 150.0
	54	ptcrit = 500.0
	55	cIM cf. FH epmax = 0.993
	56
	57	omtrain = 45.0 ! used also for snow (no disctinction rain/snow)
	58
	59	c -- misc:
	60
	61	dtovsh = -0.2 ! dT for overshoot
	62	dpbase = -40. ! definition cloud base (400m above LCL)
	63	dttrig = 5. ! (loose) condition for triggering
	64
	65	c -- rate of approach to quasi-equilibrium:
	66
	67	dtcrit = -2.0
	68	tau = 8000.
	69	beta = 1.0 - delt/tau
	70	alpha = 1.5E-3 * delt/tau
	71	c increase alpha to compensate W decrease:
	72	alpha = alpha*1.5
	73
	74	c -- interface cloud parameterization:
	75
	76	delta=0.01 ! cld
	77
	78	c -- interface with boundary-layer (gust factor): (sb)
	79
	80	betad=10.0 ! original value (from convect 4.3)
	81
	82	return
	83	end
	84
	85	SUBROUTINE cv30_prelim(len,nd,ndp1,t,q,p,ph
	86	: ,lv,cpn,tv,gz,h,hm,th)
	87	implicit none
	88
	89	!=====================================================================
	90	! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
	91	! "ori": from convect4.3 (vectorized)
	92	! "convect3": to be exactly consistent with convect3
	93	!=====================================================================
	94
	95	c inputs:
	96	integer len, nd, ndp1
	97	real t(len,nd), q(len,nd), p(len,nd), ph(len,ndp1)
	98
	99	c outputs:
	100	real lv(len,nd), cpn(len,nd), tv(len,nd)
	101	real gz(len,nd), h(len,nd), hm(len,nd)
	102	real th(len,nd)
	103
	104	c local variables:
	105	integer k, i
	106	real rdcp
	107	real tvx,tvy ! convect3
	108	real cpx(len,nd)
	109
	110	#include "cvthermo.h"
	111	#include "cv30param.h"
	112
	113
	114	c ori do 110 k=1,nlp
	115	do 110 k=1,nl ! convect3
	116	do 100 i=1,len
	117	cdebug lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
	118	lv(i,k)= lv0-clmcpv*(t(i,k)-273.15)
	119	cpn(i,k)=cpd(1.0-q(i,k))+cpvq(i,k)
	120	cpx(i,k)=cpd(1.0-q(i,k))+clq(i,k)
	121	c ori tv(i,k)=t(i,k)(1.0+q(i,k)epsim1)
	122	tv(i,k)=t(i,k)*(1.0+q(i,k)/eps-q(i,k))
	123	rdcp=(rrd(1.-q(i,k))+q(i,k)rrv)/cpn(i,k)
	124	th(i,k)=t(i,k)(1000.0/p(i,k))*rdcp
	125	100 continue
	126	110 continue
	127	c
	128	c gz = phi at the full levels (same as p).
	129	c
	130	do 120 i=1,len
	131	gz(i,1)=0.0
	132	120 continue
	133	c ori do 140 k=2,nlp
	134	do 140 k=2,nl ! convect3
	135	do 130 i=1,len
	136	tvx=t(i,k)*(1.+q(i,k)/eps-q(i,k)) !convect3
	137	tvy=t(i,k-1)*(1.+q(i,k-1)/eps-q(i,k-1)) !convect3
	138	gz(i,k)=gz(i,k-1)+0.5rrd(tvx+tvy) !convect3
	139	& *(p(i,k-1)-p(i,k))/ph(i,k) !convect3
	140
	141	c ori gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
	142	c ori & *(p(i,k-1)-p(i,k))/ph(i,k)
	143	130 continue
	144	140 continue
	145	c
	146	c h = phi + cpT (dry static energy).
	147	c hm = phi + cp(T-Tbase)+Lq
	148	c
	149	c ori do 170 k=1,nlp
	150	do 170 k=1,nl ! convect3
	151	do 160 i=1,len
	152	h(i,k)=gz(i,k)+cpn(i,k)*t(i,k)
	153	hm(i,k)=gz(i,k)+cpx(i,k)(t(i,k)-t(i,1))+lv(i,k)q(i,k)
	154	160 continue
	155	170 continue
	156
	157	return
	158	end
	159
	160	SUBROUTINE cv30_feed(len,nd,t,q,qs,p,ph,hm,gz
	161	: ,nk,icb,icbmax,iflag,tnk,qnk,gznk,plcl)
	162	implicit none
	163
	164	C================================================================
	165	C Purpose: CONVECTIVE FEED
	166	C
	167	C Main differences with cv_feed:
	168	C - ph added in input
	169	C - here, nk(i)=minorig
	170	C - icb defined differently (plcl compared with ph instead of p)
	171	C
	172	C Main differences with convect3:
	173	C - we do not compute dplcldt and dplcldr of CLIFT anymore
	174	C - values iflag different (but tests identical)
	175	C - A,B explicitely defined (!...)
	176	C================================================================
	177
	178	#include "cv30param.h"
	179
	180	c inputs:
	181	integer len, nd
	182	real t(len,nd), q(len,nd), qs(len,nd), p(len,nd)
	183	real hm(len,nd), gz(len,nd)
	184	real ph(len,nd+1)
	185
	186	c outputs:
	187	integer iflag(len), nk(len), icb(len), icbmax
	188	real tnk(len), qnk(len), gznk(len), plcl(len)
	189
	190	c local variables:
	191	integer i, k
	192	integer ihmin(len)
	193	real work(len)
	194	real pnk(len), qsnk(len), rh(len), chi(len)
	195	real A, B ! convect3
	196	cym
	197	plcl=0.0
	198	c@ !-------------------------------------------------------------------
	199	c@ ! --- Find level of minimum moist static energy
	200	c@ ! --- If level of minimum moist static energy coincides with
	201	c@ ! --- or is lower than minimum allowable parcel origin level,
	202	c@ ! --- set iflag to 6.
	203	c@ !-------------------------------------------------------------------
	204	c@
	205	c@ do 180 i=1,len
	206	c@ work(i)=1.0e12
	207	c@ ihmin(i)=nl
	208	c@ 180 continue
	209	c@ do 200 k=2,nlp
	210	c@ do 190 i=1,len
	211	c@ if((hm(i,k).lt.work(i)).and.
	212	c@ & (hm(i,k).lt.hm(i,k-1)))then
	213	c@ work(i)=hm(i,k)
	214	c@ ihmin(i)=k
	215	c@ endif
	216	c@ 190 continue
	217	c@ 200 continue
	218	c@ do 210 i=1,len
	219	c@ ihmin(i)=min(ihmin(i),nlm)
	220	c@ if(ihmin(i).le.minorig)then
	221	c@ iflag(i)=6
	222	c@ endif
	223	c@ 210 continue
	224	c@ c
	225	c@ !-------------------------------------------------------------------
	226	c@ ! --- Find that model level below the level of minimum moist static
	227	c@ ! --- energy that has the maximum value of moist static energy
	228	c@ !-------------------------------------------------------------------
	229	c@
	230	c@ do 220 i=1,len
	231	c@ work(i)=hm(i,minorig)
	232	c@ nk(i)=minorig
	233	c@ 220 continue
	234	c@ do 240 k=minorig+1,nl
	235	c@ do 230 i=1,len
	236	c@ if((hm(i,k).gt.work(i)).and.(k.le.ihmin(i)))then
	237	c@ work(i)=hm(i,k)
	238	c@ nk(i)=k
	239	c@ endif
	240	c@ 230 continue
	241	c@ 240 continue
	242
	243	!-------------------------------------------------------------------
	244	! --- Origin level of ascending parcels for convect3:
	245	!-------------------------------------------------------------------
	246
	247	do 220 i=1,len
	248	nk(i)=minorig
	249	220 continue
	250
	251	!-------------------------------------------------------------------
	252	! --- Check whether parcel level temperature and specific humidity
	253	! --- are reasonable
	254	!-------------------------------------------------------------------
	255	do 250 i=1,len
	256	if( ( ( t(i,nk(i)).lt.250.0 )
	257	& .or.( q(i,nk(i)).le.0.0 ) )
	258	c@ & .or.( p(i,ihmin(i)).lt.400.0 ) )
	259	& .and.
	260	& ( iflag(i).eq.0) ) iflag(i)=7
	261	250 continue
	262	!-------------------------------------------------------------------
	263	! --- Calculate lifted condensation level of air at parcel origin level
	264	! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
	265	!-------------------------------------------------------------------
	266
	267	A = 1669.0 ! convect3
	268	B = 122.0 ! convect3
	269
	270	do 260 i=1,len
	271
	272	if (iflag(i).ne.7) then ! modif sb Jun7th 2002
	273
	274	tnk(i)=t(i,nk(i))
	275	qnk(i)=q(i,nk(i))
	276	gznk(i)=gz(i,nk(i))
	277	pnk(i)=p(i,nk(i))
	278	qsnk(i)=qs(i,nk(i))
	279	c
	280	rh(i)=qnk(i)/qsnk(i)
	281	c ori rh(i)=min(1.0,rh(i)) ! removed for convect3
	282	c ori chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
	283	chi(i)=tnk(i)/(A-B*rh(i)-tnk(i)) ! convect3
	284	plcl(i)=pnk(i)(rh(i)*chi(i))
	285	if(((plcl(i).lt.200.0).or.(plcl(i).ge.2000.0))
	286	& .and.(iflag(i).eq.0))iflag(i)=8
	287
	288	endif ! iflag=7
	289
	290	260 continue
	291
	292	!-------------------------------------------------------------------
	293	! --- Calculate first level above lcl (=icb)
	294	!-------------------------------------------------------------------
	295
	296	c@ do 270 i=1,len
	297	c@ icb(i)=nlm
	298	c@ 270 continue
	299	c@c
	300	c@ do 290 k=minorig,nl
	301	c@ do 280 i=1,len
	302	c@ if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i)))
	303	c@ & icb(i)=min(icb(i),k)
	304	c@ 280 continue
	305	c@ 290 continue
	306	c@c
	307	c@ do 300 i=1,len
	308	c@ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
	309	c@ 300 continue
	310
	311	do 270 i=1,len
	312	icb(i)=nlm
	313	270 continue
	314	c
	315	c la modification consiste a comparer plcl a ph et non a p:
	316	c icb est defini par : ph(icb)<plcl<ph(icb-1)
	317	c@ do 290 k=minorig,nl
	318	do 290 k=3,nl-1 ! modif pour que icb soit sup/egal a 2
	319	do 280 i=1,len
	320	if( ph(i,k).lt.plcl(i) ) icb(i)=min(icb(i),k)
	321	280 continue
	322	290 continue
	323	c
	324	do 300 i=1,len
	325	c@ if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
	326	if((icb(i).eq.nlm).and.(iflag(i).eq.0))iflag(i)=9
	327	300 continue
	328
	329	do 400 i=1,len
	330	icb(i) = icb(i)-1 ! icb sup ou egal a 2
	331	400 continue
	332	c
	333	c Compute icbmax.
	334	c
	335	icbmax=2
	336	do 310 i=1,len
	337	c! icbmax=max(icbmax,icb(i))
	338	if (iflag(i).lt.7) icbmax=max(icbmax,icb(i)) ! sb Jun7th02
	339	310 continue
	340
	341	return
	342	end
	343
	344	SUBROUTINE cv30_undilute1(len,nd,t,q,qs,gz,plcl,p,nk,icb
	345	: ,tp,tvp,clw,icbs)
	346	implicit none
	347
	348	!----------------------------------------------------------------
	349	! Equivalent de TLIFT entre NK et ICB+1 inclus
	350	!
	351	! Differences with convect4:
	352	! - specify plcl in input
	353	! - icbs is the first level above LCL (may differ from icb)
	354	! - in the iterations, used x(icbs) instead x(icb)
	355	! - many minor differences in the iterations
	356	! - tvp is computed in only one time
	357	! - icbs: first level above Plcl (IMIN de TLIFT) in output
	358	! - if icbs=icb, compute also tp(icb+1),tvp(icb+1) & clw(icb+1)
	359	!----------------------------------------------------------------
	360
	361	#include "cvthermo.h"
	362	#include "cv30param.h"
	363
	364	c inputs:
	365	integer len, nd
	366	integer nk(len), icb(len)
	367	real t(len,nd), q(len,nd), qs(len,nd), gz(len,nd)
	368	real p(len,nd)
	369	real plcl(len) ! convect3
	370
	371	c outputs:
	372	real tp(len,nd), tvp(len,nd), clw(len,nd)
	373
	374	c local variables:
	375	integer i, k
	376	integer icb1(len), icbs(len), icbsmax2 ! convect3
	377	real tg, qg, alv, s, ahg, tc, denom, es, rg
	378	real ah0(len), cpp(len)
	379	real tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)
	380	real qsicb(len) ! convect3
	381	real cpinv(len) ! convect3
	382
	383	!-------------------------------------------------------------------
	384	! --- Calculates the lifted parcel virtual temperature at nk,
	385	! --- the actual temperature, and the adiabatic
	386	! --- liquid water content. The procedure is to solve the equation.
	387	! cptp+Lqp+phi=cptnk+Lqnk+gznk.
	388	!-------------------------------------------------------------------
	389
	390	do 320 i=1,len
	391	tnk(i)=t(i,nk(i))
	392	qnk(i)=q(i,nk(i))
	393	gznk(i)=gz(i,nk(i))
	394	c ori ticb(i)=t(i,icb(i))
	395	c ori gzicb(i)=gz(i,icb(i))
	396	320 continue
	397	c
	398	c * Calculate certain parcel quantities, including static energy *
	399	c
	400	do 330 i=1,len
	401	ah0(i)=(cpd(1.-qnk(i))+clqnk(i))*tnk(i)
	402	& +qnk(i)(lv0-clmcpv(tnk(i)-273.15))+gznk(i)
	403	cpp(i)=cpd(1.-qnk(i))+qnk(i)cpv
	404	cpinv(i)=1./cpp(i)
	405	330 continue
	406	c
	407	c * Calculate lifted parcel quantities below cloud base *
	408	c
	409	do i=1,len !convect3
	410	icb1(i)=MAX(icb(i),2) !convect3
	411	icb1(i)=MIN(icb(i),nl) !convect3
	412	c if icb is below LCL, start loop at ICB+1:
	413	c (icbs est le premier niveau au-dessus du LCL)
	414	icbs(i)=icb1(i) !convect3
	415	if (plcl(i).lt.p(i,icb1(i))) then
	416	icbs(i)=MIN(icbs(i)+1,nl) !convect3
	417	endif
	418	enddo !convect3
	419
	420	do i=1,len !convect3
	421	ticb(i)=t(i,icbs(i)) !convect3
	422	gzicb(i)=gz(i,icbs(i)) !convect3
	423	qsicb(i)=qs(i,icbs(i)) !convect3
	424	enddo !convect3
	425
	426	c
	427	c Re-compute icbsmax (icbsmax2): !convect3
	428	c !convect3
	429	icbsmax2=2 !convect3
	430	do 310 i=1,len !convect3
	431	icbsmax2=max(icbsmax2,icbs(i)) !convect3
	432	310 continue !convect3
	433
	434	c initialization outputs:
	435
	436	do k=1,icbsmax2 ! convect3
	437	do i=1,len ! convect3
	438	tp(i,k) = 0.0 ! convect3
	439	tvp(i,k) = 0.0 ! convect3
	440	clw(i,k) = 0.0 ! convect3
	441	enddo ! convect3
	442	enddo ! convect3
	443
	444	c tp and tvp below cloud base:
	445
	446	do 350 k=minorig,icbsmax2-1
	447	do 340 i=1,len
	448	tp(i,k)=tnk(i)-(gz(i,k)-gznk(i))*cpinv(i)
	449	tvp(i,k)=tp(i,k)*(1.+qnk(i)/eps-qnk(i)) !whole thing (convect3)
	450	340 continue
	451	350 continue
	452	c
	453	c * Find lifted parcel quantities above cloud base *
	454	c
	455	do 360 i=1,len
	456	tg=ticb(i)
	457	c ori qg=qs(i,icb(i))
	458	qg=qsicb(i) ! convect3
	459	cdebug alv=lv0-clmcpv*(ticb(i)-t0)
	460	alv=lv0-clmcpv*(ticb(i)-273.15)
	461	c
	462	c First iteration.
	463	c
	464	c ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	465	s=cpd(1.-qnk(i))+clqnk(i) ! convect3
	466	: +alvalvqg/(rrvticb(i)ticb(i)) ! convect3
	467	s=1./s
	468	c ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	469	ahg=cpdtg+(cl-cpd)qnk(i)tg+alvqg+gzicb(i) ! convect3
	470	tg=tg+s*(ah0(i)-ahg)
	471	c ori tg=max(tg,35.0)
	472	cdebug tc=tg-t0
	473	tc=tg-273.15
	474	denom=243.5+tc
	475	denom=MAX(denom,1.0) ! convect3
	476	c ori if(tc.ge.0.0)then
	477	es=6.112exp(17.67tc/denom)
	478	c ori else
	479	c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	480	c ori endif
	481	c ori qg=epses/(p(i,icb(i))-es(1.-eps))
	482	qg=epses/(p(i,icbs(i))-es(1.-eps))
	483	c
	484	c Second iteration.
	485	c
	486
	487	c ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	488	c ori s=1./s
	489	c ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	490	ahg=cpdtg+(cl-cpd)qnk(i)tg+alvqg+gzicb(i) ! convect3
	491	tg=tg+s*(ah0(i)-ahg)
	492	c ori tg=max(tg,35.0)
	493	cdebug tc=tg-t0
	494	tc=tg-273.15
	495	denom=243.5+tc
	496	denom=MAX(denom,1.0) ! convect3
	497	c ori if(tc.ge.0.0)then
	498	es=6.112exp(17.67tc/denom)
	499	c ori else
	500	c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	501	c ori end if
	502	c ori qg=epses/(p(i,icb(i))-es(1.-eps))
	503	qg=epses/(p(i,icbs(i))-es(1.-eps))
	504
	505	alv=lv0-clmcpv*(ticb(i)-273.15)
	506
	507	c ori c approximation here:
	508	c ori tp(i,icb(i))=(ah0(i)-(cl-cpd)qnk(i)ticb(i)
	509	c ori & -gz(i,icb(i))-alv*qg)/cpd
	510
	511	c convect3: no approximation:
	512	tp(i,icbs(i))=(ah0(i)-gz(i,icbs(i))-alv*qg)
	513	: /(cpd+(cl-cpd)*qnk(i))
	514
	515	c ori clw(i,icb(i))=qnk(i)-qg
	516	c ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
	517	clw(i,icbs(i))=qnk(i)-qg
	518	clw(i,icbs(i))=max(0.0,clw(i,icbs(i)))
	519
	520	rg=qg/(1.-qnk(i))
	521	c ori tvp(i,icb(i))=tp(i,icb(i))(1.+rgepsi)
	522	c convect3: (qg utilise au lieu du vrai mixing ratio rg)
	523	tvp(i,icbs(i))=tp(i,icbs(i))*(1.+qg/eps-qnk(i)) !whole thing
	524
	525	360 continue
	526	c
	527	c ori do 380 k=minorig,icbsmax2
	528	c ori do 370 i=1,len
	529	c ori tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
	530	c ori 370 continue
	531	c ori 380 continue
	532	c
	533
	534	c -- The following is only for convect3:
	535	c
	536	c * icbs is the first level above the LCL:
	537	c if plcl<p(icb), then icbs=icb+1
	538	c if plcl>p(icb), then icbs=icb
	539	c
	540	c * the routine above computes tvp from minorig to icbs (included).
	541	c
	542	c * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1)
	543	c must be known. This is the case if icbs=icb+1, but not if icbs=icb.
	544	c
	545	c * therefore, in the case icbs=icb, we compute tvp at level icb+1
	546	c (tvp at other levels will be computed in cv3_undilute2.F)
	547	c
	548
	549	do i=1,len
	550	ticb(i)=t(i,icb(i)+1)
	551	gzicb(i)=gz(i,icb(i)+1)
	552	qsicb(i)=qs(i,icb(i)+1)
	553	enddo
	554
	555	do 460 i=1,len
	556	tg=ticb(i)
	557	qg=qsicb(i) ! convect3
	558	cdebug alv=lv0-clmcpv*(ticb(i)-t0)
	559	alv=lv0-clmcpv*(ticb(i)-273.15)
	560	c
	561	c First iteration.
	562	c
	563	c ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	564	s=cpd(1.-qnk(i))+clqnk(i) ! convect3
	565	: +alvalvqg/(rrvticb(i)ticb(i)) ! convect3
	566	s=1./s
	567	c ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	568	ahg=cpdtg+(cl-cpd)qnk(i)tg+alvqg+gzicb(i) ! convect3
	569	tg=tg+s*(ah0(i)-ahg)
	570	c ori tg=max(tg,35.0)
	571	cdebug tc=tg-t0
	572	tc=tg-273.15
	573	denom=243.5+tc
	574	denom=MAX(denom,1.0) ! convect3
	575	c ori if(tc.ge.0.0)then
	576	es=6.112exp(17.67tc/denom)
	577	c ori else
	578	c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	579	c ori endif
	580	c ori qg=epses/(p(i,icb(i))-es(1.-eps))
	581	qg=epses/(p(i,icb(i)+1)-es(1.-eps))
	582	c
	583	c Second iteration.
	584	c
	585
	586	c ori s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	587	c ori s=1./s
	588	c ori ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	589	ahg=cpdtg+(cl-cpd)qnk(i)tg+alvqg+gzicb(i) ! convect3
	590	tg=tg+s*(ah0(i)-ahg)
	591	c ori tg=max(tg,35.0)
	592	cdebug tc=tg-t0
	593	tc=tg-273.15
	594	denom=243.5+tc
	595	denom=MAX(denom,1.0) ! convect3
	596	c ori if(tc.ge.0.0)then
	597	es=6.112exp(17.67tc/denom)
	598	c ori else
	599	c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	600	c ori end if
	601	c ori qg=epses/(p(i,icb(i))-es(1.-eps))
	602	qg=epses/(p(i,icb(i)+1)-es(1.-eps))
	603
	604	alv=lv0-clmcpv*(ticb(i)-273.15)
	605
	606	c ori c approximation here:
	607	c ori tp(i,icb(i))=(ah0(i)-(cl-cpd)qnk(i)ticb(i)
	608	c ori & -gz(i,icb(i))-alv*qg)/cpd
	609
	610	c convect3: no approximation:
	611	tp(i,icb(i)+1)=(ah0(i)-gz(i,icb(i)+1)-alv*qg)
	612	: /(cpd+(cl-cpd)*qnk(i))
	613
	614	c ori clw(i,icb(i))=qnk(i)-qg
	615	c ori clw(i,icb(i))=max(0.0,clw(i,icb(i)))
	616	clw(i,icb(i)+1)=qnk(i)-qg
	617	clw(i,icb(i)+1)=max(0.0,clw(i,icb(i)+1))
	618
	619	rg=qg/(1.-qnk(i))
	620	c ori tvp(i,icb(i))=tp(i,icb(i))(1.+rgepsi)
	621	c convect3: (qg utilise au lieu du vrai mixing ratio rg)
	622	tvp(i,icb(i)+1)=tp(i,icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
	623
	624	460 continue
	625
	626	return
	627	end
	628
	629	SUBROUTINE cv30_trigger(len,nd,icb,plcl,p,th,tv,tvp
	630	o ,pbase,buoybase,iflag,sig,w0)
	631	implicit none
	632
	633	!-------------------------------------------------------------------
	634	! --- TRIGGERING
	635	!
	636	! - computes the cloud base
	637	! - triggering (crude in this version)
	638	! - relaxation of sig and w0 when no convection
	639	!
	640	! Caution1: if no convection, we set iflag=4
	641	! (it used to be 0 in convect3)
	642	!
	643	! Caution2: at this stage, tvp (and thus buoy) are know up
	644	! through icb only!
	645	! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy
	646	!-------------------------------------------------------------------
	647
	648	#include "cv30param.h"
	649
	650	c input:
	651	integer len, nd
	652	integer icb(len)
	653	real plcl(len), p(len,nd)
	654	real th(len,nd), tv(len,nd), tvp(len,nd)
	655
	656	c output:
	657	real pbase(len), buoybase(len)
	658
	659	c input AND output:
	660	integer iflag(len)
	661	real sig(len,nd), w0(len,nd)
	662
	663	c local variables:
	664	integer i,k
	665	real tvpbase, tvbase, tdif, ath, ath1
	666
	667	c
	668	c *** set cloud base buoyancy at (plcl+dpbase) level buoyancy
	669	c
	670	do 100 i=1,len
	671	pbase(i) = plcl(i) + dpbase
	672	tvpbase = tvp(i,icb(i))*(pbase(i)-p(i,icb(i)+1))
	673	: /(p(i,icb(i))-p(i,icb(i)+1))
	674	: + tvp(i,icb(i)+1)*(p(i,icb(i))-pbase(i))
	675	: /(p(i,icb(i))-p(i,icb(i)+1))
	676	tvbase = tv(i,icb(i))*(pbase(i)-p(i,icb(i)+1))
	677	: /(p(i,icb(i))-p(i,icb(i)+1))
	678	: + tv(i,icb(i)+1)*(p(i,icb(i))-pbase(i))
	679	: /(p(i,icb(i))-p(i,icb(i)+1))
	680	buoybase(i) = tvpbase - tvbase
	681	100 continue
	682
	683	c
	684	c * make sure that column is dry adiabatic between the surface *
	685	c * and cloud base, and that lifted air is positively buoyant *
	686	c * at cloud base *
	687	c * if not, return to calling program after resetting *
	688	c * sig(i) and w0(i) *
	689	c
	690
	691	c oct3 do 200 i=1,len
	692	c oct3
	693	c oct3 tdif = buoybase(i)
	694	c oct3 ath1 = th(i,1)
	695	c oct3 ath = th(i,icb(i)-1) - dttrig
	696	c oct3
	697	c oct3 if (tdif.lt.dtcrit .or. ath.gt.ath1) then
	698	c oct3 do 60 k=1,nl
	699	c oct3 sig(i,k) = betasig(i,k) - 2.alphatdiftdif
	700	c oct3 sig(i,k) = AMAX1(sig(i,k),0.0)
	701	c oct3 w0(i,k) = beta*w0(i,k)
	702	c oct3 60 continue
	703	c oct3 iflag(i)=4 ! pour version vectorisee
	704	c oct3c convect3 iflag(i)=0
	705	c oct3cccc return
	706	c oct3 endif
	707	c oct3
	708	c oct3200 continue
	709
	710	c -- oct3: on reecrit la boucle 200 (pour la vectorisation)
	711
	712	do 60 k=1,nl
	713	do 200 i=1,len
	714
	715	tdif = buoybase(i)
	716	ath1 = th(i,1)
	717	ath = th(i,icb(i)-1) - dttrig
	718
	719	if (tdif.lt.dtcrit .or. ath.gt.ath1) then
	720	sig(i,k) = betasig(i,k) - 2.alphatdiftdif
	721	sig(i,k) = AMAX1(sig(i,k),0.0)
	722	w0(i,k) = beta*w0(i,k)
	723	iflag(i)=4 ! pour version vectorisee
	724	c convect3 iflag(i)=0
	725	endif
	726
	727	200 continue
	728	60 continue
	729
	730	c fin oct3 --
	731
	732	return
	733	end
	734
	735	SUBROUTINE cv30_compress( len,nloc,ncum,nd,ntra
	736	: ,iflag1,nk1,icb1,icbs1
	737	: ,plcl1,tnk1,qnk1,gznk1,pbase1,buoybase1
	738	: ,t1,q1,qs1,u1,v1,gz1,th1
	739	: ,tra1
	740	: ,h1,lv1,cpn1,p1,ph1,tv1,tp1,tvp1,clw1
	741	: ,sig1,w01
	742	o ,iflag,nk,icb,icbs
	743	o ,plcl,tnk,qnk,gznk,pbase,buoybase
	744	o ,t,q,qs,u,v,gz,th
	745	o ,tra
	746	o ,h,lv,cpn,p,ph,tv,tp,tvp,clw
	747	o ,sig,w0 )
	748	implicit none
	749
	750	#include "cv30param.h"
[1403]	751	include 'iniprint.h'
[879]	752
	753	c inputs:
	754	integer len,ncum,nd,ntra,nloc
	755	integer iflag1(len),nk1(len),icb1(len),icbs1(len)
	756	real plcl1(len),tnk1(len),qnk1(len),gznk1(len)
	757	real pbase1(len),buoybase1(len)
	758	real t1(len,nd),q1(len,nd),qs1(len,nd),u1(len,nd),v1(len,nd)
	759	real gz1(len,nd),h1(len,nd),lv1(len,nd),cpn1(len,nd)
	760	real p1(len,nd),ph1(len,nd+1),tv1(len,nd),tp1(len,nd)
	761	real tvp1(len,nd),clw1(len,nd)
	762	real th1(len,nd)
	763	real sig1(len,nd), w01(len,nd)
	764	real tra1(len,nd,ntra)
	765
	766	c outputs:
	767	c en fait, on a nloc=len pour l'instant (cf cv_driver)
	768	integer iflag(nloc),nk(nloc),icb(nloc),icbs(nloc)
	769	real plcl(nloc),tnk(nloc),qnk(nloc),gznk(nloc)
	770	real pbase(nloc),buoybase(nloc)
	771	real t(nloc,nd),q(nloc,nd),qs(nloc,nd),u(nloc,nd),v(nloc,nd)
	772	real gz(nloc,nd),h(nloc,nd),lv(nloc,nd),cpn(nloc,nd)
	773	real p(nloc,nd),ph(nloc,nd+1),tv(nloc,nd),tp(nloc,nd)
	774	real tvp(nloc,nd),clw(nloc,nd)
	775	real th(nloc,nd)
	776	real sig(nloc,nd), w0(nloc,nd)
	777	real tra(nloc,nd,ntra)
	778
	779	c local variables:
	780	integer i,k,nn,j
	781
[1403]	782	CHARACTER (LEN=20) :: modname='cv30_compress'
	783	CHARACTER (LEN=80) :: abort_message
[879]	784
[1403]	785
[879]	786	do 110 k=1,nl+1
	787	nn=0
	788	do 100 i=1,len
	789	if(iflag1(i).eq.0)then
	790	nn=nn+1
	791	sig(nn,k)=sig1(i,k)
	792	w0(nn,k)=w01(i,k)
	793	t(nn,k)=t1(i,k)
	794	q(nn,k)=q1(i,k)
	795	qs(nn,k)=qs1(i,k)
	796	u(nn,k)=u1(i,k)
	797	v(nn,k)=v1(i,k)
	798	gz(nn,k)=gz1(i,k)
	799	h(nn,k)=h1(i,k)
	800	lv(nn,k)=lv1(i,k)
	801	cpn(nn,k)=cpn1(i,k)
	802	p(nn,k)=p1(i,k)
	803	ph(nn,k)=ph1(i,k)
	804	tv(nn,k)=tv1(i,k)
	805	tp(nn,k)=tp1(i,k)
	806	tvp(nn,k)=tvp1(i,k)
	807	clw(nn,k)=clw1(i,k)
	808	th(nn,k)=th1(i,k)
	809	endif
	810	100 continue
	811	110 continue
	812
	813	c do 121 j=1,ntra
	814	c do 111 k=1,nd
	815	c nn=0
	816	c do 101 i=1,len
	817	c if(iflag1(i).eq.0)then
	818	c nn=nn+1
	819	c tra(nn,k,j)=tra1(i,k,j)
	820	c endif
	821	c 101 continue
	822	c 111 continue
	823	c 121 continue
	824
	825	if (nn.ne.ncum) then
[1403]	826	write(lunout,*)'strange! nn not equal to ncum: ',nn,ncum
	827	abort_message = ''
	828	CALL abort_gcm (modname,abort_message,1)
[879]	829	endif
	830
	831	nn=0
	832	do 150 i=1,len
	833	if(iflag1(i).eq.0)then
	834	nn=nn+1
	835	pbase(nn)=pbase1(i)
	836	buoybase(nn)=buoybase1(i)
	837	plcl(nn)=plcl1(i)
	838	tnk(nn)=tnk1(i)
	839	qnk(nn)=qnk1(i)
	840	gznk(nn)=gznk1(i)
	841	nk(nn)=nk1(i)
	842	icb(nn)=icb1(i)
	843	icbs(nn)=icbs1(i)
	844	iflag(nn)=iflag1(i)
	845	endif
	846	150 continue
	847
	848	return
	849	end
	850
	851	SUBROUTINE cv30_undilute2(nloc,ncum,nd,icb,icbs,nk
	852	: ,tnk,qnk,gznk,t,q,qs,gz
	853	: ,p,h,tv,lv,pbase,buoybase,plcl
	854	o ,inb,tp,tvp,clw,hp,ep,sigp,buoy)
	855	implicit none
	856
	857	C---------------------------------------------------------------------
	858	C Purpose:
	859	C FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
	860	C &
	861	C COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
	862	C FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
	863	C &
	864	C FIND THE LEVEL OF NEUTRAL BUOYANCY
	865	C
	866	C Main differences convect3/convect4:
	867	C - icbs (input) is the first level above LCL (may differ from icb)
	868	C - many minor differences in the iterations
	869	C - condensed water not removed from tvp in convect3
	870	C - vertical profile of buoyancy computed here (use of buoybase)
	871	C - the determination of inb is different
	872	C - no inb1, only inb in output
	873	C---------------------------------------------------------------------
	874
	875	#include "cvthermo.h"
	876	#include "cv30param.h"
	877	#include "conema3.h"
	878
	879	c inputs:
	880	integer ncum, nd, nloc
	881	integer icb(nloc), icbs(nloc), nk(nloc)
	882	real t(nloc,nd), q(nloc,nd), qs(nloc,nd), gz(nloc,nd)
	883	real p(nloc,nd)
	884	real tnk(nloc), qnk(nloc), gznk(nloc)
	885	real lv(nloc,nd), tv(nloc,nd), h(nloc,nd)
	886	real pbase(nloc), buoybase(nloc), plcl(nloc)
	887
	888	c outputs:
	889	integer inb(nloc)
	890	real tp(nloc,nd), tvp(nloc,nd), clw(nloc,nd)
	891	real ep(nloc,nd), sigp(nloc,nd), hp(nloc,nd)
	892	real buoy(nloc,nd)
	893
	894	c local variables:
	895	integer i, k
	896	real tg,qg,ahg,alv,s,tc,es,denom,rg,tca,elacrit
	897	real by, defrac, pden
	898	real ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
	899	logical lcape(nloc)
	900
	901	!=====================================================================
	902	! --- SOME INITIALIZATIONS
	903	!=====================================================================
	904
	905	do 170 k=1,nl
	906	do 160 i=1,ncum
	907	ep(i,k)=0.0
	908	sigp(i,k)=spfac
	909	160 continue
	910	170 continue
	911
	912	!=====================================================================
	913	! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
	914	!=====================================================================
	915	c
	916	c --- The procedure is to solve the equation.
	917	c cptp+Lqp+phi=cptnk+Lqnk+gznk.
	918	c
	919	c * Calculate certain parcel quantities, including static energy *
	920	c
	921	c
	922	do 240 i=1,ncum
	923	ah0(i)=(cpd(1.-qnk(i))+clqnk(i))*tnk(i)
	924	cdebug & +qnk(i)(lv0-clmcpv(tnk(i)-t0))+gznk(i)
	925	& +qnk(i)(lv0-clmcpv(tnk(i)-273.15))+gznk(i)
	926	240 continue
	927	c
	928	c
	929	c * Find lifted parcel quantities above cloud base *
	930	c
	931	c
	932	do 300 k=minorig+1,nl
	933	do 290 i=1,ncum
	934	c ori if(k.ge.(icb(i)+1))then
	935	if(k.ge.(icbs(i)+1))then ! convect3
	936	tg=t(i,k)
	937	qg=qs(i,k)
	938	cdebug alv=lv0-clmcpv*(t(i,k)-t0)
	939	alv=lv0-clmcpv*(t(i,k)-273.15)
	940	c
	941	c First iteration.
	942	c
	943	c ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
	944	s=cpd(1.-qnk(i))+clqnk(i) ! convect3
	945	: +alvalvqg/(rrvt(i,k)t(i,k)) ! convect3
	946	s=1./s
	947	c ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
	948	ahg=cpdtg+(cl-cpd)qnk(i)tg+alvqg+gz(i,k) ! convect3
	949	tg=tg+s*(ah0(i)-ahg)
	950	c ori tg=max(tg,35.0)
	951	cdebug tc=tg-t0
	952	tc=tg-273.15
	953	denom=243.5+tc
	954	denom=MAX(denom,1.0) ! convect3
	955	c ori if(tc.ge.0.0)then
	956	es=6.112exp(17.67tc/denom)
	957	c ori else
	958	c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	959	c ori endif
	960	qg=epses/(p(i,k)-es(1.-eps))
	961	c
	962	c Second iteration.
	963	c
	964	c ori s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
	965	c ori s=1./s
	966	c ori ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
	967	ahg=cpdtg+(cl-cpd)qnk(i)tg+alvqg+gz(i,k) ! convect3
	968	tg=tg+s*(ah0(i)-ahg)
	969	c ori tg=max(tg,35.0)
	970	cdebug tc=tg-t0
	971	tc=tg-273.15
	972	denom=243.5+tc
	973	denom=MAX(denom,1.0) ! convect3
	974	c ori if(tc.ge.0.0)then
	975	es=6.112exp(17.67tc/denom)
	976	c ori else
	977	c ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	978	c ori endif
	979	qg=epses/(p(i,k)-es(1.-eps))
	980	c
	981	cdebug alv=lv0-clmcpv*(t(i,k)-t0)
	982	alv=lv0-clmcpv*(t(i,k)-273.15)
	983	c print*,'cpd dans convect2 ',cpd
	984	c print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
	985	c print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
	986
	987	c ori c approximation here:
	988	c ori tp(i,k)=(ah0(i)-(cl-cpd)qnk(i)t(i,k)-gz(i,k)-alv*qg)/cpd
	989
	990	c convect3: no approximation:
	991	tp(i,k)=(ah0(i)-gz(i,k)-alvqg)/(cpd+(cl-cpd)qnk(i))
	992
	993	clw(i,k)=qnk(i)-qg
	994	clw(i,k)=max(0.0,clw(i,k))
	995	rg=qg/(1.-qnk(i))
	996	c ori tvp(i,k)=tp(i,k)(1.+rgepsi)
	997	c convect3: (qg utilise au lieu du vrai mixing ratio rg):
	998	tvp(i,k)=tp(i,k)*(1.+qg/eps-qnk(i)) ! whole thing
	999	endif
	1000	290 continue
	1001	300 continue
	1002	c
	1003	!=====================================================================
	1004	! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
	1005	! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
	1006	! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
	1007	!=====================================================================
	1008	c
	1009	c ori do 320 k=minorig+1,nl
	1010	do 320 k=1,nl ! convect3
	1011	do 310 i=1,ncum
	1012	pden=ptcrit-pbcrit
	1013	ep(i,k)=(plcl(i)-p(i,k)-pbcrit)/pden*epmax
	1014	ep(i,k)=amax1(ep(i,k),0.0)
	1015	ep(i,k)=amin1(ep(i,k),epmax)
	1016	sigp(i,k)=spfac
	1017	c ori if(k.ge.(nk(i)+1))then
	1018	c ori tca=tp(i,k)-t0
	1019	c ori if(tca.ge.0.0)then
	1020	c ori elacrit=elcrit
	1021	c ori else
	1022	c ori elacrit=elcrit*(1.0-tca/tlcrit)
	1023	c ori endif
	1024	c ori elacrit=max(elacrit,0.0)
	1025	c ori ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)
	1026	c ori ep(i,k)=max(ep(i,k),0.0 )
	1027	c ori ep(i,k)=min(ep(i,k),1.0 )
	1028	c ori sigp(i,k)=sigs
	1029	c ori endif
	1030	310 continue
	1031	320 continue
	1032	c
	1033	!=====================================================================
	1034	! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
	1035	! --- VIRTUAL TEMPERATURE
	1036	!=====================================================================
	1037	c
	1038	c dans convect3, tvp est calcule en une seule fois, et sans retirer
	1039	c l'eau condensee (~> reversible CAPE)
	1040	c
	1041	c ori do 340 k=minorig+1,nl
	1042	c ori do 330 i=1,ncum
	1043	c ori if(k.ge.(icb(i)+1))then
	1044	c ori tvp(i,k)=tvp(i,k)(1.0-qnk(i)+ep(i,k)clw(i,k))
	1045	c oric print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
	1046	c oric print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
	1047	c ori endif
	1048	c ori 330 continue
	1049	c ori 340 continue
	1050
	1051	c ori do 350 i=1,ncum
	1052	c ori tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
	1053	c ori 350 continue
	1054
	1055	do 350 i=1,ncum ! convect3
	1056	tp(i,nlp)=tp(i,nl) ! convect3
	1057	350 continue ! convect3
	1058	c
	1059	c=====================================================================
	1060	c --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
	1061	c=====================================================================
	1062
	1063	c-- this is for convect3 only:
	1064
	1065	c first estimate of buoyancy:
	1066
	1067	do 500 i=1,ncum
	1068	do 501 k=1,nl
	1069	buoy(i,k)=tvp(i,k)-tv(i,k)
	1070	501 continue
	1071	500 continue
	1072
	1073	c set buoyancy=buoybase for all levels below base
	1074	c for safety, set buoy(icb)=buoybase
	1075
	1076	do 505 i=1,ncum
	1077	do 506 k=1,nl
	1078	if((k.ge.icb(i)).and.(k.le.nl).and.(p(i,k).ge.pbase(i)))then
	1079	buoy(i,k)=buoybase(i)
	1080	endif
	1081	506 continue
	1082	cIM cf. CRio/JYG 270807 buoy(icb(i),k)=buoybase(i)
	1083	buoy(i,icb(i))=buoybase(i)
	1084	505 continue
	1085
	1086	c-- end convect3
	1087
	1088	c=====================================================================
	1089	c --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
	1090	c --- LEVEL OF NEUTRAL BUOYANCY
	1091	c=====================================================================
	1092	c
	1093	c-- this is for convect3 only:
	1094
	1095	do 510 i=1,ncum
	1096	inb(i)=nl-1
	1097	510 continue
	1098
	1099	do 530 i=1,ncum
	1100	do 535 k=1,nl-1
	1101	if ((k.ge.icb(i)).and.(buoy(i,k).lt.dtovsh)) then
	1102	inb(i)=MIN(inb(i),k)
	1103	endif
	1104	535 continue
	1105	530 continue
	1106
	1107	c-- end convect3
	1108
	1109	c ori do 510 i=1,ncum
	1110	c ori cape(i)=0.0
	1111	c ori capem(i)=0.0
	1112	c ori inb(i)=icb(i)+1
	1113	c ori inb1(i)=inb(i)
	1114	c ori 510 continue
	1115	c
	1116	c Originial Code
	1117	c
	1118	c do 530 k=minorig+1,nl-1
	1119	c do 520 i=1,ncum
	1120	c if(k.ge.(icb(i)+1))then
	1121	c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	1122	c byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	1123	c cape(i)=cape(i)+by
	1124	c if(by.ge.0.0)inb1(i)=k+1
	1125	c if(cape(i).gt.0.0)then
	1126	c inb(i)=k+1
	1127	c capem(i)=cape(i)
	1128	c endif
	1129	c endif
	1130	c520 continue
	1131	c530 continue
	1132	c do 540 i=1,ncum
	1133	c byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
	1134	c cape(i)=capem(i)+byp
	1135	c defrac=capem(i)-cape(i)
	1136	c defrac=max(defrac,0.001)
	1137	c frac(i)=-cape(i)/defrac
	1138	c frac(i)=min(frac(i),1.0)
	1139	c frac(i)=max(frac(i),0.0)
	1140	c540 continue
	1141	c
	1142	c K Emanuel fix
	1143	c
	1144	c call zilch(byp,ncum)
	1145	c do 530 k=minorig+1,nl-1
	1146	c do 520 i=1,ncum
	1147	c if(k.ge.(icb(i)+1))then
	1148	c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	1149	c cape(i)=cape(i)+by
	1150	c if(by.ge.0.0)inb1(i)=k+1
	1151	c if(cape(i).gt.0.0)then
	1152	c inb(i)=k+1
	1153	c capem(i)=cape(i)
	1154	c byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	1155	c endif
	1156	c endif
	1157	c520 continue
	1158	c530 continue
	1159	c do 540 i=1,ncum
	1160	c inb(i)=max(inb(i),inb1(i))
	1161	c cape(i)=capem(i)+byp(i)
	1162	c defrac=capem(i)-cape(i)
	1163	c defrac=max(defrac,0.001)
	1164	c frac(i)=-cape(i)/defrac
	1165	c frac(i)=min(frac(i),1.0)
	1166	c frac(i)=max(frac(i),0.0)
	1167	c540 continue
	1168	c
	1169	c J Teixeira fix
	1170	c
	1171	c ori call zilch(byp,ncum)
	1172	c ori do 515 i=1,ncum
	1173	c ori lcape(i)=.true.
	1174	c ori 515 continue
	1175	c ori do 530 k=minorig+1,nl-1
	1176	c ori do 520 i=1,ncum
	1177	c ori if(cape(i).lt.0.0)lcape(i)=.false.
	1178	c ori if((k.ge.(icb(i)+1)).and.lcape(i))then
	1179	c ori by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	1180	c ori byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	1181	c ori cape(i)=cape(i)+by
	1182	c ori if(by.ge.0.0)inb1(i)=k+1
	1183	c ori if(cape(i).gt.0.0)then
	1184	c ori inb(i)=k+1
	1185	c ori capem(i)=cape(i)
	1186	c ori endif
	1187	c ori endif
	1188	c ori 520 continue
	1189	c ori 530 continue
	1190	c ori do 540 i=1,ncum
	1191	c ori cape(i)=capem(i)+byp(i)
	1192	c ori defrac=capem(i)-cape(i)
	1193	c ori defrac=max(defrac,0.001)
	1194	c ori frac(i)=-cape(i)/defrac
	1195	c ori frac(i)=min(frac(i),1.0)
	1196	c ori frac(i)=max(frac(i),0.0)
	1197	c ori 540 continue
	1198	c
	1199	c=====================================================================
	1200	c --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
	1201	c=====================================================================
	1202	c
	1203	cym do i=1,ncum*nlp
	1204	cym hp(i,1)=h(i,1)
	1205	cym enddo
	1206
	1207	do k=1,nlp
	1208	do i=1,ncum
	1209	hp(i,k)=h(i,k)
	1210	enddo
	1211	enddo
	1212
	1213	do 600 k=minorig+1,nl
	1214	do 590 i=1,ncum
	1215	if((k.ge.icb(i)).and.(k.le.inb(i)))then
	1216	hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)t(i,k))ep(i,k)*clw(i,k)
	1217	endif
	1218	590 continue
	1219	600 continue
	1220
	1221	return
	1222	end
	1223
	1224	SUBROUTINE cv30_closure(nloc,ncum,nd,icb,inb
	1225	: ,pbase,p,ph,tv,buoy
	1226	o ,sig,w0,cape,m)
	1227	implicit none
	1228
	1229	!===================================================================
	1230	! --- CLOSURE OF CONVECT3
	1231	!
	1232	! vectorization: S. Bony
	1233	!===================================================================
	1234
	1235	#include "cvthermo.h"
	1236	#include "cv30param.h"
	1237
	1238	c input:
	1239	integer ncum, nd, nloc
	1240	integer icb(nloc), inb(nloc)
	1241	real pbase(nloc)
	1242	real p(nloc,nd), ph(nloc,nd+1)
	1243	real tv(nloc,nd), buoy(nloc,nd)
	1244
	1245	c input/output:
	1246	real sig(nloc,nd), w0(nloc,nd)
	1247
	1248	c output:
	1249	real cape(nloc)
	1250	real m(nloc,nd)
	1251
	1252	c local variables:
	1253	integer i, j, k, icbmax
	1254	real deltap, fac, w, amu
	1255	real dtmin(nloc,nd), sigold(nloc,nd)
	1256
	1257
	1258	c -------------------------------------------------------
	1259	c -- Initialization
	1260	c -------------------------------------------------------
	1261
	1262	do k=1,nl
	1263	do i=1,ncum
	1264	m(i,k)=0.0
	1265	enddo
	1266	enddo
	1267
	1268	c -------------------------------------------------------
	1269	c -- Reset sig(i) and w0(i) for i>inb and i<icb
	1270	c -------------------------------------------------------
	1271
	1272	c update sig and w0 above LNB:
	1273
	1274	do 100 k=1,nl-1
	1275	do 110 i=1,ncum
	1276	if ((inb(i).lt.(nl-1)).and.(k.ge.(inb(i)+1)))then
	1277	sig(i,k)=beta*sig(i,k)
	1278	: +2.alphabuoy(i,inb(i))*ABS(buoy(i,inb(i)))
	1279	sig(i,k)=AMAX1(sig(i,k),0.0)
	1280	w0(i,k)=beta*w0(i,k)
	1281	endif
	1282	110 continue
	1283	100 continue
	1284
	1285	c compute icbmax:
	1286
	1287	icbmax=2
	1288	do 200 i=1,ncum
	1289	icbmax=MAX(icbmax,icb(i))
	1290	200 continue
	1291
	1292	c update sig and w0 below cloud base:
	1293
	1294	do 300 k=1,icbmax
	1295	do 310 i=1,ncum
	1296	if (k.le.icb(i))then
	1297	sig(i,k)=betasig(i,k)-2.alphabuoy(i,icb(i))buoy(i,icb(i))
	1298	sig(i,k)=amax1(sig(i,k),0.0)
	1299	w0(i,k)=beta*w0(i,k)
	1300	endif
	1301	310 continue
	1302	300 continue
	1303
	1304	c! if(inb.lt.(nl-1))then
	1305	c! do 85 i=inb+1,nl-1
	1306	c! sig(i)=betasig(i)+2.alphabuoy(inb)
	1307	c! 1 abs(buoy(inb))
	1308	c! sig(i)=amax1(sig(i),0.0)
	1309	c! w0(i)=beta*w0(i)
	1310	c! 85 continue
	1311	c! end if
	1312
	1313	c! do 87 i=1,icb
	1314	c! sig(i)=betasig(i)-2.alphabuoy(icb)buoy(icb)
	1315	c! sig(i)=amax1(sig(i),0.0)
	1316	c! w0(i)=beta*w0(i)
	1317	c! 87 continue
	1318
	1319	c -------------------------------------------------------------
	1320	c -- Reset fractional areas of updrafts and w0 at initial time
	1321	c -- and after 10 time steps of no convection
	1322	c -------------------------------------------------------------
	1323
	1324	do 400 k=1,nl-1
	1325	do 410 i=1,ncum
	1326	if (sig(i,nd).lt.1.5.or.sig(i,nd).gt.12.0)then
	1327	sig(i,k)=0.0
	1328	w0(i,k)=0.0
	1329	endif
	1330	410 continue
	1331	400 continue
	1332
	1333	c -------------------------------------------------------------
	1334	c -- Calculate convective available potential energy (cape),
	1335	c -- vertical velocity (w), fractional area covered by
	1336	c -- undilute updraft (sig), and updraft mass flux (m)
	1337	c -------------------------------------------------------------
	1338
	1339	do 500 i=1,ncum
	1340	cape(i)=0.0
	1341	500 continue
	1342
	1343	c compute dtmin (minimum buoyancy between ICB and given level k):
	1344
	1345	do i=1,ncum
	1346	do k=1,nl
	1347	dtmin(i,k)=100.0
	1348	enddo
	1349	enddo
	1350
	1351	do 550 i=1,ncum
	1352	do 560 k=1,nl
	1353	do 570 j=minorig,nl
	1354	if ( (k.ge.(icb(i)+1)).and.(k.le.inb(i)).and.
	1355	: (j.ge.icb(i)).and.(j.le.(k-1)) )then
	1356	dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
	1357	endif
	1358	570 continue
	1359	560 continue
	1360	550 continue
	1361
	1362	c the interval on which cape is computed starts at pbase :
	1363
	1364	do 600 k=1,nl
	1365	do 610 i=1,ncum
	1366
	1367	if ((k.ge.(icb(i)+1)).and.(k.le.inb(i))) then
	1368
	1369	deltap = MIN(pbase(i),ph(i,k-1))-MIN(pbase(i),ph(i,k))
	1370	cape(i)=cape(i)+rrdbuoy(i,k-1)deltap/p(i,k-1)
	1371	cape(i)=AMAX1(0.0,cape(i))
	1372	sigold(i,k)=sig(i,k)
	1373
	1374	c dtmin(i,k)=100.0
	1375	c do 97 j=icb(i),k-1 ! mauvaise vectorisation
	1376	c dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
	1377	c 97 continue
	1378
	1379	sig(i,k)=betasig(i,k)+alphadtmin(i,k)*ABS(dtmin(i,k))
	1380	sig(i,k)=amax1(sig(i,k),0.0)
	1381	sig(i,k)=amin1(sig(i,k),0.01)
	1382	fac=AMIN1(((dtcrit-dtmin(i,k))/dtcrit),1.0)
	1383	w=(1.-beta)facSQRT(cape(i))+beta*w0(i,k)
	1384	amu=0.5(sig(i,k)+sigold(i,k))w
	1385	m(i,k)=amu0.007p(i,k)*(ph(i,k)-ph(i,k+1))/tv(i,k)
	1386	w0(i,k)=w
	1387	endif
	1388
	1389	610 continue
	1390	600 continue
	1391
	1392	do 700 i=1,ncum
	1393	w0(i,icb(i))=0.5*w0(i,icb(i)+1)
	1394	m(i,icb(i))=0.5*m(i,icb(i)+1)
	1395	: *(ph(i,icb(i))-ph(i,icb(i)+1))
	1396	: /(ph(i,icb(i)+1)-ph(i,icb(i)+2))
	1397	sig(i,icb(i))=sig(i,icb(i)+1)
	1398	sig(i,icb(i)-1)=sig(i,icb(i))
	1399	700 continue
	1400
	1401
	1402	c! cape=0.0
	1403	c! do 98 i=icb+1,inb
	1404	c! deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
	1405	c! cape=cape+rrdbuoy(i-1)deltap/p(i-1)
	1406	c! dcape=rrdbuoy(i-1)deltap/p(i-1)
	1407	c! dlnp=deltap/p(i-1)
	1408	c! cape=amax1(0.0,cape)
	1409	c! sigold=sig(i)
	1410
	1411	c! dtmin=100.0
	1412	c! do 97 j=icb,i-1
	1413	c! dtmin=amin1(dtmin,buoy(j))
	1414	c! 97 continue
	1415
	1416	c! sig(i)=betasig(i)+alphadtmin*abs(dtmin)
	1417	c! sig(i)=amax1(sig(i),0.0)
	1418	c! sig(i)=amin1(sig(i),0.01)
	1419	c! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
	1420	c! w=(1.-beta)facsqrt(cape)+beta*w0(i)
	1421	c! amu=0.5(sig(i)+sigold)w
	1422	c! m(i)=amu0.007p(i)*(ph(i)-ph(i+1))/tv(i)
	1423	c! w0(i)=w
	1424	c! 98 continue
	1425	c! w0(icb)=0.5*w0(icb+1)
	1426	c! m(icb)=0.5m(icb+1)(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
	1427	c! sig(icb)=sig(icb+1)
	1428	c! sig(icb-1)=sig(icb)
	1429
	1430	return
	1431	end
	1432
	1433	SUBROUTINE cv30_mixing(nloc,ncum,nd,na,ntra,icb,nk,inb
	1434	: ,ph,t,rr,rs,u,v,tra,h,lv,qnk
	1435	: ,hp,tv,tvp,ep,clw,m,sig
	1436	: ,ment,qent,uent,vent,sij,elij,ments,qents,traent)
	1437	implicit none
	1438
	1439	!---------------------------------------------------------------------
	1440	! a faire:
	1441	! - changer rr(il,1) -> qnk(il)
	1442	! - vectorisation de la partie normalisation des flux (do 789...)
	1443	!---------------------------------------------------------------------
	1444
	1445	#include "cvthermo.h"
	1446	#include "cv30param.h"
	1447
	1448	c inputs:
	1449	integer ncum, nd, na, ntra, nloc
	1450	integer icb(nloc), inb(nloc), nk(nloc)
	1451	real sig(nloc,nd)
	1452	real qnk(nloc)
	1453	real ph(nloc,nd+1)
	1454	real t(nloc,nd), rr(nloc,nd), rs(nloc,nd)
	1455	real u(nloc,nd), v(nloc,nd)
	1456	real tra(nloc,nd,ntra) ! input of convect3
	1457	real lv(nloc,na), h(nloc,na), hp(nloc,na)
	1458	real tv(nloc,na), tvp(nloc,na), ep(nloc,na), clw(nloc,na)
	1459	real m(nloc,na) ! input of convect3
	1460
	1461	c outputs:
	1462	real ment(nloc,na,na), qent(nloc,na,na)
	1463	real uent(nloc,na,na), vent(nloc,na,na)
	1464	real sij(nloc,na,na), elij(nloc,na,na)
	1465	real traent(nloc,nd,nd,ntra)
	1466	real ments(nloc,nd,nd), qents(nloc,nd,nd)
	1467	real sigij(nloc,nd,nd)
	1468
	1469	c local variables:
	1470	integer i, j, k, il, im, jm
	1471	integer num1, num2
	1472	integer nent(nloc,na)
	1473	real rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
	1474	real alt, smid, sjmin, sjmax, delp, delm
	1475	real asij(nloc), smax(nloc), scrit(nloc)
	1476	real asum(nloc,nd),bsum(nloc,nd),csum(nloc,nd)
	1477	real wgh
	1478	real zm(nloc,na)
	1479	logical lwork(nloc)
	1480
	1481	c=====================================================================
	1482	c --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
	1483	c=====================================================================
	1484
	1485	c ori do 360 i=1,ncum*nlp
	1486	do 361 j=1,nl
	1487	do 360 i=1,ncum
	1488	nent(i,j)=0
	1489	c in convect3, m is computed in cv3_closure
	1490	c ori m(i,1)=0.0
	1491	360 continue
	1492	361 continue
	1493
	1494	c ori do 400 k=1,nlp
	1495	c ori do 390 j=1,nlp
	1496	do 400 j=1,nl
	1497	do 390 k=1,nl
	1498	do 385 i=1,ncum
	1499	qent(i,k,j)=rr(i,j)
	1500	uent(i,k,j)=u(i,j)
	1501	vent(i,k,j)=v(i,j)
	1502	elij(i,k,j)=0.0
	1503	cym ment(i,k,j)=0.0
	1504	cym sij(i,k,j)=0.0
	1505	385 continue
	1506	390 continue
	1507	400 continue
	1508
	1509	cym
	1510	ment(1:ncum,1:nd,1:nd)=0.0
	1511	sij(1:ncum,1:nd,1:nd)=0.0
	1512
	1513	c do k=1,ntra
	1514	c do j=1,nd ! instead nlp
	1515	c do i=1,nd ! instead nlp
	1516	c do il=1,ncum
	1517	c traent(il,i,j,k)=tra(il,j,k)
	1518	c enddo
	1519	c enddo
	1520	c enddo
	1521	c enddo
	1522	zm(:,:)=0.
	1523
	1524	c=====================================================================
	1525	c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
	1526	c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
	1527	c --- FRACTION (sij)
	1528	c=====================================================================
	1529
	1530	do 750 i=minorig+1, nl
	1531
	1532	do 710 j=minorig,nl
	1533	do 700 il=1,ncum
	1534	if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
	1535	: (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
	1536
	1537	rti=rr(il,1)-ep(il,i)*clw(il,i)
	1538	bf2=1.+lv(il,j)lv(il,j)rs(il,j)/(rrvt(il,j)t(il,j)*cpd)
	1539	anum=h(il,j)-hp(il,i)+(cpv-cpd)t(il,j)(rti-rr(il,j))
	1540	denom=h(il,i)-hp(il,i)+(cpd-cpv)(rr(il,i)-rti)t(il,j)
	1541	dei=denom
	1542	if(abs(dei).lt.0.01)dei=0.01
	1543	sij(il,i,j)=anum/dei
	1544	sij(il,i,i)=1.0
	1545	altem=sij(il,i,j)rr(il,i)+(1.-sij(il,i,j))rti-rs(il,j)
	1546	altem=altem/bf2
	1547	cwat=clw(il,j)*(1.-ep(il,j))
	1548	stemp=sij(il,i,j)
	1549	if((stemp.lt.0.0.or.stemp.gt.1.0.or.altem.gt.cwat)
	1550	: .and.j.gt.i)then
	1551	anum=anum-lv(il,j)(rti-rs(il,j)-cwatbf2)
	1552	denom=denom+lv(il,j)*(rr(il,i)-rti)
	1553	if(abs(denom).lt.0.01)denom=0.01
	1554	sij(il,i,j)=anum/denom
	1555	altem=sij(il,i,j)rr(il,i)+(1.-sij(il,i,j))rti-rs(il,j)
	1556	altem=altem-(bf2-1.)*cwat
	1557	end if
	1558	if(sij(il,i,j).gt.0.0.and.sij(il,i,j).lt.0.95)then
	1559	qent(il,i,j)=sij(il,i,j)rr(il,i)+(1.-sij(il,i,j))rti
	1560	uent(il,i,j)=sij(il,i,j)u(il,i)+(1.-sij(il,i,j))u(il,nk(il))
	1561	vent(il,i,j)=sij(il,i,j)v(il,i)+(1.-sij(il,i,j))v(il,nk(il))
	1562	c!!! do k=1,ntra
	1563	c!!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
	1564	c!!! : +(1.-sij(il,i,j))*tra(il,nk(il),k)
	1565	c!!! end do
	1566	elij(il,i,j)=altem
	1567	elij(il,i,j)=amax1(0.0,elij(il,i,j))
	1568	ment(il,i,j)=m(il,i)/(1.-sij(il,i,j))
	1569	nent(il,i)=nent(il,i)+1
	1570	end if
	1571	sij(il,i,j)=amax1(0.0,sij(il,i,j))
	1572	sij(il,i,j)=amin1(1.0,sij(il,i,j))
	1573	endif ! new
	1574	700 continue
	1575	710 continue
	1576
	1577	c do k=1,ntra
	1578	c do j=minorig,nl
	1579	c do il=1,ncum
	1580	c if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
	1581	c : (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
	1582	c traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
	1583	c : +(1.-sij(il,i,j))*tra(il,nk(il),k)
	1584	c endif
	1585	c enddo
	1586	c enddo
	1587	c enddo
	1588
	1589	c
	1590	c * if no air can entrain at level i assume that updraft detrains *
	1591	c * at that level and calculate detrained air flux and properties *
	1592	c
	1593
	1594	c@ do 170 i=icb(il),inb(il)
	1595
	1596	do 740 il=1,ncum
	1597	if ((i.ge.icb(il)).and.(i.le.inb(il)).and.(nent(il,i).eq.0)) then
	1598	c@ if(nent(il,i).eq.0)then
	1599	ment(il,i,i)=m(il,i)
	1600	qent(il,i,i)=rr(il,nk(il))-ep(il,i)*clw(il,i)
	1601	uent(il,i,i)=u(il,nk(il))
	1602	vent(il,i,i)=v(il,nk(il))
	1603	elij(il,i,i)=clw(il,i)
	1604	cMAF sij(il,i,i)=1.0
	1605	sij(il,i,i)=0.0
	1606	end if
	1607	740 continue
	1608	750 continue
	1609
	1610	c do j=1,ntra
	1611	c do i=minorig+1,nl
	1612	c do il=1,ncum
	1613	c if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
	1614	c traent(il,i,i,j)=tra(il,nk(il),j)
	1615	c endif
	1616	c enddo
	1617	c enddo
	1618	c enddo
	1619
	1620	do 100 j=minorig,nl
	1621	do 101 i=minorig,nl
	1622	do 102 il=1,ncum
	1623	if ((j.ge.(icb(il)-1)).and.(j.le.inb(il))
	1624	: .and.(i.ge.icb(il)).and.(i.le.inb(il)))then
	1625	sigij(il,i,j)=sij(il,i,j)
	1626	endif
	1627	102 continue
	1628	101 continue
	1629	100 continue
	1630	c@ enddo
	1631
	1632	c@170 continue
	1633
	1634	c=====================================================================
	1635	c --- NORMALIZE ENTRAINED AIR MASS FLUXES
	1636	c --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
	1637	c=====================================================================
	1638
	1639	cym call zilch(asum,ncum*nd)
	1640	cym call zilch(bsum,ncum*nd)
	1641	cym call zilch(csum,ncum*nd)
	1642	call zilch(asum,nloc*nd)
	1643	call zilch(csum,nloc*nd)
	1644	call zilch(csum,nloc*nd)
	1645
	1646	do il=1,ncum
	1647	lwork(il) = .FALSE.
	1648	enddo
	1649
	1650	DO 789 i=minorig+1,nl
	1651
	1652	num1=0
	1653	do il=1,ncum
	1654	if ( i.ge.icb(il) .and. i.le.inb(il) ) num1=num1+1
	1655	enddo
	1656	if (num1.le.0) goto 789
	1657
	1658
	1659	do 781 il=1,ncum
	1660	if ( i.ge.icb(il) .and. i.le.inb(il) ) then
	1661	lwork(il)=(nent(il,i).ne.0)
	1662	qp=rr(il,1)-ep(il,i)*clw(il,i)
	1663	anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i))
	1664	: +(cpv-cpd)t(il,i)(qp-rr(il,i))
	1665	denom=h(il,i)-hp(il,i)+lv(il,i)*(rr(il,i)-qp)
	1666	: +(cpd-cpv)t(il,i)(rr(il,i)-qp)
	1667	if(abs(denom).lt.0.01)denom=0.01
	1668	scrit(il)=anum/denom
	1669	alt=qp-rs(il,i)+scrit(il)*(rr(il,i)-qp)
	1670	if(scrit(il).le.0.0.or.alt.le.0.0)scrit(il)=1.0
	1671	smax(il)=0.0
	1672	asij(il)=0.0
	1673	endif
	1674	781 continue
	1675
	1676	do 175 j=nl,minorig,-1
	1677
	1678	num2=0
	1679	do il=1,ncum
	1680	if ( i.ge.icb(il) .and. i.le.inb(il) .and.
	1681	: j.ge.(icb(il)-1) .and. j.le.inb(il)
	1682	: .and. lwork(il) ) num2=num2+1
	1683	enddo
	1684	if (num2.le.0) goto 175
	1685
	1686	do 782 il=1,ncum
	1687	if ( i.ge.icb(il) .and. i.le.inb(il) .and.
	1688	: j.ge.(icb(il)-1) .and. j.le.inb(il)
	1689	: .and. lwork(il) ) then
	1690
	1691	if(sij(il,i,j).gt.1.0e-16.and.sij(il,i,j).lt.0.95)then
	1692	wgh=1.0
	1693	if(j.gt.i)then
	1694	sjmax=amax1(sij(il,i,j+1),smax(il))
	1695	sjmax=amin1(sjmax,scrit(il))
	1696	smax(il)=amax1(sij(il,i,j),smax(il))
	1697	sjmin=amax1(sij(il,i,j-1),smax(il))
	1698	sjmin=amin1(sjmin,scrit(il))
	1699	if(sij(il,i,j).lt.(smax(il)-1.0e-16))wgh=0.0
	1700	smid=amin1(sij(il,i,j),scrit(il))
	1701	else
	1702	sjmax=amax1(sij(il,i,j+1),scrit(il))
	1703	smid=amax1(sij(il,i,j),scrit(il))
	1704	sjmin=0.0
	1705	if(j.gt.1)sjmin=sij(il,i,j-1)
	1706	sjmin=amax1(sjmin,scrit(il))
	1707	endif
	1708	delp=abs(sjmax-smid)
	1709	delm=abs(sjmin-smid)
	1710	asij(il)=asij(il)+wgh*(delp+delm)
	1711	ment(il,i,j)=ment(il,i,j)(delp+delm)wgh
	1712	endif
	1713	endif
	1714	782 continue
	1715
	1716	175 continue
	1717
	1718	do il=1,ncum
	1719	if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
	1720	asij(il)=amax1(1.0e-16,asij(il))
	1721	asij(il)=1.0/asij(il)
	1722	asum(il,i)=0.0
	1723	bsum(il,i)=0.0
	1724	csum(il,i)=0.0
	1725	endif
	1726	enddo
	1727
	1728	do 180 j=minorig,nl
	1729	do il=1,ncum
	1730	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
	1731	: .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
	1732	ment(il,i,j)=ment(il,i,j)*asij(il)
	1733	endif
	1734	enddo
	1735	180 continue
	1736
	1737	do 190 j=minorig,nl
	1738	do il=1,ncum
	1739	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
	1740	: .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
	1741	asum(il,i)=asum(il,i)+ment(il,i,j)
	1742	ment(il,i,j)=ment(il,i,j)*sig(il,j)
	1743	bsum(il,i)=bsum(il,i)+ment(il,i,j)
	1744	endif
	1745	enddo
	1746	190 continue
	1747
	1748	do il=1,ncum
	1749	if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
	1750	bsum(il,i)=amax1(bsum(il,i),1.0e-16)
	1751	bsum(il,i)=1.0/bsum(il,i)
	1752	endif
	1753	enddo
	1754
	1755	do 195 j=minorig,nl
	1756	do il=1,ncum
	1757	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
	1758	: .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
	1759	ment(il,i,j)=ment(il,i,j)asum(il,i)bsum(il,i)
	1760	endif
	1761	enddo
	1762	195 continue
	1763
	1764	do 197 j=minorig,nl
	1765	do il=1,ncum
	1766	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
	1767	: .and. j.ge.(icb(il)-1) .and. j.le.inb(il) ) then
	1768	csum(il,i)=csum(il,i)+ment(il,i,j)
	1769	endif
	1770	enddo
	1771	197 continue
	1772
	1773	do il=1,ncum
	1774	if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
	1775	: .and. csum(il,i).lt.m(il,i) ) then
	1776	nent(il,i)=0
	1777	ment(il,i,i)=m(il,i)
	1778	qent(il,i,i)=rr(il,1)-ep(il,i)*clw(il,i)
	1779	uent(il,i,i)=u(il,nk(il))
	1780	vent(il,i,i)=v(il,nk(il))
	1781	elij(il,i,i)=clw(il,i)
	1782	cMAF sij(il,i,i)=1.0
	1783	sij(il,i,i)=0.0
	1784	endif
	1785	enddo ! il
	1786
	1787	c do j=1,ntra
	1788	c do il=1,ncum
	1789	c if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
	1790	c : .and. csum(il,i).lt.m(il,i) ) then
	1791	c traent(il,i,i,j)=tra(il,nk(il),j)
	1792	c endif
	1793	c enddo
	1794	c enddo
	1795	789 continue
	1796	c
	1797	c MAF: renormalisation de MENT
	1798	do jm=1,nd
	1799	do im=1,nd
	1800	do il=1,ncum
	1801	zm(il,im)=zm(il,im)+(1.-sij(il,im,jm))*ment(il,im,jm)
	1802	end do
	1803	end do
	1804	end do
	1805	c
	1806	do jm=1,nd
	1807	do im=1,nd
	1808	do il=1,ncum
	1809	if(zm(il,im).ne.0.) then
	1810	ment(il,im,jm)=ment(il,im,jm)*m(il,im)/zm(il,im)
	1811	endif
	1812	end do
	1813	end do
	1814	end do
	1815	c
	1816	do jm=1,nd
	1817	do im=1,nd
	1818	do 999 il=1,ncum
	1819	qents(il,im,jm)=qent(il,im,jm)
	1820	ments(il,im,jm)=ment(il,im,jm)
	1821	999 continue
	1822	enddo
	1823	enddo
	1824
	1825	return
	1826	end
	1827
	1828
	1829	SUBROUTINE cv30_unsat(nloc,ncum,nd,na,ntra,icb,inb
	1830	: ,t,rr,rs,gz,u,v,tra,p,ph
	1831	: ,th,tv,lv,cpn,ep,sigp,clw
	1832	: ,m,ment,elij,delt,plcl
	1833	: ,mp,rp,up,vp,trap,wt,water,evap,b)
	1834	implicit none
	1835
	1836
	1837	#include "cvthermo.h"
	1838	#include "cv30param.h"
	1839	#include "cvflag.h"
	1840
	1841	c inputs:
	1842	integer ncum, nd, na, ntra, nloc
	1843	integer icb(nloc), inb(nloc)
	1844	real delt, plcl(nloc)
	1845	real t(nloc,nd), rr(nloc,nd), rs(nloc,nd)
	1846	real u(nloc,nd), v(nloc,nd)
	1847	real tra(nloc,nd,ntra)
	1848	real p(nloc,nd), ph(nloc,nd+1)
	1849	real th(nloc,na), gz(nloc,na)
	1850	real lv(nloc,na), ep(nloc,na), sigp(nloc,na), clw(nloc,na)
	1851	real cpn(nloc,na), tv(nloc,na)
	1852	real m(nloc,na), ment(nloc,na,na), elij(nloc,na,na)
	1853
	1854	c outputs:
	1855	real mp(nloc,na), rp(nloc,na), up(nloc,na), vp(nloc,na)
	1856	real water(nloc,na), evap(nloc,na), wt(nloc,na)
	1857	real trap(nloc,na,ntra)
	1858	real b(nloc,na)
	1859
	1860	c local variables
	1861	integer i,j,k,il,num1
	1862	real tinv, delti
	1863	real awat, afac, afac1, afac2, bfac
	1864	real pr1, pr2, sigt, b6, c6, revap, tevap, delth
	1865	real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
	1866	real ampmax
	1867	real lvcp(nloc,na)
	1868	real wdtrain(nloc)
	1869	logical lwork(nloc)
	1870
	1871
	1872	c------------------------------------------------------
	1873
	1874	delti = 1./delt
	1875	tinv=1./3.
	1876
	1877	mp(:,:)=0.
	1878
	1879	do i=1,nl
	1880	do il=1,ncum
	1881	mp(il,i)=0.0
	1882	rp(il,i)=rr(il,i)
	1883	up(il,i)=u(il,i)
	1884	vp(il,i)=v(il,i)
	1885	wt(il,i)=0.001
	1886	water(il,i)=0.0
	1887	evap(il,i)=0.0
	1888	b(il,i)=0.0
	1889	lvcp(il,i)=lv(il,i)/cpn(il,i)
	1890	enddo
	1891	enddo
	1892
	1893	c do k=1,ntra
	1894	c do i=1,nd
	1895	c do il=1,ncum
	1896	c trap(il,i,k)=tra(il,i,k)
	1897	c enddo
	1898	c enddo
	1899	c enddo
	1900
	1901	c
	1902	c * check whether ep(inb)=0, if so, skip precipitating *
	1903	c * downdraft calculation *
	1904	c
	1905
	1906	do il=1,ncum
	1907	lwork(il)=.TRUE.
	1908	if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
	1909	enddo
	1910
	1911	call zilch(wdtrain,ncum)
	1912
	1913	DO 400 i=nl+1,1,-1
	1914
	1915	num1=0
	1916	do il=1,ncum
	1917	if ( i.le.inb(il) .and. lwork(il) ) num1=num1+1
	1918	enddo
	1919	if (num1.le.0) goto 400
	1920
	1921	c
	1922	c * integrate liquid water equation to find condensed water *
	1923	c * and condensed water flux *
	1924	c
	1925
	1926	c
	1927	c * begin downdraft loop *
	1928	c
	1929
	1930	c
	1931	c * calculate detrained precipitation *
	1932	c
	1933	do il=1,ncum
	1934	if (i.le.inb(il) .and. lwork(il)) then
	1935	if (cvflag_grav) then
	1936	wdtrain(il)=gravep(il,i)m(il,i)*clw(il,i)
	1937	else
	1938	wdtrain(il)=10.0ep(il,i)m(il,i)*clw(il,i)
	1939	endif
	1940	endif
	1941	enddo
	1942
	1943	if(i.gt.1)then
	1944	do 320 j=1,i-1
	1945	do il=1,ncum
	1946	if (i.le.inb(il) .and. lwork(il)) then
	1947	awat=elij(il,j,i)-(1.-ep(il,i))*clw(il,i)
	1948	awat=amax1(awat,0.0)
	1949	if (cvflag_grav) then
	1950	wdtrain(il)=wdtrain(il)+gravawatment(il,j,i)
	1951	else
	1952	wdtrain(il)=wdtrain(il)+10.0awatment(il,j,i)
	1953	endif
	1954	endif
	1955	enddo
	1956	320 continue
	1957	endif
	1958
	1959	c
	1960	c * find rain water and evaporation using provisional *
	1961	c * estimates of rp(i)and rp(i-1) *
	1962	c
	1963
	1964	do 999 il=1,ncum
	1965
	1966	if (i.le.inb(il) .and. lwork(il)) then
	1967
	1968	wt(il,i)=45.0
	1969
	1970	if(i.lt.inb(il))then
	1971	rp(il,i)=rp(il,i+1)
	1972	: +(cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il,i)
	1973	rp(il,i)=0.5*(rp(il,i)+rr(il,i))
	1974	endif
	1975	rp(il,i)=amax1(rp(il,i),0.0)
	1976	rp(il,i)=amin1(rp(il,i),rs(il,i))
	1977	rp(il,inb(il))=rr(il,inb(il))
	1978
	1979	if(i.eq.1)then
	1980	afac=p(il,1)(rs(il,1)-rp(il,1))/(1.0e4+2000.0p(il,1)*rs(il,1))
	1981	else
	1982	rp(il,i-1)=rp(il,i)
	1983	: +(cpd*(t(il,i)-t(il,i-1))+gz(il,i)-gz(il,i-1))/lv(il,i)
	1984	rp(il,i-1)=0.5*(rp(il,i-1)+rr(il,i-1))
	1985	rp(il,i-1)=amin1(rp(il,i-1),rs(il,i-1))
	1986	rp(il,i-1)=amax1(rp(il,i-1),0.0)
	1987	afac1=p(il,i)(rs(il,i)-rp(il,i))/(1.0e4+2000.0p(il,i)*rs(il,i))
	1988	afac2=p(il,i-1)*(rs(il,i-1)-rp(il,i-1))
	1989	: /(1.0e4+2000.0p(il,i-1)rs(il,i-1))
	1990	afac=0.5*(afac1+afac2)
	1991	endif
	1992	if(i.eq.inb(il))afac=0.0
	1993	afac=amax1(afac,0.0)
	1994	bfac=1./(sigd*wt(il,i))
	1995	c
	1996	cjyg1
	1997	ccc sigt=1.0
	1998	ccc if(i.ge.icb)sigt=sigp(i)
	1999	c prise en compte de la variation progressive de sigt dans
	2000	c les couches icb et icb-1:
	2001	c pour plcl<ph(i+1), pr1=0 & pr2=1
	2002	c pour plcl>ph(i), pr1=1 & pr2=0
	2003	c pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
	2004	c sur le nuage, et pr2 est la proportion sous la base du
	2005	c nuage.
	2006	pr1=(plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
	2007	pr1=max(0.,min(1.,pr1))
	2008	pr2=(ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
	2009	pr2=max(0.,min(1.,pr2))
	2010	sigt=sigp(il,i)*pr1+pr2
	2011	cjyg2
	2012	c
	2013	b6=bfac50.sigd(ph(il,i)-ph(il,i+1))sigt*afac
	2014	c6=water(il,i+1)+bfac*wdtrain(il)
	2015	: -50.sigdbfac(ph(il,i)-ph(il,i+1))evap(il,i+1)
	2016	if(c6.gt.0.0)then
	2017	revap=0.5(-b6+sqrt(b6b6+4.*c6))
	2018	evap(il,i)=sigtafacrevap
	2019	water(il,i)=revap*revap
	2020	else
	2021	evap(il,i)=-evap(il,i+1)
	2022	: +0.02(wdtrain(il)+sigdwt(il,i)*water(il,i+1))
	2023	: /(sigd*(ph(il,i)-ph(il,i+1)))
	2024	end if
	2025	c
	2026	c * calculate precipitating downdraft mass flux under *
	2027	c * hydrostatic approximation *
	2028	c
	2029	if (i.ne.1) then
	2030
	2031	tevap=amax1(0.0,evap(il,i))
	2032	delth=amax1(0.001,(th(il,i)-th(il,i-1)))
	2033	if (cvflag_grav) then
	2034	mp(il,i)=100.ginvlvcp(il,i)sigdtevap
	2035	: *(p(il,i-1)-p(il,i))/delth
	2036	else
	2037	mp(il,i)=10.lvcp(il,i)sigdtevap(p(il,i-1)-p(il,i))/delth
	2038	endif
	2039	c
	2040	c * if hydrostatic assumption fails, *
	2041	c * solve cubic difference equation for downdraft theta *
	2042	c * and mass flux from two simultaneous differential eqns *
	2043	c
	2044	amfac=sigdsigd70.0ph(il,i)(p(il,i-1)-p(il,i))
	2045	: (th(il,i)-th(il,i-1))/(tv(il,i)th(il,i))
	2046	amp2=abs(mp(il,i+1)mp(il,i+1)-mp(il,i)mp(il,i))
	2047	if(amp2.gt.(0.1*amfac))then
	2048	xf=100.0sigdsigdsigd(ph(il,i)-ph(il,i+1))
	2049	tf=b(il,i)-5.0(th(il,i)-th(il,i-1))t(il,i)
	2050	: /(lvcp(il,i)sigdth(il,i))
	2051	af=xftf+mp(il,i+1)mp(il,i+1)*tinv
	2052	bf=2.(tinvmp(il,i+1))*3+tinvmp(il,i+1)xftf
	2053	: +50.(p(il,i-1)-p(il,i))xf*tevap
	2054	fac2=1.0
	2055	if(bf.lt.0.0)fac2=-1.0
	2056	bf=abs(bf)
	2057	ur=0.25bfbf-afafaftinvtinv*tinv
	2058	if(ur.ge.0.0)then
	2059	sru=sqrt(ur)
	2060	fac=1.0
	2061	if((0.5*bf-sru).lt.0.0)fac=-1.0
	2062	mp(il,i)=mp(il,i+1)tinv+(0.5bf+sru)**tinv
	2063	: +fac(abs(0.5bf-sru))**tinv
	2064	else
	2065	d=atan(2.*sqrt(-ur)/(bf+1.0e-28))
	2066	if(fac2.lt.0.0)d=3.14159-d
	2067	mp(il,i)=mp(il,i+1)tinv+2.sqrt(aftinv)cos(d*tinv)
	2068	endif
	2069	mp(il,i)=amax1(0.0,mp(il,i))
	2070
	2071	if (cvflag_grav) then
	2072	Cjyg : il y a vraisemblablement une erreur dans la ligne 2 suivante:
	2073	C il faut diviser par (mp(il,i)sigdgrav) et non par (mp(il,i)+sigd*0.1).
	2074	C Et il faut bien revoir les facteurs 100.
	2075	b(il,i-1)=b(il,i)+100.0(p(il,i-1)-p(il,i))tevap
	2076	2 /(mp(il,i)+sigd*0.1)
	2077	3 -10.0(th(il,i)-th(il,i-1))t(il,i)/(lvcp(il,i)sigdth(il,i))
	2078	else
	2079	b(il,i-1)=b(il,i)+100.0(p(il,i-1)-p(il,i))tevap
	2080	2 /(mp(il,i)+sigd*0.1)
	2081	3 -10.0(th(il,i)-th(il,i-1))t(il,i)/(lvcp(il,i)sigdth(il,i))
	2082	endif
	2083	b(il,i-1)=amax1(b(il,i-1),0.0)
	2084	endif
	2085	c
	2086	c * limit magnitude of mp(i) to meet cfl condition *
	2087	c
	2088	ampmax=2.0(ph(il,i)-ph(il,i+1))delti
	2089	amp2=2.0(ph(il,i-1)-ph(il,i))delti
	2090	ampmax=amin1(ampmax,amp2)
	2091	mp(il,i)=amin1(mp(il,i),ampmax)
	2092	c
	2093	c * force mp to decrease linearly to zero *
	2094	c * between cloud base and the surface *
	2095	c
	2096	if(p(il,i).gt.p(il,icb(il)))then
	2097	mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
	2098	endif
	2099
	2100	360 continue
	2101	endif ! i.eq.1
	2102	c
	2103	c * find mixing ratio of precipitating downdraft *
	2104	c
	2105
	2106	if (i.ne.inb(il)) then
	2107
	2108	rp(il,i)=rr(il,i)
	2109
	2110	if(mp(il,i).gt.mp(il,i+1))then
	2111
	2112	if (cvflag_grav) then
	2113	rp(il,i)=rp(il,i+1)mp(il,i+1)+rr(il,i)(mp(il,i)-mp(il,i+1))
	2114	: +100.ginv0.5sigd(ph(il,i)-ph(il,i+1))
	2115	: *(evap(il,i+1)+evap(il,i))
	2116	else
	2117	rp(il,i)=rp(il,i+1)mp(il,i+1)+rr(il,i)(mp(il,i)-mp(il,i+1))
	2118	: +5.sigd(ph(il,i)-ph(il,i+1))
	2119	: *(evap(il,i+1)+evap(il,i))
	2120	endif
	2121	rp(il,i)=rp(il,i)/mp(il,i)
	2122	up(il,i)=up(il,i+1)mp(il,i+1)+u(il,i)(mp(il,i)-mp(il,i+1))
	2123	up(il,i)=up(il,i)/mp(il,i)
	2124	vp(il,i)=vp(il,i+1)mp(il,i+1)+v(il,i)(mp(il,i)-mp(il,i+1))
	2125	vp(il,i)=vp(il,i)/mp(il,i)
	2126
	2127	c do j=1,ntra
	2128	c trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
	2129	ctestmaf : +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
	2130	c : +tra(il,i,j)*(mp(il,i)-mp(il,i+1))
	2131	c trap(il,i,j)=trap(il,i,j)/mp(il,i)
	2132	c end do
	2133
	2134	else
	2135
	2136	if(mp(il,i+1).gt.1.0e-16)then
	2137	if (cvflag_grav) then
	2138	rp(il,i)=rp(il,i+1)
	2139	: +100.ginv0.5sigd(ph(il,i)-ph(il,i+1))
	2140	: *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
	2141	else
	2142	rp(il,i)=rp(il,i+1)
	2143	: +5.sigd(ph(il,i)-ph(il,i+1))
	2144	: *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
	2145	endif
	2146	up(il,i)=up(il,i+1)
	2147	vp(il,i)=vp(il,i+1)
	2148
	2149	c do j=1,ntra
	2150	c trap(il,i,j)=trap(il,i+1,j)
	2151	c end do
	2152
	2153	endif
	2154	endif
	2155	rp(il,i)=amin1(rp(il,i),rs(il,i))
	2156	rp(il,i)=amax1(rp(il,i),0.0)
	2157
	2158	endif
	2159	endif
	2160	999 continue
	2161
	2162	400 continue
	2163
	2164	return
	2165	end
	2166
	2167	SUBROUTINE cv30_yield(nloc,ncum,nd,na,ntra
	2168	: ,icb,inb,delt
	2169	: ,t,rr,u,v,tra,gz,p,ph,h,hp,lv,cpn,th
	2170	: ,ep,clw,m,tp,mp,rp,up,vp,trap
	2171	: ,wt,water,evap,b
	2172	: ,ment,qent,uent,vent,nent,elij,traent,sig
	2173	: ,tv,tvp
	2174	: ,iflag,precip,VPrecip,ft,fr,fu,fv,ftra
	2175	: ,upwd,dnwd,dnwd0,ma,mike,tls,tps,qcondc,wd)
	2176	implicit none
	2177
	2178	#include "cvthermo.h"
	2179	#include "cv30param.h"
	2180	#include "cvflag.h"
	2181	#include "conema3.h"
	2182
	2183	c inputs:
	2184	integer ncum,nd,na,ntra,nloc
	2185	integer icb(nloc), inb(nloc)
	2186	real delt
	2187	real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd)
	2188	real tra(nloc,nd,ntra), sig(nloc,nd)
	2189	real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na)
	2190	real th(nloc,na), p(nloc,nd), tp(nloc,na)
	2191	real lv(nloc,na), cpn(nloc,na), ep(nloc,na), clw(nloc,na)
	2192	real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na)
	2193	real vp(nloc,na), wt(nloc,nd), trap(nloc,nd,ntra)
	2194	real water(nloc,na), evap(nloc,na), b(nloc,na)
	2195	real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na)
	2196	cym real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
	2197	real vent(nloc,na,na), elij(nloc,na,na)
	2198	integer nent(nloc,na)
	2199	real traent(nloc,na,na,ntra)
	2200	real tv(nloc,nd), tvp(nloc,nd)
	2201
	2202	c input/output:
	2203	integer iflag(nloc)
	2204
	2205	c outputs:
	2206	real precip(nloc)
	2207	real VPrecip(nloc,nd+1)
	2208	real ft(nloc,nd), fr(nloc,nd), fu(nloc,nd), fv(nloc,nd)
	2209	real ftra(nloc,nd,ntra)
	2210	real upwd(nloc,nd), dnwd(nloc,nd), ma(nloc,nd)
	2211	real dnwd0(nloc,nd), mike(nloc,nd)
	2212	real tls(nloc,nd), tps(nloc,nd)
	2213	real qcondc(nloc,nd) ! cld
	2214	real wd(nloc) ! gust
	2215
	2216	c local variables:
	2217	integer i,k,il,n,j,num1
	2218	real rat, awat, delti
	2219	real ax, bx, cx, dx, ex
	2220	real cpinv, rdcp, dpinv
	2221	real lvcp(nloc,na), mke(nloc,na)
	2222	real am(nloc), work(nloc), ad(nloc), amp1(nloc)
	2223	c!! real up1(nloc), dn1(nloc)
	2224	real up1(nloc,nd,nd), dn1(nloc,nd,nd)
	2225	real asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
	2226	real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd) ! cld
	2227	real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd) ! cld
	2228
	2229
	2230	c-------------------------------------------------------------
	2231
	2232	c initialization:
	2233
	2234	delti = 1.0/delt
	2235
	2236	do il=1,ncum
	2237	precip(il)=0.0
	2238	wd(il)=0.0 ! gust
	2239	VPrecip(il,nd+1)=0.
	2240	enddo
	2241
	2242	do i=1,nd
	2243	do il=1,ncum
	2244	VPrecip(il,i)=0.0
	2245	ft(il,i)=0.0
	2246	fr(il,i)=0.0
	2247	fu(il,i)=0.0
	2248	fv(il,i)=0.0
	2249	qcondc(il,i)=0.0 ! cld
	2250	qcond(il,i)=0.0 ! cld
	2251	nqcond(il,i)=0.0 ! cld
	2252	enddo
	2253	enddo
	2254
	2255	c do j=1,ntra
	2256	c do i=1,nd
	2257	c do il=1,ncum
	2258	c ftra(il,i,j)=0.0
	2259	c enddo
	2260	c enddo
	2261	c enddo
	2262
	2263	do i=1,nl
	2264	do il=1,ncum
	2265	lvcp(il,i)=lv(il,i)/cpn(il,i)
	2266	enddo
	2267	enddo
	2268
	2269
	2270	c
	2271	c * calculate surface precipitation in mm/day *
	2272	c
	2273	do il=1,ncum
	2274	if(ep(il,inb(il)).ge.0.0001)then
	2275	if (cvflag_grav) then
	2276	precip(il)=wt(il,1)sigdwater(il,1)86400.1000./(rowl*grav)
	2277	else
	2278	precip(il)=wt(il,1)sigdwater(il,1)*8640.
	2279	endif
	2280	endif
	2281	enddo
	2282
	2283	C *** CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg/m2/s ===
	2284	C
	2285	c MAF rajout pour lessivage
	2286	do k=1,nl
	2287	do il=1,ncum
	2288	if (k.le.inb(il)) then
	2289	if (cvflag_grav) then
	2290	VPrecip(il,k) = wt(il,k)sigdwater(il,k)/grav
	2291	else
	2292	VPrecip(il,k) = wt(il,k)sigdwater(il,k)/10.
	2293	endif
	2294	endif
	2295	end do
	2296	end do
	2297	C
	2298	c
	2299	c * Calculate downdraft velocity scale *
	2300	c * NE PAS UTILISER POUR L'INSTANT *
	2301	c
	2302	c! do il=1,ncum
	2303	c! wd(il)=betadabs(mp(il,icb(il)))0.01rrdt(il,icb(il))
	2304	c! : /(sigd*p(il,icb(il)))
	2305	c! enddo
	2306
	2307	c
	2308	c * calculate tendencies of lowest level potential temperature *
	2309	c * and mixing ratio *
	2310	c
	2311	do il=1,ncum
	2312	work(il)=1.0/(ph(il,1)-ph(il,2))
	2313	am(il)=0.0
	2314	enddo
	2315
	2316	do k=2,nl
	2317	do il=1,ncum
	2318	if (k.le.inb(il)) then
	2319	am(il)=am(il)+m(il,k)
	2320	endif
	2321	enddo
	2322	enddo
	2323
	2324	do il=1,ncum
	2325
	2326	c convect3 if((0.1dpinvam).ge.delti)iflag(il)=4
	2327	if (cvflag_grav) then
	2328	if((0.01gravwork(il)*am(il)).ge.delti)iflag(il)=1!consist vect
	2329	ft(il,1)=0.01gravwork(il)am(il)(t(il,2)-t(il,1)
	2330	: +(gz(il,2)-gz(il,1))/cpn(il,1))
	2331	else
	2332	if((0.1work(il)am(il)).ge.delti)iflag(il)=1 !consistency vect
	2333	ft(il,1)=0.1work(il)am(il)*(t(il,2)-t(il,1)
	2334	: +(gz(il,2)-gz(il,1))/cpn(il,1))
	2335	endif
	2336
	2337	ft(il,1)=ft(il,1)-0.5lvcp(il,1)sigd*(evap(il,1)+evap(il,2))
	2338
	2339	if (cvflag_grav) then
	2340	ft(il,1)=ft(il,1)-0.009gravsigd*mp(il,2)
	2341	: t(il,1)b(il,1)*work(il)
	2342	else
	2343	ft(il,1)=ft(il,1)-0.09sigdmp(il,2)t(il,1)b(il,1)*work(il)
	2344	endif
	2345
	2346	ft(il,1)=ft(il,1)+0.01sigdwt(il,1)(cl-cpd)water(il,2)*(t(il,2)
	2347	:-t(il,1))*work(il)/cpn(il,1)
	2348
	2349	if (cvflag_grav) then
	2350	Cjyg1 Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
	2351	c (sb: pour l'instant, on ne fait que le chgt concernant grav, pas evap)
	2352	fr(il,1)=0.01gravmp(il,2)(rp(il,2)-rr(il,1))work(il)
	2353	: +sigd0.5(evap(il,1)+evap(il,2))
	2354	c+tard : +sigd*evap(il,1)
	2355
	2356	fr(il,1)=fr(il,1)+0.01gravam(il)(rr(il,2)-rr(il,1))work(il)
	2357
	2358	fu(il,1)=fu(il,1)+0.01gravwork(il)(mp(il,2)(up(il,2)-u(il,1))
	2359	: +am(il)*(u(il,2)-u(il,1)))
	2360	fv(il,1)=fv(il,1)+0.01gravwork(il)(mp(il,2)(vp(il,2)-v(il,1))
	2361	: +am(il)*(v(il,2)-v(il,1)))
	2362	else ! cvflag_grav
	2363	fr(il,1)=0.1mp(il,2)(rp(il,2)-rr(il,1))*work(il)
	2364	: +sigd0.5(evap(il,1)+evap(il,2))
	2365	fr(il,1)=fr(il,1)+0.1am(il)(rr(il,2)-rr(il,1))*work(il)
	2366	fu(il,1)=fu(il,1)+0.1work(il)(mp(il,2)*(up(il,2)-u(il,1))
	2367	: +am(il)*(u(il,2)-u(il,1)))
	2368	fv(il,1)=fv(il,1)+0.1work(il)(mp(il,2)*(vp(il,2)-v(il,1))
	2369	: +am(il)*(v(il,2)-v(il,1)))
	2370	endif ! cvflag_grav
	2371
	2372	enddo ! il
	2373
	2374	c do j=1,ntra
	2375	c do il=1,ncum
	2376	c if (cvflag_grav) then
	2377	c ftra(il,1,j)=ftra(il,1,j)+0.01gravwork(il)
	2378	c : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
	2379	c : +am(il)*(tra(il,2,j)-tra(il,1,j)))
	2380	c else
	2381	c ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
	2382	c : (mp(il,2)(trap(il,2,j)-tra(il,1,j))
	2383	c : +am(il)*(tra(il,2,j)-tra(il,1,j)))
	2384	c endif
	2385	c enddo
	2386	c enddo
	2387
	2388	do j=2,nl
	2389	do il=1,ncum
	2390	if (j.le.inb(il)) then
	2391	if (cvflag_grav) then
	2392	fr(il,1)=fr(il,1)
	2393	: +0.01gravwork(il)ment(il,j,1)(qent(il,j,1)-rr(il,1))
	2394	fu(il,1)=fu(il,1)
	2395	: +0.01gravwork(il)ment(il,j,1)(uent(il,j,1)-u(il,1))
	2396	fv(il,1)=fv(il,1)
	2397	: +0.01gravwork(il)ment(il,j,1)(vent(il,j,1)-v(il,1))
	2398	else ! cvflag_grav
	2399	fr(il,1)=fr(il,1)
	2400	: +0.1work(il)ment(il,j,1)*(qent(il,j,1)-rr(il,1))
	2401	fu(il,1)=fu(il,1)
	2402	: +0.1work(il)ment(il,j,1)*(uent(il,j,1)-u(il,1))
	2403	fv(il,1)=fv(il,1)
	2404	: +0.1work(il)ment(il,j,1)*(vent(il,j,1)-v(il,1))
	2405	endif ! cvflag_grav
	2406	endif ! j
	2407	enddo
	2408	enddo
	2409
	2410	c do k=1,ntra
	2411	c do j=2,nl
	2412	c do il=1,ncum
	2413	c if (j.le.inb(il)) then
	2414
	2415	c if (cvflag_grav) then
	2416	c ftra(il,1,k)=ftra(il,1,k)+0.01gravwork(il)*ment(il,j,1)
	2417	c : *(traent(il,j,1,k)-tra(il,1,k))
	2418	c else
	2419	c ftra(il,1,k)=ftra(il,1,k)+0.1work(il)ment(il,j,1)
	2420	c : *(traent(il,j,1,k)-tra(il,1,k))
	2421	c endif
	2422
	2423	c endif
	2424	c enddo
	2425	c enddo
	2426	c enddo
	2427
	2428	c
	2429	c * calculate tendencies of potential temperature and mixing ratio *
	2430	c * at levels above the lowest level *
	2431	c
	2432	c * first find the net saturated updraft and downdraft mass fluxes *
	2433	c * through each level *
	2434	c
	2435
	2436	do 500 i=2,nl+1 ! newvecto: mettre nl au lieu nl+1?
	2437
	2438	num1=0
	2439	do il=1,ncum
	2440	if(i.le.inb(il))num1=num1+1
	2441	enddo
	2442	if(num1.le.0)go to 500
	2443
	2444	call zilch(amp1,ncum)
	2445	call zilch(ad,ncum)
	2446
	2447	do 440 k=i+1,nl+1
	2448	do 441 il=1,ncum
	2449	if (i.le.inb(il) .and. k.le.(inb(il)+1)) then
	2450	amp1(il)=amp1(il)+m(il,k)
	2451	endif
	2452	441 continue
	2453	440 continue
	2454
	2455	do 450 k=1,i
	2456	do 451 j=i+1,nl+1
	2457	do 452 il=1,ncum
	2458	if (i.le.inb(il) .and. j.le.(inb(il)+1)) then
	2459	amp1(il)=amp1(il)+ment(il,k,j)
	2460	endif
	2461	452 continue
	2462	451 continue
	2463	450 continue
	2464
	2465	do 470 k=1,i-1
	2466	do 471 j=i,nl+1 ! newvecto: nl au lieu nl+1?
	2467	do 472 il=1,ncum
	2468	if (i.le.inb(il) .and. j.le.inb(il)) then
	2469	ad(il)=ad(il)+ment(il,j,k)
	2470	endif
	2471	472 continue
	2472	471 continue
	2473	470 continue
	2474
	2475	do 1350 il=1,ncum
	2476	if (i.le.inb(il)) then
	2477	dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2478	cpinv=1.0/cpn(il,i)
	2479
	2480	c convect3 if((0.1dpinvamp1).ge.delti)iflag(il)=4
	2481	if (cvflag_grav) then
	2482	if((0.01gravdpinv*amp1(il)).ge.delti)iflag(il)=1 ! vecto
	2483	else
	2484	if((0.1dpinvamp1(il)).ge.delti)iflag(il)=1 ! vecto
	2485	endif
	2486
	2487	if (cvflag_grav) then
	2488	ft(il,i)=0.01gravdpinv(amp1(il)(t(il,i+1)-t(il,i)
	2489	: +(gz(il,i+1)-gz(il,i))*cpinv)
	2490	: -ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
	2491	: -0.5sigdlvcp(il,i)*(evap(il,i)+evap(il,i+1))
	2492	rat=cpn(il,i-1)*cpinv
	2493	ft(il,i)=ft(il,i)-0.009gravsigd(mp(il,i+1)t(il,i)*b(il,i)
	2494	: -mp(il,i)t(il,i-1)ratb(il,i-1))dpinv
	2495	ft(il,i)=ft(il,i)+0.01gravdpinvment(il,i,i)(hp(il,i)-h(il,i)
	2496	: +t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
	2497	else ! cvflag_grav
	2498	ft(il,i)=0.1dpinv(amp1(il)*(t(il,i+1)-t(il,i)
	2499	: +(gz(il,i+1)-gz(il,i))*cpinv)
	2500	: -ad(il)(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))cpinv))
	2501	: -0.5sigdlvcp(il,i)*(evap(il,i)+evap(il,i+1))
	2502	rat=cpn(il,i-1)*cpinv
	2503	ft(il,i)=ft(il,i)-0.09sigd(mp(il,i+1)t(il,i)b(il,i)
	2504	: -mp(il,i)t(il,i-1)ratb(il,i-1))dpinv
	2505	ft(il,i)=ft(il,i)+0.1dpinvment(il,i,i)*(hp(il,i)-h(il,i)
	2506	: +t(il,i)(cpv-cpd)(rr(il,i)-qent(il,i,i)))*cpinv
	2507	endif ! cvflag_grav
	2508
	2509
	2510	ft(il,i)=ft(il,i)+0.01sigdwt(il,i)(cl-cpd)water(il,i+1)
	2511	: (t(il,i+1)-t(il,i))dpinv*cpinv
	2512
	2513	if (cvflag_grav) then
	2514	fr(il,i)=0.01gravdpinv(amp1(il)(rr(il,i+1)-rr(il,i))
	2515	: -ad(il)*(rr(il,i)-rr(il,i-1)))
	2516	fu(il,i)=fu(il,i)+0.01gravdpinv(amp1(il)(u(il,i+1)-u(il,i))
	2517	: -ad(il)*(u(il,i)-u(il,i-1)))
	2518	fv(il,i)=fv(il,i)+0.01gravdpinv(amp1(il)(v(il,i+1)-v(il,i))
	2519	: -ad(il)*(v(il,i)-v(il,i-1)))
	2520	else ! cvflag_grav
	2521	fr(il,i)=0.1dpinv(amp1(il)*(rr(il,i+1)-rr(il,i))
	2522	: -ad(il)*(rr(il,i)-rr(il,i-1)))
	2523	fu(il,i)=fu(il,i)+0.1dpinv(amp1(il)*(u(il,i+1)-u(il,i))
	2524	: -ad(il)*(u(il,i)-u(il,i-1)))
	2525	fv(il,i)=fv(il,i)+0.1dpinv(amp1(il)*(v(il,i+1)-v(il,i))
	2526	: -ad(il)*(v(il,i)-v(il,i-1)))
	2527	endif ! cvflag_grav
	2528
	2529	endif ! i
	2530	1350 continue
	2531
	2532	c do k=1,ntra
	2533	c do il=1,ncum
	2534	c if (i.le.inb(il)) then
	2535	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2536	c cpinv=1.0/cpn(il,i)
	2537	c if (cvflag_grav) then
	2538	c ftra(il,i,k)=ftra(il,i,k)+0.01gravdpinv
	2539	c : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
	2540	c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
	2541	c else
	2542	c ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
	2543	c : (amp1(il)(tra(il,i+1,k)-tra(il,i,k))
	2544	c : -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
	2545	c endif
	2546	c endif
	2547	c enddo
	2548	c enddo
	2549
	2550	do 480 k=1,i-1
	2551	do 1370 il=1,ncum
	2552	if (i.le.inb(il)) then
	2553	dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2554	cpinv=1.0/cpn(il,i)
	2555
	2556	awat=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)
	2557	awat=amax1(awat,0.0)
	2558
	2559	if (cvflag_grav) then
	2560	fr(il,i)=fr(il,i)
	2561	: +0.01gravdpinvment(il,k,i)(qent(il,k,i)-awat-rr(il,i))
	2562	fu(il,i)=fu(il,i)
	2563	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
	2564	fv(il,i)=fv(il,i)
	2565	: +0.01gravdpinvment(il,k,i)(vent(il,k,i)-v(il,i))
	2566	else ! cvflag_grav
	2567	fr(il,i)=fr(il,i)
	2568	: +0.1dpinvment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))
	2569	fu(il,i)=fu(il,i)
	2570	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
	2571	fv(il,i)=fv(il,i)
	2572	: +0.1dpinvment(il,k,i)*(vent(il,k,i)-v(il,i))
	2573	endif ! cvflag_grav
	2574
	2575	c (saturated updrafts resulting from mixing) ! cld
	2576	qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat) ! cld
	2577	nqcond(il,i)=nqcond(il,i)+1. ! cld
	2578	endif ! i
	2579	1370 continue
	2580	480 continue
	2581
	2582	c do j=1,ntra
	2583	c do k=1,i-1
	2584	c do il=1,ncum
	2585	c if (i.le.inb(il)) then
	2586	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2587	c cpinv=1.0/cpn(il,i)
	2588	c if (cvflag_grav) then
	2589	c ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
	2590	c : *(traent(il,k,i,j)-tra(il,i,j))
	2591	c else
	2592	c ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
	2593	c : *(traent(il,k,i,j)-tra(il,i,j))
	2594	c endif
	2595	c endif
	2596	c enddo
	2597	c enddo
	2598	c enddo
	2599
	2600	do 490 k=i,nl+1
	2601	do 1380 il=1,ncum
	2602	if (i.le.inb(il) .and. k.le.inb(il)) then
	2603	dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2604	cpinv=1.0/cpn(il,i)
	2605
	2606	if (cvflag_grav) then
	2607	fr(il,i)=fr(il,i)
	2608	: +0.01gravdpinvment(il,k,i)(qent(il,k,i)-rr(il,i))
	2609	fu(il,i)=fu(il,i)
	2610	: +0.01gravdpinvment(il,k,i)(uent(il,k,i)-u(il,i))
	2611	fv(il,i)=fv(il,i)
	2612	: +0.01gravdpinvment(il,k,i)(vent(il,k,i)-v(il,i))
	2613	else ! cvflag_grav
	2614	fr(il,i)=fr(il,i)
	2615	: +0.1dpinvment(il,k,i)*(qent(il,k,i)-rr(il,i))
	2616	fu(il,i)=fu(il,i)
	2617	: +0.1dpinvment(il,k,i)*(uent(il,k,i)-u(il,i))
	2618	fv(il,i)=fv(il,i)
	2619	: +0.1dpinvment(il,k,i)*(vent(il,k,i)-v(il,i))
	2620	endif ! cvflag_grav
	2621	endif ! i and k
	2622	1380 continue
	2623	490 continue
	2624
	2625	c do j=1,ntra
	2626	c do k=i,nl+1
	2627	c do il=1,ncum
	2628	c if (i.le.inb(il) .and. k.le.inb(il)) then
	2629	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2630	c cpinv=1.0/cpn(il,i)
	2631	c if (cvflag_grav) then
	2632	c ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv*ment(il,k,i)
	2633	c : *(traent(il,k,i,j)-tra(il,i,j))
	2634	c else
	2635	c ftra(il,i,j)=ftra(il,i,j)+0.1dpinvment(il,k,i)
	2636	c : *(traent(il,k,i,j)-tra(il,i,j))
	2637	c endif
	2638	c endif ! i and k
	2639	c enddo
	2640	c enddo
	2641	c enddo
	2642
	2643	do 1400 il=1,ncum
	2644	if (i.le.inb(il)) then
	2645	dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2646	cpinv=1.0/cpn(il,i)
	2647
	2648	if (cvflag_grav) then
	2649	c sb: on ne fait pas encore la correction permettant de mieux
	2650	c conserver l'eau:
	2651	fr(il,i)=fr(il,i)+0.5sigd(evap(il,i)+evap(il,i+1))
	2652	: +0.01grav(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i)
	2653	: (rp(il,i)-rr(il,i-1)))dpinv
	2654
	2655	fu(il,i)=fu(il,i)+0.01grav(mp(il,i+1)*(up(il,i+1)-u(il,i))
	2656	: -mp(il,i)(up(il,i)-u(il,i-1)))dpinv
	2657	fv(il,i)=fv(il,i)+0.01grav(mp(il,i+1)*(vp(il,i+1)-v(il,i))
	2658	: -mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
	2659	else ! cvflag_grav
	2660	fr(il,i)=fr(il,i)+0.5sigd(evap(il,i)+evap(il,i+1))
	2661	: +0.1(mp(il,i+1)(rp(il,i+1)-rr(il,i))-mp(il,i)
	2662	: (rp(il,i)-rr(il,i-1)))dpinv
	2663	fu(il,i)=fu(il,i)+0.1(mp(il,i+1)(up(il,i+1)-u(il,i))
	2664	: -mp(il,i)(up(il,i)-u(il,i-1)))dpinv
	2665	fv(il,i)=fv(il,i)+0.1(mp(il,i+1)(vp(il,i+1)-v(il,i))
	2666	: -mp(il,i)(vp(il,i)-v(il,i-1)))dpinv
	2667	endif ! cvflag_grav
	2668
	2669	endif ! i
	2670	1400 continue
	2671
	2672	c sb: interface with the cloud parameterization: ! cld
	2673
	2674	do k=i+1,nl
	2675	do il=1,ncum
	2676	if (k.le.inb(il) .and. i.le.inb(il)) then ! cld
	2677	C (saturated downdrafts resulting from mixing) ! cld
	2678	qcond(il,i)=qcond(il,i)+elij(il,k,i) ! cld
	2679	nqcond(il,i)=nqcond(il,i)+1. ! cld
	2680	endif ! cld
	2681	enddo ! cld
	2682	enddo ! cld
	2683
	2684	C (particular case: no detraining level is found) ! cld
	2685	do il=1,ncum ! cld
	2686	if (i.le.inb(il) .and. nent(il,i).eq.0) then ! cld
	2687	qcond(il,i)=qcond(il,i)+(1.-ep(il,i))*clw(il,i) ! cld
	2688	nqcond(il,i)=nqcond(il,i)+1. ! cld
	2689	endif ! cld
	2690	enddo ! cld
	2691
	2692	do il=1,ncum ! cld
	2693	if (i.le.inb(il) .and. nqcond(il,i).ne.0.) then ! cld
	2694	qcond(il,i)=qcond(il,i)/nqcond(il,i) ! cld
	2695	endif ! cld
	2696	enddo
	2697
	2698	c do j=1,ntra
	2699	c do il=1,ncum
	2700	c if (i.le.inb(il)) then
	2701	c dpinv=1.0/(ph(il,i)-ph(il,i+1))
	2702	c cpinv=1.0/cpn(il,i)
	2703
	2704	c if (cvflag_grav) then
	2705	c ftra(il,i,j)=ftra(il,i,j)+0.01gravdpinv
	2706	c : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
	2707	c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
	2708	c else
	2709	c ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
	2710	c : (mp(il,i+1)(trap(il,i+1,j)-tra(il,i,j))
	2711	c : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
	2712	c endif
	2713	c endif ! i
	2714	c enddo
	2715	c enddo
	2716
	2717	500 continue
	2718
	2719
	2720	c * move the detrainment at level inb down to level inb-1 *
	2721	c * in such a way as to preserve the vertically *
	2722	c * integrated enthalpy and water tendencies *
	2723	c
	2724	do 503 il=1,ncum
	2725
	2726	ax=0.1ment(il,inb(il),inb(il))(hp(il,inb(il))-h(il,inb(il))
	2727	: +t(il,inb(il))*(cpv-cpd)
	2728	: *(rr(il,inb(il))-qent(il,inb(il),inb(il))))
	2729	: /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
	2730	ft(il,inb(il))=ft(il,inb(il))-ax
	2731	ft(il,inb(il)-1)=ft(il,inb(il)-1)+ax*cpn(il,inb(il))
	2732	: *(ph(il,inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)
	2733	: *(ph(il,inb(il)-1)-ph(il,inb(il))))
	2734
	2735	bx=0.1ment(il,inb(il),inb(il))(qent(il,inb(il),inb(il))
	2736	: -rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
	2737	fr(il,inb(il))=fr(il,inb(il))-bx
	2738	fr(il,inb(il)-1)=fr(il,inb(il)-1)
	2739	: +bx*(ph(il,inb(il))-ph(il,inb(il)+1))
	2740	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
	2741
	2742	cx=0.1ment(il,inb(il),inb(il))(uent(il,inb(il),inb(il))
	2743	: -u(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
	2744	fu(il,inb(il))=fu(il,inb(il))-cx
	2745	fu(il,inb(il)-1)=fu(il,inb(il)-1)
	2746	: +cx*(ph(il,inb(il))-ph(il,inb(il)+1))
	2747	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
	2748
	2749	dx=0.1ment(il,inb(il),inb(il))(vent(il,inb(il),inb(il))
	2750	: -v(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
	2751	fv(il,inb(il))=fv(il,inb(il))-dx
	2752	fv(il,inb(il)-1)=fv(il,inb(il)-1)
	2753	: +dx*(ph(il,inb(il))-ph(il,inb(il)+1))
	2754	: /(ph(il,inb(il)-1)-ph(il,inb(il)))
	2755
	2756	503 continue
	2757
	2758	c do j=1,ntra
	2759	c do il=1,ncum
	2760	c ex=0.1*ment(il,inb(il),inb(il))
	2761	c : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
	2762	c : /(ph(il,inb(il))-ph(il,inb(il)+1))
	2763	c ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
	2764	c ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
	2765	c : +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
	2766	c : /(ph(il,inb(il)-1)-ph(il,inb(il)))
	2767	c enddo
	2768	c enddo
	2769
	2770	c
	2771	c * homoginize tendencies below cloud base *
	2772	c
	2773	c
	2774	do il=1,ncum
	2775	asum(il)=0.0
	2776	bsum(il)=0.0
	2777	csum(il)=0.0
	2778	dsum(il)=0.0
	2779	enddo
	2780
	2781	do i=1,nl
	2782	do il=1,ncum
	2783	if (i.le.(icb(il)-1)) then
	2784	asum(il)=asum(il)+ft(il,i)*(ph(il,i)-ph(il,i+1))
	2785	bsum(il)=bsum(il)+fr(il,i)(lv(il,i)+(cl-cpd)(t(il,i)-t(il,1)))
	2786	: *(ph(il,i)-ph(il,i+1))
	2787	csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))
	2788	: *(ph(il,i)-ph(il,i+1))
	2789	dsum(il)=dsum(il)+t(il,i)*(ph(il,i)-ph(il,i+1))/th(il,i)
	2790	endif
	2791	enddo
	2792	enddo
	2793
	2794	c!!! do 700 i=1,icb(il)-1
	2795	do i=1,nl
	2796	do il=1,ncum
	2797	if (i.le.(icb(il)-1)) then
	2798	ft(il,i)=asum(il)t(il,i)/(th(il,i)dsum(il))
	2799	fr(il,i)=bsum(il)/csum(il)
	2800	endif
	2801	enddo
	2802	enddo
	2803
	2804	c
	2805	c * reset counter and return *
	2806	c
	2807	do il=1,ncum
	2808	sig(il,nd)=2.0
	2809	enddo
	2810
	2811
	2812	do i=1,nd
	2813	do il=1,ncum
	2814	upwd(il,i)=0.0
	2815	dnwd(il,i)=0.0
	2816	enddo
	2817	enddo
	2818
	2819	do i=1,nl
	2820	do il=1,ncum
	2821	dnwd0(il,i)=-mp(il,i)
	2822	enddo
	2823	enddo
	2824	do i=nl+1,nd
	2825	do il=1,ncum
	2826	dnwd0(il,i)=0.
	2827	enddo
	2828	enddo
	2829
	2830
	2831	do i=1,nl
	2832	do il=1,ncum
	2833	if (i.ge.icb(il) .and. i.le.inb(il)) then
	2834	upwd(il,i)=0.0
	2835	dnwd(il,i)=0.0
	2836	endif
	2837	enddo
	2838	enddo
	2839
	2840	do i=1,nl
	2841	do k=1,nl
	2842	do il=1,ncum
	2843	up1(il,k,i)=0.0
	2844	dn1(il,k,i)=0.0
	2845	enddo
	2846	enddo
	2847	enddo
	2848
	2849	do i=1,nl
	2850	do k=i,nl
	2851	do n=1,i-1
	2852	do il=1,ncum
	2853	if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
	2854	up1(il,k,i)=up1(il,k,i)+ment(il,n,k)
	2855	dn1(il,k,i)=dn1(il,k,i)-ment(il,k,n)
	2856	endif
	2857	enddo
	2858	enddo
	2859	enddo
	2860	enddo
	2861
	2862	do i=2,nl
	2863	do k=i,nl
	2864	do il=1,ncum
	2865	ctest if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
	2866	if (i.le.inb(il).and.k.le.inb(il)) then
	2867	upwd(il,i)=upwd(il,i)+m(il,k)+up1(il,k,i)
	2868	dnwd(il,i)=dnwd(il,i)+dn1(il,k,i)
	2869	endif
	2870	enddo
	2871	enddo
	2872	enddo
	2873
	2874
	2875	c!!! DO il=1,ncum
	2876	c!!! do i=icb(il),inb(il)
	2877	c!!!
	2878	c!!! upwd(il,i)=0.0
	2879	c!!! dnwd(il,i)=0.0
	2880	c!!! do k=i,inb(il)
	2881	c!!! up1=0.0
	2882	c!!! dn1=0.0
	2883	c!!! do n=1,i-1
	2884	c!!! up1=up1+ment(il,n,k)
	2885	c!!! dn1=dn1-ment(il,k,n)
	2886	c!!! enddo
	2887	c!!! upwd(il,i)=upwd(il,i)+m(il,k)+up1
	2888	c!!! dnwd(il,i)=dnwd(il,i)+dn1
	2889	c!!! enddo
	2890	c!!! enddo
	2891	c!!!
	2892	c!!! ENDDO
	2893
	2894	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2895	c determination de la variation de flux ascendant entre
	2896	c deux niveau non dilue mike
	2897	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2898
	2899	do i=1,nl
	2900	do il=1,ncum
	2901	mike(il,i)=m(il,i)
	2902	enddo
	2903	enddo
	2904
	2905	do i=nl+1,nd
	2906	do il=1,ncum
	2907	mike(il,i)=0.
	2908	enddo
	2909	enddo
	2910
	2911	do i=1,nd
	2912	do il=1,ncum
	2913	ma(il,i)=0
	2914	enddo
	2915	enddo
	2916
	2917	do i=1,nl
	2918	do j=i,nl
	2919	do il=1,ncum
	2920	ma(il,i)=ma(il,i)+m(il,j)
	2921	enddo
	2922	enddo
	2923	enddo
	2924
	2925	do i=nl+1,nd
	2926	do il=1,ncum
	2927	ma(il,i)=0.
	2928	enddo
	2929	enddo
	2930
	2931	do i=1,nl
	2932	do il=1,ncum
	2933	if (i.le.(icb(il)-1)) then
	2934	ma(il,i)=0
	2935	endif
	2936	enddo
	2937	enddo
	2938
	2939	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2940	c icb represente de niveau ou se trouve la
	2941	c base du nuage , et inb le top du nuage
	2942	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2943
	2944	do i=1,nd
	2945	do il=1,ncum
	2946	mke(il,i)=upwd(il,i)+dnwd(il,i)
	2947	enddo
	2948	enddo
	2949
	2950	do i=1,nd
	2951	DO 999 il=1,ncum
	2952	rdcp=(rrd(1.-rr(il,i))-rr(il,i)rrv)
	2953	: /(cpd(1.-rr(il,i))+rr(il,i)cpv)
	2954	tls(il,i)=t(il,i)(1000.0/p(il,i))*rdcp
	2955	tps(il,i)=tp(il,i)
	2956	999 CONTINUE
	2957	enddo
	2958
	2959	c
	2960	c * diagnose the in-cloud mixing ratio * ! cld
	2961	c * of condensed water * ! cld
	2962	c ! cld
	2963
	2964	do i=1,nd ! cld
	2965	do il=1,ncum ! cld
	2966	mac(il,i)=0.0 ! cld
	2967	wa(il,i)=0.0 ! cld
	2968	siga(il,i)=0.0 ! cld
	2969	sax(il,i)=0.0 ! cld
	2970	enddo ! cld
	2971	enddo ! cld
	2972
	2973	do i=minorig, nl ! cld
	2974	do k=i+1,nl+1 ! cld
	2975	do il=1,ncum ! cld
	2976	if (i.le.inb(il) .and. k.le.(inb(il)+1)) then ! cld
	2977	mac(il,i)=mac(il,i)+m(il,k) ! cld
	2978	endif ! cld
	2979	enddo ! cld
	2980	enddo ! cld
	2981	enddo ! cld
	2982
	2983	do i=1,nl ! cld
	2984	do j=1,i ! cld
	2985	do il=1,ncum ! cld
	2986	if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
	2987	: .and. j.ge.icb(il) ) then ! cld
	2988	sax(il,i)=sax(il,i)+rrd*(tvp(il,j)-tv(il,j)) ! cld
	2989	: *(ph(il,j)-ph(il,j+1))/p(il,j) ! cld
	2990	endif ! cld
	2991	enddo ! cld
	2992	enddo ! cld
	2993	enddo ! cld
	2994
	2995	do i=1,nl ! cld
	2996	do il=1,ncum ! cld
	2997	if (i.ge.icb(il) .and. i.le.(inb(il)-1) ! cld
	2998	: .and. sax(il,i).gt.0.0 ) then ! cld
	2999	wa(il,i)=sqrt(2.*sax(il,i)) ! cld
	3000	endif ! cld
	3001	enddo ! cld
	3002	enddo ! cld
	3003
	3004	do i=1,nl ! cld
	3005	do il=1,ncum ! cld
	3006	if (wa(il,i).gt.0.0) ! cld
	3007	: siga(il,i)=mac(il,i)/wa(il,i) ! cld
	3008	: rrdtvp(il,i)/p(il,i)/100./delta ! cld
	3009	siga(il,i) = min(siga(il,i),1.0) ! cld
	3010	cIM cf. FH
	3011	if (iflag_clw.eq.0) then
	3012	qcondc(il,i)=siga(il,i)clw(il,i)(1.-ep(il,i)) ! cld
	3013	: + (1.-siga(il,i))*qcond(il,i) ! cld
	3014	else if (iflag_clw.eq.1) then
	3015	qcondc(il,i)=qcond(il,i) ! cld
	3016	endif
	3017
	3018	enddo ! cld
	3019	enddo ! cld
	3020
	3021	return
	3022	end
	3023
	3024	SUBROUTINE cv30_tracer(nloc,len,ncum,nd,na,
	3025	& ment,sij,da,phi)
	3026	implicit none
	3027	c inputs:
	3028	integer ncum, nd, na, nloc,len
	3029	real ment(nloc,na,na),sij(nloc,na,na)
	3030	c ouputs:
	3031	real da(nloc,na),phi(nloc,na,na)
	3032	c local variables:
	3033	integer i,j,k
	3034	c
	3035	da(:,:)=0.
	3036	c
	3037	do j=1,na
	3038	do k=1,na
	3039	do i=1,ncum
	3040	da(i,j)=da(i,j)+(1.-sij(i,k,j))*ment(i,k,j)
	3041	phi(i,j,k)=sij(i,k,j)*ment(i,k,j)
	3042	c print *,'da',j,k,da(i,j),sij(i,k,j),ment(i,k,j)
	3043	end do
	3044	end do
	3045	end do
	3046
	3047	return
	3048	end
	3049
	3050
	3051	SUBROUTINE cv30_uncompress(nloc,len,ncum,nd,ntra,idcum
	3052	: ,iflag
	3053	: ,precip,VPrecip,sig,w0
	3054	: ,ft,fq,fu,fv,ftra
	3055	: ,inb
	3056	: ,Ma,upwd,dnwd,dnwd0,qcondc,wd,cape
	3057	: ,da,phi,mp
	3058	: ,iflag1
	3059	: ,precip1,VPrecip1,sig1,w01
	3060	: ,ft1,fq1,fu1,fv1,ftra1
	3061	: ,inb1
	3062	: ,Ma1,upwd1,dnwd1,dnwd01,qcondc1,wd1,cape1
	3063	: ,da1,phi1,mp1)
	3064	implicit none
	3065
	3066	#include "cv30param.h"
	3067
	3068	c inputs:
	3069	integer len, ncum, nd, ntra, nloc
	3070	integer idcum(nloc)
	3071	integer iflag(nloc)
	3072	integer inb(nloc)
	3073	real precip(nloc)
	3074	real VPrecip(nloc,nd+1)
	3075	real sig(nloc,nd), w0(nloc,nd)
	3076	real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
	3077	real ftra(nloc,nd,ntra)
	3078	real Ma(nloc,nd)
	3079	real upwd(nloc,nd),dnwd(nloc,nd),dnwd0(nloc,nd)
	3080	real qcondc(nloc,nd)
	3081	real wd(nloc),cape(nloc)
	3082	real da(nloc,nd),phi(nloc,nd,nd),mp(nloc,nd)
	3083
	3084	c outputs:
	3085	integer iflag1(len)
	3086	integer inb1(len)
	3087	real precip1(len)
	3088	real VPrecip1(len,nd+1)
	3089	real sig1(len,nd), w01(len,nd)
	3090	real ft1(len,nd), fq1(len,nd), fu1(len,nd), fv1(len,nd)
	3091	real ftra1(len,nd,ntra)
	3092	real Ma1(len,nd)
	3093	real upwd1(len,nd),dnwd1(len,nd),dnwd01(len,nd)
	3094	real qcondc1(nloc,nd)
	3095	real wd1(nloc),cape1(nloc)
	3096	real da1(nloc,nd),phi1(nloc,nd,nd),mp1(nloc,nd)
	3097
	3098	c local variables:
	3099	integer i,k,j
	3100
	3101	do 2000 i=1,ncum
	3102	precip1(idcum(i))=precip(i)
	3103	iflag1(idcum(i))=iflag(i)
	3104	wd1(idcum(i))=wd(i)
	3105	inb1(idcum(i))=inb(i)
	3106	cape1(idcum(i))=cape(i)
	3107	2000 continue
	3108
	3109	do 2020 k=1,nl
	3110	do 2010 i=1,ncum
	3111	VPrecip1(idcum(i),k)=VPrecip(i,k)
	3112	sig1(idcum(i),k)=sig(i,k)
	3113	w01(idcum(i),k)=w0(i,k)
	3114	ft1(idcum(i),k)=ft(i,k)
	3115	fq1(idcum(i),k)=fq(i,k)
	3116	fu1(idcum(i),k)=fu(i,k)
	3117	fv1(idcum(i),k)=fv(i,k)
	3118	Ma1(idcum(i),k)=Ma(i,k)
	3119	upwd1(idcum(i),k)=upwd(i,k)
	3120	dnwd1(idcum(i),k)=dnwd(i,k)
	3121	dnwd01(idcum(i),k)=dnwd0(i,k)
	3122	qcondc1(idcum(i),k)=qcondc(i,k)
	3123	da1(idcum(i),k)=da(i,k)
	3124	mp1(idcum(i),k)=mp(i,k)
	3125	2010 continue
	3126	2020 continue
	3127
	3128	do 2200 i=1,ncum
	3129	sig1(idcum(i),nd)=sig(i,nd)
	3130	2200 continue
	3131
	3132
	3133	c do 2100 j=1,ntra
	3134	c do 2110 k=1,nd ! oct3
	3135	c do 2120 i=1,ncum
	3136	c ftra1(idcum(i),k,j)=ftra(i,k,j)
	3137	c 2120 continue
	3138	c 2110 continue
	3139	c 2100 continue
	3140	do j=1,nd
	3141	do k=1,nd
	3142	do i=1,ncum
	3143	phi1(idcum(i),k,j)=phi(i,k,j)
	3144	end do
	3145	end do
	3146	end do
	3147
	3148	return
	3149	end
	3150

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

Download in other formats:

Original Format