module tpindex_mod
      
      implicit none
      
      contains
      
      subroutine tpindex(pw,tw,qvar,pref,tref,wrefvar,LCOEF,MT,MP,
     &     NVAR,wratio)

!==================================================================
!     
!     Purpose
!     -------
!     Interpolate K-coefficients to the given P,T and Qvar values.
!
!     Notes
!     -----
!     The interpolation is the usual one in two dimensions given
!     in "Numerical Recipes", where the "X" are P, the "Y" are
!     T, and the F(X,Y) are the N2 K-coefficients.
!
!     The interpolating box is:
!
!           (PL,TU)                        (PR,TU)
!
!                          (TW,PW)
!
!           
!           (PL,TL)                        (PR,TL)
!
!      PL  - Pressure left
!      PR  - Pressure right
!      TL  - Temperature lower
!      TU  - Temperature upper
!      PW  - Pressure wanted
!      TW  - Temperature wanted
!
!     Inputs
!     ------
!     PW                 - The pressure to interpolate to
!     TW                 - The temperature to interpolate to
!     Qvar               - The abundance to interpolate to
!     Pref               - The tabulated (log) pressure grid
!     Tref               - The tabulated temperature grid
!     WREFVAR            - The tabulated abundance grid
!   
!     Outputs
!     -------
!     MT                 - Temperature index (bottom left temperature)
!                          of bounding box
!     MP                 - Pressure index (bottom left pressure)
!                          of bounding box
!     NVAR               -  index (bottom left )
!                          of bounding box
!     LCOEF(4)           - Interpolation coefficients
!     WRATIO             - weight for the interpolation of abundance
!
!     Authors
!     -------
!     Adapted from the NASA Ames code by R. Wordsworth (2009)
!     
!==================================================================

      use radinc_h, only: L_PINT, L_NTREF, L_REFVAR

      implicit none

      REAL,INTENT(IN) :: Tref(L_NTREF)
      REAL,INTENT(IN) :: pref(L_PINT)
      REAL,INTENT(IN) :: wrefvar(L_REFVAR)
      REAL,INTENT(IN) :: PW, TW, Qvar

      INTEGER,INTENT(OUT) :: MT, MP, NVAR
      REAL,INTENT(OUT) :: LCOEF(4), wratio
      
      INTEGER :: N
      REAL :: PWL ! local value for log10(PW)
      REAL :: TWL ! local value for TW
      REAL :: T ! linear interpolation weight along log(pressure)
      REAL :: U ! linear interpolation weight along temperature

!     1. For the interpolation along temperature 
!     Get the upper and lower temperature grid indicies that bound the
!     requested temperature. If the sought temperature is below
!     the T-grid minimum, assume that lower value. If the sought temperature
!     is above the T-grid maximum, assume that maximum value.
!     TW : temperature to interpolate to
!     TREF : reference temperature array
!     MT : index of TREF for (lower) bounding temperature
!     U : weight for interpolation in temperature

      TWL = TW

      IF(TWL.LT.TREF(1)) THEN
        ! Below lowest tabulated temperature
        MT = 1
        TWL=TREF(1)*1.00
      ELSEIF(TWL.GE.TREF(L_NTREF)) THEN
       ! Above highest tabulated temperature
        MT=L_NTREF-1
        TWL = TREF(L_NTREF)*1.00
      ELSE 
        ! Regular case, find encompassing tabulated temperatures
        do n=1,L_NTREF-1
          IF (TWL.GE.TREF(N) .and. TWL.LT.TREF(N+1)) THEN
            MT = n
            EXIT
          END IF
        END DO
      
      END IF

!     weight for the temperature interpolation
      U = (TWL-TREF(MT))/(TREF(MT+1)-TREF(MT))

!     2. For interpolation along pressure 
!     Get the upper and lower pressure grid indicies that bound the
!     requested pressure. If the sought pressure is below
!     the P-grid minimum, assume that lower value. If the sought pressure
!     is above the P-grid maximum, assume that maximum value.

      PWL = log10(PW)

      IF(PWL.LT.PREF(1)) THEN
        ! Below lowest tabulated pressure
        MP = 1
        PWL=Pref(1)*1.00
      ELSEIF(PWL.GE.PREF(L_PINT)) THEN
        ! Above highest tabulated pressure
        MP = L_PINT-1
        PWL = PREF(L_PINT)*1.00
      ELSE 
        ! Regular case, find encompassing tabulated pressures
        do n=1,L_PINT-1
          IF (PWL.GE.PREF(N) .and. PWL.LT.PREF(N+1)) THEN
            MP = n
            EXIT
          END IF
        END DO
      
      END IF

!     weight for the pressure interpolation
      T = (PWL-PREF(MP))/(PREF(MP+1)-PREF(MP))

!  Build the interpolation coefficients (bilinear in temperature-log(pressure))
      LCOEF(1) = (1.0-T)*(1.0-U)
      LCOEF(2) = T*(1.0-U)
      LCOEF(3) = T*U
      LCOEF(4) = (1.0-T)*U

!  3. For interpolation along abundances 
!  Get the indicies for abundance of the varying species. There are 10 sets of 
!  k-coefficients with differing amounts of variable vs. constant gas.

      IF(QVAR.LE.WREFVAR(1)) THEN
        ! Below lowest tabulated abundance
        NVAR   = 1
        WRATIO = 0.0D0      ! put all the weight on the first point
      ELSEIF(QVAR.GE.WREFVAR(L_REFVAR)) THEN
        ! Above highest tabulated abundance
        NVAR   = L_REFVAR-1 ! TB16 in order to not oversize NVAr when doing
                                  !NVAR+1
        WRATIO = 1.00D0     ! put all the weight on the last point
      ELSE
        ! Regular case, find encompassing tabulated abundances
        ! and compute corresponding weight
        DO N=2,L_REFVAR
          IF(QVAR.GE.WREFVAR(N-1) .and. QVAR.LT.WREFVAR(N)) THEN
            NVAR   = N-1
            WRATIO = (QVAR - WREFVAR(N-1))/(WREFVAR(N) - WREFVAR(N-1))
            EXIT
          END IF
        END DO
      END IF

      end subroutine tpindex
      
      end module tpindex_mod