! ! 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) ! ! RIAM_SLW (calculate nonlinear four wave interactions by means ! of the exact FD-RIAM technique) ! ! SWINTFXNL (interface with SWAN model to compute nonlinear transfer ! with the XNL method for given action density spectrum) ! ! *** nonlinear 3 wave-wave interactions *** ! ! SWLTA (triad-wave interactions calculated with the Lumped Triad ! Approximation of Eldeberky, 1996) ! !---------------------------------------------------------------------- ! !****************************************************************** ! 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) 2004-2005 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 ! !OMP INCLUDE 'tpcomm.inc' 40.41 ! ! 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: R.C. Ris, N. Booij, | ! | IJ.G. Haagsma, A.T.M.M. Kieftenburg, | ! | M. Zijlema, E.E. Kriezi, | ! | R. Padilla-Hernandez, L.H. Holthuijsen | ! --|-----------------------------------------------------------|-- ! ! ! SWAN (Simulating WAves Nearshore); a third generation wave model ! Copyright (C) 2004-2005 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 ! !OMP INCLUDE 'tpcomm.inc' 40.41 ! ! 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 ) 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) 2004-2005 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 ! ! 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 ! ! 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 ! !OMP INCLUDE 'tpcomm.inc' 40.41 ! ! 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(*) ! 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) = -1. * DSNL(I,ID) / PI3 40.00 END IF ! 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 , & IMATDA ,IMATRA , 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) 2004-2005 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 ! ! 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 ! ! 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 ! !OMP INCLUDE 'tpcomm.inc' 40.41 ! ! 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(*) ! 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 ELSE IMATDA(ID,I) = IMATDA(ID,I) - SFNL(I,ID) / & MAX(1.E-18,AC2(ID,I,KCGRD(1))*SIGPI) 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) 2004-2005 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 ! !OMP INCLUDE 'tpcomm.inc' 40.41 ! ! 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) 2004-2005 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 ! !OMP INCLUDE 'tpcomm.inc' 40.41 ! ! 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) 2004-2005 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 ! !OMP INCLUDE 'tpcomm.inc' 40.41 ! ! 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 ) 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: R.C. Ris, N. Booij, | ! | IJ.G. Haagsma, A.T.M.M. Kieftenburg, | ! | M. Zijlema, E.E. Kriezi, | ! | R. Padilla-Hernandez, L.H. Holthuijsen | ! --|-----------------------------------------------------------|-- ! ! ! SWAN (Simulating WAves Nearshore); a third generation wave model ! Copyright (C) 2004-2005 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 ! ! 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 ! ! 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 ! !OMP INCLUDE 'tpcomm.inc' 40.41 ! ! 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) ! 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)) ELSE IMATDA(ID,IS) = IMATDA(ID,IS) - MEMNL4(ID,IS,KCGRD(1)) / & MAX(1.E-18,AC2(ID,IS,KCGRD(1))) 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 ! !******************************************************************** ! ! << Numerical Computations of the Nonlinear Energy Transfer ! of Gravity Wave Spectra in Finite Water Depth >> ! ! << developed by Noriaki Hashimoto >> ! ! References: N. Hashimoto et al. (1998) ! Komatsu and Masuda (2000) (in Japanese) ! ! ! SWAN (Simulating WAves Nearshore); a third generation wave model ! Copyright (C) 2004-2005 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 ! ! SUBROUTINE RIAM_SLW(LMAX,N,N2,G,H,DQ,DQ2,DT,DT2,W,P,ACT,SNL,MINT) 40.17 USE M_SNL4 40.41 ! LMAX ! N : number of directional bins 40.17 ! N2 ! G : gravitational acceleration 40.17 ! H : depth 40.17 ! DQ ! DQ2 ! DT : size of the directional bins (Delta Theta) 40.17 ! DT2 ! W : discretised frequency array 40.17 ! P : density of water 40.17 ! ACT : action density 40.17 ! SNL : Quadruplet source term 40.17 ! MINT ! IW4 : counter for the 4th frequency 40.17 ! W4 : frequency of the fourth quadruplet component 40.17 ! AK4 : wavenumber of the fourth quadruplet component 40.17 ! DNA : coefficient in eq. 17 of Hashimoto (98) 40.17 ! = 2 w4 k4 / Cg(k4) 40.17 ! IW3L ! II ! JJ ! DI ! DJ ! CGK4 : group velocity for the fourth quadruplet component 40.17 REAL :: W(LMAX), ACT(LMAX,N), SNL(LMAX,N) 40.17 ! Initialisation of the quadruplet source term 40.17 SNL(:,:) = 0. 40.17 ! ! ================= DO 130 IW4=1,LMAX ! ================= ! W4=W(IW4) ! WAVE converts nondimensional Sqr(w)d/g to nondimensional kd AK4=WAVE(W4**2*H/G)/H ! CGCMP computes group velocity CALL CGCMP(G,AK4,H,CGK4) ! Calculates the coefficient in equation (17) of Hashimoto (1998) DNA=2.*W4*AK4/CGK4 IW3L=MAX0(1,NINT(IW4-ALOG(3.)/DQ)) ! ------------------------------------------------------------ CALL PRESET(LMAX,N,IW3L,IW4,N2,G,H,W4,AK4,W,DQ,DQ2,DT,DT2,P) 40.41 40.17 ! ------------------------------------------------------------ ! ! ============== DO 130 IT4=1,N ! ============== ! ! ================ CURRIAM => FRIAM 40.41 DO ! (W1 - W3 - W4 - W2) 40.41 ! ================ ! K1=CURRIAM%II(1)+IW4 40.41 K2=CURRIAM%II(2)+IW4 40.41 K3=CURRIAM%II(3)+IW4 40.41 ! IF(K1.LT.1 ) GO TO 120 40.17 IF(K2.GT.LMAX) GO TO 120 40.17 ! GG=CURRIAM%SSS*DNA 40.41 ! M1P= CURRIAM%JJ(1)+IT4 40.41 M1N=-CURRIAM%JJ(1)+IT4 40.41 M2P= CURRIAM%JJ(2)+IT4 40.41 M2N=-CURRIAM%JJ(2)+IT4 40.41 M3P= CURRIAM%JJ(3)+IT4 40.41 M3N=-CURRIAM%JJ(3)+IT4 40.41 ! M1P=MOD(M1P,N) M1N=MOD(M1N,N) M2P=MOD(M2P,N) M2N=MOD(M2N,N) M3P=MOD(M3P,N) M3N=MOD(M3N,N) ! IF(M1P.LT.1) M1P=M1P+N IF(M1N.LT.1) M1N=M1N+N IF(M2P.LT.1) M2P=M2P+N IF(M2N.LT.1) M2N=M2N+N IF(M3P.LT.1) M3P=M3P+N IF(M3N.LT.1) M3N=M3N+N ! A4=ACT(IW4,IT4) ! IF(MINT.EQ.1) THEN ! K01=0 K02=0 K03=0 M01=0 M02=0 M03=0 IF(CURRIAM%DI(1).LE.1.) K01= 1 40.41 IF(CURRIAM%DI(2).LE.1.) K02= 1 40.41 IF(CURRIAM%DI(3).LE.1.) K03= 1 40.41 IF(CURRIAM%DJ(1).LE.1.) M01=-1 40.41 IF(CURRIAM%DJ(2).LE.1.) M02=-1 40.41 IF(CURRIAM%DJ(3).LE.1.) M03=-1 40.41 ! K11=K1+K01 K21=K2+K02 K31=K3+K03 IF(K11.GT.LMAX) K11=LMAX IF(K21.GT.LMAX) K21=LMAX IF(K31.GT.LMAX) K31=LMAX ! K111=K1+K01-1 K211=K2+K02-1 K311=K3+K03-1 IF(K111.LT.1) K111=1 IF(K211.LT.1) K211=1 IF(K311.LT.1) K311=1 ! M1P1=M1P+M01 M2P1=M2P+M02 M3P1=M3P+M03 IF(M1P1.LT.1) M1P1=N IF(M2P1.LT.1) M2P1=N IF(M3P1.LT.1) M3P1=N ! M1N1=M1N-M01 M2N1=M2N-M02 M3N1=M3N-M03 IF(M1N1.GT.N) M1N1=1 IF(M2N1.GT.N) M2N1=1 IF(M3N1.GT.N) M3N1=1 ! M1P11=M1P+M01+1 M2P11=M2P+M02+1 M3P11=M3P+M03+1 IF(M1P11.GT.N) M1P11=1 IF(M2P11.GT.N) M2P11=1 IF(M3P11.GT.N) M3P11=1 ! M1N11=M1N-M01-1 M2N11=M2N-M02-1 M3N11=M3N-M03-1 IF(M1N11.LT.1) M1N11=N IF(M2N11.LT.1) M2N11=N IF(M3N11.LT.1) M3N11=N ! DI1=CURRIAM%DI(1) 40.41 DI2=CURRIAM%DI(2) 40.41 DI3=CURRIAM%DI(3) 40.41 DJ1=CURRIAM%DJ(1) 40.41 DJ2=CURRIAM%DJ(2) 40.41 DJ3=CURRIAM%DJ(3) 40.41 ! A1P=(DI1*(DJ1*ACT(K11, M1P1)+ACT(K11, M1P11)) 40.41 & + DJ1*ACT(K111,M1P1)+ACT(K111,M1P11)) 40.41 & /((1.+DI1)*(1.+DJ1)) 40.41 ! A1N=(DI1*(DJ1*ACT(K11, M1N1)+ACT(K11, M1N11)) 40.41 & + DJ1*ACT(K111,M1N1)+ACT(K111,M1N11)) 40.41 & /((1.+DI1)*(1.+DJ1)) 40.41 ! A2P=(DI2*(DJ2*ACT(K21, M2P1)+ACT(K21, M2P11)) 40.41 & + DJ2*ACT(K211,M2P1)+ACT(K211,M2P11)) 40.41 & /((1.+DI2)*(1.+DJ2)) 40.41 ! A2N=(DI2*(DJ2*ACT(K21, M2N1)+ACT(K21, M2N11)) 40.41 & + DJ2*ACT(K211,M2N1)+ACT(K211,M2N11)) 40.41 & /((1.+DI2)*(1.+DJ2)) 40.41 ! A3P=(DI3*(DJ3*ACT(K31, M3P1)+ACT(K31, M3P11)) 40.41 & + DJ3*ACT(K311,M3P1)+ACT(K311,M3P11)) 40.41 & /((1.+DI3)*(1.+DJ3)) 40.41 ! A3N=(DI3*(DJ3*ACT(K31, M3N1)+ACT(K31, M3N11)) 40.41 & + DJ3*ACT(K311,M3N1)+ACT(K311,M3N11)) 40.41 & /((1.+DI3)*(1.+DJ3)) 40.41 ! ELSE ! IF(K1.LT.1) K1=1 IF(K2.LT.1) K2=1 IF(K3.LT.1) K3=1 ! A1P=ACT(K1,M1P) A1N=ACT(K1,M1N) A2P=ACT(K2,M2P) A2N=ACT(K2,M2N) A3P=ACT(K3,M3P) A3N=ACT(K3,M3N) ! ENDIF ! W1P2P=A1P+A2P W1N2N=A1N+A2N S1P2P=A1P*A2P S1N2N=A1N*A2N W3P4 =A3P+A4 W3N4 =A3N+A4 S3P4 =A3P*A4 S3N4 =A3N*A4 ! XP=S1P2P*W3P4-S3P4*W1P2P XN=S1N2N*W3N4-S3N4*W1N2N ! SNL( K1,M1P)=SNL( K1,M1P)-XP*GG SNL( K1,M1N)=SNL( K1,M1N)-XN*GG SNL( K2,M2P)=SNL( K2,M2P)-XP*GG SNL( K2,M2N)=SNL( K2,M2N)-XN*GG SNL( K3,M3P)=SNL( K3,M3P)+XP*GG SNL( K3,M3N)=SNL( K3,M3N)+XN*GG SNL(IW4,IT4)=SNL(IW4,IT4)+(XP+XN)*GG ! ! ============ 120 CONTINUE IF (.NOT.ASSOCIATED(CURRIAM%NEXTRIAM)) EXIT 40.41 CURRIAM => CURRIAM%NEXTRIAM 40.41 END DO 40.41 ! ============ ! ! ============ 130 CONTINUE ! ============ ! RETURN END SUBROUTINE PRESET(LMAX,N,IW3L,IW4,N2,G,H,W4,AK4,W,DQ,DQ2,DT,DT2,P) 40.41 40.17 ! T1A : Angle between theta_1 and theta_a ! T2A : Angle between theta_2 and theta_a ! T34 : Angle between theta_3 and theta_4 ! WA : wa = w1 + w2 = w3 + w4 (resonance conditions) ! AKA : absolute value of ka = k1 + k2 = k3 + k4 (resonance conditions) ! TA : Angle theta_a representing vector ka REAL :: W(LMAX) 40.17 ! ! ================ DO 120 IT34=1,N2 ! ================ ! T34=DT*(IT34-1) DT3=DT IF(IT34.EQ.1.OR.IT34.EQ.N2) DT3=DT2 ! ! =================== DO 110 IW3=IW3L,IW4 ! =================== ! W3=W(IW3) AK3=WAVE(W3**2*H/G)/H ! WA=W3+W4 AKA=SQRT(AK3*AK3+AK4*AK4+2.*AK3*AK4*COS(T34)) TA=ATAN2(AK3*SIN(T34),AK3*COS(T34)+AK4) R=SQRT(G*AKA*TANH(AKA*H/2.))/WA-1./SQRT(2.) ! TL=0. ACS=AKA/(2.*WAVE(WA**2*H/(4.*G))/H) ACS=SIGN(1.,ACS)*MIN(1.,ABS(ACS)) 40.41 IF(R.LT.0.) TL=ACOS(ACS) IT1S=NINT(TL/DT)+1 ! ! =================== DO 100 IT1A=IT1S,N2 ! =================== ! T1A=DT*(IT1A-1) ! DT1=DT IF(IT1A.EQ.1.OR.IT1A.EQ.N2) DT1=DT2 ! ! ----------------------------------------------------------------- CALL FINDW1W2(LMAX,IW3,W,G,H,W1,W2,W3,WA,AK1,AK2,AKA,T1A,T2A,IND) 40.17 ! ----------------------------------------------------------------- IF(IND.LT.0) GO TO 100 ! IF(W1.LE.W3.AND.W3.LE.W4.AND.W4.LE.W2) THEN CALL KERNEL(N,G,H,W1,W2,W3,W4,AK1,AK2,AK3,AK4,AKA, & T1A,T2A,T34,TA,DQ,DT,DT1,DT3,P) 40.41 40.17 ENDIF ! ! ============ 100 CONTINUE 110 CONTINUE 120 CONTINUE ! ============ ! RETURN END SUBROUTINE FINDW1W2(LMAX, 40.17 & IW3,W,G,H,W1,W2,W3,WA,AK1,AK2,AKA,T1A,T2A,IND) 40.17 REAL :: W(LMAX) 40.17 ! IND=0 EPS=0.0005 ! X1=0.005 X2=W3 ! ! --------------------------------------- 110 CALL FDW1(G,H,AKA,T1A,WA,X1,X2,X,EPS,M) ! --------------------------------------- ! IF(M.EQ.0) GO TO 999 ! W1=X ! AK1=WAVE(W1**2*H/G)/H AK2=SQRT(AKA*AKA+AK1*AK1-2.*AKA*AK1*COS(T1A)) W2=SQRT(G*AK2*TANH(AK2*H)) IF(W1.GT.W2) GO TO 999 T2A=ATAN2(-AK1*SIN(T1A),AKA-AK1*COS(T1A)) RETURN ! 999 IND=-999 RETURN END SUBROUTINE FDW1(G,H,AKA,T1A,WA,X1,X2,X,EPS,M) ! M=0 F1=FUNCW1(G,H,X1,AKA,T1A,WA) F2=FUNCW1(G,H,X2,AKA,T1A,WA) 20 IF(F1*F2) 18,19,21 18 M=M+1 X=X2-((X2-X1)/(F2-F1)*F2) IF((ABS(X-X2)/ABS(X))-EPS) 22,23,23 23 IF((ABS(X-X1)/ABS(X))-EPS) 22,24,24 22 RETURN 24 F=FUNCW1(G,H,X,AKA,T1A,WA) FM=F*F1 IF(FM) 25,22,26 25 X2=X F2=F GO TO 20 26 X1=X F1=F GO TO 20 19 M=M+1 IF(F1) 17,16,17 16 X=X1 GO TO 22 17 X=X2 GO TO 22 21 M=0 RETURN END FUNCTION FUNCW1(G,H,X,AKA,T1A,WA) ! AK1=WAVE(X**2*H/G)/H AK2=SQRT(AKA*AKA+AK1*AK1-2.*AKA*AK1*COS(T1A)) FUNCW1=WA-SQRT(G*AK1*TANH(AK1*H))-SQRT(G*AK2*TANH(AK2*H)) ! RETURN END SUBROUTINE KERNEL(N,G,H,W1,W2,W3,W4,AK1,AK2,AK3,AK4,AKA, & T1A,T2A,T34,TA,DQ,DT,DT1,DT3,P) 40.41 40.17 USE M_SNL4 40.41 ! LOGICAL, SAVE :: LRIAM = .FALSE. 40.41 TYPE(RIAMDAT), POINTER :: RIAMTMP 40.41 ! PI =4.*ATAN(1.) PI2=2.*PI ! ! ------------------------ CALL CGCMP(G,AK1,H,CGK1) CALL CGCMP(G,AK2,H,CGK2) CALL CGCMP(G,AK3,H,CGK3) ! ------------------------ ! Calculate denominator S (eq. 14, Hashimoto (1998)) SS0=ABS(1.+CGK2/CGK1*(AK1-AKA*COS(T1A))/AK2) IF(SS0.LE.1.E-7) RETURN ! k1 k3 w3 G ! Kernel function ------------- - dth3 dth1 dOm (eq. 17) ! Cg(k1) Cg(k3) S ! where S defined in eq. 14 and G defined in eq. 5 (Hashimoto, 1998) CF=AK1*AK3*W3/(CGK1*CGK3)/SS0 & *9.*PI*G*G/(4.*P*P*W1*W2*W3*W4)*DT1*DT3*DQ ! -------------------------------------- CALL NONKDD(G,H,W3,W4,-W1,AK3,AK4,AK1, & T34,0.,T1A+TA+PI,DDD1) ! -------------------------------------- GGG1=CF*DDD1*DDD1 ! ! -------------------------------------- CALL NONKDD(G,H,W3,W4,-W1,AK3,AK4,AK1, & T34,0.,-T1A+TA+PI,DDD2) ! -------------------------------------- GGG2=CF*DDD2*DDD2 ! I1=1-INTW(W1/W4,DQ) I2=1-INTW(W2/W4,DQ) I3=1-INTW(W3/W4,DQ) ! DI1=WEIGW(W1/W4,DQ) DI2=WEIGW(W2/W4,DQ) DI3=0. ! ! ---------------------------- ! TT1=T1A+TA TT2=T2A+TA TT3=T34 ! ALLOCATE(RIAMTMP) 40.41 ! RIAMTMP%SSS =GGG1 40.41 RIAMTMP%II(1)=I1 40.41 RIAMTMP%II(2)=I2 40.41 RIAMTMP%II(3)=I3 40.41 RIAMTMP%DI(1)=DI1 40.41 RIAMTMP%DI(2)=DI2 40.41 RIAMTMP%DI(3)=DI3 40.41 RIAMTMP%JJ(1)=INTT(TT1,DT,PI2,N)-1 40.41 RIAMTMP%JJ(2)=INTT(TT2,DT,PI2,N)-1 40.41 RIAMTMP%JJ(3)=INTT(TT3,DT,PI2,N)-1 40.41 RIAMTMP%DJ(1)=WEIGT(TT1,DT,PI2,N) 40.41 RIAMTMP%DJ(2)=WEIGT(TT2,DT,PI2,N) 40.41 RIAMTMP%DJ(3)=0. 40.41 ! NULLIFY(RIAMTMP%NEXTRIAM) 40.41 IF ( .NOT.LRIAM ) THEN 40.41 FRIAM = RIAMTMP 40.41 CURRIAM => FRIAM 40.41 LRIAM = .TRUE. 40.41 ELSE CURRIAM%NEXTRIAM => RIAMTMP 40.41 CURRIAM => RIAMTMP 40.41 END IF ! ! ---------------------------- ! TT1=-T1A+TA TT2=-T2A+TA TT3= T34 ! ALLOCATE(RIAMTMP) 40.41 ! RIAMTMP%SSS =GGG2 40.41 RIAMTMP%II(1)=I1 40.41 RIAMTMP%II(2)=I2 40.41 RIAMTMP%II(3)=I3 40.41 RIAMTMP%DI(1)=DI1 40.41 RIAMTMP%DI(2)=DI2 40.41 RIAMTMP%DI(3)=DI3 40.41 RIAMTMP%JJ(1)=INTT(TT1,DT,PI2,N)-1 40.41 RIAMTMP%JJ(2)=INTT(TT2,DT,PI2,N)-1 40.41 RIAMTMP%JJ(3)=INTT(TT3,DT,PI2,N)-1 40.41 RIAMTMP%DJ(1)=WEIGT(TT1,DT,PI2,N) 40.41 RIAMTMP%DJ(2)=WEIGT(TT2,DT,PI2,N) 40.41 RIAMTMP%DJ(3)=0. 40.41 ! NULLIFY(RIAMTMP%NEXTRIAM) 40.41 IF ( .NOT.LRIAM ) THEN 40.41 FRIAM = RIAMTMP 40.41 CURRIAM => FRIAM 40.41 LRIAM = .TRUE. 40.41 ELSE CURRIAM%NEXTRIAM => RIAMTMP 40.41 CURRIAM => RIAMTMP 40.41 END IF ! RETURN END FUNCTION INTW(W,DQ) AA=ALOG(W)/DQ INTW=1+NINT(-AA) RETURN END FUNCTION WEIGW(W,DQ) AA=ALOG(W)/DQ L=1+INT(-AA) A=DQ*(1-L)-ALOG(W) IF(A.NE.0.) THEN B=ALOG(W)+DQ*L WEIGW=B/A ELSE WEIGW=1000. ENDIF RETURN END FUNCTION INTT(T,DT,PI2,MM) T=AMOD(T,PI2) IF(T.LT.0.) T=T+PI2 INTT=NINT(T/DT)+1 IF(INTT.GT.MM) INTT=INTT-MM IF(INTT.LT.1) INTT=INTT+MM RETURN END FUNCTION WEIGT(T,DT,PI2,MM) T=AMOD(T,PI2) IF(T.LT.0.) T=T+PI2 N=INT(T/DT)+1 IF(N.GT.MM) N=N-MM IF(N.LT.1) N=N+MM C=T-DT*(N-1) IF(C.NE.0.) THEN D=DT*N-T WEIGT=D/C ELSE WEIGT=1000. ENDIF RETURN END SUBROUTINE NONKDD(G,H,W1,W2,W3,AK1,AK2,AK3,T1,T2,T3,DDD) ! See Herterich and Hasselmann (1980; p. 223, eq. B1) CALL NONKD(G,H,W1,W2,W3,AK1,AK2,AK3,T1,T2,T3,DD1) CALL NONKD(G,H,W2,W3,W1,AK2,AK3,AK1,T2,T3,T1,DD2) CALL NONKD(G,H,W3,W1,W2,AK3,AK1,AK2,T3,T1,T2,DD3) ! DDD=(DD1+DD2+DD3)/3. ! RETURN END SUBROUTINE NONKD(G,H,W1,W2,W3,AK1,AK2,AK3,T1,T2,T3,DDD) ! GG : square of acceleration of gravity (=g**2) ! G2 : twice the square of acceleration of gravity (=2g**2) ! WP23 : sum of w2 and w3 (=w2+w3) ! W123 : sum of w1, w2 and w3 (=w1+w2+w3) ! AKX1 : size of the x-component of wavenumber k1 ! AKY1 : size of the y-component of wavenumber k1 ! AKX2 : size of the x-component of wavenumber k2 ! AKY2 : size of the y-component of wavenumber k2 ! AKX3 : size of the x-component of wavenumber k3 ! AKY3 : size of the y-component of wavenumber k3 ! AK23 : dot product of wavenumbers k2 and k3 ! W23 : frequency corresponding to wavenumber k2+k3 DDD=0. GG=G*G G2=2.*G*G ! AKX1=AK1*COS(T1) AKY1=AK1*SIN(T1) ! AKX2=AK2*COS(T2) AKY2=AK2*SIN(T2) ! AKX3=AK3*COS(T3) AKY3=AK3*SIN(T3) ! AK23=SQRT((AKX2+AKX3)**2+(AKY2+AKY3)**2) W23=SQRT(G*AK23*TANH(AK23*H)) ! WP23=W2+W3 CF0=W23**2-WP23**2 IF(CF0.EQ.0.) RETURN ! W123=W1+W2+W3 AKXY23=AKX2*AKX3+AKY2*AKY3 AKXY123=AKX1*(AKX2+AKX3)+AKY1*(AKY2+AKY3) CF1=0. IF(AK1*H.LE.10.) CF1=AK1*AK1/COSH(AK1*H)**2 ! CF2=0. CF3=0. IF(AK3*H.LE.10.) CF2=W2*AK3*AK3/COSH(AK3*H)**2 IF(AK2*H.LE.10.) CF3=W3*AK2*AK2/COSH(AK2*H)**2 ! See Herterich and Hasselmann (1980), first eq. after eq. B2 ! Equation is divided by I, to avoid complex numbers DD=WP23*(AK2*AK3*TANH(AK2*H)*TANH(AK3*H)-AKXY23)-(CF2+CF3)/2. ! See Herterich and Hasselmann (1980), second eq. after eq. B2 EE=(AKXY23-W2*W3*(W2**2+W3**2+W2*W3)/GG)/(2.*G) CF4=0. IF(AK23*H.LE.10.) CF4=W1*AK23**2/COSH(AK23*H)**2 ! Equation B1 of Herterich and Hasselmann (1980), where DD is ! multiplied by I, to compensate for the earlier division DDD=-DD/CF0*(2.*W123*(W1**2*W23**2/GG-AKXY123)-CF4-WP23*CF1) & +DD*W1*(W1*W1+W23*W23)/GG & +EE*(W1**3*WP23/G-G*AKXY123-G*CF1) & +W1*AKXY23*(W123*(W2**2+W3**2)+W2*W3*WP23)/G2 & -W1*W2*W2*AK3*AK3*(W123+W3)/G2 & -W1*W3*W3*AK2*AK2*(W123+W2)/G2 ! RETURN END SUBROUTINE CGCMP(G,AK,H,CG) ! Calculates group velocity Cg based on depth and wavenumber ! Includes a deep water limit for kd > 10. ! G : gravitational acceleration ! AK : wave number ! H : depth ! CG : group velocity ! AKH : depth x Wave number (kd) ! RGK : square root of gk ! SECH2 : the square of the secant Hyperbolic of kd ! = 1 - TANH(kd)**2 ! Calculation of group velocity Cg AKH=AK*H IF(AKH.LE.10.) THEN 40.33 ! Shallow water: ! Cg = (1/2) (1 + (2 kd) / Sinh (2 kd)) (w / k) ! = (1/2) (1 + (kd SECH2) / Tanh (kd)) (w / k) ! = (1/2) (1 + (kd SECH2) / Tanh (kd)) (Sqrt (gk Tanh(kd)) / k) ! = (1/2) (Sqrt (gk) (Sqrt(Tanh (kd)) + kd SECH2 / Sqrt (Tanh (kd)) / k ! = (1/2) (Sqrt (k) g (Tanh (kd) + kd SECH2) / (k Sqrt (g Tanh (kd))) ! = (g Tanh(kd) + gkd SECH2) / (2 Sqrt(gk Tanh(kd)) SECH2=1.-TANH(AKH)**2 CG=(G*AKH*SECH2+G*TANH(AKH))/(2.*SQRT(G*AK*TANH(AKH))) ELSE 40.17 ! Deep water: ! Cg = w / (2 k) ! = g / (2 Sqrt (gk)) RGK=SQRT(G*AK) CG=G/(2.*RGK) ENDIF 40.17 ! RETURN END SUBROUTINE CGCMP SUBROUTINE DCGCMP(G,AK,H,DCG) ! AKH=AK*H IF(AKH.GT.10.) GO TO 100 SECH2=1-TANH(AKH)**2 SECH3=SECH2*SQRT(SECH2) DCG=G*H*SECH3*(COSH(AKH)-AKH*SINH(AKH))/SQRT(G*AK*TANH(AKH)) & -(G*AKH*SECH2+G*TANH(AKH))**2/(4.*SQRT(G*AK*TANH(AKH))**3) RETURN ! 100 RGK=SQRT(G*AK) DCG=-G*G/(4.*RGK**3) ! RETURN END REAL FUNCTION WAVE(D) ! Transforms nondimensional Sqr(w)d/g into nondimensional kd using the ! dispersion relation and an iterative method IF(D-10.) 2,2,1 1 XX=D GO TO 6 2 IF(D-1.) 3,4,4 3 X=SQRT(D) GO TO 5 4 X=D 5 COTHX=1.0/TANH(X) XX=X-(X-D*COTHX)/(1.+D*(COTHX**2-1.)) E=1.-XX/X X=XX IF(ABS(E)-0.0005) 6,5,5 6 WAVE=XX RETURN 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) 2004-2005 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 SWLTA ( AC2 , DEP2 , CGO , SPCSIG, & KWAVE , IMATRA, IMATDA, & IDDLOW, IDDTOP, ISSTOP, IDCMIN, IDCMAX, & HS , SMEBRK, PLTRI , URSELL ) ! !**************************************************************** ! USE OCPCOMM4 USE SWCOMM3 USE SWCOMM4 ! IMPLICIT NONE ! !OMP INCLUDE 'tpcomm.inc' ! ! ! --|-----------------------------------------------------------|-- ! | Delft University of Technology | ! | Faculty of Civil Engineering | ! | Environmental Fluid Mechanics Section | ! | P.O. Box 5048, 2600 GA Delft, The Netherlands | ! | | ! | Programmers: R.C. Ris, N. Booij, | ! | IJ.G. Haagsma, A.T.M.M. Kieftenburg, | ! | M. Zijlema, E.E. Kriezi, | ! | R. Padilla-Hernandez, L.H. Holthuijsen | ! --|-----------------------------------------------------------|-- ! ! ! SWAN (Simulating WAves Nearshore); a third generation wave model ! Copyright (C) 2004-2005 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.55: Marcel Zijlema ! ! 1. Updates ! ! 40.55, Feb. 06: New subroutine ! ! 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) ! ! 6. Local variables ! ! AUX1 : auxiliary real ! AUX2 : auxiliary real ! BIPH : parameterized bi-phase 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 ! 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, EM, & FT, RINT, SIGPI, SINBPH, STRI, WISM, WISM1, WISP, WISP1, & W0, WM, WN0, WNM, XIS, XISLN REAL, ALLOCATABLE :: E(:), 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 (SA(1:MDC,1:MSC+ISP1)) E = 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) ) ! 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) W0 = SPCSIG(IS) WN0 = KWAVE(IS,1) C0 = W0 / WN0 IF ( IS.GT.-ISM1 ) THEN EM = WISM * E(IS+ISM1) + WISM1 * E(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. 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 SA(ID,IS) = MAX(0., FT * ( EM * EM - 2. * EM * E0 )) 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 ELSE IMATDA(ID,IS) = IMATDA(ID,IS) - STRI / & MAX(1.E-18,AC2(ID,IS,KCGRD(1))*SIGPI) 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,SA) RETURN END SUBROUTINE SWLTA