-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathwwm_mod_xnl4v5.F90
8240 lines (8159 loc) · 265 KB
/
wwm_mod_xnl4v5.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#include "wwm_functions.h"
!------------------------------------------------------------------------------
module m_xnldata
USE DATAPOOL, only : rkind
!------------------------------------------------------------------------------
! module for computing the quadruplet interaction
! Created by Gerbrant van Vledder
!
! version 1.01 16/02/1999 Initial version
! 2.01 01/10/2001 various extensions added
! 3.1.01 01/10/2001 Array's for k4 -locus added
! 3.2 12/05/2002 Triplet data added
! 4.00 08/08/2002 Upgrade to version 4.0
! 4.01 19/08/2002 Various modifications for consistency reasons
! 5.01 9/09/2002 Length of strings aqname and bqname modified
! q_dstep added, step for BQF files
! 11/09/2002 Filtering variables added
! 5.02 12/04/2003 Switch for triplet variables corrected
! 5.03 26/05/2003 Switch for lumping along locus added
! 04/06/2003 Switch for Gauss-Legendre integration added
! 06/06/2003 Switch iq_xdia added and NXDIA removed
! 12/06/2003 Loop indices ik1,ia1,ik3,ia1 added
! 16/06/2003 Switch IQ_SYM introduced
! 04/09/2003 Version string set in subroutine q_version
! 09/09/2003 Parameter id_facmax introduced
! 5.04 24/12/2003 Tail factors for k2 and k4 always in BQF
! 30/12/2003 Parameters IQ_TAIL & FF_TAIL added
!------------------------------------------------------------------------------------
implicit none
!
character(len=60) q_version ! version string
!
character(len=20) sub_name ! Name of active subroutine
character(len=20) qbase ! base name for I/O files
character(len=20) qf_error ! name of file with error messages
!
integer iufind ! Specifies handling of unit numbers, see Z_FILEIO
integer iscreen ! identifier for screen, set in XNL_INIT
!
! unit numbers for I/O
!
integer luq_bqf ! binary file storing and retrieving precomputed loci
integer luq_cfg ! user defined configuration
integer luq_err ! file with error messages
integer luq_fil ! test output for filtering
integer luq_grd ! ASCII file storing and retrieving precomputed loci
integer luq_int ! test file for test output of integration
integer luq_loc ! statistics about computed loci
integer luq_log ! logging
integer luq_prt ! general print file for quadruplets
integer luq_trf ! testing transformation of loci
integer luq_tst ! test file for quadruplets
integer luq_txt ! reading (error) text file
integer luq_t13 ! test of basis integration
!------------------------------------------------------------------------------
! physical coefficients, to be obtained through interface XNL_INIT
!------------------------------------------------------------------------------
real(rkind) q_grav ! gravitational acceleration (Earth = 9.81 m/s^2)
real(rkind) qf_tail ! power of spectral tail of E(f), e.g. -4,, -4.5, -5
! ! these values must be set in the interface routine
!------------------------------------------------------------------------------
! filtering coefficients
!------------------------------------------------------------------------------
real(rkind) qf_krat ! maximum ratio of the interacting wave numbers k1 and k3
real(rkind) qf_dmax ! maximum directional difference between k1 and k3
real(rkind) qf_frac ! fraction of maximum action density to filter
!
! program switches, optionally to be reset in routine Q_SETCONFIG
!
integer iq_compact ! switch to compact data
! == 0, do not compact
! == 1, compact data by elimiting zero contribution along locus
!
integer iq_cple ! type of coupling coefficient
! == 1, deep water coefficient of Webb
! == 2, deep water coefficient of Zakharov
! == 3, finite depth coefficient of Hasselmann & Herterich
! == 4, finite depth coefficient of Zakharov
! == 5, finite depth coefficient of Lin & Perrie
!
integer iq_disp ! type of dispersion relation, viz. depth dependency
! == 1, deep water, possibly with geometric scaling
! == 2, linear dispersion relation, w^2 = g.k.tanh(kd)
! == 3, nonlinear dispersion relation
!
integer iq_dscale ! switch to activate depth scaling according to
! Herterich and Hasselmann
! ! == 0, No depth scaling
! ! == 1, depth scaling activated
!
integer iq_filt ! switch to activate filtering in wave number space
! ! ==0, no filtering
! ! ==1, filtering activated
!
integer iq_gauleg ! switch for Gauss-Legendre interpolation
! ! == 0, No Gauss-Legendre, default
! ! > 0 Gauss-Legendre, iq_gauleg is number of points
!
integer iq_geom ! type of scaling
! == 0, no geometric scaling, only directional scaling of loci
! == 1, geometric scaling using Resio/Tracy method
! only possible in the case IQ_DISP=1
!
integer iq_grid ! type of spectral grid
! == 1, sector & symmetric around zero
! == 2, sector & symmetric around zero & non-symmetric
! == 3, full circle & non-symmetric
!
integer iq_integ ! option to output integration results
! ! ==0 no output of integration
! ! ==1 only sum per locus
! ! ==2 also information per point on locus
! ! ==3 only basic line integrals
!
integer iq_interp ! type of interpolation to retrieve action density
! ! == 1, bi-linear interpolation in discrete spectrum (default)
! ! == 2, take nearest bins, on the basis of maximum weight
!
integer iq_locus ! Option for computation of locus
! ! ==1, explicit polar method with fixed k-step
! ! ==2, explicit polar method with adpative k-stepping
! ! ==3, explicit polar method with geometric k-spacing
!
integer iq_log ! switch to activate logging to file QBASE//.LOG
! ! == 0, No print output
! ! == 1, print output
!
integer iq_lump ! switch to activate lumping on locus
! ! == 0, No lumping
! ! == 1, Lumping along locus
!
integer iq_make ! option to make quadruplet grid
! == 1, make when needed (default)
! == 2, always make quadruplet grid
! == 3, only make grid file
!
integer iq_mod ! option to redistribute points on locus
! ! == 0, Points will be used as computed by tracing algortihm
! ! == 1, Equi-distant spacing on points along locus (NLOC1)
!
integer iq_prt ! switch to activate print output, to file QBASE//.PRT
! ! == 0, No print output
! ! == 1, print output
!
integer iq_search ! switch to determine search for a proper grid
! == 0, no search is carried out
! == 1, search nearest (relative) interaction grid
!
integer iq_screen ! option to send output to the screen
! ! == 0, no output is send to screen
! ! == 1, output is send to screen
!
integer iq_sym ! switch to activate use of symmetry reduction
! ! == 0, no symmetries are used
! ! == 1, symmetry activated (default)
!
integer iq_tail ! add parametric tail to transfer rate and diagnonal term
! ! == 0, no tail is added
! ! == 1, parametric tail is added
!
integer iq_test ! test level, output is directed to unit luqtst
! ! == 0, no test output
! ! == 1, output of basic I/O
! ! == 2, extensive test output
!
integer iq_trace ! trace option
! ! == 0, no trace of subroutine calls
! ! > 0, maximum number of traces per subroutine
! ! < 0, as for >0 but now output is send to the screen
!
integer iq_trf ! option to print transformed loci to special output file
! ! == 0, no output to data file unit luqtrf
! ! == 1, test output from routine Q_GETLOCUS
!
integer iq_t13 ! option to output T13 integration
! ! ==0, no output
! ! ==1, test output of T13 per locus
!
integer iq_xdia ! switch to activate output to extended DIA data file
! == 0, no output
! > 0, output to data file, but only when lumping is also
! activated
!---------------------------------------------------------------------------------------
!
!
! grid administration
!
character(len=13) aqname ! name of ASCII grid file
character(len=13) bqname ! name of binary quadruplet grid file
character(len=13) lastquadfile ! name of last retrieved BQF file
character(len=21) q_header ! header of Binary Quadruplet File as intended in BQF-file
character(len=21) r_header ! header of Binary Quadruplet File as exists in BQF-file
logical lq_grid ! flag to make (new) interaction grid
!
integer nkq ! number of wave numbers of quad-grid
integer naq ! number of angles of quad-grad
integer ncirc ! number of angles on a full circle
!
integer ia_k1,ik_k1 ! indices of main loop variables
integer ia_k3,ik_k3 ! indices of main loop variables
!
real(rkind) fqmin ! lowest frequency in Hz
real(rkind) fqmax ! highest frequency in Hz
real(rkind) q_sector ! half plane width in degrees (for iq_grid=1,2)
real(rkind) q_dstep ! step size for generating BQF files
!
integer, parameter :: mq_stack=10 ! maximum number of elements in stack
!
integer mlocus ! maximum number of points on locus for defining arrays
integer nlocus0 ! preferred number of points on locus
integer nlocus1 ! number of points on locus as computed in Q_CMPLOCUS
integer klocus ! number of points on locus as stored in quadruplet database
! based on nlocus0, iq_gauleg and iq_lump (without compacting)
! used in Q_ALLOCATE to define size of data arrays
integer nlocus ! number of points on locus, equal to klocus
integer nlocusx ! number of points on locus for use in computation (nlocusx <= nlocus)
!
real(rkind) kqmin ! lowest wave number
real(rkind) kqmax ! highest wave number
real(rkind) wk_max ! maximum weight for wave number interpolation, set in Q_INIT
!
real(rkind) k0x,k0y,dk0 ! components of initial wave number of locus,
real(rkind) krefx,krefy ! components of reference wave number for quad-grid
real(rkind) k1x,k1y ! components of k1 wave number
real(rkind) k2x,k2y ! components of k2 wave number
real(rkind) k3x,k3y ! components of k3 wave number
real(rkind) k4x,k4y ! components of k4 wave number
real(rkind) px,py ! components of difference k1-k3 wave number
real(rkind) pmag ! magnitude of P-vector
real(rkind) pang ! angle related of P-vector, Pang = atan2(py,px), (radians)
real(rkind) sang ! angle of symmytry axis of locus, SANG = PANG +/ pi° (radians)
real(rkind) xang ! angle of locus for the case that w1=w3, Xang=atan2(-px,py), (radians)
real(rkind) q ! difference of radian frequencies, used in Resio-Tracy method
real(rkind) kmin_loc ! minimum wave number of locus along symmetry axis
real(rkind) kmax_loc ! maximum wave number of locus along symmetry axis
real(rkind) kmid ! wave number at midpoint of locus along symmetry axis
real(rkind) kmidx ! x-component of wave number at midpoint of locus along symmetry axis
real(rkind) kmidy ! y-component of wave number at midpoint of locus along symmetry axis
real(rkind) loc_crf ! circumference of locus in (kx,ky)-space
real(rkind) loc_area ! area of locus, measured in (kx-ky)- space
real(rkind) loc_xz ! x-coordinate of center of gravity of locus in (kx,ky)-space
real(rkind) loc_yz ! y-coordinate of center of gravity of locus in (kx,ky)-space
!
! data for extended input k-grid, necessary when input grid is smaller than
! internal k-grid.
!
! real(rkind) fackx ! geometric spacing factor of input grid
! integer nkx ! new number of k-rings of extended input grid
! real(rkind), allocatable :: kx(:) ! extended k-grid
! real(rkind), allocatable :: nspecx(:,:) ! extended action density spectrum
!
! information about pre_computed locus, only half the angles need to be saved
!
!
integer, allocatable :: quad_nloc(:,:) ! number of points on locus
integer, allocatable :: quad_ik2(:,:,:) ! lower wave number index of k2
integer, allocatable :: quad_ia2(:,:,:) ! lower direction index of k2
integer, allocatable :: quad_ik4(:,:,:) ! lower wave number index of k4
integer, allocatable :: quad_ia4(:,:,:) ! lower direction index of k4
real(rkind), allocatable :: quad_w1k2(:,:,:) ! weight 1 of k2
real(rkind), allocatable :: quad_w2k2(:,:,:) ! weight 2 of k2
real(rkind), allocatable :: quad_w3k2(:,:,:) ! weight 3 of k2
real(rkind), allocatable :: quad_w4k2(:,:,:) ! weight 4 of k2
real(rkind), allocatable :: quad_w1k4(:,:,:) ! weight 1 of k4
real(rkind), allocatable :: quad_w2k4(:,:,:) ! weight 2 of k4
real(rkind), allocatable :: quad_w3k4(:,:,:) ! weight 3 of k4
real(rkind), allocatable :: quad_w4k4(:,:,:) ! weight 4 of k4
real(rkind), allocatable :: quad_zz (:,:,:) ! compound product of cple*ds*sym/jac
real(rkind), allocatable :: quad_t2(:,:,:) ! tail factors for k2 wave number
real(rkind), allocatable :: quad_t4(:,:,:) ! tail factors for k4 wave number
real(rkind), allocatable :: quad_sym(:,:,:) ! symmetry factor btween k1-k3 and k1-k4
real(rkind), allocatable :: quad_jac(:,:,:) ! jacobian term
real(rkind), allocatable :: quad_cple(:,:,:) ! coupling coefficient
real(rkind), allocatable :: quad_ws (:,:,:) ! (weighted) step size
!
! characteristic of computed locus
!
real(rkind), allocatable :: x2_loc(:) ! k2x coordinates around locus
real(rkind), allocatable :: y2_loc(:) ! k2y coordinates around locus
real(rkind), allocatable :: z_loc(:) ! data value around locus
real(rkind), allocatable :: s_loc(:) ! coordinate along locus
real(rkind), allocatable :: x4_loc(:) ! k4x coordinates around locus
real(rkind), allocatable :: y4_loc(:) ! k4y coordinates around locus
real(rkind), allocatable :: ds_loc(:) ! step size around locus
real(rkind), allocatable :: jac_loc(:) ! jacobian term around locus
real(rkind), allocatable :: cple_loc(:) ! coupling coefficient around locus
real(rkind), allocatable :: sym_loc(:) ! factor for symmetry between k3 and k4
!
real(rkind), allocatable :: k_pol(:) ! wave numbers during polar generation of locus
real(rkind), allocatable :: c_pol(:) ! cosines during polar generation of locus
real(rkind), allocatable :: a_pol(:) ! angles of polar locus
!
! characteristics of modified locus, result
!
real(rkind), allocatable :: x2_mod(:) ! k2x coordinates along locus
real(rkind), allocatable :: y2_mod(:) ! k2y coordinates along locus
real(rkind), allocatable :: x4_mod(:) ! k4x coordinates along locus
real(rkind), allocatable :: y4_mod(:) ! k4y coordinates along locus
real(rkind), allocatable :: z_mod(:) ! data value around locus
real(rkind), allocatable :: s_mod(:) ! coordinate along locus
real(rkind), allocatable :: ds_mod(:) ! step size around locus
real(rkind), allocatable :: jac_mod(:) ! jacobian term around locus
real(rkind), allocatable :: cple_mod(:) ! coupling coefficient around locus
real(rkind), allocatable :: sym_mod(:) ! factor for symmetry between k3 and k4
!
real(rkind), allocatable :: k2m_mod(:) ! k2 magnitude around locus
real(rkind), allocatable :: k2a_mod(:) ! k2 angle around locus
real(rkind), allocatable :: k4m_mod(:) ! k4 magnitude around locus
real(rkind), allocatable :: k4a_mod(:) ! k4 angle around locus
!
! result of subroutine Q_weight
!
real(rkind), allocatable :: wk_k2(:) ! position of k2 and k4 wave number
real(rkind), allocatable :: wk_k4(:) ! w.r.t. discrete k-grid
real(rkind), allocatable :: wa_k2(:) ! position of k2 and k4 wave number
real(rkind), allocatable :: wa_k4(:) ! w.r.t. discrete a-grid
real(rkind), allocatable :: wt_k2(:) ! weight factor in tail,
real(rkind), allocatable :: wt_k4(:) ! wt==1 for wave numbers inside k-grid
!
integer, allocatable :: t_ik2(:) ! transformed weight for k2-magnitude
integer, allocatable :: t_ia2(:) ! transformed direction for k2
integer, allocatable :: t_ik4(:) ! transformed tail factor for k2
integer, allocatable :: t_ia4(:) ! transformed weight for k4
real(rkind), allocatable :: t_w1k2(:) ! transformed weight 1 for k2
real(rkind), allocatable :: t_w2k2(:) ! transformed weight 2 for k2
real(rkind), allocatable :: t_w3k2(:) ! transformed weight 3 for k2
real(rkind), allocatable :: t_w4k2(:) ! transformed weight 4 for k2
real(rkind), allocatable :: t_w1k4(:) ! transformed weight 1 for k4
real(rkind), allocatable :: t_w2k4(:) ! transformed weight 2 for k4
real(rkind), allocatable :: t_w3k4(:) ! transformed weight 3 for k4
real(rkind), allocatable :: t_w4k4(:) ! transformed weight 4 for k4
real(rkind), allocatable :: t_zz(:) ! product term
real(rkind), allocatable :: t_tail2(:) ! tail factor for k2
real(rkind), allocatable :: t_tail4(:) ! tail factor for k4
real(rkind), allocatable :: t_sym(:) ! transformed symetry factor
real(rkind), allocatable :: t_cple(:) ! transformed coupling coefficient
real(rkind), allocatable :: t_jac(:) ! tranformed jacobian term
real(rkind), allocatable :: t_ws(:) ! transformed weighted/step size
!
! corresponding declarations
!
integer, allocatable :: r_ik2(:)
integer, allocatable :: r_ia2(:)
integer, allocatable :: r_ik4(:)
integer, allocatable :: r_ia4(:)
real(rkind), allocatable :: r_w1k2(:),r_w2k2(:),r_w3k2(:),r_w4k2(:)
real(rkind), allocatable :: r_w1k4(:),r_w2k4(:),r_w3k4(:),r_w4k4(:)
real(rkind), allocatable :: r_zz(:),r_jac(:),r_cple(:),r_sym(:),r_ws(:)
real(rkind), allocatable :: r_tail2(:),r_tail4(:)
!
real(rkind), allocatable :: dt13(:) ! increment along locus
!
real(rkind), allocatable :: q_xk(:) ! extended wave number array starting at index 0
real(rkind), allocatable :: q_sk(:) ! step size of extended wave number array
real(rkind) sk_max ! maximum wave number in extended array
!
real(rkind), allocatable :: q_k(:) ! wave number grid [1/m]
real(rkind), allocatable :: q_dk(:) ! width of wave number bins [1/m]
real(rkind), allocatable :: q_kpow(:) ! wave number to a certain power, used in filtering
real(rkind), allocatable :: q_f(:) ! frequencies accociated to wave number/depth
real(rkind), allocatable :: q_df(:) ! step size of frequency grid
real(rkind), allocatable :: q_sig(:) ! radian frequencies associated to wave number/depth
real(rkind), allocatable :: q_dsig(:) ! step size of radian frequency grid
real(rkind), allocatable :: q_cg(:) ! group velocity (m/s)
real(rkind), allocatable :: q_a(:) ! directions of quadruplet grid in radians
real(rkind), allocatable :: q_ad(:) ! directions of quadruplet grid in degrees
real(rkind), allocatable :: a(:,:) ! Action density on wave number grid A(sigma,theta)
real(rkind), allocatable :: nspec(:,:) ! Action density on wave number grid N(kx,ky)
real(rkind), allocatable :: nk1d(:) ! Internal 1d action density spectrum N(k)
real(rkind), allocatable :: qnl(:,:) ! Nonlinear energy transfer Snl(k,theta)
!
integer id_facmax ! Factor for determining range of depth search (Q_SEARCHGRID)
real(rkind) q_dird1,q_dird2 ! first and last direction of host model (via XNL_INIT) degrees
real(rkind) q_depth ! local water depth in m
real(rkind) q_maxdepth ! maximum water depth, set in XNL_INIT, used in Q_CTRGRID
real(rkind) q_mindepth ! minimum water depth, set in XNL_INIT, used in Q_CTRGRID
real(rkind) q_lambda ! geometric scaling factor for 'deep' water loci
real(rkind) q_scale ! additional scale factor resulting from SEARCH for neasrest grid
!
real(rkind) eps_q ! absolute accuracy for check of Q
real(rkind) eps_k ! absolute accuracy for equality check of k
real(rkind) rel_k ! relative accuracy for equality check of k
!
integer iq_stack ! Sequence number of stack with subroutine calls
character(len=21) cstack(mq_stack) ! Stack with module names
!
! characteristics of locus
!
real(rkind) crf1 ! estimated circumference of locus
!---------------------------------------------------------------------------------
!
! information about type of grid
!
integer iaref ! index of first angle of reference wave numbers
integer iamax ! maximum difference in indices for sector grids
integer iaq1,iaq2 ! indices of do-loop for directions
integer iag1,iag2 ! range of directions for precomputed interaction grid
real(rkind) q_ang1,q_ang2 ! lower and upper angle of grid in degrees
real(rkind) q_delta ! directional spacing of angular grid in radians
real(rkind) q_deltad ! directional spacing of angular grid in degrees
!
real(rkind) q_ffac ! geometric factor between subsequent frequencies
real(rkind) q_kfac ! geometric factor between subsequent wave numbers
! (only valid for IQ_IDISP==1)
real(rkind) qk_tail ! power of spectral tail of N(k), computed from qf_tail
real(rkind) ff_tail ! fraction of maximum frequency where parametric tail starts
!
!-----------------------------------------------------------------------------
!
!!/R real(rkind) wq2(4) ! interpolation weights for k2
!!/R real(rkind) wq4(4) ! interpolation weights for k4
!!/R real(rkind) wqw ! overall weight of contribution
!!/R real(rkind) wtriq(40) ! triplet weights
!!/R integer ikq2(4) ! wave number index for k2
!!/R integer idq2(4) ! angle index for k2
!!/R integer ikq4(4) ! wave number index for k4
!!/R integer idq4(4) ! angle index for k4
!!/R integer iktriq(40,3) ! k-indices of triplets
!!/R integer idtriq(40,3) ! direction indices of triplets
!
!============== General settings =================
!
integer iq_type ! method for computing the nonlinear interactions
! depending on the value of iq_type a number of settings
! for other processes or schematizations are set in Q_COMPU
! iq_type==1: deep water, symmetric spectrum, Webb coupling coefficient
! 2: deep water computation with WAM depth scaling based on Herterich
! and Hasselmann (1980)
! 3: finite depth transfer
!
integer iq_err ! counts the number of errors
! if no error occurred, IQ_ERR = 0
! for each occuring error, iq_err is incremented
! errors are always terminating
! routine Q_ERROR handles the reporting on the error
!
integer iq_warn ! counts the number of warnings
!
! indices for test output of actual integration
! these values are set and optionally modified in Q_SETCONFIG
!
integer mk1a,mk1b ! indices of k1 when test output is needed
integer mk3a,mk3b ! indices of k3 when test output is needed
contains
!----------------------------------------------------------------------------------
!------------------------------------------------------------------------------
subroutine xnl_init(sigma,dird,nsigma,ndir,pftail,x_grav,depth,ndepth, &
& iquad,iqgrid,iproc,ierror)
!------------------------------------------------------------------------------
!
! +-------+ ALKYON Hydraulic Consultancy & Research
! | | Gerbrant van Vledder
! | +---+
! | | +---+ Last update: 6 May 2004
! +---+ | | Release: 5.03
! +---+
!
use datapool,only:stat
use m_fileio
use m_constants
! do not use m_xnldata
!
implicit none
!
! 0. Update history
!
! 10/01/2001 Initial version
! 14/02/2001 Release 3
! 06/11/2001 Depth forced to 1000 when iquad ==1,2
! 08/08/2002 Upgrade to version 4
! 19/08/2002 Extra test output
! 22/08/2002 User defined directions stored in Quad-system
! 26/08/2002 Minimum water depth in variable q_mindepth
! 09/09/2002 Initialized LASTQUADFILE
! 10/09/2002 Initialisation of Q_DSTEP
! 11/09/2002 Called of Q_ALLOC moved to location after Q_SETCFG
! Output unit luq_fil added
! 16/09/2002 Parameter IPROC added to take are of MPI
! 25/09/2002 Check added for directions of sector grid
! 25/04/2003 name q_alloc changed to q_allocate
! 04/06/2003 variable IQ_INT renamed IQ_INTEG
! 11/06/2003 Call to Q_SETCFG changed into Q_SETCONFIG
! Call to Q_CHKCFG changed into Q_CHKCONFIG
! Call to subroutine Q_SUMMARY added
! Compute size of points on locus, stored in KLOCUS
! 13/06/2003 Test parameters moved to Q_SETCONFIG
! 04/09/2003 Routine Q_SETVERSION added
! 27/04/2004 Check added for directions of sector grid
! 06/05/2004 Initialisation of IQ_INTERP removed
!
! 1. Purpose:
!
! Initialize coefficients, integration space, file i/o for computation
! nonlinear quadruplet wave-wave interaction
!
! 2. Method
!
! Set version number
! Set unit unit numbers
! Open quad related files
! Optionally reset configuration by a back door option
! Compute integration spaces for given water depths
!
! 3. Parameter list:
!
!Type I/O Name Description
!------------------------------------------------------------------------------
integer, intent(in) :: nsigma ! Number of sigma values
integer, intent(in) :: ndir ! Number of directions
integer, intent(in) :: ndepth ! Number of water depths
real(rkind), intent(in) :: sigma(nsigma) ! Radian frequencies
real(rkind), intent(in) :: dird(ndir) ! Directions (degrees)
real(rkind), intent(in) :: pftail ! power of spectral tail, e.g. -4 or -5
real(rkind), intent(in) :: depth(ndepth) ! depths for which integration space must be computed
real(rkind), intent(in) :: x_grav ! gravitational acceleration
integer, intent(in) :: iquad ! Type of method for computing nonlinear interactions
integer, intent(in) :: iqgrid ! Type of grid for computing nonlinear interactions
integer, intent(in) :: iproc ! Processor number, controls output file for MPI
integer, intent(out) :: ierror ! Error indicator. If no errors are detected IERR=0
!
! 4. Error messages
!
! An error message is produced within the QUAD system.
! If no errors are detected IERROR=0
! otherwise IERROR > 0
!
! ierr Description of error
! -------------------------------
! 1 Invalid value of iquad
! 2 Invalid value of iq_grid
! 31 Incompatability between iq_grid and input directions for circle grid
! 32 Incompatability between iq_grid and input directions for sector grid grid
! 4 Error in deleting *.ERR file
! 5 Error generated by Q_SETCONFIG
! 6 Error generated by Q_CHKCFG
! 7 Error generated by Q_CTRGRID
!
! 5. Called by:
!
! host program, e.g. SWANQUAD5
!
! 6. Subroutines used:
!
! Q_SETVERSION
! Q_SETCONFIG
! Q_CHKCFG
! Q_SUMMARY
! Q_ALLOCATE
! Q_CTRGRID
! Q_INIT
!
! 7. Remarks
!
! 8. Structure
!
! 9. Switches
!
! 10. Source code
!---------------------------------------------------------------------------------------
!
! Local parameters
!
integer iuerr ! error indicator
integer idepth ! index over water depths
integer igrid ! status of quadruplet grid file
integer ia,ik ! counters
!
real(rkind) depmin ! minimum water depth
real(rkind) dstep ! directional step
real(rkind) dgap ! directional gap between first and last direction
!
call q_setversion ! set version number
!------------------------------------------------------------------------------
! user defined settings
!------------------------------------------------------------------------------
q_mindepth = 0.1 ! Set minimum water depth
q_maxdepth = 2000 ! Set maximum water depth
q_dstep = 0.1 ! Set minimum step for coding depth files
iscreen = STAT%FHNDL ! Put into STAT&FHNDL was screen before makes problems in parallel
iufind = 1 ! search for unit numbers automatically
!----------------------------------------------------------------------------
! Initialisations
!
ierror = 0 ! set error condition
iq_stack = 0 ! initialize stack for tracing subroutines
qbase = 'xnl4v5' ! Base name for quadruplet files
qf_error = 'xnl5_errors.txt' ! Text file with error messages
lastquadfile = 'quad?????.bqf' ! Initialize name of last retrieved quad file
!
! set values of physical quantities
! and store them in quad data area
!
q_grav = x_grav ! gravitational acceleration
qf_tail = pftail ! Power of parametric spectral tail
iq_type = iquad ! Type of method to compute transfer
iq_prt = 0 ! Print output on, to file /qbase/.prt
iq_test = 0 ! test level
iq_trace = 0 ! level of subroutine trace
iq_log = 0 ! Set logging of q_routines off
iq_grid = iqgrid ! Grid type for computation of nonlinear transfer
iq_screen = 0 ! enable output to screen
!
!------------------------------------------------------------------------------
! Check input
!------------------------------------------------------------------------------
if(iq_type<1 .or. iq_type>3) then
ierror = 1
goto 9999
end if
!
if(iq_grid<1 .or. iq_grid>3) then
ierror = 2
goto 9999
end if
!
! Retrieve size of spectral grid from input
!
fqmin = sigma(1)/(2.*pi) ! minimum frequency in Hz
fqmax = sigma(nsigma)/(2.*pi) ! maximum frequency in Hz
nkq = nsigma ! number of frequencies/wave numbers
naq = ndir ! number of directions
!
! check if directions are given on full circle or in a symmetric sector
!----------------------------------------------------------------------------
! 1: compute directional step
! 2: compute gap between first and last
! 3: compare gap with step
!
dstep = dird(2)-dird(1) ! directional step
dgap = 180.- abs(180.- abs(dird(1)-dird(ndir))) ! directional gap
!
!----------------------------------------------------------------------
if(iq_grid==1 .or. iq_grid==2) then
!
! check if gap equal to step in the case of full circle
!
if(abs(dstep-dgap) < 0.001) then
ierror = 31
write(iscreen,'(a,i4,2f10.2)') 'XNL_INIT iq_grid dstep dgap:',iq_grid,dstep,dgap
goto 9999
end if
!
! check if sector is symmetric around zero in the case of sector grid
!
if(abs(dird(1)+dird(ndir)) > 0.01) then
write(iscreen,'(a,i4,2f10.2)') 'XNL_INIT iq_grid dird(1) dird(n):',iq_grid,dird(1),dird(ndir)
ierror = 32
goto 9999
end if
end if
!
q_dird1 = dird(1)
q_dird2 = dird(ndir)
!
! assign unit numbers for QUAD related files
!
! If IUFIND=0, fixed prespecified unit numbers must be given
! IUFIND=1, the numbers are searched automatically
!
if(iufind==0) then
luq_err = 103
luq_tst = 104
luq_int = 105
luq_log = 106
luq_prt = 107
luq_cfg = 108
luq_bqf = 109
luq_grd = 110
luq_txt = 111
luq_loc = 112
luq_trf = 113
luq_t13 = 114
luq_fil = 117
end if
!
! delete old Error file, if it exists
!
call z_fileio(trim(qbase)//'.err','DF',iufind,luq_err,iuerr)
if(iuerr/=0) then
call q_error('e','FILEIO','Problem in deleting error file *.ERR')
ierror = 4
goto 9999
end if
!
! create new files, first create logging file
!
call z_fileio(trim(qbase)//'.log','UF',iufind,luq_log,iuerr) ! logging
call z_fileio(trim(qbase)//'.prt','UF',iufind,luq_prt,iuerr) ! general print file
call z_fileio(trim(qbase)//'.tst','UF',iufind,luq_tst,iuerr) ! test output
!
call z_fileio(trim(qbase)//'.int','UF',iufind,luq_int,iuerr) ! logging for locus integration
call z_fileio(trim(qbase)//'.trf','UF',iufind,luq_trf,iuerr) ! transformation test
call z_fileio(trim(qbase)//'.t13','UF',iufind,luq_t13,iuerr) ! test output for integration of T13
!
write(luq_log,'(2a,i4)') 'XNL_INIT: ',trim(qbase)//'.log connected to :',luq_log
write(luq_log,'(2a,i4)') 'XNL_INIT: ',trim(qbase)//'.prt connected to :',luq_prt
write(luq_log,'(2a,i4)') 'XNL_INIT: ',trim(qbase)//'.tst connected to :',luq_tst
!
write(luq_log,'(2a,i4)') 'XNL_INIT: ',trim(qbase)//'.int connected to :',luq_int
write(luq_log,'(2a,i4)') 'XNL_INIT: ',trim(qbase)//'.trf connected to :',luq_trf
write(luq_log,'(2a,i4)') 'XNL_INIT: ',trim(qbase)//'.t13 connected to :',luq_t13
!
write(luq_prt,'(a)') '---------------------------------------------------------------'
write(luq_prt,'(a)') trim(q_version)
write(luq_prt,'(a)') 'Solution of Boltzmann integral using Webb/Resio/Tracy method'
write(luq_prt,'(a)') '---------------------------------------------------------------'
write(luq_prt,*)
write(luq_prt,'(a)') 'Initialisation'
write(luq_prt,*)
!
if(iproc >=0) write(luq_prt,'(a,i5)') '(MPI) processor number:',iproc
!---------------------------------------------------------------------------------
! Reset configuration from file, using a backdoor
!---------------------------------------------------------------------------------
call q_setconfig(iquad)
if (iq_err /=0) then
ierror = 5
goto 9999
end if
!---------------------------------------------------------------------------------
! check settings for inconsistencies
!---------------------------------------------------------------------------------
call q_chkconfig
if (iq_err /=0) then
ierror = 6
goto 9999
end if
!---------------------------------------------------------------------------------
! determine minimum size of number of points on locus as stored in database
!---------------------------------------------------------------------------------
klocus = nlocus0
if(iq_gauleg > 0) klocus = min(iq_gauleg,klocus)
if(iq_lump > 0) klocus = min(iq_lump,klocus)
!---------------------------------------------------------------------------------
! write summary of program settings
!---------------------------------------------------------------------------------
call q_summary
!----------------------------------------------------------------------------------
! allocate data arrays
!-----------------------------------------------------------------------------
call q_allocate
!------------------------------------------------------------------------------
! Generate interaction grid and coefficients for each valid water depth
! Q_CTRGRID controls grid generation
!------------------------------------------------------------------------------
do idepth=1,ndepth
q_depth = depth(idepth)
!
if(iquad==1 .or. iquad==2) q_depth = q_maxdepth
!
if(q_depth < q_mindepth) then
call q_error('w','DEPTH','Invalid depth')
write(luq_err,'(a,e12.5,f10.2)') 'Incorrect depth & minimum:',q_depth,q_mindepth
else
call q_init
call q_ctrgrid(2,igrid)
if(iq_err /= 0) then
ierror = 7
goto 9999
end if
end if
!
if(iquad==1 .and. ndepth > 0) then
write(luq_prt,'(a)') 'XNL_INIT: For deep water only one grid suffices'
exit
end if
end do
!
!
! Create or open triplet output data file if iq_triq > 0
!
!
9999 continue
!
!! if (iq_log ==0) call z_fileio(trim(qbase)//'.log','DF',iufind,luq_log,iuerr)
!! if (iq_prt ==0) call z_fileio(trim(qbase)//'.prt','DF',iufind,luq_prt,iuerr)
!
if (iq_integ==0) call z_fileio(trim(qbase)//'.int','DF',iufind,luq_int,iuerr)
if (iq_test ==0) call z_fileio(trim(qbase)//'.tst','DF',iufind,luq_tst,iuerr)
if (iq_trf ==0) call z_fileio(trim(qbase)//'.trf','DF',iufind,luq_trf,iuerr)
if (iq_t13 ==0) call z_fileio(trim(qbase)//'.t13','DF',iufind,luq_t13,iuerr)
!
return
end subroutine
!-----------------------------------------------------------------------------!
subroutine xnl_main(aspec,sigma,angle,nsig,ndir,depth,iquad,xnl,diag, &
& iproc, ierror)
!-----------------------------------------------------------------------------!
!
! +-------+ ALKYON Hydraulic Consultancy & Research
! | | Gerbrant Ph. van Vledder
! | +---+
! | | +---+ Last update: 7 May 2004
! +---+ | | Release: 5.04
! +---+
!
! do not use m_xnldata
use serv_xnl4v5
!
implicit none
!
! 0. Update history
!
! 25/02/1999 Initial version
! 25/07/1999 Various restructurations
! 12/10/1999 Error handling improved
! 15/10/1999 Existing test file deleted
! 25/10/1999 Parameter iq_filt added
! 29/10/1999 iq_call renamed to i_qmain and save statement added
! 08/12/1999 Call to CHKLAW included
! 16/12/1999 Extra output and check for wave number range
! 27/12/1999 Expansion of k-grid now in new subroutine Q_EXPAND
! 06/01/2000 Deallocate of KX and NSPECX removed
! 08/01/2000 Recontructed, subroutine Q_WRTVV splitted off
! 12/01/2000 Diagonal term added to interface
! 26/06/2000 Name changed to XNL_MAIN
! 01/02/2001 Interface change, spectrum must be given as Action density(sigma,theta)
! 06/11/2001 Water depth forced to 1000 when IQUAD==1,2
! 13/08/2002 Upgrade to release 4.0
! 20/08/2002 Action density array copied to internal A-array
! 09/09/2002 Upgrade to release 5
! 16/09/2002 Parameter IPROC included to take care of MPI processors
! 27/09/2002 Description of input argument SIGMA corrected
! 30/12/2003 Call to subroutine Q_ADDTAIL added
! 07/05/2004 Depth set to Q_MAXDEPTH for iquad=1,2
!
! 1. Purpose:
!
! Compute nonlinear transfer for a given action density spectrum
! on a given sigma and direction grid
!
! 2. Method
!
! Webb/Resio/Tracy/Van Vledder
!
! 3. Parameter list:
!
! Type I/O Name Description
!---------------------------------------------------------------------------------------------
integer,intent(in) :: nsig ! number of frequencies (sigma)
integer,intent(in) :: ndir ! number of directions
integer,intent(in) :: iquad ! method of computing nonlinear quadruplet interactions
integer, intent(in) :: iproc ! MPI processor number
!
real(rkind), intent(in) :: aspec(nsig,ndir) ! Action density spectrum as a function of (sigma,theta)
real(rkind), intent(in) :: sigma(nsig) ! radian frequencies
real(rkind), intent(in) :: angle(ndir) ! directions in radians (sector or full circle)
real(rkind), intent(in) :: depth ! water depth in m
real(rkind), intent(out) :: xnl(nsig,ndir) ! nonlinear quadruplet interaction computed with
! a certain exact method (k,theta)
real(rkind), intent(out) :: diag(nsig,ndir) ! diagonal term for semi-implicit integration
integer, intent(out) :: ierror ! error indicator
!
!--------------------------------------------------------------------------------
!
! 4. Error messages
!
! 5. Called by:
!
! host program
!
! 6. Subroutines used
!
! Q_DSCALE WAM depth scaling
! Q_XNL4V4 Xnl using Webb/Resio/Tracy/VanVledder
! Q_STACK Stack administration
! Q_CHKCONS Check conservation of energy, action and momentum
!
! 7. Remarks
!
! 8. Structure
!
! 9. Switches
!
! 10. Source code:
!--------------------------------------------------------------------------------
! local variables
!
integer, save :: i_qmain ! counter number of calls of XNL_MAIN
integer i_qlast ! value of iquad in last call
!
integer isig ! counter for sigma values
integer idir ! counter of directions
real(rkind) q_dfac ! depth scale factor for nonlinear transfer
!
real(rkind) sum_e ! sum of energy
real(rkind) sum_a ! sum of action
real(rkind) sum_mx ! sum of momentum in x-direction
real(rkind) sum_my ! sum of momentum in y-direction
!
data i_qmain /0/ ! keep track of number of calls of XNL_MAIN
data i_qlast /0/ ! keep track of last call with IQUAD
!
!--------------------------------------------------------------------------------
!
iq_stack =0 ! initialize stack order every time qmain is called
!
call q_stack('+xnl_main')
!
i_qmain = i_qmain + 1
!
if(iq_prt>=1) then
write(luq_prt,*)
write(luq_prt,'(a,i4,f16.3,i4)') 'XNL_MAIN: Input arguments: iquad depth iproc:',&
& iquad,depth,iproc
end if
!
! initialisations for error handling
!
iq_err = 0 ! No errors detected at start
q_depth = depth ! water depth to be used in computation
!
if(iquad==1 .or. iquad==2) q_depth=q_maxdepth
! !
! check water depth to be used in computation
!
if(q_depth < q_mindepth) then
xnl = 0.
call q_error('w','DEPTH','Zero transfer returned')
goto 9999
end if
!
! check if iquad has changed since last call, this is no more allowed
!
!!if (iquad /= i_qlast .and. i_qmain/=1) then
!! call q_error('e','IQUAD','Value of IQUAD differs from initial value')
!! ierror = 1
!! goto 9999
!!end if
!-----------------------------------------------------------------------------+
! main choice between various options |
!-----------------------------------------------------------------------------+
!
if(iquad>=1 .and. iquad <=3) then
!
a = aspec
call q_xnl4v4(aspec,sigma,angle,nsig,ndir,depth,xnl,diag,ierror)
!
if(ierror/=0) then
call q_error('e','wrtvv','Problem in Q_XNL4V4')
goto 9999
end if
!------------------------------------------------------------------------------
! add parametric tail to transfer rate and diagonal term
!------------------------------------------------------------------------------
!
if(iq_tail==1) call q_addtail(xnl,diag,nsig,ndir,qf_tail)
!
!------------------------------------------------------------------------------
! compute scale factor to include WAM depth scaling
!------------------------------------------------------------------------------
!
if(iq_dscale ==1) then
call q_dscale(aspec,sigma,angle,nsig,ndir,depth,q_grav,q_dfac)
!
xnl = xnl*q_dfac
!
if(iq_prt >=1) write(luq_prt,'(a,f7.4)') 'XNL_MAIN depth scale factor:',q_dfac
end if
end if
!
! check conservation laws
!
call q_chkcons(xnl,nsig,ndir,sum_e,sum_a,sum_mx,sum_my)
!
if(iq_prt >= 1) then
write(luq_prt,'(a)') 'XNL_MAIN: Conservation checks'
write(luq_prt,'(a,4e13.5)') 'XNL_MAIN: E/A/MOMX/MOMY:',sum_e,sum_a,sum_mx,sum_my
end if
!
9999 continue
!
ierror = iq_err
!
if(iq_log >= 1) then
write(luq_log,*)
write(luq_log,'(a,i4)') 'XNL_MAIN: Number of warnings:',iq_warn
write(luq_log,'(a,i4)') 'XNL_MAIN: Number of errors :',iq_err
end if
!
!!i_qlast = iquad
!
call q_stack('-xnl_main')
!
return
end subroutine
!------------------------------------------------------------------------------
subroutine q_addtail(xnl,diag,nsig,na,pf_tail)
!------------------------------------------------------------------------------
!
! +-------+ ALKYON Hydraulic Consultancy & Research
! | | Gerbrant van Vledder
! | +---+
! | | +---+ Last update: 27 April 2004
! +---+ | | Release: 5.04
! +---+
!
! do not use m_xnldata
implicit none