flexdoxygen/20151015/convect43c_8f90_source.html

 !**********************************************************************

 ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *

 ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *

 ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *

 !                                                                     *

 ! This file is part of FLEXPART.                                      *

 !                                                                     *

 ! FLEXPART is free software: you can redistribute it and/or modify    *

 ! it under the terms of the GNU General Public License as published by*

 ! the Free Software Foundation, either version 3 of the License, or   *

 ! (at your option) any later version.                                 *

 !                                                                     *

 ! FLEXPART is distributed in the hope that it will be useful,         *

 ! but WITHOUT ANY WARRANTY; without even the implied warranty of      *

 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *

 ! GNU General Public License for more details.                        *

 !                                                                     *

 ! You should have received a copy of the GNU General Public License   *

 ! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *

 !**********************************************************************


 !**************************************************************************

 !****                       SUBROUTINE CONVECT                        *****

 !****                          VERSION 4.3c                           *****

 !****                          20 May, 2002                           *****

 !****                          Kerry Emanuel                          *****

 !**************************************************************************

 !

   SUBROUTINE convect &

          (nd,  nl,   delt, iflag, &

          precip, wd,   tprime, qprime, cbmf    )

   !

   !-cv *************************************************************************

   !-cv C. Forster, November 2003 - May 2004:

   !-cv

   !-cv The subroutine has been downloaded from Kerry Emanuel's homepage,

   !-cv where further infos on the convection scheme can be found

   !-cv http://www-paoc.mit.edu/~emanuel/home.html

   !-cv

   !-cv The following changes have been made to integrate this subroutine

   !-cv into FLEXPART

   !-cv

   !-cv Putting most of the variables in a new common block

   !-cv renaming eps to eps0 because there is some eps already in includepar

   !-cv

   !-cv removing the arrays U,V,TRA and related arrays

   !-cv

   !-cv renaming the original arrays T,Q,QS,P,PH to

   !-cv TCONV,QCONV,QSCONV,PCONV_HPA,PHCONV_HPA

   !-cv

   !-cv Initialization of variables has been put into parameter statements

   !-cv instead of assignment of values at each call, in order to save

   !-cv computation time.

   !***************************************************************************

   !

   !-----------------------------------------------------------------------------

   !    *** On input:      ***

   !

   !T:   Array of absolute temperature (K) of dimension ND, with first

   !      index corresponding to lowest model level. Note that this array

   !      will be altered by the subroutine if dry convective adjustment

   !      occurs and if IPBL is not equal to 0.

   !

   !Q:   Array of specific humidity (gm/gm) of dimension ND, with first

   !       index corresponding to lowest model level. Must be defined

   !       at same grid levels as T. Note that this array will be altered

   !       if dry convective adjustment occurs and if IPBL is not equal to 0.

   !

   !QS:  Array of saturation specific humidity of dimension ND, with first

   !       index corresponding to lowest model level. Must be defined

   !       at same grid levels as T. Note that this array will be altered

   !       if dry convective adjustment occurs and if IPBL is not equal to 0.

   !

   !U:   Array of zonal wind velocity (m/s) of dimension ND, witth first

   !       index corresponding with the lowest model level. Defined at

   !       same levels as T. Note that this array will be altered if

   !       dry convective adjustment occurs and if IPBL is not equal to 0.

   !

   !V:   Same as U but for meridional velocity.

   !

   !TRA: Array of passive tracer mixing ratio, of dimensions (ND,NTRA),

   !       where NTRA is the number of different tracers. If no

   !       convective tracer transport is needed, define a dummy

   !       input array of dimension (ND,1). Tracers are defined at

   !       same vertical levels as T. Note that this array will be altered

   !       if dry convective adjustment occurs and if IPBL is not equal to 0.

   !

   !P:   Array of pressure (mb) of dimension ND, with first

   !       index corresponding to lowest model level. Must be defined

   !       at same grid levels as T.

   !

   !PH:  Array of pressure (mb) of dimension ND+1, with first index

   !       corresponding to lowest level. These pressures are defined at

   !       levels intermediate between those of P, T, Q and QS. The first

   !       value of PH should be greater than (i.e. at a lower level than)

   !       the first value of the array P.

   !

   !ND:  The dimension of the arrays T,Q,QS,P,PH,FT and FQ

   !

   !NL:  The maximum number of levels to which convection can

   !       penetrate, plus 1.

   !       NL MUST be less than or equal to ND-1.

   !

   !NTRA:The number of different tracers. If no tracer transport

   !       is needed, set this equal to 1. (On most compilers, setting

   !       NTRA to 0 will bypass tracer calculation, saving some CPU.)

   !

   !DELT: The model time step (sec) between calls to CONVECT

   !

   !----------------------------------------------------------------------------

   !    ***   On Output:         ***

   !

   !IFLAG: An output integer whose value denotes the following:

   !

   !           VALUE                        INTERPRETATION

   !           -----                        --------------

   !             0               No moist convection; atmosphere is not

   !                             unstable, or surface temperature is less

   !                             than 250 K or surface specific humidity

   !                             is non-positive.

   !

   !             1               Moist convection occurs.

   !

   !             2               No moist convection: lifted condensation

   !                             level is above the 200 mb level.

   !

   !             3               No moist convection: cloud base is higher

   !                             then the level NL-1.

   !

   !             4               Moist convection occurs, but a CFL condition

   !                             on the subsidence warming is violated. This

   !                             does not cause the scheme to terminate.

   !

   !FT:   Array of temperature tendency (K/s) of dimension ND, defined at same

   !        grid levels as T, Q, QS and P.

   !

   !FQ:   Array of specific humidity tendencies ((gm/gm)/s) of dimension ND,

   !        defined at same grid levels as T, Q, QS and P.

   !

   !FU:   Array of forcing of zonal velocity (m/s^2) of dimension ND,

   !        defined at same grid levels as T.

   !

   !FV:   Same as FU, but for forcing of meridional velocity.

   !

   !FTRA: Array of forcing of tracer content, in tracer mixing ratio per

   !        second, defined at same levels as T. Dimensioned (ND,NTRA).

   !

   !PRECIP: Scalar convective precipitation rate (mm/day).

   !

   !WD:    A convective downdraft velocity scale. For use in surface

   !        flux parameterizations. See convect.ps file for details.

   !

   !TPRIME: A convective downdraft temperature perturbation scale (K).

   !         For use in surface flux parameterizations. See convect.ps

   !         file for details.

   !

   !QPRIME: A convective downdraft specific humidity

   !         perturbation scale (gm/gm).

   !         For use in surface flux parameterizations. See convect.ps

   !         file for details.

   !

   !CBMF:   The cloud base mass flux ((kg/m**2)/s). THIS SCALAR VALUE MUST

   !         BE STORED BY THE CALLING PROGRAM AND RETURNED TO CONVECT AT

   !         ITS NEXT CALL. That is, the value of CBMF must be "remembered"

   !         by the calling program between calls to CONVECT.

   !

   !-----------------------------------------------------------------------------

   !

   !    ***  THE PARAMETER NA SHOULD IN GENERAL BE GREATER THAN   ***

   !    ***                OR EQUAL TO  ND + 1                    ***

   !

   !

   use par_mod

   use conv_mod


   implicit none

   !

   !-cv====>Begin Module CONVECT    File convect.f      Undeclared variables

   !

   !Argument variables

   !

   integer :: iflag, nd, nl

   !

   real :: cbmf, delt, precip, qprime, tprime, wd

   !

   !Local variables

   !

   integer :: i, icb, ihmin, inb, inb1, j, jtt, k

   integer :: nk

   !

   real :: ad, afac, ahmax, ahmin, alt, altem

   real :: am, amp1, anum, asij, awat, b6, bf2, bsum, by

   real :: byp, c6, cape, capem, cbmfold, chi, coeff

   real :: cpinv, cwat, damps, dbo, dbosum

   real :: defrac, dei, delm, delp, delt0, delti, denom, dhdp

   real :: dpinv, dtma, dtmin, dtpbl, elacrit, ents

   real :: epmax, fac, fqold, frac, ftold

   real :: plcl, qp1, qsm, qstm, qti, rat

   real :: rdcp, revap, rh, scrit, sigt, sjmax

   real :: sjmin, smid, smin, stemp, tca

   real :: tvaplcl, tvpplcl, tvx, tvy, wdtrain


   !integer jc,jn

   !real alvnew,a2,ahm,alv,rm,sum,qnew,dphinv,tc,thbar,tnew,x


   real :: fup(na),fdown(na)

   !

   !-cv====>End Module   CONVECT    File convect.f


   INTEGER :: nent(na)

   REAL :: m(na),mp(na),ment(na,na),qent(na,na),elij(na,na)

   REAL :: sij(na,na),tvp(na),tv(na),water(na)

   REAL :: qp(na),ep(na),th(na),wt(na),evap(na),clw(na)

   REAL :: sigp(na),tp(na),cpn(na)

   REAL :: lv(na),lvcp(na),h(na),hp(na),gz(na),hm(na)

   !REAL TOLD(NA)

   !

   ! -----------------------------------------------------------------------

   !

   !   ***                     Specify Switches                         ***

   !

   !   ***   IPBL: Set to zero to bypass dry adiabatic adjustment       ***

   !   ***    Any other value results in dry adiabatic adjustment       ***

   !   ***     (Zero value recommended for use in models with           ***

   !   ***                   boundary layer schemes)                    ***

   !

   !   ***   MINORIG: Lowest level from which convection may originate  ***

   !   ***     (Should be first model level at which T is defined       ***

   !   ***      for models using bulk PBL schemes; otherwise, it should ***

   !   ***      be the first model level at which T is defined above    ***

   !   ***                      the surface layer)                      ***

   !

     INTEGER,PARAMETER :: ipbl=0

     INTEGER,PARAMETER :: minorig=1

   !

   !------------------------------------------------------------------------------

   !

   !   ***                    SPECIFY PARAMETERS                        ***

   !

   !   *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***

   !   ***  TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-        ***

   !   ***       CONVERSION THRESHOLD IS ASSUMED TO BE ZERO             ***

   !   ***     (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY            ***

   !   ***               BETWEEN 0 C AND TLCRIT)                        ***

   !   ***   ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT       ***

   !   ***                       FORMULATION                            ***

   !   ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***

   !   ***  SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE       ***

   !   ***                        OF CLOUD                              ***

   !   ***        OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN       ***

   !   ***     OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW          ***

   !   ***  COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***

   !   ***                          OF RAIN                             ***

   !   ***  COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***

   !   ***                          OF SNOW                             ***

   !   ***     CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM      ***

   !   ***                         TRANSPORT                            ***

   !   ***    DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION    ***

   !   ***        A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC      ***

   !   ***    ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF    ***

   !   ***                 APPROACH TO QUASI-EQUILIBRIUM                ***

   !   ***   (THEIR STANDARD VALUES ARE  0.20 AND 0.1, RESPECTIVELY)    ***

   !   ***                   (DAMP MUST BE LESS THAN 1)                 ***

   !

     REAL,PARAMETER :: elcrit=.0011

     REAL,PARAMETER :: tlcrit=-55.0

     REAL,PARAMETER :: entp=1.5

     REAL,PARAMETER :: sigd=0.05

     REAL,PARAMETER :: sigs=0.12

     REAL,PARAMETER :: omtrain=50.0

     REAL,PARAMETER :: omtsnow=5.5

     REAL,PARAMETER :: coeffr=1.0

     REAL,PARAMETER :: coeffs=0.8

     REAL,PARAMETER :: cu=0.7

     REAL,PARAMETER :: beta=10.0

     REAL,PARAMETER :: dtmax=0.9

     REAL,PARAMETER :: alpha=0.025  !original 0.2

     REAL,PARAMETER :: damp=0.1

   !

   !   ***        ASSIGN VALUES OF THERMODYNAMIC CONSTANTS,        ***

   !   ***            GRAVITY, AND LIQUID WATER DENSITY.           ***

   !   ***             THESE SHOULD BE CONSISTENT WITH             ***

   !   ***              THOSE USED IN CALLING PROGRAM              ***

   !   ***     NOTE: THESE ARE ALSO SPECIFIED IN SUBROUTINE TLIFT  ***

   !

   REAL,PARAMETER :: cpd=1005.7

   REAL,PARAMETER :: cpv=1870.0

   REAL,PARAMETER :: cl=2500.0

   REAL,PARAMETER :: rv=461.5

   REAL,PARAMETER :: rd=287.04

   REAL,PARAMETER :: lv0=2.501e6

   REAL,PARAMETER :: g=9.81

   REAL,PARAMETER :: rowl=1000.0

   !

   REAL,PARAMETER :: cpvmcl=cl-cpv

   REAL,PARAMETER :: eps0=rd/rv

   REAL,PARAMETER :: epsi=1./eps0

   REAL,PARAMETER :: ginv=1.0/g

   REAL,PARAMETER :: epsilon=1.e-20


   ! EPSILON IS A SMALL NUMBER USED TO EXCLUDE MASS FLUXES OF ZERO

   !

   delti=1.0/delt

   !

   !      ***  INITIALIZE OUTPUT ARRAYS AND PARAMETERS  ***

   !


     DO i=1,nl+1

      ft(i)=0.0

      fq(i)=0.0

      fdown(i)=0.0

      sub(i)=0.0

      fup(i)=0.0

      m(i)=0.0

      mp(i)=0.0

     DO j=1,nl+1

      fmass(i,j)=0.0

      ment(i,j)=0.0

     END DO

     END DO

     DO i=1,nl+1

      rdcp=(rd*(1.-qconv(i))+qconv(i)*rv)/ &

           (cpd*(1.-qconv(i))+qconv(i)*cpv)

      th(i)=tconv(i)*(1000.0/pconv_hpa(i))**rdcp

     END DO

     precip=0.0

     wd=0.0

     tprime=0.0

     qprime=0.0

     iflag=0

   !

   !  IF(IPBL.NE.0)THEN

   !

   !***            PERFORM DRY ADIABATIC ADJUSTMENT            ***

   !

   !  JC=0

   !  DO 30 I=NL-1,1,-1

   !   JN=0

   !    SUM=TH(I)*(1.+QCONV(I)*EPSI-QCONV(I))

   !   DO 10 J=I+1,NL

   !    SUM=SUM+TH(J)*(1.+QCONV(J)*EPSI-QCONV(J))

   !    THBAR=SUM/REAL(J+1-I)

   !    IF((TH(J)*(1.+QCONV(J)*EPSI-QCONV(J))).LT.THBAR)JN=J

   !  10    CONTINUE

   !   IF(I.EQ.1)JN=MAX(JN,2)

   !   IF(JN.EQ.0)GOTO 30

   !  12    CONTINUE

   !   AHM=0.0

   !   RM=0.0

   !   DO 15 J=I,JN

   !    AHM=AHM+(CPD*(1.-QCONV(J))+QCONV(J)*CPV)*TCONV(J)*

   !    +   (PHCONV_HPA(J)-PHCONV_HPA(J+1))

   !    RM=RM+QCONV(J)*(PHCONV_HPA(J)-PHCONV_HPA(J+1))

   !  15    CONTINUE

   !   DPHINV=1./(PHCONV_HPA(I)-PHCONV_HPA(JN+1))

   !   RM=RM*DPHINV

   !   A2=0.0

   !   DO 20 J=I,JN

   !    QCONV(J)=RM

   !    RDCP=(RD*(1.-QCONV(J))+QCONV(J)*RV)/

   !    1     (CPD*(1.-QCONV(J))+QCONV(J)*CPV)

   !    X=(0.001*PCONV_HPA(J))**RDCP

   !    TOLD(J)=TCONV(J)

   !    TCONV(J)=X

   !    A2=A2+(CPD*(1.-QCONV(J))+QCONV(J)*CPV)*X*

   !    1    (PHCONV_HPA(J)-PHCONV_HPA(J+1))

   !  20    CONTINUE

   !   DO 25 J=I,JN

   !    TH(J)=AHM/A2

   !    TCONV(J)=TCONV(J)*TH(J)

   !    TC=TOLD(J)-273.15

   !    ALV=LV0-CPVMCL*TC

   !    QSCONV(J)=QSCONV(J)+QSCONV(J)*(1.+QSCONV(J)*(EPSI-1.))*ALV*

   !    1    (TCONV(J)- TOLD(J))/(RV*TOLD(J)*TOLD(J))

   ! if (qslev(j) .lt. 0.) then

   !   write(*,*) 'qslev.lt.0 ',j,qslev

   ! endif

   !  25    CONTINUE

   !   IF((TH(JN+1)*(1.+QCONV(JN+1)*EPSI-QCONV(JN+1))).LT.

   !    1    (TH(JN)*(1.+QCONV(JN)*EPSI-QCONV(JN))))THEN

   !    JN=JN+1

   !    GOTO 12

   !   END IF

   !   IF(I.EQ.1)JC=JN

   !  30   CONTINUE

   !

   !   ***   Remove any supersaturation that results from adjustment ***

   !

   !IF(JC.GT.1)THEN

   ! DO 38 J=1,JC

   !    IF(QSCONV(J).LT.QCONV(J))THEN

   !     ALV=LV0-CPVMCL*(TCONV(J)-273.15)

   !     TNEW=TCONV(J)+ALV*(QCONV(J)-QSCONV(J))/(CPD*(1.-QCONV(J))+

   !    1      CL*QCONV(J)+QSCONV(J)*(CPV-CL+ALV*ALV/(RV*TCONV(J)*TCONV(J))))

   !     ALVNEW=LV0-CPVMCL*(TNEW-273.15)

   !     QNEW=(ALV*QCONV(J)-(TNEW-TCONV(J))*(CPD*(1.-QCONV(J))

   !    1     +CL*QCONV(J)))/ALVNEW

   !     PRECIP=PRECIP+24.*3600.*1.0E5*(PHCONV_HPA(J)-PHCONV_HPA(J+1))*

   !    1      (QCONV(J)-QNEW)/(G*DELT*ROWL)

   !     TCONV(J)=TNEW

   !     QCONV(J)=QNEW

   !     QSCONV(J)=QNEW

   !    END IF

   !  38  CONTINUE

   !END IF

   !

   !END IF

   !

   !  *** CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY AND STATIC ENERGY

   !

     gz(1)=0.0

     cpn(1)=cpd*(1.-qconv(1))+qconv(1)*cpv

     h(1)=tconv(1)*cpn(1)

     lv(1)=lv0-cpvmcl*(tconv(1)-273.15)

     hm(1)=lv(1)*qconv(1)

     tv(1)=tconv(1)*(1.+qconv(1)*epsi-qconv(1))

     ahmin=1.0e12

     ihmin=nl

     DO i=2,nl+1

       tvx=tconv(i)*(1.+qconv(i)*epsi-qconv(i))

       tvy=tconv(i-1)*(1.+qconv(i-1)*epsi-qconv(i-1))

       gz(i)=gz(i-1)+0.5*rd*(tvx+tvy)*(pconv_hpa(i-1)-pconv_hpa(i))/ &

            phconv_hpa(i)

       cpn(i)=cpd*(1.-qconv(i))+cpv*qconv(i)

       h(i)=tconv(i)*cpn(i)+gz(i)

       lv(i)=lv0-cpvmcl*(tconv(i)-273.15)

       hm(i)=(cpd*(1.-qconv(i))+cl*qconv(i))*(tconv(i)-tconv(1))+ &

            lv(i)*qconv(i)+gz(i)

       tv(i)=tconv(i)*(1.+qconv(i)*epsi-qconv(i))

   !

   !  ***  Find level of minimum moist static energy    ***

   !

       IF(i.GE.minorig.AND.hm(i).LT.ahmin.AND.hm(i).LT.hm(i-1))THEN

        ahmin=hm(i)

        ihmin=i

       END IF

     END DO

     ihmin=min(ihmin, nl-1)

   !

   !  ***     Find that model level below the level of minimum moist       ***

   !  ***  static energy that has the maximum value of moist static energy ***

   !

     ahmax=0.0

   !  ***  bug fixed: need to assign an initial value to NK

   !  HSO, 05.08.2009

     nk=minorig

     DO i=minorig,ihmin

      IF(hm(i).GT.ahmax)THEN

       nk=i

       ahmax=hm(i)

      END IF

     END DO

   !

   !  ***  CHECK WHETHER PARCEL LEVEL TEMPERATURE AND SPECIFIC HUMIDITY   ***

   !  ***                          ARE REASONABLE                         ***

   !  ***      Skip convection if HM increases monotonically upward       ***

   !

     IF(tconv(nk).LT.250.0.OR.qconv(nk).LE.0.0.OR.ihmin.EQ.(nl-1)) &

          THEN

      iflag=0

      cbmf=0.0

      RETURN

     END IF

   !

   !   ***  CALCULATE LIFTED CONDENSATION LEVEL OF AIR AT PARCEL ORIGIN LEVEL ***

   !   ***       (WITHIN 0.2% OF FORMULA OF BOLTON, MON. WEA. REV.,1980)      ***

   !

     rh=qconv(nk)/qsconv(nk)

     chi=tconv(nk)/(1669.0-122.0*rh-tconv(nk))

     plcl=pconv_hpa(nk)*(rh**chi)

     IF(plcl.LT.200.0.OR.plcl.GE.2000.0)THEN

      iflag=2

      cbmf=0.0

      RETURN

     END IF

   !

   !   ***  CALCULATE FIRST LEVEL ABOVE LCL (=ICB)  ***

   !

     icb=nl-1

     DO i=nk+1,nl

      IF(pconv_hpa(i).LT.plcl)THEN

       icb=min(icb,i)

      END IF

     END DO

     IF(icb.GE.(nl-1))THEN

      iflag=3

      cbmf=0.0

      RETURN

     END IF

   !

   !   *** FIND TEMPERATURE UP THROUGH ICB AND TEST FOR INSTABILITY           ***

   !

   !   *** SUBROUTINE TLIFT CALCULATES PART OF THE LIFTED PARCEL VIRTUAL      ***

   !   ***  TEMPERATURE, THE ACTUAL TEMPERATURE AND THE ADIABATIC             ***

   !   ***                   LIQUID WATER CONTENT                             ***

   !

     CALL tlift(gz,icb,nk,tvp,tp,clw,nd,nl,1)

     DO i=nk,icb

      tvp(i)=tvp(i)-tp(i)*qconv(nk)

     END DO

   !

   !   ***  If there was no convection at last time step and parcel    ***

   !   ***       is stable at ICB then skip rest of calculation        ***

   !

     IF(cbmf.EQ.0.0.AND.tvp(icb).LE.(tv(icb)-dtmax))THEN

      iflag=0

      RETURN

     END IF

   !

   !   ***  IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY ***

   !

     IF(iflag.NE.4)iflag=1

   !

   !   ***  FIND THE REST OF THE LIFTED PARCEL TEMPERATURES          ***

   !

     CALL tlift(gz,icb,nk,tvp,tp,clw,nd,nl,2)

   !

   !   ***  SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF   ***

   !   ***          PRECIPITATION FALLING OUTSIDE OF CLOUD           ***

   !   ***      THESE MAY BE FUNCTIONS OF TP(I), PCONV_HPA(I) AND CLW(I)     ***

   !

     DO i=1,nk

      ep(i)=0.0

      sigp(i)=sigs

     END DO

     DO i=nk+1,nl

      tca=tp(i)-273.15

      IF(tca.GE.0.0)THEN

       elacrit=elcrit

      ELSE

       elacrit=elcrit*(1.0-tca/tlcrit)

      END IF

      elacrit=max(elacrit,0.0)

      epmax=0.999

      ep(i)=epmax*(1.0-elacrit/max(clw(i),1.0e-8))

      ep(i)=max(ep(i),0.0)

      ep(i)=min(ep(i),epmax)

      sigp(i)=sigs

     END DO

   !

   !   ***       CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL     ***

   !   ***                    VIRTUAL TEMPERATURE                    ***

   !

     DO i=icb+1,nl

      tvp(i)=tvp(i)-tp(i)*qconv(nk)

     END DO

     tvp(nl+1)=tvp(nl)-(gz(nl+1)-gz(nl))/cpd

   !

   !   ***        NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS       ***

   !

     DO i=1,nl+1

      hp(i)=h(i)

      nent(i)=0

      water(i)=0.0

      evap(i)=0.0

      wt(i)=omtsnow

      lvcp(i)=lv(i)/cpn(i)

      DO j=1,nl+1

       qent(i,j)=qconv(j)

       elij(i,j)=0.0

       sij(i,j)=0.0

      END DO

     END DO

     qp(1)=qconv(1)

     DO i=2,nl+1

      qp(i)=qconv(i-1)

     END DO

   !

   !  ***  FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S      ***

   !  ***          HIGHEST LEVEL OF NEUTRAL BUOYANCY                 ***

   !  ***     AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)           ***

   !

     cape=0.0

     capem=0.0

     inb=icb+1

     inb1=inb

     byp=0.0

     DO i=icb+1,nl-1

      by=(tvp(i)-tv(i))*(phconv_hpa(i)-phconv_hpa(i+1))/pconv_hpa(i)

      cape=cape+by

      IF(by.GE.0.0)inb1=i+1

      IF(cape.GT.0.0)THEN

       inb=i+1

       byp=(tvp(i+1)-tv(i+1))*(phconv_hpa(i+1)-phconv_hpa(i+2))/ &

            pconv_hpa(i+1)

       capem=cape

      END IF

     END DO

     inb=max(inb,inb1)

     cape=capem+byp

     defrac=capem-cape

     defrac=max(defrac,0.001)

     frac=-cape/defrac

     frac=min(frac,1.0)

     frac=max(frac,0.0)

   !

   !   ***   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL   ***

   !

     DO i=icb,inb

      hp(i)=h(nk)+(lv(i)+(cpd-cpv)*tconv(i))*ep(i)*clw(i)

     END DO

   !

   !   ***  CALCULATE CLOUD BASE MASS FLUX AND RATES OF MIXING, M(I),  ***

   !   ***                   AT EACH MODEL LEVEL                       ***

   !

     dbosum=0.0

   !

   !   ***     INTERPOLATE DIFFERENCE BETWEEN LIFTED PARCEL AND      ***

   !   ***  ENVIRONMENTAL TEMPERATURES TO LIFTED CONDENSATION LEVEL  ***

   !

     tvpplcl=tvp(icb-1)-rd*tvp(icb-1)*(pconv_hpa(icb-1)-plcl)/ &

          (cpn(icb-1)*pconv_hpa(icb-1))

     tvaplcl=tv(icb)+(tvp(icb)-tvp(icb+1))*(plcl-pconv_hpa(icb))/ &

          (pconv_hpa(icb)-pconv_hpa(icb+1))

     dtpbl=0.0

     DO i=nk,icb-1

      dtpbl=dtpbl+(tvp(i)-tv(i))*(phconv_hpa(i)-phconv_hpa(i+1))

     END DO

     dtpbl=dtpbl/(phconv_hpa(nk)-phconv_hpa(icb))

     dtmin=tvpplcl-tvaplcl+dtmax+dtpbl

     dtma=dtmin

   !

   !   ***  ADJUST CLOUD BASE MASS FLUX   ***

   !

   cbmfold=cbmf

   ! *** C. Forster: adjustment of CBMF is not allowed to depend on FLEXPART timestep

   delt0=delt/3.

   damps=damp*delt/delt0

   cbmf=(1.-damps)*cbmf+0.1*alpha*dtma

   cbmf=max(cbmf,0.0)

   !

   !   *** If cloud base mass flux is zero, skip rest of calculation  ***

   !

   IF(cbmf.EQ.0.0.AND.cbmfold.EQ.0.0)THEN

    RETURN

   END IF


   !

   !   ***   CALCULATE RATES OF MIXING,  M(I)   ***

   !

   m(icb)=0.0

   DO i=icb+1,inb

    k=min(i,inb1)

    dbo=abs(tv(k)-tvp(k))+ &

         entp*0.02*(phconv_hpa(k)-phconv_hpa(k+1))

    dbosum=dbosum+dbo

    m(i)=cbmf*dbo

   END DO

   DO i=icb+1,inb

    m(i)=m(i)/dbosum

   END DO

   !

   !   ***  CALCULATE ENTRAINED AIR MASS FLUX (MENT), TOTAL WATER MIXING  ***

   !   ***     RATIO (QENT), TOTAL CONDENSED WATER (ELIJ), AND MIXING     ***

   !   ***                        FRACTION (SIJ)                          ***

   !

     DO i=icb+1,inb

      qti=qconv(nk)-ep(i)*clw(i)

      DO j=icb,inb

       bf2=1.+lv(j)*lv(j)*qsconv(j)/(rv*tconv(j)*tconv(j)*cpd)

       anum=h(j)-hp(i)+(cpv-cpd)*tconv(j)*(qti-qconv(j))

       denom=h(i)-hp(i)+(cpd-cpv)*(qconv(i)-qti)*tconv(j)

       dei=denom

       IF(abs(dei).LT.0.01)dei=0.01

       sij(i,j)=anum/dei

       sij(i,i)=1.0

       altem=sij(i,j)*qconv(i)+(1.-sij(i,j))*qti-qsconv(j)

       altem=altem/bf2

       cwat=clw(j)*(1.-ep(j))

       stemp=sij(i,j)

       IF((stemp.LT.0.0.OR.stemp.GT.1.0.OR. &

            altem.GT.cwat).AND.j.GT.i)THEN

        anum=anum-lv(j)*(qti-qsconv(j)-cwat*bf2)

        denom=denom+lv(j)*(qconv(i)-qti)

        IF(abs(denom).LT.0.01)denom=0.01

        sij(i,j)=anum/denom

        altem=sij(i,j)*qconv(i)+(1.-sij(i,j))*qti-qsconv(j)

        altem=altem-(bf2-1.)*cwat

       END IF

       IF(sij(i,j).GT.0.0.AND.sij(i,j).LT.0.9)THEN

        qent(i,j)=sij(i,j)*qconv(i)+(1.-sij(i,j))*qti

        elij(i,j)=altem

        elij(i,j)=max(0.0,elij(i,j))

        ment(i,j)=m(i)/(1.-sij(i,j))

        nent(i)=nent(i)+1

       END IF

       sij(i,j)=max(0.0,sij(i,j))

       sij(i,j)=min(1.0,sij(i,j))

      END DO

   !

   !   ***   IF NO AIR CAN ENTRAIN AT LEVEL I ASSUME THAT UPDRAFT DETRAINS  ***

   !   ***   AT THAT LEVEL AND CALCULATE DETRAINED AIR FLUX AND PROPERTIES  ***

   !

      IF(nent(i).EQ.0)THEN

       ment(i,i)=m(i)

       qent(i,i)=qconv(nk)-ep(i)*clw(i)

       elij(i,i)=clw(i)

       sij(i,i)=1.0

      END IF

     END DO

     sij(inb,inb)=1.0

   !

   !   ***  NORMALIZE ENTRAINED AIR MASS FLUXES TO REPRESENT EQUAL  ***

   !   ***              PROBABILITIES OF MIXING                     ***

   !

     DO i=icb+1,inb

     IF(nent(i).NE.0)THEN

      qp1=qconv(nk)-ep(i)*clw(i)

      anum=h(i)-hp(i)-lv(i)*(qp1-qsconv(i))

      denom=h(i)-hp(i)+lv(i)*(qconv(i)-qp1)

      IF(abs(denom).LT.0.01)denom=0.01

      scrit=anum/denom

      alt=qp1-qsconv(i)+scrit*(qconv(i)-qp1)

      IF(alt.LT.0.0)scrit=1.0

      scrit=max(scrit,0.0)

      asij=0.0

      smin=1.0

      DO j=icb,inb

       IF(sij(i,j).GT.0.0.AND.sij(i,j).LT.0.9)THEN

        IF(j.GT.i)THEN

         smid=min(sij(i,j),scrit)

         sjmax=smid

         sjmin=smid

         IF(smid.LT.smin.AND.sij(i,j+1).LT.smid)THEN

          smin=smid

          sjmax=min(sij(i,j+1),sij(i,j),scrit)

          sjmin=max(sij(i,j-1),sij(i,j))

          sjmin=min(sjmin,scrit)

         END IF

        ELSE

         sjmax=max(sij(i,j+1),scrit)

         smid=max(sij(i,j),scrit)

         sjmin=0.0

         IF(j.GT.1)sjmin=sij(i,j-1)

         sjmin=max(sjmin,scrit)

        END IF

        delp=abs(sjmax-smid)

        delm=abs(sjmin-smid)

        asij=asij+(delp+delm)*(phconv_hpa(j)-phconv_hpa(j+1))

        ment(i,j)=ment(i,j)*(delp+delm)* &

             (phconv_hpa(j)-phconv_hpa(j+1))

       END IF

      END DO

      asij=max(1.0e-21,asij)

      asij=1.0/asij

      DO j=icb,inb

       ment(i,j)=ment(i,j)*asij

      END DO

      bsum=0.0

      DO j=icb,inb

       bsum=bsum+ment(i,j)

      END DO

      IF(bsum.LT.1.0e-18)THEN

       nent(i)=0

       ment(i,i)=m(i)

       qent(i,i)=qconv(nk)-ep(i)*clw(i)

       elij(i,i)=clw(i)

       sij(i,i)=1.0

      END IF

     END IF

     END DO

   !

   !   ***  CHECK WHETHER EP(INB)=0, IF SO, SKIP PRECIPITATING    ***

   !   ***             DOWNDRAFT CALCULATION                      ***

   !

     IF(ep(inb).LT.0.0001)goto 405

   !

   !   ***  INTEGRATE LIQUID WATER EQUATION TO FIND CONDENSED WATER   ***

   !   ***                AND CONDENSED WATER FLUX                    ***

   !

     jtt=2

   !

   !    ***                    BEGIN DOWNDRAFT LOOP                    ***

   !

     DO i=inb,1,-1

   !

   !    ***              CALCULATE DETRAINED PRECIPITATION             ***

   !

     wdtrain=g*ep(i)*m(i)*clw(i)

     IF(i.GT.1)THEN

      DO j=1,i-1

      awat=elij(j,i)-(1.-ep(i))*clw(i)

      awat=max(0.0,awat)

        wdtrain=wdtrain+g*awat*ment(j,i)

      END DO

     END IF

   !

   !    ***    FIND RAIN WATER AND EVAPORATION USING PROVISIONAL   ***

   !    ***              ESTIMATES OF QP(I)AND QP(I-1)             ***

   !

   !

   !  ***  Value of terminal velocity and coefficient of evaporation for snow   ***

   !

     coeff=coeffs

     wt(i)=omtsnow

   !

   !  ***  Value of terminal velocity and coefficient of evaporation for rain   ***

   !

     IF(tconv(i).GT.273.0)THEN

      coeff=coeffr

      wt(i)=omtrain

     END IF

     qsm=0.5*(qconv(i)+qp(i+1))

     afac=coeff*phconv_hpa(i)*(qsconv(i)-qsm)/ &

          (1.0e4+2.0e3*phconv_hpa(i)*qsconv(i))

     afac=max(afac,0.0)

     sigt=sigp(i)

     sigt=max(0.0,sigt)

     sigt=min(1.0,sigt)

     b6=100.*(phconv_hpa(i)-phconv_hpa(i+1))*sigt*afac/wt(i)

     c6=(water(i+1)*wt(i+1)+wdtrain/sigd)/wt(i)

     revap=0.5*(-b6+sqrt(b6*b6+4.*c6))

     evap(i)=sigt*afac*revap

     water(i)=revap*revap

   !

   !    ***  CALCULATE PRECIPITATING DOWNDRAFT MASS FLUX UNDER     ***

   !    ***              HYDROSTATIC APPROXIMATION                 ***

   !

     IF(i.EQ.1)goto 360

     dhdp=(h(i)-h(i-1))/(pconv_hpa(i-1)-pconv_hpa(i))

     dhdp=max(dhdp,10.0)

     mp(i)=100.*ginv*lv(i)*sigd*evap(i)/dhdp

     mp(i)=max(mp(i),0.0)

   !

   !   ***   ADD SMALL AMOUNT OF INERTIA TO DOWNDRAFT              ***

   !

     fac=20.0/(phconv_hpa(i-1)-phconv_hpa(i))

     mp(i)=(fac*mp(i+1)+mp(i))/(1.+fac)

   !

   !    ***      FORCE MP TO DECREASE LINEARLY TO ZERO                 ***

   !    ***      BETWEEN ABOUT 950 MB AND THE SURFACE                  ***

   !

       IF(pconv_hpa(i).GT.(0.949*pconv_hpa(1)))THEN

        jtt=max(jtt,i)

        mp(i)=mp(jtt)*(pconv_hpa(1)-pconv_hpa(i))/(pconv_hpa(1)- &

             pconv_hpa(jtt))

       END IF

   360   CONTINUE

   !

   !    ***       FIND MIXING RATIO OF PRECIPITATING DOWNDRAFT     ***

   !

     IF(i.EQ.inb)goto 400

     IF(i.EQ.1)THEN

      qstm=qsconv(1)

     ELSE

      qstm=qsconv(i-1)

     END IF

     IF(mp(i).GT.mp(i+1))THEN

       rat=mp(i+1)/mp(i)

       qp(i)=qp(i+1)*rat+qconv(i)*(1.0-rat)+100.*ginv* &

            sigd*(phconv_hpa(i)-phconv_hpa(i+1))*(evap(i)/mp(i))

      ELSE

       IF(mp(i+1).GT.0.0)THEN

         qp(i)=(gz(i+1)-gz(i)+qp(i+1)*(lv(i+1)+tconv(i+1)*(cl-cpd))+ &

              cpd*(tconv(i+1)-tconv(i)))/(lv(i)+tconv(i)*(cl-cpd))

       END IF

     END IF

     qp(i)=min(qp(i),qstm)

     qp(i)=max(qp(i),0.0)

 400 CONTINUE

     END DO

   !

   !   ***  CALCULATE SURFACE PRECIPITATION IN MM/DAY     ***

   !

     precip=precip+wt(1)*sigd*water(1)*3600.*24000./(rowl*g)

   !

   405   CONTINUE

   !

   !   ***  CALCULATE DOWNDRAFT VELOCITY SCALE AND SURFACE TEMPERATURE AND  ***

   !   ***                    WATER VAPOR FLUCTUATIONS                      ***

   !

   wd=beta*abs(mp(icb))*0.01*rd*tconv(icb)/(sigd*pconv_hpa(icb))

   qprime=0.5*(qp(1)-qconv(1))

   tprime=lv0*qprime/cpd

   !

   !   ***  CALCULATE TENDENCIES OF LOWEST LEVEL POTENTIAL TEMPERATURE  ***

   !   ***                      AND MIXING RATIO                        ***

   !

     dpinv=0.01/(phconv_hpa(1)-phconv_hpa(2))

     am=0.0

     IF(nk.EQ.1)THEN

      DO k=2,inb

        am=am+m(k)

      END DO

     END IF

   ! save saturated upward mass flux for first level

     fup(1)=am

     IF((2.*g*dpinv*am).GE.delti)iflag=4

     ft(1)=ft(1)+g*dpinv*am*(tconv(2)-tconv(1)+(gz(2)-gz(1))/cpn(1))

     ft(1)=ft(1)-lvcp(1)*sigd*evap(1)

     ft(1)=ft(1)+sigd*wt(2)*(cl-cpd)*water(2)*(tconv(2)- &

          tconv(1))*dpinv/cpn(1)

     fq(1)=fq(1)+g*mp(2)*(qp(2)-qconv(1))* &

          dpinv+sigd*evap(1)

     fq(1)=fq(1)+g*am*(qconv(2)-qconv(1))*dpinv

     DO j=2,inb

      fq(1)=fq(1)+g*dpinv*ment(j,1)*(qent(j,1)-qconv(1))

     END DO

   !

   !   ***  CALCULATE TENDENCIES OF POTENTIAL TEMPERATURE AND MIXING RATIO  ***

   !   ***               AT LEVELS ABOVE THE LOWEST LEVEL                   ***

   !

   !   ***  FIRST FIND THE NET SATURATED UPDRAFT AND DOWNDRAFT MASS FLUXES  ***

   !   ***                      THROUGH EACH LEVEL                          ***

   !

     DO i=2,inb

     dpinv=0.01/(phconv_hpa(i)-phconv_hpa(i+1))

     cpinv=1.0/cpn(i)

     amp1=0.0

     ad=0.0

     IF(i.GE.nk)THEN

      DO k=i+1,inb+1

        amp1=amp1+m(k)

      END DO

     END IF

     DO k=1,i

     DO j=i+1,inb+1

      amp1=amp1+ment(k,j)

     END DO

     END DO

   ! save saturated upward mass flux

     fup(i)=amp1

     IF((2.*g*dpinv*amp1).GE.delti)iflag=4

     DO k=1,i-1

     DO j=i,inb

      ad=ad+ment(j,k)

     END DO

     END DO

   ! save saturated downward mass flux

     fdown(i)=ad

     ft(i)=ft(i)+g*dpinv*(amp1*(tconv(i+1)-tconv(i)+(gz(i+1)-gz(i))* &

          cpinv)-ad*(tconv(i)-tconv(i-1)+(gz(i)-gz(i-1))*cpinv)) &

          -sigd*lvcp(i)*evap(i)

     ft(i)=ft(i)+g*dpinv*ment(i,i)*(hp(i)-h(i)+ &

          tconv(i)*(cpv-cpd)*(qconv(i)-qent(i,i)))*cpinv

     ft(i)=ft(i)+sigd*wt(i+1)*(cl-cpd)*water(i+1)* &

          (tconv(i+1)-tconv(i))*dpinv*cpinv

     fq(i)=fq(i)+g*dpinv*(amp1*(qconv(i+1)-qconv(i))- &

          ad*(qconv(i)-qconv(i-1)))

     DO k=1,i-1

      awat=elij(k,i)-(1.-ep(i))*clw(i)

      awat=max(awat,0.0)

      fq(i)=fq(i)+g*dpinv*ment(k,i)*(qent(k,i)-awat-qconv(i))

     END DO

     DO k=i,inb

      fq(i)=fq(i)+g*dpinv*ment(k,i)*(qent(k,i)-qconv(i))

     END DO

     fq(i)=fq(i)+sigd*evap(i)+g*(mp(i+1)* &

          (qp(i+1)-qconv(i))-mp(i)*(qp(i)-qconv(i-1)))*dpinv

     END DO

   !

   !   *** Adjust tendencies at top of convection layer to reflect  ***

   !   ***       actual position of the level zero CAPE             ***

   !

     fqold=fq(inb)

     fq(inb)=fq(inb)*(1.-frac)

     fq(inb-1)=fq(inb-1)+frac*fqold*((phconv_hpa(inb)- &

          phconv_hpa(inb+1))/ &

          (phconv_hpa(inb-1)-phconv_hpa(inb)))*lv(inb)/lv(inb-1)

     ftold=ft(inb)

     ft(inb)=ft(inb)*(1.-frac)

     ft(inb-1)=ft(inb-1)+frac*ftold*((phconv_hpa(inb)- &

          phconv_hpa(inb+1))/ &

          (phconv_hpa(inb-1)-phconv_hpa(inb)))*cpn(inb)/cpn(inb-1)

   !

   !   ***   Very slightly adjust tendencies to force exact   ***

   !   ***     enthalpy, momentum and tracer conservation     ***

   !

     ents=0.0

     DO i=1,inb

      ents=ents+(cpn(i)*ft(i)+lv(i)*fq(i))* &

           (phconv_hpa(i)-phconv_hpa(i+1))

     END DO

     ents=ents/(phconv_hpa(1)-phconv_hpa(inb+1))

     DO i=1,inb

      ft(i)=ft(i)-ents/cpn(i)

     END DO


   ! ************************************************

   ! **** DETERMINE MASS DISPLACEMENT MATRIX

   ! ***** AND COMPENSATING SUBSIDENCE

   ! ************************************************


   ! mass displacement matrix due to saturated up-and downdrafts

   ! inside the cloud and determine compensating subsidence

   ! FUP(I) (saturated updrafts), FDOWN(I) (saturated downdrafts) are assumed to be

   ! balanced by  compensating subsidence (SUB(I))

   ! FDOWN(I) and SUB(I) defined positive downwards


   ! NCONVTOP IS THE TOP LEVEL AT WHICH CONVECTIVE MASS FLUXES ARE DIAGNOSED

   ! EPSILON IS A SMALL NUMBER


    sub(1)=0.

    nconvtop=1

    do i=1,inb+1

    do j=1,inb+1

     if (j.eq.nk) then

      fmass(j,i)=fmass(j,i)+m(i)

     endif

      fmass(j,i)=fmass(j,i)+ment(j,i)

      IF (fmass(j,i).GT.epsilon) nconvtop=max(nconvtop,i,j)

    end do

    if (i.gt.1) then

     sub(i)=fup(i-1)-fdown(i)

    endif

    end do

    nconvtop=nconvtop+1


     RETURN

   !

 END SUBROUTINE convect

 !

 ! ---------------------------------------------------------------------------

 !

 SUBROUTINE tlift(GZ,ICB,NK,TVP,TPK,CLW,ND,NL,KK)

   !

   !-cv

   use par_mod

   use conv_mod


   implicit none

   !-cv

   !====>Begin Module TLIFT      File convect.f      Undeclared variables

   !

   !Argument variables

   !

   integer :: icb, kk, nd, nk, nl

   !

   !Local variables

   !

   integer :: i, j, nsb, nst

   !

   real :: ah0, ahg, alv, cpinv, cpp, denom

   real :: es, qg, rg, s, tc, tg

   !

   !====>End Module   TLIFT      File convect.f


     REAL :: gz(nd),tpk(nd),clw(nd)

     REAL :: tvp(nd)

   !

   !   ***   ASSIGN VALUES OF THERMODYNAMIC CONSTANTS     ***

   !

     REAL,PARAMETER :: cpd=1005.7

     REAL,PARAMETER :: cpv=1870.0

     REAL,PARAMETER :: cl=2500.0

     REAL,PARAMETER :: rv=461.5

     REAL,PARAMETER :: rd=287.04

     REAL,PARAMETER :: lv0=2.501e6

   !

     REAL,PARAMETER :: cpvmcl=cl-cpv

     REAL,PARAMETER :: eps0=rd/rv

     REAL,PARAMETER :: epsi=1./eps0

   !

   !   ***  CALCULATE CERTAIN PARCEL QUANTITIES, INCLUDING STATIC ENERGY   ***

   !

     ah0=(cpd*(1.-qconv(nk))+cl*qconv(nk))*tconv(nk)+qconv(nk)* &

          (lv0-cpvmcl*( &

          tconv(nk)-273.15))+gz(nk)

     cpp=cpd*(1.-qconv(nk))+qconv(nk)*cpv

     cpinv=1./cpp

   !

     IF(kk.EQ.1)THEN

   !

   !   ***   CALCULATE LIFTED PARCEL QUANTITIES BELOW CLOUD BASE   ***

   !

     DO i=1,icb-1

      clw(i)=0.0

     END DO

     DO i=nk,icb-1

      tpk(i)=tconv(nk)-(gz(i)-gz(nk))*cpinv

      tvp(i)=tpk(i)*(1.+qconv(nk)*epsi)

     END DO

     END IF

   !

   !    ***  FIND LIFTED PARCEL QUANTITIES ABOVE CLOUD BASE    ***

   !

     nst=icb

     nsb=icb

     IF(kk.EQ.2)THEN

      nst=nl

      nsb=icb+1

     END IF

     DO i=nsb,nst

      tg=tconv(i)

      qg=qsconv(i)

      alv=lv0-cpvmcl*(tconv(i)-273.15)

      DO j=1,2

       s=cpd+alv*alv*qg/(rv*tconv(i)*tconv(i))

       s=1./s

       ahg=cpd*tg+(cl-cpd)*qconv(nk)*tconv(i)+alv*qg+gz(i)

       tg=tg+s*(ah0-ahg)

       tg=max(tg,35.0)

       tc=tg-273.15

       denom=243.5+tc

       IF(tc.GE.0.0)THEN

        es=6.112*exp(17.67*tc/denom)

       ELSE

        es=exp(23.33086-6111.72784/tg+0.15215*log(tg))

       END IF

       qg=eps0*es/(pconv_hpa(i)-es*(1.-eps0))

      END DO

      alv=lv0-cpvmcl*(tconv(i)-273.15)

      tpk(i)=(ah0-(cl-cpd)*qconv(nk)*tconv(i)-gz(i)-alv*qg)/cpd

      clw(i)=qconv(nk)-qg

      clw(i)=max(0.0,clw(i))

      rg=qg/(1.-qconv(nk))

      tvp(i)=tpk(i)*(1.+rg*epsi)

     END DO

     RETURN

 END SUBROUTINE tlift

convect
subroutine convect(ND, NL, DELT, IFLAG, PRECIP, WD, TPRIME, QPRIME, CBMF)
Definition: convect43c.f90:29

par_mod
Definition: par_mod.F90:34

tlift
subroutine tlift(GZ, ICB, NK, TVP, TPK, CLW, ND, NL, KK)
Definition: convect43c.f90:1015

conv_mod
Definition: conv_mod.f90:11