swancom5.f90

Go to the documentation of this file.

!
!****************************************************************
!
  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
  USE mod_action_im
  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 
  USE mod_action_im
  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
   USE mod_action_im

   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 sproxy2 (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 
   USE mod_action_im
   
   IMPLICIT NONE                                              

   INTEGER  IC,IS,ID,I1,NICMAX
   REAL     CAXL(MDC,MSC),CAYL(MDC,MSC),CG0L(MDC,MSC),ECOSL(MDC),ESINL(MDC),UX2L,UY2L

   DO id=1,mdc
     caxl(id,:) = cg0l(id,:) * ecosl(id)
     cayl(id,:) = cg0l(id,:) * esinl(id)
   END DO

!
!    --- 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 sproxy2
!
!
!****************************************************************
!
   SUBROUTINE sproxy3 (CAXL   ,CAYLA ,CAYLB,  &          !yzhang_w3
              CG0L   ,ECOSL  ,ESINL  ,UX2L   ,UY2L, DLTYETMPP, DLTXETMPP ,DLTXEA , DLTXEB)
!
!****************************************************************
!
!        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
   USE mod_action_ex

   IMPLICIT NONE
   INTEGER  IC,IS,ID,I1,NICMAX
   REAL   CAXL(MDC,MSC),CAYLA(MDC,MSC),CAYLB(MDC,MSC),CG0L(MDC,MSC),ECOSL(MDC),ESINL(MDC),UX2L,UY2L
   REAL     DLTXETMPP,DLTXEA,DLTXEB,DLTYETMPP

   DO id=1,mdc
     caxl(id,:) = cg0l(id,:) * ecosl(id) * dltyetmpp
     cayla(id,:) = cg0l(id,:) * esinl(id) * dltxea
     caylb(id,:) = cg0l(id,:) * esinl(id) * dltxeb
   END DO

!
!    --- 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*dltyetmpp
     cayla = cayla + uy2l*dltxetmpp
     caylb = caylb + uy2l*dltxetmpp
   END IF

   RETURN
   END SUBROUTINE sproxy3
!
!
!****************************************************************

!****************************************************************
!
  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   
  USE mod_action_im
  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  
      IF(abs(dloc2 - dloc1) > 100.0)THEN
        dloc2 = dloc1
      ELSE                                                           
!       no limitation                                                
        dloc2 = dloc2
      END IF                                           
    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