#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 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 , 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 !/ ------------------------------------------------------------------- /