pkg/fizhi/fizhi_moist.F

C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_moist.F,v 1.24 2004/12/11 00:34:24 molod Exp $
C $Name:  $

#include "FIZHI_OPTIONS.h"
      subroutine moistio (ndmoist,istrip,npcs,
     .   lowlevel,midlevel,nltop,nsubmin,nsubmax,Lup,
     .   pz,plz,plze,dpres,pkht,pkl,uz,vz,tz,qz,bi,bj,ntracerin,ptracer,
     .   qqz,dumoist,dvmoist,dtmoist,dqmoist,cumfric,
     .   im,jm,lm,ptop,
     .   iras,rainlsp,rainconv,snowfall,
     .   nswcld,cldtot_sw,cldras_sw,cldlsp_sw,nswlz,swlz,
     .   nlwcld,cldtot_lw,cldras_lw,cldlsp_lw,nlwlz,lwlz,
     .   lpnt,myid)

       implicit none

#ifdef ALLOW_DIAGNOSTICS
#include "SIZE.h"
#include "DIAGNOSTICS_SIZE.h"
#include "DIAGNOSTICS.h"
#endif

c Input Variables
c ---------------
      integer im,jm,lm
      integer ndmoist,istrip,npcs
      integer bi,bj,ntracerin,ptracer         
      integer lowlevel,midlevel,nltop,nsubmin,nsubmax,Lup
      _RL pz(im,jm),plz(im,jm,lm),plze(im,jm,lm+1),dpres(im,jm,lm)
      _RL pkht(im,jm,lm+1),pkl(im,jm,lm)
      _RL tz(im,jm,lm),qz(im,jm,lm,ntracerin)       
      _RL uz(im,jm,lm),vz(im,jm,lm)       
      _RL qqz(im,jm,lm)
      _RL dumoist(im,jm,lm),dvmoist(im,jm,lm)
      _RL dtmoist(im,jm,lm),dqmoist(im,jm,lm,ntracerin)
      logical cumfric
      _RL ptop
      integer iras
      _RL rainlsp(im,jm),rainconv(im,jm),snowfall(im,jm)
      integer nswcld,nswlz
      _RL cldlsp_sw(im,jm,lm),cldras_sw(im,jm,lm)
      _RL cldtot_sw(im,jm,lm),swlz(im,jm,lm)
      integer nlwcld,nlwlz
      _RL  cldlsp_lw(im,jm,lm),cldras_lw(im,jm,lm)
      _RL  cldtot_lw(im,jm,lm),lwlz(im,jm,lm)
      logical lpnt
      integer myid

c Local Variables
c ---------------
      integer    ncrnd,nsecf

      _RL       fracqq, dum
      integer    snowcrit
      parameter (fracqq = 0.1)
      _RL one
      parameter (one=1.)

      _RL   cldsr(im,jm,lm)
      _RL   srcld(istrip,lm)

      _RL plev
      _RL cldnow,cldlsp_mem,cldlsp,cldras_mem,cldras
      _RL watnow,watmin,cldmin
      _RL cldprs(im,jm),cldtmp(im,jm)
      _RL cldhi (im,jm),cldlow(im,jm)
      _RL cldmid(im,jm),totcld(im,jm)

      _RL   CLDLS(im,jm,lm)  , CPEN(im,jm,lm)
      _RL    tmpimjm(im,jm)
      _RL    lsp_new(im,jm)
      _RL   conv_new(im,jm)
      _RL   snow_new(im,jm)

      _RL  qqcolmin(im,jm)
      _RL  qqcolmax(im,jm)
      integer levpbl(im,jm)

c Gathered Arrays for Variable Cloud Base
c ---------------------------------------
      _RL    raincgath(im*jm)
      _RL     pigather(im*jm)
      _RL     thgather(im*jm,lm)
      _RL     shgather(im*jm,lm)
      _RL    pkzgather(im*jm,lm)
      _RL    pkegather(im*jm,lm+1)
      _RL    plzgather(im*jm,lm)
      _RL    plegather(im*jm,lm+1)
      _RL     dpgather(im*jm,lm)
      _RL    tmpgather(im*jm,lm)
      _RL   deltgather(im*jm,lm)
      _RL   delqgather(im*jm,lm)
      _RL      ugather(im*jm,lm,ntracerin+2-ptracer)
      _RL   delugather(im*jm,lm,ntracerin+2-ptracer)
      _RL     deltrnev(im*jm,lm)
      _RL     delqrnev(im*jm,lm)

      integer  nindeces(lm)
      integer  pblindex(im*jm)
      integer levgather(im*jm)

c Stripped Arrays
c ---------------
      _RL saveth (istrip,lm)
      _RL saveq  (istrip,lm)
      _RL saveu  (istrip,lm,ntracerin+2-ptracer)
      _RL usubcl (istrip,   ntracerin+2-ptracer)

      _RL     ple(istrip,lm+1)
      _RL      dp(istrip,lm)
      _RL      TL(ISTRIP,lm)  , SHL(ISTRIP,lm)
      _RL      PL(ISTRIP,lm)  , PLK(ISTRIP,lm)
      _RL    PLKE(ISTRIP,lm+1)
      _RL      TH(ISTRIP,lm)  ,CVTH(ISTRIP,lm)
      _RL   CVQ(ISTRIP,lm)
      _RL      UL(ISTRIP,lm,ntracerin+2-ptracer)
      _RL     cvu(istrip,lm,ntracerin+2-ptracer)
      _RL  CLMAXO(ISTRIP,lm),CLBOTH(ISTRIP,lm)
      _RL  CLSBTH(ISTRIP,lm)
      _RL    TMP1(ISTRIP,lm),  TMP2(ISTRIP,lm)
      _RL    TMP3(ISTRIP,lm),  TMP4(ISTRIP,lm+1)
      _RL    TMP5(ISTRIP,lm+1)
      integer   ITMP1(ISTRIP,lm), ITMP2(ISTRIP,lm)

      _RL   PRECIP(ISTRIP), PCNET(ISTRIP)
      _RL   SP(ISTRIP),  PREP(ISTRIP)
      _RL   PCPEN (ISTRIP,lm)
      integer pbl(istrip),depths(lm)

      _RL   cldlz(istrip,lm), cldwater(im,jm,lm)
      _RL   rhfrac(istrip), rhmin, pup, ppbl, rhcrit(istrip,lm)
      _RL   offset, alpha, rasmax

      logical first
      logical lras
      _RL    clfrac (istrip,lm)
      _RL    cldmas (istrip,lm)
      _RL    detrain(istrip,lm)
      _RL    psubcld    (istrip), psubcldg (im,jm)
      _RL    psubcld_cnt(istrip), psubcldgc(im,jm)
      _RL rnd(lm/2)
      DATA      FIRST /.TRUE./

      integer imstp,nsubcl,nlras
      integer i,j,iloop,indx,indgath,l,nn,num,numdeps,nt
      _RL tmstp,tminv,sday,grav,alhl,cp,elocp,gamfac
      _RL rkappa,p0kappa,p0kinv,ptopkap,pcheck
      _RL tice,getcon,pi
      integer ntracer,ntracedim, ntracex
 
C **********************************************************************
C ****                     INITIALIZATION                           ****
C **********************************************************************

C Add U and V components to tracer array for cumulus friction

      if(cumfric) then
       ntracer = ntracerin + 2
      else
       ntracer = ntracerin
      endif
      ntracedim= max(ntracer-ptracer,1)
      ntracex= ntracer-ptracer
      IMSTP  = nsecf(NDMOIST)
      TMSTP  = FLOAT(IMSTP)
      TMINV  = 1. /  TMSTP

C Minimum Large-Scale Cloud Fraction at rhcrit
      alpha  = 0.80         
C Difference in fraction between SR and LS Threshold
      offset = 0.10         
C Large-Scale Relative Humidity Threshold in PBL
      rhmin  = 0.90         
C Maximum Cloud Fraction associated with RAS
      rasmax = 1.00         

          nn = 3*3600.0/tmstp + 1
C Threshold for Cloud Fraction Memory
      cldmin = rasmax*(1.0-tmstp/3600.)**nn   
C Threshold for Cloud Liquid Water Memory
      watmin = 1.0e-8                         

      SDAY   = GETCON('SDAY')
      GRAV   = GETCON('GRAVITY')
      ALHL   = GETCON('LATENT HEAT COND')
      CP     = GETCON('CP')
      ELOCP  = GETCON('LATENT HEAT COND') / GETCON('CP')
      GAMFAC = GETCON('LATENT HEAT COND') * GETCON('EPS') * ELOCP
     .       / GETCON('RGAS')
      RKAPPA = GETCON('KAPPA')
      P0KAPPA  =  1000.0**RKAPPA
      P0KINV   =  1. / P0KAPPA
      PTOPKAP  =  PTOP**RKAPPA
      tice     = getcon('FREEZING-POINT')
      PI       = 4.*atan(1.)

c Determine Total number of Random Clouds to Check
c ---------------------------------------------
      ncrnd = (lm-nltop+1)/2

      if(first .and. myid.eq.1 .and. bi.eq.1 ) then
       print *
       print *,'Top Level Allowed for Convection : ',nltop
       print *,'          Highest Sub-Cloud Level: ',nsubmax
       print *,'          Lowest  Sub-Cloud Level: ',nsubmin
       print *,'    Total Number of Random Clouds: ',ncrnd
       print *
       first = .false.
      endif

c And now find PBL depth - the level where qq = fracqq * qq at surface
c --------------------------------------------------------------------
      do j = 1,jm
      do i = 1,im
       qqcolmin(i,j) = qqz(i,j,lm)*fracqq
       qqcolmax(i,j) = qqz(i,j,lm)
         levpbl(i,j) = lm+1
      enddo
      enddo

      DO L = lm-1,1,-1
       DO j = 1,jm
       DO i = 1,im
        IF((qqz(i,j,l).gt.qqcolmax(i,j))
     1                   .and.(levpbl(i,j).eq.lm+1))then
         qqcolmax(i,j) = qqz(i,j,l)
         qqcolmin(i,j) = fracqq*qqcolmax(i,j)
        endif
        if((qqz(i,j,l).lt.qqcolmin(i,j))
     1                   .and.(levpbl(i,j).eq.lm+1))levpbl(i,j)=L+1
       enddo
       enddo
      enddo

      do j = 1,jm
      do i = 1,im
        if(levpbl(i,j).gt.nsubmin) levpbl(i,j) = nsubmin
        if(levpbl(i,j).lt.nsubmax) levpbl(i,j) = nsubmax
      enddo
      enddo


c Set up the array of indeces of subcloud levels for the gathering
c ----------------------------------------------------------------
      indx = 0
      do L = nsubmin,nltop,-1
       do j = 1,jm
       do i = 1,im
        if(levpbl(i,j).eq.L) then
         indx = indx + 1
         pblindex(indx) = (j-1)*im + i
        endif
       enddo
       enddo
      enddo

      do indx = 1,im*jm
       levgather(indx) = levpbl(pblindex(indx),1)
        pigather(indx) =     pz(pblindex(indx),1)
        pkegather(indx,lm+1) = pkht(pblindex(indx),1,lm+1)
        plegather(indx,lm+1) = plze(pblindex(indx),1,lm+1)
      enddo

      do L = 1,lm
       do indx = 1,im*jm
         thgather(indx,L) = tz(pblindex(indx),1,L)
         shgather(indx,L) = qz(pblindex(indx),1,L,1)
        pkegather(indx,L) = pkht(pblindex(indx),1,L)
        pkzgather(indx,L) = pkl(pblindex(indx),1,L)
        plegather(indx,L) = plze(pblindex(indx),1,L)
        plzgather(indx,L) = plz(pblindex(indx),1,L)
         dpgather(indx,L) = dpres(pblindex(indx),1,L)
       enddo
      enddo
C General Tracers
C----------------
      do nt = 1,ntracerin-ptracer
      do L = 1,lm
       do indx = 1,im*jm
        ugather(indx,L,nt) = qz(pblindex(indx),1,L,nt+ptracer)
       enddo
      enddo
      enddo

      if(cumfric) then
C Cumulus Friction - load u and v wind components into tracer array
C------------------------------------------------------------------
      do L = 1,lm
       do indx = 1,im*jm
        ugather(indx,L,ntracerin-ptracer+1) = uz(pblindex(indx),1,L)
        ugather(indx,L,ntracerin-ptracer+2) = vz(pblindex(indx),1,L)
       enddo
      enddo

      endif

c bump the counter for number of calls to convection
c --------------------------------------------------
                        iras = iras + 1
      if( iras.ge.1e9 ) iras = 1 

c select the 'random' cloud detrainment levels for RAS
c ----------------------------------------------------
      call rndcloud(iras,ncrnd,rnd,myid)
 
      do l=1,lm
      do j=1,jm
      do i=1,im
      dumoist(i,j,l) = 0.
      dvmoist(i,j,l) = 0.
      dtmoist(i,j,l) = 0.
        do nt = 1,ntracerin
        dqmoist(i,j,l,nt) = 0.
        enddo
      enddo
      enddo
      enddo

C***********************************************************************
C ****                     LOOP OVER NPCS PEICES                    ****
C **********************************************************************

      DO 1000 NN = 1,NPCS

C **********************************************************************
C ****                   VARIABLE INITIALIZATION                    ****
C **********************************************************************

       CALL STRIP (  pigather, SP      ,im*jm,ISTRIP,1 ,NN )
       CALL STRIP ( pkzgather, PLK     ,im*jm,ISTRIP,lm,NN )
       CALL STRIP ( pkegather, PLKE    ,im*jm,ISTRIP,lm+1,NN )
       CALL STRIP ( plzgather, PL      ,im*jm,ISTRIP,lm,NN )
       CALL STRIP ( plegather, PLE     ,im*jm,ISTRIP,lm+1,NN )
       CALL STRIP (  dpgather, dp      ,im*jm,ISTRIP,lm,NN )
       CALL STRIP (  thgather, TH      ,im*jm,ISTRIP,lm,NN )
       CALL STRIP (  shgather, SHL     ,im*jm,ISTRIP,lm,NN )
       CALL STRINT( levgather, pbl     ,im*jm,ISTRIP,1 ,NN )

       do nt = 1,ntracer-ptracer
        call strip ( ugather(1,1,nt), ul(1,1,nt),im*jm,istrip,lm,nn )
       enddo

C **********************************************************************
C ****        SETUP FOR RAS CUMULUS PARAMETERIZATION                ****
C **********************************************************************

      DO L = 1,lm
      DO I = 1,ISTRIP
             TH(I,L) = TH(I,L) * P0KAPPA
         CLMAXO(I,L) = 0.
         CLBOTH(I,L) = 0.
         cldmas(I,L) = 0.
        detrain(I,L) = 0.
      ENDDO
      ENDDO

      do L = 1,lm
      depths(L) = 0
      enddo

      numdeps = 0
      do L = nsubmin,nltop,-1
                       nindeces(L) = 0
       do i = 1,istrip
       if(pbl(i).eq.L) nindeces(L) = nindeces(L) + 1
       enddo
       if(nindeces(L).gt.0) then
       numdeps = numdeps + 1
       depths(numdeps) = L
       endif
      enddo

C Initiate a do-loop around RAS for the number of different
C    sub-cloud layer depths in this strip
C --If all subcloud depths are the same, execute loop once
C   Otherwise loop over different subcloud layer depths

      num = 1
      DO iloop = 1,numdeps

       nsubcl = depths(iloop)

c Compute sub-cloud values for Temperature and Spec.Hum.
c ------------------------------------------------------
       DO 600 I=num,num+nindeces(nsubcl)-1
        TMP1(I,2) = 0.
        TMP1(I,3) = 0.
 600   CONTINUE

       NLRAS = NSUBCL - NLTOP + 1
       DO 601 L=NSUBCL,lm
       DO 602 I=num,num+nindeces(nsubcl)-1
        TMP1(I,2) = TMP1(I,2) + (PLE(I,L+1)-PLE(I,L))*TH (I,L)/sp(i)
        TMP1(I,3) = TMP1(I,3) + (PLE(I,L+1)-PLE(I,L))*SHL(I,L)/sp(i)
 602   CONTINUE
 601   CONTINUE
       DO 603 I=num,num+nindeces(nsubcl)-1
        TMP1(I,4) = 1. / ( (PLE(I,lm+1)-PLE(I,NSUBCL))/sp(I) )
         TH(I,NSUBCL) = TMP1(I,2)*TMP1(I,4)
        SHL(I,NSUBCL) = TMP1(I,3)*TMP1(I,4)
 603   CONTINUE

c Save initial value of tracers and compute sub-cloud value
c ---------------------------------------------------------
       DO NT = 1,ntracer-ptracer
          do  L = 1,lm
          do  i = num,num+nindeces(nsubcl)-1
          saveu(i,L,nt) = ul(i,L,nt)
          enddo
          enddo
          DO I=num,num+nindeces(nsubcl)-1
          TMP1(I,2) = 0.
          ENDDO
          DO L=NSUBCL,lm
          DO I=num,num+nindeces(nsubcl)-1
           TMP1(I,2) = TMP1(I,2)+(PLE(I,L+1)-PLE(I,L))*UL(I,L,NT)/sp(i)
          ENDDO
          ENDDO
          DO I=num,num+nindeces(nsubcl)-1
          UL(I,NSUBCL,NT) = TMP1(I,2)*TMP1(I,4)
             usubcl(i,nt) = ul(i,nsubcl,nt)
          ENDDO
       ENDDO

c Compute Pressure Arrays for RAS
c -------------------------------
       DO 111 L=1,lm
       DO 112 I=num,num+nindeces(nsubcl)-1
        TMP4(I,L) = PLE(I,L)
 112   CONTINUE
 111   CONTINUE
       DO I=num,num+nindeces(nsubcl)-1
        TMP5(I,1) = PTOPKAP / P0KAPPA
       ENDDO
       DO L=2,lm
       DO I=num,num+nindeces(nsubcl)-1
        TMP5(I,L) = PLKE(I,L)*P0KINV
       ENDDO
       ENDDO
       DO  I=num,num+nindeces(nsubcl)-1
        TMP4(I,lm+1) = PLE (I,lm+1)
        TMP5(I,lm+1) = PLKE(I,lm+1)*P0KINV
       ENDDO
       DO 113 I=num,num+nindeces(nsubcl)-1
        TMP4(I,NSUBCL+1) = PLE (I,lm+1)
        TMP5(I,NSUBCL+1) = PLKE(I,lm+1)*P0KINV
 113   CONTINUE

      do i=num,num+nindeces(nsubcl)-1
C Temperature at top of sub-cloud layer
      tmp2(i,1) = TH(i,NSUBCL) * PLKE(i,NSUBCL)/P0KAPPA   
C Pressure    at top of sub-cloud layer
      tmp2(i,2) = tmp4(i,nsubcl)                          
      enddo

C  CHANGED THIS: no RH requirement for RAS
c     call vqsat ( tmp2(num,1),tmp2(num,2),tmp2(num,3),
c    .             dum,.false.,nindeces(nsubcl) )
c     do i=num,num+nindeces(nsubcl)-1
c     rh =  SHL(I,NSUBCL) / tmp2(i,3)
c       if (rh .le. 0.85) then
c         rhfrac(i) = 0.
c       else if (rh .ge. 0.95) then
c         rhfrac(i) = 1.
c       else
c         rhfrac(i) = (rh-0.85)*10.
c       endif
c     enddo
      do i=num,num+nindeces(nsubcl)-1
        rhfrac(i) = 1.
      enddo

C  Compute RH threshold for Large-scale condensation
C  Used in Slingo-Ritter clouds as well - define offset between SR and LS

C Top level of atan func above this rh_threshold = rhmin
      pup = 600.                
      do i=num,num+nindeces(nsubcl)-1
       do L = nsubcl, lm
         rhcrit(i,L) = 1.
       enddo
       do L = 1, nsubcl-1
        pcheck = pl(i,L)
        if (pcheck .le. pup) then
         rhcrit(i,L) = rhmin
        else
         ppbl = pl(i,nsubcl)
         rhcrit(i,L) = rhmin + (1.-rhmin)/(19.) * 
     .      ((atan( (2.*(pcheck-pup)/(ppbl-pup)-1.) *
     .       tan(20.*pi/21.-0.5*pi) ) 
     .        + 0.5*pi) * 21./pi - 1.)
        endif
       enddo
      enddo

c Save Initial Values of Temperature and Specific Humidity
c --------------------------------------------------------
      do L = 1,lm
      do i = num,num+nindeces(nsubcl)-1
        saveth(i,L) = th (i,L)
        saveq (i,L) = shl(i,L)
        PCPEN (i,L) = 0.
        CLFRAC(i,L) = 0.
      enddo
      enddo

      CALL RAS ( NN,istrip,nindeces(nsubcl),NLRAS,NLTOP,lm,TMSTP
     1, UL(num,1,1),ntracedim,ntracex,TH(num,NLTOP),SHL(num,NLTOP)
     2, TMP4(num,NLTOP), TMP5(num,NLTOP),rnd, ncrnd, PCPEN(num,NLTOP)
     3, CLBOTH(num,NLTOP), CLFRAC(num,NLTOP)
     4, cldmas(num,nltop), detrain(num,nltop)
     8, cp,grav,rkappa,alhl,rhfrac(num),rasmax )

c Compute Diagnostic CLDMAS in RAS Subcloud Layers
c ------------------------------------------------
       do L=nsubcl,lm
       do I=num,num+nindeces(nsubcl)-1
        dum = dp(i,L)/(ple(i,lm+1)-ple(i,nsubcl))
        cldmas(i,L) = cldmas(i,L-1) - dum*cldmas(i,nsubcl-1)
       enddo
       enddo

c Update Theta and Moisture due to RAS
c ------------------------------------
       DO L=1,nsubcl
       DO I=num,num+nindeces(nsubcl)-1
        CVTH(I,L) =  (TH (I,L) - saveth(i,l))
        CVQ (I,L) =  (SHL(I,L) - saveq (i,l))
       ENDDO
       ENDDO
       DO L=nsubcl+1,lm
       DO I=num,num+nindeces(nsubcl)-1
        CVTH(I,L) = cvth(i,nsubcl)
        CVQ (I,L) = cvq (i,nsubcl)
       ENDDO
       ENDDO

       DO L=nsubcl+1,lm
       DO I=num,num+nindeces(nsubcl)-1
        TH (I,L) = saveth(i,l) + cvth(i,l)
        SHL(I,L) = saveq (i,l) + cvq (i,l)
       ENDDO
       ENDDO
       DO L=1,lm
       DO I=num,num+nindeces(nsubcl)-1
        CVTH(I,L) =  CVTH(I,L) *P0KINV*SP(I)*tminv
        CVQ (I,L) =  CVQ (I,L) *SP(I)*tminv
       ENDDO
       ENDDO

c Compute Tracer Tendency due to RAS
c ----------------------------------
       do nt = 1,ntracer-ptracer
        DO L=1,nsubcl-1
        DO I=num,num+nindeces(nsubcl)-1
         CVU(I,L,nt) = ( UL(I,L,nt)-saveu(i,l,nt) )*sp(i)*tminv
        ENDDO
        ENDDO
        DO L=nsubcl,lm
        DO I=num,num+nindeces(nsubcl)-1
         if( usubcl(i,nt).ne.0.0 ) then
          cvu(i,L,nt) = ( ul(i,nsubcl,nt)-usubcl(i,nt) ) * 
     .                     ( saveu(i,L,nt)/usubcl(i,nt) )*sp(i)*tminv
         else
          cvu(i,L,nt) = 0.0
         endif
        ENDDO
        ENDDO
       enddo

c Compute Diagnostic PSUBCLD (Subcloud Layer Pressure)
c ----------------------------------------------------
       do i=num,num+nindeces(nsubcl)-1
                             lras = .false.
       do L=nltop,nsubcl
       if( cvq(i,L).ne.0.0 ) lras = .true.
       enddo
       psubcld    (i) = 0.0
       psubcld_cnt(i) = 0.0
       if( lras ) then 
       psubcld    (i) = sp(i)+ptop-ple(i,nsubcl)
       psubcld_cnt(i) = 1.0
       endif
       enddo
      

C End of subcloud layer depth loop (iloop)

       num = num+nindeces(nsubcl)

      ENDDO

C **********************************************************************
C ****                     TENDENCY UPDATES                         ****
C ****    (Keep 'Gathered' tendencies in 'gather' arrays now)       ****
C **********************************************************************

      call paste( CVTH,deltgather,istrip,im*jm,lm,NN )
      call paste(  CVQ,delqgather,istrip,im*jm,lm,NN )
      do nt = 1,ntracer-ptracer
      call paste( CVU(1,1,nt),delugather(1,1,nt),istrip,im*jm,lm,NN )
      enddo

C **********************************************************************
C     And now paste some arrays for filling diagnostics
C (use pkegather to hold detrainment and tmpgather for cloud mass flux)
C **********************************************************************

      if(icldmas .gt.0) call paste( cldmas,tmpgather,istrip,im*jm,lm,NN)
      if(idtrain .gt.0) call paste(detrain,pkegather,istrip,im*jm,lm,NN)
      if(ipsubcld.gt.0) then 
      call paste(psubcld    ,psubcldg ,istrip,im*jm,1,NN)
      call paste(psubcld_cnt,psubcldgc,istrip,im*jm,1,NN)
      endif

C *********************************************************************
C ****         RE-EVAPORATION OF PENETRATING CONVECTIVE RAIN       ****
C *********************************************************************
 
      CALL STRIP ( thgather,TH ,im*jm,ISTRIP,lm,NN)
      CALL STRIP ( shgather,SHL,im*jm,ISTRIP,lm,NN)
      DO L=1,lm
      DO I=1,ISTRIP
           TH(I,L) =  TH(I,L) + CVTH(I,L)*tmstp/SP(I)
          SHL(I,L) = SHL(I,L) +  CVQ(I,L)*tmstp/SP(I)
           TL(I,L) =  TH(I,L)*PLK(I,L)
       saveth(I,L) =  th(I,L)
       saveq (I,L) = SHL(I,L)
      ENDDO
      ENDDO
 
       CALL RNEVP (NN,ISTRIP,lm,TL,SHL,PCPEN,PL,CLFRAC,SP,DP,PLKE,
     .  PLK,TH,TMP1,TMP2,TMP3,ITMP1,ITMP2,PCNET,PRECIP,
     .  CLSBTH,TMSTP,one,cp,grav,alhl,gamfac,cldlz,rhcrit,offset,alpha)

C **********************************************************************
C ****                     TENDENCY UPDATES                         ****
C **********************************************************************

      DO L=1,lm

       DO I =1,ISTRIP
       TMP1(I,L) = sp(i) * (SHL(I,L)-saveq(I,L)) * tminv
       ENDDO
       CALL PSTBMP(TMP1(1,L),delqgather(1,L),ISTRIP,im*jm,1,NN)

       DO I =1,ISTRIP
       TMP1(I,L) = sp(i) * ((TL(I,L)/PLK(I,L))-saveth(i,l)) * tminv
       ENDDO
       CALL PSTBMP(TMP1(1,L),deltgather(1,L),ISTRIP,im*jm,1,NN)

C Paste rain evap tendencies into arrays for diagnostic output
c ------------------------------------------------------------
       if(idtls.gt.0)then
       DO I =1,ISTRIP
        TMP1(I,L) = ((TL(I,L)/PLK(I,L))-saveth(i,l))*plk(i,l)*sday*tminv
       ENDDO
       call paste(tmp1(1,L),deltrnev(1,L),istrip,im*jm,1,NN)
       endif

       if(idqls.gt.0)then
       DO I =1,ISTRIP
        TMP1(I,L) = (SHL(I,L)-saveq(I,L)) * 1000. * sday * tminv
       ENDDO
       call paste(tmp1(1,L),delqrnev(1,L),istrip,im*jm,1,NN)
       endif

      ENDDO

C *********************************************************************
C          Add Non-Precipitating Clouds where the relative 
C                     humidity is less than 100%
C               Apply Cloud Top Entrainment Instability
C *********************************************************************

      do L=1,lm
      do i=1,istrip
      srcld(i,L) = -clsbth(i,L)
      enddo
      enddo

      call srclouds (saveth,saveq,plk,pl,plke,clsbth,cldlz,istrip,lm,
     .                                              rhcrit,offset,alpha)

      do L=1,lm
      do i=1,istrip
      srcld(i,L) = srcld(i,L)+clsbth(i,L)
      enddo
      enddo

C *********************************************************************
C ****                PASTE CLOUD AMOUNTS                          ****
C *********************************************************************

      call paste (  srcld,   cldsr,istrip,im*jm,lm,nn )
      call paste (  cldlz,cldwater,istrip,im*jm,lm,nn )
      call paste ( clsbth,   cldls,istrip,im*jm,lm,nn )
      call paste ( clboth,   cpen ,istrip,im*jm,lm,nn )
 
c compute Total Accumulated Precip for Landsurface Model
c ------------------------------------------------------
      do i = 1,istrip
C  Initialize Rainlsp, Rainconv and Snowfall
      tmp1(i,1) = 0.0   
      tmp1(i,2) = 0.0 
      tmp1(i,3) = 0.0
      enddo

      do i = 1,istrip
      prep(i)   = PRECIP(I) + PCNET(I)
      tmp1(i,1) = PRECIP(I)
      tmp1(i,2) = pcnet(i)
      enddo
c
c  check whether there is snow
c-------------------------------------------------------
c  snow algorthm:
c  if temperature profile from the surface level to 700 mb 
c  uniformaly c  below zero, then precipitation (total) is 
c  snowfall.  Else there is no snow.
c-------------------------------------------------------

        do i = 1,istrip
          snowcrit=0
          do l=lup,lm
            if (saveth(i,l)*plk(i,l).le. tice ) then
             snowcrit=snowcrit+1
            endif
          enddo
          if (snowcrit .eq. (lm-lup+1)) then
            tmp1(i,3) = prep(i)
            tmp1(i,1)=0.0
            tmp1(i,2)=0.0
          endif
        enddo
  
      CALL paste (tmp1(1,1), lsp_new,ISTRIP,im*jm,1,NN)
      CALL paste (tmp1(1,2),conv_new,ISTRIP,im*jm,1,NN)
      CALL paste (tmp1(1,3),snow_new,ISTRIP,im*jm,1,NN)

      if(iprecon.gt.0) then
      CALL paste (pcnet,raincgath,ISTRIP,im*jm,1,NN)
      endif

C *********************************************************************
C ****               End Major Stripped Region                     ****
C *********************************************************************

 1000 CONTINUE

C Large Scale Rainfall, Conv rain, and snowfall
c ---------------------------------------------
      call back2grd ( lsp_new,pblindex, lsp_new,im*jm)
      call back2grd (conv_new,pblindex,conv_new,im*jm)
      call back2grd (snow_new,pblindex,snow_new,im*jm)

      if(iprecon.gt.0) then
      call back2grd (raincgath,pblindex,raincgath,im*jm)
      endif

c Subcloud Layer Pressure
c -----------------------
      if(ipsubcld.gt.0) then
      call back2grd (psubcldg ,pblindex,psubcldg ,im*jm)
      call back2grd (psubcldgc,pblindex,psubcldgc,im*jm)
      endif

      do L = 1,lm
C Delta theta,q, convective, max and ls clouds
c --------------------------------------------
       call back2grd (deltgather(1,L),pblindex, dtmoist(1,1,L)  ,im*jm)
       call back2grd (delqgather(1,L),pblindex, dqmoist(1,1,L,1),im*jm)
       call back2grd (    cpen(1,1,L),pblindex,    cpen(1,1,L)  ,im*jm)
       call back2grd (   cldls(1,1,L),pblindex,   cldls(1,1,L)  ,im*jm)
       call back2grd (cldwater(1,1,L),pblindex,cldwater(1,1,L)  ,im*jm)
       call back2grd ( pkzgather(1,L),pblindex, pkzgather(1,L)  ,im*jm)

C Diagnostics:
c ------------
       if(icldmas.gt.0)call back2grd(tmpgather(1,L),pblindex,
     .                                            tmpgather(1,L),im*jm)
       if(idtrain.gt.0)call back2grd(pkegather(1,L),pblindex,
     .                                            pkegather(1,L),im*jm)
       if(idtls.gt.0)call back2grd(deltrnev(1,L),pblindex,
     .                                             deltrnev(1,L),im*jm)
       if(idqls.gt.0)call back2grd(delqrnev(1,L),pblindex,
     .                                             delqrnev(1,L),im*jm)
       if(icldnp.gt.0)call back2grd(cldsr(1,1,L),pblindex,
     .                                              cldsr(1,1,L),im*jm)
      enddo

c General Tracers
c ---------------
      do nt = 1,ntracerin-ptracer
       do L = 1,lm
       call back2grd (delugather(1,L,nt),pblindex,
     .                                 dqmoist(1,1,L,ptracer+nt),im*jm)
       enddo
      enddo

      if(cumfric) then

c U and V for cumulus friction
c ----------------------------
      do L = 1,lm
       call back2grd (delugather(1,L,ntracerin-ptracer+1),pblindex,
     .                                 dumoist(1,1,L),im*jm)
       call back2grd (delugather(1,L,ntracerin-ptracer+2),pblindex,
     .                                 dvmoist(1,1,L),im*jm)
      enddo

C Remove pi-weighting for u and v tendencies
      do j = 1,jm
      do i = 1,im
       tmpimjm(i,j) = 1./pz(i,j)
      enddo
      enddo
      do L = 1,lm
      do j = 1,jm
      do i = 1,im
       dumoist(i,j,L) = dumoist(i,j,L) * tmpimjm(i,j)
       dvmoist(i,j,L) = dvmoist(i,j,L) * tmpimjm(i,j)
      enddo
      enddo
      enddo


      endif

C **********************************************************************
C                          BUMP DIAGNOSTICS
C **********************************************************************


c Sub-Cloud Layer
c -------------------------
      if( ipsubcld.ne.0 ) then
      do j = 1,jm
      do i = 1,im
      qdiag(i,j,ipsubcld,bi,bj) = qdiag(i,j,ipsubcld,bi,bj) + 
     .                                           psubcldg (i,j)
      qdiag(i,j,ipsubcldc,bi,bj) = qdiag(i,j,ipsubcldc,bi,bj) + 
     .                                           psubcldgc(i,j)
      enddo
      enddo
      endif
 
c Non-Precipitating Cloud Fraction
c --------------------------------
      if( icldnp.ne.0 ) then
      do L = 1,lm
      do j = 1,jm
      do i = 1,im
      qdiag(i,j,icldnp+L-1,bi,bj) = qdiag(i,j,icldnp+L-1,bi,bj) + 
     .                                                  cldsr(i,j,L)
      enddo
      enddo
      enddo
      ncldnp = ncldnp + 1
      endif

c Moist Processes Heating Rate
c ----------------------------
      if(imoistt.gt.0) then
      do L = 1,lm
      do j = 1,jm
      do i = 1,im
       indgath = (j-1)*im + i
      qdiag(i,j,imoistt+L-1,bi,bj) = qdiag(i,j,imoistt+L-1,bi,bj) + 
     .    (dtmoist(i,j,L)*sday*pkzgather(indgath,L)/pz(i,j))
      enddo
      enddo
      enddo
      endif

c Moist Processes Moistening Rate
c -------------------------------
      if(imoistq.gt.0) then
      do L = 1,lm
      do j = 1,jm
      do i = 1,im
      qdiag(i,j,imoistq+L-1,bi,bj) = qdiag(i,j,imoistq+L-1,bi,bj) + 
     .                           (dqmoist(i,j,L,1)*sday*1000.0/pz(i,j))
      enddo
      enddo
      enddo
      endif

c Moist Processes Change in U-Momentum (Cumulus Friction)
c ------------------------------------------------------
      if(iudiag1.gt.0) then
      do L = 1,lm
      do j = 1,jm
      do i = 1,im
       indgath = (j-1)*im + i
      qdiag(i,j,iudiag1+L-1,bi,bj) = qdiag(i,j,iudiag1+L-1,bi,bj) + 
     .    dumoist(i,j,L)*sday
      enddo
      enddo
      enddo
      endif

c Moist Processes Change in V-Momentum (Cumulus Friction)
c ------------------------------------------------------
      if(iudiag2.gt.0) then
      do L = 1,lm
      do j = 1,jm
      do i = 1,im
       indgath = (j-1)*im + i
      qdiag(i,j,iudiag2+L-1,bi,bj) = qdiag(i,j,iudiag2+L-1,bi,bj) + 
     .    dvmoist(i,j,L)*sday
      enddo
      enddo
      enddo
      endif

c Cloud Mass Flux
c ---------------
      if(icldmas.gt.0) then
      do L = 1,lm
      do j = 1,jm
      do i = 1,im
       indgath = (j-1)*im + i
      qdiag(i,j,icldmas+L-1,bi,bj) = qdiag(i,j,icldmas+L-1,bi,bj) + 
     .                               tmpgather(indgath,L)
      enddo
      enddo
      enddo
      endif

c Detrained Cloud Mass Flux
c -------------------------
      if(idtrain.gt.0) then
      do L = 1,lm
      do j = 1,jm
      do i = 1,im
       indgath = (j-1)*im + i
      qdiag(i,j,idtrain+L-1,bi,bj) = qdiag(i,j,idtrain+L-1,bi,bj) + 
     .                                pkegather(indgath,L)
      enddo
      enddo
      enddo
      endif

c Grid-Scale Condensational Heating Rate
c --------------------------------------
      if(idtls.gt.0) then
      do L = 1,lm
      do j = 1,jm
      do i = 1,im
       indgath = (j-1)*im + i
      qdiag(i,j,idtls+L-1,bi,bj) = qdiag(i,j,idtls+L-1,bi,bj) + 
     .                               deltrnev(indgath,L)
      enddo
      enddo
      enddo
      endif

c Grid-Scale Condensational Moistening Rate
c -----------------------------------------
      if(idqls.gt.0) then
      do L = 1,lm
      do j = 1,jm
      do i = 1,im
       indgath = (j-1)*im + i
      qdiag(i,j,idqls+L-1,bi,bj) = qdiag(i,j,idqls+L-1,bi,bj) + 
     .                                delqrnev(indgath,L)
      enddo
      enddo
      enddo
      endif

c Total Precipitation
c -------------------
      if(ipreacc.gt.0) then
      do j = 1,jm
      do i = 1,im
      qdiag(i,j,ipreacc,bi,bj) = qdiag(i,j,ipreacc,bi,bj)
     .                   +  (  lsp_new(I,j)
     .                      + snow_new(I,j)
     .                      + conv_new(i,j) ) *sday*tminv
      enddo
      enddo
      endif

c Convective Precipitation
c ------------------------
      if(iprecon.gt.0) then
      do j = 1,jm
      do i = 1,im
       indgath = (j-1)*im + i
      qdiag(i,j,iprecon,bi,bj) = qdiag(i,j,iprecon,bi,bj) + 
     .                      raincgath(indgath)*sday*tminv
      enddo
      enddo
      endif

C **********************************************************************
C ****   Fill Rainfall and Snowfall Arrays for Land Surface Model   ****
C ****        Note:  Precip Rates work when DT(turb)<DT(moist)      ****
C **********************************************************************

      do j = 1,jm
      do i = 1,im
       rainlsp (i,j) = rainlsp (i,j) +  lsp_new(i,j)*tminv
       rainconv(i,j) = rainconv(i,j) + conv_new(i,j)*tminv
       snowfall(i,j) = snowfall(i,j) + snow_new(i,j)*tminv
      enddo
      enddo

C **********************************************************************
C ***        Compute Time-averaged Quantities for Radiation          ***
C ***               CPEN  => Cloud Fraction from RAS                 ***
C ***               CLDLS => Cloud Fraction from RNEVP               ***
C **********************************************************************

      do j = 1,jm
      do i = 1,im
          cldhi (i,j) = 0.
          cldmid(i,j) = 0.
          cldlow(i,j) = 0.
          cldtmp(i,j) = 0.
          cldprs(i,j) = 0.
         tmpimjm(i,j) = 0.
      enddo
      enddo

c Set Moist-Process Memory Coefficient
c ------------------------------------
      cldras_mem = 1.0-tmstp/ 3600.0
      cldlsp_mem = 1.0-tmstp/(3600.0*3)

      do L = 1,lm
      do i = 1,im*jm
       plev = pl(i,L)
 
c Compute Time-averaged Cloud and Water Amounts for Longwave Radiation
c --------------------------------------------------------------------
       watnow = cldwater(i,1,L)
       if( plev.le.500.0 ) then
        cldras = min( max( cldras_lw(i,1,L)*cldras_mem,cpen(i,1,L)),1.0)
       else
        cldras = 0.0
       endif
       cldlsp = min( max( cldlsp_lw(i,1,L)*cldlsp_mem,cldls(i,1,L)),1.0)

       if( cldras.lt.cldmin ) cldras = 0.0
       if( cldlsp.lt.cldmin ) cldlsp = 0.0

       cldnow = max( cldlsp,cldras )

       lwlz(i,1,L) = ( nlwlz*lwlz(i,1,L) + watnow)/(nlwlz +1)
       cldtot_lw(i,1,L) = (nlwcld*cldtot_lw(i,1,L) + cldnow)/(nlwcld+1)
       cldlsp_lw(i,1,L) = (nlwcld*cldlsp_lw(i,1,L) + cldlsp)/(nlwcld+1)
       cldras_lw(i,1,L) = (nlwcld*cldras_lw(i,1,L) + cldras)/(nlwcld+1)

 
c Compute Time-averaged Cloud and Water Amounts for Shortwave Radiation
c ---------------------------------------------------------------------
       watnow = cldwater(i,1,L)
       if( plev.le.500.0 ) then
        cldras = min( max(cldras_sw(i,1,L)*cldras_mem, cpen(i,1,L)),1.0)
       else
        cldras = 0.0
       endif
       cldlsp = min( max(cldlsp_sw(i,1,L)*cldlsp_mem,cldls(i,1,L)),1.0)

       if( cldras.lt.cldmin ) cldras = 0.0
       if( cldlsp.lt.cldmin ) cldlsp = 0.0

       cldnow = max( cldlsp,cldras )

       swlz(i,1,L) = ( nswlz*swlz(i,1,L) + watnow)/(nswlz +1)
       cldtot_sw(i,1,L) = (nswcld*cldtot_sw(i,1,L) + cldnow)/(nswcld+1)
       cldlsp_sw(i,1,L) = (nswcld*cldlsp_sw(i,1,L) + cldlsp)/(nswcld+1)
       cldras_sw(i,1,L) = (nswcld*cldras_sw(i,1,L) + cldras)/(nswcld+1)

 
c Compute Instantaneous Low-Mid-High Maximum Overlap Cloud Fractions
c ----------------------------------------------------------------------

       if( L.lt.midlevel ) cldhi (i,1) = max( cldnow,cldhi (i,1) )
       if( L.ge.midlevel  .and.
     .      L.lt.lowlevel ) cldmid(i,1) = max( cldnow,cldmid(i,1) )
       if( L.ge.lowlevel ) cldlow(i,1) = max( cldnow,cldlow(i,1) )

c Compute Cloud-Top Temperature and Pressure
c ------------------------------------------
       cldtmp(i,1) =  cldtmp(i,1) + cldnow*pkzgather(i,L)
     .             *  ( tz(i,1,L) + dtmoist(i,1,L)*tmstp/pz(i,1) )
       cldprs(i,1) =  cldprs(i,1) + cldnow*plev
       tmpimjm(i,1) = tmpimjm(i,1) + cldnow

      enddo
      enddo
 
c Compute Instantaneous Total 2-D Cloud Fraction
c ----------------------------------------------
      do j = 1,jm
      do i = 1,im
      totcld(i,j) = 1.0 - (1.-cldhi (i,j))
     .                  * (1.-cldmid(i,j))
     .                  * (1.-cldlow(i,j))
      enddo
      enddo
       

C **********************************************************************
C ***       Fill Cloud Top Pressure and Temperature Diagnostic       ***
C **********************************************************************

      if(icldtmp.gt.0) then
      do j = 1,jm
      do i = 1,im
         if( cldtmp(i,j).gt.0.0 ) then
         qdiag(i,j,icldtmp,bi,bj) = qdiag(i,j,icldtmp,bi,bj) + 
     .                       cldtmp(i,j)*totcld(i,j)/tmpimjm(i,j)
         qdiag(i,j,icttcnt,bi,bj) = qdiag(i,j,icttcnt,bi,bj) + 
     .                                                totcld(i,j)
         endif
      enddo
      enddo
      endif

      if(icldprs.gt.0) then
      do j = 1,jm
      do i = 1,im
         if( cldprs(i,j).gt.0.0 ) then
         qdiag(i,j,icldprs,bi,bj) = qdiag(i,j,icldprs,bi,bj) + 
     .                       cldprs(i,j)*totcld(i,j)/tmpimjm(i,j)
         qdiag(i,j,ictpcnt,bi,bj) = qdiag(i,j,ictpcnt,bi,bj) + 
     .                                                totcld(i,j)
         endif
      enddo
      enddo
      endif
   
C **********************************************************************
C ****                      INCREMENT COUNTERS                      ****
C **********************************************************************
 
       nlwlz    = nlwlz    + 1
       nswlz    = nswlz    + 1

       nlwcld   = nlwcld   + 1
       nswcld   = nswcld   + 1

#ifdef ALLOW_DIAGNOSTICS
       if( (bi.eq.1) .and. (bj.eq.1) ) then
       nmoistt  = nmoistt  + 1
       nmoistq  = nmoistq  + 1
       npreacc  = npreacc  + 1
       nprecon  = nprecon  + 1

       ncldmas  = ncldmas  + 1
       ndtrain  = ndtrain  + 1

       ndtls  = ndtls  + 1
       ndqls  = ndqls  + 1

       nudiag1  = nudiag1  + 1
       nudiag2  = nudiag2  + 1

       endif
#endif
 
      RETURN
      END
      SUBROUTINE RAS( NN, LNG, LENC, K, NLTOP, nlayr, DT
     *,      UOI, ntracedim, ntracer, POI, QOI, PRS, PRJ, rnd, ncrnd
     *,      RAINS, CLN, CLF, cldmas, detrain
     *,      cp,grav,rkappa,alhl,rhfrac,rasmax )
C
C*********************************************************************
C********************* SUBROUTINE  RAS   *****************************
C********************** 16 MARCH   1988 ******************************
C*********************************************************************
C
      implicit none

C Argument List
      integer nn,lng,lenc,k,nltop,nlayr
      integer ntracedim, ntracer
      integer ncrnd
      _RL dt
      _RL UOI(lng,nlayr,ntracedim),   POI(lng,K)
      _RL QOI(lng,K), PRS(lng,K+1), PRJ(lng,K+1)
      _RL rnd(ncrnd)
      _RL RAINS(lng,K), CLN(lng,K), CLF(lng,K)
      _RL cldmas(lng,K), detrain(lng,K)
      _RL cp,grav,rkappa,alhl,rhfrac(lng),rasmax

C Local Variables
      _RL TCU(lng,K), QCU(lng,K)
      _RL ucu(lng,K,ntracedim)
      _RL ALF(lng,K), BET(lng,K), GAM(lng,K)
     *,         ETA(lng,K), HOI(lng,K)
     *,         PRH(lng,K), PRI(lng,K)
      _RL HST(lng,K), QOL(lng,K), GMH(lng,K)
 
      _RL TX1(lng), TX2(lng), TX3(lng), TX4(lng), TX5(lng)
     *,         TX6(lng), TX7(lng), TX8(lng), TX9(lng)
     *,         TX11(lng), TX12(lng), TX13(lng), TX14(lng,ntracedim)
     *,         TX15(lng)
     *,         WFN(lng) 
      integer IA1(lng), IA2(lng), IA3(lng)
      _RL cloudn(lng), pcu(lng)

      integer krmin,icm
      _RL rknob, cmb2pa
      PARAMETER (KRMIN=01)
      PARAMETER (ICM=1000)
      PARAMETER (CMB2PA=100.0)
      PARAMETER (rknob = 10.)

      integer IC(ICM),   IRND(icm)
      _RL cmass(lng,K)
      LOGICAL SETRAS
      integer ifound
      _RL temp
      _RL thbef(lng,K)

      integer i,L,nc,ib,nt
      integer km1,kp1,kprv,kcr,kfx,ncmx
      _RL p00, crtmsf, frac, rasblf

      do L = 1,k
      do I = 1,LENC
       rains(i,l) = 0.
      enddo
      enddo
 
      p00 = 1000.
      crtmsf = 0.

C The numerator here is the fraction of the subcloud layer mass flux
C      allowed to entrain into the cloud

CCC   FRAC = 1./dt
CCC   FRAC = 0.5/dt
      FRAC = 0.5/dt

      KM1    = K  - 1
      KP1    = K  + 1
C we want the ras adjustment time scale to be one hour (indep of dt)
      RASBLF = 1./3600.
C
       KPRV  = KM1
C Removed KRMAX parameter
       KCR   = MIN(KM1,nlayr-2)     
       KFX   = KM1 - KCR
       NCMX  = KFX + NCRND
C
       IF (KFX .GT. 0) THEN
        DO NC=1,KFX
         IC(NC) = K - NC
        ENDDO
       ENDIF
C
      IF (NCRND .GT. 0) THEN
       DO I=1,ncrnd
        IRND(I) = (RND(I)-0.0005)*(KCR-KRMIN+1)
        IRND(I) = IRND(I) + KRMIN
       ENDDO
C
       DO NC=1,NCRND
        IC(KFX+NC) = IRND(NC)
       ENDDO
      ENDIF
C
      DO 100 NC=1,NCMX
C
       IF (NC .EQ. 1 ) THEN
        SETRAS = .TRUE.
       ELSE
        SETRAS = .FALSE.
       ENDIF
       IB = IC(NC)

c Initialize Cloud Fraction Array
c -------------------------------
      do i = 1,lenc
      cloudn(i) = 0.0
      enddo

       CALL CLOUD(nn,lng, LENC, K, NLTOP, nlayr, IB, RASBLF,SETRAS,FRAC
     *,           CP,  ALHL, RKAPPA, GRAV, P00, CRTMSF
     *,           POI, QOI, UOI, ntracedim, Ntracer, PRS, PRJ
     *,           PCU, CLOUDN, TCU, QCU, UCU, CMASS
     *,           ALF, BET, GAM, PRH, PRI, HOI, ETA
     *,           HST, QOL, GMH
     *,           TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, TX9
     *,           WFN, TX11, TX12, TX13, TX14, TX15
     *,           IA1,IA2,IA3,rhfrac)

C Compute fraction of grid box into which rain re-evap occurs (clf)
c -----------------------------------------------------------------
      do i = 1,lenc

c mass in detrainment layer
c -------------------------
       tx1(i) = cmb2pa *  (prs(i,ib+1) - prs(i,ib))/(grav*dt)

c ratio of detraining cloud mass to mass in detrainment layer
c -----------------------------------------------------------
       tx1(i) = rhfrac(i)*rknob * cmass(i,ib) / tx1(i)
       if(cmass(i,K).gt.0.) clf(i,ib) = clf(i,ib) + tx1(i)
       if( clf(i,ib).gt.1.) clf(i,ib) = 1.
      enddo
 
c Compute Total Cloud Mass Flux
c *****************************
      do L=ib,k
      do i=1,lenc
       cmass(i,L) = rhfrac(i)*cmass(i,L) * dt
      enddo
      enddo
         
      do L=ib,k
      do i=1,lenc
       cldmas(i,L) = cldmas(i,L) + cmass(i,L)
      enddo
      enddo

      do i=1,lenc
      detrain(i,ib) = detrain(i,ib) + cmass(i,ib)
      enddo

      DO L=IB,K
       DO I=1,LENC
        thbef(I,L) = POI(I,L)
        POI(I,L) = POI(I,L) + TCU(I,L) * DT * rhfrac(i)
        QOI(I,L) = QOI(I,L) + QCU(I,L) * DT * rhfrac(i)
       ENDDO
      ENDDO
      DO NT=1,Ntracer
      DO L=IB,K
       DO I=1,LENC
        UOI(I,L+nltop-1,NT)=UOI(I,L+nltop-1,NT)+UCU(I,L,NT)*DT*rhfrac(i)
       ENDDO
      ENDDO
      ENDDO
      DO I=1,LENC
       rains(I,ib) = rains(I,ib) + PCU(I)*dt * rhfrac(i)
      ENDDO

      do i = 1,lenc
       ifound = 0
       do L = 1,k
        if(tcu(i,L).ne.0.)ifound = ifound + 1
       enddo
       if(ifound.ne.0) then
c       print *,i,' made a cloud detraining at ',ib
        do L = 1,k
         temp = TCU(I,L) * DT * rhfrac(i)
c        write(6,122)L,thbef(i,L),poi(i,L),temp
        enddo
       endif
      enddo

  100 CONTINUE

 122  format(' ',i3,' TH B ',e10.3,' TH A ',e10.3,' DTH ',e10.3)
 
c Fill Convective Cloud Fractions based on 3-D Rain Amounts
c ---------------------------------------------------------
      do L=k-1,1,-1
      do i=1,lenc
         tx1(i)   = 100*(prs(i,L+1)-prs(i,L))/grav
         cln(i,L) = min(1600*rains(i,L)/tx1(i),rasmax )
      enddo
      enddo

      RETURN
      END
      subroutine rndcloud (iras,nrnd,rnd,myid)
      implicit none
      integer n,iras,nrnd,myid
      _RL random_numbx
      _RL rnd(nrnd)
      integer irm
      parameter (irm = 1000)
      _RL random(irm)
      integer i,mcheck,numrand,iseed,indx
      logical first
      data    first /.true./
      integer iras0
      data    iras0 /0/
      save random, iras0
 
      if(nrnd.eq.0.)then
       do i = 1,nrnd
        rnd(i) = 0
       enddo
       if(first .and. myid.eq.1) print *,' NO RANDOM CLOUDS IN RAS '
       go to 100
      endif

      mcheck = mod(iras-1,irm/nrnd)

c First Time In From a Continuing RESTART (IRAS.GT.1) or Reading a New RESTART
c ----------------------------------------------------------------------------
      if( first.and.(iras.gt.1) .or. iras.ne.iras0+1 )then
       print *,' first ',first,' iras ',iras,' iras0 ',iras0
       if( myid.eq.1 ) print *, 'Recreating Rand Numb Array in RNDCLOUD'
       if( myid.eq.1 ) print *, 'IRAS: ',iras,'  IRAS0: ',iras0
       numrand = mod(iras,irm/nrnd) * nrnd
       iseed   = iras * nrnd - numrand
       call random_seedx(iseed)
       do i = 1,irm
        random(i) = random_numbx(iseed)
       enddo
       indx = (iras-1)*nrnd

c Multiple Time In But have Used Up all 1000 numbers (MCHECK.EQ.0)
c ----------------------------------------------------------------
      else if (mcheck.eq.0) then
          iseed = (iras-1)*nrnd
          call random_seedx(iseed)
          do i = 1,irm
           random(i) = random_numbx(iseed)
          enddo
          indx = iseed

c Multiple Time In But have NOT Used Up all 1000 numbers (MCHECK.NE.0)
c --------------------------------------------------------------------
      else
          indx = (iras-1)*nrnd
      endif

          indx = mod(indx,irm)
      if( indx+nrnd.gt.1000 ) indx=1000-nrnd
 
      do n = 1,nrnd
       rnd(n) = random(indx+n)
      enddo
 
 100  continue
      first = .false.
      iras0 = iras
      return
      end
      function random_numbx(iseed)
      implicit none
      integer iseed
      real *8 seed,port_rand
      _RL random_numbx
      random_numbx = 0
#ifdef CRAY
      _RL ranf
      random_numbx = ranf()
#else 
#ifdef SGI
      _RL rand
      random_numbx = rand()
#endif
      random_numbx = port_rand(seed)
#endif
      return
      end
      subroutine random_seedx (iseed)
      implicit none
      integer  iseed
#ifdef CRAY
      call ranset (iseed)
#endif
#ifdef SGI
      integer*4   seed
                  seed = iseed
      call srand (seed)
#endif
      return
      end
      SUBROUTINE CLOUD(nn,lng, LENC, K, NLTOP, nlayr, IC, RASALF
     *,                 SETRAS, FRAC
     *,                 CP,  ALHL, RKAP, GRAV, P00, CRTMSF
     *,                 POI, QOI, UOI, ntracedim, Ntracer, PRS,  PRJ
     *,                 PCU, CLN, TCU, QCU, UCU, CMASS
     *,                 ALF, BET, GAM, PRH, PRI, HOL, ETA
     *,                 HST, QOL, GMH
     *,                 TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, ALM
     *,                 WFN, AKM, QS1, CLF, UHT, WLQ
     *,                 IA, I1, I2,rhfrac)
C
C*********************************************************************
C******************** Relaxed Arakawa-Schubert ***********************
C********************* Plug Compatible Version  **********************
C************************ SUBROUTINE CLOUD ***************************
C************************* 23 JULY    1992 ***************************
C*********************************************************************
C*********************************************************************
C*********************************************************************
C************************** Developed By *****************************
C**************************              *****************************
C************************ Shrinivas Moorthi **************************
C************************       and         **************************
C************************  Max J. Suarez *****************************
C************************                *****************************
C******************** Laboratory for Atmospheres *********************
C****************** NASA/GSFC, Greenbelt, MD 20771 *******************
C*********************************************************************
C*********************************************************************
C
C   The calculations of Moorthi and Suarez (1992, MWR) are 
C   contained in the CLOUD routine.
C   It is probably advisable, at least initially, to treat CLOUD
C   as a black box that computes the single cloud adjustments. RAS,
C   on the other hand, can be tailored to each GCMs configuration
C   (ie, number and placement of levels, nature of boundary layer,
C    time step and frequency with which RAS is called).
C
C
C  Input:
C  ------
C
C     lng     : The inner dimension of update and input arrays.
C
C     LENC    : The run: the number of soundings processes in a single call.
C               RAS works on the first LENC of the lng soundings
C               passed. This allows working on pieces of the world
C               say for multitasking, without declaring temporary arrays
C               and copying the data to and from them.  This is an f77
C               version. An F90 version would have to allow more
C               flexibility in the argument declarations.  Obviously
C               (LENC<=lng).  
C 
C     K       : Number of vertical layers (increasing downwards).
C               Need not be the same as the number of layers in the
C               GCM, since it is the outer dimension. The bottom layer
C               (K) is the subcloud layer. 
C
C     IC      : Detrainment level to check for presence of convection
C
C     RASALF  : Relaxation parameter (< 1.) for present cloud-type
C
C     SETRAS  : Logical parameter to control re-calculation of
C               saturation specific humidity and mid level P**kappa
C
C     FRAC    : Fraction of the PBL (layer K) mass allowed to be used
C               by a cloud-type in time DT
C
C     CP      : Specific heat at constant pressure
C
C     ALHL    : Latent Heat of condensation
C
C     RKAP    : R/Cp, where R is the gas constant
C
C     GRAV    : Acceleration due to gravity
C
C     P00     : A reference pressure in hPa, useually 1000 hPa
C
C     CRTMSF  : Critical value of mass flux above which cloudiness at 
C               the detrainment layer of that cloud-type is assumed.
C               Affects only cloudiness calculation.
C
C     POI     : 2D array of dimension (lng,K) containing potential
C               temperature. Updated but not initialized by RAS.
C
C     QOI     : 2D array of dimension (lng,K) containing specific
C               humidity. Updated but not initialized by RAS.
C
C     UOI     : 3D array of dimension (lng,K,NTRACER) containing tracers
C               Updated but not initialized by RAS.
C
C     PRS     : 2D array of dimension (lng,K+1) containing pressure
C               in hPa at the interfaces of K-layers from top of the 
C               atmosphere to the bottom. Not modified.
C
C     PRJ     : 2D array of dimension (lng,K+1) containing (PRS/P00) **
C               RKAP.  i.e. Exner function at layer edges. Not modified.
C
C     rhfrac  : 1D array of dimension (lng) containing a rel.hum. scaling
C               fraction. Not modified.
C
C  Output:
C  -------
C
C     PCU     : 1D array of length lng containing accumulated
C               precipitation in mm/sec.
C
C     CLN     : 2D array of dimension (lng,K) containing cloudiness
C               Note:  CLN is bumped but NOT initialized
C
C     TCU     : 2D array of dimension (lng,K) containing accumulated 
C               convective heating (K/sec). 
C
C     QCU     : 2D array of dimension (lng,K) containing accumulated 
C               convective drying (kg/kg/sec).
C
C     CMASS   : 2D array of dimension (lng,K) containing the
C               cloud mass flux (kg/sec). Filled from cloud top
C               to base.
C
C  Temporaries:
C
C      ALF, BET, GAM, ETA, PRH, PRI, HOI, HST, QOL, GMH are temporary
C               2D real arrays of dimension of at least (LENC,K) where LENC is 
C               the horizontal dimension over which convection is invoked.
C
C
C    TX1, TX2, TX3, TX4, TX5, TX6, TX7, TX8, TX9, AKM, QS1, CLF, UHT
C    VHT, WLQ  WFN are temporary real arrays of length at least LENC
C
C    IA, I1, and I2 are temporary integer arrays of length LENC 
C
C
C************************************************************************
      implicit none
C Argument List declarations
      integer nn,lng,LENC,K,NLTOP,nlayr,ic,ntracedim, ntracer
      _RL rasalf
      LOGICAL SETRAS
      _RL frac, cp,  alhl, rkap, grav, p00, crtmsf
      _RL POI(lng,K),QOI(lng,K),PRS(lng,K+1),PRJ(lng,K+1)
      _RL uoi(lng,nlayr,ntracedim)
      _RL PCU(LENC), CLN(lng)
      _RL TCU(lng,K),QCU(lng,K),ucu(lng,k,ntracedim),CMASS(lng,K)
      _RL ALF(lng,K), BET(lng,K),  GAM(lng,K), PRH(lng,K), PRI(lng,K)
      _RL HOL(LENC,K), ETA(LENC,K), HST(LENC,K), QOL(LENC,K)
      _RL GMH(LENC,K)
      _RL TX1(LENC), TX2(LENC), TX3(LENC), TX4(LENC)
      _RL TX5(LENC), TX6(LENC), TX7(LENC), TX8(LENC)
      _RL ALM(LENC), WFN(LENC), AKM(LENC), QS1(LENC)
      _RL WLQ(LENC), CLF(LENC)
      _RL uht(lng,ntracedim)
      integer IA(LENC), I1(LENC),I2(LENC)
      _RL      rhfrac(lng)

C Local Variables
      _RL daylen,half,one,zero,cmb2pa,rhmax
      PARAMETER (DAYLEN=86400.0,  HALF=0.5,  ONE=1.0, ZERO=0.0)
      PARAMETER (CMB2PA=100.0)
      PARAMETER (RHMAX=0.9999)
      _RL rkapp1,onebcp,albcp,onebg,cpbg,twobal
C
      integer nt,km1,ic1,i,L,len1,len2,isav,len11,ii
      integer lena,lena1,lenb
      _RL tem,tem1
 
c Explicit Inline Directives
c --------------------------
#ifdef CRAY
#ifdef f77
cfpp$ expand (qsat)
#endif
#endif

      RKAPP1 = 1.0  + RKAP
      ONEBCP = 1.0  / CP
      ALBCP  = ALHL * ONEBCP
      ONEBG  = 1.0  / GRAV
      CPBG   = CP   * ONEBG
      TWOBAL = 2.0 / ALHL

C
      KM1 = K  - 1
      IC1 = IC + 1
C
C      SETTING ALF, BET, GAM, PRH, AND PRI : DONE ONLY WHEN SETRAS=.T.
C

      IF (SETRAS) THEN

         DO 2050 L=1,K
         DO 2030 I=1,LENC
          PRH(I,L) = (PRJ(I,L+1)*PRS(I,L+1) - PRJ(I,L)*PRS(I,L))
     *                            / ((PRS(I,L+1)-PRS(I,L)) * RKAPP1)
 2030    CONTINUE
 2050    CONTINUE

         DO 2070 L=1,K
         DO 2060 I=1,LENC
          TX5(I) = POI(I,L) * PRH(I,L)
          TX1(I) = (PRS(I,L) + PRS(I,L+1)) * 0.5
          TX3(I) = TX5(I)
          CALL QSAT(TX3(I), TX1(I), TX2(I), TX4(I), .TRUE.)
          ALF(I,L) = TX2(I) - TX4(I) * TX5(I)
          BET(I,L) = TX4(I) * PRH(I,L)
          GAM(I,L) = 1.0 / ((1.0 + TX4(I)*ALBCP) * PRH(I,L))
          PRI(I,L) = (CP/CMB2PA) / (PRS(I,L+1) - PRS(I,L))
 2060    CONTINUE
 2070    CONTINUE

      ENDIF
C
C
      DO 10 L=1,K
      DO 10 I=1,lng
      TCU(I,L) = 0.0
      QCU(I,L) = 0.0
      CMASS(I,L) = 0.0
   10 CONTINUE

      do nt = 1,ntracer
      do L=1,K
      do I=1,LENC
      ucu(I,L,nt) = 0.0
      enddo
      enddo
      enddo
C
      DO 30 I=1,LENC
      TX1(I)   = PRJ(I,K+1) * POI(I,K)
      QS1(I)   = ALF(I,K) + BET(I,K)*POI(I,K)
      QOL(I,K) = MIN(QS1(I)*RHMAX,QOI(I,K))

      HOL(I,K) = TX1(I)*CP + QOL(I,K)*ALHL
      ETA(I,K) = ZERO
      TX2(I)   = (PRJ(I,K+1) - PRJ(I,K)) * POI(I,K) * CP
   30 CONTINUE
C
      IF (IC .LT. KM1) THEN
         DO 3703 L=KM1,IC1,-1
         DO 50 I=1,LENC
         QS1(I)   = ALF(I,L) + BET(I,L)*POI(I,L)
         QOL(I,L) = MIN(QS1(I)*RHMAX,QOI(I,L))
C
         TEM1 = TX2(I) + PRJ(I,L+1) * POI(I,L) * CP 
         HOL(I,L) = TEM1 + QOL(I,L )* ALHL
         HST(I,L) = TEM1 + QS1(I)   * ALHL

         TX1(I)   = (PRJ(I,L+1) - PRJ(I,L)) * POI(I,L)
         ETA(I,L) = ETA(I,L+1) + TX1(I)*CPBG
         TX2(I)   = TX2(I)     + TX1(I)*CP
   50    CONTINUE
C
 3703    CONTINUE
      ENDIF


      DO 70 I=1,LENC
      HOL(I,IC) = TX2(I)
      QS1(I)    = ALF(I,IC) + BET(I,IC)*POI(I,IC)
      QOL(I,IC) = MIN(QS1(I)*RHMAX,QOI(I,IC)) 
c
      TEM1      = TX2(I) + PRJ(I,IC1) * POI(I,IC) * CP 
      HOL(I,IC) = TEM1 + QOL(I,IC) * ALHL
      HST(I,IC) = TEM1 + QS1(I)    * ALHL
C
      TX3(I   ) = (PRJ(I,IC1) - PRH(I,IC)) * POI(I,IC)
      ETA(I,IC) = ETA(I,IC1) + CPBG * TX3(I)
   70 CONTINUE
C
      DO 130 I=1,LENC
      TX2(I) = HOL(I,K)  - HST(I,IC)
      TX1(I) = ZERO

  130 CONTINUE
C
C     ENTRAINMENT PARAMETER
C
      DO 160 L=IC,KM1
      DO 160 I=1,LENC
      TX1(I) = TX1(I) + (HST(I,IC) - HOL(I,L)) * (ETA(I,L) - ETA(I,L+1))
  160 CONTINUE
C
      LEN1 = 0
      LEN2 = 0
      ISAV = 0
      DO 195 I=1,LENC
      IF (TX1(I) .GT. ZERO .AND. TX2(I) .GT. ZERO
     .                     .AND. rhfrac(i).ne.0.0 ) THEN
         LEN1      = LEN1 + 1
         IA(LEN1)  = I
         ALM(LEN1) = TX2(I) / TX1(I)
      ENDIF
  195 CONTINUE
C
      LEN2 = LEN1
      if (IC1 .lt. K) then
         DO 196 I=1,LENC
         IF (TX2(I) .LE. 0.0 .AND. (HOL(I,K) .GT. HST(I,IC1))
     .                       .AND. rhfrac(i).ne.0.0 ) THEN
            LEN2      = LEN2 + 1
            IA(LEN2)  = I
            ALM(LEN2) = 0.0
         ENDIF
  196    CONTINUE
      endif
C
      IF (LEN2 .EQ. 0) THEN 
         DO 5010 I=1,LENC*K
         HST(I,1) = 0.0
         QOL(I,1) = 0.0
 5010    CONTINUE
         DO 5020 I=1,LENC
         PCU(I) = 0.0
 5020    CONTINUE
         RETURN
      ENDIF
      LEN11 = LEN1 + 1
C
C     NORMALIZED MASSFLUX
C
      DO 250 I=1,LEN2
      ETA(I,K) = 1.0
      II       = IA(I)
      TX2(I)   = 0.5 * (PRS(II,IC) + PRS(II,IC1))
      TX4(I)   = PRS(II,K)
  250 CONTINUE
C
      DO 252 I=LEN11,LEN2
      WFN(I)   = 0.0
      II       = IA(I)
      IF (HST(II,IC1) .LT. HST(II,IC)) THEN
         TX6(I) = (HST(II,IC1)-HOL(II,K))/(HST(II,IC1)-HST(II,IC))
      ELSE
         TX6(I) = 0.0
      ENDIF
      TX2(I) = 0.5 * (PRS(II,IC1)+PRS(II,IC1+1)) * (1.0-TX6(I))
     *                             + TX2(I)      * TX6(I)
  252 CONTINUE
C
      CALL ACRITN(LEN2, TX2, TX4, TX3)
C
      DO 260 L=KM1,IC,-1
      DO 255 I=1,LEN2
      TX1(I) = ETA(IA(I),L)
  255 CONTINUE
      DO 260 I=1,LEN2
      ETA(I,L) = 1.0 + ALM(I) * TX1(I)
  260 CONTINUE
C
C     CLOUD WORKFUNCTION
C
      IF (LEN1 .GT. 0) THEN
         DO 270 I=1,LEN1
         II = IA(I)
         WFN(I) = - GAM(II,IC) * (PRJ(II,IC1) - PRH(II,IC))
     *                         *  HST(II,IC) * ETA(I,IC1)
  270    CONTINUE
      ENDIF
C
      DO 290 I=1,LEN2
      II = IA(I)
      TX1(I) = HOL(II,K)
  290 CONTINUE
C
      IF (IC1 .LE. KM1) THEN

         DO 380 L=KM1,IC1,-1
         DO 380 I=1,LEN2
         II = IA(I)
         TEM = TX1(I) + (ETA(I,L) - ETA(I,L+1)) * HOL(II,L)
C
         PCU(I) = PRJ(II,L+1) - PRH(II,L)
         TEM1   = ETA(I,L+1) * PCU(I)
         TX1(I) = TX1(I)*PCU(I)
C
         PCU(I) = PRH(II,L) - PRJ(II,L)
         TEM1   = (TEM1 + ETA(I,L) * PCU(I)) * HST(II,L)
         TX1(I) = TX1(I) + TEM*PCU(I)
C
         WFN(I) = WFN(I) + (TX1(I) - TEM1) * GAM(II,L)
         TX1(I) = TEM
  380    CONTINUE
      ENDIF
C
      LENA = 0
      IF (LEN1 .GT. 0) THEN
         DO 512 I=1,LEN1
         II = IA(I)
         WFN(I) = WFN(I) + TX1(I) * GAM(II,IC)*(PRJ(II,IC1)-PRH(II,IC))
     *                   - TX3(I)
         IF (WFN(I) .GT. 0.0) THEN
            LENA = LENA + 1
            I1(LENA) = IA(I)
            I2(LENA) = I
            TX1(LENA) = WFN(I)
            TX2(LENA) = QS1(IA(I))
            TX6(LENA) = 1.0
         ENDIF
  512    CONTINUE
      ENDIF
      LENB = LENA
      DO 515 I=LEN11,LEN2
      WFN(I) = WFN(I) - TX3(I)
      IF (WFN(I) .GT. 0.0 .AND. TX6(I) .GT. 0.0) THEN
         LENB = LENB + 1
         I1(LENB)  = IA(I)
         I2(LENB)  = I
         TX1(LENB) = WFN(I)
         TX2(LENB) = QS1(IA(I))
         TX4(LENB) = TX6(I)
      ENDIF
  515 CONTINUE
C
      IF (LENB .LE. 0) THEN
         DO 5030 I=1,LENC*K
         HST(I,1) = 0.0
         QOL(I,1) = 0.0
 5030    CONTINUE
         DO 5040 I=1,LENC
         PCU(I) = 0.0
 5040    CONTINUE
         RETURN
      ENDIF

C
      DO 516 I=1,LENB
      WFN(I) = TX1(I)
      QS1(I) = TX2(I)
  516 CONTINUE
C
      DO 520 L=IC,K
      DO 517 I=1,LENB
      TX1(I) = ETA(I2(I),L)
  517 CONTINUE
      DO 520 I=1,LENB
      ETA(I,L) = TX1(I)
  520 CONTINUE
C
      LENA1 = LENA + 1
C
      DO 510 I=1,LENA
      II = I1(I)
      TX8(I) = HST(II,IC) - HOL(II,IC)
  510 CONTINUE
      DO 530 I=LENA1,LENB
      II = I1(I)
      TX6(I) = TX4(I)
      TEM    = TX6(I) * (HOL(II,IC)-HOL(II,IC1)) + HOL(II,IC1)
      TX8(I) = HOL(II,K) - TEM

      TEM1   = TX6(I) * (QOL(II,IC)-QOL(II,IC1)) + QOL(II,IC1)
      TX5(I) = TEM    - TEM1 * ALHL
      QS1(I) = TEM1   + TX8(I)*(ONE/ALHL)
      TX3(I) = HOL(II,IC)
  530 CONTINUE
C
C
      DO 620 I=1,LENB
      II = I1(I)
      WLQ(I) = QOL(II,K) - QS1(I)     * ETA(I,IC)
      TX7(I) = HOL(II,K)
  620 CONTINUE
      DO NT=1,Ntracer
      DO 621 I=1,LENB
      II = I1(I)
      UHT(I,NT) = UOI(II,K+nltop-1,NT)-UOI(II,IC+nltop-1,NT) * ETA(I,IC)
  621 CONTINUE
      ENDDO
C
      DO 635 L=KM1,IC,-1
      DO 630 I=1,LENB
      II = I1(I)
      TEM    = ETA(I,L) - ETA(I,L+1)
      WLQ(I) = WLQ(I) + TEM * QOL(II,L)
  630 CONTINUE
  635 CONTINUE
      DO NT=1,Ntracer
      DO L=KM1,IC,-1
      DO I=1,LENB
      II = I1(I)
      TEM    = ETA(I,L) - ETA(I,L+1)
      UHT(I,NT) = UHT(I,NT) + TEM *  UOI(II,L+nltop-1,NT)
      ENDDO
      ENDDO
      ENDDO
C
C     CALCULATE GS AND PART OF AKM (THAT REQUIRES ETA)
C
      DO 690 I=1,LENB
      II = I1(I)
c     TX7(I)     = HOL(II,K)
      TEM        = (POI(II,KM1) - POI(II,K)) / (PRH(II,K) - PRH(II,KM1))
      HOL(I,K)   = TEM * (PRJ(II,K)-PRH(II,KM1))*PRH(II,K)*PRI(II,K)
      HOL(I,KM1) = TEM * (PRH(II,K)-PRJ(II,K))*PRH(II,KM1)*PRI(II,KM1)
      AKM(I)     = ZERO
      TX2(I)     = 0.5 * (PRS(II,IC) + PRS(II,IC1))
  690 CONTINUE

      IF (IC1 .LE. KM1) THEN
         DO 750 L=KM1,IC1,-1
         DO 750 I=1,LENB
         II = I1(I)
         TEM      = (POI(II,L-1) - POI(II,L)) * ETA(I,L)
     *                            / (PRH(II,L) - PRH(II,L-1))
C
         HOL(I,L)   = TEM * (PRJ(II,L)-PRH(II,L-1)) * PRH(II,L)
     *                    *  PRI(II,L)  + HOL(I,L)
         HOL(I,L-1) = TEM * (PRH(II,L)-PRJ(II,L)) * PRH(II,L-1)
     *                                              * PRI(II,L-1)
C
         AKM(I)   = AKM(I) - HOL(I,L)
     *            * (ETA(I,L)   * (PRH(II,L)-PRJ(II,L)) +
     *               ETA(I,L+1) * (PRJ(II,L+1)-PRH(II,L))) / PRH(II,L)
  750    CONTINUE
      ENDIF
C
C
      CALL RNCL(LENB, TX2, TX1, CLF)

      DO 770 I=1,LENB
      TX2(I) = (ONE - TX1(I)) * WLQ(I)
      WLQ(I) = TX1(I) * WLQ(I)
C
      TX1(I) = HOL(I,IC)
  770 CONTINUE
      DO 790 I=LENA1, LENB
      II = I1(I)
      TX1(I) = TX1(I) + (TX5(I)-TX3(I)+QOL(II,IC)*ALHL)*(PRI(II,IC)/CP)
  790 CONTINUE

      DO 800 I=1,LENB
      HOL(I,IC) = TX1(I) - TX2(I) * ALBCP * PRI(I1(I),IC)
  800 CONTINUE

      IF (LENA .GT. 0) THEN
         DO 810 I=1,LENA
         II = I1(I)
         AKM(I) = AKM(I) - ETA(I,IC1) * (PRJ(II,IC1) - PRH(II,IC)) 
     *                                * TX1(I) / PRH(II,IC)
  810    CONTINUE
      ENDIF
c
C     CALCULATE GH
C
      DO 830 I=1,LENB
      II = I1(I)
      TX3(I)   =  QOL(II,KM1) - QOL(II,K)
      GMH(I,K) = HOL(I,K) + TX3(I) * PRI(II,K) * (ALBCP)

      AKM(I)   = AKM(I) + GAM(II,KM1)*(PRJ(II,K)-PRH(II,KM1)) 
     *                               * GMH(I,K)
      TX3(I)   =  zero
  830 CONTINUE
C
      IF (IC1 .LE. KM1) THEN
         DO 840 L=KM1,IC1,-1
         DO 840 I=1,LENB
         II = I1(I)
         TX2(I) = TX3(I)
         TX3(I) = (QOL(II,L-1) - QOL(II,L)) * ETA(I,L)
         TX2(I) = TX2(I) + TX3(I)
C
         GMH(I,L) = HOL(I,L) + TX2(I)   * PRI(II,L) * (ALBCP*HALF)
  840    CONTINUE
C
C
      ENDIF
      DO 850 I=LENA1,LENB
      TX3(I) = TX3(I) + TWOBAL
     *       * (TX7(I) - TX8(I) - TX5(I) - QOL(I1(I),IC)*ALHL)
  850 CONTINUE
      DO 860 I=1,LENB
      GMH(I,IC) = TX1(I) + PRI(I1(I),IC) * ONEBCP
     *          * (TX3(I)*(ALHL*HALF) + ETA(I,IC) * TX8(I))
  860 CONTINUE
C
C     CALCULATE HC PART OF AKM
C
      IF (IC1 .LE. KM1) THEN
         DO 870 I=1,LENB
         TX1(I) = GMH(I,K)
  870    CONTINUE
         DO 3725 L=KM1,IC1,-1
         DO 880 I=1,LENB
         II = I1(I)
         TX1(I) = TX1(I) + (ETA(I,L) - ETA(I,L+1)) * GMH(I,L)
         TX2(I) = GAM(II,L-1) * (PRJ(II,L) - PRH(II,L-1))
  880    CONTINUE
C
         IF (L .EQ. IC1) THEN
            DO 890 I=LENA1,LENB
            TX2(I) = ZERO
  890       CONTINUE
         ENDIF
         DO 900 I=1,LENB
         II = I1(I)
         AKM(I) = AKM(I) + TX1(I) * 
     *          (TX2(I) + GAM(II,L)*(PRH(II,L)-PRJ(II,L)))
  900    CONTINUE
 3725    CONTINUE
      ENDIF
C
      DO 920 I=LENA1,LENB
      II = I1(I)
      TX2(I) = 0.5 * (PRS(II,IC) + PRS(II,IC1))
     *       + 0.5*(PRS(II,IC+2) - PRS(II,IC)) * (ONE-TX6(I))
c
      TX1(I) = PRS(II,IC1)
      TX5(I) = 0.5 * (PRS(II,IC1) + PRS(II,IC+2))
C
      IF ((TX2(I) .GE. TX1(I)) .AND. (TX2(I) .LT. TX5(I))) THEN
         TX6(I)     = ONE - (TX2(I) - TX1(I)) / (TX5(I) - TX1(I))
C
         TEM        = PRI(II,IC1) / PRI(II,IC)
         HOL(I,IC1) = HOL(I,IC1) + HOL(I,IC) * TEM
         HOL(I,IC)  = ZERO
C
         GMH(I,IC1) = GMH(I,IC1) + GMH(I,IC) * TEM
         GMH(I,IC)  = ZERO
      ELSEIF (TX2(I) .LT. TX1(I)) THEN
         TX6(I) = 1.0
      ELSE
         TX6(I) = 0.0
      ENDIF
  920 CONTINUE
C
C
      DO I=1,LENC
      PCU(I) = 0.0
      ENDDO

      DO 970 I=1,LENB
      II = I1(I)
      IF (AKM(I) .LT. ZERO .AND. WLQ(I) .GE. 0.0) THEN
         WFN(I) = - TX6(I) * WFN(I) * RASALF / AKM(I)
      ELSE
         WFN(I) = ZERO
      ENDIF
      TEM       = (PRS(II,K+1)-PRS(II,K))*(CMB2PA*FRAC)
      WFN(I)    = MIN(WFN(I), TEM)
C
C     compute cloud amount
C
CC    TX1(I) = CLN(II)
CC    IF (WFN(I) .GT. CRTMSF)  TX1(I) = TX1(I) + CLF(I)
CC    IF (TX1(I) .GT. ONE)  TX1(I) = ONE
C
C     PRECIPITATION
C
      PCU(II) =  WLQ(I) * WFN(I) * ONEBG
C
C     CUMULUS FRICTION AT THE BOTTOM LAYER
C
      TX4(I)   = WFN(I) * (1.0/ALHL)
      TX5(I)   = WFN(I) * ONEBCP
  970 CONTINUE
C
C     compute cloud mass flux for diagnostic output
C
      DO L = IC,K
      DO I=1,LENB
       II = I1(I)
       if(L.lt.K)then
       CMASS(II,L) =  ETA(I,L+1) * WFN(I) * ONEBG
       else
       CMASS(II,L) =  WFN(I) * ONEBG
       endif
      ENDDO
      ENDDO
C
CC    DO 975 I=1,LENB
CC    II = I1(I)
CC    CLN(II) = TX1(I)
CC975 CONTINUE
C
C     THETA AND Q CHANGE DUE TO CLOUD TYPE IC
C
 
c     TEMA = 0.0
c     TEMB = 0.0
      DO 990 L=IC,K
      DO 980 I=1,LENB
      II = I1(I)
      TEM      = (GMH(I,L) - HOL(I,L)) * TX4(I)
      TEM1     =  HOL(I,L) * TX5(I)
C
      TCU(II,L) = TEM1 / PRH(II,L)
      QCU(II,L) = TEM
  980 CONTINUE

c     I = I1(IP1)
c
c     TEM  = (PRS(I,L+1)-PRS(I,L)) * (ONEBG*100.0)
c     TEMA = TEMA +  TCU(I,L) * PRH(I,L) * TEM * (CP/ALHL)
c     TEMB = TEMB +  QCU(I,L)            * TEM
C
  990 CONTINUE
C
c Compute Tracer Tendencies
c -------------------------
      do nt = 1,ntracer
c
c Tracer Tendency at the Bottom Layer
c -----------------------------------
      DO 995 I=1,LENB
      II = I1(I)
      TEM    = half*TX5(I) * PRI(II,K)
      TX1(I) = ( UOI(II,KM1+nltop-1,nt) - UOI(II,K+nltop-1,nt))
      ucu(II,K,nt) = TEM * TX1(I)
  995 CONTINUE
c
c Tracer Tendency at all other Levels
c -----------------------------------
      DO 1020 L=KM1,IC1,-1
      DO 1010 I=1,LENB
      II = I1(I)
      TEM = half*TX5(I) * PRI(II,L)
      TEM1   = TX1(I)
      TX1(I) = (UOI(II,L-1+nltop-1,nt)-UOI(II,L+nltop-1,nt)) * ETA(I,L)
      TX3(I) = (TX1(I) + TEM1) * TEM
 1010 CONTINUE
      DO 1020 I=1,LENB
      II = I1(I)
      ucu(II,L,nt) = TX3(I)
 1020 CONTINUE

      DO 1030 I=1,LENB
      II = I1(I)
      IF (TX6(I) .GE. 1.0) THEN
         TEM    = half*TX5(I) * PRI(II,IC)
      ELSE
         TEM = 0.0
      ENDIF
      TX1(I) = (TX1(I) + UHT(I,nt) + UHT(I,nt)) * TEM
 1030 CONTINUE
      DO 1040 I=1,LENB
      II = I1(I)
      ucu(II,IC,nt) = TX1(I)
 1040 CONTINUE
 
      enddo
C
C     PENETRATIVE CONVECTION CALCULATION OVER
C

      RETURN
      END
      SUBROUTINE RNCL(lng, PL, RNO, CLF)
C
C*********************************************************************
C********************** Relaxed Arakawa-Schubert *********************
C************************   SUBROUTINE  RNCL  ************************
C**************************** 23 July 1992 ***************************
C*********************************************************************
      implicit none
C Argument List declarations
      integer lng
      _RL PL(lng),  RNO(lng), CLF(lng)

C Local Variables
      _RL p5,p8,pt8,pt2,pfac,p4,p6,p7,p9,cucld,cfac
      PARAMETER (P5=500.0,  P8=800.0, PT8=0.8, PT2=0.2)
      PARAMETER (PFAC=PT2/(P8-P5))
      PARAMETER (P4=400.0,    P6=401.0)
      PARAMETER (P7=700.0,    P9=900.0)
      PARAMETER (CUCLD=0.5,CFAC=CUCLD/(P6-P4))

      integer i
C
      DO 10 I=1,lng
                           rno(i) = 1.0
ccc   if( pl(i).le.400.0 ) rno(i) = max( 0.75, 1.0-0.0025*(400.0-pl(i)) )

ccc   IF ( PL(I).GE.P7 .AND. PL(I).LE.P9 ) THEN
ccc     RNO(I) = ((P9-PL(I))/(P9-P7)) **2
ccc   ELSE IF (PL(I).GT.P9) THEN
ccc     RNO(I) = 0.
ccc   ENDIF

      CLF(I) = CUCLD
C
CARIESIF (PL(I) .GE. P5 .AND. PL(I) .LE. P8) THEN
CARIES    RNO(I) = (P8-PL(I))*PFAC + PT8
CARIESELSEIF (PL(I) .GT. P8 ) THEN
CARIES    RNO(I) = PT8
CARIESENDIF
CARIES
      IF (PL(I) .GE. P4 .AND. PL(I) .LE. P6) THEN
         CLF(I) = (P6-PL(I))*CFAC
      ELSEIF (PL(I) .GT. P6 ) THEN
         CLF(I) = 0.0
      ENDIF
   10 CONTINUE
C
      RETURN
      END
      SUBROUTINE ACRITN ( lng,PL,PLB,ACR )
 
C*********************************************************************
C********************** Relaxed Arakawa-Schubert *********************
C************************** SUBROUTINE  ACRIT    *********************
C******************  modified August 28, 1996   L.Takacs  ************
C****                                                            *****
C****  Note:  Data obtained from January Mean After-Analysis     *****
C****         from 4x5 46-layer GEOS Assimilation                *****
C****                                                            *****
C*********************************************************************
      implicit none
C Argument List declarations
      integer lng
      _RL PL(lng), PLB(lng), ACR(lng)

C Local variables
      integer lma
      parameter  (lma=18)
      _RL p(lma)
      _RL a(lma)
      integer i,L
      _RL temp

      data p  / 93.81, 111.65, 133.46, 157.80, 186.51,
     .         219.88, 257.40, 301.21, 352.49, 409.76,
     .         471.59, 535.04, 603.33, 672.79, 741.12,
     .         812.52, 875.31, 930.20/

      data a  / 3.35848, 3.13645, 2.48072, 2.08277, 1.53364,
     .          1.01971,  .65846,  .45867,  .38687,  .31002,
     .           .25574,  .20347,  .17254,  .15260,  .16756,
     .           .09916,  .10360,  .05880/

 
      do L=1,lma-1
      do i=1,lng
         if( pl(i).ge.p(L)   .and.
     .       pl(i).le.p(L+1)) then
         temp = ( pl(i)-p(L) )/( p(L+1)-p(L) )
         acr(i) = a(L+1)*temp + a(L)*(1-temp)
         endif
      enddo
      enddo

      do i=1,lng
      if( pl(i).lt.p(1)   ) acr(i) = a(1)
      if( pl(i).gt.p(lma) ) acr(i) = a(lma)
      enddo

      do i=1,lng
      acr(i) = acr(i) * (plb(i)-pl(i))
      enddo

      RETURN
      END
       SUBROUTINE RNEVP(NN,IRUN,NLAY,TL,QL,RAIN,PL,CLFRAC,SP,DP,PLKE,
     1   PLK,TH,TEMP1,TEMP2,TEMP3,ITMP1,ITMP2,RCON,RLAR,CLSBTH,tmscl,
     2   tmfrc,cp,gravity,alhl,gamfac,cldlz,RHCRIT,offset,alpha)
 
      implicit none
C Argument List declarations
      integer nn,irun,nlay
      _RL TL(IRUN,NLAY),QL(IRUN,NLAY),RAIN(IRUN,NLAY),
     . PL(IRUN,NLAY),CLFRAC(IRUN,NLAY),SP(IRUN),TEMP1(IRUN,NLAY),
     . TEMP2(IRUN,NLAY),PLKE(IRUN,NLAY+1),
     . RCON(IRUN),RLAR(IRUN),DP(IRUN,NLAY),PLK(IRUN,NLAY),TH(IRUN,NLAY),
     . TEMP3(IRUN,NLAY)
      integer ITMP1(IRUN,NLAY),ITMP2(IRUN,NLAY)
      _RL CLSBTH(IRUN,NLAY)
      _RL tmscl,tmfrc,cp,gravity,alhl,gamfac,offset,alpha
      _RL cldlz(irun,nlay)
      _RL rhcrit(irun,nlay)
C
C Local Variables
      _RL zm1p04,zero,two89,zp44,zp01,half,zp578,one,thousand,z3600
      _RL zp1,zp001
      PARAMETER (ZM1P04 = -1.04E-4 )
      PARAMETER (ZERO = 0.)
      PARAMETER (TWO89= 2.89E-5)
      PARAMETER ( ZP44= 0.44)
      PARAMETER ( ZP01= 0.01)
      PARAMETER ( ZP1 = 0.1 )
      PARAMETER ( ZP001= 0.001)
      PARAMETER ( HALF= 0.5)
      PARAMETER ( ZP578 = 0.578 )
      PARAMETER ( ONE = 1.)
      PARAMETER ( THOUSAND = 1000.)
      PARAMETER ( Z3600 = 3600.)
C
      _RL EVP9(IRUN,NLAY)
      _RL water(irun),crystal(irun)
      _RL watevap(irun),iceevap(irun)
      _RL fracwat,fracice, tice,rh,fact,dum
      _RL rainmax(irun)
      _RL getcon,rphf,elocp,cpog,relax
      _RL exparg,arearat,rpow

      integer i,L,n,nlaym1,irnlay,irnlm1

c Explicit Inline Directives
c --------------------------
#ifdef CRAY
#ifdef f77
cfpp$ expand (qsat)
#endif
#endif

       tice = getcon('FREEZING-POINT')

       fracwat = 0.70
       fracice = 0.01

       NLAYM1 = NLAY - 1
       IRNLAY = IRUN*NLAY
       IRNLM1 = IRUN*(NLAY-1)
 
       RPHF = Z3600/tmscl
 
       ELOCP   = alhl/cp
       CPOG    = cp/gravity

       DO I = 1,IRUN
          RLAR(I) = 0.
         water(i) = 0.
       crystal(i) = 0.
       ENDDO

       do L = 1,nlay
       do i = 1,irun
          EVP9(i,L) = 0.
         TEMP1(i,L) = 0.
         TEMP2(i,L) = 0.
         TEMP3(i,L) = 0.
        CLSBTH(i,L) = 0.
         cldlz(i,L) = 0.
       enddo
       enddo
 
C RHO(ZERO) / RHO FOR TERMINAL VELOCITY APPROX.
c ---------------------------------------------
      DO L = 1,NLAY
      DO I = 1,IRUN
        TEMP2(I,L) = PL(I,L)*ZP001
        TEMP2(I,L) = SQRT(TEMP2(I,L))
      ENDDO
      ENDDO

C INVERSE OF MASS IN EACH LAYER
c -----------------------------
      DO L = 1,NLAY
      DO I = 1,IRUN
      TEMP3(I,L) = GRAVITY*ZP01 / DP(I,L)
      ENDDO
      ENDDO
 
C DO LOOP FOR MOISTURE EVAPORATION ABILITY AND CONVEC EVAPORATION.
c ----------------------------------------------------------------
      DO 100 L=1,NLAY

      DO I = 1,IRUN
        TEMP1(I,3) = TL(I,L)
        TEMP1(I,4) = QL(I,L)
      ENDDO

      DO 50 N=1,2
       IF(N.EQ.1)RELAX=HALF
       IF(N.GT.1)RELAX=ONE

      DO I = 1,IRUN
       call qsat ( temp1(i,3),pl(i,L),temp1(i,2),temp1(i,6),.true. )
       TEMP1(I,5)=TEMP1(I,2)-TEMP1(I,4)
       TEMP1(I,6)=TEMP1(I,6)*ELOCP
       TEMP1(I,5)=TEMP1(I,5)/(ONE+TEMP1(I,6))
       TEMP1(I,4)=TEMP1(I,4)+TEMP1(I,5)*RELAX
       TEMP1(I,3)=TEMP1(I,3)-TEMP1(I,5)*ELOCP*RELAX
      ENDDO
  50  CONTINUE

      DO I = 1,IRUN
          EVP9(I,L) = (TEMP1(I,4) -  QL(I,L))/TEMP3(I,L)
C convective detrained water
         cldlz(i,L) = rain(i,L)*temp3(i,L)  
       if(  tl(i,L).gt.tice-20.) then
           water(i) =   water(i) + rain(i,L)
       else
         crystal(i) = crystal(i) + rain(i,L)
       endif
      ENDDO
 
C**********************************************************************
C  FOR CONVECTIVE PRECIP, FIND THE "EVAPORATION EFFICIENCY" USING     *
C                   KESSLERS PARAMETERIZATION                         *
C**********************************************************************
 
      DO 20 I=1,IRUN

        iceevap(i) = 0.
        watevap(i) = 0.

       if( (evp9(i,L).gt.0.) .and. (crystal(i).gt.0.) ) then
       iceevap(I) = EVP9(I,L)*fracice
       IF(iceevap(i).GE.crystal(i)) iceevap(i) = crystal(i)
       EVP9(I,L)=EVP9(I,L)-iceevap(I)
       crystal(i) = crystal(i) - iceevap(i)
       endif

C and now warm precipitate
       if( (evp9(i,L).gt.0.) .and. (water(i).gt.0.) ) then
       exparg    = ZM1P04*tmscl*((water(i)*RPHF*TEMP2(I,L))**ZP578)
       AREARAT = ONE-(EXP(EXPARG))
       watevap(I) = EVP9(I,L)*AREARAT*fracwat
       IF(watevap(I).GE.water(i)) watevap(I) = water(i)
       EVP9(I,L)=EVP9(I,L)-watevap(I)
       water(i) = water(i) - watevap(i)
       endif

       QL(I,L) = QL(I,L)+(iceevap(i)+watevap(i))*TEMP3(I,L)
       TL(I,L) = TL(I,L)-(iceevap(i)+watevap(i))*TEMP3(I,L)*ELOCP

  20  CONTINUE

  100 CONTINUE

      do i = 1,irun
      rcon(i) = water(i) + crystal(i)
      enddo

C**********************************************************************
C Large Scale Precip
C**********************************************************************

      DO 200 L=1,NLAY
      DO I = 1,IRUN
        rainmax(i) = rhcrit(i,L)*evp9(i,L) +
     .               ql(i,L)*(rhcrit(i,L)-1.)/temp3(i,L)

         if (rainmax(i).LE.0.0) then
           call qsat( tl(i,L),pl(i,L),rh,dum,.false.)
           rh = ql(i,L)/rh

           if( rhcrit(i,L).eq.1.0 ) then
           fact = 1.0
           else
           fact = min( 1.0, alpha + (1.0-alpha)*( rh-rhcrit(i,L)) /
     1                                          (1.0-rhcrit(i,L)) )
           endif

C Do not allow clouds above 10 mb
           if( pl(i,L).ge.10.0 ) CLSBTH(I,L) = fact  
           RLAR(I) = RLAR(I)-rainmax(I)
           QL(I,L) = QL(I,L)+rainmax(I)*TEMP3(I,L)
           TL(I,L) = TL(I,L)-rainmax(I)*TEMP3(I,L)*ELOCP
C Large-scale water
        cldlz(i,L) = cldlz(i,L) - rainmax(i)*temp3(i,L)  
         ENDIF
      ENDDO

      DO I=1,IRUN
         IF((RLAR(I).GT.0.0).AND.(rainmax(I).GT.0.0))THEN
          RPOW=(RLAR(I)*RPHF*TEMP2(I,L))**ZP578
          EXPARG  = ZM1P04*tmscl*RPOW
          AREARAT = ONE-(EXP(EXPARG))
          TEMP1(I,7) = rainmax(I)*AREARAT
          IF(TEMP1(I,7).GE.RLAR(I)) TEMP1(I,7) = RLAR(I)
          RLAR(I) =    RLAR(I)-TEMP1(I,7)
          QL(I,L) =    QL(I,L)+TEMP1(I,7)*TEMP3(I,L)
          TL(I,L) =    TL(I,L)-TEMP1(I,7)*TEMP3(I,L)*ELOCP
         ENDIF
      ENDDO

  200 CONTINUE

      RETURN
      END
      subroutine srclouds (th,q,plk,pl,plke,cloud,cldwat,irun,irise,
     1                                                rhc,offset,alpha)
C***********************************************************************
C
C  PURPOSE:
C  ========
C    Compute non-precipitating cloud fractions
C    based on Slingo and Ritter (1985).
C    Remove cloudiness where conditionally unstable.
C
C  INPUT:
C  ======
C    th ......... Potential Temperature (irun,irise)
C    q .......... Specific  Humidity    (irun,irise)
C    plk ........ P**Kappa at mid-layer (irun,irise)
C    pl ......... Pressure at mid-layer (irun,irise)
C    plke ....... P**Kappa at edge      (irun,irise+1)
C    irun ....... Horizontal   dimension
C    irise ...... Vertical     dimension
C
C  OUTPUT:
C  =======
C    cloud ...... Cloud Fraction        (irun,irise)
C
C***********************************************************************

      implicit none
      integer  irun,irise

      _RL   th(irun,irise)
      _RL    q(irun,irise)
      _RL  plk(irun,irise)
      _RL   pl(irun,irise)
      _RL plke(irun,irise+1)

      _RL  cloud(irun,irise)
      _RL cldwat(irun,irise)
      _RL     qs(irun,irise)

      _RL cp, alhl, getcon, akap
      _RL ratio, temp, elocp
      _RL rhcrit,rh,dum
      integer i,L

      _RL rhc(irun,irise)
      _RL offset,alpha

c Explicit Inline Directives
c --------------------------
#ifdef CRAY
#ifdef f77
cfpp$ expand (qsat)
#endif
#endif

      cp     = getcon('CP')
      alhl   = getcon('LATENT HEAT COND')
      elocp  = alhl/cp
      akap   = getcon('KAPPA')

      do L = 1,irise
      do i = 1,irun
                  temp = th(i,L)*plk(i,L)
      call qsat ( temp,pl(i,L),qs(i,L),dum,.false. )
      enddo
      enddo

      do L  = 2,irise
      do i  = 1,irun
      rh    = q(i,L)/qs(i,L)

      rhcrit = rhc(i,L) - offset
      ratio  = alpha*(rh-rhcrit)/offset
 
      if(cloud(i,L).eq.  0.0 .and. ratio.gt.0.0 ) then
         cloud(i,L) = min( ratio,1.0 )
      endif
 
      enddo
      enddo
 
c Reduce clouds from conditionally unstable layer
c -----------------------------------------------
      call ctei ( th,q,cloud,cldwat,pl,plk,plke,irun,irise )

      return
      end

      subroutine ctei ( th,q,cldfrc,cldwat,pl,plk,plke,im,lm )
      implicit none
      integer im,lm
      _RL  th(im,lm),q(im,lm),plke(im,lm+1),cldwat(im,lm)
      _RL plk(im,lm),pl(im,lm),cldfrc(im,lm)
      integer i,L
      _RL    getcon,cp,alhl,elocp,cpoel,t,p,s,qs,dqsdt,dq
      _RL    k,krd,kmm,f

      cp     = getcon('CP')
      alhl   = getcon('LATENT HEAT COND')
      cpoel  = cp/alhl
      elocp  = alhl/cp

      do L=lm,2,-1
      do i=1,im
          dq = q(i,L)+cldwat(i,L)-q(i,L-1)-cldwat(i,L-1)
      if( dq.eq.0.0 ) dq = 1.0e-20
      k = 1.0 + cpoel*plke(i,L)*( th(i,L)-th(i,L-1) ) / dq

      t = th(i,L)*plk(i,L)
      p = pl(i,L)
      call qsat ( t,p,qs,dqsdt,.true. )
   
      krd = ( cpoel*t*(1+elocp*dqsdt) )/( 1 + 1.608*dqsdt*t )

      kmm = ( 1+elocp*dqsdt )*( 1 + 0.392*cpoel*t )
     .    / ( 2+(1+1.608*cpoel*t)*elocp*dqsdt )

      s = ( (k-krd)/(kmm-krd) )
      f = 1.0 - min( 1.0, max(0.0,1.0-exp(-s)) )

      cldfrc(i,L) = cldfrc(i,L)*f
      cldwat(i,L) = cldwat(i,L)*f

      enddo
      enddo

      return
      end

      subroutine back2grd(gathered,indeces,scattered,irun)
      implicit none
      integer i,irun,indeces(irun)
      _RL gathered(irun),scattered(irun)
      _RL temp(irun)
      do i = 1,irun
       temp(indeces(i)) = gathered(i)
      enddo
      do i = 1,irun
       scattered(i) = temp(i)
      enddo
      return
      end