FMS  2024.03
Flexible Modeling System
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
334 900 format (/,1x,10('-'),' output from horiz_interp ',10('-'))
335 901 format (' input: min=',f16.9,' max=',f16.9,' avg=',f22.15)
336 902 format (' output: min=',f16.9,' max=',f16.9,' avg=',f22.15)
337 903 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_
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
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_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
mpp_pe
integer function mpp_pe()
Returns processor ID.
Definition: mpp_util.inc:407