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
For Fortran-90 declarations and interface checking:USE M3UTILIOBILIN() and BMATVEC() apply 4-band sparse matrix
<NU,CU>
; (e.g., as computed by subroutine UNGRIDB()) to arrayV
, to generate output arrayC
, according to the following equation for output vector value at rowR
and layerL
"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 variablesV
are subscripted by(COLUMN,ROW,LAYER)
, whereas the output interpolated variablesC
are subscripted by(LAYER,STACKNUMBER)
.For bilinear interpolation of gridded data
V
having dimensionsV(NC,NR)
(as specified in the Fortran storage order) to a list of locations having grid-normal coordinates<X(S),Y(S)>
, letC(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 )ThenNU
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 + 1andAX
has the bilinear-interpolation coefficientsAX(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 UNGRIDBand programs
MTXBLEND, MTXBUILD, MTXCALC, MTXCPLE.
!! 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 ...
Up: Coordinate and Grid Related Routines
To: Models-3/EDSS I/O API: The Help Pages