!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! 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
!-----------------------------------------------------------------------
! 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
!-----------------------------------------------------------------------
!
! Zhi Liang
!
!
! This program remap data from logically rectangular grid to logically rectangular grid.
!
!
! This program expects to read data from a netcdf file, which is specfied
! by the namelist variable "src_data". The number of data to be remapped is
! specified by num_flds. The name of field to be remapped is specified
! by the namelist variable "fld_name". The source data should be on the source
! grid which is specified by namelist variable src_grid. The destination grid is
! specified by nml dst_grd. The output file is a netcdf file specified
! by the namelist variable "dst_data". Each field can be a scalar variable or
! a vector field, which is specified by vector_fld. The vector field should
! always be paired together. A laplacian extrapolation will be performed when
! there is any missing value in the source data to interpolate data onto missing points.
!
use mpp_mod, only : mpp_npes, mpp_pe, mpp_root_pe, mpp_error
use mpp_mod, only : FATAL, WARNING, stdout, mpp_chksum
use mpp_domains_mod, only : mpp_define_layout, mpp_define_domains, domain2d
use mpp_domains_mod, only : mpp_global_field, mpp_domains_set_stack_size, mpp_get_compute_domain
use mpp_domains_mod, only : mpp_update_domains, CYCLIC_GLOBAL_DOMAIN, FOLD_NORTH_EDGE
use fms_mod, only : fms_init, open_namelist_file, close_file, fms_end
use fms_mod, only : check_nml_error, write_version_number, lowercase, uppercase
use fms_io_mod, only : fms_io_exit
use constants_mod, only : constants_init, PI
use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_type
use axis_utils_mod, only : interp_1d
implicit none
#include "netcdf.inc"
integer, parameter :: max_flds = 20
integer, parameter :: max_iter = 2000
integer, parameter :: max_atts = 20
real, parameter :: rel_coef = 0.6
real, parameter :: epsln = 1.0e-10
!--- namelist interface ----------------------------------------------
!
!
! Number of fields.
!
!
! Name of input file containing to-be-remapped data.
!
!
! Name of output file containing after-remapped data.
!
!
! Name of grid descriptor file containing target grid information.
!
!
! Name of grid descriptor file containing source grid information.
!
!
! Name of runoff field in input file.
!
!
! Name of grid where the field located. Valid choices are (T)racer, (C)orner, (E)ast and (N)orth.
!
!
! True if fields are vector components. All the vector field should be paired together.
! That is, if vector_field(n) is .true., then vector_field(n+1) should be true.
!
!
! The stopping criteria when extrapping data onto missing points.
!
!
! Number of nearest neighbors for regridding.
!
!
! Maximum region of influence around destination grid points.
!
!
! 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.
!
!
! 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.
!
!
! 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". Since bilinear interpolation suppose the source grid
! is a lat-lon grid, "spherical" need to be used if the source grid is not a lat-lon grid.
!
!
! For Debugging. Set true to print out chksum information for debugging reproducing ability
! accross processors. default is false.
!
!
integer :: num_flds = 0
character(len=128) :: src_data = ''
character(len=128) :: src_grid = ''
character(len=128) :: dst_data = ''
character(len=128) :: dst_grid = ''
character(len=16) :: fld_name(max_flds) = ''
character(len=1) :: fld_pos(max_flds) = '' ! with value 'T','C','E','N'
logical :: vector_fld(max_flds) = .false. ! True if fields are vector components.
real :: stop_crit(max_flds) = 0.001
integer :: num_nbrs = 5
real :: max_dist = 0.1
logical :: use_source_vertical_grid = .FALSE.
logical :: apply_mask = .TRUE.
character(len=32) :: interp_method = "bilinear"
logical :: debug = .false.
namelist /regrid_nml/ num_flds, src_data, src_grid, dst_data, dst_grid, fld_name, &
fld_pos, vector_fld, num_nbrs, max_dist, stop_crit, interp_method, &
apply_mask, use_source_vertical_grid, debug
!--- data type for easy data management ------------------------------
type cell_type
real, dimension(:,:), pointer :: geolon=>NULL(), geolat=>NULL() ! geographical grid
real, dimension(:,:,:), pointer :: mask =>NULL() ! mask at each vertical level
real, dimension(:,:), pointer :: sinrot=>NULL(), cosrot=>NULL() ! rotation
end type cell_type
type regrid_type
type(cell_type) :: T,C,E,N ! T,C,E,N-cell
integer :: ni, nj,nk ! grid size
integer :: isc, iec, jsc, jec ! compute domain decompsition src grid
integer :: isd, ied, jsd, jed ! data domain decompsition on dst grid
type(domain2d) :: domain ! domain of the grid ( only destination grid will use it).
logical :: is_cyclic ! indicate if the grid is cyclic
logical :: is_tripolar ! indicate if the grid is tripolar
real, pointer :: zt(:) =>NULL() ! vertical grid
real, pointer :: grid_xt(:)=>NULL() ! T-cell longitude grid
real, pointer :: grid_yt(:)=>NULL() ! T-cell latitude grid
real, pointer :: grid_xc(:)=>NULL() ! C-cell longitude grid
real, pointer :: grid_yc(:)=>NULL() ! C-cell latitude grid
end type regrid_type
!--- version information ---------------------------------------------
character(len=128) :: version = '$Id: regrid.F90,v 16.0 2008/07/30 22:49:37 fms Exp $'
character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $'
!--- other variables
integer :: ntime_src ! number of time levels of src data
logical :: do_vertical_interp ! indicate if vertical interpolation is needed.
logical :: is_root_pe ! if current pe is root pe.
type(regrid_type) :: Src ! regrid_type data of source
type(regrid_type) :: Dst ! regrid_type data of destination
type(horiz_interp_type), target :: Interp_T, Interp_C ! horiz_interp data at T,C-cell
type(horiz_interp_type), target :: Interp_E, Interp_N ! horiz_interp data at E,N-cell
real, dimension(:), allocatable :: missing_value ! missing value for source field.
real, dimension(:), allocatable :: time_value ! time axis value
real, dimension(:), allocatable :: t1_avg, t2_avg, dt_avg ! time average data
real, dimension(:,:), allocatable :: ht, hc, he, hn ! topography
logical, dimension(:), allocatable :: has_zaxis ! indicate if the field has z axis
logical :: has_taxis,tavg_exist ! indicate if all the field has time axis
integer :: ncid_dst ! ncid corresponding to
integer :: id_t1_dst, id_t2_dst, id_dt_dst
integer :: id_time_dst, id_dst_fld(max_flds)
!--- begin of the program -------------------------------------------
call regrid_init()
call read_grid()
call setup_meta()
call process_data ()
call regrid_end ()
contains
!#####################################################################
subroutine regrid_init
integer :: unit, io_status, ierr, n
character(len=1) :: direction
call fms_init()
call constants_init()
!--- reading namelist ------------------------------------------------
unit = open_namelist_file()
read (unit, regrid_nml,iostat=io_status)
write (stdout(),'(/)')
write (stdout(), regrid_nml)
ierr = check_nml_error(io_status,'regrid_nml')
call close_file (unit)
!--- write out version information ---------------------------------
call write_version_number( version, tagname )
if(num_flds == 0) call error_handler('regrid: nml num_fiels = 0 should be a positive number')
if(num_flds .gt. max_flds) call error_handler('regrid: nml num_fiels is greater than maximum'// &
' allowed fields max_flds. Modify num_flds or increase "max_flds" at top of the program' )
do n = 1, num_flds
if(fld_pos(n) .ne. 'T' .and. fld_pos(n) .ne. 'C' .and. fld_pos(n) .ne. 'E' .and. fld_pos(n) .ne. 'N' ) &
call error_handler('regrid: fld_pos should be "T", "C", "E" or "N" ')
enddo
!--- number of vector_fld is true should be even and should be paired adjacent.
!--- The first one of the pair should be u-component and the second one
!--- should be v-component.
n = 1
do while ( n .le. num_flds )
if(vector_fld(n)) then
n = n + 1
if( n .gt. num_flds .or. ( .NOT. vector_fld(n) ) ) then
call error_handler('regrid: vector field should be paired together')
endif
endif
n = n+1
enddo
!--- write out field information to be remapped.
Write(stdout(),*) '*****************************************************************'
write(stdout(),*) 'number of fields to be remapped are: ', num_flds
write(stdout(),*) 'The source data file is: ', trim(src_data)
write(stdout(),*) 'The source grid file is: ', trim(src_grid)
write(stdout(),*) 'The destination grid file is: ', trim(dst_grid)
if(use_source_vertical_grid) then
write(stdout(),*) 'use_source_vertical_grid is set to true. The destination data will have '// &
'same vertical level as the source data '
else
write(stdout(),*) 'use_source_vertical_grid is set to false. vertical interpolation will ', &
'be performed if the source and destination has different vertical grid. '
endif
direction = 'x'
do n = 1, num_flds
if(vector_fld(n)) then
write(stdout(),*) '** field ',n,': ',trim(fld_name(n)), ' at '// &
trim(fld_pos(n))//'-cell center is the '//direction//'-component of a vector field'
if(direction == 'x') then
direction = 'y'
else
direction = 'x'
endif
else
write(stdout(),*) '**field ',n,': ',trim(fld_name(n)), ' at '// &
trim(fld_pos(n))//'-cell center is a scalar variable'
endif
enddo
!--- memory allocation
allocate(has_zaxis(num_flds), missing_value(num_flds) )
has_zaxis = .false.
missing_value = -1e20
is_root_pe = mpp_pe() == mpp_root_pe()
end subroutine regrid_init
!#################################################################
!--- read the src grid and dst grid and land/sea mask of dst grid.
subroutine read_grid
integer :: isc, iec, jsc, jec
integer :: k, npes, layout(2)
real :: D2R
logical :: use_src_vgrid
D2R = PI/180.0
!--- read some general grid informaiton for source grid and destination grid
call get_grid_info(src_grid, Src, .TRUE.)
call get_grid_info(dst_grid, Dst, .FALSE.)
do_vertical_interp = .false.
if(use_source_vertical_grid) then !--- modify the destination vertical information for the future use.
!--- In this case the destination vertical grid should be the same as source.
Dst%nk = Src%nk
deallocate(Dst%zt)
allocate(Dst%zt(Dst%nk))
Dst%zt = Src%zt
else !--- when use_source_vertical_grid is false, vertical interpolation will be done if the source vertical
!--- grid didn't match destination vertical grid.
if(Src%nk .ne. Dst%nk) then
do_vertical_interp = .true.
else
do k = 1, Src%nk
if( abs(Src%zt(k) - Dst%zt(k)) .gt. epsln ) then
do_vertical_interp = .true.
exit
endif
enddo
endif
if(do_vertical_interp) then
write(stdout(),*)'NOTE: Since source vertical grid does not match ' //&
'destination vertical grid and nml use_source_vertical_grid is set ' //&
'to false, vertical interpolation will be performed. '
else
write(stdout(),*)'NOTE: Even though nml use_source_vertical_grid is set to false, ', &
'but since source vertical grid match destination vertical grid, ', &
'no vertical interpolation will be performed. '
endif
endif
isc= Dst%isc; iec= Dst%iec; jsc= Dst%jsc; jec= Dst%jec;
!--- when source is tripolar grid, 'spherical' interpolation will be used
!--- no matter your choice of "interp_method".
if(Src%is_tripolar .and. trim(interp_method) == "bilinear") then
write(stdout(),*) "WARNING: Since src grid is tripolar grid, even though you chose ", &
"bilinear option, spherical option will be used."
endif
!--- call horizinterp_new to calculate the reampping weight in horizontal direction
if( any(fld_pos == 'T') ) then
if(Src%is_tripolar) then
call horiz_interp_new(Interp_T, Src%T%geolon*D2R, Src%T%geolat*D2R, &
Dst%T%geolon(isc:iec,jsc:jec)*D2R, Dst%T%geolat(isc:iec,jsc:jec)*D2R, &
interp_method='spherical', num_nbrs=num_nbrs, max_dist=max_dist, &
src_modulo = Src%is_cyclic)
else
call horiz_interp_new(Interp_T, Src%T%geolon(:,1)*D2R, Src%T%geolat(1,:)*D2R, &
Dst%T%geolon(isc:iec,jsc:jec)*D2R, Dst%T%geolat(isc:iec,jsc:jec)*D2R, &
interp_method=trim(interp_method), num_nbrs=num_nbrs, max_dist=max_dist, &
src_modulo = Src%is_cyclic, grid_at_center = .true.)
endif
endif
if( any(fld_pos == 'C') ) then
if(Src%is_tripolar) then
call horiz_interp_new(Interp_C, Src%C%geolon*D2R, Src%C%geolat*D2R, &
Dst%C%geolon(isc:iec,jsc:jec)*D2R, Dst%C%geolat(isc:iec,jsc:jec)*D2R, &
interp_method='spherical', num_nbrs=num_nbrs, max_dist=max_dist, &
src_modulo = Src%is_cyclic )
else
call horiz_interp_new(Interp_C, Src%C%geolon(:,1)*D2R, Src%C%geolat(1,:)*D2R, &
Dst%C%geolon(isc:iec,jsc:jec)*D2R, Dst%C%geolat(isc:iec,jsc:jec)*D2R, &
interp_method=trim(interp_method), num_nbrs=num_nbrs, max_dist=max_dist, &
src_modulo = Src%is_cyclic, grid_at_center = .true.)
endif
endif
if( any(fld_pos == 'E') ) then
if(Src%is_tripolar) then
call horiz_interp_new(Interp_E, Src%E%geolon*D2R, Src%E%geolat*D2R, &
Dst%E%geolon(isc:iec,jsc:jec)*D2R, Dst%E%geolat(isc:iec,jsc:jec)*D2R, &
interp_method='spherical', num_nbrs=num_nbrs, max_dist=max_dist, &
src_modulo = Src%is_cyclic)
else
call horiz_interp_new(Interp_E, Src%E%geolon(:,1)*D2R, Src%E%geolat(1,:)*D2R, &
Dst%E%geolon(isc:iec,jsc:jec)*D2R, Dst%E%geolat(isc:iec,jsc:jec)*D2R, &
interp_method=trim(interp_method), num_nbrs=num_nbrs, max_dist=max_dist, &
src_modulo = Src%is_cyclic, grid_at_center = .true.)
endif
endif
if( any(fld_pos == 'N') ) then
if(Src%is_tripolar) then
call horiz_interp_new(Interp_N, Src%N%geolon*D2R, Src%N%geolat*D2R, &
Dst%N%geolon(isc:iec,jsc:jec)*D2R, Dst%N%geolat(isc:iec,jsc:jec)*D2R, &
interp_method='spherical', num_nbrs=num_nbrs, max_dist=max_dist, &
src_modulo = Src%is_cyclic)
else
call horiz_interp_new(Interp_N, Src%N%geolon(:,1)*D2R, Src%N%geolat(1,:)*D2R, &
Dst%N%geolon(isc:iec,jsc:jec)*D2R, Dst%N%geolat(isc:iec,jsc:jec)*D2R, &
interp_method=trim(interp_method), num_nbrs=num_nbrs, max_dist=max_dist, &
src_modulo = Src%is_cyclic, grid_at_center = .true.)
endif
endif
end subroutine read_grid
!#####################################################################
subroutine get_grid_info(gridfile, Regrid, is_src_grid)
character(len=*), intent(in) :: gridfile
type(regrid_type), intent(inout) :: Regrid
logical, intent(in) :: is_src_grid
integer :: rcode, ncid, natt, id_zt, id_kmt, dims(2)
integer :: id_grid_xt, id_grid_yt, id_grid_xc, id_grid_yc
integer :: npes, i, j, k, layout(2)
integer :: isc, iec, jsc, jec, isd, ied, jsd, jed
logical :: is_new_grid
character(len=128) :: name, catt
integer,allocatable,dimension(:,:) :: kmt, kmc, kme, kmn
real, allocatable,dimension(:,:) :: tmp
rcode = nf_open(trim(gridfile), NF_NOWRITE, ncid)
call error_handler('error in open file '//trim(gridfile), rcode)
Regrid%ni = 0; Regrid%nj = 0; Regrid%nk = 0
!--- get vertical grid information
rcode = nf_inq_varid(ncid, 'zt', id_zt)
call error_handler('error in inquring id of field zt from file '//trim(gridfile), rcode)
rcode = nf_inq_vardimid(ncid, id_zt, dims)
call error_handler('error in inquring dims of field zt from file '//trim(gridfile), rcode)
rcode = nf_inq_dimlen(ncid, dims(1), Regrid%nk)
call error_handler('error in inquring dimlen of field zt from file '//trim(gridfile), rcode)
allocate(Regrid%zt(Regrid%nk) )
rcode = nf_get_var_double(ncid, id_zt, Regrid%zt)
call error_handler('error in getting data of zt from file '//trim(gridfile), rcode)
!--- get horizontal grid size
rcode = nf_inq_varid(ncid, 'num_levels', id_kmt)
if(rcode == 0) then
is_new_grid = .true.
else
is_new_grid = .false.
rcode = nf_inq_varid(ncid, 'kmt', id_kmt)
endif
call error_handler('Can not find kmt/num_levels in file '//trim(gridfile), rcode)
rcode = nf_inq_vardimid(ncid, id_kmt, dims)
call error_handler('error in inquring dims of field num_levels/kmt from file '//trim(gridfile), rcode)
rcode = nf_inq_dimlen(ncid, dims(1), Regrid%ni )
call error_handler('error in inquring dimlen ni of field num_levels/kmt from file '//trim(gridfile), rcode)
rcode = nf_inq_dimlen(ncid, dims(2), Regrid%nj )
call error_handler('error in inquring dimlen nj of field num_levels/kmt from file '//trim(gridfile), rcode)
!--- get global attributes ( cyclic or tripolar )
rcode = nf_inq_natts(ncid, natt)
call error_handler('error in inquring natts of file '//trim(gridfile), rcode )
do i = 1, natt
rcode = nf_inq_attname(ncid,NF_GLOBAL,i,name )
call error_handler('error in inquring att name of file '//trim(gridfile), rcode )
catt = ''
select case( trim(name) )
case('x_boundary_type')
rcode = nf_get_att_text(ncid,NF_GLOBAL,name,catt)
call error_handler('error in inquring x_boundary_type value of file '//trim(gridfile), rcode )
if(trim(catt) == 'cyclic') then
Regrid%is_cyclic = .true.
Write(stdout(),*)' NOTE: x_boundary_type of grid '//trim(gridfile)//' is cyclic'
else
Write(stdout(),*)' NOTE: x_boundary_type of '//trim(gridfile)//' is solid_walls'
endif
case ('y_boundary_type')
rcode = nf_get_att_text(ncid,NF_GLOBAL,name,catt)
call error_handler('error in inquring y_boundary_type value of file '//trim(gridfile), rcode )
if (trim(catt) == 'fold_north_edge') then
Regrid%is_tripolar = .true.
Write(stdout(),*)' NOTE: y_boundary_type of '//trim(gridfile)//' is tripolar'
else
Write(stdout(),*)' NOTE: y_boundary_type of '//trim(gridfile)//' is solid_walls'
endif
end select
end do
!--- define domain decompsition for destination grid ---------------
if(.not. is_src_grid) then
npes = mpp_npes()
layout = (/1,0/)
call mpp_define_layout((/1,Regrid%ni,1,Regrid%nj/),npes,layout)
if(Regrid%is_tripolar) then
call mpp_define_domains((/1,Regrid%ni,1,Regrid%nj/),layout, Regrid%domain, xflags = CYCLIC_GLOBAL_DOMAIN, &
yflags = FOLD_NORTH_EDGE, xhalo=1, yhalo=1)
else if (Regrid%is_cyclic) then
call mpp_define_domains((/1,Regrid%ni,1,Regrid%nj/),layout, Regrid%domain, xflags = cyclic_global_domain, &
xhalo=1, yhalo=1)
else
call mpp_define_domains((/1,Regrid%ni,1,Regrid%nj/),layout, Regrid%domain, xhalo=1, yhalo=1)
endif
call mpp_get_compute_domain(Regrid%domain,Regrid%isc,Regrid%iec,Regrid%jsc,Regrid%jec)
endif
!--- when the grid is source grid, need to get global grid information.
!--- when the grid is destination grid, only need to get local grid information.
if(is_src_grid) then
isc = 1; iec = Regrid%ni
jsc = 1; jec = Regrid%nj
else
isc = Regrid%isc ; iec = Regrid%iec
jsc = Regrid%jsc ; jec = Regrid%jec
endif
isd = isc-1 ; ied = iec+1
jsd = jsc-1 ; jed = jec+1
! get number vertical level
allocate(kmt(isd:ied,jsd:jed), kmc(isc:iec,jsc:jec), kme(isc:iec,jsc:jec), kmn(isc:iec,jsc:jec) )
allocate(tmp(Regrid%ni,Regrid%nj) )
rcode = nf_get_var_double(ncid, id_kmt, tmp)
call error_handler('error in getting value of kmt from file '//trim(gridfile), rcode)
kmt = 0
kmt(isc:iec,jsc:jec) = tmp(isc:iec,jsc:jec)
if(is_src_grid) then
kmt(isd,jsc:jec) = kmt(iec,jsc:jec)
kmt(ied,jsc:jec) = kmt(isc,jsc:jec)
else
call mpp_update_domains(kmt,Dst%domain)
endif
!--- calculate the mask at T,C,E,N-cell if needed
if( any(fld_pos == 'T') ) then
allocate(Regrid%T%mask(isc:iec,jsc:jec,Regrid%nk))
do k=1,Regrid%nk
do j=jsc,jec
do i=isc,iec
if (kmt(i,j) .ge. k) then
Regrid%T%mask(i,j,k) = 1.0
else
Regrid%T%mask(i,j,k) = 0.0
endif
enddo
enddo
enddo
endif
if( any(fld_pos == 'C') ) then
do j = jsc, jec
do i = isc, iec
kmc(i,j) = min(kmt(i,j), kmt(i+1,j), kmt(i,j+1), kmt(i+1,j+1))
enddo
enddo
allocate(Regrid%C%mask(isc:iec,jsc:jec,Regrid%nk))
do k=1,Regrid%nk
do j=jsc,jec
do i=isc,iec
if (kmc(i,j) .ge. k) then
Regrid%C%mask(i,j,k) = 1.0
else
Regrid%C%mask(i,j,k) = 0.0
endif
enddo
enddo
enddo
endif
if( any(fld_pos == 'E') ) then
do j = jsc, jec
do i = isc, iec
kme(i,j) = min(kmt(i,j), kmt(i+1,j))
enddo
enddo
allocate(Regrid%E%mask(isc:iec,jsc:jec,Regrid%nk))
do k=1,Regrid%nk
do j=jsc,jec
do i=isc,iec
if (kme(i,j) .ge. k) then
Regrid%E%mask(i,j,k) = 1.0
else
Regrid%E%mask(i,j,k) = 0.0
endif
enddo
enddo
enddo
endif
if( any(fld_pos == 'N') ) then
do j = jsc, jec
do i = isc, iec
kmn(i,j) = min(kmt(i,j), kmt(i,j+1))
enddo
enddo
allocate(Regrid%N%mask(isc:iec,jsc:jec,Regrid%nk))
do k=1,Regrid%nk
do j=jsc,jec
do i=isc,iec
if (kmn(i,j) .ge. k) then
Regrid%N%mask(i,j,k) = 1.0
else
Regrid%N%mask(i,j,k) = 0.0
endif
enddo
enddo
enddo
endif
!--- get the destination axis grid information ---------------------------------
if(.not. is_src_grid) then
allocate(Regrid%grid_xt(Regrid%ni),Regrid%grid_yt(Regrid%nj) )
allocate(Regrid%grid_xc(Regrid%ni),Regrid%grid_yc(Regrid%nj) )
if(is_new_grid) then
rcode = nf_inq_varid(ncid, 'grid_x_T', id_grid_xt)
call error_handler('error in inquring variable id of field grid_x_T', rcode )
rcode = nf_inq_varid(ncid, 'grid_y_T', id_grid_yt)
call error_handler('error in inquring variable id of field grid_y_T', rcode )
rcode = nf_inq_varid(ncid, 'grid_x_C', id_grid_xc)
call error_handler('error in inquring variable id of field grid_x_C', rcode )
rcode = nf_inq_varid(ncid, 'grid_y_C', id_grid_yc)
call error_handler('error in inquring variable id of field grid_y_C', rcode )
else
rcode = nf_inq_varid(ncid, 'gridlon_t', id_grid_xt)
call error_handler('error in inquring variable id of field gridlon_t', rcode )
rcode = nf_inq_varid(ncid, 'gridlat_t', id_grid_yt)
call error_handler('error in inquring variable id of field gridlat_t', rcode )
rcode = nf_inq_varid(ncid, 'gridlon_c', id_grid_xc)
call error_handler('error in inquring variable id of field gridlon_c', rcode )
rcode = nf_inq_varid(ncid, 'gridlat_c', id_grid_yc)
call error_handler('error in inquring variable id of field gridlat_c', rcode )
endif
rcode = nf_get_var_double(ncid, id_grid_xt, Regrid%grid_xt)
call error_handler('error in getting data of variable grid_x_T/gridlon_t', rcode )
rcode = nf_get_var_double(ncid, id_grid_yt, Regrid%grid_yt)
call error_handler('error in getting data of variable grid_y_T/gridlat_t', rcode )
rcode = nf_get_var_double(ncid, id_grid_xc, Regrid%grid_xc)
call error_handler('error in getting data of variable grid_x_C/gridlon_c', rcode )
rcode = nf_get_var_double(ncid, id_grid_yc, Regrid%grid_yc)
call error_handler('error in getting data of variable grid_y_C/gridlat_c', rcode )
endif
if(any(fld_pos == 'T') ) then
call get_cell_info(Regrid%T,'T',ncid,Regrid%ni,Regrid%nj,isc,iec,jsc,jec,is_new_grid)
endif
if(any(fld_pos == 'C') ) then
call get_cell_info(Regrid%C,'C',ncid,Regrid%ni,Regrid%nj,isc,iec,jsc,jec,is_new_grid)
endif
if(any(fld_pos == 'E') ) then
call get_cell_info(Regrid%E,'E',ncid,Regrid%ni,Regrid%nj,isc,iec,jsc,jec,is_new_grid)
endif
if(any(fld_pos == 'N') ) then
call get_cell_info(Regrid%N,'N',ncid,Regrid%ni,Regrid%nj,isc,iec,jsc,jec,is_new_grid)
endif
deallocate(kmt, kmc, kme, kmn,tmp)
rcode = nf_close(ncid)
end subroutine get_grid_info
!#####################################################################
subroutine get_cell_info(Cell, type, ncid, ni, nj, isc, iec,jsc,jec,is_new_grid )
type(cell_type), intent(inout) :: Cell
character(len=1), intent(in) :: type
integer, intent(in) :: ncid, ni, nj
integer, intent(in) :: isc, iec, jsc, jec
logical, intent(in) :: is_new_grid
character(len=1) :: lc_type, uc_type
integer :: rcode, id_x, id_y, id_angle, id_sinrot, id_cosrot
real, allocatable :: tmp(:,:)
allocate(Cell%geolon(isc:iec, jsc:jec), Cell%geolat(isc:iec, jsc:jec) )
allocate(Cell%sinrot(isc:iec, jsc:jec), Cell%cosrot(isc:iec, jsc:jec) )
allocate(tmp(ni,nj) )
lc_type = lowercase(type)
uc_type = uppercase(type)
if(is_new_grid) then
rcode = nf_inq_varid(ncid, 'x_'//uc_type, id_x)
call error_handler('error in inquring variable id of field x_'//uc_type, rcode )
rcode = nf_inq_varid(ncid, 'y_'//uc_type, id_y)
call error_handler('error in inquring variable id of field y_'//uc_type, rcode )
rcode = nf_inq_varid(ncid, 'angle_'//uc_type, id_angle)
call error_handler('error in inquring variable id of field angle_'//uc_type, rcode )
else
rcode = nf_inq_varid(ncid, 'geolon_'//lc_type, id_x)
call error_handler('error in inquring variable id of field geolon_'//lc_type, rcode )
rcode = nf_inq_varid(ncid, 'geolat_'//lc_type, id_y)
call error_handler('error in inquring variable id of field geolat_'//lc_type, rcode )
if(lc_type == 'c') then
rcode = nf_inq_varid(ncid, 'sin_rot', id_sinrot)
call error_handler('error in inquring variable id of field sin_rot', rcode )
rcode = nf_inq_varid(ncid, 'cos_rot', id_cosrot)
call error_handler('error in inquring variable id of field cos_rot', rcode )
else
write(stdout(),*)'==> NOTE: No rotation information avaible for '// &
uc_type//'-cell, will set sinrot = 0 and cosrot = 1'
endif
endif
rcode = nf_get_var_double(ncid, id_x, tmp)
call error_handler('error in getting data of variable x_'//uc_type//'/geolon_'//lc_type, rcode )
Cell%geolon(isc:iec, jsc:jec) = tmp(isc:iec,jsc:jec)
rcode = nf_get_var_double(ncid, id_y, tmp)
call error_handler('error in getting data of variable y_'//uc_type//'/geolat_'//lc_type, rcode )
Cell%geolat(isc:iec, jsc:jec) = tmp(isc:iec,jsc:jec)
if(is_new_grid) then
rcode = nf_get_var_double(ncid, id_angle, tmp)
call error_handler('error in getting data of variable angle_'//uc_type, rcode )
Cell%sinrot(isc:iec,jsc:jec) = sin(tmp(isc:iec,jsc:jec))
Cell%cosrot(isc:iec,jsc:jec) = cos(tmp(isc:iec,jsc:jec))
else
if(lc_type == 'c') then
rcode = nf_get_var_double(ncid, id_sinrot, tmp)
call error_handler('error in getting data of variable sin_rot', rcode )
Cell%sinrot(isc:iec,jsc:jec) = tmp(isc:iec,jsc:jec)
rcode = nf_get_var_double(ncid, id_cosrot, tmp)
call error_handler('error in getting data of variable cos_rot', rcode )
Cell%cosrot(isc:iec,jsc:jec) = tmp(isc:iec,jsc:jec)
else
Cell%sinrot = 0
Cell%cosrot = 1
endif
endif
deallocate(tmp)
end subroutine get_cell_info
!#####################################################################
!--- set up axis and field meta data for destination data file
subroutine setup_meta
integer :: ncid_src, rcode, id_fld, id_dim, id_time, dims(4)
integer :: type, len, ndims, natt_dim, natt, dimids(4), ndims_dst
integer :: n, m, i, j, num_has_taxis, start(4), nwrite(4)
integer :: id_t1_src, id_t2_src, id_dt_src, dim_time_dst
integer, parameter :: max_dim = 10
integer :: dims_dst(max_dim), id_axes(max_dim), fld_dims(4), dims_id(max_dim)
character(len=1) :: dims_pos(max_dim), dims_cart(max_dim)
character(len=128) :: name, att_name, cart, dims_name(max_dim)
character(len=512) :: catt
logical :: found_dim = .false.
logical :: has_missing(max_flds)
!--- open the dst_data file from root_pe.
if(is_root_pe) then
rcode = nf_create(trim(dst_data),NF_WRITE, ncid_dst)
call error_handler('error in creating file '//trim(dst_data), rcode)
endif
!--- get src_data information ( ntime, global attributes, axis and field meta information )
rcode = nf_open(trim(src_data), NF_NOWRITE, ncid_src)
call error_handler('error in opening file '//trim(src_data), rcode)
!--- get missing value, tavg_information, has_zaxis or has_taxis
num_has_taxis = 0
ndims_dst = 0
tavg_exist = .false.
ntime_src = 1
do n = 1, num_flds
rcode = nf_inq_varid(ncid_src,fld_name(n), id_fld)
call error_handler('error in inquring id of field '//trim(fld_name(n))// ' of file '//trim(src_data), rcode)
rcode = nf_inq_varndims( ncid_src, id_fld, ndims )
call error_handler('error in inquring ndims of field '//trim(fld_name(n))// ' of file '//trim(src_data), rcode)
rcode = nf_inq_vardimid( ncid_src, id_fld, dimids )
call error_handler('error in inquring dimids of field '//trim(fld_name(n))// ' of file '//trim(src_data), rcode)
rcode = nf_inq_varnatts(ncid_src, id_fld, natt )
call error_handler('error in inquring natt of field '//trim(fld_name(n))// ' of file '//trim(src_data), rcode)
!--- to check if each field has z or time axis
do m = 1, ndims
rcode = nf_inq_dimname ( ncid_src, dimids(m), name )
call error_handler('error in inquring dimension name of dimid', rcode)
rcode = nf_inq_varid(ncid_src, name, id_dim)
call error_handler('error in inquring variable id of '//trim(name), rcode)
rcode = nf_inq_varnatts(ncid_src, id_dim, natt_dim)
call error_handler('error in inquring number of attributes of variable '//trim(name), rcode)
if(lowercase(trim(name)) == 'time') then
cart = 'T'
if(num_has_taxis == 0) then
rcode = nf_inq_dimlen(ncid_src, dimids(m), ntime_src)
call error_handler('error in inquring time dimension from file '//trim(src_data), rcode)
allocate(time_value(ntime_src) )
rcode = nf_get_var_double(ncid_src, id_dim, time_value)
call error_handler('error in inquring time value from file '//trim(src_data), rcode)
endif
num_has_taxis = num_has_taxis + 1
else
do i = 1, natt_dim
rcode = nf_inq_attname(ncid_src, id_dim,i,att_name )
call error_handler('error in inquring attribute name ', rcode)
if(trim(att_name) == 'cartesian_axis' ) then
rcode = nf_get_att_text(ncid_src,id_dim,att_name,cart)
call error_handler('error in get attribute value of attribute '//trim(att_name), rcode)
select case(cart(1:1))
case('X','Y')
! do nothing
case('Z')
has_zaxis(n) = .true.
case default
call error_handler('The cartesian_axis of each axis should be X, Y or Z ' )
end select
exit
endif
enddo
endif
found_dim = .false.
do j = 1, ndims_dst
if(trim(dims_name(j)) == trim(name)) then
found_dim = .true.
exit
endif
enddo
if(.not. found_dim) then
ndims_dst = ndims_dst + 1
dims_name(ndims_dst) = trim(name)
dims_pos(ndims_dst) = fld_pos(n)
dims_cart(ndims_dst) = cart(1:1)
dims_id(ndims_dst) = id_dim
endif
enddo
has_missing(n) = .false.
do i = 1, natt
rcode = nf_inq_attname(ncid_src, id_fld,i,name )
call error_handler('error in inquring attribute name', rcode)
if(trim(name) == 'time_avg_info') then
tavg_exist = .true.
else if(trim(name) == 'missing_value') then
has_missing(n) = .true.
rcode = nf_get_att_double(ncid_src, id_fld, name, missing_value(n) )
endif
enddo
enddo
if(num_has_taxis == 0) then
has_taxis = .false.
else if(num_has_taxis == num_flds) then
has_taxis = .true.
else
call error_handler('Either all or none of the fields has time axis ' )
endif
!--- suppose the time_avg_info is specified by average_T1, average_T2, average_DT.
if(tavg_exist) then
allocate( t1_avg(ntime_src), t2_avg(ntime_src), dt_avg(ntime_src) )
rcode = nf_inq_varid(ncid_src, 'average_T1', id_t1_src)
call error_handler('error in inquring variable id of average_T1', rcode)
rcode = nf_get_var_double(ncid_src,id_t1_src, t1_avg )
call error_handler('error in getting average_T1 data', rcode)
rcode = nf_inq_varid(ncid_src, 'average_T2', id_t2_src)
call error_handler('error in inquring variable id of average_T2', rcode)
rcode = nf_get_var_double(ncid_src,id_t2_src, t2_avg )
call error_handler('error in getting average_T2 data', rcode)
rcode = nf_inq_varid(ncid_src, 'average_DT', id_dt_src)
call error_handler('error in inquring variable id of average_Dt', rcode)
rcode = nf_get_var_double(ncid_src,id_dt_src, dt_avg )
call error_handler('error in getting average_DT data', rcode)
endif
!--- read the global attributes from source data and write it to output file
rcode = nf_inq_natts(ncid_src, natt )
call error_handler('error in inquring natts of file '//trim(src_data), rcode)
do i = 1, natt
rcode = nf_inq_attname(ncid_src,NF_GLOBAL,i,name )
call error_handler('error in inquring attribute name of file '//trim(src_data), rcode)
rcode = nf_inq_att(ncid_src,NF_GLOBAL,name,type,len)
call error_handler('error in inquring attribute of file '//trim(src_data), rcode)
if(type == NF_CHAR ) then ! only keep character global attributes
if(len .le. 512) then
if(trim(name) == 'filename') then
catt = trim(dst_data)
else
rcode = nf_get_att_text(ncid_src,NF_GLOBAL,name,catt)
call error_handler('error in getting attribute text of file '//trim(src_data), rcode)
endif
if(is_root_pe) then
rcode = nf_put_att_text(ncid_dst, NF_GLOBAL, name, len, catt )
call error_handler('error in putting attribute to file '//trim(dst_data), rcode)
endif
else
write(stdout(),*)'GLOBAL ATT '//trim(name)//' too long - not reading this metadata'
endif
endif
enddo
!--- read dimension information from src_data and write them to dst_data file.
if( .not. is_root_pe) then
rcode = nf_close(ncid_src)
return
endif
rcode = nf_inq_ndims(ncid_src, ndims)
call error_handler('error in inquring ndims of file '//trim(src_data), rcode)
do j = 1, ndims_dst
if(lowercase(trim(dims_name(j)))=='time' ) then
rcode = nf_def_dim(ncid_dst, dims_name(j), NF_UNLIMITED, dims_dst(j))
call error_handler('error in defining time dimension of file '//trim(dst_data), rcode)
rcode = nf_def_var(ncid_dst, dims_name(j), NF_DOUBLE, 1, dims_dst(j), id_axes(j))
call error_handler('error in defining time variable of file '//trim(dst_data), rcode)
id_time_dst = id_axes(j)
dim_time_dst = dims_dst(j)
else
select case(dims_cart(j))
case('X')
len = Dst%ni
case('Y')
len = Dst%nj
case('Z')
len = Dst%nk
end select
rcode = nf_def_dim(ncid_dst, dims_name(j), len, dims_dst(j))
call error_handler('error in defining dimension '//trim(dims_name(j))//' of file '//trim(dst_data), rcode)
rcode = nf_def_var(ncid_dst, dims_name(j), NF_DOUBLE, 1, dims_dst(j), id_axes(j) )
call error_handler('error in defining axis var '//trim(dims_name(j))//' of file '//trim(dst_data), rcode)
endif
!--- write out axis attribute
rcode = nf_inq_varnatts(ncid_src, dims_id(j), natt_dim)
call error_handler('error in inquring number of attributes', rcode)
do i = 1, natt_dim
rcode = nf_inq_attname(ncid_src, dims_id(j), i, name)
rcode = nf_copy_att(ncid_src,dims_id(j), name, ncid_dst, id_axes(j))
call error_handler('error in copy attribute '//trim(name), rcode)
enddo
enddo
!--- read attribue of each field and write them to dst_data file
do n = 1, num_flds
rcode = nf_inq_varid(ncid_src,fld_name(n), id_fld)
call error_handler('error in inquring id of field '//trim(fld_name(n))// ' of file '//trim(src_data), rcode)
rcode = nf_inq_varndims( ncid_src, id_fld, ndims )
call error_handler('error in inquring ndims of field '//trim(fld_name(n))// ' of file '//trim(src_data), rcode)
rcode = nf_inq_vardimid( ncid_src, id_fld, dimids )
call error_handler('error in inquring dimids of field '//trim(fld_name(n))// ' of file '//trim(src_data), rcode)
do i = 1, ndims
rcode = nf_inq_dimname ( ncid_src, dimids(i), name )
call error_handler('error in inquring dimension name of dimid', rcode)
do j = 1, ndims_dst
if(trim(name) == trim(dims_name(j))) then
fld_dims(i) = dims_dst(j)
exit
endif
enddo
enddo
rcode = nf_def_var(ncid_dst, fld_name(n), NF_DOUBLE, ndims, fld_dims(1:ndims), id_dst_fld(n) )
!--- copy the source field attribute to destination field.
rcode = nf_inq_varnatts(ncid_src, id_fld, natt )
call error_handler('error in inquring natt of field '//trim(fld_name(n))// ' of file '//trim(src_data), rcode)
do i = 1, natt
rcode = nf_inq_attname(ncid_src, id_fld, 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_fld, trim(name), ncid_dst, id_dst_fld(n))
call error_handler('error in copy attribute '//trim(name)//' of variable '//trim(fld_name(n)), rcode)
endif
enddo
if(.not. has_missing(n)) then
rcode = nf_put_att_double(ncid_dst, id_dst_fld(n), 'missing_value', NF_DOUBLE, 1, missing_value(n:n) )
call error_handler('error in put missing_value attribute of var '//trim(fld_name(n))// &
' of file '//trim(dst_data), rcode)
endif
call error_handler('error in defining var '//trim(fld_name(n))//' of file '//trim(dst_data), rcode)
enddo
!--- define time avg info, suppose the time_avg_info is specified by average_T1, average_T2, average_DT.
if(tavg_exist) then
rcode = nf_def_var(ncid_dst, 'average_T1', NF_DOUBLE, 1, dim_time_dst, id_t1_dst)
rcode = nf_inq_varnatts(ncid_src, id_t1_src, natt )
call error_handler('error in inquring natt of average_T1 of file '//trim(src_data), rcode)
do i = 1, natt
rcode = nf_inq_attname(ncid_src, id_t1_src, i, name)
call error_handler('error in inquring average_T1 attribute '//trim(name), rcode)
rcode = nf_copy_att(ncid_src,id_t1_src, name, ncid_dst, id_t1_dst)
call error_handler('error in copy average_T1 attribute '//trim(name), rcode)
enddo
rcode = nf_def_var(ncid_dst, 'average_T2', NF_DOUBLE, 1, dim_time_dst, id_t2_dst)
rcode = nf_inq_varnatts(ncid_src, id_t2_src, natt )
call error_handler('error in inquring natt of average_T2 of file '//trim(src_data), rcode)
do i = 1, natt
rcode = nf_inq_attname(ncid_src, id_t2_src, i, name)
call error_handler('error in inquring average_T2 attribute '//trim(name), rcode)
rcode = nf_copy_att(ncid_src,id_t2_src, name, ncid_dst, id_t2_dst)
call error_handler('error in copy average_T2 attribute '//trim(name), rcode)
enddo
rcode = nf_def_var(ncid_dst, 'average_DT', NF_DOUBLE, 1, dim_time_dst, id_dt_dst)
rcode = nf_inq_varnatts(ncid_src, id_dt_src, natt )
call error_handler('error in inquring natt of average_DT of file '//trim(src_data), rcode)
do i = 1, natt
rcode = nf_inq_attname(ncid_src, id_dt_src, i, name)
call error_handler('error in inquring average_DT attribute '//trim(name), rcode)
rcode = nf_copy_att(ncid_src,id_dt_src, name, ncid_dst, id_dt_dst)
call error_handler('error in copy average_DT attribute '//trim(name), rcode)
enddo
endif
rcode = nf_enddef(ncid_dst)
rcode = nf_close(ncid_src)
!--- write axis data to dst_data
start = 1; nwrite = 1
do j = 1, ndims_dst
select case(dims_cart(j))
case('X')
nwrite(1) = Dst%ni
if(dims_pos(j) == 'T' .or. dims_pos(j) == 'N') then
rcode = nf_put_vara_double(ncid_dst,id_axes(j), start, nwrite, Dst%grid_xt )
else
rcode = nf_put_vara_double(ncid_dst,id_axes(j), start, nwrite, Dst%grid_xc )
endif
call error_handler('error in putting data of '//trim(dims_name(j))//' to file '//trim(dst_data), rcode)
case('Y')
nwrite(1) = Dst%nj
if(dims_pos(j) == 'T' .or. dims_pos(j) == 'E') then
rcode = nf_put_vara_double(ncid_dst,id_axes(j), start, nwrite, Dst%grid_yt )
else
rcode = nf_put_vara_double(ncid_dst,id_axes(j), start, nwrite, Dst%grid_yc )
endif
call error_handler('error in putting data of '//trim(dims_name(j))//' to file '//trim(dst_data), rcode)
case('Z')
nwrite(1) = Dst%nk
rcode = nf_put_vara_double(ncid_dst,id_axes(j), start, nwrite, Dst%zt )
call error_handler('error in putting data of '//trim(dims_name(j))//' to file '//trim(dst_data), rcode)
end select
enddo
end subroutine setup_meta
!#####################################################################
subroutine process_data
integer :: ncid_src, rcode, i, j, n, m, k, l, nf, nk_src, nk_dst
integer :: start(4), nread(4), nwrite(4), id_fld(num_flds)
integer :: isc, iec, jsc, jec
real :: temp1, temp2
real, dimension(:,:,:,:), allocatable :: data_src, data_dst
real, dimension(:,:,:), allocatable :: tmp1, tmp2, tmp_src
real, dimension(:,:,:), allocatable :: depth_src_3d, depth_dst_3d
real, dimension(:,:), allocatable :: mask_in
type(horiz_interp_type), pointer :: Interp => NULL()
real, dimension(:,:), pointer :: sinrot => NULL(), cosrot => NULL()
real, dimension(:,:,:), pointer :: mask => NULL()
logical :: is_first = .true.
isc = Dst%isc; iec = Dst%iec; jsc = Dst%jsc; jec = Dst%jec;
!-------------------------------------------------------------------
!--- read src data, remap onto current grid and write output -------
rcode = nf_open(trim(src_data), NF_NOWRITE, ncid_src)
do n = 1, num_flds
rcode = nf_inq_varid(ncid_src, fld_name(n), id_fld(n) )
call error_handler('error in inquring varid of '//trim(fld_name(n))//' of file '//trim(src_data), rcode)
enddo
write(stdout(),*)'***************************************'
write(stdout(),*)' ntime_src is ', ntime_src
if(do_vertical_interp) then
allocate(depth_src_3d(isc:iec,jsc:jec,Src%nk), depth_dst_3d(isc:iec,jsc:jec,Dst%nk))
do k=1, Src%nk
depth_src_3d(:,:,k) = Src%zt(k)
enddo
do k=1, Dst%nk
depth_dst_3d(:,:,k) = Dst%zt(k)
enddo
endif
allocate(mask_in(Src%ni, Src%nj) )
do m = 1, ntime_src
write(stdout(),*)'******************* at time level', m
n = 1
!--- write out time information
if(is_root_pe) then
if(has_taxis) then
rcode = nf_put_var1_double(ncid_dst, id_time_dst, m, time_value(m))
call error_handler('error in put time data to file '//trim(src_data), rcode )
if(tavg_exist) then
rcode = nf_put_var1_double(ncid_dst, id_t1_dst, m, t1_avg(m))
call error_handler('error in put average_T1 data to file '//trim(src_data), rcode )
rcode = nf_put_var1_double(ncid_dst, id_t2_dst, m, t2_avg(m))
call error_handler('error in put average_T2 data to file '//trim(src_data), rcode )
rcode = nf_put_var1_double(ncid_dst, id_dt_dst, m, dt_avg(m))
call error_handler('error in put average_DT data to file '//trim(src_data), rcode )
endif
endif
endif
do while ( n .le. num_flds )
nf = 1
if(vector_fld(n) .and. Src%is_tripolar) nf = 2
nk_src = 1
if(has_zaxis(n)) nk_src = Src%nk
nk_dst = 1
if(has_zaxis(n)) nk_dst = Dst%nk
write(stdout(),*)'************************************************************'
write(stdout(),*)'The following are for the # ',n,' field '//trim(fld_name(n))
write(stdout(),*)'has_zaxis is ', has_zaxis(n)
write(stdout(),*)'nk_dst is ', nk_dst
allocate(data_src(Src%ni, Src%nj, nk_src,nf), data_dst(Dst%ni, Dst%nj, nk_dst,nf) )
allocate(tmp_src(Src%ni, Src%nj, nk_src) )
allocate(tmp1(isc:iec,jsc:jec,nk_src), tmp2(isc:iec,jsc:jec,nk_dst) )
select case(fld_pos(n))
case('T')
Interp => Interp_T
case('C')
Interp => Interp_C
case('E')
Interp => Interp_E
case('N')
Interp => Interp_N
end select
start = 1; nread = 1
nread(1) = Src%ni; nread(2) = Src%nj;
if(has_zaxis(n)) then
nread(3) = nk_src; start(4) = m
else
start(3) = m
endif
rcode = nf_get_vara_double(ncid_src,id_fld(n), start, nread, data_src(:,:,:,1) )
call error_handler('error in getting '//trim(fld_name(n))//' from file '//trim(src_data), rcode )
!--- if the field is a vector field and src_grid is tripolar, read the next field and
!--- rotate data onto spherical grid.
if( vector_fld(n) .and. Src%is_tripolar ) then
select case(fld_pos(n))
case('T')
sinrot => Src%T%sinrot; cosrot => Src%T%cosrot
case('C')
sinrot => Src%C%sinrot; cosrot => Src%C%cosrot
case('E')
sinrot => Src%E%sinrot; cosrot => Src%E%cosrot
case('N')
sinrot => Src%N%sinrot; cosrot => Src%N%cosrot
end select
rcode = nf_get_vara_double(ncid_src,id_fld(n+1), start, nread, data_src(:,:,:,2) )
do k = 1, nk_src
do j = 1, Src%nj
do i = 1, Src%ni
if(data_src(i,j,k,1) /= missing_value(n)) then
temp1 = data_src(i,j,k,1)*cosrot(i,j)-data_src(i,j,k,2)*sinrot(i,j)
temp2 = data_src(i,j,k,1)*sinrot(i,j)+data_src(i,j,k,2)*cosrot(i,j)
data_src(i,j,k,1) = temp1
data_src(i,j,k,2) = temp2
endif
enddo
enddo
enddo
endif
!--- extrap data onto land grid
do l = 1, nf
if(use_source_vertical_grid) then !--- the source grid mask will decide the destination data mask.
if(apply_mask) then
do k = 1, size(data_src,3)
where(data_src(:,:,k,l) == missing_value(n) )
mask_in = 0.0
elsewhere
mask_in = 1.0
end where
call horiz_interp(Interp, data_src(:,:,k,l), tmp2(:,:,k), mask_in=mask_in, &
missing_value=missing_value(n) )
enddo
else
if(any(data_src(:,:,:,l) == missing_value(n))) then
if(Src%is_tripolar .and. is_first) then
is_first = .false.
write(stdout(),*)'WARNING: when the source is tripolar, ', &
'the laplacian extrapolation will not very accurate at the tripolar region. ', &
'Even the source grid is lat-lon grid, the extrapolation will not also not ', &
'very accurate if the lat-lon grid is not uniform. '
endif
call extrap(data_src(:,:,:,l), tmp_src, stop_crit(n), missing_value(n), fld_pos(n) )
!--- remap data onto destination horizontal grid
do k = 1, size(tmp_src,3)
call horiz_interp(Interp, tmp_src(:,:,k), tmp2(:,:,k))
enddo
else
!--- remap data onto destination horizontal grid
do k = 1, size(data_src,3)
call horiz_interp(Interp, data_src(:,:,k,l), tmp2(:,:,k))
enddo
endif
endif
else
if(any(data_src(:,:,:,l) == missing_value(n))) then
if(Src%is_tripolar .and. is_first) then
is_first = .false.
write(stdout(),*)'WARNING: when the source is tripolar, ', &
'the laplacian extrapolation will not very accurate at the tripolar region. ', &
'Even the source grid is lat-lon grid, the extrapolation will not also not ', &
'very accurate if the lat-lon grid is not uniform. '
endif
call extrap(data_src(:,:,:,l), tmp_src, stop_crit(n), missing_value(n), fld_pos(n) )
!--- remap data onto destination horizontal grid
do k = 1, size(tmp_src,3)
call horiz_interp(Interp, tmp_src(:,:,k), tmp1(:,:,k))
enddo
else
!--- remap data onto destination horizontal grid
do k = 1, size(data_src,3)
call horiz_interp(Interp, data_src(:,:,k,l), tmp1(:,:,k))
enddo
endif
if(has_zaxis(n) .and. do_vertical_interp ) then
call interp_1d(depth_src_3d,depth_dst_3d,tmp1,tmp2)
else
tmp2 = tmp1
endif
!--- apply mask if needed,
if(apply_mask) then
select case(fld_pos(n))
case('T')
mask => Dst%T%mask
case('C')
mask => Dst%C%mask
case('E')
mask => Dst%E%mask
case('N')
mask => Dst%N%mask
end select
do k = 1, nk_dst
do j = jsc, jec
do i = isc, iec
if(mask(i,j,k) < 0.5) tmp2(i,j,k) = missing_value(n)
enddo
enddo
enddo
endif
endif
!--- get global data ------------------------------
call mpp_domains_set_stack_size(2*Dst%ni*Dst%nj*nk_dst)
call mpp_global_field(Dst%domain,tmp2, data_dst(:,:,:,l) )
enddo
!--- rotate the data if the field is a vector field and dst_grid is tripolar.
if( vector_fld(n) .and. Dst%is_tripolar ) then
select case(fld_pos(n))
case('T')
sinrot => Dst%T%sinrot; cosrot => Dst%T%cosrot
case('C')
sinrot => Dst%C%sinrot; cosrot => Dst%C%cosrot
case('E')
sinrot => Dst%E%sinrot; cosrot => Dst%E%cosrot
case('N')
sinrot => Dst%C%sinrot; cosrot => Dst%C%cosrot
end select
do k = 1, nk_dst
do j = 1, Dst%nj
do i = 1, Dst%ni
if(data_dst(i,j,k,1) /= missing_value(n)) then
temp1 = data_dst(i,j,k,1)*cosrot(i,j)+data_dst(i,j,k,2)*sinrot(i,j)
temp2 = -data_dst(i,j,k,1)*sinrot(i,j)+data_dst(i,j,k,2)*cosrot(i,j)
data_dst(i,j,k,1) = temp1
data_dst(i,j,k,2) = temp2
endif
enddo
enddo
enddo
select case(fld_pos(n))
case('C')
do k = 1, nk_dst
do i=1,Dst%ni/2
if(data_dst(Dst%ni-i,Dst%nj,k,1) /= missing_value(n)) then
data_dst(i,Dst%nj,k,1) = -1.*data_dst(Dst%ni-i,Dst%nj,k,1)
data_dst(i,Dst%nj,k,2) = -1.*data_dst(Dst%ni-i,Dst%nj,k,2)
endif
enddo
enddo
case('N')
do k = 1, nk_dst
do i=1,Dst%ni/2
if(data_dst(Dst%ni-i+1,Dst%nj,k,1) /= missing_value(n)) then
data_dst(i,Dst%nj,k,1) = -1.*data_dst(Dst%ni-i+1,Dst%nj,k,1)
data_dst(i,Dst%nj,k,2) = -1.*data_dst(Dst%ni-i+1,Dst%nj,k,2)
endif
enddo
enddo
end select
endif
!--- write out data from root_pe
if(is_root_pe) then
if(debug) write(stdout(),*)'NOTE: after regrid data chksum :', mpp_chksum(data_dst, (/mpp_root_pe()/) )
do l = 1, nf
start = 1; nwrite = 1
nwrite(1) = Dst%ni; nwrite(2) = Dst%nj;
if(has_zaxis(n)) then
nwrite(3) = nk_dst
start(4) = m
rcode = nf_put_vara_double(ncid_dst, id_dst_fld(n+l-1), start, nwrite, data_dst(:,:,:,l) )
call error_handler('error in putting '//trim(fld_name(n+l-1))//' in file '//trim(dst_data), rcode )
else
start(3) = m
rcode = nf_put_vara_double(ncid_dst, id_dst_fld(n+l-1), start, nwrite, data_dst(:,:,1,l) )
call error_handler('error in putting '//trim(fld_name(n+l-1))//' in file '//trim(dst_data), rcode )
endif
! rcode = nf_put_vara_double(ncid_dst, id_dst_fld(n+l-1), start, nwrite, data_dst(:,:,:,l) )
! call error_handler('error in putting '//trim(fld_name(n+l-1))//' in file '//trim(dst_data), rcode )
enddo
endif
n = n+nf
deallocate(data_src, data_dst, tmp_src, tmp1, tmp2)
enddo
enddo
rcode = nf_close(ncid_src)
if(is_root_pe) rcode = nf_close(ncid_dst)
if(do_vertical_interp) deallocate(depth_src_3d, depth_dst_3d)
deallocate(mask_in)
end subroutine process_data
!#####################################################################
subroutine regrid_end
call fms_io_exit
call fms_end()
end subroutine regrid_end
!#####################################################################
subroutine extrap(data_in, data_out, crit, missing, pos)
real, dimension(:,:,:), intent(in) :: data_in
real, dimension(:,:,:), intent(out) :: data_out
real, intent(in) :: crit, missing
character(len=1), intent(in) :: pos
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
ni = size(data_in,1)
nj = size(data_in,2)
nk = size(data_in,3)
tmp = 0.0
do j= 1, nj
do i= 1, ni
if(data_in(i,j,1) == missing) then
tmp(i,j) = initial_guess
endif
enddo
enddo
do k=1,nk
do j= 1, nj
do i= 1, ni
if(data_in(i,j,k) == missing) 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,pos)
! iterate
n=1
do
resmax=0.0
do j= 1, nj
do i= 1, ni
res(i,j) = 0.25*(tmp(i-1,j)+tmp(i+1,j)+tmp(i,j-1)+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,pos)
n=n+1
enddo
write(stdout(),'(a,i4,a,i4)') 'At k-level= ',k
write(stdout(),'(a,i6,a)') 'Stopped after ',n,' iterations'
write(stdout(),'(a,f10.4)') 'maxres= ',resmax
enddo
end subroutine extrap
!##########################################################
subroutine fill_boundaries(data, pos)
real, dimension(0:,0:), intent(inout) :: data
character(len=1), intent(in) :: pos
integer :: i,j, ni, nj
ni = size(data,1) - 2
nj = size(data,2) - 2
if(Dst%is_cyclic) then
data(0,1:nj) = data(ni,1:nj)
data(ni+1,1:nj) = data(1,1:nj)
endif
if(Dst%is_tripolar) then
do i = 1, ni
select case(pos)
case('T')
data(i,nj+1) = data(ni-i+1,nj)
case('C')
data(i,nj+1) = data(ni-i,nj-1)
case('E')
data(i,nj+1) = data(ni-i+1,nj-1)
case('N')
data(i,nj+1) = data(ni-i,nj)
end select
enddo
endif
return
end subroutine fill_boundaries
!#####################################################################
! 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
call mpp_error(FATAL,trim(msg) )
end subroutine error_handler
!#####################################################################
end program regrid