c======================================================================
c
c !IROUTINE: WETBULBCALC -- Calculate wet bulb temperature (C)
c
c !DESCRIPTION:
c
c   Calculates wet bulb temperature in C, given pressure in 
c      temperature in K and mixing ratio in kg/kg.
c
c !INPUT:
c    nx     - index for x dimension
c    ny     - index for y dimension
c    nz     - index for z dimension
c    prs    - pressure (mb)
c    tmk    - temperature (K)
c    qvp    - water vapor mixing ratio (kg/kg)
c
c !OUTPUT:
c    twb    - Wet bulb temperature (C)
c
c !ASSUMPTIONS:
c
c !REVISION HISTORY:
c     2009-March  - Mark T. Stoelinga - from RIP4.5
c     2010-August - J. Schramm
c     2014-March - A. Jaye - modified to run with NCL and ARW wrf output
c
c !INTERFACE:
c ------------------------------------------------------------------
C NCLFORTSTART
      subroutine wetbulbcalc(prs,tmk,qvp,twb,nx,ny,nz,psafile)
      implicit none
      integer nx, ny, nz
      double precision prs(nz,ny,nx)
      double precision tmk(nz,ny,nx)
      double precision qvp(nz,ny,nx)
      double precision twb(nz,ny,nx)
      character*(*) psafile
C NCLEND
      integer i,j,k
      integer jtch,jt,ipch,ip
      double precision q, t, p, e, tlcl, eth
      double precision fracip,fracip2,fracjt,fracjt2
      double precision PSADITHTE(150),PSADIPRS(150),PSADITMK(150,150)
      double precision tonpsadiabat
      double precision eps,tlclc1,tlclc2,tlclc3,tlclc4,gamma
      double precision gammamd,thtecon1,thtecon2,thtecon3,celkel
      double precision rgas,rgasmd,cp,cpmd

c
c  Before looping, set lookup table for getting temperature on
c  a pseudoadiabat.
c
      CALL DLOOKUP_TABLE(PSADITHTE,PSADIPRS,PSADITMK,psafile)

c  Define constants

      rgas=287.04  !J/K/kg
      rgasmd=.608   ! rgas_moist=rgas*(1.+rgasmd*qvp)
      cp=1004.     ! J/K/kg  Note: not using Bolton's value of 1005.7
      cpmd=.887   ! cp_moist=cp*(1.+cpmd*qvp)
      eps=0.622
      tlclc1=2840.
      tlclc2=3.5
      tlclc3=4.805
      tlclc4=55.
      gamma=rgas/cp
      gammamd=rgasmd-cpmd  ! gamma_moist=gamma*(1.+gammamd*qvp)
      thtecon1=3376. ! K
      thtecon2=2.54
      thtecon3=.81
      celkel=273.15

      DO k=1,nx
        DO j=1,ny
          DO i=1,nz
            q=dmax1(qvp(i,j,k),1.d-15)
            t=tmk(i,j,k)
            p=prs(i,j,k)/100.
            e=q*p/(eps+q)
            tlcl=tlclc1/(dlog(t**tlclc2/e)-tlclc3)+tlclc4
            eth=t*(1000./p)**(gamma*(1.+gammamd*q))*
     &         exp((thtecon1/tlcl-thtecon2)*q*(1.+thtecon3*q))


c
c   Now we need to find the temperature (in K) on a moist adiabat
c   (specified by eth in K) given pressure in hPa.  It uses a
c   lookup table, with data that was generated by the Bolton (1980)
c   formula for theta_e.
c
c     First check if pressure is less than min pressure in lookup table.
c     If it is, assume parcel is so dry that the given theta-e value can
c     be interpretted as theta, and get temperature from the simple dry
c     theta formula.
c

            if (p.le.psadiprs(150)) then
              tonpsadiabat=eth*(p/1000.)**gamma
            else   
c
c   Otherwise, look for the given thte/prs point in the lookup table.
c
            do jtch=1,150-1
              if (eth.ge.psadithte(jtch).and.eth.lt.
     &              psadithte(jtch+1)) then
                 jt=jtch
                 goto 213
              endif
            enddo
            jt=-1
 213        continue
            do ipch=1,150-1
              if (p.le.psadiprs(ipch).and.p.gt.psadiprs(ipch+1)) then
                 ip=ipch
                 goto 215
              endif
            enddo
            ip=-1
 215        continue
            if (jt.eq.-1.or.ip.eq.-1) then
               print*,
     &           'Outside of lookup table bounds. prs,thte=',p,eth
              stop
            endif
            fracjt=(eth-psadithte(jt))/(psadithte(jt+1)-psadithte(jt))
            fracjt2=1.-fracjt
            fracip=(psadiprs(ip)-p)/(psadiprs(ip)-psadiprs(ip+1))
            fracip2=1.-fracip
            if (psaditmk(ip,jt).gt.1e9.or.psaditmk(ip+1,jt).gt.1e9.or.
     &          psaditmk(ip,jt+1).gt.1e9.or.
     &          psaditmk(ip+1,jt+1).gt.1e9) then
                print*,
     &            'Tried to access missing tmperature in lookup table.'
                print*,
     &            'Prs and Thte probably unreasonable. prs,thte='
     &                   ,p,eth
               stop
            endif
            tonpsadiabat=fracip2*fracjt2*psaditmk(ip  ,jt  )+
     &             fracip *fracjt2*psaditmk(ip+1,jt  )+
     &             fracip2*fracjt *psaditmk(ip  ,jt+1)+
     &             fracip *fracjt *psaditmk(ip+1,jt+1)
            endif

            twb(i,j,k)=tonpsadiabat

          ENDDO
        ENDDO
      ENDDO

c
      return
      end
c======================================================================
c
c !IROUTINE: omgcalc -- Calculate omega (dp/dt)
c
c !DESCRIPTION:
c
c   Calculate approximate omega, based on vertical velocity w (dz/dt).
c   It is approximate because it cannot take into account the vertical
c   motion of pressure surfaces.
c
c !INPUT:
c    mx - index for x dimension
c    my - index for y dimension
c    mx -  index for vertical dimension
c    qvp - water vapor mixing ratio (kg/kg)
c    tmk - temperature (K)
c    www - vertical velocity (m/s)
c    prs -  pressure (Pa)
c
c !OUTPUT:
c    omg - omega (Pa/sec)
c
c !ASSUMPTIONS:
c
c !REVISION HISTORY:
c     2009-March  - Mark T. Stoelinga - from RIP4.5
c     2010-August - J. Schramm
c     2014-March - A. Jaye - modified to run with NCL and ARW wrf output
c
c ------------------------------------------------------------------
c NCLFORTSTART
      subroutine omgcalc(qvp,tmk,www,prs,omg,mx,my,mz)
      implicit none
      integer mx, my, mz
      double precision qvp(mz,my,mx)
      double precision tmk(mz,my,mx)
      double precision www(mz,my,mx)
      double precision prs(mz,my,mx) 
      double precision omg(mz,my,mx)
c NCLEND
c Local variables
      integer i, j, k
      double precision grav,rgas,eps 
c
c Constants
c
      grav=9.81           ! m/s**2
      rgas=287.04  !J/K/kg
      eps=0.622

      do k=1,mx
        do j=1,my
          do i=1,mz
            omg(i,j,k)=-grav*prs(i,j,k)/
     &        (rgas*((tmk(i,j,k)*(eps+qvp(i,j,k)))/
     &        (eps*(1.+qvp(i,j,k)))))*www(i,j,k)
          enddo
        enddo
      enddo
c
      return
      end
c======================================================================
c
c !IROUTINE: VIRTUAL_TEMP -- Calculate virtual temperature (K)
c
c !DESCRIPTION:
c
c   Calculates virtual temperature in K, given temperature
c      in K and mixing ratio in kg/kg.
c
c !INPUT:
c    NX     - index for x dimension
c    NY     - index for y dimension
c    NZ     - index for z dimension
c    RATMIX - water vapor mixing ratio (kg/kg)
c    TEMP   - temperature (K)
c
c !OUTPUT:
c    TV     - Virtual temperature (K)
c
c !ASSUMPTIONS:
c
c !REVISION HISTORY:
c     2009-March  - Mark T. Stoelinga - from RIP4.5
c     2010-August - J. Schramm
c     2014-March - A. Jaye - modified to run with NCL and ARW wrf output
c
c ------------------------------------------------------------------
C NCLFORTSTART
      SUBROUTINE VIRTUAL_TEMP(TEMP,RATMIX,TV,NX,NY,NZ)
      IMPLICIT NONE
      INTEGER NX,NY,NZ
      DOUBLE PRECISION TEMP(NZ,NY,NX)
      DOUBLE PRECISION RATMIX(NZ,NY,NX)
      DOUBLE PRECISION TV(NZ,NY,NX)
C NCLEND
      INTEGER I,J,K
      DOUBLE PRECISION EPS
      EPS = 0.622D0
      DO K=1,NX
        DO J=1,NY
          DO I=1,NZ
            TV(I,J,K) = TEMP(I,J,K)* (EPS+RATMIX(I,J,K))/ 
     &                     (EPS* (1.D0+RATMIX(I,J,K)))
          ENDDO
        ENDDO
      ENDDO
      RETURN
      END
