Skip to content

Commit

Permalink
Bugfix memory leak in ST6 source terms. Ticket NOAA-EMC#1340
Browse files Browse the repository at this point in the history 
It appears that changing the file extension
from .ftn in v6.07 to .F90 in v7 triggered a memory leak in
Intel's Fortran compilers. The issue occurs when mapping
from a dynamically allocated array to a static array.
GCC gfortran on the other hand will make sure that allocated
arrays will be freed on exit of a subroutine.

The leak is about 250B per grid point and subroutine call.
This translates to 5MB per time step on a 2.5M proint grid.

Updates
-------
 [x] Change array type from static to allocatalbe.
 [x] Explicitely add DEALLOCATE statements.
 [x] Move arrays with constantes to MODULE level.
 [x] Add checks for allocation status.
  • Loading branch information
@stefanzieger
stefanzieger committed Jan 16, 2025
1 parent 4c3973a commit 631ec0f
Show file tree
Hide file tree
Showing 2 changed files with 199 additions and 43 deletions.
206 changes: 168 additions & 38 deletions model/src/w3src6md.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,10 @@ MODULE W3SRC6MD
!/
PUBLIC :: W3SPR6, W3SIN6, W3SDS6
PRIVATE :: LFACTOR, TAUWINDS, IRANGE
PRIVATE
INTEGER :: MK10Hz
INTEGER, ALLOCATABLE, SAVE :: IKN(:)
REAL, ALLOCATABLE, SAVE :: IK10Hz(:), SIG10Hz(:), DSII10Hz(:)
CONTAINS
!/ ------------------------------------------------------------------- /

Expand Down Expand Up @@ -397,7 +401,7 @@ SUBROUTINE W3SIN6 (A, CG, WN2, UABS, USTAR, USDIR, CD, DAIR, &
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
#endif
INTEGER :: IK, ITH, IKN(NK)
INTEGER :: IK, ITH
REAL :: COSU, SINU, UPROXY
REAL, DIMENSION(NSPEC) :: CG2, ECOS2, ESIN2, DSII2
REAL, DIMENSION(NK) :: DSII, SIG, WN
Expand Down Expand Up @@ -435,8 +439,14 @@ SUBROUTINE W3SIN6 (A, CG, WN2, UABS, USTAR, USDIR, CD, DAIR, &
ECOS2 = ECOS(1:NSPEC) ! Only indices from 1 to NSPEC
ESIN2 = ESIN(1:NSPEC) ! are requested.
!
IKN = IRANGE(1,NSPEC,NTH) ! Index vector for elements of 1 ... NK
! ! such that e.g. SIG(1:NK) = SIG2(IKN).
IF (.NOT.ALLOCATED(IKN)) THEN
IKN = IRANGE(1,NSPEC,NTH) ! Index vector for elements of 1 ... NK
! ! such that e.g. SIG(1:NK) = SIG2(IKN).
IF (SIZE(IKN).NE.NK) THEN
WRITE(NDSE,601) SIZE(IKN), NK
CALL EXTCDE(2)
END IF
END IF
DSII2 = DDEN2 / DTH / SIG2 ! Frequency bandwidths (int.) (rad)
DSII = DSII2(IKN)
SIG = SIG2(IKN)
Expand Down Expand Up @@ -514,6 +524,9 @@ SUBROUTINE W3SIN6 (A, CG, WN2, UABS, USTAR, USDIR, CD, DAIR, &
END IF
!
!/
601 FORMAT (' *** WAVEWATCH III ERROR IN W3SIN6 : '/ &
' INDEX ARRAY DIM MISMATCH SIZE(IKN)=',I4,' .NE. NK=',I4)
!/
!/ End of W3SIN6 ----------------------------------------------------- /
!/
END SUBROUTINE W3SIN6
Expand Down Expand Up @@ -648,7 +661,7 @@ SUBROUTINE W3SDS6 (A, CG, WN, S, D)
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
#endif
INTEGER :: IK, ITH, IKN(NK)
INTEGER :: IK, ITH
REAL :: FREQ(NK) ! frequencies [Hz]
REAL :: DFII(NK) ! frequency bandwiths [Hz]
REAL :: ANAR(NK) ! directional narrowness
Expand All @@ -673,9 +686,16 @@ SUBROUTINE W3SDS6 (A, CG, WN, S, D)
#endif
!
!/ 0) --- Initialize essential parameters ---------------------------- /
IKN = IRANGE(1,NSPEC,NTH) ! Index vector for elements of 1,
! ! 2,..., NK such that for example
! ! SIG(1:NK) = SIG2(IKN).
IF (.NOT.ALLOCATED(IKN)) THEN
IKN = IRANGE(1,NSPEC,NTH) ! Index vector for elements of 1,
! ! 2,..., NK such that for example
! ! SIG(1:NK) = SIG2(IKN).
IF (SIZE(IKN).NE.NK) THEN
WRITE(NDSE,601) SIZE(IKN), NK
CALL EXTCDE(2)
END IF
END IF
!
FREQ = SIG2(IKN)/TPI
ANAR = 1.0
BNT = 0.035**2
Expand Down Expand Up @@ -735,6 +755,8 @@ SUBROUTINE W3SDS6 (A, CG, WN, S, D)
270 FORMAT (' TEST W3SDS6 : ',A,'(',A,')',':',70E11.3)
271 FORMAT (' TEST W3SDS6 : Total SDS =',E13.5)
#endif
601 FORMAT (' *** WAVEWATCH III ERROR IN W3SDS6 : '/ &
' INDEX ARRAY DIM MISMATCH SIZE(IKN)=',I4,' .NE. NK=',I4)
!/
!/ End of W3SDS6 ----------------------------------------------------- /
!/
Expand All @@ -745,8 +767,8 @@ END SUBROUTINE W3SDS6
!>
!> @brief Numerical approximation for the reduction factor.
!>
!> @details Numerical approximation for the reduction factor LFACTOR(f) to
!> reduce energy in the high-frequency part of the resolved part
!> @details Numerical approximation for the reduction factor LFACTOR(f)
!> to reduce energy in the high-frequency part of the resolved part
!> of the spectrum to meet the constraint on total stress (TAU).
!> The constraint is TAU <= TAU_TOT (TAU_TOT = TAU_WAV + TAU_VIS),
!> thus the wind input is reduced to match our constraint.
Expand Down Expand Up @@ -877,14 +899,13 @@ SUBROUTINE LFACTOR(S, CINV, U10, USTAR, USDIR, SIG, DSII, &
INTEGER :: IK, NK10Hz, SIGN_NEW, SIGN_OLD
!
REAL :: ECOS2(NSPEC), ESIN2(NSPEC)
REAL, ALLOCATABLE :: IK10Hz(:), LF10Hz(:), SIG10Hz(:), CINV10Hz(:)
REAL, ALLOCATABLE :: LF10Hz(:), UCINV10Hz(:), CINV10Hz(:)
REAL, ALLOCATABLE :: SDENS10Hz(:), SDENSX10Hz(:), SDENSY10Hz(:)
REAL, ALLOCATABLE :: DSII10Hz(:), UCINV10Hz(:)
REAL :: TAU_TOT, TAU, TAU_VIS, TAU_WAV
REAL :: TAUVX, TAUVY, TAUX, TAUY
REAL :: TAU_NND, TAU_INIT(2)
REAL :: UPROXY, RTAU, DRTAU, ERR
LOGICAL :: OVERSHOT
LOGICAL :: OVERSHOT, FLGSET10Hz
CHARACTER(LEN=23) :: IDTIME
!
!/ ------------------------------------------------------------------- /
Expand All @@ -895,16 +916,29 @@ SUBROUTINE LFACTOR(S, CINV, U10, USTAR, USDIR, SIG, DSII, &
!/ 0) --- Find the number of frequencies required to extend arrays
!/ up to f=10Hz and allocate arrays --------------------------- /
!/ ALOG is the same as LOG
FLGSET10Hz = .FALSE.
NK10Hz = CEILING(ALOG(FRQMAX/(SIG(1)/TPI))/ALOG(XFR))+1
NK10Hz = MAX(NK,NK10Hz)
!
ALLOCATE(IK10Hz(NK10Hz))
IK10Hz = REAL( IRANGE(1,NK10Hz,1) )
IF (.NOT.ALLOCATED(IK10Hz)) THEN
IK10Hz = RRANGE(1,NK10Hz,1)
MK10Hz = NK10Hz
IF (SIZE(IK10Hz).NE.NK10Hz) THEN
WRITE(NDSE,601) SIZE(IK10Hz), NK10Hz
CALL EXTCDE(2)
END IF
END IF
IF (.NOT.ALLOCATED(SIG10Hz)) THEN
ALLOCATE(SIG10Hz(NK10Hz))
FLGSET10Hz = .TRUE.
END IF
IF (.NOT.ALLOCATED(DSII10Hz)) THEN
ALLOCATE(DSII10Hz(NK10Hz))
FLGSET10Hz = .TRUE.
END IF
!
ALLOCATE(SIG10Hz(NK10Hz))
ALLOCATE(CINV10Hz(NK10Hz))
ALLOCATE(DSII10Hz(NK10Hz))
ALLOCATE(LF10Hz(NK10Hz))
ALLOCATE(CINV10Hz(NK10Hz))
ALLOCATE(SDENS10Hz(NK10Hz))
ALLOCATE(SDENSX10Hz(NK10Hz))
ALLOCATE(SDENSY10Hz(NK10Hz))
Expand All @@ -913,29 +947,39 @@ SUBROUTINE LFACTOR(S, CINV, U10, USTAR, USDIR, SIG, DSII, &
ECOS2 = ECOS(1:NSPEC)
ESIN2 = ESIN(1:NSPEC)
!
!/ --- Check array size --- /
IF (MK10Hz.NE.NK10Hz) THEN
WRITE(NDSE,602) MK10Hz, NK10Hz
CALL EXTCDE(3)
END IF
!
!/ 1) --- Either extrapolate arrays up to 10Hz or use discrete spectral
! grid per se. Limit the constraint to the positive part of the
! wind input only. ---------------------------------------------- /
IF (NK .LT. NK10Hz) THEN
SDENS10Hz(1:NK) = SUM(S,1) * DTH
SDENSX10Hz(1:NK) = SUM(MAX(0.,S)*RESHAPE(ECOS2,(/NTH,NK/)),1) * DTH
SDENSY10Hz(1:NK) = SUM(MAX(0.,S)*RESHAPE(ESIN2,(/NTH,NK/)),1) * DTH
SIG10Hz = SIG(1)*XFR**(IK10Hz-1.0)
IF (FLGSET10Hz) THEN
SIG10Hz = SIG(1)*XFR**(IK10Hz-1.0)
DSII10Hz = 0.5 * SIG10Hz * (XFR-1.0/XFR)
! The first and last frequency bin:
DSII10Hz(1) = 0.5 * SIG10Hz(1) * (XFR-1.0)
DSII10Hz(NK10Hz) = 0.5 * SIG10Hz(NK10Hz) * (XFR-1.0) / XFR
END IF
!
CINV10Hz(1:NK) = CINV
CINV10Hz(NK+1:NK10Hz) = SIG10Hz(NK+1:NK10Hz)*0.101978 ! 1/c=σ/g
DSII10Hz = 0.5 * SIG10Hz * (XFR-1.0/XFR)
! The first and last frequency bin:
DSII10Hz(1) = 0.5 * SIG10Hz(1) * (XFR-1.0)
DSII10Hz(NK10Hz) = 0.5 * SIG10Hz(NK10Hz) * (XFR-1.0) / XFR
!
! --- Spectral slope for S_IN(F) is proportional to F**(-2) ------ /
SDENS10Hz(NK+1:NK10Hz) = SDENS10Hz(NK) * (SIG10Hz(NK)/SIG10Hz(NK+1:NK10Hz))**2
SDENSX10Hz(NK+1:NK10Hz) = SDENSX10Hz(NK) * (SIG10Hz(NK)/SIG10Hz(NK+1:NK10Hz))**2
SDENSY10hz(NK+1:NK10Hz) = SDENSY10Hz(NK) * (SIG10Hz(NK)/SIG10Hz(NK+1:NK10Hz))**2
ELSE
SIG10Hz = SIG
CINV10Hz = CINV
DSII10Hz = DSII
IF (FLGSET10Hz) THEN
SIG10Hz = SIG
DSII10Hz = DSII
END IF
CINV10Hz = CINV
SDENS10Hz(1:NK) = SUM(S,1) * DTH
SDENSX10Hz(1:NK) = SUM(MAX(0.,S)*RESHAPE(ECOS2,(/NTH,NK/)),1) * DTH
SDENSY10Hz(1:NK) = SUM(MAX(0.,S)*RESHAPE(ESIN2,(/NTH,NK/)),1) * DTH
Expand Down Expand Up @@ -1037,9 +1081,13 @@ SUBROUTINE LFACTOR(S, CINV, U10, USTAR, USDIR, SIG, DSII, &
274 FORMAT (' TEST W3SIN6 : Total ',A,' =', E13.5 )
280 FORMAT (' WARNING LFACTOR (TIME,U10,TAU,TAU_TOT,ERR,TAUW_XY,' &
'TAUV_XY,TAU_SCALAR): ',A,F6.1,2F7.4,E10.3,4F7.4,F7.3 )
601 FORMAT (' *** WAVEWATCH III ERROR IN W3SIN6 LFACTOR: '/ &
' INDEX ARRAY DIM MISMATCH SIZE(IKN)=',I4,' .NE. NK=',I4)
602 FORMAT (' *** WAVEWATCH III ERROR IN W3SIN6 LFACTOR: '/ &
' ARRAY SIZE MISMATCH MK=',I4,' .NE. NK=',I4)
!
DEALLOCATE(IK10Hz,SIG10Hz,CINV10Hz,DSII10Hz,LF10Hz)
DEALLOCATE(SDENS10Hz,SDENSX10Hz,SDENSY10Hz,UCINV10Hz)
DEALLOCATE(LF10Hz,CINV10Hz,UCINV10Hz)
DEALLOCATE(SDENS10Hz,SDENSX10Hz,SDENSY10Hz)
!/
END SUBROUTINE LFACTOR
!/ ------------------------------------------------------------------- /
Expand Down Expand Up @@ -1117,6 +1165,8 @@ SUBROUTINE TAU_WAVE_ATMOS(S, CINV, SIG, DSII, TAUNWX, TAUNWY )
!/
USE CONSTANTS, ONLY: GRAV, TPI
USE W3GDATMD, ONLY: NK, NTH, NSPEC, DTH, XFR, ECOS, ESIN
USE W3ODATMD, ONLY: NDSE
USE W3SERVMD, ONLY: EXTCDE
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
Expand All @@ -1137,9 +1187,9 @@ SUBROUTINE TAU_WAVE_ATMOS(S, CINV, SIG, DSII, TAUNWX, TAUNWY )
INTEGER :: NK10Hz
!
REAL :: ECOS2(NSPEC), ESIN2(NSPEC)
REAL, ALLOCATABLE :: IK10Hz(:), SIG10Hz(:), CINV10Hz(:)
REAL, ALLOCATABLE :: SDENSX10Hz(:), SDENSY10Hz(:)
REAL, ALLOCATABLE :: DSII10Hz(:), UCINV10Hz(:)
REAL, ALLOCATABLE :: CINV10Hz(:), UCINV10Hz(:)
LOGICAL :: FLGSET10Hz
!
!/ ------------------------------------------------------------------- /
#ifdef W3_S
Expand All @@ -1148,29 +1198,54 @@ SUBROUTINE TAU_WAVE_ATMOS(S, CINV, SIG, DSII, TAUNWX, TAUNWY )
!
!/ 0) --- Find the number of frequencies required to extend arrays
!/ up to f=10Hz and allocate arrays --------------------------- /
FLGSET10Hz = .FALSE.
NK10Hz = CEILING(ALOG(FRQMAX/(SIG(1)/TPI))/ALOG(XFR))+1
NK10Hz = MAX(NK,NK10Hz)
!
ALLOCATE(IK10Hz(NK10Hz))
IK10Hz = REAL( IRANGE(1,NK10Hz,1) )
IF (.NOT.ALLOCATED(IK10Hz)) THEN
IK10Hz = RRANGE(1,NK10Hz,1)
MK10Hz = NK10hz
IF (SIZE(IK10Hz).NE.NK10Hz) THEN
WRITE(NDSE,601) SIZE(IK10Hz), NK10Hz
CALL EXTCDE(2)
END IF
END IF
IF (.NOT.ALLOCATED(SIG10Hz)) THEN
ALLOCATE(SIG10Hz(NK10Hz))
FLGSET10Hz = .TRUE.
END IF
IF (.NOT.ALLOCATED(DSII10Hz)) THEN
ALLOCATE(DSII10Hz(NK10Hz))
FLGSET10Hz = .TRUE.
END IF
!
ALLOCATE(SIG10Hz(NK10Hz))
ALLOCATE(CINV10Hz(NK10Hz))
ALLOCATE(DSII10Hz(NK10Hz))
ALLOCATE(SDENSX10Hz(NK10Hz))
ALLOCATE(SDENSY10Hz(NK10Hz))
ALLOCATE(UCINV10Hz(NK10Hz))
!
ECOS2 = ECOS(1:NSPEC)
ESIN2 = ESIN(1:NSPEC)
!
!/ --- Check array size --- /
IF (MK10Hz.NE.NK10Hz) THEN
WRITE(NDSE,602) MK10Hz, NK10Hz
CALL EXTCDE(3)
END IF
!/ 1) --- Either extrapolate arrays up to 10Hz or use discrete spectral
! grid per se. Limit the constraint to the positive part of the
! wind input only. ---------------------------------------------- /
IF (NK .LT. NK10Hz) THEN
SDENSX10Hz(1:NK) = SUM(ABS(MIN(0.,S))*RESHAPE(ECOS2,(/NTH,NK/)),1) * DTH
SDENSY10Hz(1:NK) = SUM(ABS(MIN(0.,S))*RESHAPE(ESIN2,(/NTH,NK/)),1) * DTH
SIG10Hz = SIG(1)*XFR**(IK10Hz-1.0)
IF (FLGSET10Hz) THEN
SIG10Hz = SIG(1)*XFR**(IK10Hz-1.0)
DSII10Hz = 0.5 * SIG10Hz * (XFR-1.0/XFR)
! The first and last frequency bin:
DSII10Hz(1) = 0.5 * SIG10Hz(1) * (XFR-1.0)
DSII10Hz(NK10Hz) = 0.5 * SIG10Hz(NK10Hz) * (XFR-1.0) / XFR
END IF
!
CINV10Hz(1:NK) = CINV
CINV10Hz(NK+1:NK10Hz) = SIG10Hz(NK+1:NK10Hz)*0.101978
DSII10Hz = 0.5 * SIG10Hz * (XFR-1.0/XFR)
Expand All @@ -1182,9 +1257,11 @@ SUBROUTINE TAU_WAVE_ATMOS(S, CINV, SIG, DSII, TAUNWX, TAUNWY )
SDENSX10Hz(NK+1:NK10Hz) = SDENSX10Hz(NK) * (SIG10Hz(NK)/SIG10Hz(NK+1:NK10Hz))**2
SDENSY10hz(NK+1:NK10Hz) = SDENSY10Hz(NK) * (SIG10Hz(NK)/SIG10Hz(NK+1:NK10Hz))**2
ELSE
SIG10Hz = SIG
IF (FLGSET10Hz) THEN
SIG10Hz = SIG
DSII10Hz = DSII
END IF
CINV10Hz = CINV
DSII10Hz = DSII
SDENSX10Hz(1:NK) = SUM(ABS(MIN(0.,S))*RESHAPE(ECOS2,(/NTH,NK/)),1) * DTH
SDENSY10Hz(1:NK) = SUM(ABS(MIN(0.,S))*RESHAPE(ESIN2,(/NTH,NK/)),1) * DTH
END IF
Expand All @@ -1194,6 +1271,14 @@ SUBROUTINE TAU_WAVE_ATMOS(S, CINV, SIG, DSII, TAUNWX, TAUNWY )
TAUNWX = TAUWINDS(SDENSX10Hz,CINV10Hz,DSII10Hz) ! x-component
TAUNWY = TAUWINDS(SDENSY10Hz,CINV10Hz,DSII10Hz) ! y-component
!/
DEALLOCATE(CINV10Hz,UCINV10Hz)
DEALLOCATE(SDENSX10Hz,SDENSY10Hz)
!/
601 FORMAT (' *** WAVEWATCH III ERROR IN W3SIN6 TAU_WAVE_ATMOS : '/ &
' INDEX ARRAY DIM MISMATCH SIZE(IKN)=',I4,' .NE. NK=',I4)
602 FORMAT (' *** WAVEWATCH III ERROR IN W3SIN6 TAU_WAV_ATMOS : '/ &
' ARRAY SIZE MISMATCH MK=',I4,' .NE. NK=',I4)
!/
END SUBROUTINE TAU_WAVE_ATMOS
!/ ------------------------------------------------------------------- /
!/
Expand All @@ -1206,7 +1291,7 @@ END SUBROUTINE TAU_WAVE_ATMOS
!> @param X0
!> @param X1
!> @param DX
!> @returns IX
!> @returns IX (I.A.)
!>
!> @author S. Zieger
!> @date 15-Feb-2011
Expand Down Expand Up @@ -1243,6 +1328,51 @@ FUNCTION IRANGE(X0,X1,DX) RESULT(IX)
END FUNCTION IRANGE
!/ ------------------------------------------------------------------- /
!/
!>
!> @brief Generate a sequence of linear-spaced integer numbers.
!>
!> @details Used for instance array addressing (indexing).
!>
!> @param X0
!> @param X1
!> @param DX
!> @returns IX (R.A.)
!>
!> @author S. Zieger
!> @date 15-Feb-2011
!>
FUNCTION RRANGE(X0,X1,DX) RESULT(RX)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | S. Zieger |
!/ | FORTRAN 90 |
!/ | Last update : 15-Feb-2011 |
!/ +-----------------------------------+
!/
!/ 15-Feb-2011 : Origination ( version 4.04 )
!/ (S. Zieger)
!/
! 1. Purpose :
! Generate a sequence of linear-spaced numbers.
! Used for instance array addressing (indexing).
!
!/
IMPLICIT NONE
INTEGER, INTENT(IN) :: X0, X1, DX
REAL, ALLOCATABLE :: RX(:)
INTEGER :: N
INTEGER :: I
!
N = INT(REAL(X1-X0)/REAL(DX))+1
ALLOCATE(RX(N))
DO I = 1, N
RX(I) = REAL(X0+ (I-1)*DX)
END DO
!/
END FUNCTION RRANGE
!/ ------------------------------------------------------------------- /
!/

!>
!> @brief Wind stress (tau) computation from wind-momentum-input.
Expand Down
Loading

0 comments on commit 631ec0f

Please sign in to comment.