From 5bce5e68e5989062100044e9ddec3395b14c62ec Mon Sep 17 00:00:00 2001 From: Benoit Pouliot Date: Thu, 18 Feb 2021 16:03:08 +0000 Subject: [PATCH 1/5] Bugfix/negative_peak_frequency: Add IS2 flag to avoid uninitialized calculation in ww3_outp --- model/src/ww3_outp.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/model/src/ww3_outp.F90 b/model/src/ww3_outp.F90 index 5f2a65bf3..6fb92b4db 100644 --- a/model/src/ww3_outp.F90 +++ b/model/src/ww3_outp.F90 @@ -1439,10 +1439,10 @@ SUBROUTINE W3EXPO RHOAIR = MAX ( 0. , DAIRO(J)) #endif CDIR = MOD ( 270. - CDO(J)*RADE , 360. ) - ICEDMAX = MAX ( 0., ICEFO(J)) - ICEF = ICEDMAX - ICETHICK = MAX (0., ICEHO(J)) - ICECON = MAX (0., ICEO(J)) +!/IS2 ICEDMAX = MAX ( 0., ICEFO(J)) +!/IS2 ICEF = ICEDMAX +!/IS2 ICETHICK = MAX (0., ICEHO(J)) +!/IS2 ICECON = MAX (0., ICEO(J)) ! #ifdef W3_STAB2 STAB0 = ZWIND * GRAV / 273. From c42cb854a9ef27a8d7be1b82c86792a64948b781 Mon Sep 17 00:00:00 2001 From: Benoit Pouliot Date: Thu, 18 Feb 2021 16:22:29 +0000 Subject: [PATCH 2/5] Bugfix/negative_peak_frequency: Fix peak frequency calculation in w3iogomd, ww3_outp and ww3_ounp (#304) --- model/src/w3iogomd.F90 | 43 ++++++++++++++++++++++++++---------------- model/src/ww3_ounp.F90 | 20 +++++++++++++++----- model/src/ww3_outp.F90 | 32 +++++++++++++++++++++---------- 3 files changed, 64 insertions(+), 31 deletions(-) diff --git a/model/src/w3iogomd.F90 b/model/src/w3iogomd.F90 index 82f0e4a97..cc3c6fa5e 100644 --- a/model/src/w3iogomd.F90 +++ b/model/src/w3iogomd.F90 @@ -57,9 +57,11 @@ MODULE W3IOGOMD !/ (Roberto Padilla-Hernandez & J.H. Alves) !/ 03-Nov-2020 : Factored out NAME matching into ( version 7.12 ) !/ seperate subroutine. (C. Bunney) -!/ 15-Jan-2020 : Added TP output based on exsiting ( version 7.12 ) +!/ 15-Jan-2021 : Added TP output based on exsiting ( version 7.12 ) !/ FP internal field. (C. Bunney) !/ 22-Mar-2021 : Add extra coupling fields as output ( version 7.13 ) +!/ 21-Jul-2022 : Correct FP0 calc for peak energy in ( version 7.14 ) +!/ min/max freq band (B. Pouliot, CMC) !/ !/ Copyright 2009-2014 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights @@ -2013,7 +2015,7 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 ) CALL INIT_GET_ISEA(ISEA, JSEA) EC (JSEA) = EBD(NK,JSEA) FP0 (JSEA) = UNDEF - IKP0(JSEA) = 0 + IKP0(JSEA) = NK THP0(JSEA) = UNDEF #ifdef W3_ST0 FP1 (JSEA) = UNDEF @@ -2062,7 +2064,7 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 ) ! ! 4.b Discrete peak frequencies ! - DO IK=NK-1, 2, -1 + DO IK=NK-1, 1, -1 ! #ifdef W3_OMPG !$OMP PARALLEL DO PRIVATE(JSEA,ISEA) @@ -2075,7 +2077,7 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 ) IKP0(JSEA) = IK END IF #ifdef W3_ST1 - IF ( IKP1(JSEA).EQ.0 & + IF ( IKP1(JSEA).EQ.0 .AND. IK .GT. 1 & .AND. EBD(IK-1,JSEA).LT.EBD(IK,JSEA) & .AND. EBD(IK-1,JSEA).LT.EBD(IK+1,JSEA) & .AND. SIG(IK).GT.FXPMC/UST(ISEA) & @@ -2083,7 +2085,7 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 ) IKP1(JSEA) = IK #endif #ifdef W3_ST3 - IF ( IKP1(JSEA).EQ.0 & + IF ( IKP1(JSEA).EQ.0 .AND. IK .GT. 1 & .AND. EBD(IK-1,JSEA).LT.EBD(IK,JSEA) & .AND. EBD(IK-1,JSEA).LT.EBD(IK+1,JSEA) & .AND. SIG(IK).GT.FXPMC/MAX(1.E-4,UST(ISEA)) & @@ -2091,7 +2093,7 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 ) IKP1(JSEA) = IK #endif #ifdef W3_ST4 - IF ( IKP1(JSEA).EQ.0 & + IF ( IKP1(JSEA).EQ.0 .AND. IK .GT. 1 & .AND. EBD(IK-1,JSEA).LT.EBD(IK,JSEA) & .AND. EBD(IK-1,JSEA).LT.EBD(IK+1,JSEA) & .AND. SIG(IK).GT.FXPMC/MAX(1.E-4,UST(ISEA)) & @@ -2112,7 +2114,7 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 ) ! DO JSEA=1, NSEAL CALL INIT_GET_ISEA(ISEA, JSEA) - IF ( IKP0(JSEA) .NE. 0 ) FP0(JSEA) = SIG(IKP0(JSEA)) * TPIINV + IF ( EC(JSEA) .GT. 0 ) FP0(JSEA) = SIG(IKP0(JSEA)) * TPIINV #ifdef W3_ST1 IF ( IKP1(JSEA) .NE. 0 ) FP1(JSEA) = SIG(IKP1(JSEA)) * TPIINV #endif @@ -2141,14 +2143,23 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 ) ! DO JSEA=1, NSEAL CALL INIT_GET_ISEA(ISEA, JSEA) - ILOW = MAX ( 1 , IKP0(JSEA)-1 ) - ICEN = MAX ( 1 , IKP0(JSEA) ) - IHGH = MIN ( NK , IKP0(JSEA)+1 ) - EL = EBD(ILOW,JSEA) - EBD(ICEN,JSEA) - EH = EBD(IHGH,JSEA) - EBD(ICEN,JSEA) - DENOM = XL*EH - XH*EL - FP0(JSEA) = FP0 (JSEA) * ( 1. + 0.5 * ( XL2*EH - XH2*EL ) & - / SIGN ( MAX(ABS(DENOM),1.E-15) , DENOM ) ) + IF ( EC(JSEA) .GT. 0 ) THEN + IF ( IKP0(JSEA) .EQ. 1 ) THEN + EL = - EBD(IKP0(JSEA), JSEA) + ELSE + EL = EBD(IKP0(JSEA)-1, JSEA) - EBD(IKP0(JSEA), JSEA) + END IF + + IF ( IKP0(JSEA) .EQ. NK ) THEN + EH = - EBD(IKP0(JSEA), JSEA) + ELSE + EH = EBD(IKP0(JSEA)+1, JSEA) - EBD(IKP0(JSEA), JSEA) + END IF + + DENOM = XL*EH - XH*EL + FP0(JSEA) = FP0 (JSEA) * ( 1. + 0.5 * ( XL2*EH - XH2*EL ) & + / SIGN ( MAX(ABS(DENOM),1.E-15) , DENOM ) ) + END IF #ifdef W3_ST1 ILOW = MAX ( 1 , IKP1(JSEA)-1 ) ICEN = MAX ( 1 , IKP1(JSEA) ) @@ -2208,7 +2219,7 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 ) ! DO JSEA=1, NSEAL CALL INIT_GET_ISEA(ISEA, JSEA) - IF (IKP0(JSEA).NE.0) THEN + IF ( EC(JSEA) .GT. 0 ) THEN ETX(JSEA) = ETX(JSEA) + A(ITH,IKP0(JSEA),JSEA)*ECOS(ITH) ETY(JSEA) = ETY(JSEA) + A(ITH,IKP0(JSEA),JSEA)*ESIN(ITH) END IF diff --git a/model/src/ww3_ounp.F90 b/model/src/ww3_ounp.F90 index 6bbc45b45..44896a6fe 100644 --- a/model/src/ww3_ounp.F90 +++ b/model/src/ww3_ounp.F90 @@ -75,6 +75,8 @@ PROGRAM W3OUNP !/ 19-Jul-2021 : Momentum and air density support ( version 7.14 ) !/ 06-Sep-2021 : scale factor on spectra output ( version 7.12 ) !/ 05-Jan-2022 : Added TIMESPLIT=0 (nodate) support ( version 7.14 ) +!/ 21-Jul-2022 : Correct FP0 calc for peak energy in ( version 7.14 ) +!/ min/max freq band (B. Pouliot, CMC) !/ !/ Copyright 2009 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights @@ -1595,7 +1597,7 @@ SUBROUTINE W3EXNC(I,NCID,NREQ,INDREQ,ORDER) !/ Local parameters !/ - INTEGER :: J, J1, I1, I2, ISP, IKM, IKL, IKH, & + INTEGER :: J, J1, I1, I2, ISP, IKM, & ITH, IK, ITT, NPART, IX, IY, ISEA INTEGER :: CURDATE(8), REFDATE(8) #ifdef W3_S @@ -2031,10 +2033,18 @@ SUBROUTINE W3EXNC(I,NCID,NREQ,INDREQ,ORDER) END IF END DO ! - IKL = MAX ( 1 , IKM-1 ) - IKH = MIN ( NK , IKM+1 ) - EL = E1(IKL) - E1(IKM) - EH = E1(IKH) - E1(IKM) + IF ( IKM .EQ. 1 ) THEN + EL = - E1(IKM) + ELSE + EL = E1(IKM-1) - E1(IKM) + END IF + + IF ( IKM .EQ. NK ) THEN + EH = - E1(IKM) + ELSE + EH = E1(IKM+1) - E1(IKM) + END IF + DENOM = XL*EH - XH*EL ! IF ( HSIG .GE. HSMIN ) THEN diff --git a/model/src/ww3_outp.F90 b/model/src/ww3_outp.F90 index 6fb92b4db..6689955ae 100644 --- a/model/src/ww3_outp.F90 +++ b/model/src/ww3_outp.F90 @@ -60,6 +60,8 @@ PROGRAM W3OUTP !/ 27-Jun-2017 : Expanding WMO table to 2 digits JHA ( version 6.02 ) !/ 18-Aug-2018 : S_{ice} IC5 (Q. Liu) ( version 6.06 ) !/ 19-Jul-2021 : Momentum and air density support ( version 7.14 ) +!/ 21-Jul-2022 : Correct FP0 calc for peak energy in ( version 7.14 ) +!/ min/max freq band (B. Pouliot, CMC) !/ !/ Copyright 2009-2014 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights @@ -1224,7 +1226,7 @@ SUBROUTINE W3EXPO !/ ------------------------------------------------------------------- / !/ Local parameters !/ - INTEGER :: J, I1, I2, ISP, IKM, IKL, IKH, ITH, & + INTEGER :: J, I1, I2, ISP, IKM, ITH, & IK, IH, IM, IS, IYR, IMTH, IDY, ITT, & I, NPART, IP, IX, IY, ISEA INTEGER, SAVE :: IPASS = 0 @@ -1439,10 +1441,12 @@ SUBROUTINE W3EXPO RHOAIR = MAX ( 0. , DAIRO(J)) #endif CDIR = MOD ( 270. - CDO(J)*RADE , 360. ) -!/IS2 ICEDMAX = MAX ( 0., ICEFO(J)) -!/IS2 ICEF = ICEDMAX -!/IS2 ICETHICK = MAX (0., ICEHO(J)) -!/IS2 ICECON = MAX (0., ICEO(J)) +#ifdef W3_IS2 + ICEDMAX = MAX ( 0., ICEFO(J)) + ICEF = ICEDMAX + ICETHICK = MAX (0., ICEHO(J)) + ICECON = MAX (0., ICEO(J)) +#endif ! #ifdef W3_STAB2 STAB0 = ZWIND * GRAV / 273. @@ -1568,11 +1572,19 @@ SUBROUTINE W3EXPO IKM = IK END IF END DO -! - IKL = MAX ( 1 , IKM-1 ) - IKH = MIN ( NK , IKM+1 ) - EL = E1(IKL) - E1(IKM) - EH = E1(IKH) - E1(IKM) + + IF ( IKM .EQ. 1 ) THEN + EL = - E1(IKM) + ELSE + EL = E1(IKM-1) - E1(IKM) + END IF + + IF ( IKM .EQ. NK ) THEN + EH = - E1(IKM) + ELSE + EH = E1(IKM+1) - E1(IKM) + END IF + DENOM = XL*EH - XH*EL ! IF ( HSIG .GE. HSMIN ) THEN From 8d8fef2b9891ca48ee06208f6f6abc4cb31c248f Mon Sep 17 00:00:00 2001 From: Benoit Pouliot Date: Mon, 8 Aug 2022 19:29:13 +0000 Subject: [PATCH 3/5] Forbid highest frequency for FP --- model/src/w3iogomd.F90 | 12 +++++------- model/src/ww3_ounp.F90 | 18 +++++++----------- model/src/ww3_outp.F90 | 24 ++++++++++++------------ 3 files changed, 24 insertions(+), 30 deletions(-) diff --git a/model/src/w3iogomd.F90 b/model/src/w3iogomd.F90 index cc3c6fa5e..823a7b99f 100644 --- a/model/src/w3iogomd.F90 +++ b/model/src/w3iogomd.F90 @@ -2114,7 +2114,9 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 ) ! DO JSEA=1, NSEAL CALL INIT_GET_ISEA(ISEA, JSEA) - IF ( EC(JSEA) .GT. 0 ) FP0(JSEA) = SIG(IKP0(JSEA)) * TPIINV + IF ( EC(JSEA) .GT. 0 .AND. IKP0(JSEA) .NE. NK ) THEN + FP0(JSEA) = SIG(IKP0(JSEA)) * TPIINV + END IF #ifdef W3_ST1 IF ( IKP1(JSEA) .NE. 0 ) FP1(JSEA) = SIG(IKP1(JSEA)) * TPIINV #endif @@ -2143,18 +2145,14 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 ) ! DO JSEA=1, NSEAL CALL INIT_GET_ISEA(ISEA, JSEA) - IF ( EC(JSEA) .GT. 0 ) THEN + IF ( EC(JSEA) .GT. 0 .AND. IKP0(JSEA) .NE. NK ) THEN IF ( IKP0(JSEA) .EQ. 1 ) THEN EL = - EBD(IKP0(JSEA), JSEA) ELSE EL = EBD(IKP0(JSEA)-1, JSEA) - EBD(IKP0(JSEA), JSEA) END IF - IF ( IKP0(JSEA) .EQ. NK ) THEN - EH = - EBD(IKP0(JSEA), JSEA) - ELSE - EH = EBD(IKP0(JSEA)+1, JSEA) - EBD(IKP0(JSEA), JSEA) - END IF + EH = EBD(IKP0(JSEA)+1, JSEA) - EBD(IKP0(JSEA), JSEA) DENOM = XL*EH - XH*EL FP0(JSEA) = FP0 (JSEA) * ( 1. + 0.5 * ( XL2*EH - XH2*EL ) & diff --git a/model/src/ww3_ounp.F90 b/model/src/ww3_ounp.F90 index 44896a6fe..794e31f05 100644 --- a/model/src/ww3_ounp.F90 +++ b/model/src/ww3_ounp.F90 @@ -2033,21 +2033,17 @@ SUBROUTINE W3EXNC(I,NCID,NREQ,INDREQ,ORDER) END IF END DO ! - IF ( IKM .EQ. 1 ) THEN - EL = - E1(IKM) - ELSE - EL = E1(IKM-1) - E1(IKM) - END IF + IF ( HSIG .GE. HSMIN .AND. IKM .NE. NK ) THEN + IF ( IKM .EQ. 1 ) THEN + EL = - E1(IKM) + ELSE + EL = E1(IKM-1) - E1(IKM) + END IF - IF ( IKM .EQ. NK ) THEN - EH = - E1(IKM) - ELSE EH = E1(IKM+1) - E1(IKM) - END IF - DENOM = XL*EH - XH*EL + DENOM = XL*EH - XH*EL ! - IF ( HSIG .GE. HSMIN ) THEN FP = SIG(IKM) * ( 1. + 0.5 * ( XL2*EH - XH2*EL ) & / SIGN ( MAX(ABS(DENOM),1.E-15) , DENOM ) ) THP = THBND(IKM) diff --git a/model/src/ww3_outp.F90 b/model/src/ww3_outp.F90 index 6689955ae..45135e08e 100644 --- a/model/src/ww3_outp.F90 +++ b/model/src/ww3_outp.F90 @@ -1573,21 +1573,21 @@ SUBROUTINE W3EXPO END IF END DO - IF ( IKM .EQ. 1 ) THEN - EL = - E1(IKM) - ELSE - EL = E1(IKM-1) - E1(IKM) - END IF + IF ( HSIG .GE. HSMIN .AND. IKM .NE. NK ) THEN + IF ( IKM .EQ. 1 ) THEN + EL = - E1(IKM) + ELSE + EL = E1(IKM-1) - E1(IKM) + END IF - IF ( IKM .EQ. NK ) THEN - EH = - E1(IKM) - ELSE - EH = E1(IKM+1) - E1(IKM) - END IF + IF ( IKM .EQ. NK ) THEN + EH = - E1(IKM) + ELSE + EH = E1(IKM+1) - E1(IKM) + END IF - DENOM = XL*EH - XH*EL + DENOM = XL*EH - XH*EL ! - IF ( HSIG .GE. HSMIN ) THEN FP = SIG(IKM) * ( 1. + 0.5 * ( XL2*EH - XH2*EL ) & / SIGN ( MAX(ABS(DENOM),1.E-15) , DENOM ) ) THP = THBND(IKM) From 2d8edc599660fe77160a02160b0ab62373d0ba99 Mon Sep 17 00:00:00 2001 From: Benoit Pouliot Date: Fri, 12 Aug 2022 13:30:17 +0000 Subject: [PATCH 4/5] Remove redundant conditions --- model/src/w3iogomd.F90 | 6 +++--- model/src/ww3_outp.F90 | 6 +----- 2 files changed, 4 insertions(+), 8 deletions(-) diff --git a/model/src/w3iogomd.F90 b/model/src/w3iogomd.F90 index 823a7b99f..6ceafe781 100644 --- a/model/src/w3iogomd.F90 +++ b/model/src/w3iogomd.F90 @@ -2114,7 +2114,7 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 ) ! DO JSEA=1, NSEAL CALL INIT_GET_ISEA(ISEA, JSEA) - IF ( EC(JSEA) .GT. 0 .AND. IKP0(JSEA) .NE. NK ) THEN + IF ( IKP0(JSEA) .NE. NK ) THEN FP0(JSEA) = SIG(IKP0(JSEA)) * TPIINV END IF #ifdef W3_ST1 @@ -2145,7 +2145,7 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 ) ! DO JSEA=1, NSEAL CALL INIT_GET_ISEA(ISEA, JSEA) - IF ( EC(JSEA) .GT. 0 .AND. IKP0(JSEA) .NE. NK ) THEN + IF ( IKP0(JSEA) .NE. NK ) THEN IF ( IKP0(JSEA) .EQ. 1 ) THEN EL = - EBD(IKP0(JSEA), JSEA) ELSE @@ -2217,7 +2217,7 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 ) ! DO JSEA=1, NSEAL CALL INIT_GET_ISEA(ISEA, JSEA) - IF ( EC(JSEA) .GT. 0 ) THEN + IF ( IKP0(JSEA) .NE. NK ) THEN ETX(JSEA) = ETX(JSEA) + A(ITH,IKP0(JSEA),JSEA)*ECOS(ITH) ETY(JSEA) = ETY(JSEA) + A(ITH,IKP0(JSEA),JSEA)*ESIN(ITH) END IF diff --git a/model/src/ww3_outp.F90 b/model/src/ww3_outp.F90 index 45135e08e..cf6a2b1f5 100644 --- a/model/src/ww3_outp.F90 +++ b/model/src/ww3_outp.F90 @@ -1580,11 +1580,7 @@ SUBROUTINE W3EXPO EL = E1(IKM-1) - E1(IKM) END IF - IF ( IKM .EQ. NK ) THEN - EH = - E1(IKM) - ELSE - EH = E1(IKM+1) - E1(IKM) - END IF + EH = E1(IKM+1) - E1(IKM) DENOM = XL*EH - XH*EL ! From 087e6ecf9d87d60804d84a87a9eb5d8424f24f4f Mon Sep 17 00:00:00 2001 From: Benoit Pouliot Date: Tue, 27 Sep 2022 14:13:56 +0000 Subject: [PATCH 5/5] Remove unused variables --- model/src/w3iogomd.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/model/src/w3iogomd.F90 b/model/src/w3iogomd.F90 index 1aa0e8696..cdc58ee19 100644 --- a/model/src/w3iogomd.F90 +++ b/model/src/w3iogomd.F90 @@ -1317,8 +1317,7 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 ) !/ INTEGER :: IK, ITH, JSEA, ISEA, IX, IY, & IKP0(NSEAL), NKH(NSEAL), & - ILOW, ICEN, IHGH, I, J, LKMS, HKMS, & - ITL + I, J, LKMS, HKMS, ITL #ifdef W3_S INTEGER, SAVE :: IENT = 0 #endif