rgrd4 Subroutine

private subroutine rgrd4(nx, ny, nz, nt, x, y, z, t, p, mx, my, mz, mt, xx, yy, zz, tt, q, intpol, w, lw, iw, liw, ier)

subroutine rgrd4 interpolates the values p(i,j,k,l) on the orthogonal grid (x(i),y(j),z(k),t(l)) for i=1,...,nx;j=1,...,ny;k=1,...,nz;l=1,...,nt onto q(ii,jj,kk,ll) on the orthogonal grid (xx(ii),yy(jj),zz(kk),tt(ll)) for ii=1,...,mx;jj=1,...,my;kk=1,...,mz;ll=1,...,mt

method

linear or cubic interpolation is used (independently) in each direction (see argument intpol).

requirements

each of the x,y,z,t grids must be strictly montonically increasing and each of the xx,yy,zz,tt grids must be montonically increasing (see ier = 4). in addition the (X,Y,Z,T) region of the q grid

[xx(1),xx(mx)] X [yy(1),yy(my)] X [zz(1),zz(mz)] X [tt(1),tt(my)]

must lie within the (X,Y,Z,T) region of the p grid

[x(1),x(nx)] X [y(1),y(ny)] X [z(1),z(nz)] X [t(1),t(nt)].

extrapolation is not allowed (see ier=3). if these (X,Y,Z,T) regions are identical and the orthogonal grids are UNIFORM in each direction then subroutine rgrd4u should be used instead of rgrd4.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx

the integer dimension of the grid vector x and the first dimension of p. nx > 1 if intpol(1) = 1 or nx > 3 if intpol(1) = 3 is required.
integer, intent(in) :: ny

the integer dimension of the grid vector y and the second dimension of p. ny > 1 if intpol(2) = 1 or ny > 3 if intpol(2) = 3 is required.
integer, intent(in) :: nz

the integer dimension of the grid vector z and the third dimension of p. nz > 1 if intpol(3) = 1 or nz > 3 if intpol(3) = 3 is required.
integer, intent(in) :: nt

the integer dimension of the grid vector t and the fourth dimension of p. nt > 1 if intpol(4) = 1 or nt > 3 if intpol(4) = 3 is required.
real(kind=wp), intent(in) :: x(nx)

a real(wp) nx vector of strictly increasing values which defines the x portion of the orthogonal grid on which p is given
real(kind=wp), intent(in) :: y(ny)

a real(wp) ny vector of strictly increasing values which defines the y portion of the orthogonal grid on which p is given
real(kind=wp), intent(in) :: z(nz)

a real(wp) nz vector of strictly increasing values which defines the z portion of the orthogonal grid on which p is given
real(kind=wp), intent(in) :: t(nt)

a real(wp) nt vector of strictly increasing values which defines the t portion of the orthogonal grid on which p is given
real(kind=wp), intent(in) :: p(nx,ny,nz,nt)

a real(wp) nx by ny by nz by nt array of values given on the (x,y,z,t) grid
integer, intent(in) :: mx

the integer dimension of the grid vector xx and the first dimension of q. mx > 0 is required.
integer, intent(in) :: my

the integer dimension of the grid vector yy and the second dimension of q. my > 0 is required.
integer, intent(in) :: mz

the integer dimension of the grid vector zz and the third dimension of q. mz > 0 is required.
integer, intent(in) :: mt

the integer dimension of the grid vector tt and the fourth dimension of q. mt > 0 is required.
real(kind=wp), intent(in) :: xx(mx)

a real(wp) mx vector of increasing values which defines the x portion of the orthogonal grid on which q is defined. xx(1) < x(1) or xx(mx) > x(nx) is not allowed (see ier = 3)
real(kind=wp), intent(in) :: yy(my)

a real(wp) my vector of increasing values which defines the y portion of the orthogonal grid on which q is defined. yy(1) < y(1) or yy(my) > y(ny) is not allowed (see ier = 3)
real(kind=wp), intent(in) :: zz(mz)

a real(wp) mz vector of increasing values which defines the z portion of the orthogonal grid on which q is defined. zz(1) < z(1) or zz(mz) > z(nz) is not allowed (see ier = 3)
real(kind=wp), intent(in) :: tt(mt)

a real(wp) mt vector of increasing values which defines the t portion of the orthogonal grid on which q is defined. tt(1) < t(1) or tt(mt) > t(nt) is not allowed (see ier = 3)
real(kind=wp), intent(out) :: q(mx,my,mz,mt)

a real(wp) mx by my by mz by mt array of values on the (xx,yy,zz,tt) grid which are interpolated from p on the (x,y,z,t) grid
integer, intent(in) :: intpol(4)

an integer vector of dimension 4 which sets linear or cubic interpolation in each of the x,y,z,t directions as follows:

  • intpol(1) = 1 sets linear interpolation in the x direction
  • intpol(1) = 3 sets cubic interpolation in the x direction.
  • intpol(2) = 1 sets linear interpolation in the y direction
  • intpol(2) = 3 sets cubic interpolation in the y direction.
  • intpol(3) = 1 sets linear interpolation in the z direction
  • intpol(3) = 3 sets cubic interpolation in the z direction.
  • intpol(4) = 1 sets linear interpolation in the t direction
  • intpol(4) = 3 sets cubic interpolation in the t direction.

values other than 1 or 3 in intpol are not allowed (ier = 6).
real(kind=wp), intent(inout) :: w(lw)

a real(wp) work space of length at least lw which must be provided in the routine calling rgrd4
integer, intent(in) :: lw

the integer length of the real(wp) work space w. let

  • lwx = mx if intpol(1) = 1
  • lwx = 4*mx if intpol(1) = 3
  • lwy = my+2*mx if intpol(2) = 1
  • lwy = 4*(my+mx) if intpol(2) = 3
  • lwz = 2mxmy+mz if intpol(3) = 1
  • lwz = 4(mxmy+mz) if intpol(3) = 3
  • lwt = 2mxmy*mz+mt if intpol(4) = 1
  • lwt = 4(mxmy*mz+mt) if intpol(4) = 3

then lw must be greater than or equal to lwx+lwy+lwz+lwt
integer, intent(inout) :: iw(liw)

an integer work space of length at least liw which must be provided in the routine calling rgrd4
integer, intent(in) :: liw

the integer length of the integer work space iw. liw must be at least mx+my+mz+mt
integer, intent(out) :: ier

an integer error flag set as follows:

  • ier = 0 if no errors in input arguments are detected
  • ier = 1 if min(mx,my,mz,mt) < 1
  • ier = 2 if nx < 2 when intpol(1)=1 or nx < 4 when intpol(1)=3 (or)
  • ny < 2 when intpol(2)=1 or ny < 4 when intpol(2)=3 (or)
  • nz < 2 when intpol(3)=1 or nz < 4 when intpol(3)=3 (or)
  • nt < 2 when intpol(4)=1 or nt < 4 when intpol(4)=3
  • ier = 3 if xx(1) < x(1) or x(nx) < xx(mx) (or)
  • yy(1) < y(1) or y(ny) < yy(my) (or)
  • zz(1) < z(1) or z(nz) < zz(mz) (or)

  • tt(1) < t(1) or t(nt) < tt(mt) to avoid this flag when end points are intended to be the same but may differ slightly due to roundoff error, they should be set exactly in the calling routine (e.g., if both grids have the same y boundaries then yy(1)=y(1) and yy(my)=y(ny) should be set before calling rgrd4)

  • ier = 4 if one of the grids x,y,z,t is not strictly monotonically increasing or if one of the grids xx,yy,zz,tt is not montonically increasing. more precisely if:

  • x(i+1) <= x(i) for some i such that 1 <= i < nx (or)

  • y(j+1) <= y(j) for some j such that 1 <= j < ny (or)
  • z(k+1) <= z(k) for some k such that 1 <= k < nz (or)
  • t(l+1) <= t(l) for some l such that 1 <= l < nt (or)
  • xx(ii+1) < xx(ii) for some ii such that 1 <= ii < mx (or)
  • yy(jj+1) < yy(jj) for some jj such that 1 <= jj < my (or)
  • zz(kk+1) < zz(k) for some kk such that 1 <= kk < mz (or)
  • tt(ll+1) < tt(l) for some ll such that 1 <= ll < mt
  • ier = 5 if lw or liw is too small (insufficient work space)
  • ier = 6 if any of intpol(1),intpol(2),intpol(3),intpol(4) is not equal to 1 or 3

Calls

proc~~rgrd4~~CallsGraph proc~rgrd4 regridpack_module::rgrd4 proc~cubnmx regridpack_module::cubnmx proc~rgrd4->proc~cubnmx proc~cubt4 regridpack_module::cubt4 proc~rgrd4->proc~cubt4 proc~linmx regridpack_module::linmx proc~rgrd4->proc~linmx proc~lint4 regridpack_module::lint4 proc~rgrd4->proc~lint4 proc~cubt3 regridpack_module::cubt3 proc~cubt4->proc~cubt3 proc~lint3 regridpack_module::lint3 proc~cubt4->proc~lint3 proc~lint4->proc~cubt3 proc~lint4->proc~lint3 proc~cubt2 regridpack_module::cubt2 proc~cubt3->proc~cubt2 proc~lint2 regridpack_module::lint2 proc~cubt3->proc~lint2 proc~lint3->proc~cubt2 proc~lint3->proc~lint2 proc~cubt1 regridpack_module::cubt1 proc~cubt2->proc~cubt1 proc~lint1 regridpack_module::lint1 proc~cubt2->proc~lint1 proc~lint2->proc~cubt1 proc~lint2->proc~lint1

Called by

proc~~rgrd4~~CalledByGraph proc~rgrd4 regridpack_module::rgrd4 interface~regrid regridpack_module::regrid interface~regrid->proc~rgrd4 proc~rgrd4_wrapper regridpack_module::rgrd4_wrapper interface~regrid->proc~rgrd4_wrapper proc~rgrd4_wrapper->proc~rgrd4

Source Code

    subroutine rgrd4(nx,ny,nz,nt,x,y,z,t,p,mx,my,mz,mt,xx,yy,zz,tt,q,intpol,w,lw,iw,liw,ier)

    implicit none

    integer,intent(in)  :: nx   !! the integer dimension of the grid vector x and the first dimension of p.
                                !! nx > 1 if intpol(1) = 1 or nx > 3 if intpol(1) = 3 is required.
    integer,intent(in)  :: ny   !! the integer dimension of the grid vector y and the second dimension of p.
                                !! ny > 1 if intpol(2) = 1 or ny > 3 if intpol(2) = 3 is required.
    integer,intent(in)  :: nz   !! the integer dimension of the grid vector z and the third dimension of p.
                                !! nz > 1 if intpol(3) = 1 or nz > 3 if intpol(3) = 3 is required.
    integer,intent(in)  :: nt   !! the integer dimension of the grid vector t and the fourth dimension of p.
                                !! nt > 1 if intpol(4) = 1 or nt > 3 if intpol(4) = 3 is required.
    integer,intent(in)  :: mx   !! the integer dimension of the grid vector xx and the first dimension
                                !! of q.  mx > 0 is required.
    integer,intent(in)  :: my   !! the integer dimension of the grid vector yy and the second dimension
                                !! of q.  my > 0 is required.
    integer,intent(in)  :: mz   !! the integer dimension of the grid vector zz and the third dimension
                                !! of q.  mz > 0 is required.
    integer,intent(in)  :: mt   !! the integer dimension of the grid vector tt and the fourth dimension
                                !! of q.  mt > 0 is required.
    integer,intent(in)  :: lw   !! the integer length of the real(wp) work space w.  let
                                !!
                                !! * lwx = mx                  if intpol(1) = 1
                                !! * lwx = 4*mx                if intpol(1) = 3
                                !! * lwy = my+2*mx             if intpol(2) = 1
                                !! * lwy = 4*(my+mx)           if intpol(2) = 3
                                !! * lwz = 2*mx*my+mz          if intpol(3) = 1
                                !! * lwz = 4*(mx*my+mz)        if intpol(3) = 3
                                !! * lwt = 2*mx*my*mz+mt       if intpol(4) = 1
                                !! * lwt = 4*(mx*my*mz+mt)     if intpol(4) = 3
                                !!
                                !! then lw must be greater than or equal to lwx+lwy+lwz+lwt
    integer,intent(in)  :: liw  !! the integer length of the integer work space iw.  liw must be at least mx+my+mz+mt
    integer,intent(out)  :: ier     !! an integer error flag set as follows:
                                    !!
                                    !! * ier = 0 if no errors in input arguments are detected
                                    !! * ier = 1 if  min(mx,my,mz,mt) < 1
                                    !! * ier = 2 if nx < 2 when intpol(1)=1 or nx < 4 when intpol(1)=3 (or)
                                    !!    * ny < 2 when intpol(2)=1 or ny < 4 when intpol(2)=3 (or)
                                    !!    * nz < 2 when intpol(3)=1 or nz < 4 when intpol(3)=3 (or)
                                    !!    * nt < 2 when intpol(4)=1 or nt < 4 when intpol(4)=3
                                    !! * ier = 3 if xx(1) < x(1) or x(nx) < xx(mx) (or)
                                    !!    * yy(1) < y(1) or y(ny) < yy(my) (or)
                                    !!    * zz(1) < z(1) or z(nz) < zz(mz) (or)
                                    !!    * tt(1) < t(1) or t(nt) < tt(mt)
                                    !!   to avoid this flag when end points are intended to be the
                                    !!   same but may differ slightly due to roundoff error, they
                                    !!   should be set exactly in the calling routine (e.g., if both
                                    !!   grids have the same y boundaries then yy(1)=y(1) and yy(my)=y(ny)
                                    !!   should be set before calling rgrd4)
                                    !!
                                    !! * ier = 4 if one of the grids x,y,z,t is not strictly monotonically
                                    !!   increasing or if one of the grids xx,yy,zz,tt is not
                                    !!   montonically increasing.  more precisely if:
                                    !!
                                    !!    * x(i+1) <= x(i) for some i such that 1 <= i < nx (or)
                                    !!    * y(j+1) <= y(j) for some j such that 1 <= j < ny (or)
                                    !!    * z(k+1) <= z(k) for some k such that 1 <= k < nz (or)
                                    !!    * t(l+1) <= t(l) for some l such that 1 <= l < nt (or)
                                    !!    * xx(ii+1) < xx(ii) for some ii such that 1 <= ii < mx (or)
                                    !!    * yy(jj+1) < yy(jj) for some jj such that 1 <= jj < my (or)
                                    !!    * zz(kk+1) < zz(k)  for some kk such that 1 <= kk < mz (or)
                                    !!    * tt(ll+1) < tt(l)  for some ll such that 1 <= ll < mt
                                    !! * ier = 5 if lw or liw is too small (insufficient work space)
                                    !! * ier = 6 if any of intpol(1),intpol(2),intpol(3),intpol(4)
                                    !!                 is not equal to 1 or 3

    integer,intent(inout)  :: iw(liw)   !! an integer work space of length at least liw
                                        !! which must be provided in the routine calling rgrd4
    integer,intent(in)  :: intpol(4)    !! an integer vector of dimension 4 which sets linear or cubic
                                        !! interpolation in each of the x,y,z,t directions as follows:
                                        !!
                                        !! * intpol(1) = 1 sets linear interpolation in the x direction
                                        !! * intpol(1) = 3 sets cubic interpolation in the x direction.
                                        !! * intpol(2) = 1 sets linear interpolation in the y direction
                                        !! * intpol(2) = 3 sets cubic interpolation in the y direction.
                                        !! * intpol(3) = 1 sets linear interpolation in the z direction
                                        !! * intpol(3) = 3 sets cubic interpolation in the z direction.
                                        !! * intpol(4) = 1 sets linear interpolation in the t direction
                                        !! * intpol(4) = 3 sets cubic interpolation in the t direction.
                                        !!
                                        !! values other than 1 or 3 in intpol are not allowed (ier = 6).
    real(wp),intent(in) :: x(nx)    !! a real(wp) nx vector of strictly increasing values which defines the x
                                    !! portion of the orthogonal grid on which p is given
    real(wp),intent(in) :: y(ny)    !! a real(wp) ny vector of strictly increasing values which defines the y
                                    !! portion of the orthogonal grid on which p is given
    real(wp),intent(in) :: z(nz)    !! a real(wp) nz vector of strictly increasing values which defines the z
                                    !! portion of the orthogonal grid on which p is given
    real(wp),intent(in) :: t(nt)    !! a real(wp) nt vector of strictly increasing values which defines the t
                                    !! portion of the orthogonal grid on which p is given
    real(wp),intent(in) :: p(nx,ny,nz,nt)   !! a real(wp) nx by ny by nz by nt array of values given on the (x,y,z,t) grid
    real(wp),intent(inout) :: w(lw) !! a real(wp) work space of length at least lw
                                    !! which must be provided in the routine calling rgrd4
    real(wp),intent(in) :: xx(mx)   !! a real(wp) mx vector of increasing values which defines the x portion of the
                                    !! orthogonal grid on which q is defined.  xx(1) < x(1) or xx(mx) > x(nx)
                                    !! is not allowed (see ier = 3)
    real(wp),intent(in) :: yy(my)   !! a real(wp) my vector of increasing values which defines the y portion of the
                                    !! orthogonal grid on which q is defined.  yy(1) < y(1) or yy(my) > y(ny)
                                    !! is not allowed (see ier = 3)
    real(wp),intent(in) :: zz(mz)   !! a real(wp) mz vector of increasing values which defines the z portion of the
                                    !! orthogonal grid on which q is defined.  zz(1) < z(1) or zz(mz) > z(nz)
                                    !! is not allowed (see ier = 3)
    real(wp),intent(in) :: tt(mt)   !! a real(wp) mt vector of increasing values which defines the t portion of the
                                    !! orthogonal grid on which q is defined.  tt(1) < t(1) or tt(mt) > t(nt)
                                    !! is not allowed (see ier = 3)
    real(wp),intent(out) :: q(mx,my,mz,mt)  !! a real(wp) mx by my by mz by mt array of values on the (xx,yy,zz,tt) grid
                                            !! which are interpolated from p on the (x,y,z,t) grid

    integer  :: l2,l3,l4,l5,l6,l7,l8,l9
    integer  :: k2,k3,k4,k5,k6,k7,k8,k9
    integer  :: j2,j3,j4,j5,j6,j7,j8,j9
    integer  :: i2,i3,i4,i5
    integer  :: lwx,lwy,lwz,lwt,mxmy,mxmymz
    integer  :: ii,jj,kk,ll,i,j,k,l
    integer  :: jy,kz,lt

    ! check input arguments

    ! check (xx,yy,zz,tt) grid resolution
    ier = 1
    if (min(mx,my,mz,mt) < 1) return

    ! check intpol
    ier = 6
    if (intpol(1)/=1 .and. intpol(1)/=3) return
    if (intpol(2)/=1 .and. intpol(2)/=3) return
    if (intpol(3)/=1 .and. intpol(3)/=3) return
    if (intpol(4)/=1 .and. intpol(4)/=3) return

    ! check (x,y,z,t) grid resolution
    ier = 2
    if (intpol(1)==1 .and. nx<2) return
    if (intpol(1)==3 .and. nx<4) return
    if (intpol(2)==1 .and. ny<2) return
    if (intpol(2)==3 .and. ny<4) return
    if (intpol(3)==1 .and. nz<2) return
    if (intpol(3)==3 .and. nz<4) return
    if (intpol(4)==1 .and. nt<2) return
    if (intpol(4)==3 .and. nt<4) return

    ! check work space length input and set minimum
    ier = 5
    mxmy = mx*my
    mxmymz = mxmy*mz
    if (intpol(1)==1) then
        lwx = mx
    else
        lwx = 4*mx
    end if
    if (intpol(2)==1) then
        lwy = (my+2*mx)
    else
        lwy = 4*(mx+my)
    end if
    if (intpol(3)==1) then
        lwz = (2*mxmy+mz)
    else
        lwz = 4*(mxmy+mz)
    end if
    if (intpol(4)==1) then
        lwt = (2*mxmymz+mt)
    else
        lwt = 4*(mxmymz+mt)
    end if
    if (lw < lwx+lwy+lwz+lwt) return
    if (liw < mx+my+mz+mt) return

    ! check (xx,yy,zz,tt) grid contained in (x,y,z,t) grid
    ier = 3
    if (xx(1)<x(1) .or. xx(mx)>x(nx)) return
    if (yy(1)<y(1) .or. yy(my)>y(ny)) return
    if (zz(1)<z(1) .or. zz(mz)>z(nz)) return
    if (tt(1)<t(1) .or. tt(mt)>t(nt)) return

    ! check montonicity of grids
    ier = 4
    do i=2,nx
        if (x(i-1)>=x(i)) return
    end do
    do j=2,ny
        if (y(j-1)>=y(j)) return
    end do
    do k=2,nz
        if (z(k-1)>=z(k)) return
    end do
    do l=2,nt
        if (t(l-1)>=t(l)) return
    end do
    do ii=2,mx
        if (xx(ii-1)>xx(ii)) return
    end do
    do jj=2,my
        if (yy(jj-1)>yy(jj)) return
    end do
    do kk=2,mz
        if (zz(kk-1)>zz(kk)) return
    end do
    do ll=2,mt
        if (tt(ll-1)>tt(ll)) return
    end do

    ! arguments o.k.

    ier = 0

    ! set pointers for integer work space iw

    jy = mx+1
    kz = mx+my+1
    lt = mx+my+mz+1

    if (intpol(4)==1) then

        ! linearly interpolate in nt, set work space pointers and scales

        l2 = 1
        l3 = l2
        l4 = l3+mt
        l5 = l4
        l6 = l5
        l7 = l6
        l8 = l7+mxmymz
        l9 = l8+mxmymz
        call linmx(nt,t,mt,tt,iw(lt),w(l3))
        k2 = l9

        if (intpol(3)==1) then
            ! linear in z
            k3 = k2
            k4 = k3+mz
            k5 = k4
            k6 = k5
            k7 = k6
            k8 = k7+mxmy
            k9 = k8+mxmy
            call linmx(nz,z,mz,zz,iw(kz),w(k3))
            j2 = k9
        else
            ! cubic in z
            k3 = k2+mz
            k4 = k3+mz
            k5 = k4+mz
            k6 = k5+mz
            k7 = k6+mxmy
            k8 = k7+mxmy
            k9 = k8+mxmy
            call cubnmx(nz,z,mz,zz,iw(kz),w(k2),w(k3),w(k4),w(k5))
            j2 = k9+mxmy
        end if

        if (intpol(2) == 1) then
            ! linear in y
            j3 = j2
            j4 = j3+my
            j5 = j4
            j6 = j5
            j7 = j6
            j8 = j7+mx
            j9 = j8+mx
            call linmx(ny,y,my,yy,iw(jy),w(j3))
            i2 = j9
        else
            ! cubic in y
            j3 = j2+my
            j4 = j3+my
            j5 = j4+my
            j6 = j5+my
            j7 = j6+mx
            j8 = j7+mx
            j9 = j8+mx
            call cubnmx(ny,y,my,yy,iw(jy),w(j2),w(j3),w(j4),w(j5))
            i2 = j9+mx
        end if

        if (intpol(1) == 1) then
            ! linear in x
            i3 = i2
            i4 = i3
            i5 = i4
            call linmx(nx,x,mx,xx,iw,w(i3))
        else
            ! cubic in x
            i3 = i2+mx
            i4 = i3+mx
            i5 = i4+mx
            call cubnmx(nx,x,mx,xx,iw,w(i2),w(i3),w(i4),w(i5))
        end if

        ! linearly interpolate in t

        call lint4(nx,ny,nz,nt,p,mx,my,mz,mt,mxmy,mxmymz,q,intpol,&
                    iw(lt),w(l3),w(l7),w(l8),&
                    iw(kz),w(k2),w(k3),w(k4),w(k5),w(k6),w(k7),w(k8),w(k9),&
                    iw(jy),w(j2),w(j3),w(j4),w(j5),w(j6),w(j7),w(j8),w(j9),&
                    iw,w(i2),w(i3),w(i4),w(i5))

    else

        ! cubically interpolate in t

        l2 = 1
        l3 = l2+mt
        l4 = l3+mt
        l5 = l4+mt
        l6 = l5+mt
        l7 = l6+mxmymz
        l8 = l7+mxmymz
        l9 = l8+mxmymz
        call cubnmx(nt,t,mt,tt,iw(lt),w(l2),w(l3),w(l4),w(l5))
        k2 = l9+mxmymz

        if (intpol(3)==1) then
            ! linear in z
            k3 = k2
            k4 = k3+mz
            k5 = k4
            k6 = k5
            k7 = k6
            k8 = k7+mxmy
            k9 = k8+mxmy
            call linmx(nz,z,mz,zz,iw(kz),w(k3))
            j2 = k9
        else
            ! cubic in z
            k3 = k2+mz
            k4 = k3+mz
            k5 = k4+mz
            k6 = k5+mz
            k7 = k6+mxmy
            k8 = k7+mxmy
            k9 = k8+mxmy
            call cubnmx(nz,z,mz,zz,iw(kz),w(k2),w(k3),w(k4),w(k5))
            j2 = k9+mxmy
        end if

        if (intpol(2) == 1) then
            j3 = j2
            j4 = j3+my
            j5 = j4
            j6 = j5
            j7 = j6
            j8 = j7+mx
            j9 = j8+mx
            call linmx(ny,y,my,yy,iw(jy),w(j3))
            i2 = j9
        else
            j3 = j2+my
            j4 = j3+my
            j5 = j4+my
            j6 = j5+my
            j7 = j6+mx
            j8 = j7+mx
            j9 = j8+mx
            call cubnmx(ny,y,my,yy,iw(jy),w(j2),w(j3),w(j4),w(j5))
            i2 = j9+mx
        end if

        ! set work space portion and indices which depend on x interpolation

        if (intpol(1) == 1) then
            i3 = i2
            i4 = i3
            i5 = i4
            call linmx(nx,x,mx,xx,iw,w(i3))
        else
            i3 = i2+mx
            i4 = i3+mx
            i5 = i4+mx
            call cubnmx(nx,x,mx,xx,iw,w(i2),w(i3),w(i4),w(i5))
        end if

        ! cubically interpolate in t

        call cubt4(nx,ny,nz,nt,p,mx,my,mz,mt,mxmy,mxmymz,q,intpol,&
                    iw(lt),w(l2),w(l3),w(l4),w(l5),w(l6),w(l7),w(l8),w(l9),&
                    iw(kz),w(k2),w(k3),w(k4),w(k5),w(k6),w(k7),w(k8),w(k9),&
                    iw(jy),w(j2),w(j3),w(j4),w(j5),w(j6),w(j7),w(j8),w(j9),&
                    iw,w(i2),w(i3),w(i4),w(i5))

    end if

    end subroutine rgrd4