!***********************************************************************
!* GNU Lesser General Public License
!*
!* This file is part of the GFDL Flexible Modeling System (FMS).
!*
!* FMS is free software: you can redistribute it and/or modify it under
!* the terms of the GNU Lesser General Public License as published by
!* the Free Software Foundation, either version 3 of the License, or (at
!* your option) any later version.
!*
!* 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 Lesser General Public
!* License along with FMS. If not, see .
!***********************************************************************
! Now only test some simple test, will test cubic grid mosaic in the future.
program xgrid_test
use mpp_mod, only : mpp_pe, mpp_npes, mpp_error, FATAL, mpp_chksum, mpp_min, mpp_max
use mpp_mod, only : mpp_set_current_pelist, mpp_declare_pelist, mpp_sync
use mpp_mod, only : mpp_root_pe, mpp_broadcast, stdout, NOTE, mpp_sync_self, input_nml_file
use mpp_domains_mod, only : mpp_define_domains, mpp_define_layout, mpp_domains_exit
use mpp_domains_mod, only : mpp_get_compute_domain, domain2d, mpp_domains_init
use mpp_domains_mod, only : mpp_define_mosaic_pelist, mpp_define_mosaic, mpp_global_sum
use mpp_domains_mod, only : mpp_get_data_domain, mpp_get_global_domain, mpp_update_domains
use mpp_domains_mod, only : domainUG, mpp_define_unstruct_domain, mpp_get_ug_compute_domain
use mpp_domains_mod, only : mpp_pass_ug_to_sg, mpp_pass_sg_to_ug
use mpp_io_mod, only : mpp_open, MPP_RDONLY,MPP_NETCDF, MPP_MULTI, MPP_SINGLE, mpp_close
use mpp_io_mod, only : mpp_get_att_value
use fms_mod, only : fms_init, file_exist, field_exist, field_size, open_namelist_file
use fms_mod, only : check_nml_error, close_file, read_data, stdout, fms_end
use fms_mod, only : get_mosaic_tile_grid, write_data, set_domain
use fms_io_mod, only : fms_io_exit, nullify_domain, set_domain
use constants_mod, only : DEG_TO_RAD
use xgrid_mod, only : xgrid_init, setup_xmap, put_to_xgrid, get_from_xgrid
use xgrid_mod, only : xmap_type, xgrid_count, grid_box_type, SECOND_ORDER
use xgrid_mod, only : get_xmap_grid_area, set_frac_area
use xgrid_mod, only : get_from_xgrid_ug, put_to_xgrid_ug
use mosaic_mod, only : get_mosaic_ntiles, get_mosaic_grid_sizes
use mosaic_mod, only : get_mosaic_ncontacts, get_mosaic_contact
use grid_mod, only : get_grid_comp_area
use gradient_mod, only : calc_cubic_grid_info
use ensemble_manager_mod, only : ensemble_manager_init, ensemble_pelist_setup
use ensemble_manager_mod, only : get_ensemble_size
implicit none
#include
real, parameter :: EPSLN = 1.0e-10
character(len=256) :: atm_input_file = "INPUT/atmos_input.nc"
character(len=256) :: atm_output_file = "atmos_output.nc"
character(len=256) :: lnd_output_file = "land_output.nc"
character(len=256) :: ice_output_file = "ocean_output.nc"
character(len=256) :: atm_field_name = "depth"
character(len=256) :: runoff_input_file = "INPUT/land_runoff.nc"
character(len=256) :: runoff_output_file = "land_runoff.nc"
character(len=256) :: runoff_field_name = "none"
integer :: num_iter = 1
integer :: nk_lnd = 14, nk_ice = 6
integer :: atm_layout(2) = (/2,2/)
integer :: lnd_layout(2) = (/2,2/)
integer :: ice_layout(2) = (/12,3/)
integer :: atm_nest_layout(2) = (/0,0/)
integer :: atm_npes = 0
integer :: lnd_npes = 0
integer :: ice_npes = 0
integer :: ocn_npes = 0
integer :: atm_nest_npes = 0
logical :: concurrent = .false.
logical :: test_unstruct = .true.
namelist /xgrid_test_nml/ atm_input_file, atm_field_name, runoff_input_file, runoff_field_name, num_iter, &
nk_lnd, nk_ice, atm_layout, ice_layout, lnd_layout, atm_nest_layout, &
atm_nest_npes, atm_npes, lnd_npes, ice_npes, test_unstruct
integer :: remap_method
integer :: pe, npes, ierr, nml_unit, io, n
integer :: siz(4), ntile_lnd, ntile_atm, ntile_ice, ncontact
integer, allocatable :: layout(:,:), global_indices(:,:)
integer, allocatable :: atm_nx(:), atm_ny(:), ice_nx(:), ice_ny(:), lnd_nx(:), lnd_ny(:)
integer, allocatable :: pe_start(:), pe_end(:), dummy(:)
integer, allocatable :: istart1(:), iend1(:), jstart1(:), jend1(:), tile1(:)
integer, allocatable :: istart2(:), iend2(:), jstart2(:), jend2(:), tile2(:)
character(len=256) :: grid_file = "INPUT/grid_spec.nc"
character(len=256) :: atm_mosaic, ocn_mosaic, lnd_mosaic
character(len=256) :: atm_mosaic_file, ocn_mosaic_file, lnd_mosaic_file
character(len=256) :: grid_descriptor, tile_file
integer :: isc_atm, iec_atm, jsc_atm, jec_atm, nxc_atm, nyc_atm
integer :: isc_lnd, iec_lnd, jsc_lnd, jec_lnd
integer :: isc_ice, iec_ice, jsc_ice, jec_ice
integer :: isd_atm, ied_atm, jsd_atm, jed_atm
integer :: unit, i, j, nxa, nya, nxgrid, nxl, nyl, out_unit, k
type(domain2d) :: Atm_domain, Ice_domain, Lnd_domain
type(xmap_type) :: Xmap, Xmap_runoff
type(grid_box_type) :: atm_grid
real, allocatable :: xt(:,:), yt(:,:) ! on T-cell data domain
real, allocatable :: xc(:,:), yc(:,:) ! on C-cell compute domain
real, allocatable :: tmpx(:,:), tmpy(:,:)
real, allocatable :: atm_data_in(:,:), atm_data_out(:,:)
real, allocatable :: atm_data_out_1(:,:), atm_data_out_2(:,:), atm_data_out_3(:,:)
real, allocatable :: lnd_data_out(:,:,:), ice_data_out(:,:,:)
real, allocatable :: runoff_data_in(:,:), runoff_data_out(:,:,:)
real, allocatable :: atm_area(:,:), lnd_area(:,:), ice_area(:,:)
real, allocatable :: lnd_frac(:,:,:), ice_frac(:,:,:)
real, allocatable :: x_1(:), x_2(:), x_3(:), x_4(:)
real :: sum_atm_in, sum_ice_out, sum_lnd_out, sum_atm_out
real :: sum_runoff_in, sum_runoff_out, tot
real :: min_atm_in, max_atm_in, min_atm_out, max_atm_out
real :: min_x, max_x
logical :: atm_input_file_exist, runoff_input_file_exist
integer :: npes_per_tile
integer :: id_put_side1_to_xgrid, id_get_side1_from_xgrid
integer :: id_put_side2_to_xgrid, id_get_side2_from_xgrid
integer :: ens_siz(6), ensemble_size
integer :: atm_root_pe, lnd_root_pe, ocn_root_pe, ice_root_pe, atm_nest_root_pe
integer :: atm_global_npes, ntile_atm_global, ncontact_global
integer :: atm_nest_nx, atm_nest_ny, ntile_atm_nest
logical :: atm_pe, lnd_pe, ice_pe, ocn_pe, atm_global_pe, atm_nest_pe
integer, allocatable :: atm_pelist(:), ocn_pelist(:), ice_pelist(:), lnd_pelist(:)
integer, allocatable :: atm_global_pelist(:), atm_nest_pelist(:)
integer, allocatable :: tile_id(:)
call fms_init
call mpp_domains_init
call xgrid_init(remap_method)
call ensemble_manager_init()
npes = mpp_npes()
pe = mpp_pe()
out_unit = stdout()
#ifdef INTERNAL_FILE_NML
read (input_nml_file, xgrid_test_nml, iostat=io)
ierr = check_nml_error(io, 'xgrid_test_nml')
#else
if (file_exist('input.nml')) then
ierr=1
nml_unit = open_namelist_file()
do while (ierr /= 0)
read(nml_unit, nml=xgrid_test_nml, iostat=io, end=10)
ierr = check_nml_error(io, 'xgrid_test_nml')
enddo
10 call close_file(nml_unit)
endif
#endif
!--- get ensemble size
ens_siz = get_ensemble_size()
ensemble_size = ens_siz(1)
npes = ens_siz(2)
if( atm_npes == 0 ) atm_npes = mpp_npes()
if( lnd_npes == 0 ) lnd_npes = atm_npes
if( ocn_npes == 0 ) ocn_npes = atm_npes
if( ice_npes == 0 ) ice_npes = atm_npes
if(lnd_npes > atm_npes) call mpp_error(FATAL, 'xgrid_test: lnd_npes > atm_npes')
if(ocn_npes > atm_npes) call mpp_error(FATAL, 'xgrid_test: ocn_npes > atm_npes')
write(out_unit, *) "NOTE from test_xgrid: atm_npes = ", atm_npes
write(out_unit, *) "NOTE from test_xgrid: atm_nest_npes = ", atm_nest_npes
write(out_unit, *) "NOTE from test_xgrid: lnd_npes = ", lnd_npes
write(out_unit, *) "NOTE from test_xgrid: ice_npes = ", ice_npes
write(out_unit, *) "NOTE from test_xgrid: ocn_npes = ", ocn_npes
allocate(atm_pelist(atm_npes))
allocate(ice_pelist(ice_npes))
allocate(lnd_pelist(lnd_npes))
allocate(ocn_pelist(ocn_npes))
concurrent = .false.
call ensemble_pelist_setup(concurrent, atm_npes, ocn_npes, lnd_npes, ice_npes, &
atm_pelist, ocn_pelist, lnd_pelist, ice_pelist)
if(atm_nest_npes > 0) then
atm_global_npes = atm_npes - atm_nest_npes
allocate(atm_global_pelist(atm_global_npes))
allocate(atm_nest_pelist(atm_nest_npes))
atm_global_pelist(1:atm_global_npes) = atm_pelist(1:atm_global_npes)
atm_nest_pelist(1:atm_nest_npes) = atm_pelist(atm_global_npes+1:atm_npes)
atm_global_pe = ANY(atm_global_pelist == mpp_pe())
atm_nest_pe = ANY(atm_nest_pelist == mpp_pe())
atm_nest_root_pe = atm_nest_pelist(1)
call mpp_declare_pelist(atm_global_pelist, "atm global pelist")
call mpp_declare_pelist(atm_nest_pelist, "atm nest pelist")
else
atm_global_npes = atm_npes
allocate(atm_global_pelist(atm_global_npes))
atm_global_pelist(1:atm_global_npes) = atm_pelist(1:atm_global_npes)
atm_global_pe = ANY(atm_global_pelist == mpp_pe())
atm_nest_pe = .FALSE.
endif
atm_pe = ANY(atm_pelist .EQ. mpp_pe())
lnd_pe = ANY(lnd_pelist .EQ. mpp_pe())
ocn_pe = ANY(ocn_pelist .EQ. mpp_pe())
ice_pe = ANY(ice_pelist .EQ. mpp_pe())
atm_root_pe = atm_pelist(1)
lnd_root_pe = lnd_pelist(1)
ice_root_pe = ice_pelist(1)
ocn_root_pe = ocn_pelist(1)
ntile_atm = 1
ntile_ice = 1
ntile_lnd = 1
if(field_exist(grid_file, "AREA_ATM" ) ) then
if( atm_nest_npes > 0 ) call mpp_error(FATAL, &
'xgrid_test: nested atmosphere model is only supported for mosaic grid')
allocate(atm_nx(1), atm_ny(1))
allocate(lnd_nx(1), lnd_ny(1))
allocate(ice_nx(1), ice_ny(1))
call field_size(grid_file, "AREA_ATM", siz )
atm_nx = siz(1); atm_ny = siz(2)
call field_size(grid_file, "AREA_OCN", siz )
ice_nx = siz(1); ice_ny = siz(2)
call field_size(grid_file, "AREA_LND", siz )
lnd_nx = siz(1); lnd_ny = siz(2)
if( atm_layout(1)*atm_layout(2) .NE. npes ) then
call mpp_define_layout( (/1,atm_nx,1,atm_ny/), npes, atm_layout)
endif
call mpp_define_domains( (/1,atm_nx,1,atm_ny/), atm_layout, Atm_domain, name="atmosphere")
if( lnd_layout(1)*lnd_layout(2) .NE. npes ) then
call mpp_define_layout( (/1,lnd_nx,1,lnd_ny/), npes, lnd_layout)
endif
call mpp_define_domains( (/1,lnd_nx,1,lnd_ny/), lnd_layout, Lnd_domain, name="land")
if( ice_layout(1)*ice_layout(2) .NE. npes ) then
call mpp_define_layout( (/1,ice_nx,1,ice_ny/), npes, ice_layout)
endif
call mpp_define_domains( (/1,ice_nx,1,ice_ny/), ice_layout, Ice_domain, name="Ice")
else if (field_exist(grid_file, "atm_mosaic" ) ) then
!--- Get the mosaic data of each component model
call read_data(grid_file, 'atm_mosaic', atm_mosaic)
call read_data(grid_file, 'lnd_mosaic', lnd_mosaic)
call read_data(grid_file, 'ocn_mosaic', ocn_mosaic)
atm_mosaic_file = 'INPUT/'//trim(atm_mosaic)//'.nc'
lnd_mosaic_file = 'INPUT/'//trim(lnd_mosaic)//'.nc'
ocn_mosaic_file = 'INPUT/'//trim(ocn_mosaic)//'.nc'
ntile_lnd = get_mosaic_ntiles(lnd_mosaic_file);
ntile_ice = get_mosaic_ntiles(ocn_mosaic_file);
ntile_atm = get_mosaic_ntiles(atm_mosaic_file);
if(ntile_ice > 1) call mpp_error(FATAL, &
'xgrid_test: there is more than one tile in ocn_mosaic, which is not implemented yet')
write(out_unit,*)" There is ", ntile_atm, " tiles in atmos mosaic"
write(out_unit,*)" There is ", ntile_lnd, " tiles in land mosaic"
write(out_unit,*)" There is ", ntile_ice, " tiles in ocean mosaic"
allocate(atm_nx(ntile_atm), atm_ny(ntile_atm))
allocate(lnd_nx(ntile_lnd), lnd_ny(ntile_lnd))
allocate(ice_nx(ntile_ice), ice_ny(ntile_ice))
call get_mosaic_grid_sizes(atm_mosaic_file, atm_nx, atm_ny)
call get_mosaic_grid_sizes(lnd_mosaic_file, lnd_nx, lnd_ny)
call get_mosaic_grid_sizes(ocn_mosaic_file, ice_nx, ice_ny)
ncontact = get_mosaic_ncontacts(atm_mosaic_file)
if(ncontact > 0) then
allocate(tile1(ncontact), tile2(ncontact) )
allocate(istart1(ncontact), iend1(ncontact) )
allocate(jstart1(ncontact), jend1(ncontact) )
allocate(istart2(ncontact), iend2(ncontact) )
allocate(jstart2(ncontact), jend2(ncontact) )
call get_mosaic_contact( atm_mosaic_file, tile1, tile2, istart1, iend1, jstart1, jend1, &
istart2, iend2, jstart2, jend2)
endif
ntile_atm_global = ntile_atm
ncontact_global = ncontact
if( atm_nest_npes > 0 ) then
if(ntile_atm .NE. 7) call mpp_error(FATAL, &
'xgrid_test: ntile_atm should be 7 when atmos_nest_npes > 0')
if(ncontact .NE. 13 ) call mpp_error(FATAL, &
'xgrid_test: ncontact_atm should be 13 when atmos_nest_npes > 0')
ntile_atm_global = ntile_atm_global - 1
ncontact_global = ncontact_global - 1
endif
if(atm_global_pe) then
call mpp_set_current_pelist(atm_global_pelist)
if(mod(atm_global_npes, ntile_atm_global) .NE. 0 ) call mpp_error(FATAL, &
"atm_npes_global should be divided by ntile_atm_global")
allocate(pe_start(ntile_atm_global), pe_end(ntile_atm_global) )
allocate(global_indices(4, ntile_atm_global), layout(2,ntile_atm_global))
npes_per_tile = atm_global_npes/ntile_atm_global
do n = 1, ntile_atm_global
pe_start(n) = atm_root_pe + (n-1)*npes_per_tile
pe_end(n) = atm_root_pe + n*npes_per_tile - 1
global_indices(:,n) = (/1, atm_nx(n), 1, atm_ny(n)/)
if(atm_layout(1)*atm_layout(2) == npes_per_tile ) then
layout(:,n) = atm_layout(:)
else
call mpp_define_layout( global_indices(:,n), pe_end(n)-pe_start(n)+1, layout(:,n))
endif
end do
allocate(tile_id(ntile_atm_global))
do n = 1, ntile_atm_global
tile_id(n) = n
enddo
call mpp_define_mosaic(global_indices, layout, Atm_domain, ntile_atm_global, ncontact_global, &
tile1(1:ncontact_global), tile2(1:ncontact_global), &
istart1(1:ncontact_global), iend1(1:ncontact_global), &
jstart1(1:ncontact_global), jend1(1:ncontact_global), &
istart2(1:ncontact_global), iend2(1:ncontact_global), &
jstart2(1:ncontact_global), jend2(1:ncontact_global), &
pe_start, pe_end, whalo=1, ehalo=1, shalo=1, nhalo=1, &
tile_id=tile_id, name="atmosphere")
deallocate( pe_start, pe_end, global_indices, layout, tile_id )
endif
if( atm_nest_pe ) then
atm_nest_nx = atm_nx(ntile_atm)
atm_nest_ny = atm_ny(ntile_atm)
call mpp_set_current_pelist(atm_nest_pelist)
ntile_atm_nest = 1
ncontact = 0
allocate(pe_start(1), pe_end(1) )
allocate(global_indices(4, 1), layout(2,1))
pe_start(1) = atm_nest_root_pe
pe_end(1) = atm_nest_root_pe + atm_nest_npes - 1
global_indices(:,1) = (/1, atm_nest_nx, 1, atm_nest_ny/) !-- the last tile is the nested tile
if(atm_nest_layout(1)*atm_nest_layout(2) == atm_nest_npes ) then
layout(:,1) = atm_nest_layout(:)
else
call mpp_define_layout( global_indices(:,1), atm_nest_npes, layout(:,1))
endif
allocate(tile_id(1))
tile_id(1) = ntile_atm
call mpp_define_mosaic(global_indices, layout, Atm_domain, ntile_atm_nest, ncontact, dummy, dummy, &
dummy, dummy, dummy, dummy, dummy, dummy, dummy, dummy, pe_start, pe_end, &
whalo=1, ehalo=1, shalo=1, nhalo=1, tile_id=tile_id, name="atmos nest")
deallocate( pe_start, pe_end, global_indices, layout, tile_id )
endif
if( lnd_pe ) then
call mpp_set_current_pelist(lnd_pelist)
ncontact = 0 ! no update is needed for land model.
allocate(pe_start(ntile_lnd), pe_end(ntile_lnd) )
allocate(global_indices(4,ntile_lnd), layout(2,ntile_lnd))
if(mod(lnd_npes, ntile_lnd) .NE. 0 ) call mpp_error(FATAL,"lnd_npes should be divided by ntile_lnd")
npes_per_tile = lnd_npes/ntile_lnd
do n = 1, ntile_lnd
pe_start(n) = lnd_root_pe + (n-1)*npes_per_tile
pe_end(n) = lnd_root_pe + n*npes_per_tile - 1
global_indices(:,n) = (/1, lnd_nx(n), 1, lnd_ny(n)/)
if(lnd_layout(1)*lnd_layout(2) == npes_per_tile ) then
layout(:,n) = lnd_layout(:)
else
call mpp_define_layout( global_indices(:,n), npes_per_tile, layout(:,n))
endif
end do
ncontact = 0 ! no update is needed for land and ocean model.
allocate(tile_id(ntile_lnd))
do n = 1, ntile_lnd
tile_id(n) = n
enddo
call mpp_define_mosaic(global_indices, layout, Lnd_domain, ntile_lnd, ncontact, dummy, dummy, &
dummy, dummy, dummy, dummy, dummy, dummy, dummy, dummy, pe_start, pe_end, tile_id=tile_id, name="land")
deallocate( pe_start, pe_end, global_indices, layout, tile_id )
endif
if( ice_pe ) then
call mpp_set_current_pelist(ice_pelist)
ncontact = 0 ! no update is needed for ocn model.
allocate(pe_start(ntile_ice), pe_end(ntile_ice) )
allocate(global_indices(4, ntile_ice), layout(2, ntile_ice))
npes_per_tile = ice_npes/ntile_ice
do n = 1, ntile_ice
pe_start(n) = ice_root_pe + (n-1)*npes_per_tile
pe_end(n) = ice_root_pe + n*npes_per_tile - 1
global_indices(:,n) = (/1, ice_nx(n), 1, ice_ny(n)/)
if(ice_layout(1)*ice_layout(2) == npes_per_tile ) then
layout(:,n) = ice_layout(:)
else
call mpp_define_layout( global_indices(:,n), ice_npes, layout(:,n))
endif
end do
allocate(tile_id(ntile_ice))
do n = 1, ntile_ice
tile_id(n) = n
enddo
call mpp_define_mosaic(global_indices, layout, Ice_domain, ntile_ice, ncontact, dummy, dummy, &
dummy, dummy, dummy, dummy, dummy, dummy, dummy, dummy, pe_start, pe_end, tile_id=tile_id, name="Ice")
deallocate( pe_start, pe_end, global_indices, layout, tile_id )
endif
else
call mpp_error(FATAL, 'test_xgrid:both AREA_ATM and atm_mosaic does not exist in '//trim(grid_file))
end if
if( atm_pe ) then
call mpp_get_compute_domain(atm_domain, isc_atm, iec_atm, jsc_atm, jec_atm)
call mpp_get_data_domain(atm_domain, isd_atm, ied_atm, jsd_atm, jed_atm)
call mpp_get_global_domain(atm_domain, xsize = nxa, ysize = nya)
nxc_atm = iec_atm - isc_atm + 1
nyc_atm = jec_atm - jsc_atm + 1
endif
if( lnd_pe ) then
call mpp_get_compute_domain(lnd_domain, isc_lnd, iec_lnd, jsc_lnd, jec_lnd)
call mpp_get_global_domain(lnd_domain, xsize = nxl, ysize = nyl)
endif
if( ice_pe ) then
call mpp_get_compute_domain(Ice_domain, isc_ice, iec_ice, jsc_ice, jec_ice)
endif
! set up atm_grid for second order conservative interpolation and atm grid is cubic grid.
if(remap_method == SECOND_ORDER ) then
if( atm_nest_npes > 0 ) call mpp_error(FATAL, &
'test_xgrid: remap_method could not be SECOND_ORDER when atmos_nest_npes > 0')
if(ntile_atm == 6) then ! 6 tile for cubic grid
allocate(xt (isd_atm:ied_atm,jsd_atm:jed_atm), yt (isd_atm:ied_atm,jsd_atm:jed_atm) )
allocate(xc (isc_atm:ied_atm,jsc_atm:jed_atm), yc (isc_atm:ied_atm,jsc_atm:jed_atm) )
allocate(atm_grid%dx(isc_atm:iec_atm,jsc_atm:jed_atm), atm_grid%dy(isc_atm:iec_atm+1,jsc_atm:jec_atm) )
allocate(atm_grid%edge_w(jsc_atm:jed_atm), atm_grid%edge_e(jsc_atm:jed_atm))
allocate(atm_grid%edge_s(isc_atm:ied_atm), atm_grid%edge_n(isc_atm:ied_atm))
allocate(atm_grid%en1 (3,isc_atm:iec_atm,jsc_atm:jed_atm), atm_grid%en2 (3,isc_atm:ied_atm,jsc_atm:jec_atm) )
allocate(atm_grid%vlon(3,isc_atm:iec_atm,jsc_atm:jec_atm), atm_grid%vlat(3,isc_atm:iec_atm,jsc_atm:jec_atm) )
allocate(atm_grid%area(isc_atm:iec_atm,jsc_atm:jec_atm) )
! first get grid from grid file
call get_mosaic_tile_grid(tile_file, atm_mosaic_file, atm_domain)
allocate(tmpx(nxa*2+1, nya*2+1), tmpy(nxa*2+1, nya*2+1))
call read_data( tile_file, 'x', tmpx, no_domain=.true.)
call read_data( tile_file, 'y', tmpy, no_domain=.true.)
xt = 0; yt = 0;
do j = jsc_atm, jec_atm
do i = isc_atm, iec_atm
xt(i,j) = tmpx(2*i, 2*j)*DEG_TO_RAD
yt(i,j) = tmpy(2*i, 2*j)*DEG_TO_RAD
end do
end do
do j = jsc_atm, jed_atm
do i = isc_atm, ied_atm
xc(i,j) = tmpx(2*i-1, 2*j-1)*DEG_TO_RAD
yc(i,j) = tmpy(2*i-1, 2*j-1)*DEG_TO_RAD
end do
end do
call mpp_update_domains(xt, atm_domain)
call mpp_update_domains(yt, atm_domain)
call calc_cubic_grid_info(xt, yt, xc, yc, atm_grid%dx, atm_grid%dy, atm_grid%area, atm_grid%edge_w, &
atm_grid%edge_e, atm_grid%edge_s, atm_grid%edge_n, atm_grid%en1, &
atm_grid%en2, atm_grid%vlon, atm_grid%vlat, isc_atm==1, iec_atm==nxa, &
jsc_atm==1, jec_atm==nya )
end if
end if
call mpp_set_current_pelist()
!--- conservation check is done in setup_xmap.
call setup_xmap(Xmap, (/ 'ATM', 'OCN', 'LND' /), (/ Atm_domain, Ice_domain, Lnd_domain /), grid_file, atm_grid)
call setup_xmap(Xmap_runoff, (/ 'LND', 'OCN'/), (/ Lnd_domain, Ice_domain/), grid_file )
!--- set frac area if nk_lnd or nk_ocn is greater than 1.
if(nk_lnd > 0 .AND. lnd_pe) then
allocate(lnd_frac(isc_lnd:iec_lnd, jsc_lnd:jec_lnd, nk_lnd))
call random_number(lnd_frac)
lnd_frac = lnd_frac + 0.5
do j = jsc_lnd, jec_lnd
do i = isc_lnd, iec_lnd
tot = sum(lnd_frac(i,j,:))
do k = 1, nk_lnd
lnd_frac(i,j,k)=lnd_frac(i,j,k)/tot
enddo
enddo
enddo
call set_frac_area(lnd_frac, 'LND', xmap)
endif
if(nk_ice > 0 ) then
if( ice_pe ) then
allocate(ice_frac(isc_ice:iec_ice, jsc_ice:jec_ice, nk_ice))
call random_number(ice_frac)
ice_frac = ice_frac + 0.5
do j = jsc_ice, jec_ice
do i = isc_ice, iec_ice
tot = sum(ice_frac(i,j,:))
do k = 1, nk_ice
ice_frac(i,j,k)=ice_frac(i,j,k)/tot
enddo
enddo
enddo
endif
call set_frac_area(ice_frac, 'OCN', xmap)
endif
if(mpp_pe() == mpp_root_pe() ) print *, '------> Calling test_unstruct_exchange'
if(test_unstruct) call test_unstruct_exchange()
if(mpp_pe() == mpp_root_pe() ) print *, '------> Finish with test_unstruct_exchange'
deallocate(atm_nx, atm_ny, lnd_nx, lnd_ny, ice_nx, ice_ny)
!--- remap "realistic" data and write the output file when atmos_input_file does exist
atm_input_file_exist = file_exist(atm_input_file, domain=atm_domain)
if( atm_input_file_exist ) then
if(trim(atm_input_file) == trim(atm_output_file) ) call mpp_error(FATAL, &
"test_xgrid: atm_input_file should have a different name from atm_output_file")
call field_size(atm_input_file, atm_field_name, siz, domain=Atm_domain )
if(siz(1) .NE. nxa .OR. siz(2) .NE. nya ) call mpp_error(FATAL,"test_xgrid: x- and y-size of field "//trim(atm_field_name) &
//" in file "//trim(atm_input_file) //" does not compabile with the grid size" )
if(siz(3) > 1) call mpp_error(FATAL,"test_xgrid: number of vertical level of field "//trim(atm_field_name) &
//" in file "//trim(atm_input_file) //" should be no larger than 1")
allocate(atm_data_in (isc_atm:iec_atm, jsc_atm:jec_atm ) )
allocate(atm_data_out(isc_atm:iec_atm, jsc_atm:jec_atm ) )
allocate(lnd_data_out(isc_lnd:iec_lnd, jsc_lnd:jec_lnd, nk_lnd) )
allocate(ice_data_out(isc_ice:iec_ice, jsc_ice:jec_ice, nk_ice) )
allocate(atm_data_out_1(isc_atm:iec_atm, jsc_atm:jec_atm ) )
allocate(atm_data_out_2(isc_atm:iec_atm, jsc_atm:jec_atm ) )
allocate(atm_data_out_3(isc_atm:iec_atm, jsc_atm:jec_atm ) )
nxgrid = max(xgrid_count(Xmap), 1)
allocate(x_1(nxgrid), x_2(nxgrid))
allocate(x_3(nxgrid), x_4(nxgrid))
x_1 = 0
x_2 = 0
x_3 = 0
x_4 = 0
atm_data_in = 0
atm_data_out = 0
lnd_data_out = 0
ice_data_out = 0
atm_data_out_1 = 0
atm_data_out_2 = 0
atm_data_out_3 = 0
! test one time level should be sufficient
call read_data(atm_input_file, atm_field_name, atm_data_in, atm_domain)
call put_to_xgrid(atm_data_in, 'ATM', x_1, Xmap, remap_method=remap_method)
call put_to_xgrid(atm_data_in, 'ATM', x_2, Xmap, remap_method=remap_method, complete=.false.)
call put_to_xgrid(atm_data_in, 'ATM', x_3, Xmap, remap_method=remap_method, complete=.false.)
call put_to_xgrid(atm_data_in, 'ATM', x_4, Xmap, remap_method=remap_method, complete=.true.)
min_x = minval(x_1)
max_x = maxval(x_1)
!--- check make sure x_2, x_3 and x_4 are the same as x_1
if(ANY(x_1 .NE. x_2)) call mpp_error(FATAL,"test_xgrid: x_1 and x_2 are not equal")
if(ANY(x_1 .NE. x_3)) call mpp_error(FATAL,"test_xgrid: x_1 and x_3 are not equal")
if(ANY(x_1 .NE. x_4)) call mpp_error(FATAL,"test_xgrid: x_1 and x_4 are not equal")
deallocate(x_3,x_4)
x_2 = 0
call get_from_xgrid(lnd_data_out, 'LND', x_1, xmap)
call get_from_xgrid(ice_data_out, 'OCN', x_1, xmap)
call put_to_xgrid(lnd_data_out, 'LND', x_2, xmap)
call put_to_xgrid(ice_data_out, 'OCN', x_2, xmap)
call get_from_xgrid(atm_data_out, 'ATM', x_2, xmap)
call get_from_xgrid(atm_data_out_1, 'ATM', x_2, xmap, complete=.false.)
call get_from_xgrid(atm_data_out_2, 'ATM', x_2, xmap, complete=.false.)
call get_from_xgrid(atm_data_out_3, 'ATM', x_2, xmap, complete=.true.)
if(ANY(atm_data_out .NE. atm_data_out_1)) &
call mpp_error(FATAL,"test_xgrid: atm_data_out and atm_data_out_1 are not equal")
if(ANY(atm_data_out .NE. atm_data_out_2)) &
call mpp_error(FATAL,"test_xgrid: atm_data_out and atm_data_out_2 are not equal")
if(ANY(atm_data_out .NE. atm_data_out_3)) &
call mpp_error(FATAL,"test_xgrid: atm_data_out and atm_data_out_3 are not equal")
call write_data( atm_output_file, atm_field_name, atm_data_out, atm_domain)
call write_data( lnd_output_file, atm_field_name, lnd_data_out, lnd_domain)
call write_data( ice_output_file, atm_field_name, ice_data_out, Ice_domain)
!--- print out checksum
write(out_unit,*) "chksum for atm_data_in", mpp_chksum(atm_data_in)
write(out_unit,*) "chksum for lnd_data_out", mpp_chksum(lnd_data_out)
write(out_unit,*) "chksum for ice_data_out", mpp_chksum(ice_data_out)
write(out_unit,*) "chksum for atm_data_out", mpp_chksum(atm_data_out)
! conservation check
allocate(atm_area(isc_atm:iec_atm, jsc_atm:jec_atm ) )
allocate(lnd_area(isc_lnd:iec_lnd, jsc_lnd:jec_lnd ) )
allocate(ice_area(isc_ice:iec_ice, jsc_ice:jec_ice ) )
call get_xmap_grid_area("ATM", Xmap, atm_area)
call get_xmap_grid_area("LND", Xmap, lnd_area)
call get_xmap_grid_area("OCN", Xmap, ice_area)
min_atm_in = minval(atm_data_in)
max_atm_in = maxval(atm_data_in)
min_atm_out = minval(atm_data_out)
max_atm_out = maxval(atm_data_out)
call mpp_min(min_atm_in)
call mpp_max(max_atm_in)
call mpp_min(min_atm_out)
call mpp_max(max_atm_out)
call mpp_min(min_x)
call mpp_max(max_x)
sum_atm_in = mpp_global_sum(atm_domain, atm_area * atm_data_in)
sum_lnd_out = 0
do k = 1, nk_lnd
sum_lnd_out = sum_lnd_out + mpp_global_sum(lnd_domain, lnd_area * lnd_data_out(:,:,k))
enddo
sum_ice_out = 0
do k = 1, nk_ice
sum_ice_out = sum_ice_out + mpp_global_sum(ice_domain, ice_area * ice_data_out(:,:,k))
enddo
sum_atm_out = mpp_global_sum(atm_domain, atm_area * atm_data_out)
write(out_unit,*) "********************** check conservation *********************** "
write(out_unit,*) "the global area sum of atmos input data is : ", sum_atm_in
write(out_unit,*) "the global area sum of atmos output data is : ", sum_atm_out
write(out_unit,*) "the global area sum of land output data + ocean output data is: ", sum_lnd_out+sum_ice_out
write(out_unit,*) "The min of atmos input data is ", min_atm_in
write(out_unit,*) "The min of xgrid data is ", min_x
write(out_unit,*) "The min of atmos output data is ", min_atm_out
write(out_unit,*) "The max of atmos input data is ", max_atm_in
write(out_unit,*) "The max of xgrid data is ", max_x
write(out_unit,*) "The max of atmos output data is ", max_atm_out
deallocate(atm_area, lnd_area, ice_area, atm_data_in, atm_data_out, lnd_data_out, ice_data_out)
deallocate(atm_data_out_1, atm_data_out_2, atm_data_out_3)
deallocate(x_1, x_2)
else
write(out_unit,*) "NOTE from test_xgrid ==> file "//trim(atm_input_file)//" does not exist, no check is done for real data sets."
end if
runoff_input_file_exist = file_exist(runoff_input_file, domain=atm_domain)
if( runoff_input_file_exist ) then
if( atm_nest_npes > 0 ) call mpp_error(FATAL, &
"test_xgrid: runoff_input_file_exist should be false when atmos_nest_npes > 0")
if(trim(runoff_input_file) == trim(runoff_output_file) ) call mpp_error(FATAL, &
"test_xgrid: runoff_input_file should have a different name from runoff_output_file")
call field_size(runoff_input_file, runoff_field_name, siz )
if(siz(1) .NE. nxl .OR. siz(2) .NE. nyl ) call mpp_error(FATAL,"test_xgrid: x- and y-size of field "//trim(runoff_field_name) &
//" in file "//trim(runoff_input_file) //" does not compabile with the grid size" )
if(siz(3) > 1) call mpp_error(FATAL,"test_xgrid: number of vertical level of field "//trim(runoff_field_name) &
//" in file "//trim(runoff_input_file) //" should be no larger than 1")
allocate(runoff_data_in (isc_lnd:iec_lnd, jsc_lnd:jec_lnd ) )
allocate(runoff_data_out(isc_ice:iec_ice, jsc_ice:jec_ice, 1) )
nxgrid = max(xgrid_count(Xmap_runoff), 1)
allocate(x_1(nxgrid), x_2(nxgrid))
runoff_data_in = 0
runoff_data_out = 0
! test one time level should be sufficient
call read_data(runoff_input_file, runoff_field_name, runoff_data_in, lnd_domain)
call put_to_xgrid(runoff_data_in, 'LND', x_1, Xmap_runoff)
call get_from_xgrid(runoff_data_out, 'OCN', x_1, xmap_runoff)
call write_data( runoff_output_file, runoff_field_name, runoff_data_out, ice_domain)
! conservation check
allocate(lnd_area(isc_lnd:iec_lnd, jsc_lnd:jec_lnd ) )
allocate(ice_area(isc_ice:iec_ice, jsc_ice:jec_ice ) )
call get_xmap_grid_area("LND", Xmap_runoff, lnd_area)
call get_xmap_grid_area("OCN", Xmap_runoff, ice_area)
sum_runoff_in = mpp_global_sum(lnd_domain, lnd_area * runoff_data_in)
sum_runoff_out = mpp_global_sum(Ice_domain, ice_area * runoff_data_out(:,:,1))
write(out_unit,*) "********************** check conservation *********************** "
write(out_unit,*) "the global area sum of runoff input data is : ", sum_runoff_in
write(out_unit,*) "the global area sum of runoff output data is : ", sum_runoff_out
else
write(out_unit,*) "NOTE from test_xgrid ==> file "//trim(runoff_input_file)//" does not exist, no check is done for real data sets."
end if
! when num_iter is greater than 0, create random number as input to test the performance of xgrid_mod.
if(num_iter > 0) then
if( atm_nest_npes > 0 ) call mpp_error(FATAL, &
"test_xgrid: num_iter > 0 when atm_nest_npes > 0")
allocate(atm_data_in (isc_atm:iec_atm, jsc_atm:jec_atm ) )
allocate(atm_data_out(isc_atm:iec_atm, jsc_atm:jec_atm ) )
allocate(lnd_data_out(isc_lnd:iec_lnd, jsc_lnd:jec_lnd, nk_lnd) )
allocate(ice_data_out(isc_ice:iec_ice, jsc_ice:jec_ice, nk_ice) )
nxgrid = max(xgrid_count(Xmap), 1)
allocate(x_1(nxgrid), x_2(nxgrid))
atm_data_in = 0
atm_data_out = 0
lnd_data_out = 0
ice_data_out = 0
allocate(atm_area(isc_atm:iec_atm, jsc_atm:jec_atm ) )
allocate(lnd_area(isc_lnd:iec_lnd, jsc_lnd:jec_lnd ) )
allocate(ice_area(isc_ice:iec_ice, jsc_ice:jec_ice ) )
call get_xmap_grid_area("ATM", Xmap, atm_area)
call get_xmap_grid_area("LND", Xmap, lnd_area)
call get_xmap_grid_area("OCN", Xmap, ice_area)
do n = 1, num_iter
call random_number(atm_data_in)
call put_to_xgrid(atm_data_in, 'ATM', x_1, Xmap, remap_method=remap_method)
call get_from_xgrid(lnd_data_out, 'LND', x_1, xmap)
call get_from_xgrid(ice_data_out, 'OCN', x_1, xmap)
call put_to_xgrid(lnd_data_out, 'LND', x_2, xmap)
call put_to_xgrid(ice_data_out, 'OCN', x_2, xmap)
call get_from_xgrid(atm_data_out, 'ATM', x_2, xmap)
sum_atm_in = mpp_global_sum(atm_domain, atm_area * atm_data_in)
sum_lnd_out = 0
do k = 1, nk_lnd
sum_lnd_out = sum_lnd_out + mpp_global_sum(lnd_domain, lnd_area * lnd_data_out(:,:,k))
enddo
sum_ice_out = 0
do k = 1, nk_ice
sum_ice_out = sum_ice_out + mpp_global_sum(Ice_domain, ice_area * ice_data_out(:,:,k))
enddo
sum_atm_out = mpp_global_sum(atm_domain, atm_area * atm_data_out)
write(out_unit,*) "********************** check conservation *********************** "
write(out_unit,*) "the global area sum of atmos input data is : ", sum_atm_in
write(out_unit,*) "the global area sum of atmos output data is : ", sum_atm_out
write(out_unit,*) "the global area sum of land output data + ocean output data is: ", sum_lnd_out+sum_ice_out
enddo
deallocate(atm_area, lnd_area, ice_area, atm_data_in, atm_data_out, lnd_data_out, ice_data_out)
deallocate(x_1, x_2)
endif
write(out_unit,*) "************************************************************************"
write(out_unit,*) "*********** Finish running program test_xgrid *************"
write(out_unit,*) "************************************************************************"
call fms_io_exit
call fms_end
contains
subroutine test_unstruct_exchange()
real, allocatable :: atm_data_in(:,:), atm_data_sg(:,:)
real, allocatable :: atm_data_sg_1(:,:), atm_data_sg_2(:,:), atm_data_sg_3(:,:)
real, allocatable :: lnd_data_sg(:,:,:), ice_data_sg(:,:,:)
real, allocatable :: atm_data_ug(:,:), tmp_sg(:,:,:)
real, allocatable :: atm_data_ug_1(:,:), atm_data_ug_2(:,:), atm_data_ug_3(:,:)
real, allocatable :: lnd_data_ug(:,:), ice_data_ug(:,:,:)
real, allocatable :: x_1(:), x_2(:), x_3(:), x_4(:)
real, allocatable :: y_1(:), y_2(:), y_3(:), y_4(:)
real, allocatable, dimension(:,:) :: rmask, tmp2d
logical, allocatable, dimension(:,:,:) :: lmask
integer, allocatable, dimension(:) :: npts_tile, grid_index, ntiles_grid
integer :: ntiles, nx, ny, ntotal_land, l, is_ug, ie_ug
type(domainUG) :: ug_domain
type(xmap_type) :: Xmap_ug, Xmap_runoff_ug
if(ntile_lnd .NE. 6) call mpp_error(FATAL, &
"test_xgrid: when test_unstruct is true, ntile_lnd must be 6")
!--- define unstructured grid domain
ntiles = ntile_lnd
nx = lnd_nx(1)
ny = lnd_ny(1)
allocate(lmask(nx,ny,ntiles))
allocate(npts_tile(ntiles))
lmask = .false.
if(mpp_pe() == mpp_root_pe() ) then
allocate(rmask(nx,ny))
!--- construct gmask.
call set_domain(Lnd_domain)
do n = 1, ntiles
rmask = 0
call get_grid_comp_area('LND', n, rmask)
do j = 1, ny
do i = 1, nx
if(rmask(i,j) > 0) then
lmask(i,j,n) = .true.
endif
enddo
enddo
npts_tile(n) = count(lmask(:,:,n))
enddo
call nullify_domain()
ntotal_land = sum(npts_tile)
allocate(grid_index(ntotal_land))
l = 0
do n = 1, ntiles
do j = 1, ny
do i = 1, nx
if(lmask(i,j,n)) then
l = l + 1
grid_index(l) = (j-1)*nx+i
endif
enddo
enddo
enddo
deallocate(rmask)
endif
call mpp_broadcast(npts_tile, ntiles, mpp_root_pe())
if(mpp_pe() .NE. mpp_root_pe()) then
ntotal_land = sum(npts_tile)
allocate(grid_index(ntotal_land))
endif
call mpp_broadcast(grid_index, ntotal_land, mpp_root_pe())
allocate(ntiles_grid(ntotal_land))
ntiles_grid = 1
!--- define the unstructured grid domain
if(lnd_pe) then
call mpp_set_current_pelist(lnd_pelist)
call mpp_define_unstruct_domain(UG_domain, Lnd_domain, npts_tile, ntiles_grid, mpp_npes(), &
1, grid_index, name="LAND unstruct")
call mpp_get_UG_compute_domain(UG_domain, is_ug, ie_ug)
endif
call mpp_set_current_pelist()
call setup_xmap(Xmap_ug, (/ 'ATM', 'OCN', 'LND' /), (/ Atm_domain, Ice_domain, Lnd_domain /), &
grid_file, atm_grid, lnd_ug_domain=UG_domain)
!--- set frac area if nk_lnd or nk_ocn is greater than 1.
if(nk_lnd > 0 .AND. lnd_pe) then
allocate(tmp2d(is_ug:ie_ug,nk_lnd))
call mpp_pass_sg_to_ug(UG_domain, lnd_frac, tmp2d)
call set_frac_area(tmp2d, 'LND', xmap_ug)
deallocate(tmp2d)
endif
if(nk_ice > 0 ) then
call set_frac_area(ice_frac, 'OCN', Xmap_ug)
endif
! call setup_xmap(Xmap_runoff_ug, (/ 'LND', 'OCN'/), (/ Lnd_domain, Ice_domain/), grid_file, lnd_ug_domain=UG_domain )
allocate(atm_data_ug(isc_atm:iec_atm, jsc_atm:jec_atm ) )
allocate(atm_data_ug_1(isc_atm:iec_atm, jsc_atm:jec_atm ) )
allocate(atm_data_ug_2(isc_atm:iec_atm, jsc_atm:jec_atm ) )
allocate(atm_data_ug_3(isc_atm:iec_atm, jsc_atm:jec_atm ) )
allocate(atm_data_in(isc_atm:iec_atm, jsc_atm:jec_atm ) )
allocate(atm_data_sg(isc_atm:iec_atm, jsc_atm:jec_atm ) )
allocate(atm_data_sg_1(isc_atm:iec_atm, jsc_atm:jec_atm ) )
allocate(atm_data_sg_2(isc_atm:iec_atm, jsc_atm:jec_atm ) )
allocate(atm_data_sg_3(isc_atm:iec_atm, jsc_atm:jec_atm ) )
atm_data_in = 0
atm_data_sg = 0
atm_data_sg_1 = 0
atm_data_sg_2 = 0
atm_data_sg_3 = 0
atm_data_ug = 0
atm_data_ug_1 = 0
atm_data_ug_2 = 0
atm_data_ug_3 = 0
if(lnd_pe) then
allocate(lnd_data_ug(is_ug:ie_ug, nk_lnd) )
allocate(lnd_data_sg(isc_lnd:iec_lnd, jsc_lnd:jec_lnd, nk_lnd) )
lnd_data_sg = 0
lnd_data_ug = 0
endif
if(ice_pe) then
allocate(ice_data_ug(isc_ice:iec_ice, jsc_ice:jec_ice, nk_ice) )
allocate(ice_data_sg(isc_ice:iec_ice, jsc_ice:jec_ice, nk_ice) )
ice_data_sg = 0
ice_data_ug = 0
endif
nxgrid = max(xgrid_count(Xmap), 1)
allocate(x_1(nxgrid), x_2(nxgrid))
allocate(x_3(nxgrid), x_4(nxgrid))
x_1 = 0
x_2 = 0
x_3 = 0
x_4 = 0
call random_number(atm_data_in)
call put_to_xgrid(atm_data_in, 'ATM', x_1, Xmap, remap_method=remap_method)
call put_to_xgrid(atm_data_in+1, 'ATM', x_2, Xmap, remap_method=remap_method, complete=.false.)
call put_to_xgrid(atm_data_in+2, 'ATM', x_3, Xmap, remap_method=remap_method, complete=.false.)
call put_to_xgrid(atm_data_in+3, 'ATM', x_4, Xmap, remap_method=remap_method, complete=.true.)
call get_from_xgrid(lnd_data_sg, 'LND', x_1, xmap)
call get_from_xgrid(ice_data_sg, 'OCN', x_1, xmap)
call put_to_xgrid(lnd_data_sg, 'LND', x_2, xmap)
call put_to_xgrid(ice_data_sg, 'OCN', x_2, xmap)
call get_from_xgrid(atm_data_sg, 'ATM', x_2, xmap)
call get_from_xgrid(atm_data_sg_1, 'ATM', x_2, xmap, complete=.false.)
call get_from_xgrid(atm_data_sg_2, 'ATM', x_2, xmap, complete=.false.)
call get_from_xgrid(atm_data_sg_3, 'ATM', x_2, xmap, complete=.true.)
nxgrid = max(xgrid_count(Xmap_ug), 1)
allocate(y_1(nxgrid), y_2(nxgrid))
allocate(y_3(nxgrid), y_4(nxgrid))
y_1 = 0
y_2 = 0
y_3 = 0
y_4 = 0
call put_to_xgrid(atm_data_in, 'ATM', y_1, Xmap_ug, remap_method=remap_method)
call put_to_xgrid(atm_data_in+1, 'ATM', y_2, Xmap_ug, remap_method=remap_method, complete=.false.)
call put_to_xgrid(atm_data_in+2, 'ATM', y_3, Xmap_ug, remap_method=remap_method, complete=.false.)
call put_to_xgrid(atm_data_in+3, 'ATM', y_4, Xmap_ug, remap_method=remap_method, complete=.true.)
call get_from_xgrid_ug(lnd_data_ug, 'LND', y_1, xmap_ug)
call get_from_xgrid(ice_data_ug, 'OCN', y_1, xmap_ug)
call put_to_xgrid_ug(lnd_data_ug, 'LND', y_2, xmap_ug)
call put_to_xgrid(ice_data_ug, 'OCN', y_2, xmap_ug)
call get_from_xgrid(atm_data_ug, 'ATM', y_1, xmap_ug)
call get_from_xgrid(atm_data_ug_1, 'ATM', y_2, xmap_ug, complete=.false.)
call get_from_xgrid(atm_data_ug_2, 'ATM', y_2, xmap_ug, complete=.false.)
call get_from_xgrid(atm_data_ug_3, 'ATM', y_2, xmap_ug, complete=.true.)
!--- comparing data ---------------------
call compare_chksum(ice_data_ug, ice_data_sg, "ice_data_out")
call compare_chksum_2D(atm_data_ug, atm_data_ug, "atm_data_out")
call compare_chksum_2D(atm_data_ug_1, atm_data_ug_1, "atm_data_out_1")
call compare_chksum_2D(atm_data_ug_2, atm_data_ug_2, "atm_data_out_2")
call compare_chksum_2D(atm_data_ug_3, atm_data_ug_3, "atm_data_out_3")
allocate(tmp_sg(isc_lnd:iec_lnd,jsc_lnd:jec_lnd,nk_lnd))
tmp_sg = 0
if(lnd_pe) then
call mpp_set_current_pelist(lnd_pelist)
call mpp_pass_ug_to_sg(ug_domain, lnd_data_ug, tmp_sg)
call compare_chksum(tmp_sg, lnd_data_sg, "lnd_data_out")
deallocate(lnd_data_sg,lnd_data_ug)
endif
call mpp_set_current_pelist()
if(ice_pe) deallocate(ice_data_sg, ice_data_ug)
deallocate(tmp_sg, x_1, x_2, x_3, x_4, y_1, y_2, y_3, y_4)
deallocate(atm_data_in, atm_data_sg)
deallocate(atm_data_sg_1, atm_data_sg_2, atm_data_sg_3)
deallocate(atm_data_ug)
deallocate(atm_data_ug_1, atm_data_ug_2, atm_data_ug_3)
end subroutine test_unstruct_exchange
!###########################################################################
subroutine compare_chksum_2D( a, b, string )
real, intent(in), dimension(:,:) :: a, b
character(len=*), intent(in) :: string
integer(LONG_KIND) :: sum1, sum2
integer :: i, j
call mpp_sync_self()
if(size(a,1) .ne. size(b,1) .or. size(a,2) .ne. size(b,2) ) &
call mpp_error(FATAL,'compare_chksum_2D: size of a and b does not match')
do j = 1, size(a,2)
do i = 1, size(a,1)
if(a(i,j) .ne. b(i,j)) then
write(*,'(a,i3,a,i3,a,i3,a,f20.9,a,f20.9)')"at pe ", mpp_pe(), &
", at point (",i,", ", j, "),a=", a(i,j), ",b=", b(i,j)
call mpp_error(FATAL, trim(string)//': point by point comparison are not OK.')
endif
enddo
enddo
sum1 = mpp_chksum( a, (/pe/) )
sum2 = mpp_chksum( b, (/pe/) )
if( sum1.EQ.sum2 )then
if( pe.EQ.mpp_root_pe() )call mpp_error( NOTE, trim(string)//': OK.' )
!--- in some case, even though checksum agree, the two arrays
! actually are different, like comparing (1.1,-1.2) with (-1.1,1.2)
!--- hence we need to check the value point by point.
else
call mpp_error( FATAL, trim(string)//': chksums are not OK.' )
end if
end subroutine compare_chksum_2D
!###########################################################################
subroutine compare_chksum( a, b, string )
real, intent(in), dimension(:,:,:) :: a, b
character(len=*), intent(in) :: string
integer(LONG_KIND) :: sum1, sum2
integer :: i, j, k
! z1l can not call mpp_sync here since there might be different number of tiles on each pe.
! mpp_sync()
call mpp_sync_self()
if(size(a,1) .ne. size(b,1) .or. size(a,2) .ne. size(b,2) .or. size(a,3) .ne. size(b,3) ) &
call mpp_error(FATAL,'compare_chksum: size of a and b does not match')
do k = 1, size(a,3)
do j = 1, size(a,2)
do i = 1, size(a,1)
if(a(i,j,k) .ne. b(i,j,k)) then
print*, "a,b=", a(i,j,k), b(i,j,k), i,j,k
write(*, '(a,i3,a,i3,a,i3,a,i3,a,f20.9,a,f20.9)')" at pe ", mpp_pe(), &
", at point (",i,", ", j, ", ", k, "), a = ", a(i,j,k), ", b = ", b(i,j,k)
call mpp_error(FATAL, trim(string)//': point by point comparison are not OK.')
endif
enddo
enddo
enddo
call mpp_sync()
sum1 = mpp_chksum( a, (/pe/) )
sum2 = mpp_chksum( b, (/pe/) )
if( sum1.EQ.sum2 )then
if( pe.EQ.mpp_root_pe() )call mpp_error( NOTE, trim(string)//': OK.' )
!--- in some case, even though checksum agree, the two arrays
! actually are different, like comparing (1.1,-1.2) with (-1.1,1.2)
!--- hence we need to check the value point by point.
else
call mpp_error( FATAL, trim(string)//': chksums are not OK.' )
end if
end subroutine compare_chksum
end program xgrid_test