rgrd3u Subroutine

private subroutine rgrd3u(nx, ny, nz, p, mx, my, mz, q, intpol, w, lw, iw, liw, ier)

subroutine rgrd3u interpolates the nx by ny by nz array p onto the mx by my by mz array q. it is assumed that p and q are values on uniform nx by ny by nz and mx by my by mz grids which are superimposed on the same box region (INCLUDING BOUNDARIES). if p and q are values on nonuniform orthogonal grids and/or if the grid on which q is defined lies within the p grid then subroutine rgrd3 should be used.

method

linear or cubic interpolation (see intpol) is used in each direction for which the q grid is not a subgrid of the p grid. [the mx (my,mz) uniform grid is a subgrid of the nx (ny,nz) uniform grid if and only if mx-1 (my-1,nz-1) divides nx-1 (ny-1,nz-1)]. Values are set directly without (the need for) interpolation in subgrid directions.

Arguments

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

the integer first dimension of p. nx > 1 if intpol(1) = 1 or nx > 3 if intpol(1) = 3 is required (see ier = 2).
integer, intent(in) :: ny

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

the integer third dimension of p. nz > 1 if intpol(3) = 1 or nz > 3 if intpol(3) = 3 is required (see ier = 2)
real(kind=wp), intent(in) :: p(nx,ny,nz)

a real(wp) nx by ny by nz array of given values
integer, intent(in) :: mx

the integer first dimension of q. mx > 1 is required (see ier = 1)
integer, intent(in) :: my

the integer second dimension of q. my > 1 is required (see ier = 1)
integer, intent(in) :: mz

the integer third dimension of q. mz > 1 is required (see ier = 1)
real(kind=wp), intent(out) :: q(mx,my,mz)

a real(wp) mx by my by mz array of values which are interpolated from p.
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 = 3).
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 rgrd3u
integer, intent(in) :: lw

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

  • let lwx = 1 if mx-1 divides nx-1; otherwise
  • let lwx = mx if intpol(1) = 1 or
  • let lwx = 4*mx if intpol(1) = 3
  • let lwy = 0 if my-1 divides ny-1; otherwise
  • let lwy = my+2*mx if intpol(2) = 1 or
  • let lwy = 4*(mx+my) if intpol(2) = 3
  • let lwz = 0 if mz-1 divides nz-1; otherwise
  • let lwz = 2mxmy+mz if intpol(3) = 1 or
  • let 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 rgrd3u
integer, intent(in) :: liw

the integer length of the integer work space iw. liw must be greater than or equal to 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) < 2
  • 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 any of intpol(1),intpol(2),intpol(3) is not equal to 1 or 3
  • ier = 4 if lw or liw is too small (insufficient work space)

Calls

proc~~rgrd3u~~CallsGraph proc~rgrd3u regridpack_module::rgrd3u proc~cubnmxu regridpack_module::cubnmxu proc~rgrd3u->proc~cubnmxu proc~cubt3u regridpack_module::cubt3u proc~rgrd3u->proc~cubt3u proc~linmxu regridpack_module::linmxu proc~rgrd3u->proc~linmxu proc~lint3u regridpack_module::lint3u proc~rgrd3u->proc~lint3u proc~cubt2u regridpack_module::cubt2u proc~cubt3u->proc~cubt2u proc~lint2u regridpack_module::lint2u proc~cubt3u->proc~lint2u proc~lint3u->proc~cubt2u proc~lint3u->proc~lint2u proc~cubt1u regridpack_module::cubt1u proc~cubt2u->proc~cubt1u proc~lint1u regridpack_module::lint1u proc~cubt2u->proc~lint1u proc~lint2u->proc~cubt1u proc~lint2u->proc~lint1u

Called by

proc~~rgrd3u~~CalledByGraph proc~rgrd3u regridpack_module::rgrd3u interface~regrid regridpack_module::regrid interface~regrid->proc~rgrd3u

Source Code

    subroutine rgrd3u(nx,ny,nz,p,mx,my,mz,q,intpol,w,lw,iw,liw,ier)

    implicit none

    integer,intent(in)     :: nx    !! the integer first dimension of p.  nx > 1 if intpol(1) = 1 or
                                    !! nx > 3 if intpol(1) = 3 is required (see ier = 2).
    integer,intent(in)     :: ny    !! the integer second dimension of p.  ny > 1 if intpol(2) = 1 or
                                    !! ny > 3 if intpol(2) = 3 is required (see ier = 2).
    integer,intent(in)     :: nz    !! the integer third dimension of p.  nz > 1 if intpol(3) = 1 or
                                    !! nz > 3 if intpol(3) = 3 is required (see ier = 2)
    integer,intent(in)     :: mx    !! the integer first dimension of q.  mx > 1 is required (see ier = 1)
    integer,intent(in)     :: my    !! the integer second dimension of q. my > 1 is required (see ier = 1)
    integer,intent(in)     :: mz    !! the integer third dimension of q. mz > 1 is required (see ier = 1)
    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 = 3).
    integer,intent(in)     :: lw        !! the integer length of the real(wp) work space w.
                                        !!
                                        !! * let lwx = 1 if mx-1 divides nx-1; otherwise
                                        !! * let lwx = mx if intpol(1) = 1 or
                                        !! * let lwx = 4*mx if intpol(1) = 3
                                        !! * let lwy = 0 if my-1 divides ny-1; otherwise
                                        !! * let lwy = my+2*mx if intpol(2) = 1 or
                                        !! * let lwy = 4*(mx+my) if intpol(2) = 3
                                        !! * let lwz = 0 if mz-1 divides nz-1; otherwise
                                        !! * let lwz = 2*mx*my+mz if intpol(3) = 1 or
                                        !! * let 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 greater than or equal to mx+my+mz
    integer,intent(inout)  :: iw(liw)   !! an integer work space of length at least liw
                                        !! which must be provided in the
                                        !! routine calling rgrd3u
    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) < 2
                                        !! * 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 any of intpol(1),intpol(2),intpol(3) is not equal to 1 or 3
                                        !! * ier = 4 if lw or liw is too small (insufficient work space)
    real(wp),intent(in)    :: p(nx,ny,nz)   !! a real(wp) nx by ny by nz array of given values
    real(wp),intent(out)   :: q(mx,my,mz)   !! a real(wp) mx by my by mz array of values which are interpolated from p.
    real(wp),intent(inout) :: w(lw) !! a real(wp) work space of length at least lw
                                    !! which must be provided in the
                                    !! routine calling rgrd3u

    integer  :: inmx,jnmy,knmz,isubx,jsuby,ksubz
    integer  :: lwx,lwy,lwz,mxmy,jy,kz
    integer  :: i2,i3,i4,i5
    integer  :: j2,j3,j4,j5,j6,j7,j8,j9
    integer  :: k2,k3,k4,k5,k6,k7,k8,k9

    ! check input arguments

    ! check mx,my,mz
    ier = 1
    if (min(mx,my,mz) < 1) return

    ! check intpol
    ier = 3
    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 nx,ny,nz
    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

    ! set subgrid indicators
    inmx = (nx-1)/(mx-1)
    jnmy = (ny-1)/(my-1)
    knmz = (nz-1)/(mz-1)
    isubx = nx - inmx*(mx-1)
    jsuby = ny - jnmy*(my-1)
    ksubz = nz - knmz*(mz-1)

    ! check work space lengths
    ier = 4
    mxmy = mx*my
    lwx = 1
    if (isubx/=1) then
        if (intpol(1)==1) then
            lwx = mx
        else
            lwx = 4*mx
        end if
    end if
    lwy = 0
    if (jsuby/=1) then
        if (intpol(2)==1) then
            lwy = (2*mx+my)
        else
            lwy = 4*(mx+my)
        end if
    end if
    lwz = 0
    if (ksubz/=1) then
        if (intpol(3)==1) then
            lwz = (2*mxmy+mz)
        else
            lwz = 4*(mxmy+mz)
        end if
    end if
    if (lw < lwx+lwy+lwz) return
    if (liw < mx+my+mz) return

    ! arguments o.k.

    ier = 0
    jy = mx+1
    kz = mx+my+1

    ! preset work space pointers

    k2 = 1
    k3 = 1
    k4 = 1
    k5 = 1
    k6 = 1
    k7 = 1
    k8 = 1
    k9 = 1
    j2 = 1
    j3 = 1
    j4 = 1
    j5 = 1
    j6 = 1
    j7 = 1
    j8 = 1
    j9 = 1
    i2 = 1
    i3 = 1
    i4 = 1
    i5 = 1

    if (intpol(3)==1) then
        if (ksubz/=1) then
            ! linearly interpolate in nz, set work space pointers
            k2 = 1
            k3 = k2
            k4 = k3+mz
            k5 = k4
            k6 = k5
            k7 = k6
            k8 = k7+mxmy
            k9 = k8+mxmy
            ! set z interpolation indices and scales
            call linmxu(nz,mz,iw(kz),w(k3))
            j2 = k9
            i2 = k9
        end if

        if (jsuby/=1) then
            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 linmxu(ny,my,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 cubnmxu(ny,my,iw(jy),w(j2),w(j3),w(j4),w(j5))
                i2 = j9+mx
            end if
        end if

        if (isubx/=1) then
            if (intpol(1) == 1) then
                ! linear in x
                i3 = i2
                i4 = i3
                i5 = i4
                call linmxu(nx,mx,iw,w(i3))
            else
                ! cubic in x
                i3 = i2+mx
                i4 = i3+mx
                i5 = i4+mx
                call cubnmxu(nx,mx,iw,w(i2),w(i3),w(i4),w(i5))
            end if
        end if

        ! linearly interpolate p onto q in z

        call lint3u(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),&
                    inmx,jnmy,knmz,isubx,jsuby,ksubz)

    else

        ! cubically interpolate in z

        if (ksubz/=1) then
            k2 = 1
            k3 = k2+mz
            k4 = k3+mz
            k5 = k4+mz
            k6 = k5+mz
            k7 = k6+mxmy
            k8 = k7+mxmy
            k9 = k8+mxmy
            call cubnmxu(nz,mz,iw(kz),w(k2),w(k3),w(k4),w(k5))
            j2 = k9+mxmy
            i2 = j2
        end if

        if (jsuby/=1) then
            if (intpol(2) == 1) then
                j3 = j2
                j4 = j3+my
                j5 = j4
                j6 = j5
                j7 = j6
                j8 = j7+mx
                j9 = j8+mx
                call linmxu(ny,my,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 cubnmxu(ny,my,iw(jy),w(j2),w(j3),w(j4),w(j5))
                i2 = j9+mx
            end if
        end if

        if (isubx/=1) then
            if (intpol(1) == 1) then
                i3 = i2
                i4 = i3
                i5 = i4
                call linmxu(nx,mx,iw,w(i3))
            else
                i3 = i2+mx
                i4 = i3+mx
                i5 = i4+mx
                call cubnmxu(nx,mx,iw,w(i2),w(i3),w(i4),w(i5))
            end if
        end if

        ! cubically interpolate p onto q in z

        call cubt3u(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),&
                    inmx,jnmy,knmz,isubx,jsuby,ksubz)

    end if

    end subroutine rgrd3u