Context Navigation

convect43c.f90 @ 3481cc1

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug

Last change on this file since 3481cc1 was 3481cc1, checked in by Ignacio Pisso <ip@…>, 4 years ago
move license from headers to a different file
Property mode set to `100644`
File size: 35.0 KB

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

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

Download in other formats:

Original Format