#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
!