swancom5.f90 File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	swapar (IG, NICMAX, DEP2, KWAVE, CGO)
subroutine	swapar1 (I, IS, ID, DEP2, KWAVEL, CGOL)
subroutine	sproxy (I1, IS, ID, CAXL, CAYL, CG0L, ECOSL, ESINL, UX2L, UY2L)
subroutine	sproxy2 (CAXL, CAYL, CG0L, ECOSL, ESINL, UX2L, UY2L)
subroutine	sproxy3 (CAXL, CAYLA, CAYLB, CG0L, ECOSL, ESINL, UX2L, UY2L, DLTYETMPP, DLTXETMPP, DLTXEA, DLTXEB)
subroutine	sprosd (KWAVE, CAS, CAD, CGO, DEP2, DEP1, ECOS, ESIN, UX2, UY2, IDCMIN, IDCMAX, COSCOS, SINSIN, SINCOS, RDX, RDY, CAX, CAY, IG)
subroutine	dspher (CAD, CAX, CAY, IG, ECOS, ESIN)
subroutine	adddis (DISSXY, LEAKXY, AC2, ANYBIN, DISC0, DISC1, LEAKC1, SPCSIG)
subroutine	spredt (AC2, CAX, CAY, IDCMIN, IDCMAX, ISSTOP, RDX, RDY)
subroutine	swpsel (IDCMIN, IDCMAX, CAX, CAY, ISCMIN, ISCMAX, IDTOT, ISTOT, IDDLOW, IDDTOP, ISSTOP, DEP2, UX2, UY2, SPCDIR)

Function/Subroutine Documentation

adddis()

subroutine adddis	(	real, dimension(mt)	DISSXY,
		real, dimension(mt)	LEAKXY,
		real, dimension(mdc,msc,0:mt)	AC2,
		logical, dimension(mdc,msc)	ANYBIN,
		real, dimension(mdc,msc)	DISC0,
		real, dimension(mdc,msc)	DISC1,
		real, dimension(mdc,msc)	LEAKC1,
		real, dimension(msc)	SPCSIG
	)

Definition at line 875 of file swancom5.f90.

!    (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

Here is the call graph for this function:

dspher()

subroutine dspher	(	real, dimension(mdc,msc,micmax)	CAD,
		real, dimension(mdc,msc,micmax)	CAX,
		real, dimension(mdc,msc,micmax)	CAY,
		integer	IG,
		real, dimension(mdc)	ECOS,
		real, dimension(mdc)	ESIN
	)

Definition at line 804 of file swancom5.f90.

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

spredt()

subroutine spredt	(	real, dimension(mdc,msc,0:mt)	AC2,
		real, dimension(mdc,msc,micmax)	CAX,
		real, dimension(mdc,msc,micmax)	CAY,
		integer, dimension(msc)	IDCMIN,
		integer, dimension(msc)	IDCMAX,
		integer	ISSTOP,
		real, dimension(10)	RDX,
		real, dimension(10)	RDY
	)

Definition at line 1007 of file swancom5.f90.

!  (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

sprosd()

subroutine sprosd	(	real, dimension(msc,micmax)	KWAVE,
		real, dimension(mdc,msc,micmax)	CAS,
		real, dimension(mdc,msc,micmax)	CAD,
		real, dimension(msc,micmax)	CGO,
		real, dimension(mt)	DEP2,
		real, dimension(mt)	DEP1,
		real, dimension(mdc)	ECOS,
		real, dimension(mdc)	ESIN,
		real, dimension(mt)	UX2,
		real, dimension(mt)	UY2,
		integer, dimension(msc), intent(in)	IDCMIN,
		integer, dimension(msc), intent(in)	IDCMAX,
		real, dimension(mdc)	COSCOS,
		real, dimension(mdc)	SINSIN,
		real, dimension(mdc)	SINCOS,
		real, dimension(10)	RDX,
		real, dimension(10)	RDY,
		real, dimension(mdc,msc,micmax)	CAX,
		real, dimension(mdc,msc,micmax)	CAY,
		integer	IG
	)

Definition at line 501 of file swancom5.f90.

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

sproxy()

subroutine sproxy	(	integer	I1,
		integer	IS,
		integer	ID,
		real	CAXL,
		real	CAYL,
		real	CG0L,
		real	ECOSL,
		real	ESINL,
		real	UX2L,
		real	UY2L
	)

Definition at line 146 of file swancom5.f90.

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

Here is the caller graph for this function:

sproxy2()

subroutine sproxy2	(	real, dimension(mdc,msc)	CAXL,
		real, dimension(mdc,msc)	CAYL,
		real, dimension(mdc,msc)	CG0L,
		real, dimension(mdc)	ECOSL,
		real, dimension(mdc)	ESINL,
		real	UX2L,
		real	UY2L
	)

Definition at line 260 of file swancom5.f90.

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

Here is the caller graph for this function:

sproxy3()

subroutine sproxy3	(	real, dimension(mdc,msc)	CAXL,
		real, dimension(mdc,msc)	CAYLA,
		real, dimension(mdc,msc)	CAYLB,
		real, dimension(mdc,msc)	CG0L,
		real, dimension(mdc)	ECOSL,
		real, dimension(mdc)	ESINL,
		real	UX2L,
		real	UY2L,
		real	DLTYETMPP,
		real	DLTXETMPP,
		real	DLTXEA,
		real	DLTXEB
	)

Definition at line 376 of file swancom5.f90.

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

swapar()

subroutine swapar	(	integer	IG,
		integer	NICMAX,
		real, dimension(mt)	DEP2,
		real, dimension(msc,nicmax)	KWAVE,
		real, dimension(msc,nicmax)	CGO
	)

Definition at line 16 of file swancom5.f90.

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

Here is the call graph for this function:

swapar1()

subroutine swapar1	(	integer	I,
		integer	IS,
		integer	ID,
		real, dimension(mt)	DEP2,
		real	KWAVEL,
		real	CGOL
	)

Definition at line 92 of file swancom5.f90.

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

Here is the call graph for this function:

Here is the caller graph for this function:

swpsel()

subroutine swpsel	(	integer, dimension(msc)	IDCMIN,
		integer, dimension(msc)	IDCMAX,
		real, dimension(mdc,msc,micmax)	CAX,
		real, dimension(mdc,msc,micmax)	CAY,
		integer, dimension(mdc)	ISCMIN,
		integer, dimension(mdc)	ISCMAX,
		integer	IDTOT,
		integer	ISTOT,
		integer	IDDLOW,
		integer	IDDTOP,
		integer	ISSTOP,
		real, dimension(mt)	DEP2,
		real, dimension(mt)	UX2,
		real, dimension(mt)	UY2,
		real, dimension(mdc,6)	SPCDIR
	)

Definition at line 1090 of file swancom5.f90.

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