SUBROUTINE gscond
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  ! Code converted using TO_F90 by Alan Miller
  ! Date: 2000-02-17  Time: 08:34:14
  ! SUBPROGRAM DOCUMENTATION BLOCK
  !                .      .    .
  ! SUBPROGRAM:    GSCOND      GRID SCALE CONDENSATION AND EVAPORATION
  !   PRGRMMR: ZHAO            ORG: W/NP2      DATE: ??-??-??
  !
  ! ABSTRACT:
  !     GSCOND COMPUTES THE GRID SCALE EVAPORATION AND CONDENSATION
  !
  ! PROGRAM HISTORY LOG:
  !   94-??-??  ZHAO       - ORIGINATOR
  !   95-03-25  BLACK      - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL
  !   95-03-28  BLACK      - ADDED EXTERNAL EDGE
  !
  ! USAGE: CALL GSCOND FROM MAIN PROGRAM EBU
  !
  !   INPUT ARGUMENT LIST:
  !       NONE
  !
  !   OUTPUT ARGUMENT LIST:
  !     NONE
  !
  !   OUTPUT FILES:
  !     NONE
  !
  !   SUBPROGRAMS CALLED:
  !
  !     UNIQUE: NONE
  !
  !     LIBRARY: NONE
  !
  !   COMMON BLOCKS: CTLBLK
  !                  LOOPS
  !                  MASKS
  !                  PHYS
  !                  SCRCH
  !                  VRBLS
  !                  CLDWTR
  !                  TEMPV
  !                  PVRBLS
  !                  ACMCLH
  !
  IMPLICIT NONE
  REAL, PARAMETER	:: h1=1.e0,h2=2.e0  &
       , d00=0.e0,d125=.125E0,d5=0.5E0  &
       , a1=610.78E0,a2=17.2693882E0,a3=273.16E0,a4=35.86E0  &
       , pq0=379.90516E0,tresh=.95E0  &
       , cp=1004.6E0,elwv=2.50E6,eliv=2.834E6,row=1.e3,g=9.8E0  &
       , epsq=2.e-12,dldt=2274.e0,tm10=263.16E0,r=287.0E0 , cpr=cp*r,rcpr=h1/(cpr)
  REAL, PARAMETER  :: arcp=a2*(a3-a4)/cp,rcp=h1/cp,pq0c=pq0*tresh,rrog=h1/(row*g)
  !----------------------------------------------------------------------
  INCLUDE "parmeta.f90"
  INCLUDE "parm.tbl.f90"
  !----------------------------------------------------------------------
  INTEGER, PARAMETER	:: la=13
  INTEGER, PARAMETER	:: imjm=im*jm-jm/2,jam=6+2*(jm-10)  &
       , lp1=lm+1,ltop=1,lbot=lm , lscrch=4*lm+1+la +1  &
       , l2=la+lm+1,l4=la+3*lm+1
  !
  INCLUDE "ctlblk.f90"
  INCLUDE "loops.f90"
  INCLUDE "masks.f90"
  INCLUDE "phys.f90"
  INCLUDE "vrbls.f90"
  INCLUDE "cldwtr.f90"
  INCLUDE "tempv.f90"
  INCLUDE "pvrbls.f90"
  INCLUDE "acmclh.f90"
  !
  REAL, DIMENSION(im,jm)	:: pdsl,qi,qint,ccr,e0,cond
  REAL, DIMENSION(im,jm,lm)	:: iw
  !----------------------------------------------------------------------
  REAL	:: mr,ke
  !
  !
  !--------------PREPARATORY CALCULATIONS---------------------------------
  REAL   :: dtph,rdtph,twodt,rtwodt,c0,c1,c2,us,eps,cclimit,climit
  INTEGER  :: i,j,l,lml,iwkl
  REAL   :: ctim1,ctim2,ctim,utim,hh,tmt0,tmt15,ai,bi,qw,u00kl
  REAL   :: fi,thh,pp,pp0,at,aq,ap,tkl,qkl,cwmkl,fiw,elv,qc
  REAL   :: rqkl,rqkll,ccrkl,ec,us00,ccrkl1,aa,ab,ac,ad,ae
  REAL   :: af,ag,condk,qtemp,rqtmp,cone0
  !
  dtph =nphs*dt
  rdtph=1./dtph
  twodt=dtph
  rtwodt=h1/twodt
  c0=1.5E-4
  c1=300.
  c2=0.5
  mr=3.0E-4
  ke=2.0E-5
  us=h1
  eps=0.622E0
  cclimit=1.0E-3
  climit =1.0E-20
  !------------------INITIALIZE T0, Q0 & P0 FOR THE FIRST TIME STEP-------
  IF(ntsd == 1)THEN
     DO  j=1,jm
        DO  i=1,im
           p0(i,j)=pd(i,j)
        END DO
     END DO
     DO  l=1,lm
        DO  j=1,jm
           DO  i=1,im
              t0(i,j,l)=t(i,j,l)
              q0(i,j,l)=q(i,j,l)
           END DO
        END DO
     END DO
  END IF
  !------------------PADDING SPECIFIC HUMIDITY & CWM IF TOO SMALL---------
  DO   l=1,lm
     DO  j=1,jm
        DO  i=1,im
           IF(q(i,j,l) < epsq)q(i,j,l)=epsq*htm(i,j,l)
           IF(cwm(i,j,l) < climit)cwm(i,j,l)=climit*htm(i,j,l)
        END DO
     END DO
  END DO
  DO  j=3,jm-2
     DO  i=1,im
        pdsl(i,j)=res(i,j)*pd(i,j)
        iw(i,j,1)=0
     END DO
  END DO
  !----------------------------COMPUTE UTIM-------------------------------
  ctim1=0.
  ctim2=24.*3600.
  ctim=ntsd*dt
  IF(ctim < ctim1)THEN
     utim=0.
  ELSE
     IF(ctim <= ctim2)THEN
        utim=(ctim-ctim1)/(ctim2-ctim1)
     ELSE
        utim=1.
     END IF
  END IF

  !*************BEGINING OF  GRID-SCALE CONDENSATION/EVAP. LOOP***********

  DO  l=2,lm
     !
     DO  j=3,jm-2
        DO  i=1,im
           cond(i,j)=d00
           e0(i,j)=d00
        END DO
     END DO
     !------------------QW, QI AND QINT--------------------------------------
     DO  j=3,jm-2
        DO  i=1,im
           lml=lm-lmh(i,j)
           hh=htm(i,j,l)*hbm2(i,j)
           tmt0=(t(i,j,l)-273.16)*hh
           tmt15=AMIN1(tmt0,-15.)*hh
           ai=0.008855
           bi=1.
           IF(tmt0 < -20.)THEN
              ai=0.007225
              bi=0.9674
           END IF
           qw=hh*pq0/(pdsl(i,j)*aeta(l)+pt) *EXP(hh*a2*(t(i,j,l)-a3)/(t(i,j,l)-a4))
           qi(i,j)=qw*(bi+ai*AMIN1(tmt0,0.))
           qint(i,j)=qw*(1.-0.00032*tmt15*(tmt15+15.))
           IF(tmt0 <= -40.)qint(i,j)=qi(i,j)
           !-------------------ICE-WATER ID NUMBER IW------------------------------
           IF(tmt0 < -15.)THEN
              u00kl=u00(i,j)+ul(l+lml)*(0.95-u00(i,j))*utim
              fi=q(i,j,l)-u00kl*qi(i,j)
              IF(fi > d00.OR.cwm(i,j,l) > climit)THEN
                 iw(i,j,l)=1
              ELSE
                 iw(i,j,l)=0
              END IF
           END IF

           IF(tmt0 >= 0.)THEN
              iw(i,j,l)=0
           END IF

           IF(tmt0 < 0.0.AND.tmt0 >= -15.)THEN
              iw(i,j,l)=0
              IF(iw(i,j,l-1) == 1.AND.cwm(i,j,l) > climit)iw(i,j,l)=1
           END IF
        END DO
     END DO
     !--------------CONDENSATION AND EVAPORATION OF CLOUD--------------------
     DO  j=3,jm-2
        DO  i=1,im
           !------------------------AT, AQ AND DP/DT-------------------------------
           lml=lm-lmh(i,j)
           hh=htm(i,j,l)*hbm2(i,j)
           thh=twodt*hh
           pp=pdsl(i,j)*aeta(l)+pt
           pp0=p0(i,j)*res(i,j)*aeta(l)+pt
           at=(t(i,j,l)-t0(i,j,l))*rtwodt
           aq=(q(i,j,l)-q0(i,j,l))*rtwodt
           ap=(pp-pp0)*rtwodt
           iwkl=iw(i,j,l)
           u00kl=u00(i,j)+ul(l+lml)*(0.95-u00(i,j))*utim
           tkl=t(i,j,l)
           qkl=q(i,j,l)
           cwmkl=cwm(i,j,l)
           !----------------THE SATUATION SPECIFIC HUMIDITY------------------------
           fiw=FLOAT(iwkl)
           elv=(h1-fiw)*elwv+fiw*eliv
           qc =(h1-fiw)*qint(i,j)+fiw*qi(i,j)
           !----------------THE RELATIVE HUMIDITY----------------------------------
           IF(qc <= d00)THEN
              rqkl=d00
           ELSE
              rqkl=q(i,j,l)/qc
           END IF
           !----------------CLOUD COVER RATIO CCR----------------------------------
           IF(rqkl <= u00kl)THEN
              ccr(i,j)=d00
           ELSE
              rqkll=AMIN1(us,rqkl)
              ccr(i,j)=h1-SQRT((us-rqkll)/(us-u00kl))
           END IF
           !-----------CORRECT CCR IF IT IS TOO SMALL IN LARGE CWM REGIONS--------
           IF(ccr(i,j) >= 0.01.AND.ccr(i,j) <= 0.2.AND.cwm(i,j,l) >= 0.2E-3)THEN
              ccr(i,j)=AMIN1(1.,cwm(i,j,l)*1.e3)
           END IF
           !
           ccrkl=ccr(i,j)
           !-------GIVE UP THIS POINT  IF NO CLOUD NOR CONDENSATION EXIST---------
           IF(ccrkl <= cclimit.AND.cwmkl <= climit)CYCLE
           !----------------EVAPORATION OF CLOUD WATER-----------------------------
           ec=d00
           IF(ccrkl <= cclimit.AND.cwmkl > climit)THEN
              ec=qc*(u00kl-rqkl)*rtwodt
              e0(i,j)=AMAX1(ec,0.0)
              e0(i,j)=AMIN1(cwmkl*rtwodt,e0(i,j))*hh
              e0(i,j)=AMAX1(0.,e0(i,j))
           END IF
           !----------------CONDENSATION OF CLOUD----------------------------------
           !      IF(ccrkl <= 0.20.OR.qc <= epsq)THEN
           !        cond(i,j)=0.
           !        GO TO 80
           !      END IF
           !-----------THE EQS. FOR COND. HAS BEEN REORGANIZED TO REDUCE CPU------
           IF(ccrkl > 0.20.AND.qc > epsq)THEN
              us00=us-u00kl
              ccrkl1=1.-ccrkl
              aa=eps*elv*pp*qkl
              ab=ccrkl*ccrkl1*qc*us00
              ac=ab+0.5*cwmkl
              ad=ab*ccrkl1
              ae=cpr*tkl*tkl
              af=ae*pp
              ag=aa*elv
              ai=cp*aa
              cond(i,j)=(ac-ad)*(af*aq-ai*at+ae*qkl*ap)/(ac*(af+ag))
              !-----------CHECK & CORRECT IF OVER CONDENSATION OCCURS-----------------
              condk=(qkl-u00kl*qc*0.1)*rtwodt
              !     CONDK=(QKL-U00KL*QC*0.6)*RTWODT
              IF(cond(i,j) > condk)cond(i,j)=condk
              !---------CHECK & CORRECT IF SUPERSATUATION IS TOO HIGH----------------
              qtemp=qkl-AMAX1(0.,(cond(i,j)-e0(i,j)))*thh
              rqtmp=qtemp/qc
              IF(rqtmp >= 1.10)cond(i,j)=(qkl-1.10*qc)*rtwodt
              !
              IF(cond(i,j) < d00)cond(i,j)=d00
           ELSE
              cond(i,j)=0.
           END IF
           !-------------------UPDATE OF T, Q AND CWM------------------------------
           !   80 CONTINUE
           cone0=cond(i,j)-e0(i,j)
           cwm(i,j,l)=cone0*thh+cwmkl
           !
           !     ACCUMULATE LATENT HEATING DUE TO GRID-SCALE PRECIP/EVAP.
           !     SCALE BY THE RECIPROCAL OF THE PERIOD AT WHICH THIS ROUTINE
           !     IS CALLED.  THIS PERIOD IS THE PHYSICS TIMESTEP.
           !     
           t(i,j,l)=elv*rcp*cone0*thh+tkl
           train(i,j,l)=elv*rcp*cone0*thh*rdtph+train(i,j,l)
           q(i,j,l)=-cone0*thh+qkl
           IF(cwm(i,j,l) <= 0.)cwm(i,j,l)=0.
        END DO
     END DO
  END DO
  !
  !-------------------SAVE T, Q AND P FOR THIS STEP-----------------------
  DO  l=1,lm
     DO  j=1,jm
        DO  i=1,im
           q0(i,j,l)=q(i,j,l)
           t0(i,j,l)=t(i,j,l)
        END DO
     END DO
  END DO
  DO  j=1,jm
     DO  i=1,im
        p0(i,j)=pd(i,j)
     END DO
  END DO
END SUBROUTINE gscond
