ecearth3/sources/tm5mp/base/phys_convec_tmpp.F90

!
! Subgrid stuff from TMPP.
!
! Copied from:
!   TMPP/SOURCE/tmpp_sub_subg.f90
!   TMPP/SOURCE/tmpp_sub_work.f90
!   TMPP/SOURCE/tmpp_geometry.f90
!   TMPP/SOURCE/phys_geopot.f90
!   

module phys_convec_tmpp

  implicit none
  
  ! --- in/out -------------------------------

  private
  
  public   ::  subscal
  public   ::  subscal_2d

  public   ::  mid2bound_uv
  public   ::  mid2bound_w
  public   ::  mid2bound_t
  public   ::  mid2bound_q
  public   ::  mid2bound_z
  public   ::  mid2bound_p
  
  public   ::  phlev_ec_rev

  public   ::  geopot

contains


  ! ==========================================================
  ! ===
  ! ===  TMPP/SOURCE/tmpp_sub_subg.f90
  ! ===
  ! ==========================================================
  

  !-----------------------------------------------------------------------
  ! calculate subscal parameters at a given horizontal location.
  ! Three different pressure levels are used in this subroutine:
  ! ppress (boundaries ECMWF levels)
  ! ppresg (boundaries TM3)
  ! zpl (centers ECMWF levels)
  ! zplg (centers TM3)
  !
  !-----------------------------------------------------------------------

  !subroutine subscal(pz,ppress,pu,pv,pw,pt,pq,pqac,pqam,pevap, &
  !                   pentug,pdetug,pentdg,pdetdg,pdkg,pzg)

  subroutine subscal_2d( np, lme, at, bt, &
                     pz,ppress,pw,pt,pq,pqac,pqam,pevap, &
                     pentug,pdetug,pentdg,pdetdg)

    ! --- in/out -----------------------------------
    
    integer, intent(in)  ::  np
    integer, intent(in)  ::  lme
    real, intent(in)     ::  at(0:lme), bt(0:lme)
    real, intent(in)     ::  pz(np,0:lme)
    real, intent(in)     ::  ppress(np,0:lme)
    real, intent(in)     ::  pw(np,0:lme)
    real, intent(in)     ::  pt(np,0:lme)
    real, intent(in)     ::  pq(np,0:lme)
    real, intent(in)     ::  pqac(np,0:lme)
    real, intent(in)     ::  pqam(np,0:lme)
    real, intent(in)     ::  pevap(np)
    real, intent(out)    ::  pentug(np,lme)
    real, intent(out)    ::  pdetug(np,lme)
    real, intent(out)    ::  pentdg(np,lme)
    real, intent(out)    ::  pdetdg(np,lme)
    
    ! --- local -------------------------------------
    
    integer     :: i
    
    ! --- begin ------------------------------------

    !   18 Apr 2012 - P. Le Sager - Commented openMP loop, since it gives
    !   wrong results: large differences with serial case results in eu/ed/du
    !   /dd fields
    
!testPLS    !$OMP PARALLEL &
!testPLS    !$OMP       default ( none ) &
!testPLS    !$OMP       shared  ( np, lme, at, bt ) &
!testPLS    !$OMP       shared  ( pz, ppress, pw, pt, pq ) &
!testPLS    !$OMP       shared  ( pqac, pqam, pevap ) &
!testPLS    !$OMP       shared  ( pentug, pdetug, pentdg, pdetdg ) &
!testPLS    !$OMP       private ( i )
!testPLS    !$OMP   DO
    do i = 1, np
      call subscal( lme, at, bt, &
                     pz(i,:), ppress(i,:), pw(i,:), pt(i,:), pq(i,:), &
                     pqac(i,:), pqam(i,:), pevap(i), &
                     pentug(i,:), pdetug(i,:), pentdg(i,:), pdetdg(i,:) )
    end do
!testPLS    !$OMP   END DO
!testPLS    !$OMP END PARALLEL
    
  end subroutine subscal_2d
 
  ! *
  
  subroutine subscal( lme, at, bt, &
                     pz,ppress,pw,pt,pq,pqac,pqam,pevap, &
                     pentug,pdetug,pentdg,pdetdg) !,pdkg,pzg)

    ! >>> adhoc: not from ECMWF to TM levels >>>>>>>>>>>>>>>>>>
    !     (set tm stuff to ec stuff)
    !use Geometry, only : lm
    !use Geometry, only : at => a_tm, bt => b_tm
!    use Geometry, only : lm => lme
!    use Geometry, only : at => a_ec, bt => b_ec
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

!    use Geometry, only : lme

    use Num, only : IntervalQuad_Lin
    use Num, only : interp_muherm

    ! --- const ------------------------------------

!    integer, parameter  ::  jplm = lme

    !integer, parameter  ::  jplmg = lm
    !integer, parameter  ::  jplmgp1 = lm+1
!    integer, parameter  ::  jplmg = lme
!    integer, parameter  ::  jplmgp1 = lme+1

    real, parameter   ::  ppg = 9.80665

    ! --- in/out ------------------------------------
    
    integer              ::  status

    integer, intent(in)  ::  lme
    real, intent(in)     ::  at(0:lme), bt(0:lme)

    real, intent(in)   :: pz(0:lme)
    real, intent(in)   :: ppress(0:lme)
    real, intent(in)   :: pw(0:lme)
    real, intent(in)   :: pt(0:lme)
    real, intent(in)   :: pq(0:lme)
    real, intent(in)   :: pqac(0:lme)
    real, intent(in)   :: pqam(0:lme)
    real, intent(in)   :: pevap

    !real, intent(in)   :: pu(0:lme)
    !real, intent(in)   :: pv(0:lme)

    ! entrainment, detrainment rates, diffusion coefficient

    real, intent(out)  :: pentug(lme)
    real, intent(out)  :: pdetug(lme)
    real, intent(out)  :: pentdg(lme)
    real, intent(out)  :: pdetdg(lme)

!    real, intent(out)  :: pdkg(lme)

!    real, intent(out)  :: pzg(0:lme)

    ! --- local ------------------------------------

    integer :: jl,ilast

    real :: sumup, sumdown


    ! output variables on TM vertical grid

    real ::  zamu(0:lme),zptc(0:lme),zpqc(0:lme),zplc(0:lme), &
         zpqp(0:lme),ppresg(0:lme)
!    real ::  zpdk(0:lme)

    real ::  zpam(lme),zped(lme),zpdd(lme),zpeu(lme),zpdu(lme)
    real ::  zpl(lme),zplg(lme)
    real ::  zam

    integer ::  lct

    ! --- begin -------------------------------------------

    ! calculate convection by clouds

    call cloud( lme, pz,ppress,pt,pq,pw,pqac,pqam,pevap, &
               zamu,zptc,zpqc,zplc,zpqp, &
               zpam,zpeu,zpdu,zped,zpdd)


!    ! calculate turbulent diffusion coefficient
!
!    call louis(pz,pt,pu,pv,zpdk)

    ! generate pressures on TM3 level boundaries and centers respectively

    do jl=0,lme
      ppresg(jl)=at(jl)+ppress(lme)*bt(jl)
    end do

    do jl=1,lme
      zplg(jl)=0.5*(ppresg(jl)+ppresg(jl-1))
    end do

    ! --- Vertical averaging of variables defined on level centers
    
    ! * interpolate entrainment/detrainment rates on TM3 vertical coordinate:

    ! Dirk Olivie, 12 May 2004
    ! No interpolation needed if lme = lme

    if ( lme .ne. lme ) then 

      do jl=1,lme
        zpl(jl)=0.5*(ppress(jl)+ppress(jl-1))
        zpeu(jl)=zpeu(jl)/(ppress(jl)-ppress(jl-1))
        zpdu(jl)=zpdu(jl)/(ppress(jl)-ppress(jl-1))
        zped(jl)=zped(jl)/(ppress(jl)-ppress(jl-1))
        zpdd(jl)=zpdd(jl)/(ppress(jl)-ppress(jl-1))
      end do
      ilast=1
      do jl=1,lme
        call IntervalQuad_Lin(zpl,zpeu,ppresg(jl-1),ppresg(jl),pentug(jl), ilast, status )
        if (status/=0) stop 'ERROR in subscal'
        call IntervalQuad_Lin(zpl,zpdu,ppresg(jl-1),ppresg(jl),pdetug(jl), ilast, status )
        if (status/=0) stop 'ERROR in subscal'
        call IntervalQuad_Lin(zpl,zped,ppresg(jl-1),ppresg(jl),pentdg(jl), ilast, status )
        if (status/=0) stop 'ERROR in subscal'
        call IntervalQuad_Lin(zpl,zpdd,ppresg(jl-1),ppresg(jl),pdetdg(jl), ilast, status )
        if (status/=0) stop 'ERROR in subscal'
      end do

    else
    
      pentug = zpeu
      pdetug = zpdu
      pentdg = zped
      pdetdg = zpdd
      
    endif

!    ! * Interpolate diffusion coefficients from upper boundaries of ECMWF layers to
!    !   upper boundaries of TM3 layers and 
!    !   reorganize the vertical index of the diffusion coefficient
!    !   pdkg(1) is diffusion coefficient at the top of layer 1 (TOA)
!    !   pdkg(lme) is diffusion coefficient at the top of layer lme
!
!    call interp_muherm( ppress, zpdk, ppresg, pdkg )
!    pdkg(1)=0.
!
!    ! * Interpolate geopotential height from ECMWF layer boundaries to TM3 layer 
!    !  boundaries 
!
!    call interp_muherm( ppress, pz, ppresg, pzg )
!    pzg(0)=pz(0)
!    pzg(lme)=pz(lme)

    !--------------------------------------------------------------------------
    ! correct massfluxes on TM3 grid to conserve mass
    !
    ! updraft

    ! correction level is uppermost level with nonzero detrainment
    lct=lme
    do jl=1,lme
      if ( pdetug(jl) > 0.0 ) then
        lct=jl
        exit
      endif
    end do

    zam=0.
    do jl=lme,1,-1
      zam=zam+pentug(jl)-pdetug(jl)
    end do

    pdetug(lct)=pdetug(lct)+zam
    ! downdraft

    ! correction level is lowermost level

    lct=lme
    zam=0.
    do jl=lme,1,-1
      zam=zam+pentdg(jl)-pdetdg(jl)
    end do


    pdetdg(lct)=pdetdg(lct)+zam
    !
    ! check conservation
    !
    sumup=0.
    sumdown=0.
    do jl=1,lme
      sumup=sumup+pentug(jl)-pdetug(jl)
      sumdown=sumdown+pentdg(jl)-pdetdg(jl)
    enddo
    if ( abs(sumup) > 1.0e-5 ) then
      write(*,*)' ERROR no massconserv in updraft'
      stop
    endif
    if ( abs(sumdown) > 1.0e-5 ) then
      write(*,*)' ERROR no massconserv in downdraft'
      stop
    endif

!!pentug = 0.0
!!pdetug = 0.0
!!pentdg = 0.0
!!pdetdg = 0.0

  end subroutine subscal


  ! =================================================
      

      !***********************************************************************
      !****     cloud   - routine to calculate cloud properties
      !
      ! m.heimann       mpi HH  Nov-12-1990
      !
      ! Purpose
      ! -------
      ! cloud calculates all properties associated with subgridscale
      ! cloud mixing, i.e. massflux in updraft and downdraft, entrainment
      ! and detrainment rates per level, and cloud properties: temperature,
      ! moisture and liquid water and precipitation rate on each level
      ! boundary.
      !
      !**       interface
      ! ---------
      !
      !       call cloud(pz,ppress,pt,pq,pw,pqac,pqam,
      !     .            pqtur,pamu,ptc,pqc,plc,pgp,
      !     .            pam,peu,pdu,ped,pdd,
      !     .            klc,klt)
      !
      ! input: (all variables are defined on level boundaries)
      !         pz      geopotential height (m)
      !         ppress  pressure        (Pa)
      !         pt      temperature     (K)
      !         pq      moisture        (kg water/kg air)
      !         pw      vertcal velocity (negative upward)   (Pa s**-1)
      !         pqac    div(q v)        (kg water/kg air s**-1)
      !         pqam    div(v)          (s**-1)
      !         pqtur   Fq surf-air     (kg water/m**2 s**-1)
      !
      ! output: variables defined on level boundaries:
      !         pamu    massflux                (kg m**-2 s**-1)
      !         ptc     updraft temperature     (K)
      !         pqc     updraft moisture        (kg water/kg air)
      !         plc     updraft liquid water    (kg water/kg air)
      !         pgp     updraft precipitation   (kg water m**-2 s**-1)
      !
      !         variables defined on each model level:
      !         pam     air mass                        (kg m**-2)
      !         peu     entrainment in updraft          (kg m**-2 s**-1)
      !         pdu     detrainment in updraft          (same)
      !         ped     entrainment in downdraft        (kg m**-2 s**-1)
      !         pdd     detrainment in downdraft        (same)
      !         klc     lowest level in cloud
      !         klt     first level above cloud
      !
      ! method
      ! ------
      !
      ! updraft cloud properties are calculated according to Tiedke scheme
      !
      ! externals
      ! ---------
      !
      ! needs functions qsat and dqsat_dt
      !
      ! references
      ! ----------
      !
      ! Tiedke, Mon. Wea. Rev., 117, 1779-1800, 1989.
      !
      ! revisions
      ! ---------
      !   26-jan-1995 , sr
      !
      !   including now cumulus downdraft and midlevel convection.
      !   Calculates corresponding entrainment and detrainment rates
      !
      !-------------------------------------------------------------------------

      subroutine cloud(lme, pz,ppress,pt,pq,pw,pqac,pqam, pqtur,&
                       pamu,ptc,pqc,plc,pgp, &
                       pam,peu,pdu,ped,pdd)
                       
        use phys_humidity, only : QSat, dQSat_dt

        ! --- flags -----------------------------------------

        ! parameter llcudo=.true.  : cumulus downdraft turned on
        ! parameter llcudo=.false. : cumulus downdraft turned off
        ! parameter llmilc=.true.  : midlevel convection turned on
        ! parameter llmilc=.false. : midlevel convection turned off

        logical       :: llcudo = .true.
        logical       :: llmilc = .true.

        ! --- const (dims) -------------------------------------
        
        !! vertical resolution (no of model layers)
        !integer   ::  jplm = lme
        !integer   ::  jplmm1 = jplm-1

        ! --- in/out ----------------------------------------

        integer, intent(in) ::  lme

        real, intent(in)  :: pz(0:lme)
        real, intent(in)  :: pt(0:lme)
        real, intent(in)  :: pq(0:lme)
        real, intent(in)  :: pw(0:lme)
        real, intent(in)  :: ppress(0:lme)
        real, intent(in)  :: pqac(0:lme),pqam(0:lme)
        real, intent(in)  :: pqtur

        real, intent(out) :: pam(lme),peu(lme),pdu(lme),ped(lme),pdd(lme)
        real, intent(out) :: pamu(0:lme),ptc(0:lme),pqc(0:lme),plc(0:lme),pgp(0:lme)
 
        ! --- const (phys) -------------------------------------

        ! physical constants

        real, parameter  ::  pprgasd = 287.05
        real, parameter  ::  pprgasv = 461.51
        real, parameter  ::  ppeps = pprgasd/pprgasv

        real, parameter  ::  ppg = 9.80665
        real, parameter  ::  ppcpd = 1005.46
        real, parameter  ::  ppalv = 2.5008e6
        real, parameter  ::  ppzeta = ppalv/ppcpd
        real, parameter  ::  ppvtcf = (1.0-ppeps)/ppeps

        ! * constants for turbulent entrainment and detrainment rates

        ! penetrative and midlevel convection
        real, parameter   ::  ppepsu = 1.e-4
        real, parameter   ::  ppdeltu = 1.e-4

        ! shallow convection
        real, parameter   ::  ppepsus = 3.e-4
        real, parameter   ::  ppdeltus = 3.e-4

        ! downdraft
        real, parameter   ::  ppepsd = 2.e-4
        real, parameter   ::  ppdeltd = 2.e-4

        ! * constants for precipitation parametrization
        real, parameter   ::  ppkmin = 1500.0
        real, parameter   ::  ppkval = 2.e-3

        ! * parameter beta determines the overshoot of the detrainment plumes
        ! above level of no buoyancy (beta=0 : no overshoot)

        ! penetrative and midlevel convection
        real, parameter   ::  ppbeta = 0.0
        ! shallow convection
        real, parameter   ::  ppbetas = 0.15

        ! * parameter gamma determines downward massflux at the level of free
        ! sinking (LFS)

        real, parameter   ::  ppgamma = 0.3

        ! parameter alpha determines specific humidity of air parcel at surface
        ! before starting the dry adiabatic ascent
        ! (alpha = 0 : air parcel has the specific humidity of the environment,
        !  alpha = 1 : air parcel has saturation specific humidity of the
        !              environment)
        real, parameter   ::  ppalpha = 0.2

        ! no of iterations for saturation calculation
        integer, parameter   ::   jpitermax = 5
 
        ! --- local --------------------------------------

        integer :: klc,klt,klfs

        real :: zfck,zamub,zamdlfs
        real :: zamu,zamd
        real :: zlc,zqc,ztc
        real :: zlcklc,zqcklc,ztcklc
        real :: ztd,zqd
 
        real :: zdq1,zdq2,zsv,zscv,zgp
        real :: zpgp(lme)

        real :: zepsu,zdeltu
        real :: zbeta
        integer :: jl,jjl,jiter
        integer :: imlc

        real  :: zdqcmin,zdqdmin
        real  :: ztenwb,zqenwb
        real  :: zttest,zqtest,zstv

        ! extra
        character(len=10) :: convection_type
 
        ! --- begin ----------------------------------------------

        convection_type = 'none'

        ! default cloud properties on each layer boundary

        do jl = 0, lme
          pamu(jl) = 0.0
          ptc(jl) = pt(jl)
          pqc(jl) = pq(jl)
          plc(jl) = 0.0
          pgp(jl) = 0.0
        end do

        ! air-masses (kg/m**2) in each layer, default entrainment/detrainment
        ! and precipitation rates

        do jl = 1, lme
          pam(jl)=(ppress(jl)-ppress(jl-1))/ppg
          peu(jl)=0.
          pdu(jl)=0.
          ped(jl)=0.
          pdd(jl)=0.
          zpgp(jl)=0.
        end do

        ! find condensation level by lifting an air parcel until supersaturation occurs

        ! we start the ascent with moisture and temperature of the upper boundary of
        ! the surface layer

        ztc = pt(lme-1)
        zqc = pq(lme-1) + ppalpha * ( qsat(pt(lme-1),ppress(lme-1)) - pq(lme-1) )

        do jl = lme-1-1, 1, -1

          ! preliminary value of parcel temperature (dry adiabatic ascent)
          ztc = ztc - ppg*(pz(jl)-pz(jl+1))/ppcpd

          ! check for supersaturation
          if ( zqc > qsat(ztc,ppress(jl)) ) then

            ! if supersaturation is detected we adjust moisture and
            ! temperature by condensation
            ! and set liquid water content equal to the condensate

            zdq2 = 0.0
            do jiter=1,jpitermax
              zdq1=(zqc-QSat(ztc,ppress(jl))) &
                   /(1.+ppzeta*dQSat_dt(ztc,ppress(jl)))
              zqc=zqc-zdq1
              ztc=ztc+zdq1*ppzeta
              zdq2=zdq2+zdq1
            end do
            zlc = zdq2

            ! check if parcel is buoyant:
            !   virtual temperature of parcel
            zscv = ztc*( 1.0 + ppvtcf*zqc - zlc )

            ! virtual temperature of environment
            zsv = pt(jl) * ( 1.0 + ppvtcf*pq(jl) )

            if ( zscv > zsv ) then
              ! if parcel is still buoyant then we have detected the cloud base
              klc=jl
              goto 20
            else
              !  if parcel is not buoyant then there is no penetrative or shallow
              !  convection
              if (llmilc) then
                goto 700
              else
                goto 3000
              endif
            endif
 
          endif

        end do

        ! no condensation level found
        goto 3000

20      continue

        ! if we arrive here a cloud base is detected:
        ! klc is cloud base level, ztc is cloud temperature, zqc moisture (at
        ! saturation)
        ! and zlc the liquid water content in the updraft at base level

        ztcklc = ztc
        zqcklc = zqc
        zlcklc = zlc

        ! calculate large scale moisture convergence below cloud base
        ! (use trapezoidal integration formula)

        !zfck=
        !    ((pq(klc)*pqam(klc)-pqac(klc))*pam(klc)+
        ! Correction Zoe Stockwell - Peter van Velthoven, 5 January 2000 &
        zfck = ((pq(klc)*pqam(klc) -pqac(klc) )*pam(klc+1)+ &
                (pq(klc)*pqam(lme)-pqac(lme))*pam(lme)      )

        do jl=klc+1,lme-1
          zfck=zfck+(pq(klc)*pqam(jl)-pqac(jl))*(pam(jl)+pam(jl+1))
        end do

        ! check for shallow or penetrative convection, set proper parameter
        ! values

        if (zfck.gt.0.) then
          ! penetrative and midlevel convection
          zepsu=ppepsu
          zdeltu=ppdeltu
          zbeta=ppbeta
          convection_type = 'deep'
        else
          ! shallow convection
          zepsu=ppepsus
          zdeltu=ppdeltus
          zbeta=ppbetas
          convection_type = 'shallow'
        endif

        zfck=pqtur+0.5*zfck


        ! if no moisture convergence then no penetrative or shallow
        ! convection is present

        if (zfck <= 0.0 ) then
          if (llmilc) then
             goto 700
          else
             goto 3000
          endif
        else
          goto 900
        endif

 700    continue

        ! check possibility for midlevel convection

        imlc = 0

        ! Check atmosphere from 2 layers above the surface to the middle of
        ! the atmosphere to see if midlevel convection might occur

        do jl=lme-1-1,lme-int(lme/2.),-1

          ! large scale ascent and an environmental relative humidity of more than
          ! 90% are needed for midlevel convection to occur

          if ( (pq(jl).gt.(0.9*qsat(pt(jl),ppress(jl))) ).and. &
              (pw(jl).lt.0.)) then

            if (imlc.eq.0) then
              ztc = pt(jl+1)
              zqc = pq(jl+1)
              zlc = 0.
              imlc = jl
            else if (imlc.gt.0) then
              if ((imlc-jl).eq.1) then
                imlc = jl
                goto 720
              else
                ztc = pt(jl+1)
                zqc = pq(jl+1)
                zlc = 0.
                imlc = jl
              end if
            end if

 720        continue

            ! do adiabatic ascent

            ztc = ztc - ppg*(pz(jl)-pz(jl+1))/ppcpd

            ! check for supersaturation

            if (zqc.gt.qsat(ztc,ppress(jl))) then
 
              ! if supersaturation is detected we adjust moisture and
              ! temperature by condensation
              ! and set liquid water content equal to the condensate

              zdq2=0.
              do jiter=1,jpitermax
                zdq1 = (zqc-qsat(ztc,ppress(jl)))/(1.+ppzeta*dqsat_dt(ztc,ppress(jl)))
                zqc=zqc-zdq1
                ztc=ztc+zdq1*ppzeta
                zdq2=zdq2+zdq1
              end do
              zlc=zdq2

            endif

            ! check if parcel is buoyant

            ! virtual temperature of parcel

            zscv = ztc*(1.+ppvtcf*zqc-zlc)

            ! virtual temperature of environment

            zsv = pt(jl)*(1.+ppvtcf*pq(jl))

            if (zscv.gt.zsv) then

              ! parcel is buoyant and we have detected the cloud base of midlevel
              ! convection

              klc = jl
              zamub = -pw(klc)/ppg
              zepsu = ppepsu
              zdeltu = ppdeltu
              zbeta = ppbeta
              llcudo = .false.
              convection_type = 'mid-level'
              goto 1000

            endif

          endif

        end do
 
        goto 3000

 900    continue

        ! massflux at base of cloud

        ! limit specific humidity difference between cloud and environment at
        ! cloud base

        zdqcmin = max(0.01*pq(klc),1.e-10)
        zdqcmin = max(zdqcmin,zqc+zlc-pq(klc))

        zamub=zfck/zdqcmin

 1000   continue

        ! limit mass flux at cloud base
        zamub=min(zamub,1.0)

        ! set updraft entrainment rates below cloud base proportional
        ! to layer air masses
        ! set updraft detrainment rates below cloud base to zero

        do jl=lme,klc+1,-1
          peu(jl) = zamub*pam(jl)*ppg/(ppress(lme)-ppress(klc))
          pdu(jl) = 0.0
        end do

        ! calculate now parcel ascent within cloud updraft
        !    cloud mass flux               zamu,
        !    cloud moisture                zqc,
        !    cloud temperature             ztc,
        !    cloud liquid water            zlc

        zamu = zamub

        do jl = klc,2,-1

          ! mass entrainment and detrainment
          peu(jl)=zepsu*zamu*(pz(jl-1)-pz(jl))
              pdu(jl)=zdeltu*zamu*(pz(jl-1)-pz(jl))
 
          ! preliminary values of cloud temperature, moisture and cloud liquid water
          ztc=ztc &
                -ppg*(pz(jl-1)-pz(jl))/ppcpd &
                +zepsu*(pz(jl-1)-pz(jl))*(pt(jl)-ztc)

              zqc=zqc &
                +zepsu*(pz(jl-1)-pz(jl))*(pq(jl)-zqc)

          zlc=zlc &
                +zepsu*(pz(jl-1)-pz(jl))*(-zlc)

          ! adjust moisture and temperature by condensation

          zdq2=0.
          do jiter=1,jpitermax
            zdq1=(zqc-qsat(ztc,ppress(jl))) &
                   /(1.+ppzeta*dqsat_dt(ztc,ppress(jl)))
            zqc=zqc-zdq1
            ztc=ztc+zdq1*ppzeta
            zdq2=zdq2+zdq1
          end do

          ! precipitation rate out of layer jl

          if ((pz(jl)+pz(jl-1))*0.5-pz(klc) .gt. ppkmin) then
            zgp=pam(jl)*ppkval/zamu
          else
            zgp=0.
          endif
 
          ! adjust liquid cloud water in updraft (use implicit scheme to prevent
          ! instability)

          zgp = zgp*zlc/(1+zgp)
 
          zpgp(jl) = zgp*zamu
 
          zlc = zlc-zgp+zdq2

          ! check for level of free sinking (LFS) where cumulus downdraft starts

          if (.not.llcudo) then
 
            ! downdraft calculation already done or turned off
 
            goto 800

          end if

          if ( zpgp(jl) == 0.0 ) then
 
            ! no downdraft exists since downdrafts are associated with convective
            ! precipitation from the updrafts
 
            goto 800

          end if

          if (jl.lt.3) then

            ! no downdraft since maximum possible cloud height is reached

            goto 800

          end if

          ! The LFS is assumed to be the highest model level where a mixture of equal
          ! parts of cloud air and environmental air (at wet bulb temperature) becomes
          ! negative buoyant with respect to the environmental air

          ! first step :
          ! calculate wet bulb temperature and moisture of the environmental air

          ztenwb = pt(jl-1)
          zqenwb = pq(jl-1)

          ! adjust temperature and moisture by evaporation
          ! zdq1 must be less equal 0 (zdq1=0 : environmental air is saturated,
          ! i.e. zqenwb = pq)

          do jiter = 1,jpitermax
            zdq1 = (zqenwb-qsat(ztenwb,ppress(jl-1)))/ &
                   (1.+ppzeta*dqsat_dt(ztenwb,ppress(jl-1)))
            zqenwb = zqenwb-zdq1
            ztenwb = ztenwb+zdq1*ppzeta
          end do

          ! second step :
          ! do mixing with cloud air

          zttest = 0.5*(ztc+ztenwb)
          zqtest = 0.5*(zqc+zqenwb)

          ! third step :
          ! check for negative buoyancy

          ! virtual temperature of the air mixture

          zstv = zttest*(1.+ppvtcf*zqtest)

          ! virtual temperature of the environmental air

          zsv = pt(jl-1)*(1.+ppvtcf*pq(jl-1))

          if (zstv.lt.zsv) then

            ! level of free sinking (LFS) is found, start downdraft calculation

            ! massflux at LFS is assumed to be directly proportional to the (preliminary)
            ! upward massflux at cloud base

            klfs = jl
            zamdlfs = -ppgamma*zamub
            zamd = zamdlfs
            ztd = zttest
            zqd = zqtest

            ped(klfs) = (-zamd)
            pdd(klfs) = 0.

            if (klfs.eq.klc) goto 45

            do jjl = klfs+1,klc,1

              ! mass entrainment and detrainment

              ped(jjl) = ppepsd*zamd*(pz(jjl)-pz(jjl-1))
              pdd(jjl) = ppdeltd*zamd*(pz(jjl)-pz(jjl-1))

              ! preliminary values of cloud temperature and moisture in downdraft

              ztd = ztd &
                      -ppg*(pz(jjl)-pz(jjl-1))/ppcpd &
                      +ppepsd*(pz(jjl)-pz(jjl-1))*(ztd-pt(jjl-1))

              zqd = zqd &
                      +ppepsd*(pz(jjl)-pz(jjl-1))*(zqd-pq(jjl-1))

              ! adjust moisture and temperature by evaporation

              do jiter=1,jpitermax
                zdq1 = (zqd-qsat(ztd,ppress(jjl-1)))/ &
                       (1.+ppzeta*dqsat_dt(ztd,ppress(jjl-1)))
                zqd = zqd-zdq1
                ztd = ztd+zdq1*ppzeta
              end do

              ! downdraft massflux at lower layer boundary

              zamd = zamd - ped(jjl) + pdd(jjl)

            end do
 45         continue

            zamd = min(0.,zamd)

            ! set downdraft detrainment rates below cloud base proportional to layer
            ! air masses
            ! set downdraft entrainment rates below cloud base to zero

            do jjl = lme,klc+1,-1
              ped(jjl) = 0.
              pdd(jjl) = zamd*pam(jjl)*ppg/(ppress(klc)-ppress(lme))
            end do

            ! recalculate updraft massflux at cloud base,
            ! now allowing for the downdraft massflux

            if (zamd.lt.0.) then
              zdqdmin = zqd-pq(klc)
              zamub = (zfck-zamd*zdqdmin)/zdqcmin
              if (zamub.le.0.) then
                do jjl=1,lme
                  peu(jjl)=0.
                  pdu(jjl)=0.
                  ped(jjl)=0.
                  pdd(jjl)=0.
                end do
                goto 3000
              endif

              ! go back to cloud base and start updraft calculation again

              ztc = ztcklc
              zqc = zqcklc
              zlc = zlcklc
              llcudo = .false.
              goto 1000

            else
              goto 800
            endif

          else

            goto 800

          endif

 800      continue

          ! check for buoyancy (in updraft)
          !   virtual temperature in updraft at upper boundary of layer jl:
          zscv=ztc*(1.+ppvtcf*zqc-zlc)

          !   virtual temperature outside of cloud
          zsv=pt(jl-1)*(1.+ppvtcf*pq(jl-1))

          if ( zscv <= zsv ) then
            klt=jl
            goto 400
          endif

          ! updraft massflux at upper layer boundary
          zamu=zamu+peu(jl)-pdu(jl)

          ! store cloud properties on upper layer boundary
          ptc(jl-1)=ztc
          pqc(jl-1)=zqc
          plc(jl-1)=zlc
 
        end do

        klt=2

400     continue

        ! set detrainment in two layers above cloud

        pdu(klt-1)=zbeta*zamu
        peu(klt-1)=0.
        pdu(klt)=(1-zbeta)*zamu
        peu(klt)=0.

        ! add up rainrate on each of the layer boundaries

        do jl=klt+1,lme
          pgp(jl)=pgp(jl-1)+zpgp(jl)
        end do

        ! determine net mass flux on each of the layer boundaries

        do jl=lme,1,-1
          pamu(jl-1)=pamu(jl)+(peu(jl)-pdu(jl))-(ped(jl)-pdd(jl))
        end do

        llcudo = .true.
        llmilc = .true.

        return

        ! no cloud present, set cloud base and top to 0 and return

3000    continue
        klc=0
        klt=0
        llcudo = .true.
        llmilc = .true.
        convection_type = 'none'  

      end subroutine cloud


  ! ==========================================================
  ! ===
  ! ===  TMPP/SOURCE/tmpp_sub_work.f90
  ! ===
  ! ==========================================================
  

  !---------------------------------------------------------------------
  ! calculate en/detrainment rates and diffusion coefficient on TM grid
  !---------------------------------------------------------------------
  ! History:
  ! Increased vertical dimension of z,t,q,u,v,w from lme to lme + 1
  ! in order to be able to use the same memory location in worksub
  ! for u and wu, for t and wt, etc.
  ! Added subroutine cen2bound for the same reason
  ! Removed dummy fields for geopotential height and zonal means
  ! Program just fits into memory of SGI machines (max stacksize = 524288) if 
  !  TM and ECMWF both have 1x1 degree resolution and 60 levels
  !                                                         pvv, 5-2-2000
  !---------------------------------------------------------------------

  ! =========================================================
  
  ! interpolate variables from the center of parent-model layers to the
  ! boundaries of parent-model layers and save result in same memory location

  
  subroutine mid2bound_uv( lme, npe, u, v, ps, mask, a, b )
  
    ! --- in/out ----------------------------------
    
    integer, intent(in)             ::  lme, npe
    real, intent(inout)             ::  u(npe,0:lme)
    real, intent(inout)             ::  v(npe,0:lme)
    real, intent(in)                ::  ps(npe)
    logical, intent(in)             ::  mask(npe)
    real, intent(in)                ::  a(0:lme)
    real, intent(in)                ::  b(0:lme)
    
    ! --- local ---------------------------------
    
    integer              ::  i
    real                 ::  wcol(0:lme)
    
    ! --- begin -------------------------
    
    do i = 1, npe
      if ( mask(i) ) then

        call cen2bound_col( lme, u(i,1:lme), ps(i), 1, wcol, a, b )
        u(i,:) = wcol

        call cen2bound_col( lme, v(i,1:lme), ps(i), 1, wcol, a, b )
        v(i,:) = wcol

      end if
    end do
    
  end subroutine mid2bound_uv
    
    
  ! ===
  
  
  subroutine mid2bound_w( lme, npe, w, ps, mask, a, b )
  
    ! --- in/out ----------------------------------
    
    integer, intent(in)             ::  lme, npe
    real, intent(inout)             ::  w(npe,0:lme)
    real, intent(in)                ::  ps(npe)
    logical, intent(in)             ::  mask(npe)
    real, intent(in)                ::  a(0:lme)
    real, intent(in)                ::  b(0:lme)
    
    ! --- local ---------------------------------
    
    integer              ::  i
    real                 ::  wcol(0:lme)
    
    ! --- begin -------------------------
    
    do i = 1, npe
      if ( mask(i) ) then

        call cen2bound_col( lme, w(i,1:lme), ps(i), 1, wcol, a, b )
        w(i,:) = wcol

        ! set to zero at top
        w(i,0) = 0.0

      end if
    end do
    
  end subroutine mid2bound_w
    
    
  ! ===
  
  
  subroutine mid2bound_t( lme, npe, t, ps, mask, a, b )
  
    ! --- in/out ----------------------------------
    
    integer, intent(in)             ::  lme, npe
    real, intent(inout)             ::  t(npe,0:lme)
    real, intent(in)                ::  ps(npe)
    logical, intent(in)             ::  mask(npe)
    real, intent(in)                ::  a(0:lme)
    real, intent(in)                ::  b(0:lme)
    
    ! --- local ---------------------------------
    
    integer              ::  i
    real                 ::  wcol(0:lme)
    
    ! --- begin -------------------------
    
    do i = 1, npe
      if ( mask(i) ) then

        call cen2bound_col( lme, t(i,1:lme), ps(i), 1, wcol, a, b )
        t(i,:) = wcol

      end if
    end do
    
  end subroutine mid2bound_t
    
    
  ! ===
  
  
  subroutine mid2bound_q( lme, npe, q, ps, mask, a, b, t )
  
    use phys_humidity, only : qsat
  
    ! --- in/out ----------------------------------
    
    integer, intent(in)             ::  lme, npe
    real, intent(inout)             ::  q(npe,0:lme)
    real, intent(in)                ::  ps(npe)
    logical, intent(in)             ::  mask(npe)
    real, intent(in)                ::  a(0:lme)
    real, intent(in)                ::  b(0:lme)
    real, intent(in)                ::  t(npe,0:lme)
    
    ! --- local ---------------------------------
    
    integer              ::  i, l
    real                 ::  wcol(0:lme)
    real            :: tmpress(0:lme)
    
    ! --- begin -------------------------
    
    do i = 1, npe
      if ( mask(i) ) then

        call cen2bound_col( lme, q(i,1:lme), ps(i), 1, wcol, a, b )
        q(i,:) = wcol

        ! limit specific humidity at 0 and qsat(t,p)
        ! first establish hybrid vertical coordinate at i,j ;  
        ! note that ps is expressed in Pa
        do l = 0, lme
          tmpress(l) = a(l) + ps(i)*b(l)
        end do
        do l = 0, lme
          q(i,l) = min( qsat(t(i,l),tmpress(l)), max(0.0,q(i,l)) )
        end do

      end if
    end do
    
  end subroutine mid2bound_q
    
    
  ! ===
  
  
  subroutine mid2bound_z( lme, npe, z, ps, mask, a, b, zsurf )
  
    use Binas, only : g => grav

    ! --- in/out ----------------------------------
    
    integer, intent(in)             ::  lme, npe
    real, intent(inout)             ::  z(npe,0:lme)
    real, intent(in)                ::  ps(npe)
    logical, intent(in)             ::  mask(npe)
    real, intent(in)                ::  a(0:lme)
    real, intent(in)                ::  b(0:lme)
    real, intent(in)                ::  zsurf(npe)
    
    ! --- local ---------------------------------
    
    integer              ::  i
    real                 ::  wcol(0:lme)

    ! --- begin -------------------------
    
    do i = 1, npe
      if ( mask(i) ) then

        call cen2bound_col( lme, z(i,1:lme), ps(i), 1, wcol, a, b )
        z(i,:) = wcol

        ! set to known value at surface:
        z(i,lme) = zsurf(i)/g

      end if
    end do
    
  end subroutine mid2bound_z
    
    
  ! ===
  
  
  subroutine mid2bound_p( lme, npe, p, ps, mask, a, b )
  
    ! --- in/out ----------------------------------
    
    integer, intent(in)             ::  lme, npe
    real, intent(inout)             ::  p(npe,0:lme)
    real, intent(in)                ::  ps(npe)
    logical, intent(in)             ::  mask(npe)
    real, intent(in)                ::  a(0:lme)
    real, intent(in)                ::  b(0:lme)
    
    ! --- local ---------------------------------
    
    integer              ::  i, j
    real                 ::  wcol(0:lme)
    
    ! --- begin -------------------------
    
    do i = 1, npe
      if ( mask(i) ) then

        call cen2bound_col( lme, p(i,:), ps(i), 0, wcol, a, b )
        p(i,:) = wcol

      end if
    end do
    
  end subroutine mid2bound_p
  
  

  ! interpolate from mid-levels of parent model to
  ! level boundaries
  !
  ! Peter van Velthoven  - 4 January 2000
  ! This subroutine is included to save memory
  ! field and field2 use the same space in memory
  !
  ! iopt = 0    : no field as input : fill wfield with pressure
  ! iopt = other: use field as input
  !

  subroutine cen2bound_col( lme, field, ps, iopt, wfield, a, b )

    use Num, only : interp_muherm

    ! --- in/out -------------------------------

    integer, intent(in)  :: lme
    real, intent(in)     :: field(lme)    ! input on (mid-)levels
    real, intent(in)     :: ps            ! surface pressure
    integer, intent(in)  :: iopt
    real, intent(out)    :: wfield(0:lme) ! output on vertical level boundaries
    real, intent(in)     ::  a(0:lme)
    real, intent(in)     ::  b(0:lme)
    
    ! --- begin -------------------------------

    integer       ::  status
    real          ::  ztemp(lme)
    real          ::  tmpress(0:lme) ! pressure at ECMWF vertical level boundaries
    real          ::  tcmpress(lme)  ! pressure at ECMWF (mid-)levels
    integer       ::  l
    
    ! --- begin --------------------------------

    ! establish hybrid vertical coordinate at i,j
    ! note that ps is expressed in Pa
    tmpress = a + ps * b

    ! calculate pressure at model layer center
    do l=1,lme
      tcmpress(l) = (tmpress(l-1)+tmpress(l))/2.
    end do

    if ( iopt == 0 ) then
      wfield = tmpress
    else
      call interp_muherm( tcmpress, field, tmpress, wfield, status )
      if (status/=0) stop 'ERROR in cen2bound_col'
    end if


  end subroutine cen2bound_col
  
  
  
  ! ==========================================================
  ! ===
  ! ===  TMPP/SOURCE/tmpp_geometry.f90
  ! ===
  ! ==========================================================
  

  ! pressure at half leves from bottom to top
  
  subroutine phlev_ec_rev( lme, a_ec, b_ec, ps, pb )
  
    ! --- in/out --------------------------
    
    integer, intent(in)        ::  lme
    real, intent(in)           ::  a_ec(0:lme)
    real, intent(in)           ::  b_ec(0:lme)
    real, intent(in)           ::  ps
    real, intent(out)          ::  pb(0:lme)
    
    ! --- local --------------------------
    
    integer          ::  l
    
    ! --- in/out -------------------------
    
    do l = 0, lme
      pb(lme-l) = a_ec(l) + b_ec(l) * ps
    end do
    
  end subroutine phlev_ec_rev
  
  
  ! ==========================================================
  ! ===
  ! ===  TMPP/SOURCE/phys_geopot.f90
  ! ===
  ! ==========================================================
  

  !
  ! NAME
  !   GeoPot_col  -  calculate geopotential height
  !
  ! DESCRIPTION
  !   Calculate geopotential height from halflevel pressures
  !   and full level virtual temperature.
  !
  ! USAGE
  !  
  !   call GeoPot( z, zsurf, pt, pq, pb, lm )
  !
  !     integer, intent(in)  ::  lm        ! number of levels 
  !
  !       (levels numbered downwards (top -> down) )
  !
  !     real, intent(out)    ::  z(lm)     ! geopotential height (m ?).
  !     real, intent(in)     ::  zsurf     ! orography (m ?)
  !     real, intent(in)     ::  pt(lm)    ! temperature (K ?)
  !     real, intent(in)     ::  pq(lm)    ! specific humidity (??)
  !
  !       (levels numbered upwards (bottom -> up) )
  !
  !     real, intent(in)     ::  pb(0:lm)  ! pressure at half levels
  !     
  ! HISTORY
  !
  !   06-11-2001, Arjo Segers
  !     Extracted from original routines 'geopot' and 'auxhyb'
  !     by Ad Jeuken
  !       

  subroutine GeoPot( lm, zsurf, pt, pq, pb, z )

    ! --- in/out -------------------------

    integer, intent(in)  ::  lm
    real, intent(out)    ::  z(lm)
    real, intent(in)     ::  zsurf
    real, intent(in)     ::  pt(lm)
    real, intent(in)     ::  pq(lm)
    real, intent(in)     ::  pb(0:lm)

    ! --- const ------------------------

    real, parameter      ::  rd = 287.05
    real, parameter      ::  g0 = 9.801

    ! --- local ------------------------------

    integer  ::  linv
    real     ::  pdelp, prdelp

    real     ::  palfa(lm)
    real     ::  plnr(lm)

    real     ::  tv(lm)

    integer  ::  l, lp1

    ! --- begin ---------------------------------
    
    ! >>> former routine 'auxhyb' >>>>>>>>>>>>>>>>>>>>>>

    ! loop from top to bottom: 
    do l = 1, lm
      linv = lm - l   ! lm-1, 0
      pdelp = pb(linv) - pb(linv+1)
      prdelp = 1.0 / pdelp
      if ( l == 1 ) then
        plnr(l) = rd * 1.3862944
      else
        plnr(l) = rd * log( pb(linv)/pb(linv+1) )
      end if
      palfa(l) = rd - pb(linv+1) * plnr(l) * prdelp
    end do

    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

    ! loop from bottom to top:
    do l = lm, 1, -1
      tv(l) = pt(l) * ( 1.0 + 0.608*pq(l) )
      if ( l == lm ) then
        z(l) = palfa(l)*tv(l)/g0 + zsurf/g0
      else
        lp1=l+1
        z(l) = z(lp1) + ( palfa(l)*tv(l) + (plnr(lp1)-palfa(lp1))*tv(lp1) )/g0
      end if
    end do

  end subroutine GeoPot



end module phys_convec_tmpp