diff --git a/model/bin/switch_E3SM_wavice b/model/bin/switch_E3SM_wavice new file mode 100644 index 0000000000..88b91fd77a --- /dev/null +++ b/model/bin/switch_E3SM_wavice @@ -0,0 +1 @@ +F90 NOGRB NC4 DIST MPI PR3 UQ FLX0 LN1 ST4 STAB0 NL1 BT1 DB1 MLIM TR0 BS0 IC4 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/w3fld1md.F90 b/model/src/w3fld1md.F90 index 970e207069..a7af8a179f 100644 --- a/model/src/w3fld1md.F90 +++ b/model/src/w3fld1md.F90 @@ -1134,7 +1134,11 @@ SUBROUTINE APPENDTAIL(INSPC, WN2, NKT, KA1, KA2, KA3, WNDDIR,SAT) ENDDO AVG=SUM(NORMSPC)/MAX(REAL(NTH),1.) DO T=1, NTH - INSPC(K,T) = SAT * NORMSPC(T)/TPI/(WN2(K)**3.0)/AVG + IF (AVG /= 0.0) THEN + INSPC(K,T)=BT(K)*INSPC(K,T)/TPI/(WN2(K)**3.0)/AVG + ELSE + INSPC(K,T)=0.0 + ENDIF ENDDO ENDDO DO T=1, NTH diff --git a/model/src/w3gdatmd.F90 b/model/src/w3gdatmd.F90 index 9331a70e33..126d5300c6 100644 --- a/model/src/w3gdatmd.F90 +++ b/model/src/w3gdatmd.F90 @@ -680,6 +680,7 @@ MODULE W3GDATMD LOGICAL :: GINIT, FLDRY, FLCX, FLCY, FLCTH, FLCK, FLSOU, IICEDISP,& IICESMOOTH + LOGICAL :: IC_NUMERICS LOGICAL :: FLAGLL LOGICAL :: CMPRTRCK LOGICAL, POINTER :: FLAGST(:) @@ -1182,6 +1183,7 @@ MODULE W3GDATMD LOGICAL, POINTER :: GINIT, FLDRY, FLCX, FLCY, FLCTH, FLCK, FLSOU, IICEDISP,& IICESMOOTH + LOGICAL, POINTER :: IC_NUMERICS LOGICAL, POINTER :: FLAGLL LOGICAL, POINTER :: CMPRTRCK LOGICAL, POINTER :: FLAGST(:) @@ -2359,7 +2361,8 @@ SUBROUTINE W3SETG ( IMOD, NDSE, NDST ) FLSOU => GRIDS(IMOD)%FLSOU IICEDISP => GRIDS(IMOD)%IICEDISP IICESMOOTH => GRIDS(IMOD)%IICESMOOTH -! + IC_NUMERICS => GRIDS(IMOD)%IC_NUMERICS + ! GNAME => GRIDS(IMOD)%GNAME FILEXT => GRIDS(IMOD)%FILEXT TRIGP => GRIDS(IMOD)%TRIGP diff --git a/model/src/w3gridmd.F90 b/model/src/w3gridmd.F90 index 967cc25726..b22adaa678 100644 --- a/model/src/w3gridmd.F90 +++ b/model/src/w3gridmd.F90 @@ -808,6 +808,7 @@ MODULE W3GRIDMD ! REAL(8) :: GSHIFT ! see notes in WMGHGH LOGICAL :: FLC, ICEDISP, TRCKCMPR + LOGICAL :: ICNUMERICS INTEGER :: PTM ! Partitioning method REAL :: PTFC ! Part. cut off freq (for method 5) REAL :: AIRCMIN, AIRGB @@ -1101,7 +1102,7 @@ MODULE W3GRIDMD STDX, STDY, STDT, ICEHMIN, ICEHINIT, ICEDISP, & ICESLN, ICEWIND, ICESNL, ICESDS, ICEHFAC, & ICEHDISP, ICEDDISP, ICEFDISP, CALTYPE, & - TRCKCMPR, PTM, PTFC, BTBET + TRCKCMPR, PTM, PTFC, BTBET, ICNUMERICS NAMELIST /OUTS/ P2SF, I1P2SF, I2P2SF, & US3D, I1US3D, I2US3D, & USSP, IUSSP, STK_WN, & @@ -2744,6 +2745,7 @@ SUBROUTINE W3GRID(MDS) STDY = -1. STDT = -1. ICEDISP = .FALSE. + ICNUMERICS=.FALSE. CALTYPE = 'standard' ! Variables for 3D array output E3D=0 @@ -3017,6 +3019,7 @@ SUBROUTINE W3GRID(MDS) IICEHDISP = ICEHDISP IICEDDISP = ICEDDISP IICEFDISP = ICEFDISP + IC_NUMERICS=ICNUMERICS PMOVE = MAX ( 0. , PMOVE ) PFMOVE = PMOVE ! @@ -3404,7 +3407,7 @@ SUBROUTINE W3GRID(MDS) ICEHINIT, ICEDISP, ICEHDISP, & ICESLN, ICEWIND, ICESNL, ICESDS, & ICEDDISP,ICEFDISP, CALTYPE, TRCKCMPR, & - BTBETA + BTBETA,ICNUMERICS ELSE WRITE (NDSO,2966) CICE0, CICEN, LICE, PMOVE, XSEED, FLAGTR, & XP, XR, XFILT, IHMAX, HSPMIN, WSMULT, & @@ -3414,7 +3417,7 @@ SUBROUTINE W3GRID(MDS) ICEHINIT, ICEDISP, ICEHDISP, & ICESLN, ICEWIND, ICESNL, ICESDS, & ICEDDISP, ICEFDISP, CALTYPE, TRCKCMPR,& - BTBETA + BTBETA,ICNUMERICS END IF ! #ifdef W3_FLD1 diff --git a/model/src/w3iogrmd.F90 b/model/src/w3iogrmd.F90 index 2775bcfb43..3c1ff36f10 100644 --- a/model/src/w3iogrmd.F90 +++ b/model/src/w3iogrmd.F90 @@ -787,7 +787,7 @@ SUBROUTINE W3IOGR ( INXOUT, NDSM, IMOD, FEXT ) FLCK, FLSOU, FLBPI, FLBPO, CLATS, CLATIS, CTHG0S, & STEXU, STEYU, STEDU, IICEHMIN, IICEHINIT, IICEDISP, & ICESCALES(1:4), CALTYPE, CMPRTRCK, IICEHFAC, IICEHDISP,& - IICEDDISP, IICEFDISP, BTBETA, & + IICEDDISP, IICEFDISP, BTBETA, IC_NUMERICS, & AAIRCMIN, AAIRGB WRITE(NDSM)GRIDSHIFT @@ -980,7 +980,7 @@ SUBROUTINE W3IOGR ( INXOUT, NDSM, IMOD, FEXT ) FLCTH, FLCK, FLSOU, FLBPI, FLBPO, CLATS, CLATIS, & CTHG0S, STEXU, STEYU, STEDU, IICEHMIN, IICEHINIT, & IICEDISP, ICESCALES(1:4), CALTYPE, CMPRTRCK, IICEHFAC, & - IICEDDISP, IICEHDISP, IICEFDISP, BTBETA, & + IICEDDISP, IICEHDISP, IICEFDISP, BTBETA, IC_NUMERICS, & AAIRCMIN, AAIRGB #ifdef W3_DEBUGIOGR WRITE(740+IAPROC,*) 'W3IOGR, step 7.14' diff --git a/model/src/w3sic4md.F90 b/model/src/w3sic4md.F90 index 8661b42af8..1597be9d8a 100644 --- a/model/src/w3sic4md.F90 +++ b/model/src/w3sic4md.F90 @@ -322,7 +322,8 @@ SUBROUTINE W3SIC4 (A, DEPTH, CG, IX, IY, S, D) REAL, ALLOCATABLE :: ALPHA(:) ! exponential decay rate for energy REAL, ALLOCATABLE :: MARG1(:), MARG2(:) ! Arguments for M2 REAL, ALLOCATABLE :: KARG1(:), KARG2(:), KARG3(:) !Arguments for M3 - + REAL :: x1,x2,x3,x1sqr,x2sqr,x3sqr !Arguments for M10 + REAL :: perfour,amhb,bmhb !Arguments for M10 !/ !/ ------------------------------------------------------------------- / !/ @@ -513,6 +514,43 @@ SUBROUTINE W3SIC4 (A, DEPTH, CG, IX, IY, S, D) ALPHA(IK) = 0.2*(FREQ**2.13)*HICE END DO WN_I= 0.5 * ALPHA + + CASE (8) + ! Cubic fit to Meylan, Horvat & Bitz 2021 + ! ICECOEF1 is thickness + ! ICECOEF5 is floe size + ! TPI/SIG is period + x3=min(ICECOEF1,3.5) ! limit thickness to 3.5 m + x3=max(x3,0.1) ! limit thickness >0.1 m since I make fit below + x2=min(ICECOEF5*0.5,100.0) ! convert dia to radius, limit to 100m + x2=max(2.5,x2) + x2sqr=x2*x2 + x3sqr=x3*x3 + amhb = 2.12e-3 + bmhb = 4.59e-2 + + DO IK=1, NK + x1=TPI/SIG(IK) ! period + x1sqr=x1*x1 + KARG1(ik)=-0.26982 + 1.5043*x3 - 0.70112*x3sqr + 0.011037*x2 + & + (-0.0073178)*x2*x3 + 0.00036604*x2*x3sqr + & + (-0.00045789)*x2sqr + 1.8034e-05*x2sqr*x3 + & + (-0.7246)*x1 + 0.12068*x1*x3 + & + (-0.0051311)*x1*x3sqr + 0.0059241*x1*x2 + & + 0.00010771*x1*x2*x3 - 1.0171e-05*x1*x2sqr + & + 0.0035412*x1sqr - 0.0031893*x1sqr*x3 + & + (-0.00010791)*x1sqr*x2 + & + 0.00031073*x1**3 + 1.5996e-06*x2**3 + 0.090994*x3**3 + KARG1(IK)=min(KARG1(IK),0.0) + ALPHA(IK) = 10.0**KARG1(IK) + perfour=x1sqr*x1sqr + if ((x1.gt.5.0) .and. (x1.lt.20.0)) then + ALPHA(IK) = ALPHA(IK) + amhb/x1sqr+bmhb/perfour + else if (x1.gt.20.0) then + ALPHA(IK) = amhb/x1sqr+bmhb/perfour + endif + WN_I(IK) = ALPHA(IK) * 0.5 + end do CASE DEFAULT WN_I = ICECOEF1 !Default to IC1: Uniform in k diff --git a/model/src/w3srcemd.F90 b/model/src/w3srcemd.F90 index 475f2715f4..1eb5ee0468 100644 --- a/model/src/w3srcemd.F90 +++ b/model/src/w3srcemd.F90 @@ -366,6 +366,7 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, & XFC, XFLT, XREL, XFT, FXFM, FXPM, DDEN, & FTE, FTF, FHMAX, ECOS, ESIN, IICEDISP, & ICESCALES, IICESMOOTH + USE W3GDATMD, ONLY: IC_NUMERICS USE W3GDATMD, ONLY: FSSOURCE, optionCall USE W3GDATMD, ONLY: B_JGS_NLEVEL, B_JGS_SOURCE_NONLINEAR #ifdef W3_REF1 @@ -1227,6 +1228,32 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, & CALL UOST_SRCTRMCOMPUTE(IX, IY, SPEC, CG1, DT, & U10ABS, U10DIR, VSUO, VDUO) #endif +! Sea Ice Source Terms if IC_NUMERICS namelist flag = True + IF (IC_NUMERICS) THEN +#ifdef W3_IC1 + IF (ICE .GT. 0) CALL W3SIC1 ( SPEC,DEPTH, CG1, IX, IY, VSIC, VDIC ) +#endif +#ifdef W3_IS2 + IF (ICE .GT. 0) CALL W3SIS2 ( SPEC, DEPTH, ICE, ICEH, ICEF, ICEDMAX, IX, IY, & + VSIR, VDIR, VDIR2, WN1, CG1, WN_R, CG_ICE, R ) +#endif +#ifdef W3_IC2 + IF (ICE .GT. 0) CALL W3SIC2 ( SPEC, DEPTH, ICEH, ICEF, CG1, WN1,& + IX, IY, VSIC, VDIC, WN_R, CG_ICE, ALPHA_LIU, R) +#endif +#ifdef W3_IC3 + IF (ICE .GT. 0) CALL W3SIC3 ( SPEC,DEPTH, CG1, WN1, IX, IY, VSIC, VDIC ) +#endif +#ifdef W3_IC4 + IF (ICE .GT. 0) CALL W3SIC4 ( SPEC,DEPTH, CG1, IX, IY, VSIC, VDIC ) +#endif +#ifdef W3_IC5 + IF (ICE .GT. 0) CALL W3SIC5 ( SPEC,DEPTH, CG1, WN1, IX, IY, VSIC, VDIC ) +#endif +#ifdef W3_IS1 + IF (ICE .GT. 0) CALL W3SIS1 ( SPEC, ICE, VSIR ) +#endif + ENDIF ! ! 2.g Dump training data if necessary ! @@ -1288,6 +1315,12 @@ 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) + IF(IC_NUMERICS) THEN +#if defined(W3_IC1) || defined(W3_IC2) || defined(W3_IC3) || defined(W3_IC4) || defined(W3_IC5) + VSIC(1:NSPECH) = ICE * VSIC(1:NSPECH) ! (see Rogers et al 2016) + VDIC(1:NSPECH) = ICE * VDIC(1:NSPECH) +#endif + ENDIF END IF ! VS = 0 @@ -1307,6 +1340,11 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, & #ifdef W3_UOST VS(IS) = VS(IS) + VSUO(IS) #endif + IF ( IC_NUMERICS .AND. ICE.GT.0. ) THEN +#if defined(W3_IC1) || defined(W3_IC2) || defined(W3_IC3) || defined(W3_IC4) || defined(W3_IC5) + VS(IS) = VS(IS) + VSIC(IS) +#endif + ENDIF VD(IS) = VDIN(IS) + VDNL(IS) & + VDDS(IS) + VDBT(IS) #ifdef W3_ST6 @@ -1321,6 +1359,11 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, & #ifdef W3_UOST VD(IS) = VD(IS) + VDUO(IS) #endif + IF ( IC_NUMERICS .AND. ICE.GT.0. ) THEN +#if defined(W3_IC1) || defined(W3_IC2) || defined(W3_IC3) || defined(W3_IC4) || defined(W3_IC5) + VD(IS) = VD(IS) + VDIC(IS) +#endif + ENDIF DAMAX = MIN ( DAM(IS) , MAX ( XREL*SPECINIT(IS) , AFILT ) ) AFAC = 1. / MAX( 1.E-10 , ABS(VS(IS)/DAMAX) ) #ifdef W3_NL5 @@ -1603,6 +1646,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 + IF ( IC_NUMERICS .AND. ICE.GT.0 ) THEN +#if defined(W3_IC1) || defined(W3_IC2) || defined(W3_IC3) || defined(W3_IC4) || defined(W3_IC5) + 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))) +#endif + ENDIF IF (VSIN(IS).GT.0.) WHITECAP(3) = WHITECAP(3) + SPEC(IS) * FACTOR HSTOT = HSTOT + SPEC(IS) * FACTOR END DO @@ -1867,6 +1918,13 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, & ! TAUOX=(GRAV*MWXFINISH+TAUWIX-TAUBBL(1))/DTG TAUOY=(GRAV*MWYFINISH+TAUWIY-TAUBBL(2))/DTG + IF (IC_NUMERICS) THEN +#if defined(W3_IC1) || defined(W3_IC2) || defined(W3_IC3) || defined(W3_IC4) || defined(W3_IC5) + TAUICE(:)=TAUICE(:)/DTG + TAUOX = TAUOX - TAUICE(1) + TAUOY = TAUOY - TAUICE(2) +#endif + ENDIF TAUWIX=TAUWIX/DTG TAUWIY=TAUWIY/DTG TAUWNX=TAUWNX/DTG @@ -1881,6 +1939,11 @@ 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 + IF (IC_NUMERICS) THEN +#if defined(W3_IC1) || defined(W3_IC2) || defined(W3_IC3) || defined(W3_IC4) || defined(W3_IC5) + PHICE =-1.*DWAT*GRAV*PHICE/DTG +#endif + ENDIF ! ! 10.1 Adds ice scattering and dissipation: implicit integration---------------- * ! INFLAGS2(4) is true if ice concentration was ever read during @@ -1893,7 +1956,7 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, & #endif IF ( INFLAGS2(4).AND.ICE.GT.0 ) THEN - + IF (.NOT. IC_NUMERICS ) THEN IF (IICEDISP) THEN ICECOEF2 = 1E-6 CALL LIU_FORWARD_DISPERSION (ICEH,ICECOEF2,DEPTH, & @@ -2021,6 +2084,7 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, & END DO PHICE =-1.*DWAT*GRAV*PHICE /DTG TAUICE(:)=TAUICE(:)/DTG + ENDIF ! end if IC_NUMERICS ELSE #ifdef W3_IS2 IF (IS2PARS(10).LT.0.5) THEN