SUBROUTINE peakfr (f, ijs, ijl, peakf)
! ----------------------------------------------------------------------
!**** *PEAKFR* - COMPUTATION OF PEAK FREQUENCY AT EACH GRID POINT.
!     H. GUNTHER      ECMWF            DECEMBER 1989
!     (CODE REMOVED FROM SUB. FEMEAN)
!*    PURPOSE.
!     --------
!       COMPUTE PEAK FREQUENCY AT EACH GRID POINT.
!**   INTERFACE.
!     ----------
!       *CALL* *PEAKFR (F, IJS, IJL, PEAKF)*
!          *F*     - SPECTRUM.
!          *IJS*   - INDEX OF FIRST GRIDPOINT
!          *IJL*   - INDEX OF LAST GRIDPOINT
!          *PEAKF* - PEAKFREQUENCY
!     METHOD.
!     -------
!       NONE.
!     EXTERNALS.
!     ----------
!       NONE.
!     REFERENCE.
!     ----------
!       NONE.
! ----------------------------------------------------------------------
use params
use fredir
implicit none

REAL   , INTENT(IN) :: f(0:niblo,nang,nfre)
REAL   , INTENT(OUT):: peakf(niblo)
INTEGER, INTENT(IN) :: ijs, ijl
REAL,ALLOCATABLE :: temp1(:), temp2(:)
INTEGER :: K, M, IJ
! ----------------------------------------------------------------------
ALLOCATE(temp1(ijs:ijl)); ALLOCATE(temp2(ijs:ijl))
!*    1. INITIALIZE ARRAYS
!        -----------------
DO IJ = IJS, IJL
   TEMP1 (IJ) = 0.
   PEAKF (IJ) = 0.
END DO
!*    2. LOOP OVER FREQUENCIES
!        ---------------------
DO M = 1, NFRE
!*    2.1 COMPUTE 1-D SPECTRUM (WITHOUT DELTA THETA)
!        -------------------------------------------
   DO IJ = IJS, IJL
      TEMP2 (IJ) = 0.
   END DO
   DO K = 1, NANG
      DO IJ = IJS, IJL
         TEMP2 (IJ) = TEMP2 (IJ) + F (IJ, K, M)
      END DO
   END DO
!*    2.2 DEFINE PEAK FREQUENCY
!         ---------------------
   DO IJ = IJS, IJL
      IF (TEMP1 (IJ) .LT.TEMP2 (IJ) ) THEN
         TEMP1 (IJ) = TEMP2 (IJ)
         PEAKF (IJ) = FR (M)
      ENDIF
   END DO
END DO
DEALLOCATE(temp2, temp1)

RETURN
END SUBROUTINE peakfr
