BILIN() and BMATVEC()

Fortran version:

For I/O API Version 3.2: BILIN() and BMATVEC() are Fortran-90 generic routines, with INTERFACEs defined in MODULE M3UTILIO. The prior versions.

For I/O API Versions 3.1 and before, BILIN() and BMATVEC() look like BILIN11L() and BMATVEC11() below, and rely upon single-indexing and F77/"C-void pointer" argument passing to achieve the required polymorphism.

    INTERFACE BILIN( IX, AX, ... )

        SUBROUTINE  BILIN11L( M, N, P, IX, AX, V, Y )
        INTEGER, INTENT(IN   ) :: M               ! length of input  vector
        INTEGER, INTENT(IN   ) :: N               ! length of output vector
        INTEGER, INTENT(IN   ) :: P               ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,N )       ! index array
        REAL   , INTENT(IN   ) :: AX( 4,N )       ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( M,P )        ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( N,P )        ! P-layered output vector
        END SUBROUTINE  BILIN11L

        SUBROUTINE  BILIN12L( M, NC, NR, P, IX, AX, V, Y )
        INTEGER, INTENT(IN   ) :: M               ! length of input  vector
        INTEGER, INTENT(IN   ) :: NC, NR          ! dims of output grid
        INTEGER, INTENT(IN   ) :: P               ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,NC*NR )   ! index array
        REAL   , INTENT(IN   ) :: AX( 4,NC*NR )   ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( M,P )        ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( NC,NR,P )    ! P-layered output vector
        END SUBROUTINE  BILIN12L

        SUBROUTINE  BILIN21L( MC, MR, N, P, IX, AX, V, Y )
        INTEGER, INTENT(IN   ) :: MC, MR          ! length of input  vector
        INTEGER, INTENT(IN   ) :: N               ! dims of output grid
        INTEGER, INTENT(IN   ) :: P               ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,N )       ! index array
        REAL   , INTENT(IN   ) :: AX( 4,N )       ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( MC,MR,P )    ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( N,P )        ! P-layered output vector
        END SUBROUTINE  BILIN21L

        SUBROUTINE  BILIN22L( MC, MR, NC, NR, P, IX, AX, V, Y )
        INTEGER, INTENT(IN   ) :: MC, MR          ! length of input  grid
        INTEGER, INTENT(IN   ) :: NC, NR          ! dims of output grid
        INTEGER, INTENT(IN   ) :: P               ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,NC*NR )   ! index array
        REAL   , INTENT(IN   ) :: AX( 4,NC*NR )   ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( MC,MR,P )    ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( NC,NR,P )    ! P-layered output vector
        END SUBROUTINE  BILIN22L

        SUBROUTINE  BILIN11( M, N, IX, AX, V, Y )
        INTEGER, INTENT(IN   ) :: M               ! length of input  vector
        INTEGER, INTENT(IN   ) :: N               ! length of output vector
        INTEGER, INTENT(IN   ) :: IX( 4,N )       ! index array
        REAL   , INTENT(IN   ) :: AX( 4,N )       ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( M )          ! nonlayered input  vector
        REAL   , INTENT(  OUT) :: Y( N )          ! nonlayered output vector
        END SUBROUTINE  BILIN11

        SUBROUTINE  BILIN12( M, NC, NR, IX, AX, V, Y )
        INTEGER, INTENT(IN   ) :: M               ! length of input  vector
        INTEGER, INTENT(IN   ) :: NC, NR          ! dims of output grid
        INTEGER, INTENT(IN   ) :: IX( 4,NC*NR )   ! index array
        REAL   , INTENT(IN   ) :: AX( 4,NC*NR )   ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( M )          ! nonlayered input  vector
        REAL   , INTENT(  OUT) :: Y( NC,NR )      ! nonlayered output vector
        END SUBROUTINE  BILIN12

        SUBROUTINE  BILIN21( MC, MR, N, IX, AX, V, Y )
        INTEGER, INTENT(IN   ) :: MC, MR          ! length of input  vector
        INTEGER, INTENT(IN   ) :: N               ! dims of output grid
        INTEGER, INTENT(IN   ) :: IX( 4,N )       ! index array
        REAL   , INTENT(IN   ) :: AX( 4,N )       ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( MC,MR )      ! nonlayered input  vector
        REAL   , INTENT(  OUT) :: Y( N )          ! nonlayered output vector
        END SUBROUTINE  BILIN21

        SUBROUTINE  BILIN22( MC, MR, NC, NR, IX, AX, V, Y )
        INTEGER, INTENT(IN   ) :: MC, MR          ! length of input  grid
        INTEGER, INTENT(IN   ) :: NC, NR          ! dims of output grid
        INTEGER, INTENT(IN   ) :: IX( 4,NC*NR )   ! index array
        REAL   , INTENT(IN   ) :: AX( 4,NC*NR )   ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( MC,MR )      ! nonlayered input  vector
        REAL   , INTENT(  OUT) :: Y( NC,NR )      ! nonlayered output vector
        END SUBROUTINE  BILIN22

    INTERFACE BMATVEC( ..., IX, AX, ... )

        !!  non-layered cases:  you really should be using BILIN(), but...

        SUBROUTINE  BMATVEC01( M, N, P, IX, AX, V, Y )
        INTEGER, INTENT(IN   ) :: M               ! length of input  vector
        INTEGER, INTENT(IN   ) :: N               ! length of output vector
        INTEGER, INTENT(IN   ) :: P               ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,N )       ! index array
        REAL   , INTENT(IN   ) :: AX( 4,N )       ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( M )          ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( N )          ! P-layered output vector
        END SUBROUTINE  BMATVEC01

        SUBROUTINE  BMATVEC02( M, NC, NR, P, IX, AX, V, Y )
        INTEGER, INTENT(IN   ) :: M               ! length of input  vector
        INTEGER, INTENT(IN   ) :: NC, NR          ! length of output vector
        INTEGER, INTENT(IN   ) :: P               ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,NC*NR )   ! index array
        REAL   , INTENT(IN   ) :: AX( 4,NC*NR )   ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( M )          ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( NC,NR )      ! P-layered output vector
        END SUBROUTINE  BMATVEC02

        SUBROUTINE  BMATVEC021( MC, MR, N, P, IX, AX, V, Y )
        INTEGER, INTENT(IN   ) :: MC, MR          ! length of input  vector
        INTEGER, INTENT(IN   ) :: N               ! length of output vector
        INTEGER, INTENT(IN   ) :: P               ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,N )       ! index array
        REAL   , INTENT(IN   ) :: AX( 4,N )       ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( MC,MR )      ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( N )          ! P-layered output vector
        END SUBROUTINE  BMATVEC021

        SUBROUTINE  BMATVEC022( MC,MR, NC,NR, P, IX, AX, V, Y )
        INTEGER, INTENT(IN   ) :: MC, MR          ! length of input  vector
        INTEGER, INTENT(IN   ) :: NC, NR          ! length of output vector
        INTEGER, INTENT(IN   ) :: P               ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,NC*NR )   ! index array
        REAL   , INTENT(IN   ) :: AX( 4,NC*NR )   ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( MC,MR )      ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( NC,NR )      ! P-layered output vector
        END SUBROUTINE  BMATVEC022

        !!  layered cases:

        SUBROUTINE  BMATVEC11( M, N, P, IX, AX, V, Y )
        INTEGER, INTENT(IN   ) :: M               ! length of input  vector
        INTEGER, INTENT(IN   ) :: N               ! length of output vector
        INTEGER, INTENT(IN   ) :: P               ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,N )       ! index array
        REAL   , INTENT(IN   ) :: AX( 4,N )       ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( M,P )        ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( P,N )        ! P-layered output vector
        END SUBROUTINE  BMATVEC11

        SUBROUTINE  BMATVEC12( M, NC, NR, P, IX, AX, V, Y )
        INTEGER, INTENT(IN   ) :: M               ! length of input  vector
        INTEGER, INTENT(IN   ) :: NC, NR          ! length of output vector
        INTEGER, INTENT(IN   ) :: P               ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,NC*NR )   ! index array
        REAL   , INTENT(IN   ) :: AX( 4,NC*NR )   ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( M,P )        ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( P,NC,NR )    ! P-layered output vector
        END SUBROUTINE  BMATVEC12

        SUBROUTINE  BMATVEC21( MC, MR, N, P, IX, AX, V, Y )
        INTEGER, INTENT(IN   ) :: MC, MR          ! length of input  vector
        INTEGER, INTENT(IN   ) :: N               ! length of output vector
        INTEGER, INTENT(IN   ) :: P               ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,N )       ! index array
        REAL   , INTENT(IN   ) :: AX( 4,N )       ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( MC,MR,P )    ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( P,N )        ! P-layered output vector
        END SUBROUTINE  BMATVEC21

        SUBROUTINE  BMATVEC22( MC, MR, NC, NR, P, IX, AX, V, Y )
        INTEGER, INTENT(IN   ) :: MC, MR          ! length of input  vector
        INTEGER, INTENT(IN   ) :: NC, NR          ! length of output vector
        INTEGER, INTENT(IN   ) :: P               ! number of layers
        INTEGER, INTENT(IN   ) :: IX( 4,NC*NR )   ! index array
        REAL   , INTENT(IN   ) :: AX( 4,NC*NR )   ! 4-band coeff matrix
        REAL   , INTENT(IN   ) :: V( MC,MR,P )    ! P-layered input  vector
        REAL   , INTENT(  OUT) :: Y( P,NC,NR )    ! P-layered output vector
        END SUBROUTINE  BMATVEC22

C version: none

Summary:

For Fortran-90 declarations and interface checking:
    USE M3UTILIO
    

BILIN() and BMATVEC() apply 4-band sparse matrix <IX,AX>; (e.g., as computed by subroutine UNGRIDB()) to array V, to generate output array C, according to the following equation for output vector value at row R and layer L"

V( R,L ) = SUMJ = 1...4 [ AX(J,R) U( IX(J,R),L ) ]
Note that BMATVEC simultaneously does a transpose on the output subscripting, as is used by the SMOKE layer-fractions processor: the input met variables V are subscripted by (COLUMN,ROW,LAYER), whereas the output interpolated variables C are subscripted by (LAYER,STACKNUMBER).

For single-indexed bilinear interpolation of gridded data V having dimensions V(NC,NR) (as specified in the Fortran storage order) to a list of locations having grid-normal coordinates <X(S),Y(S)>, let

    C(S) = 1 + INT( X(S) )
    R(S) = 1 + INT( Y(S) )
    P(S) = AMOD( X(S), 1.0 )
    Q(S) = AMOD( Y(S), 1.0 )
    
Then IX has the following single-indexing subscripts into the grid for the cell-corners surrounding the <X(S),Y(S)>:
    IX(1,S) = C + NC * ( R - 1 ) + 1
    IX(2,S) = C + NC * ( R - 1 ) 
    IX(3,S) = C + NC *   R
    IX(4,S) = C + NC *   R       + 1
    
and AX has the bilinear-interpolation coefficients
    AX(1,S) = ( 1.0 - P( S ) )*( 1.0 - Q( S ) )
    AX(2,S) =         P( S )  *( 1.0 - Q( S ) )
    AX(3,S) = ( 1.0 - P( S ) )*        Q( S )
    AX(4,S) =         P( S )  *        Q( S )
    

See also

GCTP coordinate transformation routine from USGS

MODULE MODGCTP

GRID2INDX()/PNTS2INDX()/INDXMULT(): "new" bilinear-interpolation package in MODULE MODGCTP

and subroutines
BMATVEC and BILIN, DMATVEC, PMATVEC, SMATVEC, and UNGRIDB
and programs
mtxblend, mtxbuild, mtxcalc, mtxcple.

Fortran Usage:

!! under construction !!

Here are (F90 free-source-form) examples of how to do grid-to-points and grid-to-grid interpolation, assuming that you want to interpolate 1-layer variable FOO1, single-indexed 1-layer variablae FOO2, and multilayer variable BAR1 from the grid G1 to the grid G2 (which may be Lambert or something) or to point-sites P. Either BILIN() or BMATVEC() may be used for the single layer interpolation; BILIN() preserves layered subscript order for 3-D variables, whereas BMATVEC() transposes data from I/O API order in which LAYER is the outermost subscript to a computational order in which LAYER is innermost subscript (e.g., as used in SMOKE program LAYPOINT).

Note that UNGRIDB() needs the X and Y coordinates for the G2 cells expressed in terms of the underlying coordinate system for G1. Note that this is backwards of what one might naively expect: In order to compute the matrix to transform data from G1 to G2, one needs to transform the grid-cell coordinates for G2 grid-nodes or point-locations into the coordinate system for G1. For I/O API v3.2 or later, this task of computing G2 locations relative to the coordinate system of G1 may be performed by MODULE MODGCTP subroutines GRID2XY() for grid-node locations, or XY2XY() for point-locations.

    ...
    USE M3UTILIO
    ...
    INTEGER NLAYS
    INTEGER NPNTS
    INTEGER NCOLS1, NROWS1              !  for G1
    REAL*8  XORIG1, YORIG1              !  for G1
    REAL*8  XCELL1, YCELL1              !  for G1
    REAL    FOO1( NCOLS1,NROWS1 )       !  to be interpolated
    REAL    FOO2( NCOLS1*NROWS1 )       !  to be interpolated (single-indexed)
    REAL    BAR1( NCOLS1,NROWS1,NLAYS ) !  to be interpolated
    ...
    REAL*8  XLOCP( NPNTS )      ! PNTS X-values in G1 coords
    REAL*8  YLOCP( NPNTS )      ! PNTS Y-values in G1 coords
    INTEGER NX( 4, NPNTS )      ! interpolation indexes
    REAL    CX( 4, NPNTS )      ! interpolation coeffs
    REAL    FOOP( NPNTS )       ! PNTS-interpolated result
    REAL    BARP( NLAYS,NPNTS ) ! PNTS-transpose/interpolated result
    ...
    INTEGER NCOLS2, NROWS2              !  for G2
    REAL    XLOC2( NCOLS2, NROWS2 )     ! G2 X-values in G1 coords
    REAL    YLAT2( NCOLS2, NROWS2 )     ! G2 Y-values in G1 coords
    INTEGER IX( 4, NCOLS2*NROWS2 )      ! interpolation indexes
    REAL    AX( 4, NCOLS2*NROWS2 )      ! interpolation coeffs
    REAL    FOO3( NCOLS2, NROWS2 )      ! G2-interpolated result
    REAL    FOO4( NCOLS2*NROWS2 )       ! G2-interpolated result
    REAL    BAR2( NCOLS2, NROWS2, NLAYS )   ! G2-interpolated result
    REAL    QUX2( NLAYS, NCOLS2, NROWS2 )   ! G2-transpose-interpolated result
    ...
    !![ read in, calculate, etc., XLOCP, YLOCP, then...]
    ...
    CALL UNGRIDB( NCOLS1, NROWS1,                   &
                  XORIG1, YORIG1, XCELL1, YCELL1,   &
                  NPNTS, XLOCP, YLOCP, NX, CX )
    ...
    CALL BMATVEC( NCOLS1, NROWS1,   &   !  input  size
                  NPNTS,            &   !  output size
                  NX, CX,           &   !  matrix from UNGRIDB()
                  FOO1,             &   !  grid to be interpolated
                  FOOP )                !  PNTS-interpolated result
    ...
    CALL BMATVEC( NCOLS1, NROWS1,   &   !  input  size
                  NPNTS,            &   !  output size
                  NLAYS,            &   !  number of layers
                  NX, CX,           &   !  matrix from UNGRIDB()
                  BAR2,             &   !  to be interpolated
                  BARP )                !  PNTS-transposed/interpolated result
    ...
    CALL GRID2XY( GDTYP1, P_ALP1, P_BET1, P_GAM1, XCENT1, YCENT1,   &
                  GDTYP2, P_ALP2, P_BET2, P_GAM2, XCENT2, YCENT2,   &
                  NCOLS2, NROWS2, XORIG2, YORIG2, XCELL2, YCELL2,   &
                  XLOC2, YLOC2 )

    CALL UNGRIDB( NCOLS1, NROWS1,                   &
                  XORIG1, YORIG1, XCELL1, YCELL1,   &
                  NCOLS2, NROWS2, XLOC2, YLOC2, NX, AX )
    ...
    CALL BMATVEC( NCOLS1, NROWS1,   &   !  input  size
                  NCOLS2, NROWS2,   &   !  output size
                  IX, AX,           &   !  matrix from UNGRIDB()
                  FOO1,             &   !  to be interpolated
                  FOO3 )                !  G2-interpolated result
    ...
    CALL BMATVEC( NCOLS1,NROWS1,    &   !  input  size
                  NCOLS2*NROWS2,    &   !  single-indexed output size
                  IX, AX,           &   !  matrix from UNGRIDB()
                  FOO1,             &   !  to be interpolated
                  FOO4 )                !  G2-interpolated result
    ...
    CALL BILIN( NCOLS1, NROWS1,     &   !  input  size
                NCOLS2,NROWS2,      &   !  output size
                NLAYS,              &   !  number of layers
                IX, AX,             &   !  matrix from UNGRIDB()
                BAR1,               &   !  to be interpolated
                BAR2 )                  !  G2-interpolated result
    ...
    
    CALL BMATVEC( NCOLS1, NROWS1,   &   !  input  size
                  NCOLS2, NROWS2,   &   !  output size
                  NLAYS,            &   !  number of layers
                  IX, AX,           &   !  matrix from UNGRIDB()
                  BAR1,             &   !  to be interpolated
                  QUX2 )                !  G2-transpose-interpolated result
    ...
    

Next: C/Fortran string conversion

Up: Coordinate and Grid Related Routines

Up: Utility Routines

To: Models-3/EDSS I/O API: The Help Pages