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