w3sic3md Module Reference

Functions/Subroutines
subroutine, public	w3sic3 (A, DEPTH, CG, WN, IX, IY, S, D)
subroutine, public	w3ic3wncg_v1 (WN_R, WN_I, CG, ICE1, ICE2, ICE3, ICE4, DPT)
real function, dimension(:), allocatable	delta (X)
subroutine, public	w3ic3wncg_cheng (WN_R, WN_I, CG, ICE1, ICE2, ICE3, ICE4, DPT)
subroutine, public	ic3table_cheng (ICE2, ICE3, ICE4)
subroutine	wn_precalc_cheng (WN, SIGMA, SIGMAM1, WN_O, WNM1, WNM2, ES_MOD, NU, DICE, HICE, DEPTH, SWITCHID, IK, KU)
subroutine	smooth_k (WN_R, WN_I, SIGMA, N, SWITCHID)
complex(8) function	drfun_dble_cheng (WN, SIGMA, ES, NU, DICE, HICE, DEPTH, JUDGE)
complex(8) function	drfun_quad_cheng (WN, SIGMA, ES, NU, DICE, HICE, DEPTH)

Variables
integer, save	calledic3table = 0

Function/Subroutine Documentation

delta()

real function, dimension(:), allocatable w3sic3md::delta

(

real, dimension(:), intent(in)

X

)

Definition at line 1733 of file w3sic3md.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           E. Rogers               |
    !/                  |           S. Zieger               |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         22-Oct-2013 |
    !/                  +-----------------------------------+
    !/
    !/    06-May-2013 : Origination (port from Clarkson.f90)( version 4.12 )
    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)
    !/
    !  1. Purpose :
    !
    !     This function calculates bin withs for any discretized function.
    !     May be used for numerical integration and differentiation.
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !      Parameter list
    !     ----------------------------------------------------------------
    !      X     R.A.      I  Array type of REAL
    !      DX    R.A       O  Bin widths if X
    !     ----------------------------------------------------------------
    !
    !  4. Remarks :
    !
    !     This function does not get used in update by S. Cheng.
    !     It should be removed if/when "V1" routines are removed.
    !
    !  5. Called by :
    !     W3IC3WNCG_V1
    !
    !  6. Source code :
    !/
    IMPLICIT NONE
    !/
    REAL, INTENT(IN)  :: X(:)
    REAL, ALLOCATABLE :: DX(:)
    INTEGER           :: IX, NX
    !/
    nx = SIZE(x,1)
    ALLOCATE(dx(nx))
    dx = 0.
    !/
    DO ix = 1,nx
      IF (ix==1) THEN
        dx(ix) = (x(ix+1)-x(ix  ))
      ELSE IF (ix==nx) THEN
        dx(ix) = (x(ix  )-x(ix-1))
      ELSE
        dx(ix) = (x(ix+1)-x(ix-1)) / 2.
      END IF
    END DO
    !/

Referenced by w3ic3wncg_v1().

drfun_dble_cheng()

complex(8) function w3sic3md::drfun_dble_cheng	(	complex(8), intent(in)	WN,
		real(8), intent(in)	SIGMA,
		real(8), intent(in)	ES,
		real(8), intent(in)	NU,
		real(8), intent(in)	DICE,
		real(8), intent(in)	HICE,
		real(8), intent(in)	DEPTH,
		integer	JUDGE
	)

Definition at line 2814 of file w3sic3md.F90.

    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           E. Rogers               |
    !/                  |           S. Zieger               |
    !/                  |           X. Zhao                 |
    !/                  |           S. Cheng                |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         13-Jan-2016 |
    !/                  +-----------------------------------+
    !/
    !/    06-May-2013 : Origination (port from Clarkson.f90)( version 4.10 )
    !/                                                        (E. Rogers)
    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)
    !/    30-Oct-2013 : Clarkson.f90 update added             (S. Zieger)
    !/
    !  1. Purpose :
    !
    !     Return dispersion relation function value for root finding
    !
    !  2. Method :
    !
    !     function based on dispersion relation derived by Wang and Shen
    !     2010
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !      FUNC1 CMPLX DBL O  Result (COMPLEX(8))
    !      WN    CMPLX DBL I  Wavenumber (COMPLEX(8))
    !      W     REAL  DBL I  Wave angular frequency
    !      ES    REAL  DBL I  Effective shear modulus on ice
    !      NU    REAL  DBL I  Effective viscosity
    !      DICE  REAL  DBL I  Density of ice
    !      HICE  REAL  DBL I  Thickness of ice
    !      DEPTH REAL  DBL I  Water depth
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name    Type  Module    Description
    !     ----------------------------------------------------------------
    !       None.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name    Type  Module   Description
    !     ----------------------------------------------------------------
    !      F_ZHAO_CHENG   xxx    xxxx      xxxxx
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !   Updated authors: Cheng and Shen.
    !   This code is based on Fortran code provided by Hayley Shen (Clarkson
    !     University) to Erick Rogers (NRL) on Aug 25 2015
    !   Original authors: Zhao and Shen.
    !   This code is based on Fortran code provided by Hayley Shen (Clarkson
    !     University) to Erick Rogers (NRL) on April 19 2013.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !/
    !/ ------------------------------------------------------------------- /
    USE constants, ONLY: grav, dwat
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    COMPLEX(8), INTENT(IN)   :: WN
    REAL(8), INTENT(IN)      :: SIGMA, ES, NU, DICE, HICE, DEPTH
    COMPLEX(8)               :: FUNC1,AA(4,4)                   ! RESULT
    INTEGER                  :: JUDGE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    COMPLEX(8)   :: VE,ALPHA,N,M,L,SK,CK,SA,CA,TH,THH,TEMP,J1,J2
    !/
    !/ ------------------------------------------------------------------- /
    ve    = cmplx( nu, es/dice/sigma )
    alpha = sqrt( wn**2. - sigma/ve * cmplx(0.,1.d0) )
    n     = sigma + 2. * ve * wn**2. * cmplx(0.,1.d0)
 
    temp  = exp(wn*hice)
    sk    = (temp - 1.d0/temp)/2.d0
    ck    = (temp + 1.d0/temp)/2.d0
 
    temp  = exp(alpha*hice)
    sa    = (temp - 1.d0/temp)/2.d0
    ca    = (temp + 1.d0/temp)/2.d0
    !
    temp  = (wn*depth)
    IF ( real(temp).LT.18.d0 ) THEN
      temp = exp(temp)
      th = (temp - 1./temp)/(temp + 1./temp)
    ELSE
      th = 1.d0
    ENDIF
    !      JUDGE==3 is not used yet
    IF ((es>=1.e5.AND.judge/=2).or.judge==3) THEN
      l     = 2 * wn * alpha * sigma * ve
      m     = (dble(dwat)/dice - 1) * dble(grav) * wn             &
           - dble(dwat) / dice * sigma**2 / th
      aa(1,1) = 0.
      aa(1,2) = 2 * cmplx(0.,1.) * wn**2.
      aa(1,3) = alpha**2. + wn**2.
      aa(1,4) = 0.
      !
      aa(2,1) = n * sigma
      aa(2,2) = -wn * dble(grav)
      aa(2,3) = cmplx(0.,1.) * wn * dble(grav)
      aa(2,4) = l
      !
      aa(3,1) = -2. * cmplx(0.,1.) * wn**2. * sk
      aa(3,2) =  2. * cmplx(0.,1.) * wn**2. * ck
      aa(3,3) =  (alpha**2. + wn**2.) * ca
      aa(3,4) = -(alpha**2. + wn**2.) * sa
      !
      aa(4,1) =   n * sigma * ck - m * sk
      aa(4,2) = - n * sigma * sk + m * ck
      aa(4,3) = -cmplx(0.,1.) * m * ca - l * sa
      aa(4,4) =  cmplx(0.,1.) * m * sa + l * ca
      !
      func1   = bsdet(aa,4)
    ELSE
      j1    = dice/dble(dwat)*(wn**2.*dble(grav)**2.*sk*sa - (n**4. &
           + 16.* ve**4.*wn**6.*alpha**2.)*sk*sa - 8. &
           *wn**3.*alpha*ve**2.*n**2.*(ck*ca-1.))
      j2    = (4.*wn**3.*alpha*ve**2.*sk*ca+n**2.*sa*ck &
           -dble(grav)*wn*sk*sa)
      IF (judge==2)THEN
        func1 = (sigma**2. - th*wn*dble(grav)) - th*j1/(j2+1.e-20)
      ELSEIF (judge==1)THEN
        func1 = (sigma**2. - th*wn*dble(grav))*j2 - th*j1
      ENDIF
    ENDIF
    !/

References constants::dwat, and constants::grav.

Referenced by smooth_k().

drfun_quad_cheng()

complex(8) function w3sic3md::drfun_quad_cheng	(	complex(8), intent(in)	WN,
		real(8), intent(in)	SIGMA,
		real(8), intent(in)	ES,
		real(8), intent(in)	NU,
		real(8), intent(in)	DICE,
		real(8), intent(in)	HICE,
		real(8), intent(in)	DEPTH
	)

Definition at line 2969 of file w3sic3md.F90.

    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           E. Rogers               |
    !/                  |           S. Zieger               |
    !/                  |           X. Zhao                 |
    !/                  |           S. Cheng                |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         13-Jan-2016 |
    !/                  +-----------------------------------+
    !/
    !/    06-May-2013 : Origination (port from Clarkson.f90)( version 4.10 )
    !/                                                        (E. Rogers)
    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)
    !/    30-Oct-2013 : Clarkson.f90 update added             (S. Zieger)
    !/
    !  1. Purpose :
    !
    !     Return dispersion relation function value for root finding
    !
    !  2. Method :
    !
    !     Use quadruple precision for computation
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !      FUNC1 CMPLX DBL O  Result (COMPLEX(8))
    !      WN    CMPLX DBL I  Wavenumber (COMPLEX(8))
    !      W     REAL  DBL I  Wave angular frequency
    !      ES    REAL  DBL I  Effective shear modulus on ice
    !      NU    REAL  DBL I  Effective viscosity
    !      DICE  REAL  DBL I  Density of ice
    !      HICE  REAL  DBL I  Thickness of ice
    !      DEPTH REAL  DBL I  Water depth
    !
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name    Type  Module    Description
    !     ----------------------------------------------------------------
    !      None.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name    Type  Module   Description
    !     ----------------------------------------------------------------
    !      F_ZHAO_CHENG   xxx    xxxx      xxxxx
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !   Updated authors: Cheng and Shen.
    !   This code is based on Fortran code provided by Hayley Shen (Clarkson
    !     University) to Erick Rogers (NRL) on Aug 25 2015
    !   Original authors: Zhao and Shen.
    !   This code is based on Fortran code provided by Hayley Shen (Clarkson
    !     University) to Erick Rogers (NRL) on April 19 2013.
    !   ER: S. Cheng had "COMPLEX(16)" for the local parameters. This is not
    !     supported by my compiler (gfortran on linux machine), so I changed
    !     it to "COMPLEX(8)"
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !/
    !/ ------------------------------------------------------------------- /
    USE constants, ONLY: grav, dwat
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    COMPLEX(8), INTENT(IN)   :: WN
    REAL(8), INTENT(IN)      :: SIGMA, ES, NU, DICE, HICE, DEPTH
    COMPLEX(8)               :: FUNC1                   ! RESULT
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    COMPLEX(8)   :: VE,ALPHA,N,M,L,SK,CK,SA,CA,TH,THH,TEMP,J1,J2
    !/
    !/ ------------------------------------------------------------------- /
    ve    = cmplx( nu, es/dice/sigma )
    alpha = sqrt( wn**2. - sigma/ve * cmplx(0.,1.d0) )
    n     = sigma + 2. * ve * wn**2. * cmplx(0.,1.d0)
 
    temp  = exp(wn*hice)
    sk    = (temp - 1.d0/temp)/2.d0
    ck    = (temp + 1.d0/temp)/2.d0
 
    temp  = exp(alpha*hice)
    sa    = (temp - 1.d0/temp)/2.d0
    ca    = (temp + 1.d0/temp)/2.d0
    !
    temp  = (wn*depth)
    IF ( real(temp).LT.18.d0 ) THEN
      temp = exp(temp)
      th = (temp - 1./temp)/(temp + 1./temp)
    ELSE
      th = 1.d0
    ENDIF
    !
    j1    = dice/dble(dwat)*(wn**2.*dble(grav)**2.*sk*sa - (n**4. &
         + 16.* ve**4.*wn**6.*alpha**2.)*sk*sa - 8. &
         *wn**3.*alpha*ve**2.*n**2.*(ck*ca-1.))
    j2    = (4.*wn**3.*alpha*ve**2.*sk*ca+n**2.*sa*ck &
         -dble(grav)*wn*sk*sa)
    func1 = (sigma**2. - th*wn*dble(grav))*j2 - th*j1
    !/

References constants::dwat, and constants::grav.

Referenced by smooth_k().

ic3table_cheng()

subroutine, public w3sic3md::ic3table_cheng	(	real	ICE2,
		real	ICE3,
		real	ICE4
	)

Definition at line 1963 of file w3sic3md.F90.

    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           E. Rogers               |
    !/                  |           S. Zieger               |
    !/                  |           X. Zhao                 |
    !/                  |           S. Cheng                |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         13-Jan-2016 |
    !/                  +-----------------------------------+
    !/
    !  1. Purpose :
    !
    !     It's a Preprocess part to create a table of wave nubmer,
    !        attenuation and group velocity
    !     for all ice thickness in deep water situation for main
    !        computation.
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !  4. Subroutines used :
    !
    !      Name               Type  Module   Description
    !     ----------------------------------------------------------------
    !      CMPLX_ROOT_MULLER_CHENG  Func. W3SIC3MD Find root for complex
    !                                         wavenumbers for waves in ice.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name              Type   Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE             Subr. W3WAVEMD
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !   Original authors: Cheng and Shen.
    !   This code is based on Fortran code provided by Hayley Shen (Clarkson
    !     University) to Erick Rogers (NRL) on Aug 25 2015
    !
    !   **UNRESOLVED BUG** This routine should be called again if ice
    !     rheology (visc., elast.) changes (either time or space).
    !     It doesn't. We need to set CALLEDIC3TABLE=0 if either parameter is
    !     changed.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nk,sig
    USE w3adatmd, ONLY: ic3wn_r, ic3wn_i, ic3cg
    USE w3idatmd, ONLY: inflags2
    USE w3odatmd, ONLY: ndse
    USE w3servmd, ONLY: extcde
    USE constants, ONLY: grav
    !
    IMPLICIT NONE
    REAL       :: ICE1, ICE2, ICE3, ICE4, DPT
    INTEGER ::  I1, I2, JITK, ITKNUM, IK
 
    dpt    = 999.
 
    itknum = ceiling(ic3_maxitk/ic3_ditk)
    ic3wn_r(:,0) = sig**2/grav
    ic3wn_i(:,0) = 0
    DO jitk = 1,itknum
      ice1 = jitk*ic3_ditk  !HICE
      CALL ic3precalc_cheng(ic3wn_r(:,jitk), ic3wn_i(:,jitk), &
           ic3cg(:,jitk), ice1, ice2, ice3, ice4, dpt)
    ENDDO
 
    RETURN

References w3servmd::extcde(), constants::grav, w3adatmd::ic3cg, w3adatmd::ic3wn_i, w3adatmd::ic3wn_r, w3idatmd::inflags2, w3odatmd::ndse, w3gdatmd::nk, w3gdatmd::sig, smooth_k(), w3dispmd::wavnu1(), and wn_precalc_cheng().

Referenced by w3wavemd::w3wave().

smooth_k()

subroutine w3sic3md::smooth_k	(	real, dimension(n)	WN_R,
		real, dimension(n)	WN_I,
		real, dimension(n), intent(in)	SIGMA,
		integer	N,
		integer	SWITCHID
	)

Definition at line 2498 of file w3sic3md.F90.

    REAL, INTENT(IN)  :: SIGMA(N)
    REAL              :: WN_R(N), WN_I(N),DIFF(N),REMOVEID(N)
    INTEGER           :: N,I,J,SWITCHID
    !
    diff = 0
    ! remove kr in mode switch zone,
    ! if it is a local extremum or suddenly increasing
    DO j = 1,3 ! 3 times to guarantee wavenumber increases monotonically
      removeid = 0
      DO i = 2,n
        diff(i) = wn_r(i) - wn_r(i-1)
      ENDDO
      DO i = 3,n
        IF(diff(i)<=0.OR.diff(i)>3*diff(i-1).OR.switchid==i)THEN
          removeid(i) = 1
          removeid(i-1) = 1
        ENDIF
      ENDDO
      ! fill removed location with kr,ki by interpolation
      DO i = 2,n-1
        IF (removeid(i) ==1) THEN
          wn_r(i) = wn_r(i-1) + (wn_r(i+1)-wn_r(i-1))/    &
               (sigma(i+1)-sigma(i-1))*(sigma(i)-sigma(i-1))
          wn_i(i) = wn_i(i-1) + (wn_i(i+1)-wn_i(i-1))/    &
               (sigma(i+1)-sigma(i-1))*(sigma(i)-sigma(i-1))
        ENDIF
      ENDDO
    ENDDO
    ! mode switch upward at the last frequencies
    IF (diff(n)>3*diff(n-1))THEN
      i = n
      wn_r(i) = wn_r(i-1) + (wn_r(i-1)-wn_r(i-2))/   &
           (sigma(i-1)-sigma(i-2))*(sigma(i)-sigma(i-1))
      wn_i(i) = wn_i(i-1) + (wn_i(i-1)-wn_i(i-2))/   &
           (sigma(i-1)-sigma(i-2))*(sigma(i)-sigma(i-1))
    ENDIF
    RETURN

References drfun_dble_cheng(), drfun_quad_cheng(), w3servmd::extcde(), and w3odatmd::ndse.

Referenced by ic3table_cheng(), and w3ic3wncg_cheng().

w3ic3wncg_cheng()

subroutine, public w3sic3md::w3ic3wncg_cheng	(	real, dimension(:), intent(inout)	WN_R,
		real, dimension(:), intent(inout)	WN_I,
		real, dimension(:), intent(inout)	CG,
		real, intent(in)	ICE1,
		real, intent(in)	ICE2,
		real, intent(in)	ICE3,
		real, intent(in)	ICE4,
		real, intent(in)	DPT
	)

Definition at line 1801 of file w3sic3md.F90.

    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           E. Rogers               |
    !/                  |           S. Zieger               |
    !/                  |           X. Zhao                 |
    !/                  |           S. Cheng                |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         13-Jan-2016 |
    !/                  +-----------------------------------+
    !/
    !  1. Purpose :
    !
    !     Calculation of complex wavenumber for waves in ice. Outsourced
    !     from W3SIC3 to allow update on wavenumbers and  group
    !     velocities at each time step an ice parameter is updated.
    !
    !  2. Method :
    !
    !     <text here>
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       WN_R   R. A.  I/O  Wave number (real part)
    !       WN_I   R. A.  I/O  Wave number (imag. part=wave attenuation)
    !       CG     R. A.  I/O  Group velocity
    !       ICE1   REAL    I   Thickness of ice                 [in m]
    !       ICE2   REAL    I   Effective viscosity of ice       [in m2/s]
    !       ICE3   REAL    I   Density of ice                  [in kg/m3]
    !       ICE4   REAL    I   Effective shear modulus of ice   [in Pa]
    !       DPT    REAL    I   Water depth                      [in m]
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name    Type  Module    Description
    !     ----------------------------------------------------------------
    !     WAVNU1   Subr. W3DISPMD   Wavenumber for waves in open water.
    !     WN_CMPLX Func. W3SIC3MD   Complex wavenumber for waves in ice.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name    Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3SIC3  Subr. W3SIC3MD Ice source term.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !     On pre-calculation table, S. Cheng says:
    !     Instead of interpolation, WN_R of arbitrary HICE is approximated
    !     by IC3WN_R in the look up table. IC3WN_R is related to an ice
    !     thickness in the table, which is closest to HICE and less than
    !     HICE
    !
    !     Fix submitted by Sukun March 2017:
    !         replace :
    !            I = MIN(INT(HICE/IC3_DITK)+1,ITKNUM)
    !         with :
    !            I = MIN(NINT(HICE/IC3_DITK),ITKNUM)
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Source code :
    !/
    !/ ------------------------------------------------------------------- /
    !/
    USE w3gdatmd, ONLY: nk,sig, ic3pars
    USE w3adatmd, ONLY: ic3wn_r, ic3wn_i, ic3cg
    USE w3odatmd, ONLY: ndse
    USE w3servmd, ONLY: extcde
    USE w3dispmd, ONLY: wavnu1
    !/
    IMPLICIT NONE
    !/
    REAL, INTENT(IN)   :: ICE1, ICE2, ICE3, ICE4, DPT
    REAL, INTENT(INOUT):: WN_R(:),WN_I(:),CG(:)
    REAL, ALLOCATABLE  :: SIGMA(:)
    !
    INTEGER     :: I, I1, I2, IK, KL,KU, ITKNUM
    COMPLEX(8)  :: WNCOMPLEX, X0,X1,X2, WNR, WNL
    REAL(8)     :: DEPTH, HICE, NU, DICE, ES_MOD, RR, K_OCEAN, &
         CG_OCEAN
    REAL        :: IC3HILIM,IC3KILIM
 
    ic3hilim=ic3pars(10)
    ic3kilim=ic3pars(11)
 
    ALLOCATE(  sigma( SIZE(cg) ) )
    sigma = 0.
    IF (SIZE(wn_r,1).EQ.nk) THEN
      kl    = 1
      ku    = nk
      i1    = 1
      i2    = nk
      sigma = sig(1:nk)
    ELSE IF (SIZE(wn_r,1).EQ.nk+2) THEN
      kl    = 1
      ku    = nk+2
      i1    = 0
      i2    = nk+1
      sigma = sig(0:nk+1)
    ELSE
      WRITE(ndse,900)
      CALL extcde(3)
    END IF
    depth   =   dpt     ! water depth
    hice    =   ice1    ! ice thickness
 
    ! Optional: limit ice thickness
    hice=min(dble(ic3hilim),hice)
    nu      =   ice2    ! "effective viscosity" parameter
    dice    =   ice3    ! density of ice
    es_mod  =   ice4    ! effective shear modulus of ice
    !
    itknum  = ceiling(ic3_maxitk/ic3_ditk)
    i       = min(nint(hice/ic3_ditk),itknum)
 
    ! Find values in pre-calculated look-up table
    ! See Remarks section.
    wn_r    =   ic3wn_r(i1:i2, i)
    wn_i    =   ic3wn_i(i1:i2, i)
    cg      =   ic3cg(i1:i2, i)
    rr      =   0.01
    !/ --- CHECK If it's shallow water situation, then it needs recalculate
    !              kr,ki,cg
    DO ik   =   kl,ku
      IF (wn_r(ik)*depth>4.0)THEN ! exit do-loop
        EXIT  ! assume kr is proportional to frequency
      ELSE
        x1  = cmplx(wn_r(ik),wn_i(ik))
        x0  = x1*(1-rr)
        x2  = x1*(1+rr)
        wncomplex = cmplx_root_muller_cheng(x0,x1,x2,1, &
             dble(sigma(ik)),es_mod,nu,dice,hice,depth)
        wn_i(ik)  = real(aimag(wncomplex))  !  ki
        wn_r(ik)  = real(wncomplex)         !  kr
      ENDIF
    ENDDO
 
    CALL smooth_k(wn_r,wn_i,sigma,ku-kl+ 1,0)
    ! Optional : limit ki
    DO ik   =   kl,ku
      wn_i(ik)  = min(wn_i(ik),ic3kilim)
    ENDDO
 
!!!     --- Update group velocitiy ----
    CALL cginic3_cheng(cg,sigma,wn_r,ku-kl+1)
 
    DEALLOCATE(sigma)
 
900 FORMAT (/' *** WAVEWATCH III ERROR IN W3SIC3_W3IC3WNCG : '/&
         '     CANNOT DETERMINE BOUNDS OF WAVENUMBER ARRAY.')

References w3servmd::extcde(), w3adatmd::ic3cg, w3gdatmd::ic3pars, w3adatmd::ic3wn_i, w3adatmd::ic3wn_r, w3odatmd::ndse, w3gdatmd::nk, w3gdatmd::sig, smooth_k(), and w3dispmd::wavnu1().

Referenced by w3sic3(), and w3wavemd::w3wave().

w3ic3wncg_v1()

subroutine, public w3sic3md::w3ic3wncg_v1	(	real, dimension(:), intent(inout)	WN_R,
		real, dimension(:), intent(inout)	WN_I,
		real, dimension(:), intent(inout)	CG,
		real, intent(in)	ICE1,
		real, intent(in)	ICE2,
		real, intent(in)	ICE3,
		real, intent(in)	ICE4,
		real, intent(in)	DPT
	)

Definition at line 658 of file w3sic3md.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           E. Rogers               |
    !/                  |           S. Zieger               |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         25-Oct-2013 |
    !/                  +-----------------------------------+
    !/
    !/    06-May-2013 : Origination (port from Clarkson.f90)( version 4.10 )
    !/                                                        (E. Rogers)
    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)
    !/
    !  1. Purpose :
    !
    !     Calculation of complex wavenumber for waves in ice. Outsourced
    !     from W3SIC3 to allow update on wavenumbers and  group
    !     velocities at each time step an ice parameter is updated.
    !
    !  2. Method :
    !
    !     <text here>
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       WN_R   R. A.  I/O  Wave number (real part)
    !       WN_I   R. A.  I/O  Wave number (imag. part=wave attenuation)
    !       CG     R. A.  I/O  Group velocity
    !       ICE1   REAL    I   Thickness of ice                 [in m]
    !       ICE2   REAL    I   Effective viscosity of ice       [in m2/s]
    !       ICE3   REAL    I   Density of ice                  [in kg/m3]
    !       ICE4   REAL    I   Effective shear modulus of ice   [in Pa]
    !       DPT    REAL    I   Water depth                      [in m]
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name    Type  Module    Description
    !     ----------------------------------------------------------------
    !     WAVNU1   Subr. W3DISPMD   Wavenumber for waves in open water.
    !     WN_CMPLX Func. W3SIC3MD   Complex wavenumber for waves in ice.
    !     WN_CMPLX_HF Func. W3SIC3MD   Like WN_CMPLX, but for h-f
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name    Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3SIC3  Subr. W3SIC3MD Ice source term.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !     See FORMAT 900.
    !
    !  7. Remarks :
    !
    !     Optional: Cap WN_I at 2.0E-4, since in simple tests with "normal
    !     resolution" (not finer than 1 km), WW3 has trouble resolving the
    !     dissipation if k_i>2e-4. Also, very large values of dissipation
    !     (e.g. k_i=100e-4) with IC3 occurs more in the higher frequencies
    !     which makes WW3 slow down quite a bit. This is done via ICEKILIM
    !     in the namelists.
    !
    !     This function does not get used in update by S. Cheng.
    !     It should be removed if/when "V1" routines are removed.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Source code :
    !/
    !/ ------------------------------------------------------------------- /
    !/
    USE w3gdatmd, ONLY: nk, sig, ic3pars
    USE w3dispmd, ONLY: wavnu1
    USE w3odatmd, ONLY: ndse
    USE w3servmd, ONLY: extcde
    USE constants, ONLY: tpi
    !/
    IMPLICIT NONE
    !/
    REAL, INTENT(INOUT):: WN_R(:),WN_I(:),CG(:)
    REAL, INTENT(IN)   :: ICE1, ICE2, ICE3, ICE4, DPT
 
    INTEGER            :: IK, KL,KU
    REAL, ALLOCATABLE  :: SIGMA(:),CG_IC3(:)
    REAL               :: K_OCEAN, CG_OCEAN
    DOUBLE PRECISION   :: KH, K_NOICE, HWAT, HICE, NU, DICE, ES_MOD
    DOUBLE PRECISION,PARAMETER :: KHMAX = 18.0d0 ! 18=OK, 19=fails
    DOUBLE COMPLEX     :: WNCOMPLEX,WNCOMPLEX_OLD
    REAL               :: STENSEC
    REAL               :: IC3HILIM,IC3KILIM
    !
    ALLOCATE( cg_ic3( SIZE(cg) ) )
    ALLOCATE(  sigma( SIZE(cg) ) )
    cg_ic3  = 0.
    sigma   = 0.
    stensec=tpi/10.0 ! sigma for T=10 sec
 
    ic3hilim=ic3pars(10)
    ic3kilim=ic3pars(11)
    !
    !     --- Input to routine (part 1): set 6 double precision variables
    !        using single precision variables. -------------------------- /
    hwat    = dble(dpt)     ! water depth
    hice    = dble(ice1)    ! ice thickness
    nu      = dble(ice2)    ! "effective viscosity" parameter
    dice    = dble(ice3)    ! density of ice
    es_mod  = dble(ice4)    ! effective shear modulus of ice
 
    ! Optional: limit ice thickness
    hice=min(dble(ic3hilim),hice)
 
    IF (SIZE(wn_r,1).EQ.nk) THEN
      kl    = 1
      ku    = nk
      sigma = sig(1:nk)
    ELSE IF (SIZE(wn_r,1).EQ.nk+2) THEN
      kl    = 1
      ku    = nk+2
      sigma = sig(0:nk+1)
    ELSE
      WRITE(ndse,900)
      CALL extcde(3)
    END IF
    !
    wncomplex_old=cmplx(0.0d0,0.0d0)
 
    DO ik = kl,ku
      !     --- Input to routine (part 2): set 2 double precision variables
      !         using single precision variable. --------------------------- /
      CALL wavnu1(sigma(ik),dpt,k_ocean,cg_ocean)
      k_noice = dble(k_ocean)
      !
      !     --- Muller Method fails for deep water: workaround follows ----- /
      kh        = k_noice * hwat ! kh w/out ice
      IF (kh.GT.khmax) THEN
        hwat   = khmax / k_noice
      ENDIF
      !     --- Calculate complex wavenumber ------------------------------- /
 
      IF((ik.GT.kl).AND.(sigma(ik).GT.stensec))THEN
 
        wncomplex = wn_cmplx_hf(dble(sigma(ik)),k_noice,es_mod,nu, &
             dice,hice,hwat,dble(sigma(ik-1)),wncomplex_old)
        wncomplex_old=wncomplex
 
      ELSE
 
        wncomplex = wn_cmplx_v1(dble(sigma(ik)),k_noice,es_mod,nu, &
             dice,hice,hwat)
        wncomplex_old=wncomplex
 
      ENDIF
 
      !     --- Output from function is type of DOUBLE COMPLEX. Set
      !         precision of imaginary to single precision array element --- /
      wn_i(ik)  = real(aimag(wncomplex))
 
      ! Optional : limit ki
      wn_i(ik)  = min(wn_i(ik),ic3kilim) ! see Remarks above
      wn_r(ik)  = real(wncomplex)
 
    END DO
    !     --- Update group velocitiy ----
    cg_ic3 = delta(sigma) / delta(wn_r)
    cg     = cg_ic3
 
    DEALLOCATE(cg_ic3)
    !
900 FORMAT (/' *** WAVEWATCH III ERROR IN W3SIC3_W3IC3WNCG : '/&
         '     CANNOT DETERMINE BOUNDS OF WAVENUMBER ARRAY.')
    !
    !/

References delta(), constants::dwat, w3servmd::extcde(), constants::grav, w3gdatmd::ic3pars, w3odatmd::ndse, w3gdatmd::nk, w3gdatmd::sig, constants::tpi, and w3dispmd::wavnu1().

Referenced by w3sic3(), and w3wavemd::w3wave().

w3sic3()

subroutine, public w3sic3md::w3sic3	(	real, dimension(nspec), intent(in)	A,
		real, intent(in)	DEPTH,
		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
	)

Definition at line 97 of file w3sic3md.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           E. Rogers               |
    !/                  |           S. Zieger               |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         11-Oct-2013 |
    !/                  +-----------------------------------+
    !/
    !/    06-May-2013 : Origination (copied from SICE1)     ( version 4.10 )
    !/                                                        (E. Rogers)
    !/    09-Oct-2013 : Update to meet WW3 coding standard    (S. Zieger)
    !/
    !/        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 :
    !
    !     Calculate ice source term S_{ice} according to a viscoelastic sea
    !     ice model (Wang and Shen 2010).
    !
    !     Reference: Wang, R., and H. H. Shen (2010), Gravity waves
    !     propagating into an ice‐covered ocean: A viscoelastic model, J.
    !     Geophys. Res., 115, C06024, doi:10.1029/2009JC005591 .
    !
    !/ ------------------------------------------------------------------- /
    !
    !  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)
    !
    !     Following W3SBT1 as a guide, equation 1 of W3SBT1 says:
    !         S = D * E
    !     However, the code of W3SBT1 has
    !         S = D * A
    !     This leads me to believe that the calling routine is
    !         expecting "S/sigma" not "S"
    !     Thus we will use D = S/E = -2 * Cg * k_i
    !
    !     The calling routine is expecting "S/sigma" not "S"
    !        Thus we will use D = S/E = -2 * Cg * k_i
    !        (see also documentation of W3SIC1)
    !
    !     Notes regarding numerics:
    !
    !     Experiments with constant k_i values suggest that results may be
    !        dependent on resolution if insufficient resolution is used.
    !        For detailed information, see documentation of W3SIC1.
    !
    !     Note regarding applicability/validity:
    !
    !     The Wang and Shen model is intended as a generalized model for
    !     various types of ice cover. It is a "continuum" model for
    !     which the same model is used from the ice edge to the ice
    !     interior. Though the ice types are expected to be very different
    !     from the edge to the interior, this is accomodated by the relative
    !     importance of the "effective viscosity" and the "modulus of
    !     elasticity". At the ice edge, where one finds frazil ice, pancake
    !     ice, or ice floes much smaller than the wave length, the "viscous"
    !     component of the model is believed to be most appropriate. At the
    !     interior, where one finds a continuous ice sheet, the "elastic
    !     model" component of the generalized visco-elastic model is
    !     expected to be appropriate. In addition to the case of continuous
    !     ice, Wang and Shen argue that the elastic model is also applicable
    !     to ice floes when the floe sizes are large relative to the
    !     wavelength. So to summarize,
    !     * frazil ice, pancake ice, and floes smaller than wavelength :
    !          viscosity dominates
    !     * continuous ice, and floes larger than wavelength :
    !          elasticity dominates
    !     * intermediate conditions: neither dominates
    !     All this is accomodated in WW3 by using non-uniform specification
    !     of viscosity and elasticity.
    !
    !     In the case where a user wishes to utilize only the "viscous
    !     model" aspect of Wang and Shen, and use an alternative scheme for
    !     continous ice and large ice floes, we allow this through the use
    !     of a user-defined namelist parameter "IC3MAXTHK". Floe size is
    !     not (at time of writing) an output available from ice model CICE,
    !     so we use the ice thickness as a way to anticipate the floe size.
    !     When ice thickness exceeds IC3MAXTHK, WW3 will use another model
    !     in place of the Wang and Shen formulation :
    !
    !     S_ice by F.A., an estimation of dissipation by turbulence
    !     at the ice-water interface. It uses only namelists for input, and
    !     no space/time varying input (though of course ice concentration is
    !     space/time varying). Unlike Liu et al. (IC2), it does not use
    !     ice thickness and does not yield a new C|Cg|k (i.e. it is non-
    !     dispersive), but it has the very nice feature of not requiring
    !     an eddy viscosity, which is a major drawback of the Liu et al.
    !     model. That is why we use it here, vs. Liu et al.
    !     (S_ice by Liu et al. and S_ice by F.A. are the two options
    !     available in IC2, i.e. w3sic2md.ftn)
    !
    !     At time of writing (March 23 2015), there may be some problems
    !     with the root-selection of IC3. For example with these settings:
    !     hice=1   ; rho_ice=917.0 ; emod=4.9e+12 ; visc=5e+7 ;
    !     h_water=deep (without using the IC3MAXTHK feature) the solution
    !     is rather irregular:
    !     T=11.11    k_i = 215.0e-5
    !     T=10       k_i = 266.0e-5
    !     T=9        k_i =   1.4e-5
    !     T=8.1      k_i =   1.7e-5
    !     With hice=0.1, ki is monotonically increasing in that range:
    !     T=11.11    k_i = 0.97e-5
    !     T=10       k_i = 2.64e-5
    !     T=9        k_i = 3.90e-5
    !     T=8.1      k_i = 4.47e-5
    !     Of course, when using IC3MAXTHK=0.1, the first example (hice=1)
    !     would switch to the "S_ice by F.A." model, and so this problem
    !     is circumvented.
    !
    !  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.
    !       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).
    !     ----------------------------------------------------------------
    !
    !  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.
    !
    !     Code by S. Cheng sets NOICE=.TRUE. if ISNAN(ICECOEF1).
    !     Comments are "maps may not be compatible"
    !     This feature is not understood by me (ER) and so omitted.
    !
    !/ ------------------------------------------------------------------- /
    !     On array size, S. Cheng says:
    !     Upon checking the origin of CG, I would say CG = CG_IC3, this is
    !     the topic ‘call W3IC3WNCG twice’. Recently I find they are not
    !     exactly the same due to calculation and smoothing of cg using
    !     several neighbor points. For different input array size, results
    !     are slightly different. In subr. w3wave, the size of input arrays
    !     is 0:NK+1, while in w3sic3, the size of all input arrays is 1:NK.
    !     This array size  difference is reflected in the resulting cg. The
    !     small difference in cg between calling IC3 twice or just calling
    !     once produces small difference in SWH. To eliminate this small
    !     difference, I suggest to keep CG instead of CG_IC3, as well as WN,
    !     WN_I, because other source terms use CG. I confirmed this change
    !     would make results of ICE the same whether calling twice or once
    !     by defining dimension WN_R, WN_I, CG_IC3 as 0:NK+1 instead of
    !     1:NK. Then CG_IC3 = CG.
    !/ ------------------------------------------------------------------- /
    !     On optimization, S. Cheng says:
    !     For Wang and Shen’s model, D does not change in the loop
    !     corresponding to NSTEPS in subr. W3SRCE. I find the most efficient
    !     and easy way to speed up is that Add D and NSTEPS as inputs of
    !     W3SIC3
    !     If NSTEPS==1
    !          Current Wang and Shen’s model code above.
    !     ELSE
    !          S = D*A
    !     Endif
    !/ ------------------------------------------------------------------- /
    !
    !  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, dwat, abmin, delab, sizefwtable,      &
         fwtable, grav
    USE w3odatmd, ONLY: ndse, iaproc, naproc, naperr
    !     USE WMMDATMD, ONLY: IMPROC, NMPERR ! WMMDATMD unavailable to outp
    USE w3servmd, ONLY: extcde
    USE w3gdatmd, ONLY: nk, nth, nspec, sig, mapwn, ic3pars, dden,  &
         flagll, ygrd, gtype, rlgtype
    USE w3idatmd, ONLY: icep1, icep2, icep3, icep4, icep5, icei,    &
         inflags2
#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)        :: CG(NK),  WN(NK), A(NSPEC), DEPTH
    REAL, INTENT(OUT)       :: S(NSPEC), D(NSPEC)
    INTEGER, INTENT(IN)     :: IX, IY
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    INTEGER                    :: ITH
#ifdef W3_T0
    REAL                   :: DOUT(NK,NTH)
#endif
    INTEGER                    :: IKTH, IK
    REAL                       :: ICECOEF1, ICECOEF2, ICECOEF3, &
         ICECOEF4, ICECOEF5, ICECONC
    REAL, DIMENSION(NK)        :: D1D, WN_I, WN_R, CG_IC3, CG_TMP
    LOGICAL                    :: NOICE
    REAL                       :: VISCM=1.83e-6
    REAL                       :: FREQ
    !  ............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,  &
         MAXTHK, MAXCNC, USE_CHENG,    &
         USE_CGICE, FIXEDHICE,   &
         FIXEDVISC,FIXEDDENS,FIXEDELAS
    INTEGER                 :: IND, IS, NUMIN
    !
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3SIC3')
#endif
    !
    ! 0.  Initializations ------------------------------------------------ /
    !
    !
    d        = 0.0
    d1d      = 0.0
    !
    wn_r     = wn
    wn_i     = 0.0
    cg_ic3   = 0.0
    cg_tmp   = 0.0
    !
    icecoef1 = 0.0
    icecoef2 = 0.0
    icecoef3 = 0.0
    icecoef4 = 0.0
    icecoef5 = 0.0
    iceconc  = 0.0
    !
    !     Rename variables to make code easier to read.
    maxthk=ic3pars(1)
    maxcnc=ic3pars(8)
    use_cheng=ic3pars(9)
    use_cgice=ic3pars(12)
    fixedhice=ic3pars(13)
    fixedvisc=ic3pars(14)
    fixeddens=ic3pars(15)
    fixedelas=ic3pars(16)
 
    !     --- Error checking for input ----------------------------------- /
    !     --- Allow one and only one input option for each variable ------ /
    numin=0
    IF (inflags2(-7)) numin=numin+1
    IF (fixedhice.GE.0.0) numin=numin+1
    IF (numin.NE.1) THEN
      IF ( iaproc .EQ. naperr )                       &
           WRITE (ndse,1001) 'ICE PARAMETER 1 (HICE)',numin
      CALL extcde(2)
    ENDIF
 
    numin=0
    IF (inflags2(-6)) numin=numin+1
    IF (fixedvisc.GE.0.0) numin=numin+1
    IF (numin.NE.1) THEN
      IF ( iaproc .EQ. naperr )                       &
           WRITE (ndse,1001) 'ICE PARAMETER 2 (VISC)',numin
      CALL extcde(2)
    ENDIF
 
    numin=0
    IF (inflags2(-5)) numin=numin+1
    IF (fixeddens.GE.0.0) numin=numin+1
    IF (numin.NE.1) THEN
      IF ( iaproc .EQ. naperr )                       &
           WRITE (ndse,1001) 'ICE PARAMETER 3 (DENS)',numin
      CALL extcde(2)
    ENDIF
 
    numin=0
    IF (inflags2(-4)) numin=numin+1
    IF (fixedelas.GE.0.0) numin=numin+1
    IF (numin.NE.1) THEN
      IF ( iaproc .EQ. naperr )                       &
           WRITE (ndse,1001) 'ICE PARAMETER 4 (ELAS)',numin
      CALL extcde(2)
    ENDIF
 
    !     --- Set local value to be used subsequently (ICEPx variables
    !         are not used beyond this point). --------------------------- /
    IF (inflags2(-7)) THEN
      icecoef1 = icep1(ix,iy) ! ice thickness
    ELSE
      icecoef1 = fixedhice
    ENDIF
 
    IF (inflags2(-6)) THEN
      icecoef2 = icep2(ix,iy) ! effective viscosity of ice cover
    ELSE
      icecoef2 = fixedvisc
    ENDIF
 
    IF (inflags2(-5)) THEN
      icecoef3 = icep3(ix,iy) ! density of ice
    ELSE
      icecoef3 = fixeddens
    ENDIF
 
    IF (inflags2(-4)) THEN
      icecoef4 = icep4(ix,iy) ! effective shear modulus of ice
    ELSE
      icecoef4 = fixedelas
    ENDIF
 
    !     ICECOEF5 = ICEP5(IX,IY) ! ICEP5 is inactive in W3SIC3
 
    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
 
      d1d=0.0
      !
      ! 2.  Ice ------------------------------------------------------------ /
    ELSEIF ( use_cheng==1.0 .AND.  &
         ((icecoef1.LE.maxthk).OR.(iceconc.LE.maxcnc)) ) THEN
 
      ! 2.a Write test output ---------------------------------------------- /
#ifdef W3_T38
      WRITE (ndst,9000) depth,icecoef1,icecoef2,icecoef3,icecoef4
#endif
 
      ! 2.b Make calculations using Cheng routines ------------------------- /
 
      !     --- Input to routine (part 1): 6 ice parameters from single
      !         precision variables. ---------------------------------------
 
      CALL w3ic3wncg_cheng(wn_r, wn_i, cg_ic3, icecoef1, icecoef2, &
           icecoef3, icecoef4, depth)
      !
      !    --- calculate source term --------------------------------------- /
      !    --- see Remarks section re: array size -------------------------- /
      IF ( use_cgice==1.0 ) THEN
        cg_tmp=cg_ic3
      ELSE
        cg_tmp=cg
      ENDIF
      DO ik=1, nk
        !            recall that D=S/E=-2*Cg*k_i
        d1d(ik)= -2.0 * cg_tmp(ik) * wn_i(ik)
 
      END DO
 
    ELSEIF ( (icecoef1 .LE. maxthk) .OR. (iceconc .LE. maxcnc) ) THEN
      !.......... e.g. if ice thickness is .le. 10 cm
      !...............or concentration is .le. 1.0
      !
      ! 2.a Write test output ---------------------------------------------- /
#ifdef W3_T38
      WRITE (ndst,9000) depth,icecoef1,icecoef2,icecoef3,icecoef4
#endif
      !
      ! 2.b Make calculations using original routines ---------------------- /
      !     --- Input to routine (part 1): 6 ice parameters from single
      !         precision variables. ---------------------------------------
 
      CALL w3ic3wncg_v1(wn_r, wn_i, cg_ic3, icecoef1, icecoef2, &
           icecoef3, icecoef4, depth     )
      !
      !    --- calculate source term --------------------------------------- /
      IF ( use_cgice==1.0 ) THEN
        cg_tmp=cg_ic3
      ELSE
        cg_tmp=cg
      ENDIF
      DO ik=1, nk
        !            recall that D=S/E=-2*Cg*k_i
        d1d(ik)= -2.0 * cg_tmp(ik) * wn_i(ik)
      END DO
      !
    ELSE ! .. e.g. if ice thickness is .gt. 10 cm
      ! Alternative by F.A., see Remarks section.
      IF (ic3pars(2).GT.0.) THEN
        uorb=0.
        aorb=0.
        fturb = ic3pars(2)
        IF (ic3pars(7).GT.0) THEN
          IF (ygrd(iy,ix).LT.0.AND.gtype.EQ.rlgtype.AND.flagll) &
               fturb = ic3pars(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
          !
          uorb = uorb + eb *sig(ik)**2 * dden(ik) / cg(ik)
          aorb = aorb + eb             * dden(ik) / cg(ik)
          !deep water only
        END DO
        !
        aorb = 2*sqrt(aorb)  ! significant amplitude
        uorb = 2*sqrt(uorb)  ! significant amplitude
 
        re = uorb*aorb / viscm
        smooth = 0.5*tanh((re-ic3pars(4))/ic3pars(5))
        pturb=(0.5+smooth)
        pvisc=(0.5-smooth)
 
        xi=(alog10(max(aorb/ic3pars(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=-1.* fturb*fw*uorb/grav
      ELSE ! so case of IC3PARS(2).LE.0.
        dturb = 0.
      END IF ! IF (IC3PARS(2).GT.0.)
 
      DO ik=1, nk
        dvisc = -1. *ic3pars(6) * wn(ik) * sqrt(viscm* sig(ik) / 2.)
        d1d(ik) = pturb*dturb*sig(ik)**2 +  pvisc*dvisc
      END DO
 
    END IF !   IF ( NOICE ) THEN
 
    ! 2.c Fill diagional matrix ------------------------------------------ /
    !
    DO ikth=1, nspec
      d(ikth) = d1d(mapwn(ikth))
    END DO
 
    !
    ! sign convention (example):
    !        S is from -10e-3 to 0
    !        A is from 0 to 10
    ! See Remarks section re: optimization
    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 W3SIC3 : '/               &
         '     ',a,' REQUIRED ONCE, BUT WAS PROVIDED BY USER '/ &
         '     ',i4,' TIMES.'/)
    !
#ifdef W3_T
9000 FORMAT (' TEST W3SIC3 : depth and 4 ice coef. : ',5e10.3)
#endif
    !/
    !/ End of W3SIC3 ----------------------------------------------------- /
    !/

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

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

wn_precalc_cheng()

subroutine w3sic3md::wn_precalc_cheng	(	complex(8), intent(out)	WN,
		real(8)	SIGMA,
		real(8)	SIGMAM1,
		complex(8)	WN_O,
		complex(8)	WNM1,
		complex(8)	WNM2,
		real(8)	ES_MOD,
		real(8)	NU,
		real(8)	DICE,
		real(8)	HICE,
		real(8)	DEPTH,
		integer	SWITCHID,
		integer	IK,
		integer	KU
	)

Definition at line 2166 of file w3sic3md.F90.

    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           S. Cheng                |
    !/                  |           E. Rogers               |
    !/                  |           S. Zieger               |
    !/                  |           X. Zhao                 |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         13-Jan-2016 |
    !/                  +-----------------------------------+
    !/
    !  1. Purpose :
    !
    !     Calculate complex wavenumber for waves in ice.
    !
    !  2. Method :
    !
    !     Wang and Shen (JGR 2010)
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       WN     CMPLX  O  Wave number (imag. part=wave attenuation)
    !       SIGMA  REAL   I  Wave angular frequency           [in rad]
    !       WN_O   REAL   I  Wave number (open water)
    !       ES     REAL   I  Effective shear modulus of ice   [in Pa]
    !       NU     REAL   I  Effective viscosity of ice       [in m2/s]
    !       DICE   REAL   I  Density of ice                  [in kg/m3]
    !       HICE   REAL   I  Thickness of ice                 [in m]
    !       DEPTH  REAL   I  Water depth                      [in m]
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name               Type  Module   Description
    !     ----------------------------------------------------------------
    !      CMPLX_ROOT_MULLER_CHENG  Func. W3SIC3MD Find root for complex
    !                                         wavenumbers for waves in ice.
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      IC3PRECALC_CHENG    xxx  xxx  xxxx
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !   Updated authors: Cheng and Shen.
    !   This code is based on Fortran code provided by Hayley Shen (Clarkson
    !     University) to Erick Rogers (NRL) on Aug 25 2015
    !   Original authors: Zhao and Shen.
    !   This code is based on Fortran code provided by Hayley Shen (Clarkson
    !     University) to Erick Rogers (NRL) on April 19 2013.
    !
    !   Hayley Shen says,
    !     We have determined that it may not be necessary to use curve
    !     fitting or lookup tables to get the group velocity and the
    !     attenuation coefficient. Attached is a short report with some
    !     sample numerical solutions. To implement the viscoelastic model,
    !     there are 4 fortran programs. According to Xin Zhao, the graduate
    !     student, it is very fast to find roots. I suggest that perhaps you
    !     try the pure viscous case by setting G=0 to start with. nu can be
    !     set at 0.05*ice concentration (m^2/s) to begin with, because for
    !     grease ice Newyear’s data showed nu to be about 0.02-0.03 m^2/s.
    !     By setting G=0 in you get exactly the Keller model for pure
    !     viscous layer.
    !
    !   This routine provides the initial guess according to the parameters
    !   of the present case. T>10s use open water, T<10s cases, calculate
    !   T=10s first using open water as the initial guess.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants, ONLY: tpi
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    REAL(8)       :: SIGMA,SIGMAM1,ES_MOD,NU,DICE,HICE,DEPTH
    COMPLEX(8)    :: WN_O, WNM1, WNM2, WN0, WN1,WN2
    COMPLEX(8),INTENT(OUT) :: WN ! RESULT
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER       :: I,IX,NUM,SWITCHID,IK,KU
    REAL(8)       :: RR, R2, EPS, SIGMA0, KAPPA, DIS
    COMPLEX(8)    :: X0, X1, X2, kp, ks, Gv
    !/
    !/ ------------------------------------------------------------------- /
    !   compute shear wave and pressure wave modes
    gv = cmplx(es_mod,-sigma*nu*dice)
    kp = sigma/sqrt(4*gv/dice)
    ks = 2*kp
    !   RR,R2 are empirical coefficients
    rr  = 0.2
    r2  = 4
    eps = 1.d-10   ! assuming variable =0, if it is less than this value
    !   compute root 1
    !   initial guesses from wave number of open water
    !
    x1  = wn_o    ! initial guess
    x0  = x1*(1-rr)
    x2  = x1*(1+rr)
 
    wn1 = cmplx_root_muller_cheng(x0,x1,x2,1,sigma, &
         es_mod,nu,dice,hice,depth)
    !   if root finder failed, or found shear wave mode,
    !   redo searching using simplified dispersion relation form, index 2
    IF (real(wn1)<0.or.abs(wn1-ks)/abs(ks)<0.03 .or. &
         abs(wn1-kp)/abs(kp)<0.03.and.es_mod>1.e7*nu)THEN
      wn1 = cmplx_root_muller_cheng(x0,x1,x2,2,sigma, &
           es_mod,nu,dice,hice,depth)
    ENDIF
    !   similar as wn1, but search with opposite order of inital guesses
    wn2 = cmplx_root_muller_cheng(x2,x1,x0,1,sigma, &
         es_mod,nu,dice,hice,depth)
    IF (real(wn2)<0.or.abs(wn2-ks)/abs(ks)<0.03)THEN
      wn2 = cmplx_root_muller_cheng(x2,x1,x0,2,sigma, &
           es_mod,nu,dice,hice,depth)
    ENDIF
    IF(abs(real(wn1)-real(wn_o))<abs(real(wn2)-real(wn_o)))THEN
      wn0 = wn1
    ELSE
      wn0 = wn2
    ENDIF
 
    IF(sigmam1==0.)THEN
      wn = wn0
    ELSE
      !      compute root 2
      !      Calculate the other wave number based on last frequency
      !      in the frequency array
      r2    = max(abs(wnm1-wnm2)/abs(wnm2), abs((sigma-sigmam1)/sigma))
      DO i  = 1,7,2
        IF(i<3.AND.sigma>4)THEN
          cycle
        ENDIF
        num = 2**(i-1)
        rr  = min(max(r2/real(num,8),0.001d0),0.5d0)
        x1  = wnm1
        DO ix  = 1,num
          x0 = x1*(1-rr)
          x2 = x1*(1+rr)
          sigma0 = sigmam1 + (1.0*ix)/num*(sigma-sigmam1)
          wn = cmplx_root_muller_cheng(x0,x1,x2,1,sigma0, &
               es_mod,nu,dice,hice,depth)
          IF(real(wn)<0)THEN   ! try another searching direction
            wn = cmplx_root_muller_cheng(x2,x1,x0,1,sigma0, &
                 es_mod,nu,dice,hice,depth)
          ENDIF
          IF(real(wn)<0)THEN   ! try another dispersion relation form
            wn = cmplx_root_muller_cheng(x0,x1,x2,3,sigma0, &
                 es_mod,nu,dice,hice,depth)
          ENDIF
          kp = sigma0/sqrt(4*gv/dice)
          ks = 2*kp
          !               set 3 means simple dispersion relation form
          IF(abs(wn-ks)/abs(ks)<0.03.OR.real(wn)<0)THEN
            wn = cmplx_root_muller_cheng(x0,x1,x2,2,sigma0, &
                 es_mod,nu,dice,hice,depth)
          ENDIF
          x1 = wn
          IF( real(wn)<0.99*real(wnm1).OR. abs(wn-wn0)>eps .AND. &
               (abs(x1-wn)/abs(wn)>0.3 .OR.                  &
               imag(wn)/(imag(x1)+eps)>10.OR. real(wn)<0))THEN
            EXIT  ! redo with smaller intervals
          ENDIF
          x0 = x1
          x1 = wn
        ENDDO !            DO IX  = 1,NUM
 
        IF(ix==num+1)THEN
          EXIT
        ENDIF
        wn = x1  !/ --- if exit of inner loop: give an approximate wn
      ENDDO !         DO I  = 1,7,2
 
      ! part 2
      ! assume found two roots, choose from the two condidates.
      kp = sigma/sqrt(4*gv/dice)
      IF(real(wn0)>0.AND.imag(wn0)>=0.AND.abs(wn - wn0)>eps)THEN
        !do switch at last 3 points is not worth numerically.
        !For v>5.e-2, it is low chance to be wrong at the last 3 points
        ! we suppose to use two mode in the future
        IF(nu>0.AND.ik>=ku-1.OR.   &
             nu>0.AND.switchid/=0)THEN
          ! assume one mode switch for general viscoelastic model
        ELSE
          dis   = abs(real(wn)-real(wn_o))/abs(real(wn0)-real(wn_o))
          kappa = (imag(wn)+ eps)/(imag(wn0) + eps)
          ! wn0 has smaller attenuation and closer to k0
          IF ((dis >= 1 .AND. kappa>=1 .AND.  &
               imag(wn0)>=0.1*imag(wnm1).AND. &
               abs(wn-kp)<abs(wn0-kp)) .OR.  &
               ! wn0 has larger attenuation and closer to k0
               ( dis>=1 .AND. kappa<1 .AND. &
               ((kappa> 0.2 .AND. imag(wn0)/real(wn0)<0.5).or. &
               abs(real(wn)-real(kp))<abs(real(wn0)-real(kp)))).OR.&
               ! wn0 has smaller attenuation and farther to k0,
               ! wn0 could be wrong at high G
               ! (kappa>1 .and. dis<1 .and. dis>  0.8 ))then
               ( kappa>1 .AND. dis<1 .AND. &
               abs(real(wn)-real(wnm1))>   &
               abs(real(wn0)-real(wnm1)) ))THEN
            wn = wn0
            switchid = ik
            ! wn0 has lager attenuation and farther to k0
          ELSEIF(dis<1 .AND. kappa<1) THEN
            ! keep wn without change
          ENDIF !  IF ((DIS >= 1 .AND. KAPPA>=1 .AND.  &
        ENDIF !   IF(NU>0.AND.IK>=KU-2.OR.   &
 
        ! choose dominant mode is farther than pressure wave.
        !but it doens't work for high viscosity.
        if (abs(wn-kp)/abs(kp)<0.03) then
          wn = wn0
          switchid = ik
        endif
        !if (real(wn)<=real(wnm1))then
        !    wn = wn0
        !    switchid = ik
        !endif
 
        IF (real(wn0)>real(wnm1).AND.real(wn0)<real(wn).and. &
             abs(imag(wn0)-imag(wnm1))<abs(imag(wn)-imag(wnm1)))THEN
          ! mainly for pure elastic case has multiple branchs
          wn = wn0
          switchid = ik
        ENDIF
      ENDIF ! IF(REAL(WN0)>0.AND.IMAG(WN0)>=0.AND.ABS(WN - WN0)....
    ENDIF   !      IF(SIGMAM1==0.)THEN
 
    IF(real(wn)<0)THEN
      print*, "MULLER METHOD FAILED, ES_MOD,NU,HICE:",es_mod,nu,hice
    ENDIF
    RETURN
 

References constants::tpi.

Referenced by ic3table_cheng().

Variable Documentation

calledic3table

integer, save w3sic3md::calledic3table = 0

Definition at line 88 of file w3sic3md.F90.

  INTEGER,SAVE ::  CALLEDIC3TABLE = 0

Referenced by w3wavemd::w3wave().