!> @file w3partmd.F90
!> @brief Spectral partitioning module
!>
!> @authors Barbara Tracy, H. L. Tolman, M. Szyszka, Chris Bunney
!> @date 23 Jul 2018
!>
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
!> @brief Spectral partitioning according to the watershed method
!>
!> @details Multiple partitioning methods are provided in this
!> module that can be selected via the \c PTMETH namelist variable.
!> Please see the w3part() subroutine for details
!>
!> @author Barbara Tracey, H. L. Tolman, M. Szyszka, Chris Bunney
!> @date 23 Jul 2018
!>
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
!
!> Nearest neighbour array frequency dimension size
INTEGER, PRIVATE :: MK = -1
!> Nearest neighbour array direction dimension size
INTEGER, PRIVATE :: MTH = -1
!> Nearest neighbour array
INTEGER, ALLOCATABLE, PRIVATE :: NEIGH(:,:)
!/
CONTAINS
!/ ------------------------------------------------------------------- /
!> @brief Interface to watershed partitioning routines.
!>
!> @details 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 \c PTMETH namlist variable:
!> -# Standard WW3 partitioning, as per original method described
!> by Barbary Tracy.
!> -# Met Office extended partitioning using split-partitions
!> (removes the wind sea part of any swell partiton and combines
!> with total wind sea partition).
!> -# Met Office "wave systems" - no classification or combining of
!> wind sea partitions. All partitions output and ordered simply
!> by wave height.
!> -# Classic, simple wave age based partitioning generating
!> a single wind sea and swell partition.
!> -# 2-band partitioning; produces hi and low freqency band partitions
!> using a user-defined cutoff frequency (\c PTFCUT).
!>
!> @remarks
!> - \c DIMXP will always be of size 2 when using \c PTMETH 4 or 5.
!>
!> - To achieve minimum storage but guaranteed storage of all
!> partitions <tt>DIMXP = ((NK+1)/2) * ((NTH-1)/2)</tt>
!> unless specified otherwise below.
!>
!> @param[in] SPEC 2-D spectrum E(f,theta)
!> @param[in] UABS Wind speed
!> @param[in] UDIR Wind direction
!> @param[in] DEPTH Water depth
!> @param[in] WN Wavenumebers for each frequency
!> @param[out] NP Number of partitions found
!> (-1=Spectrum without minumum energy;
!> 0=Spectrum with minumum energy but no partitions)
!> @param[out] XP Parameters describing partitions.
!> Entry '0' contains entire spectrum
!> @param[in] DIMXP Second dimension of XP
!>
!> @author Barbara Tracey, H. L. Tolman, M. Szyszka, Chris Bunney
!> @date 23 Jul 2018
!>
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
!/ ------------------------------------------------------------------- /
!>
!> @brief Sorts the image data in ascending order.
!>
!> @details This sort original to F. T. Tracy (2006)
!>
!> @param[in] IMI Input discretized spectrum
!> @param[out] IND Sorted data
!> @param[in] IHMAX Number of integer levels
!>
!> @author Barbara Tracy
!> @date 19 Oct 2006
!>
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
!/ ------------------------------------------------------------------- /
!>
!> @brief Nearest neighbour calculation
!>
!> @details
!> Computes the nearest neighbors for each grid point. Wrapping of
!> directional distribution (0 to 360) is taken care of using the
!> nearest neighbor system
!>
!> @param[in] IMI Input discretized spectrum
!> @param[out] IMD Sorted data
!> @param[in] IHMAX Number of integer levels
!>
!> @author Barbara Tracy
!> @date 20 Oct 2006
!>
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
#endif
!/
#ifdef W3_S
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
!/ ------------------------------------------------------------------- /
!>
!> @brief Image watersheding
!>
!> @details
!> This subroutine does incremental flooding of the image to
!> determine the watershed image.
!>
!> @param[in] IMI Input discretized spectrum
!> @param[in] IND Sorted addresses
!> @param[out] IMO Output partitioned spectrum
!> @param[in] ZP Spectral array
!> @param[out] NPART Number of partitions found
!>
!> @author H.L. Tolman
!> @date 01 Nov 2006
!>
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
!/ ------------------------------------------------------------------- /
!> @brief Add point to FIFO queue
!>
!> @param[in] IV Point to add
!>
!> @author Barbara Tracy
!> @date 01 Nov 2006
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
!/ ------------------------------------------------------------------- /
!> @brief Check if queue is empty.
!>
!> @param[out] IEMPTY Set to 1 if queue is empty, else 0
!>
!> @author Barbara Tracy
!> @date 01 Nov 2006
!>
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
!/ ------------------------------------------------------------------- /
!> @brief Get point out of queue.
!>
!> @param[out] IV Returned point
!>
!> @author Barbara Tracy
!> @date 01 Nov 2006
!>
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
!/
!/ End of PT_FLD ----------------------------------------------------- /
!/
END SUBROUTINE PT_FLD
!/ ------------------------------------------------------------------- /
!> @brief Compute mean parameters per partition
!>
!> @param[in] NPI Number of partitions found.
!> @param[in] IMO Partition map.
!> @param[in] ZP Input spectrum.
!> @param[in] DEPTH Water depth.
!> @param[in] UABS Wind speed.
!> @param[in] UDIR Wind direction.
!> @param[in] WN Wavenumebers for each frequency.
!> @param[out] NPO Number of partitions with mean parameters.
!> @param[out] XP Array with output parameters.
!> @param[in] DIMXP Second dimension of XP.
!> @param[out] PMAP Mapping between orig. and combined partitions
!>
!> @author Barbara Tracy, H. L. Tolman, M. Szyszka, C. Bunney
!> @date 02 Dec 2010
!>
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