rgrd2 Subroutine

private subroutine rgrd2(nx, ny, x, y, p, mx, my, xx, yy, q, intpol, w, lw, iw, liw, ier)

subroutine rgrd2 interpolates the values p(i,j) on the orthogonal grid (x(i),y(j)) for i=1,...,nx and j=1,...,ny onto q(ii,jj) on the orthogonal grid (xx(ii),yy(jj)) for ii=1,...,mx and jj=1,...,my.

method

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

requirements

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

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

must lie within the (X,Y) region

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

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

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

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

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

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

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

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

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

the integer length of the integer work space iw. liw must be at least mx+my
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) < 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
  • ier = 3 if xx(1) < x(1) or x(nx) < xx(mx) (or) yy(1) < y(1) or y(ny) < yy(my) (or) 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 rgrd2)

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

Calls

proc~~rgrd2~~CallsGraph proc~rgrd2 regridpack_module::rgrd2 proc~cubnmx regridpack_module::cubnmx proc~rgrd2->proc~cubnmx proc~cubt2 regridpack_module::cubt2 proc~rgrd2->proc~cubt2 proc~linmx regridpack_module::linmx proc~rgrd2->proc~linmx proc~lint2 regridpack_module::lint2 proc~rgrd2->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~~rgrd2~~CalledByGraph proc~rgrd2 regridpack_module::rgrd2 interface~regrid regridpack_module::regrid interface~regrid->proc~rgrd2 proc~rgrd2_wrapper regridpack_module::rgrd2_wrapper interface~regrid->proc~rgrd2_wrapper proc~rgrd2_wrapper->proc~rgrd2

Source Code

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

    integer  :: i,ii,j,jj,j2,j3,j4,j5,j6,j7,j8,j9,i2,i3,i4,i5
    integer  :: jy,lwx,lwy

    ! check input arguments

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

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

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

    ! check work space lengths
    ier = 5
    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 (lw < lwx+lwy) return
    if (liw < mx+my) return

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

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

    ! arguments o.k.

    ier = 0

    ! set pointer in integer work space

    jy = mx+1
    if (intpol(2) ==1) then

        ! linearly interpolate in y

        j2 = 1
        j3 = j2
        j4 = j3+my
        j5 = j4
        j6 = j5
        j7 = j6
        j8 = j7+mx
        j9 = j8+mx

        ! set y interpolation indices and scales and linearly interpolate

        call linmx(ny,y,my,yy,iw(jy),w(j3))
        i2 = j9

        ! 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 lint2(nx,ny,p,mx,my,q,intpol,iw(jy),w(j3),w(j7),w(j8),iw,w(i2),w(i3),w(i4),w(i5))

    else

        ! cubically interpolate in y, set indice pointers

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

        ! 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 cubt2(nx,ny,p,mx,my,q,intpol,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 rgrd2