BILIN() and BMATVEC()

Fortran version:

    SUBROUTINE BILIN( M, N, P, NU, CU, V, C )
      INTEGER  M         !  input: in-vector length  dimension (single-indexed)
      INTEGER  N         !  input: out-vector length dimension (single-indexed)
      INTEGER  P         !  input:  "layer"          dimension
      INTEGER  NU( 4,N ) !  input:  4-band matrix
      REAL     CU( 4,N ) !  input:  4-band matrix
      REAL     V( M,P )  !  input:  vector to be interpolated
      REAL     C( N,P )  !  output: vector

    SUBROUTINE BMATVEC( M, N, P, NU, CU, V, C )
      INTEGER  M         !  input: in-vector length  dimension (single-indexed)
      INTEGER  N         !  input: out-vector length dimension (single-indexed)
      INTEGER  P         !  input:  "layer"          dimension
      INTEGER  NU( 4,N ) !  input:  4-band matrix
      REAL     CU( 4,N ) !  input:  4-band matrix
      REAL     V( M,P )  !  input:  vector to be interpolated
      REAL     C( P,N )  !  output: *transposed* vector

C version: none

Summary:

For Fortran-90 declarations and interface checking:
    USE M3UTILIO
    

BILIN() and BMATVEC() apply 4-band sparse matrix <NU,CU>; (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"

C( R,L ) = SUMJ = 1...4 [ CU(J,R) V( NU(R,J),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 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 NU has the following single-indexing subscripts into the grid for the cell-corners surrounding the <X(S),Y(S)>:
    NU(1,S) = C + NC * ( R - 1 ) + 1
    NU(2,S) = C + NC * ( R - 1 ) 
    NU(3,S) = C + NC *   R
    NU(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 subroutines

BMATVEC and BILIN, DMATVEC, PMATVEC, SMATVEC, and UNGRIDB
and programs
MTXBLEND, MTXBUILD, MTXCALC, MTXCPLE.

Fortran Usage:

!! under construction !!

Here is an example of how to do grid-to-grid interpolation, assuming that you want to interpolate 1-layer variable FOO1 and multilayer variable BAR1 from the lat-lon-based grid G1 to the grid G2 (which may be Lambert or something). 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.

Note that UNGRIDB needs the X and Y coordinates for the G2 cells expressed in terms of the underlying coordinate system for G1 (in this case, since G1 is lat-lon, this means you need the G2-cell latitudes and longitudes as inputs for UNGRIDB).

    ...
    USE M3UTILIO
    ...
    INTEGER	NCOLS1, NROWS1			!  for G1
    INTEGER     P				!  number of layers is P==1
    INTEGER     NLAYS
    REAL*8	XORIG1, YORIG1			!  for G1
    REAL*8	XCELL1, YCELL1			!  for G1
    REAL	FOO1( NCOLS1, NROWS1 )		!  to be interpolated
    REAL	BAR1( NCOLS1, NROWS1, NLAYS )	!  to be interpolated
    ...
    INTEGER	NCOLS2, NROWS2	!  for G2
    REAL	ALON2( NCOLS2, NROWS2 )		! G2 X-values in G1 coords
    REAL	ALAT2( NCOLS2, NROWS2 )		! G2 Y-values in G1 coords
    REAL	CU( 4, NCOLS2*NROWS2 )  	! interpolation coeffs
    INTEGER	NU( 4, NCOLS2*NROWS2 )  	! interpolation indexes
    REAL	FOO2( NCOLS2, NROWS2 )		! G2-interpolated result
    REAL	BAR2( NCOLS2, NROWS2, NLAYS )	! G2-interpolated result
    REAL	QUX3( NLAYS, NCOLS2, NROWS2 )	! transpose-interpolated result
    ...
    CALL UNGRIDB( NCOLS1, NROWS1, XORIG1, YORIG1, XCELL1, YCELL1,
   &              NCOLS2*NROWS2, ALON2, ALAT2, NU, CU )
    ...
    P = 1	!  only one layer
    CALL BMATVEC( NCOLS1*NROWS1,	!  single-indexed input  size
   &              NCOLS2*NROWS2,	!  single-indexed output size
   &              P,			!  P=1 layer
   &              NU, CU,		!  matrix from UNGRIDB()
   &              FOO1,			!  to be interpolated
   &		  FOO2 )		!  G2-interpolated result
    ...		
    CALL BILIN( NCOLS1*NROWS1,	        !  single-indexed input  size
   &            NCOLS2*NROWS2,	        !  single-indexed output size
   &            NLAYS,		        !  number of layers
   &            NU, CU,		        !  matrix from UNGRIDB()
   &            BAR1,			!  to be interpolated
   &		BAR2 )		        !  G2-interpolated result
    ...		
    
    CALL BMATVEC( NCOLS1*NROWS1,	!  single-indexed input  size
   &              NCOLS2*NROWS2,	!  single-indexed output size
   &              NLAYS,		!  number of layers
   &              NU, CU,		!  matrix from UNGRIDB()
   &              BAR1,			!  to be interpolated
   &		  QUX3 )	        !  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