sltint.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

                           + 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