!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! 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
!
!
! Since most analyses are placed on a spherical latitude-longitude grid and
! most global ocean models configured from mom4 are run with tripolar grids.
! A regridding tool is needed to regrid data from tripolar grid onto latitude-longitude
! grid. This program map data from any logical rectangular grid
! (tripolar or latitude-longitude) to a spherical latitude-longitude grid.
! When the data is on Tracer cell (T-cell), the interpolation will be a conservative
! scheme. When the data is on other position (C-cell, E-cell or N-cell), regridding
! is accomplished non-conservatively using a nearest neighbor distance weighting algorithm.
!
!
! Before using this regridding tool, you need to use preprocessing tool
! make_xgrids to generate the grid file which contains source grid, destination
! grid and exchange grid information between source grid and destination grid.
! These exchange grid information is needed when doing conservative remapping.
! Suppose the file name of your source grid is src_grid.nc and the file name of
! your destination grid is dst_grid.nc (dst grid should be spherical
! latitude-lontitude grid), the command will be
! "make_xgrids -o src_grid.nc -a dst_grid.nc".
! Before using make_xgrids, you need to make sure src_grid.nc does not contain
! any exchange grid information. Otherwise you are not going to get the
! desired exchange grid information. If your src_grid.nc do contains exchange
! grid information, you can remove all those fields by using ncks. Those fields are
! AREA_ATMxOCN, DI_ATMxOCN, DJ_ATMxOCN, I_ATM_ATMxOCN, J_ATM_ATMxOCN, I_OCN_ATMxOCN,
! J_OCN_ATMxOCN, AREA_ATMxLND, DI_ATMxLND, DJ_ATMxLND, I_ATM_ATMxLND, J_ATM_ATMxLND,
! I_LND_ATMxLND, J_LND_ATMxLND, AREA_LNDxOCN, DI_LNDxOCN, DJ_LNDxOCN,
! I_LND_LNDxOCN, J_LND_LNDxOCN, I_OCN_LNDxOCN, J_OCN_LNDxOCN, xba, yba, xta, yta,
! AREA_ATM, xbl, ybl, xtl, ytl, AREA_LND, AREA_LND_CELL, xto, yto and AREA_OCN.
! After the grid file is generated, the grid file should be passed into the program
! through nml option "grid_spec_file".
!
! 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 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. The data will be always on the source vertical grid.
! Previous experiences show that linear vertical interpolation will create lots of noises.
! If we find better vertical interpolation algorithm, we may implement it in the
! future.
!
use mpp_mod, only : mpp_npes, mpp_pe, mpp_root_pe, mpp_error, FATAL, NOTE, WARNING
use mpp_mod, only : stdout, stdlog, mpp_chksum
use mpp_domains_mod, only : mpp_define_layout, mpp_define_domains, domain2d
use mpp_domains_mod, only : mpp_global_field, mpp_get_compute_domain
use mpp_domains_mod, only : mpp_domains_set_stack_size
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, epsln
use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_type
implicit none
#include "netcdf.inc"
integer, parameter :: max_flds = 20
integer, parameter :: max_dims = 20
integer, parameter :: max_atts = 20
real, parameter :: small = 1.0e-6
!--- 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 source and target grid information,
! exchange grid information between source and target grid. This grid file can
! be created using preprocessing tool make_xgrids.
!
!
! Name of field to be regridded 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.
!
!
! Number of nearest neighbors for regridding.
!
!
! Maximum region of influence around destination grid points.
!
!
! 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 = 'INPUT/src_data.nc'
character(len=128) :: grid_spec_file = 'INPUT/grid_spec.nc'
character(len=128) :: dst_data = 'OUTPUT/dst_data.nc'
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.
integer :: num_nbrs = 5
real :: max_dist = 0.05
logical :: debug = .false.
namelist /regrid_nml/ num_flds, src_data, dst_data, grid_spec_file, fld_name, &
fld_pos, vector_fld, num_nbrs, max_dist, debug
!--- version information ---------------------------------------------
character(len=128) :: version = '$Id: regrid.F90,v 14.0 2007/03/15 22:45:28 fms Exp $'
character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $'
!--- variables about source data or on source grid
integer :: ncid_src ! ncid corresponding to source file.
integer :: ni_src, nj_src, nk_src ! source grid size
integer :: ntime_src ! number of time levels of src data
logical :: src_is_cyclic ! indicate if the source is cyclic
logical :: src_is_tripolar ! indicate if the source grid is tripolar
real, dimension(max_flds) :: missing_value = -1e20 ! missing value for source field.
logical :: has_taxis,tavg_exist ! indicate if all the field has time axis
logical, dimension(max_flds) :: has_zaxis = .false. ! indicate if the field has z axis
real, dimension(:), allocatable :: time_value ! time axis value
real, dimension(:), allocatable :: t1_avg, t2_avg, dt_avg ! time average data
real, dimension(:,:), allocatable :: area_ocn ! ocean area
real, dimension(:), allocatable :: zt_src ! source vertical grid
real, dimension(:,:), allocatable :: geolon_t, geolat_t ! geographic grid location on T-cell
real, dimension(:,:), allocatable :: geolon_e, geolat_e ! geographic grid location on E-cell
real, dimension(:,:), allocatable :: geolon_c, geolat_c ! geographic grid location on C-cell
real, dimension(:,:), allocatable :: geolon_n, geolat_n ! geographic grid location on N-cell
real, dimension(:,:), allocatable :: sinrot, cosrot ! rotation
!--- variables about destination data
integer :: ni_dst, nj_dst ! destination grid size
integer :: ncid_dst ! ncid corresponding to destination file.
integer :: id_t1_dst, id_t2_dst, id_dt_dst
integer :: id_time_dst, id_dst_fld(max_flds)
character(len=128) :: dims_name(max_dims) ! output file dimension name
integer :: dims_id(max_dims)
character(len=1) :: dims_cart(max_dims) ! 'X', 'Y', 'Z' or 'T'
character(len=1) :: dims_pos(max_dims) ! 'T', 'C', 'E' or 'N'
real, dimension(:), allocatable :: xt_dst, yt_dst, xb_dst, yb_dst !lat-lon destination grid
type(domain2D), save :: Domain ! destination grid domain.
integer :: isc, iec, jsc, jec ! compute domain decompsition of dst grid.
integer :: ndims_dst ! number of dimensions in the output file
!--- exchange grid information
integer, dimension(:), allocatable :: i_atm_atmxocn
integer, dimension(:), allocatable :: j_atm_atmxocn
integer, dimension(:), allocatable :: i_ocn_atmxocn
integer, dimension(:), allocatable :: j_ocn_atmxocn
real, dimension(:), allocatable :: area_atmxocn
!--- other variables
logical :: is_root_pe ! if current pe is root pe.
type(horiz_interp_type), target :: Interp_c, Interp_e, Interp_n ! horiz_interp data at C,E,N-cell
integer :: jstart ! starting index to do interpolation
!--- begin of the program -------------------------------------------
call regrid_init()
call read_src_grid()
call read_dst_grid()
call setup_interp()
call read_xgrid()
call get_src_info()
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
!--- when all the data are on T-cell, should always use 1 processor. Since no parallization is
!--- implemeneted.
if(all(fld_pos(1:num_flds) == 'T') .and. mpp_npes() > 1) then
write(stdout(),*)'WARNING: All the fields are on T-cell but npes is greater than 1, '//&
'should change processor count to 1'
endif
!--- 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
if(fld_pos(n) .ne. 'C') call error_handler('regrid:When the field is a vector field, we assume '// &
'this field is located at C-cell center. ')
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 grid spec file that contains exchange grid is: ', trim(grid_spec_file)
Write(stdout(),*) 'The output data will be located on source vertical grid, '
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
is_root_pe = mpp_pe() == mpp_root_pe()
end subroutine regrid_init
!#################################################################
!--- read the src grid information from file grid_spec_file. The grid information
!--- will include boundary condition, horizontal and vertical grid information.
!--- No need to read source grid land/sea mask, since we can always obtain source grid land/sea
!--- mask from source data (the value of the point is equal to missing value will be
!--- land point).
subroutine read_src_grid
integer :: i
integer :: rcode, ncid, natt
character(len=64) :: name, catt
real, allocatable :: tmp(:,:)
logical :: is_new_grid
real :: D2R
D2R = PI/180.
rcode = nf_open(trim(grid_spec_file), NF_NOWRITE, ncid)
call error_handler('error in open file '//trim(grid_spec_file), rcode)
!--- get boundary condition ( cyclic or tripolar )
!--- the boundary condition are specified by global attribute
!--- 'x_boundary_type' and 'y_boundary_type'
rcode = nf_inq_natts(ncid, natt)
call error_handler('error in inquring ngatts of file '//trim(grid_spec_file), 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(grid_spec_file), 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(grid_spec_file), rcode )
if(trim(catt) == 'cyclic') then
src_is_cyclic = .true.
Write(stdout(),*)' NOTE: x_boundary_type of grid '//trim(grid_spec_file)//' is cyclic'
src_is_cyclic = .false.
else
Write(stdout(),*)' NOTE: x_boundary_type of '//trim(grid_spec_file)//' 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(grid_spec_file), rcode )
if (trim(catt) == 'fold_north_edge') then
src_is_tripolar = .true.
Write(stdout(),*)' NOTE: y_boundary_type of '//trim(grid_spec_file)//' is tripolar'
else
src_is_tripolar = .false.
Write(stdout(),*)' NOTE: y_boundary_type of '//trim(grid_spec_file)//' is solid_walls'
endif
end select
end do
!--- To determine the grid is in new grid format or old grid format
is_new_grid = check_is_new_grid(ncid)
!--- get the vertical grid information
nk_src = get_dimlen(ncid, 'zt')
allocate(zt_src(nk_src) )
call get_var_real_1d(ncid, 'zt', zt_src )
if(is_new_grid) then
ni_src = get_dimlen(ncid, 'grid_x_T')
nj_src = get_dimlen(ncid, 'grid_y_T')
else
ni_src = get_dimlen(ncid, 'gridlon_t')
nj_src = get_dimlen(ncid, 'gridlat_t')
endif
allocate(geolon_t(ni_src,nj_src), geolat_t(ni_src,nj_src) )
allocate(geolon_c(ni_src,nj_src), geolat_c(ni_src,nj_src) )
allocate(geolon_e(ni_src,nj_src), geolat_e(ni_src,nj_src) )
allocate(geolon_n(ni_src,nj_src), geolat_n(ni_src,nj_src) )
allocate(sinrot (ni_src,nj_src), cosrot (ni_src,nj_src) )
allocate(tmp (ni_src,nj_src) )
if(is_new_grid) then
call get_var_real_2d(ncid, 'x_T', geolon_t)
call get_var_real_2d(ncid, 'y_T', geolat_t)
call get_var_real_2d(ncid, 'x_C', geolon_c)
call get_var_real_2d(ncid, 'y_C', geolat_c)
call get_var_real_2d(ncid, 'x_T', geolon_e)
call get_var_real_2d(ncid, 'y_T', geolat_e)
call get_var_real_2d(ncid, 'x_C', geolon_n)
call get_var_real_2d(ncid, 'y_C', geolat_n)
call get_var_real_2d(ncid, 'angle_C', tmp)
sinrot = sin(tmp*D2R)
cosrot = cos(tmp*D2R)
else
call get_var_real_2d(ncid, 'geolon_t', geolon_t)
call get_var_real_2d(ncid, 'geolat_t', geolat_t)
call get_var_real_2d(ncid, 'geolon_c', geolon_c)
call get_var_real_2d(ncid, 'geolat_c', geolat_c)
call get_var_real_2d(ncid, 'geolon_e', geolon_e)
call get_var_real_2d(ncid, 'geolat_e', geolat_e)
call get_var_real_2d(ncid, 'geolon_n', geolon_n)
call get_var_real_2d(ncid, 'geolat_n', geolat_n)
call get_var_real_2d(ncid, 'sin_rot', sinrot )
call get_var_real_2d(ncid, 'cos_rot', cosrot )
endif
rcode = nf_close(ncid)
return
end subroutine read_src_grid
!#################################################################
!--- read the destination horizontal grid information from file "grid_spec_file".
subroutine read_dst_grid
integer :: ncid, rcode
!--- first reading grid information from file "grid_spec_file"
rcode = nf_open(trim(grid_spec_file), NF_NOWRITE, ncid)
call error_handler('error in open file '//trim(grid_spec_file), rcode)
!--- get the destination grid size
ni_dst = get_dimlen(ncid, 'xta')
nj_dst = get_dimlen(ncid, 'yta')
allocate(xt_dst(ni_dst), yt_dst(nj_dst), xb_dst(ni_dst+1), yb_dst(nj_dst+1) )
call get_var_real_1d(ncid, 'xta', xt_dst)
call get_var_real_1d(ncid, 'yta', yt_dst)
call get_var_real_1d(ncid, 'xba', xb_dst)
call get_var_real_1d(ncid, 'yba', yb_dst)
rcode = nf_close(ncid)
!--- For the conservative interpolation, it is very complicated to implmeent parallization
!--- also the conservative interpolation is efficiency enough for postprocessing purpose.
!--- That's the reason we didn't implement parallization for the conservative scheme.
return
end subroutine read_dst_grid
!#####################################################################
!--- In many situation, the source grid match the destination grid in the
!--- region southern of some latitude, indexed with jstart.
!--- Interpolation is necessary only for the region northen of jstart.
!--- this will save lots of interpolation time. When some field is located
!--- on C, E, or N-cell, call horiz_interp_new to set up the interpolation.
!--- in this case, parallization is implemented.
subroutine setup_interp
integer :: i, j, layout(2)
real :: D2R
D2R = PI/180.
!--- find the starting point to do interpolatioon jstart
if(ni_src == ni_dst) then
J_LOOP: do j = 1, nj_dst
do i = 1, ni_dst
if( abs( geolon_t(i,j) - xt_dst(i) ) .gt. small ) exit J_LOOP
if( abs( geolat_t(i,j) - yt_dst(j) ) .gt. small ) exit J_LOOP
if( abs( geolon_c(i,j) - xb_dst(i+1) ) .gt. small ) exit J_LOOP
if( abs( geolat_c(i,j) - yb_dst(j+1) ) .gt. small ) exit J_LOOP
enddo
enddo J_LOOP
jstart = j
else
jstart = 1
endif
if(jstart .ne. 1) write(stdout(),*)'NOTE: horizontal interpolation will be done ', &
'starting at lat(',jstart,')= ',yt_dst(jstart)
!--- set up domain for parallization if needed
layout = (/1,0/)
call mpp_define_layout((/1,ni_dst,jstart,nj_dst/), mpp_npes(), layout)
call mpp_define_domains((/1,ni_dst,jstart,nj_dst/),layout, Domain )
call mpp_get_compute_domain(Domain,isc,iec,jsc,jec)
!--- call horiz_interp_new to calculate the reampping weight in horizontal direction
if( any(fld_pos == 'C') ) then
call horiz_interp_new(Interp_c, geolon_c*D2R, geolat_c*D2R, xb_dst(isc+1:iec+1)*D2R, &
yb_dst(jsc+1:jec+1)*D2R, interp_method="spherical", &
num_nbrs=num_nbrs, max_dist=max_dist, src_modulo = src_is_cyclic)
endif
if( any(fld_pos == 'E') ) then
call horiz_interp_new(Interp_e, geolon_e*D2R, geolat_e*D2R, xb_dst(isc+1:iec+1)*D2R, &
yt_dst(jsc:jec)*D2R, interp_method="spherical", &
num_nbrs=num_nbrs, max_dist=max_dist, src_modulo = src_is_cyclic)
endif
if( any(fld_pos == 'N') ) then
call horiz_interp_new(Interp_n, geolon_n*D2R, geolat_n*D2R, xt_dst(isc:iec)*D2R, &
yb_dst(jsc+1:jec+1)*D2R, interp_method="spherical", &
num_nbrs=num_nbrs, max_dist=max_dist, src_modulo = src_is_cyclic)
endif
return
end subroutine setup_interp
!#################################################################
! read needed exchange grid information from file grid_spec_file
! for the purpose of first-order conservative interpolation for tracer data
! If needed, we may adding the option to do second-order conservative
! interpolation in the future. The algorithm will be very similiar to
! the algorithm used in coupler flux exchange, which is implemented in
! shared code xgrid.f90.
subroutine read_xgrid()
integer :: rcode, ncid, ni_lon, nj_lat, ncells, i, j
!--- It is needed only when doing conservative interpolation, which means
!--- there is a field located on T-cell.
if(.not. any(fld_pos == 'T')) return
!--- read xgrid from the grid_spec.nc file --------------------------
rcode = nf_open(grid_spec_file, NF_NOWRITE, ncid)
call error_handler('error in open file '//grid_spec_file, rcode)
!--- get exchange grid information
ncells = get_dimlen(ncid, 'I_ATM_ATMxOCN')
allocate(i_atm_atmxocn(ncells), j_atm_atmxocn(ncells), &
i_ocn_atmxocn(ncells), j_ocn_atmxocn(ncells), &
area_atmxocn(ncells) )
call get_var_int_1d (ncid, 'I_ATM_ATMxOCN', i_atm_atmxocn)
call get_var_int_1d (ncid, 'J_ATM_ATMxOCN', j_atm_atmxocn )
call get_var_int_1d (ncid, 'I_OCN_ATMxOCN', i_ocn_atmxocn )
call get_var_int_1d (ncid, 'J_OCN_ATMxOCN', j_ocn_atmxocn )
call get_var_real_1d(ncid, 'AREA_ATMxOCN', area_atmxocn )
rcode = nf_close(ncid)
!--- this variable will be used for the conservative interpolation.
allocate(area_ocn(ni_dst, jstart:nj_dst))
end subroutine read_xgrid
!#####################################################################
!--- This routine will read source data information from file src_data.
!--- This information includes missing_value for each field, if any field
!--- contains attribute "time_avg_info". Get the time information.
!--- Since the metadata of destination data will be determined by the
!--- source data, get the metadata information for the dst_data.
subroutine get_src_info
integer :: rcode, ndims, natt, natt_dim, id_fld, dimids(4), id_dim
integer :: num_has_taxis, i, j, m, n
character(len=128) :: name, cart, att_name
logical :: found_dim
rcode = nf_open(trim(src_data), NF_NOWRITE, ncid_src)
call error_handler('error in opening file '//trim(src_data), rcode)
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)
!--- get the missing value and check if the field has tag_exist attribute.
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
rcode = nf_get_att_double(ncid_src, id_fld, name, missing_value(n) )
endif
enddo
!--- 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)
if(lowercase(trim(name)) == 'time') then
cart = 'T'
if(num_has_taxis == 0) then
ntime_src = get_dimlen(ncid_src, trim(name) )
allocate(time_value(ntime_src) )
call get_var_real_1d(ncid_src, trim(name), time_value)
endif
num_has_taxis = num_has_taxis + 1
else
rcode = nf_inq_varnatts(ncid_src, id_dim, natt_dim)
call error_handler('error in inquring number of attributes of variable '//trim(name), rcode)
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
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
!--- some files don't have 'average_T1', 'average_T2' and 'average_DT' even though
!--- some variable do have time_avg_info attribute. We suppose 'average_T1',
!--- 'average_T2' and 'average_DT' always co-exist.
rcode = nf_inq_varid(ncid_src, 'average_T1', id_fld)
if(rcode == 0) then
allocate( t1_avg(ntime_src), t2_avg(ntime_src), dt_avg(ntime_src) )
call get_var_real_1d(ncid_src, 'average_T1', t1_avg )
call get_var_real_1d(ncid_src, 'average_T2', t2_avg )
call get_var_real_1d(ncid_src, 'average_DT', dt_avg )
else
tavg_exist = .false.
endif
endif
end subroutine get_src_info
!#####################################################################
!--- This routine will read source information.
!--- set up axis and field meta data for destination data file
!--- This routine is kind of long. We may change the design in the future
!--- to have shorter subroutine for easy reading.
subroutine setup_meta
integer :: rcode, natt, type, ndims, natt_dim, id_fld, dimids(4)
integer :: dim_time_dst, dims_dst(max_dims), id_axes(max_dims)
integer :: i, j, n, len, fld_dims(4), start(4), nwrite(4)
character(len=128) :: name
character(len=512) :: catt
!--- simply return if current pe is not root pe.
if(.not. is_root_pe) return
rcode = nf_create(trim(dst_data),NF_WRITE, ncid_dst)
call error_handler('error in creating file '//trim(dst_data), rcode)
!--- 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
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 = ni_dst
case('Y')
len = nj_dst
case('Z')
len = nk_src
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
call define_field_metadata(ncid_src, ncid_dst, fld_name(n), ndims, fld_dims, id_dst_fld(n))
enddo
!--- define time avg info, suppose the time_avg_info is specified by average_T1, average_T2, average_DT.
if(tavg_exist) then
fld_dims(1) = dim_time_dst
call define_field_metadata(ncid_src, ncid_dst, 'average_T1', 1, fld_dims, id_t1_dst)
call define_field_metadata(ncid_src, ncid_dst, 'average_T2', 1, fld_dims, id_t2_dst)
call define_field_metadata(ncid_src, ncid_dst, 'average_DT', 1, fld_dims, id_dt_dst)
endif
rcode = nf_enddef(ncid_dst)
!--- write axis data to dst_data
start = 1; nwrite = 1
do j = 1, ndims_dst
select case(dims_cart(j))
case('X')
nwrite(1) = ni_dst
if(dims_pos(j) == 'T' .or. dims_pos(j) == 'N') then
rcode = nf_put_vara_double(ncid_dst,id_axes(j), start, nwrite, xt_dst )
else
rcode = nf_put_vara_double(ncid_dst,id_axes(j), start, nwrite, xb_dst(2:) )
endif
call error_handler('error in putting data of '//trim(dims_name(j))//' to file '//trim(dst_data), rcode)
case('Y')
nwrite(1) = nj_dst
if(dims_pos(j) == 'T' .or. dims_pos(j) == 'E') then
rcode = nf_put_vara_double(ncid_dst,id_axes(j), start, nwrite, yt_dst )
else
rcode = nf_put_vara_double(ncid_dst,id_axes(j), start, nwrite, yb_dst(2:) )
endif
call error_handler('error in putting data of '//trim(dims_name(j))//' to file '//trim(dst_data), rcode)
case('Z')
nwrite(1) = nk_src
rcode = nf_put_vara_double(ncid_dst,id_axes(j), start, nwrite, zt_src )
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
!#####################################################################
!-------------------------------------------------------------------
!--- read src data, remap onto current grid and write output -------
subroutine process_data
integer :: n, m, l, rcode, nk, nf
integer :: start(4), nread(4), nwrite(4), id_fld(max_flds)
real, dimension(:,:,:,:), allocatable :: data_src, data_dst
real, dimension(:,:,:), allocatable :: wrk1, wrk2
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
do m = 1, ntime_src
write(stdout(),*)'************************************************'
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 = 1
if(has_zaxis(n)) nk = nk_src
allocate(data_src(ni_src, nj_src, nk,nf), data_dst(ni_dst, nj_dst, nk,nf) )
allocate(wrk1(isc:iec,jsc:jec,nk), wrk2(1:ni_dst,jstart:nj_dst,nk) )
start = 1; nread = 1
nread(1) = ni_src; nread(2) = nj_src
if(has_zaxis(n)) then
nread(3) = nk; 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 get data 1 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
rcode = nf_get_vara_double(ncid_src,id_fld(n+1), start, nread, data_src(:,:,:,2) )
call error_handler('error in get data 2 from file '//trim(src_data), rcode )
call rotate_data(data_src(:,jstart:,:,:), cosrot(:,jstart:), sinrot(:,jstart:), missing_value(n) )
endif
do l = 1, nf
!---- doing horizontal interpolaiton
if(fld_pos(n) == 'T' ) then !--- this case data is already global data
call conserve_interp(data_src(:,:,:,1), wrk2(:,:,:), missing_value(1))
else
call mpp_domains_set_stack_size(2*ni_dst*nj_dst*nk)
call nonconserve_interp(data_src(:,:,:,l), wrk1(:,:,:), fld_pos(n), missing_value(1) )
!---- get the global data
call mpp_global_field(Domain, wrk1, wrk2 )
endif
!--- assign data to global region
if(jstart == 1) then
data_dst(:,:,:,l) = wrk2
else
data_dst(:,1:jstart-1,:,l) = data_src(:,1:jstart-1,:,l)
data_dst(:,jstart:,:,l) = wrk2(:,:,:)
endif
enddo
!--- write out data from root_pe
if(is_root_pe) then
do l = 1, nf
start = 1; nwrite = 1
nwrite(1) = ni_dst; nwrite(2) = nj_dst
if(has_zaxis(n)) then
nwrite(3) = nk; 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
!--- for debugging reproducing ability accross processor count.
if(debug) write(stdout(),*)'NOTE:after regrid, chksum for '//trim(fld_name(n+l-1))// &
' at time level ',m,' ====> ',mpp_chksum(data_dst(:,:,:,l), (/mpp_root_pe()/) )
enddo
endif
n = n+nf
deallocate(data_src, data_dst, wrk1, wrk2)
enddo
enddo
rcode = nf_close(ncid_src)
if(is_root_pe) rcode = nf_close(ncid_dst)
end subroutine process_data
!#####################################################################
subroutine regrid_end
call fms_io_exit
call fms_end()
end subroutine regrid_end
!#####################################################################
! 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
!#####################################################################
! get the dimension length of any one dimensional variable
function get_dimlen(ncid, name)
integer, intent(in) :: ncid
character(len=*), intent(in) :: name
integer :: get_dimlen
integer :: varid, rcode, dims(1)
get_dimlen = 0
rcode = nf_inq_varid(ncid, trim(name), varid)
call error_handler('error in inquiring variable id of '//trim(name), rcode)
rcode = nf_inq_vardimid(ncid, varid, dims)
call error_handler('error in inquiring dimension id of '//trim(name), rcode)
rcode = nf_inq_dimlen(ncid, dims(1), get_dimlen)
call error_handler('error in inquiring dimension length of '//trim(name), rcode)
end function get_dimlen
!#####################################################################
! read the 1d integer data from netcdf file.
subroutine get_var_int_1d(ncid, name, data)
integer, intent(in) :: ncid
character(len=*), intent(in) :: name
integer, dimension(:), intent(out) :: data
integer :: rcode, varid
rcode = nf_inq_varid(ncid, name, varid)
call error_handler('error in inquiring variable id of '//trim(name), rcode)
rcode = nf_get_var_int(ncid, varid, data)
call error_handler('error in reading data of '//trim(name), rcode)
end subroutine get_var_int_1d
!#####################################################################
! read the 1d real data from netcdf file.
subroutine get_var_real_1d(ncid, name, data)
integer, intent(in) :: ncid
character(len=*), intent(in) :: name
real, dimension(:), intent(out) :: data
integer :: rcode, varid
rcode = nf_inq_varid(ncid, name, varid)
call error_handler('error in inquiring variable id of '//trim(name), rcode)
rcode = nf_get_var_double(ncid, varid, data)
call error_handler('error in reading data of '//trim(name), rcode)
end subroutine get_var_real_1d
!#####################################################################
! read the 2d real data from netcdf file.
subroutine get_var_real_2d(ncid, name, data)
integer, intent(in) :: ncid
character(len=*), intent(in) :: name
real, dimension(:,:), intent(out) :: data
integer :: rcode, varid
rcode = nf_inq_varid(ncid, name, varid)
call error_handler('error in inquiring variable id of '//trim(name), rcode)
rcode = nf_get_var_double(ncid, varid, data)
call error_handler('error in reading data of '//trim(name), rcode)
end subroutine get_var_real_2d
!#####################################################################
!--- rotate the data when the field is a vector field.
subroutine rotate_data(data, rotcos, rotsin, missing)
real, dimension(:,:,:,:), intent(inout) :: data
real, dimension(:,:), intent(in) :: rotsin, rotcos
real, intent(in) :: missing
integer :: i, j, k
real :: temp1, temp2
if(size(data,4).ne. 2 ) call error_handler('number of data should be paired when do rotation. ')
do k = 1, size(data,3)
do j = 1, size(data,2)
do i = 1, size(data,1)
if(data(i,j,k,1) /= missing ) then ! assume data(:,:,:,1) is not missing, then
! data(:,:,:,2) is not missing
temp1 = data(i,j,k,1)*rotcos(i,j)-data(i,j,k,2)*rotsin(i,j)
temp2 = data(i,j,k,1)*rotsin(i,j)+data(i,j,k,2)*rotcos(i,j)
data(i,j,k,1) = temp1
data(i,j,k,2) = temp2
endif
enddo
enddo
enddo
end subroutine rotate_data
!#####################################################################
subroutine conserve_interp(input, output, missing)
real, dimension(:,:,:), intent(in) :: input
real, dimension(:,jstart:,:), intent(out) :: output
real, intent(in) :: missing
integer :: i, j, k
!--- put data on exchange grid
output = 0.0
area_ocn = 0.0
do k = 1, size(input,3)
area_ocn = 0.0
do i = 1, size(i_atm_atmxocn)
if(j_atm_atmxocn(i) >= jstart ) then
if(input(i_ocn_atmxocn(i),j_ocn_atmxocn(i),k) /= missing) then
output(i_atm_atmxocn(i),j_atm_atmxocn(i),k) = output(i_atm_atmxocn(i),j_atm_atmxocn(i),k) &
+ area_atmxocn(i)*input(i_ocn_atmxocn(i),j_ocn_atmxocn(i),k)
area_ocn(i_atm_atmxocn(i),j_atm_atmxocn(i)) = area_ocn(i_atm_atmxocn(i),j_atm_atmxocn(i)) &
+ area_atmxocn(i)
endif
endif
enddo
do j = jstart, nj_dst
do i = 1, ni_dst
if(area_ocn(i,j) > epsln ) then
output(i,j,k) = output(i,j,k)/area_ocn(i,j)
else
output(i,j,k) = missing
endif
enddo
enddo
enddo
end subroutine conserve_interp
!#####################################################################
subroutine nonconserve_interp(input, output, position, missing)
real, dimension(:,:,:), intent(in) :: input
real, dimension(:,:,:), intent(out) :: output
character(len=1), intent(in) :: position
real, intent(in) :: missing
integer :: i, j, k
real, dimension(size(input,1), size(input,2)) :: mask
type(horiz_interp_type), pointer :: Interp=>NULL()
select case(position)
case('C')
Interp => Interp_c
case('E')
Interp => Interp_e
case('N')
Interp => Interp_n
end select
do k = 1, size(input,3)
!--- define land/sea mask
do j = 1, size(input,2)
do i = 1, size(input,1)
if(input(i,j,k) == missing) then
mask(i,j) = 0.0
else
mask(i,j) = 1.0
endif
enddo
enddo
!--- remap data onto destination horizontal grid through horiz_interp
call horiz_interp(Interp, input(:,:,k), output(:,:,k), &
mask_in=mask, missing_value=missing )
enddo
return
end subroutine nonconserve_interp
!#####################################################################
!--- To determine the grid is in new grid format or old grid format
function check_is_new_grid(ncid)
integer, intent(in) :: ncid
logical :: check_is_new_grid
integer :: rcode, id_fld
rcode = nf_inq_varid(ncid, 'x_T', id_fld)
if(rcode == 0) then
check_is_new_grid = .true.
else
rcode = nf_inq_varid(ncid, 'geolon_t', id_fld)
call error_handler('check_is_new_grid: error in inquring id of geolon_t ', rcode)
check_is_new_grid = .false.
endif
end function check_is_new_grid
!#####################################################################
!--- define a field metadata and copy the attribute from source data to destination data
subroutine define_field_metadata(ncid_in, ncid_out, name, ndim, dims, id_fld_out)
integer, intent(in) :: ncid_in, ncid_out, ndim, dims(:)
character(len=*), intent(in) :: name
integer, intent(out) :: id_fld_out
integer :: rcode, id_fld_in, natt, i
character(len=128) :: attname
rcode = nf_def_var(ncid_out, trim(name), NF_DOUBLE, ndim, dims(1:ndim), id_fld_out)
call error_handler('define_field_metadata: error in defining field '//trim(name), rcode )
rcode = nf_inq_varid(ncid_in, trim(name), id_fld_in)
call error_handler('define_field_metadata: error in inquring field of '//trim(name), rcode)
rcode = nf_inq_varnatts(ncid_in, id_fld_in, natt)
call error_handler('define_field_metadata: error in inquring natt of '//trim(name), rcode)
do i = 1, natt
rcode = nf_inq_attname(ncid_in, id_fld_in, i, attname)
call error_handler('error in inquring '//trim(name)//' attribute '//trim(attname), rcode)
if(trim(attname) /= '_FillValue') then
rcode = nf_copy_att(ncid_in,id_fld_in, attname, ncid_out, id_fld_out)
call error_handler('error in copy '//trim(name)//' attribute '//trim(attname), rcode)
endif
enddo
end subroutine define_field_metadata
!#####################################################################
end program regrid