w3sic2md Module Reference

Calculate ice dissipation source term S_{ice}. More...

Functions/Subroutines
subroutine, public	w3sic2 (A, DEPTH, ICEH, ICEF, CG, WN, IX, IY, S, D, WN_R, CG_ICE, ALPHA, R)
	S_{ice} source term using 5 parameters read from input files. More...

Detailed Description

Calculate ice dissipation source term S_{ice}.

Author

E. Rogers

S. Zieger

F. Ardhuin

G. Boutin

Date

05-Jan-2018

Copyright

Copyright 2009-2022 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.

Function/Subroutine Documentation

w3sic2()

subroutine, public w3sic2md::w3sic2	(	real, dimension(nspec), intent(in)	A,
		real, intent(in)	DEPTH,
		real, intent(in)	ICEH,
		real, intent(in)	ICEF,
		real, dimension(nk), intent(in)	CG,
		real, dimension(nk), intent(in)	WN,
		integer, intent(in)	IX,
		integer, intent(in)	IY,
		real, dimension(nspec), intent(out)	S,
		real, dimension(nspec), intent(out)	D,
		real, dimension(nk), intent(in)	WN_R,
		real, dimension(nk), intent(in)	CG_ICE,
		real, dimension(nk), intent(in)	ALPHA,
		real, dimension(nk), intent(in)	R
	)

S_{ice} source term using 5 parameters read from input files.

Parameters

[in]	A	Action density spectrum (1-D).
[in]	DEPTH	Local water depth.
[in]	ICEH	Ice thickness.
[in]	ICEF	Ice Floe diameter.
[in]	CG	Group velocities.
[in]	WN	Wavenumbers.
[in]	IX	Grid index.
[in]	IY	Grid index.
[out]	S	Source term (1-D version).
[out]	D	Diagonal term of derivative (1-D version).
[in]	WN_R	Wavenumbers in ice.
[in]	CG_ICE	Group velocities in ice.
[in]	ALPHA	Exponential decay rate of energy.
[in]	R	Ratio of energy to wave energy without ice.

Author

E. Rogers

S. Zieger

F. Ardhuin

G. Boutin

Date

04-Jan-2018

Definition at line 122 of file w3sic2md.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | 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,icei
#ifdef W3_T
    USE w3odatmd, ONLY: ndst
#endif
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
#ifdef W3_T0
    USE w3arrymd, ONLY: prt2ds
#endif
#ifdef W3_T1
    USE w3arrymd, ONLY: outmat
#endif
    !
    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
    !/
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
#ifdef W3_T0
    REAL                    :: DOUT(NK,NTH)
#endif
    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 ?
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3SIC2')
#endif
    !
    ! 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)
      !
#ifdef W3_T38
      WRITE (ndst,9000) depth,icecoef1,icecoef2
#endif
      !
      ! 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:
            !
#ifdef W3_IS2
            IF (ic2pars(8).GT.0) THEN
              wlg_r(ik)=tpi/wn_r(ik)
              smooth_dmax(ik)= (0.5*(1+tanh((icef-ic2pars(8)*wlg_r(ik))/(icef*0.5))))**2
            END IF
#endif
            !
            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
    !
#ifdef W3_T0
    DO ik=1, nk
      DO ith=1, nth
        dout(ik,ith) = d(ith+(ik-1)*nth)
      END DO
    END DO
    CALL prt2ds (ndst, nk, nk, nth, dout, sig(1:), '  ', 1.,    &
         0.0, 0.001, 'Diag Sice', ' ', 'NONAME')
#endif
    !
#ifdef W3_T1
    CALL outmat (ndst, d, nth, nth, nk, 'diag Sice')
#endif
    !
    ! Formats
    !
1001 FORMAT (/' *** WAVEWATCH III ERROR IN W3SIC2 : '/               &
         '     ',a,' REQUIRED BUT NOT SELECTED'/)
    !
#ifdef W3_T38
9000 FORMAT (' TEST W3SIC2 : DEPTH,ICECOEF1  : ',2e10.3)
#endif
    !/
    !/ End of W3SIC2 ----------------------------------------------------- /
    !/

References constants::abmin, w3gdatmd::dden, constants::delab, w3servmd::extcde(), w3gdatmd::flagll, constants::fwtable, constants::grav, w3gdatmd::gtype, w3gdatmd::ic2pars, w3idatmd::inflags2, w3gdatmd::mapwn, w3odatmd::ndse, w3odatmd::ndst, w3gdatmd::nk, w3gdatmd::nspec, w3gdatmd::nth, w3arrymd::outmat(), w3arrymd::prt2ds(), w3gdatmd::rlgtype, w3gdatmd::sig, constants::sizefwtable, w3servmd::strace(), constants::tpi, and w3gdatmd::ygrd.

Referenced by gxexpo(), and w3srcemd::w3srce().