subroutine tpindex(pw,tw,pref,tref,LCOEF,MT,MP)

!==================================================================
!     
!     Purpose
!     -------
!     Interpolate K-coefficients to the given P,T.
!
!     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 CO2 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
!     Pref(NP)           - The pressure grid array
!     Tref(NT)           - The temperature grid array
!   
!     Outputs
!     -------
!     TI                 - Interpolation term (pressure)
!     UI                 - Interpolation term (temperature)
!     MT                 - Temperature index (bottom left temperature)
!                          of bounding box
!     MP                 - Pressure index (bottom left pressure)
!                          of bounding box
!
!     Authors
!     -------
!     Adapted from the NASA Ames code by R. Wordsworth (2009)
!     Last update : B. de Batz de Trenquelléon (2023)
!     
!==================================================================

      use radinc_h

      implicit none

      ! INPUTS :
      !=========
      REAL,INTENT(IN) :: Tref(L_NTREF)
      REAL,INTENT(IN) :: pref(L_PINT)
      REAL,INTENT(IN) :: PW, TW

      ! OUTPUTS :
      !==========
      INTEGER,INTENT(OUT) :: MT, MP
      REAL,INTENT(OUT) :: LCOEF(4)

      ! LOCAL :
      !========
      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
      ENDIF

!     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

      end subroutine tpindex