Ignore:
Timestamp:
Mar 5, 2014, 2:19:12 PM (11 years ago)
Author:
lguez
Message:

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/convect1.F90

    r1988 r1992  
    1 !
     1
    22! $Header$
    3 !
    4       subroutine convect1(len,nd,ndp1,noff,minorig,
    5      &                   t,q,qs,u,v,
    6      &                   p,ph,iflag,ft,
    7      &                   fq,fu,fv,precip,cbmf,delt,Ma)
    8 C.............................START PROLOGUE............................
    9 C
    10 C SCCS IDENTIFICATION:  @(#)convect1.f  1.1 04/21/00
    11 C                       19:40:52 /h/cm/library/nogaps4/src/sub/fcst/convect1.f_v
    12 C
    13 C CONFIGURATION IDENTIFICATION:  None
    14 C
    15 C MODULE NAME:  convect1
    16 C
    17 C DESCRIPTION:
    18 C
    19 C convect1     The Emanuel Cumulus Convection Scheme
    20 C
    21 C CONTRACT NUMBER AND TITLE:  None
    22 C
    23 C REFERENCES: Programmers  K. Emanuel (MIT), Timothy F. Hogan, M. Peng (NRL)
    24 C
    25 C CLASSIFICATION:  Unclassified
    26 C
    27 C RESTRICTIONS: None
    28 C
    29 C COMPILER DEPENDENCIES: FORTRAN 77, FORTRAN 90
    30 C
    31 C COMPILE OPTIONS: Fortran 77: -Zu -Wf"-ei -o aggress"
    32 C                  Fortran 90: -O vector3,scalar3,task1,aggress,overindex  -ei -r 2
    33 C
    34 C LIBRARIES OF RESIDENCE: /a/ops/lib/libfcst159.a
    35 C
    36 C USAGE: call convect1(len,nd,noff,minorig,
    37 C    &                   t,q,qs,u,v,
    38 C    &                   p,ph,iflag,ft,
    39 C    &                   fq,fu,fv,precip,cbmf,delt)
    40 C
    41 C PARAMETERS:
    42 C      Name            Type         Usage            Description
    43 C   ----------      ----------     -------  ----------------------------
    44 C
    45 C      len           Integer        Input        first (i) dimension
    46 C      nd            Integer        Input        vertical (k) dimension
    47 C      ndp1          Integer        Input        nd + 1
    48 C      noff          Integer        Input        integer limit for convection (nd-noff)
    49 C      minorig       Integer        Input        First level of convection
    50 C      t             Real           Input        temperature
    51 C      q             Real           Input        specific hum
    52 C      qs            Real           Input        sat specific hum
    53 C      u             Real           Input        u-wind
    54 C      v             Real           Input        v-wind
    55 C      p             Real           Input        full level pressure
    56 C      ph            Real           Input        half level pressure
    57 C      iflag         Integer        Output       iflag on latitude strip
    58 C      ft            Real           Output       temp tend
    59 C      fq            Real           Output       spec hum tend
    60 C      fu            Real           Output       u-wind tend
    61 C      fv            Real           Output       v-wind tend
    62 C      cbmf          Real           In/Out       cumulus mass flux
    63 C      delt          Real           Input        time step
    64 C      iflag         Integer        Output       integer flag for Emanuel conditions
    65 C
    66 C COMMON BLOCKS:
    67 C      Block      Name     Type    Usage              Notes
    68 C     --------  --------   ----    ------   ------------------------
    69 C
    70 C FILES: None
    71 C
    72 C DATA BASES: None
    73 C
    74 C NON-FILE INPUT/OUTPUT: None
    75 C
    76 C ERROR CONDITIONS: None
    77 C
    78 C ADDITIONAL COMMENTS: None
    79 C
    80 C.................MAINTENANCE SECTION................................
    81 C
    82 C MODULES CALLED:
    83 C         Name           Description
    84 C         convect2        Emanuel cumulus convection tendency calculations
    85 C        -------     ----------------------
    86 C LOCAL VARIABLES AND
    87 C          STRUCTURES:
    88 C Name     Type    Description
    89 C -------  ------  -----------
    90 C See Comments Below
    91 C
    92 C i        Integer loop index
    93 C k        Integer loop index
    94 c
    95 C METHOD:
    96 C
    97 C See Emanuel, K. and M. Zivkovic-Rothman, 2000: Development and evaluation of a
    98 C       convective scheme for use in climate models.
    99 C
    100 C FILES: None
    101 C
    102 C INCLUDE FILES: None
    103 C
    104 C MAKEFILE: /a/ops/met/nogaps/src/sub/fcst/fcst159lib.mak
    105 C
    106 C..............................END PROLOGUE.............................
    107 c
    108 c
    109        USE dimphy
    110       implicit none
    111 c
    112 cym#include "dimensions.h"
    113 cym#include "dimphy.h"
    114 c
    115       integer len
    116       integer nd
    117       integer ndp1
    118       integer noff
    119       real t(len,nd)
    120       real q(len,nd)
    121       real qs(len,nd)
    122       real u(len,nd)
    123       real v(len,nd)
    124       real p(len,nd)
    125       real ph(len,ndp1)
    126       integer iflag(len)
    127       real ft(len,nd)
    128       real fq(len,nd)
    129       real fu(len,nd)
    130       real fv(len,nd)
    131       real precip(len)
    132       real cbmf(len)
    133       real Ma(len,nd)
    134       integer minorig
    135       real delt,cpd,cpv,cl,rv,rd,lv0,g
    136       real sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp
    137       real alpha,entp,coeffs,coeffr,omtrain,cu
    138 c
    139 !-------------------------------------------------------------------
    140 ! --- ARGUMENTS
    141 !-------------------------------------------------------------------
    142 ! --- On input:
    143 !
    144 !  t:   Array of absolute temperature (K) of dimension ND, with first
    145 !       index corresponding to lowest model level. Note that this array
    146 !       will be altered by the subroutine if dry convective adjustment
    147 !       occurs and if IPBL is not equal to 0.
    148 !
    149 !  q:   Array of specific humidity (gm/gm) of dimension ND, with first
    150 !       index corresponding to lowest model level. Must be defined
    151 !       at same grid levels as T. Note that this array will be altered
    152 !       if dry convective adjustment occurs and if IPBL is not equal to 0.
    153 !
    154 !  qs:  Array of saturation specific humidity of dimension ND, with first
    155 !       index corresponding to lowest model level. Must be defined
    156 !       at same grid levels as T. Note that this array will be altered
    157 !       if dry convective adjustment occurs and if IPBL is not equal to 0.
    158 !
    159 !  u:   Array of zonal wind velocity (m/s) of dimension ND, witth first
    160 !       index corresponding with the lowest model level. Defined at
    161 !       same levels as T. Note that this array will be altered if
    162 !       dry convective adjustment occurs and if IPBL is not equal to 0.
    163 !
    164 !  v:   Same as u but for meridional velocity.
    165 !
    166 !  tra: Array of passive tracer mixing ratio, of dimensions (ND,NTRA),
    167 !       where NTRA is the number of different tracers. If no
    168 !       convective tracer transport is needed, define a dummy
    169 !       input array of dimension (ND,1). Tracers are defined at
    170 !       same vertical levels as T. Note that this array will be altered
    171 !       if dry convective adjustment occurs and if IPBL is not equal to 0.
    172 !
    173 !  p:   Array of pressure (mb) of dimension ND, with first
    174 !       index corresponding to lowest model level. Must be defined
    175 !       at same grid levels as T.
    176 !
    177 !  ph:  Array of pressure (mb) of dimension ND+1, with first index
    178 !       corresponding to lowest level. These pressures are defined at
    179 !       levels intermediate between those of P, T, Q and QS. The first
    180 !       value of PH should be greater than (i.e. at a lower level than)
    181 !       the first value of the array P.
    182 !
    183 !  nl:  The maximum number of levels to which convection can penetrate, plus 1.
    184 !       NL MUST be less than or equal to ND-1.
    185 !
    186 !  delt: The model time step (sec) between calls to CONVECT
    187 !
    188 !----------------------------------------------------------------------------
    189 ! ---   On Output:
    190 !
    191 !  iflag: An output integer whose value denotes the following:
    192 !       VALUE   INTERPRETATION
    193 !       -----   --------------
    194 !         0     Moist convection occurs.
    195 !         1     Moist convection occurs, but a CFL condition
    196 !               on the subsidence warming is violated. This
    197 !               does not cause the scheme to terminate.
    198 !         2     Moist convection, but no precip because ep(inb) lt 0.0001
    199 !         3     No moist convection because new cbmf is 0 and old cbmf is 0.
    200 !         4     No moist convection; atmosphere is not
    201 !               unstable
    202 !         6     No moist convection because ihmin le minorig.
    203 !         7     No moist convection because unreasonable
    204 !               parcel level temperature or specific humidity.
    205 !         8     No moist convection: lifted condensation
    206 !               level is above the 200 mb level.
    207 !         9     No moist convection: cloud base is higher
    208 !               then the level NL-1.
    209 !
    210 !  ft:   Array of temperature tendency (K/s) of dimension ND, defined at same
    211 !        grid levels as T, Q, QS and P.
    212 !
    213 !  fq:   Array of specific humidity tendencies ((gm/gm)/s) of dimension ND,
    214 !        defined at same grid levels as T, Q, QS and P.
    215 !
    216 !  fu:   Array of forcing of zonal velocity (m/s^2) of dimension ND,
    217 !        defined at same grid levels as T.
    218 !
    219 !  fv:   Same as FU, but for forcing of meridional velocity.
    220 !
    221 !  ftra: Array of forcing of tracer content, in tracer mixing ratio per
    222 !        second, defined at same levels as T. Dimensioned (ND,NTRA).
    223 !
    224 !  precip: Scalar convective precipitation rate (mm/day).
    225 !
    226 !  wd:   A convective downdraft velocity scale. For use in surface
    227 !        flux parameterizations. See convect.ps file for details.
    228 !
    229 !  tprime: A convective downdraft temperature perturbation scale (K).
    230 !          For use in surface flux parameterizations. See convect.ps
    231 !          file for details.
    232 !
    233 !  qprime: A convective downdraft specific humidity
    234 !          perturbation scale (gm/gm).
    235 !          For use in surface flux parameterizations. See convect.ps
    236 !          file for details.
    237 !
    238 !  cbmf: The cloud base mass flux ((kg/m**2)/s). THIS SCALAR VALUE MUST
    239 !        BE STORED BY THE CALLING PROGRAM AND RETURNED TO CONVECT AT
    240 !        ITS NEXT CALL. That is, the value of CBMF must be "remembered"
    241 !        by the calling program between calls to CONVECT.
    242 !
    243 !  det:   Array of detrainment mass flux of dimension ND.
    244 !
    245 !-------------------------------------------------------------------
    246 c
    247 c  Local arrays
    248 c
    249       integer nl
    250       integer nlp
    251       integer nlm
    252       integer i,k,n
    253       real delti
    254       real rowl
    255       real clmcpv
    256       real clmcpd
    257       real cpdmcp
    258       real cpvmcpd
    259       real eps
    260       real epsi
    261       real epsim1
    262       real ginv
    263       real hrd
    264       real prccon1
    265       integer icbmax
    266       real lv(klon,klev)
    267       real cpn(klon,klev)
    268       real cpx(klon,klev)
    269       real tv(klon,klev)
    270       real gz(klon,klev)
    271       real hm(klon,klev)
    272       real h(klon,klev)
    273       real work(klon)
    274       integer ihmin(klon)
    275       integer nk(klon)
    276       real rh(klon)
    277       real chi(klon)
    278       real plcl(klon)
    279       integer icb(klon)
    280       real tnk(klon)
    281       real qnk(klon)
    282       real gznk(klon)
    283       real pnk(klon)
    284       real qsnk(klon)
    285       real ticb(klon)
    286       real gzicb(klon)
    287       real tp(klon,klev)
    288       real tvp(klon,klev)
    289       real clw(klon,klev)
    290 c
    291       real ah0(klon),cpp(klon)
    292       real tg,qg,s,alv,tc,ahg,denom,es,rg
    293 c
    294       integer ncum
    295       integer idcum(klon)
    296 c
    297       cpd=1005.7
    298       cpv=1870.0
    299       cl=4190.0
    300       rv=461.5
    301       rd=287.04
    302       lv0=2.501E6
    303       g=9.8
    304 C
    305 C   *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***
    306 C   ***  TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-        ***
    307 C   ***       CONVERSION THRESHOLD IS ASSUMED TO BE ZERO             ***
    308 C   ***     (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY            ***
    309 C   ***               BETWEEN 0 C AND TLCRIT)                        ***
    310 C   ***   ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT       ***
    311 C   ***                       FORMULATION                            ***
    312 C   ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
    313 C   ***  SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE       ***
    314 C   ***                        OF CLOUD                              ***
    315 C   ***        OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN       ***
    316 C   ***     OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW          ***
    317 C   ***  COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
    318 C   ***                          OF RAIN                             ***
    319 C   ***  COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
    320 C   ***                          OF SNOW                             ***
    321 C   ***     CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM      ***
    322 C   ***                         TRANSPORT                            ***
    323 C   ***    DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION    ***
    324 C   ***        A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC      ***
    325 C   ***    ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF    ***
    326 C   ***                 APPROACH TO QUASI-EQUILIBRIUM                ***
    327 C   ***   (THEIR STANDARD VALUES ARE  0.20 AND 0.1, RESPECTIVELY)    ***
    328 C   ***                   (DAMP MUST BE LESS THAN 1)                 ***
    329 c
    330       sigs=0.12
    331       sigd=0.05
    332       elcrit=0.0011
    333       tlcrit=-55.0
    334       omtsnow=5.5
    335       dtmax=0.9
    336       damp=0.1
    337       alpha=0.2
    338       entp=1.5
    339       coeffs=0.8
    340       coeffr=1.0
    341       omtrain=50.0
    342 c
    343       cu=0.70
    344       damp=0.1
    345 c
    346 c
    347 c Define nl, nlp, nlm, and delti
    348 c
    349       nl=nd-noff
    350       nlp=nl+1
    351       nlm=nl-1
    352       delti=1.0/delt
    353 !
    354 !-------------------------------------------------------------------
    355 ! --- SET CONSTANTS
    356 !-------------------------------------------------------------------
    357 !
    358       rowl=1000.0
    359       clmcpv=cl-cpv
    360       clmcpd=cl-cpd
    361       cpdmcp=cpd-cpv
    362       cpvmcpd=cpv-cpd
    363       eps=rd/rv
    364       epsi=1.0/eps
    365       epsim1=epsi-1.0
    366       ginv=1.0/g
    367       hrd=0.5*rd
    368       prccon1=86400.0*1000.0/(rowl*g)
    369 !
    370 ! dtmax is the maximum negative temperature perturbation.
    371 !
    372 !=====================================================================
    373 ! --- INITIALIZE OUTPUT ARRAYS AND PARAMETERS
    374 !=====================================================================
    375 !
    376       do 20 k=1,nd
    377         do 10 i=1,len
    378          ft(i,k)=0.0
    379          fq(i,k)=0.0
    380          fu(i,k)=0.0
    381          fv(i,k)=0.0
    382          tvp(i,k)=0.0
    383          tp(i,k)=0.0
    384          clw(i,k)=0.0
    385          gz(i,k) = 0.
    386  10     continue
    387  20   continue
    388       do 60 i=1,len
    389         precip(i)=0.0
    390         iflag(i)=0
    391  60   continue
    392 c
    393 !=====================================================================
    394 ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
    395 !=====================================================================
    396       do 110 k=1,nl+1
    397         do 100 i=1,len
    398           lv(i,k)= lv0-clmcpv*(t(i,k)-273.15)
    399           cpn(i,k)=cpd*(1.0-q(i,k))+cpv*q(i,k)
    400           cpx(i,k)=cpd*(1.0-q(i,k))+cl*q(i,k)
    401           tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1)
    402  100    continue
    403  110  continue
    404 c
    405 c gz = phi at the full levels (same as p).
    406 c
    407       do 120 i=1,len
    408         gz(i,1)=0.0
    409  120  continue
    410       do 140 k=2,nlp
    411         do 130 i=1,len
    412           gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
    413      &         *(p(i,k-1)-p(i,k))/ph(i,k)
    414  130    continue
    415  140  continue
    416 c
    417 c h  = phi + cpT (dry static energy).
    418 c hm = phi + cp(T-Tbase)+Lq
    419 c
    420       do 170 k=1,nlp
    421         do 160 i=1,len
    422           h(i,k)=gz(i,k)+cpn(i,k)*t(i,k)
    423           hm(i,k)=gz(i,k)+cpx(i,k)*(t(i,k)-t(i,1))+lv(i,k)*q(i,k)
    424  160    continue
    425  170  continue
    426 c
    427 !-------------------------------------------------------------------
    428 ! --- Find level of minimum moist static energy
    429 ! --- If level of minimum moist static energy coincides with
    430 ! --- or is lower than minimum allowable parcel origin level,
    431 ! --- set iflag to 6.
    432 !-------------------------------------------------------------------
    433       do 180 i=1,len
    434        work(i)=1.0e12
    435        ihmin(i)=nl
    436  180  continue
    437       do 200 k=2,nlp
    438         do 190 i=1,len
    439          if((hm(i,k).lt.work(i)).and.
    440      &      (hm(i,k).lt.hm(i,k-1)))then
    441            work(i)=hm(i,k)
    442            ihmin(i)=k
    443          endif
    444  190    continue
    445  200  continue
    446       do 210 i=1,len
    447         ihmin(i)=min(ihmin(i),nlm)
    448         if(ihmin(i).le.minorig)then
    449           iflag(i)=6
    450         endif
    451  210  continue
    452 c
    453 !-------------------------------------------------------------------
    454 ! --- Find that model level below the level of minimum moist static
    455 ! --- energy that has the maximum value of moist static energy
    456 !-------------------------------------------------------------------
    457  
    458       do 220 i=1,len
    459        work(i)=hm(i,minorig)
    460        nk(i)=minorig
    461  220  continue
    462       do 240 k=minorig+1,nl
    463         do 230 i=1,len
    464          if((hm(i,k).gt.work(i)).and.(k.le.ihmin(i)))then
    465            work(i)=hm(i,k)
    466            nk(i)=k
    467          endif
    468  230     continue
    469  240  continue
    470 !-------------------------------------------------------------------
    471 ! --- Check whether parcel level temperature and specific humidity
    472 ! --- are reasonable
    473 !-------------------------------------------------------------------
    474        do 250 i=1,len
    475        if(((t(i,nk(i)).lt.250.0).or.
    476      &      (q(i,nk(i)).le.0.0).or.
    477      &      (p(i,ihmin(i)).lt.400.0)).and.
    478      &      (iflag(i).eq.0))iflag(i)=7
    479  250   continue
    480 !-------------------------------------------------------------------
    481 ! --- Calculate lifted condensation level of air at parcel origin level
    482 ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
    483 !-------------------------------------------------------------------
    484        do 260 i=1,len
    485         tnk(i)=t(i,nk(i))
    486         qnk(i)=q(i,nk(i))
    487         gznk(i)=gz(i,nk(i))
    488         pnk(i)=p(i,nk(i))
    489         qsnk(i)=qs(i,nk(i))
    490 c
    491         rh(i)=qnk(i)/qsnk(i)
    492         rh(i)=min(1.0,rh(i))
    493         chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
    494         plcl(i)=pnk(i)*(rh(i)**chi(i))
    495         if(((plcl(i).lt.200.0).or.(plcl(i).ge.2000.0))
    496      &   .and.(iflag(i).eq.0))iflag(i)=8
    497  260   continue
    498 !-------------------------------------------------------------------
    499 ! --- Calculate first level above lcl (=icb)
    500 !-------------------------------------------------------------------
    501       do 270 i=1,len
    502        icb(i)=nlm
    503  270  continue
    504 c
    505       do 290 k=minorig,nl
    506         do 280 i=1,len
    507           if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i)))
    508      &    icb(i)=min(icb(i),k)
    509  280    continue
    510  290  continue
    511 c
    512       do 300 i=1,len
    513         if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
    514  300  continue
    515 c
    516 c Compute icbmax.
    517 c
    518       icbmax=2
    519       do 310 i=1,len
    520         icbmax=max(icbmax,icb(i))
    521  310  continue
    522 !
    523 !-------------------------------------------------------------------
    524 ! --- Calculates the lifted parcel virtual temperature at nk,
    525 ! --- the actual temperature, and the adiabatic
    526 ! --- liquid water content. The procedure is to solve the equation.
    527 !     cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
    528 !-------------------------------------------------------------------
    529 !
    530       do 320 i=1,len
    531         tnk(i)=t(i,nk(i))
    532         qnk(i)=q(i,nk(i))
    533         gznk(i)=gz(i,nk(i))
    534         ticb(i)=t(i,icb(i))
    535         gzicb(i)=gz(i,icb(i))
    536  320  continue
    537 c
    538 c   ***  Calculate certain parcel quantities, including static energy   ***
    539 c
    540       do 330 i=1,len
    541         ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i)
    542      &         +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i)
    543         cpp(i)=cpd*(1.-qnk(i))+qnk(i)*cpv
    544  330  continue
    545 c
    546 c   ***   Calculate lifted parcel quantities below cloud base   ***
    547 c
    548         do 350 k=minorig,icbmax-1
    549           do 340 i=1,len
    550            tp(i,k)=tnk(i)-(gz(i,k)-gznk(i))/cpp(i)
    551            tvp(i,k)=tp(i,k)*(1.+qnk(i)*epsi)
    552   340     continue
    553   350   continue
    554 c
    555 c    ***  Find lifted parcel quantities above cloud base    ***
    556 c
    557         do 360 i=1,len
    558          tg=ticb(i)
    559          qg=qs(i,icb(i))
    560          alv=lv0-clmcpv*(ticb(i)-273.15)
    561 c
    562 c First iteration.
    563 c
    564           s=cpd+alv*alv*qg/(rv*ticb(i)*ticb(i))
    565           s=1./s
    566           ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
    567           tg=tg+s*(ah0(i)-ahg)
    568           tg=max(tg,35.0)
    569           tc=tg-273.15
    570           denom=243.5+tc
    571           if(tc.ge.0.0)then
    572            es=6.112*exp(17.67*tc/denom)
    573           else
    574            es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
    575           endif
    576           qg=eps*es/(p(i,icb(i))-es*(1.-eps))
    577 c
    578 c Second iteration.
    579 c
    580           s=cpd+alv*alv*qg/(rv*ticb(i)*ticb(i))
    581           s=1./s
    582           ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
    583           tg=tg+s*(ah0(i)-ahg)
    584           tg=max(tg,35.0)
    585           tc=tg-273.15
    586           denom=243.5+tc
    587           if(tc.ge.0.0)then
    588            es=6.112*exp(17.67*tc/denom)
    589           else
    590            es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
    591           end if
    592           qg=eps*es/(p(i,icb(i))-es*(1.-eps))
    593 c
    594          alv=lv0-clmcpv*(ticb(i)-273.15)
    595          tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
    596      &   -gz(i,icb(i))-alv*qg)/cpd
    597          clw(i,icb(i))=qnk(i)-qg
    598          clw(i,icb(i))=max(0.0,clw(i,icb(i)))
    599          rg=qg/(1.-qnk(i))
    600          tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
    601   360   continue
    602 c
    603       do 380 k=minorig,icbmax
    604        do 370 i=1,len
    605          tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
    606  370   continue
    607  380  continue
    608 c
    609 !-------------------------------------------------------------------
    610 ! --- Test for instability.
    611 ! --- If there was no convection at last time step and parcel
    612 ! --- is stable at icb, then set iflag to 4.
    613 !-------------------------------------------------------------------
    614  
    615       do 390 i=1,len
    616         if((cbmf(i).eq.0.0) .and.(iflag(i).eq.0).and.
    617      &  (tvp(i,icb(i)).le.(tv(i,icb(i))-dtmax)))iflag(i)=4
    618  390  continue
    619  
    620 !=====================================================================
    621 ! --- IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY
    622 !=====================================================================
    623 c
    624       ncum=0
    625       do 400 i=1,len
    626         if(iflag(i).eq.0)then
    627            ncum=ncum+1
    628            idcum(ncum)=i
    629         endif
    630  400  continue
    631 c
    632 c Call convect2, which compresses the points and computes the heating,
    633 c moistening, velocity mixing, and precipiation.
    634 c
    635 c     print*,'cpd avant convect2 ',cpd
    636       if(ncum.gt.0)then
    637       call convect2(ncum,idcum,len,nd,ndp1,nl,minorig,
    638      &              nk,icb,
    639      &              t,q,qs,u,v,gz,tv,tp,tvp,clw,h,
    640      &              lv,cpn,p,ph,ft,fq,fu,fv,
    641      &              tnk,qnk,gznk,plcl,
    642      &              precip,cbmf,iflag,
    643      &              delt,cpd,cpv,cl,rv,rd,lv0,g,
    644      &              sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp,
    645      &              alpha,entp,coeffs,coeffr,omtrain,cu,Ma)
    646       endif
    647 c
    648       return
    649       end
     3
     4SUBROUTINE convect1(len, nd, ndp1, noff, minorig, t, q, qs, u, v, p, ph, &
     5    iflag, ft, fq, fu, fv, precip, cbmf, delt, ma)
     6  ! .............................START PROLOGUE............................
     7
     8  ! SCCS IDENTIFICATION:  @(#)convect1.f        1.1 04/21/00
     9  ! 19:40:52 /h/cm/library/nogaps4/src/sub/fcst/convect1.f_v
     10
     11  ! CONFIGURATION IDENTIFICATION:  None
     12
     13  ! MODULE NAME:  convect1
     14
     15  ! DESCRIPTION:
     16
     17  ! convect1     The Emanuel Cumulus Convection Scheme
     18
     19  ! CONTRACT NUMBER AND TITLE:  None
     20
     21  ! REFERENCES: Programmers  K. Emanuel (MIT), Timothy F. Hogan, M. Peng
     22  ! (NRL)
     23
     24  ! CLASSIFICATION:  Unclassified
     25
     26  ! RESTRICTIONS: None
     27
     28  ! COMPILER DEPENDENCIES: FORTRAN 77, FORTRAN 90
     29
     30  ! COMPILE OPTIONS: Fortran 77: -Zu -Wf"-ei -o aggress"
     31  ! Fortran 90: -O vector3,scalar3,task1,aggress,overindex  -ei -r 2
     32
     33  ! LIBRARIES OF RESIDENCE: /a/ops/lib/libfcst159.a
     34
     35  ! USAGE: call convect1(len,nd,noff,minorig,
     36  ! &                   t,q,qs,u,v,
     37  ! &                   p,ph,iflag,ft,
     38  ! &                   fq,fu,fv,precip,cbmf,delt)
     39
     40  ! PARAMETERS:
     41  ! Name            Type         Usage            Description
     42  ! ----------      ----------     -------  ----------------------------
     43
     44  ! len           Integer        Input        first (i) dimension
     45  ! nd            Integer        Input        vertical (k) dimension
     46  ! ndp1          Integer        Input        nd + 1
     47  ! noff          Integer        Input        integer limit for convection
     48  ! (nd-noff)
     49  ! minorig       Integer        Input        First level of convection
     50  ! t             Real           Input        temperature
     51  ! q             Real           Input        specific hum
     52  ! qs            Real           Input        sat specific hum
     53  ! u             Real           Input        u-wind
     54  ! v             Real           Input        v-wind
     55  ! p             Real           Input        full level pressure
     56  ! ph            Real           Input        half level pressure
     57  ! iflag         Integer        Output       iflag on latitude strip
     58  ! ft            Real           Output       temp tend
     59  ! fq            Real           Output       spec hum tend
     60  ! fu            Real           Output       u-wind tend
     61  ! fv            Real           Output       v-wind tend
     62  ! cbmf          Real           In/Out       cumulus mass flux
     63  ! delt          Real           Input        time step
     64  ! iflag         Integer        Output       integer flag for Emanuel
     65  ! conditions
     66
     67  ! COMMON BLOCKS:
     68  ! Block      Name     Type    Usage              Notes
     69  ! --------  --------   ----    ------   ------------------------
     70
     71  ! FILES: None
     72
     73  ! DATA BASES: None
     74
     75  ! NON-FILE INPUT/OUTPUT: None
     76
     77  ! ERROR CONDITIONS: None
     78
     79  ! ADDITIONAL COMMENTS: None
     80
     81  ! .................MAINTENANCE SECTION................................
     82
     83  ! MODULES CALLED:
     84  ! Name           Description
     85  ! convect2        Emanuel cumulus convection tendency calculations
     86  ! -------     ----------------------
     87  ! LOCAL VARIABLES AND
     88  ! STRUCTURES:
     89  ! Name     Type    Description
     90  ! -------  ------  -----------
     91  ! See Comments Below
     92
     93  ! i        Integer loop index
     94  ! k        Integer loop index
     95
     96  ! METHOD:
     97
     98  ! See Emanuel, K. and M. Zivkovic-Rothman, 2000: Development and evaluation
     99  ! of a
     100  ! convective scheme for use in climate models.
     101
     102  ! FILES: None
     103
     104  ! INCLUDE FILES: None
     105
     106  ! MAKEFILE: /a/ops/met/nogaps/src/sub/fcst/fcst159lib.mak
     107
     108  ! ..............................END PROLOGUE.............................
     109
     110
     111  USE dimphy
     112  IMPLICIT NONE
     113
     114  ! ym#include "dimensions.h"
     115  ! ym#include "dimphy.h"
     116
     117  INTEGER len
     118  INTEGER nd
     119  INTEGER ndp1
     120  INTEGER noff
     121  REAL t(len, nd)
     122  REAL q(len, nd)
     123  REAL qs(len, nd)
     124  REAL u(len, nd)
     125  REAL v(len, nd)
     126  REAL p(len, nd)
     127  REAL ph(len, ndp1)
     128  INTEGER iflag(len)
     129  REAL ft(len, nd)
     130  REAL fq(len, nd)
     131  REAL fu(len, nd)
     132  REAL fv(len, nd)
     133  REAL precip(len)
     134  REAL cbmf(len)
     135  REAL ma(len, nd)
     136  INTEGER minorig
     137  REAL delt, cpd, cpv, cl, rv, rd, lv0, g
     138  REAL sigs, sigd, elcrit, tlcrit, omtsnow, dtmax, damp
     139  REAL alpha, entp, coeffs, coeffr, omtrain, cu
     140
     141  ! -------------------------------------------------------------------
     142  ! --- ARGUMENTS
     143  ! -------------------------------------------------------------------
     144  ! --- On input:
     145
     146  ! t:   Array of absolute temperature (K) of dimension ND, with first
     147  ! index corresponding to lowest model level. Note that this array
     148  ! will be altered by the subroutine if dry convective adjustment
     149  ! occurs and if IPBL is not equal to 0.
     150
     151  ! q:   Array of specific humidity (gm/gm) of dimension ND, with first
     152  ! index corresponding to lowest model level. Must be defined
     153  ! at same grid levels as T. Note that this array will be altered
     154  ! if dry convective adjustment occurs and if IPBL is not equal to 0.
     155
     156  ! qs:  Array of saturation specific humidity of dimension ND, with first
     157  ! index corresponding to lowest model level. Must be defined
     158  ! at same grid levels as T. Note that this array will be altered
     159  ! if dry convective adjustment occurs and if IPBL is not equal to 0.
     160
     161  ! u:   Array of zonal wind velocity (m/s) of dimension ND, witth first
     162  ! index corresponding with the lowest model level. Defined at
     163  ! same levels as T. Note that this array will be altered if
     164  ! dry convective adjustment occurs and if IPBL is not equal to 0.
     165
     166  ! v:   Same as u but for meridional velocity.
     167
     168  ! tra: Array of passive tracer mixing ratio, of dimensions (ND,NTRA),
     169  ! where NTRA is the number of different tracers. If no
     170  ! convective tracer transport is needed, define a dummy
     171  ! input array of dimension (ND,1). Tracers are defined at
     172  ! same vertical levels as T. Note that this array will be altered
     173  ! if dry convective adjustment occurs and if IPBL is not equal to 0.
     174
     175  ! p:   Array of pressure (mb) of dimension ND, with first
     176  ! index corresponding to lowest model level. Must be defined
     177  ! at same grid levels as T.
     178
     179  ! ph:  Array of pressure (mb) of dimension ND+1, with first index
     180  ! corresponding to lowest level. These pressures are defined at
     181  ! levels intermediate between those of P, T, Q and QS. The first
     182  ! value of PH should be greater than (i.e. at a lower level than)
     183  ! the first value of the array P.
     184
     185  ! nl:  The maximum number of levels to which convection can penetrate, plus
     186  ! 1.
     187  ! NL MUST be less than or equal to ND-1.
     188
     189  ! delt: The model time step (sec) between calls to CONVECT
     190
     191  ! ----------------------------------------------------------------------------
     192  ! ---   On Output:
     193
     194  ! iflag: An output integer whose value denotes the following:
     195  ! VALUE   INTERPRETATION
     196  ! -----   --------------
     197  ! 0     Moist convection occurs.
     198  ! 1     Moist convection occurs, but a CFL condition
     199  ! on the subsidence warming is violated. This
     200  ! does not cause the scheme to terminate.
     201  ! 2     Moist convection, but no precip because ep(inb) lt 0.0001
     202  ! 3     No moist convection because new cbmf is 0 and old cbmf is 0.
     203  ! 4     No moist convection; atmosphere is not
     204  ! unstable
     205  ! 6     No moist convection because ihmin le minorig.
     206  ! 7     No moist convection because unreasonable
     207  ! parcel level temperature or specific humidity.
     208  ! 8     No moist convection: lifted condensation
     209  ! level is above the 200 mb level.
     210  ! 9     No moist convection: cloud base is higher
     211  ! then the level NL-1.
     212
     213  ! ft:   Array of temperature tendency (K/s) of dimension ND, defined at
     214  ! same
     215  ! grid levels as T, Q, QS and P.
     216
     217  ! fq:   Array of specific humidity tendencies ((gm/gm)/s) of dimension ND,
     218  ! defined at same grid levels as T, Q, QS and P.
     219
     220  ! fu:   Array of forcing of zonal velocity (m/s^2) of dimension ND,
     221  ! defined at same grid levels as T.
     222
     223  ! fv:   Same as FU, but for forcing of meridional velocity.
     224
     225  ! ftra: Array of forcing of tracer content, in tracer mixing ratio per
     226  ! second, defined at same levels as T. Dimensioned (ND,NTRA).
     227
     228  ! precip: Scalar convective precipitation rate (mm/day).
     229
     230  ! wd:   A convective downdraft velocity scale. For use in surface
     231  ! flux parameterizations. See convect.ps file for details.
     232
     233  ! tprime: A convective downdraft temperature perturbation scale (K).
     234  ! For use in surface flux parameterizations. See convect.ps
     235  ! file for details.
     236
     237  ! qprime: A convective downdraft specific humidity
     238  ! perturbation scale (gm/gm).
     239  ! For use in surface flux parameterizations. See convect.ps
     240  ! file for details.
     241
     242  ! cbmf: The cloud base mass flux ((kg/m**2)/s). THIS SCALAR VALUE MUST
     243  ! BE STORED BY THE CALLING PROGRAM AND RETURNED TO CONVECT AT
     244  ! ITS NEXT CALL. That is, the value of CBMF must be "remembered"
     245  ! by the calling program between calls to CONVECT.
     246
     247  ! det:   Array of detrainment mass flux of dimension ND.
     248
     249  ! -------------------------------------------------------------------
     250
     251  ! Local arrays
     252
     253  INTEGER nl
     254  INTEGER nlp
     255  INTEGER nlm
     256  INTEGER i, k, n
     257  REAL delti
     258  REAL rowl
     259  REAL clmcpv
     260  REAL clmcpd
     261  REAL cpdmcp
     262  REAL cpvmcpd
     263  REAL eps
     264  REAL epsi
     265  REAL epsim1
     266  REAL ginv
     267  REAL hrd
     268  REAL prccon1
     269  INTEGER icbmax
     270  REAL lv(klon, klev)
     271  REAL cpn(klon, klev)
     272  REAL cpx(klon, klev)
     273  REAL tv(klon, klev)
     274  REAL gz(klon, klev)
     275  REAL hm(klon, klev)
     276  REAL h(klon, klev)
     277  REAL work(klon)
     278  INTEGER ihmin(klon)
     279  INTEGER nk(klon)
     280  REAL rh(klon)
     281  REAL chi(klon)
     282  REAL plcl(klon)
     283  INTEGER icb(klon)
     284  REAL tnk(klon)
     285  REAL qnk(klon)
     286  REAL gznk(klon)
     287  REAL pnk(klon)
     288  REAL qsnk(klon)
     289  REAL ticb(klon)
     290  REAL gzicb(klon)
     291  REAL tp(klon, klev)
     292  REAL tvp(klon, klev)
     293  REAL clw(klon, klev)
     294
     295  REAL ah0(klon), cpp(klon)
     296  REAL tg, qg, s, alv, tc, ahg, denom, es, rg
     297
     298  INTEGER ncum
     299  INTEGER idcum(klon)
     300
     301  cpd = 1005.7
     302  cpv = 1870.0
     303  cl = 4190.0
     304  rv = 461.5
     305  rd = 287.04
     306  lv0 = 2.501E6
     307  g = 9.8
     308
     309  ! *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***
     310  ! ***  TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-        ***
     311  ! ***       CONVERSION THRESHOLD IS ASSUMED TO BE ZERO             ***
     312  ! ***     (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY            ***
     313  ! ***               BETWEEN 0 C AND TLCRIT)                        ***
     314  ! ***   ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT       ***
     315  ! ***                       FORMULATION                            ***
     316  ! ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
     317  ! ***  SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE       ***
     318  ! ***                        OF CLOUD                              ***
     319  ! ***        OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN       ***
     320  ! ***     OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW          ***
     321  ! ***  COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
     322  ! ***                          OF RAIN                             ***
     323  ! ***  COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
     324  ! ***                          OF SNOW                             ***
     325  ! ***     CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM      ***
     326  ! ***                         TRANSPORT                            ***
     327  ! ***    DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION    ***
     328  ! ***        A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC      ***
     329  ! ***    ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF    ***
     330  ! ***                 APPROACH TO QUASI-EQUILIBRIUM                ***
     331  ! ***   (THEIR STANDARD VALUES ARE  0.20 AND 0.1, RESPECTIVELY)    ***
     332  ! ***                   (DAMP MUST BE LESS THAN 1)                 ***
     333
     334  sigs = 0.12
     335  sigd = 0.05
     336  elcrit = 0.0011
     337  tlcrit = -55.0
     338  omtsnow = 5.5
     339  dtmax = 0.9
     340  damp = 0.1
     341  alpha = 0.2
     342  entp = 1.5
     343  coeffs = 0.8
     344  coeffr = 1.0
     345  omtrain = 50.0
     346
     347  cu = 0.70
     348  damp = 0.1
     349
     350
     351  ! Define nl, nlp, nlm, and delti
     352
     353  nl = nd - noff
     354  nlp = nl + 1
     355  nlm = nl - 1
     356  delti = 1.0/delt
     357
     358  ! -------------------------------------------------------------------
     359  ! --- SET CONSTANTS
     360  ! -------------------------------------------------------------------
     361
     362  rowl = 1000.0
     363  clmcpv = cl - cpv
     364  clmcpd = cl - cpd
     365  cpdmcp = cpd - cpv
     366  cpvmcpd = cpv - cpd
     367  eps = rd/rv
     368  epsi = 1.0/eps
     369  epsim1 = epsi - 1.0
     370  ginv = 1.0/g
     371  hrd = 0.5*rd
     372  prccon1 = 86400.0*1000.0/(rowl*g)
     373
     374  ! dtmax is the maximum negative temperature perturbation.
     375
     376  ! =====================================================================
     377  ! --- INITIALIZE OUTPUT ARRAYS AND PARAMETERS
     378  ! =====================================================================
     379
     380  DO k = 1, nd
     381    DO i = 1, len
     382      ft(i, k) = 0.0
     383      fq(i, k) = 0.0
     384      fu(i, k) = 0.0
     385      fv(i, k) = 0.0
     386      tvp(i, k) = 0.0
     387      tp(i, k) = 0.0
     388      clw(i, k) = 0.0
     389      gz(i, k) = 0.
     390    END DO
     391  END DO
     392  DO i = 1, len
     393    precip(i) = 0.0
     394    iflag(i) = 0
     395  END DO
     396
     397  ! =====================================================================
     398  ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
     399  ! =====================================================================
     400  DO k = 1, nl + 1
     401    DO i = 1, len
     402      lv(i, k) = lv0 - clmcpv*(t(i,k)-273.15)
     403      cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k)
     404      cpx(i, k) = cpd*(1.0-q(i,k)) + cl*q(i, k)
     405      tv(i, k) = t(i, k)*(1.0+q(i,k)*epsim1)
     406    END DO
     407  END DO
     408
     409  ! gz = phi at the full levels (same as p).
     410
     411  DO i = 1, len
     412    gz(i, 1) = 0.0
     413  END DO
     414  DO k = 2, nlp
     415    DO i = 1, len
     416      gz(i, k) = gz(i, k-1) + hrd*(tv(i,k-1)+tv(i,k))*(p(i,k-1)-p(i,k))/ph(i, &
     417        k)
     418    END DO
     419  END DO
     420
     421  ! h  = phi + cpT (dry static energy).
     422  ! hm = phi + cp(T-Tbase)+Lq
     423
     424  DO k = 1, nlp
     425    DO i = 1, len
     426      h(i, k) = gz(i, k) + cpn(i, k)*t(i, k)
     427      hm(i, k) = gz(i, k) + cpx(i, k)*(t(i,k)-t(i,1)) + lv(i, k)*q(i, k)
     428    END DO
     429  END DO
     430
     431  ! -------------------------------------------------------------------
     432  ! --- Find level of minimum moist static energy
     433  ! --- If level of minimum moist static energy coincides with
     434  ! --- or is lower than minimum allowable parcel origin level,
     435  ! --- set iflag to 6.
     436  ! -------------------------------------------------------------------
     437  DO i = 1, len
     438    work(i) = 1.0E12
     439    ihmin(i) = nl
     440  END DO
     441  DO k = 2, nlp
     442    DO i = 1, len
     443      IF ((hm(i,k)<work(i)) .AND. (hm(i,k)<hm(i,k-1))) THEN
     444        work(i) = hm(i, k)
     445        ihmin(i) = k
     446      END IF
     447    END DO
     448  END DO
     449  DO i = 1, len
     450    ihmin(i) = min(ihmin(i), nlm)
     451    IF (ihmin(i)<=minorig) THEN
     452      iflag(i) = 6
     453    END IF
     454  END DO
     455
     456  ! -------------------------------------------------------------------
     457  ! --- Find that model level below the level of minimum moist static
     458  ! --- energy that has the maximum value of moist static energy
     459  ! -------------------------------------------------------------------
     460
     461  DO i = 1, len
     462    work(i) = hm(i, minorig)
     463    nk(i) = minorig
     464  END DO
     465  DO k = minorig + 1, nl
     466    DO i = 1, len
     467      IF ((hm(i,k)>work(i)) .AND. (k<=ihmin(i))) THEN
     468        work(i) = hm(i, k)
     469        nk(i) = k
     470      END IF
     471    END DO
     472  END DO
     473  ! -------------------------------------------------------------------
     474  ! --- Check whether parcel level temperature and specific humidity
     475  ! --- are reasonable
     476  ! -------------------------------------------------------------------
     477  DO i = 1, len
     478    IF (((t(i,nk(i))<250.0) .OR. (q(i,nk(i))<=0.0) .OR. (p(i,ihmin(i))< &
     479      400.0)) .AND. (iflag(i)==0)) iflag(i) = 7
     480  END DO
     481  ! -------------------------------------------------------------------
     482  ! --- Calculate lifted condensation level of air at parcel origin level
     483  ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
     484  ! -------------------------------------------------------------------
     485  DO i = 1, len
     486    tnk(i) = t(i, nk(i))
     487    qnk(i) = q(i, nk(i))
     488    gznk(i) = gz(i, nk(i))
     489    pnk(i) = p(i, nk(i))
     490    qsnk(i) = qs(i, nk(i))
     491
     492    rh(i) = qnk(i)/qsnk(i)
     493    rh(i) = min(1.0, rh(i))
     494    chi(i) = tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
     495    plcl(i) = pnk(i)*(rh(i)**chi(i))
     496    IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag(i &
     497      ) = 8
     498  END DO
     499  ! -------------------------------------------------------------------
     500  ! --- Calculate first level above lcl (=icb)
     501  ! -------------------------------------------------------------------
     502  DO i = 1, len
     503    icb(i) = nlm
     504  END DO
     505
     506  DO k = minorig, nl
     507    DO i = 1, len
     508      IF ((k>=(nk(i)+1)) .AND. (p(i,k)<plcl(i))) icb(i) = min(icb(i), k)
     509    END DO
     510  END DO
     511
     512  DO i = 1, len
     513    IF ((icb(i)>=nlm) .AND. (iflag(i)==0)) iflag(i) = 9
     514  END DO
     515
     516  ! Compute icbmax.
     517
     518  icbmax = 2
     519  DO i = 1, len
     520    icbmax = max(icbmax, icb(i))
     521  END DO
     522
     523  ! -------------------------------------------------------------------
     524  ! --- Calculates the lifted parcel virtual temperature at nk,
     525  ! --- the actual temperature, and the adiabatic
     526  ! --- liquid water content. The procedure is to solve the equation.
     527  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
     528  ! -------------------------------------------------------------------
     529
     530  DO i = 1, len
     531    tnk(i) = t(i, nk(i))
     532    qnk(i) = q(i, nk(i))
     533    gznk(i) = gz(i, nk(i))
     534    ticb(i) = t(i, icb(i))
     535    gzicb(i) = gz(i, icb(i))
     536  END DO
     537
     538  ! ***  Calculate certain parcel quantities, including static energy   ***
     539
     540  DO i = 1, len
     541    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
     542      273.15)) + gznk(i)
     543    cpp(i) = cpd*(1.-qnk(i)) + qnk(i)*cpv
     544  END DO
     545
     546  ! ***   Calculate lifted parcel quantities below cloud base   ***
     547
     548  DO k = minorig, icbmax - 1
     549    DO i = 1, len
     550      tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))/cpp(i)
     551      tvp(i, k) = tp(i, k)*(1.+qnk(i)*epsi)
     552    END DO
     553  END DO
     554
     555  ! ***  Find lifted parcel quantities above cloud base    ***
     556
     557  DO i = 1, len
     558    tg = ticb(i)
     559    qg = qs(i, icb(i))
     560    alv = lv0 - clmcpv*(ticb(i)-273.15)
     561
     562    ! First iteration.
     563
     564    s = cpd + alv*alv*qg/(rv*ticb(i)*ticb(i))
     565    s = 1./s
     566    ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
     567    tg = tg + s*(ah0(i)-ahg)
     568    tg = max(tg, 35.0)
     569    tc = tg - 273.15
     570    denom = 243.5 + tc
     571    IF (tc>=0.0) THEN
     572      es = 6.112*exp(17.67*tc/denom)
     573    ELSE
     574      es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
     575    END IF
     576    qg = eps*es/(p(i,icb(i))-es*(1.-eps))
     577
     578    ! Second iteration.
     579
     580    s = cpd + alv*alv*qg/(rv*ticb(i)*ticb(i))
     581    s = 1./s
     582    ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
     583    tg = tg + s*(ah0(i)-ahg)
     584    tg = max(tg, 35.0)
     585    tc = tg - 273.15
     586    denom = 243.5 + tc
     587    IF (tc>=0.0) THEN
     588      es = 6.112*exp(17.67*tc/denom)
     589    ELSE
     590      es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
     591    END IF
     592    qg = eps*es/(p(i,icb(i))-es*(1.-eps))
     593
     594    alv = lv0 - clmcpv*(ticb(i)-273.15)
     595    tp(i, icb(i)) = (ah0(i)-(cl-cpd)*qnk(i)*ticb(i)-gz(i,icb(i))-alv*qg)/cpd
     596    clw(i, icb(i)) = qnk(i) - qg
     597    clw(i, icb(i)) = max(0.0, clw(i,icb(i)))
     598    rg = qg/(1.-qnk(i))
     599    tvp(i, icb(i)) = tp(i, icb(i))*(1.+rg*epsi)
     600  END DO
     601
     602  DO k = minorig, icbmax
     603    DO i = 1, len
     604      tvp(i, k) = tvp(i, k) - tp(i, k)*qnk(i)
     605    END DO
     606  END DO
     607
     608  ! -------------------------------------------------------------------
     609  ! --- Test for instability.
     610  ! --- If there was no convection at last time step and parcel
     611  ! --- is stable at icb, then set iflag to 4.
     612  ! -------------------------------------------------------------------
     613
     614  DO i = 1, len
     615    IF ((cbmf(i)==0.0) .AND. (iflag(i)==0) .AND. (tvp(i, &
     616      icb(i))<=(tv(i,icb(i))-dtmax))) iflag(i) = 4
     617  END DO
     618
     619  ! =====================================================================
     620  ! --- IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY
     621  ! =====================================================================
     622
     623  ncum = 0
     624  DO i = 1, len
     625    IF (iflag(i)==0) THEN
     626      ncum = ncum + 1
     627      idcum(ncum) = i
     628    END IF
     629  END DO
     630
     631  ! Call convect2, which compresses the points and computes the heating,
     632  ! moistening, velocity mixing, and precipiation.
     633
     634  ! print*,'cpd avant convect2 ',cpd
     635  IF (ncum>0) THEN
     636    CALL convect2(ncum, idcum, len, nd, ndp1, nl, minorig, nk, icb, t, q, qs, &
     637      u, v, gz, tv, tp, tvp, clw, h, lv, cpn, p, ph, ft, fq, fu, fv, tnk, &
     638      qnk, gznk, plcl, precip, cbmf, iflag, delt, cpd, cpv, cl, rv, rd, lv0, &
     639      g, sigs, sigd, elcrit, tlcrit, omtsnow, dtmax, damp, alpha, entp, &
     640      coeffs, coeffr, omtrain, cu, ma)
     641  END IF
     642
     643  RETURN
     644END SUBROUTINE convect1
Note: See TracChangeset for help on using the changeset viewer.