hinterp.F90 Source File


Source Code

      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