-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathw3pro3md.ftn
1858 lines (1847 loc) · 65.1 KB
/
w3pro3md.ftn
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 "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE W3PRO3MD
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 27-May-2014 |
!/ +-----------------------------------+
!/
!/ 27-Feb-2000 : Origination. ( version 2.08 )
!/ 17-Sep-2000 : Clean-up. ( version 2.13 )
!/ 10-Dec-2001 : Sub-grid obstructions. ( version 2.14 )
!/ 16-Oct-2002 : Change INTENT for ATRN in W3XYP3. ( version 3.00 )
!/ 26-Dec-2002 : Moving grid version. ( version 3.02 )
!/ 01-Aug-2003 : Moving grid GSE correction. ( version 3.03 )
!/ 17-Dec-2004 : Multiple grid version. ( version 3.06 )
!/ 07-Sep-2005 : Upgrade XY boundary conditions. ( version 3.08 )
!/ 09-Nov-2005 : Removing soft boundary option. ( version 3.08 )
!/ 05-Mar-2008 : Added NEC sxf90 compiler directives.
!/ (Chris Bunney, UK Met Office) ( version 3.13 )
!/ 01-Apr-2008 : Bug fix W3MAP3 MAPSTA range check. ( version 3.13 )
!/ 29-May-2009 : Preparing distribution version. ( version 3.14 )
!/ 30-Oct-2009 : Implement curvilinear grid type. ( version 3.14 )
!/ (W. E. Rogers & T. J. Campbell, NRL)
!/ 17-Aug-2010 : Add test output W3XYP3. ( version 3.14.5 )
!/ 06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
!/ specify index closure for a grid. ( version 3.14 )
!/ (T. J. Campbell, NRL)
!/ 26-Dec-2012 : More initializations. ( version 4.11 )
!/ 01-Jul-2013 : Adding UQ and UNO switches to chose between third
!/ and second order schemes. ( version 4.12 )
!/ 12-Sep-2013 : Add documentation for global clos. ( version 4.12 )
!/ 27-May-2014 : Adding OMPH switch. ( version 5.02 )
!/
!/ Copyright 2009-2014 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 :
!
! Bundles routines for third order propagation scheme in single
! module.
!
! 2. Variables and types :
!
! Name Type Scope Description
! ----------------------------------------------------------------
! TRNMIN R.P. Private Minimum transparancy for local
! switching off of averaging.
! ----------------------------------------------------------------
!
! Also work arrays for W3KTP3 (private).
!
! 3. Subroutines and functions :
!
! Name Type Scope Description
! ----------------------------------------------------------------
! W3MAP3 Subr. Public Set up auxiliary maps.
! W3MAPT Subr. Public Set up transparency map for GSE.
! W3XYP3 Subr. Public Third order spatial propagation.
! W3KTP3 Subr. Public Third order spectral propagation.
! ----------------------------------------------------------------
!
! 4. Subroutines and functions used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! W3QCK1 Subr. W3UQCKMD Regular grid UQ scheme.
! W3QCK2 Subr. Id. Irregular grid UQ scheme.
! W3QCK3 Subr. Id. Regular grid UQ scheme + obstructions.
! W3UNO2 Subr. W3UNO2MD UNO2 scheme for irregular grid.
! W3UNO2r Subr. Id. UNO2 scheme reduced to regular grid.
! W3UNO2s Subr. Id. UNO2 regular grid with subgrid
! obstruction.
! ----------------------------------------------------------------
!
! 5. Remarks :
!
! - The averaging is not performed around semi-transparent grid
! points to avoid that leaking through barriers occurs.
!
! 6. Switches :
!
! !/UQ 3rd order UQ propagation scheme.
! !/UNO 2nd order UNO propagation scheme.
!
! !/C90 Cray FORTRAN 90 compiler directives.
! !/NEC NEC SXF90 compiler directives.
!
! !/MGP Moving grid corrections.
! !/MGG Moving grid corrections.
!
! !/OMPH Hybrid OpenMP directives.
!
! !/S Enable subroutine tracing.
! !/Tn Enable test output.
!
! 7. Source code :
!
!/ ------------------------------------------------------------------- /
!/
!/ Public variables
!/
PUBLIC
!/
!/ Private data
!/
REAL, PRIVATE, PARAMETER:: TRNMIN = 0.95
!/
CONTAINS
!/ ------------------------------------------------------------------- /
SUBROUTINE W3MAP3
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 01-Apr-2008 |
!/ +-----------------------------------+
!/
!/ 27-Feb-2000 : Origination. ( version 2.08 )
!/ 10-Dec-2001 : Sub-grid obstructions. ( version 2.14 )
!/ (array allocation only.)
!/ 17-Dec-2004 : Multiple grid version. ( version 3.06 )
!/ 09-Nov-2005 : Removing soft boundary option. ( version 3.08 )
!/ 01-Apr-2008 : Bug fix sec. 4 MAPSTA range check. ( version 3.13 )
!/
! 1. Purpose :
!
! Generate 'map' arrays for the ULTIMATE QUICKEST scheme.
!
! 2. Method :
!
! MAPX2, MAPY2, MAPTH2 and MAPWN2 contain consecutive 1-D spatial
! grid counters (e.g., IXY = (IX-1)*MY + IY). The arrays are
! devided in three parts. For MAPX2, these ranges are :
!
! 1 - NMX0 Counters IXY for which grid point (IX,IY) and
! (IX+1,IY) both are active grid points.
! NMX0+1 - NMX1 Id. only (IX,IY) active.
! NMX1+1 - NMX2 Id. only (IX+1,IY) active.
!
! The array MAPY2 has a similar layout varying IY instead of IX.
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! See module documentation.
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3WAVE Subr. W3WAVEMD Wave model routine.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! 7. Remarks :
!
! 8. Structure :
!
! ------------------------------------------------------
! 1. Map MAPX2
! a Range 1 to NMX0
! b Range NMX0+1 to NMX1
! c Range NMX1+1 to NMX2
! 2. Map MAPY2
! a Range 1 to NMY0
! b Range NMY0+1 to NMY1
! c Range NMY1+1 to NMY2
! 3. Map MAPAXY
! 4. Map MAPCXY
! 5. Maps for intra-spectral propagation
! a MAPTH2, MAPATK
! b MAPWN2
! ------------------------------------------------------
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
! !/T Enable test output.
!
! 10. Source code :
!/ ------------------------------------------------------------------- /
USE W3GDATMD, ONLY: NK, NTH, NSPEC, NX, NY, NSEA, MAPSTA, MAPSF,&
GTYPE
USE W3ADATMD, ONLY: NMX0, NMX1, NMX2, NMY0, NMY1, NMY2, NACT, &
NCENT, MAPX2, MAPY2, MAPAXY, MAPCXY, &
MAPTH2, MAPWN2
USE W3ODATMD, ONLY: NDST
!/S USE W3SERVMD, ONLY: STRACE
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
INTEGER :: IX, IY, IXY0, IX2, IY2, IX0, IY0, &
ISEA, IK, ITH, ISP, ISP0, ISP2, NCENTC
!/S INTEGER, SAVE :: IENT = 0
!/T INTEGER :: MAPTXY(NY,NX), I, IXY
!/T INTEGER :: MAPTST(NK+2,NTH)
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'W3MAP3')
!
IF (GTYPE .LT. 3) THEN
! 1. Map MAPX2 ------------------------------------------------------ *
! 1.a Range 1 to NMX0
!
!/T MAPTXY = 0.
!
NMX0 = 0
DO IX=1, NX
IXY0 = (IX-1)*NY
IX2 = 1 + MOD(IX,NX)
DO IY=2, NY-1
IF ( MAPSTA(IY,IX).EQ.1 .AND. MAPSTA(IY,IX2).EQ.1 ) THEN
NMX0 = NMX0 + 1
MAPX2(NMX0) = IXY0 + IY
!/T MAPTXY(IY,IX) = MAPTXY(IY,IX) + 1
END IF
END DO
END DO
!
! 1.b Range NMX0+1 to NMX1
!
NMX1 = NMX0
DO IX=1, NX
IXY0 = (IX-1)*NY
IX2 = 1 + MOD(IX,NX)
DO IY=2, NY-1
IF ( MAPSTA(IY,IX).EQ.1 .AND. MAPSTA(IY,IX2).NE.1 ) THEN
NMX1 = NMX1 + 1
MAPX2(NMX1) = IXY0 + IY
!/T MAPTXY(IY,IX) = MAPTXY(IY,IX) + 2
END IF
END DO
END DO
!
! 1.c Range NMX1+1 to NMX2
!
NMX2 = NMX1
DO IX=1, NX
IXY0 = (IX-1)*NY
IX2 = 1 + MOD(IX,NX)
DO IY=2, NY-1
IF ( MAPSTA(IY,IX).NE.1 .AND. MAPSTA(IY,IX2).EQ.1 ) THEN
NMX2 = NMX2 + 1
MAPX2(NMX2) = IXY0 + IY
!/T MAPTXY(IY,IX) = MAPTXY(IY,IX) + 4
END IF
END DO
END DO
!
!/T WRITE (NDST,9000) 'MAPX2', NMX0, NMX1-NMX0, &
!/T NMX2-NMX1, NMX2
!/T DO IY=NY, 1, -1
!/T WRITE (NDST,9001) (MAPTXY(IY,IX),IX=1, NX)
!/T END DO
!
! 2. Map MAPY2 ------------------------------------------------------ *
! 2.a Range 1 to NMY0
!
!/T MAPTXY = 0.
!
NMY0 = 0
DO IX=1, NX
IXY0 = (IX-1)*NY
DO IY=1, NY-1
IY2 = IY + 1
IF ( MAPSTA(IY,IX).EQ.1 .AND. MAPSTA(IY2,IX).EQ.1 ) THEN
NMY0 = NMY0 + 1
MAPY2(NMY0) = IXY0 + IY
!/T MAPTXY(IY,IX) = MAPTXY(IY,IX) + 1
END IF
END DO
END DO
!
! 2.b Range NMY0+1 to NMY1
!
NMY1 = NMY0
DO IX=1, NX
IXY0 = (IX-1)*NY
DO IY=1, NY-1
IY2 = IY + 1
IF ( MAPSTA(IY,IX).EQ.1 .AND. MAPSTA(IY2,IX).NE.1 ) THEN
NMY1 = NMY1 + 1
MAPY2(NMY1) = IXY0 + IY
!/T MAPTXY(IY,IX) = MAPTXY(IY,IX) + 2
END IF
END DO
END DO
!
! 2.c Range NMY1+1 to NMY2
!
NMY2 = NMY1
DO IX=1, NX
IXY0 = (IX-1)*NY
DO IY=1, NY-1
IY2 = IY + 1
IF ( MAPSTA(IY,IX).NE.1 .AND. MAPSTA(IY2,IX).EQ.1 ) THEN
NMY2 = NMY2 + 1
MAPY2(NMY2) = IXY0 + IY
!/T MAPTXY(IY,IX) = MAPTXY(IY,IX) + 4
END IF
END DO
END DO
!
!/T WRITE (NDST,9000) 'MAPY2', NMY0, NMY1-NMY0, &
!/T NMY2-NMY1, NMY2
!/T DO IY=NY, 1, -1
!/T WRITE (NDST,9001) (MAPTXY(IY,IX),IX=1, NX)
!/T END DO
!
! 3. Map MAPAXY ----------------------------------------------------- *
!
NACT = 0
DO IX=1, NX
IY0 = (IX-1)*NY
DO IY=2, NY-1
IF ( MAPSTA(IY,IX).EQ.1 ) THEN
NACT = NACT + 1
MAPAXY(NACT) = IY0 + IY
END IF
END DO
END DO
!
! 4. Map MAPCXY ----------------------------------------------------- *
!
NCENT = 0
NCENTC = NSEA
MAPCXY = 0
!
DO ISEA=1, NSEA
IX = MAPSF(ISEA,1)
IX0 = IX-1
IX2 = IX+1
IY = MAPSF(ISEA,2)
IY0 = IY-1
IY2 = IY+1
IF ( IX .EQ. NX ) IX2 = 1
IF ( IX .EQ. 1 ) IX0 = NX
IF ( MAPSTA(IY,IX).EQ.2 .OR. MAPSTA(IY,IX).LT.0 ) THEN
MAPCXY(NCENTC) = ISEA
NCENTC = NCENTC - 1
ELSE IF ( MAPSTA(IY0,IX0).GE.1 .AND. &
MAPSTA(IY0,IX ).GE.1 .AND. &
MAPSTA(IY0,IX2).GE.1 .AND. &
MAPSTA(IY ,IX0).GE.1 .AND. &
MAPSTA(IY ,IX2).GE.1 .AND. &
MAPSTA(IY2,IX0).GE.1 .AND. &
MAPSTA(IY2,IX ).GE.1 .AND. &
MAPSTA(IY2,IX2).GE.1 ) THEN
NCENT = NCENT + 1
MAPCXY(NCENT) = ISEA
ELSE
MAPCXY(NCENTC) = ISEA
NCENTC = NCENTC - 1
END IF
END DO
END IF
!
! 5. Maps for intra-spectral propagation ---------------------------- *
!
IF ( MAPTH2(1) .NE. 0 ) RETURN
!
!/T MAPTST = 0
!
! 5.a MAPTH2 and MAPBTK
!
DO IK=1, NK
DO ITH=1, NTH
ISP = ITH + (IK-1)*NTH
ISP2 = (IK+1) + (ITH-1)*(NK+2)
MAPTH2(ISP) = ISP2
!/T MAPTST(IK+1,ITH) = MAPTST(IK+1,ITH) + 1
END DO
END DO
!
!/T WRITE (NDST,9000) 'MAPTH2', ISP, 0, 0, ISP
!/T DO IK=NK+2, 1, -1
!/T WRITE (NDST,9001) (MAPTST(IK,ITH),ITH=1, NTH)
!/T END DO
!
!/T MAPTST = 0
!
! 5.b MAPWN2
!
ISP0 = 0
DO IK=1, NK-1
DO ITH=1, NTH
ISP0 = ISP0 + 1
ISP2 = (IK+1) + (ITH-1)*(NK+2)
MAPWN2(ISP0) = ISP2
!/T MAPTST(IK+1,ITH) = MAPTST(IK+1,ITH) + 1
END DO
END DO
!
DO ITH=1, NTH
ISP0 = ISP0 + 1
ISP2 = NK+1 + (ITH-1)*(NK+2)
MAPWN2(ISP0) = ISP2
!/T MAPTST(NK+1,ITH) = MAPTST(NK+1,ITH) + 2
END DO
!
DO ITH=1, NTH
ISP0 = ISP0 + 1
ISP2 = 1 + (ITH-1)*(NK+2)
MAPWN2(ISP0) = ISP2
!/T MAPTST(1,ITH) = MAPTST(1,ITH) + 4
END DO
!
!/T WRITE (NDST,9000) 'MAPWN2', NSPEC-NTH, NTH, NTH, NSPEC+NTH
!/T DO IK=NK+2, 1, -1
!/T WRITE (NDST,9001) (MAPTST(IK,ITH),ITH=1, NTH)
!/T END DO
!
RETURN
!
! Formats
!
!/T 9000 FORMAT (/' TEST W3MAP3 : TEST MAP FOR PROPAGATION'/ &
!/T ' MAP : ',A/ &
!/T ' CENTRAL : ',I6/ &
!/T ' ABOVE : ',I6/ &
!/T ' BELOW : ',I6/ &
!/T ' TOTAL : ',I6/)
!/T 9001 FORMAT (1X,130I1)
!
!/T 9010 FORMAT (' TEST W3MAP3 : COMPOSITE MAPS TH2, WN2 AND BTK')
!/T 9011 FORMAT (2X,60I2)
!/
!/ End of W3MAP3 ----------------------------------------------------- /
!/
END SUBROUTINE W3MAP3
!/ ------------------------------------------------------------------- /
SUBROUTINE W3MAPT
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 17-Dec-2004 |
!/ +-----------------------------------+
!/
!/ 10-Dec-2001 : Origination. ( version 2.14 )
!/ 17-Dec-2004 : Multiple grid version. ( version 3.06 )
!/
! 1. Purpose :
!
! Generate 'map' arrays for the ULTIMATE QUICKEST scheme to combine
! GSE alleviation with obstructions.
!
! 2. Method :
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! See module documentation.
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3WAVE Subr. W3WAVEMD Wave model routine.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!/ ------------------------------------------------------------------- /
USE W3GDATMD, ONLY: NX, NY, NSEA, MAPSF
USE W3ADATMD, ONLY: ATRNX, ATRNY, MAPTRN
!/S USE W3SERVMD, ONLY: STRACE
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
INTEGER :: ISEA, IXY
!/S INTEGER, SAVE :: IENT = 0
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'W3MAPT')
!
! 1. Map MAPTRN ----------------------------------------------------- *
!
DO ISEA=1, NSEA
IXY = MAPSF(ISEA,3)
!notes: Oct 22 2012: I changed this because it *looks* like a bug.
! I have not confirmed that it is a bug.
! Old code is given (2 lines). Only the first line is
! changed.
!old MAPTRN(IXY) = MIN( ATRNX(IXY,1) ,ATRNY(IXY,-1) , &
!old ATRNY(IXY,1), ATRNY(IXY,-1) ) .LT. TRNMIN
MAPTRN(IXY) = MIN( ATRNX(IXY,1) ,ATRNX(IXY,-1) , &
ATRNY(IXY,1), ATRNY(IXY,-1) ) .LT. TRNMIN
END DO
!
RETURN
!
! Formats
!/
!/ End of W3MAPT ----------------------------------------------------- /
!/
END SUBROUTINE W3MAPT
!/ ------------------------------------------------------------------- /
SUBROUTINE W3XYP3 ( ISP, DTG, MAPSTA, MAPFS, VQ, VGX, VGY )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 27-May-2014 |
!/ +-----------------------------------+
!/
!/ 27-Feb-2000 : Origination. ( version 2.08 )
!/ 17-Sep-2000 : Clean-up. ( version 2.13 )
!/ 10-Dec-2001 : Sub-grid obstructions. ( version 2.14 )
!/ 16-Oct-2002 : Change INTENT for ATRNRX/Y. ( version 3.00 )
!/ 26-Dec-2002 : Moving grid version. ( version 3.02 )
!/ 01-Aug-2003 : Moving grid GSE correction. ( version 3.03 )
!/ 17-Dec-2004 : Multiple grid version. ( version 3.06 )
!/ 07-Sep-2005 : Upgrade XY boundary conditions. ( version 3.08 )
!/ 09-Nov-2005 : Removing soft boundary option. ( version 3.08 )
!/ 05-Mar-2008 : Added NEC sxf90 compiler directives.
!/ (Chris Bunney, UK Met Office) ( version 3.13 )
!/ 30-Oct-2009 : Implement curvilinear grid type. ( version 3.14 )
!/ (W. E. Rogers & T. J. Campbell, NRL)
!/ 17-Aug-2010 : Add test output. ( version 3.14.5 )
!/ 30-Oct-2010 : Implement unstructured grid ( version 3.14 )
!/ 06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
!/ specify index closure for a grid. ( version 3.14 )
!/ (T. J. Campbell, NRL)
!/ 01-Jul-2013 : Adding UQ and UNO switches to chose between third
!/ and second order schemes. ( version 4.12 )
!/ 12-Sep-2013 : Add documentation for global clos. ( version 4.12 )
!/ 27-May-2014 : Adding OMPH switch. ( version 5.02 )
!/
! 1. Purpose :
!
! Propagation in phyiscal space for a given spectral component.
!
! 2. Method :
!
! Third-order ULTIMATE QUICKEST scheme with averaging.
! Curvilinear grid implementation: Fluxes are computed in index space
! and then transformed back into physical space. The diffusion term
! is handled in physical space. The averaging scheme is adapted for
! curvilinear grids by applying the appropriate local rotations and
! adjustments to interpolation weights which control the strength of
! the averaging in axial directions.
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! ISP Int. I Number of spectral bin (IK-1)*NTH+ITH
! DTG Real I Total time step.
! MAPSTA I.A. I Grid point status map.
! MAPFS I.A. I Storage map.
! VQ R.A. I/O Field to propagate.
! VGX/Y Real I Speed of grid.
! ----------------------------------------------------------------
!
! Local variables.
! ----------------------------------------------------------------
! NTLOC Int Number of local time steps.
! DTLOC Real Local propagation time step.
! VCFL0X R.A. Local courant numbers for absolute group vel.
! using local X-grid step.
! VCFL0Y R.A. Id. in Y.
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! W3QCK3 Actual propagation algorithm
!
! STRACE Service routine.
!
! 5. Called by :
!
! W3WAVE Wave model routine.
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! - Note that the ULTIMATE limiter does not guarantee non-zero
! energies.
! - The present scheme shows a strong distorsion when propaga-
! ting a field under an angle with the grid in a truly 2-D
! fashion. Propagation is therefore split along the two
! axes.
! - Two boundary treatments are available. The first uses real
! boundaries in each space. In this case waves will not
! penetrate in narrow straights under an angle with the grid.
! This behavior is improved by using a 'soft' option, in
! which the 'X' or 'Y' sweep allows for energy to go onto
! the land. This improves the above behavior, but implies
! that X-Y connenctions are required in barriers for them
! to become inpenetrable.
! - Note that unlike in W3XYP2, isotropic diffusion is never
! used for growth.
! - Curvilinear grid implementation. Variables FACX, FACY, CCOS, CSIN,
! CCURX, CCURY are not needed and have been removed. FACX is accounted
! for as approriate in this subroutine. FACX is also accounted for in
! the case of .NOT.FLCX. Since FACX is removed, there is now a check for
! .NOT.FLCX in this subroutine. In CFL calcs dx and dy are omitted,
! since dx=dy=1 in index space. Curvilinear grid derivatives
! (DPDY, DQDX, etc.) and metric (GSQRT) are brought in via W3GDATMD.
! - The strength of the averaging scheme is dependent on grid resolution.
! Since grid resolution is non-uniform for curvilinear grids, this means
! that the strength of the averaging is also non-uniform. This may not be
! a desirable effect. A potential future upgrade would be to add an
! additional term/factor that balances the effect of the spatial
! variation of grid resolution.
!
! 8. Structure :
!
! ---------------------------------------------
! 1. Preparations
! a Set constants
! b Initialize arrays
! 2. Prepare arrays
! a Velocities and 'Q'
! 3. Loop over sub-steps
! ----------------------------------------
! a Average
! b Propagate
! c Update boundaries.
! ----------------------------------------
! 4. Store Q field in spectra
! ---------------------------------------------
!
! 9. Switches :
!
! !/NEC Enable NEC SXF90 compiler directives.
!
! !/S Enable subroutine tracing.
!
! !/OMPH Hybrid OpenMP directives.
!
! !/MGP Moving grid corrections.
! !/MGG Moving grid corrections.
!
! !/T Enable general test output.
! !/T0 Dump of precalcaulted interpolation data.
! !/T1 Dump of input field and fluxes.
! !/T2 Dump of output field (before boundary update).
! !/T3 Dump of output field (final).
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE CONSTANTS
!
USE W3TIMEMD, ONLY: DSEC21
!
USE W3GDATMD, ONLY: NX, NY, NSEA, MAPSF, DTCFL, CLATS, &
ICLOSE, FLCX, FLCY, NK, NTH, DTH, XFR, &
ICLOSE_NONE, ICLOSE_SMPL, ICLOSE_TRPL, &
ECOS, ESIN, SIG, WDCG, WDTH, PFMOVE, &
FLAGLL, DPDX, DPDY, DQDX, DQDY, GSQRT
USE W3WDATMD, ONLY: TIME
USE W3ADATMD, ONLY: NMX0, NMX1, NMX2, NMY0, NMY1, NMY2, NACT, &
NCENT, MAPX2, MAPY2, MAPAXY, MAPCXY, &
MAPTRN, CG, CX, CY, ATRNX, ATRNY, ITIME
USE W3IDATMD, ONLY: FLCUR
USE W3ODATMD, ONLY: NDSE, NDST, FLBPI, NBI, TBPI0, TBPIN, &
ISBPI, BBPI0, BBPIN, IAPROC, NAPERR
USE W3SERVMD, ONLY: EXTCDE
!/S USE W3SERVMD, ONLY: STRACE
!/UQ USE W3UQCKMD
!/UNO USE W3UNO2MD
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
INTEGER, INTENT(IN) :: ISP, MAPSTA(NY*NX), MAPFS(NY*NX)
REAL, INTENT(IN) :: DTG, VGX, VGY
REAL, INTENT(INOUT) :: VQ(1-NY:NY*(NX+2))
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
INTEGER :: ITH, IK, NTLOC, ITLOC, ISEA, IXY, IP
INTEGER :: IX, IY, IXC, IYC, IBI
INTEGER :: IIXY1(NSEA), IIXY2(NSEA), &
IIXY3(NSEA), IIXY4(NSEA)
!/S INTEGER, SAVE :: IENT = 0
REAL :: CG0, CGA, CGN, CGX, CGY, CXC, CYC, &
CXMIN, CXMAX, CYMIN, CYMAX
REAL :: CGC, FGSE = 1.
REAL :: FTH, FTHX, FTHY, FCG, FCGX, FCGY
REAL :: DTLOC, DTRAD, &
DXCGN, DYCGN, DXCGS, DYCGS, DXCGC, &
DYCGC
REAL :: RDI1(NSEA), RDI2(NSEA), &
RDI3(NSEA), RDI4(NSEA)
REAL :: TMPX, TMPY, RD1, RD2, RD3, RD4
LOGICAL :: YFIRST
LOGICAL :: GLOBAL
REAL :: CP, CQ
!/
!/ Automatic work arrays
!/
INTEGER :: MAPSTX(1-2*NY:NY*(NX+2))
REAL :: VLCFLX((NX+1)*NY), VLCFLY((NX+1)*NY),&
AQ(1-NY:NY*(NX+2))
REAL :: CXTOT((NX+1)*NY), CYTOT(NX*NY)
!/
!/ ------------------------------------------------------------------- /
!/
!/S CALL STRACE (IENT, 'W3XYP3')
!
! 1. Preparations --------------------------------------------------- *
IF ( ICLOSE .EQ. ICLOSE_TRPL ) THEN
!/ ------------------------------------------------------------------- /
IF (IAPROC .EQ. NAPERR) &
WRITE(NDSE,*)'SUBROUTINE W3XYP3 IS NOT YET ADAPTED FOR '// &
'TRIPOLE GRIDS. STOPPING NOW.'
CALL EXTCDE ( 1 )
END IF
! 1.a Set constants
!
GLOBAL = ICLOSE.NE.ICLOSE_NONE
ITH = 1 + MOD(ISP-1,NTH)
IK = 1 + (ISP-1)/NTH
!
CG0 = 0.575 * GRAV / SIG(1)
CGA = 0.575 * GRAV / SIG(IK)
CGX = CGA * ECOS(ITH)
CGY = CGA * ESIN(ITH)
!/MGP CGX = CGX - VGX
!/MGP CGY = CGY - VGY
CGC = SQRT ( CGX**2 + CGY**2 )
!/MGG FGSE = ( CGA / MAX(0.001*CGA,CGC) )**PFMOVE
!
IF ( FLCUR ) THEN
CXMIN = MINVAL ( CX(1:NSEA) )
CXMAX = MAXVAL ( CX(1:NSEA) )
CYMIN = MINVAL ( CY(1:NSEA) )
CYMAX = MAXVAL ( CY(1:NSEA) )
IF ( ABS(CGX+CXMIN) .GT. ABS(CGX+CXMAX) ) THEN
CGX = CGX + CXMIN
ELSE
CGX = CGX + CXMAX
END IF
IF ( ABS(CGY+CYMIN) .GT. ABS(CGY+CYMAX) ) THEN
CGY = CGY + CYMIN
ELSE
CGY = CGY + CYMAX
END IF
CXC = MAX ( ABS(CXMIN) , ABS(CXMAX) )
CYC = MAX ( ABS(CYMIN) , ABS(CYMAX) )
!/MGP CXC = MAX ( ABS(CXMIN-VGX) , ABS(CXMAX-VGX) )
!/MGP CYC = MAX ( ABS(CYMIN-VGY) , ABS(CYMAX-VGY) )
ELSE
CXC = 0.
CYC = 0.
END IF
!
CGN = MAX ( ABS(CGX) , ABS(CGY) , CXC, CYC, 0.001*CG0 )
!
NTLOC = 1 + INT(DTG/(DTCFL*CG0/CGN))
DTLOC = DTG / REAL(NTLOC)
DTRAD = DTLOC
IF ( FLAGLL ) DTRAD=DTRAD/(DERA*RADIUS)
!
YFIRST = MOD(ITIME,2) .EQ. 0
!
!/T WRITE (NDST,9000) YFIRST
!/T WRITE (NDST,9001) ISP, ITH, IK, ECOS(ITH), ESIN(ITH)
!
! 1.b Initialize arrays
!
!/T WRITE (NDST,9010)
!
VLCFLX = 0.
VLCFLY = 0.
CXTOT = 0.
CYTOT = 0.
!
MAPSTX(1:NX*NY) = MAPSTA(1:NX*NY)
!
IF ( GLOBAL ) THEN
MAPSTX(1-2*NY:0) = MAPSTA((NX-2)*NY+1:NX*NY)
MAPSTX(NX*NY+1:NX*NY+2*NY) = MAPSTA(1:2*NY)
ELSE
MAPSTX(1-2*NY:0) = 0
MAPSTX(NX*NY+1:NX*NY+2*NY) = 0
END IF
!
! 1.c Pre-calculate interpolation info
!
FTH = FGSE * WDTH * DTH * DTLOC
FCG = FGSE * WDCG * 0.5 * (XFR-1./XFR) * DTLOC
IF ( FLAGLL ) THEN
FTH = FTH / RADIUS / DERA
FCG = FCG / RADIUS / DERA
END IF
FCG = FCG / REAL(NTLOC) !TJC: required to match original (is this correct?)
FTHX = - FTH * ESIN(ITH)
FTHY = FTH * ECOS(ITH)
FCGX = FCG * ECOS(ITH)
FCGY = FCG * ESIN(ITH)
!
!/T0 WRITE (NDST,9011)
!
!/OMPH/!$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, TMPX, TMPY, &
!/OMPH/!$OMP& DXCGN, DYCGN, DXCGS, DYCGS, DXCGC, DYCGC, &
!/OMPH/!$OMP& IXC, IYC)
!
DO ISEA=1, NSEA
!
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
!
! 1.c.1 Normal and parallel width ...
!
TMPX = FTHX / CLATS(ISEA)
TMPY = FTHY
DXCGN = DPDX(IY,IX)*TMPX + DPDY(IY,IX)*TMPY
DYCGN = DQDX(IY,IX)*TMPX + DQDY(IY,IX)*TMPY
TMPX = FCGX / CLATS(ISEA)
TMPY = FCGY
DXCGS = DPDX(IY,IX)*TMPX + DPDY(IY,IX)*TMPY
DYCGS = DQDX(IY,IX)*TMPX + DQDY(IY,IX)*TMPY
!
! 1.c.2 "Sum" corner (and mirror image) ...
!
DXCGC = DXCGN + DXCGS
DYCGC = DYCGN + DYCGS
!
IXC = NY
IF ( DXCGC .LT. 0. ) IXC = - IXC
IYC = 1
IF ( DYCGC .LT. 0. ) IYC = - IYC
!
IIXY1(ISEA) = IXC + IYC
IF ( ABS(DXCGC) .GT. ABS(DYCGC) ) THEN
IIXY2(ISEA) = IXC
RDI1 (ISEA) = ABS(DYCGC/DXCGC)
RDI2 (ISEA) = ABS(DXCGC)
ELSE
IIXY2(ISEA) = IYC
IF ( ABS(DYCGC) .GT. 1.E-5 ) THEN
RDI1(ISEA) = ABS(DXCGC/DYCGC)
ELSE
RDI1(ISEA) = 1.
END IF
RDI2(ISEA) = ABS(DYCGC)
END IF
!
!/T0 WRITE (NDST,9012) ISEA, ITH, IIXY1(ISEA), IIXY2(ISEA), &
!/T0 RDI1(ISEA), RDI2(ISEA)*CG(IK,1)
!
! 1.c.2 "Difference" corner (and mirror image) ...
!
DXCGC = DXCGN - DXCGS
DYCGC = DYCGN - DYCGS
!
IXC = NY
IF ( DXCGC .LT. 0. ) IXC = - IXC
IYC = 1
IF ( DYCGC .LT. 0. ) IYC = - IYC
!
IIXY3(ISEA) = IXC + IYC
IF ( ABS(DXCGC) .GT. ABS(DYCGC) ) THEN
IIXY4(ISEA) = IXC
RDI3 (ISEA) = ABS(DYCGC/DXCGC)
RDI4 (ISEA) = ABS(DXCGC)
ELSE
IIXY4(ISEA) = IYC
IF ( ABS(DYCGC) .GT. 1.E-5 ) THEN
RDI3(ISEA) = ABS(DXCGC/DYCGC)
ELSE
RDI3(ISEA) = 1.
END IF
RDI4(ISEA) = ABS(DYCGC)
END IF
!
!/T0 WRITE (NDST,9013) IIXY3(ISEA), IIXY4(ISEA), RDI3(ISEA), &
!/T0 RDI4(ISEA)*CG(IK,1)
!
END DO
!
!/OMPH/!$OMP END PARALLEL DO
!
! 2. Calculate velocities and diffusion coefficients ---------------- *
! 2.a Velocities
!
! Q = ( A / CG * CLATS )
! LCFLX = ( COS*CG / CLATS ) * DT / DX
! LCFLY = ( SIN*CG ) * DT / DY
!
!/T WRITE (NDST,9020) NSEA
!
!/OMPH/!$OMP PARALLEL DO PRIVATE (ISEA, IXY)
!
DO ISEA=1, NSEA
IXY = MAPSF(ISEA,3)
VQ (IXY) = VQ(IXY) / CG(IK,ISEA) * CLATS(ISEA)
CXTOT(IXY) = ECOS(ITH) * CG(IK,ISEA) / CLATS(ISEA)
CYTOT(IXY) = ESIN(ITH) * CG(IK,ISEA)
!/MGP CXTOT(IXY) = CXTOT(IXY) - VGX/CLATS(ISEA)
!/MGP CYTOT(IXY) = CYTOT(IXY) - VGY
!/T1 IF ( .NOT. FLCUR ) &
!/T1 WRITE (NDST,9021) ISEA, MAPSF(ISEA,1), MAPSF(ISEA,2), &
!/T1 VQ(IXY), CXTOT(IXY), CYTOT(IXY)
END DO
!
!/OMPH/!$OMP END PARALLEL DO
!
IF ( FLCUR ) THEN
!/T WRITE (NDST,9022)
!
!/OMPH/!$OMP PARALLEL DO PRIVATE (IXY)
!
DO ISEA=1, NSEA
IXY = MAPSF(ISEA,3)
CXTOT(IXY) = CXTOT(IXY) + CX(ISEA)/CLATS(ISEA)
CYTOT(IXY) = CYTOT(IXY) + CY(ISEA)
!/T1 WRITE (NDST,9021) ISEA, MAPSF(ISEA,1), MAPSF(ISEA,2), &
!/T1 VQ(IXY), CXTOT(IXY), CYTOT(IXY)
END DO
!
!/OMPH/!$OMP END PARALLEL DO
!
END IF
!
!/OMPH/!$OMP PARALLEL DO PRIVATE (IX, IY, IXY, CP, CQ)
!
DO ISEA=1, NSEA
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
IXY = MAPSF(ISEA,3)
CP = CXTOT(IXY)*DPDX(IY,IX) + CYTOT(IXY)*DPDY(IY,IX)
CQ = CXTOT(IXY)*DQDX(IY,IX) + CYTOT(IXY)*DQDY(IY,IX)
VLCFLX(IXY) = CP*DTRAD
VLCFLY(IXY) = CQ*DTRAD
END DO
!
!/OMPH/!$OMP END PARALLEL DO
!
! 3. Loop over sub-steps -------------------------------------------- *
!
DO ITLOC=1, NTLOC
!
! 3.a Average
!
AQ = VQ
VQ = 0.
!
! 3.a.1 Central points