sltint.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

              if(kdp (i,k) .eq. 1) then
                 tmptop = ( finty(i,k,3) - finty(i,k,2) )*rdz(i,k)
                 tmpbot = 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)
                 fdp(i,k) = finty(i,k,2)*ht (i,k) + tmptop      *dht(i,k) &
                          + finty(i,k,3)*hb (i,k) + tmpbot      *dhb(i,k)
              else if(kdp (i,k) .eq. kdimm1) then
                 tmptop = 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)
!!!!!            tmpbot = 0.
                 fdp(i,k) = finty(i,k,2)*ht (i,k) + tmptop       *dht(i,k) &
                          + finty(i,k,3)*hb (i,k)
!!!!!                     + tmpbot      *dhb(i,k)
              end if
           end do
        end do
     end if
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
!  Horizontal interpolation only
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
  elseif(lhrzint) then
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
!    50XX loops: an optimized Lagrange cubic/linear algorithm (no
!                Hermite interpolator available)
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
     if(limdrh) then
!
! PART 1:  x-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 = kdp(i,k)
!
! 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)
           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,2
                 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,2
           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
              fdp(i,k) = 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)
           end do
        end do
     endif
!
     if( .not. limdrh ) then
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
!    60XX loops: Hermite cubic/linear interpolation in the horizontal
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
        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 = kdp(i,k)
!
! x-interpolants for the 4 latitudes
!
              f1 = fb(ii1+1,kk,jj-1)*xr   (i,k,1) &
                 + fb(ii1  ,kk,jj-1)*xl   (i,k,1)
              f2 = 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)
              f3 = 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)
              f4 = fb(ii4+1,kk,jj+2)*xr   (i,k,4) &
                 + fb(ii4  ,kk,jj+2)*xl   (i,k,4)
!
! y-interpolant
!
              fdp(i,k) = f1*wgt1y(i,k) + f2*wgt2y(i,k) + &
                         f3*wgt3y(i,k) + f4*wgt4y(i,k)
           end do
        end do
     end if
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
!  Vertical interpolation only
!
!    70XX loops: an optimized Lagrange cubic/linear algorithm (no
!                Hermite interpolator available)
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
  else if(lvrtint) then
     if(limdrh .or. limdrv) then
        write(6,*) 'SLTINT:  ERROR:  this routine does not provide '
        write(6,*) 'shape preserving capability for vertical-only'
        write(6,*) ' interpolation'
        call endrun
     end if

     do k=1,plev
        do i=1,nlon
           kk = kkdp(i,k)
           ii = i1+i-1
           fdp(i,k) = fb(ii,kk-1,jcen)*wgt1z(i,k) &
                    + fb(ii,kk  ,jcen)*wgt2z(i,k) &
                    + fb(ii,kk+1,jcen)*wgt3z(i,k) &
                    + fb(ii,kk+2,jcen)*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 following
! overwrites some results from the previous loop.
!
     do k=1,plev
        do i=1,nlon
           ii = i1+i-1
           if(kdp (i,k) .eq. 1) then
!!!!!         tmptop   = 0.
              tmpbot   = wdz(1,1,     2)*fb(ii,     1,jcen) &
                       + wdz(2,1,     2)*fb(ii,     2,jcen) &
                       + wdz(3,1,     2)*fb(ii,     3,jcen) &
                       + wdz(4,1,     2)*fb(ii,     4,jcen)
              fdp(i,k) = fb(ii,1     ,jcen)*ht (i,k) &
                    &  + fb(ii,2     ,jcen)*hb (i,k) &
                    &  + tmpbot            *dhb(i,k)
!!!!!               &  + tmptop            *dht(i,k)
           else if(kdp (i,k) .eq. kdimm1) then
              tmptop = wdz(1,2,kdimm2)*fb(ii,kdimm3,jcen) &
                     + wdz(2,2,kdimm2)*fb(ii,kdimm2,jcen) &
                     + wdz(3,2,kdimm2)*fb(ii,kdimm1,jcen) &
                     + wdz(4,2,kdimm2)*fb(ii,kdim  ,jcen)
!!!!!         tmpbot = 0.
              fdp(i,k) = fb(ii,kdimm1,jcen)*ht (i,k) &
                       + tmptop            *dht(i,k) &
                       + fb(ii,kdim  ,jcen)*hb (i,k)
!!!!!                  + tmpbot            *dhb(i,k)
           end if
        end do
     end do
  else
     write(6,*) 'SLTINT:  Error: must specify at least one of "lhr', &
          'zint" or "lvrtint" to be ".true."'
     call endrun
  end if
!
  return
end subroutine sltint