Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bugfix/negative peak frequency #741

Merged
merged 6 commits into from
Oct 3, 2022
Merged

Bugfix/negative peak frequency #741

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
5bce5e6
Bugfix/negative_peak_frequency: Add IS2 flag to avoid uninitialized c…
Feb 18, 2021
c42cb85
Bugfix/negative_peak_frequency: Fix peak frequency calculation in w3i…
Feb 18, 2021
8d8fef2
Forbid highest frequency for FP
benoitp-cmc Aug 8, 2022
2d8edc5
Remove redundant conditions
benoitp-cmc Aug 12, 2022
fd70804
Merge branch 'develop' into Bugfix/negative_peak_frequency
benoitp-cmc Sep 27, 2022
087e6ec
Remove unused variables
benoitp-cmc Sep 27, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 23 additions & 17 deletions model/src/w3iogomd.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,9 +69,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
Expand Down Expand Up @@ -1315,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
Expand Down Expand Up @@ -2051,7 +2052,7 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 )
DO JSEA=1, NSEAL
EC (JSEA) = EBD(NK,JSEA)
FP0 (JSEA) = UNDEF
IKP0(JSEA) = 0
IKP0(JSEA) = NK
THP0(JSEA) = UNDEF
END DO
!
Expand All @@ -2061,7 +2062,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)
Expand All @@ -2085,7 +2086,7 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 )
#endif
!
DO JSEA=1, NSEAL
IF ( IKP0(JSEA) .NE. 0 ) FP0(JSEA) = SIG(IKP0(JSEA)) * TPIINV
IF ( IKP0(JSEA) .NE. NK ) FP0(JSEA) = SIG(IKP0(JSEA)) * TPIINV
END DO
!
#ifdef W3_OMPG
Expand All @@ -2100,19 +2101,24 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 )
XH2 = XH**2
!
#ifdef W3_OMPG
!$OMP PARALLEL DO PRIVATE(JSEA,ILOW,ICEN,IHGH,EL,EH,DENOM)
!$OMP PARALLEL DO PRIVATE(JSEA,EL,EH,DENOM)
#endif
!
DO JSEA=1, NSEAL
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 ) )
END DO
IF ( 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

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 ) &
/ SIGN ( MAX(ABS(DENOM),1.E-15) , DENOM ) )
END IF
END DO
!
#ifdef W3_OMPG
!$OMP END PARALLEL DO
Expand Down Expand Up @@ -2140,7 +2146,7 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 )
#endif
!
DO JSEA=1, NSEAL
IF (IKP0(JSEA).NE.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
Expand Down
20 changes: 13 additions & 7 deletions model/src/ww3_ounp.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -2031,13 +2033,17 @@ 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)
DENOM = XL*EH - XH*EL
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

EH = E1(IKM+1) - E1(IKM)

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)
Expand Down
24 changes: 16 additions & 8 deletions model/src/ww3_outp.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -1439,10 +1441,12 @@ SUBROUTINE W3EXPO
RHOAIR = MAX ( 0. , DAIRO(J))
#endif
CDIR = MOD ( 270. - CDO(J)*RADE , 360. )
#ifdef W3_IS2
benoitp-cmc marked this conversation as resolved.
Show resolved Hide resolved
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.
Expand Down Expand Up @@ -1568,14 +1572,18 @@ SUBROUTINE W3EXPO
IKM = IK
END IF
END DO

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

EH = E1(IKM+1) - E1(IKM)

DENOM = XL*EH - XH*EL
!
IKL = MAX ( 1 , IKM-1 )
IKH = MIN ( NK , IKM+1 )
EL = E1(IKL) - E1(IKM)
EH = E1(IKH) - E1(IKM)
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)
Expand Down