Context Navigation

convect3.F90 @ 5473

Last change on this file since 5473 was 5160, checked in by abarral, 6 months ago
Put .h into modules
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 45.9 KB

Rev	Line
[524]	1	! $Header$
[766]	2
[1992]	3	SUBROUTINE convect3(dtime, epmax, ok_adj, t1, r1, rs, u, v, tra, p, ph, nd, &
[5144]	4	ndp1, nl, ntra, delt, iflag, ft, fr, fu, fv, ftra, precip, icb, inb, &
	5	upwd, dnwd, dnwd0, sig, w0, mike, mke, ma, ments, qents, tps, tls, sigij, &
	6	cape, tvp, pbase, buoybase, & ! ccc *
	7	! DTVPDT1,DTVPDQ1,DPLCLDT,DPLCLDR)
	8	dtvpdt1, dtvpdq1, dplcldt, dplcldr, & ! sbl
	9	ft2, fr2, fu2, fv2, wd, qcond, qcondc) ! sbl
[524]	10
[1992]	11	! * THE PARAMETER NA SHOULD IN GENERAL EQUAL ND *
[524]	12
[1992]	13	! #################################################################
	14	! Fleur Introduction des traceurs dans convect3 le 6 juin 200
	15	! #################################################################
	16	USE dimphy
[2320]	17	USE infotrac_phy, ONLY: nbtr
[5144]	18	USE lmdz_yomcst
	19
[2197]	20	IMPLICIT NONE
[1992]	21	INTEGER na
[5144]	22	PARAMETER (na = 60)
[524]	23
[1992]	24	REAL deltac ! cld
[5144]	25	PARAMETER (deltac = 0.01) ! cld
[524]	26
[1992]	27	INTEGER nent(na)
	28	INTEGER nd, ndp1, nl, ntra, iflag, icb, inb
	29	REAL dtime, epmax, delt, precip, cape
	30	REAL dplcldt, dplcldr
	31	REAL t1(nd), r1(nd), rs(nd), u(nd), v(nd), tra(nd, ntra)
	32	REAL p(nd), ph(ndp1)
	33	REAL ft(nd), fr(nd), fu(nd), fv(nd), ftra(nd, ntra)
	34	REAL sig(nd), w0(nd)
	35	REAL uent(na, na), vent(na, na), traent(na, na, nbtr), tratm(na)
	36	REAL up(na), vp(na), trap(na, nbtr)
	37	REAL m(na), mp(na), ment(na, na), qent(na, na), elij(na, na)
	38	REAL sij(na, na), tvp(na), tv(na), water(na)
	39	REAL rp(na), ep(na), th(na), wt(na), evap(na), clw(na)
	40	REAL sigp(na), b(na), tp(na), cpn(na)
	41	REAL lv(na), lvcp(na), h(na), hp(na), gz(na)
	42	REAL t(na), rr(na)
[524]	43
[1992]	44	REAL ft2(nd), fr2(nd), fu2(nd), fv2(nd) ! added sbl
	45	REAL u1(nd), v1(nd) ! added sbl
[524]	46
[1992]	47	REAL buoy(na) ! Lifted parcel buoyancy
	48	REAL dtvpdt1(nd), dtvpdq1(nd) ! Derivatives of parcel virtual
	49	! temperature wrt T1 and Q1
	50	REAL clw_new(na), qi(na)
[524]	51
[1992]	52	REAL wd, betad ! for gust factor (sb)
	53	REAL qcondc(nd) ! interface cld param (sb)
	54	REAL qcond(nd), nqcond(na), wa(na), maa(na), siga(na), axc(na) ! cld
[524]	55
[1992]	56	LOGICAL ice_conv, ok_adj
[5144]	57	PARAMETER (ice_conv = .TRUE.)
[524]	58
[1992]	59	! ccccccccccccccccccccccccccccccccccccccccccccc
	60	! declaration des variables a sortir
	61	! cccccccccccccccccccccccccccccccccccccccccccc
	62	REAL mke(nd)
	63	REAL mike(nd)
	64	REAL ma(nd)
	65	REAL tps(nd) !temperature dans les ascendances non diluees
	66	REAL tls(nd) !temperature potentielle
	67	REAL ments(nd, nd)
	68	REAL qents(nd, nd)
	69	REAL sigij(klev, klev)
	70	REAL pbase ! pressure at the cloud base level
	71	REAL buoybase ! buoyancy at the cloud base level
	72	! ccccccccccccccccccccccccccccccccccccccccccccc
[524]	73
[5144]	74	REAL :: cpv, cl, cpvmcl, eps, alv0, rdcp, pbcrit, ptcrit, sigd, spfac
	75	REAL :: tau, beta, alpha, dtcrit, dtovsh, ahm, rm, um, vm, dphinv
	76	REAL :: a2, x, tvx, tvy, plcl, pden, dpbase, tvpbase, tvbase, tdif
	77	REAL :: ath1, ath, delti, deltap, dcape, dlnp, sigold, dtmin, fac, w
	78	REAL :: amu, rti, cpd, bf2, anum, denom, dei, altem, cwat, stemp, qp
	79	REAL :: scrit, alt, smax, asij, wgh, sjmax, sjmin, smid, delp, delm
	80	REAL :: asum, bsum, csum, wflux, tinv, wdtrain, awat, afac, afac1, afac2
	81	REAL :: bfac, pr1, pr2, sigt, b6, c6, revap, tevap, delth, amfac, amp2
	82	REAL :: xf, tf, af, bf, fac2, ur, sru, d, ampmax, dpinv, am, amde, cpinv
	83	REAL :: amp1, ad, rat, ax, bx, cx, dx, ex, dsum
	84	INTEGER :: nk, i, j, nopt, jn, k, im, jm, n
[1992]	85
	86	REAL dnwd0(nd) ! precipitation driven unsaturated downdraft flux
	87	REAL dnwd(nd), dn1 ! in-cloud saturated downdraft mass flux
	88	REAL upwd(nd), up1 ! in-cloud saturated updraft mass flux
	89
	90	! * ASSIGN VALUES OF THERMODYNAMIC CONSTANTS *
	91	! * THESE SHOULD BE CONSISTENT WITH *
	92	! * THOSE USED IN CALLING PROGRAM *
	93	! * NOTE: THESE ARE ALSO SPECIFIED IN SUBROUTINE TLIFT *
	94
	95	! sb CPD=1005.7
	96	! sb CPV=1870.0
	97	! sb CL=4190.0
	98	! sb CPVMCL=CL-CPV
	99	! sb RV=461.5
	100	! sb RD=287.04
	101	! sb EPS=RD/RV
	102	! sb ALV0=2.501E6
	103	! cccccccccccccccccccccc
	104	! constantes coherentes avec le modele du Centre Europeen
	105	! sb RD = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 28.9644
	106	! sb RV = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 18.0153
	107	! sb CPD = 3.5 * RD
	108	! sb CPV = 4.0 * RV
	109	! sb CL = 4218.0
	110	! sb CPVMCL=CL-CPV
	111	! sb EPS=RD/RV
	112	! sb ALV0=2.5008E+06
	113	! ccccccccccccccccccccc
	114	! on utilise les constantes thermo du Centre Europeen: (SB)
	115
	116	cpd = rcpd
	117	cpv = rcpv
	118	cl = rcw
	119	cpvmcl = cl - cpv
[5144]	120	eps = rd / rv
[1992]	121	alv0 = rlvtt
	122
	123	nk = 1 ! origin level of the lifted parcel
	124
	125	! ccccccccccccccccccccc
	126
	127	! * INITIALIZE OUTPUT ARRAYS AND PARAMETERS *
	128
	129	DO i = 1, nd
	130	ft(i) = 0.0
	131	fr(i) = 0.0
	132	fu(i) = 0.0
	133	fv(i) = 0.0
	134
	135	ft2(i) = 0.0
	136	fr2(i) = 0.0
	137	fu2(i) = 0.0
	138	fv2(i) = 0.0
	139
	140	DO j = 1, ntra
	141	ftra(i, j) = 0.0
	142	END DO
	143
	144	qcondc(i) = 0.0 ! cld
	145	qcond(i) = 0.0 ! cld
	146	nqcond(i) = 0.0 ! cld
	147
	148	t(i) = t1(i)
	149	rr(i) = r1(i)
	150	u1(i) = u(i) ! added sbl
	151	v1(i) = v(i) ! added sbl
	152	END DO
	153	DO i = 1, nl
[5144]	154	rdcp = (rd * (1. - rr(i)) + rr(i) * rv) / (cpd * (1. - rr(i)) + rr(i) * cpv)
	155	th(i) = t(i) * (1000.0 / p(i))**rdcp
[1992]	156	END DO
	157
	158	! ************************************************************
	159	! * CALCUL DES TEMPERATURES POTENTIELLES A SORTIR
	160	! ************************************************************
	161	DO i = 1, nd
[5144]	162	rdcp = (rd * (1. - rr(i)) + rr(i) * rv) / (cpd * (1. - rr(i)) + rr(i) * cpv)
[1992]	163
[5144]	164	tls(i) = t(i) * (1000.0 / p(i))**rdcp
[1992]	165	END DO
	166
	167
	168
	169
	170	! ***********************************************************
	171
	172	precip = 0.0
	173	wd = 0.0 ! sb
	174	iflag = 1
	175
	176	! * SPECIFY PARAMETERS *
	177	! * PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE *
	178	! * PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO *
	179	! * PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. *
	180	! * EFFICIENCY IS ASSUMED TO BE UNITY *
	181	! * SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *
	182	! * SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *
	183	! * OF CLOUD *
	184	! * ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF *
	185	! * APPROACH TO QUASI-EQUILIBRIUM *
	186	! * (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) *
	187	! * (BETA MUST BE LESS THAN OR EQUAL TO 1) *
	188	! * DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE *
	189	! * APPROACH TO QUASI-EQUILIBRIUM *
	190	! * IT MUST BE LESS THAN 0 *
	191
	192	pbcrit = 150.0
	193	ptcrit = 500.0
	194	sigd = 0.01
	195	spfac = 0.15
	196	! sb:
	197	! EPMAX=0.993 ! precip efficiency less than unity
	198	! EPMAX=1. ! precip efficiency less than unity
	199
	200	! jyg
	201	! CC BETA=0.96
	202	! Beta is now expressed as a function of the characteristic time
	203	! of the convective process.
	204	! CC Old value : TAU = 15000. !(for dtime = 600.s)
	205	! CC Other value (inducing little change) :TAU = 8000.
	206	tau = 8000.
[5144]	207	beta = 1. - dtime / tau
[1992]	208	! jyg
	209	! CC ALPHA=1.0
[5144]	210	alpha = 1.5E-3 * dtime / tau
[1992]	211	! Increase alpha in order to compensate W decrease
[5144]	212	alpha = alpha * 1.5
[1992]	213
	214	! jyg (voir CONVECT 3)
	215	! CC DTCRIT=-0.2
	216	dtcrit = -2.
	217	! gf&jyg
	218	! CC DT pour l'overshoot.
	219	dtovsh = -0.2
	220
	221
	222	! * INCREMENT THE COUNTER *
	223
	224	sig(nd) = sig(nd) + 1.0
	225	sig(nd) = amin1(sig(nd), 12.1)
	226
	227	! * IF NOPT IS AN INTEGER OTHER THAN 0, CONVECT *
	228	! * RETURNS ARRAYS T AND R THAT MAY HAVE BEEN *
	229	! * ALTERED BY DRY ADIABATIC ADJUSTMENT; OTHERWISE *
	230	! * THE RETURNED ARRAYS ARE UNALTERED. *
	231
	232	nopt = 0
[5113]	233	! NOPT=1 ! sbl
[1992]	234
	235	! * PERFORM DRY ADIABATIC ADJUSTMENT *
	236
	237	! * DO NOT BYPASS THIS EVEN IF THE CALLING PROGRAM HAS A *
	238	! * BOUNDARY LAYER SCHEME !!! *
	239
	240	IF (ok_adj) THEN ! added sbl
	241
	242	DO i = nl - 1, 1, -1
	243	jn = 0
	244	DO j = i + 1, nl
	245	IF (th(j)<th(i)) jn = j
	246	END DO
	247	IF (jn==0) GO TO 30
	248	ahm = 0.0
	249	rm = 0.0
	250	um = 0.0
	251	vm = 0.0
	252	DO k = 1, ntra
	253	tratm(k) = 0.0
	254	END DO
	255	DO j = i, jn
[5144]	256	ahm = ahm + (cpd * (1. - rr(j)) + rr(j) * cpv) * t(j) * (ph(j) - ph(j + 1))
	257	rm = rm + rr(j) * (ph(j) - ph(j + 1))
	258	um = um + u(j) * (ph(j) - ph(j + 1))
	259	vm = vm + v(j) * (ph(j) - ph(j + 1))
[1992]	260	DO k = 1, ntra
[5144]	261	tratm(k) = tratm(k) + tra(j, k) * (ph(j) - ph(j + 1))
[1992]	262	END DO
	263	END DO
[5144]	264	dphinv = 1. / (ph(i) - ph(jn + 1))
	265	rm = rm * dphinv
	266	um = um * dphinv
	267	vm = vm * dphinv
[1992]	268	DO k = 1, ntra
[5144]	269	tratm(k) = tratm(k) * dphinv
[1992]	270	END DO
	271	a2 = 0.0
	272	DO j = i, jn
	273	rr(j) = rm
	274	u(j) = um
	275	v(j) = vm
	276	DO k = 1, ntra
	277	tra(j, k) = tratm(k)
	278	END DO
[5144]	279	rdcp = (rd * (1. - rr(j)) + rr(j) * rv) / (cpd * (1. - rr(j)) + rr(j) * cpv)
	280	x = (0.001 * p(j))**rdcp
[1992]	281	t(j) = x
[5144]	282	a2 = a2 + (cpd * (1. - rr(j)) + rr(j) * cpv) * x * (ph(j) - ph(j + 1))
[1992]	283	END DO
	284	DO j = i, jn
[5144]	285	th(j) = ahm / a2
	286	t(j) = t(j) * th(j)
[1992]	287	END DO
[5144]	288	30 END DO
[1992]	289
	290	END IF ! added sbl
	291
	292	! * RESET INPUT ARRAYS IF ok_adj 0 *
	293
	294	IF (ok_adj) THEN
	295	DO i = 1, nd
	296
[5144]	297	ft2(i) = (t(i) - t1(i)) / delt ! sbl
	298	fr2(i) = (rr(i) - r1(i)) / delt ! sbl
	299	fu2(i) = (u(i) - u1(i)) / delt ! sbl
	300	fv2(i) = (v(i) - v1(i)) / delt ! sbl
[1992]	301
[5113]	302	! T1(I)=T(I) ! commente sbl
	303	! R1(I)=RR(I) ! commente sbl
[1992]	304	END DO
	305	END IF
	306
	307	! *** CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY AND STATIC ENERGY
	308
	309	gz(1) = 0.0
[5144]	310	cpn(1) = cpd * (1. - rr(1)) + rr(1) * cpv
	311	h(1) = t(1) * cpn(1)
[1992]	312	DO i = 2, nl
[5144]	313	tvx = t(i) * (1. + rr(i) / eps - rr(i))
	314	tvy = t(i - 1) * (1. + rr(i - 1) / eps - rr(i - 1))
	315	gz(i) = gz(i - 1) + 0.5 * rd * (tvx + tvy) * (p(i - 1) - p(i)) / ph(i)
	316	cpn(i) = cpd * (1. - rr(i)) + cpv * rr(i)
	317	h(i) = t(i) * cpn(i) + gz(i)
[1992]	318	END DO
	319
	320	! * CALCULATE LIFTED CONDENSATION LEVEL OF AIR AT LOWEST MODEL LEVEL *
	321	! * (WITHIN 0.2% OF FORMULA OF BOLTON, MON. WEA. REV.,1980) *
	322
	323	IF (t(1)<250.0 .OR. rr(1)<=0.0) THEN
	324	iflag = 0
[5103]	325	! sb3d PRINT*,'je suis passe par 366'
[1992]	326	RETURN
	327	END IF
	328
[5103]	329	! jyg1 Utilisation de la SUBROUTINE CLIFT
[1992]	330	! C RH=RR(1)/RS(1)
	331	! C CHI=T(1)/(1669.0-122.0*RH-T(1))
	332	! C PLCL=P(1)(RH*CHI)
	333	CALL clift(p(1), t(1), rr(1), rs(1), plcl, dplcldt, dplcldr)
	334	! jyg2
	335	! sb3d PRINT *,' em_plcl,p1,t1,r1,rs1,rh '
	336	! sb3d $ ,PLCL,P(1),T(1),RR(1),RS(1),RH
	337
	338	IF (plcl<200.0 .OR. plcl>=2000.0) THEN
	339	iflag = 2
	340	RETURN
	341	END IF
	342	! jyg1
	343	! Essais de modification de ICB
	344
	345	! * CALCULATE FIRST LEVEL ABOVE LCL (=ICB) *
	346
	347	! C ICB=NL-1
	348	! C DO 50 I=2,NL-1
	349	! C IF(P(I).LT.PLCL)THEN
	350	! C ICB=MIN(ICB,I) ! ICB sup ou egal a 2
	351	! C END IF
	352	! C 50 CONTINUE
	353	! C IF(ICB.EQ.(NL-1))THEN
	354	! C IFLAG=3
	355	! C RETURN
	356	! C END IF
	357
	358	! * CALCULATE LAYER CONTAINING LCL (=ICB) *
	359
	360	icb = nl - 1
	361	! sb DO 50 I=2,NL-1
	362	DO i = 3, nl - 1 ! modif sb pour que ICB soit sup/egal a 2
	363	! la modification consiste a comparer PLCL a PH et non a P:
	364	! ICB est defini par : PH(ICB)<PLCL<PH(ICB-!)
	365	IF (ph(i)<plcl) THEN
	366	icb = min(icb, i)
	367	END IF
	368	END DO
[5144]	369	IF (icb==(nl - 1)) THEN
[1992]	370	iflag = 3
	371	RETURN
	372	END IF
	373	icb = icb - 1 ! ICB sup ou egal a 2
	374	! jyg2
	375
	376
	377
	378	! *** SUBROUTINE TLIFT CALCULATES PART OF THE LIFTED PARCEL VIRTUAL
	379	! ***
	380	! *** TEMPERATURE, THE ACTUAL TEMPERATURE AND THE ADIABATIC
	381	! ***
	382	! *** LIQUID WATER CONTENT
	383	! ***
	384
	385
	386	! jyg1
	387	! make sure that "Cloud base" seen by TLIFT is actually the
	388	! fisrt level where adiabatic ascent is saturated
	389	IF (plcl>p(icb)) THEN
	390	! sb CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB,TVP,TP,CLW,ND,NL)
	391	CALL tlift(p, t, rr, rs, gz, plcl, icb, nk, tvp, tp, clw, nd, nl, &
[5144]	392	dtvpdt1, dtvpdq1)
[1992]	393	ELSE
	394	! sb CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB+1,TVP,TP,CLW,ND,NL)
[5144]	395	CALL tlift(p, t, rr, rs, gz, plcl, icb + 1, nk, tvp, tp, clw, nd, nl, &
	396	dtvpdt1, dtvpdq1)
[1992]	397	END IF
	398	! jyg2
	399
	400	! *****************************************************************************
	401	! *** SORTIE DE LA TEMPERATURE DE L ASCENDANCE NON DILUE
	402	! *****************************************************************************
	403	DO i = 1, nd
	404	tps(i) = tp(i)
	405	END DO
	406
	407
	408	! *****************************************************************************
	409
	410
	411	! * SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF *
	412	! * PRECIPITATION FALLING OUTSIDE OF CLOUD *
	413	! * THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I) *
	414
	415	DO i = 1, nl
	416	pden = ptcrit - pbcrit
	417
	418	! jyg
	419	! cc EP(I)=(P(ICB)-P(I)-PBCRIT)/PDEN
	420	! sb EP(I)=(PLCL-P(I)-PBCRIT)/PDEN
[5144]	421	ep(i) = (plcl - p(i) - pbcrit) / pden * epmax ! sb
[1992]	422
	423	ep(i) = amax1(ep(i), 0.0)
	424	! sb EP(I)=AMIN1(EP(I),1.0)
	425	ep(i) = amin1(ep(i), epmax) ! sb
	426	sigp(i) = spfac
	427
	428	! * CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL *
	429	! * VIRTUAL TEMPERATURE *
	430
[5144]	431	tv(i) = t(i) * (1. + rr(i) / eps - rr(i))
[1992]	432	! cd1
	433	! . Keep all liquid water in lifted parcel (-> adiabatic CAPE)
	434
	435	! cc TVP(I)=TVP(I)-TP(I)(RR(1)-EP(I)CLW(I))
	436	! !!! sb TVP(I)=TVP(I)-TP(I)*RR(1) ! calcule dans tlift
	437	! cd2
	438
	439	! *** Calculate first estimate of buoyancy
	440
	441	buoy(i) = tvp(i) - tv(i)
	442	END DO
	443
	444	! *** Set Cloud Base Buoyancy at (Plcl+DPbase) level buoyancy
	445
	446	dpbase = -40. !That is 400m above LCL
	447	pbase = plcl + dpbase
[5144]	448	tvpbase = tvp(icb) * (pbase - p(icb + 1)) / (p(icb) - p(icb + 1)) + &
	449	tvp(icb + 1) * (p(icb) - pbase) / (p(icb) - p(icb + 1))
	450	tvbase = tv(icb) * (pbase - p(icb + 1)) / (p(icb) - p(icb + 1)) + &
	451	tv(icb + 1) * (p(icb) - pbase) / (p(icb) - p(icb + 1))
[1992]	452
	453	! test sb:
[5116]	454	! @ WRITE(,) '++++++++++++++++++++++++++++++++++++++++'
	455	! @ WRITE(,) 'plcl,dpbas,tvpbas,tvbas,tvp(icb),tvp(icb1)
[1992]	456	! @ : ,tv(icb),tv(icb1)'
[5116]	457	! @ WRITE(,) plcl,dpbase,tvpbase,tvbase,tvp(icb)
[1992]	458	! @ L ,tvp(icb+1),tv(icb),tv(icb+1)
[5116]	459	! @ WRITE(,) '++++++++++++++++++++++++++++++++++++++++'
[1992]	460	! fin test sb
	461	buoybase = tvpbase - tvbase
	462
	463	! C Set buoyancy = BUOYBASE for all levels below BASE.
	464	! C For safety, set : BUOY(ICB) = BUOYBASE
	465	DO i = icb, nl
	466	IF (p(i)>=pbase) THEN
	467	buoy(i) = buoybase
	468	END IF
	469	END DO
	470	buoy(icb) = buoybase
	471
[5160]	472	! sb3d PRINT *,'buoybase,tvp_tv,tvpbase,tvbase,pbase,plcl'
[1992]	473	! sb3d $, buoybase,tvp(icb)-tv(icb),tvpbase,tvbase,pbase,plcl
[5160]	474	! sb3d PRINT *,'TVP ',(tvp(i),i=1,nl)
	475	! sb3d PRINT *,'TV ',(tv(i),i=1,nl)
	476	! sb3d PRINT *, 'P ',(p(i),i=1,nl)
	477	! sb3d PRINT *,'ICB ',icb
[1992]	478	! test sb:
[5116]	479	! @ WRITE(,) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'
	480	! @ WRITE(,) 'icb,icbs,inb,buoybase:'
	481	! @ WRITE(,) icb,icb+1,inb,buoybase
	482	! @ WRITE(,) 'k,tvp,tv,tp,buoy,ep: '
[1992]	483	! @ do k=1,nl
[5116]	484	! @ WRITE(,) k,tvp(k),tv(k),tp(k),buoy(k),ep(k)
[1992]	485	! @ enddo
[5116]	486	! @ WRITE(,) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'
[1992]	487	! fin test sb
	488
	489
	490
	491	! * MAKE SURE THAT COLUMN IS DRY ADIABATIC BETWEEN THE SURFACE *
	492	! * AND CLOUD BASE, AND THAT LIFTED AIR IS POSITIVELY BUOYANT *
	493	! * AT CLOUD BASE *
	494	! * IF NOT, RETURN TO CALLING PROGRAM AFTER RESETTING *
	495	! * SIG(I) AND W0(I) *
	496
	497	! jyg
	498	! CC TDIF=TVP(ICB)-TV(ICB)
	499	tdif = buoy(icb)
	500	ath1 = th(1)
	501	! jyg
	502	! CC ATH=TH(ICB-1)-1.0
[5144]	503	ath = th(icb - 1) - 5.0
[1992]	504	! ATH=0. ! ajout sb
	505	! IF (ICB.GT.1) ATH=TH(ICB-1)-5.0 ! modif sb
	506	IF (tdif<dtcrit .OR. ath>ath1) THEN
	507	DO i = 1, nl
[5144]	508	sig(i) = beta * sig(i) - 2. * alpha * tdif * tdif
[1992]	509	sig(i) = amax1(sig(i), 0.0)
[5144]	510	w0(i) = beta * w0(i)
[1992]	511	END DO
	512	iflag = 0
	513	RETURN
	514	END IF
	515
	516
	517
	518	! *** IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY
	519	! ***
	520	! *** NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
	521	! ***
	522
	523	DO i = 1, nl
	524	hp(i) = h(i)
	525	wt(i) = 0.001
	526	rp(i) = rr(i)
	527	up(i) = u(i)
	528	vp(i) = v(i)
	529	DO j = 1, ntra
	530	trap(i, j) = tra(i, j)
	531	END DO
	532	nent(i) = 0
	533	water(i) = 0.0
	534	evap(i) = 0.0
	535	b(i) = 0.0
	536	mp(i) = 0.0
	537	m(i) = 0.0
[5144]	538	lv(i) = alv0 - cpvmcl * (t(i) - 273.15)
	539	lvcp(i) = lv(i) / cpn(i)
[1992]	540	DO j = 1, nl
	541	qent(i, j) = rr(j)
	542	elij(i, j) = 0.0
	543	ment(i, j) = 0.0
	544	sij(i, j) = 0.0
	545	uent(i, j) = u(j)
	546	vent(i, j) = v(j)
	547	DO k = 1, ntra
	548	traent(i, j, k) = tra(j, k)
	549	END DO
	550	END DO
	551	END DO
	552
[5144]	553	delti = 1.0 / delt
[1992]	554
	555	! * FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S *
	556	! * LEVEL OF NEUTRAL BUOYANCY *
	557
	558	inb = nl - 1
	559	DO i = icb, nl - 1
	560	! jyg
	561	! CC IF((TVP(I)-TV(I)).LT.DTCRIT)THEN
	562	IF (buoy(i)<dtovsh) THEN
	563	inb = min(inb, i)
	564	END IF
	565	END DO
	566
	567
	568
	569
	570
	571	! * RESET SIG(I) AND W0(I) FOR I>INB AND I<ICB *
	572
[5144]	573	IF (inb<(nl - 1)) THEN
[1992]	574	DO i = inb + 1, nl - 1
	575	! jyg
	576	! CC SIG(I)=BETASIG(I)-2.0E-4ALPHA(TV(INB)-TVP(INB))
	577	! CC 1 ABS(TV(INB)-TVP(INB))
[5144]	578	sig(i) = beta * sig(i) + 2. * alpha * buoy(inb) * abs(buoy(inb))
[1992]	579	sig(i) = amax1(sig(i), 0.0)
[5144]	580	w0(i) = beta * w0(i)
[1992]	581	END DO
	582	END IF
	583	DO i = 1, icb
	584	! jyg
	585	! CC SIG(I)=BETASIG(I)-2.0E-4ALPHA(TV(ICB)-TVP(ICB))
	586	! CC 1 (TV(ICB)-TVP(ICB))
[5144]	587	sig(i) = beta * sig(i) - 2. * alpha * buoy(icb) * buoy(icb)
[1992]	588	sig(i) = amax1(sig(i), 0.0)
[5144]	589	w0(i) = beta * w0(i)
[1992]	590	END DO
	591
	592	! * RESET FRACTIONAL AREAS OF UPDRAFTS AND W0 AT INITIAL TIME *
	593	! * AND AFTER 10 TIME STEPS OF NO CONVECTION *
	594
	595	IF (sig(nd)<1.5 .OR. sig(nd)>12.0) THEN
	596	DO i = 1, nl - 1
	597	sig(i) = 0.0
	598	w0(i) = 0.0
	599	END DO
	600	END IF
	601
	602	! * CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL *
	603
	604	DO i = icb, inb
[5144]	605	hp(i) = h(1) + (lv(i) + (cpd - cpv) * t(i)) * ep(i) * clw(i)
[1992]	606	END DO
	607
	608	! * CALCULATE CONVECTIVE AVAILABLE POTENTIAL ENERGY (CAPE), *
	609	! * VERTICAL VELOCITY (W), FRACTIONAL AREA COVERED BY *
	610	! * UNDILUTE UPDRAFT (SIG), AND UPDRAFT MASS FLUX (M) *
	611
	612	cape = 0.0
	613
	614	DO i = icb + 1, inb
	615	! jyg1
	616	! CC CAPE=CAPE+RD(TVP(I-1)-TV(I-1))(PH(I-1)-PH(I))/P(I-1)
	617	! CC DCAPE=RDBUOY(I-1)(PH(I-1)-PH(I))/P(I-1)
	618	! CC DLNP=(PH(I-1)-PH(I))/P(I-1)
	619	! The interval on which CAPE is computed starts at PBASE :
[5144]	620	deltap = min(pbase, ph(i - 1)) - min(pbase, ph(i))
	621	cape = cape + rd * buoy(i - 1) * deltap / p(i - 1)
	622	dcape = rd * buoy(i - 1) * deltap / p(i - 1)
	623	dlnp = deltap / p(i - 1)
[1992]	624	! jyg2
[5160]	625	! sb3d PRINT *,'buoy,dlnp,dcape,cape',buoy(i-1),dlnp,dcape,cape
[1992]	626	! test sb:
[5116]	627	! @ WRITE(,) '############################################'
	628	! @ WRITE(,) 'cape,rrd,buoy,deltap,p,pbase,ph:'
[1992]	629	! @ : ,cape,rd,buoy(i-1),deltap,p(i-1),pbase,ph(i)
[5116]	630	! @ WRITE(,) '############################################'
[1992]	631
	632	! fin test sb
	633	cape = amax1(0.0, cape)
	634
	635	sigold = sig(i)
	636	dtmin = 100.0
	637	DO j = icb, i - 1
	638	! jyg
	639	! CC DTMIN=AMIN1(DTMIN,(TVP(J)-TV(J)))
	640	dtmin = amin1(dtmin, buoy(j))
	641	END DO
[5160]	642	! sb3d PRINT *, 'DTMIN, BETA, ALPHA, SIG = ',DTMIN,BETA,ALPHA,SIG(I)
[5144]	643	sig(i) = beta * sig(i) + alpha * dtmin * abs(dtmin)
[1992]	644	sig(i) = amax1(sig(i), 0.0)
	645	sig(i) = amin1(sig(i), 0.01)
[5144]	646	fac = amin1(((dtcrit - dtmin) / dtcrit), 1.0)
[1992]	647	! jyg
	648	! C Essais de reduction de la vitesse
	649	! C FAC = FAC*.5
	650
[5144]	651	w = (1. - beta) * fac * sqrt(cape) + beta * w0(i)
	652	amu = 0.5 * (sig(i) + sigold) * w
	653	m(i) = amu * 0.007 * p(i) * (ph(i) - ph(i + 1)) / tv(i)
[1992]	654
	655	! --------- test sb:
[5116]	656	! WRITE(,) '############################################'
	657	! WRITE(,) 'k,amu,buoy(k-1),deltap,w,beta,fac,cape,w0(k)'
	658	! WRITE(,) i,amu,buoy(i-1),deltap
[1992]	659	! : ,w,beta,fac,cape,w0(i)
[5116]	660	! WRITE(,) '############################################'
[1992]	661	! ---------
	662
	663	w0(i) = w
	664	END DO
[5144]	665	w0(icb) = 0.5 * w0(icb + 1)
	666	m(icb) = 0.5 * m(icb + 1) * (ph(icb) - ph(icb + 1)) / (ph(icb + 1) - ph(icb + 2))
	667	sig(icb) = sig(icb + 1)
	668	sig(icb - 1) = sig(icb)
[1992]	669	! jyg1
[5160]	670	! sb3d PRINT *, 'Cloud base, c. top, CAPE',ICB,INB,cape
	671	! sb3d PRINT *, 'SIG ',(sig(i),i=1,inb)
	672	! sb3d PRINT *, 'W ',(w0(i),i=1,inb)
	673	! sb3d PRINT *, 'M ',(m(i), i=1,inb)
	674	! sb3d PRINT *, 'Dt1 ',(tvp(i)-tv(i),i=1,inb)
	675	! sb3d PRINT *, 'Dt_vrai ',(buoy(i),i=1,inb)
[1992]	676	! jyg2
	677
	678	! * CALCULATE ENTRAINED AIR MASS FLUX (MENT), TOTAL WATER MIXING *
	679	! * RATIO (QENT), TOTAL CONDENSED WATER (ELIJ), AND MIXING *
	680	! * FRACTION (SIJ) *
	681
	682	DO i = icb, inb
[5144]	683	rti = rr(1) - ep(i) * clw(i)
[1992]	684	DO j = icb - 1, inb
[5144]	685	bf2 = 1. + lv(j) * lv(j) * rs(j) / (rv * t(j) * t(j) * cpd)
	686	anum = h(j) - hp(i) + (cpv - cpd) * t(j) * (rti - rr(j))
	687	denom = h(i) - hp(i) + (cpd - cpv) * (rr(i) - rti) * t(j)
[1992]	688	dei = denom
	689	IF (abs(dei)<0.01) dei = 0.01
[5144]	690	sij(i, j) = anum / dei
[1992]	691	sij(i, i) = 1.0
[5144]	692	altem = sij(i, j) * rr(i) + (1. - sij(i, j)) * rti - rs(j)
	693	altem = altem / bf2
	694	cwat = clw(j) * (1. - ep(j))
[1992]	695	stemp = sij(i, j)
	696	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
[5144]	697	anum = anum - lv(j) * (rti - rs(j) - cwat * bf2)
	698	denom = denom + lv(j) * (rr(i) - rti)
[1992]	699	IF (abs(denom)<0.01) denom = 0.01
[5144]	700	sij(i, j) = anum / denom
	701	altem = sij(i, j) * rr(i) + (1. - sij(i, j)) * rti - rs(j)
	702	altem = altem - (bf2 - 1.) * cwat
[524]	703	END IF
	704
[5144]	705	IF (sij(i, j)>0.0 .AND. sij(i, j)<0.95) THEN
	706	qent(i, j) = sij(i, j) * rr(i) + (1. - sij(i, j)) * rti
	707	uent(i, j) = sij(i, j) * u(i) + (1. - sij(i, j)) * u(nk)
	708	vent(i, j) = sij(i, j) * v(i) + (1. - sij(i, j)) * v(nk)
[1992]	709	DO k = 1, ntra
[5144]	710	traent(i, j, k) = sij(i, j) * tra(i, k) + (1. - sij(i, j)) * tra(nk, k)
[1992]	711	END DO
	712	elij(i, j) = altem
[5144]	713	elij(i, j) = amax1(0.0, elij(i, j))
	714	ment(i, j) = m(i) / (1. - sij(i, j))
[1992]	715	nent(i) = nent(i) + 1
[524]	716	END IF
[5144]	717	sij(i, j) = amax1(0.0, sij(i, j))
	718	sij(i, j) = amin1(1.0, sij(i, j))
[1992]	719	END DO
[524]	720
[1992]	721	! *** IF NO AIR CAN ENTRAIN AT LEVEL I ASSUME THAT UPDRAFT DETRAINS
	722	! ***
	723	! *** AT THAT LEVEL AND CALCULATE DETRAINED AIR FLUX AND PROPERTIES
	724	! ***
[524]	725
[1992]	726	IF (nent(i)==0) THEN
	727	ment(i, i) = m(i)
[5144]	728	qent(i, i) = rr(nk) - ep(i) * clw(i)
[1992]	729	uent(i, i) = u(nk)
	730	vent(i, i) = v(nk)
	731	DO j = 1, ntra
	732	traent(i, i, j) = tra(nk, j)
	733	END DO
	734	elij(i, i) = clw(i)
	735	sij(i, i) = 1.0
	736	END IF
[524]	737
[1992]	738	DO j = icb - 1, inb
	739	sigij(i, j) = sij(i, j)
	740	END DO
	741
	742	END DO
	743
	744	! * NORMALIZE ENTRAINED AIR MASS FLUXES TO REPRESENT EQUAL *
	745	! * PROBABILITIES OF MIXING *
	746
	747	DO i = icb, inb
	748	IF (nent(i)/=0) THEN
[5144]	749	qp = rr(1) - ep(i) * clw(i)
	750	anum = h(i) - hp(i) - lv(i) * (qp - rs(i)) + (cpv - cpd) * t(i) * (qp - rr(i))
	751	denom = h(i) - hp(i) + lv(i) * (rr(i) - qp) + (cpd - cpv) * t(i) * (rr(i) - qp)
[1992]	752	IF (abs(denom)<0.01) denom = 0.01
[5144]	753	scrit = anum / denom
	754	alt = qp - rs(i) + scrit * (rr(i) - qp)
[1992]	755	IF (scrit<=0.0 .OR. alt<=0.0) scrit = 1.0
	756	smax = 0.0
	757	asij = 0.0
	758	DO j = inb, icb - 1, -1
[5144]	759	IF (sij(i, j)>1.0E-16 .AND. sij(i, j)<0.95) THEN
[1992]	760	wgh = 1.0
	761	IF (j>i) THEN
[5144]	762	sjmax = amax1(sij(i, j + 1), smax)
[1992]	763	sjmax = amin1(sjmax, scrit)
[5144]	764	smax = amax1(sij(i, j), smax)
	765	sjmin = amax1(sij(i, j - 1), smax)
[1992]	766	sjmin = amin1(sjmin, scrit)
[5144]	767	IF (sij(i, j)<(smax - 1.0E-16)) wgh = 0.0
	768	smid = amin1(sij(i, j), scrit)
[1992]	769	ELSE
[5144]	770	sjmax = amax1(sij(i, j + 1), scrit)
	771	smid = amax1(sij(i, j), scrit)
[1992]	772	sjmin = 0.0
[5144]	773	IF (j>1) sjmin = sij(i, j - 1)
[1992]	774	sjmin = amax1(sjmin, scrit)
	775	END IF
[5144]	776	delp = abs(sjmax - smid)
	777	delm = abs(sjmin - smid)
	778	asij = asij + wgh * (delp + delm)
	779	ment(i, j) = ment(i, j) * (delp + delm) * wgh
[524]	780	END IF
[1992]	781	END DO
	782	asij = amax1(1.0E-16, asij)
[5144]	783	asij = 1.0 / asij
[1992]	784	DO j = icb - 1, inb
[5144]	785	ment(i, j) = ment(i, j) * asij
[1992]	786	END DO
	787	asum = 0.0
	788	bsum = 0.0
	789	DO j = icb - 1, inb
	790	asum = asum + ment(i, j)
[5144]	791	ment(i, j) = ment(i, j) * sig(j)
[1992]	792	bsum = bsum + ment(i, j)
	793	END DO
	794	bsum = amax1(bsum, 1.0E-16)
[5144]	795	bsum = 1.0 / bsum
[1992]	796	DO j = icb - 1, inb
[5144]	797	ment(i, j) = ment(i, j) * asum * bsum
[1992]	798	END DO
	799	csum = 0.0
	800	DO j = icb - 1, inb
	801	csum = csum + ment(i, j)
	802	END DO
	803
	804	IF (csum<m(i)) THEN
	805	nent(i) = 0
	806	ment(i, i) = m(i)
[5144]	807	qent(i, i) = rr(1) - ep(i) * clw(i)
[1992]	808	uent(i, i) = u(nk)
	809	vent(i, i) = v(nk)
	810	DO j = 1, ntra
	811	traent(i, i, j) = tra(nk, j)
	812	END DO
	813	elij(i, i) = clw(i)
	814	sij(i, i) = 1.0
[524]	815	END IF
[1992]	816	END IF
	817	END DO
[524]	818
[1992]	819
	820
	821	! **************************************************************
	822	! * CALCUL DES MENTS(I,J) ET DES QENTS(I,J)
	823	! *************************************************************
	824
	825	DO im = 1, nd
	826	DO jm = 1, nd
	827
	828	qents(im, jm) = qent(im, jm)
	829	ments(im, jm) = ment(im, jm)
	830	END DO
	831	END DO
	832
	833	! **********************************************************
	834	! --- test sb:
[5116]	835	! @ WRITE(,) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
	836	! @ WRITE(,) 'inb,m(inb),ment(inb,inb),sigij(inb,inb):'
	837	! @ WRITE(,) inb,m(inb),ment(inb,inb),sigij(inb,inb)
	838	! @ WRITE(,) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
[1992]	839	! ---
	840
	841
	842
	843
	844
	845	! * CHECK WHETHER EP(INB)=0, IF SO, SKIP PRECIPITATING *
	846	! * DOWNDRAFT CALCULATION *
	847
	848	IF (ep(inb)<0.0001) GO TO 405
	849
	850	! * INTEGRATE LIQUID WATER EQUATION TO FIND CONDENSED WATER *
	851	! * AND CONDENSED WATER FLUX *
	852
	853	wflux = 0.0
[5144]	854	tinv = 1. / 3.
[1992]	855
	856	! * BEGIN DOWNDRAFT LOOP *
	857
	858	DO i = inb, 1, -1
	859
	860	! * CALCULATE DETRAINED PRECIPITATION *
	861
[5144]	862	wdtrain = 10.0 * ep(i) * m(i) * clw(i)
[1992]	863	IF (i>1) THEN
	864	DO j = 1, i - 1
[5144]	865	awat = elij(j, i) - (1. - ep(i)) * clw(i)
[1992]	866	awat = amax1(awat, 0.0)
[5144]	867	wdtrain = wdtrain + 10.0 * awat * ment(j, i)
[1992]	868	END DO
	869	END IF
	870
	871	! * FIND RAIN WATER AND EVAPORATION USING PROVISIONAL *
	872	! * ESTIMATES OF RP(I)AND RP(I-1) *
	873
	874	wt(i) = 45.0
	875	IF (i<inb) THEN
[5144]	876	rp(i) = rp(i + 1) + (cpd * (t(i + 1) - t(i)) + gz(i + 1) - gz(i)) / lv(i)
	877	rp(i) = 0.5 * (rp(i) + rr(i))
[1992]	878	END IF
	879	rp(i) = amax1(rp(i), 0.0)
	880	rp(i) = amin1(rp(i), rs(i))
	881	rp(inb) = rr(inb)
	882	IF (i==1) THEN
[5144]	883	afac = p(1) * (rs(1) - rp(1)) / (1.0E4 + 2000.0 * p(1) * rs(1))
[1992]	884	ELSE
[5144]	885	rp(i - 1) = rp(i) + (cpd * (t(i) - t(i - 1)) + gz(i) - gz(i - 1)) / lv(i)
	886	rp(i - 1) = 0.5 * (rp(i - 1) + rr(i - 1))
	887	rp(i - 1) = amin1(rp(i - 1), rs(i - 1))
	888	rp(i - 1) = amax1(rp(i - 1), 0.0)
	889	afac1 = p(i) * (rs(i) - rp(i)) / (1.0E4 + 2000.0 * p(i) * rs(i))
	890	afac2 = p(i - 1) * (rs(i - 1) - rp(i - 1)) / (1.0E4 + 2000.0 * p(i - 1) * rs(i - 1))
	891	afac = 0.5 * (afac1 + afac2)
[1992]	892	END IF
	893	IF (i==inb) afac = 0.0
	894	afac = amax1(afac, 0.0)
[5144]	895	bfac = 1. / (sigd * wt(i))
[1992]	896
	897	! jyg1
	898	! CC SIGT=1.0
	899	! CC IF(I.GE.ICB)SIGT=SIGP(I)
	900	! Prise en compte de la variation progressive de SIGT dans
	901	! les couches ICB et ICB-1:
	902	! Pour PLCL<PH(I+1), PR1=0 & PR2=1
	903	! Pour PLCL>PH(I), PR1=1 & PR2=0
	904	! Pour PH(I+1)<PLCL<PH(I), PR1 est la proportion a cheval
	905	! sur le nuage, et PR2 est la proportion sous la base du
	906	! nuage.
[5144]	907	pr1 = (plcl - ph(i + 1)) / (ph(i) - ph(i + 1))
	908	pr1 = max(0., min(1., pr1))
	909	pr2 = (ph(i) - plcl) / (ph(i) - ph(i + 1))
	910	pr2 = max(0., min(1., pr2))
	911	sigt = sigp(i) * pr1 + pr2
[5160]	912	! sb3d PRINT *,'i,sigt,pr1,pr2', i,sigt,pr1,pr2
[1992]	913	! jyg2
	914
[5144]	915	b6 = bfac * 50. * sigd * (ph(i) - ph(i + 1)) * sigt * afac
	916	c6 = water(i + 1) + bfac * wdtrain - 50. * sigd * bfac * (ph(i) - ph(i + 1)) * evap(i + 1)
[1992]	917	IF (c6>0.0) THEN
[5144]	918	revap = 0.5 * (-b6 + sqrt(b6 * b6 + 4. * c6))
	919	evap(i) = sigt * afac * revap
	920	water(i) = revap * revap
[1992]	921	ELSE
[5144]	922	evap(i) = -evap(i + 1) + 0.02 * (wdtrain + sigd * wt(i) * water(i + 1)) / (sigd * (ph(i &
	923	) - ph(i + 1)))
[1992]	924	END IF
	925
	926
	927
	928	! * CALCULATE PRECIPITATING DOWNDRAFT MASS FLUX UNDER *
	929	! * HYDROSTATIC APPROXIMATION *
	930
	931	IF (i==1) GO TO 360
	932	tevap = amax1(0.0, evap(i))
[5144]	933	delth = amax1(0.001, (th(i) - th(i - 1)))
	934	mp(i) = 10. * lvcp(i) * sigd * tevap * (p(i - 1) - p(i)) / delth
[1992]	935
	936	! * IF HYDROSTATIC ASSUMPTION FAILS, *
	937	! * SOLVE CUBIC DIFFERENCE EQUATION FOR DOWNDRAFT THETA *
	938	! * AND MASS FLUX FROM TWO SIMULTANEOUS DIFFERENTIAL EQNS *
	939
[5144]	940	amfac = sigd * sigd * 70.0 * ph(i) * (p(i - 1) - p(i)) * (th(i) - th(i - 1)) / (tv(i) * th(i))
	941	amp2 = abs(mp(i + 1) * mp(i + 1) - mp(i) * mp(i))
	942	IF (amp2>(0.1 * amfac)) THEN
	943	xf = 100.0 * sigd * sigd * sigd * (ph(i) - ph(i + 1))
	944	tf = b(i) - 5.0 * (th(i) - th(i - 1)) * t(i) / (lvcp(i) * sigd * th(i))
	945	af = xf * tf + mp(i + 1) * mp(i + 1) * tinv
	946	bf = 2. * (tinv * mp(i + 1))*3 + tinv mp(i + 1) * xf * tf + &
	947	50. * (p(i - 1) - p(i)) * xf * tevap
[1992]	948	fac2 = 1.0
	949	IF (bf<0.0) fac2 = -1.0
	950	bf = abs(bf)
[5144]	951	ur = 0.25 * bf * bf - af * af * af * tinv * tinv * tinv
[1992]	952	IF (ur>=0.0) THEN
	953	sru = sqrt(ur)
	954	fac = 1.0
[5144]	955	IF ((0.5 * bf - sru)<0.0) fac = -1.0
	956	mp(i) = mp(i + 1) * tinv + (0.5 * bf + sru)**tinv + &
	957	fac * (abs(0.5 * bf - sru))**tinv
[524]	958	ELSE
[5144]	959	d = atan(2. * sqrt(-ur) / (bf + 1.0E-28))
[1992]	960	IF (fac2<0.0) d = 3.14159 - d
[5144]	961	mp(i) = mp(i + 1) * tinv + 2. * sqrt(af * tinv) * cos(d * tinv)
[524]	962	END IF
[1992]	963	mp(i) = amax1(0.0, mp(i))
[5144]	964	b(i - 1) = b(i) + 100.0 * (p(i - 1) - p(i)) * tevap / (mp(i) + sigd * 0.1) - &
	965	10.0 * (th(i) - th(i - 1)) * t(i) / (lvcp(i) * sigd * th(i))
	966	b(i - 1) = amax1(b(i - 1), 0.0)
[1992]	967	END IF
[524]	968
	969
	970
[1992]	971	! * LIMIT MAGNITUDE OF MP(I) TO MEET CFL CONDITION *
	972
[5144]	973	ampmax = 2.0 * (ph(i) - ph(i + 1)) * delti
	974	amp2 = 2.0 * (ph(i - 1) - ph(i)) * delti
[1992]	975	ampmax = amin1(ampmax, amp2)
	976	mp(i) = amin1(mp(i), ampmax)
	977
	978	! * FORCE MP TO DECREASE LINEARLY TO ZERO *
	979	! * BETWEEN CLOUD BASE AND THE SURFACE *
	980
	981	IF (p(i)>p(icb)) THEN
[5144]	982	mp(i) = mp(icb) * (p(1) - p(i)) / (p(1) - p(icb))
[1992]	983	END IF
[5144]	984	360 CONTINUE
[1992]	985
	986	! * FIND MIXING RATIO OF PRECIPITATING DOWNDRAFT *
	987
	988	IF (i==inb) GO TO 400
	989	rp(i) = rr(i)
[5144]	990	IF (mp(i)>mp(i + 1)) THEN
	991	rp(i) = rp(i + 1) * mp(i + 1) + rr(i) * (mp(i) - mp(i + 1)) + &
	992	5. * sigd * (ph(i) - ph(i + 1)) * (evap(i + 1) + evap(i))
	993	rp(i) = rp(i) / mp(i)
	994	up(i) = up(i + 1) * mp(i + 1) + u(i) * (mp(i) - mp(i + 1))
	995	up(i) = up(i) / mp(i)
	996	vp(i) = vp(i + 1) * mp(i + 1) + v(i) * (mp(i) - mp(i + 1))
	997	vp(i) = vp(i) / mp(i)
[1992]	998	DO j = 1, ntra
[5144]	999	trap(i, j) = trap(i + 1, j) * mp(i + 1) + trap(i, j) * (mp(i) - mp(i + 1))
	1000	trap(i, j) = trap(i, j) / mp(i)
[1992]	1001	END DO
	1002	ELSE
[5144]	1003	IF (mp(i + 1)>1.0E-16) THEN
	1004	rp(i) = rp(i + 1) + 5.0 * sigd * (ph(i) - ph(i + 1)) * (evap(i + 1) + evap(i)) / mp(i + 1 &
	1005	)
	1006	up(i) = up(i + 1)
	1007	vp(i) = vp(i + 1)
[1992]	1008	DO j = 1, ntra
[5144]	1009	trap(i, j) = trap(i + 1, j)
[524]	1010	END DO
[1992]	1011	END IF
	1012	END IF
	1013	rp(i) = amin1(rp(i), rs(i))
	1014	rp(i) = amax1(rp(i), 0.0)
[5144]	1015	400 END DO
[524]	1016
[1992]	1017	! * CALCULATE SURFACE PRECIPITATION IN MM/DAY *
[524]	1018
[5144]	1019	precip = wt(1) * sigd * water(1) * 8640.0
[524]	1020
[1992]	1021	! sb *** Calculate downdraft velocity scale and surface temperature and
	1022	! ***
	1023	! sb *** water vapor fluctuations
	1024	! ***
[5158]	1025	! sb (inspire de convect 4.3)
[524]	1026
[1992]	1027	! BETAD=10.0
	1028	betad = 5.0
[5144]	1029	wd = betad * abs(mp(icb)) * 0.01 * rd * t(icb) / (sigd * p(icb))
[524]	1030
[5144]	1031	405 CONTINUE
[524]	1032
[1992]	1033	! * CALCULATE TENDENCIES OF LOWEST LEVEL POTENTIAL TEMPERATURE *
	1034	! * AND MIXING RATIO *
[524]	1035
[5144]	1036	dpinv = 1.0 / (ph(1) - ph(2))
[1992]	1037	am = 0.0
	1038	DO k = 2, inb
	1039	am = am + m(k)
	1040	END DO
[5144]	1041	IF ((0.1 * dpinv * am)>=delti) iflag = 4
	1042	ft(1) = 0.1 * dpinv * am * (t(2) - t(1) + (gz(2) - gz(1)) / cpn(1))
	1043	ft(1) = ft(1) - 0.5 * lvcp(1) * sigd * (evap(1) + evap(2))
	1044	ft(1) = ft(1) - 0.09 * sigd * mp(2) * t(1) * b(1) * dpinv
	1045	ft(1) = ft(1) + 0.01 * sigd * wt(1) * (cl - cpd) * water(2) * (t(2) - t(1)) * dpinv / cpn(1)
	1046	fr(1) = 0.1 * mp(2) * (rp(2) - rr(1)) * & ! correction bug conservation eau
	1047	! 1 DPINV+SIGD0.5(EVAP(1)+EVAP(2))
	1048	dpinv + sigd * 0.5 * (evap(1) + evap(2))
[1992]	1049	! IM cf. SBL
	1050	! 1 DPINV+SIGD*EVAP(1)
[5144]	1051	fr(1) = fr(1) + 0.1 * am * (rr(2) - rr(1)) * dpinv
	1052	fu(1) = fu(1) + 0.1 * dpinv * (mp(2) * (up(2) - u(1)) + am * (u(2) - u(1)))
	1053	fv(1) = fv(1) + 0.1 * dpinv * (mp(2) * (vp(2) - v(1)) + am * (v(2) - v(1)))
[1992]	1054	DO j = 1, ntra
[5144]	1055	ftra(1, j) = ftra(1, j) + 0.1 * dpinv * (mp(2) * (trap(2, j) - tra(1, &
	1056	j)) + am * (tra(2, j) - tra(1, j)))
[1992]	1057	END DO
	1058	amde = 0.0
	1059	DO j = 2, inb
[5144]	1060	fr(1) = fr(1) + 0.1 * dpinv * ment(j, 1) * (qent(j, 1) - rr(1))
	1061	fu(1) = fu(1) + 0.1 * dpinv * ment(j, 1) * (uent(j, 1) - u(1))
	1062	fv(1) = fv(1) + 0.1 * dpinv * ment(j, 1) * (vent(j, 1) - v(1))
[1992]	1063	DO k = 1, ntra
[5144]	1064	ftra(1, k) = ftra(1, k) + 0.1 * dpinv * ment(j, 1) * (traent(j, 1, k) - tra(1, k))
[1992]	1065	END DO
	1066	END DO
	1067
	1068	! * CALCULATE TENDENCIES OF POTENTIAL TEMPERATURE AND MIXING RATIO *
	1069	! * AT LEVELS ABOVE THE LOWEST LEVEL *
	1070
	1071	! * FIRST FIND THE NET SATURATED UPDRAFT AND DOWNDRAFT MASS FLUXES *
	1072	! * THROUGH EACH LEVEL *
	1073
	1074	DO i = 2, inb
[5144]	1075	dpinv = 1.0 / (ph(i) - ph(i + 1))
	1076	cpinv = 1.0 / cpn(i)
[1992]	1077	amp1 = 0.0
	1078	DO k = i + 1, inb + 1
	1079	amp1 = amp1 + m(k)
	1080	END DO
	1081	DO k = 1, i
	1082	DO j = i + 1, inb + 1
	1083	amp1 = amp1 + ment(k, j)
	1084	END DO
	1085	END DO
[5144]	1086	IF ((0.1 * dpinv * amp1)>=delti) iflag = 4
[1992]	1087	ad = 0.0
	1088	DO k = 1, i - 1
	1089	DO j = i, inb
	1090	ad = ad + ment(j, k)
	1091	END DO
	1092	END DO
[5144]	1093	ft(i) = 0.1 * dpinv * (amp1 * (t(i + 1) - t(i) + (gz(i + 1) - gz(i)) * cpinv) - ad * (t(i) - t(i - &
	1094	1) + (gz(i) - gz(i - 1)) * cpinv)) - 0.5 * sigd * lvcp(i) * (evap(i) + evap(i + 1))
	1095	rat = cpn(i - 1) * cpinv
	1096	ft(i) = ft(i) - 0.09 * sigd * (mp(i + 1) * t(i) * b(i) - mp(i) * t(i - 1) * rat * b(i - 1)) * &
	1097	dpinv
	1098	ft(i) = ft(i) + 0.1 * dpinv * ment(i, i) * (hp(i) - h(i) + t(i) * (cpv - cpd) * (rr(i) - &
	1099	qent(i, i))) * cpinv
	1100	ft(i) = ft(i) + 0.01 * sigd * wt(i) * (cl - cpd) * water(i + 1) * (t(i + 1) - t(i)) * dpinv * &
	1101	cpinv
	1102	fr(i) = 0.1 * dpinv * (amp1 * (rr(i + 1) - rr(i)) - ad * (rr(i) - rr(i - 1)))
	1103	fu(i) = fu(i) + 0.1 * dpinv * (amp1 * (u(i + 1) - u(i)) - ad * (u(i) - u(i - 1)))
	1104	fv(i) = fv(i) + 0.1 * dpinv * (amp1 * (v(i + 1) - v(i)) - ad * (v(i) - v(i - 1)))
[1992]	1105	DO k = 1, ntra
[5144]	1106	ftra(i, k) = ftra(i, k) + 0.1 * dpinv * (amp1 * (tra(i + 1, k) - tra(i, &
	1107	k)) - ad * (tra(i, k) - tra(i - 1, k)))
[1992]	1108	END DO
	1109	DO k = 1, i - 1
[5144]	1110	awat = elij(k, i) - (1. - ep(i)) * clw(i)
[1992]	1111	awat = amax1(awat, 0.0)
[5144]	1112	fr(i) = fr(i) + 0.1 * dpinv * ment(k, i) * (qent(k, i) - awat - rr(i))
	1113	fu(i) = fu(i) + 0.1 * dpinv * ment(k, i) * (uent(k, i) - u(i))
	1114	fv(i) = fv(i) + 0.1 * dpinv * ment(k, i) * (vent(k, i) - v(i))
[1992]	1115	! (saturated updrafts resulting from mixing) ! cld
[5144]	1116	qcond(i) = qcond(i) + (elij(k, i) - awat) ! cld
[1992]	1117	nqcond(i) = nqcond(i) + 1. ! cld
	1118	DO j = 1, ntra
[5144]	1119	ftra(i, j) = ftra(i, j) + 0.1 * dpinv * ment(k, i) * (traent(k, i, j) - tra(i, j &
	1120	))
[1992]	1121	END DO
	1122	END DO
	1123	DO k = i, inb
[5144]	1124	fr(i) = fr(i) + 0.1 * dpinv * ment(k, i) * (qent(k, i) - rr(i))
	1125	fu(i) = fu(i) + 0.1 * dpinv * ment(k, i) * (uent(k, i) - u(i))
	1126	fv(i) = fv(i) + 0.1 * dpinv * ment(k, i) * (vent(k, i) - v(i))
[1992]	1127	DO j = 1, ntra
[5144]	1128	ftra(i, j) = ftra(i, j) + 0.1 * dpinv * ment(k, i) * (traent(k, i, j) - tra(i, j &
	1129	))
[1992]	1130	END DO
	1131	END DO
[5144]	1132	fr(i) = fr(i) + 0.5 * sigd * (evap(i) + evap(i + 1)) + 0.1 * (mp(i + 1) * (rp(i + &
	1133	1) - rr(i)) - mp(i) * (rp(i) - rr(i - 1))) * dpinv
	1134	fu(i) = fu(i) + 0.1 * (mp(i + 1) * (up(i + 1) - u(i)) - mp(i) * (up(i) - u(i - 1))) * dpinv
	1135	fv(i) = fv(i) + 0.1 * (mp(i + 1) * (vp(i + 1) - v(i)) - mp(i) * (vp(i) - v(i - 1))) * dpinv
[1992]	1136	DO j = 1, ntra
[5144]	1137	ftra(i, j) = ftra(i, j) + 0.1 * dpinv * (mp(i + 1) * (trap(i + 1, j) - tra(i, &
	1138	j)) - mp(i) * (trap(i, j) - trap(i - 1, j)))
[1992]	1139	END DO
	1140	! (saturated downdrafts resulting from mixing) ! cld
	1141	DO k = i + 1, inb ! cld
	1142	qcond(i) = qcond(i) + elij(k, i) ! cld
	1143	nqcond(i) = nqcond(i) + 1. ! cld
	1144	END DO ! cld
	1145	! (particular case: no detraining level is found) ! cld
	1146	IF (nent(i)==0) THEN ! cld
[5144]	1147	qcond(i) = qcond(i) + (1 - ep(i)) * clw(i) ! cld
[1992]	1148	nqcond(i) = nqcond(i) + 1. ! cld
	1149	END IF ! cld
	1150	IF (nqcond(i)/=0.) THEN ! cld
[5144]	1151	qcond(i) = qcond(i) / nqcond(i) ! cld
[1992]	1152	END IF ! cld
	1153	END DO
	1154
	1155
	1156
	1157
	1158	! * MOVE THE DETRAINMENT AT LEVEL INB DOWN TO LEVEL INB-1 *
	1159	! * IN SUCH A WAY AS TO PRESERVE THE VERTICALLY *
	1160	! * INTEGRATED ENTHALPY AND WATER TENDENCIES *
	1161
	1162	! test sb:
[5116]	1163	! @ WRITE(,) '--------------------------------------------'
	1164	! @ WRITE(,) 'inb,ft,hp,h,t,rr,qent,ment,water,waterp,wt,mp,b'
	1165	! @ WRITE(,) inb,ft(inb),hp(inb),h(inb)
[1992]	1166	! @ : ,t(inb),rr(inb),qent(inb,inb)
	1167	! @ : ,ment(inb,inb),water(inb)
	1168	! @ : ,water(inb+1),wt(inb),mp(inb),b(inb)
[5116]	1169	! @ WRITE(,) '--------------------------------------------'
[1992]	1170	! fin test sb:
	1171
[5144]	1172	ax = 0.1 * ment(inb, inb) * (hp(inb) - h(inb) + t(inb) * (cpv - cpd) * (rr(inb) - qent(inb, &
	1173	inb))) / (cpn(inb) * (ph(inb) - ph(inb + 1)))
[1992]	1174	ft(inb) = ft(inb) - ax
[5144]	1175	ft(inb - 1) = ft(inb - 1) + ax * cpn(inb) * (ph(inb) - ph(inb + 1)) / (cpn(inb - 1) * (ph(inb &
	1176	- 1) - ph(inb)))
	1177	bx = 0.1 * ment(inb, inb) * (qent(inb, inb) - rr(inb)) / (ph(inb) - ph(inb + 1))
[1992]	1178	fr(inb) = fr(inb) - bx
[5144]	1179	fr(inb - 1) = fr(inb - 1) + bx * (ph(inb) - ph(inb + 1)) / (ph(inb - 1) - ph(inb))
	1180	cx = 0.1 * ment(inb, inb) * (uent(inb, inb) - u(inb)) / (ph(inb) - ph(inb + 1))
[1992]	1181	fu(inb) = fu(inb) - cx
[5144]	1182	fu(inb - 1) = fu(inb - 1) + cx * (ph(inb) - ph(inb + 1)) / (ph(inb - 1) - ph(inb))
	1183	dx = 0.1 * ment(inb, inb) * (vent(inb, inb) - v(inb)) / (ph(inb) - ph(inb + 1))
[1992]	1184	fv(inb) = fv(inb) - dx
[5144]	1185	fv(inb - 1) = fv(inb - 1) + dx * (ph(inb) - ph(inb + 1)) / (ph(inb - 1) - ph(inb))
[1992]	1186	DO j = 1, ntra
[5144]	1187	ex = 0.1 * ment(inb, inb) * (traent(inb, inb, j) - tra(inb, j)) / &
	1188	(ph(inb) - ph(inb + 1))
[1992]	1189	ftra(inb, j) = ftra(inb, j) - ex
[5144]	1190	ftra(inb - 1, j) = ftra(inb - 1, j) + ex * (ph(inb) - ph(inb + 1)) / (ph(inb - 1) - ph(&
	1191	inb))
[1992]	1192	END DO
	1193
	1194	! * HOMOGINIZE TENDENCIES BELOW CLOUD BASE *
	1195
	1196	asum = 0.0
	1197	bsum = 0.0
	1198	csum = 0.0
	1199	dsum = 0.0
	1200	DO i = 1, icb - 1
[5144]	1201	asum = asum + ft(i) * (ph(i) - ph(i + 1))
	1202	bsum = bsum + fr(i) * (lv(i) + (cl - cpd) * (t(i) - t(1))) * (ph(i) - ph(i + 1))
	1203	csum = csum + (lv(i) + (cl - cpd) * (t(i) - t(1))) * (ph(i) - ph(i + 1))
	1204	dsum = dsum + t(i) * (ph(i) - ph(i + 1)) / th(i)
[1992]	1205	END DO
	1206	DO i = 1, icb - 1
[5144]	1207	ft(i) = asum * t(i) / (th(i) * dsum)
	1208	fr(i) = bsum / csum
[1992]	1209	END DO
	1210
	1211	! * RESET COUNTER AND RETURN *
	1212
	1213	sig(nd) = 2.0
	1214
	1215	DO i = 1, nd
	1216	upwd(i) = 0.0
	1217	dnwd(i) = 0.0
	1218	! sb dnwd0(i) = - mp(i)
	1219	END DO
	1220
	1221	DO i = 1, nl
	1222	dnwd0(i) = -mp(i)
	1223	END DO
	1224	DO i = nl + 1, nd
	1225	dnwd0(i) = 0.
	1226	END DO
	1227
	1228	DO i = icb, inb
	1229	upwd(i) = 0.0
	1230	dnwd(i) = 0.0
	1231
	1232	DO k = i, inb
	1233	up1 = 0.0
	1234	dn1 = 0.0
	1235	DO n = 1, i - 1
	1236	up1 = up1 + ment(n, k)
	1237	dn1 = dn1 - ment(k, n)
	1238	END DO
	1239	upwd(i) = upwd(i) + m(k) + up1
	1240	dnwd(i) = dnwd(i) + dn1
	1241	END DO
	1242	END DO
	1243
	1244	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1245	! DETERMINATION DE LA VARIATION DE FLUX ASCENDANT ENTRE
	1246	! DEUX NIVEAU NON DILUE Mike
	1247	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1248
	1249
	1250	! sb do i=1,ND
	1251	! sb Mike(i)=M(i)
	1252	! sb enddo
	1253
	1254	DO i = 1, nl
	1255	mike(i) = m(i)
	1256	END DO
	1257	DO i = nl + 1, nd
	1258	mike(i) = 0.
	1259	END DO
	1260
	1261	DO i = 1, nd
	1262	ma(i) = 0
	1263	END DO
	1264
	1265	! sb do i=1,nd
	1266	! sb do j=i,nd
	1267	! sb Ma(i)=Ma(i)+M(j)
	1268	! sb enddo
	1269	! sb enddo
	1270
	1271	DO i = 1, nl
	1272	DO j = i, nl
	1273	ma(i) = ma(i) + m(j)
	1274	END DO
	1275	END DO
	1276
	1277	DO i = nl + 1, nd
	1278	ma(i) = 0.
	1279	END DO
	1280
	1281	DO i = 1, icb - 1
	1282	ma(i) = 0
	1283	END DO
	1284
	1285
	1286
	1287	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1288	! ICB REPRESENTE DE NIVEAU OU SE TROUVE LA
	1289	! BASE DU NUAGE , ET INB LE TOP DU NUAGE
	1290	! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1291
	1292	DO i = 1, nd
	1293	mke(i) = upwd(i) + dnwd(i)
	1294	END DO
	1295
	1296
	1297	! * Diagnose the in-cloud mixing ratio * ! cld
	1298	! * of condensed water * ! cld
[5113]	1299	! cld
[1992]	1300	DO i = 1, nd ! cld
	1301	maa(i) = 0.0 ! cld
	1302	wa(i) = 0.0 ! cld
	1303	siga(i) = 0.0 ! cld
	1304	END DO ! cld
	1305	DO i = nk, inb ! cld
	1306	DO k = i + 1, inb + 1 ! cld
	1307	maa(i) = maa(i) + m(k) ! cld
	1308	END DO ! cld
	1309	END DO ! cld
	1310	DO i = icb, inb - 1 ! cld
	1311	axc(i) = 0. ! cld
	1312	DO j = icb, i ! cld
[5144]	1313	axc(i) = axc(i) + rd * (tvp(j) - tv(j)) * (ph(j) - ph(j + 1)) / p(j) ! cld
[1992]	1314	END DO ! cld
	1315	IF (axc(i)>0.0) THEN ! cld
[5144]	1316	wa(i) = sqrt(2. * axc(i)) ! cld
[1992]	1317	END IF ! cld
	1318	END DO ! cld
	1319	DO i = 1, nl ! cld
	1320	IF (wa(i)>0.0) & ! cld
[5144]	1321	siga(i) = maa(i) / wa(i) * rd * tvp(i) / p(i) / 100. / deltac ! cld
[1992]	1322	siga(i) = min(siga(i), 1.0) ! cld
[5144]	1323	qcondc(i) = siga(i) * clw(i) * (1. - ep(i)) & ! cld
	1324	+ (1. - siga(i)) * qcond(i) ! cld
[1992]	1325	END DO ! cld
	1326
	1327
	1328	! @$$cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
[5101]	1329	! @$$ CALL writeg1d(1,klev,ma,'ma ','ma ')
	1330	! @$$ CALL writeg1d(1,klev,upwd,'upwd ','upwd ')
	1331	! @$$ CALL writeg1d(1,klev,dnwd,'dnwd ','dnwd ')
	1332	! @$$ CALL writeg1d(1,klev,dnwd0,'dnwd0 ','dnwd0 ')
	1333	! @$$ CALL writeg1d(1,klev,tvp,'tvp ','tvp ')
	1334	! @$$ CALL writeg1d(1,klev,tra(1:klev,3),'tra3 ','tra3 ')
	1335	! @$$ CALL writeg1d(1,klev,tra(1:klev,4),'tra4 ','tra4 ')
	1336	! @$$ CALL writeg1d(1,klev,tra(1:klev,5),'tra5 ','tra5 ')
	1337	! @$$ CALL writeg1d(1,klev,tra(1:klev,6),'tra6 ','tra6 ')
	1338	! @$$ CALL writeg1d(1,klev,tra(1:klev,7),'tra7 ','tra7 ')
	1339	! @$$ CALL writeg1d(1,klev,tra(1:klev,8),'tra8 ','tra8 ')
	1340	! @$$ CALL writeg1d(1,klev,tra(1:klev,9),'tra9 ','tra9 ')
	1341	! @$$ CALL writeg1d(1,klev,tra(1:klev,10),'tra10','tra10')
	1342	! @$$ CALL writeg1d(1,klev,tra(1:klev,11),'tra11','tra11')
	1343	! @$$ CALL writeg1d(1,klev,tra(1:klev,12),'tra12','tra12')
	1344	! @$$ CALL writeg1d(1,klev,tra(1:klev,13),'tra13','tra13')
	1345	! @$$ CALL writeg1d(1,klev,tra(1:klev,14),'tra14','tra14')
	1346	! @$$ CALL writeg1d(1,klev,tra(1:klev,15),'tra15','tra15')
	1347	! @$$ CALL writeg1d(1,klev,tra(1:klev,16),'tra16','tra16')
	1348	! @$$ CALL writeg1d(1,klev,tra(1:klev,17),'tra17','tra17')
	1349	! @$$ CALL writeg1d(1,klev,tra(1:klev,18),'tra18','tra18')
	1350	! @$$ CALL writeg1d(1,klev,tra(1:klev,19),'tra19','tra19')
	1351	! @$$ CALL writeg1d(1,klev,tra(1:klev,20),'tra20','tra20 ')
	1352	! @$$ CALL writeg1d(1,klev,trap(1:klev,1),'trp1','trp1')
	1353	! @$$ CALL writeg1d(1,klev,trap(1:klev,2),'trp2','trp2')
	1354	! @$$ CALL writeg1d(1,klev,trap(1:klev,3),'trp3','trp3')
	1355	! @$$ CALL writeg1d(1,klev,trap(1:klev,4),'trp4','trp4')
	1356	! @$$ CALL writeg1d(1,klev,trap(1:klev,5),'trp5','trp5')
	1357	! @$$ CALL writeg1d(1,klev,trap(1:klev,10),'trp10','trp10')
	1358	! @$$ CALL writeg1d(1,klev,trap(1:klev,12),'trp12','trp12')
	1359	! @$$ CALL writeg1d(1,klev,trap(1:klev,15),'trp15','trp15')
	1360	! @$$ CALL writeg1d(1,klev,trap(1:klev,20),'trp20','trp20')
	1361	! @$$ CALL writeg1d(1,klev,ftra(1:klev,1),'ftr1 ','ftr1 ')
	1362	! @$$ CALL writeg1d(1,klev,ftra(1:klev,2),'ftr2 ','ftr2 ')
	1363	! @$$ CALL writeg1d(1,klev,ftra(1:klev,3),'ftr3 ','ftr3 ')
	1364	! @$$ CALL writeg1d(1,klev,ftra(1:klev,4),'ftr4 ','ftr4 ')
	1365	! @$$ CALL writeg1d(1,klev,ftra(1:klev,5),'ftr5 ','ftr5 ')
	1366	! @$$ CALL writeg1d(1,klev,ftra(1:klev,6),'ftr6 ','ftr6 ')
	1367	! @$$ CALL writeg1d(1,klev,ftra(1:klev,7),'ftr7 ','ftr7 ')
	1368	! @$$ CALL writeg1d(1,klev,ftra(1:klev,8),'ftr8 ','ftr8 ')
	1369	! @$$ CALL writeg1d(1,klev,ftra(1:klev,9),'ftr9 ','ftr9 ')
	1370	! @$$ CALL writeg1d(1,klev,ftra(1:klev,10),'ftr10','ftr10')
	1371	! @$$ CALL writeg1d(1,klev,ftra(1:klev,11),'ftr11','ftr11')
	1372	! @$$ CALL writeg1d(1,klev,ftra(1:klev,12),'ftr12','ftr12')
	1373	! @$$ CALL writeg1d(1,klev,ftra(1:klev,13),'ftr13','ftr13')
	1374	! @$$ CALL writeg1d(1,klev,ftra(1:klev,14),'ftr14','ftr14')
	1375	! @$$ CALL writeg1d(1,klev,ftra(1:klev,15),'ftr15','ftr15')
	1376	! @$$ CALL writeg1d(1,klev,ftra(1:klev,16),'ftr16','ftr16')
	1377	! @$$ CALL writeg1d(1,klev,ftra(1:klev,17),'ftr17','ftr17')
	1378	! @$$ CALL writeg1d(1,klev,ftra(1:klev,18),'ftr18','ftr18')
	1379	! @$$ CALL writeg1d(1,klev,ftra(1:klev,19),'ftr19','ftr19')
	1380	! @$$ CALL writeg1d(1,klev,ftra(1:klev,20),'ftr20','ftr20 ')
	1381	! @$$ CALL writeg1d(1,klev,mp,'mp ','mp ')
	1382	! @$$ CALL writeg1d(1,klev,Mke,'Mke ','Mke ')
[1992]	1383
	1384
	1385
	1386	! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1387
	1388	END SUBROUTINE convect3
	1389	! ---------------------------------------------------------------------------

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

Download in other formats:

Original Format