#include "w3macros.h"
!/ ------------------------------------------------------------------- /
      PROGRAM W3OUNF
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           F. Ardhuin              |
!/                  |           M. Accensi              |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         28-Mar-2019 |
!/                  +-----------------------------------+
!/
!/    17-Mar-2010 : Creation                            ( version 3.14_SHOM )
!/    07-Nov-2011 : Debug for spectral output on UNST   ( version 4.04 )
!/    13-Mar-2012 : Update of NC attributes             ( version 4.04 )
!/    02-Apr-2013 : New structure of output fields.     ( version 4.10 )
!/    02-Jul-2013 : Bug correction for lat in unst grid ( version 4.11 )
!/    02-Nov-2013 : Removes unnecessary IDFM            ( version 4.12 )
!/    30-Apr-2014 : Correct group3 freq dim.            ( version 5.00 )
!/    23-May-2014 : Adding ice fluxes to W3SRCE         ( version 5.01 )
!/    14-Oct-2014 : Keep the output files opened        ( version 5.01 )
!/    27-Aug-2015 : ICEH and ICEF added as output       ( version 5.10 )
!/    10-Jan-2017 : Changes for US3D and USSP output    ( version 6.01 )
!/    01-May-2017 : Adds directional MSS parameters     ( version 6.04 )
!/    01-Mar-2018 : RTD option add variable de-rotation,( version 6.02 )
!/                  standard lat-lons and rotated grid
!/                  metadata
!/    15-May-2018 : Add namelist feature                ( version 6.05 )
!/    06-Jun-2018 : Add DEBUG/SETUP                     ( version 6.04 )
!/    27-Jun-2018 : Updated to handle SMC output.       ( version 6.05 )
!/    26-Jul-2018 : Changed reading of TABIPART         ( version 6.05 )
!/    12-Sep-2018 : Added extra partitioned fields      ( version 6.06 )
!/    25-Sep-2018 : Add WBT parameter                   ( version 6.06 )
!/    28-Mar-2019 : Bugfix to NBIPART check.            ( version 6.07 )
!/    21-May-2019 : Feature fixed filename              ( version 7.01 )
!/
!/    Copyright 2009-2013 National Weather Service (NWS),
!/       National Oceanic and Atmospheric Administration.  All rights
!/       reserved.  WAVEWATCH III is a trademark of the NWS.
!/       No unauthorized use without permission.
!/
!  1. Purpose :
!
!     Post-processing of grid output to NetCDF files
!
!  2. Method :
!
!     Data is read from the grid output file out_grd.ww3 (raw data)
!     and from the file ww3_ounf.nml or ww3_ounf.inp (NDSI)
!     Model definition and raw data files are read using WAVEWATCH III
!     subroutines. Extra global NetCDF attributes may be read from
!     ASCII file NC_globatt.inp.
!
!     Output types :
!      4 : NetCDF files
!
!  3. Parameters :
!
!  4. Subroutines used :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!      W3NMOD    Subr. W3GDATMD Set number of model.
!      W3SETG    Subr.   Id.    Point to selected model.
!      W3NDAT    Subr. W3WDATMD Set number of model for wave data.
!      W3SETW    Subr.   Id.    Point to selected model for wave data.
!      W2NAUX    Subr. W3ADATMD Set number of model for aux data.
!      W3SETA    Subr.   Id.    Point to selected model for aux data.
!      ITRACE    Subr. W3SERVMD Subroutine tracing initialization.
!      STRACE    Subr.   Id.    Subroutine tracing.
!      NEXTLN    Subr.   Id.    Get next line from input filw
!      EXTCDE    Subr.   Id.    Abort program as graceful as possible.
!      STME21    Subr. W3TIMEMD Convert time to string.
!      TICK21    Subr.   Id.    Advance time.
!      DSEC21    Func.   Id.    Difference between times.
!      W3IOGR    Subr. W3IOGRMD Reading/writing model definition file.
!      W3IOGO    Subr. W3IOGOMD Reading/writing raw gridded data file.
!      W3EXNC    Subr. Internal Execute grid netcdf output.
!     ----------------------------------------------------------------
!
!  5. Called by :
!
!     None, stand-alone program.
!
!  6. Error messages :
!
!     Checks on input, checks in W3IOxx.
!
!  7. Remarks :
!
!  8. Structure :
!
!     See source code.
!
!  9. Switches :
!
!     !/S     Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE CONSTANTS
 
!/
      USE W3WDATMD, ONLY: W3NDAT, W3SETW
      USE W3ADATMD, ONLY: W3NAUX, W3SETA
      USE W3ODATMD, ONLY: W3NOUT, W3SETO
      USE W3SERVMD, ONLY : ITRACE, NEXTLN, EXTCDE
      USE W3TIMEMD
      USE W3IOGRMD, ONLY: W3IOGR
      USE W3IOGOMD, ONLY: W3IOGO, W3READFLGRD, W3FLGRDFLAG
      USE W3INITMD, ONLY: WWVER, SWITCHES
      USE W3ODATMD, ONLY: NAPROC, NOSWLL, PTMETH, PTFCUT
!/
      USE W3GDATMD
      USE W3WDATMD, ONLY: TIME, WLV, ICE, ICEH, ICEF, BERG, UST, USTDIR
      USE W3ADATMD, ONLY: DW, UA, UD, AS, CX, CY, HS, WLM, T0M1, THM,  &
                          THS, FP0, THP0, DTDYN, FCUT,                 &
                          ABA, ABD, UBA, UBD, SXX, SYY, SXY, USERO,    &
                          PHS, PTP, PLP, PDIR, PSI, PWS, PWST, PNR,    &
                          PTM1, PT1, PT2, PEP,                         &
                          PTHP0, PQP, PSW, PPE, PGW, QP,               &
                          TAUOX, TAUOY, TAUWIX,                        &
                          TAUWIY, PHIAW, PHIOC, TUSX, TUSY, PRMS, TPMS,&
                          USSX, USSY, MSSX, MSSY, MSSD, MSCX, MSCY,    &
                          MSCD, CHARN, TWS,                            &
                          TAUWNX, TAUWNY, BHD, T02, HSIG, CGE,         &
                          T01, BEDFORMS, WHITECAP, TAUBBL, PHIBBL,     &
                          CFLTHMAX, CFLXYMAX, CFLKMAX, TAUICE, PHICE,  &
                          STMAXE, STMAXD, HMAXE, HCMAXE, HMAXD, HCMAXD,&
                          P2SMS, EF, US3D, TH1M, STH1M, TH2M, STH2M,   &
                          WN, USSP, WBT
      USE W3ODATMD, ONLY: NDSO, NDSE, SCREEN, NOGRP, NGRPP, IDOUT,     &
                          UNDEF, FLOGRD, FNMPRE, NOSWLL, NOGE
!
      USE W3NMLOUNFMD
      USE NETCDF
 
 
      IMPLICIT NONE
 
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      TYPE(NML_FIELD_T)       :: NML_FIELD
      TYPE(NML_FILE_T)        :: NML_FILE
      TYPE(NML_SMC_T)         :: NML_SMC
!
      INTEGER                 :: NDSI, NDSM, NDSOG,                    &
                                 NDSTRC, NTRACE, IERR, I, I1F, I2F,    &
                                 IOTEST, NOUT,                         &
                                 IFI, IFJ, NCTYPE, IX1, IXN, IY1, IYN, &
                                 IOUT, S3, IRET, HASNC4,               &
                                 NBIPART, CNTIPART, NCVARTYPE, IPART,  &
                                 RTDNX, RTDNY
      INTEGER                 :: TOUT(2), TDUM(2), STOPDATE(8)
!
      INTEGER, ALLOCATABLE    :: TABIPART(:), NCIDS(:,:,:)
!
      REAL                    :: DTREQ, DTEST
!
      CHARACTER*30            :: STRSTOPDATE, FILEPREFIX, STRINGIPART
      CHARACTER*1024          :: FLDOUT
      CHARACTER               :: COMSTR*1, IDTIME*23, IDDDAY*11
!
      LOGICAL                 :: FLG2D(NOGRP,NGRPP), FLG1D(NOGRP),     &
                                 VECTOR, TOGETHER, FLGNML
 
      LOGICAL                 :: SMCGRD = .FALSE.
!/
!/ ------------------------------------------------------------------- /
!/
! 1.  IO set-up.
!
      CALL W3NMOD ( 1, 6, 6 )
      CALL W3SETG ( 1, 6, 6 )
      CALL W3NDAT (    6, 6 )
      CALL W3SETW ( 1, 6, 6 )
      CALL W3NAUX (    6, 6 )
      CALL W3SETA ( 1, 6, 6 )
      CALL W3NOUT (    6, 6 )
      CALL W3SETO ( 1, 6, 6 )
!
      NDSI   = 10
      NDSM   = 20
      NDSOG  = 20
!
      NDSTRC =  6
      NTRACE = 10
      CALL ITRACE ( NDSTRC, NTRACE )
!
      WRITE (NDSO,900)
 
!
!--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! 2.  Read model definition file.
!
      CALL W3IOGR ( 'READ', NDSM )
      WRITE (NDSO,920) GNAME
!
!--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! 3.  Read general data and first fields from file
!
      CALL W3IOGO ( 'READ', NDSOG, IOTEST )
!
      WRITE (NDSO,930)
      DO IFI=1, NOGRP
        DO IFJ=1, NGRPP
          IF ( FLOGRD(IFI,IFJ) ) WRITE (NDSO,931) IDOUT(IFI,IFJ)
        END DO
      END DO
!
!--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! 4.  Read requests from input file.
!
 
!
! process ww3_ounf namelist
!
      INQUIRE(FILE=TRIM(FNMPRE)//"ww3_ounf.nml", EXIST=FLGNML)
      IF (FLGNML) THEN
        ! Read namelist
        CALL W3NMLOUNF (NDSI, TRIM(FNMPRE)//'ww3_ounf.nml', NML_FIELD, &
                        NML_FILE, NML_SMC, IERR)
 
! 4.1 Time setup
        READ(NML_FIELD%TIMESTRIDE, *)  DTREQ
        READ(NML_FIELD%TIMECOUNT, *)   NOUT
        READ(NML_FIELD%TIMESTART, *)   TOUT(1), TOUT(2)
 
! 4.2 Output fields
        FLDOUT = NML_FIELD%LIST
        CALL W3FLGRDFLAG ( NDSO, SCREEN, NDSE, FLDOUT, FLG1D,       &
                           FLG2D, 1, 1, IERR )
        IF (IERR.NE.0) GOTO 800
 
! 4.3 Output type
        NCTYPE = NML_FILE%NETCDF
        NCVARTYPE = NML_FIELD%TYPE
        STRINGIPART = NML_FIELD%PARTITION
        TOGETHER = NML_FIELD%SAMEFILE
        FILEPREFIX = NML_FILE%PREFIX
        S3 = NML_FIELD%TIMESPLIT
        IF(SMCGRD) THEN
        ELSE
          IX1 = NML_FILE%IX0
          IXN = NML_FILE%IXN
          IY1 = NML_FILE%IY0
          IYN = NML_FILE%IYN
        ENDIF ! SMCGRD
      END IF ! FLGNML
 
!
! process old ww3_ounf.inp format
!
      IF (.NOT. FLGNML) THEN
        OPEN (NDSI,FILE=TRIM(FNMPRE)//'ww3_ounf.inp',STATUS='OLD',ERR=800,IOSTAT=IERR)
        REWIND (NDSI)
 
        READ (NDSI,'(A)',END=801,ERR=802,IOSTAT=IERR) COMSTR
        IF (COMSTR.EQ.' ') COMSTR = '$'
        WRITE (NDSO,901) COMSTR
        CALL NEXTLN ( COMSTR , NDSI , NDSE )
 
! 4.1 Time setup
        READ (NDSI,*,END=801,ERR=802) TOUT, DTREQ, NOUT
 
! 4.2 Output fields
        CALL W3READFLGRD ( NDSI, NDSO, SCREEN, NDSE, COMSTR, FLG1D,      &
                           FLG2D, 1, 1, IERR )
        IF (IERR.NE.0) GOTO 800
 
! 4.3 Output type
        CALL NEXTLN ( COMSTR , NDSI , NDSE )
        READ (NDSI,*,END=801,ERR=802) NCTYPE, NCVARTYPE
        CALL NEXTLN ( COMSTR , NDSI , NDSE )
        READ (NDSI,'(A)',END=801,ERR=802) STRINGIPART
        CALL NEXTLN ( COMSTR , NDSI , NDSE )
        READ (NDSI,*,END=801,ERR=802) TOGETHER
        CALL NEXTLN ( COMSTR , NDSI , NDSE )
        FILEPREFIX= 'ww3.'
        READ (NDSI,*,END=801,ERR=802) FILEPREFIX
        CALL NEXTLN ( COMSTR , NDSI , NDSE )
        READ (NDSI,*,END=801,ERR=802) S3
        CALL NEXTLN ( COMSTR , NDSI , NDSE )
 
        IF( .NOT. SMCGRD) THEN
           READ (NDSI,*,END=801,ERR=802) IX1, IXN, IY1, IYN
        ENDIF
 
        CLOSE(NDSI,ERR=800,IOSTAT=IERR)
 
      END IF ! .NOT. FLGNML
 
!
 
! 4.1 Time setup
      DTREQ  = MAX ( 0. , DTREQ )
      IF ( DTREQ.EQ.0. ) NOUT = 1
      NOUT   = MAX ( 1 , NOUT )
      CALL STME21 ( TOUT , IDTIME )
      WRITE (NDSO,940) IDTIME
      TDUM = 0
      CALL TICK21 ( TDUM , DTREQ )
      CALL STME21 ( TDUM , IDTIME )
      IF ( DTREQ .GE. 86400. ) THEN
        WRITE (IDDDAY,'(I10,1X)') INT(DTREQ/86400.)
      ELSE
        IDDDAY = '           '
      END IF
      IDTIME(1:11) = IDDDAY
      IDTIME(21:23) = '   '
      WRITE (NDSO,941) IDTIME, NOUT
 
 
! 4.2 Output fields
      DO IFI=1, NOGRP
        DO IFJ=1, NGRPP
          IF ( FLG2D(IFI,IFJ) ) THEN
            IF ( FLOGRD(IFI,IFJ) ) THEN
              WRITE (NDSO,946) IDOUT(IFI,IFJ), ' '
            ELSE
              WRITE (NDSO,946) IDOUT(IFI,IFJ), '*** NOT AVAILABLE ***'
              FLG2D(IFI,IFJ) = .FALSE.
            END IF
          END IF
        END DO
      END DO
 
 
! 4.3 Output type
!!      NBIPART=0
!!      DO I=1,29
!!        IF ((STRINGIPART(I:I+1).EQ.'0').OR.(STRINGIPART(I:I+1).EQ.'1')      &
!!            .OR.(STRINGIPART(I:I+1).EQ.'2').OR.(STRINGIPART(I:I+1).EQ.'3')  &
!!            .OR.(STRINGIPART(I:I+1).EQ.'4').OR.(STRINGIPART(I:I+1).EQ.'5')) THEN
!!          NBIPART=NBIPART+1
!!        END IF
!!      END DO
!!      ALLOCATE(TABIPART(NBIPART))
!!      CNTIPART=0
!!      DO I=1,29
!!        IF ((STRINGIPART(I:I+1).EQ.'0').OR.(STRINGIPART(I:I+1).EQ.'1')      &
!!            .OR.(STRINGIPART(I:I+1).EQ.'2').OR.(STRINGIPART(I:I+1).EQ.'3')  &
!!            .OR.(STRINGIPART(I:I+1).EQ.'4').OR.(STRINGIPART(I:I+1).EQ.'5')) THEN
!!          CNTIPART=CNTIPART+1
!!          READ(STRINGIPART(I:I+1),'(I1)') TABIPART(CNTIPART)
!!        END IF
!!      END DO
 
      ! Alternative processing of TABIPART to capture requests
      ! greater than NOSWLL (C.Bunney):
      ALLOCATE(TABIPART(NOSWLL + 1))
      ALLOCATE(NCIDS(NOGRP,NGRPP,NOSWLL + 1))
      NBIPART=0
      DO I=1,30
        IF(STRINGIPART(I:I) .EQ. ' ') CYCLE
        READ(STRINGIPART(I:I),'(I1)') IPART
        IF(IPART .GT. NOSWLL) THEN
           WRITE(NDSO, 1500) IPART, NOSWLL
           CYCLE
        ENDIF
        NBIPART = NBIPART + 1
        IF(NBIPART .GT. NOSWLL + 1) THEN
           GOTO 803
        ENDIF
        TABIPART(NBIPART) = IPART
      ENDDO
!
      IF ( NCTYPE.LT.3 .OR. NCTYPE.GT.4 ) THEN
        WRITE (NDSE,1010) NCTYPE
        CALL EXTCDE ( 1 )
      END IF
      ! if NCTYPE = 4 checking that it is compiled with NC4
      HASNC4=0
      HASNC4=1
      IF ((HASNC4 .EQ. 0).AND.(NCTYPE.EQ.4)) THEN
        WRITE (NDSE,1012)
        CALL EXTCDE ( 1 )
      END IF
 
 
      IF(.NOT. SMCGRD) THEN
        IX1    = MAX ( IX1 , 1 )
        IXN    = MIN ( IXN , NX )
        IY1    = MAX ( IY1 , 1 )
        IYN    = MIN ( IYN , NY )
        WRITE (NDSO,3940) IX1, IXN, IY1, IYN
      ENDIF ! SMCGRD
 
      VECTOR = .TRUE.
!
!--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! 5.  Time management.
!
      IOUT = 0
      NCIDS(:,:,:) = 0
      WRITE (NDSO,970)
 
 
! 5.1 Loops on out_grd.ww3 to read the time and data
      DO
        DTEST  = DSEC21 ( TIME , TOUT )
        IF ( DTEST .GT. 0. ) THEN
          CALL W3IOGO ( 'READ', NDSOG, IOTEST )
            IF ( IOTEST .EQ. -1 ) THEN
              WRITE (NDSO,944)
              EXIT
            END IF
          CYCLE
        END IF
        IF ( DTEST .LT. 0. ) THEN
          CALL TICK21 ( TOUT , DTREQ )
          CYCLE
        END IF
 
 
! 5.1.1 Increments the time counter IOUT
        IOUT   = IOUT + 1
        CALL STME21 ( TOUT , IDTIME )
        WRITE (NDSO,971) IDTIME
 
 
! 5.1.2  Processes the variable value for the time step IOUT
        CALL W3EXNC ( NX, NY, IX1, IXN, IY1, IYN, NSEA, FILEPREFIX,   &
                      E3DF, P2MSF, US3DF, USSPF, NCTYPE, TOGETHER, NCVARTYPE,&
                      FLG2D, NCIDS, S3, STRSTOPDATE )
 
! 5.1.3 Defines the stop date
        CALL T2D(TOUT,STOPDATE,IERR)
        WRITE(STRSTOPDATE,'(I4.4,A,4(I2.2,A),I2.2)') STOPDATE(1),'-',STOPDATE(2), &
              '-',STOPDATE(3),' ',STOPDATE(5),':',STOPDATE(6),':',STOPDATE(7)
 
        CALL TICK21 ( TOUT , DTREQ )
        IF ( IOUT .GE. NOUT ) EXIT
      END DO
 
 
! 5.2 Closes the netCDF file
      IF (TOGETHER .AND. NCIDS(1,1,1).NE.0) THEN
        IRET = NF90_REDEF(NCIDS(1,1,1))
        CALL CHECK_ERR(IRET)
        IRET=NF90_PUT_ATT(NCIDS(1,1,1),NF90_GLOBAL,'stop_date',STRSTOPDATE)
        CALL CHECK_ERR(IRET)
        IRET=NF90_CLOSE(NCIDS(1,1,1))
        CALL CHECK_ERR(IRET)
      END IF
!
      DO IFI=1, NOGRP
        DO IFJ=1, NGRPP
          IF ( FLG2D(IFI,IFJ) ) THEN
            IF ( FLOGRD(IFI,IFJ) ) THEN
              IF (.NOT. TOGETHER) THEN
                IF (NCIDS(IFI,IFJ,1).NE.0) THEN
                  IRET = NF90_REDEF(NCIDS(IFI,IFJ,1))
                  CALL CHECK_ERR(IRET)
                  IRET=NF90_PUT_ATT(NCIDS(IFI,IFJ,1),NF90_GLOBAL,'stop_date',STRSTOPDATE)
                  CALL CHECK_ERR(IRET)
                  IRET=NF90_CLOSE(NCIDS(IFI,IFJ,1))
                  CALL CHECK_ERR(IRET)
                END IF ! NCIDS
                ! close partition files (except part 0 which is already closed by (IFI,IFJ,1)
                IF ((IFI.EQ.4).AND.(IFJ.LE.NOGE(IFI))) THEN
                  DO IPART=1,NOSWLL
                    IF (NCIDS(IFI,IFJ,IPART+1).NE.0) THEN
                      IRET = NF90_REDEF(NCIDS(IFI,IFJ,IPART+1))
                      CALL CHECK_ERR(IRET)
                      IRET=NF90_PUT_ATT(NCIDS(IFI,IFJ,IPART+1),NF90_GLOBAL,'stop_date',STRSTOPDATE)
                      CALL CHECK_ERR(IRET)
 
                      IRET=NF90_CLOSE(NCIDS(IFI,IFJ,IPART+1))
                      CALL CHECK_ERR(IRET)
                    END IF ! NCIDS
                  END DO ! IPART
                END IF ! partition
              ! else if together
              ELSE
                ! close frequency file
                IF ( ((IFI.EQ.6).AND.(IFJ.EQ.8)) .OR.                 &
                     ((IFI.EQ.6).AND.(IFJ.EQ.9)) .OR.                 &
                     (IFI.EQ.3) ) THEN
                  IF (NCIDS(IFI,IFJ,1).NE.0) THEN
                    IRET = NF90_REDEF(NCIDS(IFI,IFJ,1))
                    CALL CHECK_ERR(IRET)
                    IRET=NF90_PUT_ATT(NCIDS(IFI,IFJ,1),NF90_GLOBAL,'stop_date',STRSTOPDATE)
                    CALL CHECK_ERR(IRET)
                    IRET=NF90_CLOSE(NCIDS(IFI,IFJ,1))
                    CALL CHECK_ERR(IRET)
                  END IF ! NCIDS
                END IF ! IFI
              END IF ! TOGETHER
            END IF ! FLOGRD
          END IF ! FLG2D
        END DO ! IFJ
      END DO ! IFI
 
!
      GOTO 888
!
! Escape locations read errors :
!
  800 CONTINUE
      WRITE (NDSE,1000) IERR
      CALL EXTCDE ( 10 )
!
  801 CONTINUE
      WRITE (NDSE,1001)
      CALL EXTCDE ( 11 )
!
  802 CONTINUE
      WRITE (NDSE,1002) IERR
      CALL EXTCDE ( 12 )
!
  803 CONTINUE
      WRITE (NDSE,1003) NBIPART, NOSWLL
      CALL EXTCDE (13)
!
  888 CONTINUE
      WRITE (NDSO,999)
!
! Formats
!
  900 FORMAT (/15X,'   *** WAVEWATCH III Field output postp. ***   '/ &
               15X,'==============================================='/)
  901 FORMAT ( '  Comment character is ''',A,''''/)
!
  920 FORMAT ( '  Grid name : ',A/)
!
  930 FORMAT ( '  Fields in file : '/                                 &
               ' --------------------------')
  931 FORMAT ( '      ',A)
!
  940 FORMAT (/'  Output time data : '/                               &
               ' --------------------------------------------------'/ &
               '      First time         : ',A)
  941 FORMAT ( '      Interval           : ',A/                       &
               '      Number of requests : ',I10)
  944 FORMAT (/'      End of file reached '/)
  946 FORMAT ( '      ',A,2X,A)
!
 3940 FORMAT ( '      X range : ',2I7/                                &
               '      Y range : ',2I7)
!
  970 FORMAT (/'  Generating files '/                                 &
               ' --------------------------------------------------')
  971 FORMAT ( '      Files for ',A)
!
  999 FORMAT (/'  End of program '/                                   &
               ' ========================================='/          &
               '         WAVEWATCH III Field output '/)
!
 1000 FORMAT (/' *** WAVEWATCH III ERROR IN W3OUNF : '/               &
               '     ERROR IN OPENING INPUT FILE'/                    &
               '     IOSTAT =',I5/)
!
 1001 FORMAT (/' *** WAVEWATCH III ERROR IN W3OUNF : '/               &
               '     PREMATURE END OF INPUT FILE'/)
!
 1002 FORMAT (/' *** WAVEWATCH III ERROR IN W3OUNF : '/               &
               '     ERROR IN READING FROM INPUT FILE'/               &
               '     IOSTAT =',I5/)
!
 1003 FORMAT (/' *** WAVEWATCH III WERROR IN W3EXNC : '/              &
               '     OUT OF RANGE REQUEST FOR NBIPART =',I2, /        &
               '     MAX SWELL PARTITIONS (NOSW) =',I2 /)
!
 1010 FORMAT (/' *** WAVEWATCH III ERROR IN W3OUNF : '/               &
               '     ILLEGAL TYPE, NCTYPE =',I4/)
 1012 FORMAT (/' *** WAVEWATCH III ERROR IN W3OUNF : '/               &
               '     NCTYPE = 4 BUT NOT COMPILED WITH NC4'/)
!
 1500 FORMAT (/' *** WAVEWATCH III WARNING IN W3EXNC : '/             &
               '     IGNORING REQUEST FOR IPART =',I2, /              &
               '     MAX SWELL PARTITIONS (NOSW) =',I2 /)
!
!/
!/ Internal subroutine W3EXNC ---------------------------------------- /
!/
      CONTAINS
!/ ------------------------------------------------------------------- /
      SUBROUTINE W3EXNC ( NX, NY, IX1, IXN, IY1, IYN, NSEA,             &
                          FILEPREFIX, E3DF, P2MSF, US3DF, USSPF,NCTYPE, &
                          TOGETHER, NCVARTYPEI, FLG2D, NCIDS, S3, STRSTOPDATE )
!/
!/                  +-----------------------------------+
!/                  |           F. Ardhuin              |
!/                  |           M. Accensi              |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         14-Oct-2014 |
!/                  +-----------------------------------+
!/
!/    17-Mar-2010 : Creation                            ( version 3.14_SHOM )
!/    28-Feb-2013 : New option for float output         ( version 4.08 )
!/    02-Apr-2013 : New structure of output fields.     ( version 4.09 )
!/    12-Apr-2013 : Allows curvilinear grids            ( version 4.10 )
!/    30-Apr-2014 : Correct group3 freq dim.            ( version 5.00 )
!/    23-May-2014 : Adding ice fluxes to W3SRCE         ( version 5.01 )
!/    14-Oct-2014 : Keep the output files opened        ( version 5.01 )
!/
!  1. Purpose :
!
!     Perform actual grid output in NetCDF file.
!
!  3. Parameters :
!
!     Parameter list
!     ----------------------------------------------------------------
!       NX/Y    Int.  I  Grid dimensions.
!       IX1/IXN Int.  I  Grid indexes along X
!       IY1/IYN Int.  I  Grid indexes along Y
!       NSEA    Int.  I  Number of sea points.
!     ----------------------------------------------------------------
!
!     Internal parameters
!     ----------------------------------------------------------------
!       FLTWO   Log.  Flags for two-dimensional field X Y.
!       FLDIR   Log.  Flags for two-dimensional, directional field.
!       FLFRQ   Log.  Flags for frequency array (3D field)
!       X1, X2, XX, XY
!               R.A.  Output fields
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!      STRACE    Subr. W3SERVMD Subroutine tracing.
!      EXTCDE    Subr.   Id.    Abort program as graceful as possible.
!      W3S2XY    Subr.   Id.    Convert from storage to spatial grid.
!      PRTBLK    Subr. W3ARRYMD Print plot of array.
!      OUTA2I    Subr.   Id.    Print array of INTEGERS.
!     ----------------------------------------------------------------
!
!  5. Called by :
!
!     Main program in which it is contained.
!
!  6. Error messages :
!
!       None.
!
!  7. Remarks :
!
!     - Note that arrays CX and CY of the main program now contain
!       the absolute current speed and direction respectively.
!
!  8. Structure :
!
!     See source code.
!
!  9. Switches :
!
!       !/S  Enable subroutine tracing.
!       !/T  Enable test output.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE W3SERVMD, ONLY : W3S2XY
      USE W3ARRYMD, ONLY : OUTA2I, PRTBLK
      USE W3GDATMD, ONLY : SIG, GTYPE, FLAGLL, MAPSTA, MAPST2
      USE NETCDF
      IMPLICIT NONE
 
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
      INTEGER, INTENT(IN)     :: NX, NY, IX1, IXN, IY1, IYN, NSEA,     &
                                 E3DF(3,5), P2MSF(3), US3DF(3),        &
                                 USSPF(2), NCTYPE, NCVARTYPEI
      CHARACTER(30)           :: FILEPREFIX
      LOGICAL, INTENT(IN)     :: TOGETHER
      LOGICAL, INTENT(IN)     :: FLG2D(NOGRP,NGRPP)
      INTEGER, INTENT(INOUT) :: NCIDS(NOGRP,NGRPP,NOSWLL + 1), S3
      CHARACTER*30,INTENT(IN) :: STRSTOPDATE
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      INTEGER                 :: IFI, IFJ, MFILL, I, J, ISEA, IX, IY,  &
                                 I1, J1, IPART, INDEXIPART, COORDTYPE
      INTEGER                 :: S1, S2, S4, S5, NCID, OLDNCID,        &
                                 MAPSTAOUT, NDSDAT, VMIN, VMAX,        &
                                 NFIELD, N, IRET, IK, EXTRADIM, IVAR,  &
                                 IVAR1, IDVAROUT, NCVARTYPE
      INTEGER                 :: DIMID(6), VARID(300), START(4),       &
                                 COUNT(4), DIMLN(6),START1D(2),        &
                                 COUNT1D(2), DIMFIELD(3),              &
                                 STARTDATE(8), CURDATE(8), REFDATE(8), &
                                 MAP(NX+1,NY), MP2(NX+1,NY)
!
    INTEGER                  :: DEFLATE=1
!
      INTEGER, ALLOCATABLE    :: TRIGP2(:,:)
      ! Make the below allocatable to avoid stack overflow on some machines
      INTEGER(KIND=2), ALLOCATABLE    :: MX1(:,:), MXX(:,:), MYY(:,:), &
                                         MXY(:,:), MAPOUT(:,:)
!
      REAL                    :: FSC, FSC3, CABS, UABS, MFILLR
   REAL, PARAMETER            :: LOG2=LOG(2.)
!
      REAL,DIMENSION(:),  ALLOCATABLE    :: LON, LAT, FREQ
      REAL,DIMENSION(:,:),  ALLOCATABLE  :: LON2D, LAT2D, ANGLD2D
      ! Make the below allocatable to avoid stack overflow on some machines
      REAL, ALLOCATABLE       :: X1(:,:), X2(:,:), XX(:,:), XY(:,:),   &
                                 XK(:,:,:), XXK(:,:,:), XYK(:,:,:),    &
                                 MX1R(:,:), MXXR(:,:), MYYR(:,:),      &
                                 MXYR(:,:)
!
      DOUBLE PRECISION        :: OUTJULDAY
      DOUBLE PRECISION        :: SXD, SYD, X0D, Y0D
!
      CHARACTER*80            :: VARNM(8),VARNL(8)
      CHARACTER*120           :: VARNS(8),VARNG(8), VARND(8), STR2
      CHARACTER*512           :: VARNC(8), PARTCOM
      CHARACTER*30            :: UNITVAR(3),FORMAT1
      CHARACTER*30            :: STRSTARTDATE
      CHARACTER               :: FNAMENC*50,                           &
                                 ENAME*10,FORMF*11, UNITS*24, UNITS2*24
      CHARACTER, SAVE         :: OLDTIMEID*16 = '0000000000000000'
      CHARACTER, SAVE         :: TIMEID*16 = '0000000000000000'
!
      LOGICAL                 :: FLFRQ, FLDIR, FEXIST, FREMOVE
      LOGICAL                 :: CUSTOMFRQ=.FALSE.
 
!/
!/ ------------------------------------------------------------------- /
!/
!
!--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! 1.  Preparations
!
      ! Allocate output storage. This is required with the introduction
      ! of the SMC grid output as the regridded output grid dimensions could
      ! conceivably be larger than the NX and NY values. Making these (large)
      ! arrays allocatable also moves them to the heap and avoids stack
      ! overflow issues that can occur on some architectures. (Chris Bunney)
      IF(SMCGRD) THEN
      ELSE
        ALLOCATE(X1(NX+1,NY),X2(NX+1,NY),XX(NX+1,NY),XY(NX+1,NY))
        ALLOCATE(XK(NX+1,NY,NK), XXK(NX+1,NY,NK), XYK(NX+1,NY,NK))
        ALLOCATE(MX1(NX,NY), MXX(NX,NY), MYY(NX,NY), MXY(NX,NY), MAPOUT(NX,NY))
        ALLOCATE(MX1R(NX,NY), MXXR(NX,NY), MYYR(NX,NY), MXYR(NX,NY))
      ENDIF ! SMCGRD
 
      X1     = UNDEF
      X2     = UNDEF
      XX     = UNDEF
      XY     = UNDEF
      ! CB: Dont output MAPSTA for SMC grid - it does not make sense
      IF( SMCGRD ) THEN
        MAPSTAOUT = 0
      ELSE
        MAPSTAOUT = 1
      ENDIF
      NCVARTYPE  = NCVARTYPEI
      NDSDAT=30
      NCID = 0
!
      CALL U2D('days since 1990-01-01 00:00:00',REFDATE,IERR)
 
! 1.1 Set-up transfer files
      MFILL  = NF90_FILL_SHORT
      MFILLR  = NF90_FILL_FLOAT
      IF (GTYPE.NE.UNGTYPE) THEN
        COORDTYPE=1
      ELSE
        COORDTYPE=2
      ENDIF
 
! 1.2 Sets the date as ISO8601 convention
      ! S3 defines the number of characters in the date for the filename
      ! S3=0 -> field, S3=4-> YYYY, S3=6 -> YYYYMM, S3=10 -> YYYYMMDDHH
      ! Setups min and max date format
      IF (S3.GT.0 .AND. S3.LT.4) S3=4
      IF (S3.GT.10) S3=10
!
      ! Defines the format of FILETIME
      S5=S3-8
      S4=S3
      OLDTIMEID=TIMEID
      ! if S3=>nodate then filetime='field'
      IF (S3.EQ.0) THEN
        S4=5
        TIMEID="field"
      ! if S3=>YYYYMMDDHH then filetime='YYYYMMDDTHHZ'
      ELSE IF (S3.EQ.10) THEN
        S4=S4+2 ! add chars for ISO8601 : day T hours Z
        WRITE(FORMAT1,'(A,I1,A,I1,A)') '(I8.8,A1,I',S5,'.',S5,',A1)'
        WRITE (TIMEID,FORMAT1) TIME(1), 'T', &
               FLOOR(REAL(TIME(2))/NINT(10.**(6-S5))), 'Z'
      ! if S3=>YYYYMMDD then filetime='YYYYMMDD'
      ELSE IF (S3.EQ.8) THEN
        WRITE(FORMAT1,'(A,I1,A,I1,A)') '(I',S3,'.',S3,')'
        WRITE (TIMEID,FORMAT1) TIME(1)
      ! if S3=>YYYYMM then filetime='YYYYMM'
      ! or S3=>YYYY then filetime='YYYY'
      ELSE
        WRITE(FORMAT1,'(A,I1,A,I1,A)') '(I',S3,'.',S3,')'
        WRITE (TIMEID,FORMAT1) FLOOR(REAL(TIME(1))/NINT(10.**(8-S3)))
      END IF
      ! redefines filename with updated date format
      S1=LEN_TRIM(FILEPREFIX)
      FNAMENC=''
      FNAMENC(1:S1)=FILEPREFIX(1:S1)
      FNAMENC(S1+1:S1+S4) = TIMEID(1:S4)
 
      ! Create partitioning method comment string:
      IF ( PTMETH .EQ. 1 ) THEN
        PARTCOM = "Wind sea and swells defined using topographic " //  &
                  "partitions and partition wave-age cut-off "     //  &
                  "(WWIII default scheme)"
      ELSE IF ( PTMETH .EQ. 2 ) THEN
        PARTCOM = "Wind sea and swells defined using topographic " //  &
                  "partitions and spectral wave-age cut-off"
      ELSE IF ( PTMETH .EQ. 3 ) THEN
        PARTCOM = "Wave components defined using topographic "     //  &
                  "partitions only"
      ELSE IF ( PTMETH .EQ. 4 ) THEN
        PARTCOM = "Wind sea and swell defined using spectral wave-age cut-off"
      ELSE IF ( PTMETH .EQ. 5 ) THEN
        WRITE(PARTCOM, '("Wave components defined using ",F5.3,'   //  &
                  '"Hz spectral frequency cutoff")') PTFCUT
      ELSE
        WRITE(PARTCOM, '("PTM_",I1,"_Unknown")') PTMETH
      ENDIF
 
!
!--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! 2.  Loop over output fields.
!
 
      ! Instanciates the field and group indexes
      I1=0
      J1=0
!
      DO IFI=1, NOGRP
        DO IFJ=1, NGRPP
          ENAME='      '
          ! If the flag for the variable IFI of the group IFJ is .TRUE.
          IF ( FLG2D(IFI,IFJ) ) THEN
            ! Instanciates the partition array
            INDEXIPART=1
            IPART=TABIPART(INDEXIPART)
            VARNS(:)=''
 
 
!  Loop over IPART for partition variables
555         CONTINUE
 
            ! Initializes the index of field and group at the first flag FLG2D at .TRUE.
            IF (I1.EQ.0) I1=IFI
            IF (J1.EQ.0) J1=IFJ
            FORMF  = '(1X,32I5)'
!
! 2.1 Set output arrays and parameters
!
            ! Initializes the flags for freq and direction dimensions
            FLFRQ = .FALSE.
            FLDIR = .FALSE.
            UNITS2='SAME'
            IF (NCVARTYPEI.EQ.3) NCVARTYPE=2
!
            IF ( IFI .EQ. 1 .AND. IFJ .EQ. 1 ) THEN
              FSC    = 0.5
              UNITS  = 'm'
              NFIELD = 1
              ENAME  = '.dpt'
              IF( SMCGRD ) THEN
              ELSE ! IF(SMCGRD)
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, DW(1:NSEA)     &
                                                            , MAPSF, X1 )
              ENDIF
              VARNM(1)='dpt'
              VARNL(1)='depth'
              VARNS(1)='depth'
              VARNG(1)='depth'
              VARNC(1)=''
              VARND(1)=''
              VMIN = -90000
              VMAX = 140000
!
            ELSE IF ( IFI .EQ. 1 .AND. IFJ .EQ. 2 ) THEN
              FSC    = 0.01
              ENAME  = '.cur'
              UNITS  = 'm s-1'
              !! Note - CX and CY read in from .ww3 file are X-Y vectors
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, CX(1:NSEA)        &
                                                        , MAPSF, XX )
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, CY(1:NSEA)        &
                                                        , MAPSF, XY )
              ENDIF
              !! Commented out unnecessary statements below for time being
              !! CX,CY are in north-east convention and X1,X2
              !! are not actually written out below
              !DO ISEA=1, NSEA
              !  CABS   = SQRT(CX(ISEA)**2+CY(ISEA)**2)
              !  IF ( CABS .GT. 0.05 ) THEN
              !    CY(ISEA) = MOD ( 630. -                         &
              !              RADE*ATAN2(CY(ISEA),CX(ISEA)) , 360. )
              !  ELSE
              !    CY(ISEA) = UNDEF
              !    END IF
              !  CX(ISEA) = CABS
              !  END DO
              !CALL W3S2XY ( NSEA, NSEA, NX+1, NY, CX(1:NSEA)        &
              !                                          , MAPSF, X1 )
              !CALL W3S2XY ( NSEA, NSEA, NX+1, NY, CY(1:NSEA)        &
              !                                          , MAPSF, X2 )
              NFIELD=2
              VARNM(1)='ucur'
              VARNM(2)='vcur'
              VARNL(1)='eastward current'
              VARNL(2)='northward current'
              VARNS(1)='eastward_sea_water_velocity'
              VARNS(2)='northward_sea_water_velocity'
              VARNG(1)='eastward_sea_water_velocity'
              VARNG(2)='northward_sea_water_velocity'
              VARNC(1)='cur=sqrt(U**2+V**2)'
              VARNC(2)='cur=sqrt(U**2+V**2)'
              VARND(1)=''
              VARND(2)=''
              VMIN = -990
              VMAX =  990
!
            ELSE IF ( IFI .EQ. 1 .AND. IFJ .EQ. 3 ) THEN
              FSC    = 0.1
              ENAME  = '.wnd'
              UNITS  = 'm s-1'
              !! Note - UA and UD read in from .ww3 file are UX,UY
 
              IF( SMCGRD ) THEN
              ELSE ! IF(SMCGRD)
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, UA(1:NSEA)        &
                                                        , MAPSF, XX )
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, UD(1:NSEA)        &
                                                        , MAPSF, XY )
              ENDIF
              !! Commented out unnecessary statements below for time being
              !! UA,UD are  in north-east convention and X1,X2
              !! are not actually written out below
              !DO ISEA=1, NSEA
              !  UABS   = SQRT(UA(ISEA)**2+UD(ISEA)**2)
              !  IF ( UABS .GT. 1.0 ) THEN
              !    UD(ISEA) = MOD ( 630. -                         &
              !              RADE*ATAN2(UD(ISEA),UA(ISEA)) , 360. )
              !  ELSE
              !    UD(ISEA) = UNDEF
              !    END IF
              !  UA(ISEA) = UABS
              !END DO
              !CALL W3S2XY ( NSEA, NSEA, NX+1, NY, UA(1:NSEA)        &
              !                                          , MAPSF, X1 )
              !CALL W3S2XY ( NSEA, NSEA, NX+1, NY, UD(1:NSEA)        &
              !                                          , MAPSF, X2 )
              NFIELD=2
              VARNM(1)='uwnd'
              VARNM(2)='vwnd'
              VARNL(1)='eastward_wind'
              VARNL(2)='northward_wind'
              VARNS(1)='eastward_wind'
              VARNS(2)='northward_wind'
              VARNG(1)='eastward_wind'
              VARNG(2)='northward_wind'
              VARNC(1)='wind=sqrt(U10**2+V10**2)'
              VARNC(2)='wind=sqrt(U10**2+V10**2)'
              VARND(1)=''
              VARND(2)=''
              VMIN = -990
              VMAX =  990
!
            ELSE IF ( IFI .EQ. 1 .AND. IFJ .EQ. 4 ) THEN
              FSC    = 0.1
              ENAME  = '.ast'
              UNITS  = 'K'
              IF( SMCGRD ) THEN
              ELSE ! IF(SMCGRD)
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, AS(1:NSEA)        &
                                                        , MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='ast'
              VARNL(1)='air sea temperature difference'
              VARNS(1)='air_sea_temperature_difference'
              VARNG(1)='air_sea_temperature_difference'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 4000
!
            ELSE IF ( IFI .EQ. 1 .AND. IFJ .EQ. 5 ) THEN
              FSC    = 0.01
              UNITS  = 'm'
              ENAME  = '.wlv'
              IF( SMCGRD ) THEN
              ELSE ! IF(SMCGRD)
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, WLV   , MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='wlv'
              VARNL(1)='sea surface height above sea level'
              VARNS(1)='sea_surface_height_above_sea_level'
              VARNG(1)='sea_surface_height_above_sea_level'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 10000
!
            ELSE IF ( IFI .EQ. 1 .AND. IFJ .EQ. 6 ) THEN
              FSC    = 0.001
              UNITS  = '1'
              ENAME  = '.ice'
 
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, ICE(1:NSEA), MAPSF, X1 )
              ENDIF
 
              NFIELD=1
              VARNM(1)='ice'
              VARNL(1)='sea ice area fraction'
              VARNS(1)='sea_ice_area_fraction'
              VARNG(1)='sea_ice_area_fraction'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 1000
!
            ELSE IF ( IFI .EQ. 1 .AND. IFJ .EQ. 7 ) THEN
              FSC    = 0.0001
              UNITS  = 'km-1'
              ENAME  = '.ibg'
 
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, BERG   , MAPSF, X1 )
              ENDIF
              WHERE ( X1.NE.UNDEF) X1 = X1*0.1
              NFIELD=1
              VARNM(1)='ibg'
              VARNL(1)='icebergs_damping'
              VARNS(1)='icebergs_induced_attenuation_scale_for_waves'
              VARNG(1)='icebergs_damping'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
 ELSE IF ( IFI .EQ. 1 .AND. IFJ .EQ. 8 ) THEN
              FSC    = 0.001
              UNITS  = 'Krumbein phi scale'
              ENAME  = '.d50'
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, SED_D50   , MAPSF, X1 )
              WHERE ( X1.NE.UNDEF) X1 = -LOG(X1/0.001)/LOG2
              NFIELD=1
              VARNM(1)='d50'
              VARNL(1)='grain_size'
              VARNS(1)='sediment_grain_size'
              VARNG(1)='sediment_grain_size'
              VARNC(1)=''
              VARND(1)=''
              VMIN = -10000
              VMAX = 32000
!
 ELSE IF (IFI .EQ. 1 .AND. IFJ .EQ. 9 ) THEN
              FSC = 0.001
              UNITS = 'm'
              ENAME = '.ic1'
              CALL W3S2XY (NSEA, NSEA, NX+1, NY, ICEH(1:NSEA), MAPSF, X1 )
              NFIELD=1
              VARNM(1)='ic1'
              VARNL(1)='ice thickness'
              VARNS(1)='ice_thickness'
              VARNG(1)='ice_thickness'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 30000
!
 ELSE IF (IFI .EQ. 1 .AND. IFJ .EQ. 10 ) THEN
              FSC = 0.05
              UNITS = 'm'
              ENAME = '.ic5'
              CALL W3S2XY (NSEA, NSEA, NX+1, NY, ICEF(1:NSEA), MAPSF, X1 )
              NFIELD=1
              VARNM(1)='ic5'
              VARNL(1)='maximum floe diameter'
              VARNS(1)='maximum_ice_floe_diameter'
              VARNG(1)='maximum_ice_floe_diameter'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 30000
 
            ELSE IF ( IFI .EQ. 2 .AND. IFJ .EQ. 1 ) THEN
!
              NFIELD=1
              IF (NCVARTYPEI.EQ.3) NCVARTYPE=2
              FSC    = 0.002
              UNITS  = 'm'
              ENAME  = '.hs'
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, HS    , MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='hs'
              VARNL(1)='significant height of wind and swell waves'
              VARNS(1)='sea_surface_wave_significant_height'
              VARNG(1)='significant_wave_height'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
 
            ELSE IF ( IFI .EQ. 2 .AND. IFJ .EQ. 2 ) THEN
              FSC    = 1.
              UNITS  = 'm'
              ENAME  = '.lm'
 
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, WLM, MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='lm'
              VARNL(1)='mean wave length'
              VARNS(1)='mean_wave_length'
              VARNG(1)='mean_wave_length'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 3200
!
            ELSE IF ( IFI .EQ. 2 .AND. IFJ .EQ. 3 ) THEN
              FSC    = 0.01
              UNITS  = 's'
              ENAME  = '.t02'
 
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, T02, MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='t02'
              VARNL(1)='mean period T02'
              VARNS(1)='sea_surface_wind_wave_mean_period_from_variance_' // &
                       'spectral_density_second_frequency_moment'
              VARNG(1)='mean_period_t02'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 5000
!
            ELSE IF ( IFI .EQ. 2 .AND. IFJ .EQ. 4 ) THEN
              FSC    = 0.01
              UNITS  = 's'
              ENAME  = '.t0m1'
 
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, T0M1, MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='t0m1'
              VARNL(1)='mean period T0m1'
              VARNS(1)='sea_surface_wind_wave_mean_period_from_variance' // &
                       '_spectral_density_inverse_frequency_moment'
              VARNG(1)='mean_period_t0m1'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 5000
!
            ELSE IF ( IFI .EQ. 2 .AND. IFJ .EQ. 5 ) THEN
              FSC    = 0.01
              UNITS  = 's'
              ENAME  = '.t01'
 
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, T01   , MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='t01'
              VARNL(1)='mean period T01'
              VARNS(1)='sea_surface_wind_wave_mean_period_from_variance' // &
                       '_spectral_density_first_frequency_moment'
              VARNG(1)='mean_period_t01'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 5000
!
            ELSE IF ( IFI .EQ. 2 .AND. IFJ .EQ. 6 ) THEN
              FSC    = 0.001
              UNITS  = 's-1'
              ENAME  = '.fp'
 
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, FP0   , MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='fp'
              VARNL(1)='wave peak frequency'
              VARNS(1)='sea_surface_wave_peak_frequency'
              VARNG(1)='dominant_wave_frequency'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 10000
!
            ELSE IF ( IFI .EQ. 2 .AND. IFJ .EQ. 7 ) THEN
              FSC    = 0.1
              UNITS  = 'degree'
              ENAME  = '.dir'
 
              IF( SMCGRD ) THEN
              ELSE
                 DO ISEA=1, NSEA
                   IF ( THM(ISEA) .NE. UNDEF )  THEN
                     THM(ISEA) = MOD ( 630. - RADE*THM(ISEA) , 360. )
                     END IF
                   END DO
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, THM, MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='dir'
              VARNL(1)='wave mean direction'
              VARNS(1)='sea_surface_wave_from_direction'
              VARNG(1)='wave_from_direction'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 3600
!
            ELSE IF ( IFI .EQ. 2 .AND. IFJ .EQ. 8 ) THEN
              FSC    = 0.1
              UNITS  = 'degree'
              ENAME  = '.spr'
 
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, THS   , MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='spr'
              VARNL(1)='directional spread'
              VARNS(1)='sea_surface_wave_directional_spread'
              VARNG(1)='directional_spread'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 900
!
            ELSE IF ( IFI .EQ. 2 .AND. IFJ .EQ. 9 ) THEN
              FSC    = 1.
              UNITS  = 'degree'
              ENAME  = '.dp'
 
              IF( SMCGRD ) THEN
              ELSE
                 DO ISEA=1, NSEA
                   IF ( THP0(ISEA) .NE. UNDEF ) THEN
                     THP0(ISEA) = MOD ( 630-RADE*THP0(ISEA) , 360. )
                     END IF
                   END DO
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, THP0  , MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='dp'
              VARNL(1)='peak direction'
              VARNS(1)='sea_surface_wave_peak_direction'
              VARNG(1)='dominant_wave_direction'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 360
!
            ELSE IF ( IFI .EQ. 2 .AND. IFJ .EQ. 10 ) THEN
              FSC    = 0.0002
              UNITS  = 'm'
              ENAME  = '.hig'
 
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, HSIG, MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='hig'
              VARNL(1)='infragravity_wave_height'
              VARNS(1)='sea_surface_wave_infragravity_significant_height'
              VARNG(1)='infragravity_significant_wave_height'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 5000
!
            ELSE IF ( IFI .EQ. 2 .AND. IFJ .EQ. 11 ) THEN
              FSC    = 0.002
              UNITS  = 'm'
              ENAME  = '.mxe'
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, STMAXE, MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='stmaxe'
              VARNL(1)='expected maximum sea surface elevation (nonlinear,2nd order)'
              VARNS(1)='expected maximum sea surface elevation (nonlinear,2nd order)'
              VARNG(1)='expected maximum sea surface elevation (nonlinear,2nd order)'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 2 .AND. IFJ .EQ. 12 ) THEN
              FSC    = 0.002
              UNITS  = 'm'
              ENAME  = '.mxes'
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, STMAXD, MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='stmaxd'
              VARNL(1)='standard deviation of maximum sea surface elevation (nonlinear,2nd order)'
              VARNS(1)='std of expected maximum sea surface elevation (nonlinear,2nd order)'
              VARNG(1)='standard deviation of maximum sea surface elevation (nonlinear,2nd order)'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 2 .AND. IFJ .EQ. 13 ) THEN
              FSC    = 0.002
              UNITS  = 'm'
              ENAME  = '.mxh'
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, HMAXE, MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='hmaxe'
              VARNL(1)='expected maximum wave height (linear, 1st order)'
              VARNS(1)='expected maximum wave height (linear, 1st order)'
              VARNG(1)='expected maximum wave height (linear, 1st order)'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 2 .AND. IFJ .EQ. 14 ) THEN
              FSC    = 0.002
              UNITS  = 'm'
              ENAME  = '.mxhc'
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, HCMAXE, MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='hcmaxe'
              VARNL(1)='expected maximum wave height from crest (linear, 1st order)'
              VARNS(1)='expected maximum wave height from crest (linear, 1st order)'
              VARNG(1)='expected maximum wave height from crest (linear, 1st order)'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 2 .AND. IFJ .EQ. 15 ) THEN
              FSC    = 0.002
              UNITS  = 'm'
              ENAME  = '.sdmh'
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, HMAXD, MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='hmaxd'
              VARNL(1)='STD of maximum wave height (linear, 1st order)'
              VARNS(1)='STD of maximum wave height (linear, 1st order)'
              VARNG(1)='STD of maximum wave height (linear, 1st order)'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 2 .AND. IFJ .EQ. 16 ) THEN
              FSC    = 0.002
              UNITS  = 'm'
              ENAME  = '.sdmhc'
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, HCMAXD, MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='hcmaxd'
              VARNL(1)='STD of maximum wave height from crest (linear, 1st order)'
              VARNS(1)='STD of maximum wave height from crest (linear, 1st order)'
              VARNG(1)='STD of maximum wave height from crest (linear, 1st order)'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 2 .AND. IFJ .EQ. 17 ) THEN
              FSC    = 0.001
              UNITS  = '1'
              ENAME  = '.wbt'
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, WBT, MAPSF, X1 )
              END IF
              NFIELD=1
              VARNM(1)='wbt'
              VARNL(1)='domiant wave breaking probability'
              VARNS(1)='domiant_wave_breaking_probability'
              VARNG(1)='domiant_wave_breaking_probability'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 1000
!
            ELSE IF ( IFI .EQ. 3 .AND. IFJ .EQ. 1 ) THEN
              ! Information for spectral
              FLFRQ  = .TRUE.
              NFIELD = 1
              VARNM(1)='ef'
              VARNL(1)='wave_elevation_spectrum'
              VARNG(1)=VARNS(1)
              VARNC(1)=''
              VARND(1)=''
              IF (NCVARTYPE.LE.3) THEN
                UNITS  = 'log10(m2 s+1E-12)'
                VARNS(1)='base_ten_logarithm_of_power_spectral_density_of_surface_elevation'
                FSC     = 0.0004
              ELSE
                UNITS  = 'm2 s'
                VARNS(1)='power_spectral_density_of_surface_elevation'
                FSC     = 1.
                END IF
              ENAME  = '.ef'
              VARNG(1)=VARNS(1)
              VARNG(1)=VARNS(1)
              I1F=E3DF(2,1)
              I2F=E3DF(3,1)
              DO IK=I1F,I2F
                IF( SMCGRD ) THEN
                ELSE
                   CALL W3S2XY ( NSEA, NSEA, NX+1,NY,EF(:,IK),MAPSF, XX )
                ENDIF
                IF (NCVARTYPE.EQ.2) WHERE ( XX.GE.0.) XX = ALOG10(XX+1E-12)
                XK(:,:,IK)=XX
                END DO
              VMIN = -30000
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 3 .AND. IFJ .EQ. 2 ) THEN
              ! Information for spectral
              FLFRQ  = .TRUE.
              NFIELD = 1
              FSC     = 0.1
              VARNM(1)='th1m'
              VARNL(1)='mean wave direction frequency spectrum'
              VARNS(1)='sea_surface_wave_from_direction_frequency_spectrum'
              VARNG(1)=VARNS(1)
              VARNC(1)=''
              VARND(1)=''
              UNITS  = 'degree'
              ENAME  = '.th1m'
              I1F=E3DF(2,2)
              I2F=E3DF(3,2)
              DO IK=I1F,I2F
                IF( SMCGRD ) THEN
                ELSE
                   CALL W3S2XY ( NSEA, NSEA, NX+1,NY,TH1M(:,IK),MAPSF, XX )
                ENDIF
                XK(:,:,IK)=XX
                END DO
              VMIN = 0
              VMAX = 3600
!
            ELSE IF ( IFI .EQ. 3 .AND. IFJ .EQ. 3 ) THEN
              ! Information for spectral
              FLFRQ  = .TRUE.
              NFIELD = 1
              FSC     = 0.01
              VARNM(1)='sth1m'
              VARNL(1)='spreading frequency spectrum'
              VARNS(1)='sea_surface_wave_spreading_spectrum'
              VARNG(1)=VARNS(1)
              VARNC(1)=''
              VARND(1)=''
              UNITS  = 'degree'
              ENAME  = '.sth1m'
              I1F=E3DF(2,3)
              I2F=E3DF(3,3)
              DO IK=I1F,I2F
                IF( SMCGRD ) THEN
                ELSE
                   CALL W3S2XY ( NSEA, NSEA, NX+1,NY,STH1M(:,IK),MAPSF, XX )
                ENDIF
                XK(:,:,IK)=XX
                END DO
              VMIN = 0
              VMAX = 9000
!
            ELSE IF ( IFI .EQ. 3 .AND. IFJ .EQ. 4 ) THEN
              ! Information for spectral
              FLFRQ  = .TRUE.
              NFIELD = 1
              FSC     = 0.1
              VARNM(1)='th2m'
              VARNL(1)='second mean wave direction frequency spectrum'
              VARNS(1)='sea_surface_wave_from_direction_frequency_spectrum_from_second_moments'
              VARNG(1)=VARNS(1)
              VARNC(1)=''
              VARND(1)=''
              UNITS  = 'degree'
              ENAME  = '.th2m'
              I1F=E3DF(2,4)
              I2F=E3DF(3,4)
              DO IK=I1F,I2F
                IF( SMCGRD ) THEN
                ELSE
                   CALL W3S2XY ( NSEA, NSEA, NX+1,NY,TH2M(:,IK),MAPSF, XX )
                ENDIF
                XK(:,:,IK)=XX
                END DO
              VMIN = 0
              VMAX = 3600
!
            ELSE IF ( IFI .EQ. 3 .AND. IFJ .EQ. 5 ) THEN
              ! Information for spectral
              FLFRQ  = .TRUE.
              NFIELD = 1
              FSC     = 0.01
              VARNM(1)='sth2m'
              VARNL(1)='second spreading frequency spectrum'
              VARNS(1)='sea_surface_wave_spreading_spectrum_from_second_moments'
              VARNG(1)=VARNS(1)
              VARNC(1)=''
              VARND(1)=''
              UNITS  = 'degree'
              ENAME  = '.sth2m'
              I1F=E3DF(2,5)
              I2F=E3DF(3,5)
              DO IK=I1F,I2F
                IF( SMCGRD ) THEN
                ELSE
                   CALL W3S2XY ( NSEA, NSEA, NX+1,NY,STH2M(:,IK),MAPSF, XX )
                ENDIF
                XK(:,:,IK)=XX
                END DO
              VMIN = 0
              VMAX = 9000
!
            ELSE IF ( IFI .EQ. 3 .AND. IFJ .EQ. 6 ) THEN
              ! Information for spectral
              FLFRQ  = .TRUE.
              FSC    = 0.001
              UNITS  = 'm-1'
              ENAME  = '.wn'
              I1F=1
              I2F=NK
              DO IK=1,NK
                IF( SMCGRD ) THEN
                ELSE
                   CALL W3S2XY ( NSEA, NSEA, NX+1, NY, WN(IK,:), MAPSF, XX )
                ENDIF
                XK(:,:,IK)=XX
                END DO
              NFIELD=1
              VARNM(1)='wn'
              VARNL(1)='wave numbers'
              VARNS(1)='wave_numbers'
              VARNG(1)='wave_numbers'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 4 .AND. IFJ .EQ. 1 ) THEN
              FSC    = 0.002
              UNITS  = 'm'
              WRITE(ENAME,'(A4,I1)') '.phs',IPART
 
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PHS(:,IPART), MAPSF, X1 )
              ENDIF
 
              NFIELD=1
              WRITE(VARNM(1),'(A3,I1)') 'phs',IPART
              WRITE(VARNL(1),'(A,I1)') 'wave significant height partition ',IPART
              WRITE(VARNS(1),'(A,I1)') 'sea_surface_wave_significant_height_partition_',IPART
              WRITE(VARNG(1),'(A,I1)') 'significant_wave_height_partition_',IPART
              VARNC(1)=PARTCOM
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 4 .AND. IFJ .EQ. 2 ) THEN
              FSC    = 0.01
              UNITS  = 's'
              WRITE(ENAME,'(A4,I1)') '.ptp',IPART
 
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PTP(:,IPART), MAPSF, X1 )
              ENDIF
 
              NFIELD=1
              WRITE(VARNM(1),'(A3,I1)') 'ptp',IPART
              WRITE(VARNL(1),'(A,I1)') 'peak period partition ',IPART
              WRITE(VARNS(1),'(A,I1)') 'sea_surface_wave_period_at_' // &
                                       'variance_spectral_density_maximum_partition_',IPART
              WRITE(VARNG(1),'(A,I1)') 'dominant_wave_period_partition_',IPART
              VARNC(1)=PARTCOM
              VARND(1)=''
              VMIN = 0
              VMAX = 10000
!
            ELSE IF ( IFI .EQ. 4 .AND. IFJ .EQ. 3 ) THEN
              FSC    = 1.
              UNITS  = 'm'
              WRITE(ENAME,'(A4,I1)') '.plp',IPART
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PLP(:,IPART), MAPSF, X1 )
              ENDIF
 
              NFIELD=1
              WRITE(VARNM(1),'(A3,I1)') 'plp',IPART
              WRITE(VARNL(1),'(A,I1)') 'peak wave length partition ',IPART
              WRITE(VARNS(1),'(A,I1)') 'peak_wave_length_partition_',IPART
              WRITE(VARNG(1),'(A,I1)') 'peak_wave_length_partition_',IPART
              VARNC(1)=PARTCOM
              VARND(1)=''
              VMIN = 0
              VMAX = 10000
!
            ELSE IF ( IFI .EQ. 4 .AND. IFJ .EQ. 4 ) THEN
              FSC    = 0.1
              UNITS  = 'degree'
              WRITE(ENAME,'(A5,I1)') '.pdir',IPART
              IF( SMCGRD ) THEN
              ELSE
                DO ISEA=1, NSEA
                  IF ( PDIR(ISEA,IPART) .NE. UNDEF ) THEN
                     PDIR(ISEA,IPART) = MOD ( 630-RADE*PDIR(ISEA,IPART) , 360. )
                  END IF
                END DO
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PDIR(:,IPART), MAPSF, X1 )
              ENDIF
              NFIELD=1
              WRITE(VARNM(1),'(A4,I1)') 'pdir',IPART
              WRITE(VARNL(1),'(A,I1)') 'wave mean direction partition ',IPART
              WRITE(VARNS(1),'(A,I1)') 'sea_surface_wave_from_direction_partition_',IPART
              WRITE(VARNG(1),'(A,I1)') 'wave_from_direction_partition_',IPART
              VARNC(1)=PARTCOM
              VARND(1)=''
              VMIN = 0
              VMAX = 3600
!
            ELSE IF ( IFI .EQ. 4 .AND. IFJ .EQ. 5 ) THEN
              FSC    = 0.1
              UNITS  = 'degree'
              WRITE(ENAME,'(A5,I1)') '.pspr',IPART
 
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PSI(:,IPART), MAPSF, X1 )
              ENDIF
 
              NFIELD=1
              WRITE(VARNM(1),'(A4,I1)') 'pspr',IPART
              WRITE(VARNL(1),'(A,I1)') 'directional spread partition ',IPART
              WRITE(VARNS(1),'(A,I1)') 'sea_surface_wave_directional_spread_partition_',IPART
              WRITE(VARNG(1),'(A,I1)') 'directional_spread_partition_',IPART
              VARNC(1)=PARTCOM
              VARND(1)=''
              VMIN = 0
              VMAX = 900
!
            ELSE IF ( IFI .EQ. 4 .AND. IFJ .EQ. 6 ) THEN
              FSC    = 0.001
              UNITS  = '1'
              WRITE(ENAME,'(A4,I1)') '.pws',IPART
 
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PWS(:,IPART), MAPSF, X1 )
              ENDIF
 
              NFIELD=1
              WRITE(VARNM(1),'(A3,I1)') 'pws',IPART
              WRITE(VARNL(1),'(A,I1)') 'wind sea fraction in partition ',IPART
              WRITE(VARNS(1),'(A,I1)') 'wind_sea_fraction_in_partition_',IPART
              WRITE(VARNG(1),'(A,I1)') 'wind_sea_fraction_in_partition_',IPART
              VARNC(1)=PARTCOM
              VARND(1)=''
              VMIN = 0
              VMAX = 1000
!
            ELSE IF ( IFI .EQ. 4 .AND. IFJ .EQ. 7 ) THEN
              FSC    = 0.1
              UNITS  = 'degree'
              WRITE(ENAME,'(A4,I1)') '.pdp',IPART
              IF( SMCGRD ) THEN
              ELSE
                DO ISEA=1, NSEA
                  IF ( PTHP0(ISEA,IPART) .NE. UNDEF ) THEN
                    PTHP0(ISEA,IPART) = MOD ( 630-RADE*PTHP0(ISEA,IPART) , 360. )
                    END IF
                  END DO
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PTHP0(:,IPART), MAPSF, X1 )
              END IF
              NFIELD=1
              WRITE(VARNM(1),'(A3,I1)') 'pdp',IPART
              WRITE(VARNL(1),'(A,I1)') &
                 'peak direction partition ',IPART
              WRITE(VARNS(1),'(A,I1)') &
                 'sea_surface_wave_peak_from_direction_partition_',IPART
              WRITE(VARNG(1),'(A,I1)') &
                  'dominant_wave_from_direction_partition_',IPART
              VARNC(1)=PARTCOM
              VARND(1)=''
              VMIN = 0
              VMAX = 3600
!
            ELSE IF ( IFI .EQ. 4 .AND. IFJ .EQ. 8 ) THEN
              FSC    = 0.01
              UNITS  = '1'
              WRITE(ENAME,'(A4,I1)') '.pqp',IPART
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PQP(:,IPART), MAPSF, X1 )
              ENDIF
              NFIELD=1
              WRITE(VARNM(1),'(A3,I1)') 'pqp',IPART
              WRITE(VARNL(1),'(A,I1)') 'peakedness partition ',IPART
              WRITE(VARNS(1),'(A,I1)') &
                  'sea_surface_wave_peakedness_partition_',IPART
              WRITE(VARNG(1),'(A,I1)') &
                  'wave_peakedness_partition_',IPART
              VARNC(1)=PARTCOM
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 4 .AND. IFJ .EQ. 9 ) THEN
              FSC    = 0.01
              UNITS  = '1'
              WRITE(ENAME,'(A4,I1)') '.ppe',IPART
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PPE(:,IPART), MAPSF, X1 )
              ENDIF
              NFIELD=1
              WRITE(VARNM(1),'(A3,I1)') 'ppe',IPART
              WRITE(VARNL(1),'(A,I1)') &
                  'peak enhancement factor partition ',IPART
              WRITE(VARNS(1),'(A,I1)') &
                  'wave_peak_enhancement_factor_partition_',IPART
              WRITE(VARNG(1),'(A,I1)') &
                  'wave_peak_enhancement_factor_partition_',IPART
              VARNC(1)='JONSWAP peak enhancement factor; ' // PARTCOM
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 4 .AND. IFJ .EQ. 10 ) THEN
              FSC    = 0.0001
              UNITS  = 's-1'
              WRITE(ENAME,'(A4,I1)') '.pgw',IPART
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PGW(:,IPART), MAPSF, X1 )
              ENDIF
              NFIELD=1
              WRITE(VARNM(1),'(A3,I1)') 'pgw',IPART
              WRITE(VARNL(1),'(A,I1)') &
                  'frequency width partition ',IPART
              WRITE(VARNS(1),'(A,I1)') &
                  'Gaussian_frequency_spread_partition_',IPART
              WRITE(VARNG(1),'(A,I1)') &
                  'Gaussian_frequency_spread_partition_',IPART
              VARNC(1)='Gaussian least-square fit to ' // &
                       'omni-directional spectral partition; ' // PARTCOM
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 4 .AND. IFJ .EQ. 11 ) THEN
              FSC    = 0.0001
              UNITS  = '1'
              WRITE(ENAME,'(A4,I1)') '.psw',IPART
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PSW(:,IPART), MAPSF, X1 )
              ENDIF
              NFIELD=1
              WRITE(VARNM(1),'(A3,I1)') 'psw',IPART
              WRITE(VARNL(1),'(A,I1)') 'spectral width partition ',IPART
              WRITE(VARNS(1),'(A,I1)') &
                  'sea_surface_wave_spectral_width_partition_',IPART
              WRITE(VARNG(1),'(A,I1)') &
                  'wave_spectral_width_partition_',IPART
              VARNC(1)=PARTCOM
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 4 .AND. IFJ .EQ. 12 ) THEN
              FSC    = 0.01
              UNITS  = 's'
              WRITE(ENAME,'(A7,I1)') '.ptm10c',IPART
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PTM1(:,IPART), MAPSF, X1 )
              ENDIF
              NFIELD=1
              WRITE(VARNM(1),'(A6,I1)') 'ptm10c',IPART
              WRITE(VARNL(1),'(A,I1)') 'mean period Tm10 partition ',IPART
              WRITE(VARNS(1),'(A,I1)') &
                  'sea_surface_wave_mean_period_tm10_partition_',IPART
              WRITE(VARNG(1),'(A,I1)') &
                  'mean_wave_period_Tm10_partition_',IPART
              VARNC(1)=PARTCOM
              VARND(1)=''
              VMIN = 0
              VMAX = 10000
!
            ELSE IF ( IFI .EQ. 4 .AND. IFJ .EQ. 13 ) THEN
              FSC    = 0.01
              UNITS  = 's'
              WRITE(ENAME,'(A6,I1)') '.pt01c',IPART
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PT1(:,IPART), MAPSF, X1 )
              ENDIF
              NFIELD=1
              WRITE(VARNM(1),'(A5,I1)') 'pt01c',IPART
              WRITE(VARNL(1),'(A,I1)') 'mean period T01 partition ',IPART
              WRITE(VARNS(1),'(A,I1)') &
                  'sea_surface_wave_mean_period_t01_partition_',IPART
              WRITE(VARNG(1),'(A,I1)') &
                  'mean_wave_period_T01_partition_',IPART
              VARNC(1)=PARTCOM
              VARND(1)=''
              VMIN = 0
              VMAX = 10000
!
            ELSE IF ( IFI .EQ. 4 .AND. IFJ .EQ. 14 ) THEN
              FSC    = 0.01
              UNITS  = 's'
              WRITE(ENAME,'(A6,I1)') '.pt02c',IPART
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PT2(:,IPART), MAPSF, X1 )
              ENDIF
              NFIELD=1
              WRITE(VARNM(1),'(A5,I1)') 'pt02c',IPART
              WRITE(VARNL(1),'(A,I1)') 'mean period T02 partition ',IPART
              WRITE(VARNS(1),'(A,I1)') &
                  'sea_surface_wave_mean_period_t02_partition_',IPART
              WRITE(VARNG(1),'(A,I1)') &
                  'mean_wave_period_T02_partition_',IPART
              VARNC(1)=PARTCOM
              VARND(1)=''
              VMIN = 0
              VMAX = 10000
!
            ELSE IF ( IFI .EQ. 4 .AND. IFJ .EQ. 15 ) THEN
              FSC    = 0.02
              UNITS  = 'm2 s rad-1'
              WRITE(ENAME,'(A4,I1)') '.pep',IPART
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PEP(:,IPART), MAPSF, X1 )
              ENDIF
              NFIELD=1
              WRITE(VARNM(1),'(A3,I1)') 'pep',IPART
              WRITE(VARNL(1),'(A,I1)') 'energy at peak frequency partition ',IPART
              WRITE(VARNS(1),'(A,I1)') &
                  'sea_surface_wave_energy_at_variance_spectral_density_maximum_partition_',IPART
              WRITE(VARNG(1),'(A,I1)') &
                  'wave_energy_at_variance_spectral_density_maximum_partition_',IPART
              VARNC(1)=PARTCOM
              VARND(1)=''
              VMIN = 0
              VMAX = 10000
!
            ELSE IF ( IFI .EQ. 4 .AND. IFJ .EQ. 16 ) THEN
              FSC    = 0.001
              UNITS  = '1'
              ENAME  = '.tws'
 
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PWST(:), MAPSF, X1 )
              ENDIF
 
              NFIELD=1
              VARNM(1)='tws'
              VARNL(1)='wind sea fraction'
              VARNS(1)='wind_sea_fraction'
              VARNG(1)='wind_sea_fraction'
              VARNC(1)=PARTCOM
              VARND(1)=''
              VMIN = 0
              VMAX = 1000
!
            ELSE IF ( IFI .EQ. 4 .AND. IFJ .EQ. 17 ) THEN
              FSC    = 1.
              UNITS  = '1'
              ENAME  = '.pnr'
 
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PNR(:), MAPSF, X1 )
              ENDIF
 
              NFIELD=1
              VARNM(1)='pnr'
              VARNL(1)='number of wave partitions'
              VARNS(1)='number_of_wave_partitions'
              VARNG(1)='number_of_wave_partitions'
              VARNC(1)=PARTCOM
              VARND(1)=''
              VMIN = 0
              VMAX = 100
!
            ELSE IF ( IFI .EQ. 5 .AND. IFJ .EQ. 1 ) THEN
              FSC    = 0.01
              ENAME  = '.ust'
              UNITS  = 'm s-1'
              !! Note - UST and USTDIR read in from .ww3 file are X-Y vectors
              DO ISEA=1, NSEA
                UABS   = SQRT(UST(ISEA)**2+USTDIR(ISEA)**2)
                 IF (UABS.GE.10.) THEN
                  UST(ISEA)=UNDEF
                  USTDIR(ISEA)=UNDEF
                  END IF
                 END DO
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY (NSEA,NSEA,NX+1,NY, UST   (1:NSEA), MAPSF, XX )
                CALL W3S2XY (NSEA,NSEA,NX+1,NY, USTDIR(1:NSEA), MAPSF, XY )
              ENDIF ! SMCGRD
              !! Commented out unnecessary statements below for time being
              !! UST,USTDIR are in north-east convention and X1,X2
              !! are not actually written out below
              !DO ISEA=1, NSEA
              !  UABS   = SQRT(UST(ISEA)**2+USTDIR(ISEA)**2)
              !  IF ( UST(ISEA) .EQ. UNDEF ) THEN
              !      USTDIR(ISEA) = UNDEF
              !      UABS         = UNDEF
              !    ELSE IF ( UABS .GT. 0.05 ) THEN
              !      USTDIR(ISEA) = MOD ( 630. -                     &
              !        RADE*ATAN2(USTDIR(ISEA),UST(ISEA)) , 360. )
              !    ELSE
              !      USTDIR(ISEA) = UNDEF
              !    END IF
              !  UST(ISEA) = UABS
              !  END DO
              !CALL W3S2XY (NSEA,NSEA,NX+1,NY, UST   (1:NSEA) , MAPSF, X1 )
              !CALL W3S2XY (NSEA,NSEA,NX+1,NY, USTDIR(1:NSEA) , MAPSF, X2 )
              NFIELD=2
              VARNM(1)='uust'
              VARNM(2)='vust'
              VARNL(1)='eastward friction velocity'
              VARNL(2)='northward friction velocity'
              VARNS(1)='eastward_friction_velocity'
              VARNS(2)='northward_friction_velocity'
              VARNG(1)='eastward_friction_velocity'
              VARNG(2)='northward_friction_velocity'
              VARNC(1)='ust=sqrt(uust**2+vust**2)'
              VARNC(2)='ust=sqrt(uust**2+vust**2)'
              VARND(1)=''
              VARND(2)=''
              VMIN = -9900
              VMAX =  9900
!
            ELSE IF ( IFI .EQ. 5 .AND. IFJ .EQ. 2 ) THEN
              FSC    = 1.E-5
              UNITS  = '1'
              ENAME  = '.cha'
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, CHARN(1:NSEA), MAPSF, X1 )
              ENDIF ! SMCGRD
              NFIELD=1
              VARNM(1)='cha'
              VARNL(1)='charnock coefficient for surface roughness length for momentum in air'
              VARNS(1)='charnock_coefficient_for_surface_roughness_length_for_momentum_in_air'
              VARNG(1)='charnock_coefficient'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 32700
!
            ELSE IF ( IFI .EQ. 5 .AND. IFJ .EQ. 3 ) THEN
              FSC    = 0.1           !0.01
              UNITS  = 'kW m-1'
              ENAME  = '.cge'
              CGE=CGE*0.001  ! from W / m to kW / m
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, CGE(1:NSEA), MAPSF, X1 )
              ENDIF ! SMCGRD
              NFIELD=1
              VARNM(1)='cge'
              VARNL(1)='wave energy flux'
              VARNS(1)='sea_surface_wind_wave_energy_flux'
              VARNG(1)='wave_energy_flux'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 9990
!
            ELSE IF ( IFI .EQ. 5 .AND. IFJ .EQ. 4 ) THEN
              IF (NCVARTYPEI.EQ.3) NCVARTYPE=4
              FSC    = 0.1
              UNITS  = 'W m-2'
              ENAME  = '.faw'
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PHIAW(1:NSEA) , MAPSF, X1 )
              ENDIF ! SMCGRD
              NFIELD=1
              VARNM(1)='faw'
              VARNL(1)='wind to wave energy flux'
              VARNS(1)='wind_mixing_energy_flux_into_sea_water'
              VARNG(1)='wind_to_wave_energy_flux'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 9990
!
            ELSE IF ( IFI .EQ. 5 .AND. IFJ .EQ. 5 ) THEN
              FSC    = 0.000001
              UNITS  = 'm2 s-2'
              ENAME  = '.taw'
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, TAUWIX(1:NSEA)      &
                                                        , MAPSF, XX )
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, TAUWIY(1:NSEA)      &
                                                        , MAPSF, XY )
              ENDIF ! SMCGRD
              !! Commented out unnecessary statements below for time being
              !! TAUWIX, TAUWIY are in north-east convention and X1,X2
              !! are not actually written out below
              !DO ISEA=1, NSEA
              !  CABS   = SQRT(TAUWIX(ISEA)**2+TAUWIY(ISEA)**2)
              !  IF ( CABS .NE. UNDEF ) THEN
              !      TAUWIY(ISEA) = MOD ( 630. -                         &
              !            RADE*ATAN2(TAUWIY(ISEA),TAUWIX(ISEA)) , 360. )
              !    ELSE
              !      TAUWIY(ISEA) = UNDEF
              !    END IF
              !  TAUWIX(ISEA) = CABS
              !  END DO
              !CALL W3S2XY ( NSEA, NSEA, NX+1, NY, TAUWIX, MAPSF, X1 )
              !CALL W3S2XY ( NSEA, NSEA, NX+1, NY, TAUWIY, MAPSF, X2 )
              NFIELD=2
              VARNM(1)='utaw'
              VARNM(2)='vtaw'
              VARNL(1)='eastward wave supported wind stress'
              VARNL(2)='northward wave supported wind stress'
              VARNS(1)='eastward_wave_supported_wind_stress'
              VARNS(2)='northward_wave_supported_wind_stress'
              VARNG(1)='eastward_wave_supported_wind_stress'
              VARNG(2)='northward_wave_supported_wind_stress'
              VARNC(1)='taw=sqrt(utaw**2+vtaw**2)'
              VARNC(2)='taw=sqrt(utaw**2+vtaw**2)'
              VARND(1)=''
              VARND(2)=''
              VMIN = -32000
              VMAX =  32000
!
            ELSE IF ( IFI .EQ. 5 .AND. IFJ .EQ. 6 ) THEN
              FSC    = 0.0001
              ENAME  = '.twa'
              UNITS  = 'm2 s-2'
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, TAUWNX(1:NSEA)   &
                                                        , MAPSF, XX )
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, TAUWNY(1:NSEA)   &
                                                        , MAPSF, XY )
              ENDIF ! SMCGRD
              NFIELD=2
              VARNM(1)='utwa'
              VARNM(2)='vtwa'
              VARNL(1)='eastward wave to wind stress'
              VARNL(2)='northward wave to wind stress'
              VARNS(1)='eastward_wave_to_wind_stress'
              VARNS(2)='northward_wave_to_wind_stress'
              VARNG(1)='eastward_wave_to_wind_stress'
              VARNG(2)='northward_wave_to_wind_stress'
              VARNC(1)='twa=sqrt(utwa**2+vtwa**2)'
              VARNC(2)='twa=sqrt(utwa**2+vtwa**2)'
              VARND(1)=''
              VARND(2)=''
              VMIN = -32000
              VMAX =  32000
!
            ELSE IF ( IFI .EQ. 5 .AND. IFJ .EQ. 7 ) THEN
              FSC    = 0.0001
              UNITS  = '1'
              ENAME  = '.wcc'
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, WHITECAP(1:NSEA,1) &
                                                        , MAPSF, X1 )
              ENDIF ! SMCGRD
              NFIELD=1
              VARNM(1)='wcc'
              VARNL(1)='whitecap coverage'
              VARNS(1)='whitecap_coverage'
              VARNG(1)='whitecap_coverage'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 10000
!
            ELSE IF ( IFI .EQ. 5 .AND. IFJ .EQ. 8 ) THEN
              FSC    = 0.001
              UNITS  = 'm'
              ENAME  = '.wcf'
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, WHITECAP(1:NSEA,2) &
                                                        , MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='wcf'
              VARNL(1)='whitecap foam thickness'
              VARNS(1)='whitecap_foam_thickness'
              VARNG(1)='whitecap_foam_thickness'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 10000
!
            ELSE IF ( IFI .EQ. 5 .AND. IFJ .EQ. 9 ) THEN
              FSC    = 0.002
              UNITS  = 'm'
              ENAME  = '.wch'
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, WHITECAP(1:NSEA,3) &
                                                        , MAPSF, X1 )
              ENDIF ! SMCGRD
              NFIELD=1
              VARNM(1)='wch'
              VARNL(1)='significant breaking wave height'
              VARNS(1)='significant_breaking_wave_height'
              VARNG(1)='significant_breaking_wave_height'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 5 .AND. IFJ .EQ. 10 ) THEN
              FSC    = 0.0001
              UNITS  = '1'
              ENAME  = '.wcm'
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, WHITECAP(1:NSEA,4) &
                                                        , MAPSF, X1 )
              ENDIF ! SMCGRD
              NFIELD=1
              VARNM(1)='wcm'
              VARNL(1)='whitecap moment'
              VARNS(1)='whitecap_moment'
              VARNG(1)='whitecap_moment'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 10000
!
            ELSE IF ( IFI .EQ. 5 .AND. IFJ .EQ. 11 ) THEN
              FSC    = 0.002
              UNITS  = 's'
              ENAME  = '.fws'
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, TWS(1:NSEA), MAPSF, X1 )
              ENDIF ! SMCGRD
              NFIELD=1
              VARNM(1)='fws'
              VARNL(1)='Wind_sea_mean_period_T0M1'
              VARNS(1)='Wind_sea_mean_period_T0M1'
              VARNG(1)='Wind_sea_mean_period_T0M1'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 6 .AND. IFJ .EQ. 1 ) THEN
              FSC    = 10.
              FSC3    = 1.
              UNITS  = 'N m-1'
              ENAME  = '.sxy'
 
              IF( SMCGRD ) THEN
              ELSE
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, SXX(1:NSEA)       &
                                                        , MAPSF, X1 )
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, SYY(1:NSEA)       &
                                                        , MAPSF, X2 )
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, SXY(1:NSEA)       &
                                                        , MAPSF, XY )
              ENDIF ! SMCGRD
              NFIELD=3
              VARNM(1)='sxx'
              VARNM(2)='syy'
              VARNM(3)='sxy'
              VARNL(1)='Radiation stress component Sxx'
              VARNL(2)='Radiation stress component Syy'
              VARNL(3)='Radiation stress component Sxy'
              VARNS(1)='Radiation_stress_component_Sxx'
              VARNS(2)='Radiation_stress_component_Syy'
              VARNS(3)='Radiation_stress_component_Sxy'
              VARNS(1)='radiation_stress_component_sxx'
              VARNS(2)='radiation_stress_component_syy'
              VARNS(3)='radiation_stress_component_sxy'
              VARNC(1)=''
              VARNC(2)=''
              VARNC(3)=''
              VARND(1)=''
              VARND(2)=''
              VARND(3)=''
              VMIN = -3000
              VMAX = 3000
!
            ELSE IF ( IFI .EQ. 6 .AND. IFJ .EQ. 2 ) THEN
              FSC    = 0.000001
              UNITS  = 'm2 s-2'
              ENAME  = '.two'
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, TAUOX(1:NSEA)      &
                                                        , MAPSF, XX )
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, TAUOY(1:NSEA)      &
                                                        , MAPSF, XY )
              ENDIF ! SMCGRD
              NFIELD=2
              VARNM(1)='utwo'
              VARNM(2)='vtwo'
              VARNL(1)='eastward wave to ocean stress'
              VARNL(2)='northward wave to ocean stress'
              VARNS(1)='eastward_wave_to_ocean_stress'
              VARNS(2)='northward_wave_to_ocean_stress'
              VARNG(1)='eastward_wave_to_ocean_stress'
              VARNG(2)='northward_wave_to_ocean_stress'
              VARNC(1)='two=sqrt(utwo**2+vtwo**2)'
              VARNC(2)='two=sqrt(utwo**2+vtwo**2)'
              VARND(1)=''
              VARND(2)=''
              VMIN = -32000
              VMAX =  32000
!
            ELSE IF ( IFI .EQ. 6 .AND. IFJ .EQ. 3 ) THEN
              FSC    = 0.1
              UNITS  = 'm2 s-2'
              ENAME  = '.bhd'
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, BHD(1:NSEA)  &
                                                        , MAPSF, X1 )
              ENDIF ! SMCGRD
              NFIELD=1
              VARNM(1)='bhd'
              VARNL(1)='radiation pressure (Bernouilli Head)'
              VARNS(1)='radiation_pressure'
              VARNG(1)='radiation_pressure'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 1000
!
            ELSE IF ( IFI .EQ. 6 .AND. IFJ .EQ. 4 ) THEN
              IF (NCVARTYPEI.EQ.3) NCVARTYPE=4
              FSC    = 0.1
              UNITS  = 'W m-2'
              ENAME  = '.foc'
              DO ISEA=1, NSEA
                 PHIOC(ISEA)=MIN(3000.,PHIOC(ISEA))
                 END DO
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PHIOC(1:NSEA)  &
                                                        , MAPSF, X1 )
              ENDIF ! SMCGRD
              NFIELD=1
              VARNM(1)='foc'
              VARNL(1)='wave to ocean energy flux'
              VARNS(1)='wave_to_ocean_energy_flux'
              VARNG(1)='wave_to_ocean_energy_flux'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 9990
!
            ELSE IF ( IFI .EQ. 6 .AND. IFJ .EQ. 5 ) THEN
              FSC    = 0.001
              UNITS  = 'm2 s-1'
              ENAME  = '.tus'
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, TUSX(1:NSEA)      &
                                                        , MAPSF, XX )
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, TUSY(1:NSEA)      &
                                                        , MAPSF, XY )
              ENDIF ! SMCGRD
              DO ISEA=1, NSEA
                CABS   = SQRT(TUSX(ISEA)**2+TUSY(ISEA)**2)
                IF ( CABS .NE. UNDEF ) THEN
                    TUSY(ISEA) = MOD ( 630. -                         &
                          RADE*ATAN2(TUSY(ISEA),TUSX(ISEA)) , 360. )
                  ELSE
                    TUSY(ISEA) = UNDEF
                  END IF
                TUSX(ISEA) = CABS
                END DO
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY,TUSX,MAPSF, X1 )
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY,TUSY,MAPSF, X2 )
              ENDIF ! SMCGRD
              NFIELD=2
              VARNM(1)='utus'
              VARNM(2)='vtus'
              VARNL(1)='eastward stokes transport'
              VARNL(2)='northward stokes transport'
              VARNS(1)='eastward_stokes_transport'
              VARNS(2)='northward_stokes_transport'
              VARNG(1)='eastward_stokes_transport'
              VARNG(2)='northward_stokes_transport'
              VARNC(1)='tus=sqrt(utus**2+vtus**2)'
              VARNC(2)='tus=sqrt(utus**2+vtus**2)'
              VARND(1)=''
              VARND(2)=''
              VMIN = -9900
!
            ELSE IF ( IFI .EQ. 6 .AND. IFJ .EQ. 6 ) THEN
              FSC    = 0.0005
              UNITS  = 'm s-1'
              ENAME  = '.uss'
              DO ISEA=1, NSEA
                 USSX(ISEA)=MAX(-0.9998,MIN(0.9998,USSX(ISEA)))
                 USSY(ISEA)=MAX(-0.9998,MIN(0.9998,USSY(ISEA)))
                 END DO
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, USSX(1:NSEA)      &
                                                        , MAPSF, XX )
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, USSY(1:NSEA)      &
                                                        , MAPSF, XY )
              ENDIF ! SMCGRD
              !! Commented out unnecessary statements below for time being
              !! TAUWIX, TAUWIY are in north-east convention and X1,X2
              !! are not actually written out below
              !DO ISEA=1, NSEA
              !  CABS   = SQRT(USSX(ISEA)**2+USSY(ISEA)**2)
              !  IF ( CABS .NE. UNDEF ) THEN
              !      USSY(ISEA) = MOD ( 630. -                         &
              !            RADE*ATAN2(USSY(ISEA),USSX(ISEA)) , 360. )
              !    ELSE
              !      USSY(ISEA) = UNDEF
              !    END IF
              !  USSX(ISEA) = CABS
              !  END DO
              !CALL W3S2XY ( NSEA, NSEA, NX+1, NY,USSX,MAPSF, X1 )
              !CALL W3S2XY ( NSEA, NSEA, NX+1, NY,USSY,MAPSF, X2 )
              NFIELD=2
              VARNM(1)='uuss'
              VARNM(2)='vuss'
              VARNL(1)='eastward surface stokes drift'
              VARNL(2)='northward surface stokes drift'
              VARNS(1)='eastward_surface_stokes_drift'
              VARNS(2)='northward_surface_stokes_drift'
              VARNG(1)='eastward_surface_stokes_drift'
              VARNG(2)='northward_surface_stokes_drift'
              VARNC(1)='uss=sqrt(uuss**2+vuss**2)'
              WRITE(VARNC(2),'(A,F8.4,A,F8.4,A)') 'Frequency range ',SIG(1)*TPIINV,' to ',SIG(NK)*TPIINV,' Hz'
              VARND(1)=''
              VARND(2)=''
              VMIN = -9900
              VMAX =  9900
!
            ELSE IF ( IFI .EQ. 6 .AND. IFJ .EQ. 7 ) THEN
              NFIELD=2
              FSC    = 0.01
              ENAME  = '.p2s'
              UNITS  = 'm4'
              UNITS2 = 's-1'
              VARNM(1)='fp2s'
              VARNM(2)='pp2s'
              VARNL(1)='power spectral density of equivalent surface pressure'
              VARNL(2)='peak period of power spectral density of equivalent surface pressure'
              VARNS(1)='power_spectral_density_of_equivalent_surface_pressure'
              VARNS(2)='peak_period_of_power_spectral_density_of_equivalent_surface_pressure'
              VARNG(1)='power_spectral_density_of_equivalent_surface_pressure'
              VARNG(2)='peak_period_of_power_spectral_density_of_equivalent_surface_pressure'
              VARNC(1)=''
              VARNC(2)=''
              VARND(1)=''
              VARND(2)=''
              VMIN = -15000
              VMAX = 32000
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PRMS(1:NSEA)      &
                                                        , MAPSF, XX )
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, TPMS(1:NSEA)      &
                                                        , MAPSF, XY )
              ENDIF ! SMCGRD
!
            ELSE IF ( IFI .EQ. 6 .AND. IFJ .EQ. 8 ) THEN
              ! Information for spectral distribution of surface Stokes drift (2nd file)
              FLFRQ=.TRUE.
              NFIELD=2
              VARNM(1)='uusf'
              VARNM(2)='vusf'
              VARNL(1)='eastward spectral variance of surface stokes drift'
              VARNL(2)='northward spectral variance of surface stokes drift'
              VARNS(1)='eastward_spectral_variance_of_surface_stokes_drift'
              VARNS(2)='northward_spectral_variance_of_surface_stokes_drift'
              VARNG(1)='eastward_spectral_variance_of_surface_stokes_drift'
              VARNG(2)='northward_spectral_variance_of_surface_stokes_drift'
              VARNC(1)='usf=sqrt(uusf**2+vusf**2)'
              VARNC(2)='usf=sqrt(uusf**2+vusf**2)'
              VARND(1)=''
              VARND(2)=''
              UNITS   = 'm s-1 Hz-1'
              I1F=US3DF(2)
              I2F=US3DF(3)
              DO IK= I1F,I2F
                IF( SMCGRD ) THEN
                ELSE
                  CALL W3S2XY ( NSEA, NSEA, NX+1,NY,US3D(:,IK), MAPSF, XX )
                  CALL W3S2XY ( NSEA, NSEA, NX+1,NY,US3D(:,NK+IK), MAPSF, XY )
                ENDIF ! SMCGRD
                XXK(:,:,IK)=XX
                XYK(:,:,IK)=XY
              END DO
              FSC    = 0.0005
              ENAME  = '.usf'
 
              VMIN = -9900
              VMAX =  9900
!
            ELSE IF ( IFI .EQ. 6 .AND. IFJ .EQ.  9 ) THEN
                ! Information for spectral microseismic generation data (2nd file)
                FLFRQ=.TRUE.
                NFIELD=1
                FSC    = 0.0004
                VARNM(1)='p2l'
                VARNL(1)='base ten logarithm of power spectral density of equivalent surface pressure'
                VARNS(1)='base_ten_logarithm_of_power_spectral_density_of_equivalent_surface_pressure'
                VARNG(1)='base_ten_logarithm_of_power_spectral_density_of_equivalent_surface_pressure'
                VARNC(1)=''
                VARND(1)=''
                ENAME  = '.p2l'
                I1F=P2MSF(2)
                I2F=P2MSF(3)
                DO IK=I1F,I2F
                  IF( SMCGRD ) THEN
                  ELSE
                    CALL W3S2XY ( NSEA, NSEA, NX+1,NY,P2SMS(:,IK),MAPSF, XX )
                  ENDIF ! SMCGRD
                  IF (NCVARTYPE.EQ.2) THEN
                     WHERE ( XX.GE.0.) XX = ALOG10(XX*(DWAT*GRAV)**2+1E-12)
                     UNITS  = 'log10(Pa2 m2 s+1E-12)'
                  ELSE
                     WHERE ( XX.GE.0.) XX = XX*(DWAT*GRAV)**2
                     UNITS  = 'Pa2 m2 s'
                  END IF
 
                  XK(:,:,IK)=XX
                  END DO
              VMIN = -30000
              VMAX = 30000
!
            ELSE IF ( IFI .EQ. 6 .AND. IFJ .EQ. 10 ) THEN
              FSC    = 0.000001
              UNITS  = 'm2 s-2'
              ENAME  = '.tic'
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, TAUICE(1:NSEA,1)        &
                                                        , MAPSF, XX )
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, TAUICE(1:NSEA,2)        &
                                                        , MAPSF, XY )
              ENDIF ! SMCGRD
              NFIELD=2
              VARNM(1)='utic'
              VARNM(2)='vtic'
              VARNL(1)='eastward wave to sea ice stress'
              VARNL(2)='northward wave to sea ice stress'
              VARNS(1)='eastward_wave_to_sea_ice_stress'
              VARNS(2)='northward_wave_to_sea_ice_stress'
              VARNG(1)='eastward_wave_to_sea_ice_stress'
              VARNG(2)='northward_wave_to_sea_ice_stress'
              VARNC(1)='two=sqrt(utwo**2+vtwo**2)'
              VARNC(2)='two=sqrt(utwo**2+vtwo**2)'
              VARND(1)=''
              VARND(2)=''
              VMIN = -32000
              VMAX =  32000
!
            ELSE IF ( IFI .EQ. 6 .AND. IFJ .EQ. 11 ) THEN
              IF (NCVARTYPEI.EQ.3) NCVARTYPE=4
              FSC    = 0.1
              UNITS  = 'W m-2'
              ENAME  = '.fic'
              DO ISEA=1, NSEA
                 PHIOC(ISEA)=MIN(3000.,PHIOC(ISEA))
                 END DO
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PHICE(1:NSEA)  &
                                                        , MAPSF, X1 )
              ENDIF ! SMCGRD
              NFIELD=1
              VARNM(1)='fic'
              VARNL(1)='wave to sea ice energy flux'
              VARNS(1)='wave_to_sea_ice_energy_flux'
              VARNG(1)='wave_to_sea_ice_energy_flux'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 9990
 
           ELSE IF ( IFI .EQ. 6 .AND. IFJ .EQ. 12 ) THEN
              ! Information for spectral distribution of surface Stokes drift (2nd file)
              FLFRQ=.TRUE.
              IF (USSPF(1)==1) THEN
                 CUSTOMFRQ=.TRUE.
              ENDIF
              NFIELD=2
              VARNM(1)='ussp'
              VARNM(2)='vssp'
              VARNL(1)='eastward partitioned surface stokes drift'
              VARNL(2)='northward partitioned surface stokes drift'
              VARNS(1)='eastward_partitioned_surface_stokes_drift'
              VARNS(2)='northward_partitioned_surface_stokes_drift'
              VARNG(1)='eastward_partitioned_surface_stokes_drift'
              VARNG(2)='northward_partitioned_surface_stokes_drift'
              VARNC(1)='usp=sqrt(ussp**2+vssp**2)'
              VARNC(2)='usp=sqrt(ussp**2+vssp**2)'
              VARND(1)=''
              VARND(2)=''
              UNITS   = 'm s-1'
              I1F=1
              I2F=USSPF(2)
              DO IK= I1F,I2F
                IF( SMCGRD ) THEN
                ELSE
                  CALL W3S2XY ( NSEA, NSEA, NX+1,NY,USSP(:,IK), MAPSF, XX )
                  CALL W3S2XY ( NSEA, NSEA, NX+1,NY,USSP(:,NK+IK), MAPSF, XY )
                ENDIF ! SMCGRD
                XXK(:,:,IK)=XX
                XYK(:,:,IK)=XY
              END DO
              FSC    = 0.0005
              ENAME  = '.usp'
!
            ELSE IF ( IFI .EQ. 7 .AND. IFJ .EQ. 1 ) THEN
              FSC    = 0.01
              ENAME  = '.abr'
              UNITS  = 'm'
              ! NB: ABA and ABD are the X and Y components of the bottom displacement
              IF( SMCGRD ) THEN
              ELSE
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, ABA(1:NSEA)     &
                                                        , MAPSF, XX )
                CALL W3S2XY ( NSEA, NSEA, NX+1, NY, ABD(1:NSEA)     &
                                                        , MAPSF, XY )
              ENDIF ! SMCGRD
              NFIELD=2
              VARNM(1)='uabr'
              VARNM(2)='vabr'
              VARNL(1)='rms of bottom displacement amplitude zonal'
              VARNL(2)='rms of bottom displacement amplitude meridional'
              VARNS(1)='rms_of_bottom_displacement_amplitude_zonal'
              VARNS(2)='rms_of_bottom_displacement_amplitude_meridional'
              VARNG(1)='rms_of_bottom_displacement_amplitude_zonal'
              VARNG(2)='rms_of_bottom_displacement_amplitude_meridional'
              VARNC(1)='abr=sqrt(uabr**2+vabr**2)'
              VARNC(2)='abr=sqrt(uabr**2+vabr**2)'
              VARND(1)=''
              VARND(2)=''
              VMIN = -18000
              VMAX = 18000
!
            ELSE IF ( IFI .EQ. 7 .AND. IFJ .EQ. 2 ) THEN
              FSC    = 0.01
              ENAME  = '.ubr'
              UNITS  = 'm s-1'
              ! NB: UBA and UBD are the X and Y components of the bottom velocity
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, UBA(1:NSEA)       &
                                                        , MAPSF, XX )
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, UBD(1:NSEA)       &
                                                        , MAPSF, XY )
              NFIELD=2
              VARNM(1)='uubr'
              VARNM(2)='vubr'
              VARNL(1)='rms of bottom velocity amplitude zonal'
              VARNL(2)='rms of bottom velocity amplitude meridional'
              VARNS(1)='rms_of_bottom_velocity_amplitude_zonal'
              VARNS(2)='rms_of_bottom_velocity_amplitude_meridional'
              VARNG(1)='rms_of_bottom_velocity_amplitude_zonal'
              VARNG(2)='rms_of_bottom_velocity_amplitude_meridional'
              VARNC(1)='ubr=sqrt(uubr**2+vubr**2)'
              VARNC(2)='ubr=sqrt(uubr**2+vubr**2)'
              VARND(1)=''
              VARND(2)=''
              VMIN = -18000
              VMAX = 18000
!
            ELSE IF ( IFI .EQ. 7 .AND. IFJ .EQ. 3 ) THEN
              FSC    = 0.001
              UNITS  = 'm'
              ENAME  = '.bed'
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, BEDFORMS(1:NSEA,1)    &
                                                        , MAPSF, X1 )
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, BEDFORMS(1:NSEA,2)    &
                                                        , MAPSF, X2 )
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, BEDFORMS(1:NSEA,3)    &
                                                        , MAPSF, XY )
              NFIELD=3
              VARNM(1)='bed'
              VARNM(2)='ripplex'
              VARNM(3)='rippley'
              VARNL(1)='bottom roughness'
              VARNL(2)='eastward sea bottom ripple wavelength'
              VARNL(3)='northward sea bottom ripple wavelength'
              VARNS(1)='sea bottom roughness length'
              VARNS(2)='eastward_ripple_wavelength'
              VARNS(3)='northward_ripple_wavelength'
              VARNG(1)='ripple_wavelength'
              VARNG(2)='eastward_ripple_wavelength'
              VARNG(3)='northward_ripple_wavelength'
              VARNC(1)='ripple_length=sqrt(ripplex**2+rippley**2)'
              VARNC(2)='ripple_length=sqrt(ripplex**2+rippley**2)'
              VARNC(3)='ripple_length=sqrt(ripplex**2+rippley**2)'
              VARND(1)=''
              VARND(2)=''
              VARND(3)=''
              VMIN = 0
              VMAX = 30000
!
            ELSE IF ( IFI .EQ. 7 .AND. IFJ .EQ. 4 ) THEN
              FSC    = 0.1
              UNITS  = 'W m-2'
              ENAME  = '.fbb'
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, PHIBBL(1:NSEA)    &
                                                        , MAPSF, X1 )
              NFIELD=1
              VARNM(1)='fbb'
              VARNL(1)='wave dissipation in bbl'
              VARNS(1)='wave_energy_dissipation_in_bottom_boundary_layer'
              VARNG(1)='wave_dissipation_in_bbl'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 9990
!
            ELSE IF ( IFI .EQ. 7 .AND. IFJ .EQ. 5 ) THEN
              FSC    = 0.000001
              UNITS  = 'm2 s-2'
              ENAME  = '.tbb'
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, TAUBBL(1:NSEA,1) &
                                                        , MAPSF, XX )
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, TAUBBL(1:NSEA,2) &
                                                        , MAPSF, XY )
              NFIELD=2
              VARNM(1)='utbb'
              VARNM(2)='vtbb'
              VARNL(1)='eastward wave to bbl stress'
              VARNL(2)='northward wave to bbl stress'
              VARNS(1)='eastward_wave_to_bottom_boundary_layer_stress'
              VARNS(2)='northward_wave_to_bottom_boundary_layer_stress'
              VARNG(1)='eastward_wave_to_bbl_stress'
              VARNG(2)='northward_wave_to_bbl_stress'
              VARNC(1)='tbb=sqrt(utbb**2+vtbb**2)'
              VARNC(2)='tbb=sqrt(utbb**2+vtbb**2)'
              VARND(1)=''
              VARND(2)=''
              VMIN = -32000
              VMAX =  32000
!
            ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 1 ) THEN
              FSC    = 0.00001
              ENAME  = '.mss'
              UNITS  = '1'
              CALL W3S2XY ( NSEA, NSEA, NX+1,NY,MSSX,MAPSF, XX )
              CALL W3S2XY ( NSEA, NSEA, NX+1,NY,MSSY,MAPSF, XY )
              NFIELD=2
              VARNM(1)='mssu'
              VARNM(2)='mssc'
              VARNL(1)='downwave mean square slope'
              VARNL(2)='crosswave mean square slope'
              VARNS(1)='x_mean_square_slope'
              VARNS(2)='y_mean_square_slope'
              VARNG(1)='x_mean_square_slope'
              VARNG(2)='y_mean_square_slope'
              VARNC(1)='mss=mssu+mssc'
              WRITE(VARNC(1),'(A,F8.4,A,F8.4,A)') 'Frequency range ',SIG(1)*TPIINV,' to ',SIG(NK)*TPIINV,' Hz'
              VARND(1)=''
              VARND(2)=''
              VMIN = 0
              VMAX = 30000
!
            ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 2 ) THEN
              FSC    = 1E-7
              ENAME  = '.msc'
              UNITS  = '1'
              CALL W3S2XY ( NSEA, NSEA, NX+1,NY,MSCX,MAPSF, XX )
              CALL W3S2XY ( NSEA, NSEA, NX+1,NY,MSCY,MAPSF, XY )
              NFIELD=2
              VARNM(1)='mscx'
              VARNM(2)='mscy'
              VARNL(1)='eastward phillips constant'
              VARNL(2)='northward phillips constant'
              VARNS(1)='eastward_phillips_constant'
              VARNS(2)='northward_phillips_constant'
              VARNG(1)='eastward_phillips_constant'
              VARNG(2)='northward_phillips_constant'
              VARNC(1)='msc=mscx+mscy'
              VARNC(2)='msc=mscx+mscy'
              VARND(1)=''
              VARND(2)=''
              VMIN = 0
              VMAX = 30000
!
            ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 3 ) THEN
              FSC    = 0.1
              UNITS  = 'degree'
              ENAME  = '.msd'
              DO ISEA=1, NSEA
                IF ( MSSD(ISEA) .NE. UNDEF )  THEN
                  MSSD(ISEA) = MOD ( 630. - RADE*MSSD(ISEA) , 180. )
                  END IF
                END DO
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, MSSD   , MAPSF, X1 )
              NFIELD=1
              VARNM(1)='mssd'
              VARNL(1)='u direction for mss'
              VARNS(1)='sea_surface_wave_dominant_mean_square_slope_direction'
              VARNG(1)='sea_surface_wave_dominant_mean_square_slope_direction'
              WRITE(VARNC(1),'(A,F8.4,A,F8.4,A)') 'Frequency range ',SIG(1)*TPIINV,' to ',SIG(NK)*TPIINV,' Hz'
              VARND(1)=''
              VMIN = 0
              VMAX = 3600
!
            ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 4 ) THEN
              FSC    = 0.1
              UNITS  = 'degree'
              ENAME  = '.mcd'
              DO ISEA=1, NSEA
                IF ( MSCD(ISEA) .NE. UNDEF )  THEN
                  MSCD(ISEA) = MOD ( 630. - RADE*MSCD(ISEA) , 180. )
                  END IF
                END DO
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, MSCD   , MAPSF, X1 )
              NFIELD=1
              VARNM(1)='mscd'
              VARNL(1)='x direction for msc'
              VARNS(1)='sea_surface_wave_dominant_mean_square_slope_direction_in_highest_frequency'
              VARNG(1)='sea_surface_wave_dominant_mean_square_slope_direction_in_highest_frequency'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 3600
!
             ELSE IF ( IFI .EQ. 8 .AND. IFJ .EQ. 5 ) THEN
              FSC    = 0.001
              UNITS  = '1'
              ENAME  = '.qp'
              IF( SMCGRD ) THEN
              ELSE
                 CALL W3S2XY ( NSEA, NSEA, NX+1, NY, QP, MAPSF, X1 )
              ENDIF
              NFIELD=1
              VARNM(1)='qp'
              VARNL(1)='peakedness'
              VARNS(1)='sea_surface_wave_peakedness'
              VARNG(1)='wave_peakedness'
              VARNC(1)='Goda wave peakedness parameter'
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 9 .AND. IFJ .EQ. 1 ) THEN
              FSC    = 0.1
              UNITS  = 'min.'
              ENAME  = '.dtd'
              DO ISEA=1, NSEA
                IF ( DTDYN(ISEA) .NE. UNDEF ) THEN
                  DTDYN(ISEA) = DTDYN(ISEA) / 60.
                  END IF
                END DO
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, DTDYN , MAPSF, X1 )
              NFIELD=1
              VARNM(1)='dtd'
              VARNL(1)='dynamic time step'
              VARNS(1)='dynamic_time_step'
              VARNG(1)='dynamic_time_step'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 9 .AND. IFJ .EQ. 2 ) THEN
              FSC    = 0.001
              UNITS  = 's-1'
              ENAME  = '.fc'
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, FCUT  , MAPSF, X1 )
              NFIELD=1
              VARNM(1)='fc'
              VARNL(1)='cut off frequency'
              VARNS(1)='cut_off_frequency'
              VARNG(1)='cut_off_frequency'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 8000
!
            ELSE IF ( IFI .EQ. 9 .AND. IFJ .EQ. 3 ) THEN
              FSC    = 0.01
              UNITS  = '1'
              ENAME  = '.cfx'
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, CFLXYMAX  , MAPSF, X1 )
              NFIELD=1
              VARNM(1)='cfx'
              VARNL(1)='maximum cfl for spatial advection'
              VARNS(1)='maximum_cfl_for_spatial_advection'
              VARNG(1)='maximum_cfl_for_spatial_advection'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 9 .AND. IFJ .EQ. 4 ) THEN
              FSC    = 0.01
              UNITS  = '1'
              ENAME  = '.cfd'
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, CFLTHMAX  , MAPSF, X1 )
              NFIELD=1
              VARNM(1)='cfd'
              VARNL(1)='maximum cfl for direction advection'
              VARNS(1)='maximum_cfl_for_direction_advection'
              VARNG(1)='maximum_cfl_for_direction_advection'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 9 .AND. IFJ .EQ. 5 ) THEN
              FSC    = 0.01
              UNITS  = '1'
              ENAME  = '.cfk'
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, CFLKMAX  , MAPSF, X1 )
              NFIELD=1
              VARNM(1)='cfk'
              VARNL(1)='maximum cfl for frequency advection'
              VARNS(1)='maximum_cfl_for_frequency_advection'
              VARNG(1)='maximum_cfl_for_frequency_advection'
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 32000
!
            ELSE IF ( IFI .EQ. 10 ) THEN
              FSC    = 0.1
              UNITS  = 'm'
              WRITE (ENAME,'(A2,I2.2)') '.u', IFJ
              CALL W3S2XY ( NSEA, NSEA, NX+1, NY, USERO(:,IFJ)        &
                                                        , MAPSF, X1 )
              NFIELD=1
              WRITE (VARNM(1),'(A1,I2.2)') 'u', IFJ
              WRITE (VARNL(1),'(A12,I2.2)') 'User_defined', IFJ
              WRITE (VARNS(1),'(A12,I2.2)') 'User_defined', IFJ
              WRITE (VARNG(1),'(A12,I2.2)') 'user_defined', IFJ
              VARNC(1)=''
              VARND(1)=''
              VMIN = 0
              VMAX = 9990
!
            ELSE
              WRITE (NDSE,999) IFI, IFJ
              CALL EXTCDE ( 1 )
!
            END IF ! IFI AND IFJ
 
            ! Defines a differents for the second variable if specified
            UNITVAR(:)=UNITS
            IF (UNITS2.NE.'SAME') UNITVAR(2)=UNITS2
 
 
! 2.2 Make map
 
            ! CB: TODO - need to handle MAPSTA differently for SMC grid output.
            IF( .NOT. SMCGRD ) THEN
            DO IX=1, NX
              DO IY=1, NY
                MAPOUT(IX,IY)=INT2(MAPSTA(IY,IX) + 8*MAPST2(IY,IX))
                IF ( MAPSTA(IY,IX) .EQ. 0 ) THEN
                  X1(IX,IY) = UNDEF
                  X2(IX,IY) = UNDEF
                  XX(IX,IY) = UNDEF
                  XY(IX,IY) = UNDEF
                END IF
                IF ( X1(IX,IY) .EQ. UNDEF ) THEN
                  MAP(IX,IY) = 0
                ELSE
                  MAP(IX,IY) = 1
                END IF
                IF ( X2(IX,IY) .EQ. UNDEF ) THEN
                  MP2(IX,IY) = 0
                ELSE
                  MP2(IX,IY) = 1
                END IF
              END DO
            END DO
            ENDIF ! CB
 
 
! 2.3 Setups the output type 4 ( NetCDF file )
 
            S2=LEN_TRIM(ENAME)
            S1=LEN_TRIM(FILEPREFIX)+S4
            FNAMENC(S1+1:50)='       '
            FNAMENC(S1+1:S1+1) = '_'
 
            ! If flag TOGETHER and not variable with freq dim &
            ! (ef, p2l, ...), no variable name in file name
            IF (TOGETHER.AND.(.NOT.FLFRQ)) THEN
              S2=0
            ! If NOT flag TOGETHER or variable with freq dim &
            ! (ef, p2l, ...), add variable name in file name
            ELSE
              FNAMENC(S1+2:S1+S2) = ENAME(2:S2)
            ENDIF
            ! Defines the netcdf extension
            FNAMENC(S1+S2+1:S1+S2+3) = '.nc'
            FNAMENC(S1+S2+4:S1+S2+6) = '   '
            ! If the flag frequency is .TRUE., defines the fourth dimension
            IF (FLFRQ) THEN
              UNITVAR(:)=UNITS
              DIMLN(4)=I2F-I1F+1
              EXTRADIM=1
            ELSE
              DIMLN(4)=0
              EXTRADIM=0
            END IF
 
            ! If regular grid, initializes the lat/lon or x/y dimension lengths
            IF (GTYPE.NE.UNGTYPE) THEN
              IF( SMCGRD ) THEN
              ELSE ! SMCGRD
                DIMLN(2)=IXN-IX1+1
                DIMLN(3)=IYN-IY1+1
              ENDIF ! SMCGRD
            ! If unstructured mesh, initializes the nelem,tri dimension lengths
            ELSE
              DIMLN(2)=IXN-IX1+1
              DIMLN(3)=NTRI
            ENDIF
 
            ! Defines index of first field variable
            IVAR1=3+EXTRADIM+(COORDTYPE-1)+MAPSTAOUT
 
 
! 2.4.1 Save the id of the previous file
 
            IF (TOGETHER.AND.(.NOT.FLFRQ)) THEN
              OLDNCID = NCIDS(1,1,1)
            ELSE
              OLDNCID = NCIDS(IFI,IFJ,IPART+1)
            END IF
 
 
! 2.4.2 Remove the new file (if not created by the run)
 
            INQUIRE(FILE=FNAMENC, EXIST=FEXIST)
            IF (FEXIST) THEN
              FREMOVE = .FALSE.
              ! time splitted condition
              IF (INDEX(TIMEID,OLDTIMEID).EQ.0) THEN
                ! all variables in the samefile
                IF (TOGETHER.AND.(.NOT.FLFRQ).AND.NCID.EQ.0) FREMOVE = .TRUE.
                ! a file per variable
                IF (.NOT.TOGETHER.OR.FLFRQ) FREMOVE = .TRUE.
              END IF
 
              IF (FREMOVE) THEN
                OPEN(UNIT=1234, IOSTAT=IRET, FILE=FNAMENC, STATUS='old')
                IF (IRET == 0) CLOSE(1234, STATUS='delete')
                FEXIST=.FALSE.
              ELSE
                NCID = OLDNCID
              END IF
            END IF
 
! 2.4.3 Finalize the previous file (if a new one will be created)
 
            IF (.NOT.FEXIST) THEN
              IF (INDEX('0000000000000000',OLDTIMEID).EQ.0 .AND. INDEX(TIMEID,OLDTIMEID).EQ.0) THEN
                IRET = NF90_REDEF(OLDNCID)
                CALL CHECK_ERR(IRET)
                IRET=NF90_PUT_ATT(OLDNCID,NF90_GLOBAL,'stop_date',STRSTOPDATE)
                CALL CHECK_ERR(IRET)
                IRET=NF90_CLOSE(OLDNCID)
                CALL CHECK_ERR(IRET)
              END IF
            END IF
 
 
! 2.5 Creates the netcdf file
 
            IF (.NOT.FEXIST) THEN
 
              ! Initializes the time dimension length
              DIMLN(1)=1
 
              ! If NOT unstructure mesh (i.e. regular grid)
              IF (GTYPE.NE.UNGTYPE) THEN
                ! If spherical coordinate
                IF (FLAGLL) THEN
                  VARNM(NFIELD+1)='Longitude'
                  VARNM(NFIELD+2)='Latitude'
                ! If cartesian coordinate
                ELSE
                  VARNM(NFIELD+1)='x'
                  VARNM(NFIELD+2)='y'
                END IF
              END IF
 
 
 
 
 
              ! Initializes the time iteration counter n
              N=1
 
! 2.5.1 Creates the NetCDF file
              CALL W3CRNC(FNAMENC,NCID,DIMID,DIMLN,VARID, &
                          EXTRADIM,NCTYPE,MAPSTAOUT)
 
              ! Saves the NCID to keep the file opened to write all the variables
              ! and open/close at each time step
              IF (TOGETHER.AND.(.NOT.FLFRQ)) THEN
                NCIDS(1,1,1)=NCID
              ELSE
                NCIDS(IFI,IFJ,IPART+1)=NCID
              END IF
 
              ! If curvilinear grid, instanciates lat / lon
              IF (GTYPE.EQ.CLGTYPE) THEN
                IF (.NOT.ALLOCATED(LON2D)) ALLOCATE(LON2D(NX,NY),LAT2D(NX,NY))
                LON2D=TRANSPOSE(XGRD)
                LAT2D=TRANSPOSE(YGRD)
                IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL, &
                                     'latitude_resolution','n/a')
                CALL CHECK_ERR(IRET)
                IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL, &
                                     'longitude_resolution','n/a')
                CALL CHECK_ERR(IRET)
              ! If NOT curvilinear grid,
              ELSE
                IF( SMCGRD ) THEN
                ELSE ! SMCGRD
                  ! instanciates lon with x/lon for regular grid or nodes for unstructured mesh
                  IF (.NOT.ALLOCATED(LON)) ALLOCATE(LON(NX))
                  IF (.NOT.ALLOCATED(LAT)) THEN
                    ! If regular grid, instanciates lat with y/lat
                    IF (GTYPE.EQ.RLGTYPE) THEN
                      ALLOCATE(LAT(NY))
                    ! If unstructured mesh, instanciates lat with nodes
                    ELSE
                      ALLOCATE(LAT(NX))
                    END IF
                  END IF
                END IF ! SMCGRD
              END IF
 
 
! 2.5.2 Generates Lat-Lon arrays
 
              ! If regular grid
              IF (GTYPE.EQ.RLGTYPE) THEN
                IF( SMCGRD ) THEN
                ELSE ! SMCGRD
                  SXD=DBLE(0.000001d0*DNINT(1d6*(DBLE(SX)) ))
                  SYD=DBLE(0.000001d0*DNINT(1d6*(DBLE(SY)) ))
                  X0D=DBLE(0.000001d0*DNINT(1d6*(DBLE(X0)) ))
                  Y0D=DBLE(0.000001d0*DNINT(1d6*(DBLE(Y0)) ))
                  DO I=1,NX
                    LON(I)=REAL(X0D+SXD*DBLE(I-1))
                  END DO
                  DO I=1,NY
                    LAT(I)=REAL(Y0D+SYD*DBLE(I-1))
                  END DO
                  WRITE(STR2,'(F12.0)') SY
                  STR2=ADJUSTL(STR2)
                  IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,   &
                           'latitude_resolution', TRIM(STR2))
                  CALL CHECK_ERR(IRET)
                  WRITE(STR2,'(F12.0)') SX
                  STR2=ADJUSTL(STR2)
                  IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,   &
                           'longitude_resolution',TRIM(STR2))
                  CALL CHECK_ERR(IRET)
                END IF ! SMCGRD
              END IF
 
              ! If unstructured mesh
              IF (GTYPE.EQ.UNGTYPE) THEN
                LON(:)=XYB(:,1)
                LAT(:)=XYB(:,2)
                IF (.NOT.ALLOCATED(TRIGP2)) ALLOCATE(TRIGP2(3,NTRI))
                DIMLN(2)=NX
                DIMLN(3)=NTRI
                TRIGP2=TRANSPOSE(TRIGP)
                IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL, &
                                     'latitude_resolution','n/a')
                CALL CHECK_ERR(IRET)
                IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL, &
                                     'longitude_resolution','n/a')
                CALL CHECK_ERR(IRET)
              END IF
 
              ! Finishes declaration part in file by adding geographical bounds
              IF(SMCGRD) THEN
                WRITE(STR2,'(F12.0)') MINVAL(LAT)
              ELSE
                WRITE(STR2,'(F12.0)') MINVAL(YGRD)
              ENDIF
              STR2=ADJUSTL(STR2)
              IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,  &
                   'southernmost_latitude',TRIM(STR2))
              CALL CHECK_ERR(IRET)
 
              IF(SMCGRD) THEN
                WRITE(STR2,'(F12.0)') MAXVAL(LAT)
              ELSE
                WRITE(STR2,'(F12.0)') MAXVAL(YGRD)
              ENDIF
              STR2=ADJUSTL(STR2)
              IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,  &
                   'northernmost_latitude',TRIM(STR2))
              CALL CHECK_ERR(IRET)
 
              IF(SMCGRD) THEN
                WRITE(STR2,'(F12.0)') MINVAL(LON)
              ELSE
                WRITE(STR2,'(F12.0)') MINVAL(XGRD)
              ENDIF
              STR2=ADJUSTL(STR2)
              IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,  &
                  'westernmost_longitude',TRIM(STR2))
              CALL CHECK_ERR(IRET)
 
 
              IF(SMCGRD) THEN
                WRITE(STR2,'(F12.0)') MAXVAL(LON)
              ELSE
                WRITE(STR2,'(F12.0)') MAXVAL(XGRD)
              ENDIF
              STR2=ADJUSTL(STR2)
              IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,  &
                  'easternmost_longitude',TRIM(STR2))
              CALL CHECK_ERR(IRET)
              IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,  &
                  'minimum_altitude','-12000 m')
              CALL CHECK_ERR(IRET)
              IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,  &
                  'maximum_altitude','9000 m')
              CALL CHECK_ERR(IRET)
              IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,  &
                  'altitude_resolution','n/a')
              CALL CHECK_ERR(IRET)
 
 
              CALL T2D(TIME,STARTDATE,IERR)
              WRITE(STRSTARTDATE,'(I4.4,A,4(I2.2,A),I2.2)') STARTDATE(1),'-',STARTDATE(2),'-', &
                    STARTDATE(3),' ',STARTDATE(5),':',STARTDATE(6),':',STARTDATE(7)
 
              ! End of define mode of NetCDF file
              IRET = NF90_ENDDEF(NCID)
              CALL CHECK_ERR(IRET)
 
! 2.5.3 Writes longitudes, latitudes, triangles, frequency and status map (mapsta) to netcdf file
 
              ! If regular grid
              IF (GTYPE.EQ.RLGTYPE) THEN
                IF(SMCGRD) THEN ! CB: shelter original code from SMC grid
                ELSE ! SMCGRD
                  IRET=NF90_PUT_VAR(NCID,VARID(1),LON(IX1:IXN))
                  CALL CHECK_ERR(IRET)
                  IRET=NF90_PUT_VAR(NCID,VARID(2),LAT(IY1:IYN))
                  CALL CHECK_ERR(IRET)
                ENDIF ! SMCGRD
              END IF
 
              ! If spherical grid
              IF (GTYPE.EQ.CLGTYPE) THEN
                IRET=NF90_PUT_VAR(NCID,VARID(1),LON2D(IX1:IXN,IY1:IYN))
                CALL CHECK_ERR(IRET)
                IRET=NF90_PUT_VAR(NCID,VARID(2),LAT2D(IX1:IXN,IY1:IYN))
                CALL CHECK_ERR(IRET)
              END IF
 
              ! If unstructured mesh
              IF (GTYPE.EQ.UNGTYPE) THEN
                IRET=NF90_PUT_VAR(NCID,VARID(1),LON(IX1:IXN))
                CALL CHECK_ERR(IRET)
                IRET=NF90_PUT_VAR(NCID,VARID(2),LAT(IX1:IXN))
                CALL CHECK_ERR(IRET)
              END IF
 
              ! Writes frequencies to netcdf file
              IF (EXTRADIM.EQ.1) THEN
                ALLOCATE(FREQ(I2F-I1F+1))
                !BGR Here is where we should tell it what frequencies are.
                IF (CUSTOMFRQ) THEN
                   DO i=1,usspf(2)
                      FREQ(i)=sqrt(GRAV*USSP_WN(i))*TPIINV
                   ENDDO
                ELSE
                   DO i=1,I2F-I1F+1
                      FREQ(i)=SIG(I1F-1+i)*TPIINV
                   END DO
                ENDIF
                IRET=NF90_PUT_VAR(NCID,VARID(3),FREQ)
                CALL CHECK_ERR(IRET)
                DEALLOCATE(FREQ)
              END IF
 
              ! Writes triangles to netcdf file
              IF (GTYPE.EQ.UNGTYPE) THEN
                IRET=NF90_PUT_VAR(NCID,VARID(4+EXTRADIM),TRIGP2)
                CALL CHECK_ERR(IRET)
              END IF
 
              ! Writes status map array at variable index 2+1+coordtype+idim-4
              IF (MAPSTAOUT.EQ.1) THEN
                IDVAROUT=VARID(3+EXTRADIM+(COORDTYPE-1)+MAPSTAOUT)
                START(1)=1
                START(2)=1
                COUNT(1)=IXN-IX1+1
                COUNT(2)=IYN-IY1+1
                IF (GTYPE.NE.UNGTYPE) THEN
                  IRET=NF90_PUT_VAR(NCID,IDVAROUT,MAPOUT(IX1:IXN,IY1:IYN), &
                                     (/START(1:2)/),(/COUNT(1:2)/))
                ELSE
                  IRET=NF90_PUT_VAR(NCID,IDVAROUT,MAPOUT(IX1:IXN,1),(/START(1)/),(/COUNT(1)/))
                ENDIF
                CALL CHECK_ERR(IRET)
              END IF
 
              WRITE (NDSO,973) FNAMENC
 
! 2.5.4  Defines the field(LON,LAT,time) of the variable (i.e. ucur,vcur for current variable)
 
              IRET = NF90_REDEF(NCID)
              CALL CHECK_ERR(IRET)
              DO I=1,NFIELD
                IVAR=IVAR1+I
                IF (COORDTYPE.EQ.1) THEN
                  IF (NCVARTYPE.EQ.2) THEN
                    IF( SMCGRD ) THEN
                    ELSE ! SMCGRD
                      IRET=NF90_DEF_VAR(NCID,VARNM(I), NF90_SHORT, DIMID(2:4+EXTRADIM), VARID(IVAR))
                           IF (NCTYPE.EQ.4) IRET = NF90_DEF_VAR_DEFLATE(NCID, VARID(IVAR), 1, 1, DEFLATE)
                    ENDIF ! SMCGRD
                  ELSE
                    IF( SMCGRD ) THEN
                    ELSE ! SMCGRD
                      IRET=NF90_DEF_VAR(NCID,VARNM(I), NF90_FLOAT, DIMID(2:4+EXTRADIM), VARID(IVAR))
                          IF (NCTYPE.EQ.4) IRET = NF90_DEF_VAR_DEFLATE(NCID, VARID(IVAR), 1, 1, DEFLATE)
                    ENDIF ! SMCGRD
                  END IF
                  CALL CHECK_ERR(IRET)
                ELSE
                  DIMFIELD(1)=DIMID(2)
                  DIMFIELD(2)=DIMID(4)
                  DIMFIELD(3)=DIMID(5)
                  IF (NCVARTYPE.EQ.2) THEN
                    IRET = NF90_DEF_VAR(NCID,VARNM(I), NF90_SHORT, DIMFIELD(1:2+EXTRADIM), VARID(IVAR))
                        IF (NCTYPE.EQ.4) IRET = NF90_DEF_VAR_DEFLATE(NCID, VARID(IVAR), 1, 1, DEFLATE)
                  ELSE
                    IRET = NF90_DEF_VAR(NCID,VARNM(I), NF90_FLOAT, DIMFIELD(1:2+EXTRADIM), VARID(IVAR))
                        IF (NCTYPE.EQ.4) IRET = NF90_DEF_VAR_DEFLATE(NCID, VARID(IVAR), 1, 1, DEFLATE)
                  END IF
                  CALL CHECK_ERR(IRET)
                END IF
!
                IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'long_name',VARNL(I))
                CALL CHECK_ERR(IRET)
                IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'standard_name',VARNS(I))
                CALL CHECK_ERR(IRET)
                IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'globwave_name',VARNG(I))
                CALL CHECK_ERR(IRET)
                IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'units',UNITVAR(I))
                CALL CHECK_ERR(IRET)
                IF (NCVARTYPE.EQ.2) THEN
                  IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'_FillValue',NF90_FILL_SHORT)
                ELSE
                  IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'_FillValue',NF90_FILL_FLOAT)
                  FSC=1.
                  FSC3=1.
                END IF
                CALL CHECK_ERR(IRET)
                IF (I.LE.2) THEN
                  IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'scale_factor',FSC)
                ELSE
                  IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'scale_factor',FSC3)
                ENDIF
                CALL CHECK_ERR(IRET)
                IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'add_offset',0.)
                CALL CHECK_ERR(IRET)
                IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'valid_min',VMIN)
                CALL CHECK_ERR(IRET)
                IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'valid_max',VMAX)
                CALL CHECK_ERR(IRET)
                IF (VARNC(I).NE.'') THEN
                  IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'comment',VARNC(I))
                  CALL CHECK_ERR(IRET)
                END IF
                IF (VARND(I).NE.'') THEN
                  IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'direction_reference',VARND(I))
                  CALL CHECK_ERR(IRET)
                END IF
              END DO
!
              ! put START date in global attribute
              IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,'start_date',STRSTARTDATE)
              CALL CHECK_ERR(IRET)
!
              IRET = NF90_ENDDEF(NCID)
              CALL CHECK_ERR(IRET)
 
 
! 2.6 Append data to the existing file
 
            ELSE  ! FEXIST
 
! 2.6.1 Get the dimensions from the netcdf header
 
              ! If it is an unstructured mesh
              IF (GTYPE.EQ.UNGTYPE) THEN
                IRET=NF90_INQ_VARID (NCID, 'tri', VARID(4+EXTRADIM))
                CALL CHECK_ERR(IRET)
              ! If it is a regular grid
              ELSE
                ! If it is spherical coordinate
                IF (FLAGLL) THEN
                  IF(SMCGRD) THEN
                  ELSE
                    IRET=NF90_INQ_DIMID (NCID, 'longitude', DIMID(2))
                    IRET=NF90_INQ_DIMID (NCID, 'latitude', DIMID(3))
                  ENDIF ! SMCGRD
                  IRET=NF90_INQ_VARID (NCID, 'longitude', VARID(1))
                  IRET=NF90_INQ_VARID (NCID, 'latitude', VARID(2))
                ! If it is cartesian coordinate
                ELSE
                  IRET=NF90_INQ_DIMID (NCID, 'x', DIMID(2))
                  IRET=NF90_INQ_VARID (NCID, 'x', VARID(1))
                  IRET=NF90_INQ_DIMID (NCID, 'y', DIMID(3))
                  IRET=NF90_INQ_VARID (NCID, 'y', VARID(1))
                END IF
                CALL CHECK_ERR(IRET)
              END IF
              ! Get the dimension time
              IRET=NF90_INQ_DIMID (NCID, 'time', DIMID(4+EXTRADIM))
              IRET=NF90_INQUIRE_DIMENSION (NCID, DIMID(4+EXTRADIM),len=N)
              CALL CHECK_ERR(IRET)
              IRET=NF90_INQ_VARID (NCID, 'time', VARID(3+EXTRADIM))
              ! Get the dimension f
              IF (EXTRADIM.EQ.1) IRET=NF90_INQ_DIMID (NCID, 'f', DIMID(4))
 
! 2.6.2 Increments the time step for existing file
 
              ! If it is the first field of the file in mode together
              ! or NOT together or variable with freq dim (ef or p2l)
              ! ChrisBunney: Also - check IPART=TABIPART in case first
              ! requested output is a partitioned field.
              IF((TOGETHER .AND. IFI.EQ.I1 .AND. IFJ.EQ.J1 .AND. IPART.EQ.TABIPART(1))  &
                 .OR.(.NOT.TOGETHER).OR.FLFRQ) n=n+1
 
! 2.6.3 Defines or gets the variables identifiers
 
              ! If it is the first time step, define all the variables and attributes
              IF (N.EQ.1) THEN
                IRET = NF90_REDEF(NCID)
                CALL CHECK_ERR(IRET)
 
                ! Loops on all the fields of the variable (i.e. ucur/vcur for current)
                DO I=1,NFIELD
                  IVAR=IVAR1+I
                  IF (COORDTYPE.EQ.1) THEN
                    IF (NCVARTYPE.EQ.2) THEN
                      IF( SMCGRD ) THEN
                      ELSE
                         IRET = NF90_DEF_VAR(NCID,VARNM(I), NF90_SHORT, DIMID(2:4+EXTRADIM), VARID(IVAR))
                      ENDIF ! SMCGRD
                      IF (NCTYPE.EQ.4) IRET = NF90_DEF_VAR_DEFLATE(NCID, VARID(IVAR), 1, 1, DEFLATE)
                    ELSE
                      IF( SMCGRD ) THEN
                      ELSE
                        IRET = NF90_DEF_VAR(NCID,VARNM(I), NF90_FLOAT, DIMID(2:4+EXTRADIM), VARID(IVAR))
                        CALL CHECK_ERR(IRET)
                      ENDIF ! SMCGRD
                      IF (NCTYPE.EQ.4) IRET = NF90_DEF_VAR_DEFLATE(NCID, VARID(IVAR), 1, 1, DEFLATE)
                      CALL CHECK_ERR(IRET)
                    END IF
                  ELSE
                    DIMFIELD(1)=DIMID(2)
                    DIMFIELD(2)=DIMID(4)
                    DIMFIELD(3)=DIMID(5)
                    IF (NCVARTYPE.EQ.2) THEN
                      IRET = NF90_DEF_VAR(NCID,VARNM(I), NF90_SHORT, DIMFIELD(1:2+EXTRADIM), VARID(IVAR))
                      IF (NCTYPE.EQ.4) IRET = NF90_DEF_VAR_DEFLATE(NCID, VARID(IVAR), 1, 1, DEFLATE)
                    ELSE
                      IRET = NF90_DEF_VAR(NCID,VARNM(I), NF90_FLOAT, DIMFIELD(1:2+EXTRADIM), VARID(IVAR))
                      IF (NCTYPE.EQ.4) IRET = NF90_DEF_VAR_DEFLATE(NCID, VARID(IVAR), 1, 1, DEFLATE)
                    END IF
                  END IF
                  CALL CHECK_ERR(IRET)
!
                  IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'long_name',VARNL(I))
                  IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'standard_name',VARNS(I))
                  CALL CHECK_ERR(IRET)
                  IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'globwave_name',VARNG(I))
                  IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'units',UNITVAR(I))
                  CALL CHECK_ERR(IRET)
!
                  IF (NCVARTYPE.EQ.2) THEN
                    IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'_FillValue',NF90_FILL_SHORT)
                  ELSE
                    IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'_FillValue',NF90_FILL_FLOAT)
                    FSC=1.
                    FSC3=1.
                  END IF
                  CALL CHECK_ERR(IRET)
!
                  IF (I.LE.2) THEN
                    IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'scale_factor',FSC)
                  ELSE
                    IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'scale_factor',FSC3)
                  ENDIF
                  CALL CHECK_ERR(IRET)
!
                  IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'add_offset',0.)
                  CALL CHECK_ERR(IRET)
                  IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'valid_min',VMIN)
                  CALL CHECK_ERR(IRET)
                  IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'valid_max',VMAX)
                  CALL CHECK_ERR(IRET)
                  IF (VARNC(I).NE.'') THEN
                    IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'comment',VARNC(I))
                    CALL CHECK_ERR(IRET)
                  END IF
                  IF (VARND(I).NE.'') THEN
                    IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'direction_reference',VARND(I))
                    CALL CHECK_ERR(IRET)
                  END IF
                END DO
                IRET = NF90_ENDDEF(NCID)
                CALL CHECK_ERR(IRET)
 
              ! If it is not the first time step, get all VARID from the netcdf file opened
              ELSE
                IRET=NF90_REDEF(NCID)
                CALL CHECK_ERR(IRET)
                DO I=1,NFIELD
                  IVAR=IVAR1+I
                  IRET=NF90_INQ_VARID (NCID, VARNM(I), VARID(IVAR))
                  CALL CHECK_ERR(IRET)
                END DO
                IRET=NF90_ENDDEF(NCID)
                CALL CHECK_ERR(IRET)
              END IF !   N.EQ.1
            END IF  ! IERR.EQ.0
 
! 2.6.4 Defines the current time step and index
 
            CALL T2D(TIME,CURDATE,IERR)
            OUTJULDAY=TSUB(REFDATE,CURDATE)
            WRITE(NDSO,'(3A,I6,A,I4,A,I2.2,A,I2.2,A,I2.2,A,I2.2,A,I2.2,2A)')        &
                    'Writing new record ', ENAME(2:) ,'number ',N,                  &
                    ' for ',CURDATE(1),':',CURDATE(2),':',CURDATE(3),'T',CURDATE(5),&
                    ':',CURDATE(6),':',CURDATE(7),' in file ',TRIM(FNAMENC)
 
 
 
            ! Defines starting point and size of arrays to be written
            START(1)=1
            START(2)=1
            START(3)=1
            START(4)=1
 
            ! Sets time index
            START(3+1-COORDTYPE+EXTRADIM)=N
            COUNT(1)=IXN-IX1+1
            COUNT(2)=IYN-IY1+1
            COUNT(3)=1
            COUNT(4)=1
            START1D(1)=1
            START1D(2)=N
            COUNT1D(1)=IXN-IX1+1
            COUNT1D(2)=1
 
            ! Puts time in NetCDF file
            IF((IFI.EQ.I1.AND.IFJ.EQ.J1.AND.TOGETHER)  &
               .OR.(.NOT.TOGETHER).OR.FLFRQ) THEN
              IVAR1=3+EXTRADIM+(COORDTYPE-1)+MAPSTAOUT
              IRET=NF90_PUT_VAR(NCID,VARID(3+EXTRADIM),OUTJULDAY,(/N/))
              CALL CHECK_ERR(IRET)
            END IF
!
! 2.6.5 Puts field(s) in NetCDF file
 
! NFIELD=3
            IF (NCVARTYPE.EQ.2) THEN
              IF ( NFIELD.EQ.3 ) THEN
                IF (SMCGRD) THEN
                ELSE ! IF(SMCGRD)
                  DO IX=IX1, IXN
                    DO IY=IY1, IYN
                      IF ( MAPSTA(IY,IX) .LE. 0 .OR. X1(IX,IY) .EQ. UNDEF ) THEN
                        MXX(IX,IY) = MFILL
                        MYY(IX,IY) = MFILL
                        MXY(IX,IY) = MFILL
                      ELSE
                        MXX(IX,IY) = NINT(X1(IX,IY)/FSC)
                        MYY(IX,IY) = NINT(X2(IX,IY)/FSC)
                        MXY(IX,IY) = NINT(XY(IX,IY)/FSC3)
                      END IF
                    END DO
                  END DO
 
                  IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+1),              &
                          MXX(IX1:IXN,IY1:IYN),(/START(1:3)/),(/COUNT(1:3)/))
                  CALL CHECK_ERR(IRET)
                  IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+2),            &
                          MYY(IX1:IXN,IY1:IYN),(/START(1:3)/),(/COUNT(1:3)/))
                  CALL CHECK_ERR(IRET)
                  IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+3),            &
                          MXY(IX1:IXN,IY1:IYN),(/START(1:3)/),(/COUNT(1:3)/))
                  CALL CHECK_ERR(IRET)
                ENDIF ! SMCGRD
! NFIELD=2
              ELSE IF (NFIELD.EQ.2 ) THEN
! EXTRADIM=0
                IF (EXTRADIM.EQ.0) THEN
                  IF (SMCGRD) THEN
                  ELSE ! IF(SMCGRD)
                    DO IX=IX1, IXN
                      DO IY=IY1, IYN
                        IF ( MAPSTA(IY,IX) .LE. 0 .OR. XX(IX,IY) .EQ. UNDEF ) THEN
                          MXX(IX,IY) = MFILL
                          MYY(IX,IY) = MFILL
                        ELSE
                          MXX(IX,IY) = NINT(XX(IX,IY)/FSC)
                          MYY(IX,IY) = NINT(XY(IX,IY)/FSC)
                        END IF
                      END DO
                    END DO
                    IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+1),             &
                              MXX(IX1:IXN,IY1:IYN),(/START(1:3)/),(/COUNT(1:3)/))
                    CALL CHECK_ERR(IRET)
                    IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+2),           &
                            MYY(IX1:IXN,IY1:IYN),(/START(1:3)/),(/COUNT(1:3)/))
                    CALL CHECK_ERR(IRET)
                  ENDIF ! SMCGRD
! EXTRADIM=1
                ELSE
                  START(3+1-COORDTYPE)=0
                  DO IK=I1F,I2F
                    START(3+1-COORDTYPE)=START(3+1-COORDTYPE)+1
 
                    IF (SMCGRD) THEN
                    ELSE ! IF(SMCGRD)
                      DO IX=IX1, IXN
                        DO IY=IY1, IYN
                          IF ( MAPSTA(IY,IX) .LE. 0 .OR.XXK(IX,IY,IK) .EQ. UNDEF ) THEN
                            MXX(IX,IY) = MFILL
                            MYY(IX,IY) = MFILL
                          ELSE
                            MXX(IX,IY) = NINT(XXK(IX,IY,IK)/FSC)
                            MYY(IX,IY) = NINT(XYK(IX,IY,IK)/FSC)
                          END IF
                        END DO
                      END DO
                      IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+1),               &
                              MXX(IX1:IXN,IY1:IYN),(/START(1:4)/),(/COUNT(1:4)/))
                      IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+2),             &
                              MYY(IX1:IXN,IY1:IYN),(/START(1:4)/),(/COUNT(1:4)/))
                    ENDIF ! SMCGRD
                  END DO
                END IF  ! EXTRADIM
! NFIELD=1
              ELSE
! EXTRADIM=0
                IF (EXTRADIM.EQ.0) THEN
                  IF (SMCGRD) THEN
                  ELSE ! IF(SMCGRD)
                    DO IX=IX1, IXN
                      DO IY=IY1, IYN
                        IF ( MAPSTA(IY,IX) .LE. 0 .OR.X1(IX,IY) .EQ. UNDEF ) THEN
                          MX1(IX,IY) = MFILL
                        ELSE
                          MX1(IX,IY) = NINT(X1(IX,IY)/FSC)
                        END IF
                      END DO
                    END DO
                    IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+1),               &
                            MX1(IX1:IXN,IY1:IYN),(/START(1:3)/),(/COUNT(1:3)/))
                    CALL CHECK_ERR(IRET)
                  ENDIF ! SMCGRD
! EXTRADIM=1
                ELSE
                  START(3+1-COORDTYPE)=0
                  DO IK=I1F,I2F
                    START(3+1-COORDTYPE)=START(3+1-COORDTYPE)+1
 
                    IF (SMCGRD) THEN
                    ELSE ! IF(SMCGRD)
                      DO IX=IX1, IXN
                        DO IY=IY1, IYN
                          IF ( MAPSTA(IY,IX) .LE. 0 .OR.XK(IX,IY,IK) .EQ. UNDEF ) THEN
                            MX1(IX,IY) = MFILL
                          ELSE
                            MX1(IX,IY) = NINT(XK(IX,IY,IK)/FSC)
                          END IF
                        END DO
                      END DO
                      IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+1),               &
                          MX1(IX1:IXN,IY1:IYN),(/START(1:4)/),(/COUNT(1:4)/))
                      CALL CHECK_ERR(IRET)
                    ENDIF ! SMCGRD
                  END DO
                END IF   ! EXTRADIM
              END IF   ! NFIELD
!
! Real output (NCVARTYPE.GE.3)
!
            ELSE
              IF ( NFIELD.EQ.3 ) THEN
                IF (SMCGRD) THEN
                ELSE ! IF(SMCGRD)
                  DO IX=IX1, IXN
                    DO IY=IY1, IYN
                      IF ( MAPSTA(IY,IX) .LE. 0 .OR. X1(IX,IY) .EQ. UNDEF ) THEN
                        MXXR(IX,IY) = MFILLR
                        MYYR(IX,IY) = MFILLR
                        MXYR(IX,IY) = MFILLR
                      ELSE
                        MXXR(IX,IY) = X1(IX,IY)
                        MYYR(IX,IY) = X2(IX,IY)
                        MXYR(IX,IY) = XY(IX,IY)
                      END IF
                    END DO
                  END DO
 
                  IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+1),              &
                          MXXR(IX1:IXN,IY1:IYN),(/START(1:3)/),(/COUNT(1:3)/))
                  CALL CHECK_ERR(IRET)
                  IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+2),            &
                          MYYR(IX1:IXN,IY1:IYN),(/START(1:3)/),(/COUNT(1:3)/))
                  CALL CHECK_ERR(IRET)
                  IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+3),            &
                          MXYR(IX1:IXN,IY1:IYN),(/START(1:3)/),(/COUNT(1:3)/))
                  CALL CHECK_ERR(IRET)
                ENDIF ! SMCGRD
! NFIELD=2
              ELSE IF (NFIELD.EQ.2 ) THEN
! EXTRADIM=0
                IF (EXTRADIM.EQ.0) THEN
                  IF (SMCGRD) THEN
                  ELSE ! IF SMCGRD
                    DO IX=IX1, IXN
                      DO IY=IY1, IYN
                        IF ( MAPSTA(IY,IX) .LE. 0 .OR. XX(IX,IY) .EQ. UNDEF ) THEN
                          MXXR(IX,IY) = MFILLR
                          MYYR(IX,IY) = MFILLR
                        ELSE
                          MXXR(IX,IY) = XX(IX,IY)
                          MYYR(IX,IY) = XY(IX,IY)
                        END IF
                      END DO
                    END DO
                    IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+1),             &
                              MXXR(IX1:IXN,IY1:IYN),(/START(1:3)/),(/COUNT(1:3)/))
                    CALL CHECK_ERR(IRET)
                    IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+2),           &
                            MYYR(IX1:IXN,IY1:IYN),(/START(1:3)/),(/COUNT(1:3)/))
                    CALL CHECK_ERR(IRET)
                  ENDIF ! SMCGRD
  ! EXTRADIM=1
                ELSE
                  START(4-COORDTYPE)=0
                  DO IK=I1F,I2F
                    START(4-COORDTYPE)=START(4-COORDTYPE)+1
 
                    IF (SMCGRD) THEN
                    ELSE ! IF SMCGRD
                      DO IX=IX1, IXN
                        DO IY=IY1, IYN
                          IF ( MAPSTA(IY,IX) .LE. 0 .OR.XXK(IX,IY,IK) .EQ. UNDEF ) THEN
                            MXXR(IX,IY) = MFILLR
                            MYYR(IX,IY) = MFILLR
                          ELSE
                            MXXR(IX,IY) = XXK(IX,IY,IK)
                            MYYR(IX,IY) = XYK(IX,IY,IK)
                          END IF
                        END DO
                      END DO
                      IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+1),               &
                              MXXR(IX1:IXN,IY1:IYN),(/START(1:4)/),(/COUNT(1:4)/))
                      IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+2),             &
                              MYYR(IX1:IXN,IY1:IYN),(/START(1:4)/),(/COUNT(1:4)/))
                    ENDIF ! SMCGRD
                  END DO
                END IF  ! EXTRADIM
! NFIELD=1
              ELSE
! EXTRADIM=0
                IF (EXTRADIM.EQ.0) THEN
                  IF (SMCGRD) THEN
                  ELSE ! IF SMCGRD
                    DO IX=IX1, IXN
                      DO IY=IY1, IYN
                        IF ( MAPSTA(IY,IX) .LE. 0 .OR.X1(IX,IY) .EQ. UNDEF ) THEN
                          MX1R(IX,IY) = MFILLR
                        ELSE
                          MX1R(IX,IY) = X1(IX,IY)
                        END IF
                      END DO
                    END DO
                    IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+1),               &
                            MX1R(IX1:IXN,IY1:IYN),(/START(1:3)/),(/COUNT(1:3)/))
                    CALL CHECK_ERR(IRET)
                  ENDIF ! SMCGRD
! EXTRADIM=1
                ELSE
                  START(4-COORDTYPE)=0
                  DO IK=I1F,I2F
                    START(4-COORDTYPE)=START(4-COORDTYPE)+1
                    IF (SMCGRD) THEN
                    ELSE ! IF SMCGRD
                      DO IX=IX1, IXN
                        DO IY=IY1, IYN
                          IF ( MAPSTA(IY,IX) .LE. 0 .OR.XK(IX,IY,IK) .EQ. UNDEF ) THEN
                            MX1R(IX,IY) = MFILLR
                          ELSE
                            MX1R(IX,IY) = XK(IX,IY,IK)
                          END IF
                        END DO
                      END DO
                      IRET=NF90_PUT_VAR(NCID,VARID(IVAR1+1),               &
                          MX1R(IX1:IXN,IY1:IYN),(/START(1:4)/),(/COUNT(1:4)/))
                      CALL CHECK_ERR(IRET)
                    END IF ! SMCGRD
                  END DO
                END IF   ! EXTRADIM
              END IF   ! NFIELD
            END IF   ! NCVARTYPE
 
            ! updates the variable index
            IVAR1=IVAR1+NFIELD
 
 
            ! Loops over IPART for partition variables
            ! ChrisBunney: Don't loop IPART for last two entries in section 4
            ! (16: total wind sea fraction, 17: number of parts) as these fields
            ! do not have partitions.
            IF (IFI .EQ. 4 .AND. IFJ .LE. NOGE(IFI) - 2) THEN
560           CONTINUE
              IF (INDEXIPART.LT.NBIPART) THEN
                INDEXIPART=INDEXIPART+1
                IF (TABIPART(INDEXIPART).EQ.-1) GOTO 560
                IPART=TABIPART(INDEXIPART)
                GOTO 555
              END IF
            ELSE
              INDEXIPART=1
            END IF
!
          END IF  ! FLG2D(IFI,IFJ)
        END DO  ! IFI=1, NOGRP
      END DO  ! IFJ=1, NGRPP
!
! Clean up
      DEALLOCATE(X1, X2, XX, XY, XK, XXK, XYK)
      DEALLOCATE(MX1, MXX, MYY, MXY, MAPOUT)
      DEALLOCATE(MX1R, MXXR, MYYR, MXYR)
      IF (ALLOCATED(LON)) DEALLOCATE(LON, LAT)
      IF (ALLOCATED(LON2D)) DEALLOCATE(LON2D, LAT2D)
!
      RETURN
!
! Error escape locations
!
 
!
! Formats
!
  973 FORMAT ( 'NEW NETCDF FILE WAS CREATED ',A)
  999 FORMAT (/' *** WAVEWATCH III ERROR IN W3EXNC :'/                &
               '     PLEASE UPDATE FIELDS !!! '/                      &
               '     IFI = ',I2, '- IFJ = ',I2/)
!
!/
 
 
 
!/ End of W3EXNC ----------------------------------------------------- /
!/
      END SUBROUTINE W3EXNC
 
 
 
 
!--------------------------------------------------------------------------
      SUBROUTINE W3CRNC (NCFILE, NCID, DIMID, DIMLN, VARID,  &
                         EXTRADIM, NCTYPE, MAPSTAOUT )
      USE W3GDATMD
      USE NETCDF
      USE W3TIMEMD
 
      IMPLICIT NONE
 
 
 
      INTEGER, INTENT(IN)               :: EXTRADIM
      INTEGER, INTENT(IN)               :: NCTYPE
      CHARACTER*(*), INTENT(IN)         :: NCFILE
      INTEGER, INTENT(OUT)              :: NCID
      INTEGER, INTENT(OUT)              :: DIMID(6)
      INTEGER, INTENT(IN)               :: DIMLN(6)
      INTEGER, INTENT(OUT)              :: VARID(300)
      INTEGER, INTENT(IN)               :: MAPSTAOUT
!
!/ ------------------------------------------------------------------- /
!   Local parameters
!
      INTEGER                           :: IVAR,IRET,ICODE,STRL,STRL2
      INTEGER                           :: DIMTRI(2)
    INTEGER                           :: DEFLATE=1
!
      CHARACTER                         :: ATTNAME*120,ATTVAL*120
 
 
!
! Creation in netCDF3 or netCDF4
!
      IF(NCTYPE.EQ.3)  IRET = NF90_CREATE(TRIM(NCFILE), NF90_CLOBBER, NCID)
      IF(NCTYPE.EQ.4) IRET = NF90_CREATE(TRIM(NCFILE), NF90_NETCDF4, NCID)
      CALL CHECK_ERR(IRET)
!
! Define dimensions
!
      IRET = NF90_DEF_DIM(NCID, 'level', DIMLN(1), DIMID(1))
 
!
! Regular structured case
!
      IF (GTYPE.NE.UNGTYPE) THEN
        IF (FLAGLL) THEN
          IF (SMCGRD) THEN
          ELSE
            IRET = NF90_DEF_DIM(NCID, 'longitude', DIMLN(2), DIMID(2))
            IRET = NF90_DEF_DIM(NCID, 'latitude', DIMLN(3), DIMID(3))
          ENDIF ! SMCGRD
        ELSE
          IRET = NF90_DEF_DIM(NCID, 'x', DIMLN(2), DIMID(2))
          IRET = NF90_DEF_DIM(NCID, 'y', DIMLN(3), DIMID(3))
          END IF
        CALL CHECK_ERR(IRET)
!
! Unstructured case
!
      ELSE
        IRET = NF90_DEF_DIM(NCID, 'node', DIMLN(2), DIMID(2))
        IRET = NF90_DEF_DIM(NCID, 'element', DIMLN(3), DIMID(3))
        CALL CHECK_ERR(IRET)
      ENDIF
!
 
 
      IF (EXTRADIM.EQ.1) THEN
        IRET = NF90_DEF_DIM(NCID, 'f', DIMLN(4), DIMID(4))
        CALL CHECK_ERR(IRET)
      ENDIF
 
      IRET = NF90_DEF_DIM(NCID, 'time',NF90_UNLIMITED, DIMID(4+EXTRADIM))
      CALL CHECK_ERR(IRET)
 
      IF (GTYPE.EQ.UNGTYPE) THEN
        IRET = NF90_DEF_DIM(NCID, 'noel',3, DIMID(5+EXTRADIM))
        CALL CHECK_ERR(IRET)
      ENDIF
 
 
!
!     define variables
!
      IF (FLAGLL) THEN
!longitude
        IF (GTYPE.EQ.RLGTYPE) THEN
          IF (SMCGRD) THEN
          ELSE
            IRET = NF90_DEF_VAR(NCID, 'longitude', NF90_FLOAT, DIMID(2), VARID(1))
            IRET = NF90_DEF_VAR(NCID, 'latitude', NF90_FLOAT, DIMID(3), VARID(2))
          ENDIF ! SMCGRD
        ELSE IF (GTYPE.EQ.CLGTYPE) THEN
          IRET = NF90_DEF_VAR(NCID, 'longitude', NF90_FLOAT, (/ DIMID(2), DIMID(3)/), &
                                                                            VARID(1))
          IRET = NF90_DEF_VAR(NCID, 'latitude', NF90_FLOAT, (/ DIMID(2), DIMID(3)/), &
                                                                            VARID(2))
        ELSE
          IRET = NF90_DEF_VAR(NCID, 'longitude', NF90_FLOAT, DIMID(2), VARID(1))
          IRET = NF90_DEF_VAR(NCID, 'latitude', NF90_FLOAT, DIMID(2), VARID(2))
          END IF
        IRET=NF90_PUT_ATT(NCID,VARID(1),'units','degree_east')
        IRET=NF90_PUT_ATT(NCID,VARID(1),'long_name','longitude')
        IRET=NF90_PUT_ATT(NCID,VARID(1),'standard_name','longitude')
        IRET=NF90_PUT_ATT(NCID,VARID(1),'valid_min',-180.0)
        IRET=NF90_PUT_ATT(NCID,VARID(1),'valid_max',360.)
!
        IRET=NF90_PUT_ATT(NCID,VARID(2),'units','degree_north')
        IRET=NF90_PUT_ATT(NCID,VARID(2),'long_name','latitude')
        IRET=NF90_PUT_ATT(NCID,VARID(2),'standard_name','latitude')
        IRET=NF90_PUT_ATT(NCID,VARID(2),'valid_min',-90.0)
        IRET=NF90_PUT_ATT(NCID,VARID(2),'valid_max',180.)
!
        IF(SMCGRD) THEN
        ELSE
        ENDIF ! SMCGRD
!
      ELSE
        IF (GTYPE.EQ.RLGTYPE) THEN
          IRET = NF90_DEF_VAR(NCID, 'x', NF90_FLOAT, DIMID(2), VARID(1))
          IRET = NF90_DEF_VAR(NCID, 'y', NF90_FLOAT, DIMID(3), VARID(2))
        ELSE IF (GTYPE.EQ.CLGTYPE) THEN
          IRET = NF90_DEF_VAR(NCID, 'x', NF90_FLOAT, (/ DIMID(2), DIMID(3)/), &
                                                                            VARID(1))
          IRET = NF90_DEF_VAR(NCID, 'y', NF90_FLOAT, (/ DIMID(2), DIMID(3)/), &
                                                                            VARID(2))
        ELSE
          IRET = NF90_DEF_VAR(NCID, 'x', NF90_FLOAT, DIMID(2), VARID(1))
          IRET = NF90_DEF_VAR(NCID, 'y', NF90_FLOAT, DIMID(2), VARID(2))
          END IF
!
        IRET=NF90_PUT_ATT(NCID,VARID(1),'units','m')
        IRET=NF90_PUT_ATT(NCID,VARID(1),'long_name','x')
        IRET=NF90_PUT_ATT(NCID,VARID(2),'units','m')
        IRET=NF90_PUT_ATT(NCID,VARID(2),'long_name','y')
!
      END IF  ! FLAGLL
!
      IRET=NF90_PUT_ATT(NCID,VARID(1),'axis','X')
      IRET=NF90_PUT_ATT(NCID,VARID(2),'axis','Y')
        IF (NCTYPE.EQ.4) IRET = NF90_DEF_VAR_DEFLATE(NCID, VARID(1), 1, 1, DEFLATE)
        IF (NCTYPE.EQ.4) IRET = NF90_DEF_VAR_DEFLATE(NCID, VARID(2), 1, 1, DEFLATE)
 
!
! frequency
!
      if (EXTRADIM.EQ.1) THEN
        IRET = NF90_DEF_VAR(NCID, 'f', NF90_FLOAT, DIMID(4), VARID(3))
        IF (NCTYPE.EQ.4) IRET = NF90_DEF_VAR_DEFLATE(NCID, VARID(3), 1, 1, DEFLATE)
        CALL CHECK_ERR(IRET)
        IRET=NF90_PUT_ATT(NCID,VARID(3),'long_name','wave_frequency')
        CALL CHECK_ERR(IRET)
        IRET=NF90_PUT_ATT(NCID,VARID(3),'standard_name','wave_frequency')
        CALL CHECK_ERR(IRET)
        IRET=NF90_PUT_ATT(NCID,VARID(3),'units','s-1')
        CALL CHECK_ERR(IRET)
        IRET=NF90_PUT_ATT(NCID,VARID(3),'axis','Hz')
        CALL CHECK_ERR(IRET)
      END IF
 
 
!
!  time
!
      IRET = NF90_DEF_VAR(NCID, 'time', NF90_DOUBLE, DIMID(4+EXTRADIM), VARID(3+EXTRADIM))
      IF (NCTYPE.EQ.4) IRET = NF90_DEF_VAR_DEFLATE(NCID, VARID(3+EXTRADIM), 1, 1, DEFLATE)
      CALL CHECK_ERR(IRET)
      IRET=NF90_PUT_ATT(NCID,VARID(3+EXTRADIM),'long_name','julian day (UT)')
      CALL CHECK_ERR(IRET)
      IRET=NF90_PUT_ATT(NCID,VARID(3+EXTRADIM),'standard_name','time')
      CALL CHECK_ERR(IRET)
      IRET=NF90_PUT_ATT(NCID,VARID(3+EXTRADIM),'calendar','standard')
      CALL CHECK_ERR(IRET)
      IRET=NF90_PUT_ATT(NCID,VARID(3+EXTRADIM),'units','days since 1990-01-01 00:00:00')
      CALL CHECK_ERR(IRET)
      IRET=NF90_PUT_ATT(NCID,VARID(3+EXTRADIM),'conventions', &
        'relative julian days with decimal part (as parts of the day )')
      IRET=NF90_PUT_ATT(NCID,VARID(3+EXTRADIM),'axis','T')
      CALL CHECK_ERR(IRET)
!
! triangles for irregular grids
!
      IF (GTYPE.EQ.UNGTYPE) THEN
        IVAR=3+EXTRADIM+1
        DIMTRI(1)=DIMID(4+EXTRADIM+1)
        DIMTRI(2)=DIMID(3)
        IRET = NF90_DEF_VAR(NCID, 'tri', NF90_INT, DIMTRI, VARID(IVAR))
        IF (NCTYPE.EQ.4) IRET = NF90_DEF_VAR_DEFLATE(NCID, VARID(IVAR), 1, 1, DEFLATE)
      END IF
!
!  Status map: useful for grid combination
!
      IF (MAPSTAOUT.EQ.1) THEN
        IF (GTYPE.EQ.UNGTYPE) THEN
          IVAR=5+EXTRADIM
          IRET = NF90_DEF_VAR(NCID,'MAPSTA', NF90_SHORT,(/ DIMID(2) /), VARID(IVAR))
        ELSE
          IVAR=4+EXTRADIM
          IRET = NF90_DEF_VAR(NCID,'MAPSTA', NF90_SHORT,(/ DIMID(2) , DIMID(3) /), &
                                                                       VARID(IVAR))
          ENDIF
          IF (NCTYPE.EQ.4) IRET = NF90_DEF_VAR_DEFLATE(NCID, VARID(IVAR), 1, 1, DEFLATE)
!
        IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'long_name','status map')
        IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'standard_name','status map')
        IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'units','1')
        CALL CHECK_ERR(IRET)
        IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'valid_min',-32)
        CALL CHECK_ERR(IRET)
        IRET=NF90_PUT_ATT(NCID,VARID(IVAR),'valid_max',32)
        CALL CHECK_ERR(IRET)
      END IF
 
!
! Global attributes
!
      IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,'WAVEWATCH_III_version_number' ,TRIM(WWVER))
      CALL CHECK_ERR(IRET)
      IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,'WAVEWATCH_III_switches',TRIM(SWITCHES))
      CALL CHECK_ERR(IRET)
    IF (ZZWND.NE.10)      IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,'SIN4 namelist parameter ZWD',ZZWND)
    IF (AALPHA.NE.0.0095) IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,'SIN4 namelist parameter ALPHA0',AALPHA)
    IF (BBETA.NE.1.43)    IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,'SIN4 namelist parameter BETAMAX',BBETA)
    IF(SSDSC(7).NE.0.3)   IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,'SDS4 namelist parameter WHITECAPWIDTH', SSDSC(7))
! ... TO BE CONTINUED ...
 
      IF(SMCGRD) THEN
      ENDIF
 
      open(unit=994,file='NC_globatt.inp',status='old',iostat=ICODE)
      IF (ICODE.EQ.0) THEN
        DO WHILE (ICODE.EQ.0)
          read(994,'(a)',iostat=ICODE) ATTNAME
          read(994,'(a)',iostat=ICODE) ATTVAL
          IF (ICODE.EQ.0) THEN
            STRL=LEN_TRIM(ATTNAME)
            STRL2=LEN_TRIM(ATTVAL)
            IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,ATTNAME(1:STRL),ATTVAL(1:STRL2))
            CALL CHECK_ERR(IRET)
          END IF
        END DO
      ENDIF
      CLOSE(994)
      IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,'product_name' ,TRIM(NCFILE))
      CALL CHECK_ERR(IRET)
      IRET=NF90_PUT_ATT(NCID,NF90_GLOBAL,'area',TRIM(GNAME))
      CALL CHECK_ERR(IRET)
 
      RETURN
 
      END SUBROUTINE W3CRNC
 
!==============================================================================
 
      SUBROUTINE CHECK_ERR(IRET)
 
      USE NETCDF
      USE W3ODATMD, ONLY: NDSE
      USE W3SERVMD, ONLY: EXTCDE
 
      IMPLICIT NONE
 
      INTEGER IRET
 
      IF (IRET .NE. NF90_NOERR) THEN
        WRITE(NDSE,*) ' *** WAVEWATCH III ERROR IN OUNF :'
        WRITE(NDSE,*) ' NETCDF ERROR MESSAGE: '
        WRITE(NDSE,*) NF90_STRERROR(IRET)
        CALL EXTCDE ( 59 )
      END IF
      RETURN
 
      END SUBROUTINE CHECK_ERR
 
!==============================================================================
 
 
!/
!/ End of W3OUNF ----------------------------------------------------- /
!/
      END PROGRAM W3OUNF
 
 
 
 
 
