C **********************************************************************
C *                                                                    *
C *                                SOFTWARE LICENSING                  *
C *                                                                    *
C *                This program is free software; you can redistribute *
C *                it and/or modify it under the terms of the GNU      *
C *                General Public License as published by the Free     *
C *                Software Foundation, either Version 2 of the        *
C *                license, or (at your option) any later version.     *
C *                                                                    *
C *                This program is distributed in the hope that it     *
C *                will be useful, but without any warranty; without   *
C *                even the implied warranty of merchantability or     *
C *                fitness for a particular purpose. See the GNU       *
C *                General Public License for more details.            *
C *                                                                    *
C *                A copy of the GNU General Public License is         *
C *                available at http://www.gnu.org/copyleft/gpl.html   *
C *                or by writing to the Free Software Foundation, Inc.,*
C *                59 Temple Place - Suite 330, Boston, MA 02111, USA. *
C *                                                                    *
C **********************************************************************

      SUBROUTINE WAVESMB
C
C  CALC. BOTTOM SHEAR STRESS DUE TO WIND WAVES AND CURRENTS
C  USES GRANT/MADSEN/GLENN WAVE MODEL
C  REVISION DATE:  SEPTEMBER 26, 1997
C**************************************************************
C
      INCLUDE 'comdeck'
      SAVE
C
csfan moved to comdeck      dimension FSMW(IM,JM)

      PI=  3.141593
      TPI= 6.283185
      HPI= 1.570796
C
C    The EXTERNAL should be taken care of later
cqa      IF (WAVEDYN.EQ.'EXTERNAL') GOTO 15   
C
C  CALM WIND ( < 1 MPH = 0.45 m/s) ==> PURE CURRENTS
C
CPL      IF (WINDSP.LT.0.45) GOTO 99
C
C  CALC. SIGNIFICANT WAVE HEIGHT AND PERIOD
C  USE SHALLOW WATER SMB THEORY
C
C  WVHT = SIG. WAVE HEIGHT (m)
C  WVPD = SIG. WAVE PERIOD (s)
C
C  WINDSP = WIND SPEED (m/s)
C  WINDIR = WIND DIRECTION (degrees)
C  HMEAN(N) = MEAN DEPTH ALONG FETCH FOR DIRECTION n (m)
C  FETCH(N) = FETCH FOR DIRECTION n (m)
C
C  Previously WINDSP and AINDIR was passed into this subroutine calculated
C  in bcond.f was wrong. We introduce spatially variable wind field and
C  incorporate correct interpolation to compute WSX2D and WSY2D at time THOUR
C  by Quamrul QA 10/27/98
C
chli
      DO J=1,JM
         DO I=1,IM
            FSMW(I,J)=FSM(I,J)
         ENDDO
      ENDDO
      DO N=1,NUMQBC
         IC=IQC(N)
         JC=JQC(N)
         FSMW(IC,JC)=0.0
      ENDDO
chli
      DO 10 J=2,JM-1
        DO 10 I=2,IM-1
          IF (FSMW(I,J).GT.0.0) THEN
C
CPL 4/20/99 ADD SMALL VALUE
C
           UWIND=(WSX2D(I,J)*COS(ANG(I,J))+WSY2D(I,J)*SIN(ANG(I,J)))*
     &        FSMW(I,J)
           VWIND=(-WSX2D(I,J)*SIN(ANG(I,J))+WSY2D(I,J)*COS(ANG(I,J)))*
     &        FSMW(I,J)
           WSPDSQ=UWIND*UWIND+VWIND*VWIND
           WINDSP=SQRT(WSPDSQ)+1.E-10
C
C     WINDIR is w.r.t. zeta1-direction
C
           WINDIR=ATAN2(VWIND,UWIND)*180./PI
           IF(WINDIR .LT. 0.) WINDIR = WINDIR + 360
           IF(WINDIR .GT. 360.) WINDIR = WINDIR - 360.
C
C     Direction with respect to east (x-direction)
C
           WINDIR = WINDIR + ANG(I,J)*180./PI
           WINDIR = AMOD(WINDIR,360.)
C
C     IN stress.f and in SMB Model WAVE DIR IS IN w.r.t North Clock Wise
C
           WINDIR= 90. - WINDIR
           IF(WINDIR.LE.0.) WINDIR = WINDIR + 360.
           NDIR=NINT(WINDIR/10.)
           IF (NDIR.EQ.0) NDIR=36
           WRAT0=GRAV/(WINDSP*WINDSP)
           HRAT0=0.283/WRAT0
           TRAT0=2.4*PI*WINDSP/GRAV
           WRAT1=WRAT0*HMEAN(NDIR,I,J)
           WRAT2=WRAT0*FETCH(NDIR,I,J)
           HA1=TANH(0.53*(WRAT1**0.75))
           TA1=TANH(0.833*(WRAT1**0.375))
           WVHT(I,J)=HRAT0*HA1*TANH(0.0125*(WRAT2**0.42)/HA1)
           WVPD(I,J)=TRAT0*TA1*TANH(0.077*(WRAT2**0.25)/TA1)
           WVDR(I,J)=WINDIR
C
C  LIMIT WAVE HEIGHT TO BREAKING HEIGHT
C
            TBREAK=TIME*24.0
            HBREAK=0.78*DT(I,J)
            WVHT(I,J)=AMIN1(HBREAK,WVHT(I,J))
          ENDIF
 10   CONTINUE
C
C  OUTPUT WAVE PARAMETERS IF WAVEDYN IS SMBMODEL
C
      WRITE(1001)TIME
      WRITE(1001)WVHT
      WRITE(1001)WVPD
      WRITE(1001)WVDR
C
      RETURN
      END