SUBROUTINE sthq (f3, ijs, ijl)  
! ----------------------------------------------------------------------
!**** *STHQ* - COMPUTATION OF MEAN WAVE DIRECTION AT EACH GRID POINT.
!     S.D. HASSELMANN
!     OPTIMIZED BY L. ZAMBRESKY
!*    PURPOSE.
!     --------
!       TO COMPUTE MEAN WAVE DIRECTION AT EACH GRID POINT.
!**   INTERFACE.
!     ----------
!       *CALL* *STHQ (F3, IJS, IJL)*
!          *F3*  - SPECTRUM.
!          *IJS* - INDEX OF FIRST GRIDPOINT
!          *IJL* - INDEX OF LAST GRIDPOINT
!     METHOD.
!     -------
!       INTEGRATION OF SPECTRUM TIMES SIN AND COS OVER DIRECTION.
!     EXTERNALS.
!     ----------
!       NONE.
!     REFERENCE.
!     ----------
!       NONE.
! ----------------------------------------------------------------------
use params
use fredir
use meanpa
implicit none

REAL   , INTENT(IN):: f3(0:niblo,nang,nfre)
INTEGER, INTENT(IN):: ijs, ijl
REAL    :: TEMP(NIBLO), SI(NIBLO), CI(NIBLO)
REAL    :: CIH, SIH
INTEGER :: K, M, IJ
INCLUDE'globals.h'
! ----------------------------------------------------------------------
!*    1. INITIALISE SIN AND COS ARRAYS.
!        ------------------------------
! 1000 CONTINUE
DO IJ = IJS, IJL
   SI (IJ) = 0.
   CI (IJ) = 0.
END DO
! ----------------------------------------------------------------------
!*    2. INTEGRATE OVER FREQUENCIES AND DIRECTIONS.
!        ------------------------------------------
! 2000 CONTINUE
DO K = 1, NANG
   CIH = COSTH (K)
   SIH = SINTH (K)
   DO IJ = IJS, IJL
      TEMP (IJ) = 0.
   END DO
   DO M = 1, NFRE
      DO IJ = IJS, IJL
         TEMP (IJ) = TEMP (IJ) + F3 (IJ, K, M) * DFIM (M)
      END DO
   END DO
   DO IJ = IJS, IJL
      SI (IJ) = SI (IJ) + SIH * TEMP (IJ)
      CI (IJ) = CI (IJ) + CIH * TEMP (IJ)
   END DO
END DO
! ----------------------------------------------------------------------
!*    3. COMPUTE MEAN DIRECTION.
!        -----------------------
! 3000 CONTINUE
DO IJ = IJS, IJL
!        IF (CI(IJ).EQ.0.) CI(IJ) = 0.1E-99
   IF (CI (IJ) .EQ.0.) CI (IJ) = 0.1E-32
 END DO
DO IJ = IJS, IJL
   THQ (IJ) = ATAN2 (SI (IJ), CI (IJ) )
END DO
DO IJ = IJS, IJL
   IF (THQ (IJ) .LT.0.) THQ (IJ) = THQ (IJ) + ZPI
END DO

RETURN
END SUBROUTINE sthq
