source: LMDZ5/trunk/libf/phylmd/convect1.F90 @ 2250

! $Header$

SUBROUTINE convect1(len, nd, ndp1, noff, minorig, t, q, qs, u, v, p, ph, &
    iflag, ft, fq, fu, fv, precip, cbmf, delt, ma)
  ! .............................START PROLOGUE............................

  ! SCCS IDENTIFICATION:  @(#)convect1.f        1.1 04/21/00
  ! 19:40:52 /h/cm/library/nogaps4/src/sub/fcst/convect1.f_v

  ! CONFIGURATION IDENTIFICATION:  None

  ! MODULE NAME:  convect1

  ! DESCRIPTION:

  ! convect1     The Emanuel Cumulus Convection Scheme

  ! CONTRACT NUMBER AND TITLE:  None

  ! REFERENCES: Programmers  K. Emanuel (MIT), Timothy F. Hogan, M. Peng
  ! (NRL)

  ! CLASSIFICATION:  Unclassified

  ! RESTRICTIONS: None

  ! COMPILER DEPENDENCIES: FORTRAN 77, FORTRAN 90

  ! COMPILE OPTIONS: Fortran 77: -Zu -Wf"-ei -o aggress"
  ! Fortran 90: -O vector3,scalar3,task1,aggress,overindex  -ei -r 2

  ! LIBRARIES OF RESIDENCE: /a/ops/lib/libfcst159.a

  ! USAGE: call convect1(len,nd,noff,minorig,
  ! &                   t,q,qs,u,v,
  ! &                   p,ph,iflag,ft,
  ! &                   fq,fu,fv,precip,cbmf,delt)

  ! PARAMETERS:
  ! Name            Type         Usage            Description
  ! ----------      ----------     -------  ----------------------------

  ! len           Integer        Input        first (i) dimension
  ! nd            Integer        Input        vertical (k) dimension
  ! ndp1          Integer        Input        nd + 1
  ! noff          Integer        Input        integer limit for convection
  ! (nd-noff)
  ! minorig       Integer        Input        First level of convection
  ! t             Real           Input        temperature
  ! q             Real           Input        specific hum
  ! qs            Real           Input        sat specific hum
  ! u             Real           Input        u-wind
  ! v             Real           Input        v-wind
  ! p             Real           Input        full level pressure
  ! ph            Real           Input        half level pressure
  ! iflag         Integer        Output       iflag on latitude strip
  ! ft            Real           Output       temp tend
  ! fq            Real           Output       spec hum tend
  ! fu            Real           Output       u-wind tend
  ! fv            Real           Output       v-wind tend
  ! cbmf          Real           In/Out       cumulus mass flux
  ! delt          Real           Input        time step
  ! iflag         Integer        Output       integer flag for Emanuel
  ! conditions

  ! COMMON BLOCKS:
  ! Block      Name     Type    Usage              Notes
  ! --------  --------   ----    ------   ------------------------

  ! FILES: None

  ! DATA BASES: None

  ! NON-FILE INPUT/OUTPUT: None

  ! ERROR CONDITIONS: None

  ! ADDITIONAL COMMENTS: None

  ! .................MAINTENANCE SECTION................................

  ! MODULES CALLED:
  ! Name           Description
  ! convect2        Emanuel cumulus convection tendency calculations
  ! -------     ----------------------
  ! LOCAL VARIABLES AND
  ! STRUCTURES:
  ! Name     Type    Description
  ! -------  ------  -----------
  ! See Comments Below

  ! i        Integer loop index
  ! k        Integer loop index

  ! METHOD:

  ! See Emanuel, K. and M. Zivkovic-Rothman, 2000: Development and evaluation
  ! of a
  ! convective scheme for use in climate models.

  ! FILES: None

  ! INCLUDE FILES: None

  ! MAKEFILE: /a/ops/met/nogaps/src/sub/fcst/fcst159lib.mak

  ! ..............................END PROLOGUE.............................


  USE dimphy
  IMPLICIT NONE

  ! ym#include "dimensions.h"
  ! ym#include "dimphy.h"

  INTEGER len
  INTEGER nd
  INTEGER ndp1
  INTEGER noff
  REAL t(len, nd)
  REAL q(len, nd)
  REAL qs(len, nd)
  REAL u(len, nd)
  REAL v(len, nd)
  REAL p(len, nd)
  REAL ph(len, ndp1)
  INTEGER iflag(len)
  REAL ft(len, nd)
  REAL fq(len, nd)
  REAL fu(len, nd)
  REAL fv(len, nd)
  REAL precip(len)
  REAL cbmf(len)
  REAL ma(len, nd)
  INTEGER minorig
  REAL delt, cpd, cpv, cl, rv, rd, lv0, g
  REAL sigs, sigd, elcrit, tlcrit, omtsnow, dtmax, damp
  REAL alpha, entp, coeffs, coeffr, omtrain, cu

  ! -------------------------------------------------------------------
  ! --- ARGUMENTS
  ! -------------------------------------------------------------------
  ! --- 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.

  ! nl:  The maximum number of levels to which convection can penetrate, plus
  ! 1.
  ! NL MUST be less than or equal to ND-1.

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

  ! ----------------------------------------------------------------------------
  ! ---   On Output:

  ! iflag: An output integer whose value denotes the following:
  ! VALUE   INTERPRETATION
  ! -----   --------------
  ! 0     Moist convection occurs.
  ! 1     Moist convection occurs, but a CFL condition
  ! on the subsidence warming is violated. This
  ! does not cause the scheme to terminate.
  ! 2     Moist convection, but no precip because ep(inb) lt 0.0001
  ! 3     No moist convection because new cbmf is 0 and old cbmf is 0.
  ! 4     No moist convection; atmosphere is not
  ! unstable
  ! 6     No moist convection because ihmin le minorig.
  ! 7     No moist convection because unreasonable
  ! parcel level temperature or specific humidity.
  ! 8     No moist convection: lifted condensation
  ! level is above the 200 mb level.
  ! 9     No moist convection: cloud base is higher
  ! then the level NL-1.

  ! 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.

  ! det:   Array of detrainment mass flux of dimension ND.

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

  ! Local arrays

  INTEGER nl
  INTEGER nlp
  INTEGER nlm
  INTEGER i, k, n
  REAL delti
  REAL rowl
  REAL clmcpv
  REAL clmcpd
  REAL cpdmcp
  REAL cpvmcpd
  REAL eps
  REAL epsi
  REAL epsim1
  REAL ginv
  REAL hrd
  REAL prccon1
  INTEGER icbmax
  REAL lv(klon, klev)
  REAL cpn(klon, klev)
  REAL cpx(klon, klev)
  REAL tv(klon, klev)
  REAL gz(klon, klev)
  REAL hm(klon, klev)
  REAL h(klon, klev)
  REAL work(klon)
  INTEGER ihmin(klon)
  INTEGER nk(klon)
  REAL rh(klon)
  REAL chi(klon)
  REAL plcl(klon)
  INTEGER icb(klon)
  REAL tnk(klon)
  REAL qnk(klon)
  REAL gznk(klon)
  REAL pnk(klon)
  REAL qsnk(klon)
  REAL ticb(klon)
  REAL gzicb(klon)
  REAL tp(klon, klev)
  REAL tvp(klon, klev)
  REAL clw(klon, klev)

  REAL ah0(klon), cpp(klon)
  REAL tg, qg, s, alv, tc, ahg, denom, es, rg

  INTEGER ncum
  INTEGER idcum(klon)

  cpd = 1005.7
  cpv = 1870.0
  cl = 4190.0
  rv = 461.5
  rd = 287.04
  lv0 = 2.501E6
  g = 9.8

  ! *** 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)                 ***

  sigs = 0.12
  sigd = 0.05
  elcrit = 0.0011
  tlcrit = -55.0
  omtsnow = 5.5
  dtmax = 0.9
  damp = 0.1
  alpha = 0.2
  entp = 1.5
  coeffs = 0.8
  coeffr = 1.0
  omtrain = 50.0

  cu = 0.70
  damp = 0.1


  ! Define nl, nlp, nlm, and delti

  nl = nd - noff
  nlp = nl + 1
  nlm = nl - 1
  delti = 1.0/delt

  ! -------------------------------------------------------------------
  ! --- SET CONSTANTS
  ! -------------------------------------------------------------------

  rowl = 1000.0
  clmcpv = cl - cpv
  clmcpd = cl - cpd
  cpdmcp = cpd - cpv
  cpvmcpd = cpv - cpd
  eps = rd/rv
  epsi = 1.0/eps
  epsim1 = epsi - 1.0
  ginv = 1.0/g
  hrd = 0.5*rd
  prccon1 = 86400.0*1000.0/(rowl*g)

  ! dtmax is the maximum negative temperature perturbation.

  ! =====================================================================
  ! --- INITIALIZE OUTPUT ARRAYS AND PARAMETERS
  ! =====================================================================

  DO k = 1, nd
    DO i = 1, len
      ft(i, k) = 0.0
      fq(i, k) = 0.0
      fu(i, k) = 0.0
      fv(i, k) = 0.0
      tvp(i, k) = 0.0
      tp(i, k) = 0.0
      clw(i, k) = 0.0
      gz(i, k) = 0.
    END DO
  END DO
  DO i = 1, len
    precip(i) = 0.0
    iflag(i) = 0
  END DO

  ! =====================================================================
  ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
  ! =====================================================================
  DO k = 1, nl + 1
    DO i = 1, len
      lv(i, k) = lv0 - clmcpv*(t(i,k)-273.15)
      cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k)
      cpx(i, k) = cpd*(1.0-q(i,k)) + cl*q(i, k)
      tv(i, k) = t(i, k)*(1.0+q(i,k)*epsim1)
    END DO
  END DO

  ! gz = phi at the full levels (same as p).

  DO i = 1, len
    gz(i, 1) = 0.0
  END DO
  DO k = 2, nlp
    DO i = 1, len
      gz(i, k) = gz(i, k-1) + hrd*(tv(i,k-1)+tv(i,k))*(p(i,k-1)-p(i,k))/ph(i, &
        k)
    END DO
  END DO

  ! h  = phi + cpT (dry static energy).
  ! hm = phi + cp(T-Tbase)+Lq

  DO k = 1, nlp
    DO i = 1, len
      h(i, k) = gz(i, k) + cpn(i, k)*t(i, k)
      hm(i, k) = gz(i, k) + cpx(i, k)*(t(i,k)-t(i,1)) + lv(i, k)*q(i, k)
    END DO
  END DO

  ! -------------------------------------------------------------------
  ! --- Find level of minimum moist static energy
  ! --- If level of minimum moist static energy coincides with
  ! --- or is lower than minimum allowable parcel origin level,
  ! --- set iflag to 6.
  ! -------------------------------------------------------------------
  DO i = 1, len
    work(i) = 1.0E12
    ihmin(i) = nl
  END DO
  DO k = 2, nlp
    DO i = 1, len
      IF ((hm(i,k)<work(i)) .AND. (hm(i,k)<hm(i,k-1))) THEN
        work(i) = hm(i, k)
        ihmin(i) = k
      END IF
    END DO
  END DO
  DO i = 1, len
    ihmin(i) = min(ihmin(i), nlm)
    IF (ihmin(i)<=minorig) THEN
      iflag(i) = 6
    END IF
  END DO

  ! -------------------------------------------------------------------
  ! --- Find that model level below the level of minimum moist static
  ! --- energy that has the maximum value of moist static energy
  ! -------------------------------------------------------------------

  DO i = 1, len
    work(i) = hm(i, minorig)
    nk(i) = minorig
  END DO
  DO k = minorig + 1, nl
    DO i = 1, len
      IF ((hm(i,k)>work(i)) .AND. (k<=ihmin(i))) THEN
        work(i) = hm(i, k)
        nk(i) = k
      END IF
    END DO
  END DO
  ! -------------------------------------------------------------------
  ! --- Check whether parcel level temperature and specific humidity
  ! --- are reasonable
  ! -------------------------------------------------------------------
  DO i = 1, len
    IF (((t(i,nk(i))<250.0) .OR. (q(i,nk(i))<=0.0) .OR. (p(i,ihmin(i))< &
      400.0)) .AND. (iflag(i)==0)) iflag(i) = 7
  END DO
  ! -------------------------------------------------------------------
  ! --- Calculate lifted condensation level of air at parcel origin level
  ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
  ! -------------------------------------------------------------------
  DO i = 1, len
    tnk(i) = t(i, nk(i))
    qnk(i) = q(i, nk(i))
    gznk(i) = gz(i, nk(i))
    pnk(i) = p(i, nk(i))
    qsnk(i) = qs(i, nk(i))

    rh(i) = qnk(i)/qsnk(i)
    rh(i) = min(1.0, rh(i))
    chi(i) = tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
    plcl(i) = pnk(i)*(rh(i)**chi(i))
    IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag(i &
      ) = 8
  END DO
  ! -------------------------------------------------------------------
  ! --- Calculate first level above lcl (=icb)
  ! -------------------------------------------------------------------
  DO i = 1, len
    icb(i) = nlm
  END DO

  DO k = minorig, nl
    DO i = 1, len
      IF ((k>=(nk(i)+1)) .AND. (p(i,k)<plcl(i))) icb(i) = min(icb(i), k)
    END DO
  END DO

  DO i = 1, len
    IF ((icb(i)>=nlm) .AND. (iflag(i)==0)) iflag(i) = 9
  END DO

  ! Compute icbmax.

  icbmax = 2
  DO i = 1, len
    icbmax = max(icbmax, icb(i))
  END DO

  ! -------------------------------------------------------------------
  ! --- Calculates the lifted parcel virtual temperature at nk,
  ! --- the actual temperature, and the adiabatic
  ! --- liquid water content. The procedure is to solve the equation.
  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
  ! -------------------------------------------------------------------

  DO i = 1, len
    tnk(i) = t(i, nk(i))
    qnk(i) = q(i, nk(i))
    gznk(i) = gz(i, nk(i))
    ticb(i) = t(i, icb(i))
    gzicb(i) = gz(i, icb(i))
  END DO

  ! ***  Calculate certain parcel quantities, including static energy   ***

  DO i = 1, len
    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
      273.15)) + gznk(i)
    cpp(i) = cpd*(1.-qnk(i)) + qnk(i)*cpv
  END DO

  ! ***   Calculate lifted parcel quantities below cloud base   ***

  DO k = minorig, icbmax - 1
    DO i = 1, len
      tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))/cpp(i)
      tvp(i, k) = tp(i, k)*(1.+qnk(i)*epsi)
    END DO
  END DO

  ! ***  Find lifted parcel quantities above cloud base    ***

  DO i = 1, len
    tg = ticb(i)
    qg = qs(i, icb(i))
    alv = lv0 - clmcpv*(ticb(i)-273.15)

    ! First iteration.

    s = cpd + alv*alv*qg/(rv*ticb(i)*ticb(i))
    s = 1./s
    ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
    tg = tg + s*(ah0(i)-ahg)
    tg = max(tg, 35.0)
    tc = tg - 273.15
    denom = 243.5 + tc
    IF (tc>=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 = eps*es/(p(i,icb(i))-es*(1.-eps))

    ! Second iteration.

    s = cpd + alv*alv*qg/(rv*ticb(i)*ticb(i))
    s = 1./s
    ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
    tg = tg + s*(ah0(i)-ahg)
    tg = max(tg, 35.0)
    tc = tg - 273.15
    denom = 243.5 + tc
    IF (tc>=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 = eps*es/(p(i,icb(i))-es*(1.-eps))

    alv = lv0 - clmcpv*(ticb(i)-273.15)
    tp(i, icb(i)) = (ah0(i)-(cl-cpd)*qnk(i)*ticb(i)-gz(i,icb(i))-alv*qg)/cpd
    clw(i, icb(i)) = qnk(i) - qg
    clw(i, icb(i)) = max(0.0, clw(i,icb(i)))
    rg = qg/(1.-qnk(i))
    tvp(i, icb(i)) = tp(i, icb(i))*(1.+rg*epsi)
  END DO

  DO k = minorig, icbmax
    DO i = 1, len
      tvp(i, k) = tvp(i, k) - tp(i, k)*qnk(i)
    END DO
  END DO

  ! -------------------------------------------------------------------
  ! --- Test for instability.
  ! --- If there was no convection at last time step and parcel
  ! --- is stable at icb, then set iflag to 4.
  ! -------------------------------------------------------------------

  DO i = 1, len
    IF ((cbmf(i)==0.0) .AND. (iflag(i)==0) .AND. (tvp(i, &
      icb(i))<=(tv(i,icb(i))-dtmax))) iflag(i) = 4
  END DO

  ! =====================================================================
  ! --- IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY
  ! =====================================================================

  ncum = 0
  DO i = 1, len
    IF (iflag(i)==0) THEN
      ncum = ncum + 1
      idcum(ncum) = i
    END IF
  END DO

  ! Call convect2, which compresses the points and computes the heating,
  ! moistening, velocity mixing, and precipiation.

  ! print*,'cpd avant convect2 ',cpd
  IF (ncum>0) THEN
    CALL convect2(ncum, idcum, len, nd, ndp1, nl, minorig, nk, icb, t, q, qs, &
      u, v, gz, tv, tp, tvp, clw, h, lv, cpn, p, ph, ft, fq, fu, fv, tnk, &
      qnk, gznk, plcl, precip, cbmf, iflag, delt, cpd, cpv, cl, rv, rd, lv0, &
      g, sigs, sigd, elcrit, tlcrit, omtsnow, dtmax, damp, alpha, entp, &
      coeffs, coeffr, omtrain, cu, ma)
  END IF

  RETURN
END SUBROUTINE convect1