cubt2 Subroutine

    subroutine cubt2(nx,ny,p,mx,my,q,intpol,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)

    implicit none

    integer  :: nx,ny,mx,my,intpol(2),jy(my),ix(mx)
    integer  :: jsave,j,jj,ii
    real(wp) :: p(nx,ny),q(mx,my)
    real(wp) :: pjm(mx),pj(mx),pjp(mx),pjpp(mx)
    real(wp) :: dym(my),dy(my),dyp(my),dypp(my)
    real(wp) :: dxm(mx),dx(mx),dxp(mx),dxpp(mx)

    if (intpol(1)==1) then

        ! linear in x

        jsave = -3
        do jj=1,my

            ! load closest four j lines containing interpolate on xx mesh
            ! for j-1,j,j+1,j+2 in pjm,pj,pjp,pjpp

            j = jy(jj)
            if (j==jsave) then
                ! j pointer has not moved since last pass (no updates or interpolation)
            else if (j==jsave+1) then
                ! update j-1,j,j+1 and interpolate j+2
                do ii=1,mx
                    pjm(ii) = pj(ii)
                    pj(ii) = pjp(ii)
                    pjp(ii) = pjpp(ii)
                end do
                call lint1(nx,p(1,j+2),mx,pjpp,ix,dx)
            else if (j==jsave+2) then
                ! update j-1,j and interpolate j+1,j+2
                do ii=1,mx
                    pjm(ii) = pjp(ii)
                    pj(ii) = pjpp(ii)
                end do
                call lint1(nx,p(1,j+1),mx,pjp,ix,dx)
                call lint1(nx,p(1,j+2),mx,pjpp,ix,dx)
            else if (j==jsave+3) then
                ! update j-1 and interpolate j,j+1,j+2
                do ii=1,mx
                    pjm(ii) = pjpp(ii)
                end do
                call lint1(nx,p(1,j),mx,pj,ix,dx)
                call lint1(nx,p(1,j+1),mx,pjp,ix,dx)
                call lint1(nx,p(1,j+2),mx,pjpp,ix,dx)
            else
                ! interpolate all four j-1,j,j+1,j+2
                call lint1(nx,p(1,j-1),mx,pjm,ix,dx)
                call lint1(nx,p(1,j),mx,pj,ix,dx)
                call lint1(nx,p(1,j+1),mx,pjp,ix,dx)
                call lint1(nx,p(1,j+2),mx,pjpp,ix,dx)
            end if

            ! save j pointer for next pass

            jsave = j

            ! cubically interpolate q(ii,jj) from pjm,pj,pjp,pjpp in y direction

            do ii=1,mx
                q(ii,jj) = dym(jj)*pjm(ii)+dy(jj)*pj(ii)+dyp(jj)*pjp(ii)+dypp(jj)*pjpp(ii)
            end do
        end do
        return

    else

        ! cubic in x

        jsave = -3
        do jj=1,my

            ! load closest four j lines containing interpolate on xx mesh
            ! for j-1,j,j+1,j+2 in pjm,pj,pjp,pjpp

            j = jy(jj)
            if (j==jsave) then
                ! j pointer has not moved since last pass (no updates or interpolation)
            else if (j==jsave+1) then
                ! update j-1,j,j+1 and interpolate j+2
                do ii=1,mx
                    pjm(ii) = pj(ii)
                    pj(ii) = pjp(ii)
                    pjp(ii) = pjpp(ii)
                end do
                call cubt1(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp)
            else if (j==jsave+2) then
                ! update j-1,j and interpolate j+1,j+2
                do ii=1,mx
                    pjm(ii) = pjp(ii)
                    pj(ii) = pjpp(ii)
                end do
                call cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp)
                call cubt1(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp)
            else if (j==jsave+3) then
                ! update j-1 and interpolate j,j+1,j+2
                do ii=1,mx
                    pjm(ii) = pjpp(ii)
                end do
                call cubt1(nx,p(1,j),mx,pj,ix,dxm,dx,dxp,dxpp)
                call cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp)
                call cubt1(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp)
            else
                ! interpolate all four j-1,j,j+1,j+2
                call cubt1(nx,p(1,j-1),mx,pjm,ix,dxm,dx,dxp,dxpp)
                call cubt1(nx,p(1,j),mx,pj,ix,dxm,dx,dxp,dxpp)
                call cubt1(nx,p(1,j+1),mx,pjp,ix,dxm,dx,dxp,dxpp)
                call cubt1(nx,p(1,j+2),mx,pjpp,ix,dxm,dx,dxp,dxpp)
            end if

            ! save j pointer for next pass

            jsave = j

            ! cubically interpolate q(ii,jj) from pjm,pj,pjp,pjpp in y direction

            do ii=1,mx
                q(ii,jj) = dym(jj)*pjm(ii)+dy(jj)*pj(ii)+dyp(jj)*pjp(ii)+dypp(jj)*pjpp(ii)
            end do
        end do

    end if

    end subroutine cubt2