For I/O API Version 3.2: BILIN() and BMATVEC() are Fortran-90 generic routines, withINTERFACEs 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 rowRand 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 variablesVare subscripted by(COLUMN,ROW,LAYER), whereas the output interpolated variablesCare subscripted by(LAYER,STACKNUMBER).For single-indexed bilinear interpolation of gridded data
Vhaving 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 )ThenIXhas 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 + 1andAXhas 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 variableBAR1from 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 whichLAYERis the outermost subscript to a computational order in whichLAYERis innermost subscript (e.g., as used in SMOKE program LAYPOINT).Note that
UNGRIDB()needs theXandYcoordinates 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 MODGCTPsubroutinesGRID2XY()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