!!   Latitude 3:  Cubic interpolation!              fxl = (   - 2.*fb (ii3-1,kk+1,jj+1) &                        - 3.*fb (ii3  ,kk+1,jj+1) &                        + 6.*fb (ii3+1,kk+1,jj+1) &                        -    fb (ii3+2,kk+1,jj+1) )*rdx6(jj+1)              fxr = (        fb (ii3-1,kk+1,jj+1) &                        - 6.*fb (ii3  ,kk+1,jj+1) &                        + 3.*fb (ii3+1,kk+1,jj+1) &                        + 2.*fb (ii3+2,kk+1,jj+1) )*rdx6(jj+1)!              deli = (       fb (ii3+1,kk+1,jj+1) - &                             fb (ii3  ,kk+1,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,3) = fb (ii3  ,kk+1,jj+1)*hl (i,k,3) &                             + fb (ii3+1,kk+1,jj+1)*hr (i,k,3) &                             + fxl*dhl(i,k,3) + fxr*dhr(i,k,3)!!   Latitude 4:  Linear interpolation!              fintx(i,k,4,3) = fb (ii4  ,kk+1,jj+2)*xl (i,k,4) &                             + fb (ii4+1,kk+1,jj+2)*xr (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  ,kk+2,jj  )*xl (i,k,2) &                             + fb (ii2+1,kk+2,jj  )*xr (i,k,2)              fintx(i,k,3,4) = fb (ii3  ,kk+2,jj+1)*xl (i,k,3) &                             + fb (ii3+1,kk+2,jj+1)*xr (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):!!   Latitude 1:  Linear interpolation!                 fintx(i,k,1,2) = fb (ii1  ,1,jj-1)*xl (i,k,1) &                                + fb (ii1+1,1,jj-1)*xr (i,k,1)!!   Latitude 2:  Cubic interpolation!                 fxl = (   - 2.*fb (ii2-1,1,jj  ) &                           - 3.*fb (ii2  ,1,jj  ) &                           + 6.*fb (ii2+1,1,jj  ) &                           -    fb (ii2+2,1,jj  ) )*rdx6(jj)                 fxr = (        fb (ii2-1,1,jj  ) &                           - 6.*fb (ii2  ,1,jj  ) &                           + 3.*fb (ii2+1,1,jj  ) &                           + 2.*fb (ii2+2,1,jj  ) )*rdx6(jj)!                                                       deli = (       fb (ii2+1,1,jj  ) - &                                fb (ii2  ,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,2) = fb (ii2  ,1,jj  )*hl (i,k,2) &                                + fb (ii2+1,1,jj  )*hr (i,k,2) &                                + fxl*dhl(i,k,2) + fxr*dhr(i,k,2)!!   Latitude 3:  Cubic interpolation!                 fxl = (   - 2.*fb (ii3-1,1,jj+1) &                           - 3.*fb (ii3  ,1,jj+1) &                           + 6.*fb (ii3+1,1,jj+1) &                           -    fb (ii3+2,1,jj+1) )*rdx6(jj+1)                 fxr = (        fb (ii3-1,1,jj+1) &                           - 6.*fb (ii3  ,1,jj+1) &                           + 3.*fb (ii3+1,1,jj+1) &                           + 2.*fb (ii3+2,1,jj+1) )*rdx6(jj+1)!                                                       deli = (       fb (ii3+1,1,jj+1) - &                                fb (ii3  ,1,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  ,1,jj+1)*hl (i,k,3) &                                + fb (ii3+1,1,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  ,1,jj+2)*xl (i,k,4) &                                + fb (ii4+1,1,jj+2)*xr (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  ,3,jj  )*xl (i,k,2) &                                + fb (ii2+1,3,jj  )*xr (i,k,2)                 fintx(i,k,3,4) = fb (ii3  ,3,jj+1)*xl (i,k,3) &                                + fb (ii3+1,3,jj+1)*xr (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  ,kdimm2,jj  )*xl (i,k,2) &                                + fb (ii2+1,kdimm2,jj  )*xr (i,k,2)                 fintx(i,k,3,1) = fb (ii3  ,kdimm2,jj+1)*xl (i,k,3) &                                + fb (ii3+1,kdimm2,jj+1)*xr (i,k,3)!!!              fintx(i,k,4,1) =  not used!! Height level 4 (placed in level 3 of stencil):!!   Latitude 1:  Linear interpolation!                 fintx(i,k,1,3) = fb (ii1  ,kdim,jj-1)*xl (i,k,1) &                                + fb (ii1+1,kdim,jj-1)*xr (i,k,1)!!   Latitude 2:  Cubic interpolation!                 fxl = (   - 2.*fb (ii2-1,kdim,jj  ) &                           - 3.*fb (ii2  ,kdim,jj  ) &                           + 6.*fb (ii2+1,kdim,jj  ) &                           -    fb (ii2+2,kdim,jj  ) )*rdx6(jj)                 fxr = (        fb (ii2-1,kdim,jj  ) &                           - 6.*fb (ii2  ,kdim,jj  ) &                           + 3.*fb (ii2+1,kdim,jj  ) &                           + 2.*fb (ii2+2,kdim,jj  ) )*rdx6(jj)!                                                           deli = (       fb (ii2+1,kdim,jj  ) - &                                fb (ii2  ,kdim,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  ,kdim,jj  )*hl (i,k,2) &                                + fb (ii2+1,kdim,jj  )*hr (i,k,2) &                                + fxl*dhl(i,k,2) + fxr*dhr(i,k,2)!!   Latitude 3:  Cubic interpolation!                 fxl = (   - 2.*fb (ii3-1,kdim,jj+1) &                           - 3.*fb (ii3  ,kdim,jj+1) &                           + 6.*fb (ii3+1,kdim,jj+1) &                           -    fb (ii3+2,kdim,jj+1) )*rdx6(jj+1)                 fxr = (        fb (ii3-1,kdim,jj+1) &                           - 6.*fb (ii3  ,kdim,jj+1) &                           + 3.*fb (ii3+1,kdim,jj+1) &                           + 2.*fb (ii3+2,kdim,jj+1) )*rdx6(jj+1)!                                                           deli = (       fb (ii3+1,kdim,jj+1) - &                                fb (ii3  ,kdim,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,3) = fb (ii3  ,kdim,jj+1)*hl (i,k,3) &                                + fb (ii3+1,kdim,jj+1)*hr (i,k,3) &                                + fxl*dhl(i,k,3) + fxr*dhr(i,k,3)!!   Latitude 4:  Linear interpolation!                 fintx(i,k,4,3) = fb (ii4  ,kdim,jj+2)*xl (i,k,4) &                                + fb (ii4+1,kdim,jj+2)*xr (i,k,4)              end if           end do        end do!! PART 2:  y-derivatives!        jmin =  1000000        jmax = -1000000        do k=1,plev           do i=1,nlon              if(jdp(i,k) .lt. jmin) jmin = jdp(i,k)              if(jdp(i,k) .gt. jmax) jmax = jdp(i,k)           end do        end do!! Loop over departure latitudes!        icount = 0        do jdpval = jmin,jmax           do k=1,plev              call wheneq(nlon    ,jdp(1,k),1       ,jdpval  , &                          indx    ,nval    )              icount = icount + nval!! y derivatives at the inner height levels (kk = 2,3) needed for! z-interpolation!              do kk  = 2,3                 do ii = 1,nval                    i = indx(ii)                    fbot(i,k,kk) = lbasdy(1,1,jdpval)*fintx(i,k,1,kk) &                                 + lbasdy(2,1,jdpval)*fintx(i,k,2,kk) &                                 + lbasdy(3,1,jdpval)*fintx(i,k,3,kk) &                                 + lbasdy(4,1,jdpval)*fintx(i,k,4,kk)                    ftop(i,k,kk) = lbasdy(1,2,jdpval)*fintx(i,k,1,kk) &                                 + lbasdy(2,2,jdpval)*fintx(i,k,2,kk) &                                 + lbasdy(3,2,jdpval)*fintx(i,k,3,kk) &                                 + lbasdy(4,2,jdpval)*fintx(i,k,4,kk)                 end do              end do           end do        end do        if (icount.ne.nlon*plev) then           write(*,*)'SLTINT:  Did not complete computations for all departure points'           call endrun        end if!! Apply SCM0 limiter to derivative estimates.!        do kk  = 2,3           do k=1,plev              do i=1,nlon                 deli = ( fintx(i,k,3,kk) - fintx(i,k,2,kk) )*rdphi(i,k)                 tmp1 = fac*deli                 tmp2 = abs( tmp1 )                 if( deli*fbot(i,k,kk)   .le. 0.0  ) fbot(i,k,kk) = 0.                 if( deli*ftop(i,k,kk)   .le. 0.0  ) ftop(i,k,kk) = 0.                 if( abs( fbot(i,k,kk) ) .gt. tmp2 ) fbot(i,k,kk) = tmp1                 if( abs( ftop(i,k,kk) ) .gt. tmp2 ) ftop(i,k,kk) = tmp1              end do           end do        end do!! PART 3:  y-interpolants!        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,2,2)*hs (i,k) + fbot (i,k  ,2)*dhs(i,k) &                           + fintx(i,k,3,2)*hn (i,k) + ftop (i,k  ,2)*dhn(i,k)              finty(i,k,3) = fintx(i,k,2,3)*hs (i,k) + fbot (i,k  ,3)*dhs(i,k) &                           + fintx(i,k,3,3)*hn (i,k) + ftop (i,k  ,3)*dhn(i,k)              finty(i,k,4) = fintx(i,k,2,4)*ys (i,k) &