WAVEWATCH III  beta 0.0.1
Functions/Subroutines
w3sic5md Module Reference

Calculate ice source term S_{ice} according to different ice models: More...

Functions/Subroutines

subroutine, public w3sic5 (A, DEPTH, CG, WN, IX, IY, S, D)
 Calculate ice source term S_{ice} according to 3 different sea ice models. More...
 
subroutine, public w3ic5wncg (WN_R, WN_I, CG, HICE, IVISC, RHOI, ISMODG, HWAT)
 Calculation of complex wavenumber arrays for ice-coupled waves. More...
 
subroutine, public fsdisp (HICE, IVISC, RHOI, ISMODG, HWAT, WT, WNR, WNI)
 Calculate the complex wavenumber for waves in ice according to three different sea ice models. More...
 

Detailed Description

Calculate ice source term S_{ice} according to different ice models:

Author
Q. Liu
E. Rogers
Date
19-May-2021
Copyright
Copyright 2009-2022 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.

Function/Subroutine Documentation

◆ fsdisp()

subroutine, public w3sic5md::fsdisp ( real, intent(in)  HICE,
real, intent(in)  IVISC,
real, intent(in)  RHOI,
real, intent(in)  ISMODG,
real, intent(in)  HWAT,
real, intent(in)  WT,
real, intent(out)  WNR,
real, intent(out)  WNI 
)

Calculate the complex wavenumber for waves in ice according to three different sea ice models.

Parameters
[in]HICEThickness of ice (m).
[in]IVISCViscosity parameter of ice (m^2 s^{-1}).
[in]RHOIThe density of ice (kg m^{-3}).
[in]ISMODGEffecitive shear modulus G of ice (Pa).
[in]HWATWater depth.
[in]WTWave period (s; 1/freq).
[out]WNRThe real part of the wave number.
[out]WNIThe imaginary part of the wave number.
Author
Q. Liu
Date
19-May-2021

Definition at line 653 of file w3sic5md.F90.

653  !/
654  !/ +-----------------------------------+
655  !/ | WAVEWATCH III NOAA/NCEP |
656  !/ | Q. Liu |
657  !/ | FORTRAN 90 |
658  !/ | Last update : 19-May-2021 |
659  !/ +-----------------------------------+
660  !/
661  !/ 17-Mar-2016 : Origination ( version 5.10)
662  !/ ( Q. Liu)
663  !/ 17-Mar-2016 : Start from the Matlab code `FoxSquire.m` (provided
664  !/ by Prof. Vernon Squire from University of Otago)
665  !/ ( Q. Liu)
666  !/ 25-Apr-2017 : Add more filters ( Q. Liu)
667  !/
668  !/ 19-May-2021 : Incl. RP and M2 ice models ( Q. Liu)
669  !/
670  ! 1. Purpose :
671  !
672  ! Calculate the complex wavenumber for waves in ice according to
673  ! three different sea ice models, i.e., FS, RP and M2 (see Liu et
674  ! al. 2020)
675  !
676  ! 2. Method :
677  ! Mainly solving the dispersion relations of FS and RP models (
678  ! Eqs. (20, 24, 25)) in Mosig et al. (2015))
679  !
680  ! 3. Parameters :
681  !
682  ! Parameter list
683  ! ----------------------------------------------------------------
684  ! Name Type Intent Description
685  ! ----------------------------------------------------------------
686  ! HICE Real. IN thickness of ice (m)
687  ! IVISC Real. IN viscosity parameter of ice (m^2 s^{-1})
688  ! RHOI Real. IN the density of ice (kg m^{-3})
689  ! ISMODG Real. IN effecitive shear modulus G of ice (Pa)
690  ! HWAT Real. IN water depth
691  ! WT Real. IN wave period (s; 1/freq)
692  ! WNR Real. Out the real. part of the wave number
693  ! WNI Real. Out the imag. part of the wave number
694  ! ----------------------------------------------------------------
695  !
696  ! 4. Subroutines used :
697  !
698  ! Name Type Module Description
699  ! ----------------------------------------------------------------
700  ! STRACE Subr. W3SERVMD Subroutine tracing.
701  ! POLYROOTS Subr. / Find the roots of a general polynomial
702  ! NR_ROOT Func. / Newton-Raphson root finding
703  ! ----------------------------------------------------------------
704  !
705  ! 5. Called by :
706  !
707  ! Name Type Module Description
708  ! ----------------------------------------------------------------
709  ! W3IC5WNCG Subr. / Wavenumber and group velocity of ice-
710  ! coupled waves
711  ! ----------------------------------------------------------------
712  !
713  ! 6. Error messages :
714  !
715  ! See Format 1000, 1001, 1002
716  !
717  ! 7. Remarks :
718  !
719  ! 8. Structure :
720  !
721  ! See source code.
722  !
723  ! 9. Switches :
724  !
725  ! !/S Enable subroutine tracing.
726  !
727  ! 10. Source code :
728  !
729  !/ ------------------------------------------------------------------- /
730 #ifdef W3_S
731  USE w3servmd, ONLY: strace
732 #endif
733  !/
734  USE constants, ONLY: grav, tpi
735  USE w3dispmd, ONLY: wavnu1
736  USE w3servmd, ONLY: extcde
737  USE w3odatmd, ONLY: ndse, iaproc, naproc, naperr
738  USE w3gdatmd, ONLY: ic5pars
739  USE w3gsrumd, ONLY: w3inan
740  !/
741  IMPLICIT NONE
742  !/
743  !/ ------------------------------------------------------------------- /
744  !/ Parameter list
745  REAL, INTENT(IN) :: HICE, IVISC, RHOI, ISMODG, HWAT, WT
746  REAL, INTENT(OUT) :: WNR, WNI
747  !/
748  !/ ------------------------------------------------------------------- /
749  !/ Local parameters
750  !
751  REAL :: IC5MINIG, IC5MINWT, IC5MAXKRATIO, &
752  IC5MAXKI, IC5MINHW, IC5VEMOD
753  REAL :: TISMODG, TWT, TRATIO, THW
754  REAL, PARAMETER :: NU = 0.3, rhow = 1025.
755  ! COMPLEX :: GV, C1
756  ! REAL :: SIGMA, C2, WNO, CGO, THKH, &
757  COMPLEX :: GV, C1, C2
758  REAL :: SIGMA, WNO, CGO, THKH, &
759  RTRL(5), RTIM(5), RTANG(5)
760  INTEGER :: IREAL
761  ! COMPLEX(KDPC) :: GUESS, CROOT, C1D
762  ! REAL(KDP) :: C2D, HWATD
763  COMPLEX(KDPC) :: GUESS, CROOT, C1D, C2D
764  REAL(KDP) :: HWATD
765  !/
766 #ifdef W3_S
767  INTEGER, SAVE :: IENT = 0
768 #endif
769  !/
770  !/ ------------------------------------------------------------------- /
771  !/
772 #ifdef W3_S
773  CALL strace (ient, 'FSDISP')
774 #endif
775  ! Note, same as W3IC3WNCG_xx in w3sic3md :
776  ! HICE → ICE1
777  ! IVISC → ICE2
778  ! RHOI → ICE3
779  ! ISMODG → ICE4
780  ! 0. Initializations ------------------------------------------------ *
781  ! Set limiters
782  !
783  ! When G = 0, the FS method does not provide a solution. It is not
784  ! unexpected because the FS model is originally devised as a
785  ! thin elastic plate model in which elasticity is necessary.
786  !
787  ! The FS algorithm may also have issues for very short wave periods,
788  ! shallow waters and low G (e.g., T~3 s, d~10 m, hi~0.5 m, G<10^6 Pa)
789  !
790  ic5minig = ic5pars(1) ! Minimum G
791  ic5minwt = ic5pars(2) ! Minimum T
792  ic5maxkratio = ic5pars(3) ! Maximum k_{ow}/k_r
793  ic5maxki = ic5pars(4) ! Maximum k_i
794  ic5minhw = ic5pars(5) ! Minimum d
795  ic5vemod = ic5pars(9) ! Model selected 1: EFS, 2: RP, 3: M2
796  !
797  tismodg = max(ic5minig, ismodg)
798  twt = max(ic5minwt, wt)
799  thw = max(ic5minhw, hwat)
800  !
801  ! G <= 0. is not allowed
802  IF (abs(tismodg) < errtol) THEN
803  IF ( iaproc .EQ. naperr ) WRITE(ndse, 1000) 'FSDISP'
804  CALL extcde(1)
805  END IF
806  !
807  ! σ = 2π / T
808  sigma = tpi / twt
809  !
810  IF (abs(ic5vemod - 1.) < errtol) THEN
811  ! Complex shear modulus Gv = G - i σ ρ η (EFS model)
812  gv = cmplx(tismodg, -1. * sigma * rhoi * ivisc)
813  !
814  ! -------------------------------------------------------------------- *
815  ! Note that Eq. (24) in Mosig et al. (2015) can be written like below:
816  ! (c1 * k^5 + c2 * k) * tanh(HWAT*k) - 1 = 0
817  ! Most Important part of this module --------------------------------- *
818  c1 = gv * hice**3. / (6. * rhow * sigma**2.)
819  !
820  ! To be divided by (1-NU) or multiplied by (1+NU) ??
821  ! Beam model: then multiplied by (1+ν)
822  ! Plate model: then divided by (1-ν)
823  ! The beam version is more theoretically (J.E.M. Mosig, personal
824  ! communication, 2016), although there is only very marginal difference
825  ! between this two version as (1+NU = 1.3 and 1/(1-NU) ~ 1.4)
826  c1 = c1 * (1+nu)
827  ! C1 = C1 / (1-NU)
828  !
829  ! C2
830  ! C2 = GRAV / SIGMA**2. - RHOI * HICE / RHOW
831  c2 = cmplx(grav / sigma**2. - rhoi * hice / rhow, 0.)
832  !
833  ELSE IF (abs(ic5vemod - 2.) < errtol) THEN
834  ! See Appendix of Liu et al. (2020) - RP model
835  c1 = cmplx(tismodg * hice**3. * (1+nu) / (6. * rhow * sigma**2.), 0.)
836  c2 = cmplx(grav/sigma**2. - rhoi * hice / rhow, &
837  -1. * ivisc / (rhow * sigma))
838  !
839  ELSE IF (abs(ic5vemod - 3.) > errtol) THEN
840  WRITE(ndse, 1003) 'FSDISP', ic5vemod
841  CALL extcde(4)
842  END IF
843  ! Use the dispersion in open water to get an approximation of
844  ! tanh(HWAT * k). We can also roughly use the dispersion in deep
845  ! water case, that is tanh(HWAT*k) ~ 1.
846  ! Wavenumber in the open water
847  ! WAVNU1(SI, H, K, CG)
848  CALL wavnu1(sigma, thw, wno, cgo)
849  thkh = tanh(wno * thw)
850  !
851  IF (abs(ic5vemod - 1.) < errtol .OR. abs(ic5vemod - 2.) < errtol) THEN
852  ! Get the first guess of the complex wavenumber
853  CALL polyroots(6, &
854  (/real(real(c1))*thkh, 0., 0., 0., real(real(c2))*thkh, -1./),&
855  rtrl, rtim)
856  rtang = atan2(rtim, rtrl)
857  !
858  ! There should only be one real root in RTRL + i * RTIM because in
859  ! this case (ivisc=0) the original viscoelastic-type model reduced to
860  ! the thin elastic plate model which has only one real solution.
861  ! Find its index ...
862  !
863  ireal = minloc(abs(rtang), dim=1)
864  IF (rtrl(ireal) <= 0. .OR. abs(rtim(ireal)) > errtol) THEN
865  IF ( iaproc .EQ. naperr ) WRITE(ndse, 1001) 'FSDISP'
866  CALL extcde(2)
867  END IF
868  !
869  ! Get the first guess for iteration
870  guess = rtrl(ireal) * exp(cmplx(0., 1e-6))
871  !
872  ! Newton-Raphson method
873  ! Turn c1, c2, hwat to be double
874  c1d = c1; c2d = c2; hwatd = thw
875  croot = nr_root(c1d, c2d, hwatd, guess)
876  wnr = real(real(croot))
877  wni = real(aimag(croot))
878  !
879  ELSE IF (abs(ic5vemod - 3.) < errtol) THEN ! M2
880  ! Model with Order 3 Power Law (section 6.2 in Meylan et al. (2018, JGR-Ocean))
881  ! Based on my understanding, the wavelength does not change because
882  ! the elasticity is not considered in this model
883  wnr = wno ! Open-water wavenumber
884  ! Eq. (53) in Meylan et al. (2018)
885  wni = hice * ivisc * sigma**3. / (rhow * grav**2.)
886  END IF
887  !
888  ! RATIO Check
889  ! Using the ratio k0 / kr as a basic check for the reliability of
890  ! FSDISP. The FS dispersion relation can give a very different kr from
891  ! k0, especially for small wave periods (k0/kr is as high as 100).
892  ! From my tests, using IC5MAXKRATIO = 1000. can basically detect most
893  ! spurious solutions (although not all of them)
894  !
895  ! ISNAN Check
896  ! Common ways used are:
897  ! NAN = SQRT(-1.) or
898  ! a /= a then a is NaN or
899  ! ISNAN func (supported by gfortran & ifort)
900  ! --- ISNAN -> W3INAN because ISNAN is not supported by pgi
901  ! For very few cases, we can get nan | negative ki | kr
902  !
903  ! (N.B.) NaN problem solved by using CMPLX_TANH2
904  !
905  tratio = wno / wnr
906  IF (w3inan(wnr) .OR. w3inan(wni) .OR. wnr <= 0 .OR. wni <= 0. &
907  .OR. tratio >= ic5maxkratio) THEN
908  IF ( iaproc .EQ. naperr ) &
909  WRITE(ndse, 1002) 'FSDISP', hice, ivisc, tismodg, hwat, twt, &
910  wno, wnr, wni
911  CALL extcde(3)
912  END IF
913  !
914  ! Filter high ki
915  wni = min(ic5maxki, wni)
916  !
917  ! FORMAT
918 1000 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/ &
919  ' Subr. ', a, ': Zero shear modulus G is not allowed&
920  & in the FS viscoelastic model'/)
921  !
922 1001 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/ &
923  ' Subr. ', a, ': get a bad first guess'/)
924  !
925 1002 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/ &
926  ' -----------------------------------------------------'/&
927  ' Subr. ', a,' : get NaN/NeG/Huge kr or ki for' /&
928  ' -----------------------------------------------------'/&
929  ' Ice thickness : ', f9.1, ' m'/ &
930  ' Ice viscosity : ', e9.2, ' m2/s'/ &
931  ' Ice shear modulus : ', e9.2, ' Pa' / &
932  ' Water depth : ', f9.1, ' m'/ &
933  ' Wave period : ', f10.2, ' s'/ &
934  ' Wave number (Ko) : ', f11.3, ' rad/m'/ &
935  ' Wave number (Kr) : ', f11.3, ' rad/m'/ &
936  ' Attenu. Rate (Ki) : ', e9.2, ' /m'/)
937  !
938 1003 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/ &
939  ' Subr. ', a, ': Unknown VE model (', f5.0, ')'/)
940  !/
941  !/ End of FSDISP ----------------------------------------------------- /
942  !/

References w3servmd::extcde(), constants::grav, w3odatmd::iaproc, w3gdatmd::ic5pars, w3odatmd::naperr, w3odatmd::naproc, w3odatmd::ndse, w3servmd::strace(), constants::tpi, and w3dispmd::wavnu1().

Referenced by w3ic5wncg().

◆ w3ic5wncg()

subroutine, public w3sic5md::w3ic5wncg ( real, dimension(:), intent(inout)  WN_R,
real, dimension(:), intent(inout)  WN_I,
real, dimension(:), intent(in)  CG,
real, intent(in)  HICE,
real, intent(in)  IVISC,
real, intent(in)  RHOI,
real, intent(in)  ISMODG,
real, intent(in)  HWAT 
)

Calculation of complex wavenumber arrays for ice-coupled waves.

Using the Fox-Squire dispersion relations to get (kr, ki) and then get cg by cg = dσ / dk (here dk uses kr).

Parameters
[in,out]WN_RThe real part of the wave number.
[in,out]WN_IThe imaginary part of the wave number.
[in,out]CGGroup velocity (m s^{-1}).
[in,out]HICEThickness of ice (m).
[in,out]IVISCViscosity parameter of ice (m^2 s^{-1}).
[in,out]RHOIThe density of ice (kg m^{-3}).
[in,out]ISMODGEffecitive shear modulus G of ice (Pa).
[in,out]HWATWater depth.
Author
Q. Liu
E. Rogers
Date
25-Apr-2017

Definition at line 495 of file w3sic5md.F90.

495  !/
496  !/ +-----------------------------------+
497  !/ | WAVEWATCH III NOAA/NCEP |
498  !/ | Q. Liu |
499  !/ | E. Rogers |
500  !/ | FORTRAN 90 |
501  !/ | Last update : 25-Apr-2017 |
502  !/ +-----------------------------------+
503  !/
504  !/ 17-Apr-2016 : Origination ( version 5.10)
505  !/ ( Q. Liu )
506  !/ 17-Apr-2016 : Start from W3IC3WNCG_CHENG ( Q. Liu )
507  !/
508  !/ 1. Purpose:
509  ! Calculation of complex wavenumber arrays for ice-coupled waves.
510  !
511  ! This also allows us to use Cg_ice in the advection part of the
512  ! radiative transfer energy equation (RTE). --- abandoned in the end
513  !
514  ! 2. Method:
515  ! Using the Fox-Squire dispersion relations to get (kr, ki) and
516  ! then get cg by cg = dσ / dk (here dk uses kr)
517  !
518  ! 3. Parameters:
519  !
520  ! Parameter list:
521  ! ----------------------------------------------------------------
522  ! Name Type Intent Description
523  ! ----------------------------------------------------------------
524  ! WN_R R.A. I/O the real. part of the wave number
525  ! WN_I R.A. I/O the imag. part of the wave number
526  ! CG R.A. I group velocity (m s^{-1})
527  ! HICE Real. I thickness of ice (m)
528  ! IVISC Real. I viscosity parameter of ice (m^2 s^{-1})
529  ! RHOI Real. I the density of ice (kg m^{-3})
530  ! ISMODG Real. I effecitive shear modulus G of ice (Pa)
531  ! HWAT Real. I water depth
532  ! ----------------------------------------------------------------
533  ! * the intent of WN_R/I must be inout
534  ! * CG is unchanged but still kept here because some legacy reasons.
535  !
536  ! 4. Subroutines used:
537  !
538  ! Name Type Module Description
539  ! ----------------------------------------------------------------
540  ! FSDISP Subr. / dispersion relations for ice-coupled waves
541  ! CGINICE5 Subr. / group velocity for given (σ, kr) array
542  ! ----------------------------------------------------------------
543  !
544  ! 5. Called by :
545  !
546  ! Name Type Module Description
547  ! ----------------------------------------------------------------
548  ! W3SIC5 Subr. Public Ice source term
549  ! W3WAVE Subr. W3WAVEMD WW3 integration
550  ! ----------------------------------------------------------------
551  !
552  ! 6. Error messages :
553  !
554  ! 7. Remarks :
555  !
556  ! 8. Structure :
557  !
558  ! See source code.
559  !
560  ! 9. Switches :
561  !
562  ! !/S Enable subroutine tracing.
563  !
564  ! 10. Source code :
565  !
566  !/ ------------------------------------------------------------------- /
567 #ifdef W3_S
568  USE w3servmd, ONLY: strace
569 #endif
570  USE constants, ONLY: tpi
571  USE w3gdatmd, ONLY: nk, sig
572  USE w3odatmd, ONLY: ndse, iaproc, naproc, naperr
573  USE w3servmd, ONLY: extcde
574  !/
575  IMPLICIT NONE
576  !/
577  !/ ------------------------------------------------------------------- /
578  !/ Parameter list
579  REAL, INTENT(INOUT) :: WN_R(:), WN_I(:)
580  REAL, INTENT(IN) :: CG(:)
581  REAL, INTENT(IN) :: HICE, IVISC, RHOI, ISMODG, HWAT
582  !/
583  !/ ------------------------------------------------------------------- /
584  !/ Local parameters
585  REAL, ALLOCATABLE :: SIGMA(:)
586  INTEGER :: KL, KU, IK
587  REAL :: TWN_R, TWN_I
588  !/
589 #ifdef W3_S
590  INTEGER, SAVE :: IENT = 0
591 #endif
592  !/
593  !/ ------------------------------------------------------------------- /
594  !/
595 #ifdef W3_S
596  CALL strace (ient, 'W3IC5WNCG')
597 #endif
598  !/
599  ! Initialize SIGMA {in w3gdatmd: SIG (0: NK+1)}
600  IF (ALLOCATED(sigma)) DEALLOCATE(sigma); ALLOCATE(sigma(SIZE(cg)))
601  sigma = 0.
602 
603  IF (SIZE(wn_r, 1) .EQ. nk) THEN
604  kl = 1
605  ku = nk
606  sigma = sig(1:nk)
607  ELSE IF (SIZE(wn_r,1) .EQ. nk+2) THEN
608  kl = 1
609  ku = nk+2
610  sigma = sig(0:nk+1)
611  ELSE
612  IF ( iaproc .EQ. naperr ) WRITE(ndse,900) 'W3IC5WNCG'
613  CALL extcde(3)
614  END IF
615  !
616  ! Fox-Squire dispersion
617  DO ik = kl, ku
618  ! FSDISP(HICE, IVISC, RHOI, ISMODG, HWAT, WT, WNR, WNI)
619  CALL fsdisp(hice, ivisc, rhoi, ismodg, hwat, tpi/sigma(ik), &
620  twn_r, twn_i)
621  wn_r(ik) = twn_r
622  wn_i(ik) = twn_i
623  END DO
624  !
625  DEALLOCATE(sigma)
626  !
627 900 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/ &
628  ' Subr. ', a, ': Cannot determine bounds of&
629  & wavenumber array.'/)
630  !/
631  !/ End of W3IC5WNCG -------------------------------------------------- /
632  !/

References w3servmd::extcde(), fsdisp(), w3odatmd::iaproc, w3odatmd::naperr, w3odatmd::naproc, w3odatmd::ndse, w3gdatmd::nk, w3gdatmd::sig, w3servmd::strace(), and constants::tpi.

Referenced by w3sic5().

◆ w3sic5()

subroutine, public w3sic5md::w3sic5 ( real, dimension(nspec), intent(in)  A,
real, intent(in)  DEPTH,
real, dimension(nk), intent(in)  CG,
real, dimension(nk), intent(in)  WN,
integer, intent(in)  IX,
integer, intent(in)  IY,
real, dimension(nspec), intent(out)  S,
real, dimension(nspec), intent(out)  D 
)

Calculate ice source term S_{ice} according to 3 different sea ice models.

(Mosig et al. 2015, Meylan et al. 2018, Liu et al. 2020).

Parameters
[in]AAction density spectrum (1-D).
[in]DEPTHLocal water depth.
[in]CGGroup velocities.
[in]WNWavenumbers.
[in]IXGrid index.
[in]IYGrid index.
[out]SSource term (1-D version).
[out]DDiagonal term of derivative (1-D version).
Author
Q. Liu
E. Rogers
Date
19-May-2021

Definition at line 169 of file w3sic5md.F90.

169  !/
170  !/ +-----------------------------------+
171  !/ | WAVEWATCH III NOAA/NCEP |
172  !/ | Q. Liu |
173  !/ | E. Rogers |
174  !/ | FORTRAN 90 |
175  !/ | Last update : 19-May-2021 |
176  !/ +-----------------------------------+
177  !/
178  !/ 23-Mar-2016 : Origination ( version 5.10 )
179  !/ ( Q. Liu)
180  !/ 23-Mar-2016 : Started from w3sic1/2/3/4 subr. ( Q. Liu)
181  !/ 05-Apr-2016 : Options for Cg_{ice} or Cg ( Q. Liu)
182  !/ 25-Apr-2017 : Add more filters ( Q. Liu)
183  !/ 20-Aug-2018 : Ready to be merged to master (v6.06)( Q. Liu)
184  !/
185  !/ Copyright 2009 National Weather Service (NWS),
186  !/ National Oceanic and Atmospheric Administration. All rights
187  !/ reserved. WAVEWATCH III is a trademark of the NWS.
188  !/ No unauthorized use without permission.
189  !/
190  !/ 1. Purpose :
191  ! Calculate ice source term S_{ice} according to 3 different sea ice
192  ! models (Mosig et al. 2015, Meylan et al. 2018, Liu et al. 2020)
193  !
194  ! 2. Method :
195  ! Regarding i/o (general to all Sice modules): S_{ice} source term
196  ! is calculated using up to 5 parameters read from input files.
197  ! These parameters are allowed to vary in space and time.
198  ! The parameters control the exponential decay rate k_i.
199  ! Since there are 5 parameters, this permits description of
200  ! dependence of k_i on frequency or wavenumber.
201  !
202  ! Sea ice affects the wavenumber k of wind-generated ocean waves.
203  ! The ice-modified wavenumber can be expressed as a complex number
204  ! k = k_r + i * k_i, with the real part k_r representing impact of
205  ! the sea ice on the physical wavelength and propagation speeds,
206  ! producing something analogous to shoaling and refraction by
207  ! bathymetry, whereas the imaginary part of the complex
208  ! wavenumber, k_i, is an exponential decay coefficient
209  ! k_i(x,y,t,sigma) (depending on location, time and frequency,
210  ! respectively), representing wave attenuation, and can be
211  ! introduced in a wave model such as WW3 as S_ice/E=-2*Cg*k_i,
212  ! where S_ice is one of several dissipation mechanisms, along
213  ! with whitecapping, for example, S_ds=S_wc+S_ice+⋯. The k_r -
214  ! modified by ice would enter the model via the C calculations
215  ! on the left-hand side of the governing equation.The fundamentals
216  ! are straightforward, e.g. Rogers and Holland (2009 and
217  ! subsequent unpublished work) modified a similar model, SWAN
218  ! (Booij et al. 1999) to include the effects of a viscous mud
219  ! layer using the same approach (k = k_r + i*k_i) previously.
220  !
221  ! General approach is analogous to Rogers and Holland (2009)
222  ! approach for mud.
223  ! See text near their eq. 1 :
224  ! k = k_r + i * k_i
225  ! eta(x,t) = Real( a * exp( i * ( k * x - sigma * t ) ) )
226  ! a = a0 * exp( -k_i * x )
227  ! S / E = -2 * Cg * k_i (see also Komen et al. (1994, pg. 170)
228  !
229  ! Following W3SBT1 as a guide, equation 1 of W3SBT1 says:
230  ! S = D * E
231  ! However, the code of W3SBT1 has
232  ! S = D * A
233  ! This leads me to believe that the calling routine is
234  ! expecting "S/sigma" not "S"
235  ! Thus we will use D = S/E = -2 * Cg * k_i
236  ! (see also documentation of W3SIC1)
237  !
238  ! Notes regarding numerics:
239  ! -------------------------
240  ! Experiments with constant k_i values suggest that results may be
241  ! dependent on resolution if insufficient resolution is used.
242  ! For detailed information, see documentation of W3SIC1.
243  !
244  ! Note regarding applicability/validity:
245  ! --------------------------------------
246  ! Similar to the Wang and Shen model used in w3sic3md, the 3 models
247  ! used here are empirical medium models as well which treat the sea
248  ! ice cover as a continuum and use 1/2 empirical rheological para-
249  ! meters, i.e., the effective shear modulus of ice G and the effec-
250  ! tive viscosity η to characterize sea ices of various type. Please
251  ! see the documentation of w3sic3md for a detailed discussion of
252  ! this kind of model.
253  !
254  ! 3. Parameters :
255  !
256  ! Parameter list
257  ! ----------------------------------------------------------------
258  ! A R.A. I Action density spectrum (1-D).
259  ! DEPTH Real I Local water depth.
260  ! CG R.A. I Group velocities.
261  ! WN R.A. I Wavenumbers
262  ! IX,IY I.S. I Grid indices.
263  ! S R.A. O Source term (1-D version).
264  ! D R.A. O Diagonal term of derivative (1-D version).
265  ! ----------------------------------------------------------------
266  !
267  ! 4. Subroutines used :
268  !
269  ! Name Type Module Description
270  ! ----------------------------------------------------------------
271  ! STRACE Subr. W3SERVMD Subroutine tracing (!/S switch).
272  ! PRT2DS Subr. W3ARRYMD Print plot output (!/T0 switch).
273  ! OUTMAT Subr. W3ARRYMD Matrix output (!/T1 switch).
274  ! W3IC5WNCG Subr. / Wavenumber and group velocity of ice-
275  ! coupled waves
276  ! ----------------------------------------------------------------
277  ! * / means this module
278  !
279  ! 5. Called by :
280  !
281  ! Name Type Module Description
282  ! ----------------------------------------------------------------
283  ! W3SRCE Subr. W3SRCEMD Source term integration.
284  ! W3EXPO Subr. N/A ASCII Point output post-processor.
285  ! W3EXNC Subr. N/A NetCDF Point output post-processor.
286  ! GXEXPO Subr. N/A GrADS point output post-processor.
287  ! ----------------------------------------------------------------
288  !
289  ! 6. Error messages :
290  !
291  ! None.
292  !
293  ! 7. Remarks :
294  !
295  ! If ice parameter 1 is zero, no calculations are made.
296  !
297  ! 8. Structure :
298  !
299  ! See source code.
300  !
301  ! 9. Switches :
302  !
303  ! !/S Enable subroutine tracing.
304  ! !/T Enable general test output.
305  ! !/T0 2-D print plot of source term.
306  ! !/T1 Print arrays.
307  !
308  ! 10. Source code :
309  !/ ------------------------------------------------------------------- /
310  !/
311 #ifdef W3_T
312  USE w3odatmd, ONLY: ndst
313 #endif
314 #ifdef W3_S
315  USE w3servmd, ONLY: strace
316 #endif
317 #ifdef W3_T0
318  USE w3arrymd, ONLY: prt2ds
319 #endif
320 #ifdef W3_T1
321  USE w3arrymd, ONLY: outmat
322 #endif
323  !/
324  USE constants, ONLY: tpi
325  USE w3servmd, ONLY: extcde
326  USE w3odatmd, ONLY: ndse, iaproc, naproc, naperr
327  USE w3gdatmd, ONLY: nk, nth, nspec, sig, mapwn, ic5pars
328  USE w3idatmd, ONLY: inflags2, icep1, icep2, icep3, icep4, icei
329  !
330  IMPLICIT NONE
331  !
332  !/
333  !/ ------------------------------------------------------------------- /
334  !/ Parameter list
335  !/
336  REAL, INTENT(IN) :: CG(NK), WN(NK), A(NSPEC), DEPTH
337  REAL, INTENT(OUT) :: S(NSPEC), D(NSPEC)
338  INTEGER, INTENT(IN) :: IX, IY
339  !/ ------------------------------------------------------------------- /
340  !/ Local parameters
341  !/
342 #ifdef W3_S
343  INTEGER, SAVE :: IENT = 0
344 #endif
345 #ifdef W3_T0
346  INTEGER :: ITH
347  REAL :: DOUT(NK,NTH)
348 #endif
349  !/
350  REAL :: ICECOEF1, ICECOEF2, ICECOEF3, &
351  ICECOEF4, ICECONC
352  REAL, DIMENSION(NK) :: D1D, WN_R, WN_I
353  ! REAL :: TWN_R, TWN_I
354  INTEGER :: IK, IKTH
355  LOGICAL :: NOICE
356  !/ ------------------------------------------------------------------- /
357  !/
358 #ifdef W3_S
359  CALL strace (ient, 'W3SIC5')
360 #endif
361  !
362  ! 0. Initializations ------------------------------------------------ /
363  d = 0.
364  d1d = 0.
365  wn_r = 0.
366  wn_i = 0.
367  !
368  icecoef1 = 0.
369  icecoef2 = 0.
370  icecoef3 = 0.
371  icecoef4 = 0.
372  iceconc = 0.
373  !
374  ! Set the ice parameters from input
375  IF (inflags2(-7)) THEN
376  icecoef1 = icep1(ix, iy) ! ice thickness h_i
377  ELSE
378  IF ( iaproc .EQ. naperr ) &
379  WRITE (ndse,1001) 'ICE PARAMETER 1 (HICE)'
380  CALL extcde(2)
381  ENDIF
382  !
383  IF (inflags2(-6)) THEN
384  icecoef2 = icep2(ix, iy) ! effective viscosity of ice η
385  ELSE
386  IF ( iaproc .EQ. naperr ) &
387  WRITE (ndse,1001) 'ICE PARAMETER 2 (VISC)'
388  CALL extcde(2)
389  ENDIF
390  !
391  IF (inflags2(-5)) THEN
392  icecoef3 = icep3(ix, iy) ! density of ice ρ_i
393  ELSE
394  IF ( iaproc .EQ. naperr ) &
395  WRITE (ndse,1001) 'ICE PARAMETER 3 (DENS)'
396  CALL extcde(2)
397  ENDIF
398  !
399  IF (inflags2(-4)) THEN
400  icecoef4 = icep4(ix, iy) ! effective shear modulus of ice G
401  ELSE
402  IF ( iaproc .EQ. naperr ) &
403  WRITE (ndse,1001) 'ICE PARAMETER 4 (ELAS)'
404  CALL extcde(2)
405  ENDIF
406  !
407  IF (inflags2(4)) iceconc = icei(ix, iy) ! ice concentration
408  !
409  ! 1. No ice --------------------------------------------------------- /
410  noice = .false.
411  ! Zero ice thickness
412  ! Very small ice thickness may cause problems in POLYROOTS because
413  ! the first coefficient C1 may be very close to zero. So we regard
414  ! cases where hice is less than 0.0001 as no ice.
415  ! IF (ICECOEF1 < ERRTOL) NOICE = .TRUE.
416  IF (icecoef1 < 0.0001) noice = .true.
417  ! zero ice concentration
418  IF (inflags2(4) .AND. iceconc < errtol) noice = .true.
419  !
420  ! Calculate the decay rate k_i
421  IF ( noice ) THEN
422  d1d = 0.
423  !
424  ! 2. Ice ------------------------------------------------------------- /
425  ELSE
426  ! W3IC5WNCG(WN_R, WN_I, CG, HICE, IVISC, RHOI, ISMODG, HWAT)
427  CALL w3ic5wncg(wn_r, wn_i, cg, icecoef1, icecoef2, &
428  icecoef3, icecoef4, depth)
429  ! recall that D=S/E=-2*Cg_{ice}*k_i
430  ! In some cases, the FS model yields very large Cg_{ice}, which
431  ! subquently may result in numerical failure due to the violation of CFL
432  ! conditions, therefore we still use ice-free group velocity to advect
433  ! wave packets.
434  !
435  DO ik = 1, nk
436  d1d(ik) = -2.0 * cg(ik) * wn_i(ik)
437  END DO
438  END IF
439  !
440  ! 2.1 Fill diagonal matrix
441  DO ikth = 1, nspec
442  d(ikth) = d1d(mapwn(ikth))
443  END DO
444 
445  s = d * a
446  !
447  ! ... Test output of arrays
448  !
449 #ifdef W3_T0
450  DO ik=1, nk
451  DO ith=1, nth
452  dout(ik,ith) = d(ith+(ik-1)*nth)
453  END DO
454  END DO
455  CALL prt2ds (ndst, nk, nk, nth, dout, sig(1:), ' ', 1., &
456  0.0, 0.001, 'Diag Sice', ' ', 'NONAME')
457 #endif
458  !
459 #ifdef W3_T1
460  CALL outmat (ndst, d, nth, nth, nk, 'diag Sice')
461 #endif
462  !
463  ! Formats
464  !
465 1001 FORMAT(/' *** WAVEWATCH III ERROR IN W3SIC5MD : '/ &
466  ' ',a,' IS NOT DEFINED IN ww3_shel.inp.')
467 
468  !/
469  !/ End of W3SIC5------------------------------------------------------ /
470  !/

References w3servmd::extcde(), w3odatmd::iaproc, w3gdatmd::ic5pars, w3idatmd::inflags2, w3gdatmd::mapwn, w3odatmd::naperr, w3odatmd::naproc, w3odatmd::ndse, w3odatmd::ndst, w3gdatmd::nk, w3gdatmd::nspec, w3gdatmd::nth, w3arrymd::outmat(), w3arrymd::prt2ds(), w3gdatmd::sig, w3servmd::strace(), constants::tpi, and w3ic5wncg().

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

w3gdatmd::nk
integer, pointer nk
Definition: w3gdatmd.F90:1230
w3gdatmd::nspec
integer, pointer nspec
Definition: w3gdatmd.F90:1230
w3idatmd::inflags2
logical, dimension(:), pointer inflags2
Definition: w3idatmd.F90:260
w3gsrumd
Definition: w3gsrumd.F90:17
w3gdatmd::sig
real, dimension(:), pointer sig
Definition: w3gdatmd.F90:1234
w3gdatmd::ic5pars
real, dimension(:), pointer ic5pars
Definition: w3gdatmd.F90:1158
w3odatmd::iaproc
integer, pointer iaproc
Definition: w3odatmd.F90:457
m_constants::nu
real nu
kinematic viscosity of water
Definition: mod_constants.f90:22
w3arrymd::outmat
subroutine outmat(NDS, A, MX, NX, NY, MNAME)
Definition: w3arrymd.F90:988
w3odatmd::ndse
integer, pointer ndse
Definition: w3odatmd.F90:456
w3odatmd::naperr
integer, pointer naperr
Definition: w3odatmd.F90:457
w3servmd
Definition: w3servmd.F90:3
w3gdatmd::nth
integer, pointer nth
Definition: w3gdatmd.F90:1230
w3odatmd
Definition: w3odatmd.F90:3
w3odatmd::naproc
integer, pointer naproc
Definition: w3odatmd.F90:457
constants::tpi
real, parameter tpi
TPI 2*Pi.
Definition: constants.F90:72
w3servmd::strace
subroutine strace(IENT, SNAME)
Definition: w3servmd.F90:148
w3idatmd
Define data structures to set up wave model input data for several models simultaneously.
Definition: w3idatmd.F90:16
w3gdatmd::mapwn
integer, dimension(:), pointer mapwn
Definition: w3gdatmd.F90:1231
w3arrymd
Definition: w3arrymd.F90:3
w3odatmd::ndst
integer, pointer ndst
Definition: w3odatmd.F90:456
constants
Define some much-used constants for global use (all defined as PARAMETER).
Definition: constants.F90:20
w3gdatmd
Definition: w3gdatmd.F90:16
w3dispmd::wavnu1
subroutine wavnu1(SI, H, K, CG)
Definition: w3dispmd.F90:85
w3sic5md::fsdisp
subroutine, public fsdisp(HICE, IVISC, RHOI, ISMODG, HWAT, WT, WNR, WNI)
Calculate the complex wavenumber for waves in ice according to three different sea ice models.
Definition: w3sic5md.F90:653
w3servmd::extcde
subroutine extcde(IEXIT, UNIT, MSG, FILE, LINE, COMM)
Definition: w3servmd.F90:736
w3arrymd::prt2ds
subroutine prt2ds(NDS, NFR0, NFR, NTH, E, FR, UFR, FACSP, FSC, RRCUT, PRVAR, PRUNIT, PNTNME)
Definition: w3arrymd.F90:1943
w3sic5md::w3ic5wncg
subroutine, public w3ic5wncg(WN_R, WN_I, CG, HICE, IVISC, RHOI, ISMODG, HWAT)
Calculation of complex wavenumber arrays for ice-coupled waves.
Definition: w3sic5md.F90:495
w3dispmd
Definition: w3dispmd.F90:3
constants::grav
real, parameter grav
GRAV Acc.
Definition: constants.F90:61