WW3/w3strkmd_8F90_source.html

#include "w3macros.h"

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

MODULE w3strkmd

  !/

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

  !/                  | WAVEWATCH III           NOAA/NCEP |

  !/                  |     A. J. van der Westhuysen      |

  !/                  |            Jeff Hanson            |

  !/                  |        Eve-Marie Devaliere        |

  !/                  |                        FORTRAN 95 |

  !/                  | Last update :         03-Mar-2016 |

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

  !/

  !/    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 )

  !/    26-Sep-2016 : Optimization updates (A. van der Westhuysen)

  !/                                                      ( version 5.15 )

  !/    03-Mar-2016 : Optimization updates for INTERSECT,

  !/                  UNION, UNIQUE, SORT, SETDIFF, FINDIJ

  !/                  (S. Zieger, BoM Australia)          ( 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 constants, only: file_endian

  IMPLICIT NONE

  !

  !  1. Purpose :

  !

  !     Module containing data structures and subroutines for spatial and

  !     temporal tracking (part of wave partitioning).

  !

  !  2. Method :

  !

  !     Read raw partitioning data.

  !     Perform tracking in space.

  !     Perform tracking in time.

  !

  !  3. Variables and types :

  !

  !     NOTE: In Fortran 90/95 derived types cannot contain allocatable arrays.

  !           The same functionality is achieved here using pointers (pointing

  !           to unnamed allocatable arrays). Can be replaced by allocatable arrays

  !           when transitioning to the Fortran 2003 standard.

  !

  !     Name    Type       Description

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

  !     param   Der. type  structure of basic spectrally partitioned results at a geo point

  !     hs      Real arr   array of sign. wave height partitions

  !     tp      Real arr   array of peak period partitions

  !     dir     Real arr   array of mean direction partitions

  !     dspr    Real arr   array of mean directional spread (one-sided) of partitions

  !     wf      Real arr   array of wind fraction

  !     ipart   Int arr    array of partition indices

  !     sys     Int arr    array of system indices to which a given partition has been assigned

  !     ngbrSys Int arr    array of system indices of neighboring grid points

  !     checked    Int      0 = geo point not checked yet (in SUBROUTINE findSys)

  !                         1 = geo point has been checked

  !                        -1 = geo land point (i.e. no partitioning data found for this point)

  !                        -2 = geo land point, second passing.

  !

  TYPE param

    REAL    :: hs(10)

    REAL    :: tp(10)

    REAL    :: dir(10)

    REAL    :: dspr(10)

    !         REAL    :: wf(10)

    INTEGER :: ipart(10)

    INTEGER :: sys(10)

    INTEGER :: ngbrsys(50)

    INTEGER :: checked

  END TYPE param

  !

  !     wind    Der. type     structure containing wind-related parameters

  !     wdir       Real       wind direction at grid point (Nautical or Cartesian, invariant)

  !     wspd       Real       wind speed at grid point

  !

  TYPE wind

    REAL    :: wdir

    REAL    :: wspd

  END TYPE wind

  !

  !     dat2d   Der. type     2d data structure for storing raw partitioned data

  !     lat     Real arr      2d array of latitudes of input partitioned data

  !     lon     Real arr      2d array of longitudes of input partitioned data

  !     par  type(param) arr  2d array of partitioned parameter structures

  !     wnd  type(wind) arr   2d array of wind parameter structures

  !     maxi    Int           size of 2d array of raw partitioned data in i dimension

  !     maxj    Int           size of 2d array of raw partitioned data in j dimension

  !

  TYPE dat2d

    real*8               :: date

    REAL, POINTER        :: lat(:,:)

    REAL, POINTER        :: lon(:,:)

    TYPE(param), POINTER :: par(:,:)

    TYPE(wind), POINTER  :: wnd(:,:)

    INTEGER              :: maxi

    INTEGER              :: maxj

  END TYPE dat2d

  !

  !     neighbr  Der. type    structure for storing data of neighboring grid point

  !     par     type(param)   partitioned parameter structure at neighboring grid point

  !     i          Int        i index of neighboring grid point

  !     j          Int        j index of neighboring grid point

  !

  TYPE neighbr

    TYPE(param) :: par

    INTEGER     :: i

    INTEGER     :: j

  END TYPE neighbr

  !

  !     mtchsys  Der. type    structure for storing data of matched systems

  !     sysVal   Int arr      array of indices of matched systems

  !     tpVal    Real arr     array of peak period values of matched systems

  !     wfVal    Real arr     array of wind fraction values of matched systems

  !

  TYPE mtchsys

    INTEGER     :: sysval(50)

    REAL        :: tpval(50)

    REAL        :: dirval(50)

    REAL        :: hsval(50)

    !         REAL        :: wfVal(50)

  END TYPE mtchsys

  !

  !     system   Der. type    structure for storing spatially tracked systems (one time level)

  !     hs       Real arr     sign wave height field assoc with wave system (in 1d array)

  !     tp       Real arr     peak period field assoc with wave system (in 1d array)

  !     dir      Real arr     mean direction field assoc with wave system (in 1d array)

  !     dspr     Real arr     mean directional spread field assoc with wave system (in 1d array)

  !     wf       Real arr     wind fraction assoc with wave system (in 1d array)

  !     i        Int arr      i index of geo grid point in wave system (in 1d array)

  !     j        Int arr      j index of geo grid point in wave system (in 1d array)

  !     lat      Real arr     latitudes of grid point in wave system (in 1d array)

  !     lon      Real arr     longitudes of grid point in wave system (in 1d array)

  !     sysInd   Int          index of current wave system

  !     hsMean   Real         spatial mean sign wave height of current wave system

  !     tpMean   Real         spatial mean peak period of current wave system

  !     dirMean  Real         spatial mean wave direction of current wave system

  !     wfMean   Real         spatial mean wind fraction of current wave system

  !     nPoints  Int          total number of grid points in current wave system

  !     ngbr     Int arr      indices of neighboring wave systems

  !     grp      Int          time-tracked group that system is assigned to

  !

  TYPE system

    REAL, POINTER    :: hs(:)

    REAL, POINTER    :: tp(:)

    REAL, POINTER    :: dir(:)

    REAL, POINTER    :: dspr(:)

    !         REAL, POINTER    :: wf(:)

    INTEGER, POINTER :: i(:)

    INTEGER, POINTER :: j(:)

    INTEGER, POINTER :: indx(:)

    REAL, POINTER    :: lat(:)

    REAL, POINTER    :: lon(:)

    INTEGER          :: sysind

    REAL             :: hsmean

    REAL             :: tpmean

    REAL             :: dirmean

    !         REAL             :: wfMean

    INTEGER          :: npoints

    INTEGER          :: ngbr(1000)

    INTEGER          :: grp

  END TYPE system

  !

  !     timsys   Der. type     structure for storing time-tracked systems (all time levels)

  !     sys  type(system) arr  array of all spatially+temporally tracked systems at given

  !                            time level

  !

  TYPE timsys

    TYPE(system), POINTER :: sys(:)

  END TYPE timsys

  !

  !     sysmemory  Der. type   Structure to store key characteristics of systems over multiple

  !                            time levels. Used during the time tracking routine.

  !

  TYPE sysmemory

    INTEGER :: grp

    INTEGER :: npoints

    INTEGER, POINTER :: indx(:)

    INTEGER :: updated

    INTEGER :: length

    REAL    :: lonmean

    REAL    :: latmean

    REAL    :: tpmean

    REAL    :: dirmean

  END TYPE sysmemory

  !

  !  4. Subroutines and functions used :

  !

  !     a. Main subroutines for spatial/temporal tracking:

  !

  !     waveTracking_NWS_V2  main subroutine of spatial and temporal tracking algorithm

  !     spiralTrackV3        performs the spatial spiral tracking for a given time step

  !     timeTrackingV2       performs the time tracking of all wave systems

  !

  !     b. Auxiliary subroutines and functions for tracking:

  !

  !     findWay              find direction and no. steps in spatial search spiral

  !     findNext             find next point on spatial search spiral

  !     findSys              find all neighboring wave systems for given grid point

  !     combineWaveSystems   combine wave systems, then remove small and low-energy systems

  !     printFinalSys        output the final output systems for this time step

  !     combineSys           combine wave systems

  !     combinePartitionsV2  combine two partitions that have been assigned to the same system

  !     func. mean_angleV2   compute the mean direction from array of directions

  !     findIJV4             Find indices of system "a" that lie over or next to system "b"

  !

  !     c. Simple data manipulation (based on Matlab intrinsic functions):

  !

  !     UNIQUE               removes duplicate reals from an vector

  !     SORT                 sorts the vector in ascending or descending order

  !     SETDIFF              returns elements in vector1 that are not in vector2

  !     INTERSECT            returns elements that are mutual in vector1 and vector2

  !     UNION                returns the union of vector1 and vector2

  !     func. LENGTH         finds no. of indices in vector not filled with blank entries.

  !     func. FINDFIRST      returns index of first instance of a search value in vector

  !     func. STD            computes standard deviation

  !

  !  5. Called by :

  !

  !     WW3_SYSTRK (main program)

  !

  !  6. Error messages :

  !

  !  7. Remarks :

  !

  !  8. Structure :

  !

  !     The structure of the tracking algorithm is the following

  !     (parentheses indicate minor subroutines and functions):

  !

  !     +---- SUBROUTINE waveTracking_NWS_V2       main subroutine of spatial and temporal tracking algorithm

  !     |        (CALL UNIQUE)                     removes duplicate reals from an vector

  !     |        CALL spiralTrackV3                See below

  !     |        CALL timeTrackingV2               See below

  !     |

  !     +---+--- SUBROUTINE spiralTrackV3          performs the spatial spiral tracking for a given time step

  !     |   |       CALL findWay                   find direction and no. steps in spatial search spiral

  !     |   |       CALL findNext                  find next point on spatial search spiral

  !     |   |       CALL findSys                   See below

  !     |   |       CALL combineWaveSystems        See below

  !     |   |

  !     |   +------ SUBROUTINE findSys             find all neighboring wave systems for given grid point

  !     |   |          (CALL UNIQUE)

  !     |   |          CALL combinePartitionsV2    combine two partitions that have been assigned to the same system

  !     |   |

  !     |   +---+-- SUBROUTINE combineWaveSystems  combine wave systems, then remove small and low-energy systems

  !     |       |      CALL printFinalSys          See below

  !     |       |      CALL combineSys             See below

  !     |       |

  !     |       +----- SUBROUTINE printFinalSys    output the final output systems for this time step

  !     |       |         (CALL UNIQUE)

  !     |       |         (CALL SETDIFF)           returns elements in vector1 that are not in vector2

  !     |       |         (CALL SORT)              sorts the vector in ascending or descending order

  !     |       |

  !     |       +----- SUBROUTINE combineSys       combine wave systems

  !     |                 (CALL SORT)

  !     |                 (CALL UNIQUE)

  !     |                 (CALL UNION)             returns the union of vector1 and vector2

  !     |                 (CALL SETDIFF)

  !     |                 CALL findIJV4            Find indices of system "a" that lie over or next to system "b"

  !     |                 CALL combinePartitionsV2

  !     |

  !     +------- SUBROUTINE timeTrackingV2         performs the time tracking of all wave systems

  !                 (CALL SORT)

  !                 (CALL SETDIFF)

  !

  !  9. Switches :

  !

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

  !       !/MPI   Id.

  !

  ! 10. Source code :

  !

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

  !/

CONTAINS

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

  SUBROUTINE 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   )

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          4-Jan-2013 |

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

    !/

    !/    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 )

    !/

    !/    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.

    !/

    IMPLICIT NONE

#ifdef W3_MPI


    include "mpif.h"

#endif

    !

    !  1. Purpose :

    !

    !     Main subroutine of spatial and temporal tracking algorithm

    !

    !  2. Method

    !

    !     (1) Read the raw partitioning output from one of two file formats:

    !         (a) "partRes" format of IFP-SWAN (intype=1), or

    !         (b) WW3 spectral bulletin format (intype=2).

    !         If intype=0, the partition data is read from memory (not activated yet).

    !

    !     (2) Perform tracking in space by calling subroutine spiralTrackV3

    !     (3) Perform tracking in time by calling subroutine timeTrackingV2

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     intype         Int       input  For coupling: Type of input (0 = memory; 1 = partRes file; 2 = WW3 part file)

    !     tmax           Int       input  For coupling: Value of maxTs to apply (1 or 2)

    !     tcur           Int       input  For coupling: Index of current time step (1 or 2)

    !     filename       Char      input  File name of locally partitioned data output

    !     tstart         Char      input  Start time in raw partition file (if used)

    !     tend           Char      input  End time in raw partition file (if used)

    !     minlon         Real      input  Lower lon boundary of domain to be processed

    !     maxlon         Real      input  Upper lon boundary of domain to be processed

    !     minlat         Real      input  Lower lat boundary of domain to be processed

    !     maxlat         Real      input  Upper lat boundary of domain to be processed

    !     dirKnob        Real      input  Parameter in direction for combining fields in space

    !     perKnob        Real      input  Parameter in period for combining fields in space

    !     hsKnob         Real      input  Parameter in wave height for purging fields

    !     wetPts         Real      input  Percentage of wet points for purging fields (fraction)

    !     seedLat        Real      input  Start Lat for tracking spiral (if =0 centre of field is used)

    !     seedLon        Real      input  Start Lon for tracking spiral (if =0 centre of field is used)

    !     dirTimeKnob    Real      input  Parameter in direction for combining fields in time

    !     tpTimeKnob     Real      input  Parameter in period for combining fields in time

    !     paramFile      Char      input  File name of partitioning parameters                           Is this used???

    !     sys         TYPE(timsys) output Final set of spatially and temporally tracked systems

    !     wsdat       TYPE(dat2d)  output Final version of 2D (gridded) partition data

    !     maxGroup       Int       output Maximum number of wave systems ("groups") tracked in time

    !

    CHARACTER    :: filename*50, paramFile*32

    REAL         :: dirKnob, perKnob, hsKnob, wetPts, seedLat, &

         seedLon, dirTimeKnob, tpTimeKnob

    real*8       :: tstart, tend

    INTEGER      :: maxGroup, intype, tmax, tcur, ntint

    INTEGER, POINTER :: maxSys(:)

    TYPE(dat2d), POINTER :: wsdat(:)

    TYPE(timsys), POINTER :: sysA(:), sysAA(:)

    INTEGER      :: NumConsSys, iConsSys

    REAL         :: dt

    REAL         :: minlon, maxlon, minlat, maxlat

    INTEGER      :: mxcwt, mycwt


    !     Note: Variables wsdat, sysA and maxSys have IN/OUT intent so that they

    !     can be manipulated outside of this subroutine, e.g. re-indexing of

    !     systems and groups during the simulation.

    INTENT (IN)  intype, tmax, tcur, filename, paramfile, &

         minlon, maxlon, minlat, maxlat, &

         hsknob, wetpts, seedlat, seedlon, &

         dirknob, perknob, dirtimeknob, tptimeknob

    INTENT (OUT) maxgroup

    !      INTENT (IN OUT) wsdat, sysA, maxSys

    !

    !     Local variables

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

    !     llat     Real    Latitude of partition point, from input file

    !     llon     Real    Longitude of partition point, from input file

    !     ts       Real    Time step of partition, from input file

    !     hs0      Real    Wave height of partition, from input file

    !     tp0      Real    Peak period of partition, from input file

    !     dir0     Real    Mean direction of partition, from input file

    !     dspr0    Real    Mean directional spread of partition, from input file

    !     wf0      Real    wind fraction of partition, from input file (removed)

    !     wndSpd0  Real    Wind speed of partition, from input file

    !     wndDir0  Real    Wind direction of partition, from input file

    !     wndFce0  Real    Wind force of partition, from input file (not, used; removed)

    !     tss      Int.    Time step counter

    !     t0       Int     Index of first time step to compute

    !

    LOGICAL                :: file_exists, FLFORM, LOOP

    LOGICAL                :: testout

    parameter(testout = .false.)

    CHARACTER              :: dummy*10, dummyc*12

    CHARACTER(LEN=10)      :: VERPRT

    CHARACTER(LEN=35)      :: IDSTR

    CHARACTER(LEN=78)      :: headln1

    CHARACTER(LEN=51)      :: headln2

    INTEGER                :: line

    INTEGER, ALLOCATABLE   :: ts(:), tmp_i4(:)

    REAL, ALLOCATABLE      :: llat(:),llon(:),hs0(:), &

         tp0(:),dir0(:),dspr0(:),&

         wndSpd0(:),wndDir0(:)

    real*8, ALLOCATABLE :: date0(:),tmp_r8(:)

    INTEGER    :: maxTs, t0, nout1, nout2, maxI, maxJ

    REAL, ALLOCATABLE :: mlon(:,:), mlat(:,:), tmp_r4(:)

    REAL, POINTER :: uniqueTim(:),uniqueLatraw(:),uniqueLonraw(:), &

         uniqueLat(:),uniqueLon(:)

    INTEGER    :: ioerr,ierr, i, j, k, l, alreadyIn, ok, tss, tsA

    INTEGER    :: maxPart, DATETIME(2)

    INTEGER    :: tstep, iline, numpart, skipln, readln, filesize

    REAL       :: x,y,wnd,wnddir

    REAL       :: invar1, invar2, invar3, invar4

    REAL       :: invar5, invar6, invar7

    REAL, ALLOCATABLE :: phs(:),ptp(:),pdir(:),pspr(:),pwf(:) ! current partition values

    real*8     :: date1, date2, ttest, ttemp

    INTEGER    :: ic, leng, maxpartout                                  ! Remove?

    REAL       :: dx

    INTEGER    :: latind1, latind2, lonind1, lonind2

    REAL       :: lonext, latext


#ifdef W3_MPI

    INTEGER    :: rank, irank, nproc, EXTENT, DOMSIZE, tag1, tag2

    !      INTEGER    :: MPI_INT_DOMARR, MPI_REAL_DOMARR

    INTEGER    :: MPI_STATUS(MPI_STATUS_SIZE)

    INTEGER    :: REQ(16)

    !    INTEGER    :: ISTAT(MPI_STATUS_SIZE,16)

    REAL       :: COMMARR1(44)

    INTEGER    :: COMMARR2(11)

#endif

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      UNIQUE

    !      spiralTrackV3

    !      timeTrackingV2

    !

    !  5. Subroutines calling

    !

    !      WW3_SYSTRK (main program)

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See above

    !

    !  9. Switches :

    !

    !     None defined yet.

    !

    ! 10. Source code :

    !

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


#ifdef W3_MPI

    CALL mpi_comm_rank(mpi_comm_world, rank, ierr)

    CALL mpi_comm_size(mpi_comm_world, nproc, ierr)

#endif

    NULLIFY( sysa )

    NULLIFY( maxsys )


    !     Select input type for raw partitioning data

    IF ((intype.EQ.1).OR.(intype.EQ.2)) THEN

      !        Raw partitioning data is coming from an input file.

      !        Read file here, and set up 2d array wsdat with the data.

      t0 = 1

      IF (intype.EQ.1) THEN

#ifdef W3_MPI

        IF (rank.EQ.0) THEN

#endif

          !           Read partRes format file

          WRITE(20,*) 'Reading partRes partitioning file...'

          filesize = 7500000

          ALLOCATE(ts(filesize))

          ALLOCATE(llat(filesize))

          ALLOCATE(llon(filesize))

          ALLOCATE(hs0(filesize))

          ALLOCATE(tp0(filesize))

          ALLOCATE(dir0(filesize))

          ALLOCATE(dspr0(filesize))

          !            ALLOCATE(wf0(filesize))

          ALLOCATE(wndspd0(filesize))

          ALLOCATE(wnddir0(filesize))

          ALLOCATE(date0(filesize))

          WRITE(20,*) '*** Max number of lines read from "partRes" ', &

               'input file is = ',filesize,'!'

          WRITE(6,*) 'Reading partRes file...'

          INQUIRE(file=filename, exist=file_exists)

          IF (.NOT.file_exists) THEN

            WRITE(20,2001)

            WRITE(6,2001)

            stop 1

          END IF

          OPEN(unit=11,file=filename,status='old')

          line = 1

          DO WHILE (.true.)

            READ (11, *, END=113) dummyc,llat(line),llon(line),   &

                 ts(line),hs0(line),tp0(line),dir0(line), &

                 wndspd0(line),wnddir0(line),invar7

            !partRes file does not contain the dspr variable

            dspr0(line) = 9999.

            !               wf0(line) = 9999.

            line = line+1

          ENDDO

113       ierr = -1

          CLOSE(11)

          line = line-1

          WRITE(6,*) '... finished'

          !            DEALLOCATE(date0)

#ifdef W3_MPI

        END IF

#endif

      ELSE IF (intype.EQ.2) THEN

#ifdef W3_MPI

        IF (rank.EQ.0) THEN

#endif

          !           Read WW3 Spectral Partition format file

          !           Query input file to determine required array sizes

          INQUIRE(file=filename, exist=file_exists)

          IF (.NOT.file_exists) THEN

            WRITE(20,2001)

            WRITE(6,2001)

            stop 1

          END IF

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

          !/          Test unformatted read

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

          OPEN(unit=11,file=filename,form='UNFORMATTED', convert=file_endian,status='OLD',access='STREAM')

          READ(11,err=802,iostat=ioerr) i

          CLOSE(11)

          !/          --- First four-byte integer could possibly be byte-swapped,

          !               if ww3_shel was compiled on a different architecture. ---

          k = swapi4(i)

          flform = .NOT.(i.EQ.(len(idstr)+len(verprt)).OR.&

               k.EQ.(len(idstr)+len(verprt))    )

          !        ======== COUNT LOOP ===========

          IF (flform) THEN

            !              Input file in formatted ASCII

            WRITE(6,*) 'Reading formatted ASCII file...'

            OPEN(unit=11,file=filename,status='old')

            READ(11,'(78A)') headln1

            idstr = headln1(1:len(idstr))

            READ(11,'(78A)') headln1

            READ(11,'(51A)') headln2

          ELSE

            IF (k.EQ.(len(idstr)+len(verprt))) THEN

              !/           --- Stop here. The file appears to be endian encoded

              !                different from the native machine format. And, the

              !                compiler option will override support for FORTRAN

              !                convert statements convert=file_endian or

              !                convert='little_endian'. ---

              WRITE(20,1200)

              WRITE(6,1200)

              stop 1

            ELSE

              !                Input file in unformatted binary

              WRITE(6,*) 'Reading binary formatted file...'

              OPEN(unit=11,file=filename,form='UNFORMATTED', convert=file_endian, &

                   status='OLD')

            ENDIF

            rewind(11)

            READ(11,err=802,iostat=ioerr) idstr,verprt

            READ(11,err=802,iostat=ioerr) headln1

            READ(11,err=802,iostat=ioerr) headln2

          END IF

          !/

          IF (idstr(1:9).ne.'WAVEWATCH') THEN

            CLOSE(11)

            WRITE(20,1300)

            WRITE(6,1300)

            stop 1

          ENDIF

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

          !/          Skip to start time

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

          skipln = 3

          ttest = 0

          DO WHILE (ttest.LT.tstart)

            IF (flform) THEN

              READ (11,1000,err=802,END=112) date1,date2,x,y, &

                   numpart,wnd,wnddir,invar6,invar7

#ifdef W3_del

              write(*,*) '0:',x,y,numpart

#endif

              skipln = skipln+1

            ELSE

              READ (11,err=802,iostat=ioerr) datetime,x,y, &

                   dummy,numpart,invar1,wnd,wnddir, &

                   invar5,invar6

              ! write(*,*) '0:',DATETIME,numpart

              date1=dble(datetime(1))

              date2=dble(datetime(2))

            END IF

            ttest = date1 + date2*1.0e-6

            IF (flform) THEN

              DO line = 1,numpart+1

                READ(11,1010,END=111,ERR=802,IOSTAT=IOERR) &

                     invar1,invar2,invar3,invar4

                ! write(*,*) '0+:',line,numpart+1,invar1,invar2,invar3,invar4

                skipln = skipln+1

              END DO

            ELSE

              DO line = 1,numpart+1

                READ (11,err=802,iostat=ioerr) iline,invar1, &

                     invar2,invar3,invar4,invar5,invar6

                ! write(*,*) '0+:',line,iline,invar1,invar2,invar3,invar4,invar5,invar6

              END DO

            END IF

          END DO

          skipln = skipln-numpart-1-1

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

          !           Read file for ntint time levels

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

          readln = numpart

          tstep = 1

          ttemp = tstart

          maxpart = numpart

          DO WHILE (tstep.LE.ntint)

            IF (readln.GT.0) THEN

              IF (flform) THEN

                READ (11,1000,err=802,END=111) date1,date2,x,y, &

                     numpart,wnd,wnddir,invar6,invar7

              ELSE

                READ (11,END=111,ERR=802,IOSTAT=IOERR) DATETIME,  &

                     x,y,dummy,numpart,wnd,wnddir,invar5,invar6,invar7

                ! write(*,*) '1:',numpart,x,y

                date1=dble(datetime(1))

                date2=dble(datetime(2))

              END IF

              maxpart = max(maxpart,numpart)

            END IF


            ttest = date1 + date2*1.e-6

            IF (ttest.GT.ttemp) THEN

              tstep = tstep+1

              ttemp = ttest

              IF (tstep.GT.ntint) EXIT

            END IF

            IF (flform) THEN

              DO line = 1,numpart+1

                READ (11,1010,END=111,ERR=802,IOSTAT=IOERR)      &

                     invar1,invar2,invar3,invar4

#ifdef W3_del

                write(*,'(A,2I6,4F7.2)') '1+:',line,numpart+1,invar1,invar2,invar3,invar4

#endif

                readln = readln+1

              END DO

            ELSE

              DO line = 1,numpart+1

                READ (11,END=111,ERR=802,IOSTAT=IOERR) iline,invar1,&

                     invar2,invar3,invar4,invar5,invar6

                readln = readln+1

              END DO

            END IF

          ENDDO

111       CONTINUE

          CLOSE(11)

          !        ===== END COUNT LOOP =====

          !        ===== START READ LOOP =====

          ALLOCATE(ts(readln))

          ALLOCATE(llat(readln))

          ALLOCATE(llon(readln))

          ALLOCATE(hs0(readln))

          ALLOCATE(tp0(readln))

          ALLOCATE(dir0(readln))

          ALLOCATE(dspr0(readln))

          !            ALLOCATE(wf0(readln))

          ALLOCATE(wndspd0(readln))

          ALLOCATE(wnddir0(readln))

          ALLOCATE(date0(readln))

          ts(1:readln)   = -1

          llat(1:readln) = 9999.

          llon(1:readln) = 9999.

          hs0(1:readln)  = 9999.

          tp0(1:readln)  = 9999.

          dir0(1:readln)  = 9999.

          dspr0(1:readln)  = 9999.


          IF (flform) THEN

            OPEN(unit=11,file=filename,status='old')

          ELSE

            OPEN(unit=11,file=filename,status='old', &

                 form='unformatted', convert=file_endian)

          END IF

          line = 1

          tstep = 1

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

          !/          Skip to start time

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

          IF (flform) THEN

            DO i = 1,skipln

              READ (11, *)

            END DO

          ELSE

            ! --- Repeat from above since access='DIRECT'

            !     does not support fseek and ftell. ---

            READ(11,END=112,ERR=802,IOSTAT=IOERR) IDSTR,verprt

            READ(11,END=112,ERR=802,IOSTAT=IOERR) headln1

            READ(11,END=112,ERR=802,IOSTAT=IOERR) headln2

            !/ --- allocate buffer for all partition parameters

            !/     for a single grid point  ---

            IF (.NOT.ALLOCATED(phs)) ALLOCATE(phs(maxpart))

            IF (.NOT.ALLOCATED(ptp)) ALLOCATE(ptp(maxpart))

            IF (.NOT.ALLOCATED(pdir)) ALLOCATE(pdir(maxpart))

            IF (.NOT.ALLOCATED(pspr)) ALLOCATE(pspr(maxpart))

            IF (.NOT.ALLOCATED(pwf)) ALLOCATE(pwf(maxpart))


            ttest = 0


            DO WHILE (ttest.LT.tstart)

              READ (11,END=112,ERR=802,IOSTAT=IOERR) DATETIME, &

                   invar1,invar2,dummy,numpart,invar3,   &

                   invar4,invar5,invar6,invar7

              date1=dble(datetime(1))

              date2=dble(datetime(2))

              ttest = date1 + date2*1.0e-6

              !/ --- reset buffer ---

              phs(:) = 0.

              ptp(:) = 0.

              pdir(:) = 0.

              pspr(:) = 0.

              pwf(:) = 0.


              !/ --- fill buffer with partition data ---

              READ (11,END=112,ERR=802,IOSTAT=IOERR) iline,invar1, &

                   invar2,invar3,invar4,invar5,invar6

              DO i = 1,numpart

                READ (11,END=112,ERR=802,IOSTAT=IOERR) iline,      &

                     phs(i),ptp(i),invar3,pdir(i),pspr(i),pwf(i)

              END DO

            END DO

            !/ --- move buffer content to data array ---

            DO i=1,numpart

              hs0(line)   = phs(i)

              tp0(line)   = ptp(i)

              dir0(line)  = pdir(i)

              dspr0(line) = pspr(i)

              date0(line) = date1 + date2*1.0e-6

              ts(line) = tstep

              llat(line) = x

              llon(line) = y

              wndspd0(line) = wnd

              wnddir0(line) = wnddir


              line = line + 1

            END DO


          END IF

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

          !           Read file for ntint time levels

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

          ttemp = tstart

          DO WHILE (line.LE.readln)

            IF (flform) THEN

              READ (11,1000,END=112) date1,date2,x,y,numpart, &

                   wnd,wnddir,invar6,invar7

            ELSE

              READ (11,err=802,iostat=ioerr) datetime,x,y,    &

                   dummy,numpart,wnd,wnddir,invar5,invar6,invar7

              date1=dble(datetime(1))

              date2=dble(datetime(2))

            END IF


            ttest = date1 + date2*1.0e-6

            IF (ttest.GT.ttemp) THEN

              tstep = tstep+1

              ttemp = ttest

              IF (tstep.GT.ntint) EXIT

            END IF


            IF (flform) THEN

              READ (11,1010,END=112) invar1,invar2,invar3,invar4 ! Skip total integral parameters

              DO i = 1,numpart

                IF (line.LE.readln) THEN

                  READ (11,1010,END=112) hs0(line),tp0(line), &

                       dir0(line),dspr0(line)

                  date0(line) = ttest


                  ts(line) = tstep

                  llat(line) = x

                  llon(line) = y

                  wndspd0(line) = wnd

                  wnddir0(line) = wnddir


                  line = line+1

                END IF

              END DO

            ELSE

              READ (11,err=802,iostat=ioerr) k,invar1,invar2,  &

                   invar3,invar4,invar5

              DO i = 1,numpart

                IF (line.LE.readln) THEN

                  READ (11,END=112,ERR=802,IOSTAT=IOERR) k,         &

                       hs0(line),tp0(line),invar3,dir0(line),       &

                       dspr0(line)

                  date0(line) = ttest


                  ts(line) = tstep

                  llat(line) = x

                  llon(line) = y

                  wndspd0(line) = wnd

                  wnddir0(line) = wnddir


                  line = line+1

                END IF

              END DO

            END IF

          END DO

110       ierr = -1

          CLOSE(11)


112       CONTINUE

          IF (line.EQ.1) THEN

            WRITE(20,2002)

            WRITE(6,2002)

            stop 1

          END IF

          CLOSE(11)

          ! ===== READ LOOP FINISHED =====

          line=line-1


          WRITE(6,*) '... finished'


          IF (ttest.LT.tstart) THEN

            WRITE(20,2003) tstart

            WRITE(6,2003) tstart

            stop 1

          END IF


          IF (ALLOCATED(phs)) DEALLOCATE(phs)

          IF (ALLOCATED(ptp)) DEALLOCATE(ptp)

          IF (ALLOCATED(pdir)) DEALLOCATE(pdir)

          IF (ALLOCATED(pspr)) DEALLOCATE(pspr)

          IF (ALLOCATED(pwf)) DEALLOCATE(pwf)


#ifdef W3_MPI

        END IF

#endif

      END IF


#ifdef W3_MPI

      IF (rank.EQ.0) THEN

#endif


        !        Find unique time steps (and sort in ascending order)

        CALL unique(real(ts(1:line)),line,uniquetim,maxts)


        !        Find unique lat and lon values (and sort in ascending order)

        CALL unique(llat(1:line),SIZE(llat(1:line)),uniquelatraw,nout1)

        CALL unique(llon(1:line),SIZE(llon(1:line)),uniquelonraw,nout2)


        !--042916-----------------------

        !

        !        Redefine uniqueLatraw and uniqueLonrawto based on domain extent

        WRITE(20,*) 'uniqueLatraw(:) =', uniquelatraw(:)

        WRITE(20,*) 'uniqueLonraw(:) =', uniquelonraw(:)


        WRITE(20,*) 'No. increments: Longitude, Latitue =', mxcwt, mycwt

        DEALLOCATE(uniquelatraw)

        DEALLOCATE(uniquelonraw)

        ALLOCATE(uniquelatraw(mycwt+1))

        ALLOCATE(uniquelonraw(mxcwt+1))

        DO i = 1,(mycwt+1)

          uniquelatraw(i) =  minlat + &

               (real(i)-1)/real(mycwt)*(maxlat-minlat)

        END DO

        DO i = 1,(mxcwt+1)

          uniquelonraw(i) =  minlon + &

               (real(i)-1)/real(mxcwt)*(maxlon-minlon)

        END DO

        WRITE(20,*) 'uniqueLatraw(:) =', uniquelatraw(:)

        WRITE(20,*) 'uniqueLonraw(:) =', uniquelonraw(:)

        !

        !--042916-----------------------


        !        Filter out lats and lons outside of domain of interest

        DO latind1 = 1,SIZE(uniquelatraw)

          IF (uniquelatraw(latind1).GE.minlat) EXIT

        END DO

        DO latind2 = SIZE(uniquelatraw),1,-1

          IF (uniquelatraw(latind2).LE.maxlat) EXIT

        END DO

        DO lonind1 = 1,SIZE(uniquelonraw)

          IF (uniquelonraw(lonind1).GE.minlon) EXIT

        END DO

        DO lonind2 = SIZE(uniquelonraw),1,-1

          IF (uniquelonraw(lonind2).LE.maxlon) EXIT

        END DO

        WRITE(20,*) 'latind1, latind2, lonind1, lonind2 =', &

             latind1, latind2, lonind1, lonind2

        IF ((latind1.GE.latind2).OR.(lonind1.GE.lonind2)) THEN

          WRITE(20,1400)

          WRITE(6,1400)

          stop 1

        END IF

        NULLIFY(uniquelat)

        NULLIFY(uniquelon)

        ALLOCATE(uniquelat(latind2-latind1+1))

        ALLOCATE(uniquelon(lonind2-lonind1+1))

        uniquelat = uniquelatraw(latind1:latind2)

        uniquelon = uniquelonraw(lonind1:lonind2)

        WRITE(20,*) 'In waveTracking_NWS_V2: Longitude range =', &

             uniquelon(1), uniquelon(SIZE(uniquelon))

        WRITE(20,*) '                        Latitude range  =', &

             uniquelat(1), uniquelat(SIZE(uniquelat))


        !        Map is transposed (rotated by 90 deg), so that:

        !          I (matrix row) represents Longitute

        !          J (matrix column) represents Latitude

        !          i.e. from this point onwards the indices (i,j) represents Cart. coordinates

        ALLOCATE( mlon(SIZE(uniquelon),SIZE(uniquelat)) )

        ALLOCATE( mlat(SIZE(uniquelon),SIZE(uniquelat)) )

        !

        maxi = SIZE(uniquelon)

        maxj = SIZE(uniquelat)

        DO i = 1,maxi

          DO j = 1,maxj

            mlon(i,j) = uniquelon(i)

            mlat(i,j) = uniquelat(j)

          END DO

        END DO

#ifdef W3_MPI

      END IF

#endif


#ifdef W3_MPI

      CALL mpi_bcast(maxi,1,mpi_integer,0,mpi_comm_world,ierr)

      CALL mpi_bcast(maxj,1,mpi_integer,0,mpi_comm_world,ierr)

      CALL mpi_bcast(maxts,1,mpi_integer,0,mpi_comm_world,ierr)

#endif


      !        Allocate the wsdat structure

#ifdef W3_MPI

      IF (rank.EQ.0) THEN

#endif

        WRITE(20,*) 'Allocating wsdat...'

#ifdef W3_MPI

      END IF

#endif

      NULLIFY(wsdat)

      ALLOCATE(wsdat(maxts))

#ifdef W3_MPI

      IF (rank.EQ.0) THEN

#endif

        WRITE(20,*) 'SIZE(wsdat) = ',SIZE(wsdat)

#ifdef W3_MPI

      END IF

#endif


      !        Allocate and initialize the wsdat array

#ifdef W3_MPI

      IF (rank.EQ.0) THEN

#endif

        DO tsa = 1,maxts

          ALLOCATE(wsdat(tsa)%lat(maxi,maxj))

          ALLOCATE(wsdat(tsa)%lon(maxi,maxj))

          ALLOCATE(wsdat(tsa)%par(maxi,maxj))

          ALLOCATE(wsdat(tsa)%wnd(maxi,maxj))


          DO j = 1,maxj

            DO i = 1,maxi

              wsdat(tsa)%lat(i,j)=mlat(i,j)

              wsdat(tsa)%lon(i,j)=mlon(i,j)

              wsdat(tsa)%maxi=maxi

              wsdat(tsa)%maxj=maxj

              wsdat(tsa)%par(i,j)%hs(1:10)=9999.

              wsdat(tsa)%par(i,j)%tp(1:10)=9999.

              wsdat(tsa)%par(i,j)%dir(1:10)=9999.

              wsdat(tsa)%par(i,j)%dspr(1:10)=9999.

              !                  wsdat(tsA)%par(i,j)%wf(1:10)=9999.

              wsdat(tsa)%par(i,j)%ipart(1:10)=0

              wsdat(tsa)%par(i,j)%sys(1:10)=9999                      ! 40.PAR Increase this array, or make allocatable

              wsdat(tsa)%par(i,j)%ngbrSys(1:50)=9999

              wsdat(tsa)%wnd(i,j)%wdir=9999.

              wsdat(tsa)%wnd(i,j)%wspd=9999.

              wsdat(tsa)%par(i,j)%checked=-1

            END DO

          END DO

        END DO


        !        Assign to each line in partition file an entry in wsdat

        !        At each time step each point contains all numpart partitions.

        !        Only store the first 10 partitions.

        l = 1


        DO WHILE (l.LE.line)

          DO j = 1,maxj

            DO i = 1,maxi

              IF ( (abs(llat(l)-mlat(i,j)).LT.1.e-2).AND. &

                   (abs(llon(l)-mlon(i,j)).LT.1.e-2) ) THEN

                !                    WRITE(20,*) 'MATCHED! ',l,&

                !                       llat(l),mlat(i,j),ABS(llat(l)-mlat(i,j)),&

                !                       llon(l),mlon(i,j),ABS(llon(l)-mlon(i,j))

                wsdat(ts(l))%lat(i,j) = llat(l)

                wsdat(ts(l))%lon(i,j) = llon(l)

                !                   --- Find ALL partition values associated with

                !                       lat(i,j) and lon(i,j). Keep list index l

                !                       fixed and recycle iline as variable index. ---

                iline = l

                k = 1

                DO WHILE (                                         &

                     abs(wsdat(ts(l))%lat(i,j)-llat(iline)).LT.1.e-3 &

                     .AND.abs(wsdat(ts(l))%lon(i,j)-llon(iline)).LT.1.e-3 )

                  IF (k.LE.10) THEN

                    wsdat(ts(iline))%par(i,j)%ipart(k) = k

                    wsdat(ts(iline))%par(i,j)%hs(k)    = hs0(iline)

                    wsdat(ts(iline))%par(i,j)%tp(k)    = tp0(iline)

                    wsdat(ts(iline))%par(i,j)%dir(k)   = dir0(iline)

                    wsdat(ts(iline))%par(i,j)%dspr(k)  = dspr0(iline)

                    !                       wsdat(ts(k))%par(i,j)%wf(k) = wf0(l)

                    IF (k.EQ.1) THEN

                      wsdat(ts(iline))%date = date0(iline)

                      wsdat(ts(iline))%wnd(i,j)%wdir = wnddir0(iline)

                      wsdat(ts(iline))%wnd(i,j)%wspd = wndspd0(iline)

                      wsdat(ts(iline))%par(i,j)%checked = 0

                    END IF

                  END IF

                  k = k + 1

                  iline = iline + 1

                  if (iline.GT.line) EXIT

                END DO

                !                   --- Account for increment at the end of loop (400 CONTINUE)

                !                       and go one element back in list because of increment. ---

                l = iline-1

                GOTO 400

              END IF

            END DO

          END DO

400       CONTINUE

          IF (l+1.le.line) THEN

            IF (ts(l).LT.ts(l+1)) THEN

              k = line-l

              !             --- With each time step completed, deallocate processed 1:l

              !                 elements from 1d array. Create a temporary array size of

              !                 (l+1:line) with k elements and reallocate original array. ---

              IF (ALLOCATED(tmp_i4)) DEALLOCATE(tmp_i4)

              !             --- REALLOCATE(integer arrays) ---

              ALLOCATE(tmp_i4(k))

              tmp_i4(1:k) = ts((l+1):line)

              DEALLOCATE(ts)

              ALLOCATE(ts(k))

              ts(1:k) = tmp_i4(1:k)

              DEALLOCATE(tmp_i4)

              !             --- REALLOCATE(double precision arrays) ---

              IF (ALLOCATED(tmp_r8)) DEALLOCATE(tmp_r8)

              ALLOCATE(tmp_r8(k))

              tmp_r8(1:k) = date0((l+1):line)

              DEALLOCATE(date0)

              ALLOCATE(date0(k))

              date0(1:k) = tmp_r8(1:k)

              DEALLOCATE(tmp_r8)

              !             --- REALLOCATE(single precision arrays) ---

              IF (ALLOCATED(tmp_r4)) DEALLOCATE(tmp_r4)

              ALLOCATE(tmp_r4(k))

              tmp_r4(1:k) = llat((l+1):line)

              DEALLOCATE(llat)

              ALLOCATE(llat(k))

              llat(1:k) = tmp_r4(1:k)

              tmp_r4(1:k) = llon((l+1):line)

              DEALLOCATE(llon)

              ALLOCATE(llon(k))

              llon(1:k) = tmp_r4(1:k)

              tmp_r4(1:k) = hs0((l+1):line)

              DEALLOCATE(hs0)

              ALLOCATE(hs0(k))

              hs0(1:k) = tmp_r4(1:k)

              tmp_r4(1:k) = tp0((l+1):line)

              DEALLOCATE(tp0)

              ALLOCATE(tp0(k))

              tp0(1:k) = tmp_r4(1:k)

              tmp_r4(1:k) = dir0((l+1):line)

              DEALLOCATE(dir0)

              ALLOCATE(dir0(k))

              dir0(1:k) = tmp_r4(1:k)

              tmp_r4(1:k) = dspr0((l+1):line)

              DEALLOCATE(dspr0)

              ALLOCATE(dspr0(k))

              dspr0(1:k) = tmp_r4(1:k)

              tmp_r4(1:k) = wndspd0((l+1):line)

              DEALLOCATE(wndspd0)

              ALLOCATE(wndspd0(k))

              wndspd0(1:k) = tmp_r4(1:k)

              tmp_r4(1:k) = wnddir0((l+1):line)

              DEALLOCATE(wnddir0)

              ALLOCATE(wnddir0(k))

              wnddir0(1:k) = tmp_r4(1:k)

              DEALLOCATE(tmp_r4)

              line = k

              l = 0

            END IF

          END IF

          l = l + 1

        END DO

        !

        IF (ALLOCATED(ts)) DEALLOCATE(ts)

        IF (ALLOCATED(llat)) DEALLOCATE(llat)

        IF (ALLOCATED(llon)) DEALLOCATE(llon)

        IF (ALLOCATED(mlat)) DEALLOCATE(mlat)

        IF (ALLOCATED(mlon)) DEALLOCATE(mlon)

        IF (ALLOCATED(date0)) DEALLOCATE(date0)

        IF (ALLOCATED(hs0)) DEALLOCATE(hs0)

        IF (ALLOCATED(tp0)) DEALLOCATE(tp0)

        IF (ALLOCATED(dir0)) DEALLOCATE(dir0)

        IF (ALLOCATED(dspr0)) DEALLOCATE(dspr0)

        !         IF (ALLOCATED(wf0)) DEALLOCATE(wf0)

        IF (ALLOCATED(wndspd0)) DEALLOCATE(wndspd0)

        IF (ALLOCATED(wnddir0)) DEALLOCATE(wnddir0)

#ifdef W3_MPI

      END IF

#endif


#ifdef W3_MPI

      !         Communicate the wsdat entries from rank=0 to other ranks

      DO tsa = t0,maxts

        irank = mod((tsa-t0),min(nproc,maxts))

        !             WRITE(20,*) 'Rank,irank=',rank,irank

        IF (irank.NE.0) THEN

          !               WRITE(20,*) 'Communicating for Rank,irank=',rank,irank


          IF (rank.EQ.irank) THEN

            ALLOCATE(wsdat(tsa)%lat(maxi,maxj))

            ALLOCATE(wsdat(tsa)%lon(maxi,maxj))

            ALLOCATE(wsdat(tsa)%par(maxi,maxj))

            ALLOCATE(wsdat(tsa)%wnd(maxi,maxj))


            DO j = 1,maxj

              DO i = 1,maxi

                wsdat(tsa)%maxi=maxi

                wsdat(tsa)%maxj=maxj

                wsdat(tsa)%par(i,j)%hs(1:10)=9999.

                wsdat(tsa)%par(i,j)%tp(1:10)=9999.

                wsdat(tsa)%par(i,j)%dir(1:10)=9999.

                wsdat(tsa)%par(i,j)%dspr(1:10)=9999.

                wsdat(tsa)%par(i,j)%ipart(1:10)=0

                wsdat(tsa)%par(i,j)%sys(1:10)=9999                      ! 40.PAR Increase this array, or make allocatable

                wsdat(tsa)%par(i,j)%ngbrSys(1:50)=9999

                wsdat(tsa)%wnd(i,j)%wdir=9999.

                wsdat(tsa)%wnd(i,j)%wspd=9999.

                wsdat(tsa)%par(i,j)%checked=-1

              END DO

            END DO

          END IF


          DO j = 1,maxj

            DO i = 1,maxi

              tag1 = ((j-1)*maxi+i)*10


              IF (rank.EQ.0) THEN

                !                        WRITE(6,*) '>> Sending: rank,irank,tag1=', &

                !                                rank,irank,(tag1+1)

                commarr1 = (/wsdat(tsa)%par(i,j)%hs(:), &

                     wsdat(tsa)%par(i,j)%tp(:), &

                     wsdat(tsa)%par(i,j)%dir(:), &

                     wsdat(tsa)%par(i,j)%dspr(:), &

                     wsdat(tsa)%wnd(i,j)%wdir, &

                     wsdat(tsa)%wnd(i,j)%wspd, &

                     wsdat(tsa)%lat(i,j), &

                     wsdat(tsa)%lon(i,j)/)

                CALL mpi_send(commarr1,44,mpi_real,irank, &

                     (tag1+1),mpi_comm_world,ierr)

              END IF

              IF (rank.EQ.irank) THEN

                !                         WRITE(6,*) '<< Receiving: rank,irank,tag1=', &

                !                                rank,irank,(tag1+1)

                CALL mpi_recv(commarr1,44,mpi_real,0,(tag1+1), &

                     mpi_comm_world,mpi_status,ierr)

                wsdat(tsa)%par(i,j)%hs = commarr1(1:10)

                wsdat(tsa)%par(i,j)%tp = commarr1(11:20)

                wsdat(tsa)%par(i,j)%dir = commarr1(21:30)

                wsdat(tsa)%par(i,j)%dspr = commarr1(31:40)

                wsdat(tsa)%wnd(i,j)%wdir = commarr1(41)

                wsdat(tsa)%wnd(i,j)%wspd = commarr1(42)

                wsdat(tsa)%lat(i,j) = commarr1(43)

                wsdat(tsa)%lon(i,j) = commarr1(44)

              END IF


              IF (rank.EQ.0) THEN

                CALL mpi_send(wsdat(tsa)%date,1, &

                     mpi_double_precision,irank, &

                     (tag1+2),mpi_comm_world,ierr)

              END IF

              IF (rank.EQ.irank) THEN

                CALL mpi_recv(wsdat(tsa)%date,1, &

                     mpi_double_precision,0,(tag1+2), &

                     mpi_comm_world,mpi_status,ierr)

              END IF


              IF (rank.EQ.0) THEN

                !                        WRITE(6,*) '>> Sending: rank,irank,tag1=', &

                !                                rank,irank,(tag1+3)

                commarr2 = (/wsdat(tsa)%par(i,j)%ipart(:), &

                     wsdat(tsa)%par(i,j)%checked/)

                CALL mpi_send(commarr2,11, &

                     mpi_integer,irank,(tag1+3),mpi_comm_world,ierr)

              END IF

              IF (rank.EQ.irank) THEN

                !                         WRITE(6,*) '<< Receiving: rank,irank,tag1=', &

                !                                rank,irank,(tag1+3)

                CALL mpi_recv(commarr2,11, &

                     mpi_integer,0,(tag1+3), &

                     mpi_comm_world,mpi_status,ierr)

                wsdat(tsa)%par(i,j)%ipart(:) = commarr2(1:10)

                wsdat(tsa)%par(i,j)%checked = commarr2(11)

              END IF


            END DO

          END DO

        END IF

      END DO


      CALL mpi_barrier(mpi_comm_world,ierr)

#endif


#ifdef W3_MPI

      IF (rank.EQ.0) THEN

#endif

        ! ----*** Test Output *** --------------------------------------------------

        IF (testout) THEN

          !-----RAW PARTITION output: Coordinates

          OPEN(unit=31,file='PART_COORD.OUT',status='unknown')


          WRITE(31,*) 'Longitude ='

          DO j = maxj,1,-1

            DO i = 1,maxi

              WRITE(31,'(F7.2)',advance='NO') wsdat(1)%lon(i,j)

            END DO

            WRITE(31,'(A)',advance='YES') ''

          END DO


          WRITE(31,*) 'Latitude = '

          DO j = maxj,1,-1

            DO i = 1,maxi

              WRITE(31,'(F7.2)',advance='NO') wsdat(1)%lat(i,j)

            END DO

            WRITE(31,'(A)',advance='YES') ''

          END DO


          CLOSE(31)


          !-----RAW PARTITION output: hs

          OPEN(unit=32, file='PART_HSIGN.OUT', &

               status='unknown')


          maxpartout = 5

          DO tsa = 1,SIZE(wsdat)

            WRITE(32,'(I4,71x,A)') tsa,'Time step'

            WRITE(32,'(I4,71x,A)') maxpartout,'Tot number of raw partitions'

            DO k = 1,maxpartout

              WRITE(32,'(I4,71x,A)') k,'System number'

              WRITE(32,'(I4,71x,A)') 9999,'Number of points in system'

              DO j = maxj,1,-1

                DO i = 1,maxi

                  WRITE(32,'(F8.2)',advance='NO') wsdat(tsa)%par(i,j)%hs(k)

                END DO

                WRITE(32,'(A)',advance='YES') ''

              END DO

            END DO

          END DO


          CLOSE(32)


          !-----RAW PARTITION output: tp

          !      OPEN(unit=33,recl=2147483646, file='PART_TP.OUT', &

          !           status='unknown')

          OPEN(unit=33, file='PART_TP.OUT', &

               status='unknown')


          DO tsa = 1,SIZE(wsdat)

            WRITE(33,'(I4,71x,A)') tsa,'Time step'

            WRITE(33,'(I4,71x,A)') maxpartout,'Tot number of raw partitions'

            DO k = 1,maxpartout

              WRITE(33,'(I4,71x,A)') k,'System number'

              WRITE(33,'(I4,71x,A)') 9999,'Number of points in system'

              DO j = maxj,1,-1

                DO i = 1,maxi

                  WRITE(33,'(F8.2)',advance='NO') wsdat(tsa)%par(i,j)%tp(k)

                END DO

                WRITE(33,'(A)',advance='YES') ''

              END DO

            END DO

          END DO


          CLOSE(33)


          !-----RAW PARTITION output: dir

          OPEN(unit=34, file='PART_DIR.OUT', &

               status='unknown')


          DO tsa = 1,SIZE(wsdat)

            WRITE(34,'(I4,71x,A)') tsa,'Time step'

            WRITE(34,'(I4,71x,A)') maxpartout,'Tot number of raw partitions'

            DO k = 1,maxpartout

              WRITE(34,'(I4,71x,A)') k,'System number'

              WRITE(34,'(I4,71x,A)') 9999,'Number of points in system'

              DO j = maxj,1,-1

                DO i = 1,maxi

                  WRITE(34,'(F8.2)',advance='NO') wsdat(tsa)%par(i,j)%dir(k)

                END DO

                WRITE(34,'(A)',advance='YES') ''

              END DO

            END DO

          END DO


          CLOSE(34)


          !-----RAW PARTITION output: dspr

          OPEN(unit=35, file='PART_DSPR.OUT', &

               status='unknown')


          DO tsa = 1,SIZE(wsdat)

            WRITE(35,'(I4,71x,A)') tsa,'Time step'

            WRITE(35,'(I4,71x,A)') maxpartout,'Tot number of raw partitions'

            DO k = 1,maxpartout

              WRITE(35,'(I4,71x,A)') k,'System number'

              WRITE(35,'(I4,71x,A)') 9999,'Number of points in system'

              DO j = maxj,1,-1

                DO i = 1,maxi

                  WRITE(35,'(F8.2)',advance='NO') &

                       wsdat(tsa)%par(i,j)%dspr(k)

                END DO

                WRITE(35,'(A)',advance='YES') ''

              END DO

            END DO

          END DO


          CLOSE(35)

        END IF

#ifdef W3_MPI

      END IF

#endif


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


      !        Allocate the sysA structure

#ifdef W3_MPI

      IF (rank.EQ.0) THEN

#endif

        WRITE(20,*) 'Allocating sysA...'

#ifdef W3_MPI

      END IF

#endif

      ALLOCATE( sysa(maxts) )

#ifdef W3_MPI

      IF (rank.EQ.0) THEN

#endif

        WRITE(20,*) 'SIZE(sysA) = ',SIZE(sysa)

        WRITE(6,1020) ' Number of time levels being processed:',SIZE(sysa)

1020    FORMAT(a,i4)

#ifdef W3_MPI

      END IF

#endif


      !        Allocate maxSys

      ALLOCATE( maxsys(maxts) )

    ELSE

      !        Raw partitioning data from wave model memory, via the array wsdat.

      !        Set maxTs to the time step to compute: 1=first time step, 2=otherwise

      maxts = tmax

      t0 = tcur


      !        Allocate the sysA structure

      ALLOCATE( sysa(1) )                                              !Change to sysA(2)?

      !        Allocate maxSys

      ALLOCATE( maxsys(1) )                                            !Change to maxSys(2)?

    END IF


    !     Big loop over all time levels

#ifdef W3_MPI

    IF (rank.EQ.0) THEN

#endif

      WRITE(6,*) 'Performing spatial tracking...'

#ifdef W3_MPI

    END IF

    !      WRITE(20,*) 'rank,t0,maxTs,nproc =',rank,t0,maxTs,nproc

    DO tsa = (t0+rank),maxts,min(nproc,maxts)

      !         WRITE(20,*) 'Computing: Rank, tsA =',rank,tsA

#endif

#ifdef W3_SHRD

      DO tsa = t0,maxts

#endif


        WRITE(20,*) 'Call spiralTrackV3, tsA=',tsa,'...'

        CALL spiraltrackv3 ( wsdat(tsa), dirknob, perknob, wetpts, &

             hsknob, seedlat, seedlon, &

             maxsys(tsa), sysa(tsa)%sys )


        WRITE(20,*) '*** SIZE(sysA(1:tsA)%sys) at end of time step', &

             tsa,':'

        WRITE(20,*) SIZE(sysa(tsa)%sys)

#ifdef W3_SHRD

      END DO

#endif

#ifdef W3_MPI

    END DO

#endif


#ifdef W3_MPI

    CALL mpi_barrier(mpi_comm_world,ierr)


    !!     Define communicator for array of integers in structure "system"

    !      DOMSIZE = maxI*maxJ

    !      WRITE(20,*) 'Rank',rank,'DOMSIZE =',DOMSIZE

    !      CALL MPI_TYPE_CONTIGUOUS(DOMSIZE,MPI_INTEGER,MPI_INT_DOMARR,IERR)

    !      CALL MPI_TYPE_COMMIT(MPI_INT_DOMARR,IERR)

    !      CALL MPI_TYPE_EXTENT(MPI_INT_DOMARR,EXTENT,IERR)

    !      WRITE(20,*) 'Rank',rank,'has set up communicator MPI_INT_DOMARR, &

    !                  size =',EXTENT


    !!     Define communicator for array of reals in structure "system"

    !      CALL MPI_TYPE_CONTIGUOUS(DOMSIZE,MPI_REAL,MPI_REAL_DOMARR,IERR)

    !      CALL MPI_TYPE_COMMIT(MPI_REAL_DOMARR,IERR)

    !      CALL MPI_TYPE_EXTENT(MPI_REAL_DOMARR,EXTENT,IERR)

    !      WRITE(20,*) 'Rank',rank,'has set up communicator MPI_REAL_DOMARR, &

    !                  size =',EXTENT


    !     Communicate results back to rank 0

    DO tsa = t0,maxts

      irank = mod((tsa-t0),min(nproc,maxts))

      !         WRITE(20,*) 'Rank,irank=',rank,irank

      IF (irank.NE.0) THEN

        !            WRITE(20,*) 'Communicating for Rank,irank=',rank,irank


        !           Send maxSys(tsA) at each time level to rank 0

        tag1 = tsa

        IF (rank.EQ.irank) THEN

          !              Send results from current rank to rank 0 (blocking)

          !               WRITE(20,*) '>> Sending: rank,tsA,tag1=',rank,tsA,tag1

          CALL mpi_send(maxsys(tsa),1,mpi_integer,0,tag1, &

               mpi_comm_world,ierr)

          !               WRITE(20,*) 'Rank, IERR=',rank,IERR

        END IF

        IF (rank.EQ.0) THEN

          !               WRITE(20,*) '<< Receiving: rank,tsA,tag1=',rank,tsA,tag1

          CALL mpi_recv(maxsys(tsa),1,mpi_integer, &

               irank,tag1,mpi_comm_world,mpi_status,ierr)

          !              Allocate structure at this time level

          ALLOCATE( sysa(tsa)%sys(maxsys(tsa)) )

          DO ic = 1,maxsys(tsa)

            NULLIFY( sysa(tsa)%sys(ic)%i )

            NULLIFY( sysa(tsa)%sys(ic)%j )

            NULLIFY( sysa(tsa)%sys(ic)%lon )

            NULLIFY( sysa(tsa)%sys(ic)%lat )

            NULLIFY( sysa(tsa)%sys(ic)%hs )

            NULLIFY( sysa(tsa)%sys(ic)%tp )

            NULLIFY( sysa(tsa)%sys(ic)%dir)

            NULLIFY( sysa(tsa)%sys(ic)%dspr)

            ALLOCATE( sysa(tsa)%sys(ic)%i(maxi*maxj) )

            ALLOCATE( sysa(tsa)%sys(ic)%j(maxi*maxj) )

            ALLOCATE( sysa(tsa)%sys(ic)%lon(maxi*maxj) )

            ALLOCATE( sysa(tsa)%sys(ic)%lat(maxi*maxj) )

            ALLOCATE( sysa(tsa)%sys(ic)%hs(maxi*maxj) )

            ALLOCATE( sysa(tsa)%sys(ic)%tp(maxi*maxj) )

            ALLOCATE( sysa(tsa)%sys(ic)%dir(maxi*maxj) )

            ALLOCATE( sysa(tsa)%sys(ic)%dspr(maxi*maxj) )

            sysa(tsa)%sys(ic)%i(:) = 9999

            sysa(tsa)%sys(ic)%j(:) = 9999

            sysa(tsa)%sys(ic)%lon(:) = 9999.

            sysa(tsa)%sys(ic)%lat(:) = 9999.

            sysa(tsa)%sys(ic)%hs(:) = 9999.

            sysa(tsa)%sys(ic)%tp(:) = 9999.

            sysa(tsa)%sys(ic)%dir(:) = 9999.

            sysa(tsa)%sys(ic)%dspr(:) = 9999.

            sysa(tsa)%sys(ic)%hsMean = 9999.

            sysa(tsa)%sys(ic)%tpMean = 9999.

            sysa(tsa)%sys(ic)%dirMean = 9999.

            sysa(tsa)%sys(ic)%sysInd = 9999

            sysa(tsa)%sys(ic)%nPoints = 9999

            sysa(tsa)%sys(ic)%grp = 9999

          END DO

        END IF


        !           Send data fields at each (tsA,ic) combination

        IF ((rank.EQ.0).OR.(rank.EQ.irank)) THEN

          DO ic = 1, maxsys(tsa)

            !              Construct a unique tag for each message

            tag2 = tsa*10000 + ic*100

            domsize = maxi*maxj


            IF (rank.EQ.irank) THEN

              !                 WRITE(20,*) '>> Sending: rank,irank,tag2=', &

              !                             rank,irank,(tag2+1)

              CALL mpi_send(sysa(tsa)%sys(ic)%i(:),domsize, &

                   mpi_integer,0,(tag2+1),mpi_comm_world,req(1),ierr)

            END IF

            IF (rank.EQ.0) THEN

              !                 WRITE(20,*) '<< Receiving: rank,irank,tag2=', &

              !                             rank,irank,(tag2+1)

              CALL mpi_recv(sysa(tsa)%sys(ic)%i(:),domsize, &

                   mpi_integer,irank,(tag2+1), &

                   mpi_comm_world,mpi_status,req(2),ierr)

            END IF

            !               CALL MPI_WAITALL(2,REQ,ISTAT,IERR)


            IF (rank.EQ.irank) THEN

              !                 WRITE(20,*) '>> Sending: rank,irank,tag2=', &

              !                             rank,irank,(tag2+2)

              CALL mpi_send(sysa(tsa)%sys(ic)%j(:),domsize, &

                   mpi_integer,0,(tag2+2),mpi_comm_world,req(1),ierr)

            END IF

            IF (rank.EQ.0) THEN

              !                 WRITE(20,*) '<< Receiving: rank,irank,tag2=', &

              !                             rank,irank,(tag2+2)

              CALL mpi_recv(sysa(tsa)%sys(ic)%j(:),domsize, &

                   mpi_integer,irank,(tag2+2), &

                   mpi_comm_world,mpi_status,req(2),ierr)

            END IF

            !               CALL MPI_WAITALL(2,REQ,ISTAT,IERR)


            IF (rank.EQ.irank) THEN

              !                 WRITE(20,*) '>> Sending: rank,tag2=',rank,(tag2+3)

              CALL mpi_send(sysa(tsa)%sys(ic)%lon(:),domsize, &

                   mpi_real,0,(tag2+3),mpi_comm_world,req(1),ierr)

            END IF

            IF (rank.EQ.0) THEN

              !                 WRITE(20,*) '<< Receiving: rank,tag2=',rank,(tag2+3)

              CALL mpi_recv(sysa(tsa)%sys(ic)%lon(:),domsize, &

                   mpi_real,irank,(tag2+3), &

                   mpi_comm_world,mpi_status,req(2),ierr)

            END IF

            !               CALL MPI_WAITALL(2,REQ,ISTAT,IERR)


            IF (rank.EQ.irank) THEN

              !                 WRITE(20,*) '>> Sending: rank,tag2=',rank,(tag2+4)

              CALL mpi_send(sysa(tsa)%sys(ic)%lat(:),domsize, &

                   mpi_real,0,(tag2+4),mpi_comm_world,req(1),ierr)

            END IF

            IF (rank.EQ.0) THEN

              !                 WRITE(20,*) '<< Receiving: rank,tag2=',rank,(tag2+4)

              CALL mpi_recv(sysa(tsa)%sys(ic)%lat(:),domsize, &

                   mpi_real,irank,(tag2+4), &

                   mpi_comm_world,mpi_status,req(2),ierr)

            END IF

            !               CALL MPI_WAITALL(2,REQ,ISTAT,IERR)


            IF (rank.EQ.irank) THEN

              !                 WRITE(20,*) '>> Sending: rank,tag2=',rank,(tag2+5)

              CALL mpi_send(sysa(tsa)%sys(ic)%hs(:),domsize, &

                   mpi_real,0,(tag2+5),mpi_comm_world,req(1),ierr)

            END IF

            IF (rank.EQ.0) THEN

              !                 WRITE(20,*) '<< Receiving: rank,tag2=',rank,(tag2+5)

              CALL mpi_recv(sysa(tsa)%sys(ic)%hs(:),domsize, &

                   mpi_real,irank,(tag2+5), &

                   mpi_comm_world,mpi_status,req(2),ierr)

            END IF

            !               CALL MPI_WAITALL(2,REQ,ISTAT,IERR)


            IF (rank.EQ.irank) THEN

              !                 WRITE(20,*) '>> Sending: rank,tag2=',rank,(tag2+6)

              CALL mpi_send(sysa(tsa)%sys(ic)%tp(:),domsize, &

                   mpi_real,0,(tag2+6),mpi_comm_world,req(1),ierr)

            END IF

            IF (rank.EQ.0) THEN

              !                 WRITE(20,*) '<< Receiving: rank,tag2=',rank,(tag2+6)

              CALL mpi_recv(sysa(tsa)%sys(ic)%tp(:),domsize, &

                   mpi_real,irank,(tag2+6), &

                   mpi_comm_world,mpi_status,req(2),ierr)

            END IF

            !               CALL MPI_WAITALL(2,REQ,ISTAT,IERR)


            IF (rank.EQ.irank) THEN

              !                 WRITE(20,*) '>> Sending: rank,tag2=',rank,(tag2+7)

              CALL mpi_send(sysa(tsa)%sys(ic)%dir(:),domsize, &

                   mpi_real,0,(tag2+7),mpi_comm_world,req(1),ierr)

            END IF

            IF (rank.EQ.0) THEN

              !                 WRITE(20,*) '<< Receiving: rank,tag2=',rank,(tag2+7)

              CALL mpi_recv(sysa(tsa)%sys(ic)%dir(:),domsize, &

                   mpi_real,irank,(tag2+7), &

                   mpi_comm_world,mpi_status,req(2),ierr)

            END IF

            !               CALL MPI_WAITALL(2,REQ,ISTAT,IERR)


            IF (rank.EQ.irank) THEN

              !                 WRITE(20,*) '>> Sending: rank,tag2=',rank,(tag2+8)

              CALL mpi_send(sysa(tsa)%sys(ic)%dspr(:),domsize, &

                   mpi_real,0,(tag2+8),mpi_comm_world,req(1),ierr)

            END IF

            IF (rank.EQ.0) THEN

              !                 WRITE(20,*) '<< Receiving: rank,tag2=',rank,(tag2+8)

              CALL mpi_recv(sysa(tsa)%sys(ic)%dspr(:),domsize, &

                   mpi_real,irank,(tag2+8), &

                   mpi_comm_world,mpi_status,req(2),ierr)

            END IF

            !               CALL MPI_WAITALL(2,REQ,ISTAT,IERR)


            IF (rank.EQ.irank) THEN

              !                 WRITE(20,*) '>> Sending: rank,irank,tag2=', &

              !                             rank,irank,(tag2+9)

              CALL mpi_send(sysa(tsa)%sys(ic)%hsMean,1,mpi_real, &

                   0,(tag2+9),mpi_comm_world,ierr)

            END IF

            IF (rank.EQ.0) THEN

              !                 WRITE(20,*) '<< Receiving: rank,irank,tag2=', &

              !                             rank,irank,(tag2+9)

              CALL mpi_recv(sysa(tsa)%sys(ic)%hsMean,1,mpi_real, &

                   irank,(tag2+9),mpi_comm_world,mpi_status,ierr)

            END IF


            IF (rank.EQ.irank) THEN

              !                 WRITE(20,*) '>> Sending: rank,irank,tag2=', &

              !                             rank,irank,(tag2+10)

              CALL mpi_send(sysa(tsa)%sys(ic)%tpMean,1,mpi_real, &

                   0,(tag2+10),mpi_comm_world,ierr)

            END IF

            IF (rank.EQ.0) THEN

              !                 WRITE(20,*) '<< Receiving: rank,irank,tag2=', &

              !                             rank,irank,(tag2+10)

              CALL mpi_recv(sysa(tsa)%sys(ic)%tpMean,1,mpi_real, &

                   irank,(tag2+10),mpi_comm_world,mpi_status,ierr)

            END IF


            IF (rank.EQ.irank) THEN

              !                 WRITE(20,*) '>> Sending: rank,irank,tag2=', &

              !                             rank,irank,(tag2+11)

              CALL mpi_send(sysa(tsa)%sys(ic)%dirMean,1,mpi_real, &

                   0,(tag2+11),mpi_comm_world,ierr)

            END IF

            IF (rank.EQ.0) THEN

              !                 WRITE(20,*) '<< Receiving: rank,irank,tag2=', &

              !                             rank,irank,(tag2+11)

              CALL mpi_recv(sysa(tsa)%sys(ic)%dirMean,1,mpi_real, &

                   irank,(tag2+11),mpi_comm_world,mpi_status,ierr)

            END IF


            IF (rank.EQ.irank) THEN

              !                 WRITE(20,*) '>> Sending: rank,irank,tag2=', &

              !                             rank,irank,(tag2+12)

              CALL mpi_send(sysa(tsa)%sys(ic)%sysInd,1,mpi_integer,&

                   0,(tag2+12),mpi_comm_world,ierr)

            END IF

            IF (rank.EQ.0) THEN

              !                 WRITE(20,*) '<< Receiving: rank,irank,tag2=', &

              !                             rank,irank,(tag2+12)

              CALL mpi_recv(sysa(tsa)%sys(ic)%sysInd,1,mpi_integer,&

                   irank,(tag2+12),mpi_comm_world,mpi_status,ierr)

            END IF


            IF (rank.EQ.irank) THEN

              !                 WRITE(20,*) '>> Sending: rank,irank,tag2=', &

              !                             rank,irank,(tag2+13)

              CALL mpi_send(sysa(tsa)%sys(ic)%nPoints,1,mpi_integer,&

                   0,(tag2+13),mpi_comm_world,ierr)

            END IF

            IF (rank.EQ.0) THEN

              !                 WRITE(20,*) '<< Receiving: rank,irank,tag2=', &

              !                             rank,irank,(tag2+13)

              CALL mpi_recv(sysa(tsa)%sys(ic)%nPoints,1,mpi_integer,&

                   irank,(tag2+13),mpi_comm_world,mpi_status,ierr)

            END IF


            IF (rank.EQ.irank) THEN

              !                 WRITE(20,*) '>> Sending: rank,irank,tag2=', &

              !                             rank,irank,(tag2+14)

              CALL mpi_send(sysa(tsa)%sys(ic)%grp,1,mpi_integer,&

                   0,(tag2+14),mpi_comm_world,ierr)

            END IF

            IF (rank.EQ.0) THEN

              !                 WRITE(20,*) '<< Receiving: rank,irank,tag2=', &

              !                             rank,irank,(tag2+14)

              CALL mpi_recv(sysa(tsa)%sys(ic)%grp,1,mpi_integer,&

                   irank,(tag2+14),mpi_comm_world,mpi_status,ierr)

            END IF

          END DO

        END IF

      END IF

    END DO


    CALL mpi_barrier(mpi_comm_world,ierr)


    !      CALL MPI_TYPE_FREE(MPI_INT_DOMARR,IERR)

    !      CALL MPI_TYPE_FREE(MPI_REAL_DOMARR,IERR)


#endif


#ifdef W3_MPI

    IF (rank.EQ.0) THEN

#endif

      WRITE(6,*) 'Performing temporal tracking...'

      WRITE(20,*) 'Calling timeTrackingV2...'

      lonext = wsdat(1)%lon(maxi,1)-wsdat(1)%lon(1,1)

      latext = wsdat(1)%lat(1,maxj)-wsdat(1)%lat(1,1)


      CALL timetrackingv2 (sysa, maxsys, tptimeknob, dirtimeknob, 1, &

           maxgroup, dt, lonext, latext, maxi, maxj)

      !

#ifdef W3_MPI

    END IF

#endif

    !

    RETURN

    !

802 CONTINUE

    WRITE (6,990) ioerr

    stop 1


990 FORMAT (/' *** WAVEWATCH III ERROR IN W3STRKMD : '/            &

         '     ERROR IN READING FROM PARTITION FILE'/          &

         '     IOSTAT =',i5/)

1000 FORMAT (f9.0,f7.0,f8.3,f8.3,14x,i3,7x,f5.1,f6.1,f5.1,f6.1)

1010 FORMAT (3x,f8.2,f8.2,8x,f9.2,f9.2)

1200 FORMAT (/' *** WAVEWATCH III ERROR IN W3STRKMD : '/            &

         '     ERROR IN READING PARTITION FILE   '/            &

         '     INCOMPATIBLE ENDIANESS'/                        )

1300 FORMAT (/' *** WAVEWATCH III ERROR IN W3STRKMD : '/            &

         '     ERROR IN READING PARTITION FILE   '/                   &

         '     EXPECTED IDSTR "WAVEWATCH III PARTITIONED DATA FILE"'/ )

1400 FORMAT (/' *** WAVEWATCH III ERROR IN W3STRKMD : '/            &

         '     ERROR IN FINDING DOMAIN TO PROCESS - '/         &

         '     SPECIFIED LAT/LON LIMITS WITHIN DOMAIN '/       &

         '     OF RAW PARTITION FILE?'/                        )

2001 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/            &

         '     ERROR IN OPENING INPUT FILE'/                   )

2002 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/            &

         '     PREMATURE END OF INPUT FILE'/                   )

2003 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/            &

         '     PREMATURE END OF PARTITION FILE - '/            &

         '     TSTART=',f13.4/                                 )

    !

    !

  END SUBROUTINE wavetracking_nws_v2

  !/ End of waveTracking_NWS_V2 ---------------------------------------- /

  !/

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

  SUBROUTINE spiraltrackv3 (wsdat   ,dirKnob ,perKnob ,wetPts  , &

       hsKnob  ,seedLat ,seedLon  , &

       maxSys  ,sys   )

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          4-Jan-2013 |

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

    !/

    !/    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 )

    !/

    !/    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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     Performs the spatial spiral tracking for a given time step

    !

    !  2. Method

    !

    !     Index convention on grid:

    !

    !      j

    !      ^

    !      |+(1,maxJ)          +(maxI,maxJ)

    !      |

    !      |

    !      |

    !      |

    !      |

    !      |

    !      |(1,1)              +(maxI,1)

    !      +----------------------> i

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     dirKnob     Real     input  Parameter in direction for combining fields in space

    !     perKnob     Real     input  Parameter in period for combining fields in space

    !     wetPts      Real     input  Percentage of wet points for purging fields (fraction)

    !     hsKnob      Real     input  Parameter in wave height for purging fields

    !     seedLat     Real     input  Start Lat for tracking spiral (if =0 centre of field is used)

    !     seedLon     Real     input  Start Lon for tracking spiral (if =0 centre of field is used)

    !     wsdat       Real arr output Input 2d (gridded) data structure to be spiral tracked

    !     maxSys      Int      output Maximum number of partition systems

    !     sys     Type(system) output Final set of tracked systems, for one time level

    !

    TYPE(dat2d)  :: wsdat

    REAL         :: dirKnob,perKnob,wetPts,hsKnob,seedLat,seedLon

    INTEGER      :: maxSys

    TYPE(system), POINTER :: sys(:)


    INTENT (IN) wetpts,dirknob,perknob,hsknob,seedlat,seedlon

    INTENT (IN OUT) wsdat

    !      INTENT (OUT) maxSys,sys

    !

    !     Local variables

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

    !     ngbrExt    Int    How far do we want the neighbour to be considered

    !     combine    Int    Toggle (1=combine wave systems; 0=do not combine)

    !     maxI,MaxJ  Int    Dimensions of the 2d (gridded) data wsdat

    !     deltaLat   Real   Delta in kilometers between 2 pts (in latitude)

    !

    LOGICAL   :: first

    CHARACTER :: way *1

    INTEGER   :: ngbrExt, combine, maxI, maxJ, i, j, oldJ

    INTEGER   :: horizStepCount, vertStepCount, checkCount, sc, &

         maxPts, landPts, horizBorder, vertBorder, m, k, &

         stepCount

    REAL      :: deltaLat, minLat, maxLat, minLon, maxLon

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      findWay

    !      findNext

    !      findSys

    !      combineWaveSystems

    !

    !  5. Subroutines calling

    !

    !     waveTracking_NWS_V2

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     -

    !

    !  9. Switches :

    !

    !     None defined yet.

    !

    ! 10. Source code :

    !

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


    ! Routine starts by identifying the starting point. Choose the 'center' of the domain

    !     Set the search distance for neighbors:

    !     1: 1 row and column out, i.e. the 8 neighbors around the current point

    !     2: 2 rows and columns out... etc.

    ngbrext=1

    combine=1

    WRITE(20,*) 'In spiralTrackV3: combine = ',combine


    maxi = wsdat%maxi

    maxj = wsdat%maxj

    IF ( (seedlat.EQ.0).OR.(seedlon.EQ.0) ) THEN

      i=nint(real(maxi)/2.)

      j=nint(real(maxj)/2.)

      WRITE(20,*) 'In spiralTrackV3, i=NINT(maxI/2.) =',i

      WRITE(20,*) 'In spiralTrackV3, j=NINT(maxJ/2.) =',j

    ELSE

      i=1

      j=1

      DO WHILE ( (wsdat%lat(1,j).LT.seedlat).AND.(j.LT.wsdat%maxj) )     !40.PAR !Improve with SWAN's indice identification

        j=j+1

      END DO

      DO WHILE ( (wsdat%lon(i,1).LT.seedlon).AND.(i.LT.wsdat%maxi) )

        i=i+1

      END DO

    END IF

    !     In case center point is land point...

    IF (wsdat%par(i,j)%checked.EQ.-1) THEN

      oldj=j

      DO WHILE (wsdat%par(i,j)%checked.EQ.-1)

        j=j+1

        IF (j.EQ.maxj) THEN

          j=oldj

          i=i+1

          oldj=oldj+1

        END IF

      END DO

    END IF

    !     Compute distance in km between 2 grid points (at equator)

    deltalat=(wsdat%lat(i,j)-wsdat%lat(i,j-1))*111.18


    !     Starts the spiral

    !     Intitiate variables

    horizstepcount=0

    vertstepcount=0

    way='R'

    first=.true.

    checkcount=1

    maxsys=0

    landpts=0


    minlat=minval(wsdat%lat)

    maxlat=maxval(wsdat%lat)

    minlon=minval(wsdat%lon)

    maxlon=maxval(wsdat%lon)


    horizborder=0

    vertborder=0

    DO WHILE (checkcount.LE.(maxi*maxj-3) )

      !        From the direction (way) we were going before, find which direction we

      !        are going now and how many 'step' we need to take

      CALL findway(way, horizstepcount, vertstepcount, &

           vertborder, horizborder, stepcount)

      IF (first) THEN

        m=0

        DO k=1,length(wsdat%par(i,j)%hs, &

             SIZE(wsdat%par(i,j)%hs),9999.)

          IF ( (wsdat%par(i,j)%hs(k).EQ.0.).AND. &

               (wsdat%par(i,j)%tp(k).EQ.0.) ) THEN

            wsdat%par(i,j)%sys(k)=-1

          ELSE

            m=m+1

            wsdat%par(i,j)%sys(k)=m

          END IF

        END DO


        wsdat%par(i,j)%checked=1

        checkcount=checkcount+1

        first=.false.

      END IF

      DO sc = 1, stepcount

        CALL findnext (i,j,maxi,maxj,way,vertborder,horizborder)

        IF ( wsdat%par(i,j)%checked.EQ.-1 ) THEN

          !            Land point is one of our grid points, so we need to update counter

          checkcount=checkcount+1

          landpts=landpts+1

          !                So that we don't count the land points twice....

          wsdat%par(i,j)%checked=-2

        ELSE IF ( wsdat%par(i,j)%checked.EQ.0 ) THEN

          !            Hasn't been checked yet and is not land point

          checkcount=checkcount+1

          CALL findsys(i, j, wsdat, maxsys, ngbrext, maxi, maxj, &

               perknob, dirknob, hsknob)

        END IF

      END DO

    END DO

    !     wetPts% of wet points

    maxpts=nint(wetpts*(maxi*maxj-1))

    !

    WRITE(20,*) 'Call combineWaveSystems...'

    CALL combinewavesystems(wsdat,maxsys,maxpts,maxi,maxj, &

         perknob,dirknob,hsknob,combine,sys)


    RETURN

  END SUBROUTINE spiraltrackv3

  !/ End of spiralTrackV3 ---------------------------------------------- /

  !/

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

  SUBROUTINE timetrackingv2 (sysA       ,maxSys     ,tpTimeKnob  , &

       dirTimeKnob,ts0        ,maxGroup    , &

       dt         ,lonext     ,latext      , &

       maxI       ,maxJ                    )

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          4-Jan-2013 |

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

    !/

    !/    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 )

    !/

    !/    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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     Performs the time tracking of the systems identified within

    !     the subroutine spiralTrackV3.

    !

    !  2. Method

    !

    !     -

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     Note: perKnob, dirKnob in Matlab version replaced by tpTimeKnob, dirTimeKnob!

    !

    !     sysA        TYPE(timsys) in/out Final set of spatially and temporally tracked systems

    !     dirTimeKnob    Real      input  Parameter in direction for combining fields in time

    !     tpTimeKnob     Real      input  Parameter in period for combining fields in time

    !     ts0            Int       input  Time step to which default grp values are associated

    !     maxSys         Int arr   input  Total number of systems per time level

    !     maxGroup       Int       output Maximum number of wave systems ("groups") tracked in time

    !     lonext         Real      input  Longitudinal extent of domain

    !     latext         Real      input  Latitudinal extent of domain

    !     maxI, maxJ     Int       input  Maximum indices of wave field

    !

    TYPE(timsys), POINTER :: sysA(:)

    INTEGER, POINTER :: maxSys(:)

    REAL         :: dirTimeKnob, tpTimeKnob

    INTEGER      :: ts0, maxGroup

    REAL         :: dt

    REAL         :: lonext, latext

    INTEGER      :: maxI, maxJ


    INTENT (IN) tptimeknob, dirtimeknob, ts0, maxi, maxj

    !      INTENT (IN OUT) sysA

    INTENT (OUT) maxgroup

    !

    !     Local variables

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

    !     ic        Int   Counter for wave systems

    !     ts1       Int   Adjusted initial time step in case ts0 has only empty systems

    !

    LOGICAL :: file_exists

    CHARACTER :: dummy*23

    TYPE(sysmemory) :: sysMem(50)                                             !!! 50 memory spaces should be enough Check!!!

    INTEGER :: leng, l, i, ii, j, k, kk, idir, numSys, &

         counter, new, DIFSIZE, tpMinInd, dirMinInd, used, ok

    REAL    :: Tb,  deltaPer, deltaDir, tpMinVal, dirMinVal, &

         dirForTpMin, tpForDirMin

    REAL, ALLOCATABLE :: sysOrdered(:), TEMP(:), dirs(:)

    REAL, POINTER :: DIFARR(:)

    INTEGER, ALLOCATABLE :: indSorted(:), alreadyUsed(:), allInd(:)

    INTEGER, ALLOCATABLE :: ind(:), ind2(:)

    INTEGER :: ts1

    REAL, ALLOCATABLE :: GOF(:,:), GOFMinVal(:), GOFMinInd(:), &

         Tbsysmem(:), deltaDirsysmem(:), &

         deltaPersysmem(:),m1sysmem(:),m2sysmem(:)

    REAL    :: m1, m2

    REAL    :: lonmean, latmean, dmndiag

    INTEGER :: npnts, npnts2

    REAL, ALLOCATABLE :: mnlonlist(:), mnlatlist(:), mndist(:)

    REAL, POINTER     :: dummy1(:),dummy2(:),dummy3(:)

    INTEGER, ALLOCATABLE :: olsize(:)

    REAL    :: TEMP1, TEMP2

    INTEGER :: iii, jj, ll, idup

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      SORT

    !      SETDIFF

    !

    !  5. Subroutines calling

    !

    !     waveTracking_NWS_V2

    !

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     -

    !

    !  9. Switches :

    !

    !     None defined yet.

    !

    ! 10. Source code :

    !

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


    !     Associate default grp value to time step 1

    WRITE(20,*) 'TIME TRACKING'

    WRITE(20,*) 'Inside timeTrackingV2: SIZE(sysA(1)%sys) =', &

         SIZE(sysa(1)%sys)

    WRITE(20,*) 'Inside timeTrackingV2: maxSys(1) =',maxsys(1)

    WRITE(20,*) 'ts0 = ',ts0


    ts1 = ts0


    !     Skip initial time steps with empty systems (e.g. when starting from rest)

    DO i = ts1, SIZE(sysa)

      IF (SIZE(sysa(ts1)%sys).EQ.0) ts1 = ts1+1

      !        No non-empty systems found

      IF (ts1.GT.SIZE(sysa)) THEN

        maxgroup = 0

        GOTO 2000

      END IF

    END DO

    WRITE(20,*) 'TS = ',ts1


    IF (SIZE(sysa(ts1)%sys).GT.0) THEN

      !        Initialize system memory groups

      sysa(ts1)%sys(:)%grp = 9999

      sysmem(:)%grp = 9999

      sysmem(:)%nPoints = 0

      sysmem(:)%lonMean = 9999.

      sysmem(:)%latMean = 9999.

      sysmem(:)%tpMean = 9999.

      sysmem(:)%dirMean = 9999.

      sysmem(:)%updated = -9999

      sysmem(:)%length = 0

      DO iii = 1,50

        ALLOCATE(sysmem(iii)%indx(maxi*maxj))

        sysmem(iii)%indx = 9999

      END DO


      INQUIRE(file="sys_restart.ww3", exist=file_exists)

      IF (file_exists) THEN

        !           Use groups from wave tracking hotfile

        WRITE(20,*) '*** Using group memory hotfile'

        OPEN(unit=12,file='sys_restart.ww3',status='old')

        READ(12,'(A23,I10)') dummy,maxgroup

        WRITE(20,*) 'Reading ',maxgroup,' systems'

        DO k = 1,maxgroup

          READ(12,'(A23,I10)') dummy,sysmem(k)%grp

          READ(12,'(A23,I10)') dummy,sysmem(k)%nPoints

          READ(12,'(A23,F10.4)') dummy,sysmem(k)%lonMean

          READ(12,'(A23,F10.4)') dummy,sysmem(k)%latMean

          READ(12,'(A23,F10.3)') dummy,sysmem(k)%tpMean

          READ(12,'(A23,F10.3)') dummy,sysmem(k)%dirMean

          READ(12,'(A23,I10)') dummy,sysmem(k)%updated

          READ(12,'(A23,I10)') dummy,sysmem(k)%length

          DO j = maxj,1,-1

            READ(12,*) (sysmem(k)%indx((j-1)*maxi+i), i = 1,maxi)

          END DO

          !Reset update counter

          sysmem(k)%updated = 0

        END DO

        CLOSE(12)

        ts1 = ts1-1

      ELSE

        !        Set up the group number array for the first time level to be tracked

        ALLOCATE( sysordered(maxsys(ts1)) )

        ALLOCATE( indsorted(maxsys(ts1)) )

        CALL sort (real(sysa(ts1)%sys(1:maxsys(ts1))%nPoints), &

             maxsys(ts1),sysordered,indsorted,'D')

        sysa(ts1)%sys(1:maxsys(ts1)) = sysa(ts1)%sys(indsorted)

        IF (ALLOCATED(sysordered)) DEALLOCATE(sysordered)

        IF (ALLOCATED(indsorted)) DEALLOCATE(indsorted)


        !        Set the initial long-term system memory

        DO i = 1, maxsys(ts1)

          sysa(ts1)%sys(i)%grp = i

          !           Set initial values of long-term system memory

          sysmem(i)%grp = i

          sysmem(i)%nPoints = sysa(ts1)%sys(i)%nPoints

          sysmem(i)%lonMean = &

               sum(sysa(ts1)%sys(i)%lon(1:sysmem(i)%nPoints))/&

               sysmem(i)%nPoints

          sysmem(i)%latMean = &

               sum(sysa(ts1)%sys(i)%lat(1:sysmem(i)%nPoints))/&

               sysmem(i)%nPoints

          !070512----------- Weight averages with Hm0 ---------------------

          temp1 = 0.

          temp2 = 0.

          DO iii = 1,sysmem(i)%nPoints

            temp1 = temp1 + &

                 (sysa(ts1)%sys(i)%hs(iii)**2)*sysa(ts1)%sys(i)%lon(iii)

            temp2 = temp2 + &

                 (sysa(ts1)%sys(i)%hs(iii)**2)*sysa(ts1)%sys(i)%lat(iii)

          END DO

          sysmem(i)%lonMean = temp1/&

               max(sum(sysa(ts1)%sys(i)%hs(1:sysmem(i)%nPoints)**2),&

               0.001)

          sysmem(i)%latMean = temp2/&

               max(sum(sysa(ts1)%sys(i)%hs(1:sysmem(i)%nPoints)**2),&

               0.001)

          !070512----------- Weight averages with Hm0 ---------------------

          sysmem(i)%tpMean = sysa(ts1)%sys(i)%tpMean

          sysmem(i)%dirMean = sysa(ts1)%sys(i)%dirMean

          sysmem(i)%updated = ts1

          sysmem(i)%length = 1

          !071012----------- Grid point indexing --------------------------

          DO iii = 1,sysmem(i)%nPoints

            sysmem(i)%indx(iii) = (sysa(ts1)%sys(i)%j(iii)-1)*maxi +&

                 sysa(ts1)%sys(i)%i(iii)

          END DO

          !071012----------- Grid point indexing --------------------------

        END DO

        maxgroup = maxsys(ts1)

        !         i = ts1

      END IF


      !******** Test output ***********************

      DO i = 1, maxgroup

        WRITE(20,*) 'sysMem(',i,')%grp =',sysmem(i)%grp

        WRITE(20,*) 'sysMem(',i,')%nPoints =',sysmem(i)%nPoints

        WRITE(20,*) 'sysMem(',i,')%lonMean =',sysmem(i)%lonMean

        WRITE(20,*) 'sysMem(',i,')%latMean =',sysmem(i)%latMean

        WRITE(20,*) 'sysMem(',i,')%tpMean =',sysmem(i)%tpMean

        WRITE(20,*) 'sysMem(',i,')%dirMean =',sysmem(i)%dirMean

        WRITE(20,*) 'sysMem(',i,')%updated =',sysmem(i)%updated

        WRITE(20,*) 'sysMem(',i,')%length =',sysmem(i)%length

      END DO

      !********************************************

    END IF


    !     Loop over all time levels to track systems in time

    WRITE(20,*) 'Number of time levels = ',SIZE(sysa)

    DO i = (ts1+1), SIZE(sysa)

      WRITE(20,*) 'TS = ',i


      IF (SIZE(sysa(i)%sys).GT.0) THEN

        !           *** Added: 02/29/12 *************************************

        !           Sort groups, so that larger systems get associated first

        ALLOCATE( sysordered(maxsys(i)) )

        ALLOCATE( indsorted(maxsys(i)) )

        CALL sort (real(sysa(i)%sys(1:maxsys(i))%nPoints), &

             maxsys(i),sysordered,indsorted,'D')

        sysa(i)%sys(1:maxsys(i)) = sysa(i)%sys(indsorted)

        IF (ALLOCATED(sysordered)) DEALLOCATE(sysordered)

        IF (ALLOCATED(indsorted)) DEALLOCATE(indsorted)

        !           *** Added: 02/29/12 *************************************


        !           Initialize groups                                             ! Optimize?

        sysa(i)%sys(:)%grp = 9999                                     ! Optimize?

        counter = 0

        leng = length(real(sysmem(:)%grp), &

             SIZE(sysmem(:)%grp),real(9999))

        ALLOCATE( alreadyused(leng+10) )                              !Make space for 10 new potential entries. Improve!!!

        WRITE(20,*) 'sysMem(1:leng)%grp =', &

             sysmem(1:leng)%grp

        ALLOCATE( allind(leng) )

        alreadyused(:) = 0

        allind(:) = sysmem(1:leng)%grp


        !071212-----GoF 2D-------------------------------

        ALLOCATE( ind(SIZE(allind)) )

        ind(:) = allind

        ALLOCATE( ind2(SIZE(ind)) )

        DO ii = 1, SIZE(ind)

          ind2(ii) = findfirst(real(allind),SIZE(allind), &

               REAL(ind(ii)))

        END DO

        !        Define 2D array for evaluating best fit for systems

        ALLOCATE( gof(maxsys(i),maxgroup) )

        ALLOCATE( gofminval(maxgroup) )

        ALLOCATE( gofminind(maxgroup) )

        ALLOCATE( tbsysmem(maxgroup) )

        ALLOCATE( deltadirsysmem(maxgroup) )

        ALLOCATE( deltapersysmem(maxgroup) )

        ALLOCATE( m1sysmem(maxgroup) )

        ALLOCATE( m2sysmem(maxgroup) )

        !071212-----GoF 2D-------------------------------

        DO j = 1, maxsys(i)

          npnts = sysa(i)%sys(j)%nPoints

          lonmean = sum(sysa(i)%sys(j)%lon(1:npnts))/npnts

          latmean = sum(sysa(i)%sys(j)%lat(1:npnts))/npnts

          !070512----------- Weight averages with Hm0 ---------------------

          temp1 = 0.

          temp2 = 0.

          DO iii = 1,npnts

            temp1 = temp1 + &

                 (sysa(i)%sys(j)%hs(iii)**2)*sysa(i)%sys(j)%lon(iii)

            temp2 = temp2 + &

                 (sysa(i)%sys(j)%hs(iii)**2)*sysa(i)%sys(j)%lat(iii)

          END DO

          lonmean=temp1/max(sum(sysa(i)%sys(j)%hs(1:npnts)**2),0.001)

          latmean=temp2/max(sum(sysa(i)%sys(j)%hs(1:npnts)**2),0.001)

          !070512----------- Weight averages with Hm0 ---------------------

          !071012----------- Grid point indexing --------------------------

          ALLOCATE(sysa(i)%sys(j)%indx(maxi*maxj))

          sysa(i)%sys(j)%indx = 9999

          DO iii = 1,sysa(i)%sys(j)%nPoints

            sysa(i)%sys(j)%indx(iii) = &

                 (sysa(i)%sys(j)%j(iii)-1)*maxi + &

                 sysa(i)%sys(j)%i(iii)

          END DO

          !071012----------- Grid point indexing --------------------------

          WRITE(20,*) 'System no. ',j,' of ',maxsys(i)

          WRITE(20,*) 'Size =', npnts

          WRITE(20,*) 'lonMean =', lonmean

          WRITE(20,*) 'latMean =', latmean

          WRITE(20,*) 'tpMean =', sysa(i)%sys(j)%tpMean

          WRITE(20,*) 'dirMean =', sysa(i)%sys(j)%dirMean

          sysa(i)%sys(j)%grp = 9999                                     !Now redundant?


          !           Compute deltas

          tbsysmem = sysmem(1:maxgroup)%tpMean

          WRITE(20,*) 'Tbsysmem(:) =       ', tbsysmem(:)

          !           Compute deltas the same way as for field combining - they should

          !           be of the same degree of strictness as the latter, otherwise

          !           the time combining will lose track!

          !3stddev            m1 = -3.645*Tb + 63.211

          !3stddev            m1 = MAX(m1,10.)

          !3stddev            m2 = -0.346*Tb + 3.686

          !3stddev            m2 = MAX(m2,0.6)

          !1stddev            m1 = -2.219*Tb + 35.734

          !1stddev            m1 = MAX(m1,5.)

          !1stddev            m2 = -0.226*Tb + 2.213

          !1stddev            m2 = MAX(m2,0.35)

          !071412            m1 = -5.071*Tb + 90.688

          !071412            m1 = MAX(m1,16.)

          !071412            m2 = -0.467*Tb + 5.161

          !071412            m2 = MAX(m2,1.0)

          !071412            deltaDir = (m1*1. + dirTimeKnob)*1.

          !071412            deltaPer = (m2*1. + tpTimeKnob)*1.

          DO ii = 1,SIZE(ind2)

            m1sysmem(ii) = max((-3.645*tbsysmem(ii)+63.211),10.)

            m2sysmem(ii) = max((-0.346*tbsysmem(ii)+3.686),0.6)

          END DO

          deltadirsysmem = m1sysmem(:)*1. + dirtimeknob

          deltapersysmem = m2sysmem(:)*1. + tptimeknob

          WRITE(20,*) 'deltaDirsysmem(:) = ',deltadirsysmem

          WRITE(20,*) 'deltaPersysmem(:) = ',deltapersysmem


          !                 Criterion 1: Mean period

          ALLOCATE( temp(SIZE(ind2)) )

          temp = abs( sysa(i)%sys(j)%tpMean - &

               sysmem(ind2(:))%tpMean )

          WRITE(20,*) 'tpMean list =', &

               sysmem(ind2(:))%tpMean

          WRITE(20,*) 'tpMinVal list =', temp

          tpminval = minval(temp)

          tpminind = findfirst(temp,SIZE(temp),tpminval)


          !                 Criterion 2: Mean direction

          ALLOCATE( dirs(SIZE(ind2)) )

          dirs(:)=abs( sysa(i)%sys(j)%dirMean - &

               sysmem(ind2(:))%dirMean )

          !                 Deal with wrap around

          DO idir = 1, SIZE(dirs)

            IF (dirs(idir).GE.180.) dirs(idir)=360-dirs(idir)

          END DO

          WRITE(20,*) 'dirMean list =', &

               sysmem(ind2(:))%dirMean

          WRITE(20,*) 'dirMinVal list =', dirs


          !                 Criterion 3: Size

          WRITE(20,*) 'Size list =', &

               sysmem(ind2(:))%nPoints


          !                 Criterion 4: Distance between systems

          ALLOCATE (mnlonlist(SIZE(ind2)))

          ALLOCATE (mnlatlist(SIZE(ind2)))

          ALLOCATE (mndist(SIZE(ind2)))

          DO ii = 1,SIZE(ind2)

            mnlonlist(ii) = sysmem(ind2(ii))%lonMean

            mnlatlist(ii) = sysmem(ind2(ii))%latMean

            mndist(ii) = sqrt((lonmean-mnlonlist(ii))**2 + &

                 (latmean-mnlatlist(ii))**2)

          END DO

          dmndiag = sqrt(lonext**2+latext**2)

          WRITE(20,*) 'Distance list =',mndist(:)

          WRITE(20,*) 'Domain diagonal =',dmndiag


          !                 Criterion 5: Overlap of systems

          ALLOCATE (olsize(SIZE(ind2)))

          DO ii = 1,SIZE(ind2)


            IF (sysmem(ind2(ii))%nPoints.GT.0) THEN

              CALL intersect(real(sysa(i)%sys(j)%indx(1:npnts)),npnts, &

                   REAL(sysMem(ind2(ii))%indx(1:sysMem(ind2(ii))%nPoints)),&

                   sysMem(ind2(ii))%nPoints,dummy1,olsize(ii),dummy2,dummy3)

            ELSE

              olsize(ii) = 0

            END IF

          END DO


          gof(j,1:SIZE(ind2)) = (temp/deltapersysmem(:))**2 + &

               (dirs/deltadirsysmem(:))**2 + &

               !                    (4*mndist(:)/dmndiag)**2

               ( (real(olsize(:)) - &

               REAL(sysMem(ind2(:))%nPoints) )/&

               (0.50*MAX(REAL(sysMem(ind2(:))%nPoints),0.001)) )**2

          ! Remove GoF entries which exceed predifined tolerances

          DO ii = 1,SIZE(ind2)

            WRITE(20,*) 'Testing: ii,olsize(ii),size,frac =',&

                 ii,olsize(ii),sysmem(ind2(ii))%nPoints,&

                 REAL(olsize(ii))/&

                 MAX(REAL(sysMem(ind2(ii))%nPoints),0.001)

            IF ( real(olsize(ii)).LT.&

                 0.50*real(sysmem(ind2(ii))%nPoints) ) THEN

              gof(j,ii) = 9999.

            END IF

            IF ( (temp(ii).GT.deltapersysmem(ii)).OR.&

                 (dirs(ii).GT.deltadirsysmem(ii)) ) THEN

              gof(j,ii) = 9999.

            END IF

          END DO

          WRITE(20,*) 'GOF(j,:) =',gof(j,:)


          IF (ALLOCATED(temp)) DEALLOCATE(temp)

          IF (ALLOCATED(dirs)) DEALLOCATE(dirs)

          IF (ALLOCATED(mnlonlist)) DEALLOCATE(mnlonlist)

          IF (ALLOCATED(mnlatlist)) DEALLOCATE(mnlatlist)

          IF (ALLOCATED(mndist)) DEALLOCATE(mndist)

          IF (ALLOCATED(olsize)) DEALLOCATE(olsize)


          !071212-----------GoF 2D-------------

        END DO

        IF (ALLOCATED(tbsysmem)) DEALLOCATE(tbsysmem)

        IF (ALLOCATED(deltadirsysmem)) DEALLOCATE(deltadirsysmem)

        IF (ALLOCATED(deltapersysmem)) DEALLOCATE(deltapersysmem)

        IF (ALLOCATED(m1sysmem)) DEALLOCATE(m1sysmem)

        IF (ALLOCATED(m2sysmem)) DEALLOCATE(m2sysmem)


        WRITE(20,*) 'GoF3:'

        DO jj = 1,maxsys(i)

          WRITE(20,*) gof(jj,:)

        END DO


        !        Find minima in GoF

        DO k = 1,maxgroup

          gofminval(k) = minval(gof(:,k))

          gofminind(k) = findfirst(gof(:,k),SIZE(gof,1),gofminval(k))

          IF (gofminval(k).EQ.9999) THEN

            gofminind(k) = 0

          END IF

        END DO


        IF (ALLOCATED(gof)) DEALLOCATE(gof)


        DO j = 1, maxsys(i)

          new = 0

          !                 Look up sysMem match for this current system. If no match

          !                 is found, the index value 0 is returned.

          tpminind = 0

          temp1 = 9999.

          DO jj = 1, SIZE(gofminind)

            IF (gofminind(jj).EQ.j) THEN

              IF (gofminval(jj).LT.temp1) THEN

                tpminind = jj

                temp1 = gofminval(jj)

              END IF

            END IF

          END DO

          dirminind = tpminind

          WRITE(20,*) 'System, GOFMinInd: ',j,tpminind


          IF (tpminind.NE.0) THEN

            !                    Success

            !071212-----------GoF 2D-------------


            counter = counter+1

            sysa(i)%sys(j)%grp = &

                 sysmem(ind2(dirminind))%grp

            alreadyused(counter) = sysa(i)%sys(j)%grp


            WRITE(20,*) 'Case 1: matched this ts (',i, &

                 ') sys ',sysa(i)%sys(j)%sysInd,' (tp=', &

                 sysa(i)%sys(j)%tpMean,' dir=', &

                 sysa(i)%sys(j)%dirMean,') with grp ', &

                 sysmem(ind2(dirminind))%grp

            WRITE(20,*) 'Added ',alreadyused(counter), &

                 ' in array *alreadyUsed*'

          ELSE

            new = 1

          END IF


          IF (new.EQ.1) THEN

            used = 0

            DO k = 1, maxgroup

              ok = 1

              WRITE(20,*) 'maxGroup,k,ok,used =', &

                   maxgroup,k,ok,used

              !                       Make sure it hasn't been used yet (at current time level)

              IF ((i.GT.2).AND. &

                   (.NOT.any(alreadyused(:).EQ.k))) THEN

                !                           Make sure it hasn't been used yet (at previous time level)

                DO l = 1, maxgroup

                  !                              If last update of system was more that *6* time steps

                  !                              ago, system can be released (TO CALIBRATE)

                  IF ( (sysmem(l)%grp.EQ.k).AND. &

                       ((i-sysmem(l)%updated).LT.6) ) ok = 0

                  WRITE(20,*) 'l, ok = ',l,ok

                END DO

                IF (ok.EQ.1) THEN

                  sysa(i)%sys(j)%grp = k

                  counter = counter+1;

                  alreadyused(counter) = k

                  used = 1

                  WRITE(20,*) 'k,used,counter =', &

                       k,used,counter

                  EXIT

                END IF

              END IF

            END DO

            IF (used.EQ.0) THEN

              maxgroup = maxgroup+1

              sysa(i)%sys(j)%grp = maxgroup

              !                       Increase sysMem by one slot

              sysmem(maxgroup)%grp = maxgroup

              counter = counter+1

              alreadyused(counter) = maxgroup

            END IF

            WRITE(20,*) 'counter,maxGroup,sysA(i)%sys(j)%grp =',&

                 counter,maxgroup,sysa(i)%sys(j)%grp

            WRITE(20,*) 'NO GRP MATCH case 2'

          END IF


        END DO

        IF (ALLOCATED(ind)) DEALLOCATE(ind)                  !071212 Shifted

        IF (ALLOCATED(ind2)) DEALLOCATE(ind2)                !071212 Shifted

        IF (ALLOCATED(gofminval)) DEALLOCATE(gofminval)

        IF (ALLOCATED(gofminind)) DEALLOCATE(gofminind)


        IF (ALLOCATED(alreadyused)) DEALLOCATE(alreadyused)

        IF (ALLOCATED(allind)) DEALLOCATE(allind)


        !        Update sysMem

        DO k = 1, maxgroup

          DO kk = 1, maxsys(i)

            IF (sysa(i)%sys(kk)%grp.EQ.sysmem(k)%grp) THEN

              sysmem(k)%nPoints = sysa(i)%sys(kk)%nPoints

              sysmem(k)%lonMean = &

                   sum(sysa(i)%sys(kk)%lon(1:sysmem(k)%nPoints))/&

                   sysmem(k)%nPoints

              sysmem(k)%latMean = &

                   sum(sysa(i)%sys(kk)%lat(1:sysmem(k)%nPoints))/&

                   sysmem(k)%nPoints

              !070512----------- Weight averages with Hm0 ---------------------

              temp1 = 0.

              temp2 = 0.

              DO iii = 1,sysmem(k)%nPoints

                temp1 = temp1 + &

                     (sysa(i)%sys(kk)%hs(iii)**2)*sysa(i)%sys(kk)%lon(iii)

                temp2 = temp2 + &

                     (sysa(i)%sys(kk)%hs(iii)**2)*sysa(i)%sys(kk)%lat(iii)

              END DO

              sysmem(k)%lonMean = temp1/&

                   max(sum(sysa(i)%sys(kk)%hs(1:sysmem(k)%nPoints)**2),&

                   0.001)

              sysmem(k)%latMean = temp2/&

                   max(sum(sysa(i)%sys(kk)%hs(1:sysmem(k)%nPoints)**2),&

                   0.001)

              !070512----------- Weight averages with Hm0 ---------------------

              sysmem(k)%tpMean = sysa(i)%sys(kk)%tpMean

              sysmem(k)%dirMean = sysa(i)%sys(kk)%dirMean

              !071012----------- Grid point indexing --------------------------

              sysmem(k)%indx(:) = 9999

              DO iii = 1,sysmem(k)%nPoints

                sysmem(k)%indx(iii) = &

                     (sysa(i)%sys(kk)%j(iii)-1)*maxi + &

                     sysa(i)%sys(kk)%i(iii)

              END DO

              !071012----------- Grid point indexing --------------------------

              sysmem(k)%updated = i

              sysmem(k)%length = sysmem(k)%length + 1

            END IF

          END DO

          !Test for expired groups

          IF ((i-sysmem(k)%updated).GE.6) THEN

            sysmem(k)%nPoints = 0

            sysmem(k)%lonMean = 9999.

            sysmem(k)%latMean = 9999.

            sysmem(k)%tpMean = 9999.

            sysmem(k)%dirMean = 9999.

            sysmem(k)%indx(:) = 9999

            sysmem(k)%updated = -9999

            sysmem(k)%length = 0

          END IF

        END DO

        !083012  !Filter out duplicates groups that can develop

        DO l = 1, maxgroup

          DO ll = (l+1), maxgroup


            deltadir = max((-3.645*sysmem(l)%tpMean+63.211),10.)*1.

            deltaper = max((-0.346*sysmem(l)%tpMean+3.686),0.6)*1.


            IF ( (abs(sysmem(l)%tpMean-sysmem(ll)%tpMean).LT.&

                 deltaper).AND. &

                 (abs(sysmem(l)%dirMean-sysmem(ll)%dirMean).LT.&

                 deltadir).AND. &

                 (sysmem(l)%updated.NE.sysmem(ll)%updated).AND. &

                 (sysmem(ll)%nPoints.NE.0) ) THEN

              !Find the more recent entry, and delete from group

              IF (sysmem(ll)%length.LT.sysmem(l)%length) THEN

                idup = ll

                WRITE(20,*) 'Deleting memgroup ',ll, &

                     '(updated',sysmem(ll)%updated,', length', &

                     sysmem(ll)%length,'), duplicate of memgroup', &

                     l,'(updated',sysmem(l)%updated,', length', &

                     sysmem(l)%length,'):'

              ELSE

                idup = l

                WRITE(20,*) 'Deleting memgroup ',l, &

                     '(updated',sysmem(l)%updated,', length', &

                     sysmem(l)%length,'), duplicate of memgroup', &

                     ll,'(updated',sysmem(ll)%updated,', length', &

                     sysmem(ll)%length,'):'

              END IF

              WRITE(20,*) 'deltaPer, diff Per:',deltaper,&

                   abs(sysmem(l)%tpMean-sysmem(ll)%tpMean)

              WRITE(20,*) 'deltaDir, diff Dir:',deltadir,&

                   abs(sysmem(l)%dirMean-sysmem(ll)%dirMean)

              sysmem(idup)%nPoints = 0

              sysmem(idup)%lonMean = 9999.

              sysmem(idup)%latMean = 9999.

              sysmem(idup)%tpMean = 9999.

              sysmem(idup)%dirMean = 9999.

              sysmem(idup)%indx(:) = 9999

              sysmem(idup)%updated = -9999

              sysmem(idup)%length = 0

            END IF

          END DO

        END DO

      ELSE

        WRITE(20,*) '*** No systems at this time level. ', &

             'No. systems =',SIZE(sysa(i)%sys)

        !Test for expired groups

        DO k = 1, maxgroup

          IF ((i-sysmem(k)%updated).GE.6) THEN

            sysmem(k)%nPoints = 0

            sysmem(k)%lonMean = 9999.

            sysmem(k)%latMean = 9999.

            sysmem(k)%tpMean = 9999.

            sysmem(k)%dirMean = 9999.

            sysmem(k)%indx(:) = 9999

            sysmem(k)%updated = -9999

            sysmem(k)%length = 0

          END IF

        END DO

      END IF

      !        ******** Test output ***********************

      DO k = 1, maxgroup

        WRITE(20,*) 'sysMem(',k,')%grp =',sysmem(k)%grp

        WRITE(20,*) 'sysMem(',k,')%nPoints =',sysmem(k)%nPoints

        WRITE(20,*) 'sysMem(',k,')%lonMean =',sysmem(k)%lonMean

        WRITE(20,*) 'sysMem(',k,')%latMean =',sysmem(k)%latMean

        WRITE(20,*) 'sysMem(',k,')%tpMean =',sysmem(k)%tpMean

        WRITE(20,*) 'sysMem(',k,')%dirMean =',sysmem(k)%dirMean

        WRITE(20,*) 'sysMem(',k,')%updated =',sysmem(k)%updated

        WRITE(20,*) 'sysMem(',k,')%length =',sysmem(k)%length

      END DO

      !        ********************************************

    END DO


    !     Write hotfile of wave groups

    OPEN(unit=27,file='sys_restart1.ww3',status='unknown')

    WRITE(27,'(A23,I10)') 'maxGroup             =',maxgroup

    DO k = 1, maxgroup

      WRITE(27,'(A8,I3,A12,I10)') 'sysMem( ',k, &

           ' )%grp     =',sysmem(k)%grp

      WRITE(27,'(A8,I3,A12,I10)') 'sysMem( ',k, &

           ' )%nPoints =',sysmem(k)%nPoints

      WRITE(27,'(A8,I3,A12,F10.4)') 'sysMem( ',k, &

           ' )%lonMean =',sysmem(k)%lonMean

      WRITE(27,'(A8,I3,A12,F10.4)') 'sysMem( ',k, &

           ' )%latMean =',sysmem(k)%latMean

      WRITE(27,'(A8,I3,A12,F10.3)') 'sysMem( ',k, &

           ' )%tpMean  =',sysmem(k)%tpMean

      WRITE(27,'(A8,I3,A12,F10.3)') 'sysMem( ',k, &

           ' )%dirMean =',sysmem(k)%dirMean

      WRITE(27,'(A8,I3,A12,I10)') 'sysMem( ',k, &

           ' )%updated =',sysmem(k)%updated

      WRITE(27,'(A8,I3,A12,I10)') 'sysMem( ',k, &

           ' )%length  =',sysmem(k)%length

      DO j = maxj,1,-1

        DO i = 1,maxi

          WRITE(27,'(I8)',advance='NO') sysmem(k)%indx((j-1)*maxi+i)

        END DO

        WRITE(27,'(A)',advance='YES') ''

      END DO

    END DO

    CLOSE(27)


2000 CONTINUE

    RETURN

  END SUBROUTINE timetrackingv2

  !/ End of timeTrackingV2 --------------------------------------------- /

  !/

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

  SUBROUTINE findway (way         ,horizStepCount,vertStepCount , &

       vertBorder  ,horizBorder   ,stepCount     )

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          4-Jan-2013 |

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

    !/

    !/    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 )

    !/

    !/    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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     From the direction (way) we were going before, find which direction we

    !     are going now and how many 'steps' we need to take

    !

    !  2. Method

    !

    !     -

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     way         Char     in/out Direction of spiral search

    !     vertBorder  Int      input

    !     horizBorder Int      input

    !     stepCount   Int      output Number of steps to go in the selected direction (way)

    !

    CHARACTER   :: way *1

    INTEGER     :: horizStepCount, vertStepCount, &

         vertBorder, horizBorder, stepCount


    INTENT (IN) vertborder, horizborder

    INTENT (OUT) stepcount

    INTENT (IN OUT) way

    !

    !     Local variables

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

    !     -

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !     -

    !

    !  5. Subroutines calling

    !

    !     spiralTrackV3

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See above

    !

    !  9. Switches :

    !

    !     None defined yet.

    !

    ! 10. Source code :

    !

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

    SELECT CASE (way)

    CASE ('R')

      way='D'

      vertstepcount=vertstepcount+1

      IF (horizborder.EQ.1) THEN

        horizstepcount=horizstepcount-1

      END IF

      stepcount=vertstepcount

    CASE ('D')

      way='L'

      horizstepcount=horizstepcount+1

      IF (vertborder.EQ.1) THEN

        vertstepcount=vertstepcount-1

      END IF

      stepcount=horizstepcount

    CASE ('L')

      way='U'

      vertstepcount=vertstepcount+1

      IF (horizborder.EQ.1) THEN

        horizstepcount=horizstepcount-1

      END IF

      stepcount=vertstepcount

    CASE ('U')

      way='R'

      horizstepcount=horizstepcount+1

      IF (vertborder.EQ.1) THEN

        vertstepcount=vertstepcount-1

      END IF

      stepcount=horizstepcount

    CASE DEFAULT

      WRITE(20,*) 'In spaTack:findWay should NOT go here!'

    END SELECT


    RETURN

  END SUBROUTINE findway

  !/ End of findWay ---------------------------------------------------- /

  !/

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

  SUBROUTINE findnext (i      ,j          ,maxI       ,maxJ     , &

       way    ,vertBorder ,horizBorder          )

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           H. L. Tolman            |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          4-Jan-2013 |

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

    !/

    !/    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 )

    !/

    !/    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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     Find next point on spatial search spiral

    !

    !  2. Method

    !

    !     -

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     i,j         Int      in/out  Current grid indices

    !     maxI, maxJ  Int      input   Maximum indices of wave field

    !     way         Char     input   Direction of spiral search

    !     vertBorder  Int      output  Flag indicating that vert domain edge has been hit

    !     horizBorder Int      output  Flag indicating that hor domain edge has been hit

    !

    CHARACTER   :: way

    INTEGER     :: i, j, maxI, maxJ, vertBorder, horizBorder


    INTENT (IN) maxi, maxj, way

    INTENT (IN OUT) i, j

    INTENT (OUT) vertborder, horizborder

    !

    !     Local variables

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

    !     -

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !     -

    !

    !  5. Subroutines calling

    !

    !     spiralTrackV3

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     -

    !

    !  9. Switches :

    !

    !     None defined yet.

    !

    ! 10. Source code :

    !

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

    vertborder=0

    horizborder=0

    SELECT CASE (way)

    CASE ('R')

      IF (i.LT.maxi) THEN

        i=i+1

      ELSE

        !                 Need to tell findWay that if we hit the border we don't

        !                 increment stepCount...

        horizborder=1

      END IF

    CASE ('D')

      IF (j.GT.1) THEN

        j=j-1

      ELSE

        vertborder=1

      END IF

    CASE ('L')

      IF (i.GT.1) THEN

        i=i-1

      ELSE

        horizborder=1

      END IF

    CASE ('U')

      IF (j.LT.maxj) THEN

        j=j+1

      ELSE

        vertborder=1

      END IF

    END SELECT


    RETURN

  END SUBROUTINE findnext

  !/ End of findNext --------------------------------------------------- /

  !/

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

  SUBROUTINE findsys (i        ,j        ,wsdat    ,maxSys    , &

       ngbrExt  ,maxI     ,maxJ     ,perKnob   , &

       dirKnob  ,hsKnob                           )

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          4-Jan-2013 |

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

    !/

    !/    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 )

    !/

    !/    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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     Find all wave systems that neighbour the grid point (i,j), and

    !     match these with the systems at (i,j).

    !

    !  2. Method

    !

    !     For the given point (i,j), find all wave systems at neighbouring grid

    !     points within the reach specified by ngbrExt.

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     i,j         Int      input   Current grid indices

    !     maxI, maxJ  Int      input   Maximum indices of wave field

    !     wsdat   Type(dat2d)  in/out  Input data structure to be spiral tracked

    !     maxSys      Int      in/out  Maximum number of systems identified

    !

    TYPE(dat2d) :: wsdat

    INTEGER     :: i, j, maxI, maxJ, ngbrExt, maxSys

    REAL        :: perKnob ,dirKnob, hsKnob


    INTENT (IN) i, j, maxi, maxj, ngbrext, perknob ,dirknob

    INTENT (IN OUT) wsdat, maxsys

    !

    !     Local variables

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

    !     tmpsys  TYPE(system)   Temporary instance of the wave system variable

    !     nngbr       Int        Number of neighbours found

    !

    TYPE(system), ALLOCATABLE :: tmpsys(:)

    TYPE(neighbr) :: ngbr(50)

    TYPE(mtchsys) :: match

    LOGICAL       :: found

    INTEGER       :: counter, ii, jj, nngbr, startCount, endCount, l,&

         nout, maxS, s, p, n, countAll, ind, minInd, &

         npart, pp, leng

    INTEGER       :: allFullSys(50)

    REAL, POINTER :: realarr(:)

    INTEGER, ALLOCATABLE :: allSys(:)

    REAL         :: hsAll(50),tpAll(50),dirAll(50),GOF(50)

    REAL         :: absDir,absPer,absHs,T,&

         deltaPer,deltaDir,deltaHs,temp

    REAL         :: dx, m1, m2

    REAL         :: GOFMinVal

    INTEGER      :: GOFMinInd

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      UNIQUE

    !      combinePartitionsV2

    !

    !  5. Subroutines calling

    !

    !     spiralTrackV3

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     -

    !

    !  9. Switches :

    !

    !     None defined yet.

    !

    ! 10. Source code :

    !

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

    NULLIFY(realarr)

    !      WRITE(20,*) 'findSys: i,j,maxSys =',i,j,maxSys


    !     First find the checked neighbor

    counter=1

    DO ii=(i-ngbrext), (i+ngbrext)

      DO jj=(j-ngbrext), (j+ngbrext)

        IF ( (ii.GT.0).AND.(jj.GT.0).AND. &

             (jj.LE.maxj).AND.(ii.LE.maxi) ) THEN

          IF ( wsdat%par(ii,jj)%checked.EQ.1 ) THEN

            ngbr(counter)%par = wsdat%par(ii,jj)                !Added the par field to maintain the data structure

            ngbr(counter)%i = ii

            ngbr(counter)%j = jj

            counter=counter+1

          END IF

        END IF

      END DO

    END DO

    !     New variable nngbr

    nngbr=counter-1


    IF (nngbr.GT.0) THEN

      allfullsys(:) = 0

      startcount=1

      l=1

      DO WHILE (l.LE.nngbr)

        leng = length(real(ngbr(l)%par%sys), &

             SIZE(ngbr(l)%par%sys),real(9999))

        endcount = startcount+leng-1

        allfullsys(startcount:endcount) = ngbr(l)%par%sys(1:leng)

        startcount=endcount+1

        l=l+1

      END DO


      IF (endcount.EQ.0) WRITE(20,*) '***1.Calling UNIQUE w. len=0!'

      CALL unique (real(allfullsys),endcount,realarr,nout)            !Can one do this?

      ALLOCATE(allsys(nout))

      allsys = int(realarr)                                           !Can one do this?

      IF (ASSOCIATED(realarr)) DEALLOCATE(realarr)

      maxs = maxval(allsys)


      IF (maxsys.LT.maxs) THEN

        maxsys=maxs

      END IF

      !         Initiate sys num

      ALLOCATE( tmpsys(SIZE(allsys)) )

      !         Clear the wsdat%par(i,j)%sys field, new values assigned below.

      !         System info temporarily stored in allSys

      wsdat%par(i,j)%sys(1:10) = 9999


      DO s=1, SIZE(allsys)

        hsall(:) = 0.

        tpall(:) = 0.

        dirall(:) = 0.

        !              wfAll(:) = 0.

        n=1

        countall=0

        DO WHILE (n.LE.nngbr)

          !                Calculate mean of common neighbor wave system

          !                for every neigbor wave system

          found = .false.

          DO ind = 1, SIZE(ngbr(n)%par%sys)                       !Optimize this?

            IF ( ngbr(n)%par%sys(ind).EQ.allsys(s) ) THEN        !Put sys under par to maintain structure

              found = .true.

              EXIT

            END IF

          END DO


          IF (found) THEN

            countall=countall+1

            hsall(countall)=ngbr(n)%par%hs(ind)

            tpall(countall)=ngbr(n)%par%tp(ind)

            dirall(countall)=ngbr(n)%par%dir(ind)

            !                     wfAll(countAll)=ngbr(n)%par%wf(ind)

          ELSE

            n=n+1

            cycle

          END IF

          n=n+1

        END DO

        tmpsys(s)%hsMean = sum(hsall(1:countall))/countall

        tmpsys(s)%tpMean = sum(tpall(1:countall))/countall

        tmpsys(s)%dirMean = &

             mean_anglev2(dirall(1:countall),countall)

        !              tmpsys(s)%wfMean = SUM(wfAll(1:countAll))/countAll

      END DO


      !         Find the partition at current (i,j) point that matches previously

      !         identified wave systems if any...

      wsdat%par(i,j)%ngbrSys(1:SIZE(allsys)) = allsys


      npart = length(real(wsdat%par(i,j)%ipart), &

           SIZE(wsdat%par(i,j)%ipart),real(0))

      DO p = 1, npart

        IF ( (wsdat%par(i,j)%hs(p).LT.hsknob).OR. &

             (wsdat%par(i,j)%tp(p).EQ.0.) ) THEN

          wsdat%par(i,j)%sys(p)=-1

          cycle

        END IF


        ind=0                                                       !Replaced 'index' by 'ind'

        match%sysVal(:) = 9999

        match%tpVal(:) = 9999.

        match%dirVal(:) = 9999.

        !              match%wfVal(:) = 9999.

        !             Cycle through the neighbouring systems identified above

        DO s=1,SIZE(allsys)

          abshs = abs(wsdat%par(i,j)%hs(p)-tmpsys(s)%hsMean)

          absper = abs(wsdat%par(i,j)%tp(p)-tmpsys(s)%tpMean)

          absdir = abs(wsdat%par(i,j)%dir(p)-tmpsys(s)%dirMean)

          !                  absWf = ABS(wsdat%par(i,j)%wf(p)-tmpsys(s)%wfMean)

          IF (absdir.GT.180) THEN

            absdir = 360 - absdir

            IF (absdir.LT.0) THEN

              WRITE(20,*) '*** WARNING: absDir negative!'

              WRITE(20,*) 'wsdat%par(i,j)%dir(p) =', &

                   wsdat%par(i,j)%dir(p)

              WRITE(20,*) 'tmpsys(s)%dirMean) =', &

                   tmpsys(s)%dirMean

            END IF

          END IF

          !                 Calculate delta dir and freq as a function of the partition

          !                 dir and freq

          t = tmpsys(s)%tpMean

          dx = 0.5*( (wsdat%lon(2,1)-wsdat%lon(1,1)) + &

               (wsdat%lat(1,2)-wsdat%lat(1,1)) )

          m1 = -3.645*t + 63.211

          m1 = max(m1,10.)

          m2 = -0.346*t + 3.686

          m2 = max(m2,0.6)

          !1stddev                  m1 = -2.219*T + 35.734

          !1stddev                  m1 = MAX(m1,5.)

          !1stddev                  m2 = -0.226*T + 2.213

          !1stddev                  m2 = MAX(m2,0.35)

          !5stddev                  m1 = -5.071*T + 90.688

          !5stddev                  m1 = MAX(m1,16.)

          !5stddev                  m2 = -0.467*T + 5.161

          !5stddev                  m2 = MAX(m2,1.0)

          deltadir = m1*dx + dirknob

          deltaper = m2*dx + perknob

          deltahs = 0.25*tmpsys(s)%hsMean

          IF ((absper.LT.deltaper).AND.(absdir.LT.deltadir)) THEN

            ind=ind+1

            match%sysVal(ind) = allsys(s)

            match%tpVal(ind) = absper

            match%dirVal(ind) = absdir

            match%hsVal(ind) = abshs

            !                      match%wfVal(ind) = absWf

          END IF

        END DO


        IF (ind.GT.0) THEN

          IF (ind.EQ.1) THEN

            wsdat%par(i,j)%sys(p) = match%sysVal(1)

          ELSE

            !                    Take the closest match, using GoF function

            gof(:) = 9999.

            gof(1:ind) = (match%tpVal(1:ind)/deltaper)**2 + &

                 (match%dirVal(1:ind)/deltadir)**2 + &

                 (match%hsVal(1:ind)/deltahs)**2

            gofminval = minval(gof(1:ind))

            gofminind = findfirst(gof(1:ind),ind,gofminval)

            wsdat%par(i,j)%sys(p) = match%sysVal(gofminind)         !The index of the system is swapped - the remaining info stays the same!

          END IF

        END IF

      END DO

    END IF


    !     Now check if 2 partitions have been associated to the same wave system, if

    !     so combine them

    npart = length(real(wsdat%par(i,j)%ipart), &

         SIZE(wsdat%par(i,j)%ipart),real(0))

    DO p = 1, (npart-1)                                                 !Could probably be optimized!

      DO pp = (p+1), npart

        IF (wsdat%par(i,j)%sys(p).EQ.wsdat%par(i,j)%sys(pp)) THEN

          !              There is at least one duplicate, so combine systems

          CALL combinepartitionsv2(wsdat%par(i,j))

        END IF

      END DO

    END DO

    !     Now that we have associated any possible partition to an existing

    !     wave system, we check if any wave system is free. If so give it a

    !     new wave system number

    npart = length(real(wsdat%par(i,j)%ipart), &

         SIZE(wsdat%par(i,j)%ipart),real(0))


    DO p = 1, npart

      IF (wsdat%par(i,j)%sys(p).EQ.9999) THEN

        maxsys = maxsys + 1

        wsdat%par(i,j)%sys(p) = maxsys

      END IF

    END DO

    wsdat%par(i,j)%checked=1


    IF (ALLOCATED(allsys)) DEALLOCATE(allsys)

    IF (ALLOCATED(tmpsys)) DEALLOCATE(tmpsys)


    RETURN

  END SUBROUTINE findsys

  !/ End of findSys ---------------------------------------------------- /

  !/

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

  SUBROUTINE combinewavesystems (wsdat   ,maxSys  ,maxPts   , &

       maxI    ,maxJ    ,perKnob  , &

       dirKnob ,hsKnob  ,combine  , &

       sys   )

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          4-Jan-2013 |

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

    !/

    !/    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 )

    !/

    !/    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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     Combine wave systems. Then remove small and low-energy systems from set,

    !     based on the parameters maxPts and maxHgt.

    !

    !  2. Method

    !

    !     -

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     wsdat   Type(dat2d)  output  Combined wave system data structure

    !     sys     Type(system) output  Final set of tracked systems, for one time level

    !     maxI, maxJ  Int      input   Maximum indices of wave field

    !     maxSys      Int      input   Maximum number of systems identified

    !     maxPts      Int      input   Number of points req for valid system

    !     hsKnob      Real     input   Parameter for identifying valid system

    !     combine     Int      input   Toggle: 1=combine systems; 0=do not combine


    TYPE(dat2d) :: wsdat

    TYPE(system), POINTER :: sys(:), systemp(:)

    INTEGER      :: maxSys, maxPts, maxI, maxJ, combine

    REAL         :: perKnob ,dirKnob, hsKnob


    INTENT (IN) maxpts, maxi, maxj, hsknob, combine

    INTENT (IN OUT) wsdat, maxsys                                       !In the Matlab code maxSys is only input ???

    !      INTENT (OUT) sys

    !

    !     Local variables

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

    !     nSys       Int   Number of wave systems (for checking iterative combining loop)

    !

    LOGICAL   :: found

    INTEGER, ALLOCATABLE :: sysOut(:)

    INTEGER, ALLOCATABLE :: actSysInd(:)

    INTEGER   :: iter, ok, nSys, mS, s, so, ss, ind, leng, &

         iw, jw, iloop

    INTEGER   :: actSys

    REAL      :: dev, hsCmp, maxHgt, temp(5)

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      printFinalSys

    !      combineSys

    !

    !  5. Subroutines calling

    !

    !     spiralTrackV3

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     -

    !

    !  9. Switches :

    !

    !     None defined yet.

    !

    ! 10. Source code :

    !

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


    !012912      WRITE(20,*) 'maxSys,maxPts,maxI,maxJ,hsKnob,combine =', &

    !012912                  maxSys,maxPts,maxI,maxJ,hsKnob,combine


    !     Set up initial index array of active systems

    IF (.NOT.ALLOCATED(actsysind)) ALLOCATE( actsysind(maxsys) )

    actsysind(1:maxsys) = (/ (ind, ind = 1, maxsys) /)

    !opt      WRITE(20,*) 'actSysInd =',actSysInd


    IF (combine.EQ.1) THEN

      !         Combine wave systems

      WRITE(20,*) 'Calling printFinalSys...'

      CALL printfinalsys (wsdat,maxsys,actsysind,maxi,maxj,1,sys)

      iter=0

      ok=0

      !         Keep on combining wave systems until all possible combining

      !         has been carried out (based on the combining criteria)

      DO WHILE (ok.EQ.0)

        iter = iter+1

        !             No of systems before combining

        IF (ALLOCATED(actsysind)) THEN

          nsys = SIZE(actsysind)

        ELSE

          nsys = maxsys

        END IF

        WRITE(20,'(A,A,I3,A,I5,A)') 'Calling combineSys for ', &

             'iteration',iter,' (maxSys =',nsys,').'


        !opt              WRITE(20,*) 'SIZE(sys)=',SIZE(sys)

        CALL combinesys (wsdat,sys,maxsys,maxi,maxj, &

             actsysind,perknob,dirknob)

        !             No of systems after combining

        !opt              WRITE(20,*) 'maxSys,nSys,SIZE(actSysInd) =', &

        !opt                           maxSys,nSys,SIZE(actSysInd)

        !              IF (maxSys.EQ.nSys) ok = 1

        IF (SIZE(actsysind).EQ.nsys) ok = 1

      END DO

    ELSE

      !         Do not combine wave systems

      CALL printfinalsys (wsdat,maxsys,actsysind,maxi,maxj,3,sys)

    END IF


    !     Remove small and low-energy systems from set, based on

    !     the parameters maxPts and maxHgt.

    !      ALLOCATE( sysOut(maxSys) )

    !      sysOut = sys(1:maxSys)%sysInd

    !      mS = maxSys

    ms = SIZE(actsysind)

    ss = 1

    WRITE(20,*) 'Filtering the set of',ms,'systems on size and mag.'


    DO so = 1, ms

      s = actsysind(so)

      !opt      NOTE: if we deallocate the individual records without

      !opt      compressing sys, then s and sysInd will remain the same

      ss = s


      leng = length(sys(ss)%hs,SIZE(sys(ss)%hs),9999.)

      dev  = std(sys(ss)%hs(1:leng),leng)

      hscmp = sys(ss)%hsMean + 2.*dev

      maxhgt = hsknob


      IF ( (hscmp.LT.maxhgt).OR.(sys(ss)%nPoints.LT.maxpts) ) THEN

        !             Remove system, and shift up indices to fill the gap

        DO ind = 1, maxsys

          !                Find index to remove

          IF (ind.EQ.ss) THEN

            !                   Shift up entries, deleting the duplicate partition

            !                   REPLACE WITH CSHIFT(ARRAY, SHIFT, dim)?

            !                    IF (ind.LT.maxSys) &

            !                       sys( ind:(maxSys-1) ) = sys( (ind+1):maxSys )

            IF (ind.LE.maxsys) THEN

              !                      Since we use pointers, we have to copy each index and

              !                      field individually. Otherwise memory corruption occurs.

              DO iloop = ind,ind

                sys(iloop)%sysInd = 9999

                sys(iloop)%nPoints = 0

                sys(iloop)%grp = 9999

                DEALLOCATE( sys(iloop)%hs )

                DEALLOCATE( sys(iloop)%tp )

                DEALLOCATE( sys(iloop)%dir )

                DEALLOCATE( sys(iloop)%dspr )

                !                          DEALLOCATE( sys(iloop)%wf )

                DEALLOCATE( sys(iloop)%i )

                DEALLOCATE( sys(iloop)%j )

                DEALLOCATE( sys(iloop)%lat )

                DEALLOCATE( sys(iloop)%lon )

                !                          DEALLOCATE( sys(iloop)%hsMean )

                !                          DEALLOCATE( sys(iloop)%tpMean )

                !                          DEALLOCATE( sys(iloop)%dirMean )

                !                          DEALLOCATE( sys(iloop)%ngbr )

              END DO

            END IF

          END IF

        END DO


        !             Update wsdat as well

        DO iw = 1, maxi

          DO jw = 1, maxj

            leng = length(real(wsdat%par(iw,jw)%sys), &

                 SIZE(wsdat%par(iw,jw)%sys),real(9999))

            ind = 1

            found = .false.

            !                    Identify system index (there are no duplicate

            !                    systems at this point.

            DO WHILE (ind.LE.leng)

              IF ( wsdat%par(iw,jw)%sys(ind).EQ.s ) THEN

                found = .true.

                EXIT

              END IF

              ind = ind + 1

            END DO

            IF (found) THEN

              !                       Blank out used record

              wsdat%par(iw,jw)%sys(ind) = 9999

              wsdat%par(iw,jw)%ipart(ind) = 9999

            END IF

          END DO

        END DO

      END IF

    END DO


    !     Compile array index of active systems in sys

    actsys = 0

    DO so = 1,maxsys

      IF (sys(so)%nPoints>0) actsys = actsys + 1

    END DO

    IF (ALLOCATED(actsysind)) DEALLOCATE(actsysind)

    ALLOCATE( actsysind(actsys) )

    actsys = 0

    DO so = 1,maxsys

      IF (sys(so)%nPoints>0) THEN

        actsys = actsys + 1

        actsysind(actsys) = sys(so)%sysInd

      END IF

    END DO


    !opt      WRITE(20,*) 'actSysInd =',actSysInd

    DO so = 1,SIZE(actsysind)

      s = actsysind(so)

      !opt         WRITE(20,*) 'sys(',s,')%sysInd =',sys(s)%sysInd

    END DO


    CALL printfinalsys (wsdat,maxsys,actsysind,maxi,maxj,1,sys)

    !opt      WRITE(20,*) 'actSysInd =',actSysInd

    !opt      DO so = 1,maxSys

    !opt         WRITE(20,*) 'sys(',so,')%sysInd =',sys(so)%sysInd, &

    !opt                     ',  sys(',so,')%nPoints =',sys(so)%nPoints

    !opt      END DO


    RETURN

  END SUBROUTINE combinewavesystems

  !/ End of combineWaveSystems ----------------------------------------- /

  !/

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

  SUBROUTINE printfinalsys (wsdat    ,maxSys    ,actSysInd , &

       maxI     ,maxJ      ,flag      ,sys    )

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          4-Jan-2013 |

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

    !/

    !/    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 )

    !/

    !/    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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     Output (print) the final output systems for this time step.

    !

    !  2. Method

    !

    !     -

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     wsdat   Type(dat2d)  input   Combined data structure

    !     maxI, maxJ  Int      input   Maximum indices of wave field

    !     maxSys      Int      input   Maximum number of systems identified

    !     flag        Int      input   Flag for printing system

    !     sys     Type(system) output  Final set of tracked systems, for one time level

    !

    TYPE(dat2d)  :: wsdat

    TYPE(system), POINTER :: sys(:)

    INTEGER      :: maxSys, maxI, maxJ, flag

    INTEGER, ALLOCATABLE :: actSysInd(:)


    INTENT (IN) wsdat, actsysind, maxi, maxj, flag

    INTENT (OUT) maxsys

    !      INTENT (IN OUT) sys

    !

    !     Local variables

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

    !     ic        Int   Counter for wave systems

    !

    INTEGER   :: ic, nGuys, startInd, endInd, i, j, ind, leng, leng2

    INTEGER   :: UNISIZE, DIFSIZE

    REAL, ALLOCATABLE :: sysOrdered(:)

    REAL, POINTER :: UNIARR(:), DIFARR(:)

    INTEGER, ALLOCATABLE :: ngbrSysAll(:), sysSortedInd(:)

    REAL      :: TEMP(2), TEMP1, TEMP2

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      UNIQUE

    !      SETDIFF

    !      SORT

    !

    !  5. Subroutines calling

    !

    !     combineWaveSystems

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     -

    !

    !  9. Switches :

    !

    !     None defined yet.

    !

    ! 10. Source code :

    !

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


    !     Initialize sys structure

    IF (flag.NE.2) THEN

      !         Allocate data structure with the final wave systems

      WRITE(20,*) 'In printFinalSys...'

      maxsys = SIZE(actsysind)

      NULLIFY(sys)

      ALLOCATE( sys(maxsys) )

      WRITE(20,*) 'Allocated sys okay, SIZE(sys) =',SIZE(sys)


      ALLOCATE( ngbrsysall(50*maxi*maxj) )                            !Large enough?

      DO ic = 1, maxsys

        NULLIFY( sys(ic)%hs )

        NULLIFY( sys(ic)%tp )

        NULLIFY( sys(ic)%dir )

        NULLIFY( sys(ic)%dspr )

        !              NULLIFY( sys(ic)%wf )

        NULLIFY( sys(ic)%i )

        NULLIFY( sys(ic)%j )

        NULLIFY( sys(ic)%lat )

        NULLIFY( sys(ic)%lon )

        ALLOCATE( sys(ic)%hs(maxi*maxj) )

        ALLOCATE( sys(ic)%tp(maxi*maxj) )

        ALLOCATE( sys(ic)%dir(maxi*maxj) )

        ALLOCATE( sys(ic)%dspr(maxi*maxj) )

        !              ALLOCATE( sys(ic)%wf(maxI*maxJ) )

        ALLOCATE( sys(ic)%i(maxi*maxj) )

        ALLOCATE( sys(ic)%j(maxi*maxj) )

        ALLOCATE( sys(ic)%lat(maxi*maxj) )

        ALLOCATE( sys(ic)%lon(maxi*maxj) )

        sys(ic)%hs(:) = 9999.                                       !Optimize this further?

        sys(ic)%tp(:) = 9999.

        sys(ic)%dir(:) = 9999.

        sys(ic)%dspr(:) = 9999.

        !              sys(ic)%wf(:) = 9999.

        sys(ic)%i(:) = 9999

        sys(ic)%j(:) = 9999

        sys(ic)%lat(:) = 9999.

        sys(ic)%lon(:) = 9999.

        sys(ic)%sysInd = 9999

        sys(ic)%hsMean = 9999.

        sys(ic)%tpMean = 9999.

        sys(ic)%dirMean = 9999.

        sys(ic)%nPoints = 0

        sys(ic)%ngbr(:) = 9999

        sys(ic)%grp = 9999

        ngbrsysall(:) = 0

        startind=1

        nguys=0


        DO i = 1, maxi

          DO j = 1, maxj

            !                     ind=wsdat.par(i,j).sys==ic;

            DO ind = 1, SIZE(wsdat%par(i,j)%sys)                !40.81 !Optimize this?

              IF (wsdat%par(i,j)%sys(ind).EQ.actsysind(ic)) &

                   THEN

                nguys=nguys+1

                sys(ic)%hs(nguys)=wsdat%par(i,j)%hs(ind)

                sys(ic)%tp(nguys)=wsdat%par(i,j)%tp(ind)

                sys(ic)%dir(nguys)=wsdat%par(i,j)%dir(ind)

                sys(ic)%dspr(nguys)=wsdat%par(i,j)%dspr(ind)

                !                           sys(ic)%wf(nGuys)=wsdat%par(i,j)%wf(ind)

                sys(ic)%i(nguys)=i

                sys(ic)%j(nguys)=j

                sys(ic)%lat(nguys)=wsdat%lat(i,j)

                sys(ic)%lon(nguys)=wsdat%lon(i,j)

                leng = length(real(wsdat%par(i,j)%ngbrSys), &

                     SIZE(wsdat%par(i,j)%ngbrSys),real(9999))

                endind = startind + leng-1

                ngbrsysall(startind:endind) = &

                     wsdat%par(i,j)%ngbrSys(1:leng)

                startind=endind+1

              END IF

            END DO

          END DO

        END DO


        !              if ~isempty(sys)

        IF (nguys.GT.0) THEN

          sys(ic)%sysInd=ic

          sys(ic)%hsMean = sum(sys(ic)%hs(1:nguys))/nguys

          sys(ic)%tpMean = sum(sys(ic)%tp(1:nguys))/nguys

          !                 sys(ic)%dirMean=mean_angle_single(sys(ic).dir)          40.81 Replaced with two-argument mean_angleV2

          sys(ic)%dirMean = &

               mean_anglev2(sys(ic)%dir(1:nguys),nguys)

          !070512----------- Weight averages with Hm0 ---------------------

          temp1 = 0.

          temp2 = 0.

          DO i = 1,nguys

            temp1 = temp1 + (sys(ic)%hs(i)**2)*sys(ic)%hs(i)

            temp2 = temp2 + (sys(ic)%hs(i)**2)*sys(ic)%tp(i)

          END DO

          sys(ic)%hsMean = &

               temp1/max(sum(sys(ic)%hs(1:nguys)**2),0.001)

          sys(ic)%tpMean = &

               temp2/max(sum(sys(ic)%hs(1:nguys)**2),0.001)

          sys(ic)%dirMean = mean_anglev3(sys(ic)%dir(1:nguys), &

               sys(ic)%hs(1:nguys),nguys)

          !070512----------- Weight averages with Hm0 ---------------------

          sys(ic)%nPoints = nguys

          IF (endind.GT.0) THEN

            CALL unique(real(ngbrsysall(1:endind)),endind, &

                 uniarr,unisize)

            temp = (/real(sys(ic)%sysInd),real(sys(ic)%sysInd)/)

            CALL setdiff(real(uniarr),unisize, &

                 temp,2,difarr,difsize)

            difsize = min(difsize,SIZE(sys(ic)%ngbr))

            sys(ic)%ngbr(1:difsize) = nint(difarr(1:difsize))

            IF (ASSOCIATED(uniarr)) DEALLOCATE(uniarr)

            IF (ASSOCIATED(difarr)) DEALLOCATE(difarr)

          END IF

        ELSE

          cycle

        END IF

      END DO

      IF (ALLOCATED(ngbrsysall)) DEALLOCATE(ngbrsysall)

    END IF


    !     Print the sorted field to the screen

    leng = length(real(sys(:)%nPoints), &

         SIZE(sys(:)%nPoints),real(9999))

    ALLOCATE( sysordered(leng) )

    ALLOCATE( syssortedind(leng) )

    CALL sort (real(sys(:)%nPoints),leng, &

         sysordered,syssortedind,'D')

    leng = length(real(sysordered), &

         SIZE(sysordered),real(0))


    DO ic = 1, leng

      leng2 = length(real(sys(syssortedind(ic))%ngbr), &

           SIZE(sys(syssortedind(ic))%ngbr),real(9999))

    END DO

    IF (ALLOCATED(sysordered)) DEALLOCATE(sysordered)

    IF (ALLOCATED(syssortedind)) DEALLOCATE(syssortedind)


    RETURN

  END SUBROUTINE printfinalsys

  !/ End of printFinalSys ---------------------------------------------- /

  !/

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

  SUBROUTINE combinesys (wsdat    ,sys      ,maxSys   ,maxI     , &

       maxJ     ,actSysInd,perKnob  ,dirKnob  )

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          4-Jan-2013 |

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

    !/

    !/    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 )

    !/

    !/    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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     Combine wave systems

    !

    !  2. Method

    !

    !     -

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     wsdat   Type(dat2d)  input   Combined data structure

    !     maxI, maxJ  Int      input   Maximum indices of wave field

    !     sys     Type(system) output  Final set of tracked systems, for one time level

    !     maxSys      Int      input  Number of systems

    !     dirKnob     Real     input  Parameter in direction for combining fields in space

    !     perKnob     Real     input  Parameter in period for combining fields in space

    !

    TYPE(dat2d)  :: wsdat                                               !40.PAR

    TYPE(system), POINTER :: sys(:)                                     !40.PAR

    INTEGER      :: maxSys, maxI, maxJ                                  !40.PAR

    INTEGER, ALLOCATABLE :: actSysInd(:)

    REAL         :: perKnob ,dirKnob

    REAL         :: dx, m1, m2


    INTENT (IN) maxi, maxj, perknob, dirknob                           !40.PAR

    !      INTENT (IN OUT) wsdat, sys, maxSys                                  !40.PAR

    !

    !     Local variables

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

    !     ngbIndex    Int   Arr   Array of neighbours

    !

    INTEGER, ALLOCATABLE :: sysSortedInd(:), sysOut(:)

    INTEGER, POINTER :: indSys1(:), indSys2(:)

    REAL, ALLOCATABLE    :: sysOrdered(:), rounded(:)

    REAL, POINTER    :: uniarr(:), difarr(:), allngbr(:)

    INTEGER   :: leng, leng2, s, ss, so, ngb, lsys, lsys2, hh, i, j, &

         ii, jj, ind, ind2, nn, nbr, icEnd,ic,iii,iloop

    INTEGER   :: myngbr, indMatch, matchSys, keep, replacedInd, &

         hhForIndMatch, lMatch, tot, outsize

    INTEGER   :: ngbIndex(10000), keepInd(maxI*maxJ), oneLess(1000)     !Array large enough?

    !      REAL      :: Tb,deltaPerB,deltaDirB,absDir,absPer,absHs,absWf

    REAL      :: Tb,deltaPerB,deltaDirB,deltaHsB,absDir,absPer,absHs

    LOGICAL   :: file_exists

    INTEGER   :: MASK(maxI,maxJ)

    REAL      :: lonmean, latmean, DIST

    !061512 -----------------------------------------------

    LOGICAL   :: ZIPMATCH

    INTEGER   :: counter, count2, izp, izp2, in, jn, icnt, ngbrExt

    REAL      :: T, ngb_tp, ngb_dir

    REAL      :: ngbmatch(maxI*maxJ)

    TYPE(neighbr) :: ngbr(50)

    !061512 -----------------------------------------------

    REAL      :: TEMP1, TEMP2

    INTEGER   :: actSys


    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      SORT

    !      findIJV4

    !      UNIQUE

    !      combinePartitionsV2

    !      UNION

    !      SETDIFF

    !

    !  5. Subroutines calling

    !

    !     combineWaveSystems

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     -

    !

    !  9. Switches :

    !

    !     None defined yet.

    !

    ! 10. Source code :

    !

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

    !     Initialize pointer (first use)

    NULLIFY(indsys1)

    NULLIFY(indsys2)

    !     Flag to combine systems on a point-by-point basis along boundary,

    !     instead of using mean values.

    zipmatch = .false.

    ngbrext = 1

    !     Combine systems on the basis of tpMean

    ALLOCATE( sysordered(maxsys) )

    ALLOCATE( syssortedind(maxsys) )

    ALLOCATE( sysout(maxsys) )

    ALLOCATE( rounded(maxsys) )

    !     Sort in descending Tp: the following improves the iterative combining in

    !     the special case that the wave period is constant over the domain, but

    !     tpMean is not because of truncation errors at very high decimals.

    rounded = real(int(sys(1:maxsys)%tpMean*1.e4))*1.e-4

    CALL sort(rounded,maxsys,sysordered,syssortedind,'D')

    sysout=sys(syssortedind)%sysInd

    IF (ALLOCATED(rounded)) DEALLOCATE(rounded)


    !051612 --- Land mask addition

    mask(:,:) = 0

    INQUIRE(file="sys_mask.ww3", exist=file_exists)

    IF (file_exists) THEN

      WRITE(20,*) '*** Using land mask'

      OPEN(unit=13,file='sys_mask.ww3',status='old')

      DO j = maxj,1,-1

        READ(13,*) (mask(i,j), i=1,maxi)

      END DO

      CLOSE(13)

    END IF

    !051612 --- Land mask addition


    !opt      WRITE(20,*) 'SIZE(sysOut)=',SIZE(sysOut)

    DO so = 1, SIZE(sysout)

      !          WRITE(20,*) 'so =',so

      s = sysout(so)

      ss = findfirst(real(sys(:)%sysInd),SIZE(sys(:)%sysInd), &

           REAL(s))

      !opt          WRITE(20,*) 's,ss=',s,ss

      ngbindex(:) = 0

      ii = 1

      leng = length(real(sys(ss)%ngbr),SIZE(sys(ss)%ngbr), &

           REAL(9999))

      !         Identify the indices of all the systems that neighbour the current system s,

      !         store in ngbIndex(:)

      DO ngb = 1, leng

        IF ( sys(ss)%ngbr(ngb).NE.s ) THEN

          myngbr = 1

          DO WHILE (myngbr.LE.SIZE(sysout))

            IF (sys(myngbr)%sysInd.EQ.sys(ss)%ngbr(ngb)) THEN

              ngbindex(ii) = myngbr

              ii = ii+1

              IF (ii.GT.1000) &

                   WRITE(20,*) '*** WARNING: ngbIndex(:) exceeded!'

            END IF

            myngbr = myngbr+1

          END DO

        END IF

      END DO

      ii = ii-1

      !opt          WRITE(20,*) so,'.   sys =',s,', Tp =',sys(s)%tpMean, &

      !opt                      ', size=',sys(s)%nPoints,', #neighbours=',ii


      IF ( ii.GT.0 ) THEN

        DO ngb = 1, ii

          !                 We first need to find the (i,j) points that are either common

          !                 to both these systems, or at the boundary of the two systems. Here

          !                 sys 1 will carry the 'ss' index and sys 2 the ngbIndex(ngb) index.

          CALL findijv4 (sys(ss),sys(ngbindex(ngb)), &

               maxi,maxj,indsys1,indsys2)

          IF ((SIZE(indsys1)>10).AND.(SIZE(indsys2)>10).AND. &

               (sys(ss)%nPoints.GT.sys(ngbindex(ngb))%nPoints)) &

               THEN

            lsys = SIZE(indsys1)

            lsys2 = SIZE(indsys2)


            !061512---------------Add zipper compare

            IF (zipmatch) THEN

              !                       Omit small systems to save time

              IF ((sys(ss)%nPoints.LT.5).OR. &

                   (sys(ngbindex(ngb))%nPoints.LT.5)) THEN

                cycle

              END IF

              dx=0.5*((wsdat%lon(2,1)-wsdat%lon(1,1)) + &

                   (wsdat%lat(1,2)-wsdat%lat(1,1)))

              ngbmatch(:)=0.

              DO izp = 1,lsys

                !                          Find neighbors of this point

                counter=0

                DO in=(sys(ss)%i(indsys1(izp))-ngbrext), &

                     (sys(ss)%i(indsys1(izp))+ngbrext)

                  DO jn=(sys(ss)%j(indsys1(izp))-ngbrext), &

                       (sys(ss)%j(indsys1(izp))+ngbrext)

                    counter=counter+1

                    ngbr(counter)%i = in

                    ngbr(counter)%j = jn

                  END DO

                END DO

                !                          Find these points in neighboring system

                ngb_tp = 0.

                ngb_dir = 0.

                count2 = 0

                DO izp2 = 1,lsys2

                  DO icnt = 1,counter

                    IF ((sys(ngbindex(ngb))%i(indsys2(izp2)) &

                         .EQ.ngbr(icnt)%i).AND. &

                         (sys(ngbindex(ngb))%j(indsys2(izp2)) &

                         .EQ.ngbr(icnt)%j)) THEN

                      count2 = count2+1

                      ngb_tp = ngb_tp + &

                           sys(ngbindex(ngb))%tp(indsys2(izp2))

                      ngb_dir = ngb_dir + &

                           sys(ngbindex(ngb))%dir(indsys2(izp2))

                    END IF

                  END DO

                END DO

                IF (count2.GT.0) THEN

                  absper = abs(sys(ss)%tp(indsys1(izp))-ngb_tp/count2)

                  absdir = abs(sys(ss)%dir(indsys1(izp))-ngb_dir/count2)

                  t = sys(ss)%tp(indsys1(izp))

                  m1 = -3.645*t + 63.211

                  m1 = max(m1,10.)

                  m2 = -0.346*t + 3.686

                  m2 = max(m2,0.6)

                  deltadirb = (m1*dx + dirknob)*1.

                  deltaperb = (m2*dx + perknob)*1.

                  IF ( (absper.LT.deltaperb).AND. &

                       (absdir.LT.deltadirb) ) THEN

                    ngbmatch(izp)=1.

                  END IF

                END IF

              END DO

              !                       If >80% of neighbors fall within criteria, system is matched

              IF ((sum(ngbmatch(1:lsys))/lsys).GT.0.50) THEN

                indmatch = ngbindex(ngb)

                matchsys = sys(indmatch)%sysInd

              ELSE

                cycle

              END IF

            ELSE

              !061512---------------------------------


              tb = max(sum(sys(ss)%tp(indsys1))/lsys, &

                   sum(sys(ngbindex(ngb))%tp(indsys2))/lsys2)

              !                      deltaPerB = (-0.06*Tb+2+perKnob)*1.5

              !                      deltaDirB = (-Tb+(25+10*dirKnob))*1.5

              !                      deltaPerB = (-0.06*Tb+2+2)*1.5

              !                      deltaDirB = (-Tb+(25+10*2))*1.5

              dx=0.5*((wsdat%lon(2,1)-wsdat%lon(1,1)) + &

                   (wsdat%lat(1,2)-wsdat%lat(1,1)))

              m1 = -3.523*tb + 64.081

              m1 = max(m1,10.)

              m2 = -0.337*tb + 3.732

              m2 = max(m2,0.6)

              !1stddev                      m1 = -2.219*Tb + 35.734

              !1stddev                      m1 = MAX(m1,5.)

              !1stddev                      m2 = -0.226*Tb + 2.213

              !1stddev                      m2 = MAX(m2,0.35)

              !5stddev                      m1 = -5.071*Tb + 90.688

              !5stddev                      m1 = MAX(m1,16.)

              !5stddev                      m2 = -0.467*Tb + 5.161

              !5stddev                      m2 = MAX(m2,1.0)

              deltadirb = (m1*1. + dirknob)*1.

              deltaperb = (m2*1. + perknob)*1.

              deltahsb = 0.50*sum(sys(ss)%hs(indsys1))/lsys

              !                      deltaHsB = 0.25*SUM(sys(ss)%hs(indSys1))/lsys


              !051612               --- Land mask addition

              !                     Option 1: If system centroid is near a land mask (e.g. 3 arc-deg),

              !                               increase the tolerances

              IF (any(mask.EQ.1)) THEN

                lonmean = sum(sys(ss)%lon(indsys1))/lsys

                latmean = sum(sys(ss)%lat(indsys1))/lsys

                DO j = 1,maxj

                  DO i = 1,maxi

                    IF (mask(i,j).EQ.1) THEN

                      !                             Land point found. Compute distance to system centroid

                      dist = sqrt((lonmean-wsdat%lon(i,j))**2 +&

                           (latmean-wsdat%lat(i,j))**2)

                      IF (dist.LT.3.) THEN

                        !                               System assumed to be influenced by land,

                        !                               increase tolerances to deltaDirB=30,deltaPerB=3

                        !                                deltaDirB = (m1*1. + 30)*1.

                        !                                deltaPerB = (m2*1. + 3)*1.

                        deltadirb = (m1*1. + 30)*1.

                        deltaperb = (m2*1. + 3)*1.

                        !Remove dHs limitation from criteria

                        deltahsb = 9999.

                        GOTO 500

                      END IF

                    END IF

                  END DO

                END DO

              END IF

500           CONTINUE

              !051612               --- Land mask addition


              abshs = abs( sum(sys(ss)%hs(indsys1))/lsys - &

                   sum(sys(ngbindex(ngb))%hs(indsys2))/lsys2 )

              absper = abs( sum(sys(ss)%tp(indsys1))/lsys - &

                   sum(sys(ngbindex(ngb))%tp(indsys2))/lsys2 )

              absdir = abs( &

                   mean_anglev2(sys(ss)%dir(indsys1),lsys) - &

                   mean_anglev2(sys(ngbindex(ngb))%dir(indsys2), &

                   lsys2) )

              IF (absdir.GT.180) absdir = 360.-absdir

              !                      absWf = ABS( SUM(sys(ss)%wf(indSys1))/lsys - &

              !                        SUM(sys(ngbIndex(ngb))%wf(indSys2))/lsys2 )


              IF ( (absper.LT.deltaperb).AND. &

                   (absdir.LT.deltadirb).AND. &

                   (abshs.LT.deltahsb) ) THEN

                indmatch = ngbindex(ngb)

                matchsys = sys(indmatch)%sysInd

                !opt                         WRITE(20,*) '-> Matched sys',s, &

                !opt                                     'with neighbor sys',matchSys

              ELSE

                cycle

              END IF

              !061512---------------------------------

            END IF

            !061512---------------------------------


            keep = 0

            keepind(:) = 0


            DO hh = 1, sys(ss)%nPoints

              ii = sys(ss)%i(hh)

              jj = sys(ss)%j(hh)

              ind = 0

              ind = findfirst(real(wsdat%par(ii,jj)%sys), &

                   SIZE(wsdat%par(ii,jj)%sys),real(s))        !Shouldn't REAL(s) be matchSys...

              IF (ind.NE.0) THEN

                wsdat%par(ii,jj)%sys(ind)=matchsys            !...and matchSys be s, (i.e. add the matching neigbour to the base?)

              END IF

              !                        Remove the "-1" system from the set

              ind2 = 1

              oneless(:) = 9999                                !Streamline this?

              leng = length(real(wsdat%par(ii,jj)%sys), &

                   SIZE(wsdat%par(ii,jj)%sys),real(9999))

              DO ind = 1, leng

                IF ( wsdat%par(ii,jj)%sys(ind).NE.-1 ) THEN

                  oneless(ind2) = wsdat%par(ii,jj)%sys(ind)

                  ind2 = ind2+1

                END IF

              END DO

              ind2 = ind2-1

              !                        Combine any partitions assigned to the same systems

              !                        Check for duplicates

              IF (ind2.EQ.0) &

                   WRITE(20,*) '***2.Calling UNIQUE w. len=0!'

              CALL unique(real(oneless(1:ind2)),ind2, &

                   uniarr,outsize)

              IF (ASSOCIATED(uniarr)) DEALLOCATE(uniarr)

              IF (ind2.GT.outsize) THEN

                !                          There is at least one duplicate, so combine systems

                CALL combinepartitionsv2(wsdat%par(ii,jj))

                !                          Update the combined partitions values into the system we are keeping.

                !                          Since partitions have been combined we don't know if the index is the same

                replacedind = &

                     findfirst(real(wsdat%par(ii,jj)%sys(:)), &

                     SIZE(wsdat%par(ii,jj)%sys(:)), &

                     REAL(matchSys))

                hhforindmatch = 1

                DO WHILE (hhforindmatch.LE. &

                     sys(indmatch)%nPoints)

                  IF ( (sys(indmatch)%i(hhforindmatch) &

                       .EQ.ii).AND. &

                       (sys(indmatch)%j(hhforindmatch) &

                       .EQ.jj) ) EXIT

                  hhforindmatch = hhforindmatch + 1

                END DO

                sys(indmatch)%hs(hhforindmatch) = &

                     wsdat%par(ii,jj)%hs(replacedind)

                sys(indmatch)%tp(hhforindmatch) = &

                     wsdat%par(ii,jj)%tp(replacedind)

                sys(indmatch)%dir(hhforindmatch) = &

                     wsdat%par(ii,jj)%dir(replacedind)

                sys(indmatch)%dspr(hhforindmatch) = &

                     wsdat%par(ii,jj)%dspr(replacedind)

                !                           sys(indMatch)%wf(hhForIndMatch) = &

                !                                  wsdat%par(ii,jj)%wf(replacedInd)

              ELSE

                keep = keep+1

                keepind(keep) = hh

              END IF

            END DO

            leng = length(real(sys(indmatch)%hs), &

                 SIZE(sys(indmatch)%hs),real(9999.))


            !                     Update system info

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

            !                     First need to find which points were common to both systems =>

            !                     keepInd since that means partitions have not been combined for those

            !                     points as a result of the combination of those 2 systems =>

            !                     distinct points

            !                      keepInd = keepInd(1:keep)

            lmatch = length(real(sys(indmatch)%hs), &

                 SIZE(sys(indmatch)%hs),real(9999.))

            tot = lmatch + keep

            CALL union (real(sys(indmatch)%ngbr), &

                 SIZE(sys(indmatch)%ngbr), &

                 REAL(sys(ss)%ngbr), &

                 SIZE(sys(ss)%ngbr), &

                 allngbr,outsize)

            CALL setdiff(allngbr,SIZE(allngbr), &

                 REAL((/sys(indMatch)%sysInd, &

                 sys(ss)%sysInd/)), &

                 SIZE((/sys(indMatch)%sysInd, &

                 sys(ss)%sysInd/)),difarr,outsize)

            sys(indmatch)%ngbr(:) = 9999

            outsize = min(outsize,size(sys(indmatch)%ngbr))

            sys(indmatch)%ngbr(1:outsize) = nint(difarr(1:outsize))

            IF (ASSOCIATED(allngbr)) DEALLOCATE(allngbr)

            IF (ASSOCIATED(difarr)) DEALLOCATE(difarr)


            leng = length(real(sys(indmatch)%i), &

                 SIZE(sys(indmatch)%i),real(9999))

            sys(indmatch)%hsMean = sum((/ &

                 sys(ss)%hs(keepind(1:keep)), &

                 sys(indmatch)%hs(1:leng) /))/tot

            sys(indmatch)%tpMean = sum((/ &

                 sys(ss)%tp(keepind(1:keep)), &

                 sys(indmatch)%tp(1:leng) /))/tot

            sys(indmatch)%dirMean = &

                 mean_anglev2((/ sys(ss)%dir(keepind(1:keep)), &

                 sys(indmatch)%dir(1:leng) /),tot)

            !070512----------- Weight averages with Hm0 ---------------------

            temp1 = 0.

            temp2 = 0.

            DO iii = 1,keep

              temp1 = temp1 + (sys(ss)%hs(keepind(iii))**2)*&

                   sys(ss)%hs(keepind(iii))

              temp2 = temp2 + (sys(ss)%hs(keepind(iii))**2)*&

                   sys(ss)%tp(keepind(iii))

            END DO

            DO iii = 1,leng

              temp1 = temp1 + (sys(indmatch)%hs(iii)**2)*&

                   sys(indmatch)%hs(iii)

              temp2 = temp2 + (sys(indmatch)%hs(iii)**2)*&

                   sys(indmatch)%tp(iii)

            END DO

            sys(indmatch)%hsMean = temp1/max(sum((/ &

                 sys(ss)%hs(keepind(1:keep))**2, &

                 sys(indmatch)%hs(1:leng)**2 /)),0.001)

            sys(indmatch)%tpMean = temp2/max(sum((/ &

                 sys(ss)%hs(keepind(1:keep))**2, &

                 sys(indmatch)%hs(1:leng)**2 /)),0.001)

            sys(indmatch)%dirMean = &

                 mean_anglev3((/ sys(ss)%dir(keepind(1:keep)), &

                 sys(indmatch)%dir(1:leng) /), &

                 (/ sys(ss)%hs(keepind(1:keep)), &

                 sys(indmatch)%hs(1:leng) /),tot)

            !070512----------- Weight averages with Hm0 ---------------------


            sys(indmatch)%i(1:(keep+leng))= &

                 (/sys(ss)%i(keepind(1:keep)), &

                 sys(indmatch)%i(1:leng)/)

            sys(indmatch)%j(1:(keep+leng))= &

                 (/sys(ss)%j(keepind(1:keep)), &

                 sys(indmatch)%j(1:leng)/)

            sys(indmatch)%lat(1:(keep+leng)) = &

                 (/sys(ss)%lat(keepind(1:keep)), &

                 sys(indmatch)%lat(1:leng)/)

            sys(indmatch)%lon(1:(keep+leng)) = &

                 (/sys(ss)%lon(keepind(1:keep)), &

                 sys(indmatch)%lon(1:leng)/)

            sys(indmatch)%dir(1:(keep+leng)) = &

                 (/sys(ss)%dir(keepind(1:keep)), &

                 sys(indmatch)%dir(1:leng)/)

            sys(indmatch)%dspr(1:(keep+leng)) = &

                 (/sys(ss)%dspr(keepind(1:keep)), &

                 sys(indmatch)%dspr(1:leng)/)

            !                      sys(indMatch)%wf(1:(keep+leng)) = &

            !                          (/sys(ss)%wf(keepInd(1:keep)), &

            !                            sys(indMatch)%wf(1:leng)/)

            sys(indmatch)%hs(1:(keep+leng)) = &

                 (/sys(ss)%hs(keepind(1:keep)), &

                 sys(indmatch)%hs(1:leng)/)

            sys(indmatch)%tp(1:(keep+leng)) = &

                 (/sys(ss)%tp(keepind(1:keep)), &

                 sys(indmatch)%tp(1:leng)/)

            sys(indmatch)%nPoints = &

                 length(real(sys(indmatch)%i), &

                 SIZE(sys(indmatch)%i),real(9999))

            !                     Clear array of system that has just been combined with another

            sys(ss)%nPoints = 0

            sys(ss)%ngbr(:) = 9999

            WRITE(20,*) 'Deallocating sys',s

            DEALLOCATE( sys(ss)%hs )                            !opt

            DEALLOCATE( sys(ss)%tp )                            !opt

            DEALLOCATE( sys(ss)%dir )                           !opt

            DEALLOCATE( sys(ss)%dspr )                          !opt

            !                      DEALLOCATE( sys(ss)%wf )                           !opt

            DEALLOCATE( sys(ss)%i )                             !opt

            DEALLOCATE( sys(ss)%j )                             !opt

            DEALLOCATE( sys(ss)%lat )                           !opt

            DEALLOCATE( sys(ss)%lon )                           !opt

            !                      DEALLOCATE( sys(ss)%hsMean )                        !opt

            !                      DEALLOCATE( sys(ss)%tpMean )                        !opt

            !                      DEALLOCATE( sys(ss)%dirMean )                       !opt


            !                     Loop through wsdat to update neighbouring system values

            DO i = 1, maxi

              DO j = 1, maxj

                ind = findfirst(real(wsdat%par(i,j)%ngbrSys), &

                     SIZE(wsdat%par(i,j)%ngbrSys),real(s))

                IF (ind.NE.0) THEN

                  wsdat%par(i,j)%ngbrSys(ind)=matchsys

                END IF

                leng = length(real(wsdat%par(i,j)%ngbrSys), &

                     SIZE(wsdat%par(i,j)%ngbrSys),real(9999))

                IF (leng.GT.0) THEN

                  CALL unique( &

                       REAL(wsdat%par(i,j)%ngbrSys(1:leng)), &

                       leng,uniarr,outsize)

                  wsdat%par(i,j)%ngbrSys(:) = 9999

                  wsdat%par(i,j)%ngbrSys(1:outsize) = &

                       nint(uniarr)

                  IF (ASSOCIATED(uniarr)) DEALLOCATE(uniarr)

                ELSE

                  wsdat%par(i,j)%ngbrSys(:) = 9999

                END IF

              END DO

            END DO


            !                     Update neigbors in sys structure

            DO nn = 1, maxsys

              nbr = findfirst(real(sys(nn)%ngbr), &

                   SIZE(sys(nn)%ngbr),real(s))

              IF (nbr.NE.0) THEN

                !                             WRITE(20,*) 'update'

                sys(nn)%ngbr(nbr)=matchsys

              END IF

              leng2 = length(real(sys(nn)%ngbr), &

                   SIZE(sys(nn)%ngbr),real(9999))

              IF (leng2.GT.0) THEN

                CALL unique(real(sys(nn)%ngbr(1:leng2)), &

                     leng2,uniarr,outsize)

                sys(nn)%ngbr(:) = 9999

                sys(nn)%ngbr(1:outsize) = nint(uniarr)

                IF (ASSOCIATED(uniarr)) DEALLOCATE(uniarr)

                !                             WRITE(20,*) 'has now ngbr: ', &

                !                                sys(nn)%ngbr(1:outsize)

              END IF

            END DO

            EXIT

          END IF


          IF (ASSOCIATED(indsys1)) DEALLOCATE(indsys1)

          IF (ASSOCIATED(indsys2)) DEALLOCATE(indsys2)

        END DO

      END IF

    END DO


    IF (ALLOCATED(sysordered)) DEALLOCATE(sysordered)

    IF (ALLOCATED(syssortedind)) DEALLOCATE(syssortedind)

    IF (ALLOCATED(sysout)) DEALLOCATE(sysout)


    !     Compile array index of active systems in sys

    actsys = 0

    DO ic = 1,maxsys

      IF (sys(ic)%nPoints>0) actsys = actsys + 1

    END DO

    IF (ALLOCATED(actsysind)) DEALLOCATE(actsysind)

    ALLOCATE( actsysind(actsys) )

    actsys = 0

    DO ic = 1,maxsys

      IF (sys(ic)%nPoints>0) THEN

        actsys = actsys + 1

        actsysind(actsys) = sys(ic)%sysInd

      END IF

    END DO


    !opt      WRITE(20,*) 'actSys =',actSys

    !opt      WRITE(20,*) 'actSysInd =',actSysInd

    !opt      DO ic = 1,SIZE(actSysInd)

    !opt         s = actSysInd(ic)

    !opt         WRITE(20,*) 'sys(',s,')%sysInd =',sys(s)%sysInd

    !opt      END DO

    WRITE(20,*) 'Leaving combineSys...'


    RETURN

  END SUBROUTINE combinesys

  !/ End of combineSys ------------------------------------------------- /

  !/

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

  SUBROUTINE combinepartitionsv2 (dat)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          4-Jan-2013 |

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

    !/

    !/    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 )

    !/

    !/    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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     Combine two partitions that have been assigned to the same system

    !

    !  2. Method

    !

    !     Of all the partitions associated with a certain common system,

    !     add all the Hs values to the partition with the largest Hs,

    !     and delete the rest. NOTE that the tp and dir values of this

    !     maximum partition is not adjusted!

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     dat    TYPE(param)   in/out  Input data structure (partitions set)

    !                                  to combine

    !

    TYPE(param) :: dat


    INTENT (IN OUT) dat

    !

    !     Local variables

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

    TYPE duplicate

      INTEGER :: val

      INTEGER :: ndup

      INTEGER :: ind(50)

    END TYPE duplicate


    TYPE(duplicate) :: dup(100)                                         !40.PAR

    LOGICAL :: found

    INTEGER :: nsys, ndup, p, pp, maxInd, npart, s, ss, ppp

    REAL :: temp

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      -

    !

    !  5. Subroutines calling

    !

    !     findSys

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     -

    !

    !  9. Switches :

    !

    !     None defined yet.

    !

    ! 10. Source code :

    !

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


    !     Find indices in dat%sys(:) of all partition associated with

    !     the same wave system, and store them in the data structure

    !     dup(1:nsys). Here nsys is the number of systems for which duplicates

    !     were found, and dup(s)%ndup the number of partitions assigned

    !     to the same system s.

    nsys = 0

    dup(:)%ndup = 0

    dup(:)%val = 9999

    DO s = 1,100

      dup(s)%ind(:) = 0

    END DO


    npart = length(real(dat%ipart),SIZE(dat%ipart),real(0))

    DO p = 1, npart-1

      found = .false.

      IF (any(dat%sys(p).EQ.dup(:)%val)) cycle                         !found = .TRUE.

      DO pp = (p+1), npart

        IF (dat%sys(p).EQ.dat%sys(pp)) THEN

          !              First value

          IF (.NOT.found) THEN

            nsys=nsys+1

            dup(nsys)%val = dat%sys(p)

            dup(nsys)%ndup = 1

            dup(nsys)%ind(dup(nsys)%ndup) = p

            found = .true.

          END IF

          !              Subsequent duplicates

          IF (.NOT.any(pp.EQ.dup(nsys)%ind(:))) THEN

            dup(nsys)%ndup = dup(nsys)%ndup+1

            dup(nsys)%ind(dup(nsys)%ndup) = pp

          END IF

        END IF

      END DO

    END DO


    !     Now go through array of duplicates for each of n systems

    !     to add all the wave energy to the most energetic of the

    !     duplicates, and then remove the rest.

    maxind = 0

    temp = -9999.

    DO s = 1, nsys

      !        Find duplicate partition with the largest Hs (most energy)

      DO p = 1, dup(s)%ndup

        IF ( temp.LT.dat%hs(dup(s)%ind(p)) ) THEN

          temp = dat%hs(dup(s)%ind(p))

          maxind = p

        END IF

      END DO


      !        Add all energy (Hs) to this partition

      dat%hs(dup(s)%ind(maxind)) = &

           sqrt( sum(dat%hs(dup(s)%ind(1:dup(s)%ndup))**2) )


      !        Remove duplicate partitions which did not have the maximum Hs,

      !        and shift up indices to fill the gap

      DO p = 1, dup(s)%ndup

        !           Find index to remove

        IF (p.NE.maxind) THEN

          !              Shift up entries, deleting the duplicate partition

          !              REPLACE WITH CSHIFT(ARRAY, SHIFT, dim) ?

          dat%hs( dup(s)%ind(p):(npart-1) ) = &

               dat%hs( (dup(s)%ind(p)+1):npart)

          dat%tp( dup(s)%ind(p):(npart-1) ) = &

               dat%tp( (dup(s)%ind(p)+1):npart)

          dat%dir( dup(s)%ind(p):(npart-1) ) = &

               dat%dir( (dup(s)%ind(p)+1):npart)

          dat%dspr( dup(s)%ind(p):(npart-1) ) = &

               dat%dspr( (dup(s)%ind(p)+1):npart)

          !               dat%wf( dup(s)%ind(p):(npart-1) ) = &

          !                       dat%wf( (dup(s)%ind(p)+1):npart)

          dat%sys( dup(s)%ind(p):(npart-1) ) = &

               dat%sys( (dup(s)%ind(p)+1):npart)

          dat%ipart( dup(s)%ind(p):(npart-1) ) = &

               dat%ipart( (dup(s)%ind(p)+1):npart)

          !              Shift up indices

          DO ss = 1, nsys

            DO ppp = 1, dup(ss)%ndup

              IF (dup(ss)%ind(ppp).GT.dup(s)%ind(p)) &

                   dup(ss)%ind(ppp) = dup(ss)%ind(ppp)-1

            END DO

          END DO

          !              Add blank to end

          dat%hs(npart) = 9999.

          dat%tp(npart) = 9999.

          dat%dir(npart) = 9999.

          dat%dspr(npart) = 9999.

          !               dat%wf(npart) = 9999.

          dat%sys(npart) = 9999

          dat%ipart(npart) = 0

        END IF

      END DO

    END DO


    RETURN

  END SUBROUTINE combinepartitionsv2

  !/ End of combinePartitionsV2 ---------------------------------------- /

  !/

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

  REAL FUNCTION mean_angleV2(ang,ll)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          4-Jan-2013 |

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

    !/

    !/    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 )

    !/

    !/    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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     Compute the mean direction from array of directions

    !

    !  2. Method

    !

    !     ang is a column vector of angles

    !     m_ang is the mean from a unit-vector average of ang

    !     Assumes clockwise rotation from North = 0.

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     ang     Real   input    Array of angles to average

    !     ll      Int    input    Length of ang

    !

    REAL :: ang(ll)

    INTEGER :: ll

    !

    !     Local variables

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

    !     u,v       Real   Arrays of u,v dir components to average

    !     um,vm     Real   Mean u,v dir components

    !     theta     Real   Mean direction relative to North

    !

    REAL      :: pi

    parameter(pi = 3.1416)

    REAL      :: u(ll), v(ll), vm, um, theta

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      -

    !

    !  5. Subroutines calling

    !

    !     findSys

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     -

    !

    !  9. Switches :

    !

    !     None defined yet.

    !

    ! 10. Source code :

    !

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


    !     North and East components

    v(:) = cos(ang(:)*(pi/180.))

    u(:) = sin(ang(:)*(pi/180.))

    vm = sum(v)/ll

    um = sum(u)/ll


    !     Compute mean magnitude and direction relative to North (from Upolar.m)

    theta = (atan2(um,vm))*(180/pi)


    !     Convert inputs to radians, the to the -pi to pi range

    !     (incorporated from original function xunwrapV2.m)


    !     Convert to radians

    theta = theta*(pi/180)


    theta = pi*((abs(theta)/pi) - &

         2*ceiling(((abs(theta)/pi)-1)/2))*sign(1.,theta)


    !     Shift the points in the -pi to 0 range to the pi to 2pi range

    IF (theta.LT.0.) theta = theta + 2*pi


    !     Convert back to degrees and return value

    mean_anglev2 = theta*(180/pi)


    RETURN

  END FUNCTION mean_anglev2

  !/ End of mean_angleV2 ----------------------------------------------- /

  !/

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

  REAL function mean_anglev3(ang,hsign,ll)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          4-Jan-2013 |

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

    !/

    !/    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 )

    !/

    !/    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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     Compute the mean direction from array of directions,

    !     INCLUDING WEIGHTING WITH HMO

    !

    !  2. Method

    !

    !     ang is a column vector of angles

    !     m_ang is the mean from a unit-vector average of ang

    !     Assumes clockwise rotation from North = 0.

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     ang     Real   input    Array of angles to average

    !     ll      Int    input    Length of ang

    !

    REAL :: ang(ll), hsign(ll)

    REAL :: temp1, temp2

    INTEGER :: ll

    !

    !     Local variables

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

    !     u,v       Real   Arrays of u,v dir components to average

    !     um,vm     Real   Mean u,v dir components

    !     theta     Real   Mean direction relative to North

    !

    REAL      :: pi

    parameter(pi = 3.1416)

    REAL      :: u(ll), v(ll), vm, um, theta

    INTEGER   :: i

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      -

    !

    !  5. Subroutines calling

    !

    !     findSys

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     -

    !

    !  9. Switches :

    !

    !     None defined yet.

    !

    ! 10. Source code :

    !

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


    !     North and East components

    v(:) = cos(ang(:)*(pi/180.))

    u(:) = sin(ang(:)*(pi/180.))

    temp1 = 0.

    temp2 = 0.

    DO i = 1,ll

      temp1 = temp1 + (hsign(i)**2)*v(i)

      temp2 = temp2 + (hsign(i)**2)*u(i)

    END DO

    vm = temp1/max(sum(hsign**2),0.001)

    um = temp2/max(sum(hsign**2),0.001)


    !     Compute mean magnitude and direction relative to North (from Upolar.m)

    theta = (atan2(um,vm))*(180/pi)


    !     Convert inputs to radians, the to the -pi to pi range

    !     (incorporated from original function xunwrapV2.m)


    !     Convert to radians

    theta = theta*(pi/180)


    theta = pi*((abs(theta)/pi) - &

         2*ceiling(((abs(theta)/pi)-1)/2))*sign(1.,theta)


    !     Shift the points in the -pi to 0 range to the pi to 2pi range

    IF (theta.LT.0.) theta = theta + 2*pi


    !     Convert back to degrees and return value

    mean_anglev3 = theta*(180/pi)


    RETURN

  END FUNCTION mean_anglev3

  !/ End of mean_angleV3 ----------------------------------------------- /

  !/

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

  SUBROUTINE unique (INARRAY,INSIZE,OUTARRAY,OUTSIZE)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :         22-Dec-2016 |

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

    !/

    !/    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 )

    !/    12-Dec-2016 : Change algorithm from N*N to N*log(N)

    !/                  (S. Zieger BoM Australia)           ( 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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     Returns the sorted elements that are unique in INARRAY.

    !

    !  2. Method

    !

    !     1. Sort input array with quicksort

    !     2. Copy sequential-elements if the 'current' element

    !        is not equal to the 'previous' element in array.

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     INARRAY   REAL    ARR  input    Input array

    !     INSIZE    INTEGER      input    Size of input array

    !     OUTARRAY  REAL    ARR  output   Output array (sorted)

    !     OUTSIZE   INTEGER      output   Size of output array (number of unique elements)

    !

    INTEGER, INTENT(IN)  :: INSIZE

    INTEGER, INTENT(OUT) :: OUTSIZE

    REAL, INTENT(IN)     :: INARRAY(INSIZE)

    REAL, POINTER        :: OUTARRAY(:)

    !

    !     Local variables

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

    INTEGER          :: I, K

    REAL             :: ARRAY(INSIZE), TEMP(INSIZE)

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Scope    Description

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

    !      QSORT     Subr. Private  Quicksort algorithm

    !

    !  5. Subroutines calling

    !

    !     waveTracking_NWS_V2

    !     findSys

    !     printFinalSys

    !     combineSys

    !     findIJV4

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     -

    !

    !  9. Switches :

    !

    !     None defined yet.

    !

    ! 10. Source code :

    !

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

    k = 1


    IF ( insize.EQ.0 ) THEN

      WRITE(20,*) '*** In Subr. UNIQUE: Input array has length=0!'

    ELSE

      !/    --- Setup input arrays and temporary arrays. ---

      DO i=1,insize

        array(i) = inarray(i)

        temp(i)  = real(i)

      END DO

      !/

      !/      --- Sort input arrays (use temporary array to store indices). ---

      CALL qsort(array,temp,1,insize)

      !/

      !/      --- Reset temporary array. ---

      temp(:) = 9999.

      !/

      !/      --- Initialise first values and array index. ---

      k = 1

      temp(k) = array(k)

      k = k + 1

      !/      --- Iterate over elements in array. ---

      DO i=2,insize

        !/        --- Compare sequential array values ('previous' less than 'next')

        !/            and test against the last list element check in. ---

        IF ( array(i).GT.array(i-1) .AND. &

             array(i).GT.temp(k-1)        ) THEN

          temp(k) = array(i)

          k = k + 1

        END IF

      END DO


      !/      --- Allocate output array ---

      outsize = k - 1

      ALLOCATE(outarray(outsize))

      !/      --- Transfer output from temporary array to output array. ---

      IF ( outsize.GE.1 ) THEN

        DO i=1,outsize

          outarray(i) = temp(i)

        END DO

      END IF

    END IF

    !/

    RETURN

    !/

  END SUBROUTINE unique

  !/ End of UNIQUE ----------------------------------------------------- /

  !/

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

  SUBROUTINE sort (INARRAY,INSIZE,OUTARRAY,IY,DIRECTION)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :         20-Dec-2016 |

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

    !/

    !/    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 )

    !/    20-Dec-2016 : Add quicksort algorithm (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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     Sorts the array INARRAY in ascending (Direction = 'A') or

    !     descending (Direciton = 'D') order. The sorted array is

    !     stored in OUTARRAY, and the sorted array of the original

    !     indices is stored in IY.

    !

    !  2. Method

    !

    !     Sort algorithm based on quicksort.

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     INARRAY   REAL    ARR  input    Input array

    !     INSIZE    INTEGER      input    Size of input array

    !     OUTARRAY  REAL    ARR  output   Sorted output array

    !     IY        INTEGER ARR  output   Sorted array of the original indices


    CHARACTER         :: DIRECTION *1

    INTEGER           :: INSIZE

    INTEGER           :: IY(INSIZE)

    REAL              :: INARRAY(INSIZE), OUTARRAY(INSIZE)


    INTENT (IN)  inarray, insize, direction

    INTENT (OUT) outarray, iy

    !

    !     Local variables

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

    !      INARRAY - array of values to be sorted

    !      IY -      array to be carried with X (all swaps of X elements are       ??? EDIT!

    !                matched in IY .  After the sort IY(J) contains the original

    !                postition of the value X(J) in the unsorted X array.

    !      N -       number of values in array X to be sorted


    INTEGER           :: I

    REAL              :: IND(INSIZE)

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      -

    !

    !  5. Subroutines calling

    !

    !     printFinalSys

    !     combineSys

    !     timeTrackingV2

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     -

    !

    !  9. Switches :

    !

    !     None defined yet.

    !

    ! 10. Source code :

    !

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


    !     Sort OUTARRAY in as/decending order


    IF (insize.EQ.0) THEN

      WRITE(20,*) '*** In Subr. SORT: Input array has length=0 !!!'

    ELSE


      DO i = 1, insize

        outarray(i) = inarray(i)

        ind(i) = real(i)

      END DO


      IF (direction .EQ. 'A') THEN

        CALL qsort(outarray,ind,1,insize)

      ELSE IF (direction .EQ. 'D') THEN

        CALL qsort_desc(outarray,ind,1,insize)

      END IF


    END IF

    !

    !/    --- Cast index array to integer. ---

    DO i = 1, insize

      iy(i) = int(ind(i))

    END DO

    !

    RETURN

    !

  END SUBROUTINE sort

  !/ End of SORT ------------------------------------------------------- /

  !/

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

  SUBROUTINE setdiff (INARRAY1,   INSIZE1,  INARRAY2,   INSIZE2, &

       OUTARRAY,   OUTSIZE)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :         20-Dec-2016 |

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

    !/

    !/    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 )

    !/    20-Dec-2016 : Add quicksort algorithm (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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     (i)  Returns the elements in INARRAY1 that are not in INARRAY2.

    !     (ii) Sort the resulting array in ascending order.

    !

    !  2. Method

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     INARRAY1  REAL    ARR  input    Input array

    !     INSIZE1   INTEGER      input    Size of input array

    !     INARRAY2  REAL    ARR  input    Input array

    !     INSIZE2   INTEGER      input    Size of input array

    !     OUTARRAY  REAL    ARR  output   Output array

    !     OUTSIZE   INTEGER      output   Size of output array (number of unique elements)


    INTEGER           :: INSIZE1, INSIZE2, OUTSIZE

    REAL              :: INARRAY1(INSIZE1), INARRAY2(INSIZE2)

    REAL, POINTER     :: OUTARRAY(:)


    INTENT (IN)  inarray1, insize1, inarray2, insize2

    INTENT (OUT) outsize

    !

    !     Local variables

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

    INTEGER          :: I,J,K

    REAL             :: TEMP(INSIZE1)

    REAL             :: ARRAY1(INSIZE1),ARRAY2(INSIZE2)

    REAL             :: ID1(INSIZE1),ID2(INSIZE2)

    LOGICAL          :: LOOP

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Scope    Description

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

    !      QSORT     Subr. Private  Quicksort algorithm

    !

    !  5. Subroutines calling

    !

    !     printFinalSys

    !     combineSys

    !     timeTrackingV2

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     -

    !

    !  9. Switches :

    !

    !     None defined yet.

    !

    ! 10. Source code :

    !

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

    IF ( (insize1).EQ.0 ) THEN

      outsize = 0

      ALLOCATE(outarray(outsize))

    ELSE IF ( insize2.EQ.0 ) THEN

      CALL unique(inarray1,insize1,outarray,outsize)

    ELSE

      !/      --- Setup input arrays. ---

      DO i=1,insize1

        array1(i) = inarray1(i)

        id1(i)    = real(i)

      END DO

      DO i=1,insize2

        array2(i) = inarray2(i)

        id2(i)    = real(i)

      END DO

      !/

      !/      --- Sort input arrays. ---

      CALL qsort(array1,id1,1,insize1)

      CALL qsort(array2,id2,1,insize2)

      !/

      !/      --- Initialise indices. ---

      i = 1

      j = 1

      k = 1

      !/

      !/      --- Allocate and initialize temporary output ---

      temp(:) = 9999.

      !/

      !/      --- Loop though both arrays by incrementing I,J. ---

      loop = .true.

      DO WHILE ( loop )

        !/

        IF ( array1(i).LT.array2(j) .OR.   &

             array1(i).GT.array2(insize2)  ) THEN

          !/            --- Populate output array. Check for dumplicates

          !/                in output array. ---

          IF ( k.EQ.1 ) THEN

            temp(k) = array1(i)

            k = k + 1

          ELSE IF ( temp(k-1).LT.array1(i) ) THEN

            temp(k) = array1(i)

            k = k + 1

          END IF

          i = i + 1

        ELSE IF ( array2(j).LT.array1(i) ) THEN

          j = j + 1

        ELSE

          i = i + 1

          j = j + 1

        END IF

        !/          --- Check for exit the loop. ---

        IF ( i.GT.insize1 ) THEN

          loop = .false.

        END IF

        !/          --- Make sure array pointer I,J are within array bounds. ---

        i = min(i,insize1)

        j = min(j,insize2)

        !/

      END DO

      !/

      !/      --- Allocate output array ---

      outsize = k-1

      ALLOCATE(outarray(outsize))

      !/        --- Transfer output from temporary array to output array. ---

      DO i=1,outsize

        outarray(i) = temp(i)

      END DO

    END IF

    !/

    RETURN

    !/

  END SUBROUTINE setdiff

  !/ End of SETDIFF ---------------------------------------------------- /

  !/

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

  SUBROUTINE intersect (INARRAY1 ,INSIZE1  ,INARRAY2 ,INSIZE2  , &

       OUTARRAY ,OUTSIZE  ,IND1     ,IND2     )

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :         20-Dec-2016 |

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

    !/

    !/    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 )

    !/    20-Dec-2016 : Add count-histogram method based on

    !/                  algorithm from Mirko Velic (BoM)

    !/                  (S. Zieger BoM, Australia)          ( 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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     (i)  Returns the elements that are mutual in INARRAY1 and INARRAY2.

    !     (ii) Sort the resulting array in ascending order.

    !

    !  2. Method

    !

    !     Sort with counting/histogram method with input array being

    !     cast as integer.

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     INARRAY1  REAL    ARR  input    Input array

    !     INSIZE1   INTEGER      input    Size of input array

    !     INARRAY2  REAL    ARR  input    Input array

    !     INSIZE2   INTEGER      input    Size of input array

    !     OUTARRAY  REAL    ARR  output   Output array

    !     OUTSIZE   INTEGER      output   Size of output array (number of

    !     intersects)

    !

    INTEGER           :: INSIZE1, INSIZE2, OUTSIZE

    REAL              :: INARRAY1(INSIZE1), INARRAY2(INSIZE2)

    REAL, POINTER :: OUTARRAY(:)

    REAL, POINTER :: IND1(:), IND2(:)

    !

    INTENT (IN)  inarray1, insize1, inarray2, insize2

    INTENT (OUT) outsize

    !

    !     Local variables

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

    !     VIDX1, VIDX2 - array(s) in which the value is represented by

    !                    its index (i.e. histogram with frequency 1)

    !     N            - data range and size of possible intersections.

    !

    LOGICAL,ALLOCATABLE :: VIDX1(:),VIDX2(:)

    INTEGER,ALLOCATABLE :: IPOS1(:),IPOS2(:)

    !

    INTEGER             :: I, J

    INTEGER             :: N, IMIN, IMAX

    INTEGER             :: MINV1,MAXV1, MINV2, MAXV2

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Scope    Description

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

    !

    !  5. Subroutines calling

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     -

    !

    !  9. Switches :

    !

    !     None defined yet.

    !

    ! 10. Source code :

    !

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

    !

    outsize = 0


    !/    --- Calculate the range of the two sets. ---

    minv1  = int(minval(inarray1))

    maxv1  = int(maxval(inarray1))

    minv2  = int(minval(inarray2))

    maxv2  = int(maxval(inarray2))


    !/    --- Check if ranges overlap. ---

    IF ( maxv1.LT.minv2.OR.insize1.EQ.0.OR.insize2.EQ.0 ) THEN

      ALLOCATE(outarray(outsize))

      ALLOCATE(ind1(outsize))

      ALLOCATE(ind2(outsize))

    ELSE

      !/      --- Calculate size of temporary output arrays. Allow

      !/          extra elements: ZERO, and make sure index is 1:N. ---

      imin = min(minv1,minv2)-1

      imax = max(maxv1,maxv2)+1


      n = imax-imin


      ALLOCATE(vidx1(n),vidx2(n))

      ALLOCATE(ipos1(n),ipos2(n))


      vidx1(1:n) = .false.

      vidx2(1:n) = .false.


      DO i=1,insize1

        j = int(inarray1(i)-imin)

        vidx1(j) = .true.

        ipos1(j) = i

      END DO


      DO i=1,insize2

        j = int(inarray2(i)-imin)

        !/        --- Intersect arrays and check for

        !/            duplicate elements in array2. ---

        IF  ( vidx1(j).AND..NOT.vidx2(j) ) THEN

          outsize = outsize + 1

          vidx2(j) = .true.

          ipos2(j) = i

        END IF

      END DO

      !/      --- Allocate output arrays. ---

      ALLOCATE(outarray(outsize))

      ALLOCATE(ind1(outsize))

      ALLOCATE(ind2(outsize))

      !/      --- Transfer contents. ---

      i = 1

      DO j=1,n

        IF ( vidx1(j).AND.vidx2(j).AND.i.LE.outsize ) THEN

          outarray(i) = inarray1(ipos1(j))

          ind1(i)     = ipos1(j)

          ind2(i)     = ipos2(j)

          i = i + 1

        END IF

      END DO

      !/      --- Free memory. ---

      DEALLOCATE(vidx1,vidx2)

      DEALLOCATE(ipos1,ipos2)

    END IF

    !/

    RETURN

    !/

  END SUBROUTINE intersect

  !/ End of INTERSECT -------------------------------------------------- /

  !/

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

  SUBROUTINE union (INARRAY1,   INSIZE1,  INARRAY2,   INSIZE2, &

       OUTARRAY,   OUTSIZE)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          4-Jan-2013 |

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

    !/

    !/    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 )

    !/    20-Dec-2016 : Add count-histogram method similarly

    !/                  to INTERSECT (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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     (i)  Returns the union of INARRAY1 and INARRAY2.

    !     (ii) Sort the resulting array in ascending order.

    !

    !  2. Method

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     INARRAY   REAL    ARR  input    Input array

    !     INSIZE    INTEGER      input    Size of input array

    !     OUTARRAY  REAL    ARR  output   Output array (sorted)

    !     OUTSIZE   INTEGER      output   Size of output array (number of

    !                                     unique elements)

    INTEGER           :: INSIZE1, INSIZE2, OUTSIZE

    REAL              :: INARRAY1(INSIZE1), INARRAY2(INSIZE2)

    REAL, POINTER     :: OUTARRAY(:)

    !

    INTENT (IN)  inarray1, insize1, inarray2, insize2

    INTENT (OUT) outsize

    !

    !     Local variables

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

    !     VIDX1, VIDX2 - array(s) in which the value is represented by

    !                    its index (i.e. histogram with frequency 1)

    !     N            - data range and size of possible intersections.

    !

    LOGICAL,ALLOCATABLE :: VIDX1(:),VIDX2(:)

    INTEGER,ALLOCATABLE :: IPOS1(:),IPOS2(:)

    REAL,ALLOCATABLE    :: TEMP(:)

    !

    INTEGER             :: I, J

    INTEGER             :: N, IMIN, IMAX

    INTEGER             :: MINV1,MAXV1, MINV2, MAXV2

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      QSORT     Subr. Private  Quicksort algorithm

    !

    !  5. Subroutines calling

    !

    !     combineSys

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     -

    !

    !  9. Switches :

    !

    !     None defined yet.

    !

    ! 10. Source code :

    !

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

    !/    --- Setup input arrays. ---

    IF ( (insize1+insize2).EQ.0 ) THEN

      outsize = 0

      ALLOCATE(outarray(outsize))


    ELSEIF ( insize1.EQ.0 ) THEN

      outsize = insize2

      ALLOCATE(outarray(outsize))

      ALLOCATE(temp(outsize))


      DO i=1,outsize

        outarray(i) = inarray2(i)

        temp(i)     = real(i)

      END DO

      CALL qsort(outarray,temp,1,outsize)


    ELSEIF ( insize2.EQ.0 ) THEN

      outsize = insize1

      ALLOCATE(outarray(outsize),temp(outsize))


      DO i=1,outsize

        outarray(i) = inarray1(i)

        temp(i)     = real(i)

      END DO

      CALL qsort(outarray,temp,1,outsize)


    ELSE

      outsize = 0

      !/      --- Calculate the range of the two sets. ---

      minv1  = int(minval(inarray1))

      maxv1  = int(maxval(inarray1))

      minv2  = int(minval(inarray2))

      maxv2  = int(maxval(inarray2))

      !

      !/  --- Allow extra elementes: ZERO, and make sure index is 1:N. ---

      imin = min(minv1,minv2)-1

      imax = max(maxv1,maxv2)+1


      n = imax-imin


      ALLOCATE(vidx1(n),vidx2(n))

      ALLOCATE(ipos1(n),ipos2(n))


      vidx1(1:n) = .false.

      vidx2(1:n) = .false.

      ipos1(1:n) = -9999

      ipos2(1:n) = -9999

      !/

      DO i=1,insize1

        j = int(inarray1(i)-imin)

        IF  ( .NOT.vidx1(j) ) THEN

          outsize = outsize + 1

          vidx1(j) = .true.

          ipos1(j) = i

        END IF

      END DO


      DO i=1,insize2

        j = int(inarray2(i)-imin)

        IF  ( .NOT.vidx1(j).AND..NOT.vidx2(j) ) THEN

          outsize = outsize + 1

          vidx2(j) = .true.

          ipos2(j) = i

        END IF

      END DO


      ALLOCATE(outarray(outsize))


      i = 1

      DO j=1,n

        IF ( vidx1(j).AND.i.LE.outsize ) THEN

          outarray(i) = inarray1(ipos1(j))

          i = i + 1

        ELSEIF ( vidx2(j).AND.i.LE.outsize ) THEN

          outarray(i) = inarray2(ipos2(j))

          i = i + 1

        END IF

      END DO


      DEALLOCATE(vidx1,vidx2)

      DEALLOCATE(ipos1,ipos2)


    END IF

    !/

    RETURN

    !/

  END SUBROUTINE union

  !/ End of UNION ------------------------------------------------------ /

  !/

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

  INTEGER FUNCTION length(ARRAY,ARRSIZE,VAL)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          4-Jan-2013 |

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

    !/

    !/    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 )

    !/

    !/    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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     Find largest index in ARRAY with a value not equal to the

    !     filler value VAL.

    !     E.g. If VAL = 9999. and ARRAY = [X X X X 9999. 9999. 9999.],

    !     the function returns 4.

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    INTEGER :: arrsize

    REAL :: array(arrsize)

    REAL :: val

    !

    !     Local variables

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

    REAL :: field

    INTEGER :: i

    !

    ! 10. Source code :

    !

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

    IF (arrsize.GT.0) THEN

      i = 1

      field = array(i)

      DO WHILE (field.NE.val)

        i = i+1

        IF (i.GT.SIZE(array)) EXIT

        field = array(i)

      END DO

      length = i-1

    ELSE

      length = 0

    END IF


    RETURN

  END FUNCTION length

  !/ End of LENGTH ----------------------------------------------------- /

  !/

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

  INTEGER FUNCTION findfirst(ARRAY,ARRSIZE,VAL)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          4-Jan-2013 |

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

    !/

    !/    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 )

    !/

    !/    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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     Fast algorithm to find the *first* index IND in ARRAY

    !     for which ARRAY(IND) = VAL. Use only when there are

    !     no duplicates in ARRAY!

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    INTEGER :: arrsize

    REAL :: array(arrsize)

    REAL :: val

    !

    !     Local variables

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

    INTEGER :: ind

    !

    ! 10. Source code :

    !

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

    ind = 1

    DO WHILE (ind.LE.arrsize)

      IF ( array(ind).EQ.val ) EXIT

      ind = ind + 1

    END DO

    IF (ind.GT.arrsize) THEN

      findfirst = 0

    ELSE

      findfirst = ind

    ENDIF


    RETURN

  END FUNCTION findfirst

  !/ End of FINDFIRST -------------------------------------------------- /

  !/

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

  REAL function std(array,n)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          4-Jan-2013 |

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

    !/

    !/    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 )

    !/

    !/    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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     Computes standard deviation.

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     ARRAY    REAL    Input array for which to compute the std dev.

    !     N        INT     Size of ARRAY

    !

    REAL    :: array(n)

    INTEGER :: n

    !

    !     Local variables

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

    REAL    :: mn

    !

    ! 10. Source code :

    !

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

    IF (n.GT.1) THEN

      mn = sum(array)/n

      std = sqrt( 1/(real(n)-1)*sum( (array(:)-mn)**2 ) )

    ELSE

      std = 0.

    END IF


    RETURN

  END FUNCTION std

  !/ End of STD -------------------------------------------------------- /

  !/

  RECURSIVE SUBROUTINE qsort(ARRAY,IDX,LO,HI)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |            Stefan Zieger          |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          6-Sep-2016 |

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

    !/

    !/    06-Sep-2016 : Origination, based on code by Mirko ( version 5.16 )

    !/                  Velic (BoM, Australia)

    !/

    !/    Copyright 2009-2013 National Weather Service (NWS),

    !/       National Oceanic and Atmospheric Administration.  All rights

    !/       reserved.  WAVEWATCH III is a trademark of the NWS.

    !/       No unauthorized use without permission.

    !/

    !

    !  1. Purpose :

    !

    !     Quicksort algorithm.

    !

    !  2. Method

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     ARRAY   REAL    ARR  in/out   Input array

    !     IDX     REAL    ARR  in/out   Original indices of input array

    !     LO      INTEGER      input    First element

    !     HI      INTEGER      input    Last element

    !

    IMPLICIT NONE

    !/

    INTEGER, INTENT(IN) :: lo,hi

    REAL,INTENT(INOUT)  :: array(:),idx(:)

    !/

    !     Local variables

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

    LOGICAL :: loop

    INTEGER :: top, bot

    REAL    :: val, tmp

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Scope    Description

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

    !

    !  5. Subroutines calling

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !  9. Switches :

    !

    ! 10. Source code :

    !

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

    !

    !/    --- Check array size and bounds. ---

    IF ( SIZE(array).EQ. 0 ) THEN

      WRITE(6,199)

      CALL abort

    ELSE IF ( SIZE(array).NE.SIZE(idx) ) THEN

      WRITE(6,201)

      CALL abort

    ELSE IF ( lbound(array,1).GT.lo ) THEN

      WRITE(6,203)

      CALL abort

    ELSE IF ( ubound(array,1).LT.hi ) THEN

      WRITE(6,205)

      CALL abort

    END IF

    !

    top = lo

    bot = hi

    val = array(int((lo+hi)/2))

    !

    loop = .true.

    DO WHILE ( loop )

      DO WHILE ( array(top).LT.val )

        top = top + 1

      END DO

      DO WHILE ( val.LT.array(bot) )

        bot = bot - 1

      END DO

      IF ( top.LT.bot ) THEN

        !/        --- Swap values at indices TOP and BOT ---

        tmp = array(top)

        array(top) = array(bot)

        array(bot) = tmp

        !/        --- Swap index values at indices TOP and BOT ---

        tmp = idx(top)

        idx(top) = idx(bot)

        idx(bot) = tmp

        !

        top = top + 1

        bot = bot - 1

      ELSE

        loop = .false.

      END IF


    END DO


    !/    --- Recursive call quicksort ---

    IF (lo.LT.top-1) CALL qsort(array,idx,lo,top-1)

    IF (bot+1.LT.hi) CALL qsort(array,idx,bot+1,hi)

    !

    RETURN

    !/

199 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/            &

         '     QSORT ARRAY IS EMPTY'                           )

201 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/            &

         '     QSORT ARRAY SIZE AND INDEX ARRAY SIZE MISMATCH' )

203 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/            &

         '     QSORT ARRAY INDEX OUT OF LOWER BOUND'           )

205 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/            &

         '     QSORT ARRAY INDEX OUT OF UPPER BOUND'           )

    !/

  END SUBROUTINE qsort

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

  !/

  RECURSIVE SUBROUTINE qsort_desc(ARRAY,IDX,LO,HI)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |            Stefan Zieger          |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :          6-Sep-2016 |

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

    !/

    !/    06-Sep-2016 : Origination, based on code by Mirko ( version 5.16 )

    !/                  Velic (BoM, AUstralia)

    !/

    !/    Copyright 2009-2013 National Weather Service (NWS),

    !/       National Oceanic and Atmospheric Administration.  All rights

    !/       reserved.  WAVEWATCH III is a trademark of the NWS.

    !/       No unauthorized use without permission.

    !/

    !

    !  1. Purpose :

    !

    !     Quicksort algorithm with descending sort order.

    !

    !  2. Method

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     ARRAY   REAL    ARR  in/out   Input array

    !     LO      INTEGER      input    First element

    !     HI      INTEGER      input    Last element

    !

    IMPLICIT NONE

    !/

    INTEGER, INTENT(IN) :: lo,hi

    REAL,INTENT(INOUT)  :: array(:),idx(:)

    !/

    !     Local variables

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

    INTEGER :: top, bot, i

    REAL    :: val, tmp

    LOGICAL :: loop

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Scope    Description

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

    !

    !  5. Subroutines calling

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !  9. Switches :

    !

    ! 10. Source code :

    !

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

    !

    !/    --- Check array size and bounds. ---

    IF ( SIZE(array).EQ. 0 ) THEN

      WRITE(6,199)

      CALL abort

    ELSE IF ( SIZE(array).NE.SIZE(idx) ) THEN

      WRITE(6,201)

      CALL abort

    ELSE IF ( lbound(array,1).GT.lo ) THEN

      WRITE(6,203)

      CALL abort

    ELSE IF ( ubound(array,1).LT.hi ) THEN

      WRITE(6,205)

      CALL abort

    END IF

    !

    top = lo

    bot = hi

    val = array(int((lo+hi)/2))

    !

    loop = .true.

    DO WHILE ( loop )

      DO WHILE ( array(top).GT.val )

        top = top + 1

      END DO

      DO WHILE ( val.GT.array(bot) )

        bot = bot - 1

      END DO

      IF ( top.LT.bot ) THEN

        !/        --- Swap values at indices TOP and BOT ---

        tmp = array(top)

        array(top) = array(bot)

        array(bot) = tmp

        !/        --- Swap index values at indices TOP and BOT ---

        tmp = idx(top)

        idx(top) = idx(bot)

        idx(bot) = tmp

        !

        top = top + 1

        bot = bot - 1

      ELSE

        loop = .false.

      END IF


    END DO

    !/    --- Recursive call quicksort ---

    IF (lo.LT.top-1) CALL qsort_desc(array,idx,lo,top-1)

    IF (bot+1.LT.hi) CALL qsort_desc(array,idx,bot+1,hi)

    !

    RETURN

    !/

199 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/            &

         '     QSORT ARRAY IS EMPTY'                           )

201 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/            &

         '     QSORT ARRAY SIZE AND INDEX ARRAY SIZE MISMATCH' )

203 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/            &

         '     QSORT ARRAY INDEX OUT OF LOWER BOUND'           )

205 FORMAT (/' *** WAVEWATCH III ERROR IN W3SYSTRK : '/            &

         '     QSORT ARRAY INDEX OUT OF UPPER BOUND'           )

    !/

  END SUBROUTINE qsort_desc

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

  !/

  FUNCTION swapi4(INT4) RESULT(INT4SWP)

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           S. Zieger               |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         03-Jan-2017 |

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

    !/

    !/    03-Jan-2017 : Origination                         ( version 5.16 )

    !/                                                        (S. Zieger)

    !/

    !  1. Purpose :

    !

    !     Return a Byte-swapped integer (size of 4 bytes)

    !

    !  2. Source code :

    !

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

    IMPLICIT NONE

    INTEGER(KIND=4), INTENT(IN)  :: int4

    INTEGER(KIND=4)              :: int4swp

    !/

    !     Local variables

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

    INTEGER(KIND=1), DIMENSION(4) :: bytein, byteout

    !/

    bytein  = transfer(int4, bytein)

    byteout = (/bytein(4),bytein(3),bytein(2),bytein(1)/)

    int4swp = transfer(byteout, int4swp)

    !/

    RETURN

    !/

  END FUNCTION swapi4

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

  SUBROUTINE findijv4 (a   ,b   ,maxI,  maxJ ,indA ,indB )

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |     A. J. van der Westhuysen      |

    !/                  |            Jeff Hanson            |

    !/                  |        Eve-Marie Devaliere        |

    !/                  |                        FORTRAN 95 |

    !/                  | Last update :         03-Mar-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 )

    !/    03-Mar-2017 : Calls to INTERSECT and UNION        ( version 5.16 )

    !/                  replaced  (S. Zieger, BoM, Australia)

    !/

    !/    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.

    !/

    IMPLICIT NONE

    !

    !  1. Purpose :

    !

    !     Find a(i,j) indices of system "a" that lie over or along the

    !     fringes of system "b".

    !

    !  2. Method

    !

    !     (i)  Use an index matrix to map locations of wave systems in B

    !     (ii) Avoid multiple use of INTERSECT and UNION as in findIJV3

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !     a, b    Type(system)  input   Final set of tracked systems, for one time level

    !     maxI         Int      input   Number rows indices of wave field

    !     maxJ         Int      input   Number column indices of wave field

    !     indA*, indB* Int.A.   output  Pointer array of indices for combining systems

    !

    TYPE(system) :: a, b

    INTEGER      :: maxI, maxJ

    INTEGER, POINTER :: indA(:), indB(:)

    !

    INTENT (IN) a, b, maxi,maxj

    !

    !     Local variables

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

    !     posB         Int              Neighbour index

    !     posB_MM      Int              Neighbour index (-1,-1)

    !     posB_MP      Int              Neighbour index (-1,+1)

    !     posB_PM      Int              Neighbour index (+1,-1)

    !     posB_PP      Int              Neighbour index (+1,+1)

    !     tmpA*, tmpB* Int.A.           Array of indices for combining

    !     systems

    !

    INTEGER      :: LENG_AI,LENG_BI

    INTEGER      :: OUTA,OUTB,I,J,IND,OUTDUMB

    INTEGER      :: POSB,POSB_MM,POSB_PM,POSB_MP,POSB_PP

    INTEGER      :: IND_B2(maxI,maxJ)

    REAL,ALLOCATABLE :: TMPA(:),DUMA(:),TMPB(:)

    REAL,POINTER     :: DUMB(:)

    LOGICAL      :: FOUND

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Scope    Description

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

    !      QSORT     Subr. Private  Quicksort algorithm

    !      UNIQUE    Subr. Private  Return sorted unique numbers of an array

    !

    !  5. Subroutines calling

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !  9. Switches :

    !

    ! 10. Source code :

    !

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

    !

    NULLIFY(dumb)

    !

    IF (ASSOCIATED(inda)) DEALLOCATE(inda)

    IF (ASSOCIATED(indb)) DEALLOCATE(indb)

    !

    leng_ai = length(real(a%i),SIZE(a%i),real(9999))

    leng_bi = length(real(b%i),SIZE(b%i),real(9999))

    !

    ALLOCATE(tmpa(leng_ai))

    ALLOCATE(tmpb(5*leng_ai))

    !

    tmpa(:) = 9999.

    tmpb(:) = 9999.

    !

    outa = 0

    outb = 0

    ind_b2(:,:) = 0

    !

    DO ind=1,leng_bi

      i = b%I(ind)

      j = b%J(ind)

      IF (ind_b2(i,j).EQ.0) ind_b2(i,j) = ind

    END DO

    !

    DO ind=1,leng_ai

      i = a%I(ind)

      j = a%J(ind)

      posb = ind_b2(i,j)

      posb_mm = 0

      posb_pp = 0

      posb_mp = 0

      posb_pm = 0

      IF (i.GT.1.AND.j.GT.1)       posb_mm = ind_b2(i-1,j-1)

      IF (i.GT.1.AND.j.LT.maxj)    posb_mp = ind_b2(i-1,j+1)

      IF (i.LT.maxi.AND.j.LT.maxj) posb_pp = ind_b2(i+1,j+1)

      IF (i.LT.maxi.AND.j.GT.1)    posb_pm = ind_b2(i+1,j-1)


      found = .false.

      IF (posb.NE.0) THEN

        outb = outb + 1

        tmpb(outb) = real(posb)

        IF (.NOT.found) THEN

          outa = outa + 1

          tmpa(outa) = real(ind)

          found = .true.

        END IF

      END IF

      IF (posb_mm.NE.0) THEN

        outb = outb + 1

        tmpb(outb) = real(posb_mm)

        IF (.NOT.found) THEN

          outa = outa + 1

          tmpa(outa) = real(ind)

          found = .true.

        END IF

      END IF

      IF (posb_mp.NE.0) THEN

        outb = outb + 1

        tmpb(outb) = real(posb_mp)

        IF (.NOT.found) THEN

          outa = outa + 1

          tmpa(outa) = real(ind)

          found = .true.

        END IF

      END IF

      IF (posb_pm.NE.0) THEN

        outb = outb + 1

        tmpb(outb) = real(posb_pm)

        IF (.NOT.found) THEN

          outa = outa + 1

          tmpa(outa) = real(ind)

          found = .true.

        END IF

      END IF

      IF (posb_pp.NE.0) THEN

        outb = outb + 1

        tmpb(outb) = real(posb_pp)

        IF (.NOT.found) THEN

          outa = outa + 1

          tmpa(outa) = real(ind)

          found = .true.

        END IF

      END IF


    END DO

    !

    !/    Compact indices for wave systems in B.

    !/    Check for empty arrays first.

    IF (outb.GT.0) THEN

      CALL unique(tmpb,outb,dumb,outdumb)

      outb = outdumb

    END IF

    ALLOCATE(indb(outb))

    IF (outb.GT.0) indb(1:outb) = int(dumb(1:outb))

    IF (ASSOCIATED(dumb)) DEALLOCATE(dumb)

    !

    !/    Allocate output array and transfer content

    !/    for wave systems in A.

    ALLOCATE(inda(outa))

    IF (outa.GT.0) THEN

      ALLOCATE(duma(outa))

      duma(:) = 0

      CALL qsort(tmpa(1:outa),duma(1:outa),1,outa)

      IF (ALLOCATED(duma)) DEALLOCATE(duma)

      inda(1:outa) = int(tmpa(1:outa))

    END IF

    !/

    IF (ALLOCATED(tmpa)) DEALLOCATE(tmpa)

    IF (ALLOCATED(tmpb)) DEALLOCATE(tmpb)

    !/

    RETURN

    !/

  END SUBROUTINE findijv4

  !/ End of findIJV4 --------------------------------------------------- /

  !/

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

  !/

END MODULE w3strkmd

!/

!/ End of module W3STRKMD -------------------------------------------- /

!/