swancom4.f90 File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	fac4ww (XIS, SNLC1, DAL1, DAL2, DAL3, SPCSIG, WWINT, WWAWG, WWSWG)
subroutine	swlta (IG, DEP2, CGO, SPCSIG, KWAVE, IMATRA, IMATDA, IDDLOW, IDDTOP, ISSTOP, IDCMIN, IDCMAX, HS, SMEBRK, PLTRI, URSELL)
subroutine	range4 (WWINT, IDDLOW, IDDTOP)
subroutine	filnl3 (IDCMIN, IDCMAX, IMATRA, IMATDA, MEMNL4, PLNL4S, ISSTOP, IGC)
subroutine	swsnl8 (WWINT, UE, SA1, SA2, SPCSIG, SNLC1, DAL1, DAL2, DAL3, SFNL, DEP2, KMESPC, MEMNL4, FACHFR)
subroutine	riam_slw (LMAX, N, N2, G, H, DQ, DQ2, DT, DT2, W, P, ACT, SNL, MINT)
subroutine	swsnl4 (WWINT, WWAWG, SPCSIG, SNLC1, DAL1, DAL2, DAL3, DEP2, KMESPC, MEMNL4, FACHFR, IDIA, ITER)
subroutine	swsnl3 (WWINT, WWAWG, UE, SA1, SA2, SPCSIG, SNLC1, DAL1, DAL2, DAL3, SFNL, DEP2, KMESPC, MEMNL4, FACHFR)
subroutine	swsnl2 (IDDLOW, IDDTOP, WWINT, WWAWG, UE, SA1, ISSTOP, SA2, SPCSIG, SNLC1, DAL1, DAL2, DAL3, SFNL, DEP2, KMESPC, IMATDA, IMATRA, FACHFR, PLNL4S, IDCMIN, IDCMAX, IG, CGO, WWSWG)
subroutine	swsnl1 (WWINT, WWAWG, WWSWG, IDCMIN, IDCMAX, UE, SA1, SA2, DA1C, DA1P, DA1M, DA2C, DA2P, DA2M, SPCSIG, SNLC1, KMESPC, FACHFR, ISSTOP, DAL1, DAL2, DAL3, SFNL, DSNL, DEP2, AC2, IMATDA, IMATRA, PLNL4S, PLNL4D, IDDLOW, IDDTOP)
real function	wave (D)
subroutine	cgcmp (G, AK, H, CG)
subroutine	preset (LMAX, N, IW3L, IW4, N2, G, H, W4, AK4, W, DQ, DQ2, DT, DT2, P)
subroutine	findw1w2 (LMAX, IW3, W, G, H, W1, W2, W3, WA, AK1, AK2, AKA, T1A, T2A, IND)
subroutine	fdw1 (G, H, AKA, T1A, WA, X1, X2, X, EPS, M)
function	funcw1 (G, H, X, AKA, T1A, WA)

Function/Subroutine Documentation

cgcmp()

subroutine cgcmp	(		G,
			AK,
			H,
			CG
	)

Definition at line 2739 of file swancom4.f90.

!  Calculates group velocity Cg based on depth and wavenumber
!  Includes a deep water limit for kd > 10.

!  G     : gravitational acceleration
!  AK    : wave number
!  H     : depth
!  CG    : group velocity
!  AKH   : depth x Wave number (kd)
!  RGK   : square root of gk
!  SECH2 : the square of the secant Hyperbolic of kd
!          = 1 - TANH(kd)**2

!  Calculation of group velocity Cg

   akh=ak*h
   IF(akh <= 10.)THEN                                                 

!    Shallow water:

!    Cg = (1/2) (1 + (2 kd) / Sinh (2 kd)) (w / k)
!       = (1/2) (1 + (kd SECH2) / Tanh (kd)) (w / k)
!       = (1/2) (1 + (kd SECH2) / Tanh (kd)) (Sqrt (gk Tanh(kd)) / k)
!       = (1/2) (Sqrt (gk) (Sqrt(Tanh (kd)) + kd SECH2 / Sqrt (Tanh (kd)) / k
!       = (1/2) (Sqrt (k) g (Tanh (kd) + kd SECH2) / (k Sqrt (g Tanh (kd)))
!       = (g Tanh(kd) + gkd SECH2) / (2 Sqrt(gk Tanh(kd))

     sech2=1.-tanh(akh)**2
     cg=(g*akh*sech2+g*tanh(akh))/(2.*sqrt(g*ak*tanh(akh)))

   ELSE                                                                

!    Deep water:

!    Cg = w / (2 k)
!       = g / (2 Sqrt (gk))

     rgk=sqrt(g*ak)
     cg=g/(2.*rgk)

   ENDIF                                                               
!
   RETURN

Here is the caller graph for this function:

fac4ww()

subroutine fac4ww	(	real	XIS,
		real	SNLC1,
		real	DAL1,
		real	DAL2,
		real	DAL3,
		real, dimension(msc)	SPCSIG,
		integer, dimension(*)	WWINT,
		real, dimension(*)	WWAWG,
		real, dimension(*)	WWSWG
	)

Definition at line 17 of file swancom4.f90.

!(This subroutine has not been tested yet)
!
!******************************************************************

!
   USE swcomm3                                                         
   USE swcomm4                                                         
   USE ocpcomm4                                                        
   USE m_snl4                                                          
!
!  Purpose :
!
!     Calculate interpolation constants for Snl.
!
   REAL    SPCSIG(MSC)                                                 
   INTEGER     MSC2  ,MSC1  ,IS    ,IDP   ,IDP1  ,                &
               IDM   ,IDM1  ,ISP   ,ISP1  ,ISM   ,ISM1  ,         &
           IDPP  ,IDMM  ,ISPP  ,ISMM  ,                       &
           ISLOW ,ISHGH ,ISCLW ,ISCHG ,IDLOW ,IDHGH ,         &
               MSCMAX,MDCMAX                     

   REAL        SNLC1 ,LAMM2 ,LAMP2 ,DELTH3,                       &
               AUX1  ,DELTH4,DAL1  ,DAL2  ,DAL3  ,CIDP  ,WIDP  ,  &
           WIDP1 ,CIDM  ,WIDM  ,WIDM1 ,XIS   ,XISLN ,WISP  ,  &
           WISP1 ,WISM  ,WISM1 ,AWG1  ,AWG2  ,AWG3  ,AWG4  ,  &
           AWG5  ,AWG6  ,AWG7  ,AWG8  ,SWG1  ,SWG2  ,SWG3  ,  &
           SWG4  ,SWG5  ,SWG6  ,SWG7  ,SWG8  ,FREQ  ,         &
           RADE                                            
!
   REAL       WWAWG(*),WWSWG(*)
!
   INTEGER    WWINT(*)
!

   IF(ALLOCATED(af11)) DEALLOCATE(af11)                               

!  *** Compute frequency indices                               ***
!  *** XIS is the relative increment of the relative frequency ***
!
   msc2   = int( float(msc) / 2.0 )
   msc1   = msc2 - 1

   xis    = spcsig(msc2) / spcsig(msc1)  

!
!  *** set values for the nonlinear four-wave interactions ***
!
   snlc1  = 1. / grav_w**4                                               
!
   lamm2  = (1.-pquad(1))**2                                           
   lamp2  = (1.+pquad(1))**2                                           
   delth3 = acos( (lamm2**2+4.-lamp2**2) / (4.*lamm2) )
   aux1   = sin(delth3)
   delth4 = asin(-aux1*lamm2/lamp2)
!
   dal1   = 1. / (1.+pquad(1))**4                                      
   dal2   = 1. / (1.-pquad(1))**4                                      
   dal3   = 2. * dal1 * dal2
!
!  *** Compute directional indices in sigma and theta space ***
!
   cidp   = abs(delth4/ddir)                                           
   idp   = int(cidp)
   idp1  = idp + 1
   widp   = cidp - real(idp)
   widp1  = 1.- widp
!
   cidm   = abs(delth3/ddir)                                           
   idm   = int(cidm)
   idm1  = idm + 1
   widm   = cidm - real(idm)
   widm1  = 1.- widm
   xisln  = log( xis )
!
   isp    = int( log(1.+pquad(1)) / xisln )                            
   isp1   = isp + 1
   wisp   = (1.+pquad(1) - xis**isp) / (xis**isp1 - xis**isp)          
   wisp1  = 1. - wisp
!
   ism    = int( log(1.-pquad(1)) / xisln )                            
   ism1   = ism - 1
   wism   = (xis**ism -(1.-pquad(1))) / (xis**ism - xis**ism1)         
   wism1  = 1. - wism
!
!  *** Range of calculations ***
!
   islow =  1  + ism1
   ishgh = msc + isp1 - ism1
   isclw =  1
   ischg = msc - ism1
   idlow = 1 - mdc - max(idm1,idp1)
   idhgh = mdc + mdc + max(idm1,idp1)
!
   msc4mi = islow
   msc4ma = ishgh
   mdc4mi = idlow
   mdc4ma = idhgh
   mscmax = msc4ma - msc4mi + 1
   mdcmax = mdc4ma - mdc4mi + 1
!
!  *** Interpolation weights ***
!
   awg1   = widp  * wisp
   awg2   = widp1 * wisp
   awg3   = widp  * wisp1
   awg4   = widp1 * wisp1
!
   awg5   = widm  * wism
   awg6   = widm1 * wism
   awg7   = widm  * wism1
   awg8   = widm1 * wism1
!
!  *** quadratic interpolation ***
!
   swg1   = awg1**2
   swg2   = awg2**2
   swg3   = awg3**2
   swg4   = awg4**2
!
   swg5   = awg5**2
   swg6   = awg6**2
   swg7   = awg7**2
   swg8   = awg8**2
!
!  --- determine discrete counters for piecewise                       
!      constant interpolation                                          
!
   IF(awg1 < awg2)THEN
     IF(awg2 < awg3)THEN
       IF(awg3 < awg4)THEN
         ispp=isp
         idpp=idp
       ELSE
         ispp=isp
         idpp=idp1
       END IF
     ELSE IF(awg2 < awg4)THEN
       ispp=isp
       idpp=idp
     ELSE
       ispp=isp1
       idpp=idp
     END IF
   ELSE IF(awg1 < awg3)THEN
     IF(awg3 < awg4)THEN
       ispp=isp
       idpp=idp
     ELSE
       ispp=isp
       idpp=idp1
     END IF
   ELSE IF(awg1 < awg4)THEN
     ispp=isp
     idpp=idp
   ELSE
     ispp=isp1
     idpp=idp1
   END IF
   IF(awg5 < awg6)THEN
     IF(awg6 < awg7)THEN
       IF(awg7 < awg8)THEN
         ismm=ism
         idmm=idm
       ELSE
         ismm=ism
         idmm=idm1
       END IF
     ELSE IF(awg6 < awg8)THEN
       ismm=ism
       idmm=idm
     ELSE
       ismm=ism1
       idmm=idm
     END IF
   ELSE IF(awg5 < awg7)THEN
     IF(awg7 < awg8)THEN
       ismm=ism
       idmm=idm
     ELSE
       ismm=ism
       idmm=idm1
     END IF
   ELSE IF(awg5 < awg8)THEN
     ismm=ism
     idmm=idm
   ELSE
     ismm=ism1
     idmm=idm1
   END IF
!
!  *** fill the arrays *
!
   wwint(1) = idp
   wwint(2) = idp1
   wwint(3) = idm
   wwint(4) = idm1
   wwint(5) = isp
   wwint(6) = isp1
   wwint(7) = ism
   wwint(8) = ism1
   wwint(9) = islow
   wwint(10)= ishgh
   wwint(11)= isclw
   wwint(12)= ischg
   wwint(13)= idlow
   wwint(14)= idhgh
   wwint(15)= msc4mi
   wwint(16)= msc4ma
   wwint(17)= mdc4mi
   wwint(18)= mdc4ma
   wwint(19)= mscmax
   wwint(20)= mdcmax
   wwint(21)= idpp                                                     
   wwint(22)= idmm                                                     
   wwint(23)= ispp                                                     
   wwint(24)= ismm                                                     
!
   wwawg(1) = awg1
   wwawg(2) = awg2
   wwawg(3) = awg3
   wwawg(4) = awg4
   wwawg(5) = awg5
   wwawg(6) = awg6
   wwawg(7) = awg7
   wwawg(8) = awg8
!
   wwswg(1) = swg1
   wwswg(2) = swg2
   wwswg(3) = swg3
   wwswg(4) = swg4
   wwswg(5) = swg5
   wwswg(6) = swg6
   wwswg(7) = swg7
   wwswg(8) = swg8

   ALLOCATE (af11(msc4mi:msc4ma))                                      

!  *** Fill scaling array (f**11)                     ***
!  *** compute the radian frequency**11 for IS=1, MSC ***
!
   DO is=1, msc
     af11(is) = ( spcsig(is) / ( 2. * pi_w ) )**11                       
   END DO  
!
!  *** compute the radian frequency for the IS = MSC+1, ISHGH ***
!
   freq   = spcsig(msc) / ( 2. * pi_w )                                  
   DO is = msc+1, ishgh
     freq   = freq * xis
     af11(is) = freq**11
   END DO  
!
!  *** compute the radian frequency for IS = 0, ISLOW ***
!
   freq   = spcsig(1) / ( 2. * pi_w )                                    
   DO is = 0, islow, -1
     freq   = freq / xis
     af11(is) = freq**11
   END DO  

   RETURN

fdw1()

subroutine fdw1	(		G,
			H,
			AKA,
			T1A,
			WA,
			X1,
			X2,
			X,
			EPS,
			M
	)

Definition at line 2881 of file swancom4.f90.

!
   m=0
   f1=funcw1(g,h,x1,aka,t1a,wa)
   f2=funcw1(g,h,x2,aka,t1a,wa)
20 IF(f1*f2) 18,19,21
18 m=m+1
   x=x2-((x2-x1)/(f2-f1)*f2)
   IF((abs(x-x2)/abs(x))-eps) 22,23,23
23 IF((abs(x-x1)/abs(x))-eps) 22,24,24
22 RETURN
24 f=funcw1(g,h,x,aka,t1a,wa)
   fm=f*f1
   IF(fm) 25,22,26
25 x2=x
   f2=f
   GO TO 20
26 x1=x
   f1=f
   GO TO 20
19 m=m+1
   IF(f1) 17,16,17
16 x=x1
   GO TO 22
17 x=x2
   GO TO 22
21 m=0
   RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

filnl3()

subroutine filnl3	(	integer, dimension(msc)	IDCMIN,
		integer, dimension(msc)	IDCMAX,
		real, dimension(mdc,msc)	IMATRA,
		real, dimension(mdc,msc)	IMATDA,
		real, dimension(mdc,msc,mt)	MEMNL4,
		real, dimension(mdc,msc,nptst)	PLNL4S,
		integer	ISSTOP,
		integer	IGC
	)

Definition at line 597 of file swancom4.f90.

!(This subroutine has not been tested yet)
!
!*******************************************************************
!
   USE swcomm3                                                         
   USE swcomm4                                                         
   USE ocpcomm4 
!   USE ALL_VARS, ONLY : MT,AC2                                                       
   USE vars_wave, ONLY : mt,ac2                                                       
   
   IMPLICIT NONE
!
!  2. Purpose
!
!     Fill the IMATRA/IMATDA arrays with the nonlinear wave-wave interaction
!     source term for a gridpoint ix,iy per sweep direction
!
!*******************************************************************
!
   INTEGER   IS,ID,ISSTOP,IDDUM,IGC
!
   REAL      IMATRA(MDC,MSC)           ,                          &
             IMATDA(MDC,MSC)           ,                          &
!        AC2(MDC,MSC,0:MT)        ,                          &
         plnl4s(mdc,msc,nptst)     ,                          &
         memnl4(mdc,msc,mt)                                
!
   INTEGER   IDCMIN(MSC),IDCMAX(MSC)
!
!   SAVE IENT
!   DATA IENT/0/
!   IF (LTRACE) CALL STRACE (IENT,'FILNL3')
!
   DO is=1, isstop
     DO iddum = idcmin(is), idcmax(is)
       id = mod( iddum - 1 + mdc , mdc ) + 1
       IF(testfl) plnl4s(id,is,iptst) = memnl4(id,is,igc)         
       IF(memnl4(id,is,igc) > 0.)THEN                          
          imatra(id,is) = imatra(id,is) + memnl4(id,is,igc)
       ELSE
          imatda(id,is) = imatda(id,is) - memnl4(id,is,igc) /     &
                      max(1.e-18,ac2(id,is,igc))
       END IF
     END DO
   END DO    
!
!   IF(TESTFL .AND. ITEST >= 50)THEN
!     WRITE(PRINTF,9000) IDCMIN(1),IDCMAX(1),MSC,ISSTOP
!9000 FORMAT(' FILNL3: ID_MIN ID_MAX MSC ISTOP :',4I6)
!     IF(ITEST >= 100)THEN
!       DO IS=1, ISSTOP
!         DO IDDUM = IDCMIN(IS), IDCMAX(IS)
!           ID = MOD ( IDDUM - 1 + MDC , MDC ) + 1
!           WRITE(PRINTF,6001) IS,ID,MEMNL4(ID,IS,KCGRD(1))             
!6001       FORMAT(' FILNL3: IS ID MEMNL()          :',2I6,E12.4)
!         ENDDO
!       ENDDO
!     ENDIF
!   ENDIF
!
   RETURN

findw1w2()

subroutine findw1w2	(		LMAX,
			IW3,
		real, dimension(lmax)	W,
			G,
			H,
			W1,
			W2,
			W3,
			WA,
			AK1,
			AK2,
			AKA,
			T1A,
			T2A,
			IND
	)

Definition at line 2853 of file swancom4.f90.

   REAL  :: W(LMAX)                                                    
!
   ind=0
   eps=0.0005
!
   x1=0.005
   x2=w3
!
!     ---------------------------------------
110 CALL fdw1(g,h,aka,t1a,wa,x1,x2,x,eps,m)
!     ---------------------------------------
!
   IF(m.EQ.0) GO TO 999
!
   w1=x
!
   ak1=wave(w1**2*h/g)/h
   ak2=sqrt(aka*aka+ak1*ak1-2.*aka*ak1*cos(t1a))
   w2=sqrt(g*ak2*tanh(ak2*h))
   IF(w1.GT.w2) GO TO 999
   t2a=atan2(-ak1*sin(t1a),aka-ak1*cos(t1a))
   RETURN
!
999 ind=-999
   RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

funcw1()

function funcw1	(		G,
			H,
			X,
			AKA,
			T1A,
			WA
	)

Definition at line 2911 of file swancom4.f90.

!
   ak1=wave(x**2*h/g)/h
   ak2=sqrt(aka*aka+ak1*ak1-2.*aka*ak1*cos(t1a))
   funcw1=wa-sqrt(g*ak1*tanh(ak1*h))-sqrt(g*ak2*tanh(ak2*h))
!
   RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

preset()

subroutine preset	(		LMAX,
			N,
			IW3L,
			IW4,
			N2,
			G,
			H,
			W4,
			AK4,
		real, dimension(lmax)	W,
			DQ,
			DQ2,
			DT,
			DT2,
			P
	)

Definition at line 2786 of file swancom4.f90.

!(This subroutine has not been tested yet)

!  T1A   : Angle between theta_1 and theta_a
!  T2A   : Angle between theta_2 and theta_a
!  T34   : Angle between theta_3 and theta_4
!  WA    : wa = w1 + w2 = w3 + w4 (resonance conditions)
!  AKA   : absolute value of ka = k1 + k2 = k3 + k4 (resonance conditions)
!  TA    : Angle theta_a representing vector ka

   REAL :: W(LMAX)                                                     
!
!  ================
   DO 120 it34=1,n2
!  ================
!
     t34=dt*(it34-1)
     dt3=dt
     IF(it34 == 1 .OR. it34 == n2) dt3=dt2
!
!    ===================
     DO 110 iw3=iw3l,iw4
!    ===================
!
       w3=w(iw3)
       ak3=wave(w3**2*h/g)/h
!
       wa=w3+w4
       aka=sqrt(ak3*ak3+ak4*ak4+2.*ak3*ak4*cos(t34))
       ta=atan2(ak3*sin(t34),ak3*cos(t34)+ak4)
       r=sqrt(g*aka*tanh(aka*h/2.))/wa-1./sqrt(2.)
!
       tl=0.
       acs=aka/(2.*wave(wa**2*h/(4.*g))/h)
       acs=sign(1.,acs)*min(1.,abs(acs))                                   
       IF(r < 0.) tl=acos(acs)
       it1s=nint(tl/dt)+1
!
!     ===================
       DO 100 it1a=it1s,n2
!     ===================
!
         t1a=dt*(it1a-1)
!
         dt1=dt
         IF(it1a == 1 .OR. it1a == n2) dt1=dt2
!
!     -----------------------------------------------------------------
         CALL findw1w2(lmax,iw3,w,g,h,w1,w2,w3,wa,ak1,ak2,aka,t1a,t2a,ind)   
!     -----------------------------------------------------------------
         IF(ind < 0) GO TO 100
!
         IF(w1 <= w3 .AND. w3 <= w4 .AND. w4 <= w2)THEN
!JQI        CALL KERNEL(N,G,H,W1,W2,W3,W4,AK1,AK2,AK3,AK4,AKA,           &
!JQI                T1A,T2A,T34,TA,DQ,DT,DT1,DT3,P)                       
         ENDIF
!
! ============
  100 CONTINUE
  110 CONTINUE
  120 CONTINUE
! ============
!
   RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

range4()

subroutine range4	(	integer, dimension(*)	WWINT,
		integer	IDDLOW,
		integer	IDDTOP
	)

Definition at line 507 of file swancom4.f90.

!
!******************************************************************
!
!     calculate the minimum and maximum counters in frequency and
!     directional space which fall with the calculation for the
!     nonlinear wave-wave interactions.
!
!     Method :  review for the counters :
!
!                            Frequencies -->
!                 +---+---------------------+---------+- IDHGH
!              d  | 3 :          2          :    2    |
!              i  + - + - - - - - - - - - - + - - - - +- MDC
!              r  |   :                     :         |
!              e  | 3 :  original spectrum  :    1    |
!              c  |   :                     :         |
!              t. + - + - - - - - - - - - - + - - - - +- 1
!                 | 3 :          2          :    2    |
!                 +---+---------------------+---------+- IDLOW
!                 |   |                     |    ^    |
!             ISLOW   1                     MSC  |  ISHGH
!                     ^                          |
!                     |                          |
!                    ISCLW                     ISCHG
!              lowest discrete               highest discrete
!                central bin                   central bin
!
!
!       The directional counters depend on the numerical method that
!       is used.
!****************************************************************
!
   USE swcomm3                                                         
   USE swcomm4                                                         
   USE ocpcomm4 
   IMPLICIT NONE

   INTEGER     IDDLOW,IDDTOP                                           
!
   INTEGER     WWINT(*)
!
!   *** Range in directional domain ***
!
   IF(iquad < 3 .AND. iquad > 0)THEN                         
!   *** counters based on bins which fall within a sweep ***
     wwint(13) = iddlow - max( wwint(4), wwint(2) )
     wwint(14) = iddtop + max( wwint(4), wwint(2) )
   ELSE
!   *** counters initially based on full circle ***
     wwint(13) = 1   - max( wwint(4), wwint(2) )
     wwint(14) = mdc + max( wwint(4), wwint(2) )
   END IF
!
!  *** error message ***
!
!   IF(WWINT(9) < WWINT(15) .OR. WWINT(10) > WWINT(16) .OR.          &
!      WWINT(13) < WWINT(17) .OR. WWINT(14) > WWINT(18))THEN
!     WRITE (PRINTF,900) IXCGRD(1), IYCGRD(1),                       &
!                        WWINT(9) ,WWINT(10) ,WWINT(13) ,WWINT(14),  &
!           WWINT(15),WWINT(16) ,WWINT(17) ,WWINT(18)
!900  FORMAT ( ' ** Error : array bounds and maxima in subr RANGE4, ', &
!              ' point ', 2I5,                                         &
!         /,'            ISL,ISH : ',2I4, '   IDL,IDH : ',2I4,    &
!         /,'            SMI,SMA : ',2I4, '   DMI,DMA : ',2I4)
!     IF(ITEST > 50) WRITE (PRTEST, 901) MSC, MDC, IDDLOW, IDDTOP
!901  FORMAT (' MSC, MDC, IDDLOW, IDDTOP: ', 4I5)
!   ENDIF
!
!  test output
!
!   IF(TESTFL .AND. ITEST >= 60)THEN
!     WRITE(PRTEST,911) WWINT(4), WWINT(2), WWINT(8), WWINT(6)
!911  FORMAT (' RANGE4: IDM1 IDP1 ISM1 ISP1    :',4I5)
!     WRITE(PRTEST,916) WWINT(11), WWINT(12), IQUAD
!916  FORMAT (' RANGE4: ISCLW ISCHG IQUAD      :',3I5)
!     WRITE (PRTEST,917) WWINT(9), WWINT(10), WWINT(13), WWINT(14)
!917  FORMAT (' RANGE4: ISLOW ISHGH IDLOW IDHGH:',4I5)
!     WRITE (PRTEST,919) WWINT(15), WWINT(16), WWINT(17), WWINT(18)
!919  FORMAT (' RANGE4: MS4MI MS4MA MD4MI MD4MA:',4I5)
!     WRITE(PRINTF,*)
!   END IF
!
   RETURN

riam_slw()

subroutine riam_slw	(		LMAX,
			N,
			N2,
			G,
			H,
			DQ,
			DQ2,
			DT,
			DT2,
		real, dimension(lmax)	W,
			P,
		real, dimension(lmax,n)	ACT,
		real, dimension(lmax,n)	SNL,
			MINT
	)

Definition at line 896 of file swancom4.f90.

!     (This subroutine has not been tested yet)

      USE m_snl4                                                          

!     LMAX
!     N     : number of directional bins                                  
!     N2
!     G     : gravitational acceleration                                  
!     H     : depth                                                       
!     DQ
!     DQ2
!     DT    : size of the directional bins (Delta Theta)                  
!     DT2
!     W     : discretised frequency array                                 
!     P     : density of water                                            
!     ACT   : action density                                              
!     SNL   : Quadruplet source term                                      
!     MINT

!     IW4   : counter for the 4th frequency                               
!     W4    : frequency of the fourth quadruplet component                
!     AK4   : wavenumber of the fourth quadruplet component               
!     DNA   : coefficient in eq. 17 of Hashimoto (98)                     
!             = 2 w4 k4 / Cg(k4)                                          
!     IW3L
!     II
!     JJ
!     DI
!     DJ
!     CGK4  : group velocity for the fourth quadruplet component          

      REAL :: W(LMAX), ACT(LMAX,N), SNL(LMAX,N)                           

!     Initialisation of the quadruplet source term                        

      snl(:,:) = 0.                                                       
!
!     =================
      DO 130 iw4=1,lmax
!     =================
!
      w4=w(iw4)

!     WAVE converts nondimensional Sqr(w)d/g to nondimensional kd

      ak4=wave(w4**2*h/g)/h

!     CGCMP computes group velocity

      CALL cgcmp(g,ak4,h,cgk4)

!     Calculates the coefficient in equation (17) of Hashimoto (1998)

      dna=2.*w4*ak4/cgk4

      iw3l=max0(1,nint(iw4-alog(3.)/dq))

!     ------------------------------------------------------------
      CALL preset(lmax,n,iw3l,iw4,n2,g,h,w4,ak4,w,dq,dq2,dt,dt2,p)        
!     ------------------------------------------------------------
!
!     ==============
      DO 130 it4=1,n
!     ==============
!
!     ================
      curriam => friam                                                    
      DO                !    (W1 - W3 - W4 - W2)                          
!     ================
!
      k1=curriam%II(1)+iw4                                                
      k2=curriam%II(2)+iw4                                                
      k3=curriam%II(3)+iw4                                                
!
      IF(k1.LT.1   ) GO TO 120                                            
      IF(k2.GT.lmax) GO TO 120                                            
!
      gg=curriam%SSS*dna                                                  
!
      m1p= curriam%JJ(1)+it4                                              
      m1n=-curriam%JJ(1)+it4                                              
      m2p= curriam%JJ(2)+it4                                              
      m2n=-curriam%JJ(2)+it4                                              
      m3p= curriam%JJ(3)+it4                                              
      m3n=-curriam%JJ(3)+it4                                              
!
      m1p=mod(m1p,n)
      m1n=mod(m1n,n)
      m2p=mod(m2p,n)
      m2n=mod(m2n,n)
      m3p=mod(m3p,n)
      m3n=mod(m3n,n)
!
      IF(m1p.LT.1) m1p=m1p+n
      IF(m1n.LT.1) m1n=m1n+n
      IF(m2p.LT.1) m2p=m2p+n
      IF(m2n.LT.1) m2n=m2n+n
      IF(m3p.LT.1) m3p=m3p+n
      IF(m3n.LT.1) m3n=m3n+n
!
      a4=act(iw4,it4)
!
      IF(mint.EQ.1) THEN
!
        k01=0
        k02=0
        k03=0
        m01=0
        m02=0
        m03=0
        IF(curriam%DI(1).LE.1.) k01= 1                                    
        IF(curriam%DI(2).LE.1.) k02= 1                                    
        IF(curriam%DI(3).LE.1.) k03= 1                                    
        IF(curriam%DJ(1).LE.1.) m01=-1                                    
        IF(curriam%DJ(2).LE.1.) m02=-1                                    
        IF(curriam%DJ(3).LE.1.) m03=-1                                    
!
        k11=k1+k01
        k21=k2+k02
        k31=k3+k03
        IF(k11.GT.lmax) k11=lmax
        IF(k21.GT.lmax) k21=lmax
        IF(k31.GT.lmax) k31=lmax
!
        k111=k1+k01-1
        k211=k2+k02-1
        k311=k3+k03-1
        IF(k111.LT.1) k111=1
        IF(k211.LT.1) k211=1
        IF(k311.LT.1) k311=1
!
        m1p1=m1p+m01
        m2p1=m2p+m02
        m3p1=m3p+m03
        IF(m1p1.LT.1) m1p1=n
        IF(m2p1.LT.1) m2p1=n
        IF(m3p1.LT.1) m3p1=n
!
        m1n1=m1n-m01
        m2n1=m2n-m02
        m3n1=m3n-m03
        IF(m1n1.GT.n) m1n1=1
        IF(m2n1.GT.n) m2n1=1
        IF(m3n1.GT.n) m3n1=1
!
        m1p11=m1p+m01+1
        m2p11=m2p+m02+1
        m3p11=m3p+m03+1
        IF(m1p11.GT.n) m1p11=1
        IF(m2p11.GT.n) m2p11=1
        IF(m3p11.GT.n) m3p11=1
!
        m1n11=m1n-m01-1
        m2n11=m2n-m02-1
        m3n11=m3n-m03-1
        IF(m1n11.LT.1) m1n11=n
        IF(m2n11.LT.1) m2n11=n
        IF(m3n11.LT.1) m3n11=n
!
        di1=curriam%DI(1)                                                 
        di2=curriam%DI(2)                                                 
        di3=curriam%DI(3)                                                 
        dj1=curriam%DJ(1)                                                 
        dj2=curriam%DJ(2)                                                 
        dj3=curriam%DJ(3)                                                 
!
        a1p=(di1*(dj1*act(k11, m1p1)+act(k11, m1p11))                 &
            + dj1*act(k111,m1p1)+act(k111,m1p11))                 &
        /((1.+di1)*(1.+dj1))                                      
!
        a1n=(di1*(dj1*act(k11, m1n1)+act(k11, m1n11))                 &
            + dj1*act(k111,m1n1)+act(k111,m1n11))                 &
        /((1.+di1)*(1.+dj1))                                  
!
        a2p=(di2*(dj2*act(k21, m2p1)+act(k21, m2p11))                 &
            + dj2*act(k211,m2p1)+act(k211,m2p11))                 &
        /((1.+di2)*(1.+dj2))                                  
!
        a2n=(di2*(dj2*act(k21, m2n1)+act(k21, m2n11))                 &
            + dj2*act(k211,m2n1)+act(k211,m2n11))                 &
        /((1.+di2)*(1.+dj2))                                
!
        a3p=(di3*(dj3*act(k31, m3p1)+act(k31, m3p11))                 &
            + dj3*act(k311,m3p1)+act(k311,m3p11))                 &
        /((1.+di3)*(1.+dj3))                                 
!
        a3n=(di3*(dj3*act(k31, m3n1)+act(k31, m3n11))                 &
            + dj3*act(k311,m3n1)+act(k311,m3n11))                 &
        /((1.+di3)*(1.+dj3))                                 
!
      ELSE
!
        IF(k1.LT.1) k1=1
        IF(k2.LT.1) k2=1
        IF(k3.LT.1) k3=1
!
        a1p=act(k1,m1p)
        a1n=act(k1,m1n)
        a2p=act(k2,m2p)
        a2n=act(k2,m2n)
        a3p=act(k3,m3p)
        a3n=act(k3,m3n)
!
      ENDIF
!
      w1p2p=a1p+a2p
      w1n2n=a1n+a2n
      s1p2p=a1p*a2p
      s1n2n=a1n*a2n
      w3p4 =a3p+a4
      w3n4 =a3n+a4
      s3p4 =a3p*a4
      s3n4 =a3n*a4
!
      xp=s1p2p*w3p4-s3p4*w1p2p
      xn=s1n2n*w3n4-s3n4*w1n2n
!
      snl( k1,m1p)=snl( k1,m1p)-xp*gg
      snl( k1,m1n)=snl( k1,m1n)-xn*gg
      snl( k2,m2p)=snl( k2,m2p)-xp*gg
      snl( k2,m2n)=snl( k2,m2n)-xn*gg
      snl( k3,m3p)=snl( k3,m3p)+xp*gg
      snl( k3,m3n)=snl( k3,m3n)+xn*gg
      snl(iw4,it4)=snl(iw4,it4)+(xp+xn)*gg
!
! ============
  120 CONTINUE
      IF (.NOT.ASSOCIATED(curriam%NEXTRIAM)) EXIT                         
      curriam => curriam%NEXTRIAM                                         
      END DO                                                              
! ============
!
! ============
  130 CONTINUE
! ============
!
      RETURN

Here is the call graph for this function:

swlta()

subroutine swlta	(	integer	IG,
		real, dimension(mt)	DEP2,
		real, dimension(msc,micmax)	CGO,
		real, dimension(msc)	SPCSIG,
		real, dimension(msc,micmax)	KWAVE,
		real, dimension(mdc,msc)	IMATRA,
		real, dimension(mdc,msc)	IMATDA,
		integer	IDDLOW,
		integer	IDDTOP,
		integer	ISSTOP,
		integer, dimension(msc)	IDCMIN,
		integer, dimension(msc)	IDCMAX,
		real	HS,
		real	SMEBRK,
		real, dimension(mdc,msc,nptst)	PLTRI,
		real, dimension(mt)	URSELL
	)

Definition at line 287 of file swancom4.f90.

! (This subroutine has not been tested yet)           
!
!****************************************************************
!
   USE ocpcomm4
   USE swcomm3
   USE swcomm4
   USE vars_wave, ONLY : mt, ac2
!
   IMPLICIT NONE
!
!  2. Purpose
!
!     In this subroutine the triad-wave interactions are calculated
!     with the Lumped Triad Approximation of Eldeberky (1996). His
!     expression is based on a parametrization of the biphase (as
!     function of the Ursell number), is directionally uncoupled and
!     takes into account for the self-self interactions only.
!
!     For a full description of the equations reference is made
!     to PhD thesis of Eldeberky (1996). Here only the main expressions
!     are given.
!
!  3. Method
!
!     The parametrized biphase is given by (see eq. 3.19):
!
!                                  0.2
!     beta = - pi/2 + pi/2 tanh ( ----- )
!                                   Ur
!
!     The Ursell number is calculated in routine SINTGRL
!
!     The source term as function of frequency p is (see eq. 7.25):
!
!             +      -
!     S(p) = S(p) + S(p)
!
!     in which
!
!      +
!     S(p) = alpha Cp Cg,p (R(p/2,p/2))**2 sin (|beta|) ( E(p/2)**2 -2 E(p) E(p/2) )
!
!      -          +
!     S(p) = - 2 S(2p)
!
!     with alpha a tunable coefficient and R(p/2,p/2) is the interaction
!     coefficient of which the expression can be found in Eldeberky (1996);
!     see eq. 7.26.
!
!     Note that a slightly adapted formulation of the LTA is used in
!     in the SWAN model:
!
!     - Only positive contributions to higher harmonics are considered
!       here (no energy is transferred to lower harmonics).
!
!     - The mean frequency in the expression of the Ursell number
!       is calculated according to the first order moment over the
!       zeroth order moment (personal communication, Y.Eldeberky, 1997).
!
!     - The interactions are calculated up to 2.5 times the mean
!       frequency only.
!
!     - Since the spectral grid is logarithmically distributed in frequency
!       space, the interactions between central bin and interacting bin
!       are interpolated such that the distance between these bins is
!       factor 2 (nearly).
!
!     - The interactions are calculated in terms of energy density
!       instead of action density. So the action density spectrum
!       is firstly converted to the energy density grid, then the
!       interactions are calculated and then the spectrum is converted
!       to the action density spectrum back.
!
!     - To ensure numerical stability the Patankar rule is used.
!
   INTEGER IDDLOW, IDDTOP, ISSTOP
   INTEGER IDCMIN(MSC), IDCMAX(MSC)

   REAL :: HS, SMEBRK
!  REAL :: AC2(MDC,MSC,0:MT)

   REAL :: CGO(MSC,MICMAX)
   REAL :: DEP2(MT)
   REAL :: IMATDA(MDC,MSC), IMATRA(MDC,MSC)
   REAL :: SPCSIG(MSC)
   REAL :: KWAVE(MSC,MICMAX)
   REAL :: PLTRI(MDC,MSC,NPTST)
   REAL :: URSELL(MT)

   INTEGER I1, I2, ID, IG, IDDUM, IENT, II, IS, ISM, ISM1, ISMAX, ISP, ISP1
   REAL    AUX1, AUX2, BIPH, C0, CM, DEP, DEP_2, DEP_3, E0, EM,     &
           FT, RINT, SIGPI, SINBPH, STRI, WISM, WISM1, WISP, WISP1, &
       W0, WM, WN0, WNM, XIS, XISLN
   REAL, ALLOCATABLE :: E(:), SA(:,:)

!   DEP   = DEP2(IGC)
   dep   = dep2(ig)
   dep_2 = dep**2
   dep_3 = dep**3
!
!  --- compute some indices in sigma space
!
   i2     = int(float(msc) / 2.)
   i1     = i2 - 1
   xis    = spcsig(i2) / spcsig(i1)
   xisln  = log( xis )

   isp    = int( log(2.) / xisln )
   isp1   = isp + 1
   wisp   = (2. - xis**isp) / (xis**isp1 - xis**isp)
   wisp1  = 1. - wisp

   ism    = int( log(0.5) / xisln )
   ism1   = ism - 1
   wism   = (xis**ism -0.5) / (xis**ism - xis**ism1)
   wism1  = 1. - wism

   ALLOCATE (e(1:msc))
   ALLOCATE (sa(1:mdc,1:msc+isp1))
   e  = 0.
   sa = 0.
!
!  --- compute maximum frequency for which interactions are calculated
!
   ismax = 1
   DO is = 1, msc
     IF(spcsig(is) < ( ptriad(2) * smebrk) )THEN
       ismax = is
     ENDIF
   ENDDO
   ismax = max( ismax , isp1 )
!
!  --- compute 3 wave-wave interactions
!
!   IF( URSELL(IGC) >= PTRIAD(5) )THEN
   IF( ursell(ig) >= ptriad(5) )THEN
!
!    --- calculate biphase
!
!     BIPH   = (0.5*PI_W)*(TANH(PTRIAD(4)/URSELL(IGC))-1.)
     biph   = (0.5*pi_w)*(tanh(ptriad(4)/ursell(ig))-1.)
     sinbph = abs( sin(biph) )
!
     DO ii = iddlow, iddtop
       id = mod( ii - 1 + mdc , mdc ) + 1
!
!      --- initialize array with E(f) for the direction considered
!
       DO is = 1, msc
!         E(IS) = AC2(ID,IS,IGC) * 2. * PI_W * SPCSIG(IS)
         e(is) = ac2(id,is,ig) * 2. * pi_w * spcsig(is)
       END DO
!
       DO is = 1, ismax

         e0  = e(is)
         w0  = spcsig(is)
         wn0 = kwave(is,1)
         c0  = w0 / wn0

         IF( is > -ism1 )THEN
           em  = wism * e(is+ism1)       + wism1 * e(is+ism)
           wm  = wism * spcsig(is+ism1)  + wism1 * spcsig(is+ism)
           wnm = wism * kwave(is+ism1,1) + wism1 * kwave(is+ism,1)
           cm  = wm / wnm
         ELSE
           em  = 0.
           wm  = 0.
           wnm = 0.
           cm  = 0.
         END IF

         aux1 = wnm**2 * ( grav_w * dep + 2.*cm**2 )
         aux2 = wn0 * dep * ( grav_w * dep +                            &
            (2./15.) * grav_w * dep_3 * wn0**2 -                    &
        (2./ 5.) * w0**2 * dep_2 )
         rint = aux1 / aux2
         ft = ptriad(1) * c0 * cgo(is,1) * rint**2 * sinbph

         sa(id,is) = max(0., ft * ( em * em - 2. * em * e0 ))

       END DO
     END DO
!
!    ---  put source term together
!
     DO is = 1, isstop
       sigpi = spcsig(is) * 2. * pi_w
       DO iddum = idcmin(is), idcmax(is)
         id = mod( iddum - 1 + mdc , mdc ) + 1
!
         stri = sa(id,is) - 2.*(wisp  * sa(id,is+isp1) +     &
            wisp1 * sa(id,is+isp ))
!
!        --- store results in rhs and main diagonal according
!            to Patankar-rules
!
         IF(testfl) pltri(id,is,iptst) = stri / sigpi
         IF(stri > 0.)THEN
           imatra(id,is) = imatra(id,is) + stri / sigpi
         ELSE
           imatda(id,is) = imatda(id,is) - stri /     &
!                      MAX(1.E-18,AC2(ID,IS,IGC)*SIGPI)
                       max(1.e-18,ac2(id,is,ig)*sigpi)
         END IF
       END DO
     END DO

   END IF

   DEALLOCATE(e,sa)

   RETURN

swsnl1()

subroutine swsnl1	(	integer, dimension(*)	WWINT,
		real, dimension(*)	WWAWG,
		real, dimension(*)	WWSWG,
		integer, dimension(msc)	IDCMIN,
		integer, dimension(msc)	IDCMAX,
		real, dimension(msc4mi:msc4ma , mdc4mi:mdc4ma )	UE,
		real, dimension(msc4mi:msc4ma , mdc4mi:mdc4ma )	SA1,
		real, dimension(msc4mi:msc4ma , mdc4mi:mdc4ma )	SA2,
		real, dimension(msc4mi:msc4ma , mdc4mi:mdc4ma )	DA1C,
		real, dimension(msc4mi:msc4ma , mdc4mi:mdc4ma )	DA1P,
		real, dimension(msc4mi:msc4ma , mdc4mi:mdc4ma )	DA1M,
		real, dimension(msc4mi:msc4ma , mdc4mi:mdc4ma )	DA2C,
		real, dimension(msc4mi:msc4ma , mdc4mi:mdc4ma )	DA2P,
		real, dimension(msc4mi:msc4ma , mdc4mi:mdc4ma )	DA2M,
		real, dimension(msc)	SPCSIG,
		real	SNLC1,
		real	KMESPC,
		real	FACHFR,
			ISSTOP,
		real	DAL1,
		real	DAL2,
		real	DAL3,
		real, dimension(msc4mi:msc4ma , mdc4mi:mdc4ma )	SFNL,
		real, dimension(msc4mi:msc4ma , mdc4mi:mdc4ma )	DSNL,
		real, dimension(mt)	DEP2,
		real, dimension(mdc,msc,0:mt)	AC2,
		real, dimension(mdc,msc)	IMATDA,
		real, dimension(mdc,msc)	IMATRA,
		real, dimension(mdc,msc,nptst)	PLNL4S,
		real, dimension(mdc,msc,nptst)	PLNL4D,
		integer	IDDLOW,
		integer	IDDTOP
	)

Definition at line 2254 of file swancom4.f90.

!(This subroutine has not been tested yet and not useful any more)
!
!********************************************************************
!
   USE swcomm3                                                         
   USE swcomm4                                                         
   USE ocpcomm4                                                        
   USE m_snl4   
   USE all_vars, ONLY : mt                                                       
!
!   --|-----------------------------------------------------------|--
!     | Delft University of Technology                            |
!     | Faculty of Civil Engineering                              |
!     | Fluid Mechanics Section                                   |
!     | P.O. Box 5048, 2600 GA  Delft, The Netherlands            |
!     |                                                           |
!     | Programmers: H.L. Tolman, R.C. Ris                        |
!   --|-----------------------------------------------------------|--
!
!
!     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
!     40.13: Nico Booij
!     40.17: IJsbrand Haagsma
!     40.23: Marcel Zijlema
!     40.41: Marcel Zijlema
!
!  1. Updates
!
!     30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN
!     40.17, Dec. 01: Implentation of Multiple DIA
!     40.23, Aug. 02: some corrections
!     40.41, Oct. 04: common blocks replaced by modules, include files removed
!
!  2. Purpose
!
!     Calculate non-linear interaction using the discrete interaction
!     approximation (Hasselmann and Hasselmann 1985; WAMDI group 1988),
!     including the diagonal term for the implicit integration.
!
!     The interactions are calculated for all bin's that fall
!     within a sweep. No additional auxiliary array is required (see
!     SWSNL3)
!
!  3. Method
!
!     Discrete interaction approximation.
!
!     Since the domain in directional domain is by definition not
!     periodic, the spectral space can not beforehand
!     folded to the side angles. This can only be done if the
!     full circle has to be calculated
!
!
!                            Frequencies -->
!                 +---+---------------------+---------+- IDHGH
!              d  | 3 :          2          :    2    |
!              i  + - + - - - - - - - - - - + - - - - +- MDC
!              r  |   :                     :         |
!              e  | 3 :  original spectrum  :    1    |
!              c  |   :                     :         |
!              t. + - + - - - - - - - - - - + - - - - +- 1
!                 | 3 :          2          :    2    |
!                 +---+---------------------+---------+- IDLOW
!                 |   |                     |    ^    |
!             ISLOW   1                     MSC  |    ISHGH
!                     ^                          |
!                     |                          |
!                    ISCLW                     ISCHG
!              lowest discrete               highest discrete
!                central bin                   central bin
!
!                            1 : Extra tail added beyond MSC
!                            2 : Spectrum copied outside ID range
!                            3 : Empty bins at low frequencies
!
!     ISLOW =  1  + ISM1
!     ISHGH = MSC + ISP1 - ISM1
!     ISCLW =  1
!     ISCHG = MSC - ISM1
!     IDLOW =  IDDLOW - MAX(IDM1,IDP1)
!     IDHGH =  IDDTOP + MAX(IDM1,IDP1)
!
!     For the meaning of the counters on the right hand side of the
!     above equations see section 4.
!
!  4. Argument variables
!
!     SPCSIG: Relative frequencies in computational domain in sigma-space 
!
   REAL    SPCSIG(MSC)                                                 
!
!     Data in PARAMETER statements :
!     ----------------------------------------------------------------
!       DAL1    Real  LAMBDA dependend weight factors (see FAC4WW)
!       DAL2    Real
!       DAL3    Real
!       ITHP, ITHP1, ITHM, ITHM1, IFRP, IFRP1, IFRM, IFRM1
!               Int.  Counters of interpolation point relative to
!                     central bin, see figure below (set in FAC4WW).
!       NFRLOW, NFRHGH, NFRCHG, NTHLOW, NTHHGH
!               Int.  Range of calculations, see section 2.
!       AF11    R.A.  Scaling array (Freq**11).
!       AWGn    Real  Interpolation weights, see numbers in fig.
!       SWGn    Real  Id. squared.
!       UE      R.A.  "Unfolded" spectrum.
!       SA1     R.A.  Interaction constribution of first and second
!       SA2     R.A.    quadr. respectively (unfolded space).
!       DA1C, DA1P, DA1M, DA2C, DA2P, DA2M
!               R.A.  Idem for diagonal matrix.
!       PERCIR        full circle or sector
!     ----------------------------------------------------------------
!
!       Relative offsets of interpolation points around central bin
!       "#" and corresponding numbers of AWGn :
!
!               ISM1  ISM
!                5        7    T |
!          IDM1   +------+     H +
!                 |      |     E |      ISP      ISP1
!                 |   \  |     T |       3           1
!           IDM   +------+     A +        +---------+  IDP1
!                6       \8      |        |         |
!                                |        |  /      |
!                           \    +        +---------+  IDP
!                                |      /4           2
!                              \ |  /
!          -+-----+------+-------#--------+---------+----------+
!                                |           FREQ.
!
!  7. Common blocks used
!
!
!  8. Subroutines used
!
!     ---
!
!  9. Subroutines calling
!
!     SOURCE (in SWANCOM1)
!
! 12. Structure
!
!     -------------------------------------------
!       Initialisations.
!       Calculate proportionality constant.
!       Prepare auxiliary spectrum.
!       Calculate interactions :
!       -----------------------------------------
!         Energy at interacting bins
!         Contribution to interactions
!         Fold interactions to side angles
!       -----------------------------------------
!       Put source term together
!     -------------------------------------------
!
! 13. Source text
!
!*************************************************************
!
   INTEGER   IS     ,ID     ,I      ,J      ,                       &
             ISHGH  ,IDLOW  ,ISP    ,ISP1   ,IDP    ,IDP1   ,       &
         ISM    ,ISM1   ,IDHGH  ,IDM    ,IDM1   ,ISCLW  ,       &
         ISCHG  ,IDDLOW ,IDDTOP                                    
!
   REAL      X      ,X2     ,CONS   ,FACTOR ,SNLCS1 ,SNLCS2 ,SNLCS3, &
             E00    ,EP1    ,EM1    ,EP2    ,EM2    ,SA1A   ,SA1B  , &
         SA2A   ,SA2B   ,KMESPC ,FACHFR ,AWG1   ,AWG2   ,AWG3  , &
         AWG4   ,AWG5   ,AWG6   ,AWG7   ,AWG8   ,DAL1   ,DAL2  , &
         DAL3   ,SNLC1  ,SWG1   ,SWG2   ,SWG3   ,SWG4   ,SWG5  , &
         SWG6   ,SWG7   ,SWG8           ,JACOBI ,SIGPI            
!
   REAL      AC2(MDC,MSC,0:MT)                    ,                 &
             DEP2(MT)                           ,                 &
         UE(MSC4MI:MSC4MA , MDC4MI:MDC4MA )    ,                 &
         SA1(MSC4MI:MSC4MA , MDC4MI:MDC4MA )   ,                 &
         SA2(MSC4MI:MSC4MA , MDC4MI:MDC4MA )   ,                 &
         DA1C(MSC4MI:MSC4MA , MDC4MI:MDC4MA )  ,                 &
         DA1P(MSC4MI:MSC4MA , MDC4MI:MDC4MA )  ,                 &
         DA1M(MSC4MI:MSC4MA , MDC4MI:MDC4MA )  ,                 &
         DA2C(MSC4MI:MSC4MA , MDC4MI:MDC4MA )  ,                 &
         DA2P(MSC4MI:MSC4MA , MDC4MI:MDC4MA )  ,                 &
         DA2M(MSC4MI:MSC4MA , MDC4MI:MDC4MA )  ,                 &
         SFNL(MSC4MI:MSC4MA , MDC4MI:MDC4MA )  ,                 &
         DSNL(MSC4MI:MSC4MA , MDC4MI:MDC4MA )  ,                 &
         IMATDA(MDC,MSC)                       ,                 &
         IMATRA(MDC,MSC)                       ,                 &
         PLNL4S(MDC,MSC,NPTST)                 ,                 &
         PLNL4D(MDC,MSC,NPTST)                 ,                 &
         WWAWG(*)                              ,                 &
         WWSWG(*)
!
   INTEGER   IDCMIN(MSC),IDCMAX(MSC),WWINT(*)
!
   LOGICAL   PERCIR
!
!   SAVE IENT
!   DATA IENT/0/
!   IF (LTRACE) CALL STRACE (IENT,'SWSNL1')
!
   idp    = wwint(1)
   idp1   = wwint(2)
   idm    = wwint(3)
   idm1   = wwint(4)
   isp    = wwint(5)
   isp1   = wwint(6)
   ism    = wwint(7)
   ism1   = wwint(8)
   islow  = wwint(9)
   ishgh  = wwint(10)
   isclw  = wwint(11)
   ischg  = wwint(12)
   idlow  = wwint(13)
   idhgh  = wwint(14)
!
   awg1 = wwawg(1)
   awg2 = wwawg(2)
   awg3 = wwawg(3)
   awg4 = wwawg(4)
   awg5 = wwawg(5)
   awg6 = wwawg(6)
   awg7 = wwawg(7)
   awg8 = wwawg(8)
!
   swg1 = wwswg(1)
   swg2 = wwswg(2)
   swg3 = wwswg(3)
   swg4 = wwswg(4)
   swg5 = wwswg(5)
   swg6 = wwswg(6)
   swg7 = wwswg(7)
   swg8 = wwswg(8)
!
!  *** Initialize auxiliary arrays per gridpoint ***
!
   DO id = mdc4mi, mdc4ma
     DO is = msc4mi, msc4ma
       ue(is,id)   = 0.
       sa1(is,id)  = 0.
       sa2(is,id)  = 0.
       sfnl(is,id) = 0.
       da1c(is,id) = 0.
       da1p(is,id) = 0.
       da1m(is,id) = 0.
       da2c(is,id) = 0.
       da2p(is,id) = 0.
       da2m(is,id) = 0.
       dsnl(is,id) = 0.
     ENDDO
   ENDDO
!
!  *** Calculate factor R(X) to calculate the NL wave-wave ***
!  *** interaction for shallow water                       ***
!  *** SNLC1 = 1/GRAV**4                                   ***         
!
   snlcs1 = pquad(3)                                                   
   snlcs2 = pquad(4)                                                   
   snlcs3 = pquad(5)                                                   
!   X      = MAX ( 0.75 * DEP2(IGC) * KMESPC , 0.5 )
   x      = max( 0.75 * dep2(kcgrd(1)) * kmespc , 0.5 )
   x2     = max( -1.e15, snlcs3*x)
   cons   = snlc1 * ( 1. + snlcs1/x * (1.-snlcs2*x) * exp(x2))
   jacobi = 2. * pi_w
!
!  *** check whether the spectral domain is periodic in ***
!  *** directional space and if so, modify boundaries   ***
!
   percir = .false.
   IF(iddlow == 1 .AND. iddtop == mdc)THEN
!    *** periodic in theta -> spectrum can be folded    ***
!    *** (can only be present in presence of a current) ***
     idclow = 1
     idchgh = mdc
     iiid   = 0
     percir = .true.
   ELSE
!    *** different sectors per sweep -> extend range with IIID ***
     iiid   = max( idm1 , idp1 )
     idclow = idlow
     idchgh = idhgh
   ENDIF
!
!  *** Prepare auxiliary spectrum               ***
!  *** set action original spectrum in array UE ***
!
   DO iddum = idlow - iiid, idhgh + iiid
     id = mod( iddum - 1 + mdc , mdc ) + 1
     DO is = 1, msc
!       UE(IS,IDDUM) = AC2(ID,IS,IGC) * SPCSIG(IS) * JACOBI        
       ue(is,iddum) = ac2(id,is,kcgrd(1)) * spcsig(is) * jacobi        
     ENDDO
   ENDDO
!
!  *** set values in area 2 for IS > MSC+1  ***
!
   DO is = msc+1, ishgh
     DO id = idlow - iiid , idhgh + iiid
       ue(is,id) = ue(is-1,id) * fachfr
     ENDDO
   ENDDO
!
!  *** Calculate interactions      ***
!  *** Energy at interacting bins  ***
!
   DO is = isclw, ischg
     DO id = idclow, idchgh
       e00    =        ue(is      ,id      )
       ep1    = awg1 * ue(is+isp1,id+idp1) +                  &
                awg2 * ue(is+isp1,id+idp ) +                  &
        awg3 * ue(is+isp ,id+idp1) +                  &
        awg4 * ue(is+isp ,id+idp )
       em1    = awg5 * ue(is+ism1,id-idm1) +                  &
                awg6 * ue(is+ism1,id-idm ) +                  &
        awg7 * ue(is+ism ,id-idm1) +                  &
        awg8 * ue(is+ism ,id-idm )
!
       ep2    = awg1 * ue(is+isp1,id-idp1) +                  &
                awg2 * ue(is+isp1,id-idp ) +                  &
        awg3 * ue(is+isp ,id-idp1) +                  &
        awg4 * ue(is+isp ,id-idp )
       em2    = awg5 * ue(is+ism1,id+idm1) +                  &
                awg6 * ue(is+ism1,id+idm ) +                  &
        awg7 * ue(is+ism ,id+idm1) +                  &
        awg8 * ue(is+ism ,id+idm )
!
!      *** Contribution to interactions                          ***
!      *** CONS is the shallow water factor for the NL interact. ***
!
       factor = cons * af11(is) * e00
!
       sa1a   = e00 * ( ep1*dal1 + em1*dal2 ) * pquad(2)               
       sa1b   = sa1a - ep1*em1*dal3 * pquad(2)                         
       sa2a   = e00 * ( ep2*dal1 + em2*dal2 ) * pquad(2)               
       sa2b   = sa2a - ep2*em2*dal3 * pquad(2)                         
!
       sa1(is,id) = factor * sa1b
       sa2(is,id) = factor * sa2b
!
!       IF(ITEST >= 100 .AND. TESTFL)THEN
!         WRITE(PRINTF,9002) E00,EP1,EM1,EP2,EM2
!9002     FORMAT (' E00 EP1 EM1 EP2 EM2  :',5E11.4)
!         WRITE(PRINTF,9003) SA1A,SA1B,SA2A,SA2B
!9003     FORMAT (' SA1A SA1B SA2A SA2B  :',4E11.4)
!         WRITE(PRINTF,9004) IS,ID,SA1(IS,ID),SA2(IS,ID)
!9004     FORMAT (' IS ID SA1() SA2()    :',2I4,2E12.4)
!         WRITE(PRINTF,9005) FACTOR
!9005     FORMAT (' FACTOR               : ',E12.4)
!       END IF
!
       da1c(is,id) = cons * af11(is) * ( sa1a + sa1b )
       da1p(is,id) = factor * ( dal1*e00 - dal3*em1 ) * pquad(2)       
       da1m(is,id) = factor * ( dal2*e00 - dal3*ep1 ) * pquad(2)       
!
       da2c(is,id) = cons * af11(is) * ( sa2a + sa2b )
       da2p(is,id) = factor * ( dal1*e00 - dal3*em2 ) * pquad(2)       
       da2m(is,id) = factor * ( dal2*e00 - dal3*ep2 ) * pquad(2)       
     ENDDO
   ENDDO
!
!  *** Fold interactions to side angles if spectral domain ***
!  *** is periodic in directional space                    ***
!
   IF(percir)THEN
     DO id = 1, idhgh - mdc
       id0   = 1 - id
       DO is = isclw, ischg
         sa1(is,mdc+id) = sa1(is,  id   )
         sa2(is,mdc+id) = sa2(is,  id   )
         da1c(is,mdc+id) = da1c(is,  id   )
         da1p(is,mdc+id) = da1p(is,  id   )
         da1m(is,mdc+id) = da1m(is,  id   )
         da2c(is,mdc+id) = da2c(is,  id   )
         da2p(is,mdc+id) = da2p(is,  id   )
         da2m(is,mdc+id) = da2m(is,  id   )
!
         sa1(is,  id0 ) = sa1(is, mdc+id0)
         sa2(is,  id0 ) = sa2(is, mdc+id0)
         da1c(is,  id0 ) = da1c(is, mdc+id0)
         da1p(is,  id0 ) = da1p(is, mdc+id0)
         da1m(is,  id0 ) = da1m(is, mdc+id0)
         da2c(is,  id0 ) = da2c(is, mdc+id0)
         da2p(is,  id0 ) = da2p(is, mdc+id0)
         da2m(is,  id0 ) = da2m(is, mdc+id0)
       ENDDO
     ENDDO
   ENDIF
!
!  *** Put source term together (To save space I=IS and J=ID ***
!  *** is used)                                              ***
!
   pi3   = (2. * pi_w)**3
   DO i = 1, isstop
     sigpi = spcsig(i) * jacobi                                        
     DO j = idcmin(i), idcmax(i)
       id = mod( j - 1 + mdc , mdc ) + 1
       sfnl(i,id) = - 2. * ( sa1(i,j) + sa2(i,j) )                       &
                    + awg1 * ( sa1(i-isp1,j-idp1) + sa2(i-isp1,j+idp1) ) &
            + awg2 * ( sa1(i-isp1,j-idp ) + sa2(i-isp1,j+idp ) ) &
            + awg3 * ( sa1(i-isp ,j-idp1) + sa2(i-isp ,j+idp1) ) &
            + awg4 * ( sa1(i-isp ,j-idp ) + sa2(i-isp ,j+idp ) ) &
            + awg5 * ( sa1(i-ism1,j+idm1) + sa2(i-ism1,j-idm1) ) &
            + awg6 * ( sa1(i-ism1,j+idm ) + sa2(i-ism1,j-idm ) ) &
            + awg7 * ( sa1(i-ism ,j+idm1) + sa2(i-ism ,j-idm1) ) &
            + awg8 * ( sa1(i-ism ,j+idm ) + sa2(i-ism ,j-idm ) )
!
       dsnl(i,id) = - 2. * ( da1c(i,j) + da2c(i,j) )                       &
                    + swg1 * ( da1p(i-isp1,j-idp1) + da2p(i-isp1,j+idp1) ) &
            + swg2 * ( da1p(i-isp1,j-idp ) + da2p(i-isp1,j+idp ) ) &
            + swg3 * ( da1p(i-isp ,j-idp1) + da2p(i-isp ,j+idp1) ) &
            + swg4 * ( da1p(i-isp ,j-idp ) + da2p(i-isp ,j+idp ) ) &
            + swg5 * ( da1m(i-ism1,j+idm1) + da2m(i-ism1,j-idm1) ) &
            + swg6 * ( da1m(i-ism1,j+idm ) + da2m(i-ism1,j-idm ) ) &
            + swg7 * ( da1m(i-ism ,j+idm1) + da2m(i-ism ,j-idm1) ) &
            + swg8 * ( da1m(i-ism ,j+idm ) + da2m(i-ism ,j-idm ) )
!
!      *** store results in IMATDA and IMATRA ***
!
       IF(testfl) THEN
         plnl4s(id,i,iptst) = sfnl(i,id) / sigpi                       
         plnl4d(id,i,iptst) = -1. * dsnl(i,id) / pi3                   
       END IF
!
       imatra(id,i) = imatra(id,i) + sfnl(i,id) / sigpi
       imatda(id,i) = imatda(id,i) - dsnl(i,id) / pi3
!
!       IF(ITEST >= 90 .AND. TESTFL) THEN
!         WRITE(PRINTF,9006) I,J,SFNL(I,ID),DSNL(I,ID),SPCSIG(I)                                                    30.72
!9006     FORMAT (' IS ID SFNL DSNL SPCSIG:',2I4,3E12.4)                
!       END IF
!
     ENDDO
   ENDDO
!

   RETURN

swsnl2()

subroutine swsnl2	(	integer	IDDLOW,
		integer	IDDTOP,
		integer, dimension(*)	WWINT,
		real, dimension(*)	WWAWG,
		real, dimension(msc4mi:msc4ma , mdc4mi:mdc4ma )	UE,
		real, dimension(msc4mi:msc4ma , mdc4mi:mdc4ma )	SA1,
		integer	ISSTOP,
		real, dimension(msc4mi:msc4ma , mdc4mi:mdc4ma )	SA2,
		real, dimension(msc)	SPCSIG,
		real	SNLC1,
		real	DAL1,
		real	DAL2,
		real	DAL3,
		real, dimension(msc4mi:msc4ma , mdc4mi:mdc4ma)	SFNL,
		real, dimension(mt)	DEP2,
		real	KMESPC,
		real, dimension(mdc,msc)	IMATDA,
		real, dimension(mdc,msc)	IMATRA,
		real	FACHFR,
		real, dimension(mdc,msc,nptst)	PLNL4S,
		integer, dimension(msc)	IDCMIN,
		integer, dimension(msc)	IDCMAX,
		integer	IG,
		real, dimension(msc,micmax)	CGO,
		real, dimension(*)	WWSWG
	)

Definition at line 1959 of file swancom4.f90.

!
!*******************************************************************
!
!     Calculate non-linear interaction using the discrete interaction
!     approximation (Hasselmann and Hasselmann 1985; WAMDI group 1988)
!
!*******************************************************************
   USE swcomm3                                                         
   USE swcomm4                                                         
   USE ocpcomm4                                                        
   USE m_snl4  
!   USE ALL_VARS, ONLY : MT,AC2                   
   USE vars_wave, ONLY : mt,ac2                   
   
   IMPLICIT NONE                                     
!
   REAL    :: SPCSIG(MSC)                                                 
   INTEGER :: IS     ,ID     ,I      ,J      ,IG     ,ISHGH  ,  &
              ISSTOP ,ISP    ,ISP1   ,IDP    ,IDP1   ,ISM    ,ISM1   ,  &
          IDM    ,IDM1   ,ISCLW  ,ISCHG  ,                          &
          IDLOW  ,IDHGH  ,IDDLOW ,IDDTOP ,IDCLOW ,IDCHGH    
!
   REAL    :: X      ,X2     ,CONS   ,FACTOR ,SNLCS1 ,SNLCS2 ,SNLCS3 ,  &
              E00    ,EP1    ,EM1    ,EP2    ,EM2    ,SA1A   ,SA1B   ,  &
          SA2A   ,SA2B   ,KMESPC ,FACHFR ,AWG1   ,AWG2   ,AWG3   ,  &
          AWG4   ,AWG5   ,AWG6   ,AWG7   ,AWG8   ,DAL1   ,DAL2   ,  &
          DAL3           ,JACOBI ,SIGPI                            
!
    REAL   :: DEP2(MT)                              ,                   &
          UE(MSC4MI:MSC4MA , MDC4MI:MDC4MA )    ,                   &
          SA1(MSC4MI:MSC4MA , MDC4MI:MDC4MA )   ,                   &
          SA2(MSC4MI:MSC4MA , MDC4MI:MDC4MA )   ,                   &
          DA1C(MSC4MI:MSC4MA , MDC4MI:MDC4MA )  ,                   &
          DA1P(MSC4MI:MSC4MA , MDC4MI:MDC4MA )  ,                   &
          DA1M(MSC4MI:MSC4MA , MDC4MI:MDC4MA )  ,                   &
          DA2C(MSC4MI:MSC4MA , MDC4MI:MDC4MA )  ,                   &
          DA2P(MSC4MI:MSC4MA , MDC4MI:MDC4MA )  ,                   &
          DA2M(MSC4MI:MSC4MA , MDC4MI:MDC4MA )  ,                   &
          SFNL(MSC4MI:MSC4MA , MDC4MI:MDC4MA)   ,                   &
          DFNL(MSC4MI:MSC4MA , MDC4MI:MDC4MA)   ,                   &
          IMATRA(MDC,MSC)                       ,                   &
          IMATDA(MDC,MSC)                       ,                   &
          PLNL4S(MDC,MSC,NPTST)                 ,                   &
          WWAWG(*)                              ,                   &
          CGO(MSC,MICMAX)
!
   INTEGER :: WWINT(*)         ,IDCMIN(MSC)      ,IDCMAX(MSC)
   INTEGER :: ISLOW,IIID,IDDUM,ID0
   REAL    :: SNLC1,SWG1,SWG2,SWG3,SWG4,SWG5,SWG6,SWG7,SWG8
   REAL    :: WWSWG(*),PI3
   
!
   LOGICAL   PERCIR
!
   idp    = wwint(1)
   idp1   = wwint(2)
   idm    = wwint(3)
   idm1   = wwint(4)
   isp    = wwint(5)
   isp1   = wwint(6)
   ism    = wwint(7)
   ism1   = wwint(8)
   islow  = wwint(9)
   ishgh  = wwint(10)
   isclw  = wwint(11)
   ischg  = wwint(12)
   idlow  = wwint(13)
   idhgh  = wwint(14)

   awg1 = wwawg(1)
   awg2 = wwawg(2)
   awg3 = wwawg(3)
   awg4 = wwawg(4)
   awg5 = wwawg(5)
   awg6 = wwawg(6)
   awg7 = wwawg(7)
   awg8 = wwawg(8)
!
   swg1 = wwswg(1)
   swg2 = wwswg(2)
   swg3 = wwswg(3)
   swg4 = wwswg(4)
   swg5 = wwswg(5)
   swg6 = wwswg(6)
   swg7 = wwswg(7)
   swg8 = wwswg(8)
!
!  *** Initialize auxiliary arrays per gridpoint ***
!
   DO id = mdc4mi, mdc4ma
     DO is = msc4mi, msc4ma
       ue(is,id)   = 0.
       sa1(is,id)  = 0.
       sa2(is,id)  = 0.
       
       da1c(is,id) = 0.
       da1p(is,id) = 0.
       da1m(is,id) = 0.
       da2c(is,id) = 0.
       da2p(is,id) = 0.
       da2m(is,id) = 0.
       
       sfnl(is,id) = 0.
       dfnl(is,id) = 0.
       
     ENDDO
   ENDDO
!
!  *** Calculate prop. constant.                           ***
!  *** Calculate factor R(X) to calculate the NL wave-wave ***
!  *** interaction for shallow water                       ***
!  *** SNLC1 = 1/GRAV**4                                   ***         
!
   snlcs1 = pquad(3)                                                   
   snlcs2 = pquad(4)                                                   
   snlcs3 = pquad(5)    

   x      = max( 0.75 * dep2(ig) * kmespc , 0.5 )
   x2     = max( -1.e15, snlcs3*x)
   cons   = snlc1 * ( 1. + snlcs1/x * (1.-snlcs2*x) * exp(x2))
   jacobi = 2. * pi_w
!
!  *** check whether the spectral domain is periodic in ***
!  *** direction space and if so modify boundaries      ***
!
   percir = .false.
   IF(iddlow == 1 .AND. iddtop == mdc)THEN
!    *** periodic in theta -> spectrum can be folded  ***
!    *** (can only occur in presence of a current)    ***
     idclow = 1
     idchgh = mdc
     iiid   = 0
     percir = .true.
   ELSE
!    *** different sectors per sweep -> extend range with IIID ***
     iiid   = max( idm1 , idp1 )
     idclow = idlow
     idchgh = idhgh
   ENDIF
!
!  *** Prepare auxiliary spectrum               ***
!  *** set action original spectrum in array UE ***
!
!print*,'swsnl2-1'
!print*,'IDLOW,IDHGH,IIID',IDLOW,IDHGH,IIID
!print*,'low high MDC=',IDLOW - IIID, IDHGH + IIID,MDC

   DO iddum = idlow - iiid , idhgh + iiid
     id = mod( iddum - 1 + mdc , mdc ) + 1
     DO is = 1, msc
!      print*,'ID,IS,IG=',ID,IS,IG
       ue(is,iddum) = ac2(id,is,ig) * spcsig(is) * jacobi      
     ENDDO
   ENDDO
!print*,'swsnl2-2'
!
!  *** set values in the areas 2 for IS > MSC+1 ***
!
   DO is = msc+1, ishgh
     DO id = idlow - iiid , idhgh + iiid
       ue(is,id) = ue(is-1,id) * fachfr
     ENDDO
   ENDDO
!
!  *** Calculate interactions      ***
!  *** Energy at interacting bins  ***
!
   DO is = isclw, ischg
     DO id = idclow , idchgh
       e00 =        ue(is     ,id      )
       ep1 = awg1 * ue(is+isp1,id+idp1) +                   &
             awg2 * ue(is+isp1,id+idp ) +                   &
         awg3 * ue(is+isp ,id+idp1) +                   &
         awg4 * ue(is+isp ,id+idp )
       em1 = awg5 * ue(is+ism1,id-idm1) +                   &
             awg6 * ue(is+ism1,id-idm ) +                   &
         awg7 * ue(is+ism ,id-idm1) +                   &
         awg8 * ue(is+ism ,id-idm )
!
       ep2 = awg1 * ue(is+isp1,id-idp1) +                   &
             awg2 * ue(is+isp1,id-idp ) +                   &
         awg3 * ue(is+isp ,id-idp1) +                   &
         awg4 * ue(is+isp ,id-idp )
       em2 = awg5 * ue(is+ism1,id+idm1) +                   &
             awg6 * ue(is+ism1,id+idm ) +                   &
         awg7 * ue(is+ism ,id+idm1) +                   &
         awg8 * ue(is+ism ,id+idm )
!
!      *** Contribution to interactions                          ***
!      *** CONS is the shallow water factor for the NL interact. ***
!
       factor = cons * af11(is) * e00
!
       sa1a = e00 * ( ep1*dal1 + em1*dal2 ) * pquad(2)               
       sa1b = sa1a - ep1*em1*dal3 * pquad(2)                         
       sa2a = e00 * ( ep2*dal1 + em2*dal2 ) * pquad(2)               
       sa2b = sa2a - ep2*em2*dal3 * pquad(2)                         
!
       sa1(is,id) = factor * sa1b
       sa2(is,id) = factor * sa2b
       
       da1c(is,id) = cons*af11(is)*(sa1a+sa1b)
       da1p(is,id) = factor*(dal1*e00-dal3*em1)*pquad(2)
       da1m(is,id) = factor*(dal2*e00-dal3*ep1)*pquad(2)
       
       da2c(is,id) = cons*af11(is)*(sa2a+sa2b)
       da2p(is,id) = factor*(dal1*e00-dal3*em2)*pquad(2)
       da2m(is,id) = factor*(dal2*e00-dal3*ep2)*pquad(2)
       
     ENDDO
   ENDDO
!
!  *** Fold interactions to side angles if spectral domain ***
!  *** is periodic in directional space                    ***
!
   IF(percir)THEN
     DO id = 1, idhgh - mdc
       id0   = 1 - id
       DO is = isclw, ischg
         sa1(is,mdc+id) = sa1(is ,  id    )
         sa2(is,mdc+id) = sa2(is ,  id    )
         sa1(is,  id0 ) = sa1(is , mdc+id0)
         sa2(is,  id0 ) = sa2(is , mdc+id0)
       
       da1c(is,mdc+id) = da1c(is,id)
       da1p(is,mdc+id) = da1p(is,id)
       da1m(is,mdc+id) = da1m(is,id)
       da1c(is,id0) = da1c(is,mdc+id0)
       da1p(is,id0) = da1p(is,mdc+id0)
       da1m(is,id0) = da1m(is,mdc+id0)
       
       da2c(is,mdc+id) = da2c(is,id)
       da2p(is,mdc+id) = da2p(is,id)
       da2m(is,mdc+id) = da2m(is,id)
       da2c(is,id0) = da2c(is,mdc+id0)
       da2p(is,id0) = da2p(is,mdc+id0)
       da2m(is,id0) = da2m(is,mdc+id0)
       ENDDO
     ENDDO
   ENDIF
!
!  ***  Put source term together (To save space I=IS and J=ID ***
!  ***  is used)                                              ***
!
   pi3 = (2.0*pi_w)**3
   DO i = 1, isstop
     sigpi = spcsig(i) * jacobi 
     
     DO j = idcmin(i), idcmax(i)
       id = mod( j - 1 + mdc , mdc ) + 1
       sfnl(i,id) = - 2. * ( sa1(i,j) + sa2(i,j) )                       &
                    + awg1 * ( sa1(i-isp1,j-idp1) + sa2(i-isp1,j+idp1) ) &
            + awg2 * ( sa1(i-isp1,j-idp ) + sa2(i-isp1,j+idp ) ) &
            + awg3 * ( sa1(i-isp ,j-idp1) + sa2(i-isp ,j+idp1) ) &
            + awg4 * ( sa1(i-isp ,j-idp ) + sa2(i-isp ,j+idp ) ) &
            + awg5 * ( sa1(i-ism1,j+idm1) + sa2(i-ism1,j-idm1) ) &
            + awg6 * ( sa1(i-ism1,j+idm ) + sa2(i-ism1,j-idm ) ) &
            + awg7 * ( sa1(i-ism ,j+idm1) + sa2(i-ism ,j-idm1) ) &
            + awg8 * ( sa1(i-ism ,j+idm ) + sa2(i-ism ,j-idm ) )

       dfnl(i,id) =   - 2. * ( da1c(i,j) + da2c(i,j) )                      &
                    + swg1 * ( da1p(i-isp1,j-idp1) + da2p(i-isp1,j+idp1) )  &
            + swg2 * ( da1p(i-isp1,j-idp ) + da2p(i-isp1,j+idp ) )  &
            + swg3 * ( da1p(i-isp ,j-idp1) + da2p(i-isp ,j+idp1) )  &
            + swg4 * ( da1p(i-isp ,j-idp ) + da2p(i-isp ,j+idp ) )  &
            + swg5 * ( da1m(i-ism1,j+idm1) + da2m(i-ism1,j-idm1) )  &
            + swg6 * ( da1m(i-ism1,j+idm ) + da2m(i-ism1,j-idm ) )  &
            + swg7 * ( da1m(i-ism ,j+idm1) + da2m(i-ism ,j-idm1) )  &
            + swg8 * ( da1m(i-ism ,j+idm ) + da2m(i-ism ,j-idm ) )
            
!
!      *** store results in rhv ***
!      *** store results in rhs and main diagonal according ***        
!      *** to Patankar-rules                                ***        
!
       IF(testfl) plnl4s(id,i,iptst) =  sfnl(i,id) / sigpi  
       imatra(id,i) = imatra(id,i) + sfnl(i,id) / sigpi
       imatda(id,i) = imatda(id,i) + dfnl(i,id)     
     
     ENDDO
   ENDDO

   RETURN

swsnl3()

subroutine swsnl3	(	integer, dimension(*)	WWINT,
		real, dimension(*)	WWAWG,
		real, dimension(msc4mi:msc4ma,mdc4mi:mdc4ma)	UE,
		real, dimension(msc4mi:msc4ma,mdc4mi:mdc4ma)	SA1,
		real, dimension(msc4mi:msc4ma,mdc4mi:mdc4ma)	SA2,
		real, dimension(msc)	SPCSIG,
		real	SNLC1,
		real	DAL1,
		real	DAL2,
		real	DAL3,
		real, dimension(msc4mi:msc4ma,mdc4mi:mdc4ma)	SFNL,
		real, dimension(mt)	DEP2,
		real	KMESPC,
		real, dimension(mdc,msc,mt)	MEMNL4,
		real	FACHFR
	)

Definition at line 1558 of file swancom4.f90.

!  (This subroutine has not been tested yet)              
!
!*******************************************************************
!
   USE swcomm3                                                         
   USE swcomm4                                                         
   USE ocpcomm4                                                        
   USE m_snl4 
!   USE ALL_VARS, ONLY : MT,AC2                
   USE vars_wave, ONLY : mt,ac2                
!
!   --|-----------------------------------------------------------|--
!     | Delft University of Technology                            |
!     | Faculty of Civil Engineering                              |
!     | Fluid Mechanics Section                                   |
!     | P.O. Box 5048, 2600 GA  Delft, The Netherlands            |
!     |                                                           |
!     | Programmers: H.L. Tolman, R.C. Ris                        |
!   --|-----------------------------------------------------------|--
!
!
!     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
!     40.17: IJsbrand Haagsma
!     40.41: Marcel Zijlema
!
!  1. Updates
!
!     30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN
!     40.17, Dec. 01: Implemented Multiple DIA
!     40.41, Oct. 04: common blocks replaced by modules, include files removed
!
!  2. Purpose
!
!     Calculate non-linear interaction using the discrete interaction
!     approximation (Hasselmann and Hasselmann 1985; WAMDI group 1988)
!     for the full circle (option if a current is present). Note: using
!     this subroutine requires an additional array with size
!     (MXC*MYC*MDC*MSC). This requires more internal memory but can
!     speed up the computations sigificantly if a current is present.
!
!  3. Method
!
!     Discrete interaction approximation. To make interpolation simple,
!     the interactions are calculated in a "folded" space.
!
!                            Frequencies -->
!                 +---+---------------------+---------+- IDHGH
!              d  | 3 :          2          :    2    |
!              i  + - + - - - - - - - - - - + - - - - +- MDC
!              r  |   :                     :         |
!              e  | 3 :  original spectrum  :    1    |
!              c  |   :                     :         |
!              t. + - + - - - - - - - - - - + - - - - +- 1
!                 | 3 :          2          :    2    |
!                 +---+---------------------+---------+- IDLOW
!                 |   |                     |     ^   |
!              ISLOW  1                    MSC    |   ISHGH
!                     |                           |
!                   ISCLW                        ISCHG
!              lowest discrete               highest discrete
!                central bin                   central bin
!
!                            1 : Extra tail added beyond MSC
!                            2 : Spectrum copied outside ID range
!                            3 : Empty bins at low frequencies
!
!     ISLOW =  1  + ISM1
!     ISHGH = MSC + ISP1 - ISM1
!     ISCLW =  1
!     ISCHG = MSC - ISM1
!     IDLOW =  1  - MAX(IDM1,IDP1)
!     IDHGH = MDC + MAX(IDM1,IDP1)
!
!       Relative offsets of interpolation points around central bin
!       "#" and corresponding numbers of AWGn :
!
!               ISM1  ISM
!                5        7    T |
!          IDM1   +------+     H +
!                 |      |     E |      ISP      ISP1
!                 |   \  |     T |       3           1
!           IDM   +------+     A +        +---------+  IDP1
!                6       \8      |        |         |
!                                |        |  /      |
!                           \    +        +---------+  IDP
!                                |      /4           2
!                              \ |  /
!          -+-----+------+-------#--------+---------+----------+
!                                |           FREQ.
!
!
!  4. Argument variables
!
!     MCGRD : number of wet grid points of the computational grid
!     MDC   : grid points in theta-direction of computational grid
!     MDC4MA: highest array counter in directional space (Snl4)
!     MDC4MI: lowest array counter in directional space (Snl4)
!     MSC   : grid points in sigma-direction of computational grid
!     MSC4MA: highest array counter in frequency space (Snl4)
!     MSC4MI: lowest array counter in frequency space (Snl4)
!     WWINT : counters for quadruplet interactions
!
   INTEGER WWINT(*)
!
!     AC2   : action density
!     AF11  : scaling frequency
!     DAL1  : coefficient for the quadruplet interactions
!     DAL2  : coefficient for the quadruplet interactions
!     DAL3  : coefficient for the quadruplet interactions
!     DEP2  : depth
!     FACHFR
!     KMESPC: mean average wavenumber over full spectrum
!     MEMNL4
!     PI    : circular constant
!     SA1   : interaction contribution of first quadruplet (unfolded space)
!     SA2   : interaction contribution of second quadruplet (unfolded space)
!     SFNL
!     SNLC1
!     SPCSIG: relative frequencies in computational domain in sigma-space
!     UE    : "unfolded" spectrum
!     WWAWG : weight coefficients for the quadruplet interactions
!
   REAL    DAL1, DAL2, DAL3, FACHFR, KMESPC, SNLC1                     
!   REAL    AC2(MDC,MSC,0:MT)
   REAL    DEP2(MT)
   REAL    MEMNL4(MDC,MSC,MT)
   REAL    SA1(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
   REAL    SA2(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
   REAL    SFNL(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
   REAL    SPCSIG(MSC)                                                 
   REAL    UE(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
   REAL    WWAWG(*)
!
!  7. Common blocks used
!
!
!  8. Subroutines used
!
!     ---
!
!  9. Subroutines calling
!
!     SOURCE (in SWANCOM1)
!
! 12. Structure
!
!     -------------------------------------------
!       Initialisations.
!       Calculate proportionality constant.
!       Prepare auxiliary spectrum.
!       Calculate (unfolded) interactions :
!       -----------------------------------------
!         Energy at interacting bins
!         Contribution to interactions
!         Fold interactions to side angles
!       -----------------------------------------
!       Put source term together
!     -------------------------------------------
!
! 13. Source text
!
!*******************************************************************
!
   INTEGER   IS      ,ID      ,ID0     ,I       ,J       ,          &
             ISHGH   ,IDLOW   ,IDHGH   ,ISP     ,ISP1    ,          &
         IDP     ,IDP1    ,ISM     ,ISM1    ,IDM     ,IDM1    , &
         ISCLW   ,ISCHG
!
   REAL      X       ,X2      ,CONS    ,FACTOR  ,SNLCS2  ,          &
             SNLCS3  ,E00     ,EP1     ,EM1     ,EP2     ,EM2     , &
         SA1A    ,SA1B    ,SA2A    ,SA2B    ,                   &
         AWG1    ,AWG2    ,AWG3    ,AWG4    ,AWG5    ,AWG6    , &
         AWG7    ,AWG8    ,                                     &
         JACOBI  ,SIGPI
!
!   SAVE IENT
!   DATA IENT/0/
!   IF (LTRACE) CALL STRACE (IENT,'SWSNL3')
!
   idp    = wwint(1)
   idp1   = wwint(2)
   idm    = wwint(3)
   idm1   = wwint(4)
   isp    = wwint(5)
   isp1   = wwint(6)
   ism    = wwint(7)
   ism1   = wwint(8)
   islow  = wwint(9)
   ishgh  = wwint(10)
   isclw  = wwint(11)
   ischg  = wwint(12)
   idlow  = wwint(13)
   idhgh  = wwint(14)
!
   awg1 = wwawg(1)
   awg2 = wwawg(2)
   awg3 = wwawg(3)
   awg4 = wwawg(4)
   awg5 = wwawg(5)
   awg6 = wwawg(6)
   awg7 = wwawg(7)
   awg8 = wwawg(8)
!
!  *** Initialize auxiliary arrays per gridpoint ***
!
   DO id = mdc4mi, mdc4ma
     DO is = msc4mi, msc4ma
       ue(is,id)   = 0.
       sa1(is,id)  = 0.
       sa2(is,id)  = 0.
       sfnl(is,id) = 0.
     ENDDO
   ENDDO
!
!  *** Calculate prop. constant.                           ***
!  *** Calculate factor R(X) to calculate the NL wave-wave ***
!  *** interaction for shallow water                       ***
!  *** SNLC1 = 1/GRAV**4                                   ***         
!
   snlcs1 = pquad(3)                                                   
   snlcs2 = pquad(4)                                                   
   snlcs3 = pquad(5)                                                   
!   X      = MAX ( 0.75 * DEP2(IGC) * KMESPC , 0.5 )               
   x      = max( 0.75 * dep2(kcgrd(1)) * kmespc , 0.5 )               
   x2     = max( -1.e15, snlcs3*x)
   cons   = snlc1 * ( 1. + snlcs1/x * (1.-snlcs2*x) * exp(x2))
   jacobi = 2. * pi_w
!
!  *** extend the area with action density at periodic boundaries ***
!
   DO iddum = idlow, idhgh
     id = mod( iddum - 1 + mdc , mdc ) + 1
     DO is=1, msc
!       UE (IS,IDDUM) = AC2(ID,IS,IGC) * SPCSIG(IS) * JACOBI       
       ue(is,iddum) = ac2(id,is,kcgrd(1)) * spcsig(is) * jacobi       
     ENDDO
   ENDDO
!
   DO is = msc+1, ishgh
     DO id = idlow, idhgh
       ue(is,id) = ue(is-1,id) * fachfr
     ENDDO
   ENDDO
!
!  *** Calculate (unfolded) interactions ***
!  *** Energy at interacting bins        ***
!
   DO is = isclw, ischg
     DO id = 1, mdc
       e00    =        ue(is      ,id      )
       ep1    = awg1 * ue(is+isp1,id+idp1) +                       &
                awg2 * ue(is+isp1,id+idp ) +                       &
        awg3 * ue(is+isp ,id+idp1) +                       &
        awg4 * ue(is+isp ,id+idp )
       em1    = awg5 * ue(is+ism1,id-idm1) +                       &
                awg6 * ue(is+ism1,id-idm ) +                       &
        awg7 * ue(is+ism ,id-idm1) +                       &
        awg8 * ue(is+ism ,id-idm )
       ep2    = awg1 * ue(is+isp1,id-idp1) +                       &
                awg2 * ue(is+isp1,id-idp ) +                       &
        awg3 * ue(is+isp ,id-idp1) +                       &
        awg4 * ue(is+isp ,id-idp )
       em2    = awg5 * ue(is+ism1,id+idm1) +                       &
                awg6 * ue(is+ism1,id+idm ) +                       &
        awg7 * ue(is+ism ,id+idm1) +                       &
        awg8 * ue(is+ism ,id+idm )
!
!      Contribution to interactions
!
       factor = cons * af11(is) * e00
!
       sa1a   = e00 * ( ep1*dal1 + em1*dal2 ) * pquad(2)               
       sa1b   = sa1a - ep1*em1*dal3 * pquad(2)                         
       sa2a   = e00 * ( ep2*dal1 + em2*dal2 ) * pquad(2)               
       sa2b   = sa2a - ep2*em2*dal3 * pquad(2)                         
!
       sa1(is,id) = factor * sa1b
       sa2(is,id) = factor * sa2b
!
!       IF(ITEST >= 100 .AND. TESTFL)THEN
!         WRITE(PRINTF,9002) E00,EP1,EM1,EP2,EM2
!9002     FORMAT (' E00 EP1 EM1 EP2 EM2  :',5E11.4)
!         WRITE(PRINTF,9003) SA1A,SA1B,SA2A,SA2B
!9003     FORMAT (' SA1A SA1B SA2A SA2B  :',4E11.4)
!         WRITE(PRINTF,9004) IS,ID,SA1(IS,ID),SA2(IS,ID)
!9004     FORMAT (' IS ID SA1() SA2()    :',2I4,2E12.4)
!         WRITE(PRINTF,9005) FACTOR,JACOBI
!9005     FORMAT (' FACTOR JACOBI        : ',2E12.4)
!       END IF
!
     ENDDO
   ENDDO
!
!  *** Fold interactions to side angles -> domain in theta is ***
!  *** periodic                                               ***
!
   DO id = 1, idhgh - mdc
     id0   = 1 - id
     DO is = isclw, ischg
       sa1(is,mdc+id) = sa1(is,  id   )
       sa2(is,mdc+id) = sa2(is,  id   )
       sa1(is,  id0 ) = sa1(is,mdc+id0)
       sa2(is,  id0 ) = sa2(is,mdc+id0)
     ENDDO
   ENDDO
!
!  *** Put source term together (To save space I=IS and ***
!  *** J=MDC is used)  ----                             ***
!
   DO i = 1, msc
     sigpi = spcsig(i) * jacobi                                        
     DO j = 1, mdc
       sfnl(i,j) = - 2. * ( sa1(i,j) + sa2(i,j) )                       &
                   + awg1 * ( sa1(i-isp1,j-idp1) + sa2(i-isp1,j+idp1) ) &
           + awg2 * ( sa1(i-isp1,j-idp ) + sa2(i-isp1,j+idp ) ) &
           + awg3 * ( sa1(i-isp ,j-idp1) + sa2(i-isp ,j+idp1) ) &
           + awg4 * ( sa1(i-isp ,j-idp ) + sa2(i-isp ,j+idp ) ) &
           + awg5 * ( sa1(i-ism1,j+idm1) + sa2(i-ism1,j-idm1) ) &
           + awg6 * ( sa1(i-ism1,j+idm ) + sa2(i-ism1,j-idm ) ) &
           + awg7 * ( sa1(i-ism ,j+idm1) + sa2(i-ism ,j-idm1) ) &
           + awg8 * ( sa1(i-ism ,j+idm ) + sa2(i-ism ,j-idm ) )
!
!      *** store value in auxiliary array and use values in ***
!      *** next four sweeps (see subroutine FILNL3)         ***
!
!       MEMNL4(J,I,IGC) = SFNL(I,J) / SIGPI                        
       memnl4(j,i,kcgrd(1)) = sfnl(i,j) / sigpi                        
     ENDDO
   ENDDO
!
!  *** test output ***
!
!   IF(ITEST >= 50 .AND. TESTFL)THEN
!     WRITE(PRINTF,*)
!     WRITE(PRINTF,*) ' SWSNL3 subroutine '
!     WRITE(PRINTF,9011) IDP, IDP1, IDM, IDM1
!9011 FORMAT (' IDP IDP1 IDM IDM1     :',4I5)
!     WRITE (PRINTF,9013) ISP, ISP1, ISM, ISM1
!9013 FORMAT (' ISP ISP1 ISM ISM1     :',4I5)
!     WRITE (PRINTF,9015) ISLOW, ISHGH, IDLOW, IDHGH
!9015 FORMAT (' ISLOW ISHG IDLOW IDHG :',4I5)
!     WRITE(PRINTF,9016) ISCLW, ISCHG, JACOBI
!9016 FORMAT (' ICLW ICHG JACOBI      :',2I5,E12.4)
!     WRITE (PRINTF,9017) AWG1, AWG2, AWG3, AWG4
!9017 FORMAT (' AWG1 AWG2 AWG3 AWG4   :',4E12.4)
!     WRITE (PRINTF,9018) AWG5, AWG6, AWG7, AWG8
!9018 FORMAT (' AWG5 AWG6 AWG7 AWG8   :',4E12.4)
!     WRITE (PRINTF,9019) MSC4MI, MSC4MA, MDC4MI, MDC4MA
!9019 FORMAT (' S4MI S4MA D4MI D4MA   :',4I6)
!     WRITE(PRINTF,9020) SNLC1,X,X2,CONS
!9020 FORMAT (' SNLC1  X  X2  CONS    :',4E12.4)
!    WRITE(PRINTF,9021) DEP2(KCGRD(1)),KMESPC,FACHFR,PI
!9021 FORMAT (' DEPTH KMESPC FACHFR PI:',4E12.4)
!     WRITE(PRINTF,*)
!
!    *** value source term in every bin ***
!
!     IF(ITEST >= 150)THEN
!       DO I=1, MSC
!         DO J=1, MDC
!           WRITE(PRINTF,2006) I,J,MEMNL4(J,I,KCGRD(1)),SFNL(I,J),      &
!                         SPCSIG(I)                                
!2006       FORMAT (' I J MEMNL() SFNL() SPCSIG:',2I4,3E12.4)           
!         ENDDO
!       ENDDO
!     END IF
!   END IF
!
   RETURN
!

swsnl4()

subroutine swsnl4	(	integer, dimension(*)	WWINT,
		real, dimension(*)	WWAWG,
		real, dimension(msc)	SPCSIG,
		real	SNLC1,
		real	DAL1,
		real	DAL2,
		real	DAL3,
		real, dimension(mt)	DEP2,
		real	KMESPC,
		real, dimension(mdc,msc,mt)	MEMNL4,
		real	FACHFR,
		integer	IDIA,
			ITER
	)

Definition at line 1143 of file swancom4.f90.

!  (This subroutine has not been tested yet)              
!
!*******************************************************************
!
   USE swcomm3                                                         
   USE swcomm4                                                         
   USE ocpcomm4                                                        
   USE m_snl4
!   USE ALL_VARS, ONLY : MT,AC2
   USE vars_wave, ONLY : mt,ac2
!
!   --|-----------------------------------------------------------|--
!     | Delft University of Technology                            |
!     | Faculty of Civil Engineering                              |
!     | Fluid Mechanics Section                                   |
!     | P.O. Box 5048, 2600 GA  Delft, The Netherlands            |
!     |                                                           |
!     | Programmers: H.L. Tolman, R.C. Ris                        |
!   --|-----------------------------------------------------------|--
!
!
!     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
!
!     40.17: IJsbrand Haagsma
!     40.41: Marcel Zijlema
!
!  1. Updates
!
!     40.17, Dec. 01: New Subroutine based on SWSNL3
!     40.41, Oct. 04: common blocks replaced by modules, include files removed
!
!  2. Purpose
!
!     Calculate non-linear interaction using the discrete interaction
!     approximation (Hasselmann and Hasselmann 1985; WAMDI group 1988)
!     for the full circle (option if a current is present). Note: using
!     this subroutine requires an additional array with size
!     (MXC*MYC*MDC*MSC). This requires more internal memory but can
!     speed up the computations sigificantly if a current is present.
!
!  3. Method
!
!     Discrete interaction approximation. To make interpolation simple,
!     the interactions are calculated in a "folded" space.
!
!                            Frequencies -->
!                 +---+---------------------+---------+- IDHGH
!              d  | 3 :          2          :    2    |
!              i  + - + - - - - - - - - - - + - - - - +- MDC
!              r  |   :                     :         |
!              e  | 3 :  original spectrum  :    1    |
!              c  |   :                     :         |
!              t. + - + - - - - - - - - - - + - - - - +- 1
!                 | 3 :          2          :    2    |
!                 +---+---------------------+---------+- IDLOW
!                 |   |                     |     ^   |
!              ISLOW  1                    MSC    |   ISHGH
!                     |                           |
!                   ISCLW                        ISCHG
!              lowest discrete               highest discrete
!                central bin                   central bin
!
!                            1 : Extra tail added beyond MSC
!                            2 : Spectrum copied outside ID range
!                            3 : Empty bins at low frequencies
!
!     ISLOW =  1  + ISM1
!     ISHGH = MSC + ISP1 - ISM1
!     ISCLW =  1
!     ISCHG = MSC - ISM1
!     IDLOW =  1  - MAX(IDM1,IDP1)
!     IDHGH = MDC + MAX(IDM1,IDP1)
!
!       Relative offsets of interpolation points around central bin
!       "#" and corresponding numbers of AWGn :
!
!               ISM1  ISM
!                5        7    T |
!          IDM1   +------+     H +
!                 |      |     E |      ISP      ISP1
!                 |   \  |     T |       3           1
!           IDM   +------+     A +        +---------+  IDP1
!                6       \8      |        |         |
!                                |        |  /      |
!                           \    +        +---------+  IDP
!                                |      /4           2
!                              \ |  /
!          -+-----+------+-------#--------+---------+----------+
!                                |           FREQ.
!
!
!  4. Argument variables
!
!     MCGRD : number of wet grid points of the computational grid
!     MDC   : grid points in theta-direction of computational grid
!     MDC4MA: highest array counter in directional space (Snl4)
!     MDC4MI: lowest array counter in directional space (Snl4)
!     MSC   : grid points in sigma-direction of computational grid
!     MSC4MA: highest array counter in frequency space (Snl4)
!     MSC4MI: lowest array counter in frequency space (Snl4)
!     WWINT : counters for quadruplet interactions
!
   INTEGER WWINT(*)
   INTEGER IDIA
!
!     AC2   : action density
!     AF11  : scaling frequency
!     DAL1  : coefficient for the quadruplet interactions
!     DAL2  : coefficient for the quadruplet interactions
!     DAL3  : coefficient for the quadruplet interactions
!     DEP2  : depth
!     FACHFR
!     KMESPC: mean average wavenumber over full spectrum
!     MEMNL4
!     PI    : circular constant
!     SA1   : interaction contribution of first quadruplet (unfolded space)
!     SA2   : interaction contribution of second quadruplet (unfolded space)
!     SFNL
!     SNLC1
!     SPCSIG: relative frequencies in computational domain in sigma-space
!     UE    : "unfolded" spectrum
!     WWAWG : weight coefficients for the quadruplet interactions
!
   REAL    DAL1, DAL2, DAL3, FACHFR, KMESPC, SNLC1
!   REAL    AC2(MDC,MSC,0:MT)
   REAL    DEP2(MT)
   REAL    MEMNL4(MDC,MSC,MT)
   REAL    SPCSIG(MSC)
   REAL    WWAWG(*)
!
!  6. Local variables
!
   REAL, ALLOCATABLE :: SA1(:,:), SA2(:,:), SFNL(:,:), UE(:,:)
!
!  7. Common blocks used
!
!
!  8. Subroutines used
!
!     ---
!
!  9. Subroutines calling
!
!     SOURCE (in SWANCOM1)
!
! 12. Structure
!
!     -------------------------------------------
!       Initialisations.
!       Calculate proportionality constant.
!       Prepare auxiliary spectrum.
!       Calculate (unfolded) interactions :
!       -----------------------------------------
!         Energy at interacting bins
!         Contribution to interactions
!         Fold interactions to side angles
!       -----------------------------------------
!       Put source term together
!     -------------------------------------------
!
! 13. Source text
!
!*******************************************************************
!
   INTEGER   IS      ,ID      ,ID0     ,I       ,J       ,          &
             ISHGH   ,IDLOW   ,IDHGH   ,ISP     ,ISP1    ,          &
         IDP     ,IDP1    ,ISM     ,ISM1    ,IDM     ,IDM1    , &
         ISCLW   ,ISCHG
!
   REAL      X       ,X2      ,CONS    ,FACTOR  ,SNLCS2  ,          &
             SNLCS3  ,E00     ,EP1     ,EM1     ,EP2     ,EM2     , &
         SA1A    ,SA1B    ,SA2A    ,SA2B    ,                   &
         AWG1    ,AWG2    ,AWG3    ,AWG4    ,AWG5    ,AWG6    , &
         AWG7    ,AWG8    ,                                     &
         JACOBI  ,SIGPI
!
!   SAVE IENT
!   DATA IENT/0/
!   IF (LTRACE) CALL STRACE (IENT,'SWSNL4')

   ALLOCATE(sa1(msc4mi:msc4ma,mdc4mi:mdc4ma))
   ALLOCATE(sa2(msc4mi:msc4ma,mdc4mi:mdc4ma))
   ALLOCATE(sfnl(msc4mi:msc4ma,mdc4mi:mdc4ma))
   ALLOCATE(ue(msc4mi:msc4ma,mdc4mi:mdc4ma))
!
   idp    = wwint(1)
   idp1   = wwint(2)
   idm    = wwint(3)
   idm1   = wwint(4)
   isp    = wwint(5)
   isp1   = wwint(6)
   ism    = wwint(7)
   ism1   = wwint(8)
   islow  = wwint(9)
   ishgh  = wwint(10)
   isclw  = wwint(11)
   ischg  = wwint(12)
   idlow  = wwint(13)
   idhgh  = wwint(14)
!
   awg1 = wwawg(1)
   awg2 = wwawg(2)
   awg3 = wwawg(3)
   awg4 = wwawg(4)
   awg5 = wwawg(5)
   awg6 = wwawg(6)
   awg7 = wwawg(7)
   awg8 = wwawg(8)
!
!  *** Initialize auxiliary arrays per gridpoint ***
!
   DO id = mdc4mi, mdc4ma
     DO is = msc4mi, msc4ma
       ue(is,id)   = 0.
       sa1(is,id)  = 0.
       sa2(is,id)  = 0.
       sfnl(is,id) = 0.
     ENDDO
   ENDDO
!
!  *** Calculate prop. constant.                           ***
!  *** Calculate factor R(X) to calculate the NL wave-wave ***
!  *** interaction for shallow water                       ***
!  *** SNLC1 = 1/GRAV**4                                   ***
!
   snlcs1 = pquad(3)                                                   
   snlcs2 = pquad(4)                                                   
   snlcs3 = pquad(5)                                                   
!   X      = MAX ( 0.75 * DEP2(IGC) * KMESPC , 0.5 )
   x      = max( 0.75 * dep2(kcgrd(1)) * kmespc , 0.5 )
   x2     = max( -1.e15, snlcs3*x)
   cons   = snlc1 * ( 1. + snlcs1/x * (1.-snlcs2*x) * exp(x2))
   jacobi = 2. * pi_w
!
!  *** extend the area with action density at periodic boundaries ***
!
   DO iddum = idlow, idhgh
     id = mod( iddum - 1 + mdc , mdc ) + 1
     DO is=1, msc
!       UE (IS,IDDUM) = AC2(ID,IS,IGC) * SPCSIG(IS) * JACOBI
       ue(is,iddum) = ac2(id,is,kcgrd(1)) * spcsig(is) * jacobi
     ENDDO
   ENDDO
!
   DO is = msc+1, ishgh
     DO id = idlow, idhgh
       ue(is,id) = ue(is-1,id) * fachfr
     ENDDO
   ENDDO
!
!  *** Calculate (unfolded) interactions ***
!  *** Energy at interacting bins        ***
!
   DO is = isclw, ischg
     DO id = 1, mdc
       e00    =        ue(is      ,id      )
       ep1    = awg1 * ue(is+isp1,id+idp1) +                           &
                awg2 * ue(is+isp1,id+idp ) +                           &
        awg3 * ue(is+isp ,id+idp1) +                           &
        awg4 * ue(is+isp ,id+idp )
       em1    = awg5 * ue(is+ism1,id-idm1) +                           &
                awg6 * ue(is+ism1,id-idm ) +                           &
        awg7 * ue(is+ism ,id-idm1) +                           &
        awg8 * ue(is+ism ,id-idm )
       ep2    = awg1 * ue(is+isp1,id-idp1) +                           &
                awg2 * ue(is+isp1,id-idp ) +                           &
        awg3 * ue(is+isp ,id-idp1) +                           &
        awg4 * ue(is+isp ,id-idp )
       em2    = awg5 * ue(is+ism1,id+idm1) +                           &
                awg6 * ue(is+ism1,id+idm ) +                           &
        awg7 * ue(is+ism ,id+idm1) +                           &
        awg8 * ue(is+ism ,id+idm )
!
!      Contribution to interactions
!
       factor = cons * af11(is) * e00
!
       sa1a   = e00 * ( ep1*dal1 + em1*dal2 ) * cnl4_1(idia)
       sa1b   = sa1a - ep1*em1*dal3 * cnl4_2(idia)
       sa2a   = e00 * ( ep2*dal1 + em2*dal2 ) * cnl4_1(idia)
       sa2b   = sa2a - ep2*em2*dal3 * cnl4_2(idia)
!

       sa1(is,id) = factor * sa1b
       sa2(is,id) = factor * sa2b
!
!       IF(ITEST >= 100 .AND. TESTFL)THEN
!         WRITE(PRINTF,9002) E00,EP1,EM1,EP2,EM2
!9002     FORMAT (' E00 EP1 EM1 EP2 EM2  :',5E11.4)
!         WRITE(PRINTF,9003) SA1A,SA1B,SA2A,SA2B
!9003     FORMAT (' SA1A SA1B SA2A SA2B  :',4E11.4)
!         WRITE(PRINTF,9004) IS,ID,SA1(IS,ID),SA2(IS,ID)
!9004     FORMAT (' IS ID SA1() SA2()    :',2I4,2E12.4)
!         WRITE(PRINTF,9005) FACTOR,JACOBI
!9005     FORMAT (' FACTOR JACOBI        : ',2E12.4)
!       END IF
!
     ENDDO
   ENDDO
!
!  *** Fold interactions to side angles -> domain in theta is ***
!  *** periodic                                               ***
!
   DO id = 1, idhgh - mdc
     id0   = 1 - id
     DO is = isclw, ischg
       sa1(is,mdc+id) = sa1(is,  id   )
       sa2(is,mdc+id) = sa2(is,  id   )
       sa1(is,  id0 ) = sa1(is,mdc+id0)
       sa2(is,  id0 ) = sa2(is,mdc+id0)
     ENDDO
   ENDDO
!
!  *** Put source term together (To save space I=IS and ***
!  *** J=MDC is used)                                   ***
!
   fac = 1.                                                            

   DO i = 1, msc
     sigpi = spcsig(i) * jacobi                                        
     DO j = 1, mdc
       sfnl(i,j) = - 2. * ( sa1(i,j) + sa2(i,j) )                        &
                   + awg1 * ( sa1(i-isp1,j-idp1) + sa2(i-isp1,j+idp1) )  &
           + awg2 * ( sa1(i-isp1,j-idp ) + sa2(i-isp1,j+idp ) )  &
           + awg3 * ( sa1(i-isp ,j-idp1) + sa2(i-isp ,j+idp1) )  &
           + awg4 * ( sa1(i-isp ,j-idp ) + sa2(i-isp ,j+idp ) )  &
           + awg5 * ( sa1(i-ism1,j+idm1) + sa2(i-ism1,j-idm1) )  &
           + awg6 * ( sa1(i-ism1,j+idm ) + sa2(i-ism1,j-idm ) )  &
           + awg7 * ( sa1(i-ism ,j+idm1) + sa2(i-ism ,j-idm1) )  &
           + awg8 * ( sa1(i-ism ,j+idm ) + sa2(i-ism ,j-idm ) )
!
!      *** store value in auxiliary array and use values in ***
!      *** next four sweeps (see subroutine FILNL3)         ***
!
       IF(idia == 1)THEN
!         MEMNL4(J,I,IGC) = FAC * SFNL(I,J) / SIGPI
         memnl4(j,i,kcgrd(1)) = fac * sfnl(i,j) / sigpi
       ELSE
!         MEMNL4(J,I,IGC) = MEMNL4(J,I,IGC) + FAC * SFNL(I,J) / SIGPI
         memnl4(j,i,kcgrd(1)) = memnl4(j,i,kcgrd(1)) + fac * sfnl(i,j) / sigpi
       END IF
     ENDDO
   ENDDO
!
!  *** test output ***
!
!   IF(ITEST >= 50 .AND. TESTFL)THEN
!     WRITE(PRINTF,*)
!     WRITE(PRINTF,*) ' SWSNL4 subroutine '
!     WRITE(PRINTF,9011) IDP, IDP1, IDM, IDM1
!9011 FORMAT (' IDP IDP1 IDM IDM1     :',4I5)
!     WRITE (PRINTF,9013) ISP, ISP1, ISM, ISM1
!9013 FORMAT (' ISP ISP1 ISM ISM1     :',4I5)
!     WRITE (PRINTF,9015) ISLOW, ISHGH, IDLOW, IDHGH
!9015 FORMAT (' ISLOW ISHG IDLOW IDHG :',4I5)
!     WRITE(PRINTF,9016) ISCLW, ISCHG, JACOBI
!9016 FORMAT (' ICLW ICHG JACOBI      :',2I5,E12.4)
!     WRITE (PRINTF,9017) AWG1, AWG2, AWG3, AWG4
!9017 FORMAT (' AWG1 AWG2 AWG3 AWG4   :',4E12.4)
!     WRITE (PRINTF,9018) AWG5, AWG6, AWG7, AWG8
!9018 FORMAT (' AWG5 AWG6 AWG7 AWG8   :',4E12.4)
!     WRITE (PRINTF,9019) MSC4MI, MSC4MA, MDC4MI, MDC4MA
!9019 FORMAT (' S4MI S4MA D4MI D4MA   :',4I6)
!     WRITE(PRINTF,9020) SNLC1,X,X2,CONS
!9020 FORMAT (' SNLC1  X  X2  CONS    :',4E12.4)
!     WRITE(PRINTF,9021) DEP2(KCGRD(1)),KMESPC,FACHFR,PI
!9021 FORMAT (' DEPTH KMESPC FACHFR PI:',4E12.4)
!     WRITE(PRINTF,*)
!
!    *** value source term in every bin ***
!
!     IF(ITEST >= 150)THEN
!       DO I=1, MSC
!         DO J=1, MDC
!           WRITE(PRINTF,2006) I,J,MEMNL4(J,I,KCGRD(1)),SFNL(I,J),     &
!                         SPCSIG(I)
!2006       FORMAT (' I J MEMNL() SFNL() SPCSIG:',2I4,3E12.4)
!         ENDDO
!       ENDDO
!     END IF
!   END IF
!
   DEALLOCATE (sa1, sa2, sfnl, ue)

   RETURN
!

swsnl8()

subroutine swsnl8	(	integer, dimension(*)	WWINT,
		real, dimension(msc4mi:msc4ma,mdc4mi:mdc4ma)	UE,
		real, dimension(msc4mi:msc4ma,mdc4mi:mdc4ma)	SA1,
		real, dimension(msc4mi:msc4ma,mdc4mi:mdc4ma)	SA2,
		real, dimension(msc)	SPCSIG,
		real	SNLC1,
		real	DAL1,
		real	DAL2,
		real	DAL3,
		real, dimension(msc4mi:msc4ma,mdc4mi:mdc4ma)	SFNL,
		real, dimension(mt)	DEP2,
		real	KMESPC,
		real, dimension(mdc,msc,mt)	MEMNL4,
		real	FACHFR
	)

Definition at line 665 of file swancom4.f90.

              
!  (This subroutine has not been tested yet)              
!*********************************************************************
!
   USE swcomm3                                                         
   USE swcomm4                                                         
   USE ocpcomm4                                                        
   USE m_snl4
!   USE ALL_VARS, ONLY : MT,AC2
   USE vars_wave, ONLY : mt,ac2
!
   IMPLICIT NONE
!
!  2. Purpose
!
!     Calculate non-linear interaction using the discrete interaction
!     approximation (Hasselmann and Hasselmann 1985; WAMDI group 1988)
!     for the full circle.
!
!  3. Method
!
!     Discrete interaction approximation. To make interpolation simple,
!     the interactions are calculated in a "folded" space.
!
!                            Frequencies -->
!                 +---+---------------------+---------+- IDHGH
!              d  | 3 :          2          :    2    |
!              i  + - + - - - - - - - - - - + - - - - +- MDC
!              r  |   :                     :         |
!              e  | 3 :  original spectrum  :    1    |
!              c  |   :                     :         |
!              t. + - + - - - - - - - - - - + - - - - +- 1
!                 | 3 :          2          :    2    |
!                 +---+---------------------+---------+- IDLOW
!                 |   |                     |     ^   |
!              ISLOW  1                    MSC    |   ISHGH
!                     |                           |
!                   ISCLW                        ISCHG
!              lowest discrete               highest discrete
!                central bin                   central bin
!
!                            1 : Extra tail added beyond MSC
!                            2 : Spectrum copied outside ID range
!                            3 : Empty bins at low frequencies
!
!     ISLOW =  1  + ISM1
!     ISHGH = MSC + ISP1 - ISM1
!     ISCLW =  1
!     ISCHG = MSC - ISM1
!     IDLOW =  1  - MAX(IDM1,IDP1)
!     IDHGH = MDC + MAX(IDM1,IDP1)
!
!     Note: using this subroutine requires an additional array
!           with size MXC*MYC*MDC*MSC.
!
   INTEGER WWINT(*)
!
   REAL    DAL1, DAL2, DAL3, FACHFR, KMESPC, SNLC1
!   REAL    AC2(MDC,MSC,0:MT)
   REAL    DEP2(MT)
   REAL    MEMNL4(MDC,MSC,MT)
   REAL    SA1(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
   REAL    SA2(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
   REAL    SFNL(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
   REAL    SPCSIG(MSC)
   REAL    UE(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
!
!*******************************************************************
!
   INTEGER   IS      ,ID      ,ID0     ,I       ,J       ,         &
             ISLOW   ,ISHGH   ,IDLOW   ,IDHGH   ,                  &
         IDPP    ,IDMM    ,ISPP    ,ISMM    ,                  &
         ISCLW   ,ISCHG   ,IDDUM   ,IENT
!
   REAL      X       ,X2      ,CONS    ,FACTOR  ,SNLCS1  ,SNLCS2  ,  &
             SNLCS3  ,E00     ,EP1     ,EM1     ,EP2     ,EM2     ,  &
         SA1A    ,SA1B    ,SA2A    ,SA2B    ,                    &
         JACOBI  ,SIGPI
!
!   SAVE IENT
!   DATA IENT/0/
!   IF (LTRACE) CALL STRACE (IENT,'SWSNL8')
!
   islow  = wwint(9)
   ishgh  = wwint(10)
   isclw  = wwint(11)
   ischg  = wwint(12)
   idlow  = wwint(13)
   idhgh  = wwint(14)
   idpp   = wwint(21)
   idmm   = wwint(22)
   ispp   = wwint(23)
   ismm   = wwint(24)
!
!  *** Calculate prop. constant.                           ***
!  *** Calculate factor R(X) to calculate the NL wave-wave ***
!  *** interaction for shallow water                       ***
!  *** SNLC1 = 1/GRAV**4                                   ***
!
   snlcs1 = pquad(3)
   snlcs2 = pquad(4)
   snlcs3 = pquad(5)
!   X      = MAX ( 0.75 * DEP2(IGC) * KMESPC , 0.5 )
   x      = max( 0.75 * dep2(kcgrd(1)) * kmespc , 0.5 )
   x2     = max( -1.e15, snlcs3*x)
   cons   = snlc1 * ( 1. + snlcs1/x * (1.-snlcs2*x) * exp(x2))
   jacobi = 2. * pi_w
!
!  *** extend the area with action density at periodic boundaries ***
!
   DO iddum = idlow, idhgh
     id = mod( iddum - 1 + mdc , mdc ) + 1
     DO is=1, msc
!       UE (IS,IDDUM) = AC2(ID,IS,IGC) * SPCSIG(IS) * JACOBI
       ue(is,iddum) = ac2(id,is,kcgrd(1)) * spcsig(is) * jacobi
     ENDDO
   ENDDO
!
   DO id = idlow, idhgh
     DO is = msc+1, ishgh
       ue(is,id) = ue(is-1,id) * fachfr
     ENDDO
   ENDDO
!
!  *** Calculate (unfolded) interactions ***
!  *** Energy at interacting bins        ***
!
   DO id = 1, mdc
     DO is = isclw, ischg
       e00 = ue(is     ,id     )
       ep1 = ue(is+ispp,id+idpp)
       em1 = ue(is+ismm,id-idmm)
       ep2 = ue(is+ispp,id-idpp)
       em2 = ue(is+ismm,id+idmm)
!
!      Contribution to interactions
!
       factor = cons * af11(is) * pquad(2) * e00
!
       sa1a   = e00 * ( ep1*dal1 + em1*dal2 )
       sa1b   = sa1a - ep1*em1*dal3
       sa2a   = e00 * ( ep2*dal1 + em2*dal2 )
       sa2b   = sa2a - ep2*em2*dal3
!
       sa1(is,id) = factor * sa1b
       sa2(is,id) = factor * sa2b
!
     ENDDO
   ENDDO
!
!  *** Fold interactions to side angles -> domain in theta is ***
!  *** periodic                                               ***
!
   DO id = 1, idhgh - mdc
     id0   = 1 - id
     DO is = isclw, ischg
       sa1(is,mdc+id) = sa1(is,  id   )
       sa2(is,mdc+id) = sa2(is,  id   )
       sa1(is,  id0 ) = sa1(is,mdc+id0)
       sa2(is,  id0 ) = sa2(is,mdc+id0)
     ENDDO
   ENDDO
!
!  *** Put source term together (To save space I=IS and ***
!  *** J=MDC is used)  ----                             ***
!
   DO i = 1, msc
     sigpi = spcsig(i) * jacobi
     DO j = 1, mdc
       sfnl(i,j) =   - 2. * ( sa1(i,j) + sa2(i,j) )                 &
                   + ( sa1(i-ispp,j-idpp) + sa2(i-ispp,j+idpp) )    &
           + ( sa1(i-ismm,j+idmm) + sa2(i-ismm,j-idmm) )
!
!      *** store value in auxiliary array and use values in ***
!      *** next four sweeps (see subroutine FILNL3)         ***
!
!       MEMNL4(J,I,IGC) = SFNL(I,J) / SIGPI
       memnl4(j,i,kcgrd(1)) = sfnl(i,j) / sigpi
     ENDDO
   ENDDO
!
!  *** value source term in every bin ***
!
!   IF(ITEST >= 150 .AND. TESTFL)THEN
!     DO I=1, MSC
!       DO J=1, MDC
!         WRITE(PRINTF,2006) I,J,MEMNL4(J,I,KCGRD(1)),SFNL(I,J),SPCSIG(I)
!2006     FORMAT (' I J MEMNL() SFNL() SPCSIG:',2I4,3E12.4)
!       ENDDO
!     ENDDO
!   END IF
!
   RETURN
!

wave()

real function wave

(

D

)

Definition at line 2717 of file swancom4.f90.

!  Transforms nondimensional Sqr(w)d/g into nondimensional kd using the
!  dispersion relation and an iterative method

   IF(d-10.) 2,2,1
 1 xx=d
   GO TO 6
 2 IF(d-1.) 3,4,4
 3 x=sqrt(d)
   GO TO 5
 4 x=d
 5 cothx=1.0/tanh(x)
   xx=x-(x-d*cothx)/(1.+d*(cothx**2-1.))
   e=1.-xx/x
   x=xx
   IF(abs(e)-0.0005) 6,5,5
 6 wave=xx
   RETURN

Here is the caller graph for this function: