Context Navigation

convect43c.f90 @ 3405994

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

Last change on this file since 3405994 was 92fab65, checked in by Ignacio Pisso <ip@…>, 4 years ago
add SPDX-License-Identifier to all .f90 files
Property mode set to `100644`
File size: 35.1 KB

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

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

Download in other formats:

Original Format