SUBROUTINE stresso (f, ijs, ijl, ig, igl)
! ----------------------------------------------------------------------
!**** *STRESSO* - COMPUTATION OF WAVE STRESS.
!     H. GUNTHER      GKSS/ECMWF  NOVEMBER  1989 CODE MOVED FROM SINPUT.
!     P.A.E.M. JANSSEN      KNMI  AUGUST    1990
!*    PURPOSE.
!     --------
!       COMPUTE NORMALIZED WAVE STRESS FROM INPUT SOURCE FUNCTION
!**   INTERFACE.
!     ----------
!       *CALL* *STRESSO (F, IJS, IJL, IG, IGL)*
!          *F*   - WAVE SPECTRUM.
!          *IJS* - INDEX OF FIRST GRIDPOINT.
!          *IJL* - INDEX OF LAST GRIDPOINT.
!          *IG*  - ACTUAL BLOCK NUMBER.
!          *IGL* - NUMBER OF BLOCKS.
!     METHOD.
!     -------
!       THE INPUT SOURCE FUNCTION IS INTEGRATED OVER FREQUENCY
!       AND DIRECTIONS.
!       BECAUSE ARRAY *SL* IS USED, ONLY THE INPUT SOURCE
!       HAS TO BE STORED IN *SL* (CALL FIRST SINPUT, THEN
!       STRESSO, AND THEN THE REST OF THE SOURCE FUNCTIONS)
!     EXTERNALS.
!     -----------
!       NONE.
!     REFERENCE.
!     ----------
!       R SNYDER ET AL,1981.
!       G. KOMEN, S. HASSELMANN AND K. HASSELMANN, JPO, 1984.
!       P. JANSSEN, JPO, 1985
! ----------------------------------------------------------------------
USE params ; USE coupl ; USE fredir ; USE source ; USE table ; USE wind
IMPLICIT NONE

REAL   , INTENT(IN):: f(0:niblo,nang,nfre)
INTEGER, INTENT(IN):: ig, ijs, ijl, igl
INCLUDE'globals.h'
REAL,PARAMETER :: roair = 1.225, rowater = 1000.
REAL,PARAMETER :: xeps = roair / rowater, xinveps = 1. / xeps

REAL,ALLOCATABLE :: constf(:), tauhf(:), temp(:), xstress(:), ystress(:)
REAL    :: const0, const1, const2, deli1, deli2, delj1, delj2
REAL    :: tau1, tkd, cosw, ust, xi, xj
INTEGER :: i, j, k, m, ij
! ----------------------------------------------------------------------
ALLOCATE(constf(nfre))  ; ALLOCATE(tauhf(ijs:ijl))
ALLOCATE(temp(ijs:ijl)) ; ALLOCATE(xstress(ijs:ijl)) ; ALLOCATE(ystress(ijs:ijl))
!*    1. PRECOMPUTE FREQUENCY SCALING.
!        -----------------------------
const0 = delth * (zpi) **4 * fr(nfre) **5 / g**2

DO m = 1, nfre
   constf(m) = zpi * xinveps * fr(m) * dfim(m)
ENDDO
!*    2. COMPUTE WAVE STRESS OF ACTUEL BLOCK.
!        ------------------------------------
!*    2.1 PRESET STRESS ARRAYS.
!         ---------------------
DO ij = ijs, ijl
   xstress(ij) = 0.
   ystress(ij) = 0.
ENDDO
!*    2.2 INTEGRATE INPUT SOURCE FUNCTION OVER FREQUENCY AND DIRECTIONS.
!         --------------------------------------------------------------
DO k = 1, nang
   DO m = 1, nfre
      const1 = constf(m) * sinth(k)
      const2 = constf(m) * costh(k)
      DO ij = ijs, ijl
         xstress(ij) = xstress(ij) + sl(ij,k,m) * const1
         ystress(ij) = ystress(ij) + sl(ij,k,m) * const2
      ENDDO
   ENDDO
ENDDO
!*    2.3 CALCULATE HIGH-FREQUENCY CONTRIBUTION TO STRESS.
!     ----------------------------------------------------
DO ij = ijs, ijl
   temp(ij) = 0.
ENDDO

DO k = 1, nang
   tkd = th(k)
   DO ij = ijs, ijl
      cosw     = MAX (COS (tkd-thwnew(ij) ), 0.)
      temp(ij) = temp(ij) + f(ij, k, nfre) * cosw**3
   ENDDO
ENDDO

DO ij = ijs, ijl
! KORRES MODIFICATION 1
!  ust   = usnew(ij)
   ust   = max(usnew(ij),0.000001)
   xi    = ust / delust
! KORRES MODIFICATION 2
   XI    = MIN(REAL(IUSTAR),XI)
   i     = MIN (iustar - 1, INT (xi) )
! KORRES MODIFICATION 3
   i     = MAX (0, i)
   deli1 = MIN (1., xi - FLOAT (i) )
   deli2 = 1. - deli1

   xj    = (g * z0new(ij) / ust**2 - alpha) / delalp
! KORRES MODIFICATION 4
        XJ    = MIN(REAL(IALPHA),XJ)

   j     = MIN (ialpha - 1, INT (xj) )
! KORRES MODIFICATION 5
   j     = MAX (0 , J)
! KORRES MODIFICATION 6
!  delj1 = MIN (1., xj - FLOAT (j) )
        DELJ1 = MAX(MIN (1. ,XJ-FLOAT(J)),0.)
   delj2 = 1. - delj1

   tau1 = (tauhft(i, j)     * deli2 + tauhft(i + 1, j) * deli1) * delj2    &
         +(tauhft(i, j + 1) * deli2 + tauhft(i + 1, j + 1) * deli1) * delj1
   tauhf(ij) = const0 * temp(ij) * ust**2 * tau1
ENDDO

DO ij = ijs, ijl
   xstress(ij)  = xstress(ij) + tauhf(ij) * SIN(thwnew(ij) )
   ystress(ij)  = ystress(ij) + tauhf(ij) * COS(thwnew(ij) )
   tauw(ij, ig) = SQRT (xstress(ij) **2 + ystress(ij) **2)

   tauw_x(ij, ig) = xstress(ij) 
   tauw_y(ij, ig) = ystress(ij) 
!  tauwhf_x(ij, ig) = tauhf(ij) * SIN(thwnew(ij) )
!  tauwhf_y(ij, ig) = tauhf(ij) * COS(thwnew(ij) )
ENDDO
!   print *,'WAVE STRESS: ', tauw(12150,1),usnew(12150)
DEALLOCATE(ystress, xstress, temp, tauhf,constf)

RETURN
END SUBROUTINE stresso
