w3sic2md.F90

Go to the documentation of this file.

 
 
#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,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 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3sic2
 
  !/
  !/ End of module W3SIC2MD -------------------------------------------- /
  !/
END MODULE w3sic2md