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

Updated Numerical Solution for Wave Attenuation in sea ice #1294

Open
wants to merge 4 commits into
base: develop
Choose a base branch
from
Open

Updated Numerical Solution for Wave Attenuation in sea ice #1294

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
2b8187a
IC4 Accurate Wave Ice numerics
Jul 11, 2024
9d67e87
preprocessor directive must start with "W3_" "W3_IC4_ACCURATE_NUMERICS"
Jul 12, 2024
81513b0
add numerics switch file
erinethomas Aug 14, 2024
b703df9
fix to logic around CALL W3SIC4
erinethomas Aug 14, 2024
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
1 change: 1 addition & 0 deletions model/bin/switch_IC4Numerics
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
F90 NOGRB NC4 DIST MPI PR3 UQ FLX0 LN1 ST4 STAB0 NL1 BT1 DB1 MLIM TR0 BS0 IC4 IC4_ACCURATE_NUMERICS IS0 REF0 XX0 WNT1 WNX1 CRT1 CRX1 O0 O1 O2 O3 O4 O5 O6 O7 SCRIPNC SCRIP RTD RWND UOST
58 changes: 56 additions & 2 deletions model/src/w3srcemd.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -469,6 +469,9 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
! viscoelastic sea ice model (Wang and Shen 2010).
! !/IC4 Dissipation via interaction with ice as a function of freq.
! (empirical/parametric methods)
! !/IC4_ACCURATE_NUMERICS
! Correction for numerics for wave damping in sea ice.
! only tested for use with IC4M10 (Meylan et al 2021)
! !/IC5 Dissipation via interaction with ice according to a
! viscoelastic sea ice model (Mosig et al. 2015).
! !/DB0 No depth-limited breaking. ( Choose one )
Expand Down Expand Up @@ -1324,6 +1327,15 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
! UNRESOLVED OBSTACLES
CALL UOST_SRCTRMCOMPUTE(IX, IY, SPEC, CG1, DT, &
U10ABS, U10DIR, VSUO, VDUO)
#endif
!
! Sea Ice Damping Source Term
!
#ifdef W3_IC4
#ifdef W3_IC4_ACCURATE_NUMERICS
if (ICE.GT.0) CALL W3SIC4 ( SPEC,DEPTH, CG1, &
IX, IY, VSIC, VDIC )
#endif
#endif
!
! 2.g Dump training data if necessary
Expand Down Expand Up @@ -1384,6 +1396,10 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
VDIN(1:NSPECH) = ICESCALEIN * VDIN(1:NSPECH)
VSDS(1:NSPECH) = ICESCALEDS * VSDS(1:NSPECH)
VDDS(1:NSPECH) = ICESCALEDS * VDDS(1:NSPECH)
#ifdef W3_IC4_ACCURATE_NUMERICS
VSIC(1:NSPECH) = ICE * VSIC(1:NSPECH) ! (see Rogers et al 2016)
VDIC(1:NSPECH) = ICE * VDIC(1:NSPECH) ! **************
#endif
END IF

#ifdef W3_PDLIB
Expand Down Expand Up @@ -1422,6 +1438,9 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
#endif
#ifdef W3_UOST
VS(IS) = VS(IS) + VSUO(IS)
#endif
#ifdef W3_IC4_ACCURATE_NUMERICS
IF ( ICE.GT.0. ) VS(IS) = VS(IS) + VSIC(IS)
#endif
VD(IS) = VDIN(IS) + VDNL(IS) &
+ VDDS(IS) + VDBT(IS)
Expand All @@ -1436,6 +1455,9 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
#endif
#ifdef W3_UOST
VD(IS) = VD(IS) + VDUO(IS)
#endif
#ifdef W3_IC4_ACCURATE_NUMERICS
IF ( ICE.GT.0. ) VD(IS) = VD(IS) + VDIC(IS)
#endif
DAMAX = MIN ( DAM(IS) , MAX ( XREL*SPECINIT(IS) , AFILT ) )
AFAC = 1. / MAX( 1.E-10 , ABS(VS(IS)/DAMAX) )
Expand All @@ -1455,6 +1477,9 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
!
DT = MAX ( 0.5, DT ) ! The hardcoded min. dt is a problem for certain cases e.g. laborotary scale problems.
!
#ifdef W3_IC4_ACCURATE_NUMERICS
if (ICE.gt.0.01 .and. ICE.lt.0.95) DT=DTMIN ! always use a small timestep in ice
#endif
DTDYN = DTDYN + DT
#ifdef W3_T
DTRAW = DT
Expand Down Expand Up @@ -1743,6 +1768,14 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
/ MAX ( 1. , (1.-HDT*VDBT(IS))) ! semi-implict integration scheme
PHINL = PHINL + VSNL(IS)* DT * FACTOR &
/ MAX ( 1. , (1.-HDT*VDNL(IS))) ! semi-implict integration scheme
#ifdef W3_IC4_ACCURATE_NUMERICS
IF ( ICE.GT.0 ) THEN
PHICE = PHICE + VSIC(IS) * DT * FACTOR &
/ MAX ( 1. , (1.-HDT*VDIC(IS))) ! semi-implicit integration
TAUICE(:) = TAUICE(:) - FACTOR2*COSI(:)*VSIC(IS) * DT &
/ MAX ( 1. , (1.-HDT*VDIC(IS)))
END IF
#endif
IF (VSIN(IS).GT.0.) WHITECAP(3) = WHITECAP(3) + SPEC(IS) * FACTOR
HSTOT = HSTOT + SPEC(IS) * FACTOR
END DO
Expand Down Expand Up @@ -2011,13 +2044,21 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
TAUBBL(:)=TAUBBL(:)/DTG
TAUOCX=DAIR*COEF*COEF*USTAR*USTAR*COS(USTDIR) + DWAT*(TAUOX-TAUWIX)
TAUOCY=DAIR*COEF*COEF*USTAR*USTAR*SIN(USTDIR) + DWAT*(TAUOY-TAUWIY)
#ifdef W3_IC4_ACCURATE_NUMERICS
TAUICE(:)=TAUICE(:)/DTG
TAUOX = TAUOX - TAUICE(1)
TAUOY = TAUOY - TAUICE(2)
#endif
!
! Transformation in wave energy flux in W/m^2=kg / s^3
!
PHIOC =DWAT*GRAV*(EFINISH+PHIAW-PHIBBL)/DTG
PHIAW =DWAT*GRAV*PHIAW /DTG
PHINL =DWAT*GRAV*PHINL /DTG
PHIBBL=DWAT*GRAV*PHIBBL/DTG
#ifdef W3_IC4_ACCURATE_NUMERICS
PHICE =-1.*DWAT*GRAV*PHICE/DTG
#endif
!
! 10.1 Adds ice scattering and dissipation: implicit integration---------------- *
! INFLAGS2(4) is true if ice concentration was ever read during
Expand All @@ -2028,7 +2069,7 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
WRITE(740+IAPROC,*) '3 : sum(SPEC)=', sum(SPEC)
END IF
#endif

#ifndef W3_IC4_ACCURATE_NUMERICS
IF ( INFLAGS2(4).AND.ICE.GT.0 ) THEN

IF (IICEDISP) THEN
Expand All @@ -2037,6 +2078,7 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
SIG,WN_R,CG_ICE,ALPHA_LIU)
!
IF (IICESMOOTH) THEN
#endif
#ifdef W3_IS2
DO IK=1,NK
SMOOTH_ICEDISP=0.
Expand All @@ -2046,6 +2088,7 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
WN_R(IK)=WN1(IK)*(1-SMOOTH_ICEDISP)+WN_R(IK)*(SMOOTH_ICEDISP)
END DO
#endif
#ifndef W3_IC4_ACCURATE_NUMERICS
END IF
ELSE
WN_R=WN1
Expand All @@ -2054,6 +2097,7 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
!
R(:)=1 ! In case IC2 is defined but not IS2
!
#endif
#ifdef W3_IC1
CALL W3SIC1 ( SPEC,DEPTH, CG1, IX, IY, VSIC, VDIC )
#endif
Expand All @@ -2069,15 +2113,18 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
CALL W3SIC3 ( SPEC,DEPTH, CG1, WN1, IX, IY, VSIC, VDIC )
#endif
#ifdef W3_IC4
#ifndef W3_IC4_ACCURATE_NUMERICS
CALL W3SIC4 ( SPEC,DEPTH, CG1, IX, IY, VSIC, VDIC )
#endif
#endif
#ifdef W3_IC5
CALL W3SIC5 ( SPEC,DEPTH, CG1, WN1, IX, IY, VSIC, VDIC )
#endif
!
#ifdef W3_IS1
CALL W3SIS1 ( SPEC, ICE, VSIR )
#endif
#ifndef W3_IC4_ACCURATE_NUMERICS
SPEC2 = SPEC
!
TAUICE(:) = 0.
Expand All @@ -2088,6 +2135,7 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
! First part of ice term integration: dissipation part
!
ATT=1.
#endif
#ifdef W3_IC1
ATT=EXP(ICE*VDIC(IS)*DTG)
#endif
Expand All @@ -2097,7 +2145,7 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
#ifdef W3_IC3
ATT=EXP(ICE*VDIC(IS)*DTG)
#endif
#ifdef W3_IC4
#ifndef W3_IC4_ACCURATE_NUMERICS
Copy link
Collaborator

@sbanihash sbanihash Dec 7, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In a case of IC0 for example, where there is no IC1-IC5 defined, then VDIC has not been defined (look at line 744 in this same script) and the build will give this error: This name does not have a type, and must have an explicit type. [VDIC]
ATT=EXP(ICE*VDIC(IS)*DTG)
suggested changes is:
#ifdef WW3_IC4
#ifndef WW3_IC4_ACCURATE_NUMERICS
ATT = EXP(ICE * VDIC(IS) * DTG)
#endif
#endif

ATT=EXP(ICE*VDIC(IS)*DTG)
#endif
#ifdef W3_IC5
Expand All @@ -2115,7 +2163,9 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
ATT=ATT*EXP((ICE*VDIR(IS))*DTG)
END IF
#endif
#ifndef W3_IC4_ACCURATE_NUMERICS
SPEC(1+(IK-1)*NTH:NTH+(IK-1)*NTH) = ATT*SPEC2(1+(IK-1)*NTH:NTH+(IK-1)*NTH)
#endif
!
! Second part of ice term integration: scattering including re-distribution in directions
!
Expand All @@ -2142,6 +2192,7 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
END IF
END IF
#endif
#ifndef W3_IC4_ACCURATE_NUMERICS
!
! 10.2 Fluxes of energy and momentum due to ice effects
!
Expand All @@ -2158,12 +2209,15 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
PHICE =-1.*DWAT*GRAV*PHICE /DTG
TAUICE(:)=TAUICE(:)/DTG
ELSE
#endif
#ifdef W3_IS2
IF (IS2PARS(10).LT.0.5) THEN
ICEF = 0.
ENDIF
#endif
#ifndef W3_IC4_ACCURATE_NUMERICS
END IF
#endif
!
!
! - - - - - - - - - - - - - - - - - - - - - -
Expand Down