#include <misc.h>
#include <params.h>
subroutine bilin (arrin, xin, yin, nlondin, nlonin, & 3,13
nlevdin, nlev, nlatin, arrout, xout, &
yout, nlondout, nlonout, nlevdout, nlatout)
!-----------------------------------------------------------------------
!
! Purpose:
!
! Do a bilinear interpolation from input mesh defined by xin, yin to output
! mesh defined by xout, yout. Circularity is assumed in the x-direction so
! input x-grid must be in degrees east and must start from Greenwich. When
! extrapolation is necessary in the N-S direction, values will be copied
! from the extreme edge of the input grid. Vectorization is over the
! longitude dimension. Input grid is assumed rectangular. Output grid
! is assumed ragged, i.e. xout is a 2-d mesh.
!
! Method: Interpolate first in longitude, then in latitude.
!
! Author: Jim Rosinski
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use abortutils, only: endrun
!-----------------------------------------------------------------------
implicit none
!-----------------------------------------------------------------------
!
! Input arguments
!
integer, intent(in) :: nlondin ! longitude dimension of input grid
integer, intent(in) :: nlonin ! number of real longitudes (input)
integer, intent(in) :: nlevdin ! vertical dimension of input grid
integer, intent(in) :: nlev ! number of vertical levels
integer, intent(in) :: nlatin ! number of input latitudes
integer, intent(in) :: nlatout ! number of output latitudes
integer, intent(in) :: nlondout ! longitude dimension of output grid
integer, intent(in) :: nlonout(nlatout) ! number of output longitudes per lat
integer, intent(in) :: nlevdout ! vertical dimension of output grid
real(r8), intent(in) :: arrin(nlondin,nlevdin,nlatin) ! input array of values to interpolate
real(r8), intent(in) :: xin(nlondin) ! input x mesh
real(r8), intent(in) :: yin(nlatin) ! input y mesh
real(r8), intent(in) :: xout(nlondout,nlatout) ! output x mesh
real(r8), intent(in) :: yout(nlatout) ! output y mesh
!
! Output arguments
!
real(r8), intent(out) :: arrout(nlondout,nlevdout,nlatout) ! interpolated array
!
! Local workspace
!
integer :: i, ii, iw, ie, iiprev ! longitude indices
integer :: j, jj, js, jn, jjprev ! latitude indices
integer :: k ! level index
integer :: icount ! number of bad values
real(r8) :: extrap ! percent grid non-overlap
real(r8) :: dxinwrap ! delta-x on input grid for 2-pi
real(r8) :: avgdxin ! avg input delta-x
real(r8) :: ratio ! compare dxinwrap to avgdxin
real(r8) :: sum ! sum of weights (used for testing)
!
! Dynamic
!
integer :: iim(nlondout) ! interpolation index to the left
integer :: iip(nlondout) ! interpolation index to the right
integer :: jjm(nlatout) ! interpolation index to the south
integer :: jjp(nlatout) ! interpolation index to the north
real(r8) :: wgts(nlatout) ! interpolation weight to the north
real(r8) :: wgtn(nlatout) ! interpolation weight to the north
real(r8) :: wgte(nlondout) ! interpolation weight to the north
real(r8) :: wgtw(nlondout) ! interpolation weight to the north
real(r8) :: igrid(nlonin) ! interpolation weight to the north
!
! Check validity of input coordinate arrays: must be monotonically increasing,
! and have a total of at least 2 elements
!
if (nlonin < 2 .or. nlatin < 2) then
call endrun ('BILIN: Must have at least 2 input points for interpolation')
end if
if (xin(1) < 0. .or. xin(nlonin) > 360.) then
call endrun ('BILIN: Input x-grid must be between 0 and 360')
end if
icount = 0
do i=1,nlonin-1
if (xin(i) >= xin(i+1)) icount = icount + 1
end do
do j=1,nlatin-1
if (yin(j) >= yin(j+1)) icount = icount + 1
end do
do j=1,nlatout-1
if (yout(j) >= yout(j+1)) icount = icount + 1
end do
do j=1,nlatout
do i=1,nlonout(j)-1
if (xout(i,j) >= xout(i+1,j)) icount = icount + 1
end do
end do
if (icount > 0) then
call endrun ('BILIN: Non-monotonic coordinate array(s) found')
end if
if (yout(nlatout) <= yin(1) .or. yout(1) >= yin(nlatin)) then
call endrun ('BILIN: No overlap in y-coordinates')
end if
#if ( ! defined SCAM )
do j=1,nlatout
if (xout(1,j) < 0. .or. xout(nlonout(j),j) > 360.) then
call endrun ('BILIN: Output x-grid must be between 0 and 360')
end if
if (xout(nlonout(j),j) <= xin(1) .or. &
xout(1,j) >= xin(nlonin)) then
call endrun ('BILIN: No overlap in x-coordinates')
end if
end do
#endif
!
! Initialize index arrays for later checking
!
do j=1,nlatout
jjm(j) = 0
jjp(j) = 0
end do
!
! For values which extend beyond N and S boundaries, set weights
! such that values will just be copied.
!
do js=1,nlatout
if (yout(js) > yin(1)) exit
jjm(js) = 1
jjp(js) = 1
wgts(js) = 1.
wgtn(js) = 0.
end do
do jn=nlatout,1,-1
if (yout(jn) <= yin(nlatin)) exit
jjm(jn) = nlatin
jjp(jn) = nlatin
wgts(jn) = 1.
wgtn(jn) = 0.
end do
!
! Loop though output indices finding input indices and weights
!
jjprev = 1
do j=js,jn
do jj=jjprev,nlatin-1
if (yout(j) > yin(jj) .and. yout(j) <= yin(jj+1)) then
jjm(j) = jj
jjp(j) = jj + 1
wgts(j) = (yin(jj+1) - yout(j)) / (yin(jj+1) - yin(jj))
wgtn(j) = (yout(j) - yin(jj)) / (yin(jj+1) - yin(jj))
goto 30
end if
end do
call endrun ('BILIN: Failed to find interp values')
30 jjprev = jj
end do
dxinwrap = xin(1) + 360. - xin(nlonin)
!
! Check for sane dxinwrap values. Allow to differ no more than 10% from avg
!
avgdxin = (xin(nlonin)-xin(1))/(nlonin-1.)
ratio = dxinwrap/avgdxin
if (ratio < 0.9 .or. ratio > 1.1) then
write(6,*)'BILIN: Insane dxinwrap value =',dxinwrap,' avg=', avgdxin
call endrun
end if
!
! Check grid overlap
!
extrap = 100.*((js - 1.) + float(nlatout - jn))/nlatout
if (extrap > 20.) then
write(6,*)'BILIN:',extrap,' % of N/S output grid will have to be extrapolated'
end if
!
! Check that interp/extrap points have been found for all outputs, and that
! they are reasonable (i.e. within 32-bit roundoff)
!
icount = 0
do j=1,nlatout
if (jjm(j) == 0 .or. jjp(j) == 0) icount = icount + 1
sum = wgts(j) + wgtn(j)
if (sum < 0.99999 .or. sum > 1.00001) icount = icount + 1
if (wgts(j) < 0. .or. wgts(j) > 1.) icount = icount + 1
if (wgtn(j) < 0. .or. wgtn(j) > 1.) icount = icount + 1
end do
if (icount > 0) then
call endrun ('BILIN: something bad in latitude indices or weights')
end if
!
! Do the bilinear interpolation
!
do j=1,nlatout
!
! Initialize index arrays for later checking
!
do i=1,nlondout
iim(i) = 0
iip(i) = 0
end do
!
! For values which extend beyond E and W boundaries, set weights such that
! values will be interpolated between E and W edges of input grid.
!
do iw=1,nlonout(j)
if (xout(iw,j) > xin(1)) exit
iim(iw) = nlonin
iip(iw) = 1
wgtw(iw) = (xin(1) - xout(iw,j)) /dxinwrap
wgte(iw) = (xout(iw,j)+360. - xin(nlonin))/dxinwrap
end do
do ie=nlonout(j),1,-1
if (xout(ie,j) <= xin(nlonin)) exit
iim(ie) = nlonin
iip(ie) = 1
wgtw(ie) = (xin(1)+360. - xout(ie,j)) /dxinwrap
wgte(ie) = (xout(ie,j) - xin(nlonin))/dxinwrap
end do
!
! Loop though output indices finding input indices and weights
!
iiprev = 1
do i=iw,ie
do ii=iiprev,nlonin-1
if (xout(i,j) > xin(ii) .and. xout(i,j) <= xin(ii+1)) then
iim(i) = ii
iip(i) = ii + 1
wgtw(i) = (xin(ii+1) - xout(i,j)) / (xin(ii+1) - xin(ii))
wgte(i) = (xout(i,j) - xin(ii)) / (xin(ii+1) - xin(ii))
goto 60
end if
end do
call endrun ('BILIN: Failed to find interp values')
60 iiprev = ii
end do
icount = 0
do i=1,nlonout(j)
if (iim(i) == 0 .or. iip(i) == 0) icount = icount + 1
sum = wgtw(i) + wgte(i)
if (sum < 0.99999 .or. sum > 1.00001) icount = icount + 1
if (wgtw(i) < 0. .or. wgtw(i) > 1.) icount = icount + 1
if (wgte(i) < 0. .or. wgte(i) > 1.) icount = icount + 1
end do
if (icount > 0) then
write(6,*)'BILIN: j=',j,' Something bad in longitude indices or weights'
call endrun
end if
!
! Do the interpolation, 1st in longitude then latitude
!
do k=1,nlev
do i=1,nlonin
igrid(i) = arrin(i,k,jjm(j))*wgts(j) + arrin(i,k,jjp(j))*wgtn(j)
end do
do i=1,nlonout(j)
arrout(i,k,j) = igrid(iim(i))*wgtw(i) + igrid(iip(i))*wgte(i)
end do
end do
end do
return
end subroutine bilin