lint4u Subroutine

    subroutine lint4u(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,&
                      inmx,jnmy,knmz,lnmt,isubx,jsuby,ksubz,lsubt)

    implicit none

    integer  :: nx,ny,nz,nt,mx,my,mz,mt
    integer  :: mxmy,mxmymz
    real(wp) :: p(nx,ny,nz,nt),q(mxmymz,mt)
    integer  :: inmx,jnmy,knmz,lnmt,isubx,jsuby,ksubz,lsubt
    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)
    integer  :: lt(mt),kz(mz),jy(my),ix(mx),intpol(4)
    integer  :: l,ll,lsave,iijjkk

    if (intpol(3) == 1) then
        ! linear in z
        if (lsubt == 1) then
            ! mt grid is subset of nt grid
            do ll=1,mt
                l = lnmt*(ll-1)+1
                call lint3u(nx,ny,nz,p(1,1,1,l),mx,my,mxmy,mz,q(1,ll),intpol,kz,&
                            dz,pk,pkp,jy,dym,dy,dyp,dypp,pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp,&
                            inmx,jnmy,knmz,isubx,jsuby,ksubz)
            end do
            return
        end if

        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 lint3u(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,&
                            inmx,jnmy,knmz,isubx,jsuby,ksubz)
            else
                ! interpolate l,l+1 in pt,ptp on xx,yy,zz mesh
                call lint3u(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,&
                            inmx,jnmy,knmz,isubx,jsuby,ksubz)
                call lint3u(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,&
                            inmx,jnmy,knmz,isubx,jsuby,ksubz)
            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

        if (lsubt == 1) then
            ! mt grid is subset of nt grid
            do ll=1,mt
                l = lnmt*(ll-1)+1
                call cubt3u(nx,ny,nz,p(1,1,1,l),mx,my,mxmy,mz,q(1,ll),intpol,&
                    kz,dzm,dz,dzp,dzpp,pkm,pk,pkp,pkpp,jy,dym,dy,dyp,dypp,&
                    pjm,pj,pjp,pjpp,ix,dxm,dx,dxp,dxpp,&
                    inmx,jnmy,knmz,isubx,jsuby,ksubz)
            end do
            return
        end if

        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 cubt3u(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,&
                    inmx,jnmy,knmz,isubx,jsuby,ksubz)
            else
                ! interpolate l,l+1 in pt,ptp on xx,yy,zz mesh
                call cubt3u(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,&
                    inmx,jnmy,knmz,isubx,jsuby,ksubz)
                call cubt3u(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,&
                    inmx,jnmy,knmz,isubx,jsuby,ksubz)
            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 lint4u