#include "swancpp.h" ! ! SWAN/COMPU file 3 of 5 ! ! PROGRAM SWANCOM3.FOR ! ! This file SWANCOM3 of the main program SWAN program ! includes the next subroutines (mainly subroutines for ! the source terms for generation of wave energy ) : ! ! WNDPAR (DOLPHIN-B formulations for the SWAN model for a ! first- or a second generation first guess of the spectrum) ! WINDP1 (computation of variables derived from the wind such as ! mean wind velocity, mean wind direction, minimum counter ! for the wind, maximum counter for the wind, wind friction ! velocity and the Pierson Moskowitz frequency ) ! WINDP2 (computation of wind sea energy spectrum necessary for ! the second generation wind growth model) ! WINDP3 (limit the energy spectrum in the case of a first or ! second generation wind growth model) ! SWIND0 (linear input term Cavaleri and Malanotte Rizolli (1981) ! SWIND3 (third generation wind growth model (Snyder et al. 1981; ! Komen et al., 1984) ! SWIND4 (third generation wind growth model (Janssen, 1989,1991) ! SWIND5 (third generation wind growth model according to the ! expression of Yan (1987) (especially derived for a ! frequency range that extend to the high frequencies ! !**************************************************************** ! SUBROUTINE WNDPAR (ISSTOP,IDWMIN,IDWMAX,IDCMIN,IDCMAX, 32.06 & DEP2 ,WIND10,GENC0 ,GENC1 , 40.85 32.06 & THETAW,AC2 ,KWAVE ,IMATRA,IMATDA, 32.06 & SPCSIG,CGO ,ALIMW ,GROWW ,ETOTW , 32.06 & PLWNDS,PLWNDD,SPCDIR,ITER ) 32.06 ! !**************************************************************** ! USE SWCOMM3 40.41 USE SWCOMM4 40.41 USE OCPCOMM4 40.41 ! IMPLICIT NONE 30.82 ! ! ! --|-----------------------------------------------------------|-- ! | Delft University of Technology | ! | Faculty of Civil Engineering | ! | Environmental Fluid Mechanics Section | ! | P.O. Box 5048, 2600 GA Delft, The Netherlands | ! | | ! | Programmers: The SWAN team | ! --|-----------------------------------------------------------|-- ! ! ! SWAN (Simulating WAves Nearshore); a third generation wave model ! Copyright (C) 1993-2017 Delft University of Technology ! ! This program is free software; you can redistribute it and/or ! modify it under the terms of the GNU General Public License as ! published by the Free Software Foundation; either version 2 of ! the License, or (at your option) any later version. ! ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! A copy of the GNU General Public License is available at ! http://www.gnu.org/copyleft/gpl.html#SEC3 ! or by writing to the Free Software Foundation, Inc., ! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ! ! ! 0. Authors ! ! 30.70: Nico Booij ! 30.72: IJsbrand Haagsma ! 30.75: IJsbrand Haagsma (bug fix) ! 30.82: IJsbrand Haagsma ! 32.06: Roeland Ris ! 40.00: Nico Booij (Nonstationary boundary conditions) ! 40.41: Marcel Zijlema ! 40.85: Marcel Zijlema ! ! 1. Updates ! ! Jan. 97: New subroutine (Roeland Ris) ! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN ! 30.70, Feb. 98: argument list of WINDP2 changed ! 30.75, Mar. 98: Set FPM=SIGPKD, due to change in argument list of WINDP2 ! 40.00, July 98: argument list of WINDP2 changed ! 30.82, Oct. 98: Updated description of several variables ! 30.82, Apr. 99: Dimensioning KCGRD corrected ! 32.06, June 99: Reformulated directional spreading for first guess ! 30.82, June 99: Implicit none added; all variables declared ! 40.41, Aug. 04: COS(THETA-THETAW) replaced by sumrule to make it cheaper ! 40.41, Oct. 04: common blocks replaced by modules, include files removed ! 40.85, Aug. 08: store wind input for output purposes ! ! 2. Purpose ! ! Computation of the wind input source term with formulations ! of a first-generation model (constant porportionality coefficient) ! and a second-generation model (proportionality coefficient depends ! on the energy in the wind sea part of the spectrum). The ! expressions are from Holthuijsen and de Boer (1988) and from ! the DOLPHIN-B model (Holthuijsen and Booij). During the ! implementation of the terms modifications to the code have been ! made after personal communications with Holthuijsen and Booij. ! ! 3. Method ! ! The source term of the following nature: ! ! S = A + B E for E < Elim | t - tw | < pi/2 ! ! (Elim-E) ! S = -------- for E > Elim | t - tw | < pi/2 ! TAU ! ! S = 0 for E > Elim | t - tw | > pi/2 ! ! in which the terms A and B are: ! ! [cf10] 2 2 2 4 ! A = ------ pi (1./g) [rhoaw] [cdrag] (U cos(t-tw)) ! 2 pi ! ! and: ! U cos(t-tw) s ! B = 5 [cf20] [rhoaw] f { ----------- - [cf30] } ---- ! Cph 2 pi ! ! The coefficient TAU in the relaxation model is given by: ! 2 ! / 2 pi \ g ! TAU = [cf40] | ----- | ------------ ! \ s / U cos(d-dw) ! ! The limiting spectrum is given by: ! ! -3 ! ALPHA K s -4 2 2 ! Elim = ----------- exp ( -5/4 (----) ) --- cos ( t - tw ) ! 2 Cg spm pi ! ! in which: ! ! ALPHA wind sea and/or depth dependent proportionality ! coefficient which controls the energy scale of the ! limiting spectrum. ! * In the first-generation model ALPHA is a constant ! equal to 0.0081 (fully developed) ! * In the second-generation model ALPHA depends on the ! energy in the wind sea part of the spectrum. ALPHA ! is calculated here by: ! [cf60] ! ALPHA = [cf50] * Edml ! ! spm adapted Pierson-Moskowitz (1964) peak frequency ! ! The total non-dimensional energy in the wind sea part of the ! spectrum is calculated by (see subroutine WINDP2): ! ! 2 ! grav * ETOTW ! Edml = ------------- ! 4 ! wind10 ! ! ! 4. Argument variables ! ! i SPCDIR: (*,1); spectral directions (radians) 30.82 ! (*,2); cosine of spectral directions 30.82 ! (*,3); sine of spectral directions 30.82 ! (*,4); cosine^2 of spectral directions 30.82 ! (*,5); cosine*sine of spectral directions 30.82 ! (*,6); sine^2 of spectral directions 30.82 ! i SPCSIG: Relative frequencies in computational domain in sigma-space 30.72 ! REAL SPCDIR(MDC,6) 30.82 REAL SPCSIG(MSC) 30.72 ! ! 6. Local variables ! ! IENT : Number of entries into this subroutine ! INTEGER IENT ! REAL :: FPM ! Pierson-Moskowitz frequency REAL :: SWIND_EXP, SWIND_IMP ! explicit and implicit part of wind source ! ! INTEGERS: ! --------- ! IDWMIN Minimum counter for spectral wind direction ! IDWMAX Maximum counter for spectral wind direction ! IX Counter of gridpoint in x-direction ! IY Counter of gridpoint in y-direction ! IS Counter of frequency bin ! ISSTOP Countrer for the maximum frequency of all directions ! IDDUM Dummy counter ! ID Counter of directional distribution ! IDWMIN/IDWMAX Minimum / maximum counter in wind sector (180 degrees) ! ! REALS: ! --------- ! ALPM Coefficient for overshoot at deep water ! ALPMD Coefficient for overshoot corrected for shallow ! water using expression of Bretschnieder (1973) ! ALIMW limiting spectrum in terms of action density ! ARG1, ARG2 Exponent ! CDRAG Wind drag coefficient ! DND Nondimensional depth ! DTHETA Difference in rad between wave and wind direction ! EDML Dimensionless energy ! ETOTW Total energy of the wind sea part of the spectrum ! RHOAW Density of air devided by the density of water ! SIGPK Peak frequency in terms of rad /s ! SIGPKD Adapted peak frequency for shallow water ! TAU Variable for the wind growth equation ! THETA Spectral direction ! THETAW Mean direction of the relative wind vector ! TWOPI Two times pi ! WIND10 Velocity of the relative wind vector ! ! one and more dimensional arrays: ! --------------------------------- ! AC2 4D Action density as function of D,S,X,Y and T ! ALIMW 2D Limiting action density spectrum ! DEP2 1D Depth ! CGO 2D Group velocity ! KWAVE 2D Wave number ! LOGSIG 1D Logaritmic distribution of frequency ! IMATRA 2D Coefficients of right hand side of vector ! IMATDA 2D Coefficients of the diagonal ! PLWNDS 3D Values of explicit part of wind input for test point ! PLWNDD 3D Values of implicit part of wind input for test point ! SPCDIR 1D Spectral direction of wave component ! IDCMIN 1D Minimum counter ! IDCMAX 1D Maximum counter in directional space ! GROWW 2D Aux. array to determine whether there are ! wave generation conditions ! ! PWIND(1) = CF10 188.0 ! PWIND(2) = CF20 0.59 ! PWIND(3) = CF30 0.12 ! PWIND(4) = CF40 250.0 ! PWIND(5) = CF50 0.0023 ! PWIND(6) = CF60 -0.2233 ! PWIND(7) = CF70 0. (not used) ! PWIND(8) = CF80 -0.56 (not used) ! PWIND(9) = RHOAW 0.00125 (density air / density water) ! PWIND(10) = EDMLPM 0.0036 (limit energy Pierson Moskowitz) ! PWIND(11) = CDRAG 0.0012 (drag coefficient) ! PWIND(12) = UMIN 1.0 (minimum wind velocity) ! PWIND(13) = PMLM 0.13 ( ) ! ! 7. Common blocks used ! ! ! 5. SUBROUTINES CALLING ! ! SOURCE ! ! 6. SUBROUTINES USED ! ! WINDP2 (compute the total energy in the wind sea part of ! the spectrum). Subroutine WINDP2 is called in ! SWANCOM1 in subroutine SOURCE ! ! 7. ERROR MESSAGES ! ! --- ! ! 8. REMARKS ! ! --- ! ! 9. STRUCTURE ! ! -------------------------------------------------------------- ! Calculate the adapted peak frequency ! -------------------------------------------------------------- ! If first-generation model ! alpha is constant ! else ! Calculate energy in wind sea part of spectrum ETOTW ! Calculate alpha on the basis of ETOTW ! end ! -------------------------------------------------------------- ! Take depth effects into account for alpha ! -------------------------------------------------------------- ! For each frequency and direction ! compute limiting spectrum and determine whether there is ! grow or decay ! end ! -------------------------------------------------------------- ! Do for each frequency and direction ! If wind-wave generation conditions are present ! calculate A + B E ! else if energy is larger than limiting spectrum ! calculate dissipation rate with relaxation model ! endif ! enddo ! -------------------------------------------------------------- ! Store results in matrix (IMATRA or IMATDA) ! -------------------------------------------------------------- ! End of the subroutine WNDPAR ! -------------------------------------------------------------- ! ! 10. SOURCE ! !*********************************************************************** ! INTEGER IS ,ID ,ITER , & IDWMIN,IDWMAX,IDDUM ,ISSTOP ! REAL WIND10,THETA ,THETAW,EDML ,ARG1 ,ARG2 , & ALPM ,ALPMD ,TEMP1 ,TEMP2 ,FACTA ,FACTB , & ADUM ,BDUM ,CINV ,SIGTPI,SIGMA ,TWOPI ,TAUINV, & SIGPK ,SIGPKD,DND ,ETOTW ,ALIM1D, & CTW ,STW ,COSDIF, 40.41 & DIRDIS,AC2CEN,DTHETA ! REAL :: AC2(MDC,MSC,MCGRD) REAL :: ALIMW(MDC,MSC) REAL :: IMATDA(MDC,MSC), IMATRA(MDC,MSC) ! Changed ICMAX to MICMAX, since MICMAX doesn't vary over gridpoint 40.22 REAL :: KWAVE(MSC,MICMAX) 40.22 REAL :: PLWNDS(MDC,MSC,NPTST) 40.00 REAL :: PLWNDD(MDC,MSC,NPTST) REAL :: GENC0(MDC,MSC,MGENR) 40.85 REAL :: GENC1(MDC,MSC,MGENR) 40.85 REAL :: DEP2(MCGRD) ! Changed ICMAX to MICMAX, since MICMAX doesn't vary over gridpoint 40.22 REAL :: CGO(MSC,MICMAX) 40.22 ! INTEGER IDCMIN(MSC) , & IDCMAX(MSC) ! LOGICAL GROWW(MDC,MSC) ! SAVE IENT DATA IENT/0/ IF (LTRACE) CALL STRACE (IENT,'WNDPAR') ! ! *** initialization of arrays *** ! DO IS = 1, MSC DO ID = 1, MDC GROWW(ID,IS) = .FALSE. ALIMW(ID,IS) = 0. ENDDO ENDDO ! ! *** calculate the adapted shallow water peak frequency *** ! *** according to Bretschneider (1973) using the nondimensional *** ! *** depth DND *** ! TWOPI = 2. * PI DND = MIN( 50. , GRAV * DEP2(KCGRD(1)) / WIND10**2 ) SIGPK = TWOPI * 0.13 * GRAV / WIND10 SIGPKD = SIGPK / TANH(0.833*DND**0.375) FPM = SIGPKD 30.75 CTW = COS(THETAW) 40.41 STW = SIN(THETAW) 40.41 ! IF ( IWIND .EQ. 1 ) THEN ! ! *** first generation model *** ! ALPM = 0.0081 ! ELSE IF (IWIND .EQ. 2 ) THEN ! ! *** second generation model *** ! ! *** Determine the proportionality constant alpha on the basis *** ! *** of the total energy in the wind sea part of the spectrum *** ! *** output of subroutine (WINDP2) is ETOTW *** ! CALL WINDP2 (IDWMIN ,IDWMAX ,SIGPKD ,FPM , & ETOTW , & AC2 ,SPCSIG , WIND10 ) 40.00 EDML = MIN ( PWIND(10) , (GRAV**2 * ETOTW) / WIND10**4 ) EDML = MAX ( 1.E-25 , EDML ) ! ARG1 = ABS(PWIND(6)) ALPM = MAX( 0.0081, (PWIND(5) * (1./EDML)**ARG1) ) ! ENDIF ! ! *** Take into account depth effects for proportionality *** ! *** constant alpha through the nondimensional depth DND *** ! ALPMD = 0.0081 + ( 0.013 - 0.0081 ) * EXP ( -1. * DND ) ALPM = MIN ( 0.155 , MAX ( ALPMD , ALPM ) ) ! ! *** Calculate the limiting spectrum in terms of action density *** ! *** for the wind sea part (centered around the local wind *** ! *** direction). For conversion of f^-5 --> k^-3 and coefficients *** ! *** see Kitaigorodskii et al. 1975 *** ! DO IS = 1, ISSTOP TEMP1 = ALPM / ( 2. * KWAVE(IS,1)**3 * CGO(IS,1) ) ARG2 = MIN ( 2. , SIGPKD / SPCSIG(IS) ) 30.72 TEMP2 = EXP ( (-5./4.) * ARG2**4 ) ALIM1D = TEMP1 * TEMP2 / SPCSIG(IS) 30.72 DO IDDUM = IDWMIN, IDWMAX ID = MOD ( IDDUM - 1 + MDC, MDC ) + 1 THETA = SPCDIR(ID,1) 30.82 COSDIF = SPCDIR(ID,2)*CTW + SPCDIR(ID,3)*STW 40.41 ! ! For better convergence the first guess of the directional spreading 32.06 ! is modified in third generation mode. The new formulation better 32.06 ! fits the directional spreading of the deep water growth curves. 32.06 ! IF ((ITER.EQ.1).AND.(IGEN.EQ.3)) THEN 32.06 DIRDIS = 0.434917 * (MAX(0., COSDIF))**0.6 40.41 32.06 ELSE 32.06 DIRDIS = (2./PI) * COSDIF**2 40.41 END IF 32.06 ! ALIMW(ID,IS) = ALIM1D * DIRDIS AC2CEN = AC2(ID,IS,KCGRD(1)) IF ( AC2CEN .LE. ALIMW(ID,IS) ) THEN GROWW(ID,IS) = .TRUE. ELSE GROWW(ID,IS) = .FALSE. ENDIF ENDDO ! *** test output *** IF ( TESTFL .AND. ITEST .GE. 10 ) THEN WRITE(PRINTF,2002) IS, SPCSIG(IS), KWAVE(IS,1), CGO(IS,1) 30.72 2002 FORMAT(' WNDPAR: IS SPCSIG KWAVE CGO :',I3,3E12.4) 30.72 WRITE(PRINTF,2003) TEMP1, TEMP2, ARG2 2003 FORMAT(' WNDPAR: TEMP1 TEMP2 ARG2 :',3X,3E12.4) END IF ENDDO ! ! *** Calculate the wind input (linear term A and exponential *** ! *** term B) in wave generating conditions or disspation term *** ! *** if energy in bin is larger than limiting spectrum *** ! ! FACTA = PWIND(1) * PI * PWIND(9)**2 * PWIND(11)**2 / GRAV**2 ! DO IS = 1, ISSTOP SIGMA = SPCSIG(IS) 30.72 SIGTPI = SIGMA * TWOPI CINV = KWAVE(IS,1) / SIGMA FACTB = PWIND(2) * PWIND(9) * SIGMA / TWOPI 34.00 DO IDDUM = IDCMIN(IS), IDCMAX(IS) ID = MOD ( IDDUM - 1 + MDC, MDC ) + 1 DTHETA = SPCDIR(ID,1) - THETAW 30.82 COSDIF = SPCDIR(ID,2)*CTW + SPCDIR(ID,3)*STW 40.41 AC2CEN = AC2(ID,IS,KCGRD(1)) ! SWIND_EXP = 0. 40.13 SWIND_IMP = 0. 40.13 IF ( GROWW(ID,IS) ) THEN ! *** term A *** IF ( SIGMA .GE. ( 0.7 * SIGPKD ) ) THEN ADUM = FACTA * (WIND10 * COSDIF)**4 40.41 ADUM = MAX ( 0. , ADUM / SIGTPI ) ELSE ADUM = 0. END IF ! *** term B; Note that BDUM is multiplied with a factor 5 *** ! *** as in the DOLPHIN-B model *** ! BDUM = MAX( 0., ((WIND10 * CINV) * COSDIF-PWIND(3))) 40.41 BDUM = FACTB * BDUM * 5. SWIND_EXP = ADUM + BDUM * AC2CEN 40.13 ! ELSE IF ( .NOT. GROWW(ID,IS) .AND. AC2CEN .GT. 0. ) THEN ! ! *** for no energy dissipation outside the wind field *** ! *** TAUINV is set equal zero (as in the DOLPHIN-B model) *** ! IF ( COSDIF .LT. 0. ) THEN 40.41 TAUINV = 0. ELSE TAUINV = ( SIGMA**2 * WIND10 * ABS(COSDIF) ) / 40.41 & ( PWIND(4) * GRAV * TWOPI**2 ) ENDIF SWIND_EXP = TAUINV * ALIMW(ID,IS) SWIND_IMP = TAUINV ADUM = ALIMW(ID,IS) BDUM = TAUINV END IF ! ! *** store results in IMATDA and IMATRA *** ! IMATRA(ID,IS) = IMATRA(ID,IS) + SWIND_EXP IMATDA(ID,IS) = IMATDA(ID,IS) + SWIND_IMP IF (TESTFL) PLWNDS(ID,IS,IPTST) = SWIND_EXP 40.13 IF (TESTFL) PLWNDD(ID,IS,IPTST) = -1.*SWIND_IMP 40.13 GENC0(ID,IS,1) = GENC0(ID,IS,1) + SWIND_EXP 40.85 GENC1(ID,IS,1) = GENC1(ID,IS,1) - SWIND_IMP 40.85 ! ! *** test output *** ! ! Value of ITEST changed from 10 to 110 to reduce test output 40.13 IF ( TESTFL .AND. ITEST .GE. 110 ) THEN 40.13 WRITE(PRINTF,2004) IS, ID, GROWW(ID,IS), ADUM, BDUM 40.13 2004 FORMAT(' WNDPAR: IS ID GROWW ADUM BDUM :', 40.13 & 2I3,2X,L1,2X,2E12.4) END IF ENDDO ENDDO ! ! *** test output *** ! ! Value of ITEST changed from 10 to 60 to reduce test output 40.13 IF ( TESTFL .AND. ITEST .GE. 60 ) THEN 40.13 WRITE(PRINTF,*) WRITE(PRINTF,6051) IDWMIN, IDWMAX 6051 FORMAT(' WNDPAR : IDWMIN IDWMAX :',2I5) WRITE(PRINTF,6052) THETAW,WIND10,SIGPK,SIGPKD 6052 FORMAT(' WNDPAR : Tw U10 Spk Spk,d :',4E12.4) WRITE(PRINTF,7050) ETOTW, EDML, ALPM, ALPMD 7050 FORMAT(' WNDPAR: ETOW EDML ALPM ALPMD:',4E12.4) ENDIF ! RETURN ! end of subroutine WNDPAR END ! !**************************************************************** ! SUBROUTINE WINDP1 (WIND10 ,THETAW , & IDWMIN ,IDWMAX , & FPM ,UFRIC , & WX2 ,WY2 , & ANYWND ,SPCDIR , 40.00 & UX2 ,UY2 , #ifdef SWAN_MODEL & UFF ,UCDRAG , #endif & SPCSIG ,AC2 30.70 41.33 !ADC & ,NodeNumber 41.20 & ,GENC0 ,KWAVE 40.88 & ) ! !**************************************************************** ! USE OCPCOMM1 40.41 USE OCPCOMM2 40.41 USE OCPCOMM3 40.41 USE OCPCOMM4 40.41 USE SWCOMM1 40.41 USE SWCOMM2 40.41 USE SWCOMM3 40.41 USE SWCOMM4 40.41 USE SDSBABANIN !ADC USE WIND, ONLY: WindDrag 41.20 ! IMPLICIT NONE ! ! ! --|-----------------------------------------------------------|-- ! | Delft University of Technology | ! | Faculty of Civil Engineering | ! | Environmental Fluid Mechanics Section | ! | P.O. Box 5048, 2600 GA Delft, The Netherlands | ! | | ! | Programmers: The SWAN team | ! --|-----------------------------------------------------------|-- ! ! ! SWAN (Simulating WAves Nearshore); a third generation wave model ! Copyright (C) 1993-2017 Delft University of Technology ! ! This program is free software; you can redistribute it and/or ! modify it under the terms of the GNU General Public License as ! published by the Free Software Foundation; either version 2 of ! the License, or (at your option) any later version. ! ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! A copy of the GNU General Public License is available at ! http://www.gnu.org/copyleft/gpl.html#SEC3 ! or by writing to the Free Software Foundation, Inc., ! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ! ! ! 0. Authors ! ! 30.70: Nico Booij ! 30.82: IJsbrand Haagsma ! 32.06: Roeland Ris ! 40.00: Nico Booij (Nonstationary boundary conditions) ! 40.41: Marcel Zijlema ! 41.20: Casey Dietrich ! 41.33: Marcel Zijlema ! ! 1. Updates ! ! 20.64, : limited sector is now taken into account ! argument SPCDIR is added ! 30.70, Feb. 98: relative wind velocity is determined ! arguments UX2, UY2, SPCSIG added ! full common introduced ! 40.00, July 98: argument list changed: KCGRD removed ! 30.82, Oct. 98: Updated description of several variables ! 32.06, June 99: Reformulation of wind speed in terms of friction ! velocity for first and second generation ! 40.41, Oct. 04: common blocks replaced by modules, include files removed ! 41.20, Jun. 10: use ADCIRC's new sector-based wind drag ! 41.33, Mar. 12: extension drag coefficient based on 2nd order polynomial ! 41.33, Aug. 12: extension drag coefficient based on cross swell ! ! 2. Purpose ! ! Computation of parameters derived from the wind for several ! subroutines such as : ! SWIND1, SWIND2 SWIND3 ! CUTOFF ! ! Output of this subroutine : ! ! WIND10 , THETAW, IDWMIN, IDWMAX , UFRIC, FPM ! ! 3. METHOD ! ! a. For SWIND1 and SWIND2 : ! ! SIGMA_FPM = 0.13 * GRAV * 2 * PI / WIND10 ! ! b. For SWIND3 (wind input according to Snyder (1981) *** ! ! - wind friction velocity according to Wu (1982): ! * -3 ! U = UFRIC = wind10 sqrt( (0.8 + 0.065 wind10 ) 10 ) ! ! c. For SWIND4 (wind input according to Janssen 1991) ! ! - wind friction velocity: ! ! UFRIC = sqrt ( CDRAG) * U10 ! ! for U10 < 7.5 m/s -> CDRAG = 1.2873.e-3 ! ! else wind friction velocity according to Wu (1982): ! ! * -3 ! U = UFRIC = wind10 sqrt( (0.8 + 0.065 wind10 ) 10 ) ! ! d. ! The Pierson Moskowitz radian frequency for a fully developed ! sea state spectrum for all third generation wind input ! models is equal to: ! ! grav ! SIGMA_FPM = --------- ! 28 UFRIC ! ----------------------------------------------------------- ! ! 4. Argument variables ! ! i SPCDIR: (*,1); spectral directions (radians) 30.82 ! (*,2); cosine of spectral directions 30.82 ! (*,3); sine of spectral directions 30.82 ! (*,4); cosine^2 of spectral directions 30.82 ! (*,5); cosine*sine of spectral directions 30.82 ! (*,6); sine^2 of spectral directions 30.82 ! i SPCSIG: Relative frequencies in computational domain in sigma-space 30.82 ! REAL SPCDIR(MDC,6) 30.82 REAL SPCSIG(MSC) 30.82 ! ! IDWMIN Minimum counter for spectral wind direction ! IDWMAX Maximum counter for spectral wind direction ! IX Counter of gridpoints in x-direction ! IY Counter of gridpoints in y-direction ! MXC Maximum counter of gridppoints in x-direction ! MYC Maximum counter of gridppoints in y-direction ! KCGRD int, i Point index for grid point 30.21 ! MCGRD int, i Maximum counter of gridpoints in space 30.21 ! ICMAX int, i Maximum counter for the points of the molecule 30.21 ! ! REALS: ! --------- ! ! PI (3,14) ! GRAV Gravitational acceleration ! THETAW Mean direction of the relative wind vector ! WIND10 Velocity of the relative wind vector ! U10 wind velocity from SWANPREn ! WDIC mean wind direction from SWANPREn ! FPM PM frequency ! UFRIC Wind friction velocity ! ! one and more dimensional arrays: ! --------------------------------- ! WX2,WY2 2D Wind velocity array relative to a current ! PWIND 1D coefficients for wind expressions ! ANYWND 1D Indicator if wind input has to be taken into ! account for a bin ! ! 7. Common blocks used ! ! ! 5. SUBROUTINES CALLING ! ! SWOMPU ! ! 6. SUBROUTINES USED ! ! NONE ! ! 7. ERROR MESSAGES ! ! --- ! ! 8. REMARKS ! ! ! 9. STRUCTURE ! ! ------------------------------------------------------------ ! If constant wind then ! set parameters WIND10 and WDIC equal U10 and WDIC ! else compute ! wind velocity and mean wind direction ! ---------------------------------------------------------- ! compute the minimum and maximum counters for wind source term ! for which the wid input is active ! compute the Pierson Moskowitz frequency for IWIND = 1, 2 ! compute the wind friction velocity and PM freq. for IWIND = 3 ! ------------------------------------------------------------ ! End of the subroutine WINDP1 ! ------------------------------------------------------------ ! ! 10. SOURCE ! !*********************************************************************** ! !ADC INTEGER, OPTIONAL :: NodeNumber 41.20 !ADC! INTEGER IDWMIN ,IDWMAX 30.70 INTEGER IENT ,ID ,IDDUM, IS ! REAL WIND10 ,THETAW , 30.70 & UFRIC ,FPM ,CDRAG ,SDMEAN 30.70 ! REAL UREF ,UTL ,PP, QQ, RR 41.33 PARAMETER (UREF=31.5, PP=0.55, QQ=2.97, RR=-1.49) 41.33 REAL UREF1, UTL1, NSL1, NSH1, CS, CA,CB,CC, NA,NB,NC, 41.33 & UREF2, UTL2, NSL2, NSH2, CD,CE, ND,NE 41.33 PARAMETER (UREF1=27.5,NA=1.05,NB=1.25,NC=1.4,CA=0.7,CB=1.1,CC=6., 41.33 & NSL1=30.,NSH1=80.,CS=50., 41.33 & UREF2=54.,ND=2.3,NE=10.,CD=8.2,CE=2.5, 41.33 & NSL2=45.,NSH2=55.) 41.33 REAL A, B, C, D, E 41.33 REAL ETOTS ,EEX ,EEY , EAD ,SIGMA1, 41.33 & COSDIR,SINDIR,DETOT , FAC ,DSPR 41.33 ! REAL WX2(MCGRD) , & WY2(MCGRD) , & UX2(MCGRD) , 30.70 & UY2(MCGRD) , 30.70 #ifdef SWAN_MODEL & UFF(MCGRD) , MAI & UCDRAG(MCGRD) MAI #endif REAL AC2(MDC,MSC,MCGRD) 41.33 REAL KWAVE(MSC,MICMAX) REAL GENC0(MDC,MSC,MGENR) ! LOGICAL ANYWND(MDC) ! REAL AWX, AWY, RWX, RWY 30.70 SAVE IENT DATA IENT/0/ IF (LTRACE) CALL STRACE (IENT,'WINDP1') ! ! compute absolute wind velocity 30.70 IF (VARWI) THEN AWX = WX2(KCGRD(1)) AWY = WY2(KCGRD(1)) ELSE AWX = U10 * COS(WDIC) AWY = U10 * SIN(WDIC) ENDIF ! compute relative wind velocity 30.70 IF (ICUR.EQ.0) THEN RWX = AWX RWY = AWY ELSE RWX = AWX - UX2(KCGRD(1)) RWY = AWY - UY2(KCGRD(1)) ENDIF ! compute absolute value of relative wind velocity 30.70 WIND10 = SQRT(RWX**2+RWY**2) IF (WIND10.GT.0.) THEN THETAW = ATAN2 (RWY,RWX) THETAW = MOD ( (THETAW + PI2) , PI2 ) ELSE THETAW = 0. ENDIF ! ! *** compute the minimum and maximum counter for the active *** ! *** wind field : *** ! *** . *** ! *** IDWMAX =135 degrees .<--mean wind direction *** ! *** \ o o o | o o o . (THETAW) *** ! *** \ o o | o o . o *** ! *** \ o | o . o o *** ! *** \ | . o o o *** ! *** ----------------\-------------------- *** ! *** | \ o o o o *** ! *** | \ o o o *** ! *** | \ o o *** ! *** IDWMIN = 325 degrees *** ! *** *** ! ! move ThetaW to the right interval, shifting + or - 2*PI 20.64 SDMEAN = 0.5 * (SPCDIR(1,1) + SPCDIR(MDC,1)) 30.82 IF (THETAW .LT. SDMEAN - PI) THETAW = THETAW + 2.*PI IF (THETAW .GT. SDMEAN + PI) THETAW = THETAW - 2.*PI ! IF ( (THETAW - 0.5 * PI) .LE. SPCDIR(1,1) ) THEN 30.82 IF ( (THETAW + 1.5 * PI) .GE. SPCDIR(MDC,1) ) THEN IDWMIN = 1 ELSE IDWMIN = NINT ( (THETAW + 1.5*PI - SPCDIR(1,1)) / DDIR ) + 1 30.82 ENDIF ELSE IDWMIN = NINT ( (THETAW - 0.5*PI - SPCDIR(1,1)) / DDIR ) + 1 30.82 END IF ! IF ( (THETAW + 0.5 * PI) .GE. SPCDIR(MDC,1) ) THEN 30.82 IF ( (THETAW - 1.5 * PI) .LE. SPCDIR(1,1) ) THEN 30.82 IDWMAX = MDC ELSE IDWMAX = NINT ( (THETAW - 1.5 * PI - SPCDIR(1,1)) / DDIR ) + 1 30.82 ENDIF ELSE IDWMAX = NINT ( (THETAW + 0.5 * PI - SPCDIR(1,1)) / DDIR ) + 1 30.82 ENDIF ! IF ( IDWMIN .GT. IDWMAX) IDWMAX = MDC + IDWMAX ! ! *** determine for which bin the wind input is active *** ! *** initialize array for active wind input *** ! DO ID = 1, MDC ANYWND(ID) = .FALSE. ENDDO ! IF ( TESTFL .AND. ITEST .GE. 30 ) THEN WRITE(PRINTF,500) IDWMIN, IDWMAX 500 FORMAT(' WINDP1: IDWMIN IDWMAX :',2I15) ENDIF ! DO IDDUM = IDWMIN , IDWMAX ID = MOD ( IDDUM - 1 + MDC, MDC ) + 1 ANYWND(ID) = .TRUE. ! IF ( TESTFL .AND. ITEST .GE. 40 ) THEN WRITE(PRINTF,400) IDDUM, ID, ANYWND(ID) 400 FORMAT(' WINDP1: IDDUM ID ANYWND :',2I5,L4) ENDIF ! ENDDO ! ! --- determine directional spreading for cross swell ! IF ( IDRAG.EQ.3 ) THEN EEX = 0. EEY = 0. ETOTS = 0. DO ID = 1, MDC EAD = 0. DO IS = 1, MSC SIGMA1 = SPCSIG(IS) DETOT = SIGMA1**2 * AC2(ID,IS,KCGRD(1)) EAD = EAD + DETOT ENDDO ETOTS = ETOTS + EAD EEX = EEX + EAD * SPCDIR(ID,2) EEY = EEY + EAD * SPCDIR(ID,3) ENDDO IF ( ETOTS.GT.0. ) THEN COSDIR = EEX / ETOTS SINDIR = EEY / ETOTS FAC = MIN( 1., SQRT(COSDIR**2+SINDIR**2) ) DSPR = SQRT(2.-2.*FAC) ELSE DSPR = 0. ENDIF DSPR = DSPR * 180. / PI ENDIF ! ! *** compute the Pierson Moskowitz frequency *** ! IF ( IWIND .EQ. 1 .OR. IWIND .EQ. 2 ) THEN ! ! *** first and second generation wind wave model *** ! IF ( WIND10 .LT. PWIND(12) ) WIND10 = PWIND(12) FPM = 2. * PI * PWIND(13) * GRAV / WIND10 ! ! *** determine U friction in case predictor is obtained *** ! *** with second genaration wave model *** ! !ADC! added the default option to use the ADCIRC drag formulation !ADC IF ( IDRAG.EQ.0 ) THEN !ADC IF ( WIND10 .GT. 7.5 ) THEN !ADC CDRAG = WindDrag(DBLE(WIND10),NodeNumber) !ADC ELSE !ADC CDRAG = 0.0012873 !ADC ENDIF !ADC ENDIF IF ( IDRAG.EQ.1 ) THEN ! Wu (1982) drag formulation IF ( WIND10 .GT. 7.5 ) THEN CDRAG = ( 0.8 + 0.065 * WIND10 ) * 0.001 CDRAG = MIN ( CDCAP, CDRAG ) ELSE CDRAG = 0.0012873 ENDIF ELSE IF ( IDRAG.EQ.2 ) THEN ! Zijlema et al (2012) drag formulation UTL = WIND10/UREF CDRAG = ( PP + QQ*UTL + RR*UTL*UTL ) * 0.001 CDRAG = MIN ( CDCAP, CDRAG ) ELSE IF ( IDRAG.EQ.3 ) THEN ! drag based on swell UTL1 = WIND10/UREF1 UTL2 = WIND10/UREF2 ! IF ( DSPR.NE.CS ) THEN ! no swell, opposing or following swell A = NA B = NB C = NC D = ND E = NE IF ( DSPR.GT.NSL1 .AND. DSPR.LT.CS ) THEN FAC = (DSPR-NSL1) / (CS-NSL1) A = A + FAC * (CA-NA) B = B + FAC * (CB-NB) C = C + FAC * (CC-NC) IF ( DSPR.GT.NSL2 ) THEN FAC = (DSPR-NSL2) / (CS-NSL2) D = D + FAC * (CD-ND) E = E + FAC * (CE-NE) ENDIF ELSE IF ( DSPR.GT.CS .AND. DSPR.LT.NSH1 ) THEN FAC = (DSPR-NSH1) / (CS-NSH1) A = A + FAC * (CA-NA) B = B + FAC * (CB-NB) C = C + FAC * (CC-NC) IF ( DSPR.LT.NSH2 ) THEN FAC = (DSPR-NSH2) / (CS-NSH2) D = D + FAC * (CD-ND) E = E + FAC * (CE-NE) ENDIF ENDIF ELSE ! cross swell A = CA B = CB C = CC D = CD E = CE ENDIF ! CDRAG = MIN ( A+B*UTL1**C, D*(1-UTL2**E) ) CDRAG = MAX ( 0.7 , CDRAG ) CDRAG = CDRAG * 0.001 CDRAG = MIN ( CDCAP, CDRAG ) ENDIF ! UFRIC = SQRT ( CDRAG ) * WIND10 ! ! Reformulation of the wind speed in terms of friction velocity. 32.06 ! This formulation is based on Bouws (1986) and described in Delft 32.06 ! Hydraulics report H3515 (1999) 32.06 ! WIND10 = WIND10 * SQRT(((0.8 + 0.065 * WIND10) * 0.001) / 32.06 & ((0.8 + 0.065 * 15. ) * 0.001)) 32.06 ! ELSE IF (IWIND .GE. 3 ) THEN ! ! *** Calculate the wind friction velocity *** ! *** based on wind drag formulation *** ! *** apply cd-cap if appropriate *** ! !ADC! added the default option to use the ADCIRC drag formulation !ADC IF ( IDRAG.EQ.0 ) THEN !ADC IF ( WIND10 .GT. 7.5 ) THEN !ADC CDRAG = WindDrag(DBLE(WIND10),NodeNumber) !ADC ELSE !ADC CDRAG = 0.0012873 !ADC ENDIF !ADC UFRIC = SQRT ( CDRAG ) * WIND10 !ADC ENDIF IF ( IDRAG.EQ.1 ) THEN ! Wu (1982) drag formulation IF ( WIND10 .GT. 7.5 ) THEN CDRAG = ( 0.8 + 0.065 * WIND10 ) * 0.001 CDRAG = MIN ( CDCAP, CDRAG ) #if defined DRAGLIM_DAVIS && defined SWAN_MODEL CDRAG = MIN ( 0.0024, CDRAG ) !Limitig used by Davis maitane #endif ! this call is deprecated !ADC! CDRAG = WindDrag(DBLE(WIND10),NodeNumber) ELSE CDRAG = 0.0012873 ENDIF UFRIC = SQRT ( CDRAG ) * WIND10 ELSE IF ( IDRAG.EQ.2 ) THEN ! Zijlema et al (2012) drag formulation UTL = WIND10/UREF CDRAG = ( PP + QQ*UTL + RR*UTL*UTL ) * 0.001 CDRAG = MIN ( CDCAP, CDRAG ) UFRIC = SQRT ( CDRAG ) * WIND10 ELSE IF ( IDRAG.EQ.3 ) THEN ! drag based on swell UTL1 = WIND10/UREF1 UTL2 = WIND10/UREF2 ! IF ( DSPR.NE.CS ) THEN ! no swell, opposing or following swell A = NA B = NB C = NC D = ND E = NE IF ( DSPR.GT.NSL1 .AND. DSPR.LT.CS ) THEN FAC = (DSPR-NSL1) / (CS-NSL1) A = A + FAC * (CA-NA) B = B + FAC * (CB-NB) C = C + FAC * (CC-NC) IF ( DSPR.GT.NSL2 ) THEN FAC = (DSPR-NSL2) / (CS-NSL2) D = D + FAC * (CD-ND) E = E + FAC * (CE-NE) ENDIF ELSE IF ( DSPR.GT.CS .AND. DSPR.LT.NSH1 ) THEN FAC = (DSPR-NSH1) / (CS-NSH1) A = A + FAC * (CA-NA) B = B + FAC * (CB-NB) C = C + FAC * (CC-NC) IF ( DSPR.LT.NSH2 ) THEN FAC = (DSPR-NSH2) / (CS-NSH2) D = D + FAC * (CD-ND) E = E + FAC * (CE-NE) ENDIF ENDIF ELSE ! cross swell A = CA B = CB C = CC D = CD E = CE ENDIF ! CDRAG = MIN ( A+B*UTL1**C, D*(1-UTL2**E) ) CDRAG = MAX ( 0.7 , CDRAG ) CDRAG = CDRAG * 0.001 CDRAG = MIN ( CDCAP, CDRAG ) UFRIC = SQRT ( CDRAG ) * WIND10 ELSE IF (IDRAG.EQ.4) THEN ! Hwang (2011) drag formulation ! ! Sep 1 2010 : Wu formula replaced with Hwang eq 10. ! Source Hwang, "A note on the ocean surface roughness spectrum" ! Successfully tested with Hurricane Frances and Ivan in February 2011. ! Cap is necessary or CDRAG = 0 when winds reach > or = to 70 m/s. ! Capped at maximum Ustar for winds greater than 50.33 m/s. ! Hwang formulation added to wave age calculation on 2/15/11 in swanout1.f. ! Dec 28 2016: CDFAC added to (optionally) counter bias in wind speeds ! : Important: this is applied *after* the cap on UFRIC CDRAG = (-0.016*WIND10**2 + 0.967*WIND10 + 8.058) * 0.0001 UFRIC = SQRT ( CDRAG ) * WIND10 #ifdef SWAN_MODEL UCDRAG = CDRAG #endif IF ( WIND10.GT.50.33 ) UFRIC = 2.02558 UFRIC = MIN(USCAP,UFRIC) * SQRT(CDFAC) ! ELSE IF (IDRAG.EQ.5) THEN ! Fan et al (2012) drag formulation CALL SURF_ROUGH_FAN (FPI, WIND10, UFRIC, CDRAG) FPM = GRAV / ( 28.0 * UFRIC ) ELSE IF (IDRAG.EQ.6) THEN ! ECMWF drag formulation CALL SURF_ROUGH_ECMWF(WIND10, UFRIC, GENC0, SPCSIG, & KWAVE, CDRAG) END IF ! ! *** adapted wind friction velocity and PM-frequency *** ! IF (IDRAG.LT.5) THEN UFRIC = MAX ( 1.E-15 , UFRIC) FPM = GRAV / ( 28.0 * UFRIC ) #ifdef SWAN_MODEL UFF = UFRIC #endif END IF END IF ! ! *** test output *** ! IF ( TESTFL .AND. ITEST .GE. 50 ) THEN WRITE(PRINTF,6050) KCGRD(1), MDC, MCGRD, IWIND 30.21 6050 FORMAT(' WINDP1:INDEX MDC MCGRD IWND:',4I5) WRITE(PRINTF,6052) THETAW,WIND10,WDIC,U10 6052 FORMAT(' : THAW WIND10 WDIC U10 :',4E12.4) WRITE(PRINTF,6054) GRAV, PI, DDIR, VARWI 6054 FORMAT(' : GRAV PI DDIR VARWI :',3E12.4,L6) WRITE(PRINTF,6056) IDWMIN,IDWMAX,FPM, UFRIC 6056 FORMAT(' : IDWMIN IDWMAX FPM UFR:',2I4,2E12.4) WRITE(PRINTF,*) END IF ! RETURN ! end of subroutine WINDP1 END ! !**************************************************************** ! SUBROUTINE WINDP2 (IDWMIN ,IDWMAX ,SIGPKD ,FPM , & ETOTW , & AC2 ,SPCSIG , 40.00 & WIND10 ) 30.70 ! !**************************************************************** ! USE OCPCOMM1 40.41 USE OCPCOMM2 40.41 USE OCPCOMM3 40.41 USE OCPCOMM4 40.41 USE SWCOMM1 40.41 USE SWCOMM2 40.41 USE SWCOMM3 40.41 USE SWCOMM4 40.41 ! ! ! --|-----------------------------------------------------------|-- ! | Delft University of Technology | ! | Faculty of Civil Engineering | ! | Environmental Fluid Mechanics Section | ! | P.O. Box 5048, 2600 GA Delft, The Netherlands | ! | | ! | Programmers: The SWAN team | ! --|-----------------------------------------------------------|-- ! ! ! SWAN (Simulating WAves Nearshore); a third generation wave model ! Copyright (C) 1993-2017 Delft University of Technology ! ! This program is free software; you can redistribute it and/or ! modify it under the terms of the GNU General Public License as ! published by the Free Software Foundation; either version 2 of ! the License, or (at your option) any later version. ! ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! A copy of the GNU General Public License is available at ! http://www.gnu.org/copyleft/gpl.html#SEC3 ! or by writing to the Free Software Foundation, Inc., ! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ! ! ! 0. Authors ! ! 30.72: IJsbrand Haagsma ! 30.70: Nico Booij ! 40.41: Marcel Zijlema ! ! 1. Updates ! ! 20.72, Jan. 96: Integration modified, using FRINTF, FRINTH ! and PWTAIL(6) ! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN ! 30.70, Feb. 98: full common introduced, argument list changed ! ISFPM changed (in case of very high value of FPM) ! 40.00, July 98: argument list changed: KCGRD removed ! 40.41, Oct. 04: common blocks replaced by modules, include files removed ! ! 2. Purpose ! ! Computation of wind sea energy spectrum for the second ! generation wind growth model. Output of the subroutine ! is ETOTW (total wind sea energy spectrum) ! ! 3. Method ! ! Compute the wind sea energy spectrum : ETOTW ! ! +90 inf ! | | ! ETOTW = | | E(s,d) ds dd ! -90 0.7*FPM ! ! Compute the total energy density ETOTW for F > 0.7 FPM ! ! ^ | ! E()| | * ! | * * ! | * ! | * | * / ETOTW(e) ! | | o * / ! | * | o o/o * ! | | o o o o o * ! | * | o o o o o o o o* ! 0---------------|--------------------------- ! 0 0.7*FPM SIGMA --> s ! ! SIGMA MAX ! | ! ETOTW(d) = | E(s,d)ds ISFPM = FPM ! 0.7 FPM ! ! and over the interval +/- 90 degrees (according to ! the mean wind direction as computed in WINDP1 ) ! ! IDWMAX = "135" ! o ! ETOTW = o ETOTW(d)dd ! ! IDWMIN ="325" ! ! ! 4. Argument variables ! ! SPCSIG: Relative frequencies in computational domain in sigma-space 30.72 ! REAL SPCSIG(MSC) 30.72 ! ! ISFPM Counter in point just for the Pierson Moskowitz ! frequency ! IDWMIN Minimum counter for spectral wind direction ! IDWMAX Maximum counter for spectral wind direction ! IX Counter of gridpoints in x-direction ! IY Counter of gridpoints in y-direction ! IS Counter of relative frequency band ! ID Counter of directional distribution ! MXC Maximum counter of gridppoints in x-direction ! MYC Maximum counter of gridppoints in y-direction ! MSC Maximum counter of relative frequency ! MDC Maximum counter of directional distribution ! ! REALS: ! --------- ! ! DAC2 Difference in action density between two neighbour ! points (IS and IS+1) ! DD Directional band width ! DSIG Distance between FPM and next frequency ! DUM_DS Distance betwee the relative frequency and the ! next frequency value : LOGSIG(IS+1) - FPM ! ETOT_ Dummy variables to compute the total energy ! ETOTW Total energy density over the sea spectrum, i.e. ! a > apm and a 180 degrees spectral interval ! CNETOT contribution of high frequency tail over the ! full spectrum ! FPM Pierson Moskowitz frequency ! SIGFAC Relative difference between two frequencies ! THETAW Mean direction of the relative wind vector ! WIND10 Velocity of the relative wind vector ! ! one and more dimensional arrays: ! --------------------------------- ! AC2 4D Action density as function of D,S,X,Y and T ! ! 7. Common blocks used ! ! ! 5. SUBROUTINES CALLING ! ! SWOMPU, SOURCE ! ! 6. SUBROUTINES USED ! ! NONE ! ! 7. ERROR MESSAGES ! ! --- ! ! 8. REMARKS ! ! ! 9. STRUCTURE ! ! ------------------------------------------------------------ ! Determine the counter for the FPM frequency ! integrate the wind sea energy spectrum ! add the energy of the high frequency tail to the spectrum ! ------------------------------------------------------------ ! ! 10. SOURCE ! !*********************************************************************** ! ! INTEGER IDWMIN ,IDWMAX , & IDDUM ,ID ,IS ,ISFPM ! REAL ETOTW ,FPM ,SIG ,ATOTD REAL AC2(MDC,MSC,MCGRD) ! SAVE IENT DATA IENT/0/ IF (LTRACE) CALL STRACE (IENT,'WINDP2') ! ! *** compute wind sea energy spectrum for IS > 0.7 FPM *** ! *** minimum FPM is equal : 2 * pi * 0.13 * grav / pwind(12) *** ! *** is equal 8 rad/s = 1.27 Hz *** ! ISFPM = MSC 30.70 FACINT = 0. 30.70 DO IS = 1, MSC SIG = SPCSIG(IS) 30.72 IF (FRINTH * SIG .GT. (0.7 * FPM) ) THEN ISFPM = IS FACINT = (FRINTH - 0.7*FPM/SIG) / (FRINTH - 1./FRINTH) GOTO 11 END IF ENDDO 11 CONTINUE ! ! *** calculate the energy in the wind sea part of the spectrum *** ! *** from ISFPM. *** ! ETOTW = 0. DO IS = ISFPM, MSC SIG = SPCSIG(IS) 30.72 ATOTD = 0. DO IDDUM = IDWMIN, IDWMAX ID = MOD ( IDDUM - 1 + MDC, MDC ) + 1 ATOTD = ATOTD + AC2(ID,IS,KCGRD(1)) 30.21 ENDDO IF (IS.EQ.ISFPM) THEN ETOTW = ETOTW + FACINT * FRINTF * SIG**2 * DDIR * ATOTD ELSE ETOTW = ETOTW + FRINTF * SIG**2 * DDIR * ATOTD ENDIF ENDDO ! add high-frequency tail: ETOTW = ETOTW + PWTAIL(6) * SIG**2 * DDIR * ATOTD ! ! *** test output *** ! IF ( TESTFL .AND. ITEST .GE. 70 ) THEN WRITE(PRINTF,*) WRITE(PRINTF,6050) IWIND,IDWMIN,IDWMAX, ISFPM, ETOTW 6050 FORMAT(' WINDP2: IWND IDWMIN IDWMAX ISFPM ETOTW:',4I6,1X,E12.4) END IF ! RETURN ! end of subroutine WINDP2 END ! !******************************************************************** ! SUBROUTINE WINDP3 (ISSTOP ,ALIMW ,AC2 , & GROWW ,IDCMIN ,IDCMAX ) 40.41 ! !**************************************************************** ! USE SWCOMM3 40.41 USE SWCOMM4 40.41 USE OCPCOMM4 40.41 ! ! ! --|-----------------------------------------------------------|-- ! | Delft University of Technology | ! | Faculty of Civil Engineering | ! | Environmental Fluid Mechanics Section | ! | P.O. Box 5048, 2600 GA Delft, The Netherlands | ! | | ! | Programmers: The SWAN team | ! --|-----------------------------------------------------------|-- ! ! ! SWAN (Simulating WAves Nearshore); a third generation wave model ! Copyright (C) 1993-2017 Delft University of Technology ! ! This program is free software; you can redistribute it and/or ! modify it under the terms of the GNU General Public License as ! published by the Free Software Foundation; either version 2 of ! the License, or (at your option) any later version. ! ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! A copy of the GNU General Public License is available at ! http://www.gnu.org/copyleft/gpl.html#SEC3 ! or by writing to the Free Software Foundation, Inc., ! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ! ! ! 1. UPDATE ! ! 40.41, Oct. 04: common blocks replaced by modules, include files removed ! ! 2. PURPOSE ! ! Reduce the energy density in spectral direction direct after ! solving the tri-diagonal matrix if the energy density level is ! larger then the upper bound limit given by a Pierson Moskowitz ! spectrum. This is only carried out if a particular wave ! component is 'growing'. ! ! If the energy density in a bin is larger than the upper bound ! limit (for instance when crossing wind seas are present) then ! the energy density level is a lower limit. ! ! 3. METHOD ! ! The upper bound limit is given by: ! ! 2 -4 ! ALPHA g sigma 2 2 ! A(s,t) = ----------- exp ( -5/4 (----) ) --- cos ( d - dw ) ! sigma^6 FPM pi ! ! in which ALPHA is wind sea dependent (see subroutine ! SWIND2) : ! ! / C60 \ ! | / E_windsea g^2 \ | ! ALPHA = MIN | 0.0081 , C50 | ------------ | | ! \ \ U_10^4 / / ! ! ! 4. PARAMETERLIST ! ! INTEGERS : ! ---------- ! IX IY Grid point in geographical space ! MDC ,MSC Counters in spectral space ! ISSTOP Maximum frequency that fall within a sweep ! ! ARRAYS: ! ------- ! AC2 4D Action density as funciton of x,y,s,t ! ALIMW 2D Contains the action density upper bound ! limit regardind spectral action density per ! spectral bin (A(s,t)) ! GROWW 2D Logical array which determines is there is ! a) generation ( E < E_lim -> .TRUE. ) or ! b) dissipation ( E > E_lim -> .FALSE. ) ! IDCMIN 1D Frequency dependent minimum counter ! IDCMAX 1D Frequency dependent minimum counter ! ! 5. SUBROUTINES CALLING ! ! SWOMPU ! ! 6. SUBROUTINES USED ! ! NONE ! ! 7. Common blocks used ! ! ! 8. REMARKS ! ! ! 9. STRUCTURE ! ! ------------------------------------------------------------ ! For every spectral bin ! If wave input for a bin is true (GROWW(ID,IS) = .TRUE.) and ! if wave action is larger then maximum then reduce ! wave action to limit its value to upper boundary limit ! else if no growth (see subroutine SWIND1 and SWIND2), check ! if the energy density level in the wind sea part ! of the spectrum is lower than upper bound limit. If so ! set energy level equal to upper bound limit ! endif ! ------------------------------------------------------------ ! End of the subroutine WINDP3 ! ------------------------------------------------------------ ! ! 10. SOURCE ! !*********************************************************************** ! INTEGER IS ,ID ,ISSTOP ,IDDUM ! INTEGER IDCMIN(MSC) , & IDCMAX(MSC) ! REAL AC2CEN ! REAL AC2(MDC,MSC,MCGRD) , & ALIMW(MDC,MSC) ! LOGICAL GROWW(MDC,MSC) ! SAVE IENT DATA IENT/0/ IF (LTRACE) CALL STRACE (IENT,'WINDP3') ! ! *** limit the action density spectrum *** ! DO IS = 1, ISSTOP DO IDDUM = IDCMIN(IS) , IDCMAX(IS) ID = MOD ( IDDUM - 1 + MDC, MDC ) + 1 AC2CEN = AC2(ID,IS,KCGRD(1)) IF ( GROWW(ID,IS) .AND. AC2CEN .GT. ALIMW(ID,IS) ) & AC2(ID,IS,KCGRD(1)) = ALIMW(ID,IS) IF ( .NOT. GROWW(ID,IS) .AND. AC2CEN .LT. ALIMW(ID,IS) ) & AC2(ID,IS,KCGRD(1)) = ALIMW(ID,IS) ! IF (TESTFL .AND. ITEST .GE. 50) THEN WRITE(PRINTF,300) IS,ID,GROWW(ID,IS),AC2CEN,ALIMW(ID,IS) 300 FORMAT(' WINDP3 : IS ID GROWW AC2CEN ALIM:',2I4,L4,2E12.4) END IF ! ENDDO ENDDO ! ! *** test output *** ! IF (TESTFL .AND. ITEST .GE. 50) THEN WRITE(PRINTF,4000) KCGRD(1),ISSTOP,MSC,MDC,MCGRD 4000 FORMAT(' WINDP3 : POINT ISSTOP MSC MDC MCGRD :',5I5) END IF ! RETURN ! end of subroutine WINDP3 END ! !**************************************************************** ! SUBROUTINE SWIND0 (IDCMIN ,IDCMAX ,ISSTOP , & SPCSIG ,THETAW ,ANYWND , & UFRIC ,FPM ,PLWNDS , & IMATRA ,SPCDIR ,GENC0 , & KWAVE ) ! !**************************************************************** ! USE SWCOMM3 40.41 USE SWCOMM4 40.41 USE OCPCOMM4 40.41 ! IMPLICIT NONE ! ! ! --|-----------------------------------------------------------|-- ! | Delft University of Technology | ! | Faculty of Civil Engineering | ! | Environmental Fluid Mechanics Section | ! | P.O. Box 5048, 2600 GA Delft, The Netherlands | ! | | ! | Programmers: The SWAN team | ! --|-----------------------------------------------------------|-- ! ! ! SWAN (Simulating WAves Nearshore); a third generation wave model ! Copyright (C) 1993-2017 Delft University of Technology ! ! This program is free software; you can redistribute it and/or ! modify it under the terms of the GNU General Public License as ! published by the Free Software Foundation; either version 2 of ! the License, or (at your option) any later version. ! ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! A copy of the GNU General Public License is available at ! http://www.gnu.org/copyleft/gpl.html#SEC3 ! or by writing to the Free Software Foundation, Inc., ! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ! ! ! 0. Authors ! ! 30.72: IJsbrand Haagsma ! 30.82: IJsbrand Haagsma ! 40.41: Marcel Zijlema ! 40.85: Marcel Zijlema ! ! 1. Updates ! ! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN ! 30.82, Oct. 98: Updated description of several variables ! 40.41, Aug. 04: COS(THETA-THETAW) replaced by sumrule to make it cheaper ! 40.41, Oct. 04: common blocks replaced by modules, include files removed ! 40.85, Aug. 08: store wind input for output purposes ! ! 2. PURPOSE ! ! Computation of the source term for the wind input for a ! third generation wind growth model: ! ! 1) Linear wind input term according to Cavaleri and ! Malanotte-Rizzoli (1981) ! ! 3. METHOD ! ! To ensure wave growth when no wave energy is present in the ! numerical model a linear growth term is used (see ! Cavaleri and Malanotte-Rizzoli 1981). Contributions for ! frequencies lower then FPM have been eliminated buy aa filter : ! ! -3 ! 1.5*10 * 4 sigma -4 ! A = Sin = ------- { U max[0, (cos(d - dw )] } exp{-(-----) } / Jac ! 2 FPM ! g ! ! With Jac = Jacobian = 2 pi sigma ! ! The Pierson Moskowitz radian frequency for a fully developed ! sea state spectrum is as follows (computed in WINDP2 : ! ! 1 g ! FPM = ---- --------- * 2 pi ! 2 pi 28 UFRIC ! ! ! 4. Argument variables ! ! i SPCDIR: (*,1); spectral directions (radians) 30.82 ! (*,2); cosine of spectral directions 30.82 ! (*,3); sine of spectral directions 30.82 ! (*,4); cosine^2 of spectral directions 30.82 ! (*,5); cosine*sine of spectral directions 30.82 ! (*,6); sine^2 of spectral directions 30.82 ! i SPCSIG: Relative frequencies in computational domain in sigma-space 30.72 ! REAL SPCDIR(MDC,6) 30.82 REAL SPCSIG(MSC) 30.72 ! ! IDWMIN Minimum counter for spectral wind direction ! IDWMAX Maximum counter for spectral wind direction ! IS Counter of relative frequency band ! ID Counter of directional distribution ! MSC Maximum counter of relative frequency ! MDC Maximum counter of directional distribution ! ! REALS: ! --------- ! FPM Pierson Moskowitz frequecy (WAM) ! GRAV Gravitatiuonal acceleration ! SWINEO Coefficient stored in IMATRA ! THETA Spectral direction ! THETAW Mean direction of the relative wind vector ! UFRIC Wind friction velocity ! ! one and more dimensional arrays: ! --------------------------------- ! IMATRA 2D Coefficients of right hand side of matrix ! ANYWND 1D Determines the number of directional bins ! that fall within the wind sea part of the ! spectrum ( thetaw +- 90 degrees) ! IDCMIN 1D Frequency dependent counter ! IDCMAX 1D Frequency dependent counter ! ! 7. Common blocks used ! ! ! 5. SUBROUTINES CALLING ! ! SOURCE ! ! 6. SUBROUTINES USED ! ! --- ! ! 7. ERROR MESSAGES ! ! --- ! ! 8. REMARKS ! ! E. Rogers, May 16 2012: I have noticed that for computations beyond 1 Hz, ! the stress contribution from this term could be quite large. This is because ! 1) it is a linear term, so the smallness of E(f) in the tail does not make Sin(f) smaller ! 2) the linear term is flat, plotted as a function of frequency ! 3) the stress contribution is something like tau(f)=integral of Sin(f) / C df , so tau ! increases with frequency ! And some other remarks: ! 1) Equation (5) in Cavaleri and Rizzoli (JGR 1981) does not make sense to me: ! the units in particular do not seem correct. ! 2) I cannot see how equation (4) in Cavaleri and Rizzoli (JGR 1981) follows from ! his source material. In particular, where does the "80" come from? There ! appears to be at least one step missing in getting from Luigi's source material ! to his equation (4). ! 3) I cannot see how eq (5) in Tolman (JPO 1992) follows from his source material ! which is Luigi's paper. Again, there appears to be at least one step missing, ! this time in getting from Luigi's eq (5) to Tolman's eq (10). ! 4) In particular, Luigi's eq (5) is proportional to (sigma)/(g^2*k^2), which is ! for deep water sigma^-3. This dependency is completely missing in Tolman's eq (10)! ! ! ....so what to do? ! I don't want to make drastic changes to this formulation without knowing more about the ! highly suspicious genesis of the formulation, even though I have a strong feeling now ! that it is wrong. So, for now, I'll just apply the sigma^-3 drop-off for frequencies ! beyond 1 Hz. Thus, it won't affect most simulations, which stop at 1 Hz. ! Source for the sigma^-3 drop-off: Luigi eq (5) ! Source for the 1 Hz application/bending point: mostly arbitrary, intended to have ! zero effect on typical simulations (which stop at 1 Hz) ! Alternate bending point: fpm, though this would have a much more noticeable effect on ! typical simulations. ! ! 9. STRUCTURE ! ! ------------------------------------------------------------ ! For every spectral bin that faal witin the sweep considered: ! compute the source term Sinl and store the results in array ! -------------------------------------------------------------- ! End of the subroutine SWIND0 ! -------------------------------------------------------------- ! ! 10. SOURCE ! !*********************************************************************** ! INTEGER IENT, IDDUM ,ID ,IS ,ISSTOP ! REAL FPM ,UFRIC ,THETA ,THETAW , & SWINEA ,SIGMA ,TEMP1 ,TEMP2 , & CTW ,STW ,COSDIF , 40.41 & TEMP3 ,FILTER ,ARGU ,REDUC ,FREQ1 REAL TAUX ,TAUY ,CINV2 ,CTH ,STH ! REAL IMATRA(MDC,MSC) , & PLWNDS(MDC,MSC,NPTST) 40.00 REAL GENC0(MDC,MSC,MGENR) 40.85 REAL KWAVE(MSC,MICMAX) ! INTEGER IDCMIN(MSC) , & IDCMAX(MSC) ! LOGICAL ANYWND(MDC) ! SAVE IENT DATA IENT/0/ IF (LTRACE) CALL STRACE (IENT,'SWIND0') ! ! *** calculate linear wind input term *** ! CTW = COS(THETAW) 40.41 STW = SIN(THETAW) 40.41 FPM = GRAV / ( 28.0 * UFRIC ) TEMP1 = PWIND(31) / ( GRAV**2 * 2. * PI ) 7/MAR DO IS = 1, ISSTOP SIGMA = SPCSIG(IS) 30.72 ! ! **** ARGU = FPM / SIGMA *** ! **** the value of ARGU was change for MIN () because for *** ! **** values of fpm/sigma too small could be some problems *** ! **** with some computers to handle small numbers *** ARGU = MIN (2., FPM / SIGMA) 30.00 FILTER = EXP ( - ARGU**4 ) 20.87 ! note that SIGMA below is not in eq A-2 of Ris (1997) ! thus, we are calculating "A/sigma" here, not "A" TEMP2 = TEMP1 / SIGMA DO IDDUM = IDCMIN(IS), IDCMAX(IS) ID = MOD ( IDDUM - 1 + MDC, MDC ) + 1 IF ( ANYWND(ID) .AND. SIGMA .GE. (0.7 * FPM) ) THEN THETA = SPCDIR(ID,1) 30.82 COSDIF = SPCDIR(ID,2)*CTW + SPCDIR(ID,3)*STW 40.41 TEMP3 = ( UFRIC * MAX( 0. , COSDIF))**4 40.41 SWINEA = MAX( 0. , TEMP2 * TEMP3 * FILTER ) ! see above remarks for rationale on reduction above 1 Hz FREQ1 = SIGMA/PI2 40.88 IF ( FREQ1.GT.1. ) THEN REDUC = FREQ1**(-3) ELSE REDUC = 1. ENDIF SWINEA = REDUC * SWINEA IMATRA(ID,IS) = IMATRA(ID,IS) + SWINEA IF(TESTFL) PLWNDS(ID,IS,IPTST) = SWINEA 40.85 40.00 GENC0(ID,IS,1) = GENC0(ID,IS,1) + SWINEA 40.85 ! ! *** test output *** ! IF (ITEST .GE. 80 .AND. TESTFL ) & WRITE (PRTEST, 333) ID, IS, FILTER, SWINEA 333 FORMAT (' ID IS FILTER WIND SOURCE (IMATRA)', & 2I4, 1X, 2(1X,E11.4)) END IF ENDDO ENDDO ! ! calculate stress (test point only) IF (ITEST.GE.40.AND.TESTFL) THEN TAUX=0. TAUY=0. DO IS = 1, MSC CINV2=KWAVE(IS,1)/SPCSIG(IS) DO ID = 1, MDC CTH = SPCDIR(ID,2) ! new local variable = cos(theta) STH = SPCDIR(ID,3) ! new local variable = sin(theta) ! SPCSIG(IS) replaces "EN" as used in SWIND_Donelan: Thus, factor AC2 is removed, EN=sig*ac2 TAUX =TAUX +CTH*CINV2*PLWNDS(ID,IS,IPTST) & *SPCSIG(IS)*DDIR*FRINTF*SPCSIG(IS) TAUY =TAUY +STH*CINV2*PLWNDS(ID,IS,IPTST) & *SPCSIG(IS)*DDIR*FRINTF*SPCSIG(IS) ENDDO ENDDO TAUX=TAUX*PWIND(17)*GRAV TAUY=TAUY*PWIND(17)*GRAV WRITE(PRINTF,*)'SWIND0: TAUX,TAUY = ',TAUX,TAUY ENDIF ! ! *** test output *** ! IF (ITEST.GE.60.AND.TESTFL) THEN WRITE(PRINTF,400) KCGRD(1), THETAW*180./PI 400 FORMAT(' SWIND0: POINT THETAW :',I5,E12.4) WRITE(PRINTF,500) TEMP1, FPM, UFRIC 500 FORMAT(' SWIND0: TEMP1 FPM UFRC :',3E12.4) WRITE(PRINTF,*) IF (ITEST.GE. 120.AND.TESTFL) THEN 24/MAR DO IS = 1, ISSTOP DO IDDUM = IDCMIN(IS), IDCMAX(IS) ID = MOD ( IDDUM - 1 + MDC, MDC ) + 1 WRITE(PRINTF,100) IS,ID,ANYWND(ID) 100 FORMAT(' IS ID ANYWND : ', 2I5,1X,L2) ENDDO ENDDO ENDIF END IF ! RETURN ! end of subroutine SWIND0 END ! !**************************************************************** ! SUBROUTINE SWIND3 (SPCSIG ,THETAW , & KWAVE ,IMATRA ,GENC0 , & IDCMIN ,IDCMAX ,AC2 ,UFRIC , & FPM ,PLWNDS ,ISSTOP ,SPCDIR ,ANYWND) ! !**************************************************************** ! USE SWCOMM3 40.41 USE SWCOMM4 40.41 USE OCPCOMM4 40.41 !ESMF USE M_GENARR, ONLY: SAVE_SINBAC, SINBAC ! IMPLICIT NONE ! ! ! --|-----------------------------------------------------------|-- ! | Delft University of Technology | ! | Faculty of Civil Engineering | ! | Environmental Fluid Mechanics Section | ! | P.O. Box 5048, 2600 GA Delft, The Netherlands | ! | | ! | Programmers: The SWAN team | ! --|-----------------------------------------------------------|-- ! ! ! SWAN (Simulating WAves Nearshore); a third generation wave model ! Copyright (C) 1993-2017 Delft University of Technology ! ! This program is free software; you can redistribute it and/or ! modify it under the terms of the GNU General Public License as ! published by the Free Software Foundation; either version 2 of ! the License, or (at your option) any later version. ! ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! A copy of the GNU General Public License is available at ! http://www.gnu.org/copyleft/gpl.html#SEC3 ! or by writing to the Free Software Foundation, Inc., ! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ! ! ! 0. Authors ! ! 30.72: IJsbrand Haagsma ! 30.82: IJsbrand Haagsma ! 40.41: Marcel Zijlema ! 40.85: Marcel Zijlema ! ! 1. Updates ! ! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN ! 30.82, Oct. 98: Updated description of several variables ! 40.41, Aug. 04: COS(THETA-THETAW) replaced by sumrule to make it cheaper ! 40.41, Oct. 04: common blocks replaced by modules, include files removed ! 40.85, Aug. 08: store wind input for output purposes ! ! 2. PURPOSE ! ! Computation of the source term for the wind input for a ! third generation wind growth model: ! ! 1) Exponential input term, (Snyder et al. 1981, which ! expression has been modified by Komen et al. 1984). ! ! This input term should be combined with the dissipation ! term of Komen et al. (1984) ! ! 3. METHOD ! ! The exponential term used is taken from Snyder et al. (1981) ! and Komen et al. (1984): ! ! Sin (s,d) = B*E(s,d) ! e * ! 28 U cos( d - dw ) ! B = max(0. , (0.25 rhoaw( ------------------- -1 ) )) sigma ! sigma / kwave ! ! with : ! ! * -3 ! U = UFRIC = wind10 sqrt( (0.8 + 0.065 wind10 ) 10 ) ! ! UFRIC is computed in WINDP1 ! ! ! The Pierson Moskowitz radian frequency for a fully developed ! sea state spectrum is as follows (computed in WINDP1): ! ! 1 g ! FPM = ---- --------- * 2 pi ! 2 pi 28 UFRIC ! ! ! 4. Argument variables ! ! i SPCDIR: (*,1); spectral directions (radians) 30.82 ! (*,2); cosine of spectral directions 30.82 ! (*,3); sine of spectral directions 30.82 ! (*,4); cosine^2 of spectral directions 30.82 ! (*,5); cosine*sine of spectral directions 30.82 ! (*,6); sine^2 of spectral directions 30.82 ! i SPCSIG: Relative frequencies in computational domain in sigma-space 30.72 ! REAL SPCDIR(MDC,6) 30.82 REAL SPCSIG(MSC) 30.72 ! ! IS Counter of relative frequency band ! ID Counter of directional distribution ! MSC Maximum counter of relative frequency ! MDC Maximum counter of directional distribution ! ! REALS: ! --------- ! FPM Pierson Moskowitz frequecy (WAM) ! THETA Spectral direction ! THETAW Mean direction of the relative wind vector ! UFRIC Wind friction velocity ! ! one and more dimensional arrays: ! --------------------------------- ! KWAVE 2D Wavenumber ! LOGSIG 1D Logaritmic distribution of frequency ! IMATRA 2D Coefficients of right hand side of matrix ! PWIND 1D Wind coefficients ! IDCMIN 1D Frequency dependent counter ! IDCMAX 1D Frequency dependent counter ! ANYWND 1D Wind input for bin considered ! ! 7. Common blocks used ! ! ! 5. SUBROUTINES CALLING ! ! SOURCE ! ! 6. SUBROUTINES USED ! ! --- ! ! 7. ERROR MESSAGES ! ! --- ! ! 8. REMARKS ! ! --- ! ! 9. STRUCTURE ! ! ------------------------------------------------------------ ! For every spectral bin that fall in teh sweep considered ! compute the source term and store the results in IMATRA ! -------------------------------------------------------------- ! End of the subroutine SWIND3 ! -------------------------------------------------------------- ! ! 10. SOURCE ! !*********************************************************************** ! INTEGER IENT, IDDUM ,ID ,IS ,ISSTOP ! REAL FPM ,UFRIC ,THETA ,THETAW,SIGMA ,SWINEB,TEMP1, & CTW ,STW ,COSDIF, 40.41 & TEMP2 ,TEMP3 ,CINV ! REAL AC2(MDC,MSC,MCGRD) , & IMATRA(MDC,MSC) , & KWAVE(MSC,MICMAX) , & PLWNDS(MDC,MSC,NPTST) 40.00 REAL :: GENC0(MDC,MSC,MGENR) 40.85 ! INTEGER IDCMIN(MSC) , & IDCMAX(MSC) ! LOGICAL ANYWND(MDC) ! !/T LOGICAL IMP_EXP ! SAVE IENT DATA IENT/0/ IF (LTRACE) CALL STRACE (IENT,'SWIND3') ! CTW = COS(THETAW) 40.41 STW = SIN(THETAW) 40.41 TEMP1 = 0.25 * PWIND(9) TEMP2 = 28.0 * UFRIC DO IS = 1, ISSTOP SIGMA = SPCSIG(IS) 30.72 CINV = KWAVE(IS,1) / SIGMA TEMP3 = TEMP2 * CINV DO IDDUM = IDCMIN(IS), IDCMAX(IS) ID = MOD ( IDDUM - 1 + MDC, MDC ) + 1 IF ( ANYWND(ID) ) THEN THETA = SPCDIR(ID,1) 30.82 COSDIF = SPCDIR(ID,2)*CTW + SPCDIR(ID,3)*STW 40.41 SWINEB = TEMP1 * ( TEMP3 * COSDIF - 1.0 ) 40.41 SWINEB = MAX ( 0. , SWINEB * SIGMA ) ! IMATRA(ID,IS) = IMATRA(ID,IS) + SWINEB * AC2(ID,IS,KCGRD(1)) IF (TESTFL) PLWNDS(ID,IS,IPTST) = PLWNDS(ID,IS,IPTST) + & SWINEB*AC2(ID,IS,KCGRD(1)) 40.85 40.00 GENC0(ID,IS,1) = GENC0(ID,IS,1) + SWINEB*AC2(ID,IS,KCGRD(1)) 40.85 !ESMF IF (SAVE_SINBAC) SINBAC(ID,IS,KCGRD(1)) = !ESMF & SINBAC(ID,IS,KCGRD(1)) + SWINEB*AC2(ID,IS,KCGRD(1)) ! END IF ENDDO ENDDO ! ! *** test output *** ! IF (ITEST.GE. 80.AND.TESTFL) THEN 40.00 WRITE(PRTEST,6000) KCGRD(1), THETAW*180./PI 6000 FORMAT(' SWIND3: POINT THETAW :',I5,E12.4) WRITE(PRTEST,6100) TEMP1, FPM, UFRIC 6100 FORMAT(' SWIND3: TEMP1 FPM UFRC :',3E12.4, /, & ' IS ID1 ID2 Wind source term') DO IS = 1, MSC WRITE(PRTEST,6200) IS, IDCMIN(IS), IDCMAX(IS), & (PLWNDS(ID,IS,IPTST), ID=IDCMIN(IS), IDCMAX(IS)) 6200 FORMAT(3I4, 600e12.4) ENDDO WRITE(PRTEST,*) END IF ! RETURN ! end of subroutine SWIND3 END ! !**************************************************************** ! SUBROUTINE SWIND4 (IDWMIN ,IDWMAX , & SPCSIG ,WIND10 ,THETAW ,XIS , & DD ,KWAVE ,IMATRA ,GENC0 , & IDCMIN ,IDCMAX ,AC2 ,UFRIC , & PLWNDS ,ISSTOP ,ITER ,USTAR ,ZELEN , & SPCDIR ,ANYWND ,IT ,TAUWV ) ! !****************************************************************** ! USE SWCOMM2 40.41 USE SWCOMM3 40.41 USE SWCOMM4 40.41 USE OCPCOMM4 40.41 !ESMF USE M_GENARR, ONLY: SAVE_SINBAC, SINBAC ! IMPLICIT NONE ! ! ! --|-----------------------------------------------------------|-- ! | Delft University of Technology | ! | Faculty of Civil Engineering | ! | Environmental Fluid Mechanics Section | ! | P.O. Box 5048, 2600 GA Delft, The Netherlands | ! | | ! | Programmers: The SWAN team | ! --|-----------------------------------------------------------|-- ! ! ! SWAN (Simulating WAves Nearshore); a third generation wave model ! Copyright (C) 1993-2017 Delft University of Technology ! ! This program is free software; you can redistribute it and/or ! modify it under the terms of the GNU General Public License as ! published by the Free Software Foundation; either version 2 of ! the License, or (at your option) any later version. ! ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! A copy of the GNU General Public License is available at ! http://www.gnu.org/copyleft/gpl.html#SEC3 ! or by writing to the Free Software Foundation, Inc., ! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ! ! ! 0. Authors ! ! 30.72: IJsbrand Haagsma ! 30.82: IJsbrand Haagsma ! 40.02: IJsbrand Haagsma ! 40.31: Tim Campbell and John Cazes ! 40.41: Marcel Zijlema ! 40.61: Roop Lalbeharry ! 40.85: Marcel Zijlema ! ! 1. Updates ! ! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN ! 30.82, Oct. 98: Updated description of several variables ! 40.02, Oct. 00: References to CDRAGP and TAUWP removed ! 40.31, Jul. 03: correction calculation TAUDIR in test output ! 40.41, Aug. 04: COS(THETA-THETAW) replaced by sumrule to make it cheaper ! 40.41, Oct. 04: common blocks replaced by modules, include files removed ! 40.61, Nov. 06: improvements to WAM4 based on WAM4.5 ! 40.85, Aug. 08: store wind input for output purposes ! ! 2. Purpose ! ! Computation of the source term for the wind input for a ! third generation wind growth model: ! ! 1) Computation of the exponential input term based on a ! quasi linear theory developped by Janssen (1989, 1991). ! This formulation should be used in combination with the ! whitecapping dissipation source term according to ! Janssen (1991) and Mastenbroek et al. (1993) ! ! 3. Method ! ! The exponential term for the wind input used is taken ! from Janssen (1991): ! ! Sin (s,d) = B*E(s,d) ! e ! * ! / kwave U \ 2 ! B = max(0. , (beta rhoaw | --------- | cos( d - dw ) sigma)) ! \ sigma / ! ! The friction velocity is a fucntion of the roughness ! length Ze. A first gues for U* is given by Wu (1982), ! which is computed in subroutine WINDP2. ! ! 4. Argument variables ! ! i SPCDIR: (*,1); spectral directions (radians) 30.82 ! (*,2); cosine of spectral directions 30.82 ! (*,3); sine of spectral directions 30.82 ! (*,4); cosine^2 of spectral directions 30.82 ! (*,5); cosine*sine of spectral directions 30.82 ! (*,6); sine^2 of spectral directions 30.82 ! i SPCSIG: Relative frequencies in computational domain in sigma-space 30.72 ! REAL SPCDIR(MDC,6) 30.82 REAL SPCSIG(MSC) 30.72 ! ! IDWMIN Minimum counter for spectral wind direction ! IDWMAX Maximum counter for spectral wind direction ! IS Counter of relative frequency band ! ID Counter of directional distribution ! MSC Maximum counter of relative frequency ! MDC Maximum counter of directional distribution ! ! REALS: ! --------- ! DD Directional band width ! GRAV Gravitational acceleration ! SWINEO Coefficient stored in IMATRA ! THETA Spectral direction ! THETAW Mean direction of the relative wind vector ! WIND10 Velocity of the relative wind vector ! UFRIC Wind friction velocity ! ZALP Wave growth parameter used in WAM4 40.61 ! ! one and more dimensional arrays: ! --------------------------------- ! KWAVE 2D Wavenumber ! IMATRA 2D Coefficients of right hand side of matrix ! PWIND 1D Wind coefficients ! IDCMIN 1D Frequency dependent counter ! IDCMAX 1D Frequency dependent counter ! USTAR 2D Friction velocity at previous iteration level ! ZELEN 2D Roughness length at previous iteration level ! ANYWND 1D Determine if wind input is active for bin ! ! PWIND(9) 1D Rho air / rho water ! ------------ ! PWIND(14) 1D alfa (which is tuned at 0.01) ! PWIND(15) 1D Kappa ( 0.41) ! PWIND(16) 1D Rho air ! PWIND(17) 1D Rho water ! ! 7. Common blocks used ! ! ! 5. SUBROUTINES CALLING ! ! SOURCE ! ! 6. SUBROUTINES USED ! ! --- ! ! 7. ERROR MESSAGES ! ! --- ! ! 8. REMARKS ! ! --- ! ! 9. STRUCTURE ! ! ------------------------------------------------------------ ! For ! -------------------------------------------------------------- ! End of the subroutine SWIND4 ! -------------------------------------------------------------- ! ! 10. SOURCE ! ! !*********************************************************************** ! INTEGER IDWMAX ,IDWMIN ,IDDUM ,ID ,ISSTOP ,IS ! REAL THETA ,THETAW ,DD ,SWINEB ,WIND10 , & ZO ,ZE ,BETA1 ,BETA2 ,UFRIC ,UFRIC2 ,DS , & ZARG ,ZLOG1 ,ZLOG2 ,ZCN1 ,ZCN2 ,ZCN ,XIS , & SIGMA ,SIGMA1 ,SIGMA2 ,WAVEN ,WAVEN1 ,WAVEN2 ,TAUW , & TAUTOT ,TAUDIR ,COS1 ,COS2 ,CW1 ,RHOA ,RHOW , & RHOAW ,ALPHA ,XKAPPA ,F1 ,TAUWX ,TAUWY ,SE1 , & CTW ,STW ,COSDIF , 40.41 & SE2 ,SINWAV ,COSWAV REAL ZALP 40.61 ! REAL AC2(MDC,MSC,MCGRD) , 30.21 & IMATRA(MDC,MSC) , & KWAVE(MSC,MICMAX) , & PLWNDS(MDC,MSC,NPTST), & USTAR(MCGRD) , & ZELEN(MCGRD) , & TAUWV(MCGRD) REAL GENC0(MDC,MSC,MGENR) 40.85 ! INTEGER IDCMIN(MSC) , & IDCMAX(MSC) ! LOGICAL ANYWND(MDC) 40.41 ! INTEGER IENT , ITER, IT, J, II REAL ZTEN, RATIO, BETAMX, TXHFR, TYHFR, CW2, X1, X2, & ZARG1, ZARG2, GAMHF, SIGMAX, SIGHF1, SIGHF2, ZCNHF1, & ZCNHF2, AUX, COS3, ZAHF1, ZAHF2, FACHFR, FA, FB, FC, FD, FE, & FCEN, FF1, FF2, FF3, DCEN, TAUNEW, XFAC2, BETA, ZLOG ! SAVE IENT DATA IENT/0/ IF (LTRACE) CALL STRACE (IENT,'SWIND4') ! ! ! *** initialization *** ! ALPHA = PWIND(14) XKAPPA = PWIND(15) RHOA = PWIND(16) RHOW = PWIND(17) RHOAW = RHOA / RHOW ZTEN = 10. RATIO = 0.75 BETAMX = 1.2 F1 = BETAMX / XKAPPA**2 ZALP = 0.011 40.61 CTW = COS(THETAW) 40.41 STW = SIN(THETAW) 40.41 ! IF ( NSTATC.EQ.1 .AND. IT.EQ.1 ) THEN 40.41 40.00 ! ! *** nonstationary and first time step (the number of *** ! *** iterations however still can increase per time step *** ! ZO = ALPHA * UFRIC * UFRIC / GRAV ZE = ZO / SQRT( 1. - RATIO ) USTAR(KCGRD(1)) = UFRIC 30.21 ZELEN(KCGRD(1)) = ZE 30.21 ELSE IF ( NSTATC.EQ.0 .AND. ICOND.EQ.4 .AND. ITER .EQ. 1 ) THEN 40.41 40.00 ! ! *** non-first stationary computations and first iteration *** ! ZO = ALPHA * UFRIC * UFRIC / GRAV ZE = ZO / SQRT( 1. - RATIO ) USTAR(KCGRD(1)) = UFRIC 30.21 ZELEN(KCGRD(1)) = ZE 30.21 ELSE IF ( NSTATC.EQ.0 .AND. ICOND.NE.4 .AND. ITER .EQ. 2 ) THEN 40.41 40.00 ! ! *** first stationary computation (this subroutine is never *** ! *** excecuted anyway, this subroutine in entered after 1 *** ! *** iteration) and thus calculate ZO and ZE as a first *** ! *** prediction only and only in the second sweep *** ! ZO = ALPHA * UFRIC * UFRIC / GRAV ZE = ZO / SQRT( 1. - RATIO ) USTAR(KCGRD(1)) = UFRIC ZELEN(KCGRD(1)) = ZE ELSE ! ! *** calculate wave stress using the value of the *** ! *** velocity U* and roughness length Ze from the *** ! *** previous iteration *** ! UFRIC = USTAR(KCGRD(1)) ZE = ZELEN(KCGRD(1)) ! TAUW = 0. TAUWX = 0. TAUWY = 0. TXHFR = 0. TYHFR = 0. ! ! *** use old friction velocity to calculate wave stress *** ! UFRIC2 = UFRIC * UFRIC ! DO IS = 1, MSC-1 SIGMA1 = SPCSIG(IS) 30.72 SIGMA2 = SPCSIG(IS+1) 30.72 WAVEN1 = KWAVE(IS,1) WAVEN2 = KWAVE(IS+1,1) DS = SIGMA2 - SIGMA1 CW1 = SIGMA1 / WAVEN1 CW2 = SIGMA2 / WAVEN2 ZCN1 = ALOG ( GRAV * ZE / CW1**2 ) ZCN2 = ALOG ( GRAV * ZE / CW2**2 ) X1 = (UFRIC/CW1 + ZALP)**2 40.61 X2 = (UFRIC/CW2 + ZALP)**2 40.61 DO IDDUM = IDWMIN, IDWMAX ID = MOD ( IDDUM - 1 + MDC, MDC ) + 1 THETA = SPCDIR(ID,1) 30.82 COSDIF = SPCDIR(ID,2)*CTW + SPCDIR(ID,3)*STW 40.41 SINWAV = SPCDIR(ID,3) 40.41 COSWAV = SPCDIR(ID,2) 40.41 COS1 = MAX ( 0. , COSDIF ) 40.41 COS2 = COS1 * COS1 BETA1 = 0. BETA2 = 0. ! ! *** Miles constant Beta *** ! IF ( COS1 .GT. 0.01 ) THEN ZARG1 = XKAPPA / ( (UFRIC / CW1 + ZALP ) * COS1 ) 40.61 ZARG2 = XKAPPA / ( (UFRIC / CW2 + ZALP ) * COS1 ) 40.61 ZLOG1 = ZCN1 + ZARG1 ZLOG2 = ZCN2 + ZARG2 IF (ZLOG1.LT.0.) BETA1 = F1 * X1 * EXP (ZLOG1) * ZLOG1**4 40.61 IF (ZLOG2.LT.0.) BETA2 = F1 * X2 * EXP (ZLOG2) * ZLOG2**4 40.61 ENDIF ! ! *** calculate wave stress by integrating input source *** ! *** term in x- and y direction respectively *** ! SE1 = BETA1 * SIGMA1**3 * AC2(ID,IS ,KCGRD(1)) 40.61 SE2 = BETA2 * SIGMA2**3 * AC2(ID,IS+1,KCGRD(1)) 40.61 ! TAUWX = TAUWX + 0.5 * ( SE1 + SE2 ) * DS * COSWAV * COS2 TAUWY = TAUWY + 0.5 * ( SE1 + SE2 ) * DS * SINWAV * COS2 ! ! *** test output *** ! IF (ITEST.GE. 40 .AND. TESTFL) THEN WRITE(PRINTF,105) IS, ID, UFRIC, ZE 105 FORMAT(' SW4: IS ID UFRIC ZE :',2I4,2E12.4) WRITE(PRINTF,106) ZLOG1, ZLOG2, BETA1, BETA2 106 FORMAT(' SW4: ZOLG1-2 BETA1 BETA2:',4E12.4) IF (ABS(TAUWX).GT.0. .OR. ABS(TAUWY).GT.0.) THEN 40.31 TAUDIR = ATAN2 ( TAUWX, TAUWY ) 40.31 ELSE 40.31 TAUDIR = 0. 40.31 ENDIF 40.31 TAUDIR = MOD ( (TAUDIR + 2. * PI) , (2. * PI) ) 40.31 WRITE(PRINTF,107) TAUWX, TAUWY, TAUDIR*180./PI 107 FORMAT(' SW4: TAUWX TAUWY TAUDIR :',3E12.4) ENDIF ! ENDDO ENDDO ! ! *** determine effect of high frequency tail to wave stress *** ! *** assuming deep water conditions *** ! GAMHF = XKAPPA * GRAV / UFRIC SIGMAX = SPCSIG(MSC) 30.72 SIGHF1 = SIGMAX DO J=1, 50 SIGHF2 = XIS * SIGHF1 DS = SIGHF2 - SIGHF1 ZCNHF1 = ALOG ( ZE * SIGHF1**2 / GRAV ) ZCNHF2 = ALOG ( ZE * SIGHF2**2 / GRAV ) AUX = 0.0 SIGHF1 = SIGMAX 40.61 CW1 = GRAV/SIGHF1 40.61 CW2 = GRAV/SIGHF2 40.61 DO IDDUM = IDWMIN, IDWMAX ID = MOD ( IDDUM - 1 + MDC, MDC ) + 1 THETA = SPCDIR(ID,1) 30.82 COSDIF = SPCDIR(ID,2)*CTW + SPCDIR(ID,3)*STW 40.41 SINWAV = SPCDIR(ID,3) 40.41 COSWAV = SPCDIR(ID,2) 40.41 COS1 = MAX ( 0. , COSDIF ) 40.41 COS2 = COS1 * COS1 COS3 = COS2 * COS1 40.61 BETA1 = 0.0 BETA2 = 0.0 ! IF ( COS1 .GT. 0.01 ) THEN ! *** beta is independent of direction ! *** ZAHF1 = XKAPPA / ( UFRIC / CW1 + ZALP ) 40.61 ZAHF2 = XKAPPA / ( UFRIC / CW2 + ZALP ) 40.61 ZLOG1 = ZCNHF1 + ZAHF1 ZLOG2 = ZCNHF2 + ZAHF2 IF ( ZLOG1 .LT. 0. ) BETA1 = F1 * EXP (ZLOG1) * ZLOG1**4 IF ( ZLOG2 .LT. 0. ) BETA2 = F1 * EXP (ZLOG2) * ZLOG2**4 AUX = AUX + BETA1 + BETA2 ENDIF ! ! *** calculate contribution of high frequency tail to *** ! *** wave stress by integrating input source term in *** ! *** x- and y direction respectively *** ! FACHFR = SIGMAX**6 * AC2(ID,MSC,KCGRD(1)) * COS2 / GRAV**2 30.21 SE1 = FACHFR * BETA1 / SIGHF1 SE2 = FACHFR * BETA2 / SIGHF2 ! TXHFR = TXHFR + 0.5 * ( SE1 + SE2 ) * DS * COSWAV TYHFR = TYHFR + 0.5 * ( SE1 + SE2 ) * DS * SINWAV ! ! *** if coeffcient BETA = 0. for a frequency over *** ! *** all directions is zero skip loop *** ! IF ( AUX .EQ. 0. ) GOTO 5000 ! ENDDO ! IF (ITEST.GE. 45 ) THEN WRITE(PRINTF,407) XIS, SIGHF1, SIGHF2, J 407 FORMAT(' SW4: XIS SIGHF1 SIGHF2 J :',3E12.4,I4) WRITE(PRINTF,437) TXHFR, TYHFR, BETA1, BETA2 437 FORMAT(' SW4: TXHFR TYHFR BETA1,2 :',4E12.4) ENDIF ! SIGHF1 = SIGHF2 ENDDO 5000 CONTINUE ! IF ( ITEST .GE. 45 ) THEN WRITE(PRINTF,321) TAUWX, TAUWY, TXHFR, TYHFR 321 FORMAT(' SW4: Twx Twy Thfx Thfy:',4E12.4) ENDIF ! TAUTOT = RHOA * UFRIC2 ! *** wave stress *** TAUWX = TAUWX + TXHFR * UFRIC2 40.61 TAUWY = TAUWY + TYHFR * UFRIC2 40.61 IF (ABS(TAUWX).GT.0. .OR. ABS(TAUWY).GT.0.) THEN TAUDIR = ATAN2 ( TAUWX, TAUWY ) ELSE TAUDIR = 0. ENDIF TAUDIR = MOD ( (TAUDIR + 2. * PI) , (2. * PI) ) TAUW = RHOA * DD * SQRT ( TAUWX**2 + TAUWY**2 ) 40.61 TAUW = MIN ( TAUW , 0.999 * TAUTOT ) ! IF ( ITEST .GE. 45 ) THEN RATIO = TAUW / TAUTOT WRITE(PRINTF,301) TAUW, TAUTOT, RATIO, KCGRD(1) 30.21 301 FORMAT(' SW4: Tauw Taut ratio :',3E12.4,' in ',I5) ENDIF ! DO II = 1, 20 ! *** start iteration process *** FA = SQRT ( 1. - TAUW / TAUTOT ) FB = ZTEN * RHOA * GRAV / ALPHA FC = FA * ( FB / TAUTOT - 1. ) FD = SQRT ( TAUTOT ) FE = ALOG ( FC + 1. ) ! ! *** calculate function value and derivative in *** ! *** numerical point considered *** ! FCEN = FD * FE - SQRT(RHOA) * WIND10 * XKAPPA FF1 = 0.5 * FE / FD FF2 = 0.5 * TAUW * FC / FA - FA * FB FF3 = TAUTOT**1.5 * ( FC + FA ) DCEN = FF1 + FF2 / FF3 ! ! *** new total stress *** ! TAUNEW = TAUTOT - FCEN / DCEN ! IF ( ITEST .GT. 30 .AND. TESTFL ) THEN WRITE(PRINTF,440) TAUTOT, TAUNEW, FCEN, DCEN, II 440 FORMAT(' SW4: Tt Tnew Fcn DFcn II:',4E12.4,I2) WRITE(PRINTF,450) FA, FB, FC, FD 450 FORMAT(' SW4: FA FB FC FD :',4E12.4) WRITE(PRINTF,460) FE, FF1, FF2, FF3 460 FORMAT(' SW4: FE FF1 FF2 FF3 :',4E12.4) ENDIF ! IF ( TAUNEW .LE. TAUW ) TAUNEW = .5 * (TAUTOT + TAUW) 20.81 IF ( ABS ( TAUNEW - TAUTOT ) .LE. 1.E-5 ) GOTO 3000 ! TAUTOT = TAUNEW ENDDO 3000 CONTINUE ! UFRIC = SQRT ( TAUTOT / RHOA ) ! IF ( ITEST .GE. 20 .AND. TESTFL ) THEN WRITE(PRINTF,200) KCGRD(1) 200 FORMAT(' SW4: Values after Newton-Raphson in point:',I5) WRITE(PRINTF,206) TAUW, TAUTOT, TAUW/TAUTOT, UFRIC 206 FORMAT(' SW4: Tauw Taut rat Us :',4E12.4) WRITE(PRINTF,*) ENDIF ! ZO = ALPHA * UFRIC * UFRIC / GRAV ZE = ZO / SQRT ( 1. - TAUW / TAUTOT ) ! USTAR(KCGRD(1)) = UFRIC ZELEN(KCGRD(1)) = ZE TAUWV(KCGRD(1)) = TAUW ! ENDIF ! ! -----> ! ! *** calculate critical height and Miles parameter and *** ! *** calculate input source term B for with the updated *** ! *** values of UFRIC and ZE *** ! UFRIC2 = UFRIC * UFRIC ! DO IS = 1, ISSTOP SIGMA = SPCSIG(IS) 30.72 WAVEN = KWAVE(IS,1) CW1 = SIGMA / WAVEN ZCN = ALOG ( GRAV * ZE / CW1**2 ) DO IDDUM = IDCMIN(IS), IDCMAX(IS) ID = MOD ( IDDUM - 1 + MDC, MDC ) + 1 IF ( ANYWND(ID) ) THEN THETA = SPCDIR(ID,1) 30.82 COSDIF = SPCDIR(ID,2)*CTW + SPCDIR(ID,3)*STW 40.41 COS1 = MAX ( 0. , COSDIF ) 40.41 COS2 = COS1 * COS1 XFAC2 = ( UFRIC / CW1 + ZALP )**2 40.61 BETA = 0. IF ( COS1 .GT. 0.01 ) THEN ZARG = XKAPPA / ( (UFRIC / CW1 + ZALP ) * COS1 ) 40.61 ZLOG = ZCN + ZARG IF ( ZLOG .LT. 0. ) BETA = F1 * EXP (ZLOG) * ZLOG**4 ENDIF ! ! *** compute the factor B and store result in array *** ! SWINEB = RHOAW * BETA * XFAC2 * COS2 * SIGMA IMATRA(ID,IS) = IMATRA(ID,IS) + SWINEB * AC2(ID,IS,KCGRD(1)) 30.21 IF (TESTFL) PLWNDS(ID,IS,IPTST) = PLWNDS(ID,IS,IPTST) + & SWINEB*AC2(ID,IS,KCGRD(1)) 40.85 40.00 GENC0(ID,IS,1) = GENC0(ID,IS,1) + SWINEB*AC2(ID,IS,KCGRD(1)) 40.85 !ESMF IF (SAVE_SINBAC) SINBAC(ID,IS,KCGRD(1)) = !ESMF & SINBAC(ID,IS,KCGRD(1)) + SWINEB*AC2(ID,IS,KCGRD(1)) ! ! *** test output *** ! IF (ITEST.GE. 30 .AND. TESTFL) THEN WRITE(PRINTF,101) ZARG,ZLOG,ID,IS 101 FORMAT(' SW4: ZARG ZLOG ID IS :',2E12.4,' in',2I3) WRITE(PRINTF,102) COS1, COS2, XFAC2, BETA 102 FORMAT(' SW4: COS1/2 XFAC BETA :',4E12.4) ENDIF ENDIF ENDDO ! IF (ITEST.GE. 20 .AND. TESTFL) THEN WRITE(PRINTF,1102) SIGMA, WAVEN, CW1, ZCN 1102 FORMAT(' SW4: SIG WAV CW1 ZCN :',4E12.4) ENDIF ENDDO ! ! *** test output *** ! IF (ITEST.GE. 20.AND.TESTFL) THEN WRITE(PRINTF,9001) KCGRD(1), IDWMIN, IDWMAX 9001 FORMAT(' SW4: POINT IDWMIN IDWMAX :',3I5) WRITE(PRINTF,6053) WIND10,UFRIC,THETAW*180./PI 6053 FORMAT(' SW4: WIND10 UFRIC THETAW :',3E12.4) WRITE(PRINTF,7136) PWIND(9), PWIND(16), PWIND(17) 7136 FORMAT(' SW4: RHOAW RHOA RHOW :',3E12.4) WRITE(PRINTF,7126) PWIND(14), PWIND(15), ZTEN 7126 FORMAT(' SW4: ALPHA XKAPPA ZTEN :',3E12.4) END IF ! RETURN ! end of subroutine SWIND4 END ! !**************************************************************** ! SUBROUTINE SWIND5 (SPCSIG ,THETAW ,ISSTOP , & UFRIC ,KWAVE ,IMATRA ,IDCMIN , & IDCMAX ,AC2 ,ANYWND ,PLWNDS , & SPCDIR ,GENC0 ) ! !**************************************************************** ! USE SWCOMM3 40.41 USE SWCOMM4 40.41 USE OCPCOMM4 40.41 !ESMF USE M_GENARR, ONLY: SAVE_SINBAC, SINBAC ! IMPLICIT NONE ! ! ! --|-----------------------------------------------------------|-- ! | Delft University of Technology | ! | Faculty of Civil Engineering | ! | Environmental Fluid Mechanics Section | ! | P.O. Box 5048, 2600 GA Delft, The Netherlands | ! | | ! | Programmers: The SWAN team | ! --|-----------------------------------------------------------|-- ! ! ! SWAN (Simulating WAves Nearshore); a third generation wave model ! Copyright (C) 1993-2017 Delft University of Technology ! ! This program is free software; you can redistribute it and/or ! modify it under the terms of the GNU General Public License as ! published by the Free Software Foundation; either version 2 of ! the License, or (at your option) any later version. ! ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! A copy of the GNU General Public License is available at ! http://www.gnu.org/copyleft/gpl.html#SEC3 ! or by writing to the Free Software Foundation, Inc., ! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ! ! ! 0. Authors ! ! 30.72: IJsbrand Haagsma ! 30.82: IJsbrand Haagsma ! 40.41: Marcel Zijlema ! 40.53: Andre van der Westhuysen ! 40.85: Marcel Zijlema ! ! 1. Updates ! ! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN ! 30.82, Oct. 98: Updated description of several variables ! 40.41, Aug. 04: COS(THETA-THETAW) replaced by sumrule to make it cheaper ! 40.41, Oct. 04: common blocks replaced by modules, include files removed ! 40.53, Aug. 04: changes parameters of Yan formulae in case of Alves and ! Banner whitecapping method ! 40.85, Aug. 08: store wind input for output purposes ! ! 2. Purpose ! ! Computation of the source term for the wind input for a ! third generation wind growth model: ! ! The exponential input term is according to Yan (1987). This ! input term is valid for the higher frequency part of the ! spectrum (strongly forced wave components). The expression ! reduces to the Snyder (1982) expression form for spectral ! wave components with weak wind forcing and to the Plant (1982) ! form for more strongly forced wave components. ! ! 3. Method ! ! The expression reads --> with X = Ustar / C ! ! / / 2 \ ! Sin = | | 0.04 X + 0.00544 X + 0.000055 | * cos (theta) ! \ \ / ! \ ! - 0.00031 | sigma * AC2(d,s,x,y) ! / ! ! 4. Argument variables ! ! i SPCDIR: (*,1); spectral directions (radians) 30.82 ! (*,2); cosine of spectral directions 30.82 ! (*,3); sine of spectral directions 30.82 ! (*,4); cosine^2 of spectral directions 30.82 ! (*,5); cosine*sine of spectral directions 30.82 ! (*,6); sine^2 of spectral directions 30.82 ! i SPCSIG: Relative frequencies in computational domain in sigma-space 30.72 ! REAL SPCDIR(MDC,6) 30.82 REAL SPCSIG(MSC) 30.72 ! ! IS Counter of relative frequency band ! ID Counter of directional distribution ! MSC Maximum counter of relative frequency ! MDC Maximum counter of directional distribution ! ! REALS: ! --------- ! FPM Pierson Moskowitz frequecy (WAM) ! THETA Spectral direction ! THETAW Mean direction of the relative wind vector ! UFRIC Wind friction velocity ! ! one and more dimensional arrays: ! --------------------------------- ! KWAVE 2D Wavenumber ! IMATRA 2D Coefficients of right hand side of matrix ! IDCMIN 1D Frequency dependent counter ! IDCMAX 1D Frequency dependent counter ! ANYWND 1D Wind input for bin considered ! ! 7. Common blocks used ! ! ! 5. SUBROUTINES CALLING ! ! SOURCE ! ! 6. SUBROUTINES USED ! ! --- ! ! 7. ERROR MESSAGES ! ! --- ! ! 8. REMARKS ! ! --- ! ! 9. STRUCTURE ! ! ------------------------------------------------------------ ! For every spectral bin that fall within a sweep considered ! compute the source term and store the results in IMATRA ! -------------------------------------------------------------- ! End of the subroutine SWIND5 ! -------------------------------------------------------------- ! ! 10. SOURCE ! !*********************************************************************** ! INTEGER IENT, IDDUM ,ID ,IS ,ISSTOP ! REAL UFRIC ,THETA ,THETAW ,SIGMA ,SWINEB , TEMP3, & CTW ,STW ,COSDIF , 40.41 & USTAC1 ,USTAC2 ,COF1 ,COF2 ,COF3 ,COF4 ! REAL AC2(MDC,MSC,MCGRD) , & IMATRA(MDC,MSC) , & KWAVE(MSC,MICMAX) , & PLWNDS(MDC,MSC,NPTST) REAL GENC0(MDC,MSC,MGENR) 40.85 ! INTEGER IDCMIN(MSC) , & IDCMAX(MSC) ! LOGICAL ANYWND(MDC) ! SAVE IENT DATA IENT/0/ IF (LTRACE) CALL STRACE (IENT,'SWIND5') ! ! *** input according to Yan (1987) *** ! COF1 = 0.04 COF2 = 0.00544 COF3 = 0.000055 COF4 = 0.00031 ! ! adapted Yan fit for use of Alves and Banner method 40.53 ! IF (IWCAP.EQ.7) THEN 40.53 COF1 = 0.04 40.53 COF2 = 0.00552 40.53 COF3 = 0.000052 40.53 COF4 = 0.000302 40.53 END IF 40.53 ! CTW = COS(THETAW) 40.41 STW = SIN(THETAW) 40.41 DO IS = 1, ISSTOP SIGMA = SPCSIG(IS) USTAC1 = ( UFRIC * KWAVE(IS,1) ) / SIGMA USTAC2 = USTAC1 * USTAC1 TEMP3 = ( COF1 * USTAC2 + COF2 * USTAC1 + COF3) DO IDDUM = IDCMIN(IS), IDCMAX(IS) ID = MOD ( IDDUM - 1 + MDC, MDC ) + 1 IF ( ANYWND(ID) ) THEN THETA = SPCDIR(ID,1) 30.82 COSDIF = SPCDIR(ID,2)*CTW + SPCDIR(ID,3)*STW 40.41 SWINEB = TEMP3 * COSDIF - COF4 40.41 SWINEB = MAX ( 0. , SWINEB * SIGMA ) IMATRA(ID,IS) = IMATRA(ID,IS) + SWINEB * AC2(ID,IS,KCGRD(1)) IF (TESTFL) PLWNDS(ID,IS,IPTST) = PLWNDS(ID,IS,IPTST) + & SWINEB*AC2(ID,IS,KCGRD(1)) 40.85 40.00 GENC0(ID,IS,1) = GENC0(ID,IS,1) + SWINEB*AC2(ID,IS,KCGRD(1)) 40.85 !ESMF IF (SAVE_SINBAC) SINBAC(ID,IS,KCGRD(1)) = !ESMF & SINBAC(ID,IS,KCGRD(1)) + SWINEB*AC2(ID,IS,KCGRD(1)) END IF ENDDO ENDDO ! ! *** test output *** ! IF (ITEST.GE. 60.AND.TESTFL) THEN WRITE(PRINTF,6000) KCGRD(1), THETAW*180./PI, UFRIC 6000 FORMAT(' SWIND5: POINT THETAW UFRIC :',I5,2E12.4) WRITE(PRINTF,6100) COF1, COF2, COF3, COF4 6100 FORMAT(' SWIND5: COF1 COF2 COF3 COF4 :',4E12.4) WRITE(PRINTF,*) END IF ! RETURN ! end of subroutine SWIND5 END !