SUBROUTINE precpd
  !
  ! $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:58:35
  !
  ! subprogram documentation block
  !                .      .    .
  ! subprogram:    precpd      large scale precipitation
  !   prgrmmr: zhao            org: w/np2      date: ??-??-??
  !
  ! abstract:
  !     precpd computes the grid scale precipitation.
  !
  ! program history log:
  !   94-??-??  zhao       - originator
  !   95-03-25  black      - conversion from 1-d to 2-d in horizontal
  !   95-11-20  abeles     - parallel optimization
  !   96-03-29  black      - removed scrch common
  !   96-07-18  zhao       - new wmin calculation
  !   96-09-25  baldwin    - new sr calculation
  !
  ! usage: call precpd 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
  !                  vrbls
  !                  cldwtr
  !                  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,r=287.0e0,c0=0.15  &
       , 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,eliw=eliv-elwv
  REAL, PARAMETER	:: arcp=a2*(a3-a4)/cp,rcp=h1/cp, &
       pq0c=pq0*tresh,rrog=h1/(row*g) , rrow=h1/row
  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  &
       , l1=la+lm+1,l2=la+2*lm+1 , l3=la+3*lm+1,l4=la+4*lm+1
  INCLUDE "ctlblk.f90"
  INCLUDE "loops.f90"
  INCLUDE "masks.f90"
  INCLUDE "phys.f90"
  INCLUDE "vrbls.f90"
  INCLUDE "cldwtr.f90"
  INCLUDE "pvrbls.f90"
  INCLUDE "acmclh.f90"
  REAL, DIMENSION(imjm)	:: pdsl,precrs,precss,tt,qq,ww,htmk,hbm2k
  INTEGER, DIMENSION(imjm)	:: iprec,jprec
  REAL, DIMENSION(im,jm)	:: aprec
  REAL, DIMENSION(imjm)       :: praut,pracw,err,psaut,psaci,ers,psm, &
       precrl,precsl,qi,precrl1,precsl1,qint,psm1,psm2,const,rq,ccr, &
       hh,conde,rconde,qc,tmt0,u00kl,ppr,pps,cpdr,pid,wmink
  INTEGER, DIMENSION(imjm)  	:: iwl,iwl1
  !
  INTEGER	:: i,j,l,n,npre,lml,iwk
  REAL  :: ke,mi0
  REAL  :: dtph,rdtph,twodt,rtwodt,us,eps,cclimit,climit,cws,csm1,crs1,&
       crs2,cr,aa2,ctim1,ctim2,ctim,utim,ttemp,wfix,wmin,detal,ull,aetal,&
       precrk,precsk,tk,qk,tmt15,ai,bi,qw,fi,fiw,rqkll,cwmk,expf,aa1, &
       amaxcm,tmt0k,u00klt,amaxrq,hht,qtemp,tmt0t,tmt15t,qwt,qit,qintt, &
       qct,rqt,rqtt,erk,rprs,errt,erst,amaxps,totppt
  !
  !--------------preparatory calculations---------------------------------
  dtph  = nphs * dt
  rdtph = 1. / dtph
  twodt= dtph
  rtwodt=h1/twodt
  ke=2.0e-5
  us=h1
  eps=0.622e0
  cclimit=1.0e-3
  climit=1.0e-20
  cws=0.025
  csm1=5.0000e-8
  crs1=5.00000e-6
  crs2=6.66600e-10
  cr=5.0e-4
  mi0=5.0e-4
  aa2=1.25e-3
  avrain=avrain+1.
  !-------------------padding cloud mixing ratio if too small-------------
  DO  l=1,lm
     DO  j=1,jm
        DO  i=1,im
           cwm(i,j,l)=cwm(i,j,l)*htm(i,j,l)*hbm2(i,j)
           IF(cwm(i,j,l) < 0.)cwm(i,j,l)=0.
           !------------------padding specific humidity if too small---------------
           IF(q(i,j,l) < epsq)q(i,j,l)=epsq*htm(i,j,l)
        END DO
     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
  !-----------------------------------------------------------------------
  DO j=1,jm
     DO i=1,im
        aprec(i,j)=0.
     END DO
  END DO

  DO n=1,imjm
     iwl1(n)=0
     precrl1(n)=d00
     precsl1(n)=d00
  END DO
  !------------choose the columns where prec can be produced--------------
  npre=0
  DO  j=3,jm-2
     DO  i=1,im
        DO  l=2,lm
           !zhao
           ttemp=0.025*(t(i,j,l)-273.16)
           wfix=0.9814*EXP(0.01873*l)
           wmin=0.1e-3*EXP(ttemp)*wfix
           !zhao
           IF(cwm(i,j,l) > wmin)THEN
              npre=npre+1
              iprec(npre)=i
              jprec(npre)=j
              EXIT
           END IF
        END DO
     END DO
  END DO
  !------------------------------------------------------------------------
  DO  n=1,npre
     i=iprec(n)
     j=jprec(n)
     hbm2k(n)=hbm2(i,j)
     pdsl(n)=res(i,j)*pd(i,j)
     const(n)=pdsl(n)/g*twodt
  END DO
  !-----------------begining of precipitation calculation-----------------
  !cmic$ do all autoscope
  !cmic$* private(ccr,conde,cpdr,err,ers,hh,htmk,iwl,pid,ppr,pps)
  !cmic$* private(pracw,praut,precl,precrl,precsl,psaci,psaut,psm)
  !cmic$* private(psm1,psm2,qc,qi,qint,qq,rconde,rq,tmt0,tt,u00kl,ww)
  !cmic$* private(wmink)
  DO  l=2,lm
     detal=deta(l)
     ull=ul(l)
     aetal=aeta(l)
     !zhao
     wfix=0.9814*EXP(0.01873*l)
     !zhao
     DO  n=1,npre
        precrl(n)=d00
        precsl(n)=d00
        praut(n)=d00
        psaut(n)=d00
        pracw(n)=d00
        psaci(n)=d00
        err  (n)=d00
        ers  (n)=d00
        psm  (n)=d00
        psm1 (n)=d00
        psm2 (n)=d00
        ppr  (n)=d00
        pps  (n)=d00
        cpdr (n)=d00
        hh   (n)=d00
        pid  (n)=d00
        iwl  (n)=d00
        conde(n)=d00
        rconde(n)=d00
     END DO
     !-----------------------------------------------------------------------
     DO  n=1,npre
        i=iprec(n)
        j=jprec(n)
        lml=lm-lmh(i,j)
        tt(n)=t(i,j,l)
        qq(n)=q(i,j,l)
        ww(n)=cwm(i,j,l)
        htmk(n)=htm(i,j,l)
        u00kl(n)=u00(i,j)+ul(l+lml)*(0.95-u00(i,j))*utim
        !zhao
        ttemp=0.025*(tt(n)-273.16)
        wmink(n)=0.1e-3*EXP(ttemp)*wfix
        !zhao
     END DO
     !----------choose the points where precipitation can be produced--------
     DO  n=1,npre
        precrk=amax1(0.,precrl1(n))
        precsk=amax1(0.,precsl1(n))
        hh(n)=htmk(n)*hbm2k(n)
        IF(ww(n) < wmink(n).AND.(precrk+precsk) == d00)THEN
           pid(n)=0.
        ELSE
           pid(n)=hh(n)
        END IF
     END DO
     !-------------------qw, qi and qint-------------------------------------
     DO  n=1,npre
        IF(pid(n) == 1.)THEN
           conde(n)=const(n)*detal
           rconde(n)=h1/conde(n)
           tk=tt(n)
           qk=qq(n)
           tmt0(n)=(tk-273.16)*hh(n)
           tmt15=amin1(tmt0(n),-15.)*hh(n)
           ai=0.008855
           bi=1.
           IF(tmt0(n) < -20.)THEN
              ai=0.007225
              bi=0.9674
           END IF
           qw=hh(n)*pq0/(pdsl(n)*aetal+pt) *EXP(hh(n)*a2*(tk-a3)/(tk-a4))
           qi(n)=qw*(bi+ai*amin1(tmt0(n),0.))
           qint(n)=qw*(1.-0.00032*tmt15*(tmt15+15.))
           IF(tmt0(n) <= -40.)qint(n)=qi(n)
           !-------------------ice-water id number iw------------------------------
           IF(tmt0(n) < -15.)THEN
              fi=qk-u00kl(n)*qi(n)
              IF(fi > d00.OR.ww(n) > climit) THEN
                 iwl(n)=1
              ELSE
                 iwl(n)=0
              END IF
           END IF

           IF(tmt0(n) < 0.0.AND.tmt0(n) >= -15.0)THEN
              iwl(n)=0
              IF(iwl1(n) == 1.AND.ww(n) > climit)iwl(n)=1
           END IF

           IF(tmt0(n) >= 0.)THEN
              iwl(n)=0
           END IF
           !----------------the satuation specific humidity------------------------
           fiw=float(iwl(n))
           qc(n)=(h1-fiw)*qint(n)+fiw*qi(n)
           !----------------the relative humidity----------------------------------
           IF(qc(n) <= d00)THEN
              rq(n)=1.e-10
           ELSE
              rq(n)=qk/qc(n)
           END IF
           !----------------cloud cover ratio ccr----------------------------------
           IF(rq(n) <= u00kl(n))THEN
              ccr(n)=d00
           ELSE
              rqkll=amin1(us,rq(n))
              ccr(n)=h1-SQRT((us-rqkll)/(us-u00kl(n)))
           END IF
           !-----------correct ccr if it is too small in large cwm regions--------
           IF(ccr(n) >= 0.01.AND.ccr(n) <= 0.2.AND.ww(n) >= 0.2e-3)THEN
              ccr(n)=amin1(1.0,ww(n)*1.0e3)
           END IF
        END IF
     END DO
     !------------------precipitation production rates-----------------------
     !------------------auto-convert rates-----------------------------------
     DO  n=1,npre
        IF(pid(n) == 1.)THEN
           iwk=iwl(n)
           cwmk=amax1(0.,ww(n)-climit)
           !zhao   ttn=tt(n)
           !zhao   if(ttn.ge.-10.)then
           !zhao     mi0=1.e-4
           !zhao   elseif(ttn.le.-25.)then
           !zhao     mi0=5.e-4
           !zhao   else
           !zhao     mi0=4.e-4*(-10.-ttn)/15.+1.e-4
           !zhao   endif
           mi0=wmink(n)
           !zhao
           IF(iwk == 1)THEN
              expf=EXP(0.025*tmt0(n))
              aa1=1.e-3*expf
              psaut(n)=aa1*amax1(0.,cwmk-mi0)
              cpdr(n)=-psaut(n)*twodt
              IF(-cpdr(n) >= cwmk)THEN
                 cpdr(n)=-cwmk
                 psaut(n)=-cpdr(n)*rtwodt
              END IF
           ELSE
              amaxcm=amax1(0.,cwmk-mi0)
              praut(n)=c0*amaxcm*amaxcm
              cpdr(n)=-praut(n)*twodt
              IF(-cpdr(n) >= cwmk)THEN
                 cpdr(n)=-cwmk
                 praut(n)=-cpdr(n)*rtwodt
              END IF
           END IF
           ppr(n)=praut(n)*conde(n)
           pps(n)=psaut(n)*conde(n)
        END IF
     END DO
     DO  n=1,npre
        IF(pid(n) == 1.)THEN
           ww(n)=cpdr(n)*hh(n)+ww(n)
           precrl(n)=precrl1(n)+ppr(n)*hh(n)
           precsl(n)=precsl1(n)+pps(n)*hh(n)
        END IF
     END DO
     !-----------------------accretions--------------------------------------
     DO  n=1,npre
        IF(pid(n) == 1.)THEN
           iwk=iwl(n)
           cwmk=ww(n)
           precrk=amax1(0.,precrl1(n))
           precsk=amax1(0.,precsl1(n))
           IF(iwk == 1)THEN
              expf=EXP(0.025*tmt0(n))
              cs=aa2*expf
              psaci(n)=cs*amax1(0.,cwmk)*precsk
              cpdr(n)=-psaci(n)*twodt
              IF(-cpdr(n) >= cwmk)THEN
                 cpdr(n)=-cwmk
                 psaci(n)=-cpdr(n)*rtwodt
              END IF
           ELSE
              pracw(n)=cr*amax1(0.,cwmk)*(precrk+precsk)
              cpdr(n)=-pracw(n)*twodt
              IF(-cpdr(n) >= cwmk)THEN
                 cpdr(n)=-cwmk
                 pracw(n)=-cpdr(n)*rtwodt
              END IF
           END IF
           ppr(n)=pracw(n)*conde(n)
           pps(n)=psaci(n)*conde(n)
        END IF
     END DO
     DO  n=1,npre
        IF(pid(n) == 1.)THEN
           ww(n)=cpdr(n)*hh(n)+ww(n)
           precrl(n)=precrl(n)+ppr(n)*hh(n)
           precsl(n)=precsl(n)+pps(n)*hh(n)
        END IF
     END DO
     !-----evaporation/condensation of precipitation-------------------------
     !***** err & ers positive--evaporation
     !***** err & ers negtive---condensation
     !*****
     DO  n=1,npre
        IF(pid(n) == 1.0)THEN
           qk=qq(n)
           tmt0k=tmt0(n)
           IF(tmt0k < -30.)tmt0k=-30.
           precrk=amax1(0.,precrl(n))
           precsk=amax1(0.,precsl(n))
           !
           ! increase the evaporation/condensation for strong/light prec
           !
           u00klt=u00kl(n)
           amaxrq=amax1(0.,u00kl(n)-rq(n))
           err(n)=ke*amaxrq*precrk**0.5
           IF(tmt0(n) >= 0.)THEN
              ers(n)=0.
           ELSE
              ers(n)=(crs1+crs2*tmt0k)*amaxrq*precsk/u00klt
           END IF
           !---------------correct if over-evapo./cond. occurs--------------------
           IF(err(n)+ers(n) > 1.e-20) THEN
              hht=hh(n)*twodt
              ttemp=tt(n)-rcp*(elwv*err(n)+eliv*ers(n))*hht
              qtemp=qq(n)+hht*(err(n)+ers(n))
              tmt0t=(ttemp-273.16)*hh(n)
              IF(tmt0t < -30.)tmt0t=-30.
              tmt15t=amin1(tmt0t,-15.)*hh(n)
              ai=0.008855
              bi=1.
              IF(tmt0t < -20.)THEN
                 ai=0.007225
                 bi=0.9674
              END IF
              qwt=hh(n)*pq0/(pdsl(n)*aetal+pt) *EXP(hh(n)*a2*(ttemp-a3)/(ttemp-a4))
              qit=qwt*(bi+ai*amin1(tmt0t,0.))
              qintt=qwt*(1.-0.00032*tmt15t*(tmt15t+15.))
              IF(tmt0t <= -40.)qintt=qit
              fiw=float(iwl(n))
              qct=(h1-fiw)*qintt+fiw*qit
              IF(qct <= 1.e-10) THEN
                 rqt=1.e-10
                 rqtt=1.e-10
              ELSE
                 rqt=qtemp/qct
                 rqtt=qq(n)/qct
              END IF
              IF(rqt > u00kl(n))THEN 
                 erk=(u00kl(n)-rqtt)*qct*rtwodt
                 rprs=erk/(precrk+precsk)
                 errt=precrk*rprs
                 erst=precsk*rprs
                 err(n)=amax1(d00,0.5*(err(n)+errt))
                 ers(n)=amax1(d00,0.5*(ers(n)+erst))
              END IF
           END IF
           ppr(n)=-err(n)*conde(n)
           pps(n)=-ers(n)*conde(n)
           IF(-ppr(n) >= precrk)THEN
              ppr(n)=-precrk
              err(n)=-ppr(n)*rconde(n)
           END IF
           IF(-pps(n) >= precsk)THEN
              pps(n)=-precsk
              ers(n)=-pps(n)*rconde(n)
           END IF
        END IF
     END DO
     DO  n=1,npre
        IF(pid(n) == 1.)THEN
           precrl(n)=precrl(n)+ppr(n)*hh(n)
           precsl(n)=precsl(n)+pps(n)*hh(n)
        END IF
     END DO
     !--------------------melting of the snow--------------------------------
     DO  n=1,npre
        IF(pid(n) == 1.)THEN
           cwmk=ww(n)
           amaxps=amax1(0.,precsl(n))
           IF(tmt0(n) > 0.)THEN
              psm1(n)=csm1*tmt0(n)*tmt0(n)*amaxps
              psm2(n)=cws*cr*cwmk*amaxps
              psm(n)=psm1(n)+psm2(n)
           ELSE
              psm1(n)=d00
              psm2(n)=d00
              psm(n)=d00
           END IF
           ppr(n)=psm(n)*conde(n)
           pps(n)=-psm(n)*conde(n)
           IF(-pps(n) >= amaxps)THEN
              pps(n)=-amaxps
              ppr(n)=amaxps
              psm1(n)=-pps(n)*rconde(n)
              psm2(n)=d00
              psm(n)=psm1(n)
           END IF
        END IF
     END DO
     DO  n=1,npre
        IF(pid(n) == 1.)THEN
           precrl(n)=precrl(n)+ppr(n)*hh(n)
           precsl(n)=precsl(n)+pps(n)*hh(n)
        END IF
     END DO
     !---------------update t and q------------------------------------------
     DO  n=1,npre
        IF(pid(n) == 1.)THEN
           hht=hh(n)*twodt
           tt(n)=-rcp*(elwv*err(n)+eliv*ers(n)+eliw*psm1(n)) *hht+tt(n)
           qq(n)=(err(n)+ers(n))*hht+qq(n)
        END IF
     END DO
     DO  n=1,npre
        IF(hh(n) == 1.)THEN
           iwl1(n)=iwl(n)
           precrl1(n)=precrl(n)
           precsl1(n)=precsl(n)
        END IF
     END DO
!CDIR LOOPCNT=30000
     DO  n=1,npre
        i=iprec(n)
        j=jprec(n)

        !     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.

        train(i,j,l)=train(i,j,l)+(tt(n)-t(i,j,l))*rdtph
        t(i,j,l)=tt(n)
        q(i,j,l)=qq(n)
        cwm(i,j,l)=ww(n)
     END DO
  END DO
  !-------------------the precipitation on sfc----------------------------
  DO  n=1,npre
     precrs(n)=precrl1(n)*rrow
     precss(n)=precsl1(n)*rrow
  END DO
!CDIR LOOPCNT=30000
  DO  n=1,npre
     i=iprec(n)
     j=jprec(n)
     aprec(i,j)=precrs(n)+precss(n)
     prec(i,j)=prec(i,j)+ precrs(n)+precss(n)
     acprec(i,j)=acprec(i,j)+aprec(i,j)
  END DO
  !---------------the snow and rain ratio of sfc prec---------------------
  !----sr is the ratio of snow to the total precip------------------------
  !----if total precip is zero, sr is zero--------------------------------
!CDIR LOOPCNT=30000
  DO  n=1,npre
     i=iprec(n)
     j=jprec(n)
     totppt=precrs(n)+precss(n)
     IF (totppt > 1.e-8) THEN
        sr(i,j)=precss(n)/totppt
     ELSE
        sr(i,j)=0.
     END IF
  END DO
END SUBROUTINE precpd
