#include "w3macros.h" !/ ------------------------------------------------------------------- / MODULE W3SIC4MD !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | C. Collins | !/ | E. Rogers | !/ | FORTRAN 90 | !/ | Last update : 21-Jan-2015 | !/ +-----------------------------------+ !/ !/ For updates see W3SIC4 documentation. !/ ! 1. Purpose : ! ! Calculate ice source term S_{ice} according to simple methods. ! Attenuation is a function of frequency and specified directly ! by the user. Example: a function is based on an exponential fit to ! the empirical data of Wadhams et al. (1988). ! ! 2. Variables and types : ! ! 3. Subroutines and functions : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! W3SIC4 Subr. Public ice source term. ! ---------------------------------------------------------------- ! ! 4. Subroutines and functions used : ! ! See subroutine documentation. ! ! 5. Remarks : ! ! Source material : ! 1) Wadhams et al. JGR 1988 ! 2) Meylan et al. GRL 2014 ! 3) Kohout & Meylan JGR 2008 in Horvat & Tziperman Cryo. 2015 ! 4) Kohout et al. Nature 2014 ! 5) Doble et al. GRL 2015 ! 6) Rogers et al. JGR 2016 ! Documentation of IC4: ! 1) Collins and Rogers, NRL Memorandum report 2017 ! ---> "A Source Term for Wave Attenuation by Sea ! Ice in WAVEWATCH III® : IC4" ! ---> describes original IC4 methods, 1 to 6 ! 2) Rogers et al., NRL Memorandum report 2018a ! ---> "Forecasting and hindcasting waves in and near the ! marginal ice zone: wave modeling and the ONR “Sea ! State” Field Experiment" ! ---> IC4 method 7 added ! 2) Rogers et al., NRL Memorandum report 2018b ! ---> "Frequency Distribution of Dissipation of Energy of ! Ocean Waves by Sea Ice Using Data from Wave Array 3 of ! the ONR “Sea State” Field Experiment" ! ---> New recommendations for IC4 Method 2 (polynomial fit) ! and IC4 Method 6 (step function via namelist) ! ! 6. Switches : ! ! See subroutine documentation. ! ! 7. Source code : !/ !/ ------------------------------------------------------------------- / !/ PUBLIC :: W3SIC4 !/ CONTAINS !/ ------------------------------------------------------------------- / SUBROUTINE W3SIC4 (A, DEPTH, CG, IX, IY, S, D) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | C. Collins | !/ | E. Rogers | !/ | FORTRAN 90 | !/ | Last update : 24-Feb-2017 | !/ +-----------------------------------+ !/ !/ 03-Dec-2015 : Origination ( version 5.09 ) !/ (starting from IC1) (C. Collins) !/ 03-Dec-2015 : W3SIC4 created, Methods 1,2,3,4 (C. Collins) !/ 21-Jan-2016 : IC4 added to NCEP repository (E. Rogers) !/ 27-Jan-2016 : Method 5 added (step function) (E. Rogers) !/ 08-Apr-2016 : Method 6 added (namelist step funct.) (E. Rogers) !/ 24-Feb-2017 : Corrections to Methods 1,2,3,4 (E. Rogers) !/ 13-Apr-2017 : Method 7 added (Doble et al. 2015) (E. Rogers) !/ !/ FIXME : Move field input to W3SRCE and provide !/ (S.Zieger) input parameter to W3SIC1 to make the subroutine !/ : versatile for point output processors ww3_outp !/ and ww3_ounp. !/ !/ Copyright 2009 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights !/ reserved. WAVEWATCH III is a trademark of the NWS. !/ No unauthorized use without permission. !/ ! 1. Purpose : ! ! S_{ice} source term using 5 parameters read from input files. ! These parameters are allowed to vary in space and time. ! The parameters control the exponential decay rate k_i ! Since there are 5 parameters, this permits description of ! dependence of k_i on frequency or wavenumber. ! !/ ------------------------------------------------------------------- / ! ! 2. Method : ! ! Apply parametric/empirical functions, e.g. from the literature. ! 1) Exponential fit to Wadhams et al. 1988, Table 2 ! 2) Polynomial fit, Eq. 3 from Meylan et al. 2014 ! 3) Quadratic fit to Kohout & Meylan'08 in Horvat & Tziperman'15 ! Here, note that their eqn is given as ln(alpha)=blah, so we ! have alpha=exp(blah) ! 4) Eq. 1 from Kohout et al. 2014 ! ! 5) Simple step function for ki as a function of frequency ! with up to 4 "steps". Controlling parameters KIx and FCx are ! read in as input fields, so they may be nonstationary and ! non-uniform in the same manner that ice concentration and ! water levels may be nonstationary and non-uniform. ! 444444444444 ! 33333333333 ! 222222222222 ! 1111111111111 ! ^ ^ ^ ! | | | ! 5 6 7 ! Here, 1 indicates ki=KI1=ICECOEF1 (ICEP1) ! 2 indicates ki=KI2=ICECOEF2 (ICEP2) ! 3 indicates ki=KI3=ICECOEF3 (ICEP3) ! 4 indicates ki=KI4=ICECOEF4 (ICEP4) ! 5 indicates freq cutoff #1 =FC5=ICECOEF5 (ICEP5) ! 6 indicates freq cutoff #2 =FC6=ICECOEF6 (MUDD) ! 7 indicates freq cutoff #3 =FC7=ICECOEF7 (MUDT) ! freq cutoff is given in Hz, freq=1/T (not sigma) ! Examples using hindcast, inversion with uniform ki: ! 5.1) Beaufort Sea, AWAC mooring, 2012, Aug 17 to 20 ! 0.0418 Hz to 0.15 Hz : ki=10e-6 ! 0.15 Hz to 0.175 Hz : ki=11e-6 ! 0.175 Hz to 0.25 Hz : ki=15e-6 ! 0.25 Hz to 0.5 Hz : ki=25e-6 ! 5.2) Beaufort Sea, AWAC mooring, 2012, Oct 27 to 30 ! 0.0418 Hz to 0.1 Hz : ki=5e-6 ! 0.1 Hz to 0.12 Hz : ki=7e-6 ! 0.12 Hz to 0.16 Hz : ki=15e-6 ! 0.16 Hz to 0.5 Hz : ki=100e-6 ! ICEP1=KI1=5.0e-6 ! ICEP2=KI2=7.0e-6 ! ICEP3=KI3=15.0e-6 ! ICEP4=KI4=100.0e-6 ! ICEP5=FC5=0.10 ! MUDD=FC6=0.12 ! MUDT=FC7=0.16 ! In terms of the 3-character IDs for "Homogeneous field ! data" in ww3_shel.inp, these are, respectively, IC1, IC2, ! IC3, IC4, IC5, MDN, MTH, and so this might look like: ! 'IC1' 19680606 000000 5.0e-6 ! 'IC2' 19680606 000000 7.0e-6 ! 'IC3' 19680606 000000 15.0e-6 ! 'IC4' 19680606 000000 100.0e-6 ! 'IC5' 19680606 000000 0.10 ! 'MDN' 19680606 000000 0.12 ! 'MTH' 19680606 000000 0.16 ! ! 6) Simple step function for ki as a function of frequency ! with up to 10 "steps". Controlling parameters KIx and FCx are ! read in as namelist parameters, so they are stationary and ! uniform. ! The last non-zero FCx value should be a large number, e.g. 99 Hz ! ! 4444444444 <--- ki=ic4_ki(4) ! 3333333333 <--- ki=ic4_ki(3) ! 2222222222 <--- ki=ic4_ki(2) ! 11111111111 <--- ki=ic4_ki(1) ! ^ ^ ^ ^ ! | | | | ! ic4_fc(1) ic4_fc(2) ic4_fc(3) ic4_fc(4)=large number ! Example: Beaufort Sea, AWAC mooring, 2012, Oct 27 to 30 ! &SIC4 IC4METHOD = 6, ! IC4KI = 0.50E-05, 0.70E-05, 0.15E-04, ! 0.10E+00, 0.00E+00, 0.00E+00, ! 0.00E+00, 0.00E+00, 0.00E+00, ! 0.00E+00, ! IC4FC = 0.100, 0.120, 0.160, ! 99.00, 0.000, 0.000, ! 0.000, 0.000, 0.000, ! 0.000 ! / ! ! 7) Doble et al. (GRL 2015), eq. 3. This is a function of ice ! thickness and wave period. ! ALPHA = 0.2*(T^(-2.13)*HICE or ! ALPHA = 0.2*(FREQ^2.13)*HICE ! ! More verbose description of implementation of Sice in WW3: ! See documentation for IC1 ! ! Notes regarding numerics: ! See documentation for IC1 ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! A R.A. I Action density spectrum (1-D) ! DEPTH Real I Local water depth ! CG R.A. I Group velocities. ! IX,IY I.S. I Grid indices. ! S R.A. O Source term (1-D version). ! D R.A. O Diagonal term of derivative (1-D version). ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing (!/S switch). ! PRT2DS Subr. W3ARRYMD Print plot output (!/T1 switch). ! OUTMAT Subr. W3ARRYMD Matrix output (!/T2 switch). ! ---------------------------------------------------------------- ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! W3SRCE Subr. W3SRCEMD Source term integration. ! W3EXPO Subr. N/A ASCII Point output post-processor. ! W3EXNC Subr. N/A NetCDF Point output post-processor. ! GXEXPO Subr. N/A GrADS point output post-processor. ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! None. ! ! 7. Remarks : ! ! If ice parameter 1 is zero, no calculations are made. ! For questions, comments and/or corrections, please refer to: ! Method 1 : C. Collins ! Method 2 : C. Collins ! Method 3 : C. Collins ! Method 4 : C. Collins ! Method 5 : E. Rogers ! Method 6 : E. Rogers ! Method 7 : E. Rogers ! ! ALPHA = 2 * WN_I ! Though it may seem redundant/unnecessary to have *both* in the ! code, we do it this way to make the code easier to read and ! relate to other codes and source material, and hopefully avoid ! mistakes. !/ ------------------------------------------------------------------- / ! ! 8. Structure : ! ! See source code. ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! !/T Enable general test output. ! !/T0 2-D print plot of source term. ! !/T1 Print arrays. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE CONSTANTS, ONLY: TPI USE W3ODATMD, ONLY: NDSE USE W3SERVMD, ONLY: EXTCDE USE W3GDATMD, ONLY: NK, NTH, NSPEC, SIG, MAPWN, IC4PARS, DDEN, & IC4_KI, IC4_FC, NIC4 USE W3IDATMD, ONLY: ICEP1, ICEP2, ICEP3, ICEP4, ICEP5, & MUDT, MUDV, MUDD, INFLAGS2 !/T USE W3ODATMD, ONLY: NDST !/S USE W3SERVMD, ONLY: STRACE !/T0 USE W3ARRYMD, ONLY: PRT2DS !/T1 USE W3ARRYMD, ONLY: OUTMAT ! IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list REAL, INTENT(IN) :: CG(NK), A(NSPEC), DEPTH REAL, INTENT(OUT) :: S(NSPEC), D(NSPEC) INTEGER, INTENT(IN) :: IX, IY !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ !/S INTEGER, SAVE :: IENT = 0 !/T0 INTEGER :: ITH !/T0 REAL :: DOUT(NK,NTH) INTEGER :: IKTH, IK, ITH, IC4METHOD, IFC REAL :: D1D(NK), EB(NK) REAL :: ICECOEF1, ICECOEF2, ICECOEF3, & ICECOEF4, ICECOEF5, ICECOEF6, & ICECOEF7, ICECOEF8 REAL :: KI1,KI2,KI3,KI4,FC5,FC6,FC7,FREQ REAL :: HS, EMEAN, HICE REAL, ALLOCATABLE :: WN_I(:) ! exponential decay rate for amplitude REAL, ALLOCATABLE :: ALPHA(:) ! exponential decay rate for energy REAL, ALLOCATABLE :: MARG1(:), MARG2(:) ! Arguments for M2 REAL, ALLOCATABLE :: KARG1(:), KARG2(:), KARG3(:) !Arguments for M3 !/ !/ ------------------------------------------------------------------- / !/ !/S CALL STRACE (IENT, 'W3SIC4') ! ! 0. Initializations ------------------------------------------------ * ! D = 0.0 ! ALLOCATE(WN_I(0:NK+1)) ALLOCATE(ALPHA(0:NK+1)) ALLOCATE(MARG1(0:NK+1)) ALLOCATE(MARG2(0:NK+1)) ALLOCATE(KARG1(0:NK+1)) ALLOCATE(KARG2(0:NK+1)) ALLOCATE(KARG3(0:NK+1)) MARG1 = 0.0 MARG2 = 0.0 KARG1 = 0.0 KARG2 = 0.0 KARG3 = 0.0 WN_I = 0.0 ALPHA = 0.0 ICECOEF1 = 0.0 ICECOEF2 = 0.0 ICECOEF3 = 0.0 ICECOEF4 = 0.0 ICECOEF5 = 0.0 ICECOEF6 = 0.0 ICECOEF7 = 0.0 ICECOEF8 = 0.0 HS = 0.0 HICE = 0.0 EMEAN = 0.0 ! ! IF (.NOT.INFLAGS2(-7))THEN ! WRITE (NDSE,1001) 'ICE PARAMETER 1' ! CALL EXTCDE(201) ! ENDIF ! ! We cannot remove the other use of INFLAGS below, ! because we would get 'array not allocated' error for the methods ! that don't use MUDV, etc. and don't have MUDV allocated. IF (INFLAGS2(-7)) ICECOEF1 = ICEP1(IX,IY) ! a.k.a. IC1 IF (INFLAGS2(-6)) ICECOEF2 = ICEP2(IX,IY) ! etc. IF (INFLAGS2(-5)) ICECOEF3 = ICEP3(IX,IY) IF (INFLAGS2(-4)) ICECOEF4 = ICEP4(IX,IY) IF (INFLAGS2(-3)) ICECOEF5 = ICEP5(IX,IY) ! Borrow from Smud (error if BT8 or BT9) !/BT8 WRITE (NDSE,*) 'DUPLICATE USE OF MUD PARAMETERS' !/BT8 CALL EXTCDE(202) !/BT9 WRITE (NDSE,*) 'DUPLICATE USE OF MUD PARAMETERS' !/BT9 CALL EXTCDE(202) IF (INFLAGS2(-2)) ICECOEF6 = MUDD(IX,IY) ! a.k.a. MDN IF (INFLAGS2(-1)) ICECOEF7 = MUDT(IX,IY) ! a.k.a. MTH IF (INFLAGS2(0 )) ICECOEF8 = MUDV(IX,IY) ! a.k.a. MVS IC4METHOD = IC4PARS(1) ! ! x. No ice --------------------------------------------------------- / ! ! IF ( ICECOEF1==0. ) THEN ! D = 0. ! WRITE(*,*) '!!!No Ice!!!' ! ! x. Ice ------------------------------------------------------------ / ! ELSE ! ! x.x Set constant(s) and write test output -------------------------- / ! ! (none) ! !/T38 WRITE (NDST,9000) DEPTH,ICECOEF1,ICECOEF2,ICECOEF3,ICECOEF4 ! ! 1. Make calculations ---------------------------------------------- / ! ! 1.a Calculate WN_I SELECT CASE (IC4METHOD) CASE (1) ! IC4M1 : Exponential fit to Wadhams et al. 1988 ALPHA = EXP(-ICECOEF1 * TPI / SIG - ICECOEF2) WN_I = 0.5 * ALPHA CASE (2) ! IC4M2 : Polynomial fit, Eq. 3 from Meylan et al. 2014 !NB: Eq. 3 only includes T^2 and T^4 terms, ! which correspond to ICECOEF3, ICECOEF5, so in ! regtest: ICECOEF1=ICECOEF2=ICECOEF4=0 MARG1 = ICECOEF1 + ICECOEF2*(SIG/TPI) + ICECOEF3*(SIG/TPI)**2 MARG2 = ICECOEF4*(SIG/TPI)**3 + ICECOEF5*(SIG/TPI)**4 ALPHA = MARG1 + MARG2 WN_I = 0.5 * ALPHA CASE (3) ! IC4M3 : Quadratic fit to Kohout & Meylan'08 in Horvat & Tziperman'15 HICE=ICECOEF1 ! For this method, ICECOEF1=ice thickness KARG1 = -0.3203 + 2.058*HICE - 0.9375*(TPI/SIG) KARG2 = -0.4269*HICE**2 + 0.1566*HICE*(TPI/SIG) KARG3 = 0.0006 * (TPI/SIG)**2 ALPHA = EXP(KARG1 + KARG2 + KARG3) WN_I = 0.5 * ALPHA CASE (4) !Eq. 1 from Kohout et al. 2014 !Calculate HS DO IK=1, NK EB(IK) = 0. DO ITH=1, NTH EB(IK) = EB(IK) + A(ITH+(IK-1)*NTH) END DO END DO DO IK=1, NK EB(IK) = EB(IK) * DDEN(IK) / CG(IK) EMEAN = EMEAN + EB(IK) END DO HS = 4.*SQRT( MAX(0.,EMEAN) ) ! If Hs < 3 m then do Hs dependent calc, otherwise dH/dx is a constant IF (HS <= 3) THEN WN_I=ICECOEF1 ! from: DHDX=ICECOEF1*HS and WN_I=DHDX/HS ELSE IF (HS > 3) THEN WN_I=ICECOEF2/HS ! from: DHDX=ICECOEF2 and WN_I=DHDX/HS END IF CASE (5) ! Simple step function (time- and/or space-varying) ! rename variables for clarity KI1=ICECOEF1 KI2=ICECOEF2 KI3=ICECOEF3 KI4=ICECOEF4 FC5=ICECOEF5 FC6=ICECOEF6 FC7=ICECOEF7 IF((KI1.EQ.0.0).OR.(KI2.EQ.0.0).OR.(KI3.EQ.0.0).OR. & (KI4.EQ.0.0).OR.(FC5.EQ.0.0).OR.(FC6.EQ.0.0).OR. & (FC7.EQ.0.0))THEN WRITE (NDSE,1001)'ICE PARAMETERS' CALL EXTCDE(201) END IF DO IK=1, NK FREQ=SIG(IK)/TPI ! select ki IF(FREQ.LT.FC5)THEN WN_I(IK)=KI1 ELSEIF(FREQ.LT.FC6)THEN WN_I(IK)=KI2 ELSEIF(FREQ.LT.FC7)THEN WN_I(IK)=KI3 ELSE WN_I(IK)=KI4 ENDIF END DO CASE (6) ! Simple step function (from namelist) ! error checking: require at least 3 steps IF((IC4_KI(1).EQ.0.0).OR.(IC4_KI(2).EQ.0.0).OR. & (IC4_KI(3).EQ.0.0).OR.(IC4_FC(1).EQ.0.0).OR. & (IC4_FC(2).EQ.0.0) )THEN WRITE (NDSE,1001)'ICE PARAMETERS' CALL EXTCDE(201) END IF DO IK=1, NK FREQ=SIG(IK)/TPI ! select ki DO IFC=1,NIC4 IF(FREQ.LT.IC4_FC(IFC))THEN WN_I(IK)=IC4_KI(IFC) EXIT END IF END DO END DO CASE (7) ! Doble et al. (GRL 2015) HICE=ICECOEF1 ! For this method, ICECOEF1=ice thickness DO IK=1,NK FREQ=SIG(IK)/TPI ALPHA(IK) = 0.2*(FREQ**2.13)*HICE END DO WN_I= 0.5 * ALPHA CASE DEFAULT WN_I = ICECOEF1 !Default to IC1: Uniform in k END SELECT ! ! 1.b Calculate DID ! DO IK=1, NK ! SBT1 has: D1D(IK) = FACTOR * MAX(0., (CG(IK)*WN(IK)/SIG(IK)-0.5) ) ! recall that D=S/E=-2*Cg*k_i D1D(IK) = -2. * CG(IK) * WN_I(IK) END DO ! ! 1.c Fill diagional matrix ! DO IKTH=1, NSPEC D(IKTH) = D1D(MAPWN(IKTH)) END DO ! ! END IF ! S = D * A ! ! ... Test output of arrays ! !/T0 DO IK=1, NK !/T0 DO ITH=1, NTH !/T0 DOUT(IK,ITH) = D(ITH+(IK-1)*NTH) !/T0 END DO !/T0 END DO ! !/T0 CALL PRT2DS (NDST, NK, NK, NTH, DOUT, SIG(1:), ' ', 1., & !/T0 0.0, 0.001, 'Diag Sice', ' ', 'NONAME') ! !/T1 CALL OUTMAT (NDST, D, NTH, NTH, NK, 'diag Sice') ! ! Formats ! 1001 FORMAT (/' *** WAVEWATCH III ERROR IN W3SIC4 : '/ & ' ',A,' REQUIRED BUT NOT SELECTED'/) ! !/T 9000 FORMAT (' TEST W3SIC4 : DEPTH,ICECOEF1 : ',2E10.3) !/ !/ End of W3SIC4 --------------------------------------------------- / !/ END SUBROUTINE W3SIC4 !/ !/ End of module W3SIC4MD ------------------------------------------ / !/ END MODULE W3SIC4MD