source: branches/flexpart91_hasod/src_parallel/convect43c.f90 @ 8

Last change on this file since 8 was 8, checked in by hasod, 9 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.001*PCONV_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))/(RV*TOLD(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      CL*QCONV(J)+QSCONV(J)*(CPV-CL+ALV*ALV/(RV*TCONV(J)*TCONV(J))))
397  !     ALVNEW=LV0-CPVMCL*(TNEW-273.15)
398  !     QNEW=(ALV*QCONV(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)/(G*DELT*ROWL)
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.5*RD*(TVX+TVY)*(PCONV_HPA(I-1)-PCONV_HPA(I))/ &
425           PHCONV_HPA(I)
426      CPN(I)=CPD*(1.-QCONV(I))+CPV*QCONV(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))+CL*QCONV(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)-RD*TVP(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.1*ALPHA*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        ENTP*0.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)/(RV*TCONV(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)-CWAT*BF2)
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=G*EP(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+G*AWAT*MENT(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=COEFF*PHCONV_HPA(I)*(QSCONV(I)-QSM)/ &
807         (1.0E4+2.0E3*PHCONV_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(B6*B6+4.*C6))
815    EVAP(I)=SIGT*AFAC*REVAP
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.*GINV*LV(I)*SIGD*EVAP(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)
862400 CONTINUE
863    END DO
864  !
865  !   ***  CALCULATE SURFACE PRECIPITATION IN MM/DAY     ***
866  !
867    PRECIP=PRECIP+WT(1)*SIGD*WATER(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=BETA*ABS(MP(ICB))*0.01*RD*TCONV(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.*G*DPINV*AM).GE.DELTI)IFLAG=4
891    FT(1)=FT(1)+G*DPINV*AM*(TCONV(2)-TCONV(1)+(GZ(2)-GZ(1))/CPN(1))
892    FT(1)=FT(1)-LVCP(1)*SIGD*EVAP(1)
893    FT(1)=FT(1)+SIGD*WT(2)*(CL-CPD)*WATER(2)*(TCONV(2)- &
894         TCONV(1))*DPINV/CPN(1)
895    FQ(1)=FQ(1)+G*MP(2)*(QP(2)-QCONV(1))* &
896         DPINV+SIGD*EVAP(1)
897    FQ(1)=FQ(1)+G*AM*(QCONV(2)-QCONV(1))*DPINV
898    DO J=2,INB
899     FQ(1)=FQ(1)+G*DPINV*MENT(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.*G*DPINV*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)+G*DPINV*(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         -SIGD*LVCP(I)*EVAP(I)
936    FT(I)=FT(I)+G*DPINV*MENT(I,I)*(HP(I)-H(I)+ &
937         TCONV(I)*(CPV-CPD)*(QCONV(I)-QENT(I,I)))*CPINV
938    FT(I)=FT(I)+SIGD*WT(I+1)*(CL-CPD)*WATER(I+1)* &
939         (TCONV(I+1)-TCONV(I))*DPINV*CPINV
940    FQ(I)=FQ(I)+G*DPINV*(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)+G*DPINV*MENT(K,I)*(QENT(K,I)-AWAT-QCONV(I))
946    END DO
947    DO K=I,INB
948     FQ(I)=FQ(I)+G*DPINV*MENT(K,I)*(QENT(K,I)-QCONV(I))
949    END DO
950    FQ(I)=FQ(I)+SIGD*EVAP(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)+FRAC*FQOLD*((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)+FRAC*FTOLD*((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  !
1013END SUBROUTINE CONVECT
1014!
1015! ---------------------------------------------------------------------------
1016!
1017SUBROUTINE 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))+CL*QCONV(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+ALV*ALV*QG/(RV*TCONV(I)*TCONV(I))
1090      S=1./S
1091      AHG=CPD*TG+(CL-CPD)*QCONV(NK)*TCONV(I)+ALV*QG+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.112*EXP(17.67*TC/DENOM)
1098      ELSE
1099       ES=EXP(23.33086-6111.72784/TG+0.15215*LOG(TG))
1100      END IF
1101      QG=EPS0*ES/(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.+RG*EPSI)
1109    END DO
1110    RETURN
1111END SUBROUTINE TLIFT
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG