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