w3iopomd.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
#define nf90_err(ncerr) nf90_err_check(ncerr, __LINE__)
!/ ------------------------------------------------------------------- /
MODULE w3iopomd
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         05-Jun-2018 |
  !/                  +-----------------------------------+
  !/
  !/    25-Jan-2001 : Origination.                        ( version 2.00 )
  !/    24-Jan-2001 : Flat grid version.                  ( version 2.06 )
  !/    11-Jun-2001 : Clean-up.                           ( version 2.11 )
  !/    10-Nov-2004 : Multiple grid version.              ( version 3.06 )
  !/    27-Jun-2006 : Adding file name preamble.          ( version 3.09 )
  !/    25-Jul-2006 : Adding grid ID per point.           ( version 3.10 )
  !/    01-May-2007 : Move O7a output from W3INIT.        ( 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)
  !/    30-Oct-2009 : Implement curvilinear grid type.    ( version 3.14 )
  !/                  (W. E. Rogers & T. J. Campbell, NRL)
  !/    29-Oct-2010 : Implement unstructured grid         ( version 3.14.4 )
  !/                  (A. Roland and F. Ardhuin)
  !/    06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
  !/                  specify index closure for a grid.   ( version 3.14 )
  !/                  (T. J. Campbell, NRL)
  !/    12-Jun-2012 : Add /RTD option or rotated grid option.
  !/                  (Jian-Guo Li)                       ( version 4.06 )
  !/    02-Sep-2012 : Clean up of open BC for UG grids    ( version 4.07 )
  !/    25-Feb-2013 : ITOUT=0 bug correction for UG grids ( version 4.08 )
  !/    11-Nov-2013 : SMC and rotated grid incorporated in the main
  !/                  trunk                               ( version 4.13 )
  !/    05-Jun-2018 : Add SETUP                           ( version 6.04 )
  !/    04-Oct-2019 : Optional one file per output stride ( version 7.00 )
  !/                  (R. Padilla-Hernandez & J.H. Alves)
  !/
  !/    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 :
  !
  !     Process point output.
  !
  !  2. Variables and types :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      VEROPT    C*10  Private  Point output file version number.
  !      IDSTR     C*31  Private  Point output file ID string.
  !     ----------------------------------------------------------------
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3IOPP    Subr. Public   Preprocessing of point output req.
  !      W3IOPE    Subr. Public   Extract point data from grid.
  !      W3IOPO    Subr. Public   Point data IO.
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      W3SETO    Subr. W3ODATMD Data structure management.
  !      W3SETG    Subr. W3GDATMD Data structure management.
  !      W3SETW    Subr. W3WDATMD Data structure management.
  !      W3DMO2    Subr. W3ODATMD Data structure management.
  !      STRACE    Subr. W3SERVMD Subroutine tracing.
  !      EXTCDE    Subr. W3SERVMD Program abort with exit code.
  !      MPI_STARTALL, MPIWAITALL
  !                Subr.          MPI persistent communication routines.
  !     ----------------------------------------------------------------
  !
  !  5. Remarks :
  !
  !     - Allocation of allocatable arrays takes place at different
  !       places throughout the code, in W3IOPP on write, and in
  !       W3IOPO on read.
  !
  !  6. Switches :
  !
  !       !/S     Enable subroutine tracing.
  !       !/T     Enable test output.
  !
  !       !/SHRD  Switch for shared / distributed memory architecture.
  !       !/DIST  Id.
  !       !/MPI   MPI message passing.
  !
  !       !/O7a   Diagnostic output for output points.
  !
  !  7. Source code :
  !
  !/ ------------------------------------------------------------------- /
  PUBLIC
  !/
  !/ Private parameter statements (ID strings)
  !/
  CHARACTER(LEN=10), PARAMETER, PRIVATE :: VEROPT = '2021-04-06'
  CHARACTER(LEN=31), PARAMETER, PRIVATE ::                        &
       IDSTR = 'WAVEWATCH III POINT OUTPUT FILE'
 
  character(*), parameter, private :: DNAME_NOPTS = 'NOPTS'
 
  character(*), parameter, private :: DNAME_NSPEC = 'NSPEC'
 
  character(*), parameter, private :: DNAME_VSIZE = 'VSIZE'
 
  character(*), parameter, private :: DNAME_NAMELEN = 'NAMELEN'
 
  character(*), parameter, private :: DNAME_GRDIDLEN = 'GRDIDLEN'
 
  character(*), parameter, private :: DNAME_TIME = 'TIME'
 
  character(*), parameter, private :: DNAME_WW3TIME = 'WW3TIME'
 
  character(*), parameter, private :: VNAME_NK = 'NK'
 
  character(*), parameter, private :: VNAME_NTH = 'NTH'
 
  character(*), parameter, private :: VNAME_PTLOC = 'PTLOC'
 
  character(*), parameter, private :: VNAME_PTNME = 'PTNME'
 
  character(*), parameter, private :: VNAME_TIME = 'TIME'
 
  character(*), parameter, private :: VNAME_WW3TIME = 'WW3TIME'
 
  character(*), parameter, private :: VNAME_DPO = 'DPO'
 
  character(*), parameter, private :: VNAME_WAO = 'WAO'
 
  character(*), parameter, private :: VNAME_WDO = 'WDO'
 
  character(*), parameter, private :: VNAME_TAUAO = 'TAUAO'
 
  character(*), parameter, private :: VNAME_TAUDO = 'TAUDO'
 
  character(*), parameter, private :: VNAME_DAIRO = 'DAIRO'
 
  character(*), parameter, private :: VNAME_ZET_SETO = 'ZET_SETO'
 
  character(*), parameter, private :: VNAME_ASO = 'ASO'
 
  character(*), parameter, private :: VNAME_CAO = 'CAO'
 
  character(*), parameter, private :: VNAME_CDO = 'CDO'
 
  character(*), parameter, private :: VNAME_ICEO = 'ICEO'
 
  character(*), parameter, private :: VNAME_ICEHO = 'ICEHO'
 
  character(*), parameter, private :: VNAME_ICEFO = 'ICEFO'
 
  character(*), parameter, private :: VNAME_GRDID = 'GRDID'
 
  character(*), parameter, private :: VNAME_SPCO = 'SPCO'
 
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3iopp ( NPT, XPT, YPT, PNAMES, IMOD )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         02-Sep-2012 |
    !/                  +-----------------------------------+
    !/
    !/    14-Jan-1999 : Distributed FORTRAN 77 version.     ( version 1.18 )
    !/    30-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/                  Major changes to logistics.
    !/    24-Jan-2001 : Flat grid version.                  ( version 2.06 )
    !/    09-Nov-2004 : Multiple grid version.              ( version 3.06 )
    !/    25-Jul-2006 : Adding grid ID per point.           ( version 3.10 )
    !/    01-May-2007 : Move O7a output from W3INIT.        ( version 3.11 )
    !/    30-Oct-2009 : Implement run-time grid selection.  ( version 3.14 )
    !/                  (W. E. Rogers & T. J. Campbell, NRL)
    !/    30-Oct-2009 : Implement curvilinear grid type.    ( version 3.14 )
    !/                  (W. E. Rogers & T. J. Campbell, NRL)
    !/    06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
    !/                  specify index closure for a grid.   ( version 3.14 )
    !/                  (T. J. Campbell, NRL)
    !/    12-Jun-2012 : Add /RTD option or rotated grid option.
    !/                  (Jian-Guo Li)                       ( version 4.06 )
    !/    02-Sep-2012 : Clean up of open BC for UG grids    ( version 4.07 )
    !/    01-Mar-2018 : Add option to unrotate spectra      ( version 6.02 )
    !/                  from RTD grid models
    !/
    !  1. Purpose :
    !
    !     Preprocessing of point output.
    !
    !  2. Method :
    !
    !     Check location of points in grid and calculate interpolation
    !     factors.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       NPT     Int.   I   Number of output points in input.
    !       XPT     R.A.  I/O  X (longitude) coordinates of output points.
    !       YPT     R.A.  I/O  Id. Y.
    !       PNAMES  C*40   I   Names of output points.
    !       IMOD    Int.   I   Grid ID number.
    !     ----------------------------------------------------------------
    !
    !     Local data
    !     ----------------------------------------------------------------
    !       ACC     Real  "Accuracy" factor to determine if output point
    !                     is grid point.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3INIT    Subr. W3INITMD Wave model initialization routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !     - Warnings for points out of the grid or on land.
    !
    !  7. Remarks :
    !
    !     - The output points are obtained by bi-linear interpolation from
    !       the spectra at the grid points. Given the possibility of ice
    !       coverage, the actual interpolation factors can only be
    !       determined at the actual output time. Hence only the basic
    !       bilinear interpolation factors are stored.
    !     - Implementation of the /O7a diagnostic output section is
    !       currently incomplete and non-functional for curvilinear grids
    !       and/or tripole grids
    !
    !  8. Structure :
    !
    !     -------------------------------------------
    !      Determine grid range
    !      do for all defined points
    !      -----------------------------------------
    !        Check if point within grid
    !        Calculate interpolation data
    !        Check if point not on land
    !        Store interpolation data
    !     -------------------------------------------
    !
    !  9. Switches :
    !
    !       !/S     Enable subroutine tracing.
    !       !/T     Test output.
    !
    !       !/O7a   Diagnostic output for output points.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gsrumd
    USE w3gdatmd, ONLY: nth, nk, nspec, nx, ny, x0, y0, sx, gsu,&
         rlgtype, clgtype, ungtype, gtype, flagll,   &
         iclose,iclose_none,iclose_smpl,iclose_trpl, &
         mapsta, mapfs, filext, zb, trnx, trny
    USE w3gdatmd, ONLY: trigp,maxx, maxy, dxymax
#ifdef W3_RTD
    !!  Use rotated N-Pole lat/lon and conversion sub.  JGLi12Jun2012
    USE w3gdatmd, ONLY: polat, polon, flagunr
    USE w3servmd, ONLY: w3lltoeq
#endif
    USE w3odatmd, ONLY: w3dmo2
    USE w3odatmd, ONLY: ndse, ndst, iaproc, naperr, napout, screen, &
         nopts, ptloc, ptnme, grdid, iptint, ptifac
    USE w3servmd, ONLY: extcde
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    USE w3triamd
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)          :: NPT, IMOD
    REAL, INTENT(INOUT)          :: XPT(NPT), YPT(NPT)
    CHARACTER(LEN=40),INTENT(IN) :: PNAMES(NPT)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    LOGICAL                 :: INGRID
    INTEGER                 :: IPT, J, K
    INTEGER                 :: IX1, IY1, IXS, IYS
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    INTEGER                 :: IX(4), IY(4)   ! Indices of points used in interp.
    REAL                    :: RD(4)          ! Interpolation coefficient
    REAL, PARAMETER         :: ACC = 0.05
    REAL                    :: FACTOR
    INTEGER                 :: ITOUT          ! Triangle index in unstructured grids
#ifdef W3_O7a
    INTEGER                 :: IX0, IXN, IY0, IYN, NNX,         &
         KX, KY, JX, IIX, IX2, IY2, IS1
    REAL                    :: RD1, RD2, RDTOT, ZBOX(4), DEPTH
    CHARACTER(LEN=1)         :: SEA(5), LND(5), OUT(5)
    CHARACTER(LEN=9)         :: PARTS
    CHARACTER(LEN=1), ALLOCATABLE :: STRING(:), LINE1(:), LINE2(:)
    !
    DATA sea / ' ', 's', 'e', 'a', ' ' /
    DATA lnd / ' ', 'l', 'n', 'd', ' ' /
    DATA out / ' ', 'x', 'x', 'x', ' ' /
#endif
    !/
#ifdef W3_RTD
    !!  Declare a few temporary variables for rotated grid.  JGLi12Jun2012
    REAL, ALLOCATABLE :: EquLon(:),EquLat(:),StdLon(:),StdLat(:),AnglPT(:)
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3IOPP')
#endif
    !
    IF ( flagll ) THEN
      factor = 1.
    ELSE
      factor = 1.e-3
    END IF
    !
    CALL w3dmo2 ( imod, ndse, ndst, npt )
    grdid  = filext
    !
    nopts  = 0
    !
#ifdef W3_RTD
    !!   Convert standard lon/lat to rotated lon/lat  JGLi12Jun2012
    ALLOCATE( equlon(npt), equlat(npt),                 &
         &          stdlon(npt), stdlat(npt), anglpt(npt) )
 
    stdlon = xpt
    stdlat = ypt
 
    CALL w3lltoeq ( stdlat, stdlon, equlat, equlon,     &
         &                anglpt, polat, polon, npt )
 
    xpt = equlon
    ypt = equlat
 
#endif
    !
    !   Removed by F.A. 2011/04/04  /T      CALL W3GSUP( GSU, NDST )
    !
    ! Loop over output points
    !
    DO ipt=1, npt
      !
#ifdef W3_T
      WRITE (ndst,9010) ipt, xpt(ipt), ypt(ipt), pnames(ipt)
#endif
      !
#ifdef W3_RTD
      !!   Need to wrap rotated Elon values greater than X0.  JGLi12Jun2012
      xpt(ipt) = mod( equlon(ipt)+360.0, 360.0 )
      IF( xpt(ipt) .LT. x0 )  xpt(ipt) = xpt(ipt) + 360.0
#endif
      !
      !     Check if point within grid and compute interpolation weights
      !
      IF (gtype .NE. ungtype) THEN
        ingrid = w3grmp( gsu, xpt(ipt), ypt(ipt), ix, iy, rd )
      ELSE
        CALL is_in_ungrid(imod, dble(xpt(ipt)), dble(ypt(ipt)), itout, ix, iy, rd)
        ingrid = (itout.GT.0)
      END IF
      !
      IF ( .NOT.ingrid ) THEN
        IF ( iaproc .EQ. naperr ) THEN
          IF ( flagll ) THEN
            WRITE (ndse,1000) xpt(ipt), ypt(ipt), pnames(ipt)
          ELSE
            WRITE (ndse,1001) xpt(ipt), ypt(ipt), pnames(ipt)
          END IF
        END IF
        cycle
      END IF
      !
#ifdef W3_T
      DO k = 1,4
        WRITE (ndst,9012) ix(k), iy(k), rd(k)
      END DO
#endif
      !
      !     Check if point not on land
      !
      IF ( mapsta(iy(1),ix(1)) .EQ. 0 .AND. &
           mapsta(iy(2),ix(2)) .EQ. 0 .AND. &
           mapsta(iy(3),ix(3)) .EQ. 0 .AND. &
           mapsta(iy(4),ix(4)) .EQ. 0 ) THEN
        IF ( iaproc .EQ. naperr ) THEN
          IF ( flagll ) THEN
            WRITE (ndse,1002) xpt(ipt), ypt(ipt), pnames(ipt)
          ELSE
            WRITE (ndse,1003) xpt(ipt), ypt(ipt), pnames(ipt)
          END IF
        END IF
        cycle
      END IF
      !
      !     Store interpolation data
      !
      nopts  = nopts + 1
      !
      ptloc(1,nopts) = xpt(ipt)
      ptloc(2,nopts) = ypt(ipt)
#ifdef W3_RTD
      !!   Store the standard lon/lat in PTLOC for output purpose, assuming
      !!   they are not used for any inside calculation.  JGLi12Jun2012
      ptloc(1,nopts) = stdlon(ipt)
      ptloc(2,nopts) = stdlat(ipt)
#endif
      !
      DO k = 1,4
        iptint(1,k,nopts) = ix(k)
        iptint(2,k,nopts) = iy(k)
        ptifac(k,nopts) = rd(k)
      END DO
 
      ptnme(nopts) = pnames(ipt)
      !
    END DO ! End loop over output points (IPT).
    !
#ifdef W3_RTD
    DEALLOCATE( equlon, equlat, stdlon, stdlat, anglpt )
#endif
    !
    ! Diagnostic output
    !
#ifdef W3_O7a
    IF ( iaproc .EQ. napout ) THEN
      WRITE (screen,940) nopts
      DO j=1, nopts
        !
        WRITE (screen,941) ptnme(j), ptloc(:,j)*factor
        ix(:) = iptint(1,:,j)
        iy(:) = iptint(2,:,j)
        rd(:) = ptifac(:,j)
        WRITE (screen,942) (ix(k),iy(k),rd(k),k=1,4)
        !
        zbox   = 0.
        rdtot  = 0.
        DO k = 1,4
          IF ( mapfs(iy(k),ix(k)) .GT. 0 ) THEN
            zbox(k) = zb(ix(k))
            rdtot   = rdtot + rd(k)
          END IF
        END DO
        rdtot  = max( 1.e-7 , rdtot )
        !
        depth  = - ( rd(1)*zbox(1) + &
             rd(2)*zbox(2) + &
             rd(3)*zbox(3) + &
             rd(4)*zbox(4) ) / rdtot
        WRITE (screen,943) depth
        !
        ! *** implementation of O7a option with curvilinear grids is incomplete ***
        !
        IF ( rd1 .LT. 0.05 ) ix2 = ix1
        IF ( rd1 .GT. 0.95 ) ix1 = ix2
        IF ( rd2 .LT. 0.05 ) iy2 = iy1
        IF ( rd2 .GT. 0.95 ) iy1 = iy2
        ix0    = ix1 - 1
        ixn    = ix2 + 1
        iy0    = max( 1 , iy1 - 1 )
        iyn    = min( iy2 + 1 , ny )
        nnx    = 13 * ( ixn - ix0 + 1 )
        !
        ALLOCATE ( string(nnx), line1(nnx), line2(nnx) )
        DO kx=1, nnx
          line1(kx) = ' '
          line2(kx) = '-'
        END DO
        DO kx=7, nnx, 13
          line1(kx) = '|'
          line2(kx) = '+'
        END DO
        !
        IF ( iclose.NE.iclose_none ) THEN
          WRITE (screen,945) (1+mod(kx+nx-1,nx),kx=ix0,ixn)
        ELSE
          WRITE (screen,945) (kx,kx=ix0,ixn)
        END IF
        WRITE (screen,946) line1
        !
        DO ky=iyn, iy0, -1
          !
          string  = line1
          DO kx=ix0, ixn
            IF ( iclose.NE.iclose_none .OR. (kx.GE.1 .AND. kx.LE.nx) ) THEN
              iix    = 1 + mod(kx-1+nx,nx)
              is1    = mapfs(ky,iix)
              IF ( mapsta(ky,iix) .NE. 0 ) THEN
                WRITE (parts,'(F8.1,1X)') -zb(is1)
                nnx    = 2 + (kx-ix0)*13
                DO jx=1, 9
                  string(nnx+jx:nnx+jx) = parts(jx:jx)
                END DO
              ENDIF
            END IF
          END DO
          WRITE (screen,946) string
          !
          string = line2
          DO kx=ix0, ixn
            nnx    = 5 + (kx-ix0)*13
            IF ( iclose.EQ.iclose_none .AND. (kx.LT.1.OR.kx.GT.nx) ) THEN
              string(nnx:nnx+4) = out
            ELSE
              iix    = 1 + mod(kx-1+nx,nx)
              IF ( mapsta(ky,iix) .EQ. 0 ) THEN
                string(nnx:nnx+4) = lnd
              ELSE
                string(nnx:nnx+4) = sea
              END IF
            END IF
          END DO
          WRITE (screen,947) ky, string
          !
          string  = line1
          DO kx=ix0, ixn
            IF ( iclose.NE.iclose_none .OR. (kx.GE.1 .AND. kx.LE.nx) ) THEN
              is1    = mapfs(ky,kx)
              iix    = 1 + mod(kx-1+nx,nx)
              IF ( mapsta(ky,iix) .NE. 0 ) THEN
                WRITE (parts,'(I4,1A,I4)')               &
                     nint(1000.*trnx(ky,iix)),         &
                     '|', nint(1000.*trny(ky,iix))
                nnx    = 2 + (kx-ix0)*13
                DO jx=1, 9
                  string(nnx+jx:nnx+jx) = parts(jx:jx)
                END DO
              ENDIF
            END IF
          END DO
          WRITE (screen,946) string
          WRITE (screen,946) line1
          !
        END DO
        !
        IF ( iclose.NE.iclose_none ) THEN
          WRITE (screen,945) (1+mod(kx+nx-1,nx),kx=ix0,ixn)
        ELSE
          WRITE (screen,945) (kx,kx=ix0,ixn)
        END IF
        DEALLOCATE ( string, line1, line2 )
 
      END DO
      WRITE (screen,*)
      WRITE (screen,*)
    END IF
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_O7a
940 FORMAT (/' Diagnostic output for output points [',i3,'] :'/&
         '--------------------------------------------'/  &
         '    Bottom level in m above grid point'/        &
         '    X/Y transparency in thousands below')
941 FORMAT (/'    Point ',a,' at ',2f8.2,' (degr or km)'/    &
         '    -------------------------------------------------')
942 FORMAT ( '       Interp. cell :',4(' (',2i5,f4.2,')'))
943 FORMAT ( '       Depth (water level = 0)  :',f10.1,' m'/)
945 FORMAT ( '          IX =  ',4i13)
946 FORMAT ( '                     ',52a1)
947 FORMAT ( '          IY =',i5,2x,52a1)
#endif
    !
1000 FORMAT (/' *** WAVEWATCH-III WARNING :'/                   &
         '     OUTPUT POINT OUT OF GRID : ',2f10.3,2x,a/   &
         '     POINT SKIPPPED '/)
1001 FORMAT (/' *** WAVEWATCH-III WARNING :'/                   &
         '     OUTPUT POINT OUT OF GRID : ',2e10.3,2x,a/   &
         '     POINT SKIPPPED '/)
    !
1002 FORMAT (/' *** WAVEWATCH-III WARNING :'/                   &
         '     OUTPUT POINT ON LAND : ',2f10.3,2x,a/       &
         '     POINT SKIPPPED '/)
1003 FORMAT (/' *** WAVEWATCH-III WARNING :'/                   &
         '     OUTPUT POINT ON LAND : ',2e10.3,2x,a/       &
         '     POINT SKIPPPED '/)
    !
#ifdef W3_T
9010 FORMAT (' TEST W3IOPP : INPUT  : ',i4,2f12.2,2x,a)
9011 FORMAT ('               CORR.  :     ',2f12.2)
9012 FORMAT (' TEST W3IOPP : INT. DATA: ',2i6,1f8.2)
9013 FORMAT (' TEST W3IOPP : INT. DATA B): ',4i4,2f8.2)
9020 FORMAT (' TEST W3IOPP : PREPROCESSED DATA',i4,2x,a,2x,2f12.2, &
         4(/'             ',2i5,2f6.3))
9021 FORMAT (' TEST W3IOPP : PREPROCESSED DATA',i4,2x,a,2x,2f12.2, &
         4(/'             ',2i5,f6.3))
#endif
    !/
    !/ End of W3IOPP ----------------------------------------------------- /
    !/
  END SUBROUTINE w3iopp
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3iope ( A )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         12-Jun-2012 |
    !/                  +-----------------------------------+
    !/
    !/    12-Jan-1999 : Distributed FORTRAN 77 version.     ( version 1.18 )
    !/    25-Jan-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/                  Major changes to logistics.
    !/    11-Jun-2001 : Clean-up.                           ( version 2.11 )
    !/    09-Nov-2004 : Multiple grid version.              ( version 3.06 )
    !/    30-Oct-2009 : Implement curvilinear grid type.    ( version 3.14 )
    !/                  (W. E. Rogers & T. J. Campbell, NRL)
    !/    29-Oct-2010 : Implement unstructured grids        ( version 3.14.4 )
    !/                  (A. Roland and F. Ardhuin)
    !/    12-Jun-2012 : Add /RTD option or rotated grid option.
    !/                  (Jian-Guo Li)                       ( version 4.06 )
    !/    01-Mar-2018 : Add option to unrotate spectra      ( version 6.02 )
    !/                  from RTD grid models
    !/    19-Jul-2021 : Momentum and air density support    ( version 7.14 )
    !/
    !  1. Purpose :
    !
    !     Extract point output data and store in output COMMONs. This
    !     action is taken from an earlier version of W3IOPO so that the
    !     point output postprocessor does not need the full sea-point
    !     grid to be able to run.
    !       Note that the output spectrum is F(f,theta). Interpolation
    !     is performed for this spectrum.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       A       R.A.   I   Action spectra on storage grid.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !     - To allow for dynamic ice edges, interpolation factors are
    !       calculated for every time step separately.
    !     - Wind current and depth data are interpolated ignoring ice,
    !       spectrum is interpolated removing ice points.
    !     - Spectra are left in par list to allow for change of shape of
    !       arrays.
    !     - IMOD is not passed to this routine. Since it is used only
    !       in W3WAVE, it is assumed that the pointer are set
    !       appropriately outside this routine.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/SHRD  Switch for shared / distributed memory architecture.
    !     !/DIST  Id.
    !     !/MPI   Switch for message passing method.
    !
    !     !/S     Enable subroutine tracing.
    !     !/T     Test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3gdatmd, ONLY: nk, nth, sig, nx, ny, nsea, nseal,          &
         mapsta, mapfs
#ifdef W3_RTD
    !!   Use spectral rotation sub and angle.  JGLi12Jun2012
    USE w3gdatmd, ONLY: nspec, angld, flagunr
    USE w3servmd, ONLY: w3acturn
#endif
    USE w3wdatmd, ONLY: ice, iceh, icef
#ifdef W3_FLX5
    USE w3wdatmd, ONLY: rhoair
#endif
    USE w3adatmd, ONLY: cg, dw, ua, ud, as, cx, cy,                 &
         sp => sppnt
#ifdef W3_FLX5
    USE w3adatmd, ONLY: taua, tauadir
#endif
    USE w3odatmd, ONLY: ndst, nopts, iptint, ptifac, il, iw, ii,    &
         dpo, wao, wdo, aso, cao, cdo, iceo, iceho,  &
         icefo, spco, naproc
#ifdef W3_FLX5
    USE w3odatmd, ONLY: tauao, taudo, dairo
#endif
#ifdef W3_SETUP
    USE w3wdatmd, ONLY: zeta_setup
    USE w3odatmd, ONLY: zet_seto
#endif
#ifdef W3_MPI
    USE w3odatmd, ONLY: irqpo2
#endif
    USE w3servmd, ONLY: extcde
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
#ifdef W3_T
    USE w3arrymd, ONLY: prt2ds
#endif
    !
    IMPLICIT NONE
    !
#ifdef W3_MPI
    include "mpif.h"
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL, INTENT(IN)        :: A(NTH,NK,0:NSEAL)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: I, IX1, IY1, IX(4), IY(4), J, IS(4), &
         IM(4), IK, ITH, ISP
#ifdef W3_MPI
    INTEGER                 :: IOFF, IERR_MPI
    INTEGER                 :: STAT(MPI_STATUS_SIZE,4*NOPTS)
#endif
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: RD(4), RDS, RDI, FACRD,              &
         WNDX, WNDY, CURX, CURY, FAC1(NK),    &
         FAC2(NK), FAC3(NK), FAC4(NK)
#ifdef W3_FLX5
    REAL                    :: TAUX, TAUY
#endif
    INTEGER                 :: JSEA, ISEA
#ifdef W3_T
    REAL                    :: SPTEST(NK,NTH)
#endif
#ifdef W3_RTD
    REAL :: Spectr(NSPEC), AnglDIS
    INTEGER :: IROT
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3IOPE')
#endif
    !
    cx(0)  = 0.
    cy(0)  = 0.
    !
    ! Loop over spectra -------------------------------------------------- *
    !
    DO i=1, nopts
      !
#ifdef W3_T
      WRITE (ndst,9000) i
#endif
      !
      ! Unpack interpolation data
      !
      ix(:)  = iptint(1,:,i)
      iy(:)  = iptint(2,:,i)
      rd(:)  = ptifac(:,i)
      !
#ifdef W3_T
      !        WRITE (NDST,9001) IX1, IY1, IX(2)
#endif
      !
      !
      ! Correct for land and ice and get sea point counters
      !
      il(i)  = 0
      iw(i)  = 0
      ii(i)  = 0
      rds    = 0.
      rdi    = 0.
      !
      DO j=1, 4
        is(j)  = mapfs(iy(j),ix(j))
        im(j)  = mapsta(iy(j),ix(j))
        IF ( im(j).GT.0 ) THEN
          iw(i)  = iw(i) + 1
          rds    = rds + rd(j)
#ifdef W3_RTD
          irot   = is(j) ! For rotation angle
#endif
        ELSE
          IF ( im(j).LT.0 ) THEN
            ii(i)  = ii(i) + 1
            rdi    = rdi + rd(j)
          ELSE
            il(i)  = il(i) + 1
            rd(j)  = 0.
          END IF
        END IF
      END DO
      !
      ! Depth, wind and current, ignore ice
      !
      IF ( rds+rdi .GT. 1.e-7 ) THEN
        facrd  = 1. / (rds+rdi)
        rd     = rd * facrd
      END IF
      !
#ifdef W3_T
      WRITE (ndst,9002) (is(j),j=1,4), (im(j),j=1,4), (rd(j),j=1,4)
#endif
      !
      ! Interpolate ice depth, wind, stresses, rho air and current
      !
      IF (.NOT. lpdlib) THEN
        icefo(i) = 0
        DO j=1, 4
          isea = mapfs(iy(j),ix(j))
#ifdef W3_DIST
          jsea = 1 + (isea-1)/naproc
#endif
#ifdef W3_SHRD
          jsea = isea
#endif
          icefo(i) = icefo(i) + rd(j)*icef(jsea)
        END DO
      ELSE
        icefo(i) = rd(1)*icef(is(1)) + rd(2)*icef(is(2)) +          &
             rd(3)*icef(is(3)) + rd(4)*icef(is(4))
      END IF
 
      iceo(i) = rd(1)*ice(is(1)) + rd(2)*ice(is(2)) +               &
           rd(3)*ice(is(3)) + rd(4)*ice(is(4))
 
      iceho(i) = rd(1)*iceh(is(1)) + rd(2)*iceh(is(2)) +            &
           rd(3)*iceh(is(3)) + rd(4)*iceh(is(4))
      !
      dpo(i) = rd(1)*dw(is(1)) + rd(2)*dw(is(2)) +                  &
           rd(3)*dw(is(3)) + rd(4)*dw(is(4))
#ifdef W3_SETUP
      dpo(i) = rd(1)*zeta_setup(is(1)) +                     &
           rd(2)*zeta_setup(is(2)) +                     &
           rd(3)*zeta_setup(is(3)) +                     &
           rd(4)*zeta_setup(is(4))
#endif
      !
#ifdef W3_FLX5
      dairo(i) = rd(1)*rhoair(is(1)) + rd(2)*rhoair(is(2)) +        &
           rd(3)*rhoair(is(3)) + rd(4)*rhoair(is(4))
#endif
      !
      wndx   = rd(1) * ua(is(1)) * cos(ud(is(1))) +                 &
           rd(2) * ua(is(2)) * cos(ud(is(2))) +                 &
           rd(3) * ua(is(3)) * cos(ud(is(3))) +                 &
           rd(4) * ua(is(4)) * cos(ud(is(4)))
      wndy   = rd(1) * ua(is(1)) * sin(ud(is(1))) +                 &
           rd(2) * ua(is(2)) * sin(ud(is(2))) +                 &
           rd(3) * ua(is(3)) * sin(ud(is(3))) +                 &
           rd(4) * ua(is(4)) * sin(ud(is(4)))
      !
      wao(i) = sqrt( wndx**2 + wndy**2 )
      IF ( wao(i).GT.1.e-7 ) THEN
        wdo(i) = atan2(wndy,wndx)
#ifdef W3_RTD
        IF ( flagunr ) wdo(i) = wdo(i) - angld(is(1))*dera
#endif
      ELSE
        wdo(i) = 0.
      END IF
      !
#ifdef W3_FLX5
      taux   = rd(1) * taua(is(1)) * cos(tauadir(is(1))) +          &
           rd(2) * taua(is(2)) * cos(tauadir(is(2))) +          &
           rd(3) * taua(is(3)) * cos(tauadir(is(3))) +          &
           rd(4) * taua(is(4)) * cos(tauadir(is(4)))
      tauy   = rd(1) * taua(is(1)) * sin(tauadir(is(1))) +          &
           rd(2) * taua(is(2)) * sin(tauadir(is(2))) +          &
           rd(3) * taua(is(3)) * sin(tauadir(is(3))) +          &
           rd(4) * taua(is(4)) * sin(tauadir(is(4)))
      !
      tauao(i) = sqrt( taux**2 + tauy**2 )
      IF ( tauao(i).GT.1.e-7 ) THEN
        taudo(i) = atan2(tauy,taux)
#ifdef W3_RTD
        IF ( flagunr ) taudo(i) = taudo(i) - angld(is(1))*dera
#endif
      ELSE
        taudo(i) = 0.
      END IF
      !
#endif
      aso(i) = rd(1)*as(is(1)) + rd(2)*as(is(2)) +                  &
           rd(3)*as(is(3)) + rd(4)*as(is(4))
      !
      curx   = rd(1)*cx(is(1)) + rd(2)*cx(is(2)) +                  &
           rd(3)*cx(is(3)) + rd(4)*cx(is(4))
      cury   = rd(1)*cy(is(1)) + rd(2)*cy(is(2)) +                  &
           rd(3)*cy(is(3)) + rd(4)*cy(is(4))
      !
      cao(i) = sqrt( curx**2 + cury**2 )
      IF ( cao(i).GT.1.e-7 ) THEN
        cdo(i) = atan2(cury,curx)
#ifdef W3_RTD
        IF ( flagunr ) cdo(i) = cdo(i) - angld(is(1))*dera
#endif
      ELSE
        cdo(i) = 0.
      END IF
      !
      ! Interp. weights for spectra, no ice points (spectra by def. zero)
      !
      IF ( rds .GT. 1.e-7 ) THEN
        facrd  = (rds+rdi) / rds
        rd     = rd * facrd
      END IF
      !
#ifdef W3_T
      WRITE (ndst,9003) (rd(j),j=1,4)
#endif
      !
      ! Extract spectra, shared memory version
      !        (done in separate step for MPP compatibility)
      !
#ifdef W3_SHRD
      DO j=1, 4
        DO ik=1, nk
          DO ith=1, nth
            sp(ith,ik,j) = a(ith,ik,is(j))
          END DO
        END DO
      END DO
#endif
      !
      ! Extract spectra, distributed memory version(s)
      !
#ifdef W3_MPI
      ioff   = 1 + 4*(i-1)
      CALL mpi_startall ( 4, irqpo2(ioff), ierr_mpi )
      CALL mpi_waitall  ( 4, irqpo2(ioff), stat, ierr_mpi )
#endif
      !
      ! Interpolate spectrum
      !
      DO ik=1, nk
        fac1(ik) = tpi * sig(ik) / cg(ik,is(1))
        fac2(ik) = tpi * sig(ik) / cg(ik,is(2))
        fac3(ik) = tpi * sig(ik) / cg(ik,is(3))
        fac4(ik) = tpi * sig(ik) / cg(ik,is(4))
      END DO
      !
      DO ik=1,nk
        DO ith=1,nth
          isp    = ith + (ik-1)*nth
          spco(isp,i) = rd(1) * sp(ith,ik,1) * fac1(ik)             &
               + rd(2) * sp(ith,ik,2) * fac2(ik)             &
               + rd(3) * sp(ith,ik,3) * fac3(ik)             &
               + rd(4) * sp(ith,ik,4) * fac4(ik)
#ifdef W3_T
          sptest(ik,ith) = spco(isp,i)
#endif
        END DO
      END DO
      !
#ifdef W3_RTD
      !!  Rotate the interpolated spectrum by -AnglD(IS(1)).  JGLi12Jun2012
      !!  SPCO still holds action not energy spectrum yet.  JGLi18Jun2013
      !!  Use new index IROT rather than IS(1) as in some cases
      !!  IS(1) will be a coast point and have an index of 0. C.Bunney 15/02/2011
      IF ( flagunr ) THEN
        spectr = spco(:,i)
        angldis = -angld(irot)
        CALL  w3acturn( nth, nk, angldis, spectr )
        spco(:,i) = spectr
      END IF
 
#endif
      !
#ifdef W3_T
      WRITE (ndst,9004) dpo(i), wao(i), wdo(i)*rade,             &
           cao(i), cdo(i)*rade
#endif
 
      ! FA COMMENTED OUT: BUG
      !At line 1974 of file w3arrymd.f90
      !Fortran runtime error: Index '52' of dimension 1 of array 'pnum2' above upper bound of 51
#ifdef W3_T
      ! CALL PRT2DS (NDST, NK, NK, NTH, SPTEST, SIG(1:), ' ', 1.,0.,&
      !              0.0001, 'E(f,theta)', 'm**2s', 'TEST OUTPUT' )
#endif
      !
    END DO
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST W3IOPE : POINT NR.:',i3)
9001 FORMAT (' TEST W3IOPE :',2i8,'   (',i3,')')
9002 FORMAT (' TEST W3IOPE :',4i7,2x,4i2,2x,4f5.2)
9003 FORMAT (' TEST W3IOPE :',40x,4f5.2)
9004 FORMAT (' TEST W3IOPE :',f8.1,2(f7.2,f7.1))
#endif
    !/
    !/ End of W3IOPE ----------------------------------------------------- /
    !/
  END SUBROUTINE w3iope
 
#ifdef W3_BIN2NC
 
  integer function nf90_err_check(errcode, ILINE)
    USE netcdf
    USE w3odatmd, ONLY: ndse
    implicit none
    integer, intent(in) :: errcode, iline
    
    nf90_err_check = errcode
    if(errcode /= nf90_noerr) then
      WRITE(ndse,*) ' *** WAVEWATCH III ERROR IN W3IOPO :'
      WRITE(ndse,*) ' LINE NUMBER ', iline
      WRITE(ndse,*) ' NETCDF ERROR MESSAGE: '
      WRITE(ndse,*) 'Error: ', trim(nf90_strerror(errcode))
      return 
    endif
  end function nf90_err_check
 
  SUBROUTINE w3iopon_read(IOTST, IMOD_IN, filename, ncerr)
    USE netcdf
    USE w3odatmd, ONLY: w3dmo2
    USE w3wdatmd, ONLY: time
    USE w3gdatmd, ONLY: nth, nk, nspec, filext
    USE w3odatmd, ONLY: ndst, ndse, ipass => ipass2, nopts, iptint, &
         il, iw, ii, ptloc, ptifac, dpo, wao, wdo,   &
         aso, cao, cdo, spco, ptnme, o2init, fnmpre, &
         grdid, iceo, iceho, icefo, w3dmo2
    USE w3servmd, ONLY: extcde
#ifdef W3_FLX5
    USE w3odatmd, ONLY: tauao, taudo, dairo
#endif
#ifdef W3_SETUP
    USE w3odatmd, ONLY: zet_seto
#endif
    IMPLICIT NONE
 
    INTEGER, INTENT(OUT)          :: IOTST
    INTEGER, INTENT(IN), OPTIONAL :: IMOD_IN
    character(*), intent(in) :: filename
    integer, intent(inout) :: ncerr
    INTEGER :: IGRD,MK,MTH
    integer :: fh
    integer :: d_nopts, d_nspec, d_vsize, d_namelen, d_grdidlen, d_time, d_ww3time
    integer :: d_nopts_len, d_nspec_len, d_vsize_len, d_namelen_len, d_grdidlen_len, d_time_len, d_ww3time_len
    integer :: v_idtst, v_vertst, v_nk, v_nth, v_ptloc, v_ptnme, v_time, v_ww3time 
    integer :: v_dpo, v_wao, v_wdo
#ifdef W3_FLX5
    integer :: v_tauao,v_taudo, v_dairo
#endif
#ifdef W3_SETUP
    integer :: v_zet_seto
#endif
    integer :: v_aso, v_cao, v_cdo, v_iceo
    integer :: v_iceho, v_icefo, v_grdid, v_spco
    integer :: v_title_len, v_version_len 
    CHARACTER(LEN=31)       :: IDTST
    CHARACTER(LEN=10)       :: VERTST
 
    iotst = 0
 
    IF (PRESENT(imod_in)) THEN
      igrd = imod_in
    ELSE
      igrd = 1
    END IF
 
    ! Open the netCDF file.
    ncerr = nf90_open(filename, nf90_nowrite, fh)
    if (nf90_err(ncerr) .ne. 0) return
 
    ! Read and check the version:
    ncerr = nf90_inquire_attribute(fh, nf90_global, 'title', len = v_title_len) 
    if (nf90_err(ncerr) .ne. 0) return
    ncerr = nf90_get_att(fh, nf90_global, 'title', idtst)
    if (nf90_err(ncerr) .ne. 0) return
    ncerr = nf90_inquire_attribute(fh, nf90_global, 'version', len = v_version_len)
    if (nf90_err(ncerr) .ne. 0) return
    ncerr = nf90_get_att(fh, nf90_global, 'version', vertst)
    if (nf90_err(ncerr) .ne. 0) return
 
    IF ( idtst .NE. idstr ) THEN
      WRITE (ndse,902) idtst, idstr
      CALL extcde ( 10 )
    END IF
    IF ( vertst .NE. veropt ) THEN
      WRITE (ndse,903) vertst, veropt
      CALL extcde ( 11 )
    END IF
 
    ! Read the dimension information for NOPTS.
    ncerr = nf90_inq_dimid(fh, dname_nopts, d_nopts)
    if (nf90_err(ncerr) .ne. 0) return 
    ncerr = nf90_inquire_dimension(fh, d_nopts, len = d_nopts_len)
    if (nf90_err(ncerr) .ne. 0) return
    nopts=d_nopts_len
 
    ! Read the dimension information for NSPEC.
    ncerr = nf90_inq_dimid(fh, dname_nspec, d_nspec)
    if (nf90_err(ncerr) .ne. 0) return
    ncerr = nf90_inquire_dimension(fh, d_nspec, len = d_nspec_len)
    if (nf90_err(ncerr) .ne. 0) return
 
    ! Read the dimension information for VSIZE.
    ncerr = nf90_inq_dimid(fh, dname_vsize, d_vsize)
    if (nf90_err(ncerr) .ne. 0) return
    ncerr = nf90_inquire_dimension(fh, d_vsize, len = d_vsize_len)
    if (nf90_err(ncerr) .ne. 0) return
 
    ! Read the dimension information for NAMELEN.
    ncerr = nf90_inq_dimid(fh, dname_namelen, d_namelen)
    if (nf90_err(ncerr) .ne. 0) return
    ncerr = nf90_inquire_dimension(fh, d_namelen, len = d_namelen_len)
    if (nf90_err(ncerr) .ne. 0) return
 
    ! Read the dimension information for GRDIDLEN.
    ncerr = nf90_inq_dimid(fh, dname_grdidlen, d_grdidlen)
    if (nf90_err(ncerr) .ne. 0) return
    ncerr = nf90_inquire_dimension(fh, d_grdidlen, len = d_grdidlen_len)
    if (nf90_err(ncerr) .ne. 0) return
 
    ! Read the dimention information from time
    ncerr = nf90_inq_dimid(fh, dname_time, d_time)
    if (nf90_err(ncerr) .ne. 0) return
    ncerr = nf90_inquire_dimension(fh, d_time, len = d_time_len)
    if (nf90_err(ncerr) .ne. 0) return
 
    IF ( ipass .LE. d_time_len ) THEN 
 
      IF ( ipass.EQ.1 ) THEN
 
        ! Read scalar variables.
        ncerr = nf90_inq_varid(fh, vname_nk, v_nk)
        if (nf90_err(ncerr) .ne. 0) return
        ncerr = nf90_get_var(fh, v_nk, mk)
        if (nf90_err(ncerr) .ne. 0) return
        ncerr = nf90_inq_varid(fh, vname_nth, v_nth)
        if (nf90_err(ncerr) .ne. 0) return
        ncerr = nf90_get_var(fh, v_nth, mth)
        if (nf90_err(ncerr) .ne. 0) return
 
        !read in written variables NK, NTH as MK and MTH
        !and ensure they match 
        IF (nk.NE.mk .OR. nth.NE.mth) THEN
          WRITE (ndse,904) mk, mth, nk, nth
          CALL extcde ( 12 )
        END IF
 
        ! Allocate variables: 
        IF ( .NOT. o2init )                                     &
           CALL w3dmo2 ( igrd, ndse, ndst, nopts )
 
        ! Read vars with nopts as a dimension.
        ncerr = nf90_inq_varid(fh, vname_ptloc, v_ptloc)
        if (nf90_err(ncerr) .ne. 0) return
        ncerr = nf90_get_var(fh, v_ptloc, ptloc, start = (/ 1, 1/), &
          count = (/ d_vsize_len, d_nopts_len /)) 
        if (nf90_err(ncerr) .ne. 0) return
        ncerr = nf90_inq_varid(fh, vname_ptnme, v_ptnme)
        if (nf90_err(ncerr) .ne. 0) return
        ncerr = nf90_get_var(fh, v_ptnme, ptnme)
        if (nf90_err(ncerr) .ne. 0) return
      END IF
 
      !Variables read based on time (IPASS): 
 
      ncerr = nf90_inq_varid(fh, vname_ww3time, v_ww3time)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_get_var(fh, v_ww3time, time, start = (/ 1, ipass/), &
          count = (/ d_vsize_len, 1 /))
      if (nf90_err(ncerr) .ne. 0) return
 
      ! set IW, II and IL to 0,
      ! These values are set to 0 in binary file and have been removed 
      ! from netcdf file.  Possible can be completely removed. 
      iw = 0
      ii = 0
      il = 0
 
      ncerr = nf90_inq_varid(fh, vname_dpo, v_dpo)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_get_var(fh, v_dpo, dpo, start = (/ 1, ipass/), &
          count = (/ nopts, 1 /))
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_inq_varid(fh, vname_wao, v_wao)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_get_var(fh, v_wao, wao, start = (/ 1, ipass/), &
          count = (/ nopts, 1 /))
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_inq_varid(fh, vname_wdo, v_wdo)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_get_var(fh, v_wdo, wdo, start = (/ 1, ipass/), &
          count = (/ nopts, 1 /))
      if (nf90_err(ncerr) .ne. 0) return
#ifdef W3_FLX5
      ncerr = nf90_inq_varid(fh, vname_tauao, v_tauao)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_get_var(fh, v_tauao, tauao, start = (/ 1, ipass/), &
          count = (/ nopts, 1 /))
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_inq_varid(fh, vname_taudo, v_taudo)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_get_var(fh, v_taudo, taudo, start = (/ 1, ipass/), &
          count = (/ nopts, 1 /))
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_inq_varid(fh, vname_dairo, v_dairo)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_get_var(fh, v_dairo, dairo, start = (/ 1, ipass/), &
          count = (/ nopts, 1 /))
      if (nf90_err(ncerr) .ne. 0) return
#endif
#ifdef W3_SETUP
      ncerr = nf90_inq_varid(fh, zet_seto, v_zet_seto)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_get_var(fh, v_zet_seto, zet_seto, start = (/ 1, ipass/), &
          count = (/ nopts, 1 /))
      if (nf90_err(ncerr) .ne. 0) return
#endif
      ncerr = nf90_inq_varid(fh, vname_aso, v_aso)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_get_var(fh, v_aso, aso, start = (/ 1, ipass/), &
          count = (/ nopts, 1 /))
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_inq_varid(fh, vname_cao, v_cao)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_get_var(fh, v_cao, cao, start = (/ 1, ipass/), &
          count = (/ nopts, 1 /))
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_inq_varid(fh, vname_cdo, v_cdo)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_get_var(fh, v_cdo, cdo, start = (/ 1, ipass/), &
          count = (/ nopts, 1 /))
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_inq_varid(fh, vname_iceo, v_iceo)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_get_var(fh, v_iceo, iceo, start = (/ 1, ipass/), &
          count = (/ nopts, 1 /))
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_inq_varid(fh, vname_iceho, v_iceho)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_get_var(fh, v_iceho, iceho, start = (/ 1, ipass/), &
          count = (/ nopts, 1 /))
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_inq_varid(fh, vname_icefo, v_icefo)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_get_var(fh, v_icefo, icefo, start = (/ 1, ipass/), &
          count = (/ nopts, 1 /))
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_inq_varid(fh, vname_grdid, v_grdid)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_get_var(fh, v_grdid, grdid, start = (/ 1, 1, ipass/), &
          count = (/ 13, nopts, 1 /))
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_inq_varid(fh, vname_spco, v_spco)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_get_var(fh, v_spco, spco, start = (/ 1, 1, ipass/), &
          count = (/nspec, nopts, 1 /))
      if (nf90_err(ncerr) .ne. 0) return
 
    ELSE 
      ! Set flag to indicate IPASS > d_time_len 
      ! and are at the end of the
      iotst  = -1
    END IF 
 
    ! Close the file.
    ncerr = nf90_close(fh)
    if (nf90_err(ncerr) .ne. 0) return
 
902 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOPON :'/               &
         '     ILEGAL IDSTR, READ : ',a/                        &
         '                  CHECK : ',a/)
903 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOPON :'/               &
         '     ILEGAL VEROPT, READ : ',a/                       &
         '                   CHECK : ',a/)
904 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOPO :'/                &
         '     ERROR IN SPECTRA, MK, MTH : ',2i8/               &
         '              ARRAY DIMENSIONS : ',2i8/)
 
 
  END SUBROUTINE w3iopon_read
 
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3iopon_write(timestep_only,filename, ncerr)
    USE netcdf 
    USE w3gdatmd, ONLY: nth, nk, nspec
    USE w3wdatmd, ONLY: time
    USE w3odatmd, ONLY: ndst, ndse, ipass => ipass2, nopts, iptint, &
         ptloc, ptifac, dpo, wao, wdo,   &
         aso, cao, cdo, spco, ptnme, o2init, fnmpre, &
         grdid, iceo, iceho, icefo
  USE w3timemd, ONLY: caltype, t2d, u2d, tsub
#ifdef W3_FLX5
    USE w3odatmd, ONLY: tauao, taudo, dairo
#endif
#ifdef W3_SETUP
    USE w3odatmd, ONLY: zet_seto
#endif
 
    IMPLICIT NONE
    integer, intent(in) :: timestep_only ! 1 if only timestep should be written.
    character(*), intent(in) :: filename
    integer, intent(inout) :: ncerr
    integer :: ndim, nvar, fmt, itime, fh
    integer :: d_nopts, d_nspec, d_vsize, d_namelen, d_grdidlen, d_time
    integer :: v_idtst, v_vertst, v_nk, v_nth, v_ptloc, v_ptnme, v_time, v_ww3time
    integer :: v_dpo, v_wao, v_wdo
#ifdef W3_FLX5
    integer :: v_tauao, v_taudo, v_dairo
#endif    
#ifdef W3_SETUP
    integer :: v_zet_seto
#endif  
    integer :: v_aso, v_cao, v_cdo, v_iceo
    integer :: v_iceho, v_icefo, v_grdid, v_spco
    integer :: curdate(8), refdate(8),ierr
    double precision :: outjulday
 
    !If first pass, or if you are writting a file for every time-step: 
    IF ( ipass.EQ.1  .OR. timestep_only.EQ.1 ) THEN 
      ! Create the netCDF file.
      ncerr = nf90_create(filename, nf90_netcdf4, fh)
      if (nf90_err(ncerr) .ne. 0) return
 
      ! Define dimensions.
      ncerr = nf90_def_dim(fh, dname_nopts, nopts, d_nopts)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_def_dim(fh, dname_nspec, nspec, d_nspec)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_def_dim(fh, dname_vsize, 2, d_vsize)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_def_dim(fh, dname_namelen, 40, d_namelen)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_def_dim(fh, dname_grdidlen, 13, d_grdidlen)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_def_dim(fh, dname_time, nf90_unlimited, d_time)
      if (nf90_err(ncerr) .ne. 0) return
 
      ! Define global attributes.
      ncerr = nf90_put_att(fh, nf90_global, 'title', idstr)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_put_att(fh, nf90_global, 'version', veropt)
      if (nf90_err(ncerr) .ne. 0) return
 
      ! Define scalar variables.
      ncerr = nf90_def_var(fh, vname_nk, nf90_int, v_nk)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_def_var(fh, vname_nth, nf90_int, v_nth)
      if (nf90_err(ncerr) .ne. 0) return
 
      ! Define vars with nopts as a dimension. Point location and name
      ncerr = nf90_def_var(fh, vname_ptloc, nf90_float, (/d_vsize, d_nopts/), v_ptloc)
      if (nf90_err(ncerr) .ne. 0) return 
      ncerr = nf90_def_var(fh, vname_ptnme, nf90_char, (/d_namelen, d_nopts/), v_ptnme)
      if (nf90_err(ncerr) .ne. 0) return
 
      ! Define time for each time step 
      ncerr = nf90_def_var(fh, vname_ww3time, nf90_int, (/d_vsize, d_time/),v_ww3time)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_def_var(fh, vname_time, nf90_double, (/d_time/),v_time)
      if (nf90_err(ncerr) .ne. 0) return
      SELECT CASE (trim(caltype))
      CASE ('360_day')
        ncerr = nf90_put_att(fh, v_time, 'long_name', 'time in 360 day calendar')
        if (nf90_err(ncerr) .ne. 0) return
      CASE ('365_day')
        ncerr = nf90_put_att(fh, v_time, 'long_name', 'time in 365 day calendar') 
        if (nf90_err(ncerr) .ne. 0) return
      CASE ('standard')
        ncerr = nf90_put_att(fh, v_time, 'long_name', 'Julian day (UT)') 
        if (nf90_err(ncerr) .ne. 0) return
      END SELECT
      ncerr = nf90_put_att(fh, v_time, 'standard_name', 'time')
      if (nf90_err(ncerr) .ne. 0) return 
      ncerr = nf90_put_att(fh, v_time, 'units', 'days since 1990-01-01 00:00:00')
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_put_att(fh, v_time, 'conventions','Relative Julian days with decimal part (as parts of the day)')
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_put_att(fh, v_time, 'axis', 'T')
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_put_att(fh, v_time, 'calendar', trim(caltype))
      if (nf90_err(ncerr) .ne. 0) return
 
      ! Define vars with nopts and time as dimensions 
      ncerr = nf90_def_var(fh, vname_dpo, nf90_float, (/d_nopts, d_time/), v_dpo)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_def_var(fh, vname_wao, nf90_float, (/d_nopts, d_time/), v_wao)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_def_var(fh, vname_wdo, nf90_float, (/d_nopts, d_time/), v_wdo)
      if (nf90_err(ncerr) .ne. 0) return
#ifdef W3_FLX5
      ncerr = nf90_def_var(fh, vname_tauao, nf90_float, (/d_nopts, d_time/), v_tauao)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_def_var(fh, vname_taudo, nf90_float, (/d_nopts, d_time/), v_taudo)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_def_var(fh, vname_dairo, nf90_float, (/d_nopts, d_time/), v_dairo)
      if (nf90_err(ncerr) .ne. 0) return
#endif    
#ifdef W3_SETUP
      ncerr = nf90_def_var(fh, vname_zet_seto, nf90_float, (/d_nopts, d_time/), v_zet_seto)
      if (nf90_err(ncerr) .ne. 0) return
#endif    
      ncerr = nf90_def_var(fh, vname_aso, nf90_float, (/d_nopts, d_time/), v_aso)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_def_var(fh, vname_cao, nf90_float, (/d_nopts, d_time/), v_cao)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_def_var(fh, vname_cdo, nf90_float, (/d_nopts, d_time/), v_cdo)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_def_var(fh, vname_iceo, nf90_float, (/d_nopts, d_time/), v_iceo)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_def_var(fh, vname_iceho, nf90_float, (/d_nopts, d_time/), v_iceho)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_def_var(fh, vname_icefo, nf90_float, (/d_nopts, d_time/), v_icefo)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_def_var(fh, vname_grdid, nf90_char, (/d_grdidlen, d_nopts, d_time/), v_grdid)
      if (nf90_err(ncerr) .ne. 0) return
      
      ! Define spectral output with dimensions nspec, nopts and time
      ncerr = nf90_def_var(fh, vname_spco, nf90_float, (/d_nspec, d_nopts, d_time/), v_spco)
      if (nf90_err(ncerr) .ne. 0) return
  
      ! End of all variable definitions 
      ncerr = nf90_enddef(fh)
      if (nf90_err(ncerr) .ne. 0) return 
 
      ! Write the scalar data.
      ncerr = nf90_put_var(fh, v_nk, nk)
      if (nf90_err(ncerr) .ne. 0) return
      ncerr = nf90_put_var(fh, v_nth, nth)
      if (nf90_err(ncerr) .ne. 0) return
 
      ! Write the data with NOPTS as a dimension. (no time dimension)
      if (associated(ptloc)) then
         ncerr = nf90_put_var(fh, v_ptloc, ptloc(:,1:nopts))
         if (nf90_err(ncerr) .ne. 0) return
      endif
      if (associated(ptnme)) then
         ncerr = nf90_put_var(fh, v_ptnme, ptnme(1:nopts))
         if (nf90_err(ncerr) .ne. 0) return
      endif
      
    ELSE 
      ! If we are writing to the same file, re-open the file 
      ncerr = nf90_open(filename, nf90_write, fh)
      if (nf90_err(ncerr) .ne. 0) return
    END IF 
 
    !Determine the start for the time dimension 
    IF ( timestep_only.EQ.1 ) THEN
       itime=1
    ELSE 
       itime=ipass
    END IF
 
    ! Write Time
    IF ( itime > 1 ) THEN
       ncerr = nf90_inq_varid(fh, vname_ww3time, v_ww3time)
       if (nf90_err(ncerr) .ne. 0) return
       ncerr = nf90_inq_varid(fh, vname_time, v_time)
       if (nf90_err(ncerr) .ne. 0) return
    END IF
    ncerr = nf90_put_var(fh, v_ww3time, time, start = (/ 1, itime/), &
       count = (/ 2, 1 /))
    if (nf90_err(ncerr) .ne. 0) return
 
    CALL u2d('days since 1990-01-01 00:00:00',refdate,ierr)
    CALL t2d(time,curdate,ierr)
    outjulday=tsub(refdate,curdate)
 
    ncerr = nf90_put_var(fh, v_time, outjulday, start = (/itime/))
    if (nf90_err(ncerr) .ne. 0) return
 
 
    ! If itime > 1 need to inquire varid 
    IF ( itime > 1 ) THEN 
       ncerr = nf90_inq_varid(fh, vname_dpo, v_dpo)
       if (nf90_err(ncerr) .ne. 0) return
       ncerr = nf90_inq_varid(fh, vname_wao, v_wao)
       if (nf90_err(ncerr) .ne. 0) return
       ncerr = nf90_inq_varid(fh, vname_wdo, v_wdo)
       if (nf90_err(ncerr) .ne. 0) return
#ifdef W3_FLX5
       ncerr = nf90_inq_varid(fh, vname_tauao, v_tauao)
       if (nf90_err(ncerr) .ne. 0) return
       ncerr = nf90_inq_varid(fh, vname_taudo, v_taudo)
       if (nf90_err(ncerr) .ne. 0) return
       ncerr = nf90_inq_varid(fh, vname_dairo, v_dairo)
       if (nf90_err(ncerr) .ne. 0) return
#endif
#ifdef W3_SETUP
       ncerr = nf90_inq_varid(fh, vname_zet_seto, v_zet_seto)
       if (nf90_err(ncerr) .ne. 0) return
#endif
       ncerr = nf90_inq_varid(fh, vname_aso, v_aso)
       if (nf90_err(ncerr) .ne. 0) return
       ncerr = nf90_inq_varid(fh, vname_cao, v_cao)
       if (nf90_err(ncerr) .ne. 0) return
       ncerr = nf90_inq_varid(fh, vname_cdo, v_cdo)
       if (nf90_err(ncerr) .ne. 0) return
       ncerr = nf90_inq_varid(fh, vname_iceo, v_iceo)
       if (nf90_err(ncerr) .ne. 0) return
       ncerr = nf90_inq_varid(fh, vname_iceho, v_iceho)
       if (nf90_err(ncerr) .ne. 0) return
       ncerr = nf90_inq_varid(fh, vname_icefo, v_icefo)
       if (nf90_err(ncerr) .ne. 0) return
       ncerr = nf90_inq_varid(fh, vname_grdid, v_grdid)
       if (nf90_err(ncerr) .ne. 0) return
       ncerr = nf90_inq_varid(fh, vname_spco, v_spco)
       if (nf90_err(ncerr) .ne. 0) return
    END IF
 
    ncerr = nf90_put_var(fh, v_dpo, dpo, start = (/ 1, itime/), &
       count = (/ nopts, 1 /))
    if (nf90_err(ncerr) .ne. 0) return
 
    ncerr = nf90_put_var(fh, v_wao, wao, start = (/ 1, itime/), &
       count = (/ nopts, 1 /))
    if (nf90_err(ncerr) .ne. 0) return
 
    ncerr = nf90_put_var(fh, v_wdo, wdo, start = (/ 1, itime/), &
       count = (/ nopts, 1 /))
    if (nf90_err(ncerr) .ne. 0) return
 
#ifdef W3_FLX5
    ncerr = nf90_put_var(fh, v_tauao, tauao, start = (/ 1, itime/), &
       count = (/ nopts, 1 /))
    if (nf90_err(ncerr) .ne. 0) return
 
    ncerr = nf90_put_var(fh, v_taudo, taudo, start = (/ 1, itime/), &
       count = (/ nopts, 1 /))
    if (nf90_err(ncerr) .ne. 0) return
 
    ncerr = nf90_put_var(fh, v_dairo, dairo, start = (/ 1, itime/), &
       count = (/ nopts, 1 /))
    if (nf90_err(ncerr) .ne. 0) return
#endif
#ifdef W3_SETUP
    ncerr = nf90_put_var(fh, v_zet_seto, zet_seto, start = (/ 1, itime/), &
       count = (/ nopts, 1 /))
    if (nf90_err(ncerr) .ne. 0) return
#endif
    ncerr = nf90_put_var(fh, v_aso, aso, start = (/ 1, itime/), &
       count = (/ nopts, 1 /))
    if (nf90_err(ncerr) .ne. 0) return
 
    ncerr = nf90_put_var(fh, v_cao, cao, start = (/ 1, itime/), &
       count = (/ nopts, 1 /))
    if (nf90_err(ncerr) .ne. 0) return
 
    ncerr = nf90_put_var(fh, v_cdo, cdo, start = (/ 1, itime/), &
       count = (/ nopts, 1 /))
    if (nf90_err(ncerr) .ne. 0) return
 
    ncerr = nf90_put_var(fh, v_iceo, iceo, start = (/ 1, itime/), &
       count = (/ nopts, 1 /))
    if (nf90_err(ncerr) .ne. 0) return
 
    ncerr = nf90_put_var(fh, v_iceho, iceho, start = (/ 1, itime/), &
       count = (/ nopts, 1 /))
    if (nf90_err(ncerr) .ne. 0) return
 
    ncerr = nf90_put_var(fh, v_icefo, icefo, start = (/ 1, itime/), &
       count = (/ nopts, 1 /))
    if (nf90_err(ncerr) .ne. 0) return
 
    ncerr = nf90_put_var(fh, v_grdid, grdid, start = (/ 1, 1, itime/), &
       count = (/ 13, nopts, 1 /))
    if (nf90_err(ncerr) .ne. 0) return
 
    !write spectral output
    ncerr = nf90_put_var(fh, v_spco, spco, start = (/ 1, 1, itime/), &
       count = (/nspec, nopts, 1 /))
    if (nf90_err(ncerr) .ne. 0) return
 
    ! Close the file.
    ncerr = nf90_close(fh)
    if (nf90_err(ncerr) .ne. 0) return
 
  END SUBROUTINE w3iopon_write
 
  SUBROUTINE w3iopon ( INXOUT, NDSOP, IOTST, IMOD)
    USE w3gdatmd, ONLY: w3setg
    USE w3wdatmd, ONLY: w3setw
    USE w3odatmd, ONLY: w3seto
    USE w3gdatmd, ONLY: filext
    USE w3wdatmd, ONLY: time
    USE w3odatmd, ONLY: ndst, ndse, ipass => ipass2, fnmpre 
    USE w3odatmd, ONLY: ofiles
    USE w3servmd, ONLY: extcde
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    use netcdf
    IMPLICIT NONE
 
    CHARACTER, INTENT(IN)         :: INXOUT*(*)
    INTEGER, INTENT(IN)           :: NDSOP
    INTEGER, INTENT(OUT)          :: IOTST
    INTEGER, INTENT(IN), OPTIONAL :: IMOD
 
    CHARACTER(LEN=15) :: TIMETAG
    INTEGER :: IGRD
    character(len = 124) :: filename
    integer :: ncerr
 
#ifdef W3_S
    CALL strace (ient, 'W3IOPON')
#endif
 
    ! IPASS essentially is the time variable dimension 
    ipass  = ipass + 1
 
    ! Optimistically assume success.
    iotst = 0
 
    ! Has a model number been specified?
    IF (PRESENT(imod)) THEN
      igrd = imod
    ELSE
      igrd = 1
    END IF
 
    CALL w3seto(igrd, ndse, ndst)
    CALL w3setg(igrd, ndse, ndst)
    CALL w3setw(igrd, ndse, ndst)
 
    ! INXOUT must be 'READ' or 'WRITE'.
    IF (inxout .NE. 'READ' .AND. inxout .NE. 'WRITE') THEN
      WRITE (ndse, 900) inxout
      CALL extcde(1)
    END IF
 
    ! Determine filename.
    IF ( ofiles(2) .EQ. 1 ) THEN 
      ! Create TIMETAG for file name using YYYYMMDD.HHMMS prefix
      WRITE(timetag,"(i8.8,'.'i6.6)")time(1),time(2)
      filename = fnmpre(:len_trim(fnmpre))//timetag//'.out_pnt.'//filext(:len_trim(filext))//'.nc'
    ELSE 
      filename = fnmpre(:len_trim(fnmpre))//'out_pnt.'//filext(:len_trim(filext))//'.nc'
    END IF 
 
    ! Do a read or a write of the point file.
    IF (inxout .EQ. 'READ') THEN
      CALL w3iopon_read(iotst, imod, filename, ncerr)
    ELSE
      CALL w3iopon_write(ofiles(2), filename, ncerr)
    ENDIF
    if (nf90_err(ncerr) .ne. 0) then
      WRITE(ndse,*) ' *** WAVEWATCH III ERROR IN W3IOPO :'
      WRITE(ndse,*) 'Nonzero return at end of W3IOPON'
      WRITE(ndse,*) 'Error: ', trim(nf90_strerror(ncerr))
      CALL extcde(21)
    endif
 
    !/
    !/ End of W3IOPON ----------------------------------------------------- /
    !/
 
900 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOPO :'/                &
         '     ILEGAL INXOUT VALUE: ',a/)
  END SUBROUTINE w3iopon
#endif
  !/ ------------------------------------------------------------------- /
#ifdef W3_ASCII
 
#endif
 
  SUBROUTINE w3iopo ( INXOUT, NDSOP, IOTST, IMOD &
#ifdef W3_ASCII
      ,NDSOA &
#endif
      )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         25-Jul-2006 |
    !/                  +-----------------------------------+
    !/
    !/    07-Jan-1999 : Distributed FORTRAN 77 version.     ( version 1.18 )
    !/    30-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/                  Major changes to logistics.
    !/    10-Nov-2004 : Multiple grid version.              ( version 3.06 )
    !/    27-Jun-2006 : Adding file name preamble.          ( version 3.09 )
    !/    25-Jul-2006 : Adding grid ID per point.           ( version 3.10 )
    !/    27-Aug-2015 : Adding interpolation for the ice.   ( version 5.10 )
    !/    19-Jul-2021 : Momentum and air density support    ( version 7.14 )
    !/
    !  1. Purpose :
    !
    !     Read/write point output.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       INXOUT  C*(*)  I   Test string for read/write, valid are:
    !                          'READ' and 'WRITE'.
    !       NDSOP   Int.   I   File unit number. for binary
    !       NDSOA   Int.   I   File unit number. for ASCII
    !       IOTST   Int.   O   Test indictor for reading.
    !                           0 : Data read.
    !                          -1 : Past end of file.
    !       IMOD    I(O)   I   Model number for W3GDAT etc.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !      WW3_OUTP  Prog.   N/A    Postprocessing for point output.
    !      GX_OUTP   Prog.   N/A    Grads postprocessing for point output.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       Tests on INXOUT, file status and on array dimensions.
    !
    !  7. Remarks :
    !
    !     - The output file has the pre-defined name 'out_pnt.FILEXT'.
    !     - In MPP version of model data is supposed to be gatherd at the
    !       correct processor before the routine is called.
    !     - No error output filtering needed.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/SHRD  Switch for shared / distributed memory architecture.
    !     !/DIST  Id.
    !
    !     !/S     Enable subroutine tracing.
    !     !/T     Test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: w3setg
    USE w3wdatmd, ONLY: w3setw
    USE w3odatmd, ONLY: w3seto, w3dmo2
    !/
    USE w3gdatmd, ONLY: nth, nk, nspec, filext
    USE w3wdatmd, ONLY: time
    USE w3odatmd, ONLY: ndst, ndse, ipass => ipass2, nopts, iptint, &
         il, iw, ii, ptloc, ptifac, dpo, wao, wdo,   &
         aso, cao, cdo, spco, ptnme, o2init, fnmpre, &
         grdid, iceo, iceho, icefo
#ifdef W3_FLX5
    USE w3odatmd, ONLY: tauao, taudo, dairo
#endif
    USE w3odatmd, ONLY :  ofiles
    !/
#ifdef W3_SETUP
    USE w3odatmd, ONLY: zet_seto
#endif
    !/
    USE w3servmd, ONLY: extcde
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
 
    use constants, only: file_endian
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)           :: NDSOP
#ifdef W3_ASCII
    INTEGER, INTENT(IN), OPTIONAL :: NDSOA
#endif
    INTEGER, INTENT(OUT)          :: IOTST
    INTEGER, INTENT(IN), OPTIONAL :: IMOD
    CHARACTER, INTENT(IN)         :: INXOUT*(*)
    !/
    !/ ------------------------------------------------------------------- /
    !/ local parameters
    !/
    INTEGER                 :: IGRD, IERR, MK, MTH, I, J
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    LOGICAL,SAVE            :: WRITE
    CHARACTER(LEN=31)       :: IDTST
    CHARACTER(LEN=10)       :: VERTST
    !/
    CHARACTER(LEN=15) :: TIMETAG
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3IOPO')
#endif
    ipass  = ipass + 1
    iotst  = 0
    !
    ! test input parameters ---------------------------------------------- *
    !
    IF ( PRESENT(imod) ) THEN
      igrd   = imod
    ELSE
      igrd   = 1
    END IF
    !
    CALL w3seto ( igrd, ndse, ndst )
    CALL w3setg ( igrd, ndse, ndst )
    CALL w3setw ( igrd, ndse, ndst )
    !
    IF (inxout.NE.'READ' .AND. inxout.NE.'WRITE' ) THEN
      WRITE (ndse,900) inxout
      CALL extcde ( 1 )
    END IF
    !
    ! First pass to this file and we are only writing 1 file for all time     
    IF ( ipass.EQ.1  .AND. ofiles(2) .EQ. 0) THEN
      WRITE  = inxout.EQ.'WRITE'
    ELSE
      IF ( WRITE .AND. inxout.EQ.'READ' ) THEN
        WRITE (ndse,901) inxout
        CALL extcde ( 2 )
      END IF
    END IF
    !
    ! open file ---------------------------------------------------------- *
    !
    IF ( ipass.EQ.1 .AND. ofiles(2) .EQ. 0 ) THEN
 
      i      = len_trim(filext)
      j      = len_trim(fnmpre)
 
#ifdef W3_T
      WRITE (ndst,9001) fnmpre(:j)//'out_pnt.'//filext(:i)
#endif
      IF ( WRITE ) THEN
        OPEN (ndsop,file=fnmpre(:j)//'out_pnt.'//filext(:i),    &
             form='UNFORMATTED', convert=file_endian,err=800,iostat=ierr)
#ifdef W3_ASCII
        OPEN (ndsoa,file=fnmpre(:j)//'out_pnt.'//filext(:i)//'.txt',    &
             form='FORMATTED', err=800,iostat=ierr)
#endif
      ELSE
        OPEN (ndsop,file=fnmpre(:j)//'out_pnt.'//filext(:i),    &
             form='UNFORMATTED', convert=file_endian,err=800,iostat=ierr,status='OLD')
      END IF
      !
      rewind( ndsop )
      !
      ! test info ---------------------------------------------------------- *
      ! ( IPASS = 1 )
      !
      IF ( WRITE ) THEN
        WRITE (ndsop)                                           &
             idstr, veropt, nk, nth, nopts
#ifdef W3_ASCII
        WRITE (ndsoa,*)                                           &
             'IDSTR, VEROPT, NK, NTH, NOPTS:',                  &
             idstr, veropt, nk, nth, nopts
#endif
      ELSE
        READ (ndsop,END=801,ERR=802,IOSTAT=IERR)                &
             idtst, vertst, mk, mth, nopts
        !
        IF ( idtst .NE. idstr ) THEN
          WRITE (ndse,902) idtst, idstr
          CALL extcde ( 10 )
        END IF
        IF ( vertst .NE. veropt ) THEN
          WRITE (ndse,903) vertst, veropt
          CALL extcde ( 11 )
        END IF
        IF (nk.NE.mk .OR. nth.NE.mth) THEN
          WRITE (ndse,904) mk, mth, nk, nth
          CALL extcde ( 12 )
        END IF
        IF ( .NOT. o2init )                                     &
             CALL w3dmo2 ( igrd, ndse, ndst, nopts )
      END IF
      !
#ifdef W3_T
      WRITE (ndst,9002) idstr, veropt, nk, nth, nopts
#endif
      !
      ! Point specific info ------------------------------------------------ *
      ! ( IPASS = 1 )
      !
      IF ( WRITE ) THEN
        WRITE (ndsop)                                           &
             ((ptloc(j,i),j=1,2),i=1,nopts), (ptnme(i),i=1,nopts)
#ifdef W3_ASCII
        WRITE (ndsoa,*)                                           &
             '((PTLOC(J,I),J=1,2),I=1,NOPTS), (PTNME(I),I=1,NOPTS):', &
             ((ptloc(j,i),j=1,2),i=1,nopts), (ptnme(i),i=1,nopts)
#endif
      ELSE
        READ  (ndsop,END=801,ERR=802,IOSTAT=IERR)               &
             ((ptloc(j,i),j=1,2),i=1,nopts), (ptnme(i),i=1,nopts)
      END IF
      !
#ifdef W3_T
      WRITE (ndst,9003)
      DO i=1, nopts
        WRITE (ndst,9004) i, ptloc(1,i), ptloc(2,i), ptnme(i)
      END DO
#endif
      !
    END IF
    !
    !
    IF ( ipass.GE. 1  .AND. ofiles(2) .EQ. 1) THEN
      WRITE  = inxout.EQ.'WRITE'
    ELSE
      IF ( WRITE .AND. inxout.EQ.'READ' ) THEN
        WRITE (ndse,901) inxout
        CALL extcde ( 2 )
      END IF
    END IF
 
    ! open file ---------------------------------------------------------- *
    !
    IF ( ipass.GE.1 .AND. ofiles(2) .EQ. 1) THEN
      !
      i      = len_trim(filext)
      j      = len_trim(fnmpre)
 
      ! Create TIMETAG for file name using YYYYMMDD.HHMMS prefix
      WRITE(timetag,"(i8.8,'.'i6.6)")time(1),time(2)
      !
#ifdef W3_T
      WRITE (ndst,9001) fnmpre(:j)//timetag//'.out_pnt.'// &
           filext(:i)
#endif
      IF ( WRITE ) THEN
        OPEN (ndsop,file=fnmpre(:j)//timetag//'.out_pnt.'   &
             //filext(:i),form='UNFORMATTED', convert=file_endian,err=800,iostat=ierr)
#ifdef W3_ASCII
        OPEN (ndsoa,file=fnmpre(:j)//timetag//'.out_pnt.'   &
             //filext(:i)//'.txt',form='FORMATTED', err=800,iostat=ierr)
#endif
      END IF
      !
      rewind( ndsop )
      !
      !
      ! test info ---------------------------------------------------------- *
      ! ( IPASS GE.1 .AND. OFILES(2) .EQ. 1)
      !
      IF ( WRITE ) THEN
        WRITE (ndsop)                                           &
            idstr, veropt, nk, nth, nopts
#ifdef W3_ASCII
        WRITE (ndsoa,*)                                           &
             'IDSTR, VEROPT, NK, NTH, NOPTS:',                  &
             idstr, veropt, nk, nth, nopts
#endif
      ELSE
        READ (ndsop,END=801,ERR=802,IOSTAT=IERR)                &
             idtst, vertst, mk, mth, nopts
        !
        IF ( idtst .NE. idstr ) THEN
          WRITE (ndse,902) idtst, idstr
          CALL extcde ( 10 )
        END IF
        IF ( vertst .NE. veropt ) THEN
          WRITE (ndse,903) vertst, veropt
          CALL extcde ( 11 )
        END IF
        IF (nk.NE.mk .OR. nth.NE.mth) THEN
          WRITE (ndse,904) mk, mth, nk, nth
          CALL extcde ( 12 )
        END IF
        IF ( .NOT. o2init )                                     &
             CALL w3dmo2 ( igrd, ndse, ndst, nopts )
      END IF
      !
#ifdef W3_T
      WRITE (ndst,9002) idstr, veropt, nk, nth, nopts
#endif
      !
      ! Point specific info ------------------------------------------------ *
      ! ( IPASS GE.1 .AND. OFILES(2) .EQ. 1)
      !
      IF ( WRITE ) THEN
        WRITE (ndsop)                                           &
             ((ptloc(j,i),j=1,2),i=1,nopts), (ptnme(i),i=1,nopts)
#ifdef W3_ASCII
        WRITE (ndsoa,*)                                           &
             '((PTLOC(J,I),J=1,2),I=1,NOPTS), (PTNME(I),I=1,NOPTS):', &
             ((ptloc(j,i),j=1,2),i=1,nopts), (ptnme(i),i=1,nopts)
#endif
      ELSE
        READ  (ndsop,END=801,ERR=802,IOSTAT=IERR)               &
             ((ptloc(j,i),j=1,2),i=1,nopts), (ptnme(i),i=1,nopts)
      END IF
      !
#ifdef W3_T
      WRITE (ndst,9003)
      DO i=1, nopts
        WRITE (ndst,9004) i, ptloc(1,i), ptloc(2,i), ptnme(i)
      END DO
#endif
      !
    END IF
    !
    !
    ! TIME --------------------------------------------------------------- *
    !
    IF ( WRITE ) THEN
      WRITE (ndsop)                            time
#ifdef W3_ASCII
      WRITE (ndsoa,*) 'TIME:',                   time
#endif
    ELSE
      READ (ndsop,END=803,ERR=802,IOSTAT=IERR) time
    END IF
    !
#ifdef W3_T
    WRITE (ndst,9010) time
#endif
    !
    !
    ! Loop over spectra -------------------------------------------------- *
    !
    DO i=1, nopts
      !
      IF ( WRITE ) THEN
        ! set IW, II and IL to 0 because it is not used and gives &
        ! outlier values in out_pnt.points
        iw(i) = 0
        ii(i) = 0
        il(i) = 0
        WRITE (ndsop)                                            &
             iw(i), ii(i), il(i), dpo(i), wao(i), wdo(i),      &
#ifdef W3_FLX5
             tauao(i), taudo(i), dairo(i),                     &
#endif
#ifdef W3_SETUP
             zet_seto(i),                                      &
#endif
             aso(i), cao(i), cdo(i), iceo(i), iceho(i),        &
             icefo(i), grdid(i), (spco(j,i),j=1,nspec)
#ifdef W3_ASCII
        WRITE (ndsoa,*)                                            &
             'IW(I), II(I), IL(I), DPO(I), WAO(I), WDO(I):',      &
             iw(i), ii(i), il(i), dpo(i), wao(i), wdo(i),      &
#ifdef W3_FLX5
             'TAUAO(I), TAUDO(I), DAIRO(I):',                     &
             tauao(i), taudo(i), dairo(i),                     &
#endif
#ifdef W3_SETUP
             'ZET_SETO(I):',                                      &
             zet_seto(i),                                      &
#endif
             'ASO(I), CAO(I), CDO(I), ICEO(I), ICEHO(I):',        &
             aso(i), cao(i), cdo(i), iceo(i), iceho(i),        &
             'ICEFO(I), GRDID(I), (SPCO(J,I),J=1,NSPEC):',       &
             icefo(i), grdid(i), (spco(j,i),j=1,nspec)
#endif
      ELSE
        READ (ndsop,END=801,ERR=802,IOSTAT=IERR)                 &
             iw(i), ii(i), il(i), dpo(i), wao(i), wdo(i),      &
#ifdef W3_FLX5
             tauao(i), taudo(i), dairo(i),                     &
#endif
#ifdef W3_SETUP
             zet_seto(i),                                      &
#endif
             aso(i), cao(i), cdo(i), iceo(i), iceho(i),        &
             icefo(i), grdid(i), (spco(j,i),j=1,nspec)
      END IF
      !
    END DO
    IF (ofiles(2) .EQ. 1)  CLOSE (ndsop)
    !
    RETURN
    !
    ! Escape locations read errors
    !
800 CONTINUE
    WRITE (ndse,1000) ierr
    CALL extcde ( 20 )
    !
801 CONTINUE
    WRITE (ndse,1001)
    CALL extcde ( 21 )
    !
802 CONTINUE
    WRITE (ndse,1002) ierr
    CALL extcde ( 22 )
    !
803 CONTINUE
    iotst  = -1
#ifdef W3_T
    WRITE (ndst,9011)
#endif
    RETURN
    !
    ! Formats
    !
900 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOPO :'/                &
         '     ILEGAL INXOUT VALUE: ',a/)
901 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOPO :'/                &
         '     MIXED READ/WRITE, LAST REQUEST: ',a/)
902 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOPO :'/                &
         '     ILEGAL IDSTR, READ : ',a/                        &
         '                  CHECK : ',a/)
903 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOPO :'/                &
         '     ILEGAL VEROPT, READ : ',a/                       &
         '                   CHECK : ',a/)
904 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOPO :'/                &
         '     ERROR IN SPECTRA, MK, MTH : ',2i8/               &
         '              ARRAY DIMENSIONS : ',2i8/)
    !
1000 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOPO : '/               &
         '     ERROR IN OPENING FILE'/                          &
         '     IOSTAT =',i5/)
1001 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOPO : '/               &
         '     PREMATURE END OF FILE'/)
1002 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOPO : '/               &
         '     ERROR IN READING FROM FILE'/                     &
         '     IOSTAT =',i5/)
    !
#ifdef W3_T
9000 FORMAT (' TEST W3IOPO : IPASS =',i4,'    INXOUT = ',a,       &
         ' WRITE = ',l1,' UNIT =',i3/                         &
         '               IGRD =',i3,' FEXT = ',a)
 
9001 FORMAT (' TEST W3IOPO : OPENING NEW FILE [',a,']')
9002 FORMAT (' TEST W3IOPO : TEST PARAMETERS:'/                   &
         '       IDSTR : ',a/                                 &
         '      VEROPT : ',a/                                 &
         '      NK,NTH :',i5,i8/                              &
         '        NOPT :',i5)
9003 FORMAT (' TEST W3IOPO : POINT LOCATION AND ID')
9004 FORMAT (3x,i4,2f10.2,2x,a)
    !
9010 FORMAT (' TEST W3IOPO : TIME  :',i9.8,i7.6)
9011 FORMAT (' TEST W3IOPO : END OF FILE REACHED')
    !
9020 FORMAT (' TEST W3IOPO : POINT NR.:',i5)
9021 FORMAT (' TEST W3IOPO :',2i4,2f6.3)
9022 FORMAT (' TEST W3IOPO :',4i7,2x,4i2,2x,4f5.2)
9030 FORMAT (' TEST W3IOPO :',f8.1,2(f7.2,f7.1))
#endif
    !/
    !/ End of W3IOPO ----------------------------------------------------- /
    !/
  END SUBROUTINE w3iopo
  !/
  !/ End of module W3IOPOMD -------------------------------------------- /
  !/
END MODULE w3iopomd