source: trunk/src/convect43c.f90 @ 28

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