#include "w3macros.h" !/ ------------------------------------------------------------------- / MODULE W3SIC2MD !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | E. Rogers | !/ | S. Zieger | !/ | F. Ardhuin & G. Boutin | !/ | FORTRAN 90 | !/ | Last update : 05-Jan-2018 | !/ +-----------------------------------+ !/ !/ 10-Mar-2014 : Generalization with turbulent BL ( version 5.01 ) !/ 05-Jan-2018 : Addition of floe size effect ( version 6.04 ) !/ !/ For updates see W3SIC1 documentation. !/ ! 1. Purpose : ! ! Calculate ice dissipation source term S_{ice}. ! Exponential decay rate according to Liu et al., which ! uses as input: 1) ice thickness, and 2) an eddy ! viscosity parameter. This method is non-uniform in ! frequency. This is discussed further below, in ! subroutine "LIU_REVERSE_DISPERSION". ! ! Includes generalization by F. Ardhuin with viscous and tubulent ! boundary layers. That part is activating by setting namelist ! parameters that define the under-ice roughness and a friction ! coefficient. For example: &IC2 IC2TURB = 1. , IC2ROUGH =0.0001 ! ! References for Subtype 2: ! Liu et al. 1991: JGR 96 (C3), 4605-4621 ! Liu and Mollo 1988: JPO 18 1720-1712 ! Stopa et al. 2016: The Cryosphere ! ! 2. Variables and types : ! ! 3. Subroutines and functions : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! W3SIC2 Subr. Public Ice source term. ! ---------------------------------------------------------------- ! ! 4. Subroutines and functions used : ! ! See subroutine documentation. ! ! 5. Remarks : ! ! Reference:Rogers, W.E. and M.D. Orzech, 2013: Implementation and ! Testing of Ice and Mud Source Functions in WAVEWATCH III(R), ! NRL/MR/7320--13-9462, 31pp. ! available from http://www7320.nrlssc.navy.mil/pubs.php ! Direct link: ! http://www7320.nrlssc.navy.mil/pubs/2013/rogers2-2013.pdf ! ! 6. Switches : ! ! See subroutine documentation. ! ! 7. Source code : !/ !/ ------------------------------------------------------------------- / !/ PUBLIC :: W3SIC2 !/ CONTAINS !/ ------------------------------------------------------------------- / SUBROUTINE W3SIC2 (A, DEPTH, ICEH, ICEF, CG, WN, IX, IY, S, D, WN_R, & CG_ICE, ALPHA, R) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | E. Rogers | !/ | S. Zieger | !/ | F. Ardhuin & G. Boutin | !/ | FORTRAN 90 | !/ | Last update : 04-Jan-2018 | !/ +-----------------------------------+ !/ !/ 16-Oct-2012 : Origination. ( version 4.04 ) !/ (E. Rogers) !/ 09-Oct-2013 : W3SIC1 SUBTYPE=2 outsourced to W3SIC2 (S. Zieger) !/ 10-Mar-2014 : Generalization with turbulent BL ( version 5.01 ) !/ 16-Feb-2016 : Passes ICEH as parameter ( version 5.10 ) !/ 02-May-2016 : Call to Liu disp moved to w3srce ( version 5.10 ) !/ 04-Jan-2018 : Includes floe size dependance ( version 6.02 ) !/ 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 : ! ! Regarding i/o (general to all Sice modules): S_{ice} source term ! is calculated using up to 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. ! ! Sea ice affects the wavenumber k of wind-generated ocean waves. ! The ice-modified wavenumber can be expressed as a complex number ! k = k_r + i*k_i, with the real part k_r representing impact of ! the sea ice on the physical wavelength and propagation speeds, ! producing something analogous to shoaling and refraction by ! bathymetry, whereas the imaginary part of the complex ! wavenumber, k_i, is an exponential decay coefficient ! k_i(x,y,t,sigma) (depending on location, time and frequency, ! respectively), representing wave attenuation, and can be ! introduced in a wave model such as WW3 as S_ice/E=-2*Cg*k_i, ! where S_ice is one of several dissipation mechanisms, along ! with whitecapping, for example, S_ds=S_wc+S_ice+⋯. The k_r - ! modified by ice would enter the model via the C calculations ! on the left-hand side of the governing equation.The fundamentals ! are straightforward, e.g. Rogers and Holland (2009 and ! subsequent unpublished work) modified a similar model, SWAN ! (Booij et al. 1999) to include the effects of a viscous mud ! layer using the same approach (k = k_r + i*k_i) previously. ! ! General approach is analogous to Rogers and Holland (2009) ! approach for mud. ! See text near their eq. 1 : ! k = k_r + i * k_i ! eta(x,t) = Real( a * exp( i * ( k * x - sigma * t ) ) ) ! a = a0 * exp( -k_i * x ) ! S / E = -2 * Cg * k_i (see also Komen et al. (1994, pg. 170) ! ! Please note that S is source term for action. ! ! Notes regarding numerics: ! (Note by F. Ardhuin: these may not apply in version 5 thanks to splitting ! of ice source terms and implicit integration in W3SRCE) ! Experiments with constant k_i values suggest that : ! for dx=20.0 km, k_i should not exceed 3.5e-6 ! (assumes 2.7% Hs error in my particular test case is intolerable) ! for dx=5.0 km, k_i should not exceed 2.0e-5 ! for dx=2.5 km, k_i should not exceed 5.0e-5 ! for dx=1.0 km, k_i should not exceed 2.0e-4 ! for dx=0.35 km, error is less than 2.1% for all k_i tested ! for dx=0.10 km, error is less than 1.3% for all k_i tested ! "Ground truth" used for this is an exponential decay profile. ! ! For reference, ACNFS is 1/12th deg, so delta_latitude=9.25 km. ! ! {put more equations here} ! ! The laminar to turbulent transition is described in ! Stopa et al. (The Cryosphere, 2016). ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! A R.A. I Action density spectrum (1-D) ! DEPTH Real I Local water depth ! ICEH Real I Ice thickness ! CG R.A. I Group velocities ! WN R.A. I Wavenumbers ! 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) ! WN_R R.A. I Wavenumbers in ice ! CG_ICE R.A. I Group velocities in ice ! ALPHA R.A. I Exponential decay rate of energy ! R R.A. I Ratio of energy to wave energy without ice ! ICEF Real I Ice Floe diameter ! ! imported via module: ! ICEP2 R.A. I Eddy viscosity ! ---------------------------------------------------------------- ! ! 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. ! ! 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 USE W3ODATMD, ONLY: NDSE USE W3SERVMD, ONLY: EXTCDE USE W3DISPMD USE W3GDATMD, ONLY: NK, NTH, NSPEC, SIG, MAPWN, IC2PARS, DDEN, & FLAGLL, YGRD, GTYPE, RLGTYPE USE W3IDATMD, ONLY: INFLAGS2,ICEP1,ICEP2,ICEP3,ICEP4,ICEP5,WLEV,& ICEI !/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) :: A(NSPEC), DEPTH, ICEH REAL, INTENT(IN) :: CG(NK), WN(NK) REAL, INTENT(OUT) :: S(NSPEC), D(NSPEC) REAL, INTENT(IN) :: ALPHA(NK) ! exponential (spatial) decay rate for energy (1/m) INTEGER, INTENT(IN) :: IX, IY REAL, INTENT(IN) :: WN_R(NK), CG_ICE(NK), R(NK) REAL, INTENT(IN) :: ICEF ! Hypothesis: friction does not occur for broken ice !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ !/S INTEGER, SAVE :: IENT = 0 !/T0 INTEGER :: ITH !/T0 REAL :: DOUT(NK,NTH) INTEGER :: IKTH, IK REAL :: D1D(NK) !In SBT1: D1D was named "CBETA" REAL :: ICECOEF1, ICECOEF2, ICECONC REAL, ALLOCATABLE :: WN_I(:) ! exponential (spatial) decay rate for amplitude (1/m) REAL :: VISCM=1.83E-6 ! molecular viscosity of water at freezing REAL :: PTURB, PVISC, DTURB, DVISC, & SMOOTH, RE, UORB, AORB, EB, & DELI1, DELI2, FW, XI, FTURB, & CG_EFF(NK), WLG_R(NK), SMOOTH_DMAX(NK) INTEGER :: IND, ITH, IS LOGICAL :: NOICE=.FALSE. ! Warning, ALPHA = 2 * WN_I -> Makes WN_I useless, doesnt it ? !/ !/ ------------------------------------------------------------------- / !/ !/S CALL STRACE (IENT, 'W3SIC2') ! ! 0. Initializations ------------------------------------------------ * ! D = 0.0 ! ALLOCATE(WN_I(NK)) WN_I = 0.0 ICECOEF1 = ICEH ICECOEF2 = 0.0 ICECONC = 0.0 CG_EFF = 0. SMOOTH_DMAX(:)=1. ! IF (INFLAGS2(-7))ICECOEF1 = ICEH IF (INFLAGS2(-6))ICECOEF2 = ICEP2(IX,IY) IF (INFLAGS2(4)) ICECONC = ICEI(IX,IY) ! ! ! 1. No ice --------------------------------------------------------- / ! NOICE=.FALSE. IF (ICECOEF1==0.0) NOICE=.TRUE. IF (INFLAGS2(4).AND.(ICECONC==0.0)) NOICE=.TRUE. IF ( NOICE ) THEN D = 0.0 ! ! 2. Ice ------------------------------------------------------------ / ELSE ! ! 2.a Set constant(s) and write test output -------------------------- / ! ! (none) ! !/T38 WRITE (NDST,9000) DEPTH,ICECOEF1,ICECOEF2 ! ! 2.b Make calculations ---------------------------------------------- / ! ICECOEF1 = H_ICE ! ICECOEF2 = VISC ! ! Branches out depending on choice of dispersion relation... ! by default IC2PARS(1)=0, and attenuation computed as described in Stopa et al. 2016 ! IF (IC2PARS(1).GT.0.5) THEN IF (.NOT.INFLAGS2(-7))THEN WRITE (NDSE,1001) 'ICE PARAMETER 1' CALL EXTCDE(2) ENDIF IF (.NOT.INFLAGS2(-6))THEN WRITE (NDSE,1001) 'ICE PARAMETER 2' CALL EXTCDE(2) ENDIF ! WN_I(:) = 0.5 * ALPHA(:) ! ALPHA=2*WN_I DO IK=1, NK ! recall that D=S/E=-2*Cg*k_i ! Note: We should not use CG_ICE here unless CG_ICE is also ! used for advection in w3wavemd.ftn (see lines for IC3 ! there). D1D(IK)= -2.0 * CG(IK) * WN_I(IK) END DO ! ! Alternative by F.A.: generalization to a turbulent boundary layer ! uses the ice-free dispersion, to be updated later ! ELSE ! goes here if IC2PARS(1).LE.0.5 (this is the default behavior) IF (IC2PARS(2).GT.0.) THEN UORB=0. AORB=0. FTURB = IC2PARS(2) ! Special treatment in the southern ocean ... IF (IC2PARS(7).GT.0) THEN IF (YGRD(IY,IX).LT.0.AND.GTYPE.EQ.RLGTYPE.AND.FLAGLL) FTURB = IC2PARS(7) END IF DO IK=1, NK EB = 0. DO ITH=1, NTH IS=ITH+(IK-1)*NTH EB = EB + A(IS) END DO ! ! UORB and AORB are the variances of the orbital velocity and surface elevation ! of the water relative to the ice ... this is only correct if the ice layer ! does not move. This should is changed by taking into account DMAX when IC2DMAX > 0: ! !/IS2 IF (IC2PARS(8).GT.0) THEN !/IS2 WLG_R(IK)=TPI/WN_R(IK) !/IS2 SMOOTH_DMAX(IK)= (0.5*(1+TANH((ICEF-IC2PARS(8)*WLG_R(IK))/(ICEF*0.5))))**2 !/IS2 END IF ! IF (R(IK).GT.1.) THEN UORB = UORB + EB * SMOOTH_DMAX(IK)* SIG(IK)**2 * DDEN(IK) / CG(IK) & / (R(IK)*CG_ICE(IK)/CG(IK)) AORB = AORB + EB * SMOOTH_DMAX(IK) * DDEN(IK) / CG(IK) & / (R(IK)*CG_ICE(IK)/CG(IK)) !deep water only ELSE UORB = UORB + EB * SMOOTH_DMAX(IK) *SIG(IK)**2 * DDEN(IK) / CG(IK) AORB = AORB + EB * SMOOTH_DMAX(IK) * DDEN(IK) / CG(IK) !deep water only END IF END DO ! AORB = 2*SQRT(AORB) ! significant amplitude UORB = 2*SQRT(UORB) ! significant amplitude RE = UORB*AORB / VISCM SMOOTH = 0.5*TANH((RE-IC2PARS(4))/IC2PARS(5)) PTURB=(0.5+SMOOTH) PVISC=(0.5-SMOOTH) XI=(ALOG10(MAX(AORB/IC2PARS(3),3.))-ABMIN)/DELAB IND = MIN (SIZEFWTABLE-1, INT(XI)) DELI1= MIN (1. ,XI-FLOAT(IND)) DELI2= 1. - DELI1 FW =FWTABLE(IND)*DELI2+FWTABLE(IND+1)*DELI1 DTURB= FTURB*FW*UORB/GRAV ELSE ! so case of IC2PARS(2).LE.0. DTURB = 0. PTURB = 0. PVISC = 1. END IF ! IF (IC2PARS(2).GT.0.) ! DO IK=1, NK ! WN_R is used here but warning, this is only OK for unbroken ice DVISC = IC2PARS(6) * WN_R(IK) * SQRT(VISCM* SIG(IK) / 2.) D1D(IK) = -1.*(PTURB*MAX(DTURB*SIG(IK)**2,DVISC) + PVISC*DVISC) & *SMOOTH_DMAX(IK) END DO END IF ! IF (IC2PARS(1).GT.0.5) ! ! 2.c Fill diagional matrix ! DO IKTH=1, NSPEC D(IKTH) = D1D(MAPWN(IKTH)) END DO ! END IF ! IF ( NOICE ) THEN ! 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 W3SIC2 : '/ & ' ',A,' REQUIRED BUT NOT SELECTED'/) ! !/T38 9000 FORMAT (' TEST W3SIC2 : DEPTH,ICECOEF1 : ',2E10.3) !/ !/ End of W3SIC2 ----------------------------------------------------- / !/ END SUBROUTINE W3SIC2 !/ !/ End of module W3SIC2MD -------------------------------------------- / !/ END MODULE W3SIC2MD