-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathHydroCalc.f90
5203 lines (3431 loc) · 293 KB
/
HydroCalc.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
MODULE Waves
! This MODULE stores variables and routines associated with incident
! waves and current.
USE NWTC_Library
IMPLICIT NONE
COMPLEX(ReKi), PARAMETER :: ImagNmbr = (0.0,1.0) ! The imaginary number, SQRT(-1.0)
COMPLEX(ReKi), ALLOCATABLE :: WaveElevC0(:) ! Fourier transform of the instantaneous elevation of incident waves at the platform reference point (m-s)
REAL(ReKi), ALLOCATABLE :: DZNodes (:) ! Length of variable-length tower or platform elements for the points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed (meters)
REAL(ReKi) :: Gravity ! Gravitational acceleration (m/s^2)
REAL(ReKi), PARAMETER :: Inv2Pi = 0.15915494 ! 0.5/Pi.
REAL(ReKi), PARAMETER :: PiOvr4 = 0.78539816 ! Pi/4.
REAL(ReKi) :: RhoXg ! = WtrDens*Gravity
REAL(ReKi), ALLOCATABLE :: WaveAcc0 (:,:,:) ! Instantaneous acceleration of incident waves in the xi- (1), yi- (2), and zi- (3) directions, respectively, accounting for stretching, at each of the NWaveKin0 points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed (m/s^2)
REAL(ReKi) :: WaveDir ! Incident wave propagation heading direction (degrees)
REAL(ReKi) :: WaveDOmega ! Frequency step for incident wave calculations (rad/s)
REAL(ReKi), ALLOCATABLE :: WaveElev (:,:) ! Instantaneous elevation of incident waves at each of the NWaveElev points where the incident wave elevations can be output (meters)
REAL(ReKi), ALLOCATABLE :: WaveElev0 (:) ! Instantaneous elevation of incident waves at the platform reference point (meters)
REAL(ReKi), ALLOCATABLE :: WaveKinzi0(:) ! zi-coordinates for points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed; these are relative to the mean see level (meters)
REAL(ReKi), ALLOCATABLE :: WaveTime (:) ! Simulation times at which the instantaneous elevation of, velocity of, acceleration of, and loads associated with the incident waves are determined (sec)
REAL(ReKi), ALLOCATABLE :: WaveVel0 (:,:,:) ! Instantaneous velocity of incident waves in the xi- (1), yi- (2), and zi- (3) directions, respectively, accounting for stretching, at each of the NWaveKin0 points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed (m/s ) (The values include both the velocity of incident waves and the velocity of current.)
REAL(ReKi) :: WtrDens ! Water density (kg/m^3)
REAL(ReKi) :: WtrDpth ! Water depth (meters)
INTEGER(4) :: NStepWave ! Total number of frequency components = total number of time steps in the incident wave (-)
INTEGER(4) :: NStepWave2 ! NStepWave/2
INTEGER(4) :: NWaveElev ! Number of points where the incident wave elevation can be output (-)
INTEGER(4) :: NWaveKin0 ! Number of points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed (-)
INTEGER(4) :: WaveMod ! Incident wave kinematics model {0: none=still water, 1: plane progressive (regular), 2: JONSWAP/Pierson-Moskowitz spectrum (irregular), 3: user-defind spectrum from routine UserWaveSpctrm (irregular)}
INTEGER(4) :: WaveStMod ! Model for stretching incident wave kinematics to instantaneous free surface {0: none=no stretching, 1: vertical stretching, 2: extrapolation stretching, 3: Wheeler stretching}
CHARACTER(1024) :: DirRoot ! The name of the root file including the full path to the current working directory. This may be useful if you want this routine to write a permanent record of what it does to be stored with the simulation results: the results should be stored in a file whose name (including path) is generated by appending any suitable extension to DirRoot.
CONTAINS
!=======================================================================
SUBROUTINE InitWaves ( WtrDensIn , WtrDpthIn , WaveModIn, WaveStModIn, &
WaveTMaxIn , WaveDT , WaveHs , WaveTp , &
WavePkShp , WaveDirIn , WaveSeed , GHWvFile , &
CurrMod , CurrSSV0 , CurrSSDir, CurrNSRef , &
CurrNSV0 , CurrNSDir , CurrDIV , CurrDIDir , &
NWaveKin0In, WaveKinzi0In, DZNodesIn, NWaveElevIn, &
WaveElevxi , WaveElevyi , GravityIn, DirRootIn )
! This routine is used to initialize the variables associated with
! incident waves and current.
USE FFT_Module
IMPLICIT NONE
! Passed Variables:
INTEGER(4), INTENT(IN ) :: NWaveElevIn ! Number of points where the incident wave elevations can be output (-)
INTEGER(4), INTENT(IN ) :: NWaveKin0In ! Number of points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed (-)
REAL(ReKi), INTENT(IN ) :: CurrDIDir ! Depth-independent current heading direction (degrees)
REAL(ReKi), INTENT(IN ) :: CurrDIV ! Depth-independent current velocity (m/s)
REAL(ReKi), INTENT(IN ) :: CurrNSDir ! Near-surface current heading direction (degrees)
REAL(ReKi), INTENT(IN ) :: CurrNSRef ! Near-surface current reference depth (meters)
REAL(ReKi), INTENT(IN ) :: CurrNSV0 ! Near-surface current velocity at still water level (m/s)
REAL(ReKi), INTENT(IN ) :: CurrSSDir ! Sub-surface current heading direction (degrees)
REAL(ReKi), INTENT(IN ) :: CurrSSV0 ! Sub-surface current velocity at still water level (m/s)
REAL(ReKi), INTENT(IN ) :: DZNodesIn (NWaveKin0In) ! Length of variable-length tower or platform elements for the points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed (meters)
REAL(ReKi), INTENT(IN ) :: GravityIn ! Gravitational acceleration (m/s^2)
REAL(ReKi), INTENT(IN ) :: WaveDirIn ! Incident wave propagation heading direction (degrees)
REAL(ReKi), INTENT(IN ) :: WaveDT ! Time step for incident wave calculations (sec)
REAL(ReKi), INTENT(IN ) :: WaveElevxi (NWaveElevIn) ! xi-coordinates for points where the incident wave elevations can be output (meters)
REAL(ReKi), INTENT(IN ) :: WaveElevyi (NWaveElevIn) ! yi-coordinates for points where the incident wave elevations can be output (meters)
REAL(ReKi), INTENT(IN ) :: WaveKinzi0In(NWaveKin0In) ! zi-coordinates for points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed; these are relative to the mean see level (meters)
REAL(ReKi), INTENT(IN ) :: WaveHs ! Significant wave height of incident waves (meters)
REAL(ReKi), INTENT(IN ) :: WavePkShp ! Peak shape parameter of incident wave spectrum [1.0 for Pierson-Moskowitz] (-)
REAL(ReKi), INTENT(IN ) :: WaveTMaxIn ! Analysis time for incident wave calculations; the actual analysis time may be larger than this value in order for the maintain an effecient FFT (sec)
REAL(ReKi), INTENT(IN ) :: WaveTp ! Peak spectral period of incident waves (sec)
REAL(ReKi), INTENT(IN ) :: WtrDensIn ! Water density (kg/m^3)
REAL(ReKi), INTENT(IN ) :: WtrDpthIn ! Water depth (meters)
INTEGER(4), INTENT(IN ) :: CurrMod ! Current profile model {0: none=no current, 1: standard, 2: user-defined from routine UserCurrent}
INTEGER(4), INTENT(IN ) :: WaveStModIn ! Model for stretching incident wave kinematics to instantaneous free surface {0: none=no stretching, 1: vertical stretching, 2: extrapolation stretching, 3: Wheeler stretching}
INTEGER(4), INTENT(IN ) :: WaveModIn ! Incident wave kinematics model {0: none=still water, 1: plane progressive (regular), 2: JONSWAP/Pierson-Moskowitz spectrum (irregular), 3: user-defind spectrum from routine UserWaveSpctrm (irregular)}
INTEGER(4), INTENT(IN ) :: WaveSeed (2) ! Random seeds of incident waves [-2147483648 to 2147483647]
CHARACTER(1024), INTENT(IN ) :: DirRootIn ! The name of the root file including the full path to the current working directory. This may be useful if you want this routine to write a permanent record of what it does to be stored with the simulation results: the results should be stored in a file whose name (including path) is generated by appending any suitable extension to DirRoot.
CHARACTER(1024), INTENT(IN ) :: GHWvFile ! The root name of GH Bladed files containing wave data.
! Local Variables:
COMPLEX(ReKi) :: ImagOmega ! = ImagNmbr*Omega (rad/s)
COMPLEX(ReKi), ALLOCATABLE :: PWaveAccC0HPz0 (:) ! Partial derivative of WaveAccC0H(:) with respect to zi at zi = 0 (1/s)
COMPLEX(ReKi), ALLOCATABLE :: PWaveAccC0VPz0 (:) ! Partial derivative of WaveAccC0V(:) with respect to zi at zi = 0 (1/s)
COMPLEX(ReKi), ALLOCATABLE :: PWaveVelC0HPz0 (:) ! Partial derivative of WaveVelC0H(:) with respect to zi at zi = 0 (- )
COMPLEX(ReKi), ALLOCATABLE :: PWaveVelC0VPz0 (:) ! Partial derivative of WaveVelC0V(:) with respect to zi at zi = 0 (- )
COMPLEX(ReKi), ALLOCATABLE :: WaveAccC0H(:,:) ! Fourier transform of the instantaneous horizontal acceleration of incident waves before applying stretching at the zi-coordinates for points along a vertical line passing through the platform reference point (m/s)
COMPLEX(ReKi), ALLOCATABLE :: WaveAccC0V(:,:) ! Fourier transform of the instantaneous vertical acceleration of incident waves before applying stretching at the zi-coordinates for points along a vertical line passing through the platform reference point (m/s)
COMPLEX(ReKi), ALLOCATABLE :: WaveElevC (:,:) ! Fourier transform of the instantaneous elevation of incident waves (m-s)
COMPLEX(ReKi), ALLOCATABLE :: WaveVelC0H(:,:) ! Fourier transform of the instantaneous horizontal velocity of incident waves before applying stretching at the zi-coordinates for points along a vertical line passing through the platform reference point (meters)
COMPLEX(ReKi), ALLOCATABLE :: WaveVelC0V(:,:) ! Fourier transform of the instantaneous vertical velocity of incident waves before applying stretching at the zi-coordinates for points along a vertical line passing through the platform reference point (meters)
COMPLEX(ReKi) :: WGNC ! Fourier transform of the realization of a White Gaussian Noise (WGN) time series process with unit variance for the current frequency component (SQRT(sec))
REAL(ReKi) :: CurrVxi ! xi-component of the current velocity at the instantaneous elevation (m/s)
REAL(ReKi) :: CurrVyi ! yi-component of the current velocity at the instantaneous elevation (m/s)
REAL(ReKi) :: CurrVxi0 ! xi-component of the current velocity at zi = 0.0 meters (m/s)
REAL(ReKi) :: CurrVyi0 ! yi-component of the current velocity at zi = 0.0 meters (m/s)
REAL(ReKi) :: CurrVxiS ! xi-component of the current velocity at zi = -SmllNmbr meters (m/s)
REAL(ReKi) :: CurrVyiS ! yi-component of the current velocity at zi = -SmllNmbr meters (m/s)
REAL(ReKi) :: CWaveDir ! COS( WaveDir )
REAL(ReKi) :: GHQBar ! Unused instantaneous dynamic pressure of incident waves in GH Bladed wave data files (N/m^2)
REAL(ReKi), ALLOCATABLE :: GHWaveAcc (:,:) ! Instantaneous acceleration of incident waves in the xi- (1), yi- (2), and zi- (3) directions, respectively, at each of the GHNWvDpth vertical locations in GH Bladed wave data files (m/s^2)
REAL(ReKi) :: GHWaveTime ! Instantaneous simulation times in GH Bladed wave data files (sec)
REAL(ReKi), ALLOCATABLE :: GHWaveVel (:,:) ! Instantaneous velocity of incident waves in the xi- (1), yi- (2), and zi- (3) directions, respectively, at each of the GHNWvDpth vertical locations in GH Bladed wave data files (m/s )
REAL(ReKi), ALLOCATABLE :: GHWvDpth (:) ! Vertical locations in GH Bladed wave data files.
REAL(ReKi), PARAMETER :: n_Massel = 3.0 ! Factor used to the scale the peak spectral frequency in order to find the cut-off frequency based on the suggestion in: Massel, S. R., Ocean Surface Waves: Their Physics and Prediction, Advanced Series on Ocean Engineering - Vol. 11, World Scientific Publishing, Singapore - New Jersey - London - Hong Kong, 1996. This reference recommends n_Massel > 3.0 (higher for higher-order wave kinemetics); the ">" designation is accounted for by checking if ( Omega > OmegaCutOff ).
REAL(ReKi) :: Omega ! Wave frequency (rad/s)
REAL(ReKi) :: OmegaCutOff ! Cut-off frequency or upper frequency limit of the wave spectrum beyond which the wave spectrum is zeroed (rad/s)
REAL(ReKi) :: PCurrVxiPz0 ! Partial derivative of CurrVxi with respect to zi at zi = 0 (1/s )
REAL(ReKi) :: PCurrVyiPz0 ! Partial derivative of CurrVyi with respect to zi at zi = 0 (1/s )
REAL(ReKi), ALLOCATABLE :: PWaveAcc0HPz0 (:) ! Partial derivative of WaveAcc0H (:) with respect to zi at zi = 0 (1/s^2)
REAL(ReKi), ALLOCATABLE :: PWaveAcc0VPz0 (:) ! Partial derivative of WaveAcc0V (:) with respect to zi at zi = 0 (1/s^2)
REAL(ReKi), ALLOCATABLE :: PWaveVel0HPz0 (:) ! Partial derivative of WaveVel0H (:) with respect to zi at zi = 0 (1/s )
REAL(ReKi), ALLOCATABLE :: PWaveVel0HxiPz0(:) ! Partial derivative of WaveVel0Hxi(:) with respect to zi at zi = 0 (1/s )
REAL(ReKi), ALLOCATABLE :: PWaveVel0HyiPz0(:) ! Partial derivative of WaveVel0Hyi(:) with respect to zi at zi = 0 (1/s )
REAL(ReKi), ALLOCATABLE :: PWaveVel0VPz0 (:) ! Partial derivative of WaveVel0V (:) with respect to zi at zi = 0 (1/s )
REAL(ReKi) :: S2Sd_Fact ! Factor used to scale the magnitude of the WaveS2Sdd as required by the discrete time inverse Fourier transform (-)
REAL(ReKi) :: Slope ! Miscellanous slope used in an interpolation (-)
REAL(ReKi), PARAMETER :: SmllNmbr = 9.999E-4 ! A small number representing epsilon for taking numerical derivatives.
REAL(ReKi) :: SWaveDir ! SIN( WaveDir )
REAL(ReKi), ALLOCATABLE :: WaveAcc0H (:,:) ! Instantaneous horizontal acceleration of incident waves before applying stretching at the zi-coordinates for points along a vertical line passing through the platform reference point (m/s^2)
REAL(ReKi) :: WaveAcc0HExtrap ! Temporary value extrapolated from the WaveAcc0H (:,:) array (m/s^2)
REAL(ReKi) :: WaveAcc0HInterp ! Temporary value interpolated from the WaveAcc0H (:,:) array (m/s^2)
REAL(ReKi), ALLOCATABLE :: WaveAcc0V (:,:) ! Instantaneous vertical acceleration of incident waves before applying stretching at the zi-coordinates for points along a vertical line passing through the platform reference point (m/s^2)
REAL(ReKi) :: WaveAcc0VExtrap ! Temporary value extrapolated from the WaveAcc0V (:,:) array (m/s^2)
REAL(ReKi) :: WaveAcc0VInterp ! Temporary value interpolated from the WaveAcc0V (:,:) array (m/s^2)
REAL(ReKi) :: WaveElev_Max ! Maximum expected value of the instantaneous elevation of incident waves (meters)
REAL(ReKi) :: WaveElev_Min ! Minimum expected value of the instantaneous elevation of incident waves (meters)
REAL(ReKi), ALLOCATABLE :: WaveElevxiPrime(:) ! Locations along the wave heading direction for points where the incident wave elevations can be output (meters)
REAL(ReKi), ALLOCATABLE :: WaveKinzi0Prime(:) ! zi-coordinates for points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed before applying stretching; these are relative to the mean see level (meters)
REAL(ReKi), ALLOCATABLE :: WaveKinzi0St (:) ! Array of elevations found by stretching the elevations in the WaveKinzi0Prime(:) array using the instantaneous wave elevation; these are relative to the mean see level (meters)
REAL(ReKi) :: WaveNmbr ! Wavenumber of the current frequency component (1/meter)
REAL(ReKi) :: WaveS1Sdd ! One-sided power spectral density of the wave spectrum per unit time for the current frequency component (m^2/(rad/s))
REAL(ReKi) :: WaveS2Sdd ! Two-sided power spectral density of the wave spectrum per unit time for the current frequency component (m^2/(rad/s))
REAL(ReKi) :: WaveTMax ! Analysis time for incident wave calculations (sec)
REAL(ReKi), ALLOCATABLE :: WaveVel0H (:,:) ! Instantaneous horizontal velocity of incident waves before applying stretching at the zi-coordinates for points along a vertical line passing through the platform reference point (m/s )
REAL(ReKi), ALLOCATABLE :: WaveVel0Hxi (:,:) ! Instantaneous xi-direction velocity of incident waves before applying stretching at the zi-coordinates for points along a vertical line passing through the platform reference point (m/s )
REAL(ReKi) :: WaveVel0HxiExtrap ! Temporary value extrapolated from the WaveVel0Hxi(:,:) array (m/s )
REAL(ReKi) :: WaveVel0HxiInterp ! Temporary value interpolated from the WaveVel0Hxi(:,:) array (m/s )
REAL(ReKi), ALLOCATABLE :: WaveVel0Hyi (:,:) ! Instantaneous yi-direction velocity of incident waves before applying stretching at the zi-coordinates for points along a vertical line passing through the platform reference point (m/s )
REAL(ReKi) :: WaveVel0HyiExtrap ! Temporary value extrapolated from the WaveVel0Hyi(:,:) array (m/s )
REAL(ReKi) :: WaveVel0HyiInterp ! Temporary value interpolated from the WaveVel0Hyi(:,:) array (m/s )
REAL(ReKi), ALLOCATABLE :: WaveVel0V (:,:) ! Instantaneous vertical velocity of incident waves before applying stretching at the zi-coordinates for points along a vertical line passing through the platform reference point (m/s )
REAL(ReKi) :: WaveVel0VExtrap ! Temporary value extrapolated from the WaveVel0V (:,:) array (m/s )
REAL(ReKi) :: WaveVel0VInterp ! Temporary value interpolated from the WaveVel0V (:,:) array (m/s )
REAL(ReKi) :: WGNC_Fact ! Factor used to scale the magnitude of the WGNC as required by the discrete time inverse Fourier transform (-)
REAL(ReKi) :: zi_Max ! Maximum elevation where the wave kinematics are to be applied using stretching to the instantaneous free surface (meters)
REAL(ReKi) :: zi_Min ! Minimum elevation where the wave kinematics are to be applied using stretching to the instantaneous free surface (meters)
REAL(ReKi) :: ziPrime_Max ! Maximum elevation where the wave kinematics are computed before applying stretching to the instantaneous free surface (meters)
REAL(ReKi) :: ziPrime_Min ! Minimum elevation where the wave kinematics are computed before applying stretching to the instantaneous free surface (meters)
INTEGER(4) :: GHNStepWave ! Total number of time steps in the GH Bladed wave data files.
INTEGER(4) :: GHNWvDpth ! Number of vertical locations in GH Bladed wave data files.
INTEGER(4) :: I ! Generic index
INTEGER(4) :: I_Orig ! The index of the time step from original (input) part of data
INTEGER(4) :: I_WaveTp ! The index of the frequency component nearest to WaveTp
INTEGER(4) :: J ! Generic index
INTEGER(4) :: J_Min ! The minimum value of index J such that WaveKinzi0(J) >= -WtrDpth
INTEGER(4) :: K ! Generic index
INTEGER(4) :: LastInd = 1 ! Index into the arrays saved from the last call as a starting point for this call
INTEGER :: nSeeds ! number of seeds required to initialize the intrinsic random number generator
INTEGER(4) :: NWaveKin0Prime ! Number of points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed before applying stretching to the instantaneous free surface (-)
INTEGER(4) :: Sttus ! Status returned by an attempted allocation or READ.
INTEGER, ALLOCATABLE :: TmpWaveSeeds (:) ! A temporary array used for portability. IVF/CVF use a random number generator initialized with 2 seeds; other platforms can use different implementations (e.g. gfortran needs 8 or 12 seeds)
INTEGER(4) :: UnFA = 31 ! I/O unit number for the file needed for the GH Bladed wave data by FAST.
INTEGER(4) :: UnKi = 32 ! I/O unit number for the GH Bladed wave data file containing wave particle kinematics time history.
INTEGER(4) :: UnSu = 33 ! I/O unit number for the GH Bladed wave data file containing surface elevation time history.
LOGICAL :: Reading = .TRUE. ! Flag to say whether or not we are still reading from the GH Bladed wave data files (files not exhausted).
! Save these values for future use:
NWaveElev = NWaveElevIn
NWaveKin0 = NWaveKin0In
ALLOCATE ( DZNodes (NWaveKin0) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the DZNodes array.')
ENDIF
ALLOCATE ( WaveKinzi0(NWaveKin0) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveKinzi0 array.')
ENDIF
DirRoot = DirRootIn
DZNodes (:) = DZNodesIn (:)
Gravity = GravityIn
WaveDir = WaveDirIn
WaveKinzi0(:) = WaveKinzi0In(:)
WaveMod = WaveModIn
WaveStMod = WaveStModIn
WaveTMax = WaveTMaxIn
WtrDens = WtrDensIn
WtrDpth = WtrDpthIn
RhoXg = WtrDens*Gravity
! Initialize the variables associated with the incident wave:
SELECT CASE ( WaveMod ) ! Which incident wave kinematics model are we using?
CASE ( 0 ) ! None=still water.
! Initialize everything to zero:
NStepWave = 2 ! We must have at least two elements in order to interpolate later on
NStepWave2 = 1
ALLOCATE ( WaveTime (0:NStepWave-1 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveTime array.')
ENDIF
ALLOCATE ( WaveElevC0 (0:NStepWave2 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveElevC0 array.')
ENDIF
ALLOCATE ( WaveElev0 (0:NStepWave-1 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveElev0 array.')
ENDIF
ALLOCATE ( WaveElev (0:NStepWave-1,NWaveElev ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveElev array.')
ENDIF
ALLOCATE ( WaveVel0 (0:NStepWave-1,NWaveKin0,3) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveVel0 array.')
ENDIF
ALLOCATE ( WaveAcc0 (0:NStepWave-1,NWaveKin0,3) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveAcc0 array.')
ENDIF
WaveDOmega = 0.0
WaveTime = (/ 0.0, 1.0 /) ! We must have at least two different time steps in the interpolation
WaveElevC0 = (0.0,0.0)
WaveElev0 = 0.0
WaveElev = 0.0
WaveVel0 = 0.0
WaveAcc0 = 0.0
! Add the current velocities to the wave velocities:
DO J = 1,NWaveKin0 ! Loop through all points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed
CALL InitCurrent ( CurrMod , CurrSSV0 , CurrSSDir, CurrNSRef, &
CurrNSV0 , CurrNSDir, CurrDIV , CurrDIDir, &
WaveKinzi0(J), WtrDpth , DirRoot , CurrVxi , CurrVyi )
WaveVel0(:,J,1) = WaveVel0(:,J,1) + CurrVxi ! xi-direction
WaveVel0(:,J,2) = WaveVel0(:,J,2) + CurrVyi ! yi-direction
ENDDO ! J - All points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed
CASE ( 1, 2, 3 ) ! Plane progressive (regular) wave, JONSWAP/Pierson-Moskowitz spectrum (irregular) wave, or user-defined spectrum (irregular) wave.
! Tell our nice users what is about to happen that may take a while:
CALL WrScr ( ' Generating incident wave kinematics and current time history.' )
! Calculate the locations of the points along the wave heading direction
! where the incident wave elevations can be output:
CWaveDir = COS( D2R*WaveDir )
SWaveDir = SIN( D2R*WaveDir )
ALLOCATE ( WaveElevxiPrime (NWaveElev) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveElevxiPrime array.')
ENDIF
DO J = 1,NWaveElev ! Loop through all points where the incident wave elevations can be output
WaveElevxiPrime(J) = WaveElevxi(J)*CWaveDir + WaveElevyi(J)*SWaveDir
ENDDO ! J - All points where the incident wave elevations can be output
! Determine the number of, NWaveKin0Prime, and the zi-coordinates for,
! WaveKinzi0Prime(:), points along a vertical line passing through the
! platform reference point where the incident wave kinematics will be
! computed before applying stretching to the instantaneous free surface.
! The locations are relative to the mean see level. Also determine J_Min,
! which is the minimum value of index J such that WaveKinzi0(J) >=
! -WtrDpth. These depend on which incident wave kinematics stretching
! method is being used:
!JASON: ADD OTHER STRETCHING METHODS HERE, SUCH AS: DELTA STRETCHING (SEE ISO 19901-1) OR CHAKRABARTI STRETCHING (SEE OWTES)???
!JASON: APPLY STRETCHING TO THE DYNAMIC PRESSURE, IF YOU EVER COMPUTE THAT HERE!!!
SELECT CASE ( WaveStMod ) ! Which model are we using to extrapolate the incident wave kinematics to the instantaneous free surface?
CASE ( 0 ) ! None=no stretching.
! Since we have no stretching, NWaveKin0Prime and WaveKinzi0Prime(:) are
! equal to the number of, and the zi-coordinates for, the points in the
! WaveKinzi0(:) array between, and including, -WtrDpth and 0.0.
! Determine J_Min and NWaveKin0Prime here:
J_Min = 0
NWaveKin0Prime = 0
DO J = 1,NWaveKin0 ! Loop through all points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed
IF ( WaveKinzi0(J) >= -WtrDpth ) THEN
IF ( J_Min == 0 ) J_Min = J
IF ( WaveKinzi0(J) <= 0.0 ) THEN
NWaveKin0Prime = NWaveKin0Prime + 1
ELSE
EXIT
ENDIF
ENDIF
ENDDO ! J - All points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed
! ALLOCATE the WaveKinzi0Prime(:) array and compute its elements here:
ALLOCATE ( WaveKinzi0Prime(NWaveKin0Prime) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveKinzi0Prime array.')
ENDIF
DO J = 1,NWaveKin0Prime ! Loop through all points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed without stretching
WaveKinzi0Prime(J) = WaveKinzi0(J+J_Min-1)
ENDDO ! J - All points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed without stretching
CASE ( 1, 2 ) ! Vertical stretching or extrapolation stretching.
! Vertical stretching says that the wave kinematics above the mean sea level
! equal the wave kinematics at the mean sea level. The wave kinematics
! below the mean sea level are left unchanged.
!
! Extrapolation stretching uses a linear Taylor expansion of the wave
! kinematics (and their partial derivatives with respect to z) at the mean
! sea level to find the wave kinematics above the mean sea level. The
! wave kinematics below the mean sea level are left unchanged.
!
! Vertical stretching and extrapolation stretching do not effect the wave
! kinematics below the mean sea level; also, vertical stretching and
! extrapolation stretching say the wave kinematics above the mean sea
! level depend only on the mean sea level values. Consequently,
! NWaveKin0Prime and WaveKinzi0Prime(:) are equal to the number of, and
! the zi-coordinates for, the points in the WaveKinzi0(:) array between,
! and including, -WtrDpth and 0.0; the WaveKinzi0Prime(:) array must also
! include 0.0 even if the WaveKinzi0(:) array does not.
! Determine J_Min and NWaveKin0Prime here:
J_Min = 0
NWaveKin0Prime = 0
DO J = 1,NWaveKin0 ! Loop through all points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed
IF ( WaveKinzi0(J) >= -WtrDpth ) THEN
IF ( J_Min == 0 ) J_Min = J
NWaveKin0Prime = NWaveKin0Prime + 1
IF ( WaveKinzi0(J) >= 0.0 ) EXIT
ENDIF
ENDDO ! J - All points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed
! ALLOCATE the WaveKinzi0Prime(:) array and compute its elements here:
ALLOCATE ( WaveKinzi0Prime(NWaveKin0Prime) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveKinzi0Prime array.')
ENDIF
DO J = 1,NWaveKin0Prime ! Loop through all points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed without stretching
WaveKinzi0Prime(J) = MIN( WaveKinzi0(J+J_Min-1), 0.0 ) ! The uppermost point is always zero even if WaveKinzi0(NWaveKin0Prime+J_Min-1) > 0.0
ENDDO ! J - All points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed without stretching
CASE ( 3 ) ! Wheeler stretching.
! Wheeler stretching says that wave kinematics calculated using Airy theory
! at the mean sea level should actually be applied at the instantaneous
! free surface and that Airy wave kinematics computed at locations between
! the seabed and the mean sea level should be shifted vertically to new
! locations in proportion to their elevation above the seabed.
!
! Thus, given a range of zi(:) where we want to know the wave kinematics
! after applying Wheeler stretching, the required range of ziPrime(:)
! where the wave kinematics need to be computed before applying
! stretching, is as follows:
!
! ziPrime_Min <= ziPrime(:) <= ziPrime_Max
!
! ziPrime_Min = MAX{ ( zi_Min - WaveElev_Max )/( 1 + WaveElev_Max/WtrDpth ), -WtrDpth }
! ziPrime_Max = MIN{ ( zi_Max - WaveElev_Min )/( 1 + WaveElev_Min/WtrDpth ), 0 }
!
! where,
! zi_Max = maximum elevation where the wave kinematics are to be
! applied using stretching to the instantaneous free
! surface
! zi_Min = minimum elevation where the wave kinematics are to be
! applied using stretching to the instantaneous free
! surface
! ziPrime_Max = maximum elevation where the wave kinematics are computed
! before applying stretching to the instantaneous free
! surface
! ziPrime_Min = minimum elevation where the wave kinematics are computed
! before applying stretching to the instantaneous free
! surface
! WaveElev_Max = maximum expected value of the instantaneous elevation of
! incident waves
! WaveElev_Min = minimum expected value of the instantaneous elevation of
! incident waves
!
! Thus, in order to account for Wheeler stretching when computing the wave
! kinematics at each of the NWaveKin0 points along a vertical line passing
! through the platform reference point [defined by the zi-coordinates
! relative to the mean see level as specified in the WaveKinzi0(:) array],
! we must first compute the wave kinematics without stretching at
! alternative elevations [indicated here by the NWaveKin0Prime-element
! array WaveKinzi0Prime(:)]:
IF ( NWaveKin0 > 0 ) THEN ! .TRUE. if we have at least one point along a vertical line passing through the platform reference point where the incident wave kinematics will be computed
WaveElev_Max = WaveHs ! The maximum expected value the instantaneous wave elevation will most likely not exceed the value of WaveHs above the MSL since 99.99366% of all instantaneous wave elevations fall within WaveHs = 4*( the standard deviation of the incident waves ) assuming that the instanteous wave elevation is Gaussian distributed with zero mean (as implemented).
WaveElev_Min = -WaveHs ! The mimimum expected value the instantaneous wave elevation will most likely not exceed the value of WaveHs below the MSL since 99.99366% of all instantaneous wave elevations fall within WaveHs = 4*( the standard deviation of the incident waves ) assuming that the instanteous wave elevation is Gaussian distributed with zero mean (as implemented).
ziPrime_Min = MAX( WheelerStretching ( WaveKinzi0( 1), WaveElev_Max, WtrDpth, 'B' ), -WtrDpth )
ziPrime_Max = MIN( WheelerStretching ( WaveKinzi0(NWaveKin0), WaveElev_Min, WtrDpth, 'B' ), 0.0 )
IF ( MIN( ziPrime_Min, 0.0) == MAX( ziPrime_Max, -WtrDpth ) ) THEN ! .TRUE. only when all of the WaveKinzi0(:) elevations are lower than the seabed or higher than the maximum wave (inclusive); thus, we will not have any points to compute wave kinematics without stretching
NWaveKin0Prime = 0.0
ELSE ! At least one of the WaveKinzi0(:) elevations must lie within the water
! Determine NWaveKin0Prime here; no reason to compute J_Min here, so don't:
! NOTE: See explanation of stretching above for more information.
DO J = 1,NWaveKin0 ! Loop through all points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed
IF ( WheelerStretching ( WaveKinzi0(J), WaveElev_Max, WtrDpth, 'B' ) <= ziPrime_Min ) THEN
I = J
NWaveKin0Prime = 1
ELSEIF ( WheelerStretching ( WaveKinzi0(J), WaveElev_Min, WtrDpth, 'B' ) >= ziPrime_Max ) THEN
NWaveKin0Prime = NWaveKin0Prime + 1
DO K = J+1,NWaveKin0 ! Loop through all remaining points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed
IF ( WaveKinzi0(K) >= WaveElev_Max ) THEN
NWaveKin0Prime = NWaveKin0Prime + 1
EXIT ! EXIT this DO...LOOP
ELSE
NWaveKin0Prime = NWaveKin0Prime + 1
ENDIF
ENDDO ! K - All remaining points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed
EXIT ! EXIT this DO...LOOP
ELSE
NWaveKin0Prime = NWaveKin0Prime + 1
ENDIF
ENDDO ! J - All points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed
IF ( NWaveKin0Prime > 1 ) THEN ! .TRUE. if zi_Min /= zi_Max
zi_Min = MAX( WaveKinzi0( I ), -WtrDpth )
zi_Max = MIN( WaveKinzi0(NWaveKin0Prime+I-1), WaveElev_Max )
Slope = ( ziPrime_Max - ziPrime_Min )/( zi_Max - zi_Min )
ELSE ! we must have zi_Min == Zi_Max, but we still have ziPrime_Min /= ziPrime_Max
NWaveKin0Prime = 2
ENDIF
! ALLOCATE the WaveKinzi0Prime(:) array and compute its elements here:
! NOTE: See explanation of stretching above for more information.
ALLOCATE ( WaveKinzi0Prime(NWaveKin0Prime) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveKinzi0Prime array.')
ENDIF
WaveKinzi0Prime( 1) = ziPrime_Min ! First point - lowermost
DO J = 2,NWaveKin0Prime-1 ! Loop through all but the first and last points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed without stretching
WaveKinzi0Prime( J) = ( WaveKinzi0(J+I-1) - zi_Min )*Slope + ziPrime_Min ! Interpolate to find the middle points using the elevations of the WaveKinzi0(:) array
ENDDO ! J - All but the first and last points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed without stretching
WaveKinzi0Prime(NWaveKin0Prime) = ziPrime_Max ! Last point - uppermost
ENDIF
ELSE ! We must not have any point along a vertical line passing through the platform reference point where the incident wave kinematics will be computed; thus, neither will we have any points to compute wave kinematics without stretching
NWaveKin0Prime = 0
ENDIF
ENDSELECT
! Perform some initialization computations including initializing the
! pseudorandom number generator, calculating the total number of frequency
! components = total number of time steps in the incident wave,
! calculating the frequency step, calculating the index of the frequency
! component nearest to WaveTp, and ALLOCATing the arrays:
! NOTE: WaveDOmega = 2*Pi/WaveTMax since, in the FFT:
! Omega = (K-1)*WaveDOmega
! Time = (J-1)*WaveDT
! and therefore:
! Omega*Time = (K-1)*(J-1)*WaveDOmega*WaveDT
! = (K-1)*(J-1)*2*Pi/NStepWave [see FFT_Module]
! or:
! WaveDOmega = 2*Pi/(NStepWave*WaveDT)
! = 2*Pi/WaveTMax
CALL RANDOM_SEED ( SIZE = nSeeds )
IF ( nSeeds /= 2 ) THEN
CALL ProgWarn( ' The random number generator in use differs from the original code provided by NREL. This pRNG uses ' &
//TRIM(Int2LStr(nSeeds))//' seeds instead of the 2 in the HydroDyn input file.')
END IF
ALLOCATE ( TmpWaveSeeds ( nSeeds ), STAT=Sttus)
IF (Sttus/= 0 ) THEN
CALL ProgAbort( ' Error allocating space for TmpWaveSeeds array.' )
RETURN
END IF
! We'll just populate this with odd seeds = Seed(1) and even seeds = Seed(2)
DO I = 1,nSeeds,2
TmpWaveSeeds(I) = WaveSeed(1)
END DO
DO I = 2,nSeeds,2
TmpWaveSeeds(I) = WaveSeed(2)
END DO
CALL RANDOM_SEED ( PUT=TmpWaveSeeds )
DEALLOCATE(TmpWaveSeeds, STAT=Sttus)
IF (Sttus/= 0 ) THEN
CALL ProgWarn( ' Error deallocating space for TmpWaveSeeds array.' )
END IF
NStepWave = CEILING ( WaveTMax/WaveDT ) ! Set NStepWave to an even integer
IF ( MOD(NStepWave,2) == 1 ) NStepWave = NStepWave + 1 ! larger or equal to WaveTMax/WaveDT.
NStepWave2 = MAX( NStepWave/2, 1 ) ! Make sure that NStepWave is an even product of small factors (PSF) that is
NStepWave = 2*PSF ( NStepWave2, 9 ) ! greater or equal to WaveTMax/WaveDT to ensure that the FFT is efficient.
NStepWave2 = NStepWave/2 ! Update the value of NStepWave2 based on the value needed for NStepWave.
WaveTMax = NStepWave*WaveDT ! Update the value of WaveTMax based on the value needed for NStepWave.
WaveDOmega = TwoPi/WaveTMax ! Compute the frequency step for incident wave calculations.
I_WaveTp = NINT ( TwoPi/(WaveDOmega*WaveTp) ) ! Compute the index of the frequency component nearest to WaveTp.
IF ( WaveMod == 2 ) OmegaCutOff = n_Massel*TwoPi/WaveTp ! Compute the cut-off frequency or upper frequency limit of the wave spectrum beyond which the wave spectrum is zeroed. The TwoPi/WaveTp is the peak spectral frequency in rad/s; the cut-off frequency is a factor of N_Massel above this value based on the suggestion in: Massel, S. R., Ocean Surface Waves: Their Physics and Prediction, Advanced Series on Ocean Engineering - Vol. 11, World Scientific Publishing, Singapore - New Jersey - London - Hong Kong, 1996.
ALLOCATE ( WaveTime (0:NStepWave-1 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveTime array.')
ENDIF
ALLOCATE ( WaveElevC0 (0:NStepWave2 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveElevC0 array.')
ENDIF
ALLOCATE ( WaveElevC (0:NStepWave2 ,NWaveElev ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveElevC array.')
ENDIF
ALLOCATE ( WaveVelC0H (0:NStepWave2 ,NWaveKin0Prime ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveVelC0H array.')
ENDIF
ALLOCATE ( WaveVelC0V (0:NStepWave2 ,NWaveKin0Prime ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveVelC0V array.')
ENDIF
ALLOCATE ( WaveAccC0H (0:NStepWave2 ,NWaveKin0Prime ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveAccC0H array.')
ENDIF
ALLOCATE ( WaveAccC0V (0:NStepWave2 ,NWaveKin0Prime ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveAccC0V array.')
ENDIF
ALLOCATE ( PWaveVelC0HPz0(0:NStepWave2 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the PWaveVelC0HPz0 array.')
ENDIF
ALLOCATE ( PWaveVelC0VPz0(0:NStepWave2 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the PWaveVelC0VPz0 array.')
ENDIF
ALLOCATE ( PWaveAccC0HPz0(0:NStepWave2 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the PWaveAccC0HPz0 array.')
ENDIF
ALLOCATE ( PWaveAccC0VPz0(0:NStepWave2 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the PWaveAccC0VPz0 array.')
ENDIF
ALLOCATE ( WaveElev0 (0:NStepWave-1 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveElev0 array.')
ENDIF
ALLOCATE ( WaveElev (0:NStepWave-1,NWaveElev ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveElev array.')
ENDIF
ALLOCATE ( WaveVel0H (0:NStepWave-1,NWaveKin0Prime ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveVel0H array.')
ENDIF
ALLOCATE ( WaveVel0Hxi (0:NStepWave-1,NWaveKin0Prime ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveVel0Hxi array.')
ENDIF
ALLOCATE ( WaveVel0Hyi (0:NStepWave-1,NWaveKin0Prime ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveVel0Hyi array.')
ENDIF
ALLOCATE ( WaveVel0V (0:NStepWave-1,NWaveKin0Prime ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveVel0V array.')
ENDIF
ALLOCATE ( WaveAcc0H (0:NStepWave-1,NWaveKin0Prime ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveAcc0H array.')
ENDIF
ALLOCATE ( WaveAcc0V (0:NStepWave-1,NWaveKin0Prime ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveAcc0V array.')
ENDIF
ALLOCATE ( PWaveVel0HPz0 (0:NStepWave-1 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the PWaveVel0HPz0 array.')
ENDIF
ALLOCATE ( PWaveVel0HxiPz0 (0:NStepWave-1 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the PWaveVel0HxiPz0 array.')
ENDIF
ALLOCATE ( PWaveVel0HyiPz0 (0:NStepWave-1 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the PWaveVel0HyiPz0 array.')
ENDIF
ALLOCATE ( PWaveVel0VPz0 (0:NStepWave-1 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the PWaveVel0VPz0 array.')
ENDIF
ALLOCATE ( PWaveAcc0HPz0 (0:NStepWave-1 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the PWaveAcc0HPz0 array.')
ENDIF
ALLOCATE ( PWaveAcc0VPz0 (0:NStepWave-1 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the PWaveAcc0VPz0 array.')
ENDIF
ALLOCATE ( WaveVel0 (0:NStepWave-1,NWaveKin0,3 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveVel0 array.')
ENDIF
ALLOCATE ( WaveAcc0 (0:NStepWave-1,NWaveKin0,3 ) , STAT=Sttus )
IF ( Sttus /= 0 ) THEN
CALL ProgAbort(' Error allocating memory for the WaveAcc0 array.')
ENDIF
! Calculate the factors needed by the discrete time inverse Fourier
! transform in the calculations of the White Gaussian Noise (WGN) and
! the two-sided power spectral density of the wave spectrum per unit time:
WGNC_Fact = SQRT( Pi/(WaveDOmega*WaveDT) ) ! This factor is needed by the discrete time inverse Fourier transform to ensure that the time series WGN process has unit variance
S2Sd_Fact = 1.0/WaveDT ! This factor is also needed by the discrete time inverse Fourier transform
! Compute the positive-frequency components (including zero) of the Fourier
! transforms of the wave kinematics:
DO I = 0,NStepWave2 ! Loop through the positive frequency components (including zero) of the Fourier transforms
! Compute the frequency of this component and its imaginary value:
Omega = I* WaveDOmega
ImagOmega = ImagNmbr*Omega
! Compute the Fourier transform of the realization of a White Gaussian Noise
! (WGN) time series process with unit variance:
! NOTE: For the time series process to be real with zero mean, the values at
! Omega == 0.0 and Omega == NStepWave2*WaveDOmega (= WaveOmegaMax)
! must be zero.
IF ( ( I == 0 ) .OR. ( I == NStepWave2 ) ) THEN ! .TRUE. if ( Omega == 0.0 ) or ( Omega == NStepWave2*WaveDOmega (= WaveOmegaMax) )
WGNC = (0.0,0.0)
ELSE ! All other Omega
WGNC = BoxMuller ( )
IF ( ( WaveMod == 1 ) .AND. ( I == I_WaveTp ) ) WGNC = WGNC*( SQRT(2.0)/ABS(WGNC) ) ! This scaling of WGNC is used to ensure that the Box-Muller method is only providing a random phase, not a magnitude change, at the frequency of the plane progressive wave. The SQRT(2.0) is used to ensure that the time series WGN process has unit variance (i.e. sinusoidal with amplitude SQRT(2.0)). NOTE: the denominator here will never equal zero since U1 cannot equal 1.0, and thus, C1 cannot be 0.0 in the Box-Muller method.
ENDIF
! Compute the one-sided power spectral density of the wave spectrum per unit
! time; zero-out the wave spectrum above the cut-off frequency:
SELECT CASE ( WaveMod ) ! Which incident wave kinematics model are we using?
CASE ( 1 ) ! Plane progressive (regular) wave; the wave spectrum is an impulse function centered on frequency component closest to WaveTp.
IF ( I == I_WaveTp ) THEN ! .TRUE. if we are at the Omega closest to WaveTp.
WaveS1Sdd = 0.5*(WaveHs/2.0)*(WaveHs/2.0)/WaveDOmega
ELSE ! All other Omega
WaveS1Sdd = 0.0
ENDIF
CASE ( 2 ) ! JONSWAP/Pierson-Moskowitz spectrum (irregular) wave.
IF ( Omega > OmegaCutOff ) THEN ! .TRUE. if Omega is above the cut-off frequency
WaveS1Sdd = 0.0 ! Zero-out the wave spectrum above the cut-off frequency. We must cut-off the frequency in order to avoid nonphysical wave forces. Waves that have wavelengths much smaller than the platform diameter (high frequency) do not contribute to the net force because regions of positive and negative velocity/acceleration are experienced by the platform at the same time and cancel out. !JASON: OTHER FREQUENCY CUT-OFF CONDITIONS ARE USED THROUGHOUT THE INDUSTRY. SHOULD YOU USE ONE OF THEM INSTEAD? SEE, FOR EXAMPLE, MY E-MAIL EXCHANGES WITH PAUL SCLAVOUNOS DATED 5/26/2006 OR: "GH Bladed Thoery Manual" OR: Trumars, Jenny M. V.; Tarp-Johansen, Niels Jacob; Krogh, Thomas; "The Effect of Wave Modelling on Offshore Wind Turbine Fatigue Loads," 2005 Copenhagen Offshore Wind Conference and Exhibition, 26-28 October 2005, Copenhagen, Denmark [CD-ROM].
ELSE ! All other Omega
WaveS1Sdd = JONSWAP ( Omega, WaveHs, WaveTp, WavePkShp )
ENDIF
CASE ( 3 ) ! User-defined spectrum (irregular) wave.
CALL UserWaveSpctrm ( Omega, WaveDir, DirRoot, WaveS1Sdd )
ENDSELECT
! Compute the two-sided power spectral density of the wave spectrum per unit
! time:
WaveS2Sdd = 0.5*WaveS1Sdd
! Compute the wavenumber:
WaveNmbr = WaveNumber ( Omega, Gravity, WtrDpth )
! Compute the Fourier transform of the instantaneous elevation of incident
! waves at the platform reference point:
WaveElevC0 (I ) = ( WGNC_Fact*WGNC )*SQRT( TwoPi*( S2Sd_Fact*WaveS2Sdd ) )
! Compute the Fourier transform of the instantaneous elevation of incident
! waves at each desired point on the still water level plane where it can
! be output:
DO J = 1,NWaveElev ! Loop through all points where the incident wave elevations can be output
WaveElevC (I,J) = WaveElevC0 (I )*EXP( -ImagNmbr*WaveNmbr*WaveElevxiPrime(J) )
ENDDO ! J - All points where the incident wave elevations can be output
! Compute the Fourier transform of the incident wave kinematics before
! applying stretching at the zi-coordinates for points along a vertical
! line passing through the platform reference point:
DO J = 1,NWaveKin0Prime ! Loop through all points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed without stretching
WaveVelC0H (I,J) = Omega* WaveElevC0 (I )*COSHNumOvrSIHNDen ( WaveNmbr, WtrDpth, WaveKinzi0Prime(J) )
WaveVelC0V (I,J) = ImagOmega* WaveElevC0 (I )*SINHNumOvrSIHNDen ( WaveNmbr, WtrDpth, WaveKinzi0Prime(J) )
WaveAccC0H (I,J) = ImagOmega* WaveVelC0H (I,J)
WaveAccC0V (I,J) = ImagOmega* WaveVelC0V (I,J)
ENDDO ! J - All points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed without stretching
PWaveVelC0HPz0(I ) = Omega* WaveElevC0 (I )*WaveNmbr
PWaveVelC0VPz0(I ) = ImagOmega* WaveElevC0 (I )*WaveNmbr*COTH ( WaveNmbr*WtrDpth )
PWaveAccC0HPz0(I ) = ImagOmega*PWaveVelC0HPz0(I )
PWaveAccC0VPz0(I ) = ImagOmega*PWaveVelC0VPz0(I )
ENDDO ! I - The positive frequency components (including zero) of the Fourier transforms
! Calculate the array of simulation times at which the instantaneous
! elevation of, velocity of, acceleration of, and loads associated with
! the incident waves are to be determined:
DO I = 0,NStepWave-1 ! Loop through all time steps
WaveTime(I) = I*WaveDT
ENDDO ! I - All time steps
! Compute the inverse Fourier transforms to find the time-domain
! representations of the wave kinematics without stretcing:
CALL InitFFT ( NStepWave, .TRUE. )
CALL ApplyFFT_cx ( WaveElev0 (: ), WaveElevC0 (: ) )
DO J = 1,NWaveElev ! Loop through all points where the incident wave elevations can be output
CALL ApplyFFT_cx ( WaveElev (:,J), WaveElevC (:,J) )
ENDDO ! J - All points where the incident wave elevations can be output
DO J = 1,NWaveKin0Prime ! Loop through all points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed without stretching
CALL ApplyFFT_cx ( WaveVel0H (:,J), WaveVelC0H (:,J) )
CALL ApplyFFT_cx ( WaveVel0V (:,J), WaveVelC0V (:,J) )
CALL ApplyFFT_cx ( WaveAcc0H (:,J), WaveAccC0H (:,J) )
CALL ApplyFFT_cx ( WaveAcc0V (:,J), WaveAccC0V (:,J) )
ENDDO ! J - All points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed without stretching
CALL ApplyFFT_cx ( PWaveVel0HPz0(: ), PWaveVelC0HPz0(: ) )
CALL ApplyFFT_cx ( PWaveVel0VPz0(: ), PWaveVelC0VPz0(: ) )
CALL ApplyFFT_cx ( PWaveAcc0HPz0(: ), PWaveAccC0HPz0(: ) )
CALL ApplyFFT_cx ( PWaveAcc0VPz0(: ), PWaveAccC0VPz0(: ) )
CALL ExitFFT
! Add the current velocities to the wave velocities:
! NOTE: Both the horizontal velocities and the partial derivative of the
! horizontal velocities with respect to zi at zi = 0 are found here.
DO J = 1,NWaveKin0Prime ! Loop through all points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed without stretching
CALL InitCurrent ( CurrMod , CurrSSV0 , CurrSSDir, CurrNSRef, &
CurrNSV0 , CurrNSDir, CurrDIV , CurrDIDir, &
WaveKinzi0Prime(J), WtrDpth , DirRoot , CurrVxi , CurrVyi )
WaveVel0Hxi (:,J) = WaveVel0H (:,J)*CWaveDir + CurrVxi ! xi-direction
WaveVel0Hyi (:,J) = WaveVel0H (:,J)*SWaveDir + CurrVyi ! yi-direction
ENDDO ! J - All points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed without stretching
CALL InitCurrent ( CurrMod , CurrSSV0 , CurrSSDir, CurrNSRef, &
CurrNSV0 , CurrNSDir, CurrDIV , CurrDIDir, &
0.0 , WtrDpth , DirRoot , CurrVxi0 , CurrVyi0 )
CALL InitCurrent ( CurrMod , CurrSSV0 , CurrSSDir, CurrNSRef, &
CurrNSV0 , CurrNSDir, CurrDIV , CurrDIDir, &
-SmllNmbr , WtrDpth , DirRoot , CurrVxiS , CurrVyiS )
PCurrVxiPz0 = ( CurrVxi0 - CurrVxiS )/SmllNmbr ! xi-direction
PCurrVyiPz0 = ( CurrVyi0 - CurrVyiS )/SmllNmbr ! yi-direction
PWaveVel0HxiPz0(: ) = PWaveVel0HPz0(: )*CWaveDir + PCurrVxiPz0 ! xi-direction
PWaveVel0HyiPz0(: ) = PWaveVel0HPz0(: )*SWaveDir + PCurrVyiPz0 ! yi-direction
! Apply stretching to obtain the wave kinematics, WaveVel0 and WaveAcc0, at
! the desired locations from the wave kinematics at alternative locations,
! WaveVel0Hxi, WaveVel0Hyi, WaveVel0V, WaveAcc0H, WaveAcc0V, if the elevation
! of the point defined by WaveKinzi0(J) lies between the seabed and the
! instantaneous free surface, else set WaveVel0 and WaveAcc0 to zero.
! This depends on which incident wave kinematics stretching method is
! being used:
SELECT CASE ( WaveStMod ) ! Which model are we using to extrapolate the incident wave kinematics to the instantaneous free surface?
CASE ( 0 ) ! None=no stretching.
! Since we have no stretching, the wave kinematics between the seabed and
! the mean sea level are left unchanged; below the seabed or above the
! mean sea level, the wave kinematics are zero:
DO I = 0,NStepWave-1 ! Loop through all time steps
DO J = 1,NWaveKin0 ! Loop through all points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed
IF ( ( WaveKinzi0(J) < -WtrDpth ) .OR. ( WaveKinzi0(J) > 0.0 ) ) THEN ! .TRUE. if the elevation of the point defined by WaveKinzi0(J) lies below the seabed or above mean sea level (exclusive)
WaveVel0(I,J,:) = 0.0
WaveAcc0(I,J,:) = 0.0
ELSE ! The elevation of the point defined by WaveKinzi0(J) must lie between the seabed and the mean sea level (inclusive)
WaveVel0(I,J,1) = WaveVel0Hxi(I,J-J_Min+1 )
WaveVel0(I,J,2) = WaveVel0Hyi(I,J-J_Min+1 )
WaveVel0(I,J,3) = WaveVel0V (I,J-J_Min+1 )
WaveAcc0(I,J,1) = WaveAcc0H (I,J-J_Min+1 )*CWaveDir
WaveAcc0(I,J,2) = WaveAcc0H (I,J-J_Min+1 )*SWaveDir
WaveAcc0(I,J,3) = WaveAcc0V (I,J-J_Min+1 )
ENDIF
ENDDO ! J - All points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed
ENDDO ! I - All time steps
CASE ( 1 ) ! Vertical stretching.
! Vertical stretching says that the wave kinematics above the mean sea level
! equal the wave kinematics at the mean sea level. The wave kinematics
! below the mean sea level are left unchanged:
DO I = 0,NStepWave-1 ! Loop through all time steps
DO J = 1,NWaveKin0 ! Loop through all points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed
IF ( ( WaveKinzi0(J) < -WtrDpth ) .OR. ( WaveKinzi0(J) > WaveElev0(I) ) ) THEN ! .TRUE. if the elevation of the point defined by WaveKinzi0(J) lies below the seabed or above the instantaneous free surface (exclusive)
WaveVel0(I,J,:) = 0.0
WaveAcc0(I,J,:) = 0.0
ELSEIF ( WaveKinzi0(J) > 0.0 ) THEN ! .TRUE. if the elevation of the point devined by WaveKinzi0(J) lies between the mean sea level (exclusive) and the instantaneous free surface (inclusive)
WaveVel0(I,J,1) = WaveVel0Hxi(I,NWaveKin0Prime)
WaveVel0(I,J,2) = WaveVel0Hyi(I,NWaveKin0Prime)
WaveVel0(I,J,3) = WaveVel0V (I,NWaveKin0Prime)
WaveAcc0(I,J,1) = WaveAcc0H (I,NWaveKin0Prime)*CWaveDir
WaveAcc0(I,J,2) = WaveAcc0H (I,NWaveKin0Prime)*SWaveDir
WaveAcc0(I,J,3) = WaveAcc0V (I,NWaveKin0Prime)
ELSE ! The elevation of the point defined by WaveKinzi0(J) must lie between the seabed and the mean sea level (inclusive)
WaveVel0(I,J,1) = WaveVel0Hxi(I,J-J_Min+1 )
WaveVel0(I,J,2) = WaveVel0Hyi(I,J-J_Min+1 )
WaveVel0(I,J,3) = WaveVel0V (I,J-J_Min+1 )
WaveAcc0(I,J,1) = WaveAcc0H (I,J-J_Min+1 )*CWaveDir
WaveAcc0(I,J,2) = WaveAcc0H (I,J-J_Min+1 )*SWaveDir
WaveAcc0(I,J,3) = WaveAcc0V (I,J-J_Min+1 )
ENDIF
ENDDO ! J - All points along a vertical line passing through the platform reference point where the incident wave kinematics will be computed
ENDDO ! I - All time steps