Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

convect43c.f90 @ 8

Last change on this file since 8 was 8, checked in by hasod, 11 years ago
Added parallel version of Flexpart91
File size: 36.5 KB

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

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

Download in other formats:

Original Format