Back to OASIS3 home
Modules used for the scrip
library (in oasis3/lib/scrip/src) :
remap_distance_weight : contains necessary routines for
performing an
interpolation using a distance-weighted average of n nearest neighbors
used in scrip.F,
use
kinds_mod !
defines common data types
use constants
! defines common constants
use grids
! module containing grid info
use remap_vars
! module containing remap info
implicit none
!-----------------------------------------------------------------------
!
! module variables
!
!-----------------------------------------------------------------------
real (kind=dbl_kind), dimension(:),
allocatable, save ::
& coslat, sinlat,
! cosine, sine of grid lats (for distance)
& coslon, sinlon,
! cosine, sine of grid lons (for distance)
&
wgtstmp ! an array to
hold the link weight
!***********************************************************************
contains
!***********************************************************************
subroutine
remap_distwgt (lextrapdone, num_neighbors)
!-----------------------------------------------------------------------
!
! this routine computes the inverse-distance
weights for a
! nearest-neighbor interpolation.
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! input variables
!
!-----------------------------------------------------------------------
LOGICAL ::
&
lextrapdone ! logical, true if EXTRAP done on field
INTEGER (kind=int_kind) ::
&
num_neighbors ! number of neighbours
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
logical (kind=log_kind),
dimension(num_neighbors) ::
&
nbr_mask
! mask at nearest neighbors
logical (kind=log_kind) ::
ll_debug ! for debug outputs
integer (kind=int_kind) :: n,
&
dst_add, ! destination address
&
nmap, !
index of current map being computed
&
icount ! number
of non masked nearest neighbour
integer (kind=int_kind),
dimension(num_neighbors) ::
&
nbr_add ! source
address at nearest neighbors
real (kind=dbl_kind),
dimension(num_neighbors) ::
&
nbr_dist ! angular distance
four nearest neighbors
real (kind=dbl_kind) ::
&
coslat_dst, ! cos(lat) of destination grid point
&
coslon_dst, ! cos(lon) of destination grid point
&
sinlat_dst, ! sin(lat) of destination grid point
&
sinlon_dst, ! sin(lon) of destination grid point
&
dist_tot ! sum of neighbor
distances (for normalizing)
real (kind=dbl_kind) :: dl_test
integer (kind=int_kind) :: src_addnn,
il_debug_add
!-----------------------------------------------------------------------
!
IF (nlogprt .GE. 2) THEN
WRITE (UNIT =
nulou,FMT = *)' '
WRITE (UNIT =
nulou,FMT = *)'Entering routine remap_distwgt'
CALL FLUSH(nulou)
ENDIF
!
!-----------------------------------------------------------------------
!
! compute mappings from grid1 to grid2
!
!-----------------------------------------------------------------------
dl_test = 0.0
nmap = 1
!***
!*** allocate wgtstmp to be consistent
with store_link interface
!***
allocate (wgtstmp(num_wts))
!***
!*** compute cos, sin of lat/lon on
source grid for distance
!*** calculations
!***
allocate (coslat(grid1_size),
coslon(grid1_size),
&
sinlat(grid1_size), sinlon(grid1_size))
coslat = cos(grid1_center_lat)
coslon = cos(grid1_center_lon)
sinlat = sin(grid1_center_lat)
sinlon = sin(grid1_center_lon)
!***
!*** loop over destination grid
!***
grid_loop1: do dst_add = 1, grid2_size
if (.not.
grid2_mask(dst_add)) cycle grid_loop1
coslat_dst =
cos(grid2_center_lat(dst_add))
coslon_dst =
cos(grid2_center_lon(dst_add))
sinlat_dst =
sin(grid2_center_lat(dst_add))
sinlon_dst =
sin(grid2_center_lon(dst_add))
!***
!*** find nearest grid
points on source grid
!*** or non masked nearest
neighbour
!*** and distances to each
point
!***
call
grid_search_nbr(src_addnn,
&
nbr_add, nbr_dist,
&
grid2_center_lat(dst_add),
&
grid2_center_lon(dst_add),
&
coslat_dst, coslon_dst,
&
sinlat_dst, sinlon_dst,
&
bin_addr1, num_neighbors,lextrapdone,
&
dst_add)
ll_debug = .false.
IF (ll_debug) THEN
il_debug_add = 15248
IF
(dst_add == il_debug_add) THEN
WRITE(nulou,*) 'nbr_add = ', nbr_add(:)
WRITE(nulou,*) 'nbr_dist = ', nbr_dist(:)
CALL FLUSH(nulou)
ENDIF
ENDIF
!***
!*** compute weights based
on inverse distance
!*** if mask is false,
eliminate those points
!***
dist_tot = zero
do n=1,num_neighbors
if
(grid1_mask(nbr_add(n)) .or. lextrapdone) then
nbr_dist(n) = one/(nbr_dist(n) + epsilon(dl_test))
dist_tot = dist_tot + nbr_dist(n)
nbr_mask(n) = .true.
else
nbr_mask(n) = .false.
endif
end DO
IF (ll_debug) THEN
IF
(dst_add == il_debug_add) THEN
WRITE(nulou,*) 'nbr_mask = ', nbr_mask(:)
WRITE(nulou,*) 'nbr_dist = ', nbr_dist(:)
WRITE(nulou,*) 'dist_tot = ', dist_tot
endif
ENDIF
!***
!*** normalize weights and
store the link
!***
icount = 0
do n=1,num_neighbors
if (nbr_mask(n))
then
wgtstmp(1) = nbr_dist(n)/dist_tot
IF
(ll_debug) THEN
IF (dst_add == il_debug_add) THEN
WRITE(nulou,*) 'wgtstmp = ', wgtstmp(1)
WRITE(nulou,*) 'nbr_dist = ', nbr_dist(:)
ENDIF
ENDIF
call
store_link_nbr(nbr_add(n), dst_add, wgtstmp, nmap)
grid2_frac(dst_add) = one
icount = icount + 1
endif
end do
#ifndef NOT_NNEIGHBOUR
if (icount == 0) THEN
IF
(nlogprt .ge. 2) then
WRITE(nulou,*) ' '
WRITE(nulou,*)
& 'Using the nearest non-masked neighbour
because all 4
& surrounding source points are
masked for target point ',dst_add
CALL FLUSH(nulou)
ENDIF
wgtstmp(1) = 1.
grid2_frac(dst_add) = one
call
store_link_nbr(src_addnn, dst_add, wgtstmp, nmap)
ENDIF
#endif
end do grid_loop1
deallocate (coslat, coslon, sinlat,
sinlon)
deallocate(wgtstmp)
!
IF (nlogprt .GE. 2) THEN
WRITE (UNIT =
nulou,FMT = *)' '
WRITE (UNIT =
nulou,FMT = *)'Leaving routine remap_distwgt'
CALL FLUSH(nulou)
ENDIF
!
!-----------------------------------------------------------------------
end subroutine remap_distwgt
!***********************************************************************
subroutine grid_search_nbr(src_addnn,
&
nbr_add, nbr_dist, plat, plon,
&
coslat_dst, coslon_dst, sinlat_dst, sinlon_dst,
&
src_bin_add, num_neighbors, lextrapdone, dst_add)
!-----------------------------------------------------------------------
!
! this routine finds the closest num_neighbor
points to a search
! point and computes a distance to each of the
neighbors.
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! input variables
!
!-----------------------------------------------------------------------
integer (kind=int_kind), dimension(:,:),
intent(in) ::
&
src_bin_add
! search bins for restricting search
integer (kind=int_kind), intent(in) ::
dst_add
real (kind=dbl_kind), intent(in) ::
&
plat, ! latitude
of the search point
&
plon, ! longitude of
the search point
& coslat_dst,
! cos(lat) of the search point
& coslon_dst,
! cos(lon) of the search point
& sinlat_dst,
! sin(lat) of the search point
&
sinlon_dst ! sin(lon) of the search point
INTEGER (kind=int_kind) ::
&
num_neighbors ! number of neighbours
LOGICAL ::
&
lextrapdone, ! logical, true if EXTRAP done on field
&
ll_debug
!-----------------------------------------------------------------------
!
! output variables
!
!-----------------------------------------------------------------------
integer (kind=int_kind),
dimension(num_neighbors), intent(out) ::
& nbr_add ! address
of each of the closest points
real (kind=dbl_kind),
dimension(num_neighbors), intent(out) ::
& nbr_dist ! distance to
each of the closest points
integer (kind=int_kind), intent(out) ::
src_addnn
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (kind=int_kind) :: n, nmax,
nadd, nchk, ! dummy indices
& min_add, max_add, nm1,
np1, i, j, ip1, im1, jp1, jm1,
& il_debug_add
real (kind=dbl_kind) ::
& distance, arg,
src_latsnn ! angular distance
!-----------------------------------------------------------------------
!
! loop over source grid and find nearest
neighbors
!
!-----------------------------------------------------------------------
!***
!*** restrict the search using search
bins
!*** expand the bins to catch neighbors
!***
select case (restrict_type)
case('LATITUDE')
do n=1,num_srch_bins
if (plat >=
bin_lats(1,n) .and. plat <= bin_lats(2,n)) then
min_add = src_bin_add(1,n)
max_add = src_bin_add(2,n)
nm1
= max(n-1,1)
np1
= min(n+1,num_srch_bins)
min_add = min(min_add,src_bin_add(1,nm1))
max_add = max(max_add,src_bin_add(2,nm1))
min_add = min(min_add,src_bin_add(1,np1))
max_add = max(max_add,src_bin_add(2,np1))
endif
end do
case('LATLON')
n = 0
nmax =
nint(sqrt(real(num_srch_bins)))
do j=1,nmax
jp1 = min(j+1,nmax)
jm1 = max(j-1,1)
do i=1,nmax
ip1 =
min(i+1,nmax)
im1 = max(i-1,1)
n = n+1
if (plat >=
bin_lats(1,n) .and. plat <= bin_lats(2,n) .and.
& plon >=
bin_lons(1,n) .and. plon <= bin_lons(2,n)) then
min_add = src_bin_add(1,n)
max_add = src_bin_add(2,n)
nm1
= (jm1-1)*nmax + im1
np1
= (jp1-1)*nmax + ip1
nm1
= max(nm1,1)
np1
= min(np1,num_srch_bins)
min_add = min(min_add,src_bin_add(1,nm1))
max_add = max(max_add,src_bin_add(2,nm1))
min_add = min(min_add,src_bin_add(1,np1))
max_add = max(max_add,src_bin_add(2,np1))
endif
end do
end do
end select
!***
!*** initialize distance and address
arrays
!***
nbr_add = 0
nbr_dist = bignum
src_latsnn = bignum
src_addnn = 0
do nadd=min_add,max_add
!***
!*** find distance to this
point
!***
arg =
sinlat_dst*sinlat(nadd) +
& coslat_dst*coslat(nadd)*
&
(coslon_dst*coslon(nadd) + sinlon_dst*sinlon(nadd))
IF (arg < -1.0d0) THEN
arg
= -1.0d0
ELSE IF (arg > 1.0d0) THEN
arg
= 1.0d0
END IF
distance = acos(arg)
ll_debug = .false.
IF (ll_debug) THEN
il_debug_add = 15248
IF
(dst_add == il_debug_add) THEN
WRITE(nulou,1009) nadd, distance
ENDIF
ENDIF
!***
!*** store the address and
distance if this is one of the
!*** smallest four so far
!***
check_loop: do
nchk=1,num_neighbors
if (distance
.lt. nbr_dist(nchk)) then
do
n=num_neighbors,nchk+1,-1
nbr_add(n) = nbr_add(n-1)
nbr_dist(n) = nbr_dist(n-1)
end
do
nbr_add(nchk) = nadd
nbr_dist(nchk) = distance
exit
check_loop
endif
end do check_loop
IF (ll_debug) THEN
IF
(dst_add == il_debug_add) THEN
WRITE(nulou,1010) nadd, distance
ENDIF
ENDIF
#ifndef NOT_NNEIGHBOUR
!***
!*** store the non-masked
closest neighbour
!***
if (distance <
src_latsnn) THEN
if
(grid1_mask(nadd) .or. lextrapdone) THEN
src_addnn = nadd
src_latsnn = distance
ENDIF
IF
(ll_debug) THEN
IF (dst_add == il_debug_add) THEN
WRITE(nulou,1010) nadd, distance
ENDIF
ENDIF
ENDIF
#endif
end DO
1009 FORMAT (1X, 'nadd0 =', 1X, I6, 2X, 'distance0 =', 1X, F18.16)
1010 FORMAT (1X, 'nadd1 =', 1X, I6, 2X, 'distance0 =', 1X, F18.16)
1011 FORMAT (1X, 'nadd2 =', 1X, I6, 2X, 'distance2 =', 1X, F18.16)
!-----------------------------------------------------------------------
end subroutine grid_search_nbr
!***********************************************************************
subroutine store_link_nbr(add1, add2,
weights, nmap)
!-----------------------------------------------------------------------
!
! this routine stores the address and weight
for this link in
! the appropriate address and weight arrays and
resizes those
! arrays if necessary.
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! input variables
!
!-----------------------------------------------------------------------
integer (kind=int_kind), intent(in) ::
& add1, ! address
on grid1
& add2, ! address
on grid2
& nmap !
identifies which direction for mapping
real (kind=dbl_kind), dimension(:),
intent(in) ::
& weights ! array of
remapping weights for this link
!-----------------------------------------------------------------------
!
! increment number of links and check to see if
remap arrays need
! to be increased to accomodate the new
link. then store the
! link.
!
!-----------------------------------------------------------------------
select case (nmap)
case(1)
num_links_map1 =
num_links_map1 + 1
if (num_links_map1 >
max_links_map1)
& call
resize_remap_vars(1,resize_increment)
grid1_add_map1(num_links_map1) = add1
grid2_add_map1(num_links_map1) = add2
wts_map1
(:,num_links_map1) = weights
case(2)
num_links_map2 =
num_links_map2 + 1
if (num_links_map2 >
max_links_map2)
& call
resize_remap_vars(2,resize_increment)
grid1_add_map2(num_links_map2) = add1
grid2_add_map2(num_links_map2) = add2
wts_map2
(:,num_links_map2) = weights
end select
!-----------------------------------------------------------------------
end subroutine store_link_nbr