!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                   !!
!!                   GNU General Public License                      !!
!!                                                                   !!
!! This file is part of the Flexible Modeling System (FMS).          !!
!!                                                                   !!
!! FMS is free software; you can redistribute it and/or modify       !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation.             !!
!!                                                                   !!
!! FMS is distributed in the hope that it will be useful,            !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of    !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the     !!
!! GNU General Public License for more details.                      !!
!!                                                                   !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to:                                 !!
!!          Free Software Foundation, Inc.                           !!
!!          59 Temple Place, Suite 330                               !!
!!          Boston, MA  02111-1307  USA                              !!
!! or see:                                                           !!
!!          http://www.gnu.org/licenses/gpl.txt                      !!
!!                                                                   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
!---------------------------------------------------------------------
!------------ FMS version number and tagname for this file -----------
        
! $Id: array_lib.f90,v 1.1.2.1.2.1 2009/08/10 10:48:13 rsh Exp $
! $Name: mom4p1_pubrel_dec2009_nnz $

! ARRAY_LIB: Array procedures for F90
! Compiled/Modified:
!   07/01/06  John Haynes (haynes@atmos.colostate.edu)
!
! infind (function)
! lin_interpolate (function)
  
  module array_lib
  implicit none

  contains

! ----------------------------------------------------------------------------
! function INFIND
! ----------------------------------------------------------------------------
  function infind(list,val,sort,dist)
  use m_mrgrnk
  implicit none
!
! Purpose:
!   Finds the index of an array that is closest to a value, plus the
!   difference between the value found and the value specified
!
! Inputs:
!   [list]   an array of sequential values
!   [val]    a value to locate
! Optional input:
!   [sort]   set to 1 if [list] is in unknown/non-sequential order
!
! Returns:
!   index of [list] that is closest to [val]
!
! Optional output:
!   [dist]   set to variable containing [list([result])] - [val]
!
! Requires:
!   mrgrnk library
!
! Created:
!   10/16/03  John Haynes (haynes@atmos.colostate.edu)
! Modified:
!   01/31/06  IDL to Fortran 90
 
! ----- INPUTS -----
  real*8, dimension(:), intent(in) :: list
  real*8, intent(in) :: val  
  integer, intent(in), optional :: sort
  
! ----- OUTPUTS -----
  integer*4 :: infind
  real*8, intent(out), optional :: dist

! ----- INTERNAL -----
  real*8, dimension(size(list)) :: lists
  integer*4 :: nlist, result, tmp(1), sort_list
  integer*4, dimension(size(list)) :: mask, idx

  if (present(sort)) then
    sort_list = sort
  else
    sort_list = 0
  endif  

  nlist = size(list)
  if (sort_list == 1) then
    call mrgrnk(list,idx)
    lists = list(idx)
  else
    lists = list
  endif

  if (val >= lists(nlist)) then
    result = nlist
  else if (val <= lists(1)) then
    result = 1
  else
    mask(:) = 0
    where (lists < val) mask = 1
      tmp = minloc(mask,1)
      if (abs(lists(tmp(1)-1)-val) < abs(lists(tmp(1))-val)) then
        result = tmp(1) - 1
      else
        result = tmp(1)
      endif
  endif
  if (present(dist)) dist = lists(result)-val
  if (sort_list == 1) then
    infind = idx(result)
  else
    infind = result
  endif

  end function infind

! ----------------------------------------------------------------------------
! function LIN_INTERPOLATE
! ----------------------------------------------------------------------------  
  subroutine lin_interpolate(yarr,xarr,yyarr,xxarr,tol)
  use m_mrgrnk
  implicit none
!
! Purpose:
!   linearly interpolate a set of y2 values given a set of y1,x1,x2
!
! Inputs:
!   [yarr]    an array of y1 values
!   [xarr]    an array of x1 values
!   [xxarr]   an array of x2 values
!   [tol]     maximum distance for a match
!
! Output:
!   [yyarr]   interpolated array of y2 values
!
! Requires:
!   mrgrnk library
!
! Created:
!   06/07/06  John Haynes (haynes@atmos.colostate.edu)

! ----- INPUTS -----
  real*8, dimension(:), intent(in) :: yarr, xarr, xxarr
  real*8, intent(in) :: tol

! ----- OUTPUTS -----
  real*8, dimension(size(xxarr)), intent(out) :: yyarr

! ----- INTERNAL -----
  real*8, dimension(size(xarr)) :: ysort, xsort
  integer*4, dimension(size(xarr)) :: ist
  integer*4 :: nx, nxx, i, iloc
  real*8 :: d, m

  nx = size(xarr)
  nxx = size(xxarr)

! // xsort, ysort are sorted versions of xarr, yarr  
  call mrgrnk(xarr,ist)
  ysort = yarr(ist)
  xsort = xarr(ist)
  
  do i=1,nxx
    iloc = infind(xsort,xxarr(i),dist=d)
    if (d > tol) then
      print *, 'interpolation error'
      stop
    endif
    if (iloc == nx) then
!     :: set to the last value
      yyarr(i) = ysort(nx)
    else
!     :: is there another closeby value?
      if (abs(xxarr(i)-xsort(iloc+1)) < 2*tol) then
!       :: yes, do a linear interpolation      
        m = (ysort(iloc+1)-ysort(iloc))/(xsort(iloc+1)-xsort(iloc))
        yyarr(i) = ysort(iloc) + m*(xxarr(i)-xsort(iloc))
      else
!       :: no, set to the only nearby value
        yyarr(i) = ysort(iloc)
      endif
    endif
  enddo
  
  end subroutine lin_interpolate

  end module array_lib