Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

convect43c.f90 @ 13

Last change on this file since 13 was 13, checked in by saeck, 11 years ago
update to wetdepo.f90
File size: 36.4 KB

Rev	Line
[13]	1	!**********************************************************************
	2	! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
	3	! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
	4	! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
	5	! *
	6	! This file is part of FLEXPART. *
	7	! *
	8	! FLEXPART is free software: you can redistribute it and/or modify *
	9	! it under the terms of the GNU General Public License as published by*
	10	! the Free Software Foundation, either version 3 of the License, or *
	11	! (at your option) any later version. *
	12	! *
	13	! FLEXPART is distributed in the hope that it will be useful, *
	14	! but WITHOUT ANY WARRANTY; without even the implied warranty of *
	15	! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
	16	! GNU General Public License for more details. *
	17	! *
	18	! You should have received a copy of the GNU General Public License *
	19	! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
	20	!**********************************************************************
	21
	22	!**************************************************************************
	23	!** SUBROUTINE CONVECT ***
	24	!** VERSION 4.3c ***
	25	!** 20 May, 2002 ***
	26	!** Kerry Emanuel ***
	27	!**************************************************************************
	28	!
	29	SUBROUTINE CONVECT &
	30	(ND, NL, DELT, IFLAG, &
	31	PRECIP, WD, TPRIME, QPRIME, CBMF )
	32	!
	33	!-cv *************************************************************************
	34	!-cv C. Forster, November 2003 - May 2004:
	35	!-cv
	36	!-cv The subroutine has been downloaded from Kerry Emanuel's homepage,
	37	!-cv where further infos on the convection scheme can be found
	38	!-cv http://www-paoc.mit.edu/~emanuel/home.html
	39	!-cv
	40	!-cv The following changes have been made to integrate this subroutine
	41	!-cv into FLEXPART
	42	!-cv
	43	!-cv Putting most of the variables in a new common block
	44	!-cv renaming eps to eps0 because there is some eps already in includepar
	45	!-cv
	46	!-cv removing the arrays U,V,TRA and related arrays
	47	!-cv
	48	!-cv renaming the original arrays T,Q,QS,P,PH to
	49	!-cv TCONV,QCONV,QSCONV,PCONV_HPA,PHCONV_HPA
	50	!-cv
	51	!-cv Initialization of variables has been put into parameter statements
	52	!-cv instead of assignment of values at each call, in order to save
	53	!-cv computation time.
	54	!***************************************************************************
	55	!
	56	!-----------------------------------------------------------------------------
	57	! * On input: *
	58	!
	59	!T: Array of absolute temperature (K) of dimension ND, with first
	60	! index corresponding to lowest model level. Note that this array
	61	! will be altered by the subroutine if dry convective adjustment
	62	! occurs and if IPBL is not equal to 0.
	63	!
	64	!Q: Array of specific humidity (gm/gm) of dimension ND, with first
	65	! index corresponding to lowest model level. Must be defined
	66	! at same grid levels as T. Note that this array will be altered
	67	! if dry convective adjustment occurs and if IPBL is not equal to 0.
	68	!
	69	!QS: Array of saturation specific humidity of dimension ND, with first
	70	! index corresponding to lowest model level. Must be defined
	71	! at same grid levels as T. Note that this array will be altered
	72	! if dry convective adjustment occurs and if IPBL is not equal to 0.
	73	!
	74	!U: Array of zonal wind velocity (m/s) of dimension ND, witth first
	75	! index corresponding with the lowest model level. Defined at
	76	! same levels as T. Note that this array will be altered if
	77	! dry convective adjustment occurs and if IPBL is not equal to 0.
	78	!
	79	!V: Same as U but for meridional velocity.
	80	!
	81	!TRA: Array of passive tracer mixing ratio, of dimensions (ND,NTRA),
	82	! where NTRA is the number of different tracers. If no
	83	! convective tracer transport is needed, define a dummy
	84	! input array of dimension (ND,1). Tracers are defined at
	85	! same vertical levels as T. Note that this array will be altered
	86	! if dry convective adjustment occurs and if IPBL is not equal to 0.
	87	!
	88	!P: Array of pressure (mb) of dimension ND, with first
	89	! index corresponding to lowest model level. Must be defined
	90	! at same grid levels as T.
	91	!
	92	!PH: Array of pressure (mb) of dimension ND+1, with first index
	93	! corresponding to lowest level. These pressures are defined at
	94	! levels intermediate between those of P, T, Q and QS. The first
	95	! value of PH should be greater than (i.e. at a lower level than)
	96	! the first value of the array P.
	97	!
	98	!ND: The dimension of the arrays T,Q,QS,P,PH,FT and FQ
	99	!
	100	!NL: The maximum number of levels to which convection can
	101	! penetrate, plus 1.
	102	! NL MUST be less than or equal to ND-1.
	103	!
	104	!NTRA:The number of different tracers. If no tracer transport
	105	! is needed, set this equal to 1. (On most compilers, setting
	106	! NTRA to 0 will bypass tracer calculation, saving some CPU.)
	107	!
	108	!DELT: The model time step (sec) between calls to CONVECT
	109	!
	110	!----------------------------------------------------------------------------
	111	! * On Output: *
	112	!
	113	!IFLAG: An output integer whose value denotes the following:
	114	!
	115	! VALUE INTERPRETATION
	116	! ----- --------------
	117	! 0 No moist convection; atmosphere is not
	118	! unstable, or surface temperature is less
	119	! than 250 K or surface specific humidity
	120	! is non-positive.
	121	!
	122	! 1 Moist convection occurs.
	123	!
	124	! 2 No moist convection: lifted condensation
	125	! level is above the 200 mb level.
	126	!
	127	! 3 No moist convection: cloud base is higher
	128	! then the level NL-1.
	129	!
	130	! 4 Moist convection occurs, but a CFL condition
	131	! on the subsidence warming is violated. This
	132	! does not cause the scheme to terminate.
	133	!
	134	!FT: Array of temperature tendency (K/s) of dimension ND, defined at same
	135	! grid levels as T, Q, QS and P.
	136	!
	137	!FQ: Array of specific humidity tendencies ((gm/gm)/s) of dimension ND,
	138	! defined at same grid levels as T, Q, QS and P.
	139	!
	140	!FU: Array of forcing of zonal velocity (m/s^2) of dimension ND,
	141	! defined at same grid levels as T.
	142	!
	143	!FV: Same as FU, but for forcing of meridional velocity.
	144	!
	145	!FTRA: Array of forcing of tracer content, in tracer mixing ratio per
	146	! second, defined at same levels as T. Dimensioned (ND,NTRA).
	147	!
	148	!PRECIP: Scalar convective precipitation rate (mm/day).
	149	!
	150	!WD: A convective downdraft velocity scale. For use in surface
	151	! flux parameterizations. See convect.ps file for details.
	152	!
	153	!TPRIME: A convective downdraft temperature perturbation scale (K).
	154	! For use in surface flux parameterizations. See convect.ps
	155	! file for details.
	156	!
	157	!QPRIME: A convective downdraft specific humidity
	158	! perturbation scale (gm/gm).
	159	! For use in surface flux parameterizations. See convect.ps
	160	! file for details.
	161	!
	162	!CBMF: The cloud base mass flux ((kg/m**2)/s). THIS SCALAR VALUE MUST
	163	! BE STORED BY THE CALLING PROGRAM AND RETURNED TO CONVECT AT
	164	! ITS NEXT CALL. That is, the value of CBMF must be "remembered"
	165	! by the calling program between calls to CONVECT.
	166	!
	167	!-----------------------------------------------------------------------------
	168	!
	169	! * THE PARAMETER NA SHOULD IN GENERAL BE GREATER THAN *
	170	! * OR EQUAL TO ND + 1 *
	171	!
	172	!
	173	use par_mod
	174	use conv_mod
	175
	176	implicit none
	177	!
	178	!-cv====>Begin Module CONVECT File convect.f Undeclared variables
	179	!
	180	!Argument variables
	181	!
	182	integer :: iflag, nd, nl
	183	!
	184	real :: cbmf, delt, precip, qprime, tprime, wd
	185	!
	186	!Local variables
	187	!
	188	integer :: i, icb, ihmin, inb, inb1, j, jtt, k
	189	integer :: nk
	190	!
	191	real :: ad, afac, ahmax, ahmin, alt, altem
	192	real :: am, amp1, anum, asij, awat, b6, bf2, bsum, by
	193	real :: byp, c6, cape, capem, cbmfold, chi, coeff
	194	real :: cpinv, cwat, damps, dbo, dbosum
	195	real :: defrac, dei, delm, delp, delt0, delti, denom, dhdp
	196	real :: dpinv, dtma, dtmin, dtpbl, elacrit, ents
	197	real :: epmax, fac, fqold, frac, ftold
	198	real :: plcl, qp1, qsm, qstm, qti, rat
	199	real :: rdcp, revap, rh, scrit, sigt, sjmax
	200	real :: sjmin, smid, smin, stemp, tca
	201	real :: tvaplcl, tvpplcl, tvx, tvy, wdtrain
	202
	203	!integer jc,jn
	204	!real alvnew,a2,ahm,alv,rm,sum,qnew,dphinv,tc,thbar,tnew,x
	205
	206	real :: FUP(NA),FDOWN(NA)
	207	!
	208	!-cv====>End Module CONVECT File convect.f
	209
	210	INTEGER :: NENT(NA)
	211	REAL :: M(NA),MP(NA),MENT(NA,NA),QENT(NA,NA),ELIJ(NA,NA)
	212	REAL :: SIJ(NA,NA),TVP(NA),TV(NA),WATER(NA)
	213	REAL :: QP(NA),EP(NA),TH(NA),WT(NA),EVAP(NA),CLW(NA)
	214	REAL :: SIGP(NA),TP(NA),CPN(NA)
	215	REAL :: LV(NA),LVCP(NA),H(NA),HP(NA),GZ(NA),HM(NA)
	216	!REAL TOLD(NA)
	217	!
	218	! -----------------------------------------------------------------------
	219	!
	220	! * Specify Switches *
	221	!
	222	! * IPBL: Set to zero to bypass dry adiabatic adjustment *
	223	! * Any other value results in dry adiabatic adjustment *
	224	! * (Zero value recommended for use in models with *
	225	! * boundary layer schemes) *
	226	!
	227	! * MINORIG: Lowest level from which convection may originate *
	228	! * (Should be first model level at which T is defined *
	229	! * for models using bulk PBL schemes; otherwise, it should *
	230	! * be the first model level at which T is defined above *
	231	! * the surface layer) *
	232	!
	233	INTEGER,PARAMETER :: IPBL=0
	234	INTEGER,PARAMETER :: MINORIG=1
	235	!
	236	!------------------------------------------------------------------------------
	237	!
	238	! * SPECIFY PARAMETERS *
	239	!
	240	! * ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) *
	241	! * TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO- *
	242	! * CONVERSION THRESHOLD IS ASSUMED TO BE ZERO *
	243	! * (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY *
	244	! * BETWEEN 0 C AND TLCRIT) *
	245	! * ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT *
	246	! * FORMULATION *
	247	! * SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *
	248	! * SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *
	249	! * OF CLOUD *
	250	! * OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN *
	251	! * OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW *
	252	! * COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *
	253	! * OF RAIN *
	254	! * COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *
	255	! * OF SNOW *
	256	! * CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM *
	257	! * TRANSPORT *
	258	! * DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION *
	259	! * A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC *
	260	! * ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF *
	261	! * APPROACH TO QUASI-EQUILIBRIUM *
	262	! * (THEIR STANDARD VALUES ARE 0.20 AND 0.1, RESPECTIVELY) *
	263	! * (DAMP MUST BE LESS THAN 1) *
	264	!
	265	REAL,PARAMETER :: ELCRIT=.0011
	266	REAL,PARAMETER :: TLCRIT=-55.0
	267	REAL,PARAMETER :: ENTP=1.5
	268	REAL,PARAMETER :: SIGD=0.05
	269	REAL,PARAMETER :: SIGS=0.12
	270	REAL,PARAMETER :: OMTRAIN=50.0
	271	REAL,PARAMETER :: OMTSNOW=5.5
	272	REAL,PARAMETER :: COEFFR=1.0
	273	REAL,PARAMETER :: COEFFS=0.8
	274	REAL,PARAMETER :: CU=0.7
	275	REAL,PARAMETER :: BETA=10.0
	276	REAL,PARAMETER :: DTMAX=0.9
	277	REAL,PARAMETER :: ALPHA=0.025 !original 0.2
	278	REAL,PARAMETER :: DAMP=0.1
	279	!
	280	! * ASSIGN VALUES OF THERMODYNAMIC CONSTANTS, *
	281	! * GRAVITY, AND LIQUID WATER DENSITY. *
	282	! * THESE SHOULD BE CONSISTENT WITH *
	283	! * THOSE USED IN CALLING PROGRAM *
	284	! * NOTE: THESE ARE ALSO SPECIFIED IN SUBROUTINE TLIFT *
	285	!
	286	REAL,PARAMETER :: CPD=1005.7
	287	REAL,PARAMETER :: CPV=1870.0
	288	REAL,PARAMETER :: CL=2500.0
	289	REAL,PARAMETER :: RV=461.5
	290	REAL,PARAMETER :: RD=287.04
	291	REAL,PARAMETER :: LV0=2.501E6
	292	REAL,PARAMETER :: G=9.81
	293	REAL,PARAMETER :: ROWL=1000.0
	294	!
	295	REAL,PARAMETER :: CPVMCL=CL-CPV
	296	REAL,PARAMETER :: EPS0=RD/RV
	297	REAL,PARAMETER :: EPSI=1./EPS0
	298	REAL,PARAMETER :: GINV=1.0/G
	299	REAL,PARAMETER :: EPSILON=1.e-20
	300
	301	! EPSILON IS A SMALL NUMBER USED TO EXCLUDE MASS FLUXES OF ZERO
	302	!
	303	DELTI=1.0/DELT
	304	!
	305	! * INITIALIZE OUTPUT ARRAYS AND PARAMETERS *
	306	!
	307
	308	DO I=1,NL+1
	309	FT(I)=0.0
	310	FQ(I)=0.0
	311	FDOWN(I)=0.0
	312	SUB(I)=0.0
	313	FUP(I)=0.0
	314	M(I)=0.0
	315	MP(I)=0.0
	316	DO J=1,NL+1
	317	FMASS(I,J)=0.0
	318	MENT(I,J)=0.0
	319	END DO
	320	END DO
	321	DO I=1,NL+1
	322	RDCP=(RD(1.-QCONV(I))+QCONV(I)RV)/ &
	323	(CPD(1.-QCONV(I))+QCONV(I)CPV)
	324	TH(I)=TCONV(I)(1000.0/PCONV_HPA(I))*RDCP
	325	END DO
	326	PRECIP=0.0
	327	WD=0.0
	328	TPRIME=0.0
	329	QPRIME=0.0
	330	IFLAG=0
	331	!
	332	! IF(IPBL.NE.0)THEN
	333	!
	334	!* PERFORM DRY ADIABATIC ADJUSTMENT *
	335	!
	336	! JC=0
	337	! DO 30 I=NL-1,1,-1
	338	! JN=0
	339	! SUM=TH(I)(1.+QCONV(I)EPSI-QCONV(I))
	340	! DO 10 J=I+1,NL
	341	! SUM=SUM+TH(J)(1.+QCONV(J)EPSI-QCONV(J))
	342	! THBAR=SUM/REAL(J+1-I)
	343	! IF((TH(J)(1.+QCONV(J)EPSI-QCONV(J))).LT.THBAR)JN=J
	344	! 10 CONTINUE
	345	! IF(I.EQ.1)JN=MAX(JN,2)
	346	! IF(JN.EQ.0)GOTO 30
	347	! 12 CONTINUE
	348	! AHM=0.0
	349	! RM=0.0
	350	! DO 15 J=I,JN
	351	! AHM=AHM+(CPD(1.-QCONV(J))+QCONV(J)CPV)TCONV(J)
	352	! + (PHCONV_HPA(J)-PHCONV_HPA(J+1))
	353	! RM=RM+QCONV(J)*(PHCONV_HPA(J)-PHCONV_HPA(J+1))
	354	! 15 CONTINUE
	355	! DPHINV=1./(PHCONV_HPA(I)-PHCONV_HPA(JN+1))
	356	! RM=RM*DPHINV
	357	! A2=0.0
	358	! DO 20 J=I,JN
	359	! QCONV(J)=RM
	360	! RDCP=(RD(1.-QCONV(J))+QCONV(J)RV)/
	361	! 1 (CPD(1.-QCONV(J))+QCONV(J)CPV)
	362	! X=(0.001PCONV_HPA(J))*RDCP
	363	! TOLD(J)=TCONV(J)
	364	! TCONV(J)=X
	365	! A2=A2+(CPD(1.-QCONV(J))+QCONV(J)CPV)X
	366	! 1 (PHCONV_HPA(J)-PHCONV_HPA(J+1))
	367	! 20 CONTINUE
	368	! DO 25 J=I,JN
	369	! TH(J)=AHM/A2
	370	! TCONV(J)=TCONV(J)*TH(J)
	371	! TC=TOLD(J)-273.15
	372	! ALV=LV0-CPVMCL*TC
	373	! QSCONV(J)=QSCONV(J)+QSCONV(J)(1.+QSCONV(J)(EPSI-1.))ALV
	374	! 1 (TCONV(J)- TOLD(J))/(RVTOLD(J)TOLD(J))
	375	! if (qslev(j) .lt. 0.) then
	376	! write(,) 'qslev.lt.0 ',j,qslev
	377	! endif
	378	! 25 CONTINUE
	379	! IF((TH(JN+1)(1.+QCONV(JN+1)EPSI-QCONV(JN+1))).LT.
	380	! 1 (TH(JN)(1.+QCONV(JN)EPSI-QCONV(JN))))THEN
	381	! JN=JN+1
	382	! GOTO 12
	383	! END IF
	384	! IF(I.EQ.1)JC=JN
	385	! 30 CONTINUE
	386	!
	387	! * Remove any supersaturation that results from adjustment *
	388	!
	389	!IF(JC.GT.1)THEN
	390	! DO 38 J=1,JC
	391	! IF(QSCONV(J).LT.QCONV(J))THEN
	392	! ALV=LV0-CPVMCL*(TCONV(J)-273.15)
	393	! TNEW=TCONV(J)+ALV(QCONV(J)-QSCONV(J))/(CPD(1.-QCONV(J))+
	394	! 1 CLQCONV(J)+QSCONV(J)(CPV-CL+ALVALV/(RVTCONV(J)*TCONV(J))))
	395	! ALVNEW=LV0-CPVMCL*(TNEW-273.15)
	396	! QNEW=(ALVQCONV(J)-(TNEW-TCONV(J))(CPD*(1.-QCONV(J))
	397	! 1 +CL*QCONV(J)))/ALVNEW
	398	! PRECIP=PRECIP+24.3600.1.0E5(PHCONV_HPA(J)-PHCONV_HPA(J+1))
	399	! 1 (QCONV(J)-QNEW)/(GDELTROWL)
	400	! TCONV(J)=TNEW
	401	! QCONV(J)=QNEW
	402	! QSCONV(J)=QNEW
	403	! END IF
	404	! 38 CONTINUE
	405	!END IF
	406	!
	407	!END IF
	408	!
	409	! *** CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY AND STATIC ENERGY
	410	!
	411	GZ(1)=0.0
	412	CPN(1)=CPD(1.-QCONV(1))+QCONV(1)CPV
	413	H(1)=TCONV(1)*CPN(1)
	414	LV(1)=LV0-CPVMCL*(TCONV(1)-273.15)
	415	HM(1)=LV(1)*QCONV(1)
	416	TV(1)=TCONV(1)(1.+QCONV(1)EPSI-QCONV(1))
	417	AHMIN=1.0E12
	418	IHMIN=NL
	419	DO I=2,NL+1
	420	TVX=TCONV(I)(1.+QCONV(I)EPSI-QCONV(I))
	421	TVY=TCONV(I-1)(1.+QCONV(I-1)EPSI-QCONV(I-1))
	422	GZ(I)=GZ(I-1)+0.5RD(TVX+TVY)*(PCONV_HPA(I-1)-PCONV_HPA(I))/ &
	423	PHCONV_HPA(I)
	424	CPN(I)=CPD(1.-QCONV(I))+CPVQCONV(I)
	425	H(I)=TCONV(I)*CPN(I)+GZ(I)
	426	LV(I)=LV0-CPVMCL*(TCONV(I)-273.15)
	427	HM(I)=(CPD(1.-QCONV(I))+CLQCONV(I))*(TCONV(I)-TCONV(1))+ &
	428	LV(I)*QCONV(I)+GZ(I)
	429	TV(I)=TCONV(I)(1.+QCONV(I)EPSI-QCONV(I))
	430	!
	431	! * Find level of minimum moist static energy *
	432	!
	433	IF(I.GE.MINORIG.AND.HM(I).LT.AHMIN.AND.HM(I).LT.HM(I-1))THEN
	434	AHMIN=HM(I)
	435	IHMIN=I
	436	END IF
	437	END DO
	438	IHMIN=MIN(IHMIN, NL-1)
	439	!
	440	! * Find that model level below the level of minimum moist *
	441	! * static energy that has the maximum value of moist static energy *
	442	!
	443	AHMAX=0.0
	444	! *** bug fixed: need to assign an initial value to NK
	445	! HSO, 05.08.2009
	446	NK=MINORIG
	447	DO I=MINORIG,IHMIN
	448	IF(HM(I).GT.AHMAX)THEN
	449	NK=I
	450	AHMAX=HM(I)
	451	END IF
	452	END DO
	453	!
	454	! * CHECK WHETHER PARCEL LEVEL TEMPERATURE AND SPECIFIC HUMIDITY *
	455	! * ARE REASONABLE *
	456	! * Skip convection if HM increases monotonically upward *
	457	!
	458	IF(TCONV(NK).LT.250.0.OR.QCONV(NK).LE.0.0.OR.IHMIN.EQ.(NL-1)) &
	459	THEN
	460	IFLAG=0
	461	CBMF=0.0
	462	RETURN
	463	END IF
	464	!
	465	! * CALCULATE LIFTED CONDENSATION LEVEL OF AIR AT PARCEL ORIGIN LEVEL *
	466	! * (WITHIN 0.2% OF FORMULA OF BOLTON, MON. WEA. REV.,1980) *
	467	!
	468	RH=QCONV(NK)/QSCONV(NK)
	469	CHI=TCONV(NK)/(1669.0-122.0*RH-TCONV(NK))
	470	PLCL=PCONV_HPA(NK)(RH*CHI)
	471	IF(PLCL.LT.200.0.OR.PLCL.GE.2000.0)THEN
	472	IFLAG=2
	473	CBMF=0.0
	474	RETURN
	475	END IF
	476	!
	477	! * CALCULATE FIRST LEVEL ABOVE LCL (=ICB) *
	478	!
	479	ICB=NL-1
	480	DO I=NK+1,NL
	481	IF(PCONV_HPA(I).LT.PLCL)THEN
	482	ICB=MIN(ICB,I)
	483	END IF
	484	END DO
	485	IF(ICB.GE.(NL-1))THEN
	486	IFLAG=3
	487	CBMF=0.0
	488	RETURN
	489	END IF
	490	!
	491	! * FIND TEMPERATURE UP THROUGH ICB AND TEST FOR INSTABILITY *
	492	!
	493	! * SUBROUTINE TLIFT CALCULATES PART OF THE LIFTED PARCEL VIRTUAL *
	494	! * TEMPERATURE, THE ACTUAL TEMPERATURE AND THE ADIABATIC *
	495	! * LIQUID WATER CONTENT *
	496	!
	497	CALL TLIFT(GZ,ICB,NK,TVP,TP,CLW,ND,NL,1)
	498	DO I=NK,ICB
	499	TVP(I)=TVP(I)-TP(I)*QCONV(NK)
	500	END DO
	501	!
	502	! * If there was no convection at last time step and parcel *
	503	! * is stable at ICB then skip rest of calculation *
	504	!
	505	IF(CBMF.EQ.0.0.AND.TVP(ICB).LE.(TV(ICB)-DTMAX))THEN
	506	IFLAG=0
	507	RETURN
	508	END IF
	509	!
	510	! * IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY *
	511	!
	512	IF(IFLAG.NE.4)IFLAG=1
	513	!
	514	! * FIND THE REST OF THE LIFTED PARCEL TEMPERATURES *
	515	!
	516	CALL TLIFT(GZ,ICB,NK,TVP,TP,CLW,ND,NL,2)
	517	!
	518	! * SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF *
	519	! * PRECIPITATION FALLING OUTSIDE OF CLOUD *
	520	! * THESE MAY BE FUNCTIONS OF TP(I), PCONV_HPA(I) AND CLW(I) *
	521	!
	522	DO I=1,NK
	523	EP(I)=0.0
	524	SIGP(I)=SIGS
	525	END DO
	526	DO I=NK+1,NL
	527	TCA=TP(I)-273.15
	528	IF(TCA.GE.0.0)THEN
	529	ELACRIT=ELCRIT
	530	ELSE
	531	ELACRIT=ELCRIT*(1.0-TCA/TLCRIT)
	532	END IF
	533	ELACRIT=MAX(ELACRIT,0.0)
	534	EPMAX=0.999
	535	EP(I)=EPMAX*(1.0-ELACRIT/MAX(CLW(I),1.0E-8))
	536	EP(I)=MAX(EP(I),0.0)
	537	EP(I)=MIN(EP(I),EPMAX)
	538	SIGP(I)=SIGS
	539	END DO
	540	!
	541	! * CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL *
	542	! * VIRTUAL TEMPERATURE *
	543	!
	544	DO I=ICB+1,NL
	545	TVP(I)=TVP(I)-TP(I)*QCONV(NK)
	546	END DO
	547	TVP(NL+1)=TVP(NL)-(GZ(NL+1)-GZ(NL))/CPD
	548	!
	549	! * NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS *
	550	!
	551	DO I=1,NL+1
	552	HP(I)=H(I)
	553	NENT(I)=0
	554	WATER(I)=0.0
	555	EVAP(I)=0.0
	556	WT(I)=OMTSNOW
	557	LVCP(I)=LV(I)/CPN(I)
	558	DO J=1,NL+1
	559	QENT(I,J)=QCONV(J)
	560	ELIJ(I,J)=0.0
	561	SIJ(I,J)=0.0
	562	END DO
	563	END DO
	564	QP(1)=QCONV(1)
	565	DO I=2,NL+1
	566	QP(I)=QCONV(I-1)
	567	END DO
	568	!
	569	! * FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S *
	570	! * HIGHEST LEVEL OF NEUTRAL BUOYANCY *
	571	! * AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB) *
	572	!
	573	CAPE=0.0
	574	CAPEM=0.0
	575	INB=ICB+1
	576	INB1=INB
	577	BYP=0.0
	578	DO I=ICB+1,NL-1
	579	BY=(TVP(I)-TV(I))*(PHCONV_HPA(I)-PHCONV_HPA(I+1))/PCONV_HPA(I)
	580	CAPE=CAPE+BY
	581	IF(BY.GE.0.0)INB1=I+1
	582	IF(CAPE.GT.0.0)THEN
	583	INB=I+1
	584	BYP=(TVP(I+1)-TV(I+1))*(PHCONV_HPA(I+1)-PHCONV_HPA(I+2))/ &
	585	PCONV_HPA(I+1)
	586	CAPEM=CAPE
	587	END IF
	588	END DO
	589	INB=MAX(INB,INB1)
	590	CAPE=CAPEM+BYP
	591	DEFRAC=CAPEM-CAPE
	592	DEFRAC=MAX(DEFRAC,0.001)
	593	FRAC=-CAPE/DEFRAC
	594	FRAC=MIN(FRAC,1.0)
	595	FRAC=MAX(FRAC,0.0)
	596	!
	597	! * CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL *
	598	!
	599	DO I=ICB,INB
	600	HP(I)=H(NK)+(LV(I)+(CPD-CPV)TCONV(I))EP(I)*CLW(I)
	601	END DO
	602	!
	603	! * CALCULATE CLOUD BASE MASS FLUX AND RATES OF MIXING, M(I), *
	604	! * AT EACH MODEL LEVEL *
	605	!
	606	DBOSUM=0.0
	607	!
	608	! * INTERPOLATE DIFFERENCE BETWEEN LIFTED PARCEL AND *
	609	! * ENVIRONMENTAL TEMPERATURES TO LIFTED CONDENSATION LEVEL *
	610	!
	611	TVPPLCL=TVP(ICB-1)-RDTVP(ICB-1)(PCONV_HPA(ICB-1)-PLCL)/ &
	612	(CPN(ICB-1)*PCONV_HPA(ICB-1))
	613	TVAPLCL=TV(ICB)+(TVP(ICB)-TVP(ICB+1))*(PLCL-PCONV_HPA(ICB))/ &
	614	(PCONV_HPA(ICB)-PCONV_HPA(ICB+1))
	615	DTPBL=0.0
	616	DO I=NK,ICB-1
	617	DTPBL=DTPBL+(TVP(I)-TV(I))*(PHCONV_HPA(I)-PHCONV_HPA(I+1))
	618	END DO
	619	DTPBL=DTPBL/(PHCONV_HPA(NK)-PHCONV_HPA(ICB))
	620	DTMIN=TVPPLCL-TVAPLCL+DTMAX+DTPBL
	621	DTMA=DTMIN
	622	!
	623	! * ADJUST CLOUD BASE MASS FLUX *
	624	!
	625	CBMFOLD=CBMF
	626	! *** C. Forster: adjustment of CBMF is not allowed to depend on FLEXPART timestep
	627	DELT0=DELT/3.
	628	DAMPS=DAMP*DELT/DELT0
	629	CBMF=(1.-DAMPS)CBMF+0.1ALPHA*DTMA
	630	CBMF=MAX(CBMF,0.0)
	631	!
	632	! * If cloud base mass flux is zero, skip rest of calculation *
	633	!
	634	IF(CBMF.EQ.0.0.AND.CBMFOLD.EQ.0.0)THEN
	635	RETURN
	636	END IF
	637
	638	!
	639	! * CALCULATE RATES OF MIXING, M(I) *
	640	!
	641	M(ICB)=0.0
	642	DO I=ICB+1,INB
	643	K=MIN(I,INB1)
	644	DBO=ABS(TV(K)-TVP(K))+ &
	645	ENTP0.02(PHCONV_HPA(K)-PHCONV_HPA(K+1))
	646	DBOSUM=DBOSUM+DBO
	647	M(I)=CBMF*DBO
	648	END DO
	649	DO I=ICB+1,INB
	650	M(I)=M(I)/DBOSUM
	651	END DO
	652	!
	653	! * CALCULATE ENTRAINED AIR MASS FLUX (MENT), TOTAL WATER MIXING *
	654	! * RATIO (QENT), TOTAL CONDENSED WATER (ELIJ), AND MIXING *
	655	! * FRACTION (SIJ) *
	656	!
	657	DO I=ICB+1,INB
	658	QTI=QCONV(NK)-EP(I)*CLW(I)
	659	DO J=ICB,INB
	660	BF2=1.+LV(J)LV(J)QSCONV(J)/(RVTCONV(J)TCONV(J)*CPD)
	661	ANUM=H(J)-HP(I)+(CPV-CPD)TCONV(J)(QTI-QCONV(J))
	662	DENOM=H(I)-HP(I)+(CPD-CPV)(QCONV(I)-QTI)TCONV(J)
	663	DEI=DENOM
	664	IF(ABS(DEI).LT.0.01)DEI=0.01
	665	SIJ(I,J)=ANUM/DEI
	666	SIJ(I,I)=1.0
	667	ALTEM=SIJ(I,J)QCONV(I)+(1.-SIJ(I,J))QTI-QSCONV(J)
	668	ALTEM=ALTEM/BF2
	669	CWAT=CLW(J)*(1.-EP(J))
	670	STEMP=SIJ(I,J)
	671	IF((STEMP.LT.0.0.OR.STEMP.GT.1.0.OR. &
	672	ALTEM.GT.CWAT).AND.J.GT.I)THEN
	673	ANUM=ANUM-LV(J)(QTI-QSCONV(J)-CWATBF2)
	674	DENOM=DENOM+LV(J)*(QCONV(I)-QTI)
	675	IF(ABS(DENOM).LT.0.01)DENOM=0.01
	676	SIJ(I,J)=ANUM/DENOM
	677	ALTEM=SIJ(I,J)QCONV(I)+(1.-SIJ(I,J))QTI-QSCONV(J)
	678	ALTEM=ALTEM-(BF2-1.)*CWAT
	679	END IF
	680	IF(SIJ(I,J).GT.0.0.AND.SIJ(I,J).LT.0.9)THEN
	681	QENT(I,J)=SIJ(I,J)QCONV(I)+(1.-SIJ(I,J))QTI
	682	ELIJ(I,J)=ALTEM
	683	ELIJ(I,J)=MAX(0.0,ELIJ(I,J))
	684	MENT(I,J)=M(I)/(1.-SIJ(I,J))
	685	NENT(I)=NENT(I)+1
	686	END IF
	687	SIJ(I,J)=MAX(0.0,SIJ(I,J))
	688	SIJ(I,J)=MIN(1.0,SIJ(I,J))
	689	END DO
	690	!
	691	! * IF NO AIR CAN ENTRAIN AT LEVEL I ASSUME THAT UPDRAFT DETRAINS *
	692	! * AT THAT LEVEL AND CALCULATE DETRAINED AIR FLUX AND PROPERTIES *
	693	!
	694	IF(NENT(I).EQ.0)THEN
	695	MENT(I,I)=M(I)
	696	QENT(I,I)=QCONV(NK)-EP(I)*CLW(I)
	697	ELIJ(I,I)=CLW(I)
	698	SIJ(I,I)=1.0
	699	END IF
	700	END DO
	701	SIJ(INB,INB)=1.0
	702	!
	703	! * NORMALIZE ENTRAINED AIR MASS FLUXES TO REPRESENT EQUAL *
	704	! * PROBABILITIES OF MIXING *
	705	!
	706	DO I=ICB+1,INB
	707	IF(NENT(I).NE.0)THEN
	708	QP1=QCONV(NK)-EP(I)*CLW(I)
	709	ANUM=H(I)-HP(I)-LV(I)*(QP1-QSCONV(I))
	710	DENOM=H(I)-HP(I)+LV(I)*(QCONV(I)-QP1)
	711	IF(ABS(DENOM).LT.0.01)DENOM=0.01
	712	SCRIT=ANUM/DENOM
	713	ALT=QP1-QSCONV(I)+SCRIT*(QCONV(I)-QP1)
	714	IF(ALT.LT.0.0)SCRIT=1.0
	715	SCRIT=MAX(SCRIT,0.0)
	716	ASIJ=0.0
	717	SMIN=1.0
	718	DO J=ICB,INB
	719	IF(SIJ(I,J).GT.0.0.AND.SIJ(I,J).LT.0.9)THEN
	720	IF(J.GT.I)THEN
	721	SMID=MIN(SIJ(I,J),SCRIT)
	722	SJMAX=SMID
	723	SJMIN=SMID
	724	IF(SMID.LT.SMIN.AND.SIJ(I,J+1).LT.SMID)THEN
	725	SMIN=SMID
	726	SJMAX=MIN(SIJ(I,J+1),SIJ(I,J),SCRIT)
	727	SJMIN=MAX(SIJ(I,J-1),SIJ(I,J))
	728	SJMIN=MIN(SJMIN,SCRIT)
	729	END IF
	730	ELSE
	731	SJMAX=MAX(SIJ(I,J+1),SCRIT)
	732	SMID=MAX(SIJ(I,J),SCRIT)
	733	SJMIN=0.0
	734	IF(J.GT.1)SJMIN=SIJ(I,J-1)
	735	SJMIN=MAX(SJMIN,SCRIT)
	736	END IF
	737	DELP=ABS(SJMAX-SMID)
	738	DELM=ABS(SJMIN-SMID)
	739	ASIJ=ASIJ+(DELP+DELM)*(PHCONV_HPA(J)-PHCONV_HPA(J+1))
	740	MENT(I,J)=MENT(I,J)(DELP+DELM) &
	741	(PHCONV_HPA(J)-PHCONV_HPA(J+1))
	742	END IF
	743	END DO
	744	ASIJ=MAX(1.0E-21,ASIJ)
	745	ASIJ=1.0/ASIJ
	746	DO J=ICB,INB
	747	MENT(I,J)=MENT(I,J)*ASIJ
	748	END DO
	749	BSUM=0.0
	750	DO J=ICB,INB
	751	BSUM=BSUM+MENT(I,J)
	752	END DO
	753	IF(BSUM.LT.1.0E-18)THEN
	754	NENT(I)=0
	755	MENT(I,I)=M(I)
	756	QENT(I,I)=QCONV(NK)-EP(I)*CLW(I)
	757	ELIJ(I,I)=CLW(I)
	758	SIJ(I,I)=1.0
	759	END IF
	760	END IF
	761	END DO
	762	!
	763	! * CHECK WHETHER EP(INB)=0, IF SO, SKIP PRECIPITATING *
	764	! * DOWNDRAFT CALCULATION *
	765	!
	766	IF(EP(INB).LT.0.0001)GOTO 405
	767	!
	768	! * INTEGRATE LIQUID WATER EQUATION TO FIND CONDENSED WATER *
	769	! * AND CONDENSED WATER FLUX *
	770	!
	771	JTT=2
	772	!
	773	! * BEGIN DOWNDRAFT LOOP *
	774	!
	775	DO I=INB,1,-1
	776	!
	777	! * CALCULATE DETRAINED PRECIPITATION *
	778	!
	779	WDTRAIN=GEP(I)M(I)*CLW(I)
	780	IF(I.GT.1)THEN
	781	DO J=1,I-1
	782	AWAT=ELIJ(J,I)-(1.-EP(I))*CLW(I)
	783	AWAT=MAX(0.0,AWAT)
	784	WDTRAIN=WDTRAIN+GAWATMENT(J,I)
	785	END DO
	786	END IF
	787	!
	788	! * FIND RAIN WATER AND EVAPORATION USING PROVISIONAL *
	789	! * ESTIMATES OF QP(I)AND QP(I-1) *
	790	!
	791	!
	792	! * Value of terminal velocity and coefficient of evaporation for snow *
	793	!
	794	COEFF=COEFFS
	795	WT(I)=OMTSNOW
	796	!
	797	! * Value of terminal velocity and coefficient of evaporation for rain *
	798	!
	799	IF(TCONV(I).GT.273.0)THEN
	800	COEFF=COEFFR
	801	WT(I)=OMTRAIN
	802	END IF
	803	QSM=0.5*(QCONV(I)+QP(I+1))
	804	AFAC=COEFFPHCONV_HPA(I)(QSCONV(I)-QSM)/ &
	805	(1.0E4+2.0E3PHCONV_HPA(I)QSCONV(I))
	806	AFAC=MAX(AFAC,0.0)
	807	SIGT=SIGP(I)
	808	SIGT=MAX(0.0,SIGT)
	809	SIGT=MIN(1.0,SIGT)
	810	B6=100.(PHCONV_HPA(I)-PHCONV_HPA(I+1))SIGT*AFAC/WT(I)
	811	C6=(WATER(I+1)*WT(I+1)+WDTRAIN/SIGD)/WT(I)
	812	REVAP=0.5(-B6+SQRT(B6B6+4.*C6))
	813	EVAP(I)=SIGTAFACREVAP
	814	WATER(I)=REVAP*REVAP
	815	!
	816	! * CALCULATE PRECIPITATING DOWNDRAFT MASS FLUX UNDER *
	817	! * HYDROSTATIC APPROXIMATION *
	818	!
	819	IF(I.EQ.1)GOTO 360
	820	DHDP=(H(I)-H(I-1))/(PCONV_HPA(I-1)-PCONV_HPA(I))
	821	DHDP=MAX(DHDP,10.0)
	822	MP(I)=100.GINVLV(I)SIGDEVAP(I)/DHDP
	823	MP(I)=MAX(MP(I),0.0)
	824	!
	825	! * ADD SMALL AMOUNT OF INERTIA TO DOWNDRAFT *
	826	!
	827	FAC=20.0/(PHCONV_HPA(I-1)-PHCONV_HPA(I))
	828	MP(I)=(FAC*MP(I+1)+MP(I))/(1.+FAC)
	829	!
	830	! * FORCE MP TO DECREASE LINEARLY TO ZERO *
	831	! * BETWEEN ABOUT 950 MB AND THE SURFACE *
	832	!
	833	IF(PCONV_HPA(I).GT.(0.949*PCONV_HPA(1)))THEN
	834	JTT=MAX(JTT,I)
	835	MP(I)=MP(JTT)*(PCONV_HPA(1)-PCONV_HPA(I))/(PCONV_HPA(1)- &
	836	PCONV_HPA(JTT))
	837	END IF
	838	360 CONTINUE
	839	!
	840	! * FIND MIXING RATIO OF PRECIPITATING DOWNDRAFT *
	841	!
	842	IF(I.EQ.INB)GOTO 400
	843	IF(I.EQ.1)THEN
	844	QSTM=QSCONV(1)
	845	ELSE
	846	QSTM=QSCONV(I-1)
	847	END IF
	848	IF(MP(I).GT.MP(I+1))THEN
	849	RAT=MP(I+1)/MP(I)
	850	QP(I)=QP(I+1)RAT+QCONV(I)(1.0-RAT)+100.GINV &
	851	SIGD(PHCONV_HPA(I)-PHCONV_HPA(I+1))(EVAP(I)/MP(I))
	852	ELSE
	853	IF(MP(I+1).GT.0.0)THEN
	854	QP(I)=(GZ(I+1)-GZ(I)+QP(I+1)(LV(I+1)+TCONV(I+1)(CL-CPD))+ &
	855	CPD(TCONV(I+1)-TCONV(I)))/(LV(I)+TCONV(I)(CL-CPD))
	856	END IF
	857	END IF
	858	QP(I)=MIN(QP(I),QSTM)
	859	QP(I)=MAX(QP(I),0.0)
	860	400 CONTINUE
	861	END DO
	862	!
	863	! * CALCULATE SURFACE PRECIPITATION IN MM/DAY *
	864	!
	865	PRECIP=PRECIP+WT(1)SIGDWATER(1)3600.24000./(ROWL*G)
	866	!
	867	405 CONTINUE
	868	!
	869	! * CALCULATE DOWNDRAFT VELOCITY SCALE AND SURFACE TEMPERATURE AND *
	870	! * WATER VAPOR FLUCTUATIONS *
	871	!
	872	WD=BETAABS(MP(ICB))0.01RDTCONV(ICB)/(SIGD*PCONV_HPA(ICB))
	873	QPRIME=0.5*(QP(1)-QCONV(1))
	874	TPRIME=LV0*QPRIME/CPD
	875	!
	876	! * CALCULATE TENDENCIES OF LOWEST LEVEL POTENTIAL TEMPERATURE *
	877	! * AND MIXING RATIO *
	878	!
	879	DPINV=0.01/(PHCONV_HPA(1)-PHCONV_HPA(2))
	880	AM=0.0
	881	IF(NK.EQ.1)THEN
	882	DO K=2,INB
	883	AM=AM+M(K)
	884	END DO
	885	END IF
	886	! save saturated upward mass flux for first level
	887	FUP(1)=AM
	888	IF((2.GDPINV*AM).GE.DELTI)IFLAG=4
	889	FT(1)=FT(1)+GDPINVAM*(TCONV(2)-TCONV(1)+(GZ(2)-GZ(1))/CPN(1))
	890	FT(1)=FT(1)-LVCP(1)SIGDEVAP(1)
	891	FT(1)=FT(1)+SIGDWT(2)(CL-CPD)WATER(2)(TCONV(2)- &
	892	TCONV(1))*DPINV/CPN(1)
	893	FQ(1)=FQ(1)+GMP(2)(QP(2)-QCONV(1))* &
	894	DPINV+SIGD*EVAP(1)
	895	FQ(1)=FQ(1)+GAM(QCONV(2)-QCONV(1))*DPINV
	896	DO J=2,INB
	897	FQ(1)=FQ(1)+GDPINVMENT(J,1)*(QENT(J,1)-QCONV(1))
	898	END DO
	899	!
	900	! * CALCULATE TENDENCIES OF POTENTIAL TEMPERATURE AND MIXING RATIO *
	901	! * AT LEVELS ABOVE THE LOWEST LEVEL *
	902	!
	903	! * FIRST FIND THE NET SATURATED UPDRAFT AND DOWNDRAFT MASS FLUXES *
	904	! * THROUGH EACH LEVEL *
	905	!
	906	DO I=2,INB
	907	DPINV=0.01/(PHCONV_HPA(I)-PHCONV_HPA(I+1))
	908	CPINV=1.0/CPN(I)
	909	AMP1=0.0
	910	AD=0.0
	911	IF(I.GE.NK)THEN
	912	DO K=I+1,INB+1
	913	AMP1=AMP1+M(K)
	914	END DO
	915	END IF
	916	DO K=1,I
	917	DO J=I+1,INB+1
	918	AMP1=AMP1+MENT(K,J)
	919	END DO
	920	END DO
	921	! save saturated upward mass flux
	922	FUP(I)=AMP1
	923	IF((2.GDPINV*AMP1).GE.DELTI)IFLAG=4
	924	DO K=1,I-1
	925	DO J=I,INB
	926	AD=AD+MENT(J,K)
	927	END DO
	928	END DO
	929	! save saturated downward mass flux
	930	FDOWN(I)=AD
	931	FT(I)=FT(I)+GDPINV(AMP1(TCONV(I+1)-TCONV(I)+(GZ(I+1)-GZ(I)) &
	932	CPINV)-AD(TCONV(I)-TCONV(I-1)+(GZ(I)-GZ(I-1))CPINV)) &
	933	-SIGDLVCP(I)EVAP(I)
	934	FT(I)=FT(I)+GDPINVMENT(I,I)*(HP(I)-H(I)+ &
	935	TCONV(I)(CPV-CPD)(QCONV(I)-QENT(I,I)))*CPINV
	936	FT(I)=FT(I)+SIGDWT(I+1)(CL-CPD)WATER(I+1) &
	937	(TCONV(I+1)-TCONV(I))DPINVCPINV
	938	FQ(I)=FQ(I)+GDPINV(AMP1*(QCONV(I+1)-QCONV(I))- &
	939	AD*(QCONV(I)-QCONV(I-1)))
	940	DO K=1,I-1
	941	AWAT=ELIJ(K,I)-(1.-EP(I))*CLW(I)
	942	AWAT=MAX(AWAT,0.0)
	943	FQ(I)=FQ(I)+GDPINVMENT(K,I)*(QENT(K,I)-AWAT-QCONV(I))
	944	END DO
	945	DO K=I,INB
	946	FQ(I)=FQ(I)+GDPINVMENT(K,I)*(QENT(K,I)-QCONV(I))
	947	END DO
	948	FQ(I)=FQ(I)+SIGDEVAP(I)+G(MP(I+1)* &
	949	(QP(I+1)-QCONV(I))-MP(I)(QP(I)-QCONV(I-1)))DPINV
	950	END DO
	951	!
	952	! * Adjust tendencies at top of convection layer to reflect *
	953	! * actual position of the level zero CAPE *
	954	!
	955	FQOLD=FQ(INB)
	956	FQ(INB)=FQ(INB)*(1.-FRAC)
	957	FQ(INB-1)=FQ(INB-1)+FRACFQOLD((PHCONV_HPA(INB)- &
	958	PHCONV_HPA(INB+1))/ &
	959	(PHCONV_HPA(INB-1)-PHCONV_HPA(INB)))*LV(INB)/LV(INB-1)
	960	FTOLD=FT(INB)
	961	FT(INB)=FT(INB)*(1.-FRAC)
	962	FT(INB-1)=FT(INB-1)+FRACFTOLD((PHCONV_HPA(INB)- &
	963	PHCONV_HPA(INB+1))/ &
	964	(PHCONV_HPA(INB-1)-PHCONV_HPA(INB)))*CPN(INB)/CPN(INB-1)
	965	!
	966	! * Very slightly adjust tendencies to force exact *
	967	! * enthalpy, momentum and tracer conservation *
	968	!
	969	ENTS=0.0
	970	DO I=1,INB
	971	ENTS=ENTS+(CPN(I)FT(I)+LV(I)FQ(I))* &
	972	(PHCONV_HPA(I)-PHCONV_HPA(I+1))
	973	END DO
	974	ENTS=ENTS/(PHCONV_HPA(1)-PHCONV_HPA(INB+1))
	975	DO I=1,INB
	976	FT(I)=FT(I)-ENTS/CPN(I)
	977	END DO
	978
	979	! ************************************************
	980	! **** DETERMINE MASS DISPLACEMENT MATRIX
	981	! ***** AND COMPENSATING SUBSIDENCE
	982	! ************************************************
	983
	984	! mass displacement matrix due to saturated up-and downdrafts
	985	! inside the cloud and determine compensating subsidence
	986	! FUP(I) (saturated updrafts), FDOWN(I) (saturated downdrafts) are assumed to be
	987	! balanced by compensating subsidence (SUB(I))
	988	! FDOWN(I) and SUB(I) defined positive downwards
	989
	990	! NCONVTOP IS THE TOP LEVEL AT WHICH CONVECTIVE MASS FLUXES ARE DIAGNOSED
	991	! EPSILON IS A SMALL NUMBER
	992
	993	SUB(1)=0.
	994	NCONVTOP=1
	995	do i=1,INB+1
	996	do j=1,INB+1
	997	if (j.eq.NK) then
	998	FMASS(j,i)=FMASS(j,i)+M(i)
	999	endif
	1000	FMASS(j,i)=FMASS(j,i)+MENT(j,i)
	1001	IF (FMASS(J,I).GT.EPSILON) NCONVTOP=MAX(NCONVTOP,I,J)
	1002	end do
	1003	if (i.gt.1) then
	1004	SUB(i)=FUP(i-1)-FDOWN(i)
	1005	endif
	1006	end do
	1007	NCONVTOP=NCONVTOP+1
	1008
	1009	RETURN
	1010	!
	1011	END SUBROUTINE CONVECT
	1012	!
	1013	! ---------------------------------------------------------------------------
	1014	!
	1015	SUBROUTINE TLIFT(GZ,ICB,NK,TVP,TPK,CLW,ND,NL,KK)
	1016	!
	1017	!-cv
	1018	use par_mod
	1019	use conv_mod
	1020
	1021	implicit none
	1022	!-cv
	1023	!====>Begin Module TLIFT File convect.f Undeclared variables
	1024	!
	1025	!Argument variables
	1026	!
	1027	integer :: icb, kk, nd, nk, nl
	1028	!
	1029	!Local variables
	1030	!
	1031	integer :: i, j, nsb, nst
	1032	!
	1033	real :: ah0, ahg, alv, cpinv, cpp, denom
	1034	real :: es, qg, rg, s, tc, tg
	1035	!
	1036	!====>End Module TLIFT File convect.f
	1037
	1038	REAL :: GZ(ND),TPK(ND),CLW(ND)
	1039	REAL :: TVP(ND)
	1040	!
	1041	! * ASSIGN VALUES OF THERMODYNAMIC CONSTANTS *
	1042	!
	1043	REAL,PARAMETER :: CPD=1005.7
	1044	REAL,PARAMETER :: CPV=1870.0
	1045	REAL,PARAMETER :: CL=2500.0
	1046	REAL,PARAMETER :: RV=461.5
	1047	REAL,PARAMETER :: RD=287.04
	1048	REAL,PARAMETER :: LV0=2.501E6
	1049	!
	1050	REAL,PARAMETER :: CPVMCL=CL-CPV
	1051	REAL,PARAMETER :: EPS0=RD/RV
	1052	REAL,PARAMETER :: EPSI=1./EPS0
	1053	!
	1054	! * CALCULATE CERTAIN PARCEL QUANTITIES, INCLUDING STATIC ENERGY *
	1055	!
	1056	AH0=(CPD(1.-QCONV(NK))+CLQCONV(NK))TCONV(NK)+QCONV(NK) &
	1057	(LV0-CPVMCL*( &
	1058	TCONV(NK)-273.15))+GZ(NK)
	1059	CPP=CPD(1.-QCONV(NK))+QCONV(NK)CPV
	1060	CPINV=1./CPP
	1061	!
	1062	IF(KK.EQ.1)THEN
	1063	!
	1064	! * CALCULATE LIFTED PARCEL QUANTITIES BELOW CLOUD BASE *
	1065	!
	1066	DO I=1,ICB-1
	1067	CLW(I)=0.0
	1068	END DO
	1069	DO I=NK,ICB-1
	1070	TPK(I)=TCONV(NK)-(GZ(I)-GZ(NK))*CPINV
	1071	TVP(I)=TPK(I)(1.+QCONV(NK)EPSI)
	1072	END DO
	1073	END IF
	1074	!
	1075	! * FIND LIFTED PARCEL QUANTITIES ABOVE CLOUD BASE *
	1076	!
	1077	NST=ICB
	1078	NSB=ICB
	1079	IF(KK.EQ.2)THEN
	1080	NST=NL
	1081	NSB=ICB+1
	1082	END IF
	1083	DO I=NSB,NST
	1084	TG=TCONV(I)
	1085	QG=QSCONV(I)
	1086	ALV=LV0-CPVMCL*(TCONV(I)-273.15)
	1087	DO J=1,2
	1088	S=CPD+ALVALVQG/(RVTCONV(I)TCONV(I))
	1089	S=1./S
	1090	AHG=CPDTG+(CL-CPD)QCONV(NK)TCONV(I)+ALVQG+GZ(I)
	1091	TG=TG+S*(AH0-AHG)
	1092	TG=MAX(TG,35.0)
	1093	TC=TG-273.15
	1094	DENOM=243.5+TC
	1095	IF(TC.GE.0.0)THEN
	1096	ES=6.112EXP(17.67TC/DENOM)
	1097	ELSE
	1098	ES=EXP(23.33086-6111.72784/TG+0.15215*LOG(TG))
	1099	END IF
	1100	QG=EPS0ES/(PCONV_HPA(I)-ES(1.-EPS0))
	1101	END DO
	1102	ALV=LV0-CPVMCL*(TCONV(I)-273.15)
	1103	TPK(I)=(AH0-(CL-CPD)QCONV(NK)TCONV(I)-GZ(I)-ALV*QG)/CPD
	1104	CLW(I)=QCONV(NK)-QG
	1105	CLW(I)=MAX(0.0,CLW(I))
	1106	RG=QG/(1.-QCONV(NK))
	1107	TVP(I)=TPK(I)(1.+RGEPSI)
	1108	END DO
	1109	RETURN
	1110	END SUBROUTINE TLIFT

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

Download in other formats:

Original Format