!/ ------------------------------------------------------------------- /
      PROGRAM WW3_SYSTRK
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |     A. J. van der Westhuysen      |
!/                  |            Jeff Hanson            |
!/                  |        Eve-Marie Devaliere        |
!/                  |                        FORTRAN 95 |
!/                  | Last update :         16-Jan-2017 |
!/                  +-----------------------------------+
!/
!/    03-Feb-2012 : Origination, based on Matlab code   ( version 4.05 )
!/                  by Jeff Hanson & Eve-Marie Devaliere
!/    04-Jan-2013 : Inclusion in trunk                  ( version 4.08 )
!/    29-Nov-2013 : Remove DOC control characters,
!/                  update MPI! to MPI/! (H.L. Tolman). ( version 4.15 )
!/    11-Feb-2014 : Add NetCDF output option. Both NetCDF-3 and
!/                  NetCDF-4 are available. (B. Li).    ( version 4.18 )
!/    26-Sep-2016 : Optimization updates (A. van der Westhuysen)
!/                                                      ( version 5.15 )
!/    20-Sep-2016 : Add support for unformatted partition file.
!/                  (S.Zieger BoM, Australia)           ( version 5.16 )
!/    20-Dec-2016 : Optimized search algorithms and
!/                  set functions. (S.Zieger)           ( version 5.16 )
!/
!/    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.
!/
      USE W3STRKMD
!/F90      USE W3TIMEMD, ONLY: TDIFF
      IMPLICIT NONE
!/MPI
!/MPI      INCLUDE "mpif.h"
!
!  1. Purpose :
!
!     Perform spatial and temporal tracking of wave systems, based 
!     on spectral partition (bulletin) output.
!
!  2. Method :
!
!     This is a controller program. It reads the input parameter file 
!     ww3_systrk.inp and calls subroutine waveTracking_NWS_V2 to 
!     perform the actual tracking procedure. Write output (fields and 
!     point output).
!
!  3. Parameters :
!
      LOGICAL      :: testout
      PARAMETER (testout = .FALSE.)
      CHARACTER    :: filename*80, paramFile*32
      REAL         :: dirKnob, perKnob, hsKnob, wetPts, seedLat, &
                      seedLon, dirTimeKnob, tpTimeKnob, tint
      REAL         :: lonout(100), latout(100)                            !Increase dimension?
      INTEGER      :: maxGroup, ntint, noutp
!/F90      INTEGER      :: CLKDT0(8),CLKDT1(8)
!/F90      REAL         :: CLKFEL
      TYPE(dat2d), POINTER :: wsdat(:)
      TYPE(timsys), POINTER :: sysA(:)
      INTEGER, POINTER :: maxSys(:)
!
!     Local parameters.
!     ----------------------------------------------------------------
!     intype         Int       input  Type of input (0 = from memory; 1 = from file)
!     tmax           Int       input  Value of maxTs to apply (1 or 2, used for model coupling)
!     tcur           Int       input  Index of current time step (1 or 2, used for model coupling)
!     ulimGroup      Int       input  Upper limit of number of wave systems to output
!
      LOGICAL           :: file_exists
      CHARACTER         :: inpstr*72
      INTEGER           :: intype, tmax, tcur, maxI, maxJ
      INTEGER           :: it, igrp, sysmatch, ind, ip
      INTEGER           :: i, j, leng, ulimGroup
      REAL, ALLOCATABLE :: dum(:,:)
!/TRKNC      REAL, ALLOCATABLE :: dum2nc(:,:,:,:)
!/TRKNC      REAL, ALLOCATABLE ::  hsprt_nc(:,:,:)
!/TRKNC      REAL, ALLOCATABLE ::  tpprt_nc(:,:,:)
!/TRKNC      REAL, ALLOCATABLE ::  dirprt_nc(:,:,:)
!/TRKNC      REAL, ALLOCATABLE ::  longitude_nc(:),latitude_nc(:)
!/TRKNC      REAL, ALLOCATABLE ::  lonprt_nc(:),latprt_nc(:)
      INTEGER NTIME_NC
      INTEGER           :: outputType
      LOGICAL           :: outputCheck1,outputCheck2
      DOUBLE PRECISION  :: date1, date2, tstart, tend
      REAL              :: dlon, dlat, lonprt, latprt
      REAL              :: dt
      REAL              :: minlon, maxlon, minlat, maxlat
      INTEGER           :: mxcwt, mycwt
!/MPI      INTEGER           :: rank, nproc, ierr
!/MPI      CHARACTER         :: rankstr*4

!     For point output (bilinear interpolation)
      REAL :: hsprt(10),tpprt(10),dirprt(10)
      REAL :: BL_hsprt(10),BR_hsprt(10),TR_hsprt(10),TL_hsprt(10), &
              BL_tpprt(10),BR_tpprt(10),TR_tpprt(10),TL_tpprt(10), &
              BL_dirprt(10),BR_dirprt(10),TR_dirprt(10),TL_dirprt(10)
      REAL :: BL_dirx,BR_dirx,TR_dirx,TL_dirx, &
              BL_diry,BR_diry,TR_diry,TL_diry
      REAL :: BL_lonprt,BR_lonprt,TR_lonprt,TL_lonprt, &
              BL_latprt,BR_latprt,TR_latprt,TL_latprt
      REAL :: t, u, BL_W, BR_W, TR_W, TL_W
      REAL      :: PI
      PARAMETER  (PI = 3.1416) 
!
!  4. Subroutines used :
!
!     waveTracking_NWS_V2
!
!  5. Called by :
!
!     None, stand-alone program.
!
!  6. Error messages :
!
!  7. Remarks :
!
!  8. Structure :
!
!     Calls subroutine waveTracking_NWS_V2 in trackmd.95 - see that
!     file for structure.
!
!  9. Switches :
!
!       !/SHRD  Switch for shared / distributed memory architecture.
!       !/MPI   Id.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
!
!/MPI/!     Start of parallel region
!/MPI      CALL MPI_INIT(ierr)
!/MPI    
!/MPI      CALL MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
!/MPI      CALL MPI_COMM_SIZE(MPI_COMM_WORLD, nproc, ierr)
!/MPI
!     Open log file
!/MPI      WRITE(rankstr,'(i4.4)')  rank
!/MPI      OPEN(unit=20,file='sys_log'//rankstr//'.ww3',status='unknown')
!/SHRD      OPEN(unit=20,file='sys_log.ww3',status='unknown')

!     Print code version
!/MPI      IF (rank.EQ.0) THEN
         WRITE(6,900)
!/MPI      END IF
      WRITE(20,900)
  900 FORMAT (/15X,'    *** WAVEWATCH III Wave system tracking ***  '/ &
               15X,'==============================================='/)

!     Since this program reads the raw partitioning input from file,
!     we set intype=1 or 2, and tmax and tcur to dummy values (not used).
      intype = 2
!      intype = 1
      IF (intype.EQ.1) WRITE(6,*) &
         '*** WARNING: partRes format input used!'
      tmax = 0
      tcur = 0

!     Read input parameter file
!/MPI      IF (rank.EQ.0) THEN
      INQUIRE(FILE='ww3_systrk.inp', EXIST=file_exists)
      IF (.NOT.file_exists) THEN
         WRITE(20,2000)
         WRITE(6,2000)
         CALL ABORT
      END IF
      OPEN(unit=10,file='ww3_systrk.inp',status='old')

      READ(10,'(A72)') inpstr
      DO WHILE (inpstr(1:1).EQ.'$')
         READ(10,'(A72)') inpstr
      END DO
      BACKSPACE(10)
      READ(10,*) filename

      READ(10,'(A72)') inpstr
      DO WHILE (inpstr(1:1).EQ.'$')
         READ(10,'(A72)') inpstr
      END DO
      BACKSPACE(10)
      READ(10,*) date1, date2, dt, ntint
      tstart = date1 + date2/1000000

      READ(10,'(A72)') inpstr
      DO WHILE (inpstr(1:1).EQ.'$')
         READ(10,'(A72)') inpstr
      END DO
      BACKSPACE(10)
      READ(10,*) outputType

      !Check for correct outputType option:
      IF (outputType.EQ.1) THEN
         !ASCII output
      ELSEIF (outputType.EQ.3) THEN
         !NetCDF 3 - requrires !/TRKNC switch
         outputCheck1 = .TRUE.
!/TRKNC     outputCheck1 = .FALSE.
         IF(outputCheck1)  THEN
            WRITE(6,993)
            STOP 
         END IF
      ELSEIF (outputType.EQ.4) THEN
         !NetCDF 4 - requrires !/TRKNC and !/NC4 switch
         outputCheck1 = .TRUE.
         outputCheck2 = .TRUE.
!/TRKNC     outputCheck1 = .FALSE.
!/NC4       outputCheck2 = .FALSE.
         IF(outputCheck1.OR.outputCheck2) THEN
            WRITE(6,994)
            STOP
         END IF
      ELSE
         !Not a valid outputType
         WRITE(6,995) outputType
         STOP
      ENDIF

      READ(10,'(A72)') inpstr
      DO WHILE (inpstr(1:1).EQ.'$')
         READ(10,'(A72)') inpstr
      END DO
      BACKSPACE(10)
      READ(10,*) minlon, maxlon, mxcwt

      READ(10,'(A72)') inpstr
      DO WHILE (inpstr(1:1).EQ.'$')
         READ(10,'(A72)') inpstr
      END DO
      BACKSPACE(10)
      READ(10,*) minlat, maxlat, mycwt

      READ(10,'(A72)') inpstr
      DO WHILE (inpstr(1:1).EQ.'$')
         READ(10,'(A72)') inpstr
      END DO
      BACKSPACE(10)
      READ(10,*) dirKnob, perKnob, hsKnob, wetPts, & 
                dirTimeKnob, tpTimeKnob

      READ(10,'(A72)') inpstr
      DO WHILE (inpstr(1:1).EQ.'$')
         READ(10,'(A72)') inpstr
      END DO
      BACKSPACE(10)
      READ(10,*) seedLat, seedLon

      READ(10,'(A72)') inpstr
      DO WHILE (inpstr(1:1).EQ.'$')
         READ(10,'(A72)') inpstr
      END DO
      BACKSPACE(10)
      noutp = 1
      lonout(:) = 9999.
      latout(:) = 9999.
      DO WHILE (.TRUE.)
         READ(10,*) lonout(noutp),latout(noutp)
         IF ((lonout(noutp).EQ.0.).AND.(latout(noutp).EQ.0.)) EXIT
         noutp = noutp + 1
      END DO
      noutp = noutp - 1 

      CLOSE(10)

      WRITE(20,*) 'Raw partition file = ',filename
      WRITE(20,'(A,F15.6)') 'Start time = ',tstart
      WRITE(20,*) 'dt = ',dt
      WRITE(20,*) 'No. time levels = ',ntint
      WRITE(20,'(A,2F7.2)') 'Domain limits: Longitude =',minlon, maxlon
      WRITE(20,'(A,2F7.2)') '               Latitude  =',minlat, maxlat
      WRITE(20,*) 'No. increments: Long, Lat  =',mxcwt, mycwt
      WRITE(20,*) 'dirKnob, perKnob, hsKnob, wetPts, & 
                dirTimeKnob, tpTimeKnob, seedLat, seedLon ='
      WRITE(20,'(8F6.2)') dirKnob, perKnob, hsKnob, wetPts, & 
                dirTimeKnob, tpTimeKnob, seedLat, seedLon
      WRITE(20,*) 'No. output points =',noutp
      DO i = 1,noutp
         WRITE(20,*) lonout(i), latout(i)
      END DO

      INQUIRE(FILE=filename, EXIST=file_exists)
      IF (.NOT.file_exists) THEN
         WRITE(20,2200) filename
         WRITE(6,2200) filename
         CALL EXIT(1)
      END IF


!/F90         CALL DATE_AND_TIME ( VALUES=CLKDT0 )

!/MPI      END IF

!/MPI/!     MPI communication block 
!/MPI      CALL MPI_BCAST(filename,80,MPI_CHARACTER,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(tstart,1,MPI_DOUBLE_PRECISION,0, &
!/MPI                     MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(tend,1,MPI_DOUBLE_PRECISION,0, &
!/MPI                     MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(dt,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(ntint,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(minlon,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(maxlon,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(minlat,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(maxlat,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(mxcwt,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(mycwt,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(dirKnob,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(perKnob,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(hsKnob,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(wetPts,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(dirTimeKnob,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(tpTimeKnob,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(seedLon,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(seedLat,1,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(noutp,1,MPI_INTEGER,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(lonout,100,MPI_REAL,0,MPI_COMM_WORLD,IERR)
!/MPI      CALL MPI_BCAST(latout,100,MPI_REAL,0,MPI_COMM_WORLD,IERR)

!/MPI      CALL MPI_Barrier(MPI_COMM_WORLD,IERR)

      CALL waveTracking_NWS_V2 (intype     ,tmax       , &
                                tcur       ,filename   , &
                                tstart     ,tend       , &
                                dt         ,ntint      , &
                                minlon     ,maxlon     , & 
                                minlat     ,maxlat     , &
                                mxcwt      ,mycwt      , &
                                dirKnob    ,             &
                                perKnob    ,hsKnob     , &
                                wetPts     ,seedLat    , &
                                seedLon    ,dirTimeKnob, &
                                tpTimeKnob ,paramFile  , &
                                sysA       ,wsdat      , &
                                maxSys     ,maxGroup   )

!/MPI      IF (rank.EQ.0) THEN

!/F90         CALL DATE_AND_TIME ( VALUES=CLKDT1 )
!/F90         CLKFEL = TDIFF ( CLKDT0,CLKDT1 )
!/F90         WRITE (6,998) CLKFEL
         WRITE (6,*) 'Final system output...'

!     Set upper limit for wave systems to output (limited by AWIPS display)
      ulimGroup = 9

!-----Output systems as plain text----------------------------------------

      maxI = SIZE(wsdat(1)%lon,1)
      maxJ = SIZE(wsdat(1)%lon,2)
      dlon = wsdat(1)%lon(2,2)-wsdat(1)%lon(1,1)
      dlat = wsdat(1)%lat(2,2)-wsdat(1)%lat(1,1)
      WRITE(20,*) 'dlon, dlat =',dlon,dlat

!-----Final SYSTEM output: Coordinates
      OPEN(unit=21,file='sys_coord.ww3', status='unknown')

      WRITE(21,'(I6,69X,A)') maxJ,'Number of rows'
      WRITE(21,'(I6,69X,A)') maxI,'Number of cols'
!/TRKNC      ALLOCATE( longitude_nc(maxI) )
!/TRKNC      ALLOCATE( latitude_nc(maxJ) )

      WRITE(21,*) 'Longitude ='
      DO j = maxJ,1,-1
         DO i = 1,maxI
            WRITE(21,'(F7.2)',ADVANCE='NO') wsdat(1)%lon(i,j)
!/TRKNC        longitude_nc(i)=wsdat(1)%lon(i,1)
         END DO
         WRITE(21,'(A)',ADVANCE='YES') ''
      END DO

      WRITE(21,*) 'Latitude = '
      DO j = maxJ,1,-1
         DO i = 1,maxI
            WRITE(21,'(F7.2)',ADVANCE='NO') wsdat(1)%lat(i,j)
!/TRKNC        latitude_nc(j)=wsdat(1)%lat(1,j)
         END DO
         WRITE(21,'(A)',ADVANCE='YES') ''
      END DO

      CLOSE(21)

!-----Final SYSTEM output: hs
      IF(outputType == 1) THEN
      OPEN(unit=22,file='sys_hs.ww3', status='unknown')

      WRITE(22,'(I6,69X,A)') maxJ,'Number of rows'
      WRITE(22,'(I6,69X,A)') maxI,'Number of cols'
      ENDIF

      NTIME_NC=SIZE(sysA)
      ALLOCATE( dum(maxI,maxJ) )
!/TRKNC      IF(outputType == 3 .OR. outputType == 4) THEN
!/TRKNC      ALLOCATE( dum2nc(maxI,maxJ,maxGroup,NTIME_NC) )
!/TRKNC      ENDIF

      DO it = 1,SIZE(sysA)
!        Loop through identified groups, limiting the output in file to ulimGroup    
         IF(outputType == 1) THEN
            WRITE(22,'(F15.6,60x,A)') wsdat(it)%date,'Time'
            WRITE(22,'(I6,69x,A)') MIN(ulimGroup,maxGroup), &
               'Tot number of systems'
         ENDIF 
         DO igrp = 1,MIN(ulimGroup,maxGroup)
            dum(1:maxI,1:maxJ) = 9999.00
!           Find system with this group tag
            sysmatch = 1
            DO WHILE (sysmatch.LE.maxSys(it))
               IF (sysA(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
               sysmatch = sysmatch+1
            END DO
            IF (sysmatch.LE.maxSys(it)) THEN
!              Match found: fill the output matrix with this data
               leng = sysA(it)%sys(sysmatch)%nPoints
               DO ind = 1, leng
                  dum(sysA(it)%sys(sysmatch)%i(ind), &
                      sysA(it)%sys(sysmatch)%j(ind)) = &
                      sysA(it)%sys(sysmatch)%hs(ind)
               END DO
            ELSE
               leng = 0
            END IF
            
      IF(outputType == 1) THEN
            WRITE(22,'(I6,69x,A)') igrp,'System number'  
            WRITE(22,'(I6,69x,A)') leng,'Number of points in system'

            DO J = maxJ,1,-1
               DO i = 1,maxI
                  WRITE(22,'(F8.2)',ADVANCE='NO') dum(i,j)
               END DO
               WRITE(22,'(A)',ADVANCE='YES') ''
            END DO
      ELSE
!/TRKNC            DO J = maxJ,1,-1
!/TRKNC               DO i = 1,maxI
!/TRKNC               dum2nc(i,j,igrp,it)=dum(i,j)
!/TRKNC               END DO
!/TRKNC            END DO
      ENDIF

         END DO
      END DO
        
!/TRKNC      IF(outputType == 3 .OR. outputType == 4 ) THEN
!/TRKNC      call t2netcdf(longitude_nc,latitude_nc,dum2nc,maxI,maxJ,&
!/TRKNC         maxGroup,date1,date2,dt,NTIME_NC,1,outputType)
!/TRKNC      ENDIF

      IF(outputType.EQ.1) CLOSE(22)

!-----Final SYSTEM output: tp
      IF(outputType == 1) THEN
      OPEN(unit=23,file='sys_tp.ww3',status='unknown')

      WRITE(23,'(I6,69X,A)') maxJ,'Number of rows'
      WRITE(23,'(I6,69X,A)') maxI,'Number of cols'
      ENDIF

      DO it = 1,SIZE(sysA)
!        Loop through identified groups, limiting the output in file to ulimGroup
         IF(outputType == 1) THEN
            WRITE(23,'(F15.6,60x,A)') wsdat(it)%date,'Time'
            WRITE(23,'(I6,69X,A)') MIN(ulimGroup,maxGroup), &
               'Tot number of systems'
         ENDIF 
         DO igrp = 1,MIN(ulimGroup,maxGroup)
            dum(1:maxI,1:maxJ) = 9999.00
!           Find system with this group tag
            sysmatch = 1
            DO WHILE (sysmatch.LE.maxSys(it))
               IF (sysA(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
               sysmatch = sysmatch+1
            END DO
            IF (sysmatch.LE.maxSys(it)) THEN
!              Match found: fill the output matrix with this data
               leng = sysA(it)%sys(sysmatch)%nPoints
               DO ind = 1, leng
                  dum(sysA(it)%sys(sysmatch)%i(ind), &
                      sysA(it)%sys(sysmatch)%j(ind)) = &
                      sysA(it)%sys(sysmatch)%tp(ind)
               END DO
            ELSE
               leng = 0
            END IF
            
      IF(outputType == 1) THEN
            WRITE(23,'(I6,69X,A)') igrp,'System number'  
            WRITE(23,'(I6,69X,A)') leng,'Number of points in system'
            DO J = maxJ,1,-1
               DO i = 1,maxI
                  WRITE(23,'(F8.2)',ADVANCE='NO') dum(i,j)
               END DO
               WRITE(23,'(A)',ADVANCE='YES') ''
            END DO
      ELSE

!/TRKNC            DO J = maxJ,1,-1
!/TRKNC               DO i = 1,maxI
!/TRKNC               dum2nc(i,j,igrp,it)=dum(i,j)
!/TRKNC            END DO
!/TRKNC            END DO
      ENDIF

         END DO
      END DO

!/TRKNC      IF(outputType.EQ.3 .OR. outputType.EQ. 4 ) THEN
!/TRKNC      call t2netcdf(longitude_nc,latitude_nc,dum2nc,maxI,maxJ,&
!/TRKNC         maxGroup,date1,date2,dt,NTIME_NC,2,outputType)
!/TRKNC      ENDIF

      IF(outputType.EQ.1) CLOSE(23)

!-----Final SYSTEM output: dir
      IF(outputType == 1) THEN
      OPEN(unit=24,file='sys_dir.ww3',status='unknown')

      WRITE(24,'(I6,69X,A)') maxJ,'Number of rows'
      WRITE(24,'(I6,69X,A)') maxI,'Number of cols'
      ENDIF

      DO it = 1,SIZE(sysA)
!        Loop through identified groups, limiting the output in file to
!        ulimGroup
         IF(outputType == 1) THEN
             WRITE(24,'(F15.6,60x,A)') wsdat(it)%date,'Time'
             WRITE(24,'(I6,69X,A)') MIN(ulimGroup,maxGroup), &
               'Tot number of systems'
         ENDIF 
         DO igrp = 1,MIN(ulimGroup,maxGroup)
            dum(1:maxI,1:maxJ) = 9999.00
!           Find system with this group tag
            sysmatch = 1
            DO WHILE (sysmatch.LE.maxSys(it))
               IF (sysA(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
               sysmatch = sysmatch+1
            END DO
            IF (sysmatch.LE.maxSys(it)) THEN
!              Match found: fill the output matrix with this data
               leng = sysA(it)%sys(sysmatch)%nPoints
               DO ind = 1, leng
                  dum(sysA(it)%sys(sysmatch)%i(ind), &
                      sysA(it)%sys(sysmatch)%j(ind)) = & 
                      sysA(it)%sys(sysmatch)%dir(ind)
               END DO
            ELSE
               leng = 0
            END IF
            
      IF(outputType == 1) THEN
            WRITE(24,'(I6,69X,A)') igrp,'System number'  
            WRITE(24,'(I6,69X,A)') leng,'Number of points in system'
            DO J = maxJ,1,-1
               DO i = 1,maxI
                  WRITE(24,'(F8.2)',ADVANCE='NO') dum(i,j)
               END DO
               WRITE(24,'(A)',ADVANCE='YES') ''
            END DO
      ELSE
!/TRKNC            DO J = maxJ,1,-1
!/TRKNC               DO i = 1,maxI
!/TRKNC               dum2nc(i,j,igrp,it)=dum(i,j)
!/TRKNC               END DO
!/TRKNC            END DO
      END IF

         END DO
      END DO

!/TRKNC      IF(outputType.EQ.3 .OR. outputType.EQ.4 ) THEN
!/TRKNC      call t2netcdf(longitude_nc,latitude_nc,dum2nc,maxI,maxJ,&
!/TRKNC         maxGroup,date1,date2,dt,NTIME_NC,3,outputType)
!/TRKNC      ENDIF
      IF(outputType.EQ.1) CLOSE(24)

!-----Final SYSTEM output: dspr
      IF(outputType == 1) THEN
      OPEN(unit=25,file='sys_dspr.ww3',status='unknown')

      WRITE(25,'(I6,69X,A)') maxJ,'Number of rows'
      WRITE(25,'(I6,69X,A)') maxI,'Number of cols'
      ENDIF

      DO it = 1,SIZE(sysA)
!        Loop through identified groups, limiting the output in file to ulimGroup
         IF(outputType == 1) THEN
            WRITE(25,'(F15.6,60x,A)') wsdat(it)%date,'Time'
            WRITE(25,'(I6,69X,A)') MIN(ulimGroup,maxGroup), &
               'Tot number of systems'
         ENDIF
         DO igrp = 1,MIN(ulimGroup,maxGroup)
            dum(1:maxI,1:maxJ) = 9999.00
!           Find system with this group tag
            sysmatch = 1
            DO WHILE (sysmatch.LE.maxSys(it))
               IF (sysA(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
               sysmatch = sysmatch+1
            END DO
            IF (sysmatch.LE.maxSys(it)) THEN
!              Match found: fill the output matrix with this data
               leng = sysA(it)%sys(sysmatch)%nPoints
               DO ind = 1, leng
                  dum(sysA(it)%sys(sysmatch)%i(ind), &
                      sysA(it)%sys(sysmatch)%j(ind)) = & 
                      sysA(it)%sys(sysmatch)%dspr(ind)
               END DO
            ELSE
               leng = 0
            END IF
            
      IF(outputType == 1) THEN
            WRITE(25,'(I6,69X,A)') igrp,'System number'  
            WRITE(25,'(I6,69X,A)') leng,'Number of points in system'
            DO J = maxJ,1,-1
               DO i = 1,maxI
                  WRITE(25,'(F8.2)',ADVANCE='NO') dum(i,j)
               END DO
               WRITE(25,'(A)',ADVANCE='YES') ''
            END DO
       ELSE
!/TRKNC            DO J = maxJ,1,-1
!/TRKNC               DO i = 1,maxI
!/TRKNC               dum2nc(i,j,igrp,it)=dum(i,j)
!/TRKNC               END DO
!/TRKNC            END DO
       ENDIF

         END DO
      END DO

!/TRKNC      IF(outputType.EQ.3 .OR. outputType.EQ.4 ) THEN
!/TRKNC      call t2netcdf(longitude_nc,latitude_nc,dum2nc,maxI,maxJ,&
!/TRKNC         maxGroup,date1,date2,dt,NTIME_NC,4,outputType)
!/TRKNC      ENDIF
      IF(outputType.EQ.1) CLOSE(25)

      IF (ALLOCATED(DUM)) DEALLOCATE(dum)
!/TRKNC      IF (ALLOCATED(dum2nc)) DEALLOCATE(dum2nc)

!/TRKNC      IF(outputType.EQ.3.OR.outputType.EQ.4) THEN
!/TRKNC      ALLOCATE( hsprt_nc(10,noutp,NTIME_NC) )
!/TRKNC      ALLOCATE( tpprt_nc(10,noutp,NTIME_NC) )
!/TRKNC      ALLOCATE( dirprt_nc(10,noutp,NTIME_NC) )
!/TRKNC      ALLOCATE( lonprt_nc(noutp) )
!/TRKNC      ALLOCATE( latprt_nc(noutp) )
!/TRKNC      ENDIF

!-----Final SYSTEM output: point output
      IF(outputType == 1) THEN
      OPEN(unit=26,file='sys_pnt.ww3',status='unknown')
      WRITE(26,'(A)') '%'
      WRITE(26,'(A)') '%'
      WRITE(26,'(A)') '% WW3 Wave tracking point output'
      WRITE(26,'(A)') '%'
      WRITE(26,'(10A)') '%       Xp            Yp            ', &
           'HsSY01        HsSY02        HsSY03        HsSY04        ', &
           'HsSY05        HsSY06        HsSY07        HsSY08        ', &
           'HsSY09        HsSY10        ', &
           'TpSY01        TpSY02        TpSY03        TpSY04        ', &
           'TpSY05        TpSY06        TpSY07        TpSY08        ', &
           'TpSY09        TpSY10        ', &
           'DrSY01        DrSY02        DrSY03        DrSY04        ', &
           'DrSY05        DrSY06        DrSY07        DrSY08        ', &
           'DrSY09        DrSY10'   
      WRITE(26,'(10A)') '%       [degr]        [degr]        ', &
           '[m]           [m]           [m]           [m]           ', &
           '[m]           [m]           [m]           [m]           ', &
           '[m]           [m]           ', &
           '[sec]         [sec]         [sec]         [sec]         ', &
           '[sec]         [sec]         [sec]         [sec]         ', &
           '[sec]         [sec]         ', &
           '[degr]        [degr]        [degr]        [degr]        ', &
           '[degr]        [degr]        [degr]        [degr]        ', &
           '[degr]        [degr]'
      WRITE(26,'(A)') '%'
       ENDIF

      DO it = 1,SIZE(sysA)
      IF(outputType == 1) THEN
         WRITE(26,'(A,F15.6)') 'Time : ',wsdat(it)%date
       ENDIF

         DO ip = 1,noutp
            hsprt(1:10) = 999.9999
            tpprt(1:10) = 999.9999
            dirprt(1:10) = 999.9999
            lonprt = 999.9999
            latprt = 999.9999
            BL_hsprt(1:10) = 999.9999
            BL_tpprt(1:10) = 999.9999
            BL_dirprt(1:10) = 999.9999
            BR_hsprt(1:10) = 999.9999
            BR_tpprt(1:10) = 999.9999
            BR_dirprt(1:10) = 999.9999
            TL_hsprt(1:10) = 999.9999
            TL_tpprt(1:10) = 999.9999
            TL_dirprt(1:10) = 999.9999
            TR_hsprt(1:10) = 999.9999
            TR_tpprt(1:10) = 999.9999
            TR_dirprt(1:10) = 999.9999
            BL_lonprt = 999.9999
            BL_latprt = 999.9999
            BR_lonprt = 999.9999
            BR_latprt = 999.9999
            TL_lonprt = 999.9999
            TL_latprt = 999.9999
            TR_lonprt = 999.9999
            TR_latprt = 999.9999
            BL_W = 999
            BR_W = 999
            TR_W = 999
            TL_W = 999

            DO j = 1, (maxJ-1)
               DO i = 1, (maxI-1)
                  IF ( ( ((lonout(ip).GE. &
                           wsdat(1)%lon(i,j)).AND. &
                          (lonout(ip).LT. &
                           wsdat(1)%lon(i+1,j))).OR. &
                         ((lonout(ip).GT. &
                          wsdat(1)%lon(i,j)).AND. &
                         (lonout(ip).LE. &
                          wsdat(1)%lon(i+1,j))) ).AND. &
                       ( ((latout(ip).GE. &
                          wsdat(1)%lat(i,j)).AND. &
                          (latout(ip).LT. &
                          wsdat(1)%lat(i,j+1))).OR. &
                         ((latout(ip).GT. &
                          wsdat(1)%lat(i,j)).AND. &
                          (latout(ip).LE. &
                          wsdat(1)%lat(i,j+1))) ) ) &
                  THEN
                     BL_lonprt = wsdat(1)%lon(i,j)
                     BL_latprt = wsdat(1)%lat(i,j)
                     BR_lonprt = wsdat(1)%lon(i+1,j)
                     BR_latprt = wsdat(1)%lat(i+1,j)
                     TL_lonprt = wsdat(1)%lon(i,j+1)
                     TL_latprt = wsdat(1)%lat(i,j+1)
                     TR_lonprt = wsdat(1)%lon(i+1,j+1)
                     TR_latprt = wsdat(1)%lat(i+1,j+1)
!                    Compute weights for this point
                     t = (lonout(ip)-BL_lonprt)/(BR_lonprt-BL_lonprt)
                     u = (latout(ip)-BL_latprt)/(TL_latprt-BL_latprt)
                     BL_W = (1-t)*(1-u)
                     BR_W = t*(1-u)
                     TR_W = t*u
                     TL_W = (1-t)*u
!                    Compute output values using weights
                     lonprt = BL_W*BL_lonprt + BR_W*BR_lonprt + &
                              TL_W*TL_lonprt + TR_W*TR_lonprt
                     latprt = BL_W*BL_latprt + BR_W*BR_latprt + &
                              TL_W*TL_latprt + TR_W*TR_latprt
                  END IF
               END DO
            END DO
!           Loop through identified groups, limiting the output in file to 10
            DO igrp = 1,MIN(10,maxGroup)
!              Find system with this group tag
               sysmatch = 1
               DO WHILE (sysmatch.LE.maxSys(it))
                  IF (sysA(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
                  sysmatch = sysmatch+1
               END DO
               IF (sysmatch.LE.maxSys(it)) THEN
!                Match found: fill the output matrix with this data
                 leng = sysA(it)%sys(sysmatch)%nPoints
                 DO ind = 1, leng
!                  Write output point data with bilinear interpolation
                   IF ( (sysA(it)%sys(sysmatch)%lon(ind).EQ.&
                         BL_lonprt).AND.&
                        (sysA(it)%sys(sysmatch)%lat(ind).EQ.&
                         BL_latprt) ) THEN
                      BL_hsprt(igrp) = sysA(it)%sys(sysmatch)%hs(ind)
                      BL_tpprt(igrp) = sysA(it)%sys(sysmatch)%tp(ind)
                      BL_dirprt(igrp) = sysA(it)%sys(sysmatch)%dir(ind)
                   ELSE IF ( (sysA(it)%sys(sysmatch)%lon(ind).EQ.&
                              BR_lonprt).AND.&
                             (sysA(it)%sys(sysmatch)%lat(ind).EQ.&
                              BR_latprt)) THEN
                      BR_hsprt(igrp) = sysA(it)%sys(sysmatch)%hs(ind)
                      BR_tpprt(igrp) = sysA(it)%sys(sysmatch)%tp(ind)
                      BR_dirprt(igrp) = sysA(it)%sys(sysmatch)%dir(ind)
                   ELSE IF ( (sysA(it)%sys(sysmatch)%lon(ind).EQ.&
                              TL_lonprt).AND.&
                             (sysA(it)%sys(sysmatch)%lat(ind).EQ.&
                              TL_latprt)) THEN
                      TL_hsprt(igrp) = sysA(it)%sys(sysmatch)%hs(ind)
                      TL_tpprt(igrp) = sysA(it)%sys(sysmatch)%tp(ind)
                      TL_dirprt(igrp) = sysA(it)%sys(sysmatch)%dir(ind)
                   ELSE IF ( (sysA(it)%sys(sysmatch)%lon(ind).EQ.&
                              TR_lonprt).AND.&
                             (sysA(it)%sys(sysmatch)%lat(ind).EQ.&
                              TR_latprt)) THEN
                      TR_hsprt(igrp) = sysA(it)%sys(sysmatch)%hs(ind)
                      TR_tpprt(igrp) = sysA(it)%sys(sysmatch)%tp(ind)
                      TR_dirprt(igrp) = sysA(it)%sys(sysmatch)%dir(ind)
                   END IF
                 END DO
!                  Compute output value using weights
!                  (only if output point is surrounded by valid points)
                   IF ( (BL_hsprt(igrp).NE.999.9999).AND. &
                        (BR_hsprt(igrp).NE.999.9999).AND. &
                        (TL_hsprt(igrp).NE.999.9999).AND. &
                        (TR_hsprt(igrp).NE.999.9999) ) THEN
                      hsprt(igrp) = BL_W * BL_hsprt(igrp) + &
                                    BR_W * BR_hsprt(igrp) + &
                                    TL_W * TL_hsprt(igrp) + &
                                    TR_W * TR_hsprt(igrp)
                      tpprt(igrp) = BL_W * BL_tpprt(igrp) + &
                                    BR_W * BR_tpprt(igrp) + &
                                    TL_W * TL_tpprt(igrp) + &
                                    TR_W * TR_tpprt(igrp)
                      BL_dirx = COS((270-BL_dirprt(igrp))*PI/180.)
                      BR_dirx = COS((270-BR_dirprt(igrp))*PI/180.)
                      TR_dirx = COS((270-TR_dirprt(igrp))*PI/180.)
                      TL_dirx = COS((270-TL_dirprt(igrp))*PI/180.)
                      BL_diry = SIN((270-BL_dirprt(igrp))*PI/180.)
                      BR_diry = SIN((270-BR_dirprt(igrp))*PI/180.)
                      TR_diry = SIN((270-TR_dirprt(igrp))*PI/180.)
                      TL_diry = SIN((270-TL_dirprt(igrp))*PI/180.)
                      dirprt(igrp)=270 - 180./PI* &
                         ATAN2(BL_W*BL_diry+BR_W*BR_diry+ &
                               TL_W*TL_diry+TR_W*TR_diry, &
                               BL_W*BL_dirx+BR_W*BR_dirx+ &
                               TL_W*TL_dirx+TR_W*TR_dirx)
                      IF (dirprt(igrp).GT.360.) THEN
                         dirprt(igrp) = dirprt(igrp) - 360.
                      END IF
                   ELSE
                      hsprt(igrp) = 999.9999
                      tpprt(igrp) = 999.9999
                      dirprt(igrp) = 999.9999
                   END IF
               END IF  
            END DO
      IF(outputType == 1) THEN
            WRITE(26,'(32F14.4)') lonprt,latprt, &
                       hsprt(1:10),tpprt(1:10),dirprt(1:10)  
        ENDIF
!/TRKNC      IF(outputType.EQ.3.OR.outputType.EQ.4) THEN
!/TRKNC        lonprt_nc(ip)=lonprt
!/TRKNC        latprt_nc(ip)=latprt
!/TRKNC        do igrp=1,10
!/TRKNC        hsprt_nc(igrp,ip,it)=hsprt(igrp)
!/TRKNC        tpprt_nc(igrp,ip,it)=tpprt(igrp)
!/TRKNC        dirprt_nc(igrp,ip,it)=dirprt(igrp)
!/TRKNC        enddo
!/TRKNC      ENDIF

         END DO
      END DO
!/TRKNC      IF(outputType.EQ.3.OR.outputType.EQ.4) THEN
!/TRKNC      call pt2netcdf(lonprt_nc,latprt_nc,hsprt_nc,tpprt_nc, &
!/TRKNC       dirprt_nc,noutp,date1,date2,dt,NTIME_NC,outputType)
!/TRKNC      ENDIF

      IF(outputType.EQ.1) CLOSE(26)

!-----Final SYSTEM output: point output (Nearest neighbor, as a double check)
      IF (testout) THEN
      OPEN(unit=28,file='sys_pnt_nn.ww3',status='unknown')
      WRITE(28,'(A)') '%'
      WRITE(28,'(A)') '%'
      WRITE(28,'(A)') '% WW3 Wave tracking point output'
      WRITE(28,'(A)') '%'
      WRITE(28,'(10A)') '%       Xp            Yp            ', &
           'HsSY01        HsSY02        HsSY03        HsSY04        ', &
           'HsSY05        HsSY06        HsSY07        HsSY08        ', &
           'HsSY09        HsSY10        ', &
           'TpSY01        TpSY02        TpSY03        TpSY04        ', &
           'TpSY05        TpSY06        TpSY07        TpSY08        ', &
           'TpSY09        TpSY10        ', &
           'DrSY01        DrSY02        DrSY03        DrSY04        ', &
           'DrSY05        DrSY06        DrSY07        DrSY08        ', &
           'DrSY09        DrSY10'   
      WRITE(28,'(10A)') '%       [degr]        [degr]        ', &
           '[m]           [m]           [m]           [m]           ', &
           '[m]           [m]           [m]           [m]           ', &
           '[m]           [m]           ', &
           '[sec]         [sec]         [sec]         [sec]         ', &
           '[sec]         [sec]         [sec]         [sec]         ', &
           '[sec]         [sec]         ', &
           '[degr]        [degr]        [degr]        [degr]        ', &
           '[degr]        [degr]        [degr]        [degr]        ', &
           '[degr]        [degr]'
      WRITE(28,'(A)') '%'

      DO it = 1,SIZE(sysA)
         WRITE(28,'(A,F15.6)') 'Time : ',wsdat(it)%date

         DO ip = 1,noutp
            hsprt(1:10) = 999.9999
            tpprt(1:10) = 999.9999
            dirprt(1:10) = 999.9999
            lonprt = 999.9999
            latprt = 999.9999

            DO j = 1, maxJ
               DO i = 1, maxI
!                 Write nearest nearbor output (no bilinear interpolation)
                  IF ( (lonout(ip).GE. &
                        (wsdat(1)%lon(i,j)-dlon/2)).AND. &
                        (lonout(ip).LT. &
                        (wsdat(1)%lon(i,j)+dlon/2)).AND. &
                        (latout(ip).GE. &
                        (wsdat(1)%lat(i,j)-dlat/2)).AND. &
                        (latout(ip).LT. &
                        (wsdat(1)%lat(i,j)+dlat/2)) ) &
                  THEN
                         lonprt = wsdat(1)%lon(i,j)
                         latprt = wsdat(1)%lat(i,j)
                  END IF
               END DO
            END DO
!           Loop through identified groups, limiting the output in file to 10
            DO igrp = 1,MIN(10,maxGroup)
!              Find system with this group tag
               sysmatch = 1
               DO WHILE (sysmatch.LE.maxSys(it))
                  IF (sysA(it)%sys(sysmatch)%grp.EQ.igrp) EXIT
                  sysmatch = sysmatch+1
               END DO
               IF (sysmatch.LE.maxSys(it)) THEN
!                Match found: fill the output matrix with this data
                 leng = sysA(it)%sys(sysmatch)%nPoints
                 DO ind = 1, leng
!                  Write nearest nearbor output (no bilinear interpolation)
                   IF ( (lonout(ip).GE. &
                        (sysA(it)%sys(sysmatch)%lon(ind)-dlon/2)).AND. &
                        (lonout(ip).LT. &
                        (sysA(it)%sys(sysmatch)%lon(ind)+dlon/2)).AND. &
                        (latout(ip).GE. &
                        (sysA(it)%sys(sysmatch)%lat(ind)-dlat/2)).AND. &
                        (latout(ip).LT. &
                        (sysA(it)%sys(sysmatch)%lat(ind)+dlat/2)) ) &
                      THEN
                         hsprt(igrp) = sysA(it)%sys(sysmatch)%hs(ind)
                         tpprt(igrp) = sysA(it)%sys(sysmatch)%tp(ind)
                         dirprt(igrp) = sysA(it)%sys(sysmatch)%dir(ind)
                   END IF
                 END DO
               END IF  
            END DO
            WRITE(28,'(32F14.4)') lonprt,latprt, &
                       hsprt(1:10),tpprt(1:10),dirprt(1:10)  
         END DO
      END DO

      CLOSE(28)
      END IF

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

      WRITE(20,*) 'In ww3_systrk: Deallocating wsdat ...'
      DO it=1,size(wsdat)
        IF (ASSOCIATED(wsdat(it)%lat)) DEALLOCATE(wsdat(it)%lat)
        IF (ASSOCIATED(wsdat(it)%lon)) DEALLOCATE(wsdat(it)%lon)
        IF (ASSOCIATED(wsdat(it)%par)) DEALLOCATE(wsdat(it)%par)
        IF (ASSOCIATED(wsdat(it)%wnd)) DEALLOCATE(wsdat(it)%wnd)
      END DO
      IF (ASSOCIATED(wsdat)) DEALLOCATE(wsdat)
      WRITE(20,*) '               Deallocating sysA ...'
      DO i=1,size(sysA)
       DO j=1,size(sysA(i)%sys)
         IF (ASSOCIATED(sysA(i)%sys(j)%i)) DEALLOCATE(sysA(i)%sys(j)%i)
         IF (ASSOCIATED(sysA(i)%sys(j)%j)) DEALLOCATE(sysA(i)%sys(j)%j)
         IF (ASSOCIATED(sysA(i)%sys(j)%lon)) &
                                         DEALLOCATE(sysA(i)%sys(j)%lon)
         IF (ASSOCIATED(sysA(i)%sys(j)%lat)) &
                                         DEALLOCATE(sysA(i)%sys(j)%lat)
         IF (ASSOCIATED(sysA(i)%sys(j)%hs))  &
                                          DEALLOCATE(sysA(i)%sys(j)%hs)
         IF (ASSOCIATED(sysA(i)%sys(j)%tp))  &
                                          DEALLOCATE(sysA(i)%sys(j)%tp)
         IF (ASSOCIATED(sysA(i)%sys(j)%dir)) &
                                         DEALLOCATE(sysA(i)%sys(j)%dir)
         IF (ASSOCIATED(sysA(i)%sys(j)%dspr)) &
                                        DEALLOCATE(sysA(i)%sys(j)%dspr)
       END DO
      END DO
      IF (ASSOCIATED(sysA)) DEALLOCATE(sysA)
      WRITE(20,*) '               Deallocating maxSys ...'
      IF (ASSOCIATED(maxSys)) DEALLOCATE(maxSys)
      CLOSE(20)

      WRITE(6,*) '... ww3_systrk completed successfully.'

      WRITE(6,999)

!/MPI      END IF   !/IF (rank.EQ.0)

!/MPI      CALL MPI_FINALIZE(IERR)
!/MPI/!     End of parallel region

!/F90  998 FORMAT ( ' ... finished. Elapsed time : ',F10.2,' s')
  993 FORMAT (/' *** WAVEWATCH III ERROR IN WW3_SYSTRK : '/           &
               '     OutputType=3 needs TRKNC switch ')
  994 FORMAT (/' *** WAVEWATCH III ERROR IN WW3_SYSTRK : '/           &
               '     OutputType=4 needs TRKNC and NC4 switch ')
  995 FORMAT (/' *** WAVEWATCH III ERROR IN WW3_SYSTRK : '/           &
               '     OutputType,',I3,'not valid. Options: 1,3,4')

  999 FORMAT (/15X,'End of program '/ &
               15X,'==============================================='/ &
               15X,'     *** WAVEWATCH III Wave system tracking ***  ')

 2000 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/    &
               '     ERROR IN OPENING INPUT FILE')
 2200 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/    &
               '     ERROR IN OPENING PARTITION FILE : ',A)

      END PROGRAM WW3_SYSTRK
!
!/TRKNC      subroutine t2netcdf(lons,lats,data_in,nlons,nlats,nsys,date1,date2,&
!/TRKNC         dt,ntime,ivar, outputType)
!/TRKNC      USE W3TIMEMD
!/TRKNC      use netcdf
!/TRKNC      implicit none
!/TRKNC      character (len = 15) :: file_name
!/TRKNC      integer, parameter :: ndims = 4
!/TRKNC      integer, parameter :: deflate = 1
!/TRKNC      integer :: outputType, ncid, oldMode
!/TRKNC      integer :: nlons,nlats,nsys,rec,ntime,ivar
!/TRKNC      double precision :: date1,date2,timenc
!/TRKNC      real :: data_in(nlons, nlats, nsys,ntime)
!/TRKNC      real :: lats(nlats), lons(nlons),dt
!/TRKNC      double precision :: times(ntime)
!/TRKNC      integer :: iyc,imc,idc,ihc,iminc,isc,Jday,Jday0
!/TRKNC      integer :: iret
!
!/TRKNC      integer :: lon_varid, lat_varid, rec_varid
!/TRKNC      character (len = *), parameter :: lsys_name = "system_index"
!/TRKNC      character (len = *), parameter :: lat_name = "latitude"
!/TRKNC      character (len = *), parameter :: lon_name = "longitude"
!/TRKNC      character (len = *), parameter :: time_name = "time"
!/TRKNC      integer :: sys_dimid, lon_dimid, lat_dimid, rec_dimid
!/TRKNC      integer :: start(ndims), count(ndims)
!
!/TRKNC      character (len = *), parameter :: var1_name="hs"
!/TRKNC      character (len = *), parameter :: var2_name="tp"
!/TRKNC      character (len = *), parameter :: var3_name="dir"
!/TRKNC      character (len = *), parameter :: var4_name="dspr"
!/TRKNC      integer :: var1_varid, var2_varid, var3_varid,var4_varid
!/TRKNC      integer :: dimids(ndims)
!
!/TRKNC      character (len = *), parameter :: units = "units"
!/TRKNC      character (len = *), parameter :: var1_units = "m"
!/TRKNC      character (len = *), parameter :: var2_units = "s"
!/TRKNC      character (len = *), parameter :: var3_units = "degrees"
!/TRKNC      character (len = *), parameter :: var4_units = "degrees"
!/TRKNC      character (len = *), parameter :: lat_units = "degrees_north"
!/TRKNC      character (len = *), parameter :: lon_units = "degrees_east"
!/TRKNC      iyc=date1/10000
!/TRKNC      imc=(date1-iyc*10000)/100
!/TRKNC      idc=int(date1-DBLE(iyc*10000)-DBLE(imc*100))
!/TRKNC      ihc=date2/10000
!/TRKNC      iminc=(date2-ihc*10000)/100
!/TRKNC      isc=date2-ihc*10000-100*iminc
!/TRKNC      timenc=DBLE(julday(idc,imc,iyc))+(DBLE(ihc)+(DBLE(iminc)+ &
!/TRKNC         (DBLE(isc)/60.0D0))/60.0D0)/24.0D0
!/TRKNC      Jday0=julday(1,1,1990)
!/TRKNC      timenc=timenc-Jday0
!/TRKNC      do rec=1,ntime
!/TRKNC       times(rec)=timenc+DBLE( (rec-1)*dt)/3600.0D0/24.0D0 
!/TRKNC      enddo
!/TRKNC      if( ivar == 1) then
!/TRKNC       file_name = "sys_hs.ww3.nc"
!/TRKNC      else if( ivar == 2) then
!/TRKNC       file_name = "sys_tp.ww3.nc"
!/TRKNC      else if( ivar == 3) then
!/TRKNC       file_name = "sys_dir.ww3.nc"
!/TRKNC      else
!/TRKNC      file_name = "sys_dspr.ww3.nc"
!/TRKNC      endif
!
!/TRKNC! create the netcdf file. 
!/TRKNC      if (outputType.EQ.3) then
!/TRKNC      call check( nf90_create(file_name, NF90_CLOBBER, ncid) )
!/TRKNC      endif
!/NC4        if(outputType.EQ.4) call check( nf90_create(file_name,NF90_NETCDF4,ncid))
!/TRKNC      call check ( nf90_set_fill(ncid,nf90_nofill,oldMode) )
!/TRKNC      call check( nf90_def_dim(ncid, lsys_name, nsys, sys_dimid) )
!/TRKNC      call check( nf90_def_dim(ncid, lat_name, nlats, lat_dimid) )
!/TRKNC      call check( nf90_def_dim(ncid, lon_name, nlons, lon_dimid) )
!/TRKNC      call check( nf90_def_dim(ncid, time_name, ntime, rec_dimid) )
!/TRKNC      call check( nf90_def_var(ncid, lat_name, NF90_REAL, lat_dimid,lat_varid))
!/NC4        call check( nf90_def_var_deflate(ncid,lat_varid,1,1,deflate) )
!/TRKNC      call check( nf90_def_var(ncid, lon_name, NF90_REAL, lon_dimid,lon_varid))
!/NC4        call check( nf90_def_var_deflate(ncid,lon_varid,1,1,deflate) )
!/TRKNC      call check( nf90_def_var(ncid,time_name,NF90_DOUBLE,rec_dimid,rec_varid))
!/NC4        call check( nf90_def_var_deflate(ncid,rec_varid,1,1,deflate) )
!
!/TRKNC      call check( nf90_put_att(ncid, lat_varid, units, lat_units) )
!/TRKNC      call check( nf90_put_att(ncid, lat_varid, 'long_name', 'latitude') )
!/TRKNC      call check( nf90_put_att(ncid, lat_varid, 'standard_name', 'latitude') )
!/TRKNC      call check( nf90_put_att(ncid, lat_varid, 'axis','Y'))
!/TRKNC      call check( nf90_put_att(ncid, lon_varid, units, lon_units) )
!/TRKNC      call check( nf90_put_att(ncid, lon_varid, 'long_name', 'longitude') )
!/TRKNC      call check( nf90_put_att(ncid, lon_varid, 'standard_name', 'longitude') )
!/TRKNC      call check( nf90_put_att(ncid, lon_varid, 'axis','X'))
!/TRKNC      call check(nf90_put_att(ncid,rec_varid,units,& 
!/TRKNC                  'days since 1990-01-01 00:00:00'))
!/TRKNC      call check(nf90_put_att(ncid,rec_varid,'long_name','julian day (UT)'))
!/TRKNC      call check( nf90_put_att(ncid, rec_varid,'standard_name', 'time') )
!/TRKNC      call check( nf90_put_att(ncid, rec_varid, 'conventions',&
!/TRKNC      'relative julian day with decimal part (as part of the day)' ) )
!/TRKNC      call check( nf90_put_att(ncid, rec_varid, 'axis','T'))
!
!/TRKNC      dimids = (/ lon_dimid, lat_dimid, sys_dimid, rec_dimid /)
!/TRKNC      if( ivar == 1) then
!/TRKNC      call check( nf90_def_var(ncid, var1_name, NF90_REAL, dimids,var1_varid) )
!/NC4        call check( nf90_def_var_deflate(ncid,var1_varid,1,1,deflate) )
!/TRKNC      call check( nf90_put_att(ncid, var1_varid, units, var1_units) )
!/TRKNC      call check( nf90_put_att(ncid, var1_varid,'long_name','significant_wave_height') )
!/TRKNC      call check( nf90_put_att(ncid, var1_varid,'missing_value','9999.00'))
!/TRKNC       else if( ivar == 2) then
!/TRKNC      call check( nf90_def_var(ncid, var2_name, NF90_REAL, dimids, var2_varid) )
!/NC4        call check( nf90_def_var_deflate(ncid,var2_varid,1,1,deflate) )
!/TRKNC      call check( nf90_put_att(ncid, var2_varid, units, var2_units) )
!/TRKNC      call check( nf90_put_att(ncid, var2_varid,'long_name','peak_period') )
!/TRKNC      call check( nf90_put_att(ncid, var2_varid,'missing_value','9999.00') )
!/TRKNC       else if ( ivar ==3 ) then
!/TRKNC      call check( nf90_def_var(ncid, var3_name, NF90_REAL, dimids, var3_varid) )
!/NC4        call check( nf90_def_var_deflate(ncid,var3_varid,1,1,deflate) )
!/TRKNC      call check( nf90_put_att(ncid, var3_varid, units, var3_units) )
!/TRKNC      call check( nf90_put_att(ncid, var3_varid,'long_name','peak_direction') )
!/TRKNC      call check( nf90_put_att(ncid, var3_varid,'missing_value','9999.00') )
!/TRKNC       else
!/TRKNC      call check( nf90_def_var(ncid, var4_name, NF90_REAL, dimids, var4_varid) )
!/NC4        call check( nf90_def_var_deflate(ncid,var4_varid,1,1,deflate) )
!/TRKNC      call check( nf90_put_att(ncid, var4_varid, units, var4_units) )
!/TRKNC      call check( nf90_put_att(ncid,var4_varid,'long_name','directional_spread') )
!/TRKNC      call check( nf90_put_att(ncid, var4_varid,'missing_value','9999.00') )
!/TRKNC      endif
!/TRKNC      call check( nf90_enddef(ncid) )
!  
!/TRKNC      call check( nf90_put_var(ncid, lat_varid, lats) )
!/TRKNC      call check( nf90_put_var(ncid, lon_varid, lons) )
!/TRKNC      call check( nf90_put_var(ncid, rec_varid, times) )
!  
!/TRKNC      count = (/ nlons, nlats, nsys, ntime /)
!/TRKNC      start = (/ 1, 1, 1, 1 /)
!/TRKNC      if( ivar == 1) then
!/TRKNC      call check( nf90_put_var(ncid, var1_varid, data_in, start = start, &
!/TRKNC                              count = count) )
!/TRKNC       else if( ivar == 2) then
!/TRKNC      call check( nf90_put_var(ncid, var2_varid, data_in, start = start, &
!/TRKNC                              count = count) )
!/TRKNC       else if( ivar == 3) then
!/TRKNC      call check( nf90_put_var(ncid, var3_varid, data_in, start = start, &
!/TRKNC                              count = count) )
!/TRKNC       else
!/TRKNC      call check( nf90_put_var(ncid, var4_varid, data_in, start = start, &
!/TRKNC                              count = count) )
!/TRKNC      endif
!/TRKNC      call check( nf90_close(ncid) )
!/TRKNC      end subroutine t2netcdf
!
!/TRKNC      subroutine check(status)
!/TRKNC      use netcdf
!/TRKNC      integer, intent ( in) :: status
!/TRKNC      if(status /= nf90_noerr) then 
!/TRKNC        write(6,996)
!/TRKNC          996 FORMAT (/' *** WAVEWATCH III ERROR IN WW3_SYSTRK:'/           &
!/TRKNC                       'netCDF error:')
!/TRKNC        print *, trim(nf90_strerror(status))
!/TRKNC        stop "Stopped in netcdf output part"
!/TRKNC      endif
!/TRKNC      end subroutine check  
!
!/TRKNC      subroutine pt2netcdf(longitude,latitude,hs,tp,&
!/TRKNC       dir,npoints,date1,date2,dt,ntime,outputType)
!/TRKNC      USE W3TIMEMD
!/TRKNC      use netcdf
!/TRKNC      implicit none
!/TRKNC      integer :: ntime,npoints,outputType
!/TRKNC      integer, parameter :: deflate = 1
!/TRKNC      integer :: iret, oldMode
!/TRKNC      integer :: ncid
!/TRKNC      integer :: system_index_dim
!/TRKNC      integer :: point_dim,rec_dim
!/TRKNC      integer :: nsys
!/TRKNC      integer :: start(3), count(3)
!/TRKNC      parameter (nsys = 10)
!/TRKNC      integer :: latitude_id
!/TRKNC      integer :: longitude_id
!/TRKNC      integer :: time_id
!/TRKNC      integer :: hs_id
!/TRKNC      integer :: tp_id
!/TRKNC      integer :: dir_id
!/TRKNC      integer :: time_rank
!/TRKNC      integer :: hs_rank
!/TRKNC      integer :: tp_rank
!/TRKNC      integer :: dir_rank
!/TRKNC      parameter (time_rank = 1)
!/TRKNC      parameter (hs_rank = 3)
!/TRKNC      parameter (tp_rank = 3)
!/TRKNC      parameter (dir_rank = 3)
!
!/TRKNC      integer :: hs_dims(hs_rank)
!/TRKNC      integer :: tp_dims(tp_rank)
!/TRKNC      integer :: dir_dims(dir_rank)
!/TRKNC      real :: latitude(npoints),dt
!/TRKNC      real :: longitude(npoints)
!/TRKNC      real :: hs(nsys, npoints, ntime)
!/TRKNC      real :: tp(nsys, npoints, ntime)
!/TRKNC      real :: dir(nsys, npoints, ntime)
!/TRKNC      integer :: iyc,imc,idc,ihc,iminc,isc,Jday,Jday0,rec
!/TRKNC      double precision date1,date2,timenc
!/TRKNC      double precision times(ntime)
!
!/TRKNC      iyc=date1/10000
!/TRKNC      imc=(date1-iyc*10000)/100
!/TRKNC      idc=int(date1-DBLE(iyc*10000)-DBLE(imc*100))
!/TRKNC      ihc=date2/10000
!/TRKNC      iminc=(date2-ihc*10000)/100
!/TRKNC      isc=date2-ihc*10000-100*iminc
!/TRKNC      timenc=DBLE(julday(idc,imc,iyc))+(DBLE(ihc)+(DBLE(iminc)+ &
!/TRKNC       (DBLE(isc)/60.0D0))/60.0D0)/24.0D0
!/TRKNC      Jday0=julday(1,1,1990)
!/TRKNC      timenc=timenc-Jday0
!/TRKNC      do rec=1,ntime
!/TRKNC       times(rec)=timenc+DBLE( (rec-1)*dt)/3600.0D0/24.0D0 
!/TRKNC      enddo
!
!/TRKNC      if(outputType.EQ.3) then
!/TRKNC      iret = nf90_create('sys_pnt.ww3.nc', NF90_CLOBBER, ncid)
!/TRKNC      endif
!/TRKNC      if (outputType.EQ.4) iret = nf90_create('sys_pnt.ww3.nc',NF90_NETCDF4, ncid)
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_set_fill(ncid,nf90_nofill,oldMode) 
!/TRKNC      call check(iret)
!/TRKNC! define dimensions
!/TRKNC      iret = nf90_def_dim(ncid, 'system_index', nsys, system_index_dim)
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_def_dim(ncid, 'point', npoints, point_dim)
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_def_dim(ncid, 'time', ntime, rec_dim)
!/TRKNC      call check(iret)
!/TRKNC! define variables
!/TRKNC      iret = nf90_def_var(ncid, 'latitude', NF90_REAL, point_dim, &
!/TRKNC          latitude_id)
!/TRKNC      call check(iret)
!/NC4        call check( nf90_def_var_deflate(ncid,latitude_id,1,1,deflate))
!/TRKNC      iret = nf90_def_var(ncid, 'longitude', NF90_REAL, point_dim, &
!/TRKNC        longitude_id)
!/TRKNC      call check(iret)
!/NC4        call check( nf90_def_var_deflate(ncid,longitude_id,1,1,deflate))
!/TRKNC      iret = nf90_def_var(ncid, 'time', NF90_DOUBLE, rec_dim, &
!/TRKNC         time_id)
!/TRKNC      call check(iret)
!/NC4        call check( nf90_def_var_deflate(ncid,time_id,1,1,deflate) )
!/TRKNC      hs_dims(3) = rec_dim
!/TRKNC      hs_dims(2) = point_dim
!/TRKNC      hs_dims(1) = system_index_dim
!/TRKNC      iret = nf90_def_var(ncid, 'hs', NF90_REAL,  &
!/TRKNC        hs_dims, hs_id)
!/TRKNC      call check(iret)
!/NC4        call check( nf90_def_var_deflate(ncid,hs_id,1,1,deflate))
!/TRKNC      tp_dims(3) = rec_dim
!/TRKNC      tp_dims(2) = point_dim
!/TRKNC      tp_dims(1) = system_index_dim
!/TRKNC      iret = nf90_def_var(ncid, 'tp', NF90_REAL, &
!/TRKNC        tp_dims, tp_id)
!/TRKNC      call check(iret)
!/NC4        call check( nf90_def_var_deflate(ncid,tp_id,1,1,deflate))
!/TRKNC      dir_dims(3) = rec_dim
!/TRKNC      dir_dims(2) = point_dim
!/TRKNC      dir_dims(1) = system_index_dim
!/TRKNC      iret = nf90_def_var(ncid, 'dir', NF90_REAL, & 
!/TRKNC        dir_dims, dir_id)
!/TRKNC      call check(iret)
!/NC4        call check( nf90_def_var_deflate(ncid,dir_id,1,1,deflate))
!/TRKNC! assign attributes
!/TRKNC      iret = nf90_put_att(ncid, latitude_id, 'units', 'degrees_north')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, latitude_id, 'long_name', 'latitude')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, latitude_id, 'standard_name', 'latitude')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, latitude_id, 'axis', 'Y')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, longitude_id, 'units', 'degrees_east')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, longitude_id,'long_name','longitude')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, longitude_id,'standard_name','longitude')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, longitude_id, 'axis', 'X')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, time_id, 'units', &
!/TRKNC        'days since 1990-01-01 00:00:00')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, time_id, 'long_name','julian day(UT)')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, time_id, 'standard_name','time')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, time_id, 'conventions', &
!/TRKNC       'relative julian day with decimal part (as part of the day)')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, time_id, 'axis', 'T')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, hs_id, 'units', 'm')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, hs_id,'long_name','significant_wave_height')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, hs_id, 'missing_value', & 
!/TRKNC        '999.9999')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, tp_id, 'units', 's')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, tp_id,'long_name','peak_period')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, tp_id, 'missing_value', &
!/TRKNC       '999.9999')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, dir_id, 'units', 'degrees')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, dir_id,'long_name','peak_direction')
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_att(ncid, dir_id, 'missing_value',& 
!/TRKNC         '999.9999')
!/TRKNC      call check(iret)
!/TRKNC! leave define mode
!/TRKNC      iret = nf90_enddef(ncid)
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_var(ncid, latitude_id, latitude)
!/TRKNC      call check(iret)
!
!/TRKNC      iret = nf90_put_var(ncid, longitude_id, longitude)
!/TRKNC      call check(iret)
! 
!/TRKNC      iret = nf90_put_var(ncid, time_id, times)
!/TRKNC      call check(iret)
! 
!/TRKNC      start = (/ 1, 1, 1 /)
!/TRKNC      count = (/ nsys,npoints,ntime /)
!
!/TRKNC      iret = nf90_put_var(ncid, hs_id, hs,&
!/TRKNC       start = start, count = count )
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_var(ncid, tp_id, tp, &
!/TRKNC       start = start, count = count )
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_put_var(ncid, dir_id, dir,&
!/TRKNC       start = start, count = count )
!/TRKNC      call check(iret)
!/TRKNC      iret = nf90_close(ncid)
!/TRKNC      call check(iret)
!/TRKNC      return
!/TRKNC      end subroutine pt2netcdf

