rgrd4u Subroutine

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

subroutine rgrd4u interpolates the nx by ny by nz by nt array p onto the mx by my by mz by mt array q. it is assumed that p and q are values on uniform nx by ny by nz by nt and mx by my by mz by mt 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 rgrd4 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,mt) uniform grid is a subgrid of the nx (ny,nz,nt) uniform grid if and only if mx-1 (my-1,nz-1,nt-1) divides nx-1 (ny-1,nz-1,nt-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)
integer, intent(in) :: nt

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

a real(wp) nx by ny by nz by nt 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)
integer, intent(in) :: mt

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

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

the integer length of the 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+2mx 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
  • let lwt = 0 if mt-1 divides nt-1; otherwise let lwt = 2mxmymz+mt if intpol(4) = 1 or let lwt = 4(mxmymz+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 rgrd4u
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) < 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 (or) nt < 2 when intpol(4)=1 or nt < 4 when intpol(4)=3.
  • ier = 3 if any of intpol(1),intpol(2),intpol(3),intpol(4) is not equal to 1 or 3.
  • ier = 4 if lw or liw is too small (insufficient work space)

Calls

proc~~rgrd4u~~CallsGraph proc~rgrd4u regridpack_module::rgrd4u proc~cubnmxu regridpack_module::cubnmxu proc~rgrd4u->proc~cubnmxu proc~cubt4u regridpack_module::cubt4u proc~rgrd4u->proc~cubt4u proc~linmxu regridpack_module::linmxu proc~rgrd4u->proc~linmxu proc~lint4u regridpack_module::lint4u proc~rgrd4u->proc~lint4u proc~cubt3u regridpack_module::cubt3u proc~cubt4u->proc~cubt3u proc~lint3u regridpack_module::lint3u proc~cubt4u->proc~lint3u proc~lint4u->proc~cubt3u proc~lint4u->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~~rgrd4u~~CalledByGraph proc~rgrd4u regridpack_module::rgrd4u interface~regrid regridpack_module::regrid interface~regrid->proc~rgrd4u

Source Code

    subroutine rgrd4u(nx,ny,nz,nt,p,mx,my,mz,mt,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)      :: nt   !! the integer fourth dimension of p.  nt > 1 if intpol(4) = 1 or
                                    !! nt > 3 if intpol(4) = 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)      :: mt   !! the integer fourth dimension of q. mt > 1 is required (see ier = 1)
    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 = 3).
    integer,intent(in)      :: liw          !! the integer length of the integer work space iw.
                                            !! liw must be at least mx+my+mz+mt
    integer,intent(inout)   :: iw(liw)      !! an integer work space of length at least liw
                                            !! which must be provided in the routine calling rgrd4u
    integer,intent(in)      :: lw           !! the integer length of the 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
                                            !! * let lwt = 0 if mt-1 divides nt-1; otherwise
                                            !!   let lwt = 2*mx*my*mz+mt if intpol(4) = 1 or
                                            !!   let 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(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) < 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 (or)
                                            !!                     nt < 2 when intpol(4)=1 or nt < 4 when intpol(4)=3.
                                            !! * ier = 3 if any of intpol(1),intpol(2),intpol(3),intpol(4)  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,nt)   !! a real(wp) nx by ny by nz by nt array of given values
    real(wp),intent(out)    :: q(mx,my,mz,mt)   !! a real(wp) mx by my by mz by mt 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 rgrd4u


    integer :: inmx,jnmy,knmz,lnmt,isubx,jsuby,ksubz,lsubt
    integer :: mxmy,mxmymz,lwx,lwy,lwz,lwt,jy,kz,lt
    integer :: i2,i3,i4,i5
    integer :: j2,j3,j4,j5,j6,j7,j8,j9
    integer :: k2,k3,k4,k5,k6,k7,k8,k9
    integer :: l2,l3,l4,l5,l6,l7,l8,l9

    ! check input arguments

    ! check mx,my,mz,mt
    ier = 1
    if (min(mx,my,mz,mt) < 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
    if (intpol(4)/=1 .and. intpol(4)/=3) return

    ! check nx,ny,nz,nt
    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

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

    ! check work space length input
    ier = 4
    mxmy = mx*my
    mxmymz = mxmy*mz
    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*my+4*mx
        end if
    end if
    lwz = 0
    if (ksubz/=1) then
        if (intpol(3)==1) then
            lwz = (2*mxmy+mz)
        else
            lwz = 4*mxmy+4*mz
        end if
    end if
    lwt = 0
    if (lsubt/=1) then
        if (intpol(4)==1) then
            lwt = (2*mxmymz+mt)
        else
            lwt = 4*mxmymz+4*mt
        end if
    end if

    if (lw < lwx+lwy+lwz+lwt) return
    if (liw < mx+my+mz+mt) return

    ! arguments o.k.

    ier = 0
    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 linmxu(nt,mt,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 linmxu(nz,mz,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 cubnmxu(nz,mz,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 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

        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

        ! linearly interpolate in t

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

    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 cubnmxu(nt,mt,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 linmxu(nz,mz,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 cubnmxu(nz,mz,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 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

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

        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

        ! cubically interpolate in t

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

    end if

    end subroutine rgrd4u