!
!****************************************************************
!
  SUBROUTINE SWAPAR(IG,NICMAX,DEP2,KWAVE,CGO)
!
!****************************************************************
!
!     computes the wave parameters K and CGO in the nearby
!     points, depending of the sweep direction.
!     The nearby points are indicated with the index IC (see
!     FUNCTION ICODE(_,_)
!
!  Method
!
!     The wave number K(IS,iC) is computed with the dispersion relation:
!
!     S = GRAV K(IS,IC)tanh(K(IS,IC)DEP(IX,IY))
!
!     where S = is logarithmic distributed via LOGSIG
!
!     The group velocity CGO in the case without current is equal to
!
!                    1       K(IS,IC)DEP(IX,IY)          S
!     CGO(IS,IC) = ( - + --------------------------) -----------
!                    2   2 sinh 2K(IS,IC)DEP(IX,IY)  |k(IS,IC)|
!
!
  USE M_GENARR
  USE SWCOMM3                                                         
  USE SWCOMM4                                                         
  USE OCPCOMM4
# if defined (EXPLICIT)
  USE MOD_ACTION_EX
# else         
  USE MOD_ACTION_IM
# endif    
  USE ALL_VARS 
   
  IMPLICIT NONE                                                     

  INTEGER :: IC,IS,ID,IG,NICMAX,INDX
  REAL    :: NN(1:MSC), ND(1:MSC)                                    
  REAL    :: DEP2(MT),KWAVE(MSC,NICMAX),CGO(MSC,NICMAX)
  REAL    :: DEPLOC

  DO IC = 1, NICMAX
    IF(IC == 1)THEN
      INDX  = IG
      DEPLOC = DEP2(INDX)
      IF(DEPLOC <= DEPMIN)THEN
!     *** depth is negative ***
        DO IS = 1, MSC
          KWAVE(IS,IC) = -1.                                           
          CGO(IS,IC)   = 0.                                            
        END DO	 
      ELSE
!     *** call KSCIP1 to compute KWAVE and CGO ***
        CALL KSCIP1(MSC,SPCSIG,DEPLOC,KWAVE(1,IC),CGO(1,IC),NN,ND)                                 
      ENDIF
    ELSE                                                    
      INDX  = NBVE(IG,IC-1)
      DEPLOC = DEP2(NV(INDX,1))+DEP2(NV(INDX,2))+DEP2(NV(INDX,3))
      DEPLOC = DEPLOC/3.0
      IF(DEPLOC <= DEPMIN)THEN
!     *** depth is negative ***
        DO IS = 1, MSC
          KWAVE(IS,IC) = -1.                                           
          CGO(IS,IC)   = 0.                                            
        END DO	 
      ELSE
!     *** call KSCIP1 to compute KWAVE and CGO ***
        CALL KSCIP1(MSC,SPCSIG,DEPLOC,KWAVE(1,IC),CGO(1,IC),NN,ND)
      ENDIF
    END IF  
  ENDDO                                                              

  RETURN
  END SUBROUTINE SWAPAR
 
 !
!
!****************************************************************
!
  SUBROUTINE SWAPAR1(I,IS,ID,DEP2,KWAVEL,CGOL)
!
!****************************************************************
!
!     computes the wave parameters K and CGO in the nearby
!     points, depending of the sweep direction.
!     The nearby points are indicated with the index IC (see
!     FUNCTION ICODE(_,_)
!
!  Method
!
!     The wave number K(IS,iC) is computed with the dispersion relation:
!
!     S = GRAV K(IS,IC)tanh(K(IS,IC)DEP(IX,IY))
!
!     where S = is logarithmic distributed via LOGSIG
!
!     The group velocity CGO in the case without current is equal to
!
!                    1       K(IS,IC)DEP(IX,IY)          S
!     CGO(IS,IC) = ( - + --------------------------) -----------
!                    2   2 sinh 2K(IS,IC)DEP(IX,IY)  |k(IS,IC)|
!
  USE M_GENARR
  USE SWCOMM3                                                         
  USE SWCOMM4                                                         
  USE OCPCOMM4 
# if defined (EXPLICIT)
  USE MOD_ACTION_EX
# else        
  USE MOD_ACTION_IM
# endif    
  USE ALL_VARS, ONLY : NV,NTVE,MT
   
  IMPLICIT NONE                                                     
!
  INTEGER :: IC,IS,ID,I,INDX
  REAL    :: NN(1:MSC), ND(1:MSC)                                    
  REAL    :: DEP2(MT),KWAVEL,CGOL
  REAL    :: DEPLOC,SPCSIGL

  DEPLOC = (DEP2(NV(I,1))+DEP2(NV(I,2))+DEP2(NV(I,3)))/3.0
  IF(DEPLOC <= DEPMIN)THEN
!   *** depth is negative ***
    KWAVEL = -1.                                           
    CGOL   = 0.                                            
  ELSE
!   *** call KSCIP1 to compute KWAVE and CGO ***
    SPCSIGL = SPCSIG(IS)
    CALL KSCIP1(1,SPCSIGL,DEPLOC,KWAVEL,CGOL,NN,ND)                                 
  END IF

  RETURN
  END SUBROUTINE SWAPAR1
!    
!****************************************************************
!
   SUBROUTINE SPROXY (I1     ,IS     ,ID     ,CAXL   ,CAYL   ,  &
                      CG0L   ,ECOSL  ,ESINL  ,UX2L   ,UY2L   )
!
!****************************************************************
!
!        computes the propagation velocities of energy in X-, Y-
!        -space, i.e., CAX, CAY, in the presence or absence of
!        currents, for the action balance equation.
!
!        The propagation velocities are computed for the fully 360
!        degrees sector.
!
!     METHOD
!
!        The next equation are calculated:
!
!              @X     _
!        CAX = -- = n C cos (id) + Ux  = CGO cos(id) + Ux
!              @T
!
!              @Y     _
!        CAY = -- = n C sin(id)  + Uy  = CGO sin(id) + Uy
!              @T
!                                                         _
!
!       ******************************************************************
!       *  attention! in the action balance equation the term            *
!       *  dx                                                            *
!       *  -- = CGO + U = CX  with x, CGO, U and CX vectors              *
!       *  dt                                                            *
!       *  is in the literature the term dx/dt often indicated           *
!       *  with CX and CY in the action balance equation.                *
!       *  In this program we use:    CAX = CGO + U                      *
!       ******************************************************************
!
!   ------------------------------------------------------------
!   If depth is negative ( DEP(IX,IY) <= 0), then,
!     For every point in S and D-direction do,
!       Give propagation velocities default values :
!       CAX(ID,IS,IC)     = 0.   {propagation velocity of energy in X-dir.}
!       CAY(ID,IS,IC)     = 0.   {propagation velocity of energy in Y-dir.}
!     ---------------------------------------------------------
!   Else if current is on (ICUR > 0) then,
!     For every point in S and D-direction do,  {using the output of SWAPAR}
!       S = logaritmic distributed via LOGSIG
!       Compute propagation velocity in X-direction:
!
!               1    K(IS,IC)DEP2(IX,IY)      S cos(D)
!       CAX = ( - + ------------------------) --------- + UX2(IX,IY)
!               2   sinh 2K(IS,IC)DEP2(IX,IY) |K(IS,IC)|
!
!       ------------------------------------------------------
!       Compute propagation velocity in Y-direction:
!
!               1    K(IS,IC)DEP2(IX,IY)      S sin(D)
!       CAY = ( - + ------------------------) -------- + UY2(IX,IY)
!               2   sinh 2K(IS,IC)DEP2(IX,IY) |K(IS,IC)|
!
!       ------------------------------------------------------
!   Else if current is not on (ICUR = 0)
!     For every point in S and D-direction do
!       S = logarithmic distributed via LOGSIG
!       Compute propagation velocity in X-direction:
!
!               1    K(IS,IC)DEP2(IX,IY)        S cos(D)
!       CAX = ( - + ------------------------) ----------
!               2   sinh 2K(IS,IC)DEP2(IX,IY)  |K(IS,IC)|
!
!       ------------------------------------------------------
!       Compute propagation velocity in Y-direction:
!
!               1    K(IS,IC)DEP2(IX,IY)        S sin(D)
!       CAY = ( - + ------------------------) ----------
!               2   sinh 2K(IS,IC)DEP2(IX,IY)  |K(IS,IC)|
!
!     ----------------------------------------------------------
!   End IF
!   ------------------------------------------------------------
!****************************************************************
   USE SWCOMM3                                                         
   USE SWCOMM4                                                         
   USE OCPCOMM4                                                        
   USE M_DIFFR 
#  if defined (EXPLICIT)
   USE MOD_ACTION_EX
#  else      
   USE MOD_ACTION_IM
#  endif   
   
   IMPLICIT NONE                                              

   INTEGER  IC,IS,ID,I1,NICMAX
   REAL     CAXL,CAYL,CG0L,ECOSL,ESINL,UX2L,UY2L

   CAXL = CG0L * ECOSL
   CAYL = CG0L * ESINL
!
!    --- adapt the velocities in case of diffraction
!
   IF(IDIFFR == 1 .AND. PDIFFR(3) /= 0.)THEN 
!     CAXL = CAXL*DIFPARAM(I1)      
!     CAYL = CAYL*DIFPARAM(I1)      
   END IF
!
!    --- ambient currents added
!
   IF(ICUR == 1)THEN 
     CAXL = CAXL + UX2L
     CAYL = CAYL + UY2L
   END IF

   RETURN
   END SUBROUTINE SPROXY
!
!
!****************************************************************
!
  SUBROUTINE SPROSD (KWAVE    ,CAS      ,CAD      ,          &
                     CGO      ,DEP2     ,DEP1     ,          &
		     ECOS     ,ESIN     ,UX2      ,          &
		     UY2      ,IDCMIN   ,IDCMAX   ,          &
		     COSCOS   ,SINSIN   ,SINCOS   ,          &
		     RDX      ,RDY      ,CAX      ,          &
		     CAY      ,IG       )
!
!****************************************************************
!
!     computes the propagation velocities of energy in S- and
!     D-space, i.e., CAS, CAD, in the presence or absence of
!     currents, for the action balance equation.
!
!  Method
!
!     The next equation are solved numerically
!
!           @S   @S   @D   _     @D   @D          _   @U
!     CAS = -- = -- [ -- + U . ( -- + --) ] - CGO K . --
!           @T   @D   @T         @X   @Y              @s
!
!           with:   @S       KS
!                   -- =  ---------
!                   @D    sinh(2KD)
!
!           @D      S      @D         @D           @Ux   @Uy
!     CAD = -- = ------- [ --sin(D) - --cos(D)] + [--- - ---] *
!           @T  sinh(2KD)  @X         @Y            @X   @Y
!
!                        @Uy               @Ux
!     * sin(D)cos(D) +   ---sin(D)sin(D) - ---cos(D)cos(D)
!                        @X                @Y
!   ------------------------------------------------------------
!   For current sweep and two adjacent sweeps do
!       determine interpolation factors RDXL and RDYL
!       determine depth and current gradients
!   ------------------------------------------------------------
!   For each frequency do
!       determine auxiliary quantities depending on sigma
!       For each direction in the sweep and two neighbouring
!           directions do
!           If IREFR=-1
!           Then compute reduction factor for contribution due
!                to depth gradient
!           ----------------------------------------------------
!           determine sweep in which this direction is located
!           using gradients of the proper sweep determine
!           Csigma (CAS) and Ctheta (CAD)
!   ------------------------------------------------------------
!   If ITFRE=0
!   Then make values of CAS=0
!   ------------------------------------------------------------
!   If IREFR=0
!   Then make values of CAD=0
!   ------------------------------------------------------------
!
!****************************************************************
!
  USE M_GENARR
  USE SWCOMM2                                           
  USE SWCOMM3                                                         
  USE SWCOMM4                                                         
  USE TIMECOMM                                                        
  USE OCPCOMM4                                                        
  USE M_DIFFR   
# if defined (EXPLICIT)
  USE MOD_ACTION_EX
# else            
  USE MOD_ACTION_IM
# endif      
  USE ALL_VARS, ONLY : NTVE,NBVE,NV,VX,VY,XC,YC,ART1,MT,ISONB_W,NBSN,NTSN                                                      
!
  IMPLICIT NONE
!
  INTEGER, INTENT(IN) :: IDCMIN(MSC), IDCMAX(MSC)
  REAL  :: CAS(MDC,MSC,MICMAX)                                        
  REAL  :: CAD(MDC,MSC,MICMAX)                                        
  REAL  :: CAX(MDC,MSC,MICMAX)                                        
  REAL  :: CAY(MDC,MSC,MICMAX)                                        
  REAL  :: CGO(MSC,MICMAX)                                            
  REAL  :: DEP2(MT) ,DEP1(MT) ,ECOS(MDC) ,ESIN(MDC) ,COSCOS(MDC) ,    &
	   SINSIN(MDC) ,SINCOS(MDC)
  REAL  :: KWAVE(MSC,MICMAX)                                          
  REAL  :: UX2(MT) ,UY2(MT) ,RDX(10) ,RDY(10)  
  INTEGER  IENT ,IS ,ID ,II ,SWPNGB ,IDDUM ,ID1 ,ID2 ,ISWEEP   
  INTEGER :: ISWP                                              
  INTEGER :: IC, KCGI                                             
  LOGICAL :: VALSWP                                            
  REAL    :: VLSINH ,KD1   ,COEF
  REAL    :: RDXL(2) ,RDYL(2) ,DET ,DX2 ,DY2 ,DX3 ,DY3
  REAL    :: DPDX   ,DPDY   ,DUXDX ,DUXDY ,DUYDX ,DUYDY
  REAL    :: CAST1    ,CAST2    ,CAST3(3) ,CAST4(3) ,               &
             CAST5    ,CAST6(3) ,CAST7(3) ,CAST8(3) ,CAST9(3) ,     &
	     CADT1    ,CADT2(3) ,CADT3(3) ,                         &
	     CADT4(3) ,CADT5(3) ,CADT6(3) ,CADT7(3)
  REAL    :: DLOC1, DLOC2, DLOC3

  INTEGER, PARAMETER :: SWP_ARRAY(1:3) = (/2,1,3/)

  REAL    :: SUMHDX,SUMHDY,SUMUXDX,SUMUXDY,SUMUYDX,SUMUYDY
  REAL    :: SUMUXHDY,SUMUYHDX
  INTEGER :: IG,I
  REAL    :: UX,UY
  REAL    :: X1,X2,X3,Y1,Y2,Y3,DX1,DY1
      
  CAST1 = 0.
  CAST2 = 0.
  CAST5 = 0.
  CADT1 = 0.
  DLOC1 = DEP2(IG)

  CAS(:,:,1) = 0.
  CAD(:,:,1) = 0.
!
! *** coefficients for CAS -> function of IX and IY only ***
!
  IF(NSTATC == 0 .OR. .NOT.DYNDEP)THEN                            
!   *** stationary calculation ***
    CAST2 = 0.
  ELSE
!   nonstationary depth, CAST2 is @D/@t
    CAST2 = (DLOC1-DEP1(IG))*RDTIM
  END IF
      
  SUMHDX = 0.0
  SUMHDY = 0.0
  IF(ICUR /= 0)THEN
    SUMUXDX = 0.0
    SUMUXDY = 0.0
    SUMUYDX = 0.0
    SUMUYDY = 0.0
    SUMUXHDY = 0.0
    SUMUYHDX = 0.0
  END IF
      
  DO I = 1,NTVE(IG)
    UX = UX2(NV(NBVE(IG,I),1))+UX2(NV(NBVE(IG,I),2))+UX2(NV(NBVE(IG,I),3))
    UX = UX/3.0     
    UY = UY2(NV(NBVE(IG,I),1))+UY2(NV(NBVE(IG,I),2))+UY2(NV(NBVE(IG,I),3))
    UY = UY/3.0     
	
    DLOC2 = DEP2(NV(NBVE(IG,I),1)) +         &
            DEP2(NV(NBVE(IG,I),2)) +         &
            DEP2(NV(NBVE(IG,I),3)) 
    DLOC2 = DLOC2/3.0

    IF(IREFR == -1)THEN                                          
!     limitation of depths in neighbouring grid points
      DLOC2 = MIN(DLOC2, PNUMS(17)*DLOC1)
    ELSE                                                           
!     no limitation                                                
      DLOC2 = DLOC2                                           
    END IF                                                         
 
    X1 = VX(IG)
    Y1 = VY(IG)
	
    IF(NV(NBVE(IG,I),1) == IG)THEN
      X2 = VX(NV(NBVE(IG,I),2))	
      Y2 = VY(NV(NBVE(IG,I),2))	
      X3 = VX(NV(NBVE(IG,I),3))	
      Y3 = VY(NV(NBVE(IG,I),3))	
    ELSE IF(NV(NBVE(IG,I),2) == IG)THEN
      X2 = VX(NV(NBVE(IG,I),3))	
      Y2 = VY(NV(NBVE(IG,I),3))	
      X3 = VX(NV(NBVE(IG,I),1))	
      Y3 = VY(NV(NBVE(IG,I),1))	
    ELSE  
      X2 = VX(NV(NBVE(IG,I),1))	
      Y2 = VY(NV(NBVE(IG,I),1))	
      X3 = VX(NV(NBVE(IG,I),2))	
      Y3 = VY(NV(NBVE(IG,I),2))	
    END IF
	
    DX1 = 0.5*(X1+X2)
    DY1 = 0.5*(Y1+Y2)
    DX1 = XC(NBVE(IG,I))-DX1
    DY1 = YC(NBVE(IG,I))-DY1
	
    DX2 = 0.5*(X1+X3)
    DY2 = 0.5*(Y1+Y3)
    DX2 = DX2-XC(NBVE(IG,I))
    DY2 = DY2-YC(NBVE(IG,I))
	
    DX = DX1+DX2
    DY = DY1+DY2
	
    SUMHDX = SUMHDX - DLOC2*DX
    SUMHDY = SUMHDY - DLOC2*DY
	
    IF(ICUR /= 0)THEN
      SUMUXDX = SUMUXDX - UX*DX
      SUMUXDY = SUMUXDY - UX*DY
      SUMUYDX = SUMUYDX - UY*DX
      SUMUYDY = SUMUYDY - UY*DY
      SUMUXHDY = SUMUXHDY - UX*DLOC2*DY
      SUMUYHDX = SUMUYHDX - UY*DLOC2*DX
    END IF
  END DO	    

  IF(ISONB_W(IG) == 1)THEN
    DX = VX(NBSN(IG,2))-VX(NBSN(IG,NTSN(IG)-1))
    DX = 0.5*DX
    DY = VY(NBSN(IG,2))-VY(NBSN(IG,NTSN(IG)-1))
    DY = 0.5*DY
	
    SUMHDX = SUMHDX - DEP2(IG)*DX
    SUMHDY = SUMHDY - DEP2(IG)*DY

    IF(ICUR /= 0)THEN
      SUMUXDX = SUMUXDX - UX2(IG)*DX
      SUMUXDY = SUMUXDY - UX2(IG)*DY
      SUMUYDX = SUMUYDX - UY2(IG)*DX
      SUMUYDY = SUMUYDY - UY2(IG)*DY
      SUMUXHDY = SUMUXHDY - UX2(IG)*DEP2(IG)*DY
      SUMUYHDX = SUMUYHDX - UY2(IG)*DEP2(IG)*DX
    END IF
  END IF	

  DO IS =1,MSC
    KD1 = KWAVE(IS,1)*DLOC1
    IF(KD1 > 30.0)KD1 = 30.
    VLSINH = SINH(2.*KD1)
    COEF   = SPCSIG(IS)/VLSINH
!
!   *** coefficients for CAS -> function of IS only ***
!
    CAST1 = KWAVE(IS,1)*COEF
    CAST5 = CGO(IS,1)*KWAVE(IS,1)
!
!   *** coefficients for CAD -> function of IS only ***
!
    CADT1 =  COEF
!
!   loop over spectral directions
!
    DO IDDUM = IDCMIN(IS)-1, IDCMAX(IS)+1  
      ID = MOD(IDDUM-1+MDC, MDC)+1
      IF(ICUR == 0)THEN
        CAS(ID,IS,1) = CAST1*CAST2
        CAD(ID,IS,1) = CADT1*(ESIN(ID)*SUMHDY+ECOS(ID)*SUMHDX)
        CAD(ID,IS,1) = CAD(ID,IS,1)/ART1(IG)
	    
!       --- adapt the velocity in case of diffraction                
        IF(IDIFFR == 1)THEN                                        
!          CAD(ID,IS,1) = DIFPARAM(IG)*CAD(ID,IS,1)          &
!                       - DIFPARDX(IG)*CGO(IS,1)*ESIN(ID)    &
!	                + DIFPARDY(IG)*CGO(IS,1)*ECOS(ID)          
          PRINT*,'NOT FINISH YET. SEE SPROSD 001'
          STOP
        END IF
      ELSE
        IF(IDIFFR == 0)THEN                                        
          CAS(ID,IS,1)= CAST1*(CAST2*ART1(IG)+                            &
	 	        SUMUXHDY-SUMUYHDX-DLOC1*SUMUXDY+DLOC1*SUMUYDX)-   &
		        CAST5 *                                           &
		       (COSCOS(ID)*SUMUXDY-SINCOS(ID)*(SUMUXDX-SUMUYDY)-  &
		        SINSIN(ID)*SUMUYDX)
          CAS(ID,IS,1)=CAS(ID,IS,1)/ART1(IG)		     

          CAD(ID,IS,1) = CADT1*(ESIN(ID)*SUMHDY+ECOS(ID)*SUMHDX)      +   &
		         SINCOS(ID)*(SUMUXDY+SUMUYDX) +                   &
		         SINSIN(ID)*SUMUYDY+COSCOS(ID)*SUMUXDX
          CAD(ID,IS,1) = CAD(ID,IS,1)/ART1(IG)		   
        ELSE IF(IDIFFR == 1)THEN                                    
!          CAS(ID,IS,1)= CAST1 * (CAST2 + CAST3(ISWEEP) + CAST4(ISWEEP)) - &
!		         DIFPARAM(IG)*CAST5 *                              &
!		        (COSCOS(ID) * CAST6(ISWEEP) +                      &
!		         SINCOS(ID) * (CAST7(ISWEEP) + CAST8(ISWEEP)) +    &
!		         SINSIN(ID) * CAST9(ISWEEP) )                  

!          CAD(ID,IS,1) = DIFPARAM(IG)*CADT1 * (ESIN(ID) * CADT2(ISWEEP) - &
!			  ECOS(ID) * CADT3(ISWEEP)) -                      &
!	    	          DIFPARDX(IG)*CGO(IS,1)*ESIN(ID) +                &
!			  DIFPARDY(IG)*CGO(IS,1)*ECOS(ID) +                &
!		          SINCOS(ID) * (CADT4(ISWEEP) - CADT5(ISWEEP)) +   &
!		          SINSIN(ID) *  CADT6(ISWEEP) -                    &
!		          COSCOS(ID) *  CADT7(ISWEEP)                         
          PRINT*,'NOT FINISH YET. SEE SPROSD 002'
	  STOP
        END IF                                                        
      ENDIF
    END DO    !IDDUM
  END DO      !IS	  

!
! *** for most cases CAS and CAD will be activated. Therefore ***
! *** for IREFR is set 0 (no refraction) or ITFRE = 0 (no     ***
! *** frequency shift) we have put the IF statement outside   ***
! *** the internal loop above                                 ***
!
  IF(ITFRE == 0)THEN
    CAS(:,:,1) = 0.0
  ENDIF
!
  IF(IREFR == 0)THEN
    CAD(:,:,1) = 0.0
  ENDIF
!
  RETURN
  END SUBROUTINE SPROSD
!****************************************************************
!
  SUBROUTINE DSPHER (CAD, CAX, CAY, IG, ECOS, ESIN)       
!
!****************************************************************
!
!     computes the propagation velocities of energy in Theta-
!     space, i.e., CAD, due to use of spherical coordinates
!
!  Method
!
!     References:
!     W. E. Rogers, J. M. Kaihatu, H. A. H. Petit, N. Booij and L. H. Holthuijsen,
!     "Multiple-scale Propagation in a Third-Generation Wind Wave Model"
!     in preparation
!
!             Cg Cos(theta) Tan(latitude)
!     CAD = - ---------------------------
!                    Rearth2
!
!     The group velocity CG in the direction of the wave propagation
!     in case with a current is equal to:
!
!                     1      K(IS,IC)DEP(IX,IY)        S
!     CG(ID,IS,IC)= ( - + -----------------------) --------- +
!                     2  sinh 2K(IS,IC)DEP(IX,IY)  |k(IS,IC)|
!
!                     + (UX2(IX,IY)cos(D) + UY2(IX,IY)sin(D))
!
!     which is equivalent with CAX*cos(D) + CAY*sin(D)
!
!****************************************************************
!
  USE SWCOMM2                                                         
  USE SWCOMM3                                                         
  USE SWCOMM4                                                         
  USE OCPCOMM4  
  USE ALL_VARS, ONLY : VY                                                      
!
  IMPLICIT NONE
!
  REAL  :: CAD(MDC,MSC,MICMAX)                                        
  REAL  :: CAX(MDC,MSC,MICMAX)                                        
  REAL  :: CAY(MDC,MSC,MICMAX)                                        
  REAL  :: ECOS(MDC)
  REAL  :: ESIN(MDC)
  INTEGER :: IS, ID, IG
  REAL     TANLAT, CTTMP
!
!************************************************************************
!
!
! *** TANLAT is Tan of Latitude
!
  TANLAT = TAN(DEGRAD*(VY(IG)+YOFFS))
!
  DO ID = 1, MDC
    CTTMP = ECOS(ID) * TANLAT / REARTH2
    DO IS = 1, MSC
      CAD(ID,IS,1) = CAD(ID,IS,1) -                &
	        (CAX(ID,IS,1)*ECOS(ID) + CAY(ID,IS,1)*ESIN(ID)) * CTTMP  
    END DO
  END DO

  RETURN
  END SUBROUTINE DSPHER
 
!
!*******************************************************************
!
      SUBROUTINE ADDDIS (DISSXY     ,LEAKXY     ,                  &
                         AC2        ,ANYBIN     ,                  &
			 DISC0      ,DISC1      ,                  &
			 LEAKC1     ,SPCSIG     )               
!    (This subroutine has not been tested yet)
!
!*******************************************************************
!
      USE SWCOMM3 
      USE ALL_VARS, ONLY : MT                                                
!
!
!   --|-----------------------------------------------------------|--
!     | Delft University of Technology                            |
!     | Faculty of Civil Engineering                              |
!     | Environmental Fluid Mechanics Section                     |
!     | P.O. Box 5048, 2600 GA  Delft, The Netherlands            |
!     |                                                           |
!     | Programmers: R.C. Ris, N. Booij,                          |
!     |              IJ.G. Haagsma, A.T.M.M. Kieftenburg,         |
!     |              M. Zijlema, E.E. Kriezi,                     |
!     |              R. Padilla-Hernandez, L.H. Holthuijsen       |
!   --|-----------------------------------------------------------|--
!
!
!     SWAN (Simulating WAves Nearshore); a third generation wave model
!     Copyright (C) 2004-2005  Delft University of Technology
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation; either version 2 of
!     the License, or (at your option) any later version.
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     A copy of the GNU General Public License is available at
!     http://www.gnu.org/copyleft/gpl.html#SEC3
!     or by writing to the Free Software Foundation, Inc.,
!     59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
!
!
!  0. Authors
!
!     30.72: IJsbrand Haagsma
!
!  1. Updates
!
!     20.53, Aug. 95: New subroutine
!     30.74, Nov. 97: Prepared for version with INCLUDE statements
!     30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN
!
!  2. Purpose
!
!     Adds dissipation and leak
!
!  3. Method
!
!     ---
!
!  4. Argument variables
!
!     SPCSIG: Relative frequencies in computational domain in sigma-space 30.72
!
      REAL    SPCSIG(MSC)                                              
!
!     IX          Counter of gridpoints in x-direction
!     IY          Counter of gridpoints in y-direction
!     MXC         Maximum counter of gridppoints in x-direction
!     MYC         Maximum counter of gridppoints in y-direction
!     MSC         Maximum counter of relative frequency
!     MDC         Maximum counter of directional distribution
!
!     REALS:
!     ---------
!
!     SP          Dummy variable
!     TEMP        Dummy variable
!
!     one and more dimensional arrays:
!     ---------------------------------
!     AC2       4D    Action density as function of D,S,X,Y and T
!
!  7. Common blocks used
!
!
!  8. Subroutines used
!
!     ---
!
!  9. Subroutines calling
!
!     SWOMPU
!
! 11. Remarks
!
!     DISSXY and LEAKXY are dissipation and leak integrated over the
!     spectrum for each point in the computational grid
!     DISSC0 and DISSC1 give the dissipation distributed over the
!     spectral space in one point of the computational grid
!
! 12. Structure
!
!     -------------------------------------------------------------
!     -------------------------------------------------------------
!
! 13. Source text
!
      REAL     DISSXY(MT)    ,LEAKXY(MT)      ,          &
               DISC0(MDC,MSC)   ,DISC1(MDC,MSC)     ,          &
	       LEAKC1(MDC,MSC)  ,AC2(MDC,MSC,0:MT)            
!
      LOGICAL  ANYBIN(MDC,MSC)
      INTEGER, SAVE :: IENT=0
      CALL STRACE (IENT, 'ADDDIS')
!
      DO 100 ISC = 1, MSC
        DSDD = DDIR * FRINTF * SPCSIG(ISC)**2
        DO 90 IDC = 1, MDC
          IF (ANYBIN(IDC,ISC)) THEN
            DISSXY(KCGRD(1)) = DISSXY(KCGRD(1)) + DSDD*(DISC0(IDC,ISC) +  &
	                       DISC1(IDC,ISC) * AC2(IDC,ISC,KCGRD(1)))
            LEAKXY(KCGRD(1)) = LEAKXY(KCGRD(1)) + DSDD *                  &
	                       LEAKC1(IDC,ISC) * AC2(IDC,ISC,KCGRD(1))
          ENDIF
  90    CONTINUE
 100  CONTINUE
      RETURN
      END SUBROUTINE ADDDIS
!
!****************************************************************
!
   SUBROUTINE SPREDT (AC2     ,CAX     ,CAY     ,IDCMIN  ,IDCMAX  ,    &
	              ISSTOP  ,RDX     ,RDY     )     
!  (This subroutine has not been tested yet)
!
!****************************************************************
!
!        to estimate the action density depending of the sweep
!        direction during the first iteration of a stationary
!        computation. The reason for this is that AC2 is zero
!        at first iteration and no initialisation is given in
!        case of stationarity (NSTATC=0). Action density should
!        be nonzero because of the computation of the source
!        terms. The estimate is based on solving the equation
!
!            dN       dN
!        CAX -- + CAY -- = 0
!            dx       dy
!
!        in an explicit manner.
!
!     METHOD
!
!
!          [RDX1*CAX + RDY1*CAY]*N(i-1,j) + [RDX2*CAX + RDY2*CAY]*N(i,j-1)
! N(i,j) = ---------------------------------------------------------------
!                      (RDX1+RDX2) * CAX  +  (RDY1+RDY2) * CAY
!
!   ------------------------------------------------------------
!   For every sweep direction do,
!     For every point in S and D direction in sweep direction do,
!       predict values for action density at new point from values
!       of neighbour gridpoints taking into account spectral propagation
!       direction (with currents !!) and the boundary conditions.
!       --------------------------------------------------------
!       If wave action AC2 is negative, then
!         Give wave action initial value 1.E-10
!     ---------------------------------------------------------
!****************************************************************
!
   USE SWCOMM3                                                         
   USE SWCOMM4                                                         
   USE OCPCOMM4 
   USE ALL_VARS, ONLY : MT                                                       

   IMPLICIT NONE
   
   INTEGER :: IS ,ID ,IDDUM ,ISSTOP                                       
   REAL    :: FAC_A ,FAC_B
   REAL    :: AC2(MDC,MSC,0:MT)                                         
   REAL    :: CAX(MDC,MSC,MICMAX)                                        
   REAL    :: CAY(MDC,MSC,MICMAX)                                        
   REAL    :: RDX(10),  RDY(10)                                          
   INTEGER :: IDCMIN(MSC) ,IDCMAX(MSC)
   REAL    :: CDEN ,CNUM ,WEIG1, WEIG2
!
   DO IS = 1, ISSTOP
     DO IDDUM = IDCMIN(IS), IDCMAX(IS)
       ID = MOD ( IDDUM - 1 + MDC , MDC ) + 1
!
!      *** Computation of weighting coefs WEIG1 AND WEIG2 ***
!
       CDEN = RDX(1) * CAX(ID,IS,1) + RDY(1) * CAY(ID,IS,1)
       CNUM = (RDX(1) + RDX(2)) * CAX(ID,IS,1)              &
	    + (RDY(1) + RDY(2)) * CAY(ID,IS,1)
       WEIG1 = CDEN/CNUM
       WEIG2 = 1. - WEIG1
!
       FAC_A = WEIG1 * AC2(ID,IS,KCGRD(2))                         
       FAC_B = WEIG2 * AC2(ID,IS,KCGRD(3))                         
!
       IF (ACUPDA) AC2(ID,IS,KCGRD(1)) = MAX ( 0. , (FAC_A + FAC_B))      
!
     END DO   !IDDUM	    
   END DO     !IS  
!
   RETURN
   END SUBROUTINE SPREDT
!
!******************************************************************
!
  SUBROUTINE SWPSEL(IDCMIN    ,IDCMAX    ,CAX       ,             &
                    CAY       ,ISCMIN    ,ISCMAX    ,             &
		    IDTOT     ,ISTOT     ,IDDLOW    ,             &
		    IDDTOP    ,ISSTOP    ,DEP2      ,             &
		    UX2       ,UY2       ,SPCDIR    )
!
!******************************************************************
!
!     compute the frequency dependent counters in directional space
!     in a situation with a current and without a current.
!     The counters are only computed for the gridpoint
!     considered. This means IC = 1 
!
!******************************************************************
  USE SWCOMM1                                                         
  USE SWCOMM2                                                         
  USE SWCOMM3                                                         
  USE SWCOMM4                                                         
  USE OCPCOMM4                                                        
  USE ALL_VARS, ONLY : MT                                                       

  IMPLICIT NONE
!
  REAL    :: SPCDIR(MDC,6)                                               
  INTEGER :: IS ,ID ,IDSUM ,IDCLOW ,IDCHGH ,IDTOT ,ISTOT ,    & 
             IDDLOW ,IDDTOP ,ISSLOW ,ISSTOP ,IENT ,IDDUM ,ISCLOW ,    &
	     ISCHGH ,IX ,IY ,IC         
  REAL    :: CAXMID ,CAYMID ,GROUP ,UABS ,THDIR    
  INTEGER :: IDCMIN(MSC) ,IDCMAX(MSC) ,ISCMIN(MDC) ,ISCMAX(MDC) ,SECTOR(MSC)
  REAL    :: CAX(MDC,MSC,MICMAX) ,CAY(MDC,MSC,MICMAX) ,DEP2(MT) ,     &
             UX2(MT) ,UY2(MT) ,RDX(10) ,RDY(10)    
  LOGICAL :: LOWEST ,LOWBIN ,HGHBIN
!
! *** initialize arrays in frequency direction ***
!
  DO ID = 1, MDC
    ISCMIN(ID) = 1
    ISCMAX(ID) = 1
  END DO  
!
! *** initialize array's in theta direction ***
!
  DO IS = 1, MSC
    IDCMIN(IS) = 1
    IDCMAX(IS) = MDC
    SECTOR(IS) = 1
  END DO  
!
!     *** calculate minimum and maximum counters in frequency ***
!     *** space if a current is present: ISCMIN and ISCMAX    ***
!
  IDDLOW =  9999
  IDDTOP = -9999
  DO IS = 1 , MSC
    IF(SECTOR(IS) > 0)THEN
      IDDLOW = MIN ( IDDLOW , IDCMIN(IS) )
      IDDTOP = MAX ( IDDTOP , IDCMAX(IS) )
    END IF
  END DO
!
!     *** Determine counters ***
!
  DO IDDUM = IDDLOW, IDDTOP
    ID = MOD ( IDDUM - 1 + MDC , MDC ) + 1
    LOWEST = .TRUE.
    DO IS = 1, MSC
      IF(LOWEST)THEN
        ISCLOW = IS
        LOWEST = .FALSE.
      END IF
      ISCHGH = IS
    END DO
!
!       *** set the minimum and maximum counters in arrays ***
!
    IF(.NOT.LOWEST)THEN
      ISCMIN(ID) = ISCLOW
      ISCMAX(ID) = ISCHGH
    ELSE
    END IF
!
  END DO   !IDDUM
!
!     *** calculate the maximum number of counters in both ***
!     *** directional space and frequency space            ***
!
  IF(IDDLOW /= 9999)THEN
    IDTOT = ( IDDTOP - IDDLOW ) + 1
    IF(ICUR == 1)THEN
      IF(IDTOT < 3)THEN
        IDDTOP = IDDTOP + 1
        IF(IDTOT == 1) IDDLOW = IDDLOW - 1
        IDTOT = 3
      END IF
    END IF
  ELSE
    IDTOT = 0
  END IF
!
! *** set variables ***
!
  IDTOT  =     1
  ISTOT  =     1
  ISSLOW =  9999                                                     
  ISSTOP = -9999                                                     

  DO IS = 1, MSC
    IDCLOW  = 0
    IDCHGH  = 0
    IDSUM   = 0
    DO ID = 1, MDC
      IDSUM = IDSUM + 1
      ISSLOW = MIN(IS,ISSLOW)                                   
      ISSTOP = MAX(IS,ISSTOP)                                   
    END DO
  END DO  
!
  IF(ISSLOW /= 9999)THEN
    ISSLOW = 1
!   minimal value of ISSTOP is 4 (or MSC if MSC<4)                    
    IF(ICUR > 0) ISSTOP = MAX(MIN(4,MSC),ISSTOP)                    
    ISTOT = ( ISSTOP - ISSLOW ) + 1
  ELSE
    ISTOT = 0
  END IF
!
!     *** check if IDTOT is less then MDC ***
!
  IF(IDTOT > MDC)THEN
    IDDLOW = 1
    IDDTOP = MDC
    IDTOT  = MDC
  END IF
!
!     *** check if the lowest frequency is not blocked !    ***
!     *** this can occur in real cases if the depth is very ***
!     *** small and the current velocity is large           ***
!     *** the propagation velocity Cg = sqrt (gd) < U       ***
!
  IF(ICUR == 1 .AND. FULCIR .AND.                      &
     ISSLOW /= 1 .AND. ISSLOW /= 9999)THEN                          
!        CALL MSGERR (2,'The lowest freqency is blocked')                  
7002 FORMAT (I4, 1X, I4, 1X, I2)                                       
    IC = 1
    GROUP = SQRT ( GRAV_W * DEP2(KCGRD(IC)) )
  END IF

  RETURN
  END SUBROUTINE SWPSEL
