ww3_grib.F90

Go to the documentation of this file.

 
!
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
!
PROGRAM w3grib
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |            A. Chawla              |
  !/                  |           J.-H. Alves             |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         22-Mar-2021 |
  !/                  +-----------------------------------+
  !/
  !/    01-Nov-1999 : Final FORTRAN 77        ( version 1.18 + error fix )
  !/    24-Jan-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
  !/    25-Jan-2001 : Flat grid error exit added          ( version 2.06 )
  !/    29-Apr-2002 : Adding output fields 17-18.         ( version 2.20 )
  !/    08-May-2002 : Replace XLF switch with NCEP1.      ( version 2.21 )
  !/    13-Nov-2002 : Add stress vector.                  ( version 3.00 )
  !/    24-Dec-2004 : Multiple grid version.              ( version 3.06 )
  !/    20-Jul-2005 : Additional output parameters.       ( version 3.07 )
  !/    11-Apr-2007 : Additional output parameters.       ( version 3.11 )
  !/    18-May-2007 : Update GRIB1 for partitioning.      ( version 3.11 )
  !/    16-Jul-2007 : Adding GRIB2 capability.            ( version 3.11 )
  !/                  (A. Chawla)
  !/    01-Aug-2007 : Update FLGRIB for GRIB2.            ( version 3.11 )
  !/    29-May-2009 : Preparing distribution version.     ( version 3.14 )
  !/    30-Oct-2009 : Implement run-time grid selection.  ( version 3.14 )
  !/                  (W. E. Rogers & T. J. Campbell, NRL)
  !/    05-Oct-2011 : Updating to the 53 output parameter ( version 4.05 )
  !/                  (Arun Chawla)
  !/    01-Mar-2013 : Adding double-index output fields   ( version 4.11 )
  !/                  (J-Henrique Alves)
  !/    01-Dec-2016 : Adding lambert conformal grid       ( version 6.01 )
  !/                  (J.H. Alves)
  !/    26-Jul-2018 : Adding polar stereographic grid     ( version 6.05 )
  !/                  (J.H. Alves)
  !/    22-Mar-2021 : New coupling fields output          ( version 7.13 )
  !/    09-Jun-2021 : remove grib1 support (NCEP1)        ( version 7.14 )
  !/
  !/    Copyright 2009 National Weather Service (NWS),
  !/       National Oceanic and Atmospheric Administration.  All rights
  !/       reserved.  WAVEWATCH III is a trademark of the NWS.
  !/       No unauthorized use without permission.
  !/
  !  1. Purpose :
  !
  !     Post-processing of grid output.
  !
  !  2. Method :
  !
  !     Data is read from the grid output file out_grd.ww3 (raw data)
  !     and from the file ww3_grib.inp ( NDSI, output requests ).
  !     Model definition and raw data files are read using WAVEWATCH III
  !     subroutines.
  !     GRIB packing is performed using NCEP's W3 library (not supplied).
  !
  !     When adding new parameters to GRIB packing, keep in mind that
  !     packing is done differently for scalar and vector quantities
  !
  !  3. Parameters :
  !
  !  4. Subroutines used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      W3NMOD    Subr. W3GDATMD Set number of model.
  !      W3SETG    Subr.   Id.    Point to selected model.
  !      W3NDAT    Subr. W3WDATMD Set number of model for wave data.
  !      W3SETW    Subr.   Id.    Point to selected model for wave data.
  !      W3NAUX    Subr. W3ADATMD Set number of model for aux data.
  !      W3SETA    Subr.   Id.    Point to selected model for aux data.
  !      ITRACE    Subr. W3SERVMD Subroutine tracing initialization.
  !      STRACE    Subr.   Id.    Subroutine tracing.
  !      NEXTLN    Subr.   Id.    Get next line from input filw
  !      EXTCDE    Subr.   Id.    Abort program as graceful as possible.
  !      STME21    Subr. W3TIMEMD Convert time to string.
  !      TICK21    Subr.   Id.    Advance time.
  !      DSEC21    Func.   Id.    Difference between times.
  !      W3IOGR    Subr. W3IOGRMD Reading/writing model definition file.
  !      W3IOGO    Subr. W3IOGOMD Reading/writing raw gridded data file.
  !      W3READFLGRD Subr. W3IOGOMD Reading output fields flags.
  !      W3EXGB    Subr. Internal Execute grib output.
  !      BAOPEN    Subr.          NCEP library routine.
  !      BAOPENW    Subr.          NCEP library routine.
  !     ----------------------------------------------------------------
  !
  !  5. Called by :
  !
  !     None, stand-alone program.
  !
  !  6. Error messages :
  !
  !     Checks on input, checks in W3IOxx.
  !
  !  7. Remarks :
  !
  !  8. Structure :
  !
  !     See source code.
  !
  !  9. Switches :
  !
  !       !/S     Enable subroutine tracing.
  !
  !       !/NCO   NCEP NCO modifications for operational implementation.
  !
  !       !/NOGRB No GRIB package included.
  !       !/NCEP2 NCEP IBM links to GRIB2 packing routines.
  !
  ! 10. Source code :
  !
  !/ ------------------------------------------------------------------- /
  USE constants
  !
  !     USE W3GDATMD, ONLY: W3NMOD, W3SETG
  USE w3wdatmd, ONLY: w3ndat, w3setw
  !      USE W3ADATMD, ONLY: W3NAUX, W3SETA
  USE w3odatmd, ONLY: w3nout, w3seto
  USE w3iogrmd, ONLY: w3iogr
  USE w3iogomd, ONLY: w3readflgrd, w3iogo
  USE w3servmd, ONLY : itrace, nextln, extcde
#ifdef W3_S
  USE w3servmd, ONLY : strace
#endif
  USE w3timemd, ONLY: stme21, tick21, dsec21
  !
  USE w3gdatmd
  USE w3wdatmd, ONLY: time, wlv, ice, ust, ustdir, rhoair
  USE w3adatmd
  USE w3odatmd, ONLY: ndse, ndst, ndso, nogrp, ngrpp, idout, undef,&
       flogrd, fnmpre, noswll, noge, flogd
  !
  IMPLICIT NONE
  !/
  !/ ------------------------------------------------------------------- /
  !/ Local variables
  !/
  INTEGER                 :: ndsi, ndsm, ndsog, ndsdat, ndstrc,   &
       ntrace, ierr, iotest, i,j,k, ifi,ifj,&
       isea, ix, iy, tout(2), nout, tdum(2),&
       ftime(2), cid, pid, gid, gds, iout,  &
       gdtn
  INTEGER, ALLOCATABLE    :: ifia(:),ifja(:)
#ifdef W3_NOGRB
  INTEGER                 :: kpds(1), kgds(1)
#endif
  ! GRIB2 specific variables
#ifdef W3_NCEP2
  INTEGER                 :: kpds(200), kgds(200), idrs(200)
  INTEGER                 :: listsec0(3), listsec1(13),igds(5)
  INTEGER                 :: ideflist, idefnum, kpdsnum, numcoord
  INTEGER                 :: ibmp, lcgrib, lengrib, idrsnum
  REAL                    :: coordlist, xn
  CHARACTER(LEN=1), ALLOCATABLE  :: cgrib(:)
  INTEGER                 :: latan1, lonv, scnmod, latin1, &
       latin2, latsp, lonsp
  REAL                    :: dsx, dsy
  REAL                    :: yn, x0n, y0n
#endif
#ifdef W3_S
  INTEGER, SAVE           :: ient = 0
#endif
  REAL                    :: dtreq, dtest, rftime
  LOGICAL                 :: flreq(nogrp,ngrpp), flgrib(nogrp,ngrpp)
  CHARACTER               :: comstr*1, idtime*23, iddday*11
  CHARACTER(LEN=80)       :: linein
  CHARACTER(LEN=8)        :: words(5)
  INTEGER                 :: gen_pro
 
  !/
  !/ ------------------------------------------------------------------- /
  !/
#ifdef W3_NCO
  !     CALL W3TAGB('WAVEGRIB',1998,0007,0050,'NP21   ')
#endif
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 1.  IO set-up.
  !
  CALL w3nmod ( 1, 6, 6 )
  CALL w3setg ( 1, 6, 6 )
  CALL w3ndat (    6, 6 )
  CALL w3setw ( 1, 6, 6 )
  CALL w3naux (    6, 6 )
  CALL w3seta ( 1, 6, 6 )
  CALL w3nout (    6, 6 )
  CALL w3seto ( 1, 6, 6 )
  !
  ndsi   = 10
  ndsm   = 20
  ndsog  = 20
  ndsdat = 50
  !
  ndstrc =  6
  ntrace = 10
  !
#ifdef W3_NCO
  !
  ! Redo according to NCO
  !
  ndsi   = 11
  ndso   =  6
  ndse   = ndso
  ndst   = ndso
  ndsm   = 12
  ndsog  = 13
  ndsdat = 51
  ndstrc = ndso
#endif
  !
  WRITE (ndso,900)
  !
  CALL itrace ( ndstrc, ntrace )
#ifdef W3_S
  CALL strace (ient, 'W3GRIB')
#endif
  !
  OPEN (ndsi,file='ww3_grib.inp',status='OLD',err=800,iostat=ierr)
  READ (ndsi,'(A)',END=801,ERR=802) comstr
  IF (comstr.EQ.' ') comstr = '$'
  WRITE (ndso,901) comstr
  !
#ifdef W3_NOGRB
  WRITE (ndse,902)
#endif
#ifdef W3_NCEP2
  CALL baopenw (ndsdat,'gribfile',ierr)
#endif
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 2.  Read model definition file.
  !
  CALL w3iogr ( 'READ', ndsm )
  WRITE (ndso,920) gname
  !
  IF ( .NOT. flagll ) GOTO 810
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 3.  Read requests from input file.
  !     Output times
  !
  CALL nextln ( comstr , ndsi , ndse )
  READ (ndsi,'(A)') linein
  WRITE(ndso,*)' LINEIN:  ',linein
  READ(linein,*,iostat=ierr) words
  WRITE (ndso,*) words
  READ(words( 1 ), * ) tout(1)
  READ(words( 2 ), * ) tout(2)
  READ(words( 3 ), * ) dtreq
  READ(words( 4 ), * ) nout
  IF (words(5) .NE. '0' .AND. words(5) .NE. '1') THEN
    gen_pro=-99999
  ELSE
    READ(words( 5 ), * ) gen_pro
  ENDIF
  WRITE(ndso,*) 'GEN_PRO  ',gen_pro
  dtreq  = max( 0. , dtreq )
  IF ( dtreq.EQ.0 ) nout = 1
  nout   = max( 1 , nout )
  !
  CALL stme21 ( tout , idtime )
  WRITE (ndso,940) idtime
  !
  tdum(1) = 0
  tdum(2) = 0
  CALL tick21 ( tdum , dtreq )
  CALL stme21 ( tdum , idtime )
  IF ( dtreq .GE. 86400. ) THEN
    WRITE (iddday,'(I10,1X)') int(dtreq/86400.)
  ELSE
    iddday = '           '
  END IF
  idtime(1:11) = iddday
  idtime(21:23) = '   '
  WRITE (ndso,941) idtime, nout
  !
  ! ... Initialize FLGRD array
  !
  flreq(:,:)=.false.
  !
  ! ... Call to interface for reading flags or namelists
  !
  CALL w3readflgrd ( ndsi, ndso, 9, ndse, comstr, flogd, flreq,   &
       1, 1, ierr )
  !
  ! Inform user of parameters that were requested but failed to make the
  ! grade, as they are not available for grib encoding, or are not
  ! included presently
  !
  WRITE (ndso,944)
  ! Reset flags for variables not yet implemented in grib output
  ! interface
  !
  !
  ifi = 3 ! Entire group Frequency-dependent parameters
  DO ifj = 1,noge(ifi)
    IF ( flreq(ifi,ifj) ) THEN
      WRITE (ndso,946) idout(ifi,ifj),                           &
           '*** NOT YET CODED INTO WW3_GRIB ***'
      flreq(ifi,ifj) = .false.
    END IF
  END DO
  !
  ifi = 5 ! Atm-waves layer, all except for friction velocity
  DO ifj = 2,10
    IF ( flreq(ifi,ifj) ) THEN
      WRITE (ndso,946) idout(ifi,ifj),                           &
           '*** NOT YET CODED INTO WW3_GRIB ***'
      flreq(ifi,ifj) = .false.
    END IF
  END DO
  DO ifi = 6,8 ! Entire groups wave-ocean interaction, wave-bottom
    ! layer and spectrum parameters
    DO ifj = 1,noge(ifi)
      IF ( flreq(ifi,ifj) ) THEN
        WRITE (ndso,946) idout(ifi,ifj),                           &
             '*** NOT YET CODED INTO WW3_GRIB ***'
        flreq(ifi,ifj) = .false.
      END IF
    END DO
  END DO
  IF ( flreq(9,5) ) THEN ! CFL number for K advection
    WRITE (ndso,946) idout(9,5),'*** NOT YET CODED INTO WW3_GRIB ***'
    flreq(9,5) = .false.
  END IF
  ifi = 10 ! User defined parameters
  DO ifj = 1,noge(ifi)
    IF ( flreq(ifi,ifj) ) THEN
      WRITE (ndso,946) idout(ifi,ifj),                           &
           '*** NOT YET CODED INTO WW3_GRIB ***'
      flreq(ifi,ifj) = .false.
    END IF
  END DO
  !
  ! Compatibility with NCEP operational codes, same effect as old FLGRIB
  !  lists variables that have no code for variable names (not 100%
  !  correct in old codes... )
  !
  ! Chage this as parameters become available in grib2 tables
  !
  ALLOCATE ( ifia(13), ifja(13) )
 
  ifia = (/ 1, 2, 2, 4, 4, 4, 4, 4, 5, 9, 9, 9, 9 /)
  ifja = (/ 4, 2, 8, 3, 5, 6, 7, 8, 1, 1, 2, 3, 4 /)
  DO i = 1, 13
    IF ( flreq(ifia(i),ifja(i)) ) THEN
      flreq(ifia(i),ifja(i)) = .false.
      WRITE(ndso,946) idout(ifia(i),ifja(i)),                      &
           '*** EXCLUDED FROM GRIB OUTPUT ***'
    END IF
  END DO
  !
  ! Write to stdout parameters that have successfully been requested
  !
  WRITE (ndso,945)
  DO i=1, nogrp
    DO j=1, ngrpp
      IF ( flreq(i,j) ) WRITE (ndso,931) idout(i,j)
    END DO
  END DO
  !
  !
  !
  ! ... GRIB specific parameters
  !
  CALL nextln ( comstr , ndsi , ndse )
  READ (ndsi,*,END=801,ERR=802) FTIME, CID, PID, GID, GDS, gdtn
  !
  ! Check if grid type is curvilinear, and only go on if Lambert conformal
  ! or PolarStereo
  !
  IF ( gtype .EQ. clgtype ) THEN
#ifdef W3_NCEP2
    ! Allowing code to work with Lambert conformal grids
    IF ( gdtn .NE. 30 .AND. gdtn .NE. 20 ) THEN
#endif
      WRITE(ndse,*)'PROGRAM W3GRIB: CURVILINEAR GRID SUPPORT '// &
           'FOR GRIB OUTPUT IS NOT YET IMPLEMENTED. NOW STOPPING'
      CALL extcde ( 1 )
#ifdef W3_NCEP2
    ENDIF
#endif
  END IF
  !
  !
  ! Coded up to now only for Lamber conformal grids (GDTN=30) or
  ! PolarStereo (GDTN=20). For regular grids use GDTN=0
  !
#ifdef W3_NCEP2
  IF ( gdtn .EQ. 30 ) THEN
    ! This is a Lambert conformal grid, read projection parameters
    CALL nextln ( comstr , ndsi , ndse )
    READ (ndsi,*,END=801,ERR=802) LATAN1, LONV, DSX, DSY,          &
         scnmod, latin1, latin2, latsp, lonsp
  ELSEIF ( gdtn .EQ. 20 ) THEN
    CALL nextln ( comstr , ndsi , ndse )
    READ (ndsi,*,END=801,ERR=802) LATAN1, LONV, DSX, DSY,   &
         scnmod
#endif
 
#ifdef W3_NCEP2
  ENDIF
#endif
  !
  CALL stme21 ( ftime , idtime )
  WRITE (ndso,948) idtime, cid, pid, gid, gds
  !
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 4.  Read general data and first fields from file
  ! 4.a Read file.
  !
  CALL w3iogo ( 'READ', ndsog, iotest )
  !
  ! 4.b Output fields in file
  !
  !
  WRITE (ndso,930)
  DO i=1, nogrp
    DO j=1, ngrpp
      IF ( flogrd(i,j) ) WRITE (ndso,931) idout(i,j)
    END DO
  END DO
  !
#ifdef W3_NCEP2
  !
  IF ( gdtn .EQ. 0 ) THEN
    !
#endif
    ! 4.c Flip MAPSF for REGULAR/RECTILINEAR grids
    !
    DO isea=1, nsea
      ix            = mapsf(isea,1)
      iy            = mapsf(isea,2)
      mapsf(isea,2) = ny + 1 - iy
      mapsf(isea,3) = iy +( ix-1)*ny
    END DO
#ifdef W3_NCEP2
    !
  ENDIF
#endif
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 5. Set grib encoding parameter Sections
  !
  ! ... Initialize KPDS and KGDS (for NCEP2)
  !
  kpds     =   0
  kgds     =   0
  !
  ! ... Set GRIB2 packing arrays
  !
#ifdef W3_NCEP2
  lcgrib = 4*nx*ny
  ALLOCATE(cgrib(lcgrib))
#endif
  !
  ! ... Set GRIB2 Indicator Section
  !      ( 1)  Discipline-GRIB Master Table Number (see Code Table 0.0)
  !            0 = Metereological; 10 = Oceanographic
  !      ( 2)  GRIB Edition Number
  !      ( 3)
#ifdef W3_NCEP2
  listsec0      =   0
  listsec0(1)   =   10
  listsec0(2)   =   2
#endif
  !
  ! ... Set GRIB2 Identification Section
  !       ( 1) ID OF CENTER
  !       ( 2) ID OF SUB-CENTER
  !       ( 3) GRIB Master Tables Version Number (Code Table 1.0)
  !       ( 4) GRIB Local Tables Version Number (Code Table 1.0)
  !       ( 5) Significance of Reference Time (Code Table 1.2)
  !     * ( 6) YEAR (4 digits)
  !     * ( 7) MONTH OF YEAR
  !     * ( 8) DAY OF MONTH
  !     * ( 9) HOUR OF DAY
  !       (10) MINUTE OF HOUR
  !       (11) SECOND OF MINUTE
  !       (12) Production status of data (Code Table 1.3)
  !       (13) Type of processed data (Code Table 1.4)
  !
#ifdef W3_NCEP2
  listsec1     =  0
  listsec1(1)  = cid
  listsec1(3)  = 2
  listsec1(4)  = 1
  listsec1(5)  = 1
  listsec1(13) = 1
#endif
  !
  ! ... Set GRIB2 IGDS elements
  !       ( 1) Source of grid definition (Code Table 3.0)
  !       ( 2) Number of grid points
  !       ( 3) Number of octets needed for each additional grid points definition
  !       ( 4) Interpretation of list for optional points definition (Code Table 3.11)
  !       ( 5) Grid definition template number (Code Table 3.1)
  !
#ifdef W3_NCEP2
  igds     = 0 ! Defined in code
  igds(2)  = nx*ny
  idefnum  = 0
  ideflist = 0
  igds(5)=gdtn
  IF ( gdtn .EQ. 30 .AND. gtype .EQ. clgtype ) THEN
    idefnum  = 1
    WRITE (ndso,1011) 'LAMBERTCONF'
  ELSEIF ( gdtn .EQ. 20 .AND. gtype .EQ. clgtype ) THEN
    WRITE (ndso,1011) 'POLARSTEREO'
  ELSEIF ( gdtn .EQ. 0 ) THEN
    WRITE (ndso,1011) 'LLRECTILINEAR'
  ELSE
    WRITE(ndse,*)'PROGRAM WAVEGRIB2: SUPPORT FOR CHOSEN '// &
         'GRIB2 GRID DEFINITION TEMPLATE NOT YET IMPLEMENTED'
    CALL extcde ( 2 )
  ENDIF
#endif
  !
  ! ... Set GRIB2 KGDS elements
  !
  ! General parameters for all grids
  !       ( 1) Coordinate system (6 = spherical coordinate system with radius of 6,371,229 m)
  !       ( 2)
  !       ( 3)
  !       ( 4)
  !       ( 5)
  !       ( 6)
  !       ( 7)
  !       ( 8) Number of points along parallel
  !       ( 9) Number of points along meridian
#ifdef W3_NCEP2
  kgds( 1) = 6
  kgds( 8) = nx
  kgds( 9) = ny
#endif
  !
#ifdef W3_NCEP2
  IF ( gdtn .EQ. 30 ) THEN
#endif
    !
    ! Lambert Conformal grid
    !       (10) Latitude of first grid point
    !       (11) Longitude of first grid point
    !       (12) Resolution and component flags
    !       (13) Latitude where DX and DY are specified
    !       (14) Longitude of orientation
    !       (15) Increment of longitude
    !       (16) Increment of latitude
    !       (17) Projection center flag
    !       (18) Scanning mode
    !       (19) First latitude of secant cone
    !       (20) Second latitude of secant cone
    !       (21) Latitude of southern pole
    !       (22) Longitude of southern pole
    !
#ifdef W3_NCEP2
    x0 = mod(xgrd(1,1) + 360.,360.)
    xn = mod(xgrd(ny,nx) + 360., 360.)
    x0n = mod(xgrd(ny,1) + 360., 360.)
    kgds(11)=nint(1000000.*x0)
    y0 = ygrd(1,1)
    yn = ygrd(ny,nx)
    y0n = ygrd(ny,1)
    kgds(10)=nint(1000000.*y0)
    kgds(12)=0
    kgds(13)=dble(1000000.*latan1)
    kgds(14)=dble(1000000.*lonv)
    kgds(15)=nint(1000000*dsx)
    kgds(16)=nint(1000000*dsy)
    kgds(17)=0
    kgds(18)=scnmod
    kgds(19)=dble(1000000.*latin1)
    kgds(20)=dble(1000000.*latin2)
    kgds(21)=dble(1000000.*latsp)
    kgds(22)=dble(1000000.*lonsp)
#endif
    !
#ifdef W3_NCEP2
  ELSEIF (gdtn .EQ. 20 ) THEN
#endif
    !
    ! PolarStereo grid
    !       (10) Latitude of first grid point
    !       (11) Longitude of first grid point
    !       (12) Res and component flags
    !       (13) Latitude where DX and DY are specified
    !       (14) Longitude of orientation
    !       (15) Increment of longitude
    !       (16) Increment of latitude
    !       (17) Projection center flag
    !       (18) Scanning mode
    !
    !  Projection for PolarStereo grid was changed from
    !  KGDS( 1) = 6 to KGDS( 1) = 5 (Earth assumed represented by WGS84 -
    !  Octet No 15 Table 3.2)
#ifdef W3_NCEP2
    kgds( 1) = 5
    x0 = mod(xgrd(1,1) + 360.,360.)
    xn = mod(xgrd(ny,nx) + 360., 360.)
    x0n = mod(xgrd(ny,1) + 360., 360.)
    kgds(11)=nint(1000000.*x0)
    y0 = ygrd(1,1)
    yn = ygrd(ny,nx)
    y0n = ygrd(ny,1)
    kgds(10)=nint(1000000.*y0)
    kgds(12)=0
    kgds(13)=dble(1000000.*latan1)
    kgds(14)=dble(1000000.*lonv)
    kgds(15)=nint(1000000*dsx)
    kgds(16)=nint(1000000*dsy)
    kgds(17)=0
    kgds(18)=scnmod
#endif
    !
#ifdef W3_NCEP2
  ELSEIF (gdtn .EQ. 0 ) THEN
#endif
    !
    ! Lat Lon rectilinear grid
    !       (10)
    !       (11)
    !       (12) Latitude of first grid point
    !       (13) Longitude of first grid point
    !       (14) Res and component flags
    !       (15) Latitude of last grid point
    !       (16) Longitude of last grid point
    !       (17) Increment of longitude
    !       (18) Increment of latitude
    !       (19) Scanning mode
    !
#ifdef W3_NCEP2
    kgds(12) = nint(1000000.*(y0+(real(ny-1)*sy)))
    x0 = mod(x0 + 360.,360.)
    kgds(13) = nint(1000000.*x0)
    kgds(14) = 48
    kgds(15) = nint(1000000.*y0)
    xn = mod(x0+real(nx-1)*sx + 360., 360.)
    kgds(16) = nint(1000000.*xn)
    kgds(17) = nint(1000000.*sx)
    kgds(18) = nint(1000000.*sy)
  ENDIF
#endif
  !
  ! ... Set GRIB2 PDS elements
  !     KPDSNUM   (0 indicates forecast at a horizontal level)
  !     ( 1) Parameter category (Code Table 4.1)
  !            For oceanographic products -- 0 = waves; 1 = currents; 2 = ice
  !            For atmospheric products   -- 2 = momentum
  !     ( 2) Parameter number (Code Table 4.2)
  !     ( 3) Generating process (Code Table 4.3)
  !     ( 4) Background generating process identifier (center specific)
  !     ( 5) Process or model number
  !     ( 6) Hours of observational data cutoff after reference time
  !     ( 7) Minutes of observational data cutoff after reference time
  !     ( 8) Indicator of forecast time unit (Code Table 4.4)
  !     ( 9) Time range
  !     (10) Type of level (Code Table 4.5) 1st level
  !     (11) Scaled factor of (10)
  !     (12) Scaled value of (10)
  !     (13) Type of level (Code Table 4.5) 2nd level
  !     (14) Scaled factor of (13)
  !     (15) Scaled value of (13)
  !
  !
  !      KPDS(3)=4 ensemble forecast:ww3_grib.inp has gen_pro set to 1
  !             =2 deterministic forecast: ww3_grib.inp gen_pro set to 0
  !             =2 legacy :with no gen_pro set in ww3_grib.inp
  !              (in the case of legacy the params revert back to old names)
#ifdef W3_NCEP2
  kpdsnum = 0
  if ( gen_pro.eq.1 ) then
    kpds( 3) = 4
  else
    kpds(3)=2
  endif
  kpds( 4) = 0
  kpds( 5) = pid
  kpds( 8) = 1
  kpds(10) = 1
  kpds(12) = 1
  kpds(13) = 255
#endif
  !
  ! ... Set GRIB2 vertical layer information
  !
#ifdef W3_NCEP2
  numcoord  = 0
  coordlist = 0.0
#endif
  !
  ! ... Set GRIB2 bitmap information
  !      0 Bitmap is provided
  !
#ifdef W3_NCEP2
  ibmp = gds
#endif
  !
  ! ... Set GRIB2 Data Representation Template Number (Code Table 5.0)
  !
#ifdef W3_NCEP2
  idrsnum = 40 !jpeg2000 *** SEGFAULTS in some linux
#endif
  !                            clusters with Intel compiler ***
#ifdef W3_NCEP2
  !IDRSNUM = 0 !simple packing
  !IDRSNUM = 41 !png packing
  !IDRSNUM = 2 !Complex Packing (Grid Point Data)
#endif
  !
  ! ... Set GRIB2 IDRS elements
  !     ( 1) Reference value (R) (IEEE 32-bit floating-point value)
  !     ( 2) Binary Scale Factor (E)
  !     ( 3) Decimal Scale Factor (D)
  !     ( 4) Number of bits used for each packed value
  !     ( 5) Type of original field values (Code Table 5.1)
  !
#ifdef W3_NCEP2
  idrs    = 0
  idrs(3) = 2
#endif
  !
#ifdef W3_T
  WRITE (ndst,9050) kpds
  WRITE (ndst,9051) kgds
#endif
  !
  !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  ! 6.  Time management.
  !
  iout   = 0
  WRITE (ndso,970)
  !
  DO
    dtest  = dsec21( time , tout )
    IF ( dtest .GT. 0. ) THEN
      CALL w3iogo ( 'READ', ndsog, iotest )
      IF ( iotest .EQ. -1 ) THEN
        WRITE (ndso,942)
        GOTO 888
      END IF
      cycle
    END IF
    IF ( dtest .LT. 0. ) THEN
      CALL tick21 ( tout , dtreq )
      cycle
    END IF
    !
    iout   = iout + 1
    CALL stme21 ( tout , idtime )
    !
    rftime = dsec21( ftime , time ) / 3600.
    IF ( rftime .LT. 0. ) THEN
#ifdef W3_NCEP2
      listsec1( 6) = time(1)/10000
      listsec1( 7) = mod(time(1),10000) / 100
      listsec1( 8) = mod(time(1),100)
      listsec1( 9) = time(2) / 10000
      kpds( 9)     = 0
#endif
      WRITE (ndso,972) idtime
    ELSE
#ifdef W3_NCEP2
      listsec1( 6) = ftime(1)/10000
      listsec1( 7) = mod(ftime(1),10000) / 100
      listsec1( 8) = mod(ftime(1),100)
      listsec1( 9) = ftime(2) / 10000
      kpds( 9)     = nint(rftime)
#endif
      WRITE (ndso,971) idtime, nint(rftime)
    END IF
    !
    CALL w3exgb ( nx, ny, nsea )
    CALL tick21 ( tout , dtreq )
    IF ( iout .GE. nout ) EXIT
  END DO
  !
  GOTO 888
  !
  ! Escape locations read errors :
  !
800 CONTINUE
  WRITE (ndse,1000) ierr
  CALL extcde ( 3 )
  !
801 CONTINUE
  WRITE (ndse,1001)
  CALL extcde ( 4 )
  !
802 CONTINUE
  WRITE (ndse,1002) ierr
  CALL extcde ( 5 )
  !
810 CONTINUE
  IF ( .NOT. flagll ) THEN
    WRITE (ndse,1010)
    CALL extcde ( 10 )
  END IF
  !
888 CONTINUE
  WRITE (ndso,999)
  !
#ifdef W3_NCO
  !     CALL W3TAGE('WAVEGRIB')
#endif
  !
  ! Formats
  !
900 FORMAT (/15x,'   *** WAVEWATCH III GRIB output postp. ***   '/ &
       15x,'=============================================='/)
901 FORMAT ( '  Comment character is ''',a,''''/)
902 FORMAT (/'  *** WARNING : NO GRIB PACKAGE LINKED ***'/)
  !
920 FORMAT ( '  Grid name : ',a/)
  !
930 FORMAT ( '  Fields in file : '/                                 &
       ' --------------------------')
931 FORMAT ( '      ',a)
  !
940 FORMAT (/'  Output time data : '/                               &
       ' -----------------------------------------------------'/ &
       '      First time         : ',a)
941 FORMAT ( '      Interval           : ',a/                       &
       '      Number of requests : ',i4)
942 FORMAT (/'      End of file reached '/)
  !
944 FORMAT (/'  Requested output fields not yet available: '/        &
       ' -----------------------------------------------------')
  !
945 FORMAT (/'  Successfully requested output fields : '/            &
       ' -----------------------------------------------------')
946 FORMAT ( '      ',a,1x,a)
  !
948 FORMAT (/'  Additional GRIB parameters : '/                     &
       ' -----------------------------------------------------'/ &
       '      Run time           : ',a/                       &
       '      GRIB center ID     : ',i4/                      &
       '      GRIB gen. proc. ID : ',i4/                      &
       '      GRIB grid ID       : ',i4/                      &
       '      GRIB GDS parameter : ',i4)
  !
970 FORMAT (//'  Generating file '/                                 &
       ' -----------------------------------------------------')
971 FORMAT ( '      Data for ',a,'  ',i3,'H forecast.')
972 FORMAT ( '      Data for ',a,'       hindcast.')
  !
999 FORMAT (/'  End of program '/                                   &
       ' ========================================='/          &
       '         WAVEWATCH III GRIB output '/)
  !
#ifdef W3_T
9050 FORMAT ( ' TEST W3GRIB : KPDS : ',13i4/                      &
       '                      ',12i4)
9051 FORMAT ( ' TEST W3GRIB : KGDS : ',8i6/                       &
       '                      ',8i6/                       &
       '                      ',6i6)
#endif
  !
1000 FORMAT (/' *** WAVEWATCH III ERROR IN W3GRIB : '/               &
       '     ERROR IN OPENING INPUT FILE'/                    &
       '     IOSTAT =',i5/)
  !
1001 FORMAT (/' *** WAVEWATCH III ERROR IN W3GRIB : '/               &
       '     PREMATURE END OF INPUT FILE'/)
  !
1002 FORMAT (/' *** WAVEWATCH III ERROR IN W3GRIB : '/               &
       '     ERROR IN READING FROM INPUT FILE'/               &
       '     IOSTAT =',i5/)
  !
1010 FORMAT (/' *** WAVEWATCH-III ERROR IN W3GRIB : '/          &
       '     GRIB REQUIRES SPHERICAL GRID'/)
#ifdef W3_NCEP2
1011 FORMAT (/' CHOSEN GRID TYPE: : ',a/)
#endif
  !/
  !/ Internal subroutine W3EXGB ---------------------------------------- /
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3exgb ( NX, NY, NSEA )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |             A. Chawla             |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         22-Mar-2021 |
    !/                  +-----------------------------------+
    !/
    !/    10-Jun-1999 : Final FORTRAN 77                    ( version 1.18 )
    !/    24-Jan-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/                  Massive changes to logistics.
    !/    29-Apr-2002 : Adding output fields 17-18.         ( version 2.20 )
    !/    24-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    18-May-2007 : Update GRIB1 for partitioning.      ( version 3.11 )
    !/    16-Jul-2007 : Adding GRIB2 capability             ( version 3.11 )
    !/                  (A. Chawla)
    !/    22-Mar-2021 : New coupling fields output          ( version 7.13 )
    !/
    !  1. Purpose :
    !
    !     Perform actual GRIB output.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       NX, NY, NSEA
    !               Int.  I  Array dimensions.
    !     ----------------------------------------------------------------
    !
    !     Internal parameters
    !     ----------------------------------------------------------------
    !       X1, X2, XX, XY
    !               R.A.  Output fields
    !       BITMAP  L.A.  Data / no data bitmap
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE     Subr. W3SERVMD Subroutine tracing.
    !      EXTCDE     Subr.   Id.    Abort program as graceful as possible.
    !      W3S2XY     Subr.   Id.    Convert from storage to spatial grid.
    !      PUTGB      Subr.          NCEP GRIB1 library routine.
    !      GRIBCREATE Subr.          NCEP GRIB2 library routine.
    !      ADDGRID    Subr.          NCEP GRIB2 library routine.
    !      ADDFIELD   Subr.          NCEP GRIB2 library routine.
    !      GRIBEND    Subr.          NCEP GRIB2 library routine.
    !      WRYTE      Subr.          NCEP GRIB2 library routine.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !     Program in which it is contained.
    !
    !  6. Error messages :
    !
    !     None.
    !
    !  7. Remarks :
    !
    !     - Note that arrays CX and CY of the main program now contain
    !       the absolute current speed and direction respectively.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S      Enable subroutine tracing.
    !     !/T      Enable test output.
    !     !/NCEP2  NCEP IBM calls to GRIB2 packer (follows updated grib2
    !     tables under verification as of 02/10/2012).
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3servmd, ONLY : w3s2xy
#ifdef W3_RTD
    USE w3servmd, ONLY : w3thrtn, w3xyrtn
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: NX, NY, NSEA
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: J, IXY, NDATA
    INTEGER                 :: IO
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    ::  X1(NX*NY), X2(NX*NY), XX(NX*NY),    &
         XY(NX*NY), CABS, UABS,              &
         YY(NX*NY,0:NOSWLL), KPDS5A, KPDS5B, &
         KPDS5A1(3)
    LOGICAL*1               ::  BITMAP(NX*NY)
    LOGICAL                 :: FLONE, FLTWO, FLDIR, FLTRI, FLPRT
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3EXGB')
#endif
    !
#ifdef W3_T
    WRITE (ndst,9000) ((flreq(ifi,ifj),ifj=1,ngrpp), ifi=1,nogrp)
    WRITE (ndst,9001) ndsdat, kpds, kgds
#endif
    !
    !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ! 1.  Preparations
    !
    x1     = undef
    x2     = undef
    xx     = undef
    xy     = undef
    yy     = undef
    !
    !
    !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ! 2.  Loop over output fields.
    !
    DO ifi=1, nogrp
      DO ifj=1, ngrpp
        IF ( flreq(ifi,ifj) ) THEN
 
 
          !
          ! Initialize array dimension flags
          !
          flone = .false.
          fltwo = .false.
          fldir = .false.
          fltri = .false.
          flprt = .false.
          !
#ifdef W3_T
          WRITE (ndst,9020) idout(ifi,ifj)
#endif
          !
          ! 2.a Set output arrays and parameters
          !
          !     Water depth
          !
          IF ( ifi .EQ. 1 .AND. ifj .EQ.  1 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 14
            kpds(1) = 4
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, dw(1:nsea)          &
                 , mapsf, x1 )
            !
            !     Current
            !
          ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ.  2 ) THEN
            fltwo = .true.
#ifdef W3_NCEP2
            kpds(2) = 1
            kpds(1) = 1
#endif
#ifdef W3_RTD
            ! Rotate x,y vector back to standard pole
            IF ( flagunr ) CALL w3xyrtn(nsea, cx, cy, angld)
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, cx(1:nsea)          &
                 , mapsf, xx )
            CALL w3s2xy ( nsea, nsea, nx, ny, cy(1:nsea)          &
                 , mapsf, xy )
            DO isea=1, nsea
              IF (cx(isea) .NE. undef) THEN
                cabs   = sqrt(cx(isea)**2+cy(isea)**2)
                IF ( cabs .GT. 0.001 ) THEN
                  cy(isea) = mod( 630. -                         &
                       rade*atan2(cy(isea),cx(isea)) , 360. )
                ELSE
                  cy(isea) = 0.
                END IF
              ELSE
                cabs = undef
                cy(isea) = undef
              END IF
              cx(isea) = cabs
            END DO
            CALL w3s2xy ( nsea, nsea, nx, ny, cx(1:nsea)          &
                 , mapsf, x1 )
            CALL w3s2xy ( nsea, nsea, nx, ny, cy(1:nsea)          &
                 , mapsf, x2 )
            !
            !     Wind speed
            !
          ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ.  3 ) THEN
            fltwo = .true.
#ifdef W3_NCEP2
            kpds(2) = 1
            kpds(1) = 2
            listsec0(1) = 0
#endif
#ifdef W3_RTD
            ! Rotate x,y vector back to standard pole
            IF ( flagunr ) CALL w3xyrtn(nsea, ua, ud, angld)
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, ua(1:nsea)          &
                 , mapsf, xx )
            CALL w3s2xy ( nsea, nsea, nx, ny, ud(1:nsea)          &
                 , mapsf, xy )
            DO isea=1, nsea
              IF (ua(isea) .NE. undef) THEN
                uabs   = sqrt(ua(isea)**2+ud(isea)**2)
                IF ( uabs .GT. 0.001 ) THEN
                  ud(isea) = mod( 630. -                         &
                       rade*atan2(ud(isea),ua(isea)) , 360. )
                ELSE
                  ud(isea) = 0.
                END IF
              ELSE
                uabs = undef
                ud(isea) = undef
              END IF
              ua(isea) = uabs
            END DO
            CALL w3s2xy ( nsea, nsea, nx, ny, ua(1:nsea)          &
                 , mapsf, x1 )
            CALL w3s2xy ( nsea, nsea, nx, ny, ud(1:nsea)          &
                 , mapsf, x2 )
            !
            !     Air-sea temp. dif.
            !
          ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ.  4 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 255
            kpds(1) = 3
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, as(1:nsea)          &
                 , mapsf, x1 )
            !
            !     Water level
            !
          ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 5 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 1
            kpds(1) = 3
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, wlv   , mapsf, x1 )
            !
            !     Ice concentration
            !
          ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 6 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 0
            kpds(1) = 2
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, ice   , mapsf, x1 )
            !
            !     Atmospheric momentum
            !
          ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ.  8 ) THEN
            fltwo = .true.
#ifdef W3_NCEP2
            kpds(2) = 1
            kpds(1) = 2
            listsec0(1) = 0
#endif
#ifdef W3_RTD
            ! Rotate x,y vector back to standard pole
            IF ( flagunr ) CALL w3xyrtn(nsea, taua, tauadir, angld)
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, taua(1:nsea)        &
                 , mapsf, xx )
            CALL w3s2xy ( nsea, nsea, nx, ny, tauadir(1:nsea)     &
                 , mapsf, xy )
            DO isea=1, nsea
              IF (taua(isea) .NE. undef) THEN
                uabs   = sqrt(taua(isea)**2+tauadir(isea)**2)
                IF ( uabs .GT. 0.001 ) THEN
                  tauadir(isea) = mod( 630. -                 &
                       rade*atan2(tauadir(isea),taua(isea)) , 360. )
                ELSE
                  tauadir(isea) = 0.
                END IF
              ELSE
                uabs = undef
                tauadir(isea) = undef
              END IF
              taua(isea) = uabs
            END DO
            CALL w3s2xy ( nsea, nsea, nx, ny, taua(1:nsea)        &
                 , mapsf, x1 )
            CALL w3s2xy ( nsea, nsea, nx, ny, tauadir(1:nsea)     &
                 , mapsf, x2 )
            !
            !     Air density
            !
          ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 9 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 0
            kpds(1) = 2
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, rhoair, mapsf, x1 )
            !
            !     Significant wave height
            !
          ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ.  1 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 3
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, hs    , mapsf, x1 )
            !
            !     Mean wave length
            !
          ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ.  2 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 193
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, wlm   , mapsf, x1 )
            !
            !     Mean wave period (based on second moment)
            !
          ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ.  3 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            if ((gen_pro.eq.1) .or. (gen_pro.eq.0)) then
              kpds(2) = 28
            else
              kpds(2) = 25
            endif
#endif
 
            CALL w3s2xy ( nsea, nsea, nx, ny, t02   , mapsf, x1 )
            !
            !     Mean wave period (based on first moment)
            !
          ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ.  4 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 15
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, t0m1   , mapsf, x1 )
            !
            !     Mean wave period (based on first inverse moment)
            !
          ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ.  5 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            if ((gen_pro.eq.1) .or. (gen_pro.eq.0)) then
              kpds(2) = 34
            else
              kpds(2) = 15
            endif
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, t01   , mapsf, x1 )
            !
            !     Peak frequency
            !
          ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 6 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 11
#endif
            DO isea=1, nsea
              IF ( fp0(isea) .NE. undef .AND. fp0(isea) .NE. 0 ) THEN
                fp0(isea) = 1. / max(fr1,fp0(isea)) ! Limit FP to lowest discrete frequency
              END IF
            END DO
            CALL w3s2xy ( nsea, nsea, nx, ny, fp0   , mapsf, x1 )
            !
            !
            !     Mean wave direction
            !
          ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ.  7 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 14
#endif
#ifdef W3_RTD
            ! Rotate direction back to standard pole
            IF ( flagunr ) CALL w3thrtn(nsea, thm, angld, .false.)
#endif
            DO isea=1, nsea
              IF ( thm(isea) .NE. undef )                         &
                   thm(isea) = mod( 630. - rade*thm(isea) , 360. )
            END DO
            CALL w3s2xy ( nsea, nsea, nx, ny, thm   , mapsf, x1 )
            !
            !     Directional spread
            !
          ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 8 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 31
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, ths   , mapsf, x1 )
            !
            !     Peak direction
            !
          ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 9 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            if ((gen_pro.eq.1) .or. (gen_pro.eq.0)) then
              kpds(2) = 46
            else
              kpds(2) = 10
            endif
#endif
#ifdef W3_RTD
            ! Rotate direction back to standard pole
            IF ( flagunr ) CALL w3thrtn(nsea, thp0, angld, .false.)
#endif
            DO isea=1, nsea
              IF ( thp0(isea) .NE. undef ) THEN
                thp0(isea) = mod( 630-rade*thp0(isea) , 360. )
              END IF
            END DO
            CALL w3s2xy ( nsea, nsea, nx, ny, thp0  , mapsf, x1 )
            !
            !     Mean wave number
            !
          ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 19 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 255
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, wnmean, mapsf, x1 )
            !
            !     Partitioned wave height
            !
          ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 1 ) THEN
            flprt = .true.
#ifdef W3_NCEP2
            kpds5a  = 5
            kpds5b  = 8
            if ((gen_pro.eq.1) .or. (gen_pro.eq.0)) then
              kpds5a1(1)  =47
              kpds5a1(2)  =48
              kpds5a1(3)  =49
            else
              kpds5b  = 8
            endif
#endif
            CALL w3s2xy                                           &
                 ( nsea, nsea, nx, ny, phs(:,0), mapsf, yy(:,0) )
            DO i=1, noswll
              CALL w3s2xy                                         &
                   ( nsea, nsea, nx, ny, phs(:,i), mapsf, yy(:,i) )
            END DO
            !
            !     Partitioned peak period
            !
          ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 2 ) THEN
            flprt = .true.
#ifdef W3_NCEP2
            kpds5a  = 6
            kpds5b = 9
            if ((gen_pro.eq.1) .or. (gen_pro.eq.0)) then
              kpds5a1(1)  = 50
              kpds5a1(2)  = 51
              kpds5a1(3)  = 52
            else
              kpds5b = 9
            endif
#endif
            CALL w3s2xy                                           &
                 ( nsea, nsea, nx, ny, ptp(:,0), mapsf, yy(:,0) )
            DO i=1, noswll
              CALL w3s2xy                                         &
                   ( nsea, nsea, nx, ny, ptp(:,i), mapsf, yy(:,i) )
            END DO
            !
            !     Partitioned peak wave length
            !
          ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 3 ) THEN
            flprt = .true.
#ifdef W3_NCEP2
            kpds5a  =  193
            kpds5b  =  193
#endif
            CALL w3s2xy                                           &
                 ( nsea, nsea, nx, ny, plp(:,0), mapsf, yy(:,0) )
            DO i=1, noswll
              CALL w3s2xy                                         &
                   ( nsea, nsea, nx, ny, plp(:,i), mapsf, yy(:,i) )
            END DO
            !
            !     Partitioned mean direction
            !
          ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 4 ) THEN
            flprt = .true.
#ifdef W3_NCEP2
            kpds5a  = 4
            kpds5b = 7
            if ((gen_pro.eq.1) .or. (gen_pro.eq.0)) then
              kpds5a1(1)  = 53
              kpds5a1(2)  = 54
              kpds5a1(3)  = 55
            else
              kpds5b = 7
            endif
#endif
#ifdef W3_RTD
            DO i = 0,noswll
              ! Rotate direction back to standard pole
              IF ( flagunr ) CALL w3thrtn(nsea, pdir(:,i), angld, .false.)
            END DO
#endif
            DO isea = 1,nsea
              DO i = 0,noswll
                IF ( pdir(isea,i) .NE. undef ) THEN
                  pdir(isea,i) = mod( 630 - rade*pdir(isea,i) , 360. )
                END IF
              END DO
            END DO
            CALL w3s2xy                                           &
                 ( nsea, nsea, nx, ny, pdir(:,0), mapsf, yy(:,0) )
            DO i=1, noswll
              CALL w3s2xy                                         &
                   ( nsea, nsea, nx, ny, pdir(:,i), mapsf, yy(:,i) )
            END DO
            !
            !     Partitioned Directional spread
            !
          ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 5 ) THEN
            flprt = .true.
#ifdef W3_NCEP2
            kpds5a  =  32
            kpds5b  =  33
#endif
            CALL w3s2xy                                           &
                 ( nsea, nsea, nx, ny, psi(:,0), mapsf, yy(:,0) )
            DO i=1, noswll
              CALL w3s2xy                                         &
                   ( nsea, nsea, nx, ny, psi(:,i), mapsf, yy(:,i) )
            END DO
            !
            !     Partitioned wind sea fraction
            !
          ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 6 ) THEN
            flprt = .true.
#ifdef W3_NCEP2
            kpds5a  =  255
            kpds5b  =  255
#endif
            CALL w3s2xy                                           &
                 ( nsea, nsea, nx, ny, pws(:,0), mapsf, yy(:,0) )
            DO i=1, noswll
              CALL w3s2xy                                         &
                   ( nsea, nsea, nx, ny, pws(:,i), mapsf, yy(:,i) )
            END DO
            !
            !     Total wind sea fraction
            !
          ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 16 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 255
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, pwst  , mapsf, x1 )
            !
            !     Number of fields in partition
            !
          ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 17 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 255
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, pnr   , mapsf, x1 )
            !
            !     Friction velocity
            !
          ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ.  1 ) THEN
            fltwo = .true.
#ifdef W3_NCEP2
            kpds(2) = 17
            kpds(1) = 1
#endif
#ifdef W3_RTD
            ! Rotate x,y vector back to standard pole
            IF ( flagunr ) CALL w3xyrtn(nsea, ust, ustdir, angld)
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, ust(1:nsea)         &
                 , mapsf, x1 )
            CALL w3s2xy ( nsea, nsea, nx, ny, ustdir(1:nsea)      &
                 , mapsf, x2 )
            !
            !     Average source term time step
            !
          ELSE IF ( ifi .EQ. 9 .AND. ifj .EQ. 1 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 255
#endif
            DO isea=1, nsea
              IF ( dtdyn(isea) .NE. undef )                       &
                   dtdyn(isea) = dtdyn(isea) / 60.
            END DO
            CALL w3s2xy ( nsea, nsea, nx, ny, dtdyn , mapsf, x1 )
            !
            !     Cut-off frequency
            !
          ELSE IF ( ifi .EQ. 9 .AND. ifj .EQ. 2 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 255
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, fcut  , mapsf, x1 )
            !
            !     CFL Maximum (in spatial space)
            !
          ELSE IF ( ifi .EQ. 9 .AND. ifj .EQ. 3 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 255
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, cflxymax  , mapsf, x1 )
            !
            !     CFL Maximum (in spectral space)
            !
          ELSE IF ( ifi .EQ. 9 .AND. ifj .EQ. 4 ) THEN
            flone = .true.
#ifdef W3_NCEP2
            kpds(2) = 255
#endif
            CALL w3s2xy ( nsea, nsea, nx, ny, cflthmax  , mapsf, x1 )
            !
          ELSE
            WRITE (ndse,999)
            CALL extcde ( 1 )
            !
          END IF
          !
          ! 3   Perform output
          !
          ndata  = nx*ny
          !
          ! 3.a Partitioned data
          !
          IF ( flprt ) THEN
            !
#ifdef W3_NCEP2
            kpds(2) = kpds5a
#endif
            DO ixy=1, nx*ny
              bitmap(ixy) = yy(ixy,0) .NE. undef
            END DO
#ifdef W3_NCEP2
            CALL gribcreate (cgrib,lcgrib,listsec0,listsec1,io)
            IF (io .NE. 0) GOTO 810
            CALL addgrid (cgrib,lcgrib,igds,kgds,200,ideflist,   &
                 idefnum, io)
            IF (io .NE. 0) GOTO 820
            CALL addfield (cgrib,lcgrib,kpdsnum,kpds,200,        &
                 coordlist, numcoord, idrsnum, idrs,   &
                 200,yy(:,0), ndata, ibmp, bitmap, io)
            IF (io .NE. 0) GOTO 820
            CALL gribend (cgrib, lcgrib, lengrib, io)
            IF (io .NE. 0) GOTO 830
            CALL wryte (ndsdat, lengrib, cgrib)
#endif
            !
#ifdef W3_NCEP2
            if ((gen_pro.eq.0) .or. (gen_pro.eq.1)) then
              kpds(10) = 241
#endif
              DO i=1, noswll
#ifdef W3_NCEP2
                kpds(2) = kpds5a1(i)
                kpds(12) = i
#endif
                DO ixy=1, nx*ny
                  bitmap(ixy) = yy(ixy,i) .NE. undef
                END DO
#ifdef W3_NCEP2
                CALL gribcreate (cgrib,lcgrib,listsec0,listsec1,io)
                IF (io .NE. 0) GOTO 810
                CALL addgrid (cgrib,lcgrib,igds,kgds,200,ideflist,   &
                     idefnum, io)
                IF (io .NE. 0) GOTO 820
                CALL addfield (cgrib,lcgrib,kpdsnum,kpds,200,        &
                     coordlist, numcoord, idrsnum, idrs,   &
                     200,yy(:,i), ndata, ibmp, bitmap, io)
                IF (io .NE. 0) GOTO 820
                CALL gribend (cgrib, lcgrib, lengrib, io)
                IF (io .NE. 0) GOTO 830
                CALL wryte (ndsdat, lengrib, cgrib)
#endif
              END DO
#ifdef W3_NCEP2
            ELSE
              kpds(2) = kpds5b
              kpds(10) = 241
#endif
              DO i=1, noswll
#ifdef W3_NCEP2
                kpds(12) = i
#endif
                DO ixy=1, nx*ny
                  bitmap(ixy) = yy(ixy,i) .NE. undef
                END DO
#ifdef W3_NCEP2
                CALL gribcreate (cgrib,lcgrib,listsec0,listsec1,io)
                IF (io .NE. 0) GOTO 810
                CALL addgrid (cgrib,lcgrib,igds,kgds,200,ideflist,   &
                     idefnum, io)
                IF (io .NE. 0) GOTO 820
                CALL addfield (cgrib,lcgrib,kpdsnum,kpds,200,        &
                     coordlist, numcoord, idrsnum, idrs,   &
                     200,yy(:,i), ndata, ibmp, bitmap, io)
                IF (io .NE. 0) GOTO 820
                CALL gribend (cgrib, lcgrib, lengrib, io)
                IF (io .NE. 0) GOTO 830
                CALL wryte (ndsdat, lengrib, cgrib)
#endif
              END DO
#ifdef W3_NCEP2
            ENDIF
            kpds(10) = 1
            kpds(12) = 1
#endif
            !
            ! 3.b Other data
            !
          ELSE IF (flone) THEN
            !
            DO ixy=1, nx*ny
              bitmap(ixy) = x1(ixy) .NE. undef
            END DO
            !
#ifdef W3_NCEP2
            CALL gribcreate (cgrib,lcgrib,listsec0,listsec1,io)
            IF (io .NE. 0) GOTO 810
            CALL addgrid (cgrib,lcgrib,igds,kgds,200,ideflist,   &
                 idefnum, io)
            IF (io .NE. 0) GOTO 820
            CALL addfield (cgrib,lcgrib,kpdsnum,kpds,200,        &
                 coordlist, numcoord, idrsnum, idrs,   &
                 200,x1, ndata, ibmp, bitmap, io)
            IF (io .NE. 0) GOTO 820
            CALL gribend (cgrib, lcgrib, lengrib, io)
            IF (io .NE. 0) GOTO 830
            CALL wryte (ndsdat, lengrib, cgrib)
#endif
            !
          ELSE IF ( fltwo ) THEN
            !
            DO ixy=1, nx*ny
              bitmap(ixy) = x1(ixy) .NE. undef
            END DO
#ifdef W3_NCEP2
            CALL gribcreate (cgrib,lcgrib,listsec0,listsec1,io)
            IF (io .NE. 0) GOTO 810
            CALL addgrid (cgrib,lcgrib,igds,kgds,200,ideflist, &
                 idefnum, io)
            IF (io .NE. 0) GOTO 820
            CALL addfield (cgrib,lcgrib,kpdsnum,kpds,200, &
                 coordlist, numcoord, idrsnum, idrs, &
                 200,x1, ndata, ibmp, bitmap, io)
            IF (io .NE. 0) GOTO 820
            CALL gribend (cgrib, lcgrib, lengrib, io)
            IF (io .NE. 0) GOTO 830
            CALL wryte (ndsdat, lengrib, cgrib)
#endif
 
#ifdef W3_NCEP2
            kpds(2) = 0
            CALL gribcreate (cgrib,lcgrib,listsec0,listsec1,io)
            IF (io .NE. 0) GOTO 810
            CALL addgrid (cgrib,lcgrib,igds,kgds,200,ideflist, &
                 idefnum, io)
            IF (io .NE. 0) GOTO 820
            CALL addfield (cgrib,lcgrib,kpdsnum,kpds,200,      &
                 coordlist, numcoord, idrsnum, idrs, &
                 200,x2, ndata, ibmp, bitmap, io)
            IF (io .NE. 0) GOTO 820
            CALL gribend (cgrib, lcgrib, lengrib, io)
            IF (io .NE. 0) GOTO 830
            CALL wryte (ndsdat, lengrib, cgrib)
            kpds(2) = 2
            CALL gribcreate (cgrib,lcgrib,listsec0,listsec1,io)
            IF (io .NE. 0) GOTO 810
            CALL addgrid (cgrib,lcgrib,igds,kgds,200,ideflist, &
                 idefnum, io)
            IF (io .NE. 0) GOTO 820
            CALL addfield (cgrib,lcgrib,kpdsnum,kpds,200,      &
                 coordlist, numcoord, idrsnum, idrs, &
                 200,xx, ndata, ibmp, bitmap, io)
            IF (io .NE. 0) GOTO 820
            CALL gribend (cgrib, lcgrib, lengrib, io)
            IF (io .NE. 0) GOTO 830
            CALL wryte (ndsdat, lengrib, cgrib)
            kpds(2) = 3
            CALL gribcreate (cgrib,lcgrib,listsec0,listsec1,io)
            IF (io .NE. 0) GOTO 810
            CALL addgrid (cgrib,lcgrib,igds,kgds,200,ideflist, &
                 idefnum, io)
            IF (io .NE. 0) GOTO 820
            CALL addfield (cgrib,lcgrib,kpdsnum,kpds,200,      &
                 coordlist, numcoord, idrsnum, idrs, &
                 200,xy, ndata, ibmp, bitmap, io)
            IF (io .NE. 0) GOTO 820
            CALL gribend (cgrib, lcgrib, lengrib, io)
            IF (io .NE. 0) GOTO 830
            CALL wryte (ndsdat, lengrib, cgrib)
#endif
            !
          END IF
#ifdef W3_NCEP2
          listsec0(1) = 10
          kpds(1)     = 0
#endif
          !
          ! ... End of fields loop
          !
        END IF
      END DO
    END DO
    !
    RETURN
    !
    ! Error escape locations
    !
#ifdef W3_NCEP2
810 CONTINUE
    WRITE (ndse,1010) io
    CALL extcde ( 20 )
820 CONTINUE
    WRITE (ndse,1020) io
    CALL extcde ( 30 )
830 CONTINUE
    WRITE (ndse,1030) io
    CALL extcde ( 40 )
#endif
    !
    ! Formats
    !
999 FORMAT (/' *** WAVEWATCH III ERROR IN W3EXGB :'/                &
         '     PLEASE UPDATE FIELDS !!! '/)
    !
#ifdef W3_NCEP2
1000 FORMAT (/' *** WAVEWATCH III ERROR IN W3EXGB : '/               &
         '     ERROR IN OPENING OUTPUT FILE'/                   &
         '     IOSTAT =',i5/)
#endif
    !
#ifdef W3_NCEP2
1010 FORMAT (/' *** WAVEWATCH III ERROR IN W3EXGB : '/               &
         '     ERROR CREATING NEW GRIB2 FIELD'/                 &
         '     IOSTAT =',i5/)
#endif
    !
#ifdef W3_NCEP2
1020 FORMAT (/' *** WAVEWATCH III ERROR IN W3EXGB : '/               &
         '     ERROR ADDING GRIB2 FIELD'/                       &
         '     IOSTAT =',i5/)
#endif
    !
#ifdef W3_NCEP2
1030 FORMAT (/' *** WAVEWATCH III ERROR IN W3EXGB : '/               &
         '     ERROR ENDING GRIB2 MESSAGE'/                     &
         '     IOSTAT =',i5/)
#endif
    !
#ifdef W3_T
9000 FORMAT (' TEST W3EXGB : FLAGS  :',40l2)
9001 FORMAT (' TEST W3EXGB : NDSDAT :',i4/                        &
         '               KPDS   :',13i4/                      &
         '                       ',12i4/                      &
         '               KGDS   :',8i6/                       &
         '                       ',8i6/                       &
         '                       ',6i6)
#endif
    !
#ifdef W3_T
9012 FORMAT (' TEST W3EXGB : BLOK PARS    : ',3i4)
9014 FORMAT ('           BASE NAME : ',a)
#endif
    !
#ifdef W3_T
9020 FORMAT (' TEST W3EXGB : OUTPUT FIELD : ',a)
#endif
    !/
    !/ End of W3EXGB ----------------------------------------------------- /
    !/
  END SUBROUTINE w3exgb
  !/
  !/ End of W3GRIB ----------------------------------------------------- /
  !/
END PROGRAM w3grib