SUBROUTINE stokes_drift (f, ijs, ijl, ust, vst)  
! ----------------------------------------------------------------------
!*    PURPOSE.
!     --------
!       COMPUTE STOKES DRIFT VELOCITY COMPONENTS AT EACH GRID POINT.
!**   INTERFACE.
!     ----------
!       *CALL* *STOKES_DRIFT (F, IJS, IJL, UST, VST)*
!          *F*     - SPECTRUM.
!          *IJS*   - INDEX OF FIRST GRIDPOINT
!          *IJL*   - INDEX OF LAST GRIDPOINT
!          *UST* - X-COMPONENT OF STOKES DRIFT VELOCITY
!          *VST* - Y-COMPONENT OF STOKES DRIFT VELOCITY
!     METHOD.
!     -------
!       NONE.
!     EXTERNALS.
!     ----------
!       NONE.
!     REFERENCE.
!     ----------
!       NONE.
! ----------------------------------------------------------------------
use params
use fredir
implicit none

REAL   , INTENT(IN) :: f(0:niblo,nang,nfre)
REAL   , INTENT(OUT):: ust(niblo), vst(niblo)
INTEGER, INTENT(IN) :: ijs, ijl
INCLUDE'globals.h'!NIKOS ONLY G,PI
REAL    :: TEMPU(NIBLO), TEMPV(NIBLO)
REAL    :: FAK, FAKU, FAKV, FAKT, TAILFAC
INTEGER :: K, M, IJ
! ----------------------------------------------------------------------
!*    1. INITIALIZE ARRAYS
!        -----------------
FAK = 16 * PI**3 / G  
TAILFAC=FAK*FR(NFRE)**4


DO IJ = IJS, IJL  
   UST (IJ) = 0.  
   VST (IJ) = 0.  
END DO  
DO M = 1, NFRE  
   FAKT = FAK * FR (M) **3 * DFIM (M)  
   DO IJ = IJS, IJL  
      TEMPU (IJ) = 0.  
      TEMPV (IJ) = 0.  
   END DO  
   DO K = 1, NANG  
      FAKU = FAKT * SIN (TH (K) )  
      FAKV = FAKT * COS (TH (K) )  
      DO 4002 IJ = IJS, IJL  
         TEMPU (IJ) = TEMPU (IJ) + FAKU * F (IJ, K, M)  
 4002       TEMPV (IJ) = TEMPV (IJ) + FAKV * F (IJ, K, M)  
 4004    END DO  
   DO 4005 IJ = IJS, IJL  
      UST (IJ) = UST (IJ) + TEMPU (IJ)  
      VST (IJ) = VST (IJ) + TEMPV (IJ)  
 4005    END DO  
END DO  

!*    3. ADD CONTRIBUTION FROM TAIL.
DO K = 1, NANG
            FAKU = DELTH*SIN(TH(K))
            FAKV = DELTH*COS(TH(K))
DO IJ = IJS, IJL  
              UST(IJ)=UST(IJ)+TAILFAC*FAKU*F(IJ,K,NFRE)
              VST(IJ)=VST(IJ)+TAILFAC*FAKV*F(IJ,K,NFRE)
END DO  
END DO  

RETURN  
END SUBROUTINE stokes_drift
