WW3/w3iogomd_8F90_source.html

#include "w3macros.h"


!/ ------------------------------------------------------------------- /

MODULE w3iogomd

  !/

  !/                  +-----------------------------------+

  !/                  | WAVEWATCH III           NOAA/NCEP |

  !/                  |           H. L. Tolman            |

  !/                  |                        FORTRAN 90 |

  !/                  | Last update :         02-Mar-2024 |

  !/                  +-----------------------------------+

  !/

  !/    04-Jan-2001 : Origination.                        ( version 2.00 )

  !/    23-Apr-2002 : Clean up.                           ( version 2.19 )

  !/    29-Apr-2002 : Add output parameters 17-18.        ( version 2.20 )

  !/    30-May-2002 : Switch clean up.                    ( version 2.21 )

  !/    13-Nov-2002 : Add stress vector.                  ( version 3.00 )

  !/    25-Oct-2004 : Multiple grid version.              ( version 3.06 )

  !/    27-Jun-2005 : Adding MAPST2.                      ( version 3.07 )

  !/    21-Jul-2005 : Adding output fields 19-21.         ( version 3.07 )

  !/    23-Apr-2006 : Filter for directional spread.      ( version 3.09 )

  !/    27-Jun-2006 : Adding file name preamble.          ( version 3.09 )

  !/    05-Jul-2006 : Consolidate stress arrays.          ( version 3.09 )

  !/    02-Apr-2007 : Adding partitioned output.          ( version 3.11 )

  !/                  Adding user slots for outputs.

  !/    08-Oct-2007 : Adding ST3 source term option.      ( version 3.13 )

  !/                  ( F. Ardhuin )

  !/    05-Mar-2008 : Added NEC sxf90 compiler directives

  !/                  (Chris Bunney, UK Met Office)       ( version 3.13 )

  !/    29-May-2009 : Preparing distribution version.     ( version 3.14 )

  !/    13-Sep-2009 : Add coupling option                 ( version 3.14 )

  !/    10-Mar-2009 : Add second order pressure           ( version 3.14 )

  !/    15-Sep-2010 : Adding ST4 source term option.      ( version 3.14 )

  !/    30-Oct-2009 : Implement curvilinear grid type.    ( version 3.14 )

  !/                  (W. E. Rogers & T. J. Campbell, NRL)

  !/    05-Feb-2011 : Implement unstructured grid        ( version 3.14.3 )

  !/                  (A. Roland and F. Ardhuin)

  !/    12-Jun-2012 : Add /RTD option or rotated grid option.

  !/                  (Jian-Guo Li)                       ( version 4.06 )

  !/    25-Dec-2012 : New output structure and smaller    ( version 4.11 )

  !/                  memory footprint.

  !/    15-Apr-2013 : New subroutine to read param. names ( version 4.11 )

  !/    21-Aug-2013 : Bug correction in W3IOGO: UBR, ABR  ( version 4.11 )

  !/    11-Nov-2013 : SMC and rotated grid incorporated in the main

  !/                  trunk                               ( version 4.13 )

  !/    31-Jan-2014 : Bug fix warning output (Tolman).    ( version 4.18 )

  !/    10-Feb-2014 : Bug correction for US3D: div. by df ( version 4.18 )

  !/    30-Apr-2014 : Add th2m and sth2m calculation      ( version 5.01 )

  !/    27-May-2014 : Switch to OMPG switch.              ( version 5.02 )

  !/    27-Aug-2015 : Add ICEF,ICEH as output fields      ( version 5.10 )

  !/    01-Mar-2018 : Removed RTD code (now used in post  ( version 6.02 )

  !/                  processing code)

  !/    05-Jun-2018 : Add DEBUGSTP/SETUP                  ( version 6.04 )

  !/    22-Aug-2018 : Add WBT output parameter            ( version 6.06 )

  !/    25-Sep-2019 : Corrected th2m and sth2m            ( version 6.07 )

  !/                  calculations. (J Dykes, NRL)

  !/    04-Oct-2019 : Optional one file per output stride ( version 7.00 )

  !/                  (Roberto Padilla-Hernandez & J.H. Alves)

  !/    03-Nov-2020 : Factored out NAME matching into     ( version 7.12 )

  !/                  seperate subroutine. (C. Bunney)

  !/    15-Jan-2021 : Added TP output based on exsiting   ( version 7.12 )

  !/                  FP internal field. (C. Bunney)

  !/    22-Mar-2021 : Add extra coupling fields as output ( version 7.13 )

  !/    21-Jul-2022 : Correct FP0 calc for peak energy in ( version 7.14 )

  !/                  min/max freq band (B. Pouliot, CMC)

  !/    02-Mar-2024 : Add skweness and EM bias varaible   ( version 7.xx )

  !/

  !/    Copyright 2009-2024 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 :

  !

  !     Gridded output of mean wave parameters.

  !

  !  2. Variables and types :

  !

  !      Name      Type  Scope    Description

  !     ----------------------------------------------------------------

  !      VEROGR    C*10  Private  Gridded output file version number.

  !      IDSTR     C*30  Private  Gridded output file ID string.

  !     ----------------------------------------------------------------

  !

  !  3. Subroutines and functions :

  !

  !      Name      Type  Scope    Description

  !     ----------------------------------------------------------------

  !      W3OUTG    Subr. Public   Calculate mean parameters.

  !      W3IOGO    Subr. Public   IO to raw gridded fields file.

  !     ----------------------------------------------------------------

  !

  !  4. Subroutines and functions used :

  !

  !      Name      Type  Module   Description

  !     ----------------------------------------------------------------

  !      W3SETO    Subr. W3ODATMD Point to data structure.

  !      W3SETG    Subr. W3GDATMD Point to data structure.

  !      W3SETW    Subr. W3WDATMD Point to data structure.

  !      W3SETA    Subr. W3ADATMD Point to data structure.

  !      W3XETA    Subr. W3ADATMD Point to data structure.

  !      W3DIMW    Subr. W3WDATMD Allocate data structure.

  !      W3DIMA    Subr. W3ADATMD Allocate data structure.

  !      STRACE    Subr. W3SERVMD Subroutine tracing.           ( !/S )

  !      EXTCDE    Subr. W3SERVMD Program abort with exit code.

  !     ----------------------------------------------------------------

  !

  !  5. Remarks :

  !

  !     - The different output fields are not folded in with this module

  !       due to the different requirements for a element '0' in some of

  !       the fields.

  !

  !  6. Switches :

  !

  !       !/SHRD  Switch for shared / distributed memory architecture.

  !       !/DIST  Id.

  !

  !       !/OMPG  OpenMP compiler directive for loop splitting.

  !

  !       !/O8    Filter for low wave heights ( HSMIN )

  !       !/O9    Negative wave height alowed, other mean parameters will

  !             not be correct.

  !

  !       !/ST0   No source terms.

  !       !/ST1   Source term set 1 (WAM equiv.)

  !       !/ST2   Source term set 2 (Tolman and Chalikov)

  !       !/ST3   Source term set 3 (WAM 4+)

  !       !/ST4   Source term set 4 (Ardhuin et al. 2009, 2010)

  !       !/ST6   Source term set 6 (BYDRZ)

  !

  !       !/S     Enable subroutine tracing.

  !       !/T     Test output.

  !

  !  7. Source code :

  !

  !/ ------------------------------------------------------------------- /

#ifdef W3_S

  USE w3servmd, ONLY : strace

#endif

  !/

  PUBLIC

  CHARACTER(LEN=1024)                   :: fldout

  !/

  !/ Private parameter statements (ID strings)

  !/

  CHARACTER(LEN=10), PARAMETER, PRIVATE :: verogr = '2019-10-04'

  CHARACTER(LEN=30), PARAMETER, PRIVATE ::                        &

       idstr = 'WAVEWATCH III GRID OUTPUT FILE'

  !/

CONTAINS

  !/ ------------------------------------------------------------------- /

  SUBROUTINE w3flgrdupdt ( NDSO, NDSEN, FLGRD, FLGR2, FLGD, FLG2 )

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           F. Ardhuin              |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         15-Apr-2013 |

    !/                  +-----------------------------------+

    !/

    !/    15-Apr-2013 : Origination.                        ( version 4.10 )

    !/

    !  1. Purpose :

    !

    !     Updates the flags for output parameters based on the mod_def file

    !     this is to prevent the allocation of big 3D arrays when not requested

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !       NDSO      Int.   I   Output file logical unit number

    !       NDSEN     R.A.   I   Error output file logical unit number

    !       FLGD,FLG2 L.A.   O   1D array of flags for groups

    !       FLGRD     L.A.   O   2D array of flags

    !       FLGR2     L.A.   O   2D array of flags

    !     ----------------------------------------------------------------

    !

    !

    !  4. Subroutines used :

    !

    !     None

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      W3INIT    Subr.   N/A

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !     None.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S     Enable subroutine tracing.

    !     !/T     Test output.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

    USE constants

    USE w3gdatmd, ONLY: e3df, p2msf, us3df, usspf

    USE w3odatmd, ONLY: nogrp, ngrpp

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    INTEGER, INTENT(IN)     :: NDSO, NDSEN

    LOGICAL, INTENT(INOUT)  :: FLGRD(NOGRP,NGRPP), FLGD(NOGRP),     &

         FLGR2(NOGRP,NGRPP), FLG2(NOGRP)

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

    INTEGER             :: I

    CHARACTER(LEN=10)  :: VARNAME1(5),VARNAME2(5)

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    !/

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'W3FLGRDUPDT')

#endif

    !

    varname1(1) = 'EF';    varname2(1) = 'E3D'

    varname1(2) = 'TH1M';  varname2(2) = 'TH1MF'

    varname1(3) = 'STH1M'; varname2(3) = 'STH1MF'

    varname1(4) = 'TH2M';  varname2(4) = 'TH2MF'

    varname1(5) = 'STH2M'; varname2(5) = 'STH2MF'


    DO i=1,5

      IF (e3df(1,i).LE.0.OR.e3df(3,i).LT.e3df(2,i)) THEN

        IF (flgrd(3,i).OR.flgr2(3,i)) THEN

          WRITE(ndsen,1008) varname1(i),varname2(i)

        END IF

        flgrd(3,i)=.false.

        flgr2(3,i)=.false.

      END IF

    END DO

    IF (us3df(1).LE.0.OR.us3df(3).LT.us3df(2)) THEN

      IF (flgrd(6,8).OR.flgr2(6,8)) THEN

        WRITE(ndsen,1008) 'USF','US3D'

      END IF

      flgrd(6,8)=.false.

      flgr2(6,8)=.false.

    END IF

    IF (usspf(1).LE.0.OR.usspf(2).LE.0) THEN

      IF (flgrd(6,12).OR.flgr2(6,12)) THEN

        WRITE(ndsen,1008) 'USP','USSP'

      END IF

      flgrd(6,12)=.false.

      flgr2(6,12)=.false.

    END IF

    IF (p2msf(1).LE.0.OR.p2msf(3).LT.p2msf(2)) THEN

      IF (flgrd(6,9).OR.flgr2(6,9)) THEN

        WRITE(ndsen,1008) 'P2L','P2SF'

      END IF

      flgrd(6,9)=.false.

      flgr2(6,9)=.false.

    END IF

    !

    flgd(3) = .false.

    flg2(3) = .false.

    IF(any(flgrd(3,:))) flgd(3)=.true.

    IF(any(flgr2(3,:))) flg2(3)=.true.

    flgd(6) = .false.

    flg2(6) = .false.

    IF(any(flgrd(6,:))) flgd(6)=.true.

    IF(any(flgr2(6,:))) flg2(6)=.true.

    !

    RETURN

    !

1008 FORMAT (/' *** WAVEWATCH III WARNING  : '/                       &

         '     PARAMETER ',a,' not allowed: need to set',        &

         ' parameter ',a,' in OUTS namelist (in ww3_grid.inp)'   &

         ' with proper bounds' )

    !

  END SUBROUTINE w3flgrdupdt

  !/ ------------------------------------------------------------------- /

  SUBROUTINE w3readflgrd ( NDSI , NDSO, NDSS, NDSEN, COMSTR,      &

       FLG1D, FLG2D, IAPROC, NAPOUT, IERR)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           F. Ardhuin              |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         25-Sep-2020 |

    !/                  +-----------------------------------+

    !/

    !/    15-Apr-2013 : Origination.                        ( version 4.10 )

    !/    31-Jan-2014 : Bug fix warning output (Tolman).    ( version 4.18 )

    !/    30-Apr-2014 : Add th2m and sth2m calculation      ( version 5.01 )

    !/    25-Sep-2020 : Calculate FLG1D for any processor   ( version 7.10 )

    !/    03-Nov-2020 : Factored out NAME matching into     ( version 7.12 )

    !/                  seperate subroutine (C. Bunney)

    !/

    !  1. Purpose :

    !

    !     Fills in FLG1D and FLG2D arrays from ASCII input file

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !       NDSI    Int.   I   Input file logical unit number

    !       NDSO    Int.   I   Output file logical unit number

    !       NDSS    Int.   I   Screen file logical unit number

    !       NDSEN   R.A.   I   Error output file logical unit number

    !       COMSTR  Char   I   Comment string, usually '$'

    !       FLG1D   L.A.   O   1D array of flags for groups

    !       FLG2D   L.A.   O   2D array of flags

    !       IAPROC  Int.   I   index of current processor

    !       NAPOUT  Int.   I   index of processor for output (screen)

    !       IERR    Int.   O   Error message number

    !     ----------------------------------------------------------------

    !

    !

    !  4. Subroutines used :

    !

    !     None

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      WW3_SHEL  Prog.   N/A    Actual wave model program

    !      WW3_OUTF  Prog.   N/A    Output postprocessor.

    !      WW3_OUNF  Prog.   N/A    NetCDF output postprocessor.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !     None.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S     Enable subroutine tracing.

    !     !/T     Test output.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

    USE constants

    USE w3gdatmd, ONLY: us3df, usspf

    USE w3odatmd, ONLY: nogrp, ngrpp, noge, idout

    USE w3servmd, ONLY: nextln, strsplit, str_to_upper

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    INTEGER, INTENT(IN)     :: NDSI, NDSO, NDSS, NDSEN, IAPROC, NAPOUT

    INTEGER, INTENT(OUT)    :: IERR

    CHARACTER(LEN=1)        :: COMSTR

    LOGICAL, INTENT(OUT)    :: FLG2D(NOGRP,NGRPP), FLG1D(NOGRP)

    CHARACTER(LEN=100)      :: OUT_NAMES(100), TESTSTR

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

    INTEGER             :: IFI, IFJ, IOUT

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    CHARACTER(LEN=1)    :: AFLG

    LOGICAL             :: FLT, NAMES

    !/

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'W3READFLGRD')

#endif

    !

    !

    ! 1.  Initialize flags -------------------------------------- *

    !

    ierr=0

    flg2d(:,:)=.false. ! Initialize FLG2D

    flg1d(:)=.false. ! Initialize FLOG

    names =.false.

    !

    DO ifi=1,nogrp ! Loop over field output groups

      !

      CALL nextln ( comstr , ndsi , ndsen )

      READ (ndsi,*,END=2001,ERR=2002) aflg

      IF (aflg.EQ.'T') THEN

        flg1d(ifi)=.true.

      ELSE IF (aflg.EQ.'F') THEN

        flg1d(ifi)=.false.

      ELSE IF (aflg.EQ.'N') THEN

        names=.true.

        EXIT

      ELSE

        ierr=1

        GOTO 2005

      END IF

      IF ( flg1d(ifi) ) THEN ! Skip if group not requested

        CALL nextln ( comstr , ndsi , ndsen )

        READ (ndsi,'(A)',END=2001,ERR=2006,IOSTAT=IERR)              &

             fldout

        out_names(:)=''

        CALL strsplit(fldout,out_names)

        ifj=0

        DO WHILE (len_trim(out_names(ifj+1)).NE.0)

          ifj=ifj+1

          IF ( out_names(ifj) .EQ. 'T' )                            &

               flg2d(ifi,ifj)=.true.

        ENDDO

        IF ( iaproc .EQ. napout .AND. ifj .LT. noge(ifi) ) WRITE(ndsen,1007) ifi

      ENDIF

    END DO

    !

    IF (names) THEN

      !

      ! 2. Reads and splits list of output field names

      !

      CALL nextln ( comstr , ndsi , ndsen )

      READ (ndsi,'(A)',END=2001,ERR=2003,IOSTAT=IERR) fldout

      out_names(:)=''

      CALL strsplit(fldout,out_names)

      iout=0

      DO WHILE (len_trim(out_names(iout+1)).NE.0)

        CALL str_to_upper(out_names(iout+1))

        !

        ! 2. Matches names with expected ...

        !

        teststr=out_names(iout+1)

        CALL w3fldtoij(teststr, ifi, ifj, iaproc, napout, ndsen)


        IF(ifi .NE. -1) THEN

          flg2d(ifi, ifj) = .true.

        ENDIF

        !

        iout=iout+1

        !

      END DO

      !

    END IF

    !

    flt    = .true.

    DO ifi=1, nogrp

      IF ( iaproc .EQ. napout ) THEN

        DO ifj=1, ngrpp

          IF ( flg2d(ifi,ifj) ) THEN

            IF ( flt ) THEN

              WRITE (ndso,1945) idout(ifi,ifj)

              flt    = .false.

            ELSE

              WRITE (ndso,1946) idout(ifi,ifj)

            END IF

          END IF

        END DO

      END IF

      IF(any(flg2d(ifi,:))) flg1d(ifi)=.true. !Update FLG1D

    END DO

    IF ( iaproc .EQ. napout ) THEN

      IF ( flt ) WRITE (ndso,1945) 'no fields defined'

    END IF

    !

    RETURN

    !

2001 CONTINUE

    IF ( iaproc .EQ. napout ) WRITE (ndsen,1001)

    RETURN

2002 CONTINUE

    IF ( iaproc .EQ. napout ) WRITE (ndsen, 1002) ifi, ierr

    RETURN

2003 CONTINUE

    IF ( iaproc .EQ. napout ) WRITE (ndsen, 1003) ierr

    RETURN

    !2004 CONTINUE ! replaced by warning in code ....

2005 CONTINUE

    IF ( iaproc .EQ. napout ) WRITE (ndsen, 1005) aflg

    RETURN

2006 CONTINUE

    IF ( iaproc .EQ. napout ) WRITE (ndsen, 1006) ifi,ierr

    RETURN

    !

1945 FORMAT ( '            Fields   : ',a)

1946 FORMAT ( '                       ',a)

    !

1001 FORMAT (/' *** WAVEWATCH III ERROR  : '/                         &

         '     PREMATURE END OF INPUT FILE'/)

    !

1002 FORMAT (/' *** WAVEWATCH III ERROR  : '/                         &

         '     ERROR IN READING OUTPUT FIELDS GROUP FLAGS ',     &

         i2, /, '     IOSTAT =',i5/)

    !

1003 FORMAT (/' *** WAVEWATCH III ERROR  : '/                         &

         '     ERROR READING OUTPUT FIELD NAMES FROM INPUT FILE'/&

         '     IOSTAT =',i5/)

    !

1005 FORMAT (/' *** WAVEWATCH III ERROR    : '/                       &

         '     WAS EXPECTING "T" "F" or "N", but found "',a,'".'/)

    !

1006 FORMAT (/' *** WAVEWATCH III ERROR  : '/                         &

         '     ERROR IN READING OUTPUT FIELDS FLAGS FOR GROUP ', &

         i2, /, '     IOSTAT =',i5/)

    !

1007 FORMAT (/' *** WAVEWATCH III WARNING  : '/                       &

         '     NUMBER OF REQUESTED OUTPUT FIELD FLAGS IN GROUP ',&

         i2, /,' LESS THAN AVAILABLE, CHECK DOCS FOR MORE OPTIONS')

    !

  END SUBROUTINE w3readflgrd


  !/ ------------------------------------------------------------------- /

  SUBROUTINE w3flgrdflag ( NDSO, NDSS, NDSEN, FLDOUT,      &

       FLG1D, FLG2D, IAPROC, NAPOUT, IERR)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           F. Ardhuin              |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         25-Sep-2020 |

    !/                  +-----------------------------------+

    !/

    !/    15-Apr-2013 : Origination.                        ( version 4.10 )

    !/    31-Jan-2014 : Bug fix warning output (Tolman).    ( version 4.18 )

    !/    30-Apr-2014 : Add th2m and sth2m calculation      ( version 5.01 )

    !/    17-Feb-2016 : New version for namelist use        ( version 5.11 )

    !/    25-Sep-2020 : Calculate FLG1D for any processor   ( version 7.10 )

    !/    03-Nov-2020 : Factored out NAME matching into     ( version 7.12 )

    !/                  seperate subroutine (C. Bunney)

    !/

    !  1. Purpose :

    !

    !     Fills in FLG1D and FLG2D arrays from ASCII input file

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !       NDSO    Int.   I   Output file logical unit number

    !       NDSS    Int.   I   Screen file logical unit number

    !       NDSEN   R.A.   I   Error output file logical unit number

    !       FLDOUT  Cha.   I   List of field names

    !       FLG1D   L.A.   O   1D array of flags for groups

    !       FLG2D   L.A.   O   2D array of flags

    !       IAPROC  Int.   I   index of current processor

    !       NAPOUT  Int.   I   index of processor for output (screen)

    !       IERR    Int.   O   Error message number

    !     ----------------------------------------------------------------

    !

    !

    !  4. Subroutines used :

    !

    !     None

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      WW3_SHEL  Prog.   N/A    Actual wave model program

    !      WW3_OUTF  Prog.   N/A    Output postprocessor.

    !      WW3_OUNF  Prog.   N/A    NetCDF output postprocessor.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !     None.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S     Enable subroutine tracing.

    !     !/T     Test output.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

    USE constants

    USE w3odatmd, ONLY: nogrp, ngrpp, idout

    USE w3servmd, ONLY: strsplit, str_to_upper

    USE w3gdatmd, ONLY: us3df, usspf

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    INTEGER, INTENT(IN)     :: NDSO, NDSS, NDSEN, IAPROC, NAPOUT

    CHARACTER(1024), INTENT(IN)   :: FLDOUT

    INTEGER, INTENT(OUT)    :: IERR

    LOGICAL, INTENT(OUT)    :: FLG2D(NOGRP,NGRPP), FLG1D(NOGRP)

    CHARACTER(LEN=100)      :: OUT_NAMES(100), TESTSTR

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

    INTEGER             :: I, IFI, IFJ, IOUT

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    LOGICAL             :: FLT

    !/

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'W3FLGRDFLAG')

#endif

    !

    !

    ! 1.  Initialize flags -------------------------------------- *

    !

    ierr=0

    flg2d(:,:)=.false. ! Initialize FLG2D

    flg1d(:)=.false. ! Initialize FLOG

    !

    ! 2. Splits list of output field names

    !

    out_names(:)=''

    CALL strsplit(fldout,out_names)

    iout=0

    DO WHILE (len_trim(out_names(iout+1)).NE.0)

      CALL str_to_upper(out_names(iout+1))

      !

      ! 2. Matches names with expected ...

      !

      teststr=out_names(iout+1)

      CALL w3fldtoij(teststr, ifi, ifj, iaproc, napout, ndsen)


      IF(ifi .NE. -1) THEN

        flg2d(ifi, ifj) = .true.

      ENDIF

      !

      iout=iout+1

      !

    END DO

    !

    flt    = .true.

    DO ifi=1, nogrp

      IF ( iaproc .EQ. napout ) THEN

        DO ifj=1, ngrpp

          IF ( flg2d(ifi,ifj) ) THEN

            IF ( flt ) THEN

              WRITE (ndso,1945) idout(ifi,ifj)

              flt    = .false.

            ELSE

              WRITE (ndso,1946) idout(ifi,ifj)

            END IF

          END IF

        END DO

      ENDIF

      IF(any(flg2d(ifi,:))) flg1d(ifi)=.true. !Update FLG1D

    END DO

    IF ( iaproc .EQ. napout ) THEN

      IF ( flt ) WRITE (ndso,1945) 'no fields defined'

    ENDIF

    !

    RETURN

    !

1945 FORMAT ( '            Fields   : ',a)

1946 FORMAT ( '                       ',a)

    !

    ! 1004 FORMAT (/' *** WAVEWATCH III WARNING  : '/                       &

    !               '     REQUESTED OUTPUT FIELD ',A,' WAS NOT RECOGNIZED.'/)

    !!

    ! 1008 FORMAT (/' *** WAVEWATCH III WARNING  : '/                       &

    !               '     PARAMETER ',A,' not allowed: need to set',        &

    !               ' parameter ',A,' in OUTS namelist (in ww3_grid.inp)')

    !

  END SUBROUTINE w3flgrdflag


  !/ ------------------------------------------------------------------- /

  SUBROUTINE w3fldtoij(FLD, I, J, IAPROC, NAPOUT, NDSEN)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |            C. Bunney              |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         22-Mar-2021 |

    !/                  +-----------------------------------+

    !/

    !/    03-Nov-2020 : Origination.                        ( version 7.12 )

    !/    22-Mar-2021 : Add extra coupling fields as output ( version 7.13 )

    !

    !  1. Purpose :

    !

    !     Returns the group/field (I/J) indices for a named output field.

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !       FLD     Cha.   I   Field names

    !       I       Int.   O   Output group number (IFI)

    !       J       Int.   O   Output field number (IFJ)

    !       IAPROC  Int.   I   index of current processor

    !       NAPOUT  Int.   I   index of processor for output (screen)

    !       NDSEN   R.A.   I   Error output file logical unit number

    !     ----------------------------------------------------------------

    !

    !/ ------------------------------------------------------------------- /

    USE w3gdatmd, ONLY: us3df, usspf

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

    CHARACTER(LEN=*), INTENT(IN) :: FLD

    INTEGER, INTENT(IN) :: IAPROC, NAPOUT, NDSEN

    INTEGER, INTENT(OUT) :: I, J


    i = -1

    j = -1


    SELECT CASE(trim(fld(1:6)))

      !

      ! Group 1

      !

    CASE('DPT')

      i = 1

      j = 1

    CASE('CUR')

      i = 1

      j = 2

    CASE('WND')

      i = 1

      j = 3

    CASE('AST')

      i = 1

      j = 4

    CASE('WLV')

      i = 1

      j = 5

    CASE('ICE')

      i = 1

      j = 6

    CASE('IBG')

      i = 1

      j = 7

    CASE('TAU')

      i = 1

      j = 8

    CASE('RHO')

      i = 1

      j = 9

#ifdef W3_BT4

    CASE('D50')

      i = 1

      j = 10

#endif

#ifdef W3_IS2

    CASE('IC1')

      i = 1

      j = 11

    CASE('IC5')

      i = 1

      j = 12

#endif

      ! Group 2

      !

#ifdef W3_OASACM

    CASE('AHS')

      i = 2

      j = 1

#endif

#ifdef W3_OASOCM

    CASE('OHS')

      i = 2

      j = 1

#endif

    CASE('HS')

      i = 2

      j = 1

    CASE('LM')

      i = 2

      j = 2

    CASE('T02')

      i = 2

      j = 3

    CASE('T0M1')

      i = 2

      j = 4

    CASE('T01')

      i = 2

      j = 5

    CASE('FP')

      i = 2

      j = 6

    CASE('DIR')

      i = 2

      j = 7

    CASE('SPR')

      i = 2

      j = 8

    CASE('DP')

      i = 2

      j = 9

    CASE('HIG')

      i = 2

      j = 10

    CASE('MXE')

      i = 2

      j = 11

    CASE('MXES')

      i = 2

      j = 12

    CASE('MXH')

      i = 2

      j = 13

    CASE('MXHC')

      i = 2

      j = 14

    CASE('SDMH')

      i = 2

      j = 15

    CASE('SDMHC')

      i = 2

      j = 16

    CASE('WBT')

      i = 2

      j = 17

    CASE('TP') ! Uses FP0 internally, as per FP

      i = 2

      j = 18

    CASE('WNM')

      i = 2

      j = 19

#ifdef W3_OASOCM

    CASE('THM')

      i = 2

      j = 20

#endif

      !

      ! Group 3

      !

    CASE('EF')

      i = 3

      j = 1

    CASE('TH1M')

      i = 3

      j = 2

    CASE('STH1M')

      i = 3

      j = 3

    CASE('TH2M')

      i = 3

      j = 4

    CASE('STH2M')

      i = 3

      j = 5

    CASE('WN')

      i = 3

      j = 6

      !

      ! Group 4

      !

    CASE('PHS')

      i = 4

      j = 1

    CASE('PTP')

      i = 4

      j = 2

    CASE('PLP')

      i = 4

      j = 3

    CASE('PDIR')

      i = 4

      j = 4

    CASE('PSPR')

      i = 4

      j = 5

    CASE('PWS')

      i = 4

      j = 6

    CASE('PDP')

      i = 4

      j = 7

    CASE('PQP')

      i = 4

      j = 8

    CASE('PPE')

      i = 4

      j = 9

    CASE('PGW')

      i = 4

      j = 10

    CASE('PSW')

      i = 4

      j = 11

    CASE('PTM10')

      i = 4

      j = 12

    CASE('PT01')

      i = 4

      j = 13

    CASE('PT02')

      i = 4

      j = 14

    CASE('PEP')

      i = 4

      j = 15

    CASE('TWS')

      i = 4

      j = 16

    CASE('PNR')

      i = 4

      j = 17

      !

      ! Group 5

      !

    CASE('UST')

      i = 5

      j = 1

#ifdef W3_OASACM

    CASE('ACHA')

      i = 5

      j = 2

#endif

#ifdef W3_OASOCM

    CASE('OCHA')

      i = 5

      j = 2

#endif

    CASE('CHA')

      i = 5

      j = 2

    CASE('CGE')

      i = 5

      j = 3

    CASE('FAW')

      i = 5

      j = 4

    CASE('TAW')

      i = 5

      j = 5

    CASE('TWA')

      i = 5

      j = 6

    CASE('WCC')

      i = 5

      j = 7

    CASE('WCF')

      i = 5

      j = 8

    CASE('WCH')

      i = 5

      j = 9

    CASE('WCM')

      i = 5

      j = 10

    CASE('FWS')

      i = 5

      j = 11

      !

      ! Group 6

      !

    CASE('SXY')

      i = 6

      j = 1

    CASE('TWO')

      i = 6

      j = 2

    CASE('BHD')

      i = 6

      j = 3

    CASE('FOC')

      i = 6

      j = 4

    CASE('TUS')

      i = 6

      j = 5

    CASE('USS')

      i = 6

      j = 6

    CASE('P2S')

      i = 6

      j = 7

    CASE('USF')

      IF (us3df(1).GE.1) THEN

        i = 6

        j = 8

      ELSE

        IF ( iaproc .EQ. napout ) WRITE(ndsen,1008) 'USF','US3D'

      END IF

    CASE('P2L')

      i = 6

      j = 9

    CASE('TWI')

      i = 6

      j = 10

    CASE('FIC')

      i = 6

      j = 11

    CASE('USP')

      IF (usspf(1).GE.1) THEN

        i = 6

        j = 12

      ELSE

        IF ( iaproc .EQ. napout ) WRITE(ndsen,1008) 'USP','USSP'

      END IF

    CASE('TOC')

      i = 6

      j = 13

      !

      ! Group 7

      !

    CASE('ABR')

      i = 7

      j = 1

    CASE('UBR')

      i = 7

      j = 2

    CASE('BED')

      i = 7

      j = 3

    CASE('FBB')

      i = 7

      j = 4

    CASE('TBB')

      i = 7

      j = 5

      !

      ! Group 8

      !

    CASE('MSS')

      i = 8

      j = 1

    CASE('MSC')

      i = 8

      j = 2

    CASE('MSD')

      i = 8

      j = 3

    CASE('MCD')

      i = 8

      j = 4

    CASE('QP')

      i = 8

      j = 5

    CASE('QKK')

      i = 8

      j = 6

    CASE('SKW')

      i = 8

      j = 7

    CASE('EMB')

      i = 8

      j = 8

    CASE('EMC')

      i = 8

      j = 9

      !

      ! Group 9

      !

    CASE('DTD')

      i = 9

      j = 1

    CASE('FC')

      i = 9

      j = 2

    CASE('CFX')

      i = 9

      j = 3

    CASE('CFD')

      i = 9

      j = 4

    CASE('CFK')

      i = 9

      j = 5

      !

      ! Group 10

      !

    CASE('U1')

      i = 10

      j = 1

    CASE('U2')

      i = 10

      j = 1

      ! Not found:

#ifdef W3_COU

    CASE('DRY')

#endif

    CASE('UNSET')

    CASE DEFAULT

      i = -1

      j = -1

      IF ( iaproc .EQ. napout ) WRITE (ndsen,1004) trim(fld)

    END SELECT


1004 FORMAT (/' *** WAVEWATCH III WARNING  : '/                       &

         '     REQUESTED OUTPUT FIELD ',a,' WAS NOT RECOGNIZED.'/)

    !

1008 FORMAT (/' *** WAVEWATCH III WARNING  : '/                       &

         '     PARAMETER ',a,' not allowed: need to set',        &

         ' parameter ',a,' in OUTS namelist (in ww3_grid.inp)')

    !

  END SUBROUTINE w3fldtoij


  !/ ------------------------------------------------------------------- /

  SUBROUTINE w3outg ( A, FLPART, FLOUTG, FLOUTG2 )

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           H. L. Tolman            |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         10-Apr-2015 |

    !/                  +-----------------------------------+

    !/

    !/    10-Dec-1998 : Distributed FORTRAN 77 version.     ( version 1.18 )

    !/    04-Jan-2000 : Upgrade to FORTRAN 90               ( version 2.00 )

    !/                  Major changes to logistics.

    !/    09-May-2002 : Switch clean up.                    ( version 2.21 )

    !/    19-Oct-2004 : Multiple grid version.              ( version 3.06 )

    !/    21-Jul-2005 : Adding output fields 19-21.         ( version 3.07 )

    !/    23-Apr-2006 : Filter for directional spread.      ( version 3.09 )

    !/    02-Apr-2007 : Adding partitioned output.          ( version 3.11 )

    !/                  Adding user slots for outputs.

    !/    08-Oct-2007 : Adding ST3 source term option.      ( version 3.13 )

    !/                  ( F. Ardhuin )

    !/    05-Mar-2008 : Added NEC sxf90 compiler directives

    !/                  (Chris Bunney, UK Met Office)       ( version 3.13 )

    !/    25-Dec-2012 : New output structure and smaller    ( version 4.11 )

    !/                  memory footprint.

    !/    10-Feb-2014 : Bug correction for US3D: div. by df ( version 4.18 )

    !/    30-Apr-2014 : Add th2m and sth2m calculation      ( version 5.01 )

    !/    27-May-2014 : Switch to OMPG switch.              ( version 5.02 )

    !/    10-Apr-2015 : Remove unused variables             ( version 5.08 )

    !/    10-Jan-2017 : Separate Stokes drift calculation   ( version 6.01 )

    !/    01-Mar-2018 : Removed RTD code (now used in post  ( version 6.02 )

    !/                  processing code)

    !/    22-Aug-2018 : Add WBT parameter                   ( version 6.06 )

    !/    25-Sep-2019 : Corrected th2m and sth2m            ( version 6.07 )

    !/                  calculations. (J Dykes, NRL)

    !/

    !  1. Purpose :

    !

    !     Fill necessary arrays with gridded data for output.

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !       A       R.A.   I   Input spectra. Left in par list to change

    !                          shape.

    !       FLPART  Log.   I   Flag for filling fields with part. data.

    !       FLOUTG  Log.   I   Flag for file field output

    !       FLOUTG2 Log.   I   Flag for coupling field output

    !     ----------------------------------------------------------------

    !

    !     Locally saved parameters

    !     ----------------------------------------------------------------

    !       HSMIN   Real  Filter level in Hs for calculation of mean

    !                     wave parameters.

    !     ----------------------------------------------------------------

    !

    !  4. Subroutines used :

    !

    !     See module documentation.

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !     None.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/SHRD  Switch for shared / distributed memory architecture.

    !     !/DIST  Id.

    !

    !     !/OMPG  OpenMP compiler directive for loop splitting.

    !

    !     !/O8    Filter for low wave heights ( HSMIN )

    !     !/O9    Negative wave height alowed, other mean parameters will

    !             not be correct.

    !

    !     !/ST0   No source terms.

    !     !/ST1   Source term set 1 (WAM equiv.)

    !     !/ST2   Source term set 2 (Tolman and Chalikov)

    !     !/ST3   Source term set 3 (WAM 4+)

    !     !/ST6   Source term set 6 (BYDRZ)

    !

    !     !/S     Enable subroutine tracing.

    !     !/T     Test output.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

    USE constants

    USE w3gdatmd

    USE w3wdatmd, ONLY: ust, fpis

    USE w3adatmd, ONLY: cg, wn, dw

    USE w3adatmd, ONLY: hs, wlm, t02, t0m1, t01, fp0,               &

         thm, ths, thp0

    USE w3adatmd, ONLY: aba, abd, uba, ubd, fcut, sxx,              &

         syy, sxy, phs, ptp, plp, pdir, psi, pws,    &

         pwst, pnr, usero, tusx, tusy, prms, tpms,   &

         ussx, ussy, mssx, mssy, mssd, mscx, mscy,   &

         mscd, charn,                                &

         bhd, cge, p2sms, us3d, ef, th1m, sth1m,     &

         th2m, sth2m, hsig, stmaxe, stmaxd,          &

         hcmaxe, hmaxe, hcmaxd, hmaxd, ussp, qp, pqp,&

         pthp0, ppe, pgw, psw, ptm1, pt1, pt2, pep,  &

         wbt, qkk

    USE w3odatmd, ONLY: ndst, undef, iaproc, naproc, napfld,        &

         icprt, dtprt, wscut, noswll, flogrd, flogr2,&

         nogrp, ngrpp

    USE w3adatmd, ONLY: nsealm

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    USE w3parall, ONLY : init_get_isea

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    REAL, INTENT(IN)        :: A(NTH,NK,0:NSEAL)

    LOGICAL, INTENT(IN)     :: FLPART, FLOUTG, FLOUTG2

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

    INTEGER                 :: IK, ITH, JSEA, ISEA, IX, IY,         &

         IKP0(NSEAL), NKH(NSEAL),             &

         I, J, LKMS, HKMS, ITL

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    REAL                    :: FXPMC, FACTOR, FACTOR2, EBAND, FKD,  &

         AABS, UABS,                          &

         XL, XH, XL2, XH2, EL, EH, DENOM, KD, &

         M1, M2, MA, MB, MC, STEX, STEY, STED

    REAL                    :: ET(NSEAL), EWN(NSEAL), ETR(NSEAL),   &

         ETX(NSEAL), ETY(NSEAL), AB(NSEAL),   &

         ETXX(NSEAL), ETYY(NSEAL), ETXY(NSEAL),&

         ABX(NSEAL), ABY(NSEAL),ET02(NSEAL),  &

         EBD(NK,NSEAL), EC(NSEAL),            &

         ABR(NSEAL), UBR(NSEAL), UBS(NSEAL),  &

         ABX2(NSEAL), ABY2(NSEAL),            &

         AB2X(NSEAL), AB2Y(NSEAL),            &

         ABST(NSEAL), ABXX(NSEAL),            &

         ABYY(NSEAL), ABXY(NSEAL),            &

         ABYX(NSEAL), EET1(NSEAL),            &

         ETUSCX(NSEAL), ETUSCY(NSEAL),        &

         ETMSSL(NSEAL), ETMSSCL(NSEAL),       &

         ETTPMM(NSEAL), ETF(NSEAL),           &

         ET1(NSEAL), ABX2M(NSEAL),            &

         ABY2M(NSEAL), ABXM(NSEAL),           &

         ABYM(NSEAL), ABXYM(NSEAL),           &

         MSSXM(NSEAL), MSSYM(NSEAL),          &

         MSSXTM(NSEAL), MSSYTM(NSEAL),        &

         MSSXYM(NSEAL), THMP(NSEAL),          &

         T02P(NSEAL), NV(NSEAL), NS(NSEAL),   &

         NB(NSEAL), MODE(NSEAL),              &

         MU(NSEAL), NI(NSEAL), STMAXEL(NSEAL),&

         PHI(21,NSEAL),PHIST(NSEAL),         &

         EBC(NK,NSEAL), ABP(NSEAL),           &

         STMAXDL(NSEAL), TLPHI(NSEAL),        &

         WL02X(NSEAL), WL02Y(NSEAL),          &

         ALPXT(NSEAL), ALPYT(NSEAL),          &

         ALPXY(NSEAL), SCREST(NSEAL),         &

         QK1(NSEAL), QK2(NSEAL)

    REAL                       USSCO, FT1

    REAL, SAVE              :: HSMIN = 0.01

    LOGICAL                 :: FLOLOC(NOGRP,NGRPP)

    !/

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'W3OUTG')

#endif

    DO i=1,nogrp

      DO j=1,ngrpp

        floloc(i,j) =   &

             ((floutg.AND.flogrd(i,j)).OR.(floutg2.AND.flogr2(i,j)))

      END DO

    END DO

    !

    fxpmc  = 0.66 * grav / 28.

    hsmin  = hsmin

    ft1    =  0.3333 * sig(nk)**2 * dth * sig(nk)

    !

    ! 1.  Initialize storage arrays -------------------------------------- *

    !

    et     = 0.

    et02   = 0.

    ewn    = 0.

    etr    = 0.

    et1    = 0.

    eet1   = 0.

    etx    = 0.

    ety    = 0.

    etxx   = 0.

    etyy   = 0.

    etxy   = 0.

    abr    = 0.

    aba    = 0.

    abd    = 0.

    ubr    = 0.

    uba    = 0.

    ubd    = 0.

    ubs    = 0.

    sxx    = 0.

    syy    = 0.

    sxy    = 0.

    ussx   = 0.

    ussy   = 0.

    tusx   = 0.

    tusy   = 0.

    mssx   = 0.

    mssy   = 0.

    mssd   = 0.

    mscx   = 0.

    mscy   = 0.

    mscd   = 0.

    prms   = 0.

    tpms   = 0.

    etuscy = 0.

    etuscy = 0.

    etmssl = 0.

    etmsscl= 0.

    ettpmm = 0.

    ebd    = 0.

    ec     = 0.

    etf    = 0.

    ebc    = 0.

    bhd = 0.

    mssxm = 0.

    mssym = 0.

    mssxtm = 0.

    mssytm = 0.

    mssxym = 0.

    phi    = 0.

    phist  = 0.

    tlphi  = 0.

    stmaxel = 0.

    stmaxdl = 0.

    qk2 = 0.

    !

    hs     = undef

    wlm    = undef

    t0m1   = undef

    t01    = undef

    t02    = undef

    fp0    = undef

    thm    = undef

    ths    = undef

    thp0   = undef

    hsig   = undef

    wl02x  = undef

    wl02y  = undef

    alpxy  = undef

    alpxt  = undef

    alpyt  = undef

    qkk    = undef

    thmp = undef

    t02p = undef

    screst = undef

    nv = undef

    ns = undef

    nb = undef

    mu = undef

    ni = undef

    mode = undef

    stmaxe = undef

    stmaxd = undef

    hcmaxe = undef

    hmaxe = undef

    hcmaxd = undef

    hmaxd = undef

    qp    = undef

    wbt    = undef

    !

    ! 2.  Integral over discrete part of spectrum ------------------------ *

    !

    DO ik=1, nk

      !

      ! 2.a Initialize energy in band

      !

      ab     = 0.

      abx    = 0.

      aby    = 0.

      abx2   = 0.

      aby2   = 0.

      ab2x   = 0.

      ab2y   = 0.

      abxx   = 0.

      abyy   = 0.

      abxy   = 0.

      abyx   = 0.

      abst   = 0.

      qk1    = 0.

      !

      ! 2.b Integrate energy in band

      !

      DO ith=1, nth

        !

#ifdef W3_OMPG

        !$OMP PARALLEL DO PRIVATE(JSEA,ISEA,FACTOR)

#endif

        !

        DO jsea=1, nseal

          nkh(jsea)  = min( nk ,   &

               int(facti2+facti1*log(max(1.e-7,fcut(jsea)))) )

          ab(jsea)  = ab(jsea) + a(ith,ik,jsea)

          abx(jsea)  = abx(jsea) + a(ith,ik,jsea)*ecos(ith)

          aby(jsea)  = aby(jsea) + a(ith,ik,jsea)*esin(ith)

          ! These are the integrals with cos^2 and sin^2

          abx2(jsea) = abx2(jsea) + a(ith,ik,jsea)*ec2(ith)

          aby2(jsea) = aby2(jsea) + a(ith,ik,jsea)*es2(ith)

          ! Using trig identities to represent cos2theta and sin2theta.

          ab2x(jsea) = ab2x(jsea) + a(ith,ik,jsea)*(2*ec2(ith) - 1)

          ab2y(jsea) = ab2y(jsea) + a(ith,ik,jsea)*(2*esc(ith))

          abyx(jsea) = abyx(jsea) + a(ith,ik,jsea)*esc(ith)

          IF (ith.LE.nth/2) THEN

            abst(jsea) = abst(jsea) +                               &

                 a(ith,ik,jsea)*a(ith+nth/2,ik,jsea)

            qk1(jsea) = qk1(jsea) + (a(ith,ik,jsea)+a(ith+nth/2,ik,jsea))**2

          END IF

          CALL init_get_isea(isea, jsea)

          factor     = max( 0.5 , cg(ik,isea)/sig(ik)*wn(ik,isea) )

          abxx(jsea) = abxx(jsea) + ((1.+ec2(ith))*factor-0.5) *    &

               a(ith,ik,jsea)

          abyy(jsea) = abyy(jsea) + ((1.+es2(ith))*factor-0.5) *    &

               a(ith,ik,jsea)

          abxy(jsea) = abxy(jsea) + esc(ith)*factor * a(ith,ik,jsea)

        END DO

        !

#ifdef W3_OMPG

        !$OMP END PARALLEL DO

#endif

        !

      END DO

      !

      ! 2.c Finalize integration over band and update mean arrays

      !

      !

#ifdef W3_OMPG

      !$OMP PARALLEL DO PRIVATE(JSEA,ISEA,FACTOR,FACTOR2,MA,MC,MB,KD,FKD,USSCO,M1,M2)

#endif

      !

      DO jsea=1, nseal

        CALL init_get_isea(isea, jsea)

        factor       = dden(ik) / cg(ik,isea)

        ebd(ik,jsea) = ab(jsea) * factor            ! this is E(f)*df

        et(jsea)    = et(jsea) + ebd(ik,jsea)

#ifdef W3_IG1

        IF (ik.EQ.nint(igpars(5))) hsig(jsea) = 4*sqrt(et(jsea))

#endif

        etf(jsea)  = etf(jsea) + ebd(ik,jsea) * cg(ik,isea)

        ewn(jsea)  = ewn(jsea) + ebd(ik,jsea) / wn(ik,isea)

        etr(jsea)  = etr(jsea) + ebd(ik,jsea) / sig(ik)

        et1(jsea)  = et1(jsea) + ebd(ik,jsea) * sig(ik)

        !          EET1(JSEA) = EET1(JSEA)+ EBD(IK,JSEA)**2 * SIG(IK)

        eet1(jsea) = eet1(jsea)+ ebd(ik,jsea)**2 * sig(ik)/dsii(ik)

        et02(jsea) = et02(jsea)+ ebd(ik,jsea) * sig(ik)**2

        etx(jsea)  = etx(jsea) + abx(jsea) * factor

        ety(jsea)  = ety(jsea) + aby(jsea) * factor

        tusx(jsea) = tusx(jsea) + abx(jsea)*factor               &

             *grav*wn(ik,isea)/sig(ik)

        tusy(jsea)  = tusy(jsea) + aby(jsea)*factor               &

             *grav*wn(ik,isea)/sig(ik)

        etxx(jsea) = etxx(jsea) + abx2(jsea) * factor* wn(ik,isea)**2

        ! NB:     QK1 (JSEA) = QK1(JSEA) + A(ITH,IK,JSEA)**2

        qk2(jsea) = qk2(jsea) + qk1(jsea)  * factor* sig(ik) /wn(ik,isea)

        etyy(jsea) = etyy(jsea) + aby2(jsea) * factor* wn(ik,isea)**2

        etxy(jsea) = etxy(jsea) + abyx(jsea) * factor* wn(ik,isea)**2

        IF (sig(ik)*0.5*(1+xfr).LT.0.4*tpi) THEN

          etmssl(jsea)  = etmssl(jsea) + ab(jsea)*factor           &

               *wn(ik,isea)**2

        ELSE

          IF (sig(max(ik-1,1))*0.5*(1+xfr).LT.0.4*tpi) THEN

            etmssl(jsea)  = etmssl(jsea) + ab(jsea)*factor         &

                 *(sig(ik)*0.5*(1+1/xfr)-(0.4*tpi))/dsii(ik)     &

                 *wn(ik,isea)**2

            factor2       = sig(ik)**5/(grav**2)/dsii(ik)

            etmsscl(jsea) = ab(jsea)*factor*factor2

          END IF

        END IF

        !

        ubs(jsea) = ubs(jsea) + ab(jsea) * sig(ik)**2

        !

        !   2nd order equivalent surface pressure spectral density at K=0

        !   this is used for microseismic or microbarom sources

        !   Finite water depth corrections (Ardhuin & Herbers 2013) are not

        !   included here.

        !

        factor2 = dth*2/(tpi**2)                        &

             * sig(ik)                             &

             * (tpi*sig(ik)/cg(ik,isea))**2        &  ! Jacobian^2 to get E(f,th) from A(k,th)

             * abst(jsea)

        !

        !   Integration over seismic radian frequency : *2*dsigma

        !

        prms(jsea)  = prms(jsea) + factor2 * 2 * dsii(ik)

        IF ( floloc(6, 9).AND.(ik.GE.p2msf(2).AND.ik.LE.p2msf(3)))   &

             p2sms(jsea,ik) = factor2 * 2 * tpi

        IF (factor2 .GT. ettpmm(jsea)) THEN

          ettpmm(jsea) = factor2

          tpms(jsea) = tpi/sig(ik)

        END IF


        !

        ! Directional moments in the last freq. band

        !

        IF (ik.EQ.nk) THEN

          factor2       = sig(ik)**5/(grav**2)/dsii(ik)

          etuscx(jsea)  = abx(jsea)*factor*factor2

          etuscy(jsea)  = aby(jsea)*factor*factor2

          !

          !     NB: the slope PDF is proportional to ell1=ETYY*EC2-2*ETXY*ECS+ETYY*ES2 = A*EC2-2*B*ECS+C*ES2

          !     This is an ellipse equation with axis direction given by dir=0.5*ATAN2(-2.*ETXY,ETYY-ETXX)

          !

          ma  = abx2(jsea) * factor * factor2

          mc  = aby2(jsea) * factor * factor2

          mb  = abyx(jsea) * factor * factor2

          !

          ! Old definitions:  MSCX(JSEA)  = ABX2(JSEA) * FACTOR * FACTOR2

          !                   MSCY(JSEA)  = ABY2(JSEA) * FACTOR * FACTOR2

          mscd(jsea)=0.5*atan2(2*mb,ma-mc)


          mscx(jsea)= ma*cos(mscd(jsea))**2   &

               +2*mb*sin(mscd(jsea))*cos(mscd(jsea))+ma*sin(mscd(jsea))**2

          mscy(jsea)= mc*cos(mscd(jsea))**2   &

               -2*mb*sin(mscd(jsea))*cos(mscd(jsea))+ma*sin(mscd(jsea))**2

        END IF

        !

        ! Deep water limits

        !

        kd    = max( 0.001 , wn(ik,isea) * dw(isea) )

        IF ( kd .LT. 6. ) THEN

          fkd       = factor / sinh(kd)**2

          abr(jsea) = abr(jsea) + ab(jsea) * fkd

          aba(jsea) = aba(jsea) + abx(jsea) * fkd

          abd(jsea) = abd(jsea) + aby(jsea) * fkd

          ubr(jsea) = ubr(jsea) + ab(jsea) * sig(ik)**2 * fkd

          uba(jsea) = uba(jsea) + abx(jsea) * sig(ik)**2 * fkd

          ubd(jsea) = ubd(jsea) + aby(jsea) * sig(ik)**2 * fkd

          ussco=fkd*sig(ik)*wn(ik,isea)*cosh(2.*kd)

          bhd(jsea) = bhd(jsea) +                             &

               grav*wn(ik,isea) * ebd(ik,jsea) / (sinh(2.*kd))

        ELSE

          ussco=factor*sig(ik)*2.*wn(ik,isea)

        END IF

        !

        abxx(jsea)   = max( 0. , abxx(jsea) ) * factor

        abyy(jsea)   = max( 0. , abyy(jsea) ) * factor

        abxy(jsea)   = abxy(jsea) * factor

        sxx(jsea)    = sxx(jsea)  + abxx(jsea)

        syy(jsea)    = syy(jsea)  + abyy(jsea)

        sxy(jsea)    = sxy(jsea)  + abxy(jsea)

        ebd(ik,jsea) = ebd(ik,jsea) / dsii(ik)

        !

        IF ( floloc( 3, 1).AND.(ik.GE.e3df(2,1).AND.ik.LE.e3df(3,1)))   &

             ef(jsea,ik)  = ebd(ik,jsea) * tpi

        !

        ussx(jsea)  = ussx(jsea) + abx(jsea)*ussco

        ussy(jsea)  = ussy(jsea) + aby(jsea)*ussco

        !

        ! Fills the 3D Stokes drift spectrum array

        !  ! The US3D Stokes drift specrum array is now calculated in a

        !  subroutine and called at the end of this subroutine

        !          IF ( FLOLOC( 6, 8).AND.(IK.GE.US3DF(2).AND.IK.LE.US3DF(3) ))   THEN

        !            US3D(JSEA,IK)    =  ABX(JSEA)*USSCO/(DSII(IK)*TPIINV)

        !            US3D(JSEA,NK+IK) =  ABY(JSEA)*USSCO/(DSII(IK)*TPIINV)

        !          END IF

        IF ( floloc( 3, 2).AND.(ik.GE.e3df(2,2).AND.ik.LE.e3df(3,2)))  &

             th1m(jsea,ik)= mod( 630. - rade*atan2(aby(jsea),abx(jsea)) , 360. )

        m1 = sqrt(abx(jsea)**2+aby(jsea)**2)/max(1e-20,ab(jsea))

        IF ( floloc( 3, 3).AND.(ik.GE.e3df(2,3).AND.ik.LE.e3df(3,3)))  &

             sth1m(jsea,ik)= sqrt(abs(2.*(1-m1)))*rade

        IF ( floloc( 3, 4).AND.(ik.GE.e3df(2,4).AND.ik.LE.e3df(3,4)))  &

             th2m(jsea,ik)= mod( 270. - rade*0.5*atan2(aby2(jsea),ab2x(jsea)) , 180. )

        m2 = sqrt(ab2x(jsea)**2+ab2y(jsea)**2)/max(1e-20,ab(jsea))

        IF ( floloc( 3, 5).AND.(ik.GE.e3df(2,5).AND.ik.LE.e3df(3,5)))  &

             sth2m(jsea,ik)= sqrt(abs(0.5*(1-m2)))*rade

      END DO

      !

#ifdef W3_OMPG

      !$OMP END PARALLEL DO

#endif

      !

    END DO

    !

    ! Start of Space-Time Extremes Section

    IF ( ( stexu .GT. 0. .AND. steyu .GT. 0. ) &

         .OR. ( stedu .GT. 0. ) ) THEN

      !  Space-Time extremes

      !    (for references:

      !     - Krogstad et al, OMAE 2004

      !     - Baxevani and Rychlik, OE 2006

      !     - Adler and Taylor, 2007

      !     - Fedele, JPO 2012

      !     - Fedele et al, OM 2013

      !     - Benetazzo et al, JPO 2015)

      !

      !  Compute spectral parameters wrt the mean wave direction

      !  (no tail contribution - Prognostic)

      DO jsea=1, nseal

        CALL init_get_isea(isea, jsea)

        ix     = mapsf(isea,1)

        iy     = mapsf(isea,2)

        IF ( mapsta(iy,ix) .GT. 0 ) THEN

          IF ( abs(etx(jsea))+abs(ety(jsea)) .GT. 1.e-7 ) THEN

            thmp(jsea) = atan2(ety(jsea),etx(jsea))

          END IF

        END IF

      END DO

      !

      DO ik=1, nk

        !

        abx2m = 0.

        aby2m = 0.

        abxm = 0.

        abym = 0.

        abxym = 0.

        !

        DO ith=1, nth

          !

#ifdef W3_OMPG

          !$OMP PARALLEL DO PRIVATE(JSEA)

#endif

          !

          DO jsea=1, nseal

            abx2m(jsea) = abx2m(jsea) + a(ith,ik,jsea)*                &

                 (ecos(ith)*cos(thmp(jsea))+esin(ith)*sin(thmp(jsea)))**2

            aby2m(jsea) = aby2m(jsea) + a(ith,ik,jsea)*                &

                 (esin(ith)*cos(thmp(jsea))-ecos(ith)*sin(thmp(jsea)))**2

            abxm(jsea)  = abxm(jsea) + a(ith,ik,jsea)*                 &

                 (ecos(ith)*cos(thmp(jsea))+esin(ith)*sin(thmp(jsea)))

            abym(jsea)  = abym(jsea) + a(ith,ik,jsea)*                 &

                 (esin(ith)*cos(thmp(jsea))-ecos(ith)*sin(thmp(jsea)))

            abxym(jsea) = abxym(jsea) + a(ith,ik,jsea)*                &

                 (ecos(ith)*cos(thmp(jsea))+esin(ith)*sin(thmp(jsea)))*   &

                 (esin(ith)*cos(thmp(jsea))-ecos(ith)*sin(thmp(jsea)))

          END DO

          !

#ifdef W3_OMPG

          !$OMP END PARALLEL DO

#endif

          !

        END DO

        !

#ifdef W3_OMPG

        !$OMP PARALLEL DO PRIVATE(JSEA,ISEA,FACTOR)

#endif

        !

        DO jsea=1, nseal

          CALL init_get_isea(isea, jsea)

          factor       = dden(ik) / cg(ik,isea)

          mssxm(jsea)  = mssxm(jsea) + abx2m(jsea)*factor*             &

               wn(ik,isea)**2

          mssym(jsea)  = mssym(jsea) + aby2m(jsea)*factor*             &

               wn(ik,isea)**2

          mssxtm(jsea)  = mssxtm(jsea) + abxm(jsea)*factor*wn(ik,isea)* &

               sig(ik)

          mssytm(jsea)  = mssytm(jsea) + abym(jsea)*factor*wn(ik,isea)* &

               sig(ik)

          mssxym(jsea)  = mssxym(jsea) + abxym(jsea)*factor*           &

               wn(ik,isea)**2

        END DO

        !

#ifdef W3_OMPG

        !$OMP END PARALLEL DO

#endif

        !

      END DO


#ifdef W3_OMPG

      !$OMP PARALLEL DO PRIVATE(JSEA,STEX,STEY,STED,ITL,IK)

#endif

      !

      DO jsea=1, nseal

        !

        !  Mean wave period (no tail contribution - Prognostic)

        IF ( et02(jsea) .GT. 1.e-7 ) THEN

          t02p(jsea) = tpi * sqrt(et(jsea) / et02(jsea) )

        END IF

        !

        !  Mean wavelength and mean crest length (02) for space-time extremes

        IF ( mssxm(jsea) .GT. 1.e-7 ) THEN

          wl02x(jsea) = tpi * sqrt(et(jsea) / mssxm(jsea))

        END IF

        IF ( mssym(jsea) .GT. 1.e-7 ) THEN

          wl02y(jsea) = tpi * sqrt(et(jsea) / mssym(jsea))

        END IF

        !

        !  Irregularity parameters for space-time extremes

        IF ((mssxm(jsea) .GT. 1.e-7) .AND. (et02(jsea) .GT. 1.e-7)) THEN

          alpxt(jsea) = mssxtm(jsea) / (sqrt(mssxm(jsea) * et02(jsea)))

        ENDIF

        IF ((mssym(jsea) .GT. 1.e-7) .AND. (et02(jsea) .GT. 1.e-7)) THEN

          alpyt(jsea) = mssytm(jsea) / (sqrt(mssym(jsea) * et02(jsea)))

        ENDIF

        IF ((mssxm(jsea) .GT. 1.e-7) .AND. (mssym(jsea) .GT. 1.e-7)) THEN

          alpxy(jsea) = mssxym(jsea) / (sqrt(mssxm(jsea) * mssym(jsea)))

        ENDIF

        !

        !  Short-crestedness parameter

        IF (mssxm(jsea) .GT. 1.e-7)  THEN

          screst(jsea) = sqrt(mssym(jsea)/mssxm(jsea))

        END IF

        !

        !  Space domain size (user-defined or default)

        IF ( stexu .GT. 0 .AND. steyu .GT. 0 ) THEN

          stex = stexu

          stey = steyu

        ELSE

          stex = 0.

          stey = 0.

        END IF

        !

        !  Time domain size (user-defined or default)

        IF ( stedu .GT. 0 ) THEN

          sted = stedu

        ELSE

          sted = 0.

        END IF

        !

        !  Average numbers of waves in the space-time domain (Volume+Sides+Borders)

        IF ((wl02x(jsea) .GT. 1.e-7) .AND. (wl02y(jsea) .GT. 1.e-7)    &

             .AND. (t02p(jsea) .GT. 1.e-7)) THEN

          nv(jsea) = tpi*(stex*stey*sted)/                             &

               (wl02x(jsea)*wl02y(jsea)*t02p(jsea))  *                    &

               sqrt(1-alpxt(jsea)**2-alpyt(jsea)**2  -                    &

               alpxy(jsea)**2+2*alpxt(jsea)*alpyt(jsea)*alpxy(jsea))

          ns(jsea) = sqrt(tpi)*((stex*sted)/(wl02x(jsea)*t02p(jsea)) * &

               sqrt(1-alpxt(jsea)**2) +                                   &

               (stey*sted)/(wl02y(jsea)*t02p(jsea)) *                     &

               sqrt(1-alpyt(jsea)**2) +                                   &

               (stex*stey)/(wl02x(jsea)*wl02y(jsea)) *                    &

               sqrt(1-alpxy(jsea)**2))

          nb(jsea) = stex/wl02x(jsea) + stey/wl02y(jsea) +             &

               sted/t02p(jsea)

        END IF

        !

        ! Integral measure of wave steepness (Fedele & Tayfun, 2009) MU, as a

        ! function of the spectral width parameter NI (Longuet-Higgins, 1985)

        IF (et1(jsea) .GT. 1.e-7) THEN

          ni(jsea) = sqrt(et(jsea)*et02(jsea)/et1(jsea)**2 - 1)

        ENDIF

        IF (et(jsea) .GT. 1.e-7) THEN

          mu(jsea) = et1(jsea)**2/grav * (et(jsea))**(-1.5) *          &

               (1-ni(jsea)+ni(jsea)**2)

        ENDIF

        !

        ! Mode of the Adler&Taylor distribution

        ! (normalized on the standard deviation = Hs/4)

        ! Time extremes

        IF ((stex .EQ. 0) .AND. (stey .EQ. 0)) THEN

          mode(jsea) = sqrt(2.*log(nb(jsea)))

          ! Space extremes (strictly for STEX*STEY >> WL02X*WL02Y)

        ELSEIF (sted .EQ. 0) THEN

          mode(jsea) = sqrt(2.*log(ns(jsea))+log(2.*log(ns(jsea))+     &

               log(2.*log(ns(jsea)))))

          ! Space-time extremes (strictly for STEX*STEY >> WL02X*WL02Y)

        ELSEIF ((wl02x(jsea) .GT. 1.e-7) .AND. (wl02y(jsea) .GT. 1.e-7) &

             .AND. (t02p(jsea) .GT. 1.e-7)) THEN

          mode(jsea) = sqrt(2.*log(nv(jsea))+2.*log(2.*log(nv(jsea))+  &

               2.*log(2.*log(nv(jsea)))))

        ENDIF

        !

        ! Expected maximum sea surface elevation in the ST domain - nonlinear

        ! (in meters, Hs/4=SQRT(ET(JSEA)))

        stmaxe(jsea) = sqrt(et(jsea)) *                                &

             ( mode(jsea)+0.5*mu(jsea)*mode(jsea)**2 +                  &

             0.5772*(1+mu(jsea)*mode(jsea)) /                           &

             (mode(jsea)-(2*nv(jsea)*mode(jsea)+ns(jsea)) /             &

             (nv(jsea)*mode(jsea)**2+ns(jsea)*mode(jsea)+nb(jsea))) )

        !

        ! Standard deviation of the maximum sea surface elevation in ST domain

        !  - nonlinear (in meters, Hs/4=SQRT(ET(JSEA)))

        stmaxd(jsea) =  sqrt(et(jsea)) *                               &

             ( pi*(1+mu(jsea)*mode(jsea))/sqrt(6.) /                    &

             (mode(jsea)-(2*nv(jsea)*mode(jsea)+ns(jsea)) /             &

             (nv(jsea)*mode(jsea)**2+ns(jsea)*mode(jsea)+nb(jsea))) )

        !

        ! Autocovariance (time) function (normalized on the maximum, i.e. total

        ! variance)

        IF (t02p(jsea) .GT. 1.e-7) THEN

          tlphi(jsea) = 0.3*t02p(jsea)

          DO itl = 1, 21

            DO ik = 1, nk-3, 4

              phi(itl,jsea) = phi(itl,jsea) +                          &

                   (xfr**3*ebd(ik+3,jsea)*cos(xfr**3*sig(ik)*tlphi(jsea))+    &

                   xfr**2*ebd(ik+2,jsea)*cos(xfr**2*sig(ik)*tlphi(jsea))+     &

                   xfr*ebd(ik+1,jsea)*cos(xfr*sig(ik)*tlphi(jsea)) +          &

                   ebd(ik,jsea)*cos(sig(ik)*tlphi(jsea)))*dsii(ik)

            ENDDO

            tlphi(jsea) = tlphi(jsea) + t02p(jsea)/20.

          ENDDO

          phi(:,jsea) = phi(:,jsea)/et(jsea)

          !

          ! First minimum of the autocovariance function (absolute value)

          phist(jsea) = abs(minval(phi(:,jsea),1))

        ENDIF

        !

        ! Wave height of the wave with the maximum expected crest height

        ! and corresponding standard deviation

        ! (according to Boccotti Quasi-Determinism theory - linear)

        stmaxel(jsea) = sqrt(et(jsea)) * ( mode(jsea)+0.5772 /         &

             (mode(jsea)-(2*nv(jsea)*mode(jsea)+ns(jsea)) /             &

             (nv(jsea)*mode(jsea)**2+ns(jsea)*mode(jsea)+nb(jsea))) )

        stmaxdl(jsea) = sqrt(et(jsea)) *                               &

             ( pi/sqrt(6.) /                                            &

             (mode(jsea)-(2*nv(jsea)*mode(jsea)+ns(jsea)) /             &

             (nv(jsea)*mode(jsea)**2+ns(jsea)*mode(jsea)+nb(jsea))) )

        hcmaxe(jsea) = stmaxel(jsea)*(1+phist(jsea))

        hcmaxd(jsea) = stmaxdl(jsea)*(1+phist(jsea))

        ! Maximum expected wave height and corresponding standard deviation

        ! (according to Boccotti Quasi-Determinism theory - linear)

        hmaxe(jsea) = stmaxel(jsea)*sqrt(2*(1+phist(jsea)))

        hmaxd(jsea) = stmaxdl(jsea)*sqrt(2*(1+phist(jsea)))

      ENDDO

      !

#ifdef W3_OMPG

      !$OMP END PARALLEL DO

#endif

      !


      ! End of Space-Time Extremes Section

    ENDIF

    !

    ! 3.  Finalize computation of mean parameters ------------------------ *

    !

#ifdef W3_OMPG

    !$OMP PARALLEL DO PRIVATE(JSEA,ISEA,EBAND)

#endif

    !

    DO jsea=1, nseal

      CALL init_get_isea(isea, jsea)

      !

      ! 3.a Directional mss parameters

      !     NB: the slope PDF is proportional to ell1=ETYY*EC2-2*ETXY*ECS+ETXX*ES2 = C*EC2-2*B*ECS+A*ES2

      !     This is an ellipse equation with axis direction given by dir=0.5*ATAN2(2.*ETXY,ETXX-ETYY)

      !     From matlab script: t0=0.5*(atan2(2.*B,A-C));

      !     From matlab script: A2=A.*cos(t0).^2+2.*B.*sin(t0).*cos(t0)+A.*cos(t0).^2+C.*sin(t0)^2;

      !     From matlab script: C2=C.*cos(t0)^2-2.*B.*sin(t0).*cos(t0)+A.*sin(t0).^2;

      mssd(jsea)=0.5*(atan2(2*etxy(jsea),etxx(jsea)-etyy(jsea)))

      mssx(jsea)  = etxx(jsea)*cos(mssd(jsea))**2   &

           +2*etxy(jsea)*sin(mssd(jsea))*cos(mssd(jsea))+etyy(jsea)*sin(mssd(jsea))**2

      mssy(jsea)  = etyy(jsea)*cos(mssd(jsea))**2   &

           -2*etxy(jsea)*sin(mssd(jsea))*cos(mssd(jsea))+etxx(jsea)*sin(mssd(jsea))**2

      !

      ! 3.b Add tail

      !     ( DTH * SIG absorbed in FTxx )


      eband     = ab(jsea) / cg(nk,isea)           ! EBAND is E(sigma)/sigma for the last frequency band

      et(jsea) = et(jsea) + fte  * eband

      ewn(jsea) = ewn(jsea) + ftwl * eband

      etf(jsea) = etf(jsea) + grav * fttr * eband  ! this is the integral of CgE in deep water

      etr(jsea) = etr(jsea) + fttr * eband

      et1(jsea) = et1(jsea) + ft1  * eband

      !        EET1(JSEA)= EET1(JSEA) + FT1  * EBAND**2 : this was not correct. Actually tail may not be needed for Qp.

      et02(jsea)= et02(jsea)+ eband* 0.5 * sig(nk)**4 * dth

      etx(jsea) = etx(jsea) + fte * abx(jsea) / cg(nk,isea)

      ety(jsea) = ety(jsea) + fte * aby(jsea) / cg(nk,isea)

      sxx(jsea) = sxx(jsea) + fte * abxx(jsea) / cg(nk,isea)

      syy(jsea) = syy(jsea) + fte * abyy(jsea) / cg(nk,isea)

      sxy(jsea) = sxy(jsea) + fte * abxy(jsea) / cg(nk,isea)

      !

      ! Tail for surface stokes drift is commented out: very sensitive to tail power

      !

      !       USSX(JSEA)  = USSX(JSEA) + 2*GRAV*ETUSCX(JSEA)/SIG(NK)

      !       USSY(JSEA)  = USSY(JSEA) + 2*GRAV*ETUSCY(JSEA)/SIG(NK)

      ubs(jsea) = ubs(jsea) + ftwl * eband/grav

    END DO

    !

#ifdef W3_OMPG

    !$OMP END PARALLEL DO

#endif

    !

    sxx    = sxx * dwat * grav

    syy    = syy * dwat * grav

    sxy    = sxy * dwat * grav

    !

#ifdef W3_OMPG

    !$OMP PARALLEL DO PRIVATE(JSEA,ISEA,IX,IY)

#endif

    !

    DO jsea=1, nseal

      CALL init_get_isea(isea, jsea)

      ix     = mapsf(isea,1)

      iy     = mapsf(isea,2)

      IF ( mapsta(iy,ix) .GT. 0 ) THEN

#ifdef W3_O9

        IF ( et(jsea) .GE. 0. ) THEN

#endif

          hs(jsea) = 4. * sqrt( et(jsea) )

#ifdef W3_O9

        ELSE

          hs(jsea) = - 4. * sqrt( -et(jsea) )

        END IF

#endif

        IF ( et(jsea) .GT. 1.e-7 ) THEN

          qp(jsea) = ( 2. / et(jsea)**2 ) * eet1(jsea)

          wlm(jsea) = ewn(jsea) / et(jsea) * tpi

          t0m1(jsea) = etr(jsea) / et(jsea) * tpi

          ths(jsea) = rade * sqrt( max( 0. , 2. * ( 1. - sqrt( &

               max(0.,(etx(jsea)**2+ety(jsea)**2)/et(jsea)**2) ) ) ) )

          IF ( ths(jsea) .LT. 0.01*rade*dth ) ths(jsea) = 0.

          ! NB:           QK1 (JSEA) = QK1(JSEA) + A(ITH,IK,JSEA)**2

          !               QK2 (JSEA) = QK2 (JSEA) + QK1(JSEA)  * FACTOR* SIG(IK) /WN(IK,ISEA)

          qkk(jsea) = sqrt(0.5*qk2(jsea))/et(jsea)

        ELSE

          wlm(jsea) = 0.

          t0m1(jsea) = tpi / sig(nk)

          ths(jsea) = 0.

        END IF

        IF ( abs(etx(jsea))+abs(ety(jsea)) .GT. 1.e-7 ) THEN

          thm(jsea) = atan2(ety(jsea),etx(jsea))

        ELSE

          thm(jsea) = 0.

        END IF

        abr(jsea) = sqrt( 2. * max( 0. , abr(jsea) ) )

        IF ( abr(jsea) .GE. 1.e-7 ) THEN

          abd(jsea) = atan2(abd(jsea),aba(jsea))

        ELSE

          abd(jsea) = 0.

        ENDIF

        aba(jsea) = abr(jsea)

        ubr(jsea) = sqrt( 2. * max( 0. , ubr(jsea) ) )

        IF ( ubr(jsea) .GE. 1.e-7 ) THEN

          ubd(jsea) = atan2(ubd(jsea),uba(jsea))

        ELSE

          ubd(jsea) = 0.

        ENDIF

        uba(jsea) = ubr(jsea)

        cge(jsea) = dwat*grav*etf(jsea)

        IF ( et02(jsea) .GT. 1.e-7  .AND.  et(jsea) .GT. 0 ) THEN

          t02(jsea) = tpi * sqrt(et(jsea) / et02(jsea) )

          t01(jsea) = tpi * et(jsea) / et1(jsea)

        ELSE

          t02(jsea) = tpi / sig(nk)

          t01(jsea)= t02(jsea)

        ENDIF

        !

        !  Add here USERO(JSEA,1) ...

        !

      END IF

    END DO

    !

#ifdef W3_OMPG

    !$OMP END PARALLEL DO

#endif

    !

    ! 3.b Clean-up small values if !/O8 switch selected

    !

#ifdef W3_O8

    DO jsea=1, nseal

      IF ( hs(jsea).LE.hsmin .AND. hs(jsea).NE.undef) THEN

        wlm(jsea) = undef

        t02(jsea) = undef

        t0m1(jsea) = undef

        thm(jsea) = undef

        ths(jsea) = undef

      END IF

    END DO

#endif

    !

    ! 4.  Peak frequencies and directions -------------------------------- *

    ! 4.a Initialize

    !

#ifdef W3_OMPG

    !$OMP PARALLEL DO PRIVATE(JSEA)

#endif

    !

    DO jsea=1, nseal

      ec(jsea) = ebd(nk,jsea)

      fp0(jsea) = undef

      ikp0(jsea) = nk

      thp0(jsea) = undef

    END DO

    !

#ifdef W3_OMPG

    !$OMP END PARALLEL DO

#endif

    !

    ! 4.b Discrete peak frequencies

    !

    DO ik=nk-1, 1, -1

      !

#ifdef W3_OMPG

      !$OMP PARALLEL DO PRIVATE(JSEA)

#endif

      !

      DO jsea=1, nseal

        IF ( ec(jsea) .LT. ebd(ik,jsea) ) THEN

          ec(jsea) = ebd(ik,jsea)

          ikp0(jsea) = ik

        END IF

      END DO

      !

#ifdef W3_OMPG

      !$OMP END PARALLEL DO

#endif

      !

    END DO

    !

#ifdef W3_OMPG

    !$OMP PARALLEL DO PRIVATE(JSEA)

#endif

    !

    DO jsea=1, nseal

      IF ( ikp0(jsea) .NE. nk ) fp0(jsea) = sig(ikp0(jsea)) * tpiinv

    END DO

    !

#ifdef W3_OMPG

    !$OMP END PARALLEL DO

#endif

    !

    ! 4.c Continuous peak frequencies

    !

    xl     = 1./xfr - 1.

    xh     =  xfr - 1.

    xl2    = xl**2

    xh2    = xh**2

    !

#ifdef W3_OMPG

    !$OMP PARALLEL DO PRIVATE(JSEA,EL,EH,DENOM)

#endif

    !

    DO jsea=1, nseal

      IF ( ikp0(jsea) .NE. nk ) THEN

        IF ( ikp0(jsea) .EQ. 1 ) THEN

          el = - ebd(ikp0(jsea), jsea)

        ELSE

          el = ebd(ikp0(jsea)-1, jsea) - ebd(ikp0(jsea), jsea)

        END IF


        eh = ebd(ikp0(jsea)+1, jsea) - ebd(ikp0(jsea), jsea)


        denom  = xl*eh - xh*el

        fp0(jsea) = fp0(jsea) * ( 1. + 0.5 * ( xl2*eh - xh2*el )   &

             / sign( max(abs(denom),1.e-15) , denom ) )

      END IF

    END DO

    !

#ifdef W3_OMPG

    !$OMP END PARALLEL DO

#endif

    !

    ! 4.d Peak directions

    !

#ifdef W3_OMPG

    !$OMP PARALLEL DO PRIVATE(JSEA)

#endif

    !

    DO jsea=1, nseal

      etx(jsea) = 0.

      ety(jsea) = 0.

    END DO

    !

#ifdef W3_OMPG

    !$OMP END PARALLEL DO

#endif

    !

    DO ith=1, nth

      !

#ifdef W3_OMPG

      !$OMP PARALLEL DO PRIVATE(JSEA)

#endif

      !

      DO jsea=1, nseal

        IF ( ikp0(jsea) .NE. nk ) THEN

          etx(jsea) = etx(jsea) + a(ith,ikp0(jsea),jsea)*ecos(ith)

          ety(jsea) = ety(jsea) + a(ith,ikp0(jsea),jsea)*esin(ith)

        END IF

      END DO

      !

#ifdef W3_OMPG

      !$OMP END PARALLEL DO

#endif

      !

    END DO

    !

#ifdef W3_OMPG

    !$OMP PARALLEL DO PRIVATE(JSEA)

#endif

    !

    DO jsea=1, nseal

      IF ( abs(etx(jsea))+abs(ety(jsea)) .GT. 1.e-7 .AND.           &

           fp0(jsea).NE.undef )                                     &

           thp0(jsea) = atan2(ety(jsea),etx(jsea))

      etx(jsea) = 0.

      ety(jsea) = 0.

    END DO

    !

#ifdef W3_OMPG

    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO PRIVATE(JSEA,ISEA,IX,IY)

#endif

    !

    DO jsea =1, nseal

      CALL init_get_isea(isea, jsea)

      ix          = mapsf(isea,1)

      iy          = mapsf(isea,2)

      IF ( mapsta(iy,ix) .LE. 0 ) THEN

        fp0(jsea) = undef

        thp0(jsea) = undef

      END IF

    END DO

    !

#ifdef W3_OMPG

    !$OMP END PARALLEL DO

#endif

    !

    ! 5.  Test output (local to MPP only)

    !

#ifdef W3_T

    WRITE (ndst,9050)

    DO jsea =1, nseal

      CALL init_get_isea(isea, jsea)

      ix     = mapsf(isea,1)

      iy     = mapsf(isea,2)

      IF ( hs(jsea) .EQ. undef ) THEN

        WRITE (ndst,9051) isea, ix, iy

      ELSE IF ( wlm(jsea) .EQ. undef ) THEN

        WRITE (ndst,9052) isea, ix, iy, hs(jsea)

      ELSE IF ( fp0(jsea) .EQ. undef ) THEN

        WRITE (ndst,9053) isea, ix, iy, hs(jsea), wlm(jsea),   &

             t0m1(jsea), rade*thm(jsea), ths(jsea)

      ELSE

        WRITE (ndst,9054) isea, ix, iy, hs(jsea), wlm(jsea),   &

             t0m1(jsea), rade*thm(jsea), ths(jsea), fp0(jsea),&

             thp0(jsea)

      END IF

    END DO

#endif

    !

    ! 6.  Fill arrays wth partitioned data

    !

    IF ( flpart ) THEN

      !

      ! 6.a Initializations

      !

      phs    = undef

      ptp    = undef

      plp    = undef

      pdir   = undef

      psi    = undef

      pws    = undef

      pwst   = undef

      pnr    = undef

      pthp0  = undef

      pqp    = undef

      ppe    = undef

      pgw    = undef

      psw    = undef

      ptm1   = undef

      pt1    = undef

      pt2    = undef

      pep    = undef

      !

      ! 6.b Loop over local sea points

      !

#ifdef W3_OMPG

      !$OMP PARALLEL DO PRIVATE(ISEA,JSEA,IX,IY,I,J)

#endif

      !

      DO jsea=1, nseal

        CALL init_get_isea(isea, jsea)

        ix          = mapsf(isea,1)

        iy          = mapsf(isea,2)

        !

        IF ( mapsta(iy,ix).GT.0 ) THEN

          i         = icprt(jsea,2)

          pnr(jsea) = max( 0. , real(icprt(jsea,1)-1) )

          IF ( icprt(jsea,1).GE.1 ) pwst(jsea) = dtprt(6,i)

        END IF

        !

        IF ( mapsta(iy,ix).GT.0 .AND. icprt(jsea,1).GT.1 ) THEN

          i      = icprt(jsea,2) + 1

          IF ( dtprt(6,i) .GE. wscut ) THEN

            phs(jsea,0) = dtprt(1,i)

            ptp(jsea,0) = dtprt(2,i)

            plp(jsea,0) = dtprt(3,i)

            ! (PDIR is already in degrees nautical - convert back to

            !  Cartesian in radians to maintain internal convention)

            IF(dtprt(4,i) .NE. undef) THEN

              pdir(jsea,0) = (270. - dtprt(4,i)) * dera

            ENDIF

            psi(jsea,0) = dtprt(5,i)

            pws(jsea,0) = dtprt(6,i)

            ! (PTHP0 is already in degrees nautical - convert back to

            !  Cartesian in radians to maintain internal convention)

            IF(dtprt(7,i) .NE. undef) THEN

              pthp0(jsea,0) = (270. - dtprt(7,i)) * dera

            ENDIF

            psw(jsea,0) = dtprt(8,i)

            ppe(jsea,0) = dtprt(9,i)

            pqp(jsea,0) = dtprt(10,i)

            pgw(jsea,0) = dtprt(11,i)

            ptm1(jsea,0) = dtprt(12,i)

            pt1(jsea,0) = dtprt(13,i)

            pt2(jsea,0) = dtprt(14,i)

            pep(jsea,0) = dtprt(15,i)

            i      = i + 1

          END IF

          DO j=1, noswll

            IF ( i .GT.  icprt(jsea,2)+icprt(jsea,1)-1 ) EXIT

            phs(jsea,j) = dtprt(1,i)

            ptp(jsea,j) = dtprt(2,i)

            plp(jsea,j) = dtprt(3,i)

            ! (PDIR is already in degrees nautical - convert back to

            !  Cartesian in radians to maintain internal convention)

            IF(dtprt(4,i) .NE. undef) THEN

              pdir(jsea,j) = (270. - dtprt(4,i)) * dera

            ENDIF

            psi(jsea,j) = dtprt(5,i)

            pws(jsea,j) = dtprt(6,i)

            ! (PTHP0 is already in degrees nautical - convert back to

            !  Cartesian in radians to maintain internal convention)

            IF(dtprt(7,i) .NE. undef) THEN

              pthp0(jsea,j) = (270. - dtprt(7,i)) * dera

            ENDIF

            psw(jsea,j) = dtprt(8,i)

            ppe(jsea,j) = dtprt(9,i)

            pqp(jsea,j) = dtprt(10,i)

            pgw(jsea,j) = dtprt(11,i)

            ptm1(jsea,j) = dtprt(12,i)

            pt1(jsea,j) = dtprt(13,i)

            pt2(jsea,j) = dtprt(14,i)

            pep(jsea,j) = dtprt(15,i)

            i      = i + 1

          END DO

        END IF

        !

      END DO

      !

#ifdef W3_OMPG

      !$OMP END PARALLEL DO

#endif

      !


    END IF


    IF (floloc( 6, 8)) THEN

      CALL calc_u3stokes(a,1)

    END IF

    !

    IF (floloc( 6, 12)) THEN

      CALL calc_u3stokes(a,2)

    ENDIF

    !

    IF (floloc( 8, 7).OR.floloc( 8, 8).OR.floloc( 8, 9)) THEN

      CALL skewness(a)

    END IF


    !

    ! Dominant wave breaking probability

    !

    IF (floloc(2, 17)) CALL calc_wbt(a)

    !

    RETURN

    !

    ! Formats

    !

#ifdef W3_T

9050 FORMAT (' TEST W3OUTG : ISEA, IX, IY, HS, L, Tm, THm, THs',     &

         ', FP0, THP0')

9051 FORMAT (2x,i8,2i8)

9052 FORMAT (2x,i8,2i8,f6.2)

9053 FORMAT (2x,i8,2i8,f6.2,f7.1,f6.2,2f6.1)

9054 FORMAT (2x,i8,2i8,f6.2,f7.1,f6.2,2f6.1,f6.3,f6.0)

#endif


    !/

    !/ End of W3OUTG ----------------------------------------------------- /

    !/

  END SUBROUTINE w3outg

  !/ ------------------------------------------------------------------- /

  !/

  SUBROUTINE w3iogo ( INXOUT, NDSOG, IOTST, IMOD &

#ifdef W3_ASCII

                      ,NDSOA &

#endif

          )

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           H. L. Tolman            |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         22-Mar-2021 |

    !/                  +-----------------------------------+

    !/

    !/    17-Mar-1999 : Distributed FORTRAN 77 version.     ( version 1.18 )

    !/    04-Jan-2000 : Upgrade to FORTRAN 90               ( version 2.00 )

    !/                  Major changes to logistics.

    !/    24-Jan-2001 : Flat grid version (formats only)    ( version 2.06 )

    !/    23-Apr-2002 : Clean up                            ( version 2.19 )

    !/    29-Apr-2002 : Add output types 17-18.             ( version 2.20 )

    !/    13-Nov-2002 : Add stress vector.                  ( version 3.00 )

    !/    25-Oct-2004 : Multiple grid version.              ( version 3.06 )

    !/    27-Jun-2005 : Adding MAPST2.                      ( version 3.07 )

    !/    21-Jul-2005 : Adding output fields 19-21.         ( version 3.07 )

    !/    27-Jun-2006 : Adding file name preamble.          ( version 3.09 )

    !/    05-Jul-2006 : Consolidate stress arrays.          ( version 3.09 )

    !/    02-Apr-2007 : Adding partitioned output.          ( version 3.11 )

    !/                  Adding user slots for outputs.

    !/    30-Oct-2009 : Implement curvilinear grid type.    ( version 3.14 )

    !/                  (W. E. Rogers & T. J. Campbell, NRL)

    !/    31-Oct-2010 : Implement unstructured grids        ( version 3.14 )

    !/                  (A. Roland and F. Ardhuin)

    !/    05-Feb-2011 : Renumbering of output fields        ( version 3.14 )

    !/                  (F. Ardhuin)

    !/    25-Dec-2012 : New output structure and smaller    ( version 4.11 )

    !/                  memory footprint.

    !/    21-Aug-2013 : Added missing cos,sin for UBA, ABA  ( version 4.11 )

    !/    27-Nov-2013 : Management of coupling output       ( version 4.18 )

    !/    01-Mar-2018 : Removed RTD code (now used in post  ( version 6.02 )

    !/                  processing code)

    !/    25-Aug-2018 : Add WBT parameter                   ( version 6.06 )

    !/    22-Mar-2021 : Add extra coupling fields as output ( version 7.13 )

    !/    07-Mar-2024 : Add Skewness parameters             ( version 7.13 )

    !/

    !  1. Purpose :

    !

    !     Read/write gridded output.

    !

    !  2. Method :

    !

    !     Fields in file are determined by flags in FLOGRD in W3ODATMD.

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !       INXOUT  C*(*)  I   Test string for read/write, valid are:

    !                          'READ' and 'WRITE'.

    !       NDSOG   Int.   I   File unit number.

    !       IOTST   Int.   O   Test indictor for reading.

    !                           0 : Fields read.

    !                          -1 : Past end of file.

    !       IMOD    Int.   I   Model number for W3GDAT etc.

    !     ----------------------------------------------------------------

    !

    !  4. Subroutines used :

    !

    !       See module documentation above.

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.

    !      WW3_OUTF  Prog.   N/A    Ouput postprocessor.

    !      WW3_GRIB  Prog.   N/A    Ouput postprocessor.

    !      GX_OUTF   Prog.   N/A    Ouput postprocessor.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !       Tests on INXOUT, file status and on array dimensions.

    !

    !  7. Remarks :

    !

    !     - MAPSTA is dumped as it contains information on the ice edge.

    !       Dynamic ice edges require MAPSTA to be dumped every time step.

    !     - The output file has the pre-defined name 'out_grd.FILEXT'.

    !     - The current components CX and CY are written to out_grd as

    !       components, but converted to magnitude and direction in most

    !       gridded and point output post-processors (except gx_outf).

    !     - All written direction are in degrees, nautical convention,

    !       but in reading, all is convered back to radians and cartesian

    !       conventions.

    !     - Before writing, wind and current directions are converted,

    !       wave directions are already in correct convention (see W3OUTG).

    !     - In MPP version of model data is supposed to be gatherd at the

    !       correct processor before the routine is called.

    !     - In MPP version routine is called by only one process, therefore

    !       no test on process for error messages is needed.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/ST1   First source term package (WAM3).

    !     !/ST2   Second source term package (TC96).

    !     !/S     Enable subroutine tracing.

    !     !/T     Test output.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

    USE constants

    USE w3gdatmd

    !/

    USE w3wdatmd, ONLY: w3setw, w3dimw

    USE w3adatmd, ONLY: w3seta, w3dima, w3xeta

    USE w3odatmd, ONLY: w3seto

    !/

    USE w3wdatmd, ONLY: time, dinit, wlv, ice, icef, iceh, berg,    &

         ust,  ustdir, asf, rhoair

    USE w3adatmd, ONLY: ainit, dw, ua, ud, as, cx, cy, wn,          &

         taua, tauadir

    USE w3adatmd, ONLY: hs, wlm, t02, t0m1, t01, fp0, thm, ths, thp0,&

         wbt, wnmean

    USE w3adatmd, ONLY: dtdyn, fcut, aba, abd, uba, ubd, sxx, syy, sxy,&

         phs, ptp, plp, pdir, psi, pws, pwst, pnr,    &

         pthp0, pqp, ppe, pgw, psw, ptm1, pt1, pt2,  &

         pep, usero, tauox, tauoy, tauwix, tauwiy,    &

         phiaw, phioc, tusx, tusy, prms, tpms,        &

         ussx, ussy, mssx, mssy, mssd, mscx, mscy,    &

         mscd, qp, tauwnx, tauwny, charn, tws, bhd,   &

         phibbl, taubbl, whitecap, bedforms, cge, ef, &

         cflxymax, cflthmax, cflkmax, p2sms, us3d,    &

         th1m, sth1m, th2m, sth2m, hsig, phice, tauice,&

         stmaxe, stmaxd, hmaxe, hcmaxe, hmaxd, hcmaxd,&

         ussp, tauocx, tauocy, qkk, skew, embia1, embia2

    !/

    USE w3odatmd, ONLY: nogrp, ngrpp, idout, undef, ndst, ndse,     &

         flogrd, ipass => ipass1, write => write1,   &

         fnmpre, noswll, noextr

    !/

    USE w3servmd, ONLY: extcde

    USE w3odatmd, only : iaproc

    USE w3odatmd, ONLY :  ofiles

#ifdef W3_SETUP

    USE w3wdatmd, ONLY: zeta_setup

#endif

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    INTEGER, INTENT(INOUT)        :: IOTST

    INTEGER, INTENT(IN)           :: NDSOG

    INTEGER, INTENT(IN), OPTIONAL :: IMOD

    CHARACTER, INTENT(IN)         :: INXOUT*(*)

    CHARACTER(LEN=15) :: TIMETAG

#ifdef W3_ASCII

    INTEGER, INTENT(IN), OPTIONAL :: NDSOA

#endif

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

    INTEGER                 :: IGRD, IERR, I, J, IX, IY, MOGRP,     &

         MGRPP, ISEA, MOSWLL, IK, IFI, IFJ    &

         ,IFILOUT

    INTEGER, ALLOCATABLE    :: MAPTMP(:,:)

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    REAL                    :: AUX1(NSEA), AUX2(NSEA),              &

         AUX3(NSEA), AUX4(NSEA)

#ifdef W3_SMC

    REAL                    :: UDARC

#endif

    CHARACTER(LEN=30)       :: IDTST, TNAME

    CHARACTER(LEN=10)       :: VERTST

    !/

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'W3IOGO')

#endif

    !

    ! 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 w3seta ( igrd, ndse, ndst )

#ifdef W3_MPI

    CALL w3xeta ( igrd, ndse, ndst )

#endif

    CALL w3setw ( igrd, ndse, ndst )

    !

    ipass  = ipass + 1

    iotst  = 0

    !

    IF (inxout.NE.'READ' .AND. inxout.NE.'WRITE' ) THEN

      WRITE (ndse,900) inxout

      CALL extcde ( 1 )

    END IF

    !

    IF ( ipass.EQ.1 .AND. ofiles(1) .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

    !

#ifdef W3_T

    WRITE (ndst,9000) ipass, inxout, WRITE, ndsog, igrd, filext

#endif

    !

    !

    ! open file ---------------------------------------------------------- *

    ! ( IPASS = 1 )

    !

    IF ( ipass.EQ.1 .AND. ofiles(1) .EQ. 0) THEN

      i      = len_trim(filext)

      j      = len_trim(fnmpre)

      !

#ifdef W3_T

      WRITE (ndst,9001) fnmpre(:j)//'out_grd.'//filext(:i)

#endif

      IF ( WRITE ) THEN

        OPEN (ndsog,file=fnmpre(:j)//'out_grd.'//filext(:i),    &

             form ='UNFORMATTED', convert=file_endian,err=800,iostat=ierr)

#ifdef W3_ASCII

        OPEN (ndsoa,file=fnmpre(:j)//'out_grd.'//filext(:i)//'.txt',    &

             form ='FORMATTED',err=800,iostat=ierr)

#endif

      ELSE

        OPEN (ndsog,file=fnmpre(:j)//'out_grd.'//filext(:i),    &

             form='UNFORMATTED', convert=file_endian,err=800,iostat=ierr,status='OLD')

      END IF

      !

      rewind( ndsog )

      !

      ! test info --------------------------------------------------------- *

      ! ( IPASS = 1 )

      !

      IF ( WRITE ) THEN

        WRITE (ndsog)                                           &

             idstr, verogr, gname, nogrp, ngrpp, nsea, nx, ny,     &

             undef, noswll

#ifdef W3_ASCII

        WRITE (ndsoa,*)                                           &

             'IDSTR, VEROGR, GNAME, NOGRP, NGRPP, NSEA, NX, NY,     &

             UNDEF, NOSWLL:',                                     &

             idstr, verogr, gname, nogrp, ngrpp, nsea, nx, ny,     &

             undef, noswll

#endif

      ELSE

        READ (ndsog,END=801,ERR=802,IOSTAT=IERR)                &

             idtst, vertst, tname, mogrp, mgrpp, nsea, nx, ny,     &

             undef, moswll

        !

        IF ( idtst .NE. idstr ) THEN

          WRITE (ndse,902) idtst, idstr

          CALL extcde ( 20 )

        END IF

        IF ( vertst .NE. verogr ) THEN

          WRITE (ndse,903) vertst, verogr

          CALL extcde ( 21 )

        END IF

        IF ( nogrp .NE. mogrp .OR. ngrpp .NE. mgrpp ) THEN

          WRITE (ndse,904) mogrp, mgrpp, nogrp, ngrpp

          CALL extcde ( 22 )

        END IF

        IF ( tname .NE. gname ) THEN

          WRITE (ndse,905) tname, gname

        END IF

        IF ( noswll .NE. moswll ) THEN

          WRITE (ndse,906) moswll, noswll

          CALL extcde ( 24 )

        END IF

        !

      END IF

      !

#ifdef W3_T

      WRITE (ndst,9002) idstr, verogr, gname, nsea, nx, ny,    &

           undef

#endif

      !

    END IF

    !

    !  IN CASE OF GENERATION OF A NEW FILE OUTPUT EVERY DELTA OUTPUT

    ! open file ---------------------------------------------------------- *

    ! ( IPASS = 1 )

    !

    IF ( ipass.GE.1 .AND. ofiles(1) .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


    !

    IF ( ipass.GE.1 .AND. ofiles(1) .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_grd.'//filext(:i)

#endif

      IF ( WRITE ) THEN

        OPEN (ndsog,file=fnmpre(:j)//timetag//'.out_grd.'  &

             //filext(:i),form='UNFORMATTED', convert=file_endian,err=800,iostat=ierr)

#ifdef W3_ASCII

        OPEN (ndsoa,file=fnmpre(:j)//timetag//'.out_grd.'  &

             //filext(:i)//'.txt',form='FORMATTED',err=800,iostat=ierr)

#endif

      ELSE

        OPEN (ndsog,file=fnmpre(:j)//'out_grd.'//filext(:i),    &

             form='UNFORMATTED', convert=file_endian,err=800,iostat=ierr,status='OLD')

      END IF

      !

      rewind( ndsog )

      !

      ! test info --------------------------------------------------------- *

      ! ( IPASS >= 1 & OFILES(1) = 1)

      !

      IF ( WRITE ) THEN

        WRITE (ndsog)                                           &

             idstr, verogr, gname, nogrp, ngrpp, nsea, nx, ny,     &

             undef, noswll

#ifdef W3_ASCII

        WRITE (ndsoa,*)                                           &

             'IDSTR, VEROGR, GNAME, NOGRP, NGRPP, NSEA, NX, NY,     &

             UNDEF, NOSWLL:',                                     &

             idstr, verogr, gname, nogrp, ngrpp, nsea, nx, ny,     &

             undef, noswll

#endif

      ELSE

        READ (ndsog,END=801,ERR=802,IOSTAT=IERR)                &

             idtst, vertst, tname, mogrp, mgrpp, nsea, nx, ny,     &

             undef, moswll

        !

        IF ( idtst .NE. idstr ) THEN

          WRITE (ndse,902) idtst, idstr

          CALL extcde ( 20 )

        END IF

        IF ( vertst .NE. verogr ) THEN

          WRITE (ndse,903) vertst, verogr

          CALL extcde ( 21 )

        END IF

        IF ( nogrp .NE. mogrp .OR. ngrpp .NE. mgrpp ) THEN

          WRITE (ndse,904) mogrp, mgrpp, nogrp, ngrpp

          CALL extcde ( 22 )

        END IF

        IF ( tname .NE. gname ) THEN

          WRITE (ndse,905) tname, gname

        END IF

        IF ( noswll .NE. moswll ) THEN

          WRITE (ndse,906) moswll, noswll

          CALL extcde ( 24 )

        END IF

        !

      END IF

      !

#ifdef W3_T

      WRITE (ndst,9002) idstr, verogr, gname, nsea, nx, ny,    &

           undef

#endif

      !

    END IF

    !

    ! TIME and flags ----------------------------------------------------- *

    !

    IF ( WRITE ) THEN

      WRITE (ndsog)                            time, flogrd

#ifdef W3_ASCII

      WRITE (ndsoa,*)                         'TIME, FLOGRD:', &

                                               time, flogrd

#endif

    ELSE

      READ (ndsog,END=803,ERR=802,IOSTAT=IERR) TIME, flogrd

    END IF

    !

#ifdef W3_T

    WRITE (ndst,9003) time, flogrd

#endif

    !

    ! MAPSTA ------------------------------------------------------------- *

    !

    ALLOCATE ( maptmp(ny,nx) )

    IF ( WRITE ) THEN

      maptmp = mapsta + 8*mapst2

      WRITE (ndsog)                                               &

           ((maptmp(iy,ix),ix=1,nx),iy=1,ny)

#ifdef W3_ASCII

      WRITE (ndsoa,*) 'MAPSTA:',                                  &

           ((maptmp(iy,ix),ix=1,nx),iy=1,ny)

#endif

    ELSE

      READ (ndsog,END=801,ERR=802,IOSTAT=IERR)                    &

           ((maptmp(iy,ix),ix=1,nx),iy=1,ny)

      mapsta = mod(maptmp+2,8) - 2

      mapst2 = (maptmp-mapsta) / 8

    END IF

    DEALLOCATE ( maptmp )

    !

    ! Fields ---------------------------------------------- *

    !

    ! Initialization ---------------------------------------------- *

    !

    IF ( WRITE ) THEN

      DO isea=1, nsea

        IF ( mapsta(mapsf(isea,2),mapsf(isea,1)) .LT. 0 ) THEN

          !

          IF ( flogrd( 2, 2) ) wlm(isea) = undef

          IF ( flogrd( 2, 3) ) t02(isea) = undef

          IF ( flogrd( 2, 4) ) t0m1(isea) = undef

          IF ( flogrd( 2, 5) ) t01(isea) = undef

          IF ( flogrd( 2, 6) .OR. flogrd( 2,18) )                 &

               fp0(isea) = undef  ! FP or TP

          IF ( flogrd( 2, 7) ) thm(isea) = undef

          IF ( flogrd( 2, 8) ) ths(isea) = undef

          IF ( flogrd( 2, 9) ) thp0(isea) = undef

          ust(isea) = undef

          ustdir(isea) = undef

          IF ( flogrd( 2,10) ) hsig(isea) = undef

          IF ( flogrd( 2,11) ) stmaxe(isea) = undef

          IF ( flogrd( 2,12) ) stmaxd(isea) = undef

          IF ( flogrd( 2,13) ) hmaxe(isea) = undef

          IF ( flogrd( 2,14) ) hcmaxe(isea) = undef

          IF ( flogrd( 2,15) ) hmaxd(isea) = undef

          IF ( flogrd( 2,16) ) hcmaxd(isea) = undef

          IF ( flogrd( 2,17) ) wbt(isea) = undef

          IF ( flogrd( 2,19) ) wnmean(isea) = undef

          !

          IF ( flogrd( 3, 1) ) ef(isea,:) = undef

          IF ( flogrd( 3, 2) ) th1m(isea,:) = undef

          IF ( flogrd( 3, 3) ) sth1m(isea,:) = undef

          IF ( flogrd( 3, 4) ) th2m(isea,:) = undef

          IF ( flogrd( 3, 5) ) sth2m(isea,:) = undef

          !

          IF ( flogrd( 4, 1) ) phs(isea,:) = undef

          IF ( flogrd( 4, 2) ) ptp(isea,:) = undef

          IF ( flogrd( 4, 3) ) plp(isea,:) = undef

          IF ( flogrd( 4, 4) ) pdir(isea,:) = undef

          IF ( flogrd( 4, 5) ) psi(isea,:) = undef

          IF ( flogrd( 4, 6) ) pws(isea,:) = undef

          IF ( flogrd( 4, 7) ) pthp0(isea,:) = undef

          IF ( flogrd( 4, 8) ) pqp(isea,:) = undef

          IF ( flogrd( 4, 9) ) ppe(isea,:)  = undef

          IF ( flogrd( 4,10) ) pgw(isea,:)  = undef

          IF ( flogrd( 4,11) ) psw(isea,:) = undef

          IF ( flogrd( 4,12) ) ptm1(isea,:) = undef

          IF ( flogrd( 4,13) ) pt1(isea,:) = undef

          IF ( flogrd( 4,14) ) pt2(isea,:) = undef

          IF ( flogrd( 4,15) ) pep(isea,:) = undef

          IF ( flogrd( 4,16) ) pwst(isea  ) = undef

          IF ( flogrd( 4,17) ) pnr(isea  ) = undef

          !

          IF ( flogrd( 5, 2) ) charn(isea) = undef

          IF ( flogrd( 5, 3) ) cge(isea) = undef

          IF ( flogrd( 5, 4) ) phiaw(isea) = undef

          IF ( flogrd( 5, 5) ) THEN

            tauwix(isea) = undef

            tauwiy(isea) = undef

          END IF

          IF ( flogrd( 5, 6) ) THEN

            tauwnx(isea) = undef

            tauwny(isea) = undef

          END IF

          IF ( flogrd( 5, 7) ) whitecap(isea,1) = undef

          IF ( flogrd( 5, 8) ) whitecap(isea,2) = undef

          IF ( flogrd( 5, 9) ) whitecap(isea,3) = undef

          IF ( flogrd( 5,10) ) whitecap(isea,4) = undef

          !

          IF ( flogrd( 6, 1) ) THEN

            sxx(isea) = undef

            syy(isea) = undef

            sxy(isea) = undef

          END IF

          IF ( flogrd( 6, 2) ) THEN

            tauox(isea) = undef

            tauoy(isea) = undef

          END IF

          IF ( flogrd( 6, 3) ) bhd(isea) = undef

          IF ( flogrd( 6, 4) ) phioc(isea) = undef

          IF ( flogrd( 6, 5) ) THEN

            tusx(isea) = undef

            tusy(isea) = undef

          END IF

          IF ( flogrd( 6, 6) ) THEN

            ussx(isea) = undef

            ussy(isea) = undef

          END IF

          IF ( flogrd( 6, 7) ) THEN

            prms(isea) = undef

            tpms(isea) = undef

          END IF

          IF ( flogrd( 6, 8) ) us3d(isea,:) = undef

          IF ( flogrd( 6, 9) ) p2sms(isea,:) = undef

          IF ( flogrd( 6, 10) ) tauice(isea,:) = undef

          IF ( flogrd( 6, 11) ) phice(isea) = undef

          IF ( flogrd( 6, 12) ) ussp(isea,:) = undef

          IF ( flogrd( 6, 13) ) THEN

            tauocx(isea) = undef

            tauocy(isea) = undef

          END IF

          !

          IF ( flogrd( 7, 1) ) THEN

            aba(isea) = undef

            abd(isea) = undef

          END IF

          IF ( flogrd( 7, 2) ) THEN

            uba(isea) = undef

            ubd(isea) = undef

          END IF

          IF ( flogrd( 7, 3) ) bedforms(isea,:) = undef

          IF ( flogrd( 7, 4) ) phibbl(isea) = undef

          IF ( flogrd( 7, 5) ) taubbl(isea,:) = undef

          !

          IF ( flogrd( 8, 1) ) THEN

            mssx(isea) = undef

            mssy(isea) = undef

          END IF

          IF ( flogrd( 8, 2) ) THEN

            mscx(isea) = undef

            mscy(isea) = undef

          END IF

          IF ( flogrd( 8, 3) ) mssd(isea) = undef

          IF ( flogrd( 8, 4) ) mscd(isea) = undef

          IF ( flogrd( 8, 5) ) qp(isea) = undef

          IF ( flogrd( 8, 6) ) qkk(isea) = undef

          IF ( flogrd( 8, 7) ) skew(isea) = undef

          IF ( flogrd( 8, 8) ) embia1(isea) = undef

          IF ( flogrd( 8, 9) ) embia2(isea) = undef

          !

          IF ( flogrd( 9, 1) ) dtdyn(isea) = undef

          IF ( flogrd( 9, 2) ) fcut(isea) = undef

          IF ( flogrd( 9, 3) ) cflxymax(isea) = undef

          IF ( flogrd( 9, 4) ) cflthmax(isea) = undef

          IF ( flogrd( 9, 5) ) cflkmax(isea) = undef

          !

        END IF

        !

        IF ( mapsta(mapsf(isea,2),mapsf(isea,1)) .EQ. 2 ) THEN

          !

          IF ( flogrd( 5, 4) ) phiaw(isea) = undef

          IF ( flogrd( 5, 5) ) THEN

            tauwix(isea) = undef

            tauwiy(isea) = undef

          END IF

          IF ( flogrd( 5, 6) ) THEN

            tauwnx(isea) = undef

            tauwny(isea) = undef

          END IF

          IF ( flogrd( 5, 7) ) whitecap(isea,1) = undef

          IF ( flogrd( 5, 8) ) whitecap(isea,2) = undef

          IF ( flogrd( 5, 9) ) whitecap(isea,3) = undef

          IF ( flogrd( 5,10) ) whitecap(isea,4) = undef

          !

          IF ( flogrd( 6, 2) )THEN

            tauox(isea) = undef

            tauoy(isea) = undef

          END IF

          IF ( flogrd( 6, 4) ) phioc(isea) = undef

          !

          IF ( flogrd( 7, 3) ) bedforms(isea,:) = undef

          IF ( flogrd( 7, 4) ) phibbl(isea) = undef

          IF ( flogrd( 7, 5) ) taubbl(isea,:) = undef

          !

        END IF

        !

      END DO

      !

    ELSE

      IF (.NOT.dinit) CALL w3dimw ( igrd, ndse, ndst, .true. )

      IF (.NOT.ainit) CALL w3dima ( igrd, ndse, ndst, .true. )

    END IF

    !

    ! Actual output  ---------------------------------------------- *

    DO ifi=1, nogrp

      DO ifj=1, ngrpp


        IF ( flogrd(ifi,ifj) ) THEN

          !

#ifdef W3_T

          WRITE (ndst,9010) flogrd(ifi,ifj), idout(ifi,ifj)

#endif

          !

          IF ( WRITE ) THEN

            !

            !     Section 1)

            !

            IF ( ifi .EQ. 1 .AND. ifj .EQ. 1 ) THEN

              WRITE ( ndsog ) dw(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'DW:', dw(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 2 ) THEN

              WRITE ( ndsog ) cx(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'CX:', cx(1:nsea)

#endif

              WRITE ( ndsog ) cy(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'CY:', cy(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 3 ) THEN

              DO isea=1, nsea

#ifdef W3_SMC

                !!Li  Rotate map-east wind in Arctic part back to local east.  JGLi02Feb2016

                IF( arctc .AND. (isea .GT. nglo) ) THEN

                  udarc = ud(isea) - angarc(isea - nglo)*dera

                  ud(isea) = mod(tpi + udarc, tpi)

                ENDIF

#endif

                IF (ua(isea) .NE.undef) THEN

                  aux1(isea) = ua(isea)*cos(ud(isea))

                  aux2(isea) = ua(isea)*sin(ud(isea))

                ELSE

                  aux1(isea) = undef

                  aux2(isea) = undef

                END IF

              END DO

              WRITE ( ndsog ) aux1

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'AUX1 (UA*cos(UD)):', aux1

#endif

              WRITE ( ndsog ) aux2

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'AUX2 (UA*sin(UD)):', aux2

#endif

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 4 ) THEN

              WRITE ( ndsog ) as(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'AS:', as(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 5 ) THEN

              WRITE ( ndsog ) wlv(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'WLV:', wlv(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 6 ) THEN

              WRITE ( ndsog ) ice(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'ICE:', ice(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 7 ) THEN

              WRITE ( ndsog ) berg(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'BERG:', berg(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 8 ) THEN

              DO isea=1, nsea

#ifdef W3_SMC

                !!Li  Rotate map-east momentum in Arctic part back to local east.  JGLi02Feb2016

                IF( arctc .AND. (isea .GT. nglo) ) THEN

                  udarc = tauadir(isea) - angarc(isea - nglo)*dera

                  tauadir(isea) = mod(tpi + udarc, tpi)

                ENDIF

#endif

                IF (taua(isea) .NE.undef) THEN

                  aux1(isea) = taua(isea)*cos(tauadir(isea))

                  aux2(isea) = taua(isea)*sin(tauadir(isea))

                ELSE

                  aux1(isea) = undef

                  aux2(isea) = undef

                END IF

              END DO

              WRITE ( ndsog ) aux1

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'AUX1 (TAUA*cos(TAUADIR)):', aux1

#endif

              WRITE ( ndsog ) aux2

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'AUX2 (TAUA*sin(TAUADIR)):', aux2

#endif

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 9 ) THEN

              WRITE ( ndsog ) rhoair(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'RHOAIR:', rhoair(1:nsea)

#endif

#ifdef W3_BT4

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 10 ) THEN

              WRITE ( ndsog ) sed_d50(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'SED_D50:', sed_d50(1:nsea)

#endif

#endif

#ifdef W3_IS2

            ELSE IF (ifi .EQ. 1 .AND. ifj .EQ. 11 ) THEN

              WRITE (ndsog ) iceh(1:nsea)

#ifdef W3_ASCII

              WRITE (ndsoa,* ) 'ICEH:', iceh(1:nsea)

#endif

            ELSE IF (ifi .EQ. 1 .AND. ifj .EQ. 12 ) THEN

              WRITE (ndsog ) icef(1:nsea)

#ifdef W3_ASCII

              WRITE (ndsoa,* ) 'ICEF:', icef(1:nsea)

#endif

#endif

#ifdef W3_SETUP

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 13 ) THEN

              WRITE ( ndsog ) zeta_setup(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'ZETA_SETUP:', zeta_setup(1:nsea)

#endif

#endif


              !

              !     Section 2)

              !

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 1 ) THEN

              WRITE ( ndsog ) hs(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'HS:', hs(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 2 ) THEN

              WRITE ( ndsog ) wlm(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'WLM:', wlm(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 3 ) THEN

              WRITE ( ndsog ) t02(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'T02:', t02(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 4 ) THEN

              WRITE ( ndsog ) t0m1(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'T0M1:', t0m1(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 5 ) THEN

              WRITE ( ndsog ) t01(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'T01:', t01(1:nsea)

#endif

            ELSE IF ( (ifi .EQ. 2 .AND. ifj .EQ. 6) .OR.         &

                 (ifi .EQ. 2 .AND. ifj .EQ. 18) ) THEN

              ! Note: TP output is derived from FP field.

              WRITE ( ndsog ) fp0(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'FP0:', fp0(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 7 ) THEN

              WRITE ( ndsog ) thm(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'THM:', thm(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 8 ) THEN

              WRITE ( ndsog ) ths(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'THS:', ths(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 9 ) THEN

              WRITE ( ndsog ) thp0(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'THP0:', thp0(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 10 ) THEN

              WRITE ( ndsog ) hsig(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'HSIG:', hsig(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 11 ) THEN

              WRITE ( ndsog ) stmaxe(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'STMAXE:', stmaxe(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 12 ) THEN

              WRITE ( ndsog ) stmaxd(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'STMAXD:', stmaxd(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 13 ) THEN

              WRITE ( ndsog ) hmaxe(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'HMAXE:', hmaxe(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 14 ) THEN

              WRITE ( ndsog ) hcmaxe(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'HCMAXE:', hcmaxe(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 15 ) THEN

              WRITE ( ndsog ) hmaxd(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'HMAXD:', hmaxd(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 16 ) THEN

              WRITE ( ndsog ) hcmaxd(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'HCMAXD:', hcmaxd(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 17 ) THEN

              WRITE ( ndsog ) wbt(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'WBT:', wbt(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 19 ) THEN

              WRITE ( ndsog ) wnmean(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'WNMEAN:', wnmean(1:nsea)

#endif

              !

              !     Section 3)

              !

            ELSE IF ( ifi .EQ. 3 .AND. ifj .EQ. 1 ) THEN

              WRITE ( ndsog ) ef(1:nsea,e3df(2,1):e3df(3,1))

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'EF:', ef(1:nsea,e3df(2,1):e3df(3,1))

#endif

            ELSE IF ( ifi .EQ. 3 .AND. ifj .EQ. 2 ) THEN

              WRITE ( ndsog ) th1m(1:nsea,e3df(2,2):e3df(3,2))

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TH1M:', th1m(1:nsea,e3df(2,2):e3df(3,2))

#endif

            ELSE IF ( ifi .EQ. 3 .AND. ifj .EQ. 3 ) THEN

              WRITE ( ndsog ) sth1m(1:nsea,e3df(2,3):e3df(3,3))

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'STH1M:', sth1m(1:nsea,e3df(2,3):e3df(3,3))

#endif

            ELSE IF ( ifi .EQ. 3 .AND. ifj .EQ. 4 ) THEN

              WRITE ( ndsog ) th2m(1:nsea,e3df(2,4):e3df(3,4))

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TH2M:', th2m(1:nsea,e3df(2,4):e3df(3,4))

#endif

            ELSE IF ( ifi .EQ. 3 .AND. ifj .EQ. 5 ) THEN

              WRITE ( ndsog ) sth2m(1:nsea,e3df(2,5):e3df(3,5))

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'STH2M:', sth2m(1:nsea,e3df(2,5):e3df(3,5))

#endif

            ELSE IF ( ifi .EQ. 3 .AND. ifj .EQ. 6) THEN

              WRITE ( ndsog ) wn(1:nk,1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'WN:', wn(1:nk,1:nsea)

#endif

              !

              !     Section 4)

              !

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 1 ) THEN

              WRITE ( ndsog ) phs(1:nsea,0:noswll)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PHS:', phs(1:nsea,0:noswll)

#endif

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 2 ) THEN

              WRITE ( ndsog ) ptp(1:nsea,0:noswll)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PTP:', ptp(1:nsea,0:noswll)

#endif

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 3 ) THEN

              WRITE ( ndsog ) plp(1:nsea,0:noswll)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PLP:', plp(1:nsea,0:noswll)

#endif

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 4 ) THEN

              WRITE ( ndsog ) pdir(1:nsea,0:noswll)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PDIR:', pdir(1:nsea,0:noswll)

#endif

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 5 ) THEN

              WRITE ( ndsog ) psi(1:nsea,0:noswll)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PSI:', psi(1:nsea,0:noswll)

#endif

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 6 ) THEN

              WRITE ( ndsog ) pws(1:nsea,0:noswll)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PWS:', pws(1:nsea,0:noswll)

#endif

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 7 ) THEN

              WRITE ( ndsog ) pthp0(1:nsea,0:noswll)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PTHP0:', pthp0(1:nsea,0:noswll)

#endif

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 8  ) THEN

              WRITE ( ndsog ) pqp(1:nsea,0:noswll)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PQP:', pqp(1:nsea,0:noswll)

#endif

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 9  ) THEN

              WRITE ( ndsog ) ppe(1:nsea,0:noswll)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PPE:', ppe(1:nsea,0:noswll)

#endif

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 10 ) THEN

              WRITE ( ndsog ) pgw(1:nsea,0:noswll)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PGW:', pgw(1:nsea,0:noswll)

#endif

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 11 ) THEN

              WRITE ( ndsog ) psw(1:nsea,0:noswll)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PSW:', psw(1:nsea,0:noswll)

#endif

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 12 ) THEN

              WRITE ( ndsog ) ptm1(1:nsea,0:noswll)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PTM1:', ptm1(1:nsea,0:noswll)

#endif

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 13 ) THEN

              WRITE ( ndsog ) pt1(1:nsea,0:noswll)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PT1:', pt1(1:nsea,0:noswll)

#endif

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 14 ) THEN

              WRITE ( ndsog ) pt2(1:nsea,0:noswll)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PT2:', pt2(1:nsea,0:noswll)

#endif

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 15 ) THEN

              WRITE ( ndsog ) pep(1:nsea,0:noswll)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PEP:', pep(1:nsea,0:noswll)

#endif

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 16 ) THEN

              WRITE ( ndsog ) pwst(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PWST:', pwst(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 17 ) THEN

              WRITE ( ndsog ) pnr(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PNR:', pnr(1:nsea)

#endif

              !

              !     Section 5)

              !

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 1 ) THEN

              DO isea=1, nsea

                ix     = mapsf(isea,1)

                iy     = mapsf(isea,2)

                IF ( mapsta(iy,ix) .EQ. 1 ) THEN

                  aux1(isea) = ust(isea) * asf(isea) *        &

                       cos(ustdir(isea))

                  aux2(isea) = ust(isea) * asf(isea) *        &

                       sin(ustdir(isea))

                ELSE

                  aux1(isea) = undef

                  aux2(isea) = undef

                END IF

              END DO

              WRITE ( ndsog ) aux1

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'AUX1 (UST*ASF*cos(USTDIR)):', aux1

#endif

              WRITE ( ndsog ) aux2

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'AUX2 (UST*ASF*sin(USTDIR)):', aux2

#endif

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 2 ) THEN

              WRITE ( ndsog ) charn(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'CHARN:', charn(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 3 ) THEN

              WRITE ( ndsog ) cge(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'CGE:', cge(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 4 ) THEN

              WRITE ( ndsog ) phiaw(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PHIAW:', phiaw(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 5 ) THEN

              WRITE ( ndsog ) tauwix(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TAUWIX:', tauwix(1:nsea)

#endif

              WRITE ( ndsog ) tauwiy(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TAUWIY:', tauwiy(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 6 ) THEN

              WRITE ( ndsog ) tauwnx(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TAUWNX:', tauwnx(1:nsea)

#endif

              WRITE ( ndsog ) tauwny(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TAUWNY:', tauwny(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 7 ) THEN

              WRITE ( ndsog ) whitecap(1:nsea,1)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'WHITECAP(1):', whitecap(1:nsea,1)

#endif

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 8 ) THEN

              WRITE ( ndsog ) whitecap(1:nsea,2)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'WHITECAP(2):', whitecap(1:nsea,2)

#endif

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 9 ) THEN

              WRITE ( ndsog ) whitecap(1:nsea,3)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'WHITECAP(3):', whitecap(1:nsea,3)

#endif

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 10 ) THEN

              WRITE ( ndsog ) whitecap(1:nsea,4)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'WHITECAP(4):', whitecap(1:nsea,4)

#endif

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 11 ) THEN

              WRITE ( ndsog ) tws(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TWS:', tws(1:nsea)

#endif

              !

              !     Section 6)

              !

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 1 ) THEN

              WRITE ( ndsog ) sxx(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'SXX:', sxx(1:nsea)

#endif

              WRITE ( ndsog ) syy(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'SYY:', syy(1:nsea)

#endif

              WRITE ( ndsog ) sxy(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'SXY:', sxy(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 2 ) THEN

              WRITE ( ndsog ) tauox(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TAUOX:', tauox(1:nsea)

#endif

              WRITE ( ndsog ) tauoy(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TAUOY:', tauoy(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 3  ) THEN

              WRITE ( ndsog ) bhd(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'BHD:', bhd(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 4 ) THEN

              WRITE ( ndsog ) phioc(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PHIOC:', phioc(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 5 ) THEN

              WRITE ( ndsog ) tusx(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TUSX:', tusx(1:nsea)

#endif

              WRITE ( ndsog ) tusy(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TUSY:', tusy(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 6 ) THEN

              WRITE ( ndsog ) ussx(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'USSX:', ussx(1:nsea)

#endif

              WRITE ( ndsog ) ussy(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'USSY:', ussy(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 7 ) THEN

              WRITE ( ndsog ) prms(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PRMS:', prms(1:nsea)

#endif

              WRITE ( ndsog ) tpms(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TPMS:', tpms(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 8 ) THEN

              WRITE ( ndsog ) us3d(1:nsea,   us3df(2):us3df(3))

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'US3D:', us3d(1:nsea,   us3df(2):us3df(3))

#endif

              WRITE ( ndsog ) us3d(1:nsea,nk+us3df(2):nk+us3df(3))

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'US3D+NK:', us3d(1:nsea,nk+us3df(2):nk+us3df(3))

#endif

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ.  9 ) THEN

              WRITE ( ndsog ) p2sms(1:nsea,p2msf(2):p2msf(3))

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'P2SMS:', p2sms(1:nsea,p2msf(2):p2msf(3))

#endif

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 10 ) THEN

              WRITE ( ndsog ) tauice(1:nsea,1)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TAUICE(1):', tauice(1:nsea,1)

#endif

              WRITE ( ndsog ) tauice(1:nsea,2)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TAUICE(2):', tauice(1:nsea,2)

#endif

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 11 ) THEN

              WRITE ( ndsog ) phice(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PHICE:', phice(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 12 ) THEN

              WRITE ( ndsog ) ussp(1:nsea,   1:usspf(2))

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'USSP:', ussp(1:nsea,   1:usspf(2))

#endif

              WRITE ( ndsog ) ussp(1:nsea,nk+1:nk+usspf(2))

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'USSP:', ussp(1:nsea,nk+1:nk+usspf(2))

#endif

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 13 ) THEN

              WRITE ( ndsog ) tauocx(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TAUOCX:', tauocx(1:nsea)

#endif

              WRITE ( ndsog ) tauocy(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TAUOCY:', tauocy(1:nsea)

#endif

              !

              !     Section 7)

              !

            ELSE IF ( ifi .EQ. 7 .AND. ifj .EQ. 1 ) THEN

              DO isea=1, nsea

                IF ( aba(isea) .NE. undef ) THEN

                  aux1(isea) = aba(isea)*cos(abd(isea))

                  aux2(isea) = aba(isea)*sin(abd(isea))

                ELSE

                  aux1(isea) = undef

                  aux2(isea) = undef

                END IF

              END DO

              WRITE ( ndsog ) aux1

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'AUX1 (ABA*cos(ABD)):', aux1

#endif

              WRITE ( ndsog ) aux2

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'AUX2 (ABA*sin(ABD)):', aux2

#endif

              !WRITE ( NDSOG ) ABA(1:NSEA)

              !WRITE ( NDSOG ) ABD(1:NSEA)

            ELSE IF ( ifi .EQ. 7 .AND. ifj .EQ. 2 ) THEN

              DO isea=1, nsea

                IF ( uba(isea) .NE. undef ) THEN

                  aux1(isea) = uba(isea)*cos(ubd(isea))

                  aux2(isea) = uba(isea)*sin(ubd(isea))

                ELSE

                  aux1(isea) = undef

                  aux2(isea) = undef

                END IF

              END DO

              WRITE ( ndsog ) aux1

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'AUX1 (UBA*cos(UBD)):', aux1

#endif

              WRITE ( ndsog ) aux2

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'AUX2 (UBA*sin(UBD)):', aux2

#endif

              !                    WRITE ( NDSOG ) UBA(1:NSEA)

              !                    WRITE ( NDSOG ) UBD(1:NSEA)

            ELSE IF ( ifi .EQ. 7 .AND. ifj .EQ. 3 ) THEN

              WRITE ( ndsog ) bedforms(1:nsea,1)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'BEDFORMS(1):', bedforms(1:nsea,1)

#endif

              WRITE ( ndsog ) bedforms(1:nsea,2)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'BEDFORMS(2):', bedforms(1:nsea,2)

#endif

              WRITE ( ndsog ) bedforms(1:nsea,3)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'BEDFORMS(3):', bedforms(1:nsea,3)

#endif

            ELSE IF ( ifi .EQ. 7 .AND. ifj .EQ. 4 ) THEN

              WRITE ( ndsog ) phibbl(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'PHIBBL:', phibbl(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 7 .AND. ifj .EQ. 5 ) THEN

              WRITE ( ndsog ) taubbl(1:nsea,1)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TAUBBL(1):', taubbl(1:nsea,1)

#endif

              WRITE ( ndsog ) taubbl(1:nsea,2)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'TAUBBL(2):', taubbl(1:nsea,2)

#endif

              !

              !     Section 8)

              !Skewness

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 1 ) THEN

              WRITE ( ndsog ) mssx(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'MSSX:', mssx(1:nsea)

#endif

              WRITE ( ndsog ) mssy(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'MSSY:', mssy(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 2 ) THEN

              WRITE ( ndsog ) mscx(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'MSCX:', mscx(1:nsea)

#endif

              WRITE ( ndsog ) mscy(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'MSCY:', mscy(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 3 ) THEN

              WRITE ( ndsog ) mssd(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'MSSD:', mssd(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 4 ) THEN

              WRITE ( ndsog ) mscd(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'MSCD:', mscd(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 5 ) THEN

              WRITE ( ndsog ) qp(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'QP:', qp(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 6 ) THEN

              WRITE ( ndsog ) qkk(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'QKK:', qkk(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 7 ) THEN

              WRITE ( ndsog ) skew(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'SKW:', skew(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 8 ) THEN

              WRITE ( ndsog ) embia1(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'EMB:', embia1(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 9 ) THEN

              WRITE ( ndsog ) embia2(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'EMC:', embia2(1:nsea)

#endif

              !

              !     Section 9)

              !

            ELSE IF ( ifi .EQ. 9 .AND. ifj .EQ. 1 ) THEN

              WRITE ( ndsog ) dtdyn(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'DTDYN:', dtdyn(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 9 .AND. ifj .EQ. 2 ) THEN

              WRITE ( ndsog ) fcut(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'FCUT:', fcut(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 9 .AND. ifj .EQ. 3 ) THEN

              WRITE ( ndsog ) cflxymax(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'CFLXYMAX:', cflxymax(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 9 .AND. ifj .EQ. 4 ) THEN

              WRITE ( ndsog ) cflthmax(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'CFLTHMAX:', cflthmax(1:nsea)

#endif

            ELSE IF ( ifi .EQ. 9 .AND. ifj .EQ. 5 ) THEN

              WRITE ( ndsog ) cflkmax(1:nsea)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'CFLMAX:', cflkmax(1:nsea)

#endif

              !

              !     Section 10)

              !

            ELSE IF ( ifi .EQ. 10 ) THEN

              WRITE ( ndsog ) usero(1:nsea,ifj)

#ifdef W3_ASCII

              WRITE ( ndsoa,* ) 'USER0:', usero(1:nsea,ifj)

#endif

              !

            END IF

            !

          ELSE

            !

            !     Start of reading ......

            !

            !     Section 1)

            !

            IF ( ifi .EQ. 1 .AND. ifj .EQ. 1 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) DW(1:NSEA)

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 2 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) CX(1:NSEA)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) CY(1:NSEA)

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 3 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) UA(1:NSEA)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) UD(1:NSEA)

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 4 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) AS(1:NSEA)

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 5 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) WLV(1:NSEA)

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 6 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) ICE(1:NSEA)

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 7 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) BERG(1:NSEA)

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 8 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) TAUA(1:NSEA)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) TAUADIR(1:NSEA)

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 9 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) RHOAIR(1:NSEA)

#ifdef W3_BT4

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 10 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) SED_D50(1:NSEA)

#endif

#ifdef W3_IS2

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 11 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) ICEH(1:NSEA)

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 12 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) ICEF(1:NSEA)

#endif

#ifdef W3_SETUP

            ELSE IF ( ifi .EQ. 1 .AND. ifj .EQ. 13 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) ZETA_SETUP(1:NSEA)

#endif

              !

              !     Section 2)

              !

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 1 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) HS(1:NSEA)

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 2 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) WLM(1:NSEA)

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 3 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) T02(1:NSEA)

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 4 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) T0M1(1:NSEA)

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 5 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) T01(1:NSEA)

            ELSE IF ( (ifi .EQ. 2 .AND. ifj .EQ. 6) .OR.       &

                 (ifi .EQ. 2 .AND. ifj .EQ. 18) ) THEN

              ! Note: TP output is derived from FP field.

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) FP0(1:NSEA)

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 7 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) THM(1:NSEA)

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 8 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) THS(1:NSEA)

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 9 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   thp0(1:nsea)

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 10 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   hsig(1:nsea)

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 11 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   stmaxe(1:nsea)

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 12 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   stmaxd(1:nsea)

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 13 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   hmaxe(1:nsea)

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 14 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   hcmaxe(1:nsea)

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 15 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   hmaxd(1:nsea)

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 16 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   hcmaxd(1:nsea)

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 17 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) WBT(1:NSEA)

            ELSE IF ( ifi .EQ. 2 .AND. ifj .EQ. 19 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   wnmean(1:nsea)

              !

              !     Section 3)

              !

            ELSE IF ( ifi .EQ. 3 .AND. ifj .EQ. 1 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   ef(1:nsea,e3df(2,1):e3df(3,1))

            ELSE IF ( ifi .EQ. 3 .AND. ifj .EQ. 2 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   th1m(1:nsea,e3df(2,2):e3df(3,2))

            ELSE IF ( ifi .EQ. 3 .AND. ifj .EQ. 3 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   sth1m(1:nsea,e3df(2,3):e3df(3,3))

            ELSE IF ( ifi .EQ. 3 .AND. ifj .EQ. 4 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   th2m(1:nsea,e3df(2,4):e3df(3,4))

            ELSE IF ( ifi .EQ. 3 .AND. ifj .EQ. 5 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   sth2m(1:nsea,e3df(2,5):e3df(3,5))

            ELSE IF ( ifi .EQ. 3 .AND. ifj .EQ. 6) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)  &

                   wn(1:nk,1:nsea)

              !

              !     Section 4)

              !

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 1 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   phs(1:nsea,0:noswll)

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 2 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   ptp(1:nsea,0:noswll)

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 3 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   plp(1:nsea,0:noswll)

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 4 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   pdir(1:nsea,0:noswll)

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 5 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   psi(1:nsea,0:noswll)

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 6 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   pws(1:nsea,0:noswll)

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 7 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   pthp0(1:nsea,0:noswll)

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 8  ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   pqp(1:nsea,0:noswll)

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 9  ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   ppe(1:nsea,0:noswll)

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 10 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   pgw(1:nsea,0:noswll)

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 11 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   psw(1:nsea,0:noswll)

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 12 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   ptm1(1:nsea,0:noswll)

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 13 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   pt1(1:nsea,0:noswll)

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 14 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   pt2(1:nsea,0:noswll)

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 15 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   pep(1:nsea,0:noswll)

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 16) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   pwst(1:nsea)

            ELSE IF ( ifi .EQ. 4 .AND. ifj .EQ. 17) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) PNR(1:NSEA)

              !

              !     Section 5)

              !

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 1 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)          &

                   ust(1:nsea)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)          &

                   ustdir(1:nsea)

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 2 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   charn(1:nsea)

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 3 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) CGE(1:NSEA)

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 4 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   phiaw(1:nsea)

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 5 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   tauwix(1:nsea)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   tauwiy(1:nsea)

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 6 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   tauwnx(1:nsea)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   tauwny(1:nsea)

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 7 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   whitecap(1:nsea,1)

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 8 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   whitecap(1:nsea,2)

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 9 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   whitecap(1:nsea,3)

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 10 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   whitecap(1:nsea,4)

            ELSE IF ( ifi .EQ. 5 .AND. ifj .EQ. 11 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   tws(1:nsea)

              !

              !     Section 6)

              !

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 1 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) SXX(1:NSEA)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) SYY(1:NSEA)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) SXY(1:NSEA)

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 2 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   tauox(1:nsea)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   tauoy(1:nsea)

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 3 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   bhd(1:nsea)

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 4 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   phioc(1:nsea)

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 5 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   tusx(1:nsea)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   tusy(1:nsea)

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 6 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   ussx(1:nsea)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   ussy(1:nsea)

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 7 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   prms(1:nsea)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   tpms(1:nsea)

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 8 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)  &

                   us3d(1:nsea,us3df(2):us3df(3))

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)  &

                   us3d(1:nsea,nk+us3df(2):nk+us3df(3))

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ.  9 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)     &

                   p2sms(1:nsea,p2msf(2):p2msf(3))

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 10 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   tauice(1:nsea,1)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   tauice(1:nsea,2)

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 11 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   phice(1:nsea)

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 12 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)  &

                   ussp(1:nsea,1:usspf(2))

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)  &

                   ussp(1:nsea,nk+1:nk+usspf(2))

            ELSE IF ( ifi .EQ. 6 .AND. ifj .EQ. 13 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   tauocx(1:nsea)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   tauocy(1:nsea)


              !

              !     Section 7)

              !

            ELSE IF ( ifi .EQ. 7 .AND. ifj .EQ. 1 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) ABA(1:NSEA)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) ABD(1:NSEA)

            ELSE IF ( ifi .EQ. 7 .AND. ifj .EQ. 2 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) UBA(1:NSEA)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) UBD(1:NSEA)

            ELSE IF ( ifi .EQ. 7 .AND. ifj .EQ. 3 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   bedforms(1:nsea,1)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   bedforms(1:nsea,2)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   bedforms(1:nsea,3)

            ELSE IF ( ifi .EQ. 7 .AND. ifj .EQ. 4 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   phibbl(1:nsea)

            ELSE IF ( ifi .EQ. 7 .AND. ifj .EQ. 5 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   taubbl(1:nsea,1)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   taubbl(1:nsea,2)

              !

              !     Section 8)

              !

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 1 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   mssx(1:nsea)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   mssy(1:nsea)

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 2 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   mscx(1:nsea)

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   mscy(1:nsea)

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 3 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   mssd(1:nsea)

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 4 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   mscd(1:nsea)

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 5 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) QP(1:NSEA)

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 6 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) QKK(1:NSEA)

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 7 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) SKEW(1:NSEA)

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 8 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) EMBIA1(1:NSEA)

            ELSE IF ( ifi .EQ. 8 .AND. ifj .EQ. 9 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR) EMBIA2(1:NSEA)

              !

              !     Section 9)

              !

            ELSE IF ( ifi .EQ. 9 .AND. ifj .EQ. 1 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   dtdyn(1:nsea)

            ELSE IF ( ifi .EQ. 9 .AND. ifj .EQ. 2 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   fcut(1:nsea)

            ELSE IF ( ifi .EQ. 9 .AND. ifj .EQ. 3 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   cflxymax(1:nsea)

            ELSE IF ( ifi .EQ. 9 .AND. ifj .EQ. 4 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   cflthmax(1:nsea)

            ELSE IF ( ifi .EQ. 9 .AND. ifj .EQ. 5 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   cflkmax(1:nsea)

              !

              !     Section 10)

              !

            ELSE IF ( ifi .EQ. 10 ) THEN

              READ (ndsog,END=801,ERR=802,IOSTAT=IERR)         &

                   usero(1:nsea,ifj)

            END IF

            !

            ! End of test on WRITE/READ:

            !

          END IF

          !

          ! End of test on  FLOGRD(IFI,IFJ):

          !

        END IF

        !

        ! End of IFI and IFJ loops

        !

      END DO

    END DO

    !

    ! Flush the buffers for write

    !

    IF ( WRITE ) CALL flush ( ndsog )

    !

    IF(ofiles(1) .EQ. 1) CLOSE(ndsog)

    !

#ifdef W3_MPI

    CALL w3seta ( igrd, ndse, ndst )

#endif

    !

    RETURN

    !

    ! Escape locations read errors

    !

800 CONTINUE

    WRITE (ndse,1000) ierr

    CALL extcde ( 41 )

    !

801 CONTINUE

    WRITE (ndse,1001)

    CALL extcde ( 42 )

    !

802 CONTINUE

    WRITE (ndse,1002) ierr

    CALL extcde ( 43 )

    !

803 CONTINUE

    iotst  = -1

#ifdef W3_T

    WRITE (ndst,9020)

#endif

    RETURN

    !

    ! Formats

    !

900 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOGO :'/                &

         '     ILEGAL INXOUT VALUE: ',a/)

901 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOGO :'/                &

         '     MIXED READ/WRITE, LAST REQUEST: ',a/)

902 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOGO :'/                &

         '     ILEGAL IDSTR, READ : ',a/                        &

         '                  CHECK : ',a/)

903 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOGO :'/                &

         '     ILEGAL VEROGR, READ : ',a/                       &

         '                   CHECK : ',a/)

904 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOGO :'/                &

         '     DIFFERENT NUMBER OF FIELDS, FILE :',i8,i8/       &

         '                              PROGRAM :',i8,i8/)

905 FORMAT (/' *** WAVEWATCH III WARNING IN W3IOGO :'/              &

         '     ILEGAL GNAME, READ : ',a/                        &

         '                  CHECK : ',a/)

906 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOGO :'/                &

         '     ILEGAL NOSWLL, READ : ',i4/                      &

         '                   CHECK : ',i4/)

    !

    !  999 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOGO :'/                &

    !               '     PLEASE UPDATE FIELDS !!! '/)

    !

1000 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOGO : '/               &

         '     ERROR IN OPENING FILE'/                          &

         '     IOSTAT =',i5/)

1001 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOGO : '/               &

         '     PREMATURE END OF FILE'/)

1002 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOGO : '/               &

         '     ERROR IN READING FROM FILE'/                     &

         '     IOSTAT =',i5/)

    !

#ifdef W3_T

9000 FORMAT (' TEST W3IOGO : IPASS =',i4,' INXOUT = ',a,          &

         ' WRITE = ',l1,' UNIT =',i3/                         &

         '               IGRD =',i3,' FEXT = ',a)

9001 FORMAT (' TEST W3IOGO : OPENING NEW FILE [',a,']')

9002 FORMAT (' TEST W3IOGO : TEST PARAMETERS:'/                   &

         '       IDSTR :  ',a/                                &

         '      VEROGR :  ',a/                                &

         '       GNAME :  ',a/                                &

         '        NSEA :',i6/                                 &

         '       NX,NY : ',i9,i12/                            &

         '       UNDEF : ',f8.2)

9003 FORMAT (' TEST W3IOGO : TIME  :',i9.8,i7.6/                  &

         '               FLAGS :',20l2,1x,20l2/               &

         '                      ',20l2,2x,20l2/               &

         '                      ',20l2,2x,20l2/               &

         '                      ',20l2,2x,20l2/               &

         '                      ',20l2,2x,20l2)

9010 FORMAT (' TEST W3IOGO : PROC = ',l1,' FOR ',a)

9020 FORMAT (' TEST W3IOGO : END OF FILE REACHED')

#endif

    !/

    !/ End of W3IOGO ----------------------------------------------------- /

    !/

  END SUBROUTINE w3iogo

  !/

  !/ ------------------------------------------------------------------- /

  SUBROUTINE calc_u3stokes ( A , USS_SWITCH )

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           H. L. Tolman            |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         10-Jan-2017 |

    !/                  +-----------------------------------+

    !/

    !/    10-Jan-2017 : Separate Stokes drift calculation  ( version 6.01 )

    !/

    !  1. Purpose :

    !

    !     This code is built for the purpose of outputting Stokes drift

    !       related parameters that can be utilized to obtain full

    !       Stokes drift profiles external to the wave model.

    !

    !     Option 1: USS_SWITCH == 1

    !               This method is for outputing the Stokes drift frequency

    !               spectrum for spectral frequency bands as defined by the

    !               WW3 computation spectral frequency grid.

    !               Output Quantity: Stokes drift frequency spectrum [m/s/Hz]

    !                                X and Y componenets.

    !

    !     Option 2: USS_SWITCH == 2

    !               This method is for outputing the surface Stokes drift

    !               for a specified frequency partition/band of the

    !               wave spectrum.  These partitions do not need to be

    !               matched to WW3's computation spectral frequency grid,

    !               and will rather sum the contributions of the WW3 bands

    !               into the output partition.  The partitions are defined

    !               in the ww3_grid.inp namelist section.

    !               Output Quantity: Stokes drift surface velocity [m/s]

    !                                X and Y components

    !                                For each partition (up to 25).

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !       A       R.A.   I   Input spectra. Left in par list to change

    !                          shape.

    !      USS_SWITCH  I   I   Switch if computing US3D (spectral) or USSP

    !                           (partitions)

    !     ----------------------------------------------------------------

    !

    !

    !  4. Subroutines used :

    !

    !     See module documentation.

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !     None.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/SHRD  Switch for shared / distributed memory architecture.

    !     !/DIST  Id.

    !

    !     !/OMPG  OpenMP compiler directive for loop splitting.

    !

    !     !/S     Enable subroutine tracing.

    !     !/T     Test output.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

    USE constants, ONLY: tpiinv, grav, tpi

    USE w3gdatmd,  ONLY: dden, dsii, xfr, sig, nk, nth, nseal,    &

         ecos, esin, us3df, usspf, ussp_wn

    USE w3adatmd,  ONLY: cg, wn, dw

    USE w3adatmd,  ONLY: ussx, ussy,  us3d, ussp

    USE w3odatmd, ONLY: iaproc, naproc

    USE w3parall, ONLY: init_get_isea

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    REAL, INTENT(IN)        :: A(NTH,NK,0:NSEAL)

    INTEGER, INTENT(IN)     :: USS_SWITCH

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

    INTEGER                 :: IK, ITH, ISEA, JSEA

    INTEGER                 :: IKST, IKFI, IB

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    REAL                    :: FACTOR, FKD,KD

    REAL                    :: ABX(NSEAL), ABY(NSEAL), USSCO

    REAL                    :: MINDIFF

    INTEGER                 :: Spc2Bnd(NK)

    !/

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'CALC_U3STOKES')

#endif

    !

    ! 1.  Initialize storage arrays -------------------------------------- *

    !

    ! 2.  Integral over discrete part of spectrum ------------------------ *

    !

    !Two options ----------------------------------------------------|

    ! USS_SWITCH == 1 -> Old option, Stokes drift integrated in same |

    !                    wavenumber bands as model integrates.       |

    ! USS_SWITCH == 2 -> New option, Stokes drift integrated in a    |

    !                    defined number (NP) of user specified       |

    !                    partitions, where NP and the frequency      |

    !                    ranges for each partition can be user       |

    !                    defined at run-time.                        |

    !----------------------------------------------------------------|


    if (uss_switch==1) then

      ikst=us3df(2)!Start at US3DF(2)

      ikfi=us3df(3)!End at US3DF(3)

    ELSEif (uss_switch==2) then

      ikst=1  ! Start at 1

      ikfi=nk ! End at NK

    ENDIF


    ! Initialize US3D/USSP

    IF (uss_switch.eq.1) then

      us3d(:,:)=0.0

    ELSEIF (uss_switch.eq.2) then

      ussp(:,:)=0.0

    ENDIF

    DO ik=ikst,ikfi   !1, NK

      !

      ! 2.a Initialize energy in band

      !

      abx    = 0.

      aby    = 0.

      !

      ! 2.b Integrate energy in band

      !

      DO ith=1, nth

        !

#ifdef W3_OMPG

        !$OMP PARALLEL DO PRIVATE(JSEA)

#endif

        !

        DO jsea=1, nseal

          abx(jsea)  = abx(jsea) + a(ith,ik,jsea)*ecos(ith)

          aby(jsea)  = aby(jsea) + a(ith,ik,jsea)*esin(ith)

        END DO

        !

#ifdef W3_OMPG

        !$OMP END PARALLEL DO

#endif

        !

      END DO

      !

      ! 2.c Finalize integration over band and update mean arrays

      !

      !

#ifdef W3_OMPG

      !$OMP PARALLEL DO PRIVATE(JSEA,ISEA,FACTOR,KD,FKD,USSCO,MINDIFF,IB)

#endif

      !

      DO jsea=1, nseal

        CALL init_get_isea(isea, jsea)

        factor       = dden(ik) / cg(ik,isea)

        !

        ! Deep water limits

        !

        kd    = max( 0.001 , wn(ik,isea) * dw(isea) )

        IF ( kd .LT. 6. ) THEN

          fkd       = factor / sinh(kd)**2

          ussco=fkd*sig(ik)*wn(ik,isea)*cosh(2.*kd)

        ELSE

          ussco=factor*sig(ik)*2.*wn(ik,isea)

        END IF

        !

        !

        !USSX(JSEA)  = USSX(JSEA) + ABX(JSEA)*USSCO

        !USSY(JSEA)  = USSY(JSEA) + ABY(JSEA)*USSCO

        !

        ! Fills the 3D Stokes drift spectrum array or surface Stokes partitions

        !

        IF (uss_switch==1) THEN

          !Old method fills into WW3 bands

          IF (ik.GE.us3df(2).and.ik.LE.us3df(3)) then

            us3d(jsea,ik)    =  abx(jsea)*ussco/(dsii(ik)*tpiinv)

            us3d(jsea,nk+ik) =  aby(jsea)*ussco/(dsii(ik)*tpiinv)

          ENDIF

        ELSEIF (uss_switch==2) THEN

          ! Match each spectral component to the nearest partition

          mindiff=1.e8

          spc2bnd(ik) = 1

          mindiff=abs(ussp_wn(1)-wn(ik,isea))

          DO ib=2,usspf(2)

            IF (mindiff .gt. abs(ussp_wn(ib)-wn(ik,isea))) then

              spc2bnd(ik) = ib

              mindiff = abs(ussp_wn(ib)-wn(ik,isea))

            ENDIF

          ENDDO

          !Put spectral energey into whichever band central wavenumber fits in

          ussp(jsea,spc2bnd(ik))    =  ussp(jsea,spc2bnd(ik)) + abx(jsea)*ussco

          ussp(jsea,nk+spc2bnd(ik)) =  ussp(jsea,nk+spc2bnd(ik)) + aby(jsea)*ussco

        ENDIF

      END DO

#ifdef W3_OMPG

      !$OMP END PARALLEL DO

#endif

    END DO

    !

    RETURN

    !

    !/ End of CALC_U3STOKES

    !----------------------------------------------------- /

    !/

  END SUBROUTINE calc_u3stokes

  !/

  !/ ------------------------------------------------------------------- /

  SUBROUTINE calc_wbt (A)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           Q. Liu                  |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         24-Aug-2018 |

    !/                  +-----------------------------------+

    !/

    !/    24-Aug-2018 : Origination.                       ( version 6.06 )

    !/

    !  1. Purpose :

    !

    !     Estimate the dominant wave breaking probability b_T based on

    !     the empirical parameterization proposed by Babanin et al. (2001).

    !     From their Fig. 12, we have

    !

    !         b_T = 85.1 * [(εp - 0.055) * (1 + H_s/d)]^2.33,

    !

    !     where ε is the significant steepness of the spectral peak, H_s is

    !     the significant wave height, d is the water depth.

    !

    !     For more details, please see

    !         Banner et al.  2000: JPO,      30,  3145 -  3160.

    !         Babanin et al. 2001: JGR, 106(C6), 11569 - 11676.

    !

    !  2. Method :

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    !       A       R.A.   I   Input wave action spectra N(j, θ, k)

    !     ----------------------------------------------------------------

    !

    !  4. Subroutines used :

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

    !     ----------------------------------------------------------------

    !      W3OUTG    Subr. Public   Calculate mean parameters.

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !     None.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S     Enable subroutine tracing.

    !     !/T     Test output.

    !

    ! 10. Source code :

    !

    !/ ------------------------------------------------------------------- /

    USE w3dispmd, ONLY: wavnu1

    USE w3adatmd, ONLY: u10, u10d, wbt

    USE w3adatmd, ONLY: cg, wn, dw

    USE w3gdatmd, ONLY: nk, nth, nseal, sig, esin, ecos, dth, dsii, &

         fte, xfr, mapsf, mapsta, dmin

    USE w3gdatmd, ONLY: btbeta

    USE w3parall, ONLY: init_get_isea

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    IMPLICIT NONE

    !

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    REAL, INTENT(IN)     :: A  (NTH, NK, 0:NSEAL)

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    !

    INTEGER              :: FPOPT = 0

    !

    INTEGER              :: IK, ITH, ISEA, JSEA, IKM, IKL, IKH, IX, IY

    REAL                 :: TDPT, TU10, TUDIR, SINU, COSU, TC, TFORCE

    REAL                 :: ESIG(NK) ! E(σ)

    REAL                 :: FACTOR, ET, HS, ETP, HSP, SIGP, KP, &

         CGP, WSTP

    REAL                 :: XL, XH, XL2, XH2, EL, EH, DENOM

    REAL                 :: TWBT

    !/

    !/ ------------------------------------------------------------------- /

    !/

#ifdef W3_S

    CALL strace (ient, 'CALC_WBT')

#endif

    !

    DO jsea = 1, nseal

      ! JSEA 2 ISEA

      CALL init_get_isea(isea, jsea)

      !

      ! check the status of this grid point [escape if this point is excluded]

      !

      ix = mapsf(isea,1)

      iy = mapsf(isea,2)

      IF ( mapsta(iy,ix) .LE. 0 ) cycle

      !

      ! Wind info. is required to select wind sea partition from the wave

      ! spectrum. Two wind velocities are availabe:

      ! - U10 & U10D   (w3adatmd)

      ! - UST & USTDIR (w3wdatmd)

      !     * U10D & USTDIR are not really the same when swell are present.

      !

      ! Following Janssen et al. (1989) and Bidlot (2001), spectral components

      ! are considered to be subject to local wind forcing when

      !

      !          c / [U cos(θ - φ)] < β,

      !

      ! where c is the phase velocity c = σ/k, φ is the wind direction, U is

      ! the wind speed U10, (sometimes approximated by U10≅ 28 * ust), β is

      ! the constant forcing parameter with β∈ [1.0, 2.0]. By default, we use

      ! β = 1.2(Bidlot 2001).

      !

      tdpt  = max(dw(isea), dmin)          ! water depth d

      tu10  = u10(isea)                    ! wind velocity U10

      tudir = u10d(isea)                   ! wind direction φ (rad)

      sinu  = sin(tudir)                   ! sinφ

      cosu  = cos(tudir)                   ! cosφ

      !

      esig  = 0.                           ! E(σ)

      et    = 0.                           ! ΣE(σ)δσ

      etp   = 0.                           ! ΣE(σ)δσ at peak only

      !

      DO ik = 1, nk

        tc     = sig(ik) / wn(ik, isea)  ! phase velocity c=σ/k

        factor = sig(ik) / cg(ik, isea)  ! σ / cg

        factor = factor * dth            ! σ / cg * δθ

        !

        DO ith = 1, nth

          tforce = tc - tu10 * (cosu*ecos(ith)+sinu*esin(ith)) &

               * btbeta


          IF (tforce .LT. 0.) THEN ! wind sea component

            esig(ik) = esig(ik) + a(ith, ik, jsea) * factor

          ENDIF

        ENDDO ! ITH

        !

      ENDDO ! IK

      !

      ! ESIG is E(σ) of the wind sea after filtration of any background swell.

      ! Now we need to get Hs & σp for the wind sea spectrum.

      ! FTE    = 0.25 * SIG(NK) * DTH * SIG(NK) [ww3_grid.ftn]

      !

      et = sum(esig * dsii)

      et = et + esig(nk) * fte / (dth * sig(nk))  ! FTE: add tail

      hs = 4. * sqrt(max(0., et))                 ! Hs

      !

      ! Get σp from E(σ)

      !

      ! Here we have tried three different ways to calculate FP:

      !

      ! FPOPT = 0: fp defined by Young (1999, p. 239)

      ! FPOPT = 1: parabolic fit around the discrete peak frequency, as used

      !            by ww3_outp

      ! FPOPT = 2: discrete peak frequency

      !

      ! When the discrete peak frequency is used:

      ! * For XFR = 1.1, the **discrete** peak region [0.7σp, 1.3σp] will be

      !     {0.75, 0.83, 0.91, 1., 1.1, 1.21, 1.33}σp,

      ! * and for XFR = 1.07, the **discrete** peak region becomes

      !     {0.71, 0.76, 0.82, 0.87, 0.93, 1., 1.07, 1.14, 1.23, 1.31}σp.

      !

      ! Thus, a good approximation to the range [0.7σp, 1.3σp] is guranteed

      ! by each XFR. I however found using the discrete peak frequency yielded

      ! step-wise results. According to my test, the smoothest results were

      ! obtained with FPOPT = 0. For simplicity, the δσ values (DSII) are

      ! not modified.

      !

      ikm   = maxloc(esig, 1)                     ! index for σp

      !

      IF (fpopt .EQ. 0) THEN

        !             FP defined in Ian's book

        sigp = sum(esig**4. * sig(1:nk) * dsii) /  &

             max(1e-10, sum(esig**4. * dsii))

        !

      ELSE IF (fpopt .EQ. 1) THEN

        !             Parabolic fit around the discrete peak (ww3_outp.ftn)

        xl    = 1./xfr - 1.

        xh    = xfr - 1.

        xl2   = xl**2.

        xh2   = xh**2.

        ikl   = max(  1 , ikm-1 )

        ikh   = min( nk , ikm+1 )

        el    = esig(ikl) - esig(ikm)

        eh    = esig(ikh) - esig(ikm)

        denom = xl*eh - xh*el

        sigp  = sig(ikm) * (1. + 0.5 * ( xl2*eh - xh2*el) &

             / sign(max(abs(denom), 1.e-15), denom)) ! σp

        !

      ELSE IF (fpopt .EQ. 2) THEN

        !             Discrete peak (Give stepwise results, not used by default)

        sigp  = sig(ikm)

      ENDIF

      !

      ! kp from σp (linear dispersion)

      !

      ! N(k, θ) at first step is zero →  σp=0 → floating divided by zero error

      IF (sigp < 1e-6) sigp = sig(nk)     ! Hsp & b_T will be still 0.

      !

      CALL wavnu1 (sigp, tdpt, kp, cgp)

      !

      !                         { /1.3σp         }1/2

      ! peak wave height Hp = 4 { |      E(σ) dσ }

      !                         { /0.7σp         }

      !

      DO ik = 1, nk

        IF ( (sig(ik) >= 0.7 * sigp) .AND. &

             (sig(ik) <= 1.3 * sigp) ) THEN

          etp = etp + esig(ik) * dsii(ik)

        ENDIF

      ENDDO ! IK

      hsp = 4. * sqrt(max(0., etp))

      !

      ! significant steepness of the peak region εp

      !

      wstp = 0.5 * kp * hsp

      !

      ! Dominant wave breaking b_T

      !

      twbt = 85.1 * (max(0.0, wstp - 0.055) * (1 + hs/tdpt))**2.33

      wbt(jsea) = min(1.0, twbt)

      !

    ENDDO ! JSEA

    !/

    !/ End of  CALC_WBT -------------------------------------------------- /

    !/

  END SUBROUTINE calc_wbt

  !/ ------------------------------------------------------------------- /

  !/

      SUBROUTINE secondhh(NKHF,FAC0,FAC1,FAC2,FAC3)

!----------------------------------------------------------------


!**** *SECONDHH* - COMPUTATION OF SECOND ORDER HARMONICS AND

!                  RELEVANT TABLES FOR THE ALTIMETER CORRECTIONS.


!     P.A.E.M. JANSSEN


!     PURPOSE.

!     ---------


!          COMPUTE SECOND HARMONICS


!**   INTERFACE.

!     ----------


!          *CALL* *SECONDHH*


!     METHOD.

!     -------


!          SEE REFERENCE.


!     EXTERNALS.

!     ----------


!         VMIN_D

!         VPLUS_D


!     REFERENCES.

!     -----------


!          V E ZAKHAROV(1967)


!-------------------------------------------------------------------


!-------------------------------------------------------------------

USE constants, ONLY: grav, tpi

USE w3gdatmd,  ONLY: nk, nth, xfr, sig, th, dth, ecos, esin

      IMPLICIT NONE

 !     REAL(KIND=4) :: VMIN_D,VPLUS_D


      INTEGER, INTENT(IN) :: NKHF

      REAL(KIND=4), dimension(nth,nth,nkhf,nkhf), INTENT(OUT)  :: fac0, fac1, fac2, fac3

      REAL(KIND=4), parameter   :: fratio = 1.1


      INTEGER :: M, K1, M1, K2, M2


      REAL(KIND=4), parameter :: del1=1.0e-8

      REAL(KIND=4), parameter :: zconst = 0.0281349


      !REAL(KIND=4) :: VMIN_D, VPLUS_D

      REAL(KIND=4) :: co1

      REAL(KIND=4) :: xk1, xk1sq, xk2, xk2sq, xk3

      REAL(KIND=4) :: cosdiff

      REAL(KIND=4) :: x12, x13, x32, om1, om2, om3, f1, f2, f3

      REAL(KIND=4) :: vm, vp

      REAL(KIND=4) :: delom1, delom2

      REAL(KIND=4) :: delom321, delom312

      REAL(KIND=4) :: c22, s22


      REAL(KIND=4), dimension(nth,nth,nkhf,nkhf) :: b

      REAL(KIND=4), dimension(:), ALLOCATABLE:: fak, sighf, dfimhf


!-----------------------------------------------------------------------


!*    1. INITIALISE RELEVANT QUANTITIES.


      ALLOCATE(fak(nkhf))

      ALLOCATE(sighf(nkhf))

      ALLOCATE(dfimhf(nkhf))


      sighf(1)  = sig(1)

      DO m=2,nkhf

        sighf(m) = xfr*sighf(m-1)

      ENDDO


      DO m=1,nkhf

         fak(m) = (sighf(m))**2/grav

      ENDDO


      co1 = 0.5*(xfr-1.)*dth

      dfimhf(1) = co1*sighf(1)

      DO m=2,nkhf-1

         dfimhf(m)=co1*(sighf(m)+sighf(m-1))

      ENDDO

      dfimhf(nkhf)=co1*sighf(nkhf-1)


      DO m2=1,nkhf

        xk2 = fak(m2)

        xk2sq = fak(m2)**2

        DO  m1=1,nkhf

          xk1 = fak(m1)

          xk1sq = fak(m1)**2

          DO k1=1,nth

            DO k2=1,nth

              cosdiff = cos(th(k1)-th(k2))

              x12 = xk1*xk2*cosdiff

              xk3 = xk1sq + xk2sq +2.0*x12 +del1

              xk3 = sqrt(xk3)

              x13 = xk1sq+x12

              x32 = x12+xk2sq

              om1 = sqrt(grav*xk1)

              om2 = sqrt(grav*xk2)

              om3 = sqrt(grav*xk3)

              f1 = sqrt(xk1/(2.0*om1))

              f2 = sqrt(xk2/(2.0*om2))

              f3 = sqrt(xk3/(2.0*om3))

              vm = tpi*vmin_d(xk3,xk1,xk2,x13,x32,x12,om3,om1,om2)

              vp = tpi*vplus_d(-xk3,xk1,xk2,-x13,-x32,x12,om3,om1,om2)

              delom1 = om3-om1-om2+del1

              delom2 = om3+om1+om2+del1

              fac0(k1,k2,m1,m2) = -f3/(f1*f2)*(vm/(delom1)+             &

     &                            vp/(delom2))

            ENDDO

          ENDDO

        ENDDO

      ENDDO


      DO m2=1,nkhf

        xk2 = fak(m2)

        xk2sq = fak(m2)**2

        DO  m1=1,nkhf

          xk1 = fak(m1)

          xk1sq = fak(m1)**2

          DO k1=1,nth

            DO k2=1,nth

              cosdiff = cos(th(k1)-th(k2))

              x12 = xk1*xk2*cosdiff

              xk3 = xk1sq + xk2sq - 2.*x12 + del1

              xk3 = sqrt(xk3)

              x13 = xk1sq-x12

              x32 = x12-xk2sq

              om1 = sqrt(grav*xk1)

              om2 = sqrt(grav*xk2)

              om3 = sqrt(grav*xk3)+del1

              f1 = sqrt(xk1/(2.0*om1))

              f2 = sqrt(xk2/(2.0*om2))

              f3 = sqrt(abs(xk3)/(2.0*om3))

              vm = tpi*vmin_d(xk1,xk3,xk2,x13,x12,x32,om1,om3,om2)

              vp = tpi*vmin_d(xk2,-xk3,xk1,-x32,x12,-x13,om2,om3,om1)

              delom321 = om3+om2-om1+del1

              delom312 = om3+om1-om2+del1

              b(k1,k2,m1,m2) = -f3/(f1*f2)*(vm/(delom321)+              &

     &                         vp/(delom312))

            ENDDO

          ENDDO

        ENDDO

      ENDDO


      DO m2=1,nkhf

        xk2sq = fak(m2)**2

        DO m1=1,nkhf

          xk1sq = fak(m1)**2

          DO k2=1,nth

            DO k1=1,nth

              c22 = fac0(k1,k2,m1,m2)+b(k1,k2,m1,m2)

              s22 = b(k1,k2,m1,m2)-fac0(k1,k2,m1,m2)

              fac1(k1,k2,m1,m2) =                                       &

     &             (xk1sq*ecos(k1)**2 + xk2sq*ecos(k2)**2)*c22        &

     &             -fak(m1)*fak(m2)*ecos(k1)*ecos(k2)*s22

              fac2(k1,k2,m1,m2) =                                       &

     &             (xk1sq*esin(k1)**2 + xk2sq*esin(k2)**2)*c22        &

     &             -fak(m1)*fak(m2)*esin(k1)*esin(k2)*s22

              fac3(k1,k2,m1,m2) =                                       &

     &             (xk1sq*esin(k1)*ecos(k1) +                         &

     &              xk2sq*esin(k2)*ecos(k2))*c22                      &

     &             -fak(m1)*fak(m2)*ecos(k1)*esin(k2)*s22

              fac0(k1,k2,m1,m2) = c22

            ENDDO

          ENDDO

        ENDDO

      ENDDO


     CONTAINS


!-----------------------------------------------------------------------


     REAL(KIND=4) function vmin_d(xi,xj,xk,xij,xik,xjk,xoi,xoj,xok)


!     PETER JANSSEN


!     PURPOSE.

!     --------


!              GIVES NONLINEAR TRANSFER COEFFICIENT FOR THREE

!              WAVE INTERACTIONS OF DEEP-WATER WAVES IN THE

!              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV)


!     INTERFACE.

!     ----------

!              *VMIN_D(XI,XJ,XK)*

!                      *XI*  - WAVE NUMBER

!                      *XJ*  - WAVE NUMBER

!                      *XK*  - WAVE NUMBER

!     METHOD.

!     -------

!              NONE


!     EXTERNALS.

!     ----------

!              NONE.


!***  1. DETERMINE NONLINEAR TRANSFER.

!     --------------------------------

      IMPLICIT NONE

      REAL, INTENT(IN) :: xi, xj, xk, xij, xik, xjk, xoi, xoj, xok

      REAL :: ri, rj, rk, oi, oj, ok, sqijk, sqikj, sqjki


      ri=abs(xi)+del1

      rj=abs(xj)+del1

      rk=abs(xk)+del1

      oi=xoi+del1

      oj=xoj+del1

      ok=xok+del1

      sqijk=sqrt(oi*oj*rk/(ok*ri*rj))

      sqikj=sqrt(oi*ok*rj/(oj*ri*rk))

      sqjki=sqrt(oj*ok*ri/(oi*rj*rk))

      vmin_d=zconst*( (xij-ri*rj)*sqijk + (xik-ri*rk)*sqikj             &

     &                + (xjk+rj*rk)*sqjki )


      END FUNCTION vmin_d


!-----------------------------------------------------------------------


      REAL(kind=4) function vplus_d(xi,xj,xk,xij,xik,xjk,xoi,xoj,xok)


!***  *VPLUS_D*  DETERMINES THE NONLINEAR TRANSFER COEFFICIENT FOR THREE

!                WAVE INTERACTIONS OF DEEP-WATER WAVES.


!     PETER JANSSEN


!     PURPOSE.

!     --------


!              GIVES NONLINEAR TRANSFER COEFFICIENT FOR THREE

!              WAVE INTERACTIONS OF GRAVITY-CAPILLARY WAVES IN THE

!              IDEAL CASE OF NO CURRENT. (CF.ZAKHAROV)


!     INTERFACE.

!     ----------

!              *VPLUS_D(XI,XJ,XK)*

!                        *XI*  - WAVE NUMBER

!                        *XJ*  - WAVE NUMBER

!                        *XK*  - WAVE NUMBER

!     METHOD.

!     -------

!              NONE


!     EXTERNALS.

!     ----------

!              NONE.


!***  1. DETERMINE NONLINEAR TRANSFER.

!     --------------------------------


      IMPLICIT NONE

      REAL, INTENT(IN) :: xi, xj, xk, xij, xik, xjk, xoi, xoj, xok

      REAL :: ri, rj, rk, oi, oj, ok, sqijk, sqikj, sqjki


      ri=abs(xi)+del1

      rj=abs(xj)+del1

      rk=abs(xk)+del1

      oi=xoi+del1

      oj=xoj+del1

      ok=xok+del1

      sqijk=sqrt(oi*oj*rk/(ok*ri*rj))

      sqikj=sqrt(oi*ok*rj/(oj*ri*rk))

      sqjki=sqrt(oj*ok*ri/(oi*rj*rk))

      vplus_d=zconst*( (xij+ri*rj)*sqijk + (xik+ri*rk)*sqikj            &

     &               + (xjk+rj*rk)*sqjki )


      END FUNCTION vplus_d

!     -----------------------------------------------------------------


      END SUBROUTINE secondhh

  !/ ------------------------------------------------------------------- /

  !/

      SUBROUTINE skewness(A)


!--------------------------------------------------------------------


!*****SKEWNESS** COMPUTES PARAMETERS OF THE NEARLY-GAUSSIAN

!             DISTRIBUTION OF OCEAN WAVES AT A FIXED GRID POINT.


!     P.JANSSEN JULY 1997


!     PURPOSE

!     -------

!             DETERMINES SKEWNESS PARAMETERS IN ORDER TO OBTAIN

!             CORRECTION ON ALTIMETER WAVE HEIGHT.


!     INTERFACE

!     ---------

!             *CALL* *SKEWNESS(IU06,F1,NCOLL,XKAPPA1,DELH_ALT)*


!     METHOD

!     ------

!             EVALUATE DEVIATIONS FROM GAUSSIANITY FOLLOWING THE WORK

!             OF SROKOSZ AND LONGUET-HIGGINS. FOR SECOND ORDER

!             CORRECTIONS TO SURFACE ELEVATION THE APPROACH OF

!             ZAKHAROV HAS BEEN USED.


!     EXTERNALS

!     ---------

!             NONE


!     REFERENCES

!     ----------

!             M.A. SROKOSZ, J.G.R.,91,995-1006(1986)

!             V.E. ZAKHAROV, HAMILTONIAN APPROACH(1967)

!--------------------------------------------------------------------


!--------------------------------------------------------------------

!      *TH*        REAL      DIRECTIONS IN RADIANS.

USE constants, ONLY: grav, tpi, tpiinv

USE w3gdatmd,  ONLY: nk, nth, xfr, sig, dth, ecos, esin, nseal

USE w3parall,  ONLY: init_get_isea

USE w3adatmd,  ONLY: cg, skew, embia1, embia2


    IMPLICIT NONE


    REAL, INTENT(IN)        :: A(NTH,NK,0:NSEAL)


    INTEGER :: NKHF

    REAL(KIND=4), dimension(:,:,:,:) , ALLOCATABLE:: fac0,fac1,fac2,fac3


    INTEGER :: M, K, M1, K1, M2, K2, I, J

    INTEGER :: MSTART, JSEA


    REAL(KIND=4) :: conx, delta

    REAL(KIND=4) :: fh, delf, xk1

    REAL(KIND=4) :: xpi, xpj, xpk, xn, xfac, co1

    REAL(KIND=4), dimension(:,:), ALLOCATABLE :: f2

    REAL(KIND=4), dimension(0:3,0:2,0:2) :: xmu, xlambda

    REAL(KIND=4), dimension(:) , ALLOCATABLE::  sighf, dfimhf, fak


! ----------------------------------------------------------------------


    nkhf=nk+13 ! same offset as in ECWAM


    ALLOCATE(fac0(nth,nth,nkhf,nkhf))

    ALLOCATE(fac1(nth,nth,nkhf,nkhf))

    ALLOCATE(fac2(nth,nth,nkhf,nkhf))

    ALLOCATE(fac3(nth,nth,nkhf,nkhf))


    CALL secondhh(nkhf,fac0,fac1,fac2,fac3)


    ALLOCATE(f2(nth,nkhf))

    ALLOCATE(sighf(nkhf), dfimhf(nkhf), fak(nkhf))


!     1. COMPUTATION OF FREQUENCY-DIRECTION INCREMENT

!     -----------------------------------------------


    mstart = 1


#ifdef W3_OMPG

        !$OMP PARALLEL DO PRIVATE(JSEA)

#endif

    DO jsea=1, nseal

      xmu(:,:,:) = 0.0

      DO k=1,nth

        DO m=1,nk

          conx = tpiinv / sig(m) * cg(m,jsea)

          f2(k,m)=a(k,m,jsea)/ conx

          END DO

        END DO


      sighf(1)  = sig(1)

      DO m=2,nkhf

        sighf(m) = xfr*sighf(m-1)

      ENDDO


      co1 = 0.5*(xfr-1.)*dth*tpiinv

      dfimhf(1) = co1*sighf(1)        ! this is DF*DTH

      DO m=2,nkhf-1

        dfimhf(m)=co1*(sighf(m)+sighf(m-1))

      ENDDO

      dfimhf(nkhf)=co1*sighf(nkhf-1)


      DO m=1,nkhf

        fak(m) = (sighf(m))**2/grav

      ENDDO


! Deals with the tail ...

      DO m=nk+1,nkhf

        fh=(sighf(nk)/sighf(m))**5

        DO k=1,nth

          f2(k,m)=f2(k,nk)*fh

        ENDDO

      ENDDO


!     2. COMPUTATION OF THE SKEWNESS COEFFICIENTS

!     --------------------------------------------


      DO m1=mstart,nkhf

        DO m2=mstart,nkhf

          DO k1=1,nth

            DO k2=1,nth

              delf = dfimhf(m1)*dfimhf(m2)*f2( k1,m1)*f2(k2,m2)

              xmu(3,0,0) = xmu(3,0,0)+3.0*fac0(k1,k2,m1,m2)*delf

              xmu(1,2,0) = xmu(1,2,0)+fac1(k1,k2,m1,m2)*delf

              xmu(1,0,2) = xmu(1,0,2)+fac2(k1,k2,m1,m2)*delf

              xmu(1,1,1) = xmu(1,1,1)+fac3(k1,k2,m1,m2)*delf

            ENDDO

          ENDDO

        ENDDO

      ENDDO


      DO k1=1,nth

        DO m1=mstart,nkhf

          xk1 = fak(m1)**2

          delf = dfimhf(m1)*f2(k1,m1)

          xmu(2,0,0) = xmu(2,0,0) + delf

          xmu(0,2,0) = xmu(0,2,0) + xk1*ecos(k1)**2*delf

          xmu(0,0,2) = xmu(0,0,2) + xk1*esin(k1)**2*delf

          xmu(0,1,1) = xmu(0,1,1) + xk1*ecos(k1)*esin(k1)*delf

        ENDDO

      ENDDO


!     3. COMPUTATION OF THE NORMALISED SKEWNESS COEFFICIENTS

!     ------------------------------------------------------


      DO i=0,3

        xpi = 0.5*float(i)

        DO j=0,2

          xpj = 0.5*float(j)

          DO k=0,2

            xpk = 0.5*float(k)

            xn = xmu(2,0,0)**xpi*xmu(0,2,0)**xpj*xmu(0,0,2)**xpk  ! denom in Srokosz eq. 11

            IF (xn .NE. 0) THEN

              xlambda(i,j,k) = xmu(i,j,k)/xn

            ELSE

              xlambda(i,j,k) = 0

            END IF

          END DO

        END DO

      END DO

      IF ( xmu(2,0,0) .GT. 1.e-7 ) THEN

        skew(jsea)=xlambda(3,0,0)

        delta = ( xlambda(1,2,0) + xlambda(1,0,2)           &

                  - 2.0*xlambda(0,1,1)*xlambda(1,1,1) )/    &

                   (1.0 - xlambda(0,1,1)**2)             ! this is called gamma eq. 20

        embia1(jsea)=-0.125*delta                             ! EM Bias coefficient

        embia2(jsea)=-0.125*xlambda(3,0,0)/3.0           ! tracker bias (least squares only)

      END IF

    END DO  ! end of loop on JSEA

        !

#ifdef W3_OMPG

        !$OMP END PARALLEL DO

#endif


    DEALLOCATE(fac0,fac1,fac2,fac3)

    DEALLOCATE(f2,sighf,dfimhf,fak)


      END SUBROUTINE skewness


END MODULE w3iogomd