w3partmd Module Reference

Spectral partitioning according to the watershed method. More...

Functions/Subroutines
subroutine	w3part (SPEC, UABS, UDIR, DEPTH, WN, NP, XP, DIMXP)
	Interface to watershed partitioning routines. More...
subroutine	ptsort (IMI, IND, IHMAX)
	Sorts the image data in ascending order. More...
subroutine	ptnghb
	Nearest neighbour calculation. More...
subroutine	pt_fld (IMI, IND, IMO, ZP, NPART)
	Image watersheding. More...
subroutine	ptmean (NPI, IMO, ZP, DEPTH, UABS, UDIR, WN, NPO, XP, DIMXP, PMAP)
	Compute mean parameters per partition. More...

Detailed Description

Spectral partitioning according to the watershed method.

Multiple partitioning methods are provided in this module that can be selected via the PTMETH namelist variable. Please see the w3part() subroutine for details

Author

Barbara Tracey, H. L. Tolman, M. Szyszka, Chris Bunney

Date

23 Jul 2018

Function/Subroutine Documentation

pt_fld()

subroutine w3partmd::pt_fld	(	integer, dimension(nspec), intent(in)	IMI,
		integer, dimension(nspec), intent(in)	IND,
		integer, dimension(nspec), intent(out)	IMO,
		real, dimension(nspec), intent(in)	ZP,
		integer, intent(out)	NPART
	)

Image watersheding.

This subroutine does incremental flooding of the image to determine the watershed image.

Parameters

[in]	IMI	Input discretized spectrum
[in]	IND	Sorted addresses
[out]	IMO	Output partitioned spectrum
[in]	ZP	Spectral array
[out]	NPART	Number of partitions found

Author

H.L. Tolman

Date

01 Nov 2006

Definition at line 841 of file w3partmd.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         01-Nov-2006 !
    !/                  +-----------------------------------+
    !/
    !/    01-Nov-2006 : Origination.                       ( version 3.10 )
    !/
    !  1. Purpose :
    !
    !     This subroutine does incremental flooding of the image to
    !     determine the watershed image.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       IMI     I.A.   I   Input discretized spectrum.
    !       IND     I.A.   I   Sorted addresses.
    !       IMO     I.A.   O   Output partitioned spectrum.
    !       ZP      R.A.   I   Spectral array.
    !       NPART   Int.   O   Number of partitions found.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Sur.  W3SERVMD Subroutine tracing.
    !     ----------------------------------------------------------------
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !
    USE w3gdatmd, ONLY: nspec
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: IMI(NSPEC), IND(NSPEC)
    INTEGER, INTENT(OUT)    :: IMO(NSPEC), NPART
    REAL, INTENT(IN)        :: ZP(NSPEC)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: MASK, INIT, IWSHED, IMD(NSPEC),      &
         IC_LABEL, IFICT_PIXEL, M, IH, MSAVE, &
         IP, I, IPP, IC_DIST, IEMPTY, IPPP,   &
         JL, JN, IPT, J
    INTEGER                 :: IQ(NSPEC), IQ_START, IQ_END
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: ZPMAX, EP1, DIFF
    !/
#ifdef W3_S
    CALL strace (ient, 'PT_FLD')
#endif
    !
    ! -------------------------------------------------------------------- /
    ! 0.  Initializations
    !
    mask        = -2
    init        = -1
    iwshed      =  0
    imo         = init
    ic_label    =  0
    imd         =  0
    ifict_pixel = -100
    !
    iq_start    =  1
    iq_end      =  1
    !
    zpmax       = maxval( zp )
    !
    ! -------------------------------------------------------------------- /
    ! 1.  Loop over levels
    !
    m      =  1
    !
    DO ih=1, ihmax
      msave  = m
      !
      ! 1.a Pixels at level IH
      !
      DO
        ip     = ind(m)
        IF ( imi(ip) .NE. ih ) EXIT
        !
        !     Flag the point, if it stays flagge, it is a separate minimum.
        !
        imo(ip) = mask
        !
        !     Consider neighbors. If there is neighbor, set distance and add
        !     to queue.
        !
        DO i=1, neigh(9,ip)
          ipp    = neigh(i,ip)
          IF ( (imo(ipp).GT.0) .OR. (imo(ipp).EQ.iwshed) ) THEN
            imd(ip) = 1
            CALL fifo_add (ip)
            EXIT
          END IF
        END DO
        !
        IF ( m+1 .GT. nspec ) THEN
          EXIT
        ELSE
          m = m + 1
        END IF
        !
      END DO
      !
      ! 1.b Process the queue
      !
      ic_dist = 1
      CALL fifo_add (ifict_pixel)
      !
      DO
        CALL fifo_first (ip)
        !
        !     Check for end of processing
        !
        IF ( ip .EQ. ifict_pixel ) THEN
          CALL fifo_empty (iempty)
          IF ( iempty .EQ. 1 ) THEN
            EXIT
          ELSE
            CALL fifo_add (ifict_pixel)
            ic_dist = ic_dist + 1
            CALL fifo_first (ip)
          END IF
        END IF
        !
        !     Process queue
        !
        DO i=1, neigh(9,ip)
          ipp = neigh(i,ip)
          !
          !     Check for labeled watersheds or basins
          !
          IF ( (imd(ipp).LT.ic_dist) .AND. ( (imo(ipp).GT.0) .OR.  &
               (imo(ipp).EQ.iwshed))) THEN
            !
            IF ( imo(ipp) .GT. 0 ) THEN
              !
              IF ((imo(ip) .EQ. mask) .OR. (imo(ip) .EQ. &
                   iwshed)) THEN
                imo(ip) = imo(ipp)
              ELSE IF (imo(ip) .NE. imo(ipp)) THEN
                imo(ip) = iwshed
              END IF
              !
            ELSE IF (imo(ip) .EQ. mask) THEN
              !
              imo(ip) = iwshed
              !
            END IF
            !
          ELSE IF ( (imo(ipp).EQ.mask) .AND. (imd(ipp).EQ.0) ) THEN
            !
            imd(ipp) = ic_dist + 1
            CALL fifo_add (ipp)
            !
          END IF
          !
        END DO
        !
      END DO
      !
      ! 1.c Check for mask values in IMO to identify new basins
      !
      m = msave
      !
      DO
        ip     = ind(m)
        IF ( imi(ip) .NE. ih ) EXIT
        imd(ip) = 0
        !
        IF (imo(ip) .EQ. mask) THEN
          !
          ! ... New label for pixel
          !
          ic_label = ic_label + 1
          CALL fifo_add (ip)
          imo(ip) = ic_label
          !
          ! ... and all connected to it ...
          !
          DO
            CALL fifo_empty (iempty)
            IF ( iempty .EQ. 1 ) EXIT
            CALL fifo_first (ipp)
            !
            DO i=1, neigh(9,ipp)
              ippp   = neigh(i,ipp)
              IF ( imo(ippp) .EQ. mask ) THEN
                CALL fifo_add (ippp)
                imo(ippp) = ic_label
              END IF
            END DO
            !
          END DO
          !
        END IF
        !
        IF ( m + 1 .GT. nspec ) THEN
          EXIT
        ELSE
          m = m + 1
        END IF
        !
      END DO
      !
    END DO
    !
    ! -------------------------------------------------------------------- /
    ! 2.  Find nearest neighbor of 0 watershed points and replace
    !     use original input to check which group to affiliate with 0
    !     Soring changes first in IMD to assure symetry in adjustment.
    !
    DO j=1, 5
      imd    = imo
      DO jl=1 , nspec
        ipt    = -1
        IF ( imo(jl) .EQ. 0 ) THEN
          ep1    = zpmax
          DO jn=1, neigh(9,jl)
            diff   = abs( zp(jl) - zp(neigh(jn,jl)))
            IF ( (diff.LE.ep1) .AND. (imo(neigh(jn,jl)).NE.0) ) THEN
              ep1    = diff
              ipt    = jn
            END IF
          END DO
          IF ( ipt .GT. 0 ) imd(jl) = imo(neigh(ipt,jl))
        END IF
      END DO
      imo    = imd
      IF ( minval(imo) .GT. 0 ) EXIT
    END DO
    !
    npart = ic_label
    !
    RETURN
    !
  CONTAINS
    !/ ------------------------------------------------------------------- /
    SUBROUTINE fifo_add ( IV )
 
      INTEGER, INTENT(IN)      :: IV
      !
      iq(iq_end) = iv
      !
      iq_end = iq_end + 1
      IF ( iq_end .GT. nspec ) iq_end = 1
      !
      RETURN
    END SUBROUTINE fifo_add
    !/ ------------------------------------------------------------------- /
    SUBROUTINE fifo_empty ( IEMPTY )
 
      INTEGER, INTENT(OUT)     :: IEMPTY
      !
      IF ( iq_start .NE. iq_end ) THEN
        iempty = 0
      ELSE
        iempty = 1
      END IF
      !
      RETURN
    END SUBROUTINE fifo_empty
    !/ ------------------------------------------------------------------- /
    SUBROUTINE fifo_first ( IV )
 
      INTEGER, INTENT(OUT)     :: IV
      !
      iv = iq(iq_start)
      !
      iq_start = iq_start + 1
      IF ( iq_start .GT. nspec ) iq_start = 1
      !
      RETURN
    END SUBROUTINE fifo_first
    !/
    !/ End of PT_FLD ----------------------------------------------------- /
    !/

References fifo_add(), fifo_empty(), fifo_first(), w3odatmd::ihmax, w3gdatmd::nspec, and w3servmd::strace().

Referenced by w3part().

ptmean()

subroutine w3partmd::ptmean	(	integer, intent(in)	NPI,
		integer, dimension(nspec), intent(in)	IMO,
		real, dimension(nspec), intent(in)	ZP,
		real, intent(in)	DEPTH,
		real, intent(in)	UABS,
		real, intent(in)	UDIR,
		real, dimension(nk), intent(in)	WN,
		integer, intent(out)	NPO,
		real, dimension(dimp,0:dimxp), intent(out)	XP,
		integer, intent(in)	DIMXP,
		integer, dimension(dimxp), intent(out)	PMAP
	)

Compute mean parameters per partition.

Parameters

[in]	NPI	Number of partitions found.
[in]	IMO	Partition map.
[in]	ZP	Input spectrum.
[in]	DEPTH	Water depth.
[in]	UABS	Wind speed.
[in]	UDIR	Wind direction.
[in]	WN	Wavenumebers for each frequency.
[out]	NPO	Number of partitions with mean parameters.
[out]	XP	Array with output parameters.
[in]	DIMXP	Second dimension of XP.
[out]	PMAP	Mapping between orig. and combined partitions

Author

Barbara Tracy, H. L. Tolman, M. Szyszka, C. Bunney

Date

02 Dec 2010

Definition at line 1178 of file w3partmd.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III          USACE/NOAA |
    !/                  |          Barbara  Tracy           |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         02-Dec-2010 !
    !/                  +-----------------------------------+
    !/
    !/    28-Oct-2006 : Origination.                       ( version 3.10 )
    !/    02-Nov-2006 : Adding tail to integration.        ( version 3.10 )
    !/    24-Mar-2007 : Adding overall field.              ( version 3.11 )
    !/    02-Dec-2010 : Adding a mapping PMAP between      ( version 3.14 )
    !/                  original and combined partitions
    !/                  ( M. Szyszka )
    !/
    !  1. Purpose :
    !
    !     Compute mean parameters per partition.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       NPI     Int.   I   Number of partitions found.
    !       IMO     I.A.   I   Partition map.
    !       ZP      R.A.   I   Input spectrum.
    !       DEPTH   Real   I   Water depth.
    !       UABS    Real   I   Wind speed.
    !       UDIR    Real   I   Wind direction.
    !       WN      R.A.   I   Wavenumebers for each frequency.
    !       NPO     Int.   O   Number of partitions with mean parameters.
    !       XP      R.A.   O   Array with output parameters.
    !       DIMXP   int.   I   Second dimension of XP.
    !       PMAP    I.A.   O   Mapping between orig. and combined partitions
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Sur.  W3SERVMD Subroutine tracing.
    !      WAVNU1    Subr. W3DISPMD Wavenumber computation.
    !     ----------------------------------------------------------------
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !
    USE constants
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    USE w3dispmd, ONLY: wavnu1
    !
    USE w3gdatmd, ONLY: nk, nth, nspec, dth, sig, dsii, dsip,       &
         ecos, esin, xfr, fachfe, th, fte
    USE w3odatmd, ONLY: iaproc, naperr, ndse, ndst
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: NPI, IMO(NSPEC), DIMXP
    INTEGER, INTENT(OUT)    :: NPO, PMAP(DIMXP)
    REAL, INTENT(IN)        :: ZP(NSPEC), DEPTH, UABS, UDIR, WN(NK)
    REAL, INTENT(OUT)       :: XP(DIMP,0:DIMXP)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IK, ITH, ISP, IP, IFPMAX(0:NPI)
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: SUMF(0:NK+1,0:NPI), SUMFW(NK,0:NPI), &
         SUMFX(NK,0:NPI), SUMFY(NK,0:NPI),    &
         SUME(0:NPI), SUMEW(0:NPI),           &
         SUMEX(0:NPI), SUMEY(0:NPI),          &
         EFPMAX(0:NPI), FCDIR(NTH)
    REAL,DIMENSION(0:NPI)   :: SUME1, SUME2, SUMEM1, SUMQP
    REAL                    :: HS, XL, XH, XL2, XH2, EL, EH, DENOM, &
         SIGP, WNP, CGP, UPAR, C(NK), RD, FACT
    REAL                    :: QP, M0, M1, M2, MM1, FSPRD, EPM_FP, ALP_PM
    REAL                    :: Y, YHAT, XHAT, SUMXY, SUMYLOGY, NUMER,&
         SUMY, SUMXXY, SUMXYLOGY, SUMEXP, SUMEYP
    REAL                    :: FTEII
    !/
#ifdef W3_S
    CALL strace (ient, 'PTMEAN')
#endif
    !
    ! -------------------------------------------------------------------- /
    ! 1.  Check on need of processing
    !
    npo    = 0
    xp     = 0.
    !
    IF ( npi .EQ. 0 ) RETURN
    !
    ! -------------------------------------------------------------------- /
    ! 2.  Initialize arrays
    !
    sumf   = 0.
    sumfw  = 0.
    sumfx  = 0.
    sumfy  = 0.
    sume   = 0.
    !
    sume1  = 0.  !/ first spectral moment
    sume2  = 0.  !/ second spectral moment
    sumem1 = 0.  !/ inverse spectral moment
    sumqp  = 0.  !/ peakedness parameter
    !
    sumew  = 0.
    sumex  = 0.
    sumey  = 0.
    ifpmax = 0
    efpmax = 0.
    !
    DO ik=1, nk
      c(ik)  = sig(ik) / wn(ik)
    END DO
    !
    DO ith=1, nth
      upar   = wsmult * uabs * max(0.,cos(th(ith)-dera*udir))
      IF ( upar .LT. c(nk) ) THEN
        fcdir(ith) = sig(nk+1)
      ELSE
        DO ik=nk-1, 2, -1
          IF ( upar .LT. c(ik) ) EXIT
        END DO
        rd     = (c(ik)-upar) / (c(ik)-c(ik+1))
        IF ( rd .LT. 0 ) THEN
          ik     = 0
          rd     = max( 0., rd+1. )
        END IF
        fcdir(ith) = rd*sig(ik+1) + (1.-rd)*sig(ik)
      END IF
    END DO
    !
    ! -------------------------------------------------------------------- /
    ! 3.  Spectral integrals and preps
    ! 3.a Integrals
    !     NOTE: Factor DTH only used in Hs computation.
    !
    DO ik=1, nk
      DO ith=1, nth
        isp    = ik + (ith-1)*nk
        ip     = imo(isp)
        fact   = max( 0. , min( 1. ,                              &
             1. - ( fcdir(ith) - 0.5*(sig(ik-1)+sig(ik)) ) / dsip(ik) ) )
        sumf(ik, 0) = sumf(ik, 0) + zp(isp)
        sumfw(ik, 0) = sumfw(ik, 0) + zp(isp) * fact
        sumfx(ik, 0) = sumfx(ik, 0) + zp(isp) * ecos(ith)
        sumfy(ik, 0) = sumfy(ik, 0) + zp(isp) * esin(ith)
        IF ( ip .EQ. 0 ) cycle
        sumf(ik,ip) = sumf(ik,ip) + zp(isp)
        sumfw(ik,ip) = sumfw(ik,ip) + zp(isp) * fact
        sumfx(ik,ip) = sumfx(ik,ip) + zp(isp) * ecos(ith)
        sumfy(ik,ip) = sumfy(ik,ip) + zp(isp) * esin(ith)
      END DO
    END DO
    sumf(nk+1,:) = sumf(nk,:) * fachfe
    !
    DO ip=0, npi
      DO ik=1, nk
        sume(ip) = sume(ip) + sumf(ik,ip) * dsii(ik)
        sumqp(ip) = sumqp(ip) + sumf(ik,ip)**2 * dsii(ik) * sig(ik)
        sume1(ip) = sume1(ip) + sumf(ik,ip) * dsii(ik) * sig(ik)
        sume2(ip) = sume2(ip) + sumf(ik,ip) * dsii(ik) * sig(ik)**2
        sumem1(ip) = sumem1(ip) + sumf(ik,ip) * dsii(ik) / sig(ik)
 
        sumew(ip) = sumew(ip) + sumfw(ik,ip) * dsii(ik)
        sumex(ip) = sumex(ip) + sumfx(ik,ip) * dsii(ik)
        sumey(ip) = sumey(ip) + sumfy(ik,ip) * dsii(ik)
        IF ( sumf(ik,ip) .GT. efpmax(ip) ) THEN
          ifpmax(ip) = ik
          efpmax(ip) = sumf(ik,ip)
        END IF
      END DO
 
      !SUME (IP) = SUME (IP) + SUMF (NK,IP) * FTE
      !SUME1(IP) = SUME1(IP) + SUMF (NK,IP) * FTE
      !SUME2(IP) = SUME2(IP) + SUMF (NK,IP) * FTE
      !SUMEM1(IP) = SUMEM1(IP) + SUMF (NK,IP) * FTE
      !SUMQP(IP) = SUMQP(IP) + SUMF (NK,IP) * FTE
      !SUMEW(IP) = SUMEW(IP) + SUMFW(NK,IP) * FTE
      !SUMEX(IP) = SUMEX(IP) + SUMFX(NK,IP) * FTE
      !SUMEY(IP) = SUMEY(IP) + SUMFY(NK,IP) * FTE
      ! Met Office: Proposed bugfix for tail calculations, previously
      !  PT1 and PT2 values were found to be too low when using the
      !  FTE scaling factor for the tail. I think there are two issues:
      !  1. energy spectrum is scaled in radian frequency space above by DSII.
      !     This needs to be consistent and FTE contains a DTH*SIG(NK)
      !     factor that is not used in the DSII scaled calcs above
      !  2. the tail fit calcs for period parameters needs to follow
      !     the form used in w3iogomd and scaling should be
      !     based on the relationship between FTE and FT1, FTTR etc.
      !     as per w3iogomd and ww3_grid
      fteii = fte / (dth * sig(nk))
      sume(ip) = sume(ip) + sumf(nk,ip) * fteii
      sume1(ip) = sume1(ip) + sumf(nk,ip) * sig(nk) * fteii * (0.3333 / 0.25)
      sume2(ip) = sume2(ip) + sumf(nk,ip) * sig(nk)**2 * fteii * (0.5 / 0.25)
      sumem1(ip) = sumem1(ip) + sumf(nk,ip) / sig(nk) * fteii * (0.2 / 0.25)
      sumqp(ip) = sumqp(ip) + sumf(nk,ip) * fteii
      sumew(ip) = sumew(ip) + sumfw(nk,ip) * fteii
      sumex(ip) = sumex(ip) + sumfx(nk,ip) * fteii
      sumey(ip) = sumey(ip) + sumfy(nk,ip) * fteii
 
    END DO
    !
    ! -------------------------------------------------------------------- /
    ! 4.  Compute pars
    !
    npo    = -1
    !
    DO ip=0, npi
      !
      sumexp = 0.
      sumeyp = 0.
      !
      m0 = sume(ip)  * dth * tpiinv
      hs     = 4. * sqrt( max( m0 , 0. ) )
      IF ( hs .LT. hspmin ) THEN
        ! For wind cutoff and 2-band partitioning methods, keep the
        ! partition, but set the integrated parameters to UNDEF
        ! for Hs values less that HSPMIN:
        IF( ptmeth .EQ. 4 .OR. ptmeth .EQ. 5 ) THEN
          npo = npo + 1
          xp(:,npo) = undef
          xp(6,npo) = 0.0 ! Set wind sea frac to zero
        ENDIF
        cycle
      ENDIF
      !
      IF ( npo .GE. dimxp ) GOTO 2000
      npo = npo + 1
      IF (ip.GT.0)THEN
        IF(npo.LT.1)cycle
        pmap(npo) = ip
      ENDIF
      !
      m1 = sume1(ip) * dth * tpiinv**2
      m2 = sume2(ip) * dth * tpiinv**3
      mm1 = sumem1(ip) * dth
      qp = sumqp(ip) *(dth * tpiinv)**2
      !       M1 = MAX( M1, 1.E-7 )
      !       M2 = MAX( M2, 1.E-7 )
      !
      xl     = 1. / xfr - 1.
      xh     = xfr - 1.
      xl2    = xl**2
      xh2    = xh**2
      el     = sumf(ifpmax(ip)-1,ip) - sumf(ifpmax(ip),ip)
      eh     = sumf(ifpmax(ip)+1,ip) - sumf(ifpmax(ip),ip)
      denom  = xl*eh - xh*el
      sigp   = sig(ifpmax(ip))
      IF (denom.NE.0.) THEN
        sigp   = sigp *( 1. + 0.5 * ( xl2*eh - xh2*el ) &
             / sign( abs(denom) , denom ) )
      END IF
      CALL wavnu1 ( sigp, depth, wnp, cgp )
      !
      !/ --- Parabolic fit around the spectral peak ---
      ik = ifpmax(ip)
      efpmax(ip) = sumf(ik,ip) * dth
      IF (ik.GT.1 .AND. ik.LT.nk) THEN
        el    = sumf(ik-1,ip) * dth
        eh    = sumf(ik+1,ip) * dth
        numer = 0.125 * ( el - eh )**2
        denom = el - 2. * efpmax(ip) + eh
        IF (denom.NE.0.) efpmax(ip) = efpmax(ip)         &
             - numer / sign( abs(denom),denom )
      END IF
      !
      !/ --- Weighted least-squares regression to estimate frequency
      !/     spread (FSPRD) to an exponential function:
      !/              E(f) = A * exp(-1/2*(f-fp)/B)**2             ,
      !/     where B is frequency spread and  E(f) is used for
      !/     weighting to avoid greater weights to smalll values
      !/     in ordinary least-square fit. ---
      fsprd     = undef
      sumy      = 0.
      sumxy     = 0.
      sumxxy    = 0.
      sumylogy  = 0.
      sumxylogy = 0.
      !
      DO ik=1, nk
        y = sumf(ik,ip)*dth
        ! --- sums for weighted least-squares ---
        IF (y.GE.1.e-15) THEN
          yhat = log(y)
          xhat = -0.5 * ( (sig(ik)-sigp)*tpiinv )**2
          sumy      = sumy + y
          sumxy     = sumxy + xhat * yhat
          sumxxy    = sumxxy + xhat * xhat * y
          sumylogy  = sumylogy + y * yhat
          sumxylogy = sumxylogy + sumxy * yhat
        END IF
      END DO
      !
      numer = sumy * sumxxy - sumxy**2
      denom = sumy * sumxylogy - sumxy * sumylogy
      IF (denom.NE.0.)  fsprd = sqrt( numer / sign(abs(denom),numer) )
      !
      sumexp = sumfx(ifpmax(ip),ip) * dsii(ifpmax(ip))
      sumeyp = sumfy(ifpmax(ip),ip) * dsii(ifpmax(ip))
      !
      !/ --- Significant wave height ---
      xp(1,npo) = hs
      !/ --- Peak wave period ---
      xp(2,npo) = tpi / sigp
      !/ --- Peak wave length ---
      xp(3,npo) = tpi / wnp
      !/ --- Mean wave direction ---
      xp(4,npo) = mod( 630.-atan2(sumey(ip),sumex(ip))*rade , 360. )
      !/ --- Mean directional spread ---
      xp(5,npo) = rade * sqrt( max( 0. , 2. * ( 1. - sqrt( &
           max(0.,(sumex(ip)**2+sumey(ip)**2)/sume(ip)**2) ) ) ) )
      !/ --- Wind sea fraction ---
      xp(6,npo) = sumew(ip) / sume(ip)
      !/ --- Peak wave direction ---
      xp(7,npo) =  mod(630.-atan2(sumeyp,sumexp)*rade , 360.)
      !/ --- Spectral width (Longuet-Higgins 1975) ---
      xp(8,npo) = sqrt( max( 1. , m2*m0 / m1**2 ) - 1. )
      !/ --- JONSWAP peak enhancement parameter (E(fp)/EPM(fp))---
      !  EPM_FMX = ALPHA_PM_FMX * GRAV**2 * TPI * SIGP**-5 * EXP(-5/4)
      alp_pm = 0.3125 * hs**2 * (sigp)**4
      epm_fp = alp_pm * tpi * (sigp**(-5)) * 2.865048e-1
      xp(9,npo) = max( efpmax(ip) / epm_fp , 1.0 )
      !/ --- peakedness parameter (Goda 1970) ---
      xp(10,npo) = 2. * qp / m0**2
      !/ --- gaussian frequency width ---
      xp(11,npo) = fsprd
      !/ --- wave energy period (inverse moment) ---
      xp(12,npo) = mm1 / m0
      !/ --- mean wave period (first moment) ---
      xp(13,npo) = m0 / m1
      !/ --- zero-upcrossing period (second moment) ---
      xp(14,npo) = sqrt( m0 / m2 )
      !/ --- peak spectral density (one-dimensional) ---
      xp(15,npo) = efpmax(ip)
      !
    END DO
    !
    RETURN
    !
    ! Escape locations read errors --------------------------------------- *
    !
2000 CONTINUE
    IF ( iaproc .EQ. naperr ) WRITE (ndse,1000) npo+1
    RETURN
    !
    ! Formats
    !
1000 FORMAT (/' *** WAVEWATCH III ERROR IN PTMEAN :'/                &
         '     XP ARRAY TOO SMALL AT PARTITION',i6/)
    !/
    !/ End of PTMEAN ----------------------------------------------------- /
    !/

References constants::dera, w3gdatmd::dsii, w3gdatmd::dsip, w3gdatmd::dth, w3gdatmd::ecos, w3gdatmd::esin, w3gdatmd::fachfe, w3gdatmd::fte, w3odatmd::iaproc, w3odatmd::naperr, w3odatmd::ndse, w3odatmd::ndst, w3gdatmd::nk, w3gdatmd::nspec, w3gdatmd::nth, constants::rade, w3gdatmd::sig, w3servmd::strace(), w3gdatmd::th, constants::tpi, constants::tpiinv, constants::undef, w3dispmd::wavnu1(), and w3gdatmd::xfr.

Referenced by w3part().

ptnghb()

subroutine w3partmd::ptnghb

Nearest neighbour calculation.

Computes the nearest neighbors for each grid point. Wrapping of directional distribution (0 to 360) is taken care of using the nearest neighbor system

Parameters

[in]	IMI	Input discretized spectrum
[out]	IMD	Sorted data
[in]	IHMAX	Number of integer levels

Author

Barbara Tracy

Date

20 Oct 2006

Definition at line 632 of file w3partmd.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III          USACE/NOAA |
    !/                  |          Barbara  Tracy           |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         20-Oct-2006 !
    !/                  +-----------------------------------+
    !/
    !/    20-Oct-2006 : Origination.                       ( version 3.10 )
    !/
    !  1. Purpose :
    !
    !     This subroutine computes the nearest neighbors for each grid
    !     point. Wrapping of directional distribution (0 to 360)is taken
    !     care of using the nearest neighbor system
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       IMI     I.A.   I   Input discretized spectrum.
    !       IMD     I.A.   O   Sorted data.
    !       IHMAX   Int.   I   Number of integer levels.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Sur.  W3SERVMD Subroutine tracing.
    !     ----------------------------------------------------------------
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !
    USE w3gdatmd, ONLY: nk, nth, nspec
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    !     INTEGER, INTENT(IN)      :: IHMAX, IMI(NSPEC)
    !     INTEGER, INTENT(IN)      :: IMD(NSPEC)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: N, J, I, K
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
    CALL strace (ient, 'PTNGHB')
#endif
    !
    ! -------------------------------------------------------------------- /
    ! 1.  Check on need of processing
    !
    IF ( mk.EQ.nk .AND. mth.EQ.nth ) RETURN
    !
    IF ( mk.GT.0 ) DEALLOCATE ( neigh )
    ALLOCATE ( neigh(9,nspec) )
    mk     = nk
    mth    = nth
    !
    ! -------------------------------------------------------------------- /
    ! 2.  Build map
    !
    neigh  = 0
    !
    ! ... Base loop
    !
    DO n = 1, nspec
      !
      j      = (n-1) / nk + 1
      i      = n - (j-1) * nk
      k      = 0
      !
      ! ... Point at the left(1)
      !
      IF ( i .NE. 1 ) THEN
        k           = k + 1
        neigh(k, n) = n - 1
      END IF
      !
      ! ... Point at the right (2)
      !
      IF ( i .NE. nk ) THEN
        k           = k + 1
        neigh(k, n) = n + 1
      END IF
      !
      ! ... Point at the bottom(3)
      !
      IF ( j .NE. 1 ) THEN
        k           = k + 1
        neigh(k, n) = n - nk
      END IF
      !
      ! ... ADD Point at bottom_wrap to top
      !
      IF ( j .EQ. 1 ) THEN
        k          = k + 1
        neigh(k,n) = nspec - (nk-i)
      END IF
      !
      ! ... Point at the top(4)
      !
      IF ( j .NE. nth ) THEN
        k           = k + 1
        neigh(k, n) = n + nk
      END IF
      !
      ! ... ADD Point to top_wrap to bottom
      !
      IF ( j .EQ. nth ) THEN
        k          = k + 1
        neigh(k,n) = n - (nth-1) * nk
      END IF
      !
      ! ... Point at the bottom, left(5)
      !
      IF ( (i.NE.1) .AND. (j.NE.1) ) THEN
        k           = k + 1
        neigh(k, n) = n - nk - 1
      END IF
      !
      ! ... Point at the bottom, left with wrap.
      !
      IF ( (i.NE.1) .AND. (j.EQ.1) ) THEN
        k          = k + 1
        neigh(k,n) = n - 1 + nk * (nth-1)
      END IF
      !
      ! ... Point at the bottom, right(6)
      !
      IF ( (i.NE.nk) .AND. (j.NE.1) ) THEN
        k           = k + 1
        neigh(k, n) = n - nk + 1
      END IF
      !
      ! ... Point at the bottom, right with wrap
      !
      IF ( (i.NE.nk) .AND. (j.EQ.1) ) THEN
        k           = k + 1
        neigh(k,n) = n + 1 + nk * (nth - 1)
      END  IF
      !
      ! ... Point at the top, left(7)
      !
      IF ( (i.NE.1) .AND. (j.NE.nth) ) THEN
        k           = k + 1
        neigh(k, n) = n + nk - 1
      END IF
      !
      ! ... Point at the top, left with wrap
      !
      IF ( (i.NE.1) .AND. (j.EQ.nth) ) THEN
        k           = k + 1
        neigh(k,n) = n - 1 - (nk) * (nth-1)
      END IF
      !
      ! ... Point at the top, right(8)
      !
      IF ( (i.NE.nk) .AND. (j.NE.nth) ) THEN
        k           = k + 1
        neigh(k, n) = n + nk + 1
      END IF
      !
      ! ... Point at top, right with wrap
      !
      !
      IF ( (i.NE.nk) .AND. (j.EQ.nth) ) THEN
        k           = k + 1
        neigh(k,n) = n + 1 - (nk) * (nth-1)
      END IF
      !
      neigh(9,n) = k
      !
    END DO
    !
    RETURN
    !/
    !/ End of PTNGHB ----------------------------------------------------- /
    !/

References w3gdatmd::nk, w3gdatmd::nspec, w3gdatmd::nth, and w3servmd::strace().

Referenced by w3part().

ptsort()

subroutine w3partmd::ptsort	(	integer, dimension(nspec), intent(in)	IMI,
		integer, dimension(nspec), intent(out)	IND,
		integer, intent(in)	IHMAX
	)

Sorts the image data in ascending order.

This sort original to F. T. Tracy (2006)

Parameters

[in]	IMI	Input discretized spectrum
[out]	IND	Sorted data
[in]	IHMAX	Number of integer levels

Author

Barbara Tracy

Date

19 Oct 2006

Definition at line 513 of file w3partmd.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III          USACE/NOAA |
    !/                  |          Barbara  Tracy           |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         19-Oct-2006 !
    !/                  +-----------------------------------+
    !/
    !/    19-Oct-2006 : Origination.                       ( version 3.10 )
    !/
    !  1. Purpose :
    !
    !     This subroutine sorts the image data in ascending order.
    !     This sort original to F.T.Tracy (2006)
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       IMI     I.A.   I   Input discretized spectrum.
    !       IND     I.A.   O   Sorted data.
    !       IHMAX   Int.   I   Number of integer levels.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Sur.  W3SERVMD Subroutine tracing.
    !     ----------------------------------------------------------------
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !
    USE w3gdatmd, ONLY: nspec
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)      :: IHMAX, IMI(NSPEC)
    INTEGER, INTENT(OUT)     :: IND(NSPEC)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: I, IN, IV
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    INTEGER                 :: NUMV(IHMAX), IADDR(IHMAX),           &
         IORDER(NSPEC)
    !/
#ifdef W3_S
    CALL strace (ient, 'PTSORT')
#endif
    !
    ! -------------------------------------------------------------------- /
    ! 1.  Occurences per height
    !
    numv   = 0
    DO i=1, nspec
      numv(imi(i)) = numv(imi(i)) + 1
    END DO
    !
    ! -------------------------------------------------------------------- /
    ! 2.  Starting address per height
    !
    iaddr(1) = 1
    DO i=1, ihmax-1
      iaddr(i+1) = iaddr(i) + numv(i)
    END DO
    !
    ! -------------------------------------------------------------------- /
    ! 3.  Order points
    !
    DO i=1, nspec
      iv        = imi(i)
      in        = iaddr(iv)
      iorder(i) = in
      iaddr(iv) = in + 1
    END DO
    !
    ! -------------------------------------------------------------------- /
    ! 4.  Sort points
    !
    DO i=1, nspec
      ind(iorder(i)) = i
    END DO
    !
    RETURN
    !/
    !/ End of PTSORT ----------------------------------------------------- /
    !/

References w3gdatmd::nspec, and w3servmd::strace().

Referenced by w3part().

w3part()

subroutine w3partmd::w3part	(	real, dimension(nk,nth), intent(in)	SPEC,
		real, intent(in)	UABS,
		real, intent(in)	UDIR,
		real, intent(in)	DEPTH,
		real, dimension(nk), intent(in)	WN,
		integer, intent(out)	NP,
		real, dimension(dimp,0:dimxp), intent(out)	XP,
		integer, intent(in)	DIMXP
	)

Interface to watershed partitioning routines.

Watershed Algorithm of Vincent and Soille, 1991, implemented by Barbara Tracy (USACE/ERDC) for NOAA/NCEP.

This version of W3PART contains alternate Met Office partitioning methods, selected at runtime using the PTMETH namlist variable:

1. Standard WW3 partitioning, as per original method described by Barbary Tracy.
2. Met Office extended partitioning using split-partitions (removes the wind sea part of any swell partiton and combines with total wind sea partition).
3. Met Office "wave systems" - no classification or combining of wind sea partitions. All partitions output and ordered simply by wave height.
4. Classic, simple wave age based partitioning generating a single wind sea and swell partition.
5. 2-band partitioning; produces hi and low freqency band partitions using a user-defined cutoff frequency (PTFCUT).

Remarks

• DIMXP will always be of size 2 when using PTMETH 4 or 5.
• To achieve minimum storage but guaranteed storage of all partitions DIMXP = ((NK+1)/2) * ((NTH-1)/2) unless specified otherwise below.

Parameters

[in]	SPEC	2-D spectrum E(f,theta)
[in]	UABS	Wind speed
[in]	UDIR	Wind direction
[in]	DEPTH	Water depth
[in]	WN	Wavenumebers for each frequency
[out]	NP	Number of partitions found (-1=Spectrum without minumum energy; 0=Spectrum with minumum energy but no partitions)
[out]	XP	Parameters describing partitions. Entry '0' contains entire spectrum
[in]	DIMXP	Second dimension of XP

Author

Barbara Tracey, H. L. Tolman, M. Szyszka, Chris Bunney

Date

23 Jul 2018

Definition at line 139 of file w3partmd.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III          USACE/NOAA |
    !/                  |          Barbara  Tracy           |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         02-Dec-2010 !
    !/                  +-----------------------------------+
    !/
    !/    28-Oct-2006 : Origination.                       ( version 3.10 )
    !/    02-Dec-2010 : Adding a mapping PMAP between      ( version 3.14 )
    !/                  original and combined partitions
    !/                  ( M. Szyszka )
    !/
    !  1. Purpose :
    !
    !     Interface to watershed partitioning routines.
    !
    !  2. Method :
    !
    !     Watershed Algorithm of Vincent and Soille, 1991, implemented by
    !     Barbara Tracy (USACE/ERDC) for NOAA/NCEP.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       SPEC    R.A.   I   2-D spectrum E(f,theta).
    !       UABS    Real   I   Wind speed.
    !       UDIR    Real   I   Wind direction.
    !       DEPTH   Real   I   Water depth.
    !       WN      R.A.   I   Wavenumebers for each frequency.
    !       NP      Int.   O   Number of partitions.
    !                           -1 : Spectrum without minumum energy.
    !                            0 : Spectrum with minumum energy.
    !                                but no partitions.
    !       XP      R.A.   O   Parameters describing partitions.
    !                          Entry '0' contains entire spectrum.
    !       DIMXP   Int.   I   Second dimension of XP.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Sur.  W3SERVMD Subroutine tracing.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !     - To achieve minimum storage but guaranteed storage of all
    !       partitions DIMXP = ((NK+1)/2) * ((NTH-1)/2) unless specified
    !       otherwise below.
    !
    !     This version of W3PART contains alternate Met Office partitioning
    !     methods, selected at runtime using the PTMETH namlist variable:
    !         1) Standard WW3 partitioning
    !         2) Met Office extended partitioning using split-partitions
    !            (removes the wind sea part of any swell partiton and combines
    !            with total wind sea partition).
    !         3) Met Office "wave systems" - no classification or combining of
    !            wind sea partitions. All partitions output and ordered simply
    !            by wave height.
    !         4) Classic, simple wave age based partitioning generating
    !            a single wind sea and swell partition. [DIMXP = 2]
    !         5) 2-band partitioning; produces hi and low freqency band partitions
    !            using a user-defined cutoff frequency (PTFCUT). [DIMXP = 2]
    !
    !     (Chris Bunney, UK Met Office, Jul 2018)
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    !     !/S    Enable subroutine tracing.
    !     !/T    Enable test output
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    !/
    USE constants
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !
    USE w3gdatmd, ONLY: nk, nth, nspec, sig, th
    USE w3odatmd, ONLY: wscut, flcomb
    !
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(OUT)          :: NP
    INTEGER, INTENT(IN)           :: DIMXP
    REAL, INTENT(IN)              :: SPEC(NK,NTH), WN(NK), UABS,    &
         UDIR, DEPTH
    REAL, INTENT(OUT)             :: XP(DIMP,0:DIMXP)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: ITH, IMI(NSPEC), IMD(NSPEC),         &
         IMO(NSPEC), IND(NSPEC), NP_MAX,      &
         IP, IT(1), INDEX(DIMXP), NWS,        &
         IPW, IPT, ISP
    INTEGER                 :: PMAP(DIMXP)
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: ZP(NSPEC), ZMIN, ZMAX, Z(NSPEC),     &
         FACT, WSMAX, HSMAX
    REAL                    :: TP(DIMP,DIMXP)
    INTEGER                 :: IK, WIND_PART    ! ChrisB; added for new
    REAL                    :: C, UPAR, SIGCUT  ! UKMO partioning methods
    !/
    !/ ------------------------------------------------------------------- /
    ! 0.  Initializations
    !
#ifdef W3_S
    CALL strace (ient, 'W3PART')
#endif
    !
    np     = 0
    xp     = 0.
    !
    ! -------------------------------------------------------------------- /
    ! 1.  Process input spectrum
    ! 1.a 2-D to 1-D spectrum
    !
    DO ith=1, nth
      zp(1+(ith-1)*nk:ith*nk) = spec(:,ith)
    END DO
 
    !
    ! PTMETH == 4 : Do simple partitioning based solely on the
    ! wave age criterion (produces one swell and one wind sea only):
    !
    IF( ptmeth .EQ. 4 ) THEN
      DO ik=1, nk
        DO ith=1, nth
          isp = ik + (ith-1) * nk ! index into partition array IMO
 
          upar = wsmult * uabs * max(0.0, cos(th(ith)-dera*udir))
          c = sig(ik) / wn(ik)
 
          IF( upar .LE. c ) THEN
            ! Is swell:
            imo(isp) = 2
          ELSE
            ! Is wind sea:
            imo(isp) = 1
          ENDIF
        ENDDO
      ENDDO
 
      ! We have a max of up to two partitions:
      np_max=2
 
      ! Calculate mean parameters:
      CALL ptmean ( np_max, imo, zp, depth, uabs, udir, wn,           &
           np, xp, dimxp, pmap )
 
      ! No more processing required, return:
      RETURN
    ENDIF ! PTMETH == 4
    !
    ! PTMETH == 5 : produce "high" and "low" band partitions
    ! using a frequency cutoff:
    !
    IF( ptmeth .EQ. 5 ) THEN
      sigcut = tpi * ptfcut
      DO ik = 1, nk
        ! If bin center <= freq cutoff then mark as "low band".
        IF(sig(ik) .LE. sigcut) THEN
          ip = 2
        ELSE
          ip = 1
        ENDIF
 
        DO ith=1, nth
          isp = ik + (ith-1) * nk ! index into partition array IMO
          imo(isp) = ip
        ENDDO
      ENDDO
 
      ! We only ever have 2 partitions:
      np_max=2
 
      ! Calculate mean parameters:
      CALL ptmean ( np_max, imo, zp, depth, uabs, udir, wn,           &
           np, xp, dimxp, pmap )
 
      ! No more processing required, return:
      RETURN
    ENDIF ! PTMETH == 5
    !
    ! 1.b Invert spectrum and 'digitize'
    !
    zmin   = minval( zp )
    zmax   = maxval( zp )
    IF ( zmax-zmin .LT. 1.e-9 ) RETURN
    !
    z      = zmax - zp
    !
    fact   = real(ihmax-1) / ( zmax - zmin )
    imi    = max( 1 , min( ihmax , nint( 1. + z*fact ) ) )
    !
    ! 1.c Sort digitized image
    !
    CALL ptsort ( imi, ind, ihmax )
    !
    ! -------------------------------------------------------------------- /
    ! 2.  Perform partitioning
    ! 2.a Update nearest neighbor info as needed.
    !
    CALL ptnghb
    !
    ! 2.b Incremental flooding
    !
    CALL pt_fld ( imi, ind, imo, zp, np_max )
    !
    ! 2.c Compute parameters per partition
    !     NP and NX initialized inside routine.
    !
    CALL ptmean ( np_max, imo, zp, depth, uabs, udir, wn,           &
         np, xp, dimxp, pmap )
    !
    ! 2.d PTMETH == 2: move the wind sea part of the partitions into a
    !     seperate partition and recalculate the mean parameters.
    !
    IF ( np .GT. 0 .AND. ptmeth .EQ. 2 ) THEN
      wind_part = np_max
 
      DO ik=1, nk
        DO ith=1, nth
          isp = ik + (ith-1) * nk ! index into partition array IMO
          upar = wsmult * uabs * max(0.0, cos(th(ith)-dera*udir))
          c = sig(ik) / wn(ik)
 
          IF( c .LT. upar ) THEN
            ! Bin is wind forced - mark as new wind partition:
            wind_part = np_max + 1
 
            ! Update status map to show new wind partition
            imo(isp) = wind_part
          ENDIF
        ENDDO
      ENDDO
 
      IF( wind_part .NE. np_max ) THEN
        ! Some bins were marked as wind sea - recalculate
        ! integrated parameters:
        np_max = wind_part
        CALL ptmean ( np_max, imo, zp, depth, uabs, udir, wn,     &
             np, xp, dimxp, pmap )
      ENDIF
    ENDIF
    !
    ! -------------------------------------------------------------------- /
    ! 3.  Sort and recombine wind seas as needed
    ! 3.a Sort by wind sea fraction
    !
    IF ( np .LE. 1 ) RETURN
 
    ! -----------------------------------------------------------------
    ! PTMETH == 3: Don't classify or combine any partitions as wind sea.
    ! Simply sort by HS and return.
    ! -----------------------------------------------------------------
    IF( ptmeth .EQ. 3 ) THEN
      tp(:,1:np)  = xp(:,1:np)
      xp(:,1:np)  = 0.
 
      DO ip=1, np
        it          = maxloc(tp(1,1:np))
        xp(:,ip)    = tp(:,it(1))
        tp(1,it(1)) = -1.
      END DO
 
      RETURN ! Don't process any further
    ENDIF ! PTMETH == 3
 
    ! -----------------------------------------------------------------
    ! PTMETH == 1: Default WW3 partitioning.
    ! -----------------------------------------------------------------
    tp(:,1:np)  = xp(:,1:np)
    xp(:,1:np)  = 0.
    index(1:np) = 0
    nws         = 0
    !
    DO ip=1, np
      it          = maxloc(tp(6,1:np))
      index(ip)   = it(1)
      xp(:,ip)    = tp(:,index(ip))
      IF ( tp(6,it(1)) .GE. wscut ) nws = nws + 1
      tp(6,it(1)) = -1.
    END DO
    !
    ! 3.b Combine wind seas as needed and resort
    !
    IF ( nws.GT.1 .AND. flcomb ) THEN
      ipw    = pmap(index(1))
      DO ip=2, nws
        ipt    = pmap(index(ip))
        DO isp=1, nspec
          IF ( imo(isp) .EQ. ipt ) imo(isp) = ipw
        END DO
      END DO
      !
      CALL ptmean ( np_max, imo, zp, depth, uabs, udir, wn,       &
           np, xp, dimxp, pmap )
      IF ( np .LE. 1 ) RETURN
      !
      tp(:,1:np)  = xp(:,1:np)
      xp(:,1:np)  = 0.
      index(1:np) = 0
      nws         = 0
      !
      DO ip=1, np
        it          = maxloc(tp(6,1:np))
        index(ip)   = it(1)
        xp(:,ip)    = tp(:,index(ip))
        IF ( tp(6,it(1)) .GE. wscut ) nws = nws + 1
        tp(6,it(1)) = -1.
      END DO
      !
    END IF
    !
    ! 3.c Sort remaining fields by wave height
    !
    nws    = min( 1 , nws )
    !
    tp(:,1:np)  = xp(:,1:np)
    xp(:,1:np)  = 0.
    !
    IF ( nws .GT. 0 ) THEN
      xp(:,1) = tp(:,1)
      tp(1,1) = -1.
      nws     = 1
    END IF
    !
    DO ip=nws+1, np
      it          = maxloc(tp(1,1:np))
      xp(:,ip)    = tp(:,it(1))
      tp(1,it(1)) = -1.
    END DO
    !
    ! -------------------------------------------------------------------- /
    ! 4.  End of routine
    !
    RETURN
    !/
    !/ End of W3PART ----------------------------------------------------- /
    !/

References constants::dera, w3odatmd::flcomb, w3gdatmd::nk, w3gdatmd::nspec, w3gdatmd::nth, pt_fld(), ptmean(), ptnghb(), ptsort(), w3gdatmd::sig, w3servmd::strace(), w3gdatmd::th, constants::tpi, and w3odatmd::wscut.

Referenced by w3iosfmd::w3cprt(), w3exnc(), and w3expo().