subroutine hinterp2 ( qin,iin,jin,qout,iout,jout,mlev,undef ) 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_in (iin) real dphi_in (jin) real lons(iout), lons_out(iout*jout) real lats(jout), lats_out(iout*jout) integer i,j,loc pi = 4.0*atan(1.0) dlin = 2*pi/ iin dpin = pi/(jin-1) dlout = 2*pi/ iout dpout = pi/(jout-1) ! Compute Input DLAM & DPHI ! ------------------------- dlam_in(:) = dlin dphi_in(:) = dpin ! Compute Output Lons & Lats ! -------------------------- lons(1) = -pi do i=2,iout lons(i) = lons(i-1) + dlout enddo lats(1) = -pi*0.5 do j=2,jout-1 lats(j) = lats(j-1) + dpout enddo lats(jout) = pi*0.5 loc = 0 do j=1,jout do i=1,iout loc = loc + 1 lons_out(loc) = lons(i) enddo enddo loc = 0 do j=1,jout do i=1,iout loc = loc + 1 lats_out(loc) = lats(j) enddo enddo call interp_hh ( qin,iin,jin,mlev,dlam_in,dphi_in, & qout,iout*jout,lons_out,lats_out,undef, -pi ) return end !....................................................................................................... subroutine hhinterp ( qin,iin,jin,qout,iout,jout,mlev,undef, & 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 lons_out(iout), lats_out(jout) real dlam(iin), lons(iout*jout) real dphi(jin), lats(iout*jout) integer i,j,loc real lon_min character(len=*), parameter :: Iam='hinterp_' integer :: imh, k real lons_save(iin) pi = 4.0*atan(1.0) dlin = 2*pi/ iin dpin = pi/(jin-1) dlout = 2*pi/ iout dpout = pi/(jout-1) ! 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 imh=1 do k = 1, mlev call myhflip2_(qin(1,1,k),iin,jin) enddo lons_in(1) = -pi do i=2,iin lons_in(i) = lons_in(i-1) + dlin enddo endif ! Set INPUT DLAM and DPHI to Uniform Grid Constants ! ------------------------------------------------- dlam(:) = dlin dphi(:) = dpin 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 ! Compute Output Lons & Lats consistent with DLAM and DPHI ! -------------------------------------------------------- lons_out(1) = -pi do i=2,iout lons_out(i) = lons_out(i-1) + dlout enddo lats_out(1) = -pi*0.5 do j=2,jout-1 lats_out(j) = lats_out(j-1) + dpout enddo lats_out(jout) = pi*0.5 loc = 0 do j=1,jout do i=1,iout loc = loc + 1 lons(loc) = lons_out(i) enddo enddo loc = 0 do j=1,jout do i=1,iout loc = loc + 1 lats(loc) = lats_out(j) 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 call myhflip2_(qin(1,1,k),iin,jin) enddo endif return contains subroutine myhflip2_ ( q,im,jm ) implicit none integer im,jm,i,j real, intent(inout) :: q(im,jm) real, allocatable :: dum(:) allocate ( dum(im) ) do j=1,jm do i=1,im/2 dum(i) = q(i+im/2,j) dum(i+im/2) = q(i,j) enddo q(:,j) = dum(:) enddo deallocate ( dum ) end subroutine myhflip2_ end subroutine interp_hh ( q_cmp,im,jm,lm,dlam,dphi, & q_geo,irun,lon_geo,lat_geo, undef, lon_min) !*********************************************************************** ! ! PURPOSE: ! ======== ! Performs a horizontal interpolation from a field on a computational grid ! to arbitrary locations. ! ! INPUT: ! ====== ! q_cmp ...... Field q_cmp(im,jm,lm) on the computational grid ! im ......... Longitudinal dimension of q_cmp ! jm ......... Latitudinal dimension of q_cmp ! lm ......... Vertical dimension of q_cmp ! dlam ....... Computational Grid Delta Lambda ! dphi ....... Computational Grid Delta Phi ! irun ....... Number of Output Locations ! lon_geo .... Longitude Location of Output ! lat_geo .... Latitude Location of Output ! ! OUTPUT: ! ======= ! q_geo ...... Field q_geo(irun,lm) at arbitrary locations ! ! !*********************************************************************** !* GODDARD LABORATORY FOR ATMOSPHERES * !*********************************************************************** implicit none ! Input Variables ! --------------- 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 ! Local Variables ! --------------- integer i,j,l integer, allocatable :: ip1(:), ip0(:), im1(:), im2(:) integer, allocatable :: jp1(:), jp0(:), jm1(:), jm2(:) ! Bi-Linear Weights ! ----------------- real, allocatable :: wl_ip0jp0 (:) real, allocatable :: wl_im1jp0 (:) real, allocatable :: wl_ip0jm1 (:) real, allocatable :: wl_im1jm1 (:) ! Bi-Cubic Weights ! ---------------- 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 ap1, ap0, am1, am2 real 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 lam_cmp real phi_cmp real undef integer im1_cmp,icmp integer jm1_cmp,jcmp ! Initialization ! -------------- pi = 4.*atan(1.) ! Allocate Memory for Weights and Index Locations ! ----------------------------------------------- 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) ) ! Compute Input Computational-Grid Latitude and Longitude Locations ! ----------------------------------------------------------------- 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 ! Compute Weights for Computational to Geophysical Grid Interpolation ! ------------------------------------------------------------------- do i=1,irun lam_cmp = lon_geo(i) phi_cmp = lat_geo(i) ! Determine Indexing Based on Computational Grid ! ---------------------------------------------- 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 ! Fix Longitude Index Boundaries ! ------------------------------ 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 ! Compute Immediate Surrounding Coordinates ! ----------------------------------------- lam = lam_cmp phi = phi_cmp ! Compute and Adjust Longitude Weights ! ------------------------------------ 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 ! Compute and Adjust Latitude Weights ! Note: Latitude Index Boundaries are Adjusted during Interpolation ! ------------------------------------------------------------------ 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 ! Bi-Linear Weights ! ----------------- 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 ! Bi-Cubic Weights ! ---------------- 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 ! Interpolate Computational-Grid Quantities to Geophysical Grid ! ------------------------------------------------------------- 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 ! 1st Order Interpolation at Poles ! -------------------------------- 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 ! Cubic Interpolation away from Poles ! ----------------------------------- 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 ! Load Temp array into Output array ! --------------------------------- 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