Back to home page

MITgcm

 
 

    


File indexing completed on 2023-03-03 06:09:58 UTC

view on githubraw file Latest commit 06d4643e on 2023-01-18 18:18:37 UTC
a456aa407c Andr*0001 #include "FIZHI_OPTIONS.h"
9d62bca8f2 Andr*0002       subroutine moistio (ndmoist,istrip,npcs,
d2aaec7e02 Andr*0003      .   lowlevel,midlevel,nltop,nsubmin,nsubmax,Lup,
689620ef36 Andr*0004      .   pz,plz,plze,dpres,pkht,pkl,uz,vz,tz,qz,bi,bj,ntracerin,ptracer,
                0005      .   qqz,dumoist,dvmoist,dtmoist,dqmoist,cumfric,
d2aaec7e02 Andr*0006      .   im,jm,lm,ptop,
1662f365b2 Andr*0007      .   iras,rainlsp,rainconv,snowfall,
                0008      .   nswcld,cldtot_sw,cldras_sw,cldlsp_sw,nswlz,swlz,
                0009      .   nlwcld,cldtot_lw,cldras_lw,cldlsp_lw,nlwlz,lwlz,
7d59649d1b Andr*0010      .   lpnt,myid)
1662f365b2 Andr*0011 
7aae1e039f Andr*0012        implicit none
                0013 
a7a27977ab Jean*0014 C Input Variables
                0015 C ---------------
7aae1e039f Andr*0016       integer im,jm,lm
9d62bca8f2 Andr*0017       integer ndmoist,istrip,npcs
a7a27977ab Jean*0018       integer bi,bj,ntracerin,ptracer
d2aaec7e02 Andr*0019       integer lowlevel,midlevel,nltop,nsubmin,nsubmax,Lup
a456aa407c Andr*0020       _RL pz(im,jm),plz(im,jm,lm),plze(im,jm,lm+1),dpres(im,jm,lm)
                0021       _RL pkht(im,jm,lm+1),pkl(im,jm,lm)
a7a27977ab Jean*0022       _RL tz(im,jm,lm),qz(im,jm,lm,ntracerin)
                0023       _RL uz(im,jm,lm),vz(im,jm,lm)
a456aa407c Andr*0024       _RL qqz(im,jm,lm)
                0025       _RL dumoist(im,jm,lm),dvmoist(im,jm,lm)
689620ef36 Andr*0026       _RL dtmoist(im,jm,lm),dqmoist(im,jm,lm,ntracerin)
                0027       logical cumfric
a456aa407c Andr*0028       _RL ptop
9d62bca8f2 Andr*0029       integer iras
a456aa407c Andr*0030       _RL rainlsp(im,jm),rainconv(im,jm),snowfall(im,jm)
9d62bca8f2 Andr*0031       integer nswcld,nswlz
a456aa407c Andr*0032       _RL cldlsp_sw(im,jm,lm),cldras_sw(im,jm,lm)
                0033       _RL cldtot_sw(im,jm,lm),swlz(im,jm,lm)
9d62bca8f2 Andr*0034       integer nlwcld,nlwlz
a456aa407c Andr*0035       _RL  cldlsp_lw(im,jm,lm),cldras_lw(im,jm,lm)
                0036       _RL  cldtot_lw(im,jm,lm),lwlz(im,jm,lm)
9d62bca8f2 Andr*0037       logical lpnt
                0038       integer myid
1662f365b2 Andr*0039 
a7a27977ab Jean*0040 C Local Variables
                0041 C ---------------
d2aaec7e02 Andr*0042       integer    ncrnd,nsecf
1662f365b2 Andr*0043 
f0525ae611 Andr*0044       _RL       fracqq, rh, dum
d2aaec7e02 Andr*0045       integer    snowcrit
1662f365b2 Andr*0046       parameter (fracqq = 0.1)
daaf4a8d7f Andr*0047       _RL one
                0048       parameter (one=1.)
1662f365b2 Andr*0049 
a456aa407c Andr*0050       _RL   cldsr(im,jm,lm)
                0051       _RL   srcld(istrip,lm)
1662f365b2 Andr*0052 
a456aa407c Andr*0053       _RL plev
                0054       _RL cldnow,cldlsp_mem,cldlsp,cldras_mem,cldras
                0055       _RL watnow,watmin,cldmin
                0056       _RL cldprs(im,jm),cldtmp(im,jm)
                0057       _RL cldhi (im,jm),cldlow(im,jm)
                0058       _RL cldmid(im,jm),totcld(im,jm)
1662f365b2 Andr*0059 
a456aa407c Andr*0060       _RL   CLDLS(im,jm,lm)  , CPEN(im,jm,lm)
                0061       _RL    tmpimjm(im,jm)
                0062       _RL    lsp_new(im,jm)
                0063       _RL   conv_new(im,jm)
                0064       _RL   snow_new(im,jm)
1662f365b2 Andr*0065 
a456aa407c Andr*0066       _RL  qqcolmin(im,jm)
                0067       _RL  qqcolmax(im,jm)
1662f365b2 Andr*0068       integer levpbl(im,jm)
                0069 
a7a27977ab Jean*0070 C Gathered Arrays for Variable Cloud Base
                0071 C ---------------------------------------
a456aa407c Andr*0072       _RL    raincgath(im*jm)
                0073       _RL     pigather(im*jm)
                0074       _RL     thgather(im*jm,lm)
                0075       _RL     shgather(im*jm,lm)
                0076       _RL    pkzgather(im*jm,lm)
                0077       _RL    pkegather(im*jm,lm+1)
                0078       _RL    plzgather(im*jm,lm)
                0079       _RL    plegather(im*jm,lm+1)
                0080       _RL     dpgather(im*jm,lm)
                0081       _RL    tmpgather(im*jm,lm)
                0082       _RL   deltgather(im*jm,lm)
                0083       _RL   delqgather(im*jm,lm)
689620ef36 Andr*0084       _RL      ugather(im*jm,lm,ntracerin+2-ptracer)
                0085       _RL   delugather(im*jm,lm,ntracerin+2-ptracer)
a456aa407c Andr*0086       _RL     deltrnev(im*jm,lm)
                0087       _RL     delqrnev(im*jm,lm)
1662f365b2 Andr*0088 
                0089       integer  nindeces(lm)
                0090       integer  pblindex(im*jm)
                0091       integer levgather(im*jm)
                0092 
a7a27977ab Jean*0093 C Stripped Arrays
                0094 C ---------------
a456aa407c Andr*0095       _RL saveth (istrip,lm)
                0096       _RL saveq  (istrip,lm)
689620ef36 Andr*0097       _RL saveu  (istrip,lm,ntracerin+2-ptracer)
                0098       _RL usubcl (istrip,   ntracerin+2-ptracer)
a456aa407c Andr*0099 
                0100       _RL     ple(istrip,lm+1)
                0101       _RL      dp(istrip,lm)
                0102       _RL      TL(ISTRIP,lm)  , SHL(ISTRIP,lm)
                0103       _RL      PL(ISTRIP,lm)  , PLK(ISTRIP,lm)
                0104       _RL    PLKE(ISTRIP,lm+1)
                0105       _RL      TH(ISTRIP,lm)  ,CVTH(ISTRIP,lm)
                0106       _RL   CVQ(ISTRIP,lm)
689620ef36 Andr*0107       _RL      UL(ISTRIP,lm,ntracerin+2-ptracer)
                0108       _RL     cvu(istrip,lm,ntracerin+2-ptracer)
a456aa407c Andr*0109       _RL  CLMAXO(ISTRIP,lm),CLBOTH(ISTRIP,lm)
                0110       _RL  CLSBTH(ISTRIP,lm)
                0111       _RL    TMP1(ISTRIP,lm),  TMP2(ISTRIP,lm)
                0112       _RL    TMP3(ISTRIP,lm),  TMP4(ISTRIP,lm+1)
                0113       _RL    TMP5(ISTRIP,lm+1)
1662f365b2 Andr*0114       integer   ITMP1(ISTRIP,lm), ITMP2(ISTRIP,lm)
                0115 
a456aa407c Andr*0116       _RL   PRECIP(ISTRIP), PCNET(ISTRIP)
                0117       _RL   SP(ISTRIP),  PREP(ISTRIP)
                0118       _RL   PCPEN (ISTRIP,lm)
1662f365b2 Andr*0119       integer pbl(istrip),depths(lm)
                0120 
a456aa407c Andr*0121       _RL   cldlz(istrip,lm), cldwater(im,jm,lm)
                0122       _RL   rhfrac(istrip), rhmin, pup, ppbl, rhcrit(istrip,lm)
                0123       _RL   offset, alpha, rasmax
1662f365b2 Andr*0124 
                0125       logical first
                0126       logical lras
a456aa407c Andr*0127       _RL    clfrac (istrip,lm)
                0128       _RL    cldmas (istrip,lm)
                0129       _RL    detrain(istrip,lm)
                0130       _RL    psubcld    (istrip), psubcldg (im,jm)
                0131       _RL    psubcld_cnt(istrip), psubcldgc(im,jm)
                0132       _RL rnd(lm/2)
1662f365b2 Andr*0133       DATA      FIRST /.TRUE./
                0134 
d2aaec7e02 Andr*0135       integer imstp,nsubcl,nlras
471f7f004b Andr*0136       integer i,j,iloop,indx,indgath,l,nn,num,numdeps,nt
a456aa407c Andr*0137       _RL tmstp,tminv,sday,grav,alhl,cp,elocp,gamfac
                0138       _RL rkappa,p0kappa,p0kinv,ptopkap,pcheck
                0139       _RL tice,getcon,pi
689620ef36 Andr*0140       integer ntracer,ntracedim, ntracex
8b150b4af9 Andr*0141 
                0142 #ifdef ALLOW_DIAGNOSTICS
                0143       logical  diagnostics_is_on
                0144       external diagnostics_is_on
                0145       _RL tmpdiag(im,jm),tmpdiag2(im,jm)
                0146 #endif
a7a27977ab Jean*0147 
1662f365b2 Andr*0148 C **********************************************************************
                0149 C ****                     INITIALIZATION                           ****
                0150 C **********************************************************************
                0151 
689620ef36 Andr*0152 C Add U and V components to tracer array for cumulus friction
                0153 
                0154       if(cumfric) then
                0155        ntracer = ntracerin + 2
                0156       else
                0157        ntracer = ntracerin
                0158       endif
dfc4060574 Andr*0159       ntracedim= max(ntracer-ptracer,1)
689620ef36 Andr*0160       ntracex= ntracer-ptracer
1662f365b2 Andr*0161       IMSTP  = nsecf(NDMOIST)
                0162       TMSTP  = FLOAT(IMSTP)
                0163       TMINV  = 1. /  TMSTP
                0164 
                0165 C Minimum Large-Scale Cloud Fraction at rhcrit
a7a27977ab Jean*0166       alpha  = 0.80
7d59649d1b Andr*0167 C Difference in fraction between SR and LS Threshold
a7a27977ab Jean*0168       offset = 0.10
f0525ae611 Andr*0169 C Large-Scale Relative Humidity Threshold above 600 mb
a7a27977ab Jean*0170       rhmin  = 0.95
1662f365b2 Andr*0171 C Maximum Cloud Fraction associated with RAS
a7a27977ab Jean*0172       rasmax = 1.00
1662f365b2 Andr*0173 
f0525ae611 Andr*0174 cc        nn = 3*3600.0/tmstp + 1
                0175           nn = 3*1800.0/tmstp + 1
1662f365b2 Andr*0176 C Threshold for Cloud Fraction Memory
a7a27977ab Jean*0177 cc    cldmin = rasmax*(1.0-tmstp/3600.)**nn
                0178       cldmin = rasmax*(1.0-tmstp/1800.)**nn
1662f365b2 Andr*0179 C Threshold for Cloud Liquid Water Memory
a7a27977ab Jean*0180       watmin = 1.0e-8
1662f365b2 Andr*0181 
                0182       SDAY   = GETCON('SDAY')
                0183       GRAV   = GETCON('GRAVITY')
                0184       ALHL   = GETCON('LATENT HEAT COND')
                0185       CP     = GETCON('CP')
                0186       ELOCP  = GETCON('LATENT HEAT COND') / GETCON('CP')
                0187       GAMFAC = GETCON('LATENT HEAT COND') * GETCON('EPS') * ELOCP
                0188      .       / GETCON('RGAS')
                0189       RKAPPA = GETCON('KAPPA')
                0190       P0KAPPA  =  1000.0**RKAPPA
                0191       P0KINV   =  1. / P0KAPPA
                0192       PTOPKAP  =  PTOP**RKAPPA
                0193       tice     = getcon('FREEZING-POINT')
                0194       PI       = 4.*atan(1.)
                0195 
a7a27977ab Jean*0196 C Determine Total number of Random Clouds to Check
                0197 C ---------------------------------------------
c74d45a9d3 Andr*0198       ncrnd = (lm-nltop+1)/2
1662f365b2 Andr*0199 
97eb405d7b Andr*0200       if(first .and. myid.eq.1 .and. bi.eq.1 ) then
1662f365b2 Andr*0201        print *
d2aaec7e02 Andr*0202        print *,'Top Level Allowed for Convection : ',nltop
                0203        print *,'          Highest Sub-Cloud Level: ',nsubmax
                0204        print *,'          Lowest  Sub-Cloud Level: ',nsubmin
1662f365b2 Andr*0205        print *,'    Total Number of Random Clouds: ',ncrnd
                0206        print *
                0207        first = .false.
                0208       endif
                0209 
a7a27977ab Jean*0210 C And now find PBL depth - the level where qq = fracqq * qq at surface
                0211 C --------------------------------------------------------------------
1662f365b2 Andr*0212       do j = 1,jm
                0213       do i = 1,im
                0214        qqcolmin(i,j) = qqz(i,j,lm)*fracqq
                0215        qqcolmax(i,j) = qqz(i,j,lm)
                0216          levpbl(i,j) = lm+1
                0217       enddo
                0218       enddo
                0219 
                0220       DO L = lm-1,1,-1
                0221        DO j = 1,jm
                0222        DO i = 1,im
                0223         IF((qqz(i,j,l).gt.qqcolmax(i,j))
                0224      1                   .and.(levpbl(i,j).eq.lm+1))then
                0225          qqcolmax(i,j) = qqz(i,j,l)
                0226          qqcolmin(i,j) = fracqq*qqcolmax(i,j)
                0227         endif
                0228         if((qqz(i,j,l).lt.qqcolmin(i,j))
                0229      1                   .and.(levpbl(i,j).eq.lm+1))levpbl(i,j)=L+1
                0230        enddo
                0231        enddo
                0232       enddo
                0233 
                0234       do j = 1,jm
                0235       do i = 1,im
                0236         if(levpbl(i,j).gt.nsubmin) levpbl(i,j) = nsubmin
                0237         if(levpbl(i,j).lt.nsubmax) levpbl(i,j) = nsubmax
                0238       enddo
                0239       enddo
                0240 
a7a27977ab Jean*0241 C Set up the array of indeces of subcloud levels for the gathering
                0242 C ----------------------------------------------------------------
471f7f004b Andr*0243       indx = 0
1662f365b2 Andr*0244       do L = nsubmin,nltop,-1
                0245        do j = 1,jm
                0246        do i = 1,im
                0247         if(levpbl(i,j).eq.L) then
471f7f004b Andr*0248          indx = indx + 1
                0249          pblindex(indx) = (j-1)*im + i
1662f365b2 Andr*0250         endif
                0251        enddo
                0252        enddo
                0253       enddo
                0254 
471f7f004b Andr*0255       do indx = 1,im*jm
a7a27977ab Jean*0256 c      levgather(indx) = levpbl(pblindex(indx),1)
                0257 c       pigather(indx) =     pz(pblindex(indx),1)
                0258 c       pkegather(indx,lm+1) = pkht(pblindex(indx),1,lm+1)
                0259 c       plegather(indx,lm+1) = plze(pblindex(indx),1,lm+1)
                0260         j = 1+INT((pblindex(indx)-1)/im)
                0261         i = 1+MOD((pblindex(indx)-1),im)
                0262        levgather(indx) = levpbl(i,j)
                0263         pigather(indx) =     pz(i,j)
                0264         pkegather(indx,lm+1) = pkht(i,j,lm+1)
                0265         plegather(indx,lm+1) = plze(i,j,lm+1)
1662f365b2 Andr*0266       enddo
                0267 
                0268       do L = 1,lm
471f7f004b Andr*0269        do indx = 1,im*jm
a7a27977ab Jean*0270 c        thgather(indx,L) = tz(pblindex(indx),1,L)
                0271 c        shgather(indx,L) = qz(pblindex(indx),1,L,1)
                0272 c       pkegather(indx,L) = pkht(pblindex(indx),1,L)
                0273 c       pkzgather(indx,L) = pkl(pblindex(indx),1,L)
                0274 c       plegather(indx,L) = plze(pblindex(indx),1,L)
                0275 c       plzgather(indx,L) = plz(pblindex(indx),1,L)
                0276 c        dpgather(indx,L) = dpres(pblindex(indx),1,L)
                0277         j = 1+INT((pblindex(indx)-1)/im)
                0278         i = 1+MOD((pblindex(indx)-1),im)
                0279          thgather(indx,L) = tz(i,j,L)
                0280          shgather(indx,L) = qz(i,j,L,1)
                0281         pkegather(indx,L) = pkht(i,j,L)
                0282         pkzgather(indx,L) = pkl(i,j,L)
                0283         plegather(indx,L) = plze(i,j,L)
                0284         plzgather(indx,L) = plz(i,j,L)
                0285          dpgather(indx,L) = dpres(i,j,L)
1662f365b2 Andr*0286        enddo
                0287       enddo
689620ef36 Andr*0288 C General Tracers
                0289 C----------------
                0290       do nt = 1,ntracerin-ptracer
                0291       do L = 1,lm
                0292        do indx = 1,im*jm
a7a27977ab Jean*0293 c       ugather(indx,L,nt) = qz(pblindex(indx),1,L,nt+ptracer)
                0294         j = 1+INT((pblindex(indx)-1)/im)
                0295         i = 1+MOD((pblindex(indx)-1),im)
                0296         ugather(indx,L,nt) = qz(i,j,L,nt+ptracer)
689620ef36 Andr*0297        enddo
                0298       enddo
                0299       enddo
                0300 
                0301       if(cumfric) then
                0302 C Cumulus Friction - load u and v wind components into tracer array
                0303 C------------------------------------------------------------------
                0304       do L = 1,lm
                0305        do indx = 1,im*jm
a7a27977ab Jean*0306 c       ugather(indx,L,ntracerin-ptracer+1) = uz(pblindex(indx),1,L)
                0307 c       ugather(indx,L,ntracerin-ptracer+2) = vz(pblindex(indx),1,L)
                0308         j = 1+INT((pblindex(indx)-1)/im)
                0309         i = 1+MOD((pblindex(indx)-1),im)
                0310         ugather(indx,L,ntracerin-ptracer+1) = uz(i,j,L)
                0311         ugather(indx,L,ntracerin-ptracer+2) = vz(i,j,L)
689620ef36 Andr*0312        enddo
                0313       enddo
                0314 
                0315       endif
1662f365b2 Andr*0316 
a7a27977ab Jean*0317 C bump the counter for number of calls to convection
                0318 C --------------------------------------------------
1662f365b2 Andr*0319                         iras = iras + 1
a7a27977ab Jean*0320       if( iras.ge.1e9 ) iras = 1
1662f365b2 Andr*0321 
a7a27977ab Jean*0322 C select the 'random' cloud detrainment levels for RAS
                0323 C ----------------------------------------------------
1662f365b2 Andr*0324       call rndcloud(iras,ncrnd,rnd,myid)
a7a27977ab Jean*0325 
1662f365b2 Andr*0326       do l=1,lm
                0327       do j=1,jm
                0328       do i=1,im
689620ef36 Andr*0329       dumoist(i,j,l) = 0.
                0330       dvmoist(i,j,l) = 0.
1662f365b2 Andr*0331       dtmoist(i,j,l) = 0.
689620ef36 Andr*0332         do nt = 1,ntracerin
1662f365b2 Andr*0333         dqmoist(i,j,l,nt) = 0.
                0334         enddo
                0335       enddo
                0336       enddo
                0337       enddo
                0338 
                0339 C***********************************************************************
                0340 C ****                     LOOP OVER NPCS PEICES                    ****
                0341 C **********************************************************************
                0342 
                0343       DO 1000 NN = 1,NPCS
                0344 
                0345 C **********************************************************************
                0346 C ****                   VARIABLE INITIALIZATION                    ****
                0347 C **********************************************************************
                0348 
                0349        CALL STRIP (  pigather, SP      ,im*jm,ISTRIP,1 ,NN )
                0350        CALL STRIP ( pkzgather, PLK     ,im*jm,ISTRIP,lm,NN )
75f4744d22 Andr*0351        CALL STRIP ( pkegather, PLKE    ,im*jm,ISTRIP,lm+1,NN )
                0352        CALL STRIP ( plzgather, PL      ,im*jm,ISTRIP,lm,NN )
                0353        CALL STRIP ( plegather, PLE     ,im*jm,ISTRIP,lm+1,NN )
                0354        CALL STRIP (  dpgather, dp      ,im*jm,ISTRIP,lm,NN )
1662f365b2 Andr*0355        CALL STRIP (  thgather, TH      ,im*jm,ISTRIP,lm,NN )
                0356        CALL STRIP (  shgather, SHL     ,im*jm,ISTRIP,lm,NN )
                0357        CALL STRINT( levgather, pbl     ,im*jm,ISTRIP,1 ,NN )
                0358 
689620ef36 Andr*0359        do nt = 1,ntracer-ptracer
                0360         call strip ( ugather(1,1,nt), ul(1,1,nt),im*jm,istrip,lm,nn )
                0361        enddo
1662f365b2 Andr*0362 
                0363 C **********************************************************************
                0364 C ****        SETUP FOR RAS CUMULUS PARAMETERIZATION                ****
                0365 C **********************************************************************
f0525ae611 Andr*0366 C RAS works with real theta - convert from model theta
1662f365b2 Andr*0367       DO L = 1,lm
                0368       DO I = 1,ISTRIP
                0369              TH(I,L) = TH(I,L) * P0KAPPA
                0370          CLMAXO(I,L) = 0.
                0371          CLBOTH(I,L) = 0.
                0372          cldmas(I,L) = 0.
                0373         detrain(I,L) = 0.
                0374       ENDDO
                0375       ENDDO
                0376 
                0377       do L = 1,lm
                0378       depths(L) = 0
                0379       enddo
                0380 
                0381       numdeps = 0
                0382       do L = nsubmin,nltop,-1
                0383                        nindeces(L) = 0
                0384        do i = 1,istrip
                0385        if(pbl(i).eq.L) nindeces(L) = nindeces(L) + 1
                0386        enddo
                0387        if(nindeces(L).gt.0) then
                0388        numdeps = numdeps + 1
                0389        depths(numdeps) = L
                0390        endif
                0391       enddo
                0392 
                0393 C Initiate a do-loop around RAS for the number of different
                0394 C    sub-cloud layer depths in this strip
                0395 C --If all subcloud depths are the same, execute loop once
                0396 C   Otherwise loop over different subcloud layer depths
                0397 
                0398       num = 1
                0399       DO iloop = 1,numdeps
                0400 
                0401        nsubcl = depths(iloop)
                0402 
a7a27977ab Jean*0403 C Compute sub-cloud values for Temperature and Spec.Hum.
                0404 C ------------------------------------------------------
1662f365b2 Andr*0405        DO 600 I=num,num+nindeces(nsubcl)-1
                0406         TMP1(I,2) = 0.
                0407         TMP1(I,3) = 0.
                0408  600   CONTINUE
                0409 
                0410        NLRAS = NSUBCL - NLTOP + 1
                0411        DO 601 L=NSUBCL,lm
                0412        DO 602 I=num,num+nindeces(nsubcl)-1
                0413         TMP1(I,2) = TMP1(I,2) + (PLE(I,L+1)-PLE(I,L))*TH (I,L)/sp(i)
                0414         TMP1(I,3) = TMP1(I,3) + (PLE(I,L+1)-PLE(I,L))*SHL(I,L)/sp(i)
                0415  602   CONTINUE
                0416  601   CONTINUE
                0417        DO 603 I=num,num+nindeces(nsubcl)-1
                0418         TMP1(I,4) = 1. / ( (PLE(I,lm+1)-PLE(I,NSUBCL))/sp(I) )
                0419          TH(I,NSUBCL) = TMP1(I,2)*TMP1(I,4)
                0420         SHL(I,NSUBCL) = TMP1(I,3)*TMP1(I,4)
                0421  603   CONTINUE
                0422 
a7a27977ab Jean*0423 C Save initial value of tracers and compute sub-cloud value
                0424 C ---------------------------------------------------------
689620ef36 Andr*0425        DO NT = 1,ntracer-ptracer
                0426           do  L = 1,lm
                0427           do  i = num,num+nindeces(nsubcl)-1
                0428           saveu(i,L,nt) = ul(i,L,nt)
                0429           enddo
                0430           enddo
                0431           DO I=num,num+nindeces(nsubcl)-1
                0432           TMP1(I,2) = 0.
                0433           ENDDO
                0434           DO L=NSUBCL,lm
                0435           DO I=num,num+nindeces(nsubcl)-1
                0436            TMP1(I,2) = TMP1(I,2)+(PLE(I,L+1)-PLE(I,L))*UL(I,L,NT)/sp(i)
                0437           ENDDO
                0438           ENDDO
                0439           DO I=num,num+nindeces(nsubcl)-1
                0440           UL(I,NSUBCL,NT) = TMP1(I,2)*TMP1(I,4)
                0441              usubcl(i,nt) = ul(i,nsubcl,nt)
                0442           ENDDO
                0443        ENDDO
1662f365b2 Andr*0444 
a7a27977ab Jean*0445 C Compute Pressure Arrays for RAS
                0446 C -------------------------------
1662f365b2 Andr*0447        DO 111 L=1,lm
                0448        DO 112 I=num,num+nindeces(nsubcl)-1
                0449         TMP4(I,L) = PLE(I,L)
                0450  112   CONTINUE
                0451  111   CONTINUE
                0452        DO I=num,num+nindeces(nsubcl)-1
                0453         TMP5(I,1) = PTOPKAP / P0KAPPA
                0454        ENDDO
                0455        DO L=2,lm
                0456        DO I=num,num+nindeces(nsubcl)-1
75f4744d22 Andr*0457         TMP5(I,L) = PLKE(I,L)*P0KINV
1662f365b2 Andr*0458        ENDDO
                0459        ENDDO
                0460        DO  I=num,num+nindeces(nsubcl)-1
                0461         TMP4(I,lm+1) = PLE (I,lm+1)
75f4744d22 Andr*0462         TMP5(I,lm+1) = PLKE(I,lm+1)*P0KINV
1662f365b2 Andr*0463        ENDDO
                0464        DO 113 I=num,num+nindeces(nsubcl)-1
                0465         TMP4(I,NSUBCL+1) = PLE (I,lm+1)
75f4744d22 Andr*0466         TMP5(I,NSUBCL+1) = PLKE(I,lm+1)*P0KINV
1662f365b2 Andr*0467  113   CONTINUE
                0468 
                0469       do i=num,num+nindeces(nsubcl)-1
                0470 C Temperature at top of sub-cloud layer
a7a27977ab Jean*0471       tmp2(i,1) = TH(i,NSUBCL) * PLKE(i,NSUBCL)/P0KAPPA
1662f365b2 Andr*0472 C Pressure    at top of sub-cloud layer
a7a27977ab Jean*0473       tmp2(i,2) = tmp4(i,nsubcl)
1662f365b2 Andr*0474       enddo
                0475 
                0476       do i=num,num+nindeces(nsubcl)-1
f0525ae611 Andr*0477        call qsat (tmp2(i,1),tmp2(i,2),tmp2(i,3),dum,.false.)
                0478        rh =  SHL(I,NSUBCL) / tmp2(i,3)
                0479        if (rh .le. 0.85) then
                0480          rhfrac(i) = 0.
                0481        else if (rh .ge. 0.95) then
                0482          rhfrac(i) = 1.
                0483        else
                0484          rhfrac(i) = (rh-0.85)*10.
                0485        endif
1662f365b2 Andr*0486       enddo
f0525ae611 Andr*0487 CC  Uncomment out the previous code, comment out this to
                0488 CC          activate the RH threshold for convection (trigger)
                0489 CC    do i=num,num+nindeces(nsubcl)-1
                0490 CC      rhfrac(i) = 1.
                0491 CC    enddo
1662f365b2 Andr*0492 
                0493 C  Compute RH threshold for Large-scale condensation
                0494 C  Used in Slingo-Ritter clouds as well - define offset between SR and LS
                0495 
                0496 C Top level of atan func above this rh_threshold = rhmin
a7a27977ab Jean*0497       pup = 600.
1662f365b2 Andr*0498       do i=num,num+nindeces(nsubcl)-1
7d59649d1b Andr*0499        do L = nsubcl, lm
                0500          rhcrit(i,L) = 1.
                0501        enddo
                0502        do L = 1, nsubcl-1
75f4744d22 Andr*0503         pcheck = pl(i,L)
7d59649d1b Andr*0504         if (pcheck .le. pup) then
                0505          rhcrit(i,L) = rhmin
                0506         else
75f4744d22 Andr*0507          ppbl = pl(i,nsubcl)
a7a27977ab Jean*0508          rhcrit(i,L) = rhmin + (1.-rhmin)/(19.) *
7d59649d1b Andr*0509      .      ((atan( (2.*(pcheck-pup)/(ppbl-pup)-1.) *
a7a27977ab Jean*0510      .       tan(20.*pi/21.-0.5*pi) )
1662f365b2 Andr*0511      .        + 0.5*pi) * 21./pi - 1.)
7d59649d1b Andr*0512         endif
                0513        enddo
1662f365b2 Andr*0514       enddo
                0515 
a7a27977ab Jean*0516 C Save Initial Values of Temperature and Specific Humidity
                0517 C --------------------------------------------------------
1662f365b2 Andr*0518       do L = 1,lm
                0519       do i = num,num+nindeces(nsubcl)-1
                0520         saveth(i,L) = th (i,L)
                0521         saveq (i,L) = shl(i,L)
                0522         PCPEN (i,L) = 0.
                0523         CLFRAC(i,L) = 0.
                0524       enddo
                0525       enddo
                0526 
                0527       CALL RAS ( NN,istrip,nindeces(nsubcl),NLRAS,NLTOP,lm,TMSTP
689620ef36 Andr*0528      1, UL(num,1,1),ntracedim,ntracex,TH(num,NLTOP),SHL(num,NLTOP)
1662f365b2 Andr*0529      2, TMP4(num,NLTOP), TMP5(num,NLTOP),rnd, ncrnd, PCPEN(num,NLTOP)
                0530      3, CLBOTH(num,NLTOP), CLFRAC(num,NLTOP)
                0531      4, cldmas(num,nltop), detrain(num,nltop)
                0532      8, cp,grav,rkappa,alhl,rhfrac(num),rasmax )
                0533 
a7a27977ab Jean*0534 C Compute Diagnostic CLDMAS in RAS Subcloud Layers
                0535 C ------------------------------------------------
1662f365b2 Andr*0536        do L=nsubcl,lm
                0537        do I=num,num+nindeces(nsubcl)-1
75f4744d22 Andr*0538         dum = dp(i,L)/(ple(i,lm+1)-ple(i,nsubcl))
1662f365b2 Andr*0539         cldmas(i,L) = cldmas(i,L-1) - dum*cldmas(i,nsubcl-1)
                0540        enddo
                0541        enddo
                0542 
a7a27977ab Jean*0543 C Update Theta and Moisture due to RAS
                0544 C ------------------------------------
1662f365b2 Andr*0545        DO L=1,nsubcl
                0546        DO I=num,num+nindeces(nsubcl)-1
                0547         CVTH(I,L) =  (TH (I,L) - saveth(i,l))
                0548         CVQ (I,L) =  (SHL(I,L) - saveq (i,l))
                0549        ENDDO
                0550        ENDDO
                0551        DO L=nsubcl+1,lm
                0552        DO I=num,num+nindeces(nsubcl)-1
                0553         CVTH(I,L) = cvth(i,nsubcl)
                0554         CVQ (I,L) = cvq (i,nsubcl)
                0555        ENDDO
                0556        ENDDO
                0557 
                0558        DO L=nsubcl+1,lm
                0559        DO I=num,num+nindeces(nsubcl)-1
                0560         TH (I,L) = saveth(i,l) + cvth(i,l)
                0561         SHL(I,L) = saveq (i,l) + cvq (i,l)
                0562        ENDDO
                0563        ENDDO
                0564        DO L=1,lm
                0565        DO I=num,num+nindeces(nsubcl)-1
                0566         CVTH(I,L) =  CVTH(I,L) *P0KINV*SP(I)*tminv
                0567         CVQ (I,L) =  CVQ (I,L) *SP(I)*tminv
                0568        ENDDO
                0569        ENDDO
                0570 
a7a27977ab Jean*0571 C Compute Tracer Tendency due to RAS
                0572 C ----------------------------------
689620ef36 Andr*0573        do nt = 1,ntracer-ptracer
                0574         DO L=1,nsubcl-1
                0575         DO I=num,num+nindeces(nsubcl)-1
                0576          CVU(I,L,nt) = ( UL(I,L,nt)-saveu(i,l,nt) )*sp(i)*tminv
                0577         ENDDO
                0578         ENDDO
                0579         DO L=nsubcl,lm
                0580         DO I=num,num+nindeces(nsubcl)-1
                0581          if( usubcl(i,nt).ne.0.0 ) then
a7a27977ab Jean*0582           cvu(i,L,nt) = ( ul(i,nsubcl,nt)-usubcl(i,nt) ) *
689620ef36 Andr*0583      .                     ( saveu(i,L,nt)/usubcl(i,nt) )*sp(i)*tminv
                0584          else
                0585           cvu(i,L,nt) = 0.0
                0586          endif
                0587         ENDDO
                0588         ENDDO
                0589        enddo
1662f365b2 Andr*0590 
a7a27977ab Jean*0591 C Compute Diagnostic PSUBCLD (Subcloud Layer Pressure)
                0592 C ----------------------------------------------------
1662f365b2 Andr*0593        do i=num,num+nindeces(nsubcl)-1
                0594                              lras = .false.
                0595        do L=nltop,nsubcl
                0596        if( cvq(i,L).ne.0.0 ) lras = .true.
                0597        enddo
                0598        psubcld    (i) = 0.0
                0599        psubcld_cnt(i) = 0.0
a7a27977ab Jean*0600        if( lras ) then
1662f365b2 Andr*0601        psubcld    (i) = sp(i)+ptop-ple(i,nsubcl)
                0602        psubcld_cnt(i) = 1.0
                0603        endif
                0604        enddo
a7a27977ab Jean*0605 
1662f365b2 Andr*0606 C End of subcloud layer depth loop (iloop)
                0607 
                0608        num = num+nindeces(nsubcl)
                0609 
                0610       ENDDO
                0611 
                0612 C **********************************************************************
                0613 C ****                     TENDENCY UPDATES                         ****
                0614 C ****    (Keep 'Gathered' tendencies in 'gather' arrays now)       ****
                0615 C **********************************************************************
                0616 
                0617       call paste( CVTH,deltgather,istrip,im*jm,lm,NN )
                0618       call paste(  CVQ,delqgather,istrip,im*jm,lm,NN )
689620ef36 Andr*0619       do nt = 1,ntracer-ptracer
                0620       call paste( CVU(1,1,nt),delugather(1,1,nt),istrip,im*jm,lm,NN )
                0621       enddo
1662f365b2 Andr*0622 
                0623 C **********************************************************************
                0624 C     And now paste some arrays for filling diagnostics
                0625 C (use pkegather to hold detrainment and tmpgather for cloud mass flux)
                0626 C **********************************************************************
                0627 
8b150b4af9 Andr*0628       call paste( cldmas,tmpgather,istrip,im*jm,lm,NN)
                0629       call paste(detrain,pkegather,istrip,im*jm,lm,NN)
1662f365b2 Andr*0630       call paste(psubcld    ,psubcldg ,istrip,im*jm,1,NN)
                0631       call paste(psubcld_cnt,psubcldgc,istrip,im*jm,1,NN)
                0632 
                0633 C *********************************************************************
                0634 C ****         RE-EVAPORATION OF PENETRATING CONVECTIVE RAIN       ****
                0635 C *********************************************************************
a7a27977ab Jean*0636 
1662f365b2 Andr*0637       CALL STRIP ( thgather,TH ,im*jm,ISTRIP,lm,NN)
                0638       CALL STRIP ( shgather,SHL,im*jm,ISTRIP,lm,NN)
                0639       DO L=1,lm
                0640       DO I=1,ISTRIP
                0641            TH(I,L) =  TH(I,L) + CVTH(I,L)*tmstp/SP(I)
                0642           SHL(I,L) = SHL(I,L) +  CVQ(I,L)*tmstp/SP(I)
                0643            TL(I,L) =  TH(I,L)*PLK(I,L)
                0644        saveth(I,L) =  th(I,L)
                0645        saveq (I,L) = SHL(I,L)
                0646       ENDDO
                0647       ENDDO
a7a27977ab Jean*0648 
75f4744d22 Andr*0649        CALL RNEVP (NN,ISTRIP,lm,TL,SHL,PCPEN,PL,CLFRAC,SP,DP,PLKE,
1662f365b2 Andr*0650      .  PLK,TH,TMP1,TMP2,TMP3,ITMP1,ITMP2,PCNET,PRECIP,
daaf4a8d7f Andr*0651      .  CLSBTH,TMSTP,one,cp,grav,alhl,gamfac,cldlz,rhcrit,offset,alpha)
1662f365b2 Andr*0652 
                0653 C **********************************************************************
                0654 C ****                     TENDENCY UPDATES                         ****
                0655 C **********************************************************************
                0656 
                0657       DO L=1,lm
                0658 
                0659        DO I =1,ISTRIP
                0660        TMP1(I,L) = sp(i) * (SHL(I,L)-saveq(I,L)) * tminv
                0661        ENDDO
                0662        CALL PSTBMP(TMP1(1,L),delqgather(1,L),ISTRIP,im*jm,1,NN)
                0663 
                0664        DO I =1,ISTRIP
                0665        TMP1(I,L) = sp(i) * ((TL(I,L)/PLK(I,L))-saveth(i,l)) * tminv
                0666        ENDDO
                0667        CALL PSTBMP(TMP1(1,L),deltgather(1,L),ISTRIP,im*jm,1,NN)
                0668 
                0669 C Paste rain evap tendencies into arrays for diagnostic output
a7a27977ab Jean*0670 C ------------------------------------------------------------
1662f365b2 Andr*0671        DO I =1,ISTRIP
                0672         TMP1(I,L) = ((TL(I,L)/PLK(I,L))-saveth(i,l))*plk(i,l)*sday*tminv
                0673        ENDDO
                0674        call paste(tmp1(1,L),deltrnev(1,L),istrip,im*jm,1,NN)
                0675 
                0676        DO I =1,ISTRIP
                0677         TMP1(I,L) = (SHL(I,L)-saveq(I,L)) * 1000. * sday * tminv
                0678        ENDDO
                0679        call paste(tmp1(1,L),delqrnev(1,L),istrip,im*jm,1,NN)
                0680 
                0681       ENDDO
                0682 
                0683 C *********************************************************************
a7a27977ab Jean*0684 C          Add Non-Precipitating Clouds where the relative
1662f365b2 Andr*0685 C                     humidity is less than 100%
                0686 C               Apply Cloud Top Entrainment Instability
                0687 C *********************************************************************
                0688 
                0689       do L=1,lm
                0690       do i=1,istrip
                0691       srcld(i,L) = -clsbth(i,L)
                0692       enddo
                0693       enddo
                0694 
                0695       call srclouds (saveth,saveq,plk,pl,plke,clsbth,cldlz,istrip,lm,
                0696      .                                              rhcrit,offset,alpha)
                0697 
                0698       do L=1,lm
                0699       do i=1,istrip
                0700       srcld(i,L) = srcld(i,L)+clsbth(i,L)
                0701       enddo
                0702       enddo
                0703 
                0704 C *********************************************************************
                0705 C ****                PASTE CLOUD AMOUNTS                          ****
                0706 C *********************************************************************
                0707 
                0708       call paste (  srcld,   cldsr,istrip,im*jm,lm,nn )
                0709       call paste (  cldlz,cldwater,istrip,im*jm,lm,nn )
                0710       call paste ( clsbth,   cldls,istrip,im*jm,lm,nn )
                0711       call paste ( clboth,   cpen ,istrip,im*jm,lm,nn )
a7a27977ab Jean*0712 
                0713 C compute Total Accumulated Precip for Landsurface Model
                0714 C ------------------------------------------------------
1662f365b2 Andr*0715       do i = 1,istrip
                0716 C  Initialize Rainlsp, Rainconv and Snowfall
a7a27977ab Jean*0717       tmp1(i,1) = 0.0
                0718       tmp1(i,2) = 0.0
1662f365b2 Andr*0719       tmp1(i,3) = 0.0
                0720       enddo
                0721 
                0722       do i = 1,istrip
                0723       prep(i)   = PRECIP(I) + PCNET(I)
                0724       tmp1(i,1) = PRECIP(I)
                0725       tmp1(i,2) = pcnet(i)
                0726       enddo
a7a27977ab Jean*0727 C
                0728 C  check whether there is snow
                0729 C-------------------------------------------------------
                0730 C  snow algorthm:
                0731 C  if temperature profile from the surface level to 700 mb
                0732 C  uniformaly c  below zero, then precipitation (total) is
                0733 C  snowfall.  Else there is no snow.
                0734 C-------------------------------------------------------
1662f365b2 Andr*0735 
                0736         do i = 1,istrip
                0737           snowcrit=0
                0738           do l=lup,lm
                0739             if (saveth(i,l)*plk(i,l).le. tice ) then
                0740              snowcrit=snowcrit+1
                0741             endif
                0742           enddo
                0743           if (snowcrit .eq. (lm-lup+1)) then
                0744             tmp1(i,3) = prep(i)
                0745             tmp1(i,1)=0.0
                0746             tmp1(i,2)=0.0
                0747           endif
                0748         enddo
a7a27977ab Jean*0749 
1662f365b2 Andr*0750       CALL paste (tmp1(1,1), lsp_new,ISTRIP,im*jm,1,NN)
                0751       CALL paste (tmp1(1,2),conv_new,ISTRIP,im*jm,1,NN)
                0752       CALL paste (tmp1(1,3),snow_new,ISTRIP,im*jm,1,NN)
                0753       CALL paste (pcnet,raincgath,ISTRIP,im*jm,1,NN)
                0754 
                0755 C *********************************************************************
                0756 C ****               End Major Stripped Region                     ****
                0757 C *********************************************************************
                0758 
                0759  1000 CONTINUE
                0760 
                0761 C Large Scale Rainfall, Conv rain, and snowfall
a7a27977ab Jean*0762 C ---------------------------------------------
1662f365b2 Andr*0763       call back2grd ( lsp_new,pblindex, lsp_new,im*jm)
                0764       call back2grd (conv_new,pblindex,conv_new,im*jm)
                0765       call back2grd (snow_new,pblindex,snow_new,im*jm)
                0766       call back2grd (raincgath,pblindex,raincgath,im*jm)
                0767 
a7a27977ab Jean*0768 C Subcloud Layer Pressure
                0769 C -----------------------
1662f365b2 Andr*0770       call back2grd (psubcldg ,pblindex,psubcldg ,im*jm)
                0771       call back2grd (psubcldgc,pblindex,psubcldgc,im*jm)
                0772 
                0773       do L = 1,lm
                0774 C Delta theta,q, convective, max and ls clouds
a7a27977ab Jean*0775 C --------------------------------------------
1662f365b2 Andr*0776        call back2grd (deltgather(1,L),pblindex, dtmoist(1,1,L)  ,im*jm)
                0777        call back2grd (delqgather(1,L),pblindex, dqmoist(1,1,L,1),im*jm)
                0778        call back2grd (    cpen(1,1,L),pblindex,    cpen(1,1,L)  ,im*jm)
                0779        call back2grd (   cldls(1,1,L),pblindex,   cldls(1,1,L)  ,im*jm)
                0780        call back2grd (cldwater(1,1,L),pblindex,cldwater(1,1,L)  ,im*jm)
                0781        call back2grd ( pkzgather(1,L),pblindex, pkzgather(1,L)  ,im*jm)
                0782 
                0783 C Diagnostics:
a7a27977ab Jean*0784 C ------------
8b150b4af9 Andr*0785        call back2grd(tmpgather(1,L),pblindex,
1662f365b2 Andr*0786      .                                            tmpgather(1,L),im*jm)
8b150b4af9 Andr*0787        call back2grd(pkegather(1,L),pblindex,
1662f365b2 Andr*0788      .                                            pkegather(1,L),im*jm)
8b150b4af9 Andr*0789        call back2grd(deltrnev(1,L),pblindex,
1662f365b2 Andr*0790      .                                             deltrnev(1,L),im*jm)
8b150b4af9 Andr*0791        call back2grd(delqrnev(1,L),pblindex,
1662f365b2 Andr*0792      .                                             delqrnev(1,L),im*jm)
8b150b4af9 Andr*0793        call back2grd(cldsr(1,1,L),pblindex,
1662f365b2 Andr*0794      .                                              cldsr(1,1,L),im*jm)
                0795       enddo
                0796 
a7a27977ab Jean*0797 C General Tracers
                0798 C ---------------
689620ef36 Andr*0799       do nt = 1,ntracerin-ptracer
                0800        do L = 1,lm
                0801        call back2grd (delugather(1,L,nt),pblindex,
                0802      .                                 dqmoist(1,1,L,ptracer+nt),im*jm)
                0803        enddo
                0804       enddo
                0805 
                0806       if(cumfric) then
                0807 
a7a27977ab Jean*0808 C U and V for cumulus friction
                0809 C ----------------------------
689620ef36 Andr*0810       do L = 1,lm
                0811        call back2grd (delugather(1,L,ntracerin-ptracer+1),pblindex,
                0812      .                                 dumoist(1,1,L),im*jm)
                0813        call back2grd (delugather(1,L,ntracerin-ptracer+2),pblindex,
                0814      .                                 dvmoist(1,1,L),im*jm)
                0815       enddo
                0816 
                0817 C Remove pi-weighting for u and v tendencies
                0818       do j = 1,jm
                0819       do i = 1,im
                0820        tmpimjm(i,j) = 1./pz(i,j)
                0821       enddo
                0822       enddo
                0823       do L = 1,lm
                0824       do j = 1,jm
                0825       do i = 1,im
                0826        dumoist(i,j,L) = dumoist(i,j,L) * tmpimjm(i,j)
889be31914 Andr*0827        dvmoist(i,j,L) = dvmoist(i,j,L) * tmpimjm(i,j)
689620ef36 Andr*0828       enddo
                0829       enddo
                0830       enddo
1662f365b2 Andr*0831 
689620ef36 Andr*0832       endif
1662f365b2 Andr*0833 
                0834 C **********************************************************************
                0835 C                          BUMP DIAGNOSTICS
                0836 C **********************************************************************
                0837 
8b150b4af9 Andr*0838 #ifdef ALLOW_DIAGNOSTICS
471f7f004b Andr*0839 
a7a27977ab Jean*0840 C Sub-Cloud Layer
                0841 C -------------------------
8b150b4af9 Andr*0842       if(diagnostics_is_on('PSUBCLD ',myid) .and.
                0843      .                 diagnostics_is_on('PSUBCLDC',myid) ) then
                0844        call diagnostics_fill(psubcldg,'PSUBCLD ',0,1,3,bi,bj,myid)
                0845        call diagnostics_fill(psubcldgc,'PSUBCLDC',0,1,3,bi,bj,myid)
1662f365b2 Andr*0846       endif
8b150b4af9 Andr*0847 
a7a27977ab Jean*0848 C Non-Precipitating Cloud Fraction
                0849 C --------------------------------
8b150b4af9 Andr*0850       if(diagnostics_is_on('CLDNP   ',myid) ) then
                0851        do L=1,lm
                0852        do j=1,jm
                0853        do i=1,im
                0854         tmpdiag(i,j) = cldsr(i,j,L)
                0855        enddo
                0856        enddo
                0857        call diagnostics_fill(tmpdiag,'CLDNP   ',L,1,3,bi,bj,myid)
1662f365b2 Andr*0858       enddo
                0859       endif
                0860 
a7a27977ab Jean*0861 C Moist Processes Heating Rate
                0862 C ----------------------------
8b150b4af9 Andr*0863       if(diagnostics_is_on('MOISTT  ',myid) ) then
                0864        do L=1,lm
                0865        do j=1,jm
                0866        do i=1,im
                0867         indgath = (j-1)*im + i
                0868         tmpdiag(i,j)=(dtmoist(i,j,L)*sday*pkzgather(indgath,L)/pz(i,j))
                0869        enddo
                0870        enddo
                0871        call diagnostics_fill(tmpdiag,'MOISTT  ',L,1,3,bi,bj,myid)
1662f365b2 Andr*0872       enddo
                0873       endif
                0874 
a7a27977ab Jean*0875 C Moist Processes Moistening Rate
                0876 C -------------------------------
8b150b4af9 Andr*0877       if(diagnostics_is_on('MOISTQ  ',myid) ) then
                0878        do L=1,lm
                0879        do j=1,jm
                0880        do i=1,im
                0881         tmpdiag(i,j)=(dqmoist(i,j,L,1)*sday*1000./pz(i,j))
                0882        enddo
                0883        enddo
                0884        call diagnostics_fill(tmpdiag,'MOISTQ  ',L,1,3,bi,bj,myid)
1662f365b2 Andr*0885       enddo
                0886       endif
                0887 
a7a27977ab Jean*0888 C Moist Processes U-tendency
                0889 C ----------------------------
f0525ae611 Andr*0890       if(diagnostics_is_on('MOISTU  ',myid) ) then
                0891        do L=1,lm
                0892        do j=1,jm
                0893        do i=1,im
                0894         tmpdiag(i,j)=dumoist(i,j,L)*sday
                0895        enddo
                0896        enddo
                0897        call diagnostics_fill(tmpdiag,'MOISTU  ',L,1,3,bi,bj,myid)
                0898       enddo
                0899       endif
                0900 
a7a27977ab Jean*0901 C Moist Processes V-tendency
                0902 C ----------------------------
f0525ae611 Andr*0903       if(diagnostics_is_on('MOISTV  ',myid) ) then
                0904        do L=1,lm
                0905        do j=1,jm
                0906        do i=1,im
                0907         tmpdiag(i,j)=dvmoist(i,j,L)*sday
                0908        enddo
                0909        enddo
                0910        call diagnostics_fill(tmpdiag,'MOISTV  ',L,1,3,bi,bj,myid)
                0911       enddo
                0912       endif
                0913 
a7a27977ab Jean*0914 C Cloud Mass Flux
                0915 C ---------------
8b150b4af9 Andr*0916       if(diagnostics_is_on('CLDMAS  ',myid) ) then
                0917        do L=1,lm
                0918        do j=1,jm
                0919        do i=1,im
                0920         indgath = (j-1)*im + i
                0921         tmpdiag(i,j)=tmpgather(indgath,L)
                0922        enddo
                0923        enddo
                0924        call diagnostics_fill(tmpdiag,'CLDMAS  ',L,1,3,bi,bj,myid)
                0925        enddo
1662f365b2 Andr*0926       endif
                0927 
a7a27977ab Jean*0928 C Detrained Cloud Mass Flux
                0929 C -------------------------
8b150b4af9 Andr*0930       if(diagnostics_is_on('DTRAIN  ',myid) ) then
                0931        do L=1,lm
                0932        do j=1,jm
                0933        do i=1,im
                0934         indgath = (j-1)*im + i
                0935         tmpdiag(i,j)=pkegather(indgath,L)
                0936        enddo
                0937        enddo
                0938        call diagnostics_fill(tmpdiag,'DTRAIN  ',L,1,3,bi,bj,myid)
                0939        enddo
1662f365b2 Andr*0940       endif
                0941 
a7a27977ab Jean*0942 C Grid-Scale Condensational Heating Rate
                0943 C --------------------------------------
8b150b4af9 Andr*0944       if(diagnostics_is_on('DTLS    ',myid) ) then
                0945        do L=1,lm
                0946        do j=1,jm
                0947        do i=1,im
                0948         indgath = (j-1)*im + i
                0949         tmpdiag(i,j)=deltrnev(indgath,L)
                0950        enddo
                0951        enddo
                0952        call diagnostics_fill(tmpdiag,'DTLS    ',L,1,3,bi,bj,myid)
                0953        enddo
1662f365b2 Andr*0954       endif
                0955 
a7a27977ab Jean*0956 C Grid-Scale Condensational Moistening Rate
                0957 C -----------------------------------------
8b150b4af9 Andr*0958       if(diagnostics_is_on('DQLS    ',myid) ) then
                0959        do L=1,lm
                0960        do j=1,jm
                0961        do i=1,im
                0962         indgath = (j-1)*im + i
                0963         tmpdiag(i,j)=delqrnev(indgath,L)
                0964        enddo
                0965        enddo
                0966        call diagnostics_fill(tmpdiag,'DQLS    ',L,1,3,bi,bj,myid)
                0967        enddo
1662f365b2 Andr*0968       endif
                0969 
a7a27977ab Jean*0970 C Total Precipitation
                0971 C -------------------
8b150b4af9 Andr*0972       if(diagnostics_is_on('PREACC  ',myid) ) then
                0973        do j=1,jm
                0974        do i=1,im
a7a27977ab Jean*0975         tmpdiag(i,j) = (lsp_new(I,j) + snow_new(I,j) + conv_new(i,j))
8b150b4af9 Andr*0976      .                                                    *sday*tminv
                0977        enddo
                0978        enddo
ab43cdb872 Jean*0979        call diagnostics_fill(tmpdiag,'PREACC  ',0,1,3,bi,bj,myid)
1662f365b2 Andr*0980       endif
                0981 
a7a27977ab Jean*0982 C Convective Precipitation
                0983 C ------------------------
8b150b4af9 Andr*0984       if(diagnostics_is_on('PRECON  ',myid) ) then
                0985        do j=1,jm
                0986        do i=1,im
                0987         indgath = (j-1)*im + i
                0988         tmpdiag(i,j) = raincgath(indgath)*sday*tminv
                0989        enddo
                0990        enddo
ab43cdb872 Jean*0991        call diagnostics_fill(tmpdiag,'PRECON  ',0,1,3,bi,bj,myid)
1662f365b2 Andr*0992       endif
                0993 
8b150b4af9 Andr*0994 #endif
                0995 
1662f365b2 Andr*0996 C **********************************************************************
                0997 C ****   Fill Rainfall and Snowfall Arrays for Land Surface Model   ****
                0998 C ****        Note:  Precip Rates work when DT(turb)<DT(moist)      ****
                0999 C **********************************************************************
                1000 
                1001       do j = 1,jm
                1002       do i = 1,im
                1003        rainlsp (i,j) = rainlsp (i,j) +  lsp_new(i,j)*tminv
                1004        rainconv(i,j) = rainconv(i,j) + conv_new(i,j)*tminv
                1005        snowfall(i,j) = snowfall(i,j) + snow_new(i,j)*tminv
                1006       enddo
                1007       enddo
                1008 
                1009 C **********************************************************************
                1010 C ***        Compute Time-averaged Quantities for Radiation          ***
                1011 C ***               CPEN  => Cloud Fraction from RAS                 ***
                1012 C ***               CLDLS => Cloud Fraction from RNEVP               ***
                1013 C **********************************************************************
                1014 
                1015       do j = 1,jm
                1016       do i = 1,im
                1017           cldhi (i,j) = 0.
                1018           cldmid(i,j) = 0.
                1019           cldlow(i,j) = 0.
                1020           cldtmp(i,j) = 0.
                1021           cldprs(i,j) = 0.
                1022          tmpimjm(i,j) = 0.
                1023       enddo
                1024       enddo
                1025 
a7a27977ab Jean*1026 C Set Moist-Process Memory Coefficient
                1027 C ------------------------------------
                1028       cldras_mem = 1.0-tmstp/ 3600.0
f0525ae611 Andr*1029       cldras_mem = 1.0-tmstp/ 1800.0
a7a27977ab Jean*1030       cldlsp_mem = 1.0-tmstp/(3600.0*3)
f0525ae611 Andr*1031       cldlsp_mem = 1.0-tmstp/(1800.0)
1662f365b2 Andr*1032 
                1033       do L = 1,lm
a7a27977ab Jean*1034       do j = 1,jm
                1035       do i = 1,im
                1036 C- to recover old code, replace the 2 lines above by the 2 commented below:
                1037 c     do j = 1,1
                1038 c     do i = 1,im*jm
                1039        indx = (j-1)*im + i
                1040        plev = pl(indx,L)
                1041 
                1042 C Compute Time-averaged Cloud and Water Amounts for Longwave Radiation
                1043 C --------------------------------------------------------------------
                1044        watnow = cldwater(i,j,L)
1662f365b2 Andr*1045        if( plev.le.500.0 ) then
a7a27977ab Jean*1046           cldras = min( max( cldras_lw(i,j,L)*cldras_mem,cpen(i,j,L)),
e7070f537c Cons*1047      $         1.0 _d 0)
1662f365b2 Andr*1048        else
                1049         cldras = 0.0
                1050        endif
a7a27977ab Jean*1051        cldlsp = min( max( cldlsp_lw(i,j,L)*cldlsp_mem,cldls(i,j,L)),
e7070f537c Cons*1052      $      1.0 _d 0)
1662f365b2 Andr*1053 
                1054        if( cldras.lt.cldmin ) cldras = 0.0
                1055        if( cldlsp.lt.cldmin ) cldlsp = 0.0
                1056 
                1057        cldnow = max( cldlsp,cldras )
                1058 
a7a27977ab Jean*1059        lwlz(i,j,L) = ( nlwlz*lwlz(i,j,L) + watnow)/(nlwlz +1)
                1060        cldtot_lw(i,j,L) = (nlwcld*cldtot_lw(i,j,L) + cldnow)/(nlwcld+1)
                1061        cldlsp_lw(i,j,L) = (nlwcld*cldlsp_lw(i,j,L) + cldlsp)/(nlwcld+1)
                1062        cldras_lw(i,j,L) = (nlwcld*cldras_lw(i,j,L) + cldras)/(nlwcld+1)
                1063 
                1064 C Compute Time-averaged Cloud and Water Amounts for Shortwave Radiation
                1065 C ---------------------------------------------------------------------
                1066        watnow = cldwater(i,j,L)
1662f365b2 Andr*1067        if( plev.le.500.0 ) then
a7a27977ab Jean*1068           cldras = min( max(cldras_sw(i,j,L)*cldras_mem, cpen(i,j,L)),
e7070f537c Cons*1069      $         1.0 _d 0)
1662f365b2 Andr*1070        else
                1071         cldras = 0.0
                1072        endif
a7a27977ab Jean*1073        cldlsp = min( max(cldlsp_sw(i,j,L)*cldlsp_mem,cldls(i,j,L)),
e7070f537c Cons*1074      $      1.0 _d 0)
1662f365b2 Andr*1075 
                1076        if( cldras.lt.cldmin ) cldras = 0.0
                1077        if( cldlsp.lt.cldmin ) cldlsp = 0.0
                1078 
                1079        cldnow = max( cldlsp,cldras )
                1080 
a7a27977ab Jean*1081        swlz(i,j,L) = ( nswlz*swlz(i,j,L) + watnow)/(nswlz +1)
                1082        cldtot_sw(i,j,L) = (nswcld*cldtot_sw(i,j,L) + cldnow)/(nswcld+1)
                1083        cldlsp_sw(i,j,L) = (nswcld*cldlsp_sw(i,j,L) + cldlsp)/(nswcld+1)
                1084        cldras_sw(i,j,L) = (nswcld*cldras_sw(i,j,L) + cldras)/(nswcld+1)
                1085 
                1086 C Compute Instantaneous Low-Mid-High Maximum Overlap Cloud Fractions
                1087 C ----------------------------------------------------------------------
1662f365b2 Andr*1088 
a7a27977ab Jean*1089        if( L.lt.midlevel ) cldhi (i,j) = max( cldnow,cldhi (i,j) )
1662f365b2 Andr*1090        if( L.ge.midlevel  .and.
a7a27977ab Jean*1091      .      L.lt.lowlevel ) cldmid(i,j) = max( cldnow,cldmid(i,j) )
                1092        if( L.ge.lowlevel ) cldlow(i,j) = max( cldnow,cldlow(i,j) )
1662f365b2 Andr*1093 
a7a27977ab Jean*1094 C Compute Cloud-Top Temperature and Pressure
                1095 C ------------------------------------------
                1096        cldtmp(i,j) =  cldtmp(i,j) + cldnow*pkzgather(i,L)
                1097      .             *  ( tz(i,j,L) + dtmoist(i,j,L)*tmstp/pz(i,j) )
                1098        cldprs(i,j) =  cldprs(i,j) + cldnow*plev
                1099        tmpimjm(i,j) = tmpimjm(i,j) + cldnow
1662f365b2 Andr*1100 
                1101       enddo
                1102       enddo
a7a27977ab Jean*1103       enddo
                1104 
                1105 C Compute Instantaneous Total 2-D Cloud Fraction
                1106 C ----------------------------------------------
1662f365b2 Andr*1107       do j = 1,jm
                1108       do i = 1,im
                1109       totcld(i,j) = 1.0 - (1.-cldhi (i,j))
                1110      .                  * (1.-cldmid(i,j))
                1111      .                  * (1.-cldlow(i,j))
                1112       enddo
                1113       enddo
a7a27977ab Jean*1114 
1662f365b2 Andr*1115 C **********************************************************************
                1116 C ***       Fill Cloud Top Pressure and Temperature Diagnostic       ***
                1117 C **********************************************************************
                1118 
8b150b4af9 Andr*1119 #ifdef ALLOW_DIAGNOSTICS
                1120       if(diagnostics_is_on('CLDTMP  ',myid) .and.
                1121      .                 diagnostics_is_on('CTTCNT  ',myid) ) then
                1122        do j=1,jm
a7a27977ab Jean*1123        do i=1,im
                1124         if( cldtmp(i,j).gt.0. ) then
8b150b4af9 Andr*1125          tmpdiag(i,j) = cldtmp(i,j)*totcld(i,j)/tmpimjm(i,j)
                1126          tmpdiag2(i,j) = totcld(i,j)
                1127         else
                1128          tmpdiag(i,j) = 0.
                1129          tmpdiag2(i,j) = 0.
                1130         endif
                1131        enddo
                1132        enddo
                1133        call diagnostics_fill(tmpdiag,'CLDTMP  ',0,1,3,bi,bj,myid)
                1134        call diagnostics_fill(tmpdiag2,'CTTCNT  ',0,1,3,bi,bj,myid)
1662f365b2 Andr*1135       endif
                1136 
8b150b4af9 Andr*1137       if(diagnostics_is_on('CLDPRS  ',myid) .and.
                1138      .                 diagnostics_is_on('CTPCNTC ',myid) ) then
                1139        do j=1,jm
a7a27977ab Jean*1140        do i=1,im
                1141         if( cldprs(i,j).gt.0. ) then
8b150b4af9 Andr*1142          tmpdiag(i,j) = cldprs(i,j)*totcld(i,j)/tmpimjm(i,j)
                1143          tmpdiag2(i,j) = totcld(i,j)
                1144         else
                1145          tmpdiag(i,j) = 0.
                1146          tmpdiag2(i,j) = 0.
                1147         endif
                1148        enddo
                1149        enddo
                1150        call diagnostics_fill(tmpdiag,'CLDPRS  ',0,1,3,bi,bj,myid)
                1151        call diagnostics_fill(tmpdiag2,'CTPCNT  ',0,1,3,bi,bj,myid)
1662f365b2 Andr*1152       endif
8b150b4af9 Andr*1153 
                1154 #endif
a7a27977ab Jean*1155 
1662f365b2 Andr*1156 C **********************************************************************
                1157 C ****                      INCREMENT COUNTERS                      ****
                1158 C **********************************************************************
a7a27977ab Jean*1159 
1662f365b2 Andr*1160        nlwlz    = nlwlz    + 1
                1161        nswlz    = nswlz    + 1
                1162 
                1163        nlwcld   = nlwcld   + 1
                1164        nswcld   = nswcld   + 1
                1165 
                1166       RETURN
                1167       END
471f7f004b Andr*1168       SUBROUTINE RAS( NN, LNG, LENC, K, NLTOP, nlayr, DT
689620ef36 Andr*1169      *,      UOI, ntracedim, ntracer, POI, QOI, PRS, PRJ, rnd, ncrnd
                1170      *,      RAINS, CLN, CLF, cldmas, detrain
                1171      *,      cp,grav,rkappa,alhl,rhfrac,rasmax )
1662f365b2 Andr*1172 C
                1173 C*********************************************************************
                1174 C********************* SUBROUTINE  RAS   *****************************
                1175 C********************** 16 MARCH   1988 ******************************
                1176 C*********************************************************************
                1177 C
7aae1e039f Andr*1178       implicit none
                1179 
36f0442dce Andr*1180 C Argument List
471f7f004b Andr*1181       integer nn,lng,lenc,k,nltop,nlayr
689620ef36 Andr*1182       integer ntracedim, ntracer
36f0442dce Andr*1183       integer ncrnd
a456aa407c Andr*1184       _RL dt
689620ef36 Andr*1185       _RL UOI(lng,nlayr,ntracedim),   POI(lng,K)
471f7f004b Andr*1186       _RL QOI(lng,K), PRS(lng,K+1), PRJ(lng,K+1)
a456aa407c Andr*1187       _RL rnd(ncrnd)
471f7f004b Andr*1188       _RL RAINS(lng,K), CLN(lng,K), CLF(lng,K)
                1189       _RL cldmas(lng,K), detrain(lng,K)
                1190       _RL cp,grav,rkappa,alhl,rhfrac(lng),rasmax
36f0442dce Andr*1191 
                1192 C Local Variables
471f7f004b Andr*1193       _RL TCU(lng,K), QCU(lng,K)
689620ef36 Andr*1194       _RL ucu(lng,K,ntracedim)
471f7f004b Andr*1195       _RL ALF(lng,K), BET(lng,K), GAM(lng,K)
                1196      *,         ETA(lng,K), HOI(lng,K)
                1197      *,         PRH(lng,K), PRI(lng,K)
                1198       _RL HST(lng,K), QOL(lng,K), GMH(lng,K)
a7a27977ab Jean*1199 
471f7f004b Andr*1200       _RL TX1(lng), TX2(lng), TX3(lng), TX4(lng), TX5(lng)
                1201      *,         TX6(lng), TX7(lng), TX8(lng), TX9(lng)
689620ef36 Andr*1202      *,         TX11(lng), TX12(lng), TX13(lng), TX14(lng,ntracedim)
471f7f004b Andr*1203      *,         TX15(lng)
a7a27977ab Jean*1204      *,         WFN(lng)
471f7f004b Andr*1205       integer IA1(lng), IA2(lng), IA3(lng)
                1206       _RL cloudn(lng), pcu(lng)
1662f365b2 Andr*1207 
7aae1e039f Andr*1208       integer krmin,icm
a456aa407c Andr*1209       _RL rknob, cmb2pa
7aae1e039f Andr*1210       PARAMETER (KRMIN=01)
                1211       PARAMETER (ICM=1000)
                1212       PARAMETER (CMB2PA=100.0)
                1213       PARAMETER (rknob = 10.)
36f0442dce Andr*1214 
                1215       integer IC(ICM),   IRND(icm)
471f7f004b Andr*1216       _RL cmass(lng,K)
36f0442dce Andr*1217       LOGICAL SETRAS
c74d45a9d3 Andr*1218       integer ifound
                1219       _RL temp
                1220       _RL thbef(lng,K)
36f0442dce Andr*1221 
                1222       integer i,L,nc,ib,nt
7aae1e039f Andr*1223       integer km1,kp1,kprv,kcr,kfx,ncmx
a456aa407c Andr*1224       _RL p00, crtmsf, frac, rasblf
7aae1e039f Andr*1225 
                1226       do L = 1,k
                1227       do I = 1,LENC
                1228        rains(i,l) = 0.
                1229       enddo
                1230       enddo
a7a27977ab Jean*1231 
1662f365b2 Andr*1232       p00 = 1000.
                1233       crtmsf = 0.
                1234 
                1235 C The numerator here is the fraction of the subcloud layer mass flux
                1236 C      allowed to entrain into the cloud
                1237 
                1238 CCC   FRAC = 1./dt
c74d45a9d3 Andr*1239 CCC   FRAC = 0.5/dt
1662f365b2 Andr*1240       FRAC = 0.5/dt
                1241 
                1242       KM1    = K  - 1
                1243       KP1    = K  + 1
                1244 C we want the ras adjustment time scale to be one hour (indep of dt)
f0525ae611 Andr*1245 C     RASBLF = 1./3600.
                1246 C we want the ras adjustment time scale to be one half hour (indep of dt)
                1247       RASBLF = 1./1800.
1662f365b2 Andr*1248 C
                1249        KPRV  = KM1
                1250 C Removed KRMAX parameter
a7a27977ab Jean*1251        KCR   = MIN(KM1,nlayr-2)
c74d45a9d3 Andr*1252        KFX   = KM1 - KCR
1662f365b2 Andr*1253        NCMX  = KFX + NCRND
                1254 C
                1255        IF (KFX .GT. 0) THEN
                1256         DO NC=1,KFX
                1257          IC(NC) = K - NC
                1258         ENDDO
                1259        ENDIF
                1260 C
                1261       IF (NCRND .GT. 0) THEN
                1262        DO I=1,ncrnd
                1263         IRND(I) = (RND(I)-0.0005)*(KCR-KRMIN+1)
                1264         IRND(I) = IRND(I) + KRMIN
                1265        ENDDO
                1266 C
                1267        DO NC=1,NCRND
                1268         IC(KFX+NC) = IRND(NC)
                1269        ENDDO
                1270       ENDIF
                1271 C
                1272       DO 100 NC=1,NCMX
                1273 C
                1274        IF (NC .EQ. 1 ) THEN
                1275         SETRAS = .TRUE.
                1276        ELSE
                1277         SETRAS = .FALSE.
                1278        ENDIF
                1279        IB = IC(NC)
                1280 
a7a27977ab Jean*1281 C Initialize Cloud Fraction Array
                1282 C -------------------------------
1662f365b2 Andr*1283       do i = 1,lenc
                1284       cloudn(i) = 0.0
                1285       enddo
                1286 
471f7f004b Andr*1287        CALL CLOUD(nn,lng, LENC, K, NLTOP, nlayr, IB, RASBLF,SETRAS,FRAC
1662f365b2 Andr*1288      *,           CP,  ALHL, RKAPPA, GRAV, P00, CRTMSF
689620ef36 Andr*1289      *,           POI, QOI, UOI, ntracedim, Ntracer, PRS, PRJ
1662f365b2 Andr*1290      *,           PCU, CLOUDN, TCU, QCU, UCU, CMASS
                1291      *,           ALF, BET, GAM, PRH, PRI, HOI, ETA
                1292      *,           HST, QOL, GMH
                1293      *,           TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, TX9
                1294      *,           WFN, TX11, TX12, TX13, TX14, TX15
                1295      *,           IA1,IA2,IA3,rhfrac)
                1296 
                1297 C Compute fraction of grid box into which rain re-evap occurs (clf)
a7a27977ab Jean*1298 C -----------------------------------------------------------------
1662f365b2 Andr*1299       do i = 1,lenc
                1300 
a7a27977ab Jean*1301 C mass in detrainment layer
                1302 C -------------------------
1662f365b2 Andr*1303        tx1(i) = cmb2pa *  (prs(i,ib+1) - prs(i,ib))/(grav*dt)
                1304 
a7a27977ab Jean*1305 C ratio of detraining cloud mass to mass in detrainment layer
                1306 C -----------------------------------------------------------
1662f365b2 Andr*1307        tx1(i) = rhfrac(i)*rknob * cmass(i,ib) / tx1(i)
                1308        if(cmass(i,K).gt.0.) clf(i,ib) = clf(i,ib) + tx1(i)
                1309        if( clf(i,ib).gt.1.) clf(i,ib) = 1.
                1310       enddo
a7a27977ab Jean*1311 
                1312 C Compute Total Cloud Mass Flux
                1313 C *****************************
1662f365b2 Andr*1314       do L=ib,k
                1315       do i=1,lenc
                1316        cmass(i,L) = rhfrac(i)*cmass(i,L) * dt
                1317       enddo
                1318       enddo
a7a27977ab Jean*1319 
1662f365b2 Andr*1320       do L=ib,k
                1321       do i=1,lenc
                1322        cldmas(i,L) = cldmas(i,L) + cmass(i,L)
                1323       enddo
                1324       enddo
                1325 
                1326       do i=1,lenc
                1327       detrain(i,ib) = detrain(i,ib) + cmass(i,ib)
                1328       enddo
                1329 
                1330       DO L=IB,K
                1331        DO I=1,LENC
c74d45a9d3 Andr*1332         thbef(I,L) = POI(I,L)
1662f365b2 Andr*1333         POI(I,L) = POI(I,L) + TCU(I,L) * DT * rhfrac(i)
                1334         QOI(I,L) = QOI(I,L) + QCU(I,L) * DT * rhfrac(i)
                1335        ENDDO
                1336       ENDDO
689620ef36 Andr*1337       DO NT=1,Ntracer
                1338       DO L=IB,K
                1339        DO I=1,LENC
                1340         UOI(I,L+nltop-1,NT)=UOI(I,L+nltop-1,NT)+UCU(I,L,NT)*DT*rhfrac(i)
                1341        ENDDO
                1342       ENDDO
                1343       ENDDO
1662f365b2 Andr*1344       DO I=1,LENC
                1345        rains(I,ib) = rains(I,ib) + PCU(I)*dt * rhfrac(i)
                1346       ENDDO
                1347 
c74d45a9d3 Andr*1348       do i = 1,lenc
                1349        ifound = 0
                1350        do L = 1,k
                1351         if(tcu(i,L).ne.0.)ifound = ifound + 1
                1352        enddo
                1353        if(ifound.ne.0) then
                1354 c       print *,i,' made a cloud detraining at ',ib
                1355         do L = 1,k
                1356          temp = TCU(I,L) * DT * rhfrac(i)
                1357 c        write(6,122)L,thbef(i,L),poi(i,L),temp
                1358         enddo
                1359        endif
                1360       enddo
                1361 
1662f365b2 Andr*1362   100 CONTINUE
c74d45a9d3 Andr*1363 
                1364  122  format(' ',i3,' TH B ',e10.3,' TH A ',e10.3,' DTH ',e10.3)
a7a27977ab Jean*1365 
                1366 C Fill Convective Cloud Fractions based on 3-D Rain Amounts
                1367 C ---------------------------------------------------------
1662f365b2 Andr*1368       do L=k-1,1,-1
                1369       do i=1,lenc
                1370          tx1(i)   = 100*(prs(i,L+1)-prs(i,L))/grav
                1371          cln(i,L) = min(1600*rains(i,L)/tx1(i),rasmax )
                1372       enddo
                1373       enddo
                1374 
                1375       RETURN
                1376       END
                1377       subroutine rndcloud (iras,nrnd,rnd,myid)
                1378       implicit none
                1379       integer n,iras,nrnd,myid
a456aa407c Andr*1380       _RL random_numbx
23d5b3245c Jean*1381 c     _RL rnd(nrnd)
                1382       _RL rnd(*)
1662f365b2 Andr*1383       integer irm
                1384       parameter (irm = 1000)
a456aa407c Andr*1385       _RL random(irm)
23d5b3245c Jean*1386       integer i,mcheck,iseed,indx
1662f365b2 Andr*1387       logical first
                1388       data    first /.true./
                1389       integer iras0
                1390       data    iras0 /0/
                1391       save random, iras0
a7a27977ab Jean*1392 
23d5b3245c Jean*1393       if(nrnd.eq.0)then
1662f365b2 Andr*1394        do i = 1,nrnd
                1395         rnd(i) = 0
                1396        enddo
9524fe64b5 Andr*1397        if(first .and. myid.eq.1) print *,' NO RANDOM CLOUDS IN RAS '
1662f365b2 Andr*1398        go to 100
                1399       endif
                1400 
                1401       mcheck = mod(iras-1,irm/nrnd)
                1402 
e0beae57f7 Jean*1403 c     print *,' RNDCLOUD: first ',first,' iras ',iras,' iras0 ',iras0
                1404 c     print *,' RNDCLOUD: irm,nrnd,mcheck=',irm,nrnd,mcheck
23d5b3245c Jean*1405 
                1406       if ( iras.eq.iras0 ) then
                1407 C-    Not the 1rst tile: we are all set (already done for the 1rst tile):
a7a27977ab Jean*1408 C -----------------------------------------------------------------------
23d5b3245c Jean*1409           indx = (iras-1)*nrnd
                1410 
a7a27977ab Jean*1411 C First Time In From a Continuing RESTART (IRAS.GT.1) or Reading a New RESTART
                1412 C   -- or --
                1413 C Multiple Time In But have Used Up all 1000 numbers (MCHECK.EQ.0)
                1414 C ----------------------------------------------------------------------------
23d5b3245c Jean*1415       elseif ( first.and.(iras.gt.1) .or. mcheck.eq.0 ) then
                1416        iseed = (iras-1-mcheck)*nrnd
1662f365b2 Andr*1417        call random_seedx(iseed)
                1418        do i = 1,irm
191d60053a Andr*1419         random(i) = random_numbx(iseed)
1662f365b2 Andr*1420        enddo
471f7f004b Andr*1421        indx = (iras-1)*nrnd
1662f365b2 Andr*1422 
23d5b3245c Jean*1423        if( myid.eq.1 ) print *, 'Creating Rand Numb Array in RNDCLOUD'
                1424      &                        ,', iseed=', iseed
                1425        if( myid.eq.1 ) print *, 'IRAS: ',iras,'  IRAS0: ',iras0,
                1426      &    ' indx: ', mod(indx,irm)
1662f365b2 Andr*1427 
a7a27977ab Jean*1428 C Multiple Time In But have NOT Used Up all 1000 numbers (MCHECK.NE.0)
                1429 C --------------------------------------------------------------------
1662f365b2 Andr*1430       else
471f7f004b Andr*1431           indx = (iras-1)*nrnd
1662f365b2 Andr*1432       endif
                1433 
471f7f004b Andr*1434           indx = mod(indx,irm)
23d5b3245c Jean*1435       if( indx+nrnd.gt.irm ) then
                1436 c       if( myid.eq.1 .AND. iras.ne.iras0 ) print *,
                1437 c    &   'reach end of Rand Numb Array in RNDCLOUD',indx,irm-nrnd
                1438         indx=irm-nrnd
                1439       endif
                1440 
1662f365b2 Andr*1441       do n = 1,nrnd
471f7f004b Andr*1442        rnd(n) = random(indx+n)
1662f365b2 Andr*1443       enddo
a7a27977ab Jean*1444 
1662f365b2 Andr*1445  100  continue
                1446       first = .false.
                1447       iras0 = iras
23d5b3245c Jean*1448 
1662f365b2 Andr*1449       return
                1450       end
191d60053a Andr*1451       function random_numbx(iseed)
1662f365b2 Andr*1452       implicit none
191d60053a Andr*1453       integer iseed
                1454       real *8 seed,port_rand
a456aa407c Andr*1455       _RL random_numbx
06d4643e1f Jean*1456 #ifdef FIZHI_CRAY
a456aa407c Andr*1457       _RL ranf
1662f365b2 Andr*1458       random_numbx = ranf()
a7a27977ab Jean*1459 #else
06d4643e1f Jean*1460 #ifdef FIZHI_SGI
a456aa407c Andr*1461       _RL rand
1662f365b2 Andr*1462       random_numbx = rand()
e0beae57f7 Jean*1463 #else
c25f76287c Jean*1464       seed = -1.d0
191d60053a Andr*1465       random_numbx = port_rand(seed)
                1466 #endif
e0beae57f7 Jean*1467 #endif
1662f365b2 Andr*1468       return
                1469       end
                1470       subroutine random_seedx (iseed)
                1471       implicit none
                1472       integer  iseed
4ec31d4dbd Andr*1473       real *8 port_rand
06d4643e1f Jean*1474 #ifdef FIZHI_CRAY
1662f365b2 Andr*1475       call ranset (iseed)
c25f76287c Jean*1476 #else
06d4643e1f Jean*1477 #ifdef FIZHI_SGI
1662f365b2 Andr*1478       integer*4   seed
                1479                   seed = iseed
                1480       call srand (seed)
c25f76287c Jean*1481 #else
                1482       real*8 tmpRdN
                1483       real*8 seed
                1484       seed = iseed
                1485       tmpRdN = port_rand(seed)
                1486 #endif
1662f365b2 Andr*1487 #endif
                1488       return
                1489       end
471f7f004b Andr*1490       SUBROUTINE CLOUD(nn,lng, LENC, K, NLTOP, nlayr, IC, RASALF
1662f365b2 Andr*1491      *,                 SETRAS, FRAC
                1492      *,                 CP,  ALHL, RKAP, GRAV, P00, CRTMSF
689620ef36 Andr*1493      *,                 POI, QOI, UOI, ntracedim, Ntracer, PRS,  PRJ
1662f365b2 Andr*1494      *,                 PCU, CLN, TCU, QCU, UCU, CMASS
                1495      *,                 ALF, BET, GAM, PRH, PRI, HOL, ETA
                1496      *,                 HST, QOL, GMH
                1497      *,                 TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, ALM
                1498      *,                 WFN, AKM, QS1, CLF, UHT, WLQ
                1499      *,                 IA, I1, I2,rhfrac)
                1500 C
                1501 C*********************************************************************
                1502 C******************** Relaxed Arakawa-Schubert ***********************
                1503 C********************* Plug Compatible Version  **********************
                1504 C************************ SUBROUTINE CLOUD ***************************
                1505 C************************* 23 JULY    1992 ***************************
                1506 C*********************************************************************
                1507 C*********************************************************************
                1508 C*********************************************************************
                1509 C************************** Developed By *****************************
                1510 C**************************              *****************************
                1511 C************************ Shrinivas Moorthi **************************
                1512 C************************       and         **************************
                1513 C************************  Max J. Suarez *****************************
                1514 C************************                *****************************
                1515 C******************** Laboratory for Atmospheres *********************
                1516 C****************** NASA/GSFC, Greenbelt, MD 20771 *******************
                1517 C*********************************************************************
                1518 C*********************************************************************
                1519 C
a7a27977ab Jean*1520 C   The calculations of Moorthi and Suarez (1992, MWR) are
1662f365b2 Andr*1521 C   contained in the CLOUD routine.
                1522 C   It is probably advisable, at least initially, to treat CLOUD
                1523 C   as a black box that computes the single cloud adjustments. RAS,
                1524 C   on the other hand, can be tailored to each GCMs configuration
                1525 C   (ie, number and placement of levels, nature of boundary layer,
                1526 C    time step and frequency with which RAS is called).
                1527 C
                1528 C
                1529 C  Input:
                1530 C  ------
                1531 C
471f7f004b Andr*1532 C     lng     : The inner dimension of update and input arrays.
1662f365b2 Andr*1533 C
                1534 C     LENC    : The run: the number of soundings processes in a single call.
471f7f004b Andr*1535 C               RAS works on the first LENC of the lng soundings
1662f365b2 Andr*1536 C               passed. This allows working on pieces of the world
                1537 C               say for multitasking, without declaring temporary arrays
4e1c053948 Jean*1538 C               and copying the data to and from them. This is a F77
1662f365b2 Andr*1539 C               version. An F90 version would have to allow more
                1540 C               flexibility in the argument declarations.  Obviously
a7a27977ab Jean*1541 C               (LENC<=lng).
                1542 C
1662f365b2 Andr*1543 C     K       : Number of vertical layers (increasing downwards).
                1544 C               Need not be the same as the number of layers in the
                1545 C               GCM, since it is the outer dimension. The bottom layer
a7a27977ab Jean*1546 C               (K) is the subcloud layer.
1662f365b2 Andr*1547 C
                1548 C     IC      : Detrainment level to check for presence of convection
                1549 C
                1550 C     RASALF  : Relaxation parameter (< 1.) for present cloud-type
                1551 C
                1552 C     SETRAS  : Logical parameter to control re-calculation of
                1553 C               saturation specific humidity and mid level P**kappa
                1554 C
                1555 C     FRAC    : Fraction of the PBL (layer K) mass allowed to be used
                1556 C               by a cloud-type in time DT
                1557 C
                1558 C     CP      : Specific heat at constant pressure
                1559 C
                1560 C     ALHL    : Latent Heat of condensation
                1561 C
                1562 C     RKAP    : R/Cp, where R is the gas constant
                1563 C
                1564 C     GRAV    : Acceleration due to gravity
                1565 C
                1566 C     P00     : A reference pressure in hPa, useually 1000 hPa
                1567 C
a7a27977ab Jean*1568 C     CRTMSF  : Critical value of mass flux above which cloudiness at
1662f365b2 Andr*1569 C               the detrainment layer of that cloud-type is assumed.
                1570 C               Affects only cloudiness calculation.
                1571 C
471f7f004b Andr*1572 C     POI     : 2D array of dimension (lng,K) containing potential
1662f365b2 Andr*1573 C               temperature. Updated but not initialized by RAS.
                1574 C
471f7f004b Andr*1575 C     QOI     : 2D array of dimension (lng,K) containing specific
1662f365b2 Andr*1576 C               humidity. Updated but not initialized by RAS.
                1577 C
471f7f004b Andr*1578 C     UOI     : 3D array of dimension (lng,K,NTRACER) containing tracers
1662f365b2 Andr*1579 C               Updated but not initialized by RAS.
                1580 C
471f7f004b Andr*1581 C     PRS     : 2D array of dimension (lng,K+1) containing pressure
a7a27977ab Jean*1582 C               in hPa at the interfaces of K-layers from top of the
1662f365b2 Andr*1583 C               atmosphere to the bottom. Not modified.
                1584 C
471f7f004b Andr*1585 C     PRJ     : 2D array of dimension (lng,K+1) containing (PRS/P00) **
1662f365b2 Andr*1586 C               RKAP.  i.e. Exner function at layer edges. Not modified.
                1587 C
471f7f004b Andr*1588 C     rhfrac  : 1D array of dimension (lng) containing a rel.hum. scaling
1662f365b2 Andr*1589 C               fraction. Not modified.
                1590 C
                1591 C  Output:
                1592 C  -------
                1593 C
471f7f004b Andr*1594 C     PCU     : 1D array of length lng containing accumulated
1662f365b2 Andr*1595 C               precipitation in mm/sec.
                1596 C
471f7f004b Andr*1597 C     CLN     : 2D array of dimension (lng,K) containing cloudiness
1662f365b2 Andr*1598 C               Note:  CLN is bumped but NOT initialized
                1599 C
a7a27977ab Jean*1600 C     TCU     : 2D array of dimension (lng,K) containing accumulated
                1601 C               convective heating (K/sec).
1662f365b2 Andr*1602 C
a7a27977ab Jean*1603 C     QCU     : 2D array of dimension (lng,K) containing accumulated
1662f365b2 Andr*1604 C               convective drying (kg/kg/sec).
                1605 C
471f7f004b Andr*1606 C     CMASS   : 2D array of dimension (lng,K) containing the
1662f365b2 Andr*1607 C               cloud mass flux (kg/sec). Filled from cloud top
                1608 C               to base.
                1609 C
                1610 C  Temporaries:
                1611 C
                1612 C      ALF, BET, GAM, ETA, PRH, PRI, HOI, HST, QOL, GMH are temporary
a7a27977ab Jean*1613 C               2D real arrays of dimension of at least (LENC,K) where LENC is
1662f365b2 Andr*1614 C               the horizontal dimension over which convection is invoked.
                1615 C
                1616 C
                1617 C    TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, TX9, AKM, QS1, CLF, UHT
                1618 C    VHT, WLQ  WFN are temporary real arrays of length at least LENC
                1619 C
a7a27977ab Jean*1620 C    IA, I1, and I2 are temporary integer arrays of length LENC
1662f365b2 Andr*1621 C
                1622 C
                1623 C************************************************************************
36f0442dce Andr*1624       implicit none
                1625 C Argument List declarations
689620ef36 Andr*1626       integer nn,lng,LENC,K,NLTOP,nlayr,ic,ntracedim, ntracer
a456aa407c Andr*1627       _RL rasalf
36f0442dce Andr*1628       LOGICAL SETRAS
a456aa407c Andr*1629       _RL frac, cp,  alhl, rkap, grav, p00, crtmsf
471f7f004b Andr*1630       _RL POI(lng,K),QOI(lng,K),PRS(lng,K+1),PRJ(lng,K+1)
689620ef36 Andr*1631       _RL uoi(lng,nlayr,ntracedim)
471f7f004b Andr*1632       _RL PCU(LENC), CLN(lng)
f266a50285 Andr*1633       _RL TCU(lng,K),QCU(lng,K),ucu(lng,k,ntracedim),CMASS(lng,K)
471f7f004b Andr*1634       _RL ALF(lng,K), BET(lng,K),  GAM(lng,K), PRH(lng,K), PRI(lng,K)
a456aa407c Andr*1635       _RL HOL(LENC,K), ETA(LENC,K), HST(LENC,K), QOL(LENC,K)
                1636       _RL GMH(LENC,K)
                1637       _RL TX1(LENC), TX2(LENC), TX3(LENC), TX4(LENC)
                1638       _RL TX5(LENC), TX6(LENC), TX7(LENC), TX8(LENC)
                1639       _RL ALM(LENC), WFN(LENC), AKM(LENC), QS1(LENC)
                1640       _RL WLQ(LENC), CLF(LENC)
689620ef36 Andr*1641       _RL uht(lng,ntracedim)
36f0442dce Andr*1642       integer IA(LENC), I1(LENC),I2(LENC)
471f7f004b Andr*1643       _RL      rhfrac(lng)
1662f365b2 Andr*1644 
36f0442dce Andr*1645 C Local Variables
a456aa407c Andr*1646       _RL daylen,half,one,zero,cmb2pa,rhmax
1662f365b2 Andr*1647       PARAMETER (DAYLEN=86400.0,  HALF=0.5,  ONE=1.0, ZERO=0.0)
                1648       PARAMETER (CMB2PA=100.0)
                1649       PARAMETER (RHMAX=0.9999)
a456aa407c Andr*1650       _RL rkapp1,onebcp,albcp,onebg,cpbg,twobal
1662f365b2 Andr*1651 C
36f0442dce Andr*1652       integer nt,km1,ic1,i,L,len1,len2,isav,len11,ii
471f7f004b Andr*1653       integer lena,lena1,lenb
                1654       _RL tem,tem1
a7a27977ab Jean*1655 
                1656 C Explicit Inline Directives
                1657 C --------------------------
06d4643e1f Jean*1658 #ifdef FIZHI_CRAY
4e1c053948 Jean*1659 #ifdef FIZHI_F77_COMPIL
1662f365b2 Andr*1660 cfpp$ expand (qsat)
                1661 #endif
                1662 #endif
                1663 
                1664       RKAPP1 = 1.0  + RKAP
                1665       ONEBCP = 1.0  / CP
                1666       ALBCP  = ALHL * ONEBCP
                1667       ONEBG  = 1.0  / GRAV
                1668       CPBG   = CP   * ONEBG
                1669       TWOBAL = 2.0 / ALHL
689620ef36 Andr*1670 
1662f365b2 Andr*1671 C
                1672       KM1 = K  - 1
                1673       IC1 = IC + 1
                1674 C
36f0442dce Andr*1675 C      SETTING ALF, BET, GAM, PRH, AND PRI : DONE ONLY WHEN SETRAS=.T.
1662f365b2 Andr*1676 C
                1677 
                1678       IF (SETRAS) THEN
                1679 
                1680          DO 2050 L=1,K
                1681          DO 2030 I=1,LENC
                1682           PRH(I,L) = (PRJ(I,L+1)*PRS(I,L+1) - PRJ(I,L)*PRS(I,L))
                1683      *                            / ((PRS(I,L+1)-PRS(I,L)) * RKAPP1)
                1684  2030    CONTINUE
                1685  2050    CONTINUE
                1686 
                1687          DO 2070 L=1,K
                1688          DO 2060 I=1,LENC
                1689           TX5(I) = POI(I,L) * PRH(I,L)
                1690           TX1(I) = (PRS(I,L) + PRS(I,L+1)) * 0.5
                1691           TX3(I) = TX5(I)
                1692           CALL QSAT(TX3(I), TX1(I), TX2(I), TX4(I), .TRUE.)
                1693           ALF(I,L) = TX2(I) - TX4(I) * TX5(I)
                1694           BET(I,L) = TX4(I) * PRH(I,L)
                1695           GAM(I,L) = 1.0 / ((1.0 + TX4(I)*ALBCP) * PRH(I,L))
                1696           PRI(I,L) = (CP/CMB2PA) / (PRS(I,L+1) - PRS(I,L))
                1697  2060    CONTINUE
                1698  2070    CONTINUE
                1699 
                1700       ENDIF
                1701 C
                1702 C
                1703       DO 10 L=1,K
471f7f004b Andr*1704       DO 10 I=1,lng
1662f365b2 Andr*1705       TCU(I,L) = 0.0
                1706       QCU(I,L) = 0.0
                1707       CMASS(I,L) = 0.0
                1708    10 CONTINUE
                1709 
689620ef36 Andr*1710       do nt = 1,ntracer
                1711       do L=1,K
                1712       do I=1,LENC
                1713       ucu(I,L,nt) = 0.0
                1714       enddo
                1715       enddo
                1716       enddo
1662f365b2 Andr*1717 C
                1718       DO 30 I=1,LENC
                1719       TX1(I)   = PRJ(I,K+1) * POI(I,K)
                1720       QS1(I)   = ALF(I,K) + BET(I,K)*POI(I,K)
                1721       QOL(I,K) = MIN(QS1(I)*RHMAX,QOI(I,K))
                1722 
                1723       HOL(I,K) = TX1(I)*CP + QOL(I,K)*ALHL
                1724       ETA(I,K) = ZERO
                1725       TX2(I)   = (PRJ(I,K+1) - PRJ(I,K)) * POI(I,K) * CP
                1726    30 CONTINUE
                1727 C
                1728       IF (IC .LT. KM1) THEN
                1729          DO 3703 L=KM1,IC1,-1
                1730          DO 50 I=1,LENC
                1731          QS1(I)   = ALF(I,L) + BET(I,L)*POI(I,L)
                1732          QOL(I,L) = MIN(QS1(I)*RHMAX,QOI(I,L))
                1733 C
a7a27977ab Jean*1734          TEM1 = TX2(I) + PRJ(I,L+1) * POI(I,L) * CP
1662f365b2 Andr*1735          HOL(I,L) = TEM1 + QOL(I,L )* ALHL
                1736          HST(I,L) = TEM1 + QS1(I)   * ALHL
                1737 
                1738          TX1(I)   = (PRJ(I,L+1) - PRJ(I,L)) * POI(I,L)
                1739          ETA(I,L) = ETA(I,L+1) + TX1(I)*CPBG
                1740          TX2(I)   = TX2(I)     + TX1(I)*CP
                1741    50    CONTINUE
                1742 C
                1743  3703    CONTINUE
                1744       ENDIF
                1745 
                1746       DO 70 I=1,LENC
                1747       HOL(I,IC) = TX2(I)
                1748       QS1(I)    = ALF(I,IC) + BET(I,IC)*POI(I,IC)
a7a27977ab Jean*1749       QOL(I,IC) = MIN(QS1(I)*RHMAX,QOI(I,IC))
                1750 C
                1751       TEM1      = TX2(I) + PRJ(I,IC1) * POI(I,IC) * CP
1662f365b2 Andr*1752       HOL(I,IC) = TEM1 + QOL(I,IC) * ALHL
                1753       HST(I,IC) = TEM1 + QS1(I)    * ALHL
                1754 C
                1755       TX3(I   ) = (PRJ(I,IC1) - PRH(I,IC)) * POI(I,IC)
                1756       ETA(I,IC) = ETA(I,IC1) + CPBG * TX3(I)
                1757    70 CONTINUE
                1758 C
                1759       DO 130 I=1,LENC
                1760       TX2(I) = HOL(I,K)  - HST(I,IC)
                1761       TX1(I) = ZERO
                1762 
                1763   130 CONTINUE
                1764 C
                1765 C     ENTRAINMENT PARAMETER
                1766 C
                1767       DO 160 L=IC,KM1
                1768       DO 160 I=1,LENC
                1769       TX1(I) = TX1(I) + (HST(I,IC) - HOL(I,L)) * (ETA(I,L) - ETA(I,L+1))
                1770   160 CONTINUE
                1771 C
                1772       LEN1 = 0
                1773       LEN2 = 0
                1774       ISAV = 0
                1775       DO 195 I=1,LENC
                1776       IF (TX1(I) .GT. ZERO .AND. TX2(I) .GT. ZERO
                1777      .                     .AND. rhfrac(i).ne.0.0 ) THEN
                1778          LEN1      = LEN1 + 1
                1779          IA(LEN1)  = I
                1780          ALM(LEN1) = TX2(I) / TX1(I)
                1781       ENDIF
                1782   195 CONTINUE
                1783 C
                1784       LEN2 = LEN1
                1785       if (IC1 .lt. K) then
                1786          DO 196 I=1,LENC
                1787          IF (TX2(I) .LE. 0.0 .AND. (HOL(I,K) .GT. HST(I,IC1))
                1788      .                       .AND. rhfrac(i).ne.0.0 ) THEN
                1789             LEN2      = LEN2 + 1
                1790             IA(LEN2)  = I
                1791             ALM(LEN2) = 0.0
                1792          ENDIF
                1793   196    CONTINUE
                1794       endif
                1795 C
a7a27977ab Jean*1796       IF (LEN2 .EQ. 0) THEN
                1797 c        DO 5010 I=1,LENC*K
                1798 c        HST(I,1) = 0.0
                1799 c        QOL(I,1) = 0.0
                1800 c5010    CONTINUE
                1801          DO L = 1,K
                1802           DO I = 1,LENC
                1803            HST(I,L) = 0.0
                1804            QOL(I,L) = 0.0
                1805           ENDDO
                1806          ENDDO
1662f365b2 Andr*1807          DO 5020 I=1,LENC
                1808          PCU(I) = 0.0
                1809  5020    CONTINUE
                1810          RETURN
                1811       ENDIF
                1812       LEN11 = LEN1 + 1
                1813 C
                1814 C     NORMALIZED MASSFLUX
                1815 C
                1816       DO 250 I=1,LEN2
                1817       ETA(I,K) = 1.0
                1818       II       = IA(I)
                1819       TX2(I)   = 0.5 * (PRS(II,IC) + PRS(II,IC1))
                1820       TX4(I)   = PRS(II,K)
                1821   250 CONTINUE
                1822 C
                1823       DO 252 I=LEN11,LEN2
                1824       WFN(I)   = 0.0
                1825       II       = IA(I)
                1826       IF (HST(II,IC1) .LT. HST(II,IC)) THEN
                1827          TX6(I) = (HST(II,IC1)-HOL(II,K))/(HST(II,IC1)-HST(II,IC))
                1828       ELSE
                1829          TX6(I) = 0.0
                1830       ENDIF
                1831       TX2(I) = 0.5 * (PRS(II,IC1)+PRS(II,IC1+1)) * (1.0-TX6(I))
                1832      *                             + TX2(I)      * TX6(I)
                1833   252 CONTINUE
                1834 C
                1835       CALL ACRITN(LEN2, TX2, TX4, TX3)
                1836 C
                1837       DO 260 L=KM1,IC,-1
                1838       DO 255 I=1,LEN2
                1839       TX1(I) = ETA(IA(I),L)
                1840   255 CONTINUE
                1841       DO 260 I=1,LEN2
                1842       ETA(I,L) = 1.0 + ALM(I) * TX1(I)
                1843   260 CONTINUE
                1844 C
                1845 C     CLOUD WORKFUNCTION
                1846 C
                1847       IF (LEN1 .GT. 0) THEN
                1848          DO 270 I=1,LEN1
                1849          II = IA(I)
                1850          WFN(I) = - GAM(II,IC) * (PRJ(II,IC1) - PRH(II,IC))
                1851      *                         *  HST(II,IC) * ETA(I,IC1)
                1852   270    CONTINUE
                1853       ENDIF
                1854 C
                1855       DO 290 I=1,LEN2
                1856       II = IA(I)
                1857       TX1(I) = HOL(II,K)
                1858   290 CONTINUE
                1859 C
                1860       IF (IC1 .LE. KM1) THEN
                1861 
                1862          DO 380 L=KM1,IC1,-1
                1863          DO 380 I=1,LEN2
                1864          II = IA(I)
                1865          TEM = TX1(I) + (ETA(I,L) - ETA(I,L+1)) * HOL(II,L)
                1866 C
                1867          PCU(I) = PRJ(II,L+1) - PRH(II,L)
                1868          TEM1   = ETA(I,L+1) * PCU(I)
                1869          TX1(I) = TX1(I)*PCU(I)
                1870 C
                1871          PCU(I) = PRH(II,L) - PRJ(II,L)
                1872          TEM1   = (TEM1 + ETA(I,L) * PCU(I)) * HST(II,L)
                1873          TX1(I) = TX1(I) + TEM*PCU(I)
                1874 C
                1875          WFN(I) = WFN(I) + (TX1(I) - TEM1) * GAM(II,L)
                1876          TX1(I) = TEM
                1877   380    CONTINUE
                1878       ENDIF
                1879 C
                1880       LENA = 0
                1881       IF (LEN1 .GT. 0) THEN
                1882          DO 512 I=1,LEN1
                1883          II = IA(I)
                1884          WFN(I) = WFN(I) + TX1(I) * GAM(II,IC)*(PRJ(II,IC1)-PRH(II,IC))
                1885      *                   - TX3(I)
                1886          IF (WFN(I) .GT. 0.0) THEN
                1887             LENA = LENA + 1
                1888             I1(LENA) = IA(I)
                1889             I2(LENA) = I
                1890             TX1(LENA) = WFN(I)
                1891             TX2(LENA) = QS1(IA(I))
                1892             TX6(LENA) = 1.0
                1893          ENDIF
                1894   512    CONTINUE
                1895       ENDIF
                1896       LENB = LENA
                1897       DO 515 I=LEN11,LEN2
                1898       WFN(I) = WFN(I) - TX3(I)
                1899       IF (WFN(I) .GT. 0.0 .AND. TX6(I) .GT. 0.0) THEN
                1900          LENB = LENB + 1
                1901          I1(LENB)  = IA(I)
                1902          I2(LENB)  = I
                1903          TX1(LENB) = WFN(I)
                1904          TX2(LENB) = QS1(IA(I))
                1905          TX4(LENB) = TX6(I)
                1906       ENDIF
                1907   515 CONTINUE
                1908 C
                1909       IF (LENB .LE. 0) THEN
877fb4cee5 Jean*1910 c        DO 5030 I=1,LENC*K
                1911 c        HST(I,1) = 0.0
                1912 c        QOL(I,1) = 0.0
                1913 c5030    CONTINUE
                1914          DO L = 1,K
                1915           DO I = 1,LENC
                1916            HST(I,L) = 0.0
                1917            QOL(I,L) = 0.0
                1918           ENDDO
                1919          ENDDO
1662f365b2 Andr*1920          DO 5040 I=1,LENC
                1921          PCU(I) = 0.0
                1922  5040    CONTINUE
                1923          RETURN
                1924       ENDIF
                1925 
                1926 C
                1927       DO 516 I=1,LENB
                1928       WFN(I) = TX1(I)
                1929       QS1(I) = TX2(I)
                1930   516 CONTINUE
                1931 C
                1932       DO 520 L=IC,K
                1933       DO 517 I=1,LENB
                1934       TX1(I) = ETA(I2(I),L)
                1935   517 CONTINUE
                1936       DO 520 I=1,LENB
                1937       ETA(I,L) = TX1(I)
                1938   520 CONTINUE
                1939 C
                1940       LENA1 = LENA + 1
                1941 C
                1942       DO 510 I=1,LENA
                1943       II = I1(I)
                1944       TX8(I) = HST(II,IC) - HOL(II,IC)
                1945   510 CONTINUE
                1946       DO 530 I=LENA1,LENB
                1947       II = I1(I)
                1948       TX6(I) = TX4(I)
                1949       TEM    = TX6(I) * (HOL(II,IC)-HOL(II,IC1)) + HOL(II,IC1)
                1950       TX8(I) = HOL(II,K) - TEM
                1951 
                1952       TEM1   = TX6(I) * (QOL(II,IC)-QOL(II,IC1)) + QOL(II,IC1)
                1953       TX5(I) = TEM    - TEM1 * ALHL
                1954       QS1(I) = TEM1   + TX8(I)*(ONE/ALHL)
                1955       TX3(I) = HOL(II,IC)
                1956   530 CONTINUE
                1957 C
                1958 C
                1959       DO 620 I=1,LENB
                1960       II = I1(I)
                1961       WLQ(I) = QOL(II,K) - QS1(I)     * ETA(I,IC)
                1962       TX7(I) = HOL(II,K)
                1963   620 CONTINUE
689620ef36 Andr*1964       DO NT=1,Ntracer
                1965       DO 621 I=1,LENB
                1966       II = I1(I)
                1967       UHT(I,NT) = UOI(II,K+nltop-1,NT)-UOI(II,IC+nltop-1,NT) * ETA(I,IC)
                1968   621 CONTINUE
                1969       ENDDO
1662f365b2 Andr*1970 C
                1971       DO 635 L=KM1,IC,-1
                1972       DO 630 I=1,LENB
                1973       II = I1(I)
                1974       TEM    = ETA(I,L) - ETA(I,L+1)
                1975       WLQ(I) = WLQ(I) + TEM * QOL(II,L)
                1976   630 CONTINUE
                1977   635 CONTINUE
689620ef36 Andr*1978       DO NT=1,Ntracer
                1979       DO L=KM1,IC,-1
                1980       DO I=1,LENB
                1981       II = I1(I)
                1982       TEM    = ETA(I,L) - ETA(I,L+1)
                1983       UHT(I,NT) = UHT(I,NT) + TEM *  UOI(II,L+nltop-1,NT)
                1984       ENDDO
                1985       ENDDO
                1986       ENDDO
1662f365b2 Andr*1987 C
                1988 C     CALCULATE GS AND PART OF AKM (THAT REQUIRES ETA)
                1989 C
                1990       DO 690 I=1,LENB
                1991       II = I1(I)
                1992 c     TX7(I)     = HOL(II,K)
                1993       TEM        = (POI(II,KM1) - POI(II,K)) / (PRH(II,K) - PRH(II,KM1))
                1994       HOL(I,K)   = TEM * (PRJ(II,K)-PRH(II,KM1))*PRH(II,K)*PRI(II,K)
                1995       HOL(I,KM1) = TEM * (PRH(II,K)-PRJ(II,K))*PRH(II,KM1)*PRI(II,KM1)
                1996       AKM(I)     = ZERO
                1997       TX2(I)     = 0.5 * (PRS(II,IC) + PRS(II,IC1))
                1998   690 CONTINUE
                1999 
                2000       IF (IC1 .LE. KM1) THEN
                2001          DO 750 L=KM1,IC1,-1
                2002          DO 750 I=1,LENB
                2003          II = I1(I)
                2004          TEM      = (POI(II,L-1) - POI(II,L)) * ETA(I,L)
                2005      *                            / (PRH(II,L) - PRH(II,L-1))
                2006 C
                2007          HOL(I,L)   = TEM * (PRJ(II,L)-PRH(II,L-1)) * PRH(II,L)
                2008      *                    *  PRI(II,L)  + HOL(I,L)
                2009          HOL(I,L-1) = TEM * (PRH(II,L)-PRJ(II,L)) * PRH(II,L-1)
                2010      *                                              * PRI(II,L-1)
                2011 C
                2012          AKM(I)   = AKM(I) - HOL(I,L)
                2013      *            * (ETA(I,L)   * (PRH(II,L)-PRJ(II,L)) +
                2014      *               ETA(I,L+1) * (PRJ(II,L+1)-PRH(II,L))) / PRH(II,L)
                2015   750    CONTINUE
                2016       ENDIF
                2017 C
                2018 C
                2019       CALL RNCL(LENB, TX2, TX1, CLF)
                2020 
                2021       DO 770 I=1,LENB
                2022       TX2(I) = (ONE - TX1(I)) * WLQ(I)
                2023       WLQ(I) = TX1(I) * WLQ(I)
                2024 C
                2025       TX1(I) = HOL(I,IC)
                2026   770 CONTINUE
                2027       DO 790 I=LENA1, LENB
                2028       II = I1(I)
                2029       TX1(I) = TX1(I) + (TX5(I)-TX3(I)+QOL(II,IC)*ALHL)*(PRI(II,IC)/CP)
                2030   790 CONTINUE
                2031 
                2032       DO 800 I=1,LENB
                2033       HOL(I,IC) = TX1(I) - TX2(I) * ALBCP * PRI(I1(I),IC)
                2034   800 CONTINUE
                2035 
                2036       IF (LENA .GT. 0) THEN
                2037          DO 810 I=1,LENA
                2038          II = I1(I)
a7a27977ab Jean*2039          AKM(I) = AKM(I) - ETA(I,IC1) * (PRJ(II,IC1) - PRH(II,IC))
1662f365b2 Andr*2040      *                                * TX1(I) / PRH(II,IC)
                2041   810    CONTINUE
                2042       ENDIF
a7a27977ab Jean*2043 C
1662f365b2 Andr*2044 C     CALCULATE GH
                2045 C
                2046       DO 830 I=1,LENB
                2047       II = I1(I)
                2048       TX3(I)   =  QOL(II,KM1) - QOL(II,K)
                2049       GMH(I,K) = HOL(I,K) + TX3(I) * PRI(II,K) * (ALBCP)
                2050 
a7a27977ab Jean*2051       AKM(I)   = AKM(I) + GAM(II,KM1)*(PRJ(II,K)-PRH(II,KM1))
1662f365b2 Andr*2052      *                               * GMH(I,K)
                2053       TX3(I)   =  zero
                2054   830 CONTINUE
                2055 C
                2056       IF (IC1 .LE. KM1) THEN
                2057          DO 840 L=KM1,IC1,-1
                2058          DO 840 I=1,LENB
                2059          II = I1(I)
                2060          TX2(I) = TX3(I)
                2061          TX3(I) = (QOL(II,L-1) - QOL(II,L)) * ETA(I,L)
                2062          TX2(I) = TX2(I) + TX3(I)
                2063 C
                2064          GMH(I,L) = HOL(I,L) + TX2(I)   * PRI(II,L) * (ALBCP*HALF)
                2065   840    CONTINUE
                2066 C
                2067 C
                2068       ENDIF
                2069       DO 850 I=LENA1,LENB
                2070       TX3(I) = TX3(I) + TWOBAL
                2071      *       * (TX7(I) - TX8(I) - TX5(I) - QOL(I1(I),IC)*ALHL)
                2072   850 CONTINUE
                2073       DO 860 I=1,LENB
                2074       GMH(I,IC) = TX1(I) + PRI(I1(I),IC) * ONEBCP
                2075      *          * (TX3(I)*(ALHL*HALF) + ETA(I,IC) * TX8(I))
                2076   860 CONTINUE
                2077 C
                2078 C     CALCULATE HC PART OF AKM
                2079 C
                2080       IF (IC1 .LE. KM1) THEN
                2081          DO 870 I=1,LENB
                2082          TX1(I) = GMH(I,K)
                2083   870    CONTINUE
                2084          DO 3725 L=KM1,IC1,-1
                2085          DO 880 I=1,LENB
                2086          II = I1(I)
                2087          TX1(I) = TX1(I) + (ETA(I,L) - ETA(I,L+1)) * GMH(I,L)
                2088          TX2(I) = GAM(II,L-1) * (PRJ(II,L) - PRH(II,L-1))
                2089   880    CONTINUE
                2090 C
                2091          IF (L .EQ. IC1) THEN
                2092             DO 890 I=LENA1,LENB
                2093             TX2(I) = ZERO
                2094   890       CONTINUE
                2095          ENDIF
                2096          DO 900 I=1,LENB
                2097          II = I1(I)
a7a27977ab Jean*2098          AKM(I) = AKM(I) + TX1(I) *
1662f365b2 Andr*2099      *          (TX2(I) + GAM(II,L)*(PRH(II,L)-PRJ(II,L)))
                2100   900    CONTINUE
                2101  3725    CONTINUE
                2102       ENDIF
                2103 C
                2104       DO 920 I=LENA1,LENB
                2105       II = I1(I)
                2106       TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1))
                2107      *       + 0.5*(PRS(II,IC+2) - PRS(II,IC)) * (ONE-TX6(I))
a7a27977ab Jean*2108 C
1662f365b2 Andr*2109       TX1(I) = PRS(II,IC1)
                2110       TX5(I) = 0.5 * (PRS(II,IC1) + PRS(II,IC+2))
                2111 C
                2112       IF ((TX2(I) .GE. TX1(I)) .AND. (TX2(I) .LT. TX5(I))) THEN
                2113          TX6(I)     = ONE - (TX2(I) - TX1(I)) / (TX5(I) - TX1(I))
                2114 C
                2115          TEM        = PRI(II,IC1) / PRI(II,IC)
                2116          HOL(I,IC1) = HOL(I,IC1) + HOL(I,IC) * TEM
                2117          HOL(I,IC)  = ZERO
                2118 C
                2119          GMH(I,IC1) = GMH(I,IC1) + GMH(I,IC) * TEM
                2120          GMH(I,IC)  = ZERO
                2121       ELSEIF (TX2(I) .LT. TX1(I)) THEN
                2122          TX6(I) = 1.0
                2123       ELSE
                2124          TX6(I) = 0.0
                2125       ENDIF
                2126   920 CONTINUE
                2127 C
                2128 C
                2129       DO I=1,LENC
                2130       PCU(I) = 0.0
                2131       ENDDO
                2132 
                2133       DO 970 I=1,LENB
                2134       II = I1(I)
                2135       IF (AKM(I) .LT. ZERO .AND. WLQ(I) .GE. 0.0) THEN
                2136          WFN(I) = - TX6(I) * WFN(I) * RASALF / AKM(I)
                2137       ELSE
                2138          WFN(I) = ZERO
                2139       ENDIF
                2140       TEM       = (PRS(II,K+1)-PRS(II,K))*(CMB2PA*FRAC)
                2141       WFN(I)    = MIN(WFN(I), TEM)
                2142 C
                2143 C     compute cloud amount
                2144 C
                2145 CC    TX1(I) = CLN(II)
                2146 CC    IF (WFN(I) .GT. CRTMSF)  TX1(I) = TX1(I) + CLF(I)
                2147 CC    IF (TX1(I) .GT. ONE)  TX1(I) = ONE
                2148 C
                2149 C     PRECIPITATION
                2150 C
                2151       PCU(II) =  WLQ(I) * WFN(I) * ONEBG
                2152 C
                2153 C     CUMULUS FRICTION AT THE BOTTOM LAYER
                2154 C
                2155       TX4(I)   = WFN(I) * (1.0/ALHL)
                2156       TX5(I)   = WFN(I) * ONEBCP
                2157   970 CONTINUE
                2158 C
                2159 C     compute cloud mass flux for diagnostic output
                2160 C
                2161       DO L = IC,K
                2162       DO I=1,LENB
                2163        II = I1(I)
                2164        if(L.lt.K)then
                2165        CMASS(II,L) =  ETA(I,L+1) * WFN(I) * ONEBG
                2166        else
                2167        CMASS(II,L) =  WFN(I) * ONEBG
                2168        endif
                2169       ENDDO
                2170       ENDDO
                2171 C
                2172 CC    DO 975 I=1,LENB
                2173 CC    II = I1(I)
                2174 CC    CLN(II) = TX1(I)
                2175 CC975 CONTINUE
                2176 C
                2177 C     THETA AND Q CHANGE DUE TO CLOUD TYPE IC
                2178 C
a7a27977ab Jean*2179 
1662f365b2 Andr*2180 c     TEMA = 0.0
                2181 c     TEMB = 0.0
                2182       DO 990 L=IC,K
                2183       DO 980 I=1,LENB
                2184       II = I1(I)
                2185       TEM      = (GMH(I,L) - HOL(I,L)) * TX4(I)
                2186       TEM1     =  HOL(I,L) * TX5(I)
                2187 C
                2188       TCU(II,L) = TEM1 / PRH(II,L)
                2189       QCU(II,L) = TEM
                2190   980 CONTINUE
                2191 
                2192 c     I = I1(IP1)
                2193 c
                2194 c     TEM  = (PRS(I,L+1)-PRS(I,L)) * (ONEBG*100.0)
                2195 c     TEMA = TEMA +  TCU(I,L) * PRH(I,L) * TEM * (CP/ALHL)
                2196 c     TEMB = TEMB +  QCU(I,L)            * TEM
                2197 C
                2198   990 CONTINUE
                2199 C
a7a27977ab Jean*2200 C Compute Tracer Tendencies
                2201 C -------------------------
689620ef36 Andr*2202       do nt = 1,ntracer
a7a27977ab Jean*2203 C
                2204 C Tracer Tendency at the Bottom Layer
                2205 C -----------------------------------
689620ef36 Andr*2206       DO 995 I=1,LENB
                2207       II = I1(I)
                2208       TEM    = half*TX5(I) * PRI(II,K)
                2209       TX1(I) = ( UOI(II,KM1+nltop-1,nt) - UOI(II,K+nltop-1,nt))
                2210       ucu(II,K,nt) = TEM * TX1(I)
                2211   995 CONTINUE
a7a27977ab Jean*2212 C
                2213 C Tracer Tendency at all other Levels
                2214 C -----------------------------------
689620ef36 Andr*2215       DO 1020 L=KM1,IC1,-1
                2216       DO 1010 I=1,LENB
                2217       II = I1(I)
                2218       TEM = half*TX5(I) * PRI(II,L)
                2219       TEM1   = TX1(I)
                2220       TX1(I) = (UOI(II,L-1+nltop-1,nt)-UOI(II,L+nltop-1,nt)) * ETA(I,L)
                2221       TX3(I) = (TX1(I) + TEM1) * TEM
                2222  1010 CONTINUE
                2223       DO 1020 I=1,LENB
                2224       II = I1(I)
                2225       ucu(II,L,nt) = TX3(I)
                2226  1020 CONTINUE
889be31914 Andr*2227 
689620ef36 Andr*2228       DO 1030 I=1,LENB
                2229       II = I1(I)
                2230       IF (TX6(I) .GE. 1.0) THEN
                2231          TEM    = half*TX5(I) * PRI(II,IC)
                2232       ELSE
                2233          TEM = 0.0
                2234       ENDIF
                2235       TX1(I) = (TX1(I) + UHT(I,nt) + UHT(I,nt)) * TEM
                2236  1030 CONTINUE
                2237       DO 1040 I=1,LENB
                2238       II = I1(I)
                2239       ucu(II,IC,nt) = TX1(I)
                2240  1040 CONTINUE
a7a27977ab Jean*2241 
689620ef36 Andr*2242       enddo
1662f365b2 Andr*2243 C
                2244 C     PENETRATIVE CONVECTION CALCULATION OVER
                2245 C
                2246 
                2247       RETURN
                2248       END
471f7f004b Andr*2249       SUBROUTINE RNCL(lng, PL, RNO, CLF)
1662f365b2 Andr*2250 C
                2251 C*********************************************************************
                2252 C********************** Relaxed Arakawa-Schubert *********************
                2253 C************************   SUBROUTINE  RNCL  ************************
                2254 C**************************** 23 July 1992 ***************************
                2255 C*********************************************************************
36f0442dce Andr*2256       implicit none
                2257 C Argument List declarations
471f7f004b Andr*2258       integer lng
                2259       _RL PL(lng),  RNO(lng), CLF(lng)
1662f365b2 Andr*2260 
36f0442dce Andr*2261 C Local Variables
a456aa407c Andr*2262       _RL p5,p8,pt8,pt2,pfac,p4,p6,p7,p9,cucld,cfac
1662f365b2 Andr*2263       PARAMETER (P5=500.0,  P8=800.0, PT8=0.8, PT2=0.2)
                2264       PARAMETER (PFAC=PT2/(P8-P5))
                2265       PARAMETER (P4=400.0,    P6=401.0)
                2266       PARAMETER (P7=700.0,    P9=900.0)
                2267       PARAMETER (CUCLD=0.5,CFAC=CUCLD/(P6-P4))
36f0442dce Andr*2268 
                2269       integer i
1662f365b2 Andr*2270 C
471f7f004b Andr*2271       DO 10 I=1,lng
1662f365b2 Andr*2272                            rno(i) = 1.0
e7070f537c Cons*2273 ccc   if( pl(i).le.400.0 ) rno(i) = max( 0.75 _d 0, 1.0-0.0025*
                2274 ccc  &                                 (400.0-pl(i)) )
1662f365b2 Andr*2275 
                2276 ccc   IF ( PL(I).GE.P7 .AND. PL(I).LE.P9 ) THEN
                2277 ccc     RNO(I) = ((P9-PL(I))/(P9-P7)) **2
                2278 ccc   ELSE IF (PL(I).GT.P9) THEN
                2279 ccc     RNO(I) = 0.
                2280 ccc   ENDIF
                2281 
                2282       CLF(I) = CUCLD
                2283 C
                2284 CARIESIF (PL(I) .GE. P5 .AND. PL(I) .LE. P8) THEN
                2285 CARIES    RNO(I) = (P8-PL(I))*PFAC + PT8
                2286 CARIESELSEIF (PL(I) .GT. P8 ) THEN
                2287 CARIES    RNO(I) = PT8
                2288 CARIESENDIF
                2289 CARIES
                2290       IF (PL(I) .GE. P4 .AND. PL(I) .LE. P6) THEN
                2291          CLF(I) = (P6-PL(I))*CFAC
                2292       ELSEIF (PL(I) .GT. P6 ) THEN
                2293          CLF(I) = 0.0
                2294       ENDIF
                2295    10 CONTINUE
                2296 C
                2297       RETURN
                2298       END
471f7f004b Andr*2299       SUBROUTINE ACRITN ( lng,PL,PLB,ACR )
a7a27977ab Jean*2300 
1662f365b2 Andr*2301 C*********************************************************************
                2302 C********************** Relaxed Arakawa-Schubert *********************
                2303 C************************** SUBROUTINE  ACRIT    *********************
                2304 C******************  modified August 28, 1996   L.Takacs  ************
                2305 C****                                                            *****
                2306 C****  Note:  Data obtained from January Mean After-Analysis     *****
                2307 C****         from 4x5 46-layer GEOS Assimilation                *****
                2308 C****                                                            *****
                2309 C*********************************************************************
36f0442dce Andr*2310       implicit none
                2311 C Argument List declarations
471f7f004b Andr*2312       integer lng
                2313       _RL PL(lng), PLB(lng), ACR(lng)
1662f365b2 Andr*2314 
36f0442dce Andr*2315 C Local variables
                2316       integer lma
1662f365b2 Andr*2317       parameter  (lma=18)
a456aa407c Andr*2318       _RL p(lma)
                2319       _RL a(lma)
36f0442dce Andr*2320       integer i,L
a456aa407c Andr*2321       _RL temp
1662f365b2 Andr*2322 
                2323       data p  / 93.81, 111.65, 133.46, 157.80, 186.51,
                2324      .         219.88, 257.40, 301.21, 352.49, 409.76,
                2325      .         471.59, 535.04, 603.33, 672.79, 741.12,
                2326      .         812.52, 875.31, 930.20/
                2327 
                2328       data a  / 3.35848, 3.13645, 2.48072, 2.08277, 1.53364,
                2329      .          1.01971,  .65846,  .45867,  .38687,  .31002,
                2330      .           .25574,  .20347,  .17254,  .15260,  .16756,
                2331      .           .09916,  .10360,  .05880/
                2332 
                2333       do L=1,lma-1
471f7f004b Andr*2334       do i=1,lng
1662f365b2 Andr*2335          if( pl(i).ge.p(L)   .and.
                2336      .       pl(i).le.p(L+1)) then
                2337          temp = ( pl(i)-p(L) )/( p(L+1)-p(L) )
                2338          acr(i) = a(L+1)*temp + a(L)*(1-temp)
                2339          endif
                2340       enddo
                2341       enddo
                2342 
471f7f004b Andr*2343       do i=1,lng
1662f365b2 Andr*2344       if( pl(i).lt.p(1)   ) acr(i) = a(1)
                2345       if( pl(i).gt.p(lma) ) acr(i) = a(lma)
                2346       enddo
                2347 
471f7f004b Andr*2348       do i=1,lng
1662f365b2 Andr*2349       acr(i) = acr(i) * (plb(i)-pl(i))
                2350       enddo
                2351 
                2352       RETURN
                2353       END
75f4744d22 Andr*2354        SUBROUTINE RNEVP(NN,IRUN,NLAY,TL,QL,RAIN,PL,CLFRAC,SP,DP,PLKE,
1662f365b2 Andr*2355      1   PLK,TH,TEMP1,TEMP2,TEMP3,ITMP1,ITMP2,RCON,RLAR,CLSBTH,tmscl,
                2356      2   tmfrc,cp,gravity,alhl,gamfac,cldlz,RHCRIT,offset,alpha)
a7a27977ab Jean*2357 
36f0442dce Andr*2358       implicit none
                2359 C Argument List declarations
                2360       integer nn,irun,nlay
a456aa407c Andr*2361       _RL TL(IRUN,NLAY),QL(IRUN,NLAY),RAIN(IRUN,NLAY),
36f0442dce Andr*2362      . PL(IRUN,NLAY),CLFRAC(IRUN,NLAY),SP(IRUN),TEMP1(IRUN,NLAY),
                2363      . TEMP2(IRUN,NLAY),PLKE(IRUN,NLAY+1),
                2364      . RCON(IRUN),RLAR(IRUN),DP(IRUN,NLAY),PLK(IRUN,NLAY),TH(IRUN,NLAY),
                2365      . TEMP3(IRUN,NLAY)
                2366       integer ITMP1(IRUN,NLAY),ITMP2(IRUN,NLAY)
a456aa407c Andr*2367       _RL CLSBTH(IRUN,NLAY)
                2368       _RL tmscl,tmfrc,cp,gravity,alhl,gamfac,offset,alpha
                2369       _RL cldlz(irun,nlay)
                2370       _RL rhcrit(irun,nlay)
36f0442dce Andr*2371 C
                2372 C Local Variables
f0525ae611 Andr*2373       _RL zm1p04,zero,two89,zp44,zp01,half,zp578,one,z3600,z1800
a456aa407c Andr*2374       _RL zp1,zp001
1662f365b2 Andr*2375       PARAMETER (ZM1P04 = -1.04E-4 )
                2376       PARAMETER (ZERO = 0.)
                2377       PARAMETER (TWO89= 2.89E-5)
                2378       PARAMETER ( ZP44= 0.44)
                2379       PARAMETER ( ZP01= 0.01)
                2380       PARAMETER ( ZP1 = 0.1 )
                2381       PARAMETER ( ZP001= 0.001)
                2382       PARAMETER ( HALF= 0.5)
                2383       PARAMETER ( ZP578 = 0.578 )
                2384       PARAMETER ( ONE = 1.)
                2385       PARAMETER ( Z3600 = 3600.)
f0525ae611 Andr*2386       PARAMETER ( Z1800 = 1800.)
1662f365b2 Andr*2387 C
a456aa407c Andr*2388       _RL EVP9(IRUN,NLAY)
                2389       _RL water(irun),crystal(irun)
                2390       _RL watevap(irun),iceevap(irun)
                2391       _RL fracwat,fracice, tice,rh,fact,dum
                2392       _RL rainmax(irun)
                2393       _RL getcon,rphf,elocp,cpog,relax
                2394       _RL exparg,arearat,rpow
1662f365b2 Andr*2395 
36f0442dce Andr*2396       integer i,L,n,nlaym1,irnlay,irnlm1
1662f365b2 Andr*2397 
a7a27977ab Jean*2398 C Explicit Inline Directives
                2399 C --------------------------
06d4643e1f Jean*2400 #ifdef FIZHI_CRAY
4e1c053948 Jean*2401 #ifdef FIZHI_F77_COMPIL
1662f365b2 Andr*2402 cfpp$ expand (qsat)
                2403 #endif
                2404 #endif
                2405 
                2406        tice = getcon('FREEZING-POINT')
                2407 
                2408        fracwat = 0.70
                2409        fracice = 0.01
                2410 
                2411        NLAYM1 = NLAY - 1
                2412        IRNLAY = IRUN*NLAY
                2413        IRNLM1 = IRUN*(NLAY-1)
a7a27977ab Jean*2414 
f0525ae611 Andr*2415 C      RPHF = Z3600/tmscl
                2416        RPHF = Z1800/tmscl
a7a27977ab Jean*2417 
1662f365b2 Andr*2418        ELOCP   = alhl/cp
                2419        CPOG    = cp/gravity
                2420 
                2421        DO I = 1,IRUN
                2422           RLAR(I) = 0.
                2423          water(i) = 0.
                2424        crystal(i) = 0.
                2425        ENDDO
                2426 
                2427        do L = 1,nlay
                2428        do i = 1,irun
                2429           EVP9(i,L) = 0.
                2430          TEMP1(i,L) = 0.
                2431          TEMP2(i,L) = 0.
                2432          TEMP3(i,L) = 0.
                2433         CLSBTH(i,L) = 0.
                2434          cldlz(i,L) = 0.
                2435        enddo
                2436        enddo
a7a27977ab Jean*2437 
1662f365b2 Andr*2438 C RHO(ZERO) / RHO FOR TERMINAL VELOCITY APPROX.
a7a27977ab Jean*2439 C ---------------------------------------------
1662f365b2 Andr*2440       DO L = 1,NLAY
                2441       DO I = 1,IRUN
                2442         TEMP2(I,L) = PL(I,L)*ZP001
                2443         TEMP2(I,L) = SQRT(TEMP2(I,L))
                2444       ENDDO
                2445       ENDDO
                2446 
                2447 C INVERSE OF MASS IN EACH LAYER
a7a27977ab Jean*2448 C -----------------------------
1662f365b2 Andr*2449       DO L = 1,NLAY
                2450       DO I = 1,IRUN
75f4744d22 Andr*2451       TEMP3(I,L) = GRAVITY*ZP01 / DP(I,L)
1662f365b2 Andr*2452       ENDDO
                2453       ENDDO
a7a27977ab Jean*2454 
1662f365b2 Andr*2455 C DO LOOP FOR MOISTURE EVAPORATION ABILITY AND CONVEC EVAPORATION.
a7a27977ab Jean*2456 C ----------------------------------------------------------------
1662f365b2 Andr*2457       DO 100 L=1,NLAY
                2458 
                2459       DO I = 1,IRUN
                2460         TEMP1(I,3) = TL(I,L)
                2461         TEMP1(I,4) = QL(I,L)
                2462       ENDDO
                2463 
                2464       DO 50 N=1,2
                2465        IF(N.EQ.1)RELAX=HALF
                2466        IF(N.GT.1)RELAX=ONE
                2467 
                2468       DO I = 1,IRUN
                2469        call qsat ( temp1(i,3),pl(i,L),temp1(i,2),temp1(i,6),.true. )
                2470        TEMP1(I,5)=TEMP1(I,2)-TEMP1(I,4)
                2471        TEMP1(I,6)=TEMP1(I,6)*ELOCP
                2472        TEMP1(I,5)=TEMP1(I,5)/(ONE+TEMP1(I,6))
                2473        TEMP1(I,4)=TEMP1(I,4)+TEMP1(I,5)*RELAX
                2474        TEMP1(I,3)=TEMP1(I,3)-TEMP1(I,5)*ELOCP*RELAX
                2475       ENDDO
                2476   50  CONTINUE
                2477 
                2478       DO I = 1,IRUN
                2479           EVP9(I,L) = (TEMP1(I,4) -  QL(I,L))/TEMP3(I,L)
                2480 C convective detrained water
a7a27977ab Jean*2481          cldlz(i,L) = rain(i,L)*temp3(i,L)
1662f365b2 Andr*2482        if(  tl(i,L).gt.tice-20.) then
                2483            water(i) =   water(i) + rain(i,L)
                2484        else
                2485          crystal(i) = crystal(i) + rain(i,L)
                2486        endif
                2487       ENDDO
a7a27977ab Jean*2488 
1662f365b2 Andr*2489 C**********************************************************************
                2490 C  FOR CONVECTIVE PRECIP, FIND THE "EVAPORATION EFFICIENCY" USING     *
                2491 C                   KESSLERS PARAMETERIZATION                         *
                2492 C**********************************************************************
a7a27977ab Jean*2493 
1662f365b2 Andr*2494       DO 20 I=1,IRUN
                2495 
                2496         iceevap(i) = 0.
                2497         watevap(i) = 0.
                2498 
                2499        if( (evp9(i,L).gt.0.) .and. (crystal(i).gt.0.) ) then
                2500        iceevap(I) = EVP9(I,L)*fracice
                2501        IF(iceevap(i).GE.crystal(i)) iceevap(i) = crystal(i)
                2502        EVP9(I,L)=EVP9(I,L)-iceevap(I)
                2503        crystal(i) = crystal(i) - iceevap(i)
                2504        endif
                2505 
                2506 C and now warm precipitate
                2507        if( (evp9(i,L).gt.0.) .and. (water(i).gt.0.) ) then
                2508        exparg    = ZM1P04*tmscl*((water(i)*RPHF*TEMP2(I,L))**ZP578)
                2509        AREARAT = ONE-(EXP(EXPARG))
                2510        watevap(I) = EVP9(I,L)*AREARAT*fracwat
                2511        IF(watevap(I).GE.water(i)) watevap(I) = water(i)
                2512        EVP9(I,L)=EVP9(I,L)-watevap(I)
                2513        water(i) = water(i) - watevap(i)
                2514        endif
                2515 
                2516        QL(I,L) = QL(I,L)+(iceevap(i)+watevap(i))*TEMP3(I,L)
                2517        TL(I,L) = TL(I,L)-(iceevap(i)+watevap(i))*TEMP3(I,L)*ELOCP
                2518 
                2519   20  CONTINUE
                2520 
                2521   100 CONTINUE
                2522 
                2523       do i = 1,irun
                2524       rcon(i) = water(i) + crystal(i)
                2525       enddo
                2526 
                2527 C**********************************************************************
                2528 C Large Scale Precip
                2529 C**********************************************************************
                2530 
                2531       DO 200 L=1,NLAY
                2532       DO I = 1,IRUN
                2533         rainmax(i) = rhcrit(i,L)*evp9(i,L) +
                2534      .               ql(i,L)*(rhcrit(i,L)-1.)/temp3(i,L)
                2535 
                2536          if (rainmax(i).LE.0.0) then
                2537            call qsat( tl(i,L),pl(i,L),rh,dum,.false.)
                2538            rh = ql(i,L)/rh
                2539 
                2540            if( rhcrit(i,L).eq.1.0 ) then
                2541            fact = 1.0
                2542            else
32362b8fa8 Cons*2543            fact = min( 1.0 _d 0, alpha + (1.0-alpha)*( rh-rhcrit(i,L)) /
1662f365b2 Andr*2544      1                                          (1.0-rhcrit(i,L)) )
                2545            endif
                2546 
                2547 C Do not allow clouds above 10 mb
a7a27977ab Jean*2548            if( pl(i,L).ge.10.0 ) CLSBTH(I,L) = fact
1662f365b2 Andr*2549            RLAR(I) = RLAR(I)-rainmax(I)
                2550            QL(I,L) = QL(I,L)+rainmax(I)*TEMP3(I,L)
                2551            TL(I,L) = TL(I,L)-rainmax(I)*TEMP3(I,L)*ELOCP
                2552 C Large-scale water
a7a27977ab Jean*2553         cldlz(i,L) = cldlz(i,L) - rainmax(i)*temp3(i,L)
1662f365b2 Andr*2554          ENDIF
                2555       ENDDO
                2556 
                2557       DO I=1,IRUN
                2558          IF((RLAR(I).GT.0.0).AND.(rainmax(I).GT.0.0))THEN
                2559           RPOW=(RLAR(I)*RPHF*TEMP2(I,L))**ZP578
                2560           EXPARG  = ZM1P04*tmscl*RPOW
                2561           AREARAT = ONE-(EXP(EXPARG))
                2562           TEMP1(I,7) = rainmax(I)*AREARAT
                2563           IF(TEMP1(I,7).GE.RLAR(I)) TEMP1(I,7) = RLAR(I)
                2564           RLAR(I) =    RLAR(I)-TEMP1(I,7)
                2565           QL(I,L) =    QL(I,L)+TEMP1(I,7)*TEMP3(I,L)
                2566           TL(I,L) =    TL(I,L)-TEMP1(I,7)*TEMP3(I,L)*ELOCP
                2567          ENDIF
                2568       ENDDO
                2569 
                2570   200 CONTINUE
                2571 
                2572       RETURN
                2573       END
                2574       subroutine srclouds (th,q,plk,pl,plke,cloud,cldwat,irun,irise,
                2575      1                                                rhc,offset,alpha)
                2576 C***********************************************************************
                2577 C
                2578 C  PURPOSE:
                2579 C  ========
                2580 C    Compute non-precipitating cloud fractions
                2581 C    based on Slingo and Ritter (1985).
                2582 C    Remove cloudiness where conditionally unstable.
                2583 C
                2584 C  INPUT:
                2585 C  ======
                2586 C    th ......... Potential Temperature (irun,irise)
                2587 C    q .......... Specific  Humidity    (irun,irise)
                2588 C    plk ........ P**Kappa at mid-layer (irun,irise)
                2589 C    pl ......... Pressure at mid-layer (irun,irise)
                2590 C    plke ....... P**Kappa at edge      (irun,irise+1)
                2591 C    irun ....... Horizontal   dimension
                2592 C    irise ...... Vertical     dimension
                2593 C
                2594 C  OUTPUT:
                2595 C  =======
                2596 C    cloud ...... Cloud Fraction        (irun,irise)
                2597 C
                2598 C***********************************************************************
                2599 
                2600       implicit none
                2601       integer  irun,irise
                2602 
a456aa407c Andr*2603       _RL   th(irun,irise)
                2604       _RL    q(irun,irise)
                2605       _RL  plk(irun,irise)
                2606       _RL   pl(irun,irise)
                2607       _RL plke(irun,irise+1)
1662f365b2 Andr*2608 
a456aa407c Andr*2609       _RL  cloud(irun,irise)
                2610       _RL cldwat(irun,irise)
                2611       _RL     qs(irun,irise)
1662f365b2 Andr*2612 
a456aa407c Andr*2613       _RL cp, alhl, getcon, akap
                2614       _RL ratio, temp, elocp
                2615       _RL rhcrit,rh,dum
36f0442dce Andr*2616       integer i,L
1662f365b2 Andr*2617 
a456aa407c Andr*2618       _RL rhc(irun,irise)
                2619       _RL offset,alpha
1662f365b2 Andr*2620 
a7a27977ab Jean*2621 C Explicit Inline Directives
                2622 C --------------------------
06d4643e1f Jean*2623 #ifdef FIZHI_CRAY
4e1c053948 Jean*2624 #ifdef FIZHI_F77_COMPIL
1662f365b2 Andr*2625 cfpp$ expand (qsat)
                2626 #endif
                2627 #endif
                2628 
                2629       cp     = getcon('CP')
                2630       alhl   = getcon('LATENT HEAT COND')
                2631       elocp  = alhl/cp
                2632       akap   = getcon('KAPPA')
                2633 
                2634       do L = 1,irise
                2635       do i = 1,irun
                2636                   temp = th(i,L)*plk(i,L)
                2637       call qsat ( temp,pl(i,L),qs(i,L),dum,.false. )
                2638       enddo
                2639       enddo
                2640 
                2641       do L  = 2,irise
                2642       do i  = 1,irun
                2643       rh    = q(i,L)/qs(i,L)
                2644 
                2645       rhcrit = rhc(i,L) - offset
                2646       ratio  = alpha*(rh-rhcrit)/offset
a7a27977ab Jean*2647 
1662f365b2 Andr*2648       if(cloud(i,L).eq.  0.0 .and. ratio.gt.0.0 ) then
32362b8fa8 Cons*2649          cloud(i,L) = min( ratio,1.0 _d 0)
1662f365b2 Andr*2650       endif
a7a27977ab Jean*2651 
1662f365b2 Andr*2652       enddo
                2653       enddo
a7a27977ab Jean*2654 
                2655 C Reduce clouds from conditionally unstable layer
                2656 C -----------------------------------------------
1662f365b2 Andr*2657       call ctei ( th,q,cloud,cldwat,pl,plk,plke,irun,irise )
                2658 
                2659       return
                2660       end
                2661 
                2662       subroutine ctei ( th,q,cldfrc,cldwat,pl,plk,plke,im,lm )
                2663       implicit none
                2664       integer im,lm
a456aa407c Andr*2665       _RL  th(im,lm),q(im,lm),plke(im,lm+1),cldwat(im,lm)
                2666       _RL plk(im,lm),pl(im,lm),cldfrc(im,lm)
1662f365b2 Andr*2667       integer i,L
a456aa407c Andr*2668       _RL    getcon,cp,alhl,elocp,cpoel,t,p,s,qs,dqsdt,dq
                2669       _RL    k,krd,kmm,f
1662f365b2 Andr*2670 
                2671       cp     = getcon('CP')
                2672       alhl   = getcon('LATENT HEAT COND')
                2673       cpoel  = cp/alhl
                2674       elocp  = alhl/cp
                2675 
                2676       do L=lm,2,-1
                2677       do i=1,im
                2678           dq = q(i,L)+cldwat(i,L)-q(i,L-1)-cldwat(i,L-1)
                2679       if( dq.eq.0.0 ) dq = 1.0e-20
                2680       k = 1.0 + cpoel*plke(i,L)*( th(i,L)-th(i,L-1) ) / dq
                2681 
                2682       t = th(i,L)*plk(i,L)
                2683       p = pl(i,L)
                2684       call qsat ( t,p,qs,dqsdt,.true. )
a7a27977ab Jean*2685 
1662f365b2 Andr*2686       krd = ( cpoel*t*(1+elocp*dqsdt) )/( 1 + 1.608*dqsdt*t )
                2687 
                2688       kmm = ( 1+elocp*dqsdt )*( 1 + 0.392*cpoel*t )
                2689      .    / ( 2+(1+1.608*cpoel*t)*elocp*dqsdt )
                2690 
                2691       s = ( (k-krd)/(kmm-krd) )
a7a27977ab Jean*2692 c     f = 1.0 - min( 1.0 _d 0, max(0.0 _d 0,1.0-exp(-s)) )
                2693 C- to avoid floating overflow in exp(-s):
                2694       f = 1.
                2695       if (s.GT.0. ) f = max( 0.0 _d 0, min( 1.0 _d 0, exp(-s)) )
1662f365b2 Andr*2696 
                2697       cldfrc(i,L) = cldfrc(i,L)*f
                2698       cldwat(i,L) = cldwat(i,L)*f
                2699 
                2700       enddo
                2701       enddo
                2702 
                2703       return
                2704       end
                2705 
                2706       subroutine back2grd(gathered,indeces,scattered,irun)
                2707       implicit none
                2708       integer i,irun,indeces(irun)
a456aa407c Andr*2709       _RL gathered(irun),scattered(irun)
                2710       _RL temp(irun)
1662f365b2 Andr*2711       do i = 1,irun
                2712        temp(indeces(i)) = gathered(i)
                2713       enddo
                2714       do i = 1,irun
                2715        scattered(i) = temp(i)
                2716       enddo
                2717       return
                2718       end