#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, 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 !/S USE W3SERVMD, ONLY: STRACE ! 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) !/S INTEGER, SAVE :: IENT = 0 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 ! !/S CALL STRACE (IENT, 'W3PART') ! 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 : ! !/ ------------------------------------------------------------------- / ! !/S USE W3SERVMD, ONLY: STRACE ! 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 !/S INTEGER, SAVE :: IENT = 0 INTEGER :: NUMV(IHMAX), IADDR(IHMAX), & IORDER(NSPEC) !/ !/S CALL STRACE (IENT, 'PTSORT') ! ! -------------------------------------------------------------------- / ! 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 : ! !/ ------------------------------------------------------------------- / ! !/S USE W3SERVMD, ONLY: STRACE ! 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 !/S INTEGER, SAVE :: IENT = 0 !/ !/S CALL STRACE (IENT, 'PTNGHB') ! ! -------------------------------------------------------------------- / ! 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 : ! !/ ------------------------------------------------------------------- / ! !/S USE W3SERVMD, ONLY: STRACE ! 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 !/S INTEGER, SAVE :: IENT = 0 REAL :: ZPMAX, EP1, DIFF !/ !/S CALL STRACE (IENT, 'PT_FLD') ! ! -------------------------------------------------------------------- / ! 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 ) ! ! Add point to FIFO queue. ! INTEGER, INTENT(IN) :: IV ! IQ(IQ_END) = IV ! IQ_END = IQ_END + 1 IF ( IQ_END .GT. NSPEC ) IQ_END = 1 ! RETURN END SUBROUTINE !/ ------------------------------------------------------------------- / SUBROUTINE FIFO_EMPTY ( IEMPTY ) ! ! Check if queue is empty. ! INTEGER, INTENT(OUT) :: IEMPTY ! IF ( IQ_START .NE. IQ_END ) THEN IEMPTY = 0 ELSE IEMPTY = 1 END IF ! RETURN END SUBROUTINE !/ ------------------------------------------------------------------- / SUBROUTINE FIFO_FIRST ( IV ) ! ! Get point out of queue. ! 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 !/ ------------------------------------------------------------------- / 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 !/S USE W3SERVMD, ONLY: STRACE 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) !/S INTEGER, SAVE :: IENT = 0 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 !/ !/S CALL STRACE (IENT, 'PTMEAN') ! ! -------------------------------------------------------------------- / ! 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