#include <misc.h>#include <params.h>subroutine sltint(kdim    ,ikcnst  ,jcen    ,fb      ,lam     , &                  rdphi   ,rdz     ,lbasdy  ,wdz     ,xl      , &                  xr      ,wgt1x   ,wgt2x   ,wgt3x   ,wgt4x   , &                  hl      ,hr      ,dhl     ,dhr     ,ys      , &                  yn      ,wgt1y   ,wgt2y   ,wgt3y   ,wgt4y   , &                  hs      ,hn      ,dhs     ,dhn     ,wgt1z   , &                  wgt2z   ,wgt3z   ,wgt4z   ,hb      ,ht      , &                  dhb     ,dht     ,idp     ,jdp     ,kdp     , &                  kkdp    ,lhrzint ,lvrtint ,limdrh  ,limdrv  , &                  fdp     ,nlon    )!!-----------------------------------------------------------------------!! Purpose:! Interpolate field to departure points using Hermite or Lagrange! Cubic interpolation!! Author:  J. Olson!!-----------------------------------------------------------------------!! $Id: sltint.F90,v 1.3 2000/06/14 19:32:12 erik Exp $! $Author: erik $!!-----------------------------------------------------------------------  use precision  use pmgrid#if (!defined CRAY)  use srchutil, only: wheneq#endif  implicit none!----------------------------- Parameters ------------------------------!  integer, parameter :: ppdy = 4 ! length of interpolation grid stencil  integer, parameter :: ppdz = 4 ! length of interp. grid stencil in z!!------------------------------Arguments--------------------------------!  integer , intent(in)   :: kdim                ! vertical dimension  integer , intent(in)   :: ikcnst              ! constituent index  integer , intent(in)   :: jcen                ! index of latitude in extended grid  real(r8), intent(in)   :: fb(plond,kdim*ikcnst,beglatex:endlatex) ! input field  real(r8), intent(in)   :: lam   (plond,platd) ! longitude coordinates of model grid  real(r8), intent(in)   :: rdphi (plon,plev)   ! reciprocal of y-interval  real(r8), intent(in)   :: rdz   (plon,plev)   ! reciprocal of z-interval  real(r8), intent(in)   :: lbasdy(4,2,platd)   ! basis functions for lat deriv est.  real(r8), intent(in)   :: wdz   (4,2,kdim)    ! basis functions for vert deriv est.  real(r8), intent(in)   :: xl    (plon,plev,4) ! weight for x-interpolants (left)  real(r8), intent(in)   :: xr    (plon,plev,4) ! weight for x-interpolants (right)  real(r8), intent(in)   :: wgt1x (plon,plev,4) ! weight for x-interpolants (Lag Cubic)  real(r8), intent(in)   :: wgt2x (plon,plev,4) ! weight for x-interpolants (Lag Cubic)  real(r8), intent(in)   :: wgt3x (plon,plev,4) ! weight for x-interpolants (Lag Cubic)  real(r8), intent(in)   :: wgt4x (plon,plev,4) ! weight for x-interpolants (Lag Cubic)  real(r8), intent(in)   :: hl    (plon,plev,4) ! weight for x-interpolants (Hermite)  real(r8), intent(in)   :: hr    (plon,plev,4) ! weight for x-interpolants (Hermite)  real(r8), intent(in)   :: dhl   (plon,plev,4) ! weight for x-interpolants (Hermite)  real(r8), intent(in)   :: dhr   (plon,plev,4) ! weight for x-interpolants (Hermite)  real(r8), intent(in)   :: ys    (plon,plev)   ! weight for y-interpolants (south)  real(r8), intent(in)   :: yn    (plon,plev)   ! weight for y-interpolants (north)  real(r8), intent(in)   :: wgt1y (plon,plev)   ! weight for y-interpolants (Lag Cubic)  real(r8), intent(in)   :: wgt2y (plon,plev)   ! weight for y-interpolants (Lag Cubic)  real(r8), intent(in)   :: wgt3y (plon,plev)   ! weight for y-interpolants (Lag Cubic)  real(r8), intent(in)   :: wgt4y (plon,plev)   ! weight for y-interpolants (Lag Cubic)  real(r8), intent(in)   :: hs    (plon,plev)   ! weight for y-interpolants (Hermite)  real(r8), intent(in)   :: hn    (plon,plev)   ! weight for y-interpolants (Hermite)  real(r8), intent(in)   :: dhs   (plon,plev)   ! weight for y-interpolants (Hermite)  real(r8), intent(in)   :: dhn   (plon,plev)   ! weight for y-interpolants (Hermite)  real(r8), intent(in)   :: wgt1z (plon,plev)   ! weight for z-interpolants (Lag Cubic)  real(r8), intent(in)   :: wgt2z (plon,plev)   ! weight for z-interpolants (Lag Cubic)  real(r8), intent(in)   :: wgt3z (plon,plev)   ! weight for z-interpolants (Lag Cubic)  real(r8), intent(in)   :: wgt4z (plon,plev)   ! weight for z-interpolants (Lag Cubic)  real(r8), intent(in)   :: hb    (plon,plev)   ! weight for z-interpolants (Hermite)  real(r8), intent(in)   :: ht    (plon,plev)   ! weight for z-interpolants (Hermite)  real(r8), intent(in)   :: dhb   (plon,plev)   ! weight for z-interpolants (Hermite)  real(r8), intent(in)   :: dht   (plon,plev)   ! weight for z-interpolants (Hermite)  integer , intent(in)   :: idp   (plon,plev,4) ! index of x-coordinate of dep pt  integer , intent(in)   :: jdp   (plon,plev)   ! index of y-coordinate of dep pt  integer , intent(in)   :: kdp   (plon,plev)   ! index of z-coordinate of dep pt  integer , intent(in)   :: kkdp  (plon,plev)   ! index of z-coordinate of dep pt (alt)  logical , intent(in)   :: lhrzint             ! flag to do horizontal interpolation  logical , intent(in)   :: lvrtint             ! flag to do vertical   interpolation  logical , intent(in)   :: limdrh              ! horizontal derivative limiter flag  logical , intent(in)   :: limdrv              ! vertical   derivative limiter flag  real(r8), intent(out)  :: fdp(plon,plev)      ! interpolant  integer , intent(in)   :: nlon                ! number of longitudes for this latitude!!---------------------------Local workspace-----------------------------!  integer i                           ! |  integer ii1,ii2,ii3,ii4             ! |  integer ii,jj,j                     ! |  integer k                           ! |  integer kk                          ! |  integer jmin                        ! |  integer jmax                        ! | -- indices  integer jdpval                      ! |  integer kdpval                      ! |  integer icount                      ! |  integer indx(plon)                  ! |  integer nval                        ! |  integer kdimm1                      ! |  integer kdimm2                      ! |  integer kdimm3                      ! |!  real(r8) fac                        ! factor applied in limiter  real(r8) tmp1                       ! derivative factor  real(r8) tmp2                       ! abs(tmp1)  real(r8) deli                       ! linear derivative  real(r8) dx                         ! delta x  real(r8) rdx (platd)                ! 1./dx  real(r8) rdx6(platd)                ! 1./(6*dx)  real(r8) fxl                        ! left  derivative estimate  real(r8) fxr                        ! right derivative estimate!  real(r8) f1                         ! |  real(r8) f2                         ! |  real(r8) f3                         ! |  real(r8) f4                         ! |  real(r8) tmptop                     ! | -- work arrays  real(r8) tmpbot                     ! |  real(r8) fintx(plon,plev,ppdy,ppdz) ! |  real(r8) finty(plon,plev     ,ppdz) ! |  real(r8) fbot (plon,plev,ppdz)      ! |  real(r8) ftop (plon,plev,ppdz)      ! |!!-----------------------------------------------------------------------!  fac  = 3.*(1. - 10.*epsilon(fac))  kdimm1 = kdim - 1  kdimm2 = kdim - 2  kdimm3 = kdim - 3!  do j = 1,platd     dx      = lam(nxpt+2,j) - lam(nxpt+1,j)     rdx (j) = 1./dx     rdx6(j) = 1./(6.*dx)  end do!!-----------------------------------------------------------------------!------------------------- Code Description ----------------------------!! Each block of code performs a specific interpolation as follows:!!  For 3-D (horizontal AND vertical) interpolation:!!    10XX loops: Hermite  cubic/linear interpolation in the horizontal!    20XX loops: Lagrange cubic/linear interpolation in the horizontal!    30XX loops: Hermite  cubic/linear interpolation in the vertical!    40XX loops: Lagrange cubic/linear interpolation in the vertical!!  For horizontal interpolation only:!!    50XX loops: an optimized Lagrange cubic/linear algorithm!    60XX loops: Hermite cubic/linear interpolation in the horizontal!!  For vertical interpolation only:!!    70XX loops: an optimized Lagrange cubic/linear algorithm (no!                Hermite interpolator available)!!-----------------------------------------------------------------------!-----------------------------------------------------------------------!  if( lhrzint .and. lvrtint ) then!!-----------------------------------------------------------------------!-----------------------------------------------------------------------!!    10XX loops: Hermite  cubic/linear interpolation in the horizontal!!-----------------------------------------------------------------------!-----------------------------------------------------------------------!     if( limdrh ) then!! PART 1:  x-interpolation!! Loop over fields.! ..x interpolation at each height needed for z interpolation.! ...x interpolation at each latitude needed for y interpolation.!        do k=1,plev           do i=1,nlon              ii1 = idp(i,k,1)              ii2 = idp(i,k,2)              ii3 = idp(i,k,3)              ii4 = idp(i,k,4)              jj = jdp(i,k)              kk = kkdp(i,k)!! Height level 1:  Linear interpolation on inner two latitudes only!!!!           fintx(i,k,1,1) = not used              fintx(i,k,2,1) = fb (ii2  ,kk-1,jj  )*xl (i,k,2) &                             + fb (ii2+1,kk-1,jj  )*xr (i,k,2)              fintx(i,k,3,1) = fb (ii3  ,kk-1,jj+1)*xl (i,k,3) &                             + fb (ii3+1,kk-1,jj+1)*xr (i,k,3)!!!           fintx(i,k,4,1) = not used!! Height level 2!!   Latitude 1:  Linear interpolation!              fintx(i,k,1,2) = fb (ii1  ,kk,jj-1)*xl (i,k,1) &                             + fb (ii1+1,kk,jj-1)*xr (i,k,1)!!   Latitude 2:  Cubic interpolation!              fxl = (   - 2.*fb (ii2-1,kk,jj) &                        - 3.*fb (ii2  ,kk,jj) &                        + 6.*fb (ii2+1,kk,jj) &                        -    fb (ii2+2,kk,jj) )*rdx6(jj)              fxr = (        fb (ii2-1,kk,jj) &                        - 6.*fb (ii2  ,kk,jj) &                        + 3.*fb (ii2+1,kk,jj) &                        + 2.*fb (ii2+2,kk,jj) )*rdx6(jj)!              deli = (       fb (ii2+1,kk,jj) - &                             fb (ii2  ,kk,jj) )*rdx(jj)              tmp1 = fac*deli              tmp2 = abs( tmp1 )              if( deli*fxl   .le. 0.0  ) fxl = 0.              if( deli*fxr   .le. 0.0  ) fxr = 0.              if( abs( fxl ) .gt. tmp2 ) fxl = tmp1              if( abs( fxr ) .gt. tmp2 ) fxr = tmp1!              fintx(i,k,2,2) = fb (ii2  ,kk,jj)*hl (i,k,2) &                             + fb (ii2+1,kk,jj)*hr (i,k,2) &                             + fxl*dhl(i,k,2) + fxr*dhr(i,k,2)!!   Latitude 3:  Cubic interpolation!              fxl = (   - 2.*fb (ii3-1,kk  ,jj+1) &                        - 3.*fb (ii3  ,kk  ,jj+1) &                        + 6.*fb (ii3+1,kk  ,jj+1) &                        -    fb (ii3+2,kk  ,jj+1) )*rdx6(jj+1)              fxr = (        fb (ii3-1,kk  ,jj+1) &                        - 6.*fb (ii3  ,kk  ,jj+1) &                        + 3.*fb (ii3+1,kk  ,jj+1) &                        + 2.*fb (ii3+2,kk  ,jj+1) )*rdx6(jj+1)!              deli = (       fb (ii3+1,kk  ,jj+1) - &                             fb (ii3  ,kk  ,jj+1) )*rdx(jj+1)              tmp1 = fac*deli              tmp2 = abs( tmp1 )              if( deli*fxl   .le. 0.0  ) fxl = 0.              if( deli*fxr   .le. 0.0  ) fxr = 0.              if( abs( fxl ) .gt. tmp2 ) fxl = tmp1              if( abs( fxr ) .gt. tmp2 ) fxr = tmp1!              fintx(i,k,3,2) = fb (ii3  ,kk  ,jj+1)*hl (i,k,3) &                             + fb (ii3+1,kk  ,jj+1)*hr (i,k,3) &                             + fxl*dhl(i,k,3) + fxr*dhr(i,k,3)!!   Latitude 4:  Linear interpolation!              fintx(i,k,4,2) = fb (ii4  ,kk,jj+2)*xl (i,k,4) &                             + fb (ii4+1,kk,jj+2)*xr (i,k,4)!! Height level 3!!   Latitude 1:  Linear interpolation!              fintx(i,k,1,3) = fb (ii1  ,kk+1,jj-1)*xl (i,k,1) &                             + fb (ii1+1,kk+1,jj-1)*xr (i,k,1)!!   Latitude 2:  Cubic interpolation!              fxl = (   - 2.*fb (ii2-1,kk+1,jj  ) &                        - 3.*fb (ii2  ,kk+1,jj  ) &                        + 6.*fb (ii2+1,kk+1,jj  ) &                        -    fb (ii2+2,kk+1,jj  ) )*rdx6(jj)              fxr = (        fb (ii2-1,kk+1,jj  ) &                        - 6.*fb (ii2  ,kk+1,jj  ) &                        + 3.*fb (ii2+1,kk+1,jj  ) &                        + 2.*fb (ii2+2,kk+1,jj  ) )*rdx6(jj)!              deli = (       fb (ii2+1,kk+1,jj  ) - &                             fb (ii2  ,kk+1,jj  ) )*rdx(jj)              tmp1 = fac*deli              tmp2 = abs( tmp1 )              if( deli*fxl   .le. 0.0  ) fxl = 0.              if( deli*fxr   .le. 0.0  ) fxr = 0.              if( abs( fxl ) .gt. tmp2 ) fxl = tmp1              if( abs( fxr ) .gt. tmp2 ) fxr = tmp1!              fintx(i,k,2,3) = fb (ii2  ,kk+1,jj  )*hl (i,k,2) &                             + fb (ii2+1,kk+1,jj  )*hr (i,k,2) &                             + fxl*dhl(i,k,2) + fxr*dhr(i,k,2)