!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! 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 !!
!! !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! xgrid_mod - implements exchange grids. An exchange grid is the grid whose
! boundary set is the union of the boundaries of the participating
! grids. The exchange grid is the coarsest grid that is a
! refinement of each of the participating grids. Every exchange
! grid cell is a subarea of one and only one cell in each of the
! participating grids. The exchange grid has two purposes:
!
! (1) The exchange cell areas are used as weights for
! conservative interpolation between model grids.
!
! (2) Computation of surface fluxes takes place on it,
! thereby using the finest scale data obtainable.
!
! The exchange cells are the 2D intersections between cells of the
! participating grids. They are computed elsewhere and are
! read here from a NetCDF grid file as a sequence of quintuples
! (i and j on each of two grids and the cell area).
!
! Each processing element (PE) computes a subdomain of each of the
! participating grids as well as a subset of the exchange cells.
! The geographic regions corresponding to these subdomains will,
! in general, not be the same so communication must occur between
! the PEs. The scheme for doing this is as follows. A distinction
! is drawn between the participating grids. There is a single
! "side 1" grid and it does not have partitions (sub-grid surface
! types). There are one or more "side 2" grids and they may have
! more than 1 partition. In standard usage, the atmosphere grid is
! on side 1 and the land and sea ice grids are on side 2. The set
! of exchange cells computed on a PE corresponds to its side 2
! geographic region(s). Communication between the PEs takes place
! on the side 1 grid. Note: this scheme does not generally allow
! reproduction of answers across varying PE counts. This is
! because, in the side 1 "get", exchange cells are first summed
! locally onto a side 1 grid, then these side 1 contributions are
! further summed after they have been communicated to their target
! PE. For the make_exchange_reproduce option, a special side 1 get
! is used. This get communicates individual exchange cells. The
! cells are summed in the order they appear in the grid spec. file.
! Michael Winton (Michael.Winton@noaa.gov) Oct 2001
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
module xgrid_mod
!
! Michael Winton
!
!
!
! xgrid_mod implements exchange grids for coupled models running on
! multiple processors. An exchange grid is formed from the union of
! the bounding lines of the two (logically rectangular) participating
! grids. The exchange grid is therefore the coarsest grid that is a
! refinement of both participating grids. Exchange grids are used for
! two purposes by coupled models: (1) conservative interpolation of fields
! between models uses the exchange grid cell areas as weights and
! (2) the surface flux calculation takes place on the exchange grid thereby
! using the finest scale data available. xgrid_mod uses a NetCDF grid
! specification file containing the grid cell overlaps in combination with
! the
! mpp_domains domain decomposition information to determine
! the grid and processor connectivities.
!
!
! xgrid_mod is initialized with a list of model identifiers (three characters
! each), a list of mpp_domains domain data structures, and a grid specification
! file name. The first element in the lists refers to the "side one" grid.
! The remaining elements are on "side two". Thus, there may only be a single
! side one grid and it is further restricted to have no partitions (sub-grid
! areal divisions). In standard usage, the atmosphere model is on side one
! and the land and sea ice models are on side two. xgrid_mod performs
! interprocessor communication on the side one grid. Exchange grid variables
! contain no data for zero sized partitions. The size and format of exchange
! grid variables change every time the partition sizes or number of partitions
! are modified with a set_frac_area call on a participating side two grid.
! Existing exchange grid variables cannot be properly interpreted after
! that time; new ones must be allocated and assigned with the put_to_xgrid
! call.
!
!
! The fields of xmap_type are all private.
!
!
! xgrid_mod reads a NetCDF grid specification file to determine the
! grid and processor connectivities. The exchange grids are defined
! by a sequence of quintuples: the i/j indices of the intersecting
! cells of the two participating grids and their areal overlap.
! The names of the five fields are generated automatically from the
! three character ids of the participating grids. For example, if
! the side one grid id is "ATM" and the side two grid id is "OCN",
! xgrid_mod expects to find the following five fields in the grid
! specification file: I_ATM_ATMxOCN, J_ATM_ATMxOCN, I_OCN_ATMxOCN,
! J_OCN_ATMxOCN, and AREA_ATMxOCN. These fields may be generated
! by the make_xgrids utility.
!
use fms_mod, only: file_exist, open_namelist_file, check_nml_error, &
error_mesg, close_file, FATAL, stdlog, &
write_version_number
use mpp_mod, only: mpp_npes, mpp_pe, mpp_root_pe, mpp_send, mpp_recv, &
mpp_sync_self, stdout
use mpp_domains_mod, only: mpp_get_compute_domain, mpp_get_compute_domains, &
Domain2d, mpp_global_sum, mpp_update_domains, &
mpp_modify_domain, mpp_get_data_domain, XUPDATE, &
YUPDATE
use mpp_io_mod, only: mpp_open, MPP_MULTI, MPP_SINGLE, MPP_OVERWR
use constants_mod, only: PI
implicit none
include 'netcdf.inc'
private
public xmap_type, setup_xmap, set_frac_area, put_to_xgrid, get_from_xgrid, &
xgrid_count, some, conservation_check, xgrid_init
!--- paramters that determine the remapping method
integer, parameter :: FIRST_ORDER = 1
integer, parameter :: SECOND_ORDER = 2
integer, parameter :: SECOND_ORDER_MERID = 3
integer, parameter :: SECOND_ORDER_ZONAL = 4
!
!
! Set to .true. to make xgrid_mod reproduce answers on different
! numbers of PEs. This option has a considerable performance impact.
!
!
! exchange grid interpolation method. It has four options:
! "first_order", "second_order", "second_order_merid", "second_order_zonal".
!
logical :: make_exchange_reproduce = .false. ! exactly same on different # PEs
character(len=64) :: interp_method = 'first_order'
namelist /xgrid_nml/ make_exchange_reproduce, interp_method
!
logical :: init = .true.
integer :: remapping_method
!
!
! Scatters data from model grid onto exchange grid.
!
!
! Scatters data from model grid onto exchange grid.
!
!
! call put_to_xgrid(d, grid_id, x, xmap, remap_order)
!
!
!
!
!
!
! exchange grid interpolation method. It has four possible values:
! FIRST_ORDER (=1), SECOND_ORDER(=2), SECOND_ORDER_MERID(=3) and
! SECOND_ORDER_ZONAL(=4). Default value is FIRST_ORDER.
!
interface put_to_xgrid
module procedure put_side1_to_xgrid
module procedure put_side2_to_xgrid
end interface
!
!
!
! Sums data from exchange grid to model grid.
!
!
! Sums data from exchange grid to model grid.
!
!
! call get_from_xgrid(d, grid_id, x, xmap)
!
!
!
!
!
interface get_from_xgrid
module procedure get_side1_from_xgrid
module procedure get_side2_from_xgrid
end interface
!
!
!
! Returns three numbers which are the global sum of a variable.
!
!
! Returns three numbers which are the global sum of a
! variable (1) on its home model grid, (2) after interpolation to the other
! side grid(s), and (3) after re_interpolation back onto its home side grid(s).
! Conservation_check must be called by all PEs to work properly.
!
!
! call conservation_check(d, grid_id, xmap,remap_order)
!
!
!
!
! The global sum of a variable.
!
!
interface conservation_check
module procedure conservation_check_side1
module procedure conservation_check_side2
end interface
!
type xcell_type
integer :: i1, j1, i2, j2 ! indices of cell in model arrays on both sides
integer :: pe ! other side pe that has this cell
real :: area ! geographic area of exchange cell
real :: di, dj ! Weight for the gradient of flux
end type xcell_type
type grid_type
character(len=3) :: id ! grid identifier
integer, pointer, dimension(:) :: is =>NULL(), ie =>NULL() ! domain - i-range (pe index)
integer, pointer, dimension(:) :: js =>NULL(), je =>NULL() ! domain - j-range (pe index)
integer, pointer :: is_me =>NULL(), ie_me =>NULL() ! my domain - i-range
integer, pointer :: js_me =>NULL(), je_me =>NULL() ! my domain - j-range
integer :: isd_me, ied_me ! my data domain - i-range
integer :: jsd_me, jed_me ! my data domain - j-range
integer :: im , jm , km ! global domain range
real, pointer, dimension(:) :: lon =>NULL(), lat =>NULL() ! center of global grids
real, pointer, dimension(:,:,:) :: frac_area =>NULL() ! partition fractions
real, pointer, dimension(:,:) :: area =>NULL() ! cell area
real, pointer, dimension(:,:) :: area_inv =>NULL() ! 1 / area for normalization
integer :: first, last ! xgrid index range
integer :: size ! # xcell patterns
type(xcell_type), pointer, dimension(:) :: x =>NULL() ! xcell patterns
integer :: size_repro ! # side 1 patterns for repro
type(xcell_type), pointer, dimension(:) :: x_repro =>NULL() ! side 1 patterns for repro
type(Domain2d) :: domain ! used for conservation checks
type(Domain2d) :: domain_with_halo ! used for second order remapping
end type grid_type
type x1_type
integer :: i, j
real :: area ! (= geographic area * frac_area)
real :: di, dj ! weight for the gradient of flux
end type x1_type
type x2_type
integer :: i, j, k
real :: area ! geographic area of exchange cell
end type x2_type
type xmap_type
private
integer :: size ! # of exchange grid cells with area > 0 on this pe
integer :: me, npes, root_pe
logical, pointer, dimension(:) :: your1my2 =>NULL()! true if side 1 domain on
! indexed pe overlaps side 2
! domain on this pe
logical, pointer, dimension(:) :: your2my1 =>NULL() ! true if a side 2 domain on
! indexed pe overlaps side 1
! domain on this pe
type (grid_type), pointer, dimension(:) :: grids =>NULL() ! 1st grid is side 1;
! rest on side 2
!
! Description of the individual exchange grid cells (index is cell #)
!
type(x1_type), pointer, dimension(:) :: x1 =>NULL() ! side 1 info
type(x2_type), pointer, dimension(:) :: x2 =>NULL() ! side 2 info
real, pointer, dimension(:) :: send_buffer =>NULL() ! for non-blocking sends
real, pointer, dimension(:) :: recv_buffer =>NULL() ! for non-blocking recv
integer, pointer, dimension(:) :: send_count_repro =>NULL(), recv_count_repro =>NULL()
end type xmap_type
!-----------------------------------------------------------------------
character(len=128) :: version = '$Id$'
character(len=128) :: tagname = '$Name$'
real, parameter :: EPS = 1.0e-10
logical :: module_is_initialized = .FALSE.
contains
!#######################################################################
logical function in_box(i, j, is, ie, js, je)
integer :: i, j, is, ie, js, je
in_box = (i>=is) .and. (i<=ie) .and. (j>=js) .and. (j<=je)
end function in_box
!#######################################################################
!
!
! Initialize the xgrid_mod.
!
!
! Initialization routine for the xgrid module. It reads the xgrid_nml,
! writes the version information and xgrid_nml to the log file.
!
!
! call xgrid_init ( )
!
!
! exchange grid interpolation method. It has four possible values:
! FIRST_ORDER (=1), SECOND_ORDER(=2), SECOND_ORDER_MERID(=3) and
! SECOND_ORDER_ZONAL(=4)
!
subroutine xgrid_init(remap_method)
integer, intent(out) :: remap_method
integer :: unit, ierr, io
if (module_is_initialized) return
module_is_initialized = .TRUE.
if ( file_exist( 'input.nml' ) ) then
unit = open_namelist_file ( )
ierr = 1
do while ( ierr /= 0 )
read ( unit, nml = xgrid_nml, iostat = io, end = 10 )
ierr = check_nml_error ( io, 'xgrid_nml' )
enddo
10 continue
call close_file ( unit )
endif
!--------- write version number and namelist ------------------
call write_version_number (version, tagname)
unit = stdlog ( )
if ( mpp_pe() == mpp_root_pe() ) write (unit,nml=xgrid_nml)
call close_file (unit)
!--------- check interp_method has suitable value
select case(trim(interp_method))
case('first_order')
remap_method = FIRST_ORDER
case('second_order')
remap_method = SECOND_ORDER
case('second_order_merid')
remap_method = SECOND_ORDER_MERID
case('second_order_zonal')
remap_method = SECOND_ORDER_ZONAL
case default
call error_mesg('xgrid_mod', ' nml interp_method = ' //trim(interp_method)// &
' is not a valid namelist option', FATAL)
end select
remapping_method = remap_method
end subroutine xgrid_init
!
!#######################################################################
subroutine load_xgrid (xmap, grid, domain, ncid, id_i1, id_j1, id_i2, id_j2, &
id_area, n_areas, id_di, id_dj )
type(xmap_type), intent(inout) :: xmap
type(grid_type), intent(inout) :: grid
type(Domain2d), intent(inout) :: domain
integer, intent(in) :: ncid, id_i1, id_j1, id_i2, id_j2, id_area, n_areas
integer, intent(in), optional :: id_di, id_dj
integer, dimension(n_areas) :: i1, j1, i2, j2 ! xgrid quintuples
real, dimension(n_areas) :: area, di, dj ! from grid file
type (grid_type), pointer, save :: grid1 =>NULL()
integer, dimension(0:xmap%npes-1) :: is_2, ie_2, js_2, je_2 ! side 2 decomp.
integer :: start(4), nread(4), rcode, l, ll, ll_repro, p
grid1 => xmap%grids(1)
start = 1; nread = 1; nread(1) = n_areas
rcode = nf_get_vara_int(ncid, id_i1, start, nread, i1)
rcode = nf_get_vara_int(ncid, id_j1, start, nread, j1)
rcode = nf_get_vara_int(ncid, id_i2, start, nread, i2)
rcode = nf_get_vara_int(ncid, id_j2, start, nread, j2)
rcode = nf_get_vara_double(ncid, id_area, start, nread, area)
di = 0.0; dj = 0.0
if(present(id_di)) rcode = nf_get_vara_double(ncid, id_di, start, nread, di)
if(present(id_dj)) rcode = nf_get_vara_double(ncid, id_dj, start, nread, dj)
do l=1,n_areas
if (in_box(i1(l), j1(l), grid1%is_me, grid1%ie_me, &
grid1%js_me, grid1%je_me) ) then
grid1%area(i1(l),j1(l)) = grid1%area(i1(l),j1(l))+area(l)
do p=0,xmap%npes-1
if (in_box(i2(l), j2(l), grid%is(p), grid%ie(p), &
grid%js(p), grid%je(p))) then
xmap%your2my1(p) = .true.
end if
end do
grid%size_repro = grid%size_repro + 1
end if
if (in_box(i2(l), j2(l), grid%is_me, grid%ie_me, &
grid%js_me, grid%je_me) ) then
grid%size = grid%size + 1
grid%area(i2(l),j2(l)) = grid%area(i2(l),j2(l))+area(l)
do p=0,xmap%npes-1
if (in_box(i1(l), j1(l), grid1%is(p), grid1%ie(p), &
grid1%js(p), grid1%je(p))) then
xmap%your1my2(p) = .true.
end if
end do
end if
end do
allocate( grid%x( grid%size ) )
if (make_exchange_reproduce) allocate ( grid%x_repro(grid%size_repro) )
ll = 0
ll_repro = 0
do l=1,n_areas
if (in_box(i2(l), j2(l), grid%is_me, grid%ie_me, grid%js_me, grid%je_me)) then
! insert in this grids cell pattern list and add area to side 2 area
ll = ll + 1
grid%x(ll)%i1 = i1(l); grid%x(ll)%i2 = i2(l)
grid%x(ll)%j1 = j1(l); grid%x(ll)%j2 = j2(l)
grid%x(ll)%area = area(l)
grid%x(ll)%di = di(l)
grid%x(ll)%dj = dj(l)
if (make_exchange_reproduce) then
do p=0,xmap%npes-1
if (in_box(i1(l), j1(l), grid1%is(p), grid1%ie(p), &
grid1%js(p), grid1%je(p))) then
grid%x(ll)%pe = p + xmap%root_pe
end if
end do
end if ! make_exchange reproduce
end if
if (in_box(i1(l),j1(l), grid1%is_me,grid1%ie_me, grid1%js_me,grid1%je_me) &
.and. make_exchange_reproduce ) then
ll_repro = ll_repro + 1
grid%x_repro(ll_repro)%i1 = i1(l); grid%x_repro(ll_repro)%i2 = i2(l)
grid%x_repro(ll_repro)%j1 = j1(l); grid%x_repro(ll_repro)%j2 = j2(l)
grid%x_repro(ll_repro)%area = area(l)
grid%x_repro(ll_repro)%di = di(l)
grid%x_repro(ll_repro)%dj = dj(l)
do p=0,xmap%npes-1
if (in_box(i2(l), j2(l), grid%is(p), grid%ie(p), &
grid%js(p), grid%je(p))) then
grid%x_repro(ll_repro)%pe = p + xmap%root_pe
end if
end do
end if ! make_exchange_reproduce
end do
grid%area_inv = 0.0;
where (grid%area>0.0) grid%area_inv = 1.0/grid%area
end subroutine load_xgrid
!#######################################################################
!
! get_grid - read the center point of the grid from grid_spec.nc.
! - only the grid at the side 1 is needed, so we only read
! - atm and land grid
!
!
subroutine get_grid(grid, grid_id, ncid)
type(grid_type), intent(inout) :: grid
integer, intent(in) :: ncid
character(len=3), intent(in) :: grid_id
integer :: rcode, start(4), nread(4), id_lon, id_lat, i, j
real, dimension(grid%im) :: lonb
real, dimension(grid%jm) :: latb
real :: d2r
d2r = PI/180.0
start = 1; nread = 1
if(grid_id == 'ATM') then
nread(1) = grid%im
rcode = nf_inq_varid(ncid,'xta',id_lon)
if (rcode/=0) call error_mesg('xgrid_mod', 'cannot find grid file field xta', FATAL)
rcode = nf_get_vara_double(ncid, id_lon, start, nread,lonb)
rcode = nf_inq_varid(ncid,'yta',id_lat)
if (rcode/=0) call error_mesg('xgrid_mod', 'cannot find grid file field yta', FATAL)
nread(1) = grid%jm
rcode = nf_get_vara_double(ncid, id_lat, start, nread,latb)
else if(grid_id == 'LND') then
nread(1) = grid%im
rcode = nf_inq_varid(ncid,'xtl',id_lon)
if (rcode/=0) call error_mesg('xgrid_mod', 'cannot find grid file field xtl', FATAL)
rcode = nf_get_vara_double(ncid, id_lon, start, nread,lonb)
nread(1) = grid%jm
rcode = nf_inq_varid(ncid,'ytl',id_lat)
if (rcode/=0) call error_mesg('xgrid_mod', 'cannot find grid file field ytl', FATAL)
rcode = nf_get_vara_double(ncid, id_lat, start, nread,latb)
endif
!--- second order remapping suppose second order
if(grid_id == 'LND' .or. grid_id == 'ATM') then
grid%lon = lonb * d2r
grid%lat = latb * d2r
endif
return
end subroutine get_grid
!#######################################################################
!
!
! Sets up exchange grid connectivity using grid specification file and
! processor domain decomposition.
!
!
! Sets up exchange grid connectivity using grid specification file and
! processor domain decomposition. Initializes xmap.
!
!
! call setup_xmap(xmap, grid_ids, grid_domains, grid_file)
!
!
!
!
!
subroutine setup_xmap(xmap, grid_ids, grid_domains, grid_file )
type (xmap_type), intent(inout) :: xmap
character(len=3), dimension(:), intent(in ) :: grid_ids
!!$ type(Domain2d), dimension(size(grid_ids(:))), intent(in ) :: grid_domains
type(Domain2d), dimension(:), intent(in ) :: grid_domains
character(len=*) , intent(in ) :: grid_file
integer :: g, l, ll, p, n_areas, send_size, recv_size
integer :: ncid, i1_id, j1_id, i2_id, j2_id, area_id, di_id, dj_id
integer :: dims(4), rcode
integer :: unit
type (grid_type), pointer, save :: grid =>NULL(), grid1 =>NULL()
real, dimension(3) :: xxx
real, dimension(:,:), allocatable :: check_data
real, dimension(:,:,:), allocatable :: check_data_3D
integer :: i, j
logical :: use_higher_order = .false.
if(interp_method .ne. 'first_order') use_higher_order = .true.
xmap%me = mpp_pe ()
xmap%npes = mpp_npes()
xmap%root_pe = mpp_root_pe()
allocate( xmap%grids(1:size(grid_ids(:))) )
allocate ( xmap%your1my2(0:xmap%npes-1), xmap%your2my1(0:xmap%npes-1) )
xmap%your1my2 = .false.; xmap%your2my1 = .false.;
rcode = nf_open(grid_file,0,ncid)
if (rcode/=0) call error_mesg ('xgrid_mod', 'cannot open grid file', FATAL)
do g=1,size(grid_ids(:))
grid => xmap%grids(g)
if (g==1) grid1 => xmap%grids(g)
grid%id = grid_ids (g)
grid%domain = grid_domains(g)
call mpp_modify_domain(grid%domain, grid%domain_with_halo, xhalo = 1, yhalo =1)
call mpp_get_data_domain(grid%domain_with_halo, grid%isd_me, grid%ied_me, grid%jsd_me, grid%jed_me)
allocate ( grid%is(0:xmap%npes-1), grid%ie(0:xmap%npes-1) )
allocate ( grid%js(0:xmap%npes-1), grid%je(0:xmap%npes-1) )
call mpp_get_compute_domains(grid%domain, xbegin=grid%is, xend=grid%ie, &
ybegin=grid%js, yend=grid%je )
grid%is_me => grid%is(xmap%me-xmap%root_pe); grid%ie_me => grid%ie(xmap%me-xmap%root_pe)
grid%js_me => grid%js(xmap%me-xmap%root_pe); grid%je_me => grid%je(xmap%me-xmap%root_pe)
grid%im = maxval(grid%ie)
grid%jm = maxval(grid%je)
grid%km = 1
allocate( grid%area (grid%is_me:grid%ie_me, grid%js_me:grid%je_me) )
allocate( grid%area_inv(grid%is_me:grid%ie_me, grid%js_me:grid%je_me) )
grid%area = 0.0
grid%size = 0
grid%size_repro = 0
if (g>1) then
allocate( grid%frac_area(grid%is_me:grid%ie_me, grid%js_me:grid%je_me, &
grid%km ) )
grid%frac_area = 1.0
rcode = nf_inq_varid(ncid, &
'I_'//grid_ids(1)//'_'//grid_ids(1)//'x'//grid_ids(g),&
i1_id)
if (rcode/=0) &
call error_mesg('xgrid_mod', 'cannot find grid file field', FATAL)
rcode = nf_inq_varid(ncid, &
'J_'//grid_ids(1)//'_'//grid_ids(1)//'x'//grid_ids(g),&
j1_id)
if (rcode/=0) &
call error_mesg('xgrid_mod', 'cannot find grid file field', FATAL)
rcode = nf_inq_varid(ncid, &
'I_'//grid_ids(g)//'_'//grid_ids(1)//'x'//grid_ids(g),&
i2_id)
if (rcode/=0) &
call error_mesg('xgrid_mod', 'cannot find grid file field', FATAL)
rcode = nf_inq_varid(ncid, &
'J_'//grid_ids(g)//'_'//grid_ids(1)//'x'//grid_ids(g),&
j2_id)
if (rcode/=0) &
call error_mesg('xgrid_mod', 'cannot find grid file field', FATAL)
rcode = nf_inq_varid(ncid, 'AREA_'//grid_ids(1)//'x'//grid_ids(g), area_id)
if (rcode/=0) &
call error_mesg('xgrid_mod', 'cannot find grid file field', FATAL)
if(use_higher_order) then
rcode = nf_inq_varid(ncid, 'DI_'//grid_ids(1)//'x'//grid_ids(g), di_id)
if (rcode/=0) &
call error_mesg('xgrid_mod', 'cannot find grid file field DI1_'//grid_ids(1)//'x'//grid_ids(g), FATAL)
rcode = nf_inq_varid(ncid, 'DJ_'//grid_ids(1)//'x'//grid_ids(g), dj_id)
if (rcode/=0) &
call error_mesg('xgrid_mod', 'cannot find grid file field DJ1'//grid_ids(1)//'x'//grid_ids(g), FATAL)
endif
rcode = nf_inq_vardimid(ncid, area_id, dims)
rcode = nf_inq_dimlen(ncid, dims(1), n_areas)
! load exchange cells, sum grid cell areas, set your1my2/your2my1
if(n_areas > 0) then
if(use_higher_order) then
call load_xgrid (xmap, grid, grid%domain, ncid, i1_id, j1_id, i2_id, j2_id, &
area_id, n_areas, di_id, dj_id)
else
call load_xgrid (xmap, grid, grid%domain, ncid, i1_id, j1_id, i2_id, j2_id, &
area_id, n_areas )
endif
endif
end if
! get the center point of the grid box
allocate(grid%lon(grid%im), grid%lat(grid%jm))
call get_grid(grid, grid_ids(g), ncid)
end do
rcode = nf_close(ncid)
grid1%area_inv = 0.0;
where (grid1%area>0.0)
grid1%area_inv = 1.0/grid1%area
end where
xmap%your1my2(xmap%me-xmap%root_pe) = .false. ! this is not necessarily true but keeps
xmap%your2my1(xmap%me-xmap%root_pe) = .false. ! a PE from communicating with itself
send_size = grid1%im*grid1%jm
recv_size = maxval((grid1%ie-grid1%is+1)*(grid1%je-grid1%js+1) )
if (make_exchange_reproduce) then
allocate( xmap%send_count_repro(0:xmap%npes-1) )
allocate( xmap%recv_count_repro(0:xmap%npes-1) )
xmap%send_count_repro = 0
xmap%recv_count_repro = 0
do g=2,size(xmap%grids(:))
do p=0,xmap%npes-1
xmap%send_count_repro(p) = xmap%send_count_repro(p) &
+count(xmap%grids(g)%x (:)%pe==p+xmap%root_pe)
xmap%recv_count_repro(p) = xmap%recv_count_repro(p) &
+count(xmap%grids(g)%x_repro(:)%pe==p+xmap%root_pe)
end do
end do
send_size = max(send_size, sum(xmap%send_count_repro))
end if
allocate (xmap%send_buffer(send_size))
allocate (xmap%recv_buffer(recv_size))
call mpp_open( unit, 'xgrid.out', action=MPP_OVERWR, threading=MPP_MULTI, &
fileset=MPP_SINGLE, nohdrs=.TRUE. )
write( unit,* )xmap%grids(:)%id, ' GRID: PE ', xmap%me, ' #XCELLS=', &
xmap%grids(2:size(xmap%grids(:)))%size, ' #COMM. PARTNERS=', &
count(xmap%your1my2), '/', count(xmap%your2my1), &
pack((/(p+xmap%root_pe,p=0,xmap%npes-1)/), xmap%your1my2), &
'/', pack((/(p+xmap%root_pe,p=0,xmap%npes-1)/), xmap%your2my1)
allocate( xmap%x1(1:sum(xmap%grids(2:size(xmap%grids(:)))%size)) )
allocate( xmap%x2(1:sum(xmap%grids(2:size(xmap%grids(:)))%size)) )
call regen(xmap)
xxx = conservation_check(grid1%area*0+1.0, grid1%id, xmap)
write(stdout(),* )"Checked data is array of constant 1"
write(stdout(),* )grid1%id,'(',xmap%grids(:)%id,')=', xxx
do g=2,size(xmap%grids(:))
xxx = conservation_check(xmap%grids(g)%frac_area*0+1.0, xmap%grids(g)%id, xmap )
write( stdout(),* )xmap%grids(g)%id,'(',xmap%grids(:)%id,')=', xxx
enddo
! create an random number 2d array
if(grid1%id == "ATM") then
allocate(check_data(size(grid1%area,1), size(grid1%area,2)))
call random_number(check_data)
!--- second order along both zonal and meridinal direction
xxx = conservation_check(check_data, grid1%id, xmap, remap_method = remapping_method )
write( stdout(),* ) &
"Checked data is array of random number between 0 and 1 using "//trim(interp_method)
write( stdout(),* )grid1%id,'(',xmap%grids(:)%id,')=', xxx
deallocate(check_data)
do g=2,size(xmap%grids(:))
allocate(check_data_3d(size(xmap%grids(g)%frac_area,1),size(xmap%grids(g)%frac_area,2), &
size(xmap%grids(g)%frac_area,3) ))
call random_number(check_data_3d)
xxx = conservation_check(check_data_3d, xmap%grids(g)%id, xmap, remap_method = remapping_method )
write( stdout(),* )xmap%grids(g)%id,'(',xmap%grids(:)%id,')=', xxx
deallocate( check_data_3d)
end do
endif
call close_file (unit)
end subroutine setup_xmap
!
!#######################################################################
subroutine regen(xmap)
type (xmap_type), intent(inout) :: xmap
integer :: g, l, i, j, k, max_size
max_size = 0;
do g=2,size(xmap%grids(:))
max_size = max_size + xmap%grids(g)%size * xmap%grids(g)%km
end do
if (max_size>size(xmap%x1(:))) then
deallocate(xmap%x1)
deallocate(xmap%x2)
allocate( xmap%x1(1:max_size) )
allocate( xmap%x2(1:max_size) )
end if
xmap%size = 0
do g=2,size(xmap%grids(:))
xmap%grids(g)%first = xmap%size + 1;
do l=1,xmap%grids(g)%size
i = xmap%grids(g)%x(l)%i2
j = xmap%grids(g)%x(l)%j2
do k=1,xmap%grids(g)%km
if (xmap%grids(g)%frac_area(i,j,k)/=0.0) then
xmap%size = xmap%size+1
xmap%x1(xmap%size)%i = xmap%grids(g)%x(l)%i1
xmap%x1(xmap%size)%j = xmap%grids(g)%x(l)%j1
xmap%x1(xmap%size)%area = xmap%grids(g)%x(l)%area &
*xmap%grids(g)%frac_area(i,j,k)
xmap%x1(xmap%size)%di = xmap%grids(g)%x(l)%di
xmap%x1(xmap%size)%dj = xmap%grids(g)%x(l)%dj
xmap%x2(xmap%size)%i = xmap%grids(g)%x(l)%i2
xmap%x2(xmap%size)%j = xmap%grids(g)%x(l)%j2
xmap%x2(xmap%size)%k = k
xmap%x2(xmap%size)%area = xmap%grids(g)%x(l)%area
end if
end do
end do
xmap%grids(g)%last = xmap%size
end do
end subroutine regen
!#######################################################################
!
!
! Changes sub-grid portion areas and/or number.
!
!
! Changes sub-grid portion areas and/or number.
!
!
! call set_frac_area(f, grid_id, xmap)
!
!
!
!
subroutine set_frac_area(f, grid_id, xmap)
real, dimension(:,:,:), intent(in ) :: f
character(len=3), intent(in ) :: grid_id
type (xmap_type), intent(inout) :: xmap
integer :: g
type(grid_type), pointer, save :: grid =>NULL()
if (grid_id==xmap%grids(1)%id) call error_mesg ('xgrid_mod', &
'set_frac_area called on side 1 grid', FATAL)
do g=2,size(xmap%grids(:))
grid => xmap%grids(g)
if (grid_id==grid%id) then
if (size(f,3)/=size(grid%frac_area,3)) then
deallocate (grid%frac_area)
grid%km = size(f,3);
allocate( grid%frac_area(grid%is_me:grid%ie_me, grid%js_me:grid%je_me, &
grid%km) )
end if
grid%frac_area = f;
call regen(xmap)
return;
end if
end do
call error_mesg ('xgrid_mod', 'set_frac_area: could not find grid id', FATAL)
end subroutine set_frac_area
!
!#######################################################################
!
!
! Returns current size of exchange grid variables.
!
!
! Returns current size of exchange grid variables.
!
!
! xgrid_count(xmap)
!
!
!
integer function xgrid_count(xmap)
type (xmap_type), intent(inout) :: xmap
xgrid_count = xmap%size
end function xgrid_count
!
!#######################################################################
!
!
!
!
!
!
subroutine put_side1_to_xgrid(d, grid_id, x, xmap, remap_method)
real, dimension(:,:), intent(in ) :: d
character(len=3), intent(in ) :: grid_id
real, dimension(:), intent(inout) :: x
type (xmap_type), intent(inout) :: xmap
integer, intent(in), optional :: remap_method
integer :: g, method
method = FIRST_ORDER ! default
if(present(remap_method)) method = remap_method
if (grid_id==xmap%grids(1)%id) then
if(method == FIRST_ORDER) then
call put_1_to_xgrid_order_1(d, x, xmap)
else
call put_1_to_xgrid_order_2(d, x, xmap, method )
endif
return;
end if
do g=2,size(xmap%grids(:))
if (grid_id==xmap%grids(g)%id) &
call error_mesg ('xgrid_mod', &
'put_to_xgrid expects a 3D side 2 grid', FATAL)
end do
call error_mesg ('xgrid_mod', 'put_to_xgrid: could not find grid id', FATAL)
end subroutine put_side1_to_xgrid
!
!#######################################################################
!
!
!
!
!
subroutine put_side2_to_xgrid(d, grid_id, x, xmap)
real, dimension(:,:,:), intent(in ) :: d
character(len=3), intent(in ) :: grid_id
real, dimension(:), intent(inout) :: x
type (xmap_type), intent(inout) :: xmap
integer :: g
if (grid_id==xmap%grids(1)%id) &
call error_mesg ('xgrid_mod', &
'put_to_xgrid expects a 2D side 1 grid', FATAL)
do g=2,size(xmap%grids(:))
if (grid_id==xmap%grids(g)%id) then
call put_2_to_xgrid(d, xmap%grids(g), x, xmap)
return;
end if
end do
call error_mesg ('xgrid_mod', 'put_to_xgrid: could not find grid id', FATAL)
end subroutine put_side2_to_xgrid
!
!#######################################################################
!
!
!
!
!
subroutine get_side1_from_xgrid(d, grid_id, x, xmap)
real, dimension(:,:), intent( out) :: d
character(len=3), intent(in ) :: grid_id
real, dimension(:), intent(in ) :: x
type (xmap_type), intent(inout) :: xmap
integer :: g
if (grid_id==xmap%grids(1)%id) then
if (make_exchange_reproduce) then
call get_1_from_xgrid_repro(d, x, xmap)
else
call get_1_from_xgrid(d, x, xmap)
end if
return;
end if
do g=2,size(xmap%grids(:))
if (grid_id==xmap%grids(g)%id) &
call error_mesg ('xgrid_mod', &
'get_from_xgrid expects a 3D side 2 grid', FATAL)
end do
call error_mesg ('xgrid_mod', 'get_from_xgrid: could not find grid id', FATAL)
end subroutine get_side1_from_xgrid
!
!#######################################################################
!
!
!
!
!
subroutine get_side2_from_xgrid(d, grid_id, x, xmap)
real, dimension(:,:,:), intent( out) :: d
character(len=3), intent(in ) :: grid_id
real, dimension(:), intent(in ) :: x
type (xmap_type), intent(in ) :: xmap
integer :: g
if (grid_id==xmap%grids(1)%id) &
call error_mesg ('xgrid_mod', &
'get_from_xgrid expects a 2D side 1 grid', FATAL)
do g=2,size(xmap%grids(:))
if (grid_id==xmap%grids(g)%id) then
call get_2_from_xgrid(d, xmap%grids(g), x, xmap)
return;
end if
end do
call error_mesg ('xgrid_mod', 'get_from_xgrid: could not find grid id', FATAL)
end subroutine get_side2_from_xgrid
!
!#######################################################################
!
!
! Returns logical associating exchange grid cells with given side two grid.
!
!
! Returns logical associating exchange grid cells with given side two grid.
!
!
! call some(xmap, some_arr, grid_id)
!
!
!
!
! logical associating exchange grid cells with given side 2 grid.
!
subroutine some(xmap, some_arr, grid_id)
type (xmap_type), intent(in) :: xmap
character(len=3), optional, intent(in) :: grid_id
logical, dimension(xmap%size), intent(out) :: some_arr
integer :: g, l
if (.not.present(grid_id)) then
some_arr = .true.
return;
end if
if (grid_id==xmap%grids(1)%id) &
call error_mesg ('xgrid_mod', 'some expects a side 2 grid id', FATAL)
do g=2,size(xmap%grids(:))
if (grid_id==xmap%grids(g)%id) then
some_arr = .false.
some_arr(xmap%grids(g)%first:xmap%grids(g)%last) = .true.;
return;
end if
end do
call error_mesg ('xgrid_mod', 'some could not find grid id', FATAL)
end subroutine some
!
!#######################################################################
subroutine put_2_to_xgrid(d, grid, x, xmap)
type (grid_type), intent(in) :: grid
real, dimension(grid%is_me:grid%ie_me, &
grid%js_me:grid%je_me, grid%km), intent(in) :: d
real, dimension(: ), intent(inout) :: x
type (xmap_type), intent(in ) :: xmap
integer :: k, l
do l=grid%first,grid%last
x(l) = d(xmap%x2(l)%i,xmap%x2(l)%j,xmap%x2(l)%k)
end do
end subroutine put_2_to_xgrid
!#######################################################################
subroutine get_2_from_xgrid(d, grid, x, xmap)
type (grid_type), intent(in ) :: grid
real, dimension(grid%is_me:grid%ie_me, &
grid%js_me:grid%je_me, grid%km), intent(out) :: d
real, dimension(:), intent(in ) :: x
type (xmap_type), intent(in ) :: xmap
integer :: l, k
d = 0.0
do l=grid%first,grid%last
d(xmap%x2(l)%i,xmap%x2(l)%j,xmap%x2(l)%k) = &
d(xmap%x2(l)%i,xmap%x2(l)%j,xmap%x2(l)%k) + xmap%x2(l)%area*x(l)
end do
!
! normalize with side 2 grid cell areas
!
do k=1,size(d,3)
d(:,:,k) = d(:,:,k) * grid%area_inv
end do
end subroutine get_2_from_xgrid
!#######################################################################
function get_side_1(pe, im, jm)
integer, intent(in) :: pe, im, jm
real, dimension(im,jm) :: get_side_1
real, dimension(im*jm) :: buf
integer :: i, j, l
! call mpp_recv(buf, im*jm, pe)
! l = 0
! do j=1,jm; do i=1,im;
! l = l + 1
! get_side_1(i,j) = buf(l)
! end do; end do
! Force use of "scalar", integer pointer mpp interface.
call mpp_recv( get_side_1(1,1), glen=im*jm, from_pe=pe )
end function get_side_1
!#######################################################################
subroutine put_1_to_xgrid_order_1(d, x, xmap)
real, dimension(:,:), intent(in ) :: d
real, dimension(: ), intent(inout) :: x
type (xmap_type), intent(inout) :: xmap
integer :: i, is, ie, im, j, js, je, jm, p, l
real, dimension(xmap%grids(1)%im,xmap%grids(1)%jm) :: dg
type (grid_type), pointer, save :: grid1 =>NULL()
grid1 => xmap%grids(1)
is = grid1%is_me; ie = grid1%ie_me;
js = grid1%js_me; je = grid1%je_me;
dg(is:ie,js:je) = d;
im = ie-is+1; jm = je-js+1;
l = 0
call mpp_sync_self() !Balaji
do j=1,jm; do i=1,im;
l = l + 1;
xmap%send_buffer(l) = d(i,j)
end do; end do;
do p=0,xmap%npes-1
if (xmap%your2my1(p)) then
! Force use of "scalar", integer pointer mpp interface.
call mpp_send(xmap%send_buffer(1), plen=im*jm, to_pe=p+xmap%root_pe);
end if
end do
do p=0,xmap%npes-1
if (xmap%your1my2(p)) then
is = grid1%is(p); ie = grid1%ie(p);
js = grid1%js(p); je = grid1%je(p);
dg(is:ie,js:je) = get_side_1(p+xmap%root_pe,ie-is+1,je-js+1);
end if
end do
do l=1,xmap%size
x(l) = dg(xmap%x1(l)%i,xmap%x1(l)%j)
end do
! call mpp_sync_self
end subroutine put_1_to_xgrid_order_1
!#######################################################################
subroutine put_1_to_xgrid_order_2(d, x, xmap, remap_method)
real, dimension(:,:), intent(in ) :: d
real, dimension(: ), intent(inout) :: x
type (xmap_type), intent(inout) :: xmap
integer, intent(in) :: remap_method
integer :: i, is, ie, im, j, js, je, jm, p, l, isd, ied, jsd, jed
real, dimension(xmap%grids(1)%im,xmap%grids(1)%jm) :: dg
real, dimension(xmap%grids(1)%im,xmap%grids(1)%jm) :: grad_x, grad_y
real, dimension(xmap%grids(1)%isd_me:xmap%grids(1)%ied_me,xmap%grids(1)%jsd_me:xmap%grids(1)%jed_me) :: tmp
type (grid_type), pointer, save :: grid1 =>NULL()
integer :: num_block, send_size, recv_size
grid1 => xmap%grids(1)
is = grid1%is_me; ie = grid1%ie_me
js = grid1%js_me; je = grid1%je_me
isd = grid1%isd_me; ied = grid1%ied_me
jsd = grid1%jsd_me; jed = grid1%jed_me
im = ie-is+1; jm = je-js+1
dg(is:ie,js:je) = d
! first get the halo of data
tmp(is:ie,js:je) = d(:,:)
if (remap_method == SECOND_ORDER_ZONAL) then
call mpp_update_domains(tmp,grid1%domain_with_halo, flags=XUPDATE)
else if ( remap_method == SECOND_ORDER_MERID ) then
call mpp_update_domains(tmp,grid1%domain_with_halo, flags=YUPDATE)
else
call mpp_update_domains(tmp,grid1%domain_with_halo)
endif
num_block = 1
if( remap_method .ne. SECOND_ORDER_ZONAL ) then ! second_order_merid or second_order
grad_y(is:ie,js:je) = grad_merid(tmp, grid1%lat, is, ie, js, je,isd, ied, jsd, jed)
num_block = num_block + 1
endif
if (remap_method .ne. SECOND_ORDER_MERID) then ! second_order_zonal or second_order
grad_x(is:ie,js:je) = grad_zonal(tmp, grid1%lon, grid1%lat, is, ie, js, je, isd, ied, jsd, jed)
num_block = num_block + 1
endif
call mpp_sync_self() !Balaji
send_size = num_block*im*jm
! if size of send_buffer is not enough, need to reallocate send_buffer
if(size(xmap%send_buffer(:)) .lt. send_size) then
deallocate(xmap%send_buffer)
allocate(xmap%send_buffer(send_size))
endif
l = 0
do j=js,je; do i=is,ie
l = l + 1
xmap%send_buffer(l) = tmp(i,j)
end do; end do
if(remap_method .ne. SECOND_ORDER_ZONAL) then
do j=js,je; do i=is,ie
l = l + 1
xmap%send_buffer(l) = grad_y(i,j)
end do; end do
endif
if (remap_method .ne. SECOND_ORDER_MERID) then
do j=js,je; do i=is,ie
l = l + 1
xmap%send_buffer(l) = grad_x(i,j)
end do; end do
endif
do p=0,xmap%npes-1
if (xmap%your2my1(p)) then
! Force use of "scalar", integer pointer mpp interface.
call mpp_send(xmap%send_buffer(1), plen=send_size, to_pe=p+xmap%root_pe);
end if
end do
do p=0,xmap%npes-1
if (xmap%your1my2(p)) then
is = grid1%is(p); ie = grid1%ie(p)
js = grid1%js(p); je = grid1%je(p)
recv_size = num_block*(ie-is+1)*(je-js+1)
if(size(xmap%recv_buffer(:)) .lt. recv_size) then
deallocate(xmap%recv_buffer)
allocate(xmap%recv_buffer(recv_size))
endif
call mpp_recv(xmap%recv_buffer(1), glen = recv_size, from_pe = p+xmap%root_pe)
l = 0
do j = js,je; do i=is,ie
l = l + 1
dg(i,j) = xmap%recv_buffer(l)
enddo; enddo
if(remap_method .ne. SECOND_ORDER_ZONAL) then
do j = js,je; do i=is,ie
l = l + 1
grad_y(i,j) = xmap%recv_buffer(l)
enddo; enddo
endif
if (remap_method .ne. SECOND_ORDER_MERID) then
do j = js,je; do i=is,ie
l = l + 1
grad_x(i,j) = xmap%recv_buffer(l)
enddo; enddo
endif
end if
end do
do l=1,xmap%size
x(l) = dg(xmap%x1(l)%i,xmap%x1(l)%j)
if(remap_method .ne. SECOND_ORDER_ZONAL) then
x(l) = x(l) + grad_y(xmap%x1(l)%i,xmap%x1(l)%j ) *xmap%x1(l)%dj
endif
if (remap_method .ne. SECOND_ORDER_MERID) then
x(l) = x(l) + grad_x(xmap%x1(l)%i,xmap%x1(l)%j ) *xmap%x1(l)%di
endif
end do
end subroutine put_1_to_xgrid_order_2
!#######################################################################
subroutine get_1_from_xgrid(d, x, xmap)
real, dimension(:,:), intent(out) :: d
real, dimension(: ), intent(in ) :: x
type (xmap_type), intent(inout) :: xmap
real, dimension(xmap%grids(1)%im,xmap%grids(1)%jm), target :: dg
integer :: i, is, ie, im, j, js, je, jm, l, le, p
real , pointer, save :: dgp =>NULL()
type (grid_type) , pointer, save :: grid1 =>NULL()
grid1 => xmap%grids(1)
dg = 0.0;
do l=1,xmap%size
dgp => dg(xmap%x1(l)%i,xmap%x1(l)%j)
dgp = dgp + xmap%x1(l)%area*x(l)
end do
le = 0;
call mpp_sync_self() !Balaji
do p=0,xmap%npes-1
if (xmap%your1my2(p)) then
l = le + 1;
is = grid1%is(p); ie = grid1%ie(p);
js = grid1%js(p); je = grid1%je(p);
do j=js,je; do i=is,ie;
le = le + 1
xmap%send_buffer(le) = dg(i,j)
end do; end do;
! Force use of "scalar", integer pointer mpp interface.
call mpp_send(xmap%send_buffer(l), plen=le-l+1, to_pe=p+xmap%root_pe);
end if
end do
d = dg(grid1%is_me:grid1%ie_me,grid1%js_me:grid1%je_me);
im = grid1%ie_me-grid1%is_me+1;
jm = grid1%je_me-grid1%js_me+1;
do p=0,xmap%npes-1
if (xmap%your2my1(p)) d = d + get_side_1(p+xmap%root_pe,im,jm)
end do
!
! normalize with side 1 grid cell areas
!
d = d * grid1%area_inv
! call mpp_sync_self
end subroutine get_1_from_xgrid
!#######################################################################
subroutine get_1_from_xgrid_repro(d, x, xmap)
type (xmap_type), intent(inout) :: xmap
real, dimension(xmap%grids(1)%is_me:xmap%grids(1)%ie_me, &
xmap%grids(1)%js_me:xmap%grids(1)%je_me), intent(out) :: d
real, dimension(: ), intent(in ) :: x
real, dimension(:), allocatable :: x_psum
integer, dimension(:), allocatable :: pe_psum
integer :: l1, l2, l3, g, i, j, k, p
integer, dimension(0:xmap%npes-1) :: pl
type (grid_type), pointer, save :: grid =>NULL()
allocate ( x_psum (sum(xmap%send_count_repro)) )
allocate ( pe_psum (sum(xmap%send_count_repro)) )
x_psum = 0.0
l1 = 0 ! index into partition summed exchange grid variable
l2 = 0 ! index into exchange grid variable
do g=2,size(xmap%grids(:))
do l3=1,xmap%grids(g)%size ! index into this side 2 grid's patterns
l1 = l1 + 1
do k=1,xmap%grids(g)%km
i = xmap%grids(g)%x(l3)%i2
j = xmap%grids(g)%x(l3)%j2
if (xmap%grids(g)%frac_area(i,j,k)/=0.0) then
l2 = l2 + 1
x_psum (l1) = x_psum(l1) + xmap%x1(l2)%area * x(l2)
pe_psum(l1) = xmap%grids(g)%x(l3)%pe
end if
end do
end do
end do
l2 = 0;
call mpp_sync_self() !Balaji
do p=0,xmap%npes-1
l1 = l2 + 1
l2 = l2 + xmap%send_count_repro(p)
if (xmap%send_count_repro(p)>0) then ! can send to myself
xmap%send_buffer(l1:l2) = pack(x_psum, pe_psum==p+xmap%root_pe)
! Force use of "scalar", integer pointer mpp interface.
call mpp_send(xmap%send_buffer(l1), plen=l2-l1+1, to_pe=p+xmap%root_pe);
end if
end do
deallocate ( x_psum, pe_psum)
allocate ( x_psum (sum(xmap%recv_count_repro)) )
l2 = 0;
do p=0,xmap%npes-1
l1 = l2 + 1
l2 = l2 + xmap%recv_count_repro(p)
if (xmap%recv_count_repro(p)>0) then ! can receive from myself
! Force use of "scalar", integer pointer mpp interface.
call mpp_recv(x_psum(l1), glen=l2-l1+1, from_pe=p+xmap%root_pe);
pl(p) = l1
end if
end do
d = 0.0
do g=2,size(xmap%grids(:))
grid => xmap%grids(g)
do l3=1,grid%size_repro ! index into side1 grid's patterns
i = grid%x_repro(l3)%i1
j = grid%x_repro(l3)%j1
d(i,j) = d(i,j) + x_psum(pl(grid%x_repro(l3)%pe-xmap%root_pe))
pl(grid%x_repro(l3)%pe-xmap%root_pe) = pl(grid%x_repro(l3)%pe-xmap%root_pe) + 1
end do
end do
deallocate ( x_psum )
!
! normalize with side 1 grid cell areas
!
d = d * xmap%grids(1)%area_inv
! call mpp_sync_self
end subroutine get_1_from_xgrid_repro
!#######################################################################
!
!
!
!
!
!
! conservation_check - returns three numbers which are the global sum of a
! variable (1) on its home model grid, (2) after interpolation to the other
! side grid(s), and (3) after re_interpolation back onto its home side grid(s).
!
function conservation_check_side1(d, grid_id, xmap,remap_method) ! this one for 1->2->1
real, dimension(:,:), intent(in ) :: d
character(len=3), intent(in ) :: grid_id
type (xmap_type), intent(inout) :: xmap
real, dimension(3) :: conservation_check_side1
integer, intent(in), optional :: remap_method
real :: gsum
real, dimension(xmap%size) :: x_over, x_back
real, dimension(size(d,1),size(d,2)) :: d1
real, dimension(:,:,:), allocatable :: d2
integer :: g
type (grid_type), pointer, save :: grid1 =>NULL(), grid2 =>NULL()
grid1 => xmap%grids(1)
conservation_check_side1(1) = mpp_global_sum(grid1%domain, grid1%area*d)
conservation_check_side1(2) = 0.0
call put_to_xgrid (d, grid1%id, x_over, xmap, remap_method) ! put from side 1
do g=2,size(xmap%grids(:))
grid2 => xmap%grids(g)
allocate (d2 (grid2%is_me:grid2%ie_me, grid2%js_me:grid2%je_me, grid2%km) )
call get_from_xgrid (d2, grid2%id, x_over, xmap) ! get onto side 2's
conservation_check_side1(2) = conservation_check_side1(2) + &
mpp_global_sum( grid2%domain, grid2%area * sum(grid2%frac_area*d2,DIM=3) )
call put_to_xgrid (d2, grid2%id, x_back, xmap) ! put from side 2's
deallocate (d2)
end do
call get_from_xgrid(d1, grid1%id, x_back, xmap) ! get onto side 1
conservation_check_side1(3) = mpp_global_sum(grid1%domain, grid1%area*d1)
end function conservation_check_side1
!
!#######################################################################
!
! conservation_check - returns three numbers which are the global sum of a
! variable (1) on its home model grid, (2) after interpolation to the other
! side grid(s), and (3) after re_interpolation back onto its home side grid(s).
!
!
!
!
!
!
function conservation_check_side2(d, grid_id, xmap,remap_method) ! this one for 2->1->2
real, dimension(:,:,:), intent(in ) :: d
character(len=3), intent(in ) :: grid_id
type (xmap_type), intent(inout) :: xmap
real, dimension(3) :: conservation_check_side2
integer, intent(in), optional :: remap_method
real :: gsum
real, dimension(xmap%size) :: x_over, x_back
real, dimension(:,: ), allocatable :: d1
real, dimension(:,:,:), allocatable :: d2
integer :: g
type (grid_type), pointer, save :: grid1 =>NULL(), grid2 =>NULL()
grid1 => xmap%grids(1)
do g = 2,size(xmap%grids(:))
grid2 => xmap%grids(g)
if (grid_id==grid2%id) then
conservation_check_side2(1) = mpp_global_sum( grid2%domain, &
grid2%area * sum(grid2%frac_area*d,DIM=3) )
call put_to_xgrid(d, grid_id, x_over, xmap) ! put from this side 2
else
call put_to_xgrid(0 * grid2%frac_area, grid2%id, x_over, xmap) ! zero rest
end if
end do
allocate ( d1(size(grid1%area,1),size(grid1%area,2)) )
call get_from_xgrid(d1, grid1%id, x_over, xmap) ! get onto side 1
conservation_check_side2(2) = mpp_global_sum(grid1%domain, grid1%area*d1)
call put_to_xgrid(d1, grid1%id, x_back, xmap,remap_method) ! put from side 1
deallocate ( d1 )
conservation_check_side2(3) = 0.0;
do g = 2,size(xmap%grids(:))
grid2 => xmap%grids(g)
allocate ( d2 ( size(grid2%frac_area, 1), size(grid2%frac_area, 2), &
size(grid2%frac_area, 3) ) )
call get_from_xgrid(d2, grid2%id, x_back, xmap) ! get onto side 2's
conservation_check_side2(3) = conservation_check_side2(3) &
+mpp_global_sum( grid2%domain, &
grid2%area * sum(grid2%frac_area*d2,DIM=3) )
deallocate ( d2 )
end do
end function conservation_check_side2
!
!#######################################################################
! This function is used to calculate the gradient along zonal direction.
! Maybe need to setup a limit for the gradient.
function grad_zonal(d, lon, lat, is, ie, js, je, isd, ied, jsd, jed)
integer, intent(in) :: isd, ied, jsd, jed
real, dimension(isd:ied,jsd:jed), intent(in) :: d
real, dimension(:), intent(in) :: lon
real, dimension(:), intent(in) :: lat
integer, intent(in) :: is, ie, js, je
real, dimension(is:ie,js:je) :: grad_zonal
real :: dx, costheta
integer :: i, j, ip1, im1
! calculate the gradient of the data on each grid
do i = is, ie
if(i == 1) then
ip1 = i+1; im1 = i
else if(i==size(lon(:)) ) then
ip1 = i; im1 = i-1
else
ip1 = i+1; im1 = i-1
endif
dx = lon(ip1) - lon(im1)
if(abs(dx).lt.EPS ) call error_mesg('xgrids_mod', 'Improper grid size in lontitude', FATAL)
if(dx .gt. PI) dx = dx - 2.0* PI
if(dx .lt. -PI) dx = dx + 2.0* PI
do j = js, je
costheta = cos(lat(j))
if(abs(costheta) .lt. EPS) call error_mesg('xgrids_mod', 'Improper latitude grid', FATAL)
grad_zonal(i,j) = (d(ip1,j)-d(im1,j))/(dx*costheta)
enddo
enddo
return
end function grad_zonal
!#######################################################################
! This function is used to calculate the gradient along meridinal direction.
! Maybe need to setup a limit for the gradient.
function grad_merid(d, lat, is, ie, js, je, isd, ied, jsd, jed)
integer, intent(in) :: isd, ied, jsd, jed
real, dimension(isd:ied,jsd:jed), intent(in) :: d
real, dimension(:), intent(in) :: lat
integer, intent(in) :: is, ie, js, je
real, dimension(is:ie,js:je) :: grad_merid
real :: dy
integer :: i, j, jp1, jm1
! calculate the gradient of the data on each grid
do j = js, je
if(j == 1) then
jp1 = j+1; jm1 = j
else if(j == size(lat(:)) ) then
jp1 = j; jm1 = j-1
else
jp1 = j+1; jm1 = j-1
endif
dy = lat(jp1) - lat(jm1)
if(abs(dy).lt.EPS) call error_mesg('xgrids_mod', 'Improper grid size in latitude', FATAL)
do i = is, ie
grad_merid(i,j) = (d(i,jp1) - d(i,jm1))/dy
enddo
enddo
return
end function grad_merid
!#######################################################################
end module xgrid_mod
!
!
! A guide to grid coupling in FMS.
!
!
! A simple xgrid example.
!
!