WETBULB.f

       SUBROUTINE wetbulb(T,Q,PMID,HTM,KARR,TWET)
!
!     FILE: WETBULB.f
!     WRITTEN: 10 SEPTEMBER 1993, MICHAEL BALDWIN
!     REVISIONS: 
!     CONVERSION TO 2-D: 12 JUNE 1998 (T BLACK)
!     MPI VERSION: 04 Jan 2000 ( JIM TUCCILLO )
!     MODIFIED FOR HYBRID: OCT 2001, H CHUANG
!     02-01-15  MIKE BALDWIN - WRF VERSION
!     21-07-26  Wen Meng  - Restrict compuation from undefined grids
!     21-09-13  Jesse Meng- 2D DECOMPOSITION
!
!-----------------------------------------------------------------------
!     ROUTINE TO COMPUTE WET BULB TEMPERATURES USING THE LOOK UP TABLE
!     APPROACH THAT IS USED IN CUCNVC
!  
!     FOR A GIVEN POINT K AND LAYER L:
!      THETA E IS COMPUTED FROM THETA AND Q BY LIFTING THE PARCEL TO
!      ITS SATURATION POINT.
!      THEN THE WET BULB TEMPERATURE IS FOUND BY FOLLOWING THE THETA E
!      DOWN TO THE ORIGINAL PRESSURE LEVEL (USING SUBROUTINE TTBLEX).     
!
!   
      use lookup_mod, only: thl, rdth, jtb, qs0, sqs, rdq, itb, ptbl, plq, ttbl,&
              pl, rdp, the0, sthe, rdthe, ttblq, itbq, jtbq, rdpq, the0q, stheq,&
              rdtheq
      use ctlblk_mod, only: jsta, jend, im, jsta_2l, jend_2u, lm, spval,        &
                            ista, iend, ista_2l, iend_2u
      use cuparm_mod, only: h10e5, capa, epsq, d00, elocp
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      implicit none
!
!-----------------------------------------------------------------------
!  LIST OF VARIABLES NEEDED
!    PARAMETERS:
!      INCLUDED IN "cuparm" AND "parm.tbl"
!    INPUT:
!      T,Q,HTM,PMID(3-D),KARR (2-D)
!    OUTPUT: 
!      TWET (3-D)
!    SUBROUTINES CALLED:
!      TTBLEX
!
      real,dimension(ista_2l:iend_2u,jsta_2l:jend_2u,LM),intent(in)  :: t,q,    &
                                                           pmid,htm
      integer,dimension(ista:iend,jsta:jend),      intent(in)  :: karr
      real,dimension(ista_2l:iend_2u,jsta_2l:jend_2u,LM),intent(out) :: twet


      real,    dimension(ista:iend,jsta:jend) :: thesp, qq, pp
      integer, dimension(ista:iend,jsta:jend) :: klres,khres,iptb,ithtb
!
      integer i,j,l,ittb1,ittbk,iqtbk,it,knuml,knumh,iq
      real tbtk,qbtk,apebtk,tthbtk,tthk,qqk,bqs00k,sqs00k,bqs10k,      &
           sqs10k,bqk,sqk,tqk,ppk,tpspk,apespk,prespk,p00k,p10k,p01k,  &
           p11k,presk
!
!--------------COMPUTE WET BULB TEMPERATURES----------------------------
!!$omp  parallel do
!!$omp& private(apebtk,apespk,bqk,bqs00k,bqs10k,iq,iqtbk,it,ittb1,ittbk,
!!$omp&         karr,khres,klres,knumh,knuml,p00k,p01k,p10k,p11k,ppk,
!!$omp&         presk,qbtk,qqk,sqk,sqs00k,sqs10k,tbtk,thesp,tpspk,
!!$omp&         tqk,tthbtk,tthk)
!-----------------------------------------------------------------------
      DO 300 l=1,lm
      DO 125 j=jsta,jend
      DO 125 i=ista,iend
        IF (htm(i,j,l)<1.0) THEN
          thesp(i,j)=273.15
          cycle    
        ENDIF
        IF(t(i,j,l)<spval)THEN
        tbtk  =t(i,j,l)
        qbtk  =q(i,j,l)
        presk =pmid(i,j,l)
        apebtk=(h10e5/presk)**capa
        IF(qbtk<epsq) qbtk=htm(i,j,l)*epsq
!--------------SCALING POTENTIAL TEMPERATURE & TABLE INDEX--------------
        tthbtk  =tbtk*apebtk
        tthk    =(tthbtk-thl)*rdth
        qqk     =tthk-aint(tthk)
        ittb1   =int(tthk)+1
!--------------KEEPING INDICES WITHIN THE TABLE-------------------------
        IF(ittb1<1) THEN
          ittb1  =1
          qqk    =d00
        ENDIF
!
        IF(ittb1>=jtb) THEN
        ittb1  =jtb-1
        qqk    =d00
        ENDIF
!--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY---------------
        ittbk=ittb1
        bqs00k=qs0(ittbk)
        sqs00k=sqs(ittbk)
        bqs10k=qs0(ittbk+1)
        sqs10k=sqs(ittbk+1)
!--------------SCALING SPEC. HUMIDITY & TABLE INDEX---------------------
        bqk=(bqs10k-bqs00k)*qqk+bqs00k
        sqk=(sqs10k-sqs00k)*qqk+sqs00k
        tqk=(qbtk-bqk)/sqk*rdq
        ppk=tqk-aint(tqk)
        iqtbk=int(tqk)+1
!--------------KEEPING INDICES WITHIN THE TABLE-------------------------
        IF(iqtbk<1) THEN
          iqtbk =1
          ppk   =d00
        ENDIF
!
        IF(iqtbk>=itb) THEN
          iqtbk=itb-1
          ppk  =d00
        ENDIF
!--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.-------
        iq=iqtbk
        it=ittb1
        p00k=ptbl(iq  ,it  )
        p10k=ptbl(iq+1,it  )
        p01k=ptbl(iq  ,it+1)
        p11k=ptbl(iq+1,it+1)
!--------------SATURATION POINT VARIABLES AT THE BOTTOM-----------------
        tpspk=p00k+(p10k-p00k)*ppk+(p01k-p00k)*qqk        &
                +(p00k-p10k-p01k+p11k)*ppk*qqk
        apespk=(h10e5/tpspk)**capa
        thesp(i,j)=tthbtk*exp(elocp*qbtk*apespk/tthbtk)
        ELSE
        thesp(i,j)=spval
        ENDIF !end t(i,j,l)<spval
!      ENDIF
  125 CONTINUE
!--------------SCALING PRESSURE & TT TABLE INDEX------------------------
      knuml=0
      knumh=0
!
      DO 280 j=jsta,jend
      DO 280 i=ista,iend
      klres(i,j)=0
      khres(i,j)=0
!
!      IF(KARR(I,J)>0)THEN
        IF(pmid(i,j,l)==spval)cycle
        presk=pmid(i,j,l)
!
        IF(presk<plq)THEN
          knuml=knuml+1
          klres(i,j)=1
        ELSE
          knumh=knumh+1
          khres(i,j)=1
        ENDIF
!      ENDIF
 280  CONTINUE
!***
!***  COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE<PL
!**
      IF(knuml>0)THEN
        CALL ttblex(twet(ista_2l,jsta_2l,l),ttbl,itb,jtb,klres     &
      ,pmid(ista_2l,jsta_2l,l),pl,qq,pp,rdp,the0,sthe              &
      ,rdthe,thesp,iptb,ithtb)
      ENDIF
!***
!***  COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PL
!**
      IF(knumh>0)THEN
       CALL ttblex(twet(ista_2l,jsta_2l,l),ttblq,itbq,jtbq,khres   &
      ,pmid(ista_2l,jsta_2l,l),plq,qq,pp,rdpq,the0q,stheq          &
      ,rdtheq,thesp,iptb,ithtb)
      ENDIF
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
  300 CONTINUE
      RETURN
      END