CALWXT_REVISED.f

       SUBROUTINE calwxt_revised_post(T,Q,PMID,PINT,HTM,LMH,PREC,ZINT,IWX)
! 
!     FILE: CALWXT.f
!     WRITTEN: 11 NOVEMBER 1993, MICHAEL BALDWIN
!     REVISIONS:
!               30 SEPT 1994-SETUP NEW DECISION TREE (M BALDWIN)
!               12 JUNE 1998-CONVERSION TO 2-D (T BLACK)
!     01-10-25  H CHUANG - MODIFIED TO PROCESS HYBRID MODEL OUTPUT
!     02-01-15  MIKE BALDWIN - WRF VERSION
!     05-07-07  BINBIN ZHOU  - ADD PREC FOR RSM
!     05-08-24  GEOFF MANIKIN - MODIFIED THE AREA REQUIREMENTS
!                TO MAKE AN ALTERNATE ALGORITHM 
!     19-10-30  Bo CUI - REMOVE "GOTO" STATEMENT
!                              
!
!     ROUTINE TO COMPUTE PRECIPITATION TYPE USING A DECISION TREE
!     APPROACH THAT USES VARIABLES SUCH AS INTEGRATED WET BULB TEMP
!     BELOW FREEZING AND LOWEST LAYER TEMPERATURE
!
!     SEE BALDWIN AND CONTORNO PREPRINT FROM 13TH WEATHER ANALYSIS
!     AND FORECASTING CONFERENCE FOR MORE DETAILS
!     (OR BALDWIN ET AL, 10TH NWP CONFERENCE PREPRINT)
!
!     SINCE THE ORIGINAL VERSION OF THE ALGORITHM HAS A HIGH BIAS
!      FOR FREEZING RAIN AND SLEET, THE GOAL IS TO BALANCE THAT BIAS
!      WITH A VERSION MORE LIKELY TO PREDICT SNOW
!
     use params_mod, only: h1m12, d00, d608, h1, rog
     use ctlblk_mod, only: jsta, jend, modelname, pthresh, im, jsta_2l, jend_2u, lm,&
              lp1, spval
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      implicit none
!
!  LIST OF VARIABLES NEEDED
!    PARAMETERS:
!      D608,ROG,H1,D00
!HC       PARAMETER(D608=0.608,ROG=287.04/9.8,H1=1.0,D00=0.0)
!
!    INPUT:
!      T,Q,PMID,HTM,LMH,PREC,ZINT
      REAL,dimension(IM,jsta_2l:jend_2u,LM), intent(in) ::  t,q,pmid,htm
      REAL,dimension(IM,jsta_2l:jend_2u,LP1),intent(in) ::  pint,zint 
      REAL,dimension(IM,jsta_2l:jend_2u),    intent(in) ::  lmh
      REAL,dimension(IM,jsta_2l:jend_2u),    intent(in) ::  prec
!    OUTPUT:
!      IWX - INSTANTANEOUS WEATHER TYPE.
!        ACTS LIKE A 4 BIT BINARY
!          1111 = RAIN/FREEZING RAIN/ICE PELLETS/SNOW
!          WHERE THE ONE'S DIGIT IS FOR SNOW
!                THE TWO'S DIGIT IS FOR ICE PELLETS
!                THE FOUR'S DIGIT IS FOR FREEZING RAIN
!            AND THE EIGHT'S DIGIT IS FOR RAIN
     integer, DIMENSION(IM,jsta:jend),intent(inout) ::  iwx
!    INTERNAL:
!
      REAL, ALLOCATABLE :: twet(:,:,:)
      integer,DIMENSION(IM,jsta:jend) :: karr,licee
      real,   dimension(IM,jsta:jend) :: tcold,twarm
!
      integer i,j,l,lmhk,lice,ifrel,iwrml,ifrzl
      real psfck,tdchk,a,tdkl,tdpre,tlmhk,twrmk,areas8,areap4,area1,  &
           surfw,surfc,dzkl,pintk1,pintk2,pm150,qkl,tkl,pkl,area0,    &
           areap0
     
      logical :: jcontinue=.true.

!    SUBROUTINES CALLED:
!     WETBULB
!     
!
!     INITIALIZE WEATHER TYPE ARRAY TO ZERO (IE, OFF).
!     WE DO THIS SINCE WE WANT IWX TO REPRESENT THE
!     INSTANTANEOUS WEATHER TYPE ON RETURN.
!     
!
!     ALLOCATE LOCAL STORAGE
!
      ALLOCATE ( twet(im,jsta_2l:jend_2u,lm) )
!
!$omp  parallel do
      DO j=jsta,jend
        DO i=1,im
          iwx(i,j) = 0
        ENDDO
      ENDDO

!
!!$omp  parallel do
!!$omp& private(a,lmhk,pkl,psfck,qkl,tdchk,tdkl,tdpre,tkl)
      DO 800 j=jsta,jend
      DO 800 i=1,im
      lmhk=nint(lmh(i,j))
!
!   SKIP THIS POINT IF NO PRECIP THIS TIME STEP 
!
      IF (prec(i,j)<=pthresh) cycle    
!
!   FIND COLDEST AND WARMEST TEMPS IN SATURATED LAYER BETWEEN
!   70 MB ABOVE GROUND AND 500 MB
!   ALSO FIND HIGHEST SATURATED LAYER IN THAT RANGE
!
!meb
      psfck=pint(i,j,lmhk+1)
!meb
      tdchk=2.0

      jcontinue=.true.
      do while (jcontinue)

  760 tcold(i,j)=t(i,j,lmhk)
      twarm(i,j)=t(i,j,lmhk)
      licee(i,j)=lmhk
!
      DO 775 l=1,lmhk
      qkl=q(i,j,l)
      qkl=amax1(h1m12,qkl)
      tkl=t(i,j,l)
      pkl=pmid(i,j,l)
!
!   SKIP PAST THIS IF THE LAYER IS NOT BETWEEN 70 MB ABOVE GROUND
!       AND 500 MB
!
      IF (pkl<50000.0.OR.pkl>psfck-7000.0) cycle   
      IF(qkl<spval)THEN
      a=alog(qkl*pkl/(610.78*(0.378*qkl+0.622)))
      tdkl=(237.3*a)/(17.269-a)+273.15
      tdpre=tkl-tdkl
      IF (tdpre<tdchk.AND.tkl<tcold(i,j)) tcold(i,j)=tkl
      IF (tdpre<tdchk.AND.tkl>twarm(i,j)) twarm(i,j)=tkl
      IF (tdpre<tdchk.AND.l<licee(i,j)) licee(i,j)=l
      ENDIF
  775 CONTINUE
!
!    IF NO SAT LAYER AT DEW POINT DEP=TDCHK, INCREASE TDCHK
!     AND START AGAIN (BUT DON'T MAKE TDCHK > 6)
!
      IF (tcold(i,j)==t(i,j,lmhk).AND.tdchk<6.0) THEN
        tdchk=tdchk+2.0
      ELSE
        jcontinue=.false.
      ENDIF
      enddo     ! enddo jcontinue
  800 CONTINUE
!
!    LOWEST LAYER T
!
      DO 850 j=jsta,jend
      DO 850 i=1,im
      karr(i,j)=0
      IF (prec(i,j)<=pthresh) cycle    
      lmhk=nint(lmh(i,j))
      tlmhk=t(i,j,lmhk)
!
!    DECISION TREE TIME
!
      IF (tcold(i,j)>269.15) THEN
          IF (tlmhk<=273.15) THEN
!             TURN ON THE FLAG FOR
!             FREEZING RAIN = 4
!             IF ITS NOT ON ALREADY
!             IZR=MOD(IWX(I,J),8)/4
!             IF (IZR<1) IWX(I,J)=IWX(I,J)+4
              iwx(i,j)=iwx(i,j)+4
            cycle   
          ELSE
!             TURN ON THE FLAG FOR
!             RAIN = 8
!             IF ITS NOT ON ALREADY
!             IRAIN=IWX(I,J)/8
!             IF (IRAIN<1) IWX(I,J)=IWX(I,J)+8
              iwx(i,j)=iwx(i,j)+8
            cycle    
          ENDIF
      ENDIF
      karr(i,j)=1
  850 CONTINUE
!
!   COMPUTE WET BULB ONLY AT POINTS THAT NEED IT
!
      CALL wetbulb(t,q,pmid,htm,karr,twet)
!
!!$omp  parallel do
!!$omp& private(area1,areap4,areap0,areas8,dzkl,ifrzl,iwrml,lice,
!!$omp&         lmhk,pintk1,pintk2,pm150,psfck,surfc,surfw,
!!$omp&         tlmhk,twrmk)
      DO 1900 j=jsta,jend
      DO 1900 i=1,im
      IF(karr(i,j)>0)THEN
        lmhk=nint(lmh(i,j))
        lice=licee(i,j)
!meb
        psfck=pint(i,j,lmhk+1)
!meb
        tlmhk=t(i,j,lmhk)
        twrmk=twarm(i,j)
!
!    TWET AREA VARIABLES
!     CALCULATE ONLY WHAT IS NEEDED
!      FROM GROUND TO 150 MB ABOVE SURFACE
!      FROM GROUND TO TCOLD LAYER
!      AND FROM GROUND TO 1ST LAYER WHERE WET BULB T < 0.0
!
!     PINTK1 IS THE PRESSURE AT THE BOTTOM OF THE LAYER
!     PINTK2 IS THE PRESSURE AT THE TOP OF THE LAYER
!
!     AREAP4 IS THE AREA OF TWET ABOVE -4 C BELOW HIGHEST SAT LYR 
!     AREAP0 IS THE AREA OF TWET ABOVE 0 C BELOW HIGHEST SAT LYR
!
        areas8=d00
        areap4=d00
    areap0=d00
        surfw =d00
        surfc =d00
        
!
        DO 1945 l=lmhk,lice,-1
        dzkl=zint(i,j,l)-zint(i,j,l+1)
        area1=(twet(i,j,l)-269.15)*dzkl
        area0=(twet(i,j,l)-273.15)*dzkl
        IF (twet(i,j,l)>=269.15) areap4=areap4+area1
        IF (twet(i,j,l)>=273.15) areap0=areap0+area0
 1945   CONTINUE
!
!        IF (AREAP4<3000.0) THEN
!             TURN ON THE FLAG FOR
!             SNOW = 1
!             IF ITS NOT ON ALREADY
!             ISNO=MOD(IWX(I,J),2)
!             IF (ISNO<1) IWX(I,J)=IWX(I,J)+1
!          IWX(I,J)=IWX(I,J)+1
!          GO TO 1900
!        ENDIF
        IF (areap0<350.0) THEN
!             TURN ON THE FLAG FOR
!             SNOW = 1
              iwx(i,j)=iwx(i,j)+1
            cycle     
       ENDIF
!
!     AREAS8 IS THE NET AREA OF TWET W.R.T. FREEZING IN LOWEST 150MB
!
        pintk1=psfck
        pm150=psfck-15000.
!
        DO 1955 l=lmhk,1,-1
        pintk2=pint(i,j,l)
        IF(pintk1<pm150) THEN
          pintk1=pintk2
        ELSE
          dzkl=zint(i,j,l)-zint(i,j,l+1)
!
!    SUM PARTIAL LAYER IF IN 150 MB AGL LAYER
!
          IF(pintk2<pm150)                                   &
            dzkl=t(i,j,l)*(q(i,j,l)*d608+h1)*rog*               &
                 alog(pintk1/pm150)
          area1=(twet(i,j,l)-273.15)*dzkl
          areas8=areas8+area1
          pintk1=pintk2
        ENDIF
 1955   CONTINUE
!
!     SURFW IS THE AREA OF TWET ABOVE FREEZING BETWEEN THE GROUND
!       AND THE FIRST LAYER ABOVE GROUND BELOW FREEZING
!     SURFC IS THE AREA OF TWET BELOW FREEZING BETWEEN THE GROUND
!       AND THE WARMEST SAT LAYER
!
        ifrzl=0
        iwrml=0
!
        DO 2050 l=lmhk,1,-1
        IF (ifrzl==0.AND.t(i,j,l)<273.15) ifrzl=1
        IF (iwrml==0.AND.t(i,j,l)>=twrmk) iwrml=1
!
        IF (iwrml==0.OR.ifrzl==0) THEN
          dzkl=zint(i,j,l)-zint(i,j,l+1)
          area1=(twet(i,j,l)-273.15)*dzkl
          IF(ifrzl==0.AND.twet(i,j,l)>=273.15)surfw=surfw+area1
          IF(iwrml==0.AND.twet(i,j,l)<=273.15)surfc=surfc+area1
        ENDIF
 2050   CONTINUE
        IF(surfc<-3000.0.OR.                                    &
     &    (areas8<-3000.0.AND.surfw<50.0)) THEN
!             TURN ON THE FLAG FOR
!             ICE PELLETS = 2
!             IF ITS NOT ON ALREADY
!             IIP=MOD(IWX(I,J),4)/2
!             IF (IIP<1) IWX(I,J)=IWX(I,J)+2
          iwx(i,j)=iwx(i,j)+2
          cycle     
        ENDIF
!
        IF(tlmhk<273.15) THEN
!             TURN ON THE FLAG FOR
!             FREEZING RAIN = 4
!             IF ITS NOT ON ALREADY
!             IZR=MOD(IWX(K),8)/4
!             IF (IZR<1) IWX(K)=IWX(K)+4
          iwx(i,j)=iwx(i,j)+4
        ELSE
!             TURN ON THE FLAG FOR
!             RAIN = 8
!             IF ITS NOT ON ALREADY
!             IRAIN=IWX(K)/8
!             IF (IRAIN<1) IWX(K)=IWX(K)+8
          iwx(i,j)=iwx(i,j)+8
        ENDIF
      ENDIF
 1900 CONTINUE
!      print *, 'revised check ', IWX(500,800)

!---------------------------------------------------------
      DEALLOCATE (twet)


      RETURN
      END