    SUBROUTINE PRECPD
!     ******************************************************************
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .
! SUBPROGRAM:    PRECPD      LARGE SCALE PRECIPITATION
!   PRGRMMR: ZHAO            ORG: W/NP22     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
!   98-11-02  BLACK      - MODIFICATION FOR DISTRIBUTED MEMORY

! 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

! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE : IBM SP
!$$$
!     ******************************************************************
    PARAMETER &
    (A1=610.78,A2=17.2693882,A3=273.16,A4=35.86 &
    , PQ0=379.90516,TRESH=.95,R=287.04,C0=0.15 &
    , CP=1004.6,ELWV=2.50E6,ELIV=2.834E6,ROW=1.E3,G=9.8 &
    , EPSQ=2.E-12,DLDT=2274.0,ELIW=ELIV-ELWV)

    PARAMETER &
    (ARCP=A2*(A3-A4)/CP,RCP=1./CP,PQ0C=PQ0*TRESH,RROG=1./(ROW*G) &
    , RROW=1./ROW)
!----------------------------------------------------------------------
    INCLUDE "parmeta.f90"
    INCLUDE "parm.tbl.f90"
    INCLUDE "mpp.h"
#include "sp.h"
!----------------------------------------------------------------------
    PARAMETER &
    (IMJM=IM*JM-JM/2,JAM=6+2*(JM-10) &
    , LP1=LM+1,LTOP=1,LBOT=LM)

    PARAMETER &
    (IMJM_LOC=IDIM2*JDIM2,LDA=(IDIM2-IDIM1+1)*(JDIM2-JDIM1+1))
!-----------------------------------------------------------------------
    LOGICAL :: &
    RUN,FIRST,RESTRT,SIGMA,NOZ
!----------------------------------------------------------------------
    INCLUDE "COMM_CTLBLK.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_LOOPS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_MASKS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_PHYS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_VRBLS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_CLDWTR.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_PVRBLS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_ACMCLH.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_PPTASM.f90"
!-----------------------------------------------------------------------
    DIMENSION &
    IPREC(IMJM_LOC),JPREC(IMJM_LOC) &
    ,T_T(LM,IDIM1:IDIM2,JDIM1:JDIM2) &
    ,Q_T(LM,IDIM1:IDIM2,JDIM1:JDIM2) &
    ,TRAIN_T(LM,IDIM1:IDIM2,JDIM1:JDIM2) &
    ,HTM_T(LM,IDIM1:IDIM2,JDIM1:JDIM2) &
    ,CWM_T(LM,IDIM1:IDIM2,JDIM1:JDIM2) &
    ,TLAT_T(LM,IDIM1:IDIM2,JDIM1:JDIM2)
!-----------------------------------------------------------------------
    DIMENSION &
    PRECRL1(IMJM_LOC),PRECSL1(IMJM_LOC),IWL1(IMJM_LOC)
!-----------------------------------------------------------------------
    REAL :: &
    KE,INIT,MI0
!-----------------------------------------------------------------------
!***********************************************************************

    print*,"entrei!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
!--------------PREPARATORY CALCULATIONS---------------------------------
    DTPH  = NPHS * DT
    RDTPH = 1. / DTPH
    TWODT= DTPH
    RTWODT=1./TWODT
    KE=2.0E-5
    US=1.
    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.
    ARATIM=ARATIM+1
!-------------------PADDING CLOUD MIXING RATIO IF TOO SMALL-------------
! omp parallel do
    DO 20 L=1,LM
        DO J=MYJS,MYJE
            DO I=MYIS,MYIE
                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)
            ENDDO
        ENDDO
    20 END DO

    UTIM=1.
!-----------------------------------------------------------------------
! omp parallel do
    DO N=1,IMJM_LOC
        IWL1(N)=0
        PRECRL1(N)=0.
        PRECSL1(N)=0.
    ENDDO
!------------CHOOSE THE COLUMNS WHERE PREC CAN BE PRODUCED--------------
    NPRE=0
    DO 35 J=MYJS2,MYJE2
        DO 35 I=MYIS,MYIE
        
            DO L=2,LM
                TTEMP=0.025*(T(I,J,L)-273.16)
                WFIX=0.9814*EXP(0.01873*L)
                WMIN=0.1E-3*EXP(TTEMP)*WFIX
                IF(CWM(I,J,L) > WMIN)GO TO 33

            ENDDO
        
            GO TO 35
            33 NPRE=NPRE+1
            IPREC(NPRE)=I
            JPREC(NPRE)=J
    35 END DO
!------------------------------------------------------------------------
!***
!***  TRANSPOSE ARRAYS
!***
! omp parallel sections
! omp section
    CALL SGETMO(T,LDA,LDA,LM,T_T,LM)
! omp section
    CALL SGETMO(Q,LDA,LDA,LM,Q_T,LM)
! omp section
    CALL SGETMO(CWM,LDA,LDA,LM,CWM_T,LM)
! omp section
    CALL SGETMO(HTM,LDA,LDA,LM,HTM_T,LM)
! omp section
    CALL SGETMO(TRAIN,LDA,LDA,LM,TRAIN_T,LM)
! omp section
    CALL SGETMO(TLATGS,LDA,LDA,LM,TLAT_T,LM)
! omp end parallel sections
!-----------------------------------------------------------------------
!-----------------BEGINING OF PRECIPITATION CALCULATION-----------------
!-----------------------------------------------------------------------
!***
!***  LOOP OVER ALL POSSIBLE PRECIPITATION POINTS
!***
! omp parallel do
! omp& private(aai,aetal,ai,amaxcm,amaxps,amaxrq,bi,ccr,conde,
! omp&         const,cpdr,cs,cwmk,detal,erk,err,errt,ers,erst,expf,fi,
! omp&         fiw,hbm2k,hh,htmk,i,iwl,j,lml,mi0,pdsl,pid,ppr,pps,
! omp&         pracw,praut,precrk,precrl,precsk,precsl,precss,psaci,
! omp&         psaut,psm,psm1,psm2,qc,qct,qi,qint,qintt,qit,qk,qq,
! omp&         qtemp,qw,qwt,rconde,rprs,rq,rqkll,rqt,rqtt,tk,tmt0,
! omp&         tmt0k,tmt0t,tmt15,tmt15t,totppt,tt,ttemp,ttlat,u00ij,
! omp&         u00kl,u00klt,ull,wfix,wmink,ww)

    DO 300 N=1,NPRE
    
        I=IPREC(N)
        J=JPREC(N)
        HBM2K=HBM2(I,J)
        PDSL=RES(I,J)*PD(I,J)
        CONST=PDSL/G*TWODT
        LML=LM-LMH(I,J)
        U00IJ=U00(I,J)
    
        DO 180 L=2,LM
        
            DETAL=DETA(L)
            ULL=UL(L)
            AETAL=AETA(L)
            WFIX=0.9814*EXP(0.01873*L)
        
            PRECRL=0.
            PRECSL=0.
            PRAUT=0.
            PSAUT=0.
            PRACW=0.
            PSACI=0.
            ERR  =0.
            ERS  =0.
            PSM  =0.
            PSM1 =0.
            PSM2 =0.
            PPR  =0.
            PPS  =0.
            CPDR =0.
            HH   =0.
            PID  =0.
            IWL  =0.
            CONDE=0.
            RCONDE=0.
        !-----------------------------------------------------------------------
            TT=T_T(L,I,J)
            QQ=Q_T(L,I,J)
            WW=CWM_T(L,I,J)
            HTMK=HTM_T(L,I,J)
            TTLAT=TLAT_T(L,I,J)
        !-----------------------------------------------------------------------
            U00KL=U00IJ+UL(L+LML)*(0.95-U00IJ)*UTIM
            TTEMP=0.025*(TT-273.16)
            WMINK=0.1E-3*EXP(TTEMP)*WFIX
        !-----------------------------------------------------------------------
        !----------CHOOSE THE POINTS WHERE PRECIPITATION CAN BE PRODUCED--------
        !-----------------------------------------------------------------------
            PRECRK=AMAX1(0.,PRECRL1(N))
            PRECSK=AMAX1(0.,PRECSL1(N))
            HH=HTMK*HBM2K
            IF(WW < WMINK .AND. (PRECRK+PRECSK) == 0.)THEN
                PID=0.
            ELSE
                PID=HH
            ENDIF
        !-----------------------------------------------------------------------
        !-------------------QW, QI AND QINT-------------------------------------
        !-----------------------------------------------------------------------
            IF(PID == 1.)THEN
                CONDE=CONST*DETAL
                RCONDE=1./CONDE
                TK=TT
                QK=QQ
                TMT0=(TK-273.16)*HH
                TMT15=AMIN1(TMT0,-15.)*HH
                AI=0.008855
                BI=1.
            
                IF(TMT0 < -20.)THEN
                    AI=0.007225
                    BI=0.9674
                ENDIF
            
                QW=HH*PQ0/(PDSL*AETAL+PT) &
                *EXP(HH*A2*(TK-A3)/(TK-A4))
                QI=QW*(BI+AI*AMIN1(TMT0,0.))
                QINT=QW*(1.-0.00032*TMT15*(TMT15+15.))
                IF(TMT0 <= -40.)QINT=QI
            !-------------------ICE-WATER ID NUMBER IW------------------------------
                IF(TMT0 < -15.)THEN
                    FI=QK-U00KL*QI
                    IF(FI > 0. .OR. WW > CLIMIT) THEN
                        IWL=1
                    ELSE
                        IWL=0
                    ENDIF
                ENDIF
            
                IF(TMT0 < 0.0 .AND. TMT0 >= -15.0)THEN
                    IWL=0
                    IF(IWL1(N) == 1 .AND. WW > CLIMIT)IWL=1
                ENDIF
            
                IF(TMT0 >= 0.)THEN
                    IWL=0
                ENDIF
            !----------------THE SATUATION SPECIFIC HUMIDITY------------------------
                FIW=FLOAT(IWL)
                QC=(1.-FIW)*QINT+FIW*QI
            !----------------THE RELATIVE HUMIDITY----------------------------------
                IF(QC <= 0.)THEN
                    RQ=1.E-10
                ELSE
                    RQ=QK/QC
                ENDIF
            !----------------CLOUD COVER RATIO CCR----------------------------------
                IF(RQ <= U00KL)THEN
                    CCR=0.
                ELSE
                    RQKLL=AMIN1(US,RQ)
                    CCR=1.-SQRT((US-RQKLL)/(US-U00KL))
                ENDIF
            !-----------CORRECT CCR IF IT IS TOO SMALL IN LARGE CWM REGIONS--------
                IF(CCR >= 0.01 .AND. CCR <= 0.2.AND.WW >= 0.2E-3)THEN
                    CCR=AMIN1(1.,WW*1.0E3)
                ENDIF
            ENDIF
            60 CONTINUE
        !-----------------------------------------------------------------------
        !------------------PRECIPITATION PRODUCTION RATES-----------------------
        !------------------AUTO-CONVERT RATES-----------------------------------
        !-----------------------------------------------------------------------
            IF(PID == 1.)THEN
                IWK=IWL
                CWMK=AMAX1(0.,WW-CLIMIT)
                MI0=WMINK
            
                IF(IWK == 1)THEN
                    EXPF=EXP(0.025*TMT0)
                    AA1=1.E-3*EXPF
                    PSAUT=AA1*AMAX1(0.,CWMK-MI0)
                    CPDR=-PSAUT*TWODT
                    IF(-CPDR >= CWMK)THEN
                        CPDR=-CWMK
                        PSAUT=-CPDR*RTWODT
                    ENDIF
                ELSE
                    AMAXCM=AMAX1(0.,CWMK-MI0)
                    PRAUT=C0*AMAXCM*AMAXCM
                    CPDR=-PRAUT*TWODT
                    IF(-CPDR >= CWMK)THEN
                        CPDR=-CWMK
                        PRAUT=-CPDR*RTWODT
                    ENDIF
                ENDIF
                PPR=PRAUT*CONDE
                PPS=PSAUT*CONDE
            ENDIF
        
            IF(PID == 1.)THEN
                WW=CPDR*HH+WW
                PRECRL=PRECRL1(N)+PPR*HH
                PRECSL=PRECSL1(N)+PPS*HH
            ENDIF
        !-----------------------------------------------------------------------
        !-----------------------ACCRETIONS--------------------------------------
        !-----------------------------------------------------------------------
            IF(PID == 1.)THEN
                IWK=IWL
                CWMK=WW
                PRECRK=AMAX1(0.,PRECRL1(N))
                PRECSK=AMAX1(0.,PRECSL1(N))
                IF(IWK == 1)THEN
                    EXPF=EXP(0.025*TMT0)
                    CS=AA2*EXPF
                    PSACI=CS*AMAX1(0.,CWMK)*PRECSK
                    CPDR=-PSACI*TWODT
                    IF(-CPDR >= CWMK)THEN
                        CPDR=-CWMK
                        PSACI=-CPDR*RTWODT
                    ENDIF
                ELSE
                    PRACW=CR*AMAX1(0.,CWMK)*(PRECRK+PRECSK)
                    CPDR=-PRACW*TWODT
                    IF(-CPDR >= CWMK)THEN
                        CPDR=-CWMK
                        PRACW=-CPDR*RTWODT
                    ENDIF
                ENDIF
                PPR=PRACW*CONDE
                PPS=PSACI*CONDE
            ENDIF
        
            IF(PID == 1.)THEN
                WW=CPDR*HH+WW
                PRECRL=PRECRL+PPR*HH
                PRECSL=PRECSL+PPS*HH
            ENDIF
        !-----------------------------------------------------------------------
        !-----EVAPORATION/CONDENSATION OF PRECIPITATION-------------------------
        !***** ERR & ERS POSITIVE--EVAPORATION
        !***** ERR & ERS NEGTIVE---CONDENSATION
        !-----------------------------------------------------------------------
            IF(PID == 1.0)THEN
                QK=QQ
                TMT0K=TMT0
                IF(TMT0K < -30.)TMT0K=-30.
                PRECRK=AMAX1(0.,PRECRL)
                PRECSK=AMAX1(0.,PRECSL)
            !----------------------------------------------------------------------
            ! INCREASE THE EVAPORATION/CONDENSATION FOR STRONG/LIGHT PREC
            !----------------------------------------------------------------------
                U00KLT=U00KL
                AMAXRQ=AMAX1(0.,U00KL-RQ)
                ERR=KE*AMAXRQ*PRECRK**0.5
            
                IF(TMT0 >= 0.)THEN
                    ERS=0.
                ELSE
                    ERS=(CRS1+CRS2*TMT0K)*AMAXRQ*PRECSK/U00KLT
                ENDIF
            
                IF(ERR+ERS <= 1.E-20) GO TO 125
            !---------------CORRECT IF OVER-EVAPO./COND. OCCURS--------------------
                HHT=HH*TWODT
                TTEMP=TT-RCP*(ELWV*ERR+ELIV*ERS)*HHT
                QTEMP=QQ+HHT*(ERR+ERS)
                TMT0T=(TTEMP-273.16)*HH
                IF(TMT0T < -30.)TMT0T=-30.
                TMT15T=AMIN1(TMT0T,-15.)*HH
                AI=0.008855
                BI=1.
            
                IF(TMT0T < -20.)THEN
                    AI=0.007225
                    BI=0.9674
                ENDIF
            
                QWT=HH*PQ0/(PDSL*AETAL+PT) &
                *EXP(HH*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)
                QCT=(1.-FIW)*QINTT+FIW*QIT
            
                IF(QCT <= 1.E-10) THEN
                    RQT=1.E-10
                    RQTT=1.E-10
                ELSE
                    RQT=QTEMP/QCT
                    RQTT=QQ/QCT
                ENDIF
            
                IF(RQT <= U00KL) GO TO 125
            
                ERK=(U00KL-RQTT)*QCT*RTWODT
                RPRS=ERK/(PRECRK+PRECSK)
                ERRT=PRECRK*RPRS
                ERST=PRECSK*RPRS
                ERR=AMAX1(0.,0.5*(ERR+ERRT))
                ERS=AMAX1(0.,0.5*(ERS+ERST))
            
                125 CONTINUE
            
                PPR=-ERR*CONDE
                PPS=-ERS*CONDE
            
                IF(-PPR >= PRECRK)THEN
                    PPR=-PRECRK
                    ERR=-PPR*RCONDE
                ENDIF
            
                IF(-PPS >= PRECSK)THEN
                    PPS=-PRECSK
                    ERS=-PPS*RCONDE
                ENDIF
            
            ENDIF
        
            IF(PID == 1.)THEN
                PRECRL=PRECRL+PPR*HH
                PRECSL=PRECSL+PPS*HH
            ENDIF
        !-----------------------------------------------------------------------
        !--------------------MELTING OF THE SNOW--------------------------------
        !-----------------------------------------------------------------------
            IF(PID == 1.)THEN
                CWMK=WW
                AMAXPS=AMAX1(0.,PRECSL)
            
                IF(TMT0 > 0.)THEN
                    PSM1=CSM1*TMT0*TMT0*AMAXPS
                    PSM2=CWS*CR*CWMK*AMAXPS
                    PSM=PSM1+PSM2
                ELSE
                    PSM1=0.
                    PSM2=0.
                    PSM=0.
                ENDIF
            
                PPR=PSM*CONDE
                PPS=-PSM*CONDE
            
                IF(-PPS >= AMAXPS)THEN
                    PPS=-AMAXPS
                    PPR=AMAXPS
                    PSM1=-PPS*RCONDE
                    PSM2=0.
                    PSM=PSM1
                ENDIF
            
            ENDIF
        
            IF(PID == 1.)THEN
                PRECRL=PRECRL+PPR*HH
                PRECSL=PRECSL+PPS*HH
            ENDIF
        !-----------------------------------------------------------------------
        !---------------UPDATE T AND Q------------------------------------------
        !-----------------------------------------------------------------------
            IF(PID == 1.)THEN
                HHT=HH*TWODT
                TT=-RCP*(ELWV*ERR+ELIV*ERS+ELIW*PSM1) &
                *HHT+TT
                QQ=(ERR+ERS)*HHT+QQ
                TTLAT=-RCP*(ELWV*ERR+ELIV*ERS+ELIW*PSM1 &
                -ELWV*(PRAUT+PRACW) &
                -ELIV*(PSAUT+PSACI))*HHT+TTLAT
            ENDIF
        
            IF(HH == 1.)THEN
                IWL1(N)=IWL
                PRECRL1(N)=PRECRL
                PRECSL1(N)=PRECSL
            ENDIF
        
        !     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_T(L,I,J)=TRAIN_T(L,I,J)+(TT-T_T(L,I,J))*RDTPH
            T_T(L,I,J)=TT
            Q_T(L,I,J)=QQ
            CWM_T(L,I,J)=WW
            TLAT_T(L,I,J)=TTLAT
        180 END DO
    !-----------------------------------------------------------------------
    !-------------------THE PRECIPITATION ON SFC----------------------------
    !-----------------------------------------------------------------------
        PRECRS=PRECRL1(N)*RROW
        PRECSS=PRECSL1(N)*RROW
    
        APREC(I,J)=PRECRS+PRECSS
        PREC(I,J)=PREC(I,J)+PRECRS+PRECSS
        ACPREC(I,J)=ACPREC(I,J)+APREC(I,J)
    ! LG      if (ACPREC(I,J).gt.0) then
    ! LG        print*,I,J,ACPREC(I,J)
    ! LG      endif
    !-----------------------------------------------------------------------
    !---------------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--------------------------------
    !-----------------------------------------------------------------------
        TOTPPT=PRECRS+PRECSS
        IF (TOTPPT > 1.E-8) THEN
            SR(I,J)=PRECSS/TOTPPT
        ELSE
            SR(I,J)=0.
        ENDIF
    !-----------------------------------------------------------------------
    300 END DO
!-----------------------------------------------------------------------
!***
!***  TRANSPOSE BACK
!***
! omp parallel sections
! omp section
    CALL SGETMO(T_T,LM,LM,LDA,T,LDA)
! omp section
    CALL SGETMO(Q_T,LM,LM,LDA,Q,LDA)
! omp section
    CALL SGETMO(CWM_T,LM,LM,LDA,CWM,LDA)
! omp section
    CALL SGETMO(TRAIN_T,LM,LM,LDA,TRAIN,LDA)
! omp section
    CALL SGETMO(TLAT_T,LM,LM,LDA,TLATGS,LDA)
! omp end parallel sections
!-----------------------------------------------------------------------
    RETURN
    END SUBROUTINE PRECPD
