module pindex_mod

implicit none

contains

subroutine pindex(pti,vti,pref,nvarlev,nlev,nv)

INTEGER, INTENT(IN) :: nvarlev
INTEGER, INTENT(IN) :: nlev
REAL, INTENT(IN) :: pti(nvarlev)
REAL, INTENT(IN) :: vti(nvarlev)
REAL, INTENT(IN) :: pref(nlev)

REAL, INTENT(OUT) :: nv(nlev)

REAL :: a,b
INTEGER :: n,l

!==================================================================
!     
!     Purpose
!     -------
!     Interpolate only data to the given P values.
!
!     Inputs
!     ------
!     pti                - The pressure to interpolate to
!     vti                - The variable to interpolate to
!     pref               - The tabulated pressure grid
!     nvarlev            - Number of levels of pti
!     nlev               - Number of levels of pref
!   
!     Outputs
!     -------
!     nv                 - Interpolated variable
!
!     Local variables
!     -------
!     n,l                - integer
!     m                  - coefficient for interpolation
!     k                  - coefficient for interpolation
!                        
!==================================================================

DO n=1,nlev
  IF(pti(1) .LT. pti(nvarlev)) THEN
    IF(pref(n).LE.pti(1)) THEN
      nv(n)=vti(1)
    ELSEIF(pref(n).GE.pti(nvarlev)) THEN
      nv(n)=vti(nvarlev)
    ELSE
      DO l=2,nvarlev
        IF(pti(l-1).LE.pref(n) .AND. pref(n).LE.pti(l)) THEN
          IF(abs(vti(l-1)-vti(l)).LT.1e-6) THEN
            nv(n) = (vti(l-1) + vti(l))/2.
          ELSE
            a = (pti(l)-pti(l-1))/(vti(l)-vti(l-1))
            b = (vti(l)*pti(l-1)-vti(l-1)*pti(l))/(vti(l)-vti(l-1))
            nv(n) = (pref(n)-b)/a
          ENDIF
        ENDIF
      ENDDO
    END IF
  ELSE
    IF(pref(n).LE.pti(nvarlev)) THEN
      nv(n)=vti(nvarlev)
    ELSEIF(pref(n).GE.pti(1)) THEN
      nv(n)=vti(1)
    ELSE
      DO l=2,nvarlev
        IF(pti(l-1).GE.pref(n) .AND. pref(n).GE.pti(l)) THEN
          IF(abs(vti(l-1)-vti(l)).LT.1e-6) THEN
            nv(n) = (vti(l-1) + vti(l))/2.
          ELSE
            a = (pti(l)-pti(l-1))/(vti(l)-vti(l-1))
            b = (vti(l)*pti(l-1)-vti(l-1)*pti(l))/(vti(l)-vti(l-1))
            nv(n) = (pref(n)-b)/a
          ENDIF
        ENDIF
      ENDDO
    END IF
  END IF
ENDDO

end subroutine pindex
end module pindex_mod