subroutine hinterp ( qin,iin,jin,qout,iout,jout,mlev,undef ) 5,2
implicit none
integer iin,jin, iout,jout, mlev
real qin(iin,jin,mlev), qout(iout,jout,mlev)
real undef,pi,dlin,dpin,dlout,dpout
real dlam(iin), lons(iout*jout), lon
real dphi(jin), lats(iout*jout), lat
integer i,j,loc
pi = 4.0*atan(1.0)
dlin = 2*pi/iin
dpin = pi/(jin-1)
dlam(:) = dlin
dphi(:) = dpin
dlout = 2*pi/iout
dpout = pi/(jout-1)
loc = 0
do j=1,jout
do i=1,iout
loc = loc + 1
lon = -pi + (i-1)*dlout
lons(loc) = lon
enddo
enddo
loc = 0
do j=1,jout
lat = -pi/2.0 + (j-1)*dpout
do i=1,iout
loc = loc + 1
lats(loc) = lat
enddo
enddo
call interp_hh
( qin,iin,jin,mlev,dlam,dphi,
. qout,iout*jout,lons,lats,undef, -pi )
return
end
!.......................................................................................................
subroutine hhinterp ( qin,iin,jin,qout,iout,jout,mlev,undef, 1,1
. lons_in, lats_in )
implicit none
integer iin,jin, iout,jout, mlev
real qin(iin,jin,mlev), qout(iout,jout,mlev)
real undef,pi,dlin,dpin,dlout,dpout
real :: lons_in(iin), lats_in(jin)
real dlam(iin), lons(iout*jout), lon
real dphi(jin), lats(iout*jout), lat
integer i,j,loc,rc
real lon_min
character(len=*), parameter :: Iam='hinterp_'
integer :: i1, i2, i3, i4, n1, n2, n3, n4, imh, k
real lons_save(iin)
pi = 4.0*atan(1.0)
! Larry's code expects longitudes in the range [-180,180],
! therefore redefine longitudes and flip fields if input
! has longitude in the range [0,360].
! -------------------------------------------------------
lons_save = lons_in
imh = -1
if ( lons_in(1) >= 0 ) then
do i = 1, iin
if ( lons_in(i) >= pi ) then
lons_in(i) = lons_in(i) - 2*pi
if ( imh < 0 ) imh = i
end if
end do
i1 = 1; i2 = imh-1; i3=imh; i4=iin
n1 = 1; n2 = i4-i3+1; n3=n2+1; n4=iin
dlam = lons_in
lons_in(n1:n2) = dlam(i3:i4)
lons_in(n3:n4) = dlam(i1:i2)
do k = 1, mlev
do j = 1, jin
dlam = qin(:,j,k)
qin(n1:n2,j,k) = dlam(i3:i4)
qin(n3:n4,j,k) = dlam(i1:i2)
end do
end do
end if
do i = 1, iin-1
dlam(i) = lons_in(i+1) - lons_in(i)
enddo
dlam(iin) = lons_in(iin) - lons_in(iin-1) ! just in case
lon_min = lons_in(1)
if ( abs(lats_in(1)+pi/2) >= 0.0001 ) then
qout = undef ! requires lat to start at South Pole
return
end if
do j = 1, jin-1
dphi(j) = lats_in(j+1) - lats_in(j)
enddo
dphi(jin) = lats_in(jin) - lats_in(jin-1) ! just in case
dlout = 2*pi/iout
dpout = pi/(jout-1)
loc = 0
do j=1,jout
do i=1,iout
loc = loc + 1
lon = -pi + (i-1)*dlout ! in [-180,180)
lons(loc) = lon
enddo
enddo
loc = 0
do j=1,jout
lat = -pi/2.0 + (j-1)*dpout
do i=1,iout
loc = loc + 1
lats(loc) = lat
enddo
enddo
call interp_hh ( qin,iin,jin,mlev,dlam,dphi,
. qout,iout*jout,lons,lats,undef, lon_min )
! Return input to original form
! -----------------------------
if ( imh > 0 ) then
lons_in = lons_save
do k = 1, mlev
do j = 1, jin
dlam = qin(:,j,k)
qin(i1:i2,j,k) = dlam(n3:n4)
qin(i3:i4,j,k) = dlam(n1:n2)
end do
end do
end if
return
end
subroutine interp_hh ( q_cmp,im,jm,lm,dlam,dphi, 2
. q_geo,irun,lon_geo,lat_geo, undef, lon_min )
C***********************************************************************
C
C PURPOSE:
C ========
C Performs a horizontal interpolation from a field on a computational grid
C to arbitrary locations.
C
C INPUT:
C ======
C q_cmp ...... Field q_cmp(im,jm,lm) on the computational grid
C im ......... Longitudinal dimension of q_cmp
C jm ......... Latitudinal dimension of q_cmp
C lm ......... Vertical dimension of q_cmp
C dlam ....... Computational Grid Delta Lambda
C dphi ....... Computational Grid Delta Phi
C irun ....... Number of Output Locations
C lon_geo .... Longitude Location of Output
C lat_geo .... Latitude Location of Output
C
C OUTPUT:
C =======
C q_geo ...... Field q_geo(irun,lm) at arbitrary locations
C
C
C***********************************************************************
C* GODDARD LABORATORY FOR ATMOSPHERES *
C***********************************************************************
implicit none
c Input Variables
c ---------------
integer im,jm,lm,irun
real q_geo(irun,lm)
real lon_geo(irun)
real lat_geo(irun)
real q_cmp(im,jm,lm)
real dlam(im)
real dphi(jm)
real :: lon_min
c Local Variables
c ---------------
integer i,j,l,m,n
integer, allocatable :: ip1(:), ip0(:), im1(:), im2(:)
integer, allocatable :: jp1(:), jp0(:), jm1(:), jm2(:)
integer ip1_for_jp1, ip0_for_jp1, im1_for_jp1, im2_for_jp1
integer ip1_for_jm2, ip0_for_jm2, im1_for_jm2, im2_for_jm2
integer jm2_for_jm2, jp1_for_jp1
c Bi-Linear Weights
c -----------------
real, allocatable :: wl_ip0jp0 (:)
real, allocatable :: wl_im1jp0 (:)
real, allocatable :: wl_ip0jm1 (:)
real, allocatable :: wl_im1jm1 (:)
c Bi-Cubic Weights
c ----------------
real, allocatable :: wc_ip1jp1 (:)
real, allocatable :: wc_ip0jp1 (:)
real, allocatable :: wc_im1jp1 (:)
real, allocatable :: wc_im2jp1 (:)
real, allocatable :: wc_ip1jp0 (:)
real, allocatable :: wc_ip0jp0 (:)
real, allocatable :: wc_im1jp0 (:)
real, allocatable :: wc_im2jp0 (:)
real, allocatable :: wc_ip1jm1 (:)
real, allocatable :: wc_ip0jm1 (:)
real, allocatable :: wc_im1jm1 (:)
real, allocatable :: wc_im2jm1 (:)
real, allocatable :: wc_ip1jm2 (:)
real, allocatable :: wc_ip0jm2 (:)
real, allocatable :: wc_im1jm2 (:)
real, allocatable :: wc_im2jm2 (:)
real ux, ap1, ap0, am1, am2
real uy, bp1, bp0, bm1, bm2
real lon_cmp(im)
real lat_cmp(jm)
real q_tmp(irun)
real pi,d
real lam,lam_ip1,lam_ip0,lam_im1,lam_im2
real phi,phi_jp1,phi_jp0,phi_jm1,phi_jm2
real dl,dp,phi_np,lam_0
real lam_geo, lam_cmp
real phi_geo, phi_cmp
real undef
integer im1_cmp,icmp
integer jm1_cmp,jcmp
c Initialization
c --------------
pi = 4.*atan(1.)
dl = 2*pi/ im ! Uniform Grid Delta Lambda
dp = pi/(jm-1) ! Uniform Grid Delta Phi
c Allocate Memory for Weights and Index Locations
c -----------------------------------------------
allocate ( wl_ip0jp0(irun) , wl_im1jp0(irun) )
allocate ( wl_ip0jm1(irun) , wl_im1jm1(irun) )
allocate ( wc_ip1jp1(irun) , wc_ip0jp1(irun) , wc_im1jp1(irun) , wc_im2jp1(irun) )
allocate ( wc_ip1jp0(irun) , wc_ip0jp0(irun) , wc_im1jp0(irun) , wc_im2jp0(irun) )
allocate ( wc_ip1jm1(irun) , wc_ip0jm1(irun) , wc_im1jm1(irun) , wc_im2jm1(irun) )
allocate ( wc_ip1jm2(irun) , wc_ip0jm2(irun) , wc_im1jm2(irun) , wc_im2jm2(irun) )
allocate ( ip1(irun) , ip0(irun) , im1(irun) , im2(irun) )
allocate ( jp1(irun) , jp0(irun) , jm1(irun) , jm2(irun) )
c Compute Input Computational-Grid Latitude and Longitude Locations
c -----------------------------------------------------------------
lon_cmp(1) = lon_min ! user supplied orign
do i=2,im
lon_cmp(i) = lon_cmp(i-1) + dlam(i-1)
enddo
lat_cmp(1) = -pi*0.5
do j=2,jm-1
lat_cmp(j) = lat_cmp(j-1) + dphi(j-1)
enddo
lat_cmp(jm) = pi*0.5
c Compute Weights for Computational to Geophysical Grid Interpolation
c -------------------------------------------------------------------
do i=1,irun
lam_cmp = lon_geo(i)
phi_cmp = lat_geo(i)
c Determine Indexing Based on Computational Grid
c ----------------------------------------------
im1_cmp = 1
do icmp = 2,im
if( lon_cmp(icmp).lt.lam_cmp ) im1_cmp = icmp
enddo
jm1_cmp = 1
do jcmp = 2,jm
if( lat_cmp(jcmp).lt.phi_cmp ) jm1_cmp = jcmp
enddo
im1(i) = im1_cmp
ip0(i) = im1(i) + 1
ip1(i) = ip0(i) + 1
im2(i) = im1(i) - 1
jm1(i) = jm1_cmp
jp0(i) = jm1(i) + 1
jp1(i) = jp0(i) + 1
jm2(i) = jm1(i) - 1
c Fix Longitude Index Boundaries
c ------------------------------
if(im1(i).eq.im) then
ip0(i) = 1
ip1(i) = 2
endif
if(im1(i).eq.1) then
im2(i) = im
endif
if(ip0(i).eq.im) then
ip1(i) = 1
endif
c Compute Immediate Surrounding Coordinates
c -----------------------------------------
lam = lam_cmp
phi = phi_cmp
c Compute and Adjust Longitude Weights
c ------------------------------------
lam_im2 = lon_cmp(im2(i))
lam_im1 = lon_cmp(im1(i))
lam_ip0 = lon_cmp(ip0(i))
lam_ip1 = lon_cmp(ip1(i))
if( lam_im2.gt.lam_im1 ) lam_im2 = lam_im2 - 2*pi
if( lam_im1.gt.lam_ip0 ) lam_ip0 = lam_ip0 + 2*pi
if( lam_im1.gt.lam_ip1 ) lam_ip1 = lam_ip1 + 2*pi
if( lam_ip0.gt.lam_ip1 ) lam_ip1 = lam_ip1 + 2*pi
c Compute and Adjust Latitude Weights
c Note: Latitude Index Boundaries are Adjusted during Interpolation
c ------------------------------------------------------------------
phi_jm1 = lat_cmp(jm1(i))
if( jm2(i).eq.0 ) then
phi_jm2 = phi_jm1 - dphi(1)
else
phi_jm2 = lat_cmp(jm2(i))
endif
if( jm1(i).eq.jm ) then
phi_jp0 = phi_jm1 + dphi(jm-1)
phi_jp1 = phi_jp0 + dphi(jm-2)
else
phi_jp0 = lat_cmp(jp0(i))
if( jp1(i).eq.jm+1 ) then
phi_jp1 = phi_jp0 + dphi(jm-1)
else
phi_jp1 = lat_cmp(jp1(i))
endif
endif
c Bi-Linear Weights
c -----------------
d = (lam_ip0-lam_im1)*(phi_jp0-phi_jm1)
wl_im1jm1(i) = (lam_ip0-lam )*(phi_jp0-phi )/d
wl_ip0jm1(i) = (lam -lam_im1)*(phi_jp0-phi )/d
wl_im1jp0(i) = (lam_ip0-lam )*(phi -phi_jm1)/d
wl_ip0jp0(i) = (lam -lam_im1)*(phi -phi_jm1)/d
c Bi-Cubic Weights
c ----------------
ap1 = ( (lam -lam_ip0)*(lam -lam_im1)*(lam -lam_im2) )
. / ( (lam_ip1-lam_ip0)*(lam_ip1-lam_im1)*(lam_ip1-lam_im2) )
ap0 = ( (lam_ip1-lam )*(lam -lam_im1)*(lam -lam_im2) )
. / ( (lam_ip1-lam_ip0)*(lam_ip0-lam_im1)*(lam_ip0-lam_im2) )
am1 = ( (lam_ip1-lam )*(lam_ip0-lam )*(lam -lam_im2) )
. / ( (lam_ip1-lam_im1)*(lam_ip0-lam_im1)*(lam_im1-lam_im2) )
am2 = ( (lam_ip1-lam )*(lam_ip0-lam )*(lam_im1-lam ) )
. / ( (lam_ip1-lam_im2)*(lam_ip0-lam_im2)*(lam_im1-lam_im2) )
bp1 = ( (phi -phi_jp0)*(phi -phi_jm1)*(phi -phi_jm2) )
. / ( (phi_jp1-phi_jp0)*(phi_jp1-phi_jm1)*(phi_jp1-phi_jm2) )
bp0 = ( (phi_jp1-phi )*(phi -phi_jm1)*(phi -phi_jm2) )
. / ( (phi_jp1-phi_jp0)*(phi_jp0-phi_jm1)*(phi_jp0-phi_jm2) )
bm1 = ( (phi_jp1-phi )*(phi_jp0-phi )*(phi -phi_jm2) )
. / ( (phi_jp1-phi_jm1)*(phi_jp0-phi_jm1)*(phi_jm1-phi_jm2) )
bm2 = ( (phi_jp1-phi )*(phi_jp0-phi )*(phi_jm1-phi ) )
. / ( (phi_jp1-phi_jm2)*(phi_jp0-phi_jm2)*(phi_jm1-phi_jm2) )
wc_ip1jp1(i) = bp1*ap1
wc_ip0jp1(i) = bp1*ap0
wc_im1jp1(i) = bp1*am1
wc_im2jp1(i) = bp1*am2
wc_ip1jp0(i) = bp0*ap1
wc_ip0jp0(i) = bp0*ap0
wc_im1jp0(i) = bp0*am1
wc_im2jp0(i) = bp0*am2
wc_ip1jm1(i) = bm1*ap1
wc_ip0jm1(i) = bm1*ap0
wc_im1jm1(i) = bm1*am1
wc_im2jm1(i) = bm1*am2
wc_ip1jm2(i) = bm2*ap1
wc_ip0jm2(i) = bm2*ap0
wc_im1jm2(i) = bm2*am1
wc_im2jm2(i) = bm2*am2
enddo
c Interpolate Computational-Grid Quantities to Geophysical Grid
c -------------------------------------------------------------
do L=1,lm
do i=1,irun
if( lat_geo(i).le.lat_cmp(2) .or.
. lat_geo(i).ge.lat_cmp(jm-1) ) then
c 1st Order Interpolation at Poles
c --------------------------------
if( q_cmp( im1(i),jm1(i),L ).ne.undef .and.
. q_cmp( ip0(i),jm1(i),L ).ne.undef .and.
. q_cmp( im1(i),jp0(i),L ).ne.undef .and.
. q_cmp( ip0(i),jp0(i),L ).ne.undef ) then
q_tmp(i) = wl_im1jm1(i) * q_cmp( im1(i),jm1(i),L )
. + wl_ip0jm1(i) * q_cmp( ip0(i),jm1(i),L )
. + wl_im1jp0(i) * q_cmp( im1(i),jp0(i),L )
. + wl_ip0jp0(i) * q_cmp( ip0(i),jp0(i),L )
else
q_tmp(i) = undef
endif
else
c Cubic Interpolation away from Poles
c -----------------------------------
if( q_cmp( ip1(i),jp0(i),L ).ne.undef .and.
. q_cmp( ip0(i),jp0(i),L ).ne.undef .and.
. q_cmp( im1(i),jp0(i),L ).ne.undef .and.
. q_cmp( im2(i),jp0(i),L ).ne.undef .and.
. q_cmp( ip1(i),jm1(i),L ).ne.undef .and.
. q_cmp( ip0(i),jm1(i),L ).ne.undef .and.
. q_cmp( im1(i),jm1(i),L ).ne.undef .and.
. q_cmp( im2(i),jm1(i),L ).ne.undef .and.
. q_cmp( ip1(i),jp1(i),L ).ne.undef .and.
. q_cmp( ip0(i),jp1(i),L ).ne.undef .and.
. q_cmp( im1(i),jp1(i),L ).ne.undef .and.
. q_cmp( im2(i),jp1(i),L ).ne.undef .and.
. q_cmp( ip1(i),jm2(i),L ).ne.undef .and.
. q_cmp( ip0(i),jm2(i),L ).ne.undef .and.
. q_cmp( im1(i),jm2(i),L ).ne.undef .and.
. q_cmp( im2(i),jm2(i),L ).ne.undef ) then
q_tmp(i) = wc_ip1jp1(i) * q_cmp( ip1(i),jp1(i),L )
. + wc_ip0jp1(i) * q_cmp( ip0(i),jp1(i),L )
. + wc_im1jp1(i) * q_cmp( im1(i),jp1(i),L )
. + wc_im2jp1(i) * q_cmp( im2(i),jp1(i),L )
. + wc_ip1jp0(i) * q_cmp( ip1(i),jp0(i),L )
. + wc_ip0jp0(i) * q_cmp( ip0(i),jp0(i),L )
. + wc_im1jp0(i) * q_cmp( im1(i),jp0(i),L )
. + wc_im2jp0(i) * q_cmp( im2(i),jp0(i),L )
. + wc_ip1jm1(i) * q_cmp( ip1(i),jm1(i),L )
. + wc_ip0jm1(i) * q_cmp( ip0(i),jm1(i),L )
. + wc_im1jm1(i) * q_cmp( im1(i),jm1(i),L )
. + wc_im2jm1(i) * q_cmp( im2(i),jm1(i),L )
. + wc_ip1jm2(i) * q_cmp( ip1(i),jm2(i),L )
. + wc_ip0jm2(i) * q_cmp( ip0(i),jm2(i),L )
. + wc_im1jm2(i) * q_cmp( im1(i),jm2(i),L )
. + wc_im2jm2(i) * q_cmp( im2(i),jm2(i),L )
elseif( q_cmp( im1(i),jm1(i),L ).ne.undef .and.
. q_cmp( ip0(i),jm1(i),L ).ne.undef .and.
. q_cmp( im1(i),jp0(i),L ).ne.undef .and.
. q_cmp( ip0(i),jp0(i),L ).ne.undef ) then
q_tmp(i) = wl_im1jm1(i) * q_cmp( im1(i),jm1(i),L )
. + wl_ip0jm1(i) * q_cmp( ip0(i),jm1(i),L )
. + wl_im1jp0(i) * q_cmp( im1(i),jp0(i),L )
. + wl_ip0jp0(i) * q_cmp( ip0(i),jp0(i),L )
else
q_tmp(i) = undef
endif
endif
enddo
c Load Temp array into Output array
c ---------------------------------
do i=1,irun
q_geo(i,L) = q_tmp(i)
enddo
enddo
deallocate ( wl_ip0jp0 , wl_im1jp0 )
deallocate ( wl_ip0jm1 , wl_im1jm1 )
deallocate ( wc_ip1jp1 , wc_ip0jp1 , wc_im1jp1 , wc_im2jp1 )
deallocate ( wc_ip1jp0 , wc_ip0jp0 , wc_im1jp0 , wc_im2jp0 )
deallocate ( wc_ip1jm1 , wc_ip0jm1 , wc_im1jm1 , wc_im2jm1 )
deallocate ( wc_ip1jm2 , wc_ip0jm2 , wc_im1jm2 , wc_im2jm2 )
deallocate ( ip1 , ip0 , im1 , im2 )
deallocate ( jp1 , jp0 , jm1 , jm2 )
return
end