FV3 Bundle
xgrid.F90
Go to the documentation of this file.
1 !***********************************************************************
2 !* GNU Lesser General Public License
3 !*
4 !* This file is part of the GFDL Flexible Modeling System (FMS).
5 !*
6 !* FMS is free software: you can redistribute it and/or modify it under
7 !* the terms of the GNU Lesser General Public License as published by
8 !* the Free Software Foundation, either version 3 of the License, or (at
9 !* your option) any later version.
10 !*
11 !* FMS is distributed in the hope that it will be useful, but WITHOUT
12 !* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 !* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14 !* for more details.
15 !*
16 !* You should have received a copy of the GNU Lesser General Public
17 !* License along with FMS. If not, see <http://www.gnu.org/licenses/>.
18 !***********************************************************************
19 
20 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
21 !
22 ! xgrid_mod - implements exchange grids. An exchange grid is the grid whose
23 ! boundary set is the union of the boundaries of the participating
24 ! grids. The exchange grid is the coarsest grid that is a
25 ! refinement of each of the participating grids. Every exchange
26 ! grid cell is a subarea of one and only one cell in each of the
27 ! participating grids. The exchange grid has two purposes:
28 !
29 ! (1) The exchange cell areas are used as weights for
30 ! conservative interpolation between model grids.
31 !
32 ! (2) Computation of surface fluxes takes place on it,
33 ! thereby using the finest scale data obtainable.
34 !
35 ! The exchange cells are the 2D intersections between cells of the
36 ! participating grids. They are computed elsewhere and are
37 ! read here from a NetCDF grid file as a sequence of quintuples
38 ! (i and j on each of two grids and the cell area).
39 !
40 ! Each processing element (PE) computes a subdomain of each of the
41 ! participating grids as well as a subset of the exchange cells.
42 ! The geographic regions corresponding to these subdomains will,
43 ! in general, not be the same so communication must occur between
44 ! the PEs. The scheme for doing this is as follows. A distinction
45 ! is drawn between the participating grids. There is a single
46 ! "side 1" grid and it does not have partitions (sub-grid surface
47 ! types). There are one or more "side 2" grids and they may have
48 ! more than 1 partition. In standard usage, the atmosphere grid is
49 ! on side 1 and the land and sea ice grids are on side 2. The set
50 ! of exchange cells computed on a PE corresponds to its side 2
51 ! geographic region(s). Communication between the PEs takes place
52 ! on the side 1 grid. Note: this scheme does not generally allow
53 ! reproduction of answers across varying PE counts. This is
54 ! because, in the side 1 "get", exchange cells are first summed
55 ! locally onto a side 1 grid, then these side 1 contributions are
56 ! further summed after they have been communicated to their target
57 ! PE. For the make_exchange_reproduce option, a special side 1 get
58 ! is used. This get communicates individual exchange cells. The
59 ! cells are summed in the order they appear in the grid spec. file.
60 ! Michael Winton (Michael.Winton@noaa.gov) Oct 2001
61 !
62 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
63 module xgrid_mod
64 
65 ! <CONTACT EMAIL="Michael.Winton@noaa.gov">
66 ! Michael Winton
67 ! </CONTACT>
68 ! <CONTACT EMAIL="Zhi.Liang@noaa.gov">
69 ! Zhi Liang
70 ! </CONTACT>
71 
72 ! <HISTORY SRC="http://www.gfdl.noaa.gov/fms-cgi-bin/cvsweb.cgi/FMS/"/>
73 
74 ! <OVERVIEW>
75 ! <TT>xgrid_mod</TT> implements exchange grids for coupled models running on
76 ! multiple processors. An exchange grid is formed from the union of
77 ! the bounding lines of the two (logically rectangular) participating
78 ! grids. The exchange grid is therefore the coarsest grid that is a
79 ! refinement of both participating grids. Exchange grids are used for
80 ! two purposes by coupled models: (1) conservative interpolation of fields
81 ! between models uses the exchange grid cell areas as weights and
82 ! (2) the surface flux calculation takes place on the exchange grid thereby
83 ! using the finest scale data available. <TT>xgrid_mod</TT> uses a NetCDF grid
84 ! specification file containing the grid cell overlaps in combination with
85 ! the <LINK SRC="ftp://ftp.gfdl.gov/pub/vb/mpp/mpp_domains.F90">
86 ! <TT>mpp_domains</TT></LINK> domain decomposition information to determine
87 ! the grid and processor connectivities.
88 ! </OVERVIEW>
89 
90 ! <DESCRIPTION>
91 ! <TT>xgrid_mod</TT> is initialized with a list of model identifiers (three characters
92 ! each), a list of <TT>mpp_domains</TT> domain data structures, and a grid specification
93 ! file name. The first element in the lists refers to the "side one" grid.
94 ! The remaining elements are on "side two". Thus, there may only be a single
95 ! side one grid and it is further restricted to have no partitions (sub-grid
96 ! areal divisions). In standard usage, the atmosphere model is on side one
97 ! and the land and sea ice models are on side two. <TT>xgrid_mod</TT> performs
98 ! interprocessor communication on the side one grid. Exchange grid variables
99 ! contain no data for zero sized partitions. The size and format of exchange
100 ! grid variables change every time the partition sizes or number of partitions
101 ! are modified with a <TT>set_frac_area</TT> call on a participating side two grid.
102 ! Existing exchange grid variables cannot be properly interpreted after
103 ! that time; new ones must be allocated and assigned with the <TT>put_to_xgrid</TT>
104 ! call.
105 ! </DESCRIPTION>
106 
107 ! <DATA NAME="xmap_type" TYPE="" >
108 ! The fields of xmap_type are all private.
109 ! </DATA>
110 
111 ! <DATASET NAME="">
112 ! <TT>xgrid_mod</TT> reads a NetCDF grid specification file to determine the
113 ! grid and processor connectivities. The exchange grids are defined
114 ! by a sequence of quintuples: the <TT>i/j</TT> indices of the intersecting
115 ! cells of the two participating grids and their areal overlap.
116 ! The names of the five fields are generated automatically from the
117 ! three character ids of the participating grids. For example, if
118 ! the side one grid id is "ATM" and the side two grid id is "OCN",
119 ! <TT>xgrid_mod</TT> expects to find the following five fields in the grid
120 ! specification file: <TT>I_ATM_ATMxOCN, J_ATM_ATMxOCN, I_OCN_ATMxOCN,
121 ! J_OCN_ATMxOCN, and AREA_ATMxOCN</TT>. These fields may be generated
122 ! by the <TT>make_xgrids</TT> utility.
123 ! </DATASET>
124 
125 #include <fms_platform.h>
126 
127 use fms_mod, only: file_exist, open_namelist_file, check_nml_error, &
128  error_mesg, close_file, fatal, note, stdlog, &
129  write_version_number, read_data, field_exist, &
130  field_size, lowercase, string, &
131  get_mosaic_tile_grid
132 use fms_io_mod, only: get_var_att_value
133 use mpp_mod, only: mpp_npes, mpp_pe, mpp_root_pe, mpp_send, mpp_recv, &
134  mpp_sync_self, stdout, mpp_max, event_recv, &
135  mpp_get_current_pelist, mpp_clock_id, mpp_min, &
136  mpp_alltoall, &
137  mpp_clock_begin, mpp_clock_end, mpp_clock_sync, &
138  comm_tag_1, comm_tag_2, comm_tag_3, comm_tag_4, &
139  comm_tag_5, comm_tag_6, comm_tag_7, comm_tag_8, &
140  comm_tag_9, comm_tag_10
141 use mpp_mod, only: input_nml_file, mpp_set_current_pelist, mpp_sum, mpp_sync
145  yupdate, mpp_get_current_ntile, mpp_get_tile_id, &
146  mpp_get_ntile_count, mpp_get_tile_list, &
149  mpp_get_domain_npes, mpp_get_domain_root_pe, &
150  mpp_domain_is_initialized, mpp_broadcast_domain, &
151  mpp_get_domain_pelist, mpp_compute_extent, &
152  domainug, mpp_get_ug_compute_domains, &
153  mpp_get_ug_domains_index, mpp_get_ug_domain_grid_index, &
154  mpp_get_ug_domain_tile_list, mpp_pass_sg_to_ug
155 use mpp_io_mod, only: mpp_open, mpp_multi, mpp_single, mpp_overwr
156 use constants_mod, only: pi, radius
160 
163 use gradient_mod, only: gradient_cubic
164 
165 implicit none
166 private
167 
170 ! AREA_ATM_SPHERE, AREA_LND_SPHERE, AREA_OCN_SPHERE, &
171 ! AREA_ATM_MODEL, AREA_LND_MODEL, AREA_OCN_MODEL, &
177 
178 !--- paramters that determine the remapping method
179 integer, parameter :: first_order = 1
180 integer, parameter :: second_order = 2
181 integer, parameter :: version1 = 1 ! grid spec file
182 integer, parameter :: version2 = 2 ! mosaic grid file
183 integer, parameter :: max_fields = 80
184 
185 ! <NAMELIST NAME="xgrid_nml">
186 ! <DATA NAME="make_exchange_reproduce" TYPE="logical" DEFAULT=".false.">
187 ! Set to .true. to make <TT>xgrid_mod</TT> reproduce answers on different
188 ! numbers of PEs. This option has a considerable performance impact.
189 ! </DATA>
190 ! <DATA NAME="interp_method" TYPE="character(len=64)" DEFAULT=" 'first_order' ">
191 ! exchange grid interpolation method. It has two options:
192 ! "first_order", "second_order".
193 ! </DATA>
194 ! <DATA NAME="xgrid_log" TYPE="logical" DEFAULT=" .false. ">
195 ! Outputs exchange grid information to xgrid.out.<pe> for debug/diag purposes.
196 ! </DATA>
197 ! <DATA NAME="nsubset" TYPE="integer" DEFAULT="0">
198 ! number of processors to read exchange grid information. Those processors that read
199 ! the exchange grid information will send data to other processors to prepare for flux exchange.
200 ! Default value is 0. When nsubset is 0, each processor will read part of the exchange grid
201 ! information. The purpose of this namelist is to improve performance of setup_xmap when running
202 ! on highr processor count and solve receiving size mismatch issue on high processor count.
203 ! Try to set nsubset = mpp_npes/MPI_rank_per_node.
204 ! </DATA>
205 logical :: make_exchange_reproduce = .false. ! exactly same on different # PEs
206 logical :: xgrid_log = .false.
207 character(len=64) :: interp_method = 'first_order'
208 logical :: debug_stocks = .false.
209 logical :: xgrid_clocks_on = .false.
210 logical :: monotonic_exchange = .false.
211 integer :: nsubset = 0 ! 0 means mpp_npes()
212 logical :: do_alltoall = .true.
213 logical :: do_alltoallv = .false.
216 ! </NAMELIST>
217 logical :: init = .true.
219 
220 ! Area elements used inside each model
221 real, allocatable, dimension(:,:) :: area_atm_model, area_lnd_model, area_ocn_model
222 ! Area elements based on a the spherical model used by the ICE layer
223 real, allocatable, dimension(:,:) :: area_atm_sphere, area_lnd_sphere, area_ocn_sphere
224 
225 ! <INTERFACE NAME="put_to_xgrid">
226 
227 ! <OVERVIEW>
228 ! Scatters data from model grid onto exchange grid.
229 ! </OVERVIEW>
230 ! <DESCRIPTION>
231 ! Scatters data from model grid onto exchange grid.
232 ! </DESCRIPTION>
233 ! <TEMPLATE>
234 ! call put_to_xgrid(d, grid_id, x, xmap, remap_order)
235 ! </TEMPLATE>
236 ! <IN NAME="d" TYPE="real" > </IN>
237 ! <IN NAME="grid_id" TYPE=" character(len=3)" > </IN>
238 ! <INOUT NAME="x" TYPE="real" > </INOUT>
239 ! <INOUT NAME="xmap" TYPE="xmap_type" > </INOUT>
240 ! <IN NAME="remap_method" TYPE="integer,optional">
241 ! exchange grid interpolation method. It has four possible values:
242 ! FIRST_ORDER (=1), SECOND_ORDER(=2). Default value is FIRST_ORDER.
243 ! </IN>
244 interface put_to_xgrid
245  module procedure put_side1_to_xgrid
246  module procedure put_side2_to_xgrid
247 end interface
248 ! </INTERFACE>
249 
250 ! <INTERFACE NAME="get_from_xgrid">
251 
252 ! <OVERVIEW>
253 ! Sums data from exchange grid to model grid.
254 ! </OVERVIEW>
255 ! <DESCRIPTION>
256 ! Sums data from exchange grid to model grid.
257 ! </DESCRIPTION>
258 ! <TEMPLATE>
259 ! call get_from_xgrid(d, grid_id, x, xmap)
260 ! </TEMPLATE>
261 ! <IN NAME="x" TYPE="real" > </IN>
262 ! <IN NAME="grid_id" TYPE=" character(len=3)" > </IN>
263 ! <OUT NAME="d" TYPE="real" > </OUT>
264 ! <INOUT NAME="xmap" TYPE="xmap_type" > </INOUT>
265 interface get_from_xgrid
266  module procedure get_side1_from_xgrid
267  module procedure get_side2_from_xgrid
268 end interface
269 ! </INTERFACE>
270 
272  module procedure put_side1_to_xgrid_ug
273  module procedure put_side2_to_xgrid_ug
274 end interface
275 
277  module procedure get_side2_from_xgrid_ug
278  module procedure get_side1_from_xgrid_ug
279 end interface
280 
281 interface set_frac_area
282  module procedure set_frac_area_sg
283  module procedure set_frac_area_ug
284 end interface
285 
286 ! <INTERFACE NAME="conservation_check">
287 
288 ! <OVERVIEW>
289 ! Returns three numbers which are the global sum of a variable.
290 ! </OVERVIEW>
291 ! <DESCRIPTION>
292 ! Returns three numbers which are the global sum of a
293 ! variable (1) on its home model grid, (2) after interpolation to the other
294 ! side grid(s), and (3) after re_interpolation back onto its home side grid(s).
295 ! Conservation_check must be called by all PEs to work properly.
296 ! </DESCRIPTION>
297 ! <TEMPLATE>
298 ! call conservation_check(d, grid_id, xmap,remap_order)
299 ! </TEMPLATE>
300 ! <IN NAME="d" TYPE="real" DIM="(:,:)" > </IN>
301 ! <IN NAME="grid_id" TYPE="character(len=3)" > </IN>
302 ! <INOUT NAME="xmap" TYPE="xmap_type" > </INOUT>
303 ! <OUT NAME="" TYPE="real" DIM="3">The global sum of a variable.</OUT>
304 ! <IN NAME="remap_method" TYPE="integer,optional">
305 ! </IN>
307  module procedure conservation_check_side1
308  module procedure conservation_check_side2
309 end interface
310 ! </INTERFACE>
312  module procedure conservation_check_ug_side1
313  module procedure conservation_check_ug_side2
314 end interface
315 
316 
318  integer :: i1, j1, i2, j2 ! indices of cell in model arrays on both sides
319  integer :: l1, l2
320  integer :: recv_pos ! position in the receive buffer.
321  integer :: pe ! other side pe that has this cell
322  integer :: tile ! tile index of side 1 mosaic.
323  real :: area ! geographic area of exchange cell
324 ! real :: area1_ratio !(= x_area/grid1_area), will be added in the future to improve efficiency
325 ! real :: area2_ratio !(= x_area/grid2_area), will be added in the future to improve efficiency
326  real :: di, dj ! Weight for the gradient of flux
327  real :: scale
328 end type xcell_type
329 
331  real, dimension(:,:), pointer :: dx => null()
332  real, dimension(:,:), pointer :: dy => null()
333  real, dimension(:,:), pointer :: area => null()
334  real, dimension(:), pointer :: edge_w => null()
335  real, dimension(:), pointer :: edge_e => null()
336  real, dimension(:), pointer :: edge_s => null()
337  real, dimension(:), pointer :: edge_n => null()
338  real, dimension(:,:,:), pointer :: en1 => null()
339  real, dimension(:,:,:), pointer :: en2 => null()
340  real, dimension(:,:,:), pointer :: vlon => null()
341  real, dimension(:,:,:), pointer :: vlat => null()
342 end type grid_box_type
343 
345  character(len=3) :: id ! grid identifier
346  integer :: npes ! number of processor on this grid.
347  logical :: on_this_pe ! indicate the domain is defined on this pe
348  integer :: root_pe ! indicate the root pe of the domain
349  integer, pointer, dimension(:) :: pelist ! pelist of the domain
350  integer :: ntile ! number of tiles in mosaic
351  integer :: ni, nj ! max of global size of all the tiles
352  integer, pointer, dimension(:) :: tile =>null() ! tile id ( pe index )
353  integer, pointer, dimension(:) :: is =>null(), ie =>null() ! domain - i-range (pe index)
354  integer, pointer, dimension(:) :: js =>null(), je =>null() ! domain - j-range (pe index)
355  integer, pointer :: is_me =>null(), ie_me =>null() ! my domain - i-range
356  integer, pointer :: js_me =>null(), je_me =>null() ! my domain - j-range
357  integer :: isd_me, ied_me ! my data domain - i-range
358  integer :: jsd_me, jed_me ! my data domain - j-range
359  integer :: nxd_me, nyd_me ! data domain size
360  integer :: nxc_me, nyc_me ! compute domain size
361  integer, pointer :: tile_me ! my tile id
362  integer :: im , jm , km ! global domain range
363  real, pointer, dimension(:) :: lon =>null(), lat =>null() ! center of global grids
364  real, pointer, dimension(:,:) :: geolon=>null(), geolat=>null() ! geographical grid center
365  real, pointer, dimension(:,:,:) :: frac_area =>null() ! partition fractions
366  real, pointer, dimension(:,:) :: area =>null() ! cell area
367  real, pointer, dimension(:,:) :: area_inv =>null() ! 1 / area for normalization
368  integer :: first, last ! xgrid index range
369  integer :: first_get, last_get ! xgrid index range for get_2_from_xgrid
370  integer :: size ! # xcell patterns
371  type(xcell_type), pointer :: x(:) =>null() ! xcell patterns
372  integer :: size_repro ! # side 1 patterns for repro
373  type(xcell_type), pointer :: x_repro(:) =>null() ! side 1 patterns for repro
374  type(domain2d) :: domain ! used for conservation checks
375  type(domain2d) :: domain_with_halo ! used for second order remapping
376  logical :: is_latlon ! indicate if the grid is lat-lon grid or not.
377  type(grid_box_type) :: box ! used for second order remapping.
378  !--- The following is for land unstruct domain
379  logical :: is_ug
380  integer :: nxl_me
381  integer, pointer :: ls_me =>null(), le_me =>null() ! unstruct domain
382  integer, pointer, dimension(:) :: ls =>null(), le =>null()
383  integer, pointer :: gs_me =>null(), ge_me =>null()
384  integer, pointer, dimension(:) :: gs =>null(), ge =>null()
385  integer, pointer, dimension(:) :: l_index =>null()
386  type(domainug) :: ug_domain
387 
388 end type grid_type
389 
391  integer :: i, j
392  real :: area ! (= geographic area * frac_area)
393 ! real :: area_ratio !(= x1_area/grid1_area) ! will be added in the future to improve efficiency
394  real :: di, dj ! weight for the gradient of flux
395  integer :: tile ! tile index of side 1 mosaic.
396  integer :: pos
397 end type x1_type
398 
400  integer :: i, j, l, k, pos
401  real :: area ! geographic area of exchange cell
402 ! real :: area_ratio !(=x2_area/grid2_area ) ! will be added in the future to improve efficiency
403 end type x2_type
404 
406  integer :: count
407  integer :: pe
408  integer :: buffer_pos
409  integer, _allocatable :: i(:) _null
410  integer, _allocatable :: j(:) _null
411  integer, _allocatable :: g(:) _null
412  integer, _allocatable :: xloc(:) _null
413  integer, _allocatable :: tile(:) _null
414  real, _allocatable :: di(:) _null
415  real, _allocatable :: dj(:) _null
416 end type overlap_type
417 
419  integer :: nsend, nrecv
420  integer :: sendsize, recvsize
421  integer, pointer, dimension(:) :: unpack_ind=>null()
422  type(overlap_type), pointer, dimension(:) :: send=>null()
423  type(overlap_type), pointer, dimension(:) :: recv=>null()
424 end type comm_type
425 
427  private
428  integer :: size ! # of exchange grid cells with area > 0 on this pe
429  integer :: size_put1 ! # of exchange grid cells for put_1_to_xgrid
430  integer :: size_get2 ! # of exchange grid cells for get_2_to_xgrid
431  integer :: me, npes, root_pe
432  logical, pointer, dimension(:) :: your1my2 =>null()! true if side 1 domain on
433  ! indexed pe overlaps side 2
434  ! domain on this pe
435  logical, pointer, dimension(:) :: your2my1 =>null() ! true if a side 2 domain on
436  ! indexed pe overlaps side 1
437  ! domain on this pe
438  integer, pointer, dimension(:) :: your2my1_size=>null() ! number of exchange grid of
439  ! a side 2 domain on
440  ! indexed pe overlaps side 1
441  ! domain on this pe
442 
443  type(grid_type), pointer, dimension(:) :: grids =>null() ! 1st grid is side 1;
444  ! rest on side 2
445  !
446  ! Description of the individual exchange grid cells (index is cell #)
447  !
448  type(x1_type), pointer, dimension(:) :: x1 =>null() ! side 1 info
449  type(x1_type), pointer, dimension(:) :: x1_put =>null() ! side 1 info
450  type(x2_type), pointer, dimension(:) :: x2 =>null() ! side 2 info
451  type(x2_type), pointer, dimension(:) :: x2_get =>null() ! side 2 info
452 
453  integer, pointer, dimension(:) :: send_count_repro =>null()
454  integer, pointer, dimension(:) :: recv_count_repro =>null()
455  integer :: send_count_repro_tot ! sum(send_count_repro)
456  integer :: recv_count_repro_tot ! sum(recv_count_repro)
457  integer :: version ! version of xgrids. version=VERSION! is for grid_spec file
458  ! and version=VERSION2 is for mosaic grid.
459  integer, pointer, dimension(:) :: ind_get1 =>null() ! indx for side1 get and side2 put.
460  integer, pointer, dimension(:) :: ind_put1 =>null() ! indx for side1 put and side 2get.
461  type(comm_type), pointer :: put1 =>null() ! for put_1_to_xgrid
462  type(comm_type), pointer :: get1 =>null() ! for get_1_from_xgrid
463  type(comm_type), pointer :: get1_repro =>null()! for get_1_from_xgrid_repro
464 end type xmap_type
465 
466 !-----------------------------------------------------------------------
467 ! Include variable "version" to be written to log file.
468 #include<file_version.h>
469 
470  real, parameter :: eps = 1.0e-10
471  real, parameter :: large_number = 1.e20
472  logical :: module_is_initialized = .false.
475  integer :: id_get_1_from_xgrid = 0
477  integer :: id_get_2_from_xgrid = 0
478  integer :: id_put_2_to_xgrid = 0
479  integer :: id_setup_xmap = 0
483 
484 
485  ! The following is for nested model
486  integer :: nnest=0, tile_nest, tile_parent
487  integer :: is_nest=0, ie_nest=0, js_nest=0, je_nest=0
488  integer :: is_parent=0, ie_parent=0, js_parent=0, je_parent=0
489 
490  ! The following is required to compute stocks of water, heat, ...
491 
492  interface stock_move
493  module procedure stock_move_3d, stock_move_2d
494  end interface
495 
496  interface stock_move_ug
497  module procedure stock_move_ug_3d
498  end interface
499 
502 
503 contains
504 
505 !#######################################################################
506 
507 logical function in_box(i, j, is, ie, js, je)
508  integer, intent(in) :: i, j, is, ie, js, je
509 
510  in_box = (i>=is) .and. (i<=ie) .and. (j>=js) .and. (j<=je)
511 end function in_box
512 
513 !#######################################################################
514 
515 ! <SUBROUTINE NAME="xgrid_init">
516 
517 ! <OVERVIEW>
518 ! Initialize the xgrid_mod.
519 ! </OVERVIEW>
520 ! <DESCRIPTION>
521 ! Initialization routine for the xgrid module. It reads the xgrid_nml,
522 ! writes the version information and xgrid_nml to the log file.
523 ! </DESCRIPTION>
524 ! <TEMPLATE>
525 ! call xgrid_init ( )
526 ! </TEMPLATE>
527 ! <OUT NAME="remap_method" TYPE="integer">
528 ! exchange grid interpolation method. It has four possible values:
529 ! FIRST_ORDER (=1), SECOND_ORDER(=2).
530 ! </OUT>
531 subroutine xgrid_init(remap_method)
532  integer, intent(out) :: remap_method
533 
534  integer :: unit, ierr, io, out_unit
535 
536  if (module_is_initialized) return
537  module_is_initialized = .true.
538 
539 
540 #ifdef INTERNAL_FILE_NML
541  read (input_nml_file, xgrid_nml, iostat=io)
542  ierr = check_nml_error( io, 'xgrid_nml' )
543 #else
544  if ( file_exist( 'input.nml' ) ) then
545  unit = open_namelist_file( )
546  ierr = 1
547  do while ( ierr /= 0 )
548  read ( unit, nml = xgrid_nml, iostat = io, end = 10 )
549  ierr = check_nml_error( io, 'xgrid_nml' )
550  enddo
551  10 continue
552  call close_file ( unit )
553  endif
554 #endif
555 
556 !--------- write version number and namelist ------------------
557  call write_version_number("XGRID_MOD", version)
558 
559  unit = stdlog( )
560  out_unit = stdout()
561  if ( mpp_pe() == mpp_root_pe() ) write (unit,nml=xgrid_nml)
562  call close_file (unit)
563 
564 !--------- check interp_method has suitable value
565 !--- when monotonic_exchange is true, interp_method must be second order.
566 
567  select case(trim(interp_method))
568  case('first_order')
569  remap_method = first_order
570  if( monotonic_exchange ) call error_mesg('xgrid_mod', &
571  'xgrid_nml monotonic_exchange must be .false. when interp_method = first_order', fatal)
572  write(out_unit,*)"NOTE from xgrid_mod: use first_order conservative exchange"
573  case('second_order')
574  if(monotonic_exchange) then
575  write(out_unit,*)"NOTE from xgrid_mod: use monotonic second_order conservative exchange"
576  else
577  write(out_unit,*)"NOTE from xgrid_mod: use second_order conservative exchange"
578  endif
579  remap_method = second_order
580  case default
581  call error_mesg('xgrid_mod', ' nml interp_method = ' //trim(interp_method)// &
582  ' is not a valid namelist option', fatal)
583  end select
584 
585  if(xgrid_clocks_on) then
586  id_put_1_to_xgrid_order_1 = mpp_clock_id("put_1_to_xgrid_order_1", flags=mpp_clock_sync)
587  id_put_1_to_xgrid_order_2 = mpp_clock_id("put_1_to_xgrid_order_2", flags=mpp_clock_sync)
588  id_get_1_from_xgrid = mpp_clock_id("get_1_from_xgrid", flags=mpp_clock_sync)
589  id_get_1_from_xgrid_repro = mpp_clock_id("get_1_from_xgrid_repro", flags=mpp_clock_sync)
590  id_get_2_from_xgrid = mpp_clock_id("get_2_from_xgrid", flags=mpp_clock_sync)
591  id_put_2_to_xgrid = mpp_clock_id("put_2_to_xgrid", flags=mpp_clock_sync)
592  id_setup_xmap = mpp_clock_id("setup_xmap", flags=mpp_clock_sync)
593  id_set_comm = mpp_clock_id("set_comm")
594  id_regen = mpp_clock_id("regen")
595  id_conservation_check = mpp_clock_id("conservation_check")
596  id_load_xgrid = mpp_clock_id("load_xgrid")
597  id_load_xgrid1 = mpp_clock_id("load_xgrid1")
598  id_load_xgrid2 = mpp_clock_id("load_xgrid2")
599  id_load_xgrid3 = mpp_clock_id("load_xgrid3")
600  id_load_xgrid4 = mpp_clock_id("load_xgrid4")
601  id_load_xgrid5 = mpp_clock_id("load_xgrid5")
602  endif
603 
604  remapping_method = remap_method
605 
606 end subroutine xgrid_init
607 ! </SUBROUTINE>
608 
609 !#######################################################################
610 
611 subroutine load_xgrid (xmap, grid, grid_file, grid1_id, grid_id, tile1, tile2, use_higher_order)
612 type(xmap_type), intent(inout) :: xmap
613 type(grid_type), intent(inout) :: grid
614 character(len=*), intent(in) :: grid_file
615 character(len=3), intent(in) :: grid1_id, grid_id
616 integer, intent(in) :: tile1, tile2
617 logical, intent(in) :: use_higher_order
618 
619  integer, pointer, dimension(:) :: i1=>null(), j1=>null()
620  integer, pointer, dimension(:) :: i2=>null(), j2=>null()
621  real, pointer, dimension(:) :: di=>null(), dj=>null()
622  real, pointer, dimension(:) :: area =>null()
623  integer, pointer, dimension(:) :: i1_tmp=>null(), j1_tmp=>null()
624  integer, pointer, dimension(:) :: i2_tmp=>null(), j2_tmp=>null()
625  real, pointer, dimension(:) :: di_tmp=>null(), dj_tmp=>null()
626  real, pointer, dimension(:) :: area_tmp =>null()
627  integer, pointer, dimension(:) :: i1_side1=>null(), j1_side1=>null()
628  integer, pointer, dimension(:) :: i2_side1=>null(), j2_side1=>null()
629  real, pointer, dimension(:) :: di_side1=>null(), dj_side1=>null()
630  real, pointer, dimension(:) :: area_side1 =>null()
631 
632  real, allocatable, dimension(:,:) :: tmp
633  real, allocatable, dimension(:) :: send_buffer, recv_buffer
634  type(grid_type), pointer, save :: grid1 =>null()
635  integer :: l, ll, ll_repro, p, siz(4), nxgrid, size_prev
636  type(xcell_type), allocatable :: x_local(:)
637  integer :: size_repro, out_unit
638  logical :: scale_exist = .false.
639  logical :: is_distribute = .false.
640  real, allocatable, dimension(:) :: scale
641  real :: garea
642  integer :: npes, isc, iec, nxgrid_local, pe, nxgrid_local_orig
643  integer :: nxgrid1, nxgrid2, nset1, nset2, ndivs, cur_ind
644  integer :: pos, nsend, nrecv, l1, l2, n, mypos, m
645  integer :: start(4), nread(4)
646  logical :: found
647  character(len=128) :: attvalue
648  integer, dimension(0:xmap%npes-1) :: pelist
649  logical, dimension(0:xmap%npes-1) :: subset_rootpe
650  integer, dimension(0:xmap%npes-1) :: nsend1, nsend2, nrecv1, nrecv2
651  integer, dimension(0:xmap%npes-1) :: send_cnt, recv_cnt
652  integer, dimension(0:xmap%npes-1) :: send_buffer_pos, recv_buffer_pos
653  integer, dimension(0:xmap%npes-1) :: ibegin, iend, pebegin, peend
654  integer, dimension(2*xmap%npes) :: ibuf1, ibuf2
655  integer, dimension(0:xmap%npes-1) :: pos_x, y2m1_size
656  integer, allocatable, dimension(:) :: y2m1_pe
657  integer, pointer, save :: iarray(:), jarray(:)
658  integer, allocatable, save :: pos_s(:)
659  integer, pointer, dimension(:) :: iarray2(:)=>null(), jarray2(:)=>null()
660  logical :: last_grid
661  integer :: nxgrid1_old
662  integer :: lll
663 
664  scale_exist = .false.
665  grid1 => xmap%grids(1)
666  out_unit = stdout()
667  npes = xmap%npes
668  pe = mpp_pe()
669  mypos = mpp_pe()-mpp_root_pe()
670 
671  call mpp_get_current_pelist(pelist)
672  !--- make sure npes = pelist(npes-1) - pelist(0) + 1
673  if( npes .NE. pelist(npes-1) - pelist(0) + 1 ) then
674  print*, "npes =", npes, ", pelist(npes-1)=", pelist(npes-1), ", pelist(0)=", pelist(0)
675  call error_mesg('xgrid_mod', .NE.'npes pelist(npes-1) - pelist(0)', fatal)
676  endif
677 
678  select case(xmap%version)
679  case(version1)
680  call field_size(grid_file, 'AREA_'//grid1_id//'x'//grid_id, siz)
681  nxgrid = siz(1);
682  if(nxgrid .LE. 0) return
683  case(version2)
684  !--- max_size is the exchange grid size between super grid.
685  nxgrid = get_mosaic_xgrid_size(grid_file)
686  if(nxgrid .LE. 0) return
687  end select
688 
689  !--- define a domain to read exchange grid.
690  if(nxgrid > npes) then
691  ndivs = npes
692  if(nsubset >0 .AND. nsubset < npes) ndivs = nsubset
693  call mpp_compute_extent( 1, nxgrid, ndivs, ibegin, iend)
694  if(npes == ndivs) then
695  p = mpp_pe()-mpp_root_pe()
696  isc = ibegin(p)
697  iec = iend(p)
698  subset_rootpe(:) = .true.
699  else
700  isc = 0; iec = -1
701  call mpp_compute_extent(pelist(0), pelist(npes-1), ndivs, pebegin, peend)
702  do n = 0, ndivs-1
703  if(pe == pebegin(n)) then
704  isc = ibegin(n)
705  iec = iend(n)
706  exit
707  endif
708  enddo
709  cur_ind = 0
710  subset_rootpe(:) = .false.
711 
712  do n = 0, npes-1
713  if(pelist(n) == pebegin(cur_ind)) then
714  subset_rootpe(n) = .true.
715  cur_ind = cur_ind+1
716  if(cur_ind == ndivs) exit
717  endif
718  enddo
719  endif
720  is_distribute = .true.
721  else
722  is_distribute = .false.
723  isc = 1; iec = nxgrid
724  endif
725 
726  nset1 = 5
727  nset2 = 5
728  if(use_higher_order) then
729  nset1 = nset1 + 2
730  nset2 = nset2 + 2
731  end if
732  if(scale_exist) nset2 = nset1 + 1
733 
734  call mpp_clock_begin(id_load_xgrid1)
735  if(iec .GE. isc) then
736  nxgrid_local = iec - isc + 1
737  allocate(i1_tmp(isc:iec), j1_tmp(isc:iec), i2_tmp(isc:iec), j2_tmp(isc:iec), area_tmp(isc:iec) )
738  if(use_higher_order) allocate(di_tmp(isc:iec), dj_tmp(isc:iec))
739 
740  start = 1; nread = 1
741 
742  select case(xmap%version)
743  case(version1)
744  start(1) = isc; nread(1) = nxgrid_local
745  allocate(tmp(nxgrid_local,1))
746  call read_data(grid_file, 'I_'//grid1_id//'_'//grid1_id//'x'//grid_id, tmp, start, nread, no_domain=.true.)
747  i1_tmp = tmp(:,1)
748  call read_data(grid_file, 'J_'//grid1_id//'_'//grid1_id//'x'//grid_id, tmp, start, nread, no_domain=.true.)
749  j1_tmp = tmp(:,1)
750  call read_data(grid_file, 'I_'//grid_id//'_'//grid1_id//'x'//grid_id, tmp, start, nread, no_domain=.true.)
751  i2_tmp = tmp(:,1)
752  call read_data(grid_file, 'J_'//grid_id//'_'//grid1_id//'x'//grid_id, tmp, start, nread, no_domain=.true.)
753  j2_tmp = tmp(:,1)
754  call read_data(grid_file, 'AREA_'//grid1_id//'x'//grid_id, tmp, start, nread, no_domain=.true.)
755  area_tmp = tmp(:,1)
756  if(use_higher_order) then
757  call read_data(grid_file, 'DI_'//grid1_id//'x'//grid_id, tmp, start, nread, no_domain=.true.)
758  di_tmp = tmp(:,1)
759  call read_data(grid_file, 'DJ_'//grid1_id//'x'//grid_id, tmp, start, nread, no_domain=.true.)
760  dj_tmp = tmp(:,1)
761  end if
762  deallocate(tmp)
763  case(version2)
764  nread(1) = 2; start(2) = isc; nread(2) = nxgrid_local
765  allocate(tmp(2, isc:iec))
766  call read_data(grid_file, "tile1_cell", tmp, start, nread, no_domain=.true.)
767  i1_tmp(isc:iec) = tmp(1, isc:iec)
768  j1_tmp(isc:iec) = tmp(2, isc:iec)
769  call read_data(grid_file, "tile2_cell", tmp, start, nread, no_domain=.true.)
770  i2_tmp(isc:iec) = tmp(1, isc:iec)
771  j2_tmp(isc:iec) = tmp(2, isc:iec)
772  if(use_higher_order) then
773  call read_data(grid_file, "tile1_distance", tmp, start, nread, no_domain=.true.)
774  di_tmp(isc:iec) = tmp(1, isc:iec)
775  dj_tmp(isc:iec) = tmp(2, isc:iec)
776  end if
777  start = 1; nread = 1
778  start(1) = isc; nread(1) = nxgrid_local
779  deallocate(tmp)
780  allocate(tmp(isc:iec,1) )
781  call read_data(grid_file, "xgrid_area", tmp(:,1:1), start, nread, no_domain=.true.)
782  ! check the units of "xgrid_area
783  call get_var_att_value(grid_file, "xgrid_area", "units", attvalue)
784  if( trim(attvalue) == 'm2' ) then
785  garea = 4.0*pi*radius*radius;
786  area_tmp = tmp(:,1)/garea
787  else if( trim(attvalue) == 'none' ) then
788  area_tmp = tmp(:,1)
789  else
790  call error_mesg('xgrid_mod', 'In file '//trim(grid_file)//', xgrid_area units = '// &
791  trim(attvalue)//' should be "m2" or "none"', fatal)
792  endif
793 
794  !--- if field "scale" exist, read this field. Normally this
795  !--- field only exist in landXocean exchange grid cell.
796  if(grid1_id == 'LND' .AND. grid_id == 'OCN') then
797  if(field_exist(grid_file, "scale")) then
798  allocate(scale(isc:iec))
799  write(out_unit, *)"NOTE from load_xgrid(xgrid_mod): field 'scale' exist in the file "// &
800  trim(grid_file)//", this field will be read and the exchange grid cell area will be multiplied by scale"
801  call read_data(grid_file, "scale", tmp, start, nread, no_domain=.true.)
802  scale = tmp(:,1)
803  scale_exist = .true.
804  endif
805  endif
806  deallocate(tmp)
807  end select
808 
809  !---z1l: The following change is for the situation that some processor is masked out.
810  !---loop through all the pe to see if side 1 and side of each exchange grid is on some processor
811  nxgrid_local_orig = nxgrid_local
812  allocate(i1(isc:iec), j1(isc:iec), i2(isc:iec), j2(isc:iec), area(isc:iec) )
813  if(use_higher_order) allocate(di(isc:iec), dj(isc:iec))
814  pos = isc-1
815  do l = isc, iec
816  found = .false.
817  !--- first check if the exchange grid is on one of side 1 processor
818  do p = 0, npes - 1
819  if(grid1%tile(p) == tile1) then
820  if(in_box_nbr(i1_tmp(l), j1_tmp(l), grid1, p)) then
821  found = .true.
822  exit
823  endif
824  endif
825  enddo
826  !--- Then check if the exchange grid is on one of side 2 processor
827  if( found ) then
828  do p = 0, npes - 1
829  if(grid%tile(p) == tile2) then
830  if (in_box_nbr(i2_tmp(l), j2_tmp(l), grid, p)) then
831  pos = pos+1
832  i1(pos) = i1_tmp(l)
833  j1(pos) = j1_tmp(l)
834  i2(pos) = i2_tmp(l)
835  j2(pos) = j2_tmp(l)
836  area(pos) = area_tmp(l)
837  if(use_higher_order) then
838  di(pos) = di_tmp(l)
839  dj(pos) = dj_tmp(l)
840  endif
841  exit
842  endif
843  endif
844  enddo
845  endif
846  enddo
847 
848  deallocate(i1_tmp, i2_tmp, j1_tmp, j2_tmp, area_tmp)
849  if(use_higher_order) deallocate( di_tmp, dj_tmp)
850  iec = pos
851  if(iec .GE. isc) then
852  nxgrid_local = iec - isc + 1
853  else
854  nxgrid_local = 0
855  endif
856  else
857  nxgrid_local = 0
858  nxgrid_local_orig = 0
859  endif
860 
861  call mpp_clock_end(id_load_xgrid1)
862 
863  if(is_distribute) then
864  !--- Since the xgrid is distributed according to side 2 grid. Send all the xgrid to its own side 2.
865  !--- Also need to send the xgrid to its own side 1 for the reproducing ability between processor count.
866  !--- first find out number of points need to send to other pe and fill the send buffer.
867  nsend1(:) = 0; nrecv1(:) = 0
868  nsend2(:) = 0; nrecv2(:) = 0
869  ibuf1(:)= 0; ibuf2(:)= 0
870 
871  call mpp_clock_begin(id_load_xgrid2)
872  if(nxgrid_local>0) then
873  allocate( send_buffer(nxgrid_local * (nset1+nset2)) )
874  pos = 0
875  do p = 0, npes - 1
876  send_buffer_pos(p) = pos
877  if(grid%tile(p) == tile2) then
878  do l = isc, iec
879  if(in_box_nbr(i2(l), j2(l), grid, p) ) then
880  nsend2(p) = nsend2(p) + 1
881  send_buffer(pos+1) = i1(l)
882  send_buffer(pos+2) = j1(l)
883  send_buffer(pos+3) = i2(l)
884  send_buffer(pos+4) = j2(l)
885  send_buffer(pos+5) = area(l)
886  if(use_higher_order) then
887  send_buffer(pos+6) = di(l)
888  send_buffer(pos+7) = dj(l)
889  endif
890  if(scale_exist) send_buffer(pos+nset2) = scale(l)
891  pos = pos + nset2
892  endif
893  enddo
894  endif
895  if(grid1%tile(p) == tile1) then
896  do l = isc, iec
897  if(in_box_nbr(i1(l), j1(l), grid1, p)) then
898  nsend1(p) = nsend1(p) + 1
899  send_buffer(pos+1) = i1(l)
900  send_buffer(pos+2) = j1(l)
901  send_buffer(pos+3) = i2(l)
902  send_buffer(pos+4) = j2(l)
903  send_buffer(pos+5) = area(l)
904  if(use_higher_order) then
905  send_buffer(pos+6) = di(l)
906  send_buffer(pos+7) = dj(l)
907  endif
908  pos = pos + nset1
909  endif
910  enddo
911  endif
912  enddo
913  endif
914  call mpp_clock_end(id_load_xgrid2)
915 
916  !--- send the size of the data on side 1 to be sent over.
917  call mpp_clock_begin(id_load_xgrid3)
918 
919  if (do_alltoall) then
920  do p = 0, npes-1
921  ibuf1(2*p+1) = nsend1(p)
922  ibuf1(2*p+2) = nsend2(p)
923  enddo
924  call mpp_alltoall(ibuf1, 2, ibuf2, 2)
925  else
926  do n = 0, npes-1
927  p = mod(mypos+npes-n, npes)
928  if(.not. subset_rootpe(p)) cycle
929  call mpp_recv( ibuf2(2*p+1), glen=2, from_pe=pelist(p), block=.false., tag=comm_tag_1)
930  enddo
931 
932  if(nxgrid_local_orig>0) then
933  do n = 0, npes-1
934  p = mod(mypos+n, npes)
935  ibuf1(2*p+1) = nsend1(p)
936  ibuf1(2*p+2) = nsend2(p)
937  call mpp_send( ibuf1(2*p+1), plen=2, to_pe=pelist(p), tag=comm_tag_1)
938  enddo
939  endif
940  call mpp_sync_self(check=event_recv)
941  endif
942  do p = 0, npes-1
943  nrecv1(p) = ibuf2(2*p+1)
944  nrecv2(p) = ibuf2(2*p+2)
945  enddo
946 
947  if(.not. do_alltoall) call mpp_sync_self()
948  call mpp_clock_end(id_load_xgrid3)
949  call mpp_clock_begin(id_load_xgrid4)
950  pos = 0
951  do p = 0, npes - 1
952  recv_buffer_pos(p) = pos
953  pos = pos + nrecv1(p) * nset1 + nrecv2(p) * nset2
954  end do
955 
956  !--- now get the data
957  nxgrid1 = sum(nrecv1)
958  nxgrid2 = sum(nrecv2)
959  if(nxgrid1>0 .OR. nxgrid2>0) allocate(recv_buffer(nxgrid1*nset1+nxgrid2*nset2))
960 
961  if (do_alltoallv) then
962  ! Construct the send and receive counters
963  send_cnt(:) = nset1 * nsend1(:) + nset2 * nsend2(:)
964  recv_cnt(:) = nset1 * nrecv1(:) + nset2 * nrecv2(:)
965 
966  call mpp_alltoall(send_buffer, send_cnt, send_buffer_pos, &
967  recv_buffer, recv_cnt, recv_buffer_pos)
968  else
969  do n = 0, npes-1
970  p = mod(mypos+npes-n, npes)
971  nrecv = nrecv1(p)*nset1+nrecv2(p)*nset2
972  if(nrecv==0) cycle
973  pos = recv_buffer_pos(p)
974  call mpp_recv(recv_buffer(pos+1), glen=nrecv, from_pe=pelist(p), &
975  block=.false., tag=comm_tag_2)
976  end do
977 
978  do n = 0, npes-1
979  p = mod(mypos+n, npes)
980  nsend = nsend1(p)*nset1 + nsend2(p)*nset2
981  if(nsend==0) cycle
982  pos = send_buffer_pos(p)
983  call mpp_send(send_buffer(pos+1), plen=nsend, to_pe=pelist(p), &
984  tag=comm_tag_2)
985  end do
986  call mpp_sync_self(check=event_recv)
987  end if
988  call mpp_clock_end(id_load_xgrid4)
989  !--- unpack buffer.
990  if( nxgrid_local>0) then
991  deallocate(i1,j1,i2,j2,area)
992  endif
993 
994  allocate(i1(nxgrid2), j1(nxgrid2))
995  allocate(i2(nxgrid2), j2(nxgrid2))
996  allocate(area(nxgrid2))
997  allocate(i1_side1(nxgrid1), j1_side1(nxgrid1))
998  allocate(i2_side1(nxgrid1), j2_side1(nxgrid1))
999  allocate(area_side1(nxgrid1))
1000  if(use_higher_order) then
1001  if(nxgrid_local>0) deallocate(di,dj)
1002  allocate(di(nxgrid2), dj(nxgrid2))
1003  allocate(di_side1(nxgrid1), dj_side1(nxgrid1))
1004  endif
1005  if(scale_exist) then
1006  if(nxgrid_local>0)deallocate(scale)
1007  allocate(scale(nxgrid2))
1008  endif
1009  pos = 0
1010  l1 = 0; l2 = 0
1011  do p = 0,npes-1
1012  do n = 1, nrecv2(p)
1013  l2 = l2+1
1014  i1(l2) = recv_buffer(pos+1)
1015  j1(l2) = recv_buffer(pos+2)
1016  i2(l2) = recv_buffer(pos+3)
1017  j2(l2) = recv_buffer(pos+4)
1018  area(l2) = recv_buffer(pos+5)
1019  if(use_higher_order) then
1020  di(l2) = recv_buffer(pos+6)
1021  dj(l2) = recv_buffer(pos+7)
1022  endif
1023  if(scale_exist)scale(l2) = recv_buffer(pos+nset2)
1024  pos = pos + nset2
1025  enddo
1026  do n = 1, nrecv1(p)
1027  l1 = l1+1
1028  i1_side1(l1) = recv_buffer(pos+1)
1029  j1_side1(l1) = recv_buffer(pos+2)
1030  i2_side1(l1) = recv_buffer(pos+3)
1031  j2_side1(l1) = recv_buffer(pos+4)
1032  area_side1(l1) = recv_buffer(pos+5)
1033  if(use_higher_order) then
1034  di_side1(l1) = recv_buffer(pos+6)
1035  dj_side1(l1) = recv_buffer(pos+7)
1036  endif
1037  pos = pos + nset1
1038  enddo
1039  enddo
1040  call mpp_sync_self()
1041  if(allocated(send_buffer)) deallocate(send_buffer)
1042  if(allocated(recv_buffer)) deallocate(recv_buffer)
1043 
1044  else
1045  nxgrid1 = nxgrid
1046  nxgrid2 = nxgrid
1047  i1_side1 => i1; j1_side1 => j1
1048  i2_side1 => i2; j2_side1 => j2
1049  area_side1 => area
1050  if(use_higher_order) then
1051  di_side1 => di
1052  dj_side1 => dj
1053  endif
1054  endif
1055 
1056  call mpp_clock_begin(id_load_xgrid5)
1057 
1058 
1059  size_prev = grid%size
1060 
1061  if(grid%tile_me == tile2) then
1062  do l=1,nxgrid2
1063  if (in_box_me(i2(l), j2(l), grid) ) then
1064  grid%size = grid%size + 1
1065  ! exclude the area overlapped with parent grid
1066  if( grid1_id .NE. "ATM" .OR. tile1 .NE. tile_parent .OR. &
1067  .NOT. in_box(i1(l), j1(l), is_parent, ie_parent, js_parent, je_parent) ) then
1068  if(grid%is_ug) then
1069  lll = grid%l_index((j2(l)-1)*grid%im+i2(l))
1070  grid%area(lll,1) = grid%area(lll,1)+area(l)
1071  else
1072  grid%area(i2(l),j2(l)) = grid%area(i2(l),j2(l))+area(l)
1073  endif
1074  endif
1075  do p=0,xmap%npes-1
1076  if(grid1%tile(p) == tile1) then
1077  if (in_box_nbr(i1(l), j1(l), grid1, p)) then
1078  xmap%your1my2(p) = .true.
1079  end if
1080  end if
1081  end do
1082  end if
1083  end do
1084  end if
1085 
1086  if(grid%size > size_prev) then
1087  if(size_prev > 0) then ! need to extend data
1088  allocate(x_local(size_prev))
1089  x_local = grid%x
1090  if(ASSOCIATED(grid%x)) deallocate(grid%x)
1091  allocate( grid%x( grid%size ) )
1092  grid%x(1:size_prev) = x_local
1093  deallocate(x_local)
1094  else
1095  allocate( grid%x( grid%size ) )
1096  grid%x%di = 0.0; grid%x%dj = 0.0
1097  end if
1098  end if
1099 
1100  ll = size_prev
1101  if( grid%tile_me == tile2 ) then ! me is tile2
1102  do l=1,nxgrid2
1103  if (in_box_me(i2(l), j2(l), grid)) then
1104  ! insert in this grids cell pattern list and add area to side 2 area
1105  ll = ll + 1
1106  grid%x(ll)%i1 = i1(l); grid%x(ll)%i2 = i2(l)
1107  grid%x(ll)%j1 = j1(l); grid%x(ll)%j2 = j2(l)
1108  if(grid%is_ug) then
1109  grid%x(ll)%l2 = grid%l_index((j2(l)-1)*grid%im + i2(l))
1110  endif
1111 ! if(grid1%is_ug) then
1112 ! grid1%x(ll)%l1 = grid1%l_index((j1(l)-1)*grid1%im + i1(l))
1113 ! endif
1114  grid%x(ll)%tile = tile1
1115  grid%x(ll)%area = area(l)
1116  if(scale_exist) then
1117  grid%x(ll)%scale = scale(l)
1118  else
1119  grid%x(ll)%scale = 1.0
1120  endif
1121  if(use_higher_order) then
1122  grid%x(ll)%di = di(l)
1123  grid%x(ll)%dj = dj(l)
1124  end if
1125 
1126  if (make_exchange_reproduce) then
1127  do p=0,xmap%npes-1
1128  if(grid1%tile(p) == tile1) then
1129  if (in_box_nbr(i1(l), j1(l), grid1, p)) then
1130  grid%x(ll)%pe = p + xmap%root_pe
1131  end if
1132  end if
1133  end do
1134  end if ! make_exchange reproduce
1135  end if
1136  end do
1137  end if
1138 
1139  if(grid%id == xmap%grids(size(xmap%grids(:)))%id) then
1140  last_grid = .true.
1141  else
1142  last_grid = .false.
1143  endif
1144 
1145  size_repro = 0
1146  if(grid1%tile_me == tile1) then
1147  if(associated(iarray)) then
1148  nxgrid1_old = size(iarray(:))
1149  else
1150  nxgrid1_old = 0
1151  endif
1152 
1153  allocate(y2m1_pe(nxgrid1))
1154  if(.not. last_grid ) allocate(pos_s(0:xmap%npes-1))
1155  y2m1_pe = -1
1156  if(nxgrid1_old > 0) then
1157  do p=0,xmap%npes-1
1158  y2m1_size(p) = xmap%your2my1_size(p)
1159  enddo
1160  else
1161  y2m1_size = 0
1162  endif
1163 
1164  do l=1,nxgrid1
1165  if (in_box_me(i1_side1(l), j1_side1(l), grid1) ) then
1166  if(grid1%is_ug) then
1167  lll = grid1%l_index((j1_side1(l)-1)*grid1%im+i1_side1(l))
1168  grid1%area(lll,1) = grid1%area(lll,1) + area_side1(l)
1169  else
1170  grid1%area(i1_side1(l),j1_side1(l)) = grid1%area(i1_side1(l),j1_side1(l))+area_side1(l)
1171  endif
1172  do p=0,xmap%npes-1
1173  if (grid%tile(p) == tile2) then
1174  if (in_box_nbr(i2_side1(l), j2_side1(l), grid, p)) then
1175  xmap%your2my1(p) = .true.
1176  y2m1_pe(l) = p
1177  y2m1_size(p) = y2m1_size(p) + 1
1178  endif
1179  endif
1180  enddo
1181  size_repro = size_repro + 1
1182  endif
1183  enddo
1184  pos_x = 0
1185  do p = 1, npes-1
1186  pos_x(p) = pos_x(p-1) + y2m1_size(p-1)
1187  enddo
1188 
1189  if(.not. last_grid) pos_s(:) = pos_x(:)
1190 
1191  if(nxgrid1_old > 0) then
1192  y2m1_size(:) = xmap%your2my1_size(:)
1193  iarray2 => iarray
1194  jarray2 => jarray
1195  allocate(iarray(nxgrid1+nxgrid1_old), jarray(nxgrid1+nxgrid1_old))
1196  ! copy the i-j index
1197  do p=0,xmap%npes-1
1198  do n = 1, xmap%your2my1_size(p)
1199  iarray(pos_x(p)+n) = iarray2(pos_s(p)+n)
1200  jarray(pos_x(p)+n) = jarray2(pos_s(p)+n)
1201  enddo
1202  enddo
1203  deallocate(iarray2, jarray2)
1204  else
1205  allocate(iarray(nxgrid1), jarray(nxgrid1))
1206  iarray(:) = 0
1207  jarray(:) = 0
1208  y2m1_size(:) = 0
1209  endif
1210 
1211  do l=1,nxgrid1
1212  p = y2m1_pe(l)
1213  if(p<0) cycle
1214  found = .false.
1215  if(y2m1_size(p) > 0) then
1216  pos = pos_x(p)+y2m1_size(p)
1217  if( i1_side1(l) == iarray(pos) .AND. j1_side1(l) == jarray(pos) ) then
1218  found = .true.
1219  else
1220  !---may need to replace with a fast search algorithm
1221  do n = 1, y2m1_size(p)
1222  pos = pos_x(p)+n
1223  if(i1_side1(l) == iarray(pos) .AND. j1_side1(l) == jarray(pos)) then
1224  found = .true.
1225  exit
1226  endif
1227  enddo
1228  endif
1229  endif
1230  if( (.NOT. found) .OR. monotonic_exchange ) then
1231  y2m1_size(p) = y2m1_size(p)+1
1232  pos = pos_x(p)+y2m1_size(p)
1233  iarray(pos) = i1_side1(l)
1234  jarray(pos) = j1_side1(l)
1235  endif
1236  end do
1237  xmap%your2my1_size(:) = y2m1_size(:)
1238  deallocate(y2m1_pe)
1239  if(last_grid) then
1240  deallocate(iarray, jarray)
1241  if(allocated(pos_s)) deallocate(pos_s)
1242  end if
1243  end if
1244 
1245  if (grid1%tile_me == tile1 .and. size_repro > 0) then
1246  ll_repro = grid%size_repro
1247  grid%size_repro = ll_repro + size_repro
1248  if(ll_repro > 0) then ! extend data
1249  allocate(x_local(ll_repro))
1250  x_local = grid%x_repro
1251  if(ASSOCIATED(grid%x_repro)) deallocate(grid%x_repro)
1252  allocate( grid%x_repro(grid%size_repro ) )
1253  grid%x_repro(1:ll_repro) = x_local
1254  deallocate(x_local)
1255  else
1256  allocate( grid%x_repro( grid%size_repro ) )
1257  grid%x_repro%di = 0.0; grid%x_repro%dj = 0.0
1258  end if
1259  do l=1,nxgrid1
1260  if (in_box_me(i1_side1(l),j1_side1(l), grid1) ) then
1261  ll_repro = ll_repro + 1
1262  grid%x_repro(ll_repro)%i1 = i1_side1(l); grid%x_repro(ll_repro)%i2 = i2_side1(l)
1263  grid%x_repro(ll_repro)%j1 = j1_side1(l); grid%x_repro(ll_repro)%j2 = j2_side1(l)
1264  if(grid1%is_ug) then
1265  grid%x_repro(ll_repro)%l1 = grid1%l_index((j1_side1(l)-1)*grid1%im+i1_side1(l))
1266  endif
1267  if(grid%is_ug) then
1268 ! grid%x_repro(ll_repro)%l2 = grid%l_index((j2_side1(l)-1)*grid%im+i2_side1(l))
1269  endif
1270  grid%x_repro(ll_repro)%tile = tile1
1271  grid%x_repro(ll_repro)%area = area_side1(l)
1272  if(use_higher_order) then
1273  grid%x_repro(ll_repro)%di = di_side1(l)
1274  grid%x_repro(ll_repro)%dj = dj_side1(l)
1275  end if
1276 
1277  do p=0,xmap%npes-1
1278  if(grid%tile(p) == tile2) then
1279  if (in_box_nbr(i2_side1(l), j2_side1(l), grid, p)) then
1280  grid%x_repro(ll_repro)%pe = p + xmap%root_pe
1281  end if
1282  end if
1283  end do
1284  end if ! make_exchange_reproduce
1285  end do
1286  end if
1287 
1288  deallocate(i1, j1, i2, j2, area)
1289  if(use_higher_order) deallocate(di, dj)
1290  if(scale_exist) deallocate(scale)
1291  if(is_distribute) then
1292  deallocate(i1_side1, j1_side1, i2_side1, j2_side1, area_side1)
1293  if(use_higher_order) deallocate(di_side1, dj_side1)
1294  endif
1295 
1296  i1=>null(); j1=>null(); i2=>null(); j2=>null()
1297  call mpp_clock_end(id_load_xgrid5)
1298 
1299 
1300 
1301 end subroutine load_xgrid
1302 
1303 !#######################################################################
1304 !
1305 ! get_grid - read the center point of the grid from grid_spec.nc.
1306 ! - only the grid at the side 1 is needed, so we only read
1307 ! - atm and land grid
1308 !
1309 !
1310 
1311 subroutine get_grid(grid, grid_id, grid_file, grid_version)
1312  type(grid_type), intent(inout) :: grid
1313  character(len=3), intent(in) :: grid_id
1314  character(len=*), intent(in) :: grid_file
1315  integer, intent(in) :: grid_version
1316 
1317  real, dimension(grid%im) :: lonb
1318  real, dimension(grid%jm) :: latb
1319  real, allocatable :: tmpx(:,:), tmpy(:,:)
1320  real :: d2r
1321  integer :: is, ie, js, je, nlon, nlat, siz(4), i, j
1322  integer :: start(4), nread(4), isc2, iec2, jsc2, jec2
1323 
1324  d2r = pi/180.0
1325 
1326  call mpp_get_compute_domain(grid%domain, is, ie, js, je)
1327 
1328  select case(grid_version)
1329  case(version1)
1330  allocate(grid%lon(grid%im), grid%lat(grid%jm))
1331  if(grid_id == 'ATM') then
1332  call read_data(grid_file, 'xta', lonb)
1333  call read_data(grid_file, 'yta', latb)
1334 
1335  if(.not. allocated(area_atm_model)) then
1336  allocate(area_atm_model(is:ie, js:je))
1337  call get_area_elements(grid_file, 'AREA_ATM_MODEL', grid%domain, area_atm_model)
1338  endif
1339  if(.not. allocated(area_atm_sphere)) then
1340  allocate(area_atm_sphere(is:ie, js:je))
1341  call get_area_elements(grid_file, 'AREA_ATM', grid%domain, area_atm_sphere)
1342  endif
1343  else if(grid_id == 'LND') then
1344  call read_data(grid_file, 'xtl', lonb)
1345  call read_data(grid_file, 'ytl', latb)
1346  if(.not. allocated(area_lnd_model)) then
1347  allocate(area_lnd_model(is:ie, js:je))
1348  call get_area_elements(grid_file, 'AREA_LND_MODEL', grid%domain, area_lnd_model)
1349  endif
1350  if(.not. allocated(area_lnd_sphere)) then
1351  allocate(area_lnd_sphere(is:ie, js:je))
1352  call get_area_elements(grid_file, 'AREA_LND', grid%domain, area_lnd_sphere)
1353  endif
1354  else if(grid_id == 'OCN' ) then
1355  if(.not. allocated(area_ocn_sphere)) then
1356  allocate(area_ocn_sphere(is:ie, js:je))
1357  call get_area_elements(grid_file, 'AREA_OCN', grid%domain, area_ocn_sphere)
1358  endif
1359  endif
1360  !--- second order remapping suppose second order
1361  if(grid_id == 'LND' .or. grid_id == 'ATM') then
1362  grid%lon = lonb * d2r
1363  grid%lat = latb * d2r
1364  endif
1365  grid%is_latlon = .true.
1366  case(version2)
1367  call field_size(grid_file, 'area', siz)
1368  nlon = siz(1); nlat = siz(2)
1369  if( mod(nlon,2) .NE. 0) call error_mesg('xgrid_mod', &
1370  'flux_exchange_mod: atmos supergrid longitude size can not be divided by 2', fatal)
1371  if( mod(nlat,2) .NE. 0) call error_mesg('xgrid_mod', &
1372  'flux_exchange_mod: atmos supergrid latitude size can not be divided by 2', fatal)
1373  nlon = nlon/2
1374  nlat = nlat/2
1375  if(nlon .NE. grid%im .OR. nlat .NE. grid%jm) call error_mesg('xgrid_mod', &
1376  'grid size in tile_file does not match the global grid size', fatal)
1377 
1378  if( grid_id == 'LND' .or. grid_id == 'ATM' .or. grid_id == 'WAV' ) then
1379  isc2 = 2*grid%is_me-1; iec2 = 2*grid%ie_me+1
1380  jsc2 = 2*grid%js_me-1; jec2 = 2*grid%je_me+1
1381  allocate(tmpx(isc2:iec2, jsc2:jec2) )
1382  allocate(tmpy(isc2:iec2, jsc2:jec2) )
1383  start = 1; nread = 1
1384  start(1) = isc2; nread(1) = iec2 - isc2 + 1
1385  start(2) = jsc2; nread(2) = jec2 - jsc2 + 1
1386  call read_data(grid_file, 'x', tmpx, start, nread, no_domain=.true.)
1387  call read_data(grid_file, 'y', tmpy, start, nread, no_domain=.true.)
1388  if(is_lat_lon(tmpx, tmpy) ) then
1389  deallocate(tmpx, tmpy)
1390  start = 1; nread = 1
1391  start(2) = 2; nread(1) = nlon*2+1
1392  allocate(tmpx(nlon*2+1, 1), tmpy(1, nlat*2+1))
1393  call read_data(grid_file, "x", tmpx, start, nread, no_domain=.true.)
1394  allocate(grid%lon(grid%im), grid%lat(grid%jm))
1395  do i = 1, grid%im
1396  grid%lon(i) = tmpx(2*i,1) * d2r
1397  end do
1398  start = 1; nread = 1
1399  start(1) = 2; nread(2) = nlat*2+1
1400  call read_data(grid_file, "y", tmpy, start, nread, no_domain=.true.)
1401  do j = 1, grid%jm
1402  grid%lat(j) = tmpy(1, 2*j) * d2r
1403  end do
1404  grid%is_latlon = .true.
1405  else
1406  allocate(grid%geolon(grid%isd_me:grid%ied_me, grid%jsd_me:grid%jed_me))
1407  allocate(grid%geolat(grid%isd_me:grid%ied_me, grid%jsd_me:grid%jed_me))
1408  grid%geolon = 1e10
1409  grid%geolat = 1e10
1410  !--- area_ocn_sphere, area_lnd_sphere, area_atm_sphere is not been defined.
1411  do j = grid%js_me,grid%je_me
1412  do i = grid%is_me,grid%ie_me
1413  grid%geolon(i, j) = tmpx(i*2,j*2)*d2r
1414  grid%geolat(i, j) = tmpy(i*2,j*2)*d2r
1415  end do
1416  end do
1417  call mpp_update_domains(grid%geolon, grid%domain)
1418  call mpp_update_domains(grid%geolat, grid%domain)
1419  grid%is_latlon = .false.
1420  end if
1421  deallocate(tmpx, tmpy)
1422  end if
1423  end select
1424 
1425  return
1426 
1427 end subroutine get_grid
1428 
1429 !#######################################################################
1430 ! Read the area elements from NetCDF file
1431 subroutine get_area_elements(file, name, domain, data)
1432  character(len=*), intent(in) :: file
1433  character(len=*), intent(in) :: name
1434  type(domain2d), intent(in) :: domain
1435  real, intent(out) :: data(:,:)
1436 
1437  if(field_exist(file, name)) then
1438  call read_data(file, name, data, domain)
1439  else
1440  call error_mesg('xgrid_mod', 'no field named '//trim(name)//' in grid file '//trim(file)// &
1441  ' Will set data to negative values...', note)
1442  ! area elements no present in grid_spec file, set to negative values....
1443  data = -1.0
1444  endif
1445 
1446 end subroutine get_area_elements
1447 
1448 !#######################################################################
1449 ! Read the OCN model area elements from NetCDF file
1450 ! <SUBROUTINE NAME="get_ocean_model_area_elements">
1451 
1452 ! <OVERVIEW>
1453 ! Read Ocean area element data.
1454 ! </OVERVIEW>
1455 ! <DESCRIPTION>
1456 ! If available in the NetCDF file, this routine will read the
1457 ! AREA_OCN_MODEL field and load the data into global AREA_OCN_MODEL.
1458 ! If not available, then the array AREA_OCN_MODEL will be left
1459 ! unallocated. Must be called by all PEs.
1460 ! </DESCRIPTION>
1461 ! <TEMPLATE>
1462 ! call get_ocean_model_area_elements(ocean_domain, grid_file)
1463 ! </TEMPLATE>
1464 
1465 ! <IN NAME="ocean_domain" TYPE="type(Domain2d)"> </IN>
1466 ! <IN NAME="grid_file" TYPE="character(len=*)" > </IN>
1467 subroutine get_ocean_model_area_elements(domain, grid_file)
1469  type(domain2d), intent(in) :: domain
1470  character(len=*), intent(in) :: grid_file
1471  integer :: is, ie, js, je
1472 
1473  if(allocated(area_ocn_model)) return
1474 
1475  call mpp_get_compute_domain(domain, is, ie, js, je)
1476  ! allocate even if ie<is, ... in which case the array will have zero size
1477  ! but will still return .T. for allocated(...)
1478  allocate(area_ocn_model(is:ie, js:je))
1479  if(ie < is .or. je < js ) return
1480 
1481 
1482  if(field_exist(grid_file, 'AREA_OCN_MODEL') )then
1483  call read_data(grid_file, 'AREA_OCN_MODEL', area_ocn_model, domain)
1484  else
1485  deallocate(area_ocn_model)
1486  endif
1487 
1488 
1489 end subroutine get_ocean_model_area_elements
1490 ! </SUBROUTINE>
1491 !#######################################################################
1492 
1493 ! <SUBROUTINE NAME="setup_xmap">
1494 
1495 ! <OVERVIEW>
1496 ! Sets up exchange grid connectivity using grid specification file and
1497 ! processor domain decomposition.
1498 ! </OVERVIEW>
1499 ! <DESCRIPTION>
1500 ! Sets up exchange grid connectivity using grid specification file and
1501 ! processor domain decomposition. Initializes xmap.
1502 ! </DESCRIPTION>
1503 ! <TEMPLATE>
1504 ! call setup_xmap(xmap, grid_ids, grid_domains, grid_file, atm_grid)
1505 ! </TEMPLATE>
1506 
1507 ! <IN NAME="grid_ids" TYPE="character(len=3)" DIM="(:)"> </IN>
1508 ! <IN NAME="grid_domains" TYPE="type(Domain2d)" DIM="(:)"> </IN>
1509 ! <IN NAME="grid_file" TYPE="character(len=*)" > </IN>
1510 ! <IN NAME="atmos_grid" TYPE="type(grid_box_type),optional" > </IN>
1511 ! <OUT NAME="xmap" TYPE="xmap_type" > </OUT>
1512 
1513 subroutine setup_xmap(xmap, grid_ids, grid_domains, grid_file, atm_grid, lnd_ug_domain)
1514  type(xmap_type), intent(inout) :: xmap
1515  character(len=3), dimension(:), intent(in ) :: grid_ids
1516  type(domain2d), dimension(:), intent(in ) :: grid_domains
1517  character(len=*), intent(in ) :: grid_file
1518  type(grid_box_type), optional, intent(in ) :: atm_grid
1519  type(domainug), optional, intent(in ) :: lnd_ug_domain
1520 
1521  integer :: g, p, send_size, recv_size, i, siz(4)
1522  integer :: unit, nxgrid_file, i1, i2, i3, tile1, tile2, j
1523  integer :: nxc, nyc, out_unit
1524  type(grid_type), pointer, save :: grid =>null(), grid1 =>null()
1525  real, dimension(3) :: xxx
1526  real, dimension(:,:), allocatable :: check_data
1527  real, dimension(:,:,:), allocatable :: check_data_3d
1528  real, allocatable :: tmp_2d(:,:), tmp_3d(:,:,:)
1529  character(len=256) :: xgrid_file, xgrid_name
1530  character(len=256) :: tile_file, mosaic_file
1531  character(len=256) :: mosaic1, mosaic2, contact
1532  character(len=256) :: tile1_name, tile2_name
1533  character(len=256), allocatable :: tile1_list(:), tile2_list(:)
1534  integer :: npes, npes2
1535  integer, allocatable :: pelist(:)
1536  type(domain2d), save :: domain2
1537  logical :: use_higher_order = .false.
1538  integer :: lnd_ug_id, l
1539  integer, allocatable :: grid_index(:)
1540 
1541  call mpp_clock_begin(id_setup_xmap)
1542 
1543  if(interp_method .ne. 'first_order') use_higher_order = .true.
1544 
1545  out_unit = stdout()
1546  xmap%me = mpp_pe()
1547  xmap%npes = mpp_npes()
1548  xmap%root_pe = mpp_root_pe()
1549 
1550  allocate( xmap%grids(1:size(grid_ids(:))) )
1551 
1552  allocate ( xmap%your1my2(0:xmap%npes-1), xmap%your2my1(0:xmap%npes-1) )
1553  allocate ( xmap%your2my1_size(0:xmap%npes-1) )
1554 
1555  xmap%your1my2 = .false.; xmap%your2my1 = .false.;
1556  xmap%your2my1_size = 0
1557 
1558 ! check the exchange grid file version to be used by checking the field in the file
1559  if(field_exist(grid_file, "AREA_ATMxOCN" ) ) then
1560  xmap%version = version1
1561  else if(field_exist(grid_file, "ocn_mosaic_file" ) ) then
1562  xmap%version = version2
1563  else
1564  call error_mesg('xgrid_mod', 'both AREA_ATMxOCN and ocn_mosaic_file does not exist in '//trim(grid_file), fatal)
1565  end if
1566 
1567  if(xmap%version==version1) then
1568  call error_mesg('xgrid_mod', 'reading exchange grid information from grid spec file', note)
1569  else
1570  call error_mesg('xgrid_mod', 'reading exchange grid information from mosaic grid file', note)
1571  end if
1572 
1573  ! check to see the id of lnd.
1574  lnd_ug_id = 0
1575  if(present(lnd_ug_domain)) then
1576  do g=1,size(grid_ids(:))
1577  if(grid_ids(g) == 'LND') lnd_ug_id = g
1578  enddo
1579  endif
1580 
1581  call mpp_clock_begin(id_load_xgrid)
1582  do g=1,size(grid_ids(:))
1583  grid => xmap%grids(g)
1584  if (g==1) grid1 => xmap%grids(g)
1585  grid%id = grid_ids(g)
1586  grid%domain = grid_domains(g)
1587  grid%on_this_pe = mpp_domain_is_initialized(grid_domains(g))
1588  allocate ( grid%is(0:xmap%npes-1), grid%ie(0:xmap%npes-1) )
1589  allocate ( grid%js(0:xmap%npes-1), grid%je(0:xmap%npes-1) )
1590  allocate ( grid%tile(0:xmap%npes-1) )
1591  grid%npes = 0
1592  grid%ni = 0
1593  grid%nj = 0
1594  grid%is = 0
1595  grid%ie = -1
1596  grid%js = 0
1597  grid%je = -1
1598  grid%tile = -1
1599 
1600  select case(xmap%version)
1601  case(version1)
1602  grid%ntile = 1
1603  case(version2)
1604  call read_data(grid_file, lowercase(grid_ids(g))//'_mosaic_file', mosaic_file)
1605  grid%ntile = get_mosaic_ntiles('INPUT/'//trim(mosaic_file))
1606  end select
1607 
1608  if( g == 1 .AND. grid_ids(1) == 'ATM' ) then
1609  if( .NOT. grid%on_this_pe ) call error_mesg('xgrid_mod', 'ATM domain is not defined on some processor' ,fatal)
1610  endif
1611  grid%npes = mpp_get_domain_npes(grid%domain)
1612  if( xmap%npes > grid%npes .AND. g == 1 .AND. grid_ids(1) == 'ATM' ) then
1613  call mpp_broadcast_domain(grid%domain, domain2)
1614  else if(xmap%npes > grid%npes) then
1615  call mpp_broadcast_domain(grid%domain)
1616  grid%npes = mpp_get_domain_npes(grid%domain)
1617  endif
1618 
1619  npes = grid%npes
1620  allocate(grid%pelist(0:npes-1))
1621  call mpp_get_domain_pelist(grid%domain, grid%pelist)
1622  grid%root_pe = mpp_get_domain_root_pe(grid%domain)
1623 
1624  call mpp_get_data_domain(grid%domain, grid%isd_me, grid%ied_me, grid%jsd_me, grid%jed_me, &
1625  xsize=grid%nxd_me, ysize=grid%nyd_me)
1626  call mpp_get_global_domain(grid%domain, xsize=grid%ni, ysize=grid%nj)
1627 
1628  if( grid%root_pe == xmap%root_pe ) then
1629  call mpp_get_compute_domains(grid%domain, xbegin=grid%is(0:npes-1), xend=grid%ie(0:npes-1), &
1630  ybegin=grid%js(0:npes-1), yend=grid%je(0:npes-1) )
1631  call mpp_get_tile_list(grid%domain, grid%tile(0:npes-1))
1632  if( xmap%npes > npes .AND. g == 1 .AND. grid_ids(1) == 'ATM' ) then
1633  call mpp_get_compute_domains(domain2, xbegin=grid%is(npes:xmap%npes-1), xend=grid%ie(npes:xmap%npes-1), &
1634  ybegin=grid%js(npes:xmap%npes-1), yend=grid%je(npes:xmap%npes-1) )
1635  call mpp_get_tile_list(domain2, grid%tile(npes:xmap%npes-1))
1636  endif
1637  else
1638  npes2 = xmap%npes-npes
1639  call mpp_get_compute_domains(domain2, xbegin=grid%is(0:npes2-1), xend=grid%ie(0:npes2-1), &
1640  ybegin=grid%js(0:npes2-1), yend=grid%je(0:npes2-1) )
1641  call mpp_get_compute_domains(grid%domain, xbegin=grid%is(npes2:xmap%npes-1), xend=grid%ie(npes2:xmap%npes-1), &
1642  ybegin=grid%js(npes2:xmap%npes-1), yend=grid%je(npes2:xmap%npes-1) )
1643  call mpp_get_tile_list(domain2, grid%tile(0:npes2-1))
1644  call mpp_get_tile_list(grid%domain, grid%tile(npes2:xmap%npes-1))
1645  endif
1646  if( xmap%npes > grid%npes .AND. g == 1 .AND. grid_ids(1) == 'ATM' ) then
1647  call mpp_deallocate_domain(domain2)
1648  endif
1649  npes = grid%npes
1650  if( g == 1 .AND. grid_ids(1) == 'ATM' ) npes = xmap%npes
1651  do p = 0, npes-1
1652  if(grid%tile(p) > grid%ntile .or. grid%tile(p) < 1) call error_mesg('xgrid_mod', &
1653  'tile id should between 1 and ntile', fatal)
1654  end do
1655 
1656  grid%im = grid%ni
1657  grid%jm = grid%nj
1658  call mpp_max(grid%ni)
1659  call mpp_max(grid%nj)
1660 
1661  grid%is_me => grid%is(xmap%me-xmap%root_pe); grid%ie_me => grid%ie(xmap%me-xmap%root_pe)
1662  grid%js_me => grid%js(xmap%me-xmap%root_pe); grid%je_me => grid%je(xmap%me-xmap%root_pe)
1663  grid%nxc_me = grid%ie_me - grid%is_me + 1
1664  grid%nyc_me = grid%je_me - grid%js_me + 1
1665  grid%tile_me => grid%tile(xmap%me-xmap%root_pe)
1666 
1667  grid%km = 1
1668  grid%is_ug = .false.
1669  !--- setup for land unstructure grid
1670  if( g == lnd_ug_id ) then
1671  if(xmap%version == version1) call error_mesg('xgrid_mod', &
1672  'does not support unstructured grid for VERSION1 grid' ,fatal)
1673  grid%is_ug = .true.
1674  grid%ug_domain = lnd_ug_domain
1675  allocate ( grid%ls(0:xmap%npes-1), grid%le(0:xmap%npes-1) )
1676  allocate ( grid%gs(0:xmap%npes-1), grid%ge(0:xmap%npes-1) )
1677  grid%ls = 0
1678  grid%le = -1
1679  grid%gs = 0
1680  grid%ge = -1
1681  if(xmap%npes > grid%npes) then
1682  call mpp_broadcast_domain(grid%ug_domain)
1683  endif
1684  call mpp_get_ug_compute_domains(grid%ug_domain, begin=grid%ls(0:npes-1), end=grid%le(0:npes-1) )
1685  call mpp_get_ug_domains_index(grid%ug_domain, grid%gs(0:npes-1), grid%ge(0:npes-1) )
1686  call mpp_get_ug_domain_tile_list(grid%ug_domain, grid%tile(0:npes-1))
1687  grid%ls_me => grid%ls(xmap%me-xmap%root_pe); grid%le_me => grid%le(xmap%me-xmap%root_pe)
1688  grid%gs_me => grid%gs(xmap%me-xmap%root_pe); grid%ge_me => grid%ge(xmap%me-xmap%root_pe)
1689  grid%tile_me => grid%tile(xmap%me-xmap%root_pe)
1690  grid%nxl_me = grid%le_me - grid%ls_me + 1
1691  allocate(grid%l_index(grid%gs_me:grid%ge_me))
1692  allocate(grid_index(grid%ls_me:grid%le_me))
1693  call mpp_get_ug_domain_grid_index(grid%ug_domain, grid_index)
1694 
1695  grid%l_index = 0
1696  do l = grid%ls_me,grid%le_me
1697  grid%l_index(grid_index(l)) = l
1698  enddo
1699 
1700  if( grid%on_this_pe ) then
1701  allocate( grid%area (grid%ls_me:grid%le_me,1) )
1702  allocate( grid%area_inv(grid%ls_me:grid%le_me,1) )
1703  grid%area = 0.0
1704  grid%size = 0
1705  grid%size_repro = 0
1706  endif
1707  else if( grid%on_this_pe ) then
1708  allocate( grid%area (grid%is_me:grid%ie_me, grid%js_me:grid%je_me) )
1709  allocate( grid%area_inv(grid%is_me:grid%ie_me, grid%js_me:grid%je_me) )
1710  grid%area = 0.0
1711  grid%size = 0
1712  grid%size_repro = 0
1713  endif
1714 
1715  ! get the center point of the grid box
1716  if(.not. grid%is_ug) then
1717  select case(xmap%version)
1718  case(version1)
1719  if( grid%npes .NE. xmap%npes ) then
1720  call error_mesg('xgrid_mod', .NE.' grid%npes xmap%npes ', fatal)
1721  endif
1722  call get_grid(grid, grid_ids(g), grid_file, xmap%version)
1723  case(version2)
1724  allocate(pelist(0:xmap%npes-1))
1725  call mpp_get_current_pelist(pelist)
1726  if( grid%on_this_pe ) then
1727  call mpp_set_current_pelist(grid%pelist)
1728  call get_mosaic_tile_grid(tile_file, 'INPUT/'//trim(mosaic_file), grid%domain)
1729  call get_grid(grid, grid_ids(g), tile_file, xmap%version)
1730  endif
1731  call mpp_set_current_pelist(pelist)
1732  deallocate(pelist)
1733  ! read the contact information from mosaic_file to check if atmosphere is nested model
1734  if( g == 1 .AND. grid_ids(1) == 'ATM' ) then
1735  nnest = get_nest_contact('INPUT/'//trim(mosaic_file), tile_nest, tile_parent, is_nest, &
1737  endif
1738  end select
1739 
1740  if( use_higher_order .AND. grid%id == 'ATM') then
1741  if( nnest > 0 ) call error_mesg('xgrid_mod', 'second_order is not supported for nested coupler', fatal)
1742  if( grid%is_latlon ) then
1743  call mpp_modify_domain(grid%domain, grid%domain_with_halo, whalo=1, ehalo=1, shalo=1, nhalo=1)
1744  call mpp_get_data_domain(grid%domain_with_halo, grid%isd_me, grid%ied_me, grid%jsd_me, grid%jed_me, &
1745  xsize=grid%nxd_me, ysize=grid%nyd_me)
1746  else
1747  if(.NOT. present(atm_grid)) call error_mesg('xgrid_mod', &
1748  'when first grid is "ATM", atm_grid should be present', fatal)
1749  if(grid%is_me-grid%isd_me .NE. 1 .or. grid%ied_me-grid%ie_me .NE. 1 .or. &
1750  grid%js_me-grid%jsd_me .NE. 1 .or. grid%jed_me-grid%je_me .NE. 1 ) call error_mesg( &
1751  'xgrid_mod', 'for non-latlon grid (cubic grid), the halo size should be 1 in all four direction', fatal)
1752  if(.NOT.( ASSOCIATED(atm_grid%dx) .AND. ASSOCIATED(atm_grid%dy) .AND. ASSOCIATED(atm_grid%edge_w) .AND. &
1753  ASSOCIATED(atm_grid%edge_e) .AND. ASSOCIATED(atm_grid%edge_s) .AND. ASSOCIATED(atm_grid%edge_n) .AND. &
1754  ASSOCIATED(atm_grid%en1) .AND. ASSOCIATED(atm_grid%en2) .AND. ASSOCIATED(atm_grid%vlon) .AND. &
1755  ASSOCIATED(atm_grid%vlat) ) ) call error_mesg( 'xgrid_mod', &
1756  'for non-latlon grid (cubic grid), all the fields in atm_grid data type should be allocated', fatal)
1757  nxc = grid%ie_me - grid%is_me + 1
1758  nyc = grid%je_me - grid%js_me + 1
1759  if(size(atm_grid%dx,1) .NE. nxc .OR. size(atm_grid%dx,2) .NE. nyc+1) &
1760  call error_mesg('xgrid_mod', 'incorrect dimension size of atm_grid%dx', fatal)
1761  if(size(atm_grid%dy,1) .NE. nxc+1 .OR. size(atm_grid%dy,2) .NE. nyc) &
1762  call error_mesg('xgrid_mod', 'incorrect dimension sizeof atm_grid%dy', fatal)
1763  if(size(atm_grid%area,1) .NE. nxc .OR. size(atm_grid%area,2) .NE. nyc) &
1764  call error_mesg('xgrid_mod', 'incorrect dimension size of atm_grid%area', fatal)
1765  if(size(atm_grid%edge_w(:)) .NE. nyc+1 .OR. size(atm_grid%edge_e(:)) .NE. nyc+1) &
1766  call error_mesg('xgrid_mod', 'incorrect dimension size of atm_grid%edge_w/edge_e', fatal)
1767  if(size(atm_grid%edge_s(:)) .NE. nxc+1 .OR. size(atm_grid%edge_n(:)) .NE. nxc+1) &
1768  call error_mesg('xgrid_mod', 'incorrect dimension size of atm_grid%edge_s/edge_n', fatal)
1769  if(size(atm_grid%en1,1) .NE. 3 .OR. size(atm_grid%en1,2) .NE. nxc .OR. size(atm_grid%en1,3) .NE. nyc+1) &
1770  call error_mesg( 'xgrid_mod', 'incorrect dimension size of atm_grid%en1', fatal)
1771  if(size(atm_grid%en2,1) .NE. 3 .OR. size(atm_grid%en2,2) .NE. nxc+1 .OR. size(atm_grid%en2,3) .NE. nyc) &
1772  call error_mesg( 'xgrid_mod', 'incorrect dimension size of atm_grid%en2', fatal)
1773  if(size(atm_grid%vlon,1) .NE. 3 .OR. size(atm_grid%vlon,2) .NE. nxc .OR. size(atm_grid%vlon,3) .NE. nyc) &
1774  call error_mesg('xgrid_mod', 'incorrect dimension size of atm_grid%vlon', fatal)
1775  if(size(atm_grid%vlat,1) .NE. 3 .OR. size(atm_grid%vlat,2) .NE. nxc .OR. size(atm_grid%vlat,3) .NE. nyc) &
1776  call error_mesg('xgrid_mod', 'incorrect dimension size of atm_grid%vlat', fatal)
1777  allocate(grid%box%dx (grid%is_me:grid%ie_me, grid%js_me:grid%je_me+1 ))
1778  allocate(grid%box%dy (grid%is_me:grid%ie_me+1, grid%js_me:grid%je_me ))
1779  allocate(grid%box%area (grid%is_me:grid%ie_me, grid%js_me:grid%je_me ))
1780  allocate(grid%box%edge_w(grid%js_me:grid%je_me+1))
1781  allocate(grid%box%edge_e(grid%js_me:grid%je_me+1))
1782  allocate(grid%box%edge_s(grid%is_me:grid%ie_me+1))
1783  allocate(grid%box%edge_n(grid%is_me:grid%ie_me+1))
1784  allocate(grid%box%en1 (3, grid%is_me:grid%ie_me, grid%js_me:grid%je_me+1 ))
1785  allocate(grid%box%en2 (3, grid%is_me:grid%ie_me+1, grid%js_me:grid%je_me ))
1786  allocate(grid%box%vlon (3, grid%is_me:grid%ie_me, grid%js_me:grid%je_me ))
1787  allocate(grid%box%vlat (3, grid%is_me:grid%ie_me, grid%js_me:grid%je_me ))
1788  grid%box%dx = atm_grid%dx
1789  grid%box%dy = atm_grid%dy
1790  grid%box%area = atm_grid%area
1791  grid%box%edge_w = atm_grid%edge_w
1792  grid%box%edge_e = atm_grid%edge_e
1793  grid%box%edge_s = atm_grid%edge_s
1794  grid%box%edge_n = atm_grid%edge_n
1795  grid%box%en1 = atm_grid%en1
1796  grid%box%en2 = atm_grid%en2
1797  grid%box%vlon = atm_grid%vlon
1798  grid%box%vlat = atm_grid%vlat
1799  end if
1800  end if
1801  end if
1802 
1803  if (g>1) then
1804  if(grid%on_this_pe) then
1805  if(grid%is_ug) then
1806  allocate( grid%frac_area(grid%ls_me:grid%le_me, 1, grid%km) )
1807  else
1808  allocate( grid%frac_area(grid%is_me:grid%ie_me, grid%js_me:grid%je_me, grid%km) )
1809  endif
1810  grid%frac_area = 1.0
1811  endif
1812 
1813  ! load exchange cells, sum grid cell areas, set your1my2/your2my1
1814  select case(xmap%version)
1815  case(version1)
1816  call load_xgrid (xmap, grid, grid_file, grid_ids(1), grid_ids(g), 1, 1, use_higher_order)
1817  case(version2)
1818  select case(grid_ids(1))
1819  case( 'ATM' )
1820  xgrid_name = 'a'
1821  case( 'LND' )
1822  xgrid_name = 'l'
1823  case( 'WAV' )
1824  xgrid_name = 'w'
1825  case default
1826  call error_mesg('xgrid_mod', 'grid_ids(1) should be ATM, LND or WAV', fatal)
1827  end select
1828  select case(grid_ids(g))
1829  case( 'LND' )
1830  xgrid_name = trim(xgrid_name)//'Xl_file'
1831  case( 'OCN' )
1832  xgrid_name = trim(xgrid_name)//'Xo_file'
1833  case( 'WAV' )
1834  xgrid_name = trim(xgrid_name)//'Xw_file'
1835  case default
1836  call error_mesg('xgrid_mod', 'grid_ids(g) should be LND, OCN or WAV', fatal)
1837  end select
1838  ! get the tile list for each mosaic
1839  call read_data(grid_file, lowercase(grid_ids(1))//'_mosaic_file', mosaic1)
1840  call read_data(grid_file, lowercase(grid_ids(g))//'_mosaic_file', mosaic2)
1841  mosaic1 = 'INPUT/'//trim(mosaic1)
1842  mosaic2 = 'INPUT/'//trim(mosaic2)
1843  allocate(tile1_list(grid1%ntile), tile2_list(grid%ntile) )
1844  do j = 1, grid1%ntile
1845  call read_data(mosaic1, 'gridtiles', tile1_list(j), level=j)
1846  end do
1847  do j = 1, grid%ntile
1848  call read_data(mosaic2, 'gridtiles', tile2_list(j), level=j)
1849  end do
1850  if(field_exist(grid_file, xgrid_name)) then
1851  call field_size(grid_file, xgrid_name, siz)
1852  nxgrid_file = siz(2)
1853  ! loop through all the exchange grid file
1854  do i = 1, nxgrid_file
1855  call read_data(grid_file, xgrid_name, xgrid_file, level = i)
1856  xgrid_file = 'INPUT/'//trim(xgrid_file)
1857  if( .NOT. file_exist(xgrid_file) )call error_mesg('xgrid_mod', &
1858  'file '//trim(xgrid_file)//' does not exist, check your xgrid file.', fatal)
1859 
1860  ! find the tile number of side 1 and side 2 mosaic, which is contained in field contact
1861  call read_data(xgrid_file, "contact", contact)
1862  i1 = index(contact, ":")
1863  i2 = index(contact, "::")
1864  i3 = index(contact, ":", back=.true. )
1865  if(i1 == 0 .OR. i2 == 0) call error_mesg('xgrid_mod', &
1866  'field contact in file '//trim(xgrid_file)//' should contains ":" and "::" ', fatal)
1867  if(i1 == i3) call error_mesg('xgrid_mod', &
1868  'field contact in file '//trim(xgrid_file)//' should contains two ":"', fatal)
1869  tile1_name = contact(i1+1:i2-1)
1870  tile2_name = contact(i3+1:len_trim(contact))
1871  tile1 = 0; tile2 = 0
1872  do j = 1, grid1%ntile
1873  if( tile1_name == tile1_list(j) ) then
1874  tile1 = j
1875  exit
1876  end if
1877  end do
1878  do j = 1, grid%ntile
1879  if( tile2_name == tile2_list(j) ) then
1880  tile2 = j
1881  exit
1882  end if
1883  end do
1884  if(tile1 == 0) call error_mesg('xgrid_mod', &
1885  trim(tile1_name)//' is not a tile of mosaic '//trim(mosaic1), fatal)
1886  if(tile2 == 0) call error_mesg('xgrid_mod', &
1887  trim(tile2_name)//' is not a tile of mosaic '//trim(mosaic2), fatal)
1888 
1889  call load_xgrid (xmap, grid, xgrid_file, grid_ids(1), grid_ids(g), tile1, tile2, &
1890  use_higher_order)
1891  end do
1892  endif
1893  deallocate(tile1_list, tile2_list)
1894  end select
1895  if(grid%on_this_pe) then
1896  grid%area_inv = 0.0;
1897  where (grid%area>0.0) grid%area_inv = 1.0/grid%area
1898  endif
1899  end if
1900  end do
1901 
1902  call mpp_clock_end(id_load_xgrid)
1903 
1904  grid1%area_inv = 0.0;
1905  where (grid1%area>0.0)
1906  grid1%area_inv = 1.0/grid1%area
1907  end where
1908 
1909  xmap%your1my2(xmap%me-xmap%root_pe) = .false. ! this is not necessarily true but keeps
1910  xmap%your2my1(xmap%me-xmap%root_pe) = .false. ! a PE from communicating with itself
1911 
1912  if (make_exchange_reproduce) then
1913  allocate( xmap%send_count_repro(0:xmap%npes-1) )
1914  allocate( xmap%recv_count_repro(0:xmap%npes-1) )
1915  xmap%send_count_repro = 0
1916  xmap%recv_count_repro = 0
1917  do g=2,size(xmap%grids(:))
1918  do p=0,xmap%npes-1
1919  if(xmap%grids(g)%size >0) &
1920  xmap%send_count_repro(p) = xmap%send_count_repro(p) &
1921  +count(xmap%grids(g)%x (:)%pe==p+xmap%root_pe)
1922  if(xmap%grids(g)%size_repro >0) &
1923  xmap%recv_count_repro(p) = xmap%recv_count_repro(p) &
1924  +count(xmap%grids(g)%x_repro(:)%pe==p+xmap%root_pe)
1925  end do
1926  end do
1927  xmap%send_count_repro_tot = sum(xmap%send_count_repro)
1928  xmap%recv_count_repro_tot = sum(xmap%recv_count_repro)
1929  else
1930  xmap%send_count_repro_tot = 0
1931  xmap%recv_count_repro_tot = 0
1932  end if
1933 
1934  if (xgrid_log) then
1935  call mpp_open( unit, 'xgrid.out', action=mpp_overwr, threading=mpp_multi, &
1936  fileset=mpp_multi, nohdrs=.true. )
1937 
1938  write( unit,* )xmap%grids(:)%id, ' GRID: PE ', xmap%me, ' #XCELLS=', &
1939  xmap%grids(2:size(xmap%grids(:)))%size, ' #COMM. PARTNERS=', &
1940  count(xmap%your1my2), '/', count(xmap%your2my1), &
1941  pack((/(p+xmap%root_pe,p=0,xmap%npes-1)/), xmap%your1my2), &
1942  '/', pack((/(p+xmap%root_pe,p=0,xmap%npes-1)/), xmap%your2my1)
1943  call close_file (unit)
1944  endif
1945 
1946  allocate( xmap%x1(1:sum(xmap%grids(2:size(xmap%grids(:)))%size)) )
1947  allocate( xmap%x2(1:sum(xmap%grids(2:size(xmap%grids(:)))%size)) )
1948  allocate( xmap%x1_put(1:sum(xmap%grids(2:size(xmap%grids(:)))%size)) )
1949  allocate( xmap%x2_get(1:sum(xmap%grids(2:size(xmap%grids(:)))%size)) )
1950 
1951  !--- The following will setup indx to be used in regen
1952  allocate(xmap%get1, xmap%put1)
1953  call mpp_clock_begin(id_set_comm)
1954 
1955  call set_comm_get1(xmap)
1956 
1957  call set_comm_put1(xmap)
1958 
1959  if(make_exchange_reproduce) then
1960  allocate(xmap%get1_repro)
1961  call set_comm_get1_repro(xmap)
1962  endif
1963 
1964  call mpp_clock_end(id_set_comm)
1965 
1966  call mpp_clock_begin(id_regen)
1967  call regen(xmap)
1968  call mpp_clock_end(id_regen)
1969 
1970  call mpp_clock_begin(id_conservation_check)
1971 
1972  if(lnd_ug_id ==0) then
1973  xxx = conservation_check(grid1%area*0.0+1.0, grid1%id, xmap)
1974  else
1975  allocate(tmp_2d(grid1%is_me:grid1%ie_me, grid1%js_me:grid1%je_me))
1976  tmp_2d = 1.0
1977  xxx = conservation_check_ug(tmp_2d, grid1%id, xmap)
1978  deallocate(tmp_2d)
1979  endif
1980  write(out_unit,* )"Checked data is array of constant 1"
1981  write(out_unit,* )grid1%id,'(',xmap%grids(:)%id,')=', xxx
1982 
1983  if(lnd_ug_id == 0) then
1984  do g=2,size(xmap%grids(:))
1985  xxx = conservation_check(xmap%grids(g)%frac_area*0.0+1.0, xmap%grids(g)%id, xmap )
1986  write( out_unit,* )xmap%grids(g)%id,'(',xmap%grids(:)%id,')=', xxx
1987  enddo
1988  else
1989  do g=2,size(xmap%grids(:))
1990  grid => xmap%grids(g)
1991  allocate(tmp_3d(grid%is_me:grid%ie_me, grid%js_me:grid%je_me,grid%km))
1992  tmp_3d = 1.0
1993  xxx = conservation_check_ug(tmp_3d, xmap%grids(g)%id, xmap )
1994  write( out_unit,* )xmap%grids(g)%id,'(',xmap%grids(:)%id,')=', xxx
1995  deallocate(tmp_3d)
1996  enddo
1997  endif
1998  ! create an random number 2d array
1999  if(grid1%id == "ATM") then
2000  allocate(check_data(size(grid1%area,1), size(grid1%area,2)))
2001  call random_number(check_data)
2002 
2003  !--- second order along both zonal and meridinal direction
2004  if(lnd_ug_id ==0) then
2005  xxx = conservation_check(check_data, grid1%id, xmap, remap_method = remapping_method )
2006  else
2007  xxx = conservation_check_ug(check_data, grid1%id, xmap, remap_method = remapping_method )
2008  endif
2009  write( out_unit,* ) &
2010  "Checked data is array of random number between 0 and 1 using "//trim(interp_method)
2011  write( out_unit,* )grid1%id,'(',xmap%grids(:)%id,')=', xxx
2012 
2013  deallocate(check_data)
2014  do g=2,size(xmap%grids(:))
2015  allocate(check_data_3d(xmap%grids(g)%is_me:xmap%grids(g)%ie_me, &
2016  xmap%grids(g)%js_me:xmap%grids(g)%je_me, grid1%km))
2017  call random_number(check_data_3d)
2018  if(lnd_ug_id ==0) then
2019  xxx = conservation_check(check_data_3d, xmap%grids(g)%id, xmap, remap_method = remapping_method )
2020  else
2021  xxx = conservation_check_ug(check_data_3d, xmap%grids(g)%id, xmap, remap_method = remapping_method )
2022  endif
2023  write( out_unit,* )xmap%grids(g)%id,'(',xmap%grids(:)%id,')=', xxx
2024  deallocate( check_data_3d)
2025  end do
2026  endif
2027  call mpp_clock_end(id_conservation_check)
2028 
2029  call mpp_clock_end(id_setup_xmap)
2030 
2031 end subroutine setup_xmap
2032 ! </SUBROUTINE>
2033 
2034 !----------------------------------------------------------------------------
2035 ! currently we are assuming there is only one nest region
2036 function get_nest_contact(mosaic_file, tile_nest_out, tile_parent_out, is_nest_out, &
2037  ie_nest_out, js_nest_out, je_nest_out, is_parent_out, &
2038  ie_parent_out, js_parent_out, je_parent_out)
2039 character(len=*), intent(in) :: mosaic_file
2040 integer, intent(out) :: tile_nest_out, tile_parent_out
2041 integer, intent(out) :: is_nest_out, ie_nest_out
2042 integer, intent(out) :: js_nest_out, je_nest_out
2043 integer, intent(out) :: is_parent_out, ie_parent_out
2044 integer, intent(out) :: js_parent_out, je_parent_out
2045 integer :: get_nest_contact
2046 !--- local variables
2047 integer :: ntiles, ncontacts, n, t1, t2
2048 integer :: nx1_contact, ny1_contact
2049 integer :: nx2_contact, ny2_contact
2050 integer, allocatable, dimension(:) :: nx, ny
2051 integer, allocatable, dimension(:) :: tile1, tile2
2052 integer, allocatable, dimension(:) :: istart1, iend1, jstart1, jend1
2053 integer, allocatable, dimension(:) :: istart2, iend2, jstart2, jend2
2054 
2055  tile_nest_out = 0; tile_parent_out = 0
2056  is_nest_out = 0; ie_nest_out = 0
2057  js_nest_out = 0; je_nest_out = 0
2058  is_parent_out = 0; ie_parent_out = 0
2059  js_parent_out = 0; je_parent_out = 0
2060  get_nest_contact = 0
2061 
2062  ! first read the contact information
2063  ntiles = get_mosaic_ntiles(mosaic_file)
2064  if( ntiles == 1 ) return
2065  allocate(nx(ntiles), ny(ntiles))
2066  call get_mosaic_grid_sizes(mosaic_file, nx, ny)
2067 
2068  ncontacts = get_mosaic_ncontacts(mosaic_file)
2069  if(ncontacts == 0) return
2070  allocate(tile1(ncontacts), tile2(ncontacts))
2071  allocate(istart1(ncontacts), iend1(ncontacts))
2072  allocate(jstart1(ncontacts), jend1(ncontacts))
2073  allocate(istart2(ncontacts), iend2(ncontacts))
2074  allocate(jstart2(ncontacts), jend2(ncontacts))
2075 
2076  call get_mosaic_contact( mosaic_file, tile1, tile2, istart1, iend1, jstart1, jend1, &
2077  istart2, iend2, jstart2, jend2)
2078 
2079  do n = 1, ncontacts
2080  if( tile1(n) == tile2(n) ) cycle ! same tile could not be nested
2081 
2082  nx1_contact = iend1(n)-istart1(n)+1
2083  ny1_contact = jend1(n)-jstart1(n)+1
2084  nx2_contact = iend2(n)-istart2(n)+1
2085  ny2_contact = jend2(n)-jstart2(n)+1
2086  t1 = tile1(n);
2087  t2 = tile2(n);
2088  ! For nesting, the contact index of one tile must match its global domain
2089  if( (nx(t1) .NE. nx1_contact .OR. ny(t1) .NE. ny1_contact ) .AND. &
2090  (nx(t2) .NE. nx2_contact .OR. ny(t2) .NE. ny2_contact ) ) cycle
2091  if(nx1_contact == nx2_contact .AND. ny1_contact == ny2_contact) then
2092  call error_mesg('xgrid_mod', 'There is no refinement for the overlapping region', fatal)
2093  endif
2094 
2096  if(get_nest_contact>1) then
2097  call error_mesg('xgrid_mod', 'only support one nest region, contact developer' ,fatal)
2098  endif
2099  if(nx2_contact*ny2_contact > nx1_contact*ny1_contact) then
2100  is_nest_out = istart2(n);
2101  ie_nest_out = iend2(n);
2102  js_nest_out = jstart2(n);
2103  je_nest_out = jend2(n);
2104  tile_nest_out = tile2(n);
2105  is_parent_out = istart1(n);
2106  ie_parent_out = iend1(n);
2107  js_parent_out = jstart1(n);
2108  je_parent_out = jend1(n);
2109  tile_parent_out = tile1(n);
2110  else
2111  is_nest_out = istart1(n);
2112  ie_nest_out = iend1(n);
2113  js_nest_out = jstart1(n);
2114  je_nest_out = jend1(n);
2115  tile_nest_out = tile1(n);
2116  is_parent_out = istart2(n);
2117  ie_parent_out = iend2(n);
2118  js_parent_out = jstart2(n);
2119  je_parent_out = jend2(n);
2120  tile_parent_out = tile2(n);
2121  endif
2122  enddo
2123 
2124  deallocate(nx, ny, tile1, tile2)
2125  deallocate(istart1, iend1, jstart1, jend1)
2126  deallocate(istart2, iend2, jstart2, jend2)
2127 
2128 
2129  return
2130 
2131 end function get_nest_contact
2132 
2133 !#######################################################################
2134 subroutine set_comm_get1_repro(xmap)
2135  type(xmap_type), intent(inout) :: xmap
2136  integer, dimension(xmap%npes) :: pe_ind, cnt
2137  integer, dimension(0:xmap%npes-1) :: send_ind, recv_ind, pl
2138  integer :: npes, nsend, nrecv, mypos
2139  integer :: m, p, pos, n, g, l, im, i, j
2140  type(comm_type), pointer, save :: comm => null()
2141 
2142  comm => xmap%get1_repro
2143  npes = xmap%npes
2144 
2145  nrecv = 0
2146  mypos = mpp_pe() - mpp_root_pe()
2147  do m=0,npes-1
2148  p = mod(mypos+npes-m, npes)
2149  if( xmap%recv_count_repro(p) > 0 ) then
2150  nrecv = nrecv + 1
2151  pe_ind(nrecv) = p
2152  endif
2153  enddo
2154 
2155  comm%nrecv = nrecv
2156  if( nrecv > 0 ) then
2157  allocate(comm%recv(nrecv))
2158  pos = 0
2159  do n = 1, nrecv
2160  p = pe_ind(n)
2161  comm%recv(n)%count = xmap%recv_count_repro(p)
2162  comm%recv(n)%pe = p + xmap%root_pe
2163  comm%recv(n)%buffer_pos = pos
2164  pos = pos + comm%recv(n)%count
2165  enddo
2166  endif
2167 
2168 
2169  ! send information
2170  nsend = 0
2171  mypos = mpp_pe() - mpp_root_pe()
2172  do m=0,xmap%npes-1
2173  p = mod(mypos+m, npes)
2174  if( xmap%send_count_repro(p) > 0 ) then
2175  nsend = nsend + 1
2176  pe_ind(nsend) = p
2177  send_ind(p) = nsend
2178  endif
2179  enddo
2180 
2181  comm%nsend = nsend
2182  if( nsend > 0 ) then
2183  allocate(comm%send(nsend))
2184  pos = 0
2185  cnt(:) = 0
2186  do n = 1, nsend
2187  p = pe_ind(n)
2188  comm%send(n)%count = xmap%send_count_repro(p)
2189  comm%send(n)%pe = p + xmap%root_pe
2190  comm%send(n)%buffer_pos = pos
2191  pos = pos + comm%send(n)%count
2192  allocate(comm%send(n)%i(comm%send(n)%count))
2193  allocate(comm%send(n)%j(comm%send(n)%count))
2194  allocate(comm%send(n)%g(comm%send(n)%count))
2195  allocate(comm%send(n)%xLoc(comm%send(n)%count))
2196  enddo
2197 
2198  do g=2,size(xmap%grids(:))
2199  im = xmap%grids(g)%im
2200  do l=1,xmap%grids(g)%size ! index into this side 2 grid's patterns
2201  p = xmap%grids(g)%x(l)%pe-xmap%root_pe
2202  n = send_ind(p)
2203  cnt(n) = cnt(n) + 1
2204  pos = cnt(n)
2205  i = xmap%grids(g)%x(l)%i2
2206  j = xmap%grids(g)%x(l)%j2
2207  if(xmap%grids(g)%is_ug) then
2208  comm%send(n)%i(pos) = xmap%grids(g)%l_index((j-1)*im+i)
2209  comm%send(n)%j(pos) = 1
2210  else
2211  comm%send(n)%i(pos) = xmap%grids(g)%x(l)%i2
2212  comm%send(n)%j(pos) = xmap%grids(g)%x(l)%j2
2213  endif
2214  comm%send(n)%g(pos) = g
2215  enddo
2216  enddo
2217  !--- make sure the count is correct
2218  do n = 1, nsend
2219  if( comm%send(n)%count .NE. cnt(n) ) call error_mesg('xgrid_mod', &
2220  .NE.'comm%send(n)%count cnt(n)', fatal)
2221  enddo
2222  endif
2223 
2224  !--- set up the recv_pos for unpack the data.
2225  pl(:) = 1
2226  do g=2,size(xmap%grids(:))
2227  do l=1,xmap%grids(g)%size_repro ! index into side1 grid's patterns
2228  p = xmap%grids(g)%x_repro(l)%pe-xmap%root_pe
2229  xmap%grids(g)%x_repro(l)%recv_pos = pl(p)
2230  pl(p) = pl(p) + 1
2231  end do
2232  end do
2233 
2234 
2235 
2236 end subroutine set_comm_get1_repro
2237 
2238 !#######################################################################
2239 subroutine set_comm_get1(xmap)
2240  type(xmap_type), intent(inout) :: xmap
2241  type(grid_type), pointer, save :: grid1 =>null()
2242  integer, allocatable :: send_size(:)
2243  integer, allocatable :: recv_size(:)
2244  integer :: max_size, g, npes, l, ll, nset, m
2245  integer :: i1, j1, tile1, p, n, pos, buffer_pos, mypos
2246  integer :: nsend, nrecv, rbuf_size, sbuf_size, msgsize
2247  logical :: found
2248  real, allocatable :: recv_buf(:), send_buf(:)
2249  real, allocatable :: diarray(:), djarray(:)
2250  integer, allocatable :: iarray(:), jarray(:), tarray(:)
2251  integer, allocatable :: pos_x(:), pelist(:), size_pe(:), pe_side1(:)
2252  integer :: recv_buffer_pos(0:xmap%npes)
2253  integer :: send_buffer_pos(0:xmap%npes)
2254  type(comm_type), pointer, save :: comm => null()
2255  integer :: i, j
2256 
2257  max_size = 0
2258  do g=2,size(xmap%grids(:))
2259  max_size = max_size + xmap%grids(g)%size
2260  enddo
2261  comm => xmap%get1
2262  grid1 => xmap%grids(1)
2263  comm%nsend = 0
2264  comm%nrecv = 0
2265  npes = xmap%npes
2266 
2267  allocate(pelist(0:npes-1))
2268  call mpp_get_current_pelist(pelist)
2269  allocate(send_size(0:npes-1))
2270  allocate(recv_size(0:npes-1))
2271  allocate(size_pe(0:npes-1))
2272  allocate(pos_x(0:npes-1))
2273  size_pe = 0
2274  send_size = 0
2275  recv_size = 0
2276 
2277  if(max_size > 0) then
2278  allocate(pe_side1(max_size))
2279  allocate(xmap%ind_get1(max_size))
2280 
2281  !--- find the recv_indx
2282  ll = 0
2283  do g=2,size(xmap%grids(:))
2284  do l=1,xmap%grids(g)%size
2285  i1 = xmap%grids(g)%x(l)%i1
2286  j1 = xmap%grids(g)%x(l)%j1
2287  tile1 = xmap%grids(g)%x(l)%tile
2288  do p=0,npes-1
2289  if(grid1%tile(p) == tile1) then
2290  if(in_box_nbr(i1, j1, grid1, p)) then
2291  size_pe(p) = size_pe(p) + 1
2292  exit
2293  endif
2294  endif
2295  enddo
2296  if( p == npes ) then
2297  call error_mesg('xgrid_mod', 'tile is not in grid1%tile(:)', fatal)
2298  endif
2299  ll = ll + 1
2300  pe_side1(ll) = p
2301  enddo
2302  enddo
2303 
2304  pos_x = 0
2305  do p = 1, npes-1
2306  pos_x(p) = pos_x(p-1) + size_pe(p-1)
2307  enddo
2308 
2309  !---find the send size for get_1_from_xgrid
2310  allocate(iarray(max_size))
2311  allocate(jarray(max_size))
2312  allocate(tarray(max_size))
2313  if(monotonic_exchange) then
2314  allocate(diarray(max_size))
2315  allocate(djarray(max_size))
2316  endif
2317 
2318  ll = 0
2319 
2320  do g=2,size(xmap%grids(:))
2321  do l=1,xmap%grids(g)%size
2322  i1 = xmap%grids(g)%x(l)%i1
2323  j1 = xmap%grids(g)%x(l)%j1
2324  tile1 = xmap%grids(g)%x(l)%tile
2325  ll = ll + 1
2326  p = pe_side1(ll)
2327 
2328  found = .false.
2329  if(send_size(p) > 0) then
2330  if( i1 == iarray(pos_x(p)+send_size(p)) .AND. j1 == jarray(pos_x(p)+send_size(p)) &
2331  .AND. tile1 == tarray(pos_x(p)+send_size(p))) then
2332  found = .true.
2333  n = send_size(p)
2334  else
2335  !---may need to replace with a fast search algorithm
2336  do n = 1, send_size(p)
2337  if(i1 == iarray(pos_x(p)+n) .AND. j1 == jarray(pos_x(p)+n) .AND. tile1 == tarray(pos_x(p)+n)) then
2338  found = .true.
2339  exit
2340  endif
2341  enddo
2342  endif
2343  endif
2344  if( (.NOT. found) .OR. monotonic_exchange ) then
2345  send_size(p) = send_size(p)+1
2346  pos = pos_x(p)+send_size(p)
2347  iarray(pos) = i1
2348  jarray(pos) = j1
2349  tarray(pos) = tile1
2350  if(monotonic_exchange) then
2351  diarray(pos) = xmap%grids(g)%x(l)%di
2352  djarray(pos) = xmap%grids(g)%x(l)%dj
2353  endif
2354  n = send_size(p)
2355  endif
2356  xmap%ind_get1(ll) = n
2357  enddo
2358  enddo
2359 
2360  pos_x = 0
2361  do p = 1, npes-1
2362  pos_x(p) = pos_x(p-1) + send_size(p-1)
2363  enddo
2364 
2365  ll = 0
2366  do g=2,size(xmap%grids(:))
2367  do l=1,xmap%grids(g)%size
2368  ll = ll + 1
2369  p = pe_side1(ll)
2370  xmap%ind_get1(ll) = pos_x(p) + xmap%ind_get1(ll)
2371  enddo
2372  enddo
2373  endif
2374 
2375  mypos = mpp_pe()-mpp_root_pe()
2376 
2377  ! send/recv for get_1_from_xgrid_recv
2378  recv_size(:) = xmap%your2my1_size(:)
2379  nsend = count( send_size> 0)
2380  comm%nsend = nsend
2381  if(nsend>0) then
2382  allocate(comm%send(nsend))
2383  comm%send(:)%count = 0
2384  endif
2385 
2386  pos = 0
2387  do p = 0, npes-1
2388  send_buffer_pos(p) = pos
2389  pos = pos + send_size(p)
2390  enddo
2391 
2392  pos = 0
2393  comm%sendsize = 0
2394  do n = 0, npes-1
2395  p = mod(mypos+n, npes)
2396  if(send_size(p)>0) then
2397  pos = pos + 1
2398  allocate(comm%send(pos)%i(send_size(p)))
2399  comm%send(pos)%buffer_pos = send_buffer_pos(p)
2400  comm%send(pos)%count = send_size(p)
2401  comm%send(pos)%pe = pelist(p)
2402  comm%sendsize = comm%sendsize + send_size(p)
2403  endif
2404  enddo
2405 
2406  nset = 3
2407  if(monotonic_exchange) nset = 5
2408  rbuf_size = sum(recv_size)*nset
2409  sbuf_size = sum(send_size)*nset
2410  if(rbuf_size>0) allocate(recv_buf(rbuf_size))
2411  if(sbuf_size>0) allocate(send_buf(sbuf_size))
2412 
2413  pos = 0
2414  do n = 0, npes-1
2415  p = mod(mypos+npes-n, npes)
2416  if(recv_size(p) ==0) cycle
2417  msgsize = recv_size(p)*nset
2418  call mpp_recv(recv_buf(pos+1), glen=msgsize, from_pe=pelist(p), block=.false., tag=comm_tag_4)
2419  pos = pos + msgsize
2420  enddo
2421 
2422  pos_x = 0
2423  do p = 1, npes-1
2424  pos_x(p) = pos_x(p-1) + size_pe(p-1)
2425  enddo
2426  ll = 0
2427  pos = 0
2428  do n = 0, npes-1
2429  p = mod(mypos+n, npes)
2430  do l = 1, send_size(p)
2431  send_buf(pos+1) = iarray(pos_x(p)+l)
2432  send_buf(pos+2) = jarray(pos_x(p)+l)
2433  send_buf(pos+3) = tarray(pos_x(p)+l)
2434  if(monotonic_exchange) then
2435  send_buf(pos+4) = diarray(pos_x(p)+l)
2436  send_buf(pos+5) = djarray(pos_x(p)+l)
2437  endif
2438  pos = pos + nset
2439  enddo
2440  enddo
2441 
2442  pos = 0
2443  do n = 0, npes-1
2444  p = mod(mypos+n, npes)
2445  if(send_size(p) ==0) cycle
2446  msgsize = send_size(p)*nset
2447  call mpp_send(send_buf(pos+1), plen=msgsize, to_pe=pelist(p), tag=comm_tag_4 )
2448  pos = pos + msgsize
2449  enddo
2450 
2451  call mpp_sync_self(check=event_recv)
2452  nrecv = count(recv_size>0)
2453  comm%nrecv = nrecv
2454  comm%recvsize = 0
2455 
2456  if(nrecv >0) then
2457  allocate(comm%recv(nrecv))
2458  comm%recv(:)%count = 0
2459  !--- set up the buffer pos for each receiving
2460  buffer_pos = 0
2461  do p = 0, npes-1
2462  recv_buffer_pos(p) = buffer_pos
2463  buffer_pos = buffer_pos + recv_size(p)
2464  enddo
2465  pos = 0
2466  buffer_pos = 0
2467  do m=0,npes-1
2468  p = mod(mypos+npes-m, npes)
2469  if(recv_size(p)>0) then
2470  pos = pos + 1
2471  allocate(comm%recv(pos)%i(recv_size(p)))
2472  allocate(comm%recv(pos)%j(recv_size(p)))
2473  allocate(comm%recv(pos)%tile(recv_size(p)))
2474  comm%recv(pos)%buffer_pos = recv_buffer_pos(p)
2475  comm%recv(pos)%pe = pelist(p)
2476  comm%recv(pos)%count = recv_size(p)
2477  comm%recvsize = comm%recvsize + recv_size(p)
2478  if(monotonic_exchange) then
2479  allocate(comm%recv(pos)%di(recv_size(p)))
2480  allocate(comm%recv(pos)%dj(recv_size(p)))
2481  endif
2482  if(grid1%is_ug) then
2483  do n = 1, recv_size(p)
2484  i = recv_buf(buffer_pos+1)
2485  j = recv_buf(buffer_pos+2)
2486  comm%recv(pos)%i(n) = grid1%l_index((j-1)*grid1%im+i)
2487  comm%recv(pos)%j(n) = 1
2488  comm%recv(pos)%tile(n) = recv_buf(buffer_pos+3)
2489  if(monotonic_exchange) then
2490  comm%recv(pos)%di(n) = recv_buf(buffer_pos+4)
2491  comm%recv(pos)%dj(n) = recv_buf(buffer_pos+5)
2492  endif
2493  buffer_pos = buffer_pos + nset
2494  enddo
2495  else
2496  do n = 1, recv_size(p)
2497  comm%recv(pos)%i(n) = recv_buf(buffer_pos+1) - grid1%is_me + 1
2498  comm%recv(pos)%j(n) = recv_buf(buffer_pos+2) - grid1%js_me + 1
2499  comm%recv(pos)%tile(n) = recv_buf(buffer_pos+3)
2500  if(monotonic_exchange) then
2501  comm%recv(pos)%di(n) = recv_buf(buffer_pos+4)
2502  comm%recv(pos)%dj(n) = recv_buf(buffer_pos+5)
2503  endif
2504  buffer_pos = buffer_pos + nset
2505  enddo
2506  endif
2507  endif
2508  enddo
2509  allocate(comm%unpack_ind(nrecv))
2510  pos = 0
2511  do p = 0, npes-1
2512  if(recv_size(p)>0) then
2513  pos = pos + 1
2514  do m = 1, nrecv
2515  if(comm%recv(m)%pe == pelist(p)) then
2516  comm%unpack_ind(pos) = m
2517  exit
2518  endif
2519  enddo
2520  endif
2521  enddo
2522  endif
2523  call mpp_sync_self()
2524 
2525  if(allocated(send_buf) ) deallocate(send_buf)
2526  if(allocated(recv_buf) ) deallocate(recv_buf)
2527  if(allocated(pelist) ) deallocate(pelist)
2528  if(allocated(pos_x) ) deallocate(pos_x)
2529  if(allocated(pelist) ) deallocate(pelist)
2530  if(allocated(iarray) ) deallocate(iarray)
2531  if(allocated(jarray) ) deallocate(jarray)
2532  if(allocated(tarray) ) deallocate(tarray)
2533  if(allocated(size_pe) ) deallocate(size_pe)
2534 
2535 end subroutine set_comm_get1
2536 
2537 !###############################################################################
2538 subroutine set_comm_put1(xmap)
2539  type(xmap_type), intent(inout) :: xmap
2540  type(grid_type), pointer, save :: grid1 =>null()
2541  integer, allocatable :: send_size(:)
2542  integer, allocatable :: recv_size(:)
2543  integer :: max_size, g, npes, l, ll, m, mypos
2544  integer :: i1, j1, tile1, p, n, pos, buffer_pos
2545  integer :: nsend, nrecv, msgsize, nset, rbuf_size, sbuf_size
2546  logical :: found
2547  real, allocatable :: recv_buf(:), send_buf(:)
2548  real, allocatable :: diarray(:), djarray(:)
2549  integer, allocatable :: iarray(:), jarray(:), tarray(:)
2550  integer, allocatable :: pos_x(:), pelist(:), size_pe(:), pe_put1(:)
2551  integer :: root_pe, recvsize, sendsize
2552  integer :: recv_buffer_pos(0:xmap%npes)
2553  type(comm_type), pointer, save :: comm => null()
2554 
2555 
2556  comm => xmap%put1
2557  if(nnest == 0 .OR. xmap%grids(1)%id .NE. 'ATM' ) then
2558  comm%nsend = xmap%get1%nrecv
2559  comm%nrecv = xmap%get1%nsend
2560  comm%sendsize = xmap%get1%recvsize
2561  comm%recvsize = xmap%get1%sendsize
2562  comm%send => xmap%get1%recv
2563  comm%recv => xmap%get1%send
2564  xmap%ind_put1 => xmap%ind_get1
2565  return
2566  endif
2567 
2568  max_size = 0
2569  do g=2,size(xmap%grids(:))
2570  max_size = max_size + xmap%grids(g)%size
2571  enddo
2572  grid1 => xmap%grids(1)
2573  comm%nsend = 0
2574  comm%nrecv = 0
2575  npes = xmap%npes
2576  allocate(pelist(0:npes-1))
2577  call mpp_get_current_pelist(pelist)
2578  allocate(send_size(0:npes-1))
2579  allocate(recv_size(0:npes-1))
2580  allocate(size_pe(0:npes-1))
2581  allocate(pos_x(0:npes-1))
2582  size_pe = 0
2583  send_size = 0
2584  recv_size = 0
2585 
2586  if(max_size > 0) then
2587  allocate(pe_put1(max_size))
2588  allocate(xmap%ind_put1(max_size))
2589 
2590  !--- find the recv_indx
2591  ll = 0
2592  do g=2,size(xmap%grids(:))
2593  do l=1,xmap%grids(g)%size
2594  i1 = xmap%grids(g)%x(l)%i1
2595  j1 = xmap%grids(g)%x(l)%j1
2596  tile1 = xmap%grids(g)%x(l)%tile
2597  do p=0,npes-1
2598  if(grid1%tile(p) == tile1) then
2599  if(in_box(i1, j1, grid1%is(p), grid1%ie(p), grid1%js(p), grid1%je(p))) then
2600  size_pe(p) = size_pe(p) + 1
2601  exit
2602  endif
2603  endif
2604  enddo
2605  ll = ll + 1
2606  pe_put1(ll) = p
2607  enddo
2608  enddo
2609 
2610  pos_x = 0
2611  do p = 1, npes-1
2612  pos_x(p) = pos_x(p-1) + size_pe(p-1)
2613  enddo
2614 
2615  !---find the send size for get_1_from_xgrid
2616  allocate(iarray(max_size))
2617  allocate(jarray(max_size))
2618  allocate(tarray(max_size))
2619  if(monotonic_exchange) then
2620  allocate(diarray(max_size))
2621  allocate(djarray(max_size))
2622  endif
2623 
2624  ll = 0
2625 
2626  do g=2,size(xmap%grids(:))
2627  do l=1,xmap%grids(g)%size
2628  i1 = xmap%grids(g)%x(l)%i1
2629  j1 = xmap%grids(g)%x(l)%j1
2630  tile1 = xmap%grids(g)%x(l)%tile
2631  ll = ll + 1
2632  p = pe_put1(ll)
2633 
2634  found = .false.
2635  if(send_size(p) > 0) then
2636  if( i1 == iarray(pos_x(p)+send_size(p)) .AND. j1 == jarray(pos_x(p)+send_size(p)) &
2637  .AND. tile1 == tarray(pos_x(p)+send_size(p))) then
2638  found = .true.
2639  n = send_size(p)
2640  else
2641  !---may need to replace with a fast search algorithm
2642  do n = 1, send_size(p)
2643  if(i1 == iarray(pos_x(p)+n) .AND. j1 == jarray(pos_x(p)+n) .AND. tile1 == tarray(pos_x(p)+n)) then
2644  found = .true.
2645  exit
2646  endif
2647  enddo
2648  endif
2649  endif
2650  if( (.NOT. found) .OR. monotonic_exchange ) then
2651  send_size(p) = send_size(p)+1
2652  pos = pos_x(p)+send_size(p)
2653  iarray(pos) = i1
2654  jarray(pos) = j1
2655  tarray(pos) = tile1
2656  if(monotonic_exchange) then
2657  diarray(pos) = xmap%grids(g)%x(l)%di
2658  djarray(pos) = xmap%grids(g)%x(l)%dj
2659  endif
2660  n = send_size(p)
2661  endif
2662  xmap%ind_put1(ll) = n
2663  enddo
2664  enddo
2665 
2666  pos_x = 0
2667  do p = 1, npes-1
2668  pos_x(p) = pos_x(p-1) + send_size(p-1)
2669  enddo
2670 
2671  ll = 0
2672  do g=2,size(xmap%grids(:))
2673  do l=1,xmap%grids(g)%size
2674  i1 = xmap%grids(g)%x(l)%i1
2675  j1 = xmap%grids(g)%x(l)%j1
2676  tile1 = xmap%grids(g)%x(l)%tile
2677  ll = ll + 1
2678  p = pe_put1(ll)
2679  xmap%ind_put1(ll) = pos_x(p) + xmap%ind_put1(ll)
2680  enddo
2681  enddo
2682  endif
2683 
2684 
2685  mypos = mpp_pe()-mpp_root_pe()
2686 
2687  if (do_alltoall) then
2688  call mpp_alltoall(send_size, 1, recv_size, 1)
2689  else
2690  do n = 0, npes-1
2691  p = mod(mypos+npes-n, npes)
2692  call mpp_recv(recv_size(p), glen=1, from_pe=pelist(p), block=.false., tag=comm_tag_5)
2693  enddo
2694 
2695  !--- send data
2696  do n = 0, npes-1
2697  p = mod(mypos+n, npes)
2698  call mpp_send(send_size(p), plen=1, to_pe=pelist(p), tag=comm_tag_5)
2699  enddo
2700 
2701  call mpp_sync_self(check=event_recv)
2702  call mpp_sync_self()
2703  endif
2704  !--- recv for put_1_to_xgrid
2705  nrecv = count( send_size> 0)
2706  comm%nrecv = nrecv
2707  if(nrecv>0) then
2708  allocate(comm%recv(nrecv))
2709  comm%recv(:)%count = 0
2710  endif
2711  pos = 0
2712  comm%recvsize = 0
2713  do p = 0, npes-1
2714  recv_buffer_pos(p) = pos
2715  pos = pos + send_size(p)
2716  enddo
2717 
2718  pos = 0
2719  do n = 0, npes-1
2720  p = mod(mypos+npes-n, npes)
2721  if(send_size(p)>0) then
2722  pos = pos + 1
2723  allocate(comm%recv(pos)%i(send_size(p)))
2724  comm%recv(pos)%buffer_pos = recv_buffer_pos(p)
2725  comm%recv(pos)%count = send_size(p)
2726  comm%recv(pos)%pe = pelist(p)
2727  comm%recvsize = comm%recvsize + send_size(p)
2728  endif
2729  enddo
2730 
2731  nset = 3
2732  if(monotonic_exchange) nset = 5
2733  rbuf_size = sum(recv_size)*nset
2734  sbuf_size = sum(send_size)*nset
2735  if(rbuf_size>0) allocate(recv_buf(rbuf_size))
2736  if(sbuf_size>0) allocate(send_buf(sbuf_size))
2737 
2738  pos = 0
2739  do n = 0, npes-1
2740  p = mod(mypos+npes-n, npes)
2741  if(recv_size(p) ==0) cycle
2742  msgsize = recv_size(p)*nset
2743  call mpp_recv(recv_buf(pos+1), glen=msgsize, from_pe=pelist(p), block=.false., tag=comm_tag_6)
2744  pos = pos + msgsize
2745  enddo
2746 
2747  pos_x = 0
2748  do p = 1, npes-1
2749  pos_x(p) = pos_x(p-1) + size_pe(p-1)
2750  enddo
2751  ll = 0
2752  pos = 0
2753  do n = 0, npes-1
2754  p = mod(mypos+n, npes)
2755  do l = 1, send_size(p)
2756  send_buf(pos+1) = iarray(pos_x(p)+l)
2757  send_buf(pos+2) = jarray(pos_x(p)+l)
2758  send_buf(pos+3) = tarray(pos_x(p)+l)
2759  if(monotonic_exchange) then
2760  send_buf(pos+4) = diarray(pos_x(p)+l)
2761  send_buf(pos+5) = djarray(pos_x(p)+l)
2762  endif
2763  pos = pos + nset
2764  enddo
2765  enddo
2766 
2767  pos = 0
2768  do n = 0, npes-1
2769  p = mod(mypos+n, npes)
2770  if(send_size(p) ==0) cycle
2771  msgsize = send_size(p)*nset
2772  call mpp_send(send_buf(pos+1), plen=msgsize, to_pe=pelist(p), tag=comm_tag_6 )
2773  pos = pos + msgsize
2774  enddo
2775 
2776  call mpp_sync_self(check=event_recv)
2777  nsend = count(recv_size>0)
2778  comm%nsend = nsend
2779  comm%sendsize = 0
2780 
2781  if(nsend >0) then
2782  allocate(comm%send(nsend))
2783  comm%send(:)%count = 0
2784  pos = 0
2785  buffer_pos = 0
2786  do m=0,npes-1
2787  p = mod(mypos+npes-m, npes)
2788  if(recv_size(p)>0) then
2789  pos = pos + 1
2790  allocate(comm%send(pos)%i(recv_size(p)))
2791  allocate(comm%send(pos)%j(recv_size(p)))
2792  allocate(comm%send(pos)%tile(recv_size(p)))
2793  comm%send(pos)%pe = pelist(p)
2794  comm%send(pos)%count = recv_size(p)
2795  comm%sendsize = comm%sendsize + recv_size(p)
2796  if(monotonic_exchange) then
2797  allocate(comm%send(pos)%di(recv_size(p)))
2798  allocate(comm%send(pos)%dj(recv_size(p)))
2799  endif
2800  do n = 1, recv_size(p)
2801  comm%send(pos)%i(n) = recv_buf(buffer_pos+1) - grid1%is_me + 1
2802  comm%send(pos)%j(n) = recv_buf(buffer_pos+2) - grid1%js_me + 1
2803  comm%send(pos)%tile(n) = recv_buf(buffer_pos+3)
2804  if(monotonic_exchange) then
2805  comm%send(pos)%di(n) = recv_buf(buffer_pos+4)
2806  comm%send(pos)%dj(n) = recv_buf(buffer_pos+5)
2807  endif
2808  buffer_pos = buffer_pos + nset
2809  enddo
2810  endif
2811  enddo
2812  endif
2813 
2814  call mpp_sync_self()
2815  if(allocated(send_buf) ) deallocate(send_buf)
2816  if(allocated(recv_buf) ) deallocate(recv_buf)
2817  if(allocated(pelist) ) deallocate(pelist)
2818  if(allocated(pos_x) ) deallocate(pos_x)
2819  if(allocated(pelist) ) deallocate(pelist)
2820  if(allocated(iarray) ) deallocate(iarray)
2821  if(allocated(jarray) ) deallocate(jarray)
2822  if(allocated(tarray) ) deallocate(tarray)
2823  if(allocated(size_pe) ) deallocate(size_pe)
2824 
2825 end subroutine set_comm_put1
2826 
2827 
2828 !###############################################################################
2829 subroutine regen(xmap)
2830 type(xmap_type), intent(inout) :: xmap
2831 
2832  integer :: g, l, k, max_size
2833  integer :: i1, j1, i2, j2, p
2834  integer :: tile1
2835  integer :: ll, lll
2836  logical :: overlap_with_nest
2837  integer :: cnt(xmap%get1%nsend)
2838  integer :: i,j,n,xloc,pos,nsend,m,npes, mypos
2839  integer :: send_ind(0:xmap%npes-1)
2840 
2841  max_size = 0
2842 
2843  do g=2,size(xmap%grids(:))
2844  max_size = max_size + xmap%grids(g)%size * xmap%grids(g)%km
2845  end do
2846 
2847  if (max_size>size(xmap%x1(:))) then
2848  deallocate(xmap%x1)
2849  deallocate(xmap%x2)
2850  allocate( xmap%x1(1:max_size) )
2851  allocate( xmap%x2(1:max_size) )
2852  endif
2853 
2854 
2855  do g=2,size(xmap%grids(:))
2856  xmap%grids(g)%first = 1
2857  xmap%grids(g)%last = 0
2858  end do
2859 
2860  xmap%size = 0
2861  ll = 0
2862  do g=2,size(xmap%grids(:))
2863  xmap%grids(g)%first = xmap%size + 1;
2864 
2865  do l=1,xmap%grids(g)%size
2866  i1 = xmap%grids(g)%x(l)%i1
2867  j1 = xmap%grids(g)%x(l)%j1
2868  i2 = xmap%grids(g)%x(l)%i2
2869  j2 = xmap%grids(g)%x(l)%j2
2870  tile1 = xmap%grids(g)%x(l)%tile
2871  ll = ll + 1
2872  if(xmap%grids(g)%is_ug) then
2873  do k=1,xmap%grids(g)%km
2874  lll = xmap%grids(g)%l_index((j2-1)*xmap%grids(g)%im+i2)
2875  if (xmap%grids(g)%frac_area(lll,1,k)/=0.0) then
2876  xmap%size = xmap%size+1
2877  xmap%x1(xmap%size)%pos = xmap%ind_get1(ll)
2878  xmap%x1(xmap%size)%i = xmap%grids(g)%x(l)%i1
2879  xmap%x1(xmap%size)%j = xmap%grids(g)%x(l)%j1
2880  xmap%x1(xmap%size)%tile = xmap%grids(g)%x(l)%tile
2881  xmap%x1(xmap%size)%area = xmap%grids(g)%x(l)%area &
2882  *xmap%grids(g)%frac_area(lll,1,k)
2883  xmap%x1(xmap%size)%di = xmap%grids(g)%x(l)%di
2884  xmap%x1(xmap%size)%dj = xmap%grids(g)%x(l)%dj
2885  xmap%x2(xmap%size)%i = xmap%grids(g)%x(l)%i2
2886  xmap%x2(xmap%size)%j = xmap%grids(g)%x(l)%j2
2887  xmap%x2(xmap%size)%l = lll
2888  xmap%x2(xmap%size)%k = k
2889  xmap%x2(xmap%size)%area = xmap%grids(g)%x(l)%area * xmap%grids(g)%x(l)%scale
2890  endif
2891  enddo
2892  else
2893  do k=1,xmap%grids(g)%km
2894  if (xmap%grids(g)%frac_area(i2,j2,k)/=0.0) then
2895  xmap%size = xmap%size+1
2896  xmap%x1(xmap%size)%pos = xmap%ind_get1(ll)
2897  xmap%x1(xmap%size)%i = xmap%grids(g)%x(l)%i1
2898  xmap%x1(xmap%size)%j = xmap%grids(g)%x(l)%j1
2899  xmap%x1(xmap%size)%tile = xmap%grids(g)%x(l)%tile
2900  xmap%x1(xmap%size)%area = xmap%grids(g)%x(l)%area &
2901  *xmap%grids(g)%frac_area(i2,j2,k)
2902  xmap%x1(xmap%size)%di = xmap%grids(g)%x(l)%di
2903  xmap%x1(xmap%size)%dj = xmap%grids(g)%x(l)%dj
2904  xmap%x2(xmap%size)%i = xmap%grids(g)%x(l)%i2
2905  xmap%x2(xmap%size)%j = xmap%grids(g)%x(l)%j2
2906  xmap%x2(xmap%size)%k = k
2907  xmap%x2(xmap%size)%area = xmap%grids(g)%x(l)%area * xmap%grids(g)%x(l)%scale
2908  end if
2909  enddo
2910  end if
2911  end do
2912  xmap%grids(g)%last = xmap%size
2913  end do
2914 
2915 
2916  if (max_size>size(xmap%x1_put(:))) then
2917  deallocate(xmap%x1_put)
2918  allocate( xmap%x1_put(1:max_size) )
2919  endif
2920  if (max_size>size(xmap%x2_get(:))) then
2921  deallocate(xmap%x2_get)
2922  allocate( xmap%x2_get(1:max_size) )
2923  endif
2924 
2925  do g=2,size(xmap%grids(:))
2926  xmap%grids(g)%first_get = 1
2927  xmap%grids(g)%last_get = 0
2928  end do
2929 
2930  xmap%size_put1 = 0
2931  xmap%size_get2 = 0
2932  ll = 0
2933  do g=2,size(xmap%grids(:))
2934  xmap%grids(g)%first_get = xmap%size_get2 + 1;
2935 
2936  do l=1,xmap%grids(g)%size
2937  i1 = xmap%grids(g)%x(l)%i1
2938  j1 = xmap%grids(g)%x(l)%j1
2939  i2 = xmap%grids(g)%x(l)%i2
2940  j2 = xmap%grids(g)%x(l)%j2
2941  tile1 = xmap%grids(g)%x(l)%tile
2942  ll = ll + 1
2943  overlap_with_nest = .false.
2944  if( xmap%grids(1)%id == "ATM" .AND. tile1 == tile_parent .AND. &
2945  in_box(i1, j1, is_parent, ie_parent, js_parent, je_parent) ) overlap_with_nest = .true.
2946  if(xmap%grids(g)%is_ug) then
2947  do k=1,xmap%grids(g)%km
2948  lll = xmap%grids(g)%l_index((j2-1)*xmap%grids(g)%im+i2)
2949  if (xmap%grids(g)%frac_area(lll,1,k)/=0.0) then
2950  xmap%size_put1 = xmap%size_put1+1
2951  xmap%x1_put(xmap%size_put1)%pos = xmap%ind_put1(ll)
2952  xmap%x1_put(xmap%size_put1)%i = xmap%grids(g)%x(l)%i1
2953  xmap%x1_put(xmap%size_put1)%j = xmap%grids(g)%x(l)%j1
2954  xmap%x1_put(xmap%size_put1)%tile = xmap%grids(g)%x(l)%tile
2955  xmap%x1_put(xmap%size_put1)%area = xmap%grids(g)%x(l)%area &
2956  *xmap%grids(g)%frac_area(lll,1,k)
2957  xmap%x1_put(xmap%size_put1)%di = xmap%grids(g)%x(l)%di
2958  xmap%x1_put(xmap%size_put1)%dj = xmap%grids(g)%x(l)%dj
2959  if( .not. overlap_with_nest) then
2960  xmap%size_get2 = xmap%size_get2+1
2961  xmap%x2_get(xmap%size_get2)%i = xmap%grids(g)%x(l)%i2
2962  xmap%x2_get(xmap%size_get2)%j = xmap%grids(g)%x(l)%j2
2963  xmap%x2_get(xmap%size_get2)%l = lll
2964  xmap%x2_get(xmap%size_get2)%k = k
2965  xmap%x2_get(xmap%size_get2)%area = xmap%grids(g)%x(l)%area * xmap%grids(g)%x(l)%scale
2966  xmap%x2_get(xmap%size_get2)%pos = xmap%size_put1
2967  endif
2968  end if
2969  end do
2970  else
2971  do k=1,xmap%grids(g)%km
2972  if (xmap%grids(g)%frac_area(i2,j2,k)/=0.0) then
2973  xmap%size_put1 = xmap%size_put1+1
2974  xmap%x1_put(xmap%size_put1)%pos = xmap%ind_put1(ll)
2975  xmap%x1_put(xmap%size_put1)%i = xmap%grids(g)%x(l)%i1
2976  xmap%x1_put(xmap%size_put1)%j = xmap%grids(g)%x(l)%j1
2977  xmap%x1_put(xmap%size_put1)%tile = xmap%grids(g)%x(l)%tile
2978  xmap%x1_put(xmap%size_put1)%area = xmap%grids(g)%x(l)%area &
2979  *xmap%grids(g)%frac_area(i2,j2,k)
2980  xmap%x1_put(xmap%size_put1)%di = xmap%grids(g)%x(l)%di
2981  xmap%x1_put(xmap%size_put1)%dj = xmap%grids(g)%x(l)%dj
2982  if( .not. overlap_with_nest) then
2983  xmap%size_get2 = xmap%size_get2+1
2984  xmap%x2_get(xmap%size_get2)%i = xmap%grids(g)%x(l)%i2
2985  xmap%x2_get(xmap%size_get2)%j = xmap%grids(g)%x(l)%j2
2986  xmap%x2_get(xmap%size_get2)%k = k
2987  xmap%x2_get(xmap%size_get2)%area = xmap%grids(g)%x(l)%area * xmap%grids(g)%x(l)%scale
2988  xmap%x2_get(xmap%size_get2)%pos = xmap%size_put1
2989  endif
2990  end if
2991  end do
2992  endif
2993  end do
2994  xmap%grids(g)%last_get = xmap%size_get2
2995  end do
2996 
2997  !---set up information for get_1_from_xgrid_repro
2998  if (make_exchange_reproduce) then
2999  if (xmap%get1_repro%nsend > 0) then
3000  xloc = 0
3001  nsend = 0
3002  npes = xmap%npes
3003  mypos = mpp_pe() - mpp_root_pe()
3004  cnt(:) = 0
3005  do m=0,npes-1
3006  p = mod(mypos+m, npes)
3007  if( xmap%send_count_repro(p) > 0 ) then
3008  nsend = nsend + 1
3009  send_ind(p) = nsend
3010  endif
3011  enddo
3012  do g=2,size(xmap%grids(:))
3013  do l=1,xmap%grids(g)%size ! index into this side 2 grid's patterns
3014  p = xmap%grids(g)%x(l)%pe-xmap%root_pe
3015  n = send_ind(p)
3016  cnt(n) = cnt(n) + 1
3017  pos = cnt(n)
3018  xmap%get1_repro%send(n)%xLoc(pos) = xloc
3019  if( xmap%grids(g)%is_ug ) then
3020  i = xmap%grids(g)%x(l)%l2
3021  xloc = xloc + count(xmap%grids(g)%frac_area(i,1,:)/=0.0)
3022  else
3023  i = xmap%grids(g)%x(l)%i2
3024  j = xmap%grids(g)%x(l)%j2
3025  xloc = xloc + count(xmap%grids(g)%frac_area(i,j,:)/=0.0)
3026  endif
3027  enddo
3028  enddo
3029  endif
3030  endif
3031 
3032 end subroutine regen
3033 
3034 !#######################################################################
3035 
3036 ! <SUBROUTINE NAME="set_frac_area">
3037 
3038 ! <OVERVIEW>
3039 ! Changes sub-grid portion areas and/or number.
3040 ! </OVERVIEW>
3041 ! <DESCRIPTION>
3042 ! Changes sub-grid portion areas and/or number.
3043 ! </DESCRIPTION>
3044 ! <TEMPLATE>
3045 ! call set_frac_area(f, grid_id, xmap)
3046 ! </TEMPLATE>
3047 
3048 ! <IN NAME="f" TYPE="real" DIM="(:,:,:)"> </IN>
3049 ! <IN NAME="grid_id" TYPE="character(len=3)" > </IN>
3050 ! <INOUT NAME="xmap" TYPE="xmap_type" > </INOUT>
3051 
3052 subroutine set_frac_area_sg(f, grid_id, xmap)
3053 real, dimension(:,:,:), intent(in ) :: f
3054 character(len=3), intent(in ) :: grid_id
3055 type(xmap_type), intent(inout) :: xmap
3056 
3057  integer :: g
3058  type(grid_type), pointer, save :: grid =>null()
3059 
3060  if (grid_id==xmap%grids(1)%id) call error_mesg ('xgrid_mod', &
3061  'set_frac_area called on side 1 grid', fatal)
3062  do g=2,size(xmap%grids(:))
3063  grid => xmap%grids(g)
3064  if (grid_id==grid%id) then
3065  if (size(f,3)/=size(grid%frac_area,3)) then
3066  deallocate (grid%frac_area)
3067  grid%km = size(f,3);
3068  allocate( grid%frac_area(grid%is_me:grid%ie_me, grid%js_me:grid%je_me, &
3069  grid%km) )
3070  end if
3071  grid%frac_area = f;
3072  call regen(xmap)
3073  return;
3074  end if
3075  end do
3076 
3077  call error_mesg ('xgrid_mod', 'set_frac_area: could not find grid id', fatal)
3078 
3079 end subroutine set_frac_area_sg
3080 ! </SUBROUTINE>
3081 
3082 !#######################################################################
3083 
3084 ! <SUBROUTINE NAME="set_frac_area_ug">
3085 
3086 ! <OVERVIEW>
3087 ! Changes sub-grid portion areas and/or number.
3088 ! </OVERVIEW>
3089 ! <DESCRIPTION>
3090 ! Changes sub-grid portion areas and/or number.
3091 ! </DESCRIPTION>
3092 ! <TEMPLATE>
3093 ! call set_frac_area_ug(f, grid_id, xmap)
3094 ! </TEMPLATE>
3095 
3096 ! <IN NAME="f" TYPE="real" DIM="(:,:,:)"> </IN>
3097 ! <IN NAME="grid_id" TYPE="character(len=3)" > </IN>
3098 ! <INOUT NAME="xmap" TYPE="xmap_type" > </INOUT>
3099 
3100 subroutine set_frac_area_ug(f, grid_id, xmap)
3101 real, dimension(:,:), intent(in ) :: f
3102 character(len=3), intent(in ) :: grid_id
3103 type(xmap_type), intent(inout) :: xmap
3104 
3105  integer :: g
3106  type(grid_type), pointer, save :: grid =>null()
3107 
3108  if (grid_id==xmap%grids(1)%id) call error_mesg ('xgrid_mod', &
3109  'set_frac_area_ug called on side 1 grid', fatal)
3110  if (grid_id .NE. 'LND' ) call error_mesg ('xgrid_mod', &
3111  .NE.'set_frac_area_ug called for grid_id LND', fatal)
3112  do g=2,size(xmap%grids(:))
3113  grid => xmap%grids(g)
3114  if (grid_id==grid%id) then
3115  if (size(f,2)/=size(grid%frac_area,3)) then
3116  deallocate (grid%frac_area)
3117  grid%km = size(f,2);
3118  allocate( grid%frac_area(grid%ls_me:grid%le_me, 1, grid%km) )
3119  end if
3120  grid%frac_area(:,1,:) = f(:,:);
3121  call regen(xmap)
3122  return;
3123  end if
3124  end do
3125 
3126  call error_mesg ('xgrid_mod', 'set_frac_area_ug: could not find grid id', fatal)
3127 
3128 end subroutine set_frac_area_ug
3129 ! </SUBROUTINE>
3130 
3131 
3132 
3133 !#######################################################################
3134 
3135 ! <FUNCTION NAME="xgrid_count">
3136 
3137 ! <OVERVIEW>
3138 ! Returns current size of exchange grid variables.
3139 ! </OVERVIEW>
3140 ! <DESCRIPTION>
3141 ! Returns current size of exchange grid variables.
3142 ! </DESCRIPTION>
3143 ! <TEMPLATE>
3144 ! xgrid_count(xmap)
3145 ! </TEMPLATE>
3146 
3147 ! <IN NAME="xmap" TYPE="xmap_type" > </IN>
3148 ! <OUT NAME="xgrid_count" TYPE="integer" > </OUT>
3149 
3150 integer function xgrid_count(xmap)
3151 type(xmap_type), intent(inout) :: xmap
3152 
3153  xgrid_count = xmap%size
3154 end function xgrid_count
3155 ! </FUNCTION>
3156 
3157 !#######################################################################
3158 
3159 ! <SUBROUTINE NAME="put_side1_to_xgrid" INTERFACE="put_to_xgrid">
3160 ! <IN NAME="d" TYPE="real" DIM="(:,:)" > </IN>
3161 ! <IN NAME="grid_id" TYPE=" character(len=3)" > </IN>
3162 ! <INOUT NAME="x" TYPE="real" DIM="(:)" > </INOUT>
3163 ! <INOUT NAME="xmap" TYPE="xmap_type" > </INOUT>
3164 ! <IN NAME="remap_method" TYPE="integer,optional"></IN>
3165 
3166 subroutine put_side1_to_xgrid(d, grid_id, x, xmap, remap_method, complete)
3167  real, dimension(:,:), intent(in ) :: d
3168  character(len=3), intent(in ) :: grid_id
3169  real, dimension(:), intent(inout) :: x
3170  type(xmap_type), intent(inout) :: xmap
3171  integer, intent(in), optional :: remap_method
3172  logical, intent(in), optional :: complete
3173 
3174  logical :: is_complete, set_mismatch
3175  integer :: g, method
3176  character(len=2) :: text
3177  integer, save :: isize=0
3178  integer, save :: jsize=0
3179  integer, save :: lsize=0
3180  integer, save :: xsize=0
3181  integer, save :: method_saved=0
3182  character(len=3), save :: grid_id_saved=""
3183  integer(LONG_KIND), dimension(MAX_FIELDS), save :: d_addrs=-9999
3184  integer(LONG_KIND), dimension(MAX_FIELDS), save :: x_addrs=-9999
3185 
3186  if (grid_id==xmap%grids(1)%id) then
3187  method = first_order ! default
3188  if(present(remap_method)) method = remap_method
3189  is_complete = .true.
3190  if(present(complete)) is_complete=complete
3191  lsize = lsize + 1
3192  if( lsize > max_fields ) then
3193  write( text,'(i2)' ) max_fields
3194  call error_mesg ('xgrid_mod', 'MAX_FIELDS='//trim(text)//' exceeded for group put_side1_to_xgrid', fatal)
3195  endif
3196  d_addrs(lsize) = loc(d)
3197  x_addrs(lsize) = loc(x)
3198 
3199  if(lsize == 1) then
3200  isize = size(d,1)
3201  jsize = size(d,2)
3202  xsize = size(x(:))
3203  method_saved = method
3204  grid_id_saved = grid_id
3205  else
3206  set_mismatch = .false.
3207  set_mismatch = set_mismatch .OR. (isize /= size(d,1))
3208  set_mismatch = set_mismatch .OR. (jsize /= size(d,2))
3209  set_mismatch = set_mismatch .OR. (xsize /= size(x(:)))
3210  set_mismatch = set_mismatch .OR. (method_saved /= method)
3211  set_mismatch = set_mismatch .OR. (grid_id_saved /= grid_id)
3212  if(set_mismatch)then
3213  write( text,'(i2)' ) lsize
3214  call error_mesg ('xgrid_mod', 'Incompatible field at count '//text//' for group put_side1_to_xgrid', fatal )
3215  endif
3216  endif
3217 
3218  if(is_complete) then
3219  !--- when exchange_monotonic is true and the side 1 ia atm, will always use monotonic second order conservative.
3220  if(monotonic_exchange .AND. grid_id == 'ATM') then
3221  call put_1_to_xgrid_order_2(d_addrs, x_addrs, xmap, isize, jsize, xsize, lsize)
3222  else if(method == first_order) then
3223  call put_1_to_xgrid_order_1(d_addrs, x_addrs, xmap, isize, jsize, xsize, lsize)
3224  else
3225  if(grid_id .NE. 'ATM') call error_mesg ('xgrid_mod', &
3226  "second order put_to_xgrid should only be applied to 'ATM' model, "//&
3227  "contact developer", fatal)
3228  call put_1_to_xgrid_order_2(d_addrs, x_addrs, xmap, isize, jsize, xsize, lsize)
3229  endif
3230 
3231  d_addrs = -9999
3232  x_addrs = -9999
3233  isize = 0
3234  jsize = 0
3235  xsize = 0
3236  lsize = 0
3237  method_saved = 0
3238  grid_id_saved = ""
3239  endif
3240  return
3241  end if
3242 
3243  do g=2,size(xmap%grids(:))
3244  if (grid_id==xmap%grids(g)%id) &
3245  call error_mesg ('xgrid_mod', &
3246  'put_to_xgrid expects a 3D side 2 grid', fatal)
3247  end do
3248 
3249  call error_mesg ('xgrid_mod', 'put_to_xgrid: could not find grid id', fatal)
3250 
3251 end subroutine put_side1_to_xgrid
3252 ! </SUBROUTINE>
3253 
3254 !#######################################################################
3255 
3256 ! <SUBROUTINE NAME="put_side2_to_xgrid" INTERFACE="put_to_xgrid">
3257 ! <IN NAME="d" TYPE="real" DIM="(:,:,:)" > </IN>
3258 ! <IN NAME="grid_id" TYPE=" character(len=3)" > </IN>
3259 ! <INOUT NAME="x" TYPE="real" DIM="(:)" > </INOUT>
3260 ! <INOUT NAME="xmap" TYPE="xmap_type" > </INOUT>
3261 
3262 subroutine put_side2_to_xgrid(d, grid_id, x, xmap)
3263 real, dimension(:,:,:), intent(in ) :: d
3264 character(len=3), intent(in ) :: grid_id
3265 real, dimension(:), intent(inout) :: x
3266 type(xmap_type), intent(inout) :: xmap
3267 
3268  integer :: g
3269 
3270  if (grid_id==xmap%grids(1)%id) &
3271  call error_mesg ('xgrid_mod', &
3272  'put_to_xgrid expects a 2D side 1 grid', fatal)
3273 
3274  do g=2,size(xmap%grids(:))
3275  if (grid_id==xmap%grids(g)%id) then
3276  call put_2_to_xgrid(d, xmap%grids(g), x, xmap)
3277  return;
3278  end if
3279  end do
3280 
3281  call error_mesg ('xgrid_mod', 'put_to_xgrid: could not find grid id', fatal)
3282 
3283 end subroutine put_side2_to_xgrid
3284 ! </SUBROUTINE>
3285 
3286 !#######################################################################
3287 
3288 ! <SUBROUTINE NAME="get_side1_from_xgrid" INTERFACE="get_from_xgrid">
3289 ! <IN NAME="x" TYPE="real" DIM="(:)" > </IN>
3290 ! <IN NAME="grid_id" TYPE=" character(len=3)" > </IN>
3291 ! <OUT NAME="d" TYPE="real" DIM="(:,:)" > </OUT>
3292 ! <INOUT NAME="xmap" TYPE="xmap_type" > </INOUT>
3293 
3294 subroutine get_side1_from_xgrid(d, grid_id, x, xmap, complete)
3295  real, dimension(:,:), intent( out) :: d
3296  character(len=3), intent(in ) :: grid_id
3297  real, dimension(:), intent(in ) :: x
3298  type(xmap_type), intent(inout) :: xmap
3299  logical, intent(in), optional :: complete
3300 
3301  logical :: is_complete, set_mismatch
3302  integer :: g
3303  character(len=2) :: text
3304  integer, save :: isize=0
3305  integer, save :: jsize=0
3306  integer, save :: lsize=0
3307  integer, save :: xsize=0
3308  character(len=3), save :: grid_id_saved=""
3309  integer(LONG_KIND), dimension(MAX_FIELDS), save :: d_addrs=-9999
3310  integer(LONG_KIND), dimension(MAX_FIELDS), save :: x_addrs=-9999
3311 
3312  if (grid_id==xmap%grids(1)%id) then
3313  is_complete = .true.
3314  if(present(complete)) is_complete=complete
3315  lsize = lsize + 1
3316  if( lsize > max_fields ) then
3317  write( text,'(i2)' ) max_fields
3318  call error_mesg ('xgrid_mod', 'MAX_FIELDS='//trim(text)//' exceeded for group get_side1_from_xgrid', fatal)
3319  endif
3320  d_addrs(lsize) = loc(d)
3321  x_addrs(lsize) = loc(x)
3322 
3323  if(lsize == 1) then
3324  isize = size(d,1)
3325  jsize = size(d,2)
3326  xsize = size(x(:))
3327  grid_id_saved = grid_id
3328  else
3329  set_mismatch = .false.
3330  set_mismatch = set_mismatch .OR. (isize /= size(d,1))
3331  set_mismatch = set_mismatch .OR. (jsize /= size(d,2))
3332  set_mismatch = set_mismatch .OR. (xsize /= size(x(:)))
3333  set_mismatch = set_mismatch .OR. (grid_id_saved /= grid_id)
3334  if(set_mismatch)then
3335  write( text,'(i2)' ) lsize
3336  call error_mesg ('xgrid_mod', 'Incompatible field at count '//text//' for group get_side1_from_xgrid', fatal )
3337  endif
3338  endif
3339 
3340  if(is_complete) then
3341  if (make_exchange_reproduce) then
3342  call get_1_from_xgrid_repro(d_addrs, x_addrs, xmap, xsize, lsize)
3343  else
3344  call get_1_from_xgrid(d_addrs, x_addrs, xmap, isize, jsize, xsize, lsize)
3345  end if
3346  d_addrs(1:lsize) = -9999
3347  x_addrs(1:lsize) = -9999
3348  isize = 0
3349  jsize = 0
3350  xsize = 0
3351  lsize = 0
3352  grid_id_saved = ""
3353  endif
3354  return;
3355  end if
3356 
3357  do g=2,size(xmap%grids(:))
3358  if (grid_id==xmap%grids(g)%id) &
3359  call error_mesg ('xgrid_mod', &
3360  'get_from_xgrid expects a 3D side 2 grid', fatal)
3361  end do
3362 
3363  call error_mesg ('xgrid_mod', 'get_from_xgrid: could not find grid id', fatal)
3364 
3365 end subroutine get_side1_from_xgrid
3366 ! </SUBROUTINE>
3367 
3368 !#######################################################################
3369 
3370 ! <SUBROUTINE NAME="get_side2_from_xgrid" INTERFACE="get_from_xgrid">
3371 ! <IN NAME="x" TYPE="real" DIM="(:)" > </IN>
3372 ! <IN NAME="grid_id" TYPE=" character(len=3)" > </IN>
3373 ! <OUT NAME="d" TYPE="real" DIM="(:,:,:)" > </OUT>
3374 ! <INOUT NAME="xmap" TYPE="xmap_type" > </INOUT>
3375 
3376 subroutine get_side2_from_xgrid(d, grid_id, x, xmap)
3377 real, dimension(:,:,:), intent( out) :: d
3378 character(len=3), intent(in ) :: grid_id
3379 real, dimension(:), intent(in ) :: x
3380 type(xmap_type), intent(in ) :: xmap
3381 
3382  integer :: g
3383 
3384  if (grid_id==xmap%grids(1)%id) &
3385  call error_mesg ('xgrid_mod', &
3386  'get_from_xgrid expects a 2D side 1 grid', fatal)
3387 
3388  do g=2,size(xmap%grids(:))
3389  if (grid_id==xmap%grids(g)%id) then
3390  call get_2_from_xgrid(d, xmap%grids(g), x, xmap)
3391  return;
3392  end if
3393  end do
3394 
3395  call error_mesg ('xgrid_mod', 'get_from_xgrid: could not find grid id', fatal)
3396 
3397 end subroutine get_side2_from_xgrid
3398 ! </SUBROUTINE>
3399 
3400 !#######################################################################
3401 
3402 ! <SUBROUTINE NAME="some">
3403 
3404 ! <OVERVIEW>
3405 ! Returns logical associating exchange grid cells with given side two grid.
3406 ! </OVERVIEW>
3407 ! <DESCRIPTION>
3408 ! Returns logical associating exchange grid cells with given side two grid.
3409 ! </DESCRIPTION>
3410 ! <TEMPLATE>
3411 ! call some(xmap, some_arr, grid_id)
3412 ! </TEMPLATE>
3413 
3414 ! <IN NAME="xmap" TYPE="xmap_type" ></IN>
3415 ! <IN NAME="grid_id" TYPE="character(len=3)" ></IN>
3416 ! <OUT NAME="some_arr" TYPE="logical" DIM="(xmap%size)" >
3417 ! logical associating exchange grid cells with given side 2 grid.
3418 ! </OUT>
3419 
3420 subroutine some(xmap, some_arr, grid_id)
3421 type(xmap_type), intent(in) :: xmap
3422 character(len=3), optional, intent(in) :: grid_id
3423 logical, dimension(:), intent(out) :: some_arr
3424 
3425  integer :: g
3426 
3427  if (.not.present(grid_id)) then
3428 
3429  if(xmap%size > 0) then
3430  some_arr = .true.
3431  else
3432  some_arr = .false.
3433  end if
3434  return;
3435  end if
3436 
3437  if (grid_id==xmap%grids(1)%id) &
3438  call error_mesg ('xgrid_mod', 'some expects a side 2 grid id', fatal)
3439 
3440  do g=2,size(xmap%grids(:))
3441  if (grid_id==xmap%grids(g)%id) then
3442  some_arr = .false.
3443  some_arr(xmap%grids(g)%first:xmap%grids(g)%last) = .true.;
3444  return;
3445  end if
3446  end do
3447 
3448  call error_mesg ('xgrid_mod', 'some could not find grid id', fatal)
3449 
3450 end subroutine some
3451 ! </SUBROUTINE>
3452 
3453 !#######################################################################
3454 
3455 subroutine put_2_to_xgrid(d, grid, x, xmap)
3456 type(grid_type), intent(in) :: grid
3457 real, dimension(grid%is_me:grid%ie_me, & grid%js_me:grid%je_me, grid%km), intent(in) :: d
3458 real, dimension(: ), intent(inout) :: x
3459 type(xmap_type), intent(in ) :: xmap
3460 
3461  integer :: l
3462  call mpp_clock_begin(id_put_2_to_xgrid)
3463 
3464  do l=grid%first,grid%last
3465  x(l) = d(xmap%x2(l)%i,xmap%x2(l)%j,xmap%x2(l)%k)
3466  end do
3467 
3468  call mpp_clock_end(id_put_2_to_xgrid)
3469 end subroutine put_2_to_xgrid
3470 
3471 !#######################################################################
3472 
3473 subroutine get_2_from_xgrid(d, grid, x, xmap)
3474 type(grid_type), intent(in ) :: grid
3475 real, dimension(grid%is_me:grid%ie_me, & grid%js_me:grid%je_me, grid%km), intent(out) :: d
3476 real, dimension(:), intent(in ) :: x
3477 type(xmap_type), intent(in ) :: xmap
3478 
3479  integer :: l, k
3480 
3481  call mpp_clock_begin(id_get_2_from_xgrid)
3482 
3483  d = 0.0
3484  do l=grid%first_get,grid%last_get
3485  d(xmap%x2_get(l)%i,xmap%x2_get(l)%j,xmap%x2_get(l)%k) = &
3486  d(xmap%x2_get(l)%i,xmap%x2_get(l)%j,xmap%x2_get(l)%k) + xmap%x2_get(l)%area*x(xmap%x2_get(l)%pos)
3487  end do
3488  !
3489  ! normalize with side 2 grid cell areas
3490  !
3491  do k=1,size(d,3)
3492  d(:,:,k) = d(:,:,k) * grid%area_inv
3493  end do
3494 
3495  call mpp_clock_end(id_get_2_from_xgrid)
3496 
3497 end subroutine get_2_from_xgrid
3498 
3499 !#######################################################################
3500 
3501 subroutine put_1_to_xgrid_order_1(d_addrs, x_addrs, xmap, isize, jsize, xsize, lsize)
3502  integer(LONG_KIND), dimension(:), intent(in) :: d_addrs
3503  integer(LONG_KIND), dimension(:), intent(in) :: x_addrs
3504  type(xmap_type), intent(inout) :: xmap
3505  integer, intent(in) :: isize, jsize, xsize, lsize
3506 
3507  integer :: i, j, p, buffer_pos, msgsize
3508  integer :: from_pe, to_pe, pos, n, l, count
3509  integer :: ibegin, istart, iend, start_pos
3510  type(comm_type), pointer, save :: comm =>null()
3511  real :: recv_buffer(xmap%put1%recvsize*lsize)
3512  real :: send_buffer(xmap%put1%sendsize*lsize)
3513  real :: unpack_buffer(xmap%put1%recvsize)
3514 
3515  real, dimension(isize, jsize) :: d
3516  real, dimension(xsize) :: x
3517  pointer(ptr_d, d)
3518  pointer(ptr_x, x)
3519 
3520  call mpp_clock_begin(id_put_1_to_xgrid_order_1)
3521 
3522  !--- pre-post receiving
3523  comm => xmap%put1
3524  do p = 1, comm%nrecv
3525  msgsize = comm%recv(p)%count*lsize
3526  from_pe = comm%recv(p)%pe
3527  buffer_pos = comm%recv(p)%buffer_pos*lsize
3528  call mpp_recv(recv_buffer(buffer_pos+1), glen=msgsize, from_pe = from_pe, block=.false., tag=comm_tag_7)
3529  enddo
3530 
3531  !--- send the data
3532  buffer_pos = 0
3533  do p = 1, comm%nsend
3534  msgsize = comm%send(p)%count*lsize
3535  to_pe = comm%send(p)%pe
3536  pos = buffer_pos
3537  do l = 1, lsize
3538  ptr_d = d_addrs(l)
3539  do n = 1, comm%send(p)%count
3540  pos = pos + 1
3541  i = comm%send(p)%i(n)
3542  j = comm%send(p)%j(n)
3543  send_buffer(pos) = d(i,j)
3544  enddo
3545  enddo
3546  call mpp_send(send_buffer(buffer_pos+1), plen=msgsize, to_pe = to_pe, tag=comm_tag_7 )
3547  buffer_pos = buffer_pos + msgsize
3548  enddo
3549 
3550  call mpp_sync_self(check=event_recv)
3551 
3552  !--- unpack the buffer
3553  if( lsize == 1) then
3554  ptr_x = x_addrs(1)
3555  do l=1,xmap%size_put1
3556  x(l) = recv_buffer(xmap%x1_put(l)%pos)
3557  end do
3558  else
3559  start_pos = 0
3560 !$OMP parallel do default(none) shared(lsize,x_addrs,comm,recv_buffer,xmap) &
3561 !$OMP private(ptr_x,count,ibegin,istart,iend,pos,unpack_buffer)
3562  do l = 1, lsize
3563  ptr_x = x_addrs(l)
3564  do p = 1, comm%nrecv
3565  count = comm%recv(p)%count
3566  ibegin = comm%recv(p)%buffer_pos*lsize + 1
3567  istart = ibegin + (l-1)*count
3568  iend = istart + count - 1
3569  pos = comm%recv(p)%buffer_pos
3570  do n = istart, iend
3571  pos = pos + 1
3572  unpack_buffer(pos) = recv_buffer(n)
3573  enddo
3574  enddo
3575  do i=1,xmap%size_put1
3576  x(i) = unpack_buffer(xmap%x1_put(i)%pos)
3577  end do
3578  enddo
3579  endif
3580 
3581  call mpp_sync_self()
3582 
3583  call mpp_clock_end(id_put_1_to_xgrid_order_1)
3584 
3585 end subroutine put_1_to_xgrid_order_1
3586 
3587 !#######################################################################
3588 
3589 
3590 subroutine put_1_to_xgrid_order_2(d_addrs, x_addrs, xmap, isize, jsize, xsize, lsize)
3591  integer(LONG_KIND), dimension(:), intent(in) :: d_addrs
3592  integer(LONG_KIND), dimension(:), intent(in) :: x_addrs
3593  type(xmap_type), intent(inout) :: xmap
3594  integer, intent(in) :: isize, jsize, xsize, lsize
3595 
3596  !: NOTE: halo size is assumed to be 1 in setup_xmap
3597  real, dimension(0:isize+1, 0:jsize+1, lsize) :: tmp
3598  real, dimension(isize, jsize, lsize) :: tmpx, tmpy
3599  real, dimension(isize, jsize, lsize) :: d_bar_max, d_bar_min
3600  real, dimension(isize, jsize, lsize) :: d_max, d_min
3601  real :: d_bar
3602  integer :: i, is, ie, im, j, js, je, jm, ii, jj
3603  integer :: p, l, ioff, joff, isd, jsd
3604  type(grid_type), pointer, save :: grid1 =>null()
3605  type(comm_type), pointer, save :: comm =>null()
3606  integer :: buffer_pos, msgsize, from_pe, to_pe, pos, n
3607  integer :: ibegin, count, istart, iend
3608  real :: recv_buffer(xmap%put1%recvsize*lsize*3)
3609  real :: send_buffer(xmap%put1%sendsize*lsize*3)
3610  real :: unpack_buffer(xmap%put1%recvsize*3)
3611  logical :: on_west_edge, on_east_edge, on_south_edge, on_north_edge
3612  real, dimension(isize, jsize) :: d
3613  real, dimension(xsize) :: x
3614  pointer(ptr_d, d)
3615  pointer(ptr_x, x)
3616 
3617  call mpp_clock_begin(id_put_1_to_xgrid_order_2)
3618  grid1 => xmap%grids(1)
3619 
3620  is = grid1%is_me; ie = grid1%ie_me
3621  js = grid1%js_me; je = grid1%je_me
3622  isd = grid1%isd_me
3623  jsd = grid1%jsd_me
3624 
3625 !$OMP parallel do default(none) shared(lsize,tmp,d_addrs,isize,jsize) private(ptr_d)
3626  do l = 1, lsize
3627  tmp(:,:,l) = large_number
3628  ptr_d = d_addrs(l)
3629  tmp(1:isize,1:jsize,l) = d(:,:)
3630  enddo
3631 
3632  if(grid1%is_latlon) then
3633  call mpp_update_domains(tmp,grid1%domain_with_halo)
3634 !$OMP parallel do default(none) shared(lsize,tmp,grid1,is,ie,js,je,isd,jsd,tmpx,tmpy)
3635  do l = 1, lsize
3636  tmpy(:,:,l) = grad_merid_latlon(tmp(:,:,l), grid1%lat, is, ie, js, je, isd, jsd)
3637  tmpx(:,:,l) = grad_zonal_latlon(tmp(:,:,l), grid1%lon, grid1%lat, is, ie, js, je, isd, jsd)
3638  enddo
3639  else
3640  call mpp_update_domains(tmp,grid1%domain)
3641  on_west_edge = (is==1)
3642  on_east_edge = (ie==grid1%im)
3643  on_south_edge = (js==1)
3644  on_north_edge = (je==grid1%jm)
3645 !$OMP parallel do default(none) shared(lsize,tmp,grid1,tmpx,tmpy, &
3646 !$OMP on_west_edge,on_east_edge,on_south_edge,on_north_edge)
3647  do l = 1, lsize
3648  call gradient_cubic(tmp(:,:,l), grid1%box%dx, grid1%box%dy, grid1%box%area, &
3649  grid1%box%edge_w, grid1%box%edge_e, grid1%box%edge_s, &
3650  grid1%box%edge_n, grid1%box%en1, grid1%box%en2, &
3651  grid1%box%vlon, grid1%box%vlat, tmpx(:,:,l), tmpy(:,:,l), &
3652  on_west_edge, on_east_edge, on_south_edge, on_north_edge)
3653  enddo
3654  end if
3655 
3656  !--- pre-post receiving
3657  buffer_pos = 0
3658  comm => xmap%put1
3659  do p = 1, comm%nrecv
3660  msgsize = comm%recv(p)%count*lsize
3661  buffer_pos = comm%recv(p)%buffer_pos*lsize
3662  if(.NOT. monotonic_exchange) then
3663  msgsize = msgsize*3
3664  buffer_pos = buffer_pos*3
3665  endif
3666  from_pe = comm%recv(p)%pe
3667  call mpp_recv(recv_buffer(buffer_pos+1), glen=msgsize, from_pe = from_pe, block=.false., tag=comm_tag_8)
3668  enddo
3669 
3670  !--- compute d_bar_max and d_bar_min.
3671  if(monotonic_exchange) then
3672 !$OMP parallel do default(none) shared(lsize,isize,jsize,d_bar_max,d_bar_min,d_max,d_min,tmp)
3673  do l = 1, lsize
3674  do j = 1, jsize
3675  do i = 1, isize
3676  d_bar_max(i,j,l) = -large_number
3677  d_bar_min(i,j,l) = large_number
3678  d_max(i,j,l) = -large_number
3679  d_min(i,j,l) = large_number
3680  do jj = j-1, j+1
3681  do ii = i-1, i+1
3682  if(tmp(i,j,l) .NE. large_number) then
3683  if(tmp(i,j,l) > d_bar_max(i,j,l)) d_bar_max(i,j,l) = tmp(i,j,l)
3684  if(tmp(i,j,l) < d_bar_min(i,j,l)) d_bar_min(i,j,l) = tmp(i,j,l)
3685  endif
3686  enddo
3687  enddo
3688  enddo
3689  enddo
3690  enddo
3691  endif
3692 
3693  !--- send the data
3694  buffer_pos = 0
3695  if(monotonic_exchange) then
3696  pos = 0
3697  do p = 1, comm%nsend
3698  msgsize = comm%send(p)%count*lsize
3699  to_pe = comm%send(p)%pe
3700  do l = 1, lsize
3701  ptr_d = d_addrs(l)
3702  do n = 1, comm%send(p)%count
3703  pos = pos + 1
3704  i = comm%send(p)%i(n)
3705  j = comm%send(p)%j(n)
3706  send_buffer(pos) = d(i,j) + tmpy(i,j,l)*comm%send(p)%dj(n) + tmpx(i,j,l)*comm%send(p)%di(n)
3707  if(send_buffer(pos) > d_max(i,j,l)) d_max(i,j,l) = send_buffer(pos)
3708  if(send_buffer(pos) < d_min(i,j,l)) d_min(i,j,l) = send_buffer(pos)
3709  enddo
3710  enddo
3711  enddo
3712 
3713  do p = 1, comm%nsend
3714  msgsize = comm%send(p)%count*lsize
3715  to_pe = comm%send(p)%pe
3716  pos = buffer_pos
3717  do l = 1, lsize
3718  ptr_d = d_addrs(l)
3719  do n = 1, comm%send(p)%count
3720  pos = pos + 1
3721  i = comm%send(p)%i(n)
3722  j = comm%send(p)%j(n)
3723  d_bar = d(i,j)
3724  if( d_max(i,j,l) > d_bar_max(i,j,l) ) then
3725  send_buffer(pos) = d_bar + ((send_buffer(pos)-d_bar)/(d_max(i,j,l)-d_bar)) * (d_bar_max(i,j,l)-d_bar)
3726  else if( d_min(i,j,l) < d_bar_min(i,j,l) ) then
3727  send_buffer(pos) = d_bar + ((send_buffer(pos)-d_bar)/(d_min(i,j,l)-d_bar)) * (d_bar_min(i,j,l)-d_bar)
3728  endif
3729  enddo
3730  enddo
3731  call mpp_send(send_buffer(buffer_pos+1), plen=msgsize, to_pe = to_pe, tag=comm_tag_8 )
3732  buffer_pos = buffer_pos + msgsize
3733  enddo
3734  else
3735  do p = 1, comm%nsend
3736  msgsize = comm%send(p)%count*lsize*3
3737  to_pe = comm%send(p)%pe
3738  pos = buffer_pos
3739  do l = 1, lsize
3740  ptr_d = d_addrs(l)
3741  do n = 1, comm%send(p)%count
3742  pos = pos + 3
3743  i = comm%send(p)%i(n)
3744  j = comm%send(p)%j(n)
3745  send_buffer(pos-2) = d(i,j)
3746  send_buffer(pos-1) = tmpy(i,j,l)
3747  send_buffer(pos ) = tmpx(i,j,l)
3748  enddo
3749  enddo
3750  call mpp_send(send_buffer(buffer_pos+1), plen=msgsize, to_pe = to_pe, tag=comm_tag_8 )
3751  buffer_pos = buffer_pos + msgsize
3752  enddo
3753  endif
3754 
3755  call mpp_sync_self(check=event_recv)
3756 
3757  !--- unpack the buffer
3758  if(monotonic_exchange) then
3759  if( lsize == 1) then
3760  ptr_x = x_addrs(1)
3761  do l=1,xmap%size_put1
3762  pos = xmap%x1_put(l)%pos
3763  x(l) = recv_buffer(pos)
3764  end do
3765  else
3766  do l = 1, lsize
3767  ptr_x = x_addrs(l)
3768  pos = 0
3769  do p = 1, comm%nsend
3770  count = comm%send(p)%count
3771  ibegin = comm%recv(p)%buffer_pos*lsize + 1
3772  istart = ibegin + (l-1)*count
3773  iend = istart + count - 1
3774  pos = comm%recv(p)%buffer_pos
3775  do n = istart, iend
3776  pos = pos + 1
3777  unpack_buffer(pos) = recv_buffer(n)
3778  enddo
3779  enddo
3780  do i=1,xmap%size_put1
3781  pos = xmap%x1_put(i)%pos
3782  x(i) = unpack_buffer(pos)
3783  end do
3784  enddo
3785  endif
3786  else
3787  if( lsize == 1) then
3788  ptr_x = x_addrs(1)
3789 !$OMP parallel do default(none) shared(xmap,recv_buffer,ptr_x) private(pos)
3790  do l=1,xmap%size_put1
3791  pos = xmap%x1_put(l)%pos
3792  x(l) = recv_buffer(3*pos-2) + recv_buffer(3*pos-1)*xmap%x1_put(l)%dj + recv_buffer(3*pos)*xmap%x1_put(l)%di
3793  end do
3794  else
3795 !$OMP parallel do default(none) shared(lsize,comm,xmap,recv_buffer,x_addrs) &
3796 !$OMP private(ptr_x,pos,ibegin,istart,iend,count,unpack_buffer)
3797  do l = 1, lsize
3798  ptr_x = x_addrs(l)
3799  pos = 0
3800  ibegin = 1
3801  do p = 1, comm%nrecv
3802  count = comm%recv(p)%count*3
3803  ibegin = comm%recv(p)%buffer_pos*lsize*3 + 1
3804  istart = ibegin + (l-1)*count
3805  iend = istart + count - 1
3806  pos = comm%recv(p)%buffer_pos*3
3807  do n = istart, iend
3808  pos = pos + 1
3809  unpack_buffer(pos) = recv_buffer(n)
3810  enddo
3811  enddo
3812  do i=1,xmap%size_put1
3813  pos = xmap%x1_put(i)%pos
3814  x(i) = unpack_buffer(3*pos-2) + unpack_buffer(3*pos-1)*xmap%x1_put(i)%dj + unpack_buffer(3*pos)*xmap%x1_put(i)%di
3815  end do
3816  enddo
3817  endif
3818  endif
3819 
3820  call mpp_sync_self()
3821  call mpp_clock_end(id_put_1_to_xgrid_order_2)
3822 
3823 end subroutine put_1_to_xgrid_order_2
3824 
3825 !#######################################################################
3826 
3827 subroutine get_1_from_xgrid(d_addrs, x_addrs, xmap, isize, jsize, xsize, lsize)
3828  integer(LONG_KIND), dimension(:), intent(in) :: d_addrs
3829  integer(LONG_KIND), dimension(:), intent(in) :: x_addrs
3830  type(xmap_type), intent(inout) :: xmap
3831  integer, intent(in) :: isize, jsize, xsize, lsize
3832 
3833  real, dimension(xmap%size), target :: dg(xmap%size, lsize)
3834  integer :: i, j, l, p, n, m
3835  integer :: msgsize, buffer_pos, pos
3836  integer :: istart, iend, count
3837  real , pointer, save :: dgp =>null()
3838  type(grid_type) , pointer, save :: grid1 =>null()
3839  type(comm_type) , pointer, save :: comm =>null()
3840  type(overlap_type), pointer, save :: send => null()
3841  type(overlap_type), pointer, save :: recv => null()
3842  real :: recv_buffer(xmap%get1%recvsize*lsize*3)
3843  real :: send_buffer(xmap%get1%sendsize*lsize*3)
3844  real :: unpack_buffer(xmap%get1%recvsize*3)
3845  real :: d(isize,jsize)
3846  real, dimension(xsize) :: x
3847  pointer(ptr_d, d)
3848  pointer(ptr_x, x)
3849 
3850  call mpp_clock_begin(id_get_1_from_xgrid)
3851 
3852  comm => xmap%get1
3853  grid1 => xmap%grids(1)
3854 
3855  do p = 1, comm%nrecv
3856  recv => comm%recv(p)
3857  msgsize = recv%count*lsize
3858  buffer_pos = recv%buffer_pos*lsize
3859  call mpp_recv(recv_buffer(buffer_pos+1), glen=msgsize, from_pe = recv%pe, block=.false., tag=comm_tag_9)
3860  enddo
3861 
3862  dg = 0.0;
3863 !$OMP parallel do default(none) shared(lsize,xmap,dg,x_addrs) private(dgp,ptr_x)
3864  do l = 1, lsize
3865  ptr_x = x_addrs(l)
3866  do i=1,xmap%size
3867  dgp => dg(xmap%x1(i)%pos,l)
3868  dgp = dgp + xmap%x1(i)%area*x(i)
3869  enddo
3870  enddo
3871 
3872 
3873  !--- send the data
3874  buffer_pos = 0
3875  istart = 1
3876  do p = 1, comm%nsend
3877  send => comm%send(p)
3878  msgsize = send%count*lsize
3879  pos = buffer_pos
3880  istart = send%buffer_pos+1
3881  iend = istart + send%count - 1
3882  do l = 1, lsize
3883  do n = istart, iend
3884  pos = pos + 1
3885  send_buffer(pos) = dg(n,l)
3886  enddo
3887  enddo
3888  call mpp_send(send_buffer(buffer_pos+1), plen=msgsize, to_pe = send%pe, tag=comm_tag_9 )
3889  buffer_pos = buffer_pos + msgsize
3890  istart = iend + 1
3891  enddo
3892 
3893  call mpp_sync_self(check=event_recv)
3894 
3895  !--- unpack the buffer
3896  do l = 1, lsize
3897  ptr_d = d_addrs(l)
3898  d = 0.0
3899  enddo
3900  !--- To bitwise reproduce old results, first copy the data onto its own pe.
3901 
3902  do p = 1, comm%nrecv
3903  recv => comm%recv(p)
3904  count = recv%count
3905  buffer_pos = recv%buffer_pos*lsize
3906  if( recv%pe == xmap%me ) then
3907 !$OMP parallel do default(none) shared(lsize,recv,recv_buffer,buffer_pos,d_addrs,count) &
3908 !$OMP private(ptr_d,i,j,pos)
3909  do l = 1, lsize
3910  pos = buffer_pos + (l-1)*count
3911  ptr_d = d_addrs(l)
3912  do n = 1,count
3913  i = recv%i(n)
3914  j = recv%j(n)
3915  pos = pos + 1
3916  d(i,j) = recv_buffer(pos)
3917  enddo
3918  enddo
3919  exit
3920  endif
3921  enddo
3922 
3923  pos = 0
3924  do m = 1, comm%nrecv
3925  p = comm%unpack_ind(m)
3926  recv => comm%recv(p)
3927  if( recv%pe == xmap%me ) then
3928  cycle
3929  endif
3930  buffer_pos = recv%buffer_pos*lsize
3931 !$OMP parallel do default(none) shared(lsize,recv,recv_buffer,buffer_pos,d_addrs) &
3932 !$OMP private(ptr_d,i,j,pos)
3933  do l = 1, lsize
3934  pos = buffer_pos + (l-1)*recv%count
3935  ptr_d = d_addrs(l)
3936  do n = 1, recv%count
3937  i = recv%i(n)
3938  j = recv%j(n)
3939  pos = pos + 1
3940  d(i,j) = d(i,j) + recv_buffer(pos)
3941  enddo
3942  enddo
3943  enddo
3944 
3945  !
3946  ! normalize with side 1 grid cell areas
3947  !
3948 !$OMP parallel do default(none) shared(lsize,d_addrs,grid1) private(ptr_d)
3949  do l = 1, lsize
3950  ptr_d = d_addrs(l)
3951  d = d * grid1%area_inv
3952  enddo
3953  call mpp_sync_self()
3954  call mpp_clock_end(id_get_1_from_xgrid)
3955 
3956 end subroutine get_1_from_xgrid
3957 
3958 !#######################################################################
3959 
3960 subroutine get_1_from_xgrid_repro(d_addrs, x_addrs, xmap, xsize, lsize)
3961  integer(LONG_KIND), dimension(:), intent(in) :: d_addrs
3962  integer(LONG_KIND), dimension(:), intent(in) :: x_addrs
3963  type(xmap_type), intent(inout) :: xmap
3964  integer, intent(in) :: xsize, lsize
3965 
3966  integer :: g, i, j, k, p, l, n, l2, m, l3
3967  integer :: msgsize, buffer_pos, pos
3968  type(grid_type), pointer, save :: grid =>null()
3969  type(comm_type), pointer, save :: comm => null()
3970  type(overlap_type), pointer, save :: send => null()
3971  type(overlap_type), pointer, save :: recv => null()
3972  integer, dimension(0:xmap%npes-1) :: pl, ml
3973  real :: recv_buffer(xmap%recv_count_repro_tot*lsize)
3974  real :: send_buffer(xmap%send_count_repro_tot*lsize)
3975  real :: d(xmap%grids(1)%is_me:xmap%grids(1)%ie_me, &
3976  xmap%grids(1)%js_me:xmap%grids(1)%je_me)
3977  real, dimension(xsize) :: x
3978  pointer(ptr_d, d)
3979  pointer(ptr_x, x)
3980 
3981  call mpp_clock_begin(id_get_1_from_xgrid_repro)
3982  comm => xmap%get1_repro
3983  !--- pre-post receiving
3984  do p = 1, comm%nrecv
3985  recv => comm%recv(p)
3986  msgsize = recv%count*lsize
3987  buffer_pos = recv%buffer_pos*lsize
3988  call mpp_recv(recv_buffer(buffer_pos+1), glen=msgsize, from_pe = recv%pe, block=.false., tag=comm_tag_10)
3989  n = recv%pe -xmap%root_pe
3990  pl(n) = buffer_pos
3991  ml(n) = recv%count
3992  enddo
3993 
3994  !pack the data
3995  send_buffer(:) = 0.0
3996 !$OMP parallel do default(none) shared(lsize,x_addrs,comm,xmap,send_buffer) &
3997 !$OMP private(ptr_x,i,j,g,l2,pos,send)
3998  do p = 1, comm%nsend
3999  pos = comm%send(p)%buffer_pos*lsize
4000  send => comm%send(p)
4001  do l = 1,lsize
4002  ptr_x = x_addrs(l)
4003  do n = 1, send%count
4004  i = send%i(n)
4005  j = send%j(n)
4006  g = send%g(n)
4007  l2 = send%xloc(n)
4008  pos = pos + 1
4009  do k =1, xmap%grids(g)%km
4010  if(xmap%grids(g)%frac_area(i,j,k)/=0.0) then
4011  l2 = l2+1
4012  send_buffer(pos) = send_buffer(pos) + xmap%x1(l2)%area *x(l2)
4013  endif
4014  enddo
4015  enddo
4016  enddo
4017  enddo
4018 
4019  do p =1, comm%nsend
4020  buffer_pos = comm%send(p)%buffer_pos*lsize
4021  msgsize = comm%send(p)%count*lsize
4022  call mpp_send(send_buffer(buffer_pos+1), plen=msgsize, to_pe=comm%send(p)%pe, tag=comm_tag_10)
4023  enddo
4024 
4025  do l = 1, lsize
4026  ptr_d = d_addrs(l)
4027  d = 0
4028  enddo
4029 
4030  call mpp_sync_self(check=event_recv)
4031 
4032 !$OMP parallel do default(none) shared(lsize,d_addrs,xmap,recv_buffer,pl,ml) &
4033 !$OMP private(ptr_d,grid,i,j,p,pos)
4034  do l = 1, lsize
4035  ptr_d = d_addrs(l)
4036  do g=2,size(xmap%grids(:))
4037  grid => xmap%grids(g)
4038  do l3=1,grid%size_repro ! index into side1 grid's patterns
4039  i = grid%x_repro(l3)%i1
4040  j = grid%x_repro(l3)%j1
4041  p = grid%x_repro(l3)%pe-xmap%root_pe
4042  pos = pl(p) + (l-1)*ml(p) + grid%x_repro(l3)%recv_pos
4043  d(i,j) = d(i,j) + recv_buffer(pos)
4044  end do
4045  end do
4046  ! normalize with side 1 grid cell areas
4047  d = d * xmap%grids(1)%area_inv
4048  enddo
4049 
4050  call mpp_sync_self()
4051 
4052  call mpp_clock_end(id_get_1_from_xgrid_repro)
4053 
4054 end subroutine get_1_from_xgrid_repro
4055 
4056 !#######################################################################
4057 
4058 ! <FUNCTION NAME="conservation_check_side1" INTERFACE="conservation_check">
4059 ! <IN NAME="d" TYPE="real" DIM="(:,:)" > </IN>
4060 ! <IN NAME="grid_id" TYPE="character(len=3)" > </IN>
4061 ! <INOUT NAME="xmap" TYPE="xmap_type" > </INOUT>
4062 ! <OUT NAME="conservation_check_side1" TYPE="real" DIM="dimension(3)" > </OUT>
4063 ! <IN NAME="remap_method" TYPE="integer,optional"></IN>
4064 ! conservation_check - returns three numbers which are the global sum of a
4065 ! variable (1) on its home model grid, (2) after interpolation to the other
4066 ! side grid(s), and (3) after re_interpolation back onto its home side grid(s).
4067 !
4068 function conservation_check_side1(d, grid_id, xmap,remap_method) ! this one for 1->2->1
4069 real, dimension(:,:), intent(in ) :: d
4070 character(len=3), intent(in ) :: grid_id
4071 type(xmap_type), intent(inout) :: xmap
4072 real, dimension(3) :: conservation_check_side1
4073 integer, intent(in), optional :: remap_method
4074 
4075 
4076  real, dimension(xmap%size) :: x_over, x_back
4077  real, dimension(size(d,1),size(d,2)) :: d1
4078  real, dimension(:,:,:), allocatable :: d2
4079  integer :: g
4080  type(grid_type), pointer, save :: grid1 =>null(), grid2 =>null()
4081 
4082  grid1 => xmap%grids(1)
4084  if(grid1%tile_me .NE. tile_nest) conservation_check_side1(1) = sum(grid1%area*d)
4085 ! if(grid1%tile_me .NE. tile_parent .OR. grid1%id .NE. "ATM") &
4086 ! conservation_check_side1(1) = sum(grid1%area*d)
4087 
4088  call put_to_xgrid (d, grid1%id, x_over, xmap, remap_method) ! put from side 1
4089  do g=2,size(xmap%grids(:))
4090  grid2 => xmap%grids(g)
4091  if(grid2%on_this_pe) then
4092  allocate (d2(grid2%is_me:grid2%ie_me, grid2%js_me:grid2%je_me, grid2%km) )
4093  endif
4094  call get_from_xgrid (d2, grid2%id, x_over, xmap) ! get onto side 2's
4095  if(grid2%on_this_pe) then
4096  conservation_check_side1(2) = conservation_check_side1(2) + sum( grid2%area * sum(grid2%frac_area*d2,dim=3) )
4097  endif
4098  call put_to_xgrid (d2, grid2%id, x_back, xmap) ! put from side 2's
4099  if(allocated(d2))deallocate (d2)
4100  end do
4101  call get_from_xgrid(d1, grid1%id, x_back, xmap) ! get onto side 1
4102  if(grid1%tile_me .NE. tile_nest) conservation_check_side1(3) = sum(grid1%area*d1)
4103 ! if(grid1%tile_me .NE. tile_parent .OR. grid1%id .NE. "ATM") &
4104 ! conservation_check_side1(3) = sum(grid1%area*d1)
4106 
4107 end function conservation_check_side1
4108 ! </FUNCTION>
4109 
4110 !#######################################################################
4111 !
4112 ! conservation_check - returns three numbers which are the global sum of a
4113 ! variable (1) on its home model grid, (2) after interpolation to the other
4114 ! side grid(s), and (3) after re_interpolation back onto its home side grid(s).
4115 !
4116 ! <FUNCTION NAME="conservation_check_side2" INTERFACE="conservation_check">
4117 ! <IN NAME="d" TYPE="real" DIM="(:,:,:)" > </IN>
4118 ! <IN NAME="grid_id" TYPE="character(len=3)" > </IN>
4119 ! <INOUT NAME="xmap" TYPE="xmap_type" > </INOUT>
4120 ! <OUT NAME="conservation_check_side2" TYPE="real" DIM="dimension(3)" > </OUT>
4121 
4122 function conservation_check_side2(d, grid_id, xmap,remap_method) ! this one for 2->1->2
4123 real, dimension(:,:,:), intent(in ) :: d
4124 character(len=3), intent(in ) :: grid_id
4125 type(xmap_type), intent(inout) :: xmap
4126 real, dimension(3) :: conservation_check_side2
4127 integer, intent(in), optional :: remap_method
4128 
4129 
4130  real, dimension(xmap%size) :: x_over, x_back
4131  real, dimension(:,: ), allocatable :: d1
4132  real, dimension(:,:,:), allocatable :: d2
4133  integer :: g
4134  type(grid_type), pointer, save :: grid1 =>null(), grid2 =>null()
4135 
4136  grid1 => xmap%grids(1)
4138  do g = 2,size(xmap%grids(:))
4139  grid2 => xmap%grids(g)
4140  if (grid_id==grid2%id) then
4141  if(grid2%on_this_pe) then
4142  conservation_check_side2(1) = sum( grid2%area * sum(grid2%frac_area*d,dim=3) )
4143  endif
4144  call put_to_xgrid(d, grid_id, x_over, xmap) ! put from this side 2
4145  else
4146  call put_to_xgrid(0.0 * grid2%frac_area, grid2%id, x_over, xmap) ! zero rest
4147  end if
4148  end do
4149 
4150  allocate ( d1(size(grid1%area,1),size(grid1%area,2)) )
4151  call get_from_xgrid(d1, grid1%id, x_over, xmap) ! get onto side 1
4152  if(grid1%tile_me .NE. tile_nest)conservation_check_side2(2) = sum(grid1%area*d1)
4153  call put_to_xgrid(d1, grid1%id, x_back, xmap,remap_method) ! put from side 1
4154  deallocate ( d1 )
4155 
4156  conservation_check_side2(3) = 0.0;
4157  do g = 2,size(xmap%grids(:))
4158  grid2 => xmap%grids(g)
4159  if(grid2%on_this_pe) then
4160  allocate ( d2( size(grid2%frac_area, 1), size(grid2%frac_area, 2), &
4161  size(grid2%frac_area, 3) ) )
4162  endif
4163  call get_from_xgrid(d2, grid2%id, x_back, xmap) ! get onto side 2's
4164  conservation_check_side2(3) = conservation_check_side2(3) + sum( grid2%area * sum(grid2%frac_area*d2,dim=3) )
4165  if(allocated(d2) )deallocate ( d2 )
4166  end do
4168 
4169 end function conservation_check_side2
4170 ! </FUNCTION>
4171 
4172 !#######################################################################
4173 
4174 ! <FUNCTION NAME="conservation_check_ug_side1" INTERFACE="conservation_check_ug">
4175 ! <IN NAME="d" TYPE="real" DIM="(:,:)" > </IN>
4176 ! <IN NAME="grid_id" TYPE="character(len=3)" > </IN>
4177 ! <INOUT NAME="xmap" TYPE="xmap_type" > </INOUT>
4178 ! <OUT NAME="conservation_check_ug_side1" TYPE="real" DIM="dimension(3)" > </OUT>
4179 ! <IN NAME="remap_method" TYPE="integer,optional"></IN>
4180 ! conservation_check_ug - returns three numbers which are the global sum of a
4181 ! variable (1) on its home model grid, (2) after interpolation to the other
4182 ! side grid(s), and (3) after re_interpolation back onto its home side grid(s).
4183 !
4184 function conservation_check_ug_side1(d, grid_id, xmap,remap_method) ! this one for 1->2->1
4185 real, dimension(:,:), intent(in ) :: d
4186 character(len=3), intent(in ) :: grid_id
4187 type(xmap_type), intent(inout) :: xmap
4188 real, dimension(3) :: conservation_check_ug_side1
4189 integer, intent(in), optional :: remap_method
4190 
4191 
4192  real, dimension(xmap%size) :: x_over, x_back
4193  real, dimension(size(d,1),size(d,2)) :: d1
4194  real, dimension(:,:,:), allocatable :: d2
4195  real, dimension(: ), allocatable :: d_ug
4196  real, dimension(:,:), allocatable :: d2_ug
4197  integer :: g
4198  type(grid_type), pointer, save :: grid1 =>null(), grid2 =>null()
4199 
4200  grid1 => xmap%grids(1)
4202 
4203 
4204  if(grid1%is_ug) then
4205  allocate(d_ug(grid1%ls_me:grid1%le_me))
4206  call mpp_pass_sg_to_ug(grid1%ug_domain, d, d_ug)
4207  if(grid1%tile_me .NE. tile_nest) conservation_check_ug_side1(1) = sum(grid1%area(:,1)*d_ug)
4208  call put_to_xgrid_ug (d_ug, grid1%id, x_over, xmap) ! put from side 1
4209  else
4210  if(grid1%tile_me .NE. tile_nest) conservation_check_ug_side1(1) = sum(grid1%area*d)
4211  call put_to_xgrid (d, grid1%id, x_over, xmap, remap_method) ! put from side 1
4212  endif
4213  do g=2,size(xmap%grids(:))
4214  grid2 => xmap%grids(g)
4215  if(grid2%is_ug) then
4216  if(grid2%on_this_pe) then
4217  allocate (d2_ug(grid2%ls_me:grid2%le_me, grid2%km) )
4218  d2_ug = 0
4219  endif
4220  call get_from_xgrid_ug (d2_ug, grid2%id, x_over, xmap) ! get onto side 2's
4221  if(grid2%on_this_pe) then
4223  sum( grid2%area(:,1) * sum(grid2%frac_area(:,1,:)*d2_ug,dim=2) )
4224  endif
4225  call put_to_xgrid_ug (d2_ug, grid2%id, x_back, xmap) ! put from side 2's
4226  if(allocated(d2_ug))deallocate (d2_ug)
4227  else
4228  if(grid2%on_this_pe) then
4229  allocate (d2(grid2%is_me:grid2%ie_me, grid2%js_me:grid2%je_me, grid2%km) )
4230  endif
4231  call get_from_xgrid (d2, grid2%id, x_over, xmap) ! get onto side 2's
4232  if(grid2%on_this_pe) then
4233  conservation_check_ug_side1(2) = conservation_check_ug_side1(2) + sum( grid2%area * sum(grid2%frac_area*d2,dim=3) )
4234  endif
4235  call put_to_xgrid (d2, grid2%id, x_back, xmap) ! put from side 2's
4236  if(allocated(d2))deallocate (d2)
4237  endif
4238  end do
4239  if(grid1%is_ug) then
4240 ! call get_from_xgrid_ug(d_ug, grid1%id, x_back, xmap) ! get onto side 1
4241  if(grid1%tile_me .NE. tile_nest) conservation_check_ug_side1(3) = sum(grid1%area(:,1)*d_ug)
4242  else
4243  call get_from_xgrid(d1, grid1%id, x_back, xmap) ! get onto side 1
4244  if(grid1%tile_me .NE. tile_nest) conservation_check_ug_side1(3) = sum(grid1%area*d1)
4245  endif
4246  if(allocated(d_ug)) deallocate(d_ug)
4248 
4249 end function conservation_check_ug_side1
4250 ! </FUNCTION>
4251 
4252 !#######################################################################
4253 !
4254 ! conservation_check_ug - returns three numbers which are the global sum of a
4255 ! variable (1) on its home model grid, (2) after interpolation to the other
4256 ! side grid(s), and (3) after re_interpolation back onto its home side grid(s).
4257 !
4258 ! <FUNCTION NAME="conservation_check_ug_side2" INTERFACE="conservation_check_ug">
4259 ! <IN NAME="d" TYPE="real" DIM="(:,:,:)" > </IN>
4260 ! <IN NAME="grid_id" TYPE="character(len=3)" > </IN>
4261 ! <INOUT NAME="xmap" TYPE="xmap_type" > </INOUT>
4262 ! <OUT NAME="conservation_check_ug_side2" TYPE="real" DIM="dimension(3)" > </OUT>
4263 
4264 function conservation_check_ug_side2(d, grid_id, xmap,remap_method) ! this one for 2->1->2
4265 real, dimension(:,:,:), intent(in ) :: d
4266 character(len=3), intent(in ) :: grid_id
4267 type(xmap_type), intent(inout) :: xmap
4268 real, dimension(3) :: conservation_check_ug_side2
4269 integer, intent(in), optional :: remap_method
4270 
4271 
4272  real, dimension(xmap%size) :: x_over, x_back
4273  real, dimension(:,: ), allocatable :: d1, d_ug
4274  real, dimension(:,:,:), allocatable :: d2
4275  integer :: g
4276  type(grid_type), pointer, save :: grid1 =>null(), grid2 =>null()
4277 
4278  grid1 => xmap%grids(1)
4280  do g = 2,size(xmap%grids(:))
4281  grid2 => xmap%grids(g)
4282  if (grid_id==grid2%id) then
4283  if(grid2%on_this_pe) then
4284  if(grid2%is_ug) then
4285  allocate(d_ug(grid2%ls_me:grid2%le_me,grid2%km))
4286  call mpp_pass_sg_to_ug(grid2%ug_domain, d, d_ug)
4287  conservation_check_ug_side2(1) = sum( grid2%area(:,1) * sum(grid2%frac_area(:,1,:)*d_ug,dim=2) )
4288  else
4289  conservation_check_ug_side2(1) = sum( grid2%area(:,:) * sum(grid2%frac_area(:,:,:)*d,dim=3) )
4290  endif
4291  endif
4292  if(grid2%is_ug) then
4293  call put_to_xgrid_ug(d_ug, grid_id, x_over, xmap) ! put from this side 2
4294  else
4295  call put_to_xgrid(d, grid_id, x_over, xmap) ! put from this side 2
4296  endif
4297  if(allocated(d_ug)) deallocate(d_ug)
4298  else
4299  if(grid2%is_ug) then
4300  call put_to_xgrid_ug(0.0 * grid2%frac_area(:,1,:), grid2%id, x_over, xmap) ! zero rest
4301  else
4302  call put_to_xgrid(0.0 * grid2%frac_area, grid2%id, x_over, xmap) ! zero rest
4303  endif
4304  end if
4305  end do
4306 
4307  allocate ( d1(size(grid1%area,1),size(grid1%area,2)) )
4308  if(grid1%is_ug) then
4309  call get_from_xgrid_ug(d1(:,1), grid1%id, x_over, xmap) ! get onto side 1
4310  else
4311  call get_from_xgrid(d1, grid1%id, x_over, xmap) ! get onto side 1
4312  endif
4313  if(grid1%tile_me .NE. tile_nest)conservation_check_ug_side2(2) = sum(grid1%area*d1)
4314  if(grid1%is_ug) then
4315  call put_to_xgrid_ug(d1(:,1), grid1%id, x_back, xmap) ! put from side 1
4316  else
4317  call put_to_xgrid(d1, grid1%id, x_back, xmap,remap_method) ! put from side 1
4318  endif
4319  deallocate ( d1 )
4320 
4321  conservation_check_ug_side2(3) = 0.0;
4322  do g = 2,size(xmap%grids(:))
4323  grid2 => xmap%grids(g)
4324  if(grid2%on_this_pe) then
4325  allocate ( d2( size(grid2%frac_area, 1), size(grid2%frac_area, 2), &
4326  size(grid2%frac_area, 3) ) )
4327  endif
4328  if(grid2%is_ug) then
4329  call get_from_xgrid_ug(d2(:,1,:), grid2%id, x_back, xmap) ! get onto side 2's
4330  else
4331  call get_from_xgrid(d2, grid2%id, x_back, xmap) ! get onto side 2's
4332  endif
4333  conservation_check_ug_side2(3) = conservation_check_ug_side2(3) + sum( grid2%area * sum(grid2%frac_area*d2,dim=3) )
4334  if(allocated(d2) )deallocate ( d2 )
4335  end do
4337 
4338 end function conservation_check_ug_side2
4339 ! </FUNCTION>
4340 
4341 
4342 !******************************************************************************
4343 ! This routine is used to get the grid area of component model with id.
4344 subroutine get_xmap_grid_area(id, xmap, area)
4345  character(len=3), intent(in ) :: id
4346  type(xmap_type), intent(inout) :: xmap
4347  real, dimension(:,:), intent(out ) :: area
4348  integer :: g
4349  logical :: found
4350 
4351  found = .false.
4352  do g = 1, size(xmap%grids(:))
4353  if (id==xmap%grids(g)%id ) then
4354  if(size(area,1) .NE. size(xmap%grids(g)%area,1) .OR. size(area,2) .NE. size(xmap%grids(g)%area,2) ) &
4355  call error_mesg("xgrid_mod", "size mismatch between area and xmap%grids(g)%area", fatal)
4356  area = xmap%grids(g)%area
4357  found = .true.
4358  exit
4359  end if
4360  end do
4361 
4362  if(.not. found) call error_mesg("xgrid_mod", id//" is not found in xmap%grids id", fatal)
4363 
4364 end subroutine get_xmap_grid_area
4365 
4366 !#######################################################################
4367 
4368 ! This function is used to calculate the gradient along zonal direction.
4369 ! Maybe need to setup a limit for the gradient. The grid is assumeed
4370 ! to be regular lat-lon grid
4371 
4372 function grad_zonal_latlon(d, lon, lat, is, ie, js, je, isd, jsd)
4373 
4374  integer, intent(in) :: isd, jsd
4375  real, dimension(isd:,jsd:), intent(in) :: d
4376  real, dimension(:), intent(in) :: lon
4377  real, dimension(:), intent(in) :: lat
4378  integer, intent(in) :: is, ie, js, je
4379  real, dimension(is:ie,js:je) :: grad_zonal_latlon
4380  real :: dx, costheta
4381  integer :: i, j, ip1, im1
4382 
4383  ! calculate the gradient of the data on each grid
4384  do i = is, ie
4385  if(i == 1) then
4386  ip1 = i+1; im1 = i
4387  else if(i==size(lon(:)) ) then
4388  ip1 = i; im1 = i-1
4389  else
4390  ip1 = i+1; im1 = i-1
4391  endif
4392  dx = lon(ip1) - lon(im1)
4393  if(abs(dx).lt.eps ) call error_mesg('xgrids_mod(grad_zonal_latlon)', 'Improper grid size in lontitude', fatal)
4394  if(dx .gt. pi) dx = dx - 2.0* pi
4395  if(dx .lt. -pi) dx = dx + 2.0* pi
4396  do j = js, je
4397  costheta = cos(lat(j))
4398  if(abs(costheta) .lt. eps) call error_mesg('xgrids_mod(grad_zonal_latlon)', 'Improper latitude grid', fatal)
4399  grad_zonal_latlon(i,j) = (d(ip1,j)-d(im1,j))/(dx*costheta)
4400  enddo
4401  enddo
4402 
4403  return
4404 
4405 end function grad_zonal_latlon
4406 
4407 !#######################################################################
4408 
4409 ! This function is used to calculate the gradient along meridinal direction.
4410 ! Maybe need to setup a limit for the gradient. regular lat-lon grid are assumed
4411 
4412 function grad_merid_latlon(d, lat, is, ie, js, je, isd, jsd)
4413  integer, intent(in) :: isd, jsd
4414  real, dimension(isd:,jsd:), intent(in) :: d
4415  real, dimension(:), intent(in) :: lat
4416  integer, intent(in) :: is, ie, js, je
4417  real, dimension(is:ie,js:je) :: grad_merid_latlon
4418  real :: dy
4419  integer :: i, j, jp1, jm1
4420 
4421  ! calculate the gradient of the data on each grid
4422  do j = js, je
4423  if(j == 1) then
4424  jp1 = j+1; jm1 = j
4425  else if(j == size(lat(:)) ) then
4426  jp1 = j; jm1 = j-1
4427  else
4428  jp1 = j+1; jm1 = j-1
4429  endif
4430  dy = lat(jp1) - lat(jm1)
4431  if(abs(dy).lt.eps) call error_mesg('xgrids_mod(grad_merid_latlon)', 'Improper grid size in latitude', fatal)
4432 
4433  do i = is, ie
4434  grad_merid_latlon(i,j) = (d(i,jp1) - d(i,jm1))/dy
4435  enddo
4436  enddo
4437 
4438  return
4439 end function grad_merid_latlon
4440 
4441 !#######################################################################
4442 subroutine get_index_range(xmap, grid_index, is, ie, js, je, km)
4443 
4444  type(xmap_type), intent(in) :: xmap
4445  integer, intent(in) :: grid_index
4446  integer, intent(out) :: is, ie, js, je, km
4447 
4448  is = xmap % grids(grid_index) % is_me
4449  ie = xmap % grids(grid_index) % ie_me
4450  js = xmap % grids(grid_index) % js_me
4451  je = xmap % grids(grid_index) % je_me
4452  km = xmap % grids(grid_index) % km
4453 
4454 end subroutine get_index_range
4455 !#######################################################################
4456 
4457 subroutine stock_move_3d(from, to, grid_index, data, xmap, &
4458  & delta_t, from_side, to_side, radius, verbose, ier)
4459 
4460  ! this version takes rank 3 data, it can be used to compute the flux on anything but the
4461  ! first grid, which typically is on the atmos side.
4462  ! note that "from" and "to" are optional, the stocks will be subtracted, resp. added, only
4463  ! if these are present.
4464 
4465  use mpp_mod, only : mpp_sum
4467 
4468  type(stock_type), intent(inout), optional :: from, to
4469  integer, intent(in) :: grid_index ! grid index
4470  real, intent(in) :: data(:,:,:) ! data array is 3d
4471  type(xmap_type), intent(in) :: xmap
4472  real, intent(in) :: delta_t
4473  integer, intent(in) :: from_side, to_side ! ISTOCK_TOP, ISTOCK_BOTTOM, or ISTOCK_SIDE
4474  real, intent(in) :: radius ! earth radius
4475  character(len=*), intent(in), optional :: verbose
4476  integer, intent(out) :: ier
4477 
4478  real :: from_dq, to_dq
4479 
4480  ier = 0
4481  if(grid_index == 1) then
4482  ! data has rank 3 so grid index must be > 1
4483  ier = 1
4484  return
4485  endif
4486 
4487  if(.not. associated(xmap%grids) ) then
4488  ier = 2
4489  return
4490  endif
4491 
4492  from_dq = delta_t * 4.0*pi*radius**2 * sum( sum(xmap%grids(grid_index)%area * &
4493  & sum(xmap%grids(grid_index)%frac_area * data, dim=3), dim=1))
4494  to_dq = from_dq
4495 
4496  ! update only if argument is present.
4497  if(present(to )) to % dq( to_side) = to % dq( to_side) + to_dq
4498  if(present(from)) from % dq(from_side) = from % dq(from_side) - from_dq
4499 
4500  if(present(verbose).and.debug_stocks) then
4501  call mpp_sum(from_dq)
4502  call mpp_sum(to_dq)
4503  from_dq = from_dq/(4.0*pi*radius**2)
4504  to_dq = to_dq /(4.0*pi*radius**2)
4505  if(mpp_pe()==mpp_root_pe()) then
4506  write(stocks_file,'(a,es19.12,a,es19.12,a)') verbose, from_dq,' [*/m^2]'
4507  endif
4508  endif
4509 
4510 end subroutine stock_move_3d
4511 
4512 !...................................................................
4513 
4514 subroutine stock_move_2d(from, to, grid_index, data, xmap, &
4515  & delta_t, from_side, to_side, radius, verbose, ier)
4516 
4517  ! this version takes rank 2 data, it can be used to compute the flux on the atmos side
4518  ! note that "from" and "to" are optional, the stocks will be subtracted, resp. added, only
4519  ! if these are present.
4520 
4521  use mpp_mod, only : mpp_sum
4523 
4524  type(stock_type), intent(inout), optional :: from, to
4525  integer, optional, intent(in) :: grid_index
4526  real, intent(in) :: data(:,:) ! data array is 2d
4527  type(xmap_type), intent(in) :: xmap
4528  real, intent(in) :: delta_t
4529  integer, intent(in) :: from_side, to_side ! ISTOCK_TOP, ISTOCK_BOTTOM, or ISTOCK_SIDE
4530  real, intent(in) :: radius ! earth radius
4531  character(len=*), intent(in) :: verbose
4532  integer, intent(out) :: ier
4533 
4534  real :: to_dq, from_dq
4535 
4536  ier = 0
4537 
4538  if(.not. associated(xmap%grids) ) then
4539  ier = 3
4540  return
4541  endif
4542 
4543  if( .not. present(grid_index) .or. grid_index==1 ) then
4544 
4545  ! only makes sense if grid_index == 1
4546  from_dq = delta_t * 4.0*pi*radius**2 * sum(sum(xmap%grids(1)%area * data, dim=1))
4547  to_dq = from_dq
4548 
4549  else
4550 
4551  ier = 4
4552  return
4553 
4554  endif
4555 
4556  ! update only if argument is present.
4557  if(present(to )) to % dq( to_side) = to % dq( to_side) + to_dq
4558  if(present(from)) from % dq(from_side) = from % dq(from_side) - from_dq
4559 
4560  if(debug_stocks) then
4561  call mpp_sum(from_dq)
4562  call mpp_sum(to_dq)
4563  from_dq = from_dq/(4.0*pi*radius**2)
4564  to_dq = to_dq /(4.0*pi*radius**2)
4565  if(mpp_pe()==mpp_root_pe()) then
4566  write(stocks_file,'(a,es19.12,a,es19.12,a)') verbose, from_dq,' [*/m^2]'
4567  endif
4568  endif
4569 
4570 end subroutine stock_move_2d
4571 
4572 !#######################################################################
4573 
4574 subroutine stock_move_ug_3d(from, to, grid_index, data, xmap, &
4575  & delta_t, from_side, to_side, radius, verbose, ier)
4576 
4577  ! this version takes rank 3 data, it can be used to compute the flux on anything but the
4578  ! first grid, which typically is on the atmos side.
4579  ! note that "from" and "to" are optional, the stocks will be subtracted, resp. added, only
4580  ! if these are present.
4581 
4582  use mpp_mod, only : mpp_sum
4584 
4585  type(stock_type), intent(inout), optional :: from, to
4586  integer, intent(in) :: grid_index ! grid index
4587  real, intent(in) :: data(:,:) ! data array is 3d
4588  type(xmap_type), intent(in) :: xmap
4589  real, intent(in) :: delta_t
4590  integer, intent(in) :: from_side, to_side ! ISTOCK_TOP, ISTOCK_BOTTOM, or ISTOCK_SIDE
4591  real, intent(in) :: radius ! earth radius
4592  character(len=*), intent(in), optional :: verbose
4593  integer, intent(out) :: ier
4594  real, dimension(size(data,1),size(data,2)) :: tmp
4595 
4596  real :: from_dq, to_dq
4597 
4598  ier = 0
4599  if(grid_index == 1) then
4600  ! data has rank 3 so grid index must be > 1
4601  ier = 1
4602  return
4603  endif
4604 
4605  if(.not. associated(xmap%grids) ) then
4606  ier = 2
4607  return
4608  endif
4609 
4610  tmp = xmap%grids(grid_index)%frac_area(:,1,:) * data
4611  from_dq = delta_t * 4.0*pi*radius**2 * sum( xmap%grids(grid_index)%area(:,1) * &
4612  & sum(tmp, dim=2))
4613  to_dq = from_dq
4614 
4615  ! update only if argument is present.
4616  if(present(to )) to % dq( to_side) = to % dq( to_side) + to_dq
4617  if(present(from)) from % dq(from_side) = from % dq(from_side) - from_dq
4618 
4619  if(present(verbose).and.debug_stocks) then
4620  call mpp_sum(from_dq)
4621  call mpp_sum(to_dq)
4622  from_dq = from_dq/(4.0*pi*radius**2)
4623  to_dq = to_dq /(4.0*pi*radius**2)
4624  if(mpp_pe()==mpp_root_pe()) then
4625  write(stocks_file,'(a,es19.12,a,es19.12,a)') verbose, from_dq,' [*/m^2]'
4626  endif
4627  endif
4628 
4629 end subroutine stock_move_ug_3d
4630 
4631 
4632 
4633 !#######################################################################
4634 subroutine stock_integrate_2d(data, xmap, delta_t, radius, res, ier)
4635 
4636  ! surface/time integral of a 2d array
4638  use mpp_mod, only : mpp_sum
4639 
4640  real, intent(in) :: data(:,:) ! data array is 2d
4641  type(xmap_type), intent(in) :: xmap
4642  real, intent(in) :: delta_t
4643  real, intent(in) :: radius ! earth radius
4644  real, intent(out) :: res
4645  integer, intent(out) :: ier
4646 
4647  ier = 0
4648  res = 0.0
4649 
4650  if(.not. associated(xmap%grids) ) then
4651  ier = 6
4652  return
4653  endif
4654 
4655  res = delta_t * 4.0*pi*radius**2 * sum(sum(xmap%grids(1)%area * data, dim=1))
4656 
4657 end subroutine stock_integrate_2d
4658 !#######################################################################
4659 
4660 !#######################################################################
4661 
4662 
4663 
4664 subroutine stock_print(stck, Time, comp_name, index, ref_value, radius, pelist)
4665 
4666  use mpp_mod, only : mpp_pe, mpp_root_pe, mpp_sum
4669 
4670  type(stock_type), intent(in) :: stck
4671  type(time_type), intent(in) :: time
4672  character(len=*) :: comp_name
4673  integer, intent(in) :: index ! to map stock element (water, heat, ..) to a name
4674  real, intent(in) :: ref_value ! the stock value returned by the component per PE
4675  real, intent(in) :: radius
4676  integer, intent(in), optional :: pelist(:)
4677 
4678  integer, parameter :: initid = -2 ! initial value for diag IDs. Must not be equal to the value
4679  ! that register_diag_field returns when it can't register the filed -- otherwise the registration
4680  ! is attempted every time this subroutine is called
4681 
4682  real :: f_value, c_value, planet_area
4683  character(len=80) :: formatstring
4684  integer :: iday, isec, hours
4685  integer :: diagid, compind
4686  integer, dimension(NELEMS,4), save :: f_valuediagid = initid
4687  integer, dimension(NELEMS,4), save :: c_valuediagid = initid
4688  integer, dimension(NELEMS,4), save :: fmc_valuediagid = initid
4689 
4690  real :: diagfield
4691  logical :: used
4692  character(len=30) :: field_name, units
4693 
4694  f_value = sum(stck % dq)
4695  c_value = ref_value - stck % q_start
4696  if(present(pelist)) then
4697  call mpp_sum(f_value, pelist=pelist)
4698  call mpp_sum(c_value, pelist=pelist)
4699  else
4700  call mpp_sum(f_value)
4701  call mpp_sum(c_value)
4702  endif
4703 
4704  if(mpp_pe() == mpp_root_pe()) then
4705  ! normalize to 1 earth m^2
4706  planet_area = 4.0*pi*radius**2
4707  f_value = f_value / planet_area
4708  c_value = c_value / planet_area
4709 
4710  if(comp_name == 'ATM') compind = 1
4711  if(comp_name == 'LND') compind = 2
4712  if(comp_name == 'ICE') compind = 3
4713  if(comp_name == 'OCN') compind = 4
4714 
4715 
4716  if(f_valuediagid(index,compind) == initid) then
4717  field_name = trim(comp_name) // trim(stock_names(index))
4718  field_name = trim(field_name) // 'StocksChange_Flux'
4719  units = trim(stock_units(index))
4720  f_valuediagid(index,compind) = register_diag_field('stock_print', field_name, time, &
4721  units=units)
4722  endif
4723 
4724  if(c_valuediagid(index,compind) == initid) then
4725  field_name = trim(comp_name) // trim(stock_names(index))
4726  field_name = trim(field_name) // 'StocksChange_Comp'
4727  units = trim(stock_units(index))
4728  c_valuediagid(index,compind) = register_diag_field('stock_print', field_name, time, &
4729  units=units)
4730  endif
4731 
4732  if(fmc_valuediagid(index,compind) == initid) then
4733  field_name = trim(comp_name) // trim(stock_names(index))
4734  field_name = trim(field_name) // 'StocksChange_Diff'
4735  units = trim(stock_units(index))
4736  fmc_valuediagid(index,compind) = register_diag_field('stock_print', field_name, time, &
4737  units=units)
4738  endif
4739 
4740  diagid=f_valuediagid(index,compind)
4741  diagfield = f_value
4742  if (diagid > 0) used = send_data(diagid, diagfield, time)
4743  diagid=c_valuediagid(index,compind)
4744  diagfield = c_value
4745  if (diagid > 0) used = send_data(diagid, diagfield, time)
4746  diagid=fmc_valuediagid(index,compind)
4747  diagfield = f_value-c_value
4748  if (diagid > 0) used = send_data(diagid, diagfield, time)
4749 
4750 
4751  call get_time(time, isec, iday)
4752  hours = iday*24 + isec/3600
4753  formatstring = '(a,a,a,i16,2x,es22.15,2x,es22.15,2x,es22.15)'
4754  write(stocks_file,formatstring) trim(comp_name),stock_names(index),stock_units(index) &
4755  ,hours,f_value,c_value,f_value-c_value
4756 
4757  endif
4758 
4759 
4760 end subroutine stock_print
4761 
4762 
4763 !###############################################################################
4764  function is_lat_lon(lon, lat)
4765  real, dimension(:,:), intent(in) :: lon, lat
4766  logical :: is_lat_lon
4767  integer :: i, j, nlon, nlat, num
4768 
4769  is_lat_lon = .true.
4770  nlon = size(lon,1)
4771  nlat = size(lon,2)
4772  loop_lat: do j = 1, nlat
4773  do i = 2, nlon
4774  if(lat(i,j) .NE. lat(1,j)) then
4775  is_lat_lon = .false.
4776  exit loop_lat
4777  end if
4778  end do
4779  end do loop_lat
4780 
4781  if(is_lat_lon) then
4782  loop_lon: do i = 1, nlon
4783  do j = 2, nlat
4784  if(lon(i,j) .NE. lon(i,1)) then
4785  is_lat_lon = .false.
4786  exit loop_lon
4787  end if
4788  end do
4789  end do loop_lon
4790  end if
4791 
4792  num = 0
4793  if(is_lat_lon) num = 1
4794  call mpp_min(num)
4795  if(num == 1) then
4796  is_lat_lon = .true.
4797  else
4798  is_lat_lon = .false.
4799  end if
4800 
4801  return
4802  end function is_lat_lon
4803 
4804 !#######################################################################
4805 
4806 ! <SUBROUTINE NAME="get_side1_from_xgrid_ug" INTERFACE="get_from_xgrid_ug">
4807 ! <IN NAME="x" TYPE="real" DIM="(:)" > </IN>
4808 ! <IN NAME="grid_id" TYPE=" character(len=3)" > </IN>
4809 ! <OUT NAME="d" TYPE="real" DIM="(:,:)" > </OUT>
4810 ! <INOUT NAME="xmap" TYPE="xmap_type" > </INOUT>
4811 
4812 subroutine get_side1_from_xgrid_ug(d, grid_id, x, xmap, complete)
4813  real, dimension(:), intent( out) :: d
4814  character(len=3), intent(in ) :: grid_id
4815  real, dimension(:), intent(in ) :: x
4816  type(xmap_type), intent(inout) :: xmap
4817  logical, intent(in), optional :: complete
4818 
4819  logical :: is_complete, set_mismatch
4820  integer :: g
4821  character(len=2) :: text
4822  integer, save :: isize=0
4823  integer, save :: lsize=0
4824  integer, save :: xsize=0
4825  character(len=3), save :: grid_id_saved=""
4826  integer(LONG_KIND), dimension(MAX_FIELDS), save :: d_addrs=-9999
4827  integer(LONG_KIND), dimension(MAX_FIELDS), save :: x_addrs=-9999
4828 
4829  if (grid_id==xmap%grids(1)%id) then
4830  is_complete = .true.
4831  if(present(complete)) is_complete=complete
4832  lsize = lsize + 1
4833  if( lsize > max_fields ) then
4834  write( text,'(i2)' ) max_fields
4835  call error_mesg ('xgrid_mod', 'MAX_FIELDS='//trim(text)//' exceeded for group get_side1_from_xgrid_ug', fatal)
4836  endif
4837  d_addrs(lsize) = loc(d)
4838  x_addrs(lsize) = loc(x)
4839 
4840  if(lsize == 1) then
4841  isize = size(d(:))
4842  xsize = size(x(:))
4843  grid_id_saved = grid_id
4844  else
4845  set_mismatch = .false.
4846  set_mismatch = set_mismatch .OR. (isize /= size(d(:)))
4847  set_mismatch = set_mismatch .OR. (xsize /= size(x(:)))
4848  set_mismatch = set_mismatch .OR. (grid_id_saved /= grid_id)
4849  if(set_mismatch)then
4850  write( text,'(i2)' ) lsize
4851  call error_mesg ('xgrid_mod', 'Incompatible field at count '//text//' for group get_side1_from_xgrid_ug', fatal )
4852  endif
4853  endif
4854 
4855  if(is_complete) then
4856  if (make_exchange_reproduce) then
4857  call get_1_from_xgrid_ug_repro(d_addrs, x_addrs, xmap, xsize, lsize)
4858  else
4859  call get_1_from_xgrid_ug(d_addrs, x_addrs, xmap, isize, xsize, lsize)
4860  end if
4861  d_addrs(1:lsize) = -9999
4862  x_addrs(1:lsize) = -9999
4863  isize = 0
4864  xsize = 0
4865  lsize = 0
4866  grid_id_saved = ""
4867  endif
4868  return;
4869  end if
4870 
4871  do g=2,size(xmap%grids(:))
4872  if (grid_id==xmap%grids(g)%id) &
4873  call error_mesg ('xgrid_mod', &
4874  'get_from_xgrid_ug expects a 3D side 2 grid', fatal)
4875  end do
4876 
4877  call error_mesg ('xgrid_mod', 'get_from_xgrid_ug: could not find grid id', fatal)
4878 
4879 end subroutine get_side1_from_xgrid_ug
4880 ! </SUBROUTINE>
4881 
4882 !#######################################################################
4883 
4884 ! <SUBROUTINE NAME="put_side1_to_xgrid_ug" INTERFACE="put_to_xgrid_ug">
4885 ! <IN NAME="d" TYPE="real" DIM="(:,:)" > </IN>
4886 ! <IN NAME="grid_id" TYPE=" character(len=3)" > </IN>
4887 ! <INOUT NAME="x" TYPE="real" DIM="(:)" > </INOUT>
4888 ! <INOUT NAME="xmap" TYPE="xmap_type" > </INOUT>
4889 ! <IN NAME="remap_method" TYPE="integer,optional"></IN>
4890 ! Currently only support first order.
4891 subroutine put_side1_to_xgrid_ug(d, grid_id, x, xmap, complete)
4892  real, dimension(:), intent(in ) :: d
4893  character(len=3), intent(in ) :: grid_id
4894  real, dimension(:), intent(inout) :: x
4895  type(xmap_type), intent(inout) :: xmap
4896  logical, intent(in), optional :: complete
4897 
4898  logical :: is_complete, set_mismatch
4899  integer :: g
4900  character(len=2) :: text
4901  integer, save :: dsize=0
4902  integer, save :: lsize=0
4903  integer, save :: xsize=0
4904  character(len=3), save :: grid_id_saved=""
4905  integer(LONG_KIND), dimension(MAX_FIELDS), save :: d_addrs=-9999
4906  integer(LONG_KIND), dimension(MAX_FIELDS), save :: x_addrs=-9999
4907 
4908  if (grid_id==xmap%grids(1)%id) then
4909  is_complete = .true.
4910  if(present(complete)) is_complete=complete
4911  lsize = lsize + 1
4912  if( lsize > max_fields ) then
4913  write( text,'(i2)' ) max_fields
4914  call error_mesg ('xgrid_mod', 'MAX_FIELDS='//trim(text)//' exceeded for group put_side1_to_xgrid_ug', fatal)
4915  endif
4916  d_addrs(lsize) = loc(d)
4917  x_addrs(lsize) = loc(x)
4918 
4919  if(lsize == 1) then
4920  dsize = size(d(:))
4921  xsize = size(x(:))
4922  grid_id_saved = grid_id
4923  else
4924  set_mismatch = .false.
4925  set_mismatch = set_mismatch .OR. (dsize /= size(d(:)))
4926  set_mismatch = set_mismatch .OR. (xsize /= size(x(:)))
4927  set_mismatch = set_mismatch .OR. (grid_id_saved /= grid_id)
4928  if(set_mismatch)then
4929  write( text,'(i2)' ) lsize
4930  call error_mesg ('xgrid_mod', 'Incompatible field at count '//text//' for group put_side1_to_xgrid_ug', fatal )
4931  endif
4932  endif
4933 
4934  if(is_complete) then
4935  call put_1_to_xgrid_ug_order_1(d_addrs, x_addrs, xmap, dsize, xsize, lsize)
4936  d_addrs(1:lsize) = -9999
4937  x_addrs(1:lsize) = -9999
4938  dsize = 0
4939  xsize = 0
4940  lsize = 0
4941  grid_id_saved = ""
4942  endif
4943  return
4944  end if
4945 
4946  do g=2,size(xmap%grids(:))
4947  if (grid_id==xmap%grids(g)%id) &
4948  call error_mesg ('xgrid_mod', &
4949  'put_to_xgrid_ug expects a 2D side 2 grid', fatal)
4950  end do
4951 
4952  call error_mesg ('xgrid_mod', 'put_to_xgrid_ug: could not find grid id', fatal)
4953 
4954 end subroutine put_side1_to_xgrid_ug
4955 ! </SUBROUTINE>
4956 
4957 !#######################################################################
4958 
4959 ! <SUBROUTINE NAME="put_side2_to_xgrid_ug" INTERFACE="put_to_xgrid_ug">
4960 ! <IN NAME="d" TYPE="real" DIM="(:,:)" > </IN>
4961 ! <IN NAME="grid_id" TYPE=" character(len=3)" > </IN>
4962 ! <INOUT NAME="x" TYPE="real" DIM="(:)" > </INOUT>
4963 ! <INOUT NAME="xmap" TYPE="xmap_type" > </INOUT>
4964 
4965 subroutine put_side2_to_xgrid_ug(d, grid_id, x, xmap)
4966  real, dimension(:,:), intent(in ) :: d
4967  character(len=3), intent(in ) :: grid_id
4968  real, dimension(:), intent(inout) :: x
4969  type(xmap_type), intent(inout) :: xmap
4970 
4971  integer :: g
4972 
4973  if (grid_id==xmap%grids(1)%id) &
4974  call error_mesg ('xgrid_mod', &
4975  'put_to_xgrid_ug expects a 2D side 1 grid', fatal)
4976 
4977  do g=2,size(xmap%grids(:))
4978  if (grid_id==xmap%grids(g)%id) then
4979  call put_2_to_xgrid_ug(d, xmap%grids(g), x, xmap)
4980  return;
4981  end if
4982  end do
4983 
4984  call error_mesg ('xgrid_mod', 'put_to_xgrid_ug: could not find grid id', fatal)
4985 
4986 end subroutine put_side2_to_xgrid_ug
4987 ! </SUBROUTINE>
4988 
4989 !#######################################################################
4990 
4991 ! <SUBROUTINE NAME="get_side2_from_xgrid_ug" INTERFACE="get_from_xgrid_ug">
4992 ! <IN NAME="x" TYPE="real" DIM="(:)" > </IN>
4993 ! <IN NAME="grid_id" TYPE=" character(len=3)" > </IN>
4994 ! <OUT NAME="d" TYPE="real" DIM="(:,:)" > </OUT>
4995 ! <INOUT NAME="xmap" TYPE="xmap_type" > </INOUT>
4996 
4997 subroutine get_side2_from_xgrid_ug(d, grid_id, x, xmap)
4998  real, dimension(:,:), intent( out) :: d
4999  character(len=3), intent(in ) :: grid_id
5000  real, dimension(:), intent(in ) :: x
5001  type(xmap_type), intent(in ) :: xmap
5002 
5003  integer :: g
5004 
5005  if (grid_id==xmap%grids(1)%id) &
5006  call error_mesg ('xgrid_mod', &
5007  'get_from_xgrid_ug expects a 2D side 1 grid', fatal)
5008 
5009  do g=2,size(xmap%grids(:))
5010  if (grid_id==xmap%grids(g)%id) then
5011  call get_2_from_xgrid_ug(d, xmap%grids(g), x, xmap)
5012  return;
5013  end if
5014  end do
5015 
5016  call error_mesg ('xgrid_mod', 'get_from_xgrid_ug: could not find grid id', fatal)
5017 
5018 end subroutine get_side2_from_xgrid_ug
5019 ! </SUBROUTINE>
5020 
5021 
5022 !#######################################################################
5023 
5024 subroutine put_1_to_xgrid_ug_order_1(d_addrs, x_addrs, xmap, dsize, xsize, lsize)
5025  integer(LONG_KIND), dimension(:), intent(in) :: d_addrs
5026  integer(LONG_KIND), dimension(:), intent(in) :: x_addrs
5027  type(xmap_type), intent(inout) :: xmap
5028  integer, intent(in) :: dsize, xsize, lsize
5029 
5030  integer :: i, j, p, buffer_pos, msgsize
5031  integer :: from_pe, to_pe, pos, n, l, count
5032  integer :: ibegin, istart, iend, start_pos
5033  type(comm_type), pointer, save :: comm =>null()
5034  real :: recv_buffer(xmap%put1%recvsize*lsize)
5035  real :: send_buffer(xmap%put1%sendsize*lsize)
5036  real :: unpack_buffer(xmap%put1%recvsize)
5037 
5038  real, dimension(dsize) :: d
5039  real, dimension(xsize) :: x
5040  pointer(ptr_d, d)
5041  pointer(ptr_x, x)
5042  integer :: lll
5043 
5044  call mpp_clock_begin(id_put_1_to_xgrid_order_1)
5045 
5046  !--- pre-post receiving
5047  comm => xmap%put1
5048  do p = 1, comm%nrecv
5049  msgsize = comm%recv(p)%count*lsize
5050  from_pe = comm%recv(p)%pe
5051  buffer_pos = comm%recv(p)%buffer_pos*lsize
5052  call mpp_recv(recv_buffer(buffer_pos+1), glen=msgsize, from_pe = from_pe, block=.false., tag=comm_tag_7)
5053  enddo
5054 
5055  !--- send the data
5056  buffer_pos = 0
5057  do p = 1, comm%nsend
5058  msgsize = comm%send(p)%count*lsize
5059  to_pe = comm%send(p)%pe
5060  pos = buffer_pos
5061  do l = 1, lsize
5062  ptr_d = d_addrs(l)
5063  do n = 1, comm%send(p)%count
5064  pos = pos + 1
5065  lll = comm%send(p)%i(n)
5066  send_buffer(pos) = d(lll)
5067  enddo
5068  enddo
5069  call mpp_send(send_buffer(buffer_pos+1), plen=msgsize, to_pe = to_pe, tag=comm_tag_7 )
5070  buffer_pos = buffer_pos + msgsize
5071  enddo
5072 
5073  call mpp_sync_self(check=event_recv)
5074 
5075  !--- unpack the buffer
5076  if( lsize == 1) then
5077  ptr_x = x_addrs(1)
5078  do l=1,xmap%size_put1
5079  x(l) = recv_buffer(xmap%x1_put(l)%pos)
5080  end do
5081  else
5082  start_pos = 0
5083 !$OMP parallel do default(none) shared(lsize,x_addrs,comm,recv_buffer,xmap) &
5084 !$OMP private(ptr_x,count,ibegin,istart,iend,pos,unpack_buffer)
5085  do l = 1, lsize
5086  ptr_x = x_addrs(l)
5087  do p = 1, comm%nrecv
5088  count = comm%recv(p)%count
5089  ibegin = comm%recv(p)%buffer_pos*lsize + 1
5090  istart = ibegin + (l-1)*count
5091  iend = istart + count - 1
5092  pos = comm%recv(p)%buffer_pos
5093  do n = istart, iend
5094  pos = pos + 1
5095  unpack_buffer(pos) = recv_buffer(n)
5096  enddo
5097  enddo
5098  do i=1,xmap%size_put1
5099  x(i) = unpack_buffer(xmap%x1_put(i)%pos)
5100  end do
5101  enddo
5102  endif
5103 
5104  call mpp_sync_self()
5105 
5106  call mpp_clock_end(id_put_1_to_xgrid_order_1)
5107 
5108 end subroutine put_1_to_xgrid_ug_order_1
5109 
5110 !#######################################################################
5111 
5112 subroutine put_2_to_xgrid_ug(d, grid, x, xmap)
5113 type(grid_type), intent(in) :: grid
5114 real, dimension(grid%ls_me:grid%le_me, grid%km), intent(in) :: d
5115 real, dimension(: ), intent(inout) :: x
5116 type(xmap_type), intent(in ) :: xmap
5117 
5118  integer :: l
5119  call mpp_clock_begin(id_put_2_to_xgrid)
5120 
5121  do l=grid%first,grid%last
5122  x(l) = d(xmap%x2(l)%l,xmap%x2(l)%k)
5123  end do
5124 
5125  call mpp_clock_end(id_put_2_to_xgrid)
5126 end subroutine put_2_to_xgrid_ug
5127 
5128 
5129 subroutine get_1_from_xgrid_ug(d_addrs, x_addrs, xmap, isize, xsize, lsize)
5130  integer(LONG_KIND), dimension(:), intent(in) :: d_addrs
5131  integer(LONG_KIND), dimension(:), intent(in) :: x_addrs
5132  type(xmap_type), intent(inout) :: xmap
5133  integer, intent(in) :: isize, xsize, lsize
5134 
5135  real, dimension(xmap%size), target :: dg(xmap%size, lsize)
5136  integer :: i, j, l, p, n, m
5137  integer :: msgsize, buffer_pos, pos
5138  integer :: istart, iend, count
5139  real , pointer, save :: dgp =>null()
5140  type(grid_type) , pointer, save :: grid1 =>null()
5141  type(comm_type) , pointer, save :: comm =>null()
5142  type(overlap_type), pointer, save :: send => null()
5143  type(overlap_type), pointer, save :: recv => null()
5144  real :: recv_buffer(xmap%get1%recvsize*lsize*3)
5145  real :: send_buffer(xmap%get1%sendsize*lsize*3)
5146  real :: unpack_buffer(xmap%get1%recvsize*3)
5147  real :: d(isize)
5148  real, dimension(xsize) :: x
5149  pointer(ptr_d, d)
5150  pointer(ptr_x, x)
5151 
5152  call mpp_clock_begin(id_get_1_from_xgrid)
5153 
5154  comm => xmap%get1
5155  grid1 => xmap%grids(1)
5156 
5157  do p = 1, comm%nrecv
5158  recv => comm%recv(p)
5159  msgsize = recv%count*lsize
5160  buffer_pos = recv%buffer_pos*lsize
5161  call mpp_recv(recv_buffer(buffer_pos+1), glen=msgsize, from_pe = recv%pe, block=.false., tag=comm_tag_9)
5162  enddo
5163 
5164  dg = 0.0;
5165 !$OMP parallel do default(none) shared(lsize,xmap,dg,x_addrs) private(dgp,ptr_x)
5166  do l = 1, lsize
5167  ptr_x = x_addrs(l)
5168  do i=1,xmap%size
5169  dgp => dg(xmap%x1(i)%pos,l)
5170  dgp = dgp + xmap%x1(i)%area*x(i)
5171  enddo
5172  enddo
5173 
5174 
5175  !--- send the data
5176  buffer_pos = 0
5177  istart = 1
5178  do p = 1, comm%nsend
5179  send => comm%send(p)
5180  msgsize = send%count*lsize
5181  pos = buffer_pos
5182  istart = send%buffer_pos+1
5183  iend = istart + send%count - 1
5184  do l = 1, lsize
5185  do n = istart, iend
5186  pos = pos + 1
5187  send_buffer(pos) = dg(n,l)
5188  enddo
5189  enddo
5190  call mpp_send(send_buffer(buffer_pos+1), plen=msgsize, to_pe = send%pe, tag=comm_tag_9 )
5191  buffer_pos = buffer_pos + msgsize
5192  istart = iend + 1
5193  enddo
5194 
5195  call mpp_sync_self(check=event_recv)
5196 
5197  !--- unpack the buffer
5198  do l = 1, lsize
5199  ptr_d = d_addrs(l)
5200  d = 0.0
5201  enddo
5202  !--- To bitwise reproduce old results, first copy the data onto its own pe.
5203 
5204  do p = 1, comm%nrecv
5205  recv => comm%recv(p)
5206  count = recv%count
5207  buffer_pos = recv%buffer_pos*lsize
5208  if( recv%pe == xmap%me ) then
5209 !$OMP parallel do default(none) shared(lsize,recv,recv_buffer,buffer_pos,d_addrs,count) &
5210 !$OMP private(ptr_d,i,pos)
5211  do l = 1, lsize
5212  pos = buffer_pos + (l-1)*count
5213  ptr_d = d_addrs(l)
5214  do n = 1,count
5215  i = recv%i(n)
5216  pos = pos + 1
5217  d(i) = recv_buffer(pos)
5218  enddo
5219  enddo
5220  exit
5221  endif
5222  enddo
5223 
5224  pos = 0
5225  do m = 1, comm%nrecv
5226  p = comm%unpack_ind(m)
5227  recv => comm%recv(p)
5228  if( recv%pe == xmap%me ) then
5229  cycle
5230  endif
5231  buffer_pos = recv%buffer_pos*lsize
5232 !$OMP parallel do default(none) shared(lsize,recv,recv_buffer,buffer_pos,d_addrs) &
5233 !$OMP private(ptr_d,i,j,pos)
5234  do l = 1, lsize
5235  pos = buffer_pos + (l-1)*recv%count
5236  ptr_d = d_addrs(l)
5237  do n = 1, recv%count
5238  i = recv%i(n)
5239  pos = pos + 1
5240  d(i) = d(i) + recv_buffer(pos)
5241  enddo
5242  enddo
5243  enddo
5244 
5245  !
5246  ! normalize with side 1 grid cell areas
5247  !
5248 !$OMP parallel do default(none) shared(lsize,d_addrs,grid1) private(ptr_d)
5249  do l = 1, lsize
5250  ptr_d = d_addrs(l)
5251  d = d * grid1%area_inv(:,1)
5252  enddo
5253  call mpp_sync_self()
5254  call mpp_clock_end(id_get_1_from_xgrid)
5255 
5256 end subroutine get_1_from_xgrid_ug
5257 
5258 !#######################################################################
5259 
5260 subroutine get_1_from_xgrid_ug_repro(d_addrs, x_addrs, xmap, xsize, lsize)
5261  integer(LONG_KIND), dimension(:), intent(in) :: d_addrs
5262  integer(LONG_KIND), dimension(:), intent(in) :: x_addrs
5263  type(xmap_type), intent(inout) :: xmap
5264  integer, intent(in) :: xsize, lsize
5265 
5266  integer :: g, i, j, k, p, l, n, l2, m, l3
5267  integer :: msgsize, buffer_pos, pos
5268  type(grid_type), pointer, save :: grid =>null()
5269  type(comm_type), pointer, save :: comm => null()
5270  type(overlap_type), pointer, save :: send => null()
5271  type(overlap_type), pointer, save :: recv => null()
5272  integer, dimension(0:xmap%npes-1) :: pl, ml
5273  real :: recv_buffer(xmap%recv_count_repro_tot*lsize)
5274  real :: send_buffer(xmap%send_count_repro_tot*lsize)
5275  real :: d(xmap%grids(1)%ls_me:xmap%grids(1)%le_me)
5276  real, dimension(xsize) :: x
5277  pointer(ptr_d, d)
5278  pointer(ptr_x, x)
5279 
5280  call mpp_clock_begin(id_get_1_from_xgrid_repro)
5281  comm => xmap%get1_repro
5282  !--- pre-post receiving
5283  do p = 1, comm%nrecv
5284  recv => comm%recv(p)
5285  msgsize = recv%count*lsize
5286  buffer_pos = recv%buffer_pos*lsize
5287  call mpp_recv(recv_buffer(buffer_pos+1), glen=msgsize, from_pe = recv%pe, block=.false., tag=comm_tag_10)
5288  n = recv%pe -xmap%root_pe
5289  pl(n) = buffer_pos
5290  ml(n) = recv%count
5291  enddo
5292 
5293  !pack the data
5294  send_buffer(:) = 0.0
5295 !$OMP parallel do default(none) shared(lsize,x_addrs,comm,xmap,send_buffer) &
5296 !$OMP private(ptr_x,i,j,g,l2,pos,send)
5297  do p = 1, comm%nsend
5298  pos = comm%send(p)%buffer_pos*lsize
5299  send => comm%send(p)
5300  do l = 1,lsize
5301  ptr_x = x_addrs(l)
5302  do n = 1, send%count
5303  i = send%i(n)
5304  j = send%j(n)
5305  g = send%g(n)
5306  l2 = send%xloc(n)
5307  pos = pos + 1
5308  do k =1, xmap%grids(g)%km
5309  if(xmap%grids(g)%frac_area(i,j,k)/=0.0) then
5310  l2 = l2+1
5311  send_buffer(pos) = send_buffer(pos) + xmap%x1(l2)%area *x(l2)
5312  endif
5313  enddo
5314  enddo
5315  enddo
5316  enddo
5317 
5318  do p =1, comm%nsend
5319  buffer_pos = comm%send(p)%buffer_pos*lsize
5320  msgsize = comm%send(p)%count*lsize
5321  call mpp_send(send_buffer(buffer_pos+1), plen=msgsize, to_pe=comm%send(p)%pe, tag=comm_tag_10)
5322  enddo
5323 
5324  do l = 1, lsize
5325  ptr_d = d_addrs(l)
5326  d = 0
5327  enddo
5328 
5329  call mpp_sync_self(check=event_recv)
5330 
5331 !$OMP parallel do default(none) shared(lsize,d_addrs,xmap,recv_buffer,pl,ml) &
5332 !$OMP private(ptr_d,grid,i,j,p,pos)
5333  do l = 1, lsize
5334  ptr_d = d_addrs(l)
5335  do g=2,size(xmap%grids(:))
5336  grid => xmap%grids(g)
5337  do l3=1,grid%size_repro ! index into side1 grid's patterns
5338  i = grid%x_repro(l3)%l1
5339  p = grid%x_repro(l3)%pe-xmap%root_pe
5340  pos = pl(p) + (l-1)*ml(p) + grid%x_repro(l3)%recv_pos
5341  d(i) = d(i) + recv_buffer(pos)
5342  end do
5343  end do
5344  ! normalize with side 1 grid cell areas
5345  d = d * xmap%grids(1)%area_inv(:,1)
5346  enddo
5347 
5348  call mpp_sync_self()
5349 
5350  call mpp_clock_end(id_get_1_from_xgrid_repro)
5351 
5352 end subroutine get_1_from_xgrid_ug_repro
5353 
5354 !#######################################################################
5355 
5356 subroutine get_2_from_xgrid_ug(d, grid, x, xmap)
5357 type(grid_type), intent(in ) :: grid
5358 real, dimension(grid%ls_me:grid%le_me, grid%km), intent(out) :: d
5359 real, dimension(:), intent(in ) :: x
5360 type(xmap_type), intent(in ) :: xmap
5361 
5362  integer :: l, k
5363 
5364  call mpp_clock_begin(id_get_2_from_xgrid)
5365 
5366  d = 0.0
5367  do l=grid%first_get,grid%last_get
5368  d(xmap%x2_get(l)%l,xmap%x2_get(l)%k) = &
5369  d(xmap%x2_get(l)%l,xmap%x2_get(l)%k) + xmap%x2_get(l)%area*x(xmap%x2_get(l)%pos)
5370  end do
5371  !
5372  ! normalize with side 2 grid cell areas
5373  !
5374  do k=1,size(d,2)
5375  d(:,k) = d(:,k) * grid%area_inv(:,1)
5376  end do
5377 
5378  call mpp_clock_end(id_get_2_from_xgrid)
5379 
5380 end subroutine get_2_from_xgrid_ug
5381 
5382 !######################################################################
5383 logical function in_box_me(i, j, grid)
5384  integer, intent(in) :: i, j
5385  type(grid_type), intent(in) :: grid
5386  integer :: g
5387 
5388  if(grid%is_ug) then
5389  g = (j-1)*grid%ni + i
5390  in_box_me = (g>=grid%gs_me) .and. (g<=grid%ge_me)
5391  else
5392  in_box_me = (i>=grid%is_me) .and. (i<=grid%ie_me) .and. (j>=grid%js_me) .and. (j<=grid%je_me)
5393  endif
5394 
5395 end function in_box_me
5396 
5397 !######################################################################
5398 logical function in_box_nbr(i, j, grid, p)
5399  integer, intent(in) :: i, j, p
5400  type(grid_type), intent(in) :: grid
5401  integer :: g
5402 
5403  if(grid%is_ug) then
5404  g = (j-1)*grid%ni + i
5405  in_box_nbr = (g>=grid%gs(p)) .and. (g<=grid%ge(p))
5406  else
5407  in_box_nbr = (i>=grid%is(p)) .and. (i<=grid%ie(p)) .and. (j>=grid%js(p)) .and. (j<=grid%je(p))
5408  endif
5409 
5410 end function in_box_nbr
5411 
5412 
5413 end module xgrid_mod
5414 
5415 
5416 ! <INFO>
5417 
5418 ! <REFERENCE>
5419 ! A <LINK SRC="http://www.gfdl.noaa.gov/~mw/docs/grid_coupling.html"> guide </LINK>to grid coupling in FMS.
5420 ! </REFERENCE>
5421 ! <REFERENCE>
5422 ! A simple xgrid <LINK SRC="http://www.gfdl.gov/~mw/docs/xgrid_example.f90.txt"> example. </LINK>
5423 ! </REFERENCE>
5424 
5425 ! </INFO>
5426 
5427 
5428 !======================================================================================
5429 ! standalone unit test
5430 
5431 
5432 #ifdef _XGRID_MAIN
5433 ! to compile on Altix:
5434 ! setenv FMS /net2/ap/regression/ia64/10-Aug-2006/CM2.1U_Control-1990_E1.k_dyn30pe/exec
5435 ! ifort -fpp -r8 -i4 -g -check all -D_XGRID_MAIN -I $FMS xgrid.f90 $FMS/stock_constants.o $FMS/fms*.o $FMS/mpp*.o $FMS/constants.o $FMS/time_manager.o $FMS/memutils.o $FMS/threadloc.o -L/usr/local/lib -lnetcdf -L/usr/lib -lmpi -lsma
5436 ! mpirun -np 30 a.out
5437 program main
5438  use mpp_mod
5440  use mpp_domains_mod
5441  use xgrid_mod, only : xmap_type, setup_xmap, stock_move, stock_type, get_index_range
5443  use constants_mod, only : pi
5444  implicit none
5445 
5446  type(xmap_type) :: xmap_sfc, xmap_runoff
5447  integer :: npes, pe, root, i, nx, ny, ier
5448  integer :: patm_beg, patm_end, pocn_beg, pocn_end
5449  integer :: is, ie, js, je, km, index_ice, index_lnd
5450  integer :: layout(2)
5451  type(stock_type), save :: atm_stock(nelems), ice_stock(nelems), &
5452  & Lnd_stock(NELEMS), Ocn_stock(NELEMS)
5453  type(domain2d) :: atm_domain, ice_domain, lnd_domain, ocn_domain
5454  logical, pointer :: maskmap(:,:)
5455  real, allocatable :: data2d(:,:), data3d(:,:,:)
5456  real :: dt, dq_tot_atm, dq_tot_ice, dq_tot_lnd, dq_tot_ocn
5457 
5458 
5459 
5460  call fms_init
5461 
5462  npes = mpp_npes()
5463  pe = mpp_pe()
5464  root = mpp_root_pe()
5465  patm_beg = 0
5466  patm_end = npes/2 - 1
5467  pocn_beg = patm_end + 1
5468  pocn_end = npes - 1
5469 
5470  if(npes /= 30) call mpp_error(fatal,'must run unit test on 30 pes')
5471 
5472  call mpp_domains_init ! (MPP_DEBUG)
5473 
5474  call mpp_declare_pelist( (/ (i, i=patm_beg, patm_end) /), 'atm_lnd_ice pes' )
5475  call mpp_declare_pelist( (/ (i, i=pocn_beg, pocn_end) /), 'ocn pes' )
5476 
5477  index_ice = 2 ! 2nd exchange grid
5478  index_lnd = 3 ! 3rd exchange grid
5479 
5480  dt = 1.0
5481 
5482  if(pe < 15) then
5483 
5484  call mpp_set_current_pelist( (/ (i, i=patm_beg, patm_end) /) )
5485 
5486  ! Lnd
5487  nx = 144
5488  ny = 90
5489  layout = (/ 5, 3 /)
5490  call mpp_define_domains( (/1,nx, 1,ny/), layout, lnd_domain, &
5491  & xflags = cyclic_global_domain, name = 'LAND MODEL' )
5492 
5493  ! Atm
5494  nx = 144
5495  ny = 90
5496  layout = (/1, 15/)
5497  call mpp_define_domains( (/1,nx, 1,ny/), layout, atm_domain)
5498 
5499  ! Ice
5500  nx = 360
5501  ny = 200
5502  layout = (/15, 1/)
5503  call mpp_define_domains( (/1,nx, 1,ny/), layout, ice_domain, name='ice_nohalo' )
5504 
5505  ! Build exchange grid
5506  call setup_xmap(xmap_sfc, (/ 'ATM', 'OCN', 'LND' /), &
5507  & (/ atm_domain, ice_domain, lnd_domain /), &
5508  & "INPUT/grid_spec.nc")
5509 
5510 ! call setup_xmap(xmap_sfc, (/ 'LND', 'OCN' /), &
5511 ! & (/ Lnd_domain, Ice_domain /), &
5512 ! & "INPUT/grid_spec.nc")
5513 
5514 
5515  ! Atm -> Ice
5516 
5517  i = index_ice
5518 
5519  call get_index_range(xmap=xmap_sfc, grid_index=i, is=is, ie=ie, js=js, je=je, km=km)
5520 
5521  allocate(data3d(is:ie, js:je, km))
5522  data3d(:,:,1 ) = 1.0/(4.0*pi)
5523  data3d(:,:,2:km) = 0.0
5524  call stock_move(from=atm_stock(istock_water), to=ice_stock(istock_water), &
5525  & grid_index=i, data=data3d, xmap=xmap_sfc, &
5526  & delta_t=dt, from_side=istock_bottom, to_side=istock_top, radius=1.0, ier=ier)
5527  deallocate(data3d)
5528 
5529  ! Atm -> Lnd
5530 
5531  i = index_lnd
5532 
5533  call get_index_range(xmap=xmap_sfc, grid_index=i, is=is, ie=ie, js=js, je=je, km=km)
5534 
5535  allocate(data3d(is:ie, js:je, km))
5536  data3d(:,:,1 ) = 1.0/(4.0*pi)
5537  data3d(:,:,2:km) = 0.0
5538  call stock_move(from=atm_stock(istock_water), to=lnd_stock(istock_water), &
5539  & grid_index=i, data=data3d, xmap=xmap_sfc, &
5540  & delta_t=dt, from_side=istock_bottom, to_side=istock_top, radius=1.0, ier=ier)
5541  deallocate(data3d)
5542 
5543  else ! pes: 15...29
5544 
5545  call mpp_set_current_pelist( (/ (i, i=pocn_beg, pocn_end) /) )
5546 
5547  ! Ocn
5548  nx = 360
5549  ny = 200
5550  layout = (/ 5, 3 /)
5551  call mpp_define_domains( (/1,nx,1,ny/), layout, ocn_domain, name='ocean model')
5552 
5553  endif
5554 
5555  ! Ice -> Ocn (same grid different layout)
5556 
5557  i = index_ice
5558 
5559  if( pe < pocn_beg ) then
5560 
5561  call get_index_range(xmap=xmap_sfc, grid_index=i, is=is, ie=ie, js=js, je=je, km=km)
5562 
5563  allocate(data3d(is:ie, js:je, km))
5564  data3d(:,:,1 ) = 1.0/(4.0*pi)
5565  data3d(:,:,2:km) = 0.0
5566  else
5567  is = 0
5568  ie = 0
5569  js = 0
5570  je = 0
5571  km = 0
5572  allocate(data3d(is:ie, js:je, km))
5573  endif
5574 
5575  call stock_move(from=ice_stock(istock_water), to=ocn_stock(istock_water), &
5576  & grid_index=i, data=data3d(:,:,1), xmap=xmap_sfc, &
5577  & delta_t=dt, from_side=istock_bottom, to_side=istock_top, radius=1.0, ier=ier)
5578  deallocate(data3d)
5579 
5580  ! Sum across sides and PEs
5581 
5582  dq_tot_atm = sum(atm_stock(istock_water)%dq)
5583  call mpp_sum(dq_tot_atm)
5584 
5585  dq_tot_lnd = sum(lnd_stock(istock_water)%dq)
5586  call mpp_sum(dq_tot_lnd)
5587 
5588  dq_tot_ice = sum(ice_stock(istock_water)%dq)
5589  call mpp_sum(dq_tot_ice)
5590 
5591  dq_tot_ocn = sum(ocn_stock(istock_water)%dq)
5592  call mpp_sum(dq_tot_ocn)
5593 
5594  if(pe==root) then
5595  write(*,'(a,4f10.7,a,e10.2)') ' Total delta_q(water) Atm/Lnd/Ice/Ocn: ', &
5596  & dq_tot_atm, dq_tot_lnd, dq_tot_ice, dq_tot_ocn, &
5597  & ' residue: ', dq_tot_atm + dq_tot_lnd + dq_tot_ice + dq_tot_ocn
5598  endif
5599 
5600 
5601 end program main
5602 ! end of _XGRID_MAIN
5603 #endif
5604 
Definition: fms.F90:20
subroutine, public get_mosaic_contact(mosaic_file, tile1, tile2, istart1, iend1, jstart1, jend1, istart2, iend2, jstart2, jend2)
Definition: mosaic.F90:352
real, parameter, public radius
Radius of the Earth [m].
Definition: constants.F90:72
integer id_load_xgrid2
Definition: xgrid.F90:480
subroutine get_1_from_xgrid_ug_repro(d_addrs, x_addrs, xmap, xsize, lsize)
Definition: xgrid.F90:5263
subroutine put_2_to_xgrid_ug(d, grid, x, xmap)
Definition: xgrid.F90:5115
integer is_nest
Definition: xgrid.F90:487
subroutine load_xgrid(xmap, grid, grid_file, grid1_id, grid_id, tile1, tile2, use_higher_order)
Definition: xgrid.F90:612
integer id_put_1_to_xgrid_order_2
Definition: xgrid.F90:474
character(len=64) interp_method
Definition: xgrid.F90:207
real function, dimension(3) conservation_check_ug_side1(d, grid_id, xmap, remap_method)
Definition: xgrid.F90:4187
integer je_nest
Definition: xgrid.F90:487
subroutine put_1_to_xgrid_order_2(d_addrs, x_addrs, xmap, isize, jsize, xsize, lsize)
Definition: xgrid.F90:3593
logical make_exchange_reproduce
Definition: xgrid.F90:205
real, dimension(:,:), allocatable, public area_atm_model
Definition: xgrid.F90:221
logical monotonic_exchange
Definition: xgrid.F90:210
subroutine check(action, status)
integer je_parent
Definition: xgrid.F90:488
integer id_load_xgrid4
Definition: xgrid.F90:481
real(kind=kind_real), parameter g
gravity (m^2 s^{-2})
integer id_put_1_to_xgrid_order_1
Definition: xgrid.F90:473
integer id_load_xgrid3
Definition: xgrid.F90:480
subroutine get_2_from_xgrid_ug(d, grid, x, xmap)
Definition: xgrid.F90:5359
subroutine stock_move_2d(from, to, grid_index, data, xmap, delta_t, from_side, to_side, radius, verbose, ier)
Definition: xgrid.F90:4518
real, dimension(:,:), allocatable, public area_ocn_sphere
Definition: xgrid.F90:223
integer js_parent
Definition: xgrid.F90:488
integer, parameter, public istock_heat
integer, parameter, public istock_bottom
integer remapping_method
Definition: xgrid.F90:218
integer function, public get_mosaic_ntiles(mosaic_file)
Definition: mosaic.F90:214
integer tile_parent
Definition: xgrid.F90:486
logical xgrid_log
Definition: xgrid.F90:206
integer id_put_2_to_xgrid
Definition: xgrid.F90:478
subroutine get_area_elements(file, name, domain, data)
Definition: xgrid.F90:1432
logical do_alltoall
Definition: xgrid.F90:212
character(len=12), dimension(nelems), parameter stock_units
subroutine, public get_ocean_model_area_elements(domain, grid_file)
Definition: xgrid.F90:1468
integer id_get_1_from_xgrid_repro
Definition: xgrid.F90:476
Definition: mpp.F90:39
subroutine put_side2_to_xgrid_ug(d, grid_id, x, xmap)
Definition: xgrid.F90:4968
integer function, public check_nml_error(IOSTAT, NML_NAME)
Definition: fms.F90:658
integer, parameter, public first_order
Definition: xgrid.F90:179
integer, parameter, public istock_top
integer, parameter, public nelems
subroutine, public get_mosaic_grid_sizes(mosaic_file, nx, ny)
Definition: mosaic.F90:279
real, dimension(:,:), allocatable area_lnd_model
Definition: xgrid.F90:221
integer, parameter, public second_order
Definition: xgrid.F90:180
integer, parameter, public istock_side
subroutine get_side1_from_xgrid(d, grid_id, x, xmap, complete)
Definition: xgrid.F90:3295
character(len=input_str_length), dimension(:), allocatable, target, public input_nml_file
Definition: mpp.F90:1378
integer function get_nest_contact(mosaic_file, tile_nest_out, tile_parent_out, is_nest_out, ie_nest_out, js_nest_out, je_nest_out, is_parent_out, ie_parent_out, js_parent_out, je_parent_out)
Definition: xgrid.F90:2039
subroutine, public get_mosaic_xgrid(xgrid_file, i1, j1, i2, j2, area, ibegin, iend)
Definition: mosaic.F90:150
integer is_parent
Definition: xgrid.F90:488
integer id_load_xgrid5
Definition: xgrid.F90:481
subroutine set_frac_area_sg(f, grid_id, xmap)
Definition: xgrid.F90:3053
subroutine put_side2_to_xgrid(d, grid_id, x, xmap)
Definition: xgrid.F90:3263
real, dimension(:,:), allocatable, public area_ocn_model
Definition: xgrid.F90:221
subroutine regen(xmap)
Definition: xgrid.F90:2830
integer, public stocks_file
subroutine get_1_from_xgrid(d_addrs, x_addrs, xmap, isize, jsize, xsize, lsize)
Definition: xgrid.F90:3830
int field_exist(const char *file, const char *name)
Definition: read_mosaic.c:69
subroutine, public fms_init(localcomm)
Definition: fms.F90:353
real function, dimension(is:ie, js:je) grad_merid_latlon(d, lat, is, ie, js, je, isd, jsd)
Definition: xgrid.F90:4415
subroutine, public stock_print(stck, Time, comp_name, index, ref_value, radius, pelist)
Definition: xgrid.F90:4667
subroutine put_side1_to_xgrid_ug(d, grid_id, x, xmap, complete)
Definition: xgrid.F90:4894
logical do_alltoallv
Definition: xgrid.F90:213
integer id_load_xgrid1
Definition: xgrid.F90:480
subroutine stock_move_3d(from, to, grid_index, data, xmap, delta_t, from_side, to_side, radius, verbose, ier)
Definition: xgrid.F90:4461
subroutine get_side1_from_xgrid_ug(d, grid_id, x, xmap, complete)
Definition: xgrid.F90:4815
subroutine get_grid
logical function in_box(i, j, is, ie, js, je)
Definition: xgrid.F90:508
logical function in_box_me(i, j, grid)
Definition: xgrid.F90:5386
integer tile_nest
Definition: xgrid.F90:486
subroutine get_2_from_xgrid(d, grid, x, xmap)
Definition: xgrid.F90:3475
integer id_get_1_from_xgrid
Definition: xgrid.F90:475
logical function is_lat_lon(lon, lat)
Definition: xgrid.F90:4767
subroutine put_1_to_xgrid_ug_order_1(d_addrs, x_addrs, xmap, dsize, xsize, lsize)
Definition: xgrid.F90:5027
subroutine, public stock_integrate_2d(data, xmap, delta_t, radius, res, ier)
Definition: xgrid.F90:4637
subroutine, public get_index_range(xmap, grid_index, is, ie, js, je, km)
Definition: xgrid.F90:4445
logical xgrid_clocks_on
Definition: xgrid.F90:209
real function, dimension(is:ie, js:je) grad_zonal_latlon(d, lon, lat, is, ie, js, je, isd, jsd)
Definition: xgrid.F90:4375
logical debug_stocks
Definition: xgrid.F90:208
integer id_setup_xmap
Definition: xgrid.F90:479
integer, parameter version1
Definition: xgrid.F90:181
character(len=5), dimension(nelems), parameter stock_names
real function, dimension(3) conservation_check_side1(d, grid_id, xmap, remap_method)
Definition: xgrid.F90:4071
integer function, public get_mosaic_ncontacts(mosaic_file)
Definition: mosaic.F90:239
integer ie_parent
Definition: xgrid.F90:488
subroutine, public setup_xmap(xmap, grid_ids, grid_domains, grid_file, atm_grid, lnd_ug_domain)
Definition: xgrid.F90:1514
integer id_regen
Definition: xgrid.F90:482
integer, parameter, public istock_water
logical init
Definition: xgrid.F90:217
integer id_get_2_from_xgrid
Definition: xgrid.F90:477
subroutine get_side2_from_xgrid_ug(d, grid_id, x, xmap)
Definition: xgrid.F90:5000
************************************************************************GNU Lesser General Public License **This file is part of the GFDL Flexible Modeling System(FMS). ! *! *FMS is free software 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 see< http:! ***********************************************************************subroutine READ_RECORD_CORE_(unit, field, nwords, data, start, axsiz) integer, intent(in) ::unit type(fieldtype), intent(in) ::field integer, intent(in) ::nwords MPP_TYPE_, intent(inout) ::data(nwords) integer, intent(in) ::start(:), axsiz(:) integer(SHORT_KIND) ::i2vals(nwords)!rab used in conjunction with transfer intrinsic to determine size of a variable integer(KIND=1) ::one_byte(8) integer ::word_sz!#ifdef __sgi integer(INT_KIND) ::ivals(nwords) real(FLOAT_KIND) ::rvals(nwords)!#else! integer ::ivals(nwords)! real ::rvals(nwords)!#endif real(DOUBLE_KIND) ::r8vals(nwords) pointer(ptr1, i2vals) pointer(ptr2, ivals) pointer(ptr3, rvals) pointer(ptr4, r8vals) if(mpp_io_stack_size< nwords) call mpp_io_set_stack_size(nwords) call mpp_error(FATAL, 'MPP_READ currently requires use_netCDF option') end subroutine READ_RECORD_CORE_ subroutine READ_RECORD_(unit, field, nwords, data, time_level, domain, position, tile_count, start_in, axsiz_in)!routine that is finally called by all mpp_read routines to perform the read!a non-netCDF record contains:! field ID! a set of 4 coordinates(is:ie, js:je) giving the data subdomain! a timelevel and a timestamp(=NULLTIME if field is static)! 3D real data(stored as 1D)!if you are using direct access I/O, the RECL argument to OPEN must be large enough for the above!in a global direct access file, record position on PE is given by %record.!Treatment of timestamp:! We assume that static fields have been passed without a timestamp.! Here that is converted into a timestamp of NULLTIME.! For non-netCDF fields, field is treated no differently, but is written! with a timestamp of NULLTIME. There is no check in the code to prevent! the user from repeatedly writing a static field. integer, intent(in) ::unit, nwords type(fieldtype), intent(in) ::field MPP_TYPE_, intent(inout) ::data(nwords) integer, intent(in), optional ::time_level type(domain2D), intent(in), optional ::domain integer, intent(in), optional ::position, tile_count integer, intent(in), optional ::start_in(:), axsiz_in(:) integer, dimension(size(field%axes(:))) ::start, axsiz integer ::tlevel !, subdomain(4) integer ::i, error, is, ie, js, je, isg, ieg, jsg, jeg type(domain2d), pointer ::io_domain=> tlevel if(PRESENT(start_in) .AND. PRESENT(axsiz_in)) then if(size(start(! the data domain and compute domain must refer to the subdomain being passed ! In this ! since that attempts to gather all data on PE size(field%axes(:)) axsiz(i)
integer id_set_comm
Definition: xgrid.F90:482
integer id_conservation_check
Definition: xgrid.F90:482
real, parameter large_number
Definition: xgrid.F90:471
real, parameter eps
Definition: xgrid.F90:470
subroutine get_side2_from_xgrid(d, grid_id, x, xmap)
Definition: xgrid.F90:3377
logical module_is_initialized
Definition: xgrid.F90:472
subroutine get_1_from_xgrid_repro(d_addrs, x_addrs, xmap, xsize, lsize)
Definition: xgrid.F90:3963
subroutine put_side1_to_xgrid(d, grid_id, x, xmap, remap_method, complete)
Definition: xgrid.F90:3167
subroutine set_comm_put1(xmap)
Definition: xgrid.F90:2539
integer js_nest
Definition: xgrid.F90:487
subroutine put_1_to_xgrid_order_1(d_addrs, x_addrs, xmap, isize, jsize, xsize, lsize)
Definition: xgrid.F90:3504
integer function, public xgrid_count(xmap)
Definition: xgrid.F90:3151
subroutine, public gradient_cubic(pin, dx, dy, area, edge_w, edge_e, edge_s, edge_n, en_n, en_e, vlon, vlat, grad_x, grad_y, on_west_edge, on_east_edge, on_south_edge, on_north_edge)
Definition: gradient.F90:69
integer function, public get_mosaic_xgrid_size(xgrid_file)
Definition: mosaic.F90:114
logical function in_box_nbr(i, j, grid, p)
Definition: xgrid.F90:5401
real function, dimension(3) conservation_check_ug_side2(d, grid_id, xmap, remap_method)
Definition: xgrid.F90:4267
integer, parameter version2
Definition: xgrid.F90:182
program main
Definition: xgrid.F90:5439
integer ie_nest
Definition: xgrid.F90:487
real, dimension(:,:), allocatable area_lnd_sphere
Definition: xgrid.F90:223
subroutine, public get_time(Time, seconds, days, ticks, err_msg)
subroutine, public xgrid_init(remap_method)
Definition: xgrid.F90:532
subroutine, public error_mesg(routine, message, level)
Definition: fms.F90:529
subroutine, public some(xmap, some_arr, grid_id)
Definition: xgrid.F90:3421
subroutine, public get_xmap_grid_area(id, xmap, area)
Definition: xgrid.F90:4347
subroutine set_comm_get1_repro(xmap)
Definition: xgrid.F90:2135
integer, parameter max_fields
Definition: xgrid.F90:183
subroutine set_comm_get1(xmap)
Definition: xgrid.F90:2240
integer nsubset
Definition: xgrid.F90:211
real function, dimension(3) conservation_check_side2(d, grid_id, xmap, remap_method)
Definition: xgrid.F90:4125
integer nnest
Definition: xgrid.F90:486
subroutine get_1_from_xgrid_ug(d_addrs, x_addrs, xmap, isize, xsize, lsize)
Definition: xgrid.F90:5132
real, dimension(:,:), allocatable, public area_atm_sphere
Definition: xgrid.F90:223
integer id_load_xgrid
Definition: xgrid.F90:482
subroutine stock_move_ug_3d(from, to, grid_index, data, xmap, delta_t, from_side, to_side, radius, verbose, ier)
Definition: xgrid.F90:4578
subroutine put_2_to_xgrid(d, grid, x, xmap)
Definition: xgrid.F90:3456
subroutine, public set_frac_area_ug(f, grid_id, xmap)
Definition: xgrid.F90:3101
real(fp), parameter, public pi