SUBROUTINE femean (f, ijs, ijl)
! ----------------------------------------------------------------------
!**** *FEMEAN* - COMPUTATION OF MEAN FREQUENCY AT EACH GRID POINT.
!     S.D. HASSELMANN
!     MODIFIED : P.JANSSEN (INTEGRATION OF F**-4 TAIL)
!     OPTIMIZED BY : L. ZAMBRESKY AND H. GUENTHER
!*    PURPOSE.
!     --------
!       COMPUTE MEAN FREQUENCY AT EACH GRID POINT.
!**   INTERFACE.
!     ----------
!       *CALL* *FEMEAN (F, IJS, IJL)*
!              *F*   - SPECTRUM.
!              *IJS* - INDEX OF FIRST GRIDPOINT
!              *IJL* - INDEX OF LAST GRIDPOINT
!     METHOD.
!     -------
!       NONE.
!     EXTERNALS.
!     ----------
!       NONE.
!     REFERENCE.
!     ----------
!       NONE.
! ----------------------------------------------------------------------
USE params ; USE fredir ; USE meanpa ; USE stat ; USE shallow
IMPLICIT NONE

REAL   , INTENT(IN):: f(0:niblo,nang,nfre)
INTEGER, INTENT(IN):: ijs, ijl
INCLUDE'globals.h'
!SHALLOW
REAL,ALLOCATABLE :: temp1(:), temp2(:)
INTEGER:: ij, m, k
REAL   :: delt25, fd, del2
! ----------------------------------------------------------------------
ALLOCATE(temp1(ijs:ijl)) ; ALLOCATE(temp2(ijs:ijl))
!*    1. INITIALISE MEAN FREQUENCY ARRAY AND TAIL FACTOR.
!        ------------------------------------------------
DO ij = ijs, ijl
!        FMEAN(IJ) = 0.1E-180
   fmean(ij) = 0.1e-32
ENDDO
delt25 = 0.20 * delth
! ----------------------------------------------------------------------
!*    2. INTEGRATE OVER FREQUENCIES AND DIRECTIONS.
!        ------------------------------------------
IF( ishallo .EQ. 1 ) THEN
!*    2.1 DEEP WATER INTEGRATION.
!         -----------------------
   DO m = 1, nfre
      fd = dfim(m) / fr(m)
      DO ij = ijs, ijl
         temp2(ij) = 0.
      ENDDO
      DO k = 1, nang
         DO ij = ijs, ijl
            temp2(ij) = temp2(ij) + f(ij,k,m)
         ENDDO
      ENDDO
      DO ij = ijs, ijl
         fmean(ij) = fmean(ij) + fd * temp2(ij)
      ENDDO
   ENDDO
!SHALLOW
ELSE
!*    2.2 SHALLOW WATER INTEGRATION.
!         --------------------------
   DO ij = ijs, ijl
!           AKMEAN(IJ) = 0.1E-180
      akmean(ij) = 0.1e-32
   ENDDO
   DO m = 1, nfre
      fd = dfim(m) / fr(m)
      DO ij = ijs, ijl
         temp1(ij) = dfim(m) / SQRT (tfak(indep(ij),m))
         temp2(ij) = 0.
      ENDDO
      DO k = 1, nang
         DO ij = ijs, ijl
            temp2(ij) = temp2(ij) + f(ij,k,m)
         ENDDO
      ENDDO
      DO ij = ijs, ijl
         fmean(ij)  = fmean(ij) + fd * temp2(ij)
         akmean(ij) = akmean(ij) + temp1(ij) * temp2(ij)
      ENDDO
   ENDDO
!        ADD TAIL TO MEAN WAVENUMBER AND NORMALIZE WITH TOTAL ENERGY.
   del2 = delt25 * SQRT (g) / zpi
   DO ij = ijs, ijl
      akmean(ij) = akmean(ij) + del2 * temp2(ij)
      akmean(ij) = (emean(ij) / akmean(ij) ) **2
   ENDDO
ENDIF
!SHALLOW
!*    3. ADD TAIL CORRECTION TO MEAN FREQUENCY AND
!*       NORMALIZE WITH TOTAL ENERGY.
!        ------------------------------------------
DO ij = ijs, ijl
   fmean(ij) = fmean(ij) + delt25 * temp2(ij)
   fmean(ij) = emean(ij) / fmean(ij)
ENDDO
DEALLOCATE(temp2, temp1)

RETURN
END SUBROUTINE femean
