module htinterp
  use shr_kind_mod ,only: r8 => shr_kind_r8
  use cam_logfile  ,only: iulog
  use cam_abortutils   ,only: endrun
  implicit none
  contains
!-----------------------------------------------------------------------
  subroutine cuthtint(fin,fht,idim1,nzp,fout,hts,nhts,logint,extrp,ier,iprnt)
    implicit none
!
! Do linear or log interpolation of fin to constant height surfaces
!   in hts(nhts). Interpolated field returned in fout.
!
! On input:
!   fin(idim1,nzp) = input field
!   fht(idim1,nzp) = heights field
!   idim1,nzp = dimensions of fin and fht
!   hts(nhts) = heights at which to interpolate
!   logint <= 0 -> do linear interpolation,
!                  otherwise do exponential (log) interpolation
!
! Nov, 2014 B.Foster: spval is no longer input, is set as parameter=zero
!
!   spval: If spval is non-zero and hts(k) is not within fht range
!          at the current grid point, then fout gets spval at the
!          current grid point.
!          If spval=0 and hts(k) is out of range, then
!          fout will be extrapolated to hts(k)
!   iprnt = 1 -> some error msgs will be printed if necessary
!
! On output:
!   fout(idim1,nhts) = contains interpolated field
!   ier = non-zero if an error has occurred
!
! Args:
      integer,intent(in) :: idim1,nzp,nhts,logint,extrp,iprnt
      real(r8),intent(in) :: fin(idim1,nzp),fht(idim1,nzp),hts(nhts)
      real(r8),intent(out) :: fout(idim1,nhts)
      integer,intent(out) :: ier
!
! Local:
      real(r8),parameter :: spval = 0.
      real(r8) :: col(nzp),htmin,htmax,exparg,f1
      integer :: n,k,k0,k1,ierr

      ier = 0
      do n=1,idim1
        col(:) = fht(n,:)
        htmin = minval(col)
        htmax = maxval(col)
        do k=1,nhts
          if (hts(k).lt.htmin.or.hts(k).gt.htmax) then
            if (extrp.ne.1) then
            ! fout(n,k) = spval

!             write(iulog,"('cuthtint: k=',i5,' hts(k)=',1pe12.4,' out of range: htmin,max=',2(1pe12.4),' fin(n,nzp)=',1pe12.4)") &
!              k,hts(k),htmin,htmax,fin(n,nzp) 

              fout(n,k) = fin(n,nzp)
              goto 100
            else
              k0 = nzp-1
              k1 = nzp
              write(*,*)'k0,k1',k0,k1,fin(n,k0),fin(n,k1)
            endif
            if (iprnt.gt.0) write(iulog,&
              "('>>> cuthtint: height ',f9.2,' out of range: htmin,max=',2f9.2,' k=',i5,' n=',i5)") &
                hts(k),htmin,htmax,k,n
          else
            call bracket(hts(k),col,nzp,1,k0,k1,ierr)
            if (ierr.ne.0) then
!             if (iprnt.gt.0) write(iulog,&
!               "('>>> cuthtint: bracket error: n=',i5,' k=',i5,' hts(k)=',f12.3,' hts col=',/(6f12.3))") &
!                 n,k,hts(k),col
            ! fout(n,k) = spval
!
! If the output coordinate cannot be bracketed in the input coordinate,
! then assign the input data. For example, time3d heights to higher (600 km) 
! than waccm heights (~400 km), so when interpolating waccm data to time3d
! heights, use the time3d data above the top of waccm.
!
              fout(n,k) = fin(n,nzp)
            ! ier = 1
              goto 100
            endif
          endif
!
! Do linear or log interpolation:
!
          if (logint.gt.0) then
            if (fin(n,k1).gt.0..and.fin(n,k1).ne.spval.and.&
                fin(n,k0).gt.0..and.fin(n,k0).ne.spval) then
              exparg = (log(fin(n,k1) / fin(n,k0)) / &
                       (col(k1)-col(k0))) * (hts(k)-col(k0))
              fout(n,k) = fin(n,k0)*exp(exparg)
            else
            ! fout(n,k) = spval
              fout(n,k) = fin(n,nzp)
            endif
          else
            if (fin(n,k0).ne.spval.and.fin(n,k1).ne.spval) then
              f1 = (hts(k)-col(k0)) / (col(k1)-col(k0))
              fout(n,k) = f1*fin(n,k1) + (1.-f1)*fin(n,k0)
            else
            ! fout(n,k) = spval
              fout(n,k) = fin(n,nzp)
            endif
          endif
 100      continue
        enddo   ! k=1,nhts
      enddo     ! n=1,idim1
      return
      end subroutine cuthtint
!-------------------------------------------------------------------
!     subroutine interp_column(fin,levin,nlevin,fout,levout,nlevout,kout,&
!       logint,iprint)
      subroutine interp_column(fin,levin,nlevin,fout,levout,nlevout,&
        logint,iprint)
!
! Interpolate column fin (at levin coord) to column fout (at levout coord).
! If logint > 0, do log interpolation, otherwise linear.
!
! Args:
      integer,intent(in) :: nlevin,nlevout,logint,iprint
!     integer,intent(out) :: kout
      real(r8),intent(in) :: fin(nlevin),levin(nlevin),levout(nlevout)
      real(r8),intent(out) :: fout(nlevout)
!
! Local:
      integer :: k,inc,icount,k0,k1,ier
      logical :: increasing
      real(r8) :: exparg,f1
!
! levin must be monotonic, but can be increasing or decreasing (1 to nlevin).
!
      increasing = .true.
      icount = 0
      do k=1,nlevin-1
        if (levin(k) > levin(k+1)) icount = icount+1
      enddo
      if (icount == nlevin-1) then
        increasing = .false.
        icount = 0
      endif
      if (icount > 0) &
        call endrun('>>> interp_column: Non-monotonic input levin')
      inc = 0
      if (increasing) inc = 1

      if (iprint > 0) then
        write(iulog,"('Enter interp_column: nlevin=',i4,' inc=',i2,' levin=',/,(6e12.4))") &
          nlevin,inc,levin
        write(iulog,"('Enter interp_column: fin=',/,(6e12.4))") fin
        write(iulog,"('Enter interp_column: nlevout=',i4,' levout=',/,(6e12.4))") &
          nlevout,levout
      endif

!     kout = nlevout
      do k=1,nlevout
        call bracket(levout(k),levin,nlevin,inc,k0,k1,ier)

        if (ier /= 0) then
          write(iulog,"('>>> interp_column: error from bracket=',i4,' k=',i4)") ier,k
          write(iulog,"('    Could not find levout(k)=',e12.4,' in levin=',/,(6e12.4))") &
            levout(k),levin
          call endrun('>>> interp_column: ier from bracket')
        endif
!
! Return k-index (kout) level at which last good interpolation occurred:
!       if (ier /= 0) then
!         if (iprint > 0) then
!           write(iulog,"('>>> interp_column: error from bracket=',i4,' k=',i4)") ier,k
!           write(iulog,"('    Could not find levout(k)=',e12.4,' in levin=',/,(6e12.4))") &
!             levout(k),levin
!         endif
!         if (k==1) then
!           call endrun('>>> interp_column: ier from bracket')
!         endif
!         kout = k-1
!         if (iprint > 0) then
!           write(iulog,"('interp_column returning: kout=',i4,' fout=',/,(6e12.4))") &
!             kout,fout(1:kout)
!         endif
!         return
!       endif

        if (logint > 0) then ! log interp
          exparg = (log(fin(k1) / fin(k0)) / &
                   (levin(k1)-levin(k0))) * (levout(k)-levin(k0))
          fout(k) = fin(k0)*exp(exparg)
        else                 ! linear interp
          f1 = (levout(k)-levin(k0)) / (levin(k1)-levin(k0))
          fout(k) = f1*fin(k1) + (1.-f1)*fin(k0)
        endif
      enddo ! k=1,nlevout

      if (iprint > 0) then
        write(iulog,"('interp_column returning: nlevout=',i4,' fout=',/,(6e12.4))") &
          nlevout,fout
      endif

      end subroutine interp_column
!-------------------------------------------------------------------
      subroutine bracket(x,xx,nx,inc,n1,n2,ier)
      implicit none
!
! Bracket x in xx(nx), returning lower index in n1 and upper index in n2
! If inc > 0 -> array increases from bottom to top,
! If inc <= 0 -> array increases from top to bottom
!
! Args:
      integer,intent(in) :: nx,inc
      integer,intent(out) :: n1,n2,ier
      real(r8),intent(in) :: x,xx(nx)
!
! Local:
      integer :: i
!
! Array increases from bottom to top:
      ier = 0
      if (inc.gt.0) then
        if (x.lt.xx(1)) then
          n1 = 1
          n2 = 2
          ier = 1
          return
        endif
        if (x.gt.xx(nx)) then
          n1 = nx-1
          n2 = nx
          ier = 2
          return
        endif
        do i=1,nx-1
          if (x.ge.xx(i).and.x.le.xx(i+1)) then
            n1 = i
            n2 = i+1
            return
          endif
        enddo
!
! Array increases from top to bottom:
      else
        if (x.gt.xx(1)) then
          n1 = 1
          n2 = 2
          ier = 3
          return
        endif
        if (x.lt.xx(nx)) then
          n1 = nx-1
          n2 = nx
          ier = 4
          return
        endif
        do i=1,nx-1
          if (x.le.xx(i).and.x.ge.xx(i+1)) then
            n1 = i
            n2 = i+1
            return
          endif
        enddo
      endif
      return
      end subroutine bracket
!-------------------------------------------------------------------
end module htinterp