diff --git a/model/bin/switch_IC4Numerics b/model/bin/switch_IC4Numerics new file mode 100644 index 000000000..32b0925c5 --- /dev/null +++ b/model/bin/switch_IC4Numerics @@ -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 diff --git a/model/src/w3srcemd.F90 b/model/src/w3srcemd.F90 index e90ba88eb..5ed1a857b 100644 --- a/model/src/w3srcemd.F90 +++ b/model/src/w3srcemd.F90 @@ -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 ) @@ -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 @@ -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 @@ -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) @@ -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) ) @@ -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 @@ -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 @@ -2011,6 +2044,11 @@ 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 ! @@ -2018,6 +2056,9 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, & 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 @@ -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 @@ -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. @@ -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 @@ -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 @@ -2069,8 +2113,10 @@ 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 @@ -2078,6 +2124,7 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, & #ifdef W3_IS1 CALL W3SIS1 ( SPEC, ICE, VSIR ) #endif +#ifndef W3_IC4_ACCURATE_NUMERICS SPEC2 = SPEC ! TAUICE(:) = 0. @@ -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 @@ -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 ATT=EXP(ICE*VDIC(IS)*DTG) #endif #ifdef W3_IC5 @@ -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 ! @@ -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 ! @@ -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 ! ! ! - - - - - - - - - - - - - - - - - - - - - -