swancom3.f90 File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	windp1 (WIND10, THETAW, IDWMIN, IDWMAX, FPM, UFRIC, WX2, WY2, SPCDIR, UX2, UY2, SPCSIG, IG)
subroutine	windp2 (IDWMIN, IDWMAX, SIGPKD, FPM, ETOTW, AC2, SPCSIG, WIND10)
subroutine	windp3 (ISSTOP, ALIMW, AC2, GROWW, IDCMIN, IDCMAX)
subroutine	swind0 (IDCMIN, IDCMAX, ISSTOP, SPCSIG, THETAW, UFRIC, FPM, PLWNDA, IMATRA, SPCDIR)
subroutine	swind3 (SPCSIG, THETAW, IMATDA, KWAVE, IMATRA, IDCMIN, IDCMAX, UFRIC, FPM, PLWNDB, ISSTOP, SPCDIR, IG)
subroutine	swind4 (IDWMIN, IDWMAX, SPCSIG, WIND10, THETAW, XIS, DD, KWAVE, IMATRA, IMATDA, IDCMIN, IDCMAX, UFRIC, PLWNDB, ISSTOP, ITER, USTAR, ZELEN, SPCDIR, IT, IG)
subroutine	swind5 (SPCSIG, THETAW, ISSTOP, UFRIC, KWAVE, IMATRA, IMATDA, IDCMIN, IDCMAX, PLWNDB, SPCDIR, IG)
subroutine	wndpar (ISSTOP, IDWMIN, IDWMAX, IDCMIN, IDCMAX, DEP2, WIND10, THETAW, AC2, KWAVE, IMATRA, IMATDA, SPCSIG, CGO, ALIMW, GROWW, ETOTW, PLWNDA, PLWNDB, SPCDIR, ITER)

Function/Subroutine Documentation

swind0()

subroutine swind0	(	integer, dimension(msc)	IDCMIN,
		integer, dimension(msc)	IDCMAX,
		integer	ISSTOP,
		real, dimension(msc)	SPCSIG,
		real	THETAW,
		real	UFRIC,
		real	FPM,
		real, dimension(mdc,msc,nptst)	PLWNDA,
		real, dimension(mdc,msc)	IMATRA,
		real, dimension(mdc,6)	SPCDIR
	)

Definition at line 471 of file swancom3.f90.

!
!****************************************************************
!
!     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)
!
!  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
!
!****************************************************************
!
   USE swcomm3                                                         
   USE swcomm4                                                         
   USE ocpcomm4                                                        
   
   IMPLICIT NONE

   REAL    :: SPCDIR(MDC,6) ,SPCSIG(MSC)                                                 
   INTEGER :: IDDUM   ,ID      ,IS      ,ISSTOP
   REAL    :: FPM     ,UFRIC   ,THETA   ,THETAW  ,SWINEA  ,SIGMA   ,  &
              TEMP1   ,TEMP2   ,CTW     ,STW     ,COSDIF  ,TEMP3   ,  &
          FILTER  ,ARGU
   REAL    :: IMATRA(MDC,MSC)      ,PLWNDA(MDC,MSC,NPTST)        
   INTEGER :: IDCMIN(MSC)          ,IDCMAX(MSC)
!
!     *** calculate linear wind input term ***
!
   ctw = cos(thetaw)                                                   
   stw = sin(thetaw)                                                   
   fpm =  grav_w / ( 28.0 * ufric )
   temp1 = pwind(31) / ( grav_w**2 * 2. * pi_w ) 
   DO is = 1, isstop
     sigma  = spcsig(is)                                               
!
!    ****            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)                                    
     filter = exp( - argu**4 )                                        
     temp2  = temp1 / sigma
     DO iddum = idcmin(is), idcmax(is)
       id = mod( iddum - 1 + mdc, mdc ) + 1
       IF(sigma >= (0.7 * fpm))THEN
         theta  = spcdir(id,1)                                         
         cosdif = spcdir(id,2)*ctw + spcdir(id,3)*stw                  
         temp3  = ( ufric *  max( 0. , cosdif))**4                     
         swinea = max( 0. , temp2 * temp3 * filter )
         imatra(id,is) = imatra(id,is) + swinea
         IF(testfl) plwnda(id,is,iptst) = swinea                       
       END IF
     ENDDO
   ENDDO

   RETURN

swind3()

subroutine swind3	(	real, dimension(msc)	SPCSIG,
		real	THETAW,
		real, dimension(mdc,msc)	IMATDA,
		real, dimension(msc,icmax)	KWAVE,
		real, dimension(mdc,msc)	IMATRA,
		integer, dimension(msc)	IDCMIN,
		integer, dimension(msc)	IDCMAX,
		real	UFRIC,
		real	FPM,
		real, dimension(mdc,msc,nptst)	PLWNDB,
		integer	ISSTOP,
		real, dimension(mdc,6)	SPCDIR,
		integer	IG
	)

Definition at line 556 of file swancom3.f90.

!
!****************************************************************
!
!     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 combinated with the dissipation
!         term of Komen et al. (1984)
!
!  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
!
!
!****************************************************************
   USE swcomm3                                                         
   USE swcomm4                                                         
   USE ocpcomm4   
!   USE ALL_VARS, ONLY : MT,AC2  
   USE vars_wave, ONLY : mt,ac2  
   
   IMPLICIT NONE                                                   

   REAL    :: SPCDIR(MDC,6),SPCSIG(MSC)                                                 
   INTEGER :: IDDUM ,ID    ,IS    ,ISSTOP  ,IG
   REAL    :: FPM   ,UFRIC ,THETA ,THETAW,SIGMA ,SWINEB,TEMP1,          &
              CTW   ,STW   ,COSDIF,TEMP2 ,TEMP3 ,CINV
   REAL    :: IMATDA(MDC,MSC),IMATRA(MDC,MSC),        &
              KWAVE(MSC,ICMAX),PLWNDB(MDC,MSC,NPTST)
   INTEGER :: IDCMIN(MSC),IDCMAX(MSC)

   ctw   = cos(thetaw)                                                 
   stw   = sin(thetaw)                                                 
   temp1 = 0.25 * pwind(9)
   temp2 = 28.0 * ufric
   DO is = 1, isstop
     sigma = spcsig(is)                                                
     cinv  = kwave(is,1) / sigma
     temp3 = temp2 * cinv
     DO iddum = idcmin(is), idcmax(is)
       id = mod( iddum - 1 + mdc, mdc ) + 1
       theta  = spcdir(id,1)                                         
       cosdif = spcdir(id,2)*ctw + spcdir(id,3)*stw                  
       swineb = temp1 * ( temp3 * cosdif - 1.0 )   
       swineb = max( 0. , swineb * sigma )

       imatra(id,is) = imatra(id,is) + swineb * ac2(id,is,ig)
       imatda(id,is) = imatda(id,is) - swineb
       IF(testfl) plwndb(id,is,iptst) = swineb                      
     ENDDO
   ENDDO

   RETURN

swind4()

subroutine swind4	(	integer	IDWMIN,
		integer	IDWMAX,
		real, dimension(msc)	SPCSIG,
		real	WIND10,
		real	THETAW,
		real	XIS,
		real	DD,
		real, dimension(msc,icmax)	KWAVE,
		real, dimension(mdc,msc)	IMATRA,
		real, dimension(mdc,msc)	IMATDA,
		integer, dimension(msc)	IDCMIN,
		integer, dimension(msc)	IDCMAX,
		real	UFRIC,
		real, dimension(mdc,msc,nptst)	PLWNDB,
		integer	ISSTOP,
		integer	ITER,
		real, dimension(mt)	USTAR,
		real, dimension(mt)	ZELEN,
		real, dimension(mdc,6)	SPCDIR,
		integer	IT,
		integer	IG
	)

Definition at line 644 of file swancom3.f90.

!
!******************************************************************
!
!     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)
!
!  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.
!
!******************************************************************
   USE swcomm2                                                         
   USE swcomm3                                                         
   USE swcomm4                                                         
   USE ocpcomm4  
!   USE ALL_VARS, ONLY : MT,AC2    
   USE vars_wave, ONLY : mt,ac2    
   
   IMPLICIT NONE                                                  
   REAL    :: SPCDIR(MDC,6),SPCSIG(MSC)                                                 
   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 ,                                   &
          SE2    ,SINWAV ,COSWAV
   REAL    :: IMATDA(MDC,MSC)      ,                                     &
              IMATRA(MDC,MSC)      ,                                      &
          KWAVE(MSC,ICMAX)     ,                                      &
          PLWNDB(MDC,MSC,NPTST),                                      &
          USTAR(MT)         ,                                      &
          ZELEN(MT)
   INTEGER :: IDCMIN(MSC)          ,IDCMAX(MSC)
   REAL    :: ZTEN ,RATIO ,BETAMX ,TXHFR ,TYHFR ,CW2 ,ZARG1 ,ZARG2
   REAL    :: GAMHF ,SIGMAX ,SIGHF1 ,SIGHF2 ,ZCNHF1 ,ZCNHF2 ,AUX
   REAL    :: ZAHF1 ,ZAHF2 ,FACHFR ,FA ,FB ,FC ,FD ,FE ,FCEN
   REAL    :: FF1 ,FF2 ,FF3 ,DCEN ,TAUNEW ,XFAC2 ,BETA ,ZLOG
   INTEGER :: IT ,IG ,ITER ,J ,II                                         

!
!  *** 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
   ctw    = cos(thetaw)                                                
   stw    = sin(thetaw)                                                
!
   IF(nstatc == 1 .AND. it == 1)THEN                               
!
!     *** nonstationary and first time step (the number of        ***
!     *** iterations however still can increase per time step     ***
!
     zo     = alpha * ufric * ufric / grav_w
     ze     = zo / sqrt( 1. - ratio )
     ustar(ig) = ufric                                           
     zelen(ig) = ze                                              
   ELSE IF(nstatc == 0 .AND. icond == 4 .AND. iter == 1)THEN     
!
!    *** non-first stationary computations and first iteration   ***
!
     zo     = alpha * ufric * ufric / grav_w
     ze     = zo / sqrt( 1. - ratio )
     ustar(ig) = ufric                                           
     zelen(ig) = ze                                              
   ELSE IF ( nstatc.EQ.0 .AND. icond.NE.4 .AND. iter .EQ. 2 ) THEN     
!
!    *** 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_w
     ze     = zo / sqrt( 1. - ratio )
     ustar(ig) = ufric
     zelen(ig) = ze
   ELSE
!
!    *** calculate wave stress using the value of the  ***
!    *** velocity U* and roughness length Ze from the  ***
!    *** previous iteration                            ***
!
     ufric = ustar(ig)
     ze    = zelen(ig)
!
     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)                                             
       sigma2 = spcsig(is+1)                                           
       waven1 = kwave(is,1)
       waven2 = kwave(is+1,1)
       ds     = sigma2 - sigma1
       cw1    = sigma1 / waven1
       cw2    = sigma2 / waven2
       zcn1   = alog( grav_w * ze / cw1**2 )
       zcn2   = alog( grav_w * ze / cw2**2 )
       DO iddum = idwmin, idwmax
         id = mod( iddum - 1 + mdc, mdc ) + 1
         theta  = spcdir(id,1)                                         
         cosdif = spcdir(id,2)*ctw + spcdir(id,3)*stw                  
         sinwav = spcdir(id,3)                                         
         coswav = spcdir(id,2)                                         
         cos1   = max( 0. , cosdif )                                  
         cos2   = cos1 * cos1
         beta1  = 0.
         beta2  = 0.
!
!        *** Miles constant Beta ***
!
         IF(cos1 > 0.01)THEN
           zarg1 = xkappa * cw1 / ( ufric * cos1 )
           zarg2 = xkappa * cw2 / ( ufric * cos1 )
           zlog1 = zcn1 + zarg1
           zlog2 = zcn2 + zarg2
           IF(zlog1 < 0.) beta1 = f1 * exp(zlog1) * zlog1**4
           IF(zlog2 < 0.) beta2 = f1 * exp(zlog2) * zlog2**4
         ENDIF
!
!        *** calculate wave stress by integrating input source ***
!        *** term in x- and y direction respectively           ***
!
         se1 = waven1**2 * beta1 * sigma1 * ac2(id,is  ,ig)
         se2 = waven2**2 * beta2 * sigma2 * ac2(id,is+1,ig)
!
         tauwx = tauwx + 0.5 * ( se1 + se2 ) * ds * coswav * cos2
         tauwy = tauwy + 0.5 * ( se1 + se2 ) * ds * sinwav * cos2
       ENDDO
     ENDDO
!
!    *** determine effect of high frequency tail to wave stress ***
!    *** assuming deep water conditions                         ***
!
     gamhf =  xkappa * grav_w / ufric
     sigmax = spcsig(msc)                                              
     sighf1 = sigmax
     DO j=1, 50
       sighf2 = xis * sighf1
       ds     = sighf2 - sighf1
       zcnhf1 = alog( ze * sighf1**2 / grav_w )
       zcnhf2 = alog( ze * sighf2**2 / grav_w )
       aux    = 0.0
       DO iddum = idwmin, idwmax
         id = mod( iddum - 1 + mdc, mdc ) + 1
         theta  = spcdir(id,1)                                         
         cosdif = spcdir(id,2)*ctw + spcdir(id,3)*stw                  
         sinwav = spcdir(id,3)                                         
         coswav = spcdir(id,2)                                         
         cos1   = max( 0. , cosdif )                                  
         cos2   = cos1 * cos1
         beta1  = 0.0
         beta2  = 0.0
!
         IF(cos1 > 0.01)THEN
!          *** beta is independent of direction ! ***
           zahf1 = gamhf / sighf1
           zahf2 = gamhf / sighf2
           zlog1 = zcnhf1 + zahf1
           zlog2 = zcnhf2 + zahf2
           IF(zlog1 < 0.) beta1 = f1 * exp(zlog1) * zlog1**4
           IF(zlog2 < 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,ig) * cos2 / grav_w**2    

         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 == 0.) GOTO 5000
!
       ENDDO

       sighf1 = sighf2
     ENDDO
5000 CONTINUE

     tautot = rhoa * ufric2
!    *** wave stress ***
     tauwx  = tauwx + txhfr
     tauwy  = tauwy + tyhfr
     IF(abs(tauwx) > 0. .OR. abs(tauwy) > 0.)THEN
       taudir = atan2( tauwx, tauwy )
     ELSE
       taudir = 0.
     ENDIF
     taudir = mod( (taudir + 2. * pi_w) , (2. * pi_w) )
     tauw   = rhoa * ufric2 * dd * sqrt( tauwx**2 + tauwy**2 )
     tauw   = min( tauw , 0.999 * tautot )

     DO ii = 1, 20
!      *** start iteration process ***
       fa = sqrt( 1. - tauw / tautot )
       fb = zten * rhoa * grav_w / 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**2 - fa * fb
       ff3  = tautot**1.5 * ( fc + 1. )
       dcen = ff1 + ff2 / ff3
!
!      *** new total stress ***
!
       taunew = tautot - fcen / dcen
!
       IF(taunew <= tauw) taunew = .5 * (tautot + tauw)           
       IF(abs( taunew - tautot ) <= 1.e-5 ) GOTO 3000
!
       tautot = taunew
     ENDDO
3000 CONTINUE
!
     ufric  = sqrt( tautot / rhoa )
!
     zo     = alpha * ufric * ufric / grav_w
     ze     = zo / sqrt( 1. - tauw / tautot )
!
     ustar(ig) = ufric
     zelen(ig) = ze
!
   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)                                               
     waven  = kwave(is,1)
     cw1    = sigma / waven
     zcn    = alog( grav_w * ze / cw1**2 )
     DO iddum = idcmin(is), idcmax(is)
       id = mod( iddum - 1 + mdc, mdc ) + 1
         theta  = spcdir(id,1)                                         
         cosdif = spcdir(id,2)*ctw + spcdir(id,3)*stw                  
         cos1   = max( 0. , cosdif )                                  
         cos2   = cos1 * cos1
         xfac2  = ( ufric / cw1 )**2
         beta   = 0.
         IF(cos1 > 0.01)THEN
           zarg = xkappa * cw1 / ( ufric * cos1 )
           zlog = zcn + zarg
           IF(zlog < 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,ig)  
         imatda(id,is) = imatda(id,is) - swineb  
         IF(testfl) plwndb(id,is,iptst) = swineb                      
     ENDDO
   ENDDO

   RETURN

swind5()

subroutine swind5	(	real, dimension(msc)	SPCSIG,
		real	THETAW,
		integer	ISSTOP,
		real	UFRIC,
		real, dimension(msc,icmax)	KWAVE,
		real, dimension(mdc,msc)	IMATRA,
		real, dimension(mdc,msc)	IMATDA,
		integer, dimension(msc)	IDCMIN,
		integer, dimension(msc)	IDCMAX,
		real, dimension(mdc,msc,nptst)	PLWNDB,
		real, dimension(mdc,6)	SPCDIR,
		integer	IG
	)

Definition at line 962 of file swancom3.f90.

!
!****************************************************************
!
!     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:
!
!  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)
!                       /
!****************************************************************
!
   USE swcomm3                                                         
   USE swcomm4                                                         
   USE ocpcomm4 
!   USE ALL_VARS, ONLY : MT,AC2   
   USE vars_wave, ONLY : mt,ac2   
   
   IMPLICIT NONE                                                    

   REAL    :: SPCDIR(MDC,6),SPCSIG(MSC)                                                 
   INTEGER :: IDDUM  ,ID     ,IS     ,ISSTOP,  IG
   REAL    :: UFRIC  ,THETA  ,THETAW ,SIGMA  ,SWINEB ,                &
              CTW    ,STW    ,COSDIF ,                                &
          USTAC1 ,USTAC2 ,COF1   ,COF2   ,COF3   ,COF4
   REAL    :: IMATDA(MDC,MSC)   ,IMATRA(MDC,MSC)      , &
              KWAVE(MSC,ICMAX)  ,PLWNDB(MDC,MSC,NPTST)
   INTEGER :: IDCMIN(MSC)       ,IDCMAX(MSC)
   REAL    :: TEMP3

!
!  *** 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                  
!
   IF(iwcap == 7)THEN                                                
     cof1 = 0.04                                                      
     cof2 = 0.00552                                                   
     cof3 = 0.000052                                                  
     cof4 = 0.000302                                                  
   END IF                                                              
!
   ctw  = cos(thetaw)                                                  
   stw  = sin(thetaw)                                                  
   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
         theta  = spcdir(id,1)                                         
         cosdif = spcdir(id,2)*ctw + spcdir(id,3)*stw                  
         swineb = temp3 * cosdif - cof4                                
         swineb = max( 0. , swineb * sigma )
!         IMATRA(ID,IS) = IMATRA(ID,IS) + SWINEB * AC2(ID,IS,IGC)
         imatra(id,is) = imatra(id,is) + swineb * ac2(id,is,ig)
         imatda(id,is) = imatda(id,is) - swineb
         IF(testfl) plwndb(id,is,iptst) = swineb                      
     ENDDO
   ENDDO

   RETURN

windp1()

subroutine windp1	(	real	WIND10,
		real	THETAW,
		integer	IDWMIN,
		integer	IDWMAX,
		real	FPM,
		real	UFRIC,
		real, dimension(mt)	WX2,
		real, dimension(mt)	WY2,
		real, dimension(mdc,6)	SPCDIR,
		real, dimension(mt)	UX2,
		real, dimension(mt)	UY2,
		real, dimension(msc)	SPCSIG,
		integer	IG
	)

Definition at line 20 of file swancom3.f90.

!
!****************************************************************
!
!        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
!
!     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
!        -----------------------------------------------------------
!****************************************************************
!
   USE ocpcomm1                                                        
   USE ocpcomm2                                                        
   USE ocpcomm3                                                        
   USE ocpcomm4                                                        
   USE swcomm1                                                         
   USE swcomm2                                                         
   USE swcomm3                                                         
   USE swcomm4  
   USE all_vars, ONLY : mt                                                       
   
   IMPLICIT NONE
!
   REAL    :: SPCDIR(MDC,6),SPCSIG(MSC)        
   INTEGER :: IDWMIN ,IDWMAX
   INTEGER :: ID,IDDUM,IG                                         
!
   REAL    :: WIND10 ,THETAW ,UFRIC ,FPM ,CDRAG ,SDMEAN       
   REAL    :: WX2(MT), WY2(MT), UX2(MT), UY2(MT)         

   REAL         AWX, AWY, RWX, RWY                     
!
!  compute absolute wind velocity                                    
   IF(varwi)THEN
     awx = wx2(ig)
     awy = wy2(ig)
   ELSE
     awx = u10 * cos(wdic)
     awy = u10 * sin(wdic)
   END IF
!  compute relative wind velocity                                      
   IF(icur == 0)THEN
     rwx = awx
     rwy = awy
   ELSE
     rwx = awx - ux2(ig)
     rwy = awy - uy2(ig)
   END IF
!  compute absolute value of relative wind velocity                   
   wind10 = sqrt(rwx**2+rwy**2)
   IF(wind10 > 0.)THEN
     thetaw = atan2(rwy,rwx)
     thetaw = mod( (thetaw + pi2_w) , pi2_w )
   ELSE
     thetaw = 0.
   END IF
!
!  *** 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             
   sdmean = 0.5 * (spcdir(1,1) + spcdir(mdc,1))                   
   IF(thetaw < sdmean - pi_w) thetaw = thetaw + 2.*pi_w
   IF(thetaw > sdmean + pi_w) thetaw = thetaw - 2.*pi_w
!
   IF((thetaw-0.5*pi_w) <= spcdir(1,1))THEN                   
     IF((thetaw+1.5*pi_w) >= spcdir(mdc,1))THEN
       idwmin = 1
     ELSE
       idwmin = nint( (thetaw + 1.5*pi_w - spcdir(1,1)) / ddir ) + 1    
     END IF
   ELSE
     idwmin = nint( (thetaw - 0.5*pi_w - spcdir(1,1)) / ddir ) + 1      
   END IF
!
   IF((thetaw + 0.5 * pi_w) >= spcdir(mdc,1))THEN                  
     IF((thetaw - 1.5 * pi_w) <= spcdir(1,1))THEN                  
       idwmax = mdc
     ELSE
       idwmax = nint( (thetaw - 1.5 * pi_w - spcdir(1,1)) / ddir ) + 1  
     END IF
   ELSE
     idwmax = nint( (thetaw + 0.5 * pi_w - spcdir(1,1)) / ddir ) + 1    
   END IF
!
   IF(idwmin > idwmax) idwmax = mdc + idwmax
!
!  *** compute the Pierson Moskowitz frequency ***
!
   IF(iwind == 1 .OR. iwind == 2)THEN
!
!    *** first and second generation wind wave model ***
!
     IF(wind10 < pwind(12) ) wind10 = pwind(12)
     fpm = 2. * pi_w * pwind(13) * grav_w / wind10
!
!    *** determine U friction in case predictor is obtained ***
!    *** with second genaration wave model                  ***
!
     IF(wind10 > 7.5)THEN
       ufric = wind10 * sqrt(( 0.8 + 0.065 * wind10 ) * 0.001 )
     ELSE
       cdrag = 0.0012873
       ufric = sqrt( cdrag ) * wind10
     END IF
!
!    Reformulation of the wind speed in terms of friction velocity.      
!    This formulation is based on Bouws (1986) and described in Delft    
!    Hydraulics report H3515 (1999)                                      
!
     wind10 = wind10 * sqrt(((0.8 + 0.065 * wind10) * 0.001) /    &
                            ((0.8 + 0.065 * 15.   ) * 0.001))          
!
   ELSE IF(iwind >= 3)THEN
!
!    *** Calculate the wind friction velocity  ***
!    *** according to Wu (1982)                ***
!
     IF(wind10 > 7.5)THEN
       ufric = wind10 * sqrt(( 0.8 + 0.065 * wind10 ) * 0.001 )
     ELSE
       cdrag = 0.0012873
       ufric = sqrt( cdrag ) * wind10
     END IF
!
!    *** Wind friction velocity and PM-frequency ***
!
     ufric = max( 1.e-15 , ufric)
     fpm =  grav_w / ( 28.0 * ufric )

   END IF

   RETURN

windp2()

subroutine windp2	(	integer	IDWMIN,
		integer	IDWMAX,
			SIGPKD,
		real	FPM,
		real	ETOTW,
		real, dimension(mdc,msc,0:mt)	AC2,
		real, dimension(msc)	SPCSIG,
			WIND10
	)

Definition at line 210 of file swancom3.f90.

!  (This subroutine has not been tested yet)
!
!****************************************************************
!
   USE ocpcomm1                                                        
   USE ocpcomm2                                                        
   USE ocpcomm3                                                        
   USE ocpcomm4                                                        
   USE swcomm1                                                         
   USE swcomm2                                                         
   USE swcomm3                                                         
   USE swcomm4 
   USE all_vars, ONLY : mt                                                        
!
!
!  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"
!
!
   REAL    SPCSIG(MSC)                                                 
!
!     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
!     ------------------------------------------------------------
!
!***********************************************************************
!
!
   INTEGER  IDWMIN  ,IDWMAX  ,IDDUM   ,ID      ,IS      ,ISFPM
!
   REAL     ETOTW   ,FPM     ,SIG     ,ATOTD

   REAL     AC2(MDC,MSC,0:MT)
!
!  *** 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                                                         
   facint = 0.                                                         
   DO is = 1, msc
     sig = spcsig(is)                                                  
     IF(frinth * sig > (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)                                                  
     atotd = 0.
     DO iddum = idwmin, idwmax
       id = mod( iddum - 1 + mdc, mdc ) + 1
       atotd = atotd + ac2(id,is,kcgrd(1))                             
     ENDDO
     IF(is == 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

   RETURN

Here is the caller graph for this function:

windp3()

subroutine windp3	(	integer	ISSTOP,
		real, dimension(mdc,msc)	ALIMW,
		real, dimension(mdc,msc,0:mt)	AC2,
		logical, dimension(mdc,msc)	GROWW,
		integer, dimension(msc)	IDCMIN,
		integer, dimension(msc)	IDCMAX
	)

Definition at line 331 of file swancom3.f90.

!    (This subroutine has not been tested yet)
!
!****************************************************************
!
      USE swcomm3                                                        
      USE swcomm4                                                        
      USE ocpcomm4  
      USE all_vars, ONLY : mt                                                     
!
!
!     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,0:MT)    ,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

Here is the call graph for this function:

wndpar()

subroutine wndpar	(	integer	ISSTOP,
		integer	IDWMIN,
		integer	IDWMAX,
		integer, dimension(msc)	IDCMIN,
		integer, dimension(msc)	IDCMAX,
		real, dimension(mt)	DEP2,
		real	WIND10,
		real	THETAW,
		real, dimension(mdc,msc,0:mt)	AC2,
		real, dimension(msc,micmax)	KWAVE,
		real, dimension(mdc,msc)	IMATRA,
		real, dimension(mdc,msc)	IMATDA,
		real, dimension(msc)	SPCSIG,
		real, dimension(msc,micmax)	CGO,
		real, dimension(mdc,msc)	ALIMW,
		logical, dimension(mdc,msc)	GROWW,
		real	ETOTW,
		real, dimension(mdc,msc,nptst)	PLWNDA,
		real, dimension(mdc,msc,nptst)	PLWNDB,
		real, dimension(mdc,6)	SPCDIR,
		integer	ITER
	)

Definition at line 1052 of file swancom3.f90.

!  (This subroutine has not been tested yet)              
!
!****************************************************************
!
   USE swcomm3                                                         
   USE swcomm4                                                         
   USE ocpcomm4 
   USE all_vars, ONLY : mt                                                       
!
   IMPLICIT NONE                                                       
!
!
!  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)                        
!             (*,2); cosine of spectral directions                        
!             (*,3); sine of spectral directions                          
!             (*,4); cosine^2 of spectral directions                      
!             (*,5); cosine*sine of spectral directions                   
!             (*,6); sine^2 of spectral directions                        
! i   SPCSIG: Relative frequencies in computational domain in sigma-space 
!
   REAL    SPCDIR(MDC,6)                                               
   REAL    SPCSIG(MSC)                                                 
!
!  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
!        PLWNDA    3D    Values of source term for test point
!        PLWNDB    3D    Values of source term 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,                                  &
        DIRDIS,AC2CEN,DTHETA
!
   REAL  :: AC2(MDC,MSC,0:MT)
   REAL  :: ALIMW(MDC,MSC)
   REAL  :: IMATDA(MDC,MSC), IMATRA(MDC,MSC)
!  Changed ICMAX to MICMAX, since MICMAX doesn't vary over gridpoint   
   REAL  :: KWAVE(MSC,MICMAX)                                          
   REAL  :: PLWNDA(MDC,MSC,NPTST)                                      
   REAL  :: PLWNDB(MDC,MSC,NPTST)
   REAL  :: DEP2(MT)
!  Changed ICMAX to MICMAX, since MICMAX doesn't vary over gridpoint   
   REAL  :: CGO(MSC,MICMAX)                                            
!
   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_w
!   DND    = MIN( 50. , GRAV_W * DEP2(IGC) / WIND10**2 )
   dnd    = min( 50. , grav_w * dep2(kcgrd(1)) / wind10**2 )
   sigpk  = twopi * 0.13 * grav_w / wind10
   sigpkd = sigpk / tanh(0.833*dnd**0.375)
   fpm    = sigpkd                                                     
   ctw    = cos(thetaw)                                                
   stw    = sin(thetaw)                                                
!
   IF(iwind == 1)THEN
!
!    *** first generation model ***
!
     alpm = 0.0081
!
   ELSE IF(iwind == 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               )    

     edml = min( pwind(10) , (grav_w**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) )                         
     temp2  = exp( (-5./4.) * arg2**4 )
     alim1d = temp1 * temp2 / spcsig(is)                               
     DO iddum = idwmin, idwmax
       id     = mod( iddum - 1 + mdc, mdc ) + 1
       theta  = spcdir(id,1)                                           
       cosdif = spcdir(id,2)*ctw + spcdir(id,3)*stw                    
!
!  For better convergence the first guess of the directional spreading 
!  is modified in third generation mode. The new formulation better    
!  fits the directional spreading of the deep water growth curves.     
!
       IF((iter == 1).AND.(igen == 3))THEN                           
         dirdis = 0.434917 * (max(0., cosdif))**0.6                    
       ELSE                                                            
         dirdis = (2./pi_w) * cosdif**2                                  
       END IF                                                          
!
       alimw(id,is) = alim1d * dirdis
!       AC2CEN       = AC2(ID,IS,IGC)
       ac2cen       = ac2(id,is,kcgrd(1))
       IF(ac2cen <= alimw(id,is))THEN
         groww(id,is) = .true.
       ELSE
         groww(id,is) = .false.
       ENDIF
     ENDDO
!    *** test output ***
!     IF(TESTFL .AND. ITEST >= 10)THEN
!       WRITE(PRINTF,2002) IS, SPCSIG(IS), KWAVE(IS,1), CGO(IS,1)       
!2002   FORMAT(' WNDPAR: IS SPCSIG KWAVE CGO :',I3,3E12.4)              
!       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_w * pwind(9)**2 * pwind(11)**2 / grav_w**2
!
   DO is = 1, isstop
     sigma   = spcsig(is)                                              
     sigtpi  = sigma * twopi
     cinv    = kwave(is,1) / sigma
     factb = pwind(2) * pwind(9) * sigma / twopi                       
     DO iddum = idcmin(is), idcmax(is)
       id     = mod( iddum - 1 + mdc, mdc ) + 1
       dtheta = spcdir(id,1) - thetaw                                  
       cosdif = spcdir(id,2)*ctw + spcdir(id,3)*stw                    
!       AC2CEN = AC2(ID,IS,IGC)
       ac2cen = ac2(id,is,kcgrd(1))
!
       swind_exp = 0.                                                  
       swind_imp = 0.                                                  

       IF(groww(id,is))THEN
!        *** term A ***
         IF(sigma >= ( 0.7 * sigpkd ))THEN
           adum = facta * (wind10 * cosdif)**4                         
           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)))          
         bdum = factb * bdum * 5.
         swind_exp = adum + bdum * ac2cen                              
!
       ELSE IF(.NOT. groww(id,is) .AND. ac2cen > 0.)THEN
!
!        *** for no energy dissipation outside the wind field     ***
!        *** TAUINV is set equal zero (as in the DOLPHIN-B model) ***
!
         IF(cosdif < 0.) THEN                                    
           tauinv = 0.
         ELSE
           tauinv = ( sigma**2 * wind10 * abs(cosdif) ) /        &
                ( pwind(4) * grav_w * 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) plwnda(id,is,iptst) = swind_exp                     
       IF (testfl) plwndb(id,is,iptst) = swind_imp                     
!
!
!      *** test output ***
!
!      Value of ITEST changed from 10 to 110 to reduce test output     
!       IF ( TESTFL .AND. ITEST .GE. 110 ) THEN                        
!         WRITE(PRINTF,2004) IS, ID, GROWW(ID,IS), ADUM, BDUM           
!2004     FORMAT(' WNDPAR: IS ID GROWW ADUM BDUM     :',            &
!                2I3,2X,L1,2X,2E12.4)
!       END IF
     ENDDO
   ENDDO
!
!  *** test output ***
!
!  Value of ITEST changed from 10 to 60 to reduce test output          
!   IF(TESTFL .AND. ITEST >= 60)THEN                              
!     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

Here is the call graph for this function: