module horiz_interp_bilinear_mod

  ! <CONTACT EMAIL="Zhi.Liang@noaa.gov"> Zhi Liang </CONTACT>

  ! <HISTORY SRC="http://www.gfdl.noaa.gov/fms-cgi-bin/cvsweb.cgi/FMS/"/>

  ! <OVERVIEW>
  !   Performs spatial interpolation between grids using bilinear interpolation
  ! </OVERVIEW>

  ! <DESCRIPTION>
  !     This module can interpolate data from regular rectangular grid
  !     to rectangular/tripolar grid. The interpolation scheme is bilinear interpolation.
  !     There is an optional mask field for missing input data.
  !     An optional output mask field may be used in conjunction with
  !     the input mask to show where output data exists.
  ! </DESCRIPTION>

  use mpp_mod,               only: mpp_error, FATAL, stdout, mpp_pe, mpp_root_pe
  use fms_mod,               only: write_version_number
  use constants_mod,         only: PI
  use horiz_interp_type_mod, only: horiz_interp_type, stats

  implicit none
  private


  public :: horiz_interp_bilinear_new, horiz_interp_bilinear, horiz_interp_bilinear_del
  public :: horiz_interp_bilinear_init

  !--- public interface
  interface horiz_interp_bilinear_new
    module procedure horiz_interp_bilinear_new_1d
    module procedure horiz_interp_bilinear_new_2d
  end interface


  real, parameter :: epsln=1.e-10
  integer, parameter :: DUMMY = -999

  !-----------------------------------------------------------------------
  character(len=128) :: version = '$Id$'
  character(len=128) :: tagname = '$Name$'
  logical            :: module_is_initialized = .FALSE.

contains

  !#######################################################################
  !  <SUBROUTINE NAME="horiz_interp_bilinear_init">
  !  <OVERVIEW>
  !     writes version number and tag name to logfile.out
  !  </OVERVIEW>
  !  <DESCRIPTION>       
  !     writes version number and tag name to logfile.out
  !  </DESCRIPTION>

  subroutine horiz_interp_bilinear_init

    if(module_is_initialized) return
    call write_version_number (version, tagname)
    module_is_initialized = .true.

  end subroutine horiz_interp_bilinear_init

  !  </SUBROUTINE>

  !########################################################################

  subroutine horiz_interp_bilinear_new_1d ( Interp, lon_in, lat_in, lon_out, lat_out, &
       verbose, src_modulo )

    !-----------------------------------------------------------------------
    type(horiz_interp_type), intent(inout) :: Interp
    real, intent(in),  dimension(:)        :: lon_in , lat_in
    real, intent(in),  dimension(:,:)      :: lon_out, lat_out
    integer, intent(in),          optional :: verbose
    logical, intent(in),          optional :: src_modulo

    logical :: src_is_modulo
    integer :: nlon_in, nlat_in, nlon_out, nlat_out, n, m
    integer :: ie, is, je, js, ln_err, lt_err, warns, unit
    real    :: wtw, wte, wts, wtn, lon, lat, tpi, hpi
    real    :: glt_min, glt_max, gln_min, gln_max, min_lon, max_lon

    warns = 0
    if(present(verbose)) warns = verbose
    src_is_modulo = .true. 
    if (present(src_modulo)) src_is_modulo = src_modulo

    hpi = 0.5*pi
    tpi = 4.0*hpi
    glt_min = hpi
    glt_max = -hpi
    gln_min = tpi
    gln_max = -tpi
    min_lon = 0.0
    max_lon = tpi
    ln_err = 0
    lt_err = 0
    !-----------------------------------------------------------------------

    allocate ( Interp % wti (size(lon_out,1),size(lon_out,2),2),   &
               Interp % wtj (size(lon_out,1),size(lon_out,2),2),   &
               Interp % i_lon (size(lon_out,1),size(lon_out,2),2), &
               Interp % j_lat (size(lon_out,1),size(lon_out,2),2))
    !-----------------------------------------------------------------------

    nlon_in = size(lon_in(:))  ; nlat_in = size(lat_in(:))
    nlon_out = size(lon_out, 1); nlat_out = size(lon_out, 2)
    Interp%nlon_src = nlon_in;  Interp%nlat_src = nlat_in
    Interp%nlon_dst = nlon_out; Interp%nlat_dst = nlat_out

    if(src_is_modulo) then
       if(lon_in(nlon_in) - lon_in(1) .gt. tpi + epsln) &
            call mpp_error(FATAL,'horiz_interp_bilinear_mod: '// & 
            'The range of source grid longitude should be no larger than tpi')

       if(lon_in(1) .lt. 0.0 .OR. lon_in(nlon_in) > tpi ) then
          min_lon = lon_in(1)
          max_lon = lon_in(nlon_in)
       endif
    endif

    do n = 1, nlat_out
       do m = 1, nlon_out
          lon = lon_out(m,n)
          lat = lat_out(m,n)

          if(src_is_modulo) then
             if(lon .lt. min_lon) then
                lon = lon + tpi
             else if(lon .gt. max_lon) then
                lon = lon - tpi
             endif
          else  ! when the input grid is in not cyclic, the output grid should located inside
             ! the input grid
             if((lon .lt. lon_in(1)) .or. (lon .gt. lon_in(nlon_in))) &
                  call mpp_error(FATAL,'horiz_interp_bilinear_mod: ' //&
                  'when input grid is not modulo, output grid should locate inside input grid')
          endif

          glt_min = min(lat,glt_min);  glt_max = max(lat,glt_max)
          gln_min = min(lon,gln_min);  gln_max = max(lon,gln_max)

          is = indp(lon, lon_in ) 
          if( lon_in(is) .gt. lon ) is = max(is-1,1)
          if( lon_in(is) .eq. lon .and. is .eq. nlon_in) is = max(is - 1,1)
          ie = min(is+1,nlon_in)
          if(lon_in(is) .ne. lon_in(ie) .and. lon_in(is) .le. lon) then
             wtw = ( lon_in(ie) - lon) / (lon_in(ie) - lon_in(is) )
          else
             !     east or west of the last data value. this could be because a
             !     cyclic condition is needed or the dataset is too small. 
             ln_err = 1
             ie = 1
             is = nlon_in
             if (lon_in(ie) .ge. lon ) then
                wtw = (lon_in(ie) -lon)/(lon_in(ie)-lon_in(is)+tpi+epsln)
             else
                wtw = (lon_in(ie) -lon+tpi+epsln)/(lon_in(ie)-lon_in(is)+tpi+epsln)
             endif
          endif
          wte = 1. - wtw

          js = indp(lat, lat_in ) 

          if( lat_in(js) .gt. lat ) js = max(js - 1, 1)
          if( lat_in(js) .eq. lat .and. js .eq. nlat_in) js = max(js - 1, 1)
          je = min(js + 1, nlat_in)

          if ( lat_in(js) .ne. lat_in(je) .and. lat_in(js) .le. lat) then
             wts = ( lat_in(je) - lat )/(lat_in(je)-lat_in(js))
          else
             !     north or south of the last data value. this could be because a
             !     pole is not included in the data set or the dataset is too small.
             !     in either case extrapolate north or south
             lt_err = 1
             wts = 1.
          endif

          wtn = 1. - wts

          Interp % i_lon (m,n,1) = is; Interp % i_lon (m,n,2) = ie
          Interp % j_lat (m,n,1) = js; Interp % j_lat (m,n,2) = je
          Interp % wti   (m,n,1) = wtw
          Interp % wti   (m,n,2) = wte
          Interp % wtj   (m,n,1) = wts
          Interp % wtj   (m,n,2) = wtn

       enddo
    enddo

    unit = stdout()

    if (ln_err .eq. 1 .and. warns > 0) then
       write (unit,'(/,(1x,a))')                                      &
            '==> Warning: the geographic data set does not extend far   ', &
            '             enough east or west - a cyclic boundary       ', &
            '             condition was applied. check if appropriate   '
       write (unit,'(/,(1x,a,2f8.4))')                                &
            '    data required between longitudes:', gln_min, gln_max,     &
            '      data set is between longitudes:', lon_in(1), lon_in(nlon_in)
       warns = warns - 1
    endif

    if (lt_err .eq. 1 .and. warns > 0) then
       write (unit,'(/,(1x,a))')                                     &
            '==> Warning: the geographic data set does not extend far   ',&
            '             enough north or south - extrapolation from    ',&
            '             the nearest data was applied. this may create ',&
            '             artificial gradients near a geographic pole   ' 
       write (unit,'(/,(1x,a,2f8.4))')                             &
            '    data required between latitudes:', glt_min, glt_max,   &
            '      data set is between latitudes:', lat_in(1), lat_in(nlat_in)
    endif

    return

  end subroutine horiz_interp_bilinear_new_1d

  !#######################################################################
  ! <SUBROUTINE NAME="horiz_interp_bilinear_new">

  !   <OVERVIEW>
  !      Initialization routine.
  !   </OVERVIEW>
  !   <DESCRIPTION>
  !      Allocates space and initializes a derived-type variable
  !      that contains pre-computed interpolation indices and weights.
  !   </DESCRIPTION>
  !   <TEMPLATE>
  !     call horiz_interp_bilinear_new ( Interp, lon_in, lat_in, lon_out, lat_out, verbose, src_modulo )

  !   </TEMPLATE>
  !   
  !   <IN NAME="lon_in" TYPE="real, dimension(:,:)" UNITS="radians">
  !      Longitude (in radians) for source data grid. 
  !   </IN>

  !   <IN NAME="lat_in" TYPE="real, dimension(:,:)" UNITS="radians">
  !      Latitude (in radians) for source data grid.
  !   </IN>

  !   <IN NAME="lon_out" TYPE="real, dimension(:,:)" UNITS="radians" >
  !      Longitude (in radians) for source data grid. 
  !   </IN>

  !   <IN NAME="lat_out" TYPE="real, dimension(:,:)" UNITS="radians" >
  !      Latitude (in radians) for source data grid. 
  !   </IN>

  !   <IN NAME="src_modulo" TYPE="logical, optional">
  !      logical variable to indicate if the boundary condition along zonal boundary
  !      is cyclic or not. When true, the zonal boundary condition is cyclic.
  !   </IN>

  !   <IN NAME="verbose" TYPE="integer, optional" >
  !      flag for the amount of print output.
  !   </IN>

  !   <INOUT NAME="Interp" TYPE="type(horiz_interp_type)" >
  !      A derived-type variable containing indices and weights used for subsequent 
  !      interpolations. To reinitialize this variable for a different grid-to-grid 
  !      interpolation you must first use the "horiz_interp_del" interface.
  !   </INOUT>

  subroutine horiz_interp_bilinear_new_2d ( Interp, lon_in, lat_in, lon_out, lat_out, &
       verbose, src_modulo, new_search, no_crash_when_not_found )

    !-----------------------------------------------------------------------
    type(horiz_interp_type), intent(inout) :: Interp
    real, intent(in),  dimension(:,:)      :: lon_in , lat_in
    real, intent(in),  dimension(:,:)      :: lon_out, lat_out
    integer, intent(in),          optional :: verbose
    logical, intent(in),          optional :: src_modulo
    logical, intent(in),          optional :: new_search
    logical, intent(in),          optional :: no_crash_when_not_found
    integer                                :: warns 
    logical                                :: src_is_modulo
    integer                                :: nlon_in, nlat_in, nlon_out, nlat_out
    integer                                :: m, n, is, ie, js, je, num_solution
    real                                   :: lon, lat, quadra, x, y, y1, y2
    real                                   :: a1, b1, c1, d1, a2, b2, c2, d2, a, b, c
    real                                   :: lon1, lat1, lon2, lat2, lon3, lat3, lon4, lat4
    real                                   :: tpi, lon_min, lon_max
    real                                   :: epsln2
    logical                                :: use_new_search, no_crash

    tpi = 2.0*pi

    warns = 0
    if(present(verbose)) warns = verbose
    src_is_modulo = .true. 
    if (present(src_modulo)) src_is_modulo = src_modulo
    use_new_search = .false.
    if (present(new_search)) use_new_search = new_search
    no_crash = .false.
    if(present(no_crash_when_not_found)) no_crash = no_crash_when_not_found

    ! make sure lon and lat has the same dimension
    if(size(lon_out,1) /= size(lat_out,1) .or. size(lon_out,2) /= size(lat_out,2) ) &
         call mpp_error(FATAL,'horiz_interp_bilinear_mod: when using bilinear ' // &
         'interplation, the output grids should be geographical grids')    

    if(size(lon_in,1) /= size(lat_in,1) .or. size(lon_in,2) /= size(lat_in,2) ) &
         call mpp_error(FATAL,'horiz_interp_bilinear_mod: when using bilinear '// &
         'interplation, the input grids should be geographical grids')  

    !--- get the grid size 
    nlon_in  = size(lon_in,1) ; nlat_in  = size(lat_in,2)
    nlon_out = size(lon_out,1); nlat_out = size(lon_out,2)
    Interp%nlon_src = nlon_in;  Interp%nlat_src = nlat_in
    Interp%nlon_dst = nlon_out; Interp%nlat_dst = nlat_out

    allocate ( Interp % wti (size(lon_out,1),size(lon_out,2),2),   &
               Interp % wtj (size(lon_out,1),size(lon_out,2),2),   &
               Interp % i_lon (size(lon_out,1),size(lon_out,2),2), &
               Interp % j_lat (size(lon_out,1),size(lon_out,2),2))

    !--- first fine the neighbor points for the destination points.
    if(use_new_search) then
       epsln2 = epsln*1e5
       call find_neighbor_new(Interp, lon_in, lat_in, lon_out, lat_out, src_is_modulo, no_crash)
    else
       epsln2 = epsln
       call find_neighbor(Interp, lon_in, lat_in, lon_out, lat_out, src_is_modulo)
    endif

    !***************************************************************************
    !         Algorithm explanation (from disscussion with Steve Garner )      *
    !                                                                          *
    !    lon(x,y) = a1*x + b1*y + c1*x*y + d1         (1)                      *
    !    lat(x,y) = a2*x + b2*y + c2*x*y + d2         (2)                      *
    !    f (x,y) = a3*x + b3*y + c3*x*y + d3          (3)                      *
    !    with x and y is between 0 and 1.                                      *
    !    lon1 = lon(0,0) = d1,          lat1 = lat(0,0) = d2                   *
    !    lon2 = lon(1,0) = a1+d1,       lat2 = lat(1,0) = a2+d2                *
    !    lon3 = lon(1,1) = a1+b1+c1+d1, lat3 = lat(1,1) = a2+b2+c2+d2          *
    !    lon4 = lon(0,1) = b1+d1,       lat4 = lat(0,1) = b2+d2                *
    !    where (lon1,lat1),(lon2,lat2),(lon3,lat3),(lon4,lat4) represents      *
    !    the four corners starting from the left lower corner of grid box      *
    !    that encloses a destination grid ( the rotation direction is          *
    !    counterclockwise ). With these conditions, we get                     *
    !    a1 = lon2-lon1,           a2 = lat2-lat1                              *
    !    b1 = lon4-lon1,           b2 = lat4-lat1                              *
    !    c1 = lon3-lon2-lon4+lon1, c2 = lat3-lat2-lat4+lat1                    *
    !    d1 = lon1                 d2 = lat1                                   *
    !    So given any point (lon,lat), from equation (1) and (2) we can        *
    !    solve (x,y).                                                          *
    !    From equation (3)                                                     *
    !    f1 = f(0,0) = d3,          f2 = f(1,0) = a3+d3                        *
    !    f3 = f(1,1) = a3+b3+c3+d3, f4 = f(0,1) = b3+d3                        *
    !    we obtain                                                             *
    !    a3 = f2-f1,       b3 = f4-f1                                          *
    !    c3 = f3-f2-f4+f1, d3 = f1                                             *
    !    at point (lon,lat) ---> (x,y)                                         *
    !    f(x,y) = (f2-f1)x + (f4-f1)y + (f3-f2-f4+f1)xy + f1                   *
    !           = f1*(1-x)*(1-y) + f2*x*(1-y) + f3*x*y + f4*y*(1-x)            *
    !    wtw=1-x; wte=x; wts=1-y; xtn=y                                        *
    !                                                                          *
    !***************************************************************************

    lon_min = minval(lon_in);
    lon_max = maxval(lon_in);
    !--- calculate the weight
    do n = 1, nlat_out
       do m = 1, nlon_out
          lon = lon_out(m,n)
          lat = lat_out(m,n)
          if(lon .lt. lon_min) then
             lon = lon + tpi
          else if(lon .gt. lon_max) then
             lon = lon - tpi
          endif
          is = Interp%i_lon(m,n,1); ie = Interp%i_lon(m,n,2)
          js = Interp%j_lat(m,n,1); je = Interp%j_lat(m,n,2)
          if( is == DUMMY) cycle
          lon1 = lon_in(is,js); lat1 = lat_in(is,js);
          lon2 = lon_in(ie,js); lat2 = lat_in(ie,js);
          lon3 = lon_in(ie,je); lat3 = lat_in(ie,je);
          lon4 = lon_in(is,je); lat4 = lat_in(is,je); 
          if(lon .lt. lon_min) then
             lon1 = lon1 -tpi; lon4 = lon4 - tpi
          else if(lon .gt. lon_max) then
             lon2 = lon2 +tpi; lon3 = lon3 + tpi
          endif                      
          a1 = lon2-lon1
          b1 = lon4-lon1
          c1 = lon1+lon3-lon4-lon2
          d1 = lon1
          a2 = lat2-lat1
          b2 = lat4-lat1
          c2 = lat1+lat3-lat4-lat2
          d2 = lat1
          !--- the coefficient of the quadratic equation
          a  = b2*c1-b1*c2
          b  = a1*b2-a2*b1+c1*d2-c2*d1+c2*lon-c1*lat
          c  = a2*lon-a1*lat+a1*d2-a2*d1
          quadra = b*b-4.*a*c
          if(abs(quadra) < epsln) quadra = 0.0
          if(quadra < 0.0) call mpp_error(FATAL, &
               "horiz_interp_bilinear_mod: No solution existed for this quadratic equation")
          if ( abs(a) .lt. epsln2) then  ! a = 0 is a linear equation
             if( abs(b) .lt. epsln) call mpp_error(FATAL, &
                  "horiz_interp_bilinear_mod: no unique solution existed for this linear equation")
             y = -c/b
          else
             y1 = 0.5*(-b+sqrt(quadra))/a
             y2 = 0.5*(-b-sqrt(quadra))/a
             if(abs(y1) < epsln2) y1 = 0.0
             if(abs(y2) < epsln2) y2 = 0.0
             if(abs(1.0-y1) < epsln2) y1 = 1.0
             if(abs(1.0-y2) < epsln2) y2 = 1.0
             num_solution = 0
             if(y1 >= 0.0 .and. y1 <= 1.0) then
                y = y1
                num_solution = num_solution +1
             endif
             if(y2 >= 0.0 .and. y2 <= 1.0) then
                y = y2
                num_solution = num_solution + 1
             endif
             if(num_solution == 0) then
                call mpp_error(FATAL, "horiz_interp_bilinear_mod: No solution found")
             else if(num_solution == 2) then
                call mpp_error(FATAL, "horiz_interp_bilinear_mod: Two solutions found")
             endif
           endif
           if(abs(a1+c1*y) < epsln) call mpp_error(FATAL, &
               "horiz_interp_bilinear_mod: the denomenator is 0")
           if(abs(y) < epsln2) y = 0.0
           if(abs(1.0-y) < epsln2) y = 1.0
           x = (lon-b1*y-d1)/(a1+c1*y)
           if(abs(x) < epsln2) x = 0.0
           if(abs(1.0-x) < epsln2) x = 1.0
           ! x and y should be between 0 and 1.
           !! Added for ECDA 
           if(use_new_search) then
             if (x < 0.0) x = 0.0 ! snz
             if (y < 0.0) y = 0.0 ! snz
             if (x > 1.0) x = 1.0
             if (y > 1.0) y = 1.0
           endif 
           if( x>1. .or. x<0. .or. y>1. .or. y < 0.) call mpp_error(FATAL, &
                "horiz_interp_bilinear_mod: weight should be between 0 and 1")
           Interp % wti(m,n,1)=1.0-x; Interp % wti(m,n,2)=x   
           Interp % wtj(m,n,1)=1.0-y; Interp % wtj(m,n,2)=y          
       enddo
    enddo

  end subroutine horiz_interp_bilinear_new_2d
  ! </SUBROUTINE>

  !#######################################################################
  ! this routine will search the source grid to fine the grid box that encloses 
  ! each destination grid.
  subroutine find_neighbor( Interp, lon_in, lat_in, lon_out, lat_out, src_modulo )
    type(horiz_interp_type), intent(inout) :: Interp
    real, intent(in),       dimension(:,:) :: lon_in , lat_in
    real, intent(in),       dimension(:,:) :: lon_out, lat_out
    logical,                 intent(in)    :: src_modulo
    integer                                :: nlon_in, nlat_in, nlon_out, nlat_out
    integer                                :: max_step, n, m, l, i, j, ip1, jp1, step
    integer                                :: is, js, jstart, jend, istart, iend, npts
    integer, allocatable, dimension(:)     :: ilon, jlat
    real                                   :: lon_min, lon_max, lon, lat, tpi
    logical                                :: found
    real                                   :: lon1, lat1, lon2, lat2, lon3, lat3, lon4, lat4

    tpi = 2.0*pi
    nlon_in  = size(lon_in,1) ; nlat_in  = size(lat_in,2)
    nlon_out = size(lon_out,1); nlat_out = size(lon_out,2)

    lon_min = minval(lon_in);
    lon_max = maxval(lon_in);

    max_step = min(nlon_in,nlat_in)/2 ! can be adjusted if needed
    allocate(ilon(8*max_step), jlat(8*max_step) )

    do n = 1, nlat_out
       do m = 1, nlon_out
          found = .false.
          lon = lon_out(m,n)
          lat = lat_out(m,n)

          if(src_modulo) then
             if(lon .lt. lon_min) then
                lon = lon + tpi
             else if(lon .gt. lon_max) then
                lon = lon - tpi
             endif
          else
             if(lon .lt. lon_min .or. lon .gt. lon_max ) &
             call mpp_error(FATAL,'horiz_interp_bilinear_mod: ' //&
                  'when input grid is not modulo, output grid should locate inside input grid')
          endif
          !--- search for the surrounding four points locatioon.
          if(m==1 .and. n==1) then
             J_LOOP: do j = 1, nlat_in-1
                do i = 1, nlon_in
                   ip1 = i+1
                   jp1 = j+1
                   if(i==nlon_in) then
                      if(src_modulo)then
                         ip1 = 1
                      else
                         cycle
                      endif
                   endif
                   lon1 = lon_in(i,  j);   lat1 = lat_in(i,j)
                   lon2 = lon_in(ip1,j);   lat2 = lat_in(ip1,j)
                   lon3 = lon_in(ip1,jp1); lat3 = lat_in(ip1,jp1)
                   lon4 = lon_in(i,  jp1); lat4 = lat_in(i,  jp1)  

                   if(lon .lt. lon_min .or. lon .gt. lon_max) then
                      if(i .ne. nlon_in) then
                         cycle
                      else
                         if(lon .lt. lon_min) then
                             lon1 = lon1 -tpi; lon4 = lon4 - tpi
                         else if(lon .gt. lon_max) then
                             lon2 = lon2 +tpi; lon3 = lon3 + tpi
                         endif
                      endif
                   endif

                   if(lat .ge. intersect(lon1,lat1,lon2,lat2,lon))then ! south
                      if(lon .le. intersect(lat2,lon2,lat3,lon3,lat))then ! east
                         if(lat .le. intersect(lon3,lat3,lon4,lat4,lon))then ! north
                            if(lon .ge. intersect(lat4,lon4,lat1,lon1,lat))then  ! west
                               found = .true.
                               Interp % i_lon (m,n,1) = i; Interp % i_lon (m,n,2) = ip1
                               Interp % j_lat (m,n,1) = j; Interp % j_lat (m,n,2) = jp1
                               exit J_LOOP
                            endif
                         endif
                      endif
                   endif
                enddo
             enddo J_LOOP
          else
             step = 0
             do while ( .not. found .and. step .lt. max_step )
                !--- take the adajcent point as the starting point
                if(m == 1) then
                   is = Interp % i_lon (m,n-1,1)
                   js = Interp % j_lat (m,n-1,1)
                else
                   is = Interp % i_lon (m-1,n,1)
                   js = Interp % j_lat (m-1,n,1)
                endif
                if(step==0) then
                   npts = 1
                   ilon(1) = is
                   jlat(1) = js
                else
                   npts = 0
                   !--- bottom and top boundary
                   jstart = max(js-step,1)
                   jend   = min(js+step,nlat_in)

                   do l = -step, step
                      i = is+l
                      if(src_modulo)then
                         if( i < 1) then
                            i = i + nlon_in
                         else if (i > nlon_in) then
                            i = i - nlon_in
                         endif
                         if( i < 1 .or. i > nlon_in) call mpp_error(FATAL, &
                              'horiz_interp_bilinear_mod: max_step is too big, decrease max_step' )
                      else
                         if( i < 1 .or. i > nlon_in) cycle
                      endif

                      npts       = npts + 1
                      ilon(npts) = i
                      jlat(npts) = jstart
                      npts       = npts + 1
                      ilon(npts) = i
                      jlat(npts) = jend                         
                   enddo

                   !--- right and left boundary -----------------------------------------------
                   istart = is - step
                   iend   = is + step
                   if(src_modulo) then
                      if( istart < 1)       istart = istart + nlon_in
                      if( iend   > nlon_in) iend   = iend   - nlon_in
                   else 
                      istart = max(istart,1)
                      iend   = min(iend, nlon_in)
                   endif
                   do l = -step, step
                      j = js+l
                         if( j < 1 .or. j > nlat_in) cycle
                         npts = npts+1
                         ilon(npts) = istart
                         jlat(npts) = j
                         npts = npts+1
                         ilon(npts) = iend
                         jlat(npts) = j
                  end do
                end if

                !--- find the surrouding points             
                do l = 1, npts
                   i = ilon(l)
                   j = jlat(l)
                   ip1 = i+1
                   if(ip1>nlon_in) then
                      if(src_modulo) then
                         ip1 = 1
                      else
                         cycle
                      endif
                   endif
                   jp1 = j+1
                   if(jp1>nlat_in) cycle
                   lon1 = lon_in(i,  j);   lat1 = lat_in(i,j)
                   lon2 = lon_in(ip1,j);   lat2 = lat_in(ip1,j)
                   lon3 = lon_in(ip1,jp1); lat3 = lat_in(ip1,jp1)
                   lon4 = lon_in(i,  jp1); lat4 = lat_in(i,  jp1)  

                   if(lon .lt. lon_min .or. lon .gt. lon_max) then
                      if(i .ne. nlon_in) then
                         cycle
                      else
                         if(lon .lt. lon_min) then
                             lon1 = lon1 -tpi; lon4 = lon4 - tpi
                         else if(lon .gt. lon_max) then
                             lon2 = lon2 +tpi; lon3 = lon3 + tpi
                         endif
                      endif
                   endif

                   if(lat .ge. intersect(lon1,lat1,lon2,lat2,lon))then ! south
                      if(lon .le. intersect(lat2,lon2,lat3,lon3,lat))then ! east
                         if(lat .le. intersect(lon3,lat3,lon4,lat4,lon))then !north
                            if(lon .ge. intersect(lat4,lon4,lat1,lon1,lat))then ! west
                               found = .true.
                               is=i; js=j  
                               Interp % i_lon (m,n,1) = i; Interp % i_lon (m,n,2) = ip1
                               Interp % j_lat (m,n,1) = j; Interp % j_lat (m,n,2) = jp1
                               exit
                            endif
                         endif
                      endif
                   endif
                enddo
                step = step + 1
             enddo
          endif
          if(.not.found) then
             call mpp_error(FATAL, &
                  'horiz_interp_bilinear_mod: the destination point is not inside the source grid' )
          endif
       enddo
    enddo

  end subroutine find_neighbor

  !#######################################################################
  ! 
  ! The function will return true if the point x,y is inside a polygon, or
  ! NO if it is not.  If the point is exactly on the edge of a polygon,
  ! the function will return .true.
  !
  ! real polyx(:) : longitude coordinates of corners 
  ! real polyx(:) : latitude  coordinates of corners
  ! real x,y      : point to be tested 
  ! ??? How to deal with truncation error.
  !
  function inside_polygon(polyx, polyy, x, y)
     real, dimension(:), intent(in) :: polyx, polyy
     real,               intent(in) :: x, y
     logical                        :: inside_polygon
     integer                        :: i, j, nedges
     real                           :: xx

     inside_polygon = .false.
     nedges = size(polyx(:))
     j = nedges
     do i = 1, nedges
        if( (polyy(i) < y .AND. polyy(j) >= y) .OR. (polyy(j) < y .AND. polyy(i) >= y) ) then
           xx = polyx(i)+(y-polyy(i))/(polyy(j)-polyy(i))*(polyx(j)-polyx(i))
           if( xx == x ) then
             inside_polygon = .true.
             return
           else if( xx < x ) then
             inside_polygon = .not. inside_polygon
           endif
        endif
        j = i
     enddo

     return

  end function inside_polygon

  !#######################################################################
  ! this routine will search the source grid to fine the grid box that encloses 
  ! each destination grid.
  subroutine find_neighbor_new( Interp, lon_in, lat_in, lon_out, lat_out, src_modulo, no_crash )
    type(horiz_interp_type), intent(inout) :: Interp
    real, intent(in),       dimension(:,:) :: lon_in , lat_in
    real, intent(in),       dimension(:,:) :: lon_out, lat_out
    logical,                 intent(in)    :: src_modulo, no_crash
    integer                                :: nlon_in, nlat_in, nlon_out, nlat_out
    integer                                :: max_step, n, m, l, i, j, ip1, jp1, step
    integer                                :: is, js, jstart, jend, istart, iend, npts
    integer, allocatable, dimension(:)     :: ilon, jlat
    real                                   :: lon_min, lon_max, lon, lat, tpi
    logical                                :: found
    real                                   :: polyx(4), polyy(4)
    real                                   :: min_lon, min_lat, max_lon, max_lat    

    tpi = 2.0*pi
    nlon_in  = size(lon_in,1) ; nlat_in  = size(lat_in,2)
    nlon_out = size(lon_out,1); nlat_out = size(lon_out,2)

    lon_min = minval(lon_in);
    lon_max = maxval(lon_in);

    max_step = min(nlon_in,nlat_in)/2 ! can be adjusted if needed
    allocate(ilon(8*max_step), jlat(8*max_step) )

    do n = 1, nlat_out
       do m = 1, nlon_out
          found = .false.
          lon = lon_out(m,n)
          lat = lat_out(m,n)

          if(src_modulo) then
             if(lon .lt. lon_min) then
                lon = lon + tpi
             else if(lon .gt. lon_max) then
                lon = lon - tpi
             endif
          else
             if(lon .lt. lon_min .or. lon .gt. lon_max ) &
             call mpp_error(FATAL,'horiz_interp_bilinear_mod: ' //&
                  'when input grid is not modulo, output grid should locate inside input grid')
          endif
          !--- search for the surrounding four points locatioon.
          if(m==1 .and. n==1) then
             J_LOOP: do j = 1, nlat_in-1
                do i = 1, nlon_in
                   ip1 = i+1
                   jp1 = j+1
                   if(i==nlon_in) then
                      if(src_modulo)then
                         ip1 = 1
                      else
                         cycle
                      endif
                   endif
                   
                   polyx(1) = lon_in(i,  j);   polyy(1) = lat_in(i,j)
                   polyx(2) = lon_in(ip1,j);   polyy(2) = lat_in(ip1,j)
                   polyx(3) = lon_in(ip1,jp1); polyy(3) = lat_in(ip1,jp1)
                   polyx(4) = lon_in(i,  jp1); polyy(4) = lat_in(i,  jp1) 
                   if(lon .lt. lon_min .or. lon .gt. lon_max) then
                      if(i .ne. nlon_in) then
                         cycle
                      else
                         if(lon .lt. lon_min) then
                             polyx(1) = polyx(1) -tpi; polyx(4) = polyx(4) - tpi
                         else if(lon .gt. lon_max) then
                             polyx(2) = polyx(2) +tpi; polyx(3) = polyx(3) + tpi
                         endif
                      endif
                   endif
 
                   min_lon = minval(polyx)
                   max_lon = maxval(polyx)
                   min_lat = minval(polyy)
                   max_lat = maxval(polyy)
!                   if( lon .GE. min_lon .AND. lon .LE. max_lon .AND. &
!                       lat .GE. min_lat .AND. lat .LE. max_lat ) then
!                      print*, 'i =', i, 'j = ', j
!                      print '(5f15.11)', lon, polyx
!                      print '(5f15.11)', lat, polyy    
!                   endif 

                   if(inside_polygon(polyx, polyy, lon, lat)) then
                      found = .true.
!                      print*, " found ", i, j
                      Interp % i_lon (m,n,1) = i; Interp % i_lon (m,n,2) = ip1
                      Interp % j_lat (m,n,1) = j; Interp % j_lat (m,n,2) = jp1
                      exit J_LOOP
                   endif
                enddo
             enddo J_LOOP
          else
             step = 0
             do while ( .not. found .and. step .lt. max_step )
                !--- take the adajcent point as the starting point
                if(m == 1) then
                   is = Interp % i_lon (m,n-1,1)
                   js = Interp % j_lat (m,n-1,1)
                else
                   is = Interp % i_lon (m-1,n,1)
                   js = Interp % j_lat (m-1,n,1)
                endif
                if(step==0) then
                   npts = 1
                   ilon(1) = is
                   jlat(1) = js
                else
                   npts = 0
                   !--- bottom and top boundary
                   jstart = max(js-step,1)
                   jend   = min(js+step,nlat_in)

                   do l = -step, step
                      i = is+l
                      if(src_modulo)then
                         if( i < 1) then
                            i = i + nlon_in
                         else if (i > nlon_in) then
                            i = i - nlon_in
                         endif
                         if( i < 1 .or. i > nlon_in) call mpp_error(FATAL, &
                              'horiz_interp_bilinear_mod: max_step is too big, decrease max_step' )
                      else
                         if( i < 1 .or. i > nlon_in) cycle
                      endif

                      npts       = npts + 1
                      ilon(npts) = i
                      jlat(npts) = jstart
                      npts       = npts + 1
                      ilon(npts) = i
                      jlat(npts) = jend                         
                   enddo

                   !--- right and left boundary -----------------------------------------------
                   istart = is - step
                   iend   = is + step
                   if(src_modulo) then
                      if( istart < 1)       istart = istart + nlon_in
                      if( iend   > nlon_in) iend   = iend   - nlon_in
                   else 
                      istart = max(istart,1)
                      iend   = min(iend, nlon_in)
                   endif
                   do l = -step, step
                      j = js+l
                         if( j < 1 .or. j > nlat_in) cycle
                         npts = npts+1
                         ilon(npts) = istart
                         jlat(npts) = j
                         npts = npts+1
                         ilon(npts) = iend
                         jlat(npts) = j
                  end do
                end if

                !--- find the surrouding points             
                do l = 1, npts
                   i = ilon(l)
                   j = jlat(l)
                   ip1 = i+1
                   if(ip1>nlon_in) then
                      if(src_modulo) then
                         ip1 = 1
                      else
                         cycle
                      endif
                   endif
                   jp1 = j+1
                   if(jp1>nlat_in) cycle
                   polyx(1) = lon_in(i,  j);   polyy(1) = lat_in(i,j)
                   polyx(2) = lon_in(ip1,j);   polyy(2) = lat_in(ip1,j)
                   polyx(3) = lon_in(ip1,jp1); polyy(3) = lat_in(ip1,jp1)
                   polyx(4) = lon_in(i,  jp1); polyy(4) = lat_in(i,  jp1)  
                   if(inside_polygon(polyx, polyy, lon, lat)) then
                      found = .true.
                      Interp % i_lon (m,n,1) = i; Interp % i_lon (m,n,2) = ip1
                      Interp % j_lat (m,n,1) = j; Interp % j_lat (m,n,2) = jp1
                      exit
                   endif
                enddo
                step = step + 1
             enddo
          endif
          if(.not.found) then
             if(no_crash) then
                Interp % i_lon (m,n,1:2) = DUMMY
                Interp % j_lat (m,n,1:2) = DUMMY
                print*,'lon,lat=',lon,lat ! snz
             else
                call mpp_error(FATAL, &
                    'horiz_interp_bilinear_mod: the destination point is not inside the source grid' )
             endif
          endif
       enddo
    enddo

  end subroutine find_neighbor_new

  !#######################################################################
  function intersect(x1, y1, x2, y2, x)
     real, intent(in) :: x1, y1, x2, y2, x
     real             :: intersect

     intersect = (y2-y1)*(x-x1)/(x2-x1) + y1

  return

  end function intersect

  !#######################################################################
  ! <SUBROUTINE NAME="horiz_interp_bilinear">

  !   <OVERVIEW>
  !      Subroutine for performing the horizontal interpolation between two grids.
  !   </OVERVIEW>
  !   <DESCRIPTION>
  !     Subroutine for performing the horizontal interpolation between two grids. 
  !     horiz_interp_bilinear_new must be called before calling this routine.
  !   </DESCRIPTION>
  !   <TEMPLATE>
  !     call horiz_interp_bilinear ( Interp, data_in, data_out, verbose, mask_in,mask_out, missing_value, missing_permit)
  !   </TEMPLATE>
  !   
  !   <IN NAME="Interp" TYPE="type(horiz_interp_type)">
  !     Derived-type variable containing interpolation indices and weights.
  !     Returned by a previous call to horiz_interp_bilinear_new.
  !   </IN>
  !   <IN NAME="data_in" TYPE="real, dimension(:,:)">
  !      Input data on source grid.
  !   </IN>
  !   <IN NAME="verbose" TYPE="integer, optional">
  !      flag for the amount of print output.
  !               verbose = 0, no output; = 1, min,max,means; = 2, still more
  !   </IN>
  !   <IN NAME="mask_in" TYPE="real, dimension(:,:),optional">
  !      Input mask, must be the same size as the input data. The real value of
  !      mask_in must be in the range (0.,1.). Set mask_in=0.0 for data points 
  !      that should not be used or have missing data. 
  !   </IN>
  !   <IN NAME="missing_value" TYPE="real, optional">
  !      Use the missing_value to indicate missing data.
  !   </IN>

  !   <IN NAME="missing_permit" TUPE="integer, optional">
  !      numbers of points allowed to miss for the bilinear interpolation. The value
  !      should be between 0 and 3.
  !   </IN>

  !   <OUT NAME="data_out" TYPE="real, dimension(:,:)">
  !      Output data on destination grid.
  !   </OUT>
  !   <OUT NAME="mask_out" TYPE="real, dimension(:,:),optional">
  !      Output mask that specifies whether data was computed.
  !   </OUT>

  subroutine horiz_interp_bilinear ( Interp, data_in, data_out, verbose, mask_in,mask_out, &
       missing_value, missing_permit, new_handle_missing )
    !-----------------------------------------------------------------------
    type (horiz_interp_type), intent(in)        :: Interp
    real, intent(in),  dimension(:,:)           :: data_in
    real, intent(out), dimension(:,:)           :: data_out
    integer, intent(in),               optional :: verbose
    real, intent(in), dimension(:,:),  optional :: mask_in
    real, intent(out), dimension(:,:), optional :: mask_out
    real, intent(in),                  optional :: missing_value
    integer, intent(in),               optional :: missing_permit
    logical, intent(in),               optional :: new_handle_missing
    !-----------------------------------------------------------------------
    integer :: nlon_in, nlat_in, nlon_out, nlat_out, n, m,         &
         is, ie, js, je, iverbose, max_missing, num_missing, &
         miss_in, miss_out, unit
    real    :: dwtsum, wtsum, min_in, max_in, avg_in, &
         min_out, max_out, avg_out, wtw, wte, wts, wtn
    real    :: mask(size(data_in,1), size(data_in,2) )
    logical :: set_to_missing, is_missing(4), new_handler
    real    :: f1, f2, f3, f4, middle, w, s

    num_missing = 0

    nlon_in  = Interp%nlon_src;  nlat_in  = Interp%nlat_src
    nlon_out = Interp%nlon_dst; nlat_out = Interp%nlat_dst

    if(present(mask_in)) then
       mask = mask_in
    else
       mask = 1.0
    endif

    if (present(verbose)) then
       iverbose = verbose
    else
       iverbose = 0
    endif

    if(present(missing_permit)) then
       max_missing = missing_permit
    else
       max_missing = 0
    endif

    if(present(new_handle_missing)) then
       new_handler = new_handle_missing
    else
       new_handler = .false.
    endif

    if(max_missing .gt. 3 .or. max_missing .lt. 0) call mpp_error(FATAL, &
         'horiz_interp_bilinear_mod: missing_permit should be between 0 and 3')

    if (size(data_in,1) /= nlon_in .or. size(data_in,2) /= nlat_in) &
         call mpp_error(FATAL,'horiz_interp_bilinear_mod: size of input array incorrect')

    if (size(data_out,1) /= nlon_out .or. size(data_out,2) /= nlat_out) &
         call mpp_error(FATAL,'horiz_interp_bilinear_mod: size of output array incorrect')

    if(new_handler) then
       if( .not. present(missing_value) )  call mpp_error(FATAL, &
            "horiz_interp_bilinear_mod: misisng_value must be present when new_handle_missing is .true.")
       if( present(mask_in) ) call mpp_error(FATAL, &
            "horiz_interp_bilinear_mod: mask_in should not be present when new_handle_missing is .true.")
       do n = 1, nlat_out
          do m = 1, nlon_out
             is = Interp % i_lon (m,n,1); ie = Interp % i_lon (m,n,2)
             js = Interp % j_lat (m,n,1); je = Interp % j_lat (m,n,2)
             wtw = Interp % wti   (m,n,1)
             wte = Interp % wti   (m,n,2)
             wts = Interp % wtj   (m,n,1)
             wtn = Interp % wtj   (m,n,2)

             is_missing = .false.
             num_missing = 0
             set_to_missing = .false.
             if(data_in(is,js) == missing_value) then
                num_missing = num_missing+1
                is_missing(1) = .true.
                if(wtw .GE. 0.5 .AND. wts .GE. 0.5) set_to_missing = .true.
             endif

             if(data_in(ie,js) == missing_value) then
                num_missing = num_missing+1
                is_missing(2) = .true.
                if(wte .GE. 0.5 .AND. wts .GE. 0.5) set_to_missing = .true.
             endif
             if(data_in(ie,je) == missing_value) then
                num_missing = num_missing+1
                is_missing(3) = .true.
                if(wte .GE. 0.5 .AND. wtn .GE. 0.5) set_to_missing = .true.
             endif
             if(data_in(is,je) == missing_value) then
                num_missing = num_missing+1
                is_missing(4) = .true.
                if(wtw .GE. 0.5 .AND. wtn .GE. 0.5) set_to_missing = .true.
             endif

             if( num_missing == 4 .OR. set_to_missing ) then
                data_out(m,n) = missing_value
                if(present(mask_out)) mask_out(m,n) = 0.0
                 cycle
             else if(num_missing == 0) then
                f1 = data_in(is,js)
                f2 = data_in(ie,js)
                f3 = data_in(ie,je)
                f4 = data_in(is,je)
                w = wtw
                s = wts
             else if(num_missing == 3) then  !--- three missing value 
                if(.not. is_missing(1) ) then
                   data_out(m,n) = data_in(is,js)
                else if(.not. is_missing(2) ) then
                   data_out(m,n) = data_in(ie,js)
                else if(.not. is_missing(3) ) then
                   data_out(m,n) = data_in(ie,je)
                else if(.not. is_missing(4) ) then
                   data_out(m,n) = data_in(is,je)
                endif
                if(present(mask_out) ) mask_out(m,n) = 1.0
                cycle
             else   !--- one or two missing value
                if( num_missing == 1) then
                   if( is_missing(1) .OR. is_missing(3) ) then
                      middle = 0.5*(data_in(ie,js)+data_in(is,je))
                   else
                      middle = 0.5*(data_in(is,js)+data_in(ie,je))
                   endif
                else ! num_missing = 2
                   if( is_missing(1) .AND. is_missing(2) ) then
                      middle = 0.5*(data_in(ie,je)+data_in(is,je))
                   else if( is_missing(1) .AND. is_missing(3) ) then
                      middle = 0.5*(data_in(ie,js)+data_in(is,je))
                   else if( is_missing(1) .AND. is_missing(4) ) then
                      middle = 0.5*(data_in(ie,js)+data_in(ie,je))
                   else if( is_missing(2) .AND. is_missing(3) ) then
                      middle = 0.5*(data_in(is,js)+data_in(is,je))
                   else if( is_missing(2) .AND. is_missing(4) ) then
                      middle = 0.5*(data_in(is,js)+data_in(ie,je))
                   else if( is_missing(3) .AND. is_missing(4) ) then
                      middle = 0.5*(data_in(is,js)+data_in(ie,js))
                   endif
                endif

                if( wtw .GE. 0.5 .AND. wts .GE. 0.5 ) then  ! zone 1
                   w = 2.0*(wtw-0.5)
                   s = 2.0*(wts-0.5)
                   f1 = data_in(is,js)
                   if(is_missing(2)) then
                      f2 = f1
                   else
                      f2 = 0.5*(data_in(is,js)+data_in(ie,js))
                   endif
                   f3 = middle
                   if(is_missing(4)) then
                      f4 = f1
                   else
                      f4 = 0.5*(data_in(is,js)+data_in(is,je))
                   endif
                else if( wte .GE. 0.5 .AND. wts .GE. 0.5 ) then  ! zone 2
                   w = 2.0*(1.0-wte)
                   s = 2.0*(wts-0.5)
                   f2 = data_in(ie,js)
                   if(is_missing(1)) then
                      f1 = f2
                   else
                      f1 = 0.5*(data_in(is,js)+data_in(ie,js))
                   endif
                   f4 = middle
                   if(is_missing(3)) then
                      f3 = f2
                   else
                      f3 = 0.5*(data_in(ie,js)+data_in(ie,je))
                   endif
                else if( wte .GE. 0.5 .AND. wtn .GE. 0.5 ) then  ! zone 3
                   w = 2.0*(1.0-wte)
                   s = 2.0*(1.0-wtn)
                   f3 = data_in(ie,je)
                   if(is_missing(2)) then
                      f2 = f3
                   else
                      f2 = 0.5*(data_in(ie,js)+data_in(ie,je))
                   endif
                   f1 = middle
                   if(is_missing(4)) then
                      f4 = f3
                   else
                      f4 = 0.5*(data_in(ie,je)+data_in(is,je))
                   endif
                else if( wtw .GE. 0.5 .AND. wtn .GE. 0.5 ) then  ! zone 4
                   w = 2.0*(wtw-0.5)
                   s = 2.0*(1.0-wtn)
                   f4 = data_in(is,je)
                   if(is_missing(1)) then
                      f1 = f4
                   else
                      f1 = 0.5*(data_in(is,js)+data_in(is,je))
                   endif
                   f2 = middle
                   if(is_missing(3)) then
                      f3 = f4
                   else
                      f3 = 0.5*(data_in(ie,je)+data_in(is,je))
                   endif
                else 
                   call mpp_error(FATAL, &
                      "horiz_interp_bilinear_mod: the point should be in one of the four zone")
                endif
             endif

             data_out(m,n) = f3 + (f4-f3)*w + (f2-f3)*s + ((f1-f2)+(f3-f4))*w*s 
            if(present(mask_out)) mask_out(m,n) = 1.0
          enddo
       enddo       
    else
       do n = 1, nlat_out
          do m = 1, nlon_out
             is = Interp % i_lon (m,n,1); ie = Interp % i_lon (m,n,2)
             js = Interp % j_lat (m,n,1); je = Interp % j_lat (m,n,2)
             wtw = Interp % wti   (m,n,1)
             wte = Interp % wti   (m,n,2)
             wts = Interp % wtj   (m,n,1)
             wtn = Interp % wtj   (m,n,2)

             if(present(missing_value) ) then
                num_missing = 0
                if(data_in(is,js) == missing_value) then
                   num_missing = num_missing+1
                   mask(is,js) = 0.0
                endif
                if(data_in(ie,js) == missing_value) then
                   num_missing = num_missing+1
                   mask(ie,js) = 0.0
                endif
                if(data_in(ie,je) == missing_value) then
                   num_missing = num_missing+1
                   mask(ie,je) = 0.0
                endif
                if(data_in(is,je) == missing_value) then
                   num_missing = num_missing+1
                   mask(is,je) = 0.0
                endif
             endif

             dwtsum = data_in(is,js)*mask(is,js)*wtw*wts &
                  + data_in(ie,js)*mask(ie,js)*wte*wts &
                  + data_in(ie,je)*mask(ie,je)*wte*wtn &
                  + data_in(is,je)*mask(is,je)*wtw*wtn 
             wtsum  = mask(is,js)*wtw*wts + mask(ie,js)*wte*wts  &
                  + mask(ie,je)*wte*wtn + mask(is,je)*wtw*wtn

             if(.not. present(mask_in) .and. .not. present(missing_value)) wtsum = 1.0

             if(num_missing .gt. max_missing ) then
                data_out(m,n) = missing_value
                if(present(mask_out)) mask_out(m,n) = 0.0
             else if(wtsum .lt. epsln) then 
                if(present(missing_value)) then
                   data_out(m,n) = missing_value
                else
                   data_out(m,n) = 0.0
                endif
                if(present(mask_out)) mask_out(m,n) = 0.0      
             else
                data_out(m,n) = dwtsum/wtsum
                if(present(mask_out)) mask_out(m,n) = wtsum
             endif
          enddo
       enddo
    endif
    !***********************************************************************
    ! compute statistics: minimum, maximum, and mean
    !-----------------------------------------------------------------------
    if (iverbose > 0) then

       ! compute statistics of input data

       call stats (data_in, min_in, max_in, avg_in, miss_in, missing_value, mask_in)

       ! compute statistics of output data
       call stats (data_out, min_out, max_out, avg_out, miss_out, missing_value, mask_out)

       !---- output statistics ----
       unit = stdout()
       write (unit,900)
       write (unit,901)  min_in ,max_in, avg_in
       if (present(mask_in))  write (unit,903)  miss_in
       write (unit,902)  min_out,max_out,avg_out
       if (present(mask_out)) write (unit,903)  miss_out

900    format (/,1x,10('-'),' output from horiz_interp ',10('-'))
901    format ('  input:  min=',f16.9,'  max=',f16.9,'  avg=',f22.15)
902    format (' output:  min=',f16.9,'  max=',f16.9,'  avg=',f22.15)
903    format ('          number of missing points = ',i6)

    endif

    return

  end subroutine horiz_interp_bilinear
  ! </SUBROUTINE>

  !#######################################################################
  ! <SUBROUTINE NAME="horiz_interp_bilinear_del">

  !   <OVERVIEW>
  !     Deallocates memory used by "horiz_interp_type" variables.
  !     Must be called before reinitializing with horiz_interp_bilinear_new.
  !   </OVERVIEW>
  !   <DESCRIPTION>
  !     Deallocates memory used by "horiz_interp_type" variables.
  !     Must be called before reinitializing with horiz_interp_bilinear_new.
  !   </DESCRIPTION>
  !   <TEMPLATE>
  !     call horiz_interp_bilinear_del ( Interp )
  !   </TEMPLATE>

  !   <INOUT NAME="Interp" TYPE="horiz_interp_type">
  !     A derived-type variable returned by previous call
  !     to horiz_interp_bilinear_new. The input variable must have
  !     allocated arrays. The returned variable will contain
  !     deallocated arrays.
  !   </INOUT>

  subroutine horiz_interp_bilinear_del( Interp )

    type (horiz_interp_type), intent(inout) :: Interp

    if(associated(Interp%wti))   deallocate(Interp%wti)
    if(associated(Interp%wtj))   deallocate(Interp%wtj)
    if(associated(Interp%i_lon)) deallocate(Interp%i_lon)
    if(associated(Interp%j_lat)) deallocate(Interp%j_lat)

  end subroutine horiz_interp_bilinear_del
  ! </SUBROUTINE>

  !#######################################################################

  function indp (value, array)
    integer                        :: indp
    real, dimension(:), intent(in) :: array
    real, intent(in)               :: value
    !
    !=======================================================================
    !
    !     indp = index of nearest data point within "array" corresponding to
    !            "value".

    !     inputs:
    !     value  = arbitrary data...same units as elements in "array"
    !     array  = array of data points  (must be monotonically increasing)

    !     output:
    !     indp =  index of nearest data point to "value"
    !             if "value" is outside the domain of "array" then indp = 1
    !             or "ia" depending on whether array(1) or array(ia) is
    !             closest to "value"
    !=======================================================================
    !
    integer i, ia, unit
    logical keep_going
    !
    ia = size(array(:))
    do i=2,ia
       if (array(i) .lt. array(i-1)) then
          unit = stdout()
          write (unit,*) &
               ' => Error: array must be monotonically increasing in "indp"' , &
               '           when searching for nearest element to value=',value
          write (unit,*) '           array(i) < array(i-1) for i=',i 
          write (unit,*) '           array(i) for i=1..ia follows:'
          call abort()
       endif
    enddo
    if (value .lt. array(1) .or. value .gt. array(ia)) then
       if (value .lt. array(1))  indp = 1
       if (value .gt. array(ia)) indp = ia
    else
       i=1
       keep_going = .true.
       do while (i .le. ia .and. keep_going)
          i = i+1
          if (value .le. array(i)) then
             indp = i
             if (array(i)-value .gt. value-array(i-1)) indp = i-1
             keep_going = .false.
          endif
       enddo
    endif
    return
  end function indp

  !######################################################################

end module horiz_interp_bilinear_mod
