forked from NOAA-EMC/WW3
-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy pathw3srcemd.F90
2669 lines (2580 loc) · 87.8 KB
/
w3srcemd.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
!> @file
!> @brief Source term integration routine.
!>
!> @author H. L. Tolman
!> @author F. Ardhuin
!> @date 22-Mar-2021
!>
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
!>
!> @brief Source term integration routine.
!>
!> @author H. L. Tolman
!> @author F. Ardhuin
!> @date 22-Mar-2021
!>
!> @copyright Copyright 2009-2022 National Weather Service (NWS),
!> National Oceanic and Atmospheric Administration. All rights
!> reserved. WAVEWATCH III is a trademark of the NWS.
!> No unauthorized use without permission.
!>
MODULE W3SRCEMD
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | F. Ardhuin |
!/ | FORTRAN 90 |
!/ | Last update : 22-Mar-2021 |
!/ +-----------------------------------+
!/
!/ For updates see subroutine.
!/
! 1. Purpose :
!
! Source term integration routine.
!
! 2. Variables and types :
!
! Name Type Scope Description
! ----------------------------------------------------------------
! OFFSET R.P. Private Offset in time integration scheme.
! 0.5 in original WAM, now 1.0
! ----------------------------------------------------------------
!
! 3. Subroutines and functions :
!
! Name Type Scope Description
! ----------------------------------------------------------------
! W3SRCE Subr. Public Calculate and integrate source terms.
! ----------------------------------------------------------------
!
! 4. Subroutines and functions used :
!
! See corresponding documentation of W3SRCE.
!
! 5. Remarks :
!
! 6. Switches :
!
! See section 9 of W3SRCE.
!
! 7. Source code :
!
!/ ------------------------------------------------------------------- /
!/
REAL, PARAMETER, PRIVATE:: OFFSET = 1.
!/
CONTAINS
!/ ------------------------------------------------------------------- /
!>
!> @brief Calculate and integrate source terms for a single grid point.
!>
!> @verbatim
!> Physics : see manual and corresponding subroutines.
!>
!> Numerics :
!>
!> Dynamic-implicit integration of the source terms based on
!> WW-II (Tolman 1992). The dynamic time step is calculated
!> given a maximum allowed change of spectral densities for
!> frequencies / wavenumbers below the usual cut-off.
!> The maximum change is given by the minimum of a parametric
!> and a relative change. The parametric change relates to a
!> PM type equilibrium range
!>
!> -1 (2pi)**4 1
!> dN(k) = Xp alpha pi ---------- ------------
!> max g**2 k**3 sigma
!>
!> 1 .
!> = FACP ------------ (1)
!> k**3 sigma .
!>
!> where
!> alpha = 0.62e-4 (set in W3GRID)
!> Xp fraction of PM shape (read in W3GRID)
!> FACP combined factor (set in W3GRID)
!>
!> The maximum relative change is given as
!>
!> / +- -+ \ .
!> dN(k) = Xr max | N(k) , max | Nx , Xfilt N(k) | | (2)
!> max \ +- max-+ / .
!>
!> where
!> Xr fraction of relative change (read in W3GRID)
!> Xfilt filter level (read in W3GRID)
!> Nx Maximum parametric change (1)
!> for largest wavenumber.
!> @endverbatim
!>
!> @param[in] srce_call
!> @param[in] IT
!> @param[in] ISEA
!> @param[in] JSEA
!> @param[in] IX Discrete grid point counters.
!> @param[in] IY Discrete grid point counters.
!> @param[in] IMOD Model number.
!> @param[in] SPECOLD
!> @param[inout] SPEC Spectrum (action) in 1-D form.
!> @param[out] VSIO
!> @param[out] VDIO
!> @param[out] SHAVEIO
!> @param[inout] ALPHA Nondimensional 1-D spectrum corresponding
!> to above full spectra (Phillip's const.).
!> @param[inout] WN1 Discrete wavenumbers.
!> @param[inout] CG1 Id. group velocities.
!> @param[in] CLATSL
!> @param[in] D_INP Depth, compared to DMIN to get DEPTH.
!> @param[in] U10ABS Wind speed at reference height.
!> @param[in] U10DIR Id. wind direction.
!> @param[inout] TAUA Magnitude of total atmospheric stress.
!> @param[inout] TAUADIR Direction of atmospheric stress.
!> @param[in] AS Air-sea temperature difference.
!> @param[inout] USTAR Friction velocity.
!> @param[inout] USTDIR Idem, direction.
!> @param[in] CX Current velocity component.
!> @param[in] CY Current velocity component.
!> @param[in] ICE Sea ice concentration.
!> @param[in] ICEH Sea ice thickness.
!> @param[inout] ICEF Sea ice floe diameter.
!> @param[in] ICEDMAX Sea ice maximum floe diameter
!> @param[in] REFLEC Reflection coefficients.
!> @param[in] REFLED Reflection direction.
!> @param[in] DELX Grid cell size in X direction.
!> @param[in] DELY Grid cell size in Y direction.
!> @param[in] DELA Grid cell area.
!> @param[in] TRNX Grid transparency in X.
!> @param[in] TRNY Grid transparency in Y.
!> @param[in] BERG Iceberg damping coefficient.
!> @param[inout] FPI Peak-input frequency.
!> @param[out] DTDYN Average dynamic time step.
!> @param[out] FCUT Cut-off frequency for tail.
!> @param[in] DTG Global time step.
!> @param[inout] TAUWX
!> @param[inout] TAUWY
!> @param[inout] TAUOX
!> @param[inout] TAUWIX
!> @param[inout] TAUWIY
!> @param[inout] TAUWNX
!> @param[inout] TAUWNY
!> @param[inout] PHIAW
!> @param[inout] CHARN
!> @param[inout] TWS
!> @param[inout] PHIOC
!> @param[inout] WHITECAP Whitecap statistics.
!> @param[in] D50 Sand grain size.
!> @param[in] PSIC Critical shields.
!> @param[inout] BEDFORM Bedform parameters.
!> @param[inout] PHIBBL Energy flux to BBL.
!> @param[inout] TAUBBL Momentum flux to BBL.
!> @param[inout] TAUICE Momentum flux to sea ice.
!> @param[inout] PHICE Energy flux to sea ice.
!> @param[inout] TAUOCX Total ocean momentum component.
!> @param[inout] TAUOCY Total ocean momentum component.
!> @param[inout] WNMEAN Mean wave number.
!> @param[in] DAIR Air density.
!> @param[in] COEF
!>
!> @author H. L. Tolman
!> @author F. Ardhuin
!> @author A. Roland
!> @author M. Dutour Sikiric
!> @date 22-Mar-2021
!>
SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
SPECOLD, SPEC, VSIO, VDIO, SHAVEIO, &
ALPHA, WN1, CG1, CLATSL, &
D_INP, U10ABS, U10DIR, &
#ifdef W3_FLX5
TAUA, TAUADIR, &
#endif
AS, USTAR, USTDIR, &
CX, CY, ICE, ICEH, ICEF, ICEDMAX, &
REFLEC, REFLED, DELX, DELY, DELA, TRNX, &
TRNY, BERG, FPI, DTDYN, FCUT, DTG, TAUWX, &
TAUWY, TAUOX, TAUOY, TAUWIX, TAUWIY, TAUWNX,&
TAUWNY, PHIAW, CHARN, TWS, PHIOC, WHITECAP, &
D50, PSIC, BEDFORM , PHIBBL, TAUBBL, TAUICE,&
PHICE, TAUOCX, TAUOCY, WNMEAN, DAIR, COEF)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | F. Ardhuin |
!/ | A. Roland |
!/ | M. Dutour Sikiric |
!/ | FORTRAN 90 |
!/ | Last update : 22-Mar-2021 |
!/ +-----------------------------------+
!/
!/ 06-Dec-1996 : Final FORTRAN 77 ( version 1.18 )
!/ 04-Feb-2000 : Upgrade to FORTRAN 90 ( version 2.00 )
!/ 14-Feb-2000 : Exact-NL added ( version 2.01 )
!/ 04-May-2000 : Non-central integration ( version 2.03 )
!/ 02-Feb-2001 : Xnl version 3.0 ( version 2.07 )
!/ 09-May-2002 : Switch clean up. ( version 2.21 )
!/ 13-Nov-2002 : Add stress vector. ( version 3.00 )
!/ 27-Nov-2002 : First version of VDIA and MDIA. ( version 3.01 )
!/ 07-Oct-2003 : Output options for NN training. ( version 3.05 )
!/ 24-Dec-2004 : Multiple model version. ( version 3.06 )
!/ 23-Jun-2006 : Linear input added. ( version 3.09 )
!/ 27-Jun-2006 : Adding file name preamble. ( version 3.09 )
!/ 04-Jul-2006 : Separation of stress computation. ( version 3.09 )
!/ 16-Apr-2007 : Miche style limiter added. ( version 3.11 )
!/ (J. H. Alves)
!/ 25-Apr-2007 : Battjes-Janssen Sdb added. ( version 3.11 )
!/ (J. H. Alves)
!/ 09-Oct-2007 : Adding WAM 4+ and SB1 options. ( version 3.13 )
!/ (F. Ardhuin)
!/ 29-May-2009 : Preparing distribution version. ( version 3.14 )
!/ 19-Aug-2010 : Making treatment of 0 water depth ( version 3.14.6 )
!/ consistent with the rest of the model.
!/ 31-Mar-2010 : Adding ice conc. and reflections ( version 3.14.4 )
!/ 15-May-2010 : Adding transparencies ( version 3.14.4 )
!/ 01-Jun-2011 : Movable bed bottom friction in BT4 ( version 4.01 )
!/ 01-Jul-2011 : Energy and momentum flux, friction ( version 4.01 )
!/ 24-Aug-2011 : Uses true depth for depth-induced ( version 4.04 )
!/ 16-Sep-2011 : Initialization of TAUWAX, TAUWAY ( version 4.04 )
!/ 1-Dec-2011 : Adding BYDRZ source term package ( version 4.04 )
!/ ST6 and optional Hwang (2011)
!/ stresses FLX4.
!/ 14-Mar-2012 : Update of BT4, passing PSIC ( version 4.04 )
!/ 13-Jul-2012 : Move GMD (SNL3) and nonlinear filter (SNLS)
!/ from 3.15 (HLT). ( version 4.08 )
!/ 28-Aug-2013 : Corrected MLIM application ( version 4.11 )
!/ 10-Sep-2013 : Special treatment for IG band ( version 4.15 )
!/ 14-Nov-2013 : Make orphaned pars in par lst local ( version 4.13 )
!/ 17-Nov-2013 : Coupling fraction of ice-free ( version 4.13 )
!/ surface to SIN and SDS. (S. Zieger)
!/ 01-Avr-2014 : Adding ice thickness and floe size ( version 4.18 )
!/ 23-May-2014 : Adding ice fluxes to W3SRCE ( version 5.01 )
!/ 27-Aug-2015 : Adding inputs to function W3SIS2 ( version 5.10 )
!/ 13-Dec-2015 : Implicit integration of Sice (F.A.) ( version 5.10 )
!/ 30-Jul-2017 : Adds TWS in interface ( version 6.04 )
!/ 07-Jan-2018 : Allows variable ice scaling (F.A.) ( version 6.04 )
!/ 01-Jan-2018 : Add implicit source term integration ( version 6.04)
!/ 01-Jan-2018 : within PDLIB (A. Roland, M. Dutour
!/ 18-Aug-2018 : S_{ice} IC5 (Q. Liu) ( version 6.06)
!/ 26-Aug-2018 : UOST (Mentaschi et al. 2015, 2018) ( version 6.06 )
!/ 22-Mar-2021 : Add extra fields used in coupling ( version 7.13 )
!/ 07-Jun-2021 : S_{nl5} GKE NL5 (Q. Liu) ( version 7.13 )
!/ 19-Jul-2021 : Momentum and air density support ( version 7.14 )
!/
!/ Copyright 2009-2013 National Weather Service (NWS),
!/ National Oceanic and Atmospheric Administration. All rights
!/ reserved. WAVEWATCH III is a trademark of the NWS.
!/ No unauthorized use without permission.
!/
! 1. Purpose :
!
! Calculate and integrate source terms for a single grid point.
!
! 2. Method :
!
! Physics : see manual and corresponding subroutines.
!
! Numerics :
!
! Dynamic-implicit integration of the source terms based on
! WW-II (Tolman 1992). The dynamic time step is calculated
! given a maximum allowed change of spectral densities for
! frequencies / wavenumbers below the usual cut-off.
! The maximum change is given by the minimum of a parametric
! and a relative change. The parametric change relates to a
! PM type equilibrium range
!
! -1 (2pi)**4 1
! dN(k) = Xp alpha pi ---------- ------------
! max g**2 k**3 sigma
!
! 1 .
! = FACP ------------ (1)
! k**3 sigma .
!
! where
! alpha = 0.62e-4 (set in W3GRID)
! Xp fraction of PM shape (read in W3GRID)
! FACP combined factor (set in W3GRID)
!
! The maximum relative change is given as
!
! / +- -+ \ .
! dN(k) = Xr max | N(k) , max | Nx , Xfilt N(k) | | (2)
! max \ +- max-+ / .
!
! where
! Xr fraction of relative change (read in W3GRID)
! Xfilt filter level (read in W3GRID)
! Nx Maximum parametric change (1)
! for largest wavenumber.
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! IX,IY Int. I Discrete grid point counters.
! IMOD Int. I Model number.
! SPEC R.A. I/O Spectrum (action) in 1-D form.
! ALPHA R.A. I/O Nondimenional 1-D spectrum corresponding
! to above full spectra (Phillip's const.).
! Calculated separately for numerical
! economy on vector machine (W3SPR2).
! WN1 R.A. I Discrete wavenumbers.
! CG1 R.A. I Id. group velocities.
! D_INP Real. I Depth. Compared to DMIN to get DEPTH.
! U10ABS Real. I Wind speed at reference height.
! U10DIR Real. I Id. wind direction.
! TAUA Real. I Magnitude of total atmospheric stress ( !/FLX5 )
! TAUADIR Real. I Direction of atmospheric stress ( !/FLX5 )
! AS Real. I Air-sea temp. difference. ( !/ST3 )
! USTAR Real. !/O Friction velocity.
! USTDIR Real !/O Idem, direction.
! CX-Y Real. I Current velocity components. ( !/BS1 )
! ICE Real I Sea ice concentration
! ICEH Real I Sea ice thickness
! ICEF Real I/O Sea ice maximum floe diameter (updated)
! ICEDMAX Real I/O Sea ice maximum floe diameter
! BERG Real I Iceberg damping coefficient ( !/BS1 )
! REFLEC R.A. I reflection coefficients ( !/BS1 )
! REFLED I.A. I reflection direction ( !/BS1 )
! TRNX-Y Real I Grid transparency in X and Y ( !/BS1 )
! DELX Real. I grid cell size in X direction ( !/BS1 )
! DELY Real. I grid cell size in Y direction ( !/BS1 )
! DELA Real. I grid cell area ( !/BS1 )
! FPI Real I/O Peak-input frequency. ( !/ST2 )
! WHITECAP R.A. O Whitecap statisics ( !/ST4 )
! DTDYN Real O Average dynamic time step.
! FCUT Real O Cut-off frequency for tail.
! DTG Real I Global time step.
! D50 Real I Sand grain size ( !/BT4 )
! BEDFORM R.A. I/O Bedform parameters ( !/BT4 )
! PSIC Real I Critical Shields ( !/BT4 )
! PHIBBL Real O Energy flux to BBL ( !/BTx )
! TAUBBL R.A. O Momentum flux to BBL ( !/BTx )
! TAUICE R.A. O Momentum flux to sea ice ( !/ICx )
! PHICE Real O Energy flux to sea ice ( !/ICx )
! TAUOCX-YReal O Total ocean momentum components
! WNMEAN Real O Mean wave number
! DAIR Real I Air density
! ----------------------------------------------------------------
! Note: several pars are set to I/O to avoid compiler warnings.
!
! 4. Subroutines used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3SPRn Subr. W3SRCnMD Mean wave parameters for use in
! source terms.
! W3FLXn Subr. W3FLXnMD Flux/stress computation.
! W3SLNn Subr. W3SLNnMD Linear input.
! W3SINn Subr. W3SRCnMD Input source term.
! W3SNLn Subr. W3SNLnMD Nonlinear interactions.
! W3SNLS Subr. W3SNLSMD Nonlinear smoother.
! W3SDSn Subr. W3SRCnMD Whitecapping source term
! W3SBTn Subr. W3SBTnMD Bottom friction source term.
! W3SDBn Subr. W3SBTnMD Depth induced breaking source term.
! W3STRn Subr. W3STRnMD Triad interaction source term.
! W3SBSn Subr. W3SBSnMD Bottom scattering source term.
! W3REFn Subr. W3REFnMD Reflexions (shore, icebergs ...).
! STRACE Subr. W3SERVMD Subroutine tracing (!/S)
! ----------------------------------------------------------------
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3WAVE Subr. W3WAVEMD Actual wave model routine.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! - No testing is performed on the status of the grid point.
!
! 8. Structure :
!
! -----------------------------------------------------------------
! 1. Preparations
! a Set maximum change and wavenumber arrays.
! b Prepare dynamic time stepping.
! c Compute mean parameters. ( W3SPRn )
! d Compute stresses (if posible).
! e Prepare cut-off
! f Test output for !/NNT option.
! --start-dynamic-integration-loop---------------------------------
! 2. Calculate source terms
! a Input. ( W3SLNx, W3SINn )
! b Nonlinear interactions. ( W3SNLn )
! c Dissipation ( W3SDSn )
! 1 as included in source terms ( W3SDSn )
! 2 optional dissipation due to different physics ( W3SWLn )
! d Bottom friction. ( W3SBTn )
! 3. Calculate cut-off frequencie(s)
! 4. Summation of source terms and diagonal term and time step.
! 5. Increment spectrum.
! 6. Add tail
! a Mean wave parameters and cut-off ( W3SPRn )
! b 'Seeding' of spectrum. ( !/SEED )
! c Add tail
! 7. Check if integration complete.
! --end-dynamic-integration-loop-----------------------------------
! 8. Save integration data.
! -----------------------------------------------------------------
!
! 9. Switches :
!
! !/FLX1 Wu (1980) stress computation. ( Choose one )
! !/FLX2 T&C (1996) stress computation.
! !/FLX3 T&C (1996) stress computation with cap.
! !/FLX4 Hwang (2011) stress computation (2nd order).
! !/FLX5 Direct use of stress from atmoshperic model.
!
! !/LN0 No linear input. ( Choose one )
!
! !/ST0 No input and dissipation. ( Choose one )
! !/ST1 WAM-3 input and dissipation.
! !/ST2 Tolman and Chalikov (1996) input and dissipation.
! !/ST3 WAM 4+ input and dissipation.
! !/ST4 Ardhuin et al. (2009, 2010)
! !/ST6 BYDB source terms after Babanin, Young, Donelan and Banner.
!
! !/NL0 No nonlinear interactions. ( Choose one )
! !/NL1 Discrete interaction approximation.
! !/NL2 Exact nonlinear interactions.
! !/NL3 Generalized Multiple DIA.
! !/NL4 Two Scale Approximation
! !/NL5 Generalized Kinetic Equation.
! !/NLS Nonlinear HF smoother.
!
! !/BT0 No bottom friction. ( Choose one )
! !/BT1 JONSWAP bottom friction.
! !/BT4 Bottom friction using movable bed roughness
! (Tolman 1994, Ardhuin & al. 2003)
! !/BT8 Muddy bed (Dalrymple & Liu).
! !/BT9 Muddy bed (Ng).
!
! !/IC1 Dissipation via interaction with ice according to simple
! methods: 1) uniform in frequency or
! !/IC2 2) Liu et al. model
! !/IC3 Dissipation via interaction with ice according to a
! viscoelastic sea ice model (Wang and Shen 2010).
! !/IC4 Dissipation via interaction with ice as a function of freq.
! (empirical/parametric methods)
! !/IC5 Dissipation via interaction with ice according to a
! viscoelastic sea ice model (Mosig et al. 2015).
! !/DB0 No depth-limited breaking. ( Choose one )
! !/DB1 Battjes-Janssen depth-limited breaking.
!
! !/TR0 No triad interactions. ( Choose one )
! !/TR1 Lumped Triad Approximation (LTA).
!
! !/BS0 No bottom scattering. ( Choose one )
! !/BS1 Scattering term by Ardhuin and Magne (2007).
!
! !/MLIM Miche style limiter for shallow water and steepness.
!
! !/SEED 'Seeding' of lowest frequency for suffuciently strong
! winds.
!
! !/NNT Write output to file test_data_NNN.ww3 for NN training.
!
! !/S Enable subroutine tracing.
! !/T Enable general test output.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE CONSTANTS, ONLY: DWAT, srce_imp_post, srce_imp_pre, &
srce_direct, GRAV, TPI, TPIINV
USE W3GDATMD, ONLY: NK, NTH, NSPEC, SIG, TH, DMIN, DTMAX, &
DTMIN, FACTI1, FACTI2, FACSD, FACHFA, FACP, &
XFC, XFLT, XREL, XFT, FXFM, FXPM, DDEN, &
FTE, FTF, FHMAX, ECOS, ESIN, IICEDISP, &
ICESCALES, IICESMOOTH
USE W3WDATMD, ONLY: TIME
USE W3ODATMD, ONLY: NDSE, NDST, IAPROC
USE W3IDATMD, ONLY: INFLAGS2
USE W3DISPMD
#ifdef W3_T
USE CONSTANTS, ONLY: RADE
#endif
#ifdef W3_REF1
USE W3GDATMD, ONLY: IOBP, IOBPD, GTYPE, UNGTYPE, REFPARS
#endif
#ifdef W3_NNT
USE W3ODATMD, ONLY: IAPROC, SCREEN, FNMPRE
#endif
#ifdef W3_FLD1
USE W3FLD1MD, ONLY: W3FLD1
USE W3GDATMD, ONLY: AALPHA
#endif
#ifdef W3_FLD2
USE W3FLD2MD, ONLY: W3FLD2
USE W3GDATMD, ONLY: AALPHA
#endif
#ifdef W3_FLX1
USE W3FLX1MD
#endif
#ifdef W3_FLX2
USE W3FLX2MD
#endif
#ifdef W3_FLX3
USE W3FLX3MD
#endif
#ifdef W3_FLX4
USE W3FLX4MD
#endif
#ifdef W3_FLX5
USE W3FLX5MD
#endif
#ifdef W3_LN1
USE W3SLN1MD
#endif
#ifdef W3_ST0
USE W3SRC0MD
#endif
#ifdef W3_ST1
USE W3SRC1MD
#endif
#ifdef W3_ST2
USE W3SRC2MD
USE W3GDATMD, ONLY : ZWIND
#endif
#ifdef W3_ST3
USE W3SRC3MD
USE W3GDATMD, ONLY : ZZWND, FFXFM, FFXPM
#endif
#ifdef W3_ST4
USE W3SRC4MD, ONLY : W3SPR4, W3SIN4, W3SDS4
USE W3GDATMD, ONLY : ZZWND, FFXFM, FFXPM, FFXFA
#endif
#ifdef W3_ST6
USE W3SRC6MD
USE W3SWLDMD, ONLY : W3SWL6
USE W3GDATMD, ONLY : SWL6S6
#endif
#ifdef W3_NL1
USE W3SNL1MD
#endif
#ifdef W3_NL2
USE W3SNL2MD
#endif
#ifdef W3_NL3
USE W3SNL3MD
#endif
#ifdef W3_NL4
USE W3SNL4MD
#endif
#ifdef W3_NL5
USE W3SNL5MD
USE W3TIMEMD, ONLY: TICK21
#endif
#ifdef W3_NLS
USE W3SNLSMD
#endif
#ifdef W3_BT1
USE W3SBT1MD
#endif
#ifdef W3_BT4
USE W3SBT4MD
#endif
#ifdef W3_BT8
USE W3SBT8MD
#endif
#ifdef W3_BT9
USE W3SBT9MD
#endif
#ifdef W3_IC1
USE W3SIC1MD
#endif
#ifdef W3_IC2
USE W3SIC2MD
#endif
#ifdef W3_IC3
USE W3SIC3MD
#endif
#ifdef W3_IC4
USE W3SIC4MD
#endif
#ifdef W3_IC5
USE W3SIC5MD
#endif
#ifdef W3_IS1
USE W3SIS1MD
#endif
#ifdef W3_IS2
USE W3SIS2MD
USE W3GDATMD, ONLY : IS2PARS
#endif
#ifdef W3_DB1
USE W3SDB1MD
#endif
#ifdef W3_TR1
USE W3STR1MD
#endif
#ifdef W3_BS1
USE W3SBS1MD
#endif
#ifdef W3_REF1
USE W3REF1MD
#endif
#ifdef W3_IG1
USE W3GDATMD, ONLY : IGPARS
#endif
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
#ifdef W3_NNT
USE W3SERVMD, ONLY: EXTCDE
#endif
#ifdef W3_UOST
USE W3UOSTMD, ONLY: UOST_SRCTRMCOMPUTE
#endif
#ifdef W3_PDLIB
USE PDLIB_W3PROFSMD, ONLY : B_JAC, ASPAR_JAC, ASPAR_DIAG_ALL
USE yowNodepool, ONLY: PDLIB_I_DIAG, PDLIB_SI
USE W3GDATMD, ONLY: B_JGS_LIMITER, FSSOURCE, optionCall
USE W3GDATMD, ONLY: IOBP_LOC, IOBPD_LOC, B_JGS_LIMITER_FUNC
USE W3WDATMD, ONLY: VA
USE W3PARALL, ONLY: IMEM, LSLOC
#endif
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
INTEGER, INTENT(IN) :: srce_call, IT, ISEA, JSEA, IX, IY, IMOD
REAL, intent(in) :: SPECOLD(NSPEC), CLATSL
REAL, INTENT(OUT) :: VSIO(NSPEC), VDIO(NSPEC)
LOGICAL, INTENT(OUT) :: SHAVEIO
REAL, INTENT(IN) :: D_INP, U10ABS, &
U10DIR, AS, CX, CY, DTG, D50,PSIC, &
ICE, ICEH
#ifdef W3_FLX5
REAL, INTENT(IN) :: TAUA, TAUADIR
#endif
INTEGER, INTENT(IN) :: REFLED(6)
REAL, INTENT(IN) :: REFLEC(4), DELX, DELY, DELA, &
TRNX, TRNY, BERG, ICEDMAX, DAIR
REAL, INTENT(INOUT) :: WN1(NK), CG1(NK), &
SPEC(NSPEC), ALPHA(NK), USTAR, &
USTDIR, FPI, TAUOX, TAUOY, &
TAUWX, TAUWY, PHIAW, PHIOC, PHICE, &
CHARN, TWS, BEDFORM(3), PHIBBL, &
TAUBBL(2), TAUICE(2), WHITECAP(4), &
TAUWIX, TAUWIY, TAUWNX, TAUWNY, &
ICEF, TAUOCX, TAUOCY, WNMEAN
REAL, INTENT(OUT) :: DTDYN, FCUT
REAL, INTENT(IN) :: COEF
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
INTEGER :: IK, ITH, IS, IS0, NSTEPS, NKH, NKH1, &
IKS1, IS1, NSPECH, IDT, IERR, ISP
REAL :: DTTOT, FHIGH, DT, AFILT, DAMAX, AFAC, &
HDT, ZWND, FP, DEPTH, TAUSCX, TAUSCY, FHIGI
! Scaling factor for SIN, SDS, SNL
REAL :: ICESCALELN, ICESCALEIN, ICESCALENL, ICESCALEDS
REAL :: EMEAN, FMEAN, AMAX, CD, Z0, SCAT, &
SMOOTH_ICEDISP
REAL :: WN_R(NK), CG_ICE(NK), ALPHA_LIU(NK), ICECOEF2, R(NK)
DOUBLE PRECISION :: ATT, ISO
REAL :: EBAND, DIFF, EFINISH, HSTOT, PHINL, &
FMEAN1, FMEANWS, &
FACTOR, FACTOR2, DRAT, TAUWAX, TAUWAY, &
MWXFINISH, MWYFINISH, A1BAND, B1BAND, &
COSI(2)
REAL :: SPECINIT(NSPEC), SPEC2(NSPEC), FRLOCAL, JAC2
REAL :: DAM (NSPEC), DAM2(NSPEC), WN2(NSPEC), &
VSLN(NSPEC), &
VSIN(NSPEC), VDIN(NSPEC), &
VSNL(NSPEC), VDNL(NSPEC), &
VSDS(NSPEC), VDDS(NSPEC), &
VSBT(NSPEC), VDBT(NSPEC)
REAL :: VS(NSPEC), VD(NSPEC), EB(NK)
LOGICAL :: SHAVE
LOGICAL :: LBREAK
LOGICAL, SAVE :: FIRST = .TRUE.
LOGICAL :: PrintDeltaSmDA
REAL :: eInc1, eInc2, eVS, eVD, JAC
REAL :: DeltaSRC(NSPEC)
REAL :: FOUT(NK,NTH), SOUT(NK,NTH), DOUT(NK,NTH)
REAL, SAVE :: TAUNUX, TAUNUY
LOGICAL, SAVE :: FLTEST = .FALSE., FLAGNN = .TRUE.
#ifdef W3_OMPG
!$omp threadprivate( TAUNUX, TAUNUY)
!$omp threadprivate( FLTEST, FLAGNN )
!$omp threadprivate( FIRST )
#endif
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters dependent on compile switch
!/
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
#endif
#ifdef W3_NNT
INTEGER, SAVE :: NDSD = 89, NDSD2 = 88, J
REAL :: QCERR = 0. !/XNL2 and !/NNT
#endif
#ifdef W3_NL5
INTEGER :: QI5TSTART(2)
REAL :: QR5KURT
INTEGER, PARAMETER :: NL5_SELECT = 1
REAL, PARAMETER :: NL5_OFFSET = 0. ! explicit dyn.
#endif
#ifdef W3_SEED
REAL :: UC, SLEV
#endif
#ifdef W3_MLIM
REAL :: HM, EM
#endif
#ifdef W3_NNT
REAL :: FACNN
#endif
#ifdef W3_T
REAL :: DTRAW
#endif
#if defined(W3_IC1) || defined(W3_IC2) || defined(W3_IC3) || defined(W3_IC4) || defined(W3_IC5)
REAL :: VSIC(NSPEC), VDIC(NSPEC)
#endif
#ifdef W3_DB1
REAL :: VSDB(NSPEC), VDDB(NSPEC)
#endif
#ifdef W3_TR1
REAL :: VSTR(NSPEC), VDTR(NSPEC)
#endif
#ifdef W3_BS1
REAL :: VSBS(NSPEC), VDBS(NSPEC)
#endif
#ifdef W3_REF1
REAL :: VREF(NSPEC)
#endif
#if defined(W3_IS1) || defined(W3_IS2)
REAL :: VSIR(NSPEC), VDIR(NSPEC)
#endif
#ifdef W3_IS2
REAL :: VDIR2(NSPEC)
DOUBLE PRECISION :: SCATSPEC(NTH)
#endif
#ifdef W3_UOST
REAL :: VSUO(NSPEC), VDUO(NSPEC)
#endif
#ifdef W3_ST1
REAL :: FH1, FH2
#endif
#ifdef W3_ST2
REAL :: FHTRAN, DFH, FACDIA, FACPAR
#endif
#ifdef W3_ST3
REAL :: FMEANS, FH1, FH2
#endif
#ifdef W3_ST4
REAL :: FMEANS, FH1, FH2, FAGE, DLWMEAN
REAL :: BRLAMBDA(NSPEC)
#endif
#if defined(W3_ST3) || defined(W3_ST4)
LOGICAL :: LLWS(NSPEC)
#endif
#ifdef W3_ST6
REAL :: VSWL(NSPEC), VDWL(NSPEC)
#endif
#ifdef W3_PDLIB
REAL :: PreVS, DVS, SIDT, FAKS, MAXDAC
#endif
#ifdef W3_NNT
CHARACTER(LEN=17), SAVE :: FNAME = 'test_data_nnn.ww3'
#endif
!
!/ -- End of variable delclarations
!
#ifdef W3_S
CALL STRACE (IENT, 'W3SRCE')
#endif
#ifdef W3_T
FLTEST = .TRUE.
#endif
!
IKS1 = 1
#ifdef W3_IG1
! Does not integrate source terms for IG band if IGPARS(12) = 0.
IF (NINT(IGPARS(12)).EQ.0) IKS1 = NINT(IGPARS(5))
#endif
IS1=(IKS1-1)*NTH+1
!! Initialise source term arrays:
VD = 0.
VS = 0.
VDIO = 0.
VSIO = 0.
VSBT = 0.
VDBT = 0.
#if defined(W3_LN0) || defined(W3_LN1) || defined(W3_SEED)
VSLN = 0.
#endif
#if defined(W3_ST0) || defined(W3_ST3) || defined(W3_ST4)
VSIN = 0.
VDIN = 0.
#endif
#if defined(W3_NL0) || defined(W3_NL1)
VSNL = 0.
VDNL = 0.
#endif
#ifdef W3_TR1
VSTR = 0.
VDTR = 0.
#endif
#if defined(W3_ST0) || defined(W3_ST4)
VSDS = 0.
VDDS = 0.
#endif
#ifdef W3_DB1
VSDB = 0.
VDDB = 0.
#endif
#if defined(W3_IC1) || defined(W3_IC2) || defined(W3_IC3) || defined(W3_IC4) || defined(W3_IC5)
VSIC = 0.
VDIC = 0.
#endif
#ifdef W3_UOST
VSUO = 0.
VDUO = 0.
#endif
#if defined(W3_IS1) || defined(W3_IS2)
VSIR = 0.
VDIR = 0.
#endif
#ifdef W3_IS2
VDIR2 = 0.
#endif
#ifdef W3_ST6
VSWL = 0.
VDWL = 0.
#endif
#if defined(W3_ST0) || defined(W3_ST1) || defined(W3_ST6)
ZWND = 10.
#endif
#if defined(W3_ST2)
ZWND = ZWIND
#endif
#if defined(W3_ST4)
ZWND = ZZWND
#endif
!
! 1. Preparations --------------------------------------------------- *
!
DEPTH = MAX ( DMIN , D_INP )
DRAT = DAIR / DWAT
ICESCALELN = MAX(0.,MIN(1.,1.-ICE*ICESCALES(1)))
ICESCALEIN = MAX(0.,MIN(1.,1.-ICE*ICESCALES(2)))
ICESCALENL = MAX(0.,MIN(1.,1.-ICE*ICESCALES(3)))
ICESCALEDS = MAX(0.,MIN(1.,1.-ICE*ICESCALES(4)))
#ifdef W3_T
WRITE (NDST,9000)
WRITE (NDST,9001) DEPTH, U10ABS, U10DIR*RADE
#endif
! 1.a Set maximum change and wavenumber arrays.
!
!XP = 0.15
!FACP = XP / PI * 0.62E-3 * TPI**4 / GRAV**2
!
DO IK=1, NK
DAM(1+(IK-1)*NTH) = FACP / ( SIG(IK) * WN1(IK)**3 )
WN2(1+(IK-1)*NTH) = WN1(IK)
END DO
!
DO IK=1, NK
IS0 = (IK-1)*NTH
DO ITH=2, NTH
DAM(ITH+IS0) = DAM(1+IS0)
WN2(ITH+IS0) = WN2(1+IS0)
END DO
END DO
!
! 1.b Prepare dynamic time stepping
!
DTDYN = 0.
DTTOT = 0.
NSTEPS = 0
PHIAW = 0.
CHARN = 0.
TWS = 0.
PHINL = 0.
PHIBBL = 0.
TAUWIX = 0.
TAUWIY = 0.
TAUWNX = 0.
TAUWNY = 0.
TAUWAX = 0.
TAUWAY = 0.
TAUSCX = 0.
TAUSCY = 0.
TAUBBL = 0.
TAUICE = 0.
PHICE = 0.
TAUOCX = 0.
TAUOCY = 0.
WNMEAN = 0.
!
! TIME is updated in W3WAVEMD prior to the call of W3SCRE, we should
! move 'TIME' one time step backward (QL)
#ifdef W3_NL5
QI5TSTART = TIME
CALL TICK21 (QI5TSTART, -1.0 * DTG)
#endif
!
#ifdef W3_DEBUGSRC
IF (IX .eq. DEBUG_NODE) THEN
WRITE(740+IAPROC,*) 'W3SRCE start sum(SPEC)=', sum(SPEC)
WRITE(740+IAPROC,*) 'W3SRCE start sum(SPECOLD)=', sum(SPECOLD)
WRITE(740+IAPROC,*) 'W3SRCE start sum(SPECINIT)=', sum(SPECINIT)
WRITE(740+IAPROC,*) 'W3SRCE start sum(VSIO)=', sum(VSIO)
WRITE(740+IAPROC,*) 'W3SRCE start sum(VDIO)=', sum(VDIO)
WRITE(740+IAPROC,*) 'W3SRCE start USTAR=', USTAR
END IF
#endif
#ifdef W3_ST4
DLWMEAN= 0.
BRLAMBDA(:)=0.
WHITECAP(:)=0.
#endif
!
! 1.c Set mean parameters
!
#ifdef W3_ST0