FMS  2025.04
Flexible Modeling System
diag_grid.F90
1 !***********************************************************************
2 !* Apache License 2.0
3 !*
4 !* This file is part of the GFDL Flexible Modeling System (FMS).
5 !*
6 !* Licensed under the Apache License, Version 2.0 (the "License");
7 !* you may not use this file except in compliance with the License.
8 !* You may obtain a copy of the License at
9 !*
10 !* http://www.apache.org/licenses/LICENSE-2.0
11 !*
12 !* FMS is distributed in the hope that it will be useful, but WITHOUT
13 !* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied;
14 !* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
15 !* PARTICULAR PURPOSE. See the License for the specific language
16 !* governing permissions and limitations under the License.
17 !***********************************************************************
18 !> @defgroup diag_grid_mod diag_grid_mod
19 !> @ingroup diag_manager
20 !> @brief diag_grid_mod is a set of procedures to work with the
21 !! model's global grid to allow regional output.
22 !!
23 !> @author Seth Underwood seth.underwood@noaa.gov
24 !!
25 !! <TT>diag_grid_mod</TT> contains useful utilities for dealing
26 !! with, mostly, regional output for grids other than the standard
27 !! lat/lon grid. This module contains three public procedures <TT>
28 !! diag_grid_init</TT>, which is shared globably in the <TT>
29 !! diag_manager_mod</TT>, <TT>diag_grid_end</TT> which will free
30 !! up memory used during the register field calls, and
31 !! <TT>get_local_indexes</TT>. The <TT>send_global_grid</TT>
32 !! procedure is called by the model that creates the global grid.
33 !! <TT>send_global_grid</TT> needs to be called before any fields
34 !! are registered that will output only regions. <TT>get_local_indexes</TT>
35 !! is to be called by the <TT>diag_manager_mod</TT> to discover the
36 !! global indexes defining a subregion on the tile.
37 
38 MODULE diag_grid_mod
39 use platform_mod
40 
41  ! <INFO>
42  ! <FUTURE>
43  ! Multi-tile regional output in the cubed sphere.
44  ! </FUTURE>
45  ! <FUTURE>
46  ! Single grid in the tri-polar grid.
47  ! </FUTURE>
48  ! <FUTURE>
49  ! Multi-tile regional output in the tri-polar grid.
50  ! </FUTURE>
51  ! <FUTURE>
52  ! Regional output using array masking. This should allow
53  ! regional output to work on any current or future grid.
54  ! </FUTURE>
55  ! </INFO>
56 
57  USE constants_mod, ONLY: deg_to_rad, rad_to_deg, radius
58  USE fms_mod, ONLY: write_version_number, error_mesg, warning, fatal,&
59  & mpp_pe
60  USE mpp_mod, ONLY: mpp_root_pe, mpp_npes, mpp_max, mpp_min
61  USE mpp_domains_mod, ONLY: domain2d, mpp_get_tile_id,&
63 
64  IMPLICIT NONE
65 
66  ! Include variable "version" to be written to log file.
67 #include<file_version.h>
68 
69  !> @brief Private type to hold the model's global grid data, and other grid information for use
70  !! in this module.
71  !> @ingroup diag_grid_mod
72  type, private :: diag_global_grid_type
73  REAL, allocatable, DIMENSION(:,:) :: glo_lat !< The latitude values on the global grid.
74  REAL, allocatable, DIMENSION(:,:) :: glo_lon !< The longitude values on the global grid.
75  REAL, allocatable, DIMENSION(:,:) :: aglo_lat !< The latitude values on the global a-grid. Here
76  !! we expect isc-1:iec+1 and
77  !! jsc=1:jec+1 to be passed in.
78  REAL, allocatable, DIMENSION(:,:) :: aglo_lon !< The longitude values on the global a-grid. Here
79  !! we expec isc-1:iec+j and
80  !! jsc-1:jec+1 to be passed in.
81  INTEGER :: myxbegin !< The starting index of the compute domain on the current PE.
82  INTEGER :: myybegin !< The starting index of the compute domain on the current PE.
83  INTEGER :: dimi !< The dimension of the global grid in the 'i' / longitudal direction.
84  INTEGER :: dimj !< The dimension of the global grid in the 'j' / latitudal direction.
85  INTEGER :: adimi !< The dimension of the global a-grid in the 'i' / longitudal direction. Again,
86  !! the expected dimension for diag_grid_mod is isc-1:iec+1.
87  INTEGER :: adimj !< The dimension of the global a-grid in the 'j' / latitudal direction. Again,
88  !! the expected dimension for diag_grid_mod is jsc-1:jec+1.
89  INTEGER :: tile_number !< The tile the <TT>glo_lat</TT> and <TT>glo_lon</TT> define.
90  INTEGER :: ntiles !< The number of tiles.
91  INTEGER :: pestart !< The starting PE number for the current tile.
92  INTEGER :: peend !< The ending PE number for the current tile.
93  CHARACTER(len=128) :: grid_type !< The global grid type.
94  END TYPE diag_global_grid_type
95 
96  !> @brief Private type to hold the corresponding (x,y,z) location for a (lat,lon)
97  !! location.
98  !> @ingroup diag_grid_mod
99  type, private :: point
100  REAL :: x !< The x value of the (x,y,z) coordinates.
101  REAL :: y !< The y value of the (x,y,z) coordinates.
102  REAL :: z !< The z value of the (x,y,z) coordinates.
103  END TYPE point
104 
105 !> @addtogroup diag_grid_mod
106 !> @{
107 
108  TYPE(diag_global_grid_type) :: diag_global_grid !< Variable to hold the global grid data
109 
110  LOGICAL :: diag_grid_initialized = .false. !< Indicates if the diag_grid_mod has been initialized.
111 
112  PRIVATE
115 
116 CONTAINS
117 
118  !> @brief Send the global grid to the <TT>diag_manager_mod</TT> for
119  !! regional output.
120  !!
121  !> In order for the diag_manager to do regional output for grids
122  !! other than the standard lat/lon grid, the <TT>
123  !! diag_manager_mod</TT> needs to know the latitude and
124  !! longitude values for the entire global grid. This procedure
125  !! is the mechanism the models will use to share their grid with
126  !! the diagnostic manager.
127  !!
128  !! This procedure needs to be called after the grid is created,
129  !! and before the first call to register the fields.
130  SUBROUTINE diag_grid_init(domain, glo_lat, glo_lon, aglo_lat, aglo_lon)
131  TYPE(domain2d), INTENT(in) :: domain !< The domain to which the grid data corresponds.
132  CLASS(*), INTENT(in), DIMENSION(:,:) :: glo_lat !< The latitude information for the grid tile.
133  CLASS(*), INTENT(in), DIMENSION(:,:) :: glo_lon !< The longitude information for the grid tile.
134  CLASS(*), INTENT(in), DIMENSION(:,:) :: aglo_lat !< The latitude information for the a-grid tile.
135  CLASS(*), INTENT(in), DIMENSION(:,:) :: aglo_lon !< The longitude information for the a-grid tile.
136 
137  INTEGER, DIMENSION(1) :: tile
138  INTEGER :: ntiles
139  INTEGER :: stat
140  INTEGER :: i_dim, j_dim
141  INTEGER :: ai_dim, aj_dim
142  INTEGER, DIMENSION(2) :: latdim, londim
143  INTEGER, DIMENSION(2) :: alatdim, alondim
144  INTEGER :: mype, npes, npespertile
145  INTEGER, ALLOCATABLE, DIMENSION(:) :: xbegin, xend, ybegin, yend
146 
147  ! Write the file version to the logfile
148  CALL write_version_number("DIAG_GRID_MOD", version)
149 
150  ! Verify all allocatable / pointers for diag_global_grid hare not
151  ! allocated / associated.
152  IF ( ALLOCATED(xbegin) ) DEALLOCATE(xbegin)
153  IF ( ALLOCATED(ybegin) ) DEALLOCATE(ybegin)
154  IF ( ALLOCATED(xend) ) DEALLOCATE(xend)
155  IF ( ALLOCATED(yend) ) DEALLOCATE(yend)
156 
157  ! What is my PE
158  mype = mpp_pe() -mpp_root_pe() + 1
159 
160  ! Get the domain/pe layout, and allocate the [xy]begin|end arrays/pointers
161  npes = mpp_npes()
162  ALLOCATE(xbegin(npes), &
163  & ybegin(npes), &
164  & xend(npes), &
165  & yend(npes), stat=stat)
166  IF ( stat .NE. 0 ) THEN
167  CALL error_mesg('diag_grid_mod::diag_grid_init',&
168  &'Could not allocate memory for the compute grid indices&
169  &.', fatal)
170  END IF
171 
172  ! Get tile information
173  ntiles = mpp_get_ntile_count(domain)
174  tile = mpp_get_tile_id(domain)
175 
176  ! Number of PEs per tile
177  npespertile = npes / ntiles
178  diag_global_grid%peEnd = npespertile * tile(1)
179  diag_global_grid%peStart = diag_global_grid%peEnd - npespertile + 1
180 
181  ! Get the compute domains
182  CALL mpp_get_compute_domains(domain,&
183  & xbegin=xbegin, xend=xend,&
184  & ybegin=ybegin, yend=yend)
185 
186  ! Module initialized
187  diag_grid_initialized = .true.
188 
189  ! Get the size of the grids
190  latdim = shape(glo_lat)
191  londim = shape(glo_lon)
192  IF ( (latdim(1) == londim(1)) .AND.&
193  &(latdim(2) == londim(2)) ) THEN
194  i_dim = latdim(1)
195  j_dim = latdim(2)
196  ELSE
197  CALL error_mesg('diag_grid_mod::diag_grid_init',&
198  &'glo_lat and glo_lon must be the same shape.', fatal)
199  END IF
200 
201  ! Same thing for the a-grid
202  alatdim = shape(aglo_lat)
203  alondim = shape(aglo_lon)
204  IF ( (alatdim(1) == alondim(1)) .AND. &
205  &(alatdim(2) == alondim(2)) ) THEN
206  IF ( tile(1) == 4 .OR. tile(1) == 5 ) THEN
207  ! These tiles need to be transposed.
208  ai_dim = alatdim(2)
209  aj_dim = alatdim(1)
210  ELSE
211  ai_dim = alatdim(1)
212  aj_dim = alatdim(2)
213  END IF
214  ELSE
215  CALL error_mesg('diag_grid_mod::diag_grid_init',&
216  & "a-grid's glo_lat and glo_lon must be the same shape.", fatal)
217  END IF
218 
219  ! Allocate the grid arrays
220  IF ( allocated(diag_global_grid%glo_lat) .OR.&
221  & allocated(diag_global_grid%glo_lon) ) THEN
222  IF ( mpp_pe() == mpp_root_pe() ) &
223  & CALL error_mesg('diag_grid_mod::diag_grid_init',&
224  &'The global grid has already been initialized', warning)
225  ELSE
226  ALLOCATE(diag_global_grid%glo_lat(i_dim,j_dim),&
227  & diag_global_grid%glo_lon(i_dim,j_dim), stat=stat)
228  IF ( stat .NE. 0 ) THEN
229  CALL error_mesg('diag_grid_mod::diag_grid_init',&
230  &'Could not allocate memory for the global grid.', fatal)
231  END IF
232  END IF
233 
234  ! Same thing for the a-grid
235  IF ( allocated(diag_global_grid%aglo_lat) .OR.&
236  & allocated(diag_global_grid%aglo_lon) ) THEN
237  IF ( mpp_pe() == mpp_root_pe() ) &
238  & CALL error_mesg('diag_grid_mod::diag_grid_init',&
239  &'The global a-grid has already been initialized', warning)
240  ELSE
241  ALLOCATE(diag_global_grid%aglo_lat(0:ai_dim-1,0:aj_dim-1),&
242  & diag_global_grid%aglo_lon(0:ai_dim-1,0:aj_dim-1), stat=stat)
243  IF ( stat .NE. 0 ) THEN
244  CALL error_mesg('diag_global_mod::diag_grid_init',&
245  &'Could not allocate memory for the global a-grid', fatal)
246  END IF
247  END IF
248 
249  ! Set the values for diag_global_grid
250 
251  ! If we are on tile 4 or 5, we need to transpose the grid to get
252  ! this to work.
253  IF ( tile(1) == 4 .OR. tile(1) == 5 ) THEN
254  SELECT TYPE (aglo_lat)
255  TYPE IS (real(kind=r4_kind))
256  diag_global_grid%aglo_lat = transpose(aglo_lat)
257  TYPE IS (real(kind=r8_kind))
258  diag_global_grid%aglo_lat = transpose(real(aglo_lat))
259  CLASS DEFAULT
260  CALL error_mesg('diag_grid_mod::diag_grid_init',&
261  & 'The a-grid latitude data is not one of the supported types of real(kind=4) or real(kind=8)', fatal)
262  END SELECT
263 
264  SELECT TYPE (aglo_lon)
265  TYPE IS (real(kind=r4_kind))
266  diag_global_grid%aglo_lon = transpose(aglo_lon)
267  TYPE IS (real(kind=r8_kind))
268  diag_global_grid%aglo_lon = transpose(real(aglo_lon))
269  CLASS DEFAULT
270  CALL error_mesg('diag_grid_mod::diag_grid_init',&
271  & 'The a-grid longitude data is not one of the supported types of real(kind=4) or real(kind=8)', fatal)
272  END SELECT
273  ELSE
274  SELECT TYPE (aglo_lat)
275  TYPE IS (real(kind=r4_kind))
276  diag_global_grid%aglo_lat = aglo_lat
277  TYPE IS (real(kind=r8_kind))
278  diag_global_grid%aglo_lat = real(aglo_lat)
279  CLASS DEFAULT
280  CALL error_mesg('diag_grid_mod::diag_grid_init',&
281  & 'The a-grid latitude data is not one of the supported types of real(kind=4) or real(kind=8)', fatal)
282  END SELECT
283 
284  SELECT TYPE (aglo_lon)
285  TYPE IS (real(kind=r4_kind))
286  diag_global_grid%aglo_lon = aglo_lon
287  TYPE IS (real(kind=r8_kind))
288  diag_global_grid%aglo_lon = real(aglo_lon)
289  CLASS DEFAULT
290  CALL error_mesg('diag_grid_mod::diag_grid_init',&
291  & 'The a-grid longitude data is not one of the supported types of real(kind=4) or real(kind=8)', fatal)
292  END SELECT
293  END IF
294 
295  SELECT TYPE (glo_lat)
296  TYPE IS (real(kind=r4_kind))
297  diag_global_grid%glo_lat = glo_lat
298  TYPE IS (real(kind=r8_kind))
299  diag_global_grid%glo_lat = real(glo_lat)
300  CLASS DEFAULT
301  CALL error_mesg('diag_grid_mod::diag_grid_init',&
302  & 'The grid latitude data is not one of the supported types of real(kind=4) or real(kind=8)', fatal)
303  END SELECT
304 
305  SELECT TYPE (glo_lon)
306  TYPE IS (real(kind=r4_kind))
307  diag_global_grid%glo_lon = glo_lon
308  TYPE IS (real(kind=r8_kind))
309  diag_global_grid%glo_lon = real(glo_lon)
310  CLASS DEFAULT
311  CALL error_mesg('diag_grid_mod::diag_grid_init',&
312  & 'The grid longitude data is not one of the supported types of real(kind=4) or real(kind=8)', fatal)
313  END SELECT
314 
315  diag_global_grid%dimI = i_dim
316  diag_global_grid%dimJ = j_dim
317  diag_global_grid%adimI = ai_dim
318  diag_global_grid%adimJ = aj_dim
319  !--- For the nested model, the nested region only has 1 tile ( ntiles = 1) but
320  !--- the tile_id is 7 for the nested region. In the routine get_local_indexes,
321  !--- local variables ijMin and ijMax have dimesnion (ntiles) and will access
322  !--- ijMin(diag_global_grid%tile_number,:). For the nested region, ntiles = 1 and
323  !--- diag_global_grid%tile_number = 7 will cause out of bounds. So need to
324  !--- set diag_global_grid%tile_number = 1 when ntiles = 1 for the nested model.
325  if(ntiles == 1) then
326  diag_global_grid%tile_number = 1
327  else
328  diag_global_grid%tile_number = tile(1)
329  endif
330  diag_global_grid%ntiles = ntiles
331  diag_global_grid%myXbegin = xbegin(mype)
332  diag_global_grid%myYbegin = ybegin(mype)
333 
334  ! Unallocate arrays used here
335  DEALLOCATE(xbegin)
336  DEALLOCATE(ybegin)
337  DEALLOCATE(xend)
338  DEALLOCATE(yend)
339  END SUBROUTINE diag_grid_init
340 
341  !> @brief Unallocate the diag_global_grid variable.
342  !!
343  !> The <TT>diag_global_grid</TT> variable is only needed during
344  !! the register field calls, and then only if there are fields
345  !! requestion regional output. Once all the register fields
346  !! calls are complete (before the first <TT>send_data</TT> call
347  !! this procedure can be called to free up memory.
348  SUBROUTINE diag_grid_end()
349 
350  IF ( diag_grid_initialized ) THEN
351  ! De-allocate grid
352  IF ( allocated(diag_global_grid%glo_lat) ) THEN
353  DEALLOCATE(diag_global_grid%glo_lat)
354  ELSE IF ( mpp_pe() == mpp_root_pe() ) THEN
355  CALL error_mesg('diag_grid_mod::diag_grid_end',&
356  &'diag_global_grid%glo_lat was not allocated.', warning)
357  END IF
358 
359  IF ( allocated(diag_global_grid%glo_lon) ) THEN
360  DEALLOCATE(diag_global_grid%glo_lon)
361  ELSE IF ( mpp_pe() == mpp_root_pe() ) THEN
362  CALL error_mesg('diag_grid_mod::diag_grid_end',&
363  &'diag_global_grid%glo_lon was not allocated.', warning)
364  END IF
365  ! De-allocate a-grid
366  IF ( allocated(diag_global_grid%aglo_lat) ) THEN
367  DEALLOCATE(diag_global_grid%aglo_lat)
368  ELSE IF ( mpp_pe() == mpp_root_pe() ) THEN
369  CALL error_mesg('diag_grid_mod::diag_grid_end',&
370  &'diag_global_grid%aglo_lat was not allocated.', warning)
371  END IF
372 
373  IF ( allocated(diag_global_grid%aglo_lon) ) THEN
374  DEALLOCATE(diag_global_grid%aglo_lon)
375  ELSE IF ( mpp_pe() == mpp_root_pe() ) THEN
376  CALL error_mesg('diag_grid_mod::diag_grid_end',&
377  &'diag_global_grid%aglo_lon was not allocated.', warning)
378  END IF
379 
380  diag_grid_initialized = .false.
381  END IF
382  END SUBROUTINE diag_grid_end
383 
384  !> @brief Find the local start and local end indexes on the local PE
385  !! for regional output.
386  !!
387  !> Given a defined region, find the local indexes on the local
388  !! PE surrounding the region.
389  SUBROUTINE get_local_indexes(latStart, latEnd, lonStart, lonEnd,&
390  & istart, iend, jstart, jend)
391  REAL, INTENT(in) :: latstart !< lat start angles
392  REAL, INTENT(in) :: lonstart !< lon start angles
393  REAL, INTENT(in) :: latend !< lat end angles
394  REAL, INTENT(in) :: lonend !< lon end angles
395  INTEGER, INTENT(out) :: istart !< i start indexes
396  INTEGER, INTENT(out) :: jstart !< j start indexes
397  INTEGER, INTENT(out) :: iend !< i end indexes
398  INTEGER, INTENT(out) :: jend !< j end indexes
399 
400  REAL, ALLOCATABLE, DIMENSION(:,:) :: delta_lat, delta_lon, grid_lon
401 
402  REAL, DIMENSION(4) :: dists_lon, dists_lat
403  REAL :: lonendadj, my_lonstart, my_lonend
404 
405  INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ijmin, ijmax
406  INTEGER :: mytile, ntiles, i, j, k, dimi, dimj, istat
407  INTEGER :: count
408 
409  LOGICAL :: onmype
410 
411  !For cfsite potential fix.
412  INTEGER :: mini
413  INTEGER :: minj
414  REAL :: minimum_distance
415  REAL :: global_min_distance
416  INTEGER :: rank_buf
417 
418  IF ( .NOT. diag_grid_initialized )&
419  & CALL error_mesg('diag_grid_mod::get_local_indexes',&
420  &'Module not initialized, first initialze module with a call &
421  &to diag_grid_init', fatal)
422 
423  ! Make adjustment for negative longitude values
424  if ( lonstart < 0. ) then
425  my_lonstart = lonstart + 360.
426  else
427  my_lonstart = lonstart
428  end if
429  if ( lonend < 0. ) then
430  my_lonend = lonend + 360.
431  else
432  my_lonend = lonend
433  end if
434 
435  IF (latstart .EQ. latend .AND. my_lonstart .EQ. my_lonend) THEN
436 
437  !For a single point, use the a-grid longitude and latitude
438  !values.
439 
440  mytile = diag_global_grid%tile_number
441  ntiles = diag_global_grid%ntiles
442 
443  allocate(ijmin(ntiles,2))
444  ijmin = 0
445 
446  !Find the i,j indices of the a-grid point nearest to the
447  !my_lonStart,latStart point.
448  CALL find_nearest_agrid_index(latstart, &
449  my_lonstart, &
450  mini, &
451  minj, &
452  minimum_distance)
453 
454  !Find the minimum distance across all ranks.
455  global_min_distance = minimum_distance
456  CALL mpp_min(global_min_distance)
457 
458  !In the case of a tie (i.e. two ranks with exactly the same
459  !minimum distance), use the i,j values from the larger rank id.
460  IF (global_min_distance .EQ. minimum_distance) THEN
461  rank_buf = mpp_pe()
462  ELSE
463  rank_buf = -1
464  ENDIF
465  CALL mpp_max(rank_buf)
466 
467  !Sanity check.
468  IF (rank_buf .EQ. -1) THEN
469  CALL error_mesg("get_local_indexes", &
470  "No rank has minimum distance.", &
471  fatal)
472  ENDIF
473 
474  IF (rank_buf .EQ. mpp_pe()) THEN
475  ijmin(mytile,1) = mini + diag_global_grid%myXbegin - 1
476  ijmin(mytile,2) = minj + diag_global_grid%myYbegin - 1
477  ENDIF
478 
479  DO i = 1,ntiles
480  CALL mpp_max(ijmin(i,1))
481  CALL mpp_max(ijmin(i,2))
482  ENDDO
483 
484  istart = ijmin(mytile,1)
485  jstart = ijmin(mytile,2)
486  iend = istart
487  jend = jstart
488 
489  DEALLOCATE(ijmin)
490  ELSE
491 
492  mytile = diag_global_grid%tile_number
493  ntiles = diag_global_grid%ntiles
494 
495  ! Arrays to home min/max for each tile
496  ALLOCATE(ijmin(ntiles,2), stat=istat)
497  IF ( istat .NE. 0 )&
498  & CALL error_mesg('diag_grid_mod::get_local_indexes',&
499  &'Cannot allocate ijMin index array', fatal)
500  ALLOCATE(ijmax(ntiles,2), stat=istat)
501  IF ( istat .NE. 0 )&
502  & CALL error_mesg('diag_grid_mod::get_local_indexes',&
503  &'Cannot allocate ijMax index array', fatal)
504  ijmin = 0
505  ijmax = 0
506 
507  ! There will be four points to define a region, find all four.
508  ! Need to call the correct function depending on if the tile is a
509  ! pole tile or not.
510  dimi = diag_global_grid%dimI
511  dimj = diag_global_grid%dimJ
512 
513  ! Build the delta array
514  ALLOCATE(delta_lat(dimi,dimj), stat=istat)
515  IF ( istat .NE. 0 )&
516  & CALL error_mesg('diag_grid_mod::get_local_indexes',&
517  &'Cannot allocate latitude delta array', fatal)
518  ALLOCATE(delta_lon(dimi,dimj), stat=istat)
519  IF ( istat .NE. 0 )&
520  & CALL error_mesg('diag_grid_mod::get_local_indexes',&
521  &'Cannot allocate longitude delta array', fatal)
522  DO j=1, dimj
523  DO i=1, dimi
524  count = 0
525  dists_lon = 0.
526  dists_lat = 0.
527  IF ( i < dimi ) THEN
528  dists_lon(1) = abs(diag_global_grid%glo_lon(i+1,j) - diag_global_grid%glo_lon(i,j))
529  dists_lat(1) = abs(diag_global_grid%glo_lat(i+1,j) - diag_global_grid%glo_lat(i,j))
530  count = count+1
531  END IF
532  IF ( j < dimj ) THEN
533  dists_lon(2) = abs(diag_global_grid%glo_lon(i,j+1) - diag_global_grid%glo_lon(i,j))
534  dists_lat(2) = abs(diag_global_grid%glo_lat(i,j+1) - diag_global_grid%glo_lat(i,j))
535  count = count+1
536  END IF
537  IF ( i > 1 ) THEN
538  dists_lon(3) = abs(diag_global_grid%glo_lon(i,j) - diag_global_grid%glo_lon(i-1,j))
539  dists_lat(3) = abs(diag_global_grid%glo_lat(i,j) - diag_global_grid%glo_lat(i-1,j))
540  count = count+1
541  END IF
542  IF ( j > 1 ) THEN
543  dists_lon(4) = abs(diag_global_grid%glo_lon(i,j) - diag_global_grid%glo_lon(i,j-1))
544  dists_lat(4) = abs(diag_global_grid%glo_lat(i,j) - diag_global_grid%glo_lat(i,j-1))
545  count = count+1
546  END IF
547 
548  ! Fix wrap around problem
549  DO k=1, 4
550  IF ( dists_lon(k) > 180.0 ) THEN
551  dists_lon(k) = 360.0 - dists_lon(k)
552  END IF
553  END DO
554  delta_lon(i,j) = sum(dists_lon)/real(count)
555  delta_lat(i,j) = sum(dists_lat)/real(count)
556  END DO
557  END DO
558 
559  ijmin = huge(1)
560  ijmax = -huge(1)
561 
562  ! Adjusted longitude array
563  ALLOCATE(grid_lon(dimi,dimj), stat=istat)
564  IF ( istat .NE. 0 )&
565  & CALL error_mesg('diag_grid_mod::get_local_indexes',&
566  &'Cannot allocate temporary longitude array', fatal)
567  grid_lon = diag_global_grid%glo_lon
568 
569  ! Make adjustments where required
570  IF ( my_lonstart > my_lonend ) THEN
571  WHERE ( grid_lon < my_lonstart )
572  grid_lon = grid_lon + 360.0
573  END WHERE
574  lonendadj = my_lonend + 360.0
575  ELSE
576  lonendadj = my_lonend
577  END IF
578 
579  DO j=1, dimj-1
580  DO i=1, dimi-1
581  onmype = .false.
582  IF ( latstart-delta_lat(i,j) <= diag_global_grid%glo_lat(i,j) .AND.&
583  & diag_global_grid%glo_lat(i,j) < latend+delta_lat(i,j) ) THEN
584  ! Short-cut for the poles
585  IF ( (abs(latstart)-delta_lat(i,j) <= 90.0 .AND.&
586  & 90.0 <= abs(latend)+delta_lat(i,j)) .AND.&
587  & abs(diag_global_grid%glo_lat(i,j)) == 90.0 ) THEN
588  onmype = .true.
589  ELSE IF ( (my_lonstart-delta_lon(i,j) <= grid_lon(i,j) .AND.&
590  & grid_lon(i,j) < lonendadj+delta_lon(i,j)) ) THEN
591  onmype = .true.
592  ELSE
593  onmype = .false.
594  END IF
595  IF ( onmype ) THEN
596  ijmin(mytile,1) = min(ijmin(mytile,1),i + diag_global_grid%myXbegin - 1)
597  ijmax(mytile,1) = max(ijmax(mytile,1),i + diag_global_grid%myXbegin - 1)
598  ijmin(mytile,2) = min(ijmin(mytile,2),j + diag_global_grid%myYbegin - 1)
599  ijmax(mytile,2) = max(ijmax(mytile,2),j + diag_global_grid%myYbegin - 1)
600  END IF
601  END IF
602  END DO
603  END DO
604  DEALLOCATE(delta_lon)
605  DEALLOCATE(delta_lat)
606  DEALLOCATE(grid_lon)
607 
608  ! Global min/max reduce
609  DO i=1, ntiles
610  CALL mpp_min(ijmin(i,1))
611  CALL mpp_max(ijmax(i,1))
612  CALL mpp_min(ijmin(i,2))
613  CALL mpp_max(ijmax(i,2))
614  END DO
615 
616  IF ( ijmin(mytile,1) == huge(1) .OR. ijmax(mytile,1) == -huge(1) ) THEN
617  ijmin(mytile,1) = 0
618  ijmax(mytile,1) = 0
619  END IF
620  IF ( ijmin(mytile,2) == huge(1) .OR. ijmax(mytile,2) == -huge(1) ) THEN
621  ijmin(mytile,2) = 0
622  ijmax(mytile,2) = 0
623  END IF
624 
625  istart = ijmin(mytile,1)
626  jstart = ijmin(mytile,2)
627  iend = ijmax(mytile,1)
628  jend = ijmax(mytile,2)
629 
630  DEALLOCATE(ijmin)
631  DEALLOCATE(ijmax)
632  END IF
633 
634  END SUBROUTINE get_local_indexes
635 
636  !> @brief Find the indices of the nearest grid point of the a-grid to the
637  !! specified (lon,lat) location on the local PE. if desired point not
638  !! within domain of local PE, return (0,0) as the indices.
639  SUBROUTINE get_local_indexes2(lat, lon, iindex, jindex)
640  REAL, INTENT(in) :: lat !< lat location
641  REAL, INTENT(in) :: lon !< lon location
642  INTEGER, INTENT(out) :: iindex !< i indexes
643  INTEGER, INTENT(out) :: jindex !< j indexes
644 
645  INTEGER :: indexes(2)
646 
647  IF ( .NOT. diag_grid_initialized )&
648  & CALL error_mesg('diag_grid_mod::get_local_indexes2',&
649  &'Module not initialized, first initialze module with a call &
650  &to diag_grid_init', fatal)
651 
652  indexes = 0
653 
654  IF ( mod(diag_global_grid%tile_number,3) == 0 ) THEN
655  IF ( lat > 30.0 .AND. diag_global_grid%tile_number == 3 ) THEN
656  indexes(:) = find_pole_index_agrid(lat,lon)
657  ELSE IF ( lat < -30.0 .AND. diag_global_grid%tile_number == 6 ) THEN
658  indexes(:) = find_pole_index_agrid(lat,lon)
659  ENDIF
660  ELSE
661  indexes(:) = find_equator_index_agrid(lat,lon)
662  END IF
663 
664  iindex = indexes(1)
665  jindex = indexes(2)
666  IF (iindex == diag_global_grid%adimI -1 .OR.&
667  jindex == diag_global_grid%adimJ -1 ) THEN
668  iindex = 0
669  jindex = 0
670  ENDIF
671 
672  END SUBROUTINE get_local_indexes2
673 
674  !> @fn pure elemental real rad2deg(real angle)
675  !> @brief Convert an angle in radian to degrees.
676  !!
677  !> Given a scalar, or an array of angles in radians this
678  !! function will return a scalar or array (of the same
679  !! dimension) of angles in degrees.
680  !! @return Scalar or array (depending on the size of angle) of angles in
681  !! degrees.
682  PURE ELEMENTAL REAL FUNCTION rad2deg(angle)
683  REAL, INTENT(in) :: angle !< Scalar or array of angles in radians.
684 
685  rad2deg = rad_to_deg * angle
686  END FUNCTION rad2deg
687 
688  !> @brief Convert an angle in degrees to radians.
689  !!
690  !> Given a scalar, or an array of angles in degrees this
691  !! function will return a scalar or array (of the same
692  !! dimension) of angles in radians.
693  !! @return Scalar or array (depending on the size of angle) of angles in
694  !! radians.
695  PURE ELEMENTAL REAL FUNCTION deg2rad(angle)
696  REAL, INTENT(in) :: angle !< Scalar or array of angles in degrees.
697 
698  deg2rad = deg_to_rad * angle
699  END FUNCTION deg2rad
700 
701  !> @brief Return the closest index (i,j) to the given (lat,lon) point.
702  !!
703  !> This function searches a pole a-grid tile looking for the grid point
704  !! closest to the give (lat, lon) location, and returns the i
705  !! and j indexes of the point.
706  !! @return The (i, j) location of the closest grid to the given (lat,
707  !! lon) location.
708  PURE FUNCTION find_pole_index_agrid(lat, lon)
709  INTEGER, DIMENSION(2) :: find_pole_index_agrid !< The (i, j) location of the closest grid to the given (lat,
710  !! lon) location.
711  REAL, INTENT(in) :: lat !< Latitude location
712  REAL, INTENT(in) :: lon !< Longitude location
713 
714  INTEGER :: indxi !< Indexes to be returned.
715  INTEGER :: indxj !< Indexes to be returned.
716  INTEGER :: dimi !< Size of the grid dimensions
717  INTEGER :: dimj !< Size of the grid dimensions
718  INTEGER :: i !< Count indexes
719  INTEGER :: j !< Count indexes
720  INTEGER :: nearestcorner !< index of the nearest corner.
721  INTEGER, DIMENSION(4,2) :: ijarray !< indexes of the cornerPts and pntDistances arrays
722  REAL :: llat, llon
723  REAL :: maxctrdist !< maximum distance to center of grid
724  REAL, DIMENSION(4) :: pntdistances !< distance from origPt to corner
725  TYPE(point) :: origpt !< Original point
726  REAL, DIMENSION(4,2) :: cornerpts !< Corner points using (lat,lon)
727  TYPE(point), DIMENSION(9) :: points !< xyz of 8 nearest neighbors
728  REAL, DIMENSION(9) :: distsqrd !< distance between origPt and points(:)
729 
730  ! Set the inital fail values for indxI and indxJ
731  indxi = 0
732  indxj = 0
733 
734  dimi = diag_global_grid%adimI
735  dimj = diag_global_grid%adimJ
736 
737  ! Since the poles have an non-unique longitude value, make a small correction if looking for one of the poles.
738  IF ( lat == 90.0 ) THEN
739  llat = lat - .1
740  ELSE IF ( lat == -90.0 ) THEN
741  llat = lat + .1
742  ELSE
743  llat = lat
744  END IF
745  llon = lon
746 
747  origpt = latlon2xyz(llat,llon)
748 
749  iloop: DO i=0, dimi-2
750  jloop: DO j = 0, dimj-2
751  cornerpts = reshape( (/ diag_global_grid%aglo_lat(i, j), diag_global_grid%aglo_lon(i, j),&
752  & diag_global_grid%aglo_lat(i+1,j+1),diag_global_grid%aglo_lon(i+1,j+1),&
753  & diag_global_grid%aglo_lat(i+1,j), diag_global_grid%aglo_lon(i+1,j),&
754  & diag_global_grid%aglo_lat(i, j+1),diag_global_grid%aglo_lon(i, j+1) /),&
755  & (/ 4, 2 /), order=(/2,1/) )
756  ! Find the maximum half distance of the corner points
757  maxctrdist = max(gcirdistance(cornerpts(1,1),cornerpts(1,2), cornerpts(2,1),cornerpts(2,2)),&
758  & gcirdistance(cornerpts(3,1),cornerpts(3,2), cornerpts(4,1),cornerpts(4,2)))
759 
760  ! Find the distance of the four corner points to the point of interest.
761  pntdistances = gcirdistance(cornerpts(:,1),cornerpts(:,2), llat,llon)
762 
763  IF ( (minval(pntdistances) <= maxctrdist) .AND. (i*j.NE.0) ) THEN
764  ! Set up the i,j index array
765  ijarray = reshape( (/ i, j, i+1, j+1, i+1, j, i, j+1 /), (/ 4, 2 /), order=(/2,1/) )
766 
767  ! the nearest point index
768  nearestcorner = minloc(pntdistances,1)
769 
770  indxi = ijarray(nearestcorner,1)
771  indxj = ijarray(nearestcorner,2)
772 
773  EXIT iloop
774  END IF
775  END DO jloop
776  END DO iloop
777 
778 
779  ! Make sure we have indexes in the correct range
780  valid: IF ( (indxi <= 0 .OR. dimi-1 <= indxi) .OR. &
781  & (indxj <= 0 .OR. dimj-1 <= indxj) ) THEN
782  indxi = 0
783  indxj = 0
784  ELSE ! indxI and indxJ are valid.
785  ! Since we are looking for the closest grid point to the
786  ! (lat,lon) point, we need to check the surrounding
787  ! points. The indexes for the variable points are as follows
788  !
789  ! 1---4---7
790  ! | | |
791  ! 2---5---8
792  ! | | |
793  ! 3---6---9
794 
795  ! Set the 'default' values for points(:) x,y,z to some large
796  ! value.
797  DO i=1, 9
798  points(i)%x = 1.0e20
799  points(i)%y = 1.0e20
800  points(i)%z = 1.0e20
801  END DO
802 
803  ! All the points around the i,j indexes
804  points(1) = latlon2xyz(diag_global_grid%aglo_lat(indxi-1,indxj+1),&
805  & diag_global_grid%aglo_lon(indxi-1,indxj+1))
806  points(2) = latlon2xyz(diag_global_grid%aglo_lat(indxi-1,indxj),&
807  & diag_global_grid%aglo_lon(indxi-1,indxj))
808  points(3) = latlon2xyz(diag_global_grid%aglo_lat(indxi-1,indxj-1),&
809  & diag_global_grid%aglo_lon(indxi-1,indxj-1))
810  points(4) = latlon2xyz(diag_global_grid%aglo_lat(indxi, indxj+1),&
811  & diag_global_grid%aglo_lon(indxi, indxj+1))
812  points(5) = latlon2xyz(diag_global_grid%aglo_lat(indxi, indxj),&
813  & diag_global_grid%aglo_lon(indxi, indxj))
814  points(6) = latlon2xyz(diag_global_grid%aglo_lat(indxi, indxj-1),&
815  & diag_global_grid%aglo_lon(indxi, indxj-1))
816  points(7) = latlon2xyz(diag_global_grid%aglo_lat(indxi+1,indxj+1),&
817  & diag_global_grid%aglo_lon(indxi+1,indxj+1))
818  points(8) = latlon2xyz(diag_global_grid%aglo_lat(indxi+1,indxj),&
819  & diag_global_grid%aglo_lon(indxi+1,indxj))
820  points(9) = latlon2xyz(diag_global_grid%aglo_lat(indxi+1,indxj-1),&
821  & diag_global_grid%aglo_lon(indxi+1,indxj-1))
822 
823 
824  ! Calculate the distance squared between the points(:) and the origPt
825  distsqrd = distancesqrd(origpt, points)
826 
827  SELECT CASE (minloc(distsqrd,1))
828  CASE ( 1 )
829  indxi = indxi-1
830  indxj = indxj+1
831  CASE ( 2 )
832  indxi = indxi-1
833  indxj = indxj
834  CASE ( 3 )
835  indxi = indxi-1
836  indxj = indxj-1
837  CASE ( 4 )
838  indxi = indxi
839  indxj = indxj+1
840  CASE ( 5 )
841  indxi = indxi
842  indxj = indxj
843  CASE ( 6 )
844  indxi = indxi
845  indxj = indxj-1
846  CASE ( 7 )
847  indxi = indxi+1
848  indxj = indxj+1
849  CASE ( 8 )
850  indxi = indxi+1
851  indxj = indxj
852  CASE ( 9 )
853  indxi = indxi+1
854  indxj = indxj-1
855  CASE DEFAULT
856  indxi = 0
857  indxj = 0
858  END SELECT
859  END IF valid
860 
861  ! Set the return value for the function
862  find_pole_index_agrid = (/indxi, indxj/)
863  END FUNCTION find_pole_index_agrid
864 
865  !> @brief Return the closest index (i,j) to the given (lat,lon) point.
866  !!
867  !> This function searches a equator grid tile looking for the grid point
868  !! closest to the give (lat, lon) location, and returns the i
869  !! and j indexes of the point.
870  !! @return The (i, j) location of the closest grid to the given (lat,
871  !! lon) location.
872  PURE FUNCTION find_equator_index_agrid(lat, lon)
873  INTEGER, DIMENSION(2) :: find_equator_index_agrid !< The (i, j) location of the closest grid to the given (lat,
874  !! lon) location.
875  REAL, INTENT(in) :: lat !< Latitude location
876  REAL, INTENT(in) :: lon !< Longitude location
877 
878  INTEGER :: indxi, indxj !< Indexes to be returned.
879  INTEGER :: indxi_tmp !< Hold the indxI value if on tile 3 or 4
880  INTEGER :: dimi, dimj !< Size of the grid dimensions
881  INTEGER :: i,j !< Count indexes
882  INTEGER :: jstart, jend, nextj !< j counting variables
883  TYPE(point) :: origpt !< Original point
884  TYPE(point), DIMENSION(4) :: points !< xyz of 8 nearest neighbors
885  REAL, DIMENSION(4) :: distsqrd !< distance between origPt and points(:)
886 
887  ! Set the inital fail values for indxI and indxJ
888  indxi = 0
889  indxj = 0
890 
891  dimi = diag_global_grid%adimI
892  dimj = diag_global_grid%adimJ
893 
894  ! check to see if the 'fix' for the latitude index is needed
895  IF ( diag_global_grid%aglo_lat(1,1) > &
896  &diag_global_grid%aglo_lat(1,2) ) THEN
897  ! reverse the j search
898  jstart = dimj-1
899  jend = 1
900  nextj = -1
901  ELSE
902  jstart = 0
903  jend = dimj-2
904  nextj = 1
905  END IF
906 
907  ! find the I index
908  iloop: DO i=0, dimi-2
909  IF ( diag_global_grid%aglo_lon(i,0) >&
910  & diag_global_grid%aglo_lon(i+1,0) ) THEN
911  ! We are at the 0 longitudal line
912  IF ( (diag_global_grid%aglo_lon(i,0) <= lon .AND. lon <= 360.) .OR.&
913  & (0. <= lon .AND. lon < diag_global_grid%aglo_lon(i+1, 0)) ) THEN
914  indxi = i
915  EXIT iloop
916  END IF
917  ELSEIF ( diag_global_grid%aglo_lon(i,0) <= lon .AND.&
918  & lon <= diag_global_grid%aglo_lon(i+1,0) ) THEN
919  indxi = i
920  EXIT iloop
921  END IF
922  END DO iloop
923 
924  ! Find the J index
925  IF ( indxi > 0 ) THEN
926  jloop: DO j=jstart, jend, nextj
927  IF ( diag_global_grid%aglo_lat(indxi,j) <= lat .AND.&
928  & lat <= diag_global_grid%aglo_lat(indxi,j+nextj) ) THEN
929  indxj = j
930  EXIT jloop
931  END IF
932  END DO jloop
933  END IF
934 
935  ! Make sure we have indexes in the correct range
936  valid: IF ( (indxi <= 0 .OR. dimi-1 < indxi) .OR. &
937  & (indxj <= 0 .OR. dimj-1 < indxj) ) THEN
938  indxi = 0
939  indxj = 0
940  ELSE ! indxI and indxJ are valid.
941  ! Since we are looking for the closest grid point to the
942  ! (lat,lon) point, we need to check the surrounding
943  ! points. The indexes for the variable points are as follows
944  !
945  ! 1---3
946  ! | |
947  ! 2---4
948 
949  ! The original point
950  origpt = latlon2xyz(lat,lon)
951 
952  ! Set the 'default' values for points(:) x,y,z to some large
953  ! value.
954  DO i=1, 4
955  points(i)%x = 1.0e20
956  points(i)%y = 1.0e20
957  points(i)%z = 1.0e20
958  END DO
959 
960  ! The original point
961  origpt = latlon2xyz(lat,lon)
962 
963  points(1) = latlon2xyz(diag_global_grid%aglo_lat(indxi,indxj),&
964  & diag_global_grid%aglo_lon(indxi,indxj))
965  points(2) = latlon2xyz(diag_global_grid%aglo_lat(indxi,indxj+nextj),&
966  & diag_global_grid%aglo_lon(indxi,indxj+nextj))
967  points(3) = latlon2xyz(diag_global_grid%aglo_lat(indxi+1,indxj+nextj),&
968  & diag_global_grid%aglo_lon(indxi+1,indxj+nextj))
969  points(4) = latlon2xyz(diag_global_grid%aglo_lat(indxi+1,indxj),&
970  & diag_global_grid%aglo_lon(indxi+1,indxj))
971 
972  ! Find the distance between the original point and the four
973  ! grid points
974  distsqrd = distancesqrd(origpt, points)
975 
976  SELECT CASE (minloc(distsqrd,1))
977  CASE ( 1 )
978  indxi = indxi;
979  indxj = indxj;
980  CASE ( 2 )
981  indxi = indxi;
982  indxj = indxj+nextj;
983  CASE ( 3 )
984  indxi = indxi+1;
985  indxj = indxj+nextj;
986  CASE ( 4 )
987  indxi = indxi+1;
988  indxj = indxj;
989  CASE DEFAULT
990  indxi = 0;
991  indxj = 0;
992  END SELECT
993 
994  ! If we are on tile 3 or 4, then the indxI and indxJ are
995  ! reversed due to the transposed grids.
996  IF ( diag_global_grid%tile_number == 4 .OR.&
997  & diag_global_grid%tile_number == 5 ) THEN
998  indxi_tmp = indxi
999  indxi = indxj
1000  indxj = indxi_tmp
1001  END IF
1002  END IF valid
1003 
1004  ! Set the return value for the function
1005  find_equator_index_agrid = (/indxi, indxj/)
1006  END FUNCTION find_equator_index_agrid
1007 
1008  !> @brief Return the (x,y,z) position of a given (lat,lon) point.
1009  !!
1010  !> Given a specific (lat, lon) point on the Earth, return the
1011  !! corresponding (x,y,z) location. The return of latlon2xyz
1012  !! will be either a scalar or an array of the same size as lat
1013  !! and lon.
1014  !! @return The return of latlon2xyz
1015  !! will be either a scalar or an array of the same size as lat
1016  !! and lon.
1017  PURE ELEMENTAL TYPE(point) function latlon2xyz(lat, lon)
1018  REAL, INTENT(in) :: lat !< The latitude of the (x,y,z) location to find. <TT>lat</TT>
1019  !! can be either a scalar or array. <TT>lat</TT> must be of the
1020  !! same rank / size as <TT>lon</TT>. This function assumes
1021  !! <TT>lat</TT> is in the range [-90,90].
1022  REAL, INTENT(in) :: lon !< The longitude of the (x,y,z) location to find. <TT>lon</TT>
1023  !! can be either a scalar or array. <TT>lon</TT> must be of the
1024  !! same rank / size as <TT>lat</TT>. This function assumes
1025  !! <TT>lon</TT> is in the range [0,360].
1026 
1027  ! lat/lon angles in radians
1028  REAL :: theta !< lat angles in radians
1029  REAL :: phi !< lon angles in radians
1030 
1031  ! Convert the lat lon values to radians The lat values passed in
1032  ! are in the range [-90,90], but we need to have a radian range
1033  ! [0,pi], where 0 is at the north pole. This is the reason for
1034  ! the subtraction from 90
1035  theta = deg2rad(90.-lat)
1036  phi = deg2rad(lon)
1037 
1038  ! Calculate the x,y,z point
1039  latlon2xyz%x = radius * sin(theta) * cos(phi)
1040  latlon2xyz%y = radius * sin(theta) * sin(phi)
1041  latlon2xyz%z = radius * cos(theta)
1042  END FUNCTION latlon2xyz
1043 
1044  !> @brief Find the distance between two points in the Cartesian
1045  !! coordinate space.
1046  !!
1047  !> <TT>distanceSqrd</TT> will find the distance squared between
1048  !! two points in the xyz coordinate space. <TT>pt1</TT> and <TT>
1049  !! pt2</TT> can either be both scalars, both arrays of the same
1050  !! size, or one a scalar and one an array. The return value
1051  !! will be a scalar or array of the same size as the input array.
1052  !! @return The return value will be a scalar or array of the same size as the input array.
1053  PURE ELEMENTAL REAL FUNCTION distancesqrd(pt1, pt2)
1054  TYPE(point), INTENT(in) :: pt1, pt2
1055 
1056  distancesqrd = (pt1%x-pt2%x)**2 +&
1057  & (pt1%y-pt2%y)**2 +&
1058  & (pt1%z-pt2%z)**2
1059  END FUNCTION distancesqrd
1060 
1061  !> @brief Find the distance, along the geodesic, between two points.
1062  !!
1063  !> <TT>gCirDistance</TT> will find the distance, along the geodesic, between two points defined
1064  !! by the (lat,lon) position of each point.
1065  !! @return real
1066  PURE ELEMENTAL REAL FUNCTION gcirdistance(lat1, lon1, lat2, lon2)
1067  REAL, INTENT(in) :: lat1, lat2, lon1, lon2
1068 
1069  REAL :: theta1, theta2
1070  REAL :: deltalambda !< Difference in longitude angles, in radians.
1071  REAL :: deltatheta !< Difference in latitude angels, in radians.
1072 
1073  theta1 = deg2rad(lat1)
1074  theta2 = deg2rad(lat2)
1075  deltalambda = deg2rad(lon2-lon1)
1076  deltatheta = deg2rad(lat2-lat1)
1077 
1078  gcirdistance = radius * 2. * asin(sqrt((sin(deltatheta/2.))**2 + cos(theta1)*cos(theta2)*(sin(deltalambda/2.))**2))
1079  END FUNCTION gcirdistance
1080 
1081  !> @brief Find the i,j indices and distance of the a-grid point nearest to
1082  !! the inputted lat,lon point.
1083  SUBROUTINE find_nearest_agrid_index(lat, &
1084  lon, &
1085  minI, &
1086  minJ, &
1087  minimum_distance)
1088 
1089  !Inputs/outputs
1090  REAL,INTENT(IN) :: lat
1091  REAL,INTENT(IN) :: lon
1092  INTEGER,INTENT(OUT) :: minI
1093  INTEGER,INTENT(OUT) :: minJ
1094  REAL,INTENT(OUT) :: minimum_distance
1095 
1096  !Local variables
1097  REAL :: llat
1098  REAL :: llon
1099  INTEGER :: j
1100  INTEGER :: i
1101  REAL :: dist
1102 
1103  !Since the poles have an non-unique longitude value, make a small
1104  !correction if looking for one of the poles.
1105  IF (lat .EQ. 90.0) THEN
1106  llat = lat - .1
1107  ELSEIF (lat .EQ. -90.0) THEN
1108  llat = lat + .1
1109  ELSE
1110  llat = lat
1111  END IF
1112  llon = lon
1113 
1114  !Loop through non-halo points. Calculate the distance
1115  !between each a-grid point and the point that we
1116  !are seeking. Store the minimum distance and its
1117  !corresponding i,j indices.
1118  mini = 0
1119  minj = 0
1120  minimum_distance = 2.0*radius*3.141592653
1121  DO j = 1,diag_global_grid%adimJ-2
1122  DO i = 1,diag_global_grid%adimI-2
1123  dist = gcirdistance(llat, &
1124  llon, &
1125  diag_global_grid%aglo_lat(i,j), &
1126  diag_global_grid%aglo_lon(i,j))
1127  IF (dist .LT. minimum_distance) THEN
1128 
1129  !These number shouldn't be hardcoded, but they have to
1130  !match the ones in diag_grid_init.
1131  if (diag_global_grid%tile_number .eq. 4 .or. &
1132  diag_global_grid%tile_number .eq. 5) then
1133 
1134  !Because of transpose in diag_grid_init.
1135  mini = j
1136  minj = i
1137 
1138  else
1139  mini = i
1140  minj = j
1141  endif
1142  minimum_distance = dist
1143  ENDIF
1144  ENDDO
1145  ENDDO
1146 
1147  !Check that valid i,j indices have been found.
1148  IF (mini .EQ. 0 .OR. minj .EQ. 0) THEN
1149  call error_mesg("find_nearest_agrid_index", &
1150  "A minimum distance was not found.", &
1151  fatal)
1152  ENDIF
1153 
1154  END SUBROUTINE find_nearest_agrid_index
1155 
1156 END MODULE diag_grid_mod
1157 !> @}
1158 ! close documentation grouping
pure elemental real function distancesqrd(pt1, pt2)
Find the distance between two points in the Cartesian coordinate space.
Definition: diag_grid.F90:1054
pure integer function, dimension(2) find_equator_index_agrid(lat, lon)
Return the closest index (i,j) to the given (lat,lon) point.
Definition: diag_grid.F90:873
subroutine, public diag_grid_init(domain, glo_lat, glo_lon, aglo_lat, aglo_lon)
Send the global grid to the diag_manager_mod for regional output.
Definition: diag_grid.F90:131
pure elemental type(point) function latlon2xyz(lat, lon)
Return the (x,y,z) position of a given (lat,lon) point.
Definition: diag_grid.F90:1018
pure elemental real function gcirdistance(lat1, lon1, lat2, lon2)
Find the distance, along the geodesic, between two points.
Definition: diag_grid.F90:1067
type(diag_global_grid_type) diag_global_grid
Variable to hold the global grid data.
Definition: diag_grid.F90:108
logical diag_grid_initialized
Indicates if the diag_grid_mod has been initialized.
Definition: diag_grid.F90:110
subroutine, public diag_grid_end()
Unallocate the diag_global_grid variable.
Definition: diag_grid.F90:349
pure elemental real function deg2rad(angle)
Convert an angle in degrees to radians.
Definition: diag_grid.F90:696
subroutine, public get_local_indexes2(lat, lon, iindex, jindex)
Find the indices of the nearest grid point of the a-grid to the specified (lon,lat) location on the l...
Definition: diag_grid.F90:640
subroutine, public get_local_indexes(latStart, latEnd, lonStart, lonEnd, istart, iend, jstart, jend)
Find the local start and local end indexes on the local PE for regional output.
Definition: diag_grid.F90:391
pure integer function, dimension(2) find_pole_index_agrid(lat, lon)
Return the closest index (i,j) to the given (lat,lon) point.
Definition: diag_grid.F90:709
subroutine find_nearest_agrid_index(lat, lon, minI, minJ, minimum_distance)
Find the i,j indices and distance of the a-grid point nearest to the inputted lat,...
Definition: diag_grid.F90:1088
pure elemental real function rad2deg(angle)
Definition: diag_grid.F90:683
Private type to hold the model's global grid data, and other grid information for use in this module.
Definition: diag_grid.F90:72
Private type to hold the corresponding (x,y,z) location for a (lat,lon) location.
Definition: diag_grid.F90:99
subroutine, public write_version_number(version, tag, unit)
Prints to the log file (or a specified unit) the version id string and tag name.
Definition: fms.F90:701
subroutine, public error_mesg(routine, message, level)
Print notes, warnings and error messages; terminates program for warning and error messages....
Definition: fms.F90:441
integer function, dimension(size(domain%tile_id(:))) mpp_get_tile_id(domain)
Returns the tile_id on current pe.
integer function mpp_get_ntile_count(domain)
Returns number of tiles in mosaic.
Retrieve the entire array of compute domain extents associated with a decomposition.
The domain2D type contains all the necessary information to define the global, compute and data domai...
integer function mpp_npes()
Returns processor count for current pelist.
Definition: mpp_util.inc:420
integer function mpp_pe()
Returns processor ID.
Definition: mpp_util.inc:406
Reduction operations. Find the max of scalar a from the PEs in pelist result is also automatically br...
Definition: mpp.F90:537
Reduction operations. Find the min of scalar a from the PEs in pelist result is also automatically br...
Definition: mpp.F90:559