SUBROUTINE f4spec (alfa, fm, gamma, sa, thetaq, ustar, spint, spt)
!----------------------------------------------------------------------
!**** *F4SPEC* - PRODUCE A WINDSEA SPECTRUM FOR ASSIMILATION.
!      P. LIONELLO    ECMWF        FEBRUARY 1990
!      (DERIVED FROM JSPEC OF PETER JANSSEN  1987.)
!**   INTERFACE.
!     ----------
!       *CALL* *F4SPEC (ALFA, FM, GAMMA, SA, THETAQ, USTAR, SPINT,SPT)*
!          *ALFA*  - TOBAS' CONSTANT
!          *FM*    - PEAK FREQUENCY
!          *GAMMA* - OVERSHOOT PARAMETER
!          *SA*    - WIDTH PARAMETER
!          *THETAQ*- MEAN WINDSEA DIRECTION
!          *USTAR* - FRICTION VELOCITY
!          *SPINT* - TOTAL ENERGY
!          *SPT*   - 2-D SPECTRUM
!      METHOD.
!      -------
!      THE ONE DIMENSIONAL PARAMETRIC SPECTRUM IS EVALUATED AND
!      DISTRIBUTED AROUND THE MEAN SEA DIRECTION.
!      THE TOTAL ENERGY IS COMPUTED.
!      EXTERNALS.
!      ----------
!         *SPR*  -  SUB TO COMPUTE COS**2 SPREADING FUNCTION.
! ---------------------------------------------------------------------
use params
use fredir
implicit none

REAL, INTENT(IN) :: alfa, gamma, thetaq, ustar, fm, sa
REAL, INTENT(OUT):: spint, spt(nang,nfre)
INCLUDE'globals.h' !NIKOS ONLY G,PI,ZPI

REAL :: delfr(nfre), et(nfre), st(nang)
!      *SPT*        REAL    2-D SPECTRUM.
!      *DELFR*      REAL    FREQUENCY INCREMENTS.
!      *ET*         REAL    1-D JONSWAP SPECTRUM.
!      *ST*         REAL    SPREADING FUNCTION.
INTEGER :: kl, ml, kl2, ml2, k, m
REAL    :: fakal, frh, earg, fjon, fmpf, fjonh
!----------------------------------------------------------------------
ml = nfre
kl = nang
!*    1. COMPUTE FREQUENCY INTERVALLS.
!        -----------------------------
DO m = 1, ml
   delfr(m) = fr(m) * 0.1
ENDDO
!----------------------------------------------------------------------
!*    2. COMPUTE TOBA SPECTRUM.
!        ----------------------
fakal = alfa * g * ustar / (zpi**3)
DO m = 1, ml
   frh   = fr(m)
   earg  = (.5 * ( (frh - fm) / (sa * fm) ) **2)
   fjon  = GAMMA**EXP ( - earg)
   fmpf  = 4. / 6. * (fm / frh) **6
   fjonh = EXP ( - fmpf)
   et(m) = fakal * fjonh * fjon / frh**4
ENDDO
!----------------------------------------------------------------------
!*  3. COMPUTATION OF SPREADING FUNCTION.
!      ----------------------------------
CALL spr (kl, nang, thetaq, th, st)
!----------------------------------------------------------------------
!*  4. COMPUTATION OF 2-D SPECTRUM.
!      ---------------------------
DO m = 1, ml
   DO k = 1, kl
      spt(k, m) = et(m) * st(k)
   ENDDO
ENDDO
! ---------------------------------------------------------------------
!*  5. COMPUTE TOTAL ENERGY.
!      --------------------
!*  5.1 INTEGRATE OVER DIRECTIONS.
!       --------------------------
kl2 = kl - 2
spint = 0.
DO m = 1, ml
   et(m) = spt(1, m)
   DO k = 2, kl2, 2
      et(m) = et(m) + 4. * spt(k, m) + 2. * spt(k + 1, m)
   ENDDO
   et(m) = (et(m) + 4. * spt(kl, m) + spt(1, m) ) * delth / 3.
ENDDO
!*  5.2 INTEGRATE 1-D SPECTRUM OVER FREQUENCIES.
!       ----------------------------------------
spint = et(1) * delfr(1)
ml2 = ml - 2
DO m = 2, ml2, 2
   spint = spint + (4. * et(m) * delfr(m) + 2. * et(m+1) * delfr(m+1))
ENDDO
IF (2 * (ML / 2) .EQ.ML) THEN
   spint = (spint + 4. * et(ml) * delfr(ml) ) / 3.
ELSE
   spint = (spint + 4. * et(ml-1) * delfr(ml-1) + et(ml) * delfr(ml)) / 3.
ENDIF

RETURN
END SUBROUTINE f4spec
