#include "swancpp.h"
!
! SWAN/COMPU file 4 of 5
!
!
! PROGRAM SWANCOM4.FOR
!
!
! This file SWANCOM4 of the main program SWAN
! include the next subroutines
!
! *** nonlinear 4 wave-wave interactions ***
!
! FAC4WW (compute the constants for the nonlinear wave
! interactions)
! RANGE4 (compute the counters for the different types of
! computations for the nonlinear wave interactions)
! SWSNL1 (nonlinear four wave interactions; semi-implicit and computed
! for all bins that fall within a sweep with DIA technique.
! Interaction are calculated per sweep)
! SWSNL2 (nonlinear four wave interactions; fully explicit and computed
! for all bins that fall within a sweep with DIA technique.
! Interaction are calculated per sweep)
! SWSNL3 (calculate nonlinear four wave interactions fully explicitly
! for the full circle per iteration by means of DIA approach
! and store results in auxiliary array MEMNL4)
! SWSNL4 (calculate nonlinear four wave interactions fully explicitly
! for the full circle per iteration by means of MDIA approach
! and store results in auxiliary array MEMNL4)
! SWSNL8 (calculate nonlinear four wave interactions fully explicitly
! for the full circle per iteration by means of DIA approach
! and store results in auxiliary array MEMNL4. Neighbouring
! interactions are interpolated in piecewise constant manner)
! FILNL3 (fill main diagonal and right-hand side of the system with
! results of array MEMNL4)
!
! SWINTFXNL (interface with SWAN model to compute nonlinear transfer
! with the XNL method for given action density spectrum)
!
! *** nonlinear 3 wave-wave interactions ***
!
! TCOEF (compute coupling coefficients for SPB of Becq-Girard et al, 1999)
! COAPP (interpolation of energy densities at sum and diff freqs for triads)
! SWLTA (triad-wave interactions calculated with the Lumped Triad
! Approximation of Eldeberky, 1996)
! SWSPB (triad-wave interactions calculated with the Stochastic
! Parametric model based on the Boussinesq eqs (SPB) of
! Becq-Girard et al, 1999)
!
!----------------------------------------------------------------------
!
!******************************************************************
!
SUBROUTINE FAC4WW (XIS ,SNLC1 , 40.41 34.00
& DAL1 ,DAL2 ,DAL3 ,SPCSIG, 34.00
& WWINT ,WWAWG ,WWSWG ) 40.17 34.00
!
!******************************************************************
!
USE SWCOMM3 40.41
USE SWCOMM4 40.41
USE OCPCOMM4 40.41
USE M_SNL4 40.17
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: H.L. Tolman, R.C. Ris |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 30.72: IJsbrand Haagsma
! 40.17: IJsbrand Haagsma
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN
! 40.17, Dec. 01: Implementation of Multiple DIA
! 40.41, Sep. 04: compute indices for interactions which will be
! interpolated in piecewise constant manner
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose :
!
! Calculate interpolation constants for Snl.
!
! 3. Method :
!
!
! 4. Argument variables
!
! SPCSIG: Relative frequencies in computational domain in sigma-space 30.72
!
REAL SPCSIG(MSC) 30.72
!
! INTEGERS:
! ---------
! MSC2,MSC1 Auxiliary variables
! MSC,MDC Maximum counters in spectral space
! IDP,IDP1 Positive range for ID
! IDM,IDM1 Negative range for ID
! ISP,ISP1 idem for IS
! ISM,ISM1 idem for IS
! ISCLW,ISCHG Minimum and maximum counter for discrete
! computations in frequency space
! ISLOW,ISHGH Minimum and maximum range in frequency space
! IDLOW,IDHGH idem in directional space
! IS Frequency counter
! MSC4MI,MSC4MA Array dimensions in frequency space
! MDC4MI,MDC4MA Array dimensions in direction space
!
! REALS:
! ------
! LAMBDA Coefficient set 0.25
! GRAV Gravitational acceleration
! SNLC1 Coefficient for the subroutines SWSNLn
! LAMM2,LAMP2
! DELTH3,DELTH4 Angles between the interacting wavenumbers
! DAL1,DAL2,DAL3 Coefficients for the non linear interactions
! CIDP,CIDM
! WIDP,WIDP1,WIDM,WIDM1 Weight factors
! WISP,WISP1,WISM,WISM1 idem
! AWGn Interpolation weight factors
! SWGn Quadratic interpolation weight factors
! XIS,XISLN Difference between succeeding frequencies
! PI 3.14
! FREQ Auxiliary frequency to fill scaling array
! DDIR,RADE band width in directional space and factor 34.00
!
! ARRAYS
! ------
! AF11 1D Scaling frequency
! WWINT 1D counters for 4WAVE interactions
! WWAWG 1D values for the interpolation
! WWSWG 1D vaules for the interpolation
!
! WWINT ( 1 = IDP WWAWG ( = AGW1 WWSWG ( = SWG1
! 2 = IDP1 = AWG2 = SWG2
! 3 = IDM = AWG3 = SWG3
! 4 = IDM1 = AWG4 = SWG4
! 5 = ISP = AWG5 = SWG5
! 6 = ISP1 = AWG6 = SWG6
! 7 = ISM = AWG7 = SWG7
! 8 = ISM1 = AWG8 ) = SWG8 )
! 9 = ISLOW
! 10= ISHGH
! 11= ISCLW
! 12= ISCHG
! 13= IDLOW
! 14= IDHGH
! 15= MSC4MI
! 16= MSC4MA
! 17= MDC4MI
! 18= MDC4MA
! 19= MSCMAX
! 20= MDCMAX
! 21= IDPP
! 22= IDMM
! 23= ISPP
! 24= ISMM )
!
! 7. Common blocks used
!
!
! 9. Source code :
!
! -----------------------------------------------------------------
! Calculate :
! 1. counters for frequency and direction for NL-interaction
! 2. weight factors
! 3. the minimum and maximum counter in IS and ID space
! 4. the interpolation weights
! 5. the quadratic interpolation rates
! 6. fill the array for the frequency**11
! ----------------------------------------------------------
!
!****************************************************************
!
INTEGER MSC2 ,MSC1 ,IS ,IDP ,IDP1 , 40.41 34.00
& IDM ,IDM1 ,ISP ,ISP1 ,ISM ,ISM1 , 34.00
& IDPP ,IDMM ,ISPP ,ISMM , 40.41
& ISLOW ,ISHGH ,ISCLW ,ISCHG ,IDLOW ,IDHGH ,
& MSCMAX,MDCMAX 34.00
!
REAL SNLC1 ,LAMM2 ,LAMP2 ,DELTH3, 40.17 34.00
& AUX1 ,DELTH4,DAL1 ,DAL2 ,DAL3 ,CIDP ,WIDP ,
& WIDP1 ,CIDM ,WIDM ,WIDM1 ,XIS ,XISLN ,WISP ,
& WISP1 ,WISM ,WISM1 ,AWG1 ,AWG2 ,AWG3 ,AWG4 ,
& AWG5 ,AWG6 ,AWG7 ,AWG8 ,SWG1 ,SWG2 ,SWG3 ,
& SWG4 ,SWG5 ,SWG6 ,SWG7 ,SWG8 ,FREQ , 34.00
& RADE 34.00
!
REAL WWAWG(*) , 40.17
& WWSWG(*)
!
INTEGER WWINT(*)
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'FAC4WW')
IF (ALLOCATED(AF11)) DEALLOCATE(AF11) 40.17
! *** Compute frequency indices ***
! *** XIS is the relative increment of the relative frequency ***
!
MSC2 = INT ( FLOAT(MSC) / 2.0 )
MSC1 = MSC2 - 1
XIS = SPCSIG(MSC2) / SPCSIG(MSC1) 30.72
!
! *** set values for the nonlinear four-wave interactions ***
!
SNLC1 = 1. / GRAV**4 40.17 34.00
!
LAMM2 = (1.-PQUAD(1))**2 40.17
LAMP2 = (1.+PQUAD(1))**2 40.17
DELTH3 = ACOS( (LAMM2**2+4.-LAMP2**2) / (4.*LAMM2) )
AUX1 = SIN(DELTH3)
DELTH4 = ASIN(-AUX1*LAMM2/LAMP2)
!
DAL1 = 1. / (1.+PQUAD(1))**4 40.17
DAL2 = 1. / (1.-PQUAD(1))**4 40.17
DAL3 = 2. * DAL1 * DAL2
!
! *** Compute directional indices in sigma and theta space ***
!
CIDP = ABS(DELTH4/DDIR) 40.00
IDP = INT(CIDP)
IDP1 = IDP + 1
WIDP = CIDP - REAL(IDP)
WIDP1 = 1.- WIDP
!
CIDM = ABS(DELTH3/DDIR) 40.00
IDM = INT(CIDM)
IDM1 = IDM + 1
WIDM = CIDM - REAL(IDM)
WIDM1 = 1.- WIDM
XISLN = LOG( XIS )
!
ISP = INT( LOG(1.+PQUAD(1)) / XISLN ) 40.17
ISP1 = ISP + 1
WISP = (1.+PQUAD(1) - XIS**ISP) / (XIS**ISP1 - XIS**ISP) 40.17
WISP1 = 1. - WISP
!
ISM = INT( LOG(1.-PQUAD(1)) / XISLN ) 40.17
ISM1 = ISM - 1
WISM = (XIS**ISM -(1.-PQUAD(1))) / (XIS**ISM - XIS**ISM1) 40.17
WISM1 = 1. - WISM
!
! *** Range of calculations ***
!
ISLOW = 1 + ISM1
ISHGH = MSC + ISP1 - ISM1
ISCLW = 1
ISCHG = MSC - ISM1
IDLOW = 1 - MDC - MAX(IDM1,IDP1)
IDHGH = MDC + MDC + MAX(IDM1,IDP1)
!
MSC4MI = ISLOW
MSC4MA = ISHGH
MDC4MI = IDLOW
MDC4MA = IDHGH
MSCMAX = MSC4MA - MSC4MI + 1
MDCMAX = MDC4MA - MDC4MI + 1
!
! *** Interpolation weights ***
!
AWG1 = WIDP * WISP
AWG2 = WIDP1 * WISP
AWG3 = WIDP * WISP1
AWG4 = WIDP1 * WISP1
!
AWG5 = WIDM * WISM
AWG6 = WIDM1 * WISM
AWG7 = WIDM * WISM1
AWG8 = WIDM1 * WISM1
!
! *** quadratic interpolation ***
!
SWG1 = AWG1**2
SWG2 = AWG2**2
SWG3 = AWG3**2
SWG4 = AWG4**2
!
SWG5 = AWG5**2
SWG6 = AWG6**2
SWG7 = AWG7**2
SWG8 = AWG8**2
!
! --- determine discrete counters for piecewise 40.41
! constant interpolation 40.41
!
IF (AWG1.LT.AWG2) THEN
IF (AWG2.LT.AWG3) THEN
IF (AWG3.LT.AWG4) THEN
ISPP=ISP
IDPP=IDP
ELSE
ISPP=ISP
IDPP=IDP1
END IF
ELSE IF (AWG2.LT.AWG4) THEN
ISPP=ISP
IDPP=IDP
ELSE
ISPP=ISP1
IDPP=IDP
END IF
ELSE IF (AWG1.LT.AWG3) THEN
IF (AWG3.LT.AWG4) THEN
ISPP=ISP
IDPP=IDP
ELSE
ISPP=ISP
IDPP=IDP1
END IF
ELSE IF (AWG1.LT.AWG4) THEN
ISPP=ISP
IDPP=IDP
ELSE
ISPP=ISP1
IDPP=IDP1
END IF
IF (AWG5.LT.AWG6) THEN
IF (AWG6.LT.AWG7) THEN
IF (AWG7.LT.AWG8) THEN
ISMM=ISM
IDMM=IDM
ELSE
ISMM=ISM
IDMM=IDM1
END IF
ELSE IF (AWG6.LT.AWG8) THEN
ISMM=ISM
IDMM=IDM
ELSE
ISMM=ISM1
IDMM=IDM
END IF
ELSE IF (AWG5.LT.AWG7) THEN
IF (AWG7.LT.AWG8) THEN
ISMM=ISM
IDMM=IDM
ELSE
ISMM=ISM
IDMM=IDM1
END IF
ELSE IF (AWG5.LT.AWG8) THEN
ISMM=ISM
IDMM=IDM
ELSE
ISMM=ISM1
IDMM=IDM1
END IF
!
! *** fill the arrays *
!
WWINT(1) = IDP
WWINT(2) = IDP1
WWINT(3) = IDM
WWINT(4) = IDM1
WWINT(5) = ISP
WWINT(6) = ISP1
WWINT(7) = ISM
WWINT(8) = ISM1
WWINT(9) = ISLOW
WWINT(10)= ISHGH
WWINT(11)= ISCLW
WWINT(12)= ISCHG
WWINT(13)= IDLOW
WWINT(14)= IDHGH
WWINT(15)= MSC4MI
WWINT(16)= MSC4MA
WWINT(17)= MDC4MI
WWINT(18)= MDC4MA
WWINT(19)= MSCMAX
WWINT(20)= MDCMAX
WWINT(21)= IDPP 40.41
WWINT(22)= IDMM 40.41
WWINT(23)= ISPP 40.41
WWINT(24)= ISMM 40.41
!
WWAWG(1) = AWG1
WWAWG(2) = AWG2
WWAWG(3) = AWG3
WWAWG(4) = AWG4
WWAWG(5) = AWG5
WWAWG(6) = AWG6
WWAWG(7) = AWG7
WWAWG(8) = AWG8
!
WWSWG(1) = SWG1
WWSWG(2) = SWG2
WWSWG(3) = SWG3
WWSWG(4) = SWG4
WWSWG(5) = SWG5
WWSWG(6) = SWG6
WWSWG(7) = SWG7
WWSWG(8) = SWG8
ALLOCATE (AF11(MSC4MI:MSC4MA)) 40.17
! *** Fill scaling array (f**11) ***
! *** compute the radian frequency**11 for IS=1, MSC ***
!
DO 100 IS=1, MSC
AF11(IS) = ( SPCSIG(IS) / ( 2. * PI ) )**11 30.72
100 CONTINUE
!
! *** compute the radian frequency for the IS = MSC+1, ISHGH ***
!
FREQ = SPCSIG(MSC) / ( 2. * PI ) 30.72
DO 110 IS = MSC+1, ISHGH
FREQ = FREQ * XIS
AF11(IS) = FREQ**11
110 CONTINUE
!
! *** compute the radian frequency for IS = 0, ISLOW ***
!
FREQ = SPCSIG(1) / ( 2. * PI ) 30.72
DO 120 IS = 0, ISLOW, -1
FREQ = FREQ / XIS
AF11(IS) = FREQ**11
120 CONTINUE
!
! *** test output ***
!
IF (ISLOW .LT. MSC4MI .OR. ISHGH .GT. MSC4MA .OR.
& IDLOW .LT. MDC4MI .OR. IDHGH .GT. MDC4MA) THEN
WRITE (PRINTF,900) IXCGRD(1), IYCGRD(1),
& ISLOW, ISHGH, IDLOW, IDHGH,
& MSC4MI,MSC4MA, MDC4MI, MDC4MA
900 FORMAT ( ' ** Error : array bounds and maxima in subr FAC4WW, ',
& ' point ', 2I5,
& /,' ISL,ISH : ',2I4, ' IDL,IDH : ',2I4,
& /,' SMI,SMA : ',2I4, ' DMI,DMA : ',2I4)
ENDIF
!
IF (ITEST .GE. 40) THEN
RADE = 360.0 / ( 2. * PI )
WRITE(PRINTF,*)
WRITE(PRINTF,*) ' FAC4WW subroutine '
WRITE(PRINTF,9000) DELTH4*RADE, DELTH3*RADE, DDIR*RADE, XIS
9000 FORMAT (' THET3 THET4 DDIR XIS :',4E12.4)
WRITE(PRINTF,9011) IDP, IDP1, IDM, IDM1
9011 FORMAT (' IDP IDP1 IDM IDM1 :',4I5)
WRITE(PRINTF,9012) WIDP, WIDP1, WIDM, WIDM1
9012 FORMAT (' WIDP WIDP1 WIDM WIDM1 :',4E12.4)
WRITE (PRINTF,9013) ISP, ISP1, ISM, ISM1
9013 FORMAT (' ISP ISP1 ISM ISM1 :',4I5)
WRITE (PRINTF,9014) WISP, WISP1, WISM, WISM1
9014 FORMAT (' WISP WISP1 WISM WISM1 :',4E12.4)
WRITE(PRINTF,9016) ISCLW, ISCHG
9016 FORMAT (' ICLW ICHG :',2I5)
WRITE (PRINTF,9017) AWG1, AWG2, AWG3, AWG4
9017 FORMAT (' AWG1 AWG2 AWG3 AWG4 :',4E12.4)
WRITE (PRINTF,9018) AWG5, AWG6, AWG7, AWG8
9018 FORMAT (' AWG5 AWG6 AWG7 AWG8 :',4E12.4)
WRITE (PRINTF,9019) MSC4MI, MSC4MA, MDC4MI, MDC4MA
9019 FORMAT (' S4MI S4MA D4MI D4MA :',4I6)
WRITE (PRINTF,9015) ISLOW, ISHGH, IDLOW,IDHGH
9015 FORMAT (' ISLOW ISHG IDLOW IDHG :',4I5)
WRITE(PRINTF,*)
END IF
!
RETURN
! End of FAC4WW
END
!
!******************************************************************
!
SUBROUTINE RANGE4 (WWINT ,IDDLOW,IDDTOP) 40.00
!
!******************************************************************
!
USE SWCOMM3 40.41
USE SWCOMM4 40.41
USE OCPCOMM4 40.41
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 40.00: Nico Booij
! 40.10: IJsbrand Haagsma
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 40.10, Mar. 00: Made modification for exact quadruplets
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose :
!
! calculate the minimum and maximum counters in frequency and
! directional space which fall with the calculation for the
! nonlinear wave-wave interactions.
!
! 3. Method : review for the counters :
!
! Frequencies -->
! +---+---------------------+---------+- IDHGH
! d | 3 : 2 : 2 |
! i + - + - - - - - - - - - - + - - - - +- MDC
! r | : : |
! e | 3 : original spectrum : 1 |
! c | : : |
! t. + - + - - - - - - - - - - + - - - - +- 1
! | 3 : 2 : 2 |
! +---+---------------------+---------+- IDLOW
! | | | ^ |
! ISLOW 1 MSC | ISHGH
! ^ |
! | |
! ISCLW ISCHG
! lowest discrete highest discrete
! central bin central bin
!
!
! The directional counters depend on the numerical method that
! is used.
!
! 4. Parameters :
!
! INTEGER
! -------
! IQUAD Counter for 4 wave interactions
! ISLOW,ISHGH Minimum and maximum counter in frequency space
! ISCLW,ISCHG idem for discrete computations
! IDLOW,IDHGH Minimum and maximum counters in directional space
! MSC,MDC Range of the original arrays
! ISM1,ISP1,
! IDM1,IDP1 see subroutine FAC4WW
! IDDLOW minimum counter of the bin that is propagated
! within a sweep
! IDDTOP minimum counter of the bin that is propagated
! within a sweep
!
! array:
! ------
! WWINT counters for the nonlinear interactions
!
! WWINT ( 1 = IDP 2 = IDP1 3 = IDM 4 = IDM1
! 5 = ISP 6 = ISP1 7 = ISM 8 = ISM1
! 9 = ISLOW 10 = ISHGH 11 = ISCLW 12 = ISCHG
! 13 = IDLOW 14 = IDHGH 15 = MSC4MI 16 = MSC4MA
! 17 = MDC4MI 18 = MDC4MA
! 19 = MSCMAX 20 = MDCMAX )
!
! 5. Subroutines used :
!
! ---
!
! 6. Called by :
!
! SOURCE
!
! 7. Common blocks used
!
!
! 9. Source code :
!
! -----------------------------------------------------------------
! Calculate :
! In absence of a current there are always four sectors
! equal 90 degrees within a sweep that are propagated
! Extend the boundaries to calculate the source term
! In presence of a current and if IDTOT .eq. MDC then calculate
! boundaries for calculation of interaction using the
! unfolded area.
! ----------------------------------------------------------
!
!****************************************************************
!
INTEGER IDDLOW,IDDTOP 40.00
!
INTEGER WWINT(*)
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'RANGE4')
!
! *** Range in directional domain ***
!
IF ( IQUAD .LT. 3 .AND. IQUAD .GT. 0 ) THEN 40.10
! *** counters based on bins which fall within a sweep ***
WWINT(13) = IDDLOW - MAX( WWINT(4), WWINT(2) )
WWINT(14) = IDDTOP + MAX( WWINT(4), WWINT(2) )
ELSE
! *** counters initially based on full circle ***
WWINT(13) = 1 - MAX( WWINT(4), WWINT(2) )
WWINT(14) = MDC + MAX( WWINT(4), WWINT(2) )
END IF
!
! *** error message ***
!
IF (WWINT(9) .LT. WWINT(15) .OR. WWINT(10) .GT. WWINT(16) .OR.
& WWINT(13) .LT. WWINT(17) .OR. WWINT(14) .GT. WWINT(18) ) THEN
WRITE (PRINTF,900) IXCGRD(1), IYCGRD(1),
& WWINT(9) ,WWINT(10) ,WWINT(13) ,WWINT(14),
& WWINT(15),WWINT(16) ,WWINT(17) ,WWINT(18)
900 FORMAT ( ' ** Error : array bounds and maxima in subr RANGE4, ',
& ' point ', 2I5,
& /,' ISL,ISH : ',2I4, ' IDL,IDH : ',2I4,
& /,' SMI,SMA : ',2I4, ' DMI,DMA : ',2I4)
IF (ITEST.GE.50) WRITE (PRTEST, 901) MSC, MDC, IDDLOW, IDDTOP
901 FORMAT (' MSC, MDC, IDDLOW, IDDTOP: ', 4I5)
ENDIF
!
! test output
!
IF (TESTFL .AND. ITEST .GE. 60) THEN
WRITE(PRTEST,911) WWINT(4), WWINT(2), WWINT(8), WWINT(6)
911 FORMAT (' RANGE4: IDM1 IDP1 ISM1 ISP1 :',4I5)
WRITE(PRTEST,916) WWINT(11), WWINT(12), IQUAD
916 FORMAT (' RANGE4: ISCLW ISCHG IQUAD :',3I5)
WRITE (PRTEST,917) WWINT(9), WWINT(10), WWINT(13), WWINT(14)
917 FORMAT (' RANGE4: ISLOW ISHGH IDLOW IDHGH:',4I5)
WRITE (PRTEST,919) WWINT(15), WWINT(16), WWINT(17), WWINT(18)
919 FORMAT (' RANGE4: MS4MI MS4MA MD4MI MD4MA:',4I5)
WRITE(PRINTF,*)
END IF
!
RETURN
! End of RANGE4
END
!
!********************************************************************
!
SUBROUTINE SWSNL1 (WWINT ,WWAWG ,WWSWG , 34.00
& IDCMIN ,IDCMAX ,UE ,SA1 , 40.17
& SA2 ,DA1C ,DA1P ,DA1M ,DA2C ,
& DA2P ,DA2M ,SPCSIG ,SNLC1 ,KMESPC , 30.72
& FACHFR ,ISSTOP ,DAL1 ,DAL2 ,DAL3 ,
& SFNL ,DSNL ,DEP2 ,AC2 ,IMATDA ,
& IMATRA ,PLNL4S ,PLNL4D , 34.00
& IDDLOW ,IDDTOP ,REDC0 ,REDC1 ) 40.85 34.00
!
!********************************************************************
!
USE SWCOMM3 40.41
USE SWCOMM4 40.41
USE OCPCOMM4 40.41
USE M_SNL4 40.17
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: H.L. Tolman, R.C. Ris |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 30.72: IJsbrand Haagsma
! 40.13: Nico Booij
! 40.17: IJsbrand Haagsma
! 40.23: Marcel Zijlema
! 40.41: Marcel Zijlema
! 40.85: Marcel Zijlema
!
! 1. Updates
!
! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN
! 40.17, Dec. 01: Implentation of Multiple DIA
! 40.23, Aug. 02: some corrections
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
! 40.85, Aug. 08: store quadruplets for output purposes
!
! 2. Purpose
!
! Calculate non-linear interaction using the discrete interaction
! approximation (Hasselmann and Hasselmann 1985; WAMDI group 1988),
! including the diagonal term for the implicit integration.
!
! The interactions are calculated for all bin's that fall
! within a sweep. No additional auxiliary array is required (see
! SWSNL3)
!
! 3. Method
!
! Discrete interaction approximation.
!
! Since the domain in directional domain is by definition not
! periodic, the spectral space can not beforehand
! folded to the side angles. This can only be done if the
! full circle has to be calculated
!
!
! Frequencies -->
! +---+---------------------+---------+- IDHGH
! d | 3 : 2 : 2 |
! i + - + - - - - - - - - - - + - - - - +- MDC
! r | : : |
! e | 3 : original spectrum : 1 |
! c | : : |
! t. + - + - - - - - - - - - - + - - - - +- 1
! | 3 : 2 : 2 |
! +---+---------------------+---------+- IDLOW
! | | | ^ |
! ISLOW 1 MSC | ISHGH
! ^ |
! | |
! ISCLW ISCHG
! lowest discrete highest discrete
! central bin central bin
!
! 1 : Extra tail added beyond MSC
! 2 : Spectrum copied outside ID range
! 3 : Empty bins at low frequencies
!
! ISLOW = 1 + ISM1
! ISHGH = MSC + ISP1 - ISM1
! ISCLW = 1
! ISCHG = MSC - ISM1
! IDLOW = IDDLOW - MAX(IDM1,IDP1)
! IDHGH = IDDTOP + MAX(IDM1,IDP1)
!
! For the meaning of the counters on the right hand side of the
! above equations see section 4.
!
! 4. Argument variables
!
! SPCSIG: Relative frequencies in computational domain in sigma-space 30.72
!
REAL SPCSIG(MSC) 30.72
!
! Data in PARAMETER statements :
! ----------------------------------------------------------------
! DAL1 Real LAMBDA dependend weight factors (see FAC4WW)
! DAL2 Real
! DAL3 Real
! ITHP, ITHP1, ITHM, ITHM1, IFRP, IFRP1, IFRM, IFRM1
! Int. Counters of interpolation point relative to
! central bin, see figure below (set in FAC4WW).
! NFRLOW, NFRHGH, NFRCHG, NTHLOW, NTHHGH
! Int. Range of calculations, see section 2.
! AF11 R.A. Scaling array (Freq**11).
! AWGn Real Interpolation weights, see numbers in fig.
! SWGn Real Id. squared.
! UE R.A. "Unfolded" spectrum.
! SA1 R.A. Interaction constribution of first and second
! SA2 R.A. quadr. respectively (unfolded space).
! DA1C, DA1P, DA1M, DA2C, DA2P, DA2M
! R.A. Idem for diagonal matrix.
! PERCIR full circle or sector
! ----------------------------------------------------------------
!
! Relative offsets of interpolation points around central bin
! "#" and corresponding numbers of AWGn :
!
! ISM1 ISM
! 5 7 T |
! IDM1 +------+ H +
! | | E | ISP ISP1
! | \ | T | 3 1
! IDM +------+ A + +---------+ IDP1
! 6 \8 | | |
! | | / |
! \ + +---------+ IDP
! | /4 2
! \ | /
! -+-----+------+-------#--------+---------+----------+
! | FREQ.
!
! 7. Common blocks used
!
!
! 8. Subroutines used
!
! ---
!
! 9. Subroutines calling
!
! SOURCE (in SWANCOM1)
!
! 12. Structure
!
! -------------------------------------------
! Initialisations.
! Calculate proportionality constant.
! Prepare auxiliary spectrum.
! Calculate interactions :
! -----------------------------------------
! Energy at interacting bins
! Contribution to interactions
! Fold interactions to side angles
! -----------------------------------------
! Put source term together
! -------------------------------------------
!
! 13. Source text
!
!*************************************************************
!
INTEGER IS ,ID ,I ,J , 34.00
& ISHGH ,IDLOW ,ISP ,ISP1 ,IDP ,IDP1 ,
& ISM ,ISM1 ,IDHGH ,IDM ,IDM1 ,ISCLW ,
& ISCHG ,IDDLOW ,IDDTOP 34.00
!
REAL X ,X2 ,CONS ,FACTOR ,SNLCS1 ,SNLCS2 ,SNLCS3,
& E00 ,EP1 ,EM1 ,EP2 ,EM2 ,SA1A ,SA1B ,
& SA2A ,SA2B ,KMESPC ,FACHFR ,AWG1 ,AWG2 ,AWG3 ,
& AWG4 ,AWG5 ,AWG6 ,AWG7 ,AWG8 ,DAL1 ,DAL2 ,
& DAL3 ,SNLC1 ,SWG1 ,SWG2 ,SWG3 ,SWG4 ,SWG5 ,
& SWG6 ,SWG7 ,SWG8 ,JACOBI ,SIGPI 34.00
!
REAL AC2(MDC,MSC,MCGRD) ,
& DEP2(MCGRD) ,
& UE(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) ,
& SA1(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) ,
& SA2(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) ,
& DA1C(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) ,
& DA1P(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) ,
& DA1M(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) ,
& DA2C(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) ,
& DA2P(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) ,
& DA2M(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) ,
& SFNL(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) ,
& DSNL(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) ,
& IMATDA(MDC,MSC) ,
& IMATRA(MDC,MSC) ,
& PLNL4S(MDC,MSC,NPTST) , 40.00
& PLNL4D(MDC,MSC,NPTST) ,
& WWAWG(*) ,
& WWSWG(*)
REAL :: REDC0 (MDC,MSC,MREDS) 40.85
REAL :: REDC1 (MDC,MSC,MREDS) 40.85
!
INTEGER IDCMIN(MSC) ,
& IDCMAX(MSC) ,
& WWINT(*)
!
LOGICAL PERCIR
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'SWSNL1')
!
IDP = WWINT(1)
IDP1 = WWINT(2)
IDM = WWINT(3)
IDM1 = WWINT(4)
ISP = WWINT(5)
ISP1 = WWINT(6)
ISM = WWINT(7)
ISM1 = WWINT(8)
ISLOW = WWINT(9)
ISHGH = WWINT(10)
ISCLW = WWINT(11)
ISCHG = WWINT(12)
IDLOW = WWINT(13)
IDHGH = WWINT(14)
!
AWG1 = WWAWG(1)
AWG2 = WWAWG(2)
AWG3 = WWAWG(3)
AWG4 = WWAWG(4)
AWG5 = WWAWG(5)
AWG6 = WWAWG(6)
AWG7 = WWAWG(7)
AWG8 = WWAWG(8)
!
SWG1 = WWSWG(1)
SWG2 = WWSWG(2)
SWG3 = WWSWG(3)
SWG4 = WWSWG(4)
SWG5 = WWSWG(5)
SWG6 = WWSWG(6)
SWG7 = WWSWG(7)
SWG8 = WWSWG(8)
!
! *** Initialize auxiliary arrays per gridpoint ***
!
DO ID = MDC4MI, MDC4MA
DO IS = MSC4MI, MSC4MA
UE(IS,ID) = 0.
SA1(IS,ID) = 0.
SA2(IS,ID) = 0.
SFNL(IS,ID) = 0.
DA1C(IS,ID) = 0.
DA1P(IS,ID) = 0.
DA1M(IS,ID) = 0.
DA2C(IS,ID) = 0.
DA2P(IS,ID) = 0.
DA2M(IS,ID) = 0.
DSNL(IS,ID) = 0.
ENDDO
ENDDO
!
! *** Calculate factor R(X) to calculate the NL wave-wave ***
! *** interaction for shallow water ***
! *** SNLC1 = 1/GRAV**4 *** 40.17
!
SNLCS1 = PQUAD(3) 34.00
SNLCS2 = PQUAD(4) 34.00
SNLCS3 = PQUAD(5) 34.00
X = MAX ( 0.75 * DEP2(KCGRD(1)) * KMESPC , 0.5 )
X2 = MAX ( -1.E15, SNLCS3*X)
CONS = SNLC1 * ( 1. + SNLCS1/X * (1.-SNLCS2*X) * EXP(X2))
JACOBI = 2. * PI
!
! *** check whether the spectral domain is periodic in ***
! *** directional space and if so, modify boundaries ***
!
PERCIR = .FALSE.
IF ( IDDLOW .EQ. 1 .AND. IDDTOP .EQ. MDC ) THEN
! *** periodic in theta -> spectrum can be folded ***
! *** (can only be present in presence of a current) ***
IDCLOW = 1
IDCHGH = MDC
IIID = 0
PERCIR = .TRUE.
ELSE
! *** different sectors per sweep -> extend range with IIID ***
IIID = MAX ( IDM1 , IDP1 )
IDCLOW = IDLOW
IDCHGH = IDHGH
ENDIF
!
! *** Prepare auxiliary spectrum ***
! *** set action original spectrum in array UE ***
!
DO IDDUM = IDLOW - IIID, IDHGH + IIID
ID = MOD ( IDDUM - 1 + MDC , MDC ) + 1
DO IS = 1, MSC
UE(IS,IDDUM) = AC2(ID,IS,KCGRD(1)) * SPCSIG(IS) * JACOBI 30.72
ENDDO
ENDDO
!
! *** set values in area 2 for IS > MSC+1 ***
!
DO IS = MSC+1, ISHGH
DO ID = IDLOW - IIID , IDHGH + IIID
UE (IS,ID) = UE(IS-1,ID) * FACHFR
ENDDO
ENDDO
!
! *** Calculate interactions ***
! *** Energy at interacting bins ***
!
DO IS = ISCLW, ISCHG
DO ID = IDCLOW, IDCHGH
E00 = UE(IS ,ID )
EP1 = AWG1 * UE(IS+ISP1,ID+IDP1) +
& AWG2 * UE(IS+ISP1,ID+IDP ) +
& AWG3 * UE(IS+ISP ,ID+IDP1) +
& AWG4 * UE(IS+ISP ,ID+IDP )
EM1 = AWG5 * UE(IS+ISM1,ID-IDM1) +
& AWG6 * UE(IS+ISM1,ID-IDM ) +
& AWG7 * UE(IS+ISM ,ID-IDM1) +
& AWG8 * UE(IS+ISM ,ID-IDM )
!
EP2 = AWG1 * UE(IS+ISP1,ID-IDP1) +
& AWG2 * UE(IS+ISP1,ID-IDP ) +
& AWG3 * UE(IS+ISP ,ID-IDP1) +
& AWG4 * UE(IS+ISP ,ID-IDP )
EM2 = AWG5 * UE(IS+ISM1,ID+IDM1) +
& AWG6 * UE(IS+ISM1,ID+IDM ) +
& AWG7 * UE(IS+ISM ,ID+IDM1) +
& AWG8 * UE(IS+ISM ,ID+IDM )
!
! *** Contribution to interactions ***
! *** CONS is the shallow water factor for the NL interact. ***
!
FACTOR = CONS * AF11(IS) * E00
!
SA1A = E00 * ( EP1*DAL1 + EM1*DAL2 ) * PQUAD(2) 40.17
SA1B = SA1A - EP1*EM1*DAL3 * PQUAD(2) 40.17
SA2A = E00 * ( EP2*DAL1 + EM2*DAL2 ) * PQUAD(2) 40.17
SA2B = SA2A - EP2*EM2*DAL3 * PQUAD(2) 40.17
!
SA1 (IS,ID) = FACTOR * SA1B
SA2 (IS,ID) = FACTOR * SA2B
!
IF(ITEST.GE.100 .AND. TESTFL) THEN
WRITE(PRINTF,9002) E00,EP1,EM1,EP2,EM2
9002 FORMAT (' E00 EP1 EM1 EP2 EM2 :',5E11.4)
WRITE(PRINTF,9003) SA1A,SA1B,SA2A,SA2B
9003 FORMAT (' SA1A SA1B SA2A SA2B :',4E11.4)
WRITE(PRINTF,9004) IS,ID,SA1(IS,ID),SA2(IS,ID)
9004 FORMAT (' IS ID SA1() SA2() :',2I4,2E12.4)
WRITE(PRINTF,9005) FACTOR
9005 FORMAT (' FACTOR : ',E12.4)
END IF
!
DA1C(IS,ID) = CONS * AF11(IS) * ( SA1A + SA1B )
DA1P(IS,ID) = FACTOR * ( DAL1*E00 - DAL3*EM1 ) * PQUAD(2) 40.23
DA1M(IS,ID) = FACTOR * ( DAL2*E00 - DAL3*EP1 ) * PQUAD(2) 40.23
!
DA2C(IS,ID) = CONS * AF11(IS) * ( SA2A + SA2B )
DA2P(IS,ID) = FACTOR * ( DAL1*E00 - DAL3*EM2 ) * PQUAD(2) 40.23
DA2M(IS,ID) = FACTOR * ( DAL2*E00 - DAL3*EP2 ) * PQUAD(2) 40.23
ENDDO
ENDDO
!
! *** Fold interactions to side angles if spectral domain ***
! *** is periodic in directional space ***
!
IF ( PERCIR ) THEN
DO ID = 1, IDHGH - MDC
ID0 = 1 - ID
DO IS = ISCLW, ISCHG
SA1 (IS,MDC+ID) = SA1 (IS, ID )
SA2 (IS,MDC+ID) = SA2 (IS, ID )
DA1C(IS,MDC+ID) = DA1C(IS, ID )
DA1P(IS,MDC+ID) = DA1P(IS, ID )
DA1M(IS,MDC+ID) = DA1M(IS, ID )
DA2C(IS,MDC+ID) = DA2C(IS, ID )
DA2P(IS,MDC+ID) = DA2P(IS, ID )
DA2M(IS,MDC+ID) = DA2M(IS, ID )
!
SA1 (IS, ID0 ) = SA1 (IS, MDC+ID0)
SA2 (IS, ID0 ) = SA2 (IS, MDC+ID0)
DA1C(IS, ID0 ) = DA1C(IS, MDC+ID0)
DA1P(IS, ID0 ) = DA1P(IS, MDC+ID0)
DA1M(IS, ID0 ) = DA1M(IS, MDC+ID0)
DA2C(IS, ID0 ) = DA2C(IS, MDC+ID0)
DA2P(IS, ID0 ) = DA2P(IS, MDC+ID0)
DA2M(IS, ID0 ) = DA2M(IS, MDC+ID0)
ENDDO
ENDDO
ENDIF
!
! *** Put source term together (To save space I=IS and J=ID ***
! *** is used) ***
!
PI3 = (2. * PI)**3
DO I = 1, ISSTOP
SIGPI = SPCSIG(I) * JACOBI 30.72
DO J = IDCMIN(I), IDCMAX(I)
ID = MOD ( J - 1 + MDC , MDC ) + 1
SFNL(I,ID) = - 2. * ( SA1(I,J) + SA2(I,J) )
& + AWG1 * ( SA1(I-ISP1,J-IDP1) + SA2(I-ISP1,J+IDP1) )
& + AWG2 * ( SA1(I-ISP1,J-IDP ) + SA2(I-ISP1,J+IDP ) )
& + AWG3 * ( SA1(I-ISP ,J-IDP1) + SA2(I-ISP ,J+IDP1) )
& + AWG4 * ( SA1(I-ISP ,J-IDP ) + SA2(I-ISP ,J+IDP ) )
& + AWG5 * ( SA1(I-ISM1,J+IDM1) + SA2(I-ISM1,J-IDM1) )
& + AWG6 * ( SA1(I-ISM1,J+IDM ) + SA2(I-ISM1,J-IDM ) )
& + AWG7 * ( SA1(I-ISM ,J+IDM1) + SA2(I-ISM ,J-IDM1) )
& + AWG8 * ( SA1(I-ISM ,J+IDM ) + SA2(I-ISM ,J-IDM ) )
!
DSNL(I,ID) = - 2. * ( DA1C(I,J) + DA2C(I,J) )
& + SWG1 * ( DA1P(I-ISP1,J-IDP1) + DA2P(I-ISP1,J+IDP1) )
& + SWG2 * ( DA1P(I-ISP1,J-IDP ) + DA2P(I-ISP1,J+IDP ) )
& + SWG3 * ( DA1P(I-ISP ,J-IDP1) + DA2P(I-ISP ,J+IDP1) )
& + SWG4 * ( DA1P(I-ISP ,J-IDP ) + DA2P(I-ISP ,J+IDP ) )
& + SWG5 * ( DA1M(I-ISM1,J+IDM1) + DA2M(I-ISM1,J-IDM1) )
& + SWG6 * ( DA1M(I-ISM1,J+IDM ) + DA2M(I-ISM1,J-IDM ) )
& + SWG7 * ( DA1M(I-ISM ,J+IDM1) + DA2M(I-ISM ,J-IDM1) )
& + SWG8 * ( DA1M(I-ISM ,J+IDM ) + DA2M(I-ISM ,J-IDM ) )
!
! *** store results in IMATDA and IMATRA ***
!
IF(TESTFL) THEN
PLNL4S(ID,I,IPTST) = SFNL(I,ID) / SIGPI 40.00
PLNL4D(ID,I,IPTST) = DSNL(I,ID) / PI3 40.85 40.00
END IF
REDC0(ID,I,1) = REDC0(ID,I,1) + SFNL(I,ID) / SIGPI 40.85
REDC1(ID,I,1) = REDC1(ID,I,1) + DSNL(I,ID) / PI3 40.85
!
IMATRA(ID,I) = IMATRA(ID,I) + SFNL(I,ID) / SIGPI
IMATDA(ID,I) = IMATDA(ID,I) - DSNL(I,ID) / PI3
!
IF(ITEST.GE.90 .AND. TESTFL) THEN
WRITE(PRINTF,9006) I,J,SFNL(I,ID),DSNL(I,ID),
& SPCSIG(I) 30.72
9006 FORMAT (' IS ID SFNL DSNL SPCSIG:',2I4,3E12.4) 30.72
END IF
!
ENDDO
ENDDO
!
! *** test output ***
!
IF (ITEST .GE. 50 .AND. TESTFL) THEN
WRITE(PRINTF,*)
WRITE(PRINTF,*) ' SWSNL1 subroutine '
WRITE(PRINTF,9011) IDP, IDP1, IDM, IDM1
9011 FORMAT (' IDP IDP1 IDM IDM1 :',4I5)
WRITE (PRINTF,9013) ISP, ISP1, ISM, ISM1
9013 FORMAT (' ISP ISP1 ISM ISM1 :',4I5)
WRITE (PRINTF,9015) ISLOW, ISHGH, IDLOW,IDHGH
9015 FORMAT (' ISLOW ISHGH IDLOW IDHG:',4I5)
WRITE(PRINTF,9016) ISCLW, ISCHG, IDDLOW, IDDTOP
9016 FORMAT (' ICLW ICHG IDDLOW IDDTO:',2I5)
WRITE (PRINTF,9017) AWG1, AWG2, AWG3, AWG4
9017 FORMAT (' AWG1 AWG2 AWG3 AWG4 :',4E12.4)
WRITE (PRINTF,9018) AWG5, AWG6, AWG7, AWG8
9018 FORMAT (' AWG5 AWG6 AWG7 AWG8 :',4E12.4)
WRITE (PRINTF,9019) MSC4MI, MSC4MA, MDC4MI, MDC4MA
9019 FORMAT (' S4MI S4MA D4MI D4MA :',4I6)
WRITE(PRINTF,9020) SNLC1,X,X2,CONS
9020 FORMAT (' SNLC1 X X2 CONS :',4E12.4)
WRITE(PRINTF,9021) DEP2(KCGRD(1)),KMESPC, FACHFR, PI
9021 FORMAT (' DEPTH KMESPC FACHFR PI:',4E12.4)
WRITE(PRINTF,9023) JACOBI
9023 FORMAT (' JACOBI :',E12.4)
WRITE(PRINTF,*)
END IF
!
RETURN
! End of the subroutine SWSNL1
END
!
!*******************************************************************
!
SUBROUTINE SWSNL2 (IDDLOW ,IDDTOP ,WWINT , 34.00
& WWAWG ,UE ,SA1 ,ISSTOP , 40.17
& SA2 ,SPCSIG ,SNLC1 ,DAL1 ,DAL2 , 30.72
& DAL3 ,SFNL ,DEP2 ,AC2 ,KMESPC ,
& REDC0 ,REDC1 ,IMATDA ,IMATRA , 40.85 40.23 34.00
& FACHFR ,PLNL4S , IDCMIN ,IDCMAX ) 34.00
!
!*******************************************************************
!
USE SWCOMM3 40.41
USE SWCOMM4 40.41
USE OCPCOMM4 40.41
USE M_SNL4 40.17
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: H.L. Tolman, R.C. Ris |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 30.72: IJsbrand Haagsma
! 40.13: Nico Booij
! 40.17: IJsbrand Haagsma
! 40.23: Marcel Zijlema
! 40.41: Marcel Zijlema
! 40.85: Marcel Zijlema
!
! 1. Updates
!
! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN
! 40.17, Dec. 01: Implemented Multiple DIA
! 40.23, Aug. 02: rhs and main diagonal adjusted according to Patankar-rules
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
! 40.85, Aug. 08: store quadruplets for output purposes
!
! 2. Purpose
!
! Calculate non-linear interaction using the discrete interaction
! approximation (Hasselmann and Hasselmann 1985; WAMDI group 1988)
!
! 3. Method
!
! Discrete interaction approximation.
!
! Frequencies -->
! +---+---------------------+---------+- IDHGH
! d | 3 : 2 : 2 |
! i + - + - - - - - - - - - - + - - - - +- MDC
! r | : : |
! e | 3 : original spectrum : 1 |
! c | : : |
! t. + - + - - - - - - - - - - + - - - - +- 1
! | 3 : 2 : 2 |
! +---+---------------------+---------+- IDLOW
! | | | ^ |
! ISLOW 1 MSC | ISHGH
! | |
! ISCLW ISCHG
! lowest discrete highest discrete
! central bin central bin
!
! 1 : Extra tail added beyond MSC
! 2 : Spectrum copied outside ID range
! 3 : Empty bins at low frequencies
!
! ISLOW = 1 + ISM1
! ISHGH = MSC + ISP1 - ISM1
! ISCLW = 1
! ISCHG = MSC - ISM1
! IDLOW = IDDLOW - MAX(IDM1,IDP1)
! IDHGH = IDDTOP + MAX(IDM1,IDP1)
!
! Relative offsets of interpolation points around central bin
! "#" and corresponding numbers of AWGn :
!
! ISM1 ISM
! 5 7 T |
! IDM1 +------+ H +
! | | E | ISP ISP1
! | \ | T | 3 1
! IDM +------+ A + +---------+ IDP1
! 6 \8 | | |
! | | / |
! \ + +---------+ IDP
! | /4 2
! \ | /
! -+-----+------+-------#--------+---------+----------+
! | FREQ.
!
!
! 4. Argument variables
!
! SPCSIG: Relative frequencies in computational domain in sigma-space 30.72
!
REAL SPCSIG(MSC) 30.72
!
! 7. Common blocks used
!
!
! 8. Subroutines used
!
! ---
!
! 9. Subroutines calling
!
! SOURCE (in SWANCOM1)
!
! 12. Structure
!
! -------------------------------------------
! Initialisations.
! Calculate proportionality constant.
! Prepare auxiliary spectrum.
! Calculate (unfolded) interactions :
! -----------------------------------------
! Energy at interacting bins
! Contribution to interactions
! Fold interactions to side angles
! -----------------------------------------
! Put source term together
! -------------------------------------------
!
! 13. Source text
!
!*******************************************************************
!
INTEGER IS ,ID ,I ,J ,ISHGH , 34.00
& ISSTOP ,ISP ,ISP1 ,IDP ,IDP1 ,ISM ,ISM1 ,
& IDM ,IDM1 ,ISCLW ,ISCHG , 34.00
& IDLOW ,IDHGH ,IDDLOW ,IDDTOP ,IDCLOW ,IDCHGH 34.00
!
REAL X ,X2 ,CONS ,FACTOR ,SNLCS1 ,SNLCS2 ,SNLCS3 ,
& E00 ,EP1 ,EM1 ,EP2 ,EM2 ,SA1A ,SA1B ,
& SA2A ,SA2B ,KMESPC ,FACHFR ,AWG1 ,AWG2 ,AWG3 ,
& AWG4 ,AWG5 ,AWG6 ,AWG7 ,AWG8 ,DAL1 ,DAL2 ,
& DAL3 ,JACOBI ,SIGPI 34.00
!
REAL AC2(MDC,MSC,MCGRD) , 30.21
& DEP2(MCGRD) , 30.21
& UE(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) ,
& SA1(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) ,
& SA2(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) ,
& SFNL(MSC4MI:MSC4MA , MDC4MI:MDC4MA) ,
& IMATRA(MDC,MSC) ,
& IMATDA(MDC,MSC) , 40.23
& PLNL4S(MDC,MSC,NPTST) , 40.00
& WWAWG(*)
REAL :: REDC0 (MDC,MSC,MREDS) 40.85
REAL :: REDC1 (MDC,MSC,MREDS) 40.85
!
INTEGER WWINT(*) ,
& IDCMIN(MSC) ,
& IDCMAX(MSC)
!
LOGICAL PERCIR
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'SWSNL2')
!
IDP = WWINT(1)
IDP1 = WWINT(2)
IDM = WWINT(3)
IDM1 = WWINT(4)
ISP = WWINT(5)
ISP1 = WWINT(6)
ISM = WWINT(7)
ISM1 = WWINT(8)
ISLOW = WWINT(9)
ISHGH = WWINT(10)
ISCLW = WWINT(11)
ISCHG = WWINT(12)
IDLOW = WWINT(13)
IDHGH = WWINT(14)
!
AWG1 = WWAWG(1)
AWG2 = WWAWG(2)
AWG3 = WWAWG(3)
AWG4 = WWAWG(4)
AWG5 = WWAWG(5)
AWG6 = WWAWG(6)
AWG7 = WWAWG(7)
AWG8 = WWAWG(8)
!
! *** Initialize auxiliary arrays per gridpoint ***
!
DO ID = MDC4MI, MDC4MA
DO IS = MSC4MI, MSC4MA
UE(IS,ID) = 0.
SA1(IS,ID) = 0.
SA2(IS,ID) = 0.
SFNL(IS,ID) = 0.
ENDDO
ENDDO
!
! *** Calculate prop. constant. ***
! *** Calculate factor R(X) to calculate the NL wave-wave ***
! *** interaction for shallow water ***
! *** SNLC1 = 1/GRAV**4 *** 40.17
!
SNLCS1 = PQUAD(3) 34.00
SNLCS2 = PQUAD(4) 34.00
SNLCS3 = PQUAD(5) 34.00
X = MAX ( 0.75 * DEP2(KCGRD(1)) * KMESPC , 0.5 )
X2 = MAX ( -1.E15, SNLCS3*X)
CONS = SNLC1 * ( 1. + SNLCS1/X * (1.-SNLCS2*X) * EXP(X2))
JACOBI = 2. * PI
!
! *** check whether the spectral domain is periodic in ***
! *** direction space and if so modify boundaries ***
!
PERCIR = .FALSE.
IF ( IDDLOW .EQ. 1 .AND. IDDTOP .EQ. MDC ) THEN
! *** periodic in theta -> spectrum can be folded ***
! *** (can only occur in presence of a current) ***
IDCLOW = 1
IDCHGH = MDC
IIID = 0
PERCIR = .TRUE.
ELSE
! *** different sectors per sweep -> extend range with IIID ***
IIID = MAX ( IDM1 , IDP1 )
IDCLOW = IDLOW
IDCHGH = IDHGH
ENDIF
!
! *** Prepare auxiliary spectrum ***
! *** set action original spectrum in array UE ***
!
DO IDDUM = IDLOW - IIID , IDHGH + IIID
ID = MOD ( IDDUM - 1 + MDC , MDC ) + 1
DO IS = 1, MSC
UE(IS,IDDUM) = AC2(ID,IS,KCGRD(1)) * SPCSIG(IS) * JACOBI 30.72
ENDDO
ENDDO
!
! *** set values in the areas 2 for IS > MSC+1 ***
!
DO IS = MSC+1, ISHGH
DO ID = IDLOW - IIID , IDHGH + IIID
UE (IS,ID) = UE(IS-1,ID) * FACHFR
ENDDO
ENDDO
!
! *** Calculate interactions ***
! *** Energy at interacting bins ***
!
DO IS = ISCLW, ISCHG
DO ID = IDCLOW , IDCHGH
E00 = UE(IS ,ID )
EP1 = AWG1 * UE(IS+ISP1,ID+IDP1) +
& AWG2 * UE(IS+ISP1,ID+IDP ) +
& AWG3 * UE(IS+ISP ,ID+IDP1) +
& AWG4 * UE(IS+ISP ,ID+IDP )
EM1 = AWG5 * UE(IS+ISM1,ID-IDM1) +
& AWG6 * UE(IS+ISM1,ID-IDM ) +
& AWG7 * UE(IS+ISM ,ID-IDM1) +
& AWG8 * UE(IS+ISM ,ID-IDM )
!
EP2 = AWG1 * UE(IS+ISP1,ID-IDP1) +
& AWG2 * UE(IS+ISP1,ID-IDP ) +
& AWG3 * UE(IS+ISP ,ID-IDP1) +
& AWG4 * UE(IS+ISP ,ID-IDP )
EM2 = AWG5 * UE(IS+ISM1,ID+IDM1) +
& AWG6 * UE(IS+ISM1,ID+IDM ) +
& AWG7 * UE(IS+ISM ,ID+IDM1) +
& AWG8 * UE(IS+ISM ,ID+IDM )
!
! *** Contribution to interactions ***
! *** CONS is the shallow water factor for the NL interact. ***
!
FACTOR = CONS * AF11(IS) * E00
!
SA1A = E00 * ( EP1*DAL1 + EM1*DAL2 ) * PQUAD(2) 40.17
SA1B = SA1A - EP1*EM1*DAL3 * PQUAD(2) 40.17
SA2A = E00 * ( EP2*DAL1 + EM2*DAL2 ) * PQUAD(2) 40.17
SA2B = SA2A - EP2*EM2*DAL3 * PQUAD(2) 40.17
!
SA1 (IS,ID) = FACTOR * SA1B
SA2 (IS,ID) = FACTOR * SA2B
!
IF(ITEST.GE.100 .AND. TESTFL) THEN
WRITE(PRINTF,9002) E00,EP1,EM1,EP2,EM2
9002 FORMAT (' E00 EP1 EM1 EP2 EM2 :',5E11.4)
WRITE(PRINTF,9003) SA1A,SA1B,SA2A,SA2B
9003 FORMAT (' SA1A SA1B SA2A SA2B :',4E11.4)
WRITE(PRINTF,9004) IS,ID,SA1(IS,ID),SA2(IS,ID)
9004 FORMAT (' IS ID SA1() SA2() :',2I4,2E12.4)
WRITE(PRINTF,9005) FACTOR ,ISLOW
9005 FORMAT (' FACTOR ISLOW : ',E12.4,I4)
END IF
!
ENDDO
ENDDO
!
! *** Fold interactions to side angles if spectral domain ***
! *** is periodic in directional space ***
!
IF ( PERCIR ) THEN
DO ID = 1, IDHGH - MDC
ID0 = 1 - ID
DO IS = ISCLW, ISCHG
SA1 (IS,MDC+ID) = SA1 (IS , ID )
SA2 (IS,MDC+ID) = SA2 (IS , ID )
SA1 (IS, ID0 ) = SA1 (IS , MDC+ID0)
SA2 (IS, ID0 ) = SA2 (IS , MDC+ID0)
ENDDO
ENDDO
ENDIF
!
! *** Put source term together (To save space I=IS and J=ID ***
! *** is used) ***
!
DO I = 1, ISSTOP
SIGPI = SPCSIG(I) * JACOBI 30.72
DO J = IDCMIN(I), IDCMAX(I)
ID = MOD ( J - 1 + MDC , MDC ) + 1
SFNL(I,ID) = - 2. * ( SA1(I,J) + SA2(I,J) )
& + AWG1 * ( SA1(I-ISP1,J-IDP1) + SA2(I-ISP1,J+IDP1) )
& + AWG2 * ( SA1(I-ISP1,J-IDP ) + SA2(I-ISP1,J+IDP ) )
& + AWG3 * ( SA1(I-ISP ,J-IDP1) + SA2(I-ISP ,J+IDP1) )
& + AWG4 * ( SA1(I-ISP ,J-IDP ) + SA2(I-ISP ,J+IDP ) )
& + AWG5 * ( SA1(I-ISM1,J+IDM1) + SA2(I-ISM1,J-IDM1) )
& + AWG6 * ( SA1(I-ISM1,J+IDM ) + SA2(I-ISM1,J-IDM ) )
& + AWG7 * ( SA1(I-ISM ,J+IDM1) + SA2(I-ISM ,J-IDM1) )
& + AWG8 * ( SA1(I-ISM ,J+IDM ) + SA2(I-ISM ,J-IDM ) )
!
! *** store results in rhv ***
! *** store results in rhs and main diagonal according *** 40.23
! *** to Patankar-rules *** 40.23
!
IF(TESTFL) PLNL4S(ID,I,IPTST) = SFNL(I,ID) / SIGPI 40.00
IF (SFNL(I,ID).GT.0.) THEN 40.23
IMATRA(ID,I) = IMATRA(ID,I) + SFNL(I,ID) / SIGPI
REDC0(ID,I,1)= REDC0(ID,I,1)+ SFNL(I,ID) / SIGPI 40.85
ELSE
IMATDA(ID,I) = IMATDA(ID,I) - SFNL(I,ID) /
& MAX(1.E-18,AC2(ID,I,KCGRD(1))*SIGPI)
REDC1(ID,I,1)= REDC1(ID,I,1)+ SFNL(I,ID) / 40.85
& MAX(1.E-18,AC2(ID,I,KCGRD(1))*SIGPI) 40.85
END IF
!
ENDDO
ENDDO
!
! *** test output ***
!
IF (ITEST .GE. 40 .AND. TESTFL) THEN
WRITE(PRINTF,*) ' SWSNL2 subroutine '
WRITE(PRINTF,9011) IDP, IDP1, IDM, IDM1
9011 FORMAT (' IDP IDP1 IDM IDM1 :',4I5)
WRITE (PRINTF,9013) ISP, ISP1, ISM, ISM1
9013 FORMAT (' ISP ISP1 ISM ISM1 :',4I5)
WRITE (PRINTF,9015) ISHGH, IDDLOW, IDDTOP
9015 FORMAT (' ISHG IDDLOW IDDTOP :',3I5)
WRITE(PRINTF,9016) ISCLW, ISCHG, IDLOW, IDHGH
9016 FORMAT (' ICLW ICHG IDLOW IDHGH :',4I5)
WRITE (PRINTF,9017) AWG1, AWG2, AWG3, AWG4
9017 FORMAT (' AWG1 AWG2 AWG3 AWG4 :',4E12.4)
WRITE (PRINTF,9018) AWG5, AWG6, AWG7, AWG8
9018 FORMAT (' AWG5 AWG6 AWG7 AWG8 :',4E12.4)
WRITE (PRINTF,9019) MSC4MI, MSC4MA, MDC4MI, MDC4MA
9019 FORMAT (' S4MI S4MA D4MI D4MA :',4I6)
WRITE(PRINTF,9020) SNLC1,X,X2,CONS
9020 FORMAT (' SNLC1 X X2 CONS :',4E12.4)
WRITE(PRINTF,9021) DEP2(KCGRD(1)),KMESPC, FACHFR,PI
9021 FORMAT (' DEPTH KMESPC FACHFR PI:',4E12.4)
WRITE(PRINTF,9023) JACOBI,ISLOW
9023 FORMAT (' JACOBI ISLOW :',E12.4,I4)
WRITE(PRINTF,*)
END IF
!
RETURN
! End of SWSNL2
END
!
!************************************************************
!
SUBROUTINE SWSNL3 ( WWINT ,WWAWG , 40.17
& UE ,SA1 ,SA2 ,SPCSIG ,SNLC1 , 40.17
& DAL1 ,DAL2 ,DAL3 ,SFNL ,DEP2 , 40.17
& AC2 ,KMESPC ,MEMNL4 ,FACHFR ) 40.17
!
!*******************************************************************
!
USE SWCOMM3 40.41
USE SWCOMM4 40.41
USE OCPCOMM4 40.41
USE M_SNL4 40.17
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: H.L. Tolman, R.C. Ris |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 30.72: IJsbrand Haagsma
! 40.17: IJsbrand Haagsma
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN
! 40.17, Dec. 01: Implemented Multiple DIA
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose
!
! Calculate non-linear interaction using the discrete interaction
! approximation (Hasselmann and Hasselmann 1985; WAMDI group 1988)
! for the full circle (option if a current is present). Note: using
! this subroutine requires an additional array with size
! (MXC*MYC*MDC*MSC). This requires more internal memory but can
! speed up the computations sigificantly if a current is present.
!
! 3. Method
!
! Discrete interaction approximation. To make interpolation simple,
! the interactions are calculated in a "folded" space.
!
! Frequencies -->
! +---+---------------------+---------+- IDHGH
! d | 3 : 2 : 2 |
! i + - + - - - - - - - - - - + - - - - +- MDC
! r | : : |
! e | 3 : original spectrum : 1 |
! c | : : |
! t. + - + - - - - - - - - - - + - - - - +- 1
! | 3 : 2 : 2 |
! +---+---------------------+---------+- IDLOW
! | | | ^ |
! ISLOW 1 MSC | ISHGH
! | |
! ISCLW ISCHG
! lowest discrete highest discrete
! central bin central bin
!
! 1 : Extra tail added beyond MSC
! 2 : Spectrum copied outside ID range
! 3 : Empty bins at low frequencies
!
! ISLOW = 1 + ISM1
! ISHGH = MSC + ISP1 - ISM1
! ISCLW = 1
! ISCHG = MSC - ISM1
! IDLOW = 1 - MAX(IDM1,IDP1)
! IDHGH = MDC + MAX(IDM1,IDP1)
!
! Relative offsets of interpolation points around central bin
! "#" and corresponding numbers of AWGn :
!
! ISM1 ISM
! 5 7 T |
! IDM1 +------+ H +
! | | E | ISP ISP1
! | \ | T | 3 1
! IDM +------+ A + +---------+ IDP1
! 6 \8 | | |
! | | / |
! \ + +---------+ IDP
! | /4 2
! \ | /
! -+-----+------+-------#--------+---------+----------+
! | FREQ.
!
!
! 4. Argument variables
!
! MCGRD : number of wet grid points of the computational grid
! MDC : grid points in theta-direction of computational grid
! MDC4MA: highest array counter in directional space (Snl4)
! MDC4MI: lowest array counter in directional space (Snl4)
! MSC : grid points in sigma-direction of computational grid
! MSC4MA: highest array counter in frequency space (Snl4)
! MSC4MI: lowest array counter in frequency space (Snl4)
! WWINT : counters for quadruplet interactions
!
INTEGER WWINT(*)
!
! AC2 : action density
! AF11 : scaling frequency
! DAL1 : coefficient for the quadruplet interactions
! DAL2 : coefficient for the quadruplet interactions
! DAL3 : coefficient for the quadruplet interactions
! DEP2 : depth
! FACHFR
! KMESPC: mean average wavenumber over full spectrum
! MEMNL4
! PI : circular constant
! SA1 : interaction contribution of first quadruplet (unfolded space)
! SA2 : interaction contribution of second quadruplet (unfolded space)
! SFNL
! SNLC1
! SPCSIG: relative frequencies in computational domain in sigma-space
! UE : "unfolded" spectrum
! WWAWG : weight coefficients for the quadruplet interactions
!
REAL DAL1, DAL2, DAL3, FACHFR, KMESPC, SNLC1 40.17
REAL AC2(MDC,MSC,MCGRD)
REAL DEP2(MCGRD)
REAL MEMNL4(MDC,MSC,MCGRD)
REAL SA1(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
REAL SA2(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
REAL SFNL(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
REAL SPCSIG(MSC) 30.72
REAL UE(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
REAL WWAWG(*)
!
! 7. Common blocks used
!
!
! 8. Subroutines used
!
! ---
!
! 9. Subroutines calling
!
! SOURCE (in SWANCOM1)
!
! 12. Structure
!
! -------------------------------------------
! Initialisations.
! Calculate proportionality constant.
! Prepare auxiliary spectrum.
! Calculate (unfolded) interactions :
! -----------------------------------------
! Energy at interacting bins
! Contribution to interactions
! Fold interactions to side angles
! -----------------------------------------
! Put source term together
! -------------------------------------------
!
! 13. Source text
!
!*******************************************************************
!
INTEGER IS ,ID ,ID0 ,I ,J ,
& ISHGH ,IDLOW ,IDHGH ,ISP ,ISP1 ,
& IDP ,IDP1 ,ISM ,ISM1 ,IDM ,IDM1 ,
& ISCLW ,ISCHG
!
REAL X ,X2 ,CONS ,FACTOR ,SNLCS2 ,
& SNLCS3 ,E00 ,EP1 ,EM1 ,EP2 ,EM2 ,
& SA1A ,SA1B ,SA2A ,SA2B ,
& AWG1 ,AWG2 ,AWG3 ,AWG4 ,AWG5 ,AWG6 ,
& AWG7 ,AWG8 ,
& JACOBI ,SIGPI
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'SWSNL3')
!
IDP = WWINT(1)
IDP1 = WWINT(2)
IDM = WWINT(3)
IDM1 = WWINT(4)
ISP = WWINT(5)
ISP1 = WWINT(6)
ISM = WWINT(7)
ISM1 = WWINT(8)
ISLOW = WWINT(9)
ISHGH = WWINT(10)
ISCLW = WWINT(11)
ISCHG = WWINT(12)
IDLOW = WWINT(13)
IDHGH = WWINT(14)
!
AWG1 = WWAWG(1)
AWG2 = WWAWG(2)
AWG3 = WWAWG(3)
AWG4 = WWAWG(4)
AWG5 = WWAWG(5)
AWG6 = WWAWG(6)
AWG7 = WWAWG(7)
AWG8 = WWAWG(8)
!
! *** Initialize auxiliary arrays per gridpoint ***
!
DO ID = MDC4MI, MDC4MA
DO IS = MSC4MI, MSC4MA
UE(IS,ID) = 0.
SA1(IS,ID) = 0.
SA2(IS,ID) = 0.
SFNL(IS,ID) = 0.
ENDDO
ENDDO
!
! *** Calculate prop. constant. ***
! *** Calculate factor R(X) to calculate the NL wave-wave ***
! *** interaction for shallow water ***
! *** SNLC1 = 1/GRAV**4 *** 40.17
!
SNLCS1 = PQUAD(3) 34.00
SNLCS2 = PQUAD(4) 34.00
SNLCS3 = PQUAD(5) 34.00
X = MAX ( 0.75 * DEP2(KCGRD(1)) * KMESPC , 0.5 ) 30.21
X2 = MAX ( -1.E15, SNLCS3*X)
CONS = SNLC1 * ( 1. + SNLCS1/X * (1.-SNLCS2*X) * EXP(X2))
JACOBI = 2. * PI
!
! *** extend the area with action density at periodic boundaries ***
!
DO IDDUM = IDLOW, IDHGH
ID = MOD ( IDDUM - 1 + MDC , MDC ) + 1
DO IS=1, MSC
UE (IS,IDDUM) = AC2(ID,IS,KCGRD(1)) * SPCSIG(IS) * JACOBI 30.72
ENDDO
ENDDO
!
DO IS = MSC+1, ISHGH
DO ID = IDLOW, IDHGH
UE(IS,ID) = UE(IS-1,ID) * FACHFR
ENDDO
ENDDO
!
! *** Calculate (unfolded) interactions ***
! *** Energy at interacting bins ***
!
DO IS = ISCLW, ISCHG
DO ID = 1, MDC
E00 = UE(IS ,ID )
EP1 = AWG1 * UE(IS+ISP1,ID+IDP1) +
& AWG2 * UE(IS+ISP1,ID+IDP ) +
& AWG3 * UE(IS+ISP ,ID+IDP1) +
& AWG4 * UE(IS+ISP ,ID+IDP )
EM1 = AWG5 * UE(IS+ISM1,ID-IDM1) +
& AWG6 * UE(IS+ISM1,ID-IDM ) +
& AWG7 * UE(IS+ISM ,ID-IDM1) +
& AWG8 * UE(IS+ISM ,ID-IDM )
EP2 = AWG1 * UE(IS+ISP1,ID-IDP1) +
& AWG2 * UE(IS+ISP1,ID-IDP ) +
& AWG3 * UE(IS+ISP ,ID-IDP1) +
& AWG4 * UE(IS+ISP ,ID-IDP )
EM2 = AWG5 * UE(IS+ISM1,ID+IDM1) +
& AWG6 * UE(IS+ISM1,ID+IDM ) +
& AWG7 * UE(IS+ISM ,ID+IDM1) +
& AWG8 * UE(IS+ISM ,ID+IDM )
!
! Contribution to interactions
!
FACTOR = CONS * AF11(IS) * E00
!
SA1A = E00 * ( EP1*DAL1 + EM1*DAL2 ) * PQUAD(2) 40.17
SA1B = SA1A - EP1*EM1*DAL3 * PQUAD(2) 40.17
SA2A = E00 * ( EP2*DAL1 + EM2*DAL2 ) * PQUAD(2) 40.17
SA2B = SA2A - EP2*EM2*DAL3 * PQUAD(2) 40.17
!
SA1 (IS,ID) = FACTOR * SA1B
SA2 (IS,ID) = FACTOR * SA2B
!
IF(ITEST.GE.100 .AND. TESTFL) THEN
WRITE(PRINTF,9002) E00,EP1,EM1,EP2,EM2
9002 FORMAT (' E00 EP1 EM1 EP2 EM2 :',5E11.4)
WRITE(PRINTF,9003) SA1A,SA1B,SA2A,SA2B
9003 FORMAT (' SA1A SA1B SA2A SA2B :',4E11.4)
WRITE(PRINTF,9004) IS,ID,SA1(IS,ID),SA2(IS,ID)
9004 FORMAT (' IS ID SA1() SA2() :',2I4,2E12.4)
WRITE(PRINTF,9005) FACTOR,JACOBI
9005 FORMAT (' FACTOR JACOBI : ',2E12.4)
END IF
!
ENDDO
ENDDO
!
! *** Fold interactions to side angles -> domain in theta is ***
! *** periodic ***
!
DO ID = 1, IDHGH - MDC
ID0 = 1 - ID
DO IS = ISCLW, ISCHG
SA1 (IS,MDC+ID) = SA1 (IS, ID )
SA2 (IS,MDC+ID) = SA2 (IS, ID )
SA1 (IS, ID0 ) = SA1 (IS,MDC+ID0)
SA2 (IS, ID0 ) = SA2 (IS,MDC+ID0)
ENDDO
ENDDO
!
! *** Put source term together (To save space I=IS and ***
! *** J=MDC is used) ---- ***
!
DO I = 1, MSC
SIGPI = SPCSIG(I) * JACOBI 30.72
DO J = 1, MDC
SFNL(I,J) = - 2. * ( SA1(I,J) + SA2(I,J) )
& + AWG1 * ( SA1(I-ISP1,J-IDP1) + SA2(I-ISP1,J+IDP1) )
& + AWG2 * ( SA1(I-ISP1,J-IDP ) + SA2(I-ISP1,J+IDP ) )
& + AWG3 * ( SA1(I-ISP ,J-IDP1) + SA2(I-ISP ,J+IDP1) )
& + AWG4 * ( SA1(I-ISP ,J-IDP ) + SA2(I-ISP ,J+IDP ) )
& + AWG5 * ( SA1(I-ISM1,J+IDM1) + SA2(I-ISM1,J-IDM1) )
& + AWG6 * ( SA1(I-ISM1,J+IDM ) + SA2(I-ISM1,J-IDM ) )
& + AWG7 * ( SA1(I-ISM ,J+IDM1) + SA2(I-ISM ,J-IDM1) )
& + AWG8 * ( SA1(I-ISM ,J+IDM ) + SA2(I-ISM ,J-IDM ) )
!
! *** store value in auxiliary array and use values in ***
! *** next four sweeps (see subroutine FILNL3) ***
!
MEMNL4(J,I,KCGRD(1)) = SFNL(I,J) / SIGPI 30.21
ENDDO
ENDDO
!
! *** test output ***
!
IF (ITEST .GE. 50 .AND. TESTFL) THEN
WRITE(PRINTF,*)
WRITE(PRINTF,*) ' SWSNL3 subroutine '
WRITE(PRINTF,9011) IDP, IDP1, IDM, IDM1
9011 FORMAT (' IDP IDP1 IDM IDM1 :',4I5)
WRITE (PRINTF,9013) ISP, ISP1, ISM, ISM1
9013 FORMAT (' ISP ISP1 ISM ISM1 :',4I5)
WRITE (PRINTF,9015) ISLOW, ISHGH, IDLOW, IDHGH
9015 FORMAT (' ISLOW ISHG IDLOW IDHG :',4I5)
WRITE(PRINTF,9016) ISCLW, ISCHG, JACOBI
9016 FORMAT (' ICLW ICHG JACOBI :',2I5,E12.4)
WRITE (PRINTF,9017) AWG1, AWG2, AWG3, AWG4
9017 FORMAT (' AWG1 AWG2 AWG3 AWG4 :',4E12.4)
WRITE (PRINTF,9018) AWG5, AWG6, AWG7, AWG8
9018 FORMAT (' AWG5 AWG6 AWG7 AWG8 :',4E12.4)
WRITE (PRINTF,9019) MSC4MI, MSC4MA, MDC4MI, MDC4MA
9019 FORMAT (' S4MI S4MA D4MI D4MA :',4I6)
WRITE(PRINTF,9020) SNLC1,X,X2,CONS
9020 FORMAT (' SNLC1 X X2 CONS :',4E12.4)
WRITE(PRINTF,9021) DEP2(KCGRD(1)),KMESPC,FACHFR,PI
9021 FORMAT (' DEPTH KMESPC FACHFR PI:',4E12.4)
WRITE(PRINTF,*)
!
! *** value source term in every bin ***
!
IF(ITEST.GE. 150 ) THEN
DO I=1, MSC
DO J=1, MDC
WRITE(PRINTF,2006) I,J,MEMNL4(J,I,KCGRD(1)),SFNL(I,J), 30.21
& SPCSIG(I) 30.72
2006 FORMAT (' I J MEMNL() SFNL() SPCSIG:',2I4,3E12.4) 30.72
ENDDO
ENDDO
END IF
END IF
!
RETURN
!
END SUBROUTINE SWSNL3
!
!*******************************************************************
!
SUBROUTINE SWSNL4 (WWINT ,WWAWG ,
& SPCSIG ,SNLC1 ,
& DAL1 ,DAL2 ,DAL3 ,DEP2 ,
& AC2 ,KMESPC ,MEMNL4 ,FACHFR ,
& IDIA ,ITER )
!
!*******************************************************************
!
USE SWCOMM3 40.41
USE SWCOMM4 40.41
USE OCPCOMM4 40.41
USE M_SNL4
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: H.L. Tolman, R.C. Ris |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 40.17: IJsbrand Haagsma
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 40.17, Dec. 01: New Subroutine based on SWSNL3
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose
!
! Calculate non-linear interaction using the discrete interaction
! approximation (Hasselmann and Hasselmann 1985; WAMDI group 1988)
! for the full circle (option if a current is present). Note: using
! this subroutine requires an additional array with size
! (MXC*MYC*MDC*MSC). This requires more internal memory but can
! speed up the computations sigificantly if a current is present.
!
! 3. Method
!
! Discrete interaction approximation. To make interpolation simple,
! the interactions are calculated in a "folded" space.
!
! Frequencies -->
! +---+---------------------+---------+- IDHGH
! d | 3 : 2 : 2 |
! i + - + - - - - - - - - - - + - - - - +- MDC
! r | : : |
! e | 3 : original spectrum : 1 |
! c | : : |
! t. + - + - - - - - - - - - - + - - - - +- 1
! | 3 : 2 : 2 |
! +---+---------------------+---------+- IDLOW
! | | | ^ |
! ISLOW 1 MSC | ISHGH
! | |
! ISCLW ISCHG
! lowest discrete highest discrete
! central bin central bin
!
! 1 : Extra tail added beyond MSC
! 2 : Spectrum copied outside ID range
! 3 : Empty bins at low frequencies
!
! ISLOW = 1 + ISM1
! ISHGH = MSC + ISP1 - ISM1
! ISCLW = 1
! ISCHG = MSC - ISM1
! IDLOW = 1 - MAX(IDM1,IDP1)
! IDHGH = MDC + MAX(IDM1,IDP1)
!
! Relative offsets of interpolation points around central bin
! "#" and corresponding numbers of AWGn :
!
! ISM1 ISM
! 5 7 T |
! IDM1 +------+ H +
! | | E | ISP ISP1
! | \ | T | 3 1
! IDM +------+ A + +---------+ IDP1
! 6 \8 | | |
! | | / |
! \ + +---------+ IDP
! | /4 2
! \ | /
! -+-----+------+-------#--------+---------+----------+
! | FREQ.
!
!
! 4. Argument variables
!
! MCGRD : number of wet grid points of the computational grid
! MDC : grid points in theta-direction of computational grid
! MDC4MA: highest array counter in directional space (Snl4)
! MDC4MI: lowest array counter in directional space (Snl4)
! MSC : grid points in sigma-direction of computational grid
! MSC4MA: highest array counter in frequency space (Snl4)
! MSC4MI: lowest array counter in frequency space (Snl4)
! WWINT : counters for quadruplet interactions
!
INTEGER WWINT(*)
INTEGER IDIA
!
! AC2 : action density
! AF11 : scaling frequency
! DAL1 : coefficient for the quadruplet interactions
! DAL2 : coefficient for the quadruplet interactions
! DAL3 : coefficient for the quadruplet interactions
! DEP2 : depth
! FACHFR
! KMESPC: mean average wavenumber over full spectrum
! MEMNL4
! PI : circular constant
! SA1 : interaction contribution of first quadruplet (unfolded space)
! SA2 : interaction contribution of second quadruplet (unfolded space)
! SFNL
! SNLC1
! SPCSIG: relative frequencies in computational domain in sigma-space
! UE : "unfolded" spectrum
! WWAWG : weight coefficients for the quadruplet interactions
!
REAL DAL1, DAL2, DAL3, FACHFR, KMESPC, SNLC1
REAL AC2(MDC,MSC,MCGRD)
REAL DEP2(MCGRD)
REAL MEMNL4(MDC,MSC,MCGRD)
REAL SPCSIG(MSC)
REAL WWAWG(*)
!
! 6. Local variables
!
REAL, ALLOCATABLE :: SA1(:,:), SA2(:,:), SFNL(:,:), UE(:,:)
!
! 7. Common blocks used
!
!
! 8. Subroutines used
!
! ---
!
! 9. Subroutines calling
!
! SOURCE (in SWANCOM1)
!
! 12. Structure
!
! -------------------------------------------
! Initialisations.
! Calculate proportionality constant.
! Prepare auxiliary spectrum.
! Calculate (unfolded) interactions :
! -----------------------------------------
! Energy at interacting bins
! Contribution to interactions
! Fold interactions to side angles
! -----------------------------------------
! Put source term together
! -------------------------------------------
!
! 13. Source text
!
!*******************************************************************
!
INTEGER IS ,ID ,ID0 ,I ,J ,
& ISHGH ,IDLOW ,IDHGH ,ISP ,ISP1 ,
& IDP ,IDP1 ,ISM ,ISM1 ,IDM ,IDM1 ,
& ISCLW ,ISCHG
!
REAL X ,X2 ,CONS ,FACTOR ,SNLCS2 ,
& SNLCS3 ,E00 ,EP1 ,EM1 ,EP2 ,EM2 ,
& SA1A ,SA1B ,SA2A ,SA2B ,
& AWG1 ,AWG2 ,AWG3 ,AWG4 ,AWG5 ,AWG6 ,
& AWG7 ,AWG8 ,
& JACOBI ,SIGPI
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'SWSNL4')
ALLOCATE(SA1(MSC4MI:MSC4MA,MDC4MI:MDC4MA))
ALLOCATE(SA2(MSC4MI:MSC4MA,MDC4MI:MDC4MA))
ALLOCATE(SFNL(MSC4MI:MSC4MA,MDC4MI:MDC4MA))
ALLOCATE(UE(MSC4MI:MSC4MA,MDC4MI:MDC4MA))
!
IDP = WWINT(1)
IDP1 = WWINT(2)
IDM = WWINT(3)
IDM1 = WWINT(4)
ISP = WWINT(5)
ISP1 = WWINT(6)
ISM = WWINT(7)
ISM1 = WWINT(8)
ISLOW = WWINT(9)
ISHGH = WWINT(10)
ISCLW = WWINT(11)
ISCHG = WWINT(12)
IDLOW = WWINT(13)
IDHGH = WWINT(14)
!
AWG1 = WWAWG(1)
AWG2 = WWAWG(2)
AWG3 = WWAWG(3)
AWG4 = WWAWG(4)
AWG5 = WWAWG(5)
AWG6 = WWAWG(6)
AWG7 = WWAWG(7)
AWG8 = WWAWG(8)
!
! *** Initialize auxiliary arrays per gridpoint ***
!
DO ID = MDC4MI, MDC4MA
DO IS = MSC4MI, MSC4MA
UE(IS,ID) = 0.
SA1(IS,ID) = 0.
SA2(IS,ID) = 0.
SFNL(IS,ID) = 0.
ENDDO
ENDDO
!
! *** Calculate prop. constant. ***
! *** Calculate factor R(X) to calculate the NL wave-wave ***
! *** interaction for shallow water ***
! *** SNLC1 = 1/GRAV**4 ***
!
SNLCS1 = PQUAD(3) 34.00
SNLCS2 = PQUAD(4) 34.00
SNLCS3 = PQUAD(5) 34.00
X = MAX ( 0.75 * DEP2(KCGRD(1)) * KMESPC , 0.5 )
X2 = MAX ( -1.E15, SNLCS3*X)
CONS = SNLC1 * ( 1. + SNLCS1/X * (1.-SNLCS2*X) * EXP(X2))
JACOBI = 2. * PI
!
! *** extend the area with action density at periodic boundaries ***
!
DO IDDUM = IDLOW, IDHGH
ID = MOD ( IDDUM - 1 + MDC , MDC ) + 1
DO IS=1, MSC
UE (IS,IDDUM) = AC2(ID,IS,KCGRD(1)) * SPCSIG(IS) * JACOBI
ENDDO
ENDDO
!
DO IS = MSC+1, ISHGH
DO ID = IDLOW, IDHGH
UE(IS,ID) = UE(IS-1,ID) * FACHFR
ENDDO
ENDDO
!
! *** Calculate (unfolded) interactions ***
! *** Energy at interacting bins ***
!
DO IS = ISCLW, ISCHG
DO ID = 1, MDC
E00 = UE(IS ,ID )
EP1 = AWG1 * UE(IS+ISP1,ID+IDP1) +
& AWG2 * UE(IS+ISP1,ID+IDP ) +
& AWG3 * UE(IS+ISP ,ID+IDP1) +
& AWG4 * UE(IS+ISP ,ID+IDP )
EM1 = AWG5 * UE(IS+ISM1,ID-IDM1) +
& AWG6 * UE(IS+ISM1,ID-IDM ) +
& AWG7 * UE(IS+ISM ,ID-IDM1) +
& AWG8 * UE(IS+ISM ,ID-IDM )
EP2 = AWG1 * UE(IS+ISP1,ID-IDP1) +
& AWG2 * UE(IS+ISP1,ID-IDP ) +
& AWG3 * UE(IS+ISP ,ID-IDP1) +
& AWG4 * UE(IS+ISP ,ID-IDP )
EM2 = AWG5 * UE(IS+ISM1,ID+IDM1) +
& AWG6 * UE(IS+ISM1,ID+IDM ) +
& AWG7 * UE(IS+ISM ,ID+IDM1) +
& AWG8 * UE(IS+ISM ,ID+IDM )
!
! Contribution to interactions
!
FACTOR = CONS * AF11(IS) * E00
!
SA1A = E00 * ( EP1*DAL1 + EM1*DAL2 ) * CNL4_1(IDIA)
SA1B = SA1A - EP1*EM1*DAL3 * CNL4_2(IDIA)
SA2A = E00 * ( EP2*DAL1 + EM2*DAL2 ) * CNL4_1(IDIA)
SA2B = SA2A - EP2*EM2*DAL3 * CNL4_2(IDIA)
!
SA1 (IS,ID) = FACTOR * SA1B
SA2 (IS,ID) = FACTOR * SA2B
!
IF(ITEST.GE.100 .AND. TESTFL) THEN
WRITE(PRINTF,9002) E00,EP1,EM1,EP2,EM2
9002 FORMAT (' E00 EP1 EM1 EP2 EM2 :',5E11.4)
WRITE(PRINTF,9003) SA1A,SA1B,SA2A,SA2B
9003 FORMAT (' SA1A SA1B SA2A SA2B :',4E11.4)
WRITE(PRINTF,9004) IS,ID,SA1(IS,ID),SA2(IS,ID)
9004 FORMAT (' IS ID SA1() SA2() :',2I4,2E12.4)
WRITE(PRINTF,9005) FACTOR,JACOBI
9005 FORMAT (' FACTOR JACOBI : ',2E12.4)
END IF
!
ENDDO
ENDDO
!
! *** Fold interactions to side angles -> domain in theta is ***
! *** periodic ***
!
DO ID = 1, IDHGH - MDC
ID0 = 1 - ID
DO IS = ISCLW, ISCHG
SA1 (IS,MDC+ID) = SA1 (IS, ID )
SA2 (IS,MDC+ID) = SA2 (IS, ID )
SA1 (IS, ID0 ) = SA1 (IS,MDC+ID0)
SA2 (IS, ID0 ) = SA2 (IS,MDC+ID0)
ENDDO
ENDDO
!
! *** Put source term together (To save space I=IS and ***
! *** J=MDC is used) ***
!
FAC = 1. 40.17
DO I = 1, MSC
SIGPI = SPCSIG(I) * JACOBI 30.72
DO J = 1, MDC
SFNL(I,J) = - 2. * ( SA1(I,J) + SA2(I,J) )
& + AWG1 * ( SA1(I-ISP1,J-IDP1) + SA2(I-ISP1,J+IDP1) )
& + AWG2 * ( SA1(I-ISP1,J-IDP ) + SA2(I-ISP1,J+IDP ) )
& + AWG3 * ( SA1(I-ISP ,J-IDP1) + SA2(I-ISP ,J+IDP1) )
& + AWG4 * ( SA1(I-ISP ,J-IDP ) + SA2(I-ISP ,J+IDP ) )
& + AWG5 * ( SA1(I-ISM1,J+IDM1) + SA2(I-ISM1,J-IDM1) )
& + AWG6 * ( SA1(I-ISM1,J+IDM ) + SA2(I-ISM1,J-IDM ) )
& + AWG7 * ( SA1(I-ISM ,J+IDM1) + SA2(I-ISM ,J-IDM1) )
& + AWG8 * ( SA1(I-ISM ,J+IDM ) + SA2(I-ISM ,J-IDM ) )
!
! *** store value in auxiliary array and use values in ***
! *** next four sweeps (see subroutine FILNL3) ***
!
IF (IDIA.EQ.1) THEN
MEMNL4(J,I,KCGRD(1)) = FAC * SFNL(I,J) / SIGPI
ELSE
MEMNL4(J,I,KCGRD(1)) = MEMNL4(J,I,KCGRD(1)) +
& FAC * SFNL(I,J) / SIGPI
END IF
ENDDO
ENDDO
!
! *** test output ***
!
IF (ITEST .GE. 50 .AND. TESTFL) THEN
WRITE(PRINTF,*)
WRITE(PRINTF,*) ' SWSNL4 subroutine '
WRITE(PRINTF,9011) IDP, IDP1, IDM, IDM1
9011 FORMAT (' IDP IDP1 IDM IDM1 :',4I5)
WRITE (PRINTF,9013) ISP, ISP1, ISM, ISM1
9013 FORMAT (' ISP ISP1 ISM ISM1 :',4I5)
WRITE (PRINTF,9015) ISLOW, ISHGH, IDLOW, IDHGH
9015 FORMAT (' ISLOW ISHG IDLOW IDHG :',4I5)
WRITE(PRINTF,9016) ISCLW, ISCHG, JACOBI
9016 FORMAT (' ICLW ICHG JACOBI :',2I5,E12.4)
WRITE (PRINTF,9017) AWG1, AWG2, AWG3, AWG4
9017 FORMAT (' AWG1 AWG2 AWG3 AWG4 :',4E12.4)
WRITE (PRINTF,9018) AWG5, AWG6, AWG7, AWG8
9018 FORMAT (' AWG5 AWG6 AWG7 AWG8 :',4E12.4)
WRITE (PRINTF,9019) MSC4MI, MSC4MA, MDC4MI, MDC4MA
9019 FORMAT (' S4MI S4MA D4MI D4MA :',4I6)
WRITE(PRINTF,9020) SNLC1,X,X2,CONS
9020 FORMAT (' SNLC1 X X2 CONS :',4E12.4)
WRITE(PRINTF,9021) DEP2(KCGRD(1)),KMESPC,FACHFR,PI
9021 FORMAT (' DEPTH KMESPC FACHFR PI:',4E12.4)
WRITE(PRINTF,*)
!
! *** value source term in every bin ***
!
IF(ITEST.GE. 150 ) THEN
DO I=1, MSC
DO J=1, MDC
WRITE(PRINTF,2006) I,J,MEMNL4(J,I,KCGRD(1)),SFNL(I,J),
& SPCSIG(I)
2006 FORMAT (' I J MEMNL() SFNL() SPCSIG:',2I4,3E12.4)
ENDDO
ENDDO
END IF
END IF
!
DEALLOCATE (SA1, SA2, SFNL, UE)
RETURN
!
END SUBROUTINE SWSNL4
!
!*********************************************************************
SUBROUTINE SWSNL8 (WWINT ,UE ,SA1 ,SA2 ,SPCSIG ,
& SNLC1 ,DAL1 ,DAL2 ,DAL3 ,SFNL ,
& DEP2 ,AC2 ,KMESPC ,MEMNL4 ,FACHFR )
!*********************************************************************
!
USE SWCOMM3 40.41
USE SWCOMM4 40.41
USE OCPCOMM4 40.41
USE M_SNL4
!
IMPLICIT NONE
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: H.L. Tolman, R.C. Ris |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 40.41: Marcel Zijlema
!
! 1. Updates
!
! 40.41, Sep. 04: piecewise constant interpolation instead
! of bi-linear one
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
!
! 2. Purpose
!
! Calculate non-linear interaction using the discrete interaction
! approximation (Hasselmann and Hasselmann 1985; WAMDI group 1988)
! for the full circle.
!
! 3. Method
!
! Discrete interaction approximation. To make interpolation simple,
! the interactions are calculated in a "folded" space.
!
! Frequencies -->
! +---+---------------------+---------+- IDHGH
! d | 3 : 2 : 2 |
! i + - + - - - - - - - - - - + - - - - +- MDC
! r | : : |
! e | 3 : original spectrum : 1 |
! c | : : |
! t. + - + - - - - - - - - - - + - - - - +- 1
! | 3 : 2 : 2 |
! +---+---------------------+---------+- IDLOW
! | | | ^ |
! ISLOW 1 MSC | ISHGH
! | |
! ISCLW ISCHG
! lowest discrete highest discrete
! central bin central bin
!
! 1 : Extra tail added beyond MSC
! 2 : Spectrum copied outside ID range
! 3 : Empty bins at low frequencies
!
! ISLOW = 1 + ISM1
! ISHGH = MSC + ISP1 - ISM1
! ISCLW = 1
! ISCHG = MSC - ISM1
! IDLOW = 1 - MAX(IDM1,IDP1)
! IDHGH = MDC + MAX(IDM1,IDP1)
!
! Note: using this subroutine requires an additional array
! with size MXC*MYC*MDC*MSC.
!
! 4. Argument variables
!
! MCGRD : number of wet grid points of the computational grid
! MDC : grid points in theta-direction of computational grid
! MDC4MA: highest array counter in directional space (Snl4)
! MDC4MI: lowest array counter in directional space (Snl4)
! MSC : grid points in sigma-direction of computational grid
! MSC4MA: highest array counter in frequency space (Snl4)
! MSC4MI: lowest array counter in frequency space (Snl4)
! WWINT : counters for quadruplet interactions
!
INTEGER WWINT(*)
!
! AC2 : action density
! AF11 : scaling frequency
! DAL1 : coefficient for the quadruplet interactions
! DAL2 : coefficient for the quadruplet interactions
! DAL3 : coefficient for the quadruplet interactions
! DEP2 : depth
! FACHFR
! KMESPC: mean average wavenumber over full spectrum
! MEMNL4
! PI : circular constant
! SA1 : interaction contribution of first quadruplet (unfolded space)
! SA2 : interaction contribution of second quadruplet (unfolded space)
! SFNL
! SNLC1
! SPCSIG: relative frequencies in computational domain in sigma-space
! UE : "unfolded" spectrum
!
REAL DAL1, DAL2, DAL3, FACHFR, KMESPC, SNLC1
REAL AC2(MDC,MSC,MCGRD)
REAL DEP2(MCGRD)
REAL MEMNL4(MDC,MSC,MCGRD)
REAL SA1(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
REAL SA2(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
REAL SFNL(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
REAL SPCSIG(MSC)
REAL UE(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
!
! 7. Common blocks used
!
!
! 8. Subroutines used
!
! ---
!
! 9. Subroutines calling
!
! SOURCE (in SWANCOM1)
!
! 12. Structure
!
! -------------------------------------------
! Initialisations.
! Calculate proportionality constant.
! Prepare auxiliary spectrum.
! Calculate (unfolded) interactions :
! -----------------------------------------
! Energy at interacting bins
! Contribution to interactions
! Fold interactions to side angles
! -----------------------------------------
! Put source term together
! -------------------------------------------
!
! 13. Source text
!
!*******************************************************************
!
INTEGER IS ,ID ,ID0 ,I ,J ,
& ISLOW ,ISHGH ,IDLOW ,IDHGH ,
& IDPP ,IDMM ,ISPP ,ISMM ,
& ISCLW ,ISCHG ,IDDUM ,IENT
!
REAL X ,X2 ,CONS ,FACTOR ,SNLCS1 ,SNLCS2 ,
& SNLCS3 ,E00 ,EP1 ,EM1 ,EP2 ,EM2 ,
& SA1A ,SA1B ,SA2A ,SA2B ,
& JACOBI ,SIGPI
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'SWSNL8')
!
ISLOW = WWINT(9)
ISHGH = WWINT(10)
ISCLW = WWINT(11)
ISCHG = WWINT(12)
IDLOW = WWINT(13)
IDHGH = WWINT(14)
IDPP = WWINT(21)
IDMM = WWINT(22)
ISPP = WWINT(23)
ISMM = WWINT(24)
!
! *** Calculate prop. constant. ***
! *** Calculate factor R(X) to calculate the NL wave-wave ***
! *** interaction for shallow water ***
! *** SNLC1 = 1/GRAV**4 ***
!
SNLCS1 = PQUAD(3)
SNLCS2 = PQUAD(4)
SNLCS3 = PQUAD(5)
X = MAX ( 0.75 * DEP2(KCGRD(1)) * KMESPC , 0.5 )
X2 = MAX ( -1.E15, SNLCS3*X)
CONS = SNLC1 * ( 1. + SNLCS1/X * (1.-SNLCS2*X) * EXP(X2))
JACOBI = 2. * PI
!
! *** extend the area with action density at periodic boundaries ***
!
DO IDDUM = IDLOW, IDHGH
ID = MOD ( IDDUM - 1 + MDC , MDC ) + 1
DO IS=1, MSC
UE (IS,IDDUM) = AC2(ID,IS,KCGRD(1)) * SPCSIG(IS) * JACOBI
ENDDO
ENDDO
!
DO ID = IDLOW, IDHGH
DO IS = MSC+1, ISHGH
UE(IS,ID) = UE(IS-1,ID) * FACHFR
ENDDO
ENDDO
!
! *** Calculate (unfolded) interactions ***
! *** Energy at interacting bins ***
!
DO ID = 1, MDC
DO IS = ISCLW, ISCHG
E00 = UE(IS ,ID )
EP1 = UE(IS+ISPP,ID+IDPP)
EM1 = UE(IS+ISMM,ID-IDMM)
EP2 = UE(IS+ISPP,ID-IDPP)
EM2 = UE(IS+ISMM,ID+IDMM)
!
! Contribution to interactions
!
FACTOR = CONS * AF11(IS) * PQUAD(2) * E00
!
SA1A = E00 * ( EP1*DAL1 + EM1*DAL2 )
SA1B = SA1A - EP1*EM1*DAL3
SA2A = E00 * ( EP2*DAL1 + EM2*DAL2 )
SA2B = SA2A - EP2*EM2*DAL3
!
SA1 (IS,ID) = FACTOR * SA1B
SA2 (IS,ID) = FACTOR * SA2B
!
ENDDO
ENDDO
!
! *** Fold interactions to side angles -> domain in theta is ***
! *** periodic ***
!
DO ID = 1, IDHGH - MDC
ID0 = 1 - ID
DO IS = ISCLW, ISCHG
SA1 (IS,MDC+ID) = SA1 (IS, ID )
SA2 (IS,MDC+ID) = SA2 (IS, ID )
SA1 (IS, ID0 ) = SA1 (IS,MDC+ID0)
SA2 (IS, ID0 ) = SA2 (IS,MDC+ID0)
ENDDO
ENDDO
!
! *** Put source term together (To save space I=IS and ***
! *** J=MDC is used) ---- ***
!
DO I = 1, MSC
SIGPI = SPCSIG(I) * JACOBI
DO J = 1, MDC
SFNL(I,J) = - 2. * ( SA1(I,J) + SA2(I,J) )
& + ( SA1(I-ISPP,J-IDPP) + SA2(I-ISPP,J+IDPP) )
& + ( SA1(I-ISMM,J+IDMM) + SA2(I-ISMM,J-IDMM) )
!
! *** store value in auxiliary array and use values in ***
! *** next four sweeps (see subroutine FILNL3) ***
!
MEMNL4(J,I,KCGRD(1)) = SFNL(I,J) / SIGPI
ENDDO
ENDDO
!
! *** value source term in every bin ***
!
IF ( ITEST.GE.150 .AND. TESTFL ) THEN
DO I=1, MSC
DO J=1, MDC
WRITE(PRINTF,2006) I,J,MEMNL4(J,I,KCGRD(1)),SFNL(I,J),
& SPCSIG(I)
2006 FORMAT (' I J MEMNL() SFNL() SPCSIG:',2I4,3E12.4)
ENDDO
ENDDO
END IF
!
RETURN
!
END SUBROUTINE SWSNL8
!
!*******************************************************************
!
SUBROUTINE FILNL3 (IDCMIN ,IDCMAX ,IMATRA ,IMATDA ,AC2 , 40.23
& MEMNL4 ,PLNL4S ,ISSTOP ,REDC0 ,REDC1 ) 40.85 40.41
!
!*******************************************************************
!
USE SWCOMM3 40.41
USE SWCOMM4 40.41
USE OCPCOMM4 40.41
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 40.23: Marcel Zijlema
! 40.41: Marcel Zijlema
! 40.85: Marcel Zijlema
!
! 1. Updates
!
! 40.23, Aug. 02: rhs and main diagonal adjusted according to Patankar-rules
! 40.41, Oct. 04: common blocks replaced by modules, include files removed
! 40.85, Aug. 08: store quadruplets for output purposes
!
! 2. Purpose
!
! Fill the IMATRA/IMATDA arrays with the nonlinear wave-wave interaction
! source term for a gridpoint ix,iy per sweep direction
!
! 3. Method
!
!
! 4. Argument variables
!
! 7. Common blocks used
!
!
! 8. Subroutines used
!
! ---
!
! 9. Subroutines calling
!
! SOURCE (in SWANCOM1)
!
! 12. Structure
!
! -------------------------------------------
! Do for every frequency and spectral direction within a sweep
! fill IMATRA/IMATDA
! -------------------------------------------
! End of FILNL3
! -------------------------------------------
!
! 13. Source text
!
!*******************************************************************
!
INTEGER IS ,ID ,ISSTOP
!
REAL IMATRA(MDC,MSC) ,
& IMATDA(MDC,MSC) , 40.23
& AC2(MDC,MSC,MCGRD) , 40.23
& PLNL4S(MDC,MSC,NPTST) , 40.00
& MEMNL4(MDC,MSC,MCGRD) 30.21
!
INTEGER IDCMIN(MSC) ,
& IDCMAX(MSC)
REAL :: REDC0 (MDC,MSC,MREDS) 40.85
REAL :: REDC1 (MDC,MSC,MREDS) 40.85
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'FILNL3')
!
DO 990 IS=1, ISSTOP
DO 980 IDDUM = IDCMIN(IS), IDCMAX(IS)
ID = MOD ( IDDUM - 1 + MDC , MDC ) + 1
IF(TESTFL) PLNL4S(ID,IS,IPTST) = MEMNL4(ID,IS,KCGRD(1)) 40.00
IF (MEMNL4(ID,IS,KCGRD(1)).GT.0.) THEN 40.23
IMATRA(ID,IS) = IMATRA(ID,IS) + MEMNL4(ID,IS,KCGRD(1))
REDC0(ID,IS,1)= REDC0(ID,IS,1)+ MEMNL4(ID,IS,KCGRD(1)) 40.85
ELSE
IMATDA(ID,IS) = IMATDA(ID,IS) - MEMNL4(ID,IS,KCGRD(1)) /
& MAX(1.E-18,AC2(ID,IS,KCGRD(1)))
REDC1(ID,IS,1)= REDC1(ID,IS,1)+ MEMNL4(ID,IS,KCGRD(1)) / 40.85
& MAX(1.E-18,AC2(ID,IS,KCGRD(1))) 40.85
END IF
980 CONTINUE
990 CONTINUE
!
IF ( TESTFL .AND. ITEST.GE.50 ) THEN
WRITE(PRINTF,9000) IDCMIN(1),IDCMAX(1),MSC,ISSTOP
9000 FORMAT(' FILNL3: ID_MIN ID_MAX MSC ISTOP :',4I6)
IF ( ITEST .GE. 100 ) THEN
DO IS=1, ISSTOP
DO IDDUM = IDCMIN(IS), IDCMAX(IS)
ID = MOD ( IDDUM - 1 + MDC , MDC ) + 1
WRITE(PRINTF,6001) IS,ID,MEMNL4(ID,IS,KCGRD(1)) 30.21
6001 FORMAT(' FILNL3: IS ID MEMNL() :',2I6,E12.4)
ENDDO
ENDDO
ENDIF
ENDIF
!
RETURN
! End of FILNL3
END
!
!----------------------------------------------------------------------
SUBROUTINE SWINTFXNL ( ASWAN,SIGMA,DIR,NDIR,NSIG,NGRID,DEPTH,
& IQTYPE,SNL,KCGRD,ICMAX,IERROR )
!----------------------------------------------------------------------
!
! +-------+ ALKYON Hydraulic Consultancy & Research
! | | Gerbrant Ph. van Vledder
! | +---+
! | | +---+ Last update: 9 September 2002
! +---+ | | Release: 5.0
! +---+
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
USE M_PARALL
USE serv_xnl4v5
USE m_xnldata
!
IMPLICIT NONE
!-----------------------------------------------------------------------------------
!
! 0. Update history
!
! Date Modification
!
! 25/02/1999 Initial version
! 16/03/1999 Interface with SWAN updated
! 24/03/1999 Interface extended with KCGRD and ICMAX
! 22/06/1999 Parameter IERR added to interface
! 20/07/1999 Call to Q_QMAIN modified
! 25/07/1999 Output files updated
! 24/09/1999 Finite depth effects included
! 25/11/1999 Bug fixed in deleting files
! 27/12/1999 Interface extended with grav,rho and ftail
! 02/02/2001 Interface modified, was N(k), now A(sigma) like SWAN
! 06/11/2001 Bug fixed in initialisation of Snl
! 08/08/2002 Version 4
! 22/08/2002 Size of direction array modified to conform with SWAN 40.11
! 09/09/2002 Release 5
! 18/05/2004 Implemented in SWAN 40.41, with adapted values of IQTYPE
!
!
! 1. Purpose:
!
! interface with SWAN model to compute nonlinear transfer with
! the XNL method for given action density spectrum
!
! 2. Method
!
! Resio/Tracy deep water geometric scaling
! Rewritten by Gerbrant van Vledder
!
! 3. Parameter list:
!
! Type I/O Name Description
!----------------------------------------------------------------------------------
INTEGER, INTENT(IN) :: NDIR ! number of directions
INTEGER, INTENT(IN) :: NSIG ! number of sigma values in array
INTEGER, INTENT(IN) :: NGRID ! number of sea points in the comp. grid
INTEGER, INTENT(IN) :: IQTYPE ! method of computing nonlinear quadruplet
! interactions
REAL , INTENT(IN) :: ASWAN(NDIR,NSIG,NGRID) ! action density spectrum as a
! ! function of (sigma,dir)
REAL , INTENT(IN) :: SIGMA(NSIG) ! Intrinsic frequencies
REAL , INTENT(IN) :: DIR(NDIR,6) ! directions in radians (when second index =1)
REAL , INTENT(IN) :: DEPTH(NGRID) ! depth array
INTEGER, INTENT(IN) :: ICMAX ! number of points of the computational stencil
INTEGER, INTENT(IN) :: KCGRD(ICMAX) ! grid addresses for points of computational stencil
REAL , INTENT(OUT):: SNL(NDIR,NSIG,NGRID) ! nonlinear quadruplet interaction computed with
! ! a certain exact method (sigma,dir)
INTEGER, INTENT(OUT):: IERROR ! Error indicator. If no errors are detected IERR=0
!--------------------------------------------------------------------------------------------------
!
! 4. Error messages
!
! An error message is produced within the QUAD system.
! If no errors are detected IERROR=0
! 1, incorrect IQUAD
! 2, depth < 0
!
! 5. Subroutines calling
!
! SOURCE
!
! 6. Subroutines used
!
! 7. Remarks
!
! The SWAN spectrum is given as an action density spectrum
! as a function of Sigma and Theta: ASWAN(itheta,isig)
!
! SWINTFXNL is called for each active grid point in a stencil
! and for each time the complete array with all grid points
! is given. Related grid points are specified in the array
! KCGRD and ICMAX.
!
! 8. Structure
!
! 9. Switches
!
! 10. Source code
!------------------------------------------------------------------------------
! Local parameters
!
INTEGER ISIG,IDIR ! counters
INTEGER IGRID ! grid index
INTEGER IQUAD ! type of computational method for Xnl
!
! --- assign arrays for intermediate storage of results
!
REAL AQUAD(NSIG,NDIR) ! action density spectrum A(sigma,dir)
REAL XNL(NSIG,NDIR) ! transfer rate dA/dt(sigma,dir)
REAL DIAG(NSIG,NDIR) ! diagonal term (dXnl/dA)
REAL DIRR(NDIR) ! single array with directions in radians
!------------------------------------------------------------------------------
!
! --- initialisations
!
IERROR = 0
!
IGRID = KCGRD(1) ! set index of current grid index
DIRR(:) = DIR(:,1) ! copy radian directions to single array
!
SNL(:,:,IGRID) = 0.
DIAG = 0.
!
! --- check value of iquad
!
IF ((IQTYPE.NE.51).AND.(IQTYPE.NE.52).AND.(IQTYPE.NE.53)) THEN
IERROR = 1
GOTO 9999
END IF
!
! N.B. Take care of different order of indices in QUAD compared to SWAN
!
! --- compute nonlinear interactions per individual spectrum
! switch order of indices
!
DO ISIG = 1, NSIG
DO IDIR = 1, NDIR
AQUAD(ISIG,IDIR) = ASWAN(IDIR,ISIG,IGRID)
END DO
END DO
!
! --- transform parameter iquad of SWAN to the parameter iq_quad as
! needed by the QUAD suite
!
IQUAD = IQTYPE - 50
!
! IQTYPE/IQUAD = 51/1 ! deep water transfer
! IQTYPE/IQUAD = 52/2 ! deep water transfer with WAM depth scaling
! IQTYPE/IQUAD = 53/3 ! finite depth transfer
!
! --- call of main subroutine to compute nonlinear quadruplet interactions
! for a given action density spectrum on a given spectral grid
!
XNL = 0.
DIAG = 0.
CALL XNL_MAIN( AQUAD,SIGMA,DIRR,NSIG,NDIR,DEPTH(IGRID),IQUAD,
& XNL,DIAG,INODE,IERROR )
!
IF (IERROR.NE.0) GOTO 9999
!
! --- convert nonlinear transfer to SWAN convention, only sequence of indices
!
DO ISIG = 1, NSIG
DO IDIR = 1, NDIR
SNL(IDIR,ISIG,IGRID) = XNL(ISIG,IDIR)
END DO
END DO
!
9999 CONTINUE
!
RETURN
!
END SUBROUTINE SWINTFXNL
!
!********************************************************************
!
SUBROUTINE TCOEF (W1, W2, W12, K1, K2, K12, DEP2, R, S, ICC)
!
!********************************************************************
!
USE SWCOMM3
IMPLICIT NONE
!
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
! 41.46: James Salmon
! 1. Updates
! 41.46, October 2013: new subroutine
! 2. Purpose
!
! Calculates coupling coefficients
!
! 3. Method
!
! R_m,p-m = (k_m + k_p-m)^2 * [0.5 + (w_m*w_p-m / g*h*k_m*k_p-m)]
!
! S_p = -2/g * ( g*h*k_p + 2*B*g*h^3*k_p^3 - (B+(1/3))*h^2*w_p^2*k_p )
!
! where: B = 1/15
!
! 4. Argument variables
!
! ICC type of coupling coefficient
! [1] - Boussinesq: Freilich & Guza (1984), Herbers & Burton (1997)
! [2] - Deterministic Boussinesq: Madsen & Sorensen (1993); suggested for SPB
!
! W1, W2, W12 w_p, w_m, w_I
! K1, K2, K12 k_p, k_m, k_I where I represents the sum or difference
! DEP2 water depth
!
INTEGER ICC
!
REAL :: W1, W2, W12
REAL :: K1, K2, K12
!
REAL :: DEP2(MCGRD)
!
REAL, INTENT(OUT) :: R
REAL, INTENT(OUT) :: S
!
! Values from common
!
! 6. Local variables
!
REAL :: DEP , DEP_2, DEP_3
REAL :: B , B2 , B3
REAL :: KPROD
!
! 7. SUBROUTINES USED
! 8. SUBROUTINES CALLING
! SWSPB
! 9. ERROR MESSAGES
! ---
! 10. REMARKS
! ---
! 11. STRUCTURE
! -----------------------------------------------------------------
! -----------------------------------------------------------------
! 13. Source text
!
DEP = DEP2(KCGRD(1))
DEP_2 = DEP**2
DEP_3 = DEP**3
!
B = 1./15.
B2 = 2.*B
B3 = (B + (1./3.))
!
IF (ICC .EQ. 1) THEN !Boussinesq
R = 0.75 * (W2 + W12)
S = -DEP * SQRT(GRAV*DEP)
!
ELSEIF (ICC .EQ. 2) THEN !Deterministic Boussinesq
! to avoid NaN when W2*W12=0=K2*K12, adapt product of K2 and K12
KPROD = SIGN(MAX(1.E-20, ABS(K2*K12)),K2*K12)
R = (0.5 + ((W2*W12)/(GRAV*DEP*KPROD))) * (K2 + K12)**2
!
S = (-2./GRAV) * ( (GRAV*DEP*K1)
& + (B2*GRAV*DEP_3*K1**3)
& - (B3*DEP_2*K1*W1**2) )
!
ENDIF
!
RETURN
END subroutine TCOEF
!
!********************************************************************
!
SUBROUTINE COAPP (AC2, SPCSIG, DEP2, IDCMIN, IDCMAX, INTPL, SX,
& KX , CX, CGX, E , EI , EX , EIX , WX)
!
!********************************************************************
!
USE OCPCOMM4
USE SWCOMM3
USE SWCOMM4
IMPLICIT NONE
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
! 41.46: James Salmon
! 1. Updates
! 41.46, Nov 2014: new subroutine
! 2. Purpose
!
! Pre-calculate directional integrations for symmetric subst.
! Interpolates energy densities and integrated energy densities at
! interacting frequencies
! Calculates corresponding wave numbers and phase velocities
!
! 3. Method
!
! Use technique from SWLTA that due to logarithmic distribution, a given factor between two freq.
! is the same number of freq. bins
!
! 4. Argument variables
!
! AC2 action density
! SPCSIG relative frequencies in computational domain in sigma-space
! IDCMIN minimum counter in directional space
! IDCMAX maximum counter in directional space
! INTPL Switch to control the interpolation
! 0 = use lumping i.e. sum at half; dif. at double
! 1 = use correct terms i.e. sum with p-m; dif. with p+m
!
REAL :: AC2(MDC,MSC,MCGRD)
REAL :: SPCSIG(MSC)
REAL :: DEP2(MCGRD)
INTEGER :: IDCMIN(MSC), IDCMAX(MSC)
INTEGER :: INTPL
!
! SX relative frequencies at sum and difference frequencies
! KX wave number at sum and difference frequencies
! CGX group velocity at sum and difference velocities
! CX phase velocity at sum and difference velocities
! E energy density as function of frequency and direction
! EI integrated E over PWDTH as funct. of dir. and freq..
! EX E at sum and difference frequencies SX
! EIX EI densities at sum and dif. freq. SX
! WX weights for logarithmic interpolation
!
REAL, INTENT(OUT) :: SX(MSC,MSC,2)
REAL, INTENT(OUT) :: KX(MSC,MSC,2)
REAL, INTENT(OUT) :: CGX(MSC,MSC,2)
REAL, INTENT(OUT) :: CX(MSC,MSC,2)
REAL, INTENT(OUT) :: E(MDC,MSC)
REAL, INTENT(OUT) :: EI(MDC,MSC)
REAL, INTENT(OUT) :: EX(MDC,MSC,MSC,2)
REAL, INTENT(OUT) :: EIX(MDC,MSC,MSC,2)
REAL, INTENT(OUT) :: WX(MSC,MSC,2)
!
! Values from common
!
! MDC : Size of array in theta-direction
! MSC : Size of array in sigma-direction
! PI : Circular constant Pi
! PTRIAD(8) : range to integrate over (in degrees)
!
! 6. Local variables
!
! B : interger inteval of bins for interpolation
! DUMMY : dummy variable
! FX : factors for interpolation to sum and dif. freqencies
! I1 : auxiliary integer
! I2 : auxiliary integer
! ID : counter
! II : loop counter
! IS : loop counter in frequency space
! PAVL : integration over PINT using sector
! PEXC : fraction of excess at outermost directional bins
! PFAC : correction using PEXC using sector
! PINT : outermost directional bin for integration
! PLIM : relative indices for directional bin limits
! PSUM : sum of E(tht,f) over PWDTH
! PWDTH : integral range in rad. / 2
! XIS : rate between two succeeding frequency counters
! XISLN : log of XIS
!
INTEGER ID, II, IS, I1, I2
INTEGER PINT, PLIM(2), B(2)
REAL PEXC, PSUM , PWDTH , PAVL(2), PFAC(2)
REAL XIS , XISLN
REAL FX(MSC,MSC,2)
REAL D , TMP(4)
!
!
! 7. SUBROUTINES USED
! 8. SUBROUTINES CALLING
! SWSPB
! 9. ERROR MESSAGES
! ---
! 10. REMARKS
! ---
! 11. STRUCTURE
! -----------------------------------------------------------------
! -----------------------------------------------------------------
! 13. Source text
!
D = DEP2(KCGRD(1))
!
! --- convert from AC(tht,f) to E(tht,f)
!
E = 0.
DO II = 1, MDC
DO IS = 1, MSC
E(II,IS) = AC2(II,IS,KCGRD(1)) * 2. * PI * SPCSIG(IS)
ENDDO
ENDDO
!
! --- integrate over range dir0-(p) <= tht < dir+(p)
!
EI = 0.
PWDTH = PTRIAD(8) * PI/180.
!IF (PWDTH.EQ.0. .OR. PWDTH.GE.DDIR*MDC) THEN
IF (.NOT.PWDTH.NE.0.) THEN
! --- Full circle integration
DO IS = 1, MSC
DO II = IDCMIN(IS), IDCMAX(IS)
ID = MOD(II - 1 + MDC, MDC) + 1
EI(ID,IS) = SUM(E(:,IS), DIM=1) * DDIR
ENDDO
ENDDO
ELSE
! --- integrate over range PWDTH
PINT = INT((PWDTH/DDIR) - 0.5) + 1
PEXC = (((PINT + 0.5) * DDIR) - PWDTH) / DDIR
IF (FULCIR) THEN
! --- integration over circle
DO IS = 1, MSC
DO II = IDCMIN(IS), IDCMAX(IS)
ID = MOD(II - 1 + MDC, MDC) + 1
PSUM = 0.
PLIM = 0
! --- define integration limits
PLIM(1) = MOD(ID - PINT + MDC - 1, MDC) + 1
PLIM(2) = MAX(MOD(ID + PINT , MDC), 1)
! --- sum values by full int.&subtract or zero and adding
PSUM = SUM(E(:,IS), DIM=1)
& * MIN(MAX(REAL(PLIM(1) - PLIM(2)), 0.), 1.)
PSUM = PSUM - SUM(E(1:PLIM(1),IS), DIM=1)
& + SUM(E(1:PLIM(2),IS), DIM=1)
PSUM = PSUM - ( (E(PLIM(1),IS)
& + E(PLIM(2),IS)) * PEXC)
!
EI(ID,IS) = PSUM * DDIR
ENDDO
ENDDO
ELSE
! --- integration over a sector
DO IS = 1, MSC
DO II = IDCMIN(IS), IDCMAX(IS)
ID = MOD(II - 1 + MDC, MDC) + 1
PSUM = 0.
PLIM = 0
! --- define integration limits
PLIM(1) = MAX(ID - PINT, 1)
PLIM(2) = MIN(ID + PINT, MDC)
PSUM = SUM(E(ID-PINT:ID+PINT,IS), DIM=1)
! --- determine range covered
PAVL(1) = (ID - 0.5) * DDIR
PAVL(2) = (MDC - ID + 0.5) * DDIR
PFAC(1) = MIN(MAX(PAVL(1) - PWDTH, 0.), 1.) * PEXC
PFAC(2) = MIN(MAX(PAVL(2) - PWDTH, 0.), 1.) * PEXC
! --- remove excess
PSUM = PSUM - ( E(PLIM(1),IS) * PFAC(1)
& +E(PLIM(2),IS) * PFAC(2))
!
EI(ID,IS) = PSUM*DDIR
ENDDO
ENDDO
ENDIF
ENDIF
!
! --- compute some indices in sigma space
! --- generate factors for logarithmic interpolation
! --- interpolate E and EI
!
I2 = INT(FLOAT(MSC) / 2.)
I1 = I2-1
XIS = SPCSIG(I2) / SPCSIG(I1) !=1 + df/f
XISLN = LOG(XIS)
!
FX = 0.
SX = 0.
KX = 0.
CGX = 0.
CX = 0.
EX = 0.
EIX = 0.
WX = 0.
IF (INTPL .EQ. 0) THEN
! --- compute only self-self interactions
! --- factors for lumping
FX(:,:,1) = 0.5
FX(:,:,2) = 2.
DO IS = 1, MSC
DO ID = 1, 2
! --- determine # freq. bins to either side of target
TMP = 0.
B = 0
B(1) = INT(LOG(FX(IS,1,ID)) / XISLN)
B(2) = B(1) + SIGN(1, B(1))
! --- determine interpolation weight
WX(IS,1,ID) = (FX(IS,1,ID) - XIS**B(2))
& / (XIS**B(1) - XIS**B(2))
! --- interpolation
! --- !Check if diff. freq. >= min(freq.) and sum < 2*max*freq(
IF (IS .GT. -B(2) .AND. MSC .GT. IS+B(2)+1) THEN
SX(IS,:,ID) = (SPCSIG(IS+B(2))*(1.-WX(IS,1,ID)))
& + (SPCSIG(IS+B(1))* WX(IS,1,ID) )
TMP(1) = SX(IS,1,ID)
CALL KSCIP1(1,TMP(1),D,TMP(2),TMP(3),TMP(4),TMP(4))
KX(IS,:,ID) = TMP(2)
CGX(IS,:,ID) = TMP(3)
CX(IS,:,ID) = SX(IS,:,ID) / KX(IS,:,ID)
EX(:,IS,IS,ID) = (E(:,IS+B(2))*(1.-WX(IS,1,ID)))
& + (E(:,IS+B(1)) * WX(IS,1,ID) )
EIX(:,IS,IS,ID) = (EI(:,IS+B(2))*(1.-WX(IS,1,ID)))
& + (EI(:,IS+B(1)) * WX(IS,1,ID) )
ELSE
SX(IS,:,ID) = ( (SPCSIG(IS)*XIS**B(2)) *
& (1.-WX(IS,II,ID)) )
& + ( (SPCSIG(IS)*XIS**B(1)) *
& (WX(IS,II,ID)) )
TMP(1) = SX(IS,1,ID)
CALL KSCIP1(1,TMP(1),D,TMP(2),TMP(3),TMP(4),TMP(4))
KX(IS,:,ID) = TMP(2)
CGX(IS,:,ID) = TMP(3)
CX(IS,:,ID) = SX(IS,1,ID) / KX(IS,1,ID)
EX(:,IS,IS,ID) = 0.
EIX(:,IS,IS,ID) = 0.
ENDIF
ENDDO
ENDDO
ELSEIF (INTPL .EQ. 1) THEN
! --- compute all interacting frequencies
DO IS = 1, MSC
DO II = 1, MSC
! --- factors for all sum and diff. interactions
FX(IS,II,1) = MAX((SPCSIG(IS) - SPCSIG(II))
& / SPCSIG(IS), 1.E-20)
FX(IS,II,2) = (SPCSIG(IS) + SPCSIG(II)) / SPCSIG(IS)
!
DO ID = 1, 2
! --- determine # freq. bins to either side of target
TMP = 0.
B = 0
B(1) = INT(LOG(FX(IS,II,ID)) / XISLN)
B(2) = B(1) + SIGN(1, B(1))
! --- determine interpolation weight
WX(IS,II,ID) = (FX(IS,II,ID) - XIS**B(2))
& / (XIS**B(1) - XIS**B(2))
! --- interpolation
! --- !Check if diff. freq. >= min(freq.) and sum < 2*max*freq(
IF (IS .GT. -B(2) .AND. MSC .GT. IS+B(2)+1) THEN
SX(IS,II,ID) = (SPCSIG(IS+B(2))*(1.-WX(IS,II,ID)))
& + (SPCSIG(IS+B(1))* WX(IS,II,ID) )
TMP(1) = SX(IS,II,ID)
CALL KSCIP1(1,TMP(1),D,TMP(2),TMP(3),TMP(4),TMP(4))
KX(IS,II,ID) = TMP(2)
CGX(IS,II,ID) = TMP(3)
CX(IS,II,ID) = SX(IS,II,ID) / KX(IS,II,ID)
EX(:,IS,II,ID) = (E(:,IS+B(2)) *(1.-WX(IS,II,ID)))
& + (E(:,IS+B(1)) * WX(IS,II,ID))
EIX(:,IS,II,ID) = (EI(:,IS+B(2))*(1.-WX(IS,II,ID)))
& + (EI(:,IS+B(1))* WX(IS,II,ID))
ELSE
SX(IS,II,ID) = ( (SPCSIG(IS)*XIS**B(2)) *
& (1.-WX(IS,II,ID)) )
& + ( (SPCSIG(IS)*XIS**B(1)) *
& (WX(IS,II,ID)) )
TMP(1) = SX(IS,II,ID)
CALL KSCIP1(1,TMP(1),D,TMP(2),TMP(3),TMP(4),TMP(4))
KX(IS,II,ID) = TMP(2)
CGX(IS,II,ID) = TMP(3)
CX(IS,II,ID) = SX(IS,II,ID) / KX(IS,II,ID)
EX(:,IS,II,ID) = 0.
EIX(:,IS,II,ID) = 0.
ENDIF
ENDDO
ENDDO
ENDDO
ENDIF
!
RETURN
END subroutine COAPP
!
!****************************************************************
!
SUBROUTINE SWLTA ( AC2 , DEP2 , CGO , SPCSIG,
& KWAVE , IMATRA, IMATDA, REDC0 , REDC1 ,
& IDDLOW, IDDTOP, ISSTOP, IDCMIN, IDCMAX,
& HS , SMEBRK, PLTRI , URSELL )
!
!****************************************************************
!
USE OCPCOMM4
USE SWCOMM3
USE SWCOMM4
!
IMPLICIT NONE
!
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
!
! 40.56: Marcel Zijlema
! 40.85: Marcel Zijlema
! 41.44: James Salmon, Pieter Smit
!
! 1. Updates
!
! 40.56, Feb. 06: New subroutine
! 40.85, Aug. 08: store triads for output purposes
! 41.44, Oct. 13: consistency fix
!
! 2. Purpose
!
! In this subroutine the triad-wave interactions are calculated
! with the Lumped Triad Approximation of Eldeberky (1996). His
! expression is based on a parametrization of the biphase (as
! function of the Ursell number), is directionally uncoupled and
! takes into account for the self-self interactions only.
!
! For a full description of the equations reference is made
! to PhD thesis of Eldeberky (1996). Here only the main expressions
! are given.
!
! 3. Method
!
! The parametrized biphase is given by (see eq. 3.19):
!
! 0.2
! beta = - pi/2 + pi/2 tanh ( ----- )
! Ur
!
! The Ursell number is calculated in routine SINTGRL
!
! The source term as function of frequency p is (see eq. 7.25):
!
! + -
! S(p) = S(p) + S(p)
!
! in which
!
! +
! S(p) = alpha Cp Cg,p (R(p/2,p/2))**2 sin (|beta|) ( E(p/2)**2 -2 E(p) E(p/2) )
!
! - +
! S(p) = - 2 S(2p)
!
! with alpha a tunable coefficient and R(p/2,p/2) is the interaction
! coefficient of which the expression can be found in Eldeberky (1996);
! see eq. 7.26.
!
! Note that a slightly adapted formulation of the LTA is used in
! in the SWAN model:
!
! - Only positive contributions to higher harmonics are considered
! here (no energy is transferred to lower harmonics).
!
! - The mean frequency in the expression of the Ursell number
! is calculated according to the first order moment over the
! zeroth order moment (personal communication, Y.Eldeberky, 1997).
!
! - The interactions are calculated up to 2.5 times the mean
! frequency only.
!
! - Since the spectral grid is logarithmically distributed in frequency
! space, the interactions between central bin and interacting bin
! are interpolated such that the distance between these bins is
! factor 2 (nearly).
!
! - The interactions are calculated in terms of energy density
! instead of action density. So the action density spectrum
! is firstly converted to the energy density grid, then the
! interactions are calculated and then the spectrum is converted
! to the action density spectrum back.
!
! - To ensure numerical stability the Patankar rule is used.
!
! 4. Argument variables
!
! AC2 action density
! CGO group velocity
! DEP2 water depth
! HS significant wave height
! IDCMIN minimum counter in directional space
! IDCMAX maximum counter in directional space
! IDDLOW minimum direction that is propagated within a sweep
! IDDTOP maximum direction that is propagated within a sweep
! IMATDA main diagonal of the linear system
! IMATRA right-hand side of system of equations
! ISSTOP maximum frequency counter in a sweep
! KWAVE wave number
! PLTRI triad contribution in TEST points
! SMEBRK average (angular) frequency
! SPCSIG relative frequencies in computational domain in sigma-space
! URSELL Ursell number
!
INTEGER IDDLOW, IDDTOP, ISSTOP
INTEGER IDCMIN(MSC), IDCMAX(MSC)
REAL :: HS, SMEBRK
REAL :: AC2(MDC,MSC,MCGRD)
REAL :: CGO(MSC,MICMAX)
REAL :: DEP2(MCGRD)
REAL :: IMATDA(MDC,MSC), IMATRA(MDC,MSC)
REAL :: SPCSIG(MSC)
REAL :: KWAVE(MSC,MICMAX)
REAL :: PLTRI(MDC,MSC,NPTST)
REAL :: URSELL(MCGRD)
REAL :: REDC0 (MDC,MSC,MREDS) 40.85
REAL :: REDC1 (MDC,MSC,MREDS) 40.85
!
! 6. Local variables
!
! AUX1 : auxiliary real
! AUX2 : auxiliary real
! BIPH : parameterized biphase of the spectrum
! C0 : phase velocity at central bin
! CM : phase velocity at interacting bin
! DEP : water depth
! DEP_2 : water depth to power 2
! DEP_3 : water depth to power 3
! E : energy density as function of frequency
! E0 : energy density at central bin
! ED : integral energy density over directions 41.44
! ED0 : integral energy density over directions at central bin 41.44
! EDM : integral energy density over directions at 41.44
! interacting bin
! EM : energy density at interacting bin
! FT : auxiliary real indicating multiplication factor
! for triad contribution
! I1 : auxiliary integer
! I2 : auxiliary integer
! ID : counter
! IDDUM : loop counter in direction space
! IENT : number of entries
! II : loop counter
! IS : loop counter in frequency space
! ISM : negative range for IS
! ISM1 : negative range for IS
! ISMAX : maximum of the counter in frequency space for
! which the triad interactions are calculated (cut-off)
! ISP : positive range for IS
! ISP1 : positive range for IS
! RINT : interaction coefficient
! SA : interaction contribution of triad
! SIGPI : frequency times 2pi
! SINBPH: absolute sine of biphase
! STRI : total triad contribution
! WISM : interpolation weight factor corresponding to lower harmonic
! WISM1 : interpolation weight factor corresponding to lower harmonic
! WISP : interpolation weight factor corresponding to higher harmonic
! WISP1 : interpolation weight factor corresponding to higher harmonic
! W0 : radian frequency of central bin
! WM : radian frequency of interacting bin
! WN0 : wave number at central bin
! WNM : wave number at interacting bin
! XIS : rate between two succeeding frequency counters
! XISLN : log of XIS
!
INTEGER I1, I2, ID, IDDUM, IENT, II, IS, ISM, ISM1, ISMAX,
& ISP, ISP1
REAL AUX1, AUX2, BIPH, C0, CM, DEP, DEP_2, DEP_3, E0, ED0, EDM,
& EM, FT, RINT, SIGPI, SINBPH, STRI, WISM, WISM1, WISP,
& WISP1, W0, WM, WN0, WNM, XIS, XISLN
REAL, ALLOCATABLE :: E(:), ED(:), SA(:,:)
!
! 9. Subroutines calling
!
! SOURCE
!
! 12. Structure
!
! Determine resonance condition and the maximum discrete freq.
! for which the interactions are calculated.
!
! If Ursell number larger than prescribed value compute interactions
! Calculate biphase
! Do for each direction
! Convert action density to energy density
! Do for all frequencies
! Calculate interaction coefficient and interaction factor
! Compute interactions and store results in matrix
!
! 13. Source text
!
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'SWLTA')
DEP = DEP2(KCGRD(1))
DEP_2 = DEP**2
DEP_3 = DEP**3
!
! --- compute some indices in sigma space
!
I2 = INT (FLOAT(MSC) / 2.)
I1 = I2 - 1
XIS = SPCSIG(I2) / SPCSIG(I1)
XISLN = LOG( XIS )
ISP = INT( LOG(2.) / XISLN )
ISP1 = ISP + 1
WISP = (2. - XIS**ISP) / (XIS**ISP1 - XIS**ISP)
WISP1 = 1. - WISP
ISM = INT( LOG(0.5) / XISLN )
ISM1 = ISM - 1
WISM = (XIS**ISM -0.5) / (XIS**ISM - XIS**ISM1)
WISM1 = 1. - WISM
ALLOCATE (E (1:MSC))
ALLOCATE (ED(1:MSC))
ALLOCATE (SA(1:MDC,1:MSC+ISP1))
E = 0.
ED = 0.
SA = 0.
!
! --- compute maximum frequency for which interactions are calculated
!
ISMAX = 1
DO IS = 1, MSC
IF ( SPCSIG(IS) .LT. ( PTRIAD(2) * SMEBRK) ) THEN
ISMAX = IS
ENDIF
ENDDO
ISMAX = MAX ( ISMAX , ISP1 )
!
! --- compute 3 wave-wave interactions
!
IF ( URSELL(KCGRD(1)).GE.PTRIAD(5) ) THEN
!
! --- calculate biphase
!
BIPH = (0.5*PI)*(TANH(PTRIAD(4)/URSELL(KCGRD(1)))-1.)
SINBPH = ABS( SIN(BIPH) )
!
! --- calculate integral of energy density over directions
!
ED(:) = SUM(AC2(:,:,KCGRD(1)),DIM=1) * 2.*PI * SPCSIG(:) * DDIR 41.44
!
DO II = IDDLOW, IDDTOP
ID = MOD ( II - 1 + MDC , MDC ) + 1
!
! --- initialize array with E(f) for the direction considered
!
DO IS = 1, MSC
E(IS) = AC2(ID,IS,KCGRD(1)) * 2. * PI * SPCSIG(IS)
END DO
!
DO IS = 1, ISMAX
E0 = E(IS)
ED0 = ED(IS)
W0 = SPCSIG(IS)
WN0 = KWAVE(IS,1)
C0 = W0 / WN0
IF ( IS.GT.-ISM1 ) THEN
EM = WISM * E(IS+ISM1) + WISM1 * E(IS+ISM)
EDM = WISM * ED(IS+ISM1) + WISM1 * ED(IS+ISM)
WM = WISM * SPCSIG(IS+ISM1) + WISM1 * SPCSIG(IS+ISM)
WNM = WISM * KWAVE(IS+ISM1,1) + WISM1 * KWAVE(IS+ISM,1)
CM = WM / WNM
ELSE
EM = 0.
EDM = 0.
WM = 0.
WNM = 0.
CM = 0.
END IF
AUX1 = WNM**2 * ( GRAV * DEP + 2.*CM**2 )
AUX2 = WN0 * DEP * ( GRAV * DEP +
& (2./15.) * GRAV * DEP_3 * WN0**2 -
& (2./ 5.) * W0**2 * DEP_2 )
RINT = AUX1 / AUX2
FT = PTRIAD(1) * C0 * CGO(IS,1) * RINT**2 * SINBPH
IF (ITRIAD.EQ.11) THEN
SA(ID,IS) = MAX(0., FT * ( EM * EM - 2. * EM * E0 ))
ELSE
SA(ID,IS) = MAX(0., FT * ( EDM * (EM - E0) - ED0 * EM )) 41.44
ENDIF
END DO
END DO
!
! --- put source term together
!
DO IS = 1, ISSTOP
SIGPI = SPCSIG(IS) * 2. * PI
DO IDDUM = IDCMIN(IS), IDCMAX(IS)
ID = MOD ( IDDUM - 1 + MDC , MDC ) + 1
!
STRI = SA(ID,IS) - 2.*(WISP * SA(ID,IS+ISP1) +
& WISP1 * SA(ID,IS+ISP ))
!
! --- store results in rhs and main diagonal according
! to Patankar-rules
!
IF(TESTFL) PLTRI(ID,IS,IPTST) = STRI / SIGPI
IF (STRI.GT.0.) THEN
IMATRA(ID,IS) = IMATRA(ID,IS) + STRI / SIGPI
REDC0(ID,IS,2)= REDC0(ID,IS,2)+ STRI / SIGPI 40.85
ELSE
IMATDA(ID,IS) = IMATDA(ID,IS) - STRI /
& MAX(1.E-18,AC2(ID,IS,KCGRD(1))*SIGPI)
REDC1(ID,IS,2)= REDC1(ID,IS,2)+ STRI / 40.85
& MAX(1.E-18,AC2(ID,IS,KCGRD(1))*SIGPI) 40.85
END IF
END DO
END DO
END IF
!
! --- test output
!
IF ( ITEST .GE. 5 .AND. TESTFL ) THEN
WRITE(PRINTF,2000) KCGRD(1), ISMAX
2000 FORMAT (' SWLTA: KCGRD ISMAX :',2I4)
WRITE(PRINTF,2001) GRAV, DEP, DEP_2, DEP_3
2001 FORMAT (' SWLTA: G DEP DEP2 DEP3 :',4E12.4)
WRITE(PRINTF,2002) PTRIAD(1), PTRIAD(2), URSELL(KCGRD(1))
2002 FORMAT (' SWLTA: P(1) P(2) P4) URSELL :',4E12.4)
WRITE(PRINTF,2003) SMEBRK, HS, BIPH, ABS(SIN(BIPH))
2003 FORMAT (' SWLTA: SMEBRK HS B |SIN(B)|:',4E12.4)
END IF
DEALLOCATE(E,ED,SA)
RETURN
END SUBROUTINE SWLTA
!
!********************************************************************
!
SUBROUTINE SWSPB ( AC2 , DEP2 , CGO , SPCSIG,
& KWAVE , IMATRA, IMATDA, REDC0 , REDC1 ,
& IDDLOW, IDDTOP, ISSTOP, IDCMIN, IDCMAX,
& PLTRI )
!
!********************************************************************
!
USE OCPCOMM4
USE SWCOMM3
USE SWCOMM4
IMPLICIT NONE
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
! 41.46: James Salmon
! 1. Updates
! 41.46, Nov 2014: new subroutine
! 2. Purpose
!
! Computes triad source term based on Becq-Girard et al. (1999)
!
! 3. Method
!
! Expression given by Becq-Girard et al. (1999) for spatial evolution
! of the continuous variance density E(fp) of component p:
!
! E(fp)' = <shoaling term> + 4 * <SUM> - 8 * <DIF.>
!
! <SUM> = (1/Sp) * INT{0,fp} Rm,p-m * Im(B(fm,fp-fm)) dfm
! <DIF> = (1/Sp) * INT{0,inf} R-m,p+m * Im(B*(-fm,fp+fm)) dfm
!
! Im(B(fm,fp-fm)) = 0.5 * [ Rm,p-m * E(fm) * E(fp-fm) / Sp
! - Rp,-m * E(fp) * E(fm) / Sp-m
! - Rp,m-p * E(fp) * E(fp-fm) / Sm ] * MU
!
! Im(B*(-fm,fp+fm)) = Im(B(fm,-(fp+fm))) = Im(fm,fp)
! Im(B(fm,fp)) = 0.5 * [ Rm,p * E(fp) * E(fm) / Sp+m
! - Rp+m,-m* E(p+m)* E(fm) / Sp
! - Rp+m,-p* E(p+m)* E(fp) / Sm ] * MU
!
! where Sx and Rx,y are coupling coeffients (see TCOEF)
!
! MU = Im(1/(delta_k - iK)) = K / (delta_k^2 + K^2)
!
! delta_k = k1,2 - k1 - k2 = kp-m + km -kp = kp+m -km - kp
! K = (C1 * k_p) + C2
! where k_p is the peak wave number; defaults C1 = 0.95, C2 = - 0.75
!
! The source term is then given by: cg,p(4*<SUM> - 8*<DIF.>)
! where cg,p is the group velocity of component p, freq = fp
!
! Note(1): in the literature occasionally the following form is given:
!
! (cg,p*Rx,y/Rx+y) * [Ia,b - Ia+b,-a - Ia+b,-b]
!
! where Ia,b = sqrt(cg|a+b|/cg|a|/cg|b|)*cg|a|*cg|b| * Ja,b
! Ja,b = Ra,b * E(a) * E(b) / Sa+b
!
! The end summation is then divided by factor: sqrt(cg,p*cg,x*cg,y)
!
! All the velocity terms cancel to simply a factor cg,p, i.e:
! (cg,p*Rx,y/Rx+y) * [Ja,b - Ja+b,-a - Ja+b,-b]
!
! Note(2): the above expression is factored by a missing Df term
! from transforming from discrete to continuous
!
! 4. Argument variables
! AC2 action density
! CGO group velocity
! DEP2 water depth
! IDCMIN minimum counter in directional space
! IDCMAX maximum counter in directional space
! IDDLOW minimum direction that is propagated within a sweep
! IDDTOP maximum direction that is propagated within a sweep
! IMATDA main diagonal of the linear system
! IMATRA right-hand side of system of equations
! ISSTOP maximum frequency counter in a sweep
! KWAVE wave number
! PLTRI triad contribution in TEST points
! SPCSIG relative frequencies in computational domain in sigma-space
!
INTEGER IDDLOW, IDDTOP, ISSTOP
INTEGER IDCMIN(MSC), IDCMAX(MSC)
!
REAL :: AC2(MDC,MSC,MCGRD)
!
REAL :: CGO(MSC,MICMAX)
REAL :: DEP2(MCGRD)
REAL :: IMATDA(MDC,MSC), IMATRA(MDC,MSC)
REAL :: SPCSIG(MSC)
REAL :: KWAVE(MSC,MICMAX)
REAL :: PLTRI(MDC,MSC,NPTST)
REAL :: REDC0(MDC,MSC,MREDS)
REAL :: REDC1(MDC,MSC,MREDS)
!
! Values from common
! MDC : Size of array in theta-direction
! MSC : Size of array in sigma-direction
! PI : Circular constant Pi
! PTRIAD : Tunable coefficients for nonlinear triad sourceterms
!
! PTRIAD(6) : value for proportionality coefficient K
! PTRIAD(7) : value for proportionality coefficient K
!
! 6. Local variables
!
INTEGER IS, ID, II
INTEGER IP ,IM, IT
!
REAL :: KP, K
!
REAL :: K1, K2, K12, W1, W2, W12, C1, C2, C12
REAL :: V0, R0, S0 , MU, DK
REAL :: R1, R2, R3 , S1, S2, S3
REAL :: I1, I2, I3 , EE1, EE2, EE3, CC1, CC2, CC3
!
REAL :: SIGPI, STRI
!
! --- for use of COAPP
REAL, ALLOCATABLE :: SX(:,:,:)
REAL, ALLOCATABLE :: KX(:,:,:)
REAL, ALLOCATABLE :: CGX(:,:,:)
REAL, ALLOCATABLE :: CX(:,:,:)
REAL, ALLOCATABLE :: E2(:,:)
REAL, ALLOCATABLE :: EI(:,:)
REAL, ALLOCATABLE :: EX(:,:,:,:)
REAL, ALLOCATABLE :: EIX(:,:,:,:)
REAL, ALLOCATABLE :: WX(:,:,:)
!
REAL, ALLOCATABLE :: CONTS(:,:) !sum contribution
REAL, ALLOCATABLE :: CONTD(:,:) !difference contribution
!
! 7. SUBROUTINES USED
!
! interaction coefficients are calculated in subroutine TCOEF
! interpolation and integration of energy densities in subroutine COAPP
! 8. SUBROUTINES CALLING
! SOURCE
! 9. ERROR MESSAGES
! ---
! 10. REMARKS
! ---
! 11. STRUCTURE
! -----------------------------------------------------------------
! -----------------------------------------------------------------
! 13. Source text
!
! --- allocate and initialize
!
ALLOCATE (SX(MSC,MSC,2))
ALLOCATE (KX(MSC,MSC,2))
ALLOCATE (CGX(MSC,MSC,2))
ALLOCATE (CX(MSC,MSC,2))
ALLOCATE (E2(MDC,MSC))
ALLOCATE (EI(MDC,MSC))
ALLOCATE (EX(MDC,MSC,MSC,2))
ALLOCATE (EIX(MDC,MSC,MSC,2))
ALLOCATE (WX(MSC,MSC,2))
!
ALLOCATE (CONTS(1:MDC,1:MSC))
ALLOCATE (CONTD(1:MDC,1:MSC))
!
CONTS = 0.
CONTD = 0.
!
! --- call COAPP for interpolation and integration of energy densities
CALL COAPP(AC2,SPCSIG,DEP2,IDCMIN,IDCMAX,1,SX,
& KX,CX,CGX,E2,EI,EX,EIX,WX)
!
! --- loop over directions propagated within sweep
DO II = IDDLOW, IDDTOP
ID = MOD ( II - 1 + MDC , MDC ) + 1
DO IP = 1, MSC
! --- define properties at IP
W1 = SPCSIG(IP)
K1 = KWAVE(IP,1)
C1 = CGO(IP,1)
!
! --- sum contribution (m, p - m)
!
IF (IP.GT.1) THEN
DO IM = 1, IP-1 !Eq. m
! --- define remaining properties in sum triad
W2 = SPCSIG(IM)
K2 = KWAVE(IM,1)
C2 = CGO(IM,1)
W12 = SX(IP,IM,1) !Eq. W1 - W2
K12 = KX(IP,IM,1)
C12 = CGX(IP,IM,1)
KP = MIN(K1,K2,K12)
K = (PTRIAD(6)*KP) + PTRIAD(7)
! --- SPB vars.
! DK = (p - m) + m - p
DK = K12 + K2 - K1
MU = K / ((DK**2) + (K**2))
! --- coupling coefficients
CALL TCOEF(W1 ,W2, W12,K1 ,K2, K12,DEP2,R0,S0,2)
CALL TCOEF(W12,W1,-W2 ,K12,K1,-K2 ,DEP2,R2,S2,2)
CALL TCOEF(W2 ,W1,-W12,K2 ,K1,-K12,DEP2,R3,S3,2)
CC1 = R0 / S0
CC2 = R2 / S2
CC3 = R3 / S3
!
! --- compute energy products
! (improved collinear approximation)
!
EE1 = 0.5 * ( E2(ID,IM) * EIX(ID,IP,IM,1)
& + EI(ID,IM) * EX(ID,IP,IM,1))
EE2 = 0.5 * ( E2(ID,IP) * EI(ID,IM)
& + EI(ID,IP) * E2(ID,IM))
EE3 = 0.5 * ( E2(ID,IP) * EIX(ID,IP,IM,1)
& + EI(ID,IP) * EX(ID,IP,IM,1))
!
! --- compute contribution
!
CONTS(ID,IP) = CONTS(ID,IP) + 0.5 * C1 * MU *
& (R0/S0) * (CC1*EE1 - CC2*EE2 - CC3*EE3) *
& PTRIAD(1) * ((FRINTF*W2)/(2.*PI))
ENDDO
ENDIF
!
! --- difference contribution (p + m, m)
!
DO IM = 1, MSC !Eq. m
! --- define rem. properties in diff. triad
W2 = SPCSIG(IM)
K2 = KWAVE(IM,1)
C2 = CGO(IM,1)
W12 = SX(IP,IM,2) !Eq. W1 + W2
K12 = KX(IP,IM,2)
C12 = CGX(IP,IM,2)
KP = MIN(K1,K2,K12)
K = (PTRIAD(6)*KP) + PTRIAD(7)
! --- SPB vars.
! DK = (p + m) - m - p
DK = K12 - K2 - K1
MU = K / ((DK**2) + (K**2))
! --- coupling coefficients
CALL TCOEF(W12,W1,W2,K12,K1,K2,DEP2,R1,S1,2)
CALL TCOEF(W2,-W1,W12,K2,-K1,K12,DEP2,R2,S2,2)
CALL TCOEF(W1,-W2,W12,K1,-K2,K12,DEP2,R0,S0,2)
CC1 = R1 / S1
CC2 = R2 / S2
CC3 = R0 / S0
!
! --- compute energy products
! (improved collinear approximation)
!
EE1 = 0.5 * ( E2(ID,IP) * EI(ID,IM)
& + EI(ID,IP) * E2(ID,IM))
EE2 = 0.5 * ( E2(ID,IP) * EIX(ID,IP,IM,2)
& + EI(ID,IP) * EX(ID,IP,IM,2))
EE3 = 0.5 * ( E2(ID,IM) * EIX(ID,IP,IM,2)
& + EI(ID,IM) * EX(ID,IP,IM,2))
!
! --- compute contributions
!
CONTD(ID,IP) = CONTD(ID,IP) + 0.5 * C1 * MU *
& (R0/S0) * (CC1*EE1 - CC2*EE2 - CC3*EE3) *
& PTRIAD(1) * ((FRINTF*W2)/(2.*PI))
ENDDO
ENDDO
ENDDO
!
! --- compute total contribution
!
DO IS = 1, ISSTOP
SIGPI = SPCSIG(IS) * 2. * PI ! Factor to convert back to N(sigma)
DO II = IDCMIN(IS), IDCMAX(IS)
ID = MOD ( II - 1 + MDC , MDC ) + 1
!
STRI = (4. * CONTS(ID,IS)) - (8. * CONTD(ID,IS))
!
IF(TESTFL) PLTRI(ID,IS,IPTST) = STRI / SIGPI
!
! --- store results in RHS and main diagonal according
! to Patankar-rules
!
IF (STRI .GT. 0.) THEN
IMATRA(ID,IS) = IMATRA(ID,IS) + STRI / SIGPI
REDC0(ID,IS,2)= REDC0(ID,IS,2)+ STRI / SIGPI
ELSE
IMATDA(ID,IS) = IMATDA(ID,IS) - STRI /
& MAX(1.E-18,AC2(ID,IS,KCGRD(1))*SIGPI)
REDC1(ID,IS,2)= REDC1(ID,IS,2)+ STRI /
& MAX(1.E-18,AC2(ID,IS,KCGRD(1))*SIGPI)
ENDIF
ENDDO
ENDDO
!
DEALLOCATE(SX,KX,CGX,CX,E2,EI,EX,EIX,WX,CONTS,CONTD)
!
RETURN
END subroutine SWSPB
!
!********************************************************************
!
SUBROUTINE STRICL (ACLOC ,DEPLOC ,SPCSIG ,KWAVE ,
& IDDLOW ,IDDTOP ,ANYBIN ,IMATDA ,IMATRA,
& CGO ,KMESPC ,ETOT ,SMEBRK ) 40.45
!
!********************************************************************
!
USE OCPCOMM4
USE SWCOMM3
USE SWCOMM4
IMPLICIT NONE
!
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 Delft University of Technology
!
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
! 40.45: Nico Booij
! 40.96: Matthijs Benit
! 41.45: James Salmon
! 1. Updates
! 40.45, July 04: new subroutine
! 41.45, October 2013: energy conservative form
! 2. Purpose
! In this subroutine the triad-wave interactions are calculated
! with the empiric distributed colinear approximation.
! 3. Method
! Transfer of action between two components under influence of a
! third one is formulated as follows:
! P+1 P+1
! M * N (Sigma N - Sigma N )
! 3 1 1 2 2
! in which:
! M is a dimensional coefficient
! N is action density
! Sigma spectral (angular) frequency
! Dimensions:
! ACLOC#: m2s2
! KW# : 1/m
! SIG# : 1/s
! SDIA# : 1/s
! 4. Argument variables
REAL, INTENT(IN) :: ACLOC(1:MDC,1:MSC) ! local action density spectrum
REAL, INTENT(IN) :: DEPLOC ! Depth at gridpoint ix,iy (obtained from SWANCOM1)
REAL, INTENT(IN) :: SPCSIG(1:MSC) ! Relative frequencies in computational domain in sigma-space
REAL :: KWAVE(1:MSC) ! Wave number in stencil points
INTEGER, INTENT(IN) :: IDDLOW ! Minimum counter in directional space
INTEGER, INTENT(IN) :: IDDTOP ! Maximum counter in directional space
LOGICAL, INTENT(IN) :: ANYBIN(1:MDC,1:MSC) ! if True this bin is going to be updated using the matrix
REAL :: IMATDA(1:MDC,1:MSC) ! IMATDA: Diagonal of matrix
REAL :: IMATRA(1:MDC,1:MSC) ! IMATRA: Right hand vector of matrix
REAL, INTENT(IN) :: CGO(1:MSC) ! Group velocities in stencil points
REAL, INTENT(IN) :: KMESPC ! Mean wave number of the spectrum
REAL, INTENT(IN) :: ETOT !
REAL, INTENT(IN) :: SMEBRK !
! Values from common
! MDC : Size of array in theta-direction
! MSC : Size of array in sigma-direction
! PI : Circular constant Pi
! PTRIAD(1) : Interaction coefficient (lambda)
! PTRIAD(2) : Power of the tail of the spectrum (p)
! PTRIAD(4) : xxx (delta)
! 6. Local variables
INTEGER, SAVE :: IENT=0 ! Number of entries into this subroutine
INTEGER :: ID ! Grid counter in spectral space (direction)
INTEGER :: II ! Counter
INTEGER :: IS1, IS2, IS3 ! Grid counter for spectral frequency
REAL :: SIG1, SIG2, SIG3 ! frequencies of 3 components
REAL :: E1, E2, E3 ! energy densities of 2 components
REAL :: KW1, KW2, KW3 ! wave numbers of 3 components
REAL :: CG1, CG2 ! wave group velocities at interacting frequencies
REAL :: RS3, SS ! aux. var. for determining SIG3
REAL :: DSIG ! frequency increment
REAL :: SIGMEAN ! mean freq of the 3 components
REAL :: KMEAN ! mean wave number
REAL :: SDIA1, SDIA2 ! source term to diagonal
REAL :: SRHS
REAL :: DISPC ! dispersion coefficient
REAL :: BETA
REAL :: SINABS
REAL :: BIPH
REAL :: URSLOC ! auxiliary coefficients
!
REAL :: ED1, ED2, ED3 ! integration of energy density over all directions
REAL :: SDIA1B, SDIA2B
REAL :: F(2)
INTEGER :: ISHAP
REAL :: SMEAN
REAL :: SHAP(2) ! coefficient for shape
!
REAL, ALLOCATABLE :: ED(:) ! integration of energy density over all directions
! 7. SUBROUTINES USED
! 8. SUBROUTINES CALLING
! SOURCE
! 9. ERROR MESSAGES
! ---
! 10. REMARKS
! ---
! 11. STRUCTURE
! -----------------------------------------------------------------
! For all active directions do
! For all first components do
! determine frequency and wave number of component
! For all second components do
! determine frequency and wave number of components
! determine frequency of third resonating component
! If this frequency is within spectral range
! Then determine source terms
! determine contributions to matrix
! -----------------------------------------------------------------
! 13. Source text
IF (LTRACE) CALL STRACE (IENT,'STRICL')
!
ISHAP = 0
!
URSLOC = (GRAV*2.*SQRT(ETOT))/(SQRT(2.)*SMEBRK**2*DEPLOC**2)
BIPH = 0.5 * PI * (TANH(PTRIAD(4)/URSLOC) - 1.)
SINABS = ABS(SIN(BIPH))
BETA = PTRIAD(1) / DEPLOC / DEPLOC * SINABS *
& KMESPC**(1.-PTRIAD(2))
ALLOCATE(ED(MSC))
ED = 0.
!
! --- calculate integral of energy density over directions
!
ED(:) = SUM(ACLOC(:,:),DIM=1) * DDIR
! first loop over all directions
DO II = IDDLOW, IDDTOP
ID = MOD (II - 1 + MDC, MDC) + 1
! first loop over all frequencies
DO IS1 = 1, MSC
SIG1 = SPCSIG(IS1)
! second loop over all higher frequencies
DO IS2 = IS1+1, MSC
SIG2 = SPCSIG(IS2)
! determine properties of 3rd component
SIG3 = SIG2 - SIG1
IF (SIG3.GT.SPCSIG(1)) THEN
SS = ALOG(SIG3/SPCSIG(1)) / FRINTF
IS3 = INT(SS) + 1
RS3 = SS - REAL(IS3 - 1)
!
! --- convert action densities to energy densities
!
E1 = ACLOC(ID,IS1) * SIG1
E2 = ACLOC(ID,IS2) * SIG2
E3 = RS3 * ACLOC(ID,IS3+1) * SPCSIG(IS3+1)
& + (1. - RS3) * ACLOC(ID,IS3 ) * SPCSIG(IS3 )
!
KW1 = KWAVE(IS1)
KW2 = KWAVE(IS2)
KW3 = RS3 * KWAVE(IS3+1) + (1. - RS3) * KWAVE(IS3)
CG1 = CGO(IS1)
CG2 = CGO(IS2)
! determine properties of triad
KMEAN = (KW1 + KW2 + KW3)/3.
DISPC = TANH(2.*KMEAN * DEPLOC) / (2.*KMEAN * DEPLOC)
F(1) = 1. !(SIG1 ) / SMEAN
F(2) = 1. !(SIG2 ) / SMEAN
SDIA1 = BETA * E3 * CG1 * KW1**(PTRIAD(2)) * DISPC
SDIA2 = BETA * E3 * CG2 * KW2**(PTRIAD(2)) * DISPC
!
! --- integrals over all directions
!
ED1 = ED(IS1) * SIG1
ED2 = ED(IS2) * SIG2
ED3 = RS3 * ED(IS3+1) * SPCSIG(IS3+1)
& + (1. - RS3) * ED(IS3 ) * SPCSIG(IS3 )
SDIA1B= BETA * ED3 * SIG1 * CG1 * KW1**(PTRIAD(2)) * DISPC
SDIA2B= BETA * ED3 * SIG2 * CG2 * KW2**(PTRIAD(2)) * DISPC
!
! --- apply shape function to relaxation
!
IF (ISHAP.EQ.1) THEN
SHAP(1) = EXP(-1.25*((SIG1/(2.*PI))/0.4)**-4)
SHAP(2) = EXP(-1.25*((SIG2/(2.*PI))/0.4)**-4)
!
SDIA1 = SDIA1 / SHAP(1)
SDIA2 = SDIA2 / SHAP(2)
SDIA1B = SDIA1B / SHAP(1)
SDIA2B = SDIA2B / SHAP(2)
ENDIF
!
IF (ANYBIN(ID,IS1)) THEN
DSIG = (FRINTF * SIG2) / SIG1
DSIG = DSIG * F(1)
SRHS = 0.5 * (SDIA2B*E2 + SDIA2*ED2 - SDIA1*ED1)
IMATDA(ID,IS1) = IMATDA(ID,IS1) + SDIA1B * DSIG * 0.5
IMATRA(ID,IS1) = IMATRA(ID,IS1) + SRHS * DSIG
ENDIF
IF (ANYBIN(ID,IS2)) THEN
DSIG = (FRINTF * SIG1) / SIG2
DSIG = DSIG * F(2)
SRHS = 0.5 * (SDIA1B*E1 + SDIA1*ED1 - SDIA2*ED2)
IMATDA(ID,IS2) = IMATDA(ID,IS2) + SDIA2B * DSIG * 0.5
IMATRA(ID,IS2) = IMATRA(ID,IS2) + SRHS * DSIG
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
DEALLOCATE(ED)
RETURN
END subroutine STRICL