w3strkmd Module Reference

Data Types
type	dat2d
type	mtchsys
type	neighbr
type	param
type	sysmemory
type	system
type	timsys
type	wind

Functions/Subroutines
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)
subroutine	spiraltrackv3 (wsdat, dirKnob, perKnob, wetPts, hsKnob, seedLat, seedLon, maxSys, sys)
subroutine	timetrackingv2 (sysA, maxSys, tpTimeKnob, dirTimeKnob, ts0, maxGroup, dt, lonext, latext, maxI, maxJ)
subroutine	findway (way, horizStepCount, vertStepCount, vertBorder, horizBorder, stepCount)
subroutine	findnext (i, j, maxI, maxJ, way, vertBorder, horizBorder)
subroutine	findsys (i, j, wsdat, maxSys, ngbrExt, maxI, maxJ, perKnob, dirKnob, hsKnob)
subroutine	combinewavesystems (wsdat, maxSys, maxPts, maxI, maxJ, perKnob, dirKnob, hsKnob, combine, sys)
subroutine	printfinalsys (wsdat, maxSys, actSysInd, maxI, maxJ, flag, sys)
subroutine	combinesys (wsdat, sys, maxSys, maxI, maxJ, actSysInd, perKnob, dirKnob)
subroutine	combinepartitionsv2 (dat)
real function	mean_anglev2 (ang, ll)
real function	mean_anglev3 (ang, hsign, ll)
subroutine	unique (INARRAY, INSIZE, OUTARRAY, OUTSIZE)
subroutine	sort (INARRAY, INSIZE, OUTARRAY, IY, DIRECTION)
subroutine	setdiff (INARRAY1, INSIZE1, INARRAY2, INSIZE2, OUTARRAY, OUTSIZE)
subroutine	intersect (INARRAY1, INSIZE1, INARRAY2, INSIZE2, OUTARRAY, OUTSIZE, IND1, IND2)
subroutine	union (INARRAY1, INSIZE1, INARRAY2, INSIZE2, OUTARRAY, OUTSIZE)
integer function	length (ARRAY, ARRSIZE, VAL)
integer function	findfirst (ARRAY, ARRSIZE, VAL)
real function	std (ARRAY, N)
recursive subroutine	qsort (ARRAY, IDX, LO, HI)
recursive subroutine	qsort_desc (ARRAY, IDX, LO, HI)
integer(kind=4) function	swapi4 (INT4)
subroutine	findijv4 (a, b, maxI, maxJ, indA, indB)

Function/Subroutine Documentation

combinepartitionsv2()

subroutine w3strkmd::combinepartitionsv2

(

type(param), intent(inout)

dat

)

Definition at line 4282 of file w3strkmd.F90.

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

References length().

Referenced by combinesys(), and findsys().

combinesys()

subroutine w3strkmd::combinesys	(	type(dat2d)	wsdat,
		type(system), dimension(:), pointer	sys,
		integer	maxSys,
		integer, intent(in)	maxI,
		integer, intent(in)	maxJ,
		integer, dimension(:), allocatable	actSysInd,
		real, intent(in)	perKnob,
		real, intent(in)	dirKnob
	)

Definition at line 3683 of file w3strkmd.F90.

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

References combinepartitionsv2(), file(), findfirst(), findijv4(), length(), mean_anglev2(), mean_anglev3(), setdiff(), sort(), union(), and unique().

Referenced by combinewavesystems().

combinewavesystems()

subroutine w3strkmd::combinewavesystems	(	type(dat2d), intent(inout)	wsdat,
		integer, intent(inout)	maxSys,
		integer, intent(in)	maxPts,
		integer, intent(in)	maxI,
		integer, intent(in)	maxJ,
		real	perKnob,
		real	dirKnob,
		real, intent(in)	hsKnob,
		integer, intent(in)	combine,
		type(system), dimension(:), pointer	sys
	)

Definition at line 3209 of file w3strkmd.F90.

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

References combinesys(), length(), printfinalsys(), and std().

Referenced by spiraltrackv3().

findfirst()

integer function w3strkmd::findfirst	(	real, dimension(arrsize)	ARRAY,
		integer	ARRSIZE,
		real	VAL
	)

Definition at line 5495 of file w3strkmd.F90.

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

Referenced by combinesys(), findsys(), and timetrackingv2().

findijv4()

subroutine w3strkmd::findijv4	(	type(system), intent(in)	a,
		type(system), intent(in)	b,
		integer, intent(in)	maxI,
		integer, intent(in)	maxJ,
		integer, dimension(:), pointer	indA,
		integer, dimension(:), pointer	indB
	)

Definition at line 5893 of file w3strkmd.F90.

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

References length(), qsort(), and unique().

Referenced by combinesys().

findnext()

subroutine w3strkmd::findnext	(	integer, intent(inout)	i,
		integer, intent(inout)	j,
		integer, intent(in)	maxI,
		integer, intent(in)	maxJ,
		character, intent(in)	way,
		integer, intent(out)	vertBorder,
		integer, intent(out)	horizBorder
	)

Definition at line 2794 of file w3strkmd.F90.

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

Referenced by spiraltrackv3().

findsys()

subroutine w3strkmd::findsys	(	integer, intent(in)	i,
		integer, intent(in)	j,
		type(dat2d), intent(inout)	wsdat,
		integer, intent(inout)	maxSys,
		integer, intent(in)	ngbrExt,
		integer, intent(in)	maxI,
		integer, intent(in)	maxJ,
		real, intent(in)	perKnob,
		real, intent(in)	dirKnob,
		real	hsKnob
	)

Definition at line 2908 of file w3strkmd.F90.

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

References combinepartitionsv2(), findfirst(), length(), mean_anglev2(), and unique().

Referenced by spiraltrackv3().

findway()

subroutine w3strkmd::findway	(	character, intent(inout)	way,
		integer	horizStepCount,
		integer	vertStepCount,
		integer, intent(in)	vertBorder,
		integer, intent(in)	horizBorder,
		integer, intent(out)	stepCount
	)

Definition at line 2678 of file w3strkmd.F90.

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

Referenced by spiraltrackv3().

intersect()

subroutine w3strkmd::intersect	(	real, dimension(insize1), intent(in)	INARRAY1,
		integer, intent(in)	INSIZE1,
		real, dimension(insize2), intent(in)	INARRAY2,
		integer, intent(in)	INSIZE2,
		real, dimension(:), pointer	OUTARRAY,
		integer, intent(out)	OUTSIZE,
		real, dimension(:), pointer	IND1,
		real, dimension(:), pointer	IND2
	)

Definition at line 5096 of file w3strkmd.F90.

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

Referenced by timetrackingv2().

length()

integer function w3strkmd::length	(	real, dimension(arrsize)	ARRAY,
		integer	ARRSIZE,
		real	VAL
	)

Definition at line 5432 of file w3strkmd.F90.

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

Referenced by combinepartitionsv2(), combinesys(), combinewavesystems(), findijv4(), findsys(), printfinalsys(), spiraltrackv3(), and timetrackingv2().

mean_anglev2()

real function w3strkmd::mean_anglev2	(	real, dimension(ll)	ang,
		integer	ll
	)

Definition at line 4463 of file w3strkmd.F90.

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

References constants::pi.

Referenced by combinesys(), findsys(), and printfinalsys().

mean_anglev3()

real function w3strkmd::mean_anglev3	(	real, dimension(ll)	ang,
		real, dimension(ll)	hsign,
		integer	ll
	)

Definition at line 4570 of file w3strkmd.F90.

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

References constants::pi.

Referenced by combinesys(), and printfinalsys().

printfinalsys()

subroutine w3strkmd::printfinalsys	(	type(dat2d), intent(in)	wsdat,
		integer, intent(out)	maxSys,
		integer, dimension(:), intent(in), allocatable	actSysInd,
		integer, intent(in)	maxI,
		integer, intent(in)	maxJ,
		integer, intent(in)	flag,
		type(system), dimension(:), pointer	sys
	)

Definition at line 3455 of file w3strkmd.F90.

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

References length(), mean_anglev2(), mean_anglev3(), setdiff(), sort(), and unique().

Referenced by combinewavesystems().

qsort()

recursive subroutine w3strkmd::qsort	(	real, dimension(:), intent(inout)	ARRAY,
		real, dimension(:), intent(inout)	IDX,
		integer, intent(in)	LO,
		integer, intent(in)	HI
	)

Definition at line 5608 of file w3strkmd.F90.

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

Referenced by findijv4(), setdiff(), sort(), union(), and unique().

qsort_desc()

recursive subroutine w3strkmd::qsort_desc	(	real, dimension(:), intent(inout)	ARRAY,
		real, dimension(:), intent(inout)	IDX,
		integer, intent(in)	LO,
		integer, intent(in)	HI
	)

Definition at line 5734 of file w3strkmd.F90.

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

Referenced by sort().

setdiff()

subroutine w3strkmd::setdiff	(	real, dimension(insize1), intent(in)	INARRAY1,
		integer, intent(in)	INSIZE1,
		real, dimension(insize2), intent(in)	INARRAY2,
		integer, intent(in)	INSIZE2,
		real, dimension(:), pointer	OUTARRAY,
		integer, intent(out)	OUTSIZE
	)

Definition at line 4937 of file w3strkmd.F90.

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

References qsort(), and unique().

Referenced by combinesys(), and printfinalsys().

sort()

subroutine w3strkmd::sort	(	real, dimension(insize), intent(in)	INARRAY,
		integer, intent(in)	INSIZE,
		real, dimension(insize), intent(out)	OUTARRAY,
		integer, dimension(insize), intent(out)	IY,
		character, intent(in)	DIRECTION
	)

Definition at line 4817 of file w3strkmd.F90.

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

References qsort(), and qsort_desc().

Referenced by combinesys(), printfinalsys(), and timetrackingv2().

spiraltrackv3()

subroutine w3strkmd::spiraltrackv3	(	type(dat2d), intent(inout)	wsdat,
		real, intent(in)	dirKnob,
		real, intent(in)	perKnob,
		real, intent(in)	wetPts,
		real, intent(in)	hsKnob,
		real, intent(in)	seedLat,
		real, intent(in)	seedLon,
		integer	maxSys,
		type(system), dimension(:), pointer	sys
	)

Definition at line 1749 of file w3strkmd.F90.

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

References combinewavesystems(), findnext(), findsys(), findway(), and length().

Referenced by wavetracking_nws_v2().

std()

real function w3strkmd::std	(	real, dimension(n)	ARRAY,
		integer	N
	)

Definition at line 5554 of file w3strkmd.F90.

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

Referenced by combinewavesystems().

swapi4()

integer(kind=4) function w3strkmd::swapi4

(

integer(kind=4), intent(in)

INT4

)

Definition at line 5858 of file w3strkmd.F90.

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

Referenced by wavetracking_nws_v2().

timetrackingv2()

subroutine w3strkmd::timetrackingv2	(	type(timsys), dimension(:), pointer	sysA,
		integer, dimension(:), pointer	maxSys,
		real, intent(in)	tpTimeKnob,
		real, intent(in)	dirTimeKnob,
		integer, intent(in)	ts0,
		integer, intent(out)	maxGroup,
		real	dt,
		real	lonext,
		real	latext,
		integer, intent(in)	maxI,
		integer, intent(in)	maxJ
	)

Definition at line 1968 of file w3strkmd.F90.

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

References file(), findfirst(), intersect(), length(), and sort().

Referenced by wavetracking_nws_v2().

union()

subroutine w3strkmd::union	(	real, dimension(insize1), intent(in)	INARRAY1,
		integer, intent(in)	INSIZE1,
		real, dimension(insize2), intent(in)	INARRAY2,
		integer, intent(in)	INSIZE2,
		real, dimension(:), pointer	OUTARRAY,
		integer, intent(out)	OUTSIZE
	)

Definition at line 5256 of file w3strkmd.F90.

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

References qsort().

Referenced by combinesys().

unique()

subroutine w3strkmd::unique	(	real, dimension(insize), intent(in)	INARRAY,
		integer, intent(in)	INSIZE,
		real, dimension(:), pointer	OUTARRAY,
		integer, intent(out)	OUTSIZE
	)

Definition at line 4686 of file w3strkmd.F90.

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

References qsort().

Referenced by combinesys(), findijv4(), findsys(), printfinalsys(), setdiff(), and wavetracking_nws_v2().

wavetracking_nws_v2()

subroutine w3strkmd::wavetracking_nws_v2	(	integer, intent(in)	intype,
		integer, intent(in)	tmax,
		integer, intent(in)	tcur,
		character, intent(in)	filename,
		real*8	tstart,
		real*8	tend,
		real	dt,
		integer	ntint,
		real, intent(in)	minlon,
		real, intent(in)	maxlon,
		real, intent(in)	minlat,
		real, intent(in)	maxlat,
		integer	mxcwt,
		integer	mycwt,
		real, intent(in)	dirKnob,
		real, intent(in)	perKnob,
		real, intent(in)	hsKnob,
		real, intent(in)	wetPts,
		real, intent(in)	seedLat,
		real, intent(in)	seedLon,
		real, intent(in)	dirTimeKnob,
		real, intent(in)	tpTimeKnob,
		character, intent(in)	paramFile,
		type(timsys), dimension(:), pointer	sysA,
		type(dat2d), dimension(:), pointer	wsdat,
		integer, dimension(:), pointer	maxSys,
		integer, intent(out)	maxGroup
	)

042916 IF ( (llat(l).EQ.mlat(i,j)).AND. & 042916 (llon(l).EQ.mlon(i,j)) ) THEN

Definition at line 298 of file w3strkmd.F90.

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

References file(), constants::file_endian, include(), spiraltrackv3(), swapi4(), timetrackingv2(), and unique().

Referenced by ww3_systrk().