program runoff_regrid
  !-----------------------------------------------------------------------
  !                   GNU General Public License                        
  !                                                                      
  ! This program 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; either version 2 of   
  ! the License, or (at your option) any later version.                 
  !                                                                      
  ! MOM 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.                                           
  !                                                                      
  ! For the full text of the GNU General Public License,                
  ! write to: Free Software Foundation, Inc.,                           
  !           675 Mass Ave, Cambridge, MA 02139, USA.                   
  ! or see:   http://www.gnu.org/licenses/gpl.html                      
  !-----------------------------------------------------------------------
  !
  ! <CONTACT EMAIL="GFDL.Climate.Model.Info@noaa.gov"> Sergey Malyshev </CONTACT>
  ! <CONTACT EMAIL="GFDL.Climate.Model.Info@noaa.gov"> Zhi Liang </CONTACT>
  !
  ! <OVERVIEW>
  !   This program works together with program make_xgrids, create_grid to 
  !   remap runoff data from a spherical grid onto any grid (spherical or tripolar) 
  !   using conservative scheme.
  ! </OVERVIEW>

  !<DESCRIPTION>
  ! This program expects to read runoff data from a netcdf file, which is specfied
  ! by the namelist variable "src_data". The name of the runoff field in input file is specified 
  ! by the namelist variable "src_fld_name". The output file is a netcdf file specified
  ! by the namelist variable "dst_data". The name of the runoff field in output file is specified 
  ! by the namelist variable "src_fld_name". Both the source grid and  destination grid 
  ! are read from netcdf file "grid_spec.nc". The grid_spec.nc is expected contains
  ! field 'x_T' and 'y_T', which is the new name convention used as a prototype for ESMF.
  ! The source grid read from variable "xta" and "yta" of "grid_spec.nc" should match 
  ! the grid in "src_data". "grid_spec.nc" is generated by 
  ! " make_xgrid -a atmos_grid.nc -l atmos_grid.nc -o dst_grid.nc ",
  ! where dst_grid is the grid which you want remap runoff data on and atmos_grid.nc 
  ! is generated through program "create_grid", which is located at the same directory
  ! as this program. make_xgrid is a c program that generates exchange grid between 
  ! atmos, ocean and land grid. make_xgrid source code is located at 
  ! "../generate_grids/make_xgrids". If there is any land point has runoff data 
  ! after remapping runoff data onto destination grid, the runoff value of that land point
  ! will be moved to the nearest ocean point.
  !
  !</DESCRIPTION>
  !

  implicit none

#include "netcdf.inc"

  !--- namelist interface ----------------------------------------------
  !<NAMELIST NAME="runoff_regrid_nml">
  ! <DATA NAME="src_data" TYPE="character(len=128)" >
  !  Name of input file containing runoff data.
  ! </DATA>
  ! <DATA NAME="src_fld_name" TYPE="character(len=128)" >
  !  Name of runoff field in input file.
  ! </DATA>
  ! <DATA NAME="dst_fld_name" TYPE="character(len=128)" >
  !  Name of runoff field in output file.
  ! </DATA>
  ! <DATA NAME="dst_data" TYPE="character(len=128)" >
  !  Name of output file containing runoff data.
  ! </DATA>
  !</NAMELIST>

  character(len=128) :: src_fld_name = 'runoff'      ! name of the input runoff field
  character(len=128) :: dst_fld_name = 'runoff'      ! name of the output runoff field
  character(len=128) :: src_data     = 'src_data.nc' ! name of the runoff source data file
  character(len=128) :: dst_data     = 'dst_data.nc' ! name of the runoff output file

  namelist /runoff_regrid_nml/ src_fld_name, dst_fld_name, src_data, dst_data


  !--- version information ---------------------------------------------
  character(len=128)   :: version= '$ID$'
  character(len=128)   :: tagname= '$Name: tikal $'
  !--- other variables -------------------------------------------------
  real,    allocatable :: time_value(:)
  real,    allocatable :: runoff(:,:,:), runoff_src(:,:,:)
  logical, allocatable :: mask(:,:)
  real,    allocatable :: lon(:), lat(:), x_t_dst(:,:), y_t_dst(:,:), grid_x_t(:), grid_y_t(:)
  integer, allocatable :: i_atm_atmxocn(:), j_atm_atmxocn(:), i_ocn_atmxocn(:), j_ocn_atmxocn(:)
  integer, allocatable :: i_atm_atmxlnd(:), j_atm_atmxlnd(:), i_lnd_atmxlnd(:), j_lnd_atmxlnd(:)
  real,    allocatable :: area_atmxocn(:), area_atmxlnd(:)
  integer              :: ni_src, nj_src, ni_dst, nj_dst, ntime
  real,      parameter :: PI = 3.141592658979
  integer              :: stdout = 6 

  !--- initialize routine 
  call runoff_regrid_init()

  !--- read source data file -------------------------------------------
  call read_src_data()

  !--- read exchange grid information ----------------------------------
  call read_xgrid()

  !--- read the runoff data and remap onto destination grid.
  call remap_data()

  !--- write runoff data to destination file
  call write_dst_data()

  !--- release memory
  call runoff_regrid_end


contains

  !#######################################################################
  !--- initialization routine: read namelist, write namelist and version information
  subroutine runoff_regrid_init
    integer :: unit=7, iostat 
    logical :: opened

    !--- read namelist 
    do 
       inquire( unit, opened = opened )
       if( .not. opened ) exit
       unit = unit + 1
       if( unit.EQ.100 )call error_handler( 'Error: unable to locate unit number.' )
    enddo

    open (unit, file='input.nml', iostat = iostat )
    if(iostat .ne. 0) call error_handler('error in open input.nml file')

    read (unit, runoff_regrid_nml, iostat=iostat)
    if(iostat .ne. 0) call error_handler('error in while reading runoff_regrid_nml ')

    close(unit)

    !--- write out namelist and version information
    write( stdout,'(/a)' )'runoff_regrid '//trim(version)//trim(tagname)
    write( stdout, runoff_regrid_nml)

  end subroutine runoff_regrid_init

  !#######################################################################
  ! read runoff data from src_data
  subroutine read_src_data

  integer :: dims(4), start(4), nread(4)
  integer :: rcode, ncid, ndims, nvar, id_time, id_runoff, i
  real    :: missing_value

  rcode = nf_open(trim(src_data), NF_NOWRITE, ncid)

  rcode = nf_inq_varid(ncid, trim(src_fld_name), id_runoff )

  rcode = nf_inq_varndims(ncid, id_runoff, ndims)

  rcode = nf_inq_vardimid(ncid, id_runoff, dims)

  rcode = nf_inq_dimlen(ncid, dims(1), ni_src)
  rcode = nf_inq_dimlen(ncid, dims(2), nj_src)

  ntime = 1
  if(ndims .gt. 2) then
     rcode = nf_inq_dimlen(ncid, dims(3), ntime)
  endif

  !--- get the time values if there are more than one time level
  if(ndims .gt. 2) then
    allocate(time_value(ntime) )
    rcode = nf_inq_varid(ncid, 'Time', id_time)
    rcode = nf_get_var_double(ncid, id_time, time_value)
  endif

  allocate(runoff_src(ni_src, nj_src, ntime) )

  !--- read the source runoff data
  start = 1; nread = 1
  nread(1) = ni_src; nread(2) = nj_src; nread(3) = ntime
  rcode = nf_get_vara_double(ncid, id_runoff, start, nread, runoff_src)

  !--- get the missing value -------------------------------------------
  missing_value = -1.0e20
  rcode = nf_get_att_double(ncid, id_runoff, 'missing_value', missing_value)
  if(rcode == 0) then
    where(runoff_src == missing_value) runoff_src = 0.0
  endif

  rcode = nf_close(ncid)

  end subroutine read_src_data

  !#######################################################################
  ! read component grid and exchange grid from grid_spec.nc
  subroutine read_xgrid
    integer              :: rcode, ncid, ncells, ni_lon, nj_lat, i, j
    real,    allocatable :: wet(:,:)
    logical, allocatable :: mask0(:,:)

  !--- read xgrid from the grid_spec.nc file --------------------------
  rcode = nf_open('grid_spec.nc', NF_NOWRITE, ncid)
  call error_handler('error in open file grid_spec.nc', rcode)

  ni_dst = get_dimlen(ncid, 'grid_x_T')
  nj_dst = get_dimlen(ncid, 'grid_y_T')

  allocate(grid_x_t(ni_dst), grid_y_t(nj_dst) )
  allocate(wet(ni_dst,nj_dst), x_t_dst(ni_dst,nj_dst), y_t_dst(ni_dst,nj_dst) )
  allocate(mask0(ni_dst,nj_dst),mask(ni_dst,nj_dst)) 

  !--- get grid axis data ----------------------------------------------
  call get_var_real_1d(ncid,'grid_x_T', grid_x_t)
  call get_var_real_1d(ncid,'grid_y_T', grid_y_t)

  call get_var_real_2d(ncid, 'wet', wet)
  
  ! calculate coastal mask, on ocn grid
  mask0 = (wet > 0.5)
  do j = 1, nj_dst
     do i = 1, ni_dst
        mask(i,j) = coastal(mask0,i,j)
     enddo
  enddo

  deallocate(mask0)

  call get_var_real_2d(ncid, 'x_T',  x_t_dst)
  call get_var_real_2d(ncid, 'y_T',  y_t_dst)

  !--- get exchange grid information
  ncells = get_dimlen(ncid, 'I_ATM_ATMxOCN')
  allocate(i_atm_atmxocn(ncells), j_atm_atmxocn(ncells), &
           i_ocn_atmxocn(ncells), j_ocn_atmxocn(ncells), &
           area_atmxocn(ncells) )

  call get_var_int_1d(ncid, 'I_ATM_ATMxOCN', i_atm_atmxocn)
  call get_var_int_1d(ncid, 'J_ATM_ATMxOCN',j_atm_atmxocn )
  call get_var_int_1d(ncid, 'I_OCN_ATMxOCN',i_ocn_atmxocn )
  call get_var_int_1d(ncid, 'J_OCN_ATMxOCN',j_ocn_atmxocn )
  call get_var_real_1d(ncid, 'AREA_ATMxOCN',area_atmxocn )

  ncells = get_dimlen(ncid, 'I_ATM_ATMxLND')
  allocate(i_atm_atmxlnd(ncells), j_atm_atmxlnd(ncells), &
           i_lnd_atmxlnd(ncells), j_lnd_atmxlnd(ncells), &
           area_atmxlnd(ncells) )

  call get_var_int_1d(ncid, 'I_ATM_ATMxLND', i_atm_atmxlnd)
  call get_var_int_1d(ncid, 'J_ATM_ATMxLND',j_atm_atmxlnd )
  call get_var_int_1d(ncid, 'I_LND_ATMxLND',i_lnd_atmxlnd )
  call get_var_int_1d(ncid, 'J_LND_ATMxLND',j_lnd_atmxlnd )
  call get_var_real_1d(ncid, 'AREA_ATMxLND',area_atmxlnd )

  !--- get the longitude and latitude of simple lon lat grid
  !--- it should has the same size as the runoff data ni_src, nj_src
  ni_lon = get_dimlen(ncid, 'xta')
  nj_lat = get_dimlen(ncid, 'yta')

  if(ni_lon .ne. ni_src .or. nj_lat .ne. nj_src)  &
      call error_handler('Error: size mismatch between file grid_spec.nc and file '//trim(src_data) )

  allocate(lon(ni_src), lat(nj_src))

  call get_var_real_1d(ncid, 'xta', lon)
  call get_var_real_1d(ncid, 'yta', lat)  

  rcode = nf_close(ncid)

  end subroutine read_xgrid

  !#######################################################################
  ! remapping runoff data onto destination grid
  subroutine remap_data
     integer              :: i, j, k, ncells_axo, ncells_axl
     real,    allocatable :: area_ocn(:,:), area_lnd(:,:), runoff_lnd(:,:), x(:), xl(:)
     integer, allocatable :: imap(:,:), jmap(:,:)
     real                 :: D2R

     D2R = PI/180.
     ncells_axo = size(i_atm_atmxocn(:))
     ncells_axl = size(i_atm_atmxlnd(:))

     allocate(area_ocn(ni_dst, nj_dst), area_lnd(ni_src, nj_src))
     allocate(runoff_lnd(ni_src, nj_src), runoff(ni_dst, nj_dst, ntime) )
     allocate(x(ncells_axo), xl(ncells_axl) )
     allocate(imap(ni_src, nj_src), jmap(ni_src, nj_src) )

     imap = 0
     jmap = 0
     runoff = 0
     do k = 1, ntime

        ! put input runoff to ocean grid
        do i = 1, ncells_axo
           x(i) = runoff_src(i_atm_atmxocn(i),j_atm_atmxocn(i),k)
        enddo

        area_ocn = 0
        do i = 1, ncells_axo
           runoff(i_ocn_atmxocn(i),j_ocn_atmxocn(i), k) = &
                runoff(i_ocn_atmxocn(i),j_ocn_atmxocn(i), k)+area_atmxocn(i)*x(i)
           area_ocn(i_ocn_atmxocn(i),j_ocn_atmxocn(i)) = &
                area_ocn(i_ocn_atmxocn(i),j_ocn_atmxocn(i))+area_atmxocn(i)
        enddo

        where (area_ocn(:,:) > 0) runoff(:,:,k) = runoff(:,:,k)/area_ocn(:,:)

        ! put input runoff to land grid -- it will be mostly zero, but in places where
        ! the input data grid and model ocean grid do not coincide
        do i = 1, ncells_axl
           xl(i) = runoff_src(i_atm_atmxlnd(i),j_atm_atmxlnd(i),k)
        enddo
        runoff_lnd = 0
        area_lnd = 0
        do i = 1, ncells_axl
           runoff_lnd(i_lnd_atmxlnd(i),j_lnd_atmxlnd(i)) = &
                runoff_lnd(i_lnd_atmxlnd(i),j_lnd_atmxlnd(i))+area_atmxlnd(i)*xl(i)
           area_lnd(i_lnd_atmxlnd(i),j_lnd_atmxlnd(i)) = &
                area_lnd(i_lnd_atmxlnd(i),j_lnd_atmxlnd(i))+area_atmxlnd(i)
        enddo
        where (area_lnd > 0) runoff_lnd = runoff_lnd/area_lnd
        write (stdout,*) 'At time level ', k,' There are ',count(runoff_lnd/=0), 'land points that has runoff value'

        ! calculate remapping
        do j = 1,nj_src
        do i = 1,ni_src
           if(runoff_lnd(i,j)>0) then
              if(imap(i,j)==0) &
                   call nearest(mask,x_t_dst*D2R,y_t_dst*D2R,lon(i)*D2R,lat(j)*D2R,imap(i,j),jmap(i,j))
              runoff(imap(i,j),jmap(i,j),k) = runoff(imap(i,j),jmap(i,j),k) + &
                   runoff_lnd(i,j)*area_lnd(i,j)/area_ocn(imap(i,j),jmap(i,j))
           endif
        enddo
        enddo

     enddo


  end subroutine remap_data

  !#####################################################################
  ! write out output
  subroutine write_dst_data
    integer           :: id_lon, id_lat, id_time, id_time_src, id_runoff_src, id_runoff
    integer           :: type, ndims, natts, rcode, ncid, ncid_src, i, natts_time
    integer           :: dims(4), start(4), nwrite(4)
    character(len=64) :: name
    integer           ::fsize = 65536, inital = 0


    rcode = nf_open(trim(src_data), NF_NOWRITE, ncid_src)
    call error_handler('error in opening file '//trim(src_data) , rcode)  
    rcode = nf_inq_varid(ncid_src, trim(src_fld_name), id_runoff_src )
    call error_handler('error in inquring varible id of variable'//trim(src_fld_name) , rcode) 
    rcode = nf_inq_varndims(ncid_src, id_runoff_src, ndims ) 
    call error_handler('error in inquring ndims of variable'//trim(src_fld_name) , rcode) 
    rcode = nf_inq_varnatts(ncid_src, id_runoff_src, natts )
    call error_handler('error in inquring natts of variable'//trim(src_fld_name) , rcode) 


#ifdef use_netCDF3
    rcode = NF__CREATE(dst_data, NF_CLOBBER, inital, fsize, ncid)
#else  
    rcode = NF__CREATE(dst_data, IOR(NF_NETCDF4,NF_CLASSIC_MODEL), inital, fsize, ncid )
#endif
    call error_handler('error in creating file '//trim(dst_data), rcode)    

    rcode = nf_def_dim(ncid, 'grid_x_T', ni_dst, dims(1))
    call error_handler('error in defining x dimension' , rcode)    
    rcode = nf_def_dim(ncid, 'grid_y_T', nj_dst, dims(2))
    call error_handler('error in defining y dimension' , rcode) 
    if(ndims > 2) then
       rcode = nf_def_dim(ncid, 'Time', NF_UNLIMITED, dims(3))
       call error_handler('error in defining time dimension' , rcode) 
    endif

    rcode = nf_def_var(ncid, 'grid_x_T', NF_DOUBLE, 1, dims(1), id_lon)
    call error_handler('error in defining variable grid_x_T' , rcode)  
    rcode = nf_put_att_text(ncid,id_lon, 'long_name', 34,'Nominal Longitude of T-cell center')
    call error_handler('error in putting longnmae of variable grid_x_T' , rcode)  
    rcode = nf_put_att_text(ncid,id_lon, 'cartesian_axis', 1,'X')
    call error_handler('error in putting  cartesian_axis of variable grid_x_T' , rcode)  
    rcode = nf_put_att_text(ncid,id_lon, 'units', 11,'degree_east')
    call error_handler('error in putting units of variable grid_x_T' , rcode)  

    rcode = nf_def_var(ncid, 'grid_y_T', NF_DOUBLE, 1, dims(2), id_lat)
    call error_handler('error in defining variable grid_y_T' , rcode)  
    rcode = nf_put_att_text(ncid,id_lat, 'long_name', 33,'Nominal latitude of T-cell center')
    call error_handler('error in putting longnmae of variable grid_x_T' , rcode)  
    rcode = nf_put_att_text(ncid,id_lat, 'cartesian_axis', 1,'Y')
    call error_handler('error in putting  cartesian_axis of variable grid_x_T' , rcode)  
    rcode = nf_put_att_text(ncid,id_lat, 'units', 12,'degree_north')
    call error_handler('error in putting units of variable grid_x_T' , rcode)  

    if(ndims > 2) then
       rcode = nf_def_var(ncid, 'Time', NF_DOUBLE, 1, dims(3), id_time)
       call error_handler('error in defining variable Time' , rcode)  
       rcode = nf_inq_varid(ncid_src,'Time',id_time_src)
       call error_handler('error in inquiring variable id of Time of file'//trim(src_data) , rcode)  
       rcode = nf_inq_varnatts(ncid_src, id_time_src, natts_time )
       call error_handler('error in inquring natts of variable Time', rcode) 
       do i = 1, natts_time
          rcode = nf_inq_attname(ncid_src, id_time_src, i, name)
          call error_handler('error in inquring attribute '//trim(name), rcode)  
          rcode = nf_copy_att(ncid_src,id_time_src, trim(name), ncid, id_time)
          call error_handler('error in copy attribute '//trim(name)//' of Time', rcode)  
       enddo
    endif

    if(ndims .le. 2) then
       rcode = nf_def_var(ncid, trim(dst_fld_name), NF_DOUBLE, 2, dims(1:2), id_runoff)
    else
       rcode = nf_def_var(ncid, trim(dst_fld_name), NF_DOUBLE, 3, dims(1:3), id_runoff)
    endif
    call error_handler('error in defining variable '//trim(dst_fld_name) , rcode)  

    !--- write the field attribute ---------------------------------------
    do i = 1, natts
       rcode = nf_inq_attname(ncid_src, id_runoff_src, i, name)
       call error_handler('error in inquring attribute '//trim(name), rcode)  
       if(trim(name) .ne. '_FillValue') then
          rcode = nf_copy_att(ncid_src,id_runoff_src, trim(name), ncid, id_runoff)
          call error_handler('error in copy attribute '//trim(name)//' of variable '//trim(src_fld_name), rcode)  
       endif
    enddo

    rcode = nf_enddef(ncid)

    !--- write out data --------------------------------------------------
    start = 1; nwrite = 1;
    nwrite(1) = ni_dst
    rcode = nf_put_vara_double(ncid, id_lon, start, nwrite, grid_x_t)
    call error_handler('error in putting grid_x_T data', rcode)
    nwrite(1) = nj_dst
    rcode = nf_put_vara_double(ncid, id_lat, start, nwrite, grid_y_t)
    call error_handler('error in putting grid_y_T data', rcode)

    do i =1, ntime
       if(ndims>2) then
          rcode = nf_put_var1_double(ncid, id_time, i, time_value(i))
          call error_handler('error in putting time data', rcode)
       endif
       nwrite(1) = ni_dst; nwrite(2) = nj_dst
       start(3)  = i
       rcode = nf_put_vara_double(ncid, id_runoff, start, nwrite, runoff(:,:,i))  
       call error_handler('error in putting data of '//trim(dst_fld_name), rcode)
    enddo

    rcode = nf_close(ncid_src)
    rcode = nf_close(ncid)

  end subroutine write_dst_data

  !#####################################################################
  ! get the dimension length of any one dimensional variable
  function get_dimlen(ncid, name)
    integer,          intent(in) :: ncid
    character(len=*), intent(in) :: name
    integer                      :: get_dimlen
    integer                      :: varid, rcode, dims(1)

    rcode = nf_inq_varid(ncid, trim(name), varid)
    call error_handler('error in inquiring variable id of '//trim(name), rcode)

    rcode = nf_inq_vardimid(ncid, varid, dims)
    call error_handler('error in inquiring dimension id of '//trim(name), rcode)

    rcode = nf_inq_dimlen(ncid, dims(1), get_dimlen)
    call error_handler('error in inquiring dimension length of '//trim(name), rcode)

  end function get_dimlen

  !#####################################################################
  ! read the 1d real data from netcdf file.
  subroutine get_var_real_1d(ncid, name, data)
    integer,             intent(in) :: ncid
    character(len=*),    intent(in) :: name
    real, dimension(:), intent(out) :: data
    integer                         :: rcode, varid

    rcode = nf_inq_varid(ncid, name, varid)
    call error_handler('error in inquiring variable id of '//trim(name), rcode)

    rcode = nf_get_var_double(ncid, varid, data)
    call error_handler('error in reading data of '//trim(name), rcode)


  end subroutine get_var_real_1d

  !#####################################################################
  ! read the 2d real data from netcdf file.
  subroutine get_var_real_2d(ncid, name, data)
    integer,               intent(in) :: ncid
    character(len=*),      intent(in) :: name
    real, dimension(:,:), intent(out) :: data
    integer                           :: rcode, varid


    rcode = nf_inq_varid(ncid, name, varid)
    call error_handler('error in inquiring variable id of '//trim(name), rcode)

    rcode = nf_get_var_double(ncid, varid, data)
    call error_handler('error in reading data of '//trim(name), rcode)

  end subroutine get_var_real_2d

  !#####################################################################
  ! read the 1d integer data from netcdf file.
  subroutine get_var_int_1d(ncid, name, data)
    integer,                intent(in) :: ncid
    character(len=*),       intent(in) :: name
    integer, dimension(:), intent(out) :: data
    integer                            :: rcode, varid

    rcode = nf_inq_varid(ncid, name, varid)
    call error_handler('error in inquiring variable id of '//trim(name), rcode)
    rcode = nf_get_var_int(ncid, varid, data)
    call error_handler('error in reading data of '//trim(name), rcode)

  end subroutine get_var_int_1d

  !#####################################################################
  ! find the nearest ocean point of a given land point.
  subroutine nearest(mask, lon, lat, plon, plat, iout, jout)
    ! finde nearest point that is not masked out; brute force approach 
    logical, intent(in) :: mask(:,:)  ! mask of valid input points
    real,    intent(in) :: lon(:,:)   ! longitudes of input grd central points
    real,    intent(in) :: lat(:,:)   ! latitudes of input grd central points
    real,    intent(in) :: plon, plat ! coordinates of destination point
    integer, intent(out):: iout, jout ! output map values

    integer :: i,j
    real    :: r, r1

    r = distance(0.0,PI/2,0.0,-PI/2)  ! some big value

    do j = 1, size(mask,2)
       do i = 1, size(mask,1)
          if (.not.mask(i,j)) cycle
          r1 = distance(plon,plat,lon(i,j),lat(i,j))
          if ( r1 < r ) then
             iout = i
             jout = j
             r = r1
          endif
       enddo
    enddo

  end subroutine nearest

  !#####################################################################
  ! returns true if (i,j) is a coastal point of ocean
  function coastal(wet, i,j)
    logical, intent(in) :: wet(:,:)
    integer, intent(in) :: i,j
    logical             :: coastal

    type delta
       integer :: i,j
    end type delta

    integer :: i0,j0,k
    type(delta):: stencil(4)
    data stencil/delta(-1,0),delta(1,0),delta(0,-1),delta(0,1)/

    coastal = .false.
    if(.not.wet(i,j)) return

    do k = 1, size(stencil(:))
       i0 = lbound(wet,1)+modulo(i+stencil(k)%i-lbound(wet,1),size(wet,1)); 
       j0 = min( max(j+stencil(k)%j,lbound(wet,2)), ubound(wet,2) )

       coastal = coastal.or..not.wet(i0,j0)
    enddo
  end function coastal

!#####################################################################
  ! calculate a distance in 3-d between points on unit square
  function distance(lon1, lat1, lon2, lat2)
    real, intent(in) :: lon1,lat1,lon2,lat2 
    real             :: distance
    real             :: x1,y1,z1, x2,y2,z2

    z1 = sin(lat1)           ;  z2 = sin(lat2)
    x1 = cos(lat1)*cos(lon1) ;  x2 = cos(lat2)*cos(lon2)
    y1 = cos(lat1)*sin(lon1) ;  y2 = cos(lat2)*sin(lon2)

    distance = sqrt((z1-z2)**2+(x1-x2)**2+(y1-y2)**2)
  end function distance

  !#####################################################################
  ! error handling routine.
  subroutine error_handler(mesg, status)
    character(len=*),  intent(in) :: mesg
    integer, optional, intent(in) :: status
    character(len=256) :: msg

    if(present(status)) then
       if(status == 0) return
       msg = nf_strerror(status)
       msg = trim(mesg)//': '// trim(msg)
    else
       msg = trim(mesg)
    endif

    write(stdout,*) msg

    call ABORT()

  end subroutine error_handler

  !#######################################################################
  ! destructor: release memory.
  subroutine runoff_regrid_end

    deallocate(runoff, runoff_src, mask, lon, lat, x_t_dst, y_t_dst)
    deallocate(i_atm_atmxocn, j_atm_atmxocn, i_ocn_atmxocn, j_ocn_atmxocn)
    deallocate(i_atm_atmxlnd, j_atm_atmxlnd, i_lnd_atmxlnd, j_lnd_atmxlnd)
    deallocate(area_atmxocn, area_atmxlnd)

  end subroutine runoff_regrid_end

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

end program runoff_regrid