FV3DYCORE  Version1.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, dist2side_latlon, &
124  spherical_linear_interpolation, big_number
125  use fv_timing_mod, only: timing_on, timing_off
126  use fv_mp_mod, only: ng, is_master, fill_corners, xdir, ydir
127  use fv_mp_mod, only: mp_gather, mp_bcst, mp_reduce_max, mp_stop
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  implicit none
154  private
155 #include <netcdf.inc>
156 
157  real(kind=R_GRID), parameter:: radius = cnst_radius
158 
159  real(kind=R_GRID) , parameter:: todeg = 180.0d0/pi
160  real(kind=R_GRID) , parameter:: torad = pi/180.0d0
161  real(kind=R_GRID) , parameter:: missing = 1.d25
162 
163  real(kind=R_GRID) :: csfac
164 
165  logical, parameter :: debug_message_size = .false.
166  logical :: write_grid_char_file = .false.
167 
168 
169  public :: todeg, missing, init_grid, spherical_to_cartesian
170 
171 contains
172 
174  subroutine read_grid(Atm, grid_file, ndims, nregions, ng)
175  type(fv_atmos_type), intent(inout), target :: Atm
176  character(len=*), intent(IN) :: grid_file
177  integer, intent(IN) :: ndims
178  integer, intent(IN) :: nregions
179  integer, intent(IN) :: ng
180 
181  real, allocatable, dimension(:,:) :: tmpx, tmpy
182  real(kind=R_GRID), pointer, dimension(:,:,:) :: grid
183  character(len=128) :: units = ""
184  character(len=256) :: atm_mosaic, atm_hgrid, grid_form
185  character(len=1024) :: attvalue
186  integer :: ntiles, i, j, stdunit
187  integer :: isc2, iec2, jsc2, jec2
188  integer :: start(4), nread(4)
189  integer :: is, ie, js, je
190  integer :: isd, ied, jsd, jed
191  integer,save :: halo=3
192 
193  is = atm%bd%is
194  ie = atm%bd%ie
195  js = atm%bd%js
196  je = atm%bd%je
197  isd = atm%bd%isd
198  ied = atm%bd%ied
199  jsd = atm%bd%jsd
200  jed = atm%bd%jed
201  grid => atm%gridstruct%grid_64
202 
203  if(.not. file_exist(grid_file)) call mpp_error(fatal, 'fv_grid_tools(read_grid): file '// &
204  trim(grid_file)//' does not exist')
205 
206  !--- make sure the grid file is mosaic file.
207  if( field_exist(grid_file, 'atm_mosaic_file') .OR. field_exist(grid_file, 'gridfiles') ) then
208  stdunit = stdout()
209  write(stdunit,*) '==>Note from fv_grid_tools_mod(read_grid): read atmosphere grid from mosaic version grid'
210  else
211  call mpp_error(fatal, 'fv_grid_tools(read_grid): neither atm_mosaic_file nor gridfiles exists in file ' &
212  //trim(grid_file))
213  endif
214 
215  if(field_exist(grid_file, 'atm_mosaic_file')) then
216  call read_data(grid_file, "atm_mosaic_file", atm_mosaic)
217  atm_mosaic = "INPUT/"//trim(atm_mosaic)
218  else
219  atm_mosaic = trim(grid_file)
220  endif
221 
222  call get_mosaic_tile_grid(atm_hgrid, atm_mosaic, atm%domain)
223 
224  grid_form = "none"
225  if( get_global_att_value(atm_hgrid, "history", attvalue) ) then
226  if( index(attvalue, "gnomonic_ed") > 0) grid_form = "gnomonic_ed"
227  endif
228  if(grid_form .NE. "gnomonic_ed") call mpp_error(fatal, &
229  "fv_grid_tools(read_grid): the grid should be 'gnomonic_ed' when reading from grid file, contact developer")
230 
231  !FIXME: Doesn't work for a nested grid
232  ntiles = get_mosaic_ntiles(atm_mosaic)
233 
234  if( .not. atm%flagstruct%regional) 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%flagstruct%regional) 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%flagstruct%regional) 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)
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 !--------------------------------------------------------
539  real(kind=R_GRID) :: xs(npx,npy)
540  real(kind=R_GRID) :: ys(npx,npy)
541 
542  real(kind=R_GRID) :: dp, dl
543  real(kind=R_GRID) :: x1,x2,y1,y2,z1,z2
544  integer :: i,j,k,n,nreg
545  integer :: filelun
546 
547  real(kind=R_GRID) :: p1(3), p2(3), p3(3), p4(3)
548  real(kind=R_GRID) :: dist,dist1,dist2, pa(2), pa1(2), pa2(2), pb(2)
549  real(kind=R_GRID) :: pt(3), pt1(3), pt2(3), pt3(3)
550  real(kind=R_GRID) :: ee1(3), ee2(3)
551 
552  real(kind=R_GRID) :: angn,angm,angav,ang
553  real(kind=R_GRID) :: aspn,aspm,aspav,asp
554  real(kind=R_GRID) :: dxn, dxm, dxav
555  real(kind=R_GRID) :: dx_local, dy_local
556 
557  real(kind=R_GRID) :: vec1(3), vec2(3), vec3(3), vec4(3)
558  real(kind=R_GRID) :: vecavg(3), vec3a(3), vec3b(3), vec4a(3), vec4b(3)
559  real(kind=R_GRID) :: xyz1(3), xyz2(3)
560 
561 ! real(kind=R_GRID) :: grid_global(1-ng:npx +ng,1-ng:npy +ng,ndims,1:nregions)
562  integer :: ios, ip, jp
563 
564  integer :: igrid
565 
566  integer :: tmplun
567  character(len=80) :: tmpfile
568 
569  real(kind=R_GRID), dimension(Atm%bd%is:Atm%bd%ie) :: sbuffer, nbuffer
570  real(kind=R_GRID), dimension(Atm%bd%js:Atm%bd%je) :: wbuffer, ebuffer
571 
572  real(kind=R_GRID), pointer, dimension(:,:,:) :: agrid, grid
573  real(kind=R_GRID), pointer, dimension(:,:) :: area, area_c
574 
575  real(kind=R_GRID), pointer, dimension(:,:) :: sina, cosa, dx, dy, dxc, dyc, dxa, dya
576  real, pointer, dimension(:,:,:) :: e1, e2
577 
578  real, pointer, dimension(:,:) :: rarea, rarea_c
579  real, pointer, dimension(:,:) :: rdx, rdy, rdxc, rdyc, rdxa, rdya
580 
581  integer, pointer, dimension(:,:,:) :: iinta, jinta, iintb, jintb
582 
583  real(kind=R_GRID), pointer, dimension(:,:,:,:) :: grid_global
584 
585  integer, pointer :: npx_g, npy_g, ntiles_g, tile
586  logical, pointer :: sw_corner, se_corner, ne_corner, nw_corner
587  logical, pointer :: latlon, cubed_sphere, have_south_pole, have_north_pole, stretched_grid
588 
589  type(domain2d), pointer :: domain
590  integer :: is, ie, js, je
591  integer :: isd, ied, jsd, jed
592  integer :: istart, iend, jstart, jend
593 
594  is = atm%bd%is
595  ie = atm%bd%ie
596  js = atm%bd%js
597  je = atm%bd%je
598  isd = atm%bd%isd
599  ied = atm%bd%ied
600  jsd = atm%bd%jsd
601  jed = atm%bd%jed
602 
603  !!! Associate pointers
604  agrid => atm%gridstruct%agrid_64
605  grid => atm%gridstruct%grid_64
606 
607  area => atm%gridstruct%area_64
608  area_c => atm%gridstruct%area_c_64
609  rarea => atm%gridstruct%rarea
610  rarea_c => atm%gridstruct%rarea_c
611 
612  sina => atm%gridstruct%sina_64
613  cosa => atm%gridstruct%cosa_64
614  dx => atm%gridstruct%dx_64
615  dy => atm%gridstruct%dy_64
616  dxc => atm%gridstruct%dxc_64
617  dyc => atm%gridstruct%dyc_64
618  dxa => atm%gridstruct%dxa_64
619  dya => atm%gridstruct%dya_64
620  rdx => atm%gridstruct%rdx
621  rdy => atm%gridstruct%rdy
622  rdxc => atm%gridstruct%rdxc
623  rdyc => atm%gridstruct%rdyc
624  rdxa => atm%gridstruct%rdxa
625  rdya => atm%gridstruct%rdya
626  e1 => atm%gridstruct%e1
627  e2 => atm%gridstruct%e2
628 
629  if (atm%neststruct%nested .or. any(atm%neststruct%child_grids)) then
630  grid_global => atm%grid_global
631  else if( trim(grid_file) .NE. 'INPUT/grid_spec.nc') then
632  allocate(grid_global(1-ng:npx +ng,1-ng:npy +ng,ndims,1:nregions))
633  endif
634 
635  iinta => atm%gridstruct%iinta
636  jinta => atm%gridstruct%jinta
637  iintb => atm%gridstruct%iintb
638  jintb => atm%gridstruct%jintb
639  npx_g => atm%gridstruct%npx_g
640  npy_g => atm%gridstruct%npy_g
641  ntiles_g => atm%gridstruct%ntiles_g
642  sw_corner => atm%gridstruct%sw_corner
643  se_corner => atm%gridstruct%se_corner
644  ne_corner => atm%gridstruct%ne_corner
645  nw_corner => atm%gridstruct%nw_corner
646  latlon => atm%gridstruct%latlon
647  cubed_sphere => atm%gridstruct%cubed_sphere
648  have_south_pole => atm%gridstruct%have_south_pole
649  have_north_pole => atm%gridstruct%have_north_pole
650  stretched_grid => atm%gridstruct%stretched_grid
651 
652  tile => atm%tile
653 
654  domain => atm%domain
655 
656  npx_g = npx
657  npy_g = npy
658  ntiles_g = nregions
659  latlon = .false.
660  cubed_sphere = .false.
661 
662  if ( atm%flagstruct%do_schmidt .and. abs(atm%flagstruct%stretch_fac-1.) > 1.e-5 ) stretched_grid = .true.
663 
664  if (atm%flagstruct%grid_type>3) then
665  if (atm%flagstruct%grid_type == 4) then
666  call setup_cartesian(npx, npy, atm%flagstruct%dx_const, atm%flagstruct%dy_const, &
667  atm%flagstruct%deglat, atm%bd)
668  else
669  call mpp_error(fatal, 'init_grid: unsupported grid type')
670  endif
671  else
672 
673  cubed_sphere = .true.
674  if (atm%neststruct%nested) then
675  call setup_aligned_nest(atm)
676  else
677  if(trim(grid_file) == 'INPUT/grid_spec.nc') then
678  call read_grid(atm, grid_file, ndims, nregions, ng)
679  else
680  if (atm%flagstruct%grid_type>=0) call gnomonic_grids(atm%flagstruct%grid_type, npx-1, xs, ys)
681  if (is_master()) then
682  if (atm%flagstruct%grid_type>=0) then
683  do j=1,npy
684  do i=1,npx
685  grid_global(i,j,1,1) = xs(i,j)
686  grid_global(i,j,2,1) = ys(i,j)
687  enddo
688  enddo
689 ! mirror_grid assumes that the tile=1 is centered on equator and greenwich meridian Lon[-pi,pi]
690  call mirror_grid(grid_global, ng, npx, npy, 2, 6)
691  do n=1,nregions
692  do j=1,npy
693  do i=1,npx
694 !---------------------------------
695 ! Shift the corner away from Japan
696 !---------------------------------
697 !--------------------- This will result in the corner close to east coast of China ------------------
698  if ( .not.atm%flagstruct%do_schmidt .and. (atm%flagstruct%shift_fac)>1.e-4 ) &
699  grid_global(i,j,1,n) = grid_global(i,j,1,n) - pi/atm%flagstruct%shift_fac
700 !----------------------------------------------------------------------------------------------------
701  if ( grid_global(i,j,1,n) < 0. ) &
702  grid_global(i,j,1,n) = grid_global(i,j,1,n) + 2.*pi
703  if (abs(grid_global(i,j,1,1)) < 1.d-10) grid_global(i,j,1,1) = 0.0
704  if (abs(grid_global(i,j,2,1)) < 1.d-10) grid_global(i,j,2,1) = 0.0
705  enddo
706  enddo
707  enddo
708  else
709  call mpp_error(fatal, "fv_grid_tools: reading of ASCII grid files no longer supported")
710  endif ! Atm%flagstruct%grid_type>=0
711 
712  grid_global( 1,1:npy,:,2)=grid_global(npx,1:npy,:,1)
713  grid_global( 1,1:npy,:,3)=grid_global(npx:1:-1,npy,:,1)
714  grid_global(1:npx,npy,:,5)=grid_global(1,npy:1:-1,:,1)
715  grid_global(1:npx,npy,:,6)=grid_global(1:npx,1,:,1)
716 
717  grid_global(1:npx, 1,:,3)=grid_global(1:npx,npy,:,2)
718  grid_global(1:npx, 1,:,4)=grid_global(npx,npy:1:-1,:,2)
719  grid_global(npx,1:npy,:,6)=grid_global(npx:1:-1,1,:,2)
720 
721  grid_global( 1,1:npy,:,4)=grid_global(npx,1:npy,:,3)
722  grid_global( 1,1:npy,:,5)=grid_global(npx:1:-1,npy,:,3)
723 
724  grid_global(npx,1:npy,:,3)=grid_global(1,1:npy,:,4)
725  grid_global(1:npx, 1,:,5)=grid_global(1:npx,npy,:,4)
726  grid_global(1:npx, 1,:,6)=grid_global(npx,npy:1:-1,:,4)
727 
728  grid_global( 1,1:npy,:,6)=grid_global(npx,1:npy,:,5)
729 
730 !------------------------
731 ! Schmidt transformation:
732 !------------------------
733  if ( atm%flagstruct%do_schmidt ) then
734  do n=1,nregions
735  call direct_transform(atm%flagstruct%stretch_fac, 1, npx, 1, npy, &
736  atm%flagstruct%target_lon, atm%flagstruct%target_lat, &
737  n, grid_global(1:npx,1:npy,1,n), grid_global(1:npx,1:npy,2,n))
738  enddo
739  endif
740  endif !is master
741  call mpp_broadcast(grid_global, size(grid_global), mpp_root_pe())
742 !--- copy grid to compute domain
743  do n=1,ndims
744  do j=js,je+1
745  do i=is,ie+1
746  grid(i,j,n) = grid_global(i,j,n,tile)
747  enddo
748  enddo
749  enddo
750  endif !(trim(grid_file) == 'INPUT/grid_spec.nc')
751 
752 !
753 ! SJL: For phys/exchange grid, etc
754 !
755  call mpp_update_domains( grid, atm%domain, position=corner)
756  if (.not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
757  call fill_corners(grid(:,:,1), npx, npy, fill=xdir, bgrid=.true.)
758  call fill_corners(grid(:,:,2), npx, npy, fill=xdir, bgrid=.true.)
759  endif
760 
761 
762  !--- dx and dy
763 
764  if( .not. atm%flagstruct%regional) then
765  istart=is
766  iend=ie
767  jstart=js
768  jend=je
769  else
770  istart=isd
771  iend=ied
772  jstart=jsd
773  jend=jed
774  endif
775 
776  do j = jstart, jend+1
777  do i = istart, iend
778  p1(1) = grid(i ,j,1)
779  p1(2) = grid(i ,j,2)
780  p2(1) = grid(i+1,j,1)
781  p2(2) = grid(i+1,j,2)
782  dx(i,j) = great_circle_dist( p2, p1, radius )
783  enddo
784  enddo
785  if( stretched_grid ) then
786  do j = jstart, jend
787  do i = istart, iend+1
788  p1(1) = grid(i,j, 1)
789  p1(2) = grid(i,j, 2)
790  p2(1) = grid(i,j+1,1)
791  p2(2) = grid(i,j+1,2)
792  dy(i,j) = great_circle_dist( p2, p1, radius )
793  enddo
794  enddo
795  else
796  call get_symmetry(dx(is:ie,js:je+1), dy(is:ie+1,js:je), 0, 1, atm%layout(1), atm%layout(2), &
797  atm%domain, atm%tile, atm%gridstruct%npx_g, atm%bd)
798  endif
799 
800  call mpp_get_boundary( dy, dx, atm%domain, ebufferx=ebuffer, wbufferx=wbuffer, sbuffery=sbuffer, nbuffery=nbuffer,&
801  flags=scalar_pair+xupdate, gridtype=cgrid_ne_param)
802  if( .not. atm%flagstruct%regional ) then
803  if(is == 1 .AND. mod(tile,2) .NE. 0) then ! on the west boundary
804  dy(is, js:je) = wbuffer(js:je)
805  endif
806  if(ie == npx-1) then ! on the east boundary
807  dy(ie+1, js:je) = ebuffer(js:je)
808  endif
809  endif
810 
811  call mpp_update_domains( dy, dx, atm%domain, flags=scalar_pair, &
812  gridtype=cgrid_ne_param, complete=.true.)
813  if (cubed_sphere .and. (.not. (atm%neststruct%nested .or. atm%flagstruct%regional))) then
814  call fill_corners(dx, dy, npx, npy, dgrid=.true.)
815  endif
816 
817  if( .not. stretched_grid ) &
818  call sorted_inta(isd, ied, jsd, jed, cubed_sphere, grid, iinta, jinta)
819 
820  agrid(:,:,:) = -1.e25
821 
822  !--- compute agrid (use same indices as for dx/dy above)
823 
824  do j=jstart,jend
825  do i=istart,iend
826  if ( stretched_grid ) then
827  call cell_center2(grid(i,j, 1:2), grid(i+1,j, 1:2), &
828  grid(i,j+1,1:2), grid(i+1,j+1,1:2), &
829  agrid(i,j,1:2) )
830  else
831  call cell_center2(grid(iinta(1,i,j),jinta(1,i,j),1:2), &
832  grid(iinta(2,i,j),jinta(2,i,j),1:2), &
833  grid(iinta(3,i,j),jinta(3,i,j),1:2), &
834  grid(iinta(4,i,j),jinta(4,i,j),1:2), &
835  agrid(i,j,1:2) )
836  endif
837  enddo
838  enddo
839 
840  call mpp_update_domains( agrid, atm%domain, position=center, complete=.true. )
841  if (.not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
842  call fill_corners(agrid(:,:,1), npx, npy, xdir, agrid=.true.)
843  call fill_corners(agrid(:,:,2), npx, npy, ydir, agrid=.true.)
844  endif
845 
846  do j=jsd,jed
847  do i=isd,ied
848  call mid_pt_sphere(grid(i, j,1:2), grid(i, j+1,1:2), p1)
849  call mid_pt_sphere(grid(i+1,j,1:2), grid(i+1,j+1,1:2), p2)
850  dxa(i,j) = great_circle_dist( p2, p1, radius )
851 !
852  call mid_pt_sphere(grid(i,j ,1:2), grid(i+1,j ,1:2), p1)
853  call mid_pt_sphere(grid(i,j+1,1:2), grid(i+1,j+1,1:2), p2)
854  dya(i,j) = great_circle_dist( p2, p1, radius )
855  enddo
856  enddo
857 ! call mpp_update_domains( dxa, dya, Atm%domain, flags=SCALAR_PAIR, gridtype=AGRID_PARAM)
858  if (cubed_sphere .and. (.not. (atm%neststruct%nested .or. atm%flagstruct%regional))) then
859  call fill_corners(dxa, dya, npx, npy, agrid=.true.)
860  endif
861 
862 
863  end if !if nested
864 
865  do j=jsd,jed
866  do i=isd+1,ied
867  dxc(i,j) = great_circle_dist(agrid(i,j,:), agrid(i-1,j,:), radius)
868  enddo
869 !xxxxxx
870  !Are the following 2 lines appropriate for the regional domain?
871 !xxxxxx
872  dxc(isd,j) = dxc(isd+1,j)
873  dxc(ied+1,j) = dxc(ied,j)
874  enddo
875 
876 ! do j=js,je+1
877 ! do i=is,ie
878  do j=jsd+1,jed
879  do i=isd,ied
880  dyc(i,j) = great_circle_dist(agrid(i,j,:), agrid(i,j-1,:), radius)
881  enddo
882  enddo
883 !xxxxxx
884  !Are the following 2 lines appropriate for the regional domain?
885 !xxxxxx
886  do i=isd,ied
887  dyc(i,jsd) = dyc(i,jsd+1)
888  dyc(i,jed+1) = dyc(i,jed)
889  end do
890 
891 
892  if( .not. stretched_grid ) &
893  call sorted_intb(isd, ied, jsd, jed, is, ie, js, je, npx, npy, &
894  cubed_sphere, agrid, iintb, jintb)
895 
896  call grid_area( npx, npy, ndims, nregions, atm%neststruct%nested, atm%gridstruct, atm%domain, atm%bd, atm%flagstruct%regional )
897 ! stretched_grid = .false.
898 
899 !----------------------------------
900 ! Compute area_c, rarea_c, dxc, dyc
901 !----------------------------------
902 
903  if ( .not. stretched_grid .and. (.not. (atm%neststruct%nested .or. atm%flagstruct%regional))) then
904 ! For symmetrical grids:
905  if ( is==1 ) then
906  i = 1
907  do j=js,je+1
908  call mid_pt_sphere(grid(i,j-1,1:2), grid(i,j, 1:2), p1)
909  call mid_pt_sphere(grid(i,j ,1:2), grid(i,j+1,1:2), p4)
910  p2(1:2) = agrid(i,j-1,1:2)
911  p3(1:2) = agrid(i,j, 1:2)
912  area_c(i,j) = 2.*get_area(p1, p4, p2, p3, radius)
913  enddo
914  do j=js,je
915  call mid_pt_sphere(grid(i,j,1:2), grid(i,j+1,1:2), p1)
916  p2(1:2) = agrid(i,j,1:2)
917  dxc(i,j) = 2.*great_circle_dist( p1, p2, radius )
918  enddo
919  endif
920  if ( (ie+1)==npx ) then
921  i = npx
922  do j=js,je+1
923  p1(1:2) = agrid(i-1,j-1,1:2)
924  call mid_pt_sphere(grid(i,j-1,1:2), grid(i,j, 1:2), p2)
925  call mid_pt_sphere(grid(i,j ,1:2), grid(i,j+1,1:2), p3)
926  p4(1:2) = agrid(i-1,j,1:2)
927  area_c(i,j) = 2.*get_area(p1, p4, p2, p3, radius)
928  enddo
929  do j=js,je
930  p1(1:2) = agrid(i-1,j,1:2)
931  call mid_pt_sphere(grid(i,j,1:2), grid(i,j+1,1:2), p2)
932  dxc(i,j) = 2.*great_circle_dist( p1, p2, radius )
933  enddo
934  endif
935  if ( js==1 ) then
936  j = 1
937  do i=is,ie+1
938  call mid_pt_sphere(grid(i-1,j,1:2), grid(i, j,1:2), p1)
939  call mid_pt_sphere(grid(i, j,1:2), grid(i+1,j,1:2), p2)
940  p3(1:2) = agrid(i, j,1:2)
941  p4(1:2) = agrid(i-1,j,1:2)
942  area_c(i,j) = 2.*get_area(p1, p4, p2, p3, radius)
943  enddo
944  do i=is,ie
945  call mid_pt_sphere(grid(i,j,1:2), grid(i+1,j,1:2), p1)
946  p2(1:2) = agrid(i,j,1:2)
947  dyc(i,j) = 2.*great_circle_dist( p1, p2, radius )
948  enddo
949  endif
950  if ( (je+1)==npy ) then
951  j = npy
952  do i=is,ie+1
953  p1(1:2) = agrid(i-1,j-1,1:2)
954  p2(1:2) = agrid(i ,j-1,1:2)
955  call mid_pt_sphere(grid(i, j,1:2), grid(i+1,j,1:2), p3)
956  call mid_pt_sphere(grid(i-1,j,1:2), grid(i, j,1:2), p4)
957  area_c(i,j) = 2.*get_area(p1, p4, p2, p3, radius)
958  enddo
959  do i=is,ie
960  p1(1:2) = agrid(i,j-1,1:2)
961  call mid_pt_sphere(grid(i,j,1:2), grid(i+1,j,1:2), p2)
962  dyc(i,j) = 2.*great_circle_dist( p1, p2, radius )
963  enddo
964  endif
965 
966  if ( sw_corner ) then
967  i=1; j=1
968  p1(1:2) = grid(i,j,1:2)
969  call mid_pt_sphere(grid(i,j,1:2), grid(i+1,j,1:2), p2)
970  p3(1:2) = agrid(i,j,1:2)
971  call mid_pt_sphere(grid(i,j,1:2), grid(i,j+1,1:2), p4)
972  area_c(i,j) = 3.*get_area(p1, p4, p2, p3, radius)
973  endif
974  if ( se_corner ) then
975  i=npx; j=1
976  call mid_pt_sphere(grid(i-1,j,1:2), grid(i,j,1:2), p1)
977  p2(1:2) = grid(i,j,1:2)
978  call mid_pt_sphere(grid(i,j,1:2), grid(i,j+1,1:2), p3)
979  p4(1:2) = agrid(i,j,1:2)
980  area_c(i,j) = 3.*get_area(p1, p4, p2, p3, radius)
981  endif
982  if ( ne_corner ) then
983  i=npx; j=npy
984  p1(1:2) = agrid(i-1,j-1,1:2)
985  call mid_pt_sphere(grid(i,j-1,1:2), grid(i,j,1:2), p2)
986  p3(1:2) = grid(i,j,1:2)
987  call mid_pt_sphere(grid(i-1,j,1:2), grid(i,j,1:2), p4)
988  area_c(i,j) = 3.*get_area(p1, p4, p2, p3, radius)
989  endif
990  if ( nw_corner ) then
991  i=1; j=npy
992  call mid_pt_sphere(grid(i,j-1,1:2), grid(i,j,1:2), p1)
993  p2(1:2) = agrid(i,j-1,1:2)
994  call mid_pt_sphere(grid(i,j,1:2), grid(i+1,j,1:2), p3)
995  p4(1:2) = grid(i,j,1:2)
996  area_c(i,j) = 3.*get_area(p1, p4, p2, p3, radius)
997  endif
998  endif
999 !-----------------
1000 
1001  call mpp_update_domains( dxc, dyc, atm%domain, flags=scalar_pair, &
1002  gridtype=cgrid_ne_param, complete=.true.)
1003  if (cubed_sphere .and. (.not. (atm%neststruct%nested .or. atm%flagstruct%regional))) then
1004  call fill_corners(dxc, dyc, npx, npy, cgrid=.true.)
1005  endif
1006 
1007  call mpp_update_domains( area, atm%domain, complete=.true. )
1008 
1009 
1010  !Handling outermost ends for area_c
1011  if (atm%neststruct%nested .or. atm%flagstruct%regional) then
1012  if (is == 1) then
1013  do j=jsd,jed
1014  area_c(isd,j) = area_c(isd+1,j)
1015  end do
1016  if (js == 1) area_c(isd,jsd) = area_c(isd+1,jsd+1)
1017  if (js == npy-1) area_c(isd,jed+1) = area_c(isd+1,jed)
1018  end if
1019  if (ie == npx-1) then
1020  do j=jsd,jed
1021  area_c(ied+1,j) = area_c(ied,j)
1022  end do
1023  if (js == 1) area_c(ied+1,jsd) = area_c(ied,jsd+1)
1024  if (js == npy-1) area_c(ied+1,jed+1) = area_c(ied,jed)
1025  end if
1026  if (js == 1) then
1027  do i=isd,ied
1028  area_c(i,jsd) = area_c(i,jsd+1)
1029  end do
1030  end if
1031  if (je == npy-1) then
1032  do i=isd,ied
1033  area_c(i,jed+1) = area_c(i,jed)
1034  end do
1035  end if
1036  end if
1037 
1038  call mpp_update_domains( area_c, atm%domain, position=corner, complete=.true.)
1039 
1040  ! Handle corner Area ghosting
1041  if (cubed_sphere .and. (.not. (atm%neststruct%nested .or. atm%flagstruct%regional))) then
1042  call fill_ghost(area, npx, npy, -big_number, atm%bd) ! fill in garbage values
1043  call fill_corners(area_c, npx, npy, fill=xdir, bgrid=.true.)
1044  endif
1045 
1046  do j=jsd,jed+1
1047  do i=isd,ied
1048  rdx(i,j) = 1.0/dx(i,j)
1049  enddo
1050  enddo
1051  do j=jsd,jed
1052  do i=isd,ied+1
1053  rdy(i,j) = 1.0/dy(i,j)
1054  enddo
1055  enddo
1056  do j=jsd,jed
1057  do i=isd,ied+1
1058  rdxc(i,j) = 1.0/dxc(i,j)
1059  enddo
1060  enddo
1061  do j=jsd,jed+1
1062  do i=isd,ied
1063  rdyc(i,j) = 1.0/dyc(i,j)
1064  enddo
1065  enddo
1066  do j=jsd,jed
1067  do i=isd,ied
1068  rarea(i,j) = 1.0/area(i,j)
1069  rdxa(i,j) = 1./dxa(i,j)
1070  rdya(i,j) = 1./dya(i,j)
1071  enddo
1072  enddo
1073  do j=jsd,jed+1
1074  do i=isd,ied+1
1075  rarea_c(i,j) = 1.0/area_c(i,j)
1076  enddo
1077  enddo
1078 
1079 200 format(a,f9.2,a,f9.2,a,f9.2)
1080 201 format(a,f9.2,a,f9.2,a,f9.2,a,f9.2)
1081 202 format(a,a,i4.4,a,i4.4,a)
1082 
1083  ! Get and print Grid Statistics
1084  dxav =0.0
1085  angav=0.0
1086  aspav=0.0
1087  dxn = missing
1088  dxm = -missing
1089  angn = missing
1090  angm = -missing
1091  aspn = missing
1092  aspm = -missing
1093  if (tile == 1) then
1094  do j=js, je
1095  do i=is, ie
1096  if(i>ceiling(npx/2.) .OR. j>ceiling(npy/2.)) cycle
1097  ang = get_angle(2, grid(i,j+1,1:2), grid(i,j,1:2), grid(i+1,j,1:2))
1098  ang = abs(90.0 - ang)
1099 
1100  if ( (i==1) .and. (j==1) ) then
1101  else
1102  angav = angav + ang
1103  angm = max(angm,ang)
1104  angn = min(angn,ang)
1105  endif
1106 
1107  dx_local = dx(i,j)
1108  dy_local = dy(i,j)
1109 
1110  dxav = dxav + 0.5 * (dx_local + dy_local)
1111  dxm = max(dxm,dx_local)
1112  dxm = max(dxm,dy_local)
1113  dxn = min(dxn,dx_local)
1114  dxn = min(dxn,dy_local)
1115 
1116  asp = abs(dx_local/dy_local)
1117  if (asp < 1.0) asp = 1.0/asp
1118  aspav = aspav + asp
1119  aspm = max(aspm,asp)
1120  aspn = min(aspn,asp)
1121  enddo
1122  enddo
1123  endif
1124  call mpp_sum(angav)
1125  call mpp_sum(dxav)
1126  call mpp_sum(aspav)
1127  call mpp_max(angm)
1128  call mpp_min(angn)
1129  call mpp_max(dxm)
1130  call mpp_min(dxn)
1131  call mpp_max(aspm)
1132  call mpp_min(aspn)
1133 
1134  if( is_master() ) then
1135 
1136  angav = angav / ( (ceiling(npy/2.0))*(ceiling(npx/2.0)) - 1 )
1137  dxav = dxav / ( (ceiling(npy/2.0))*(ceiling(npx/2.0)) )
1138  aspav = aspav / ( (ceiling(npy/2.0))*(ceiling(npx/2.0)) )
1139  write(*,* ) ''
1140 #ifdef SMALL_EARTH
1141  write(*,*) ' REDUCED EARTH: Radius is ', radius, ', omega is ', omega
1142 #endif
1143  write(*,* ) ' Cubed-Sphere Grid Stats : ', npx,'x',npy,'x',nregions
1144 !xxxxxxx write(*,201) ' Grid Length : min: ', dxN,' max: ', dxM,' avg: ', dxAV, ' min/max: ',dxN/dxM
1145  write(*,200) ' Deviation from Orthogonal : min: ',angn,' max: ',angm,' avg: ',angav
1146  write(*,200) ' Aspect Ratio : min: ',aspn,' max: ',aspm,' avg: ',aspav
1147  write(*,* ) ''
1148  endif
1149  endif!if gridtype > 3
1150 
1151  if (atm%neststruct%nested .or. any(atm%neststruct%child_grids)) then
1152  nullify(grid_global)
1153  else if( trim(grid_file) .NE. 'INPUT/grid_spec.nc') then
1154  deallocate(grid_global)
1155  endif
1156 
1157  nullify(agrid)
1158  nullify(grid)
1159 
1160  nullify( area)
1161  nullify(rarea)
1162  nullify( area_c)
1163  nullify(rarea_c)
1164 
1165  nullify(sina)
1166  nullify(cosa)
1167  nullify(dx)
1168  nullify(dy)
1169  nullify(dxc)
1170  nullify(dyc)
1171  nullify(dxa)
1172  nullify(dya)
1173  nullify(rdx)
1174  nullify(rdy)
1175  nullify(rdxc)
1176  nullify(rdyc)
1177  nullify(rdxa)
1178  nullify(rdya)
1179  nullify(e1)
1180  nullify(e2)
1181 
1182  nullify(iinta)
1183  nullify(jinta)
1184  nullify(iintb)
1185  nullify(jintb)
1186  nullify(npx_g)
1187  nullify(npy_g)
1188  nullify(ntiles_g)
1189  nullify(sw_corner)
1190  nullify(se_corner)
1191  nullify(ne_corner)
1192  nullify(nw_corner)
1193  nullify(latlon)
1194  nullify(cubed_sphere)
1195  nullify(have_south_pole)
1196  nullify(have_north_pole)
1197  nullify(stretched_grid)
1198 
1199  nullify(tile)
1200 
1201  nullify(domain)
1202 
1203  contains
1204 
1205  subroutine setup_cartesian(npx, npy, dx_const, dy_const, deglat, bd)
1206 
1207  type(fv_grid_bounds_type), intent(IN) :: bd
1208  integer, intent(in):: npx, npy
1209  real(kind=R_GRID), intent(IN) :: dx_const, dy_const, deglat
1210  real(kind=R_GRID) lat_rad, lon_rad, domain_rad
1211  integer i,j
1212  integer :: is, ie, js, je
1213  integer :: isd, ied, jsd, jed
1214 
1215  is = bd%is
1216  ie = bd%ie
1217  js = bd%js
1218  je = bd%je
1219  isd = bd%isd
1220  ied = bd%ied
1221  jsd = bd%jsd
1222  jed = bd%jed
1223 
1224  domain_rad = pi/16. ! arbitrary
1225  lat_rad = deglat * pi/180.
1226  lon_rad = 0. ! arbitrary
1227 
1228  dx(:,:) = dx_const
1229  rdx(:,:) = 1./dx_const
1230  dy(:,:) = dy_const
1231  rdy(:,:) = 1./dy_const
1232 
1233  dxc(:,:) = dx_const
1234  rdxc(:,:) = 1./dx_const
1235  dyc(:,:) = dy_const
1236  rdyc(:,:) = 1./dy_const
1237 
1238  dxa(:,:) = dx_const
1239  rdxa(:,:) = 1./dx_const
1240  dya(:,:) = dy_const
1241  rdya(:,:) = 1./dy_const
1242 
1243  area(:,:) = dx_const*dy_const
1244  rarea(:,:) = 1./(dx_const*dy_const)
1245 
1246  area_c(:,:) = dx_const*dy_const
1247  rarea_c(:,:) = 1./(dx_const*dy_const)
1248 
1249 ! The following is a hack to get pass the am2 phys init:
1250  do j=max(1,jsd),min(jed,npy)
1251  do i=max(1,isd),min(ied,npx)
1252  grid(i,j,1) = lon_rad - 0.5*domain_rad + real(i-1)/real(npx-1)*domain_rad
1253  grid(i,j,2) = lat_rad - 0.5*domain_rad + real(j-1)/real(npy-1)*domain_rad
1254  enddo
1255  enddo
1256 
1257  agrid(:,:,1) = lon_rad
1258  agrid(:,:,2) = lat_rad
1259 
1260  sina(:,:) = 1.
1261  cosa(:,:) = 0.
1262 
1263  e1(1,:,:) = 1.
1264  e1(2,:,:) = 0.
1265  e1(3,:,:) = 0.
1266 
1267  e2(1,:,:) = 0.
1268  e2(2,:,:) = 1.
1269  e2(3,:,:) = 0.
1270 
1271  end subroutine setup_cartesian
1272 
1273  subroutine setup_aligned_nest(Atm)
1274 
1275  type(fv_atmos_type), intent(INOUT), target :: Atm
1276 
1277  integer :: isd_p, ied_p, jsd_p, jed_p
1278  integer :: isg, ieg, jsg, jeg
1279  integer :: ic, jc, imod, jmod
1280 
1281 
1282  real(kind=R_GRID), allocatable, dimension(:,:,:) :: p_grid_u, p_grid_v, pa_grid, p_grid, c_grid_u, c_grid_v
1283  integer :: p_ind(1-ng:npx +ng,1-ng:npy +ng,4)
1286 
1287  integer i,j,k, p
1288  real(kind=R_GRID) sum
1289  real(kind=R_GRID) :: dist1, dist2, dist3, dist4
1290  real(kind=R_GRID), dimension(2) :: q1, q2
1291 
1292  integer, pointer :: parent_tile, refinement, ioffset, joffset
1293  integer, pointer, dimension(:,:,:) :: ind_h, ind_u, ind_v, ind_update_h
1294  real, pointer, dimension(:,:,:) :: wt_h, wt_u, wt_v
1295 
1296  integer, pointer, dimension(:,:,:) :: ind_b
1297  real, pointer, dimension(:,:,:) :: wt_b
1298 
1299  integer :: is, ie, js, je
1300  integer :: isd, ied, jsd, jed
1301 
1302  is = atm%bd%is
1303  ie = atm%bd%ie
1304  js = atm%bd%js
1305  je = atm%bd%je
1306  isd = atm%bd%isd
1307  ied = atm%bd%ied
1308  jsd = atm%bd%jsd
1309  jed = atm%bd%jed
1310 
1311 
1312 
1313  parent_tile => atm%neststruct%parent_tile
1314  refinement => atm%neststruct%refinement
1315  ioffset => atm%neststruct%ioffset
1316  joffset => atm%neststruct%joffset
1317 
1318  ind_h => atm%neststruct%ind_h
1319  ind_u => atm%neststruct%ind_u
1320  ind_v => atm%neststruct%ind_v
1321 
1322  ind_update_h => atm%neststruct%ind_update_h
1323 
1324  wt_h => atm%neststruct%wt_h
1325  wt_u => atm%neststruct%wt_u
1326  wt_v => atm%neststruct%wt_v
1327 
1328  ind_b => atm%neststruct%ind_b
1329  wt_b => atm%neststruct%wt_b
1330 
1331  call mpp_get_data_domain( atm%parent_grid%domain, &
1332  isd_p, ied_p, jsd_p, jed_p )
1333  call mpp_get_global_domain( atm%parent_grid%domain, &
1334  isg, ieg, jsg, jeg)
1335 
1336  allocate(p_grid_u(isg:ieg ,jsg:jeg+1,1:2))
1337  allocate(p_grid_v(isg:ieg+1,jsg:jeg ,1:2))
1338  allocate(pa_grid(isg:ieg,jsg:jeg ,1:2))
1339  p_ind = -1000000000
1340 
1341  allocate(p_grid( isg-ng:ieg+1+ng, jsg-ng:jeg+1+ng,1:2) )
1342  p_grid = 1.e25
1343 
1344  !Need to RECEIVE grid_global; matching mpp_send of grid_global from parent grid is in fv_control
1345  if( is_master() ) then
1346  p_ind = -1000000000
1347 
1348  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)), &
1349  atm%parent_grid%pelist(1))
1350  !Check that the grid does not lie outside its parent
1351  !3aug15: allows halo of nest to lie within halo of coarse grid.
1352  ! NOTE: will this then work with the mpp_update_nest_fine?
1353  if ( joffset + floor( real(1-ng) / real(refinement) ) < 1-ng .or. &
1354  ioffset + floor( real(1-ng) / real(refinement) ) < 1-ng .or. &
1355  joffset + floor( real(npy+ng) / real(refinement) ) > Atm%parent_grid%npy+ng .or. &
1356  ioffset + floor( real(npx+ng) / real(refinement) ) > Atm%parent_grid%npx+ng ) then
1357  call mpp_error(fatal, 'nested grid lies outside its parent')
1358  end if
1359 
1360  do j=1-ng,npy+ng
1361  jc = joffset + (j-1)/refinement !int( real(j-1) / real(refinement) )
1362  jmod = mod(j-1,refinement)
1363  if (j-1 < 0 .and. jmod /= 0) jc = jc - 1
1364  if (jmod < 0) jmod = jmod + refinement
1365 
1366  do i=1-ng,npx+ng
1367  ic = ioffset + (i-1)/refinement !int( real(i-1) / real(refinement) )
1368  imod = mod(i-1,refinement)
1369  if (i-1 < 0 .and. imod /= 0) ic = ic - 1
1370  if (imod < 0) imod = imod + refinement
1371 
1372  if (ic+1 > ieg+1 .or. ic < isg .or. jc+1 > jeg+1 .or. jc < jsg) then
1373  print*, 'p_grid:', i, j, ' OUT OF BOUNDS'
1374  print*, ic, jc
1375  print*, isg, ieg, jsg, jeg
1376  print*, imod, jmod
1377  end if
1378 
1379  if (jmod == 0) then
1380  q1 = p_grid(ic, jc, 1:2)
1381  q2 = p_grid(ic+1,jc,1:2)
1382  else
1383  call spherical_linear_interpolation( real(jmod,kind=r_grid)/real(refinement,kind=R_GRID), &
1384  p_grid(ic, jc, 1:2), p_grid(ic, jc+1, 1:2), q1)
1385  call spherical_linear_interpolation( real(jmod,kind=r_grid)/real(refinement,kind=R_GRID), &
1386  p_grid(ic+1, jc, 1:2), p_grid(ic+1, jc+1, 1:2), q2)
1387  end if
1388 
1389  if (imod == 0) then
1390  grid_global(i,j,:,1) = q1
1391  else
1392  call spherical_linear_interpolation( real(imod,kind=r_grid)/real(refinement,kind=R_GRID), &
1393  q1,q2,grid_global(i,j,:,1))
1394  end if
1395 
1396  !SW coarse-grid index; assumes grid does
1397  !not overlie other cube panels. (These indices
1398  !are also for the corners and thus need modification
1399  !to be used for cell-centered and edge-
1400  !centered variables; see below)
1401  p_ind(i,j,1) = ic
1402  p_ind(i,j,2) = jc
1403  p_ind(i,j,3) = imod
1404  p_ind(i,j,4) = jmod
1405 
1406  if (grid_global(i,j,1,1) > 2.*pi) grid_global(i,j,1,1) = grid_global(i,j,1,1) - 2.*pi
1407  if (grid_global(i,j,1,1) < 0.) grid_global(i,j,1,1) = grid_global(i,j,1,1) + 2.*pi
1408 
1409  end do
1410  end do
1411 
1412  ! Set up parent grids for interpolation purposes
1413  do j=jsg,jeg+1
1414  do i=isg,ieg
1415  call mid_pt_sphere(p_grid(i, j,1:2), p_grid(i+1, j,1:2), p_grid_u(i,j,:))
1416  !call mid_pt_sphere(p_grid(i, j,1:2), p_grid(i, j+1,1:2), p_grid_u(i,j,:))
1417  end do
1418  end do
1419  do j=jsg,jeg
1420  do i=isg,ieg+1
1421  call mid_pt_sphere(p_grid(i, j,1:2), p_grid(i, j+1,1:2), p_grid_v(i,j,:))
1422  !call mid_pt_sphere(p_grid(i, j,1:2), p_grid(i+1, j,1:2), p_grid_v(i,j,:))
1423  end do
1424  end do
1425  do j=jsg,jeg
1426  do i=isg,ieg
1427  call cell_center2(p_grid(i,j, 1:2), p_grid(i+1,j, 1:2), &
1428  p_grid(i,j+1,1:2), p_grid(i+1,j+1,1:2), &
1429  pa_grid(i,j,1:2) )
1430  end do
1431  end do
1432 
1433  end if
1434 
1435  call mpp_broadcast(grid_global(1-ng:npx+ng, 1-ng:npy+ng ,:,1), &
1436  ((npx+ng)-(1-ng)+1)*((npy+ng)-(1-ng)+1)*ndims, mpp_root_pe() )
1437  call mpp_broadcast( p_ind(1-ng:npx+ng, 1-ng:npy+ng ,1:4), &
1438  ((npx+ng)-(1-ng)+1)*((npy+ng)-(1-ng)+1)*4, mpp_root_pe() )
1439  call mpp_broadcast( pa_grid( isg:ieg , jsg:jeg , :), &
1440  ((ieg-isg+1))*(jeg-jsg+1)*ndims, mpp_root_pe())
1441  call mpp_broadcast( p_grid_u( isg:ieg , jsg:jeg+1, :), &
1442  (ieg-isg+1)*(jeg-jsg+2)*ndims, mpp_root_pe())
1443  call mpp_broadcast( p_grid_v( isg:ieg+1, jsg:jeg , :), &
1444  (ieg-isg+2)*(jeg-jsg+1)*ndims, mpp_root_pe())
1445 
1446  call mpp_broadcast( p_grid(isg-ng:ieg+ng+1, jsg-ng:jeg+ng+1, :), &
1447  (ieg-isg+2+2*ng)*(jeg-jsg+2+2*ng)*ndims, mpp_root_pe() )
1448 
1449  do n=1,ndims
1450  do j=jsd,jed+1
1451  do i=isd,ied+1
1452  grid(i,j,n) = grid_global(i,j,n,1)
1453  enddo
1454  enddo
1455  enddo
1456 
1457  ind_h = -999999999
1458  do j=jsd,jed
1459  do i=isd,ied
1460  ic = p_ind(i,j,1)
1461  jc = p_ind(i,j,2)
1462  imod = p_ind(i,j,3)
1463  jmod = p_ind(i,j,4)
1464 
1465 
1466  if (imod < refinement/2) then
1467 !!$ !!! DEBUG CODE
1468 !!$ if (ic /= ic) print*, gid, ' Bad ic ', i, j
1469 !!$ print*, i, j, ic
1470 !!$ !!! END DEBUG CODE
1471  ind_h(i,j,1) = ic - 1
1472  else
1473  ind_h(i,j,1) = ic
1474  end if
1475 
1476  if (jmod < refinement/2) then
1477  ind_h(i,j,2) = jc - 1
1478  else
1479  ind_h(i,j,2) = jc
1480  end if
1481  ind_h(i,j,3) = imod
1482  ind_h(i,j,4) = jmod
1483 
1484  end do
1485  end do
1486 
1487  ind_b = -999999999
1488  do j=jsd,jed+1
1489  do i=isd,ied+1
1490  ic = p_ind(i,j,1)
1491  jc = p_ind(i,j,2)
1492  imod = p_ind(i,j,3)
1493  jmod = p_ind(i,j,4)
1494 
1495  ind_b(i,j,1) = ic
1496  ind_b(i,j,2) = jc
1497 
1498  ind_b(i,j,3) = imod
1499  ind_b(i,j,4) = jmod
1500  enddo
1501  enddo
1502 
1503  !In a concurrent simulation, p_ind was passed off to the parent processes above, so they can create ind_update_h
1504 
1505  ind_u = -99999999
1506  !New BCs for wind components:
1507  ! For aligned grid segments (mod(j-1,R) == 0) set
1508  ! identically equal to the coarse-grid value
1509  ! Do linear interpolation in the y-dir elsewhere
1510 
1511  do j=jsd,jed+1
1512  do i=isd,ied
1513  ic = p_ind(i,j,1)
1514  jc = p_ind(i,j,2)
1515  imod = p_ind(i,j,3)
1516 
1517 #ifdef NEW_BC
1518  ind_u(i,j,1) = ic
1519 #else
1520  if (imod < refinement/2) then
1521 !!$ !!! DEBUG CODE
1522 !!$ print*, i, j, ic
1523 !!$ !!! END DEBUG CODE
1524  ind_u(i,j,1) = ic - 1
1525  else
1526  ind_u(i,j,1) = ic
1527  end if
1528 #endif
1529 
1530  ind_u(i,j,2) = jc
1531  ind_u(i,j,3) = imod
1532  ind_u(i,j,4) = p_ind(i,j,4)
1533 
1534  end do
1535  end do
1536 
1537  ind_v = -999999999
1538 
1539  do j=jsd,jed
1540  do i=isd,ied+1
1541  ic = p_ind(i,j,1)
1542  jc = p_ind(i,j,2)
1543  jmod = p_ind(i,j,4)
1544 
1545  ind_v(i,j,1) = ic
1546 
1547 #ifdef NEW_BC
1548  ind_v(i,j,2) = jc
1549 #else
1550  if (jmod < refinement/2) then
1551  ind_v(i,j,2) = jc - 1
1552  else
1553  ind_v(i,j,2) = jc
1554  end if
1555 #endif
1556 
1557  ind_v(i,j,4) = jmod
1558  ind_v(i,j,3) = p_ind(i,j,3)
1559  end do
1560  end do
1561 
1562 
1563 
1564  agrid(:,:,:) = -1.e25
1565 
1566  do j=jsd,jed
1567  do i=isd,ied
1568  call cell_center2(grid(i,j, 1:2), grid(i+1,j, 1:2), &
1569  grid(i,j+1,1:2), grid(i+1,j+1,1:2), &
1570  agrid(i,j,1:2) )
1571  enddo
1572  enddo
1573 
1574  call mpp_update_domains( agrid, atm%domain, position=center, complete=.true. )
1575 
1576  ! Compute dx
1577  do j=jsd,jed+1
1578  do i=isd,ied
1579  dx(i,j) = great_circle_dist(grid_global(i,j,:,1), grid_global(i+1,j,:,1), radius)
1580  enddo
1581  enddo
1582 
1583  ! Compute dy
1584  do j=jsd,jed
1585  do i=isd,ied+1
1586  dy(i,j) = great_circle_dist(grid_global(i,j,:,1), grid_global(i,j+1,:,1), radius)
1587  enddo
1588  enddo
1589 
1590  !We will use Michael Herzog's algorithm for computing the weights.
1591 
1592  !h points - need center (a-grid) points of parent grid
1593  do j=jsd,jed
1594  do i=isd,ied
1595 
1596  ic = ind_h(i,j,1)
1597  jc = ind_h(i,j,2)
1598 
1599  dist1 = dist2side_latlon(pa_grid(ic,jc,:) ,pa_grid(ic,jc+1,:),agrid(i,j,:))
1600  dist2 = dist2side_latlon(pa_grid(ic,jc+1,:) ,pa_grid(ic+1,jc+1,:),agrid(i,j,:))
1601  dist3 = dist2side_latlon(pa_grid(ic+1,jc+1,:),pa_grid(ic+1,jc,:),agrid(i,j,:))
1602  dist4 = dist2side_latlon(pa_grid(ic,jc,:) ,pa_grid(ic+1,jc,:),agrid(i,j,:))
1603 
1604  wt_h(i,j,1)=dist2*dist3 ! ic, jc weight
1605  wt_h(i,j,2)=dist3*dist4 ! ic, jc+1 weight
1606  wt_h(i,j,3)=dist4*dist1 ! ic+1, jc+1 weight
1607  wt_h(i,j,4)=dist1*dist2 ! ic+1, jc weight
1608 
1609  sum=wt_h(i,j,1)+wt_h(i,j,2)+wt_h(i,j,3)+wt_h(i,j,4)
1610  wt_h(i,j,:)=wt_h(i,j,:)/sum
1611 
1612  end do
1613  end do
1614 
1615 
1616 
1617  deallocate(pa_grid)
1618 
1619  do j=jsd,jed+1
1620  do i=isd,ied+1
1621 
1622  ic = ind_b(i,j,1)
1623  jc = ind_b(i,j,2)
1624 
1625  dist1 = dist2side_latlon(p_grid(ic,jc,:), p_grid(ic,jc+1,:), grid(i,j,:))
1626  dist2 = dist2side_latlon(p_grid(ic,jc+1,:), p_grid(ic+1,jc+1,:), grid(i,j,:))
1627  dist3 = dist2side_latlon(p_grid(ic+1,jc+1,:), p_grid(ic+1,jc,:), grid(i,j,:))
1628  dist4 = dist2side_latlon(p_grid(ic,jc,:), p_grid(ic+1,jc,:), grid(i,j,:))
1629 
1630  wt_b(i,j,1)=dist2*dist3 ! ic, jc weight
1631  wt_b(i,j,2)=dist3*dist4 ! ic, jc+1 weight
1632  wt_b(i,j,3)=dist4*dist1 ! ic+1, jc+1 weight
1633  wt_b(i,j,4)=dist1*dist2 ! ic+1, jc weight
1634 
1635  sum=wt_b(i,j,1)+wt_b(i,j,2)+wt_b(i,j,3)+wt_b(i,j,4)
1636  wt_b(i,j,:)=wt_b(i,j,:)/sum
1637 
1638  enddo
1639  enddo
1640 
1641  deallocate(p_grid)
1642 
1643 
1644  allocate(c_grid_u(isd:ied+1,jsd:jed,2))
1645  allocate(c_grid_v(isd:ied,jsd:jed+1,2))
1646 
1647  do j=jsd,jed
1648  do i=isd,ied+1
1649  call mid_pt_sphere(grid(i, j,1:2), grid(i, j+1,1:2), c_grid_u(i,j,:))
1650  end do
1651  end do
1652 
1653  do j=jsd,jed+1
1654  do i=isd,ied
1655  call mid_pt_sphere(grid(i,j ,1:2), grid(i+1,j ,1:2), c_grid_v(i,j,:))
1656  end do
1657  end do
1658 
1659 
1660  do j=jsd,jed
1661  do i=isd,ied
1662  dxa(i,j) = great_circle_dist(c_grid_u(i,j,:), c_grid_u(i+1,j,:), radius)
1663  end do
1664  end do
1665 
1666  do j=jsd,jed
1667  do i=isd,ied
1668  dya(i,j) = great_circle_dist(c_grid_v(i,j,:), c_grid_v(i,j+1,:), radius)
1669  end do
1670  end do
1671 
1672 
1673  !Compute interpolation weights. (Recall that the weights are defined with respect to a d-grid)
1674 
1675  !U weights
1676 
1677  do j=jsd,jed+1
1678  do i=isd,ied
1679 
1680  ic = ind_u(i,j,1)
1681  jc = ind_u(i,j,2)
1682 
1683  if (ic+1 > ieg .or. ic < isg .or. jc+1 > jeg+1 .or. jc < jsg) then
1684  print*, 'IND_U ', i, j, ' OUT OF BOUNDS'
1685  print*, ic, jc
1686  print*, isg, ieg, jsg, jeg
1687  end if
1688 
1689 
1690 #ifdef NEW_BC
1691  !New vorticity-conserving weights. Note that for the C-grid winds these
1692  ! become divergence-conserving weights!!
1693  jmod = p_ind(i,j,4)
1694  if (jmod == 0) then
1695  wt_u(i,j,1) = 1. ; wt_u(i,j,2) = 0.
1696  else
1697  dist1 = dist2side_latlon(p_grid_u(ic,jc,:), p_grid_u(ic+1,jc,:), c_grid_v(i,j,:))
1698  dist2 = dist2side_latlon(p_grid_u(ic,jc+1,:), p_grid_u(ic+1,jc+1,:), c_grid_v(i,j,:))
1699  sum = dist1+dist2
1700  wt_u(i,j,1) = dist2/sum
1701  wt_u(i,j,2) = dist1/sum
1702  endif
1703  wt_u(i,j,3:4) = 0.
1704 #else
1705  dist1 = dist2side_latlon(p_grid_u(ic,jc,:) ,p_grid_u(ic,jc+1,:), c_grid_v(i,j,:))
1706  dist2 = dist2side_latlon(p_grid_u(ic,jc+1,:) ,p_grid_u(ic+1,jc+1,:),c_grid_v(i,j,:))
1707  dist3 = dist2side_latlon(p_grid_u(ic+1,jc+1,:),p_grid_u(ic+1,jc,:), c_grid_v(i,j,:))
1708  !dist4 = dist2side_latlon(p_grid_u(ic+1,jc,:) ,p_grid_u(ic,jc,:), c_grid_v(i,j,:))
1709  dist4 = dist2side_latlon(p_grid_u(ic,jc,:) ,p_grid_u(ic+1,jc,:), c_grid_v(i,j,:))
1710 
1711  wt_u(i,j,1)=dist2*dist3 ! ic, jc weight
1712  wt_u(i,j,2)=dist3*dist4 ! ic, jc+1 weight
1713  wt_u(i,j,3)=dist4*dist1 ! ic+1, jc+1 weight
1714  wt_u(i,j,4)=dist1*dist2 ! ic+1, jc weight
1715 
1716  sum=wt_u(i,j,1)+wt_u(i,j,2)+wt_u(i,j,3)+wt_u(i,j,4)
1717  wt_u(i,j,:)=wt_u(i,j,:)/sum
1718 #endif
1719 
1720  end do
1721  end do
1722  !v weights
1723 
1724  do j=jsd,jed
1725  do i=isd,ied+1
1726 
1727  ic = ind_v(i,j,1)
1728  jc = ind_v(i,j,2)
1729 
1730  if (ic+1 > ieg .or. ic < isg .or. jc+1 > jeg+1 .or. jc < jsg) then
1731  print*, 'IND_V ', i, j, ' OUT OF BOUNDS'
1732  print*, ic, jc
1733  print*, isg, ieg, jsg, jeg
1734  end if
1735 
1736 #ifdef NEW_BC
1737  imod = p_ind(i,j,3)
1738  if (imod == 0) then
1739  wt_v(i,j,1) = 1. ; wt_v(i,j,4) = 0.
1740  else
1741  dist1 = dist2side_latlon(p_grid_v(ic,jc,:), p_grid_v(ic,jc+1,:), c_grid_u(i,j,:))
1742  dist2 = dist2side_latlon(p_grid_v(ic+1,jc,:), p_grid_v(ic+1,jc+1,:), c_grid_u(i,j,:))
1743  sum = dist1+dist2
1744  wt_v(i,j,1) = dist2/sum
1745  wt_v(i,j,4) = dist1/sum
1746  endif
1747  wt_v(i,j,2) = 0. ; wt_v(i,j,3) = 0.
1748 #else
1749  dist1 = dist2side_latlon(p_grid_v(ic,jc,:) ,p_grid_v(ic,jc+1,:), c_grid_u(i,j,:))
1750  dist2 = dist2side_latlon(p_grid_v(ic,jc+1,:) ,p_grid_v(ic+1,jc+1,:),c_grid_u(i,j,:))
1751  dist3 = dist2side_latlon(p_grid_v(ic+1,jc+1,:),p_grid_v(ic+1,jc,:), c_grid_u(i,j,:))
1752  dist4 = dist2side_latlon(p_grid_v(ic,jc,:) ,p_grid_v(ic+1,jc,:), c_grid_u(i,j,:))
1753 
1754  wt_v(i,j,1)=dist2*dist3 ! ic, jc weight
1755  wt_v(i,j,2)=dist3*dist4 ! ic, jc+1 weight
1756  wt_v(i,j,3)=dist4*dist1 ! ic+1, jc+1 weight
1757  wt_v(i,j,4)=dist1*dist2 ! ic+1, jc weight
1758 
1759  sum=wt_v(i,j,1)+wt_v(i,j,2)+wt_v(i,j,3)+wt_v(i,j,4)
1760  wt_v(i,j,:)=wt_v(i,j,:)/sum
1761 #endif
1762 
1763  end do
1764  end do
1765 
1766  deallocate(c_grid_u)
1767  deallocate(c_grid_v)
1768 
1769 
1770  deallocate(p_grid_u)
1771  deallocate(p_grid_v)
1772 
1773  if (is_master()) then
1774  if (atm%neststruct%nested) then
1775  !Nesting position information
1776  write(*,*) 'NESTED GRID ', atm%grid_number
1777  ic = p_ind(1,1,1) ; jc = p_ind(1,1,1)
1778  write(*,'(A, 2I5, 4F10.4)') 'SW CORNER: ', ic, jc, grid_global(1,1,:,1)*90./pi
1779  ic = p_ind(1,npy,1) ; jc = p_ind(1,npy,1)
1780  write(*,'(A, 2I5, 4F10.4)') 'NW CORNER: ', ic, jc, grid_global(1,npy,:,1)*90./pi
1781  ic = p_ind(npx,npy,1) ; jc = p_ind(npx,npy,1)
1782  write(*,'(A, 2I5, 4F10.4)') 'NE CORNER: ', ic, jc, grid_global(npx,npy,:,1)*90./pi
1783  ic = p_ind(npx,1,1) ; jc = p_ind(npx,1,1)
1784  write(*,'(A, 2I5, 4F10.4)') 'SE CORNER: ', ic, jc, grid_global(npx,1,:,1)*90./pi
1785  else
1786  write(*,*) 'PARENT GRID ', atm%parent_grid%grid_number, atm%parent_grid%tile
1787  ic = p_ind(1,1,1) ; jc = p_ind(1,1,1)
1788  write(*,'(A, 2I5, 4F10.4)') 'SW CORNER: ', ic, jc, atm%parent_grid%grid_global(ic,jc,:,parent_tile)*90./pi
1789  ic = p_ind(1,npy,1) ; jc = p_ind(1,npy,1)
1790  write(*,'(A, 2I5, 4F10.4)') 'NW CORNER: ', ic, jc, atm%parent_grid%grid_global(ic,jc,:,parent_tile)*90./pi
1791  ic = p_ind(npx,npy,1) ; jc = p_ind(npx,npy,1)
1792  write(*,'(A, 2I5, 4F10.4)') 'NE CORNER: ', ic, jc, atm%parent_grid%grid_global(ic,jc,:,parent_tile)*90./pi
1793  ic = p_ind(npx,1,1) ; jc = p_ind(npx,1,1)
1794  write(*,'(A, 2I5, 4F10.4)') 'SE CORNER: ', ic, jc, atm%parent_grid%grid_global(ic,jc,:,parent_tile)*90./pi
1795  endif
1796  end if
1797 
1798  end subroutine setup_aligned_nest
1799 
1800  subroutine setup_latlon(deglon_start,deglon_stop, deglat_start, deglat_stop, bd )
1801 
1802  type(fv_grid_bounds_type), intent(IN) :: bd
1803  real(kind=R_GRID), intent(IN) :: deglon_start,deglon_stop, deglat_start, deglat_stop
1804  real(kind=R_GRID) :: lon_start, lat_start, area_j
1805  integer :: is, ie, js, je
1806  integer :: isd, ied, jsd, jed
1807 
1808  is = bd%is
1809  ie = bd%ie
1810  js = bd%js
1811  je = bd%je
1812  isd = bd%isd
1813  ied = bd%ied
1814  jsd = bd%jsd
1815  jed = bd%jed
1816 
1817 
1818  dl = (deglon_stop-deglon_start)*pi/(180.*(npx-1))
1819  dp = (deglat_stop-deglat_start)*pi/(180.*(npy-1))
1820 
1821  lon_start = deglon_start*pi/180.
1822  lat_start = deglat_start*pi/180.
1823 
1824  do j=jsd,jed+1
1825  do i=isd,ied+1
1826  grid(i,j,1) = lon_start + real(i-1)*dl
1827  grid(i,j,2) = lat_start + real(j-1)*dp
1828  enddo
1829  enddo
1830 
1831  do j=jsd,jed
1832  do i=isd,ied
1833  agrid(i,j,1) = (grid(i,j,1) + grid(i+1,j,1))/2.
1834  agrid(i,j,2) = (grid(i,j,2) + grid(i,j+1,2))/2.
1835  enddo
1836  enddo
1837 
1838 
1839  do j=jsd,jed
1840  do i=isd,ied+1
1841  dxc(i,j) = dl*radius*cos(agrid(is,j,2))
1842  rdxc(i,j) = 1./dxc(i,j)
1843  enddo
1844  enddo
1845  do j=jsd,jed+1
1846  do i=isd,ied
1847  dyc(i,j) = dp*radius
1848  rdyc(i,j) = 1./dyc(i,j)
1849  enddo
1850  enddo
1851 
1852  do j=jsd,jed
1853  do i=isd,ied
1854  dxa(i,j) = dl*radius*cos(agrid(i,j,2))
1855  dya(i,j) = dp*radius
1856  rdxa(i,j) = 1./dxa(i,j)
1857  rdya(i,j) = 1./dya(i,j)
1858  enddo
1859  enddo
1860 
1861  do j=jsd,jed+1
1862  do i=isd,ied
1863  dx(i,j) = dl*radius*cos(grid(i,j,2))
1864  rdx(i,j) = 1./dx(i,j)
1865  enddo
1866  enddo
1867  do j=jsd,jed
1868  do i=isd,ied+1
1869  dy(i,j) = dp*radius
1870  rdy(i,j) = 1./dy(i,j)
1871  enddo
1872  enddo
1873 
1874  do j=jsd,jed
1875  area_j = radius*radius*dl*(sin(grid(is,j+1,2))-sin(grid(is,j,2)))
1876  do i=isd,ied
1877  area(i,j) = area_j
1878  rarea(i,j) = 1./area_j
1879  enddo
1880  enddo
1881 
1882  do j=jsd+1,jed
1883  area_j = radius*radius*dl*(sin(agrid(is,j,2))-sin(agrid(is,j-1,2)))
1884  do i=isd,ied+1
1885  area_c(i,j) = area_j
1886  rarea_c(i,j) = 1./area_j
1887  enddo
1888  enddo
1889  if (jsd==1) then
1890  j=1
1891  area_j = radius*radius*dl*(sin(agrid(is,j,2))-sin(agrid(is,j,2)-dp))
1892  do i=isd,ied+1
1893  area_c(i,j) = area_j
1894  rarea_c(i,j) = 1./area_j
1895  enddo
1896  endif
1897  if (jed+1==npy) then
1898  j=npy
1899  area_j = radius*radius*dl*(sin(agrid(is,j-1,2)+dp)-sin(agrid(is,j-1,2)))
1900  do i=isd,ied+1
1901  area_c(i,j) = area_j
1902  rarea_c(i,j) = 1./area_j
1903  enddo
1904  endif
1905  call mpp_update_domains( area_c, atm%domain, position=corner, complete=.true.)
1906 
1907  sina(:,:) = 1.
1908  cosa(:,:) = 0.
1909 
1910  e1(1,:,:) = 1.
1911  e1(2,:,:) = 0.
1912  e1(3,:,:) = 0.
1913 
1914  e2(1,:,:) = 0.
1915  e2(2,:,:) = 1.
1916  e2(3,:,:) = 0.
1917 
1918  end subroutine setup_latlon
1919 
1920  end subroutine init_grid
1921 
1922  subroutine cartesian_to_spherical(x, y, z, lon, lat, r)
1923  real(kind=R_GRID) , intent(IN) :: x, y, z
1924  real(kind=R_GRID) , intent(OUT) :: lon, lat, r
1925 
1926  r = sqrt(x*x + y*y + z*z)
1927  if ( (abs(x) + abs(y)) < 1.e-10 ) then ! poles:
1928  lon = 0.
1929  else
1930  lon = atan2(y,x) ! range: [-pi,pi]
1931  endif
1932 
1933 #ifdef RIGHT_HAND
1934  lat = asin(z/r)
1935 #else
1936  lat = acos(z/r) - pi/2.
1937 #endif
1938 
1939  end subroutine cartesian_to_spherical
1940  subroutine spherical_to_cartesian(lon, lat, r, x, y, z)
1941  real(kind=R_GRID) , intent(IN) :: lon, lat, r
1942  real(kind=R_GRID) , intent(OUT) :: x, y, z
1943 
1944  x = r * cos(lon) * cos(lat)
1945  y = r * sin(lon) * cos(lat)
1946 #ifdef RIGHT_HAND
1947  z = r * sin(lat)
1948 #else
1949  z = -r * sin(lat)
1950 #endif
1951  end subroutine spherical_to_cartesian
1952 
1955  subroutine rot_3d(axis, x1in, y1in, z1in, angle, x2out, y2out, z2out, degrees, convert)
1956  integer, intent(IN) :: axis
1957  real(kind=R_GRID) , intent(IN) :: x1in, y1in, z1in
1958  real(kind=R_GRID) , intent(INOUT) :: angle
1959  real(kind=R_GRID) , intent(OUT) :: x2out, y2out, z2out
1960  integer, intent(IN), optional :: degrees
1962  integer, intent(IN), optional :: convert
1965 
1966  real(kind=R_GRID) :: c, s
1967  real(kind=R_GRID) :: x1,y1,z1, x2,y2,z2
1968 
1969  if ( present(convert) ) then
1970  call spherical_to_cartesian(x1in, y1in, z1in, x1, y1, z1)
1971  else
1972  x1=x1in
1973  y1=y1in
1974  z1=z1in
1975  endif
1976 
1977  if ( present(degrees) ) then
1978  angle = angle*torad
1979  endif
1980 
1981  c = cos(angle)
1982  s = sin(angle)
1983 
1984  SELECT CASE(axis)
1985 
1986  CASE(1)
1987  x2 = x1
1988  y2 = c*y1 + s*z1
1989  z2 = -s*y1 + c*z1
1990  CASE(2)
1991  x2 = c*x1 - s*z1
1992  y2 = y1
1993  z2 = s*x1 + c*z1
1994  CASE(3)
1995  x2 = c*x1 + s*y1
1996  y2 = -s*x1 + c*y1
1997  z2 = z1
1998  CASE DEFAULT
1999  write(*,*) "Invalid axis: must be 1 for X, 2 for Y, 3 for Z."
2000 
2001  END SELECT
2002 
2003  if ( present(convert) ) then
2004  call cartesian_to_spherical(x2, y2, z2, x2out, y2out, z2out)
2005  else
2006  x2out=x2
2007  y2out=y2
2008  z2out=z2
2009  endif
2010 
2011  end subroutine rot_3d
2012 
2013 
2014 
2018  real(kind=R_GRID) function get_area_tri(ndims, p_1, p_2, p_3) &
2019  result(myarea)
2020 
2021  integer, intent(IN) :: ndims
2022  real(kind=R_GRID) , intent(IN) :: p_1(ndims) !
2023  real(kind=R_GRID) , intent(IN) :: p_2(ndims) !
2024  real(kind=R_GRID) , intent(IN) :: p_3(ndims) !
2025 
2026  real(kind=R_GRID) :: anga, angb, angc
2027 
2028  if ( ndims==3 ) then
2029  anga = spherical_angle(p_1, p_2, p_3)
2030  angb = spherical_angle(p_2, p_3, p_1)
2031  angc = spherical_angle(p_3, p_1, p_2)
2032  else
2033  anga = get_angle(ndims, p_1, p_2, p_3, 1)
2034  angb = get_angle(ndims, p_2, p_3, p_1, 1)
2035  angc = get_angle(ndims, p_3, p_1, p_2, 1)
2036  endif
2037 
2038  myarea = (anga+angb+angc - pi) * radius**2
2039 
2040  end function get_area_tri
2041 
2045  subroutine grid_area(nx, ny, ndims, nregions, nested, gridstruct, domain, bd, regional )
2046  type(fv_grid_bounds_type), intent(IN) :: bd
2047  integer, intent(IN) :: nx, ny, ndims, nregions
2048  logical, intent(IN) :: nested, regional
2049  type(fv_grid_type), intent(IN), target :: gridstruct
2050  type(domain2d), intent(INOUT) :: domain
2051 
2052  real(kind=R_GRID) :: p_lL(ndims)
2053  real(kind=R_GRID) :: p_uL(ndims)
2054  real(kind=R_GRID) :: p_lR(ndims)
2055  real(kind=R_GRID) :: p_uR(ndims)
2056  real(kind=R_GRID) :: a1, d1, d2, mydx, mydy, globalarea
2057 
2058  real(kind=R_GRID) :: p1(ndims), p2(ndims), p3(ndims), pi1(ndims), pi2(ndims)
2059 
2060  real(kind=R_GRID) :: maxarea, minarea
2061 
2062  integer :: i,j,n, nreg
2063  integer :: nh = 0
2064 
2065  real(kind=R_GRID), allocatable :: p_R8(:,:,:)
2066 
2067  real(kind=R_GRID), pointer, dimension(:,:,:) :: grid, agrid
2068  integer, pointer, dimension(:,:,:) :: iinta, jinta, iintb, jintb
2069  real(kind=R_GRID), pointer, dimension(:,:) :: area, area_c
2070 
2071  integer :: is, ie, js, je
2072  integer :: isd, ied, jsd, jed
2073 
2074  is = bd%is
2075  ie = bd%ie
2076  js = bd%js
2077  je = bd%je
2078  isd = bd%isd
2079  ied = bd%ied
2080  jsd = bd%jsd
2081  jed = bd%jed
2082 
2083  grid => gridstruct%grid_64
2084  agrid => gridstruct%agrid_64
2085  iinta => gridstruct%iinta
2086  jinta => gridstruct%jinta
2087  iintb => gridstruct%iintb
2088  jintb => gridstruct%jintb
2089 
2090  area => gridstruct%area_64
2091  area_c => gridstruct%area_c_64
2092 
2093  if (nested .or. regional) nh = ng
2094 
2095  maxarea = -1.e25
2096  minarea = 1.e25
2097 
2098  globalarea = 0.0
2099  do j=js-nh,je+nh
2100  do i=is-nh,ie+nh
2101  do n=1,ndims
2102  if ( gridstruct%stretched_grid .or. nested ) then
2103  p_ll(n) = grid(i ,j ,n)
2104  p_ul(n) = grid(i ,j+1,n)
2105  p_lr(n) = grid(i+1,j ,n)
2106  p_ur(n) = grid(i+1,j+1,n)
2107  else
2108  p_ll(n) = grid(iinta(1,i,j), jinta(1,i,j), n)
2109  p_ul(n) = grid(iinta(2,i,j), jinta(2,i,j), n)
2110  p_lr(n) = grid(iinta(4,i,j), jinta(4,i,j), n)
2111  p_ur(n) = grid(iinta(3,i,j), jinta(3,i,j), n)
2112  endif
2113  enddo
2114 
2115  ! Spherical Excess Formula
2116  area(i,j) = get_area(p_ll, p_ul, p_lr, p_ur, radius)
2117  !maxarea=MAX(area(i,j),maxarea)
2118  !minarea=MIN(area(i,j),minarea)
2119  !globalarea = globalarea + area(i,j)
2120  enddo
2121  enddo
2122 
2123 !!$ allocate( p_R8(nx-1,ny-1,ntiles_g) ) ! this is a "global" array
2124 !!$ do j=js,je
2125 !!$ do i=is,ie
2126 !!$ p_R8(i,j,tile) = area(i,j)
2127 !!$ enddo
2128 !!$ enddo
2129 !!$ call mp_gather(p_R8, is,ie, js,je, nx-1, ny-1, ntiles_g)
2130 !!$ if (is_master()) then
2131 !!$ globalarea = 0.0
2132 !!$ do n=1,ntiles_g
2133 !!$ do j=1,ny-1
2134 !!$ do i=1,nx-1
2135 !!$ globalarea = globalarea + p_R8(i,j,n)
2136 !!$ enddo
2137 !!$ enddo
2138 !!$ enddo
2139 !!$ endif
2140 !!$
2141 !!$ call mpp_broadcast(globalarea, mpp_root_pe())
2142 !!$
2143 !!$ deallocate( p_R8 )
2144 !!$
2145 !!$ call mp_reduce_max(maxarea)
2146 !!$ minarea = -minarea
2147 !!$ call mp_reduce_max(minarea)
2148 !!$ minarea = -minarea
2149 
2150  globalarea = mpp_global_sum(domain, area)
2151  maxarea = mpp_global_max(domain, area)
2152  minarea = mpp_global_min(domain, area)
2153 
2154  if (is_master()) write(*,209) 'MAX AREA (m*m):', maxarea, ' MIN AREA (m*m):', minarea
2155  if (is_master()) write(*,209) 'GLOBAL AREA (m*m):', globalarea, ' IDEAL GLOBAL AREA (m*m):', 4.0*pi*radius**2
2156  209 format(a,e21.14,a,e21.14)
2157 
2158  if (nested .or. regional) then
2159  nh = ng-1 !cannot get rarea_c on boundary directly
2160  area_c = 1.e30
2161  end if
2162 
2163  do j=js-nh,je+nh+1
2164  do i=is-nh,ie+nh+1
2165  do n=1,ndims
2166  if ( gridstruct%stretched_grid .or. nested ) then
2167  p_ll(n) = agrid(i-1,j-1,n)
2168  p_lr(n) = agrid(i ,j-1,n)
2169  p_ul(n) = agrid(i-1,j ,n)
2170  p_ur(n) = agrid(i ,j ,n)
2171  else
2172  p_ll(n) = agrid(iintb(1,i,j), jintb(1,i,j), n)
2173  p_lr(n) = agrid(iintb(2,i,j), jintb(2,i,j), n)
2174  p_ul(n) = agrid(iintb(4,i,j), jintb(4,i,j), n)
2175  p_ur(n) = agrid(iintb(3,i,j), jintb(3,i,j), n)
2176  endif
2177  enddo
2178  ! Spherical Excess Formula
2179  area_c(i,j) = get_area(p_ll, p_ul, p_lr, p_ur, radius)
2180  enddo
2181  enddo
2182 
2183 ! Corners: assuming triangular cells
2184  if (gridstruct%cubed_sphere .and. .not. (nested .or. regional)) then
2185 ! SW:
2186  i=1
2187  j=1
2188  if ( (is==1) .and. (js==1) ) then
2189  do n=1,ndims
2190  if ( gridstruct%stretched_grid ) then
2191  p1(n) = agrid(i-1,j ,n)
2192  p2(n) = agrid(i ,j ,n)
2193  p3(n) = agrid(i ,j-1,n)
2194  else
2195  p1(n) = agrid(iintb(1,i,j), jintb(1,i,j), n)
2196  p2(n) = agrid(iintb(2,i,j), jintb(2,i,j), n)
2197  p3(n) = agrid(iintb(3,i,j), jintb(3,i,j), n)
2198  endif
2199  enddo
2200  area_c(i,j) = get_area_tri(ndims, p1, p2, p3)
2201  endif
2202 
2203  i=nx
2204  j=1
2205  if ( (ie+1==nx) .and. (js==1) ) then
2206  do n=1,ndims
2207  if ( gridstruct%stretched_grid ) then
2208  p1(n) = agrid(i ,j ,n)
2209  p2(n) = agrid(i-1,j ,n)
2210  p3(n) = agrid(i-1,j-1,n)
2211  else
2212  p1(n) = agrid(iintb(1,i,j), jintb(1,i,j), n)
2213  p2(n) = agrid(iintb(2,i,j), jintb(2,i,j), n)
2214  p3(n) = agrid(iintb(3,i,j), jintb(3,i,j), n)
2215  endif
2216  enddo
2217  area_c(i,j) = get_area_tri(ndims, p1, p2, p3)
2218  endif
2219 
2220  i=nx
2221  j=ny
2222  if ( (ie+1==nx) .and. (je+1==ny) ) then
2223  do n=1,ndims
2224  if ( gridstruct%stretched_grid ) then
2225  p1(n) = agrid(i-1,j ,n)
2226  p2(n) = agrid(i-1,j-1,n)
2227  p3(n) = agrid(i ,j-1,n)
2228  else
2229  p1(n) = agrid(iintb(1,i,j), jintb(1,i,j), n)
2230  p2(n) = agrid(iintb(2,i,j), jintb(2,i,j), n)
2231  p3(n) = agrid(iintb(3,i,j), jintb(3,i,j), n)
2232  endif
2233  enddo
2234  area_c(i,j) = get_area_tri(ndims, p1, p2, p3)
2235  endif
2236 
2237  i=1
2238  j=ny
2239  if ( (is==1) .and. (je+1==ny) ) then
2240  do n=1,ndims
2241  if ( gridstruct%stretched_grid ) then
2242  p1(n) = agrid(i ,j ,n)
2243  p2(n) = agrid(i ,j-1,n)
2244  p3(n) = agrid(i-1,j-1,n)
2245  else
2246  p1(n) = agrid(iintb(1,i,j), jintb(1,i,j), n)
2247  p2(n) = agrid(iintb(2,i,j), jintb(2,i,j), n)
2248  p3(n) = agrid(iintb(3,i,j), jintb(3,i,j), n)
2249  endif
2250  enddo
2251  area_c(i,j) = get_area_tri(ndims, p1, p2, p3)
2252  endif
2253  endif
2254 
2255  end subroutine grid_area
2256 
2261  real(kind=R_GRID) function get_angle(ndims, p1, p2, p3, rad) result (angle)
2262 
2263  integer, intent(IN) :: ndims
2264  real(kind=R_GRID) , intent(IN) :: p1(ndims)
2265  real(kind=R_GRID) , intent(IN) :: p2(ndims)
2266  real(kind=R_GRID) , intent(IN) :: p3(ndims)
2267  integer, intent(in), optional:: rad
2268 
2269  real(kind=R_GRID) :: e1(3), e2(3), e3(3)
2270 
2271  if (ndims == 2) then
2272  call spherical_to_cartesian(p2(1), p2(2), real(1.,kind=R_GRID), e1(1), e1(2), e1(3))
2273  call spherical_to_cartesian(p1(1), p1(2), real(1.,kind=R_GRID), e2(1), e2(2), e2(3))
2274  call spherical_to_cartesian(p3(1), p3(2), real(1.,kind=R_GRID), e3(1), e3(2), e3(3))
2275  else
2276  e1 = p2; e2 = p1; e3 = p3
2277  endif
2278 
2279 ! High precision version:
2280  if ( present(rad) ) then
2281  angle = spherical_angle(e1, e2, e3)
2282  else
2283  angle = todeg * spherical_angle(e1, e2, e3)
2284  endif
2285 
2286  end function get_angle
2287 
2289  subroutine mirror_grid(grid_global,ng,npx,npy,ndims,nregions)
2290  integer, intent(IN) :: ng,npx,npy,ndims,nregions
2291  real(kind=R_GRID) , intent(INOUT) :: grid_global(1-ng:npx +ng,1-ng:npy +ng,ndims,1:nregions)
2292  integer :: i,j,n,n1,n2,nreg
2293  real(kind=R_GRID) :: x1,y1,z1, x2,y2,z2, ang
2294 
2295  nreg = 1
2296  do j=1,ceiling(npy/2.)
2297  do i=1,ceiling(npx/2.)
2298 
2299  x1 = 0.25d0 * (abs(grid_global(i ,j ,1,nreg)) + &
2300  abs(grid_global(npx-(i-1),j ,1,nreg)) + &
2301  abs(grid_global(i ,npy-(j-1),1,nreg)) + &
2302  abs(grid_global(npx-(i-1),npy-(j-1),1,nreg)))
2303  grid_global(i ,j ,1,nreg) = sign(x1,grid_global(i ,j ,1,nreg))
2304  grid_global(npx-(i-1),j ,1,nreg) = sign(x1,grid_global(npx-(i-1),j ,1,nreg))
2305  grid_global(i ,npy-(j-1),1,nreg) = sign(x1,grid_global(i ,npy-(j-1),1,nreg))
2306  grid_global(npx-(i-1),npy-(j-1),1,nreg) = sign(x1,grid_global(npx-(i-1),npy-(j-1),1,nreg))
2307 
2308  y1 = 0.25d0 * (abs(grid_global(i ,j ,2,nreg)) + &
2309  abs(grid_global(npx-(i-1),j ,2,nreg)) + &
2310  abs(grid_global(i ,npy-(j-1),2,nreg)) + &
2311  abs(grid_global(npx-(i-1),npy-(j-1),2,nreg)))
2312  grid_global(i ,j ,2,nreg) = sign(y1,grid_global(i ,j ,2,nreg))
2313  grid_global(npx-(i-1),j ,2,nreg) = sign(y1,grid_global(npx-(i-1),j ,2,nreg))
2314  grid_global(i ,npy-(j-1),2,nreg) = sign(y1,grid_global(i ,npy-(j-1),2,nreg))
2315  grid_global(npx-(i-1),npy-(j-1),2,nreg) = sign(y1,grid_global(npx-(i-1),npy-(j-1),2,nreg))
2316 
2317  ! force dateline/greenwich-meridion consitency
2318  if (mod(npx,2) /= 0) then
2319  if ( (i==1+(npx-1)/2.0d0) ) then
2320  grid_global(i,j ,1,nreg) = 0.0d0
2321  grid_global(i,npy-(j-1),1,nreg) = 0.0d0
2322  endif
2323  endif
2324 
2325  enddo
2326  enddo
2327 
2328  do nreg=2,nregions
2329  do j=1,npy
2330  do i=1,npx
2331 
2332  x1 = grid_global(i,j,1,1)
2333  y1 = grid_global(i,j,2,1)
2334  z1 = radius
2335 
2336  if (nreg == 2) then
2337  ang = -90.d0
2338  call rot_3d( 3, x1, y1, z1, ang, x2, y2, z2, 1, 1) ! rotate about the z-axis
2339  elseif (nreg == 3) then
2340  ang = -90.d0
2341  call rot_3d( 3, x1, y1, z1, ang, x2, y2, z2, 1, 1) ! rotate about the z-axis
2342  ang = 90.d0
2343  call rot_3d( 1, x2, y2, z2, ang, x1, y1, z1, 1, 1) ! rotate about the x-axis
2344  x2=x1
2345  y2=y1
2346  z2=z1
2347 
2348  ! force North Pole and dateline/greenwich-meridion consitency
2349  if (mod(npx,2) /= 0) then
2350  if ( (i==1+(npx-1)/2.0d0) .and. (i==j) ) then
2351  x2 = 0.0d0
2352  y2 = pi/2.0d0
2353  endif
2354  if ( (j==1+(npy-1)/2.0d0) .and. (i < 1+(npx-1)/2.0d0) ) then
2355  x2 = 0.0d0
2356  endif
2357  if ( (j==1+(npy-1)/2.0d0) .and. (i > 1+(npx-1)/2.0d0) ) then
2358  x2 = pi
2359  endif
2360  endif
2361 
2362  elseif (nreg == 4) then
2363  ang = -180.d0
2364  call rot_3d( 3, x1, y1, z1, ang, x2, y2, z2, 1, 1) ! rotate about the z-axis
2365  ang = 90.d0
2366  call rot_3d( 1, x2, y2, z2, ang, x1, y1, z1, 1, 1) ! rotate about the x-axis
2367  x2=x1
2368  y2=y1
2369  z2=z1
2370 
2371  ! force dateline/greenwich-meridion consitency
2372  if (mod(npx,2) /= 0) then
2373  if ( (j==1+(npy-1)/2.0d0) ) then
2374  x2 = pi
2375  endif
2376  endif
2377 
2378  elseif (nreg == 5) then
2379  ang = 90.d0
2380  call rot_3d( 3, x1, y1, z1, ang, x2, y2, z2, 1, 1) ! rotate about the z-axis
2381  ang = 90.d0
2382  call rot_3d( 2, x2, y2, z2, ang, x1, y1, z1, 1, 1) ! rotate about the y-axis
2383  x2=x1
2384  y2=y1
2385  z2=z1
2386  elseif (nreg == 6) then
2387  ang = 90.d0
2388  call rot_3d( 2, x1, y1, z1, ang, x2, y2, z2, 1, 1) ! rotate about the y-axis
2389  ang = 0.d0
2390  call rot_3d( 3, x2, y2, z2, ang, x1, y1, z1, 1, 1) ! rotate about the z-axis
2391  x2=x1
2392  y2=y1
2393  z2=z1
2394 
2395  ! force South Pole and dateline/greenwich-meridion consitency
2396  if (mod(npx,2) /= 0) then
2397  if ( (i==1+(npx-1)/2.0d0) .and. (i==j) ) then
2398  x2 = 0.0d0
2399  y2 = -pi/2.0d0
2400  endif
2401  if ( (i==1+(npx-1)/2.0d0) .and. (j > 1+(npy-1)/2.0d0) ) then
2402  x2 = 0.0d0
2403  endif
2404  if ( (i==1+(npx-1)/2.0d0) .and. (j < 1+(npy-1)/2.0d0) ) then
2405  x2 = pi
2406  endif
2407  endif
2408 
2409  endif
2410 
2411  grid_global(i,j,1,nreg) = x2
2412  grid_global(i,j,2,nreg) = y2
2413 
2414  enddo
2415  enddo
2416  enddo
2417 
2418  end subroutine mirror_grid
2419 
2420  end module fv_grid_tools_mod
2421 
fv_grid_utils_mod::mid_pt_sphere
subroutine, public mid_pt_sphere(p1, p2, pm)
Definition: fv_grid_utils.F90:1954
fv_grid_tools_mod::torad
real(kind=r_grid), parameter torad
convert to radians
Definition: fv_grid_tools.F90:160
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:120
fv_arrays_mod::fv_grid_bounds_type
Definition: fv_arrays.F90:1137
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:161
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:1801
fv_grid_tools_mod::radius
real(kind=r_grid), parameter radius
Definition: fv_grid_tools.F90:157
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:918
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:2020
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:2770
fv_grid_utils_mod::spherical_angle
real(kind=r_grid) function, public spherical_angle(p1, p2, p3)
Definition: fv_grid_utils.F90:2795
fv_grid_utils_mod::inner_prod
real function, public inner_prod(v1, v2)
Definition: fv_grid_utils.F90:973
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:1902
setup_aligned_nest
subroutine setup_aligned_nest(Atm)
Definition: fv_grid_tools.F90:1274
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:1956
fv_grid_tools_mod::todeg
real(kind=r_grid), parameter, public todeg
convert to degrees
Definition: fv_grid_tools.F90:159
fv_arrays_mod::r_grid
integer, parameter, public r_grid
Definition: fv_arrays.F90:35
fv_timing_mod
The module &#39;fv_timing&#39; contains FV3 timers.
Definition: fv_timing.F90:24
fv_grid_tools_mod::init_grid
subroutine, public init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, ng)
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_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:2013
fv_grid_tools_mod::grid_area
subroutine grid_area(nx, ny, ndims, nregions, nested, gridstruct, domain, bd, regional)
The subroutine &#39;grid_area&#39; gets the surface area on a grid in lat/lon or xyz coordinates.
Definition: fv_grid_tools.F90:2046
fv_grid_utils_mod::get_area
real(kind=r_grid) function, public get_area(p1, p4, p2, p3, radius)
Definition: fv_grid_utils.F90:2708
fv_grid_utils_mod::cell_center2
subroutine, public cell_center2(q1, q2, q3, q4, e2)
Definition: fv_grid_utils.F90:2660
fv_grid_tools_mod::csfac
real(kind=r_grid) csfac
Definition: fv_grid_tools.F90:163
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:175
fv_grid_tools_mod::cartesian_to_spherical
subroutine cartesian_to_spherical(x, y, z, lon, lat, r)
Definition: fv_grid_tools.F90:1923
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_mp_mod::ng
integer, parameter, public ng
Definition: fv_mp_mod.F90:2716
fv_grid_utils_mod::fill_ghost
Definition: fv_grid_utils.F90:124
fv_grid_tools_mod::spherical_to_cartesian
subroutine, public spherical_to_cartesian(lon, lat, r, x, y, z)
Definition: fv_grid_tools.F90:1941
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:1166
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_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:166
fv_grid_tools_mod::debug_message_size
logical, parameter debug_message_size
Definition: fv_grid_tools.F90:165
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:2290
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:2262
fv_grid_utils_mod::gnomonic_grids
subroutine, public gnomonic_grids(grid_type, im, lon, lat)
Definition: fv_grid_utils.F90:1274
setup_cartesian
subroutine setup_cartesian(npx, npy, dx_const, dy_const, deglat, bd)
Definition: fv_grid_tools.F90:1206