1module getGridPosition_mod
2 ! MODULE getGridPosition_mod (prefix='gpos' category='8. Low-level utilities and constants')
3 !
4 !:Purpose: A place to collect numerous interpolation related routines.
5 ! The main task of the module is to compute the grid XY position from a lat-lon.
6 ! This simply calls the ezsint routine gdxyfll for simple grids. For
7 ! Yin-Yan grids it calls the function gpos_xyfll_yinYangGrid
8 ! (see below in this module). There is also support for
9 ! RPN Y grids, in which case it calls the subroutine gpos_xyfll_unstructGrid.
10 !
11 use kdTree2_mod
12 use mathPhysConstants_mod
13 use physicsFunctions_mod
14 use utilities_mod
15
16 implicit none
17 save
18 private
19
20 ! public procedures
21 public :: gpos_getPositionXY, gpos_gridIsOrca
22
23 integer, parameter :: maxNumLocalGridPointsSearch = 350
24 integer, external :: get_max_rss
25
26 type(kdtree2), pointer :: tree => null()
27
28contains
29
30 !---------------------------------------------------------
31 ! gpos_getPositionXY
32 !---------------------------------------------------------
33 function gpos_getPositionXY( gdid, xpos_r4, ypos_r4, xpos2_r4, ypos2_r4, &
34 lat_deg_r4, lon_deg_r4, subGridIndex ) result(ierr)
35 !
36 !:Purpose: Compute the grid XY position from a lat-lon. This
37 ! simply calls the ezsint routine gdxyfll for simple grids. For
38 ! Yin-Yan grids it calls the function gpos_xyfll_yinYangGrid
39 ! (see below in this module). There is also support for
40 ! RPN Y grids, in which case it calls the subroutine gpos_xyfll_unstructGrid.
41 !
42 implicit none
43
44 ! Arguments:
45 integer, intent(in) :: gdid
46 integer, intent(out) :: subGridIndex
47 real(4), intent(out) :: xpos_r4
48 real(4), intent(out) :: ypos_r4
49 real(4), intent(out) :: xpos2_r4
50 real(4), intent(out) :: ypos2_r4
51 real(4), intent(in) :: lat_deg_r4
52 real(4), intent(in) :: lon_deg_r4
53 ! Result:
54 integer :: ierr ! returned value of function
55
56 ! Locals:
57 integer :: numSubGrids
58 integer :: ezget_nsubGrids, gdxyfll, ezgprm
59 character(len=1) :: grtyp
60 integer :: ni, nj, ig1, ig2, ig3, ig4
61
62 numSubGrids = ezget_nsubGrids(gdid)
63 xpos2_r4 = -999.0
64 ypos2_r4 = -999.0
65
66 if ( numSubGrids == 1 ) then
67
68 ierr = ezgprm(gdid, grtyp, ni, nj, ig1, ig2, ig3, ig4)
69
70 if (grtyp == 'Y') then
71
72 !$omp critical
73 ierr = gpos_xyfll_unstructGrid(gdid, xpos_r4, ypos_r4, lat_deg_r4, lon_deg_r4)
74 !$omp end critical
75
76 else
77
78 ! Not an unstructured grid, nor a Yin-Yang grid, call the standard ezscint routine
79 ierr = gdxyfll(gdid, xpos_r4, ypos_r4, lat_deg_r4, lon_deg_r4, 1)
80
81 end if
82
83 subGridIndex = 1
84
85 else
86
87 ! This is a Yin-Yang grid, do something different
88 !$omp critical
89 ierr = gpos_xyfll_yinYangGrid(gdid, xpos_r4, ypos_r4, xpos2_r4, ypos2_r4, &
90 lat_deg_r4, lon_deg_r4, subGridIndex)
91 !$omp end critical
92
93 end if
94
95 if ( subGridIndex /= 3 ) then
96 ! when only returning 1 position, copy values to pos2
97 xpos2_r4 = xpos_r4
98 ypos2_r4 = ypos_r4
99 end if
100
101 end function gpos_getPositionXY
102
103 !---------------------------------------------------------
104 ! gpos_xyfll_yinYangGrid
105 !---------------------------------------------------------
106 function gpos_xyfll_yinYangGrid( gdid, xpos_r4, ypos_r4, xpos2_r4, ypos2_r4, &
107 lat_deg_r4, lon_deg_r4, subGridIndex ) result(ierr)
108 !
109 !:Purpose: Compute the grid XY position from a lat-lon for a Yin-Yang grid.
110 ! It returns locations from both the Yin and Yang
111 ! subgrids when in the overlap region, depending on the logical
112 ! variable `useSingleValueOverlap`.
113 !
114 implicit none
115
116 ! Arguments:
117 integer, intent(in) :: gdid
118 integer, intent(out) :: subGridIndex
119 real(4), intent(in) :: lat_deg_r4
120 real(4), intent(in) :: lon_deg_r4
121 real(4), intent(out) :: xpos_r4
122 real(4), intent(out) :: ypos_r4
123 real(4), intent(out) :: xpos2_r4
124 real(4), intent(out) :: ypos2_r4
125 ! Result:
126 integer :: ierr ! returned value of function
127
128 ! Locals:
129 integer :: ezget_subGridids, gdgaxes, gdxyfll, ezgprm
130 integer :: EZscintIDvec(2)
131 integer, save :: EZscintIDvec1_old = -999
132 integer :: lonIndex, latIndex
133 real :: lonrot, latrot
134 real, allocatable, save :: ax_yin(:), ay_yin(:), ax_yan(:), ay_yan(:)
135 logical :: axesDifferent
136 character(len=1) :: grtyp
137 integer :: ni, nj, ig1, ig2, ig3, ig4
138 ! this controls which approach to use for interpolation within the YIN-YAN overlap
139 logical :: useSingleValueOverlap = .true.
140
141 ierr = ezget_subGridids(gdid, EZscintIDvec)
142 ! get ni nj of subGrid, assume same for both YIN and YANG
143 ierr = ezgprm(EZscintIDvec(1), grtyp, ni, nj, ig1, ig2, ig3, ig4)
144
145 ! first check YIN
146 ierr = gdxyfll(EZscintIDvec(1), xpos_r4, ypos_r4, lat_deg_r4, lon_deg_r4, 1)
147
148 ! compute rotated lon and lat at obs location
149 axesDifferent = (EZscintIDvec1_old /= EZscintIDvec(1))
150 if (axesDifferent) then
151 write(*,*) 'gpos_getPositionXY: axesDifferent, compute needed parameters'
152 if (allocated(ax_yin)) deallocate(ax_yin,ay_yin)
153 allocate(ax_yin(ni),ay_yin(nj))
154 ierr = gdgaxes(EZscintIDvec(1), ax_yin, ay_yin)
155 EZscintIDvec1_old = EZscintIDvec(1)
156 end if
157 lonIndex = floor(xpos_r4)
158 if ( lonIndex >= 1 .and. (lonIndex+1) <= ni ) then
159 lonrot = ax_yin(lonIndex) + (ax_yin(lonIndex+1) - ax_yin(lonIndex)) * &
160 (xpos_r4 - lonIndex)
161 else
162 lonrot = -999.0
163 end if
164 latIndex = floor(ypos_r4)
165 if ( latIndex >= 1 .and. (latIndex+1) <= nj ) then
166 latrot = ay_yin(latIndex) + (ay_yin(latIndex+1) - ay_yin(latIndex)) * &
167 (ypos_r4 - latIndex)
168 else
169 latrot = -999.0
170 end if
171 subGridIndex = 1
172
173 if ( useSingleValueOverlap ) then
174
175 ! this approach is most similar to how ezsint works, preferentially take YIN
176
177 if ( lonrot < 45.0 .or. lonrot > 315.0 .or. latrot < -45.0 .or. latrot > 45.0 ) then
178 ! Outside YIN, therefore use YANG (assume it is inside YANG)
179 ierr = gdxyfll(EZscintIDvec(2), xpos_r4, ypos_r4, lat_deg_r4, lon_deg_r4, 1)
180 ypos_r4 = ypos_r4 + real(nj) ! shift from YANG position to Supergrid position
181 subGridIndex = 2
182 else
183 subGridIndex = 1
184 end if
185
186 else ! not useSingleValueOverlap
187
188 ! this approach returns both the YIN and YAN locations when point is inside both
189
190 if ( lonrot < 45.0 .or. lonrot > 315.0 .or. latrot < -45.0 .or. latrot > 45.0 ) then
191 ! Outside YIN, therefore use YANG (assume it is inside YANG)
192 ierr = gdxyfll(EZscintIDvec(2), xpos_r4, ypos_r4, lat_deg_r4, lon_deg_r4, 1)
193 ypos_r4 = ypos_r4 + real(nj) ! shift from YANG position to Supergrid position
194 subGridIndex = 2
195 else
196 ! inside YIN, check if also inside YANG
197 allocate(ax_yan(ni),ay_yan(nj))
198 ierr = gdgaxes(EZscintIDvec(2), ax_yan, ay_yan)
199 ierr = gdxyfll(EZscintIDvec(2), xpos2_r4, ypos2_r4, lat_deg_r4, lon_deg_r4, 1)
200 if ( lonIndex >= 1 .and. (lonIndex+1) <= ni ) then
201 lonrot = ax_yan(lonIndex) + (ax_yan(lonIndex+1) - ax_yan(lonIndex)) * &
202 (xpos2_r4 - lonIndex)
203 else
204 lonrot = -999.0
205 end if
206 latIndex = floor(ypos2_r4)
207 if ( latIndex >= 1 .and. (latIndex+1) <= nj ) then
208 latrot = ay_yan(latIndex) + (ay_yan(latIndex+1) - ay_yan(latIndex)) * &
209 (ypos2_r4 - latIndex)
210 else
211 latrot = -999.0
212 end if
213 deallocate(ax_yan,ay_yan)
214 if ( lonrot < 45.0 .or. lonrot > 315.0 .or. latrot < -45.0 .or. latrot > 45.0 ) then
215 ! outside YANG, only inside YIN
216 xpos2_r4 = -999.0
217 ypos2_r4 = -999.0
218 subGridIndex = 1
219 else
220 ! inside both YIN and YANG
221 ypos2_r4 = ypos2_r4 + real(nj) ! shift from YANG position to Supergrid position
222 subGridIndex = 3
223 end if
224 end if
225
226 end if
227
228 end function gpos_xyfll_yinYangGrid
229
230 !---------------------------------------------------------
231 ! gpos_xyfll_unstructGrid
232 !---------------------------------------------------------
233 function gpos_xyfll_unstructGrid( gdid, xpos_r4, ypos_r4, lat_deg_r4, lon_deg_r4 ) result(ierr)
234
235 !:Purpose: This function is used to interpolate from a RPN grid to a location
236 ! specified by a latitude and longitude.
237 ! It has special treatment for the ORCA12 and ORCA025 tri-polar grids
238 ! used for sea ice and ocean models.
239 ! The kdtree2 module is used to efficiently
240 ! perform this task. The kdtree itself is constructed on the first call.
241 implicit none
242
243 ! Arguments:
244 integer, intent(in) :: gdid
245 real(4), intent(in) :: lat_deg_r4
246 real(4), intent(in) :: lon_deg_r4
247 real(4), intent(out) :: xpos_r4
248 real(4), intent(out) :: ypos_r4
249 ! Result:
250 integer :: ierr ! returned value of function
251
252 ! Locals:
253 real(8) :: lon_rad_r8, lat_rad_r8, pertGridLonRad, pertGridLatRad
254 real(8), parameter :: deltaGrid = 0.00001d0
255 real(8) :: pertPosition(3)
256 real(8) :: gridSpacing
257 real(8) :: gridSpacingSquared, lowerLeftCornerDistSquared, lowerRightCornerDistSquared, upperLeftCornerDistSquared
258 integer, save :: gdidOld = -999
259 integer :: nx, ny
260 integer, save :: ni, nj
261 real(8), allocatable, save :: grid_lon_rad(:,:), grid_lat_rad(:,:)
262 real(4), save :: maxGridSpacing
263 integer, save :: startXIndex, startYIndex, endXIndex, endYIndex
264 character(len=1) :: grtyp
265 integer :: ig1, ig2, ig3, ig4
266 real(4), allocatable :: grid_lat_deg_r4(:,:), grid_lon_deg_r4(:,:)
267 integer :: closePointsIndex
268 integer :: gridIndex
269 integer :: xIndex, yIndex, xIndexMin, xIndexMax, yIndexMin, yIndexMax
270 integer :: ezgprm, gdll ! Functions
271 integer :: numLocalGridPointsFound
272 real(kdkind), allocatable :: positionArray(:,:)
273 type(kdtree2_result) :: searchResults(maxNumLocalGridPointsSearch)
274 real(kdkind) :: maxRadiusSquared
275 real(kdkind) :: refPosition(3)
276
277 if ( gdid /= gdidOld .and. gdidOld > 0) then
278 write(*,*) 'gpos_xyfll_unstructGrid: gdid gdidOld = ',gdid, gdidOld
279 call utl_abort('gpos_xyfll_unstructGrid: only one grid expected. Change code !')
280 end if
281
282 ! create the kdtree on the first call
283 if (.not. associated(tree)) then
284 write(*,*) 'gpos_xyfll_unstructGrid: start creating kdtree'
285 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
286 ierr = ezgprm(gdid, grtyp, ni, nj, ig1, ig2, ig3, ig4)
287
288 startXIndex = 1
289 startYIndex = 1
290 endXIndex = ni
291 endYIndex = nj
292 maxGridSpacing = 10000.0
293
294 if (allocated(grid_lat_rad)) deallocate(grid_lat_rad, grid_lon_rad)
295 allocate(grid_lat_rad(ni,nj), grid_lon_rad(ni,nj))
296 allocate(grid_lat_deg_r4(ni,nj), grid_lon_deg_r4(ni,nj))
297 ierr = gdll(gdid, grid_lat_deg_r4, grid_lon_deg_r4)
298 grid_lat_rad(:,:) = real(grid_lat_deg_r4(:,:),8)*MPC_RADIANS_PER_DEGREE_R8
299 grid_lon_rad(:,:) = real(grid_lon_deg_r4(:,:),8)*MPC_RADIANS_PER_DEGREE_R8
300 deallocate(grid_lat_deg_r4, grid_lon_deg_r4)
301 gdidOld = gdid
302
303 if( gpos_gridIsOrca(ni, nj, grid_lat_rad, grid_lon_rad) ) then
304 ! The last 2 columns (i=ni-1 and ni) are repetitions of the first 2 columns (i=1 and 2)
305 ! so both ends are removed from the search area.
306 startXIndex = 2
307 endXIndex = ni - 1
308 ! The last line (j=nj) is a repetition in reverse order of the line (j=nj-2)
309 ! and is thus eliminated from the search domain.
310 endYIndex = nj - 1
311 end if
312 if(ni == 1442 .and. nj == 1021) then
313 ! This is orca025
314 maxGridSpacing = 40000.0
315 end if
316
317 nx = endXIndex - startXIndex + 1
318 ny = endYIndex - startYIndex + 1
319
320 allocate(positionArray(3,nx*ny))
321
322 gridIndex = 0
323 do yIndex = startYIndex, endYIndex
324 do xIndex = startXIndex, endXIndex
325 gridIndex = gridIndex + 1
326 positionArray(:,gridIndex) = kdtree2_3dPosition(grid_lon_rad(xIndex,yIndex), grid_lat_rad(xIndex,yIndex))
327 end do
328 end do
329 tree => kdtree2_create(positionArray, sort=.true., rearrange=.true.)
330 write(*,*) 'gpos_xyfll_unstructGrid: done creating kdtree'
331 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
332
333 else
334
335 ierr = 0
336 nx = endXIndex - startXIndex + 1
337
338 end if
339
340 ! do the search
341 maxRadiusSquared = maxGridSpacing**2
342 lon_rad_r8 = real(lon_deg_r4,8)*MPC_RADIANS_PER_DEGREE_R8
343 lat_rad_r8 = real(lat_deg_r4,8)*MPC_RADIANS_PER_DEGREE_R8
344 refPosition(:) = kdtree2_3dPosition(lon_rad_r8, lat_rad_r8)
345
346 call kdtree2_r_nearest(tp=tree, qv=refPosition, r2=maxRadiusSquared, nfound=numLocalGridPointsFound,&
347 nalloc=maxNumLocalGridPointsSearch, results=searchResults)
348 if (numLocalGridPointsFound > maxNumLocalGridPointsSearch) then
349 call utl_abort('gpos_xyfll_unstructGrid: the parameter maxNumLocalGridPointsSearch must be increased')
350 end if
351
352 if (numLocalGridPointsFound < 4) then
353 write(*,*) 'gpos_xyfll_unstructGrid: numLocalGridPointsFound = ',numLocalGridPointsFound
354 write(*,*) 'gpos_xyfll_unstructGrid: search radius = ',maxGridSpacing,' meters'
355 write(*,*) 'gpos_xyfll_unstructGrid: obs lon, lat = ', lon_deg_r4, lat_deg_r4
356 call utl_abort('gpos_xyfll_unstructGrid: the search did not find 4 close points.')
357 end if
358
359 if (searchResults(1)%dis > maxRadiusSquared) then
360 write(*,*) 'gpos_xyfll_unstructGrid: No grid point found within ',maxGridSpacing,' meters'
361 write(*,*) 'of the reference location lat-lon (degrees): ',lat_deg_r4,lon_deg_r4
362 xpos_r4 = -999.0
363 ypos_r4 = -999.0
364 return
365 end if
366
367 ! We found the closest grid point to the reference point.
368 ! Now we need to determine in which of the 4 quadrants around the grid point the reference point lies.
369
370 closePointsIndex = 1
371
372 gridIndex = searchResults(closePointsIndex)%idx
373 yIndex = (gridIndex-1)/nx + startYIndex
374 xIndex = gridIndex - (yIndex-startYIndex)*nx + startXIndex - 1
375
376 ! Perturbing the closest grid point in the x direction
377
378 if ( xIndex < ni ) then
379
380 if ( grid_lon_rad(xIndex+1,yIndex) - grid_lon_rad(xIndex,yIndex) > 5.0d0 ) then
381 pertGridLonRad = grid_lon_rad(xIndex,yIndex) + (grid_lon_rad(xIndex+1,yIndex) - (grid_lon_rad(xIndex,yIndex) + 2.d0 * MPC_PI_R8))*deltaGrid
382 else if ( grid_lon_rad(xIndex+1,yIndex) - grid_lon_rad(xIndex,yIndex) < -5.0d0 ) then
383 pertGridLonRad = grid_lon_rad(xIndex,yIndex) + ((grid_lon_rad(xIndex+1,yIndex) + 2.d0 * MPC_PI_R8) - grid_lon_rad(xIndex,yIndex))*deltaGrid
384 else
385 pertGridLonRad = grid_lon_rad(xIndex,yIndex) + (grid_lon_rad(xIndex+1,yIndex) - grid_lon_rad(xIndex,yIndex))*deltaGrid
386 end if
387 pertGridLatRad = grid_lat_rad(xIndex,yIndex) + (grid_lat_rad(xIndex+1,yIndex) - grid_lat_rad(xIndex,yIndex))*deltaGrid
388
389 ! Calculate the distance between the reference point and the perturbed grid point.
390 pertPosition(:) = kdtree2_3dPosition(pertGridLonRad, pertGridLatRad)
391
392 gridSpacingSquared = sum( (pertPosition(:) - refPosition(:))**2 )
393
394 if (gridSpacingSquared < searchResults(closePointsIndex)%dis) then
395 xIndexMin = xIndex
396 xIndexMax = xIndex + 1
397 else
398 if ( xIndex > 1 ) then
399 xIndexMin = xIndex - 1
400 xIndexMax = xIndex
401 else
402 write(*,*) 'xIndex yIndex ', xIndex, yIndex
403 write(*,*) 'Closest distance = ',sqrt(searchResults(closePointsIndex)%dis)
404 write(*,*) 'Perturbed distance = ',sqrt(gridSpacingSquared)
405 call utl_abort('gpos_xyfll_unstructGrid: 1. the reference point is outside the grid.')
406 end if
407 end if
408
409 else
410
411 if ( grid_lon_rad(xIndex-1,yIndex) - grid_lon_rad(xIndex,yIndex) > 5.0d0 ) then
412 pertGridLonRad = grid_lon_rad(xIndex,yIndex) + (grid_lon_rad(xIndex-1,yIndex) - (grid_lon_rad(xIndex,yIndex) + 2.d0 * MPC_PI_R8))*deltaGrid
413 else if ( grid_lon_rad(xIndex-1,yIndex) - grid_lon_rad(xIndex,yIndex) < -5.0d0 ) then
414 pertGridLonRad = grid_lon_rad(xIndex,yIndex) + ((grid_lon_rad(xIndex-1,yIndex) + 2.d0 * MPC_PI_R8) - grid_lon_rad(xIndex,yIndex))*deltaGrid
415 else
416 pertGridLonRad = grid_lon_rad(xIndex,yIndex) + (grid_lon_rad(xIndex-1,yIndex) - grid_lon_rad(xIndex,yIndex))*deltaGrid
417 end if
418 pertGridLatRad = grid_lat_rad(xIndex,yIndex) + (grid_lat_rad(xIndex-1,yIndex) - grid_lat_rad(xIndex,yIndex))*deltaGrid
419
420 ! Calculate the distance between the reference point and the perturbed grid point.
421 pertPosition(:) = kdtree2_3dPosition(pertGridLonRad, pertGridLatRad)
422
423 gridSpacingSquared = sum( (pertPosition(:) - refPosition(:))**2 )
424
425 if (gridSpacingSquared < searchResults(closePointsIndex)%dis) then
426 xIndexMin = xIndex - 1
427 xIndexMax = xIndex
428 else
429 write(*,*) 'xIndex yIndex ', xIndex, yIndex
430 write(*,*) 'Closest distance = ',sqrt(searchResults(closePointsIndex)%dis)
431 write(*,*) 'Perturbed distance = ',sqrt(gridSpacingSquared)
432 call utl_abort('gpos_xyfll_unstructGrid: 2. the reference point is outside the grid.')
433 end if
434
435 end if
436
437 ! Perturbing the closest grid point in the y direction
438
439 if ( yIndex < nj ) then
440
441 if ( grid_lon_rad(xIndex,yIndex+1) - grid_lon_rad(xIndex,yIndex) > 5.0d0 ) then
442 pertGridLonRad = grid_lon_rad(xIndex,yIndex) + (grid_lon_rad(xIndex,yIndex+1) - (grid_lon_rad(xIndex,yIndex) + 2.d0 * MPC_PI_R8))*deltaGrid
443 else if ( grid_lon_rad(xIndex,yIndex+1) - grid_lon_rad(xIndex,yIndex) < -5.0d0 ) then
444 pertGridLonRad = grid_lon_rad(xIndex,yIndex) + ((grid_lon_rad(xIndex,yIndex+1) + 2.d0 * MPC_PI_R8) - grid_lon_rad(xIndex,yIndex))*deltaGrid
445 else
446 pertGridLonRad = grid_lon_rad(xIndex,yIndex) + (grid_lon_rad(xIndex,yIndex+1) - grid_lon_rad(xIndex,yIndex))*deltaGrid
447 end if
448 pertGridLatRad = grid_lat_rad(xIndex,yIndex) + (grid_lat_rad(xIndex,yIndex+1) - grid_lat_rad(xIndex,yIndex))*deltaGrid
449
450 ! Calculate the distance between the reference point and the perturbed grid point.
451 pertPosition(:) = kdtree2_3dPosition(pertGridLonRad, pertGridLatRad)
452
453 gridSpacingSquared = sum( (pertPosition(:) - refPosition(:))**2 )
454
455 if (gridSpacingSquared < searchResults(closePointsIndex)%dis) then
456 yIndexMin = yIndex
457 yIndexMax = yIndex + 1
458 else
459 if ( yIndex > 1 ) then
460 yIndexMin = yIndex - 1
461 yIndexMax = yIndex
462 else
463 write(*,*) 'xIndex yIndex ', xIndex, yIndex
464 write(*,*) 'Closest distance = ',sqrt(searchResults(closePointsIndex)%dis)
465 write(*,*) 'Perturbed distance = ',sqrt(gridSpacingSquared)
466 call utl_abort('gpos_xyfll_unstructGrid: 3. the reference point is outside the grid.')
467 end if
468 end if
469
470 else
471
472 if ( grid_lon_rad(xIndex,yIndex-1) - grid_lon_rad(xIndex,yIndex) > 5.0d0 ) then
473 pertGridLonRad = grid_lon_rad(xIndex,yIndex) + (grid_lon_rad(xIndex,yIndex-1) - (grid_lon_rad(xIndex,yIndex) + 2.d0 * MPC_PI_R8))*deltaGrid
474 else if ( grid_lon_rad(xIndex,yIndex-1) - grid_lon_rad(xIndex,yIndex) < -5.0d0 ) then
475 pertGridLonRad = grid_lon_rad(xIndex,yIndex) + ((grid_lon_rad(xIndex,yIndex-1) + 2.d0 * MPC_PI_R8) - grid_lon_rad(xIndex,yIndex))*deltaGrid
476 else
477 pertGridLonRad = grid_lon_rad(xIndex,yIndex) + (grid_lon_rad(xIndex,yIndex-1) - grid_lon_rad(xIndex,yIndex))*deltaGrid
478 end if
479 pertGridLatRad = grid_lat_rad(xIndex,yIndex) + (grid_lat_rad(xIndex,yIndex-1) - grid_lat_rad(xIndex,yIndex))*deltaGrid
480
481 ! Calculate the distance between the reference point and the perturbed grid point.
482 pertPosition(:) = kdtree2_3dPosition(pertGridLonRad, pertGridLatRad)
483
484 gridSpacingSquared = sum( (pertPosition(:) - refPosition(:))**2 )
485
486 if (gridSpacingSquared < searchResults(closePointsIndex)%dis) then
487 yIndexMin = yIndex - 1
488 yIndexMax = yIndex
489 else
490 write(*,*) 'xIndex yIndex ', xIndex, yIndex
491 write(*,*) 'Closest distance = ',sqrt(searchResults(closePointsIndex)%dis)
492 write(*,*) 'Perturbed distance = ',sqrt(gridSpacingSquared)
493 call utl_abort('gpos_xyfll_unstructGrid: 4. the reference point is outside the grid.')
494 end if
495
496 end if
497
498 ! Calculate real x and y position in the grid.
499 ! kdtree returns distance in square meters !
500
501 gridSpacingSquared = phf_calcDistance(grid_lat_rad(xIndexMin,yIndexMin), grid_lon_rad(xIndexMin,yIndexMin), &
502 grid_lat_rad(xIndexMax,yIndexMin), grid_lon_rad(xIndexMax,yIndexMin))**2
503
504 lowerLeftCornerDistSquared = phf_calcDistance(grid_lat_rad(xIndexMin,yIndexMin), grid_lon_rad(xIndexMin,yIndexMin), &
505 lat_rad_r8, lon_rad_r8)**2
506
507 lowerRightCornerDistSquared = phf_calcDistance(grid_lat_rad(xIndexMax,yIndexMin), grid_lon_rad(xIndexMax,yIndexMin), &
508 lat_rad_r8, lon_rad_r8)**2
509
510 xpos_r4 = real(xIndexMin) + (lowerLeftCornerDistSquared + gridSpacingSquared - lowerRightCornerDistSquared)/(2.0*(gridSpacingSquared))
511
512 gridSpacingSquared = phf_calcDistance(grid_lat_rad(xIndexMin,yIndexMin), grid_lon_rad(xIndexMin,yIndexMin), &
513 grid_lat_rad(xIndexMin,yIndexMax), grid_lon_rad(xIndexMin,yIndexMax))**2
514
515 upperLeftCornerDistSquared = phf_calcDistance(grid_lat_rad(xIndexMin,yIndexMax), grid_lon_rad(xIndexMin,yIndexMax), &
516 lat_rad_r8, lon_rad_r8)**2
517
518 ypos_r4 = real(yIndexMin) + (lowerLeftCornerDistSquared + gridSpacingSquared - upperLeftCornerDistSquared)/(2.0*(gridSpacingSquared))
519
520 if ( abs(ypos_r4 - yIndexMin) > 2.0 .or. abs(xpos_r4 - xIndexMin) > 4.3 ) then
521 write(*,*) 'xpos_r4 = ',xpos_r4
522 write(*,*) 'ypos_r4 = ',ypos_r4
523 write(*,*) 'xIndexMin = ',xIndexMin
524 write(*,*) 'yIndexMin = ',yIndexMin
525 do closePointsIndex = 1,min(numLocalGridPointsFound,5)
526 gridIndex = searchResults(closePointsIndex)%idx
527 yIndex = (gridIndex-1)/nx + startYIndex
528 xIndex = gridIndex - (yIndex-startYIndex)*nx
529 write(*,*) 'closePointsIndex: ', closePointsIndex
530 write(*,*) 'xIndex yIndex ', xIndex, yIndex
531 write(*,*) 'idx: ',searchResults(closePointsIndex)%idx
532 write(*,*) 'distance (meters): ',sqrt(searchResults(closePointsIndex)%dis)
533 end do
534 write(*,*) 'lowerLeftCornerDist = ',sqrt(lowerLeftCornerDistSquared)
535 write(*,*) 'lowerRightCornerDist = ',sqrt(lowerRightCornerDistSquared)
536 write(*,*) 'upperLeftCornerDist = ',sqrt(upperLeftCornerDistSquared)
537 write(*,*) 'y gridSpacing = ',sqrt(gridSpacingSquared)
538 gridSpacing = phf_calcDistance(grid_lat_rad(xIndexMin,yIndexMin), grid_lon_rad(xIndexMin,yIndexMin), &
539 grid_lat_rad(xIndexMax,yIndexMin), grid_lon_rad(xIndexMax,yIndexMin))
540 write(*,*) 'x gridSpacing = ',gridSpacing
541 write(*,*) 'lon-lat (deg):'
542 write(*,*) 'reference point: ',lon_deg_r4, lat_deg_r4
543 write(*,*) xIndexMin,yIndexMin,grid_lon_rad(xIndexMin,yIndexMin)* MPC_DEGREES_PER_RADIAN_R8,grid_lat_rad(xIndexMin,yIndexMin)* MPC_DEGREES_PER_RADIAN_R8
544 write(*,*) xIndexMax,yIndexMin,grid_lon_rad(xIndexMax,yIndexMin)* MPC_DEGREES_PER_RADIAN_R8,grid_lat_rad(xIndexMax,yIndexMin)* MPC_DEGREES_PER_RADIAN_R8
545 write(*,*) xIndexMax+1,yIndexMin,grid_lon_rad(xIndexMax+1,yIndexMin)* MPC_DEGREES_PER_RADIAN_R8,grid_lat_rad(xIndexMax+1,yIndexMin)* MPC_DEGREES_PER_RADIAN_R8
546 write(*,*) xIndexMin,yIndexMax,grid_lon_rad(xIndexMin,yIndexMax)* MPC_DEGREES_PER_RADIAN_R8,grid_lat_rad(xIndexMin,yIndexMax)* MPC_DEGREES_PER_RADIAN_R8
547 write(*,*) xIndexMax,yIndexMax,grid_lon_rad(xIndexMax,yIndexMax)* MPC_DEGREES_PER_RADIAN_R8,grid_lat_rad(xIndexMax,yIndexMax)* MPC_DEGREES_PER_RADIAN_R8
548 write(*,*) xIndexMax+1,yIndexMax,grid_lon_rad(xIndexMax+1,yIndexMax)* MPC_DEGREES_PER_RADIAN_R8,grid_lat_rad(xIndexMax+1,yIndexMax)* MPC_DEGREES_PER_RADIAN_R8
549 write(*,*) xIndexMin,yIndexMax+1,grid_lon_rad(xIndexMin,yIndexMax+1)* MPC_DEGREES_PER_RADIAN_R8,grid_lat_rad(xIndexMin,yIndexMax+1)* MPC_DEGREES_PER_RADIAN_R8
550 write(*,*) xIndexMax,yIndexMax+1,grid_lon_rad(xIndexMax,yIndexMax+1)* MPC_DEGREES_PER_RADIAN_R8,grid_lat_rad(xIndexMax,yIndexMax+1)* MPC_DEGREES_PER_RADIAN_R8
551 write(*,*) xIndexMax+1,yIndexMax+1,grid_lon_rad(xIndexMax+1,yIndexMax+1)* MPC_DEGREES_PER_RADIAN_R8,grid_lat_rad(xIndexMax+1,yIndexMax+1)* MPC_DEGREES_PER_RADIAN_R8
552 call utl_abort('gpos_xyfll_unstructGrid: Check code !')
553 end if
554
555 end function gpos_xyfll_unstructGrid
556
557 !---------------------------------------------------------
558 ! gpos_gridIsOrca
559 !---------------------------------------------------------
560 function gpos_gridIsOrca( ni, nj, latitude, longitude ) result(gridIsOrca)
561
562 !:Purpose: Check if the grid is of the orca tri-polar family
563 implicit none
564
565 ! Arguments:
566 integer, intent(in) :: ni ! first dimension of the grid
567 integer, intent(in) :: nj ! second dimension of the grid
568 real(8), intent(in) :: longitude(ni,nj) ! longitude (degrees or radians)
569 real(8), intent(in) :: latitude(ni,nj) ! latitude (degrees or radians)
570 ! Result:
571 logical :: gridIsOrca ! returned value of function
572
573 ! Locals:
574 integer :: xIndex, yIndex
575
576 gridIsOrca = .true.
577
578 ! The last 2 columns (ni-1 and ni) are repetitions of the first 2 columns (i=1 and 2)
579 ! Check that is true
580 do yIndex = 1, nj
581 if(latitude(1,yIndex) /= latitude(ni-1,yIndex) .or. &
582 latitude(2,yIndex) /= latitude(ni, yIndex) .or. &
583 longitude(1,yIndex) /= longitude(ni-1,yIndex) .or. &
584 longitude(2,yIndex) /= longitude(ni, yIndex)) then
585 gridIsOrca = .false.
586 write(*,*) '1. Not an orca grid',latitude(1,yIndex),latitude(ni-1,yIndex),latitude(2,yIndex),latitude(ni, yIndex),longitude(1,yIndex),longitude(ni-1,yIndex),longitude(2,yIndex),longitude(ni, yIndex)
587 return
588 end if
589 end do
590
591 ! The last line (j=nj) is a repetition in reverse order of the line (j=nj-2)
592 ! Check that is true (does not work at (ni/2 + 1) for orca025 lats as with the letkf ocean unit test)
593 do xIndex = 2, ni
594 if(longitude(xIndex,nj) /= longitude(ni-xIndex+2,nj-2)) then
595 gridIsOrca = .false.
596 write(*,*) '2. Not an orca grid',xIndex,longitude(xIndex,nj),longitude(ni-xIndex+2,nj-2)
597 return
598 end if
599 end do
600 do xIndex = 2, ni/2
601 if(latitude(xIndex,nj) /= latitude(ni-xIndex+2,nj-2)) then
602 gridIsOrca = .false.
603 write(*,*) '3. Not an orca grid',xIndex,latitude(xIndex,nj),latitude(ni-xIndex+2,nj-2)
604 return
605 end if
606 end do
607 do xIndex = ni/2 + 2, ni
608 if(latitude(xIndex,nj) /= latitude(ni-xIndex+2,nj-2)) then
609 gridIsOrca = .false.
610 write(*,*) '3. Not an orca grid',xIndex,latitude(xIndex,nj),latitude(ni-xIndex+2,nj-2)
611 return
612 end if
613 end do
614
615 end function gpos_gridIsOrca
616
617end module getGridPosition_mod
