25 num_nbrs, max_dist, src_modulo)
27 type(horiz_interp_type),
intent(inout) :: Interp
31 real(FMS_HI_KIND_),
intent(in),
dimension(:,:) :: lon_in
32 real(FMS_HI_KIND_),
intent(in),
dimension(:,:) :: lat_in
33 real(FMS_HI_KIND_),
intent(in),
dimension(:,:) :: lon_out
34 real(FMS_HI_KIND_),
intent(in),
dimension(:,:) :: lat_out
35 logical,
intent(in),
optional :: src_modulo
37 integer,
intent(in),
optional :: num_nbrs
44 real(FMS_HI_KIND_),
optional,
intent(in) :: max_dist
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_
66 root_pe = mpp_root_pe()
68 tpi = 2.0_kindl*real(pi, fms_hi_kind_); hpi = 0.5_kindl*real(pi,fms_hi_kind_)
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')
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
78 src_is_modulo = .true.
79 if (
PRESENT(src_modulo)) src_is_modulo = src_modulo
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
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')
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(:,:)
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_)
101 where(theta_dst<0.0_kindl) theta_dst = theta_dst+real(tpi,fms_hi_kind_)
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_)
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
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))
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))
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'
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'
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')
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) )
150 map_src_dist = real(large, fms_hi_kind_)
155 select case(trim(search_method))
156 case (
"radial_search")
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)
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 )
163 call mpp_error(fatal,
"HORIZ_INTERP_SPHERICAL_NEW_: nml search_method = "// &
164 trim(search_method)//
" is not a valid namelist option")
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
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
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,:)
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
202 type(horiz_interp_type),
intent(in) :: Interp
205 real(FMS_HI_KIND_),
intent(in),
dimension(:,:) :: data_in
206 real(FMS_HI_KIND_),
intent(out),
dimension(:,:) :: data_out
207 integer,
intent(in),
optional :: verbose
208 real(FMS_HI_KIND_),
intent(in),
dimension(:,:),
optional :: mask_in
212 real(FMS_HI_KIND_),
intent(out),
dimension(:,:),
optional :: mask_out
214 real(FMS_HI_KIND_),
intent(in),
optional :: missing_value
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_
226 iverbose = 0;
if (
present(verbose)) iverbose = verbose
228 nlon_in = interp%nlon_src; nlat_in = interp%nlat_src
229 nlon_out = interp%nlon_dst; nlat_out = interp%nlat_dst
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')
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')
237 mask_src = 1.0_kindl; mask_dst = 1.0_kindl
238 if(
present(mask_in)) mask_src = mask_in
244 num_found = interp%num_found(m,n)
245 if(num_found == 0 )
then
246 mask_dst(m,n) = 0.0_kindl
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
253 if(num_found .gt. 1 )
then
254 i2 = interp%i_lon(m,n,2); j2 = interp%j_lat(m,n,2)
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
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
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)
275 wt(m,n,k) = 0.0_kindl
279 if (sum > real(epsln,fms_hi_kind_))
then
281 wt(m,n,k) = wt(m,n,k)/sum
284 mask_dst(m,n) = 0.0_kindl
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)
300 if(
present(missing_value))
then
301 data_out(m,n) = missing_value
303 data_out(m,n) = 0.0_kindl
309 if(
present(mask_out)) mask_out = mask_dst
315 if (iverbose > 0)
then
319 call stats (data_in, min_in, max_in, avg_in, miss_in, missing_value, mask=mask_src)
322 call stats (data_out, min_out, max_out, avg_out, miss_out, missing_value, mask=mask_dst)
326 if(pe == root_pe)
then
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
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)
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
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_
363 iverbose = 0;
if (
present(verbose)) iverbose = verbose
365 nlon_in = interp%nlon_src; nlat_in = interp%nlat_src
366 nlon_out = interp%nlon_dst; nlat_out = interp%nlat_dst
368 mask_src = 1.0_kindl ; mask_dst = 1.0_kindl
369 if(
present(mask_in)) mask_src = mask_in
375 num_found = interp%num_found(m,n)
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
382 if(num_found == 0 )
then
383 mask_dst(m,n) = 0.0_kindl
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
390 if(num_found .gt. 1 )
then
391 i2 = interp%i_lon(m,n,2); j2 = interp%j_lat(m,n,2)
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
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
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)
412 wt(m,n,k) = 0.0_kindl
416 if (sum > real(epsln,fms_hi_kind_))
then
418 wt(m,n,k) = wt(m,n,k)/sum
421 mask_dst(m,n) = 0.0_kindl
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
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
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
458 continue_search=.true.
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)
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 )
470 i0 = mod(step,map_src_xsize)
472 if (i0 == 0) i0 = map_src_xsize
473 res = real(step, fms_hi_kind_)/real(map_src_xsize, fms_hi_kind_)
475 continue_radial_search = .true.
476 do while (continue_radial_search)
477 continue_radial_search = .false.
479 if(n > max_nbrs)
exit
482 if (i_left <= 0)
then
483 if (src_is_modulo)
then
484 i_left = map_src_xsize + i_left
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
497 bound = jj * map_src_xsize + i_left
500 d = spherical_distance(theta_dst(i,j),phi_dst(i,j),theta_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.
511 if (i_right > map_src_xsize)
then
512 if (src_is_modulo)
then
513 i_right = i_right - map_src_xsize
515 i_right = map_src_xsize
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
527 bound = jj * map_src_xsize + i_right
530 d = spherical_distance(theta_dst(i,j),phi_dst(i,j),theta_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.
541 if( i_left > i_right)
then
545 i_right2 = map_src_xsize
553 bound_start = ( 1 - jj)*map_src_xsize - i_right1
554 bound_end = ( 1 - jj)*map_src_xsize - i_left1
556 bound_start = jj * map_src_xsize + i_left1
557 bound_end = jj * map_src_xsize + i_right1
561 do while (bound <= bound_end)
562 d = spherical_distance(theta_dst(i,j),phi_dst(i,j),theta_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.
573 if(i_left2 > 0 )
then
575 bound_start = ( 1 - jj)*map_src_xsize - i_right2
576 bound_end = ( 1 - jj)*map_src_xsize - i_left2
578 bound_start = jj * map_src_xsize + i_left2
579 bound_end = jj * map_src_xsize + i_right2
583 do while (bound <= bound_end)
584 d = spherical_distance(theta_dst(i,j),phi_dst(i,j),theta_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.
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
601 bound_start = jj * map_src_xsize + i_left1
602 bound_end = jj * map_src_xsize + i_right1
606 do while (bound <= bound_end)
607 d = spherical_distance(theta_dst(i,j),phi_dst(i,j),theta_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.
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
622 bound_start = jj * map_src_xsize + i_left2
623 bound_end = jj * map_src_xsize + i_right2
627 do while (bound <= bound_end)
628 d = spherical_distance(theta_dst(i,j),phi_dst(i,j),theta_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.
640 continue_search = .false.
646 step_size = step_size/2
653 end subroutine radial_search_
658 function update_dest_neighbors_(map_src_add, map_src_dist, src_add,d, num_found, min_nbrs)
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
667 logical :: UPDATE_DEST_NEIGHBORS_, already_exist = .false.
671 update_dest_neighbors_ = .false.
674 nloop :
do while ( n .le. num_found )
676 dist_chk :
if (d .le. map_src_dist(n))
then
678 if (src_add == map_src_add(m))
then
679 already_exist = .true.
683 if(num_found < max_neighbors)
then
684 num_found = num_found + 1
686 call mpp_error(fatal,
'UPDATE_DEST_NEIGHBORS_: '// &
687 'number of neighbor points found is greated than maxium neighbor points' )
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)
693 map_src_add(n) = src_add
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
700 if( map_src_dist(min_nbrs+1) > map_src_dist(min_nbrs) )
then
707 if(already_exist)
return
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
721 end function update_dest_neighbors_
768 function horiz_interp_spherical_distance_(theta1,phi1,theta2,phi2)
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_
774 if(theta1 == theta2 .and. phi1 == phi2)
then
775 horiz_interp_spherical_distance_ = 0.0_kindl
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)
786 end function horiz_interp_spherical_distance_
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
801 integer :: i,j,map_src_size, step
802 integer :: map_dst_xsize,map_dst_ysize
803 real(FMS_HI_KIND_) :: d
806 map_dst_xsize=
size(theta_dst,1);map_dst_ysize=
size(theta_dst,2)
807 map_src_size =
size(theta_src(:))
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 )
821 end subroutine full_search_
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_...
subroutine horiz_interp_spherical_new_(Interp, lon_in, lat_in, lon_out, lat_out, num_nbrs, max_dist, src_modulo)
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...
integer function mpp_pe()
Returns processor ID.