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 20.0 2013/12/14 00:31:05 fms Exp $'
character(len=128) :: tagname = '$Name: tikal $'
!--- 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)
integer ::fsize = 65536, inital = 0
!--- open the dst_data file from root_pe.
if(is_root_pe) then
#ifdef use_netCDF3
rcode = NF__CREATE(trim(dst_data), NF_CLOBBER, inital, fsize, ncid_dst)
#else ! netcdf 4
rcode = NF__CREATE(trim(dst_data), IOR(NF_NETCDF4,NF_CLASSIC_MODEL), inital, fsize, ncid_dst )
#endif
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
integer :: kstart, kend
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
! for vertical interpolation, Set value of levels shallower than the shallowest source level to
! be the source value at shallowest level.
! Set value of levels deeper than the deepest source level to
! be the source value at deepest level.
do kstart = 1, Dst%nk
if( Dst%zt(kstart) .GE. Src%zt(1) ) exit
enddo
do kend = Dst%nk, 1, -1
if( Dst%zt(kend) .LE. Src%zt(Src%nk) ) exit
enddo
if( kstart > 1 ) then
write(stdout(),*)"NOTE from regrid_3d: the value from level 1 to level ", kstart-1, &
" will be set to the value at the shallowest source levle."
endif
if( kend < Dst%nk ) then
write(stdout(),*)"NOTE from regrid_3d: the value from level ", kend+1, " to level ", Dst%nk, &
" will be set to the value at the deepest source levle."
endif
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
do k = 1, kstart-1
tmp2(:,:,k) = tmp1(:,:,1)
enddo
do k = kend+1, nk_dst
tmp2(:,:,k) = tmp1(:,:,nk_dst)
enddo
call interp_1d(depth_src_3d,depth_dst_3d(:,:,kstart:kend),tmp1,tmp2(:,:,kstart:kend))
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