swancom3.f90

Go to the documentation of this file.

!
!****************************************************************
!
   SUBROUTINE windp1 (WIND10     ,THETAW     ,IDWMIN     ,      &
                      IDWMAX     ,FPM        ,UFRIC      ,      &
              WX2        ,WY2        ,SPCDIR     ,      &
              UX2        ,UY2        ,SPCSIG     ,      &
              IG         )  
!
!****************************************************************
!
!        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
   END SUBROUTINE windp1
 
!
!****************************************************************
!
   SUBROUTINE windp2 (IDWMIN  ,IDWMAX  ,SIGPKD  ,FPM     ,      &
                      ETOTW   ,                                 &
              AC2     ,SPCSIG  ,                        &
              WIND10                                      )    
!  (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
   END SUBROUTINE windp2
!
!********************************************************************
!
      SUBROUTINE windp3 (ISSTOP  ,ALIMW   ,AC2     ,               &
                         GROWW   ,IDCMIN  ,IDCMAX  )                 
!    (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
      END SUBROUTINE windp3
!
!****************************************************************
!
   SUBROUTINE swind0 (IDCMIN  ,IDCMAX  ,ISSTOP  ,SPCSIG  ,THETAW  , &
              UFRIC   ,FPM     ,PLWNDA  ,IMATRA  ,SPCDIR  )
!
!****************************************************************
!
!     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
   END SUBROUTINE swind0

!
!****************************************************************
!
   SUBROUTINE swind3 (SPCSIG  ,THETAW  ,IMATDA  ,KWAVE   ,IMATRA  , &
              IDCMIN  ,IDCMAX  ,UFRIC   ,FPM     ,PLWNDB  , &
              ISSTOP  ,SPCDIR  ,IG      )
!
!****************************************************************
!
!     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
   END SUBROUTINE swind3
 
!
!****************************************************************
!
   SUBROUTINE swind4 (IDWMIN  ,IDWMAX  ,SPCSIG  ,WIND10  ,THETAW  ,  &
                      XIS     ,DD      ,KWAVE   ,IMATRA  ,IMATDA  ,  &
              IDCMIN  ,IDCMAX  ,UFRIC   ,PLWNDB  ,ISSTOP  ,  &
              ITER    ,USTAR   ,ZELEN   ,SPCDIR  ,IT      ,  &
              IG      )
!
!******************************************************************
!
!     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
   END SUBROUTINE swind4
 
!
!****************************************************************
!
   SUBROUTINE swind5 (SPCSIG  ,THETAW  ,ISSTOP  ,UFRIC   ,KWAVE   , &
                      IMATRA  ,IMATDA  ,IDCMIN  ,IDCMAX  ,PLWNDB  , &
              SPCDIR  ,IG      )
!
!****************************************************************
!
!     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
   END SUBROUTINE swind5
!
 
!
!****************************************************************
!
   SUBROUTINE wndpar (ISSTOP,IDWMIN,IDWMAX,IDCMIN,IDCMAX,       &
                      DEP2  ,WIND10,THETAW,AC2   ,KWAVE ,       &
              IMATRA,IMATDA,SPCSIG,CGO   ,ALIMW ,       &
              GROWW ,ETOTW ,PLWNDA,PLWNDB,SPCDIR,       &
              ITER  )   
!  (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

   END SUBROUTINE wndpar