+ fintx(i,k,3,4)*yn (i,k)           end do        end do     endif!!-----------------------------------------------------------------------!-----------------------------------------------------------------------!!    20XX loops: Lagrange cubic/linear interpolation in the horizontal!!-----------------------------------------------------------------------!-----------------------------------------------------------------------!     if( .not. 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+1,kk-1,jj  )*xr   (i,k,2) &                             + fb (ii2  ,kk-1,jj  )*xl   (i,k,2)              fintx(i,k,3,1) = fb (ii3+1,kk-1,jj+1)*xr   (i,k,3) &                             + fb (ii3  ,kk-1,jj+1)*xl   (i,k,3)!!!           fintx(i,k,4,1) = not used!! Height level 2:  Linear interpolation on outer two latitudes;!                  Cubic  interpolation on inner two latitudes.!              fintx(i,k,1,2) = fb (ii1+1,kk  ,jj-1)*xr   (i,k,1) &                             + fb (ii1  ,kk  ,jj-1)*xl   (i,k,1)              fintx(i,k,2,2) = fb (ii2-1,kk  ,jj  )*wgt1x(i,k,2) &                             + fb (ii2  ,kk  ,jj  )*wgt2x(i,k,2) &                             + fb (ii2+1,kk  ,jj  )*wgt3x(i,k,2) &                             + fb (ii2+2,kk  ,jj  )*wgt4x(i,k,2)              fintx(i,k,3,2) = fb (ii3-1,kk  ,jj+1)*wgt1x(i,k,3) &                             + fb (ii3  ,kk  ,jj+1)*wgt2x(i,k,3) &                             + fb (ii3+1,kk  ,jj+1)*wgt3x(i,k,3) &                             + fb (ii3+2,kk  ,jj+1)*wgt4x(i,k,3)              fintx(i,k,4,2) = fb (ii4+1,kk  ,jj+2)*xr   (i,k,4) &                             + fb (ii4  ,kk  ,jj+2)*xl   (i,k,4)!! Height level 3:  Linear interpolation on outer two latitudes;!                  Cubic  interpolation on inner two latitudes.!              fintx(i,k,1,3) = fb (ii1+1,kk+1,jj-1)*xr   (i,k,1) &                             + fb (ii1  ,kk+1,jj-1)*xl   (i,k,1)              fintx(i,k,2,3) = fb (ii2-1,kk+1,jj  )*wgt1x(i,k,2) &                             + fb (ii2  ,kk+1,jj  )*wgt2x(i,k,2) &                             + fb (ii2+1,kk+1,jj  )*wgt3x(i,k,2) &                             + fb (ii2+2,kk+1,jj  )*wgt4x(i,k,2)              fintx(i,k,3,3) = fb (ii3-1,kk+1,jj+1)*wgt1x(i,k,3) &                             + fb (ii3  ,kk+1,jj+1)*wgt2x(i,k,3) &                             + fb (ii3+1,kk+1,jj+1)*wgt3x(i,k,3) &                             + fb (ii3+2,kk+1,jj+1)*wgt4x(i,k,3)              fintx(i,k,4,3) = fb (ii4+1,kk+1,jj+2)*xr   (i,k,4) &                             + fb (ii4  ,kk+1,jj+2)*xl   (i,k,4)!! Height level 4:  Linear interpolation on inner two latitudes only!!!!           fintx(i,k,1,4) = not used              fintx(i,k,2,4) = fb (ii2+1,kk+2,jj  )*xr   (i,k,2) &                             + fb (ii2  ,kk+2,jj  )*xl   (i,k,2)              fintx(i,k,3,4) = fb (ii3+1,kk+2,jj+1)*xr   (i,k,3) &                             + fb (ii3  ,kk+2,jj+1)*xl   (i,k,3)!!!           fintx(i,k,4,4) = not used           end do        end do!! The following loop computes x-derivatives for those cases when the! departure point lies in either the top or bottom interval of the ! model grid.  In this special case, data are shifted up or down to! keep the departure point in the middle interval of the 4-point! stencil.  Therefore, some derivatives that were computed above will ! be over-written.!        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)!! TOP interval!              if(kdp (i,k) .eq. 1) then!! shift levels 4 and 2 data to levels 1 and 3, respectively!                 fintx(i,k,2,1) = fintx(i,k,2,4)                 fintx(i,k,3,1) = fintx(i,k,3,4)!                 fintx(i,k,1,3) = fintx(i,k,1,2)                 fintx(i,k,2,3) = fintx(i,k,2,2)                 fintx(i,k,3,3) = fintx(i,k,3,2)                 fintx(i,k,4,3) = fintx(i,k,4,2)!! Height level 1 (placed in level 2 of stencil):!  Linear interpolation on outer two latitudes;!  Cubic  interpolation on inner two latitudes.!                 fintx(i,k,1,2) = fb (ii1+1,1,jj-1)*xr   (i,k,1) &                                + fb (ii1  ,1,jj-1)*xl   (i,k,1)                 fintx(i,k,2,2) = fb (ii2-1,1,jj  )*wgt1x(i,k,2) &                                + fb (ii2  ,1,jj  )*wgt2x(i,k,2) &                                + fb (ii2+1,1,jj  )*wgt3x(i,k,2) &                                + fb (ii2+2,1,jj  )*wgt4x(i,k,2)                 fintx(i,k,3,2) = fb (ii3-1,1,jj+1)*wgt1x(i,k,3) &                                + fb (ii3  ,1,jj+1)*wgt2x(i,k,3) &                                + fb (ii3+1,1,jj+1)*wgt3x(i,k,3) &                                + fb (ii3+2,1,jj+1)*wgt4x(i,k,3)                 fintx(i,k,4,2) = fb (ii4+1,1,jj+2)*xr   (i,k,4) &                                + fb (ii4  ,1,jj+2)*xl   (i,k,4)!! Height level 3 (placed in level 4 of stencil):!  Linear interpolation on inner two latitudes only!!!!              fintx(i,k,1,4) = not used                 fintx(i,k,2,4) = fb (ii2+1,3,jj  )*xr   (i,k,2) &                                + fb (ii2  ,3,jj  )*xl   (i,k,2)                 fintx(i,k,3,4) = fb (ii3+1,3,jj+1)*xr   (i,k,3) &                                + fb (ii3  ,3,jj+1)*xl   (i,k,3)!!!              fintx(i,k,4,4) = not used!! BOT interval!              else if(kdp (i,k) .eq. kdimm1) then!! shift levels 1 and 3 data to levels 4 and 2, respectively!                 fintx(i,k,2,4) = fintx(i,k,2,1)                 fintx(i,k,3,4) = fintx(i,k,3,1)!                 fintx(i,k,1,2) = fintx(i,k,1,3)                 fintx(i,k,2,2) = fintx(i,k,2,3)                 fintx(i,k,3,2) = fintx(i,k,3,3)                 fintx(i,k,4,2) = fintx(i,k,4,3)!! Height level 2 (placed in level 1 of stencil):!  Linear interpolation on inner two latitudes only!!!!              fintx(i,k,1,1) = not used                 fintx(i,k,2,1) = fb (ii2+1,kdimm2,jj  )*xr   (i,k,2) &                                + fb (ii2  ,kdimm2,jj  )*xl   (i,k,2)                 fintx(i,k,3,1) = fb (ii3+1,kdimm2,jj+1)*xr   (i,k,3) &                                + fb (ii3  ,kdimm2,jj+1)*xl   (i,k,3)!!!              fintx(i,k,4,1) = not used!! Height level 4 (placed in level 3 of stencil):!  Linear interpolation on outer two latitudes;!  Cubic  interpolation on inner two latitudes.!                 fintx(i,k,1,3) = fb (ii1+1,kdim,jj-1)*xr   (i,k,1) &                                + fb (ii1  ,kdim,jj-1)*xl   (i,k,1)                 fintx(i,k,2,3) = fb (ii2-1,kdim,jj  )*wgt1x(i,k,2) &                                + fb (ii2  ,kdim,jj  )*wgt2x(i,k,2) &                                + fb (ii2+1,kdim,jj  )*wgt3x(i,k,2) &                                + fb (ii2+2,kdim,jj  )*wgt4x(i,k,2)                 fintx(i,k,3,3) = fb (ii3-1,kdim,jj+1)*wgt1x(i,k,3) &                                + fb (ii3  ,kdim,jj+1)*wgt2x(i,k,3) &                                + fb (ii3+1,kdim,jj+1)*wgt3x(i,k,3) &                                + fb (ii3+2,kdim,jj+1)*wgt4x(i,k,3)                 fintx(i,k,4,3) = fb (ii4+1,kdim,jj+2)*xr   (i,k,4) &                                + fb (ii4  ,kdim,jj+2)*xl   (i,k,4)              end if           end do        end do!! PART 2:  Y-INTERPOLATION!! Linear on outside of stencil; Lagrange cubic on inside.!        do k=1,plev           do i=1,nlon              finty(i,k,1) = fintx(i,k,2,1)*ys   (i,k) &                           + fintx(i,k,3,1)*yn   (i,k)              finty(i,k,2) = fintx(i,k,1,2)*wgt1y(i,k) &                           + fintx(i,k,2,2)*wgt2y(i,k) &                           + fintx(i,k,3,2)*wgt3y(i,k) &                           + fintx(i,k,4,2)*wgt4y(i,k)              finty(i,k,3) = fintx(i,k,1,3)*wgt1y(i,k) &                           + fintx(i,k,2,3)*wgt2y(i,k) &                           + fintx(i,k,3,3)*wgt3y(i,k) &                           + fintx(i,k,4,3)*wgt4y(i,k)              finty(i,k,4) = fintx(i,k,2,4)*ys   (i,k) &                           + fintx(i,k,3,4)*yn   (i,k)           end do        end do     endif!!-----------------------------------------------------------------------!-----------------------------------------------------------------------!!    30XX loops: Hermite  cubic/linear interpolation in the vertical!!-----------------------------------------------------------------------!-----------------------------------------------------------------------!     if( limdrv ) then        icount = 0        do kdpval = 1,kdimm1           do k=1,plev              call wheneq(nlon    ,kdp(1,k),1       ,kdpval  , &                          indx    ,nval    )              icount = icount + nval              do ii = 1,nval                 i = indx(ii)                 if(kdpval .eq. 1) then                    ftop(i,k,1) = ( finty(i,k,3) - finty(i,k,2) )*rdz(i,k)                    fbot(i,k,1) = wdz(1,1,      2)*finty(i,k,2) + &                                  wdz(2,1,      2)*finty(i,k,3) + &                                  wdz(3,1,      2)*finty(i,k,4) + &                                  wdz(4,1,      2)*finty(i,k,1)                 else if(kdpval .eq. kdimm1) then                    ftop(i,k,1) = wdz(1,2,kdimm2 )*finty(i,k,4) + &                                  wdz(2,2,kdimm2 )*finty(i,k,1) + &                                  wdz(3,2,kdimm2 )*finty(i,k,2) + &                                  wdz(4,2,kdimm2 )*finty(i,k,3)                    fbot(i,k,1) = 0.                 else                    ftop(i,k,1) = wdz(1,1,kdpval )*finty(i,k,1) + &                                  wdz(2,1,kdpval )*finty(i,k,2) + &                                  wdz(3,1,kdpval )*finty(i,k,3) + &                                  wdz(4,1,kdpval )*finty(i,k,4)                    fbot(i,k,1) = wdz(1,2,kdpval )*finty(i,k,1) + &                                  wdz(2,2,kdpval )*finty(i,k,2) + &                                  wdz(3,2,kdpval )*finty(i,k,3) + &                                  wdz(4,2,kdpval )*finty(i,k,4)                 endif              end do           end do        end do!! Apply SCM0 limiter to derivative estimates.!        do k=1,plev           do i=1,nlon              deli = ( finty(i,k,3) - finty(i,k,2) )*rdz(i,k)              tmp1 = fac*deli              tmp2 = abs( tmp1 )              if( deli*fbot(i,k,1)   .le. 0.0  ) fbot(i,k,1) = 0.              if( deli*ftop(i,k,1)   .le. 0.0  ) ftop(i,k,1) = 0.              if( abs( fbot(i,k,1) ) .gt. tmp2 ) fbot(i,k,1) = tmp1              if( abs( ftop(i,k,1) ) .gt. tmp2 ) ftop(i,k,1) = tmp1              fdp(i,k) = finty(i,k,2)*ht(i,k) + ftop(i,k,1)*dht(i,k) + &                         finty(i,k,3)*hb(i,k) + fbot(i,k,1)*dhb(i,k)           end do        end do        if (icount.ne.nlon*plev) then           write(6,*)'SLTINT:  Did not complete computations for all departure points'           call endrun        endif     endif!!-----------------------------------------------------------------------!-----------------------------------------------------------------------!!    40XX loops: Lagrange cubic/linear interpolation in the vertical!!-----------------------------------------------------------------------!-----------------------------------------------------------------------!     if( .not. limdrv ) then        do k=1,plev           do i=1,nlon              fdp(i,k) = finty(i,k,1)*wgt1z(i,k) &                       + finty(i,k,2)*wgt2z(i,k) &                       + finty(i,k,3)*wgt3z(i,k) &                       + finty(i,k,4)*wgt4z(i,k)           end do        end do!! IF the departure point is in either the top or bottom interval of the! model grid:  THEN! Perform Hermite cubic interpolation.  The data are shifted up or down! such that the departure point sits in the middle interval of the! 4 point stencil (the shift originally took place in routine "LAGXIN").! Therefore the derivative weights must be applied appropriately to! account for this shift.  The following overwrites some results from! the previous loop.!        do k=1,plev           do i=1,nlon