FMS 2025.01-dev
Flexible Modeling System
Loading...
Searching...
No Matches
horiz_interp_spherical.inc
1!***********************************************************************
2!* GNU Lesser General Public License
3!*
4!* This file is part of the GFDL Flexible Modeling System (FMS).
5!*
6!* FMS is free software: you can redistribute it and/or modify it under
7!* the terms of the GNU Lesser General Public License as published by
8!* the Free Software Foundation, either version 3 of the License, or (at
9!* your option) any later version.
10!*
11!* FMS is distributed in the hope that it will be useful, but WITHOUT
12!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13!* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14!* for more details.
15!*
16!* You should have received a copy of the GNU Lesser General Public
17!* License along with FMS. If not, see <http://www.gnu.org/licenses/>.
18!***********************************************************************
19!> @addtogroup horiz_interp_spherical_mod
20!> @{
21 !> Initialization routine.
22 !!
23 !> Allocates space and initializes a derived-type variable
24 !! that contains pre-computed interpolation indices and weights.
25 subroutine horiz_interp_spherical_new_(Interp, lon_in,lat_in,lon_out,lat_out, &
26 num_nbrs, max_dist, src_modulo)
27
28 type(horiz_interp_type), intent(inout) :: Interp !< A derived type variable containing indices
29 !! and weights for subsequent interpolations. To
30 !! reinitialize for different grid-to-grid interpolation
31 !! @ref horiz_interp_del must be used first.
32 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_in !< Latitude (radians) for source data grid
33 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lat_in !< Longitude (radians) for source data grid
34 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lon_out !< Longitude (radians) for output data grid
35 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: lat_out !< Latitude (radians) for output data grid
36 logical, intent(in), optional :: src_modulo !< indicates if the boundary condition
37 !! along zonal boundary is cyclic or not. Cyclic when true
38 integer, intent(in), optional :: num_nbrs !< Number of nearest neighbors for regridding
39 !! When number of neighbors within the radius max_dist
40 !! is less than num_nbrs, All the neighbors will be used
41 !! to interpolate onto destination grid. when number of
42 !! neighbors within the radius max_dist is greater than
43 !! num_nbrs, at least "num_nbrs"
44 ! neighbors will be used to remap onto destination grid
45 real(FMS_HI_KIND_), optional, intent(in) :: max_dist !< Maximum region of influence around
46 !! destination grid points
47
48 !------local variables ---------------------------------------
49 integer :: i, j, n
50 integer :: map_dst_xsize, map_dst_ysize, map_src_xsize, map_src_ysize
51 integer :: map_src_size, num_neighbors
52 real(FMS_HI_KIND_) :: max_src_dist, tpi, hpi
53 logical :: src_is_modulo
54 real(FMS_HI_KIND_) :: min_theta_dst, max_theta_dst, min_phi_dst, max_phi_dst
55 real(FMS_HI_KIND_) :: min_theta_src, max_theta_src, min_phi_src, max_phi_src
56 integer :: map_src_add(size(lon_out,1),size(lon_out,2),max_neighbors)
57 real(FMS_HI_KIND_) :: map_src_dist(size(lon_out,1),size(lon_out,2),max_neighbors)
58 integer :: num_found(size(lon_out,1),size(lon_out,2))
59 integer :: ilon(max_neighbors), jlat(max_neighbors)
60 real(FMS_HI_KIND_), dimension(size(lon_out,1),size(lon_out,2)) :: theta_dst, phi_dst
61 real(FMS_HI_KIND_), dimension(size(lon_in,1)*size(lon_in,2)) :: theta_src, phi_src
62 integer, parameter :: kindl = fms_hi_kind_
63
64 !--------------------------------------------------------------
65
66 pe = mpp_pe()
67 root_pe = mpp_root_pe()
68
69 tpi = 2.0_kindl*real(pi, fms_hi_kind_); hpi = 0.5_kindl*real(pi,fms_hi_kind_)
70
71 num_neighbors = num_nbrs_default
72 if(present(num_nbrs)) num_neighbors = num_nbrs
73 if (num_neighbors <= 0) call mpp_error(fatal,'horiz_interp_spherical_mod: num_neighbors must be > 0')
74
75 max_src_dist = real(max_dist_default, fms_hi_kind_)
76 if (PRESENT(max_dist)) max_src_dist = max_dist
77 interp%HI_KIND_TYPE_%max_src_dist = max_src_dist
78
79 src_is_modulo = .true.
80 if (PRESENT(src_modulo)) src_is_modulo = src_modulo
81
82 !--- check the grid size comformable
83 map_dst_xsize=size(lon_out,1);map_dst_ysize=size(lon_out,2)
84 map_src_xsize=size(lon_in,1); map_src_ysize=size(lon_in,2)
85 map_src_size = map_src_xsize*map_src_ysize
86
87 if (map_dst_xsize /= size(lat_out,1) .or. map_dst_ysize /= size(lat_out,2)) &
88 call mpp_error(fatal,'horiz_interp_spherical_mod: destination grids not conformable')
89 if (map_src_xsize /= size(lat_in,1) .or. map_src_ysize /= size(lat_in,2)) &
90 call mpp_error(fatal,'horiz_interp_spherical_mod: source grids not conformable')
91
92 theta_src = reshape(lon_in,(/map_src_size/))
93 phi_src = reshape(lat_in,(/map_src_size/))
94 theta_dst(:,:) = lon_out(:,:)
95 phi_dst(:,:) = lat_out(:,:)
96
97 min_theta_dst=real(tpi, fms_hi_kind_);max_theta_dst=0.0_kindl
98 min_phi_dst=real(pi, fms_hi_kind_);max_phi_dst=real(-pi, fms_hi_kind_)
99 min_theta_src=real(tpi, fms_hi_kind_);max_theta_src=0.0_kindl
100 min_phi_src=real(pi, fms_hi_kind_);max_phi_src=real(-pi, fms_hi_kind_)
101
102 where(theta_dst<0.0_kindl) theta_dst = theta_dst+real(tpi,fms_hi_kind_)
103
104 where(theta_dst>real(tpi,fms_hi_kind_)) theta_dst = theta_dst-real(tpi,fms_hi_kind_)
105 where(theta_src<0.0_kindl) theta_src = theta_src+real(tpi,fms_hi_kind_)
106 where(theta_src>real(tpi,fms_hi_kind_)) theta_src = theta_src-real(tpi,fms_hi_kind_)
107
108 where(phi_dst < -hpi) phi_dst = -hpi
109 where(phi_dst > hpi) phi_dst = hpi
110 where(phi_src < -hpi) phi_src = -hpi
111 where(phi_src > hpi) phi_src = hpi
112
113 do j=1,map_dst_ysize
114 do i=1,map_dst_xsize
115 min_theta_dst = min(min_theta_dst,theta_dst(i,j))
116 max_theta_dst = max(max_theta_dst,theta_dst(i,j))
117 min_phi_dst = min(min_phi_dst,phi_dst(i,j))
118 max_phi_dst = max(max_phi_dst,phi_dst(i,j))
119 enddo
120 enddo
121
122 do i=1,map_src_size
123 min_theta_src = min(min_theta_src,theta_src(i))
124 max_theta_src = max(max_theta_src,theta_src(i))
125 min_phi_src = min(min_phi_src,phi_src(i))
126 max_phi_src = max(max_phi_src,phi_src(i))
127 enddo
128
129 if (min_phi_dst < min_phi_src) print *, '=> WARNING: latitute of dest grid exceeds src'
130 if (max_phi_dst > max_phi_src) print *, '=> WARNING: latitute of dest grid exceeds src'
131 ! when src is cyclic, no need to print out the following warning.
132 if(.not. src_is_modulo) then
133 if (min_theta_dst < min_theta_src) print *, '=> WARNING : longitude of dest grid exceeds src'
134 if (max_theta_dst > max_theta_src) print *, '=> WARNING : longitude of dest grid exceeds src'
135 endif
136
137 ! allocate memory to data type
138 if(ALLOCATED(interp%i_lon)) then
139 if(size(interp%i_lon,1) .NE. map_dst_xsize .OR. &
140 size(interp%i_lon,2) .NE. map_dst_ysize ) call mpp_error(fatal, &
141 .NE..OR.'horiz_interp_spherical_mod: size(Interp%i_lon(:),1) map_dst_xsize '// &
142 .NE.'size(Interp%i_lon(:),2) map_dst_ysize')
143 else
144 allocate(interp%i_lon(map_dst_xsize,map_dst_ysize,max_neighbors), &
145 interp%j_lat(map_dst_xsize,map_dst_ysize,max_neighbors), &
146 interp%HI_KIND_TYPE_%src_dist(map_dst_xsize,map_dst_ysize,max_neighbors), &
147 interp%num_found(map_dst_xsize,map_dst_ysize) )
148 endif
149
150 map_src_add = 0
151 map_src_dist = real(large, fms_hi_kind_)
152 num_found = 0
153
154 !using radial_search to find the nearest points and corresponding distance.
155
156 select case(trim(search_method))
157 case ("radial_search") ! will be efficient, but may be not so accurate for some cases
158 call radial_search(theta_src, phi_src, theta_dst, phi_dst, map_src_xsize, map_src_ysize, &
159 map_src_add, map_src_dist, num_found, num_neighbors,max_src_dist,src_is_modulo)
160 case ("full_search") ! always accurate, but less efficient.
161 call full_search(theta_src, phi_src, theta_dst, phi_dst, map_src_add, map_src_dist, &
162 num_found, num_neighbors,max_src_dist )
163 case default
164 call mpp_error(fatal,"HORIZ_INTERP_SPHERICAL_NEW_: nml search_method = "// &
165 trim(search_method)//" is not a valid namelist option")
166 end select
167
168 do j=1,map_dst_ysize
169 do i=1,map_dst_xsize
170 do n=1,num_found(i,j)
171 if(map_src_add(i,j,n) == 0) then
172 jlat(n) = 0; ilon(n) = 0
173 else
174 jlat(n) = map_src_add(i,j,n)/map_src_xsize + 1
175 ilon(n) = map_src_add(i,j,n) - (jlat(n)-1)*map_src_xsize
176 if(ilon(n) == 0) then
177 jlat(n) = jlat(n) - 1
178 ilon(n) = map_src_xsize
179 endif
180 endif
181 enddo
182 interp%i_lon(i,j,:) = ilon(:)
183 interp%j_lat(i,j,:) = jlat(:)
184 interp%num_found(i,j) = num_found(i,j)
185 interp%HI_KIND_TYPE_%src_dist(i,j,:) = map_src_dist(i,j,:)
186 enddo
187 enddo
188
189 interp%nlon_src = map_src_xsize; interp%nlat_src = map_src_ysize
190 interp%nlon_dst = map_dst_xsize; interp%nlat_dst = map_dst_ysize
191 interp% HI_KIND_TYPE_ % is_allocated = .true.
192 interp% interp_method = spherical
193
194 return
195
196 end subroutine horiz_interp_spherical_new_
197
198 !#######################################################################
199
200 !> Subroutine for performing the horizontal interpolation between two grids.
201 !! HORIZ_INTERP_SPHERICAL_NEW_ must be called before calling this routine.
202 subroutine horiz_interp_spherical_( Interp, data_in, data_out, verbose, mask_in, mask_out, missing_value)
203 type(horiz_interp_type), intent(in) :: Interp !< A derived type variable containing indices
204 !! and weights for subsequent interpolations. Returned
205 !! by a previous call to HORIZ_INTERP_SPHERICAL_NEW_
206 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: data_in !< Input data on source grid
207 real(FMS_HI_KIND_), intent(out), dimension(:,:) :: data_out !< Output data on destination grid
208 integer, intent(in), optional :: verbose !< 0 = no output; 1 = min,max,means; 2 = most output
209 real(FMS_HI_KIND_), intent(in), dimension(:,:), optional :: mask_in !< Input mask, must be the same size as
210 !! the input data. The real(FMS_HI_KIND_) value of mask_in must be
211 !! in the range (0.,1.). Set mask_in=0.0 for data points
212 !! that should not be used or have missing data
213 real(FMS_HI_KIND_), intent(out), dimension(:,:), optional :: mask_out !< Output mask that specifies whether data
214 !! was computed.
215 real(FMS_HI_KIND_), intent(in), optional :: missing_value !< Used to indicate missing data
216
217 !--- some local variables ----------------------------------------
218 real(FMS_HI_KIND_), dimension(Interp%nlon_dst, Interp%nlat_dst,size(Interp%HI_KIND_TYPE_%src_dist,3)) :: wt
219 real(FMS_HI_KIND_), dimension(Interp%nlon_src, Interp%nlat_src) :: mask_src
220 real(FMS_HI_KIND_), dimension(Interp%nlon_dst, Interp%nlat_dst) :: mask_dst
221 integer :: nlon_in, nlat_in, nlon_out, nlat_out, num_found
222 integer :: m, n, i, j, k, miss_in, miss_out, i1, i2, j1, j2, iverbose
223 real(FMS_HI_KIND_) :: min_in, max_in, avg_in, min_out, max_out, avg_out, sum
224 integer, parameter :: kindl = fms_hi_kind_ !< compiled real kind size
225 !-----------------------------------------------------------------
226
227 iverbose = 0; if (present(verbose)) iverbose = verbose
228
229 nlon_in = interp%nlon_src; nlat_in = interp%nlat_src
230 nlon_out = interp%nlon_dst; nlat_out = interp%nlat_dst
231
232 if(size(data_in,1) .ne. nlon_in .or. size(data_in,2) .ne. nlat_in ) &
233 call mpp_error(fatal,'horiz_interp_spherical_mod: size of input array incorrect')
234
235 if(size(data_out,1) .ne. nlon_out .or. size(data_out,2) .ne. nlat_out ) &
236 call mpp_error(fatal,'horiz_interp_spherical_mod: size of output array incorrect')
237
238 mask_src = 1.0_kindl; mask_dst = 1.0_kindl
239 if(present(mask_in)) mask_src = mask_in
240
241 do n=1,nlat_out
242 do m=1,nlon_out
243 ! neighbors are sorted nearest to farthest
244 ! check nearest to see if it is a land point
245 num_found = interp%num_found(m,n)
246 if(num_found == 0 ) then
247 mask_dst(m,n) = 0.0_kindl
248 else
249 i1 = interp%i_lon(m,n,1); j1 = interp%j_lat(m,n,1)
250 if (mask_src(i1,j1) .lt. 0.5_kindl ) then
251 mask_dst(m,n) = 0.0_kindl
252 endif
253
254 if(num_found .gt. 1 ) then
255 i2 = interp%i_lon(m,n,2); j2 = interp%j_lat(m,n,2)
256 ! compare first 2 nearest neighbors -- if they are nearly
257 ! equidistant then use this mask for robustness
258 if(abs(interp%HI_KIND_TYPE_%src_dist(m,n,2)-interp%HI_KIND_TYPE_%src_dist(m,n,1)) .lt. &
259 real(epsln,FMS_HI_KIND_)) then
260 if((mask_src(i1,j1) .lt. 0.5_kindl )) mask_dst(m,n) = 0.0_kindl
261 endif
262 endif
263
264 sum=0.0_kindl
265 do k=1, num_found
266 if(mask_src(interp%i_lon(m,n,k),interp%j_lat(m,n,k)) .lt. 0.5_kindl ) then
267 wt(m,n,k) = 0.0_kindl
268 else
269 if (interp%HI_KIND_TYPE_%src_dist(m,n,k) <= real(epsln,fms_hi_kind_)) then
270 wt(m,n,k) = real(large, fms_hi_kind_)
271 sum = sum + real(large, fms_hi_kind_)
272 else if(interp%HI_KIND_TYPE_%src_dist(m,n,k) <= interp%HI_KIND_TYPE_%max_src_dist ) then
273 wt(m,n,k) = 1.0_kindl /interp%HI_KIND_TYPE_%src_dist(m,n,k)
274 sum = sum+wt(m,n,k)
275 else
276 wt(m,n,k) = 0.0_kindl
277 endif
278 endif
279 enddo
280 if (sum > real(epsln,fms_hi_kind_)) then
281 do k = 1, num_found
282 wt(m,n,k) = wt(m,n,k)/sum
283 enddo
284 else
285 mask_dst(m,n) = 0.0_kindl
286 endif
287 endif
288 enddo
289 enddo
290
291 data_out = 0.0_kindl
292 do n=1,nlat_out
293 do m=1,nlon_out
294 if(mask_dst(m,n) .gt. 0.5_kindl ) then
295 do k=1, interp%num_found(m,n)
296 i = interp%i_lon(m,n,k)
297 j = interp%j_lat(m,n,k)
298 data_out(m,n) = data_out(m,n)+data_in(i,j)*wt(m,n,k)
299 enddo
300 else
301 if(present(missing_value)) then
302 data_out(m,n) = missing_value
303 else
304 data_out(m,n) = 0.0_kindl
305 endif
306 endif
307 enddo
308 enddo
309
310 if(present(mask_out)) mask_out = mask_dst
311
312 !***********************************************************************
313 ! compute statistics: minimum, maximum, and mean
314 !-----------------------------------------------------------------------
315
316 if (iverbose > 0) then
317
318 ! compute statistics of input data
319
320 call stats (data_in, min_in, max_in, avg_in, miss_in, missing_value, mask=mask_src)
321
322 ! compute statistics of output data
323 call stats (data_out, min_out, max_out, avg_out, miss_out, missing_value, mask=mask_dst)
324
325 !---- output statistics ----
326 ! root_pe have the information of global mean, min and max
327 if(pe == root_pe) then
328 write (*,900)
329 write (*,901) min_in ,max_in, avg_in
330 if (present(mask_in)) write (*,903) miss_in
331 write (*,902) min_out,max_out,avg_out
332 if (present(mask_out)) write (*,903) miss_out
333 endif
334900 format (/,1x,10('-'),' output from horiz_interp ',10('-'))
335901 format (' input: min=',f16.9,' max=',f16.9,' avg=',f22.15)
336902 format (' output: min=',f16.9,' max=',f16.9,' avg=',f22.15)
337903 format (' number of missing points = ',i6)
338
339 endif
340
341 return
342 end subroutine horiz_interp_spherical_
343
344 !#######################################################################
345 !> This routine isn't used internally
346 !! it's similar to the routine above, just gets the weights as an out variable
347 subroutine horiz_interp_spherical_wght_( Interp, wt, verbose, mask_in, mask_out, missing_value)
348 type (horiz_interp_type), intent(in) :: Interp
349 real(FMS_HI_KIND_), intent(out), dimension(:,:,:) :: wt
350 integer, intent(in), optional :: verbose
351 real(FMS_HI_KIND_), intent(in), dimension(:,:), optional :: mask_in
352 real(FMS_HI_KIND_), intent(inout), dimension(:,:), optional :: mask_out
353 real(FMS_HI_KIND_), intent(in), optional :: missing_value
354
355 !--- some local variables ----------------------------------------
356 real(FMS_HI_KIND_), dimension(Interp%nlon_src, Interp%nlat_src) :: mask_src
357 real(FMS_HI_KIND_), dimension(Interp%nlon_dst, Interp%nlat_dst) :: mask_dst
358 integer :: nlon_in, nlat_in, nlon_out, nlat_out, num_found
359 integer :: m, n, k, i1, i2, j1, j2, iverbose
360 real(FMS_HI_KIND_) :: sum
361 integer, parameter :: kindl = fms_hi_kind_
362 !-----------------------------------------------------------------
363
364 iverbose = 0; if (present(verbose)) iverbose = verbose
365
366 nlon_in = interp%nlon_src; nlat_in = interp%nlat_src
367 nlon_out = interp%nlon_dst; nlat_out = interp%nlat_dst
368
369 mask_src = 1.0_kindl ; mask_dst = 1.0_kindl
370 if(present(mask_in)) mask_src = mask_in
371
372 do n=1,nlat_out
373 do m=1,nlon_out
374 ! neighbors are sorted nearest to farthest
375 ! check nearest to see if it is a land point
376 num_found = interp%num_found(m,n)
377
378 if (num_found > num_nbrs_default) then
379 if( iverbose .gt. 0) print *,'pe=',mpp_pe(),'num_found=',num_found
380 num_found = num_nbrs_default
381 end if
382
383 if(num_found == 0 ) then
384 mask_dst(m,n) = 0.0_kindl
385 else
386 i1 = interp%i_lon(m,n,1); j1 = interp%j_lat(m,n,1)
387 if (mask_src(i1,j1) .lt. 0.5_kindl) then
388 mask_dst(m,n) = 0.0_kindl
389 endif
390
391 if(num_found .gt. 1 ) then
392 i2 = interp%i_lon(m,n,2); j2 = interp%j_lat(m,n,2)
393 ! compare first 2 nearest neighbors -- if they are nearly
394 ! equidistant then use this mask for robustness
395 if(abs(interp%HI_KIND_TYPE_%src_dist(m,n,2)-interp%HI_KIND_TYPE_%src_dist(m,n,1)) .lt. &
396 real(epsln,FMS_HI_KIND_)) then
397 if((mask_src(i1,j1) .lt. 0.5_kindl )) mask_dst(m,n) = 0.0_kindl
398 endif
399 endif
400
401 sum=0.0_kindl
402 do k=1, num_found
403 if(mask_src(interp%i_lon(m,n,k),interp%j_lat(m,n,k)) .lt. 0.5_kindl ) then
404 wt(m,n,k) = 0.0_kindl
405 else
406 if (interp%HI_KIND_TYPE_%src_dist(m,n,k) <= real(epsln, fms_hi_kind_)) then
407 wt(m,n,k) = real(large, fms_hi_kind_)
408 sum = sum + real(large, fms_hi_kind_)
409 else if(interp%HI_KIND_TYPE_%src_dist(m,n,k) <= interp%HI_KIND_TYPE_%max_src_dist ) then
410 wt(m,n,k) = 1.0_kindl /interp%HI_KIND_TYPE_%src_dist(m,n,k)
411 sum = sum+wt(m,n,k)
412 else
413 wt(m,n,k) = 0.0_kindl
414 endif
415 endif
416 enddo
417 if (sum > real(epsln,fms_hi_kind_)) then
418 do k = 1, num_found
419 wt(m,n,k) = wt(m,n,k)/sum
420 enddo
421 else
422 mask_dst(m,n) = 0.0_kindl
423 endif
424 endif
425 enddo
426 enddo
427
428 return
429 end subroutine horiz_interp_spherical_wght_
430
431 !#######################################################################
432
433 subroutine radial_search_(theta_src,phi_src,theta_dst,phi_dst, map_src_xsize, map_src_ysize, &
434 map_src_add, map_src_dist, num_found, num_neighbors,max_src_dist,src_is_modulo)
435 real(FMS_HI_KIND_), intent(in), dimension(:) :: theta_src, phi_src
436 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: theta_dst, phi_dst
437 integer, intent(in) :: map_src_xsize, map_src_ysize
438 integer, intent(out), dimension(:,:,:) :: map_src_add
439 real(FMS_HI_KIND_), intent(out), dimension(:,:,:) :: map_src_dist
440 integer, intent(inout), dimension(:,:) :: num_found
441 integer, intent(in) :: num_neighbors
442 real(FMS_HI_KIND_), intent(in) :: max_src_dist
443 logical, intent(in) :: src_is_modulo
444
445 !---------- local variables ----------------------------------------
446 integer, parameter :: max_nbrs = 50
447 integer :: i, j, jj, i0, j0, n, l,i_left, i_right
448 integer :: map_dst_xsize, map_dst_ysize
449 integer :: i_left1, i_left2, i_right1, i_right2
450 integer :: map_src_size, step, step_size, bound, bound_start, bound_end
451 logical :: continue_search, result, continue_radial_search
452 real(FMS_HI_KIND_) :: d, res
453 !------------------------------------------------------------------
454 map_dst_xsize=size(theta_dst,1);map_dst_ysize=size(theta_dst,2)
455 map_src_size = map_src_xsize*map_src_ysize
456
457 do j=1,map_dst_ysize
458 do i=1,map_dst_xsize
459 continue_search=.true.
460 step = 1
461 step_size = int( sqrt(real(map_src_size, kind=fms_hi_kind_ )))
462 do while (continue_search .and. step_size > 0)
463 do while (step <= map_src_size .and. continue_search)
464 ! count land points as nearest neighbors
465 d = spherical_distance(theta_dst(i,j),phi_dst(i,j),theta_src(step),phi_src(step))
466 if (d <= max_src_dist) then
467 result = update_dest_neighbors_(map_src_add(i,j,:),map_src_dist(i,j,:), &
468 step,d, num_found(i,j), num_neighbors )
469 if (result) then
470 n = 0
471 i0 = mod(step,map_src_xsize)
472
473 if (i0 == 0) i0 = map_src_xsize
474 res = real(step, fms_hi_kind_)/real(map_src_xsize, fms_hi_kind_)
475 j0 = ceiling(res)
476 continue_radial_search = .true.
477 do while (continue_radial_search)
478 continue_radial_search = .false.
479 n = n+1 ! radial counter
480 if(n > max_nbrs) exit
481 ! ************** left boundary *******************************
482 i_left = i0-n
483 if (i_left <= 0) then
484 if (src_is_modulo) then
485 i_left = map_src_xsize + i_left
486 else
487 i_left = 1
488 endif
489 endif
490
491 do l = 0, 2*n
492 jj = j0 - n - 1 + l
493 if( jj < 0) then
494 bound = ( 1 - jj )*map_src_xsize - i_left
495 else if ( jj >= map_src_ysize ) then
496 bound = ( 2*map_src_ysize - jj ) * map_src_xsize - i_left
497 else
498 bound = jj * map_src_xsize + i_left
499 endif
500
501 d = spherical_distance(theta_dst(i,j),phi_dst(i,j),theta_src(bound), &
502 phi_src(bound))
503 if(d<=max_src_dist) then
504 result = update_dest_neighbors_(map_src_add(i,j,:),map_src_dist(i,j,:), &
505 bound,d, num_found(i,j), num_neighbors)
506 if (result) continue_radial_search = .true.
507 endif
508 enddo
509
510 ! ***************************right boundary *******************************
511 i_right = i0+n
512 if (i_right > map_src_xsize) then
513 if (src_is_modulo) then
514 i_right = i_right - map_src_xsize
515 else
516 i_right = map_src_xsize
517 endif
518 endif
519
520 do l = 0, 2*n
521 jj = j0 - n - 1 + l
522 if( jj < 0) then
523 bound = ( 1 - jj )*map_src_xsize - i_right
524 else if ( jj >= map_src_ysize ) then
525 bound = ( 2*map_src_ysize - jj) * map_src_xsize - i_right
526
527 else
528 bound = jj * map_src_xsize + i_right
529 endif
530
531 d = spherical_distance(theta_dst(i,j),phi_dst(i,j),theta_src(bound), &
532 phi_src(bound))
533 if(d<=max_src_dist) then
534 result = update_dest_neighbors_(map_src_add(i,j,:),map_src_dist(i,j,:), &
535 bound,d, num_found(i,j), num_neighbors)
536 if (result) continue_radial_search = .true.
537 endif
538 enddo
539
540 ! ************************* bottom boundary **********************************
541 i_left2 = 0
542 if( i_left > i_right) then
543 i_left1 = 1
544 i_right1 = i_right
545 i_left2 = i_left
546 i_right2 = map_src_xsize
547 else
548 i_left1 = i_left
549 i_right1 = i_right
550 endif
551
552 jj = j0 - n - 1
553 if( jj < 0 ) then
554 bound_start = ( 1 - jj)*map_src_xsize - i_right1
555 bound_end = ( 1 - jj)*map_src_xsize - i_left1
556 else
557 bound_start = jj * map_src_xsize + i_left1
558 bound_end = jj * map_src_xsize + i_right1
559 endif
560
561 bound = bound_start
562 do while (bound <= bound_end)
563 d = spherical_distance(theta_dst(i,j),phi_dst(i,j),theta_src(bound), &
564 phi_src(bound))
565 if(d<=max_src_dist) then
566 result = update_dest_neighbors_(map_src_add(i,j,:),map_src_dist(i,j,:), &
567 bound,d, num_found(i,j), num_neighbors)
568 if (result) continue_radial_search = .true.
569 endif
570 bound = bound + 1
571
572 enddo
573
574 if(i_left2 > 0 ) then
575 if( jj < 0 ) then
576 bound_start = ( 1 - jj)*map_src_xsize - i_right2
577 bound_end = ( 1 - jj)*map_src_xsize - i_left2
578 else
579 bound_start = jj * map_src_xsize + i_left2
580 bound_end = jj * map_src_xsize + i_right2
581 endif
582
583 bound = bound_start
584 do while (bound <= bound_end)
585 d = spherical_distance(theta_dst(i,j),phi_dst(i,j),theta_src(bound), &
586 phi_src(bound))
587 if(d<=max_src_dist) then
588 result = update_dest_neighbors_(map_src_add(i,j,:),map_src_dist(i,j,:), &
589 bound,d, num_found(i,j), num_neighbors)
590 if (result) continue_radial_search = .true.
591 endif
592 bound = bound + 1
593 enddo
594 endif
595
596 ! ************************** top boundary ************************************
597 jj = j0 + n - 1
598 if( jj >= map_src_ysize) then
599 bound_start = ( 2*map_src_ysize - jj ) * map_src_xsize - i_right1
600 bound_end = ( 2*map_src_ysize - jj ) * map_src_xsize - i_left1
601 else
602 bound_start = jj * map_src_xsize + i_left1
603 bound_end = jj * map_src_xsize + i_right1
604 endif
605
606 bound = bound_start
607 do while (bound <= bound_end)
608 d = spherical_distance(theta_dst(i,j),phi_dst(i,j),theta_src(bound), &
609 phi_src(bound))
610 if(d<=max_src_dist) then
611 result = update_dest_neighbors_(map_src_add(i,j,:),map_src_dist(i,j,:), &
612 bound,d, num_found(i,j), num_neighbors)
613 if (result) continue_radial_search = .true.
614 endif
615 bound = bound + 1
616 enddo
617
618 if(i_left2 > 0) then
619 if( jj >= map_src_ysize) then
620 bound_start = ( 2*map_src_ysize - jj ) * map_src_xsize - i_right2
621 bound_end = ( 2*map_src_ysize - jj ) * map_src_xsize - i_left2
622 else
623 bound_start = jj * map_src_xsize + i_left2
624 bound_end = jj * map_src_xsize + i_right2
625 endif
626
627 bound = bound_start
628 do while (bound <= bound_end)
629 d = spherical_distance(theta_dst(i,j),phi_dst(i,j),theta_src(bound), &
630 phi_src(bound))
631 if(d<=max_src_dist) then
632 result = update_dest_neighbors_(map_src_add(i,j,:),map_src_dist(i,j,:), &
633 bound,d, num_found(i,j), num_neighbors)
634 if (result) continue_radial_search = .true.
635 endif
636 bound = bound + 1
637 enddo
638 endif
639
640 enddo
641 continue_search = .false. ! stop looking
642 endif
643 endif
644 step=step+step_size
645 enddo ! search loop
646 step = 1
647 step_size = step_size/2
648 enddo
649 enddo
650 enddo
651
652 return
653
654 end subroutine radial_search_
655
656
657 !#####################################################################
658
659 function update_dest_neighbors_(map_src_add, map_src_dist, src_add,d, num_found, min_nbrs)
660
661 integer, intent(inout), dimension(:) :: map_src_add
662 real(FMS_HI_KIND_), intent(inout), dimension(:) :: map_src_dist
663 integer, intent(in) :: src_add
664 real(FMS_HI_KIND_), intent(in) :: d
665 integer, intent(inout) :: num_found
666 integer, intent(in) :: min_nbrs
667
668 logical :: UPDATE_DEST_NEIGHBORS_, already_exist = .false.
669
670 integer :: n,m
671
672 update_dest_neighbors_ = .false.
673
674 n = 0
675 nloop : do while ( n .le. num_found )
676 n = n + 1
677 dist_chk : if (d .le. map_src_dist(n)) then
678 do m=n,num_found
679 if (src_add == map_src_add(m)) then
680 already_exist = .true.
681 exit nloop
682 endif
683 enddo
684 if(num_found < max_neighbors) then
685 num_found = num_found + 1
686 else
687 call mpp_error(fatal,'UPDATE_DEST_NEIGHBORS_: '// &
688 'number of neighbor points found is greated than maxium neighbor points' )
689 endif
690 do m=num_found,n+1,-1
691 map_src_add(m) = map_src_add(m-1)
692 map_src_dist(m) = map_src_dist(m-1)
693 enddo
694 map_src_add(n) = src_add
695 map_src_dist(n) = d
696 update_dest_neighbors_ = .true.
697 if( num_found > min_nbrs ) then
698 if( map_src_dist(num_found) > map_src_dist(num_found-1) ) then
699 num_found = num_found - 1
700 endif
701 if( map_src_dist(min_nbrs+1) > map_src_dist(min_nbrs) ) then
702 num_found = min_nbrs
703 endif
704 endif
705 exit nloop ! n loop
706 endif dist_chk
707 end do nloop
708 if(already_exist) return
709
710 if( .not. update_dest_neighbors_ ) then
711 if( num_found < min_nbrs ) then
712 num_found = num_found + 1
713 update_dest_neighbors_ = .true.
714 map_src_add(num_found) = src_add
715 map_src_dist(num_found) = d
716 endif
717 endif
718
719
720 return
721
722 end function update_dest_neighbors_
723
724 !########################################################################
725! function HORIZ_INTERP_SPHERICAL_DISTANCE_(theta1,phi1,theta2,phi2)
726
727! real(FMS_HI_KIND_), intent(in) :: theta1, phi1, theta2, phi2
728! real(FMS_HI_KIND_) :: HORIZ_INTERP_SPHERICAL_DISTANCE_
729
730! real(FMS_HI_KIND_) :: r1(3), r2(3), cross(3), s, dot, ang
731
732 ! this is a simple, enough way to calculate distance on the sphere
733 ! first, construct cartesian vectors r1 and r2
734 ! then calculate the cross-product which is proportional to the area
735 ! between the 2 vectors. The angular distance is arcsin of the
736 ! distancealong the sphere
737 !
738 ! theta is longitude and phi is latitude
739 !
740
741
742! r1(1) = cos(theta1)*cos(phi1);r1(2)=sin(theta1)*cos(phi1);r1(3)=sin(phi1)
743! r2(1) = cos(theta2)*cos(phi2);r2(2)=sin(theta2)*cos(phi2);r2(3)=sin(phi2)
744
745! cross(1) = r1(2)*r2(3)-r1(3)*r2(2)
746! cross(2) = r1(3)*r2(1)-r1(1)*r2(3)
747! cross(3) = r1(1)*r2(2)-r1(2)*r2(1)
748
749! s = sqrt(cross(1)**2.+cross(2)**2.+cross(3)**2.)
750
751! s = min(s,real(1.0, FMS_HI_KIND_)-epsln)
752
753! dot = r1(1)*r2(1) + r1(2)*r2(2) + r1(3)*r2(3)
754
755! if (dot > 0) then
756! ang = asin(s)
757! else if (dot < 0) then
758! ang = pi + asin(s) !? original is pi - asin(s)
759! else
760! ang = pi/real(2., FMS_HI_KIND_)
761! endif
762
763! HORIZ_INTERP_SPHERICAL_DISTANCE_ = abs(ang) ! in radians
764
765! return
766
767! end function HORIZ_INTERP_SPHERICAL_DISTANCE_
768 ! The great cycle distance
769 function horiz_interp_spherical_distance_(theta1,phi1,theta2,phi2)
770
771 real(FMS_HI_KIND_), intent(in) :: theta1, phi1, theta2, phi2
772 real(FMS_HI_KIND_) :: HORIZ_INTERP_SPHERICAL_DISTANCE_, dot
773 integer, parameter :: kindl = fms_hi_kind_
774
775 if(theta1 == theta2 .and. phi1 == phi2) then
776 horiz_interp_spherical_distance_ = 0.0_kindl
777 return
778 endif
779
780 dot = cos(phi1)*cos(phi2)*cos(theta1-theta2) + sin(phi1)*sin(phi2)
781 if(dot > 1.0_kindl) dot = 1.0_kindl
782 if(dot < real(-1.0_kindl, fms_hi_kind_)) dot = -1.0_kindl
783 horiz_interp_spherical_distance_ = acos(dot)
784
785 return
786
787 end function horiz_interp_spherical_distance_
788
789
790 !#######################################################################
791
792 subroutine full_search_(theta_src,phi_src,theta_dst,phi_dst,map_src_add, map_src_dist,num_found, &
793 num_neighbors,max_src_dist)
794 real(FMS_HI_KIND_), intent(in), dimension(:) :: theta_src, phi_src
795 real(FMS_HI_KIND_), intent(in), dimension(:,:) :: theta_dst, phi_dst
796 integer, intent(out), dimension(:,:,:) :: map_src_add
797 real(FMS_HI_KIND_), intent(out), dimension(:,:,:) :: map_src_dist
798 integer, intent(out), dimension(:,:) :: num_found
799 integer, intent(in) :: num_neighbors
800 real(FMS_HI_KIND_), intent(in) :: max_src_dist
801
802 integer :: i,j,map_src_size, step
803 integer :: map_dst_xsize,map_dst_ysize
804 real(FMS_HI_KIND_) :: d
805 logical :: found
806
807 map_dst_xsize=size(theta_dst,1);map_dst_ysize=size(theta_dst,2)
808 map_src_size =size(theta_src(:))
809
810 do j=1,map_dst_ysize
811 do i=1,map_dst_xsize
812 do step = 1, map_src_size
813 d = spherical_distance(theta_dst(i,j),phi_dst(i,j),theta_src(step),phi_src(step))
814 if( d <= max_src_dist) then
815 found = update_dest_neighbors_(map_src_add(i,j,:),map_src_dist(i,j,:), &
816 step,d,num_found(i,j), num_neighbors )
817 endif
818 enddo
819 enddo
820 enddo
821
822 end subroutine full_search_
823!> @}
horiz_interp_spherical_wght_
subroutine horiz_interp_spherical_wght_(interp, wt, verbose, mask_in, mask_out, missing_value)
This routine isn't used internally it's similar to the routine above, just gets the weights as an out...
Definition horiz_interp_spherical.inc:348
horiz_interp_spherical_new_
subroutine horiz_interp_spherical_new_(interp, lon_in, lat_in, lon_out, lat_out, num_nbrs, max_dist, src_modulo)
Definition horiz_interp_spherical.inc:27
horiz_interp_spherical_
subroutine horiz_interp_spherical_(interp, data_in, data_out, verbose, mask_in, mask_out, missing_value)
Subroutine for performing the horizontal interpolation between two grids. HORIZ_INTERP_SPHERICAL_NEW_...
Definition horiz_interp_spherical.inc:203