gdswzd_mod.F90

Go to the documentation of this file.

!> @file
!! @brief Driver module for gdswzd routines.
!!
!! @date Jan 2015
!! @author George Gayno, Kyle Gerheiser
 
!> @brief Driver module for gdswzd routines.
!!
!! These routines do the following for several map projections:
!! - Convert from earth to grid coordinates or vice versa.
!! - Compute vector rotation sines and cosines.
!! - Compute map jacobians.
!! - Compute grid box area.
!!
!! Map projections include:
!! - Equidistant Cyclindrical
!! - Mercator Cylindrical
!! - Gaussian Cylindrical
!! - Polar stereographic
!! - Lambert Conformal Conic
!! - Rotated Equidistant Cyclindrical ("E" and non-"E" staggers)
!!
!! @author Mark Iredell, George Gayno, Kyle Gerheiser
!! @date Jan 2015
MODULE gdswzd_mod
  use ip_grid_descriptor_mod
  use ip_grids_mod
  use ip_grid_mod
  use ip_grid_factory_mod
  
  IMPLICIT NONE
 
  PRIVATE
 
  public :: gdswzd_2d_array_grib1, gdswzd_grib1, gdswzd
 
  INTERFACE gdswzd
     MODULE PROCEDURE gdswzd_1d_array
     MODULE PROCEDURE gdswzd_2d_array
     MODULE PROCEDURE gdswzd_scalar
     module procedure gdswzd_grib1
     module procedure gdswzd_2d_array_grib1
     module procedure gdswzd_grid
  INTERFACE gdswzd …
  END INTERFACE gdswzd
 
CONTAINS
 
  !> Returns one of the following for a grid object:
  !! - 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.
  !!
  !! 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, map jacobians and
  !! grid box areas may be returned.
  !!
  !! To compute the vector rotations, the optional arguments 'srot' and 'crot'
  !! must be present. To compute the map jacobians, the
  !! optional arguments 'xlon', 'xlat', 'ylon', 'ylat' must be present.
  !!
  !! To compute the grid box areas, the optional argument
  !! 'area' must be present.
  !!
  !! @param[in] grid Grid to call gdswzd on.
  !!
  !! @param[in] iopt Option flag.
  !! - 0 Earth coords of all the grid points.
  !! - 1 Earth coords of selected grid coords.
  !! - -1 Grid coords of selected earth coords
  !!
  !! @param[in] npts Maximum number of coordinates.
  !! @param[in] fill Fill value to set invalid output data.
  !!      Must be impossible value; suggested value: -9999.
  !! @param[inout] xpts Grid x point coordinates.
  !! @param[inout] ypts Grid y point coordinates.
  !! @param[inout] rlon Earth longitudes in degrees E.
  !! (Acceptable range: -360. to 360.)
  !!
  !! @param[inout] rlat Earth latitudes in degrees N.
  !! (Acceptable range: -90. to 90.)
  !!
  !! @param[out] nret Number of valid points computed.
  !! @param[out] crot Clockwise vector rotation cosines.
  !! @param[out] srot Clockwise vector rotation sines.
  !! ugrid=crot*uearth-srot*vearth;
  !! vgrid=srot*uearth+crot*vearth)
  !!
  !! @param[out] xlon dx/dlat in 1/degrees
  !! @param[out] xlat dy/dlon in 1/degrees
  !! @param[out] ylon dy/dlon in 1/degrees
  !! @param[out] ylat dy/dlat in 1/degrees
  !! @param[out] area Area weights in m^2.
  !! Proportional to the square of the map factor in the case of
  !! conformal projections
  !!
  !! @author Kyle Gerheiser
  !! @date July 2021
  subroutine gdswzd_grid(grid,IOPT,NPTS,FILL, &
       XPTS,YPTS,RLON,RLAT,NRET, &
       CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
 
    class(ip_grid), intent(in) :: grid
    INTEGER,        INTENT(IN   ) :: IOPT, NPTS
    INTEGER,        INTENT(  OUT) :: NRET
    !
    REAL,           INTENT(IN   ) :: FILL
    REAL,           INTENT(INOUT) :: RLON(NPTS),RLAT(NPTS)
    REAL,           INTENT(INOUT) :: XPTS(NPTS),YPTS(NPTS)
    REAL, OPTIONAL, INTENT(  OUT) :: CROT(NPTS),SROT(NPTS)
    REAL, OPTIONAL, INTENT(  OUT) :: XLON(NPTS),XLAT(NPTS)
    REAL, OPTIONAL, INTENT(  OUT) :: YLON(NPTS),YLAT(NPTS),AREA(NPTS)
 
    INTEGER                       :: IS1, IM, JM, NM, KSCAN, NSCAN, N
    INTEGER                       :: IOPF, NN, I, J
 
    !  COMPUTE GRID COORDINATES FOR ALL GRID POINTS
    IF(iopt.EQ.0) THEN
       iopf=1
 
       if(grid%descriptor%grid_num.eq.-1)then
          nm = npts
       else
          im = grid%im
          jm = grid%jm
          nm = im * jm
       endif
       nscan = grid%nscan
       kscan = grid%kscan
 
       if (nm > npts) then
          rlat=fill
          rlon=fill
          xpts=fill
          ypts=fill
          return
       end if
 
       select type(grid)
       type is(ip_rot_equid_cylind_egrid)
          if(kscan == 0) then
             is1 = (jm + 1) / 2
          else
             is1 = jm / 2
          end if
 
          DO n=1,nm
             IF(nscan.EQ.0) THEN
                j=(n-1)/im+1
                i=(n-im*(j-1))*2-mod(j+kscan,2)
             ELSE
                nn=(n*2)-1+kscan
                i = (nn-1)/jm + 1
                j = mod(nn-1,jm) + 1
                IF (mod(jm,2)==0.AND.mod(i,2)==0.AND.kscan==0) j = j + 1
                IF (mod(jm,2)==0.AND.mod(i,2)==0.AND.kscan==1) j = j - 1
             ENDIF
             xpts(n)=is1+(i-(j-kscan))/2
             ypts(n)=(i+(j-kscan))/2
          ENDDO
       type is(ip_station_points_grid)
          DO n=1,nm
             xpts(n)=fill
             ypts(n)=fill
          ENDDO
       class default
          DO n=1,nm
             IF(nscan.EQ.0) THEN
                j=(n-1)/im+1
                i=n-im*(j-1)
             ELSE
                i=(n-1)/jm+1
                j=n-jm*(i-1)
             ENDIF
             xpts(n)=i
             ypts(n)=j
          ENDDO
       end select
 
       DO n=nm+1,npts
          xpts(n)=fill
          ypts(n)=fill
       ENDDO
 
    ELSE  ! IOPT /= 0
       iopf=iopt
    ENDIF ! IOPT CHECK
 
    call grid%gdswzd(iopf,npts,fill,  &
            xpts,ypts,rlon,rlat,nret, &
            crot,srot,xlon,xlat,ylon,ylat,area)
 
  subroutine gdswzd_grid(grid,IOPT,NPTS,FILL, & …
  end subroutine gdswzd_grid
 
 
  !> Decodes the grib 2 grid definition template and returns
  !! one of the following (for scalars):
  !! - 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,
  !! where "igdtnum" is the grid definition template number:
  !! - igdtnum=00 Equidistant Cylindrical
  !! - igdtnum=01 Rotated Equidistant Cylindrical. "E" and non-"E" staggered
  !! - igdtnum=10 Mercator Cyclindrical
  !! - igdtnum=20 Polar Stereographic Azimuthal
  !! - igdtnum=30 Lambert Conformal Conical
  !! - igdtnum=40 Gaussian Equidistant Cyclindrical
  !!
  !! 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, map jacobians and
  !! grid box areas may be returned.
  !!
  !! To compute the vector rotations, the optional arguments 'srot' and 'crot'
  !! must be present. To compute the map jacobians, the
  !! optional arguments 'xlon', 'xlat', 'ylon', 'ylat' must be present.
  !!
  !! To compute the grid box areas, the optional argument
  !! 'area' must be present.
  !!
  !! @param[in] igdtnum Grid definition template number.
  !!
  !! @param[in] igdtmpl Grid definition template array.
  !! Corresponds to the gfld%igdtmpl component of the
  !! NCEPLIBS-g2 gridmod data structure
  !! See igdtmpl definition in gdswzd_1d_array() for full details.
  !!
  !! @param[in] igdtlen Number of elements of the grid definition
  !! template array. Corresponds to the gfld%igdtlen
  !! component of the ncep g2 library gridmod data structure.
  !!
  !! @param[in] iopt Option flag.
  !! - 0 Earth coords of all the grid points.
  !! - 1 Earth coords of selected grid coords.
  !! - -1 Grid coords of selected earth coords
  !!
  !! @param[in] npts Maximum number of coordinates.
  !! @param[in] fill Fill value to set invalid output data.
  !!      Must be impossible value; suggested value: -9999.
  !! @param[inout] xpts Grid x point coordinates.
  !! @param[inout] ypts Grid y point coordinates.
  !! @param[inout] rlon Earth longitudes in degrees E.
  !! (Acceptable range: -360. to 360.)
  !!
  !! @param[inout] rlat Earth latitudes in degrees N.
  !! (Acceptable range: -90. to 90.)
  !!
  !! @param[out] nret Number of valid points computed.
  !! @param[out] crot Clockwise vector rotation cosines.
  !! @param[out] srot Clockwise vector rotation sines.
  !! ugrid=crot*uearth-srot*vearth;
  !! vgrid=srot*uearth+crot*vearth)
  !!
  !! @param[out] xlon dx/dlat in 1/degrees
  !! @param[out] xlat dy/dlon in 1/degrees
  !! @param[out] ylon dy/dlon in 1/degrees
  !! @param[out] ylat dy/dlat in 1/degrees
  !! @param[out] area Area weights in m^2.
  !! Proportional to the square of the map factor in the case of
  !! conformal projections
  !!
  !! @author George Gayno, Mark Iredell
  !! @date Jan 2015
  SUBROUTINE gdswzd_scalar(IGDTNUM,IGDTMPL,IGDTLEN,IOPT,NPTS,FILL, &
       XPTS,YPTS,RLON,RLAT,NRET, &
       CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
 
    IMPLICIT NONE
    !
    INTEGER,        INTENT(IN   ) :: IGDTNUM, IGDTLEN
    INTEGER,        INTENT(IN   ) :: IGDTMPL(IGDTLEN)
    INTEGER,        INTENT(IN   ) :: IOPT, NPTS
    INTEGER,        INTENT(  OUT) :: NRET
    !
    REAL,           INTENT(IN   ) :: FILL
    REAL,           INTENT(INOUT) :: RLON, RLAT
    REAL,           INTENT(INOUT) :: XPTS, YPTS
    REAL, OPTIONAL, INTENT(  OUT) :: CROT, SROT
    REAL, OPTIONAL, INTENT(  OUT) :: XLON, XLAT
    REAL, OPTIONAL, INTENT(  OUT) :: YLON, YLAT, AREA
 
    REAL                          :: RLONA(1),RLATA(1)
    REAL                          :: XPTSA(1),YPTSA(1)
    REAL                          :: CROTA(1),SROTA(1)
    REAL                          :: XLONA(1),XLATA(1)
    REAL                          :: YLONA(1),YLATA(1),AREAA(1)
 
    rlona(1) = rlon
    rlata(1) = rlat
    xptsa(1) = xpts
    yptsa(1) = ypts
 
    nret = 0
 
    ! CALL WITHOUT EXTRA FIELDS.
 
    IF (.NOT. PRESENT(crot) .AND. &
         .NOT. PRESENT(srot) .AND. &
         .NOT. PRESENT(xlon) .AND. &
         .NOT. PRESENT(xlat) .AND. &
         .NOT. PRESENT(ylon) .AND. &
         .NOT. PRESENT(ylat) .AND. &
         .NOT. PRESENT(area) ) THEN
 
       CALL gdswzd_1d_array(igdtnum,igdtmpl,igdtlen,iopt,npts,fill, &
            xptsa,yptsa,rlona,rlata,nret)
 
       rlon = rlona(1)
       rlat = rlata(1)
       xpts = xptsa(1)
       ypts = yptsa(1)
 
    ENDIF
 
    ! MIMIC CALL TO OLD 'GDSWIZ' ROUTINES.
 
    IF (PRESENT(crot) .AND. &
         PRESENT(srot) .AND. &
         .NOT. PRESENT(xlon) .AND. &
         .NOT. PRESENT(xlat) .AND. &
         .NOT. PRESENT(ylon) .AND. &
         .NOT. PRESENT(ylat) .AND. &
         .NOT. PRESENT(area) ) THEN
 
       CALL gdswzd_1d_array(igdtnum,igdtmpl,igdtlen,iopt,npts,fill, &
            xptsa,yptsa,rlona,rlata,nret,crota,srota)
 
       rlon = rlona(1)
       rlat = rlata(1)
       xpts = xptsa(1)
       ypts = yptsa(1)
       crot = crota(1)
       srot = srota(1)
 
    ENDIF
 
    ! MIMIC CALL TO OLD 'GDSWZD' ROUTINES.
 
    IF (PRESENT(crot) .AND. &
         PRESENT(srot) .AND. &
         PRESENT(xlon) .AND. &
         PRESENT(xlat) .AND. &
         PRESENT(ylon) .AND. &
         PRESENT(ylat) .AND. &
         PRESENT(area) ) THEN
 
       CALL gdswzd_1d_array(igdtnum,igdtmpl,igdtlen,iopt,npts,fill, &
            xptsa,yptsa,rlona,rlata,nret, &
            crota,srota,xlona,xlata,ylona,ylata,areaa)
 
       rlon = rlona(1)
       rlat = rlata(1)
       xpts = xptsa(1)
       ypts = yptsa(1)
       crot = crota(1)
       srot = srota(1)
       xlon = xlona(1)
       xlat = xlata(1)
       ylon = ylona(1)
       ylat = ylata(1)
       area = areaa(1)
 
    ENDIF
 
    RETURN
 
  SUBROUTINE gdswzd_scalar(IGDTNUM,IGDTMPL,IGDTLEN,IOPT,NPTS,FILL, & …
  END SUBROUTINE gdswzd_scalar
 
  !> Decodes the grib 2 grid definition template and returns
  !! one of the following (for 2d-arrays):
  !! - 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,
  !! where "igdtnum" is the grid definition template number:
  !! - igdtnum=00 Equidistant Cylindrical
  !! - igdtnum=01 Rotated Equidistant Cylindrical. "E" and non-"E" staggered
  !! - igdtnum=10 Mercator Cyclindrical
  !! - igdtnum=20 Polar Stereographic Azimuthal
  !! - igdtnum=30 Lambert Conformal Conical
  !! - igdtnum=40 Gaussian Equidistant Cyclindrical
  !!
  !! 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, map jacobians and
  !! grid box areas may be returned.
  !!
  !! To compute the vector rotations, the optional arguments 'srot' and 'crot'
  !! must be present. To compute the map jacobians, the
  !! optional arguments 'xlon', 'xlat', 'ylon', 'ylat' must be present.
  !!
  !! To compute the grid box areas, the optional argument
  !! 'area' must be present.
  !!
  !! @param[in] igdtnum Grid definition template number.
  !!
  !! @param[in] igdtmpl Grid definition template array.
  !! Corresponds to the gfld%igdtmpl component of the
  !! NCEPLIBS-g2 gridmod data structure.
  !! See igdtmpl definition in gdswzd_1d_array() for full details.
  !!
  !! @param[in] igdtlen Number of elements of the grid definition
  !! template array. Corresponds to the gfld%igdtlen
  !! component of the ncep g2 library gridmod data structure.
  !!
  !! @param[in] iopt Option flag.
  !! - 0 Earth coords of all the grid points.
  !! - 1 Earth coords of selected grid coords.
  !! - -1 Grid coords of selected earth coords
  !!
  !! @param[in] npts Maximum number of coordinates.
  !! @param[in] fill Fill value to set invalid output data.
  !!      Must be impossible value; suggested value: -9999.
  !! @param[inout] xpts Grid x point coordinates.
  !! @param[inout] ypts Grid y point coordinates.
  !! @param[inout] rlon Earth longitudes in degrees E.
  !! (Acceptable range: -360. to 360.)
  !!
  !! @param[inout] rlat Earth latitudes in degrees N.
  !! (Acceptable range: -90. to 90.)
  !!
  !! @param[out] nret Number of valid points computed.
  !! @param[out] crot Clockwise vector rotation cosines.
  !! @param[out] srot Clockwise vector rotation sines.
  !! ugrid=crot*uearth-srot*vearth;
  !! vgrid=srot*uearth+crot*vearth)
  !!
  !! @param[out] xlon dx/dlat in 1/degrees
  !! @param[out] xlat dy/dlon in 1/degrees
  !! @param[out] ylon dy/dlon in 1/degrees
  !! @param[out] ylat dy/dlat in 1/degrees
  !! @param[out] area Area weights in m^2.
  !! Proportional to the square of the map factor in the case of
  !! conformal projections
  !!
  !! @author George Gayno, Mark Iredell
  !! @date Jan 2015
  SUBROUTINE gdswzd_2d_array(IGDTNUM,IGDTMPL,IGDTLEN,IOPT,NPTS,FILL, &
       XPTS,YPTS,RLON,RLAT,NRET, &
       CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
 
    IMPLICIT NONE
    !
    INTEGER,        INTENT(IN   ) :: IGDTNUM, IGDTLEN
    INTEGER,        INTENT(IN   ) :: IGDTMPL(IGDTLEN)
    INTEGER,        INTENT(IN   ) :: IOPT, NPTS
    INTEGER,        INTENT(  OUT) :: NRET
    !
    REAL,           INTENT(IN   ) :: FILL
    REAL,           INTENT(INOUT) :: RLON(:,:),RLAT(:,:)
    REAL,           INTENT(INOUT) :: XPTS(:,:),YPTS(:,:)
    REAL, OPTIONAL, INTENT(  OUT) :: CROT(:,:),SROT(:,:)
    REAL, OPTIONAL, INTENT(  OUT) :: XLON(:,:),XLAT(:,:)
    REAL, OPTIONAL, INTENT(  OUT) :: YLON(:,:),YLAT(:,:),AREA(:,:)
 
    CALL gdswzd_1d_array(igdtnum,igdtmpl,igdtlen,iopt,npts,fill, &
         xpts,ypts,rlon,rlat,nret, &
         crot,srot,xlon,xlat,ylon,ylat,area)
 
  SUBROUTINE gdswzd_2d_array(IGDTNUM,IGDTMPL,IGDTLEN,IOPT,NPTS,FILL, & …
  END SUBROUTINE gdswzd_2d_array
 
  !> Decodes the grib 2 grid definition template 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,
  !! where "igdtnum" is the grid definition template number:
  !! - igdtnum=00 Equidistant Cylindrical
  !! - igdtnum=01 Rotated Equidistant Cylindrical. "E" and non-"E" staggered
  !! - igdtnum=10 Mercator Cyclindrical
  !! - igdtnum=20 Polar Stereographic Azimuthal
  !! - igdtnum=30 Lambert Conformal Conical
  !! - igdtnum=40 Gaussian Equidistant Cyclindrical
  !!
  !! 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, map jacobians and
  !! grid box areas may be returned.
  !!
  !! To compute the vector rotations, the optional arguments 'srot' and 'crot'
  !! must be present. To compute the map jacobians, the
  !! optional arguments 'xlon', 'xlat', 'ylon', 'ylat' must be present.
  !!
  !! To compute the grid box areas, the optional argument
  !! 'area' must be present.
  !!
  !! @param[in] igdtnum Grid definition template number.
  !! Corresponds to the gfld%igdtnum component of the ncep g2 library
  !! gridmod data structure:
  !! - 00 - Equidistant Cylindrical
  !! - 01 - Rotated Equidistant Cylindrical. "E" and non-"E" staggered
  !! - 10 - Mercator Cyclindrical
  !! - 20 - Polar Stereographic Azimuthal
  !! - 30 - Lambert Conformal Conical
  !! - 40 - Gaussian Equidistant Cyclindrical
  !!
  !! @param[in] igdtmpl Grid definition template array.
  !! Corresponds to the gfld%igdtmpl component of the
  !! NCEPLIBS-g2 gridmod data structure
  !!
  !! Section 3 Info:
  !!
  !! All Map Projections:
  !! - 1: Shape of earth, octet 15.
  !! - 2: Scale factor of spherical earth radius, octet 16.
  !! - 3: Scaled value of radius of spherical earth, octets 17-20.
  !! - 4: Scale factor of major axis of elliptical earth, octet 21.
  !! - 5: Scaled value of major axis of elliptical earth, octets 22-25.
  !! - 6: Scale factor of minor axis of elliptical earth, octet 26.
  !! - 7: Scaled value of minor axis of elliptical earth, octets 27-30.
  !!
  !! Equidistant Cyclindrical:
  !! - 8:  Number of points along a parallel, octs 31-34.
  !! - 9:  Number of points along a meridian, octs 35-38.
  !! - 10: Basic angle of initial production domain, octets 39-42.
  !! - 11: Subdivisions of basic angle, octets 43-46.
  !! - 12: Latitude of first grid point, octets 47-50.
  !! - 13: Longitude of first grid point, octets 51-54.
  !! - 14: Resolution and component flags, octet 55.
  !! - 15: Latitude of last grid point, octets 56-59.
  !! - 16: Longitude of last grid point, octets 60-63.
  !! - 17: i-direction increment, octets 64-67.
  !! - 18: j-direction increment, octets 68-71.
  !! - 19: Scanning mode, octet 72.
  !!
  !! Mercator Cyclindrical:
  !! - 8:  Number of points along a parallel, octs 31-34.
  !! - 9:  Number of points along a meridian, octs 35-38.
  !! - 10: Latitude of first point, octets 39-42.
  !! - 11: Longitude of first point, octets 43-46.
  !! - 12: Resolution and component flags, octet 47.
  !! - 13: Tangent latitude, octets 48-51.
  !! - 14: Latitude of last point, octets 52-55.
  !! - 15: Longitude of last point, octets 56-59.
  !! - 16: Scanning mode flags, octet 60.
  !! - 17: Orientation of grid, octets 61-64.
  !! - 18: Longitudinal grid length, octets 65-68.
  !! - 19: Latitudinal grid length, octets 69-72.
  !!
  !! Lambert Conformal Conical:
  !! - 8:  Number of points along x-axis, octs 31-34.
  !! - 9:  Number of points along y-axis, octs 35-38.
  !! - 10: Latitude of first point, octets 39-42.
  !! - 11: Longitude of first point, octets 43-46.
  !! - 12: Resolution of component flag, octet 47.
  !! - 13: Latitude where grid lengths specified,octets 48-51.
  !! - 14: Longitude of meridian that is parallel to y-axis, octets 52-55.
  !! - 15: x-direction grid length, octets 56-59.
  !! - 16: y-direction grid length, octets 60-63.
  !! - 17: Projection center flag, octet 64.
  !! - 18: Scanning mode, octet 65.
  !! - 19: First tangent latitude from pole, octets 66-69.
  !! - 20: Second tangent latitude from pole, octets 70-73.
  !! - 21: Latitude of south pole of projection, octets 74-77.
  !! - 22: Longitude of south pole of projection, octets 78-81.
  !!
  !! Gaussian Cylindrical:
  !! - 8:  Number of points along a parallel, octs 31-34.
  !! - 9:  Number of points along a meridian, octs 35-38.
  !! - 10: Basic angle of initial production domain, octets 39-42.
  !! - 11: Subdivisions of basic angle, octets 43-46.
  !! - 12: Latitude of first grid point, octets 47-50.
  !! - 13: Longitude of first grid point, octets 51-54.
  !! - 14: Resolution and component flags, octet 55.
  !! - 15: Latitude of last grid point, octets 56-59.
  !! - 16: Longitude of last grid point, octets 60-63.
  !! - 17: i-direction increment, octets 64-67.
  !! - 18: Number of parallels between pole and equator, octets 68-71.
  !! - 19: Scanning mode, octet 72.
  !!
  !! Polar Stereographic Azimuthal:
  !! - 8:  Number of points along x-axis, octets 31-34.
  !! - 9:  Number of points along y-axis, octets 35-38.
  !! - 10: Latitude of first grid point, octets 39-42.
  !! - 11: Longitude of first grid point, octets 43-46.
  !! - 12: Resolution and component flags, octet 47.
  !! - 13: True latitude, octets 48-51.
  !! - 14: Orientation longitude, octets 52-55.
  !! - 15: x-direction grid length, octets 56-59.
  !! - 16: y-direction grid length, octets 60-63.
  !! - 17: Projection center flag, octet 64.
  !! - 18: Scanning mode flags, octet 65.
  !!
  !! Rotated Equidistant Cyclindrical:
  !! - 8:  Number of points along a parallel, octs 31-34.
  !! - 9:  Number of points along a meridian, octs 35-38.
  !! - 10: Basic angle of initial production domain, octets 39-42.
  !! - 11: Subdivisions of basic angle, octets 43-46.
  !! - 12: Latitude of first grid point, octets 47-50.
  !! - 13: Longitude of first grid point, octets 51-54.
  !! - 14: Resolution and component flags, octet 55.
  !! - 15: Latitude of last grid point, octets 56-59.
  !! - 16: Longitude of last grid point, octets 60-63.
  !! - 17: i-direction increment, octets 64-67.
  !! - 18: j-direction increment, octets 68-71.
  !! - 19: Scanning mode, octet 72.
  !! - 20: Latitude of southern pole of projection, octets 73-76.
  !! - 21: Longitude of southern pole of projection, octets 77-80.
  !! - 22: Angle of rotation of projection, octs 81-84.
  !!
  !! @param[in] igdtlen Number of elements of the grid definition
  !! template array. Corresponds to the gfld%igdtlen
  !! component of the ncep g2 library gridmod data structure.
  !!
  !! @param[in] iopt Option flag.
  !! - 0 Earth coords of all the grid points.
  !! - 1 Earth coords of selected grid coords.
  !! - -1 Grid coords of selected earth coords
  !!
  !! @param[in] npts Maximum number of coordinates.
  !! @param[in] fill Fill value to set invalid output data.
  !!      Must be impossible value; suggested value: -9999.
  !! @param[inout] xpts Grid x point coordinates.
  !! @param[inout] ypts Grid y point coordinates.
  !! @param[inout] rlon Earth longitudes in degrees E.
  !! (Acceptable range: -360. to 360.)
  !!
  !! @param[inout] rlat Earth latitudes in degrees N.
  !! (Acceptable range: -90. to 90.)
  !!
  !! @param[out] nret Number of valid points computed.
  !! @param[out] crot Clockwise vector rotation cosines.
  !! @param[out] srot Clockwise vector rotation sines.
  !! ugrid=crot*uearth-srot*vearth;
  !! vgrid=srot*uearth+crot*vearth)
  !!
  !! @param[out] xlon dx/dlat in 1/degrees
  !! @param[out] xlat dy/dlon in 1/degrees
  !! @param[out] ylon dy/dlon in 1/degrees
  !! @param[out] ylat dy/dlat in 1/degrees
  !! @param[out] area Area weights in m^2.
  !! Proportional to the square of the map factor in the case of
  !! conformal projections
  !!
  !! @author George Gayno, Mark Iredell
  !! @date Jan 2015
  SUBROUTINE gdswzd_1d_array(IGDTNUM,IGDTMPL,IGDTLEN,IOPT,NPTS,FILL, &
       XPTS,YPTS,RLON,RLAT,NRET, &
       CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
    INTEGER,        INTENT(IN   ) :: IGDTNUM, IGDTLEN
    INTEGER,        INTENT(IN   ) :: IGDTMPL(IGDTLEN)
    INTEGER,        INTENT(IN   ) :: IOPT, NPTS
    INTEGER,        INTENT(  OUT) :: NRET
    !
    REAL,           INTENT(IN   ) :: FILL
    REAL,           INTENT(INOUT) :: RLON(NPTS),RLAT(NPTS)
    REAL,           INTENT(INOUT) :: XPTS(NPTS),YPTS(NPTS)
    REAL, OPTIONAL, INTENT(  OUT) :: CROT(NPTS),SROT(NPTS)
    REAL, OPTIONAL, INTENT(  OUT) :: XLON(NPTS),XLAT(NPTS)
    REAL, OPTIONAL, INTENT(  OUT) :: YLON(NPTS),YLAT(NPTS),AREA(NPTS)
 
    type(grib2_descriptor) :: desc
    class(ip_grid), allocatable :: grid
 
    desc = init_descriptor(igdtnum, igdtlen, igdtmpl)
    call init_grid(grid, desc)
    
    call gdswzd_grid(grid,iopt,npts,fill, &
         xpts,ypts,rlon,rlat,nret, &
         crot,srot,xlon,xlat,ylon,ylat,area)
  
  SUBROUTINE gdswzd_1d_array(IGDTNUM,IGDTMPL,IGDTLEN,IOPT,NPTS,FILL, & …
  END SUBROUTINE gdswzd_1d_array
 
  !> Decodes the grib grid description section and
  !! returns one of the following (for 1-d arrays):
  !! - 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.
  !!
  !! 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, map jacobians and
  !! grid box areas may be returned.
  !!
  !! To compute the vector rotations, the optional arguments 'srot' and 'crot'
  !! must be present. To compute the map jacobians, the
  !! optional arguments 'xlon', 'xlat', 'ylon', 'ylat' must be present.
  !!
  !! To compute the grid box areas, the optional argument
  !! 'area' must be present.
  !!
  !! 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)=203 E-staggered Rotated Equidistant Cylindrical
  !! - kgds(1)=205 B-staggered Rotated Equidistant Cylindrical
  !!
  !! @param[in] kgds GDS parameters as decoded by w3fi63.  
  !! @param[in] iopt Option flag.
  !! - 0 Earth coords of all the grid points.
  !! - 1 Earth coords of selected grid coords.
  !! - -1 Grid coords of selected earth coords
  !!
  !! @param[in] npts Maximum number of coordinates.
  !! @param[in] fill Fill value to set invalid output data.
  !!      Must be impossible value; suggested value: -9999.
  !! @param[inout] xpts Grid x point coordinates.
  !! @param[inout] ypts Grid y point coordinates.
  !! @param[inout] rlon Earth longitudes in degrees E.
  !! (Acceptable range: -360. to 360.)
  !!
  !! @param[inout] rlat Earth latitudes in degrees N.
  !! (Acceptable range: -90. to 90.)
  !!
  !! @param[out] nret Number of valid points computed.
  !! @param[out] crot Clockwise vector rotation cosines.
  !! @param[out] srot Clockwise vector rotation sines.
  !! ugrid=crot*uearth-srot*vearth;
  !! vgrid=srot*uearth+crot*vearth)
  !!
  !! @param[out] xlon dx/dlat in 1/degrees
  !! @param[out] xlat dy/dlon in 1/degrees
  !! @param[out] ylon dy/dlon in 1/degrees
  !! @param[out] ylat dy/dlat in 1/degrees
  !! @param[out] area Area weights in m^2.
  !! Proportional to the square of the map factor in the case of
  !! conformal projections
  !!
  !! @author George Gayno, Mark Iredell
  !! @date April 1996
  SUBROUTINE gdswzd_grib1(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, &
       CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
    INTEGER,        INTENT(IN   ) :: IOPT, KGDS(200), NPTS
    INTEGER,        INTENT(  OUT) :: NRET
    !
    REAL,           INTENT(IN   ) :: FILL
    REAL,           INTENT(INOUT) :: RLON(NPTS),RLAT(NPTS)
    REAL,           INTENT(INOUT) :: XPTS(NPTS),YPTS(NPTS)
    REAL, OPTIONAL, INTENT(  OUT) :: CROT(NPTS),SROT(NPTS)
    REAL, OPTIONAL, INTENT(  OUT) :: XLON(NPTS),XLAT(NPTS)
    REAL, OPTIONAL, INTENT(  OUT) :: YLON(NPTS),YLAT(NPTS),AREA(NPTS)
 
    
    type(grib1_descriptor) :: desc
    class(ip_grid), allocatable :: grid
 
    desc = init_descriptor(kgds)
    call init_grid(grid, desc)
    
    call gdswzd_grid(grid,iopt,npts,fill, &
         xpts,ypts,rlon,rlat,nret, &
         crot,srot,xlon,xlat,ylon,ylat,area)
 
  SUBROUTINE gdswzd_grib1(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & …
  END SUBROUTINE gdswzd_grib1
 
  
  !> Decodes the grib grid description section and returns
  !! one of the following (for 2-d arrays):
  !! - 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.
  !!
  !! 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, map jacobians and
  !! grid box areas may be returned.
  !!
  !! To compute the vector rotations, the optional arguments 'srot' and 'crot'
  !! must be present. To compute the map jacobians, the
  !! optional arguments 'xlon', 'xlat', 'ylon', 'ylat' must be present.
  !!
  !! To compute the grid box areas, the optional argument
  !! 'area' must be present.
  !!
  !! 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)=203 E-staggered Rotated Equidistant Cylindrical
  !! - kgds(1)=205 B-staggered Rotated Equidistant Cylindrical
  !!
  !! @param[in] kgds GDS parameters as decoded by w3fi63.  
  !! @param[in] iopt Option flag.
  !! - 0 Earth coords of all the grid points.
  !! - 1 Earth coords of selected grid coords.
  !! - -1 Grid coords of selected earth coords
  !!
  !! @param[in] npts Maximum number of coordinates.
  !! @param[in] fill Fill value to set invalid output data.
  !!      Must be impossible value; suggested value: -9999.
  !! @param[inout] xpts Grid x point coordinates.
  !! @param[inout] ypts Grid y point coordinates.
  !! @param[inout] rlon Earth longitudes in degrees E.
  !! (Acceptable range: -360. to 360.)
  !!
  !! @param[inout] rlat Earth latitudes in degrees N.
  !! (Acceptable range: -90. to 90.)
  !!
  !! @param[out] nret Number of valid points computed.
  !! @param[out] crot Clockwise vector rotation cosines.
  !! @param[out] srot Clockwise vector rotation sines.
  !! ugrid=crot*uearth-srot*vearth;
  !! vgrid=srot*uearth+crot*vearth)
  !!
  !! @param[out] xlon dx/dlat in 1/degrees
  !! @param[out] xlat dy/dlon in 1/degrees
  !! @param[out] ylon dy/dlon in 1/degrees
  !! @param[out] ylat dy/dlat in 1/degrees
  !! @param[out] area Area weights in m^2.
  !! Proportional to the square of the map factor in the case of
  !! conformal projections
  !!
  !! @author George Gayno, Mark Iredell
  !! @date April 1996
  SUBROUTINE gdswzd_2d_array_grib1(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, &
       CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
   
    !$$$
    INTEGER,        INTENT(IN   ) :: IOPT, KGDS(200), NPTS
    INTEGER,        INTENT(  OUT) :: NRET
    !
    REAL,           INTENT(IN   ) :: FILL
    REAL,           INTENT(INOUT) :: RLON(:,:),RLAT(:,:)
    REAL,           INTENT(INOUT) :: XPTS(:,:),YPTS(:,:)
    REAL, OPTIONAL, INTENT(  OUT) :: CROT(:,:),SROT(:,:)
    REAL, OPTIONAL, INTENT(  OUT) :: XLON(:,:),XLAT(:,:)
    REAL, OPTIONAL, INTENT(  OUT) :: YLON(:,:),YLAT(:,:),AREA(:,:)
 
    
    type(grib1_descriptor) :: desc
    class(ip_grid), allocatable :: grid
 
    desc = init_descriptor(kgds)
    call init_grid(grid, desc)
    
    call gdswzd_grid(grid,iopt,npts,fill, &
         xpts,ypts,rlon,rlat,nret, &
         crot,srot,xlon,xlat,ylon,ylat,area)
 
  SUBROUTINE gdswzd_2d_array_grib1(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & …
  END SUBROUTINE gdswzd_2d_array_grib1
 
 
 
END MODULE gdswzd_mod