FV3DYCORE  Version 1.1.0
fv_grid_utils.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 
25 
26  module fv_grid_utils_mod
27 
28 ! Modules Included:
29 ! <table>
30 ! <tr>
31 ! <th>Module Name</th>
32 ! <th>Functions Included</th>
33 ! </tr>
34 ! <tr>
35 ! <td>constants_mod</td>
36 ! <td>omega, pi=>pi_8, cnst_radius=>radius</td>
37 ! </tr>
38 ! <tr>
39 ! <td>external_sst_mod</td>
40 ! <td>i_sst, j_sst, sst_ncep, sst_anom</td>
41 ! </tr>
42 ! <tr>
43 ! <td>fv_arrays_mod</td>
44 ! <td>fv_atmos_type, fv_grid_type, fv_grid_bounds_type,R_GRID</td>
45 ! </tr>
46 ! <tr>
47 ! <td>fv_diagnostics_mod</td>
48 ! <td>fv_time, prt_mxm, range_check, prt_minmax</td>
49 ! </tr>
50 ! <tr>
51 ! <td>fv_eta_mod</td>
52 ! <td>set_eta</td>
53 ! </tr>
54 ! <tr>
55 ! <td>fv_mp_mod</td>
56 ! <td>ng, is_master, mp_reduce_sum, mp_reduce_min, mp_reduce_max,fill_corners, XDir, YDir</td>
57 ! </tr>
58 ! <tr>
59 ! <td>fv_parameter_mod</td>
60 ! <td>AGRID_PARAM=>AGRID, CGRID_NE_PARAM=>CGRID_NE, CORNER, SCALAR_PAIR</td>
61 ! </tr>
62 ! <tr>
63 ! <td>fv_timing_mod</td>
64 ! <td>timing_on, timing_off</td>
65 ! </tr>
66 ! <tr>
67 ! <td>mpp_mod/td>
68 ! <td>FATAL, mpp_error, WARNING</td>
69 ! </tr>
70 ! <tr>
71 ! <td>mpp_domains_mod/td>
72 ! <td>mpp_update_domains, DGRID_NE, mpp_global_sum, BITWISE_EXACT_SUM, domain2d, BITWISE_EFP_SUM</td>
73 ! </tr>
74 ! </table>
75 
76 
77 #include <fms_platform.h>
78  use constants_mod, only: omega, pi=>pi_8, cnst_radius=>radius
79  use mpp_mod, only: fatal, mpp_error, warning
80  use external_sst_mod, only: i_sst, j_sst, sst_ncep, sst_anom
81  use mpp_domains_mod, only: mpp_update_domains, dgrid_ne, mpp_global_sum
82  use mpp_domains_mod, only: bitwise_exact_sum, domain2d, bitwise_efp_sum
83  use mpp_parameter_mod, only: agrid_param=>agrid, cgrid_ne_param=>cgrid_ne
84  use mpp_parameter_mod, only: corner, scalar_pair
85 
86  use fv_arrays_mod, only: fv_atmos_type, fv_grid_type, fv_grid_bounds_type, &
87  r_grid
88  use fv_eta_mod, only: set_eta
89  use fv_mp_mod, only: ng, is_master
90  use fv_mp_mod, only: mp_reduce_sum, mp_reduce_min, mp_reduce_max
91  use fv_mp_mod, only: fill_corners, xdir, ydir
92  use fv_timing_mod, only: timing_on, timing_off
93 
94  implicit none
95  private
96  logical:: symm_grid
97 #ifdef NO_QUAD_PRECISION
98 ! 64-bit precision (kind=8)
99  integer, parameter:: f_p = selected_real_kind(15)
100 #else
101 ! Higher precision (kind=16) for grid geometrical factors:
102  integer, parameter:: f_p = selected_real_kind(20)
103 #endif
104  real, parameter:: big_number=1.d8
105  real, parameter:: tiny_number=1.d-8
106 
107  real(kind=R_GRID) :: radius=cnst_radius
108 
109  real, parameter:: ptop_min=1.d-8
110 
111  public f_p
112  public ptop_min, big_number !CLEANUP: OK to keep since they are constants?
113  public cos_angle
114  public latlon2xyz, gnomonic_grids, &
115  global_mx, unit_vect_latlon, &
116  cubed_to_latlon, c2l_ord2, g_sum, global_qsum, great_circle_dist, &
117  v_prod, get_unit_vect2, project_sphere_v
118  public mid_pt_sphere, mid_pt_cart, vect_cross, grid_utils_init, grid_utils_end, &
119  spherical_angle, cell_center2, get_area, inner_prod, fill_ghost, direct_transform, &
120  make_eta_level, expand_cell, cart_to_latlon, intp_great_circle, normalize_vect, &
121  dist2side_latlon, spherical_linear_interpolation, get_latlon_vector
122  public symm_grid
123 
124  INTERFACE fill_ghost
125 #ifdef OVERLOAD_R4
126  MODULE PROCEDURE fill_ghost_r4
127 #endif
128  MODULE PROCEDURE fill_ghost_r8
129  END INTERFACE
130 
131  contains
132 
133  subroutine grid_utils_init(Atm, npx, npy, npz, non_ortho, grid_type, c2l_order)
134 
135  type(fv_atmos_type), intent(inout), target :: Atm
136  logical, intent(in):: non_ortho
137  integer, intent(in):: npx, npy, npz
138  integer, intent(in):: grid_type, c2l_order
139 !
140 ! Super (composite) grid:
141 
142 ! 9---4---8
143 ! | |
144 ! 1 5 3
145 ! | |
146 ! 6---2---7
147 
148  real(kind=R_GRID) grid3(3,atm%bd%isd:atm%bd%ied+1,atm%bd%jsd:atm%bd%jed+1)
149  real(kind=R_GRID) p1(3), p2(3), p3(3), p4(3), pp(3), ex(3), ey(3), e1(3), e2(3)
150  real(kind=R_GRID) pp1(2), pp2(2), pp3(2)
151  real(kind=R_GRID) sin2, tmp1, tmp2
152  integer i, j, k, n, ip
153 
154  integer :: is, ie, js, je
155  integer :: isd, ied, jsd, jed
156 
157  !Local pointers
158  real(kind=R_GRID), pointer, dimension(:,:,:) :: agrid, grid
159  real(kind=R_GRID), pointer, dimension(:,:) :: area, area_c
160  real(kind=R_GRID), pointer, dimension(:,:) :: sina, cosa, dx, dy, dxc, dyc, dxa, dya
161  real, pointer, dimension(:,:) :: del6_u, del6_v
162  real, pointer, dimension(:,:) :: divg_u, divg_v
163  real, pointer, dimension(:,:) :: cosa_u, cosa_v, cosa_s
164  real, pointer, dimension(:,:) :: sina_u, sina_v
165  real, pointer, dimension(:,:) :: rsin_u, rsin_v
166  real, pointer, dimension(:,:) :: rsina, rsin2
167  real, pointer, dimension(:,:,:) :: sin_sg, cos_sg
168  real(kind=R_GRID), pointer, dimension(:,:,:) :: ee1, ee2, ec1, ec2
169  real(kind=R_GRID), pointer, dimension(:,:,:,:) :: ew, es
170  real(kind=R_GRID), pointer, dimension(:,:,:) :: en1, en2
171 ! real(kind=R_GRID), pointer, dimension(:,:) :: eww, ess
172  logical, pointer :: sw_corner, se_corner, ne_corner, nw_corner
173 
174  is = atm%bd%is
175  ie = atm%bd%ie
176  js = atm%bd%js
177  je = atm%bd%je
178  isd = atm%bd%isd
179  ied = atm%bd%ied
180  jsd = atm%bd%jsd
181  jed = atm%bd%jed
182 
183 !--- pointers to higher-order precision quantities
184  agrid => atm%gridstruct%agrid_64
185  grid => atm%gridstruct%grid_64
186  area => atm%gridstruct%area_64
187  area_c => atm%gridstruct%area_c_64
188  dx => atm%gridstruct%dx_64
189  dy => atm%gridstruct%dy_64
190  dxc => atm%gridstruct%dxc_64
191  dyc => atm%gridstruct%dyc_64
192  dxa => atm%gridstruct%dxa_64
193  dya => atm%gridstruct%dya_64
194  sina => atm%gridstruct%sina_64
195  cosa => atm%gridstruct%cosa_64
196 
197  divg_u => atm%gridstruct%divg_u
198  divg_v => atm%gridstruct%divg_v
199 
200  del6_u => atm%gridstruct%del6_u
201  del6_v => atm%gridstruct%del6_v
202 
203  cosa_u => atm%gridstruct%cosa_u
204  cosa_v => atm%gridstruct%cosa_v
205  cosa_s => atm%gridstruct%cosa_s
206  sina_u => atm%gridstruct%sina_u
207  sina_v => atm%gridstruct%sina_v
208  rsin_u => atm%gridstruct%rsin_u
209  rsin_v => atm%gridstruct%rsin_v
210  rsina => atm%gridstruct%rsina
211  rsin2 => atm%gridstruct%rsin2
212  ee1 => atm%gridstruct%ee1
213  ee2 => atm%gridstruct%ee2
214  ec1 => atm%gridstruct%ec1
215  ec2 => atm%gridstruct%ec2
216  ew => atm%gridstruct%ew
217  es => atm%gridstruct%es
218  sin_sg => atm%gridstruct%sin_sg
219  cos_sg => atm%gridstruct%cos_sg
220  en1 => atm%gridstruct%en1
221  en2 => atm%gridstruct%en2
222 ! eww => Atm%gridstruct%eww
223 ! ess => Atm%gridstruct%ess
224 
225  sw_corner => atm%gridstruct%sw_corner
226  se_corner => atm%gridstruct%se_corner
227  ne_corner => atm%gridstruct%ne_corner
228  nw_corner => atm%gridstruct%nw_corner
229 
230  if ( atm%flagstruct%do_schmidt .and. abs(atm%flagstruct%stretch_fac-1.) > 1.e-5 ) then
231  atm%gridstruct%stretched_grid = .true.
232  symm_grid = .false.
233  else
234  atm%gridstruct%stretched_grid = .false.
235  symm_grid = .true.
236  endif
237 
238  if ( npz == 1 ) then
239  atm%ak(1) = 0.
240  atm%ak(2) = 0.
241  atm%bk(1) = 0.
242  atm%bk(2) = 1.
243  atm%ptop = 0.
244  atm%ks = 0
245  elseif ( .not. atm%flagstruct%hybrid_z ) then
246 ! Initialize (ak,bk) for cold start; overwritten with restart file
247  if (.not. atm%flagstruct%external_eta) then
248  call set_eta(npz, atm%ks, atm%ptop, atm%ak, atm%bk)
249  if ( is_master() ) then
250  write(*,*) 'Grid_init', npz, atm%ks, atm%ptop
251  tmp1 = atm%ak(atm%ks+1)
252  do k=atm%ks+1,npz
253  tmp1 = max(tmp1, (atm%ak(k)-atm%ak(k+1))/max(1.e-9, (atm%bk(k+1)-atm%bk(k))) )
254  enddo
255  write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
256  if ( tmp1 > 420.e2 ) write(*,*) 'Warning: the chosen setting in set_eta can cause instability'
257  endif
258  endif
259  endif
260 
261 ! NCEP analysis available from amip-Interp (allocate if needed)
262 #ifndef DYCORE_SOLO
263  if (.not. allocated(sst_ncep)) allocate (sst_ncep(i_sst,j_sst))
264  if (.not. allocated(sst_anom)) allocate (sst_anom(i_sst,j_sst))
265 #endif
266 
267 
268  cos_sg(:,:,:) = big_number
269  sin_sg(:,:,:) = tiny_number
270 
271  sw_corner = .false.
272  se_corner = .false.
273  ne_corner = .false.
274  nw_corner = .false.
275 
276  if (grid_type < 3 .and. .not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
277  if ( is==1 .and. js==1 ) sw_corner = .true.
278  if ( (ie+1)==npx .and. js==1 ) se_corner = .true.
279  if ( (ie+1)==npx .and. (je+1)==npy ) ne_corner = .true.
280  if ( is==1 .and. (je+1)==npy ) nw_corner = .true.
281  endif
282 
283  if ( sw_corner ) then
284  tmp1 = great_circle_dist(grid(1,1,1:2), agrid(1,1,1:2))
285  tmp2 = great_circle_dist(grid(1,1,1:2), agrid(2,2,1:2))
286  write(*,*) 'Corner interpolation coefficient=', tmp2/(tmp2-tmp1)
287  endif
288 
289  if (grid_type < 3) then
290 !xxx if ( .not. Atm%neststruct%nested ) then
291  if ( .not. atm%neststruct%nested .and. .not.atm%flagstruct%regional ) then
292  call fill_corners(grid(:,:,1), npx, npy, fill=xdir, bgrid=.true.)
293  call fill_corners(grid(:,:,2), npx, npy, fill=xdir, bgrid=.true.)
294  end if
295 
296  do j=jsd,jed+1
297  do i=isd,ied+1
298  call latlon2xyz(grid(i,j,1:2), grid3(1,i,j))
299  enddo
300  enddo
301 
302 
303  call get_center_vect( npx, npy, grid3, ec1, ec2, atm%bd )
304 
305 ! Fill arbitrary values in the non-existing corner regions:
306  if (.not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
307  do k=1,3
308  call fill_ghost(ec1(k,:,:), npx, npy, big_number, atm%bd)
309  call fill_ghost(ec2(k,:,:), npx, npy, big_number, atm%bd)
310  enddo
311  end if
312 
313 
314  do j=jsd,jed
315  do i=isd+1,ied
316  if ( ( (i<1 .and. j<1 ) .or. (i>npx .and. j<1 ) .or. &
317  (i>npx .and. j>(npy-1)) .or. (i<1 .and. j>(npy-1)) ) .and. .not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
318  ew(1:3,i,j,1:2) = 0.
319  else
320  call mid_pt_cart( grid(i,j,1:2), grid(i,j+1,1:2), pp)
321  if (i==1 .and. .not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
322  call latlon2xyz( agrid(i,j,1:2), p1)
323  call vect_cross(p2, pp, p1)
324  elseif(i==npx .and. .not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
325  call latlon2xyz( agrid(i-1,j,1:2), p1)
326  call vect_cross(p2, p1, pp)
327  else
328  call latlon2xyz( agrid(i-1,j,1:2), p3)
329  call latlon2xyz( agrid(i, j,1:2), p1)
330  call vect_cross(p2, p3, p1)
331  endif
332  call vect_cross(ew(1:3,i,j,1), p2, pp)
333  call normalize_vect(ew(1:3,i,j,1))
334 !---
335  call vect_cross(p1, grid3(1,i,j), grid3(1,i,j+1))
336  call vect_cross(ew(1:3,i,j,2), p1, pp)
337  call normalize_vect(ew(1:3,i,j,2))
338  endif
339  enddo
340  enddo
341 
342  do j=jsd+1,jed
343  do i=isd,ied
344  if ( ( (i<1 .and. j<1 ) .or. (i>(npx-1) .and. j<1 ) .or. &
345  (i>(npx-1) .and. j>npy) .or. (i<1 .and. j>npy) ) .and. .not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
346  es(1:3,i,j,1:2) = 0.
347  else
348  call mid_pt_cart(grid(i,j,1:2), grid(i+1,j,1:2), pp)
349  if (j==1 .and. .not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
350  call latlon2xyz( agrid(i,j,1:2), p1)
351  call vect_cross(p2, pp, p1)
352  elseif (j==npy .and. .not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
353  call latlon2xyz( agrid(i,j-1,1:2), p1)
354  call vect_cross(p2, p1, pp)
355  else
356  call latlon2xyz( agrid(i,j ,1:2), p1)
357  call latlon2xyz( agrid(i,j-1,1:2), p3)
358  call vect_cross(p2, p3, p1)
359  endif
360  call vect_cross(es(1:3,i,j,2), p2, pp)
361  call normalize_vect(es(1:3,i,j,2))
362 !---
363  call vect_cross(p3, grid3(1,i,j), grid3(1,i+1,j))
364  call vect_cross(es(1:3,i,j,1), p3, pp)
365  call normalize_vect(es(1:3,i,j,1))
366  endif
367  enddo
368  enddo
369 
370 ! 9---4---8
371 ! | |
372 ! 1 5 3
373 ! | |
374 ! 6---2---7
375 
376  do j=jsd,jed
377  do i=isd,ied
378 ! Testing using spherical formular: exact if coordinate lines are along great circles
379 ! SW corner:
380  cos_sg(i,j,6) = cos_angle( grid3(1,i,j), grid3(1,i+1,j), grid3(1,i,j+1) )
381 ! SE corner:
382  cos_sg(i,j,7) = -cos_angle( grid3(1,i+1,j), grid3(1,i,j), grid3(1,i+1,j+1) )
383 ! NE corner:
384  cos_sg(i,j,8) = cos_angle( grid3(1,i+1,j+1), grid3(1,i+1,j), grid3(1,i,j+1) )
385 ! NW corner:
386  cos_sg(i,j,9) = -cos_angle( grid3(1,i,j+1), grid3(1,i,j), grid3(1,i+1,j+1) )
387 ! Mid-points by averaging:
388 !!! cos_sg(i,j,1) = 0.5*( cos_sg(i,j,6) + cos_sg(i,j,9) )
389 !!! cos_sg(i,j,2) = 0.5*( cos_sg(i,j,6) + cos_sg(i,j,7) )
390 !!! cos_sg(i,j,3) = 0.5*( cos_sg(i,j,7) + cos_sg(i,j,8) )
391 !!! cos_sg(i,j,4) = 0.5*( cos_sg(i,j,8) + cos_sg(i,j,9) )
392 !!!!! cos_sg(i,j,5) = 0.25*(cos_sg(i,j,6)+cos_sg(i,j,7)+cos_sg(i,j,8)+cos_sg(i,j,9))
393 ! No averaging -----
394  call latlon2xyz(agrid(i,j,1:2), p3) ! righ-hand system consistent with grid3
395  call mid_pt3_cart(grid3(1,i,j), grid3(1,i,j+1), p1)
396  cos_sg(i,j,1) = cos_angle( p1, p3, grid3(1,i,j+1) )
397  call mid_pt3_cart(grid3(1,i,j), grid3(1,i+1,j), p1)
398  cos_sg(i,j,2) = cos_angle( p1, grid3(1,i+1,j), p3 )
399  call mid_pt3_cart(grid3(1,i+1,j), grid3(1,i+1,j+1), p1)
400  cos_sg(i,j,3) = cos_angle( p1, p3, grid3(1,i+1,j) )
401  call mid_pt3_cart(grid3(1,i,j+1), grid3(1,i+1,j+1), p1)
402  cos_sg(i,j,4) = cos_angle( p1, grid3(1,i,j+1), p3 )
403 ! Center point:
404 ! Using center_vect: [ec1, ec2]
405  cos_sg(i,j,5) = inner_prod( ec1(1:3,i,j), ec2(1:3,i,j) )
406  enddo
407  enddo
408 
409  do ip=1,9
410  do j=jsd,jed
411  do i=isd,ied
412  sin_sg(i,j,ip) = min(1.0, sqrt( max(0., 1.-cos_sg(i,j,ip)**2) ) )
413  enddo
414  enddo
415  enddo
416 
417 ! -------------------------------
418 ! For transport operation
419 ! -------------------------------
420 !xxx if (.not. Atm%neststruct%nested) then
421  if (.not. atm%neststruct%nested .and. .not. atm%flagstruct%regional) then
422  if ( sw_corner ) then
423  do i=-2,0
424  sin_sg(0,i,3) = sin_sg(i,1,2)
425  sin_sg(i,0,4) = sin_sg(1,i,1)
426  enddo
427  endif
428  if ( nw_corner ) then
429  do i=npy,npy+2
430  sin_sg(0,i,3) = sin_sg(npy-i,npy-1,4)
431  enddo
432  do i=-2,0
433  sin_sg(i,npy,2) = sin_sg(1,npx+i,1)
434  enddo
435  endif
436  if ( se_corner ) then
437  do j=-2,0
438  sin_sg(npx,j,1) = sin_sg(npx-j,1,2)
439  enddo
440  do i=npx,npx+2
441  sin_sg(i,0,4) = sin_sg(npx-1,npx-i,3)
442  enddo
443  endif
444  if ( ne_corner ) then
445  do i=npy,npy+2
446  sin_sg(npx,i,1) = sin_sg(i,npy-1,4)
447  sin_sg(i,npy,2) = sin_sg(npx-1,i,3)
448  enddo
449  endif
450  endif
451 
452 ! For AAM correction:
453  do j=js,je
454  do i=is,ie+1
455  pp1(:) = grid(i ,j ,1:2)
456  pp2(:) = grid(i,j+1 ,1:2)
457  call mid_pt_sphere(pp1, pp2, pp3)
458  call get_unit_vect2(pp1, pp2, e2)
459  call get_latlon_vector(pp3, ex, ey)
460  atm%gridstruct%l2c_v(i,j) = cos(pp3(2)) * inner_prod(e2, ex)
461  enddo
462  enddo
463  do j=js,je+1
464  do i=is,ie
465  pp1(:) = grid(i, j,1:2)
466  pp2(:) = grid(i+1,j,1:2)
467  call mid_pt_sphere(pp1, pp2, pp3)
468  call get_unit_vect2(pp1, pp2, e1)
469  call get_latlon_vector(pp3, ex, ey)
470  atm%gridstruct%l2c_u(i,j) = cos(pp3(2)) * inner_prod(e1, ex)
471  enddo
472  enddo
473 
474  else
475  cos_sg(:,:,:) = 0.
476  sin_sg(:,:,:) = 1.
477 
478  ec1(1,:,:)=1.
479  ec1(2,:,:)=0.
480  ec1(3,:,:)=0.
481 
482  ec2(1,:,:)=0.
483  ec2(2,:,:)=1.
484  ec2(3,:,:)=0.
485 
486  ew(1,:,:,1)=1.
487  ew(2,:,:,1)=0.
488  ew(3,:,:,1)=0.
489 
490  ew(1,:,:,2)=0.
491  ew(2,:,:,2)=1.
492  ew(3,:,:,2)=0.
493 
494  es(1,:,:,1)=1.
495  es(2,:,:,1)=0.
496  es(3,:,:,1)=0.
497 
498  es(1,:,:,2)=0.
499  es(2,:,:,2)=1.
500  es(3,:,:,2)=0.
501  endif
502 
503  if ( non_ortho ) then
504  cosa_u = big_number
505  cosa_v = big_number
506  cosa_s = big_number
507  sina_u = big_number
508  sina_v = big_number
509  rsin_u = big_number
510  rsin_v = big_number
511  rsina = big_number
512  rsin2 = big_number
513  cosa = big_number
514  sina = big_number
515 
516  do j=js,je+1
517  do i=is,ie+1
518 ! unit vect in X-dir: ee1
519  if (i==1 .and. .not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
520  call vect_cross(pp, grid3(1,i, j), grid3(1,i+1,j))
521  elseif(i==npx .and. .not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
522  call vect_cross(pp, grid3(1,i-1,j), grid3(1,i, j))
523  else
524  call vect_cross(pp, grid3(1,i-1,j), grid3(1,i+1,j))
525  endif
526  call vect_cross(ee1(1:3,i,j), pp, grid3(1:3,i,j))
527  call normalize_vect( ee1(1:3,i,j) )
528 
529 ! unit vect in Y-dir: ee2
530  if (j==1 .and. .not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
531  call vect_cross(pp, grid3(1:3,i,j ), grid3(1:3,i,j+1))
532  elseif(j==npy .and. .not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
533  call vect_cross(pp, grid3(1:3,i,j-1), grid3(1:3,i,j ))
534  else
535  call vect_cross(pp, grid3(1:3,i,j-1), grid3(1:3,i,j+1))
536  endif
537  call vect_cross(ee2(1:3,i,j), pp, grid3(1:3,i,j))
538  call normalize_vect( ee2(1:3,i,j) )
539 
540 ! symmetrical grid
541 #ifdef TEST_FP
542  tmp1 = inner_prod(ee1(1:3,i,j), ee2(1:3,i,j))
543  cosa(i,j) = sign(min(1., abs(tmp1)), tmp1)
544  sina(i,j) = sqrt(max(0.,1. -cosa(i,j)**2))
545 #else
546  cosa(i,j) = 0.5*(cos_sg(i-1,j-1,8)+cos_sg(i,j,6))
547  sina(i,j) = 0.5*(sin_sg(i-1,j-1,8)+sin_sg(i,j,6))
548 #endif
549  enddo
550  enddo
551 
552 ! 9---4---8
553 ! | |
554 ! 1 5 3
555 ! | |
556 ! 6---2---7
557  do j=jsd,jed
558  do i=isd+1,ied
559  cosa_u(i,j) = 0.5*(cos_sg(i-1,j,3)+cos_sg(i,j,1))
560  sina_u(i,j) = 0.5*(sin_sg(i-1,j,3)+sin_sg(i,j,1))
561 ! rsin_u(i,j) = 1. / sina_u(i,j)**2
562  rsin_u(i,j) = 1. / max(tiny_number, sina_u(i,j)**2)
563  enddo
564  enddo
565  do j=jsd+1,jed
566  do i=isd,ied
567  cosa_v(i,j) = 0.5*(cos_sg(i,j-1,4)+cos_sg(i,j,2))
568  sina_v(i,j) = 0.5*(sin_sg(i,j-1,4)+sin_sg(i,j,2))
569 ! rsin_v(i,j) = 1. / sina_v(i,j)**2
570  rsin_v(i,j) = 1. / max(tiny_number, sina_v(i,j)**2)
571  enddo
572  enddo
573 
574  do j=jsd,jed
575  do i=isd,ied
576  cosa_s(i,j) = cos_sg(i,j,5)
577 ! rsin2(i,j) = 1. / sin_sg(i,j,5)**2
578  rsin2(i,j) = 1. / max(tiny_number, sin_sg(i,j,5)**2)
579  enddo
580  enddo
581 ! Force the model to fail if incorrect corner values are to be used:
582  if (.not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
583  call fill_ghost(cosa_s, npx, npy, big_number, atm%bd)
584  end if
585 !------------------------------------
586 ! Set special sin values at edges:
587 !------------------------------------
588  do j=js,je+1
589  do i=is,ie+1
590  if ( i==npx .and. j==npy .and. .not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
591  else if ( ( i==1 .or. i==npx .or. j==1 .or. j==npy ) .and. .not. (atm%neststruct%nested.or.atm%flagstruct%regional) ) then
592  rsina(i,j) = big_number
593  else
594 ! rsina(i,j) = 1. / sina(i,j)**2
595  rsina(i,j) = 1. / max(tiny_number, sina(i,j)**2)
596  endif
597  enddo
598  enddo
599 
600  do j=jsd,jed
601  do i=is,ie+1
602  if ( (i==1 .or. i==npx) .and. .not. (atm%neststruct%nested .or. atm%flagstruct%regional) ) then
603 ! rsin_u(i,j) = 1. / sina_u(i,j)
604  rsin_u(i,j) = 1. / sign(max(tiny_number,abs(sina_u(i,j))), sina_u(i,j))
605  endif
606  enddo
607  enddo
608 
609  do j=js,je+1
610  do i=isd,ied
611  if ( (j==1 .or. j==npy) .and. .not. (atm%neststruct%nested .or. atm%flagstruct%regional) ) then
612 ! rsin_v(i,j) = 1. / sina_v(i,j)
613  rsin_v(i,j) = 1. / sign(max(tiny_number,abs(sina_v(i,j))), sina_v(i,j))
614  endif
615  enddo
616  enddo
617 
618  !EXPLANATION HERE: calling fill_ghost overwrites **SOME** of the sin_sg
619  !values along the outward-facing edge of a tile in the corners, which is incorrect.
620  !What we will do is call fill_ghost and then fill in the appropriate values
621 
622  if (.not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
623  do k=1,9
624  call fill_ghost(sin_sg(:,:,k), npx, npy, tiny_number, atm%bd) ! this will cause NAN if used
625  call fill_ghost(cos_sg(:,:,k), npx, npy, big_number, atm%bd)
626  enddo
627  end if
628 
629 ! -------------------------------
630 ! For transport operation
631 ! -------------------------------
632  if ( sw_corner ) then
633  do i=0,-2,-1
634  sin_sg(0,i,3) = sin_sg(i,1,2)
635  sin_sg(i,0,4) = sin_sg(1,i,1)
636  cos_sg(0,i,3) = cos_sg(i,1,2)
637  cos_sg(i,0,4) = cos_sg(1,i,1)
638 !!! cos_sg(0,i,7) = cos_sg(i,1,6)
639 !!! cos_sg(0,i,8) = cos_sg(i,1,7)
640 !!! cos_sg(i,0,8) = cos_sg(1,i,9)
641 !!! cos_sg(i,0,9) = cos_sg(1,i,6)
642  enddo
643 !!! cos_sg(0,0,8) = 0.5*(cos_sg(0,1,7)+cos_sg(1,0,9))
644 
645  endif
646  if ( nw_corner ) then
647  do i=npy,npy+2
648  sin_sg(0,i,3) = sin_sg(npy-i,npy-1,4)
649  cos_sg(0,i,3) = cos_sg(npy-i,npy-1,4)
650 !!! cos_sg(0,i,7) = cos_sg(npy-i,npy-1,8)
651 !!! cos_sg(0,i,8) = cos_sg(npy-i,npy-1,9)
652  enddo
653  do i=0,-2,-1
654  sin_sg(i,npy,2) = sin_sg(1,npy-i,1)
655  cos_sg(i,npy,2) = cos_sg(1,npy-i,1)
656 !!! cos_sg(i,npy,6) = cos_sg(1,npy-i,9)
657 !!! cos_sg(i,npy,7) = cos_sg(1,npy-i,6)
658  enddo
659 !!! cos_sg(0,npy,7) = 0.5*(cos_sg(1,npy,6)+cos_sg(0,npy-1,8))
660  endif
661  if ( se_corner ) then
662  do j=0,-2,-1
663  sin_sg(npx,j,1) = sin_sg(npx-j,1,2)
664  cos_sg(npx,j,1) = cos_sg(npx-j,1,2)
665 !!! cos_sg(npx,j,6) = cos_sg(npx-j,1,7)
666 !!! cos_sg(npx,j,9) = cos_sg(npx-j,1,6)
667  enddo
668  do i=npx,npx+2
669  sin_sg(i,0,4) = sin_sg(npx-1,npx-i,3)
670  cos_sg(i,0,4) = cos_sg(npx-1,npx-i,3)
671 !!! cos_sg(i,0,9) = cos_sg(npx-1,npx-i,8)
672 !!! cos_sg(i,0,8) = cos_sg(npx-1,npx-i,7)
673  enddo
674 !!! cos_sg(npx,0,9) = 0.5*(cos_sg(npx,1,6)+cos_sg(npx-1,0,8))
675  endif
676  if ( ne_corner ) then
677  do i=0,2
678  sin_sg(npx,npy+i,1) = sin_sg(npx+i,npy-1,4)
679  sin_sg(npx+i,npy,2) = sin_sg(npx-1,npy+i,3)
680  cos_sg(npx,npy+i,1) = cos_sg(npx+i,npy-1,4)
681 !!! cos_sg(npx,npy+i,6) = cos_sg(npx+i,npy-1,9)
682 !!! cos_sg(npx,npy+i,9) = cos_sg(npx+i,npy-1,8)
683  cos_sg(npx+i,npy,2) = cos_sg(npx-1,npy+i,3)
684 !!! cos_sg(npx+i,npy,6) = cos_sg(npx-1,npy+i,7)
685 !!! cos_sg(npx+i,npy,7) = cos_sg(npx-1,npy+i,8)
686  end do
687 !!! cos_sg(npx,npy,6) = 0.5*(cos_sg(npx-1,npy,7)+cos_sg(npx,npy-1,9))
688  endif
689 
690  else
691  sina = 1.
692  cosa = 0.
693  rsina = 1.
694  rsin2 = 1.
695  sina_u = 1.
696  sina_v = 1.
697  cosa_u = 0.
698  cosa_v = 0.
699  cosa_s = 0.
700  rsin_u = 1.
701  rsin_v = 1.
702  endif
703 
704  if ( grid_type < 3 ) then
705 
706 #ifdef USE_NORM_VECT
707 !-------------------------------------------------------------
708 ! Make normal vect at face edges after consines are computed:
709 !-------------------------------------------------------------
710 ! for old d2a2c_vect routines
711  if (.not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
712  do j=js-1,je+1
713  if ( is==1 ) then
714  i=1
715  call vect_cross(ew(1,i,j,1), grid3(1,i,j+1), grid3(1,i,j))
716  call normalize_vect( ew(1,i,j,1) )
717  endif
718  if ( (ie+1)==npx ) then
719  i=npx
720  call vect_cross(ew(1,i,j,1), grid3(1,i,j+1), grid3(1,i,j))
721  call normalize_vect( ew(1,i,j,1) )
722  endif
723  enddo
724 
725  if ( js==1 ) then
726  j=1
727  do i=is-1,ie+1
728  call vect_cross(es(1,i,j,2), grid3(1,i,j),grid3(1,i+1,j))
729  call normalize_vect( es(1,i,j,2) )
730  enddo
731  endif
732  if ( (je+1)==npy ) then
733  j=npy
734  do i=is-1,ie+1
735  call vect_cross(es(1,i,j,2), grid3(1,i,j),grid3(1,i+1,j))
736  call normalize_vect( es(1,i,j,2) )
737  enddo
738  endif
739  endif
740 #endif
741 
742 ! For omega computation:
743 ! Unit vectors:
744  do j=js,je+1
745  do i=is,ie
746  call vect_cross(en1(1:3,i,j), grid3(1,i,j), grid3(1,i+1,j))
747  call normalize_vect( en1(1:3,i,j) )
748  enddo
749  enddo
750  do j=js,je
751  do i=is,ie+1
752  call vect_cross(en2(1:3,i,j), grid3(1,i,j+1), grid3(1,i,j))
753  call normalize_vect( en2(1:3,i,j) )
754  enddo
755  enddo
756 !-------------------------------------------------------------
757 ! Make unit vectors for the coordinate extension:
758 !-------------------------------------------------------------
759  endif
760 
761 !xxxx
762 !!! should we insert .not.regional into the following loops alongside .not.nested ????
763 !xxxx
764  do j=jsd,jed+1
765  if ((j==1 .OR. j==npy) .and. .not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
766  do i=isd,ied
767  divg_u(i,j) = 0.5*(sin_sg(i,j,2)+sin_sg(i,j-1,4))*dyc(i,j)/dx(i,j)
768  del6_u(i,j) = 0.5*(sin_sg(i,j,2)+sin_sg(i,j-1,4))*dx(i,j)/dyc(i,j)
769  enddo
770  else
771  do i=isd,ied
772  divg_u(i,j) = sina_v(i,j)*dyc(i,j)/dx(i,j)
773  del6_u(i,j) = sina_v(i,j)*dx(i,j)/dyc(i,j)
774  enddo
775  end if
776  enddo
777  do j=jsd,jed
778  do i=isd,ied+1
779  divg_v(i,j) = sina_u(i,j)*dxc(i,j)/dy(i,j)
780  del6_v(i,j) = sina_u(i,j)*dy(i,j)/dxc(i,j)
781  enddo
782  if (is == 1 .and. .not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
783  divg_v(is,j) = 0.5*(sin_sg(1,j,1)+sin_sg(0,j,3))*dxc(is,j)/dy(is,j)
784  del6_v(is,j) = 0.5*(sin_sg(1,j,1)+sin_sg(0,j,3))*dy(is,j)/dxc(is,j)
785  endif
786  if (ie+1 == npx .and. .not. (atm%neststruct%nested .or. atm%flagstruct%regional)) then
787  divg_v(ie+1,j) = 0.5*(sin_sg(npx,j,1)+sin_sg(npx-1,j,3))*dxc(ie+1,j)/dy(ie+1,j)
788  del6_v(ie+1,j) = 0.5*(sin_sg(npx,j,1)+sin_sg(npx-1,j,3))*dy(ie+1,j)/dxc(ie+1,j)
789  endif
790  enddo
791 
792 ! Initialize cubed_sphere to lat-lon transformation:
793  call init_cubed_to_latlon( atm%gridstruct, atm%flagstruct%hydrostatic, agrid, grid_type, c2l_order, atm%bd )
794 
795  call global_mx(area, ng, atm%gridstruct%da_min, atm%gridstruct%da_max, atm%bd)
796  if( is_master() ) write(*,*) 'da_max/da_min=', atm%gridstruct%da_max/atm%gridstruct%da_min
797 
798  call global_mx_c(area_c(is:ie,js:je), is, ie, js, je, atm%gridstruct%da_min_c, atm%gridstruct%da_max_c)
799 
800  if( is_master() ) write(*,*) 'da_max_c/da_min_c=', atm%gridstruct%da_max_c/atm%gridstruct%da_min_c
801 
802 !------------------------------------------------
803 ! Initialization for interpolation at face edges
804 !------------------------------------------------
805 ! A->B scalar:
806  if (grid_type < 3 .and. .not. atm%neststruct%nested .and. .not. atm%flagstruct%regional ) then
807  call mpp_update_domains(divg_v, divg_u, atm%domain, flags=scalar_pair, &
808  gridtype=cgrid_ne_param, complete=.true.)
809  call mpp_update_domains(del6_v, del6_u, atm%domain, flags=scalar_pair, &
810  gridtype=cgrid_ne_param, complete=.true.)
811  call edge_factors (atm%gridstruct%edge_s, atm%gridstruct%edge_n, atm%gridstruct%edge_w, &
812  atm%gridstruct%edge_e, non_ortho, grid, agrid, npx, npy, atm%bd)
813  call efactor_a2c_v(atm%gridstruct%edge_vect_s, atm%gridstruct%edge_vect_n, &
814  atm%gridstruct%edge_vect_w, atm%gridstruct%edge_vect_e, &
815  non_ortho, grid, agrid, npx, npy, atm%neststruct%nested, atm%bd, atm%flagstruct%regional)
816 ! call extend_cube_s(non_ortho, grid, agrid, npx, npy, .false., Atm%neststruct%nested)
817 ! call van2d_init(grid, agrid, npx, npy)
818  else
819 
820  atm%gridstruct%edge_s = big_number
821  atm%gridstruct%edge_n = big_number
822  atm%gridstruct%edge_w = big_number
823  atm%gridstruct%edge_e = big_number
824 
825  atm%gridstruct%edge_vect_s = big_number
826  atm%gridstruct%edge_vect_n = big_number
827  atm%gridstruct%edge_vect_w = big_number
828  atm%gridstruct%edge_vect_e = big_number
829 
830  endif
831 
832 !32-bit versions of the data
833  atm%gridstruct%grid = atm%gridstruct%grid_64
834  atm%gridstruct%agrid = atm%gridstruct%agrid_64
835  atm%gridstruct%area = atm%gridstruct%area_64
836  atm%gridstruct%area_c = atm%gridstruct%area_c_64
837  atm%gridstruct%dx = atm%gridstruct%dx_64
838  atm%gridstruct%dy = atm%gridstruct%dy_64
839  atm%gridstruct%dxa = atm%gridstruct%dxa_64
840  atm%gridstruct%dya = atm%gridstruct%dya_64
841  atm%gridstruct%dxc = atm%gridstruct%dxc_64
842  atm%gridstruct%dyc = atm%gridstruct%dyc_64
843  atm%gridstruct%cosa = atm%gridstruct%cosa_64
844  atm%gridstruct%sina = atm%gridstruct%sina_64
845 
846 !--- deallocate the higher-order gridstruct arrays
847 !rab deallocate ( Atm%gridstruct%grid_64 )
848 !rab deallocate ( Atm%gridstruct%agrid_64 )
849 !rab deallocate ( Atm%gridstruct%area_64 )
850  deallocate ( atm%gridstruct%area_c_64 )
851 !rab deallocate ( Atm%gridstruct%dx_64 )
852 !rab deallocate ( Atm%gridstruct%dy_64 )
853  deallocate ( atm%gridstruct%dxa_64 )
854  deallocate ( atm%gridstruct%dya_64 )
855  deallocate ( atm%gridstruct%dxc_64 )
856  deallocate ( atm%gridstruct%dyc_64 )
857  deallocate ( atm%gridstruct%cosa_64 )
858  deallocate ( atm%gridstruct%sina_64 )
859 
860  nullify(agrid)
861  nullify(grid)
862  nullify(area)
863  nullify(area_c)
864  nullify(dx)
865  nullify(dy)
866  nullify(dxc)
867  nullify(dyc)
868  nullify(dxa)
869  nullify(dya)
870  nullify(sina)
871  nullify(cosa)
872  nullify(divg_u)
873  nullify(divg_v)
874 
875  nullify(del6_u)
876  nullify(del6_v)
877 
878  nullify(cosa_u)
879  nullify(cosa_v)
880  nullify(cosa_s)
881  nullify(sina_u)
882  nullify(sina_v)
883  nullify(rsin_u)
884  nullify(rsin_v)
885  nullify(rsina)
886  nullify(rsin2)
887  nullify(ee1)
888  nullify(ee2)
889  nullify(ec1)
890  nullify(ec2)
891  nullify(ew)
892  nullify(es)
893  nullify(sin_sg)
894  nullify(cos_sg)
895  nullify(en1)
896  nullify(en2)
897  nullify(sw_corner)
898  nullify(se_corner)
899  nullify(ne_corner)
900  nullify(nw_corner)
901 
902  end subroutine grid_utils_init
903 
904 
905  subroutine grid_utils_end
906 ! deallocate sst_ncep (if allocated)
907 #ifndef DYCORE_SOLO
908  if (allocated(sst_ncep)) deallocate( sst_ncep )
909  if (allocated(sst_anom)) deallocate( sst_anom )
910 #endif
911  end subroutine grid_utils_end
912 
917  subroutine direct_transform(c, i1, i2, j1, j2, lon_p, lat_p, n, lon, lat)
918  real(kind=R_GRID), intent(in):: c
919  real(kind=R_GRID), intent(in):: lon_p, lat_p
920  integer, intent(in):: n
921  integer, intent(in):: i1, i2, j1, j2
922 ! 0 <= lon <= 2*pi ; -pi/2 <= lat <= pi/2
923  real(kind=R_GRID), intent(inout), dimension(i1:i2,j1:j2):: lon, lat
924 !
925  real(f_p):: lat_t, sin_p, cos_p, sin_lat, cos_lat, sin_o, p2, two_pi
926  real(f_p):: c2p1, c2m1
927  integer:: i, j
928 
929  p2 = 0.5d0*pi
930  two_pi = 2.d0*pi
931 
932  if( is_master() .and. n==1 ) then
933  write(*,*) n, 'Schmidt transformation: stretching factor=', c, ' center=', lon_p, lat_p
934  endif
935 
936  c2p1 = 1.d0 + c*c
937  c2m1 = 1.d0 - c*c
938 
939  sin_p = sin(lat_p)
940  cos_p = cos(lat_p)
941 
942  do j=j1,j2
943  do i=i1,i2
944  if ( abs(c2m1) > 1.d-7 ) then
945  sin_lat = sin(lat(i,j))
946  lat_t = asin( (c2m1+c2p1*sin_lat)/(c2p1+c2m1*sin_lat) )
947  else ! no stretching
948  lat_t = lat(i,j)
949  endif
950  sin_lat = sin(lat_t)
951  cos_lat = cos(lat_t)
952  sin_o = -(sin_p*sin_lat + cos_p*cos_lat*cos(lon(i,j)))
953  if ( (1.-abs(sin_o)) < 1.d-7 ) then ! poles
954  lon(i,j) = 0.d0
955  lat(i,j) = sign( p2, sin_o )
956  else
957  lat(i,j) = asin( sin_o )
958  lon(i,j) = lon_p + atan2( -cos_lat*sin(lon(i,j)), &
959  -sin_lat*cos_p+cos_lat*sin_p*cos(lon(i,j)))
960  if ( lon(i,j) < 0.d0 ) then
961  lon(i,j) = lon(i,j) + two_pi
962  elseif( lon(i,j) >= two_pi ) then
963  lon(i,j) = lon(i,j) - two_pi
964  endif
965  endif
966  enddo
967  enddo
968 
969  end subroutine direct_transform
970 
971 
972  real function inner_prod(v1, v2)
973  real(kind=R_GRID),intent(in):: v1(3), v2(3)
974  real (f_p) :: vp1(3), vp2(3), prod16
975  integer k
976 
977  do k=1,3
978  vp1(k) = real(v1(k),kind=f_p)
979  vp2(k) = real(v2(k),kind=f_p)
980  enddo
981  prod16 = vp1(1)*vp2(1) + vp1(2)*vp2(2) + vp1(3)*vp2(3)
982  inner_prod = prod16
983 
984  end function inner_prod
985 
988  subroutine efactor_a2c_v(edge_vect_s, edge_vect_n, edge_vect_w, edge_vect_e, non_ortho, grid, agrid, npx, npy, nested, bd, regional)
989  type(fv_grid_bounds_type), intent(IN) :: bd
990  real(kind=R_GRID), intent(INOUT), dimension(bd%isd:bd%ied) :: edge_vect_s, edge_vect_n
991  real(kind=R_GRID), intent(INOUT), dimension(bd%jsd:bd%jed) :: edge_vect_w, edge_vect_e
992  logical, intent(in):: non_ortho, nested, regional
993  real(kind=R_GRID), intent(in):: grid(bd%isd:bd%ied+1,bd%jsd:bd%jed+1,2)
994  real(kind=R_GRID), intent(in):: agrid(bd%isd:bd%ied ,bd%jsd:bd%jed ,2)
995  integer, intent(in):: npx, npy
996 
997  real(kind=R_GRID) px(2,bd%isd:bd%ied+1), py(2,bd%jsd:bd%jed+1)
998  real(kind=R_GRID) p1(2,bd%isd:bd%ied+1), p2(2,bd%jsd:bd%jed+1)
999  real(kind=R_GRID) d1, d2
1000  integer i, j
1001  integer im2, jm2
1002 
1003  integer :: is, ie, js, je
1004  integer :: isd, ied, jsd, jed
1005 
1006  is = bd%is
1007  ie = bd%ie
1008  js = bd%js
1009  je = bd%je
1010  isd = bd%isd
1011  ied = bd%ied
1012  jsd = bd%jsd
1013  jed = bd%jed
1014 
1015 
1016  if ( .not. non_ortho ) then
1017  edge_vect_s = 0.
1018  edge_vect_n = 0.
1019  edge_vect_w = 0.
1020  edge_vect_e = 0.
1021  else
1022  edge_vect_s = big_number
1023  edge_vect_n = big_number
1024  edge_vect_w = big_number
1025  edge_vect_e = big_number
1026 
1027  if ( npx /= npy .and. .not. (nested .or. regional)) call mpp_error(fatal, 'efactor_a2c_v: npx /= npy')
1028  if ( (npx/2)*2 == npx ) call mpp_error(fatal, 'efactor_a2c_v: npx/npy is not an odd number')
1029 
1030  im2 = (npx-1)/2
1031  jm2 = (npy-1)/2
1032 
1033  if ( is==1 ) then
1034  i=1
1035  do j=js-2,je+2
1036  call mid_pt_sphere(agrid(i-1,j,1:2), agrid(i,j, 1:2), py(1,j))
1037  call mid_pt_sphere( grid(i, j,1:2), grid(i,j+1,1:2), p2(1,j))
1038  enddo
1039 
1040 ! west edge:
1041 !------------------------------------------------------------------
1042 ! v_sw(j) = (1.-edge_vect_w(j)) * p(j) + edge_vect_w(j) * p(j+1)
1043 !------------------------------------------------------------------
1044  do j=js-1,je+1
1045  if ( j<=jm2 ) then
1046  d1 = great_circle_dist( py(1,j ), p2(1,j) )
1047  d2 = great_circle_dist( py(1,j+1), p2(1,j) )
1048  edge_vect_w(j) = d1 / ( d1 + d2 )
1049  else
1050  d2 = great_circle_dist( py(1,j-1), p2(1,j) )
1051  d1 = great_circle_dist( py(1,j ), p2(1,j) )
1052  edge_vect_w(j) = d1 / ( d2 + d1 )
1053  endif
1054  enddo
1055  if ( js==1 ) then
1056  edge_vect_w(0) = edge_vect_w(1)
1057  endif
1058  if ( (je+1)==npy ) then
1059  edge_vect_w(npy) = edge_vect_w(je)
1060  endif
1061  do j=js-1,je+1
1062 ! if ( is_master() ) write(*,*) j, edge_vect_w(j)
1063  enddo
1064  endif
1065 
1066  if ( (ie+1)==npx ) then
1067  i=npx
1068  do j=jsd,jed
1069  call mid_pt_sphere(agrid(i-1,j,1:2), agrid(i,j, 1:2), py(1,j))
1070  call mid_pt_sphere( grid(i, j,1:2), grid(i,j+1,1:2), p2(1,j))
1071  enddo
1072 
1073  do j=js-1,je+1
1074  if ( j<=jm2 ) then
1075  d1 = great_circle_dist( py(1,j ), p2(1,j) )
1076  d2 = great_circle_dist( py(1,j+1), p2(1,j) )
1077  edge_vect_e(j) = d1 / ( d1 + d2 )
1078  else
1079  d2 = great_circle_dist( py(1,j-1), p2(1,j) )
1080  d1 = great_circle_dist( py(1,j ), p2(1,j) )
1081  edge_vect_e(j) = d1 / ( d2 + d1 )
1082  endif
1083  enddo
1084  if ( js==1 ) then
1085  edge_vect_e(0) = edge_vect_e(1)
1086  endif
1087  if ( (je+1)==npy ) then
1088  edge_vect_e(npy) = edge_vect_e(je)
1089  endif
1090  do j=js-1,je+1
1091 ! if ( is_master() ) write(*,*) j, edge_vect_e(j)
1092  enddo
1093  endif
1094 
1095  if ( js==1 ) then
1096  j=1
1097  do i=isd,ied
1098  call mid_pt_sphere(agrid(i,j-1,1:2), agrid(i, j,1:2), px(1,i))
1099  call mid_pt_sphere( grid(i,j, 1:2), grid(i+1,j,1:2), p1(1,i))
1100  enddo
1101 ! south_west edge:
1102 !------------------------------------------------------------------
1103 ! v_s(i) = (1.-edge_vect_s(i)) * p(i) + edge_vect_s(i) * p(i+1)
1104 !------------------------------------------------------------------
1105  do i=is-1,ie+1
1106  if ( i<=im2 ) then
1107  d1 = great_circle_dist( px(1,i ), p1(1,i) )
1108  d2 = great_circle_dist( px(1,i+1), p1(1,i) )
1109  edge_vect_s(i) = d1 / ( d1 + d2 )
1110  else
1111  d2 = great_circle_dist( px(1,i-1), p1(1,i) )
1112  d1 = great_circle_dist( px(1,i ), p1(1,i) )
1113  edge_vect_s(i) = d1 / ( d2 + d1 )
1114  endif
1115  enddo
1116  if ( is==1 ) then
1117  edge_vect_s(0) = edge_vect_s(1)
1118  endif
1119  if ( (ie+1)==npx ) then
1120  edge_vect_s(npx) = edge_vect_s(ie)
1121  endif
1122  do i=is-1,ie+1
1123 ! if ( is_master() ) write(*,*) i, edge_vect_s(i)
1124  enddo
1125  endif
1126 
1127 
1128  if ( (je+1)==npy ) then
1129 ! v_n(i) = (1.-edge_vect_n(i)) * p(i) + edge_vect_n(i) * p(i+1)
1130  j=npy
1131  do i=isd,ied
1132  call mid_pt_sphere(agrid(i,j-1,1:2), agrid(i, j,1:2), px(1,i))
1133  call mid_pt_sphere( grid(i,j, 1:2), grid(i+1,j,1:2), p1(1,i))
1134  enddo
1135 
1136  do i=is-1,ie+1
1137  if ( i<=im2 ) then
1138  d1 = great_circle_dist( px(1,i ), p1(1,i) )
1139  d2 = great_circle_dist( px(1,i+1), p1(1,i) )
1140  edge_vect_n(i) = d1 / ( d1 + d2 )
1141  else
1142  d2 = great_circle_dist( px(1,i-1), p1(1,i) )
1143  d1 = great_circle_dist( px(1,i ), p1(1,i) )
1144  edge_vect_n(i) = d1 / ( d2 + d1 )
1145  endif
1146  enddo
1147  if ( is==1 ) then
1148  edge_vect_n(0) = edge_vect_n(1)
1149  endif
1150  if ( (ie+1)==npx ) then
1151  edge_vect_n(npx) = edge_vect_n(ie)
1152  endif
1153  do i=is-1,ie+1
1154 ! if ( is_master() ) write(*,*) i, edge_vect_n(i)
1155  enddo
1156  endif
1157 
1158  endif
1159 
1160  end subroutine efactor_a2c_v
1161 
1164 ! Sets up edge_?
1165  subroutine edge_factors(edge_s, edge_n, edge_w, edge_e, non_ortho, grid, agrid, npx, npy, bd)
1166  type(fv_grid_bounds_type), intent(IN) :: bd
1167  real(kind=R_GRID), intent(INOUT), dimension(npx) :: edge_s, edge_n
1168  real(kind=R_GRID), intent(INOUT), dimension(npy) :: edge_w, edge_e
1169  logical, intent(in):: non_ortho
1170  real(kind=R_GRID), intent(in):: grid(bd%isd:bd%ied+1,bd%jsd:bd%jed+1,2)
1171  real(kind=R_GRID), intent(in):: agrid(bd%isd:bd%ied ,bd%jsd:bd%jed ,2)
1172  integer, intent(in):: npx, npy
1173 
1174  real(kind=R_GRID) px(2,npx), py(2,npy)
1175  real(kind=R_GRID) d1, d2
1176  integer i, j
1177 
1178  integer :: is, ie, js, je
1179  integer :: isd, ied, jsd, jed
1180 
1181  is = bd%is
1182  ie = bd%ie
1183  js = bd%js
1184  je = bd%je
1185  isd = bd%isd
1186  ied = bd%ied
1187  jsd = bd%jsd
1188 
1189 
1190  if ( .not. non_ortho ) then
1191  edge_s = 0.5d0
1192  edge_n = 0.5d0
1193  edge_w = 0.5d0
1194  edge_e = 0.5d0
1195  else
1196  edge_s = big_number
1197  edge_n = big_number
1198  edge_w = big_number
1199  edge_e = big_number
1200 
1201 ! west edge:
1202 !----------------------------------------------------------
1203 ! p_west(j) = (1.-edge_w(j)) * p(j) + edge_w(j) * p(j-1)
1204 !----------------------------------------------------------
1205  if ( is==1 ) then
1206  i=1
1207  do j=max(1,js-1), min(npy-1,je+1)
1208  call mid_pt_sphere(agrid(i-1,j,1:2), agrid(i,j,1:2), py(1,j))
1209  enddo
1210  do j=max(2,js), min(npy-1,je+1)
1211  d1 = great_circle_dist( py(1,j-1), grid(i,j,1:2) )
1212  d2 = great_circle_dist( py(1,j ), grid(i,j,1:2) )
1213  edge_w(j) = d2 / ( d1 + d2 )
1214  enddo
1215  endif
1216 
1217 ! east edge:
1218 !----------------------------------------------------------
1219 ! p_east(j) = (1.-edge_e(j)) * p(j) + edge_e(j) * p(j-1)
1220 !----------------------------------------------------------
1221  if ( (ie+1)==npx ) then
1222  i=npx
1223  do j=max(1,js-1), min(npy-1,je+1)
1224  call mid_pt_sphere(agrid(i-1,j,1:2), agrid(i,j,1:2), py(1,j))
1225  enddo
1226  do j=max(2,js), min(npy-1,je+1)
1227  d1 = great_circle_dist( py(1,j-1), grid(i,j,1:2) )
1228  d2 = great_circle_dist( py(1,j ), grid(i,j,1:2) )
1229  edge_e(j) = d2 / ( d1 + d2 )
1230 ! Check rounding difference:
1231 ! if(is_master()) write(*,*) j, edge_w(j) - edge_e(j)
1232  enddo
1233  endif
1234 
1235 
1236 ! south edge:
1237 !----------------------------------------------------------
1238 ! p_south(j) = (1.-edge_s(i)) * p(i) + edge_s(i) * p(i-1)
1239 !----------------------------------------------------------
1240  if ( js==1 ) then
1241  j=1
1242  do i=max(1,is-1), min(npx-1,ie+1)
1243  call mid_pt_sphere(agrid(i,j-1,1:2), agrid(i,j,1:2), px(1,i))
1244  enddo
1245  do i=max(2,is), min(npx-1,ie+1)
1246  d1 = great_circle_dist( px(1,i-1), grid(i,j,1:2) )
1247  d2 = great_circle_dist( px(1,i ), grid(i,j,1:2) )
1248  edge_s(i) = d2 / ( d1 + d2 )
1249  enddo
1250  endif
1251 
1252 ! North edge:
1253 !----------------------------------------------------------
1254 ! p_north(j) = (1.-edge_n(i)) * p(i) + edge_n(i) * p(i-1)
1255 !----------------------------------------------------------
1256  if ( (je+1)==npy ) then
1257  j=npy
1258  do i=max(1,is-1), min(npx-1,ie+1)
1259  call mid_pt_sphere(agrid(i,j-1,1:2), agrid(i,j,1:2), px(1,i))
1260  enddo
1261  do i=max(2,is), min(npx-1,ie+1)
1262  d1 = great_circle_dist( px(1,i-1), grid(i,j,1:2) )
1263  d2 = great_circle_dist( px(1,i ), grid(i,j,1:2) )
1264  edge_n(i) = d2 / ( d1 + d2 )
1265 ! if(is_master()) write(*,*) i, edge_s(i), edge_n(i)-edge_s(i)
1266  enddo
1267  endif
1268  endif
1269 
1270  end subroutine edge_factors
1271 
1272 
1273  subroutine gnomonic_grids(grid_type, im, lon, lat)
1274  integer, intent(in):: im, grid_type
1275  real(kind=R_GRID), intent(out):: lon(im+1,im+1)
1276  real(kind=R_GRID), intent(out):: lat(im+1,im+1)
1277  integer i, j
1278 
1279  if(grid_type==0) call gnomonic_ed( im, lon, lat)
1280  if(grid_type==1) call gnomonic_dist(im, lon, lat)
1281  if(grid_type==2) call gnomonic_angl(im, lon, lat)
1282 
1283 
1284  if(grid_type<3) then
1285  call symm_ed(im, lon, lat)
1286  do j=1,im+1
1287  do i=1,im+1
1288  lon(i,j) = lon(i,j) - pi
1289  enddo
1290  enddo
1291 ! call van2_init(lon, lat, im+1, im+1)
1292  endif
1293 
1294  end subroutine gnomonic_grids
1295 
1296 !-----------------------------------------------------
1299 !-----------------------------------------------------
1308  subroutine gnomonic_ed(im, lamda, theta)
1309 
1310  integer, intent(in):: im
1311  real(kind=R_GRID), intent(out):: lamda(im+1,im+1)
1312  real(kind=R_GRID), intent(out):: theta(im+1,im+1)
1313 
1314 ! Local:
1315  real(kind=R_GRID) pp(3,im+1,im+1)
1316  real(kind=R_GRID) p1(2), p2(2)
1317  real(f_p):: rsq3, alpha, delx, dely
1318  integer i, j, k
1319 
1320  rsq3 = 1.d0/sqrt(3.d0)
1321  alpha = asin( rsq3 )
1322 
1323 ! Ranges:
1324 ! lamda = [0.75*pi, 1.25*pi]
1325 ! theta = [-alpha, alpha]
1326 
1327  dely = 2.d0*alpha / real(im,kind=f_p)
1328 
1329 ! Define East-West edges:
1330  do j=1,im+1
1331  lamda(1, j) = 0.75d0*pi ! West edge
1332  lamda(im+1,j) = 1.25d0*pi ! East edge
1333  theta(1, j) = -alpha + dely*real(j-1,kind=f_p) ! West edge
1334  theta(im+1,j) = theta(1,j) ! East edge
1335  enddo
1336 
1337 ! Get North-South edges by symmetry:
1338 
1339  do i=2,im
1340  call mirror_latlon(lamda(1,1), theta(1,1), lamda(im+1,im+1), theta(im+1,im+1), &
1341  lamda(1,i), theta(1,i), lamda(i,1), theta(i, 1) )
1342  lamda(i,im+1) = lamda(i,1)
1343  theta(i,im+1) = -theta(i,1)
1344  enddo
1345 
1346 ! Set 4 corners:
1347  call latlon2xyz2(lamda(1 , 1), theta(1, 1), pp(1, 1, 1))
1348  call latlon2xyz2(lamda(im+1, 1), theta(im+1, 1), pp(1,im+1, 1))
1349  call latlon2xyz2(lamda(1, im+1), theta(1, im+1), pp(1, 1,im+1))
1350  call latlon2xyz2(lamda(im+1,im+1), theta(im+1,im+1), pp(1,im+1,im+1))
1351 
1352 ! Map edges on the sphere back to cube:
1353 ! Intersections at x=-rsq3
1354 
1355  i=1
1356  do j=2,im
1357  call latlon2xyz2(lamda(i,j), theta(i,j), pp(1,i,j))
1358  pp(2,i,j) = -pp(2,i,j)*rsq3/pp(1,i,j)
1359  pp(3,i,j) = -pp(3,i,j)*rsq3/pp(1,i,j)
1360  enddo
1361 
1362  j=1
1363  do i=2,im
1364  call latlon2xyz2(lamda(i,j), theta(i,j), pp(1,i,1))
1365  pp(2,i,1) = -pp(2,i,1)*rsq3/pp(1,i,1)
1366  pp(3,i,1) = -pp(3,i,1)*rsq3/pp(1,i,1)
1367  enddo
1368 
1369  do j=1,im+1
1370  do i=1,im+1
1371  pp(1,i,j) = -rsq3
1372  enddo
1373  enddo
1374 
1375  do j=2,im+1
1376  do i=2,im+1
1377 ! Copy y-z face of the cube along j=1
1378  pp(2,i,j) = pp(2,i,1)
1379 ! Copy along i=1
1380  pp(3,i,j) = pp(3,1,j)
1381  enddo
1382  enddo
1383 
1384  call cart_to_latlon( (im+1)*(im+1), pp, lamda, theta)
1385 
1386 ! Compute great-circle-distance "resolution" along the face edge:
1387  if ( is_master() ) then
1388  p1(1) = lamda(1,1); p1(2) = theta(1,1)
1389  p2(1) = lamda(2,1); p2(2) = theta(2,1)
1390  write(*,*) 'Gird distance at face edge (km)=',great_circle_dist( p1, p2, radius ) ! earth radius is assumed
1391  endif
1392 
1393  end subroutine gnomonic_ed
1394 
1404  subroutine gnomonic_ed_limited(im, in, nghost, lL, lR, uL, uR, lamda, theta)
1405  integer, intent(IN) :: im, in, nghost
1406  real(kind=R_GRID), intent(IN), dimension(2) :: lL, lR, uL, uR
1407  real(kind=R_GRID), intent(OUT) :: lamda(1-nghost:im+1+nghost,1-nghost:in+1+nghost)
1408  real(kind=R_GRID), intent(OUT) :: theta(1-nghost:im+1+nghost,1-nghost:in+1+nghost)
1409 
1410  ! Local:
1411  real(kind=R_GRID) pp(3,1-nghost:im+1+nghost,1-nghost:in+1+nghost)
1412  real(kind=R_GRID) p1(2), p2(2)
1413  real(f_p):: rsq3, alpha, delx, dely
1414  integer i, j, k, irefl
1415 
1416  rsq3 = 1.d0/sqrt(3.d0)
1417  alpha = asin( rsq3 )
1418 
1419  lamda(1,1) = ll(1); theta(1,1) = ll(2)
1420  lamda(im+1,1) = lr(1); theta(im+1,1) = lr(2)
1421  lamda(1,in+1) = ul(1); theta(1,in+1) = ul(2)
1422  lamda(im+1,in+1) = ur(1); theta(im+1,in+1) = ur(2)
1423 
1424  !Since meridians are great circles, grid spacing is equidistant in
1425  !lat-lon space along the east and west edges of the grid
1426  dely = (ul(2) - ll(2))/in
1427  do j=2,in+1+nghost
1428  theta(1,j) = theta(1,j-1) + dely
1429  theta(in+1,j) = theta(in+1,j-1) + dely
1430  lamda(1,j) = lamda(1,1)
1431  lamda(in+1,j) = lamda(in+1,1)
1432  end do
1433  do j=0,1-nghost,-1
1434  theta(1,j) = theta(1,j+1) - dely
1435  theta(in+1,j) = theta(in+1,j+1) - dely
1436  lamda(1,j) = lamda(1,1)
1437  lamda(in+1,j) = lamda(in+1,1)
1438  end do
1439 
1440  lamda(1,:) = lamda(1,1)
1441  lamda(in+1,:) = lamda(in+1,1)
1442 
1443  !Here, instead of performing a reflection (as in gnomonic_ed) to get the north and south
1444  !edges we interpolate along the great circle connecting the upper (or lower) two corners.
1445  do i=1-nghost,im+1+nghost
1446 
1447  if (i == 1) cycle
1448 
1449  call spherical_linear_interpolation(real(i-1,kind=r_grid)/real(im,kind=R_GRID), &
1450  (/lamda(1,1),theta(1,1)/), (/lamda(im+1,1),theta(im+1,1)/), p1 )
1451  call spherical_linear_interpolation(real(i-1,kind=r_grid)/real(im,kind=R_GRID), &
1452  (/lamda(1,in+1),theta(1,in+1)/), (/lamda(im+1,in+1),theta(im+1,in+1)/), p2 )
1453 
1454  lamda(i,1) = p1(1); theta(i,1) = p1(2)
1455  lamda(i,in+1) = p2(1); theta(i,in+1) = p2(2)
1456 
1457  end do
1458 
1459  !Get cartesian coordinates and project onto the cube face with x = -rsq3
1460 
1461  i=1
1462  do j=1-nghost,in+1+nghost
1463  call latlon2xyz2(lamda(i,j), theta(i,j), pp(1,i,j))
1464  pp(2,i,j) = -pp(2,i,j)*rsq3/pp(1,i,j)
1465  pp(3,i,j) = -pp(3,i,j)*rsq3/pp(1,i,j)
1466  enddo
1467 
1468  j=1
1469  do i=1-nghost,im+1+nghost
1470  call latlon2xyz2(lamda(i,j), theta(i,j), pp(1,i,1))
1471  pp(2,i,1) = -pp(2,i,1)*rsq3/pp(1,i,1)
1472  pp(3,i,1) = -pp(3,i,1)*rsq3/pp(1,i,1)
1473  enddo
1474 
1475  !We are now on the cube.
1476 
1477  do j=1-nghost,in+1+nghost
1478  do i=1-nghost,im+1+nghost
1479  pp(1,i,j) = -rsq3
1480  enddo
1481  enddo
1482 
1483  do j=1-nghost,in+1+nghost
1484  do i=1-nghost,im+1+nghost
1485  ! Copy y-z face of the cube along j=1
1486  pp(2,i,j) = pp(2,i,1)
1487  ! Copy along i=1
1488  pp(3,i,j) = pp(3,1,j)
1489  enddo
1490  enddo
1491 
1492  call cart_to_latlon( (im+1+2*nghost)*(in+1+2*nghost), &
1493  pp(:,1-nghost:im+1+nghost,1-nghost:in+1+nghost), &
1494  lamda(1-nghost:im+1+nghost,1-nghost:in+1+nghost), &
1495  theta(1-nghost:im+1+nghost,1-nghost:in+1+nghost))
1496  !call cart_to_latlon( (im+1)*(in+1), pp(:,1:im+1,1:in+1), lamda(1:im+1,1:in+1), theta(1:im+1,1:in+1))
1497 
1498  ! Compute great-circle-distance "resolution" along the face edge:
1499  if ( is_master() ) then
1500  p1(1) = lamda(1,1); p1(2) = theta(1,1)
1501  p2(1) = lamda(2,1); p2(2) = theta(2,1)
1502  write(*,*) 'Grid x-distance at face edge (km)=',great_circle_dist( p1, p2, radius ) ! earth radius is assumed
1503  p2(1) = lamda(1,2); p2(2) = theta(1,2)
1504  write(*,*) 'Grid y-distance at face edge (km)=',great_circle_dist( p1, p2, radius ) ! earth radius is assumed
1505  !print*, 'dtheta = ', dely
1506  !print*, 'dlambda = ', lamda(2,1) - lamda(1,1)
1507  endif
1508 
1509 
1510  end subroutine gnomonic_ed_limited
1511 
1512 
1513  subroutine gnomonic_angl(im, lamda, theta)
1514 ! This is the commonly known equi-angular grid
1515  integer im
1516  real(kind=R_GRID) lamda(im+1,im+1)
1517  real(kind=R_GRID) theta(im+1,im+1)
1518  real(kind=R_GRID) p(3,im+1,im+1)
1519 ! Local
1520  real(kind=R_GRID) rsq3, xf, y0, z0, y, x, z, ds
1521  real(kind=R_GRID) dy, dz
1522  integer j,k
1523  real(kind=R_GRID) dp
1524 
1525  dp = 0.5d0*pi/real(im,kind=r_grid)
1526 
1527  rsq3 = 1.d0/sqrt(3.d0)
1528  do k=1,im+1
1529  do j=1,im+1
1530  p(1,j,k) =-rsq3 ! constant
1531  p(2,j,k) =-rsq3*tan(-0.25d0*pi+(j-1)*dp)
1532  p(3,j,k) = rsq3*tan(-0.25d0*pi+(k-1)*dp)
1533  enddo
1534  enddo
1535 
1536  call cart_to_latlon( (im+1)*(im+1), p, lamda, theta)
1537 
1538  end subroutine gnomonic_angl
1539 
1540  subroutine gnomonic_dist(im, lamda, theta)
1541 ! This is the commonly known equi-distance grid
1542  integer im
1543  real(kind=R_GRID) lamda(im+1,im+1)
1544  real(kind=R_GRID) theta(im+1,im+1)
1545  real(kind=R_GRID) p(3,im+1,im+1)
1546 ! Local
1547  real(kind=R_GRID) rsq3, xf, y0, z0, y, x, z, ds
1548  real(kind=R_GRID) dy, dz
1549  integer j,k
1550 
1551 ! Face-2
1552 
1553  rsq3 = 1.d0/sqrt(3.d0)
1554  xf = -rsq3
1555  y0 = rsq3; dy = -2.d0*rsq3/im
1556  z0 = -rsq3; dz = 2.d0*rsq3/im
1557 
1558  do k=1,im+1
1559  do j=1,im+1
1560  p(1,j,k) = xf
1561  p(2,j,k) = y0 + (j-1)*dy
1562  p(3,j,k) = z0 + (k-1)*dz
1563  enddo
1564  enddo
1565  call cart_to_latlon( (im+1)*(im+1), p, lamda, theta)
1566 
1567  end subroutine gnomonic_dist
1568 
1569  subroutine symm_ed(im, lamda, theta)
1570 ! Make grid symmetrical to i=im/2+1
1571  integer im
1572  real(kind=R_GRID) lamda(im+1,im+1)
1573  real(kind=R_GRID) theta(im+1,im+1)
1574  integer i,j,ip,jp
1575  real(kind=R_GRID) avg
1576 
1577  do j=2,im+1
1578  do i=2,im
1579  lamda(i,j) = lamda(i,1)
1580  enddo
1581  enddo
1582 
1583  do j=1,im+1
1584  do i=1,im/2
1585  ip = im + 2 - i
1586  avg = 0.5d0*(lamda(i,j)-lamda(ip,j))
1587  lamda(i, j) = avg + pi
1588  lamda(ip,j) = pi - avg
1589  avg = 0.5d0*(theta(i,j)+theta(ip,j))
1590  theta(i, j) = avg
1591  theta(ip,j) = avg
1592  enddo
1593  enddo
1594 
1595 ! Make grid symmetrical to j=im/2+1
1596  do j=1,im/2
1597  jp = im + 2 - j
1598  do i=2,im
1599  avg = 0.5d0*(lamda(i,j)+lamda(i,jp))
1600  lamda(i, j) = avg
1601  lamda(i,jp) = avg
1602  avg = 0.5d0*(theta(i,j)-theta(i,jp))
1603  theta(i, j) = avg
1604  theta(i,jp) = -avg
1605  enddo
1606  enddo
1607 
1608  end subroutine symm_ed
1609 
1610  subroutine latlon2xyz2(lon, lat, p3)
1611  real(kind=R_GRID), intent(in):: lon, lat
1612  real(kind=R_GRID), intent(out):: p3(3)
1613  real(kind=R_GRID) e(2)
1614 
1615  e(1) = lon; e(2) = lat
1616  call latlon2xyz(e, p3)
1617 
1618  end subroutine latlon2xyz2
1619 
1621  subroutine latlon2xyz(p, e, id)
1622 
1623  real(kind=R_GRID), intent(in) :: p(2)
1624  real(kind=R_GRID), intent(out):: e(3)
1625  integer, optional, intent(in):: id
1626 
1627  integer n
1628  real (f_p):: q(2)
1629  real (f_p):: e1, e2, e3
1630 
1631  do n=1,2
1632  q(n) = p(n)
1633  enddo
1634 
1635  e1 = cos(q(2)) * cos(q(1))
1636  e2 = cos(q(2)) * sin(q(1))
1637  e3 = sin(q(2))
1638 !-----------------------------------
1639 ! Truncate to the desired precision:
1640 !-----------------------------------
1641  e(1) = e1
1642  e(2) = e2
1643  e(3) = e3
1644 
1645  end subroutine latlon2xyz
1646 
1650  subroutine mirror_xyz(p1, p2, p0, p)
1651 !-------------------------------------------------------------------------------
1652 ! for k=1,2,3 (x,y,z)
1653 !
1654 ! p(k) = p0(k) - 2 * [p0(k) .dot. NB(k)] * NB(k)
1655 !
1656 ! where
1657 ! NB(k) = p1(k) .cross. p2(k) ---- direction of NB is imaterial
1658 ! the normal unit vector to the "mirror" plane
1659 !-------------------------------------------------------------------------------
1660 
1661  real(kind=R_GRID), intent(in) :: p1(3), p2(3), p0(3)
1662  real(kind=R_GRID), intent(out):: p(3)
1663 !
1664  real(kind=R_GRID):: x1, y1, z1, x2, y2, z2, x0, y0, z0
1665  real(kind=R_GRID) nb(3)
1666  real(kind=R_GRID) pdot
1667  integer k
1668 
1669  call vect_cross(nb, p1, p2)
1670  pdot = sqrt(nb(1)**2+nb(2)**2+nb(3)**2)
1671  do k=1,3
1672  nb(k) = nb(k) / pdot
1673  enddo
1674 
1675  pdot = p0(1)*nb(1) + p0(2)*nb(2) + p0(3)*nb(3)
1676  do k=1,3
1677  p(k) = p0(k) - 2.d0*pdot*nb(k)
1678  enddo
1679 
1680  end subroutine mirror_xyz
1681 
1685  subroutine mirror_latlon(lon1, lat1, lon2, lat2, lon0, lat0, lon3, lat3)
1686  real(kind=R_GRID), intent(in):: lon1, lat1, lon2, lat2, lon0, lat0
1687  real(kind=R_GRID), intent(out):: lon3, lat3
1688 !
1689  real(kind=R_GRID) p0(3), p1(3), p2(3), nb(3), pp(3), sp(2)
1690  real(kind=R_GRID) pdot
1691  integer k
1692 
1693  call latlon2xyz2(lon0, lat0, p0)
1694  call latlon2xyz2(lon1, lat1, p1)
1695  call latlon2xyz2(lon2, lat2, p2)
1696  call vect_cross(nb, p1, p2)
1697 
1698  pdot = sqrt(nb(1)**2+nb(2)**2+nb(3)**2)
1699  do k=1,3
1700  nb(k) = nb(k) / pdot
1701  enddo
1702 
1703  pdot = p0(1)*nb(1) + p0(2)*nb(2) + p0(3)*nb(3)
1704  do k=1,3
1705  pp(k) = p0(k) - 2.d0*pdot*nb(k)
1706  enddo
1707 
1708  call cart_to_latlon(1, pp, sp(1), sp(2))
1709  lon3 = sp(1)
1710  lat3 = sp(2)
1711 
1712  end subroutine mirror_latlon
1713 
1714 
1715  subroutine cart_to_latlon(np, q, xs, ys)
1716 ! vector version of cart_to_latlon1
1717  integer, intent(in):: np
1718  real(kind=R_GRID), intent(inout):: q(3,np)
1719  real(kind=R_GRID), intent(inout):: xs(np), ys(np)
1720 ! local
1721  real(kind=R_GRID), parameter:: esl=1.d-10
1722  real (f_p):: p(3)
1723  real (f_p):: dist, lat, lon
1724  integer i,k
1725 
1726  do i=1,np
1727  do k=1,3
1728  p(k) = q(k,i)
1729  enddo
1730  dist = sqrt(p(1)**2 + p(2)**2 + p(3)**2)
1731  do k=1,3
1732  p(k) = p(k) / dist
1733  enddo
1734 
1735  if ( (abs(p(1))+abs(p(2))) < esl ) then
1736  lon = real(0.,kind=f_p)
1737  else
1738  lon = atan2( p(2), p(1) ) ! range [-pi,pi]
1739  endif
1740 
1741  if ( lon < 0.) lon = real(2.,kind=f_p)*pi + lon
1742 ! RIGHT_HAND system:
1743  lat = asin(p(3))
1744 
1745  xs(i) = lon
1746  ys(i) = lat
1747 ! q Normalized:
1748  do k=1,3
1749  q(k,i) = p(k)
1750  enddo
1751  enddo
1752 
1753  end subroutine cart_to_latlon
1754 
1757  subroutine vect_cross(e, p1, p2)
1758  real(kind=R_GRID), intent(in) :: p1(3), p2(3)
1759  real(kind=R_GRID), intent(out):: e(3)
1760 
1761  e(1) = p1(2)*p2(3) - p1(3)*p2(2)
1762  e(2) = p1(3)*p2(1) - p1(1)*p2(3)
1763  e(3) = p1(1)*p2(2) - p1(2)*p2(1)
1764 
1765  end subroutine vect_cross
1766 
1767  subroutine get_center_vect( npx, npy, pp, u1, u2, bd )
1768  type(fv_grid_bounds_type), intent(IN) :: bd
1769  integer, intent(in):: npx, npy
1770  real(kind=R_GRID), intent(in) :: pp(3,bd%isd:bd%ied+1,bd%jsd:bd%jed+1)
1771  real(kind=R_GRID), intent(out):: u1(3,bd%isd:bd%ied, bd%jsd:bd%jed)
1772  real(kind=R_GRID), intent(out):: u2(3,bd%isd:bd%ied, bd%jsd:bd%jed)
1773 ! Local:
1774  integer i,j,k
1775  real(kind=R_GRID) p1(3), p2(3), pc(3), p3(3)
1776 
1777  integer :: isd, ied, jsd, jed
1778 
1779  isd = bd%isd
1780  ied = bd%ied
1781  jsd = bd%jsd
1782  jed = bd%jed
1783 
1784  do j=jsd,jed
1785  do i=isd,ied
1786  if ( (i<1 .and. j<1 ) .or. (i>(npx-1) .and. j<1) .or. &
1787  (i>(npx-1) .and. j>(npy-1)) .or. (i<1 .and. j>(npy-1))) then
1788  u1(1:3,i,j) = 0.d0
1789  u2(1:3,i,j) = 0.d0
1790  else
1791 #ifdef OLD_VECT
1792  do k=1,3
1793  u1(k,i,j) = pp(k,i+1,j)+pp(k,i+1,j+1) - pp(k,i,j)-pp(k,i,j+1)
1794  u2(k,i,j) = pp(k,i,j+1)+pp(k,i+1,j+1) - pp(k,i,j)-pp(k,i+1,j)
1795  enddo
1796  call normalize_vect( u1(1,i,j) )
1797  call normalize_vect( u2(1,i,j) )
1798 #else
1799  call cell_center3(pp(1,i,j), pp(1,i+1,j), pp(1,i,j+1), pp(1,i+1,j+1), pc)
1800 ! e1:
1801  call mid_pt3_cart(pp(1,i,j), pp(1,i,j+1), p1)
1802  call mid_pt3_cart(pp(1,i+1,j), pp(1,i+1,j+1), p2)
1803  call vect_cross(p3, p2, p1)
1804  call vect_cross(u1(1,i,j), pc, p3)
1805  call normalize_vect( u1(1,i,j) )
1806 ! e2:
1807  call mid_pt3_cart(pp(1,i,j), pp(1,i+1,j), p1)
1808  call mid_pt3_cart(pp(1,i,j+1), pp(1,i+1,j+1), p2)
1809  call vect_cross(p3, p2, p1)
1810  call vect_cross(u2(1,i,j), pc, p3)
1811  call normalize_vect( u2(1,i,j) )
1812 #endif
1813  endif
1814  enddo
1815  enddo
1816 
1817  end subroutine get_center_vect
1818 
1819 
1820  subroutine get_unit_vect2( e1, e2, uc )
1821  real(kind=R_GRID), intent(in) :: e1(2), e2(2)
1822  real(kind=R_GRID), intent(out):: uc(3)
1823 ! Local:
1824  real(kind=R_GRID), dimension(3):: pc, p1, p2, p3
1825 
1826 ! RIGHT_HAND system:
1827  call latlon2xyz(e1, p1)
1828  call latlon2xyz(e2, p2)
1829 
1830  call mid_pt3_cart(p1, p2, pc)
1831  call vect_cross(p3, p2, p1)
1832  call vect_cross(uc, pc, p3)
1833  call normalize_vect( uc )
1834 
1835  end subroutine get_unit_vect2
1836 
1837  subroutine get_unit_vect3( p1, p2, uc )
1838  real(kind=R_GRID), intent(in) :: p1(3), p2(3)
1839  real(kind=R_GRID), intent(out):: uc(3)
1840 ! Local:
1841  real(kind=R_GRID), dimension(3):: pc, p3
1842 
1843  call mid_pt3_cart(p1, p2, pc)
1844  call vect_cross(p3, p2, p1)
1845  call vect_cross(uc, pc, p3)
1846  call normalize_vect( uc )
1847 
1848  end subroutine get_unit_vect3
1849 
1851  subroutine normalize_vect(e)
1852 
1853  real(kind=R_GRID), intent(inout):: e(3)
1854  real(f_p):: pdot
1855  integer k
1856 
1857  pdot = e(1)**2 + e(2)**2 + e(3)**2
1858  pdot = sqrt( pdot )
1859 
1860  do k=1,3
1861  e(k) = e(k) / pdot
1862  enddo
1863 
1864  end subroutine normalize_vect
1865 
1866 
1867  subroutine intp_great_circle(beta, p1, p2, x_o, y_o)
1868  real(kind=R_GRID), intent(in):: beta
1869  real(kind=R_GRID), intent(in):: p1(2), p2(2)
1870  real(kind=R_GRID), intent(out):: x_o, y_o
1871 !------------------------------------------
1872  real(kind=R_GRID):: pm(2)
1873  real(kind=R_GRID):: e1(3), e2(3), e3(3)
1874  real(kind=R_GRID):: s1, s2, s3, dd, alpha
1875 
1876  call latlon2xyz(p1, e1)
1877  call latlon2xyz(p2, e2)
1878 
1879  alpha = 1.d0 - beta
1880 
1881  s1 = alpha*e1(1) + beta*e2(1)
1882  s2 = alpha*e1(2) + beta*e2(2)
1883  s3 = alpha*e1(3) + beta*e2(3)
1884 
1885  dd = sqrt( s1**2 + s2**2 + s3**2 )
1886 
1887  e3(1) = s1 / dd
1888  e3(2) = s2 / dd
1889  e3(3) = s3 / dd
1890 
1891  call cart_to_latlon(1, e3, pm(1), pm(2))
1892 
1893  x_o = pm(1)
1894  y_o = pm(2)
1895 
1896  end subroutine intp_great_circle
1897 
1901  subroutine spherical_linear_interpolation(beta, p1, p2, pb)
1902 
1903  real(kind=R_GRID), intent(in):: beta
1904  real(kind=R_GRID), intent(in):: p1(2), p2(2)
1905  real(kind=R_GRID), intent(out):: pb(2)
1906 !------------------------------------------
1907  real(kind=R_GRID):: pm(2)
1908  real(kind=R_GRID):: e1(3), e2(3), eb(3)
1909  real(kind=R_GRID):: dd, alpha, omg
1910 
1911  if ( abs(p1(1) - p2(1)) < 1.d-8 .and. abs(p1(2) - p2(2)) < 1.d-8) then
1912  call mpp_error(warning, 'spherical_linear_interpolation was passed two colocated points.')
1913  pb = p1
1914  return
1915  end if
1916 
1917  call latlon2xyz(p1, e1)
1918  call latlon2xyz(p2, e2)
1919 
1920  dd = sqrt( e1(1)**2 + e1(2)**2 + e1(3)**2 )
1921 
1922  e1(1) = e1(1) / dd
1923  e1(2) = e1(2) / dd
1924  e1(3) = e1(3) / dd
1925 
1926  dd = sqrt( e2(1)**2 + e2(2)**2 + e2(3)**2 )
1927 
1928  e2(1) = e2(1) / dd
1929  e2(2) = e2(2) / dd
1930  e2(3) = e2(3) / dd
1931 
1932  alpha = 1.d0 - beta
1933 
1934  omg = acos( e1(1)*e2(1) + e1(2)*e2(2) + e1(3)*e2(3) )
1935 
1936  if ( abs(omg) < 1.d-5 ) then
1937  print*, 'spherical_linear_interpolation: ', omg, p1, p2
1938  call mpp_error(fatal, 'spherical_linear_interpolation: interpolation not well defined between antipodal points')
1939  end if
1940 
1941  eb(1) = sin( beta*omg )*e2(1) + sin(alpha*omg)*e1(1)
1942  eb(2) = sin( beta*omg )*e2(2) + sin(alpha*omg)*e1(2)
1943  eb(3) = sin( beta*omg )*e2(3) + sin(alpha*omg)*e1(3)
1944 
1945  eb(1) = eb(1) / sin(omg)
1946  eb(2) = eb(2) / sin(omg)
1947  eb(3) = eb(3) / sin(omg)
1948 
1949  call cart_to_latlon(1, eb, pb(1), pb(2))
1950 
1951  end subroutine spherical_linear_interpolation
1952 
1953  subroutine mid_pt_sphere(p1, p2, pm)
1954  real(kind=R_GRID) , intent(IN) :: p1(2), p2(2)
1955  real(kind=R_GRID) , intent(OUT) :: pm(2)
1956 !------------------------------------------
1957  real(kind=R_GRID) e1(3), e2(3), e3(3)
1958 
1959  call latlon2xyz(p1, e1)
1960  call latlon2xyz(p2, e2)
1961  call mid_pt3_cart(e1, e2, e3)
1962  call cart_to_latlon(1, e3, pm(1), pm(2))
1963 
1964  end subroutine mid_pt_sphere
1965 
1966 
1967 
1968  subroutine mid_pt3_cart(p1, p2, e)
1969  real(kind=R_GRID), intent(IN) :: p1(3), p2(3)
1970  real(kind=R_GRID), intent(OUT) :: e(3)
1971 !
1972  real (f_p):: q1(3), q2(3)
1973  real (f_p):: dd, e1, e2, e3
1974  integer k
1975 
1976  do k=1,3
1977  q1(k) = p1(k)
1978  q2(k) = p2(k)
1979  enddo
1980 
1981  e1 = q1(1) + q2(1)
1982  e2 = q1(2) + q2(2)
1983  e3 = q1(3) + q2(3)
1984 
1985  dd = sqrt( e1**2 + e2**2 + e3**2 )
1986  e1 = e1 / dd
1987  e2 = e2 / dd
1988  e3 = e3 / dd
1989 
1990  e(1) = e1
1991  e(2) = e2
1992  e(3) = e3
1993 
1994  end subroutine mid_pt3_cart
1995 
1996 
1997 
1998  subroutine mid_pt_cart(p1, p2, e3)
1999  real(kind=R_GRID), intent(IN) :: p1(2), p2(2)
2000  real(kind=R_GRID), intent(OUT) :: e3(3)
2001 !-------------------------------------
2002  real(kind=R_GRID) e1(3), e2(3)
2003 
2004  call latlon2xyz(p1, e1)
2005  call latlon2xyz(p2, e2)
2006  call mid_pt3_cart(e1, e2, e3)
2007 
2008  end subroutine mid_pt_cart
2009 
2010 
2011 
2012  real function great_circle_dist( q1, q2, radius )
2013  real(kind=R_GRID), intent(IN) :: q1(2), q2(2)
2014  real(kind=R_GRID), intent(IN), optional :: radius
2015 
2016  real (f_p):: p1(2), p2(2)
2017  real (f_p):: beta
2018  integer n
2019 
2020  do n=1,2
2021  p1(n) = q1(n)
2022  p2(n) = q2(n)
2023  enddo
2024 
2025  beta = asin( sqrt( sin((p1(2)-p2(2))/2.)**2 + cos(p1(2))*cos(p2(2))* &
2026  sin((p1(1)-p2(1))/2.)**2 ) ) * 2.
2027 
2028  if ( present(radius) ) then
2029  great_circle_dist = radius * beta
2030  else
2031  great_circle_dist = beta ! Returns the angle
2032  endif
2033 
2034  end function great_circle_dist
2035 
2040  function great_circle_dist_cart(v1, v2, radius)
2041 
2042  real(kind=R_GRID) :: great_circle_dist_cart
2043  real(kind=R_GRID), dimension(3), intent(in) :: v1, v2
2044  real(kind=R_GRID), intent(IN), optional :: radius
2045  real(kind=R_GRID) :: norm
2046 
2047  norm = (v1(1)*v1(1)+v1(2)*v1(2)+v1(3)*v1(3)) &
2048  *(v2(1)*v2(1)+v2(2)*v2(2)+v2(3)*v2(3))
2049 
2050  !if (norm <= 0.) print*, 'negative norm: ', norm, v1, v2
2051 
2052  great_circle_dist_cart=(v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)) &
2053  /sqrt(norm)
2054  great_circle_dist_cart = sign(min(1.,abs(great_circle_dist_cart)),great_circle_dist_cart)
2055  great_circle_dist_cart=acos(great_circle_dist_cart)
2056 
2057  if ( present(radius) ) then
2058  great_circle_dist_cart = radius * great_circle_dist_cart
2059  endif
2060 
2061 
2062  end function great_circle_dist_cart
2063 
2075  subroutine intersect(a1,a2,b1,b2,radius,x_inter,local_a,local_b)
2076 
2077  real(kind=R_GRID), dimension(3), intent(in) :: a1, a2, b1, b2
2078  real(kind=R_GRID), intent(in) :: radius
2079  real(kind=R_GRID), dimension(3), intent(out) :: x_inter
2080  logical, intent(out) :: local_a,local_b
2081  !------------------------------------------------------------------!
2082  ! local variables !
2083  !------------------------------------------------------------------!
2084  real(kind=R_GRID) :: a2_xy, b1_xy, b2_xy, a2_xz, b1_xz, b2_xz, &
2085  b1_xyz, b2_xyz, length
2086  !------------------------------------------------------------------!
2087  ! calculate intersection point !
2088  !------------------------------------------------------------------!
2089  a2_xy=a2(1)*a1(2)-a2(2)*a1(1)
2090  b1_xy=b1(1)*a1(2)-b1(2)*a1(1)
2091  b2_xy=b2(1)*a1(2)-b2(2)*a1(1)
2092 
2093  a2_xz=a2(1)*a1(3)-a2(3)*a1(1)
2094  b1_xz=b1(1)*a1(3)-b1(3)*a1(1)
2095  b2_xz=b2(1)*a1(3)-b2(3)*a1(1)
2096 
2097  b1_xyz=b1_xy*a2_xz-b1_xz*a2_xy
2098  b2_xyz=b2_xy*a2_xz-b2_xz*a2_xy
2099 
2100  if (b1_xyz==0.0d0) then
2101  x_inter(:)=b1(:)
2102  elseif (b2_xyz==0.0d0) then
2103  x_inter(:)=b2(:)
2104  else
2105  x_inter(:)=b2(:)-b1(:)*b2_xyz/b1_xyz
2106  length=sqrt(x_inter(1)*x_inter(1)+x_inter(2)*x_inter(2)+x_inter(3)*x_inter(3))
2107  x_inter(:)=radius/length*x_inter(:)
2108  endif
2109  !------------------------------------------------------------------!
2110  ! check if intersection is between pairs of points on sphere !
2111  !------------------------------------------------------------------!
2112  call get_nearest()
2113  call check_local(a1,a2,local_a)
2114  call check_local(b1,b2,local_b)
2115 
2116  contains
2117  !------------------------------------------------------------------!
2118  subroutine get_nearest()
2119  real(kind=R_GRID), dimension(3) :: center, dx
2120  real(kind=R_GRID) :: dist1,dist2
2121 
2122  center(:)=0.25*(a1(:)+a2(:)+b1(:)+b2(:))
2123  dx(:)=+x_inter(:)-center(:)
2124  dist1=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2125  dx(:)=-x_inter(:)-center(:)
2126  dist2=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2127 
2128  if (dist2<dist1) x_inter(:)=-x_inter(:)
2129 
2130  end subroutine get_nearest
2131  !------------------------------------------------------------------!
2132  subroutine check_local(x1,x2,local)
2133  real(kind=R_GRID), dimension(3), intent(in) :: x1,x2
2134  logical, intent(out) :: local
2135 
2136  real(kind=R_GRID), dimension(3) :: dx
2137  real(kind=R_GRID) :: dist, dist1, dist2
2138 
2139  dx(:)=x1(:)-x2(:)
2140  dist=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2141 
2142  dx(:)=x1(:)-x_inter(:)
2143  dist1=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2144  dx(:)=x2(:)-x_inter(:)
2145  dist2=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2146 
2147  if (dist1<=dist .and. dist2<=dist) then
2148  local=.true.
2149  else
2150  local=.false.
2151  endif
2152 
2153  end subroutine check_local
2154  !------------------------------------------------------------------!
2155  end subroutine intersect
2156 
2167  subroutine intersect_cross(a1,a2,b1,b2,radius,x_inter,local_a,local_b)
2168 
2169  real(kind=R_GRID), dimension(3), intent(in) :: a1, a2, b1, b2
2170  real(kind=R_GRID), intent(in) :: radius
2171  real(kind=R_GRID), dimension(3), intent(out) :: x_inter
2172  logical, intent(out) :: local_a,local_b
2173  real(kind=R_GRID), dimension(3) :: v1, v2
2174 
2188  call vect_cross(v1, a1, a2)
2189  call vect_cross(v2, b1, b2)
2190 
2191  v1 = v1/sqrt(v1(1)**2 + v1(2)**2 + v1(3)**2)
2192  v2 = v2/sqrt(v2(1)**2 + v2(2)**2 + v2(3)**2)
2193  call vect_cross(x_inter, v1, v2)
2194 
2195  !Normalize
2196  x_inter = x_inter/sqrt(x_inter(1)**2 + x_inter(2)**2 + x_inter(3)**2)
2197 
2198  ! check if intersection is between pairs of points on sphere
2199  call get_nearest()
2200  call check_local(a1,a2,local_a)
2201  call check_local(b1,b2,local_b)
2202 
2203  contains
2204  subroutine get_nearest()
2205  real(kind=R_GRID), dimension(3) :: center, dx
2206  real(kind=R_GRID) :: dist1,dist2
2207 
2208  center(:)=0.25*(a1(:)+a2(:)+b1(:)+b2(:))
2209  dx(:)=+x_inter(:)-center(:)
2210  dist1=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2211  dx(:)=-x_inter(:)-center(:)
2212  dist2=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2213 
2214  if (dist2<dist1) x_inter(:)=-x_inter(:)
2215 
2216  end subroutine get_nearest
2217 
2218  subroutine check_local(x1,x2,local)
2219  real(kind=R_GRID), dimension(3), intent(in) :: x1,x2
2220  logical, intent(out) :: local
2221 
2222  real(kind=R_GRID), dimension(3) :: dx
2223  real(kind=R_GRID) :: dist, dist1, dist2
2224 
2225  dx(:)=x1(:)-x2(:)
2226  dist=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2227 
2228  dx(:)=x1(:)-x_inter(:)
2229  dist1=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2230  dx(:)=x2(:)-x_inter(:)
2231  dist2=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2232 
2233  if (dist1<=dist .and. dist2<=dist) then
2234  local=.true.
2235  else
2236  local=.false.
2237  endif
2238 
2239  end subroutine check_local
2240  !------------------------------------------------------------------!
2241  end subroutine intersect_cross
2242 
2243 
2244 
2245  subroutine unit_vect_latlon(pp, elon, elat)
2246  real(kind=R_GRID), intent(IN) :: pp(2)
2247  real(kind=R_GRID), intent(OUT) :: elon(3), elat(3)
2248 
2249  real (f_p):: lon, lat
2250  real (f_p):: sin_lon, cos_lon, sin_lat, cos_lat
2251 
2252  lon = pp(1)
2253  lat = pp(2)
2254 
2255  sin_lon = sin(lon)
2256  cos_lon = cos(lon)
2257  sin_lat = sin(lat)
2258  cos_lat = cos(lat)
2259 
2260  elon(1) = -sin_lon
2261  elon(2) = cos_lon
2262  elon(3) = 0.d0
2263 
2264  elat(1) = -sin_lat*cos_lon
2265  elat(2) = -sin_lat*sin_lon
2266  elat(3) = cos_lat
2267 
2268  end subroutine unit_vect_latlon
2269 
2270 
2271 
2272  real(kind=R_GRID) function v_prod(v1, v2)
2273  real(kind=R_GRID) v1(3), v2(3)
2274 
2275  v_prod = v1(1)*v2(1) + v1(2)*v2(2) + v1(3)*v2(3)
2276 
2277  end function v_prod
2278 
2279 
2280  subroutine init_cubed_to_latlon( gridstruct, hydrostatic, agrid, grid_type, ord, bd )
2281  type(fv_grid_bounds_type), intent(IN) :: bd
2282  logical, intent(in):: hydrostatic
2283  real(kind=R_GRID), intent(in) :: agrid(bd%isd:bd%ied,bd%jsd:bd%jed,2)
2284  integer, intent(in) :: grid_type
2285  integer, intent(in) :: ord
2286  type(fv_grid_type), intent(INOUT), target :: gridstruct
2287  integer i, j
2288 
2289  integer :: is, ie, js, je
2290 
2291  !Local pointers
2292  real, pointer, dimension(:,:) :: a11, a12, a21, a22
2293  real, pointer, dimension(:,:) :: z11, z12, z21, z22
2294  real(kind=R_GRID), pointer, dimension(:,:,:) :: vlon, vlat
2295  real(kind=R_GRID), pointer, dimension(:,:,:) :: ee1, ee2, ec1, ec2
2296 
2297  is = bd%is
2298  ie = bd%ie
2299  js = bd%js
2300  je = bd%je
2301 
2302  if ( grid_type < 4 ) then
2303 
2304  vlon => gridstruct%vlon
2305  vlat => gridstruct%vlat
2306  a11 => gridstruct%a11
2307  a12 => gridstruct%a12
2308  a21 => gridstruct%a21
2309  a22 => gridstruct%a22
2310  z11 => gridstruct%z11
2311  z12 => gridstruct%z12
2312  z21 => gridstruct%z21
2313  z22 => gridstruct%z22
2314  ee1 => gridstruct%ee1
2315  ee2 => gridstruct%ee2
2316  ec1 => gridstruct%ec1
2317  ec2 => gridstruct%ec2
2318 
2319  do j=js-2,je+2
2320  do i=is-2,ie+2
2321  call unit_vect_latlon(agrid(i,j,1:2), vlon(i,j,1:3), vlat(i,j,1:3))
2322  enddo
2323  enddo
2324 
2325  do j=js-1,je+1
2326  do i=is-1,ie+1
2327  z11(i,j) = v_prod(ec1(1:3,i,j), vlon(i,j,1:3))
2328  z12(i,j) = v_prod(ec1(1:3,i,j), vlat(i,j,1:3))
2329  z21(i,j) = v_prod(ec2(1:3,i,j), vlon(i,j,1:3))
2330  z22(i,j) = v_prod(ec2(1:3,i,j), vlat(i,j,1:3))
2331 !-------------------------------------------------------------------------
2332  a11(i,j) = 0.5d0*z22(i,j) / gridstruct%sin_sg(i,j,5)
2333  a12(i,j) = -0.5d0*z12(i,j) / gridstruct%sin_sg(i,j,5)
2334  a21(i,j) = -0.5d0*z21(i,j) / gridstruct%sin_sg(i,j,5)
2335  a22(i,j) = 0.5d0*z11(i,j) / gridstruct%sin_sg(i,j,5)
2336 ! For 3D Coriolis force
2337 ! if(.not.hydrostatic) gridstruct%w00(i,j) = 2.d0*omega*cos(agrid(i,j,2))
2338  enddo
2339  enddo
2340  endif
2341 
2342  end subroutine init_cubed_to_latlon
2343 
2344 
2345  subroutine cubed_to_latlon(u, v, ua, va, gridstruct, npx, npy, km, mode, grid_type, domain, nested, c2l_ord, bd)
2346  type(fv_grid_bounds_type), intent(IN) :: bd
2347  integer, intent(in) :: km, npx, npy, grid_type, c2l_ord
2348  integer, intent(in) :: mode ! update if present
2349  type(fv_grid_type), intent(IN) :: gridstruct
2350  real, intent(inout):: u(bd%isd:bd%ied,bd%jsd:bd%jed+1,km)
2351  real, intent(inout):: v(bd%isd:bd%ied+1,bd%jsd:bd%jed,km)
2352  real, intent(out):: ua(bd%isd:bd%ied, bd%jsd:bd%jed,km)
2353  real, intent(out):: va(bd%isd:bd%ied, bd%jsd:bd%jed,km)
2354  type(domain2d), intent(INOUT) :: domain
2355  logical, intent(IN) :: nested
2356 
2357  if ( c2l_ord == 2 ) then
2358  call c2l_ord2(u, v, ua, va, gridstruct, km, grid_type, bd, .false.)
2359  else
2360  call c2l_ord4(u, v, ua, va, gridstruct, npx, npy, km, grid_type, domain, nested, mode, bd)
2361  endif
2362 
2363  end subroutine cubed_to_latlon
2364 
2365 
2366  subroutine c2l_ord4(u, v, ua, va, gridstruct, npx, npy, km, grid_type, domain, nested, mode, bd)
2367 
2368  type(fv_grid_bounds_type), intent(IN) :: bd
2369  integer, intent(in) :: km, npx, npy, grid_type
2370  integer, intent(in):: mode ! update if present
2371  type(fv_grid_type), intent(IN), target :: gridstruct
2372  real, intent(inout):: u(bd%isd:bd%ied,bd%jsd:bd%jed+1,km)
2373  real, intent(inout):: v(bd%isd:bd%ied+1,bd%jsd:bd%jed,km)
2374  real, intent(out):: ua(bd%isd:bd%ied, bd%jsd:bd%jed,km)
2375  real, intent(out):: va(bd%isd:bd%ied, bd%jsd:bd%jed,km)
2376  type(domain2d), intent(INOUT) :: domain
2377  logical, intent(IN) :: nested
2378 ! Local
2379 ! 4-pt Lagrange interpolation
2380  real :: a1 = 0.5625
2381  real :: a2 = -0.0625
2382  real :: c1 = 1.125
2383  real :: c2 = -0.125
2384  real utmp(bd%is:bd%ie, bd%js:bd%je+1)
2385  real vtmp(bd%is:bd%ie+1,bd%js:bd%je)
2386  real wu(bd%is:bd%ie, bd%js:bd%je+1)
2387  real wv(bd%is:bd%ie+1,bd%js:bd%je)
2388  integer i, j, k
2389 
2390  integer :: is, ie, js, je
2391 
2392 
2393  is = bd%is
2394  ie = bd%ie
2395  js = bd%js
2396  je = bd%je
2397 
2398  if ( mode > 0 ) then
2399  call timing_on('COMM_TOTAL')
2400  call mpp_update_domains(u, v, domain, gridtype=dgrid_ne)
2401  call timing_off('COMM_TOTAL')
2402  endif
2403 
2404 !$OMP parallel do default(none) shared(is,ie,js,je,km,npx,npy,grid_type,nested,c2,c1, &
2405 !$OMP u,v,gridstruct,ua,va,a1,a2) &
2406 !$OMP private(utmp, vtmp, wu, wv)
2407  do k=1,km
2408  if ( grid_type < 4 ) then
2409  if (nested) then
2410  do j=max(1,js),min(npy-1,je)
2411  do i=max(1,is),min(npx-1,ie)
2412  utmp(i,j) = c2*(u(i,j-1,k)+u(i,j+2,k)) + c1*(u(i,j,k)+u(i,j+1,k))
2413  vtmp(i,j) = c2*(v(i-1,j,k)+v(i+2,j,k)) + c1*(v(i,j,k)+v(i+1,j,k))
2414  enddo
2415  enddo
2416  else
2417  do j=max(2,js),min(npy-2,je)
2418  do i=max(2,is),min(npx-2,ie)
2419  utmp(i,j) = c2*(u(i,j-1,k)+u(i,j+2,k)) + c1*(u(i,j,k)+u(i,j+1,k))
2420  vtmp(i,j) = c2*(v(i-1,j,k)+v(i+2,j,k)) + c1*(v(i,j,k)+v(i+1,j,k))
2421  enddo
2422  enddo
2423 
2424  if ( js==1 ) then
2425  do i=is,ie+1
2426  wv(i,1) = v(i,1,k)*gridstruct%dy(i,1)
2427  enddo
2428  do i=is,ie
2429  vtmp(i,1) = 2.*(wv(i,1) + wv(i+1,1)) / (gridstruct%dy(i,1)+gridstruct%dy(i+1,1))
2430  utmp(i,1) = 2.*(u(i,1,k)*gridstruct%dx(i,1) + u(i,2,k)*gridstruct%dx(i,2)) &
2431  / ( gridstruct%dx(i,1) + gridstruct%dx(i,2))
2432 !!! vtmp(i,1) = (wv(i,1) + wv(i+1,1)) * gridstruct%rdya(i,1)
2433 !!! utmp(i,1) = (u(i,1,k)*gridstruct%dx(i,1) + u(i,2,k)*gridstruct%dx(i,2)) * gridstruct%rdxa(i,1)
2434  enddo
2435  endif
2436 
2437  if ( (je+1)==npy ) then
2438  j = npy-1
2439  do i=is,ie+1
2440  wv(i,j) = v(i,j,k)*gridstruct%dy(i,j)
2441  enddo
2442  do i=is,ie
2443  vtmp(i,j) = 2.*(wv(i,j) + wv(i+1,j)) / (gridstruct%dy(i,j)+gridstruct%dy(i+1,j))
2444  utmp(i,j) = 2.*(u(i,j,k)*gridstruct%dx(i,j) + u(i,j+1,k)*gridstruct%dx(i,j+1)) &
2445  / ( gridstruct%dx(i,j) + gridstruct%dx(i,j+1))
2446 !!! vtmp(i,j) = (wv(i,j) + wv(i+1,j)) * gridstruct%rdya(i,j)
2447 !!! utmp(i,j) = (u(i,j,k)*gridstruct%dx(i,j) + u(i,j+1,k)*gridstruct%dx(i,j+1)) * gridstruct%rdxa(i,j)
2448  enddo
2449  endif
2450 
2451  if ( is==1 ) then
2452  i = 1
2453  do j=js,je
2454  wv(1,j) = v(1,j,k)*gridstruct%dy(1,j)
2455  wv(2,j) = v(2,j,k)*gridstruct%dy(2,j)
2456  enddo
2457  do j=js,je+1
2458  wu(i,j) = u(i,j,k)*gridstruct%dx(i,j)
2459  enddo
2460  do j=js,je
2461  utmp(i,j) = 2.*(wu(i,j) + wu(i,j+1))/(gridstruct%dx(i,j)+gridstruct%dx(i,j+1))
2462  vtmp(i,j) = 2.*(wv(1,j) + wv(2,j ))/(gridstruct%dy(1,j)+gridstruct%dy(2,j))
2463 !!! utmp(i,j) = (wu(i,j) + wu(i, j+1)) * gridstruct%rdxa(i,j)
2464 !!! vtmp(i,j) = (wv(i,j) + wv(i+1,j )) * gridstruct%rdya(i,j)
2465  enddo
2466  endif
2467 
2468  if ( (ie+1)==npx) then
2469  i = npx-1
2470  do j=js,je
2471  wv(i, j) = v(i, j,k)*gridstruct%dy(i, j)
2472  wv(i+1,j) = v(i+1,j,k)*gridstruct%dy(i+1,j)
2473  enddo
2474  do j=js,je+1
2475  wu(i,j) = u(i,j,k)*gridstruct%dx(i,j)
2476  enddo
2477  do j=js,je
2478  utmp(i,j) = 2.*(wu(i,j) + wu(i, j+1))/(gridstruct%dx(i,j)+gridstruct%dx(i,j+1))
2479  vtmp(i,j) = 2.*(wv(i,j) + wv(i+1,j ))/(gridstruct%dy(i,j)+gridstruct%dy(i+1,j))
2480 !!! utmp(i,j) = (wu(i,j) + wu(i, j+1)) * gridstruct%rdxa(i,j)
2481 !!! vtmp(i,j) = (wv(i,j) + wv(i+1,j )) * gridstruct%rdya(i,j)
2482  enddo
2483  endif
2484 
2485  endif !nested
2486 
2487  !Transform local a-grid winds into latitude-longitude coordinates
2488  do j=js,je
2489  do i=is,ie
2490  ua(i,j,k) = gridstruct%a11(i,j)*utmp(i,j) + gridstruct%a12(i,j)*vtmp(i,j)
2491  va(i,j,k) = gridstruct%a21(i,j)*utmp(i,j) + gridstruct%a22(i,j)*vtmp(i,j)
2492  enddo
2493  enddo
2494  else
2495 ! Simple Cartesian Geometry:
2496  do j=js,je
2497  do i=is,ie
2498  ua(i,j,k) = a2*(u(i,j-1,k)+u(i,j+2,k)) + a1*(u(i,j,k)+u(i,j+1,k))
2499  va(i,j,k) = a2*(v(i-1,j,k)+v(i+2,j,k)) + a1*(v(i,j,k)+v(i+1,j,k))
2500  enddo
2501  enddo
2502  endif
2503  enddo
2504  end subroutine c2l_ord4
2505 
2506  subroutine c2l_ord2(u, v, ua, va, gridstruct, km, grid_type, bd, do_halo)
2507  type(fv_grid_bounds_type), intent(IN) :: bd
2508  integer, intent(in) :: km, grid_type
2509  real, intent(in) :: u(bd%isd:bd%ied,bd%jsd:bd%jed+1,km)
2510  real, intent(in) :: v(bd%isd:bd%ied+1,bd%jsd:bd%jed,km)
2511  type(fv_grid_type), intent(IN), target :: gridstruct
2512  logical, intent(in) :: do_halo
2513 !
2514  real, intent(out):: ua(bd%isd:bd%ied, bd%jsd:bd%jed,km)
2515  real, intent(out):: va(bd%isd:bd%ied, bd%jsd:bd%jed,km)
2516 !--------------------------------------------------------------
2517 ! Local
2518  real wu(bd%is-1:bd%ie+1, bd%js-1:bd%je+2)
2519  real wv(bd%is-1:bd%ie+2, bd%js-1:bd%je+1)
2520  real u1(bd%is-1:bd%ie+1), v1(bd%is-1:bd%ie+1)
2521  integer i, j, k
2522  integer :: is, ie, js, je
2523 
2524  real, dimension(:,:), pointer :: a11, a12, a21, a22
2525  real, dimension(:,:), pointer :: dx, dy, rdxa, rdya
2526 
2527  a11 => gridstruct%a11
2528  a12 => gridstruct%a12
2529  a21 => gridstruct%a21
2530  a22 => gridstruct%a22
2531 
2532  dx => gridstruct%dx
2533  dy => gridstruct%dy
2534  rdxa => gridstruct%rdxa
2535  rdya => gridstruct%rdya
2536 
2537  if (do_halo) then
2538  is = bd%is-1
2539  ie = bd%ie+1
2540  js = bd%js-1
2541  je = bd%je+1
2542  else
2543  is = bd%is
2544  ie = bd%ie
2545  js = bd%js
2546  je = bd%je
2547  endif
2548 
2549 !$OMP parallel do default(none) shared(is,ie,js,je,km,grid_type,u,dx,v,dy,ua,va,a11,a12,a21,a22) &
2550 !$OMP private(u1, v1, wu, wv)
2551  do k=1,km
2552  if ( grid_type < 4 ) then
2553  do j=js,je+1
2554  do i=is,ie
2555  wu(i,j) = u(i,j,k)*dx(i,j)
2556  enddo
2557  enddo
2558  do j=js,je
2559  do i=is,ie+1
2560  wv(i,j) = v(i,j,k)*dy(i,j)
2561  enddo
2562  enddo
2563 
2564  do j=js,je
2565  do i=is,ie
2566 ! Co-variant to Co-variant "vorticity-conserving" interpolation
2567  u1(i) = 2.*(wu(i,j) + wu(i,j+1)) / (dx(i,j)+dx(i,j+1))
2568  v1(i) = 2.*(wv(i,j) + wv(i+1,j)) / (dy(i,j)+dy(i+1,j))
2569 !!! u1(i) = (wu(i,j) + wu(i,j+1)) * rdxa(i,j)
2570 !!! v1(i) = (wv(i,j) + wv(i+1,j)) * rdya(i,j)
2571 ! Cubed (cell center co-variant winds) to lat-lon:
2572  ua(i,j,k) = a11(i,j)*u1(i) + a12(i,j)*v1(i)
2573  va(i,j,k) = a21(i,j)*u1(i) + a22(i,j)*v1(i)
2574  enddo
2575  enddo
2576  else
2577 ! 2nd order:
2578  do j=js,je
2579  do i=is,ie
2580  ua(i,j,k) = 0.5*(u(i,j,k)+u(i, j+1,k))
2581  va(i,j,k) = 0.5*(v(i,j,k)+v(i+1,j, k))
2582  enddo
2583  enddo
2584  endif
2585  enddo
2586 
2587  end subroutine c2l_ord2
2588 
2589 
2590  subroutine expand_cell(q1, q2, q3, q4, a1, a2, a3, a4, fac)
2591 ! Util for land model (for BW)
2592 !
2593 ! 4----3
2594 ! | . |
2595 ! 1----2
2596 !
2597  real(kind=R_GRID), intent(in):: q1(2), q2(2), q3(2), q4(2)
2598  real(kind=R_GRID), intent(in):: fac ! expansion factor: outside: > 1
2599  ! fac = 1: qq1 returns q1
2600  ! fac = 0: qq1 returns the center position
2601  real(kind=R_GRID), intent(out):: a1(2), a2(2), a3(2), a4(2)
2602 ! Local
2603  real(kind=R_GRID) qq1(3), qq2(3), qq3(3), qq4(3)
2604  real(kind=R_GRID) p1(3), p2(3), p3(3), p4(3)
2605  real(kind=R_GRID) ec(3)
2606  real(f_p):: dd, d1, d2, d3, d4
2607  integer k
2608 
2609 ! Transform to (x,y,z)
2610  call latlon2xyz(q1, p1)
2611  call latlon2xyz(q2, p2)
2612  call latlon2xyz(q3, p3)
2613  call latlon2xyz(q4, p4)
2614 
2615 ! Get center position:
2616  do k=1,3
2617  ec(k) = p1(k) + p2(k) + p3(k) + p4(k)
2618  enddo
2619  dd = sqrt( ec(1)**2 + ec(2)**2 + ec(3)**2 )
2620 
2621  do k=1,3
2622  ec(k) = ec(k) / dd ! cell center position
2623  enddo
2624 
2625 ! Perform the "extrapolation" in 3D (x-y-z)
2626  do k=1,3
2627  qq1(k) = ec(k) + fac*(p1(k)-ec(k))
2628  qq2(k) = ec(k) + fac*(p2(k)-ec(k))
2629  qq3(k) = ec(k) + fac*(p3(k)-ec(k))
2630  qq4(k) = ec(k) + fac*(p4(k)-ec(k))
2631  enddo
2632 
2633 !--------------------------------------------------------
2634 ! Force the points to be on the sphere via normalization
2635 !--------------------------------------------------------
2636  d1 = sqrt( qq1(1)**2 + qq1(2)**2 + qq1(3)**2 )
2637  d2 = sqrt( qq2(1)**2 + qq2(2)**2 + qq2(3)**2 )
2638  d3 = sqrt( qq3(1)**2 + qq3(2)**2 + qq3(3)**2 )
2639  d4 = sqrt( qq4(1)**2 + qq4(2)**2 + qq4(3)**2 )
2640  do k=1,3
2641  qq1(k) = qq1(k) / d1
2642  qq2(k) = qq2(k) / d2
2643  qq3(k) = qq3(k) / d3
2644  qq4(k) = qq4(k) / d4
2645  enddo
2646 
2647 !----------------------------------------
2648 ! Transform back to lat-lon coordinates:
2649 !----------------------------------------
2650 
2651  call cart_to_latlon(1, qq1, a1(1), a1(2))
2652  call cart_to_latlon(1, qq2, a2(1), a2(2))
2653  call cart_to_latlon(1, qq3, a3(1), a3(2))
2654  call cart_to_latlon(1, qq4, a4(1), a4(2))
2655 
2656  end subroutine expand_cell
2657 
2658 
2659  subroutine cell_center2(q1, q2, q3, q4, e2)
2660  real(kind=R_GRID) , intent(in ) :: q1(2), q2(2), q3(2), q4(2)
2661  real(kind=R_GRID) , intent(out) :: e2(2)
2662 ! Local
2663  real(kind=R_GRID) p1(3), p2(3), p3(3), p4(3)
2664  real(kind=R_GRID) ec(3)
2665  real(kind=R_GRID) dd
2666  integer k
2667 
2668  call latlon2xyz(q1, p1)
2669  call latlon2xyz(q2, p2)
2670  call latlon2xyz(q3, p3)
2671  call latlon2xyz(q4, p4)
2672 
2673  do k=1,3
2674  ec(k) = p1(k) + p2(k) + p3(k) + p4(k)
2675  enddo
2676  dd = sqrt( ec(1)**2 + ec(2)**2 + ec(3)**2 )
2677 
2678  do k=1,3
2679  ec(k) = ec(k) / dd
2680  enddo
2681 
2682  call cart_to_latlon(1, ec, e2(1), e2(2))
2683 
2684  end subroutine cell_center2
2685 
2687  subroutine cell_center3(p1, p2, p3, p4, ec)
2688  real(kind=R_GRID) , intent(IN) :: p1(3), p2(3), p3(3), p4(3)
2689  real(kind=R_GRID) , intent(OUT) :: ec(3)
2690 ! Local
2691  real(kind=r_grid)dd
2692  integer k
2693 
2694  do k=1,3
2695  ec(k) = p1(k) + p2(k) + p3(k) + p4(k)
2696  enddo
2697  dd = sqrt( ec(1)**2 + ec(2)**2 + ec(3)**2 )
2698 
2699  do k=1,3
2700  ec(k) = ec(k) / dd
2701  enddo
2702 
2703  end subroutine cell_center3
2704 
2705 
2706 
2707  real(kind=R_GRID) function get_area(p1, p4, p2, p3, radius)
2708 !-----------------------------------------------
2709  real(kind=R_GRID), intent(in), dimension(2):: p1, p2, p3, p4
2710  real(kind=R_GRID), intent(in), optional:: radius
2711 !-----------------------------------------------
2712  real(kind=R_GRID) e1(3), e2(3), e3(3)
2713  real(kind=R_GRID) ang1, ang2, ang3, ang4
2714 
2715 ! S-W: 1
2716  call latlon2xyz(p1, e1) ! p1
2717  call latlon2xyz(p2, e2) ! p2
2718  call latlon2xyz(p4, e3) ! p4
2719  ang1 = spherical_angle(e1, e2, e3)
2720 !----
2721 ! S-E: 2
2722 !----
2723  call latlon2xyz(p2, e1)
2724  call latlon2xyz(p3, e2)
2725  call latlon2xyz(p1, e3)
2726  ang2 = spherical_angle(e1, e2, e3)
2727 !----
2728 ! N-E: 3
2729 !----
2730  call latlon2xyz(p3, e1)
2731  call latlon2xyz(p4, e2)
2732  call latlon2xyz(p2, e3)
2733  ang3 = spherical_angle(e1, e2, e3)
2734 !----
2735 ! N-W: 4
2736 !----
2737  call latlon2xyz(p4, e1)
2738  call latlon2xyz(p3, e2)
2739  call latlon2xyz(p1, e3)
2740  ang4 = spherical_angle(e1, e2, e3)
2741 
2742  if ( present(radius) ) then
2743  get_area = (ang1 + ang2 + ang3 + ang4 - 2.*pi) * radius**2
2744  else
2745  get_area = ang1 + ang2 + ang3 + ang4 - 2.*pi
2746  endif
2747 
2748  end function get_area
2749 
2755  function dist2side(v1, v2, point)
2756  real(kind=R_GRID) :: dist2side
2757  real(kind=R_GRID), dimension(3), intent(in) :: v1, v2, point
2758 
2759  real(kind=R_GRID) :: angle, side
2760 
2761  angle = spherical_angle(v1, v2, point)
2762  side = great_circle_dist_cart(v1, point)
2763  dist2side = asin(sin(side)*sin(angle))
2764 
2765  end function dist2side
2766 
2769  function dist2side_latlon(v1,v2,point)
2770 
2771  real(kind=R_GRID) :: dist2side_latlon
2772  real(kind=R_GRID), dimension(2), intent(in) :: v1, v2, point
2773 
2774  real(kind=R_GRID),dimension(3) :: c1, c2, cpoint
2775 
2776  real(kind=R_GRID) :: angle,side
2777 
2778  !no version of spherical angle for lat-lon coords
2779  call latlon2xyz(v1,c1)
2780  call latlon2xyz(v2,c2)
2781  call latlon2xyz(point,cpoint)
2782  angle = spherical_angle(c1,c2,cpoint)
2783 
2784  side = great_circle_dist(v1,point)
2785 
2786  dist2side_latlon = asin(sin(side)*sin(angle))
2787 
2788  !!dist2side_latlon = dist2side(c1,c2,cpoint)
2789 
2790  end function dist2side_latlon
2791 
2792 
2793 
2794  real(kind=R_GRID) function spherical_angle(p1, p2, p3)
2795 
2796 ! p3
2797 ! /
2798 ! /
2799 ! p1 ---> angle
2800 ! \
2801 ! \
2802 ! p2
2803 
2804  real(kind=R_GRID) p1(3), p2(3), p3(3)
2805 
2806  real (f_p):: e1(3), e2(3), e3(3)
2807  real (f_p):: px, py, pz
2808  real (f_p):: qx, qy, qz
2809  real (f_p):: angle, ddd
2810  integer n
2811 
2812  do n=1,3
2813  e1(n) = p1(n)
2814  e2(n) = p2(n)
2815  e3(n) = p3(n)
2816  enddo
2817 
2818 !-------------------------------------------------------------------
2819 ! Page 41, Silverman's book on Vector Algebra; spherical trigonmetry
2820 !-------------------------------------------------------------------
2821 ! Vector P:
2822  px = e1(2)*e2(3) - e1(3)*e2(2)
2823  py = e1(3)*e2(1) - e1(1)*e2(3)
2824  pz = e1(1)*e2(2) - e1(2)*e2(1)
2825 ! Vector Q:
2826  qx = e1(2)*e3(3) - e1(3)*e3(2)
2827  qy = e1(3)*e3(1) - e1(1)*e3(3)
2828  qz = e1(1)*e3(2) - e1(2)*e3(1)
2829 
2830  ddd = (px*px+py*py+pz*pz)*(qx*qx+qy*qy+qz*qz)
2831 
2832  if ( ddd <= 0.0d0 ) then
2833  angle = 0.d0
2834  else
2835  ddd = (px*qx+py*qy+pz*qz) / sqrt(ddd)
2836  if ( abs(ddd)>1.d0) then
2837  angle = 2.d0*atan(1.0) ! 0.5*pi
2838  !FIX (lmh) to correctly handle co-linear points (angle near pi or 0)
2839  if (ddd < 0.d0) then
2840  angle = 4.d0*atan(1.0d0) !should be pi
2841  else
2842  angle = 0.d0
2843  end if
2844  else
2845  angle = acos( ddd )
2846  endif
2847  endif
2848 
2849  spherical_angle = angle
2850 
2851  end function spherical_angle
2852 
2853 
2854  real(kind=R_GRID) function cos_angle(p1, p2, p3)
2855 ! As spherical_angle, but returns the cos(angle)
2856 ! p3
2857 ! ^
2858 ! |
2859 ! |
2860 ! p1 ---> p2
2861 !
2862  real(kind=R_GRID), intent(in):: p1(3), p2(3), p3(3)
2863 
2864  real (f_p):: e1(3), e2(3), e3(3)
2865  real (f_p):: px, py, pz
2866  real (f_p):: qx, qy, qz
2867  real (f_p):: angle, ddd
2868  integer n
2869 
2870  do n=1,3
2871  e1(n) = p1(n)
2872  e2(n) = p2(n)
2873  e3(n) = p3(n)
2874  enddo
2875 
2876 !-------------------------------------------------------------------
2877 ! Page 41, Silverman's book on Vector Algebra; spherical trigonmetry
2878 !-------------------------------------------------------------------
2879 ! Vector P:= e1 X e2
2880  px = e1(2)*e2(3) - e1(3)*e2(2)
2881  py = e1(3)*e2(1) - e1(1)*e2(3)
2882  pz = e1(1)*e2(2) - e1(2)*e2(1)
2883 
2884 ! Vector Q: e1 X e3
2885  qx = e1(2)*e3(3) - e1(3)*e3(2)
2886  qy = e1(3)*e3(1) - e1(1)*e3(3)
2887  qz = e1(1)*e3(2) - e1(2)*e3(1)
2888 
2889 ! ddd = sqrt[ (P*P) (Q*Q) ]
2890  ddd = sqrt( (px**2+py**2+pz**2)*(qx**2+qy**2+qz**2) )
2891  if ( ddd > 0.d0 ) then
2892  angle = (px*qx+py*qy+pz*qz) / ddd
2893  else
2894  angle = 1.d0
2895  endif
2896  cos_angle = angle
2897 
2898  end function cos_angle
2899 
2901  real function g_sum(domain, p, ifirst, ilast, jfirst, jlast, ngc, area, mode, reproduce)
2902  integer, intent(IN) :: ifirst, ilast
2903  integer, intent(IN) :: jfirst, jlast, ngc
2904  integer, intent(IN) :: mode ! if ==1 divided by area
2905  logical, intent(in), optional :: reproduce
2906  real, intent(IN) :: p(ifirst:ilast,jfirst:jlast) ! field to be summed
2907  real(kind=R_GRID), intent(IN) :: area(ifirst-ngc:ilast+ngc,jfirst-ngc:jlast+ngc)
2908  type(domain2d), intent(IN) :: domain
2909  integer :: i,j
2910  real gsum
2911  logical, SAVE :: g_sum_initialized = .false.
2912  real(kind=R_GRID), SAVE :: global_area
2913  real :: tmp(ifirst:ilast,jfirst:jlast)
2914 
2915  if ( .not. g_sum_initialized ) then
2916  global_area = mpp_global_sum(domain, area, flags=bitwise_efp_sum)
2917  if ( is_master() ) write(*,*) 'Global Area=',global_area
2918  g_sum_initialized = .true.
2919  end if
2920 
2921 !-------------------------
2922 ! FMS global sum algorithm:
2923 !-------------------------
2924  if ( present(reproduce) ) then
2925  if (reproduce) then
2926  gsum = mpp_global_sum(domain, p(:,:)*area(ifirst:ilast,jfirst:jlast), &
2927  flags=bitwise_efp_sum)
2928  else
2929  gsum = mpp_global_sum(domain, p(:,:)*area(ifirst:ilast,jfirst:jlast))
2930  endif
2931  else
2932 !-------------------------
2933 ! Quick local sum algorithm
2934 !-------------------------
2935  gsum = 0.
2936  do j=jfirst,jlast
2937  do i=ifirst,ilast
2938  gsum = gsum + p(i,j)*area(i,j)
2939  enddo
2940  enddo
2941  call mp_reduce_sum(gsum)
2942  endif
2943 
2944  if ( mode==1 ) then
2945  g_sum = gsum / global_area
2946  else
2947  g_sum = gsum
2948  endif
2949 
2950  end function g_sum
2951 
2953  real function global_qsum(p, ifirst, ilast, jfirst, jlast)
2954  integer, intent(IN) :: ifirst, ilast
2955  integer, intent(IN) :: jfirst, jlast
2956  real, intent(IN) :: p(ifirst:ilast,jfirst:jlast)
2957  integer :: i,j
2958  real gsum
2959 
2960  gsum = 0.
2961  do j=jfirst,jlast
2962  do i=ifirst,ilast
2963  gsum = gsum + p(i,j)
2964  enddo
2965  enddo
2966  call mp_reduce_sum(gsum)
2967 
2968  global_qsum = gsum
2969 
2970  end function global_qsum
2971 
2972 
2973  subroutine global_mx(q, n_g, qmin, qmax, bd)
2974 
2975  type(fv_grid_bounds_type), intent(IN) :: bd
2976  integer, intent(in):: n_g
2977  real(kind=R_GRID), intent(in)::q(bd%is-n_g:bd%ie+n_g, bd%js-n_g:bd%je+n_g)
2978  real(kind=R_GRID), intent(out):: qmin, qmax
2979  integer i,j
2980 
2981  integer :: is, ie, js, je
2982 
2983  is = bd%is
2984  ie = bd%ie
2985  js = bd%js
2986  je = bd%je
2987 
2988  qmin = q(is,js)
2989  qmax = qmin
2990  do j=js,je
2991  do i=is,ie
2992  qmin = min(qmin, q(i,j))
2993  qmax = max(qmax, q(i,j))
2994  enddo
2995  enddo
2996  call mp_reduce_min(qmin)
2997  call mp_reduce_max(qmax)
2998 
2999  end subroutine global_mx
3000 
3002  subroutine global_mx_c(q, i1, i2, j1, j2, qmin, qmax)
3003  integer, intent(in):: i1, i2, j1, j2
3004  real(kind=R_GRID), intent(in) :: q(i1:i2,j1:j2)
3005  real(kind=R_GRID), intent(out) :: qmin, qmax
3006  integer i,j
3007 
3008  qmin = q(i1,j1)
3009  qmax = qmin
3010  do j=j1,j2
3011  do i=i1,i2
3012  qmin = min(qmin, q(i,j))
3013  qmax = max(qmax, q(i,j))
3014  enddo
3015  enddo
3016  call mp_reduce_min(qmin)
3017  call mp_reduce_max(qmax)
3018 
3019  end subroutine global_mx_c
3020 
3021 
3022 #ifdef OVERLOAD_R4
3023  subroutine fill_ghost_r4(q, npx, npy, value, bd)
3024 
3025  type(fv_grid_bounds_type), intent(IN) :: bd
3026  real(kind=4), intent(inout):: q(bd%isd:bd%ied,bd%jsd:bd%jed)
3027  integer, intent(in):: npx, npy
3028  real, intent(in):: value
3029  integer i,j
3030 
3031  integer :: is, ie, js, je
3032  integer :: isd, ied, jsd, jed
3033 
3034  is = bd%is
3035  ie = bd%ie
3036  js = bd%js
3037  je = bd%je
3038  isd = bd%isd
3039  ied = bd%ied
3040  jsd = bd%jsd
3041  jed = bd%jed
3042 
3043  do j=jsd,jed
3044  do i=isd,ied
3045  if ( (i<1 .and. j<1) ) then
3046  q(i,j) = value
3047  endif
3048  if ( i>(npx-1) .and. j<1 ) then
3049  q(i,j) = value
3050  endif
3051  if ( i>(npx-1) .and. j>(npy-1) ) then
3052  q(i,j) = value
3053  endif
3054  if ( i<1 .and. j>(npy-1) ) then
3055  q(i,j) = value
3056  endif
3057  enddo
3058  enddo
3059 
3060  end subroutine fill_ghost_r4
3061 #endif
3062 
3063  subroutine fill_ghost_r8(q, npx, npy, value, bd)
3064 
3065  type(fv_grid_bounds_type), intent(IN) :: bd
3066  real(kind=R_GRID), intent(inout):: q(bd%isd:bd%ied,bd%jsd:bd%jed)
3067  integer, intent(in):: npx, npy
3068  real, intent(in):: value
3069  integer i,j
3070 
3071  integer :: is, ie, js, je
3072  integer :: isd, ied, jsd, jed
3073 
3074  is = bd%is
3075  ie = bd%ie
3076  js = bd%js
3077  je = bd%je
3078  isd = bd%isd
3079  ied = bd%ied
3080  jsd = bd%jsd
3081  jed = bd%jed
3082 
3083  do j=jsd,jed
3084  do i=isd,ied
3085  if ( (i<1 .and. j<1) ) then
3086  q(i,j) = value
3087  endif
3088  if ( i>(npx-1) .and. j<1 ) then
3089  q(i,j) = value
3090  endif
3091  if ( i>(npx-1) .and. j>(npy-1) ) then
3092  q(i,j) = value
3093  endif
3094  if ( i<1 .and. j>(npy-1) ) then
3095  q(i,j) = value
3096  endif
3097  enddo
3098  enddo
3099 
3100  end subroutine fill_ghost_r8
3101 
3102 
3103 
3104  subroutine make_eta_level(km, pe, area, kks, ak, bk, ptop, domain, bd)
3105  type(fv_grid_bounds_type), intent(IN) :: bd
3106  integer, intent(in ):: km
3107  integer, intent(out):: kks
3108  real(kind=R_GRID), intent(in):: area(bd%isd:bd%ied,bd%jsd:bd%jed)
3109  real, intent(INOUT) :: ptop
3110  real, intent(inout):: pe(bd%is-1:bd%ie+1,km+1,bd%js-1:bd%je+1)
3111  real, intent(out):: ak(km+1), bk(km+1)
3112  type(domain2d), intent(IN) :: domain
3113 ! local:
3114  real ph(km+1)
3115  real, allocatable:: pem(:,:)
3116  real(kind=4) :: p4
3117  integer k, i, j
3118  integer :: is, ie, js, je
3119 
3120  is = bd%is
3121  ie = bd%ie
3122  js = bd%js
3123  je = bd%je
3124 
3125  allocate ( pem(is:ie,js:je) )
3126 
3127 ! Compute global mean values:
3128  do k=1,km+1
3129  do j=js,je
3130  do i=is,ie
3131  pem(i,j) = pe(i,k,j)
3132  enddo
3133  enddo
3134 ! Make it the same across all PEs
3135 ! ph(k) = g_sum(pem, is, ie, js, je, ng, area, 1, .true.)
3136  p4 = g_sum(domain, pem, is, ie, js, je, ng, area, 1)
3137  ph(k) = p4
3138  enddo
3139 
3140  ptop = ph(1)
3141  do j=js-1,je+1
3142  do i=is-1,ie+1
3143  pe(i,1,j) = ptop
3144  enddo
3145  enddo
3146 
3147 ! Faking a's and b's for code compatibility with hybrid sigma-p
3148  kks = 0
3149  ak(1) = ph(1)
3150  bk(1) = 0.
3151  ak(km+1) = 0.
3152  bk(km+1) = 1.
3153 
3154  do k=2,km
3155  bk(k) = (ph(k) - ph(1)) / (ph(km+1)-ph(1))
3156  ak(k) = ph(1)*(1.-bk(k))
3157  enddo
3158 
3159  if ( is_master() ) then
3160  write(*,*) 'Make_eta_level ...., ptop=', ptop
3161 #ifdef PRINT_GRID
3162  do k=1,km+1
3163  write(*,*) ph(k), ak(k), bk(k)
3164  enddo
3165 #endif
3166  endif
3167 
3168  deallocate ( pem )
3169 
3170  end subroutine make_eta_level
3171 
3172  subroutine invert_matrix(n, a, x)
3173  integer, intent (in) :: n
3174  integer :: i,j,k
3175  real(kind=R_GRID), intent (inout), dimension (n,n):: a
3176  real(kind=R_GRID), intent (out), dimension (n,n):: x
3177  real(kind=R_GRID), dimension (n,n) :: b
3178  integer indx(n)
3179 
3180  do i = 1, n
3181  do j = 1, n
3182  b(i,j) = 0.0
3183  end do
3184  end do
3185 
3186  do i = 1, n
3187  b(i,i) = 1.0
3188  end do
3189 
3190  call elgs (a,n,indx)
3191 
3192  do i = 1, n-1
3193  do j = i+1, n
3194  do k = 1, n
3195  b(indx(j),k) = b(indx(j),k) - a(indx(j),i)*b(indx(i),k)
3196  end do
3197  end do
3198  end do
3199 
3200  do i = 1, n
3201  x(n,i) = b(indx(n),i)/a(indx(n),n)
3202  do j = n-1, 1, -1
3203  x(j,i) = b(indx(j),i)
3204  do k = j+1, n
3205  x(j,i) = x(j,i)-a(indx(j),k)*x(k,i)
3206  end do
3207  x(j,i) = x(j,i)/a(indx(j),j)
3208  end do
3209  end do
3210 
3211  end subroutine invert_matrix
3212 
3216  subroutine elgs (a,n,indx)
3217  integer, intent (in) :: n
3218  integer :: i,j,k,itmp
3219  integer, intent (out), dimension (n) :: indx
3220  real(kind=R_GRID), intent (inout), dimension (n,n) :: a
3221 !
3222  real(kind=R_GRID) :: c1, pie, pi1, pj
3223  real(kind=R_GRID), dimension (n) :: c
3224 
3225  do i = 1, n
3226  indx(i) = i
3227  end do
3228 !
3229 ! find the rescaling factors, one from each row
3230 !
3231  do i = 1, n
3232  c1= 0.0
3233  do j = 1, n
3234  c1 = max(c1,abs(a(i,j)))
3235  end do
3236  c(i) = c1
3237  end do
3238 !
3239 ! search the pivoting (largest) element from each column
3240 !
3241  do j = 1, n-1
3242  pi1 = 0.0
3243  do i = j, n
3244  pie = abs(a(indx(i),j))/c(indx(i))
3245  if (pie > pi1) then
3246  pi1 = pie
3247  k = i
3248  endif
3249  end do
3250 !
3251 ! interchange the rows via indx(n) to record pivoting order
3252 !
3253  itmp = indx(j)
3254  indx(j) = indx(k)
3255  indx(k) = itmp
3256  do i = j+1, n
3257  pj = a(indx(i),j)/a(indx(j),j)
3258 !
3259 ! record pivoting ratios below the diagonal
3260 !
3261  a(indx(i),j) = pj
3262 !
3263 ! modify other elements accordingly
3264 !
3265  do k = j+1, n
3266  a(indx(i),k) = a(indx(i),k)-pj*a(indx(j),k)
3267  end do
3268  end do
3269  end do
3270 
3271  end subroutine elgs
3272 
3273  subroutine get_latlon_vector(pp, elon, elat)
3274  real(kind=R_GRID), intent(IN) :: pp(2)
3275  real(kind=R_GRID), intent(OUT) :: elon(3), elat(3)
3276 
3277  elon(1) = -sin(pp(1))
3278  elon(2) = cos(pp(1))
3279  elon(3) = 0.0
3280  elat(1) = -sin(pp(2))*cos(pp(1))
3281  elat(2) = -sin(pp(2))*sin(pp(1))
3282 !!! RIGHT_HAND
3283  elat(3) = cos(pp(2))
3284 ! Left-hand system needed to be consistent with rest of the codes
3285 ! elat(3) = -COS(pp(2))
3286 
3287  end subroutine get_latlon_vector
3288 
3289 
3290 
3291 
3292  subroutine project_sphere_v( np, f, e )
3293 !---------------------------------
3294  integer, intent(in):: np
3295  real(kind=R_GRID), intent(in):: e(3,np)
3296  real(kind=R_GRID), intent(inout):: f(3,np)
3297 ! local
3298  real(f_p):: ap
3299  integer i
3300 
3301  do i=1,np
3302  ap = f(1,i)*e(1,i) + f(2,i)*e(2,i) + f(3,i)*e(3,i)
3303  f(1,i) = f(1,i) - ap*e(1,i)
3304  f(2,i) = f(2,i) - ap*e(2,i)
3305  f(3,i) = f(3,i) - ap*e(3,i)
3306  enddo
3307 
3308  end subroutine project_sphere_v
3309 
3310 #ifdef TO_DO_MQ
3311  subroutine init_mq(phis, gridstruct, npx, npy, is, ie, js, je, ng)
3312  integer, intent(in):: npx, npy, is, ie, js, je, ng
3313  real, intent(in):: phis(is-ng:ie+ng, js-ng:je+ng)
3314  type(fv_grid_type), intent(IN), target :: gridstruct
3315 
3316 ! local:
3317  real zs(is-ng:ie+ng, js-ng:je+ng)
3318  real zb(is-ng:ie+ng, js-ng:je+ng)
3319  real pdx(3,is:ie,js:je+1)
3320  real pdy(3,is:ie+1,js:je)
3321  integer i, j, n
3322 
3323  real, pointer :: rarea(:,:)
3324  real, pointer, dimension(:,:) :: dx, dy
3325  real, pointer, dimension(:,:,:) :: en1, en2, agrid, vlon, vlat
3326 
3327  rarea => gridstruct%rarea
3328  dx => gridstruct%dx
3329  dy => gridstruct%dy
3330  en1 => gridstruct%en1
3331  en2 => gridstruct%en2
3332  agrid => gridstruct%agrid
3333  vlon => gridstruct%vlon
3334  vlat => gridstruct%vlat
3335 
3336 ! do j=js,je
3337 ! do i=is,ie
3338  do j=js-ng,je+ng
3339  do i=is-ng,ie+ng
3340  zs(i,j) = phis(i,j) / grav
3341  enddo
3342  enddo
3343 ! call mpp_update_domains( zs, domain )
3344 
3345 ! call a2b_ord2(zs, zb, gridstruct, npx, npy, is, ie, js, je, ng)
3346  call a2b_ord4(zs, zb, gridstruct, npx, npy, is, ie, js, je, ng)
3347 
3348  do j=js,je+1
3349  do i=is,ie
3350  do n=1,3
3351  pdx(n,i,j) = 0.5*(zb(i,j)+zb(i+1,j))*dx(i,j)*en1(n,i,j)
3352  enddo
3353  enddo
3354  enddo
3355  do j=js,je
3356  do i=is,ie+1
3357  do n=1,3
3358  pdy(n,i,j) = 0.5*(zb(i,j)+zb(i,j+1))*dy(i,j)*en2(n,i,j)
3359  enddo
3360  enddo
3361  enddo
3362 
3363 ! Compute "volume-mean" gradient by Green's theorem
3364  do j=js,je
3365  do i=is,ie
3366  idiag%zxg(i,j) = vlon(i,j,1)*(pdx(1,i,j+1)-pdx(1,i,j)-pdy(1,i,j)+pdy(1,i+1,j)) &
3367  + vlon(i,j,2)*(pdx(2,i,j+1)-pdx(2,i,j)-pdy(2,i,j)+pdy(2,i+1,j)) &
3368  + vlon(i,j,3)*(pdx(3,i,j+1)-pdx(3,i,j)-pdy(3,i,j)+pdy(3,i+1,j))
3369 ! dF/d(lamda) = radius*cos(agrid(i,j,2)) * dF/dx, F is a scalar
3370 ! ________________________
3371  idiag%zxg(i,j) = idiag%zxg(i,j)*rarea(i,j) * radius*cos(agrid(i,j,2))
3372 ! ^^^^^^^^^^^^^^^^^^^^^^^^
3373  enddo
3374  enddo
3375 
3376  end subroutine init_mq
3377 #endif
3378 
3379  end module fv_grid_utils_mod
fv_grid_utils_mod::mid_pt_cart
subroutine, public mid_pt_cart(p1, p2, e3)
Definition: fv_grid_utils.F90:1999
fv_grid_utils_mod::cart_to_latlon
subroutine, public cart_to_latlon(np, q, xs, ys)
Definition: fv_grid_utils.F90:1716
fv_grid_utils_mod::elgs
subroutine elgs(a, n, indx)
The subroutine &#39;elgs&#39; performs the partial-pivoting gaussian elimination.
Definition: fv_grid_utils.F90:3217
fv_grid_utils_mod::global_mx_c
subroutine global_mx_c(q, i1, i2, j1, j2, qmin, qmax)
The subroutine &#39;global_mx_c&#39; computes the global max and min at the cell corners. ...
Definition: fv_grid_utils.F90:3003
fv_grid_utils_mod::mid_pt_sphere
subroutine, public mid_pt_sphere(p1, p2, pm)
Definition: fv_grid_utils.F90:1954
fv_grid_utils_mod::cos_angle
real(kind=r_grid) function, public cos_angle(p1, p2, p3)
Definition: fv_grid_utils.F90:2855
fv_grid_utils_mod::init_cubed_to_latlon
subroutine init_cubed_to_latlon(gridstruct, hydrostatic, agrid, grid_type, ord, bd)
Definition: fv_grid_utils.F90:2281
fv_grid_utils_mod::get_unit_vect3
subroutine get_unit_vect3(p1, p2, uc)
Definition: fv_grid_utils.F90:1838
fv_grid_utils_mod::global_qsum
real function, public global_qsum(p, ifirst, ilast, jfirst, jlast)
The function &#39;global_qsum&#39; computes the quick global sum without area weighting.
Definition: fv_grid_utils.F90:2954
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_grid_utils_mod::tiny_number
real, parameter tiny_number
Definition: fv_grid_utils.F90:105
fv_arrays_mod::fv_grid_bounds_type
Definition: fv_arrays.F90:1137
fv_grid_utils_mod::normalize_vect
subroutine, public normalize_vect(e)
The subroutine &#39;normalize_vect&#39; makes &#39;e&#39; a unit vector.
Definition: fv_grid_utils.F90:1852
fv_grid_utils_mod::init_mq
subroutine init_mq(phis, gridstruct, npx, npy, is, ie, js, je, ng)
Definition: fv_grid_utils.F90:3312
fv_grid_utils_mod::intersect_cross
subroutine intersect_cross(a1, a2, b1, b2, radius, x_inter, local_a, local_b)
The subroutine &#39;intersect_cross&#39; calculates the intersection of two great circles.
Definition: fv_grid_utils.F90:2168
fv_grid_utils_mod::cell_center3
subroutine cell_center3(p1, p2, p3, p4, ec)
The subroutine &#39;cell_center3&#39; gets the center position of a cell.
Definition: fv_grid_utils.F90:2688
fv_eta_mod::set_eta
subroutine, public set_eta(km, ks, ptop, ak, bk)
This is the version of set_eta used in fvGFS and AM4.
Definition: fv_eta.F90:420
fv_grid_utils_mod::intp_great_circle
subroutine, public intp_great_circle(beta, p1, p2, x_o, y_o)
Definition: fv_grid_utils.F90:1868
fv_grid_utils_mod::invert_matrix
subroutine invert_matrix(n, a, x)
Definition: fv_grid_utils.F90:3173
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_utils_mod::latlon2xyz2
subroutine latlon2xyz2(lon, lat, p3)
Definition: fv_grid_utils.F90:1611
fv_grid_utils_mod::efactor_a2c_v
subroutine efactor_a2c_v(edge_vect_s, edge_vect_n, edge_vect_w, edge_vect_e, non_ortho, grid, agrid, npx, npy, nested, bd, regional)
The subroutine &#39;efactor_a2c_v&#39; initializes interpolation factors at face edges for interpolating vect...
Definition: fv_grid_utils.F90:989
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::g_sum
real function, public g_sum(domain, p, ifirst, ilast, jfirst, jlast, ngc, area, mode, reproduce)
The function &#39;g_sum&#39; is the fast version of &#39;globalsum&#39;.
Definition: fv_grid_utils.F90:2902
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
fv_grid_utils_mod::get_latlon_vector
subroutine, public get_latlon_vector(pp, elon, elat)
Definition: fv_grid_utils.F90:3274
fv_grid_utils_mod::get_unit_vect2
subroutine, public get_unit_vect2(e1, e2, uc)
Definition: fv_grid_utils.F90:1821
fv_grid_utils_mod::gnomonic_angl
subroutine gnomonic_angl(im, lamda, theta)
Definition: fv_grid_utils.F90:1514
fv_arrays_mod::r_grid
integer, parameter, public r_grid
Definition: fv_arrays.F90:35
fv_grid_utils_mod::make_eta_level
subroutine, public make_eta_level(km, pe, area, kks, ak, bk, ptop, domain, bd)
Definition: fv_grid_utils.F90:3105
external_sst_mod
Definition: external_sst.F90:21
fv_grid_utils_mod::c2l_ord4
subroutine c2l_ord4(u, v, ua, va, gridstruct, npx, npy, km, grid_type, domain, nested, mode, bd)
Definition: fv_grid_utils.F90:2367
fv_timing_mod
The module &#39;fv_timing&#39; contains FV3 timers.
Definition: fv_timing.F90:24
fv_grid_utils_mod::unit_vect_latlon
subroutine, public unit_vect_latlon(pp, elon, elat)
Definition: fv_grid_utils.F90:2246
fv_grid_utils_mod::mirror_latlon
subroutine mirror_latlon(lon1, lat1, lon2, lat2, lon0, lat0, lon3, lat3)
The subroutine &#39;mirror_latlon&#39; computes the mirror image of (lon0, lat0) as (lon3, lat3) given the "mirror" as defined by (lon1, lat1), (lon2, lat2), and center of the sphere.
Definition: fv_grid_utils.F90:1686
fv_grid_utils_mod::global_mx
subroutine, public global_mx(q, n_g, qmin, qmax, bd)
Definition: fv_grid_utils.F90:2974
fv_grid_utils_mod::ptop_min
real, parameter, public ptop_min
Definition: fv_grid_utils.F90:109
fv_grid_utils_mod::f_p
integer, parameter, public f_p
Definition: fv_grid_utils.F90:102
fv_grid_utils_mod::gnomonic_ed
subroutine gnomonic_ed(im, lamda, theta)
The subroutine &#39;gnomonic_ed&#39; computes the equal distances along the 4 edges of the cubed sphere...
Definition: fv_grid_utils.F90:1309
fv_grid_utils_mod::symm_grid
logical, public symm_grid
Definition: fv_grid_utils.F90:96
fv_grid_utils_mod::symm_ed
subroutine symm_ed(im, lamda, theta)
Definition: fv_grid_utils.F90:1570
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_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::gnomonic_ed_limited
subroutine gnomonic_ed_limited(im, in, nghost, lL, lR, uL, uR, lamda, theta)
The subroutine &#39;gnomonic_ed_limited&#39; creates a limited-area equidistant gnomonic grid with corners gi...
Definition: fv_grid_utils.F90:1405
fv_grid_utils_mod::cell_center2
subroutine, public cell_center2(q1, q2, q3, q4, e2)
Definition: fv_grid_utils.F90:2660
fv_eta_mod
The module &#39;fv_eta&#39; contains routine to set up the reference (Eulerian) pressure coordinate.
Definition: fv_eta.F90:25
fv_grid_utils_mod::v_prod
real(kind=r_grid) function, public v_prod(v1, v2)
Definition: fv_grid_utils.F90:2273
fv_grid_utils_mod::grid_utils_init
subroutine, public grid_utils_init(Atm, npx, npy, npz, non_ortho, grid_type, c2l_order)
Definition: fv_grid_utils.F90:134
fv_grid_utils_mod::vect_cross
subroutine, public vect_cross(e, p1, p2)
The subroutine &#39;vect_cross performs cross products of 3D vectors: e = P1 X P2.
Definition: fv_grid_utils.F90:1758
fv_timing_mod::timing_on
subroutine timing_on(blk_name)
The subroutine &#39;timing_on&#39; starts a timer.
Definition: fv_timing.F90:116
fv_grid_utils_mod::project_sphere_v
subroutine, public project_sphere_v(np, f, e)
Definition: fv_grid_utils.F90:3293
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_utils_mod::c2l_ord2
subroutine, public c2l_ord2(u, v, ua, va, gridstruct, km, grid_type, bd, do_halo)
Definition: fv_grid_utils.F90:2507
get_nearest
subroutine get_nearest()
Definition: fv_grid_utils.F90:2119
fv_grid_utils_mod::intersect
subroutine intersect(a1, a2, b1, b2, radius, x_inter, local_a, local_b)
The subroutine &#39;intersect&#39; calculates the intersection of two great circles.
Definition: fv_grid_utils.F90:2076
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_grid_utils_mod::dist2side
real(kind=r_grid) function dist2side(v1, v2, point)
The function &#39;dist2side&#39; calculates the shortest normalized distance on a sphere from point to straig...
Definition: fv_grid_utils.F90:2756
check_local
subroutine check_local(x1, x2, local)
Definition: fv_grid_utils.F90:2133
fv_grid_utils_mod::fill_ghost_r4
subroutine fill_ghost_r4(q, npx, npy, value, bd)
Definition: fv_grid_utils.F90:3024
fv_grid_utils_mod::edge_factors
subroutine edge_factors(edge_s, edge_n, edge_w, edge_e, non_ortho, grid, agrid, npx, npy, bd)
The subroutine &#39;edge_factors&#39; initializes interpolation factors at face edges for interpolation from ...
Definition: fv_grid_utils.F90:1166
fv_grid_utils_mod::gnomonic_dist
subroutine gnomonic_dist(im, lamda, theta)
Definition: fv_grid_utils.F90:1541
fv_grid_utils_mod::mirror_xyz
subroutine mirror_xyz(p1, p2, p0, p)
The subroutine &#39;mirror_xyz&#39; computes the mirror image of p0(x0, y0, z0) as p(x, y, z) given the "mirror" as defined by p1(x1, y1, z1), p2(x2, y2, z2), and the center of the sphere.
Definition: fv_grid_utils.F90:1651
fv_arrays_mod::fv_atmos_type
Definition: fv_arrays.F90:1166
fv_grid_utils_mod::fill_ghost_r8
subroutine fill_ghost_r8(q, npx, npy, value, bd)
Definition: fv_grid_utils.F90:3064
fv_grid_utils_mod::get_center_vect
subroutine get_center_vect(npx, npy, pp, u1, u2, bd)
Definition: fv_grid_utils.F90:1768
fv_grid_utils_mod::grid_utils_end
subroutine, public grid_utils_end
Definition: fv_grid_utils.F90:906
fv_grid_utils_mod::expand_cell
subroutine, public expand_cell(q1, q2, q3, q4, a1, a2, a3, a4, fac)
Definition: fv_grid_utils.F90:2591
fv_grid_utils_mod::cubed_to_latlon
subroutine, public cubed_to_latlon(u, v, ua, va, gridstruct, npx, npy, km, mode, grid_type, domain, nested, c2l_ord, bd)
Definition: fv_grid_utils.F90:2346
fv_grid_utils_mod::gnomonic_grids
subroutine, public gnomonic_grids(grid_type, im, lon, lat)
Definition: fv_grid_utils.F90:1274
fv_grid_utils_mod::radius
real(kind=r_grid) radius
Definition: fv_grid_utils.F90:107
fv_grid_utils_mod::great_circle_dist_cart
real(kind=r_grid) function great_circle_dist_cart(v1, v2, radius)
The function &#39;great_circle_dist_cart&#39; calculates the normalized great circle distance between &#39;v1&#39; an...
Definition: fv_grid_utils.F90:2041
fv_grid_utils_mod::mid_pt3_cart
subroutine mid_pt3_cart(p1, p2, e)
Definition: fv_grid_utils.F90:1969
fv_grid_utils_mod::latlon2xyz
subroutine, public latlon2xyz(p, e, id)
The subroutine &#39;latlon2xyz&#39; maps (lon, lat) to (x,y,z)
Definition: fv_grid_utils.F90:1622