!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                   !!
!!                   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 regrid_3d
  !-----------------------------------------------------------------------
  !                   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= "Bonnie.Samuels@noaa.gov">Bonnie Samuels </CONTACT>
  ! <CONTACT EMAIL="Zhi.Liang@noaa.gov">Zhi Liang</CONTACT>
  ! <REVIEWER EMAIL="Matthew.Harrison@noaa.gov">M.J. Harrison</REVIEWER>
  !
  !<DESCRIPTION>
  ! regrid 3-d lat-lon gridded data to logically rectangular grid
  ! described by grid descriptor file. Applies only to scalar fields
  ! No missing points allowed on input grid.
  !
  !</DESCRIPTION>

  use mpp_mod,          only : mpp_error, mpp_pe,  mpp_npes, mpp_root_pe
  use mpp_mod,          only : FATAL, WARNING, stdout, stdlog, mpp_chksum
  use mpp_io_mod,       only : mpp_open, mpp_close, mpp_read, mpp_write, mpp_write_meta
  use mpp_io_mod,       only : mpp_copy_meta, axistype, fieldtype, atttype
  use mpp_io_mod,       only : mpp_get_atts, mpp_get_info, mpp_get_fields, mpp_get_times
  use mpp_io_mod,       only : mpp_get_axes, mpp_get_axis_data
  use mpp_io_mod,       only : MPP_RDONLY, MPP_NETCDF, MPP_MULTI, MPP_SINGLE, MPP_OVERWR
  use mpp_domains_mod,  only : mpp_update_domains, mpp_define_domains, mpp_global_field
  use mpp_domains_mod,  only : domain2d, mpp_define_layout, mpp_get_compute_domain
  use mpp_domains_mod,  only : mpp_domains_set_stack_size
  use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_end, horiz_interp_type
  use axis_utils_mod,   only : get_axis_cart, interp_1d
  use fms_mod,          only : fms_init, fms_end, open_namelist_file, close_file, file_exist
  use fms_mod,          only : check_nml_error, write_version_number, lowercase
  use constants_mod,    only : constants_init, PI

  implicit none

  integer, parameter :: max_fields       = 10
  integer, parameter :: max_ntimes_saved = 12

  !--- namelist interface
  !<NAMELIST NAME="regrid_3d_nml">
  ! <DATA NAME="src_file" TYPE="character(len=128)" DEFAULT="src_file.nc">
  !  Name of input file containing grid and data to be regridded.
  ! </DATA>
  ! <DATA NAME="numfields"  TYPE="integer" DEFAULT="2">
  !  Number of fields.
  ! </DATA>
  ! <DATA NAME="src_field_name" TYPE="character(len=128), dimension(max_fields)" >
  !  Name of input field(s). default is (/'temp', 'salt'/)
  ! </DATA>
  ! <DATA NAME="dest_field_name" TYPE="character(len=128), dimension(max_fields)" >
  !  Name of output field(s). If it is not specified in the namelist, it will 
  !  get the value from src_field_name
  ! </DATA>
  ! <DATA NAME="dest_grid" TYPE="character(len=128)" DEFAULT="dest_grid.nc">
  !  Name of grid descriptor file containing target grid information.
  ! </DATA>
  ! <DATA NAME="dest_file"  TYPE="character(len=128)" DEFAULT="dest_file.nc">
  !  Name of output file.
  ! </DATA>
  ! <DATA NAME="num_nbrs"  TYPE="integer" DEFAULT="10">
  ! Number of nearest neighbors for regridding  
  ! </DATA>
  ! <DATA NAME="max_dist"  TYPE="integer" DEFAULT="0.17" UNITS="radians">
  !  Maximum radial influence for regridding.
  ! </DATA>
  ! <DATA NAME="scale_factor" TYPE="real ">
  ! scaling factor for data (e.g. -1 to flip sign or 0.01 to convert from centimeters)
  ! </DATA>
  ! <DATA NAME="stop_crit"  TYPE="character(len=1),dimension(2)" DEFAULT="0.001">
  !  The stopping criteria when extrapping data onto missing points.
  ! </DATA>
  ! <DATA NAME="use_source_vertical_grid" TYPE="logical" DEFAULT=".false.">
  !  when use_source_vertical_grid is set to true, the destination data will 
  !  have the same vertical level as the source data. When use_source_vertical_grid 
  !  is false, the vertical grid of destination data will come from dest_grid. 
  !  A linear vertical interpolation will be done when the source vertical is different
  !  from destination vertical grid.
  ! </DATA>
  ! <DATA NAME="apply_mask"  TYPE="logical" DEFAULT="true">
  !  flag to indicate if the land/sea mask of source/destination grid will be applied 
  !  on the output dest_file. When apply_mask is false, the destination data will be 
  !  global data, i.e. no missing value in the destination data file. When apply_mask 
  !  is true, mask will be applied to the destination data. The mask can be either 
  !  source grid or destination grid determined by nml use_source_vertical_grid. 
  !  When use_source_vertical_grid is true, source grid mask will be applied, otherwise
  !  destination grid mask will be applied.
  ! </DATA>
  ! <DATA NAME="interp_method"  TYPE= "character(len=20)" >
  !  specifying the remapping method when remampping data onto current grid.
  !  Its value can be "spherical" or " bilinear". "spherical" interpolation is a 
  !  inverse distance weighted interpolation algorithm. Default value is "bilinear". 
  !  "bilinear" interpolation is recommanded, since bilinear interpolation will provide 
  !  more smooth results than "spherical" interpolation (especially when interpolating 
  !  from coarse grid to fine grid). Plus bilinear interpolation is much more efficiency 
  !  than "spherical interpolation". 
  ! </DATA>
  ! <DATA NAME="ntimes_saved" TYPE= "integer" >
  !  number of time levels to be saved. Its value has to be less than or equal to the number 
  !  of time levels in the source data file.
  ! </DATA>
  ! <DATA NAME="timelevel_saved" TYPE= "integer(max_ntimes_saved)" >
  !  specify the selection of time levels to be saved. The number of elements to be specified 
  !  should be equal to ntimes_saved. 
  ! </DATA>
  ! <DATA NAME="debug" TYPE="logical">
  ! For Debugging. Set true to print out chksum information for debugging reproducing ability 
  ! accross processors. default is false.
  ! </DATA>
  !</NAMELIST>
  character(len=128) :: src_file                          = 'src_file.nc'
  character(len=128) :: dest_grid                         = 'dest_grid.nc'
  character(len=128) :: dest_file                         = 'dest_file.nc'
  integer            :: numfields                         = 2
  character(len=128) :: src_field_name(max_fields)        = ''
  character(len=128) :: dest_field_name(max_fields)       = ''
  real               :: stop_crit(max_fields)             = 0.005 
  integer            :: num_nbrs                          = 5
  real               :: max_dist                          = 0.1
  real               :: scale_factor(max_fields)          = 1.0
  logical            :: use_source_vertical_grid          = .FALSE.
  logical            :: apply_mask                        = .TRUE.
  character(len=32)  :: interp_method                     = "bilinear"
  logical            :: debug                             = .FALSE.
  integer            :: ntimes_saved                      = 0
  integer            :: timelevel_saved(max_ntimes_saved) = 0

  namelist /regrid_3d_nml/ src_file, src_field_name, dest_field_name, numfields, dest_file, &
                           dest_grid, scale_factor, num_nbrs, max_dist, stop_crit,          &
                           interp_method, apply_mask, use_source_vertical_grid, debug,      &
                           ntimes_saved, timelevel_saved

  !---------------------------------------------------------------------
  integer            :: ni_src, nj_src, nk_src, ni_dst, nj_dst, nk_dst, ntime_src
  type(axistype)     :: depth_axis, time_axis, axes_dst(2)
  type(fieldtype)    :: field_lon_dst, field_lat_dst
  type(fieldtype)    :: dest_field(max_fields)
  integer            :: dest_unit         
  logical            :: time_axis_exists = .false.
  real, parameter    :: tol = 1.e-10   ! tolerance for detecting missing values
  real, parameter    :: max_val=1.e20
  real, parameter    :: rel_coef = 0.9
  integer, parameter :: max_iter = 2000
  logical            :: is_cyclic = .true.  ! we suppose the source data is always global data.
  real               :: missing(max_fields)
  real               :: D2R
  type(fieldtype), dimension(:), allocatable :: src_field
  integer                                    :: src_unit
  real, dimension(:),            allocatable :: time_in
  real, dimension(:),            allocatable :: depth_src, depth_dst
  real, dimension(:),            allocatable :: lon_src, lat_src
  real, dimension(:,:),          allocatable :: lon_dst, lat_dst
  real, dimension(:,:,:),        allocatable :: mask_dst
  !--- version information variables
  character(len=128) :: version='CVS $Id: regrid_3d.f90,v 14.0.2.1 2009/10/15 16:12:00 z1l Exp $'
  character(len=128) :: tagname='Tag $Name: mom4p1_pubrel_dec2009_nnz $'

  ! --- Begin of the program

  ! --- call fms_init, which will call mpp_init, mpp_io_init, mpp_domains_init
  call fms_init
  call constants_init

  !--- call regrid_3d initialization routine
  call regrid_3d_init

  !--- read the dest_grid file
  call read_dst_grid

  !--- read src_file
  call read_src_file

  !--- set up metadata of output file
  call setup_meta()

  !--- remap data from src grid to dest grid and write out data to output file
  call process_data()

  call regrid_3d_end

  call fms_end

contains

  !#####################################################################
  ! --- read the namelist and write the version and namelist to logfile. Also
  ! --- write the namelist to standard output
  subroutine regrid_3d_init

    integer :: io_status, unit, ierr, n

    D2R = PI/180.0
    ! --- the default src_field_name is 'temp and 'salt'
    src_field_name(1) = 'temp' 
    src_field_name(2) = 'salt'

    ! --- read namelist ------------------------------------------------
    if(file_exist('input.nml')) then
       unit = open_namelist_file()
       read (unit,regrid_3d_nml,IOSTAT=io_status)
       write (stdout(),'(/)')
       write (stdout(),regrid_3d_nml)  
       ierr = check_nml_error(io_status, 'regrid_3d_nml')
       call close_file(unit)
    else
       call mpp_error(FATAL, 'regrid_3d: file input.nml does not exist' )
    endif

    if(ntimes_saved > 0) then
       if(ntimes_saved > max_ntimes_saved) call mpp_error(FATAL, 'regrid_3d: nml ntimes_saved is greater than ' &
               //'max_ntimes_saved, increase max_ntimes_saved or decrease ntimes_saved')
       do n = 1, ntimes_saved
          if(timelevel_saved(n) .le. 0) call mpp_error(FATAL, 'regrid_3d: nml timelevel_saved is not ' &
              //'specified properly, modify timelevel_saved')
       enddo
    endif

    if(numfields .gt. max_fields) call mpp_error(FATAL, &
             'regrid_3d: numfields should be less than max_fields')
    if (numfields .le. 0) call mpp_error(FATAL,'regrid_3d: No field specified')
    !--- if dest_field_name is not defined in the namelist, get it from src_field_name
    do n = 1, numfields
       if(trim(dest_field_name(n)) == '') dest_field_name(n) = trim(src_field_name(n))
    enddo 
 
    !--- write version information
    call write_version_number(version, tagname)

  end subroutine regrid_3d_init

  !#####################################################################
  !--- open grid file and store grid info
  subroutine read_dst_grid

    integer                                    :: unit, ndim, nvar, natt, ntime, i, j, k
    integer                                    :: len1, siz_in(3)
    logical                                    :: found_xt, found_yt, found_kmt
    character(len=32)                          :: name
    real, allocatable, dimension(:,:)          :: kmt
    type(axistype), allocatable, dimension(:)  :: axes
    type(fieldtype), allocatable, dimension(:) :: fields

    if(.not. file_exist(trim(dest_grid)) ) &
         call mpp_error(FATAL, 'regrid_3d: file '//trim(dest_grid)//' does not exist')

    call mpp_open(unit, trim(dest_grid),&
         action=MPP_RDONLY, form=MPP_NETCDF, threading=MPP_MULTI, fileset=MPP_SINGLE)

    call mpp_get_info(unit, ndim, nvar, natt, ntime)

    allocate(fields(nvar), axes(ndim) )
    call mpp_get_axes(unit, axes)
    call mpp_get_fields(unit,fields)

    !--------------------------------------------------------------------
    ! get output grid information
    !--------------------------------------------------------------------
    ni_dst=0; nj_dst=0; nk_dst=0
    do i=1,ndim
       call mpp_get_atts(axes(i),name=name,len=len1)
       select case (trim(name))
       case ('grid_x_T')
          ni_dst = len1
       case ('grid_y_T')
          nj_dst = len1
       case ('zt')
          nk_dst = len1
          allocate(depth_dst(nk_dst))
          call mpp_get_axis_data(axes(i),depth_dst)
          if(.not.use_source_vertical_grid) depth_axis = axes(i)         
       end select
    enddo
    if(ni_dst==0) call mpp_error(FATAL,'regrid_3d: file '//trim(dest_grid)//' does not contain axis grid_x_T')
    if(nj_dst==0) call mpp_error(FATAL,'regrid_3d: file '//trim(dest_grid)//' does not contain axis grid_y_T')
    if(nk_dst==0) call mpp_error(FATAL,'regrid_3d: file '//trim(dest_grid)//' does not contain axis zt')

    allocate(lon_dst(ni_dst,nj_dst), lat_dst(ni_dst,nj_dst), kmt(ni_dst,nj_dst) )
    found_xt = .FALSE.;  found_yt = .FALSE.; found_kmt = .false.
    do i=1,nvar
       call mpp_get_atts(fields(i),name=name,ndim=ndim)
       select case (trim(name))
       case ('x_T')
          found_xt = .true.
          call mpp_read(unit,fields(i),lon_dst)
          field_lon_dst = fields(i)
          call mpp_get_atts(fields(i),axes=axes_dst)
       case ('y_T')
          found_yt = .true.
          call mpp_read(unit,fields(i),lat_dst)
          field_lat_dst = fields(i)           
       case ('num_levels')
          found_kmt = .true.
          call mpp_read(unit,fields(i),kmt )
       end select
    enddo
    if(.not.found_kmt) call mpp_error(FATAL,'regrid_3d: field num_levels is not in the file '//trim(dest_grid) )
    if(.not.found_xt) call mpp_error(FATAL,'regrid_3d: field x_T is not in the file '//trim(dest_grid) )
    if(.not.found_yt) call mpp_error(FATAL,'regrid_3d: field y_T is not in the file '//trim(dest_grid) )
    call mpp_close(unit)

    allocate(mask_dst(ni_dst,nj_dst,nk_dst) )

    do k = 1, nk_dst
       do j = 1, nj_dst
          do i = 1, ni_dst
             if(kmt(i,j) .ge. k) then
                mask_dst(i,j,k) = 1.0
             else
                mask_dst(i,j,k) = 0.0
             endif
          enddo
       enddo
    enddo
    
    deallocate(fields, axes, kmt)    

  end subroutine read_dst_grid

  !#####################################################################
  !--- read source grid informationfrom src_file
  subroutine read_src_file

    integer                                    :: unit, ndim, nvar, natt, n
    integer                                    :: nt, i, j, k, jj, len1
    logical                                    :: found_src_field(numfields)
    character(len=1)                           :: cart
    character(len=32)                          :: name, units
    type(axistype), allocatable, dimension(:)  :: axes
    type(fieldtype), allocatable, dimension(:) :: fields
    real                                       :: scale, add

    if(.not. file_exist(trim(src_file)) ) &
         call mpp_error(FATAL, 'regrid_3d: file '//trim(src_file)//' does not exist')

    call mpp_open(src_unit, trim(src_file),&
         action=MPP_RDONLY, form=MPP_NETCDF, threading=MPP_MULTI, fileset=MPP_SINGLE)
    call mpp_get_info(src_unit, ndim, nvar, natt, ntime_src)

    allocate(fields(nvar))
    call mpp_get_fields(src_unit, fields)
    allocate(src_field(numfields))
    found_src_field = .FALSE.
    do n=1,numfields
       do i=1,nvar
          call mpp_get_atts(fields(i),name=name)
          if (lowercase(trim(src_field_name(n))) == lowercase(trim(name))) then
             src_field(n) = fields(i)
             found_src_field(n) = .TRUE.
             write(stdout(),*) 'Interpolating src field : ',trim(name), ' grid ', trim(dest_grid)
          endif
       end do
    enddo
    do n=1,numfields 
       if(.not. found_src_field(n)) call mpp_error(FATAL, 'regrid_3d: field '&
            //trim(src_field_name(n))//' is not in the file '//trim(src_file) )
    enddo
    !--- get the src grid
    call mpp_get_atts(src_field(1),ndim=ndim)
    allocate(axes(ndim))
    call mpp_get_atts(src_field(1),axes=axes)
    ni_src=0; nj_src=0; nk_src=0; ntime_src = 1
    do j=1,ndim
       call mpp_get_atts(axes(j),len=len1,units=units)
       call get_axis_cart(axes(j),cart)
       select case (cart)
       case ('X')
          ni_src = len1
          allocate(lon_src(ni_src))
          call mpp_get_axis_data(axes(j),lon_src)
       case('Y')
          nj_src = len1
          allocate(lat_src(nj_src))
          call mpp_get_axis_data(axes(j),lat_src)
       case('Z')
          nk_src = len1
          allocate(depth_src(nk_src))
          if(use_source_vertical_grid) depth_axis = axes(j)
          call mpp_get_axis_data(axes(j),depth_src)
          if (trim(units) == 'cm') then
             depth_src(:) = depth_src(:)*.01
          endif
       case ('T')
          ntime_src = len1
          time_axis_exists = .true.
          allocate(time_in(ntime_src))
          call mpp_get_times(src_unit, time_in)
          time_axis = axes(j)
       end select
    enddo
    if(ni_src==0) call mpp_error(FATAL,'regrid_3d: file ' &
         //trim(src_file)//' does not contain axis with cartesian attributes = "X" ')
    if(nj_src==0) call mpp_error(FATAL,'regrid_3d: file '&
         //trim(src_file)//' does not contain axis with cartesian attributes = "Y" ')
    if(nk_src==0) call mpp_error(FATAL,'regrid_3d: file '&
         //trim(src_file)//' does not contain axis with cartesian attributes = "Z" ')
    !--- if use selective time level, make sure the time level selection is suitable.
    if(ntimes_saved > ntime_src) call mpp_error(FATAL, 'regrid_3d: nml ntime_src is greater than ' &
             //'number of time levels in the file '// trim(src_file) )
    do n = 1, ntimes_saved
       if(timelevel_saved(n) > ntime_src) call mpp_error(FATAL, 'regrid_3d: some entry in nml timelevel_saved ' &
             // 'is greater than number of time levels in the file '// trim(src_file) )
    enddo

    !--- get the missing value
    do n=1, numfields
       call mpp_get_atts(src_field(n),missing=missing(n),scale=scale, add=add )
       if(scale .NE. 1 .OR. add .NE. 0) missing(n) = missing(n)*scale + add
    enddo

    deallocate(fields, axes)    

  end subroutine read_src_file

  !#####################################################################
  !--- setup metadata of output file
  subroutine setup_meta

    integer            :: i, j, n, nt
    character(len=32)  :: units
    character(len=128) :: longname
    !--------------------------------------------------------------------
    ! write output file metadata
    !--------------------------------------------------------------------

    call mpp_open(dest_unit, trim(dest_file),MPP_OVERWR,MPP_NETCDF,threading=MPP_SINGLE,&
         fileset=MPP_SINGLE)
    !
    ! write axis metadata
    !
    do i=1,size(axes_dst(:))
       call mpp_copy_meta(dest_unit,axes_dst(i))
    end do

    call mpp_copy_meta(dest_unit,depth_axis)

    if (time_axis_exists) then
       call mpp_copy_meta(dest_unit, time_axis)
    endif

    !
    ! write variable metadata
    !
    call mpp_copy_meta(dest_unit,field_lon_dst,axes=axes_dst)
    call mpp_copy_meta(dest_unit,field_lat_dst,axes=axes_dst)

    do n = 1, numfields
       call mpp_get_atts(src_field(n),units=units,longname=longname)
       if (time_axis_exists) then 
          call mpp_write_meta(dest_unit, dest_field(n), (/axes_dst(1),axes_dst(2), depth_axis,time_axis/), &
               trim(dest_field_name(n)), units,longname, missing=missing(n), pack=1)
       else
          call mpp_write_meta(dest_unit, dest_field(n), (/axes_dst(1),axes_dst(2), depth_axis/), &
               trim(dest_field_name(n)), units,longname, missing=missing(n), pack=1)
       endif
    enddo
    ! write axis data

    do i=1,size(axes_dst(:))
       call mpp_write(dest_unit,axes_dst(i))
    end do

    call mpp_write(dest_unit,depth_axis)

    ! write variable data

    call mpp_write(dest_unit, field_lon_dst,lon_dst)
    call mpp_write(dest_unit, field_lat_dst,lat_dst)

  end subroutine setup_meta


  !#####################################################################
  !--- remap data from src grid to destination grid.
  subroutine process_data

    integer                             :: ndivs, i, j, k, n, nt, stat, ntimelevel, l
    integer                             :: isc, iec, jsc, jec, layout(2) = (/1,0/)
    real                                :: tmp_x, tmp_y
    type(domain2d)                      :: Domain
    type(horiz_interp_type)             :: Interp
    real, dimension(:,:,:), allocatable :: tmp1, tmp2, tmp, mask_src, data_dst, data_src
    real, dimension(:,:,:), allocatable :: depth_src_3d, depth_dst_3d

    !--- decompose model grid points
    !--- mapping can get expensive so we distribute the task at this level
    ndivs = mpp_npes()
    call mpp_define_layout ((/1,ni_dst,1,nj_dst/), ndivs, layout)
    call mpp_define_domains((/1,ni_dst,1,nj_dst/),layout, Domain,xhalo=0,yhalo=0)  
    call mpp_get_compute_domain (Domain, isc, iec, jsc, jec)

    call horiz_interp_new(Interp, lon_src*D2R, lat_src*D2R, lon_dst(isc:iec,jsc:jec)*D2R, &
         lat_dst(isc:iec,jsc:jec)*D2R, interp_method = trim(interp_method),                &
         num_nbrs = num_nbrs, max_dist=max_dist, grid_at_center = .true. )

    allocate(tmp1(ni_src,nj_src, nk_src))
    allocate(data_src(ni_src,nj_src, nk_src))
    if( use_source_vertical_grid) then 
       allocate(data_dst(ni_dst,nj_dst,nk_src), tmp(isc:iec,jsc:jec, nk_src) )
       call mpp_domains_set_stack_size(2*ni_dst*nj_dst*nk_src)
    else
       call mpp_domains_set_stack_size(2*ni_dst*nj_dst*nk_dst)
       allocate(tmp2(isc:iec,jsc:jec, nk_src) )
       allocate(data_dst(ni_dst,nj_dst,nk_dst), tmp(isc:iec,jsc:jec, nk_dst) )
       allocate(depth_src_3d(isc:iec,jsc:jec, nk_src))
       allocate(depth_dst_3d(isc:iec,jsc:jec,nk_dst))
       do k=1,nk_src
          depth_src_3d(:,:,k) = depth_src(k)
       enddo
       do k=1,nk_dst
          depth_dst_3d(:,:,k) = depth_dst(k)
       enddo
    endif


    if(ntimes_saved > 0) then
       ntimelevel = ntimes_saved       
    else
       ntimelevel = ntime_src
    endif

    write(stdout(),*)' There are ', ntimelevel, ' time steps to be saved to output file.'

    do l = 1, ntimelevel 
       if(ntimes_saved > 0) then
          nt = timelevel_saved(l)
       else
          nt = l
       endif
       write(stdout(),*)'**************At time step ', l       
       do n = 1, numfields
          !--- read source data
          call mpp_read(src_unit,src_field(n),data_src,nt)
          if( n==1 .and. nt == 1) then
             allocate(mask_src(ni_src,nj_src,nk_src))
             mask_src = 1.0
             do k = 1, nk_src
                do j = 1, nj_src
                   do i = 1, ni_src
                      if (abs(data_src(i,j,k) - missing(1)) <= tol) mask_src(i,j,k) = 0.
                      if (abs(data_src(i,j,k)) > max_val) mask_src(i,j,k) = 0.0
                   enddo
                enddo
             enddo
          endif
          !--- if scale_factor is not 1, multiple scale_factor to each data
          if(scale_factor(n) .ne. 1) then
             do k =1, nk_src
                do j =1, nj_src
                   do i =1, ni_src
                      if(mask_src(i,j,k) >0.5) data_src(i,j,k) = data_src(i,j,k)*scale_factor(n)
                   enddo
                enddo
             enddo
          endif

          !--- begin to regrid data
          if( use_source_vertical_grid) then 
             if(apply_mask) then
                do k = 1, nk_src
                   call horiz_interp(Interp, data_src(:,:,k), tmp(:,:,k),  &
                        mask_in = mask_src(:,:,k), missing_value=missing(n) )
                enddo
             else
                if(any(mask_src == 0.0) ) then ! do laplace extrap if needed.
                   call extrap(data_src(:,:,:), tmp1(:,:,:), stop_crit(n), &
                      missing(n)*scale_factor(n), is_cyclic ) 
                   do k = 1, nk_src 
                      call horiz_interp(Interp, tmp1(:,:,k), tmp(:,:,k))
                   enddo
                else
                   do k = 1, nk_src
                      call horiz_interp(Interp, data_src(:,:,k), tmp(:,:,k) )
                   enddo
                endif
             endif
          else
             if(any(mask_src == 0.0) ) then ! do laplace extrap if needed.
                call extrap(data_src(:,:,:), tmp1(:,:,:), stop_crit(n), &
                      missing(n)*scale_factor(n), is_cyclic ) 
                do k = 1, nk_src
                   call horiz_interp(Interp, tmp1(:,:,k), tmp2(:,:,k) )
                enddo
             else
                do k = 1, nk_src
                   call horiz_interp(Interp, data_src(:,:,k), tmp2(:,:,k) )
                enddo
             endif
                
             call interp_1d(depth_src_3d,depth_dst_3d,tmp2,tmp)

             do k = 1, nk_dst
                if(apply_mask) then
                   do j = jsc, jec
                      do i = isc, iec
                         if(mask_dst(i,j,k) < 0.5) tmp(i,j,k) = missing(n)
                      enddo
                   enddo
                endif
             enddo
          endif

          !--- get global data
          call mpp_global_field(Domain,tmp, data_dst(:,:,:) )
          if(mpp_pe()==mpp_root_pe()) then  ! --- write out data from root pe
             if (time_axis_exists) then 
                call mpp_write(dest_unit, dest_field(n), data_dst,time_in(nt))
             else
                call mpp_write(dest_unit, dest_field(n), data_dst)  
             endif

             if(debug) then            !--- the chksum is for debugging purpose

                write(stdout(),*)'NOTE: At time level ',nt,', Chksum for ',trim(src_field_name(n)), &
                     ' after-regrid data ', mpp_chksum(data_dst, (/mpp_root_pe()/) )
             endif
          endif
       enddo
    enddo

    deallocate(tmp, tmp1, data_src, data_dst, mask_src)
    if(.not. use_source_vertical_grid) deallocate(depth_src_3d, depth_dst_3d, tmp2 )
    !--- write out chksum for parallel checking


    call mpp_close(dest_unit)
    call mpp_close(src_unit)
    call horiz_interp_end

  end subroutine process_data

  !#####################################################################
  !--- release the memory
  subroutine regrid_3d_end

  deallocate(depth_src, depth_dst, lon_src, lat_src, lon_dst, lat_dst, mask_dst)
  if(time_axis_exists) deallocate(time_in)

  end subroutine regrid_3d_end


  !#####################################################################
  subroutine extrap(data_in, data_out, crit, missing_value, is_cyclic )
    real, dimension(:,:,:),  intent(in) :: data_in
    real, dimension(:,:,:), intent(out) :: data_out
    real,                    intent(in) :: crit, missing_value 
    logical,                 intent(in) :: is_cyclic 
    real                                :: resmax, initial_guess = 0.0
    integer                             :: ni, nj, nk, i, j, k, n

    real, dimension(0:size(data_in,1)+1, 0:size(data_in,2)+1) :: tmp
    real, dimension(size(data_in,1), size(data_in,2) )        :: sor, res

    real, dimension(size(data_in,1), size(data_in,2) )        :: cfn, cfs, cfe, cfw, cfc
    real                                                      :: latp, latm, cstr, csm, csj
    real, dimension(size(data_in,1))                          :: dxu, dxt
    real, dimension(size(data_in,2))                          :: dyu, dyt

    ni = size(data_in,1)
    nj = size(data_in,2)
    nk = size(data_in,3)

    ! construct grid factors for a sphere
    do j=1,nj-1
      dyu(j) = lat_src(j+1)-lat_src(j)
    enddo
    dyu(nj) = dyu(nj-1)
    do j=2,nj
      dyt(j) = 0.5*(dyu(j)+dyu(j-1))
    enddo
    dyt(1) = dyt(2)
    do i=1,ni-1
      dxu(i) = lon_src(i+1)-lon_src(i)
    enddo
    dxu(ni) = dxu(ni-1)
    do i=2,ni
      dxt(i) = 0.5*(dxu(i)+dxu(i-1))
    enddo
    dxt(1) = dxt(2)
    do j= 1, nj
      if (j == nj) then
        latp = lat_src(j) + 0.5*(lat_src(j) - lat_src(j-1))
      else
        latp = 0.5*(lat_src(j) + lat_src(j+1))
      endif
      if (j == 1) then
        latm = lat_src(j) - 0.5*(lat_src(j+1) - lat_src(j))
      else
        latm = 0.5*(lat_src(j) + lat_src(j-1))
      endif
      csj  = cos(latp*pi/180.0)
      csm  = cos(latm*pi/180.0)
      cstr = 1.0/cos(lat_src(j)*pi/180.)       
      do i= 1,ni
        cfn(i,j) = csj*cstr/(dyt(j)*dyu(j))
        cfs(i,j) = csm*cstr/(dyt(j)*dyu(max(j-1,1)))
        cfe(i,j) = cstr**2/(dxu(i)*dxt(i))
        cfw(i,j) = cstr**2/(dxu(max(i-1,1))*dxt(i))
        cfc(i,j) = 1.0/(cfn(i,j)+cfs(i,j)+cfe(i,j)+cfw(i,j))
        cfn(i,j) = cfn(i,j)*cfc(i,j)
        cfs(i,j) = cfs(i,j)*cfc(i,j)
        cfe(i,j) = cfe(i,j)*cfc(i,j)
        cfw(i,j) = cfw(i,j)*cfc(i,j)         
      enddo
    enddo

    tmp = 0.0
    do j= 1, nj
       do i= 1, ni
          if(abs(data_in(i,j,1) - missing_value) <= tol) then
             tmp(i,j) = initial_guess
          endif
       enddo
    enddo

    do k = 1, nk
       do j= 1, nj
          do i= 1, ni
             if(abs(data_in(i,j,k) - missing_value) <= tol ) then
                sor(i,j) = rel_coef
             else
                tmp(i,j) = data_in(i,j,k)
                sor(i,j) = 0.0
             endif
          enddo
       enddo

       call fill_boundaries(tmp, is_cyclic)

       ! iterate
       n=1
       do            
          resmax=0.0
          do j= 1, nj
             do i= 1, ni
                res(i,j) = cfw(i,j)*tmp(i-1,j) + cfe(i,j)*tmp(i+1,j) + cfs(i,j)*tmp(i,j-1) + cfn(i,j)*tmp(i,j+1) - tmp(i,j)
             enddo
          enddo

          do j= 1, nj
             do i= 1, ni
                res(i,j) = res(i,j)*sor(i,j)
                tmp(i,j)=tmp(i,j)+res(i,j)
                resmax = max(abs(res(i,j)),resmax)
             enddo
          enddo

          if(resmax .le. crit .or. n > max_iter) then
             data_out(:,:,k) = tmp(1:ni,1:nj)
             exit
          endif

          !--- update boundaries

          call fill_boundaries(tmp, is_cyclic)

          n=n+1

       enddo

       if(mpp_pe() == mpp_root_pe() ) write(stdout(),'(a,i6,a)') 'Stopped after ',n,' iterations'
       if(mpp_pe() == mpp_root_pe() ) write(stdout(),'(a,f10.4)') 'maxres= ',resmax
    enddo

  end subroutine extrap

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


  subroutine fill_boundaries(data, is_cyclic)

    real, dimension(0:,0:), intent(inout) :: data
    logical,                   intent(in) :: is_cyclic
    integer :: i,j, ni, nj

    ni = size(data,1) - 2
    nj = size(data,2) - 2

    if(is_cyclic) then
       data(0,1:nj) = data(ni,1:nj)
       data(ni+1,1:nj) = data(1,1:nj)
    endif

    return

  end subroutine fill_boundaries

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

end program regrid_3d