#include "w3macros.h"
!/ ------------------------------------------------------------------- /
      PROGRAM W3TRNC
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           M. Accensi              |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         11-Apr-2016 |
!/                  +-----------------------------------+
!/
!/    17-Feb-2016 : Creation                            ( version 5.11 )
!/    11-Apr-2016 : Adapted to use more options         ( version 5.11 )
!/
!/    Copyright 2014 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 :
!
!     Convert direct access track output file to netCDF file.
!
!  2. Method :
!
!     Info read from track_o.ww3, written to track.nc
!
!  3. Parameters :
!
!  4. Subroutines used :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!      W3NMOD    Subr. W3GDATMD Set number of model.
!      W3NOUT    Subr. W3ODATMD Set number of model for output.
!     ----------------------------------------------------------------
!
!  5. Called by :
!
!     None, stand-alone program.
!
!  6. Error messages :
!
!  7. Remarks :
!
!  8. Structure :
!
!     See source code.
!
!  9. Switches :
!
!       !/S    Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE CONSTANTS

      USE W3GDATMD, ONLY : W3NMOD, W3SETG, FLAGLL, XFR
      USE W3ODATMD, ONLY : W3NOUT, W3SETO, FNMPRE
      USE W3SERVMD, ONLY : ITRACE, NEXTLN, EXTCDE
!/S      USE W3SERVMD, ONLY : STRACE
      USE W3TIMEMD
!
      USE W3ODATMD, ONLY: NDSO, NDSE, NDST
!
      USE NETCDF
!
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      CHARACTER*34, PARAMETER ::                                      &
                       IDTST  = 'WAVEWATCH III TRACK OUTPUT SPECTRA'
!
      INTEGER                 :: NDSI, NDSINP,                        &
                                 NDSOUT, NDSTRC, NTRACE,              &
                                 NSPEC, IERR, MK, MTH, IT,            &
                                 ILOC, ISPEC, TIME(2), J, S3, HasNC4, &
                                 IOUT, TOUT(2), NOUT, TDUM(2)
!/S      INTEGER, SAVE           :: IENT   = 0
      INTEGER                 :: iret, NCTYPE,NCID, ITH
      CHARACTER(30)           :: FILEPREFIX
      CHARACTER*20            :: FORMAT1
      CHARACTER               :: IDTIME*23, IDDDAY*11, stopdate*30,   &
                                 TRCKID*32
      INTEGER(kind=4)         :: iyc,imc,idc,ihc,iminc,isc

      REAL                    :: TH1, DTH, X, Y, DW, CX, CY, CAO, CDO,&
                                 WX, WY, WAO, WDO, UST, AS, DTEST,    &
                                 DTREQ, DTHD, RTH0
      REAL                    :: M2KM
      REAL, ALLOCATABLE       :: FREQ(:), FREQ1(:), FREQ2(:), DSIP(:),&
                                 SPEC(:,:), E(:,:), THD(:), DIR(:)
      CHARACTER               :: COMSTR*1, IDSTR*34, TSTSTR*3, STIME*23
      INTEGER                 :: dimid(4), varid(18), dimln(4)
!/
!/ ------------------------------------------------------------------- /
!/
!
! 1.a Initialize data structure
!
      CALL W3NMOD ( 1, 6, 6 )
      CALL W3SETG ( 1, 6, 6 )
      CALL W3NOUT (    6, 6 )
      CALL W3SETO ( 1, 6, 6 )
!
! 1.b IO set-up.
!
      NDSI   = 10
      NDSINP = 11
      NDSOUT = 51
!
      NDSTRC =  6
      NTRACE = 10
      CALL ITRACE ( NDSTRC, NTRACE )
!
!/S      CALL STRACE ( IENT, 'W3TRNC' )
!
      WRITE (NDSO,900)
!
      J      = LEN_TRIM(FNMPRE)
      OPEN (NDSI,FILE=FNMPRE(:J)//'ww3_trnc.inp',STATUS='OLD',        &
            ERR=805,IOSTAT=IERR)
      READ (NDSI,'(A)',END=806,ERR=807) COMSTR
      IF (COMSTR.EQ.' ') COMSTR = '$'
      WRITE (NDSO,901) COMSTR

!
! 1.c  Read requests from input file.
!     Output times
!
      CALL NEXTLN ( COMSTR , NDSI , NDSE )
      READ (NDSI,*,END=801,ERR=802) TOUT, DTREQ, NOUT
      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
!
      CALL NEXTLN ( COMSTR , NDSI , NDSE )
      READ (NDSI,*,END=801,ERR=802) NCTYPE

      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
!/NC4      HasNC4=1
      IF ((HasNC4 .eq. 0).and.(NCTYPE.eq.4)) THEN
        WRITE (NDSE,1012)
        CALL EXTCDE ( 1 )
      END IF

      CALL NEXTLN ( COMSTR , NDSI , NDSE )
      FILEPREFIX= 'ww3.'
      READ (NDSI,*,END=801,ERR=802) FILEPREFIX

      ! S3 defines the number of characters in the date for the filename
      ! S3=4-> YYYY, S3=6 -> YYYYMM, S3=10 -> YYYYMMDDTHHZ ...
      CALL NEXTLN ( COMSTR , NDSI , NDSE )
      READ (NDSI,*,END=801,ERR=802) S3

!
! 1.e Check consistency with input file and track_o.ww3
!
      OPEN (NDSINP,FILE=FNMPRE(:J)//'track_o.ww3',FORM='UNFORMATTED', &
            STATUS='OLD',ERR=800,IOSTAT=IERR)
      READ (NDSINP,ERR=801,IOSTAT=IERR) IDSTR, FLAGLL, MK, MTH, XFR
!
      IF ( FLAGLL ) THEN
          M2KM  = 1.
        ELSE
          M2KM  = 1.E-3
        END IF
!
      IF ( IDSTR .NE. IDTST ) GOTO 810

      WRITE (NDSO,902) MK, MTH
      NSPEC  = MK * MTH
      ALLOCATE ( FREQ(MK), FREQ1(MK), FREQ2(MK), DSIP(MK), &
                 SPEC(MK,MTH), E(MK,MTH), THD(MTH), DIR(MTH) )
!
      READ (NDSINP,ERR=801,IOSTAT=IERR) TH1, DTH, FREQ, DSIP

!
!--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! 2.  Time management.
!
      IOUT = 0
      NCID = 0
      WRITE (NDSO,970)
      READ (NDSINP,END=444, ERR=801,IOSTAT=IERR) TIME
      BACKSPACE (NDSINP)


! 2.1 Loops on track_o.ww3 to read the time and data
      DO
        DTEST  = DSEC21 ( TIME , TOUT )

        ! cycle to reach the start time of input file
        IF ( DTEST .LT. 0. ) THEN
          CALL TICK21 ( TOUT , DTREQ )
          CYCLE
        END IF

        IF ( DTEST .GE. 0. ) THEN
          TRCKID=''
          READ (NDSINP,END=444, ERR=801,IOSTAT=IERR) TIME, X, Y, TSTSTR, TRCKID
          IF ( TSTSTR .EQ. 'SEA' ) THEN
            READ (NDSINP,ERR=801,IOSTAT=IERR) DW, CX, CY, WX, WY, UST, &
                                              AS, SPEC
          END IF
          IF ( IERR .EQ. -1 ) THEN
            WRITE (NDSO,944)
            EXIT
          END IF


          IF ( TIME(1).EQ.TOUT(1) .AND. TIME(2).EQ.TOUT(2) ) THEN
            ILOC = ILOC + 1
            IF ( TSTSTR .EQ. 'SEA' ) ISPEC = ISPEC + 1
          ENDIF
          IF ( TIME(1).GT.TOUT(1) .OR. TIME(2).GT.TOUT(2) ) THEN
            CALL STME21 ( TIME , STIME )
            WRITE (NDSO,945) STIME, ILOC, ISPEC
            ILOC    = 1
            ISPEC   = 0
            IF ( TSTSTR .EQ. 'SEA' ) ISPEC = ISPEC + 1
            TOUT(1) = TIME(1)
            TOUT(2) = TIME(2)
          ENDIF
        END IF


! 2.1.1 Increments the global time counter IOUT
        IOUT = IOUT + 1
        CALL STME21 ( TOUT , IDTIME )
        WRITE (NDSO,971) IDTIME


! 2.1.2  Processes the variable value for the time step IOUT
        CALL W3EXNC ( FILEPREFIX, NCTYPE, NCID, S3, stopdate, MK, MTH )


! 2.1.3 Defines the stop date
        iyc=TOUT(1)/10000
        imc=(TOUT(1)-iyc*10000)/100
        idc=TOUT(1)-iyc*10000-100*imc
        ihc=TOUT(2)/10000
        iminc=(TOUT(2)-ihc*10000)/100
        isc=TOUT(2)-ihc*10000-100*iminc
        ! Puts the start and stop date in global attributes
        WRITE(stopdate,'(I4.4,A,4(I2.2,A),I2.2)') &
              iyc,'-',imc,'-',idc,' ',ihc,':',iminc,':',isc 


        IF ( IOUT .GE. NOUT ) EXIT
      END DO


    444 CONTINUE

! 2.2 Closes the netCDF file
      IF (ncid.NE.0) THEN
        iret = nf90_redef(ncid)
        call check_err(iret)
        iret=nf90_put_att(ncid,NF90_GLOBAL,'stop_date',stopdate)
        call check_err(iret)
        iret=nf90_close(ncid)
        call check_err(iret)
      END IF

!
      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 )
!
  805 CONTINUE
      WRITE (NDSE,1004) IERR
      CALL EXTCDE ( 14 )
!
  806 CONTINUE
      WRITE (NDSE,1005) IERR
      CALL EXTCDE ( 15 )
!
  807 CONTINUE
      WRITE (NDSE,1006) IERR
      CALL EXTCDE ( 16 )
!

  810 CONTINUE
      WRITE (NDSE,1010) IDSTR, IDTST
      CALL EXTCDE ( 20 )
!
  888 CONTINUE
      WRITE (NDSO,999)
!
! Formats
!
  900 FORMAT (/15X,'   *** WAVEWATCH III Track output postp. ***   '/ &
               15X,'==============================================='/)
  901 FORMAT ( '  Comment character is ''',A,''''/)
!
  902 FORMAT ( '      Spectral grid size : ',I3,' by ',I3//              &
                   '  Opening file : '/                               &
                   ' -----------------------------------------------')
  940 FORMAT (/'  Output time data : '/                               &
               ' --------------------------------------------------'/ &
               '      First time         : ',A)
  941 FORMAT ( '      Interval           : ',A/                       &
               '      Number of requests : ',I7)
!
  944 FORMAT (/'      End of file reached '/)
!
  945 FORMAT ( '     ',A,' :',I6,' points and',I6,'  spectra.')
!
  970 FORMAT (//'  Generating files '/                                &
               ' --------------------------------------------------')
  971 FORMAT ( '      Files for ',A)
!
  999 FORMAT (/'  End of program '/                                   &
               ' ========================================='/          &
               '         WAVEWATCH III Field output '/)
!
 1000 FORMAT (/' *** WAVEWATCH III ERROR IN W3TRNC : '/               &
               '     ERROR IN OPENING INPUT FILE'/                    &
               '     IOSTAT =',I5/)
!
 1001 FORMAT (/' *** WAVEWATCH III ERROR IN W3TRNC : '/               &
               '     PREMATURE END OF INPUT FILE'/)
!
 1002 FORMAT (/' *** WAVEWATCH III ERROR IN W3TRNC : '/               &
               '     ERROR IN READING FROM INPUT FILE'/               &
               '     IOSTAT =',I5/)
!
 1004 FORMAT (/' *** WAVEWATCH III ERROR IN W3TRNC : '/               &
               '     ERROR IN OPENING INPUT FILE'/                    &
               '     IOSTAT =',I5/)
!
 1005 FORMAT (/' *** WAVEWATCH III ERROR IN W3TRNC : '/               &
               '     ERROR IN READING FROM INPUT FILE'/               &
               '     IOSTAT =',I5/)
!
 1006 FORMAT (/' *** WAVEWATCH III ERROR IN W3TRNC : '/               &
               '     ERROR IN OPENING OUTPUT FILE'/                   &
               '     IOSTAT =',I5/)
!
 1010 FORMAT (/' *** WAVEWATCH III ERROR IN W3TRNC : '/               &
               '     ILLEGAL TYPE, NCTYPE =',I4/)
!
 1012 FORMAT (/' *** WAVEWATCH III ERROR IN W3TRNC : '/               &
               '     NCTYPE = 4 but ww3_trnc is compiled without NC4'/)
!/
!/ Internal subroutine W3EXNC ---------------------------------------- /
!/
      CONTAINS
!/ ------------------------------------------------------------------- /
      SUBROUTINE W3EXNC ( FILEPREFIX, NCTYPE, NCID, S3, stopdate, MK, MTH )
!/
!/                  +-----------------------------------+
!/                  |           M. Accensi              |
!/                  |                        FORTRAN 90 |
!/                  | Last update :          8-Apr-2016 |
!/                  +-----------------------------------+
!/
!/     8-apr-2016 : Creation                            ( version 5.11 )
!/
!  1. Purpose :
!
!     Perform actual track output in NetCDF file.
!
!  3. Parameters :
!
!     Parameter list
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!     Internal parameters
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!      STRACE    Subr. W3SERVMD Subroutine tracing.
!      EXTCDE    Subr.   Id.    Abort program as graceful as possible.
!     ----------------------------------------------------------------
!
!  5. Called by :
!
!     Main program in which it is contained.
!
!  6. Error messages :
!
!       None.
!
!  7. Remarks :
!
!       None.
!
!  8. Structure :
!
!     See source code.
!
!  9. Switches :
!
!       !/S  Enable subroutine tracing.
!       !/T  Enable test output.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE NETCDF
      IMPLICIT NONE

!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
      INTEGER, INTENT(IN)       :: NCTYPE, MK, MTH
      CHARACTER(30), INTENT(IN) :: FILEPREFIX, stopdate
      INTEGER, INTENT(INOUT)    :: NCID, S3
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      INTEGER                 :: S1, S2, S4, S5, NDSDAT, IRET,        &
                                 iyc, imc, idc, ihc, iminc, isc
!/NC4    INTEGER                  :: deflate=1

      CHARACTER*30            :: startdate
      DOUBLE PRECISION        :: timenc

      CHARACTER               :: FNAMENC*50,                          &
                                 ENAME*6
      CHARACTER, SAVE         :: OLDTIMEID*16 = '0000000000000000'
      CHARACTER, SAVE         :: TIMEID*16 = '0000000000000000'
!/S      INTEGER, SAVE           :: IENT   =   0
!/
!/ ------------------------------------------------------------------- /
!/
!
!/S      CALL STRACE (IENT, 'W3EXNC')
!

! 1.1 Sets the date as ISO8601 convention
      ! S3 defines the number of characters in the date for the filename
      ! S3=4-> YYYY, S3=6 -> YYYYMM, S3=10 -> YYYYMMDDHH
      ! Setups min and max date format
      IF (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=>YYYYMMDDHH then filetime='YYYYMMDDTHHMMSSZ'
      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)


! 1.2 Setups the output type 4 ( NetCDF file )

      ENAME='.trck'
      S2=LEN_TRIM(ENAME)
      S1=LEN_TRIM(FILEPREFIX)+S4
      FNAMENC(S1+1:50)='       '
      FNAMENC(S1+1:S1+1) = '_'

      ! add variable name in file name
      FNAMENC(S1+2:S1+S2) = ENAME(2:S2)

      ! Defines the netcdf extension
      FNAMENC(S1+S2+1:S1+S2+3) = '.nc'
      FNAMENC(S1+S2+4:S1+S2+6) = '   '

      ! Defines the dimensions
      dimln(1)=NF90_UNLIMITED ! time
      dimln(2)=MK ! frequency
      dimln(3)=MTH ! direction
      dimln(4)=32  ! string track name length


! 1.3 Gets the netcdf id

      NDSDAT=30
      OPEN (NDSDAT,FILE=FNAMENC,status='new',IOSTAT=iret)
      IF (iret.EQ.0) THEN
        ! close old file
        IF (INDEX('0000000000000000',OLDTIMEID).EQ.0 .AND. INDEX(TIMEID,OLDTIMEID).EQ.0) THEN
          iret = nf90_redef(ncid)
          call check_err(iret)
          iret=nf90_put_att(ncid,NF90_GLOBAL,'stop_date',stopdate)
          call check_err(iret)
          iret=nf90_close(ncid)
          call check_err(iret)
        END IF
        ncid=0
      ELSE
        ncid=ncid
      END IF

                        
! 1.4 Creates the netcdf file

      IF (NCID.EQ.0) THEN 

        ! Initializes the time iteration counter n
        IT      = 0
        ILOC    = 0
        ISPEC   = 0

! 1.4.1 Creates the NetCDF file 

        CALL W3CRNC(nctype,fnamenc,ncid,dimid,dimln,varid)

        ! put start date in global attribute
        iyc=TIME(1)/10000
        imc=(TIME(1)-iyc*10000)/100
        idc=TIME(1)-iyc*10000-100*imc
        ihc=TIME(2)/10000
        iminc=(TIME(2)-ihc*10000)/100
        isc=TIME(2)-ihc*10000-100*iminc
        WRITE(startdate,'(I4.4,A,4(I2.2,A),I2.2)')           &
              iyc,'-',imc,'-',idc,' ',ihc,':',iminc,':',isc
!
        iret=nf90_put_att(ncid,NF90_GLOBAL,'start_date',startdate)
        call check_err(iret)

        ! End of define mode of NetCDF file
        iret = nf90_enddef(ncid)
        call check_err(iret)

        ! Process lower band and higher band frequencies
        freq1(1:MK)=freq(1:MK)-0.5*(freq(1:MK)-(freq(1:MK)/XFR))
        freq2(1:MK)=freq(1:MK)+0.5*(-freq(1:MK)+(freq(1:MK)*XFR))
        freq1(1)=freq(1)
        freq2(MK)=freq(MK)

        ! Converts direction unit in degree
        DTHD=360./MTH
        RTH0=TH1/DTH
        DO ITH=1, MTH
          THD(ITH)=DTHD*(RTH0+REAL(ITH-1))
        END DO
        DIR(1:MTH)=MOD(360-THD(1:MTH),360.)


! 1.4.2 Adds general variables to NetCDF file
        iret=nf90_put_var(ncid,varid(2),FREQ)
        call check_err(iret)

        iret=nf90_put_var(ncid,varid(3),FREQ1)
        call check_err(iret)

        iret=nf90_put_var(ncid,varid(4),FREQ2)
        call check_err(iret)

        iret=nf90_put_var(ncid,varid(5),DSIP)
        call check_err(iret)

        iret=nf90_put_var(ncid,varid(6),DIR)
        call check_err(iret)

        WRITE (NDSO,973) FNAMENC

      END IF  ! IERR.EQ.0


! 1.5 Defines the current time step and index
     
      iyc=TIME(1)/10000
      imc=(TIME(1)-iyc*10000)/100
      idc=TIME(1)-iyc*10000-100*imc
      ihc=TIME(2)/10000
      iminc=(TIME(2)-ihc*10000)/100
      isc=TIME(2)-ihc*10000-100*iminc
      !
      timenc=DBLE(julday(idc,imc,iyc)) - DBLE(julday(1,1,1990))              &
             + DBLE((ihc + DBLE((iminc + DBLE(isc /60.)) /60.)) /24.)
      !
      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 ',IT, & 
              ' for ',iyc,':',imc,':',idc,'T',ihc,':',iminc,':',isc,         &
              ' in file ',TRIM(FNAMENC)


!
! 1.6 Exit from W3EXNC if not sea point
!
      IF ( TSTSTR .NE. 'SEA' ) GOTO 888


!
! 1.6.1 Process speed and direction components
!
      WAO = SQRT ( WX**2 + WY**2 )
      IF ( WAO.GT.1.E-7 ) THEN
        WDO = MOD(270.-ATAN2(WY,WX)*RADE,360.)
      ELSE
        WDO = 0.
      END IF

      CAO = SQRT ( CX**2 + CY**2 )
      IF ( CAO.GT.1.E-7 ) THEN
        CDO = MOD(270.-ATAN2(CY,CX)*RADE,360.)
      ELSE
        CDO = 0.
      END IF
              
!
! 1.7.1 Puts dimensions variables in NetCDF file
!
      IT=IT+1
      IF ( UST .LT. 0. ) UST = -1.0

      ! time
      iret=nf90_put_var(ncid,varid(1),timenc,start=(/IT/))
      call check_err(iret)
      ! longitude
      iret=nf90_put_var(ncid,varid(7),M2KM*X,start=(/IT/))
      call check_err(iret)
      ! latitude
      iret=nf90_put_var(ncid,varid(8),M2KM*Y,start=(/IT/))
      call check_err(iret)



! 1.7.2 Puts fields in NetCDF file


! 1.7.2.a Write spectrum

      iret=nf90_put_var(ncid,varid(9),                               &
           TRANSPOSE(SPEC),start=(/1,1,IT/), count=(/MTH,MK,1/))
       call check_err(iret)

! 1.7.2.b Write the basic stuff

      ! Write DW (depth)
      iret=nf90_put_var(ncid, varid(10),DW ,start=(/IT/))
      call check_err(iret)
      ! Write CAO (current - x direction)
      iret=nf90_put_var(ncid, varid(11),CAO ,start=(/IT/))
      call check_err(iret)   
      ! Write CDO (current - y direction)
      iret=nf90_put_var(ncid,varid(12),CDO ,start=(/IT/))
      call check_err(iret)
      ! Write WAO (wind velocity - x direction)
      iret=nf90_put_var(ncid,varid(13),WAO ,start=(/IT/))
      call check_err(iret)
      ! Write WDO (wind velocity - y direction)
      iret=nf90_put_var(ncid,varid(14),WDO ,start=(/IT/))
      call check_err(iret)
      ! Write UST (friction velocity)
      iret=nf90_put_var(ncid,varid(15),UST,start=(/IT/))
      call check_err(iret)
      ! Write AS (air sea temperature difference)
      iret=nf90_put_var(ncid,varid(16),AS ,start=(/IT/))
      call check_err(iret)
      ! Track name
      iret=nf90_put_var(ncid,varid(18),TRCKID,start=(/1,IT/),count=(/LEN_TRIM(TRCKID),1/))
      call check_err(iret)


!
  888 CONTINUE
!
      RETURN

!
! Formats
!
  973 FORMAT ( 'NEW NetCDF file was created ',A)


!/ End of W3EXNC ----------------------------------------------------- /
!/
      END SUBROUTINE W3EXNC




!--------------------------------------------------------------------------
      SUBROUTINE W3CRNC (nctype,ncfile,ncid,dimid,dimln,varid)

      USE NETCDF

      IMPLICIT NONE

      INTEGER, INTENT(IN)               :: nctype
      CHARACTER*(*), INTENT(IN)         :: ncfile
      INTEGER, INTENT(IN)               :: dimln(:)
      INTEGER, INTENT(OUT)              :: dimid(:), varid(:), ncid
      INTEGER                           :: iret, deflate=1

!
! Creation in netCDF3 or netCDF4
!
      IF(NCTYPE.EQ.3)  iret = nf90_create(trim(ncfile), NF90_CLOBBER, ncid)
!/NC4      IF(NCTYPE.EQ.4) iret = nf90_create(trim(ncfile), NF90_NETCDF4, ncid)
      call check_err(iret)

!
!     Define generals dimensions
!
      iret = nf90_def_dim(ncid, 'time', dimln(1), dimid(1))
      call check_err(iret)
      iret = nf90_def_dim(ncid, 'frequency', dimln(2), dimid(2))
      call check_err(iret)
      iret = nf90_def_dim(ncid, 'direction', dimln(3), dimid(3))
      call check_err(iret)
      iret = nf90_def_dim(ncid, 'string32', dimln(4), dimid(4))
      call check_err(iret)

!
!     define generals variables
!

!  time
      iret=nf90_def_var(ncid, 'time', NF90_DOUBLE, (/dimid(1)/), varid(1))
!/NC4      IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(1), 1, 1, deflate)
      iret=nf90_put_att(ncid,varid(1),'long_name','julian day (UT)')
      iret=nf90_put_att(ncid,varid(1),'standard_name','time')
      iret=nf90_put_att(ncid,varid(1),'units','days since 1990-01-01 00:00:00')
      iret=nf90_put_att(ncid,varid(1),'conventions',                   &
         'Relative julian days with decimal part (as parts of the day)')
      iret=nf90_put_att(ncid,varid(1),'axis','T') 

! frequency
      iret=nf90_def_var(ncid, 'frequency', NF90_FLOAT, (/dimid(2)/),varid(2))
!/NC4      IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(2), 1, 1, deflate)
      iret=nf90_put_att(ncid,varid(2),'long_name','center frequencies for spectra')
      iret=nf90_put_att(ncid,varid(2),'standard_name','frequency')
      iret=nf90_put_att(ncid,varid(2),'units','s-1')
      iret=nf90_put_att(ncid,varid(2),'scale_factor',1.)
      iret=nf90_put_att(ncid,varid(2),'add_offset',0.)
      iret=nf90_put_att(ncid,varid(2),'valid_min',0.)
      iret=nf90_put_att(ncid,varid(2),'valid_max',10.)
      iret=nf90_put_att(ncid,varid(2),'_FillValue',NF90_FILL_FLOAT)
      iret=nf90_put_att(ncid,varid(2),'axis','Y') 

!frequency1
      iret=nf90_def_var(ncid, 'frequency1', NF90_FLOAT, (/dimid(2)/), varid(3))
      call check_err(iret)
!/NC4        IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(3), 1, 1, deflate)
      iret=nf90_put_att(ncid,varid(3),'long_name','frequency of lower band')
      iret=nf90_put_att(ncid,varid(3),'standard_name','frequency_of_lower_band')
      iret=nf90_put_att(ncid,varid(3),'globwave_name','frequency_lower_band')
      iret=nf90_put_att(ncid,varid(3),'units','s-1')
      iret=nf90_put_att(ncid,varid(3),'scale_factor',1.)
      iret=nf90_put_att(ncid,varid(3),'add_offset',0.)
      iret=nf90_put_att(ncid,varid(3),'valid_min',0.)
      iret=nf90_put_att(ncid,varid(3),'valid_max',10.)
      iret=nf90_put_att(ncid,varid(3),'_FillValue',NF90_FILL_FLOAT)
      iret=nf90_put_att(ncid,varid(3),'content','Y')
      iret=nf90_put_att(ncid,varid(3),'associates','frequency')

!frequency2
      iret=nf90_def_var(ncid, 'frequency2', NF90_FLOAT, (/dimid(2)/), varid(4))
      call check_err(iret)
!/NC4        IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(4), 1, 1, deflate)
      iret=nf90_put_att(ncid,varid(4),'long_name','frequency of upper band')
      iret=nf90_put_att(ncid,varid(4),'standard_name','frequency_of_upper_band')
      iret=nf90_put_att(ncid,varid(4),'globwave_name','frequency_upper_band')
      iret=nf90_put_att(ncid,varid(4),'units','s-1')
      iret=nf90_put_att(ncid,varid(4),'scale_factor',1.)
      iret=nf90_put_att(ncid,varid(4),'add_offset',0.)
      iret=nf90_put_att(ncid,varid(4),'valid_min',0.)
      iret=nf90_put_att(ncid,varid(4),'valid_max',10.)
      iret=nf90_put_att(ncid,varid(4),'_FillValue',NF90_FILL_FLOAT)
      iret=nf90_put_att(ncid,varid(4),'content','Y')
      iret=nf90_put_att(ncid,varid(4),'associates','frequency')

! frequency area
      iret=nf90_def_var(ncid, 'frequency_area', NF90_FLOAT,(/dimid(2)/),varid(5))
!/NC4      IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(5), 1, 1, deflate)
      iret=nf90_put_att(ncid,varid(5),'long_name','frequency spectral bin width')
      iret=nf90_put_att(ncid,varid(5),'standard_name','frequency_area')
      iret=nf90_put_att(ncid,varid(5),'units','s-2')
      iret=nf90_put_att(ncid,varid(5),'scale_factor',1.)
      iret=nf90_put_att(ncid,varid(5),'add_offset',0.)
      iret=nf90_put_att(ncid,varid(5),'valid_min',0.)
      iret=nf90_put_att(ncid,varid(5),'valid_max',10.)
      iret=nf90_put_att(ncid,varid(5),'_FillValue',NF90_FILL_FLOAT)
      iret=nf90_put_att(ncid,varid(5),'content','Y')
      iret=nf90_put_att(ncid,varid(5),'associates','frequency')

! direction
      iret=nf90_def_var(ncid, 'direction', NF90_FLOAT, (/dimid(3)/),varid(6))
!/NC4      IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(6), 1, 1, deflate)
      iret=nf90_put_att(ncid,varid(6),'long_name','sea surface wave to direction')
      iret=nf90_put_att(ncid,varid(6),'standard_name','sea_surface_wave_to_direction')
      iret=nf90_put_att(ncid,varid(6),'units','degree')
      iret=nf90_put_att(ncid,varid(6),'scale_factor',1.)
      iret=nf90_put_att(ncid,varid(6),'add_offset',0.)
      iret=nf90_put_att(ncid,varid(6),'valid_min',0.)
      iret=nf90_put_att(ncid,varid(6),'valid_max',360.)
      iret=nf90_put_att(ncid,varid(6),'_FillValue',NF90_FILL_FLOAT)
      iret=nf90_put_att(ncid,varid(6),'axis','Z') 

      IF (FLAGLL) THEN 
!  longitude
        iret=nf90_def_var(ncid, 'longitude', NF90_FLOAT, (/dimid(1)/),varid(7))
!/NC4        IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(7), 1, 1, deflate)
        iret=nf90_put_att(ncid,varid(7),'long_name','longitude')
        iret=nf90_put_att(ncid,varid(7),'standard_name','longitude')
        iret=nf90_put_att(ncid,varid(7),'units','degree_east')
        iret=nf90_put_att(ncid,varid(7),'valid_min',-180.0)
        iret=nf90_put_att(ncid,varid(7),'valid_max',360.)
        iret=nf90_put_att(ncid,varid(7),'_FillValue',NF90_FILL_FLOAT)
        iret=nf90_put_att(ncid,varid(7),'content','T')
        iret=nf90_put_att(ncid,varid(7),'associates','time')


!  latitude
        iret=nf90_def_var(ncid, 'latitude', NF90_FLOAT, (/dimid(1)/),varid(8))
!/NC4        IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(8), 1, 1, deflate)
        iret=nf90_put_att(ncid,varid(8),'long_name','latitude')
        iret=nf90_put_att(ncid,varid(8),'standard_name','latitude')
        iret=nf90_put_att(ncid,varid(8),'units','degree_north')
        iret=nf90_put_att(ncid,varid(8),'valid_min',-90.0)
        iret=nf90_put_att(ncid,varid(8),'valid_max',180.)
        iret=nf90_put_att(ncid,varid(8),'_FillValue',NF90_FILL_FLOAT)
        iret=nf90_put_att(ncid,varid(8),'content','T')
        iret=nf90_put_att(ncid,varid(8),'associates','time')
      ELSE
!  longitude
        iret=nf90_def_var(ncid, 'x', NF90_FLOAT, (/dimid(1)/),varid(7))
!/NC4        IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(7), 1, 1, deflate)
        iret=nf90_put_att(ncid,varid(7),'long_name','x')
        iret=nf90_put_att(ncid,varid(7),'standard_name','x')
        iret=nf90_put_att(ncid,varid(7),'units','m')
        iret=nf90_put_att(ncid,varid(7),'_FillValue',NF90_FILL_FLOAT)
        iret=nf90_put_att(ncid,varid(7),'content','T')
        iret=nf90_put_att(ncid,varid(7),'associates','time')

!  latitude
        iret=nf90_def_var(ncid, 'y', NF90_FLOAT, (/dimid(1)/),varid(8))
!/NC4        IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(8), 1, 1, deflate)
        iret=nf90_put_att(ncid,varid(8),'long_name','y')
        iret=nf90_put_att(ncid,varid(8),'standard_name','y')
        iret=nf90_put_att(ncid,varid(8),'units','m')
        iret=nf90_put_att(ncid,varid(8),'_FillValue',NF90_FILL_FLOAT)
        iret=nf90_put_att(ncid,varid(8),'content','T')
        iret=nf90_put_att(ncid,varid(8),'associates','time')

      END IF


! Efth
      iret=nf90_def_var(ncid,'efth',NF90_FLOAT,(/dimid(3),dimid(2),dimid(1)/),varid(9))
!/NC4            IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(9), 1, 1, deflate)
      iret=nf90_put_att(ncid,varid(9),'long_name',                     &
         'sea surface wave directional variance spectral density')
      iret=nf90_put_att(ncid,varid(9),'standard_name',                 &
         'sea_surface_wave_directional_variance_spectral_density')
      iret=nf90_put_att(ncid,varid(9),'globwave_name',                 &
         'directional_variance_spectral_density')
      iret=nf90_put_att(ncid,varid(9),'units','m2 s rad-1')
      iret=nf90_put_att(ncid,varid(9),'scale_factor',1.)
      iret=nf90_put_att(ncid,varid(9),'add_offset',0.)
      iret=nf90_put_att(ncid,varid(9),'valid_min',0.)
      iret=nf90_put_att(ncid,varid(9),'valid_max',10.)
      iret=nf90_put_att(ncid,varid(9),'_FillValue',NF90_FILL_FLOAT)
      iret=nf90_put_att(ncid,varid(9),'content','TYZ')
      iret=nf90_put_att(ncid,varid(9),'associates','time frequency direction')

! DW - depth
      iret=nf90_def_var(ncid, 'dpt', NF90_FLOAT, (/dimid(1)/),varid(10))
!/NC4      IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(10), 1, 1, deflate)
      iret=nf90_put_att(ncid,varid(10),'long_name','depth')
      iret=nf90_put_att(ncid,varid(10),'standard_name','depth')
      iret=nf90_put_att(ncid,varid(10),'globwave_name','depth')
      iret=nf90_put_att(ncid,varid(10),'units','m')
      iret=nf90_put_att(ncid,varid(10),'scale_factor',1.)
      iret=nf90_put_att(ncid,varid(10),'add_offset',0.)
      iret=nf90_put_att(ncid,varid(10),'_FillValue',NF90_FILL_FLOAT)
      iret=nf90_put_att(ncid,varid(10),'content','T')
      iret=nf90_put_att(ncid,varid(10),'associates','time')

! CAO - current speed (m/s)
      iret=nf90_def_var(ncid, 'cur', NF90_FLOAT,(/dimid(1)/), varid(11))
!/NC4      IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(11), 1, 1, deflate)
      iret=nf90_put_att(ncid,varid(11),'long_name','sea water speed')
      iret=nf90_put_att(ncid,varid(11),'standard_name','sea_water_speed')
      iret=nf90_put_att(ncid,varid(11),'globwave_name','sea_water_speed')
      iret=nf90_put_att(ncid,varid(11),'units','m s-1')
      iret=nf90_put_att(ncid,varid(11),'scale_factor',1.)
      iret=nf90_put_att(ncid,varid(11),'add_offset',0.)
      iret=nf90_put_att(ncid,varid(11),'_FillValue',NF90_FILL_FLOAT)
      iret=nf90_put_att(ncid,varid(11),'content','T')
      iret=nf90_put_att(ncid,varid(11),'associates','time')

! CDO - current direction (degree)
      iret=nf90_def_var(ncid, 'curdir', NF90_FLOAT,(/dimid(1)/), varid(12))
!/NC4      IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(12), 1, 1, deflate)
      iret=nf90_put_att(ncid,varid(12),'long_name','direction from of sea water velocity')
      iret=nf90_put_att(ncid,varid(12),'standard_name','direction_of_sea_water_velocity')
      iret=nf90_put_att(ncid,varid(12),'globwave_name','direction_of_sea_water_velocity')
      iret=nf90_put_att(ncid,varid(12),'units','degree')
      iret=nf90_put_att(ncid,varid(12),'scale_factor',1.)
      iret=nf90_put_att(ncid,varid(12),'add_offset',0.)
      iret=nf90_put_att(ncid,varid(12),'_FillValue',NF90_FILL_FLOAT)
      iret=nf90_put_att(ncid,varid(12),'content','T')
      iret=nf90_put_att(ncid,varid(12),'associates','time')

! WAO - wind speed (m/s)
      iret=nf90_def_var(ncid, 'wnd', NF90_FLOAT,(/dimid(1)/), varid(13))
!/NC4      IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(13), 1, 1, deflate)
      iret=nf90_put_att(ncid,varid(13),'long_name','wind speed at 10m')
      iret=nf90_put_att(ncid,varid(13),'standard_name','wind_speed')
      iret=nf90_put_att(ncid,varid(13),'globwave_name','wind_speed')
      iret=nf90_put_att(ncid,varid(13),'units','m s-1')
      iret=nf90_put_att(ncid,varid(13),'scale_factor',1.)
      iret=nf90_put_att(ncid,varid(13),'add_offset',0.)
      iret=nf90_put_att(ncid,varid(13),'_FillValue',NF90_FILL_FLOAT)
      iret=nf90_put_att(ncid,varid(13),'content','T')
      iret=nf90_put_att(ncid,varid(13),'associates','time')

! WDO - wind direction (degree)
      iret=nf90_def_var(ncid, 'wnddir', NF90_FLOAT,(/dimid(1)/), varid(14))
!/NC4      IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(14), 1, 1, deflate)
      iret=nf90_put_att(ncid,varid(14),'long_name','wind direction')
      iret=nf90_put_att(ncid,varid(14),'standard_name','wind_from_direction')
      iret=nf90_put_att(ncid,varid(14),'globwave_name','wind_from_direction')
      iret=nf90_put_att(ncid,varid(14),'units','m s-1')
      iret=nf90_put_att(ncid,varid(14),'scale_factor',1.)
      iret=nf90_put_att(ncid,varid(14),'add_offset',0.)
      iret=nf90_put_att(ncid,varid(14),'_FillValue',NF90_FILL_FLOAT)
      iret=nf90_put_att(ncid,varid(14),'content','T')
      iret=nf90_put_att(ncid,varid(14),'associates','time')

! UST - friction velocity  (m/s)
      iret=nf90_def_var(ncid, 'ust', NF90_FLOAT,(/dimid(1)/), varid(15))
!/NC4      IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(15), 1, 1, deflate)
      iret=nf90_put_att(ncid,varid(15),'long_name','friction velocity')
      iret=nf90_put_att(ncid,varid(15),'standard_name','friction_velocity')
      iret=nf90_put_att(ncid,varid(15),'globwave_name','friction_velocity')
      iret=nf90_put_att(ncid,varid(15),'units','m s-1')
      iret=nf90_put_att(ncid,varid(15),'scale_factor',1.)
      iret=nf90_put_att(ncid,varid(15),'add_offset',0.)
      iret=nf90_put_att(ncid,varid(15),'_FillValue',NF90_FILL_FLOAT)
      iret=nf90_put_att(ncid,varid(15),'content','T')
      iret=nf90_put_att(ncid,varid(15),'associates','time')

! AS - air-sea temperature difference (deg C)
      iret=nf90_def_var(ncid, 'ast',NF90_FLOAT,(/dimid(1)/), varid(16))
!/NC4      IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(16), 1, 1, deflate)
      iret=nf90_put_att(ncid,varid(16),'long_name','air sea temperature difference')
      iret=nf90_put_att(ncid,varid(16),'standard_name','air_sea_temperature_difference')
      iret=nf90_put_att(ncid,varid(16),'globwave_name','air_sea_temperature_difference')
      iret=nf90_put_att(ncid,varid(16),'units','degree')
      iret=nf90_put_att(ncid,varid(16),'scale_factor',1.)
      iret=nf90_put_att(ncid,varid(16),'add_offset',0.)
      iret=nf90_put_att(ncid,varid(16),'_FillValue',NF90_FILL_FLOAT)
      iret=nf90_put_att(ncid,varid(16),'content','T')
      iret=nf90_put_att(ncid,varid(16),'associates','time')

!  string32
      iret=nf90_def_var(ncid, 'string32', NF90_INT, (/dimid(4)/), varid(17))
      call check_err(iret)
!/NC4      IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(17), 1, 1, deflate)
      iret=nf90_put_att(ncid,varid(17),'long_name','track_name number of characters')
      iret=nf90_put_att(ncid,varid(17),'_FillValue',NF90_FILL_INT)
      iret=nf90_put_att(ncid,varid(17),'axis','W') 

!  track_name
      iret=nf90_def_var(ncid, 'track_name', NF90_CHAR, (/dimid(4),dimid(1)/), varid(18))
      call check_err(iret)
!/NC4      IF (NCTYPE.EQ.4) iret=nf90_def_var_deflate(ncid, varid(18), 1, 1, deflate)
      iret=nf90_put_att(ncid,varid(18),'long_name','track name')
      iret=nf90_put_att(ncid,varid(18),'content','TX')
      iret=nf90_put_att(ncid,varid(18),'associates','time string16')


!
      RETURN


      END SUBROUTINE W3CRNC 


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

      subroutine check_err(iret)
   
      USE NETCDF

      implicit none

      integer iret

      if (iret .ne. NF90_NOERR) then
         print *, nf90_strerror(iret)
         stop
      endif
      return
      END subroutine check_err

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

!/
!/ End of W3TRNC ----------------------------------------------------- /
!/
      END PROGRAM W3TRNC

  


