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

WW3-MPAS-SI Coupling for E3SMv3 #1

Open
wants to merge 8 commits into
base: e3sm
Choose a base branch
from
Open

WW3-MPAS-SI Coupling for E3SMv3 #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
12c941a
add ic4metho8: Meylan et al 2021: wave attenuation due to floe scatte…
Dec 21, 2023
bf8c5fe
add wavice switch file (uses IC4M8)
Jan 23, 2024
16cef75
bug fx in IC4M8: WN_I needed to be multiplied by factor of .5
Mar 22, 2024
0c00190
Implementing Cecilia Bitz Numerics modifications for IC4. (see https:…
Jun 19, 2024
838fd81
remove v3 switch file
Jul 11, 2024
77be71d
updated comments
Jul 11, 2024
77170c5
separate the numerics fix into a separate 'switch' setting 'IC4NUM'
erinethomas Sep 17, 2024
cb5f176
bug fix to w3fld1md
erinethomas Sep 17, 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_E3SM_wavice
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 IS0 REF0 XX0 WNT1 WNX1 CRT1 CRX1 O0 O1 O2 O3 O4 O5 O6 O7 SCRIPNC SCRIP RTD RWND UOST
1 change: 1 addition & 0 deletions model/bin/switch_E3SM_waviceNUM
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 IC4NUM IS0 REF0 XX0 WNT1 WNX1 CRT1 CRX1 O0 O1 O2 O3 O4 O5 O6 O7 SCRIPNC SCRIP RTD RWND UOST
6 changes: 5 additions & 1 deletion model/ftn/w3fld1md.ftn
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1092,7 +1092,11 @@
DO K=KA1, KA2-1
AVG=SUM(INSPC(K,:))/MAX(REAL(NTH),1.)
DO T=1,NTH
INSPC(K,T)=BT(K)*INSPC(K,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
!-----------------------------------------------------------
Expand Down
39 changes: 39 additions & 0 deletions model/ftn/w3sic4md.ftn
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -310,6 +310,9 @@
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 !case 8
REAL :: perfour,amhb,bmhb !case 8

!/
!/ ------------------------------------------------------------------- /
Expand Down Expand Up @@ -493,6 +496,42 @@
ALPHA(IK) = 0.2*(FREQ**2.13)*HICE
END DO
WN_I= 0.5 * ALPHA
CASE (8)
!CC Bitz added option of 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) = 0.5 * ALPHA(IK)
end do

CASE DEFAULT
WN_I = ICECOEF1 !Default to IC1: Uniform in k
Expand Down
169 changes: 132 additions & 37 deletions model/ftn/w3srcemd.ftn
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -326,6 +326,8 @@
! viscoelastic sea ice model (Wang and Shen 2010).
! !/IC4 Dissipation via interaction with ice as a function of freq.
! (empirical/parametric methods)
! !/IC4NUM Dissipation via interaction with ice as a function of freq.
!
! !/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 @@ -408,6 +410,7 @@
!/IC2 USE W3SIC2MD
!/IC3 USE W3SIC3MD
!/IC4 USE W3SIC4MD
!/IC4NUM USE W3SIC4MD
!/IC5 USE W3SIC5MD
!/IS1 USE W3SIS1MD
!/IS2 USE W3SIS2MD
Expand Down Expand Up @@ -498,6 +501,7 @@
!/IC2 VSIC(NSPEC), VDIC(NSPEC), &
!/IC3 VSIC(NSPEC), VDIC(NSPEC), &
!/IC4 VSIC(NSPEC), VDIC(NSPEC), &
!/IC4NUM VSIC(NSPEC), VDIC(NSPEC), &
!/IC5 VSIC(NSPEC), VDIC(NSPEC), &
!/DB1 VSDB(NSPEC), VDDB(NSPEC), &
!/DBX VSDB(NSPEC), VDDB(NSPEC), &
Expand Down Expand Up @@ -575,6 +579,8 @@
!/IC3 VDIC = 0.
!/IC4 VSIC = 0.
!/IC4 VDIC = 0.
!/IC4NUM VSIC = 0.
!/IC4NUM VDIC = 0.
!/UOST VSUO = 0.
!/UOST VDUO = 0.
!/IC5 VSIC = 0.
Expand Down Expand Up @@ -872,6 +878,10 @@
!/UOST CALL UOST_SRCTRMCOMPUTE(IX, IY, SPEC, CG1, DT, &
!/UOST U10ABS, U10DIR, VSUO, VDUO)
!
! new location for dissipation due to ice if running with numerics fix
!/IC4NUM IF (ICE .GT. 0) CALL W3SIC4 ( SPEC,DEPTH, CG1, &
!/IC4NUM IX, IY, VSIC, VDIC )
!
! 2.f Additional sources.
!
!/XXX CALL W3SXXX
Expand Down Expand Up @@ -928,6 +938,8 @@
VDIN(1:NSPECH) = ICESCALEIN * VDIN(1:NSPECH)
VSDS(1:NSPECH) = ICESCALEDS * VSDS(1:NSPECH)
VDDS(1:NSPECH) = ICESCALEDS * VDDS(1:NSPECH)
!/IC4NUM VSIC(1:NSPECH) = ICE * VSIC(1:NSPECH) ! (see Rogers et al 2016)
!/IC4NUM VDIC(1:NSPECH) = ICE * VDIC(1:NSPECH) ! ^
END IF
!
VS = 0
Expand All @@ -941,16 +953,18 @@
!/BS1 VS(IS) = VS(IS) + VSBS(IS)
!/BSX VS(IS) = VS(IS) + VSBS(IS)
!/XXX VS(IS) = VS(IS) + VSXX(IS)
!/UOST VS(IS) = VS(IS) + VSUO(IS)
VD(IS) = VDIN(IS) + VDNL(IS) &
!/UOST VS(IS) = VS(IS) + VSUO(IS)
!/IC4NUM IF ( ICE.GT.0. ) VS(IS) = VS(IS) + VSIC(IS)
VD(IS) = VDIN(IS) + VDNL(IS) &
+ VDDS(IS) + VDBT(IS)
!/ST6 VD(IS) = VD(IS) + VDWL(IS)
!/TR1 VD(IS) = VD(IS) + VDTR(IS)
!/TRX VD(IS) = VD(IS) + VDTR(IS)
!/BS1 VD(IS) = VD(IS) + VDBS(IS)
!/BSX VD(IS) = VD(IS) + VDBS(IS)
!/XXX VD(IS) = VD(IS) + VDXX(IS)
!/UOST VD(IS) = VD(IS) + VDUO(IS)
!/UOST VD(IS) = VD(IS) + VDUO(IS)
!/IC4NUM IF ( ICE.GT.0. ) VD(IS) = VD(IS) + VDIC(IS)
DAMAX = MIN ( DAM(IS) , MAX ( XREL*SPECINIT(IS) , AFILT ) )
AFAC = 1. / MAX( 1.E-10 , ABS(VS(IS)/DAMAX) )
DT = MIN ( DT , AFAC / ( MAX ( 1.E-10, &
Expand All @@ -964,6 +978,7 @@
! IF (IX == DEBUG_NODE) WRITE(*,*) 'TIMINGS 1', DT
DT = MAX ( 0.5, DT ) ! Here we have a hardlimit, which is not too usefull, at least not as a fixed constant
!
!/IC4NUM IF (ICE .GT. 0.01 .and. ICE .LT. 0.95) DT=DTMIN ! always use a small timestep in ice
DTDYN = DTDYN + DT
!/T DTRAW = DT
IDT = 1 + INT ( 0.99*(DTG-DTTOT)/DT ) ! number of iterations
Expand Down Expand Up @@ -1139,6 +1154,12 @@
/ 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
!/IC4NUM IF ( ICE.GT.0 ) THEN
!/IC4NUM PHICE = PHICE + VSIC(IS) * DT * FACTOR &
!/IC4NUM / MAX ( 1. , (1.-HDT*VDIC(IS))) ! semi-implicit integration
!/IC4NUM TAUICE(:) = TAUICE(:) - FACTOR2*COSI(:)*VSIC(IS) * DT &
!/IC4NUM / MAX ( 1. , (1.-HDT*VDIC(IS)))
!/IC4NUM ENDIF
IF (VSIN(IS).GT.0.) WHITECAP(3) = WHITECAP(3) + SPEC(IS) * FACTOR
HSTOT = HSTOT + SPEC(IS) * FACTOR
END DO
Expand Down Expand Up @@ -1336,6 +1357,9 @@
!
TAUOX=(GRAV*MWXFINISH+TAUWIX-TAUBBL(1))/DTG
TAUOY=(GRAV*MWYFINISH+TAUWIY-TAUBBL(2))/DTG
!/IC4NUM TAUICE(:)=TAUICE(:)/DTG
!/IC4NUM TAUOX = TAUOX - TAUICE(1)
!/IC4NUM TAUOY = TAUOY - TAUICE(2)
TAUWIX=TAUWIX/DTG
TAUWIY=TAUWIY/DTG
TAUWNX=TAUWNX/DTG
Expand All @@ -1350,6 +1374,7 @@
PHIAW =DWAT*GRAV*PHIAW /DTG
PHINL =DWAT*GRAV*PHINL /DTG
PHIBBL=DWAT*GRAV*PHIBBL/DTG
!/IC4NUM PHICE =-1.*DWAT*GRAV*PHICE/DTG
!
! 10.1 Adds ice scattering and dissipation: implicit integration---------------- *
! INFLAGS2(4) is true if ice concentration was ever read during
Expand All @@ -1358,14 +1383,23 @@
!/DEBUGSRC IF (IX .eq. DEBUG_NODE) THEN
!/DEBUGSRC WRITE(740+IAPROC,*) '3 : sum(SPEC)=', sum(SPEC)
!/DEBUGSRC END IF

! The following section of code should only get used for any IC
! switch EXCEPT IC4NUM.
IF ( INFLAGS2(4).AND.ICE.GT.0 ) THEN

IF (IICEDISP) THEN
ICECOEF2 = 1E-6
CALL LIU_FORWARD_DISPERSION (ICEH,ICECOEF2,DEPTH, &
SIG,WN_R,CG_ICE,ALPHA_LIU)
!
!/IC1 ICECOEF2 = 1E-6
!/IC1 CALL LIU_FORWARD_DISPERSION (ICEH,ICECOEF2,DEPTH, &
!/IC1 SIG,WN_R,CG_ICE,ALPHA_LIU)
!/IC2 ICECOEF2 = 1E-6
!/IC2 CALL LIU_FORWARD_DISPERSION (ICEH,ICECOEF2,DEPTH, &
!/IC2 SIG,WN_R,CG_ICE,ALPHA_LIU)
!/IC3 ICECOEF2 = 1E-6
!/IC3 CALL LIU_FORWARD_DISPERSION (ICEH,ICECOEF2,DEPTH, &
!/IC3 SIG,WN_R,CG_ICE,ALPHA_LIU)
!/IC4 ICECOEF2 = 1E-6
!/IC4 CALL LIU_FORWARD_DISPERSION (ICEH,ICECOEF2,DEPTH, &
!/IC4 SIG,WN_R,CG_ICE,ALPHA_LIU)
IF (IICESMOOTH) THEN
!/IS2 DO IK=1,NK
!/IS2 SMOOTH_ICEDISP=0.
Expand All @@ -1374,13 +1408,22 @@
!/IS2 END IF
!/IS2 WN_R(IK)=WN1(IK)*(1-SMOOTH_ICEDISP)+WN_R(IK)*(SMOOTH_ICEDISP)
!/IS2 END DO
END IF
ELSE
WN_R=WN1
CG_ICE=CG1
END IF
!
R(:)=1 ! In case IC2 is defined but not IS2
END IF ! end icesmooth
ELSE ! iicedisp
!/IC1 WN_R=WN1
!/IC1 CG_ICE=CG1
!/IC2 WN_R=WN1
!/IC2 CG_ICE=CG1
!/IC3 WN_R=WN1
!/IC3 CG_ICE=CG1
!/IC4 WN_R=WN1
!/IC4 CG_ICE=CG1
END IF ! end icedisp
!
!/IC1 R(:)=1 ! In case IC2 is defined but not IS2
!/IC2 R(:)=1 ! In case IC2 is defined but not IS2
!/IC3 R(:)=1 ! In case IC2 is defined but not IS2
!/IC4 R(:)=1 ! In case IC2 is defined but not IS2
!
!/IC1 CALL W3SIC1 ( SPEC,DEPTH, CG1, IX, IY, VSIC, VDIC )
!/IS2 CALL W3SIS2 ( SPEC, DEPTH, ICE, ICEH, ICEF, ICEDMAX, IX, IY, &
Expand All @@ -1393,16 +1436,31 @@
!/IC5 CALL W3SIC5 ( SPEC,DEPTH, CG1, WN1, IX, IY, VSIC, VDIC )
!
!/IS1 CALL W3SIS1 ( SPEC, ICE, VSIR )
SPEC2 = SPEC
!
TAUICE(:) = 0.
PHICE = 0.
DO IK=1,NK
IS = 1+(IK-1)*NTH
!/IC1 SPEC2 = SPEC
!/IC2 SPEC2 = SPEC
!/IC3 SPEC2 = SPEC
!/IC4 SPEC2 = SPEC
!
!/IC1 TAUICE(:) = 0.
!/IC2 TAUICE(:) = 0.
!/IC3 TAUICE(:) = 0.
!/IC1 PHICE = 0.
!/IC2 PHICE = 0.
!/IC3 PHICE = 0.
!/IC4 PHICE = 0.
DO IK=1,NK
!/IC1 IS = 1+(IK-1)*NTH
!/IC2 IS = 1+(IK-1)*NTH
!/IC3 IS = 1+(IK-1)*NTH
!/IC4 IS = 1+(IK-1)*NTH
!
! First part of ice term integration: dissipation part
!
ATT=1.
!/IC1 ATT=1.
!/IC2 ATT=1.
!/IC3 ATT=1.
!/IC4 ATT=1.

!/IC1 ATT=EXP(ICE*VDIC(IS)*DTG)
!/IC2 ATT=EXP(ICE*VDIC(IS)*DTG)
!/IC3 ATT=EXP(ICE*VDIC(IS)*DTG)
Expand All @@ -1416,7 +1474,11 @@
!/IS2!
!/IS2 ATT=ATT*EXP((ICE*VDIR(IS))*DTG)
!/IS2 END IF
SPEC(1+(IK-1)*NTH:NTH+(IK-1)*NTH) = ATT*SPEC2(1+(IK-1)*NTH:NTH+(IK-1)*NTH)

!/IC1 SPEC(1+(IK-1)*NTH:NTH+(IK-1)*NTH) = ATT*SPEC2(1+(IK-1)*NTH:NTH+(IK-1)*NTH)
!/IC2 SPEC(1+(IK-1)*NTH:NTH+(IK-1)*NTH) = ATT*SPEC2(1+(IK-1)*NTH:NTH+(IK-1)*NTH)
!/IC3 SPEC(1+(IK-1)*NTH:NTH+(IK-1)*NTH) = ATT*SPEC2(1+(IK-1)*NTH:NTH+(IK-1)*NTH)
!/IC4 SPEC(1+(IK-1)*NTH:NTH+(IK-1)*NTH) = ATT*SPEC2(1+(IK-1)*NTH:NTH+(IK-1)*NTH)
!
! Second part of ice term integration: scattering including re-distribution in directions
!
Expand Down Expand Up @@ -1444,23 +1506,56 @@
!
! 10.2 Fluxes of energy and momentum due to ice effects
!
FACTOR = DDEN(IK)/CG1(IK) !Jacobian to get energy in band
FACTOR2= FACTOR*GRAV*WN1(IK)/SIG(IK) ! coefficient to get momentum
DO ITH = 1,NTH
IS = ITH+(IK-1)*NTH
PHICE = PHICE + (SPEC(IS)-SPEC2(IS)) * FACTOR
COSI(1)=ECOS(IS)
COSI(2)=ESIN(IS)
TAUICE(:) = TAUICE(:) - (SPEC(IS)-SPEC2(IS))*FACTOR2*COSI(:)
END DO
END DO
PHICE =-1.*DWAT*GRAV*PHICE /DTG
TAUICE(:)=TAUICE(:)/DTG
ELSE
!/IC1 FACTOR = DDEN(IK)/CG1(IK) !Jacobian to get energy in band
!/IC1 FACTOR2= FACTOR*GRAV*WN1(IK)/SIG(IK) ! coefficient to get momentum
!/IC1 DO ITH = 1,NTH
!/IC1 IS = ITH+(IK-1)*NTH
!/IC1 PHICE = PHICE + (SPEC(IS)-SPEC2(IS)) * FACTOR
!/IC1 COSI(1)=ECOS(IS)
!/IC1 COSI(2)=ESIN(IS)
!/IC1 TAUICE(:) = TAUICE(:) - (SPEC(IS)-SPEC2(IS))*FACTOR2*COSI(:)
!/IC1 END DO
!/IC2 FACTOR = DDEN(IK)/CG1(IK) !Jacobian to get energy in band
!/IC2 FACTOR2= FACTOR*GRAV*WN1(IK)/SIG(IK) ! coefficient to get momentum
!/IC2 DO ITH = 1,NTH
!/IC2 IS = ITH+(IK-1)*NTH
!/IC2 PHICE = PHICE + (SPEC(IS)-SPEC2(IS)) * FACTOR
!/IC2 COSI(1)=ECOS(IS)
!/IC2 COSI(2)=ESIN(IS)
!/IC2 TAUICE(:) = TAUICE(:) - (SPEC(IS)-SPEC2(IS))*FACTOR2*COSI(:)
!/IC2 END DO
!/IC3 FACTOR = DDEN(IK)/CG1(IK) !Jacobian to get energy in band
!/IC3 FACTOR2= FACTOR*GRAV*WN1(IK)/SIG(IK) ! coefficient to get momentum
!/IC3 DO ITH = 1,NTH
!/IC3 IS = ITH+(IK-1)*NTH
!/IC3 PHICE = PHICE + (SPEC(IS)-SPEC2(IS)) * FACTOR
!/IC3 COSI(1)=ECOS(IS)
!/IC3 COSI(2)=ESIN(IS)
!/IC3 TAUICE(:) = TAUICE(:) - (SPEC(IS)-SPEC2(IS))*FACTOR2*COSI(:)
!/IC3 END DO
!/IC4 FACTOR = DDEN(IK)/CG1(IK) !Jacobian to get energy in band
!/IC4 FACTOR2= FACTOR*GRAV*WN1(IK)/SIG(IK) ! coefficient to get momentum
!/IC4 DO ITH = 1,NTH
!/IC4 IS = ITH+(IK-1)*NTH
!/IC4 PHICE = PHICE + (SPEC(IS)-SPEC2(IS)) * FACTOR
!/IC4 COSI(1)=ECOS(IS)
!/IC4 COSI(2)=ESIN(IS)
!/IC4 TAUICE(:) = TAUICE(:) - (SPEC(IS)-SPEC2(IS))*FACTOR2*COSI(:)
!/IC4 END DO
END DO !ik
!/IC1 PHICE =-1.*DWAT*GRAV*PHICE /DTG
!/IC1 TAUICE(:)=TAUICE(:)/DTG
!/IC2 PHICE =-1.*DWAT*GRAV*PHICE /DTG
!/IC2 TAUICE(:)=TAUICE(:)/DTG
!/IC3 PHICE =-1.*DWAT*GRAV*PHICE /DTG
!/IC3 TAUICE(:)=TAUICE(:)/DTG
!/IC4 PHICE =-1.*DWAT*GRAV*PHICE /DTG
!/IC4 TAUICE(:)=TAUICE(:)/DTG
ELSE ! INFLAGS2(4).AND.ICE.GT.0
!/IS2 IF (IS2PARS(10).LT.0.5) THEN
!/IS2 ICEF = 0.
!/IS2 ENDIF
END IF
END IF ! INFLAGS2(4).AND.ICE.GT.0
!
!
! - - - - - - - - - - - - - - - - - - - - - -
Expand Down