lint4 Subroutine

    subroutine lint4(nx,ny,nz,nt,p,mx,my,mz,mt,mxmy,mxmymz,q,intpol,&
                     lt,dt,pt,ptp,kz,dzm,dz,dzp,dzpp,pkm,pk,pkp,pkpp,&
                     jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)

    implicit none

    integer  :: nx,ny,nz,nt,mx,my,mz,mt,mxmy,mxmymz,lsave,ll,l,iijjkk
    integer  :: lt(mt),kz(mz),jy(my),ix(mx),intpol(4)
    real(wp) :: p(nx,ny,nz,nt),q(mxmymz,mt)
    real(wp) :: dt(mt),pt(mxmymz),ptp(mxmymz)
    real(wp) :: dzm(mz),dz(mz),dzp(mz),dzpp(mz)
    real(wp) :: pkm(mxmy),pk(mxmy),pkp(mxmy),pkpp(mxmy)
    real(wp) :: dym(my),dy(my),dyp(my),dypp(my)
    real(wp) :: pjm(mx),pj(mx),pjp(mx),pjpp(mx)
    real(wp) :: dxm(mx),dx(mx),dxp(mx),dxpp(mx)

    if (intpol(3) == 1) then

        ! linear in z

        lsave = -1
        do ll=1,mt
            l = lt(ll)
            if (l==lsave) then
                ! l pointer has not moved since last pass (no updates or interpolation)
            else if (l==lsave+1) then
                ! update l and interpolate l+1
                do iijjkk=1,mxmymz
                    pt(iijjkk) = ptp(iijjkk)
                end do
                call lint3(nx,ny,nz,p(1,1,1,l+1),mx,my,mxmy,mz,ptp,intpol,&
                            kz,dz,pk,pkp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,&
                            dxm,dx,dxp,dxpp)
            else
                ! interpolate l,l+1 in pt,ptp on xx,yy,zz mesh
                call lint3(nx,ny,nz,p(1,1,1,l),mx,my,mxmy,mz,pt,intpol,kz,&
                            dz,pk,pkp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
                call lint3(nx,ny,nz,p(1,1,1,l+1),mx,my,mxmy,mz,ptp,intpol,kz,&
                            dz,pk,pkp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
            end if

            ! save l pointer for next pass

            lsave = l

            ! linearly interpolate q(ii,jj,,kk,ll) from pt,ptp in t direction

            do iijjkk=1,mxmymz
                q(iijjkk,ll) = pt(iijjkk)+dt(ll)*(ptp(iijjkk)-pt(iijjkk))
            end do
        end do

    else

        ! cubic in z

        lsave = -1
        do ll=1,mt
            l = lt(ll)
            if (l==lsave) then
                ! l pointer has not moved since last pass (no updates or interpolation)
            else if (l==lsave+1) then
                ! update l and interpolate l+1
                do iijjkk=1,mxmymz
                    pt(iijjkk) = ptp(iijjkk)
                end do
                call cubt3(nx,ny,nt,p(1,1,1,l+1),mx,my,mxmy,mz,ptp,intpol,&
                            kz,dzm,dz,dzp,dzpp,pkm,pk,pkp,pkpp,&
                            jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
            else
                ! interpolate l,l+1 in pt,ptp on xx,yy,zz mesh
                call cubt3(nx,ny,nt,p(1,1,1,l),mx,my,mxmy,mz,pt,intpol,&
                            kz,dzm,dz,dzp,dzpp,pkm,pk,pkp,pkpp,&
                            jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
                call cubt3(nx,ny,nt,p(1,1,1,l+1),mx,my,mxmy,mz,ptp,intpol,&
                            kz,dzm,dz,dzp,dzpp,pkm,pk,pkp,pkpp,&
                            jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp)
            end if

            ! save l pointer for next pass

            lsave = l

            ! linearly interpolate q(ii,jj,kk,ll) from pt,ptp in t direction

            do iijjkk=1,mxmymz
                q(iijjkk,ll) = pt(iijjkk)+dt(ll)*(ptp(iijjkk)-pt(iijjkk))
            end do

        end do

    end if

    end subroutine lint4