Back to OASIS3 home
Modules used for the scrip
library (in oasis3/lib/scrip/src) :
remap_bicubic : contains necessary routines for
performing a bicubic interpolation
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
!-----------------------------------------------------------------------
integer (kind=int_kind), parameter ::
& max_iter =
100 ! max iteration count for i,j iteration
real (kind=dbl_kind), parameter ::
& converge =
epsilon(1.0_dbl_kind) ! convergence criterion
!***********************************************************************
contains
!***********************************************************************
subroutine
remap_bicub (lextrapdone)
!-----------------------------------------------------------------------
!
! this routine computes the weights for a
bicubic interpolation.
!
!-----------------------------------------------------------------------
LOGICAL :: lextrapdone !
logical, true if EXTRAP done on field
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (kind=int_kind) :: n,icount,i,
&
dst_add, ! destination address
&
iter, !
iteration counter
&
nmap
! index of current map being computed
integer (kind=int_kind), dimension(4) ::
&
src_add ! address for
the four source points
real (kind=dbl_kind), dimension(4)
::
&
src_lats, ! latitudes of
four corners
&
src_lons ! longitudes of
four corners
real (kind=dbl_kind),
dimension(4,4) ::
&
wgts
! bicubic weights for four corners
real (kind=dbl_kind) ::
& plat,
plon, ! lat/lon coords of
destination point
& iguess,
jguess, ! current guess for bicubic coordinate
& deli,
delj, ! corrections to i,j
& dth1, dth2,
dth3, ! some latitude differences
& dph1, dph2,
dph3, ! some longitude differences
& dthp,
dphp, ! difference between point
and sw corner
& mat1, mat2,
mat3, mat4, ! matrix elements
&
determinant, ! matrix determinant
&
sum_wgts ! sum of
weights for normalization
real (kind=dbl_kind) ::
&
coslat_dst, sinlat_dst, coslon_dst, sinlon_dst,
& dist_min,
distance, ! for computing dist-weighted avg
&
src_latsnn, arg
integer (kind=int_kind) :: min_add,
max_add, srch_add, src_addnn
!-----------------------------------------------------------------------
!
! compute mappings from grid1 to grid2
!
!-----------------------------------------------------------------------
!
IF (nlogprt .GE. 2) THEN
WRITE (UNIT =
nulou,FMT = *)' '
WRITE (UNIT =
nulou,FMT = *)'Entering routine remap_bicub'
CALL FLUSH(nulou)
ENDIF
!
nmap = 1
if (grid1_rank /= 2) then
stop 'Can not do bicubic
interpolation when grid_rank /= 2'
endif
!***
!*** loop over destination grid
!***
grid_loop1: do dst_add = 1, grid2_size
if (.not.
grid2_mask(dst_add)) cycle grid_loop1
plat =
grid2_center_lat(dst_add)
plon =
grid2_center_lon(dst_add)
!***
!*** find nearest square of
grid points on source grid
!***
call
grid_search_bicub(src_add, src_lats, src_lons,
&
min_add, max_add,
&
plat, plon, grid1_dims,
&
grid1_center_lat, grid1_center_lon,
&
grid1_bound_box, bin_addr1,
&
lextrapdone)
if (src_add(1) > 0) then
!***
!*** if the 4
surrounding points have been found and are
!*** non-masked,
do bicubic interpolation
!***
grid2_frac(dst_add) = one
dth1 =
src_lats(2) - src_lats(1)
dth2 =
src_lats(4) - src_lats(1)
dth3 =
src_lats(3) - src_lats(2) - dth2
dph1 =
src_lons(2) - src_lons(1)
dph2 =
src_lons(4) - src_lons(1)
dph3 =
src_lons(3) - src_lons(2)
if (dph1
> three*pih) dph1 = dph1 - pi2
if (dph2
> three*pih) dph2 = dph2 - pi2
if (dph3
> three*pih) dph3 = dph3 - pi2
if (dph1 <
-three*pih) dph1 = dph1 + pi2
if (dph2 <
-three*pih) dph2 = dph2 + pi2
if (dph3 <
-three*pih) dph3 = dph3 + pi2
dph3 = dph3 -
dph2
iguess = half
jguess = half
iter_loop1: do
iter=1,max_iter
dthp
= plat - src_lats(1) - dth1*iguess -
&
dth2*jguess - dth3*iguess*jguess
dphp
= plon - src_lons(1)
if
(dphp > three*pih) dphp = dphp - pi2
if
(dphp < -three*pih) dphp = dphp + pi2
dphp
= dphp - dph1*iguess - dph2*jguess -
&
dph3*iguess*jguess
mat1
= dth1 + dth3*jguess
mat2
= dth2 + dth3*iguess
mat3
= dph1 + dph3*jguess
mat4
= dph2 + dph3*iguess
determinant = mat1*mat4 - mat2*mat3
deli
= (dthp*mat4 - mat2*dphp)/determinant
delj
= (mat1*dphp - dthp*mat3)/determinant
if
(abs(deli) < converge .and.
& abs(delj)
< converge) exit iter_loop1
iguess = iguess + deli
jguess = jguess + delj
end do iter_loop1
if (iter <=
max_iter) then
!***
!***
successfully found i,j - compute weights
!***
wgts(1,1) = (one - jguess**2*(three-two*jguess))*
&
(one - iguess**2*(three-two*iguess))
wgts(1,2) = (one - jguess**2*(three-two*jguess))*
&
iguess**2*(three-two*iguess)
wgts(1,3) =
jguess**2*(three-two*jguess)*
&
iguess**2*(three-two*iguess)
wgts(1,4) =
jguess**2*(three-two*jguess)*
&
(one - iguess**2*(three-two*iguess))
wgts(2,1) = (one - jguess**2*(three-two*jguess))*
&
iguess*(iguess-one)**2
wgts(2,2) = (one - jguess**2*(three-two*jguess))*
&
iguess**2*(iguess-one)
wgts(2,3) =
jguess**2*(three-two*jguess)*
&
iguess**2*(iguess-one)
wgts(2,4) =
jguess**2*(three-two*jguess)*
&
iguess*(iguess-one)**2
wgts(3,1) =
jguess*(jguess-one)**2*
&
(one - iguess**2*(three-two*iguess))
wgts(3,2) =
jguess*(jguess-one)**2*
&
iguess**2*(three-two*iguess)
wgts(3,3) =
jguess**2*(jguess-one)*
&
iguess**2*(three-two*iguess)
wgts(3,4) =
jguess**2*(jguess-one)*
&
(one - iguess**2*(three-two*iguess))
wgts(4,1) =
iguess*(iguess-one)**2*
&
jguess*(jguess-one)**2
wgts(4,2) =
iguess**2*(iguess-one)*
&
jguess*(jguess-one)**2
wgts(4,3) =
iguess**2*(iguess-one)*
&
jguess**2*(jguess-one)
wgts(4,4) =
iguess*(iguess-one)**2*
&
jguess**2*(jguess-one)
call
store_link_bicub(dst_add, src_add, wgts, nmap)
ELSE
WRITE (UNIT = nulou,FMT = *)
&
'Iteration for i,j exceed max iteration count'
STOP
endif
else if (src_add(1) < 0)
THEN
!***
!*** Search for
bicubic failed or at least one of the 4
!*** neighbours
was masked. Do distance-weighted average using
!*** the
non-masked points among the 4 closest ones.
!***
IF (nlogprt .eq.
2) then
WRITE(nulou,*) ' '
WRITE(nulou,*)
& 'WARNING: Cannot make bicubic
interpolation for target point'
&
,dst_add
WRITE(nulou,*)
& 'Using non-masked
points among 4 nearest neighbors.'
WRITE(nulou,*) ' '
ENDIF
!***
!*** Find the 4
closest points
!***
coslat_dst =
cos(plat)
sinlat_dst =
sin(plat)
coslon_dst =
cos(plon)
sinlon_dst =
sin(plon)
src_add = 0
dist_min = bignum
src_lats = bignum
do srch_add =
min_add,max_add
arg
= coslat_dst*cos(grid1_center_lat(srch_add))*
&
(coslon_dst*cos(grid1_center_lon(srch_add)) +
&
sinlon_dst*sin(grid1_center_lon(srch_add)))+
&
sinlat_dst*sin(grid1_center_lat(srch_add))
IF
(arg < -1.0d0) THEN
arg = -1.0d0
ELSE
IF (arg > 1.0d0) THEN
arg = 1.0d0
END
IF
distance=acos(arg)
if
(distance < dist_min) then
sort_loop: do n=1,4
if (distance < src_lats(n)) then
do i=4,n+1,-1
src_add (i) = src_add (i-1)
src_lats(i) = src_lats(i-1)
end do
src_add (n) = srch_add
src_lats(n) = distance
dist_min = src_lats(4)
exit sort_loop
endif
end do sort_loop
endif
end do
src_lons =
one/(src_lats + tiny)
distance =
sum(src_lons)
src_lats =
src_lons/distance
!***
!*** Among 4
closest points, keep only the non-masked ones
!***
icount = 0
do n=1,4
if
(grid1_mask(src_add(n)) .or.
& (.not.
grid1_mask(src_add(n)) .and. lextrapdone)) then
icount = icount + 1
else
src_lats(n) = zero
endif
end do
if (icount >
0) then
!*** renormalize weights
sum_wgts = sum(src_lats)
wgts(1,1) = src_lats(1)/sum_wgts
wgts(1,2) = src_lats(2)/sum_wgts
wgts(1,3) = src_lats(3)/sum_wgts
wgts(1,4) = src_lats(4)/sum_wgts
wgts(2,:) = 0.
wgts(3,:) = 0.
wgts(4,:) = 0.
grid2_frac(dst_add) = one
call store_link_bicub(dst_add, src_add, wgts, nmap)
ELSE
IF (NLOGPRT .GE. 2) THEN
WRITE(nulou,*) ' '
WRITE(nulou,*)' Using the nearest non-masked neighbour
& as all 4 surrounding source points are
masked for target point',
&
dst_add
WRITE(nulou,*) 'with longitude and latitude',
&
plon, plat
ENDIF
src_latsnn = bignum
!cdir novector
do srch_add = min_add,max_add
if (grid1_mask(srch_add) .or.
& (.not.
grid1_mask(srch_add) .and. lextrapdone)) THEN
arg = coslat_dst*cos(grid1_center_lat(srch_add))*
&
(coslon_dst*cos(grid1_center_lon(srch_add)) +
&
sinlon_dst*sin(grid1_center_lon(srch_add)))+
&
sinlat_dst*sin(grid1_center_lat(srch_add))
IF (arg < -1.0d0) THEN
arg = -1.0d0
ELSE IF (arg > 1.0d0) THEN
arg = 1.0d0
END
IF
distance=acos(arg)
if (distance < src_latsnn) then
src_addnn = srch_add
src_latsnn = distance
endif
endif
end DO
IF (NLOGPRT .GE. 2) THEN
WRITE(nulou,*) ' '
WRITE(nulou,*)
&
'Nearest non masked neighbour is source point ',
&
src_addnn
WRITE(nulou,*) 'with longitude and latitude',
&
grid1_center_lon(src_addnn), grid1_center_lat(src_addnn)
ENDIF
!
wgts(1,1) = 1.
wgts(1,2:4) = 0.
wgts(2,:) = 0.
wgts(3,:) = 0.
wgts(4,:) = 0.
src_add(1) = src_addnn
src_add(2) = 0
src_add(3) = 0
src_add(4) = 0
grid2_frac(dst_add) = one
call store_link_bicub(dst_add, src_add, wgts, nmap)
ENDIF
ENDIF
end do grid_loop1
!
IF (nlogprt .GE. 2) THEN
WRITE (UNIT =
nulou,FMT = *)' '
WRITE (UNIT =
nulou,FMT = *)'Leaving routine remap_bicub'
CALL FLUSH(nulou)
ENDIF
!
end subroutine remap_bicub
!***********************************************************************
subroutine grid_search_bicub(src_add,
src_lats, src_lons,
&
min_add, max_add,
&
plat, plon, src_grid_dims,
&
src_center_lat, src_center_lon,
&
src_bound_box,
&
src_bin_add, lextrapdone)
!-----------------------------------------------------------------------
!
! this routine finds the location of the search
point plat, plon
! in the source grid and returns the corners
needed for a bicubic
! interpolation.
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! output variables
!
!-----------------------------------------------------------------------
integer (kind=int_kind), dimension(4),
intent(out) ::
& src_add ! address
of each corner point enclosing P
real (kind=dbl_kind), dimension(4),
intent(out) ::
& src_lats, !
latitudes of the four corner points
& src_lons !
longitudes of the four corner points
integer (kind=int_kind) :: min_add,
max_add
!-----------------------------------------------------------------------
!
! input variables
!
!-----------------------------------------------------------------------
real (kind=dbl_kind), intent(in) ::
& plat, !
latitude of the search point
& plon
! longitude of the search point
integer (kind=int_kind), dimension(2),
intent(in) ::
& src_grid_dims !
size of each src grid dimension
real (kind=dbl_kind), dimension(:),
intent(in) ::
& src_center_lat, !
latitude of each src grid center
& src_center_lon !
longitude of each src grid center
real (kind=dbl_kind), dimension(:,:),
intent(in) ::
&
src_bound_box ! bounding box for src grid search
integer (kind=int_kind), dimension(:,:),
intent(in) ::
&
src_bin_add ! search bins for restricting
LOGICAL ::
&
lextrapdone ! logical, true if EXTRAP done on field
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (kind=int_kind) :: n, next_n,
srch_add, ni, ! dummy indices
& nx, ny,
ntotmask, !
dimensions of src grid
& i, j, jp1, ip1, n_add,
e_add, ne_add ! addresses
real (kind=dbl_kind) :: ! vectors
for cross-product check
& vec1_lat,
vec1_lon,
& vec2_lat,
vec2_lon, cross_product, cross_product_last
!-----------------------------------------------------------------------
!
! restrict search first using search bins.
!
!-----------------------------------------------------------------------
src_add = 0
min_add = size(src_center_lat)
max_add = 1
do n=1,num_srch_bins
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 =
min(min_add, src_bin_add(1,n))
max_add =
max(max_add, src_bin_add(2,n))
endif
end do
!-----------------------------------------------------------------------
!
! now perform a more detailed search
!
!-----------------------------------------------------------------------
nx = src_grid_dims(1)
ny = src_grid_dims(2)
srch_loop: do srch_add = min_add,max_add
!*** first check bounding box
if (plat <=
src_bound_box(2,srch_add) .and.
& plat >=
src_bound_box(1,srch_add) .and.
& plon <=
src_bound_box(4,srch_add) .and.
& plon >=
src_bound_box(3,srch_add)) then
!***
!*** we are
within bounding box so get really serious
!***
!*** find N,S
and NE points to this grid point
j = (srch_add -
1)/nx +1
i = srch_add -
(j-1)*nx
!*** find ip1
!*** Note: I do
not want to take into account the number
!*** of
overlapping grid points, as the underlying cell
!*** will be
found in all cases if the grid overlaps.
if (i < nx)
then
ip1 = i + 1
else
ip1 = 1
endif
if (j < ny)
then
jp1 = j+1
else
jp1 = 1
endif
n_add = (jp1 -
1)*nx + i
e_add = (j -
1)*nx + ip1
ne_add = (jp1 -
1)*nx + ip1
src_lats(1) =
src_center_lat(srch_add)
src_lats(2) =
src_center_lat(e_add)
src_lats(3) =
src_center_lat(ne_add)
src_lats(4) =
src_center_lat(n_add)
src_lons(1) =
src_center_lon(srch_add)
src_lons(2) =
src_center_lon(e_add)
src_lons(3) =
src_center_lon(ne_add)
src_lons(4) =
src_center_lon(n_add)
!***
!*** for
consistency, we must make sure all lons are in
!*** same 2pi
interval
!***
vec1_lon =
src_lons(1) - plon
if (vec1_lon
> pi) then
src_lons(1) = src_lons(1) - pi2
else if
(vec1_lon < -pi) then
src_lons(1) = src_lons(1) + pi2
endif
do n=2,4
vec1_lon = src_lons(n) - src_lons(1)
if
(vec1_lon > pi) then
src_lons(n) = src_lons(n) - pi2
else
if (vec1_lon < -pi) then
src_lons(n) = src_lons(n) + pi2
endif
end do
corner_loop: do
n=1,4
next_n = MOD(n,4) + 1
!***
!***
here we take the cross product of the vector making
!***
up each box side with the vector formed by the vertex
!***
and search point. if all the cross products are
!***
same sign, the point is contained in the box.
!***
vec1_lat = src_lats(next_n) - src_lats(n)
vec1_lon = src_lons(next_n) - src_lons(n)
vec2_lat = plat - src_lats(n)
vec2_lon = plon - src_lons(n)
!***
!***
check for 0,2pi crossings
!***
if
(vec1_lon > three*pih) then
vec1_lon = vec1_lon - pi2
else
if (vec1_lon < -three*pih) then
vec1_lon = vec1_lon + pi2
endif
if
(vec2_lon > three*pih) then
vec2_lon = vec2_lon - pi2
else
if (vec2_lon < -three*pih) then
vec2_lon = vec2_lon + pi2
endif
cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat
!***
!***
if cross product is less than zero, this cell
!***
doesn't work
!***
if
(n == 1) cross_product_last = cross_product
if
(cross_product*cross_product_last < zero) then
exit corner_loop
else
cross_product_last = cross_product
endif
end do
corner_loop
!***
!*** if cross
products all positive, we found the location
!***
if (n > 4)
then
src_add(1) = srch_add
src_add(2) = e_add
src_add(3) = ne_add
src_add(4) = n_add
! Check if one point is masked; IF so, nearest-neighbour
! interpolation will be used
ntotmask = 0
do ni=1,4
if (.not. grid1_mask(src_add(ni)).and.
&
.not. lextrapdone)
&
ntotmask = ntotmask + 1
end DO
IF (ntotmask .gt. 0) src_add(1) = -src_add(1)
return
endif
!***
!*** otherwise
move on to next cell
!***
endif
!bounding box check
end do srch_loop
!***
!*** if no cell found, point is likely
either in a box that straddles
!*** either pole or is outside the grid.
Put src_add = -1 so that
!*** distance-weighted average of the 4
non-masked closest points
!*** is done in calling routine.
src_add = -1
!-----------------------------------------------------------------------
end subroutine grid_search_bicub
!***********************************************************************
subroutine store_link_bicub(dst_add,
src_add, weights, nmap)
!-----------------------------------------------------------------------
!
! this routine stores the address and weight
for four links
! associated with one destination point in the
appropriate address
! and weight arrays and resizes those arrays if
necessary.
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! input variables
!
!-----------------------------------------------------------------------
integer (kind=int_kind), intent(in) ::
& dst_add, !
address on destination grid
&
nmap ! identifies which direction for
mapping
integer (kind=int_kind), dimension(4),
intent(in) ::
& src_add !
addresses on source grid
real (kind=dbl_kind), dimension(4,4),
intent(in) ::
& weights ! array of
remapping weights for these links
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (kind=int_kind) :: n, ! dummy
index
&
num_links_old !
placeholder for old link number
!-----------------------------------------------------------------------
!
! 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_old =
num_links_map1
num_links_map1 =
num_links_old + 4
if (num_links_map1 >
max_links_map1)
& call
resize_remap_vars(1,resize_increment)
do n=1,4
grid1_add_map1(num_links_old+n) = src_add(n)
grid2_add_map1(num_links_old+n) = dst_add
wts_map1 (:,num_links_old+n) = weights(:,n)
end do
case(2)
num_links_old =
num_links_map2
num_links_map2 =
num_links_old + 4
if (num_links_map2 >
max_links_map2)
& call
resize_remap_vars(2,resize_increment)
do n=1,4
grid1_add_map2(num_links_old+n) = dst_add
grid2_add_map2(num_links_old+n) = src_add(n)
wts_map2 (:,num_links_old+n) = weights(:,n)
end do
end select
!-----------------------------------------------------------------------
end subroutine store_link_bicub