SUBROUTINE sdissip (f, fl, diss,ijs, ijl)
! ----------------------------------------------------------------------
!**** *SDISSIP* - COMPUTATION OF DISSIPATION SOURCE FUNCTION.
!     S.D.HASSELMANN.
!     MODIFIED TO SHALLOW WATER : G. KOMEN , P. JANSSEN
!     OPTIMIZATION : L. ZAMBRESKY
!*    PURPOSE.
!     --------
!       COMPUTE DISSIPATION SOURCE FUNCTION AND STORE ADDITIVELY INTO
!       NET SOURCE FUNCTION ARRAY. ALSO COMPUTE FUNCTIONAL DERIVATIVE
!       OF DISSIPATION SOURCE FUNCTION.
!**   INTERFACE.
!     ----------
!       *CALL* *SDISSIP (F, FL, IJS, IJL)*
!          *F*   - SPECTRUM.
!          *FL*  - DIAGONAL MATRIX OF FUNCTIONAL DERIVATIVE
!          *IJS* - INDEX OF FIRST GRIDPOINT
!          *IJL* - INDEX OF LAST GRIDPOINT
!     METHOD.
!     -------
!       SEE REFERENCES.
!     EXTERNALS.
!     ----------
!       NONE.
!     REFERENCE.
!     ----------
!       G.KOMEN, S. HASSELMANN AND K. HASSELMANN, ON THE EXISTENCE
!          OF A FULLY DEVELOPED WINDSEA SPECTRUM, JGR, 1984.
! ----------------------------------------------------------------------
USE params ; USE fredir ; USE meanpa ; USE source ; USE stat
!SHALLOW
USE shallow
!SHALLOW
IMPLICIT NONE

REAL   , INTENT(IN) :: f(0:niblo,nang,nfre)
REAL   , INTENT(OUT):: fl(0:niblo,nang,nfre)
REAL   , INTENT(OUT):: diss(0:niblo,nang,nfre)
INTEGER, INTENT(IN) :: ijs, ijl
INCLUDE'globals.h'
REAL,ALLOCATABLE:: temp1(:), sds(:)
REAL, PARAMETER :: cdis = 4.5
REAL, PARAMETER :: consd = - 0.5 * cdis * zpi**9 / g**4
REAL, PARAMETER :: conss = - 0.5 * cdis * zpi
INTEGER :: k, m, ij
REAL    :: x
! ----------------------------------------------------------------------
ALLOCATE(temp1(ijs:ijl)) ; ALLOCATE(sds(ijs:ijl))
!*    1. ADDING DISSIPATION AND ITS FUNCTIONAL DERIVATIVE TO NET SOURCE
!*       FUNCTION AND NET SOURCE FUNCTION DERIVATIVE.
!        --------------------------------------------------------------
IF( ishallo .EQ. 1 ) THEN
   DO ij = ijs, ijl
      sds(ij) = consd * emean(ij) **2 * fmean(ij) **9
   ENDDO
   DO m = 1, nfre
      DO ij = ijs, ijl
         x         = (fr(m) / fmean(ij) ) **2
         temp1(ij) = sds(ij) * (x + x**2)
      ENDDO
      DO k = 1, nang
         DO ij = ijs, ijl
            sl(ij,k,m) = sl(ij,k,m) + temp1(ij) * f(ij,k,m)
            diss(ij,k,m) =  -temp1(ij) * f(ij,k,m)
            fl(ij,k,m) = fl(ij,k,m) + temp1(ij)
         ENDDO
      ENDDO
   ENDDO
ELSE
!SHALLOW
   DO ij = ijs, ijl
      sds(ij) = conss * fmean(ij) * emean(ij)**2 * akmean(ij) **4
   ENDDO
   DO m = 1, nfre
      DO ij = ijs, ijl
         x         = tfak(indep(ij), m) / akmean(ij)
         temp1(ij) = sds(ij) * (x + x**2)
      ENDDO
      DO  k = 1, nang
         DO ij = ijs, ijl
            sl(ij,k,m) = sl(ij,k,m) + temp1(ij) * f(ij,k,m)
            diss(ij,k,m) =  -temp1(ij) * f(ij,k,m)
            fl(ij,k,m) = fl(ij,k,m) + temp1(ij)
         ENDDO
      ENDDO
   ENDDO
!SHALLOW
ENDIF
DEALLOCATE(sds,temp1)

RETURN
END SUBROUTINE sdissip
