!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                   !!
!!                   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                      !!
!!                                                                   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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= "Sergey.Malyshev@noaa.gov"> Sergey Malyshev </CONTACT>
  ! <CONTACT EMAIL= "Zhi.Liang@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: mom4p1_pubrel_dec2009_nnz $'
  !--- 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

    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) 


    rcode = nf_create(trim(dst_data),NF_WRITE, ncid)
    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