GAUSSLAT

The GAUSSLAT routine computes cosines of colatitude and Gaussian
             weights on the Gaussian latitudes.  The Gaussian latitudes
             are at the zeroes of the legendre polynomial of the given order.

USAGE:    CALL GAUSSLAT(JMAX,SLAT,WLAT)

   Input argument list:
      JMAX     - INPUT NUMBER OF LATITUDES.
 
   Output argument list:
      SLAT     - REAL (K) COSINES OF COLATITUDE.
      WLAT     - REAL (K) GAUSSIAN WEIGHTS.


GCDIST

The GCDIST function computes great circle distance
           between two points on the earth.

USAGE:    ...GCDIST(RLAT1,RLON1,RLAT2,RLON2)

   Input argument list:
      RLAT1    - REAL LATITUDE OF POINT 1 IN DEGREES
      RLON1    - REAL LONGITUDE OF POINT 1 IN DEGREES
      RLAT2    - REAL LATITUDE OF POINT 2 IN DEGREES
      RLON2    - REAL LONGITUDE OF POINT 2 IN DEGREES
 
   Output argument list:
      GCDIST   - REAL GREAT CIRCLE DISTANCE IN METERS


GDSAWT

The GDSAWT routine decodes the grib grid description section
           and returns the earth's area covered by each grid point.
           The grid and earth coordinates are also returned. The
           edge points and corner points are given full weight now.

USAGE:    CALL GDSAWT(KGDS,KB,KA,XPTS,YPTS,RLAT,RLON,AWT)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
      KB       - INTEGER MAXIMUM NUMBER OF COORDINATES

   Output argument list:
      KA       - INTEGER NUMBER OF VALID POINTS COMPUTED
      XPTS     - REAL (KB) GRID X POINT COORDINATES
      YPTS     - REAL (KB) GRID Y POINT COORDINATES
      RLON     - REAL (KB) EARTH LONGITUDES IN DEGREES
      RLAT     - REAL (KB) EARTH LATITUDES IN DEGREES N
      AWT      - REAL (KB) AREA WEIGHTS IN M**2


GDSWIZ

The GDSWIZ routine decodes the grib grid description section
           (passed in integer form as decoded by subprogram w3fi63)
           and returns one of the following:
             (IOPT= 0) grid and earth coordinates of all grid points
             (IOPT=+1) earth coordinates of selected grid coordinates
             (IOPT=-1) grid coordinates of selected earth coordinates

           The current code recognizes the following projections:
             (KGDS(1)=000) Equidistant cylindrical
             (KGDS(1)=001) Mercator cylindrical
             (KGDS(1)=003) Lambert conformal conical
             (KGDS(1)=004) Gaussian cylindrical
             (KGDS(1)=005) Polar stereographic azimuthaL
             (KGDS(1)=201) Staggered rotated equidistant cylindrical
             (KGDS(1)=202) Rotated equidistant cylindrical
             (KGDS(1)=203) Staggered rotated equidistant cylindrical 2-D

           If the selected coordinates are more than one gridpoint
           beyond the the edges of the grid domain, then the relevant
           output elements are set to fill values.  Also if iopt=0,
           if the number of grid points exceeds the number allotted,
           then all the output elements are set to fill values.
           The actual number of valid points computed is returned too.

USAGE:    CALL GDSWIZ(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                      LROT,CROT,SROT)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
      IOPT     - INTEGER OPTION FLAG
                 ( 0 TO COMPUTE EARTH COORDS OF ALL THE GRID POINTS)
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<=0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<=0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>=0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>=0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
                 (-1 IF PROJECTION UNRECOGNIZED)
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
 

GDSWIZ00

The GDSWIZ00 routine decodes the grib grid description section
             (passed in integer form as decoded by subprogram w3fi63)
             and returns one of the following:
               (IOPT=+1) earth coordinates of selected grid coordinates
               (IOPT=-1) grid coordinates of selected earth coordinates

             For equidistant cylindrical projections.
             If the selected coordinates are more than one gridpoint
             beyond the the edges of the grid domain, then the relevant
             output elements are set to fill values.
             The actual number of valid points computed is returned too.

USAGE:    CALL GDSWIZ00(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                        LROT,CROT,SROT)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
      IOPT     - INTEGER OPTION FLAG
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
 

GDSWIZ01

The GDSWIZ01 decodes the grib grid description section
             (passed in integer form as decoded by subprogram w3fi63)
             and returns one of the following:
               (IOPT=+1) earth coordinates of selected grid coordinates
               (IOPT=-1) grid coordinates of selected earth coordinates

             For Mercator cylindrical projections.
             If the selected coordinates are more than one gridpoint
             beyond the the edges of the grid domain, then the relevant
             output elements are set to fill values.
             The actual number of valid points computed is returned too.

USAGE:    CALL GDSWIZ01(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                        LROT,CROT,SROT)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
      IOPT     - INTEGER OPTION FLAG
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
 

GDSWIZ03

The GDSWIZ03 routine decodes the grib grid description section
             (passed in integer form as decoded by subprogram w3fi63)
             and returns one of the following:
               (IOPT=+1) earth coordinates of selected grid coordinates
               (IOPT=-1) grid coordinates of selected earth coordinates

             For Lambert Conformal conical projections.
             If the selected coordinates are more than one gridpoint
             beyond the the edges of the grid domain, then the relevant
             output elements are set to fill values.
             The actual number of valid points computed is returned too.

USAGE:    CALL GDSWIZ03(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                        LROT,CROT,SROT)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
      IOPT     - INTEGER OPTION FLAG
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
 

GDSWIZ04

The GDSWIZ04 routine decodes the grib grid description section
             (passed in integer form as decoded by subprogram w3fi63)
             and returns one of the following:
               (IOPT=+1) earth coordinates of selected grid coordinates
               (IOPT=-1) grid coordinates of selected earth coordinates

             For Gaussian cylindrical projections.
             If the selected coordinates are more than one gridpoint
             beyond the the edges of the grid domain, then the relevant
             output elements are set to fill values.
             The actual number of valid points computed is returned too.

USAGE:    CALL GDSWIZ04(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                        LROT,CROT,SROT)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
      IOPT     - INTEGER OPTION FLAG
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)


GDSWIZ05

The GDSWIZ05 routine decodes the grib grid description section
             (passed in integer form as decoded by subprogram w3fi63)
             and returns one of the following:
               (IOPT=+1) earth coordinates of selected grid coordinates
               (IOPT=-1) grid coordinates of selected earth coordinates

             For Polar Stereographic azimuthal projections.
             If the selected coordinates are more than one gridpoint
             beyond the the edges of the grid domain, then the relevant
             output elements are set to fill values.
             The actual number of valid points computed is returned too.

USAGE:    CALL GDSWIZ05(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                        LROT,CROT,SROT)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
      IOPT     - INTEGER OPTION FLAG
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
 

GDSWIZC9

The GDSWIZC9 routine DECODES THE GRIB GRID DESCRIPTION SECTION
             (passed in integer form as decoded by subprogram w3fi63)
             and returns one of the following:
               (IOPT=+1) earth coordinates of selected grid coordinates
               (IOPT=-1) grid coordinates of selected earth coordinates

             For Staggered rotated equidistant cylindrical projections.
             (see under the description of kgds to determine whether
             to compute a Staggered wind grid or a staggered mass grid.)
             If the selected coordinates are more than one gridpoint
             beyond the the edges of the grid domain, then the relevant
             output elements are set to fill values.
             The actual number of valid points computed is returned too.

USAGE:    CALL GDSWIZC9(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                        LROT,CROT,SROT)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
                 IMPORTANT NOTE: IF THE 9TH BIT (FROM RIGHT) OF KGDS(11)
                                 (SCANNING MODE FLAG) IS 1, THEN THIS
                                 THE GRID IS COMPUTED FOR A WIND FIELD;
                                 OTHERWISE IT IS FOR A MASS FIELD.  THUS
                                 MOD(KGDS(11)/256,2)=0 FOR MASS GRID.
      IOPT     - INTEGER OPTION FLAG
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
 

GDSWIZCA

The GDSWIZCA routine decodes the grib grid description section
             (passed in integer form as decoded by subprogram w3fi63)
             and returns one of the following:
                (IOPT=+1) earth coordinates of selected grid coordinates
                (IOPT=-1) grid coordinates of selected earth coordinates

             For Rotated Equidistant cylindrical projections.
             If the selected coordinates are more than one gridpoint
             beyond the the edges of the grid domain, then the relevant
             output elements are set to fill values.
             The actual number of valid points computed is returned too.

USAGE:    CALL GDSWIZCA(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                        LROT,CROT,SROT)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
      IOPT     - INTEGER OPTION FLAG
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
 

GDSWIZCB

The GDSWIZCB routine decodes the grib grid description section
             (passed in integer form as decoded by subprogram w3fi63)
             and returns one of the following:
               (IOPT=+1) earth coordinates of selected grid coordinates
               (IOPT=-1) grid coordinates of selected earth coordinates

             For Staggered rotated equidistant cylindrical projections.
             (see under the description of kgds to determine whether
             to compute a Staggered wind grid or a staggered mass grid.)
             If the selected coordinates are more than one gridpoint
             beyond the the edges of the grid domain, then the relevant
             output elements are set to fill values.
             The actual number of valid points computed is returned too.

USAGE:    CALL GDSWIZCB(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                        LROT,CROT,SROT)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
                 IMPORTANT NOTE: IF THE 9TH BIT (FROM RIGHT) OF KGDS(11)
                                 (SCANNING MODE FLAG) IS 1, THEN THIS
                                 THE GRID IS COMPUTED FOR A WIND FIELD;
                                 OTHERWISE IT IS FOR A MASS FIELD.  THUS
                                 MOD(KGDS(11)/256,2)=0 FOR MASS GRID.
      IOPT     - INTEGER OPTION FLAG
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
 

GDSWZD

The GDSWZD routine decodes the grib grid description section
           (passed in integer form as decoded by subprogram w3fi63)
           and returns one of the following:
             (IOPT= 0) grid and earth coordinates of all grid points
             (IOPT=+1) earth coordinates of selected grid coordinates
             (IOPT=-1) grid coordinates of selected earth coordinates

           The current code recognizes the following projections:
             (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
             (KGDS(1)=001) MERCATOR CYLINDRICAL
             (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
             (KGDS(1)=004) GAUSSIAN CYLINDRICAL
             (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
             (KGDS(1)=201) STAGGERED ROTATED EQUIDISTANT CYLINDRICAL
             (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL
             (KGDS(1)=203) STAGGERED ROTATED EQUIDISTANT CYLINDRICAL 2-D

           If the selected coordinates are more than one gridpoint
           beyond the the edges of the grid domain, then the relevant
           output elements are set to fill values.  Also if IOPT=0,
           If the number of grid points exceeds the number allotted,
           then all the output elements are set to fill values.
           The actual number of valid points computed is returned too.
           Optionally, the vector rotations and the map jacobians
           for this grid may be returned as well.

USAGE:    CALL GDSWZD(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                      LROT,CROT,SROT,LMAP,XLON,XLAT,YLON,YLAT,AREA)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
      IOPT     - INTEGER OPTION FLAG
                 ( 0 TO COMPUTE EARTH COORDS OF ALL THE GRID POINTS)
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
      LMAP     - INTEGER FLAG TO RETURN MAP JACOBIANS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<=0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<=0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>=0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>=0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
                 (-1 IF PROJECTION UNRECOGNIZED)
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
      XLON     - REAL (NPTS) DX/DLON IN 1/DEGREES IF LMAP=1
      XLAT     - REAL (NPTS) DX/DLAT IN 1/DEGREES IF LMAP=1
      YLON     - REAL (NPTS) DY/DLON IN 1/DEGREES IF LMAP=1
      YLAT     - REAL (NPTS) DY/DLAT IN 1/DEGREES IF LMAP=1
      AREA     - REAL (NPTS) AREA WEIGHTS IN M**2 IF LMAP=1
                 (PROPORTIONAL TO THE SQUARE OF THE MAP FACTOR
                  IN THE CASE OF CONFORMAL PROJECTIONS)
 

GDSWZD00

The GDSWZD00 routine decodes the grib grid description section
             (passed in integer form as decoded by subprogram w3fi63)
             and returns one of the following:
               (IOPT=+1) earth coordinates of selected grid coordinates
               (IOPT=-1) grid coordinates of selected earth coordinates

             For Equidistant cylindrical projections.
             If the selected coordinates are more than one gridpoint
             beyond the the edges of the grid domain, then the relevant
             output elements are set to fill values.
             The actual number of valid points computed is returned too.
             optionally, the vector rotations and the map Jacobians
             for this grid may be returned as well.

USAGE:    CALL GDSWZD00(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                     LROT,CROT,SROT,LMAP,XLON,XLAT,YLON,YLAT,AREA)

   INPUT ARGUMENT LIST:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
      IOPT     - INTEGER OPTION FLAG
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
      LMAP     - INTEGER FLAG TO RETURN MAP JACOBIANS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
      XLON     - REAL (NPTS) DX/DLON IN 1/DEGREES IF LMAP=1
      XLAT     - REAL (NPTS) DX/DLAT IN 1/DEGREES IF LMAP=1
      YLON     - REAL (NPTS) DY/DLON IN 1/DEGREES IF LMAP=1
      YLAT     - REAL (NPTS) DY/DLAT IN 1/DEGREES IF LMAP=1
      AREA     - REAL (NPTS) AREA WEIGHTS IN M**2 IF LMAP=1
 

GDSWZD01

The GDSWZD01 routine decodes the grib grid description section
             (passed in integer form as decoded by subprogram w3fi63)
             and returns one of the following:
               (IOPT=+1) earth coordinates of selected grid coordinates
               (IOPT=-1) grid coordinates of selected earth coordinates

             For Mercator cylindrical projections.
             If the selected coordinates are more than one gridpoint
             beyond the the edges of the grid domain, then the relevant
             output elements are set to fill values.
             The actual number of valid points computed is returned too.
             Optionally, the vector rotations and the map Jacobians
             for this grid may be returned as well.

USAGE:    CALL GDSWZD01(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                        LROT,CROT,SROT,LMAP,XLON,XLAT,YLON,YLAT,AREA)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
      IOPT     - INTEGER OPTION FLAG
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
      LMAP     - INTEGER FLAG TO RETURN MAP JACOBIANS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
      XLON     - REAL (NPTS) DX/DLON IN 1/DEGREES IF LMAP=1
      XLAT     - REAL (NPTS) DX/DLAT IN 1/DEGREES IF LMAP=1
      YLON     - REAL (NPTS) DY/DLON IN 1/DEGREES IF LMAP=1
      YLAT     - REAL (NPTS) DY/DLAT IN 1/DEGREES IF LMAP=1
      AREA     - REAL (NPTS) AREA WEIGHTS IN M**2 IF LMAP=1
                 (PROPORTIONAL TO THE SQUARE OF THE MAP FACTOR)
 

GDSWZD03

The GDSWZD03 routine decodes the grib grid description section
             (passed in integer form as decoded by subprogram w3fi63)
             and returns one of the following:
               (IOPT=+1) earth coordinates of selected grid coordinates
               (IOPT=-1) grid coordinates of selected earth coordinates

             For Lambert Conformal conical projections.
             If the selected coordinates are more than one gridpoint
             beyond the the edges of the grid domain, then the relevant
             output elements are set to fill values.
             The actual number of valid points computed is returned too.
             optionally, the vector rotations and the map Jacobians
             For this grid may be returned as well.

USAGE:    CALL GDSWZD03(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                        LROT,CROT,SROT,LMAP,XLON,XLAT,YLON,YLAT,AREA)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
      IOPT     - INTEGER OPTION FLAG
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
      LMAP     - INTEGER FLAG TO RETURN MAP JACOBIANS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
      XLON     - REAL (NPTS) DX/DLON IN 1/DEGREES IF LMAP=1
      XLAT     - REAL (NPTS) DX/DLAT IN 1/DEGREES IF LMAP=1
      YLON     - REAL (NPTS) DY/DLON IN 1/DEGREES IF LMAP=1
      YLAT     - REAL (NPTS) DY/DLAT IN 1/DEGREES IF LMAP=1
      AREA     - REAL (NPTS) AREA WEIGHTS IN M**2 IF LMAP=1
                 (PROPORTIONAL TO THE SQUARE OF THE MAP FACTOR)
 

GDSWZD04

The GDSWZD04 routine decodes the grib grid description section
             (passed in integer form as decoded by subprogram w3fi63)
             and returns one of the following:
               (IOPT=+1) earth coordinates of selected grid coordinates
               (IOPT=-1) grid coordinates of selected earth coordinates

             For Gaussian cylindrical projections.
             If the selected coordinates are more than one gridpoint
             beyond the the edges of the grid domain, then the relevant
             output elements are set to fill values.
             the actual number of valid points computed is returned too.
             Optionally, the vector rotations and the map Jacobians
             for this grid may be returned as well.

USAGE:    CALL GDSWZD04(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                        LROT,CROT,SROT,LMAP,XLON,XLAT,YLON,YLAT,AREA)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
      IOPT     - INTEGER OPTION FLAG
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
      LMAP     - INTEGER FLAG TO RETURN MAP JACOBIANS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
      XLON     - REAL (NPTS) DX/DLON IN 1/DEGREES IF LMAP=1
      XLAT     - REAL (NPTS) DX/DLAT IN 1/DEGREES IF LMAP=1
      YLON     - REAL (NPTS) DY/DLON IN 1/DEGREES IF LMAP=1
      YLAT     - REAL (NPTS) DY/DLAT IN 1/DEGREES IF LMAP=1
      AREA     - REAL (NPTS) AREA WEIGHTS IN M**2 IF LMAP=1
 

GDSWZD05

The GDSWZD05 routine decodes the grib grid description section
             (passed in integer form as decoded by subprogram w3fi63)
             and returns one of the following:
               (IOPT=+1) earth coordinates of selected grid coordinates
               (IOPT=-1) grid coordinates of selected earth coordinates

             For Polar Stereographic Azimuthal projections.
             If the selected coordinates are more than one gridpoint
             beyond the the edges of the grid domain, then the relevant
             output elements are set to fill values.
             The actual number of valid points computed is returned too.
             Optionally, the vector rotations and the map Jacobians
             for this grid may be returned as well.

USAGE:    CALL GDSWZD05(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                        LROT,CROT,SROT,LMAP,XLON,XLAT,YLON,YLAT,AREA)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
      IOPT     - INTEGER OPTION FLAG
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
      LMAP     - INTEGER FLAG TO RETURN MAP JACOBIANS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
      XLON     - REAL (NPTS) DX/DLON IN 1/DEGREES IF LMAP=1
      XLAT     - REAL (NPTS) DX/DLAT IN 1/DEGREES IF LMAP=1
      YLON     - REAL (NPTS) DY/DLON IN 1/DEGREES IF LMAP=1
      YLAT     - REAL (NPTS) DY/DLAT IN 1/DEGREES IF LMAP=1
      AREA     - REAL (NPTS) AREA WEIGHTS IN M**2 IF LMAP=1
                 (PROPORTIONAL TO THE SQUARE OF THE MAP FACTOR)
 

GDSWZDC9

The GDSWZDC9 routine decodes the grib grid description section
             (passed in integer form as decoded by subprogram w3fi63)
             and returns one of the following:
               (IOPT=+1) earth coordinates of selected grid coordinates
               (IOPT=-1) grid coordinates of selected earth coordinates

             For Staggered Rotated Equidistant cylindrical projections.
             (See under the description of kgds to determine whether
             to compute a Staggered wind grid or a Staggered mass grid.)
             If the selected coordinates are more than one gridpoint
             beyond the the edges of the grid domain, then the relevant
             output elements are set to fill values.
             The actual number of valid points computed is returned too.
             Optionally, the vector rotations and the map Jacobians
             for this grid may be returned as well.

USAGE:    CALL GDSWZDC9(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                        LROT,CROT,SROT,LMAP,XLON,XLAT,YLON,YLAT,AREA)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
                 IMPORTANT NOTE: IF THE 9TH BIT (FROM RIGHT) OF KGDS(11)
                                 (SCANNING MODE FLAG) IS 1, THEN THIS
                                 THE GRID IS COMPUTED FOR A WIND FIELD;
                                 OTHERWISE IT IS FOR A MASS FIELD.  THUS
                                 MOD(KGDS(11)/256,2)=0 FOR MASS GRID.
      IOPT     - INTEGER OPTION FLAG
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
      LMAP     - INTEGER FLAG TO RETURN MAP JACOBIANS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
      XLON     - REAL (NPTS) DX/DLON IN 1/DEGREES IF LMAP=1
      XLAT     - REAL (NPTS) DX/DLAT IN 1/DEGREES IF LMAP=1
      YLON     - REAL (NPTS) DY/DLON IN 1/DEGREES IF LMAP=1
      YLAT     - REAL (NPTS) DY/DLAT IN 1/DEGREES IF LMAP=1
      AREA     - REAL (NPTS) AREA WEIGHTS IN M**2 IF LMAP=1


GDSWZDCA

The GDSWZDCA routine decodes the grib grid description section
             (passed in integer form as decoded by subprogram w3fi63)
             and returns one of the following:
               (IOPT=+1) earth coordinates of selected grid coordinates
               (IOPT=-1) grid coordinates of selected earth coordinates

             For Rotated Equidistant cylindrical projections.
             If the selected coordinates are more than one gridpoint
             beyond the the edges of the grid domain, then the relevant
             output elements are set to fill values.
             The actual number of valid points computed is returned too.
             Optionally, the vector rotations and the map Jacobians
             for this grid may be returned as well.

USAGE:    CALL GDSWZDCA(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                    LROT,CROT,SROT,LMAP,XLON,XLAT,YLON,YLAT,AREA)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
      IOPT     - INTEGER OPTION FLAG
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
      LMAP     - INTEGER FLAG TO RETURN MAP JACOBIANS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
      XLON     - REAL (NPTS) DX/DLON IN 1/DEGREES IF LMAP=1
      XLAT     - REAL (NPTS) DX/DLAT IN 1/DEGREES IF LMAP=1
      YLON     - REAL (NPTS) DY/DLON IN 1/DEGREES IF LMAP=1
      YLAT     - REAL (NPTS) DY/DLAT IN 1/DEGREES IF LMAP=1
      AREA     - REAL (NPTS) AREA WEIGHTS IN M**2 IF LMAP=1


GDSWZDCB

The GDSWZDCB routine decodes the grib grid description section
             (passed in integer form as decoded by subprogram w3fi63)
             and returns one of the following:
               (IOPT=+1) earth coordinates of selected grid coordinates
               (IOPT=-1) grid coordinates of selected earth coordinates

             For Staggered Rotated Equidistant cylindrical projections.
             (See under the description of kgds to determine whether
             to compute a Staggered wind grid or a Staggered mass grid.)
             If the selected coordinates are more than one gridpoint
             beyond the the edges of the grid domain, then the relevant
             output elements are set to fill values.
             The actual number of valid points computed is returned too.
             Optionally, the vector rotations and the map Jacobians
             for this grid may be returned as well.

USAGE:    CALL GDSWZDCB(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,
                        LROT,CROT,SROT,LMAP,XLON,XLAT,YLON,YLAT,AREA)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
                 IMPORTANT NOTE: IF THE 9TH BIT (FROM RIGHT) OF KGDS(11)
                                 (SCANNING MODE FLAG) IS 1, THEN THIS
                                 THE GRID IS COMPUTED FOR A WIND FIELD;
                                 OTHERWISE IT IS FOR A MASS FIELD.  THUS
                                 MOD(KGDS(11)/256,2)=0 FOR MASS GRID.
      IOPT     - INTEGER OPTION FLAG
                 (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
                 (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
      NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
      FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
                 (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
                 (ACCEPTABLE RANGE: -360. TO 360.)
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
                 (ACCEPTABLE RANGE: -90. TO 90.)
      LROT     - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1
      LMAP     - INTEGER FLAG TO RETURN MAP JACOBIANS IF 1
 
   Output argument list:
      XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
      YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
      RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
      RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
      NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
      CROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1
      SROT     - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
      XLON     - REAL (NPTS) DX/DLON IN 1/DEGREES IF LMAP=1
      XLAT     - REAL (NPTS) DX/DLAT IN 1/DEGREES IF LMAP=1
      YLON     - REAL (NPTS) DY/DLON IN 1/DEGREES IF LMAP=1
      YLAT     - REAL (NPTS) DY/DLAT IN 1/DEGREES IF LMAP=1
      AREA     - REAL (NPTS) AREA WEIGHTS IN M**2 IF LMAP=1
 

IJKGDS

The IJKGDS function decodes the grib grid description section
           and returns the field position for a given grid point.

USAGE:    ...IJKGDS(I,J,KGDS)

   Input argument list:
      I        - INTEGER X GRID POINT
      J        - INTEGER Y GRID POINT
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
 
   Output argument list:
      IJKGDS   - INTEGER POSITION IN GRIB FIELD TO LOCATE GRID POINT
                 (0 IF OUT OF BOUNDS)


IJKGDS0

The IJKGDS0 routine decodes the grib grid description section
            and returns a navigation parameter array to allow function
            IJKGDS1 to decode the field position for a given grid point.

USAGE:    CALL IJKGDS0(KGDS,IJKGDSA)

   Input argument list:
      KGDS     - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63
 
   Output argument list:
      IJKGDSA  - INTEGER (20) NAVIGATION PARAMETER ARRAY
                 IJKGDSA(1) IS NUMBER OF X POINTS
                 IJKGDSA(2) IS NUMBER OF Y POINTS
                 IJKGDSA(3) IS X WRAPAROUND INCREMENT
                            (0 IF NO WRAPAROUND)
                 IJKGDSA(4) IS Y WRAPAROUND LOWER PIVOT POINT
                            (0 IF NO WRAPAROUND)
                 IJKGDSA(5) IS Y WRAPAROUND UPPER PIVOT POINT
                            (0 IF NO WRAPAROUND)
                 IJKGDSA(6) IS SCANNING MODE
                            (0 IF X FIRST THEN Y; 1 IF Y FIRST THEN X;
                             2 IF STAGGERED DIAGONAL LIKE PROJECTION 201;
                             3 IF STAGGERED DIAGONAL LIKE PROJECTION 203)
                 IJKGDSA(7) IS MASS/WIND FLAG FOR STAGGERED DIAGONAL
                            (0 IF MASS; 1 IF WIND)
                 IJKGDSA(8:20) ARE UNUSED AT THE MOMENT
 

IJKGDS1

The IJKGDS1 function decodes the grib grid description section
            and returns the field position for a given grid point.
            call IJKGDS0 to set up the navigation parameter array.

USAGE:    ...IJKGDS1(I,J,IJKGDSA)

   Input argument list:
      I        - INTEGER X GRID POINT
      J        - INTEGER Y GRID POINT
      IJKGDSA  - INTEGER (20) NAVIGATION PARAMETER ARRAY
                 IJKGDSA(1) IS NUMBER OF X POINTS
                 IJKGDSA(2) IS NUMBER OF Y POINTS
                 IJKGDSA(3) IS X WRAPAROUND INCREMENT
                            (0 IF NO WRAPAROUND)
                 IJKGDSA(4) IS Y WRAPAROUND LOWER PIVOT POINT
                            (0 IF NO WRAPAROUND)
                 IJKGDSA(5) IS Y WRAPAROUND UPPER PIVOT POINT
                            (0 IF NO WRAPAROUND)
                 IJKGDSA(6) IS SCANNING MODE
                            (0 IF X FIRST THEN Y; 1 IF Y FIRST THEN X;
                             2 IF STAGGERED DIAGONAL LIKE PROJECTION 201;
                             3 IF STAGGERED DIAGONAL LIKE PROJECTION 203)
                 IJKGDSA(7) IS MASS/WIND FLAG FOR STAGGERED DIAGONAL
                            (0 IF MASS; 1 IF WIND)
                 IJKGDSA(8:20) ARE UNUSED AT THE MOMENT
 
   Output argument list:
      IJKGDS   - INTEGER POSITION IN GRIB FIELD TO LOCATE GRID POINT
                 (0 IF OUT OF BOUNDS)
 

IPMERGE2

The IPMERGE2 routine merges two bitmapped fields into one.
             where both input fields are valid, the first field is taken.
             where neither input field is valid, the output is invalid.

USAGE:    CALL IPMERGE2(NO,NF,M1,L1,F1,M2,L2,F2,MO,LO,FO)

   Input argument list:
      NO       - INTEGER NUMBER OF POINTS IN EACH FIELD
      NF       - INTEGER NUMBER OF FIELDS MERGES
      M1       - INTEGER FIRST DIMENSION OF FIELD 1 ARRAYS
      L1       - LOGICAL(1) (M1,NF) BITMAP FOR FIELD 1
      F1       - REAL (M1,NF) DATA FOR FIELD 1
      M2       - INTEGER FIRST DIMENSION OF FIELD 2 ARRAYS
      L2       - LOGICAL(1) (M2,NF) BITMAP FOR FIELD 2
      F2       - REAL (M2,NF) DATA FOR FIELD 2
      MO       - INTEGER FIRST DIMENSION OF OUTPUT FIELD ARRAYS
 
   Output argument list:
      LO       - LOGICAL(1) (MO,NF) BITMAP FOR OUTPUT FIELD
      FO       - REAL (MO,NF) DATA FOR OUTPUT FIELD


IPOLATES

The IPOLATES routine interpolates scalar fields
             from any grid to any grid.  Only horizontal
             interpolation is performed.

             The following interpolation methods are possible:
                (IP=0) BILINEAR
                (IP=1) BICUBIC
                (IP=2) NEIGHBOR
                (IP=3) BUDGET
                (IP=4) SPECTRAL
                (IP=6) NEIGHBOR-BUDGET

             Some of these methods have interpolation options and/or
             restrictions on the input or output grids, both of which
             are documented more fully in their respective subprograms.
             The grids are defined by their grid description sections
             (passed in integer form as decoded by subprogram w3fi63).

             The current code recognizes the following projections:
                (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
                (KGDS(1)=001) MERCATOR CYLINDRICAL
                (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
                (KGDS(1)=004) GAUSSIAN CYLINDRICAL
                (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
                (KGDS(1)=201) ROTATED EQUIDISTANT CYLINDRICAL
                (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL

             Where KGDS could be either input KGDSI or output KGDSO.
             As an added bonus the number of output grid points
             and their latitudes and longitudes are also returned.
             On the other hand, the output can be a set of station points
             If KGDSO(1)<0, in which case the number of points
             and their latitudes and longitudes must be input.
             input bitmaps will be interpolated to output bitmaps.
             output bitmaps will also be created when the output grid
             extends outside of the domain of the input grid.
             The output field is set to 0 where the output bitmap is off.
 
USAGE:    CALL IPOLATES(IP,IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
                        NO,RLAT,RLON,IBO,LO,GO,IRET)
 
   Input argument list:
      IP       - INTEGER INTERPOLATION METHOD
                 (IP=0 FOR BILINEAR;
                  IP=1 FOR BICUBIC;
                  IP=2 FOR NEIGHBOR;
                  IP=3 FOR BUDGET;
                  IP=4 FOR SPECTRAL;
                  IP=6 FOR NEIGHBOR-BUDGET)
      IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
                 (IP=0: (NO OPTIONS)
                  IP=1: CONSTRAINT OPTION
                  IP=2: (NO OPTIONS)
                  IP=3: NUMBER IN RADIUS, RADIUS WEIGHTS ...
                  IP=4: SPECTRAL SHAPE, SPECTRAL TRUNCATION
                  IP=6: NUMBER IN RADIUS, RADIUS WEIGHTS ...)
      KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
      KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
      MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF INPUT GRID FIELDS IF KM=1
      MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
      KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
      IBI      - INTEGER (KM) INPUT BITMAP FLAGS
      LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF RESPECTIVE IBI(K)=1)
      GI       - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)
      RLAT     - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
      RLON     - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
 
   Output argument list:
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
      RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
      RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
      IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
      LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
      GO       - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL INTERPOLATION
                 1    UNRECOGNIZED INTERPOLATION METHOD
                 2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
                 3    UNRECOGNIZED OUTPUT GRID
                 1X   INVALID BICUBIC METHOD PARAMETERS
                 3X   INVALID BUDGET METHOD PARAMETERS
                 4X   INVALID SPECTRAL METHOD PARAMETERS
 

IPOLATEV

The IPOLATEV routine interpolates vector fields
             from any grid to any grid.  Only horizontal
             interpolation is performed.

             The following interpolation methods are possible:
                (IP=0) BILINEAR
                (IP=1) BICUBIC
                (IP=2) NEIGHBOR
                (IP=3) BUDGET
                (IP=4) SPECTRAL
                (IP=6) NEIGHBOR-BUDGET

             Some of these methods have interpolation options and/or
             restrictions on the input or output grids, both of which
             are documented more fully in their respective subprograms.
             The grids are defined by their grid description sections
             (passed in integer form as decoded by subprogram w3fi63).

             The current code recognizes the following projections:
                (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
                (KGDS(1)=001) MERCATOR CYLINDRICAL
                (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
                (KGDS(1)=004) GAUSSIAN CYLINDRICAL
                (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
                (KGDS(1)=201) ROTATED EQUIDISTANT CYLINDRICAL
                (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL

             Where KGDS could be either input KGDSI or output KGDSO.
             As an added bonus the number of output grid points
             and their latitudes and longitudes are also returned.
             On the other hand, the output can be a set of station points
             If KGDSO(1)<0, in which case the number of points
             and their latitudes and longitudes must be input.
             input bitmaps will be interpolated to output bitmaps.
             output bitmaps will also be created when the output grid
             extends outside of the domain of the input grid.
             The output field is set to 0 where the output bitmap is off.

USAGE:    CALL IPOLATEV(IP,IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,UI,VI,
                        NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)

   Input argument list:
      IP       - INTEGER INTERPOLATION METHOD
                 (IP=0 FOR BILINEAR;
                  IP=1 FOR BICUBIC;
                  IP=2 FOR NEIGHBOR;
                  IP=3 FOR BUDGET;
                  IP=4 FOR SPECTRAL;
                  IP=6 FOR NEIGHBOR-BUDGET)
      IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
                 (IP=0: (NO OPTIONS)
                  IP=1: CONSTRAINT OPTION
                  IP=2: (NO OPTIONS)
                  IP=3: NUMBER IN RADIUS, RADIUS WEIGHTS ...
                  IP=4: SPECTRAL SHAPE, SPECTRAL TRUNCATION
                  IP=6: NUMBER IN RADIUS, RADIUS WEIGHTS ...)
      KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
                 NOTE: IF KGDSI(201)=1 THEN THE 9TH BIT OF KGDSI(11)
                       IS TEMPORARILY SET TO 1 TO ALERT THE GDS WIZARD
                       THAT THESE FIELDS ARE STAGGERED ETA WINDS.
      KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
                 NOTE: IF KGDSO(201)=1 THEN THE 9TH BIT OF KGDSO(11)
                       IS TEMPORARILY SET TO 1 TO ALERT THE GDS WIZARD
                       THAT THESE FIELDS ARE STAGGERED ETA WINDS.
      MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF INPUT GRID FIELDS IF KM=1
      MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
      KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
      IBI      - INTEGER (KM) INPUT BITMAP FLAGS
      LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF RESPECTIVE IBI(K)=1)
      UI       - REAL (MI,KM) INPUT U-COMPONENT FIELDS TO INTERPOLATE
      VI       - REAL (MI,KM) INPUT V-COMPONENT FIELDS TO INTERPOLATE
      NO       - INTEGER NUMBER OF OUTPUT POINTS (IF KGDSO(1)<0)
      RLAT     - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
      RLON     - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
      CROT     - REAL (NO) VECTOR ROTATION COSINES (IF KGDSO(1)<0)
      SROT     - REAL (NO) VECTOR ROTATION SINES (IF KGDSO(1)<0)
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
 
   Output argument list:
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
      RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
      RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
      CROT     - REAL (MO) VECTOR ROTATION COSINES (IF KGDSO(1)>=0)
      SROT     - REAL (MO) VECTOR ROTATION SINES (IF KGDSO(1)>=0)
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
      IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
      LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
      UO       - REAL (MO,KM) OUTPUT U-COMPONENT FIELDS INTERPOLATED
      VO       - REAL (MO,KM) OUTPUT V-COMPONENT FIELDS INTERPOLATED
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL INTERPOLATION
                 1    UNRECOGNIZED INTERPOLATION METHOD
                 2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
                 3    UNRECOGNIZED OUTPUT GRID
                 1X   INVALID BICUBIC METHOD PARAMETERS
                 3X   INVALID BUDGET METHOD PARAMETERS
                 4X   INVALID SPECTRAL METHOD PARAMETERS
 
REMARKS:
    EXAMPLES DEMONSTRATING RELATIVE CPU COSTS.
    THIS EXAMPLE IS INTERPOLATING 12 LEVELS OF WINDS
    FROM THE 360 X 181 GLOBAL GRID (NCEP GRID 3)
    TO THE 93 X 68 HAWAIIAN MERCATOR GRID (NCEP GRID 204).
    THE EXAMPLE TIMES ARE FOR THE C90.  AS A REFERENCE, THE CP TIME
    FOR UNPACKING THE GLOBAL 12 PAIRS OF WIND FIELDS IS 0.07 SECONDS.
 
    BILINEAR    0                   0.05
    BICUBIC     1   0               0.16
    BICUBIC     1   1               0.17
    NEIGHBOR    2                   0.02
    BUDGET      3   -1,-1           0.94
    SPECTRAL    4   0,40            0.31
    SPECTRAL    4   1,40            0.33
    SPECTRAL    4   0,-1            0.59
    N-BUDGET    6   0,-1            0.31
 
    THE SPECTRAL INTERPOLATION IS FAST FOR THE MERCATOR GRID.
    HOWEVER, FOR SOME GRIDS THE SPECTRAL INTERPOLATION IS SLOW.
    THE FOLLOWING EXAMPLE IS INTERPOLATING 12 LEVELS OF WINDS
    FROM THE 360 X 181 GLOBAL GRID (NCEP GRID 3)
    TO THE 93 X 65 CONUS LAMBERT CONFORMAL GRID (NCEP GRID 211).
 

IPSECTOR

The IPSECTOR routine creates a subsector of a grid,
            computing a new grid definition section and copying data
            from the original grid data.

USAGE:    CALL IPSECTOR(I1,I2,J1,J2,NF,M,KGDS,L,F,MS,
                        KGDSS,LS,FS,IRET)

   Input argument list:
      I1       - INTEGER FIRST X POINT OF THE SECTOR
                 OR IF 0<=I2IPSPASTE

The IPSPASTE routine pastes a subsector of a grid back
             into the original grid.

USAGE:    CALL IPSPASTE(I1,I2,J1,J2,NF,MS,LS,FS,M,KGDS,L,F,IRET)

   Input argument list:
      I1       - INTEGER FIRST X POINT OF THE SECTOR
                 OR IF 0<=I2IPXETAS

The IPXETAS routine transforms between the staggered eta grids
            as used in the eta model and for native grid transmission
            and their full expansion as used for general interpolation
            and graphics.  The eta grids are rotated latitude-longitude
            grids staggered as defined by the Arakawa e-grid, that is
            with mass data points alternating with wind data points.
            The expansion of the fields is done by 4-point averaging.

USAGE:    CALL IPXETAS(IDIR,M1,M2,KM,KGDS1,F1,KGDS2,F2,IRET)

   Input argument list:
      IDIR     - INTEGER TRANSFORM OPTION
                 (+1 TO EXPAND STAGGERED MASS FIELDS TO FULL FIELDS)
                 (+2 TO EXPAND STAGGERED WIND FIELDS TO FULL FIELDS)
                 (-1 TO CONTRACT FULL MASS FIELDS TO STAGGERED FIELDS)
                 (-2 TO CONTRACT FULL WIND FIELDS TO STAGGERED FIELDS)
      M1       - INTEGER SKIP NUMBER BETWEEN STAGGERED GRID FIELDS
      M2       - INTEGER SKIP NUMBER BETWEEN FULL GRID FIELDS
      KM       - INTEGER NUMBER OF FIELDS TO TRANSFORM
      KGDS1    - INTEGER (200) GDS PARMS OF STAGGERED GRID IF IDIR>0
      F1       - REAL (M1,KM) STAGGERED GRID FIELDS IF IDIR>0
      KGDS2    - INTEGER (200) GDS PARMS OF FULL GRID IF IDIR<0
      F2       - REAL (M2,KM) FULL GRID FIELDS IF IDIR<0

   Output argument list:
      KGDS1    - INTEGER (200) GDS PARMS OF STAGGERED GRID IF IDIR<0
      F1       - REAL (M1,KM) STAGGERED GRID FIELDS IF IDIR<0
      KGDS2    - INTEGER (200) GDS PARMS OF FULL GRID IF IDIR>0
      F2       - REAL (M2,KM) FULL GRID FIELDS IF IDIR>0
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL TRANSFORMATION
                 1    IMPROPER GRID SPECIFICATION


IPXWAFS

The IPXWAFS routine transforms between the thinned WAFS grids
            as used for transmitting to the aviation community
            and their full expansion as used for general interpolation
            and graphics.  The WAFS grids are latitude-longitude grids
            thinned only in the zonal direction as indicated
            By the PL parameters in the grib grid description section.
            The PL parameters must be supplied for contraction
            (starting at KGDS1(22)).  Otherwise (If KGDS1(22)<=0)
            grid contraction is performed specifically for
            The 1.25 degree NCEP WAFS grid identifications 37-44.
            The expansion and contraction of the fields are done
            by linear interpolation, so that they are not reversible.

USAGE:    CALL IPXWAFS(IDIR,M1,M2,KM,KGDS1,F1,KGDS2,F2,IRET)

   Input argument list:
      IDIR     - INTEGER TRANSFORM OPTION
                 (+1 TO EXPAND THINNED FIELDS TO FULL FIELDS)
                 (-1 TO CONTRACT FULL FIELDS TO THINNED FIELDS)
      M1       - INTEGER SKIP NUMBER BETWEEN THINNED GRID FIELDS
      M2       - INTEGER SKIP NUMBER BETWEEN FULL GRID FIELDS
      KM       - INTEGER NUMBER OF FIELDS TO TRANSFORM
      KGDS1    - INTEGER (200) GDS PARMS OF THINNED GRID IF IDIR>0
                 (IF IDIR<0, THEN EITHER THE PL PARAMETERS STARTING AT
                  KGDS1(22) MUST BE SUPPLIED OR IF KGDS1(22)<=0,
                  THEN THE PL PARAMETERS DEFAULT TO THOSE FOR
                  SPECIFIC NCEP WAFS GRIDS 37-44).
      F1       - REAL (M1,KM) THINNED GRID FIELDS IF IDIR>0
      KGDS2    - INTEGER (200) GDS PARMS OF FULL GRID IF IDIR<0
      F2       - REAL (M2,KM) FULL GRID FIELDS IF IDIR<0
 
   Output argument list:
      KGDS1    - INTEGER (200) GDS PARMS OF THINNED GRID IF IDIR<0
      F1       - REAL (M1,KM) THINNED GRID FIELDS IF IDIR<0
      KGDS2    - INTEGER (200) GDS PARMS OF FULL GRID IF IDIR>0
      F2       - REAL (M2,KM) FULL GRID FIELDS IF IDIR>0
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL TRANSFORMATION
                 1    IMPROPER GRID SPECIFICATION


IPXWAFS2

The IPXWAFS2 routine transforms between the thinned wafs grids
             as used for transmitting to the aviation community
             and their full expansion as used for general interpolation
             and graphics.  The WAFS grids are latitude-longitude grids
             thinned only in the zonal direction as indicated
             by the PL parameters in the grib grid description section.
             The PL parameters must be supplied for contraction
             (starting at KGDS1(22)).  Otherwise (if KGDS1(22)<=0)
             grid contraction is performed specifically for
             The 1.25 degree NCEP WAFS grid identifications 37-44.
             The expansion and contraction of the fields are done
             by linear interpolation, so that they are not reversible.
             This version allows a bitmap.

USAGE:    CALL IPXWAFS2(IDIR,M1,M2,KM,
                        KGDS1,IB1,L1,F1,KGDS2,IB2,L2,F2,IRET)

   Input argument list:
      IDIR     - INTEGER TRANSFORM OPTION
                 (+1 TO EXPAND THINNED FIELDS TO FULL FIELDS)
                 (-1 TO CONTRACT FULL FIELDS TO THINNED FIELDS)
      M1       - INTEGER SKIP NUMBER BETWEEN THINNED GRID FIELDS
      M2       - INTEGER SKIP NUMBER BETWEEN FULL GRID FIELDS
      KM       - INTEGER NUMBER OF FIELDS TO TRANSFORM
      KGDS1    - INTEGER (200) GDS PARMS OF THINNED GRID IF IDIR>0
                 (IF IDIR<0, THEN EITHER THE PL PARAMETERS STARTING AT
                  KGDS1(22) MUST BE SUPPLIED OR IF KGDS1(22)<=0,
                  THEN THE PL PARAMETERS DEFAULT TO THOSE FOR
                  SPECIFIC NCEP WAFS GRIDS 37-44).
      IB1      - INTEGER (KM) THINNED BITMAP FLAGS IF IDIR>0
      L1       - LOGICAL*1 (M1,KM) THINNED BITMAP FIELDS IF IDIR>0
      F1       - REAL (M1,KM) THINNED GRID FIELDS IF IDIR>0
      KGDS2    - INTEGER (200) GDS PARMS OF FULL GRID IF IDIR<0
      IB2      - INTEGER (KM) FULL BITMAP FLAGS IF IDIR<0
      L2       - LOGICAL*1 (M1,KM) FULL BITMAP FIELDS IF IDIR<0
      F2       - REAL (M2,KM) FULL GRID FIELDS IF IDIR<0
 
   Output argument list:
      KGDS1    - INTEGER (200) GDS PARMS OF THINNED GRID IF IDIR<0
      IB1      - INTEGER (KM) THINNED BITMAP FLAGS IF IDIR<0
      L1       - LOGICAL*1 (M1,KM) THINNED BITMAP FIELDS IF IDIR<0
      F1       - REAL (M1,KM) THINNED GRID FIELDS IF IDIR<0
      KGDS2    - INTEGER (200) GDS PARMS OF FULL GRID IF IDIR>0
      IB2      - INTEGER (KM) FULL BITMAP FLAGS IF IDIR>0
      L2       - LOGICAL*1 (M1,KM) FULL BITMAP FIELDS IF IDIR>0
      F2       - REAL (M2,KM) FULL GRID FIELDS IF IDIR>0
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL TRANSFORMATION
                 1    IMPROPER GRID SPECIFICATION
 
MAKGDS

The MAKGDS routine makes or breaks a grid description section.
           It can do one of the following:
              (IOPT=-1)    unpack a GDS into w3fi63 kgds integer form
              (IOPT=255)   pack a GDS from W3FI63 KGDS integer form
              (IOPT<255) pack a GDS from an NCEP grid identification.

USAGE:    CALL MAKGDS(IOPT,KGDS,GDS,LENGDS,IRET)

   Input argument list:
      IOPT     - INTEGER OPTION
                 IOPT=-1 TO UNPACK GDS INTO KGDS;
                 IOPT=255 TO USE KGDS TO PACK GDS;
                 IOPT<255 NCEP GRID ID TO MAKE GDS AND KGDS.
      KGDS     - INTEGER (200) W3FI63-STYLE UNPACKED GDS (IF IOPT=255)
                 (ONLY FIRST 22 VALUES ARE ACCESSED IF KGDS(20)=255.)
                 (SEE REMARKS BELOW FOR A DETAILED DESCRIPTION OF KGDS.)
      GDS      - CHARACTER (400) GRID DEFINITION SECTION (IF IOPT=-1)
 
   Output argument list:
      KGDS     - INTEGER (200) W3FI63-STYLE UNPACKED GDS (IF IOPT<255)
                 (ONLY FIRST 22 VALUES ARE ACCESSED IF KGDS(20)=255.)
                 (SEE REMARKS BELOW FOR A DETAILED DESCRIPTION OF KGDS.)
      GDS      - CHARACTER (400) GRID DEFINITION SECTION (IF IOPT>0)
      LENGDS   - INTEGER LENGTH OF THE GDS (IF IOPT>0)
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL
                 1    GRID REPRESENTATION TYPE NOT VALID
                 4    DATA REPRESENTATION TYPE NOT CURRENTLY ACCEPTABLE
 
REMARKS:
        THE KGDS PARAMETERS ARE DESCRIBED BELOW
           AS COPIED FROM THE W3FI63 DOCBLOCK.
           (1)   - DATA REPRESENTATION TYPE
           (19)  - NUMBER OF VERTICAL COORDINATE PARAMETERS
           (20)  - OCTET NUMBER OF THE LIST OF VERTICAL COORDINATE
                   PARAMETERS
                   OR
                   OCTET NUMBER OF THE LIST OF NUMBERS OF POINTS
                   IN EACH ROW
                   OR
                   255 IF NEITHER ARE PRESENT
           (21)  - FOR GRIDS WITH PL, NUMBER OF POINTS IN GRID
           (22)  - NUMBER OF WORDS IN EACH ROW
        LATITUDE/LONGITUDE GRIDS
           (2)   - N(I) NR POINTS ON LATITUDE CIRCLE
           (3)   - N(J) NR POINTS ON LONGITUDE MERIDIAN
           (4)   - LA(1) LATITUDE OF ORIGIN
           (5)   - LO(1) LONGITUDE OF ORIGIN
           (6)   - RESOLUTION FLAG (RIGHT ADJ COPY OF OCTET 17)
           (7)   - LA(2) LATITUDE OF EXTREME POINT
           (8)   - LO(2) LONGITUDE OF EXTREME POINT
           (9)   - DI LATITUDINAL DIRECTION OF INCREMENT
           (10)  - DJ LONGITUDINAL DIRECTION INCREMENT
           (11)  - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28)
        GAUSSIAN  GRIDS
           (2)   - N(I) NR POINTS ON LATITUDE CIRCLE
           (3)   - N(J) NR POINTS ON LONGITUDE MERIDIAN
           (4)   - LA(1) LATITUDE OF ORIGIN
           (5)   - LO(1) LONGITUDE OF ORIGIN
           (6)   - RESOLUTION FLAG  (RIGHT ADJ COPY OF OCTET 17)
           (7)   - LA(2) LATITUDE OF EXTREME POINT
           (8)   - LO(2) LONGITUDE OF EXTREME POINT
           (9)   - DI LATITUDINAL DIRECTION OF INCREMENT
           (10)  - N - NR OF CIRCLES POLE TO EQUATOR
           (11)  - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28)
           (12)  - NV - NR OF VERT COORD PARAMETERS
           (13)  - PV - OCTET NR OF LIST OF VERT COORD PARAMETERS
                              OR
                   PL - LOCATION OF THE LIST OF NUMBERS OF POINTS IN
                        EACH ROW (IF NO VERT COORD PARAMETERS
                        ARE PRESENT
                              OR
                   255 IF NEITHER ARE PRESENT
        POLAR STEREOGRAPHIC GRIDS
           (2)   - N(I) NR POINTS ALONG LAT CIRCLE
           (3)   - N(J) NR POINTS ALONG LON CIRCLE
           (4)   - LA(1) LATITUDE OF ORIGIN
           (5)   - LO(1) LONGITUDE OF ORIGIN
           (6)   - RESOLUTION FLAG  (RIGHT ADJ COPY OF OCTET 17)
           (7)   - LOV GRID ORIENTATION
           (8)   - DX - X DIRECTION INCREMENT
           (9)   - DY - Y DIRECTION INCREMENT
           (10)  - PROJECTION CENTER FLAG
           (11)  - SCANNING MODE (RIGHT ADJ COPY OF OCTET 28)
        SPHERICAL HARMONIC COEFFICIENTS
           (2)   - J PENTAGONAL RESOLUTION PARAMETER
           (3)   - K      "          "         "
           (4)   - M      "          "         "
           (5)   - REPRESENTATION TYPE
           (6)   - COEFFICIENT STORAGE MODE
        MERCATOR GRIDS
           (2)   - N(I) NR POINTS ON LATITUDE CIRCLE
           (3)   - N(J) NR POINTS ON LONGITUDE MERIDIAN
           (4)   - LA(1) LATITUDE OF ORIGIN
           (5)   - LO(1) LONGITUDE OF ORIGIN
           (6)   - RESOLUTION FLAG (RIGHT ADJ COPY OF OCTET 17)
           (7)   - LA(2) LATITUDE OF LAST GRID POINT
           (8)   - LO(2) LONGITUDE OF LAST GRID POINT
           (9)   - LATIT - LATITUDE OF PROJECTION INTERSECTION
           (10)  - RESERVED
           (11)  - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28)
           (12)  - LONGITUDINAL DIR GRID LENGTH
           (13)  - LATITUDINAL DIR GRID LENGTH
        LAMBERT CONFORMAL GRIDS
           (2)   - NX NR POINTS ALONG X-AXIS
           (3)   - NY NR POINTS ALONG Y-AXIS
           (4)   - LA1 LAT OF ORIGIN (LOWER LEFT)
           (5)   - LO1 LON OF ORIGIN (LOWER LEFT)
           (6)   - RESOLUTION (RIGHT ADJ COPY OF OCTET 17)
           (7)   - LOV - ORIENTATION OF GRID
           (8)   - DX - X-DIR INCREMENT
           (9)   - DY - Y-DIR INCREMENT
           (10)  - PROJECTION CENTER FLAG
           (11)  - SCANNING MODE FLAG (RIGHT ADJ COPY OF OCTET 28)
           (12)  - LATIN 1 - FIRST LAT FROM POLE OF SECANT CONE INTER
           (13)  - LATIN 2 - SECOND LAT FROM POLE OF SECANT CONE INTER
 

MOVECT

The MOVECT routine provides the rotation parameters
           to move a vector along a great circle from one
           position to another while conserving its orientation
           with respect to the great circle.  These rotation
           parameters are useful for vector interpolation.

USAGE:    CALL MOVECT(FLAT,FLON,TLAT,TLON,CROT,SROT)

   Input argument list:
      FLAT     - REAL LATITUDE IN DEGREES FROM WHICH TO MOVE THE VECTOR
      FLON     - REAL LONGITUDE IN DEGREES FROM WHICH TO MOVE THE VECTOR
      TLAT     - REAL LATITUDE IN DEGREES TO WHICH TO MOVE THE VECTOR
      TLON     - REAL LONGITUDE IN DEGREES TO WHICH TO MOVE THE VECTOR

   Output argument list:
      CROT     - REAL CLOCKWISE VECTOR ROTATION COSINE
      SROT     - REAL CLOCKWISE VECTOR ROTATION SINE
                 (UTO=CROT*UFROM-SROT*VFROM;
                  VTO=SROT*UFROM+CROT*VFROM)

REMARKS:
        THIS SUBPROGRAM IS CORRECT TO SEVEN DIGITS ON THE CRAYS.
        USE DOUBLE PRECISION IF BETTER PRECISION IS REQUIRED.


POLATEG0

The POLATEG0 routine performs bilinear interpolation
             from any grid to any grid for scalar fields,
             returning their vector gradients.  No options are
             allowed.  Only horizontal interpolation is performed.
             The grids are defined by their grid description sections
             (passed in integer form as decoded by subprogram w3fi63).

             The current code recognizes the following projections:
                (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
                (KGDS(1)=001) MERCATOR CYLINDRICAL
                (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
                (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
                (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
                (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)

             Where KGDS could be either input KGDSI or output KGDSO.
             As an added bonus the number of output grid points
             and their latitudes and longitudes are also returned.
             On the other hand, the output can be a set of station points
             If KGDSO(1)<0, in which case the number of points
             and their latitudes and longitudes must be input.
             Input bitmaps will be interpolated to output bitmaps.
             Output bitmaps will also be created when the output grid
             extends outside of the domain of the input grid.
             The output field is set to 0 where the output bitmap is off.
 

USAGE:    CALL POLATEG0(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
                        NO,RLAT,RLON,CROT,SROT,IBO,LO,XO,YO,IRET)

   INPUT ARGUMENT LIST:
      IPOPT    - INTEGER (20) INTERPOLATION OPTIONS (NO OPTIONS)
      KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
      KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
                 (KGDSO(1)<0 IMPLIES RANDOM STATION POINTS)
      MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF INPUT GRID FIELDS IF KM=1
      MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
      KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
      IBI      - INTEGER (KM) INPUT BITMAP FLAGS
      LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
      GI       - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)
      RLAT     - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
      RLON     - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
      CROT     - REAL (NO) VECTOR ROTATION COSINES (IF KGDSO(1)<0)
      SROT     - REAL (NO) VECTOR ROTATION SINES (IF KGDSO(1)<0)
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
 
   OUTPUT ARGUMENT LIST:
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
      RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
      RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
      IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
      LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
      XO       - REAL (MO,KM) OUTPUT X-GRADIENT FIELDS INTERPOLATED
      YO       - REAL (MO,KM) OUTPUT Y-GRADIENT FIELDS INTERPOLATED
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL INTERPOLATION
                 2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
                 3    UNRECOGNIZED OUTPUT GRID

REMARKS:
       THE GRADIENT COMPUTATIONS ARE NOT ROBUST NEAR THE POLES.
       IN FACT, NO GRADIENTS ARE COMPUTED POLEWARD OF 89 LATITUDE.


POLATEG1

The POLATEG1 routine performs bicubic interpolation
             from any grid to any grid for scalar fields,
             returning their vector gradients. Bitmaps are now
             allowed, but bilinear interpolation is done when any
             invalid data is within the bicubic template.
             options allow choices between Straight Bicubic (ipopt(1)=0)
             and Constrained Bicubic (ipopt(1)=1) where the value is
             confined within the range of the surrounding 4 points.
             bilinear used within one grid length of boundaries.
             only horizontal interpolation is performed.
             The grids are defined by their grid description sections
             (passed in integer form as decoded by subprogram w3fi63).

             The current code recognizes the following projections:
                (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
                (KGDS(1)=001) MERCATOR CYLINDRICAL
                (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
                (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
                (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
                (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)

             Where KGDS could be either input KGDSI or output KGDSO.
             As an added bonus the number of output grid points
             and their latitudes and longitudes are also returned.
             On the other hand, the output can be a set of station points
             If KGDSO(1)<0, in which case the number of points
             and their latitudes and longitudes must be input.
             Input bitmaps will be interpolated to output bitmaps.
             Output bitmaps will also be created when the output grid
             extends outside of the domain of the input grid.
             The output field is set to 0 where the output bitmap is off.

USAGE:    CALL POLATEG1(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
                        NO,RLAT,RLON,CROT,SROT,IBO,LO,XO,YO,IRET)

   Input argument list:
      IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
                 IPOPT(1)=0 FOR STRAIGHT BICUBIC;
                 IPOPT(1)=1 FOR CONSTRAINED BICUBIC WHERE VALUE IS
                 CONFINED WITHIN THE RANGE OF THE SURROUNDING 4 POINTS.
      KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
      KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
                 (KGDSO(1)<0 IMPLIES RANDOM STATION POINTS)
      MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF INPUT GRID FIELDS IF KM=1
      MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
      KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
      IBI      - INTEGER (KM) INPUT BITMAP FLAGS (MUST BE ALL 0)
      LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
      GI       - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)
      RLAT     - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
      RLON     - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
      CROT     - REAL (NO) VECTOR ROTATION COSINES (IF KGDSO(1)<0)
      SROT     - REAL (NO) VECTOR ROTATION SINES (IF KGDSO(1)<0)
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
 
    Output argument list:
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
      RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
      RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
      IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
      LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
      XO       - REAL (MO,KM) OUTPUT X-GRADIENT FIELDS INTERPOLATED
      YO       - REAL (MO,KM) OUTPUT Y-GRADIENT FIELDS INTERPOLATED
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL INTERPOLATION
                 2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
                 3    UNRECOGNIZED OUTPUT GRID


POLATEG4

The POLATEG4 routine PERFORMS SPECTRAL INTERPOLATION
             from any grid to any grid for scalar fields,
             returning their vector gradients.  It requires
             that the input fields be uniformly global.
             Options allow choices between triangular shape (IPOPT(1)=0)
             and rhomboidal shape (IPOPT(1)=1) which has no default;
             A second option is the truncation (IPOPT(2)) which defaults
             To a sensible truncation for the input grid (if OPT(2)=-1).
             Note that if the output grid is not found in a special list,
             then the transform back to grid is not very fast.
             This special list contains global cylindrical grids,
             Polar Stereographic grids centered at the pole and Mercator
             grids.  Only horizontal interpolation is performed.
             The grids are defined by their grid description sections
             (passed in integer form as decoded by subprogram w3fi63).

             The current code recognizes the following projections:
                (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
                (KGDS(1)=001) MERCATOR CYLINDRICAL
                (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
                (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
                (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
                (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)

             Where KGDS could be either input KGDSI or output KGDSO.
             As an added bonus the number of output grid points
             and their latitudes and longitudes are also returned.
             On the other hand, the output can be a set of station points
             If KGDSO(1)<0, in which case the number of points
             and their latitudes and longitudes must be input.
             Output bitmaps will not be created.

USAGE:    CALL POLATEG4(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
                        NO,RLAT,RLON,CROT,SROT,IBO,LO,GO,IRET)

   Input argument list:
      IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
                 IPOPT(1)=0 FOR TRIANGULAR, IPOPT(1)=1 FOR RHOMBOIDAL;
                 IPOPT(2) IS TRUNCATION NUMBER
                 (DEFAULTS TO SENSIBLE IF IPOPT(2)=-1).
      KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
      KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
                 (KGDSO(1)<0 IMPLIES RANDOM STATION POINTS)
      MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF INPUT GRID FIELDS IF KM=1
      MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
      KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
      IBI      - INTEGER (KM) INPUT BITMAP FLAGS (MUST BE ALL 0)
      LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
      GI       - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)
      RLAT     - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
      RLON     - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
      CROT     - REAL (NO) VECTOR ROTATION COSINES (IF KGDSO(1)<0)
      SROT     - REAL (NO) VECTOR ROTATION SINES (IF KGDSO(1)<0)
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
 
   Output argument list:
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
      RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
      RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
      IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
      LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
      GO       - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL INTERPOLATION
                 2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
                 3    UNRECOGNIZED OUTPUT GRID
                 41   INVALID NONGLOBAL INPUT GRID
                 42   INVALID SPECTRAL METHOD PARAMETERS
                 43   UNIMPLEMENTED OUTPUT GRID
 

POLATES0

The POLATES0 routine performs bilinear interpolation
             from any grid to any grid for scalar fields.
             options allow varying the minimum percentage for mask,
             i.e. percent valid input data required to make output data,
             (ipopt(1)) which defaults to 50 (if ipopt(1)=-1).
             Only horizontal interpolation is performed.
             The grids are defined by their grid description sections
             (passed in integer form as decoded by subprogram w3fi63).

             The current code recognizes the following projections:
                (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
                (KGDS(1)=001) MERCATOR CYLINDRICAL
                (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
                (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
                (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
                (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)

             Where KGDS could be either input KGDSI or output KGDSO.
             As an added bonus the number of output grid points
             and their latitudes and longitudes are also returned.
             On the other hand, the output can be a set of station points
             If KGDSO(1)<0, in which case the number of points
             and their latitudes and longitudes must be input.
             Input bitmaps will be interpolated to output bitmaps.
             Output bitmaps will also be created when the output grid
             extends outside of the domain of the input grid.
             The output field is set to 0 where the output bitmap is off.

USAGE:    CALL POLATES0(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
                        NO,RLAT,RLON,IBO,LO,GO,IRET)

   Input argument list:
      IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
                 IPOPT(1) IS MINIMUM PERCENTAGE FOR MASK
                 (DEFAULTS TO 50 IF IPOPT(1)=-1)
      KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
      KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
                 (KGDSO(1)<0 IMPLIES RANDOM STATION POINTS)
      MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF INPUT GRID FIELDS IF KM=1
      MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
      KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
      IBI      - INTEGER (KM) INPUT BITMAP FLAGS
      LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
      GI       - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)
      RLAT     - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
      RLON     - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
 
   Output argument list:
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
      RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
      RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
      IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
      LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
      GO       - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL INTERPOLATION
                 2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
                 3    UNRECOGNIZED OUTPUT GRID
 

POLATES1

The POLATES1 routine performs bicubic interpolation
             from any grid to any grid for scalar fields.
             bitmaps are now allowed, but bilinear interpolation is done
             when any invalid data is within the bicubic template.
             Options allow choices between straight bicubic (ipopt(1)=0)
             and constrained bicubic (ipopt(1)=1) where the value is
             confined within the range of the surrounding 4 points.
             Another option is the minimum percentage for mask,
             i.e. percent valid input data required to make output data,
             (ipopt(2)) which defaults to 50 (if ipopt(2)=-1).
             Bilinear used within one grid length of boundaries.
             Only horizontal interpolation is performed.
             The grids are defined by their grid description sections
             (passed in integer form as decoded by subprogram w3fi63).

             The current code recognizes the following projections:
                (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
                (KGDS(1)=001) MERCATOR CYLINDRICAL
                (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
                (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
                (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
                (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)

             Where KGDS could be either input KGDSI or output KGDSO.
             As an added bonus the number of output grid points
             and their latitudes and longitudes are also returned.
             On the other hand, the output can be a set of station points
             If KGDSO(1)<0, in which case the number of points
             and their latitudes and longitudes must be input.
             Input bitmaps will be interpolated to output bitmaps.
             Output bitmaps will also be created when the output grid
             extends outside of the domain of the input grid.
             The output field is set to 0 where the output bitmap is off.

USAGE:    CALL POLATES1(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
                        NO,RLAT,RLON,IBO,LO,GO,IRET)

   Input argument list:
      IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
                 IPOPT(1)=0 FOR STRAIGHT BICUBIC;
                 IPOPT(1)=1 FOR CONSTRAINED BICUBIC WHERE VALUE IS
                 CONFINED WITHIN THE RANGE OF THE SURROUNDING 4 POINTS.
                 IPOPT(2) IS MINIMUM PERCENTAGE FOR MASK
                 (DEFAULTS TO 50 IF IPOPT(2)=-1)
      KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
      KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
                 (KGDSO(1)<0 IMPLIES RANDOM STATION POINTS)
      MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF INPUT GRID FIELDS IF KM=1
      MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
      KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
      IBI      - INTEGER (KM) INPUT BITMAP FLAGS
      LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
      GI       - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)
      RLAT     - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
      RLON     - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
 
   Output argument list:
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
      RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
      RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
      IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
      LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
      GO       - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL INTERPOLATION
                 2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
                 3    UNRECOGNIZED OUTPUT GRID
 

POLATES2

The POLATES2 routine performs neighbor interpolation
             from any grid to any grid for scalar fields.
             options allow choosing the width of the grid square
             (ipopt(1)) to search for valid data, which defaults to 1
             (if ipopt(1)=-1).  Odd width squares are centered on
             the nearest input grid point; Even width squares are
             centered on the nearest four input grid points.
             Squares are searched for valid data in a spiral pattern
             starting from the center.  No searching is done where
             the output grid is outside the input grid.
             Only horizontal interpolation is performed.
             The grids are defined by their grid description sections
             (passed in integer form as decoded by subprogram w3fi63).

             The current code recognizes the following projections:
                (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
                (KGDS(1)=001) MERCATOR CYLINDRICAL
                (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
                (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
                (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
                (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)

             Where KGDS could be either input KGDSI or output KGDSO.
             As an added bonus the number of output grid points
             and their latitudes and longitudes are also returned.
             On the other hand, the output can be a set of station points
             If KGDSO(1)<0, in which case the number of points
             and their latitudes and longitudes must be input.
             Input bitmaps will be interpolated to output bitmaps.
             Output bitmaps will also be created when the output grid
             extends outside of the domain of the input grid.
             The output field is set to 0 where the output bitmap is off.

USAGE:    CALL POLATES2(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
                        NO,RLAT,RLON,IBO,LO,GO,IRET)

   Input argument list:
      IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
                 IPOPT(1) IS WIDTH OF SQUARE TO EXAMINE IN SPIRAL SEARCH
                 (DEFAULTS TO 1 IF IPOPT(1)=-1)
      KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
      KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
                 (KGDSO(1)<0 IMPLIES RANDOM STATION POINTS)
      MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF INPUT GRID FIELDS IF KM=1
      MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
      KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
      IBI      - INTEGER (KM) INPUT BITMAP FLAGS
      LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
      GI       - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)
      RLAT     - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
      RLON     - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
 
   Output argument list:
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
      RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
      RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
      IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
      LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
      GO       - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL INTERPOLATION
                 2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
                 3    UNRECOGNIZED OUTPUT GRID
 

POLATES3

The POLATES3 routine PERFORMS BUDGET INTERPOLATION
             FROM ANY GRID TO ANY GRID FOR SCALAR FIELDS.
             IT REQUIRES A GRID FOR THE OUTPUT FIELDS (KGDSO(1)>=0).
             THE ALGORITHM SIMPLY COMPUTES (WEIGHTED) AVERAGES
             OF BILINEARLY INTERPOLATED POINTS ARRANGED IN A SQUARE BOX
             CENTERED AROUND EACH OUTPUT GRID POINT AND STRETCHING
             NEARLY HALFWAY TO EACH OF THE NEIGHBORING GRID POINTS.
             OPTIONS ALLOW CHOICES OF NUMBER OF POINTS IN EACH RADIUS
             FROM THE CENTER POINT (IPOPT(1)) WHICH DEFAULTS TO 2
             (IF IPOPT(1)=-1) MEANING THAT 25 POINTS WILL BE AVERAGED;
             FURTHER OPTIONS ARE THE RESPECTIVE WEIGHTS FOR THE RADIUS
             POINTS STARTING AT THE CENTER POINT (IPOPT(2:2+IPOPT(1))
             WHICH DEFAULTS TO ALL 1 (IF IPOPT(1)=-1 OR IPOPT(2)=-1).
             A SPECIAL INTERPOLATION IS DONE IF IPOPT(2)=-2.
             IN THIS CASE, THE BOXES STRETCH NEARLY ALL THE WAY TO
             EACH OF THE NEIGHBORING GRID POINTS AND THE WEIGHTS
             ARE THE ADJOINT OF THE BILINEAR INTERPOLATION WEIGHTS.
             THIS CASE GIVES QUASI-SECOND-ORDER BUDGET INTERPOLATION.
             ANOTHER OPTION IS THE MINIMUM PERCENTAGE FOR MASK,
             I.E. PERCENT VALID INPUT DATA REQUIRED TO MAKE OUTPUT DATA,
             (IPOPT(3+IPOPT(1)) WHICH DEFAULTS TO 50 (IF -1).
             ONLY HORIZONTAL INTERPOLATION IS PERFORMED.
             The grids are defined by their grid description sections
             (passed in integer form as decoded by subprogram w3fi63).

             The current code recognizes the following projections:
                (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
                (KGDS(1)=001) MERCATOR CYLINDRICAL
                (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
                (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
                (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
                (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)

             Where KGDS could be either input KGDSI or output KGDSO.
             As an added bonus the number of output grid points
             and their latitudes and longitudes are also returned.
             On the other hand, the output can be a set of station points
             If KGDSO(1)<0, in which case the number of points
             and their latitudes and longitudes must be input.
             Input bitmaps will be interpolated to output bitmaps.
             Output bitmaps will also be created when the output grid
             extends outside of the domain of the input grid.
             The output field is set to 0 where the output bitmap is off.


USAGE:    CALL POLATES3(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
                        NO,RLAT,RLON,IBO,LO,GO,IRET)

   Input argument list:
      IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
                 IPOPT(1) IS NUMBER OF RADIUS POINTS
                 (DEFAULTS TO 2 IF IPOPT(1)=-1);
                 IPOPT(2:2+IPOPT(1)) ARE RESPECTIVE WEIGHTS
                 (DEFAULTS TO ALL 1 IF IPOPT(1)=-1 OR IPOPT(2)=-1).
                 IPOPT(3+IPOPT(1)) IS MINIMUM PERCENTAGE FOR MASK
                 (DEFAULTS TO 50 IF IPOPT(3+IPOPT(1)=-1)
      KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
      KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
      MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF INPUT GRID FIELDS IF KM=1
      MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
      KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
      IBI      - INTEGER (KM) INPUT BITMAP FLAGS
      LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
      GI       - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
 
  Output argument list:
      NO       - INTEGER NUMBER OF OUTPUT POINTS
      RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES
      RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES
      IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
      LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
      GO       - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL INTERPOLATION
                 2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
                 3    UNRECOGNIZED OUTPUT GRID
                 31   INVALID UNDEFINED OUTPUT GRID
                 32   INVALID BUDGET METHOD PARAMETERS
 

POLATES4

The POLATES4 routine performs spectral interpolation
             from any grid to any grid for scalar fields.
             It requires that the input fields be uniformly global.
             Options allow choices between triangular shape (ipopt(1)=0)
             and rhomboidal shape (ipopt(1)=1) which has no default;
             A second option is the truncation (ipopt(2)) which defaults
             to a sensible truncation for the input grid (if opt(2)=-1).
             Note that if the output grid is not found in a special list,
             then the transform back to grid is not very fast.
             This special list contains global cylindrical grids,
             Polar Stereographic grids centered at the pole and
             Mercator grids.

             The current code recognizes the following projections:
                (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
                (KGDS(1)=001) MERCATOR CYLINDRICAL
                (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
                (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
                (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
                (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)

             Where KGDS could be either input KGDSI or output KGDSO.
             As an added bonus the number of output grid points
             and their latitudes and longitudes are also returned.
             On the other hand, the output can be a set of station points
             If KGDSO(1)<0, in which case the number of points
             and their latitudes and longitudes must be input.
             Output bitmaps will not be created.

USAGE:    CALL POLATES4(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
                        NO,RLAT,RLON,IBO,LO,GO,IRET)

   Input argument list:
      IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
                 IPOPT(1)=0 FOR TRIANGULAR, IPOPT(1)=1 FOR RHOMBOIDAL;
                 IPOPT(2) IS TRUNCATION NUMBER
                 (DEFAULTS TO SENSIBLE IF IPOPT(2)=-1).
      KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
      KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
                 (KGDSO(1)<0 IMPLIES RANDOM STATION POINTS)
      MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF INPUT GRID FIELDS IF KM=1
      MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
      KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
      IBI      - INTEGER (KM) INPUT BITMAP FLAGS (MUST BE ALL 0)
      LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
      GI       - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)
      RLAT     - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
      RLON     - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
 
   Output argument list:
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
      RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
      RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
      IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
      LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
      GO       - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL INTERPOLATION
                 2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
                 3    UNRECOGNIZED OUTPUT GRID
                 41   INVALID NONGLOBAL INPUT GRID
                 42   INVALID SPECTRAL METHOD PARAMETERS
 

POLATES6

The POLATES6 routine performs budget interpolation
             from any grid to any grid for scalar fields.
             It requires a grid for the output fields (kgdso(1)>=0).
             The algorithm simply computes (weighted) averages
             of neighbor points arranged in a square box
             centered around each output grid point and stretching
             nearly halfway to each of the neighboring grid points.
             Options allow choices of number of points in each radius
             from the center point (ipopt(1)) which defaults to 2
             (if ipopt(1)=-1) meaning that 25 points will be averaged;
             Further options are the respective weights for the radius
             points starting at the center point (ipopt(2:2+ipopt(1))
             which defaults to all 1 (if ipopt(1)=-1 or ipopt(2)=-1).
             another option is the minimum percentage for mask,
             i.e. percent valid input data required to make output data,
             (ipopt(3+ipopt(1)) which defaults to 50 (if -1).
             Only horizontal interpolation is performed.
             The grids are defined by their grid description sections
             (passed in integer form as decoded by subprogram w3fi63).

             The current code recognizes the following projections:
                (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
                (KGDS(1)=001) MERCATOR CYLINDRICAL
                (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
                (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
                (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
                (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)

             Where KGDS could be either input KGDSI or output KGDSO.
             As an added bonus the number of output grid points
             and their latitudes and longitudes are also returned.
             Input bitmaps will be interpolated to output bitmaps.
             output bitmaps will also be created when the output grid
             extends outside of the domain of the input grid.
             the output field is set to 0 where the output bitmap is off.

USAGE:    CALL POLATES6(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
                        NO,RLAT,RLON,IBO,LO,GO,IRET)

   Input argument list:
     IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
                IPOPT(1) IS NUMBER OF RADIUS POINTS
                (DEFAULTS TO 2 IF IPOPT(1)=-1);
                IPOPT(2:2+IPOPT(1)) ARE RESPECTIVE WEIGHTS
                (DEFAULTS TO ALL 1 IF IPOPT(1)=-1 OR IPOPT(2)=-1).
                IPOPT(3+IPOPT(1)) IS MINIMUM PERCENTAGE FOR MASK
                (DEFAULTS TO 50 IF IPOPT(3+IPOPT(1)=-1)
     KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
     KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
     MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
                OR DIMENSION OF INPUT GRID FIELDS IF KM=1
     MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
                OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
     KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
     IBI      - INTEGER (KM) INPUT BITMAP FLAGS
     LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
     GI       - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE

  Output argument list:
     NO       - INTEGER NUMBER OF OUTPUT POINTS
     RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES
     RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES
     IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
     LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
     GO       - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
     IRET     - INTEGER RETURN CODE
                0    SUCCESSFUL INTERPOLATION
                2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
                3    UNRECOGNIZED OUTPUT GRID
                31   INVALID UNDEFINED OUTPUT GRID
                32   INVALID BUDGET METHOD PARAMETERS


POLATEV0

The POLATEV0 routine performs bilinear interpolation
             from any grid to any grid for vector fields.
             options allow varying the minimum percentage for mask,
             i.e. percent valid input data required to make output data,
             (ipopt(1)) which defaults to 50 (if ipopt(1)=-1).
             only horizontal interpolation is performed.
             The grids are defined by their grid description sections
             (passed in integer form as decoded by subprogram w3fi63).

             The current code recognizes the following projections:
                (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
                (KGDS(1)=001) MERCATOR CYLINDRICAL
                (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
                (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
                (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
                (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)

             Where KGDS could be either input KGDSI or output KGDSO.
             As an added bonus the number of output grid points
             and their latitudes and longitudes are also returned.
             Along with their vector rotation parameters.
             On the other hand, the output can be a set of station points
             If kgdso(1)<0, in which case the number of points
             and their latitudes and longitudes must be input
             Along with their vector rotation parameters.
             Input bitmaps will be interpolated to output bitmaps.
             output bitmaps will also be created when the output grid
             extends outside of the domain of the input grid.
             The output field is set to 0 where the output bitmap is off.

USAGE:    CALL POLATEV0(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,UI,VI,
                        NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)

   Input argument list:
      IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
                 IPOPT(1) IS MINIMUM PERCENTAGE FOR MASK
                 (DEFAULTS TO 50 IF IPOPT(1)=-1)
      KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
      KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
                 (KGDSO(1)<0 IMPLIES RANDOM STATION POINTS)
      MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF INPUT GRID FIELDS IF KM=1
      MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
      KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
      IBI      - INTEGER (KM) INPUT BITMAP FLAGS
      LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
      UI       - REAL (MI,KM) INPUT U-COMPONENT FIELDS TO INTERPOLATE
      VI       - REAL (MI,KM) INPUT V-COMPONENT FIELDS TO INTERPOLATE
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)
      RLAT     - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
      RLON     - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
      CROT     - REAL (NO) VECTOR ROTATION COSINES (IF KGDSO(1)<0)
      SROT     - REAL (NO) VECTOR ROTATION SINES (IF KGDSO(1)<0)
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
 
   Output argument list:
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
      RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
      RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
      CROT     - REAL (NO) VECTOR ROTATION COSINES (IF KGDSO(1)>=0)
      SROT     - REAL (NO) VECTOR ROTATION SINES (IF KGDSO(1)>=0)
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
      IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
      LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
      UO       - REAL (MO,KM) OUTPUT U-COMPONENT FIELDS INTERPOLATED
      VO       - REAL (MO,KM) OUTPUT V-COMPONENT FIELDS INTERPOLATED
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL INTERPOLATION
                 2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
                 3    UNRECOGNIZED OUTPUT GRID
 

POLATEV1

The POLATEV1 routine performs bicubic interpolation
             from any grid to any grid for scalar fields.
             Bitmaps are now allowed, but bilinear interpolation is done
             when any invalid data is within the bicubic template.
             Options allow choices between straight bicubic (ipopt(1)=0)
             and constrained bicubic (ipopt(1)=1) where the value is
             confined within the range of the surrounding 4 points.
             Another option is the minimum percentage for mask,
             i.e. percent valid input data required to make output data,
             (ipopt(2)) which defaults to 50 (if ipopt(2)=-1).
             Bilinear used within one grid length of boundaries.
             Only horizontal interpolation is performed.
             The grids are defined by their grid description sections
             (passed in integer form as decoded by subprogram w3fi63).

             The current code recognizes the following projections:
                (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
                (KGDS(1)=001) MERCATOR CYLINDRICAL
                (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
                (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
                (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
                (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)

             Where KGDS could be either input KGDSI or output KGDSO.
             As an added bonus the number of output grid points
             and their latitudes and longitudes are also returned.
             Along with their vector rotation parameters.
             On the other hand, the output can be a set of station points
             If kgdso(1)<0, in which case the number of points
             and their latitudes and longitudes must be input
             along with their vector rotation parameters.
             Output bitmaps will only be created when the output grid
             extends outside of the domain of the input grid.
             The output field is set to 0 where the output bitmap is off.

USAGE:    CALL POLATEV1(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,UI,VI,
                        NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)

   Input argument list:
      IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
                 IPOPT(1)=0 FOR STRAIGHT BICUBIC;
                 IPOPT(1)=1 FOR CONSTRAINED BICUBIC WHERE VALUE IS
                 CONFINED WITHIN THE RANGE OF THE SURROUNDING 4 POINTS.
                 IPOPT(2) IS MINIMUM PERCENTAGE FOR MASK
                 (DEFAULTS TO 50 IF IPOPT(2)=-1)
      KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
      KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
                 (KGDSO(1)<0 IMPLIES RANDOM STATION POINTS)
      MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF INPUT GRID FIELDS IF KM=1
      MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
      KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
      IBI      - INTEGER (KM) INPUT BITMAP FLAGS
      LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
      UI       - REAL (MI,KM) INPUT U-COMPONENT FIELDS TO INTERPOLATE
      VI       - REAL (MI,KM) INPUT V-COMPONENT FIELDS TO INTERPOLATE
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)
      RLAT     - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
      RLON     - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
      CROT     - REAL (NO) VECTOR ROTATION COSINES (IF KGDSO(1)<0)
      SROT     - REAL (NO) VECTOR ROTATION SINES (IF KGDSO(1)<0)
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
 
   Output argument list:
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
      RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
      RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
      CROT     - REAL (NO) VECTOR ROTATION COSINES (IF KGDSO(1)>=0)
      SROT     - REAL (NO) VECTOR ROTATION SINES (IF KGDSO(1)>=0)
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
      IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
      LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
      UO       - REAL (MO,KM) OUTPUT U-COMPONENT FIELDS INTERPOLATED
      VO       - REAL (MO,KM) OUTPUT V-COMPONENT FIELDS INTERPOLATED
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL INTERPOLATION
                 2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
                 3    UNRECOGNIZED OUTPUT GRID
 

POLATEV2

The POLATEV2 routine performs neighbor interpolation
             from any grid to any grid for vector fields.
             Options allow choosing the width of the grid square
             (ipopt(1)) to search for valid data, which defaults to 1
             (if ipopt(1)=-1).  Odd width squares are centered on
             the nearest input grid point; Even width squares are
             centered on the nearest four input grid points.
             Squares are searched for valid data in a spiral pattern
             starting from the center.  No searching is done where
             the output grid is outside the input grid.  Only horizontal
             interpolation is performed.
             The grids are defined by their grid description sections
             (passed in integer form as decoded by subprogram w3fi63).

             The current code recognizes the following projections:
                (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
                (KGDS(1)=001) MERCATOR CYLINDRICAL
                (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
                (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
                (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
                (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)

             Where KGDS could be either input KGDSI or output KGDSO.
             As an added bonus the number of output grid points
             and their latitudes and longitudes are also returned.
             Along with their vector rotation parameters.
             On the other hand, the output can be a set of station points
             If KGDSO(1)<0, in which case the number of points
             and their latitudes and longitudes must be input.
             Input bitmaps will be interpolated to output bitmaps.
             Output bitmaps will also be created when the output grid
             extends outside of the domain of the input grid.
             The output field is set to 0 where the output bitmap is off.

USAGE:    CALL POLATEV2(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,UI,VI,
                        NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)

   Input argument list:
      IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
                 IPOPT(1) IS WIDTH OF SQUARE TO EXAMINE IN SPIRAL SEARCH
                 (DEFAULTS TO 1 IF IPOPT(1)=-1)
      KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
      KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
                 (KGDSO(1)<0 IMPLIES RANDOM STATION POINTS)
      MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF INPUT GRID FIELDS IF KM=1
      MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
      KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
      IBI      - INTEGER (KM) INPUT BITMAP FLAGS
      LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
      UI       - REAL (MI,KM) INPUT U-COMPONENT FIELDS TO INTERPOLATE
      VI       - REAL (MI,KM) INPUT V-COMPONENT FIELDS TO INTERPOLATE
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)
      RLAT     - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
      RLON     - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
      CROT     - REAL (NO) VECTOR ROTATION COSINES (IF KGDSO(1)<0)
      SROT     - REAL (NO) VECTOR ROTATION SINES (IF KGDSO(1)<0)
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
 
   Output argument list:
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
      RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
      RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
      CROT     - REAL (NO) VECTOR ROTATION COSINES (IF KGDSO(1)>=0)
      SROT     - REAL (NO) VECTOR ROTATION SINES (IF KGDSO(1)>=0)
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
      IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
      LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
      UO       - REAL (MO,KM) OUTPUT U-COMPONENT FIELDS INTERPOLATED
      VO       - REAL (MO,KM) OUTPUT V-COMPONENT FIELDS INTERPOLATED
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL INTERPOLATION
                 2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
                 3    UNRECOGNIZED OUTPUT GRID
 

POLATEV3

The POLATEV3 routine performs budget interpolation
             from any grid to any grid for vector fields.
             It requires a grid for the output fields (kgdso(1)>=0).
             The algorithm simply computes (weighted) averages
             of Bilinearly interpolated points arranged in a square box
             centered around each output grid point and stretching
             nearly halfway to each of the neighboring grid points.
             Options allow choices of number of points in each radius
             from the center point (ipopt(1)) which defaults to 2
             (if ipopt(1)=-1) meaning that 25 points will be averaged;
             Further options are the respective weights for the radius
             points starting at the center point (ipopt(2:2+ipopt(1))
             which defaults to all 1 (if ipopt(1)=-1 or ipopt(2)=-1).
             A special interpolation is done if ipopt(2)=-2.
             In this case, the boxes stretch nearly all the way to
             each of the neighboring grid points and the weights
             are the adjoint of the Bilinear interpolation weights.
             This case gives Quasi-Second-Order Budget interpolation.
             Another option is the minimum percentage for mask,
             i.e. percent valid input data required to make output data,
             (ipopt(3+ipopt(1)) which defaults to 50 (if -1).
             Only horizontal interpolation is performed.
             The grids are defined by their grid description sections
             (passed in integer form as decoded by subprogram w3fi63).

             The current code recognizes the following projections:
                (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
                (KGDS(1)=001) MERCATOR CYLINDRICAL
                (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
                (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
                (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
                (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)

             Where KGDS could be either input KGDSI or output KGDSO.
             The input and output vectors are rotated so that they are
             either resolved relative to the defined grid
             in the direction of increasing x and y coordinates
             or resolved relative to easterly and northerly directions,
             as designated by their respective grid description sections.
             As an added bonus the number of output grid points
             and their latitudes and longitudes are also returned
             Along with their vector rotation parameters.
             Input bitmaps will be interpolated to output bitmaps.
             Output bitmaps will also be created when the output grid
             Extends outside of the domain of the input grid.
             The output field is set to 0 where the output bitmap is off.

USAGE:    CALL POLATEV3(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,UI,VI,
                        NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)

   Input argument list:
      IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
                 IPOPT(1) IS NUMBER OF RADIUS POINTS
                 (DEFAULTS TO 2 IF IPOPT(1)=-1);
                 IPOPT(2:2+IPOPT(1)) ARE RESPECTIVE WEIGHTS
                 (DEFAULTS TO ALL 1 IF IPOPT(1)=-1 OR IPOPT(2)=-1).
                 IPOPT(3+IPOPT(1)) IS MINIMUM PERCENTAGE FOR MASK
                 (DEFAULTS TO 50 IF IPOPT(3+IPOPT(1)=-1)
      KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
      KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
      MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF INPUT GRID FIELDS IF KM=1
      MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
      KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
      IBI      - INTEGER (KM) INPUT BITMAP FLAGS
      LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
      UI       - REAL (MI,KM) INPUT U-COMPONENT FIELDS TO INTERPOLATE
      VI       - REAL (MI,KM) INPUT V-COMPONENT FIELDS TO INTERPOLATE
 
  Output argument list:
      NO       - INTEGER NUMBER OF OUTPUT POINTS
      RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES
      RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES
      CROT     - REAL (NO) VECTOR ROTATION COSINES
      SROT     - REAL (NO) VECTOR ROTATION SINES
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
      IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
      LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
      UO       - REAL (MO,KM) OUTPUT U-COMPONENT FIELDS INTERPOLATED
      VO       - REAL (MO,KM) OUTPUT V-COMPONENT FIELDS INTERPOLATED
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL INTERPOLATION
                 2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
                 3    UNRECOGNIZED OUTPUT GRID
                 31   INVALID UNDEFINED OUTPUT GRID
                 32   INVALID BUDGET METHOD PARAMETERS
 

POLATEV4

The POLATEV4 routine performs spectral interpolation
             from any grid to any grid for vector fields.
             It requires that the input fields be uniformly global.
             Options allow choices between triangular shape (ipopt(1)=0)
             and Rhomboidal shape (ipopt(1)=1) which has no default;
             A second option is the truncation (ipopt(2)) which defaults
             to a sensible truncation for the input grid (if opt(2)=-1).
             Note that if the output grid is not found in a special list,
             then the transform back to grid is not very fast.
             This special list contains global cylindrical grids,
             Polar Stereographic grids centered at the pole and
             Mercator grids.  Only horizontal interpolation is performed.
             The grids are defined by their grid description sections
             (passed in integer form as decoded by subprogram w3fi63).

             The current code recognizes the following projections:
                (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
                (KGDS(1)=001) MERCATOR CYLINDRICAL
                (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
                (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
                (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
                (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)

             Where KGDS could be either input KGDSI or output KGDSO.
             The input and output vectors are rotated so that they are
             either resolved relative to the defined grid.
             in the direction of increasing x and y coordinates
             or resolved relative to easterly and northerly directions,
             As designated by their respective grid description sections.
             As an added bonus the number of output grid points and
             their latitudes and longitudes are also returned along with
             their vector rotation parameters.
             On the other hand, the output can be a set of station points
             if kgdso(1)<0, in which case the number of points
             and their latitudes and longitudes must be input
             along with their vector rotation parameters.
             output bitmaps will only be created when the output grid
             extends outside of the domain of the input grid.
             The output field is set to 0 where the output bitmap is off.

USAGE:    CALL POLATEV4(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,UI,VI,
    &                    NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)

   Input argument list:
      IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
                 IPOPT(1)=0 FOR TRIANGULAR, IPOPT(1)=1 FOR RHOMBOIDAL;
                 IPOPT(2) IS TRUNCATION NUMBER
                 (DEFAULTS TO SENSIBLE IF IPOPT(2)=-1).
      KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
      KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
                 (KGDSO(1)<0 IMPLIES RANDOM STATION POINTS)
      MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF INPUT GRID FIELDS IF KM=1
      MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
      KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
      IBI      - INTEGER (KM) INPUT BITMAP FLAGS (MUST BE ALL 0)
      LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
      UI       - REAL (MI,KM) INPUT U-COMPONENT FIELDS TO INTERPOLATE
      VI       - REAL (MI,KM) INPUT V-COMPONENT FIELDS TO INTERPOLATE
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)
      RLAT     - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
      RLON     - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
      CROT     - REAL (NO) VECTOR ROTATION COSINES (IF KGDSO(1)<0)
      SROT     - REAL (NO) VECTOR ROTATION SINES (IF KGDSO(1)<0)
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
 
   Output argument list:
      NO       - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
      RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
      RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
      CROT     - REAL (NO) VECTOR ROTATION COSINES (IF KGDSO(1)>=0)
      SROT     - REAL (NO) VECTOR ROTATION SINES (IF KGDSO(1)>=0)
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
      IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
      LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
      UO       - REAL (MO,KM) OUTPUT U-COMPONENT FIELDS INTERPOLATED
      VO       - REAL (MO,KM) OUTPUT V-COMPONENT FIELDS INTERPOLATED
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL INTERPOLATION
                 2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
                 3    UNRECOGNIZED OUTPUT GRID
                 41   INVALID NONGLOBAL INPUT GRID
                 42   INVALID SPECTRAL METHOD PARAMETERS
 

POLATEV6

The POLATEV6 routine performs budget interpolation
             from any grid to any grid for vector fields.
             It requires a grid for the output fields (kgdso(1)>=0).
             The algorithm simply computes (weighted) averages
             of neighbor points arranged in a square box
             centered around each output grid point and stretching
             nearly halfway to each of the neighboring grid points.
             Options allow choices of number of points in each radius
             from the center point (IPOPT(1)) which defaults to 2
             (if ipopt(1)=-1) meaning that 25 points will be averaged;
             Further options are the respective weights for the radius
             points starting at the center point (ipopt(2:2+ipopt(1))
             which defaults to all 1 (if ipopt(1)=-1 or ipopt(2)=-1).
             Another option is the minimum percentage for mask,
             i.e. percent valid input data required to make output data,
             (ipopt(3+ipopt(1)) which defaults to 50 (if -1).
             only horizontal interpolation is performed.
             The grids are defined by their grid description sections
             (passed in integer form as decoded by subprogram w3fi63).

             The current code recognizes the following projections:
                (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
                (KGDS(1)=001) MERCATOR CYLINDRICAL
                (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
                (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
                (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
                (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)

             Where KGDS could be either input KGDSI or output KGDSO.
             The input and output vectors are rotated so that they are
             either resolved relative to the defined grid 
             in the direction of increasing x and y coordinates
             or resolved relative to easterly and northerly directions,
             As designated by their respective grid description sections.
             As an added bonus the number of output grid points
             and their latitudes and longitudes are also returned
             along with their vector rotation parameters.
             input bitmaps will be interpolated to output bitmaps.
             output bitmaps will also be created when the output grid
             extends outside of the domain of the input grid.
             The output field is set to 0 where the output bitmap is off.

USAGE:    CALL POLATEV6(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,UI,VI,
                        NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)

   Input argument list:
      IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
                 IPOPT(1) IS NUMBER OF RADIUS POINTS
                 (DEFAULTS TO 2 IF IPOPT(1)=-1);
                 IPOPT(2:2+IPOPT(1)) ARE RESPECTIVE WEIGHTS
                 (DEFAULTS TO ALL 1 IF IPOPT(1)=-1 OR IPOPT(2)=-1).
                 IPOPT(3+IPOPT(1)) IS MINIMUM PERCENTAGE FOR MASK
                 (DEFAULTS TO 50 IF IPOPT(3+IPOPT(1)=-1)
      KGDSI    - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
      KGDSO    - INTEGER (200) OUTPUT GDS PARAMETERS
      MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF INPUT GRID FIELDS IF KM=1
      MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
                 OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
      KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
      IBI      - INTEGER (KM) INPUT BITMAP FLAGS
      LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
      UI       - REAL (MI,KM) INPUT U-COMPONENT FIELDS TO INTERPOLATE
      VI       - REAL (MI,KM) INPUT V-COMPONENT FIELDS TO INTERPOLATE
 
   Output argument list:
      NO       - INTEGER NUMBER OF OUTPUT POINTS
      RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES
      RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES
      CROT     - REAL (NO) VECTOR ROTATION COSINES
      SROT     - REAL (NO) VECTOR ROTATION SINES
                 (UGRID=CROT*UEARTH-SROT*VEARTH;
                  VGRID=SROT*UEARTH+CROT*VEARTH)
      IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
      LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
      UO       - REAL (MO,KM) OUTPUT U-COMPONENT FIELDS INTERPOLATED
      VO       - REAL (MO,KM) OUTPUT V-COMPONENT FIELDS INTERPOLATED
      IRET     - INTEGER RETURN CODE
                 0    SUCCESSFUL INTERPOLATION
                 2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
                 3    UNRECOGNIZED OUTPUT GRID
                 31   INVALID UNDEFINED OUTPUT GRID
                 32   INVALID BUDGET METHOD PARAMETERS
 

POLFIXS

The POLFIXS routine averages multiple pole scalar values
            on a latitude/longitude grid.  Bitmaps may be averaged too.

USAGE:    CALL POLFIXS(NM,NX,KM,RLAT,RLON,IB,LO,GO)

   Input argument list:
      NO       - INTEGER NUMBER OF GRID POINTS
      NX       - INTEGER LEADING DIMENSION OF FIELDS
      KM       - INTEGER NUMBER OF FIELDS
      RLAT     - REAL (NO) LATITUDES IN DEGREES
      RLON     - REAL (NO) LONGITUDES IN DEGREES
      IB       - INTEGER (KM) BITMAP FLAGS
      LO       - LOGICAL*1 (NX,KM) BITMAPS (IF SOME IB(K)=1)
      GO       - REAL (NX,KM) FIELDS
 
   Output argument list:
      LO       - LOGICAL*1 (NX,KM) BITMAPS (IF SOME IB(K)=1)
      GO       - REAL (NX,KM) FIELDS
 

POLFIXV

The POLFIXV routine averages multiple pole vector values
            on a latitude/longitude grid.  Bitmaps may be averaged too.
            Vectors are rotated with respect to their longitude.

USAGE:    CALL POLFIXV(NM,NX,KM,RLAT,RLON,IB,LO,UO,VO)

   Input argument list:
      NO       - INTEGER NUMBER OF GRID POINTS
      NX       - INTEGER LEADING DIMENSION OF FIELDS
      KM       - INTEGER NUMBER OF FIELDS
      RLAT     - REAL (NO) LATITUDES IN DEGREES
      RLON     - REAL (NO) LONGITUDES IN DEGREES
      IB       - INTEGER (KM) BITMAP FLAGS
      LO       - LOGICAL*1 (NX,KM) BITMAPS (IF SOME IB(K)=1)
      UO       - REAL (NX,KM) U-WINDS
      VO       - REAL (NX,KM) V-WINDS
 
   Output argument list:
      LO       - LOGICAL*1 (NX,KM) BITMAPS (IF SOME IB(K)=1)
      UO       - REAL (NX,KM) U-WINDS
      VO       - REAL (NX,KM) V-WINDS
 
 
iplib.tar Library contains subroutines to be used for interpolating almost any grids used at NCEP. (Fortran90)
Date posted: 9/05/2012