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