w3partmd.F90

Go to the documentation of this file.

 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE w3partmd
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III          USACE/NOAA |
  !/                  |          Barbara  Tracy           |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         23-Jul-2018 |
  !/                  +-----------------------------------+
  !/
  !/    01-Nov-2006 : Origination.                        ( version 3.10 )
  !/    02-Nov-2006 : Adding tail to integration.         ( version 3.10 )
  !/    24-Mar-2007 : Bug fix IMI, adding overall field   ( version 3.11 )
  !/                  and sorting.
  !/    15-Apr-2008 : Clean up for distribution.          ( version 3.14 )
  !/    02-Dec-2010 : Adding a mapping PMAP between       ( version 3.14 )
  !/                  original and combined partitions
  !/                  ( M. Szyszka )
  !/    23-Jul-2018 : Added alternative partitioning      ( version 6.05 )
  !/                  methods (C. Bunney, UKMO)
  !  1. Purpose :
  !
  !     Spectral partitioning according to the watershed method.
  !
  !  2. Variables and types :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      MK, MTH   Int.  Private  Dimensions of stored neighour array.
  !      NEIGH     I.A.  Private  Nearest Neighbor array.
  !     ----------------------------------------------------------------
  !      Note: IHMAX, HSPMIN, WSMULT, WSCUT and FLCOMB used from W3ODATMD.
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3PART    Subr. Public   Interface to watershed routines.
  !      PTSORT    Subr. Public   Sort discretized image.
  !      PTNGHB    Subr. Public   Defeine nearest neighbours.
  !      PT_FLD    Subr. Public   Incremental flooding algorithm.
  !      FIFO_ADD, FIFO_EMPTY, FIFO_FIRST
  !                Subr. PT_FLD   Queue management.
  !      PTMEAN    Subr. Public   Compute mean parameters.
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      STRACE    Subr. W3SERVMD Subroutine traceing.
  !      WAVNU1    Subr. W3DISPMD Wavenumber computation.
  !     ----------------------------------------------------------------
  !
  !  5. Remarks :
  !
  !  6. Switches :
  !
  !     !/S    Enable subroutine tracing.
  !     !/T    Enable test output
  !
  !  7. Source code :
  !
  !/ ------------------------------------------------------------------- /
  !
  USE w3odatmd, ONLY: ihmax, hspmin, wsmult, dimp, ptmeth, ptfcut
  !
  PUBLIC
  !
  INTEGER, PRIVATE              :: MK = -1
  INTEGER, PRIVATE              :: MTH = -1
  INTEGER, ALLOCATABLE, PRIVATE :: NEIGH(:,:)
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3part ( SPEC, UABS, UDIR, DEPTH, WN, NP, XP, DIMXP )
    !/
    !/                  +-----------------------------------+
    !/                  | 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 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3part
  !/ ------------------------------------------------------------------- /
  SUBROUTINE ptsort ( IMI, IND, IHMAX )
    !/
    !/                  +-----------------------------------+
    !/                  | 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 ----------------------------------------------------- /
    !/
  END SUBROUTINE ptsort
  !/ ------------------------------------------------------------------- /
  SUBROUTINE ptnghb
    !/
    !/                  +-----------------------------------+
    !/                  | 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 ----------------------------------------------------- /
    !/
  END SUBROUTINE ptnghb
  !/ ------------------------------------------------------------------- /
  SUBROUTINE pt_fld ( IMI, IND, IMO, ZP, NPART )
    !/
    !/                  +-----------------------------------+
    !/                  | 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 ----------------------------------------------------- /
    !/
  END SUBROUTINE pt_fld
  !/ ------------------------------------------------------------------- /
  SUBROUTINE ptmean ( NPI, IMO, ZP, DEPTH, UABS, UDIR, WN,        &
       NPO, XP, DIMXP, PMAP )
    !/
    !/                  +-----------------------------------+
    !/                  | 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 ----------------------------------------------------- /
    !/
  END SUBROUTINE ptmean
  !/
  !/ End of module W3PARTMD -------------------------------------------- /
  !/
END MODULE w3partmd