#include "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE W3SIC5MD
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | Q. Liu |
!/ | E. Rogers |
!/ | FORTRAN 90 |
!/ | Last update : 02-Jun-2017 |
!/ +-----------------------------------+
!/
!/ 15-Mar-2016 : Origination. ( version 5.10 )
!/ ( Q. Liu )
!/ 15-Mar-2016 : Started from w3sic1/2/3/4 module ( Q. Liu )
!/
!/ 24-Apr-2017 : Adding more filters ( Q. Liu )
!/
!/ 29-Apr-2017 : Introducing CMPLX_TANH2 ( Q. Liu )
!/
!/ 02-Jun-2017 : Update to version 5.16 ( Q. Liu )
!/
!/ 17-Jun-2017 : Remove some unnecessary lines ( Q. Liu )
!/ (cg_ice, detla function, complx_tanh etc.)
!/
!/ 20-Aug-2018 : Ready to be merged to master (v6.06)( Q. Liu)
!/
!/ 1. Purpose :
! Calculate ice source term S_{ice} according to a viscoelastic sea
! ice model (Mosig et al. 2015)
!
! Reference:
! Mosig, J. E. M., F. Montiel, and V. A. Squire (2015),
! Comparison of viscoelastic-type models for ocean wave attenuation
! in ice-covered seas, J. Geophys. Res. Oceans, 120, 6072–6090,
! doi:10.1002/2015JC010881.
!
! 2. Variables and types :
!
! Name Type Scope Description
! ----------------------------------------------------------------
! KSP Int. Private the kind parameter for single precision
! real variables
! KDP Int. Private Same as KSP but for double precision
! KSPC Int. Private the kind parameter for single precision
! complex variables
! KDPC Int. Private Same as KSPC but for double precision
! ERRTOL Real Private A real parameter used for "==" test
! ----------------------------------------------------------------
!
! 3. Subroutines and functions :
!
! Name Type Scope Description
! ----------------------------------------------------------------
! W3SIC5 Subr. Public Ice source term
! W3IC5WNCG Subr. Public Wavenumber and group velocity of ice-
! coupled waves
! FSDISP Subr. Public Solving the ice-coupled wave dispersion
! BALANCING_MATRIX
! Subr. Private Balancing the matrix before we try to
! find its eigenvalues
! EIG_HQR Subr. Private QR algorithm for real Hessenberg matrix
! (eigenvalues-finding algorithm)
! POLYROOTS Subr. Private Finding roots of a general polynomial
! NR_CORR Func. Private Get the Newton-Raphson correction term
! for iteration
! NR_ROOT Func. Private Newton-Raphson algorithm for solving
! the ice-coupled wave dispersion
! CMPLX_SINH, CMPLX_COSH, CMPLX_TANH2
! Func. Private sinh, cosh, tanh for complex inputs
! INIT_RANDOM_SEED
! Subr. Private Initialize the random seed based on
! the system's time
! ----------------------------------------------------------------
!
! 4. Subroutines and functions used :
! See subroutine documentation
!
! 5. Remarks :
!
! 6. Switches :
! See subroutine documentation
!
! 7. Source code:
!/
!/ ------------------------------------------------------------------- /
IMPLICIT NONE
!/
PUBLIC :: W3SIC5, W3IC5WNCG, FSDISP
PRIVATE :: BALANCING_MATRIX, EIG_HQR, POLYROOTS
PRIVATE :: NR_CORR, NR_ROOT
PRIVATE :: CMPLX_SINH, CMPLX_COSH, CMPLX_TANH2
PRIVATE :: INIT_RANDOM_SEED
!/
PRIVATE :: KSP, KDP, KSPC, KDPC, ERRTOL
!/ ------------------------------------------------------------------- /
!/ Parameter list
! Kind for single- and double-precision real type
INTEGER, PARAMETER :: KSP = KIND(1.0)
INTEGER, PARAMETER :: KDP = KIND(1.0D0)
!
! Kind for single- and double-precision complex type
INTEGER, PARAMETER :: KSPC = KIND((1.0, 1.0))
INTEGER, PARAMETER :: KDPC = KIND((1.0D0, 1.0D0))
REAL, PARAMETER :: ERRTOL = 1.E-12
!/
CONTAINS
!/ ------------------------------------------------------------------- /
!/
SUBROUTINE W3SIC5 (A, DEPTH, CG, WN, IX, IY, S, D)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | Q. Liu |
!/ | E. Rogers |
!/ | FORTRAN 90 |
!/ | Last update : 25-Apr-2017 |
!/ +-----------------------------------+
!/
!/ 23-Mar-2016 : Origination ( version 5.10 )
!/ ( Q. Liu)
!/ 23-Mar-2016 : Started from w3sic1/2/3/4 subr. ( Q. Liu)
!/ 05-Apr-2016 : Options for Cg_{ice} or Cg ( Q. Liu)
!/ 25-Apr-2017 : Add more filters ( Q. Liu)
!/ 20-Aug-2018 : Ready to be merged to master (v6.06)( Q. Liu)
!/
!/ 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 (Mosig et al. 2015)
!
! Reference:
! Mosig, J. E. M., F. Montiel, and V. A. Squire (2015),
! Comparison of viscoelastic-type models for ocean wave attenuation
! in ice-covered seas, J. Geophys. Res. Oceans, 120, 6072–6090,
! doi:10.1002/2015JC010881.
!
! 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:
! --------------------------------------
! Similar to the Wang and Shen model used in w3sic3md, the FS model
! used here is also a viscoelastic-type model which treats the sea
! ice cover as a continuum and uses two empirical rheological para-
! meters, i.e., the effective shear modulus of ice G and the effec-
! tive viscosity η to characterize sea ices of various type. Please
! see the documentation of w3sic3md for a detailed discussion of
! this kind of model.
!
! 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 (!/T0 switch).
! OUTMAT Subr. W3ARRYMD Matrix output (!/T1 switch).
! W3IC5WNCG Subr. / Wavenumber and group velocity of ice-
! coupled waves
! ----------------------------------------------------------------
! * / means this module
!
! 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 :
!/ ------------------------------------------------------------------- /
!/
!/T USE W3ODATMD, ONLY: NDST
!/S USE W3SERVMD, ONLY: STRACE
!/T0 USE W3ARRYMD, ONLY: PRT2DS
!/T1 USE W3ARRYMD, ONLY: OUTMAT
!/
USE CONSTANTS, ONLY: TPI
USE W3SERVMD, ONLY: EXTCDE
USE W3ODATMD, ONLY: NDSE, IAPROC, NAPROC, NAPERR
USE W3GDATMD, ONLY: NK, NTH, NSPEC, SIG, MAPWN, IC5PARS
USE W3IDATMD, ONLY: INFLAGS2, ICEP1, ICEP2, ICEP3, ICEP4, ICEI
!
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
!/
!/S INTEGER, SAVE :: IENT = 0
!/T0 INTEGER :: ITH
!/T0 REAL :: DOUT(NK,NTH)
!/
REAL :: ICECOEF1, ICECOEF2, ICECOEF3, &
ICECOEF4, ICECONC
REAL, DIMENSION(NK) :: D1D, WN_R, WN_I
REAL :: TWN_R, TWN_I
INTEGER :: IK, IKTH
LOGICAL :: NOICE
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'W3SIC5')
!
! 0. Initializations ------------------------------------------------ /
D = 0.
D1D = 0.
WN_R = 0.
WN_I = 0.
!
ICECOEF1 = 0.
ICECOEF2 = 0.
ICECOEF3 = 0.
ICECOEF4 = 0.
ICECONC = 0.
!
! Set the ice parameters from input
IF (INFLAGS2(-7)) THEN
ICECOEF1 = ICEP1(IX, IY) ! ice thickness h_i
ELSE
IF ( IAPROC .EQ. NAPERR ) &
WRITE (NDSE,1001) 'ICE PARAMETER 1 (HICE)'
CALL EXTCDE(2)
ENDIF
!
IF (INFLAGS2(-6)) THEN
ICECOEF2 = ICEP2(IX, IY) ! effective viscosity of ice η
ELSE
IF ( IAPROC .EQ. NAPERR ) &
WRITE (NDSE,1001) 'ICE PARAMETER 2 (VISC)'
CALL EXTCDE(2)
ENDIF
!
IF (INFLAGS2(-5)) THEN
ICECOEF3 = ICEP3(IX, IY) ! density of ice ρ_i
ELSE
IF ( IAPROC .EQ. NAPERR ) &
WRITE (NDSE,1001) 'ICE PARAMETER 3 (DENS)'
CALL EXTCDE(2)
ENDIF
!
IF (INFLAGS2(-4)) THEN
ICECOEF4 = ICEP4(IX, IY) ! effective shear modulus of ice G
ELSE
IF ( IAPROC .EQ. NAPERR ) &
WRITE (NDSE,1001) 'ICE PARAMETER 4 (ELAS)'
CALL EXTCDE(2)
ENDIF
!
IF (INFLAGS2(4)) ICECONC = ICEI(IX, IY) ! ice concentration
!
! 1. No ice --------------------------------------------------------- /
NOICE = .FALSE.
! Zero ice thickness
! Very small ice thickness may cause problems in POLYROOTS because
! the first coefficient C1 may be very close to zero. So we regard
! cases where hice is less than 0.0001 as no ice.
! IF (ICECOEF1 < ERRTOL) NOICE = .TRUE.
IF (ICECOEF1 < 0.0001) NOICE = .TRUE.
! zero ice concentration
IF (INFLAGS2(4) .AND. ICECONC < ERRTOL) NOICE = .TRUE.
!
! Calculate the decay rate k_i
IF ( NOICE ) THEN
D1D = 0.
!
! 2. Ice ------------------------------------------------------------- /
ELSE
! W3IC5WNCG(WN_R, WN_I, CG, HICE, IVISC, RHOI, ISMODG, HWAT)
CALL W3IC5WNCG(WN_R, WN_I, CG, ICECOEF1, ICECOEF2, &
ICECOEF3, ICECOEF4, DEPTH)
! recall that D=S/E=-2*Cg_{ice}*k_i
! In some cases, the FS model yields very larg Cg_{ice}, which
! subquently may result in numerical failure due to the violation of CFL
! conditions, therefore we still use ice-free group velocity to advect
! wave packets.
!
DO IK = 1, NK
D1D(IK) = -2.0 * CG(IK) * WN_I(IK)
END DO
END IF
!
! 2.1 Fill diagonal matrix
DO IKTH = 1, NSPEC
D(IKTH) = D1D(MAPWN(IKTH))
END DO
S = D * A
!
! ... Test output of arrays
!
!/T0 DO IK=1, NK
!/T0 DO ITH=1, NTH
!/T0 DOUT(IK,ITH) = D(ITH+(IK-1)*NTH)
!/T0 END DO
!/T0 END DO
!
!/T0 CALL PRT2DS (NDST, NK, NK, NTH, DOUT, SIG(1:), ' ', 1., &
!/T0 0.0, 0.001, 'Diag Sice', ' ', 'NONAME')
!
!/T1 CALL OUTMAT (NDST, D, NTH, NTH, NK, 'diag Sice')
!
! Formats
!
1001 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/ &
' ',A,' IS NOT DEFINED IN ww3_shel.inp.')
!/
!/ End of W3SIC5------------------------------------------------------ /
!/
END SUBROUTINE W3SIC5
!/ ------------------------------------------------------------------- /
!/
SUBROUTINE W3IC5WNCG(WN_R, WN_I, CG, HICE, IVISC, RHOI, ISMODG, &
HWAT)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | Q. Liu |
!/ | E. Rogers |
!/ | FORTRAN 90 |
!/ | Last update : 25-Apr-2017 |
!/ +-----------------------------------+
!/
!/ 17-Apr-2016 : Origination ( version 5.10)
!/ ( Q. Liu )
!/ 17-Apr-2016 : Start from W3IC3WNCG_CHENG ( Q. Liu )
!/
!/ 1. Purpose:
! Calculation of complex wavenumber arrays for ice-coupled waves.
!
! This also allows us to use Cg_ice in the advection part of the
! radiative transfer energy equation (RTE). --- abandoned in the end
!
! 2. Method:
! Using the Fox-Squire dispersion relations to get (kr, ki) and
! then get cg by cg = dσ / dk (here dk uses kr)
!
! 3. Parameters:
!
! Parameter list:
! ----------------------------------------------------------------
! Name Type Intent Description
! ----------------------------------------------------------------
! WN_R R.A. I/O the real. part of the wave number
! WN_I R.A. I/O the imag. part of the wave number
! CG R.A. I group velocity (m s^{-1})
! HICE Real. I thickness of ice (m)
! IVISC Real. I viscosity parameter of ice (m^2 s^{-1})
! RHOI Real. I the density of ice (kg m^{-3})
! ISMODG Real. I effecitive shear modulus G of ice (Pa)
! HWAT Real. I water depth
! ----------------------------------------------------------------
! * the intent of WN_R/I must be inout
! * CG is unchanged but still kept here because some legacy reasons.
!
! 4. Subroutines used:
!
! Name Type Module Description
! ----------------------------------------------------------------
! FSDISP Subr. / dispersion relations for ice-coupled waves
! CGINICE5 Subr. / group velocity for given (σ, kr) array
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3SIC5 Subr. Public Ice source term
! W3WAVE Subr. W3WAVEMD WW3 integration
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
USE CONSTANTS, ONLY: TPI
USE W3GDATMD, ONLY: NK, SIG
USE W3ODATMD, ONLY: NDSE, IAPROC, NAPROC, NAPERR
USE W3SERVMD, ONLY: EXTCDE
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
REAL, INTENT(INOUT) :: WN_R(:), WN_I(:)
REAL, INTENT(IN) :: CG(:)
REAL, INTENT(IN) :: HICE, IVISC, RHOI, ISMODG, HWAT
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
REAL, ALLOCATABLE :: SIGMA(:)
INTEGER :: KL, KU, IK
REAL :: TWN_R, TWN_I
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'W3IC5WNCG')
!/
! Initialize SIGMA {in w3gdatmd: SIG (0: NK+1)}
IF (ALLOCATED(SIGMA)) DEALLOCATE(SIGMA); ALLOCATE(SIGMA(SIZE(CG)))
SIGMA = 0.
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
IF ( IAPROC .EQ. NAPERR ) WRITE(NDSE,900) 'W3IC5WNCG'
CALL EXTCDE(3)
END IF
!
! Fox-Squire dispersion
DO IK = KL, KU
! FSDISP(HICE, IVISC, RHOI, ISMODG, HWAT, WT, WNR, WNI)
CALL FSDISP(HICE, IVISC, RHOI, ISMODG, HWAT, TPI/SIGMA(IK), &
TWN_R, TWN_I)
WN_R(IK) = TWN_R
WN_I(IK) = TWN_I
END DO
!
DEALLOCATE(SIGMA)
!
900 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/ &
' Subr. ', A, ': Cannot determine bounds of&
& wavenumber array.'/)
!/
!/ End of W3IC5WNCG -------------------------------------------------- /
!/
END SUBROUTINE W3IC5WNCG
!/ ------------------------------------------------------------------- /
!/
SUBROUTINE FSDISP(HICE, IVISC, RHOI, ISMODG, HWAT, WT, WNR, WNI)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | Q. Liu |
!/ | FORTRAN 90 |
!/ | Last update : 25-Apr-2017 |
!/ +-----------------------------------+
!/
!/ 17-Mar-2016 : Origination ( version 5.10)
!/ ( Q. Liu )
!/ 17-Mar-2016 : Start from the Matlab code `FoxSquire.m` (provided
!/ by Prof. Vernon Squire from University of Otago)
!/ ( Q. Liu )
!/ 25-Apr-2017 : Add more filters ( Q. Liu)
!/
! 1. Purpose :
!
! Calculate the complex wavenumber for waves in ice according to
! the viscoelastic model, i.e., FS model in Mosig et al. (2015)
!
! 2. Method :
! Solving the dispersion relations of FS model (Eq. (20) and (24)
! in Mosig et al. (2015))
!
! Reference:
! Mosig, J. E. M., F. Montiel, and V. A. Squire (2015),
! Comparison of viscoelastic-type models for ocean wave attenuation
! in ice-covered seas, J. Geophys. Res. Oceans, 120, 6072–6090,
! doi:10.1002/2015JC010881.
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! Name Type Intent Description
! ----------------------------------------------------------------
! HICE Real. IN thickness of ice (m)
! IVISC Real. IN viscosity parameter of ice (m^2 s^{-1})
! RHOI Real. IN the density of ice (kg m^{-3})
! ISMODG Real. IN effecitive shear modulus G of ice (Pa)
! HWAT Real. IN water depth
! WT Real. IN wave period (s; 1/freq)
! WNR Real. Out the real. part of the wave number
! WNI Real. Out the imag. part of the wave number
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! POLYROOTS Subr. / Find the roots of a general polynomial
! NR_ROOT Func. / Newton-Raphson root finding
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3IC5WNCG Subr. / Wavenumber and group velocity of ice-
! coupled waves
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! See Format 1000, 1001, 1002
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!/
USE CONSTANTS, ONLY: GRAV, TPI
USE W3DISPMD, ONLY: WAVNU1
USE W3SERVMD, ONLY: EXTCDE
USE W3ODATMD, ONLY: NDSE, IAPROC, NAPROC, NAPERR
USE W3GDATMD, ONLY: IC5PARS
USE W3GSRUMD, ONLY: W3INAN
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
REAL, INTENT(IN) :: HICE, IVISC, RHOI, ISMODG, HWAT, WT
REAL, INTENT(OUT) :: WNR, WNI
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!
REAL :: IC5MINIG, IC5MINWT, IC5MAXKRATIO, &
IC5MAXKI, IC5MINHW
REAL :: TISMODG, TWT, TRATIO, THW
REAL, PARAMETER :: NU = 0.3, RHOW = 1025.
COMPLEX :: GV, C1
REAL :: SIGMA, C2, WNO, CGO, THKH, &
RTRL(5), RTIM(5), RTANG(5)
INTEGER :: IREAL
COMPLEX(KDPC) :: GUESS, CROOT, C1D
REAL(KDP) :: C2D, HWATD
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'FSDISP')
! Note, same as W3IC3WNCG_xx in w3sic3md :
! HICE → ICE1
! IVISC → ICE2
! RHOI → ICE3
! ISMODG → ICE4
! 0. Initializations ------------------------------------------------ *
! Set limiters
!
! When G = 0, the FS method does not provide a solution. It is not
! unexpected because the FS model is originally devised as a
! thin elastic plate model in which elasticity is necessary.
!
! The FS algorithm may also have issues for very short wave periods,
! shallow waters and low G (e.g., T~3 s, d~10 m, hi~0.5 m, G<10^6 Pa)
!
IC5MINIG = IC5PARS(1) ! Minimum G
IC5MINWT = IC5PARS(2) ! Minimum T
IC5MAXKRATIO = IC5PARS(3) ! Maximum k_{ow}/k_r
IC5MAXKI = IC5PARS(4) ! Maximum k_i
IC5MINHW = IC5PARS(5) ! Minimum d
!
TISMODG = MAX(IC5MINIG, ISMODG)
TWT = MAX(IC5MINWT, WT)
THW = MAX(IC5MINHW, HWAT)
!
! G <= 0. is not allowed
IF (ABS(TISMODG) < ERRTOL) THEN
IF ( IAPROC .EQ. NAPERR ) WRITE(NDSE, 1000) 'FSDISP'
CALL EXTCDE(1)
END IF
!
! σ = 2π / T
SIGMA = TPI / TWT
!
! Complex shear modulus Gv = G - i σ ρ η
GV = CMPLX(TISMODG, -1. * SIGMA * RHOI * IVISC)
!
! -------------------------------------------------------------------- *
! Note that Eq. (24) in Mosig et al. (2015) can be written like below:
! (c1 * k^5 + c2 * k) * tanh(HWAT*k) - 1 = 0
! Most Important part of this module --------------------------------- *
C1 = GV * HICE**3. / (6. * RHOW * SIGMA**2.)
!
! To be divided by (1-NU) or multiplied by (1+NU) ??
! Beam model: then multiplied by (1+ν)
! Plate model: then divided by (1-ν)
! The beam version is more theoretically (J.E.M. Mosig, personal
! communication, 2016), although there is only very marginal difference
! between this two version as (1+NU = 1.3 and 1/(1-NU) ~ 1.4)
C1 = C1 * (1+NU)
! C1 = C1 / (1-NU)
!
! C2
C2 = GRAV / SIGMA**2. - RHOI * HICE / RHOW
!
! Use the dispersion in open water to get an approximation of
! tanh(HWAT * k). We can also roughly use the dispersion in deep
! water case, that is tanh(HWAT*k) ~ 1.
! Wavenumber in the open water
! WAVNU1(SI, H, K, CG)
CALL WAVNU1(SIGMA, THW, WNO, CGO)
THKH = TANH(WNO * THW)
!
! Get the first guess of the complex wavenumber
CALL POLYROOTS(6, &
(/REAL(REAL(C1))*THKH, 0., 0., 0., C2*THKH, -1./),&
RTRL, RTIM)
RTANG = ATAN2(RTIM, RTRL)
!
! There should only be one real root in RTRL + i * RTIM because in
! this case (ivisc=0) the original viscoelastic-type model reduced to
! the thin elastic plate model which has only one real solution.
! Find its index ...
!
IREAL = MINLOC(ABS(RTANG), DIM=1)
IF (RTRL(IREAL) <= 0. .OR. ABS(RTIM(IREAL)) > ERRTOL) THEN
IF ( IAPROC .EQ. NAPERR ) WRITE(NDSE, 1001) 'FSDISP'
CALL EXTCDE(2)
END IF
!
! Get the first guess for iteration
GUESS = RTRL(IREAL) * EXP(CMPLX(0., 1E-6))
!
! Newton-Raphson method
! Turn c1, c2, hwat to be double
C1D = C1; C2D = C2; HWATD = THW
CROOT = NR_ROOT(C1D, C2D, HWATD, GUESS)
WNR = REAL(REAL(CROOT))
WNI = REAL(AIMAG(CROOT))
!
! RATIO Check
! Using the ratio k0 / kr as a basic check for the reliability of
! FSDISP. The FS dispersion relation can give a very different kr from
! k0, especially for small wave periods (k0/kr is as high as 100).
! From my tests, using IC5MAXKRATIO = 1000. can basically detect most
! spurious solutions (although not all of them)
!
! ISNAN Check
! Common ways used are:
! NAN = SQRT(-1.) or
! a /= a then a is NaN or
! ISNAN func (supported by gfortran & ifort)
! --- ISNAN -> W3INAN because ISNAN is not supported by pgi
! For very few cases, we can get nan | negative ki | kr
!
! (N.B.) NaN problem solved by using CMPLX_TANH2
!
TRATIO = WNO / WNR
IF (W3INAN(WNR) .OR. W3INAN(WNI) .OR. WNR <= 0 .OR. WNI <= 0. &
.OR. TRATIO >= IC5MAXKRATIO) THEN
IF ( IAPROC .EQ. NAPERR ) &
WRITE(NDSE, 1002) 'FSDISP', HICE, IVISC, TISMODG, HWAT, TWT, &
WNO, WNR, WNI
CALL EXTCDE(3)
END IF
!
! Filter high ki
WNI = MIN(IC5MAXKI, WNI)
!
! FORMAT
1000 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/ &
' Subr. ', A, ': Zero shear modulus G is not allowed&
& in the FS viscoelastic model'/)
!
1001 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/ &
' Subr. ', A, ': get a bad first guess'/)
!
1002 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/ &
' -----------------------------------------------------'/&
' Subr. ', A,' : get NaN/NeG/Huge kr or ki for' /&
' -----------------------------------------------------'/&
' Ice thickness : ', F9.1, ' m'/ &
' Ice viscosity : ', E9.2, ' m2/s'/ &
' Ice shear modulus : ', E9.2, ' Pa' / &
' Water depth : ', F9.1, ' m'/ &
' Wave period : ', F10.2, ' s'/ &
' Wave number (Ko) : ', F11.3, ' rad/m'/ &
' Wave number (Kr) : ', F11.3, ' rad/m'/ &
' Attenu. Rate (Ki) : ', E9.2, ' /m'/)
!/
!/ End of FSDISP ----------------------------------------------------- /
!/
END SUBROUTINE FSDISP
!/ ------------------------------------------------------------------- /
!/
SUBROUTINE BALANCING_MATRIX(NMAT, MATRIX)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | Q. Liu |
!/ | FORTRAN 90 |
!/ | Last update : 15-Mar-2016 |
!/ +-----------------------------------+
!/
!/ 15-Mar-2016 : Origination ( version 5.10)
!/ ( Q. Liu )
!/ 15-Mar-2016 : Borrowed from Numerical Recipes in Fortran
!/ ( Q. Liu )
! 1. Purpose :
! Reducing the sensitivity of eigenvalues to rounding errors during
! the execution of some algorithms.
!
! 2. Method :
! The errors in the eigensystem found by a numerical procedure are
! generally proportional to the Euclidean norm of the matrix, that
! is, to the square root of the sum of the squares of the elements
! (sqrt(sum(a_{i, j} ** 2.)). The idea of balancing is to use
! similarity transformations to make corresponding rows and columns
! of the matrix have comparable norms, thus reducing the overall
! norm of the matrix while leaving the eigenvalues unchanged. Note
! that the symmetric matrix is already balanced.
!
! The output is matrix that is balanced in the norm given by
! summing the absolute magnitudes of the matrix elements(
! sum(abs(a_{i, j})) ). This is more efficient than using the
! Euclidean norm, and equally effective: a large reduction in
! one norm implies a large reduction in the other.
!
! For the details of this method, please refer to
! 1) Numerical Recipes in Fortran 77 (Volume 1, 2nd Edition)
! [Chapter 11.5 / subroutine balanc]
! 2) Numerical Recipes in Fortran 90 (Volume 2)
! [Chapter B11 / subroutine balanc]
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! Name Type Intent Description
! ----------------------------------------------------------------
! NMAT Int. I The size of one dimension of MATRIX
! MATRIX R.A. I/O A matrix with the shape (NMAT, NMAT)
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing (!/S switch).
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! POLYROOTS Subr. / Find the roots of polynomials
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks:
! Balancing only needs marginal computational efforts but can
! substantially improve the accuracy of the eigenvalues computed
! for a badly balanced matrix. It is therefore recommended that
! you always balance nonsymmetric matrices.
!
! Given a (NMAT, NMAT) MATRIX, this routine replaces it by a
! balanced matrix with identical eigenvalues. A symmetric matrix is
! already balanced and is unaffected by this procedure.
!
! 8. Structure :
!
! See the source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
INTEGER, INTENT(IN) :: NMAT
REAL, INTENT(INOUT) :: MATRIX(NMAT, NMAT)
!/ ------------------------------------------------------------------- /
!/ Local parameter
!/S INTEGER, SAVE :: IENT = 0
! the parameter radx is the machine's floating-point radix
REAL, PARAMETER :: RADX = RADIX(MATRIX), &
SQRADX = RADX ** 2
INTEGER :: I, LAST
REAL :: C, F, G, R, S
!/ ------------------------------------------------------------------- /
!/S CALL STRACE (IENT, 'BALANCING_MATRIX')
!
DO
LAST = 1
DO I = 1, NMAT
! Calculate row and column norms
C = SUM( ABS(MATRIX(:, I)) ) - MATRIX(I, I)
R = SUM( ABS(MATRIX(I, :)) ) - MATRIX(I, I)
! If both are non-zero
IF (C /= 0.0 .AND. R /= 0.0) THEN
! Find the integer power of the machine radix that comes closest to
! balancing the matrix (get G, F from C, R)
G = R / RADX
F = 1.0
S = C + R
DO
IF (C >= G) EXIT
F = F * RADX
C = C * SQRADX
END DO
!
G = R * RADX
DO
IF (C <= G) EXIT
F = F / RADX
C = C / SQRADX
END DO
!
IF ( (C+R)/F < 0.95*S) THEN
LAST = 0
G = 1.0 / F
! Apply similarity tranformation
MATRIX(I, :) = MATRIX(I, :) * G
MATRIX(:, I) = MATRIX(:, I) * F
END IF
END IF
END DO
IF (LAST /= 0) EXIT
END DO
!/
!/ End of subroutine BALANCING_MATRIX -------------------------------- /
!/
END SUBROUTINE BALANCING_MATRIX
!/ ------------------------------------------------------------------- /
!/
SUBROUTINE EIG_HQR (NMAT, HMAT, EIGR, EIGI)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | Q. Liu |
!/ | FORTRAN 90 |
!/ | Last update : 17-Mar-2016 |
!/ +-----------------------------------+
!/
!/ 16-Mar-2016 : Origination ( version 5.10)
!/ ( Q. Liu )
!/ 16-Mar-2016 : Borrowed from Numerical Recipes in Fortran
!/ ( Q. Liu )
!/ 17-Mar-2016 : Update the NR code v2.08 to v2.10 ( Q. Liu )
!/
! 1. Purpose :
!
! When we calculate the eigenvalues of a general matrix, we first
! reduce the matrix to a simpler form (e.g., Hessenberg form) and
! then we perform the iterative procedures.
!
! A upper Hessenberg matrix has zeros everywhere below the diagnal
! except for the first subdiagonal row. For example, in the 6x6
! case, the non-zero elements are:
! |x x x x x x|
! |x x x x x x|
! | x x x x x|
! | x x x x|
! | x x x|
! | x x|
!
! This subroutine uses QR algorithm to get the eigenvalues of a
! Hessenberg matrix. So make sure the input array HMAT is a
! Hessenberg-type matrix.
!
! 2. Method :
! QR algorithm for real Hessenberg matrices.
! (I did not understand this algorithm well, so I could not give
! any detailed explanations)
!
! For the details of this HQR method, please refer to
! 1) Numerical Recipes in Fortran 77 (Volume 1, 2nd Edition)
! [Chapter 11.6 / subroutine hqr]
! 2) Numerical Recipes in Fortran 90 (Volume 2)
! [Chapter B11 / subroutine hqr]
!
! Note that there is a bug in the `hqr` subroutine in NR v2.08.
! See http://numerical.recipes/latest-known-bugs.html. Please use
! the updated code in NR v2.10.
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! Name Type Intent Description
! ----------------------------------------------------------------
! NMAT Int. I the size of one dimension of HMAT
! HMAT R.A. I/O the Hessenberg-type matrix (NMAT, NMAT)
! EIGR R.A. O the real part of the N eigenvalues
! EIGI R.A. O the imag part of the N eigenvalues
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! POLYROOTS Subr. / Find the roots of polynomials
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! See Format 1001
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!/
USE W3SERVMD, ONLY: EXTCDE
USE W3ODATMD, ONLY: NDSE, IAPROC, NAPROC, NAPERR
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
INTEGER, INTENT(IN) :: NMAT
REAL, INTENT(INOUT) :: HMAT(NMAT, NMAT)
REAL, INTENT(OUT) :: EIGR(NMAT), EIGI(NMAT)
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
!/S INTEGER, SAVE :: IENT = 0
!/
INTEGER :: I, ITS, K, L, M, NN, MNNK, IDIAG
REAL :: ANORM, P, Q, R, S, T, U, V, W, X, Y, Z
REAL :: PP(NMAT)
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'EIG_HQR')
!
! Compute matrix norm for possible use in locating single small
! subdiagonal element.
!
! Note the speciality of Hessenberg matrix :
! Elements below the diagonal are zeros except for the first
! subdiagonal row. It might be more accurate if we use a mask array
! to mask all zero elments.
!
ANORM = SUM(ABS(HMAT))
NN = NMAT
! Gets changed only by an exceptional shift.
T = 0.0
! Begin search for next eigenvalue: "do while nn >= 1"
DO
IF (NN < 1) EXIT
ITS=0
! Begin iteration
ITERATE:DO
! Look for single small subdiagonal element.
SMALL: DO L=NN, 2, -1
S = ABS( HMAT(L-1, L-1) ) + ABS( HMAT(L, L) )
! IF (S == 0.0) S = ANORM
IF (ABS(S) < ERRTOL) S = ANORM
! IF ( ABS(HMAT(L, L-1)) + S == S ) THEN
IF ( ABS(HMAT(L, L-1)) < ERRTOL ) THEN
HMAT(L, L-1) = 0.0
EXIT SMALL
END IF
END DO SMALL
X = HMAT(NN, NN)
! One root found
IF (L == NN) THEN
EIGR(NN) = X + T
EIGI(NN) = 0.0
NN=NN-1
! Go back for next eigenvalue
EXIT ITERATE
END IF
Y = HMAT(NN-1, NN-1)
W = HMAT(NN, NN-1) * HMAT(NN-1, NN)
! Two roots found . . .
IF (L == NN-1) THEN
P = 0.5 * (Y - X)
Q = P**2 + W
Z = SQRT( ABS(Q) )
X = X + T
! . . . A real pair . . .
IF (Q >= 0.0) THEN
Z = P + SIGN(Z, P)
EIGR(NN) = X + Z
EIGR(NN-1) = EIGR(NN)
IF (Z /= 0.0) EIGR(NN) = X - W/Z
EIGI(NN) = 0.0
EIGI(NN-1) = 0.0
! . . . A complex pair
ELSE
EIGR(NN) = X + P
EIGR(NN-1) = EIGR(NN)
EIGI(NN) = Z
EIGI(NN-1) = -Z
END IF
NN=NN-2
! GO BACK FOR NEXT EIGENVALUE.
EXIT ITERATE
END IF
! NO ROOTS FOUND. CONTINUE ITERATION.
IF (ITS == 30) THEN
IF ( IAPROC .EQ. NAPERR ) WRITE(NDSE, 1001) 'EIG_HQR'
CALL EXTCDE(2)
END IF
! FORM EXCEPTIONAL SHIFT.
IF (ITS == 10 .OR. ITS == 20) THEN
T = T + X
! Add -X to the diagonal of HMAT
DO IDIAG = 1, NN
HMAT(IDIAG, IDIAG) = HMAT(IDIAG, IDIAG) + (-X)
END DO
S = ABS(HMAT(NN, NN-1)) + ABS(HMAT(NN-1, NN-2))
X = 0.75 * S
Y = X
W = -0.4375 * S**2
END IF
ITS = ITS + 1
! Form shift and then look for 2 consecutive small subdiagonal elements.
DO M = NN-2, L, -1
Z = HMAT(M, M)
R = X - Z
S = Y - Z
! Equation (11.6.23).
P = (R * S - W) / HMAT(M+1, M) + HMAT(M, M+1)
Q = HMAT(M+1, M+1) - Z - R - S
R = HMAT(M+2, M+1)
! Scale to prevent overflow or underflow
S = ABS(P) + ABS(Q) + ABS(R)
P = P / S
Q = Q / S
R = R / S
IF (M == L) EXIT
U = ABS( HMAT(M, M-1) ) * ( ABS(Q) + ABS(R) )
V = ABS(P) * ( ABS(HMAT(M-1, M-1)) + ABS(Z) + &
ABS( HMAT(M+1, M+1) ))
! Equation (11.6.26)
IF (U+V == V) EXIT
END DO
DO I= M+2, NN
HMAT(I, I-2) = 0.0
IF (I /= M+2) HMAT(I, I-3)=0.0
END DO
! Double QR step on rows L to NN and columns M to NN
DO K=M, NN-1
IF (K /= M) THEN
! Begin setup of householder vector
P = HMAT(K, K-1)
Q = HMAT(K+1, K-1)
R = 0.0
IF (K /= NN-1) R = HMAT(K+2, K-1)
X = ABS(P) + ABS(Q) + ABS(R)
IF (X /= 0.0) THEN
! Scale to prevent overflow or underflow
P = P / X
Q = Q / X
R = R / X
END IF
END IF
S = SIGN(SQRT(P**2 + Q**2 + R**2), P)
IF (S /= 0.0) THEN
IF (K == M) THEN
IF (L /= M) HMAT(K, K-1) = -HMAT(K, K-1)
ELSE
HMAT(K, K-1) = -S * X
END IF
! Equations (11.6.24).
P = P + S
X = P / S
Y = Q / S
Z = R / S
Q = Q / P
! READY FOR ROW MODIFICATION.
R = R / P
PP(K:NN) = HMAT(K, K:NN) + Q * HMAT(K+1, K:NN)
IF (K /= NN-1) THEN
PP(K:NN) = PP(K:NN) + R * HMAT(K+2, K:NN)
HMAT(K+2, K:NN) = HMAT(K+2, K:NN) - &
PP(K:NN)*Z
END IF
HMAT(K+1, K:NN) = HMAT(K+1, K:NN) - PP(K:NN) * Y
HMAT(K, K:NN) = HMAT(K, K:NN) - PP(K:NN) * X
! COLUMN MODIFICATION.
MNNK = MIN(NN, K+3)
PP(L:MNNK) = X * HMAT(L:MNNK, K) + Y * &
HMAT(L:MNNK, K+1)
IF (K /= NN-1) THEN
PP(L:MNNK) = PP(L:MNNK) + Z*HMAT(L:MNNK, K+2)
HMAT(L:MNNK, K+2) = HMAT(L:MNNK, K+2) - &
PP(L:MNNK) * R
END IF
HMAT(L:MNNK, K+1) = HMAT(L:MNNK, K+1) - &
PP(L:MNNK) * Q
HMAT(L:MNNK, K) = HMAT(L:MNNK, K) - PP(L:MNNK)
END IF
END DO
! GO BACK FOR NEXT ITERATION ON CURRENT EIGENEND DO VALUE.
END DO ITERATE
END DO
!
! Formats
1001 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/ &
' Subr. ', A, ': TOO MANY ITERATIONS'/)
!/ ------------------------------------------------------------------- /
!/
!/ End of EIG_HQR ---------------------------------------------------- /
!/
END SUBROUTINE EIG_HQR
!/ ------------------------------------------------------------------- /
!/
SUBROUTINE POLYROOTS(NPC, PCVEC, RTRL, RTIM)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | Q. Liu |
!/ | FORTRAN 90 |
!/ | Last update : 16-Mar-2016 |
!/ +-----------------------------------+
!/
!/ 16-Mar-2016 : Origination ( version 5.10)
!/ ( Q. Liu )
!/ 16-Mar-2016 : Started from Numerical Recipes in Fortran
!/ ( Q. Liu )
!/
! 1. Purpose :
!
! Find the roots of arbitrary polynomials through finding the
! eigenvalues of companion matrix.
!
! 2. Method :
! Suppose we have a general polynomial, which reads
! P(x) = c_n * x^n + c_{n-1} * x^{n-1} + ... + c_1 * x + c_0
!
! Then finding the roots of P(x) is equivalent to find the eigen-
! values of the special n x n companion matrix A
! | -c_{n-1}/c_n -c_{n-2}/c_n ... -c_1/c_n -c_0/c_n |
! | 1 0 ... 0 0 |
! A = | 0 1 ... 0 0 |
! | : : : : |
! | 0 0 1 0 |
!
! In fact, P(x) is the characteristic polynomial of matrix A, i.e.,
! P(x) = det|A-xI| and x is the eigenvalues of A (this is a
! Hessenberg matrix).
!
! In this subrountine, we will use the two subroutines above
! (BALANCING_MATRIX & EIG_HQR) to get the complex eigenvalues of
! an abitrary Hessenberg matrix
!
! For the details of this method, please refer to
! 1) Numerical Recipes in Fortran 77 (Volume 1, 2nd Edition)
! [Chapter 9.5 / Eigenvalue Methods / subroutine zrhqr]
! 2) Numerical Recipes in Fortran 90 (Volume 2)
! [Chapter B9 / subroutine zrhqr]
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! Name Type Intent Description
! ----------------------------------------------------------------
! NPC Int. I The # of the Polynomial coefficients
! (from c_n to c_0)
! PCVEC R.A. I The 1d vector for the Polynomial
! coefficients [c_n, c_{n-1}, ..., c_0]
! RTRL R.A. O The real part of all of the roots
! shape: [NPC-1]
! RTIM R.A. O The real part of all of the roots
! shape: [NPC-1]
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! Name Type Module Description
! ------------------------------------- ---------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! BALANCING_MATRIX Subr. / Balancing matrix
! EIG_HQR Subr. / Finding eigenvalues
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! FSDISP Subr. / Solving the dispersion relations
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! See Format 1001
!
! 7. Remarks :
! The built-in MATLAB function <roots> uses the same method to
! find roots of a general polynomial. But perhaps MATLAB uses
! different methods to find eigenvalues of the companion matrix.
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!/
USE W3SERVMD, ONLY: EXTCDE
USE W3ODATMD, ONLY: NDSE, IAPROC, NAPROC, NAPERR
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
INTEGER, INTENT(IN) :: NPC
REAL, INTENT(IN) :: PCVEC(NPC)
REAL, INTENT(OUT) :: RTRL(NPC-1), RTIM(NPC-1)
!/
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
!/S INTEGER, SAVE :: IENT = 0
REAL :: HESS(NPC-1, NPC-1)
INTEGER :: J
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'POLYROOTS')
!
!
IF (ABS(PCVEC(1)) < ERRTOL) THEN
IF ( IAPROC .EQ. NAPERR ) WRITE(NDSE, 1001) 'POLYROOTS'
CALL EXTCDE(2)
END IF
!
! Generate the Hessenberg matrix
HESS = 0.
HESS(1, :) = -1 * PCVEC(2:) / PCVEC(1)
DO J = 1, NPC-2
HESS(J+1, J) = 1.
END DO
! Balancing the matrix HESS
CALL BALANCING_MATRIX(NPC-1, HESS)
! Eigenvalues of the matrix HESS
CALL EIG_HQR(NPC-1, HESS, RTRL, RTIM)
! Formats
1001 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/ &
' Subr. ', A, ': the coeff. of x^n must not be 0'/)
!/
!/ End of POLYROOTS -------------------------------------------------- /
!/
END SUBROUTINE POLYROOTS
!/ ------------------------------------------------------------------- /
!/
FUNCTION NR_CORR(K, C1, C2, H)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | Q. Liu |
!/ | FORTRAN 90 |
!/ | Last update : 24-Mar-2016 |
!/ +-----------------------------------+
!/
!/ 18-Mar-2016 : Origination. ( version 5.10 )
!/ ( Q. Liu )
!/ 18-Mar-2016 : Start from the Matlab code `FoxSquire.m` (provided
!/ by Prof. Vernon Squire from University of Otago)
!/ ( Q. Liu )
!/ 24-Mar-2016 : Adding the cmplx_sinh/cosh/tanh ( Q. Liu )
!
! 1. Purpose :
!
! Calculate the corrected term in the Newton-Raphson root-finding
! method (Must use double precision)
!
! 2. Method :
! Suppose we want to find the root of f(x) = 0, then according to
! the Newton-Raphson method, the root is iteratively updated by the
! formula below:
!
! x_{i+1} = x_{i} - f(x_{i}) / f'(x_{i}),
!
! where f'(x) denotes the derivative of f(x). In this function,
! our f(x) reads (see also subr. FSDISP)
!
! f(x) = (c1 * k**4 + c2) * k * tanh(kH) -1
!
! we finally will get the Newton-Raphson correted term, i.e.,
!
! dx = f(x_{i}) / f'(x_{i})
!
! For the details of this method, please refer to
! 1) Numerical Recipes in Fortran 77 (Volume 1, 2nd Edition)
! Chapter 9.4
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! Name Type Intent Description
! ----------------------------------------------------------------
! K CMPL.(D) I complex wave number
! C1 CMPL.(D) I C1 in FSDISP
! C2 Real.(D) I C2 in FSDISP
! H Real.(D) I water depth
! NR_CORR CMPL.(D) O Newton-Raphson corrected term (DK)
! ----------------------------------------------------------------
! * (D) means double precision
!
! 4. Subroutines used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! CMPLX_SINH Func. / sinh for complex var.
! CMPLX_COSH Func. / cosh for complex var.
! CMPLX_TANH2 Func. / tanh for complex var.
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! NR_ROOT Func. / Newton-Raphson root finding.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
COMPLEX(KDPC), INTENT(IN) :: K, C1
REAL(KDP), INTENT(IN) :: C2, H
COMPLEX(KDPC) :: NR_CORR
!/
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
!/S INTEGER, SAVE :: IENT = 0
! A rough value to differentiate deep water case from finite water case
REAL(KDP), PARAMETER :: KH_LIM = 7.5
COMPLEX(KDPC) :: LAM, LAMPR, FV, DF, TKH
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'NR_CORR')
! f(k) = (c1 * k**4 + c2) * k * tanh(k*H) - 1
! = lam * k * tanh(k*H) - 1
!
TKH = K * H
LAM = C1 * K**4 + C2
! the derivative of (lam * k)
LAMPR = 5 * C1 * K**4 + C2
!
IF (REAL(REAL(TKH)) <= KH_LIM) THEN
! KH is small enough
! FV = LAM * K * SINH(K*H) - COSH(K*H)
! DF = LAM * (K*H) * COSH(K*H) + (LAMPR - H) * SINH(K*H)
FV = LAM * K * CMPLX_SINH(TKH) - CMPLX_COSH(TKH)
DF = LAM * TKH * CMPLX_COSH(TKH) + (LAMPR-H) * CMPLX_SINH(TKH)
ELSE
! FV = LAM * K * TANH(K*H) - 1
! DF = LAM * K * H + (LAMPR - H) * TANH(K*H)
! DF = LAMPR * TANH(K*H) + LAM * K * H / (COSH(K*H) **2)
FV = LAM * K * CMPLX_TANH2(TKH) - 1
DF = LAMPR * CMPLX_TANH2(TKH) + LAM * TKH * &
(1 - CMPLX_TANH2(TKH) ** 2.)
END IF
!
NR_CORR = FV / DF
!/
!/ End of NR_CORR ---------------------------------------------------- /
!/
END FUNCTION NR_CORR
!/ ------------------------------------------------------------------- /
!/
FUNCTION NR_ROOT(C1, C2, H, GUESS)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | Q. Liu |
!/ | FORTRAN 90 |
!/ | Last update : 18-Mar-2016 |
!/ +-----------------------------------+
!/
!/ 18-Mar-2016 : Origination. ( version 5.10 )
!/ ( Q. Liu )
!/ 18-Mar-2016 : Start from the Matlab code `FoxSquire.m` (provided
!/ by Prof. Vernon Squire from University of Otago)
!/ ( Q. Liu )
! 1. Purpose :
!
! The iterative procedure of the Newton-Raphson method
!
! 2. Method :
! See the document of Subr. NR_CORR (Must use double precision)
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! Name Type Intent Description
! ----------------------------------------------------------------
! C1 CMPL.(D) I C1 in FS dipsersion relations
! See the doc. of subr. NR_CORR
! C2 REAL (D) I C2 in FS dipsersion relations
! H REAL (D) I water depth
! GUESS CMPL.(D) I the first guess obtained from POLYROOTS
! NR_ROOT CMPL.(D) O the calculated complex wave number.
! ----------------------------------------------------------------
! * (D) means double precision
!
! 4. Subroutines used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! NR_CORR Func. / Newton-Raphson correction term
! INIT_RANDOM_SEED
! Subr. / Initialize the random seed based on
! the system's time
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! FSDISP Subr. / Solve FS dispersion relations
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!/
USE W3SERVMD, ONLY: EXTCDE
USE W3ODATMD, ONLY: NDSE, IAPROC, NAPROC, NAPERR
USE W3GDATMD, ONLY: IC5PARS
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
COMPLEX(KDPC), INTENT(IN) :: C1, GUESS
REAL(KDP), INTENT(IN) :: C2, H
COMPLEX(KDPC) :: NR_ROOT
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
!/S INTEGER, SAVE :: IENT = 0
COMPLEX(KDPC) :: K0, K1, DK
INTEGER :: ITER
REAL :: TRANVAL
REAL :: IC5MAXITER, IC5RKICK, IC5KFILTER
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'NR_ROOT')
!/ Set parameters
IC5MAXITER = IC5PARS(6)
IC5RKICK = IC5PARS(7) ! 0: False, 1: True
IC5KFILTER = IC5PARS(8)
!
K0 = GUESS
DK = NR_CORR(K0, C1, C2, H)
K1 = K0 - DK
ITER = 0
IF (IC5RKICK > 0.5) CALL INIT_RANDOM_SEED()
!
DO WHILE (ABS(DK) > ERRTOL)
K0 = K1
DK = NR_CORR(K0, C1, C2, H)
K1 = K0 - DK
ITER = ITER + 1
!
! Random kick to avoid converging to evanescent modes
! Note: do not use RAND(1) because it alway gives a same random no.
! The built in function of RAND is not available in <ifort>, use
! random_seed/number instead.
!
! Based on many tests, I found the random kick & the corridor excluded
! from imaginary axis are kind of helpful to avoid spurious solutions.
! However, it may also lead to no solutions returned, especially for
! high G and high T.
!
IF (IC5RKICK > 0.5 .AND. ABS(REAL(K1)) < IC5KFILTER) THEN
! K1 = K1 + 2*RAND(0)
CALL RANDOM_NUMBER(TRANVAL)
K1 = K1 + 2 * TRANVAL
END IF
!
IF (ITER >= IC5MAXITER) THEN
IF ( IAPROC .EQ. NAPERR ) WRITE(NDSE, 1001) 'NR_ROOT'
CALL EXTCDE(1)
END IF
!
END DO
!
NR_ROOT = K1
!
! Formats
1001 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/ &
' Subr. ', A, ': TOO MANY ITERATIONS'/)
!
!/
!/ End of NR_ROOT ---------------------------------------------------- /
!/
END FUNCTION NR_ROOT
!/ ------------------------------------------------------------------- /
!/
FUNCTION CMPLX_SINH(X)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | Q. Liu |
!/ | FORTRAN 90 |
!/ | Last update : 24-Mar-2016 |
!/ +-----------------------------------+
!/
!/ 24-Mar-2016 : Origination. ( version 5.10 )
!/ ( Q. Liu )
!/
! 1. Purpose :
!
! For a number of compilers, the built-in function sinh, cosh and
! tanh do not support the complex inputs. So here I write an
! external one.
!
! 2. Method :
!
! sinh(x) = (e**x - e**(-x)) / 2 (The built in function exp supports
! complex input)
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! Name Type Intent Description
! ----------------------------------------------------------------
! X CMPL(D) I a double-precision complex var.
! ----------------------------------------------------------------
! * Note, this subr. will be only called by NR_CORR,
! so for simplicity, I only use double-precision complex var.
! as input.
!
! 4. Subroutines used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! NR_CORR Subr. / Newton-Raphson correction term.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
COMPLEX(KDPC), INTENT(IN) :: X
COMPLEX(KDPC) :: CMPLX_SINH
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/S CALL STRACE (IENT, 'CMPLX_SINH')
!/
CMPLX_SINH = (EXP(X) - EXP(-X)) * 0.5
!/
!/ End of CMPLX_SINH ------------------------------------------------- /
!/
END FUNCTION CMPLX_SINH
!/ ------------------------------------------------------------------- /
!/
FUNCTION CMPLX_COSH(X)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | Q. Liu |
!/ | FORTRAN 90 |
!/ | Last update : 24-Mar-2016 |
!/ +-----------------------------------+
!/
!/ 24-Mar-2016 : Origination. ( version 5.10 )
!/ ( Q. Liu )
!/
! 1. Purpose :
!
! For a number of compilers, the built-in function sinh, cosh and
! tanh do not support the complex inputs. So here I write an
! external one.
!
! 2. Method :
!
! cosh(x) = (e**x + e**(-x)) / 2 (The built in function exp supports
! complex input)
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! Name Type Intent Description
! ----------------------------------------------------------------
! X CMPL(D) I a double-precision complex var.
! ----------------------------------------------------------------
! * Note, this subr. will be only called by NR_CORR,
! so for simplicity, I only use double-precision complex var.
! as input.
!
! 4. Subroutines used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! NR_CORR Subr. / Newton-Raphson correction term.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
COMPLEX(KDPC), INTENT(IN) :: X
COMPLEX(KDPC) :: CMPLX_COSH
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/S CALL STRACE (IENT, 'CMPLX_COSH')
!/
CMPLX_COSH = (EXP(X) + EXP(-X)) * 0.5
!/
!/ End of CMPLX_COSH ------------------------------------------------- /
!/
END FUNCTION CMPLX_COSH
!/ ------------------------------------------------------------------- /
!/
FUNCTION CMPLX_TANH2(X)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | Q. Liu |
!/ | FORTRAN 90 |
!/ | Last update : 24-Mar-2016 |
!/ +-----------------------------------+
!/
!/ 24-Mar-2016 : Origination. ( version 5.10 )
!/ ( Q. Liu )
!/
! 1. Purpose :
! We may encounter overflow error for the above tanh function as kh
! becomes huge. This is another version of tanh function
!
! 2. Method :
!
! See https://en.wikipedia.org/wiki/Hyperbolic_function
! tanh(x) = (exp(x) - exp(-x)) / (exp(x) + exp(-x))
! = (1 - exp(-2x)) / (1 + exp(-2x))
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! Name Type Intent Description
! ----------------------------------------------------------------
! X CMPL(D) I a double-precision complex var.
! ----------------------------------------------------------------
! * Note, this subr. will be only called by NR_CORR, so for
! simplicity, I only use double-precision complex var. as input.
!
! 4. Subroutines used :
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! CMPLX_SINH Func. / sinh for complex var.
! CMPLX_COSH Func. / cosh for complex var.
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! NR_CORR Subr. / Newton-Raphson correction term.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
! Calculating tanh in this way may have problems when x ->
! -inf. But in our cases x is alway >0.
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
COMPLEX(KDPC), INTENT(IN) :: X
COMPLEX(KDPC) :: CMPLX_TANH2
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/S CALL STRACE (IENT, 'CMPLX_TANH2')
!/
CMPLX_TANH2 = (1 - EXP(-2*X)) / (1 + EXP(-2*X))
!/
!/ End of CMPLX_TANH2 ------------------------------------------------ /
!/
END FUNCTION CMPLX_TANH2
!/
!/ ------------------------------------------------------------------- /
!/
SUBROUTINE INIT_RANDOM_SEED()
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | Q. Liu |
!/ | FORTRAN 90 |
!/ | Last update : 24-Mar-2016 |
!/ +-----------------------------------+
!/
!/ 24-Mar-2016 : Origination. ( version 5.10 )
!/ ( Q. Liu )
!/ 24-Mar-2016 : Borrowed from Fortran Wiki ( Q. Liu )
!
! 1. Purpose :
!
! Initialize the random seed based on the system's time.
!
! 2. Method :
!
! See http://fortranwiki.org/fortran/show/random_seed
!
! 3. Parameters :
!
! 4. Subroutines used :
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! NR_ROOT Func. / Newton-Raphson root finding.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
!/S USE W3SERVMD, ONLY: STRACE
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
!/S INTEGER, SAVE :: IENT = 0
!/
INTEGER :: I, N, CLOCK
INTEGER, DIMENSION(:), ALLOCATABLE :: SEED
!/ ------------------------------------------------------------------- /
!/S CALL STRACE (IENT, 'INIT_RANDOM_SEED')
!/
CALL RANDOM_SEED(SIZE = N)
ALLOCATE(SEED(N))
!
CALL SYSTEM_CLOCK(COUNT=CLOCK)
!
SEED = CLOCK + 37 * (/ (I - 1, I = 1, N) /)
CALL RANDOM_SEED(PUT = SEED)
!
DEALLOCATE(SEED)
!/
!/ End of INIT_RANDOM_SEED ------------------------------------------- /
!/
END SUBROUTINE INIT_RANDOM_SEED
!/
!/ ------------------------------------------------------------------- /
!/
!/ End of module W3SIC5MD -------------------------------------------- /
!/
END MODULE W3SIC5MD
!/ ------------------------------------------------------------------- /