! Contents of this file ! ! W3ODATMD parameters required for partitioning model output ! SWPARTMD spectral partitioning according to the watershed method ! MODULE W3ODATMD !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | H. L. Tolman | !/ | FORTRAN 90 | !/ | Last update : 27-Jul-2010 | !/ +-----------------------------------+ !/ !/ Copyright 2009-2010 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ ! 1. Purpose : ! ! This module considers the parameters required for model output. ! (Module truncated from original WW3 version W3ODATMD.FTN) ! ! 2. Variables and types : ! ! Elements of OUT6 are aliased to pointers with the same ! name. These pointers are defined as : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! IHMAX Int. Public Number of discrete spectral levels. ! HSPMIN Real Public Minimum significant height per part. ! WSMULT Real Public Multiplier for wind sea boundary. ! WSCUT Real Public Cut-off wind factor for wind seas. ! FLCOMB Log. Public Flag for combining wind seas. ! ---------------------------------------------------------------- ! ! 3. Subroutines and functions : ! ! 4. Subroutines and functions used : ! ! 5. Remarks : ! ! 6. Switches : ! ! 7. Source code : ! !/ ------------------------------------------------------------------- / PUBLIC INTEGER :: IHMAX REAL :: HSPMIN, WSMULT, WSCUT LOGICAL :: FLCOMB ! Values taken from ww3_grid.ftn PARAMETER (IHMAX=100, HSPMIN=0.05, WSMULT=1.7, 40.PART !Can also be moved to user input -> OUTPUT command & WSCUT=0.333, FLCOMB=.FALSE.) !/ CONTAINS !/ !/ End of module W3ODATMD -------------------------------------------- / !/ END MODULE W3ODATMD MODULE SWPARTMD !/ !/ +-----------------------------------+ !/ | WAVEWATCH III USACE/NOAA | !/ | Barbara Tracy | !/ | H. L. Tolman | !/ | FORTRAN 90 | !/ | Last update : 15-Apr-2008 | !/ +-----------------------------------+ !/ !/ 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 ) !/ ! 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 ! ---------------------------------------------------------------- ! SWPART 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 40.PART ! PUBLIC ! INTEGER, PRIVATE :: MK = -1, MTH = -1 INTEGER, ALLOCATABLE, PRIVATE :: NEIGH(:,:) !/ CONTAINS !/ ------------------------------------------------------------------- / SUBROUTINE SWPART ( SPEC, UABS, UDIR, DEPTH, WN, SPCSIG, SPCDIR, & NP, XP, DIMXP ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III USACE/NOAA | !/ | Barbara Tracy | !/ | H. L. Tolman | !/ | FORTRAN 90 | !/ | Last update : 28-Oct-2006 ! !/ +-----------------------------------+ !/ !/ 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 (Cartesian). 40.PART ! DEPTH Real I Water depth. ! WN R.A. I Wavenumbers for each frequency. ! SPCDIR REAL I Spectral directions 40.PART ! SPCSIG REAL I Spectral frequencies 40.PART ! 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 = ((MSC+1)/2) * ((MDC-1)/2) ! ! 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 USE W3ODATMD, ONLY: WSCUT, FLCOMB 40.PART USE SWCOMM3, ONLY: MSC, MDC 40.PART USE OCPCOMM4, ONLY: PRINTF 40.PART !REMOVE!! ! IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(OUT) :: NP INTEGER, INTENT(IN) :: DIMXP REAL, INTENT(IN) :: SPEC(MSC,MDC), WN(MSC), UABS, & UDIR, DEPTH REAL, INTENT(OUT) :: XP(7,0:DIMXP) REAL, INTENT(IN) :: SPCDIR(MDC,6) 40.PART REAL, INTENT(IN) :: SPCSIG(MSC) 40.PART !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: ITH, IMI(MSC*MDC), IMD(MSC*MDC), & IMO(MSC*MDC), IND(MSC*MDC), NP_MAX, & IP, IT(1), INDEX(DIMXP), NWS, & IPW, IPT, ISP INTEGER :: PMAP(DIMXP) !/S INTEGER, SAVE :: IENT = 0 REAL :: ZP(MSC*MDC), ZMIN, ZMAX, Z(MSC*MDC), & FACT, WSMAX, HSMAX REAL :: TP(7,DIMXP) 40.PART Extended from TP(6,1:DIMXP) !/ !/ ------------------------------------------------------------------- / ! 0. Initializations ! !/S CALL STRACE (IENT, 'SWPART') ! NP = 0 XP = 0. ! ! -------------------------------------------------------------------- / ! 1. Process input spectrum ! 1.a 2-D to 1-D spectrum ! DO ITH=1, MDC ZP(1+(ITH-1)*MSC:ITH*MSC) = SPEC(:,ITH) END DO ! ! 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, & SPCSIG, SPCDIR, NP, XP, DIMXP, PMAP ) 40.PART ! ! -------------------------------------------------------------------- / ! 3. Sort and recombine wind seas as needed ! 3.a Sort by wind sea fraction ! 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 ! ! 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, MSC*MDC IF ( IMO(ISP) .EQ. IPT ) IMO(ISP) = IPW END DO END DO ! CALL PTMEAN ( NP_MAX, IMO, ZP, DEPTH, UABS, UDIR, WN, & SPCSIG, SPCDIR, NP, XP, DIMXP, PMAP ) 40.PART 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 !> IF (ALLOCATED(NEIGH )) DEALLOCATE(NEIGH ) 40.PART ! ! -------------------------------------------------------------------- / ! 4. End of routine ! RETURN !/ !/ End of SWPART ----------------------------------------------------- / !/ END SUBROUTINE SWPART !/ ------------------------------------------------------------------- / 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 USE SWCOMM3, ONLY: MSC, MDC 40.PART ! IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: IHMAX, IMI(MSC*MDC) INTEGER, INTENT(OUT) :: IND(MSC*MDC) !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: I, IN, IV !/S INTEGER, SAVE :: IENT = 0 INTEGER :: NUMV(IHMAX), IADDR(IHMAX), & IORDER(MSC*MDC) !/ !/S CALL STRACE (IENT, 'PTSORT') ! ! -------------------------------------------------------------------- / ! 1. Occurences per height ! NUMV = 0 DO I=1, MSC*MDC 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, MSC*MDC IV = IMI(I) IN = IADDR(IV) IORDER(I) = IN IADDR(IV) = IN + 1 END DO ! ! -------------------------------------------------------------------- / ! 4. Sort points ! DO I=1, MSC*MDC 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 USE SWCOMM3, ONLY: MSC, MDC 40.PART USE OCPCOMM4, ONLY: PRINTF 40.PART !REMOVE!! ! IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ ! INTEGER, INTENT(IN) :: IHMAX, IMI(MSC*MDC) ! INTEGER, INTENT(IN) :: IMD(MSC*MDC) !/ !/ ------------------------------------------------------------------- / !/ 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.MSC .AND. MTH.EQ.MDC ) RETURN ! IF ( MK.GT.0 ) DEALLOCATE ( NEIGH ) ALLOCATE ( NEIGH(9,MSC*MDC) ) MK = MSC MTH = MDC ! ! -------------------------------------------------------------------- / ! 2. Build map ! NEIGH = 0 ! ! ... Base loop ! DO N = 1, MSC*MDC ! J = (N-1) / MSC + 1 I = N - (J-1) * MSC 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. MSC ) 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 - MSC END IF ! ! ... ADD Point at bottom_wrap to top ! IF ( J .EQ. 1 ) THEN K = K + 1 NEIGH(K,N) = MSC*MDC - (MSC-I) END IF ! ! ... Point at the top(4) ! IF ( J .NE. MDC ) THEN K = K + 1 NEIGH(K, N) = N + MSC END IF ! ! ... ADD Point to top_wrap to bottom ! IF ( J .EQ. MDC ) THEN K = K + 1 NEIGH(K,N) = N - (MDC-1) * MSC END IF ! ! ... Point at the bottom, left(5) ! IF ( (I.NE.1) .AND. (J.NE.1) ) THEN K = K + 1 NEIGH(K, N) = N - MSC - 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 + MSC * (MDC-1) END IF ! ! ... Point at the bottom, right(6) ! IF ( (I.NE.MSC) .AND. (J.NE.1) ) THEN K = K + 1 NEIGH(K, N) = N - MSC + 1 END IF ! ! ... Point at the bottom, right with wrap ! IF ( (I.NE.MSC) .AND. (J.EQ.1) ) THEN K = K + 1 NEIGH(K,N) = N + 1 + MSC * (MDC - 1) END IF ! ! ... Point at the top, left(7) ! IF ( (I.NE.1) .AND. (J.NE.MDC) ) THEN K = K + 1 NEIGH(K, N) = N + MSC - 1 END IF ! ! ... Point at the top, left with wrap ! IF ( (I.NE.1) .AND. (J.EQ.MDC) ) THEN K = K + 1 NEIGH(K,N) = N - 1 - (MSC) * (MDC-1) END IF ! ! ... Point at the top, right(8) ! IF ( (I.NE.MSC) .AND. (J.NE.MDC) ) THEN K = K + 1 NEIGH(K, N) = N + MSC + 1 END IF ! ! ... Point at top, right with wrap ! ! IF ( (I.NE.MSC) .AND. (J.EQ.MDC) ) THEN K = K + 1 NEIGH(K,N) = N + 1 - (MSC) * (MDC-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 USE SWCOMM3, ONLY: MSC, MDC 40.PART ! IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: IMI(MSC*MDC), IND(MSC*MDC) INTEGER, INTENT(OUT) :: IMO(MSC*MDC), NPART REAL, INTENT(IN) :: ZP(MSC*MDC) !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: MASK, INIT, IWSHED, IMD(MSC*MDC), & IC_LABEL, IFICT_PIXEL, M, IH, MSAVE, & IP, I, IPP, IC_DIST, IEMPTY, IPPP, & JL, JN, IPT, J INTEGER :: IQ(MSC*MDC), 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. MSC*MDC ) 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. MSC*MDC ) 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 , MSC*MDC 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. MSC*MDC ) 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. MSC*MDC ) IQ_START = 1 ! RETURN END SUBROUTINE !/ !/ End of PT_FLD ----------------------------------------------------- / !/ END SUBROUTINE PT_FLD !/ ------------------------------------------------------------------- / SUBROUTINE PTMEAN ( NPI, IMO, ZP, DEPTH, UABS, UDIR, WN, & SPCSIG, SPCDIR, NPO, XP, DIMXP, PMAP ) 40.PART !/ !/ +-----------------------------------+ !/ | WAVEWATCH III USACE/NOAA | !/ | Barbara Tracy | !/ | H. L. Tolman | !/ | FORTRAN 90 | !/ | Last update : 24-Mar-2007 ! !/ +-----------------------------------+ !/ !/ 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 pean 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 (Cartesian). 40.PART ! WN R.A. I Wavenumebers for each frequency. ! SPCDIR REAL I Spectral directions 40.PART ! SPCSIG REAL I Spectral frequencies 40.PART ! NPO Int. O Number of partitions with mean parameters. ! XP R.A. O Array with output parameters. ! DIMXP int. I Second dimesion 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 SWCOMM3, ONLY: MSC, MDC, DDIR, PI2, DEGRAD, FRINTF 40.PART USE SWCOMM1, ONLY: OUTPAR 40.PARTWL USE OCPCOMM4, ONLY: PRINTF 40.PART !REMOVE!! ! IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: NPI, IMO(MSC*MDC), DIMXP INTEGER, INTENT(OUT) :: NPO, PMAP(DIMXP) REAL, INTENT(IN) :: ZP(MSC*MDC), DEPTH, UABS, UDIR, WN(MSC) REAL, INTENT(OUT) :: XP(7,0:DIMXP) REAL, INTENT(IN) :: SPCDIR(MDC,6) 40.PART REAL, INTENT(IN) :: SPCSIG(MSC) 40.PART !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: IK, ITH, ISP, IP, IFPMAX(0:NPI) !/S INTEGER, SAVE :: IENT = 0 REAL :: SUMF(0:MSC+1,0:NPI), SUMFW(MSC,0:NPI), & SUMFX(MSC,0:NPI), SUMFY(MSC,0:NPI), & SUME(0:NPI), SUMEW(0:NPI), & SUMEX(0:NPI), SUMEY(0:NPI), & EFPMAX(0:NPI), FCDIR(MDC) REAL :: SUMFK(MSC,0:NPI), SUMEK(0:NPI) 40.PARTWL REAL :: HS, XL, XH, XL2, XH2, EL, EH, DENOM, & SIGP, WNP, CGP, UPAR, C(MSC), RD, FACT REAL :: DS, FTE 40.PARTWL !/ !/S CALL STRACE (IENT, 'PTMEAN') ! ! -------------------------------------------------------------------- / ! WRITE(PRINTF,*) 'Inside PTMEAN: UABS, UDIR =', UABS, UDIR 40.PART ! 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. SUMFK = 0. 40.PARTWL SUME = 0. SUMEW = 0. SUMEX = 0. SUMEY = 0. SUMEK = 0. 40.PARTWL IFPMAX = 0 EFPMAX = 0. ! DO IK=1, MSC C(IK) = SPCSIG(IK) / WN(IK) ! WRITE(PRINTF,*) 'In PTMEAN:,IK,SPCSIG(IK),WN(IK),C(IK) =', 40.PART ! & IK,SPCSIG(IK),WN(IK),C(IK) 40.PART END DO ! DO ITH=1, MDC UPAR = WSMULT * UABS * MAX(0.,COS(SPCDIR(ITH,1)-DEGRAD*UDIR)) ! WRITE(PRINTF,*) 'In PTMEAN:,ITH,SPCDIR(ITH,1),UDIR,UPAR =', 40.PART ! & ITH,SPCDIR(ITH,1)/DEGRAD,UDIR,UPAR 40.PART IF ( UPAR .LT. C(MSC) ) THEN !> FCDIR(ITH) = SPCSIG(MSC+1) 40.PART !Could be a problem to have MSC+1 frequencies!! FCDIR(ITH) = SPCSIG(MSC)*(1+FRINTF) 40.PART !Took away last frequency - CHECK REASON FOR THIS!!! ELSE DO IK=MSC-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*SPCSIG(IK+1) + (1.-RD)*SPCSIG(IK) END IF ! WRITE(PRINTF,*) 'In PTMEAN: ITH, FCDIR(ITH) =',ITH, FCDIR(ITH) 40.PART END DO ! ! -------------------------------------------------------------------- / ! 3. Spectral integrals and preps ! 3.a Integrals ! NOTE: Factor DDIR only used in Hs computation. ! DO IK=2, MSC 40.PART !For SWAN, counter is started at IK=2! DO ITH=1, MDC ISP = IK + (ITH-1)*MSC IP = IMO(ISP) DS = SPCSIG(IK)-SPCSIG(IK-1) FACT = MAX ( 0. , MIN ( 1. , & 1.-( FCDIR(ITH) - 0.5*(SPCSIG(IK-1)+SPCSIG(IK)) )/DS) ) 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) * SPCDIR(ITH,2) SUMFY(IK, 0) = SUMFY(IK, 0) + ZP(ISP) * SPCDIR(ITH,3) SUMFK(IK, 0) = SUMFK(IK, 0) + ZP(ISP) * WN(IK)**OUTPAR(3) 40.PARTWL !WLEN - ADD tail!!! 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) * SPCDIR(ITH,2) SUMFY(IK,IP) = SUMFY(IK,IP) + ZP(ISP) * SPCDIR(ITH,3) SUMFK(IK,IP) = SUMFK(IK,IP) + ZP(ISP) * WN(IK)**OUTPAR(3) 40.PARTWL !WLEN - ADD tail!!! END DO END DO !> SUMF(MSC+1,:) = SUMF(MSC,:) * FACHFE 40.PART !TAIL ADDITION - TEMPORARILY DEACTIVATED!!! ! DO IP=0, NPI DO IK=2, MSC 40.PART !For SWAN, counter is started at IK=2! DS = SPCSIG(IK)-SPCSIG(IK-1) ! SUME (IP) = SUME (IP) + SUMF (IK,IP) * DS ! SUMEW(IP) = SUMEW(IP) + SUMFW(IK,IP) * DS ! SUMEX(IP) = SUMEX(IP) + SUMFX(IK,IP) * DS ! SUMEY(IP) = SUMEY(IP) + SUMFY(IK,IP) * DS !Replaced original with more accurate trapezoidal rule 40.PART SUME (IP) = SUME (IP) + 0.5*(SPCSIG(IK)*SUMF(IK,IP) + 40.PART & SPCSIG(IK-1)*SUMF(IK-1,IP))*DS 40.PART SUMEW(IP) = SUMEW(IP) + 0.5*(SPCSIG(IK)*SUMFW(IK,IP) + 40.PART & SPCSIG(IK-1)*SUMFW(IK-1,IP))*DS 40.PART SUMEX(IP) = SUMEX(IP) + 0.5*(SPCSIG(IK)*SUMFX(IK,IP) + 40.PART & SPCSIG(IK-1)*SUMFX(IK-1,IP))*DS 40.PART SUMEY(IP) = SUMEY(IP) + 0.5*(SPCSIG(IK)*SUMFY(IK,IP) + 40.PART & SPCSIG(IK-1)*SUMFY(IK-1,IP))*DS 40.PART SUMEK(IP) = SUMEK(IP) + 0.5*(SPCSIG(IK)*SUMFK(IK,IP) + 40.PARTWL !WLEN - ADD tail!!! & SPCSIG(IK-1)*SUMFK(IK-1,IP))*DS 40.PARTWL !WLEN - ADD tail!!! IF ( SUMF(IK,IP) .GT. EFPMAX(IP) ) THEN IFPMAX(IP) = IK EFPMAX(IP) = SUMF(IK,IP) END IF END DO !> FTE = 0.25 * SPCSIG(MSC) * DDIR * SPCSIG(MSC) 40.PART !Added locally. Not so sure about this - CHECK!! Check whether DDIR is defined in Rads! FTE=0. 40.PART !REMOVE!! SUME (IP) = SUME (IP) + SUMF (MSC,IP) * FTE SUMEW(IP) = SUMEW(IP) + SUMFW(MSC,IP) * FTE SUMEX(IP) = SUMEX(IP) + SUMFX(MSC,IP) * FTE SUMEY(IP) = SUMEY(IP) + SUMFY(MSC,IP) * FTE ! WRITE(PRINTF,*) 'In PTMEAN: IP, SUME (IP) =',IP, SUME (IP) 40.PART ! WRITE(PRINTF,*) 'In PTMEAN: IP, SUMEW(IP) =',IP, SUMEW(IP) 40.PART ! WRITE(PRINTF,*) 'In PTMEAN: IP, SUMEX(IP) =',IP, SUMEX(IP) 40.PART ! WRITE(PRINTF,*) 'In PTMEAN: IP, SUMEY(IP) =',IP, SUMEY(IP) 40.PART END DO !> FROM SWANOUT1.FTN: !>! integration over [0,inf] 40.87 !> ETOT = 0. !>! trapezoidal rule is applied !> DO ID=1, MDC !> DO IS=2,MSC !> DS=SPCSIG(IS)-SPCSIG(IS-1) 30.72 !> EAD = 0.5*(SPCSIG(IS)*ACLOC(ID,IS)+ 30.72 !> & SPCSIG(IS-1)*ACLOC(ID,IS-1))*DS*DDIR 30.72 !> ETOT = ETOT + EAD !> ENDDO !> IF (MSC .GT. 3) THEN 10.20 !>! contribution of tail to total energy density !> EHFR = ACLOC(ID,MSC) * SPCSIG(MSC) 30.72 !> ETOT = ETOT + DDIR * EHFR * SPCSIG(MSC) * EFTAIL 30.72 !> ENDIF !> ENDDO !> !> Wave length ! ETOT = 0. ! EKTOT = 0. !! new integration method involving FRINTF 20.59 ! DO IS=1, MSC ! SIG2 = (SPCSIG(IS))**2 30.72 ! SKK = SIG2 * (WK(IS))**OUTPAR(3) 40.00 ! DO ID=1,MDC ! ETOT = ETOT + SIG2 * ACLOC(ID,IS) 20.59 ! EKTOT = EKTOT + SKK * ACLOC(ID,IS) 20.59 ! ENDDO ! ENDDO ! ETOT = FRINTF * ETOT ! EKTOT = FRINTF * EKTOT ! IF (MSC .GT. 3) THEN 10.20 !! contribution of tail to total energy density ! PPTAIL = PWTAIL(1) - 1. 20.59 ! CETAIL = 1. / (PPTAIL * (1. + PPTAIL * (FRINTH-1.))) 20.59 ! PPTAIL = PWTAIL(1) - 1. - 2.*OUTPAR(3) 40.00 ! IF (PPTAIL.LE.0.) THEN ! CALL MSGERR (2,'error tail computation') ! GOTO 480 ! ENDIF ! CKTAIL = 1. / (PPTAIL * (1. + PPTAIL * (FRINTH-1.))) 20.59 ! DO ID=1,MDC ! ETOT = ETOT + CETAIL * SIG2 * ACLOC(ID,MSC) 20.59 ! EKTOT = EKTOT + CKTAIL * SKK * ACLOC(ID,MSC) 20.59 ! ENDDO ! 480 CONTINUE ! ENDIF ! IF (EKTOT.GT.0.) THEN ! WLMEAN = PI2 * (ETOT / EKTOT) ** (1./OUTPAR(3)) 40.00 ! VOQ(IP,VOQR(IVTYPE)) = WLMEAN ! ELSE ! VOQ(IP,VOQR(IVTYPE)) = OVEXCV(IVTYPE) ! ENDIF ! ! -------------------------------------------------------------------- / ! 4. Compute pars ! NPO = -1 ! DO IP=0, NPI ! HS = 4. * SQRT ( SUME(IP) * DDIR ) 40.PART !Removed factor (* 1/PI2) CHECK!! Check whether DDIR is defined in Rads! IF ( HS .LT. HSPMIN ) CYCLE ! XL = 1./(1+FRINTF) - 1. XH = FRINTF 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 = SPCSIG(IFPMAX(IP)) * ( 1. + 0.5 * ( XL2*EH - XH2*EL ) & / SIGN ( MAX(ABS(DENOM),1.E-15) , DENOM ) ) !> CALL WAVNU1 ( SIGP, DEPTH, WNP, CGP ) 40.PART !Temporaliry deactivated - replace which WLEN computation acc. to SWAN ! IF ( NPO .GE. DIMXP ) GOTO 2000 NPO = NPO + 1 IF (IP.GT.0) PMAP(NPO) = IP XP(1,NPO) = HS XP(2,NPO) = PI2 / SIGP !> XP(3,NPO) = PI2 / WNP IF ( SUMEK(IP) .GT. 0 ) THEN 40.PART XP(3,NPO) = PI2 * (SUME(IP) / SUMEK(IP)) ** (1./OUTPAR(3)) 40.PART ELSE 40.PART XP(3,NPO) = 0. 40.PART END IF 40.PART XP(4,NPO) = MOD( 630.-ATAN2(SUMEY(IP),SUMEX(IP))/DEGRAD , 360. ) 40.PART !This can be improved by using DD instead of 1/DEGRAD. EXTENT to include CART settings!! XP(5,NPO) = (1/DEGRAD) * SQRT ( MAX ( 0. , 2. * ( 1. - SQRT ( & MAX(0.,(SUMEX(IP)**2+SUMEY(IP)**2)/SUME(IP)**2) ) ) ) ) ! WRITE(PRINTF,*) 'In PTMEAN: IP, SUME (IP) =',IP, SUME (IP) 40.PART ! WRITE(PRINTF,*) 'In PTMEAN: IP, SUMEW(IP) =',IP, SUMEW(IP) 40.PART XP(6,NPO) = SUMEW(IP) / SUME(IP) IF ( XP(3,NPO) .GT. 0. ) THEN 40.PART XP(7,NPO) = HS / XP(3,NPO) 40.PART ! Is this only correct if OUTPAR(3) = 1 ?? ELSE 40.PART XP(7,NPO) = 0. 40.PART END IF 40.PART ! WRITE(PRINTF,*) 'In PTMEAN: NPO, XP(1,NPO) =',NPO, XP(1,NPO) 40.PART ! WRITE(PRINTF,*) 'In PTMEAN: NPO, XP(2,NPO) =',NPO, XP(2,NPO) 40.PART ! WRITE(PRINTF,*) 'In PTMEAN: NPO, XP(3,NPO) =',NPO, XP(3,NPO) 40.PART ! WRITE(PRINTF,*) 'In PTMEAN: NPO, XP(4,NPO) =',NPO, XP(4,NPO) 40.PART ! WRITE(PRINTF,*) 'In PTMEAN: NPO, XP(5,NPO) =',NPO, XP(5,NPO) 40.PART ! WRITE(PRINTF,*) 'In PTMEAN: NPO, XP(6,NPO) =',NPO, XP(6,NPO) 40.PART ! END DO ! RETURN ! ! Escape locations read errors --------------------------------------- * ! 2000 CONTINUE !> IF ( IAPROC .EQ. NAPERR ) WRITE (NDSE,1000) NPO+1 40.PART !Replace with SWAN error message WRITE (PRINTF,1000) NPO+1 40.PART !Replace with SWAN error message RETURN ! ! Formats ! 1000 FORMAT (/' *** WAVEWATCH III ERROR IN PTMEAN :'/ 40.PART !Replace with SWAN error message & ' XP ARRAY TOO SMALL AT PARTITION',I6/) !/ !/ End of PTMEAN ----------------------------------------------------- / !/ END SUBROUTINE PTMEAN !/ !/ End of module SWPARTMD -------------------------------------------- / !/ END MODULE SWPARTMD