CALWXT.f

       SUBROUTINE calwxt_post(T,Q,PMID,PINT,HTM,LMH,PREC,ZINT,IWX,ZWET)
! 
!     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
!     19-10-30  Bo CUI - REMOVE "GOTO" STATEMENT
!     21-07-26  Wen Meng - Restrict computation from undefined grids
!                              
!
!     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)
! 
      use params_mod, only: h1m12, d00, d608, h1, rog
      use ctlblk_mod, only: jsta, jend, spval, modelname,pthresh, im,   &
                            jsta_2l, jend_2u, lm, lp1
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      implicit none
!
!    INPUT:
!      T,Q,PMID,HTM,LMH,PREC,ZINT
!
      real,dimension(IM,jsta_2l:jend_2u),intent(in)     :: lmh
      real,dimension(IM,jsta_2l:jend_2u,LM),intent(in)  :: t,q,pmid,htm
      real,dimension(IM,jsta_2l:jend_2u,LP1),intent(in) :: zint,pint
      integer,DIMENSION(IM,jsta:jend),intent(inout)     :: iwx
      real,dimension(IM,jsta_2l:jend_2u),intent(inout) :: prec
      real,DIMENSION(IM,jsta:jend),intent(inout)       :: zwet


!    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
!
!    INTERNAL:
!
      REAL, ALLOCATABLE :: twet(:,:,:)
      integer,DIMENSION(IM,jsta:jend) :: karr,licee
      real,   DIMENSION(IM,jsta:jend) :: tcold,twarm

      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
!

      integer i,j,l,lmhk,lice,ifrel,iwrml,ifrzl
      real psfck,tdchk,a,tdkl,tdpre,tlmhk,twrmk,areas8,areap4,   &
           surfw,surfc,dzkl,area1,pintk1,pintk2,pm150,pkl,tkl,qkl

      ALLOCATE ( twet(im,jsta_2l:jend_2u,lm) )

!
!!$omp  parallel do
      DO j=jsta,jend
        DO i=1,im
          iwx(i,j)  = 0
          zwet(i,j) = spval
!           if (I == 324 .and. J == 390) then
!           LMHK = NINT(LMH(I,J))
!           DO L=LMHK,1,-1
!             print *, 'tprof ', L, T(I,J,L)
!           ENDDO
!         endif
        ENDDO
      ENDDO

      IF(modelname=='RSM') THEN          !add by Binbin because of different unit
       DO j=jsta,jend
       DO i=1,im
        prec(i,j) = prec(i,j)*3*3600.0
       ENDDO
       ENDDO
      END IF


!
!!$omp  parallel do 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 = max(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)
      CALL wetfrzlvl(twet,zwet)
!
!!$omp  parallel do                                                 &
!    & private(area1,areap4,areas8,dzkl,ifrzl,iwrml,lice,          &
!    &         lmhk,pintk1,pintk2,pm150,psfck,surfc,surfw,         &
!    &         tlmhk,twrmk)
      DO 1900 j=jsta,jend
      DO 1900 i=1,im
!       IF (I == 324 .AND. J == 390) THEN
!          LMHK=NINT(LMH(I,J))
!          DO L=LMHK,1,-1          
!           print *, 'TW NCEP ', TWET(I,J,L)
!          ENDDO
!       ENDIF
      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 
!
        areas8=d00
        areap4=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
        IF (twet(i,j,l)>=269.15) areap4=areap4+area1
 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
          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
!---------------------------------------------------------
      DEALLOCATE (twet)

      IF(modelname == 'RSM') THEN    !add by Binbin, change back
!!$omp parallel do private(i,j)
        DO j=jsta,jend
          DO i=1,im
            prec(i,j) = prec(i,j)/(3*3600.0)
          ENDDO
        ENDDO
      END IF

      RETURN
      END