My Project
swancom4.f90
Go to the documentation of this file.
1 
2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 !
13 !******************************************************************
14 !
15  SUBROUTINE fac4ww (XIS ,SNLC1 ,DAL1 ,DAL2 ,DAL3 ,SPCSIG, &
16  WWINT ,WWAWG ,WWSWG )
17 !(This subroutine has not been tested yet)
18 !
19 !******************************************************************
20 
21 !
22  USE swcomm3
23  USE swcomm4
24  USE ocpcomm4
25  USE m_snl4
26 !
27 ! Purpose :
28 !
29 ! Calculate interpolation constants for Snl.
30 !
31  REAL SPCSIG(MSC)
32  INTEGER MSC2 ,MSC1 ,IS ,IDP ,IDP1 , &
33  IDM ,IDM1 ,ISP ,ISP1 ,ISM ,ISM1 , &
34  IDPP ,IDMM ,ISPP ,ISMM , &
35  ISLOW ,ISHGH ,ISCLW ,ISCHG ,IDLOW ,IDHGH , &
36  MSCMAX,MDCMAX
37 
38  REAL SNLC1 ,LAMM2 ,LAMP2 ,DELTH3, &
39  AUX1 ,DELTH4,DAL1 ,DAL2 ,DAL3 ,CIDP ,WIDP , &
40  WIDP1 ,CIDM ,WIDM ,WIDM1 ,XIS ,XISLN ,WISP , &
41  WISP1 ,WISM ,WISM1 ,AWG1 ,AWG2 ,AWG3 ,AWG4 , &
42  AWG5 ,AWG6 ,AWG7 ,AWG8 ,SWG1 ,SWG2 ,SWG3 , &
43  SWG4 ,SWG5 ,SWG6 ,SWG7 ,SWG8 ,FREQ , &
44  RADE
45 !
46  REAL WWAWG(*),WWSWG(*)
47 !
48  INTEGER WWINT(*)
49 !
50 
51  IF(ALLOCATED(af11)) DEALLOCATE(af11)
52 
53 ! *** Compute frequency indices ***
54 ! *** XIS is the relative increment of the relative frequency ***
55 !
56  msc2 = int( float(msc) / 2.0 )
57  msc1 = msc2 - 1
58 
59  xis = spcsig(msc2) / spcsig(msc1)
60 
61 !
62 ! *** set values for the nonlinear four-wave interactions ***
63 !
64  snlc1 = 1. / grav_w**4
65 !
66  lamm2 = (1.-pquad(1))**2
67  lamp2 = (1.+pquad(1))**2
68  delth3 = acos( (lamm2**2+4.-lamp2**2) / (4.*lamm2) )
69  aux1 = sin(delth3)
70  delth4 = asin(-aux1*lamm2/lamp2)
71 !
72  dal1 = 1. / (1.+pquad(1))**4
73  dal2 = 1. / (1.-pquad(1))**4
74  dal3 = 2. * dal1 * dal2
75 !
76 ! *** Compute directional indices in sigma and theta space ***
77 !
78  cidp = abs(delth4/ddir)
79  idp = int(cidp)
80  idp1 = idp + 1
81  widp = cidp - real(idp)
82  widp1 = 1.- widp
83 !
84  cidm = abs(delth3/ddir)
85  idm = int(cidm)
86  idm1 = idm + 1
87  widm = cidm - real(idm)
88  widm1 = 1.- widm
89  xisln = log( xis )
90 !
91  isp = int( log(1.+pquad(1)) / xisln )
92  isp1 = isp + 1
93  wisp = (1.+pquad(1) - xis**isp) / (xis**isp1 - xis**isp)
94  wisp1 = 1. - wisp
95 !
96  ism = int( log(1.-pquad(1)) / xisln )
97  ism1 = ism - 1
98  wism = (xis**ism -(1.-pquad(1))) / (xis**ism - xis**ism1)
99  wism1 = 1. - wism
100 !
101 ! *** Range of calculations ***
102 !
103  islow = 1 + ism1
104  ishgh = msc + isp1 - ism1
105  isclw = 1
106  ischg = msc - ism1
107  idlow = 1 - mdc - max(idm1,idp1)
108  idhgh = mdc + mdc + max(idm1,idp1)
109 !
110  msc4mi = islow
111  msc4ma = ishgh
112  mdc4mi = idlow
113  mdc4ma = idhgh
114  mscmax = msc4ma - msc4mi + 1
115  mdcmax = mdc4ma - mdc4mi + 1
116 !
117 ! *** Interpolation weights ***
118 !
119  awg1 = widp * wisp
120  awg2 = widp1 * wisp
121  awg3 = widp * wisp1
122  awg4 = widp1 * wisp1
123 !
124  awg5 = widm * wism
125  awg6 = widm1 * wism
126  awg7 = widm * wism1
127  awg8 = widm1 * wism1
128 !
129 ! *** quadratic interpolation ***
130 !
131  swg1 = awg1**2
132  swg2 = awg2**2
133  swg3 = awg3**2
134  swg4 = awg4**2
135 !
136  swg5 = awg5**2
137  swg6 = awg6**2
138  swg7 = awg7**2
139  swg8 = awg8**2
140 !
141 ! --- determine discrete counters for piecewise
142 ! constant interpolation
143 !
144  IF(awg1 < awg2)THEN
145  IF(awg2 < awg3)THEN
146  IF(awg3 < awg4)THEN
147  ispp=isp
148  idpp=idp
149  ELSE
150  ispp=isp
151  idpp=idp1
152  END IF
153  ELSE IF(awg2 < awg4)THEN
154  ispp=isp
155  idpp=idp
156  ELSE
157  ispp=isp1
158  idpp=idp
159  END IF
160  ELSE IF(awg1 < awg3)THEN
161  IF(awg3 < awg4)THEN
162  ispp=isp
163  idpp=idp
164  ELSE
165  ispp=isp
166  idpp=idp1
167  END IF
168  ELSE IF(awg1 < awg4)THEN
169  ispp=isp
170  idpp=idp
171  ELSE
172  ispp=isp1
173  idpp=idp1
174  END IF
175  IF(awg5 < awg6)THEN
176  IF(awg6 < awg7)THEN
177  IF(awg7 < awg8)THEN
178  ismm=ism
179  idmm=idm
180  ELSE
181  ismm=ism
182  idmm=idm1
183  END IF
184  ELSE IF(awg6 < awg8)THEN
185  ismm=ism
186  idmm=idm
187  ELSE
188  ismm=ism1
189  idmm=idm
190  END IF
191  ELSE IF(awg5 < awg7)THEN
192  IF(awg7 < awg8)THEN
193  ismm=ism
194  idmm=idm
195  ELSE
196  ismm=ism
197  idmm=idm1
198  END IF
199  ELSE IF(awg5 < awg8)THEN
200  ismm=ism
201  idmm=idm
202  ELSE
203  ismm=ism1
204  idmm=idm1
205  END IF
206 !
207 ! *** fill the arrays *
208 !
209  wwint(1) = idp
210  wwint(2) = idp1
211  wwint(3) = idm
212  wwint(4) = idm1
213  wwint(5) = isp
214  wwint(6) = isp1
215  wwint(7) = ism
216  wwint(8) = ism1
217  wwint(9) = islow
218  wwint(10)= ishgh
219  wwint(11)= isclw
220  wwint(12)= ischg
221  wwint(13)= idlow
222  wwint(14)= idhgh
223  wwint(15)= msc4mi
224  wwint(16)= msc4ma
225  wwint(17)= mdc4mi
226  wwint(18)= mdc4ma
227  wwint(19)= mscmax
228  wwint(20)= mdcmax
229  wwint(21)= idpp
230  wwint(22)= idmm
231  wwint(23)= ispp
232  wwint(24)= ismm
233 !
234  wwawg(1) = awg1
235  wwawg(2) = awg2
236  wwawg(3) = awg3
237  wwawg(4) = awg4
238  wwawg(5) = awg5
239  wwawg(6) = awg6
240  wwawg(7) = awg7
241  wwawg(8) = awg8
242 !
243  wwswg(1) = swg1
244  wwswg(2) = swg2
245  wwswg(3) = swg3
246  wwswg(4) = swg4
247  wwswg(5) = swg5
248  wwswg(6) = swg6
249  wwswg(7) = swg7
250  wwswg(8) = swg8
251 
252  ALLOCATE (af11(msc4mi:msc4ma))
253 
254 ! *** Fill scaling array (f**11) ***
255 ! *** compute the radian frequency**11 for IS=1, MSC ***
256 !
257  DO is=1, msc
258  af11(is) = ( spcsig(is) / ( 2. * pi_w ) )**11
259  END DO
260 !
261 ! *** compute the radian frequency for the IS = MSC+1, ISHGH ***
262 !
263  freq = spcsig(msc) / ( 2. * pi_w )
264  DO is = msc+1, ishgh
265  freq = freq * xis
266  af11(is) = freq**11
267  END DO
268 !
269 ! *** compute the radian frequency for IS = 0, ISLOW ***
270 !
271  freq = spcsig(1) / ( 2. * pi_w )
272  DO is = 0, islow, -1
273  freq = freq / xis
274  af11(is) = freq**11
275  END DO
276 
277  RETURN
278  END SUBROUTINE fac4ww
279 
280 !
281 !****************************************************************
282 !
283  SUBROUTINE swlta ( IG ,DEP2 ,CGO ,SPCSIG , &
284  KWAVE ,IMATRA ,IMATDA ,IDDLOW , &
285  IDDTOP ,ISSTOP ,IDCMIN ,IDCMAX , &
286  HS ,SMEBRK ,PLTRI ,URSELL )
287 ! (This subroutine has not been tested yet)
288 !
289 !****************************************************************
290 !
291  USE ocpcomm4
292  USE swcomm3
293  USE swcomm4
294  USE vars_wave, ONLY : mt, ac2
295 !
296  IMPLICIT NONE
297 !
298 ! 2. Purpose
299 !
300 ! In this subroutine the triad-wave interactions are calculated
301 ! with the Lumped Triad Approximation of Eldeberky (1996). His
302 ! expression is based on a parametrization of the biphase (as
303 ! function of the Ursell number), is directionally uncoupled and
304 ! takes into account for the self-self interactions only.
305 !
306 ! For a full description of the equations reference is made
307 ! to PhD thesis of Eldeberky (1996). Here only the main expressions
308 ! are given.
309 !
310 ! 3. Method
311 !
312 ! The parametrized biphase is given by (see eq. 3.19):
313 !
314 ! 0.2
315 ! beta = - pi/2 + pi/2 tanh ( ----- )
316 ! Ur
317 !
318 ! The Ursell number is calculated in routine SINTGRL
319 !
320 ! The source term as function of frequency p is (see eq. 7.25):
321 !
322 ! + -
323 ! S(p) = S(p) + S(p)
324 !
325 ! in which
326 !
327 ! +
328 ! S(p) = alpha Cp Cg,p (R(p/2,p/2))**2 sin (|beta|) ( E(p/2)**2 -2 E(p) E(p/2) )
329 !
330 ! - +
331 ! S(p) = - 2 S(2p)
332 !
333 ! with alpha a tunable coefficient and R(p/2,p/2) is the interaction
334 ! coefficient of which the expression can be found in Eldeberky (1996);
335 ! see eq. 7.26.
336 !
337 ! Note that a slightly adapted formulation of the LTA is used in
338 ! in the SWAN model:
339 !
340 ! - Only positive contributions to higher harmonics are considered
341 ! here (no energy is transferred to lower harmonics).
342 !
343 ! - The mean frequency in the expression of the Ursell number
344 ! is calculated according to the first order moment over the
345 ! zeroth order moment (personal communication, Y.Eldeberky, 1997).
346 !
347 ! - The interactions are calculated up to 2.5 times the mean
348 ! frequency only.
349 !
350 ! - Since the spectral grid is logarithmically distributed in frequency
351 ! space, the interactions between central bin and interacting bin
352 ! are interpolated such that the distance between these bins is
353 ! factor 2 (nearly).
354 !
355 ! - The interactions are calculated in terms of energy density
356 ! instead of action density. So the action density spectrum
357 ! is firstly converted to the energy density grid, then the
358 ! interactions are calculated and then the spectrum is converted
359 ! to the action density spectrum back.
360 !
361 ! - To ensure numerical stability the Patankar rule is used.
362 !
363  INTEGER IDDLOW, IDDTOP, ISSTOP
364  INTEGER IDCMIN(MSC), IDCMAX(MSC)
365 
366  REAL :: HS, SMEBRK
367 ! REAL :: AC2(MDC,MSC,0:MT)
368 
369  REAL :: CGO(MSC,MICMAX)
370  REAL :: DEP2(MT)
371  REAL :: IMATDA(MDC,MSC), IMATRA(MDC,MSC)
372  REAL :: SPCSIG(MSC)
373  REAL :: KWAVE(MSC,MICMAX)
374  REAL :: PLTRI(MDC,MSC,NPTST)
375  REAL :: URSELL(MT)
376 
377  INTEGER I1, I2, ID, IG, IDDUM, IENT, II, IS, ISM, ISM1, ISMAX, ISP, ISP1
378  REAL AUX1, AUX2, BIPH, C0, CM, DEP, DEP_2, DEP_3, E0, EM, &
379  FT, RINT, SIGPI, SINBPH, STRI, WISM, WISM1, WISP, WISP1, &
380  W0, WM, WN0, WNM, XIS, XISLN
381  REAL, ALLOCATABLE :: E(:), SA(:,:)
382 
383 ! DEP = DEP2(IGC)
384  dep = dep2(ig)
385  dep_2 = dep**2
386  dep_3 = dep**3
387 !
388 ! --- compute some indices in sigma space
389 !
390  i2 = int(float(msc) / 2.)
391  i1 = i2 - 1
392  xis = spcsig(i2) / spcsig(i1)
393  xisln = log( xis )
394 
395  isp = int( log(2.) / xisln )
396  isp1 = isp + 1
397  wisp = (2. - xis**isp) / (xis**isp1 - xis**isp)
398  wisp1 = 1. - wisp
399 
400  ism = int( log(0.5) / xisln )
401  ism1 = ism - 1
402  wism = (xis**ism -0.5) / (xis**ism - xis**ism1)
403  wism1 = 1. - wism
404 
405  ALLOCATE (e(1:msc))
406  ALLOCATE (sa(1:mdc,1:msc+isp1))
407  e = 0.
408  sa = 0.
409 !
410 ! --- compute maximum frequency for which interactions are calculated
411 !
412  ismax = 1
413  DO is = 1, msc
414  IF(spcsig(is) < ( ptriad(2) * smebrk) )THEN
415  ismax = is
416  ENDIF
417  ENDDO
418  ismax = max( ismax , isp1 )
419 !
420 ! --- compute 3 wave-wave interactions
421 !
422 ! IF( URSELL(IGC) >= PTRIAD(5) )THEN
423  IF( ursell(ig) >= ptriad(5) )THEN
424 !
425 ! --- calculate biphase
426 !
427 ! BIPH = (0.5*PI_W)*(TANH(PTRIAD(4)/URSELL(IGC))-1.)
428  biph = (0.5*pi_w)*(tanh(ptriad(4)/ursell(ig))-1.)
429  sinbph = abs( sin(biph) )
430 !
431  DO ii = iddlow, iddtop
432  id = mod( ii - 1 + mdc , mdc ) + 1
433 !
434 ! --- initialize array with E(f) for the direction considered
435 !
436  DO is = 1, msc
437 ! E(IS) = AC2(ID,IS,IGC) * 2. * PI_W * SPCSIG(IS)
438  e(is) = ac2(id,is,ig) * 2. * pi_w * spcsig(is)
439  END DO
440 !
441  DO is = 1, ismax
442 
443  e0 = e(is)
444  w0 = spcsig(is)
445  wn0 = kwave(is,1)
446  c0 = w0 / wn0
447 
448  IF( is > -ism1 )THEN
449  em = wism * e(is+ism1) + wism1 * e(is+ism)
450  wm = wism * spcsig(is+ism1) + wism1 * spcsig(is+ism)
451  wnm = wism * kwave(is+ism1,1) + wism1 * kwave(is+ism,1)
452  cm = wm / wnm
453  ELSE
454  em = 0.
455  wm = 0.
456  wnm = 0.
457  cm = 0.
458  END IF
459 
460  aux1 = wnm**2 * ( grav_w * dep + 2.*cm**2 )
461  aux2 = wn0 * dep * ( grav_w * dep + &
462  (2./15.) * grav_w * dep_3 * wn0**2 - &
463  (2./ 5.) * w0**2 * dep_2 )
464  rint = aux1 / aux2
465  ft = ptriad(1) * c0 * cgo(is,1) * rint**2 * sinbph
466 
467  sa(id,is) = max(0., ft * ( em * em - 2. * em * e0 ))
468 
469  END DO
470  END DO
471 !
472 ! --- put source term together
473 !
474  DO is = 1, isstop
475  sigpi = spcsig(is) * 2. * pi_w
476  DO iddum = idcmin(is), idcmax(is)
477  id = mod( iddum - 1 + mdc , mdc ) + 1
478 !
479  stri = sa(id,is) - 2.*(wisp * sa(id,is+isp1) + &
480  wisp1 * sa(id,is+isp ))
481 !
482 ! --- store results in rhs and main diagonal according
483 ! to Patankar-rules
484 !
485  IF(testfl) pltri(id,is,iptst) = stri / sigpi
486  IF(stri > 0.)THEN
487  imatra(id,is) = imatra(id,is) + stri / sigpi
488  ELSE
489  imatda(id,is) = imatda(id,is) - stri / &
490 ! MAX(1.E-18,AC2(ID,IS,IGC)*SIGPI)
491  max(1.e-18,ac2(id,is,ig)*sigpi)
492  END IF
493  END DO
494  END DO
495 
496  END IF
497 
498  DEALLOCATE(e,sa)
499 
500  RETURN
501  END SUBROUTINE swlta
502 
503 !
504 !******************************************************************
505 !
506  SUBROUTINE range4 (WWINT,IDDLOW,IDDTOP)
507 !
508 !******************************************************************
509 !
510 ! calculate the minimum and maximum counters in frequency and
511 ! directional space which fall with the calculation for the
512 ! nonlinear wave-wave interactions.
513 !
514 ! Method : review for the counters :
515 !
516 ! Frequencies -->
517 ! +---+---------------------+---------+- IDHGH
518 ! d | 3 : 2 : 2 |
519 ! i + - + - - - - - - - - - - + - - - - +- MDC
520 ! r | : : |
521 ! e | 3 : original spectrum : 1 |
522 ! c | : : |
523 ! t. + - + - - - - - - - - - - + - - - - +- 1
524 ! | 3 : 2 : 2 |
525 ! +---+---------------------+---------+- IDLOW
526 ! | | | ^ |
527 ! ISLOW 1 MSC | ISHGH
528 ! ^ |
529 ! | |
530 ! ISCLW ISCHG
531 ! lowest discrete highest discrete
532 ! central bin central bin
533 !
534 !
535 ! The directional counters depend on the numerical method that
536 ! is used.
537 !****************************************************************
538 !
539  USE swcomm3
540  USE swcomm4
541  USE ocpcomm4
542  IMPLICIT NONE
543 
544  INTEGER IDDLOW,IDDTOP
545 !
546  INTEGER WWINT(*)
547 !
548 ! *** Range in directional domain ***
549 !
550  IF(iquad < 3 .AND. iquad > 0)THEN
551 ! *** counters based on bins which fall within a sweep ***
552  wwint(13) = iddlow - max( wwint(4), wwint(2) )
553  wwint(14) = iddtop + max( wwint(4), wwint(2) )
554  ELSE
555 ! *** counters initially based on full circle ***
556  wwint(13) = 1 - max( wwint(4), wwint(2) )
557  wwint(14) = mdc + max( wwint(4), wwint(2) )
558  END IF
559 !
560 ! *** error message ***
561 !
562 ! IF(WWINT(9) < WWINT(15) .OR. WWINT(10) > WWINT(16) .OR. &
563 ! WWINT(13) < WWINT(17) .OR. WWINT(14) > WWINT(18))THEN
564 ! WRITE (PRINTF,900) IXCGRD(1), IYCGRD(1), &
565 ! WWINT(9) ,WWINT(10) ,WWINT(13) ,WWINT(14), &
566 ! WWINT(15),WWINT(16) ,WWINT(17) ,WWINT(18)
567 !900 FORMAT ( ' ** Error : array bounds and maxima in subr RANGE4, ', &
568 ! ' point ', 2I5, &
569 ! /,' ISL,ISH : ',2I4, ' IDL,IDH : ',2I4, &
570 ! /,' SMI,SMA : ',2I4, ' DMI,DMA : ',2I4)
571 ! IF(ITEST > 50) WRITE (PRTEST, 901) MSC, MDC, IDDLOW, IDDTOP
572 !901 FORMAT (' MSC, MDC, IDDLOW, IDDTOP: ', 4I5)
573 ! ENDIF
574 !
575 ! test output
576 !
577 ! IF(TESTFL .AND. ITEST >= 60)THEN
578 ! WRITE(PRTEST,911) WWINT(4), WWINT(2), WWINT(8), WWINT(6)
579 !911 FORMAT (' RANGE4: IDM1 IDP1 ISM1 ISP1 :',4I5)
580 ! WRITE(PRTEST,916) WWINT(11), WWINT(12), IQUAD
581 !916 FORMAT (' RANGE4: ISCLW ISCHG IQUAD :',3I5)
582 ! WRITE (PRTEST,917) WWINT(9), WWINT(10), WWINT(13), WWINT(14)
583 !917 FORMAT (' RANGE4: ISLOW ISHGH IDLOW IDHGH:',4I5)
584 ! WRITE (PRTEST,919) WWINT(15), WWINT(16), WWINT(17), WWINT(18)
585 !919 FORMAT (' RANGE4: MS4MI MS4MA MD4MI MD4MA:',4I5)
586 ! WRITE(PRINTF,*)
587 ! END IF
588 !
589  RETURN
590  END SUBROUTINE range4
591 
592 !
593 !*******************************************************************
594 !
595  SUBROUTINE filnl3 (IDCMIN ,IDCMAX ,IMATRA ,IMATDA , &
596  MEMNL4 ,PLNL4S ,ISSTOP ,IGC )
597 !(This subroutine has not been tested yet)
598 !
599 !*******************************************************************
600 !
601  USE swcomm3
602  USE swcomm4
603  USE ocpcomm4
604 ! USE ALL_VARS, ONLY : MT,AC2
605  USE vars_wave, ONLY : mt,ac2
606 
607  IMPLICIT NONE
608 !
609 ! 2. Purpose
610 !
611 ! Fill the IMATRA/IMATDA arrays with the nonlinear wave-wave interaction
612 ! source term for a gridpoint ix,iy per sweep direction
613 !
614 !*******************************************************************
615 !
616  INTEGER IS,ID,ISSTOP,IDDUM,IGC
617 !
618  REAL IMATRA(MDC,MSC) , &
619  IMATDA(MDC,MSC) , &
620 ! AC2(MDC,MSC,0:MT) , &
621  plnl4s(mdc,msc,nptst) , &
622  memnl4(mdc,msc,mt)
623 !
624  INTEGER IDCMIN(MSC),IDCMAX(MSC)
625 !
626 ! SAVE IENT
627 ! DATA IENT/0/
628 ! IF (LTRACE) CALL STRACE (IENT,'FILNL3')
629 !
630  DO is=1, isstop
631  DO iddum = idcmin(is), idcmax(is)
632  id = mod( iddum - 1 + mdc , mdc ) + 1
633  IF(testfl) plnl4s(id,is,iptst) = memnl4(id,is,igc)
634  IF(memnl4(id,is,igc) > 0.)THEN
635  imatra(id,is) = imatra(id,is) + memnl4(id,is,igc)
636  ELSE
637  imatda(id,is) = imatda(id,is) - memnl4(id,is,igc) / &
638  max(1.e-18,ac2(id,is,igc))
639  END IF
640  END DO
641  END DO
642 !
643 ! IF(TESTFL .AND. ITEST >= 50)THEN
644 ! WRITE(PRINTF,9000) IDCMIN(1),IDCMAX(1),MSC,ISSTOP
645 !9000 FORMAT(' FILNL3: ID_MIN ID_MAX MSC ISTOP :',4I6)
646 ! IF(ITEST >= 100)THEN
647 ! DO IS=1, ISSTOP
648 ! DO IDDUM = IDCMIN(IS), IDCMAX(IS)
649 ! ID = MOD ( IDDUM - 1 + MDC , MDC ) + 1
650 ! WRITE(PRINTF,6001) IS,ID,MEMNL4(ID,IS,KCGRD(1))
651 !6001 FORMAT(' FILNL3: IS ID MEMNL() :',2I6,E12.4)
652 ! ENDDO
653 ! ENDDO
654 ! ENDIF
655 ! ENDIF
656 !
657  RETURN
658  END SUBROUTINE filnl3
659 
660 !
661 !*********************************************************************
662  SUBROUTINE swsnl8 (WWINT ,UE ,SA1 ,SA2 ,SPCSIG , &
663  SNLC1 ,DAL1 ,DAL2 ,DAL3 ,SFNL , &
664  DEP2 ,KMESPC ,MEMNL4 ,FACHFR )
665 
666 ! (This subroutine has not been tested yet)
667 !*********************************************************************
668 !
669  USE swcomm3
670  USE swcomm4
671  USE ocpcomm4
672  USE m_snl4
673 ! USE ALL_VARS, ONLY : MT,AC2
674  USE vars_wave, ONLY : mt,ac2
675 !
676  IMPLICIT NONE
677 !
678 ! 2. Purpose
679 !
680 ! Calculate non-linear interaction using the discrete interaction
681 ! approximation (Hasselmann and Hasselmann 1985; WAMDI group 1988)
682 ! for the full circle.
683 !
684 ! 3. Method
685 !
686 ! Discrete interaction approximation. To make interpolation simple,
687 ! the interactions are calculated in a "folded" space.
688 !
689 ! Frequencies -->
690 ! +---+---------------------+---------+- IDHGH
691 ! d | 3 : 2 : 2 |
692 ! i + - + - - - - - - - - - - + - - - - +- MDC
693 ! r | : : |
694 ! e | 3 : original spectrum : 1 |
695 ! c | : : |
696 ! t. + - + - - - - - - - - - - + - - - - +- 1
697 ! | 3 : 2 : 2 |
698 ! +---+---------------------+---------+- IDLOW
699 ! | | | ^ |
700 ! ISLOW 1 MSC | ISHGH
701 ! | |
702 ! ISCLW ISCHG
703 ! lowest discrete highest discrete
704 ! central bin central bin
705 !
706 ! 1 : Extra tail added beyond MSC
707 ! 2 : Spectrum copied outside ID range
708 ! 3 : Empty bins at low frequencies
709 !
710 ! ISLOW = 1 + ISM1
711 ! ISHGH = MSC + ISP1 - ISM1
712 ! ISCLW = 1
713 ! ISCHG = MSC - ISM1
714 ! IDLOW = 1 - MAX(IDM1,IDP1)
715 ! IDHGH = MDC + MAX(IDM1,IDP1)
716 !
717 ! Note: using this subroutine requires an additional array
718 ! with size MXC*MYC*MDC*MSC.
719 !
720  INTEGER WWINT(*)
721 !
722  REAL DAL1, DAL2, DAL3, FACHFR, KMESPC, SNLC1
723 ! REAL AC2(MDC,MSC,0:MT)
724  REAL DEP2(MT)
725  REAL MEMNL4(MDC,MSC,MT)
726  REAL SA1(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
727  REAL SA2(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
728  REAL SFNL(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
729  REAL SPCSIG(MSC)
730  REAL UE(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
731 !
732 !*******************************************************************
733 !
734  INTEGER IS ,ID ,ID0 ,I ,J , &
735  ISLOW ,ISHGH ,IDLOW ,IDHGH , &
736  IDPP ,IDMM ,ISPP ,ISMM , &
737  ISCLW ,ISCHG ,IDDUM ,IENT
738 !
739  REAL X ,X2 ,CONS ,FACTOR ,SNLCS1 ,SNLCS2 , &
740  SNLCS3 ,E00 ,EP1 ,EM1 ,EP2 ,EM2 , &
741  SA1A ,SA1B ,SA2A ,SA2B , &
742  JACOBI ,SIGPI
743 !
744 ! SAVE IENT
745 ! DATA IENT/0/
746 ! IF (LTRACE) CALL STRACE (IENT,'SWSNL8')
747 !
748  islow = wwint(9)
749  ishgh = wwint(10)
750  isclw = wwint(11)
751  ischg = wwint(12)
752  idlow = wwint(13)
753  idhgh = wwint(14)
754  idpp = wwint(21)
755  idmm = wwint(22)
756  ispp = wwint(23)
757  ismm = wwint(24)
758 !
759 ! *** Calculate prop. constant. ***
760 ! *** Calculate factor R(X) to calculate the NL wave-wave ***
761 ! *** interaction for shallow water ***
762 ! *** SNLC1 = 1/GRAV**4 ***
763 !
764  snlcs1 = pquad(3)
765  snlcs2 = pquad(4)
766  snlcs3 = pquad(5)
767 ! X = MAX ( 0.75 * DEP2(IGC) * KMESPC , 0.5 )
768  x = max( 0.75 * dep2(kcgrd(1)) * kmespc , 0.5 )
769  x2 = max( -1.e15, snlcs3*x)
770  cons = snlc1 * ( 1. + snlcs1/x * (1.-snlcs2*x) * exp(x2))
771  jacobi = 2. * pi_w
772 !
773 ! *** extend the area with action density at periodic boundaries ***
774 !
775  DO iddum = idlow, idhgh
776  id = mod( iddum - 1 + mdc , mdc ) + 1
777  DO is=1, msc
778 ! UE (IS,IDDUM) = AC2(ID,IS,IGC) * SPCSIG(IS) * JACOBI
779  ue(is,iddum) = ac2(id,is,kcgrd(1)) * spcsig(is) * jacobi
780  ENDDO
781  ENDDO
782 !
783  DO id = idlow, idhgh
784  DO is = msc+1, ishgh
785  ue(is,id) = ue(is-1,id) * fachfr
786  ENDDO
787  ENDDO
788 !
789 ! *** Calculate (unfolded) interactions ***
790 ! *** Energy at interacting bins ***
791 !
792  DO id = 1, mdc
793  DO is = isclw, ischg
794  e00 = ue(is ,id )
795  ep1 = ue(is+ispp,id+idpp)
796  em1 = ue(is+ismm,id-idmm)
797  ep2 = ue(is+ispp,id-idpp)
798  em2 = ue(is+ismm,id+idmm)
799 !
800 ! Contribution to interactions
801 !
802  factor = cons * af11(is) * pquad(2) * e00
803 !
804  sa1a = e00 * ( ep1*dal1 + em1*dal2 )
805  sa1b = sa1a - ep1*em1*dal3
806  sa2a = e00 * ( ep2*dal1 + em2*dal2 )
807  sa2b = sa2a - ep2*em2*dal3
808 !
809  sa1(is,id) = factor * sa1b
810  sa2(is,id) = factor * sa2b
811 !
812  ENDDO
813  ENDDO
814 !
815 ! *** Fold interactions to side angles -> domain in theta is ***
816 ! *** periodic ***
817 !
818  DO id = 1, idhgh - mdc
819  id0 = 1 - id
820  DO is = isclw, ischg
821  sa1(is,mdc+id) = sa1(is, id )
822  sa2(is,mdc+id) = sa2(is, id )
823  sa1(is, id0 ) = sa1(is,mdc+id0)
824  sa2(is, id0 ) = sa2(is,mdc+id0)
825  ENDDO
826  ENDDO
827 !
828 ! *** Put source term together (To save space I=IS and ***
829 ! *** J=MDC is used) ---- ***
830 !
831  DO i = 1, msc
832  sigpi = spcsig(i) * jacobi
833  DO j = 1, mdc
834  sfnl(i,j) = - 2. * ( sa1(i,j) + sa2(i,j) ) &
835  + ( sa1(i-ispp,j-idpp) + sa2(i-ispp,j+idpp) ) &
836  + ( sa1(i-ismm,j+idmm) + sa2(i-ismm,j-idmm) )
837 !
838 ! *** store value in auxiliary array and use values in ***
839 ! *** next four sweeps (see subroutine FILNL3) ***
840 !
841 ! MEMNL4(J,I,IGC) = SFNL(I,J) / SIGPI
842  memnl4(j,i,kcgrd(1)) = sfnl(i,j) / sigpi
843  ENDDO
844  ENDDO
845 !
846 ! *** value source term in every bin ***
847 !
848 ! IF(ITEST >= 150 .AND. TESTFL)THEN
849 ! DO I=1, MSC
850 ! DO J=1, MDC
851 ! WRITE(PRINTF,2006) I,J,MEMNL4(J,I,KCGRD(1)),SFNL(I,J),SPCSIG(I)
852 !2006 FORMAT (' I J MEMNL() SFNL() SPCSIG:',2I4,3E12.4)
853 ! ENDDO
854 ! ENDDO
855 ! END IF
856 !
857  RETURN
858 !
859  END SUBROUTINE swsnl8
860 
861 !
862 !********************************************************************
863 !
864 ! << Numerical Computations of the Nonlinear Energy Transfer
865 ! of Gravity Wave Spectra in Finite Water Depth >>
866 !
867 ! << developed by Noriaki Hashimoto >>
868 !
869 ! References: N. Hashimoto et al. (1998)
870 ! Komatsu and Masuda (2000) (in Japanese)
871 !
872 !
873 ! SWAN (Simulating WAves Nearshore); a third generation wave model
874 ! Copyright (C) 2004-2005 Delft University of Technology
875 !
876 ! This program is free software; you can redistribute it and/or
877 ! modify it under the terms of the GNU General Public License as
878 ! published by the Free Software Foundation; either version 2 of
879 ! the License, or (at your option) any later version.
880 !
881 ! This program is distributed in the hope that it will be useful,
882 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
883 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
884 ! GNU General Public License for more details.
885 !
886 ! A copy of the GNU General Public License is available at
887 ! http://www.gnu.org/copyleft/gpl.html#SEC3
888 ! or by writing to the Free Software Foundation, Inc.,
889 ! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
890 !
891 !
892  SUBROUTINE riam_slw(LMAX ,N ,N2 ,G , &
893  H ,DQ ,DQ2 ,DT , &
894  DT2 ,W ,P ,ACT , &
895  SNL ,MINT )
896 
897 ! (This subroutine has not been tested yet)
898 
899  USE m_snl4
900 
901 ! LMAX
902 ! N : number of directional bins
903 ! N2
904 ! G : gravitational acceleration
905 ! H : depth
906 ! DQ
907 ! DQ2
908 ! DT : size of the directional bins (Delta Theta)
909 ! DT2
910 ! W : discretised frequency array
911 ! P : density of water
912 ! ACT : action density
913 ! SNL : Quadruplet source term
914 ! MINT
915 
916 ! IW4 : counter for the 4th frequency
917 ! W4 : frequency of the fourth quadruplet component
918 ! AK4 : wavenumber of the fourth quadruplet component
919 ! DNA : coefficient in eq. 17 of Hashimoto (98)
920 ! = 2 w4 k4 / Cg(k4)
921 ! IW3L
922 ! II
923 ! JJ
924 ! DI
925 ! DJ
926 ! CGK4 : group velocity for the fourth quadruplet component
927 
928  REAL :: W(LMAX), ACT(LMAX,N), SNL(LMAX,N)
929 
930 ! Initialisation of the quadruplet source term
931 
932  snl(:,:) = 0.
933 !
934 ! =================
935  DO 130 iw4=1,lmax
936 ! =================
937 !
938  w4=w(iw4)
939 
940 ! WAVE converts nondimensional Sqr(w)d/g to nondimensional kd
941 
942  ak4=wave(w4**2*h/g)/h
943 
944 ! CGCMP computes group velocity
945 
946  CALL cgcmp(g,ak4,h,cgk4)
947 
948 ! Calculates the coefficient in equation (17) of Hashimoto (1998)
949 
950  dna=2.*w4*ak4/cgk4
951 
952  iw3l=max0(1,nint(iw4-alog(3.)/dq))
953 
954 ! ------------------------------------------------------------
955  CALL preset(lmax,n,iw3l,iw4,n2,g,h,w4,ak4,w,dq,dq2,dt,dt2,p)
956 ! ------------------------------------------------------------
957 !
958 ! ==============
959  DO 130 it4=1,n
960 ! ==============
961 !
962 ! ================
963  curriam => friam
964  DO ! (W1 - W3 - W4 - W2)
965 ! ================
966 !
967  k1=curriam%II(1)+iw4
968  k2=curriam%II(2)+iw4
969  k3=curriam%II(3)+iw4
970 !
971  IF(k1.LT.1 ) GO TO 120
972  IF(k2.GT.lmax) GO TO 120
973 !
974  gg=curriam%SSS*dna
975 !
976  m1p= curriam%JJ(1)+it4
977  m1n=-curriam%JJ(1)+it4
978  m2p= curriam%JJ(2)+it4
979  m2n=-curriam%JJ(2)+it4
980  m3p= curriam%JJ(3)+it4
981  m3n=-curriam%JJ(3)+it4
982 !
983  m1p=mod(m1p,n)
984  m1n=mod(m1n,n)
985  m2p=mod(m2p,n)
986  m2n=mod(m2n,n)
987  m3p=mod(m3p,n)
988  m3n=mod(m3n,n)
989 !
990  IF(m1p.LT.1) m1p=m1p+n
991  IF(m1n.LT.1) m1n=m1n+n
992  IF(m2p.LT.1) m2p=m2p+n
993  IF(m2n.LT.1) m2n=m2n+n
994  IF(m3p.LT.1) m3p=m3p+n
995  IF(m3n.LT.1) m3n=m3n+n
996 !
997  a4=act(iw4,it4)
998 !
999  IF(mint.EQ.1) THEN
1000 !
1001  k01=0
1002  k02=0
1003  k03=0
1004  m01=0
1005  m02=0
1006  m03=0
1007  IF(curriam%DI(1).LE.1.) k01= 1
1008  IF(curriam%DI(2).LE.1.) k02= 1
1009  IF(curriam%DI(3).LE.1.) k03= 1
1010  IF(curriam%DJ(1).LE.1.) m01=-1
1011  IF(curriam%DJ(2).LE.1.) m02=-1
1012  IF(curriam%DJ(3).LE.1.) m03=-1
1013 !
1014  k11=k1+k01
1015  k21=k2+k02
1016  k31=k3+k03
1017  IF(k11.GT.lmax) k11=lmax
1018  IF(k21.GT.lmax) k21=lmax
1019  IF(k31.GT.lmax) k31=lmax
1020 !
1021  k111=k1+k01-1
1022  k211=k2+k02-1
1023  k311=k3+k03-1
1024  IF(k111.LT.1) k111=1
1025  IF(k211.LT.1) k211=1
1026  IF(k311.LT.1) k311=1
1027 !
1028  m1p1=m1p+m01
1029  m2p1=m2p+m02
1030  m3p1=m3p+m03
1031  IF(m1p1.LT.1) m1p1=n
1032  IF(m2p1.LT.1) m2p1=n
1033  IF(m3p1.LT.1) m3p1=n
1034 !
1035  m1n1=m1n-m01
1036  m2n1=m2n-m02
1037  m3n1=m3n-m03
1038  IF(m1n1.GT.n) m1n1=1
1039  IF(m2n1.GT.n) m2n1=1
1040  IF(m3n1.GT.n) m3n1=1
1041 !
1042  m1p11=m1p+m01+1
1043  m2p11=m2p+m02+1
1044  m3p11=m3p+m03+1
1045  IF(m1p11.GT.n) m1p11=1
1046  IF(m2p11.GT.n) m2p11=1
1047  IF(m3p11.GT.n) m3p11=1
1048 !
1049  m1n11=m1n-m01-1
1050  m2n11=m2n-m02-1
1051  m3n11=m3n-m03-1
1052  IF(m1n11.LT.1) m1n11=n
1053  IF(m2n11.LT.1) m2n11=n
1054  IF(m3n11.LT.1) m3n11=n
1055 !
1056  di1=curriam%DI(1)
1057  di2=curriam%DI(2)
1058  di3=curriam%DI(3)
1059  dj1=curriam%DJ(1)
1060  dj2=curriam%DJ(2)
1061  dj3=curriam%DJ(3)
1062 !
1063  a1p=(di1*(dj1*act(k11, m1p1)+act(k11, m1p11)) &
1064  + dj1*act(k111,m1p1)+act(k111,m1p11)) &
1065  /((1.+di1)*(1.+dj1))
1066 !
1067  a1n=(di1*(dj1*act(k11, m1n1)+act(k11, m1n11)) &
1068  + dj1*act(k111,m1n1)+act(k111,m1n11)) &
1069  /((1.+di1)*(1.+dj1))
1070 !
1071  a2p=(di2*(dj2*act(k21, m2p1)+act(k21, m2p11)) &
1072  + dj2*act(k211,m2p1)+act(k211,m2p11)) &
1073  /((1.+di2)*(1.+dj2))
1074 !
1075  a2n=(di2*(dj2*act(k21, m2n1)+act(k21, m2n11)) &
1076  + dj2*act(k211,m2n1)+act(k211,m2n11)) &
1077  /((1.+di2)*(1.+dj2))
1078 !
1079  a3p=(di3*(dj3*act(k31, m3p1)+act(k31, m3p11)) &
1080  + dj3*act(k311,m3p1)+act(k311,m3p11)) &
1081  /((1.+di3)*(1.+dj3))
1082 !
1083  a3n=(di3*(dj3*act(k31, m3n1)+act(k31, m3n11)) &
1084  + dj3*act(k311,m3n1)+act(k311,m3n11)) &
1085  /((1.+di3)*(1.+dj3))
1086 !
1087  ELSE
1088 !
1089  IF(k1.LT.1) k1=1
1090  IF(k2.LT.1) k2=1
1091  IF(k3.LT.1) k3=1
1092 !
1093  a1p=act(k1,m1p)
1094  a1n=act(k1,m1n)
1095  a2p=act(k2,m2p)
1096  a2n=act(k2,m2n)
1097  a3p=act(k3,m3p)
1098  a3n=act(k3,m3n)
1099 !
1100  ENDIF
1101 !
1102  w1p2p=a1p+a2p
1103  w1n2n=a1n+a2n
1104  s1p2p=a1p*a2p
1105  s1n2n=a1n*a2n
1106  w3p4 =a3p+a4
1107  w3n4 =a3n+a4
1108  s3p4 =a3p*a4
1109  s3n4 =a3n*a4
1110 !
1111  xp=s1p2p*w3p4-s3p4*w1p2p
1112  xn=s1n2n*w3n4-s3n4*w1n2n
1113 !
1114  snl( k1,m1p)=snl( k1,m1p)-xp*gg
1115  snl( k1,m1n)=snl( k1,m1n)-xn*gg
1116  snl( k2,m2p)=snl( k2,m2p)-xp*gg
1117  snl( k2,m2n)=snl( k2,m2n)-xn*gg
1118  snl( k3,m3p)=snl( k3,m3p)+xp*gg
1119  snl( k3,m3n)=snl( k3,m3n)+xn*gg
1120  snl(iw4,it4)=snl(iw4,it4)+(xp+xn)*gg
1121 !
1122 ! ============
1123  120 CONTINUE
1124  IF (.NOT.ASSOCIATED(curriam%NEXTRIAM)) EXIT
1125  curriam => curriam%NEXTRIAM
1126  END DO
1127 ! ============
1128 !
1129 ! ============
1130  130 CONTINUE
1131 ! ============
1132 !
1133  RETURN
1134  END SUBROUTINE riam_slw
1135 
1136 !
1137 !*******************************************************************
1138 !
1139  SUBROUTINE swsnl4 (WWINT ,WWAWG ,SPCSIG ,SNLC1 , &
1140  DAL1 ,DAL2 ,DAL3 ,DEP2 , &
1141  KMESPC ,MEMNL4 ,FACHFR ,IDIA , &
1142  ITER )
1143 
1144 ! (This subroutine has not been tested yet)
1145 !
1146 !*******************************************************************
1147 !
1148  USE swcomm3
1149  USE swcomm4
1150  USE ocpcomm4
1151  USE m_snl4
1152 ! USE ALL_VARS, ONLY : MT,AC2
1153  USE vars_wave, ONLY : mt,ac2
1154 !
1155 ! --|-----------------------------------------------------------|--
1156 ! | Delft University of Technology |
1157 ! | Faculty of Civil Engineering |
1158 ! | Fluid Mechanics Section |
1159 ! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
1160 ! | |
1161 ! | Programmers: H.L. Tolman, R.C. Ris |
1162 ! --|-----------------------------------------------------------|--
1163 !
1164 !
1165 ! SWAN (Simulating WAves Nearshore); a third generation wave model
1166 ! Copyright (C) 2004-2005 Delft University of Technology
1167 !
1168 ! This program is free software; you can redistribute it and/or
1169 ! modify it under the terms of the GNU General Public License as
1170 ! published by the Free Software Foundation; either version 2 of
1171 ! the License, or (at your option) any later version.
1172 !
1173 ! This program is distributed in the hope that it will be useful,
1174 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
1175 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1176 ! GNU General Public License for more details.
1177 !
1178 ! A copy of the GNU General Public License is available at
1179 ! http://www.gnu.org/copyleft/gpl.html#SEC3
1180 ! or by writing to the Free Software Foundation, Inc.,
1181 ! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
1182 !
1183 !
1184 ! 0. Authors
1185 !
1186 ! 40.17: IJsbrand Haagsma
1187 ! 40.41: Marcel Zijlema
1188 !
1189 ! 1. Updates
1190 !
1191 ! 40.17, Dec. 01: New Subroutine based on SWSNL3
1192 ! 40.41, Oct. 04: common blocks replaced by modules, include files removed
1193 !
1194 ! 2. Purpose
1195 !
1196 ! Calculate non-linear interaction using the discrete interaction
1197 ! approximation (Hasselmann and Hasselmann 1985; WAMDI group 1988)
1198 ! for the full circle (option if a current is present). Note: using
1199 ! this subroutine requires an additional array with size
1200 ! (MXC*MYC*MDC*MSC). This requires more internal memory but can
1201 ! speed up the computations sigificantly if a current is present.
1202 !
1203 ! 3. Method
1204 !
1205 ! Discrete interaction approximation. To make interpolation simple,
1206 ! the interactions are calculated in a "folded" space.
1207 !
1208 ! Frequencies -->
1209 ! +---+---------------------+---------+- IDHGH
1210 ! d | 3 : 2 : 2 |
1211 ! i + - + - - - - - - - - - - + - - - - +- MDC
1212 ! r | : : |
1213 ! e | 3 : original spectrum : 1 |
1214 ! c | : : |
1215 ! t. + - + - - - - - - - - - - + - - - - +- 1
1216 ! | 3 : 2 : 2 |
1217 ! +---+---------------------+---------+- IDLOW
1218 ! | | | ^ |
1219 ! ISLOW 1 MSC | ISHGH
1220 ! | |
1221 ! ISCLW ISCHG
1222 ! lowest discrete highest discrete
1223 ! central bin central bin
1224 !
1225 ! 1 : Extra tail added beyond MSC
1226 ! 2 : Spectrum copied outside ID range
1227 ! 3 : Empty bins at low frequencies
1228 !
1229 ! ISLOW = 1 + ISM1
1230 ! ISHGH = MSC + ISP1 - ISM1
1231 ! ISCLW = 1
1232 ! ISCHG = MSC - ISM1
1233 ! IDLOW = 1 - MAX(IDM1,IDP1)
1234 ! IDHGH = MDC + MAX(IDM1,IDP1)
1235 !
1236 ! Relative offsets of interpolation points around central bin
1237 ! "#" and corresponding numbers of AWGn :
1238 !
1239 ! ISM1 ISM
1240 ! 5 7 T |
1241 ! IDM1 +------+ H +
1242 ! | | E | ISP ISP1
1243 ! | \ | T | 3 1
1244 ! IDM +------+ A + +---------+ IDP1
1245 ! 6 \8 | | |
1246 ! | | / |
1247 ! \ + +---------+ IDP
1248 ! | /4 2
1249 ! \ | /
1250 ! -+-----+------+-------#--------+---------+----------+
1251 ! | FREQ.
1252 !
1253 !
1254 ! 4. Argument variables
1255 !
1256 ! MCGRD : number of wet grid points of the computational grid
1257 ! MDC : grid points in theta-direction of computational grid
1258 ! MDC4MA: highest array counter in directional space (Snl4)
1259 ! MDC4MI: lowest array counter in directional space (Snl4)
1260 ! MSC : grid points in sigma-direction of computational grid
1261 ! MSC4MA: highest array counter in frequency space (Snl4)
1262 ! MSC4MI: lowest array counter in frequency space (Snl4)
1263 ! WWINT : counters for quadruplet interactions
1264 !
1265  INTEGER WWINT(*)
1266  INTEGER IDIA
1267 !
1268 ! AC2 : action density
1269 ! AF11 : scaling frequency
1270 ! DAL1 : coefficient for the quadruplet interactions
1271 ! DAL2 : coefficient for the quadruplet interactions
1272 ! DAL3 : coefficient for the quadruplet interactions
1273 ! DEP2 : depth
1274 ! FACHFR
1275 ! KMESPC: mean average wavenumber over full spectrum
1276 ! MEMNL4
1277 ! PI : circular constant
1278 ! SA1 : interaction contribution of first quadruplet (unfolded space)
1279 ! SA2 : interaction contribution of second quadruplet (unfolded space)
1280 ! SFNL
1281 ! SNLC1
1282 ! SPCSIG: relative frequencies in computational domain in sigma-space
1283 ! UE : "unfolded" spectrum
1284 ! WWAWG : weight coefficients for the quadruplet interactions
1285 !
1286  REAL DAL1, DAL2, DAL3, FACHFR, KMESPC, SNLC1
1287 ! REAL AC2(MDC,MSC,0:MT)
1288  REAL DEP2(MT)
1289  REAL MEMNL4(MDC,MSC,MT)
1290  REAL SPCSIG(MSC)
1291  REAL WWAWG(*)
1292 !
1293 ! 6. Local variables
1294 !
1295  REAL, ALLOCATABLE :: SA1(:,:), SA2(:,:), SFNL(:,:), UE(:,:)
1296 !
1297 ! 7. Common blocks used
1298 !
1299 !
1300 ! 8. Subroutines used
1301 !
1302 ! ---
1303 !
1304 ! 9. Subroutines calling
1305 !
1306 ! SOURCE (in SWANCOM1)
1307 !
1308 ! 12. Structure
1309 !
1310 ! -------------------------------------------
1311 ! Initialisations.
1312 ! Calculate proportionality constant.
1313 ! Prepare auxiliary spectrum.
1314 ! Calculate (unfolded) interactions :
1315 ! -----------------------------------------
1316 ! Energy at interacting bins
1317 ! Contribution to interactions
1318 ! Fold interactions to side angles
1319 ! -----------------------------------------
1320 ! Put source term together
1321 ! -------------------------------------------
1322 !
1323 ! 13. Source text
1324 !
1325 !*******************************************************************
1326 !
1327  INTEGER IS ,ID ,ID0 ,I ,J , &
1328  ISHGH ,IDLOW ,IDHGH ,ISP ,ISP1 , &
1329  IDP ,IDP1 ,ISM ,ISM1 ,IDM ,IDM1 , &
1330  ISCLW ,ISCHG
1331 !
1332  REAL X ,X2 ,CONS ,FACTOR ,SNLCS2 , &
1333  SNLCS3 ,E00 ,EP1 ,EM1 ,EP2 ,EM2 , &
1334  SA1A ,SA1B ,SA2A ,SA2B , &
1335  AWG1 ,AWG2 ,AWG3 ,AWG4 ,AWG5 ,AWG6 , &
1336  AWG7 ,AWG8 , &
1337  JACOBI ,SIGPI
1338 !
1339 ! SAVE IENT
1340 ! DATA IENT/0/
1341 ! IF (LTRACE) CALL STRACE (IENT,'SWSNL4')
1342 
1343  ALLOCATE(sa1(msc4mi:msc4ma,mdc4mi:mdc4ma))
1344  ALLOCATE(sa2(msc4mi:msc4ma,mdc4mi:mdc4ma))
1345  ALLOCATE(sfnl(msc4mi:msc4ma,mdc4mi:mdc4ma))
1346  ALLOCATE(ue(msc4mi:msc4ma,mdc4mi:mdc4ma))
1347 !
1348  idp = wwint(1)
1349  idp1 = wwint(2)
1350  idm = wwint(3)
1351  idm1 = wwint(4)
1352  isp = wwint(5)
1353  isp1 = wwint(6)
1354  ism = wwint(7)
1355  ism1 = wwint(8)
1356  islow = wwint(9)
1357  ishgh = wwint(10)
1358  isclw = wwint(11)
1359  ischg = wwint(12)
1360  idlow = wwint(13)
1361  idhgh = wwint(14)
1362 !
1363  awg1 = wwawg(1)
1364  awg2 = wwawg(2)
1365  awg3 = wwawg(3)
1366  awg4 = wwawg(4)
1367  awg5 = wwawg(5)
1368  awg6 = wwawg(6)
1369  awg7 = wwawg(7)
1370  awg8 = wwawg(8)
1371 !
1372 ! *** Initialize auxiliary arrays per gridpoint ***
1373 !
1374  DO id = mdc4mi, mdc4ma
1375  DO is = msc4mi, msc4ma
1376  ue(is,id) = 0.
1377  sa1(is,id) = 0.
1378  sa2(is,id) = 0.
1379  sfnl(is,id) = 0.
1380  ENDDO
1381  ENDDO
1382 !
1383 ! *** Calculate prop. constant. ***
1384 ! *** Calculate factor R(X) to calculate the NL wave-wave ***
1385 ! *** interaction for shallow water ***
1386 ! *** SNLC1 = 1/GRAV**4 ***
1387 !
1388  snlcs1 = pquad(3)
1389  snlcs2 = pquad(4)
1390  snlcs3 = pquad(5)
1391 ! X = MAX ( 0.75 * DEP2(IGC) * KMESPC , 0.5 )
1392  x = max( 0.75 * dep2(kcgrd(1)) * kmespc , 0.5 )
1393  x2 = max( -1.e15, snlcs3*x)
1394  cons = snlc1 * ( 1. + snlcs1/x * (1.-snlcs2*x) * exp(x2))
1395  jacobi = 2. * pi_w
1396 !
1397 ! *** extend the area with action density at periodic boundaries ***
1398 !
1399  DO iddum = idlow, idhgh
1400  id = mod( iddum - 1 + mdc , mdc ) + 1
1401  DO is=1, msc
1402 ! UE (IS,IDDUM) = AC2(ID,IS,IGC) * SPCSIG(IS) * JACOBI
1403  ue(is,iddum) = ac2(id,is,kcgrd(1)) * spcsig(is) * jacobi
1404  ENDDO
1405  ENDDO
1406 !
1407  DO is = msc+1, ishgh
1408  DO id = idlow, idhgh
1409  ue(is,id) = ue(is-1,id) * fachfr
1410  ENDDO
1411  ENDDO
1412 !
1413 ! *** Calculate (unfolded) interactions ***
1414 ! *** Energy at interacting bins ***
1415 !
1416  DO is = isclw, ischg
1417  DO id = 1, mdc
1418  e00 = ue(is ,id )
1419  ep1 = awg1 * ue(is+isp1,id+idp1) + &
1420  awg2 * ue(is+isp1,id+idp ) + &
1421  awg3 * ue(is+isp ,id+idp1) + &
1422  awg4 * ue(is+isp ,id+idp )
1423  em1 = awg5 * ue(is+ism1,id-idm1) + &
1424  awg6 * ue(is+ism1,id-idm ) + &
1425  awg7 * ue(is+ism ,id-idm1) + &
1426  awg8 * ue(is+ism ,id-idm )
1427  ep2 = awg1 * ue(is+isp1,id-idp1) + &
1428  awg2 * ue(is+isp1,id-idp ) + &
1429  awg3 * ue(is+isp ,id-idp1) + &
1430  awg4 * ue(is+isp ,id-idp )
1431  em2 = awg5 * ue(is+ism1,id+idm1) + &
1432  awg6 * ue(is+ism1,id+idm ) + &
1433  awg7 * ue(is+ism ,id+idm1) + &
1434  awg8 * ue(is+ism ,id+idm )
1435 !
1436 ! Contribution to interactions
1437 !
1438  factor = cons * af11(is) * e00
1439 !
1440  sa1a = e00 * ( ep1*dal1 + em1*dal2 ) * cnl4_1(idia)
1441  sa1b = sa1a - ep1*em1*dal3 * cnl4_2(idia)
1442  sa2a = e00 * ( ep2*dal1 + em2*dal2 ) * cnl4_1(idia)
1443  sa2b = sa2a - ep2*em2*dal3 * cnl4_2(idia)
1444 !
1445 
1446  sa1(is,id) = factor * sa1b
1447  sa2(is,id) = factor * sa2b
1448 !
1449 ! IF(ITEST >= 100 .AND. TESTFL)THEN
1450 ! WRITE(PRINTF,9002) E00,EP1,EM1,EP2,EM2
1451 !9002 FORMAT (' E00 EP1 EM1 EP2 EM2 :',5E11.4)
1452 ! WRITE(PRINTF,9003) SA1A,SA1B,SA2A,SA2B
1453 !9003 FORMAT (' SA1A SA1B SA2A SA2B :',4E11.4)
1454 ! WRITE(PRINTF,9004) IS,ID,SA1(IS,ID),SA2(IS,ID)
1455 !9004 FORMAT (' IS ID SA1() SA2() :',2I4,2E12.4)
1456 ! WRITE(PRINTF,9005) FACTOR,JACOBI
1457 !9005 FORMAT (' FACTOR JACOBI : ',2E12.4)
1458 ! END IF
1459 !
1460  ENDDO
1461  ENDDO
1462 !
1463 ! *** Fold interactions to side angles -> domain in theta is ***
1464 ! *** periodic ***
1465 !
1466  DO id = 1, idhgh - mdc
1467  id0 = 1 - id
1468  DO is = isclw, ischg
1469  sa1(is,mdc+id) = sa1(is, id )
1470  sa2(is,mdc+id) = sa2(is, id )
1471  sa1(is, id0 ) = sa1(is,mdc+id0)
1472  sa2(is, id0 ) = sa2(is,mdc+id0)
1473  ENDDO
1474  ENDDO
1475 !
1476 ! *** Put source term together (To save space I=IS and ***
1477 ! *** J=MDC is used) ***
1478 !
1479  fac = 1.
1480 
1481  DO i = 1, msc
1482  sigpi = spcsig(i) * jacobi
1483  DO j = 1, mdc
1484  sfnl(i,j) = - 2. * ( sa1(i,j) + sa2(i,j) ) &
1485  + awg1 * ( sa1(i-isp1,j-idp1) + sa2(i-isp1,j+idp1) ) &
1486  + awg2 * ( sa1(i-isp1,j-idp ) + sa2(i-isp1,j+idp ) ) &
1487  + awg3 * ( sa1(i-isp ,j-idp1) + sa2(i-isp ,j+idp1) ) &
1488  + awg4 * ( sa1(i-isp ,j-idp ) + sa2(i-isp ,j+idp ) ) &
1489  + awg5 * ( sa1(i-ism1,j+idm1) + sa2(i-ism1,j-idm1) ) &
1490  + awg6 * ( sa1(i-ism1,j+idm ) + sa2(i-ism1,j-idm ) ) &
1491  + awg7 * ( sa1(i-ism ,j+idm1) + sa2(i-ism ,j-idm1) ) &
1492  + awg8 * ( sa1(i-ism ,j+idm ) + sa2(i-ism ,j-idm ) )
1493 !
1494 ! *** store value in auxiliary array and use values in ***
1495 ! *** next four sweeps (see subroutine FILNL3) ***
1496 !
1497  IF(idia == 1)THEN
1498 ! MEMNL4(J,I,IGC) = FAC * SFNL(I,J) / SIGPI
1499  memnl4(j,i,kcgrd(1)) = fac * sfnl(i,j) / sigpi
1500  ELSE
1501 ! MEMNL4(J,I,IGC) = MEMNL4(J,I,IGC) + FAC * SFNL(I,J) / SIGPI
1502  memnl4(j,i,kcgrd(1)) = memnl4(j,i,kcgrd(1)) + fac * sfnl(i,j) / sigpi
1503  END IF
1504  ENDDO
1505  ENDDO
1506 !
1507 ! *** test output ***
1508 !
1509 ! IF(ITEST >= 50 .AND. TESTFL)THEN
1510 ! WRITE(PRINTF,*)
1511 ! WRITE(PRINTF,*) ' SWSNL4 subroutine '
1512 ! WRITE(PRINTF,9011) IDP, IDP1, IDM, IDM1
1513 !9011 FORMAT (' IDP IDP1 IDM IDM1 :',4I5)
1514 ! WRITE (PRINTF,9013) ISP, ISP1, ISM, ISM1
1515 !9013 FORMAT (' ISP ISP1 ISM ISM1 :',4I5)
1516 ! WRITE (PRINTF,9015) ISLOW, ISHGH, IDLOW, IDHGH
1517 !9015 FORMAT (' ISLOW ISHG IDLOW IDHG :',4I5)
1518 ! WRITE(PRINTF,9016) ISCLW, ISCHG, JACOBI
1519 !9016 FORMAT (' ICLW ICHG JACOBI :',2I5,E12.4)
1520 ! WRITE (PRINTF,9017) AWG1, AWG2, AWG3, AWG4
1521 !9017 FORMAT (' AWG1 AWG2 AWG3 AWG4 :',4E12.4)
1522 ! WRITE (PRINTF,9018) AWG5, AWG6, AWG7, AWG8
1523 !9018 FORMAT (' AWG5 AWG6 AWG7 AWG8 :',4E12.4)
1524 ! WRITE (PRINTF,9019) MSC4MI, MSC4MA, MDC4MI, MDC4MA
1525 !9019 FORMAT (' S4MI S4MA D4MI D4MA :',4I6)
1526 ! WRITE(PRINTF,9020) SNLC1,X,X2,CONS
1527 !9020 FORMAT (' SNLC1 X X2 CONS :',4E12.4)
1528 ! WRITE(PRINTF,9021) DEP2(KCGRD(1)),KMESPC,FACHFR,PI
1529 !9021 FORMAT (' DEPTH KMESPC FACHFR PI:',4E12.4)
1530 ! WRITE(PRINTF,*)
1531 !
1532 ! *** value source term in every bin ***
1533 !
1534 ! IF(ITEST >= 150)THEN
1535 ! DO I=1, MSC
1536 ! DO J=1, MDC
1537 ! WRITE(PRINTF,2006) I,J,MEMNL4(J,I,KCGRD(1)),SFNL(I,J), &
1538 ! SPCSIG(I)
1539 !2006 FORMAT (' I J MEMNL() SFNL() SPCSIG:',2I4,3E12.4)
1540 ! ENDDO
1541 ! ENDDO
1542 ! END IF
1543 ! END IF
1544 !
1545  DEALLOCATE (sa1, sa2, sfnl, ue)
1546 
1547  RETURN
1548 !
1549  END SUBROUTINE swsnl4
1550 
1551 !
1552 !************************************************************
1553 !
1554 
1555  SUBROUTINE swsnl3 (WWINT ,WWAWG ,UE ,SA1 ,SA2 , &
1556  SPCSIG ,SNLC1 ,DAL1 ,DAL2 ,DAL3 , &
1557  SFNL ,DEP2 ,KMESPC ,MEMNL4 ,FACHFR )
1558 ! (This subroutine has not been tested yet)
1559 !
1560 !*******************************************************************
1561 !
1562  USE swcomm3
1563  USE swcomm4
1564  USE ocpcomm4
1565  USE m_snl4
1566 ! USE ALL_VARS, ONLY : MT,AC2
1567  USE vars_wave, ONLY : mt,ac2
1568 !
1569 ! --|-----------------------------------------------------------|--
1570 ! | Delft University of Technology |
1571 ! | Faculty of Civil Engineering |
1572 ! | Fluid Mechanics Section |
1573 ! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
1574 ! | |
1575 ! | Programmers: H.L. Tolman, R.C. Ris |
1576 ! --|-----------------------------------------------------------|--
1577 !
1578 !
1579 ! SWAN (Simulating WAves Nearshore); a third generation wave model
1580 ! Copyright (C) 2004-2005 Delft University of Technology
1581 !
1582 ! This program is free software; you can redistribute it and/or
1583 ! modify it under the terms of the GNU General Public License as
1584 ! published by the Free Software Foundation; either version 2 of
1585 ! the License, or (at your option) any later version.
1586 !
1587 ! This program is distributed in the hope that it will be useful,
1588 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
1589 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1590 ! GNU General Public License for more details.
1591 !
1592 ! A copy of the GNU General Public License is available at
1593 ! http://www.gnu.org/copyleft/gpl.html#SEC3
1594 ! or by writing to the Free Software Foundation, Inc.,
1595 ! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
1596 !
1597 !
1598 ! 0. Authors
1599 !
1600 ! 30.72: IJsbrand Haagsma
1601 ! 40.17: IJsbrand Haagsma
1602 ! 40.41: Marcel Zijlema
1603 !
1604 ! 1. Updates
1605 !
1606 ! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN
1607 ! 40.17, Dec. 01: Implemented Multiple DIA
1608 ! 40.41, Oct. 04: common blocks replaced by modules, include files removed
1609 !
1610 ! 2. Purpose
1611 !
1612 ! Calculate non-linear interaction using the discrete interaction
1613 ! approximation (Hasselmann and Hasselmann 1985; WAMDI group 1988)
1614 ! for the full circle (option if a current is present). Note: using
1615 ! this subroutine requires an additional array with size
1616 ! (MXC*MYC*MDC*MSC). This requires more internal memory but can
1617 ! speed up the computations sigificantly if a current is present.
1618 !
1619 ! 3. Method
1620 !
1621 ! Discrete interaction approximation. To make interpolation simple,
1622 ! the interactions are calculated in a "folded" space.
1623 !
1624 ! Frequencies -->
1625 ! +---+---------------------+---------+- IDHGH
1626 ! d | 3 : 2 : 2 |
1627 ! i + - + - - - - - - - - - - + - - - - +- MDC
1628 ! r | : : |
1629 ! e | 3 : original spectrum : 1 |
1630 ! c | : : |
1631 ! t. + - + - - - - - - - - - - + - - - - +- 1
1632 ! | 3 : 2 : 2 |
1633 ! +---+---------------------+---------+- IDLOW
1634 ! | | | ^ |
1635 ! ISLOW 1 MSC | ISHGH
1636 ! | |
1637 ! ISCLW ISCHG
1638 ! lowest discrete highest discrete
1639 ! central bin central bin
1640 !
1641 ! 1 : Extra tail added beyond MSC
1642 ! 2 : Spectrum copied outside ID range
1643 ! 3 : Empty bins at low frequencies
1644 !
1645 ! ISLOW = 1 + ISM1
1646 ! ISHGH = MSC + ISP1 - ISM1
1647 ! ISCLW = 1
1648 ! ISCHG = MSC - ISM1
1649 ! IDLOW = 1 - MAX(IDM1,IDP1)
1650 ! IDHGH = MDC + MAX(IDM1,IDP1)
1651 !
1652 ! Relative offsets of interpolation points around central bin
1653 ! "#" and corresponding numbers of AWGn :
1654 !
1655 ! ISM1 ISM
1656 ! 5 7 T |
1657 ! IDM1 +------+ H +
1658 ! | | E | ISP ISP1
1659 ! | \ | T | 3 1
1660 ! IDM +------+ A + +---------+ IDP1
1661 ! 6 \8 | | |
1662 ! | | / |
1663 ! \ + +---------+ IDP
1664 ! | /4 2
1665 ! \ | /
1666 ! -+-----+------+-------#--------+---------+----------+
1667 ! | FREQ.
1668 !
1669 !
1670 ! 4. Argument variables
1671 !
1672 ! MCGRD : number of wet grid points of the computational grid
1673 ! MDC : grid points in theta-direction of computational grid
1674 ! MDC4MA: highest array counter in directional space (Snl4)
1675 ! MDC4MI: lowest array counter in directional space (Snl4)
1676 ! MSC : grid points in sigma-direction of computational grid
1677 ! MSC4MA: highest array counter in frequency space (Snl4)
1678 ! MSC4MI: lowest array counter in frequency space (Snl4)
1679 ! WWINT : counters for quadruplet interactions
1680 !
1681  INTEGER WWINT(*)
1682 !
1683 ! AC2 : action density
1684 ! AF11 : scaling frequency
1685 ! DAL1 : coefficient for the quadruplet interactions
1686 ! DAL2 : coefficient for the quadruplet interactions
1687 ! DAL3 : coefficient for the quadruplet interactions
1688 ! DEP2 : depth
1689 ! FACHFR
1690 ! KMESPC: mean average wavenumber over full spectrum
1691 ! MEMNL4
1692 ! PI : circular constant
1693 ! SA1 : interaction contribution of first quadruplet (unfolded space)
1694 ! SA2 : interaction contribution of second quadruplet (unfolded space)
1695 ! SFNL
1696 ! SNLC1
1697 ! SPCSIG: relative frequencies in computational domain in sigma-space
1698 ! UE : "unfolded" spectrum
1699 ! WWAWG : weight coefficients for the quadruplet interactions
1700 !
1701  REAL DAL1, DAL2, DAL3, FACHFR, KMESPC, SNLC1
1702 ! REAL AC2(MDC,MSC,0:MT)
1703  REAL DEP2(MT)
1704  REAL MEMNL4(MDC,MSC,MT)
1705  REAL SA1(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
1706  REAL SA2(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
1707  REAL SFNL(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
1708  REAL SPCSIG(MSC)
1709  REAL UE(MSC4MI:MSC4MA,MDC4MI:MDC4MA)
1710  REAL WWAWG(*)
1711 !
1712 ! 7. Common blocks used
1713 !
1714 !
1715 ! 8. Subroutines used
1716 !
1717 ! ---
1718 !
1719 ! 9. Subroutines calling
1720 !
1721 ! SOURCE (in SWANCOM1)
1722 !
1723 ! 12. Structure
1724 !
1725 ! -------------------------------------------
1726 ! Initialisations.
1727 ! Calculate proportionality constant.
1728 ! Prepare auxiliary spectrum.
1729 ! Calculate (unfolded) interactions :
1730 ! -----------------------------------------
1731 ! Energy at interacting bins
1732 ! Contribution to interactions
1733 ! Fold interactions to side angles
1734 ! -----------------------------------------
1735 ! Put source term together
1736 ! -------------------------------------------
1737 !
1738 ! 13. Source text
1739 !
1740 !*******************************************************************
1741 !
1742  INTEGER IS ,ID ,ID0 ,I ,J , &
1743  ISHGH ,IDLOW ,IDHGH ,ISP ,ISP1 , &
1744  IDP ,IDP1 ,ISM ,ISM1 ,IDM ,IDM1 , &
1745  ISCLW ,ISCHG
1746 !
1747  REAL X ,X2 ,CONS ,FACTOR ,SNLCS2 , &
1748  SNLCS3 ,E00 ,EP1 ,EM1 ,EP2 ,EM2 , &
1749  SA1A ,SA1B ,SA2A ,SA2B , &
1750  AWG1 ,AWG2 ,AWG3 ,AWG4 ,AWG5 ,AWG6 , &
1751  AWG7 ,AWG8 , &
1752  JACOBI ,SIGPI
1753 !
1754 ! SAVE IENT
1755 ! DATA IENT/0/
1756 ! IF (LTRACE) CALL STRACE (IENT,'SWSNL3')
1757 !
1758  idp = wwint(1)
1759  idp1 = wwint(2)
1760  idm = wwint(3)
1761  idm1 = wwint(4)
1762  isp = wwint(5)
1763  isp1 = wwint(6)
1764  ism = wwint(7)
1765  ism1 = wwint(8)
1766  islow = wwint(9)
1767  ishgh = wwint(10)
1768  isclw = wwint(11)
1769  ischg = wwint(12)
1770  idlow = wwint(13)
1771  idhgh = wwint(14)
1772 !
1773  awg1 = wwawg(1)
1774  awg2 = wwawg(2)
1775  awg3 = wwawg(3)
1776  awg4 = wwawg(4)
1777  awg5 = wwawg(5)
1778  awg6 = wwawg(6)
1779  awg7 = wwawg(7)
1780  awg8 = wwawg(8)
1781 !
1782 ! *** Initialize auxiliary arrays per gridpoint ***
1783 !
1784  DO id = mdc4mi, mdc4ma
1785  DO is = msc4mi, msc4ma
1786  ue(is,id) = 0.
1787  sa1(is,id) = 0.
1788  sa2(is,id) = 0.
1789  sfnl(is,id) = 0.
1790  ENDDO
1791  ENDDO
1792 !
1793 ! *** Calculate prop. constant. ***
1794 ! *** Calculate factor R(X) to calculate the NL wave-wave ***
1795 ! *** interaction for shallow water ***
1796 ! *** SNLC1 = 1/GRAV**4 ***
1797 !
1798  snlcs1 = pquad(3)
1799  snlcs2 = pquad(4)
1800  snlcs3 = pquad(5)
1801 ! X = MAX ( 0.75 * DEP2(IGC) * KMESPC , 0.5 )
1802  x = max( 0.75 * dep2(kcgrd(1)) * kmespc , 0.5 )
1803  x2 = max( -1.e15, snlcs3*x)
1804  cons = snlc1 * ( 1. + snlcs1/x * (1.-snlcs2*x) * exp(x2))
1805  jacobi = 2. * pi_w
1806 !
1807 ! *** extend the area with action density at periodic boundaries ***
1808 !
1809  DO iddum = idlow, idhgh
1810  id = mod( iddum - 1 + mdc , mdc ) + 1
1811  DO is=1, msc
1812 ! UE (IS,IDDUM) = AC2(ID,IS,IGC) * SPCSIG(IS) * JACOBI
1813  ue(is,iddum) = ac2(id,is,kcgrd(1)) * spcsig(is) * jacobi
1814  ENDDO
1815  ENDDO
1816 !
1817  DO is = msc+1, ishgh
1818  DO id = idlow, idhgh
1819  ue(is,id) = ue(is-1,id) * fachfr
1820  ENDDO
1821  ENDDO
1822 !
1823 ! *** Calculate (unfolded) interactions ***
1824 ! *** Energy at interacting bins ***
1825 !
1826  DO is = isclw, ischg
1827  DO id = 1, mdc
1828  e00 = ue(is ,id )
1829  ep1 = awg1 * ue(is+isp1,id+idp1) + &
1830  awg2 * ue(is+isp1,id+idp ) + &
1831  awg3 * ue(is+isp ,id+idp1) + &
1832  awg4 * ue(is+isp ,id+idp )
1833  em1 = awg5 * ue(is+ism1,id-idm1) + &
1834  awg6 * ue(is+ism1,id-idm ) + &
1835  awg7 * ue(is+ism ,id-idm1) + &
1836  awg8 * ue(is+ism ,id-idm )
1837  ep2 = awg1 * ue(is+isp1,id-idp1) + &
1838  awg2 * ue(is+isp1,id-idp ) + &
1839  awg3 * ue(is+isp ,id-idp1) + &
1840  awg4 * ue(is+isp ,id-idp )
1841  em2 = awg5 * ue(is+ism1,id+idm1) + &
1842  awg6 * ue(is+ism1,id+idm ) + &
1843  awg7 * ue(is+ism ,id+idm1) + &
1844  awg8 * ue(is+ism ,id+idm )
1845 !
1846 ! Contribution to interactions
1847 !
1848  factor = cons * af11(is) * e00
1849 !
1850  sa1a = e00 * ( ep1*dal1 + em1*dal2 ) * pquad(2)
1851  sa1b = sa1a - ep1*em1*dal3 * pquad(2)
1852  sa2a = e00 * ( ep2*dal1 + em2*dal2 ) * pquad(2)
1853  sa2b = sa2a - ep2*em2*dal3 * pquad(2)
1854 !
1855  sa1(is,id) = factor * sa1b
1856  sa2(is,id) = factor * sa2b
1857 !
1858 ! IF(ITEST >= 100 .AND. TESTFL)THEN
1859 ! WRITE(PRINTF,9002) E00,EP1,EM1,EP2,EM2
1860 !9002 FORMAT (' E00 EP1 EM1 EP2 EM2 :',5E11.4)
1861 ! WRITE(PRINTF,9003) SA1A,SA1B,SA2A,SA2B
1862 !9003 FORMAT (' SA1A SA1B SA2A SA2B :',4E11.4)
1863 ! WRITE(PRINTF,9004) IS,ID,SA1(IS,ID),SA2(IS,ID)
1864 !9004 FORMAT (' IS ID SA1() SA2() :',2I4,2E12.4)
1865 ! WRITE(PRINTF,9005) FACTOR,JACOBI
1866 !9005 FORMAT (' FACTOR JACOBI : ',2E12.4)
1867 ! END IF
1868 !
1869  ENDDO
1870  ENDDO
1871 !
1872 ! *** Fold interactions to side angles -> domain in theta is ***
1873 ! *** periodic ***
1874 !
1875  DO id = 1, idhgh - mdc
1876  id0 = 1 - id
1877  DO is = isclw, ischg
1878  sa1(is,mdc+id) = sa1(is, id )
1879  sa2(is,mdc+id) = sa2(is, id )
1880  sa1(is, id0 ) = sa1(is,mdc+id0)
1881  sa2(is, id0 ) = sa2(is,mdc+id0)
1882  ENDDO
1883  ENDDO
1884 !
1885 ! *** Put source term together (To save space I=IS and ***
1886 ! *** J=MDC is used) ---- ***
1887 !
1888  DO i = 1, msc
1889  sigpi = spcsig(i) * jacobi
1890  DO j = 1, mdc
1891  sfnl(i,j) = - 2. * ( sa1(i,j) + sa2(i,j) ) &
1892  + awg1 * ( sa1(i-isp1,j-idp1) + sa2(i-isp1,j+idp1) ) &
1893  + awg2 * ( sa1(i-isp1,j-idp ) + sa2(i-isp1,j+idp ) ) &
1894  + awg3 * ( sa1(i-isp ,j-idp1) + sa2(i-isp ,j+idp1) ) &
1895  + awg4 * ( sa1(i-isp ,j-idp ) + sa2(i-isp ,j+idp ) ) &
1896  + awg5 * ( sa1(i-ism1,j+idm1) + sa2(i-ism1,j-idm1) ) &
1897  + awg6 * ( sa1(i-ism1,j+idm ) + sa2(i-ism1,j-idm ) ) &
1898  + awg7 * ( sa1(i-ism ,j+idm1) + sa2(i-ism ,j-idm1) ) &
1899  + awg8 * ( sa1(i-ism ,j+idm ) + sa2(i-ism ,j-idm ) )
1900 !
1901 ! *** store value in auxiliary array and use values in ***
1902 ! *** next four sweeps (see subroutine FILNL3) ***
1903 !
1904 ! MEMNL4(J,I,IGC) = SFNL(I,J) / SIGPI
1905  memnl4(j,i,kcgrd(1)) = sfnl(i,j) / sigpi
1906  ENDDO
1907  ENDDO
1908 !
1909 ! *** test output ***
1910 !
1911 ! IF(ITEST >= 50 .AND. TESTFL)THEN
1912 ! WRITE(PRINTF,*)
1913 ! WRITE(PRINTF,*) ' SWSNL3 subroutine '
1914 ! WRITE(PRINTF,9011) IDP, IDP1, IDM, IDM1
1915 !9011 FORMAT (' IDP IDP1 IDM IDM1 :',4I5)
1916 ! WRITE (PRINTF,9013) ISP, ISP1, ISM, ISM1
1917 !9013 FORMAT (' ISP ISP1 ISM ISM1 :',4I5)
1918 ! WRITE (PRINTF,9015) ISLOW, ISHGH, IDLOW, IDHGH
1919 !9015 FORMAT (' ISLOW ISHG IDLOW IDHG :',4I5)
1920 ! WRITE(PRINTF,9016) ISCLW, ISCHG, JACOBI
1921 !9016 FORMAT (' ICLW ICHG JACOBI :',2I5,E12.4)
1922 ! WRITE (PRINTF,9017) AWG1, AWG2, AWG3, AWG4
1923 !9017 FORMAT (' AWG1 AWG2 AWG3 AWG4 :',4E12.4)
1924 ! WRITE (PRINTF,9018) AWG5, AWG6, AWG7, AWG8
1925 !9018 FORMAT (' AWG5 AWG6 AWG7 AWG8 :',4E12.4)
1926 ! WRITE (PRINTF,9019) MSC4MI, MSC4MA, MDC4MI, MDC4MA
1927 !9019 FORMAT (' S4MI S4MA D4MI D4MA :',4I6)
1928 ! WRITE(PRINTF,9020) SNLC1,X,X2,CONS
1929 !9020 FORMAT (' SNLC1 X X2 CONS :',4E12.4)
1930 ! WRITE(PRINTF,9021) DEP2(KCGRD(1)),KMESPC,FACHFR,PI
1931 !9021 FORMAT (' DEPTH KMESPC FACHFR PI:',4E12.4)
1932 ! WRITE(PRINTF,*)
1933 !
1934 ! *** value source term in every bin ***
1935 !
1936 ! IF(ITEST >= 150)THEN
1937 ! DO I=1, MSC
1938 ! DO J=1, MDC
1939 ! WRITE(PRINTF,2006) I,J,MEMNL4(J,I,KCGRD(1)),SFNL(I,J), &
1940 ! SPCSIG(I)
1941 !2006 FORMAT (' I J MEMNL() SFNL() SPCSIG:',2I4,3E12.4)
1942 ! ENDDO
1943 ! ENDDO
1944 ! END IF
1945 ! END IF
1946 !
1947  RETURN
1948 !
1949  END SUBROUTINE swsnl3
1950 
1951 !
1952 !*******************************************************************
1953 !
1954  SUBROUTINE swsnl2 (IDDLOW ,IDDTOP ,WWINT ,WWAWG ,UE , &
1955  SA1 ,ISSTOP ,SA2 ,SPCSIG ,SNLC1 , &
1956  DAL1 ,DAL2 ,DAL3 ,SFNL ,DEP2 , &
1957  KMESPC ,IMATDA ,IMATRA ,FACHFR ,PLNL4S , &
1958  IDCMIN ,IDCMAX ,IG ,CGO ,WWSWG )
1959 !
1960 !*******************************************************************
1961 !
1962 ! Calculate non-linear interaction using the discrete interaction
1963 ! approximation (Hasselmann and Hasselmann 1985; WAMDI group 1988)
1964 !
1965 !*******************************************************************
1966  USE swcomm3
1967  USE swcomm4
1968  USE ocpcomm4
1969  USE m_snl4
1970 ! USE ALL_VARS, ONLY : MT,AC2
1971  USE vars_wave, ONLY : mt,ac2
1972 
1973  IMPLICIT NONE
1974 !
1975  REAL :: SPCSIG(MSC)
1976  INTEGER :: IS ,ID ,I ,J ,IG ,ISHGH , &
1977  ISSTOP ,ISP ,ISP1 ,IDP ,IDP1 ,ISM ,ISM1 , &
1978  IDM ,IDM1 ,ISCLW ,ISCHG , &
1979  IDLOW ,IDHGH ,IDDLOW ,IDDTOP ,IDCLOW ,IDCHGH
1980 !
1981  REAL :: X ,X2 ,CONS ,FACTOR ,SNLCS1 ,SNLCS2 ,SNLCS3 , &
1982  E00 ,EP1 ,EM1 ,EP2 ,EM2 ,SA1A ,SA1B , &
1983  SA2A ,SA2B ,KMESPC ,FACHFR ,AWG1 ,AWG2 ,AWG3 , &
1984  AWG4 ,AWG5 ,AWG6 ,AWG7 ,AWG8 ,DAL1 ,DAL2 , &
1985  DAL3 ,JACOBI ,SIGPI
1986 !
1987  REAL :: DEP2(MT) , &
1988  UE(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
1989  SA1(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
1990  SA2(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
1991  DA1C(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
1992  DA1P(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
1993  DA1M(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
1994  DA2C(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
1995  DA2P(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
1996  DA2M(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
1997  SFNL(MSC4MI:MSC4MA , MDC4MI:MDC4MA) , &
1998  DFNL(MSC4MI:MSC4MA , MDC4MI:MDC4MA) , &
1999  IMATRA(MDC,MSC) , &
2000  IMATDA(MDC,MSC) , &
2001  PLNL4S(MDC,MSC,NPTST) , &
2002  WWAWG(*) , &
2003  CGO(MSC,MICMAX)
2004 !
2005  INTEGER :: WWINT(*) ,IDCMIN(MSC) ,IDCMAX(MSC)
2006  INTEGER :: ISLOW,IIID,IDDUM,ID0
2007  REAL :: SNLC1,SWG1,SWG2,SWG3,SWG4,SWG5,SWG6,SWG7,SWG8
2008  REAL :: WWSWG(*),PI3
2009 
2010 !
2011  LOGICAL PERCIR
2012 !
2013  idp = wwint(1)
2014  idp1 = wwint(2)
2015  idm = wwint(3)
2016  idm1 = wwint(4)
2017  isp = wwint(5)
2018  isp1 = wwint(6)
2019  ism = wwint(7)
2020  ism1 = wwint(8)
2021  islow = wwint(9)
2022  ishgh = wwint(10)
2023  isclw = wwint(11)
2024  ischg = wwint(12)
2025  idlow = wwint(13)
2026  idhgh = wwint(14)
2027 
2028  awg1 = wwawg(1)
2029  awg2 = wwawg(2)
2030  awg3 = wwawg(3)
2031  awg4 = wwawg(4)
2032  awg5 = wwawg(5)
2033  awg6 = wwawg(6)
2034  awg7 = wwawg(7)
2035  awg8 = wwawg(8)
2036 !
2037  swg1 = wwswg(1)
2038  swg2 = wwswg(2)
2039  swg3 = wwswg(3)
2040  swg4 = wwswg(4)
2041  swg5 = wwswg(5)
2042  swg6 = wwswg(6)
2043  swg7 = wwswg(7)
2044  swg8 = wwswg(8)
2045 !
2046 ! *** Initialize auxiliary arrays per gridpoint ***
2047 !
2048  DO id = mdc4mi, mdc4ma
2049  DO is = msc4mi, msc4ma
2050  ue(is,id) = 0.
2051  sa1(is,id) = 0.
2052  sa2(is,id) = 0.
2053 
2054  da1c(is,id) = 0.
2055  da1p(is,id) = 0.
2056  da1m(is,id) = 0.
2057  da2c(is,id) = 0.
2058  da2p(is,id) = 0.
2059  da2m(is,id) = 0.
2060 
2061  sfnl(is,id) = 0.
2062  dfnl(is,id) = 0.
2063 
2064  ENDDO
2065  ENDDO
2066 !
2067 ! *** Calculate prop. constant. ***
2068 ! *** Calculate factor R(X) to calculate the NL wave-wave ***
2069 ! *** interaction for shallow water ***
2070 ! *** SNLC1 = 1/GRAV**4 ***
2071 !
2072  snlcs1 = pquad(3)
2073  snlcs2 = pquad(4)
2074  snlcs3 = pquad(5)
2075 
2076  x = max( 0.75 * dep2(ig) * kmespc , 0.5 )
2077  x2 = max( -1.e15, snlcs3*x)
2078  cons = snlc1 * ( 1. + snlcs1/x * (1.-snlcs2*x) * exp(x2))
2079  jacobi = 2. * pi_w
2080 !
2081 ! *** check whether the spectral domain is periodic in ***
2082 ! *** direction space and if so modify boundaries ***
2083 !
2084  percir = .false.
2085  IF(iddlow == 1 .AND. iddtop == mdc)THEN
2086 ! *** periodic in theta -> spectrum can be folded ***
2087 ! *** (can only occur in presence of a current) ***
2088  idclow = 1
2089  idchgh = mdc
2090  iiid = 0
2091  percir = .true.
2092  ELSE
2093 ! *** different sectors per sweep -> extend range with IIID ***
2094  iiid = max( idm1 , idp1 )
2095  idclow = idlow
2096  idchgh = idhgh
2097  ENDIF
2098 !
2099 ! *** Prepare auxiliary spectrum ***
2100 ! *** set action original spectrum in array UE ***
2101 !
2102 !print*,'swsnl2-1'
2103 !print*,'IDLOW,IDHGH,IIID',IDLOW,IDHGH,IIID
2104 !print*,'low high MDC=',IDLOW - IIID, IDHGH + IIID,MDC
2105 
2106  DO iddum = idlow - iiid , idhgh + iiid
2107  id = mod( iddum - 1 + mdc , mdc ) + 1
2108  DO is = 1, msc
2109 ! print*,'ID,IS,IG=',ID,IS,IG
2110  ue(is,iddum) = ac2(id,is,ig) * spcsig(is) * jacobi
2111  ENDDO
2112  ENDDO
2113 !print*,'swsnl2-2'
2114 !
2115 ! *** set values in the areas 2 for IS > MSC+1 ***
2116 !
2117  DO is = msc+1, ishgh
2118  DO id = idlow - iiid , idhgh + iiid
2119  ue(is,id) = ue(is-1,id) * fachfr
2120  ENDDO
2121  ENDDO
2122 !
2123 ! *** Calculate interactions ***
2124 ! *** Energy at interacting bins ***
2125 !
2126  DO is = isclw, ischg
2127  DO id = idclow , idchgh
2128  e00 = ue(is ,id )
2129  ep1 = awg1 * ue(is+isp1,id+idp1) + &
2130  awg2 * ue(is+isp1,id+idp ) + &
2131  awg3 * ue(is+isp ,id+idp1) + &
2132  awg4 * ue(is+isp ,id+idp )
2133  em1 = awg5 * ue(is+ism1,id-idm1) + &
2134  awg6 * ue(is+ism1,id-idm ) + &
2135  awg7 * ue(is+ism ,id-idm1) + &
2136  awg8 * ue(is+ism ,id-idm )
2137 !
2138  ep2 = awg1 * ue(is+isp1,id-idp1) + &
2139  awg2 * ue(is+isp1,id-idp ) + &
2140  awg3 * ue(is+isp ,id-idp1) + &
2141  awg4 * ue(is+isp ,id-idp )
2142  em2 = awg5 * ue(is+ism1,id+idm1) + &
2143  awg6 * ue(is+ism1,id+idm ) + &
2144  awg7 * ue(is+ism ,id+idm1) + &
2145  awg8 * ue(is+ism ,id+idm )
2146 !
2147 ! *** Contribution to interactions ***
2148 ! *** CONS is the shallow water factor for the NL interact. ***
2149 !
2150  factor = cons * af11(is) * e00
2151 !
2152  sa1a = e00 * ( ep1*dal1 + em1*dal2 ) * pquad(2)
2153  sa1b = sa1a - ep1*em1*dal3 * pquad(2)
2154  sa2a = e00 * ( ep2*dal1 + em2*dal2 ) * pquad(2)
2155  sa2b = sa2a - ep2*em2*dal3 * pquad(2)
2156 !
2157  sa1(is,id) = factor * sa1b
2158  sa2(is,id) = factor * sa2b
2159 
2160  da1c(is,id) = cons*af11(is)*(sa1a+sa1b)
2161  da1p(is,id) = factor*(dal1*e00-dal3*em1)*pquad(2)
2162  da1m(is,id) = factor*(dal2*e00-dal3*ep1)*pquad(2)
2163 
2164  da2c(is,id) = cons*af11(is)*(sa2a+sa2b)
2165  da2p(is,id) = factor*(dal1*e00-dal3*em2)*pquad(2)
2166  da2m(is,id) = factor*(dal2*e00-dal3*ep2)*pquad(2)
2167 
2168  ENDDO
2169  ENDDO
2170 !
2171 ! *** Fold interactions to side angles if spectral domain ***
2172 ! *** is periodic in directional space ***
2173 !
2174  IF(percir)THEN
2175  DO id = 1, idhgh - mdc
2176  id0 = 1 - id
2177  DO is = isclw, ischg
2178  sa1(is,mdc+id) = sa1(is , id )
2179  sa2(is,mdc+id) = sa2(is , id )
2180  sa1(is, id0 ) = sa1(is , mdc+id0)
2181  sa2(is, id0 ) = sa2(is , mdc+id0)
2182 
2183  da1c(is,mdc+id) = da1c(is,id)
2184  da1p(is,mdc+id) = da1p(is,id)
2185  da1m(is,mdc+id) = da1m(is,id)
2186  da1c(is,id0) = da1c(is,mdc+id0)
2187  da1p(is,id0) = da1p(is,mdc+id0)
2188  da1m(is,id0) = da1m(is,mdc+id0)
2189 
2190  da2c(is,mdc+id) = da2c(is,id)
2191  da2p(is,mdc+id) = da2p(is,id)
2192  da2m(is,mdc+id) = da2m(is,id)
2193  da2c(is,id0) = da2c(is,mdc+id0)
2194  da2p(is,id0) = da2p(is,mdc+id0)
2195  da2m(is,id0) = da2m(is,mdc+id0)
2196  ENDDO
2197  ENDDO
2198  ENDIF
2199 !
2200 ! *** Put source term together (To save space I=IS and J=ID ***
2201 ! *** is used) ***
2202 !
2203  pi3 = (2.0*pi_w)**3
2204  DO i = 1, isstop
2205  sigpi = spcsig(i) * jacobi
2206 
2207  DO j = idcmin(i), idcmax(i)
2208  id = mod( j - 1 + mdc , mdc ) + 1
2209  sfnl(i,id) = - 2. * ( sa1(i,j) + sa2(i,j) ) &
2210  + awg1 * ( sa1(i-isp1,j-idp1) + sa2(i-isp1,j+idp1) ) &
2211  + awg2 * ( sa1(i-isp1,j-idp ) + sa2(i-isp1,j+idp ) ) &
2212  + awg3 * ( sa1(i-isp ,j-idp1) + sa2(i-isp ,j+idp1) ) &
2213  + awg4 * ( sa1(i-isp ,j-idp ) + sa2(i-isp ,j+idp ) ) &
2214  + awg5 * ( sa1(i-ism1,j+idm1) + sa2(i-ism1,j-idm1) ) &
2215  + awg6 * ( sa1(i-ism1,j+idm ) + sa2(i-ism1,j-idm ) ) &
2216  + awg7 * ( sa1(i-ism ,j+idm1) + sa2(i-ism ,j-idm1) ) &
2217  + awg8 * ( sa1(i-ism ,j+idm ) + sa2(i-ism ,j-idm ) )
2218 
2219  dfnl(i,id) = - 2. * ( da1c(i,j) + da2c(i,j) ) &
2220  + swg1 * ( da1p(i-isp1,j-idp1) + da2p(i-isp1,j+idp1) ) &
2221  + swg2 * ( da1p(i-isp1,j-idp ) + da2p(i-isp1,j+idp ) ) &
2222  + swg3 * ( da1p(i-isp ,j-idp1) + da2p(i-isp ,j+idp1) ) &
2223  + swg4 * ( da1p(i-isp ,j-idp ) + da2p(i-isp ,j+idp ) ) &
2224  + swg5 * ( da1m(i-ism1,j+idm1) + da2m(i-ism1,j-idm1) ) &
2225  + swg6 * ( da1m(i-ism1,j+idm ) + da2m(i-ism1,j-idm ) ) &
2226  + swg7 * ( da1m(i-ism ,j+idm1) + da2m(i-ism ,j-idm1) ) &
2227  + swg8 * ( da1m(i-ism ,j+idm ) + da2m(i-ism ,j-idm ) )
2228 
2229 !
2230 ! *** store results in rhv ***
2231 ! *** store results in rhs and main diagonal according ***
2232 ! *** to Patankar-rules ***
2233 !
2234  IF(testfl) plnl4s(id,i,iptst) = sfnl(i,id) / sigpi
2235  imatra(id,i) = imatra(id,i) + sfnl(i,id) / sigpi
2236  imatda(id,i) = imatda(id,i) + dfnl(i,id)
2237 
2238  ENDDO
2239  ENDDO
2240 
2241  RETURN
2242  END SUBROUTINE swsnl2
2243 
2244 !
2245 !********************************************************************
2246 !
2247  SUBROUTINE swsnl1 (WWINT ,WWAWG ,WWSWG ,IDCMIN ,IDCMAX , &
2248  UE ,SA1 ,SA2 ,DA1C ,DA1P , &
2249  DA1M ,DA2C ,DA2P ,DA2M ,SPCSIG , &
2250  SNLC1 ,KMESPC ,FACHFR ,ISSTOP ,DAL1 , &
2251  DAL2 ,DAL3 ,SFNL ,DSNL ,DEP2 , &
2252  AC2 ,IMATDA ,IMATRA ,PLNL4S ,PLNL4D , &
2253  IDDLOW ,IDDTOP )
2254 
2255 !(This subroutine has not been tested yet and not useful any more)
2256 !
2257 !********************************************************************
2258 !
2259  USE swcomm3
2260  USE swcomm4
2261  USE ocpcomm4
2262  USE m_snl4
2263  USE all_vars, ONLY : mt
2264 !
2265 ! --|-----------------------------------------------------------|--
2266 ! | Delft University of Technology |
2267 ! | Faculty of Civil Engineering |
2268 ! | Fluid Mechanics Section |
2269 ! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
2270 ! | |
2271 ! | Programmers: H.L. Tolman, R.C. Ris |
2272 ! --|-----------------------------------------------------------|--
2273 !
2274 !
2275 ! SWAN (Simulating WAves Nearshore); a third generation wave model
2276 ! Copyright (C) 2004-2005 Delft University of Technology
2277 !
2278 ! This program is free software; you can redistribute it and/or
2279 ! modify it under the terms of the GNU General Public License as
2280 ! published by the Free Software Foundation; either version 2 of
2281 ! the License, or (at your option) any later version.
2282 !
2283 ! This program is distributed in the hope that it will be useful,
2284 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
2285 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2286 ! GNU General Public License for more details.
2287 !
2288 ! A copy of the GNU General Public License is available at
2289 ! http://www.gnu.org/copyleft/gpl.html#SEC3
2290 ! or by writing to the Free Software Foundation, Inc.,
2291 ! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
2292 !
2293 !
2294 ! 0. Authors
2295 !
2296 ! 30.72: IJsbrand Haagsma
2297 ! 40.13: Nico Booij
2298 ! 40.17: IJsbrand Haagsma
2299 ! 40.23: Marcel Zijlema
2300 ! 40.41: Marcel Zijlema
2301 !
2302 ! 1. Updates
2303 !
2304 ! 30.72, Feb. 98: Introduced generic names XCGRID, YCGRID and SPCSIG for SWAN
2305 ! 40.17, Dec. 01: Implentation of Multiple DIA
2306 ! 40.23, Aug. 02: some corrections
2307 ! 40.41, Oct. 04: common blocks replaced by modules, include files removed
2308 !
2309 ! 2. Purpose
2310 !
2311 ! Calculate non-linear interaction using the discrete interaction
2312 ! approximation (Hasselmann and Hasselmann 1985; WAMDI group 1988),
2313 ! including the diagonal term for the implicit integration.
2314 !
2315 ! The interactions are calculated for all bin's that fall
2316 ! within a sweep. No additional auxiliary array is required (see
2317 ! SWSNL3)
2318 !
2319 ! 3. Method
2320 !
2321 ! Discrete interaction approximation.
2322 !
2323 ! Since the domain in directional domain is by definition not
2324 ! periodic, the spectral space can not beforehand
2325 ! folded to the side angles. This can only be done if the
2326 ! full circle has to be calculated
2327 !
2328 !
2329 ! Frequencies -->
2330 ! +---+---------------------+---------+- IDHGH
2331 ! d | 3 : 2 : 2 |
2332 ! i + - + - - - - - - - - - - + - - - - +- MDC
2333 ! r | : : |
2334 ! e | 3 : original spectrum : 1 |
2335 ! c | : : |
2336 ! t. + - + - - - - - - - - - - + - - - - +- 1
2337 ! | 3 : 2 : 2 |
2338 ! +---+---------------------+---------+- IDLOW
2339 ! | | | ^ |
2340 ! ISLOW 1 MSC | ISHGH
2341 ! ^ |
2342 ! | |
2343 ! ISCLW ISCHG
2344 ! lowest discrete highest discrete
2345 ! central bin central bin
2346 !
2347 ! 1 : Extra tail added beyond MSC
2348 ! 2 : Spectrum copied outside ID range
2349 ! 3 : Empty bins at low frequencies
2350 !
2351 ! ISLOW = 1 + ISM1
2352 ! ISHGH = MSC + ISP1 - ISM1
2353 ! ISCLW = 1
2354 ! ISCHG = MSC - ISM1
2355 ! IDLOW = IDDLOW - MAX(IDM1,IDP1)
2356 ! IDHGH = IDDTOP + MAX(IDM1,IDP1)
2357 !
2358 ! For the meaning of the counters on the right hand side of the
2359 ! above equations see section 4.
2360 !
2361 ! 4. Argument variables
2362 !
2363 ! SPCSIG: Relative frequencies in computational domain in sigma-space
2364 !
2365  REAL SPCSIG(MSC)
2366 !
2367 ! Data in PARAMETER statements :
2368 ! ----------------------------------------------------------------
2369 ! DAL1 Real LAMBDA dependend weight factors (see FAC4WW)
2370 ! DAL2 Real
2371 ! DAL3 Real
2372 ! ITHP, ITHP1, ITHM, ITHM1, IFRP, IFRP1, IFRM, IFRM1
2373 ! Int. Counters of interpolation point relative to
2374 ! central bin, see figure below (set in FAC4WW).
2375 ! NFRLOW, NFRHGH, NFRCHG, NTHLOW, NTHHGH
2376 ! Int. Range of calculations, see section 2.
2377 ! AF11 R.A. Scaling array (Freq**11).
2378 ! AWGn Real Interpolation weights, see numbers in fig.
2379 ! SWGn Real Id. squared.
2380 ! UE R.A. "Unfolded" spectrum.
2381 ! SA1 R.A. Interaction constribution of first and second
2382 ! SA2 R.A. quadr. respectively (unfolded space).
2383 ! DA1C, DA1P, DA1M, DA2C, DA2P, DA2M
2384 ! R.A. Idem for diagonal matrix.
2385 ! PERCIR full circle or sector
2386 ! ----------------------------------------------------------------
2387 !
2388 ! Relative offsets of interpolation points around central bin
2389 ! "#" and corresponding numbers of AWGn :
2390 !
2391 ! ISM1 ISM
2392 ! 5 7 T |
2393 ! IDM1 +------+ H +
2394 ! | | E | ISP ISP1
2395 ! | \ | T | 3 1
2396 ! IDM +------+ A + +---------+ IDP1
2397 ! 6 \8 | | |
2398 ! | | / |
2399 ! \ + +---------+ IDP
2400 ! | /4 2
2401 ! \ | /
2402 ! -+-----+------+-------#--------+---------+----------+
2403 ! | FREQ.
2404 !
2405 ! 7. Common blocks used
2406 !
2407 !
2408 ! 8. Subroutines used
2409 !
2410 ! ---
2411 !
2412 ! 9. Subroutines calling
2413 !
2414 ! SOURCE (in SWANCOM1)
2415 !
2416 ! 12. Structure
2417 !
2418 ! -------------------------------------------
2419 ! Initialisations.
2420 ! Calculate proportionality constant.
2421 ! Prepare auxiliary spectrum.
2422 ! Calculate interactions :
2423 ! -----------------------------------------
2424 ! Energy at interacting bins
2425 ! Contribution to interactions
2426 ! Fold interactions to side angles
2427 ! -----------------------------------------
2428 ! Put source term together
2429 ! -------------------------------------------
2430 !
2431 ! 13. Source text
2432 !
2433 !*************************************************************
2434 !
2435  INTEGER IS ,ID ,I ,J , &
2436  ISHGH ,IDLOW ,ISP ,ISP1 ,IDP ,IDP1 , &
2437  ISM ,ISM1 ,IDHGH ,IDM ,IDM1 ,ISCLW , &
2438  ISCHG ,IDDLOW ,IDDTOP
2439 !
2440  REAL X ,X2 ,CONS ,FACTOR ,SNLCS1 ,SNLCS2 ,SNLCS3, &
2441  E00 ,EP1 ,EM1 ,EP2 ,EM2 ,SA1A ,SA1B , &
2442  SA2A ,SA2B ,KMESPC ,FACHFR ,AWG1 ,AWG2 ,AWG3 , &
2443  AWG4 ,AWG5 ,AWG6 ,AWG7 ,AWG8 ,DAL1 ,DAL2 , &
2444  DAL3 ,SNLC1 ,SWG1 ,SWG2 ,SWG3 ,SWG4 ,SWG5 , &
2445  SWG6 ,SWG7 ,SWG8 ,JACOBI ,SIGPI
2446 !
2447  REAL AC2(MDC,MSC,0:MT) , &
2448  DEP2(MT) , &
2449  UE(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
2450  SA1(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
2451  SA2(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
2452  DA1C(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
2453  DA1P(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
2454  DA1M(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
2455  DA2C(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
2456  DA2P(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
2457  DA2M(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
2458  SFNL(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
2459  DSNL(MSC4MI:MSC4MA , MDC4MI:MDC4MA ) , &
2460  IMATDA(MDC,MSC) , &
2461  IMATRA(MDC,MSC) , &
2462  PLNL4S(MDC,MSC,NPTST) , &
2463  PLNL4D(MDC,MSC,NPTST) , &
2464  WWAWG(*) , &
2465  WWSWG(*)
2466 !
2467  INTEGER IDCMIN(MSC),IDCMAX(MSC),WWINT(*)
2468 !
2469  LOGICAL PERCIR
2470 !
2471 ! SAVE IENT
2472 ! DATA IENT/0/
2473 ! IF (LTRACE) CALL STRACE (IENT,'SWSNL1')
2474 !
2475  idp = wwint(1)
2476  idp1 = wwint(2)
2477  idm = wwint(3)
2478  idm1 = wwint(4)
2479  isp = wwint(5)
2480  isp1 = wwint(6)
2481  ism = wwint(7)
2482  ism1 = wwint(8)
2483  islow = wwint(9)
2484  ishgh = wwint(10)
2485  isclw = wwint(11)
2486  ischg = wwint(12)
2487  idlow = wwint(13)
2488  idhgh = wwint(14)
2489 !
2490  awg1 = wwawg(1)
2491  awg2 = wwawg(2)
2492  awg3 = wwawg(3)
2493  awg4 = wwawg(4)
2494  awg5 = wwawg(5)
2495  awg6 = wwawg(6)
2496  awg7 = wwawg(7)
2497  awg8 = wwawg(8)
2498 !
2499  swg1 = wwswg(1)
2500  swg2 = wwswg(2)
2501  swg3 = wwswg(3)
2502  swg4 = wwswg(4)
2503  swg5 = wwswg(5)
2504  swg6 = wwswg(6)
2505  swg7 = wwswg(7)
2506  swg8 = wwswg(8)
2507 !
2508 ! *** Initialize auxiliary arrays per gridpoint ***
2509 !
2510  DO id = mdc4mi, mdc4ma
2511  DO is = msc4mi, msc4ma
2512  ue(is,id) = 0.
2513  sa1(is,id) = 0.
2514  sa2(is,id) = 0.
2515  sfnl(is,id) = 0.
2516  da1c(is,id) = 0.
2517  da1p(is,id) = 0.
2518  da1m(is,id) = 0.
2519  da2c(is,id) = 0.
2520  da2p(is,id) = 0.
2521  da2m(is,id) = 0.
2522  dsnl(is,id) = 0.
2523  ENDDO
2524  ENDDO
2525 !
2526 ! *** Calculate factor R(X) to calculate the NL wave-wave ***
2527 ! *** interaction for shallow water ***
2528 ! *** SNLC1 = 1/GRAV**4 ***
2529 !
2530  snlcs1 = pquad(3)
2531  snlcs2 = pquad(4)
2532  snlcs3 = pquad(5)
2533 ! X = MAX ( 0.75 * DEP2(IGC) * KMESPC , 0.5 )
2534  x = max( 0.75 * dep2(kcgrd(1)) * kmespc , 0.5 )
2535  x2 = max( -1.e15, snlcs3*x)
2536  cons = snlc1 * ( 1. + snlcs1/x * (1.-snlcs2*x) * exp(x2))
2537  jacobi = 2. * pi_w
2538 !
2539 ! *** check whether the spectral domain is periodic in ***
2540 ! *** directional space and if so, modify boundaries ***
2541 !
2542  percir = .false.
2543  IF(iddlow == 1 .AND. iddtop == mdc)THEN
2544 ! *** periodic in theta -> spectrum can be folded ***
2545 ! *** (can only be present in presence of a current) ***
2546  idclow = 1
2547  idchgh = mdc
2548  iiid = 0
2549  percir = .true.
2550  ELSE
2551 ! *** different sectors per sweep -> extend range with IIID ***
2552  iiid = max( idm1 , idp1 )
2553  idclow = idlow
2554  idchgh = idhgh
2555  ENDIF
2556 !
2557 ! *** Prepare auxiliary spectrum ***
2558 ! *** set action original spectrum in array UE ***
2559 !
2560  DO iddum = idlow - iiid, idhgh + iiid
2561  id = mod( iddum - 1 + mdc , mdc ) + 1
2562  DO is = 1, msc
2563 ! UE(IS,IDDUM) = AC2(ID,IS,IGC) * SPCSIG(IS) * JACOBI
2564  ue(is,iddum) = ac2(id,is,kcgrd(1)) * spcsig(is) * jacobi
2565  ENDDO
2566  ENDDO
2567 !
2568 ! *** set values in area 2 for IS > MSC+1 ***
2569 !
2570  DO is = msc+1, ishgh
2571  DO id = idlow - iiid , idhgh + iiid
2572  ue(is,id) = ue(is-1,id) * fachfr
2573  ENDDO
2574  ENDDO
2575 !
2576 ! *** Calculate interactions ***
2577 ! *** Energy at interacting bins ***
2578 !
2579  DO is = isclw, ischg
2580  DO id = idclow, idchgh
2581  e00 = ue(is ,id )
2582  ep1 = awg1 * ue(is+isp1,id+idp1) + &
2583  awg2 * ue(is+isp1,id+idp ) + &
2584  awg3 * ue(is+isp ,id+idp1) + &
2585  awg4 * ue(is+isp ,id+idp )
2586  em1 = awg5 * ue(is+ism1,id-idm1) + &
2587  awg6 * ue(is+ism1,id-idm ) + &
2588  awg7 * ue(is+ism ,id-idm1) + &
2589  awg8 * ue(is+ism ,id-idm )
2590 !
2591  ep2 = awg1 * ue(is+isp1,id-idp1) + &
2592  awg2 * ue(is+isp1,id-idp ) + &
2593  awg3 * ue(is+isp ,id-idp1) + &
2594  awg4 * ue(is+isp ,id-idp )
2595  em2 = awg5 * ue(is+ism1,id+idm1) + &
2596  awg6 * ue(is+ism1,id+idm ) + &
2597  awg7 * ue(is+ism ,id+idm1) + &
2598  awg8 * ue(is+ism ,id+idm )
2599 !
2600 ! *** Contribution to interactions ***
2601 ! *** CONS is the shallow water factor for the NL interact. ***
2602 !
2603  factor = cons * af11(is) * e00
2604 !
2605  sa1a = e00 * ( ep1*dal1 + em1*dal2 ) * pquad(2)
2606  sa1b = sa1a - ep1*em1*dal3 * pquad(2)
2607  sa2a = e00 * ( ep2*dal1 + em2*dal2 ) * pquad(2)
2608  sa2b = sa2a - ep2*em2*dal3 * pquad(2)
2609 !
2610  sa1(is,id) = factor * sa1b
2611  sa2(is,id) = factor * sa2b
2612 !
2613 ! IF(ITEST >= 100 .AND. TESTFL)THEN
2614 ! WRITE(PRINTF,9002) E00,EP1,EM1,EP2,EM2
2615 !9002 FORMAT (' E00 EP1 EM1 EP2 EM2 :',5E11.4)
2616 ! WRITE(PRINTF,9003) SA1A,SA1B,SA2A,SA2B
2617 !9003 FORMAT (' SA1A SA1B SA2A SA2B :',4E11.4)
2618 ! WRITE(PRINTF,9004) IS,ID,SA1(IS,ID),SA2(IS,ID)
2619 !9004 FORMAT (' IS ID SA1() SA2() :',2I4,2E12.4)
2620 ! WRITE(PRINTF,9005) FACTOR
2621 !9005 FORMAT (' FACTOR : ',E12.4)
2622 ! END IF
2623 !
2624  da1c(is,id) = cons * af11(is) * ( sa1a + sa1b )
2625  da1p(is,id) = factor * ( dal1*e00 - dal3*em1 ) * pquad(2)
2626  da1m(is,id) = factor * ( dal2*e00 - dal3*ep1 ) * pquad(2)
2627 !
2628  da2c(is,id) = cons * af11(is) * ( sa2a + sa2b )
2629  da2p(is,id) = factor * ( dal1*e00 - dal3*em2 ) * pquad(2)
2630  da2m(is,id) = factor * ( dal2*e00 - dal3*ep2 ) * pquad(2)
2631  ENDDO
2632  ENDDO
2633 !
2634 ! *** Fold interactions to side angles if spectral domain ***
2635 ! *** is periodic in directional space ***
2636 !
2637  IF(percir)THEN
2638  DO id = 1, idhgh - mdc
2639  id0 = 1 - id
2640  DO is = isclw, ischg
2641  sa1(is,mdc+id) = sa1(is, id )
2642  sa2(is,mdc+id) = sa2(is, id )
2643  da1c(is,mdc+id) = da1c(is, id )
2644  da1p(is,mdc+id) = da1p(is, id )
2645  da1m(is,mdc+id) = da1m(is, id )
2646  da2c(is,mdc+id) = da2c(is, id )
2647  da2p(is,mdc+id) = da2p(is, id )
2648  da2m(is,mdc+id) = da2m(is, id )
2649 !
2650  sa1(is, id0 ) = sa1(is, mdc+id0)
2651  sa2(is, id0 ) = sa2(is, mdc+id0)
2652  da1c(is, id0 ) = da1c(is, mdc+id0)
2653  da1p(is, id0 ) = da1p(is, mdc+id0)
2654  da1m(is, id0 ) = da1m(is, mdc+id0)
2655  da2c(is, id0 ) = da2c(is, mdc+id0)
2656  da2p(is, id0 ) = da2p(is, mdc+id0)
2657  da2m(is, id0 ) = da2m(is, mdc+id0)
2658  ENDDO
2659  ENDDO
2660  ENDIF
2661 !
2662 ! *** Put source term together (To save space I=IS and J=ID ***
2663 ! *** is used) ***
2664 !
2665  pi3 = (2. * pi_w)**3
2666  DO i = 1, isstop
2667  sigpi = spcsig(i) * jacobi
2668  DO j = idcmin(i), idcmax(i)
2669  id = mod( j - 1 + mdc , mdc ) + 1
2670  sfnl(i,id) = - 2. * ( sa1(i,j) + sa2(i,j) ) &
2671  + awg1 * ( sa1(i-isp1,j-idp1) + sa2(i-isp1,j+idp1) ) &
2672  + awg2 * ( sa1(i-isp1,j-idp ) + sa2(i-isp1,j+idp ) ) &
2673  + awg3 * ( sa1(i-isp ,j-idp1) + sa2(i-isp ,j+idp1) ) &
2674  + awg4 * ( sa1(i-isp ,j-idp ) + sa2(i-isp ,j+idp ) ) &
2675  + awg5 * ( sa1(i-ism1,j+idm1) + sa2(i-ism1,j-idm1) ) &
2676  + awg6 * ( sa1(i-ism1,j+idm ) + sa2(i-ism1,j-idm ) ) &
2677  + awg7 * ( sa1(i-ism ,j+idm1) + sa2(i-ism ,j-idm1) ) &
2678  + awg8 * ( sa1(i-ism ,j+idm ) + sa2(i-ism ,j-idm ) )
2679 !
2680  dsnl(i,id) = - 2. * ( da1c(i,j) + da2c(i,j) ) &
2681  + swg1 * ( da1p(i-isp1,j-idp1) + da2p(i-isp1,j+idp1) ) &
2682  + swg2 * ( da1p(i-isp1,j-idp ) + da2p(i-isp1,j+idp ) ) &
2683  + swg3 * ( da1p(i-isp ,j-idp1) + da2p(i-isp ,j+idp1) ) &
2684  + swg4 * ( da1p(i-isp ,j-idp ) + da2p(i-isp ,j+idp ) ) &
2685  + swg5 * ( da1m(i-ism1,j+idm1) + da2m(i-ism1,j-idm1) ) &
2686  + swg6 * ( da1m(i-ism1,j+idm ) + da2m(i-ism1,j-idm ) ) &
2687  + swg7 * ( da1m(i-ism ,j+idm1) + da2m(i-ism ,j-idm1) ) &
2688  + swg8 * ( da1m(i-ism ,j+idm ) + da2m(i-ism ,j-idm ) )
2689 !
2690 ! *** store results in IMATDA and IMATRA ***
2691 !
2692  IF(testfl) THEN
2693  plnl4s(id,i,iptst) = sfnl(i,id) / sigpi
2694  plnl4d(id,i,iptst) = -1. * dsnl(i,id) / pi3
2695  END IF
2696 !
2697  imatra(id,i) = imatra(id,i) + sfnl(i,id) / sigpi
2698  imatda(id,i) = imatda(id,i) - dsnl(i,id) / pi3
2699 !
2700 ! IF(ITEST >= 90 .AND. TESTFL) THEN
2701 ! WRITE(PRINTF,9006) I,J,SFNL(I,ID),DSNL(I,ID),SPCSIG(I) 30.72
2702 !9006 FORMAT (' IS ID SFNL DSNL SPCSIG:',2I4,3E12.4)
2703 ! END IF
2704 !
2705  ENDDO
2706  ENDDO
2707 !
2708 
2709  RETURN
2710  END SUBROUTINE swsnl1
2711 
2712 !=====================================================================
2713 ! (The functions and subroutines below have not been tested yet)
2714 !=====================================================================
2715 
2716  REAL FUNCTION WAVE(D)
2717 
2718 ! Transforms nondimensional Sqr(w)d/g into nondimensional kd using the
2719 ! dispersion relation and an iterative method
2720 
2721  IF(d-10.) 2,2,1
2722  1 xx=d
2723  GO TO 6
2724  2 IF(d-1.) 3,4,4
2725  3 x=sqrt(d)
2726  GO TO 5
2727  4 x=d
2728  5 cothx=1.0/tanh(x)
2729  xx=x-(x-d*cothx)/(1.+d*(cothx**2-1.))
2730  e=1.-xx/x
2731  x=xx
2732  IF(abs(e)-0.0005) 6,5,5
2733  6 wave=xx
2734  RETURN
2735  END FUNCTION wave
2736 
2737 
2738  SUBROUTINE cgcmp(G,AK,H,CG)
2739 
2740 ! Calculates group velocity Cg based on depth and wavenumber
2741 ! Includes a deep water limit for kd > 10.
2742 
2743 ! G : gravitational acceleration
2744 ! AK : wave number
2745 ! H : depth
2746 ! CG : group velocity
2747 ! AKH : depth x Wave number (kd)
2748 ! RGK : square root of gk
2749 ! SECH2 : the square of the secant Hyperbolic of kd
2750 ! = 1 - TANH(kd)**2
2751 
2752 ! Calculation of group velocity Cg
2753 
2754  akh=ak*h
2755  IF(akh <= 10.)THEN
2756 
2757 ! Shallow water:
2758 
2759 ! Cg = (1/2) (1 + (2 kd) / Sinh (2 kd)) (w / k)
2760 ! = (1/2) (1 + (kd SECH2) / Tanh (kd)) (w / k)
2761 ! = (1/2) (1 + (kd SECH2) / Tanh (kd)) (Sqrt (gk Tanh(kd)) / k)
2762 ! = (1/2) (Sqrt (gk) (Sqrt(Tanh (kd)) + kd SECH2 / Sqrt (Tanh (kd)) / k
2763 ! = (1/2) (Sqrt (k) g (Tanh (kd) + kd SECH2) / (k Sqrt (g Tanh (kd)))
2764 ! = (g Tanh(kd) + gkd SECH2) / (2 Sqrt(gk Tanh(kd))
2765 
2766  sech2=1.-tanh(akh)**2
2767  cg=(g*akh*sech2+g*tanh(akh))/(2.*sqrt(g*ak*tanh(akh)))
2768 
2769  ELSE
2770 
2771 ! Deep water:
2772 
2773 ! Cg = w / (2 k)
2774 ! = g / (2 Sqrt (gk))
2775 
2776  rgk=sqrt(g*ak)
2777  cg=g/(2.*rgk)
2778 
2779  ENDIF
2780 !
2781  RETURN
2782 
2783  END SUBROUTINE cgcmp
2784 
2785  SUBROUTINE preset(LMAX,N,IW3L,IW4,N2,G,H,W4,AK4,W,DQ,DQ2,DT,DT2,P)
2786 
2787 !(This subroutine has not been tested yet)
2788 
2789 ! T1A : Angle between theta_1 and theta_a
2790 ! T2A : Angle between theta_2 and theta_a
2791 ! T34 : Angle between theta_3 and theta_4
2792 ! WA : wa = w1 + w2 = w3 + w4 (resonance conditions)
2793 ! AKA : absolute value of ka = k1 + k2 = k3 + k4 (resonance conditions)
2794 ! TA : Angle theta_a representing vector ka
2795 
2796  REAL :: W(LMAX)
2797 !
2798 ! ================
2799  DO 120 it34=1,n2
2800 ! ================
2801 !
2802  t34=dt*(it34-1)
2803  dt3=dt
2804  IF(it34 == 1 .OR. it34 == n2) dt3=dt2
2805 !
2806 ! ===================
2807  DO 110 iw3=iw3l,iw4
2808 ! ===================
2809 !
2810  w3=w(iw3)
2811  ak3=wave(w3**2*h/g)/h
2812 !
2813  wa=w3+w4
2814  aka=sqrt(ak3*ak3+ak4*ak4+2.*ak3*ak4*cos(t34))
2815  ta=atan2(ak3*sin(t34),ak3*cos(t34)+ak4)
2816  r=sqrt(g*aka*tanh(aka*h/2.))/wa-1./sqrt(2.)
2817 !
2818  tl=0.
2819  acs=aka/(2.*wave(wa**2*h/(4.*g))/h)
2820  acs=sign(1.,acs)*min(1.,abs(acs))
2821  IF(r < 0.) tl=acos(acs)
2822  it1s=nint(tl/dt)+1
2823 !
2824 ! ===================
2825  DO 100 it1a=it1s,n2
2826 ! ===================
2827 !
2828  t1a=dt*(it1a-1)
2829 !
2830  dt1=dt
2831  IF(it1a == 1 .OR. it1a == n2) dt1=dt2
2832 !
2833 ! -----------------------------------------------------------------
2834  CALL findw1w2(lmax,iw3,w,g,h,w1,w2,w3,wa,ak1,ak2,aka,t1a,t2a,ind)
2835 ! -----------------------------------------------------------------
2836  IF(ind < 0) GO TO 100
2837 !
2838  IF(w1 <= w3 .AND. w3 <= w4 .AND. w4 <= w2)THEN
2839 !JQI CALL KERNEL(N,G,H,W1,W2,W3,W4,AK1,AK2,AK3,AK4,AKA, &
2840 !JQI T1A,T2A,T34,TA,DQ,DT,DT1,DT3,P)
2841  ENDIF
2842 !
2843 ! ============
2844  100 CONTINUE
2845  110 CONTINUE
2846  120 CONTINUE
2847 ! ============
2848 !
2849  RETURN
2850  END SUBROUTINE preset
2851 
2852  SUBROUTINE findw1w2(LMAX,IW3,W,G,H,W1,W2,W3,WA,AK1,AK2,AKA,T1A,T2A,IND)
2853  REAL :: W(LMAX)
2854 !
2855  ind=0
2856  eps=0.0005
2857 !
2858  x1=0.005
2859  x2=w3
2860 !
2861 ! ---------------------------------------
2862 110 CALL fdw1(g,h,aka,t1a,wa,x1,x2,x,eps,m)
2863 ! ---------------------------------------
2864 !
2865  IF(m.EQ.0) GO TO 999
2866 !
2867  w1=x
2868 !
2869  ak1=wave(w1**2*h/g)/h
2870  ak2=sqrt(aka*aka+ak1*ak1-2.*aka*ak1*cos(t1a))
2871  w2=sqrt(g*ak2*tanh(ak2*h))
2872  IF(w1.GT.w2) GO TO 999
2873  t2a=atan2(-ak1*sin(t1a),aka-ak1*cos(t1a))
2874  RETURN
2875 !
2876 999 ind=-999
2877  RETURN
2878  END SUBROUTINE findw1w2
2879 
2880  SUBROUTINE fdw1(G,H,AKA,T1A,WA,X1,X2,X,EPS,M)
2881 !
2882  m=0
2883  f1=funcw1(g,h,x1,aka,t1a,wa)
2884  f2=funcw1(g,h,x2,aka,t1a,wa)
2885 20 IF(f1*f2) 18,19,21
2886 18 m=m+1
2887  x=x2-((x2-x1)/(f2-f1)*f2)
2888  IF((abs(x-x2)/abs(x))-eps) 22,23,23
2889 23 IF((abs(x-x1)/abs(x))-eps) 22,24,24
2890 22 RETURN
2891 24 f=funcw1(g,h,x,aka,t1a,wa)
2892  fm=f*f1
2893  IF(fm) 25,22,26
2894 25 x2=x
2895  f2=f
2896  GO TO 20
2897 26 x1=x
2898  f1=f
2899  GO TO 20
2900 19 m=m+1
2901  IF(f1) 17,16,17
2902 16 x=x1
2903  GO TO 22
2904 17 x=x2
2905  GO TO 22
2906 21 m=0
2907  RETURN
2908  END SUBROUTINE fdw1
2909 
2910  FUNCTION funcw1(G,H,X,AKA,T1A,WA)
2911 !
2912  ak1=wave(x**2*h/g)/h
2913  ak2=sqrt(aka*aka+ak1*ak1-2.*aka*ak1*cos(t1a))
2914  funcw1=wa-sqrt(g*ak1*tanh(ak1*h))-sqrt(g*ak2*tanh(ak2*h))
2915 !
2916  RETURN
2917  END FUNCTION funcw1
vars_wave
Definition: mod_main_wave.f90:14
swcomm4::testfl
logical testfl
Definition: swmod1.f90:2274
range4
subroutine range4(WWINT, IDDLOW, IDDTOP)
Definition: swancom4.f90:507
mod_time::abs
Definition: mod_time.f90:64
swcomm4::nptst
integer nptst
Definition: swmod1.f90:2272
swcomm3::mdc4mi
integer mdc4mi
Definition: swmod1.f90:1672
funcw1
function funcw1(G, H, X, AKA, T1A, WA)
Definition: swancom4.f90:2911
m_snl4::af11
real, dimension(:), allocatable, save, public af11
Definition: swmod2.f90:370
swcomm3::iquad
integer iquad
Definition: swmod1.f90:2127
m_snl4::cnl4_2
real, dimension(:), allocatable, save, public cnl4_2
Definition: swmod2.f90:372
swcomm3::grav_w
real grav_w
Definition: swmod1.f90:1721
wave
real function wave(D)
Definition: swancom4.f90:2717
vars_wave::ac2
real, dimension(:,:,:), allocatable ac2
Definition: mod_main_wave.f90:24
swcomm3::msc4mi
integer msc4mi
Definition: swmod1.f90:1673
swcomm3::pquad
real, dimension(mquad) pquad
Definition: swmod1.f90:2140
m_snl4::friam
type(riamdat), target, save friam
Definition: swmod2.f90:381
swsnl2
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)
Definition: swancom4.f90:1959
swcomm3::mdc
integer mdc
Definition: swmod1.f90:1672
mod_time::mod
Definition: mod_time.f90:69
cgcmp
subroutine cgcmp(G, AK, H, CG)
Definition: swancom4.f90:2739
m_snl4::curriam
type(riamdat), pointer, save curriam
Definition: swmod2.f90:382
swcomm3::pi_w
real pi_w
Definition: swmod1.f90:1721
fdw1
subroutine fdw1(G, H, AKA, T1A, WA, X1, X2, X, EPS, M)
Definition: swancom4.f90:2881
swcomm3::mdc4ma
integer mdc4ma
Definition: swmod1.f90:1672
riam_slw
subroutine riam_slw(LMAX, N, N2, G, H, DQ, DQ2, DT, DT2, W, P, ACT, SNL, MINT)
Definition: swancom4.f90:896
swsnl1
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)
Definition: swancom4.f90:2254
swlta
subroutine swlta(IG, DEP2, CGO, SPCSIG, KWAVE, IMATRA, IMATDA, IDDLOW, IDDTOP, ISSTOP, IDCMIN, IDCMAX, HS, SMEBRK, PLTRI, URSELL)
Definition: swancom4.f90:287
swsnl3
subroutine swsnl3(WWINT, WWAWG, UE, SA1, SA2, SPCSIG, SNLC1, DAL1, DAL2, DAL3, SFNL, DEP2, KMESPC, MEMNL4, FACHFR)
Definition: swancom4.f90:1558
swcomm3::ptriad
real, dimension(mtriad) ptriad
Definition: swmod1.f90:2135
preset
subroutine preset(LMAX, N, IW3L, IW4, N2, G, H, W4, AK4, W, DQ, DQ2, DT, DT2, P)
Definition: swancom4.f90:2786
m_snl4
Definition: swmod2.f90:292
all_vars
Definition: mod_main.f90:959
filnl3
subroutine filnl3(IDCMIN, IDCMAX, IMATRA, IMATDA, MEMNL4, PLNL4S, ISSTOP, IGC)
Definition: swancom4.f90:597
swcomm3::kcgrd
integer, dimension(micmax) kcgrd
Definition: swmod1.f90:1670
m_snl4::cnl4_1
real, dimension(:), allocatable, save, public cnl4_1
Definition: swmod2.f90:371
swsnl8
subroutine swsnl8(WWINT, UE, SA1, SA2, SPCSIG, SNLC1, DAL1, DAL2, DAL3, SFNL, DEP2, KMESPC, MEMNL4, FACHFR)
Definition: swancom4.f90:665
swcomm4::iptst
integer iptst
Definition: swmod1.f90:2269
findw1w2
subroutine findw1w2(LMAX, IW3, W, G, H, W1, W2, W3, WA, AK1, AK2, AKA, T1A, T2A, IND)
Definition: swancom4.f90:2853
fac4ww
subroutine fac4ww(XIS, SNLC1, DAL1, DAL2, DAL3, SPCSIG, WWINT, WWAWG, WWSWG)
Definition: swancom4.f90:17
swcomm3::msc4ma
integer msc4ma
Definition: swmod1.f90:1673
swsnl4
subroutine swsnl4(WWINT, WWAWG, SPCSIG, SNLC1, DAL1, DAL2, DAL3, DEP2, KMESPC, MEMNL4, FACHFR, IDIA, ITER)
Definition: swancom4.f90:1143
swcomm3::ddir
real ddir
Definition: swmod1.f90:1676
Generated by   doxygen 1.8.15