FV3DYCORE  Version 2.0.0
fv_grid_tools.F90
Go to the documentation of this file.
1 !!***********************************************************************
2 !* GNU Lesser General Public License
3 !*
4 !* This file is part of the FV3 dynamical core.
5 !*
6 !* The FV3 dynamical core is free software: you can redistribute it
7 !* and/or modify it under the terms of the
8 !* GNU Lesser General Public License as published by the
9 !* Free Software Foundation, either version 3 of the License, or
10 !* (at your option) any later version.
11 !*
12 !* The FV3 dynamical core is distributed in the hope that it will be
13 !* useful, but WITHOUT ANYWARRANTY; without even the implied warranty
14 !* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
15 !* See the GNU General Public License for more details.
16 !*
17 !* You should have received a copy of the GNU Lesser General Public
18 !* License along with the FV3 dynamical core.
19 !* If not, see <http://www.gnu.org/licenses/>.
20 !***********************************************************************
21 module fv_grid_tools_mod
22 
23 ! <table>
24 ! <tr>
25 ! <th>Module Name</th>
26 ! <th>Functions Included</th>
27 ! </tr>
28 ! <tr>
29 ! <td>constants_mod</td>
30 ! <td>grav, omega, pi=>pi_8, cnst_radius=>radius</td>
31 ! </tr>
32 ! <tr>
33 ! <td>fms_mod</td>
34 ! <td>get_mosaic_tile_grid</td>
35 ! </tr>
36 ! <tr>
37 ! <td>fms_io_mod</td>
38 ! <td>file_exist, field_exist, read_data, get_global_att_value, get_var_att_value</td>
39 ! </tr>
40 ! <tr>
41 ! <td>fv_arrays_mod</td>
42 ! <td>fv_atmos_type, fv_grid_type, fv_grid_bounds_type, R_GRID</td>
43 ! </tr>
44 ! <tr>
45 ! <td>fv_control_mod</td>
46 ! <td>fv_init, fv_end, ngrids</td>
47 ! </tr>
48 ! <tr>
49 ! <td>fv_diagnostics_mod</td>
50 ! <td>prt_maxmin, prt_gb_nh_sh, prt_height</td>
51 ! </tr>
52 ! <tr>
53 ! <td>fv_eta_mod</td>
54 ! <td>set_eta, set_external_eta</td>
55 ! </tr>
56 ! <tr>
57 ! <td>fv_fill_mod</td>
58 ! <td>fillz</td>
59 ! </tr>
60 ! <tr>
61 ! <td>fv_grid_utils_mod</td>
62 ! <td>gnomonic_grids, great_circle_dist,mid_pt_sphere, spherical_angle,
63 ! cell_center2, get_area, inner_prod, fill_ghost, direct_transform,
64 ! dist2side_latlon,spherical_linear_interpolation, big_number</td>
65 ! </tr>
66 ! <tr>
67 ! <td>fv_io_mod</td>
68 ! <td>fv_io_read_tracers</td>
69 ! </tr>
70 ! <tr>
71 ! <td>fv_mp_mod</td>
72 ! <td>ng, is_master, fill_corners, XDir, YDir,mp_gather,
73 ! mp_bcst, mp_reduce_max, mp_stop </td>
74 ! </tr>
75 ! <tr>
76 ! <td>fv_timing_mod</td>
77 ! <td>timing_on, timing_off</td>
78 ! </tr>
79 ! <tr>
80 ! <td>mosaic_mod</td>
81 ! <td>get_mosaic_ntiles</td>
82 ! </tr>
83 ! <tr>
84 ! <td>mpp_mod</td>
85 ! <td>mpp_error, FATAL, get_unit, mpp_chksum, mpp_pe, stdout,
86 ! mpp_send, mpp_recv, mpp_sync_self, EVENT_RECV, mpp_npes,
87 ! mpp_sum, mpp_max, mpp_min, mpp_root_pe, mpp_broadcast, mpp_transmit</td>
88 ! </tr>
89 ! <tr>
90 ! <td>mpp_domains_mod</td>
91 ! <td>domain2d,mpp_update_domains, mpp_get_boundary,mpp_get_ntile_count,
92 ! mpp_get_pelist, mpp_get_compute_domains, mpp_global_field,
93 ! mpp_get_data_domain, mpp_get_compute_domain,mpp_get_global_domain,
94 ! mpp_global_sum, mpp_global_max, mpp_global_min</td>
95 ! </tr>
96 ! <tr>
97 ! <td>mpp_io_mod</td>
98 ! <td>mpp_get_att_value</td>
99 ! </tr>
100 ! <tr>
101 ! <td>mpp_parameter_mod</td>
102 ! <td>AGRID_PARAM=>AGRID,DGRID_NE_PARAM=>DGRID_NE,
103 ! CGRID_NE_PARAM=>CGRID_NE,CGRID_SW_PARAM=>CGRID_SW,
104 ! BGRID_NE_PARAM=>BGRID_NE,BGRID_SW_PARAM=>BGRID_SW,
105 ! SCALAR_PAIR,CORNER, CENTER, XUPDATE</td>
106 ! </tr>
107 ! <tr>
108 ! <td>sorted_index_mod</td>
109 ! <td>sorted_inta, sorted_intb</td>
110 ! </tr>
111 ! <tr>
112 ! <td>tracer_manager_mod</td>
113 ! <td>get_tracer_names, get_number_tracers, get_tracer_index, set_tracer_profile</td>
114 ! </tr>
115 ! </table>
116 
117 
118  use constants_mod, only: grav, omega, pi=>pi_8, cnst_radius=>radius
119  use fv_arrays_mod, only: fv_atmos_type, fv_grid_type, fv_grid_bounds_type, r_grid
120  use fv_grid_utils_mod, only: gnomonic_grids, great_circle_dist, &
121  mid_pt_sphere, spherical_angle, &
122  cell_center2, get_area, inner_prod, fill_ghost, &
123  direct_transform, cube_transform, dist2side_latlon, &
124  spherical_linear_interpolation, big_number
125  use fv_timing_mod, only: timing_on, timing_off
126  use fv_mp_mod, only: is_master, fill_corners, xdir, ydir
127  use fv_mp_mod, only: mp_gather, mp_bcst, mp_reduce_max, mp_stop, grids_master_procs
128  use sorted_index_mod, only: sorted_inta, sorted_intb
129  use mpp_mod, only: mpp_error, fatal, get_unit, mpp_chksum, mpp_pe, stdout, &
130  mpp_send, mpp_recv, mpp_sync_self, event_recv, mpp_npes, &
131  mpp_sum, mpp_max, mpp_min, mpp_root_pe, mpp_broadcast, mpp_transmit
132  use mpp_domains_mod, only: mpp_update_domains, mpp_get_boundary, &
133  mpp_get_ntile_count, mpp_get_pelist, &
134  mpp_get_compute_domains, mpp_global_field, &
135  mpp_get_data_domain, mpp_get_compute_domain, &
136  mpp_get_global_domain, mpp_global_sum, mpp_global_max, mpp_global_min
137  use mpp_domains_mod, only: domain2d
138  use mpp_io_mod, only: mpp_get_att_value
139 
140  use mpp_parameter_mod, only: agrid_param=>agrid, &
141  dgrid_ne_param=>dgrid_ne, &
142  cgrid_ne_param=>cgrid_ne, &
143  cgrid_sw_param=>cgrid_sw, &
144  bgrid_ne_param=>bgrid_ne, &
145  bgrid_sw_param=>bgrid_sw, &
146  scalar_pair, &
147  corner, center, xupdate
148  use fms_mod, only: get_mosaic_tile_grid
149  use fms_io_mod, only: file_exist, field_exist, read_data, &
150  get_global_att_value, get_var_att_value
151  use mosaic_mod, only : get_mosaic_ntiles
152 
153  use mpp_mod, only: mpp_transmit, mpp_recv
154  implicit none
155  private
156 #include <netcdf.inc>
157 
158  real(kind=R_GRID), parameter:: radius = cnst_radius
159 
160  real(kind=R_GRID) , parameter:: todeg = 180.0d0/pi
161  real(kind=R_GRID) , parameter:: torad = pi/180.0d0
162  real(kind=R_GRID) , parameter:: missing = 1.d25
163 
164  real(kind=R_GRID) :: csfac
165 
166  logical, parameter :: debug_message_size = .false.
167  logical :: write_grid_char_file = .false.
168 
169 
170  public :: todeg, missing, init_grid, spherical_to_cartesian
171 
172 contains
173 
175  subroutine read_grid(Atm, grid_file, ndims, nregions, ng)
176  type(fv_atmos_type), intent(inout), target :: Atm
177  character(len=*), intent(IN) :: grid_file
178  integer, intent(IN) :: ndims
179  integer, intent(IN) :: nregions
180  integer, intent(IN) :: ng
181 
182  real, allocatable, dimension(:,:) :: tmpx, tmpy
183  real(kind=R_GRID), pointer, dimension(:,:,:) :: grid
184  character(len=128) :: units = ""
185  character(len=256) :: atm_mosaic, atm_hgrid, grid_form
186  character(len=1024) :: attvalue
187  integer :: ntiles, i, j, stdunit
188  integer :: isc2, iec2, jsc2, jec2
189  integer :: start(4), nread(4)
190  integer :: is, ie, js, je
191  integer :: isd, ied, jsd, jed
192  integer,save :: halo=3 ! for regional domain external tools
193 
194  is = atm%bd%is
195  ie = atm%bd%ie
196  js = atm%bd%js
197  je = atm%bd%je
198  isd = atm%bd%isd
199  ied = atm%bd%ied
200  jsd = atm%bd%jsd
201  jed = atm%bd%jed
202  grid => atm%gridstruct%grid_64
203 
204  if(.not. file_exist(grid_file)) call mpp_error(fatal, 'fv_grid_tools(read_grid): file '// &
205  trim(grid_file)//' does not exist')
206 
207  !--- make sure the grid file is mosaic file.
208  if( field_exist(grid_file, 'atm_mosaic_file') .OR. field_exist(grid_file, 'gridfiles') ) then
209  stdunit = stdout()
210  write(stdunit,*) '==>Note from fv_grid_tools_mod(read_grid): read atmosphere grid from mosaic version grid'
211  else
212  call mpp_error(fatal, 'fv_grid_tools(read_grid): neither atm_mosaic_file nor gridfiles exists in file ' &
213  //trim(grid_file))
214  endif
215 
216  if(field_exist(grid_file, 'atm_mosaic_file')) then
217  call read_data(grid_file, "atm_mosaic_file", atm_mosaic)
218  atm_mosaic = "INPUT/"//trim(atm_mosaic)
219  else
220  atm_mosaic = trim(grid_file)
221  endif
222 
223  call get_mosaic_tile_grid(atm_hgrid, atm_mosaic, atm%domain)
224 
225  grid_form = "none"
226  if( get_global_att_value(atm_hgrid, "history", attvalue) ) then
227  if( index(attvalue, "gnomonic_ed") > 0) grid_form = "gnomonic_ed"
228  endif
229  if(grid_form .NE. "gnomonic_ed") call mpp_error(fatal, &
230  "fv_grid_tools(read_grid): the grid should be 'gnomonic_ed' when reading from grid file, contact developer")
231 
232  !FIXME: Doesn't work for a nested grid
233  ntiles = get_mosaic_ntiles(atm_mosaic)
234  if( .not. atm%gridstruct%bounded_domain) then
235  if(ntiles .NE. 6) call mpp_error(fatal, &
236  'fv_grid_tools(read_grid): ntiles should be 6 in mosaic file '//trim(atm_mosaic) )
237  if(nregions .NE. 6) call mpp_error(fatal, &
238  'fv_grid_tools(read_grid): nregions should be 6 when reading from mosaic file '//trim(grid_file) )
239  endif
240 
241  call get_var_att_value(atm_hgrid, 'x', 'units', units)
242 
243  !--- get the geographical coordinates of super-grid.
244  isc2 = 2*is-1; iec2 = 2*ie+1
245  jsc2 = 2*js-1; jec2 = 2*je+1
246  if( atm%gridstruct%bounded_domain ) then
247  isc2 = 2*(isd+halo)-1; iec2 = 2*(ied+1+halo)-1 ! For the regional domain the cell corner locations must be transferred
248  jsc2 = 2*(jsd+halo)-1; jec2 = 2*(jed+1+halo)-1 ! from the entire supergrid to the compute grid, including the halo region.
249  endif
250  allocate(tmpx(isc2:iec2, jsc2:jec2) )
251  allocate(tmpy(isc2:iec2, jsc2:jec2) )
252  start = 1; nread = 1
253  start(1) = isc2; nread(1) = iec2 - isc2 + 1
254  start(2) = jsc2; nread(2) = jec2 - jsc2 + 1
255  call read_data(atm_hgrid, 'x', tmpx, start, nread, no_domain=.true.)
256  call read_data(atm_hgrid, 'y', tmpy, start, nread, no_domain=.true.)
257 
258  !--- geographic grid at cell corner
259  grid(isd: is-1, jsd:js-1,1:ndims)=0.
260  grid(isd: is-1, je+2:jed+1,1:ndims)=0.
261  grid(ie+2:ied+1,jsd:js-1,1:ndims)=0.
262  grid(ie+2:ied+1,je+2:jed+1,1:ndims)=0.
263  if(len_trim(units) < 6) call mpp_error(fatal, &
264  "fv_grid_tools_mod(read_grid): the length of units must be no less than 6")
265  if(units(1:6) == 'degree') then
266  if( .not. atm%gridstruct%bounded_domain) then
267  do j = js, je+1
268  do i = is, ie+1
269  grid(i,j,1) = tmpx(2*i-1,2*j-1)*pi/180.
270  grid(i,j,2) = tmpy(2*i-1,2*j-1)*pi/180.
271  enddo
272  enddo
273  else
274 !
275 !*** In the regional case the halo surrounding the domain was included in the read.
276 !*** Transfer the compute and halo regions to the compute grid.
277 !
278  do j = jsd, jed+1
279  do i = isd, ied+1
280  grid(i,j,1) = tmpx(2*i+halo+2,2*j+halo+2)*pi/180.
281  grid(i,j,2) = tmpy(2*i+halo+2,2*j+halo+2)*pi/180.
282  enddo
283  enddo
284  endif
285 
286  else if(units(1:6) == 'radian') then
287  do j = js, je+1
288  do i = is, ie+1
289  grid(i,j,1) = tmpx(2*i-1,2*j-1)
290  grid(i,j,2) = tmpy(2*i-1,2*j-1)
291  enddo
292  enddo
293  else
294  print*, 'units is ' , trim(units), len_trim(units), mpp_pe()
295  call mpp_error(fatal, 'fv_grid_tools_mod(read_grid): units must start with degree or radian')
296  endif
297 
298  deallocate(tmpx, tmpy)
299  nullify(grid)
300 
301  end subroutine read_grid
302 
303 
304 
305  !#################################################################################
306  subroutine get_symmetry(data_in, data_out, ishift, jshift, npes_x, npes_y, domain, tile, npx_g, bd)
307  type(fv_grid_bounds_type), intent(IN) :: bd
308  integer, intent(in) :: ishift, jshift, npes_x, npes_y
309  real(kind=R_GRID), dimension(bd%is:bd%ie+ishift, bd%js:bd%je+jshift ), intent(in) :: data_in
310  real(kind=R_GRID), dimension(bd%is:bd%ie+jshift, bd%js:bd%je+ishift ), intent(out) :: data_out
311  real(kind=R_GRID), dimension(:), allocatable :: send_buffer
312  real(kind=R_GRID), dimension(:), allocatable :: recv_buffer
313  integer, dimension(:), allocatable :: is_recv, ie_recv, js_recv, je_recv, pe_recv
314  integer, dimension(:), allocatable :: is_send, ie_send, js_send, je_send, pe_send
315  integer, dimension(:), allocatable :: isl, iel, jsl, jel, pelist, msg1, msg2
316  integer :: msgsize, pos, ntiles, npes_per_tile, npes
317  integer :: send_buf_size, recv_buf_size, buffer_pos
318  integer :: is0, ie0, js0, je0
319  integer :: is1, ie1, js1, je1
320  integer :: is2, ie2, js2, je2
321  integer :: i, j, p, nrecv, nsend, tile_you, is3, ie3, nlist
322  integer :: start_pe, ipos, jpos, from_pe, to_pe
323  type(domain2d) :: domain
324  integer :: tile, npx_g
325  integer :: is, ie, js, je
326  integer :: isd, ied, jsd, jed
327 
328  is = bd%is
329  ie = bd%ie
330  js = bd%js
331  je = bd%je
332  isd = bd%isd
333  ied = bd%ied
334  jsd = bd%jsd
335  jed = bd%jed
336 
337  !--- This routine will be called only for cubic sphere grid. so 6 tiles will be assumed
338  !--- also number of processors on each tile will be the same.
339  ntiles = mpp_get_ntile_count(domain)
340  npes = mpp_npes()
341 
342  if(ntiles .NE. 6 ) call mpp_error(fatal, 'fv_grid_tools(get_symmetry): ntiles should be 6 ')
343  if(mod(npes,ntiles) /= 0) call mpp_error(fatal, 'fv_grid_tools(get_symmetry): npes should be divided by ntiles')
344  npes_per_tile = npes/ntiles
345 
346 ! if(npes_x == npes_y) then ! even, simple communication
347  if(npes_x == npes_y .AND. mod(npx_g-1,npes_x) == 0 ) then ! even,
348  msgsize = (ie-is+1+jshift)*(je-js+1+ishift)
349 
350  pos = mod((mpp_pe()-mpp_root_pe()), npes_x*npes_y)
351  start_pe = mpp_pe() - pos
352  ipos = mod(pos, npes_x)
353  jpos = pos/npes_x
354  from_pe = start_pe + ipos*npes_x + jpos
355  to_pe = from_pe
356  allocate(recv_buffer(msgsize))
357  call mpp_recv(recv_buffer(1), glen=msgsize, from_pe=from_pe, block=.false. )
358 
359  pos = 0
360  allocate(send_buffer(msgsize))
361  do j = js, je+jshift
362  do i = is, ie+ishift
363  pos = pos + 1
364  send_buffer(pos) = data_in(i,j)
365  enddo
366  enddo
367 
368  call mpp_send(send_buffer(1), plen=msgsize, to_pe=to_pe)
369  call mpp_sync_self(check=event_recv) ! To ensure recv is completed.
370 
371  !--unpack buffer
372  pos = 0
373  do i = is, ie+jshift
374  do j = js, je+ishift
375  pos = pos + 1
376  data_out(i,j) = recv_buffer(pos)
377  enddo
378  enddo
379 
380  call mpp_sync_self()
381  deallocate(send_buffer, recv_buffer)
382  else
383 
384  allocate(is_recv(0:npes_per_tile-1), ie_recv(0:npes_per_tile-1))
385  allocate(js_recv(0:npes_per_tile-1), je_recv(0:npes_per_tile-1))
386  allocate(is_send(0:npes_per_tile-1), ie_send(0:npes_per_tile-1))
387  allocate(js_send(0:npes_per_tile-1), je_send(0:npes_per_tile-1))
388  allocate(pe_send(0:npes_per_tile-1), pe_recv(0:npes_per_tile-1))
389  if(debug_message_size) then
390  allocate(msg1(0:npes_per_tile-1), msg2(0:npes_per_tile-1))
391  msg1 = 0
392  msg2 = 0
393  endif
394 
395  allocate(pelist(0:npes-1))
396  call mpp_get_pelist(domain, pelist)
397  allocate(isl(0:npes-1), iel(0:npes-1), jsl(0:npes-1), jel(0:npes-1) )
398  call mpp_get_compute_domains(domain, xbegin=isl, xend=iel, ybegin=jsl, yend=jel)
399  !--- pre-post receiving
400  buffer_pos = 0
401  nrecv = 0
402  nsend = 0
403  recv_buf_size = 0
404 
405  !--- first set up the receiving index
406  nlist = 0
407  do p = 0, npes-1
408  tile_you = p/(npes_x*npes_y) + 1
409  if(tile_you .NE. tile) cycle
410 
411  !--- my index for data_out after rotation
412  is1 = js; ie1 = je + ishift;
413  js1 = is; je1 = ie + jshift;
414  !--- your index for data_out
415  is2 = isl(p); ie2 = iel(p) + ishift;
416  js2 = jsl(p); je2 = jel(p) + jshift;
417  is0 = max(is1,is2); ie0 = min(ie1,ie2)
418  js0 = max(js1,js2); je0 = min(je1,je2)
419  msgsize = 0
420  if(ie0 .GE. is0 .AND. je0 .GE. js0) then
421  msgsize = (ie0-is0+1)*(je0-js0+1)
422  recv_buf_size = recv_buf_size + msgsize
423  pe_recv(nrecv) = pelist(p)
424  !--- need to rotate back the index
425  is_recv(nrecv) = js0; ie_recv(nrecv) = je0
426  js_recv(nrecv) = is0; je_recv(nrecv) = ie0
427  nrecv = nrecv+1
428  endif
429  if(debug_message_size) then
430  msg1(nlist) = msgsize
431  call mpp_recv(msg2(nlist), glen=1, from_pe=pelist(p), block=.false. )
432  nlist = nlist + 1
433  endif
434  enddo
435 
436  !--- Then setup the sending index.
437  send_buf_size = 0
438  do p = 0, npes-1
439  tile_you = p/(npes_x*npes_y) + 1
440  if(tile_you .NE. tile) cycle
441  !--- my index on data_in
442  is1 = is; ie1 = ie + ishift;
443  js1 = js; je1 = je + jshift;
444  !--- your index on data_out after rotate
445  is2 = jsl(p); ie2 = jel(p) + ishift;
446  js2 = isl(p); je2 = iel(p) + jshift;
447  is0 = max(is1,is2); ie0 = min(ie1,ie2)
448  js0 = max(js1,js2); je0 = min(je1,je2)
449  msgsize = 0
450  if(ie0 .GE. is0 .AND. je0 .GE. js0 )then
451  msgsize = (ie0-is0+1)*(je0-js0+1)
452  send_buf_size = send_buf_size + msgsize
453  pe_send(nsend) = pelist(p)
454  is_send(nsend) = is0; ie_send(nsend) = ie0
455  js_send(nsend) = js0; je_send(nsend) = je0
456  nsend = nsend+1
457  endif
458  IF(debug_message_size) call mpp_send(msgsize, plen=1, to_pe=pelist(p) )
459  enddo
460 
461  !--- check to make sure send and recv size match.
462  if(debug_message_size) then
463  call mpp_sync_self(check=event_recv) ! To ensure recv is completed.
464  do p = 0, nlist-1
465  if(msg1(p) .NE. msg2(p)) then
466  call mpp_error(fatal, "fv_grid_tools_mod(get_symmetry): mismatch on send and recv size")
467  endif
468  enddo
469  call mpp_sync_self()
470  deallocate(msg1, msg2)
471  endif
472 
473  !--- pre-post data
474  allocate(recv_buffer(recv_buf_size))
475  buffer_pos = 0
476  do p = 0, nrecv-1
477  is0 = is_recv(p); ie0 = ie_recv(p)
478  js0 = js_recv(p); je0 = je_recv(p)
479  msgsize = (ie0-is0+1)*(je0-js0+1)
480  call mpp_recv(recv_buffer(buffer_pos+1), glen=msgsize, from_pe=pe_recv(p), block=.false. )
481  buffer_pos = buffer_pos + msgsize
482  enddo
483 
484  !--- send the data
485  buffer_pos = 0
486  allocate(send_buffer(send_buf_size))
487  do p = 0, nsend-1
488  is0 = is_send(p); ie0 = ie_send(p)
489  js0 = js_send(p); je0 = je_send(p)
490  msgsize = (ie0-is0+1)*(je0-js0+1)
491  pos = buffer_pos
492  do j = js0, je0
493  do i = is0, ie0
494  pos = pos+1
495  send_buffer(pos) = data_in(i,j)
496  enddo
497  enddo
498  call mpp_send(send_buffer(buffer_pos+1), plen=msgsize, to_pe=pe_send(p) )
499  buffer_pos = buffer_pos + msgsize
500  enddo
501 
502  call mpp_sync_self(check=event_recv) ! To ensure recv is completed.
503 
504  !--- unpack buffer
505  pos = 0
506  do p = 0, nrecv-1
507  is0 = is_recv(p); ie0 = ie_recv(p)
508  js0 = js_recv(p); je0 = je_recv(p)
509 
510  do i = is0, ie0
511  do j = js0, je0
512  pos = pos + 1
513  data_out(i,j) = recv_buffer(pos)
514  enddo
515  enddo
516  enddo
517 
518  call mpp_sync_self()
519  deallocate(isl, iel, jsl, jel, pelist)
520  deallocate(is_recv, ie_recv, js_recv, je_recv, pe_recv)
521  deallocate(is_send, ie_send, js_send, je_send, pe_send)
522  deallocate(recv_buffer, send_buffer)
523  endif
524 
525  end subroutine get_symmetry
526 
529  subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, ng, tile_coarse)
530 !--------------------------------------------------------
531  type(fv_atmos_type), intent(inout), target :: Atm
532  character(len=80), intent(IN) :: grid_name
533  character(len=120),intent(IN) :: grid_file
534  integer, intent(IN) :: npx, npy, npz
535  integer, intent(IN) :: ndims
536  integer, intent(IN) :: nregions
537  integer, intent(IN) :: ng
538  integer, intent(IN) :: tile_coarse(:)
539 !--------------------------------------------------------
540  real(kind=R_GRID) :: xs(npx,npy)
541  real(kind=R_GRID) :: ys(npx,npy)
542 
543  real(kind=R_GRID) :: dp, dl
544  real(kind=R_GRID) :: x1,x2,y1,y2,z1,z2
545  integer :: i,j,k,n,nreg
546  integer :: fileLun
547 
548  real(kind=R_GRID) :: p1(3), p2(3), p3(3), p4(3)
549  real(kind=R_GRID) :: dist,dist1,dist2, pa(2), pa1(2), pa2(2), pb(2)
550  real(kind=R_GRID) :: pt(3), pt1(3), pt2(3), pt3(3)
551  real(kind=R_GRID) :: ee1(3), ee2(3)
552 
553  real(kind=R_GRID) :: angN,angM,angAV,ang
554  real(kind=R_GRID) :: aspN,aspM,aspAV,asp
555  real(kind=R_GRID) :: dxN, dxM, dxAV
556  real(kind=R_GRID) :: dx_local, dy_local
557 
558  real(kind=R_GRID) :: vec1(3), vec2(3), vec3(3), vec4(3)
559  real(kind=R_GRID) :: vecAvg(3), vec3a(3), vec3b(3), vec4a(3), vec4b(3)
560  real(kind=R_GRID) :: xyz1(3), xyz2(3)
561 
562 ! real(kind=R_GRID) :: grid_global(1-ng:npx +ng,1-ng:npy +ng,ndims,1:nregions)
563  integer :: ios, ip, jp
564 
565  integer :: igrid
566 
567  integer :: tmplun
568  character(len=80) :: tmpFile
569 
570  real(kind=R_GRID), dimension(Atm%bd%is:Atm%bd%ie) :: sbuffer, nbuffer
571  real(kind=R_GRID), dimension(Atm%bd%js:Atm%bd%je) :: wbuffer, ebuffer
572 
573  real(kind=R_GRID), pointer, dimension(:,:,:) :: agrid, grid
574  real(kind=R_GRID), pointer, dimension(:,:) :: area, area_c
575 
576  real(kind=R_GRID), pointer, dimension(:,:) :: sina, cosa, dx, dy, dxc, dyc, dxa, dya
577  real, pointer, dimension(:,:,:) :: e1, e2
578 
579  real, pointer, dimension(:,:) :: rarea, rarea_c
580  real, pointer, dimension(:,:) :: rdx, rdy, rdxc, rdyc, rdxa, rdya
581 
582  integer, pointer, dimension(:,:,:) :: iinta, jinta, iintb, jintb
583 
584  real(kind=R_GRID), pointer, dimension(:,:,:,:) :: grid_global
585 
586  integer, pointer :: npx_g, npy_g, ntiles_g, tile
587  logical, pointer :: sw_corner, se_corner, ne_corner, nw_corner
588  logical, pointer :: latlon, cubed_sphere, have_south_pole, have_north_pole, stretched_grid
589 
590  type(domain2d), pointer :: domain
591  integer :: is, ie, js, je
592  integer :: isd, ied, jsd, jed
593  integer :: istart, iend, jstart, jend
594 
595  is = atm%bd%is
596  ie = atm%bd%ie
597  js = atm%bd%js
598  je = atm%bd%je
599  isd = atm%bd%isd
600  ied = atm%bd%ied
601  jsd = atm%bd%jsd
602  jed = atm%bd%jed
603 
604  !!! Associate pointers
605  agrid => atm%gridstruct%agrid_64
606  grid => atm%gridstruct%grid_64
607 
608  area => atm%gridstruct%area_64
609  area_c => atm%gridstruct%area_c_64
610  rarea => atm%gridstruct%rarea
611  rarea_c => atm%gridstruct%rarea_c
612 
613  sina => atm%gridstruct%sina_64
614  cosa => atm%gridstruct%cosa_64
615  dx => atm%gridstruct%dx_64
616  dy => atm%gridstruct%dy_64
617  dxc => atm%gridstruct%dxc_64
618  dyc => atm%gridstruct%dyc_64
619  dxa => atm%gridstruct%dxa_64
620  dya => atm%gridstruct%dya_64
621  rdx => atm%gridstruct%rdx
622  rdy => atm%gridstruct%rdy
623  rdxc => atm%gridstruct%rdxc
624  rdyc => atm%gridstruct%rdyc
625  rdxa => atm%gridstruct%rdxa
626  rdya => atm%gridstruct%rdya
627  e1 => atm%gridstruct%e1
628  e2 => atm%gridstruct%e2
629 
630  if (atm%neststruct%nested .or. any(atm%neststruct%child_grids)) then
631  grid_global => atm%grid_global
632  else if( trim(grid_file) .EQ. 'Inline') then
633  allocate(grid_global(1-ng:npx +ng,1-ng:npy +ng,ndims,1:nregions))
634  endif
635 
636  iinta => atm%gridstruct%iinta
637  jinta => atm%gridstruct%jinta
638  iintb => atm%gridstruct%iintb
639  jintb => atm%gridstruct%jintb
640  npx_g => atm%gridstruct%npx_g
641  npy_g => atm%gridstruct%npy_g
642  ntiles_g => atm%gridstruct%ntiles_g
643  sw_corner => atm%gridstruct%sw_corner
644  se_corner => atm%gridstruct%se_corner
645  ne_corner => atm%gridstruct%ne_corner
646  nw_corner => atm%gridstruct%nw_corner
647  latlon => atm%gridstruct%latlon
648  cubed_sphere => atm%gridstruct%cubed_sphere
649  have_south_pole => atm%gridstruct%have_south_pole
650  have_north_pole => atm%gridstruct%have_north_pole
651  stretched_grid => atm%gridstruct%stretched_grid
652 
653  tile => atm%tile_of_mosaic
654 
655  domain => atm%domain
656 
657  npx_g = npx
658  npy_g = npy
659  ntiles_g = nregions
660  latlon = .false.
661  cubed_sphere = .false.
662 
663  if ( (atm%flagstruct%do_schmidt .or. atm%flagstruct%do_cube_transform) .and. abs(atm%flagstruct%stretch_fac-1.) > 1.e-5 ) then
664  stretched_grid = .true.
665  if (atm%flagstruct%do_schmidt .and. atm%flagstruct%do_cube_transform) then
666  call mpp_error(fatal, ' Cannot set both do_schmidt and do_cube_transform to .true.')
667  endif
668  endif
669 
670  if (atm%flagstruct%grid_type>3) then
671  if (atm%flagstruct%grid_type == 4) then
672  call setup_cartesian(npx, npy, atm%flagstruct%dx_const, atm%flagstruct%dy_const, &
673  atm%flagstruct%deglat, atm%bd)
674  else
675  call mpp_error(fatal, 'init_grid: unsupported grid type')
676  endif
677  else
678 
679  cubed_sphere = .true.
680 
681  if (atm%neststruct%nested) then
682  !Read grid if it exists
683  ! still need to set up
684  call setup_aligned_nest(atm)
685  else
686  if(trim(grid_file) .NE. 'Inline') then
687  call read_grid(atm, grid_file, ndims, nregions, ng)
688  else
689  if (atm%flagstruct%grid_type>=0) call gnomonic_grids(atm%flagstruct%grid_type, npx-1, xs, ys)
690  if (is_master()) then
691  if (atm%flagstruct%grid_type>=0) then
692  do j=1,npy
693  do i=1,npx
694  grid_global(i,j,1,1) = xs(i,j)
695  grid_global(i,j,2,1) = ys(i,j)
696  enddo
697  enddo
698 ! mirror_grid assumes that the tile=1 is centered on equator and greenwich meridian Lon[-pi,pi]
699  call mirror_grid(grid_global, ng, npx, npy, 2, 6)
700  do n=1,nregions
701  do j=1,npy
702  do i=1,npx
703 !---------------------------------
704 ! Shift the corner away from Japan
705 !---------------------------------
706 !--------------------- This will result in the corner close to east coast of China ------------------
707  if ( .not. ( atm%flagstruct%do_schmidt .or. atm%flagstruct%do_cube_transform) .and. (atm%flagstruct%shift_fac)>1.e-4 ) &
708  grid_global(i,j,1,n) = grid_global(i,j,1,n) - pi/atm%flagstruct%shift_fac
709 !----------------------------------------------------------------------------------------------------
710  if ( grid_global(i,j,1,n) < 0. ) &
711  grid_global(i,j,1,n) = grid_global(i,j,1,n) + 2.*pi
712  if (abs(grid_global(i,j,1,1)) < 1.d-10) grid_global(i,j,1,1) = 0.0
713  if (abs(grid_global(i,j,2,1)) < 1.d-10) grid_global(i,j,2,1) = 0.0
714  enddo
715  enddo
716  enddo
717  else
718  call mpp_error(fatal, "fv_grid_tools: reading of ASCII grid files no longer supported")
719  endif ! Atm%flagstruct%grid_type>=0
720 
721  grid_global( 1,1:npy,:,2)=grid_global(npx,1:npy,:,1)
722  grid_global( 1,1:npy,:,3)=grid_global(npx:1:-1,npy,:,1)
723  grid_global(1:npx,npy,:,5)=grid_global(1,npy:1:-1,:,1)
724  grid_global(1:npx,npy,:,6)=grid_global(1:npx,1,:,1)
725 
726  grid_global(1:npx, 1,:,3)=grid_global(1:npx,npy,:,2)
727  grid_global(1:npx, 1,:,4)=grid_global(npx,npy:1:-1,:,2)
728  grid_global(npx,1:npy,:,6)=grid_global(npx:1:-1,1,:,2)
729 
730  grid_global( 1,1:npy,:,4)=grid_global(npx,1:npy,:,3)
731  grid_global( 1,1:npy,:,5)=grid_global(npx:1:-1,npy,:,3)
732 
733  grid_global(npx,1:npy,:,3)=grid_global(1,1:npy,:,4)
734  grid_global(1:npx, 1,:,5)=grid_global(1:npx,npy,:,4)
735  grid_global(1:npx, 1,:,6)=grid_global(npx,npy:1:-1,:,4)
736 
737  grid_global( 1,1:npy,:,6)=grid_global(npx,1:npy,:,5)
738 
739 !------------------------
740 ! Schmidt transformation:
741 !------------------------
742  if ( atm%flagstruct%do_schmidt ) then
743  do n=1,nregions
744  call direct_transform(atm%flagstruct%stretch_fac, 1, npx, 1, npy, &
745  atm%flagstruct%target_lon, atm%flagstruct%target_lat, &
746  n, grid_global(1:npx,1:npy,1,n), grid_global(1:npx,1:npy,2,n))
747  enddo
748  elseif (atm%flagstruct%do_cube_transform) then
749  do n=1,nregions
750  call cube_transform(atm%flagstruct%stretch_fac, 1, npx, 1, npy, &
751  atm%flagstruct%target_lon, atm%flagstruct%target_lat, &
752  n, grid_global(1:npx,1:npy,1,n), grid_global(1:npx,1:npy,2,n))
753  enddo
754  endif
755  endif !is master
756  call mpp_broadcast(grid_global, size(grid_global), mpp_root_pe())
757  !--- copy grid to compute domain
758  do n=1,ndims
759  do j=js,je+1
760  do i=is,ie+1
761  grid(i,j,n) = grid_global(i,j,n,tile)
762  enddo
763  enddo
764  enddo
765  endif !(trim(grid_file) == 'INPUT/grid_spec.nc')
766 !
767 ! SJL: For phys/exchange grid, etc
768 !
769  call mpp_update_domains( grid, atm%domain, position=corner)
770  if (.not. (atm%gridstruct%bounded_domain)) then
771  call fill_corners(grid(:,:,1), npx, npy, fill=xdir, bgrid=.true.)
772  call fill_corners(grid(:,:,2), npx, npy, fill=xdir, bgrid=.true.)
773  endif
774 
775  !--- dx and dy
776  if( .not. atm%gridstruct%bounded_domain) then
777  istart=is
778  iend=ie
779  jstart=js
780  jend=je
781  else
782  istart=isd
783  iend=ied
784  jstart=jsd
785  jend=jed
786  endif
787 
788  do j = jstart, jend+1
789  do i = istart, iend
790  p1(1) = grid(i ,j,1)
791  p1(2) = grid(i ,j,2)
792  p2(1) = grid(i+1,j,1)
793  p2(2) = grid(i+1,j,2)
794  dx(i,j) = great_circle_dist( p2, p1, radius )
795  enddo
796  enddo
797  if( stretched_grid .or. atm%gridstruct%bounded_domain ) then
798  do j = jstart, jend
799  do i = istart, iend+1
800  p1(1) = grid(i,j, 1)
801  p1(2) = grid(i,j, 2)
802  p2(1) = grid(i,j+1,1)
803  p2(2) = grid(i,j+1,2)
804  dy(i,j) = great_circle_dist( p2, p1, radius )
805  enddo
806  enddo
807  else
808  call get_symmetry(dx(is:ie,js:je+1), dy(is:ie+1,js:je), 0, 1, atm%layout(1), atm%layout(2), &
809  atm%domain, atm%tile_of_mosaic, atm%gridstruct%npx_g, atm%bd)
810  endif
811 
812  call mpp_get_boundary( dy, dx, atm%domain, ebufferx=ebuffer, wbufferx=wbuffer, sbuffery=sbuffer, nbuffery=nbuffer,&
813  flags=scalar_pair+xupdate, gridtype=cgrid_ne_param)
814  if( .not. atm%gridstruct%bounded_domain ) then
815  if(is == 1 .AND. mod(tile,2) .NE. 0) then ! on the west boundary
816  dy(is, js:je) = wbuffer(js:je)
817  endif
818  if(ie == npx-1) then ! on the east boundary
819  dy(ie+1, js:je) = ebuffer(js:je)
820  endif
821  endif
822 
823  call mpp_update_domains( dy, dx, atm%domain, flags=scalar_pair, &
824  gridtype=cgrid_ne_param, complete=.true.)
825  if (cubed_sphere .and. (.not. (atm%gridstruct%bounded_domain))) then
826  call fill_corners(dx, dy, npx, npy, dgrid=.true.)
827  endif
828 
829  if( .not. stretched_grid ) &
830  call sorted_inta(isd, ied, jsd, jed, cubed_sphere, grid, iinta, jinta)
831 
832  agrid(:,:,:) = -1.e25
833 
834  !--- compute agrid (use same indices as for dx/dy above)
835 
836  do j=jstart,jend
837  do i=istart,iend
838  if ( stretched_grid ) then
839  call cell_center2(grid(i,j, 1:2), grid(i+1,j, 1:2), &
840  grid(i,j+1,1:2), grid(i+1,j+1,1:2), &
841  agrid(i,j,1:2) )
842  else
843  call cell_center2(grid(iinta(1,i,j),jinta(1,i,j),1:2), &
844  grid(iinta(2,i,j),jinta(2,i,j),1:2), &
845  grid(iinta(3,i,j),jinta(3,i,j),1:2), &
846  grid(iinta(4,i,j),jinta(4,i,j),1:2), &
847  agrid(i,j,1:2) )
848  endif
849  enddo
850  enddo
851 
852  call mpp_update_domains( agrid, atm%domain, position=center, complete=.true. )
853  if (.not. (atm%gridstruct%bounded_domain)) then
854  call fill_corners(agrid(:,:,1), npx, npy, xdir, agrid=.true.)
855  call fill_corners(agrid(:,:,2), npx, npy, ydir, agrid=.true.)
856  endif
857 
858  do j=jsd,jed
859  do i=isd,ied
860  call mid_pt_sphere(grid(i, j,1:2), grid(i, j+1,1:2), p1)
861  call mid_pt_sphere(grid(i+1,j,1:2), grid(i+1,j+1,1:2), p2)
862  dxa(i,j) = great_circle_dist( p2, p1, radius )
863  !
864  call mid_pt_sphere(grid(i,j ,1:2), grid(i+1,j ,1:2), p1)
865  call mid_pt_sphere(grid(i,j+1,1:2), grid(i+1,j+1,1:2), p2)
866  dya(i,j) = great_circle_dist( p2, p1, radius )
867  enddo
868  enddo
869 ! call mpp_update_domains( dxa, dya, Atm%domain, flags=SCALAR_PAIR, gridtype=AGRID_PARAM)
870  if (cubed_sphere .and. (.not. (atm%gridstruct%bounded_domain))) then
871  call fill_corners(dxa, dya, npx, npy, agrid=.true.)
872  endif
873 
874 
875 
876  end if !if nested
877 
878  do j=jsd,jed
879  do i=isd+1,ied
880  dxc(i,j) = great_circle_dist(agrid(i,j,:), agrid(i-1,j,:), radius)
881  enddo
882 !xxxxxx
883  !Are the following 2 lines appropriate for the regional domain?
884 !xxxxxx
885  dxc(isd,j) = dxc(isd+1,j)
886  dxc(ied+1,j) = dxc(ied,j)
887  enddo
888 
889 ! do j=js,je+1
890 ! do i=is,ie
891  do j=jsd+1,jed
892  do i=isd,ied
893  dyc(i,j) = great_circle_dist(agrid(i,j,:), agrid(i,j-1,:), radius)
894  enddo
895  enddo
896 !xxxxxx
897  !Are the following 2 lines appropriate for the regional domain?
898 !xxxxxx
899  do i=isd,ied
900  dyc(i,jsd) = dyc(i,jsd+1)
901  dyc(i,jed+1) = dyc(i,jed)
902  end do
903 
904 
905  if( .not. stretched_grid ) &
906  call sorted_intb(isd, ied, jsd, jed, is, ie, js, je, npx, npy, &
907  cubed_sphere, agrid, iintb, jintb)
908 
909  call grid_area( npx, npy, ndims, nregions, atm%gridstruct%bounded_domain, atm%gridstruct, atm%domain, atm%bd )
910 ! stretched_grid = .false.
911 
912 !----------------------------------
913 ! Compute area_c, rarea_c, dxc, dyc
914 !----------------------------------
915  if ( .not. stretched_grid .and. (.not. (atm%gridstruct%bounded_domain))) then
916 ! For symmetrical grids:
917  if ( is==1 ) then
918  i = 1
919  do j=js,je+1
920  call mid_pt_sphere(grid(i,j-1,1:2), grid(i,j, 1:2), p1)
921  call mid_pt_sphere(grid(i,j ,1:2), grid(i,j+1,1:2), p4)
922  p2(1:2) = agrid(i,j-1,1:2)
923  p3(1:2) = agrid(i,j, 1:2)
924  area_c(i,j) = 2.*get_area(p1, p4, p2, p3, radius)
925  enddo
926  do j=js,je
927  call mid_pt_sphere(grid(i,j,1:2), grid(i,j+1,1:2), p1)
928  p2(1:2) = agrid(i,j,1:2)
929  dxc(i,j) = 2.*great_circle_dist( p1, p2, radius )
930  enddo
931  endif
932  if ( (ie+1)==npx ) then
933  i = npx
934  do j=js,je+1
935  p1(1:2) = agrid(i-1,j-1,1:2)
936  call mid_pt_sphere(grid(i,j-1,1:2), grid(i,j, 1:2), p2)
937  call mid_pt_sphere(grid(i,j ,1:2), grid(i,j+1,1:2), p3)
938  p4(1:2) = agrid(i-1,j,1:2)
939  area_c(i,j) = 2.*get_area(p1, p4, p2, p3, radius)
940  enddo
941  do j=js,je
942  p1(1:2) = agrid(i-1,j,1:2)
943  call mid_pt_sphere(grid(i,j,1:2), grid(i,j+1,1:2), p2)
944  dxc(i,j) = 2.*great_circle_dist( p1, p2, radius )
945  enddo
946  endif
947  if ( js==1 ) then
948  j = 1
949  do i=is,ie+1
950  call mid_pt_sphere(grid(i-1,j,1:2), grid(i, j,1:2), p1)
951  call mid_pt_sphere(grid(i, j,1:2), grid(i+1,j,1:2), p2)
952  p3(1:2) = agrid(i, j,1:2)
953  p4(1:2) = agrid(i-1,j,1:2)
954  area_c(i,j) = 2.*get_area(p1, p4, p2, p3, radius)
955  enddo
956  do i=is,ie
957  call mid_pt_sphere(grid(i,j,1:2), grid(i+1,j,1:2), p1)
958  p2(1:2) = agrid(i,j,1:2)
959  dyc(i,j) = 2.*great_circle_dist( p1, p2, radius )
960  enddo
961  endif
962  if ( (je+1)==npy ) then
963  j = npy
964  do i=is,ie+1
965  p1(1:2) = agrid(i-1,j-1,1:2)
966  p2(1:2) = agrid(i ,j-1,1:2)
967  call mid_pt_sphere(grid(i, j,1:2), grid(i+1,j,1:2), p3)
968  call mid_pt_sphere(grid(i-1,j,1:2), grid(i, j,1:2), p4)
969  area_c(i,j) = 2.*get_area(p1, p4, p2, p3, radius)
970  enddo
971  do i=is,ie
972  p1(1:2) = agrid(i,j-1,1:2)
973  call mid_pt_sphere(grid(i,j,1:2), grid(i+1,j,1:2), p2)
974  dyc(i,j) = 2.*great_circle_dist( p1, p2, radius )
975  enddo
976  endif
977 
978  if ( sw_corner ) then
979  i=1; j=1
980  p1(1:2) = grid(i,j,1:2)
981  call mid_pt_sphere(grid(i,j,1:2), grid(i+1,j,1:2), p2)
982  p3(1:2) = agrid(i,j,1:2)
983  call mid_pt_sphere(grid(i,j,1:2), grid(i,j+1,1:2), p4)
984  area_c(i,j) = 3.*get_area(p1, p4, p2, p3, radius)
985  endif
986  if ( se_corner ) then
987  i=npx; j=1
988  call mid_pt_sphere(grid(i-1,j,1:2), grid(i,j,1:2), p1)
989  p2(1:2) = grid(i,j,1:2)
990  call mid_pt_sphere(grid(i,j,1:2), grid(i,j+1,1:2), p3)
991  p4(1:2) = agrid(i,j,1:2)
992  area_c(i,j) = 3.*get_area(p1, p4, p2, p3, radius)
993  endif
994  if ( ne_corner ) then
995  i=npx; j=npy
996  p1(1:2) = agrid(i-1,j-1,1:2)
997  call mid_pt_sphere(grid(i,j-1,1:2), grid(i,j,1:2), p2)
998  p3(1:2) = grid(i,j,1:2)
999  call mid_pt_sphere(grid(i-1,j,1:2), grid(i,j,1:2), p4)
1000  area_c(i,j) = 3.*get_area(p1, p4, p2, p3, radius)
1001  endif
1002  if ( nw_corner ) then
1003  i=1; j=npy
1004  call mid_pt_sphere(grid(i,j-1,1:2), grid(i,j,1:2), p1)
1005  p2(1:2) = agrid(i,j-1,1:2)
1006  call mid_pt_sphere(grid(i,j,1:2), grid(i+1,j,1:2), p3)
1007  p4(1:2) = grid(i,j,1:2)
1008  area_c(i,j) = 3.*get_area(p1, p4, p2, p3, radius)
1009  endif
1010  endif
1011 !-----------------
1012 
1013  call mpp_update_domains( dxc, dyc, atm%domain, flags=scalar_pair, &
1014  gridtype=cgrid_ne_param, complete=.true.)
1015  if (cubed_sphere .and. (.not. (atm%gridstruct%bounded_domain))) then
1016  call fill_corners(dxc, dyc, npx, npy, cgrid=.true.)
1017  endif
1018 
1019  call mpp_update_domains( area, atm%domain, complete=.true. )
1020 
1021 
1022  !Handling outermost ends for area_c
1023  if (atm%gridstruct%bounded_domain) then
1024  if (is == 1) then
1025  do j=jsd,jed
1026  area_c(isd,j) = area_c(isd+1,j)
1027  end do
1028  if (js == 1) area_c(isd,jsd) = area_c(isd+1,jsd+1)
1029  if (js == npy-1) area_c(isd,jed+1) = area_c(isd+1,jed)
1030  end if
1031  if (ie == npx-1) then
1032  do j=jsd,jed
1033  area_c(ied+1,j) = area_c(ied,j)
1034  end do
1035  if (js == 1) area_c(ied+1,jsd) = area_c(ied,jsd+1)
1036  if (js == npy-1) area_c(ied+1,jed+1) = area_c(ied,jed)
1037  end if
1038  if (js == 1) then
1039  do i=isd,ied
1040  area_c(i,jsd) = area_c(i,jsd+1)
1041  end do
1042  end if
1043  if (je == npy-1) then
1044  do i=isd,ied
1045  area_c(i,jed+1) = area_c(i,jed)
1046  end do
1047  end if
1048  end if
1049 
1050  call mpp_update_domains( area_c, atm%domain, position=corner, complete=.true.)
1051 
1052  ! Handle corner Area ghosting
1053  if (cubed_sphere .and. (.not. (atm%gridstruct%bounded_domain))) then
1054  call fill_ghost(area, npx, npy, -big_number, atm%bd) ! fill in garbage values
1055  call fill_corners(area_c, npx, npy, fill=xdir, bgrid=.true.)
1056  endif
1057 
1058  do j=jsd,jed+1
1059  do i=isd,ied
1060  rdx(i,j) = 1.0/dx(i,j)
1061  enddo
1062  enddo
1063  do j=jsd,jed
1064  do i=isd,ied+1
1065  rdy(i,j) = 1.0/dy(i,j)
1066  enddo
1067  enddo
1068  do j=jsd,jed
1069  do i=isd,ied+1
1070  rdxc(i,j) = 1.0/dxc(i,j)
1071  enddo
1072  enddo
1073  do j=jsd,jed+1
1074  do i=isd,ied
1075  rdyc(i,j) = 1.0/dyc(i,j)
1076  enddo
1077  enddo
1078  do j=jsd,jed
1079  do i=isd,ied
1080  rarea(i,j) = 1.0/area(i,j)
1081  rdxa(i,j) = 1./dxa(i,j)
1082  rdya(i,j) = 1./dya(i,j)
1083  enddo
1084  enddo
1085  do j=jsd,jed+1
1086  do i=isd,ied+1
1087  rarea_c(i,j) = 1.0/area_c(i,j)
1088  enddo
1089  enddo
1090 
1091 200 format(a,f9.2,a,f9.2,a,f9.2)
1092 201 format(a,f9.2,a,f9.2,a,f9.2,a,f9.2)
1093 202 format(a,a,i4.4,a,i4.4,a)
1094 
1095  ! Get and print Grid Statistics
1096  dxav =0.0
1097  angav=0.0
1098  aspav=0.0
1099  dxn = missing
1100  dxm = -missing
1101  angn = missing
1102  angm = -missing
1103  aspn = missing
1104  aspm = -missing
1105  !if (tile == 1) then ! doing a GLOBAL domain search on each grid
1106  do j=js, je
1107  do i=is, ie
1108  if(i>ceiling(npx/2.) .OR. j>ceiling(npy/2.)) cycle
1109  ang = get_angle(2, grid(i,j+1,1:2), grid(i,j,1:2), grid(i+1,j,1:2))
1110  ang = abs(90.0 - ang)
1111 
1112  if ( (i==1) .and. (j==1) ) then
1113  else
1114  angav = angav + ang
1115  angm = max(angm,ang)
1116  angn = min(angn,ang)
1117  endif
1118 
1119  dx_local = dx(i,j)
1120  dy_local = dy(i,j)
1121 
1122  dxav = dxav + 0.5 * (dx_local + dy_local)
1123  dxm = max(dxm,dx_local)
1124  dxm = max(dxm,dy_local)
1125  dxn = min(dxn,dx_local)
1126  dxn = min(dxn,dy_local)
1127 
1128  asp = abs(dx_local/dy_local)
1129  if (asp < 1.0) asp = 1.0/asp
1130  aspav = aspav + asp
1131  aspm = max(aspm,asp)
1132  aspn = min(aspn,asp)
1133  enddo
1134  enddo
1135  !endif
1136  call mpp_sum(angav)
1137  call mpp_sum(dxav)
1138  call mpp_sum(aspav)
1139  call mpp_max(angm)
1140  call mpp_min(angn)
1141  call mpp_max(dxm)
1142  call mpp_min(dxn)
1143  call mpp_max(aspm)
1144  call mpp_min(aspn)
1145 
1146  if( is_master() ) then
1147 
1148  angav = angav / ( (ceiling(npy/2.0))*(ceiling(npx/2.0)) - 1 )
1149  dxav = dxav / ( (ceiling(npy/2.0))*(ceiling(npx/2.0)) )
1150  aspav = aspav / ( (ceiling(npy/2.0))*(ceiling(npx/2.0)) )
1151  write(*,* ) ''
1152 #ifdef SMALL_EARTH
1153  write(*,*) ' REDUCED EARTH: Radius is ', radius, ', omega is ', omega
1154 #endif
1155  write(*,* ) ' Cubed-Sphere Grid Stats : ', npx,'x',npy,'x',nregions
1156  print*, dxn, dxm, dxav, dxn, dxm
1157  write(*,201) ' Grid Length : min: ', dxn,' max: ', dxm,' avg: ', dxav, ' min/max: ',dxn/dxm
1158  write(*,200) ' Deviation from Orthogonal : min: ',angn,' max: ',angm,' avg: ',angav
1159  write(*,200) ' Aspect Ratio : min: ',aspn,' max: ',aspm,' avg: ',aspav
1160  write(*,* ) ''
1161  endif
1162  endif!if gridtype > 3
1163 
1164  !SEND grid global if any child nests
1165  !Matching receive in setup_aligned_nest
1166  do n=1,size(atm%neststruct%child_grids)
1167  if (atm%neststruct%child_grids(n) .and. is_master()) then
1168  !need to get tile_coarse AND determine local number for tile
1169  if (ntiles_g > 1) then ! coarse grid only!!
1170 !!$ !!! DEBUG CODE
1171 !!$ print*, 'SENDING GRID_GLOBAL: ', mpp_pe(), tile_coarse(n), grids_master_procs(n), grid_global(1,npy,:,tile_coarse(n))
1172 !!$ !!! END DEBUG CODE
1173  call mpp_send(grid_global(:,:,:,tile_coarse(n)), &
1174  size(grid_global)/atm%flagstruct%ntiles,grids_master_procs(n))
1175  else
1176  call mpp_send(grid_global(:,:,:,1),size(grid_global),grids_master_procs(n))
1177  endif
1178  call mpp_sync_self()
1179  endif
1180  enddo
1181 
1182  if (atm%neststruct%nested .or. any(atm%neststruct%child_grids)) then
1183  nullify(grid_global)
1184  else if( trim(grid_file) .EQ. 'Inline') then
1185  deallocate(grid_global)
1186  endif
1187 
1188  nullify(agrid)
1189  nullify(grid)
1190 
1191  nullify( area)
1192  nullify(rarea)
1193  nullify( area_c)
1194  nullify(rarea_c)
1195 
1196  nullify(sina)
1197  nullify(cosa)
1198  nullify(dx)
1199  nullify(dy)
1200  nullify(dxc)
1201  nullify(dyc)
1202  nullify(dxa)
1203  nullify(dya)
1204  nullify(rdx)
1205  nullify(rdy)
1206  nullify(rdxc)
1207  nullify(rdyc)
1208  nullify(rdxa)
1209  nullify(rdya)
1210  nullify(e1)
1211  nullify(e2)
1212 
1213  nullify(iinta)
1214  nullify(jinta)
1215  nullify(iintb)
1216  nullify(jintb)
1217  nullify(npx_g)
1218  nullify(npy_g)
1219  nullify(ntiles_g)
1220  nullify(sw_corner)
1221  nullify(se_corner)
1222  nullify(ne_corner)
1223  nullify(nw_corner)
1224  nullify(latlon)
1225  nullify(cubed_sphere)
1226  nullify(have_south_pole)
1227  nullify(have_north_pole)
1228  nullify(stretched_grid)
1229 
1230  nullify(tile)
1231 
1232  nullify(domain)
1233 
1234  contains
1235 
1236  subroutine setup_cartesian(npx, npy, dx_const, dy_const, deglat, bd)
1237 
1238  type(fv_grid_bounds_type), intent(IN) :: bd
1239  integer, intent(in):: npx, npy
1240  real(kind=R_GRID), intent(IN) :: dx_const, dy_const, deglat
1241  real(kind=R_GRID) lat_rad, lon_rad, domain_rad
1242  integer i,j
1243  integer :: is, ie, js, je
1244  integer :: isd, ied, jsd, jed
1245 
1246  is = bd%is
1247  ie = bd%ie
1248  js = bd%js
1249  je = bd%je
1250  isd = bd%isd
1251  ied = bd%ied
1252  jsd = bd%jsd
1253  jed = bd%jed
1254 
1255  domain_rad = pi/16. ! arbitrary
1256  lat_rad = deglat * pi/180.
1257  lon_rad = 0. ! arbitrary
1258 
1259  dx(:,:) = dx_const
1260  rdx(:,:) = 1./dx_const
1261  dy(:,:) = dy_const
1262  rdy(:,:) = 1./dy_const
1263 
1264  dxc(:,:) = dx_const
1265  rdxc(:,:) = 1./dx_const
1266  dyc(:,:) = dy_const
1267  rdyc(:,:) = 1./dy_const
1268 
1269  dxa(:,:) = dx_const
1270  rdxa(:,:) = 1./dx_const
1271  dya(:,:) = dy_const
1272  rdya(:,:) = 1./dy_const
1273 
1274  area(:,:) = dx_const*dy_const
1275  rarea(:,:) = 1./(dx_const*dy_const)
1276 
1277  area_c(:,:) = dx_const*dy_const
1278  rarea_c(:,:) = 1./(dx_const*dy_const)
1279 
1280 ! The following is a hack to get pass the am2 phys init:
1281  do j=max(1,jsd),min(jed,npy)
1282  do i=max(1,isd),min(ied,npx)
1283  grid(i,j,1) = lon_rad - 0.5*domain_rad + real(i-1)/real(npx-1)*domain_rad
1284  grid(i,j,2) = lat_rad - 0.5*domain_rad + real(j-1)/real(npy-1)*domain_rad
1285  enddo
1286  enddo
1287 
1288  agrid(:,:,1) = lon_rad
1289  agrid(:,:,2) = lat_rad
1290 
1291  sina(:,:) = 1.
1292  cosa(:,:) = 0.
1293 
1294  e1(1,:,:) = 1.
1295  e1(2,:,:) = 0.
1296  e1(3,:,:) = 0.
1297 
1298  e2(1,:,:) = 0.
1299  e2(2,:,:) = 1.
1300  e2(3,:,:) = 0.
1301 
1302  end subroutine setup_cartesian
1303 
1304  !This routine currently does two things:
1305  ! 1) Create the nested grid on-the-fly from the parent
1306  ! 2) Compute the weights and indices for the boundary conditions
1307  ! We should split these into two routines in case we can
1308  ! read the nest from the input mosaic. Then we only need
1309  ! to set up the weights.
1310  ! When creating the nest on-the-fly we need the global parent grid,
1311  ! as we are doing now. For nests crossing a cube edge
1312  ! new code is needed.
1313  ! Creating the indices should be relatvely straightforward procedure
1314  ! since we will always know ioffset and joffset, which are needed
1315  ! to initialize the mpp nesting structure
1316  ! Computing the weights can be simplified by simply retreiving the
1317  ! BC agrid/grid structures?
1318 
1319  subroutine setup_aligned_nest(Atm)
1320 
1321  type(fv_atmos_type), intent(INOUT), target :: Atm
1322 
1323  integer :: isd_p, ied_p, jsd_p, jed_p
1324  integer :: isg, ieg, jsg, jeg
1325  integer :: ic, jc, imod, jmod
1326 
1327 
1328  real(kind=R_GRID), allocatable, dimension(:,:,:) :: p_grid_u, p_grid_v, pa_grid, p_grid, c_grid_u, c_grid_v
1329  integer :: p_ind(1-ng:npx +ng,1-ng:npy +ng,4)
1332 
1333  integer i,j,k, p
1334  real(kind=R_GRID) sum
1335  real(kind=R_GRID) :: dist1, dist2, dist3, dist4
1336  real(kind=R_GRID), dimension(2) :: q1, q2
1337 
1338  integer, pointer :: parent_tile, refinement, ioffset, joffset
1339  integer, pointer, dimension(:,:,:) :: ind_h, ind_u, ind_v
1340  real, pointer, dimension(:,:,:) :: wt_h, wt_u, wt_v
1341 
1342  integer, pointer, dimension(:,:,:) :: ind_b
1343  real, pointer, dimension(:,:,:) :: wt_b
1344 
1345  integer :: is, ie, js, je
1346  integer :: isd, ied, jsd, jed
1347 
1348  is = atm%bd%is
1349  ie = atm%bd%ie
1350  js = atm%bd%js
1351  je = atm%bd%je
1352  isd = atm%bd%isd
1353  ied = atm%bd%ied
1354  jsd = atm%bd%jsd
1355  jed = atm%bd%jed
1356 
1357 
1358 
1359  parent_tile => atm%neststruct%parent_tile
1360  refinement => atm%neststruct%refinement
1361  ioffset => atm%neststruct%ioffset
1362  joffset => atm%neststruct%joffset
1363 
1364  ind_h => atm%neststruct%ind_h
1365  ind_u => atm%neststruct%ind_u
1366  ind_v => atm%neststruct%ind_v
1367 
1368  wt_h => atm%neststruct%wt_h
1369  wt_u => atm%neststruct%wt_u
1370  wt_v => atm%neststruct%wt_v
1371 
1372  ind_b => atm%neststruct%ind_b
1373  wt_b => atm%neststruct%wt_b
1374 
1375  call mpp_get_data_domain( atm%parent_grid%domain, &
1376  isd_p, ied_p, jsd_p, jed_p )
1377  call mpp_get_global_domain( atm%parent_grid%domain, &
1378  isg, ieg, jsg, jeg)
1379 
1380  allocate(p_grid_u(isg:ieg ,jsg:jeg+1,1:2))
1381  allocate(p_grid_v(isg:ieg+1,jsg:jeg ,1:2))
1382  allocate(pa_grid(isg:ieg,jsg:jeg ,1:2))
1383  p_ind = -1000000000
1384 
1385  allocate(p_grid( isg-ng:ieg+1+ng, jsg-ng:jeg+1+ng,1:2) )
1386  p_grid = 1.e25
1387 
1388  !Need to RECEIVE parent grid_global;
1389  !matching mpp_send of grid_global from parent grid is in init_grid()
1390  if( is_master() ) then
1391 
1392  call mpp_recv(p_grid( isg-ng:ieg+1+ng, jsg-ng:jeg+1+ng,1:2), size(p_grid( isg-ng:ieg+1+ng, jsg-ng:jeg+1+ng,1:2)), &
1393  atm%parent_grid%pelist(1))
1394 !!$ !!!! DEBUG CODE
1395 !!$ print*, 'RECEIVING GRID GLOBAL: ', mpp_pe(), Atm%parent_grid%pelist(1), p_grid(1,jeg+1,:)
1396 !!$ !!!! END DEBUG CODE
1397 
1398  endif
1399 
1400  call mpp_broadcast( p_grid(isg-ng:ieg+ng+1, jsg-ng:jeg+ng+1, :), &
1401  (ieg-isg+2+2*ng)*(jeg-jsg+2+2*ng)*ndims, mpp_root_pe() )
1402 
1403  !NOTE : Grid now allowed to lie outside of parent
1404  !Check that the grid does not lie outside its parent
1405  !3aug15: allows halo of nest to lie within halo of coarse grid.
1406 !!$ ! NOTE: will this then work with the mpp_update_nest_fine?
1407 !!$ if ( joffset + floor( real(1-ng) / real(refinement) ) < 1-ng .or. &
1408 !!$ ioffset + floor( real(1-ng) / real(refinement) ) < 1-ng .or. &
1409 !!$ joffset + floor( real(npy+ng) / real(refinement) ) > Atm%parent_grid%npy+ng .or. &
1410 !!$ ioffset + floor( real(npx+ng) / real(refinement) ) > Atm%parent_grid%npx+ng ) then
1411 !!$ call mpp_error(FATAL, 'nested grid lies outside its parent')
1412 !!$ end if
1413 
1414  do j=1-ng,npy+ng
1415  jc = joffset + (j-1)/refinement !int( real(j-1) / real(refinement) )
1416  jmod = mod(j-1,refinement)
1417  if (j-1 < 0 .and. jmod /= 0) jc = jc - 1
1418  if (jmod < 0) jmod = jmod + refinement
1419 
1420  do i=1-ng,npx+ng
1421  ic = ioffset + (i-1)/refinement !int( real(i-1) / real(refinement) )
1422  imod = mod(i-1,refinement)
1423  if (i-1 < 0 .and. imod /= 0) ic = ic - 1
1424  if (imod < 0) imod = imod + refinement
1425 
1426  if (ic+1 > ieg+1 .or. ic < isg .or. jc+1 > jeg+1 .or. jc < jsg) then
1427  print*, 'p_grid:', i, j, ' OUT OF BOUNDS'
1428  print*, ic, jc
1429  print*, isg, ieg, jsg, jeg
1430  print*, imod, jmod
1431  end if
1432 
1433  if (jmod == 0) then
1434  q1 = p_grid(ic, jc, 1:2)
1435  q2 = p_grid(ic+1,jc,1:2)
1436  else
1437  call spherical_linear_interpolation( real(jmod,kind=r_grid)/real(refinement,kind=R_GRID), &
1438  p_grid(ic, jc, 1:2), p_grid(ic, jc+1, 1:2), q1)
1439  call spherical_linear_interpolation( real(jmod,kind=r_grid)/real(refinement,kind=R_GRID), &
1440  p_grid(ic+1, jc, 1:2), p_grid(ic+1, jc+1, 1:2), q2)
1441  end if
1442 
1443  if (imod == 0) then
1444  grid_global(i,j,:,1) = q1
1445  else
1446  call spherical_linear_interpolation( real(imod,kind=r_grid)/real(refinement,kind=R_GRID), &
1447  q1,q2,grid_global(i,j,:,1))
1448  end if
1449 
1450  !SW coarse-grid index; assumes grid does
1451  !not overlie other cube panels. (These indices
1452  !are also for the corners and thus need modification
1453  !to be used for cell-centered and edge-
1454  !centered variables; see below)
1455  p_ind(i,j,1) = ic
1456  p_ind(i,j,2) = jc
1457  p_ind(i,j,3) = imod
1458  p_ind(i,j,4) = jmod
1459 
1460  if (grid_global(i,j,1,1) > 2.*pi) grid_global(i,j,1,1) = grid_global(i,j,1,1) - 2.*pi
1461  if (grid_global(i,j,1,1) < 0.) grid_global(i,j,1,1) = grid_global(i,j,1,1) + 2.*pi
1462 
1463  end do
1464  end do
1465 
1466  ! Set up parent grids for interpolation purposes
1467  do j=jsg,jeg+1
1468  do i=isg,ieg
1469  call mid_pt_sphere(p_grid(i, j,1:2), p_grid(i+1, j,1:2), p_grid_u(i,j,:))
1470  !call mid_pt_sphere(p_grid(i, j,1:2), p_grid(i, j+1,1:2), p_grid_u(i,j,:))
1471  end do
1472  end do
1473  do j=jsg,jeg
1474  do i=isg,ieg+1
1475  call mid_pt_sphere(p_grid(i, j,1:2), p_grid(i, j+1,1:2), p_grid_v(i,j,:))
1476  !call mid_pt_sphere(p_grid(i, j,1:2), p_grid(i+1, j,1:2), p_grid_v(i,j,:))
1477  end do
1478  end do
1479  do j=jsg,jeg
1480  do i=isg,ieg
1481  call cell_center2(p_grid(i,j, 1:2), p_grid(i+1,j, 1:2), &
1482  p_grid(i,j+1,1:2), p_grid(i+1,j+1,1:2), &
1483  pa_grid(i,j,1:2) )
1484  end do
1485  end do
1486 
1487 !!$ !TODO: can we just send around ONE grid and re-calculate
1488 !!$ ! staggered grids from that??
1489 !!$ call mpp_broadcast(grid_global(1-ng:npx+ng, 1-ng:npy+ng ,:,1), &
1490 !!$ ((npx+ng)-(1-ng)+1)*((npy+ng)-(1-ng)+1)*ndims, mpp_root_pe() )
1491 !!$ call mpp_broadcast( p_ind(1-ng:npx+ng, 1-ng:npy+ng ,1:4), &
1492 !!$ ((npx+ng)-(1-ng)+1)*((npy+ng)-(1-ng)+1)*4, mpp_root_pe() )
1493 !!$ call mpp_broadcast( pa_grid( isg:ieg , jsg:jeg , :), &
1494 !!$ ((ieg-isg+1))*(jeg-jsg+1)*ndims, mpp_root_pe())
1495 !!$ call mpp_broadcast( p_grid_u( isg:ieg , jsg:jeg+1, :), &
1496 !!$ (ieg-isg+1)*(jeg-jsg+2)*ndims, mpp_root_pe())
1497 !!$ call mpp_broadcast( p_grid_v( isg:ieg+1, jsg:jeg , :), &
1498 !!$ (ieg-isg+2)*(jeg-jsg+1)*ndims, mpp_root_pe())
1499 
1500  do n=1,ndims
1501  do j=jsd,jed+1
1502  do i=isd,ied+1
1503  grid(i,j,n) = grid_global(i,j,n,1)
1504  enddo
1505  enddo
1506  enddo
1507 
1508  ind_h = -999999999
1509  do j=jsd,jed
1510  do i=isd,ied
1511  ic = p_ind(i,j,1)
1512  jc = p_ind(i,j,2)
1513  imod = p_ind(i,j,3)
1514  jmod = p_ind(i,j,4)
1515 
1516 
1517  if (imod < refinement/2) then
1518 !!$ !!! DEBUG CODE
1519 !!$ if (ic /= ic) print*, gid, ' Bad ic ', i, j
1520 !!$ print*, i, j, ic
1521 !!$ !!! END DEBUG CODE
1522  ind_h(i,j,1) = ic - 1
1523  else
1524  ind_h(i,j,1) = ic
1525  end if
1526 
1527  if (jmod < refinement/2) then
1528  ind_h(i,j,2) = jc - 1
1529  else
1530  ind_h(i,j,2) = jc
1531  end if
1532  ind_h(i,j,3) = imod
1533  ind_h(i,j,4) = jmod
1534 
1535  end do
1536  end do
1537 
1538  ind_b = -999999999
1539  do j=jsd,jed+1
1540  do i=isd,ied+1
1541  ic = p_ind(i,j,1)
1542  jc = p_ind(i,j,2)
1543  imod = p_ind(i,j,3)
1544  jmod = p_ind(i,j,4)
1545 
1546  ind_b(i,j,1) = ic
1547  ind_b(i,j,2) = jc
1548 
1549  ind_b(i,j,3) = imod
1550  ind_b(i,j,4) = jmod
1551  enddo
1552  enddo
1553 
1554  ind_u = -99999999
1555  !New BCs for wind components:
1556  ! For aligned grid segments (mod(j-1,R) == 0) set
1557  ! identically equal to the coarse-grid value
1558  ! Do linear interpolation in the y-dir elsewhere
1559 
1560  do j=jsd,jed+1
1561  do i=isd,ied
1562  ic = p_ind(i,j,1)
1563  jc = p_ind(i,j,2)
1564  imod = p_ind(i,j,3)
1565 
1566 #ifdef NEW_BC
1567  ind_u(i,j,1) = ic
1568 #else
1569  if (imod < refinement/2) then
1570 !!$ !!! DEBUG CODE
1571 !!$ print*, i, j, ic
1572 !!$ !!! END DEBUG CODE
1573  ind_u(i,j,1) = ic - 1
1574  else
1575  ind_u(i,j,1) = ic
1576  end if
1577 #endif
1578 
1579  ind_u(i,j,2) = jc
1580  ind_u(i,j,3) = imod
1581  ind_u(i,j,4) = p_ind(i,j,4)
1582 
1583  end do
1584  end do
1585 
1586  ind_v = -999999999
1587 
1588  do j=jsd,jed
1589  do i=isd,ied+1
1590  ic = p_ind(i,j,1)
1591  jc = p_ind(i,j,2)
1592  jmod = p_ind(i,j,4)
1593 
1594  ind_v(i,j,1) = ic
1595 
1596 #ifdef NEW_BC
1597  ind_v(i,j,2) = jc
1598 #else
1599  if (jmod < refinement/2) then
1600  ind_v(i,j,2) = jc - 1
1601  else
1602  ind_v(i,j,2) = jc
1603  end if
1604 #endif
1605 
1606  ind_v(i,j,4) = jmod
1607  ind_v(i,j,3) = p_ind(i,j,3)
1608  end do
1609  end do
1610 
1611 
1612 
1613  agrid(:,:,:) = -1.e25
1614 
1615  do j=jsd,jed
1616  do i=isd,ied
1617  call cell_center2(grid(i,j, 1:2), grid(i+1,j, 1:2), &
1618  grid(i,j+1,1:2), grid(i+1,j+1,1:2), &
1619  agrid(i,j,1:2) )
1620  enddo
1621  enddo
1622 
1623  call mpp_update_domains( agrid, atm%domain, position=center, complete=.true. )
1624 
1625  ! Compute dx
1626  do j=jsd,jed+1
1627  do i=isd,ied
1628  dx(i,j) = great_circle_dist(grid_global(i,j,:,1), grid_global(i+1,j,:,1), radius)
1629  enddo
1630  enddo
1631 
1632  ! Compute dy
1633  do j=jsd,jed
1634  do i=isd,ied+1
1635  dy(i,j) = great_circle_dist(grid_global(i,j,:,1), grid_global(i,j+1,:,1), radius)
1636  enddo
1637  enddo
1638 
1639  !We will use Michael Herzog's algorithm for computing the weights.
1640 
1641  !h points - need center (a-grid) points of parent grid
1642  do j=jsd,jed
1643  do i=isd,ied
1644 
1645  ic = ind_h(i,j,1)
1646  jc = ind_h(i,j,2)
1647 
1648  dist1 = dist2side_latlon(pa_grid(ic,jc,:) ,pa_grid(ic,jc+1,:),agrid(i,j,:))
1649  dist2 = dist2side_latlon(pa_grid(ic,jc+1,:) ,pa_grid(ic+1,jc+1,:),agrid(i,j,:))
1650  dist3 = dist2side_latlon(pa_grid(ic+1,jc+1,:),pa_grid(ic+1,jc,:),agrid(i,j,:))
1651  dist4 = dist2side_latlon(pa_grid(ic,jc,:) ,pa_grid(ic+1,jc,:),agrid(i,j,:))
1652 
1653  wt_h(i,j,1)=dist2*dist3 ! ic, jc weight
1654  wt_h(i,j,2)=dist3*dist4 ! ic, jc+1 weight
1655  wt_h(i,j,3)=dist4*dist1 ! ic+1, jc+1 weight
1656  wt_h(i,j,4)=dist1*dist2 ! ic+1, jc weight
1657 
1658  sum=wt_h(i,j,1)+wt_h(i,j,2)+wt_h(i,j,3)+wt_h(i,j,4)
1659  wt_h(i,j,:)=wt_h(i,j,:)/sum
1660 
1661  end do
1662  end do
1663 
1664 
1665 
1666  deallocate(pa_grid)
1667 
1668  do j=jsd,jed+1
1669  do i=isd,ied+1
1670 
1671  ic = ind_b(i,j,1)
1672  jc = ind_b(i,j,2)
1673 
1674  dist1 = dist2side_latlon(p_grid(ic,jc,:), p_grid(ic,jc+1,:), grid(i,j,:))
1675  dist2 = dist2side_latlon(p_grid(ic,jc+1,:), p_grid(ic+1,jc+1,:), grid(i,j,:))
1676  dist3 = dist2side_latlon(p_grid(ic+1,jc+1,:), p_grid(ic+1,jc,:), grid(i,j,:))
1677  dist4 = dist2side_latlon(p_grid(ic,jc,:), p_grid(ic+1,jc,:), grid(i,j,:))
1678 
1679  wt_b(i,j,1)=dist2*dist3 ! ic, jc weight
1680  wt_b(i,j,2)=dist3*dist4 ! ic, jc+1 weight
1681  wt_b(i,j,3)=dist4*dist1 ! ic+1, jc+1 weight
1682  wt_b(i,j,4)=dist1*dist2 ! ic+1, jc weight
1683 
1684  sum=wt_b(i,j,1)+wt_b(i,j,2)+wt_b(i,j,3)+wt_b(i,j,4)
1685  wt_b(i,j,:)=wt_b(i,j,:)/sum
1686 
1687  enddo
1688  enddo
1689 
1690  deallocate(p_grid)
1691 
1692 
1693  allocate(c_grid_u(isd:ied+1,jsd:jed,2))
1694  allocate(c_grid_v(isd:ied,jsd:jed+1,2))
1695 
1696  do j=jsd,jed
1697  do i=isd,ied+1
1698  call mid_pt_sphere(grid(i, j,1:2), grid(i, j+1,1:2), c_grid_u(i,j,:))
1699  end do
1700  end do
1701 
1702  do j=jsd,jed+1
1703  do i=isd,ied
1704  call mid_pt_sphere(grid(i,j ,1:2), grid(i+1,j ,1:2), c_grid_v(i,j,:))
1705  end do
1706  end do
1707 
1708 
1709  do j=jsd,jed
1710  do i=isd,ied
1711  dxa(i,j) = great_circle_dist(c_grid_u(i,j,:), c_grid_u(i+1,j,:), radius)
1712  end do
1713  end do
1714 
1715  do j=jsd,jed
1716  do i=isd,ied
1717  dya(i,j) = great_circle_dist(c_grid_v(i,j,:), c_grid_v(i,j+1,:), radius)
1718  end do
1719  end do
1720 
1721 
1722  !Compute interpolation weights. (Recall that the weights are defined with respect to a d-grid)
1723 
1724  !U weights
1725 
1726  do j=jsd,jed+1
1727  do i=isd,ied
1728 
1729  ic = ind_u(i,j,1)
1730  jc = ind_u(i,j,2)
1731 
1732  if (ic+1 > ieg .or. ic < isg .or. jc+1 > jeg+1 .or. jc < jsg) then
1733  print*, 'IND_U ', i, j, ' OUT OF BOUNDS'
1734  print*, ic, jc
1735  print*, isg, ieg, jsg, jeg
1736  end if
1737 
1738 
1739 #ifdef NEW_BC
1740  !New vorticity-conserving weights. Note that for the C-grid winds these
1741  ! become divergence-conserving weights!!
1742  jmod = p_ind(i,j,4)
1743  if (jmod == 0) then
1744  wt_u(i,j,1) = 1. ; wt_u(i,j,2) = 0.
1745  else
1746  dist1 = dist2side_latlon(p_grid_u(ic,jc,:), p_grid_u(ic+1,jc,:), c_grid_v(i,j,:))
1747  dist2 = dist2side_latlon(p_grid_u(ic,jc+1,:), p_grid_u(ic+1,jc+1,:), c_grid_v(i,j,:))
1748  sum = dist1+dist2
1749  wt_u(i,j,1) = dist2/sum
1750  wt_u(i,j,2) = dist1/sum
1751  endif
1752  wt_u(i,j,3:4) = 0.
1753 #else
1754  dist1 = dist2side_latlon(p_grid_u(ic,jc,:) ,p_grid_u(ic,jc+1,:), c_grid_v(i,j,:))
1755  dist2 = dist2side_latlon(p_grid_u(ic,jc+1,:) ,p_grid_u(ic+1,jc+1,:),c_grid_v(i,j,:))
1756  dist3 = dist2side_latlon(p_grid_u(ic+1,jc+1,:),p_grid_u(ic+1,jc,:), c_grid_v(i,j,:))
1757  !dist4 = dist2side_latlon(p_grid_u(ic+1,jc,:) ,p_grid_u(ic,jc,:), c_grid_v(i,j,:))
1758  dist4 = dist2side_latlon(p_grid_u(ic,jc,:) ,p_grid_u(ic+1,jc,:), c_grid_v(i,j,:))
1759 
1760  wt_u(i,j,1)=dist2*dist3 ! ic, jc weight
1761  wt_u(i,j,2)=dist3*dist4 ! ic, jc+1 weight
1762  wt_u(i,j,3)=dist4*dist1 ! ic+1, jc+1 weight
1763  wt_u(i,j,4)=dist1*dist2 ! ic+1, jc weight
1764 
1765  sum=wt_u(i,j,1)+wt_u(i,j,2)+wt_u(i,j,3)+wt_u(i,j,4)
1766  wt_u(i,j,:)=wt_u(i,j,:)/sum
1767 #endif
1768 
1769  end do
1770  end do
1771  !v weights
1772 
1773  do j=jsd,jed
1774  do i=isd,ied+1
1775 
1776  ic = ind_v(i,j,1)
1777  jc = ind_v(i,j,2)
1778 
1779  if (ic+1 > ieg .or. ic < isg .or. jc+1 > jeg+1 .or. jc < jsg) then
1780  print*, 'IND_V ', i, j, ' OUT OF BOUNDS'
1781  print*, ic, jc
1782  print*, isg, ieg, jsg, jeg
1783  end if
1784 
1785 #ifdef NEW_BC
1786  imod = p_ind(i,j,3)
1787  if (imod == 0) then
1788  wt_v(i,j,1) = 1. ; wt_v(i,j,4) = 0.
1789  else
1790  dist1 = dist2side_latlon(p_grid_v(ic,jc,:), p_grid_v(ic,jc+1,:), c_grid_u(i,j,:))
1791  dist2 = dist2side_latlon(p_grid_v(ic+1,jc,:), p_grid_v(ic+1,jc+1,:), c_grid_u(i,j,:))
1792  sum = dist1+dist2
1793  wt_v(i,j,1) = dist2/sum
1794  wt_v(i,j,4) = dist1/sum
1795  endif
1796  wt_v(i,j,2) = 0. ; wt_v(i,j,3) = 0.
1797 #else
1798  dist1 = dist2side_latlon(p_grid_v(ic,jc,:) ,p_grid_v(ic,jc+1,:), c_grid_u(i,j,:))
1799  dist2 = dist2side_latlon(p_grid_v(ic,jc+1,:) ,p_grid_v(ic+1,jc+1,:),c_grid_u(i,j,:))
1800  dist3 = dist2side_latlon(p_grid_v(ic+1,jc+1,:),p_grid_v(ic+1,jc,:), c_grid_u(i,j,:))
1801  dist4 = dist2side_latlon(p_grid_v(ic,jc,:) ,p_grid_v(ic+1,jc,:), c_grid_u(i,j,:))
1802 
1803  wt_v(i,j,1)=dist2*dist3 ! ic, jc weight
1804  wt_v(i,j,2)=dist3*dist4 ! ic, jc+1 weight
1805  wt_v(i,j,3)=dist4*dist1 ! ic+1, jc+1 weight
1806  wt_v(i,j,4)=dist1*dist2 ! ic+1, jc weight
1807 
1808  sum=wt_v(i,j,1)+wt_v(i,j,2)+wt_v(i,j,3)+wt_v(i,j,4)
1809  wt_v(i,j,:)=wt_v(i,j,:)/sum
1810 #endif
1811 
1812  end do
1813  end do
1814 
1815  deallocate(c_grid_u)
1816  deallocate(c_grid_v)
1817 
1818 
1819  deallocate(p_grid_u)
1820  deallocate(p_grid_v)
1821 
1822  if (is_master()) then
1823  if (atm%neststruct%nested) then
1824  !Nesting position information
1825  write(*,*) 'NESTED GRID ', atm%grid_number
1826  ic = p_ind(1,1,1) ; jc = p_ind(1,1,1)
1827  write(*,'(A, 2I5, 4F10.4)') 'SW CORNER: ', ic, jc, grid_global(1,1,:,1)*90./pi
1828  ic = p_ind(1,npy,1) ; jc = p_ind(1,npy,1)
1829  write(*,'(A, 2I5, 4F10.4)') 'NW CORNER: ', ic, jc, grid_global(1,npy,:,1)*90./pi
1830  ic = p_ind(npx,npy,1) ; jc = p_ind(npx,npy,1)
1831  write(*,'(A, 2I5, 4F10.4)') 'NE CORNER: ', ic, jc, grid_global(npx,npy,:,1)*90./pi
1832  ic = p_ind(npx,1,1) ; jc = p_ind(npx,1,1)
1833  write(*,'(A, 2I5, 4F10.4)') 'SE CORNER: ', ic, jc, grid_global(npx,1,:,1)*90./pi
1834  else
1835  write(*,*) 'PARENT GRID ', atm%parent_grid%grid_number, atm%parent_grid%global_tile
1836  ic = p_ind(1,1,1) ; jc = p_ind(1,1,1)
1837  write(*,'(A, 2I5, 4F10.4)') 'SW CORNER: ', ic, jc, atm%parent_grid%grid_global(ic,jc,:,parent_tile)*90./pi
1838  ic = p_ind(1,npy,1) ; jc = p_ind(1,npy,1)
1839  write(*,'(A, 2I5, 4F10.4)') 'NW CORNER: ', ic, jc, atm%parent_grid%grid_global(ic,jc,:,parent_tile)*90./pi
1840  ic = p_ind(npx,npy,1) ; jc = p_ind(npx,npy,1)
1841  write(*,'(A, 2I5, 4F10.4)') 'NE CORNER: ', ic, jc, atm%parent_grid%grid_global(ic,jc,:,parent_tile)*90./pi
1842  ic = p_ind(npx,1,1) ; jc = p_ind(npx,1,1)
1843  write(*,'(A, 2I5, 4F10.4)') 'SE CORNER: ', ic, jc, atm%parent_grid%grid_global(ic,jc,:,parent_tile)*90./pi
1844  endif
1845  end if
1846 
1847  end subroutine setup_aligned_nest
1848 
1849  subroutine setup_latlon(deglon_start,deglon_stop, deglat_start, deglat_stop, bd )
1850 
1851  type(fv_grid_bounds_type), intent(IN) :: bd
1852  real(kind=R_GRID), intent(IN) :: deglon_start,deglon_stop, deglat_start, deglat_stop
1853  real(kind=R_GRID) :: lon_start, lat_start, area_j
1854  integer :: is, ie, js, je
1855  integer :: isd, ied, jsd, jed
1856 
1857  is = bd%is
1858  ie = bd%ie
1859  js = bd%js
1860  je = bd%je
1861  isd = bd%isd
1862  ied = bd%ied
1863  jsd = bd%jsd
1864  jed = bd%jed
1865 
1866 
1867  dl = (deglon_stop-deglon_start)*pi/(180.*(npx-1))
1868  dp = (deglat_stop-deglat_start)*pi/(180.*(npy-1))
1869 
1870  lon_start = deglon_start*pi/180.
1871  lat_start = deglat_start*pi/180.
1872 
1873  do j=jsd,jed+1
1874  do i=isd,ied+1
1875  grid(i,j,1) = lon_start + real(i-1)*dl
1876  grid(i,j,2) = lat_start + real(j-1)*dp
1877  enddo
1878  enddo
1879 
1880  do j=jsd,jed
1881  do i=isd,ied
1882  agrid(i,j,1) = (grid(i,j,1) + grid(i+1,j,1))/2.
1883  agrid(i,j,2) = (grid(i,j,2) + grid(i,j+1,2))/2.
1884  enddo
1885  enddo
1886 
1887 
1888  do j=jsd,jed
1889  do i=isd,ied+1
1890  dxc(i,j) = dl*radius*cos(agrid(is,j,2))
1891  rdxc(i,j) = 1./dxc(i,j)
1892  enddo
1893  enddo
1894  do j=jsd,jed+1
1895  do i=isd,ied
1896  dyc(i,j) = dp*radius
1897  rdyc(i,j) = 1./dyc(i,j)
1898  enddo
1899  enddo
1900 
1901  do j=jsd,jed
1902  do i=isd,ied
1903  dxa(i,j) = dl*radius*cos(agrid(i,j,2))
1904  dya(i,j) = dp*radius
1905  rdxa(i,j) = 1./dxa(i,j)
1906  rdya(i,j) = 1./dya(i,j)
1907  enddo
1908  enddo
1909 
1910  do j=jsd,jed+1
1911  do i=isd,ied
1912  dx(i,j) = dl*radius*cos(grid(i,j,2))
1913  rdx(i,j) = 1./dx(i,j)
1914  enddo
1915  enddo
1916  do j=jsd,jed
1917  do i=isd,ied+1
1918  dy(i,j) = dp*radius
1919  rdy(i,j) = 1./dy(i,j)
1920  enddo
1921  enddo
1922 
1923  do j=jsd,jed
1924  area_j = radius*radius*dl*(sin(grid(is,j+1,2))-sin(grid(is,j,2)))
1925  do i=isd,ied
1926  area(i,j) = area_j
1927  rarea(i,j) = 1./area_j
1928  enddo
1929  enddo
1930 
1931  do j=jsd+1,jed
1932  area_j = radius*radius*dl*(sin(agrid(is,j,2))-sin(agrid(is,j-1,2)))
1933  do i=isd,ied+1
1934  area_c(i,j) = area_j
1935  rarea_c(i,j) = 1./area_j
1936  enddo
1937  enddo
1938  if (jsd==1) then
1939  j=1
1940  area_j = radius*radius*dl*(sin(agrid(is,j,2))-sin(agrid(is,j,2)-dp))
1941  do i=isd,ied+1
1942  area_c(i,j) = area_j
1943  rarea_c(i,j) = 1./area_j
1944  enddo
1945  endif
1946  if (jed+1==npy) then
1947  j=npy
1948  area_j = radius*radius*dl*(sin(agrid(is,j-1,2)+dp)-sin(agrid(is,j-1,2)))
1949  do i=isd,ied+1
1950  area_c(i,j) = area_j
1951  rarea_c(i,j) = 1./area_j
1952  enddo
1953  endif
1954  call mpp_update_domains( area_c, atm%domain, position=corner, complete=.true.)
1955 
1956  sina(:,:) = 1.
1957  cosa(:,:) = 0.
1958 
1959  e1(1,:,:) = 1.
1960  e1(2,:,:) = 0.
1961  e1(3,:,:) = 0.
1962 
1963  e2(1,:,:) = 0.
1964  e2(2,:,:) = 1.
1965  e2(3,:,:) = 0.
1966 
1967  end subroutine setup_latlon
1968 
1969  end subroutine init_grid
1970 
1971  subroutine cartesian_to_spherical(x, y, z, lon, lat, r)
1972  real(kind=R_GRID) , intent(IN) :: x, y, z
1973  real(kind=R_GRID) , intent(OUT) :: lon, lat, r
1974 
1975  r = sqrt(x*x + y*y + z*z)
1976  if ( (abs(x) + abs(y)) < 1.e-10 ) then ! poles:
1977  lon = 0.
1978  else
1979  lon = atan2(y,x) ! range: [-pi,pi]
1980  endif
1981 
1982 #ifdef RIGHT_HAND
1983  lat = asin(z/r)
1984 #else
1985  lat = acos(z/r) - pi/2.
1986 #endif
1987 
1988  end subroutine cartesian_to_spherical
1989  subroutine spherical_to_cartesian(lon, lat, r, x, y, z)
1990  real(kind=R_GRID) , intent(IN) :: lon, lat, r
1991  real(kind=R_GRID) , intent(OUT) :: x, y, z
1992 
1993  x = r * cos(lon) * cos(lat)
1994  y = r * sin(lon) * cos(lat)
1995 #ifdef RIGHT_HAND
1996  z = r * sin(lat)
1997 #else
1998  z = -r * sin(lat)
1999 #endif
2000  end subroutine spherical_to_cartesian
2001 
2004  subroutine rot_3d(axis, x1in, y1in, z1in, angle, x2out, y2out, z2out, degrees, convert)
2005  integer, intent(IN) :: axis
2006  real(kind=R_GRID) , intent(IN) :: x1in, y1in, z1in
2007  real(kind=R_GRID) , intent(INOUT) :: angle
2008  real(kind=R_GRID) , intent(OUT) :: x2out, y2out, z2out
2009  integer, intent(IN), optional :: degrees
2011  integer, intent(IN), optional :: convert
2014 
2015  real(kind=R_GRID) :: c, s
2016  real(kind=R_GRID) :: x1,y1,z1, x2,y2,z2
2017 
2018  if ( present(convert) ) then
2019  call spherical_to_cartesian(x1in, y1in, z1in, x1, y1, z1)
2020  else
2021  x1=x1in
2022  y1=y1in
2023  z1=z1in
2024  endif
2025 
2026  if ( present(degrees) ) then
2027  angle = angle*torad
2028  endif
2029 
2030  c = cos(angle)
2031  s = sin(angle)
2032 
2033  SELECT CASE(axis)
2034 
2035  CASE(1)
2036  x2 = x1
2037  y2 = c*y1 + s*z1
2038  z2 = -s*y1 + c*z1
2039  CASE(2)
2040  x2 = c*x1 - s*z1
2041  y2 = y1
2042  z2 = s*x1 + c*z1
2043  CASE(3)
2044  x2 = c*x1 + s*y1
2045  y2 = -s*x1 + c*y1
2046  z2 = z1
2047  CASE DEFAULT
2048  write(*,*) "Invalid axis: must be 1 for X, 2 for Y, 3 for Z."
2049 
2050  END SELECT
2051 
2052  if ( present(convert) ) then
2053  call cartesian_to_spherical(x2, y2, z2, x2out, y2out, z2out)
2054  else
2055  x2out=x2
2056  y2out=y2
2057  z2out=z2
2058  endif
2059 
2060  end subroutine rot_3d
2061 
2062 
2063 
2067  real(kind=R_GRID) function get_area_tri(ndims, p_1, p_2, p_3) &
2068  result(myarea)
2069 
2070  integer, intent(IN) :: ndims
2071  real(kind=R_GRID) , intent(IN) :: p_1(ndims) !
2072  real(kind=R_GRID) , intent(IN) :: p_2(ndims) !
2073  real(kind=R_GRID) , intent(IN) :: p_3(ndims) !
2074 
2075  real(kind=R_GRID) :: angA, angB, angC
2076 
2077  if ( ndims==3 ) then
2078  anga = spherical_angle(p_1, p_2, p_3)
2079  angb = spherical_angle(p_2, p_3, p_1)
2080  angc = spherical_angle(p_3, p_1, p_2)
2081  else
2082  anga = get_angle(ndims, p_1, p_2, p_3, 1)
2083  angb = get_angle(ndims, p_2, p_3, p_1, 1)
2084  angc = get_angle(ndims, p_3, p_1, p_2, 1)
2085  endif
2086 
2087  myarea = (anga+angb+angc - pi) * radius**2
2088 
2089  end function get_area_tri
2090 
2094  subroutine grid_area(nx, ny, ndims, nregions, bounded_domain, gridstruct, domain, bd )
2095  type(fv_grid_bounds_type), intent(IN) :: bd
2096  integer, intent(IN) :: nx, ny, ndims, nregions
2097  logical, intent(IN) :: bounded_domain
2098  type(fv_grid_type), intent(IN), target :: gridstruct
2099  type(domain2d), intent(INOUT) :: domain
2100 
2101  real(kind=R_GRID) :: p_lL(ndims)
2102  real(kind=R_GRID) :: p_uL(ndims)
2103  real(kind=R_GRID) :: p_lR(ndims)
2104  real(kind=R_GRID) :: p_uR(ndims)
2105  real(kind=R_GRID) :: a1, d1, d2, mydx, mydy, globalarea
2106 
2107  real(kind=R_GRID) :: p1(ndims), p2(ndims), p3(ndims), pi1(ndims), pi2(ndims)
2108 
2109  real(kind=R_GRID) :: maxarea, minarea
2110 
2111  integer :: i,j,n, nreg
2112  integer :: nh = 0
2113 
2114  real(kind=R_GRID), allocatable :: p_R8(:,:,:)
2115 
2116  real(kind=R_GRID), pointer, dimension(:,:,:) :: grid, agrid
2117  integer, pointer, dimension(:,:,:) :: iinta, jinta, iintb, jintb
2118  real(kind=R_GRID), pointer, dimension(:,:) :: area, area_c
2119 
2120  integer :: is, ie, js, je
2121  integer :: isd, ied, jsd, jed, ng
2122 
2123  is = bd%is
2124  ie = bd%ie
2125  js = bd%js
2126  je = bd%je
2127  isd = bd%isd
2128  ied = bd%ied
2129  jsd = bd%jsd
2130  jed = bd%jed
2131  ng = bd%ng
2132 
2133  grid => gridstruct%grid_64
2134  agrid => gridstruct%agrid_64
2135  iinta => gridstruct%iinta
2136  jinta => gridstruct%jinta
2137  iintb => gridstruct%iintb
2138  jintb => gridstruct%jintb
2139 
2140  area => gridstruct%area_64
2141  area_c => gridstruct%area_c_64
2142 
2143  if (bounded_domain) nh = ng
2144 
2145  maxarea = -1.e25
2146  minarea = 1.e25
2147 
2148  globalarea = 0.0
2149  do j=js-nh,je+nh
2150  do i=is-nh,ie+nh
2151  do n=1,ndims
2152  if ( gridstruct%stretched_grid .or. bounded_domain ) then
2153  p_ll(n) = grid(i ,j ,n)
2154  p_ul(n) = grid(i ,j+1,n)
2155  p_lr(n) = grid(i+1,j ,n)
2156  p_ur(n) = grid(i+1,j+1,n)
2157  else
2158  p_ll(n) = grid(iinta(1,i,j), jinta(1,i,j), n)
2159  p_ul(n) = grid(iinta(2,i,j), jinta(2,i,j), n)
2160  p_lr(n) = grid(iinta(4,i,j), jinta(4,i,j), n)
2161  p_ur(n) = grid(iinta(3,i,j), jinta(3,i,j), n)
2162  endif
2163  enddo
2164 
2165  ! Spherical Excess Formula
2166  area(i,j) = get_area(p_ll, p_ul, p_lr, p_ur, radius)
2167  !maxarea=MAX(area(i,j),maxarea)
2168  !minarea=MIN(area(i,j),minarea)
2169  !globalarea = globalarea + area(i,j)
2170  enddo
2171  enddo
2172 
2173 !!$ allocate( p_R8(nx-1,ny-1,ntiles_g) ) ! this is a "global" array
2174 !!$ do j=js,je
2175 !!$ do i=is,ie
2176 !!$ p_R8(i,j,tile) = area(i,j)
2177 !!$ enddo
2178 !!$ enddo
2179 !!$ call mp_gather(p_R8, is,ie, js,je, nx-1, ny-1, ntiles_g)
2180 !!$ if (is_master()) then
2181 !!$ globalarea = 0.0
2182 !!$ do n=1,ntiles_g
2183 !!$ do j=1,ny-1
2184 !!$ do i=1,nx-1
2185 !!$ globalarea = globalarea + p_R8(i,j,n)
2186 !!$ enddo
2187 !!$ enddo
2188 !!$ enddo
2189 !!$ endif
2190 !!$
2191 !!$ call mpp_broadcast(globalarea, mpp_root_pe())
2192 !!$
2193 !!$ deallocate( p_R8 )
2194 !!$
2195 !!$ call mp_reduce_max(maxarea)
2196 !!$ minarea = -minarea
2197 !!$ call mp_reduce_max(minarea)
2198 !!$ minarea = -minarea
2199 
2200  globalarea = mpp_global_sum(domain, area)
2201  maxarea = mpp_global_max(domain, area)
2202  minarea = mpp_global_min(domain, area)
2203 
2204  if (is_master()) write(*,209) 'MAX AREA (m*m):', maxarea, ' MIN AREA (m*m):', minarea
2205  if (is_master()) write(*,209) 'GLOBAL AREA (m*m):', globalarea, ' IDEAL GLOBAL AREA (m*m):', 4.0*pi*radius**2
2206  209 format(a,e21.14,a,e21.14)
2207 
2208  if (bounded_domain) then
2209  nh = ng-1 !cannot get rarea_c on boundary directly
2210  area_c = 1.e30
2211  end if
2212 
2213  do j=js-nh,je+nh+1
2214  do i=is-nh,ie+nh+1
2215  do n=1,ndims
2216  if ( gridstruct%stretched_grid .or. bounded_domain ) then
2217  p_ll(n) = agrid(i-1,j-1,n)
2218  p_lr(n) = agrid(i ,j-1,n)
2219  p_ul(n) = agrid(i-1,j ,n)
2220  p_ur(n) = agrid(i ,j ,n)
2221  else
2222  p_ll(n) = agrid(iintb(1,i,j), jintb(1,i,j), n)
2223  p_lr(n) = agrid(iintb(2,i,j), jintb(2,i,j), n)
2224  p_ul(n) = agrid(iintb(4,i,j), jintb(4,i,j), n)
2225  p_ur(n) = agrid(iintb(3,i,j), jintb(3,i,j), n)
2226  endif
2227  enddo
2228  ! Spherical Excess Formula
2229  area_c(i,j) = get_area(p_ll, p_ul, p_lr, p_ur, radius)
2230  enddo
2231  enddo
2232 
2233 ! Corners: assuming triangular cells
2234  if (gridstruct%cubed_sphere .and. .not. bounded_domain) then
2235 ! SW:
2236  i=1
2237  j=1
2238  if ( (is==1) .and. (js==1) ) then
2239  do n=1,ndims
2240  if ( gridstruct%stretched_grid ) then
2241  p1(n) = agrid(i-1,j ,n)
2242  p2(n) = agrid(i ,j ,n)
2243  p3(n) = agrid(i ,j-1,n)
2244  else
2245  p1(n) = agrid(iintb(1,i,j), jintb(1,i,j), n)
2246  p2(n) = agrid(iintb(2,i,j), jintb(2,i,j), n)
2247  p3(n) = agrid(iintb(3,i,j), jintb(3,i,j), n)
2248  endif
2249  enddo
2250  area_c(i,j) = get_area_tri(ndims, p1, p2, p3)
2251  endif
2252 
2253  i=nx
2254  j=1
2255  if ( (ie+1==nx) .and. (js==1) ) then
2256  do n=1,ndims
2257  if ( gridstruct%stretched_grid ) then
2258  p1(n) = agrid(i ,j ,n)
2259  p2(n) = agrid(i-1,j ,n)
2260  p3(n) = agrid(i-1,j-1,n)
2261  else
2262  p1(n) = agrid(iintb(1,i,j), jintb(1,i,j), n)
2263  p2(n) = agrid(iintb(2,i,j), jintb(2,i,j), n)
2264  p3(n) = agrid(iintb(3,i,j), jintb(3,i,j), n)
2265  endif
2266  enddo
2267  area_c(i,j) = get_area_tri(ndims, p1, p2, p3)
2268  endif
2269 
2270  i=nx
2271  j=ny
2272  if ( (ie+1==nx) .and. (je+1==ny) ) then
2273  do n=1,ndims
2274  if ( gridstruct%stretched_grid ) then
2275  p1(n) = agrid(i-1,j ,n)
2276  p2(n) = agrid(i-1,j-1,n)
2277  p3(n) = agrid(i ,j-1,n)
2278  else
2279  p1(n) = agrid(iintb(1,i,j), jintb(1,i,j), n)
2280  p2(n) = agrid(iintb(2,i,j), jintb(2,i,j), n)
2281  p3(n) = agrid(iintb(3,i,j), jintb(3,i,j), n)
2282  endif
2283  enddo
2284  area_c(i,j) = get_area_tri(ndims, p1, p2, p3)
2285  endif
2286 
2287  i=1
2288  j=ny
2289  if ( (is==1) .and. (je+1==ny) ) then
2290  do n=1,ndims
2291  if ( gridstruct%stretched_grid ) then
2292  p1(n) = agrid(i ,j ,n)
2293  p2(n) = agrid(i ,j-1,n)
2294  p3(n) = agrid(i-1,j-1,n)
2295  else
2296  p1(n) = agrid(iintb(1,i,j), jintb(1,i,j), n)
2297  p2(n) = agrid(iintb(2,i,j), jintb(2,i,j), n)
2298  p3(n) = agrid(iintb(3,i,j), jintb(3,i,j), n)
2299  endif
2300  enddo
2301  area_c(i,j) = get_area_tri(ndims, p1, p2, p3)
2302  endif
2303  endif
2304 
2305  end subroutine grid_area
2306 
2311  real(kind=R_GRID) function get_angle(ndims, p1, p2, p3, rad) result (angle)
2312 
2313  integer, intent(IN) :: ndims
2314  real(kind=R_GRID) , intent(IN) :: p1(ndims)
2315  real(kind=R_GRID) , intent(IN) :: p2(ndims)
2316  real(kind=R_GRID) , intent(IN) :: p3(ndims)
2317  integer, intent(in), optional:: rad
2318 
2319  real(kind=R_GRID) :: e1(3), e2(3), e3(3)
2320 
2321  if (ndims == 2) then
2322  call spherical_to_cartesian(p2(1), p2(2), real(1.,kind=R_GRID), e1(1), e1(2), e1(3))
2323  call spherical_to_cartesian(p1(1), p1(2), real(1.,kind=R_GRID), e2(1), e2(2), e2(3))
2324  call spherical_to_cartesian(p3(1), p3(2), real(1.,kind=R_GRID), e3(1), e3(2), e3(3))
2325  else
2326  e1 = p2; e2 = p1; e3 = p3
2327  endif
2328 
2329 ! High precision version:
2330  if ( present(rad) ) then
2331  angle = spherical_angle(e1, e2, e3)
2332  else
2333  angle = todeg * spherical_angle(e1, e2, e3)
2334  endif
2335 
2336  end function get_angle
2337 
2339  subroutine mirror_grid(grid_global,ng,npx,npy,ndims,nregions)
2340  integer, intent(IN) :: ng,npx,npy,ndims,nregions
2341  real(kind=R_GRID) , intent(INOUT) :: grid_global(1-ng:npx +ng,1-ng:npy +ng,ndims,1:nregions)
2342  integer :: i,j,n,n1,n2,nreg
2343  real(kind=R_GRID) :: x1,y1,z1, x2,y2,z2, ang
2344 
2345  nreg = 1
2346  do j=1,ceiling(npy/2.)
2347  do i=1,ceiling(npx/2.)
2348 
2349  x1 = 0.25d0 * (abs(grid_global(i ,j ,1,nreg)) + &
2350  abs(grid_global(npx-(i-1),j ,1,nreg)) + &
2351  abs(grid_global(i ,npy-(j-1),1,nreg)) + &
2352  abs(grid_global(npx-(i-1),npy-(j-1),1,nreg)))
2353  grid_global(i ,j ,1,nreg) = sign(x1,grid_global(i ,j ,1,nreg))
2354  grid_global(npx-(i-1),j ,1,nreg) = sign(x1,grid_global(npx-(i-1),j ,1,nreg))
2355  grid_global(i ,npy-(j-1),1,nreg) = sign(x1,grid_global(i ,npy-(j-1),1,nreg))
2356  grid_global(npx-(i-1),npy-(j-1),1,nreg) = sign(x1,grid_global(npx-(i-1),npy-(j-1),1,nreg))
2357 
2358  y1 = 0.25d0 * (abs(grid_global(i ,j ,2,nreg)) + &
2359  abs(grid_global(npx-(i-1),j ,2,nreg)) + &
2360  abs(grid_global(i ,npy-(j-1),2,nreg)) + &
2361  abs(grid_global(npx-(i-1),npy-(j-1),2,nreg)))
2362  grid_global(i ,j ,2,nreg) = sign(y1,grid_global(i ,j ,2,nreg))
2363  grid_global(npx-(i-1),j ,2,nreg) = sign(y1,grid_global(npx-(i-1),j ,2,nreg))
2364  grid_global(i ,npy-(j-1),2,nreg) = sign(y1,grid_global(i ,npy-(j-1),2,nreg))
2365  grid_global(npx-(i-1),npy-(j-1),2,nreg) = sign(y1,grid_global(npx-(i-1),npy-(j-1),2,nreg))
2366 
2367  ! force dateline/greenwich-meridion consitency
2368  if (mod(npx,2) /= 0) then
2369  if ( (i==1+(npx-1)/2.0d0) ) then
2370  grid_global(i,j ,1,nreg) = 0.0d0
2371  grid_global(i,npy-(j-1),1,nreg) = 0.0d0
2372  endif
2373  endif
2374 
2375  enddo
2376  enddo
2377 
2378  do nreg=2,nregions
2379  do j=1,npy
2380  do i=1,npx
2381 
2382  x1 = grid_global(i,j,1,1)
2383  y1 = grid_global(i,j,2,1)
2384  z1 = radius
2385 
2386  if (nreg == 2) then
2387  ang = -90.d0
2388  call rot_3d( 3, x1, y1, z1, ang, x2, y2, z2, 1, 1) ! rotate about the z-axis
2389  elseif (nreg == 3) then
2390  ang = -90.d0
2391  call rot_3d( 3, x1, y1, z1, ang, x2, y2, z2, 1, 1) ! rotate about the z-axis
2392  ang = 90.d0
2393  call rot_3d( 1, x2, y2, z2, ang, x1, y1, z1, 1, 1) ! rotate about the x-axis
2394  x2=x1
2395  y2=y1
2396  z2=z1
2397 
2398  ! force North Pole and dateline/greenwich-meridion consitency
2399  if (mod(npx,2) /= 0) then
2400  if ( (i==1+(npx-1)/2.0d0) .and. (i==j) ) then
2401  x2 = 0.0d0
2402  y2 = pi/2.0d0
2403  endif
2404  if ( (j==1+(npy-1)/2.0d0) .and. (i < 1+(npx-1)/2.0d0) ) then
2405  x2 = 0.0d0
2406  endif
2407  if ( (j==1+(npy-1)/2.0d0) .and. (i > 1+(npx-1)/2.0d0) ) then
2408  x2 = pi
2409  endif
2410  endif
2411 
2412  elseif (nreg == 4) then
2413  ang = -180.d0
2414  call rot_3d( 3, x1, y1, z1, ang, x2, y2, z2, 1, 1) ! rotate about the z-axis
2415  ang = 90.d0
2416  call rot_3d( 1, x2, y2, z2, ang, x1, y1, z1, 1, 1) ! rotate about the x-axis
2417  x2=x1
2418  y2=y1
2419  z2=z1
2420 
2421  ! force dateline/greenwich-meridion consitency
2422  if (mod(npx,2) /= 0) then
2423  if ( (j==1+(npy-1)/2.0d0) ) then
2424  x2 = pi
2425  endif
2426  endif
2427 
2428  elseif (nreg == 5) then
2429  ang = 90.d0
2430  call rot_3d( 3, x1, y1, z1, ang, x2, y2, z2, 1, 1) ! rotate about the z-axis
2431  ang = 90.d0
2432  call rot_3d( 2, x2, y2, z2, ang, x1, y1, z1, 1, 1) ! rotate about the y-axis
2433  x2=x1
2434  y2=y1
2435  z2=z1
2436  elseif (nreg == 6) then
2437  ang = 90.d0
2438  call rot_3d( 2, x1, y1, z1, ang, x2, y2, z2, 1, 1) ! rotate about the y-axis
2439  ang = 0.d0
2440  call rot_3d( 3, x2, y2, z2, ang, x1, y1, z1, 1, 1) ! rotate about the z-axis
2441  x2=x1
2442  y2=y1
2443  z2=z1
2444 
2445  ! force South Pole and dateline/greenwich-meridion consitency
2446  if (mod(npx,2) /= 0) then
2447  if ( (i==1+(npx-1)/2.0d0) .and. (i==j) ) then
2448  x2 = 0.0d0
2449  y2 = -pi/2.0d0
2450  endif
2451  if ( (i==1+(npx-1)/2.0d0) .and. (j > 1+(npy-1)/2.0d0) ) then
2452  x2 = 0.0d0
2453  endif
2454  if ( (i==1+(npx-1)/2.0d0) .and. (j < 1+(npy-1)/2.0d0) ) then
2455  x2 = pi
2456  endif
2457  endif
2458 
2459  endif
2460 
2461  grid_global(i,j,1,nreg) = x2
2462  grid_global(i,j,2,nreg) = y2
2463 
2464  enddo
2465  enddo
2466  enddo
2467 
2468  end subroutine mirror_grid
2469 
2470  end module fv_grid_tools_mod
2471 
fv_grid_utils_mod::mid_pt_sphere
subroutine, public mid_pt_sphere(p1, p2, pm)
Definition: fv_grid_utils.F90:2020
fv_grid_tools_mod::torad
real(kind=r_grid), parameter torad
convert to radians
Definition: fv_grid_tools.F90:161
fv_mp_mod
The module &#39;fv_mp_mod&#39; is a single program multiple data (SPMD) parallel decompostion/communication m...
Definition: fv_mp_mod.F90:24
fv_timing_mod::timing_off
subroutine timing_off(blk_name)
The subroutine &#39;timing_off&#39; stops a timer.
Definition: fv_timing.F90:180
fv_arrays_mod::fv_grid_type
The type &#39;fv_grid_type&#39; is made up of grid-dependent information from fv_grid_tools and fv_grid_utils...
Definition: fv_arrays.F90:123
fv_arrays_mod::fv_grid_bounds_type
Definition: fv_arrays.F90:1187
sorted_index_mod
The module &#39;sorted_index&#39; sorts cell corner indices in lat-lon space to ensure the same order of oper...
Definition: sorted_index.F90:28
fv_grid_tools_mod::missing
real(kind=r_grid), parameter, public missing
Definition: fv_grid_tools.F90:162
sorted_index_mod::sorted_intb
subroutine, public sorted_intb(isd, ied, jsd, jed, is, ie, js, je, npx, npy, cubed_sphere, agrid, iintb, jintb)
The subroutine &#39;sorted_intb&#39; sorts cell corner indices in latlon space based on grid locations in ind...
Definition: sorted_index.F90:235
setup_latlon
subroutine setup_latlon(deglon_start, deglon_stop, deglat_start, deglat_stop, bd)
Definition: fv_grid_tools.F90:1850
fv_grid_tools_mod::radius
real(kind=r_grid), parameter radius
Definition: fv_grid_tools.F90:158
fv_grid_utils_mod::direct_transform
subroutine, public direct_transform(c, i1, i2, j1, j2, lon_p, lat_p, n, lon, lat)
The subroutine &#39;direct_transform&#39; performs a direct transformation of the standard (symmetrical) cubi...
Definition: fv_grid_utils.F90:920
fv_grid_tools_mod::get_area_tri
real(kind=r_grid) function get_area_tri(ndims, p_1, p_2, p_3)
brief The function &#39;get_area_tri&#39; gets the surface area of a cell defined as a triangle on the sphere...
Definition: fv_grid_tools.F90:2069
fv_grid_utils_mod::dist2side_latlon
real(kind=r_grid) function, public dist2side_latlon(v1, v2, point)
The function &#39;dist2side_latlon&#39; is the version of &#39;dist2side&#39; that takes points in latitude-longitude...
Definition: fv_grid_utils.F90:2836
fv_grid_utils_mod::spherical_angle
real(kind=r_grid) function, public spherical_angle(p1, p2, p3)
Definition: fv_grid_utils.F90:2861
fv_grid_utils_mod::inner_prod
real function, public inner_prod(v1, v2)
Definition: fv_grid_utils.F90:1039
fv_grid_utils_mod::spherical_linear_interpolation
subroutine, public spherical_linear_interpolation(beta, p1, p2, pb)
The subroutine &#39;spherical_linear_interpolation&#39; interpolates along the great circle connecting points...
Definition: fv_grid_utils.F90:1968
setup_aligned_nest
subroutine setup_aligned_nest(Atm)
Definition: fv_grid_tools.F90:1320
fv_grid_tools_mod::rot_3d
subroutine rot_3d(axis, x1in, y1in, z1in, angle, x2out, y2out, z2out, degrees, convert)
The subroutine &#39;rot_3d&#39; rotates points on a sphere in xyz coordinates.
Definition: fv_grid_tools.F90:2005
fv_grid_tools_mod::todeg
real(kind=r_grid), parameter, public todeg
convert to degrees
Definition: fv_grid_tools.F90:160
fv_arrays_mod::r_grid
integer, parameter, public r_grid
Definition: fv_arrays.F90:34
fv_timing_mod
The module &#39;fv_timing&#39; contains FV3 timers.
Definition: fv_timing.F90:24
fv_arrays_mod
The module &#39;fv_arrays&#39; contains the &#39;fv_atmos_type&#39; and associated datatypes.
Definition: fv_arrays.F90:24
fv_grid_utils_mod::big_number
real, parameter, public big_number
Definition: fv_grid_utils.F90:104
fv_grid_utils_mod::great_circle_dist
real function, public great_circle_dist(q1, q2, radius)
Definition: fv_grid_utils.F90:2079
fv_grid_utils_mod::get_area
real(kind=r_grid) function, public get_area(p1, p4, p2, p3, radius)
Definition: fv_grid_utils.F90:2774
fv_grid_utils_mod::cell_center2
subroutine, public cell_center2(q1, q2, q3, q4, e2)
Definition: fv_grid_utils.F90:2726
fv_grid_tools_mod::grid_area
subroutine grid_area(nx, ny, ndims, nregions, bounded_domain, gridstruct, domain, bd)
The subroutine &#39;grid_area&#39; gets the surface area on a grid in lat/lon or xyz coordinates.
Definition: fv_grid_tools.F90:2095
fv_grid_tools_mod::csfac
real(kind=r_grid) csfac
Definition: fv_grid_tools.F90:164
fv_grid_tools_mod::get_symmetry
subroutine get_symmetry(data_in, data_out, ishift, jshift, npes_x, npes_y, domain, tile, npx_g, bd)
Definition: fv_grid_tools.F90:307
fv_grid_tools_mod::read_grid
subroutine read_grid(Atm, grid_file, ndims, nregions, ng)
The subroutine &#39;read_grid&#39; reads the grid from the mosaic grid file.
Definition: fv_grid_tools.F90:176
fv_grid_tools_mod::init_grid
subroutine, public init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, ng, tile_coarse)
The subroutine &#39;init_grid&#39; reads the grid from the input file and sets up grid descriptors.
Definition: fv_grid_tools.F90:530
fv_grid_tools_mod::cartesian_to_spherical
subroutine cartesian_to_spherical(x, y, z, lon, lat, r)
Definition: fv_grid_tools.F90:1972
fv_timing_mod::timing_on
subroutine timing_on(blk_name)
The subroutine &#39;timing_on&#39; starts a timer.
Definition: fv_timing.F90:116
fv_grid_utils_mod::fill_ghost
Definition: fv_grid_utils.F90:125
fv_grid_tools_mod::spherical_to_cartesian
subroutine, public spherical_to_cartesian(lon, lat, r, x, y, z)
Definition: fv_grid_tools.F90:1990
fv_grid_utils_mod
The module &#39;fv_grid_utils&#39; contains routines for setting up and computing grid-related quantities...
Definition: fv_grid_utils.F90:26
fv_arrays_mod::fv_atmos_type
Definition: fv_arrays.F90:1216
sorted_index_mod::sorted_inta
subroutine, public sorted_inta(isd, ied, jsd, jed, cubed_sphere, bgrid, iinta, jinta)
The subroutine &#39;sorted_inta&#39; sorts cell corner indices in latlon space based on grid locations in ind...
Definition: sorted_index.F90:47
fv_grid_utils_mod::cube_transform
subroutine, public cube_transform(c, i1, i2, j1, j2, lon_p, lat_p, n, lon, lat)
Definition: fv_grid_utils.F90:975
fv_grid_tools_mod
Definition: fv_grid_tools.F90:21
fv_grid_tools_mod::write_grid_char_file
logical write_grid_char_file
Definition: fv_grid_tools.F90:167
fv_grid_tools_mod::debug_message_size
logical, parameter debug_message_size
Definition: fv_grid_tools.F90:166
fv_grid_tools_mod::mirror_grid
subroutine mirror_grid(grid_global, ng, npx, npy, ndims, nregions)
The subroutine &#39;mirror_grid&#39; mirrors the grid across the 0-longitude line.
Definition: fv_grid_tools.F90:2340
fv_grid_tools_mod::get_angle
real(kind=r_grid) function get_angle(ndims, p1, p2, p3, rad)
The function &#39;get_angle&#39; gets the angle between 3 points on a sphere in lat/lon or xyz coordinates...
Definition: fv_grid_tools.F90:2312
fv_grid_utils_mod::gnomonic_grids
subroutine, public gnomonic_grids(grid_type, im, lon, lat)
Definition: fv_grid_utils.F90:1340
setup_cartesian
subroutine setup_cartesian(npx, npy, dx_const, dy_const, deglat, bd)
Definition: fv_grid_tools.F90:1237