rgrd3 Subroutine

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

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

method

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

requirements

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

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

must lie within the (X,Y,Z) region

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

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

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.
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) :: p(nx,ny,nz)

a real(wp) nx by ny by nz array of values given on the (x,y,z) 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.
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(out) :: q(mx,my,mz)

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

an integer vector of dimension 3 which sets linear or cubic interpolation in each of the x,y,z 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.

values other than 1 or 3 in intpol are not allowed (ier = 5).
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 rgrd3
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*(mx+my) if intpol(2) = 3
  • lwz = 2mxmy+mz if intpol(3) = 1
  • lwz = 4(mxmy+mz) if intpol(3) = 3

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

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

the integer length of the integer work space iw. liw must be at least mx+my+mz
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) < 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.
  • 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) 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 rgrd3)

  • ier = 4 if one of the grids x,y,z is not strictly monotonically increasing or if one of the grids xx,yy,zz 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)
  • 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
  • ier = 5 if lw or liw is too small (insufficient work space)
  • ier = 6 if any of intpol(1),intpol(2),intpol(3) is not equal to 1 or 3

Calls

proc~~rgrd3~~CallsGraph proc~rgrd3 regridpack_module::rgrd3 proc~cubnmx regridpack_module::cubnmx proc~rgrd3->proc~cubnmx proc~cubt3 regridpack_module::cubt3 proc~rgrd3->proc~cubt3 proc~linmx regridpack_module::linmx proc~rgrd3->proc~linmx proc~lint3 regridpack_module::lint3 proc~rgrd3->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~~rgrd3~~CalledByGraph proc~rgrd3 regridpack_module::rgrd3 interface~regrid regridpack_module::regrid interface~regrid->proc~rgrd3 proc~rgrd3_wrapper regridpack_module::rgrd3_wrapper interface~regrid->proc~rgrd3_wrapper proc~rgrd3_wrapper->proc~rgrd3

Source Code

    subroutine rgrd3(nx,ny,nz,x,y,z,p,mx,my,mz,xx,yy,zz,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)      :: 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)      :: 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*(mx+my)       if intpol(2) = 3
                                    !! * lwz = 2*mx*my+mz      if intpol(3) = 1
                                    !! * lwz = 4*(mx*my+mz)    if intpol(3) = 3
                                    !!
                                    !! then lw must be greater than or equal to lwx+lwy+lwz
    integer,intent(in)      :: liw  !! the integer length of the integer work space iw.  liw must be at least mx+my+mz
    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) < 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.
                                    !! * 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)
                                    !!    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 rgrd3)
                                    !! * ier = 4 if one of the grids x,y,z is not strictly monotonically
                                    !!   increasing or if one of the grids xx,yy,zz 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)
                                    !!    * 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
                                    !! * ier = 5 if lw or liw is too small (insufficient work space)
                                    !! * ier = 6 if any of intpol(1),intpol(2),intpol(3) is not equal to 1 or 3
    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)     :: p(nx,ny,nz)  !! a real(wp) nx by ny by nz array of values given on the (x,y,z) grid
    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(out)    :: q(mx,my,mz)  !! a real(wp) mx by my by mz array of values
                                            !! on the (xx,yy,zz) grid which are
                                            !! interpolated from p on the (x,y,z) 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 rgrd3
    integer,intent(in)      :: intpol(3)    !! an integer vector of dimension 3 which sets linear or cubic
                                            !! interpolation in each of the x,y,z 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.
                                            !!
                                            !! values other than 1 or 3 in intpol are not allowed (ier = 5).
    integer,intent(inout)   :: iw(liw)  !! an integer work space of length at least
                                        !! liw which must be provided in the
                                        !! routine calling rgrd3

    integer  :: i,ii,j,jj,k,kk
    integer  :: i2,i3,i4,i5
    integer  :: j2,j3,j4,j5,j6,j7,j8,j9
    integer  :: k2,k3,k4,k5,k6,k7,k8,k9
    integer  :: lwx,lwy,lwz,jy,kz,mxmy

    ! check input arguments

    ! check (xx,yy,zz) grid resolution
    ier = 1
    if (min(mx,my,mz) < 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

    ! check (x,y,z) 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

    ! check work space length input and set minimum
    ier = 5
    mxmy = mx*my
    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*(my+mx)
    end if
    if (intpol(3)==1) then
        lwz = (2*mxmy+mz)
    else
        lwz = 4*(mxmy+mz)
    end if
    if (lw < lwx+lwy+lwz) return
    if (liw < mx+my+mz) return

    ! check (xx,yy,zz) grid contained in (x,y,z) 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

    ! 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 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

    ! arguments o.k.

    ier = 0
    jy = mx+1
    kz = mx+my+1
    if (intpol(3)==1) then

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

        k2 = 1
        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

        ! set indices and scales which depend on y interpolation

        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

        ! set indices and scales which depend on x interpolation

        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
        call lint3(nx,ny,nz,p,mx,my,mxmy,mz,q,intpol,iw(kz),&
                    w(k3),w(k7),w(k8),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 z

        k2 = 1
        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

        ! set indices which depend on y interpolation

        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
        call cubt3(nx,ny,nz,p,mx,my,mxmy,mz,q,intpol,&
                    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 rgrd3