Ignore:
Timestamp:
Mar 5, 2014, 2:19:12 PM (10 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/cv_routines.F90

    r1988 r1992  
    1 !
     1
    22! $Id$
    3 !
    4       SUBROUTINE cv_param(nd)
    5       implicit none
    6 
    7 c------------------------------------------------------------
    8 c Set parameters for convectL
    9 c (includes microphysical parameters and parameters that
    10 c  control the rate of approach to quasi-equilibrium)
    11 c------------------------------------------------------------
    12 
    13 C   *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***
    14 C   ***  TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-        ***
    15 C   ***       CONVERSION THRESHOLD IS ASSUMED TO BE ZERO             ***
    16 C   ***     (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY            ***
    17 C   ***               BETWEEN 0 C AND TLCRIT)                        ***
    18 C   ***   ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT       ***
    19 C   ***                       FORMULATION                            ***
    20 C   ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
    21 C   ***  SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE       ***
    22 C   ***                        OF CLOUD                              ***
    23 C   ***        OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN       ***
    24 C   ***     OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW          ***
    25 C   ***  COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
    26 C   ***                          OF RAIN                             ***
    27 C   ***  COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
    28 C   ***                          OF SNOW                             ***
    29 C   ***     CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM      ***
    30 C   ***                         TRANSPORT                            ***
    31 C   ***    DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION    ***
    32 C   ***        A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC      ***
    33 C   ***    ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF    ***
    34 C   ***                 APPROACH TO QUASI-EQUILIBRIUM                ***
    35 C   ***   (THEIR STANDARD VALUES ARE  0.20 AND 0.1, RESPECTIVELY)    ***
    36 C   ***                   (DAMP MUST BE LESS THAN 1)                 ***
    37 
    38       include "cvparam.h"
    39       integer nd
    40       CHARACTER (LEN=20) :: modname='cv_routines'
    41       CHARACTER (LEN=80) :: abort_message
    42 
    43 c noff: integer limit for convection (nd-noff)
    44 c minorig: First level of convection
    45 
    46       noff = 2
    47       minorig = 2
    48 
    49       nl=nd-noff
    50       nlp=nl+1
    51       nlm=nl-1
    52 
    53       elcrit=0.0011
    54       tlcrit=-55.0
    55       entp=1.5
    56       sigs=0.12
    57       sigd=0.05
    58       omtrain=50.0
    59       omtsnow=5.5
    60       coeffr=1.0
    61       coeffs=0.8
    62       dtmax=0.9
    63 c
    64       cu=0.70
    65 c
    66       betad=10.0
    67 c
    68       damp=0.1
    69       alpha=0.2
    70 c
    71       delta=0.01  ! cld
    72 c
    73       return
    74       end
    75 
    76       SUBROUTINE cv_prelim(len,nd,ndp1,t,q,p,ph
    77      :                    ,lv,cpn,tv,gz,h,hm)
    78       implicit none
    79 
    80 !=====================================================================
    81 ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
    82 !=====================================================================
    83 
    84 c inputs:
    85       integer len, nd, ndp1
    86       real t(len,nd), q(len,nd), p(len,nd), ph(len,ndp1)
    87 
    88 c outputs:
    89       real lv(len,nd), cpn(len,nd), tv(len,nd)
    90       real gz(len,nd), h(len,nd), hm(len,nd)
    91 
    92 c local variables:
    93       integer k, i
    94       real cpx(len,nd)
    95 
    96       include "cvthermo.h"
    97       include "cvparam.h"
    98 
    99 
    100       do 110 k=1,nlp
    101         do 100 i=1,len
    102           lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
    103           cpn(i,k)=cpd*(1.0-q(i,k))+cpv*q(i,k)
    104           cpx(i,k)=cpd*(1.0-q(i,k))+cl*q(i,k)
    105           tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1)
    106  100    continue
    107  110  continue
    108 c
    109 c gz = phi at the full levels (same as p).
    110 c
    111       do 120 i=1,len
    112         gz(i,1)=0.0
    113  120  continue
    114       do 140 k=2,nlp
    115         do 130 i=1,len
    116           gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
    117      &         *(p(i,k-1)-p(i,k))/ph(i,k)
    118  130    continue
    119  140  continue
    120 c
    121 c h  = phi + cpT (dry static energy).
    122 c hm = phi + cp(T-Tbase)+Lq
    123 c
    124       do 170 k=1,nlp
    125         do 160 i=1,len
    126           h(i,k)=gz(i,k)+cpn(i,k)*t(i,k)
    127           hm(i,k)=gz(i,k)+cpx(i,k)*(t(i,k)-t(i,1))+lv(i,k)*q(i,k)
    128  160    continue
    129  170  continue
    130 
    131       return
    132       end
    133 
    134       SUBROUTINE cv_feed(len,nd,t,q,qs,p,hm,gz
    135      :                  ,nk,icb,icbmax,iflag,tnk,qnk,gznk,plcl)
    136       implicit none
    137 
    138 C================================================================
    139 C Purpose: CONVECTIVE FEED
    140 C================================================================
    141 
    142       include "cvparam.h"
    143 
    144 c inputs:
    145           integer len, nd
    146       real t(len,nd), q(len,nd), qs(len,nd), p(len,nd)
    147       real hm(len,nd), gz(len,nd)
    148 
    149 c outputs:
    150           integer iflag(len), nk(len), icb(len), icbmax
    151       real tnk(len), qnk(len), gznk(len), plcl(len)
    152 
    153 c local variables:
    154       integer i, k
    155       integer ihmin(len)
    156       real work(len)
    157       real pnk(len), qsnk(len), rh(len), chi(len)
    158 
    159 !-------------------------------------------------------------------
    160 ! --- Find level of minimum moist static energy
    161 ! --- If level of minimum moist static energy coincides with
    162 ! --- or is lower than minimum allowable parcel origin level,
    163 ! --- set iflag to 6.
    164 !-------------------------------------------------------------------
    165 
    166       do 180 i=1,len
    167        work(i)=1.0e12
    168        ihmin(i)=nl
    169  180  continue
    170       do 200 k=2,nlp
    171         do 190 i=1,len
    172          if((hm(i,k).lt.work(i)).and.
    173      &      (hm(i,k).lt.hm(i,k-1)))then
    174            work(i)=hm(i,k)
    175            ihmin(i)=k
    176          endif
    177  190    continue
    178  200  continue
    179       do 210 i=1,len
    180         ihmin(i)=min(ihmin(i),nlm)
    181         if(ihmin(i).le.minorig)then
    182           iflag(i)=6
    183         endif
    184  210  continue
    185 c
    186 !-------------------------------------------------------------------
    187 ! --- Find that model level below the level of minimum moist static
    188 ! --- energy that has the maximum value of moist static energy
    189 !-------------------------------------------------------------------
    190  
    191       do 220 i=1,len
    192        work(i)=hm(i,minorig)
    193        nk(i)=minorig
    194  220  continue
    195       do 240 k=minorig+1,nl
    196         do 230 i=1,len
    197          if((hm(i,k).gt.work(i)).and.(k.le.ihmin(i)))then
    198            work(i)=hm(i,k)
    199            nk(i)=k
    200          endif
    201  230     continue
    202  240  continue
    203 !-------------------------------------------------------------------
    204 ! --- Check whether parcel level temperature and specific humidity
    205 ! --- are reasonable
    206 !-------------------------------------------------------------------
    207        do 250 i=1,len
    208        if(((t(i,nk(i)).lt.250.0).or.
    209      &      (q(i,nk(i)).le.0.0).or.
    210      &      (p(i,ihmin(i)).lt.400.0)).and.
    211      &      (iflag(i).eq.0))iflag(i)=7
    212  250   continue
    213 !-------------------------------------------------------------------
    214 ! --- Calculate lifted condensation level of air at parcel origin level
    215 ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
    216 !-------------------------------------------------------------------
    217        do 260 i=1,len
    218         tnk(i)=t(i,nk(i))
    219         qnk(i)=q(i,nk(i))
    220         gznk(i)=gz(i,nk(i))
    221         pnk(i)=p(i,nk(i))
    222         qsnk(i)=qs(i,nk(i))
    223 c
    224         rh(i)=qnk(i)/qsnk(i)
    225         rh(i)=min(1.0,rh(i))
    226         chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
    227         plcl(i)=pnk(i)*(rh(i)**chi(i))
    228         if(((plcl(i).lt.200.0).or.(plcl(i).ge.2000.0))
    229      &   .and.(iflag(i).eq.0))iflag(i)=8
    230  260   continue
    231 !-------------------------------------------------------------------
    232 ! --- Calculate first level above lcl (=icb)
    233 !-------------------------------------------------------------------
    234       do 270 i=1,len
    235        icb(i)=nlm
    236  270  continue
    237 c
    238       do 290 k=minorig,nl
    239         do 280 i=1,len
    240           if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i)))
    241      &    icb(i)=min(icb(i),k)
    242  280    continue
    243  290  continue
    244 c
    245       do 300 i=1,len
    246         if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
    247  300  continue
    248 c
    249 c Compute icbmax.
    250 c
    251       icbmax=2
    252       do 310 i=1,len
    253         icbmax=max(icbmax,icb(i))
    254  310  continue
    255 
    256       return
    257       end
    258 
    259       SUBROUTINE cv_undilute1(len,nd,t,q,qs,gz,p,nk,icb,icbmax
    260      :                       ,tp,tvp,clw)
    261       implicit none
    262 
    263       include "cvthermo.h"
    264       include "cvparam.h"
    265 
    266 c inputs:
    267       integer len, nd
    268       integer nk(len), icb(len), icbmax
    269       real t(len,nd), q(len,nd), qs(len,nd), gz(len,nd)
    270       real p(len,nd)
    271 
    272 c outputs:
    273       real tp(len,nd), tvp(len,nd), clw(len,nd)
    274 
    275 c local variables:
    276       integer i, k
    277       real tg, qg, alv, s, ahg, tc, denom, es, rg
    278       real ah0(len), cpp(len)
    279       real tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)
    280 
    281 !-------------------------------------------------------------------
    282 ! --- Calculates the lifted parcel virtual temperature at nk,
    283 ! --- the actual temperature, and the adiabatic
    284 ! --- liquid water content. The procedure is to solve the equation.
    285 !     cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
    286 !-------------------------------------------------------------------
    287 
    288       do 320 i=1,len
    289         tnk(i)=t(i,nk(i))
    290         qnk(i)=q(i,nk(i))
    291         gznk(i)=gz(i,nk(i))
    292         ticb(i)=t(i,icb(i))
    293         gzicb(i)=gz(i,icb(i))
    294  320  continue
    295 c
    296 c   ***  Calculate certain parcel quantities, including static energy   ***
    297 c
    298       do 330 i=1,len
    299         ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i)
    300      &         +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i)
    301         cpp(i)=cpd*(1.-qnk(i))+qnk(i)*cpv
    302  330  continue
    303 c
    304 c   ***   Calculate lifted parcel quantities below cloud base   ***
    305 c
    306         do 350 k=minorig,icbmax-1
    307           do 340 i=1,len
    308            tp(i,k)=tnk(i)-(gz(i,k)-gznk(i))/cpp(i)
    309            tvp(i,k)=tp(i,k)*(1.+qnk(i)*epsi)
    310   340     continue
    311   350   continue
    312 c
    313 c    ***  Find lifted parcel quantities above cloud base    ***
    314 c
    315         do 360 i=1,len
    316          tg=ticb(i)
    317          qg=qs(i,icb(i))
    318          alv=lv0-clmcpv*(ticb(i)-t0)
    319 c
    320 c First iteration.
    321 c
    322           s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
    323           s=1./s
    324           ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
    325           tg=tg+s*(ah0(i)-ahg)
    326           tg=max(tg,35.0)
    327           tc=tg-t0
    328           denom=243.5+tc
    329           if(tc.ge.0.0)then
    330            es=6.112*exp(17.67*tc/denom)
    331           else
    332            es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
    333           endif
    334           qg=eps*es/(p(i,icb(i))-es*(1.-eps))
    335 c
    336 c Second iteration.
    337 c
    338           s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
    339           s=1./s
    340           ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
    341           tg=tg+s*(ah0(i)-ahg)
    342           tg=max(tg,35.0)
    343           tc=tg-t0
    344           denom=243.5+tc
    345           if(tc.ge.0.0)then
    346            es=6.112*exp(17.67*tc/denom)
    347           else
    348            es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
    349           end if
    350           qg=eps*es/(p(i,icb(i))-es*(1.-eps))
    351 c
    352          alv=lv0-clmcpv*(ticb(i)-273.15)
    353          tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
    354      &   -gz(i,icb(i))-alv*qg)/cpd
    355          clw(i,icb(i))=qnk(i)-qg
    356          clw(i,icb(i))=max(0.0,clw(i,icb(i)))
    357          rg=qg/(1.-qnk(i))
    358          tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
    359   360   continue
    360 c
    361       do 380 k=minorig,icbmax
    362        do 370 i=1,len
    363          tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
    364  370   continue
    365  380  continue
    366 c
    367       return
    368       end
    369 
    370       SUBROUTINE cv_trigger(len,nd,icb,cbmf,tv,tvp,iflag)
    371       implicit none
    372 
    373 !-------------------------------------------------------------------
    374 ! --- Test for instability.
    375 ! --- If there was no convection at last time step and parcel
    376 ! --- is stable at icb, then set iflag to 4.
    377 !-------------------------------------------------------------------
    378  
    379       include "cvparam.h"
    380 
    381 c inputs:
    382        integer len, nd, icb(len)
    383        real cbmf(len), tv(len,nd), tvp(len,nd)
    384 
    385 c outputs:
    386        integer iflag(len) ! also an input
    387 
    388 c local variables:
    389        integer i
    390 
    391 
    392       do 390 i=1,len
    393         if((cbmf(i).eq.0.0) .and.(iflag(i).eq.0).and.
    394      &  (tvp(i,icb(i)).le.(tv(i,icb(i))-dtmax)))iflag(i)=4
    395  390  continue
    396  
    397       return
    398       end
    399 
    400       SUBROUTINE cv_compress( len,nloc,ncum,nd
    401      :   ,iflag1,nk1,icb1
    402      :   ,cbmf1,plcl1,tnk1,qnk1,gznk1
    403      :   ,t1,q1,qs1,u1,v1,gz1
    404      :   ,h1,lv1,cpn1,p1,ph1,tv1,tp1,tvp1,clw1
    405      o   ,iflag,nk,icb
    406      o   ,cbmf,plcl,tnk,qnk,gznk
    407      o   ,t,q,qs,u,v,gz,h,lv,cpn,p,ph,tv,tp,tvp,clw
    408      o   ,dph          )
    409       implicit none
    410 
    411       include "cvparam.h"
    412 
    413 c inputs:
    414       integer len,ncum,nd,nloc
    415       integer iflag1(len),nk1(len),icb1(len)
    416       real cbmf1(len),plcl1(len),tnk1(len),qnk1(len),gznk1(len)
    417       real t1(len,nd),q1(len,nd),qs1(len,nd),u1(len,nd),v1(len,nd)
    418       real gz1(len,nd),h1(len,nd),lv1(len,nd),cpn1(len,nd)
    419       real p1(len,nd),ph1(len,nd+1),tv1(len,nd),tp1(len,nd)
    420       real tvp1(len,nd),clw1(len,nd)
    421 
    422 c outputs:
    423       integer iflag(nloc),nk(nloc),icb(nloc)
    424       real cbmf(nloc),plcl(nloc),tnk(nloc),qnk(nloc),gznk(nloc)
    425       real t(nloc,nd),q(nloc,nd),qs(nloc,nd),u(nloc,nd),v(nloc,nd)
    426       real gz(nloc,nd),h(nloc,nd),lv(nloc,nd),cpn(nloc,nd)
    427       real p(nloc,nd),ph(nloc,nd+1),tv(nloc,nd),tp(nloc,nd)
    428       real tvp(nloc,nd),clw(nloc,nd)
    429       real dph(nloc,nd)
    430 
    431 c local variables:
    432       integer i,k,nn
    433       CHARACTER (LEN=20) :: modname='cv_compress'
    434       CHARACTER (LEN=80) :: abort_message
    435 
    436       include 'iniprint.h'
    437 
    438 
    439       do 110 k=1,nl+1
    440        nn=0
    441       do 100 i=1,len
    442       if(iflag1(i).eq.0)then
    443         nn=nn+1
    444         t(nn,k)=t1(i,k)
    445         q(nn,k)=q1(i,k)
    446         qs(nn,k)=qs1(i,k)
    447         u(nn,k)=u1(i,k)
    448         v(nn,k)=v1(i,k)
    449         gz(nn,k)=gz1(i,k)
    450         h(nn,k)=h1(i,k)
    451         lv(nn,k)=lv1(i,k)
    452         cpn(nn,k)=cpn1(i,k)
    453         p(nn,k)=p1(i,k)
    454         ph(nn,k)=ph1(i,k)
    455         tv(nn,k)=tv1(i,k)
    456         tp(nn,k)=tp1(i,k)
    457         tvp(nn,k)=tvp1(i,k)
    458         clw(nn,k)=clw1(i,k)
    459       endif
    460  100    continue
    461  110  continue
    462 
    463       if (nn.ne.ncum) then
    464          write(lunout,*)'strange! nn not equal to ncum: ',nn,ncum
    465          abort_message = ''
    466          CALL abort_gcm (modname,abort_message,1)
    467       endif
    468 
    469       nn=0
    470       do 150 i=1,len
    471       if(iflag1(i).eq.0)then
    472       nn=nn+1
    473       cbmf(nn)=cbmf1(i)
    474       plcl(nn)=plcl1(i)
    475       tnk(nn)=tnk1(i)
    476       qnk(nn)=qnk1(i)
    477       gznk(nn)=gznk1(i)
    478       nk(nn)=nk1(i)
    479       icb(nn)=icb1(i)
    480       iflag(nn)=iflag1(i)
    481       endif
    482  150  continue
    483 
    484       do 170 k=1,nl
    485        do 160 i=1,ncum
    486         dph(i,k)=ph(i,k)-ph(i,k+1)
    487  160   continue
    488  170  continue
    489 
    490       return
    491       end
    492 
    493       SUBROUTINE cv_undilute2(nloc,ncum,nd,icb,nk
    494      :                       ,tnk,qnk,gznk,t,q,qs,gz
    495      :                       ,p,dph,h,tv,lv
    496      o                       ,inb,inb1,tp,tvp,clw,hp,ep,sigp,frac)
    497       implicit none
    498 
    499 C---------------------------------------------------------------------
    500 C Purpose:
    501 C     FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
    502 C     &
    503 C     COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
    504 C     FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
    505 C     &
    506 C     FIND THE LEVEL OF NEUTRAL BUOYANCY
    507 C---------------------------------------------------------------------
    508 
    509       include "cvthermo.h"
    510       include "cvparam.h"
    511 
    512 c inputs:
    513       integer ncum, nd, nloc
    514       integer icb(nloc), nk(nloc)
    515       real t(nloc,nd), q(nloc,nd), qs(nloc,nd), gz(nloc,nd)
    516       real p(nloc,nd), dph(nloc,nd)
    517       real tnk(nloc), qnk(nloc), gznk(nloc)
    518       real lv(nloc,nd), tv(nloc,nd), h(nloc,nd)
    519 
    520 c outputs:
    521       integer inb(nloc), inb1(nloc)
    522       real tp(nloc,nd), tvp(nloc,nd), clw(nloc,nd)
    523       real ep(nloc,nd), sigp(nloc,nd), hp(nloc,nd)
    524       real frac(nloc)
    525 
    526 c local variables:
    527       integer i, k
    528       real tg,qg,ahg,alv,s,tc,es,denom,rg,tca,elacrit
    529       real by, defrac
    530       real ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
    531       logical lcape(nloc)
    532 
    533 !=====================================================================
    534 ! --- SOME INITIALIZATIONS
    535 !=====================================================================
    536 
    537       do 170 k=1,nl
    538       do 160 i=1,ncum
    539        ep(i,k)=0.0
    540        sigp(i,k)=sigs
    541  160  continue
    542  170  continue
    543 
    544 !=====================================================================
    545 ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
    546 !=====================================================================
    547 c
    548 c ---       The procedure is to solve the equation.
    549 c              cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
    550 c
    551 c   ***  Calculate certain parcel quantities, including static energy   ***
    552 c
    553 c
    554       do 240 i=1,ncum
    555          ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i)
    556      &         +qnk(i)*(lv0-clmcpv*(tnk(i)-t0))+gznk(i)
    557  240  continue
    558 c
    559 c
    560 c    ***  Find lifted parcel quantities above cloud base    ***
    561 c
    562 c
    563         do 300 k=minorig+1,nl
    564           do 290 i=1,ncum
    565             if(k.ge.(icb(i)+1))then
    566               tg=t(i,k)
    567               qg=qs(i,k)
    568               alv=lv0-clmcpv*(t(i,k)-t0)
    569 c
    570 c First iteration.
    571 c
    572                s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
    573                s=1./s
    574                ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
    575                tg=tg+s*(ah0(i)-ahg)
    576                tg=max(tg,35.0)
    577                tc=tg-t0
    578                denom=243.5+tc
    579                if(tc.ge.0.0)then
    580                         es=6.112*exp(17.67*tc/denom)
    581                else
    582                         es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
    583                endif
    584                         qg=eps*es/(p(i,k)-es*(1.-eps))
    585 c
    586 c Second iteration.
    587 c
    588                s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
    589                s=1./s
    590                ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
    591                tg=tg+s*(ah0(i)-ahg)
    592                tg=max(tg,35.0)
    593                tc=tg-t0
    594                denom=243.5+tc
    595                if(tc.ge.0.0)then
    596                         es=6.112*exp(17.67*tc/denom)
    597                else
    598                         es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
    599                endif
    600                         qg=eps*es/(p(i,k)-es*(1.-eps))
    601 c
    602                alv=lv0-clmcpv*(t(i,k)-t0)
    603 c      print*,'cpd dans convect2 ',cpd
    604 c      print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
    605 c      print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
    606         tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
    607 c              if (.not.cpd.gt.1000.) then
    608 c                  print*,'CPD=',cpd
    609 c                  stop
    610 c              endif
    611                clw(i,k)=qnk(i)-qg
    612                clw(i,k)=max(0.0,clw(i,k))
    613                rg=qg/(1.-qnk(i))
    614                tvp(i,k)=tp(i,k)*(1.+rg*epsi)
    615             endif
    616   290     continue
    617   300   continue
    618 c
    619 !=====================================================================
    620 ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
    621 ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
    622 ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
    623 !=====================================================================
    624 c
    625       do 320 k=minorig+1,nl
    626         do 310 i=1,ncum
    627           if(k.ge.(nk(i)+1))then
    628             tca=tp(i,k)-t0
    629             if(tca.ge.0.0)then
    630               elacrit=elcrit
    631             else
    632               elacrit=elcrit*(1.0-tca/tlcrit)
    633             endif
    634             elacrit=max(elacrit,0.0)
    635             ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)
    636             ep(i,k)=max(ep(i,k),0.0 )
    637             ep(i,k)=min(ep(i,k),1.0 )
    638             sigp(i,k)=sigs
    639           endif
    640  310    continue
    641  320  continue
    642 c
    643 !=====================================================================
    644 ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
    645 ! --- VIRTUAL TEMPERATURE
    646 !=====================================================================
    647 c
    648       do 340 k=minorig+1,nl
    649         do 330 i=1,ncum
    650         if(k.ge.(icb(i)+1))then
    651           tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
    652 c         print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
    653 c         print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
    654         endif
    655  330    continue
    656  340  continue
    657       do 350 i=1,ncum
    658        tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
    659  350  continue
    660 c
    661 c=====================================================================
    662 c  --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
    663 c  --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
    664 c  --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
    665 c=====================================================================
    666 c
    667       do 510 i=1,ncum
    668         cape(i)=0.0
    669         capem(i)=0.0
    670         inb(i)=icb(i)+1
    671         inb1(i)=inb(i)
    672  510  continue
    673 c
    674 c Originial Code
    675 c
    676 c     do 530 k=minorig+1,nl-1
    677 c       do 520 i=1,ncum
    678 c         if(k.ge.(icb(i)+1))then
    679 c           by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
    680 c           byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
    681 c           cape(i)=cape(i)+by
    682 c           if(by.ge.0.0)inb1(i)=k+1
    683 c           if(cape(i).gt.0.0)then
    684 c             inb(i)=k+1
    685 c             capem(i)=cape(i)
    686 c           endif
    687 c         endif
    688 c520    continue
    689 c530  continue
    690 c     do 540 i=1,ncum
    691 c         byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
    692 c         cape(i)=capem(i)+byp
    693 c         defrac=capem(i)-cape(i)
    694 c         defrac=max(defrac,0.001)
    695 c         frac(i)=-cape(i)/defrac
    696 c         frac(i)=min(frac(i),1.0)
    697 c         frac(i)=max(frac(i),0.0)
    698 c540   continue
    699 c
    700 c K Emanuel fix
    701 c
    702 c     call zilch(byp,ncum)
    703 c     do 530 k=minorig+1,nl-1
    704 c       do 520 i=1,ncum
    705 c         if(k.ge.(icb(i)+1))then
    706 c           by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
    707 c           cape(i)=cape(i)+by
    708 c           if(by.ge.0.0)inb1(i)=k+1
    709 c           if(cape(i).gt.0.0)then
    710 c             inb(i)=k+1
    711 c             capem(i)=cape(i)
    712 c             byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
    713 c           endif
    714 c         endif
    715 c520    continue
    716 c530  continue
    717 c     do 540 i=1,ncum
    718 c         inb(i)=max(inb(i),inb1(i))
    719 c         cape(i)=capem(i)+byp(i)
    720 c         defrac=capem(i)-cape(i)
    721 c         defrac=max(defrac,0.001)
    722 c         frac(i)=-cape(i)/defrac
    723 c         frac(i)=min(frac(i),1.0)
    724 c         frac(i)=max(frac(i),0.0)
    725 c540   continue
    726 c
    727 c J Teixeira fix
    728 c
    729       call zilch(byp,ncum)
    730       do 515 i=1,ncum
    731         lcape(i)=.true.
    732  515  continue
    733       do 530 k=minorig+1,nl-1
    734         do 520 i=1,ncum
    735           if(cape(i).lt.0.0)lcape(i)=.false.
    736           if((k.ge.(icb(i)+1)).and.lcape(i))then
    737             by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
    738             byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
    739             cape(i)=cape(i)+by
    740             if(by.ge.0.0)inb1(i)=k+1
    741             if(cape(i).gt.0.0)then
    742               inb(i)=k+1
    743               capem(i)=cape(i)
    744             endif
    745           endif
    746  520    continue
    747  530  continue
    748       do 540 i=1,ncum
    749           cape(i)=capem(i)+byp(i)
    750           defrac=capem(i)-cape(i)
    751           defrac=max(defrac,0.001)
    752           frac(i)=-cape(i)/defrac
    753           frac(i)=min(frac(i),1.0)
    754           frac(i)=max(frac(i),0.0)
    755  540  continue
    756 c
    757 c=====================================================================
    758 c ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
    759 c=====================================================================
    760 c
    761 c initialization:
    762       do i=1,ncum*nlp
    763        hp(i,1)=h(i,1)
    764       enddo
    765 
    766       do 600 k=minorig+1,nl
    767         do 590 i=1,ncum
    768         if((k.ge.icb(i)).and.(k.le.inb(i)))then
    769           hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k)
    770         endif
    771  590    continue
    772  600  continue
    773 c
    774         return
    775         end
    776 c
    777       SUBROUTINE cv_closure(nloc,ncum,nd,nk,icb
    778      :                     ,tv,tvp,p,ph,dph,plcl,cpn
    779      :                     ,iflag,cbmf)
    780       implicit none
    781 
    782 c inputs:
    783       integer ncum, nd, nloc
    784       integer nk(nloc), icb(nloc)
    785       real tv(nloc,nd), tvp(nloc,nd), p(nloc,nd), dph(nloc,nd)
    786       real ph(nloc,nd+1) ! caution nd instead ndp1 to be consistent...
    787       real plcl(nloc), cpn(nloc,nd)
    788 
    789 c outputs:
    790       integer iflag(nloc)
    791       real cbmf(nloc) ! also an input
    792 
    793 c local variables:
    794       integer i, k, icbmax
    795       real dtpbl(nloc), dtmin(nloc), tvpplcl(nloc), tvaplcl(nloc)
    796       real work(nloc)
    797 
    798       include "cvthermo.h"
    799       include "cvparam.h"
    800 
    801 c-------------------------------------------------------------------
    802 c Compute icbmax.
    803 c-------------------------------------------------------------------
    804 
    805       icbmax=2
    806       do 230 i=1,ncum
    807        icbmax=max(icbmax,icb(i))
    808  230  continue
    809 
    810 c=====================================================================
    811 c ---  CALCULATE CLOUD BASE MASS FLUX
    812 c=====================================================================
    813 c
    814 c tvpplcl = parcel temperature lifted adiabatically from level
    815 c           icb-1 to the LCL.
    816 c tvaplcl = virtual temperature at the LCL.
    817 c
    818       do 610 i=1,ncum
    819         dtpbl(i)=0.0
    820         tvpplcl(i)=tvp(i,icb(i)-1)
    821      &  -rrd*tvp(i,icb(i)-1)*(p(i,icb(i)-1)-plcl(i))
    822      &  /(cpn(i,icb(i)-1)*p(i,icb(i)-1))
    823         tvaplcl(i)=tv(i,icb(i))
    824      &  +(tvp(i,icb(i))-tvp(i,icb(i)+1))*(plcl(i)-p(i,icb(i)))
    825      &  /(p(i,icb(i))-p(i,icb(i)+1))
    826  610  continue
    827 
    828 c-------------------------------------------------------------------
    829 c --- Interpolate difference between lifted parcel and
    830 c --- environmental temperatures to lifted condensation level
    831 c-------------------------------------------------------------------
    832 c
    833 c dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
    834 c
    835       do 630 k=minorig,icbmax
    836         do 620 i=1,ncum
    837         if((k.ge.nk(i)).and.(k.le.(icb(i)-1)))then
    838           dtpbl(i)=dtpbl(i)+(tvp(i,k)-tv(i,k))*dph(i,k)
    839         endif
    840  620    continue
    841  630  continue
    842       do 640 i=1,ncum
    843         dtpbl(i)=dtpbl(i)/(ph(i,nk(i))-ph(i,icb(i)))
    844         dtmin(i)=tvpplcl(i)-tvaplcl(i)+dtmax+dtpbl(i)
    845  640  continue
    846 c
    847 c-------------------------------------------------------------------
    848 c --- Adjust cloud base mass flux
    849 c-------------------------------------------------------------------
    850 c
    851       do 650 i=1,ncum
    852        work(i)=cbmf(i)
    853        cbmf(i)=max(0.0,(1.0-damp)*cbmf(i)+0.1*alpha*dtmin(i))
    854        if((work(i).eq.0.0).and.(cbmf(i).eq.0.0))then
    855          iflag(i)=3
    856        endif
    857  650  continue
    858 
    859        return
    860        end
    861 
    862       SUBROUTINE cv_mixing(nloc,ncum,nd,icb,nk,inb,inb1
    863      :                    ,ph,t,q,qs,u,v,h,lv,qnk
    864      :                    ,hp,tv,tvp,ep,clw,cbmf
    865      :                    ,m,ment,qent,uent,vent,nent,sij,elij)
    866       implicit none
    867 
    868       include "cvthermo.h"
    869       include "cvparam.h"
    870 
    871 c inputs:
    872       integer ncum, nd, nloc
    873       integer icb(nloc), inb(nloc), inb1(nloc), nk(nloc)
    874       real cbmf(nloc), qnk(nloc)
    875       real ph(nloc,nd+1)
    876       real t(nloc,nd), q(nloc,nd), qs(nloc,nd), lv(nloc,nd)
    877       real u(nloc,nd), v(nloc,nd), h(nloc,nd), hp(nloc,nd)
    878       real tv(nloc,nd), tvp(nloc,nd), ep(nloc,nd), clw(nloc,nd)
    879 
    880 c outputs:
    881       integer nent(nloc,nd)
    882       real m(nloc,nd), ment(nloc,nd,nd), qent(nloc,nd,nd)
    883       real uent(nloc,nd,nd), vent(nloc,nd,nd)
    884       real sij(nloc,nd,nd), elij(nloc,nd,nd)
    885 
    886 c local variables:
    887       integer i, j, k, ij
    888       integer num1, num2
    889       real dbo, qti, bf2, anum, denom, dei, altem, cwat, stemp
    890       real alt, qp1, smid, sjmin, sjmax, delp, delm
    891       real work(nloc), asij(nloc), smin(nloc), scrit(nloc)
    892       real bsum(nloc,nd)
    893       logical lwork(nloc)
    894 
    895 c=====================================================================
    896 c --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
    897 c=====================================================================
    898 c
    899         do 360 i=1,ncum*nlp
    900           nent(i,1)=0
    901           m(i,1)=0.0
    902  360    continue
    903 c
    904       do 400 k=1,nlp
    905        do 390 j=1,nlp
    906           do 385 i=1,ncum
    907             qent(i,k,j)=q(i,j)
    908             uent(i,k,j)=u(i,j)
    909             vent(i,k,j)=v(i,j)
    910             elij(i,k,j)=0.0
    911             ment(i,k,j)=0.0
    912             sij(i,k,j)=0.0
    913  385      continue
    914  390    continue
    915  400  continue
    916 c
    917 c-------------------------------------------------------------------
    918 c --- Calculate rates of mixing,  m(i)
    919 c-------------------------------------------------------------------
    920 c
    921       call zilch(work,ncum)
    922 c
    923       do 670 j=minorig+1,nl
    924         do 660 i=1,ncum
    925           if((j.ge.(icb(i)+1)).and.(j.le.inb(i)))then
    926              k=min(j,inb1(i))
    927              dbo=abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1))
    928      &       +entp*0.04*(ph(i,k)-ph(i,k+1))
    929              work(i)=work(i)+dbo
    930              m(i,j)=cbmf(i)*dbo
    931           endif
    932  660    continue
    933  670  continue
    934       do 690 k=minorig+1,nl
    935         do 680 i=1,ncum
    936           if((k.ge.(icb(i)+1)).and.(k.le.inb(i)))then
    937             m(i,k)=m(i,k)/work(i)
    938           endif
    939  680    continue
    940  690  continue
    941 c
    942 c
    943 c=====================================================================
    944 c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
    945 c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
    946 c --- FRACTION (sij)
    947 c=====================================================================
    948 c
    949 c
    950        do 750 i=minorig+1,nl
    951          do 710 j=minorig+1,nl
    952            do 700 ij=1,ncum
    953              if((i.ge.(icb(ij)+1)).and.(j.ge.icb(ij))
    954      &         .and.(i.le.inb(ij)).and.(j.le.inb(ij)))then
    955                qti=qnk(ij)-ep(ij,i)*clw(ij,i)
    956                bf2=1.+lv(ij,j)*lv(ij,j)*qs(ij,j)
    957      &         /(rrv*t(ij,j)*t(ij,j)*cpd)
    958                anum=h(ij,j)-hp(ij,i)+(cpv-cpd)*t(ij,j)*(qti-q(ij,j))
    959                denom=h(ij,i)-hp(ij,i)+(cpd-cpv)*(q(ij,i)-qti)*t(ij,j)
    960                dei=denom
    961                if(abs(dei).lt.0.01)dei=0.01
    962                sij(ij,i,j)=anum/dei
    963                sij(ij,i,i)=1.0
    964                altem=sij(ij,i,j)*q(ij,i)+(1.-sij(ij,i,j))*qti-qs(ij,j)
    965                altem=altem/bf2
    966                cwat=clw(ij,j)*(1.-ep(ij,j))
    967                stemp=sij(ij,i,j)
    968                if((stemp.lt.0.0.or.stemp.gt.1.0.or.
    969      1           altem.gt.cwat).and.j.gt.i)then
    970                  anum=anum-lv(ij,j)*(qti-qs(ij,j)-cwat*bf2)
    971                  denom=denom+lv(ij,j)*(q(ij,i)-qti)
    972                  if(abs(denom).lt.0.01)denom=0.01
    973                  sij(ij,i,j)=anum/denom
    974                  altem=sij(ij,i,j)*q(ij,i)+(1.-sij(ij,i,j))*qti-qs(ij,j)
    975                  altem=altem-(bf2-1.)*cwat
    976                endif
    977                if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then
    978                  qent(ij,i,j)=sij(ij,i,j)*q(ij,i)
    979      &                        +(1.-sij(ij,i,j))*qti
    980                  uent(ij,i,j)=sij(ij,i,j)*u(ij,i)
    981      &                        +(1.-sij(ij,i,j))*u(ij,nk(ij))
    982                  vent(ij,i,j)=sij(ij,i,j)*v(ij,i)
    983      &                        +(1.-sij(ij,i,j))*v(ij,nk(ij))
    984                  elij(ij,i,j)=altem
    985                  elij(ij,i,j)=max(0.0,elij(ij,i,j))
    986                  ment(ij,i,j)=m(ij,i)/(1.-sij(ij,i,j))
    987                  nent(ij,i)=nent(ij,i)+1
    988                endif
    989              sij(ij,i,j)=max(0.0,sij(ij,i,j))
    990              sij(ij,i,j)=min(1.0,sij(ij,i,j))
    991              endif
    992   700      continue
    993   710    continue
    994 c
    995 c   ***   If no air can entrain at level i assume that updraft detrains  ***
    996 c   ***   at that level and calculate detrained air flux and properties  ***
    997 c
    998            do 740 ij=1,ncum
    999              if((i.ge.(icb(ij)+1)).and.(i.le.inb(ij))
    1000      &       .and.(nent(ij,i).eq.0))then
    1001                ment(ij,i,i)=m(ij,i)
    1002                qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
    1003                uent(ij,i,i)=u(ij,nk(ij))
    1004                vent(ij,i,i)=v(ij,nk(ij))
    1005                elij(ij,i,i)=clw(ij,i)
    1006                sij(ij,i,i)=1.0
    1007              endif
    1008  740       continue
    1009  750   continue
    1010 c
    1011       do 770 i=1,ncum
    1012         sij(i,inb(i),inb(i))=1.0
    1013  770  continue
    1014 c
    1015 c=====================================================================
    1016 c   ---  NORMALIZE ENTRAINED AIR MASS FLUXES
    1017 c   ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
    1018 c=====================================================================
    1019 c
    1020        call zilch(bsum,ncum*nlp)
    1021        do 780 ij=1,ncum
    1022          lwork(ij)=.false.
    1023  780   continue
    1024        do 789 i=minorig+1,nl
    1025 c
    1026          num1=0
    1027          do 779 ij=1,ncum
    1028            if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))num1=num1+1
    1029  779     continue
    1030          if(num1.le.0)go to 789
    1031 c
    1032            do 781 ij=1,ncum
    1033              if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))then
    1034                 lwork(ij)=(nent(ij,i).ne.0)
    1035                 qp1=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
    1036                 anum=h(ij,i)-hp(ij,i)-lv(ij,i)*(qp1-qs(ij,i))
    1037                 denom=h(ij,i)-hp(ij,i)+lv(ij,i)*(q(ij,i)-qp1)
    1038                 if(abs(denom).lt.0.01)denom=0.01
    1039                 scrit(ij)=anum/denom
    1040                 alt=qp1-qs(ij,i)+scrit(ij)*(q(ij,i)-qp1)
    1041                 if(scrit(ij).lt.0.0.or.alt.lt.0.0)scrit(ij)=1.0
    1042                 asij(ij)=0.0
    1043                 smin(ij)=1.0
    1044              endif
    1045  781       continue
    1046          do 783 j=minorig,nl
    1047 c
    1048          num2=0
    1049          do 778 ij=1,ncum
    1050              if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
    1051      &       .and.(j.ge.icb(ij)).and.(j.le.inb(ij))
    1052      &       .and.lwork(ij))num2=num2+1
    1053  778     continue
    1054          if(num2.le.0)go to 783
    1055 c
    1056            do 782 ij=1,ncum
    1057              if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
    1058      &       .and.(j.ge.icb(ij)).and.(j.le.inb(ij)).and.lwork(ij))then
    1059                   if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then
    1060                     if(j.gt.i)then
    1061                       smid=min(sij(ij,i,j),scrit(ij))
    1062                       sjmax=smid
    1063                       sjmin=smid
    1064                         if(smid.lt.smin(ij)
    1065      &                  .and.sij(ij,i,j+1).lt.smid)then
    1066                           smin(ij)=smid
    1067                           sjmax=min(sij(ij,i,j+1),sij(ij,i,j),scrit(ij))
    1068                           sjmin=max(sij(ij,i,j-1),sij(ij,i,j))
    1069                           sjmin=min(sjmin,scrit(ij))
    1070                         endif
    1071                     else
    1072                       sjmax=max(sij(ij,i,j+1),scrit(ij))
    1073                       smid=max(sij(ij,i,j),scrit(ij))
    1074                       sjmin=0.0
    1075                       if(j.gt.1)sjmin=sij(ij,i,j-1)
    1076                       sjmin=max(sjmin,scrit(ij))
    1077                     endif
    1078                     delp=abs(sjmax-smid)
    1079                     delm=abs(sjmin-smid)
    1080                     asij(ij)=asij(ij)+(delp+delm)
    1081      &                           *(ph(ij,j)-ph(ij,j+1))
    1082                     ment(ij,i,j)=ment(ij,i,j)*(delp+delm)
    1083      &                           *(ph(ij,j)-ph(ij,j+1))
    1084                   endif
    1085               endif
    1086   782    continue
    1087   783    continue
    1088             do 784 ij=1,ncum
    1089             if((i.ge.icb(ij)+1).and.(i.le.inb(ij)).and.lwork(ij))then
    1090                asij(ij)=max(1.0e-21,asij(ij))
    1091                asij(ij)=1.0/asij(ij)
    1092                bsum(ij,i)=0.0
    1093             endif
    1094  784        continue
    1095             do 786 j=minorig,nl+1
    1096               do 785 ij=1,ncum
    1097                 if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
    1098      &          .and.(j.ge.icb(ij)).and.(j.le.inb(ij))
    1099      &          .and.lwork(ij))then
    1100                    ment(ij,i,j)=ment(ij,i,j)*asij(ij)
    1101                    bsum(ij,i)=bsum(ij,i)+ment(ij,i,j)
    1102                 endif
    1103  785     continue
    1104  786     continue
    1105              do 787 ij=1,ncum
    1106                if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
    1107      &         .and.(bsum(ij,i).lt.1.0e-18).and.lwork(ij))then
    1108                  nent(ij,i)=0
    1109                  ment(ij,i,i)=m(ij,i)
    1110                  qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
    1111                  uent(ij,i,i)=u(ij,nk(ij))
    1112                  vent(ij,i,i)=v(ij,nk(ij))
    1113                  elij(ij,i,i)=clw(ij,i)
    1114                  sij(ij,i,i)=1.0
    1115                endif
    1116   787        continue
    1117   789  continue
    1118 c
    1119        return
    1120        end
    1121 
    1122       SUBROUTINE cv_unsat(nloc,ncum,nd,inb,t,q,qs,gz,u,v,p,ph
    1123      :                  ,h,lv,ep,sigp,clw,m,ment,elij
    1124      :                  ,iflag,mp,qp,up,vp,wt,water,evap)
    1125       implicit none
    1126 
    1127 
    1128       include "cvthermo.h"
    1129       include "cvparam.h"
    1130 
    1131 c inputs:
    1132       integer ncum, nd, nloc
    1133       integer inb(nloc)
    1134       real t(nloc,nd), q(nloc,nd), qs(nloc,nd)
    1135       real gz(nloc,nd), u(nloc,nd), v(nloc,nd)
    1136       real p(nloc,nd), ph(nloc,nd+1), h(nloc,nd)
    1137       real lv(nloc,nd), ep(nloc,nd), sigp(nloc,nd), clw(nloc,nd)
    1138       real m(nloc,nd), ment(nloc,nd,nd), elij(nloc,nd,nd)
    1139 
    1140 c outputs:
    1141       integer iflag(nloc) ! also an input
    1142       real mp(nloc,nd), qp(nloc,nd), up(nloc,nd), vp(nloc,nd)
    1143       real water(nloc,nd), evap(nloc,nd), wt(nloc,nd)
    1144 
    1145 c local variables:
    1146       integer i,j,k,ij,num1
    1147       integer jtt(nloc)
    1148       real awat, coeff, qsm, afac, sigt, b6, c6, revap
    1149       real dhdp, fac, qstm, rat
    1150       real wdtrain(nloc)
    1151       logical lwork(nloc)
    1152 
    1153 c=====================================================================
    1154 c --- PRECIPITATING DOWNDRAFT CALCULATION
    1155 c=====================================================================
    1156 c
    1157 c Initializations:
    1158 c
    1159          do i = 1, ncum
    1160          do k = 1, nl+1
    1161           wt(i,k) = omtsnow
    1162           mp(i,k) = 0.0
    1163           evap(i,k) = 0.0
    1164           water(i,k) = 0.0
    1165          enddo
    1166          enddo
    1167 
    1168          do 420 i=1,ncum
    1169           qp(i,1)=q(i,1)
    1170           up(i,1)=u(i,1)
    1171           vp(i,1)=v(i,1)
    1172  420     continue
    1173 
    1174          do 440 k=2,nl+1
    1175          do 430 i=1,ncum
    1176           qp(i,k)=q(i,k-1)
    1177           up(i,k)=u(i,k-1)
    1178           vp(i,k)=v(i,k-1)
    1179  430     continue
    1180  440     continue
    1181 
    1182 
    1183 c   ***  Check whether ep(inb)=0, if so, skip precipitating    ***
    1184 c   ***             downdraft calculation                      ***
    1185 c
    1186 c
    1187 c   ***  Integrate liquid water equation to find condensed water   ***
    1188 c   ***                and condensed water flux                    ***
    1189 c
    1190 c
    1191       do 890 i=1,ncum
    1192         jtt(i)=2
    1193         if(ep(i,inb(i)).le.0.0001)iflag(i)=2
    1194         if(iflag(i).eq.0)then
    1195           lwork(i)=.true.
    1196         else
    1197           lwork(i)=.false.
    1198         endif
    1199  890  continue
    1200 c
    1201 c    ***                    Begin downdraft loop                    ***
    1202 c
    1203 c
    1204         call zilch(wdtrain,ncum)
    1205         do 899 i=nl+1,1,-1
    1206 c
    1207           num1=0
    1208           do 879 ij=1,ncum
    1209             if((i.le.inb(ij)).and.lwork(ij))num1=num1+1
    1210  879      continue
    1211           if(num1.le.0)go to 899
    1212 c
    1213 c
    1214 c    ***        Calculate detrained precipitation             ***
    1215 c
    1216           do 891 ij=1,ncum
    1217             if((i.le.inb(ij)).and.(lwork(ij)))then
    1218             wdtrain(ij)=g*ep(ij,i)*m(ij,i)*clw(ij,i)
    1219             endif
    1220  891      continue
    1221 c
    1222           if(i.gt.1)then
    1223             do 893 j=1,i-1
    1224               do 892 ij=1,ncum
    1225                 if((i.le.inb(ij)).and.(lwork(ij)))then
    1226                   awat=elij(ij,j,i)-(1.-ep(ij,i))*clw(ij,i)
    1227                   awat=max(0.0,awat)
    1228                   wdtrain(ij)=wdtrain(ij)+g*awat*ment(ij,j,i)
    1229                 endif
    1230  892          continue
    1231  893      continue
    1232           endif
    1233 c
    1234 c    ***    Find rain water and evaporation using provisional   ***
    1235 c    ***              estimates of qp(i)and qp(i-1)             ***
    1236 c
    1237 c
    1238 c  ***  Value of terminal velocity and coeffecient of evaporation for snow   ***
    1239 c
    1240           do 894 ij=1,ncum
    1241             if((i.le.inb(ij)).and.(lwork(ij)))then
    1242             coeff=coeffs
    1243             wt(ij,i)=omtsnow
    1244 c
    1245 c  ***  Value of terminal velocity and coeffecient of evaporation for rain   ***
    1246 c
    1247             if(t(ij,i).gt.273.0)then
    1248               coeff=coeffr
    1249               wt(ij,i)=omtrain
    1250             endif
    1251             qsm=0.5*(q(ij,i)+qp(ij,i+1))
    1252             afac=coeff*ph(ij,i)*(qs(ij,i)-qsm)
    1253      &       /(1.0e4+2.0e3*ph(ij,i)*qs(ij,i))
    1254             afac=max(afac,0.0)
    1255             sigt=sigp(ij,i)
    1256             sigt=max(0.0,sigt)
    1257             sigt=min(1.0,sigt)
    1258             b6=100.*(ph(ij,i)-ph(ij,i+1))*sigt*afac/wt(ij,i)
    1259             c6=(water(ij,i+1)*wt(ij,i+1)+wdtrain(ij)/sigd)/wt(ij,i)
    1260             revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
    1261             evap(ij,i)=sigt*afac*revap
    1262             water(ij,i)=revap*revap
    1263 c
    1264 c    ***  Calculate precipitating downdraft mass flux under     ***
    1265 c    ***              hydrostatic approximation                 ***
    1266 c
    1267             if(i.gt.1)then
    1268               dhdp=(h(ij,i)-h(ij,i-1))/(p(ij,i-1)-p(ij,i))
    1269               dhdp=max(dhdp,10.0)
    1270               mp(ij,i)=100.*ginv*lv(ij,i)*sigd*evap(ij,i)/dhdp
    1271               mp(ij,i)=max(mp(ij,i),0.0)
    1272 c
    1273 c   ***   Add small amount of inertia to downdraft              ***
    1274 c
    1275               fac=20.0/(ph(ij,i-1)-ph(ij,i))
    1276               mp(ij,i)=(fac*mp(ij,i+1)+mp(ij,i))/(1.+fac)
    1277 c
    1278 c    ***      Force mp to decrease linearly to zero                 ***
    1279 c    ***      between about 950 mb and the surface                  ***
    1280 c
    1281               if(p(ij,i).gt.(0.949*p(ij,1)))then
    1282                  jtt(ij)=max(jtt(ij),i)
    1283                  mp(ij,i)=mp(ij,jtt(ij))*(p(ij,1)-p(ij,i))
    1284      &           /(p(ij,1)-p(ij,jtt(ij)))
    1285               endif
    1286             endif
    1287 c
    1288 c    ***       Find mixing ratio of precipitating downdraft     ***
    1289 c
    1290             if(i.ne.inb(ij))then
    1291               if(i.eq.1)then
    1292                 qstm=qs(ij,1)
    1293               else
    1294                 qstm=qs(ij,i-1)
    1295               endif
    1296               if(mp(ij,i).gt.mp(ij,i+1))then
    1297                  rat=mp(ij,i+1)/mp(ij,i)
    1298                  qp(ij,i)=qp(ij,i+1)*rat+q(ij,i)*(1.0-rat)+100.*ginv*
    1299      &             sigd*(ph(ij,i)-ph(ij,i+1))*(evap(ij,i)/mp(ij,i))
    1300                  up(ij,i)=up(ij,i+1)*rat+u(ij,i)*(1.-rat)
    1301                  vp(ij,i)=vp(ij,i+1)*rat+v(ij,i)*(1.-rat)
    1302                else
    1303                  if(mp(ij,i+1).gt.0.0)then
    1304                    qp(ij,i)=(gz(ij,i+1)-gz(ij,i)
    1305      &               +qp(ij,i+1)*(lv(ij,i+1)+t(ij,i+1)
    1306      &               *(cl-cpd))+cpd*(t(ij,i+1)-t(ij,i)))
    1307      &               /(lv(ij,i)+t(ij,i)*(cl-cpd))
    1308                    up(ij,i)=up(ij,i+1)
    1309                    vp(ij,i)=vp(ij,i+1)
    1310                  endif
    1311               endif
    1312               qp(ij,i)=min(qp(ij,i),qstm)
    1313               qp(ij,i)=max(qp(ij,i),0.0)
    1314             endif
    1315             endif
    1316  894      continue
    1317  899    continue
    1318 c
    1319         return
    1320         end
    1321 
    1322       SUBROUTINE cv_yield(nloc,ncum,nd,nk,icb,inb,delt
    1323      :             ,t,q,u,v,gz,p,ph,h,hp,lv,cpn
    1324      :             ,ep,clw,frac,m,mp,qp,up,vp
    1325      :             ,wt,water,evap
    1326      :             ,ment,qent,uent,vent,nent,elij
    1327      :             ,tv,tvp
    1328      o             ,iflag,wd,qprime,tprime
    1329      o             ,precip,cbmf,ft,fq,fu,fv,Ma,qcondc)
    1330       implicit none
    1331 
    1332       include "cvthermo.h"
    1333       include "cvparam.h"
    1334 
    1335 c inputs
    1336       integer ncum, nd, nloc
    1337       integer nk(nloc), icb(nloc), inb(nloc)
    1338       integer nent(nloc,nd)
    1339       real delt
    1340       real t(nloc,nd), q(nloc,nd), u(nloc,nd), v(nloc,nd)
    1341       real gz(nloc,nd)
    1342       real p(nloc,nd), ph(nloc,nd+1), h(nloc,nd)
    1343       real hp(nloc,nd), lv(nloc,nd)
    1344       real cpn(nloc,nd), ep(nloc,nd), clw(nloc,nd), frac(nloc)
    1345       real m(nloc,nd), mp(nloc,nd), qp(nloc,nd)
    1346       real up(nloc,nd), vp(nloc,nd)
    1347       real wt(nloc,nd), water(nloc,nd), evap(nloc,nd)
    1348       real ment(nloc,nd,nd), qent(nloc,nd,nd), elij(nloc,nd,nd)
    1349       real uent(nloc,nd,nd), vent(nloc,nd,nd)
    1350       real tv(nloc,nd), tvp(nloc,nd)
    1351 
    1352 c outputs
    1353       integer iflag(nloc)  ! also an input
    1354       real cbmf(nloc)      ! also an input
    1355       real wd(nloc), tprime(nloc), qprime(nloc)
    1356       real precip(nloc)
    1357       real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
    1358       real Ma(nloc,nd)
    1359       real qcondc(nloc,nd)
    1360 
    1361 c local variables
    1362       integer i,j,ij,k,num1
    1363       real dpinv,cpinv,awat,fqold,ftold,fuold,fvold,delti
    1364       real work(nloc), am(nloc),amp1(nloc),ad(nloc)
    1365       real ents(nloc), uav(nloc),vav(nloc),lvcp(nloc,nd)
    1366       real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd) ! cld
    1367       real siga(nloc,nd), ax(nloc,nd), mac(nloc,nd)     ! cld
    1368 
    1369  
    1370 c -- initializations:
    1371 
    1372       delti = 1.0/delt
    1373 
    1374       do 160 i=1,ncum
    1375       precip(i)=0.0
    1376       wd(i)=0.0
    1377       tprime(i)=0.0
    1378       qprime(i)=0.0
    1379        do 170 k=1,nl+1
    1380         ft(i,k)=0.0
    1381         fu(i,k)=0.0
    1382         fv(i,k)=0.0
    1383         fq(i,k)=0.0
    1384         lvcp(i,k)=lv(i,k)/cpn(i,k)
    1385         qcondc(i,k)=0.0              ! cld
    1386         qcond(i,k)=0.0               ! cld
    1387         nqcond(i,k)=0.0              ! cld
    1388  170   continue
    1389  160  continue
    1390 
    1391 c
    1392 c   ***  Calculate surface precipitation in mm/day     ***
    1393 c
    1394         do 1190 i=1,ncum
    1395           if(iflag(i).le.1)then
    1396 cc            precip(i)=precip(i)+wt(i,1)*sigd*water(i,1)*3600.*24000.
    1397 cc     &                /(rowl*g)
    1398 cc            precip(i)=precip(i)*delt/86400.
    1399             precip(i) = wt(i,1)*sigd*water(i,1)*86400/g
    1400           endif
    1401  1190   continue
    1402 c
    1403 c
    1404 c   ***  Calculate downdraft velocity scale and surface temperature and  ***
    1405 c   ***                    water vapor fluctuations                      ***
    1406 c
    1407       do i=1,ncum
    1408        wd(i)=betad*abs(mp(i,icb(i)))*0.01*rrd*t(i,icb(i))
    1409      :           /(sigd*p(i,icb(i)))
    1410        qprime(i)=0.5*(qp(i,1)-q(i,1))
    1411        tprime(i)=lv0*qprime(i)/cpd
    1412       enddo
    1413 c
    1414 c   ***  Calculate tendencies of lowest level potential temperature  ***
    1415 c   ***                      and mixing ratio                        ***
    1416 c
    1417         do 1200 i=1,ncum
    1418           work(i)=0.01/(ph(i,1)-ph(i,2))
    1419           am(i)=0.0
    1420  1200   continue
    1421         do 1220 k=2,nl
    1422           do 1210 i=1,ncum
    1423             if((nk(i).eq.1).and.(k.le.inb(i)).and.(nk(i).eq.1))then
    1424               am(i)=am(i)+m(i,k)
    1425             endif
    1426  1210     continue
    1427  1220   continue
    1428         do 1240 i=1,ncum
    1429           if((g*work(i)*am(i)).ge.delti)iflag(i)=1
    1430           ft(i,1)=ft(i,1)+g*work(i)*am(i)*(t(i,2)-t(i,1)
    1431      &    +(gz(i,2)-gz(i,1))/cpn(i,1))
    1432           ft(i,1)=ft(i,1)-lvcp(i,1)*sigd*evap(i,1)
    1433           ft(i,1)=ft(i,1)+sigd*wt(i,2)*(cl-cpd)*water(i,2)*(t(i,2)
    1434      &     -t(i,1))*work(i)/cpn(i,1)
    1435           fq(i,1)=fq(i,1)+g*mp(i,2)*(qp(i,2)-q(i,1))*
    1436      &    work(i)+sigd*evap(i,1)
    1437           fq(i,1)=fq(i,1)+g*am(i)*(q(i,2)-q(i,1))*work(i)
    1438           fu(i,1)=fu(i,1)+g*work(i)*(mp(i,2)*(up(i,2)-u(i,1))
    1439      &    +am(i)*(u(i,2)-u(i,1)))
    1440           fv(i,1)=fv(i,1)+g*work(i)*(mp(i,2)*(vp(i,2)-v(i,1))
    1441      &    +am(i)*(v(i,2)-v(i,1)))
    1442  1240   continue
    1443         do 1260 j=2,nl
    1444            do 1250 i=1,ncum
    1445              if(j.le.inb(i))then
    1446                fq(i,1)=fq(i,1)
    1447      &                 +g*work(i)*ment(i,j,1)*(qent(i,j,1)-q(i,1))
    1448                fu(i,1)=fu(i,1)
    1449      &                 +g*work(i)*ment(i,j,1)*(uent(i,j,1)-u(i,1))
    1450                fv(i,1)=fv(i,1)
    1451      &                 +g*work(i)*ment(i,j,1)*(vent(i,j,1)-v(i,1))
    1452              endif
    1453  1250      continue
    1454  1260   continue
    1455 c
    1456 c   ***  Calculate tendencies of potential temperature and mixing ratio  ***
    1457 c   ***               at levels above the lowest level                   ***
    1458 c
    1459 c   ***  First find the net saturated updraft and downdraft mass fluxes  ***
    1460 c   ***                      through each level                          ***
    1461 c
    1462         do 1500 i=2,nl+1
    1463 c
    1464           num1=0
    1465           do 1265 ij=1,ncum
    1466             if(i.le.inb(ij))num1=num1+1
    1467  1265     continue
    1468           if(num1.le.0)go to 1500
    1469 c
    1470           call zilch(amp1,ncum)
    1471           call zilch(ad,ncum)
    1472 c
    1473           do 1280 k=i+1,nl+1
    1474             do 1270 ij=1,ncum
    1475               if((i.ge.nk(ij)).and.(i.le.inb(ij))
    1476      &            .and.(k.le.(inb(ij)+1)))then
    1477                 amp1(ij)=amp1(ij)+m(ij,k)
    1478               endif
    1479  1270         continue
    1480  1280     continue
    1481 c
    1482           do 1310 k=1,i
    1483             do 1300 j=i+1,nl+1
    1484                do 1290 ij=1,ncum
    1485                  if((j.le.(inb(ij)+1)).and.(i.le.inb(ij)))then
    1486                    amp1(ij)=amp1(ij)+ment(ij,k,j)
    1487                  endif
    1488  1290          continue
    1489  1300       continue
    1490  1310     continue
    1491           do 1340 k=1,i-1
    1492             do 1330 j=i,nl+1
    1493               do 1320 ij=1,ncum
    1494                 if((i.le.inb(ij)).and.(j.le.inb(ij)))then
    1495                    ad(ij)=ad(ij)+ment(ij,j,k)
    1496                 endif
    1497  1320         continue
    1498  1330       continue
    1499  1340     continue
    1500 c
    1501           do 1350 ij=1,ncum
    1502           if(i.le.inb(ij))then
    1503             dpinv=0.01/(ph(ij,i)-ph(ij,i+1))
    1504             cpinv=1.0/cpn(ij,i)
    1505 c
    1506             ft(ij,i)=ft(ij,i)
    1507      &       +g*dpinv*(amp1(ij)*(t(ij,i+1)-t(ij,i)
    1508      &       +(gz(ij,i+1)-gz(ij,i))*cpinv)
    1509      &       -ad(ij)*(t(ij,i)-t(ij,i-1)+(gz(ij,i)-gz(ij,i-1))*cpinv))
    1510      &       -sigd*lvcp(ij,i)*evap(ij,i)
    1511             ft(ij,i)=ft(ij,i)+g*dpinv*ment(ij,i,i)*(hp(ij,i)-h(ij,i)+
    1512      &        t(ij,i)*(cpv-cpd)*(q(ij,i)-qent(ij,i,i)))*cpinv
    1513             ft(ij,i)=ft(ij,i)+sigd*wt(ij,i+1)*(cl-cpd)*water(ij,i+1)*
    1514      &        (t(ij,i+1)-t(ij,i))*dpinv*cpinv
    1515             fq(ij,i)=fq(ij,i)+g*dpinv*(amp1(ij)*(q(ij,i+1)-q(ij,i))-
    1516      &        ad(ij)*(q(ij,i)-q(ij,i-1)))
    1517             fu(ij,i)=fu(ij,i)+g*dpinv*(amp1(ij)*(u(ij,i+1)-u(ij,i))-
    1518      &        ad(ij)*(u(ij,i)-u(ij,i-1)))
    1519             fv(ij,i)=fv(ij,i)+g*dpinv*(amp1(ij)*(v(ij,i+1)-v(ij,i))-
    1520      &        ad(ij)*(v(ij,i)-v(ij,i-1)))
    1521          endif
    1522  1350    continue
    1523          do 1370 k=1,i-1
    1524            do 1360 ij=1,ncum
    1525              if(i.le.inb(ij))then
    1526                awat=elij(ij,k,i)-(1.-ep(ij,i))*clw(ij,i)
    1527                awat=max(awat,0.0)
    1528                fq(ij,i)=fq(ij,i)
    1529      &         +g*dpinv*ment(ij,k,i)*(qent(ij,k,i)-awat-q(ij,i))
    1530                fu(ij,i)=fu(ij,i)
    1531      &         +g*dpinv*ment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
    1532                fv(ij,i)=fv(ij,i)
    1533      &         +g*dpinv*ment(ij,k,i)*(vent(ij,k,i)-v(ij,i))
    1534 c (saturated updrafts resulting from mixing)               ! cld
    1535                qcond(ij,i)=qcond(ij,i)+(elij(ij,k,i)-awat) ! cld
    1536                nqcond(ij,i)=nqcond(ij,i)+1.                ! cld
    1537              endif
    1538  1360      continue
    1539  1370    continue
    1540          do 1390 k=i,nl+1
    1541            do 1380 ij=1,ncum
    1542              if((i.le.inb(ij)).and.(k.le.inb(ij)))then
    1543                fq(ij,i)=fq(ij,i)
    1544      &                  +g*dpinv*ment(ij,k,i)*(qent(ij,k,i)-q(ij,i))
    1545                fu(ij,i)=fu(ij,i)
    1546      &                  +g*dpinv*ment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
    1547                fv(ij,i)=fv(ij,i)
    1548      &                  +g*dpinv*ment(ij,k,i)*(vent(ij,k,i)-v(ij,i))
    1549              endif
    1550  1380      continue
    1551  1390    continue
    1552           do 1400 ij=1,ncum
    1553            if(i.le.inb(ij))then
    1554              fq(ij,i)=fq(ij,i)
    1555      &                +sigd*evap(ij,i)+g*(mp(ij,i+1)*
    1556      &                (qp(ij,i+1)-q(ij,i))
    1557      &                -mp(ij,i)*(qp(ij,i)-q(ij,i-1)))*dpinv
    1558              fu(ij,i)=fu(ij,i)
    1559      &                +g*(mp(ij,i+1)*(up(ij,i+1)-u(ij,i))-mp(ij,i)*
    1560      &                (up(ij,i)-u(ij,i-1)))*dpinv
    1561              fv(ij,i)=fv(ij,i)
    1562      &               +g*(mp(ij,i+1)*(vp(ij,i+1)-v(ij,i))-mp(ij,i)*
    1563      &               (vp(ij,i)-v(ij,i-1)))*dpinv
    1564 C (saturated downdrafts resulting from mixing)               ! cld
    1565             do k=i+1,inb(ij)                                 ! cld
    1566              qcond(ij,i)=qcond(ij,i)+elij(ij,k,i)            ! cld
    1567              nqcond(ij,i)=nqcond(ij,i)+1.                    ! cld
    1568             enddo                                            ! cld
    1569 C (particular case: no detraining level is found)            ! cld
    1570             if (nent(ij,i).eq.0) then                        ! cld
    1571              qcond(ij,i)=qcond(ij,i)+(1.-ep(ij,i))*clw(ij,i) ! cld
    1572              nqcond(ij,i)=nqcond(ij,i)+1.                    ! cld
    1573             endif                                            ! cld
    1574             if (nqcond(ij,i).ne.0.) then                     ! cld
    1575              qcond(ij,i)=qcond(ij,i)/nqcond(ij,i)            ! cld
    1576             endif                                            ! cld
    1577            endif
    1578  1400     continue
    1579  1500   continue
    1580 c
    1581 c   *** Adjust tendencies at top of convection layer to reflect  ***
    1582 c   ***       actual position of the level zero cape             ***
    1583 c
    1584         do 503 ij=1,ncum
    1585         fqold=fq(ij,inb(ij))
    1586         fq(ij,inb(ij))=fq(ij,inb(ij))*(1.-frac(ij))
    1587         fq(ij,inb(ij)-1)=fq(ij,inb(ij)-1)
    1588      &   +frac(ij)*fqold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
    1589      1   (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*lv(ij,inb(ij))
    1590      &   /lv(ij,inb(ij)-1)
    1591         ftold=ft(ij,inb(ij))
    1592         ft(ij,inb(ij))=ft(ij,inb(ij))*(1.-frac(ij))
    1593         ft(ij,inb(ij)-1)=ft(ij,inb(ij)-1)
    1594      &   +frac(ij)*ftold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
    1595      1   (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*cpn(ij,inb(ij))
    1596      &   /cpn(ij,inb(ij)-1)
    1597         fuold=fu(ij,inb(ij))
    1598         fu(ij,inb(ij))=fu(ij,inb(ij))*(1.-frac(ij))
    1599         fu(ij,inb(ij)-1)=fu(ij,inb(ij)-1)
    1600      &   +frac(ij)*fuold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
    1601      1   (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
    1602         fvold=fv(ij,inb(ij))
    1603         fv(ij,inb(ij))=fv(ij,inb(ij))*(1.-frac(ij))
    1604         fv(ij,inb(ij)-1)=fv(ij,inb(ij)-1)
    1605      &  +frac(ij)*fvold*((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
    1606      1   (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
    1607  503    continue
    1608 c
    1609 c   ***   Very slightly adjust tendencies to force exact   ***
    1610 c   ***     enthalpy, momentum and tracer conservation     ***
    1611 c
    1612         do 682 ij=1,ncum
    1613         ents(ij)=0.0
    1614         uav(ij)=0.0
    1615         vav(ij)=0.0
    1616         do 681 i=1,inb(ij)
    1617          ents(ij)=ents(ij)
    1618      &  +(cpn(ij,i)*ft(ij,i)+lv(ij,i)*fq(ij,i))*(ph(ij,i)-ph(ij,i+1))   
    1619          uav(ij)=uav(ij)+fu(ij,i)*(ph(ij,i)-ph(ij,i+1))
    1620          vav(ij)=vav(ij)+fv(ij,i)*(ph(ij,i)-ph(ij,i+1))
    1621   681   continue
    1622   682   continue
    1623         do 683 ij=1,ncum
    1624         ents(ij)=ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
    1625         uav(ij)=uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
    1626         vav(ij)=vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
    1627  683    continue
    1628         do 642 ij=1,ncum
    1629         do 641 i=1,inb(ij)
    1630          ft(ij,i)=ft(ij,i)-ents(ij)/cpn(ij,i)
    1631          fu(ij,i)=(1.-cu)*(fu(ij,i)-uav(ij))
    1632          fv(ij,i)=(1.-cu)*(fv(ij,i)-vav(ij))
    1633   641   continue
    1634  642    continue
    1635 c
    1636         do 1810 k=1,nl+1
    1637           do 1800 i=1,ncum
    1638             if((q(i,k)+delt*fq(i,k)).lt.0.0)iflag(i)=10
    1639  1800     continue
    1640  1810   continue
    1641 c
    1642 c
    1643         do 1900 i=1,ncum
    1644           if(iflag(i).gt.2)then
    1645           precip(i)=0.0
    1646           cbmf(i)=0.0
    1647           endif
    1648  1900   continue
    1649         do 1920 k=1,nl
    1650          do 1910 i=1,ncum
    1651            if(iflag(i).gt.2)then
    1652              ft(i,k)=0.0
    1653              fq(i,k)=0.0
    1654              fu(i,k)=0.0
    1655              fv(i,k)=0.0
    1656              qcondc(i,k)=0.0                               ! cld
    1657            endif
    1658  1910    continue
    1659  1920   continue
    1660 
    1661         do k=1,nl+1
    1662         do i=1,ncum
    1663           Ma(i,k) = 0.
    1664         enddo
    1665         enddo
    1666         do k=nl,1,-1
    1667         do i=1,ncum
    1668           Ma(i,k) = Ma(i,k+1)+m(i,k)
    1669         enddo
    1670         enddo
    1671 
    1672 c
    1673 c   *** diagnose the in-cloud mixing ratio   ***            ! cld
    1674 c   ***           of condensed water         ***            ! cld
    1675 c                                                           ! cld
    1676       DO ij=1,ncum                                          ! cld   
    1677        do i=1,nd                                            ! cld
    1678         mac(ij,i)=0.0                                       ! cld   
    1679         wa(ij,i)=0.0                                        ! cld
    1680         siga(ij,i)=0.0                                      ! cld
    1681        enddo                                                ! cld
    1682        do i=nk(ij),inb(ij)                                  ! cld
    1683        do k=i+1,inb(ij)+1                                   ! cld
    1684         mac(ij,i)=mac(ij,i)+m(ij,k)                         ! cld
    1685        enddo                                                ! cld
    1686        enddo                                                ! cld
    1687        do i=icb(ij),inb(ij)-1                               ! cld
    1688         ax(ij,i)=0.                                         ! cld
    1689         do j=icb(ij),i                                      ! cld
    1690          ax(ij,i)=ax(ij,i)+rrd*(tvp(ij,j)-tv(ij,j))         ! cld   
    1691      :       *(ph(ij,j)-ph(ij,j+1))/p(ij,j)                 ! cld   
    1692         enddo                                               ! cld
    1693         if (ax(ij,i).gt.0.0) then                           ! cld   
    1694          wa(ij,i)=sqrt(2.*ax(ij,i))                         ! cld
    1695         endif                                               ! cld
    1696        enddo                                                ! cld
    1697        do i=1,nl                                            ! cld
    1698         if (wa(ij,i).gt.0.0)                                ! cld
    1699      :    siga(ij,i)=mac(ij,i)/wa(ij,i)                     ! cld   
    1700      :        *rrd*tvp(ij,i)/p(ij,i)/100./delta             ! cld   
    1701         siga(ij,i) = min(siga(ij,i),1.0)                    ! cld
    1702         qcondc(ij,i)=siga(ij,i)*clw(ij,i)*(1.-ep(ij,i))     ! cld   
    1703      :          + (1.-siga(ij,i))*qcond(ij,i)               ! cld   
    1704        enddo                                                ! cld
    1705       ENDDO                                                 ! cld   
    1706 
    1707         return
    1708         end
    1709 
    1710       SUBROUTINE cv_uncompress(nloc,len,ncum,nd,idcum
    1711      :         ,iflag
    1712      :         ,precip,cbmf
    1713      :         ,ft,fq,fu,fv
    1714      :         ,Ma,qcondc           
    1715      :         ,iflag1
    1716      :         ,precip1,cbmf1
    1717      :         ,ft1,fq1,fu1,fv1
    1718      :         ,Ma1,qcondc1           
    1719      :                               )
    1720       implicit none
    1721 
    1722       include "cvparam.h"
    1723 
    1724 c inputs:
    1725       integer len, ncum, nd, nloc
    1726       integer idcum(nloc)
    1727       integer iflag(nloc)
    1728       real precip(nloc), cbmf(nloc)
    1729       real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
    1730       real Ma(nloc,nd)
    1731       real qcondc(nloc,nd) !cld
    1732 
    1733 c outputs:
    1734       integer iflag1(len)
    1735       real precip1(len), cbmf1(len)
    1736       real ft1(len,nd), fq1(len,nd), fu1(len,nd), fv1(len,nd)
    1737       real Ma1(len,nd)
    1738       real qcondc1(len,nd) !cld
    1739 
    1740 c local variables:
    1741       integer i,k
    1742 
    1743         do 2000 i=1,ncum
    1744          precip1(idcum(i))=precip(i)
    1745          cbmf1(idcum(i))=cbmf(i)
    1746          iflag1(idcum(i))=iflag(i)
    1747  2000   continue
    1748 
    1749         do 2020 k=1,nl
    1750           do 2010 i=1,ncum
    1751             ft1(idcum(i),k)=ft(i,k)
    1752             fq1(idcum(i),k)=fq(i,k)
    1753             fu1(idcum(i),k)=fu(i,k)
    1754             fv1(idcum(i),k)=fv(i,k)
    1755             Ma1(idcum(i),k)=Ma(i,k)
    1756             qcondc1(idcum(i),k)=qcondc(i,k)
    1757  2010     continue
    1758  2020   continue
    1759 
    1760         return
    1761         end
    1762 
     3
     4SUBROUTINE cv_param(nd)
     5  IMPLICIT NONE
     6
     7  ! ------------------------------------------------------------
     8  ! Set parameters for convectL
     9  ! (includes microphysical parameters and parameters that
     10  ! control the rate of approach to quasi-equilibrium)
     11  ! ------------------------------------------------------------
     12
     13  ! *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***
     14  ! ***  TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-        ***
     15  ! ***       CONVERSION THRESHOLD IS ASSUMED TO BE ZERO             ***
     16  ! ***     (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY            ***
     17  ! ***               BETWEEN 0 C AND TLCRIT)                        ***
     18  ! ***   ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT       ***
     19  ! ***                       FORMULATION                            ***
     20  ! ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
     21  ! ***  SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE       ***
     22  ! ***                        OF CLOUD                              ***
     23  ! ***        OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN       ***
     24  ! ***     OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW          ***
     25  ! ***  COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
     26  ! ***                          OF RAIN                             ***
     27  ! ***  COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
     28  ! ***                          OF SNOW                             ***
     29  ! ***     CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM      ***
     30  ! ***                         TRANSPORT                            ***
     31  ! ***    DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION    ***
     32  ! ***        A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC      ***
     33  ! ***    ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF    ***
     34  ! ***                 APPROACH TO QUASI-EQUILIBRIUM                ***
     35  ! ***   (THEIR STANDARD VALUES ARE  0.20 AND 0.1, RESPECTIVELY)    ***
     36  ! ***                   (DAMP MUST BE LESS THAN 1)                 ***
     37
     38  include "cvparam.h"
     39  INTEGER nd
     40  CHARACTER (LEN=20) :: modname = 'cv_routines'
     41  CHARACTER (LEN=80) :: abort_message
     42
     43  ! noff: integer limit for convection (nd-noff)
     44  ! minorig: First level of convection
     45
     46  noff = 2
     47  minorig = 2
     48
     49  nl = nd - noff
     50  nlp = nl + 1
     51  nlm = nl - 1
     52
     53  elcrit = 0.0011
     54  tlcrit = -55.0
     55  entp = 1.5
     56  sigs = 0.12
     57  sigd = 0.05
     58  omtrain = 50.0
     59  omtsnow = 5.5
     60  coeffr = 1.0
     61  coeffs = 0.8
     62  dtmax = 0.9
     63
     64  cu = 0.70
     65
     66  betad = 10.0
     67
     68  damp = 0.1
     69  alpha = 0.2
     70
     71  delta = 0.01 ! cld
     72
     73  RETURN
     74END SUBROUTINE cv_param
     75
     76SUBROUTINE cv_prelim(len, nd, ndp1, t, q, p, ph, lv, cpn, tv, gz, h, hm)
     77  IMPLICIT NONE
     78
     79  ! =====================================================================
     80  ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
     81  ! =====================================================================
     82
     83  ! inputs:
     84  INTEGER len, nd, ndp1
     85  REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1)
     86
     87  ! outputs:
     88  REAL lv(len, nd), cpn(len, nd), tv(len, nd)
     89  REAL gz(len, nd), h(len, nd), hm(len, nd)
     90
     91  ! local variables:
     92  INTEGER k, i
     93  REAL cpx(len, nd)
     94
     95  include "cvthermo.h"
     96  include "cvparam.h"
     97
     98
     99  DO k = 1, nlp
     100    DO i = 1, len
     101      lv(i, k) = lv0 - clmcpv*(t(i,k)-t0)
     102      cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k)
     103      cpx(i, k) = cpd*(1.0-q(i,k)) + cl*q(i, k)
     104      tv(i, k) = t(i, k)*(1.0+q(i,k)*epsim1)
     105    END DO
     106  END DO
     107
     108  ! gz = phi at the full levels (same as p).
     109
     110  DO i = 1, len
     111    gz(i, 1) = 0.0
     112  END DO
     113  DO k = 2, nlp
     114    DO i = 1, len
     115      gz(i, k) = gz(i, k-1) + hrd*(tv(i,k-1)+tv(i,k))*(p(i,k-1)-p(i,k))/ph(i, &
     116        k)
     117    END DO
     118  END DO
     119
     120  ! h  = phi + cpT (dry static energy).
     121  ! hm = phi + cp(T-Tbase)+Lq
     122
     123  DO k = 1, nlp
     124    DO i = 1, len
     125      h(i, k) = gz(i, k) + cpn(i, k)*t(i, k)
     126      hm(i, k) = gz(i, k) + cpx(i, k)*(t(i,k)-t(i,1)) + lv(i, k)*q(i, k)
     127    END DO
     128  END DO
     129
     130  RETURN
     131END SUBROUTINE cv_prelim
     132
     133SUBROUTINE cv_feed(len, nd, t, q, qs, p, hm, gz, nk, icb, icbmax, iflag, tnk, &
     134    qnk, gznk, plcl)
     135  IMPLICIT NONE
     136
     137  ! ================================================================
     138  ! Purpose: CONVECTIVE FEED
     139  ! ================================================================
     140
     141  include "cvparam.h"
     142
     143  ! inputs:
     144  INTEGER len, nd
     145  REAL t(len, nd), q(len, nd), qs(len, nd), p(len, nd)
     146  REAL hm(len, nd), gz(len, nd)
     147
     148  ! outputs:
     149  INTEGER iflag(len), nk(len), icb(len), icbmax
     150  REAL tnk(len), qnk(len), gznk(len), plcl(len)
     151
     152  ! local variables:
     153  INTEGER i, k
     154  INTEGER ihmin(len)
     155  REAL work(len)
     156  REAL pnk(len), qsnk(len), rh(len), chi(len)
     157
     158  ! -------------------------------------------------------------------
     159  ! --- Find level of minimum moist static energy
     160  ! --- If level of minimum moist static energy coincides with
     161  ! --- or is lower than minimum allowable parcel origin level,
     162  ! --- set iflag to 6.
     163  ! -------------------------------------------------------------------
     164
     165  DO i = 1, len
     166    work(i) = 1.0E12
     167    ihmin(i) = nl
     168  END DO
     169  DO k = 2, nlp
     170    DO i = 1, len
     171      IF ((hm(i,k)<work(i)) .AND. (hm(i,k)<hm(i,k-1))) THEN
     172        work(i) = hm(i, k)
     173        ihmin(i) = k
     174      END IF
     175    END DO
     176  END DO
     177  DO i = 1, len
     178    ihmin(i) = min(ihmin(i), nlm)
     179    IF (ihmin(i)<=minorig) THEN
     180      iflag(i) = 6
     181    END IF
     182  END DO
     183
     184  ! -------------------------------------------------------------------
     185  ! --- Find that model level below the level of minimum moist static
     186  ! --- energy that has the maximum value of moist static energy
     187  ! -------------------------------------------------------------------
     188
     189  DO i = 1, len
     190    work(i) = hm(i, minorig)
     191    nk(i) = minorig
     192  END DO
     193  DO k = minorig + 1, nl
     194    DO i = 1, len
     195      IF ((hm(i,k)>work(i)) .AND. (k<=ihmin(i))) THEN
     196        work(i) = hm(i, k)
     197        nk(i) = k
     198      END IF
     199    END DO
     200  END DO
     201  ! -------------------------------------------------------------------
     202  ! --- Check whether parcel level temperature and specific humidity
     203  ! --- are reasonable
     204  ! -------------------------------------------------------------------
     205  DO i = 1, len
     206    IF (((t(i,nk(i))<250.0) .OR. (q(i,nk(i))<=0.0) .OR. (p(i,ihmin(i))< &
     207      400.0)) .AND. (iflag(i)==0)) iflag(i) = 7
     208  END DO
     209  ! -------------------------------------------------------------------
     210  ! --- Calculate lifted condensation level of air at parcel origin level
     211  ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
     212  ! -------------------------------------------------------------------
     213  DO i = 1, len
     214    tnk(i) = t(i, nk(i))
     215    qnk(i) = q(i, nk(i))
     216    gznk(i) = gz(i, nk(i))
     217    pnk(i) = p(i, nk(i))
     218    qsnk(i) = qs(i, nk(i))
     219
     220    rh(i) = qnk(i)/qsnk(i)
     221    rh(i) = min(1.0, rh(i))
     222    chi(i) = tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
     223    plcl(i) = pnk(i)*(rh(i)**chi(i))
     224    IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag(i &
     225      ) = 8
     226  END DO
     227  ! -------------------------------------------------------------------
     228  ! --- Calculate first level above lcl (=icb)
     229  ! -------------------------------------------------------------------
     230  DO i = 1, len
     231    icb(i) = nlm
     232  END DO
     233
     234  DO k = minorig, nl
     235    DO i = 1, len
     236      IF ((k>=(nk(i)+1)) .AND. (p(i,k)<plcl(i))) icb(i) = min(icb(i), k)
     237    END DO
     238  END DO
     239
     240  DO i = 1, len
     241    IF ((icb(i)>=nlm) .AND. (iflag(i)==0)) iflag(i) = 9
     242  END DO
     243
     244  ! Compute icbmax.
     245
     246  icbmax = 2
     247  DO i = 1, len
     248    icbmax = max(icbmax, icb(i))
     249  END DO
     250
     251  RETURN
     252END SUBROUTINE cv_feed
     253
     254SUBROUTINE cv_undilute1(len, nd, t, q, qs, gz, p, nk, icb, icbmax, tp, tvp, &
     255    clw)
     256  IMPLICIT NONE
     257
     258  include "cvthermo.h"
     259  include "cvparam.h"
     260
     261  ! inputs:
     262  INTEGER len, nd
     263  INTEGER nk(len), icb(len), icbmax
     264  REAL t(len, nd), q(len, nd), qs(len, nd), gz(len, nd)
     265  REAL p(len, nd)
     266
     267  ! outputs:
     268  REAL tp(len, nd), tvp(len, nd), clw(len, nd)
     269
     270  ! local variables:
     271  INTEGER i, k
     272  REAL tg, qg, alv, s, ahg, tc, denom, es, rg
     273  REAL ah0(len), cpp(len)
     274  REAL tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)
     275
     276  ! -------------------------------------------------------------------
     277  ! --- Calculates the lifted parcel virtual temperature at nk,
     278  ! --- the actual temperature, and the adiabatic
     279  ! --- liquid water content. The procedure is to solve the equation.
     280  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
     281  ! -------------------------------------------------------------------
     282
     283  DO i = 1, len
     284    tnk(i) = t(i, nk(i))
     285    qnk(i) = q(i, nk(i))
     286    gznk(i) = gz(i, nk(i))
     287    ticb(i) = t(i, icb(i))
     288    gzicb(i) = gz(i, icb(i))
     289  END DO
     290
     291  ! ***  Calculate certain parcel quantities, including static energy   ***
     292
     293  DO i = 1, len
     294    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
     295      273.15)) + gznk(i)
     296    cpp(i) = cpd*(1.-qnk(i)) + qnk(i)*cpv
     297  END DO
     298
     299  ! ***   Calculate lifted parcel quantities below cloud base   ***
     300
     301  DO k = minorig, icbmax - 1
     302    DO i = 1, len
     303      tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))/cpp(i)
     304      tvp(i, k) = tp(i, k)*(1.+qnk(i)*epsi)
     305    END DO
     306  END DO
     307
     308  ! ***  Find lifted parcel quantities above cloud base    ***
     309
     310  DO i = 1, len
     311    tg = ticb(i)
     312    qg = qs(i, icb(i))
     313    alv = lv0 - clmcpv*(ticb(i)-t0)
     314
     315    ! First iteration.
     316
     317    s = cpd + alv*alv*qg/(rrv*ticb(i)*ticb(i))
     318    s = 1./s
     319    ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
     320    tg = tg + s*(ah0(i)-ahg)
     321    tg = max(tg, 35.0)
     322    tc = tg - t0
     323    denom = 243.5 + tc
     324    IF (tc>=0.0) THEN
     325      es = 6.112*exp(17.67*tc/denom)
     326    ELSE
     327      es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
     328    END IF
     329    qg = eps*es/(p(i,icb(i))-es*(1.-eps))
     330
     331    ! Second iteration.
     332
     333    s = cpd + alv*alv*qg/(rrv*ticb(i)*ticb(i))
     334    s = 1./s
     335    ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
     336    tg = tg + s*(ah0(i)-ahg)
     337    tg = max(tg, 35.0)
     338    tc = tg - t0
     339    denom = 243.5 + tc
     340    IF (tc>=0.0) THEN
     341      es = 6.112*exp(17.67*tc/denom)
     342    ELSE
     343      es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
     344    END IF
     345    qg = eps*es/(p(i,icb(i))-es*(1.-eps))
     346
     347    alv = lv0 - clmcpv*(ticb(i)-273.15)
     348    tp(i, icb(i)) = (ah0(i)-(cl-cpd)*qnk(i)*ticb(i)-gz(i,icb(i))-alv*qg)/cpd
     349    clw(i, icb(i)) = qnk(i) - qg
     350    clw(i, icb(i)) = max(0.0, clw(i,icb(i)))
     351    rg = qg/(1.-qnk(i))
     352    tvp(i, icb(i)) = tp(i, icb(i))*(1.+rg*epsi)
     353  END DO
     354
     355  DO k = minorig, icbmax
     356    DO i = 1, len
     357      tvp(i, k) = tvp(i, k) - tp(i, k)*qnk(i)
     358    END DO
     359  END DO
     360
     361  RETURN
     362END SUBROUTINE cv_undilute1
     363
     364SUBROUTINE cv_trigger(len, nd, icb, cbmf, tv, tvp, iflag)
     365  IMPLICIT NONE
     366
     367  ! -------------------------------------------------------------------
     368  ! --- Test for instability.
     369  ! --- If there was no convection at last time step and parcel
     370  ! --- is stable at icb, then set iflag to 4.
     371  ! -------------------------------------------------------------------
     372
     373  include "cvparam.h"
     374
     375  ! inputs:
     376  INTEGER len, nd, icb(len)
     377  REAL cbmf(len), tv(len, nd), tvp(len, nd)
     378
     379  ! outputs:
     380  INTEGER iflag(len) ! also an input
     381
     382  ! local variables:
     383  INTEGER i
     384
     385
     386  DO i = 1, len
     387    IF ((cbmf(i)==0.0) .AND. (iflag(i)==0) .AND. (tvp(i, &
     388      icb(i))<=(tv(i,icb(i))-dtmax))) iflag(i) = 4
     389  END DO
     390
     391  RETURN
     392END SUBROUTINE cv_trigger
     393
     394SUBROUTINE cv_compress(len, nloc, ncum, nd, iflag1, nk1, icb1, cbmf1, plcl1, &
     395    tnk1, qnk1, gznk1, t1, q1, qs1, u1, v1, gz1, h1, lv1, cpn1, p1, ph1, tv1, &
     396    tp1, tvp1, clw1, iflag, nk, icb, cbmf, plcl, tnk, qnk, gznk, t, q, qs, u, &
     397    v, gz, h, lv, cpn, p, ph, tv, tp, tvp, clw, dph)
     398  IMPLICIT NONE
     399
     400  include "cvparam.h"
     401
     402  ! inputs:
     403  INTEGER len, ncum, nd, nloc
     404  INTEGER iflag1(len), nk1(len), icb1(len)
     405  REAL cbmf1(len), plcl1(len), tnk1(len), qnk1(len), gznk1(len)
     406  REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
     407  REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
     408  REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd)
     409  REAL tvp1(len, nd), clw1(len, nd)
     410
     411  ! outputs:
     412  INTEGER iflag(nloc), nk(nloc), icb(nloc)
     413  REAL cbmf(nloc), plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
     414  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
     415  REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
     416  REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd)
     417  REAL tvp(nloc, nd), clw(nloc, nd)
     418  REAL dph(nloc, nd)
     419
     420  ! local variables:
     421  INTEGER i, k, nn
     422  CHARACTER (LEN=20) :: modname = 'cv_compress'
     423  CHARACTER (LEN=80) :: abort_message
     424
     425  include 'iniprint.h'
     426
     427
     428  DO k = 1, nl + 1
     429    nn = 0
     430    DO i = 1, len
     431      IF (iflag1(i)==0) THEN
     432        nn = nn + 1
     433        t(nn, k) = t1(i, k)
     434        q(nn, k) = q1(i, k)
     435        qs(nn, k) = qs1(i, k)
     436        u(nn, k) = u1(i, k)
     437        v(nn, k) = v1(i, k)
     438        gz(nn, k) = gz1(i, k)
     439        h(nn, k) = h1(i, k)
     440        lv(nn, k) = lv1(i, k)
     441        cpn(nn, k) = cpn1(i, k)
     442        p(nn, k) = p1(i, k)
     443        ph(nn, k) = ph1(i, k)
     444        tv(nn, k) = tv1(i, k)
     445        tp(nn, k) = tp1(i, k)
     446        tvp(nn, k) = tvp1(i, k)
     447        clw(nn, k) = clw1(i, k)
     448      END IF
     449    END DO
     450  END DO
     451
     452  IF (nn/=ncum) THEN
     453    WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum
     454    abort_message = ''
     455    CALL abort_gcm(modname, abort_message, 1)
     456  END IF
     457
     458  nn = 0
     459  DO i = 1, len
     460    IF (iflag1(i)==0) THEN
     461      nn = nn + 1
     462      cbmf(nn) = cbmf1(i)
     463      plcl(nn) = plcl1(i)
     464      tnk(nn) = tnk1(i)
     465      qnk(nn) = qnk1(i)
     466      gznk(nn) = gznk1(i)
     467      nk(nn) = nk1(i)
     468      icb(nn) = icb1(i)
     469      iflag(nn) = iflag1(i)
     470    END IF
     471  END DO
     472
     473  DO k = 1, nl
     474    DO i = 1, ncum
     475      dph(i, k) = ph(i, k) - ph(i, k+1)
     476    END DO
     477  END DO
     478
     479  RETURN
     480END SUBROUTINE cv_compress
     481
     482SUBROUTINE cv_undilute2(nloc, ncum, nd, icb, nk, tnk, qnk, gznk, t, q, qs, &
     483    gz, p, dph, h, tv, lv, inb, inb1, tp, tvp, clw, hp, ep, sigp, frac)
     484  IMPLICIT NONE
     485
     486  ! ---------------------------------------------------------------------
     487  ! Purpose:
     488  ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
     489  ! &
     490  ! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
     491  ! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
     492  ! &
     493  ! FIND THE LEVEL OF NEUTRAL BUOYANCY
     494  ! ---------------------------------------------------------------------
     495
     496  include "cvthermo.h"
     497  include "cvparam.h"
     498
     499  ! inputs:
     500  INTEGER ncum, nd, nloc
     501  INTEGER icb(nloc), nk(nloc)
     502  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd)
     503  REAL p(nloc, nd), dph(nloc, nd)
     504  REAL tnk(nloc), qnk(nloc), gznk(nloc)
     505  REAL lv(nloc, nd), tv(nloc, nd), h(nloc, nd)
     506
     507  ! outputs:
     508  INTEGER inb(nloc), inb1(nloc)
     509  REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
     510  REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd)
     511  REAL frac(nloc)
     512
     513  ! local variables:
     514  INTEGER i, k
     515  REAL tg, qg, ahg, alv, s, tc, es, denom, rg, tca, elacrit
     516  REAL by, defrac
     517  REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
     518  LOGICAL lcape(nloc)
     519
     520  ! =====================================================================
     521  ! --- SOME INITIALIZATIONS
     522  ! =====================================================================
     523
     524  DO k = 1, nl
     525    DO i = 1, ncum
     526      ep(i, k) = 0.0
     527      sigp(i, k) = sigs
     528    END DO
     529  END DO
     530
     531  ! =====================================================================
     532  ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
     533  ! =====================================================================
     534
     535  ! ---       The procedure is to solve the equation.
     536  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
     537
     538  ! ***  Calculate certain parcel quantities, including static energy   ***
     539
     540
     541  DO i = 1, ncum
     542    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
     543      t0)) + gznk(i)
     544  END DO
     545
     546
     547  ! ***  Find lifted parcel quantities above cloud base    ***
     548
     549
     550  DO k = minorig + 1, nl
     551    DO i = 1, ncum
     552      IF (k>=(icb(i)+1)) THEN
     553        tg = t(i, k)
     554        qg = qs(i, k)
     555        alv = lv0 - clmcpv*(t(i,k)-t0)
     556
     557        ! First iteration.
     558
     559        s = cpd + alv*alv*qg/(rrv*t(i,k)*t(i,k))
     560        s = 1./s
     561        ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k)
     562        tg = tg + s*(ah0(i)-ahg)
     563        tg = max(tg, 35.0)
     564        tc = tg - t0
     565        denom = 243.5 + tc
     566        IF (tc>=0.0) THEN
     567          es = 6.112*exp(17.67*tc/denom)
     568        ELSE
     569          es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
     570        END IF
     571        qg = eps*es/(p(i,k)-es*(1.-eps))
     572
     573        ! Second iteration.
     574
     575        s = cpd + alv*alv*qg/(rrv*t(i,k)*t(i,k))
     576        s = 1./s
     577        ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k)
     578        tg = tg + s*(ah0(i)-ahg)
     579        tg = max(tg, 35.0)
     580        tc = tg - t0
     581        denom = 243.5 + tc
     582        IF (tc>=0.0) THEN
     583          es = 6.112*exp(17.67*tc/denom)
     584        ELSE
     585          es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
     586        END IF
     587        qg = eps*es/(p(i,k)-es*(1.-eps))
     588
     589        alv = lv0 - clmcpv*(t(i,k)-t0)
     590        ! print*,'cpd dans convect2 ',cpd
     591        ! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
     592        ! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
     593        tp(i, k) = (ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
     594        ! if (.not.cpd.gt.1000.) then
     595        ! print*,'CPD=',cpd
     596        ! stop
     597        ! endif
     598        clw(i, k) = qnk(i) - qg
     599        clw(i, k) = max(0.0, clw(i,k))
     600        rg = qg/(1.-qnk(i))
     601        tvp(i, k) = tp(i, k)*(1.+rg*epsi)
     602      END IF
     603    END DO
     604  END DO
     605
     606  ! =====================================================================
     607  ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
     608  ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
     609  ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
     610  ! =====================================================================
     611
     612  DO k = minorig + 1, nl
     613    DO i = 1, ncum
     614      IF (k>=(nk(i)+1)) THEN
     615        tca = tp(i, k) - t0
     616        IF (tca>=0.0) THEN
     617          elacrit = elcrit
     618        ELSE
     619          elacrit = elcrit*(1.0-tca/tlcrit)
     620        END IF
     621        elacrit = max(elacrit, 0.0)
     622        ep(i, k) = 1.0 - elacrit/max(clw(i,k), 1.0E-8)
     623        ep(i, k) = max(ep(i,k), 0.0)
     624        ep(i, k) = min(ep(i,k), 1.0)
     625        sigp(i, k) = sigs
     626      END IF
     627    END DO
     628  END DO
     629
     630  ! =====================================================================
     631  ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
     632  ! --- VIRTUAL TEMPERATURE
     633  ! =====================================================================
     634
     635  DO k = minorig + 1, nl
     636    DO i = 1, ncum
     637      IF (k>=(icb(i)+1)) THEN
     638        tvp(i, k) = tvp(i, k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
     639        ! print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
     640        ! print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
     641      END IF
     642    END DO
     643  END DO
     644  DO i = 1, ncum
     645    tvp(i, nlp) = tvp(i, nl) - (gz(i,nlp)-gz(i,nl))/cpd
     646  END DO
     647
     648  ! =====================================================================
     649  ! --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
     650  ! --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
     651  ! --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
     652  ! =====================================================================
     653
     654  DO i = 1, ncum
     655    cape(i) = 0.0
     656    capem(i) = 0.0
     657    inb(i) = icb(i) + 1
     658    inb1(i) = inb(i)
     659  END DO
     660
     661  ! Originial Code
     662
     663  ! do 530 k=minorig+1,nl-1
     664  ! do 520 i=1,ncum
     665  ! if(k.ge.(icb(i)+1))then
     666  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
     667  ! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
     668  ! cape(i)=cape(i)+by
     669  ! if(by.ge.0.0)inb1(i)=k+1
     670  ! if(cape(i).gt.0.0)then
     671  ! inb(i)=k+1
     672  ! capem(i)=cape(i)
     673  ! endif
     674  ! endif
     675  ! 520    continue
     676  ! 530  continue
     677  ! do 540 i=1,ncum
     678  ! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
     679  ! cape(i)=capem(i)+byp
     680  ! defrac=capem(i)-cape(i)
     681  ! defrac=max(defrac,0.001)
     682  ! frac(i)=-cape(i)/defrac
     683  ! frac(i)=min(frac(i),1.0)
     684  ! frac(i)=max(frac(i),0.0)
     685  ! 540   continue
     686
     687  ! K Emanuel fix
     688
     689  ! call zilch(byp,ncum)
     690  ! do 530 k=minorig+1,nl-1
     691  ! do 520 i=1,ncum
     692  ! if(k.ge.(icb(i)+1))then
     693  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
     694  ! cape(i)=cape(i)+by
     695  ! if(by.ge.0.0)inb1(i)=k+1
     696  ! if(cape(i).gt.0.0)then
     697  ! inb(i)=k+1
     698  ! capem(i)=cape(i)
     699  ! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
     700  ! endif
     701  ! endif
     702  ! 520    continue
     703  ! 530  continue
     704  ! do 540 i=1,ncum
     705  ! inb(i)=max(inb(i),inb1(i))
     706  ! cape(i)=capem(i)+byp(i)
     707  ! defrac=capem(i)-cape(i)
     708  ! defrac=max(defrac,0.001)
     709  ! frac(i)=-cape(i)/defrac
     710  ! frac(i)=min(frac(i),1.0)
     711  ! frac(i)=max(frac(i),0.0)
     712  ! 540   continue
     713
     714  ! J Teixeira fix
     715
     716  CALL zilch(byp, ncum)
     717  DO i = 1, ncum
     718    lcape(i) = .TRUE.
     719  END DO
     720  DO k = minorig + 1, nl - 1
     721    DO i = 1, ncum
     722      IF (cape(i)<0.0) lcape(i) = .FALSE.
     723      IF ((k>=(icb(i)+1)) .AND. lcape(i)) THEN
     724        by = (tvp(i,k)-tv(i,k))*dph(i, k)/p(i, k)
     725        byp(i) = (tvp(i,k+1)-tv(i,k+1))*dph(i, k+1)/p(i, k+1)
     726        cape(i) = cape(i) + by
     727        IF (by>=0.0) inb1(i) = k + 1
     728        IF (cape(i)>0.0) THEN
     729          inb(i) = k + 1
     730          capem(i) = cape(i)
     731        END IF
     732      END IF
     733    END DO
     734  END DO
     735  DO i = 1, ncum
     736    cape(i) = capem(i) + byp(i)
     737    defrac = capem(i) - cape(i)
     738    defrac = max(defrac, 0.001)
     739    frac(i) = -cape(i)/defrac
     740    frac(i) = min(frac(i), 1.0)
     741    frac(i) = max(frac(i), 0.0)
     742  END DO
     743
     744  ! =====================================================================
     745  ! ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
     746  ! =====================================================================
     747
     748  ! initialization:
     749  DO i = 1, ncum*nlp
     750    hp(i, 1) = h(i, 1)
     751  END DO
     752
     753  DO k = minorig + 1, nl
     754    DO i = 1, ncum
     755      IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
     756        hp(i, k) = h(i, nk(i)) + (lv(i,k)+(cpd-cpv)*t(i,k))*ep(i, k)*clw(i, k &
     757          )
     758      END IF
     759    END DO
     760  END DO
     761
     762  RETURN
     763END SUBROUTINE cv_undilute2
     764
     765SUBROUTINE cv_closure(nloc, ncum, nd, nk, icb, tv, tvp, p, ph, dph, plcl, &
     766    cpn, iflag, cbmf)
     767  IMPLICIT NONE
     768
     769  ! inputs:
     770  INTEGER ncum, nd, nloc
     771  INTEGER nk(nloc), icb(nloc)
     772  REAL tv(nloc, nd), tvp(nloc, nd), p(nloc, nd), dph(nloc, nd)
     773  REAL ph(nloc, nd+1) ! caution nd instead ndp1 to be consistent...
     774  REAL plcl(nloc), cpn(nloc, nd)
     775
     776  ! outputs:
     777  INTEGER iflag(nloc)
     778  REAL cbmf(nloc) ! also an input
     779
     780  ! local variables:
     781  INTEGER i, k, icbmax
     782  REAL dtpbl(nloc), dtmin(nloc), tvpplcl(nloc), tvaplcl(nloc)
     783  REAL work(nloc)
     784
     785  include "cvthermo.h"
     786  include "cvparam.h"
     787
     788  ! -------------------------------------------------------------------
     789  ! Compute icbmax.
     790  ! -------------------------------------------------------------------
     791
     792  icbmax = 2
     793  DO i = 1, ncum
     794    icbmax = max(icbmax, icb(i))
     795  END DO
     796
     797  ! =====================================================================
     798  ! ---  CALCULATE CLOUD BASE MASS FLUX
     799  ! =====================================================================
     800
     801  ! tvpplcl = parcel temperature lifted adiabatically from level
     802  ! icb-1 to the LCL.
     803  ! tvaplcl = virtual temperature at the LCL.
     804
     805  DO i = 1, ncum
     806    dtpbl(i) = 0.0
     807    tvpplcl(i) = tvp(i, icb(i)-1) - rrd*tvp(i, icb(i)-1)*(p(i,icb(i)-1)-plcl( &
     808      i))/(cpn(i,icb(i)-1)*p(i,icb(i)-1))
     809    tvaplcl(i) = tv(i, icb(i)) + (tvp(i,icb(i))-tvp(i,icb(i)+1))*(plcl(i)-p(i &
     810      ,icb(i)))/(p(i,icb(i))-p(i,icb(i)+1))
     811  END DO
     812
     813  ! -------------------------------------------------------------------
     814  ! --- Interpolate difference between lifted parcel and
     815  ! --- environmental temperatures to lifted condensation level
     816  ! -------------------------------------------------------------------
     817
     818  ! dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
     819
     820  DO k = minorig, icbmax
     821    DO i = 1, ncum
     822      IF ((k>=nk(i)) .AND. (k<=(icb(i)-1))) THEN
     823        dtpbl(i) = dtpbl(i) + (tvp(i,k)-tv(i,k))*dph(i, k)
     824      END IF
     825    END DO
     826  END DO
     827  DO i = 1, ncum
     828    dtpbl(i) = dtpbl(i)/(ph(i,nk(i))-ph(i,icb(i)))
     829    dtmin(i) = tvpplcl(i) - tvaplcl(i) + dtmax + dtpbl(i)
     830  END DO
     831
     832  ! -------------------------------------------------------------------
     833  ! --- Adjust cloud base mass flux
     834  ! -------------------------------------------------------------------
     835
     836  DO i = 1, ncum
     837    work(i) = cbmf(i)
     838    cbmf(i) = max(0.0, (1.0-damp)*cbmf(i)+0.1*alpha*dtmin(i))
     839    IF ((work(i)==0.0) .AND. (cbmf(i)==0.0)) THEN
     840      iflag(i) = 3
     841    END IF
     842  END DO
     843
     844  RETURN
     845END SUBROUTINE cv_closure
     846
     847SUBROUTINE cv_mixing(nloc, ncum, nd, icb, nk, inb, inb1, ph, t, q, qs, u, v, &
     848    h, lv, qnk, hp, tv, tvp, ep, clw, cbmf, m, ment, qent, uent, vent, nent, &
     849    sij, elij)
     850  IMPLICIT NONE
     851
     852  include "cvthermo.h"
     853  include "cvparam.h"
     854
     855  ! inputs:
     856  INTEGER ncum, nd, nloc
     857  INTEGER icb(nloc), inb(nloc), inb1(nloc), nk(nloc)
     858  REAL cbmf(nloc), qnk(nloc)
     859  REAL ph(nloc, nd+1)
     860  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), lv(nloc, nd)
     861  REAL u(nloc, nd), v(nloc, nd), h(nloc, nd), hp(nloc, nd)
     862  REAL tv(nloc, nd), tvp(nloc, nd), ep(nloc, nd), clw(nloc, nd)
     863
     864  ! outputs:
     865  INTEGER nent(nloc, nd)
     866  REAL m(nloc, nd), ment(nloc, nd, nd), qent(nloc, nd, nd)
     867  REAL uent(nloc, nd, nd), vent(nloc, nd, nd)
     868  REAL sij(nloc, nd, nd), elij(nloc, nd, nd)
     869
     870  ! local variables:
     871  INTEGER i, j, k, ij
     872  INTEGER num1, num2
     873  REAL dbo, qti, bf2, anum, denom, dei, altem, cwat, stemp
     874  REAL alt, qp1, smid, sjmin, sjmax, delp, delm
     875  REAL work(nloc), asij(nloc), smin(nloc), scrit(nloc)
     876  REAL bsum(nloc, nd)
     877  LOGICAL lwork(nloc)
     878
     879  ! =====================================================================
     880  ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
     881  ! =====================================================================
     882
     883  DO i = 1, ncum*nlp
     884    nent(i, 1) = 0
     885    m(i, 1) = 0.0
     886  END DO
     887
     888  DO k = 1, nlp
     889    DO j = 1, nlp
     890      DO i = 1, ncum
     891        qent(i, k, j) = q(i, j)
     892        uent(i, k, j) = u(i, j)
     893        vent(i, k, j) = v(i, j)
     894        elij(i, k, j) = 0.0
     895        ment(i, k, j) = 0.0
     896        sij(i, k, j) = 0.0
     897      END DO
     898    END DO
     899  END DO
     900
     901  ! -------------------------------------------------------------------
     902  ! --- Calculate rates of mixing,  m(i)
     903  ! -------------------------------------------------------------------
     904
     905  CALL zilch(work, ncum)
     906
     907  DO j = minorig + 1, nl
     908    DO i = 1, ncum
     909      IF ((j>=(icb(i)+1)) .AND. (j<=inb(i))) THEN
     910        k = min(j, inb1(i))
     911        dbo = abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1)) + &
     912          entp*0.04*(ph(i,k)-ph(i,k+1))
     913        work(i) = work(i) + dbo
     914        m(i, j) = cbmf(i)*dbo
     915      END IF
     916    END DO
     917  END DO
     918  DO k = minorig + 1, nl
     919    DO i = 1, ncum
     920      IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
     921        m(i, k) = m(i, k)/work(i)
     922      END IF
     923    END DO
     924  END DO
     925
     926
     927  ! =====================================================================
     928  ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
     929  ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
     930  ! --- FRACTION (sij)
     931  ! =====================================================================
     932
     933
     934  DO i = minorig + 1, nl
     935    DO j = minorig + 1, nl
     936      DO ij = 1, ncum
     937        IF ((i>=(icb(ij)+1)) .AND. (j>=icb(ij)) .AND. (i<=inb(ij)) .AND. (j<= &
     938            inb(ij))) THEN
     939          qti = qnk(ij) - ep(ij, i)*clw(ij, i)
     940          bf2 = 1. + lv(ij, j)*lv(ij, j)*qs(ij, j)/(rrv*t(ij,j)*t(ij,j)*cpd)
     941          anum = h(ij, j) - hp(ij, i) + (cpv-cpd)*t(ij, j)*(qti-q(ij,j))
     942          denom = h(ij, i) - hp(ij, i) + (cpd-cpv)*(q(ij,i)-qti)*t(ij, j)
     943          dei = denom
     944          IF (abs(dei)<0.01) dei = 0.01
     945          sij(ij, i, j) = anum/dei
     946          sij(ij, i, i) = 1.0
     947          altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j)
     948          altem = altem/bf2
     949          cwat = clw(ij, j)*(1.-ep(ij,j))
     950          stemp = sij(ij, i, j)
     951          IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
     952            anum = anum - lv(ij, j)*(qti-qs(ij,j)-cwat*bf2)
     953            denom = denom + lv(ij, j)*(q(ij,i)-qti)
     954            IF (abs(denom)<0.01) denom = 0.01
     955            sij(ij, i, j) = anum/denom
     956            altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j)
     957            altem = altem - (bf2-1.)*cwat
     958          END IF
     959          IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN
     960            qent(ij, i, j) = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti
     961            uent(ij, i, j) = sij(ij, i, j)*u(ij, i) + &
     962              (1.-sij(ij,i,j))*u(ij, nk(ij))
     963            vent(ij, i, j) = sij(ij, i, j)*v(ij, i) + &
     964              (1.-sij(ij,i,j))*v(ij, nk(ij))
     965            elij(ij, i, j) = altem
     966            elij(ij, i, j) = max(0.0, elij(ij,i,j))
     967            ment(ij, i, j) = m(ij, i)/(1.-sij(ij,i,j))
     968            nent(ij, i) = nent(ij, i) + 1
     969          END IF
     970          sij(ij, i, j) = max(0.0, sij(ij,i,j))
     971          sij(ij, i, j) = min(1.0, sij(ij,i,j))
     972        END IF
     973      END DO
     974    END DO
     975
     976    ! ***   If no air can entrain at level i assume that updraft detrains
     977    ! ***
     978    ! ***   at that level and calculate detrained air flux and properties
     979    ! ***
     980
     981    DO ij = 1, ncum
     982      IF ((i>=(icb(ij)+1)) .AND. (i<=inb(ij)) .AND. (nent(ij,i)==0)) THEN
     983        ment(ij, i, i) = m(ij, i)
     984        qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
     985        uent(ij, i, i) = u(ij, nk(ij))
     986        vent(ij, i, i) = v(ij, nk(ij))
     987        elij(ij, i, i) = clw(ij, i)
     988        sij(ij, i, i) = 1.0
     989      END IF
     990    END DO
     991  END DO
     992
     993  DO i = 1, ncum
     994    sij(i, inb(i), inb(i)) = 1.0
     995  END DO
     996
     997  ! =====================================================================
     998  ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
     999  ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
     1000  ! =====================================================================
     1001
     1002  CALL zilch(bsum, ncum*nlp)
     1003  DO ij = 1, ncum
     1004    lwork(ij) = .FALSE.
     1005  END DO
     1006  DO i = minorig + 1, nl
     1007
     1008    num1 = 0
     1009    DO ij = 1, ncum
     1010      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) num1 = num1 + 1
     1011    END DO
     1012    IF (num1<=0) GO TO 789
     1013
     1014    DO ij = 1, ncum
     1015      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) THEN
     1016        lwork(ij) = (nent(ij,i)/=0)
     1017        qp1 = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
     1018        anum = h(ij, i) - hp(ij, i) - lv(ij, i)*(qp1-qs(ij,i))
     1019        denom = h(ij, i) - hp(ij, i) + lv(ij, i)*(q(ij,i)-qp1)
     1020        IF (abs(denom)<0.01) denom = 0.01
     1021        scrit(ij) = anum/denom
     1022        alt = qp1 - qs(ij, i) + scrit(ij)*(q(ij,i)-qp1)
     1023        IF (scrit(ij)<0.0 .OR. alt<0.0) scrit(ij) = 1.0
     1024        asij(ij) = 0.0
     1025        smin(ij) = 1.0
     1026      END IF
     1027    END DO
     1028    DO j = minorig, nl
     1029
     1030      num2 = 0
     1031      DO ij = 1, ncum
     1032        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
     1033          ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) num2 = num2 + 1
     1034      END DO
     1035      IF (num2<=0) GO TO 783
     1036
     1037      DO ij = 1, ncum
     1038        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
     1039            ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
     1040          IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN
     1041            IF (j>i) THEN
     1042              smid = min(sij(ij,i,j), scrit(ij))
     1043              sjmax = smid
     1044              sjmin = smid
     1045              IF (smid<smin(ij) .AND. sij(ij,i,j+1)<smid) THEN
     1046                smin(ij) = smid
     1047                sjmax = min(sij(ij,i,j+1), sij(ij,i,j), scrit(ij))
     1048                sjmin = max(sij(ij,i,j-1), sij(ij,i,j))
     1049                sjmin = min(sjmin, scrit(ij))
     1050              END IF
     1051            ELSE
     1052              sjmax = max(sij(ij,i,j+1), scrit(ij))
     1053              smid = max(sij(ij,i,j), scrit(ij))
     1054              sjmin = 0.0
     1055              IF (j>1) sjmin = sij(ij, i, j-1)
     1056              sjmin = max(sjmin, scrit(ij))
     1057            END IF
     1058            delp = abs(sjmax-smid)
     1059            delm = abs(sjmin-smid)
     1060            asij(ij) = asij(ij) + (delp+delm)*(ph(ij,j)-ph(ij,j+1))
     1061            ment(ij, i, j) = ment(ij, i, j)*(delp+delm)*(ph(ij,j)-ph(ij,j+1))
     1062          END IF
     1063        END IF
     1064      END DO
     1065783 END DO
     1066    DO ij = 1, ncum
     1067      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. lwork(ij)) THEN
     1068        asij(ij) = max(1.0E-21, asij(ij))
     1069        asij(ij) = 1.0/asij(ij)
     1070        bsum(ij, i) = 0.0
     1071      END IF
     1072    END DO
     1073    DO j = minorig, nl + 1
     1074      DO ij = 1, ncum
     1075        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
     1076            ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
     1077          ment(ij, i, j) = ment(ij, i, j)*asij(ij)
     1078          bsum(ij, i) = bsum(ij, i) + ment(ij, i, j)
     1079        END IF
     1080      END DO
     1081    END DO
     1082    DO ij = 1, ncum
     1083      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (bsum(ij, &
     1084          i)<1.0E-18) .AND. lwork(ij)) THEN
     1085        nent(ij, i) = 0
     1086        ment(ij, i, i) = m(ij, i)
     1087        qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
     1088        uent(ij, i, i) = u(ij, nk(ij))
     1089        vent(ij, i, i) = v(ij, nk(ij))
     1090        elij(ij, i, i) = clw(ij, i)
     1091        sij(ij, i, i) = 1.0
     1092      END IF
     1093    END DO
     1094789 END DO
     1095
     1096  RETURN
     1097END SUBROUTINE cv_mixing
     1098
     1099SUBROUTINE cv_unsat(nloc, ncum, nd, inb, t, q, qs, gz, u, v, p, ph, h, lv, &
     1100    ep, sigp, clw, m, ment, elij, iflag, mp, qp, up, vp, wt, water, evap)
     1101  IMPLICIT NONE
     1102
     1103
     1104  include "cvthermo.h"
     1105  include "cvparam.h"
     1106
     1107  ! inputs:
     1108  INTEGER ncum, nd, nloc
     1109  INTEGER inb(nloc)
     1110  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd)
     1111  REAL gz(nloc, nd), u(nloc, nd), v(nloc, nd)
     1112  REAL p(nloc, nd), ph(nloc, nd+1), h(nloc, nd)
     1113  REAL lv(nloc, nd), ep(nloc, nd), sigp(nloc, nd), clw(nloc, nd)
     1114  REAL m(nloc, nd), ment(nloc, nd, nd), elij(nloc, nd, nd)
     1115
     1116  ! outputs:
     1117  INTEGER iflag(nloc) ! also an input
     1118  REAL mp(nloc, nd), qp(nloc, nd), up(nloc, nd), vp(nloc, nd)
     1119  REAL water(nloc, nd), evap(nloc, nd), wt(nloc, nd)
     1120
     1121  ! local variables:
     1122  INTEGER i, j, k, ij, num1
     1123  INTEGER jtt(nloc)
     1124  REAL awat, coeff, qsm, afac, sigt, b6, c6, revap
     1125  REAL dhdp, fac, qstm, rat
     1126  REAL wdtrain(nloc)
     1127  LOGICAL lwork(nloc)
     1128
     1129  ! =====================================================================
     1130  ! --- PRECIPITATING DOWNDRAFT CALCULATION
     1131  ! =====================================================================
     1132
     1133  ! Initializations:
     1134
     1135  DO i = 1, ncum
     1136    DO k = 1, nl + 1
     1137      wt(i, k) = omtsnow
     1138      mp(i, k) = 0.0
     1139      evap(i, k) = 0.0
     1140      water(i, k) = 0.0
     1141    END DO
     1142  END DO
     1143
     1144  DO i = 1, ncum
     1145    qp(i, 1) = q(i, 1)
     1146    up(i, 1) = u(i, 1)
     1147    vp(i, 1) = v(i, 1)
     1148  END DO
     1149
     1150  DO k = 2, nl + 1
     1151    DO i = 1, ncum
     1152      qp(i, k) = q(i, k-1)
     1153      up(i, k) = u(i, k-1)
     1154      vp(i, k) = v(i, k-1)
     1155    END DO
     1156  END DO
     1157
     1158
     1159  ! ***  Check whether ep(inb)=0, if so, skip precipitating    ***
     1160  ! ***             downdraft calculation                      ***
     1161
     1162
     1163  ! ***  Integrate liquid water equation to find condensed water   ***
     1164  ! ***                and condensed water flux                    ***
     1165
     1166
     1167  DO i = 1, ncum
     1168    jtt(i) = 2
     1169    IF (ep(i,inb(i))<=0.0001) iflag(i) = 2
     1170    IF (iflag(i)==0) THEN
     1171      lwork(i) = .TRUE.
     1172    ELSE
     1173      lwork(i) = .FALSE.
     1174    END IF
     1175  END DO
     1176
     1177  ! ***                    Begin downdraft loop                    ***
     1178
     1179
     1180  CALL zilch(wdtrain, ncum)
     1181  DO i = nl + 1, 1, -1
     1182
     1183    num1 = 0
     1184    DO ij = 1, ncum
     1185      IF ((i<=inb(ij)) .AND. lwork(ij)) num1 = num1 + 1
     1186    END DO
     1187    IF (num1<=0) GO TO 899
     1188
     1189
     1190    ! ***        Calculate detrained precipitation             ***
     1191
     1192    DO ij = 1, ncum
     1193      IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
     1194        wdtrain(ij) = g*ep(ij, i)*m(ij, i)*clw(ij, i)
     1195      END IF
     1196    END DO
     1197
     1198    IF (i>1) THEN
     1199      DO j = 1, i - 1
     1200        DO ij = 1, ncum
     1201          IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
     1202            awat = elij(ij, j, i) - (1.-ep(ij,i))*clw(ij, i)
     1203            awat = max(0.0, awat)
     1204            wdtrain(ij) = wdtrain(ij) + g*awat*ment(ij, j, i)
     1205          END IF
     1206        END DO
     1207      END DO
     1208    END IF
     1209
     1210    ! ***    Find rain water and evaporation using provisional   ***
     1211    ! ***              estimates of qp(i)and qp(i-1)             ***
     1212
     1213
     1214    ! ***  Value of terminal velocity and coeffecient of evaporation for snow
     1215    ! ***
     1216
     1217    DO ij = 1, ncum
     1218      IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
     1219        coeff = coeffs
     1220        wt(ij, i) = omtsnow
     1221
     1222        ! ***  Value of terminal velocity and coeffecient of evaporation for
     1223        ! rain   ***
     1224
     1225        IF (t(ij,i)>273.0) THEN
     1226          coeff = coeffr
     1227          wt(ij, i) = omtrain
     1228        END IF
     1229        qsm = 0.5*(q(ij,i)+qp(ij,i+1))
     1230        afac = coeff*ph(ij, i)*(qs(ij,i)-qsm)/(1.0E4+2.0E3*ph(ij,i)*qs(ij,i))
     1231        afac = max(afac, 0.0)
     1232        sigt = sigp(ij, i)
     1233        sigt = max(0.0, sigt)
     1234        sigt = min(1.0, sigt)
     1235        b6 = 100.*(ph(ij,i)-ph(ij,i+1))*sigt*afac/wt(ij, i)
     1236        c6 = (water(ij,i+1)*wt(ij,i+1)+wdtrain(ij)/sigd)/wt(ij, i)
     1237        revap = 0.5*(-b6+sqrt(b6*b6+4.*c6))
     1238        evap(ij, i) = sigt*afac*revap
     1239        water(ij, i) = revap*revap
     1240
     1241        ! ***  Calculate precipitating downdraft mass flux under     ***
     1242        ! ***              hydrostatic approximation                 ***
     1243
     1244        IF (i>1) THEN
     1245          dhdp = (h(ij,i)-h(ij,i-1))/(p(ij,i-1)-p(ij,i))
     1246          dhdp = max(dhdp, 10.0)
     1247          mp(ij, i) = 100.*ginv*lv(ij, i)*sigd*evap(ij, i)/dhdp
     1248          mp(ij, i) = max(mp(ij,i), 0.0)
     1249
     1250          ! ***   Add small amount of inertia to downdraft              ***
     1251
     1252          fac = 20.0/(ph(ij,i-1)-ph(ij,i))
     1253          mp(ij, i) = (fac*mp(ij,i+1)+mp(ij,i))/(1.+fac)
     1254
     1255          ! ***      Force mp to decrease linearly to zero
     1256          ! ***
     1257          ! ***      between about 950 mb and the surface
     1258          ! ***
     1259
     1260          IF (p(ij,i)>(0.949*p(ij,1))) THEN
     1261            jtt(ij) = max(jtt(ij), i)
     1262            mp(ij, i) = mp(ij, jtt(ij))*(p(ij,1)-p(ij,i))/ &
     1263              (p(ij,1)-p(ij,jtt(ij)))
     1264          END IF
     1265        END IF
     1266
     1267        ! ***       Find mixing ratio of precipitating downdraft     ***
     1268
     1269        IF (i/=inb(ij)) THEN
     1270          IF (i==1) THEN
     1271            qstm = qs(ij, 1)
     1272          ELSE
     1273            qstm = qs(ij, i-1)
     1274          END IF
     1275          IF (mp(ij,i)>mp(ij,i+1)) THEN
     1276            rat = mp(ij, i+1)/mp(ij, i)
     1277            qp(ij, i) = qp(ij, i+1)*rat + q(ij, i)*(1.0-rat) + &
     1278              100.*ginv*sigd*(ph(ij,i)-ph(ij,i+1))*(evap(ij,i)/mp(ij,i))
     1279            up(ij, i) = up(ij, i+1)*rat + u(ij, i)*(1.-rat)
     1280            vp(ij, i) = vp(ij, i+1)*rat + v(ij, i)*(1.-rat)
     1281          ELSE
     1282            IF (mp(ij,i+1)>0.0) THEN
     1283              qp(ij, i) = (gz(ij,i+1)-gz(ij,i)+qp(ij,i+1)*(lv(ij,i+1)+t(ij, &
     1284                i+1)*(cl-cpd))+cpd*(t(ij,i+1)-t(ij, &
     1285                i)))/(lv(ij,i)+t(ij,i)*(cl-cpd))
     1286              up(ij, i) = up(ij, i+1)
     1287              vp(ij, i) = vp(ij, i+1)
     1288            END IF
     1289          END IF
     1290          qp(ij, i) = min(qp(ij,i), qstm)
     1291          qp(ij, i) = max(qp(ij,i), 0.0)
     1292        END IF
     1293      END IF
     1294    END DO
     1295899 END DO
     1296
     1297  RETURN
     1298END SUBROUTINE cv_unsat
     1299
     1300SUBROUTINE cv_yield(nloc, ncum, nd, nk, icb, inb, delt, t, q, u, v, gz, p, &
     1301    ph, h, hp, lv, cpn, ep, clw, frac, m, mp, qp, up, vp, wt, water, evap, &
     1302    ment, qent, uent, vent, nent, elij, tv, tvp, iflag, wd, qprime, tprime, &
     1303    precip, cbmf, ft, fq, fu, fv, ma, qcondc)
     1304  IMPLICIT NONE
     1305
     1306  include "cvthermo.h"
     1307  include "cvparam.h"
     1308
     1309  ! inputs
     1310  INTEGER ncum, nd, nloc
     1311  INTEGER nk(nloc), icb(nloc), inb(nloc)
     1312  INTEGER nent(nloc, nd)
     1313  REAL delt
     1314  REAL t(nloc, nd), q(nloc, nd), u(nloc, nd), v(nloc, nd)
     1315  REAL gz(nloc, nd)
     1316  REAL p(nloc, nd), ph(nloc, nd+1), h(nloc, nd)
     1317  REAL hp(nloc, nd), lv(nloc, nd)
     1318  REAL cpn(nloc, nd), ep(nloc, nd), clw(nloc, nd), frac(nloc)
     1319  REAL m(nloc, nd), mp(nloc, nd), qp(nloc, nd)
     1320  REAL up(nloc, nd), vp(nloc, nd)
     1321  REAL wt(nloc, nd), water(nloc, nd), evap(nloc, nd)
     1322  REAL ment(nloc, nd, nd), qent(nloc, nd, nd), elij(nloc, nd, nd)
     1323  REAL uent(nloc, nd, nd), vent(nloc, nd, nd)
     1324  REAL tv(nloc, nd), tvp(nloc, nd)
     1325
     1326  ! outputs
     1327  INTEGER iflag(nloc) ! also an input
     1328  REAL cbmf(nloc) ! also an input
     1329  REAL wd(nloc), tprime(nloc), qprime(nloc)
     1330  REAL precip(nloc)
     1331  REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
     1332  REAL ma(nloc, nd)
     1333  REAL qcondc(nloc, nd)
     1334
     1335  ! local variables
     1336  INTEGER i, j, ij, k, num1
     1337  REAL dpinv, cpinv, awat, fqold, ftold, fuold, fvold, delti
     1338  REAL work(nloc), am(nloc), amp1(nloc), ad(nloc)
     1339  REAL ents(nloc), uav(nloc), vav(nloc), lvcp(nloc, nd)
     1340  REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
     1341  REAL siga(nloc, nd), ax(nloc, nd), mac(nloc, nd) ! cld
     1342
     1343
     1344  ! -- initializations:
     1345
     1346  delti = 1.0/delt
     1347
     1348  DO i = 1, ncum
     1349    precip(i) = 0.0
     1350    wd(i) = 0.0
     1351    tprime(i) = 0.0
     1352    qprime(i) = 0.0
     1353    DO k = 1, nl + 1
     1354      ft(i, k) = 0.0
     1355      fu(i, k) = 0.0
     1356      fv(i, k) = 0.0
     1357      fq(i, k) = 0.0
     1358      lvcp(i, k) = lv(i, k)/cpn(i, k)
     1359      qcondc(i, k) = 0.0 ! cld
     1360      qcond(i, k) = 0.0 ! cld
     1361      nqcond(i, k) = 0.0 ! cld
     1362    END DO
     1363  END DO
     1364
     1365
     1366  ! ***  Calculate surface precipitation in mm/day     ***
     1367
     1368  DO i = 1, ncum
     1369    IF (iflag(i)<=1) THEN
     1370      ! c            precip(i)=precip(i)+wt(i,1)*sigd*water(i,1)*3600.*24000.
     1371      ! c     &                /(rowl*g)
     1372      ! c            precip(i)=precip(i)*delt/86400.
     1373      precip(i) = wt(i, 1)*sigd*water(i, 1)*86400/g
     1374    END IF
     1375  END DO
     1376
     1377
     1378  ! ***  Calculate downdraft velocity scale and surface temperature and  ***
     1379  ! ***                    water vapor fluctuations                      ***
     1380
     1381  DO i = 1, ncum
     1382    wd(i) = betad*abs(mp(i,icb(i)))*0.01*rrd*t(i, icb(i))/(sigd*p(i,icb(i)))
     1383    qprime(i) = 0.5*(qp(i,1)-q(i,1))
     1384    tprime(i) = lv0*qprime(i)/cpd
     1385  END DO
     1386
     1387  ! ***  Calculate tendencies of lowest level potential temperature  ***
     1388  ! ***                      and mixing ratio                        ***
     1389
     1390  DO i = 1, ncum
     1391    work(i) = 0.01/(ph(i,1)-ph(i,2))
     1392    am(i) = 0.0
     1393  END DO
     1394  DO k = 2, nl
     1395    DO i = 1, ncum
     1396      IF ((nk(i)==1) .AND. (k<=inb(i)) .AND. (nk(i)==1)) THEN
     1397        am(i) = am(i) + m(i, k)
     1398      END IF
     1399    END DO
     1400  END DO
     1401  DO i = 1, ncum
     1402    IF ((g*work(i)*am(i))>=delti) iflag(i) = 1
     1403    ft(i, 1) = ft(i, 1) + g*work(i)*am(i)*(t(i,2)-t(i,1)+(gz(i,2)-gz(i, &
     1404      1))/cpn(i,1))
     1405    ft(i, 1) = ft(i, 1) - lvcp(i, 1)*sigd*evap(i, 1)
     1406    ft(i, 1) = ft(i, 1) + sigd*wt(i, 2)*(cl-cpd)*water(i, 2)*(t(i,2)-t(i,1))* &
     1407      work(i)/cpn(i, 1)
     1408    fq(i, 1) = fq(i, 1) + g*mp(i, 2)*(qp(i,2)-q(i,1))*work(i) + &
     1409      sigd*evap(i, 1)
     1410    fq(i, 1) = fq(i, 1) + g*am(i)*(q(i,2)-q(i,1))*work(i)
     1411    fu(i, 1) = fu(i, 1) + g*work(i)*(mp(i,2)*(up(i,2)-u(i,1))+am(i)*(u(i, &
     1412      2)-u(i,1)))
     1413    fv(i, 1) = fv(i, 1) + g*work(i)*(mp(i,2)*(vp(i,2)-v(i,1))+am(i)*(v(i, &
     1414      2)-v(i,1)))
     1415  END DO
     1416  DO j = 2, nl
     1417    DO i = 1, ncum
     1418      IF (j<=inb(i)) THEN
     1419        fq(i, 1) = fq(i, 1) + g*work(i)*ment(i, j, 1)*(qent(i,j,1)-q(i,1))
     1420        fu(i, 1) = fu(i, 1) + g*work(i)*ment(i, j, 1)*(uent(i,j,1)-u(i,1))
     1421        fv(i, 1) = fv(i, 1) + g*work(i)*ment(i, j, 1)*(vent(i,j,1)-v(i,1))
     1422      END IF
     1423    END DO
     1424  END DO
     1425
     1426  ! ***  Calculate tendencies of potential temperature and mixing ratio  ***
     1427  ! ***               at levels above the lowest level                   ***
     1428
     1429  ! ***  First find the net saturated updraft and downdraft mass fluxes  ***
     1430  ! ***                      through each level                          ***
     1431
     1432  DO i = 2, nl + 1
     1433
     1434    num1 = 0
     1435    DO ij = 1, ncum
     1436      IF (i<=inb(ij)) num1 = num1 + 1
     1437    END DO
     1438    IF (num1<=0) GO TO 1500
     1439
     1440    CALL zilch(amp1, ncum)
     1441    CALL zilch(ad, ncum)
     1442
     1443    DO k = i + 1, nl + 1
     1444      DO ij = 1, ncum
     1445        IF ((i>=nk(ij)) .AND. (i<=inb(ij)) .AND. (k<=(inb(ij)+1))) THEN
     1446          amp1(ij) = amp1(ij) + m(ij, k)
     1447        END IF
     1448      END DO
     1449    END DO
     1450
     1451    DO k = 1, i
     1452      DO j = i + 1, nl + 1
     1453        DO ij = 1, ncum
     1454          IF ((j<=(inb(ij)+1)) .AND. (i<=inb(ij))) THEN
     1455            amp1(ij) = amp1(ij) + ment(ij, k, j)
     1456          END IF
     1457        END DO
     1458      END DO
     1459    END DO
     1460    DO k = 1, i - 1
     1461      DO j = i, nl + 1
     1462        DO ij = 1, ncum
     1463          IF ((i<=inb(ij)) .AND. (j<=inb(ij))) THEN
     1464            ad(ij) = ad(ij) + ment(ij, j, k)
     1465          END IF
     1466        END DO
     1467      END DO
     1468    END DO
     1469
     1470    DO ij = 1, ncum
     1471      IF (i<=inb(ij)) THEN
     1472        dpinv = 0.01/(ph(ij,i)-ph(ij,i+1))
     1473        cpinv = 1.0/cpn(ij, i)
     1474
     1475        ft(ij, i) = ft(ij, i) + g*dpinv*(amp1(ij)*(t(ij,i+1)-t(ij, &
     1476          i)+(gz(ij,i+1)-gz(ij,i))*cpinv)-ad(ij)*(t(ij,i)-t(ij, &
     1477          i-1)+(gz(ij,i)-gz(ij,i-1))*cpinv)) - sigd*lvcp(ij, i)*evap(ij, i)
     1478        ft(ij, i) = ft(ij, i) + g*dpinv*ment(ij, i, i)*(hp(ij,i)-h(ij,i)+t(ij &
     1479          ,i)*(cpv-cpd)*(q(ij,i)-qent(ij,i,i)))*cpinv
     1480        ft(ij, i) = ft(ij, i) + sigd*wt(ij, i+1)*(cl-cpd)*water(ij, i+1)*(t( &
     1481          ij,i+1)-t(ij,i))*dpinv*cpinv
     1482        fq(ij, i) = fq(ij, i) + g*dpinv*(amp1(ij)*(q(ij,i+1)-q(ij, &
     1483          i))-ad(ij)*(q(ij,i)-q(ij,i-1)))
     1484        fu(ij, i) = fu(ij, i) + g*dpinv*(amp1(ij)*(u(ij,i+1)-u(ij, &
     1485          i))-ad(ij)*(u(ij,i)-u(ij,i-1)))
     1486        fv(ij, i) = fv(ij, i) + g*dpinv*(amp1(ij)*(v(ij,i+1)-v(ij, &
     1487          i))-ad(ij)*(v(ij,i)-v(ij,i-1)))
     1488      END IF
     1489    END DO
     1490    DO k = 1, i - 1
     1491      DO ij = 1, ncum
     1492        IF (i<=inb(ij)) THEN
     1493          awat = elij(ij, k, i) - (1.-ep(ij,i))*clw(ij, i)
     1494          awat = max(awat, 0.0)
     1495          fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-awat-q &
     1496            (ij,i))
     1497          fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i &
     1498            ))
     1499          fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i &
     1500            ))
     1501          ! (saturated updrafts resulting from mixing)               ! cld
     1502          qcond(ij, i) = qcond(ij, i) + (elij(ij,k,i)-awat) ! cld
     1503          nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
     1504        END IF
     1505      END DO
     1506    END DO
     1507    DO k = i, nl + 1
     1508      DO ij = 1, ncum
     1509        IF ((i<=inb(ij)) .AND. (k<=inb(ij))) THEN
     1510          fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-q(ij,i &
     1511            ))
     1512          fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i &
     1513            ))
     1514          fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i &
     1515            ))
     1516        END IF
     1517      END DO
     1518    END DO
     1519    DO ij = 1, ncum
     1520      IF (i<=inb(ij)) THEN
     1521        fq(ij, i) = fq(ij, i) + sigd*evap(ij, i) + g*(mp(ij,i+1)*(qp(ij, &
     1522          i+1)-q(ij,i))-mp(ij,i)*(qp(ij,i)-q(ij,i-1)))*dpinv
     1523        fu(ij, i) = fu(ij, i) + g*(mp(ij,i+1)*(up(ij,i+1)-u(ij, &
     1524          i))-mp(ij,i)*(up(ij,i)-u(ij,i-1)))*dpinv
     1525        fv(ij, i) = fv(ij, i) + g*(mp(ij,i+1)*(vp(ij,i+1)-v(ij, &
     1526          i))-mp(ij,i)*(vp(ij,i)-v(ij,i-1)))*dpinv
     1527        ! (saturated downdrafts resulting from mixing)               ! cld
     1528        DO k = i + 1, inb(ij) ! cld
     1529          qcond(ij, i) = qcond(ij, i) + elij(ij, k, i) ! cld
     1530          nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
     1531        END DO ! cld
     1532        ! (particular case: no detraining level is found)            ! cld
     1533        IF (nent(ij,i)==0) THEN ! cld
     1534          qcond(ij, i) = qcond(ij, i) + (1.-ep(ij,i))*clw(ij, i) ! cld
     1535          nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
     1536        END IF ! cld
     1537        IF (nqcond(ij,i)/=0.) THEN ! cld
     1538          qcond(ij, i) = qcond(ij, i)/nqcond(ij, i) ! cld
     1539        END IF ! cld
     1540      END IF
     1541    END DO
     15421500 END DO
     1543
     1544  ! *** Adjust tendencies at top of convection layer to reflect  ***
     1545  ! ***       actual position of the level zero cape             ***
     1546
     1547  DO ij = 1, ncum
     1548    fqold = fq(ij, inb(ij))
     1549    fq(ij, inb(ij)) = fq(ij, inb(ij))*(1.-frac(ij))
     1550    fq(ij, inb(ij)-1) = fq(ij, inb(ij)-1) + frac(ij)*fqold*((ph(ij, &
     1551      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, &
     1552      inb(ij))))*lv(ij, inb(ij))/lv(ij, inb(ij)-1)
     1553    ftold = ft(ij, inb(ij))
     1554    ft(ij, inb(ij)) = ft(ij, inb(ij))*(1.-frac(ij))
     1555    ft(ij, inb(ij)-1) = ft(ij, inb(ij)-1) + frac(ij)*ftold*((ph(ij, &
     1556      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, &
     1557      inb(ij))))*cpn(ij, inb(ij))/cpn(ij, inb(ij)-1)
     1558    fuold = fu(ij, inb(ij))
     1559    fu(ij, inb(ij)) = fu(ij, inb(ij))*(1.-frac(ij))
     1560    fu(ij, inb(ij)-1) = fu(ij, inb(ij)-1) + frac(ij)*fuold*((ph(ij, &
     1561      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
     1562    fvold = fv(ij, inb(ij))
     1563    fv(ij, inb(ij)) = fv(ij, inb(ij))*(1.-frac(ij))
     1564    fv(ij, inb(ij)-1) = fv(ij, inb(ij)-1) + frac(ij)*fvold*((ph(ij, &
     1565      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
     1566  END DO
     1567
     1568  ! ***   Very slightly adjust tendencies to force exact   ***
     1569  ! ***     enthalpy, momentum and tracer conservation     ***
     1570
     1571  DO ij = 1, ncum
     1572    ents(ij) = 0.0
     1573    uav(ij) = 0.0
     1574    vav(ij) = 0.0
     1575    DO i = 1, inb(ij)
     1576      ents(ij) = ents(ij) + (cpn(ij,i)*ft(ij,i)+lv(ij,i)*fq(ij,i))*(ph(ij,i)- &
     1577        ph(ij,i+1))
     1578      uav(ij) = uav(ij) + fu(ij, i)*(ph(ij,i)-ph(ij,i+1))
     1579      vav(ij) = vav(ij) + fv(ij, i)*(ph(ij,i)-ph(ij,i+1))
     1580    END DO
     1581  END DO
     1582  DO ij = 1, ncum
     1583    ents(ij) = ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
     1584    uav(ij) = uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
     1585    vav(ij) = vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
     1586  END DO
     1587  DO ij = 1, ncum
     1588    DO i = 1, inb(ij)
     1589      ft(ij, i) = ft(ij, i) - ents(ij)/cpn(ij, i)
     1590      fu(ij, i) = (1.-cu)*(fu(ij,i)-uav(ij))
     1591      fv(ij, i) = (1.-cu)*(fv(ij,i)-vav(ij))
     1592    END DO
     1593  END DO
     1594
     1595  DO k = 1, nl + 1
     1596    DO i = 1, ncum
     1597      IF ((q(i,k)+delt*fq(i,k))<0.0) iflag(i) = 10
     1598    END DO
     1599  END DO
     1600
     1601
     1602  DO i = 1, ncum
     1603    IF (iflag(i)>2) THEN
     1604      precip(i) = 0.0
     1605      cbmf(i) = 0.0
     1606    END IF
     1607  END DO
     1608  DO k = 1, nl
     1609    DO i = 1, ncum
     1610      IF (iflag(i)>2) THEN
     1611        ft(i, k) = 0.0
     1612        fq(i, k) = 0.0
     1613        fu(i, k) = 0.0
     1614        fv(i, k) = 0.0
     1615        qcondc(i, k) = 0.0 ! cld
     1616      END IF
     1617    END DO
     1618  END DO
     1619
     1620  DO k = 1, nl + 1
     1621    DO i = 1, ncum
     1622      ma(i, k) = 0.
     1623    END DO
     1624  END DO
     1625  DO k = nl, 1, -1
     1626    DO i = 1, ncum
     1627      ma(i, k) = ma(i, k+1) + m(i, k)
     1628    END DO
     1629  END DO
     1630
     1631
     1632  ! *** diagnose the in-cloud mixing ratio   ***            ! cld
     1633  ! ***           of condensed water         ***            ! cld
     1634  ! ! cld
     1635  DO ij = 1, ncum ! cld
     1636    DO i = 1, nd ! cld
     1637      mac(ij, i) = 0.0 ! cld
     1638      wa(ij, i) = 0.0 ! cld
     1639      siga(ij, i) = 0.0 ! cld
     1640    END DO ! cld
     1641    DO i = nk(ij), inb(ij) ! cld
     1642      DO k = i + 1, inb(ij) + 1 ! cld
     1643        mac(ij, i) = mac(ij, i) + m(ij, k) ! cld
     1644      END DO ! cld
     1645    END DO ! cld
     1646    DO i = icb(ij), inb(ij) - 1 ! cld
     1647      ax(ij, i) = 0. ! cld
     1648      DO j = icb(ij), i ! cld
     1649        ax(ij, i) = ax(ij, i) + rrd*(tvp(ij,j)-tv(ij,j)) & ! cld
     1650          *(ph(ij,j)-ph(ij,j+1))/p(ij, j) ! cld
     1651      END DO ! cld
     1652      IF (ax(ij,i)>0.0) THEN ! cld
     1653        wa(ij, i) = sqrt(2.*ax(ij,i)) ! cld
     1654      END IF ! cld
     1655    END DO ! cld
     1656    DO i = 1, nl ! cld
     1657      IF (wa(ij,i)>0.0) &          ! cld
     1658        siga(ij, i) = mac(ij, i)/wa(ij, i) & ! cld
     1659        *rrd*tvp(ij, i)/p(ij, i)/100./delta ! cld
     1660      siga(ij, i) = min(siga(ij,i), 1.0) ! cld
     1661      qcondc(ij, i) = siga(ij, i)*clw(ij, i)*(1.-ep(ij,i)) & ! cld
     1662        +(1.-siga(ij,i))*qcond(ij, i) ! cld
     1663    END DO ! cld
     1664  END DO ! cld
     1665
     1666  RETURN
     1667END SUBROUTINE cv_yield
     1668
     1669SUBROUTINE cv_uncompress(nloc, len, ncum, nd, idcum, iflag, precip, cbmf, ft, &
     1670    fq, fu, fv, ma, qcondc, iflag1, precip1, cbmf1, ft1, fq1, fu1, fv1, ma1, &
     1671    qcondc1)
     1672  IMPLICIT NONE
     1673
     1674  include "cvparam.h"
     1675
     1676  ! inputs:
     1677  INTEGER len, ncum, nd, nloc
     1678  INTEGER idcum(nloc)
     1679  INTEGER iflag(nloc)
     1680  REAL precip(nloc), cbmf(nloc)
     1681  REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
     1682  REAL ma(nloc, nd)
     1683  REAL qcondc(nloc, nd) !cld
     1684
     1685  ! outputs:
     1686  INTEGER iflag1(len)
     1687  REAL precip1(len), cbmf1(len)
     1688  REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
     1689  REAL ma1(len, nd)
     1690  REAL qcondc1(len, nd) !cld
     1691
     1692  ! local variables:
     1693  INTEGER i, k
     1694
     1695  DO i = 1, ncum
     1696    precip1(idcum(i)) = precip(i)
     1697    cbmf1(idcum(i)) = cbmf(i)
     1698    iflag1(idcum(i)) = iflag(i)
     1699  END DO
     1700
     1701  DO k = 1, nl
     1702    DO i = 1, ncum
     1703      ft1(idcum(i), k) = ft(i, k)
     1704      fq1(idcum(i), k) = fq(i, k)
     1705      fu1(idcum(i), k) = fu(i, k)
     1706      fv1(idcum(i), k) = fv(i, k)
     1707      ma1(idcum(i), k) = ma(i, k)
     1708      qcondc1(idcum(i), k) = qcondc(i, k)
     1709    END DO
     1710  END DO
     1711
     1712  RETURN
     1713END SUBROUTINE cv_uncompress
     1714
Note: See TracChangeset for help on using the changeset viewer.