For I/O API Version 3.2: BILIN() and BMATVEC() are Fortran-90 generic routines, withINTERFACE
s defined inMODULE M3UTILIO
. The prior versions.For I/O API Versions 3.1 and before, BILIN() and BMATVEC() look like
BILIN11L()
andBMATVEC11()
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
For Fortran-90 declarations and interface checking:USE M3UTILIOBILIN() and BMATVEC() apply 4-band sparse matrix
<IX,AX>
; (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
"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 variablesV
are subscripted by(COLUMN,ROW,LAYER)
, whereas the output interpolated variablesC
are subscripted by(LAYER,STACKNUMBER)
.For single-indexed 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 )ThenIX
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 + 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
GCTP coordinate transformation routine from USGSand subroutines
GRID2INDX()/PNTS2INDX()/INDXMULT()
: "new" bilinear-interpolation package inMODULE MODGCTP
and programsBMATVEC and BILIN
,DMATVEC
,PMATVEC
,SMATVEC
, andUNGRIDB
mtxblend, mtxbuild, mtxcalc, mtxcple.
!! 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 variableFOO1
, single-indexed 1-layer variablaeFOO2
, and multilayer variableBAR1
from the grid G1 to the grid G2 (which may be Lambert or something) or to point-sites P. EitherBILIN()
orBMATVEC()
may be used for the single layer interpolation;BILIN()
preserves layered subscript order for 3-D variables, whereasBMATVEC()
transposes data from I/O API order in whichLAYER
is the outermost subscript to a computational order in whichLAYER
is innermost subscript (e.g., as used in SMOKE program LAYPOINT).Note that
UNGRIDB()
needs theX
andY
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 byMODULE MODGCTP
subroutinesGRID2XY()
for grid-node locations, orXY2XY()
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 ...
Up: Coordinate and Grid Related Routines
To: Models-3/EDSS I/O API: The Help Pages