module smooth_polcarf 4,1
!$$$ module documentation block
! . . . .
! module: smooth_polcarf smooth polar cascade interpolation stuff
! prgmmr: parrish org: np23 date: 2005-04-28
!
! abstract: contains everything needed to implement polar smoothing
! spline interpolation, as an alternative to the
! existing polcas polar cascade routines which have
! noise problems near the pole due to magnification
! of normally small errors in the longitude direction
! very close to the pole.
! B-splines are used in smoothing form, not as a direct
! fit.
!
! program history log:
! 2005-04-28 parrish
!
! subroutines included:
! sub smooth_polcas - interpolate from x-y to lat-lon using smooth cascade interpolation.
!
! Variable Definitions:
! def norsp - order of smooth cascade interpolation
! if =0, then old cascade interpolation is used
! def nlon_m89 - 1st longitude index after leaving east side of exclusion zone around y < 0
! def nlon_p89 - last longitude index before entering west side of exclusion zone around y > 0
! def nlon_p91 - 1st longitude index after leaving east side of exclusion zone around y > 0
! def nlon_p269 - last longitude index before entering west side of exclusion zone around y < 0
! def nlon_p1 - 1st longitude index after leaving east side of exclusion zone around x > 0
! def nlon_p179 - last longitude index before entering west side of exclusion zone around x < 0
! def nlon_m179 - 1st longitude index after leaving east side of exclusion zone around x < 0
! def nlon_m1 - last longitude index before entering west side of exclusion zone around x > 0
! def nlon_m89b - lower bound on nlon_m89
! def nlon_p89b - upper bound on nlon_p89
! def nlon_p91b - lower bound on nlon_p91
! def nlon_p269b - upper bound on nlon_p269
! def nlon_p1b - lower bound on nlon_p1b
! def nlon_p179b - upper bound on nlon_p179
! def nlon_m179b - lower bound on nlon_m179
! def nlon_m1b - upper bound on nlon_m1
!
! def wxp - weights for interpolating from y to longitude in right half plane
! def iwxp - indices for wxp
! def norwxp - interpolation order for wxp
! def wxm - weights for interpolating from y to longitude in left half plane
! def iwxm - indices for wxm
! def norwxm - interpolation order for wxm
! def wyp - weights for interpolating from x to longitude in upper half plane
! def iwyp - indices for wyp
! def norwyp - interpolation order for wyp
! def wym - weights for interpolating from x to longitude in lower half plane
! def iwym - indices for wym
! def norwym - interpolation order for wym
! def wx2lat - weights for interpolating to latitude along longitude lines for y axis excluded.
! def iwx2lat - indices for wx2lat
! def nwx2lat - interpolation order for wx2lat
! def wy2lat - weights for interpolating to latitude along longitude lines for x axis excluded.
! def iwy2lat - indices for wy2lat
! def nwy2lat - interpolation order for wy2lat
! def blendxp - blend weights for right half plane
! def blendxm - blend weights for left half plane
! def blendyp - blend weights for upper half plane
! def blendym - blend weights for lower half plane
! def nlon_p90 - =0 unless nlon is divisible by 4, in which case nlon_p90=nlon/4
! def nlon_p270 - =0 unless nlon is divisible by 4, in which case nlon_p270=3*nlon/4
!
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds
, only: r_kind,i_kind
implicit none
integer(i_kind) norsp
real(r_kind),allocatable,dimension(:,:,:):: xwtxys,ywtxys
integer(i_kind),allocatable,dimension(:,:,:):: ixwtxys,iywtxys
integer(i_kind),allocatable,dimension(:,:):: nxwtxys,nywtxys
contains
subroutine init_smooth_polcas 1
norsp=0
end subroutine init_smooth_polcas
subroutine destroy_smooth_polcas 1
deallocate(xwtxys,ywtxys,ixwtxys,iywtxys,nxwtxys,nywtxys)
end subroutine destroy_smooth_polcas
subroutine setup_smooth_polcas 1,7
!$$$ subprogram documentation block
! . . . .
! subprogram: setup_smooth_polcas create info for smooth_polcas
! prgmmr: parrish org: np23 date: 2005-04-28
!
! abstract: create various information needed for using smooth_polcas.
!
use gridmod, only: nlon,nlat,rlats
use constants, only: zero,half,one,two,pi
use berror, only: nf,nr
integer(i_kind) i,j,jj,ir,iy,nlon4,ix,ilon,iord,k,nin,nxgrid,nygrid
real(r_kind) dlon,df,dr,pi2
real(r_kind) xgrid(-nf:nf),ygrid(-nf:nf),rgrid(-nf:nf),rs(0:nr)
real(r_kind) xlon,ylon
integer(i_kind) nbord
real(r_kind) arg,angle0
integer(i_kind) nor1,iwgt1(0:norsp),nor2,iwgt2(0:norsp)
real(r_kind) wgt1(0:norsp),wgt2(0:norsp)
integer(i_kind) ieven1,ieven2,iodd1,iodd2,nord_evenmax1,nord_oddmax1
integer(i_kind) ii,iord1,iord2,nord_evenmax2,nord_oddmax2
real(r_kind) tin,xin,yin
real(r_kind) slon(nlon),clon(nlon)
! define lat-lon grid lats in polar stereographic units (distance from pole)
do i=1,nr
rs(i)=cos(rlats(nlat-i))/(one+sin(rlats(nlat-i)))
end do
rs(0)=zero
pi2=two*pi
! define polar stereographic grid increment (same in x and y directions)
df=tan(pi2/nlon*half)
dlon=two*pi/nlon
! compute interpolation weights:
! first define xgrid, ygrid, and various other things
nxgrid=2*nf+1
nygrid=2*nf+1
do i=-nf,nf
xgrid(i)=df*i
ygrid(i)=xgrid(i)
end do
allocate(xwtxys(0:norsp,nlon,0:nr),ywtxys(0:norsp,nlon,0:nr))
allocate(ixwtxys(0:norsp,nlon,0:nr),iywtxys(0:norsp,nlon,0:nr))
allocate(nxwtxys(nlon,0:nr),nywtxys(nlon,0:nr))
do j=1,nlon
clon(j)=cos((j-one)*dlon)
slon(j)=sin((j-one)*dlon)
end do
do i=0,nr
do j=1,nlon
!----------------------------x coordinate first
xin=rs(i)*clon(j)
ieven1=nint(xin/df)
iodd1=int(xin/df)
if(xin.lt.zero) iodd1=iodd1-1
nord_oddmax1=min(2*(iodd1+nf)+1,2*(nf-iodd1-1)+1)
nord_evenmax1=2*min(ieven1+nf,nf-ieven1)
iord1=min(norsp,nord_oddmax1,nord_evenmax1)
if(iord1.le.0) then
nor1=0 ; wgt1(0)=one ; iwgt1(0)=min(nf,max(-nf,nint(xin/df)))
end if
if(iord1.gt.0.and.mod(iord1,2).eq.0) then
tin=half+(xin-xgrid(ieven1))/df
call bspline(tin,iord1+1,wgt1)
nor1=iord1
do ii=-iord1/2,iord1/2
iwgt1(ii+iord1/2)=ieven1+ii
end do
end if
if(iord1.gt.0.and.mod(iord1,2).ne.0) then
tin=(xin-xgrid(iodd1))/df
call bspline(tin,iord1+1,wgt1)
nor1=iord1
do ii=-(iord1-1)/2,1+(iord1-1)/2
iwgt1(ii+(iord1-1)/2)=iodd1+ii
end do
end if
!----------------------------y coordinate next
yin=rs(i)*slon(j)
ieven2=nint(yin/df)
iodd2=int(yin/df)
if(yin.lt.zero) iodd2=iodd2-1
nord_oddmax2=min(2*(iodd2+nf)+1,2*(nf-iodd2-1)+1)
nord_evenmax2=2*min(ieven2+nf,nf-ieven2)
iord2=min(norsp,nord_oddmax2,nord_evenmax2)
if(iord2.le.0) then
nor2=0 ; wgt2(0)=one ; iwgt2(0)=min(nf,max(-nf,nint(yin/df)))
end if
if(iord2.gt.0.and.mod(iord2,2).eq.0) then
tin=half+(yin-ygrid(ieven2))/df
call bspline(tin,iord2+1,wgt2)
nor2=iord2
do ii=-iord2/2,iord2/2
iwgt2(ii+iord2/2)=ieven2+ii
end do
end if
if(iord2.gt.0.and.mod(iord2,2).ne.0) then
tin=(yin-ygrid(iodd2))/df
call bspline(tin,iord2+1,wgt2)
nor2=iord2
do ii=-(iord2-1)/2,1+(iord2-1)/2
iwgt2(ii+(iord2-1)/2)=iodd2+ii
end do
end if
!---------------------------now consolidate and get final weights, addresses for this point
nxwtxys(j,i)=nor1
nywtxys(j,i)=nor2
do jj=0,nor2
iywtxys(jj,j,i)=iwgt2(jj)
ywtxys(jj,j,i)=wgt2(jj)
end do
do ii=0,nor1
ixwtxys(ii,j,i)=iwgt1(ii)
xwtxys(ii,j,i)=wgt1(ii)
end do
end do
end do
end subroutine setup_smooth_polcas
subroutine smooth_polcas(fxy,hlatlon) 2,4
!$$$ subprogram documentation block
! . . . .
! subprogram: smooth_polcas 2d smoothing spline interpolation
! prgmmr: parrish org: np22 date: 2005-04-22
!
! abstract: Interpolate polar stereo field fxy(-nf:nf,-nf:nf) to
! corresponding lat-lon field hlatlon(0:nr,nlon), where
! fxy(0,0) is pole point, and hlatlon(0,:) is also nlon copies of
! pole point, using smoothing splines.
!
! program history log:
! 2005-05-14 parrish
!
! input argument list:
! fxy - input data on cartesian grid, dimensions [-nf:nf,-nf:nf].
!
! output argument list:
! hlatlon_out - output data on polar grid, dimensions [0:nlon,0:nr]
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!$$$
use kinds, only: r_kind,i_kind
use constants, only: zero
use berror, only: nf,nr
use gridmod, only: nlon
implicit none
real(r_kind),intent(in),dimension(-nf:nf,-nf:nf):: fxy
real(r_kind),intent(out),dimension(nlon+1,0:nr):: hlatlon
! Declare local arrays variables:
integer(i_kind) i,ii,j,jj,jjj
real(r_kind) sum,ywgt,xywgt
hlatlon=zero
do i=0,nr
do j=1,nlon
sum=zero
do jj=0,nywtxys(j,i)
ywgt=ywtxys(jj,j,i)
jjj=iywtxys(jj,j,i)
do ii=0,nxwtxys(j,i)
xywgt=ywgt*xwtxys(ii,j,i)
sum=sum+xywgt*fxy(ixwtxys(ii,j,i),jjj)
end do
end do
hlatlon(j,i)=sum
end do
end do
end subroutine smooth_polcas
subroutine smooth_polcasa(fxy,hlatlon) 2,4
!$$$ subprogram documentation block
! . . . .
! subprogram: smooth_polcasa adjoint of smooth_polcas
! prgmmr: parrish org: np22 date: 2005-04-22
!
! abstract: adjoint of smooth_polcas
!
! program history log:
! 2005-05-14 parrish
!
! input argument list:
! fxy - input data on cartesian grid, dimensions [-nf:nf,-nf:nf].
!
! output argument list:
! hlatlon_out - output data on polar grid, dimensions [0:nlon,0:nr]
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!$$$
use kinds, only: r_kind,i_kind
use constants, only: zero
use berror, only: nf,nr
use gridmod, only: nlon
implicit none
real(r_kind),intent(out),dimension(-nf:nf,-nf:nf):: fxy
real(r_kind),intent(in),dimension(nlon+1,0:nr):: hlatlon
! Declare local arrays variables:
integer(i_kind) i,ii,j,jj,jjj
real(r_kind) sum,ywgt,xywgt
fxy=zero
do i=0,nr
do j=1,nlon
sum=hlatlon(j,i)
do jj=0,nywtxys(j,i)
ywgt=ywtxys(jj,j,i)
jjj=iywtxys(jj,j,i)
do ii=0,nxwtxys(j,i)
xywgt=ywgt*xwtxys(ii,j,i)
fxy(ixwtxys(ii,j,i),jjj)=fxy(ixwtxys(ii,j,i),jjj)+sum*xywgt
end do
end do
end do
end do
end subroutine smooth_polcasa
subroutine bspline(tin,k,wout) 4,1
! compute weights for a bspline of order k (degree k-1, continuous k-2)
! k >= 1
! for unit spaced grid.
! 0 <= tin <= 1
use constants, only: zero,one
integer(i_kind) k
real(r_kind) tin,wout(0:k-1)
integer(i_kind) i,m
real(r_kind) t,w(0:k),rmi(k)
do m=1,k
rmi(m)=one/m
end do
t=tin+k-1
w=zero
w(k-1)=one
if(k.gt.1) then
do m=2,k
do i=0,k-1
w(i)=w(i)*(t-i)*rmi(m-1) + w(i+1)*(i+m-t)*rmi(m-1)
end do
end do
end if
do i=0,k-1
wout(i)=w(i)
end do
end subroutine bspline
end module smooth_polcarf