#include "w3macros.h"
!/ ------------------------------------------------------------------- /
      MODULE W3SIS2MD
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           F. Ardhuin              |
!/                  |           P. Nicot                |
!/                  |           C. Sevigny              |
!/                  |           G. Boutin               |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         21-Jan-2018 |
!/                  +-----------------------------------+
!/
!/    For updates see W3SID1 documentation.
!/
!  1. Purpose :
!
!     Floe-size dependant scattering of waves in the marginal ice zone based on tabulated
!     scattering coefficients for a semi-infinite ice sheet. See papers
!     by Dumont et al. (JGR 2011) and Williams et al. (OM 2013)
!     combined with flexural dissipation and ice break-up.
!
!  2. Variables and types :
!
!  3. Subroutines and functions :
!
!      Name      Type  Scope    Description
!     ----------------------------------------------------------------
!      W3SIS2    Subr. Public   Ice scattering term.
!     ----------------------------------------------------------------
!
!  4. Subroutines and functions used :
!
!     See subroutine documentation.
!
!  5. Remarks :
!
!  6. Switches :
!
!     See subroutine documentation.
!
!  7. Source code :
 
!/
!/ ------------------------------------------------------------------- /
!/
      INTEGER ,    PARAMETER             :: NTHICK = 20, NICED = 500
      REAL, PARAMETER                    :: FRAGILITY = 0.9
      REAL                               :: THICK1 = 0.1, DTHICK = 0.25
      REAL                               :: ICEDMIN      ! minimum floe diameter
      REAL                               :: ICEDAVETAB(NICED)
      REAL, DIMENSION(:)   , ALLOCATABLE :: SIS2ALPHAS(:,:),SIS2ALPHA2(:,:)
      DOUBLE PRECISION, ALLOCATABLE,DIMENSION(:,:)     :: IS2EIGVEC, IS2SCATMAT
      DOUBLE PRECISION    , ALLOCATABLE,DIMENSION(:)   :: IS2EIGVAL
      PRIVATE :: SIS2ALPHAS, SIS2ALPHA2
      PUBLIC :: IS2EIGVEC, IS2EIGVAL
!
      PUBLIC :: W3SIS2, INSIS2, W3RPWNICE
      PRIVATE :: FINDROOTS_NR, W3FSD_DAVE, FUNCD_FVAL, FUNCD_FDERIV
!/
      CONTAINS
!/ ------------------------------------------------------------------- /
      SUBROUTINE INSIS2
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |      P. Nicot & F. Ardhuin        |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         21-Jan-2018 |
!/                  +-----------------------------------+
!/
!/    01-Apr-2014 : Creation                            ( version 4.18 )
!/    13-Dec-2015 : Adds diagonalization of scat. matrix( version 5.10 )
!/    21-Jan-2018 : Implements non-isotropic example    ( version 6.04 )
!/
!  1. Purpose :
!
!     Fill tables used for scattering
!
!  2. Method :
!
!     Linear interpolation
!
!  3. Parameters :
!
!     See module documentation.
!
!  4. Error messages :
!
!     - None.
!
!  5. Called by :
!
!     - W3IOGR (initialization after reading mod_def.ww3)
!
!  6. Subroutines used :
!
!     - None
!
!  7. Remarks :
!
!  8. Structure
!
!  9. Switches :
!
!       !/S   Enable subroutine tracing.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
!/
      USE W3GDATMD,   ONLY: SIG, DSIP, NK, NTH, IS2PARS,  &
                            EC2, ES2, ESC, ESIN, ECOS
      USE CONSTANTS,  ONLY: TPI, TPIINV
      USE W3SERVMD, ONLY: DIAGONALIZE
!
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      INTEGER                 :: I, J, K, IND, NFTAB, NROT
      REAL                    :: SIS1HTABLE(20), SIS1FTABLE(25)
      REAL                    :: SIS1ALPHATABLE(NTHICK,25), X
      REAL                    :: SIS1ALPHATABLE2(NTHICK,25)
!/
!/ ------------------------------------------------------------------- /
!/
!
! -------------------------------------------------------------------- /
! 1.  Fills array of reflection as a function of frequency and ice thickness
!
      ALLOCATE(SIS2ALPHAS(NTHICK,NK))
      ALLOCATE(SIS2ALPHA2(NTHICK,NK))
!
! Table of ice thickness for which the reflection was computed
!
      SIS1HTABLE = (/ 0.1, 0.35, 0.6, 0.85, 1.1, 1.35, 1.6, 1.85, 2.1,  2.35, &
                      2.6, 2.85, 3.1, 3.35, 3.6, 3.85, 4.1, 4.35, 4.6, 4.85 /)
      NFTAB = 25
 
!
! Table of frequencies for which the reflection was computed
!
      SIS1FTABLE = (/ 0.0420, 0.04620, 0.050820, 0.0559020, 0.06149220, 0.06764142,0.0744055620000000, &
               0.0818461182, 0.09003073002, 0.099033803022, 0.1089371833242, 0.11983090, 0.13181399, &
               0.144995391, 0.15949493 ,  0.175444423115458, 0.192988865427003, 0.212287751969704, &
               0.233516527, 0.256868, 0.28255499787, 0.310810497658843, 0.341891547424728, 0.376080702167200, 0.413688772383920 /)
IF (IS2PARS(18).LT.0.5) THEN
      SIS1ALPHATABLE = reshape((/  &
        1.78E-007, 2.21E-006, 6.57E-006, 1.31E-005, 2.28E-005, 3.60E-005, 5.23E-005, 7.19E-005, 9.60E-005, 0.0001260665, &
        0.0001621645, 0.0002032202, 0.0002483457, 0.0002987263, 0.0003571903, &
        0.0004247358, 0.00049714, 0.0005689017, 0.0006469729, 0.0007470985, &
        2.74E-007, 3.40E-006, 1.04E-005, 2.10E-005, 3.67E-005, 5.82E-005, &
        8.52E-005, 0.0001179843, 0.0001581719, 0.0002075378,    &
        0.0002663323, 0.0003330095, 0.0004063158, 0.000487845, 0.0005815828, &
        0.0006893214, 0.0008061048, 0.0009252626, 0.0010565619, 0.0012218782, &
        4.16E-007, 5.21E-006, 1.62E-005, 3.32E-005, 5.91E-005, 9.56E-005, 0.0001423817, &
        0.0002001783, 0.0002716622, 0.0003599526, 0.0004656353, &
        0.0005861828, 0.0007194785, 0.0008683691, 0.0010401368, 0.0012386623, 0.0014562639, &
        0.0016819935, 0.0019335276, 0.0022507523, &
        6.25E-007, 7.95E-006, 2.50E-005, 5.28E-005, 9.71E-005, 0.0001615969, 0.0002467311, &
        0.0003540015, 0.0004891301, 0.0006588358, 0.0008648602, &
        0.0011024698, 0.0013675698, 0.0016664349, 0.0020150216, 0.0024225002, 0.0028732348, &
        0.0033437743, 0.0038725363, 0.0045480433, &
        9.32E-007, 1.21E-005, 3.91E-005, 8.60E-005, 0.000165062, 0.0002848805, 0.0004473356,  &
        0.0006565411, 0.0009257279, 0.0012707025, 0.0016966237, &
        0.0021939892, 0.002754291, 0.0033926056, 0.0041466541, 0.0050385513, 0.0060311467, &
        0.0070669538, 0.0082363302, 0.0097527074, &
        1.38E-006, 1.86E-005, 6.30E-005, 0.0001464051, 0.0002950074, 0.0005287398, 0.0008535728, &
        0.0012801516, 0.0018391215, 0.0025674934, 0.0034790152,&
        0.0045539812, 0.0057742581, 0.0071757314, 0.0088467499, 0.0108382999, &
        0.0130588683, 0.0153657347, 0.0179697664, 0.0213785773, &
        2.05E-006, 2.91E-005, 0.0001064725, 0.0002646727, 0.0005601397, 0.0010380532, &
        0.0017148531, 0.0026159282, 0.0038087163, 0.0053744112, &
        0.0073439386, 0.0096752268, 0.0123295445, 0.0153861547, 0.0190373907, 0.0233881203, &
        0.0282220195, 0.0332084504, 0.0388053869, 0.0461221517, &
        3.03E-006, 4.71E-005, 0.0001918323, 0.000513363, 0.0011329164,  &
        0.00214921, 0.0036015483, 0.0055437802, 0.0081104965, 0.0114597841, 0.0156453149, &
        0.0205751599, 0.0261681604, 0.0325784506, 0.040175062, 0.0491364185, &
        0.0589976693, 0.0690831851, 0.0802738154, 0.0946821772, &
        4.53E-006, 8.02E-005, 0.0003720589, 0.001067721, 0.0024215478, &
        0.0046335423, 0.0077820302, 0.0119609444, 0.0173982353, 0.0243482166, &
        0.0328666235, 0.0427569128, 0.0538584883, 0.0664227317, 0.0810452442, &
        0.0979544205, 0.1162816559, 0.1348568677, 0.1551243874, 0.1804646861, &
        6.84E-006, 0.0001465372, 0.0007768361, 0.0023492877, 0.0053624751, &
        0.0101694228, 0.0168827981, 0.0256087065, 0.0366354159, 0.0502611381, &
        0.0664627541, 0.0848577616, 0.1051596929, 0.1276991046, 0.1532607056, &
        0.1820303259, 0.2126339282, 0.2433779188, 0.276250064, 0.3157667291, &
        1.06E-005, 0.0002908265, 0.0017246725, 0.0053239059, 0.0119340853, &
        0.0220417832, 0.035687347, 0.0528302869, 0.073649217, 0.098327901, &
        0.1266340957, 0.1579155214, 0.1917144685, 0.2283923343, 0.2688390962, &
        0.3131039168, 0.359282913, 0.4052034867, 0.4533028057, 0.5090275342, &
        1.68E-005, 0.0006285069, 0.0039687874, 0.0119764661, 0.0257390989, &
        0.0456577018, 0.0713257381, 0.1021551099, 0.1379813612, 0.1787934975, &
        0.2240723105, 0.2727534775, 0.3241066558, 0.3786049119, 0.4374575502, &
        0.5006384418, 0.5653760852, 0.6286093824, 0.6937583927, 0.7684424911, &
        2.82E-005, 0.0014590465, 0.0091291531, 0.0257006983, 0.051998859, &
        0.0879406127, 0.1317806874, 0.181756664, 0.2375759087, 0.2995488732, &
        0.3668698493, 0.4374715634, 0.5099661328, 0.5855233305, 0.6668652222, &
        0.754201342, 0.8422038741, 0.925009729, 1.0096608603, 1.1112763605, &
        5.06E-005, 0.0035218608, 0.0201532593, 0.0509283911, 0.0962607638, &
        0.1553638291, 0.2233707746, 0.2967358894, 0.3765436929, 0.4650074273, &
        0.5607827923, 0.6590228472, 0.756685343, 0.8572097524, 0.9677298862,&
        1.0892101997, 1.2094355521, 1.3151711317, 1.4236588573, 1.5698199322, &
        9.99E-005, 0.0084186275, 0.0412091805, 0.0918067318, 0.1629311719, &
        0.2530998556, 0.3511295005, 0.4513397894, 0.5593550205, 0.6817564883, &
        0.815653431, 0.9499632903, 1.0782538468, 1.2090808019, 1.3588213156, &
        1.5298071338, 1.6953922745, 1.8273584757, 1.9646959304, 2.1824935098, &
        0.0002197702, 0.0188788997, 0.076522568, 0.1517352519, 0.2567566667, &
        0.3885311336, 0.5247441812, 0.6573850155, 0.8007012438, 0.9681737964, &
        1.1534672992, 1.3345908371, 1.4996817772, 1.666222509, 1.8655897771, &
        2.1018122262, 2.3245234997, 2.4819251275, 2.6492217053, 2.9647668115, &
        0.0005376597, 0.0381101752, 0.1293683129, 0.2359311434, 0.3874916537, &
        0.576671054, 0.7645781239, 0.9404856491, 1.1305780655, 1.356723506, &
        1.6068566622, 1.844140891, 2.0502598744, 2.2550088512, 2.5085040228, &
        2.815733869, 3.096422752, 3.2715199709, 3.462031842, 3.8799520309, &
        0.0014123358, 0.0685490628, 0.203577783, 0.354894301, 0.5714749709, &
        0.8380065625, 1.0957796855, 1.3307089678, 1.5810360789, 1.8762692254,&
        2.197628874, 2.4940515596, 2.7425260367, 2.9852393578, 3.2886471617, &
        3.6568112069, 3.9833963075, 4.1693103781, 4.3750188868, 4.868451483, &
        0.0037199175, 0.1126937266, 0.3066036571, 0.5225755275, 0.8226663225, &
        1.1819668289, 1.5247740065, 1.8332140124, 2.152798111, 2.5177593097, &
        2.9057350396, 3.2594397982, 3.5553142762, 3.8420910788, 4.1929464318,&
         4.6097641365, 4.975238806, 5.1850270882, 5.4153321104, 5.9516719058, &
        0.009037082, 0.1761561171, 0.4486599083, 0.747307698, 1.1406376531, &
        1.5966776005, 2.0299434584, 2.4186982894, 2.8096987013, 3.2399866752, &
        3.6899433755, 4.105747098, 4.4658248514, 4.8172532951, 5.2300901821, &
        5.7046450426, 6.1270633187, 6.3975380116, 6.6851833028, 7.2647998139, &
        0.0191118494, 0.2660897811, 0.6359174284, 1.0316197503, 1.5313672569, &
        2.0961171527, 2.6294074335, 3.1054178673, 3.5736395046, 4.0762156975, &
        4.5979660806, 5.0883063988, 5.5262800229, 5.9567260476, 6.4472029911, &
        6.9971641707, 7.4947435133, 7.841785986, 8.2023388524, 8.8466621683, &
        0.0352335589, 0.3840118688, 0.8727819679, 1.3979077636, 2.0496249898, &
        2.7683726012, 3.4325758193, 4.0133546648, 4.572769913, 5.1624475081, &
        5.7663362962, 6.3283500032, 6.8264493734, 7.3114532863, 7.8576036504, &
        8.4626783134, 9.0041860776, 9.378581111, 9.7670030485, 10.46036883, &
        0.0588123086, 0.5407217302, 1.1866725247, 1.8744101474, 2.7029832853, &
        3.586993579, 4.3802119773, 5.0573372156, 5.6991092734, 6.3694458919, &
        7.0521552001, 7.6857035561, 8.2469182713, 8.7926955021, 9.4042992778, &
        10.0780089999, 10.6798375161, 11.0993158393, 11.5359791079, 12.3054058274, &
        0.0912925253, 0.7587053794, 1.5933132286, 2.4467380656, 3.4467575836, &
        4.4845019618, 5.3910117536, 6.1520041613, 6.8749329856, 7.6395860278, &
        8.4275368684, 9.1660844666, 9.8268448897, 10.4746437466, 11.2031746957, &
        12.0078539018, 12.7326579385, 13.2484048884, 13.7868789237, 14.7164075718, &
        0.1320035456, 1.258347597, 2.2697363962, 2.8626691529, 3.4800285532, &
        4.1591450976, 4.8036627318, 5.4044338329, 6.0517560386, 6.8415021436, &
        7.8080659794, 8.9257124242, 10.1640694041, 11.5476843252, 13.1494718506, &
        14.9965430946, 16.972215146, 18.903468102, 20.9125176371, 23.3776351255  &
        /) ,(/NTHICK,NFTAB/))
 
ELSE
! May be changed, but according to T. Williams, wim1 is okay from 0.25 only
       SIS1HTABLE = (/ 0.25, 0.35, 0.6, 0.85, 1.1, 1.35, 1.6, 1.85, 2.1,  2.35, &
                      2.6, 2.85, 3.1, 3.35, 3.6, 3.85, 4.1, 4.35, 4.6, 4.85 /)
       SIS1ALPHATABLE = reshape((/  &
       3.80373e-06 ,  6.02822e-06 , 1.12121e-05 , 2.24588e-05 , &
3.05165e-05 ,  3.85142e-05 , 5.29712e-05 , 7.55453e-05 , &
1.05531e-04 ,  1.38952e-04 , 1.73378e-04 , 2.09421e-04 , &
2.47299e-04 ,  2.87669e-04 , 3.34304e-04 , 3.91094e-04 , &
4.57710e-04 ,  5.31901e-04 , 6.14490e-04 , 7.08299e-04 , &
4.76980e-06 ,  6.95433e-06 , 1.54179e-05 , 2.46138e-05 , &
3.79955e-05 ,  5.96725e-05 , 8.81668e-05 , 1.20881e-04 , &
1.57364e-04 ,  1.99666e-04 , 2.50723e-04 , 3.13346e-04 , &
3.89907e-04 ,  4.81606e-04 , 5.88164e-04 , 7.08707e-04 , &
8.42536e-04 ,  9.88622e-04 , 1.14498e-03 , 1.30931e-03 , &
5.84556e-06 ,  8.58738e-06 , 2.00947e-05 , 3.18410e-05 , &
5.67385e-05 ,  9.37456e-05 , 1.37304e-04 , 1.87832e-04 , &
2.49218e-04 ,  3.28900e-04 , 4.31988e-04 , 5.59364e-04 , &
7.11302e-04 ,  8.86736e-04 , 1.08074e-03 , 1.28808e-03 , &
1.50809e-03 ,  1.74230e-03 , 1.98930e-03 , 2.24607e-03 , &
6.90305e-06 ,  1.12155e-05 , 2.59514e-05 , 5.15063e-05 , &
9.18902e-05 ,  1.45703e-04 , 2.18304e-04 , 3.18476e-04 , &
4.51992e-04 ,  6.21037e-04 , 8.24428e-04 , 1.05881e-03 , &
1.32112e-03 ,  1.60970e-03 , 1.92396e-03 , 2.26484e-03 , &
2.63577e-03 ,  3.04238e-03 , 3.49116e-03 , 3.98911e-03 , &
8.46004e-06 ,  1.46791e-05 , 3.69837e-05 , 8.32870e-05 , &
1.49151e-04 ,  2.46175e-04 , 3.92165e-04 , 5.95871e-04 , &
8.57938e-04 ,  1.17280e-03 , 1.53668e-03 , 1.95109e-03 , &
2.42048e-03 ,  2.95341e-03 , 3.56468e-03 , 4.27150e-03 , &
5.08843e-03 ,  6.02927e-03 , 7.11153e-03 , 8.35459e-03 , &
1.11513e-05 ,  1.95090e-05 , 5.78215e-05 , 1.32267e-04 , &
2.57932e-04 ,  4.61051e-04 , 7.55044e-04 , 1.14193e-03 , &
1.62478e-03 ,  2.21208e-03 , 2.92101e-03 , 3.77684e-03 , &
4.80845e-03 ,  6.04635e-03 , 7.52256e-03 , 9.26818e-03 , &
1.13101e-02 ,  1.36712e-02 , 1.63732e-02 , 1.94360e-02 , &
1.52571e-05 ,  2.78128e-05 , 9.43665e-05 , 2.30636e-04 , &
4.90057e-04 ,  9.02105e-04 , 1.48012e-03 , 2.24822e-03 , &
3.24754e-03 ,  4.53444e-03 , 6.16979e-03 , 8.21171e-03 , &
1.07145e-02 ,  1.37253e-02 , 1.72805e-02 , 2.14091e-02 , &
2.61367e-02 ,  3.14840e-02 , 3.74632e-02 , 4.40798e-02 , &
2.13731e-05 ,  4.29771e-05 , 1.64050e-04 , 4.49510e-04 , &
9.84576e-04 ,  1.82921e-03 , 3.06425e-03 , 4.79802e-03 , &
7.14520e-03 ,  1.02139e-02 , 1.40940e-02 , 1.88529e-02 , &
2.45396e-02 ,  3.11875e-02 , 3.88141e-02 , 4.74255e-02 , &
5.70233e-02 ,  6.76049e-02 , 7.91611e-02 , 9.16776e-02 , &
3.17688e-05 ,  7.05749e-05 , 3.13449e-04 , 9.26652e-04 , &
2.08059e-03 ,  3.98338e-03 , 6.87070e-03 , 1.09528e-02 , &
1.63868e-02 ,  2.32726e-02 , 3.16649e-02 , 4.15851e-02 , &
5.30293e-02 ,  6.59789e-02 , 8.04085e-02 , 9.62863e-02 , &
1.13575e-01 ,  1.32236e-01 , 1.52233e-01 , 1.73533e-01 , &
5.14019e-05 ,  1.24975e-04 , 6.49063e-04 , 2.01553e-03 , &
4.70028e-03 ,  9.17913e-03 , 1.58010e-02 , 2.47576e-02 , &
3.61186e-02 ,  4.98708e-02 , 6.59568e-02 , 8.42994e-02 , &
1.04811e-01 ,  1.27404e-01 , 1.51996e-01 , 1.78511e-01 , &
2.06873e-01 ,  2.37011e-01 , 2.68862e-01 , 3.02365e-01 , &
9.05394e-05 ,  2.44433e-04 , 1.43884e-03 , 4.66981e-03 , &
1.09558e-02 ,  2.08947e-02 , 3.46681e-02 , 5.22109e-02 , &
7.33489e-02 ,  9.78738e-02 , 1.25579e-01 , 1.56274e-01 , &
1.89790e-01 ,  2.25973e-01 , 2.64687e-01 , 3.05808e-01 , &
3.49223e-01 ,  3.94830e-01 , 4.42530e-01 , 4.92227e-01 , &
1.75074e-04 ,  5.25080e-04 , 3.39573e-03 , 1.10036e-02 , &
2.46312e-02 ,  4.44096e-02 , 6.99511e-02 , 1.00742e-01 , &
1.36302e-01 ,  1.76219e-01 , 2.20147e-01 , 2.67802e-01 , &
3.18937e-01 ,  3.73339e-01 , 4.30813e-01 , 4.91183e-01 , &
5.54282e-01 ,  6.19954e-01 , 6.88048e-01 , 7.58418e-01 , &
3.74908e-04 ,  1.22443e-03 , 8.17252e-03 , 2.46994e-02 , &
5.11207e-02 ,  8.63422e-02 , 1.29173e-01 , 1.78641e-01 , &
2.34006e-01 ,  2.94691e-01 , 3.60235e-01 , 4.30258e-01 , &
5.04426e-01 ,  5.82438e-01 , 6.64019e-01 , 7.48914e-01 , &
8.36872e-01 ,  9.27649e-01 , 1.02102e+00 , 1.11676e+00 , &
8.82800e-04 ,  3.03043e-03 , 1.88248e-02 , 5.09125e-02 , &
9.68861e-02 ,  1.54064e-01 , 2.20524e-01 , 2.94949e-01 , &
3.76394e-01 ,  4.64138e-01 , 5.57588e-01 , 6.56231e-01 , &
7.59605e-01 ,  8.67281e-01 , 9.78849e-01 , 1.09393e+00 , &
1.21214e+00 ,  1.33315e+00 , 1.45662e+00 , 1.58224e+00 , &
2.22521e-03 ,  7.51676e-03 , 3.98357e-02 , 9.57376e-02 , &
1.68983e-01 ,  2.55527e-01 , 3.52939e-01 , 4.59644e-01 , &
5.74490e-01 ,  6.96566e-01 , 8.25086e-01 , 9.59331e-01 , &
1.09864e+00 ,  1.24238e+00 , 1.38997e+00 , 1.54085e+00 , &
1.69451e+00 ,  1.85049e+00 , 2.00833e+00 , 2.16765e+00 , &
5.66833e-03 ,  1.75464e-02 , 7.66453e-02 , 1.65727e-01 , &
2.74951e-01 ,  3.99380e-01 , 5.36249e-01 , 6.83721e-01 , &
8.40370e-01 ,  1.00497e+00 , 1.17642e+00 , 1.35368e+00 , &
1.53580e+00 ,  1.72190e+00 , 1.91118e+00 , 2.10291e+00 , &
2.29645e+00 ,  2.49123e+00 , 2.68675e+00 , 2.88259e+00 , &
1.36539e-02 ,  3.72456e-02 , 1.35140e-01 , 2.67780e-01 , &
4.22901e-01 ,  5.95111e-01 , 7.81333e-01 , 9.79307e-01 , &
1.18712e+00 ,  1.40304e+00 , 1.62548e+00 , 1.85300e+00 , &
2.08433e+00 ,  2.31832e+00 , 2.55400e+00 , 2.79055e+00 , &
3.02727e+00 ,  3.26361e+00 , 3.49910e+00 , 3.73337e+00 , &
2.98722e-02 ,  7.15314e-02 , 2.21423e-01 , 4.09271e-01 , &
6.21810e-01 ,  8.53369e-01 , 1.10028e+00 , 1.35949e+00 , &
1.62829e+00 ,  1.90425e+00 , 2.18522e+00 , 2.46937e+00 , &
2.75516e+00 ,  3.04136e+00 , 3.32695e+00 , 3.61118e+00 , &
3.89345e+00 ,  4.17331e+00 , 4.45047e+00 , 4.72471e+00 , &
5.88578e-02 ,  1.25612e-01 , 3.41821e-01 , 5.98292e-01 , &
8.81793e-01 ,  1.18598e+00 , 1.50596e+00 , 1.83748e+00 , &
2.17675e+00 ,  2.52062e+00 , 2.86649e+00 , 3.21233e+00 , &
3.55662e+00 ,  3.89823e+00 , 4.23636e+00 , 4.57049e+00 , &
4.90029e+00 ,  5.22559e+00 , 5.54632e+00 , 5.86250e+00 , &
1.05468e-01 ,  2.04724e-01 , 5.03158e-01 , 8.44026e-01 , &
1.21416e+00 ,  1.60557e+00 , 2.01139e+00 , 2.42573e+00 , &
2.84383e+00 ,  3.26205e+00 , 3.67776e+00 , 4.08914e+00 , &
4.49503e+00 ,  4.89472e+00 , 5.28788e+00 , 5.67440e+00 , &
6.05432e+00 ,  6.42782e+00 , 6.79514e+00 , 7.15656e+00 , &
1.74545e-01 ,  3.14137e-01 , 7.13063e-01 , 1.15691e+00 , &
1.63126e+00 ,  2.12499e+00 , 2.62866e+00 , 3.13493e+00 , &
3.63860e+00 ,  4.13623e+00 , 4.62576e+00 , 5.10610e+00 , &
5.57680e+00 ,  6.03787e+00 , 6.48955e+00 , 6.93227e+00 , &
7.36652e+00 ,  7.79283e+00 , 8.21174e+00 , 8.62378e+00 , &
2.70834e-01 ,  4.59350e-01 , 9.80330e-01 , 1.54856e+00 , &
2.14570e+00 ,  2.75617e+00 , 3.36813e+00 , 3.97376e+00 , &
4.56847e+00 ,  5.14997e+00 , 5.71743e+00 , 6.27093e+00 , &
6.81101e+00 ,  7.33850e+00 , 7.85432e+00 , 8.35942e+00 , &
8.85471e+00 ,  9.34106e+00 , 9.81925e+00 , 1.02900e+01 , &
3.99120e-01 ,  6.46435e-01 , 1.31507e+00 , 2.03136e+00 , &
2.76957e+00 ,  3.50946e+00 , 4.23832e+00 , 4.94958e+00 , &
5.64069e+00 ,  6.31137e+00 , 6.96253e+00 , 7.59564e+00 , &
8.21234e+00 ,  8.81426e+00 , 9.40295e+00 , 9.97983e+00 , &
1.05461e+01 ,  1.11030e+01 , 1.16515e+01 , 1.21924e+01 , &
5.64538e-01 ,  8.82405e-01 , 1.72870e+00 , 2.61770e+00 , &
3.51364e+00 ,  4.39362e+00 , 5.24695e+00 , 6.07052e+00 , &
6.86499e+00 ,  7.63269e+00 , 8.37647e+00 , 9.09919e+00 , &
9.80351e+00 ,  1.04918e+01 , 1.11660e+01 , 1.18280e+01 , &
1.24792e+01 ,  1.31209e+01 , 1.37543e+01 , 1.43803e+01 , &
7.72889e-01 ,  1.17552e+00 , 2.23335e+00 , 3.31894e+00 , &
4.38704e+00 ,  5.41672e+00 , 6.40295e+00 , 7.34788e+00 , &
8.25618e+00 ,  9.13289e+00 , 9.98272e+00 , 1.08097e+01 , &
1.16174e+01 ,  1.24086e+01 , 1.31857e+01 , 1.39506e+01 , &
1.47049e+01 ,  1.54501e+01 , 1.61872e+01 , 1.69173e+01   &
        /) ,(/NTHICK,NFTAB/))
        END IF
      DO I=1,NK
        DO J=1,NTHICK
          IF (SIG(I)*TPIINV.LT.SIS1FTABLE(1)) THEN
            SIS2ALPHAS(J,I) = SIS1ALPHATABLE(J,1)
          ELSE IF (SIG(I)*TPIINV.GT. SIS1FTABLE(NFTAB)) THEN
            SIS2ALPHAS(J,I) = SIS1ALPHATABLE(J,NFTAB)
          ELSE
            IND = 1
            DO K = 1, NFTAB-1
              IF (SIS1FTABLE(K).LT.SIG(I)*TPIINV) IND = K
              END DO
            X=(SIG(I)*TPIINV-SIS1FTABLE(IND))/(SIS1FTABLE(IND+1)-SIS1FTABLE(IND))
            SIS2ALPHAS(J,I)=SIS1ALPHATABLE(J,IND)*(1-X)+SIS1ALPHATABLE(J,IND+1)*X
            END IF
! WRITE(998,*) I, J, SIG(I)*TPIINV,SIS1FTABLE(NFTAB), X, IND, SIS2ALPHAS(J,I),SIS1ALPHATABLE(J,NFTAB)
          END DO
        END DO
!
 
!
      SIS1ALPHATABLE2 = reshape((/  &
0.000001693306, 0.000001694535, 0.000001709985, 0.000001718371, 0.000001718689, &
0.000001715291, 0.000001711765, 0.000001709562, 0.000001708573, 0.000001708047, &
0.000001707216, 0.000001705508, 0.000001702489, 0.000001697757, 0.000001690926, &
0.000001681715, 0.000001670039, 0.000001655970, 0.000001639540, 0.000001620666, &
0.000002345944, 0.000002373942, 0.000002352788, 0.000002336058, 0.000002340157, &
0.000002361578, 0.000002389818, 0.000002414910, 0.000002430787, 0.000002435510, &
0.000002429671, 0.000002414451, 0.000002390516, 0.000002358094, 0.000002317788, &
0.000002271221, 0.000002220674, 0.000002167641, 0.000002111561, 0.000002051905, &
0.000003223057, 0.000003267049, 0.000003224977, 0.000003204685, 0.000003222416, &
0.000003261346, 0.000003299552, 0.000003321243, 0.000003319410, 0.000003293875, &
0.000003247777, 0.000003184750, 0.000003107742, 0.000003019246, 0.000002922032, &
0.000002819435, 0.000002714739, 0.000002609966, 0.000002505271, 0.000002400993, &
0.000004390462, 0.000004433012, 0.000004403102, 0.000004405652, 0.000004433604, &
0.000004455373, 0.000004447247, 0.000004400212, 0.000004316137, 0.000004201667, &
0.000004064075, 0.000003909816, 0.000003744650, 0.000003573915, 0.000003402281, &
0.000003233187, 0.000003068613, 0.000002909686, 0.000002757781, 0.000002614371, &
0.000005932907, 0.000005957402, 0.000005982223, 0.000006017129, 0.000006011597, &
0.000005931834, 0.000005775734, 0.000005561133, 0.000005309860, 0.000005038831, &
0.000004758982, 0.000004477929, 0.000004202314, 0.000003938073, 0.000003689194, &
0.000003456677, 0.000003239230, 0.000003035853, 0.000002848219, 0.000002677864, &
0.000007962387, 0.000007965182, 0.000008068497, 0.000008081519, 0.000007918056, &
0.000007586175, 0.000007143591, 0.000006652694, 0.000006157274, 0.000005680186, &
0.000005230854, 0.000004813199, 0.000004429817, 0.000004082375, 0.000003770294, &
0.000003490093, 0.000003236746, 0.000003006720, 0.000002800108, 0.000002616762, &
0.000010630745, 0.000010631503, 0.000010754695, 0.000010555056, 0.000009993969, &
0.000009206201, 0.000008341755, 0.000007500204, 0.000006728370, 0.000006039825, &
0.000005432987, 0.000004901489, 0.000004438118, 0.000004035401, 0.000003685333, &
0.000003379621, 0.000003110711, 0.000002873062, 0.000002663457, 0.000002478682, &
0.000014147503, 0.000014179502, 0.000014068022, 0.000013255509, 0.000011955952, &
0.000010517681, 0.000009161184, 0.000007970672, 0.000006956434, 0.000006102184, &
0.000005386202, 0.000004787198, 0.000004285040, 0.000003861250, 0.000003500220, &
0.000003190247, 0.000002923212, 0.000002692594, 0.000002491216, 0.000002312444, &
0.000018803447, 0.000018843681, 0.000017888469, 0.000015846997, 0.000013469781, &
0.000011291688, 0.000009484334, 0.000008030183, 0.000006862685, 0.000005921675, &
0.000005161810, 0.000004546889, 0.000004045234, 0.000003629478, 0.000003278778, &
0.000002980118, 0.000002726294, 0.000002510734, 0.000002323048, 0.000002154027, &
0.000024999293, 0.000024767047, 0.000021866213, 0.000017907275, 0.000014286707, &
0.000011453674, 0.000009339677, 0.000007754652, 0.000006540404, 0.000005593636, &
0.000004848620, 0.000004258144, 0.000003783667, 0.000003393700, 0.000003065741, &
0.000002787091, 0.000002551635, 0.000002353114, 0.000002179924, 0.000002022176, &
0.000033276027, 0.000031812088, 0.000025413911, 0.000019081783, 0.000014360014, &
0.000011113013, 0.000008878307, 0.000007287536, 0.000006111165, 0.000005216269, &
0.000004523259, 0.000003978625, 0.000003542278, 0.000003183888, 0.000002882748, &
0.000002627060, 0.000002410531, 0.000002226889, 0.000002066156, 0.000001920400, &
0.000044336466, 0.000039340417, 0.000027862847, 0.000019248910, 0.000013854771, &
0.000010491793, 0.000008293280, 0.000006775936, 0.000005680752, 0.000004862822, &
0.000004234499, 0.000003739067, 0.000003338786, 0.000003008719, 0.000002732485, &
0.000002498817, 0.000002299071, 0.000002126144, 0.000001974716, 0.000001841606, &
0.000059030780, 0.000046157443, 0.000028771759, 0.000018580063, 0.000013047269, &
0.000009812168, 0.000007739586, 0.000006322672, 0.000005313968, 0.000004570787, &
0.000004001034, 0.000003546046, 0.000003172357, 0.000002862735, 0.000002606289, &
0.000002391441, 0.000002205202, 0.000002038837, 0.000001893870, 0.000001773775, &
0.000078242681, 0.000050872198, 0.000028192943, 0.000017445281, 0.000012194039, &
0.000009215299, 0.000007293336, 0.000005968851, 0.000005030058, 0.000004344187, &
0.000003817214, 0.000003389423, 0.000003031921, 0.000002735020, 0.000002492733, &
0.000002292508, 0.000002116550, 0.000001954358, 0.000001814339, 0.000001706341, &
0.000102553944, 0.000052628108, 0.000026654291, 0.000016235235, 0.000011448810, &
0.000008738110, 0.000006951009, 0.000005700488, 0.000004811678, 0.000004163305, &
0.000003662892, 0.000003251855, 0.000002905015, 0.000002617308, 0.000002385271, &
0.000002195475, 0.000002027393, 0.000001869747, 0.000001734845, 0.000001635748, &
0.000131528803, 0.000051648919, 0.000024847696, 0.000015214276, 0.000010847813, &
0.000008336757, 0.000006656766, 0.000005471182, 0.000004622914, 0.000003999733, &
0.000003516266, 0.000003119051, 0.000002785020, 0.000002508708, 0.000002285704, &
0.000002102851, 0.000001941273, 0.000001790864, 0.000001662390, 0.000001566768, &
0.000162640998, 0.000049040647, 0.000023269162, 0.000014460197, 0.000010345625, &
0.000007943173, 0.000006352063, 0.000005238097, 0.000004434583, 0.000003834908, &
0.000003367217, 0.000002987629, 0.000002673539, 0.000002414062, 0.000002200682, &
0.000002022244, 0.000001866432, 0.000001726382, 0.000001605605, 0.000001507441, &
0.000190593962, 0.000046018523, 0.000022036249, 0.000013884361, 0.000009879926, &
0.000007534440, 0.000006029990, 0.000004998037, 0.000004246561, 0.000003673278, &
0.000003223583, 0.000002865466, 0.000002575963, 0.000002336404, 0.000002132787, &
0.000001957123, 0.000001806342, 0.000001677844, 0.000001564855, 0.000001460377, &
0.000208777708, 0.000043254960, 0.000020953050, 0.000013327645, 0.000009433705, &
0.000007161759, 0.000005746018, 0.000004790209, 0.000004085598, 0.000003536085, &
0.000003102390, 0.000002762488, 0.000002493229, 0.000002269788, 0.000002073862, &
0.000001900259, 0.000001753553, 0.000001634519, 0.000001527844, 0.000001418102, &
0.000213643676, 0.000040801760, 0.000019820045, 0.000012719494, 0.000009042103, &
0.000006888517, 0.000005551025, 0.000004645394, 0.000003968666, 0.000003433279, &
0.000003008733, 0.000002677964, 0.000002418184, 0.000002202481, 0.000002011120, &
0.000001840100, 0.000001696853, 0.000001583290, 0.000001480832, 0.000001370962, &
0.000208202335, 0.000038538552, 0.000018774435, 0.000012168047, 0.000008723021, &
0.000006681551, 0.000005399851, 0.000004523630, 0.000003863808, 0.000003339344, &
0.000002922742, 0.000002598268, 0.000002343776, 0.000002132781, 0.000001945837, &
0.000001779030, 0.000001639702, 0.000001529688, 0.000001430567, 0.000001324016, &
0.000198860865, 0.000036596393, 0.000018094381, 0.000011782298, 0.000008423914, &
0.000006432368, 0.000005194822, 0.000004356699, 0.000003726655, 0.000003224518, &
0.000002825485, 0.000002515982, 0.000002274472, 0.000002073961, 0.000001894595, &
0.000001733161, 0.000001598541, 0.000001493403, 0.000001398055, 0.000001292778, &
0.000188086831, 0.000035005432, 0.000017449034, 0.000011378985, 0.000008128438, &
0.000006206255, 0.000005018705, 0.000004216819, 0.000003612753, 0.000003129463, &
0.000002744757, 0.000002446723, 0.000002214548, 0.000002021421, 0.000001847651, &
0.000001690537, 0.000001559668, 0.000001458051, 0.000001365539, 0.000001262003, &
0.000176803511, 0.000033479878, 0.000016713766, 0.000010962586, 0.000007891418, &
0.000006057305, 0.000004907363, 0.000004122061, 0.000003528516, 0.000003054204, &
0.000002676573, 0.000002382977, 0.000002153342, 0.000001962633, 0.000001792438, &
0.000001639642, 0.000001512208, 0.000001412404, 0.000001322051, 0.000001223019, &
0.000157501278, 0.000031061181, 0.000013781133, 0.000009601049, 0.000007947466, &
0.000006760411, 0.000005655833, 0.000004680347, 0.000003892156, 0.000003277589, &
0.000002790517, 0.000002392986, 0.000002067558, 0.000001808761, 0.000001609460, &
0.000001454171, 0.000001323348, 0.000001206327, 0.000001111024, 0.000001045460, &
0.000007782240, 0.000002755736, 0.000000005024, 0.000000101927, 0.000088216814, &
0.009792270789, 0.030135654570, 0.009293514348, 0.001549599316, 0.000298014001, &
0.000055965309, 0.000007063938, 0.000000637549, 0.000000072256, 0.000000017854, &
0.000000008901, 0.000000004238, 0.000000001169, 0.000000000492, 0.000000001113 /), &
(/NTHICK,NFTAB/))
 
      DO I=1,NK
        DO J=1,NTHICK
          IF (SIG(I)*TPIINV.LT. SIS1FTABLE(1)) THEN
            SIS2ALPHA2(J,I) = SIS1ALPHATABLE2(J,1)
          ELSE IF (SIG(I)*TPIINV.GT. SIS1FTABLE(NFTAB)) THEN
            SIS2ALPHA2(J,I) = SIS1ALPHATABLE2(J,NFTAB)
          ELSE
            IND = 1
            DO K = 1, NFTAB-1
              IF (SIS1FTABLE(K).LT.SIG(I)*TPIINV) IND = K
              END DO
            X=(SIG(I)*TPIINV-SIS1FTABLE(IND))/(SIS1FTABLE(IND+1)-SIS1FTABLE(IND))
            SIS2ALPHA2(J,I)=SIS1ALPHATABLE2(J,IND)*(1-X)+SIS1ALPHATABLE2(J,IND+1)*X
            END IF
          END DO
        END DO
!
! -------------------------------------------------------------------- /
! 2.  Fills array of ICEDMAX to ICEDAVE
!
      DO I=1,NICED
        ICEDAVETAB(I) = W3FSD_DAVE(IS2PARS(9),REAL(I),IS2PARS(8))
        ENDDO
!
! -------------------------------------------------------------------- /
! 2.  Defines and diagonalizes the scattering matrix
!
      ALLOCATE(IS2SCATMAT(NTH,NTH))
      ALLOCATE(IS2EIGVEC(NTH,NTH))
      ALLOCATE(IS2EIGVAL(NTH))
!
      DO I=1,NTH
        DO J=1,NTH
! This is for isotropic back-scatter
!         IS2SCATMAT(I,J)=-1./DBLE(NTH)
! Other example that looks like figure 12 in Masson & LeBlond
          IS2SCATMAT(I,J)=-1./DBLE(NTH)*(2.*EC2(ABS(I-J)+1)**2+0.8*ECOS(ABS(I-J)+1)**3)
          IF (ECOS(ABS(I-J)+1).LT.0.001) IS2SCATMAT(I,J)=IS2SCATMAT(I,J)-1./DBLE(NTH)*0.8*ES2(ABS(I-J)+1)
          END DO
!WRITE(997,'(36G16.8)') IS2SCATMAT(I,:)
! Now removes sum from diagonal to enforce energy conservation ...
          IS2SCATMAT(I,I)=IS2SCATMAT(I,I)-SUM(IS2SCATMAT(I,1:NTH))
        END DO
      CALL DIAGONALIZE(IS2SCATMAT,IS2EIGVAL,IS2EIGVEC,nrot)
      DO I=1,NTH
        IS2EIGVAL(I)=MAX(0.d0,IS2EIGVAL(I))
!WRITE(994,'(36G16.8)') I,IS2EIGVAL(I)
!WRITE(995,'(36G16.8)') IS2EIGVEC(I,:)
!WRITE(996,'(36G16.8)') IS2SCATMAT(I,:)
        END DO
!CLOSE(994)
!CLOSE(995)
!CLOSE(996)
!CLOSE(997)
 
      RETURN
!/
!/ End of INSIS2 ----------------------------------------------------- /
!/
      END SUBROUTINE INSIS2
!/ ------------------------------------------------------------------- /
      SUBROUTINE W3SIS2 (A, DEPTH, CICE, ICEH, ICEF, ICEDMAX, IX, IY, &
                         S, D, DISSIP, WN, CG, WN_R, CG_ICE, R)
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  | P. Nicot & F. Ardhuin & G. Boutin |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         04-May-2016 |
!/                  +-----------------------------------+
!/
!/    16-Mar-2014 : Origination.                        ( version 4.18 )
!/    19-Sep-2014 : Correcting group speed factor       ( version 5.03 )
!/    20-Sep-2014 : Adding back-scattered energy        ( version 5.03 )
!/    27-Aug-2015 : Add breaking criterion, WIM1d       ( version 5.05 )
!/                  (ref. Williams, 2012)
!/    02-Nov-2015 : Integration of strain over bandwidth( version 5.05 )
!/    13-Jan-2016 : Changed initialization of ICEDMAX   ( version 5.10 )
!/    06-Feb-2016 : Added IICEHMIN and creep dissipation( version 5.10 )
!/    10-Mar-2016 : Added depth and call to Liu disp.   ( version 5.10 )
!/    02-May-2016 : Call to Liu disp moved to w3srce    ( version 5.10 )
!
!  1. Purpose :
!     Wave scattering in the MIZ, adapted from Dumont et al.
!
!/ ------------------------------------------------------------------- /
!
!  2. Method :
!     This scattering routine allows the estimation of the maximum floe
!     size and an estimate of the creep-induced dissipation.
!     For the scattering, it is based on the normal incidence results of
!     Kohout and Meylan which are provided in a table.
!
!  3. Parameters :
!
!     Parameter list
!     ----------------------------------------------------------------
!       A       R.A.  I     Action density spectrum (1-D)
!       DEPTH   Real  I     Water depth
!       CICE    Real  I     Sea ice concentration
!       ICEH    Real  I     ice thickness
!       ICEF    Real  I/O   Maximum floe size (updated)
!       ICEDMAX Real  I     Maximum floe size
!       IX,IY   Int   I     Not used
!       S       R.A.  O     Source term (1-D version)
!       D       R.A.  O     Diagonal part of scattering (1-D version)
!       DISSIP  R.A.  O     Diagonal dissipation term (1-D version)
!       WN      R.A.  I     Wave number
!       CG      R.A.  I     Group speed
!       WN_R    R.A.  I     Wave number in ice
!       CG_ICE  R.A.  I     Group speed in ice
!       R       R.A.  O     Ratio of energy to wave energy without ice
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!  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 concentration is zero, no calculations are made.
!
!  8. Structure :
!
!     See source code.
!
!  9. Switches :
!
!     !/S  Enable subroutine tracing.
!     !/T   Enable general test output.
!           2-D print plot of source term.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE W3ODATMD, ONLY: NDSE
      USE CONSTANTS, ONLY: TPIINV, PI, TPI, GRAV, DWAT
      USE W3SERVMD, ONLY: EXTCDE
      USE W3GDATMD, ONLY: NK, NTH, NSPEC, SIG, SIG2, DDEN, IS2PARS, XFR, &
                          IICEHMIN,IICESMOOTH
      USE W3DISPMD
!
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
      REAL, INTENT(IN)        :: A(NSPEC), DEPTH, CICE, ICEH, ICEDMAX
      INTEGER, INTENT(IN)     :: IX, IY
      REAL, INTENT(OUT)       :: S(NSPEC), D(NSPEC), DISSIP(NSPEC), R(NK)
      REAL, INTENT(INOUT)     :: ICEF
      REAL, INTENT(IN)        :: WN(NK), CG(NK), WN_R(NK), CG_ICE(NK)
 
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      INTEGER                 :: IK, IKP1, IKM1, ITH, ITH2, IS, IS2, IND1, IND2
      REAL                    :: W
      INTEGER                 :: IKBREAK, NSUM
      LOGICAL                 :: BRK_CRIT_W(NK), BRK_CRIT
      REAL                    :: ALPHA, STRAIN_C, WAMP(NK), D_FLEX_FAIL,             &
                                 SMOOTHD, TAU_D, S_D(NK), ALPHA_D, DELTA_D,B_COLE,    &
                                 DMAX, S_ATT, FACTOR, BETA
      REAL                    :: ICEDAVE(NK), CURVTOSTRAIN, CREEPFAC, MP2, B, ICEF_CREEP
      REAL                    :: SUMALLDIR, SUMA, SUME, CURVSPEC(NK), ESPEC(NK),STRAIN
      REAL, PARAMETER         :: YOUNG = 5.49E+9   ! Young modulus
      REAL, PARAMETER         :: POISSON = 0.3     ! Poisson Ratio
      REAL                    :: SIGMA_C
      REAL, PARAMETER         :: DENS = 1025.0     ! ice density
      REAL                    :: GAMMA_TOY
      REAL, DIMENSION(NK)        :: WN_I, WN_RP, WSQ, WLG, WLG_I, CG_I,   &
                                    CURV, CGRATIO, CG_EFF, DUMMY, ALPHA_DISP
!/
!/ ------------------------------------------------------------------- /
!/
!
! 0.  Initializations ------------------------------------------------ *
!
      S = 0.
      D = 0.
      DISSIP = 0.
      DUMMY = WN
      WN_I = 0.
      WN_RP = 0.
      CG_I = 0.
      CG_EFF = 0.
      SIGMA_C = IS2PARS(19)
      MP2=(1-POISSON**2)
      GAMMA_TOY =  2 + log(0.9)/log(2.)
      ! Variables from Cole et al. 1995
      ALPHA_D=0.54
      !B_COLE=1.205E-9 * EXP(0.55*1.60218E-19/(1.38064852e-23*263.15))
      B_COLE=1.205E-9 * EXP(IS2PARS(24)*1.60218E-19/(1.38064852e-23*268.15))
      TAU_D=B_COLE/0.07
      S_D=LOG(SIG(1:NK)*TAU_D)
      DELTA_D=IS2PARS(21)
!
      IF (IS2PARS(9).GT.0) ICEDMIN = IS2PARS(9)
      IF (IS2PARS(12).GT.0) THEN
        B=IS2PARS(12)
! 2 is the ratio of <H/2^4> to (Hrms/2)^4 for a Rayleigh distribution
! 0.375 is the average of cos^4
! 0.4 is 2/5
        CREEPFAC = -2*(0.25/(IS2PARS(15)+2))*0.375*B*ICEH**(IS2PARS(15)+2)  &
                        *(YOUNG/(2*B*MP2))**(IS2PARS(15)+1)/(DWAT*GRAV)
      ELSE
        CREEPFAC=1.
        ENDIF
 
      WLG = TPI/WN
      ICEF_CREEP=ICEF
      DMAX = ICEF
      BRK_CRIT = .FALSE.
      NSUM = NINT(0.3/(XFR-1.))
!
      STRAIN_C = SIGMA_C*MP2/YOUNG
 
! Minimum floe size that can break
      D_FLEX_FAIL =  0.5* ( (PI**4*YOUNG*ICEH**3)/( 48*DENS*GRAV*MP2 ) )**.25
! Estimates mean floe diameter from max floe diameter
      IF (DMAX.GT.NICED) THEN
        ICEDAVE=DMAX
      ELSE
        IF (IS2PARS(9).GT.0) THEN
          ICEDAVE= ICEDAVETAB(MAX(1,NINT(DMAX)))
        ELSE
          DO IK=1,NK
            IF (IS2PARS(14)*(TPI/WN_R(IK)).LT.DMAX) THEN
              ICEDAVE(IK)=W3FSD_DAVE( (TPI/WN_R(IK))*IS2PARS(14),DMAX,IS2PARS(8))
            ELSE
              ICEDAVE(IK)=DMAX
              ENDIF
            END DO
          ENDIF
        ENDIF
!
      IF (CICE .GT. 0) THEN
!
!  1. Calculate wavelength, Robinson-Palmer dispersion relation
!     (should be tabulated )
!     Warning: it only applies to unbroken ice.
!
        IF (IS2PARS(6).GT.0.5) THEN
          CALL W3RPWNICE(ICEH,WN_I,WN_RP,CG_I)
        ELSE
          WN_I  = WN_R
          CG_I  = CG_ICE
          WN_RP=WN*0.
          END IF
!
        WLG_I = TPI/WN_I  ! Ice wavelength
        WSQ = (WLG/WLG_I)
!
        IF (IS2PARS(16).GT.0.5) THEN
          CGRATIO(:)=CG_I(:)/CG(:)
        ELSE
          CGRATIO(:)=1
          END IF
!
!  2. gets reflection coefficient from table
!
        IND1 = 1+FLOOR((ICEH-THICK1)/DTHICK)
        IND2 = IND1+1
! defines weight for interpolation of ice thickness
        W = (ICEH-THICK1)/DTHICK - (IND1-1)
        IF (IND1.LT.1) THEN
          IND1 = 1
          W = 0
        ELSE IF (IND2.GT.NTHICK) THEN
          IND2 = NTHICK
          IND1 = NTHICK
          W = 0
          END IF
!
        DO IK = 1,NK
!
! Spatial decay scale taken from table. This corresponds the the values shown
! in Dumont et al. (JGR 2011, fig 3).
!
          CG_EFF(IK) = CICE*CG_I(IK) + (1-CICE)*CG(IK)
 
! Note by FA: dissipation should be done by ICx not by ISx
! the RP damping is thus defined by an optional IS2PARS(7), which is 0 by default
          ALPHA = -1.*(  (SIS2ALPHAS(IND1,IK)*(1-W)+SIS2ALPHAS(IND2,IK)*W)/ICEDAVE(IK) &
                        +2.*IS2PARS(7)*WN_RP(IK) )*CICE
!
! Additional scattering for pack ice defined by IS2PARS(4:5)  (see Squire et al. GRL 2009)
!
          ALPHA = IS2PARS(1) * ALPHA -2.*IS2PARS(4)*EXP(-1.*IS2PARS(5)/SIG(IK))
          IF (IS2PARS(11).GT.0) THEN
            IF (CICE.LT.0.2) ALPHA = 0.
            IF (CICE.GT.0.8) ALPHA = 0.
            IF (CICE.GE.0.2.AND.CICE.LE.0.8) ALPHA = ALPHA*(CICE-0.2)*(0.8-CICE)
            END IF
!
! time decay
!
          BETA = ALPHA * CICE * CG_EFF(IK)
!
!  3. attenuation due to scattering for all spectral components
!     with added backscattering for energy conservation ( if IS2PARS(2).EQ.1)
!
          SUMALLDIR= 0.
          SUMA = 0.
          CURVTOSTRAIN =  (0.25*MAX(ICEH,IICEHMIN)**2)
          DO ITH = 1,NTH
            IS  = ITH+(IK-1)*NTH
            D(IS) = BETA
            S(IS) = BETA * A(IS)
            SUMALLDIR = SUMALLDIR + S(IS)
            SUMA = SUMA + A(IS)
            END DO ! loop over directions
!
! R is the ratio of energy (including bending of ice) to wave energy without ice
! Wadhams 1973 eq. 34, warning, his ice thickness is 2*h
! Warning : R uses DMAX=ICEF, even if IS2DUPDATE=F
!
          IF (IICESMOOTH) THEN
            IF (IS2PARS(14)*WLG_I(IK).LT.DMAX) THEN
              SMOOTHD=TANH((DMAX-IS2PARS(14)*WLG_I(IK))/(DMAX*IS2PARS(13)))
            ELSE
              SMOOTHD=0.
              END IF
          ELSE
            SMOOTHD=1.
            END IF
!
          R(IK) =1+IS2PARS(16)*SMOOTHD*4*YOUNG*ICEH**3*(PI/WLG_I(IK))**4/(3*DWAT*GRAV*MP2)
!
! Converting action to surface elevation variance SUME with units m^2
!
          SUME = SUMA*DDEN(IK) / CG(IK)  / (R(IK)*CGRATIO(IK))
!
! CURVSPEC is the curvature variance = elevation variance * k^4
!
          CURVSPEC (IK) = SUME * (2*PI/ WLG_I(IK))**4
          ESPEC (IK) = SUME
          SUMALLDIR = SUMALLDIR / REAL(NTH)
!
! Adds the scattered energy isotropically to conserve the energy
! This may not be a very good scheme numerically. Another possible
! approach is the matrix inversion used for bottom scattering (w3sbs1md.ftn)
!
          S(1+(IK-1)*NTH:IK*NTH)=S(1+(IK-1)*NTH:IK*NTH)-SUMALLDIR*IS2PARS(2)
          END DO ! loop over wavenumbers IK
!
!  4. update of floe size
!
        IF (IS2PARS(10).LT.0.5) THEN ! resets max floe size to the last forcing or initial value
          DMAX = ICEDMAX
          ICEF = ICEDMAX
          END IF
!
        DO IK = 1, NK
! CURV is the variance of the curvature integrated over a finite bandwidth
          CURV(IK) = SUM(CURVSPEC(MAX(1,IK-NSUM):MIN(NK,IK+NSUM)))
! Now converts curvature variance to strain variance
          END DO ! end of loop on IK
!
! If IS2PARS(3)=IS2BREAK is set to true in ww3_grid, then activates ice break-up
!
        IF (IS2PARS(3).GT.0.5) THEN
          IKBREAK=0
          DO IK = 1, NK
            STRAIN = CURV(IK)*CURVTOSTRAIN
            IF (D_FLEX_FAIL .LT. DMAX) THEN
! Note that Williams et al. used IS2PARS(17)=SQRT(2), Here our default is 3.6
              WAMP(IK)= IS2PARS(17)*SQRT(STRAIN)
!
              IF (IS2PARS(9).EQ.0) THEN
                ICEDMIN=(TPI/WN_R(IK))*IS2PARS(14)
                END IF
              BRK_CRIT_W(IK) = WAMP(IK) .GT. STRAIN_C .AND. WLG_I(IK)/2 .GT. ICEDMIN .AND. WLG_I(IK)/2 .LT. DMAX &
                               .AND. WLG_I(IK)/2 .GT. D_FLEX_FAIL
 
!
              IF (BRK_CRIT_W(IK)) THEN
                IKBREAK=IK
                BRK_CRIT = .TRUE.
                END IF
              END IF
            END DO ! end of loop on IK
!
! 4.b Correction for bias introduced by the finite bandwidth sum .. .
!
          IF (BRK_CRIT) THEN
            DO IK=MAX(IKBREAK-NSUM,1),IKBREAK,1
! Modified by F.A. on Jan. 31, 2017: uses the maximum of CURVSPEC instead of CURV.
! this is better for very narrow spectra.
            IF (CURVSPEC(IK).GE.CURVSPEC(IKBREAK).AND.DMAX.GE.(WLG_I(IK)/2)) THEN
              IKBREAK = IK
              END IF
            END DO
!
            DMAX = WLG_I(IKBREAK)/2
!
!  Uses  a weighting by CURVSPEC to have a continuous shift of DMAX ..
!
            IKP1=MIN(IKBREAK+1,NK)
            IKM1=MAX(IKBREAK-1,1)
            IF (BRK_CRIT_W(IKP1).AND.BRK_CRIT_W(IKM1)) THEN
              DMAX = (WLG_I(IKBREAK)*CURVSPEC(IKBREAK)  &
                    +WLG_I(IKP1)*CURVSPEC(IKP1)+WLG_I(IKM1)*CURVSPEC(IKM1)) &
                    /(2.*(CURVSPEC(IKBREAK)+CURVSPEC(IKM1)+CURVSPEC(IKP1)))
              END IF
!	
            ICEF = DMAX
            END IF
          END IF  !end of test  (IS2PARS(3).GT.0.5)
!
!  5. inelastic or anelastic dissipation
!
        IF (IS2PARS(12).GT.0) THEN
          DO IK = 1, NK
!
! The TANH((DMAX-D*WLG_I(IK))/DMAX*C)
! is an ad hoc factor that goes to zero for WLG << DMAX and 1 for WLG >> DMAX
! this should probably be adjusted.
!
            IF (IS2PARS(14)*WLG_I(IK).LT.DMAX) THEN
              SMOOTHD=TANH((DMAX-IS2PARS(14)*WLG_I(IK))/(DMAX*IS2PARS(13)))
              IF (IS2PARS(23).LE.0.5) THEN  ! this is the inelastic option
                DISSIP(1+(IK-1)*NTH:IK*NTH)=CREEPFAC*4*CURV(IK)                &
                                            *((2*PI)/WLG_I(IK))**(IS2PARS(15)+1) &
                                            /(CGRATIO(IK)**1*R(IK)**2)           &
                                            *SMOOTHD
              ELSE ! this is the inelastic option
                DISSIP(1+(IK-1)*NTH:IK*NTH) =-4*4/3*SIG(IK)* DELTA_D*ALPHA_D *WN_I(IK)**4 * (YOUNG/MP2)**2         &
                                             * (ICEH/2)**3/3 * 1/( EXP(ALPHA_D*S_D(IK)) + EXP(-ALPHA_D*S_D(IK))) &
                                             * SMOOTHD /(R(IK)**2*CGRATIO(IK)) / (DWAT*GRAV) *TPIINV
              END IF
            END IF
            S=S+DISSIP*CICE*A
            END DO ! end of loop on IK
          ENDIF ! end of test (IS2PARS(12).GT.0)
!
!  6. Case of no scattering nor dissipation
!
      ELSE
        DMAX = 0.
        ICEF = 0.
        END IF   ! end of test  (CICE .GT. 0 .AND. ICEDAVE .GT. 0)
!
! Formats
      8000 FORMAT (' TEST W3SIS2 : ALPHA :',E10.3)
!
!/
!/ End of W3SIS2 ----------------------------------------------------- /
!/
      END SUBROUTINE W3SIS2
 
 
!/ ------------------------------------------------------------------- /
      SUBROUTINE W3RPWNICE(ICEH,WN_I,DAMPING,CG_I)
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |   C. Sevigny                      |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         27-Aug-2015 |
!/                  +-----------------------------------+
!/
!/    27-Aug-2015 : Origination.                        ( version 5.10 )
!
!  1. Purpose :
!
!/ ------------------------------------------------------------------- /
!
!  2. Method :
!
!  3. Parameters :
!
!     Parameter list
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!  5. Called by :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!  6. Error messages :
!
!     None.
!
!  7. Remarks :
!
!  8. Structure :
!
!     See source code.
!
!  9. Switches :
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE CONSTANTS, ONLY: GRAV
      USE W3GDATMD, ONLY: NK, SIG
 
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
      REAL, INTENT(IN)        :: ICEH
      REAL, INTENT(INOUT)  :: WN_I(:), DAMPING(:), CG_I(:)
 
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      DOUBLE COMPLEX :: WN_ROOT, GS1
      INTEGER :: IK
      REAL :: FLEX_RIGID
      REAL, PARAMETER :: VISC_RP = 10
      REAL, PARAMETER :: DENS = 1025.0
      REAL, PARAMETER :: DENS_ICE = 922.5
      REAL, PARAMETER :: POISSON = 0.3
      REAL, PARAMETER :: YOUNG = 5.49E+9
 
 
      FLEX_RIGID =  YOUNG * ICEH**3 /(12 *(1-POISSON**2) )
      ! Guess value for roots
      GS1 = CMPLX(SIG(1)**2/GRAV,0.)
 
      DO IK=1,NK
        CALL FINDROOTS_NR(GS1,0.,WN_ROOT,ICEH,SIG(IK))
        WN_I(IK) = REAL(WN_ROOT)
        CALL FINDROOTS_NR(GS1,VISC_RP,WN_ROOT,ICEH,SIG(IK))
        DAMPING(IK) = AIMAG(WN_ROOT)
        GS1 = WN_I(IK)
        CG_I(IK) = (5* FLEX_RIGID*WN_I(IK)**4 + DENS*GRAV - DENS_ICE*ICEH*SIG(IK)**2) &
              /(2*SIG(IK)*(DENS+DENS_ICE*WN_I(IK)*ICEH))
        END DO
 
!/
      END SUBROUTINE W3RPWNICE
 
 
!/ ------------------------------------------------------------------- /
!     SUBROUTINE FINDROOTS_NR(FUNCD,X0,VISC_RP,WN_ROOT)
      SUBROUTINE FINDROOTS_NR(GUESS,VISC_RP,X,ICEH,SIGMA)
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |   C. Sevigny                      |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         27-Aug-2015 |
!/                  +-----------------------------------+
!/
!/    27-Aug-2015 : Origination.                        ( version 5.10 )
!
!  1. Purpose :
!
!/ ------------------------------------------------------------------- /
!
!  2. Method :
!
!  3. Parameters :
!
!     Parameter list
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!  5. Called by :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!  6. Error messages :
!
!     None.
!
!  7. Remarks :
!
!  8. Structure :
!
!     See source code.
!
!  9. Switches :
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE W3ODATMD, ONLY: NDSE
      USE W3SERVMD, ONLY: EXTCDE
!/
!/ ------------------------------------------------------------------- /
      IMPLICIT NONE
 
!/ ------------------------------------------------------------------- /
!/ Parameter list
      REAL, INTENT(IN) :: VISC_RP, ICEH, SIGMA
      double COMPLEX, INTENT(IN) :: GUESS
      double COMPLEX, intent(INOUT) :: X
 
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      double complex :: FVAL, FDERIV, DX, X0
      INTEGER, PARAMETER :: MAXIT = 300
      REAL, PARAMETER :: TOL = 1E-9
      INTEGER :: J
      LOGICAL ::  UNFINISHED_ROOTS
 
      X0 = GUESS
      UNFINISHED_ROOTS = .TRUE.
      J = 1
      DO WHILE (J .LT. MAXIT .AND. UNFINISHED_ROOTS)
        FVAL = FUNCD_FVAL(X0,VISC_RP,ICEH,SIGMA)
        FDERIV = FUNCD_FDERIV(X0,VISC_RP,ICEH,SIGMA)
        DX = FVAL/FDERIV
        X = X0-DX
        X0 = X
        IF (ABS(DX) .GT. TOL .OR. ABS(FVAL) .GT. 1) THEN
          J = J+1
        ELSE
            UNFINISHED_ROOTS = .FALSE.
          END IF
 
        IF (J .GT. MAXIT) THEN
            WRITE (NDSE,1000)
            CALL EXTCDE ( 1 )
            END IF
        END DO
 
!
! Formats
      1000 FORMAT (/' *** ERROR FINDROOTS_NR *** '/                              &
               '     ROOT NOT CONVERGED'/)
 
      END SUBROUTINE FINDROOTS_NR
 
 
!/ ------------------------------------------------------------------- /
      FUNCTION W3FSD_DAVE(ICEDMIN,ICEDMAX,FRAGILITY)
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |   C. Sevigny & F. Ardhuin         |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         6-Nov-2015  |
!/                  +-----------------------------------+
!/
!/    27-Aug-2015 : Origination.                        ( version 5.10 )
!/     6-Nov-2015 : Uses a continuous DMAX->DAVE function(version 5.10 )
!
!  1. Purpose : Computes the mean flow size from the max floe size
!
!/ ------------------------------------------------------------------- /
!
!  2. Method :
!
!  3. Parameters :
!
!     Parameter list
!     ----------------------------------------------------------------
!       ICEDMIN    Real   I   Minimum floe diameter
!       FRAGILITY  Real   I   Parameter that gives the power in the FSD power law
!       ICEDMAX    R.A.  I/O  Maximum floe diameter
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!      None
!     ----------------------------------------------------------------
!
!  5. Called by :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!      INSIS2    Proc. W3SIS2MD Initialisation of parameters for IS2
!     ----------------------------------------------------------------
!
!  6. Error messages :
!
!     None.
!
!  7. Remarks :
!
!  8. Structure :
!
!     See source code.
!
!  9. Switches :
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
!
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
      REAL, INTENT(IN) :: ICEDMIN, ICEDMAX, FRAGILITY
      REAL             :: W3FSD_DAVE,  W3FSD_DAVE2
 
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
!
! analytic solution, if the FSD is given by P(x) = x^(-1-gam) for icedmin <= x <= icedmax and 0 elsewhere
!
       REAL                    :: GAM, MR, R, DENOM, ICEDMAXL
       REAL, PARAMETER :: xi = 2
!
       ICEDMAXL=MAX(ICEDMIN*1.01,ICEDMAX)
       GAM = 2 + log(FRAGILITY)/log(2.)
       IF (ICEDMIN.EQ.0) THEN
         W3FSD_DAVE =ICEDMAXL
       ELSE
         W3FSD_DAVE = GAM/(GAM-1) * ( (ICEDMAXL**(-GAM+1)-ICEDMIN**(-GAM+1))&
                      /(ICEDMAXL**(-GAM)-ICEDMIN**(-GAM)) )
      END IF
!
! Other possibility: analytical solution to Toyota algorithm (F. Arduin)
!       MR=log(ICEDMAXL/ICEDMIN)/log(2.)
!       R=(1-FRAGILITY*xi**2)/(1-FRAGILITY*xi)
!       DENOM = (1-(FRAGILITY*xi**2)**(MR+1.))
!       W3FSD_DAVE2 = MAX(ICEDMAXL*R* (1-(FRAGILITY*xi)**(MR+1.))/DENOM,ICEDMIN)
!  WRITE(991,*) ICEDMAX, GAM, W3FSD_DAVE, W3FSD_DAVE2
       END FUNCTION W3FSD_DAVE
!
!/ ------------------------------------------------------------------- /
      FUNCTION FUNCD_FVAL(WN_GUESS, VISC_RP,ICEH,SIGMA)  RESULT(FVAL)
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |   C. Sevigny                      |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         27-Aug-2015 |
!/                  +-----------------------------------+
!/
!/    27-Aug-2015 : Origination.                        ( version 5.10 )
!
!  1. Purpose :
!
!/ ------------------------------------------------------------------- /
!
!  2. Method :
!
!  3. Parameters :
!
!     Parameter list
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!  5. Called by :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!  6. Error messages :
!
!     None.
!
!  7. Remarks :
!
!  8. Structure :
!
!     See source code.
!
!  9. Switches :
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE CONSTANTS, ONLY: GRAV
 
      IMPLICIT NONE
 
      REAL, INTENT(IN) :: VISC_RP, ICEH, SIGMA
      DOUBLE COMPLEX, INTENT(IN) :: WN_GUESS
 
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      REAL :: ALP, FLEX_RIGID, GAMMA, C5, DRAFT, DENS_GRAV
      double complex :: C1, FVAL
 
      REAL, PARAMETER :: DENS = 1025.0
      REAL, PARAMETER :: DENS_ICE = 922.5
      REAL, PARAMETER :: POISSON = 0.3
      ! Effective modulus or strain modulus [Pa], YOUNG
      REAL, PARAMETER :: YOUNG = 5.49E+9
 
!/
!/ ------------------------------------------------------------------- /
 
      ! Length below which flexural felure cannot occur, flexural rigidity
      FLEX_RIGID =  YOUNG * ICEH**3 /(12 *(1-POISSON**2) )
      DRAFT = (DENS_ICE/DENS)*ICEH
      DENS_GRAV = DENS*GRAV
 
      ALP = SIGMA**2/GRAV
      ! Artificial viscosity (Robinson & Palmer, 1990)
      GAMMA = SIGMA * VISC_RP/DENS_GRAV
      C1 = CMPLX(1 - ALP*DRAFT, - 1.*GAMMA)
      C5 = FLEX_RIGID/DENS_GRAV
 
      FVAL = C5*WN_GUESS**5 + C1*WN_GUESS - ALP
 
      END FUNCTION FUNCD_FVAL
 
 
!/ ------------------------------------------------------------------- /
!     FUNCTION FUNCD(WN_GUESS, VISC_RP,ICEH)
      FUNCTION FUNCD_FDERIV(WN_GUESS, VISC_RP,ICEH,SIGMA)  RESULT(FDERIV)
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |   C. Sevigny                      |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         27-Aug-2015 |
!/                  +-----------------------------------+
!/
!/    27-Aug-2015 : Origination.                        ( version 5.10 )
!
!  1. Purpose :
!
!/ ------------------------------------------------------------------- /
!
!  2. Method :
!
!  3. Parameters :
!
!     Parameter list
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!  5. Called by :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!     ----------------------------------------------------------------
!
!  6. Error messages :
!
!     None.
!
!  7. Remarks :
!
!  8. Structure :
!
!     See source code.
!
!  9. Switches :
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE CONSTANTS, ONLY: GRAV
 
      IMPLICIT NONE
 
      REAL, INTENT(IN) :: VISC_RP, ICEH, SIGMA
      DOUBLE COMPLEX, INTENT(IN) :: WN_GUESS
 
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      REAL :: ALP, FLEX_RIGID, GAMMA, C5, DRAFT, DENS_GRAV
      double complex :: C1, FDERIV
 
      REAL, PARAMETER :: DENS = 1025.0
      REAL, PARAMETER :: DENS_ICE = 922.5
      REAL, PARAMETER :: POISSON = 0.3
      REAL, PARAMETER :: YOUNG = 5.49E+9
 
!/
!/ ------------------------------------------------------------------- /
!
! Length below which flexural felure cannot occur, flexural rigidity
!
      FLEX_RIGID =  YOUNG * ICEH**3 /(12 * (1-POISSON**2))
      DRAFT = (DENS_ICE/DENS)*ICEH
      DENS_GRAV = DENS*GRAV
 
      ALP = SIGMA**2/GRAV
!
! Artificial viscosity (Robinson & Palmer, 1990)
!
      GAMMA = SIGMA * VISC_RP/DENS_GRAV
      C1 = CMPLX(1 - ALP*DRAFT, -1.*GAMMA)
      C5 = FLEX_RIGID/DENS_GRAV
!
      FDERIV = 5*C5*WN_GUESS**4 +C1
!
      END FUNCTION FUNCD_FDERIV
!/
!/ End of module W3SIS1MD -------------------------------------------- /
!/
      END MODULE W3SIS2MD
 
