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