-
Notifications
You must be signed in to change notification settings - Fork 0
/
SdsBabanin.ftn90
2708 lines (2462 loc) · 101 KB
/
SdsBabanin.ftn90
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
! This file contains subroutines for the Babanin physics according to Rogers et al (JTECH 2012)
! based on work of Babanin, Young, Tsagareli, Ardhuin and others
MODULE SDSBABANIN
CONTAINS
SUBROUTINE CALC_SDS(NFREQ,EDENS,F,KDS,ANAR_IN,TESTFL,KWAVE,CG)
!
USE SWCOMM1, ONLY: CHTIME
USE SWCOMM3, ONLY: A1SDS,A2SDS,P1SDS,P2SDS,UPWARDS,GRAV,PI
! REAL A1SDS : coefficient on T1
! REAL A2SDS : coefficient on T2
! REAL P1SDS : power on T1
! REAL P2SDS : power on T2
! P1SDS,P2SDS is the power on (F-Fth)/Fth or (F-Fth)/F which is part of
! how I formulate "X" given by Young and Babanin eq 27 :
! P1SDS: power on threshold exceedence in T1 (precise definition depends
! on whether using "U" vs "D" method)
! P2SDS: power on threshold exceedence in T2 (precise definition depends on
! whether using "U" vs "D" method)
! LOGICAL UPWARDS : true if concave up
IMPLICIT NONE
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2024 Delft University of Technology
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program. If not, see <http://www.gnu.org/licenses/>.
!
!
! INPUT: EDENS(f), f, OUTPUT: calculate Kds(f) [Kds(f)=Sds(f)/EDENS(f)]
! EDENS is given wave spectrum m2/Hz
! f is frequency in Hz example: f=0.01:0.01:0.5
! @EDENS/@t=Sin+Snl+Sds
!
! Kds=T1+T2
!
! T1=a1*f*A*DEDENS
! (inherent breaking dissipation)
! T2=a2*cumsum(A*DEDENS) an integration of A*DEDENS from fp to f
! (induced breaking dissipation)
!
! A is directional narrowness, calculated in calling subroutine, but set
! to unity in this version
!
! DEDENS=EDENS-EDENST
! EDENST=(2*g.^2)*((2*pi).^-4)*(f.^-5).*(A.^-1)*Bnt
! a1 and a2 are function of U/cp given by Babanin
! Bnt is an empirical coefficient related to the spectral density in the
! saturation range (Babanin et al 2007)
!
! Purpose: to calculate the spectral dissipation Sds based on given wave
! spectrum
! References:
! 1) Babanin and Young ("WAVES" conf. 2005)
! 2) Young and Babanin (JPO 2006)
! 3) Babanin et al. (H/F conf. 2007), "Implementation of new ..."
! 4) WISE presentations by Rogers in 2009, 2010
! 5) Babanin et al. (JPO 2010,) paper re: Sds term
! 6) Tsagareli et al. (JPO 2010), paper re: Sin term
! 7) Rogers et al. (JTECH 2012)
! The code is "unified" insofar as depending on the setting of the
! logical "UPWARDS", one gets something similar to the method of :
! =true: Fabrice et al. (more similar to the early, nondirectional
! version of Fabrice's Sds published in a conference paper, as
! opposed to the Ardhuin et al. JPO 2010 directional Sds). With L=2
! and large exceedence, the dissipation will get very strong
! NDEDENS(is)=DEDENS(is)/EDENST(is)
! =false: Tsagareli, Babanin, Young. (but it differs from theirs
! insofar as our normalization for the cumulative term is more correct
! mathemetically, so we do not have to created a set of a1, b1 for each
! wave age...with the normalization used here, one calibration works
! for all wave age values. Another improvement is in the exponent L :
! in the Tsagereli formula, the only dimensionally correct form of the
! equation is for L=1. This means that the normalized value is always
! less than 1, which means that using L=2 instead of L=1 will
! actually make the dissipation weaker (though since a1 a2 would need
! to be recalibrated, this relation is a bit more complicated)
!
! Notes on calibration, last updated Sep 19 2013.
! Because these physics are "observation-based", and in Rogers et al (2012)
! we apply a number of rules and physical constraints based on
! observations, this means that the new physics are *less* tunable than
! prior source term packages, like KHH1984. This is a bit of a paradox,
! since we probably have more free parameters than KHH. Most important
! example, we can easily control SWH (total energy) but have little
! control over Tp, Tm02, Tm01, Tm-1,0
! Possibilities for calibrations are noted in the SWAN user's manual.
!
! Further documentation with more up-to-date information than what you get
! from Rogers et al. (2012) :
! If you are inside the .navy.mil domain, see these wiki pages:
! https://www7320.nrlssc.navy.mil/Alvin/index.php/SWAN
! (especially the pages about * Babanin physics * Improvements * and
! * Suggested Settings * )
! If you are not inside the .navy.mil domain, you can ask us for
! this info: erick.rogers@nrlssc.navy.mil
!
! Subroutine arguments:
LOGICAL , INTENT(IN) :: TESTFL
! TESTFL: true if output requested for this point (use sparingly, can result
! in large files)
REAL , INTENT(IN) :: EDENS(:) ! E(f) m2/Hz
REAL , INTENT(IN) :: F(:) ! frequency in Hz
REAL , INTENT(IN) :: ANAR_IN(:)
! ANAR_IN: directional narrowness as defined in Babanin publications
! (input value)
INTEGER , INTENT(IN) :: NFREQ ! # freqs
REAL , INTENT(IN) :: KWAVE(:,:) ! KWAVE(MSC,MICMAX)
REAL , INTENT(IN) :: CG(:,:) ! CG(MSC,MICMAX)
REAL , INTENT(OUT) :: KDS(:) ! Kds(f)=Sds(f)/E(f)
! Local variables:
REAL , ALLOCATABLE :: SDS(:) ! Sds(f), the source term
REAL , ALLOCATABLE :: NDEDENS(:)
! NDEDENS(f)=DEDENS(f)/EDENST(f)
REAL , ALLOCATABLE :: DEDENS(:) ! DEDENS(f)=EDENS(f)-EDENST(f)
REAL , ALLOCATABLE :: EDENST(:) ! E(f) threshold for breaking
REAL , ALLOCATABLE :: T1(:) ! inherent dissipation/E(f)
REAL , ALLOCATABLE :: T2(:) ! induced dissipation/E(f)
REAL , ALLOCATABLE :: ST1(:) ! inherent dissipation
REAL , ALLOCATABLE :: ST2(:) ! induced dissipation
REAL , ALLOCATABLE :: ANAR(:)
! ANAR = directional narrowness as defined in Babanin publications
REAL , ALLOCATABLE :: XFF(:),ADF(:) ! temporary arrays
REAL :: ASUM ! temporary variable for integration
REAL :: BNT
! BNT is an empirical coefficient related to the spectral density in the
! saturation range (Babanin et al 2007)
INTEGER :: II,IS ! counters
! *_INT: integrated values (for test output only)
REAL :: ST1_INT,ST2_INT,SDS_INT
! IMAX,I3FP, fp, FD(:) re: T1 and T2 at 3fp (for test output only)
INTEGER :: IMAX,I3FP
REAL :: FP
REAL, ALLOCATABLE :: FD(:)
REAL :: ELIM ! needed for UPWARDS=.FALSE.
REAL :: CTMP1 ! temporary variable
! ------- START SUBROUTINE -------------------------------------------------
ALLOCATE(SDS(NFREQ))
ALLOCATE(XFF(NFREQ))
ALLOCATE(NDEDENS(NFREQ))
ALLOCATE(DEDENS(NFREQ))
ALLOCATE(EDENST(NFREQ))
ALLOCATE(T1(NFREQ))
ALLOCATE(T2(NFREQ))
ALLOCATE(ANAR(NFREQ))
ALLOCATE(ST1(NFREQ))
ALLOCATE(ST2(NFREQ))
ALLOCATE(ADF(NFREQ))
ALLOCATE(FD(NFREQ))
BNT=(0.035**2) ! Bnt value given by Babanin et al (2007)
! --------------------------------------------------------------------------
! get a1 and a2 as a function of U/cp
! --------------------------------------------------------------------------
! needed: A(f), see Young and Babanin eq 19.
! Originally, it was read in and applied, but per recommendation by Alex
! that A(f) is unnecessary/obsolete, I set it to 1.0 here:
ANAR=1.0 ! option 1
! ANAR=ANAR_IN ! option 2
! (point output write location 1)
!IF(TESTFL)WRITE(*,*)'calc_Sds : a1,a2 = ',A1SDS,A2SDS
! new calculation, depends on k and Cg instead of f :
CTMP1=2.0*PI*BNT
DO IS=1,NFREQ
EDENST(IS)=CTMP1/(CG(IS,1)*KWAVE(IS,1)**3)
END DO
! potential optimization: if ANAR is always unity, then EDENST does not vary
! with time step and only needs to be calculated once, rather than every
! time this routine is called.
DEDENS=EDENS-EDENST ! matrix operation
! if DEDENS < 0 set to zero
DEDENS=MAX(0.0,DEDENS)
! notes: Mar 22 2011, Stefan has noticed sensitivity to ELIM. (this variable
! only applies to concave down case). I have experimented with ELIM=0.0 and
! noticed no sensitivity for the U10=12 m/s point model case.
! ELIM is needed for "concave down"
! ELIM=maxval(Edens)*1.0e-5 ! option 1
ELIM=0.0 ! option 2
! See notes above re: UPWARDS
IF(UPWARDS)THEN
NDEDENS(1:NFREQ)=DEDENS(1:NFREQ)/EDENST(1:NFREQ)
ELSE
DO IS=1,NFREQ
IF(EDENS(IS).GT.ELIM)THEN
NDEDENS(IS)=DEDENS(IS)/EDENS(IS)
ELSE
NDEDENS(IS)=0.0
END IF
ENDDO
ENDIF
NDEDENS=MAX(0.0,NDEDENS)
! --------------------------------------------------------------------------
! calculate T1
! -------------------------------------------------------------------------
DO IS=1,NFREQ
T1(IS)=A1SDS*F(IS)*ANAR(IS)*NDEDENS(IS)**P1SDS
END DO
! --------------------------------------------------------------------------
! calculate T2 an integration from fp to f
! --------------------------------------------------------------------------
! Note that we do not worry about starting at fp (which can be difficult to
! define), since stuff below fp is typically not breaking, so it is not
! part of the calculation anyway.
IF(A2SDS.GT.0.)THEN
XFF=0.0
DO IS=1,NFREQ
ASUM=0.0
DO II=1,IS
XFF(II)=F(II)
ADF(II)=ANAR(II)*NDEDENS(II)**P2SDS
ENDDO
! THE "INTEGRATE" ROUTINE IS BASED ON FINITE DIFFERENCING, BUT WE COULD
! REPLACE WITH AN OPERATION THAT USES FRINTF
CALL INTEGRATE(ASUM,XFF,ADF,IS)
T2(IS)=A2SDS*ASUM
ENDDO
ELSE
T2=0.
ENDIF
T1=MAX(0.0,T1)
IF(A2SDS.GT.0.)THEN
T2=MAX(0.0,T2)
KDS=(T1+T2)
ELSE
KDS=T1
ENDIF
! (test output for NRL purposes)
!NRL IF(TESTFL)THEN
!NRL WRITE(411,*)CHTIME,' % CHTIME'
!NRL WRITE(412,*)'% f(is),EDENS(is),EDENST(is),T1(is),T2(is)'
!NRL DO IS=1,NFREQ
!NRL WRITE(412,205)F(IS),EDENS(IS),EDENST(IS),T1(IS),T2(IS)
!NRL END DO
!NRL205 FORMAT(5(1X,E11.5))
!NRL
!NRL ! calculate integrated T1 and T2 (for output purposes only)
!NRL
!NRL ! Since we have already non-dimensionalized by using NDEDENS instead of
!NRL ! DEDENS, we do it backwards: Sds=Kds*Edens.
!NRL ! Note that in most cases, the calling routine will not use Sds.
!NRL DO IS=1,NFREQ
!NRL SDS(IS)=KDS(IS)*EDENS(IS)
!NRL ST1(IS)=T1(IS)*EDENS(IS)
!NRL ST2(IS)=T2(IS)*EDENS(IS)
!NRL END DO
!NRL
!NRL CALL INTEGRATE(ST1_INT,F,ST1,NFREQ)
!NRL CALL INTEGRATE(ST2_INT,F,ST2,NFREQ)
!NRL CALL INTEGRATE(SDS_INT,F,SDS,NFREQ)
!NRL
!NRL ! calculate T1 and T2 at 3fp
!NRL IMAX=MAXLOC(EDENS,1)
!NRL FP=F(IMAX)
!NRL FD=ABS(F-FP*3.0)
!NRL I3FP=MINLOC(FD,1)
!NRL
!NRL ! (POINT OUTPUT WRITE LOCATION 3)
!NRL WRITE(*,208)ST1_INT,ST2_INT,SDS_INT,ST1(I3FP),ST2(I3FP)
!NRL208 FORMAT('integral of T1,T2,Sds = ',5(1X,E14.8))
!NRL ENDIF
! Deallocate and return
DEALLOCATE(SDS)
DEALLOCATE(XFF)
DEALLOCATE(NDEDENS)
DEALLOCATE(DEDENS)
DEALLOCATE(EDENST)
DEALLOCATE(T1)
DEALLOCATE(T2)
DEALLOCATE(ANAR)
DEALLOCATE(ST1)
DEALLOCATE(ST2)
DEALLOCATE(ADF)
DEALLOCATE(FD)
END SUBROUTINE CALC_SDS
!****************************************************************************
SUBROUTINE SWIND_DBYB ( SPCSIG,THETAW,KWAVE,MEMSINA,MEMSINB, &
AC2,UFRIC,WIND10,SPCDIR,ANYWND,CG &
,ZELEN )
!****************************************************************************
USE SWCOMM1, ONLY: CHTIME
USE SWCOMM3
USE SWCOMM4 ! includes TESTFL
USE OCPCOMM4
!ESMF USE M_GENARR, ONLY: SAVE_SINBAC, SINBAC
IMPLICIT NONE
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2024 Delft University of Technology
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program. If not, see <http://www.gnu.org/licenses/>.
!
!
! SUBROUTINE SWIND_DBYB follows the structure of SWIND3 of SWAN.
! As such, it contains features that are not general to all wave models.
! In particular, the quadrant sweeping creates special challenges that
! are not general.
!
! 0. Authors
!
! W.E. Rogers, NRL 7320
!
! 1. Updates
!
! 2. PURPOSE
!
! Computation of the source term for the wind input for a
! third generation wind growth model:
!
! 1) Exponential input term, (Donelan et al. JPO 2006)
!
! 3. METHOD
!
! Sin (s,d) = B*E(s,d)
!
! 4. Argument variables
!
! i SPCDIR: (*,1); spectral directions (radians)
! (*,2); cosine of spectral directions
! (*,3); sine of spectral directions
! (*,4); cosine^2 of spectral directions
! (*,5); cosine*sine of spectral directions
! (*,6); sine^2 of spectral directions
! i SPCSIG: Relative frequencies in computational domain in sigma-space
!
! IS Counter of relative frequency band
! ID Counter of directional distribution
! MSC Maximum counter of relative frequency
! MDC Maximum counter of directional distribution
!
! REALS:
! ---------
! THETA Spectral direction
! THETAW Mean direction of the relative wind vector
! UFRIC Wind friction velocity
!
! one and more dimensional arrays:
! ---------------------------------
! KWAVE 2D Wavenumber
! PWIND 1D Wind coefficients
! ANYWND 1D Wind input for bin considered
!
! 5. SUBROUTINES CALLING
!
! SOURCE
!
! 6. SUBROUTINES USED
!
! calc_Lfactor
!
! 7. ERROR MESSAGES
!
! ---
!
! 8. REMARKS
!
! Here is what I had prior to Oct 4 2011:
! U10PROXY=28.0*UFRIC
!
! U10=28U* implies a fixed Cd. In order for the model to scale with U*,
! we must use a fixed Cd to convert the Donelan formula from U10 to U*.
! If we use a proper, realistic Cd to do this, then we are basically
! keeping the Donelan formula in terms of U10. The 28U* comes from Komen
! et al. (1984). The physical constraint, which is in terms of U*, make
! the resulting scaling behavior less obvious.
!
! Updated Aug 9 2012: When the new logical variable "TRUE_U10" is true,
! SWAN uses U10 instead of 28Ustar in SWIND_DBYB. Ustar is used only
! for the physical constraint in this case.
!
! Updated Feb 7 2014: "28" in the formula U10PROXY=28*UFRIC is made a
! user-defined variable. Note that "28" corresponds to C10=1.28e-3,
! which one might see for U10~5 m/s (pertinent to Snyder measurements
! and therefore used by Komen et al. (1984)) whereas a typical Lake
! George / DBYB value is C10=1.41e-3 to 1.43e-3. So U10=26.5*Ustar,
! and the max value of Cd would give U10~23.5*Ustar. Or, if we go by
! wind speed, a typical U10 during AUSWEX was 11 m/s, so with Hwang Cd,
! that comes out to a factor of 24 to 25. Via numerical experiments, I
! have found that using a factor (which I call "WINDSCALING") smaller
! than 28, and recalibrating the model, gives a stronger tendency to
! overpredict high frequency energy (already a problem mentioned in
! Rogers et al. JTECH 2012). Using WINDSCALING>28 and recalibrating
! the model improves the problem with overprediction of high frequency
! energy, but unfortunately, is making our manipulation of the DBYB
! formula more severe (noting that if we want to use DBYB without
! manipulation, we should use the "TRUE_U10" setting). Using larger factor
! also tends to reduce waveheights during tropical cyclones, e.g.
! maximum SWH during my Ivan simulation at 42040 is reduced by ~60 cm,
! which is a slight improvement. Suprisingly, using the larger factor
! does not make the "tail reduction" physical constraint more active.
! This is because energy in the tail is reduced, so Sin in the tail is
! reduced, so "first guess" of tau_wave is not significantly increased.
!
! To summarize, the two options are as follows:
! a) U10PROXY=WINDSCALING*USTAR (recommended)
! (default WINDSCALING=28)
! b) U10=U10 (activated by using "TRUE_U10" setting)
!
! Additional remarks about the tendency to overpredict high frequency
! energy (particularly with WINDSCALING=28) :
! This is evidenced by tendency to underpredict Tm02 even when
! SWH is overpredicted (in simulations of windsea-only situations). It
! can also be seen in direct comparison of E(f), of course. Another clear
! indicator is the tendency of E/ET to exceed 5, for example at 0.3 Hz
! when U10~12 m/s (see my slides from NRL External Review July 2009).
! E/ET, according to Babanin and Soloviev (1998) should not exceed 5, and
! this tendency of the "real ocean" is confirmed from buoy analyses by
! DW in these slides. Using larger values of p1 p2 also improves the
! overprediction, but not much. For example changing [p1,p2] from [4,4]
! to [8,8] improves the problem less than changing WINDSCALING from 28
! 32 (with recalibration). Essentially, the Sds is like a spring that pulls
! E(f) toward ET(f). The values a1 a2 p1 p2 control the tension in the
! spring. Going from 28 to 32 requires a recalibration, an increase in
! a1 a2, i.e. an increase in the tension of the spring. Changing [p1,p2]
! from [4,4] to [8,8] is another type of increase in the tension of the
! spring. Another piece of evidence of overprediction of high frequency
! energy: it can be shown that the slope of the high-frequency tail in
! the saturation region is between -4 and -5, where we expect it to be -5.
!
! Notes on method of calculating AINV:
!
! 1) Using single point model, U10=12 m/s
! CIRCLE 36 0.0271 10.0 62
! GEN3 BABANIN 6.24E-7 8.74E-6 4.0 4.0 1.2 0.0000 UP VECTAU AGROW
!
! Using Babanin AINV as before
! max(HSKn1),max(HSKn2),max(HS) = 2.29 ; 2.77 ; 2.55
!
! 2) Using AINV=4.0112*((RMSDIR(IS)-0.2) ** (0.81)) ! ER FORMULA
! (used in swancom2 in 2009, but precise origin unknown)
! (tends to give a high estimate of AINV)
! max(HSKn1),max(HSKn2),max(HS) = 2.29 ; 2.77 ; 2.52
!
! 3) Using AINV=(RMSDIR(IS)/0.38) ** (1.0/0.94) ; ! DW FORMULA
! (used in swancom2 in 2009, but precise origin unknown)
! max(HSKn1),max(HSKn2),max(HS) = 2.29 ; 2.77 ; 2.60
!
! 4) Using AINV=(RMSDIR(IS)*3.4)-0.25 ! eyeball from Kuik_vs_AINV.m
! max(HSKn1),max(HSKn2),max(HS) = 2.29 ; 2.77 ; 2.55
!
! 5) Using AINV=(RMSDIR(IS)*3.5333)-0.34862 ! robustfit from
! Kuik_vs_AINV.m
! max(HSKn1),max(HSKn2),max(HS) = 2.29 ; 2.77 ; 2.55
!
! 6) Using AINV=RMSDIR(IS)*2.51 ! from experiments with
! test_DSPR_normal_sigma.m
! (tends to give a low estimate of AINV)
! max(HSKn1),max(HSKn2),max(HS) = 2.29 ; 2.77 ; 2.63
!
! 10. SOURCE
!
!***********************************************************************
!
! Subroutine arguments:
REAL SPCDIR(MDC,6)
REAL SPCSIG(MSC)
REAL THETAW, KWAVE(MSC,ICMAX)
REAL MEMSINA(MDC,MSC,MCGRD), MEMSINB(MDC,MSC,MCGRD)
REAL UFRIC ,AC2(MDC,MSC,MCGRD) , WIND10
REAL CG(MSC,ICMAX)
REAL ZELEN(MCGRD)
LOGICAL ANYWND(MDC)
! Local variables:
INTEGER IDDUM ,ID ,IS
REAL THETA ,SIGMA ,SWINEB, CTW ,STW ,COSDIF
REAL TEMP2 ,UoverC
REAL S_IN(MDC,MSC)
REAL DMAX,AINV
! Ktheta is like D(theta), except max value at each freq is unity
REAL KTHETA(MSC,MDC)
REAL ANAR(MSC),SIGDENS(MSC),SQRTBN(MSC),CINV(MSC)
REAL GAMMAD,GDONEL,TEMP4,WPSI,TEMP5,TEMP6,BN
REAL SIN1D(MSC),SWND,FREQ,STRESS
INTEGER IENT
REAL RMSDIR(MSC)
REAL EDENS2D ! for test calcs
REAL CINV2,CTH,STH
REAL TAUX_linear,TAUY_linear
!NRL REAL TAUX_tmp,TAUY_tmp,ENCHECK
REAL ZE !roughness length
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'SWIND_DBYB')
CTW = COS(THETAW)
STW = SIN(THETAW)
! Calculate 1d spectrum E(sigma)
DO IS = 1, MSC
SIGDENS(IS) = 0.
DO ID = 1, MDC
SIGDENS(IS) = SIGDENS(IS) + SPCSIG(IS) * AC2(ID,IS,KCGRD(1))
KTHETA(IS,ID)=AC2(ID,IS,KCGRD(1))
END DO
SIGDENS(IS)=SIGDENS(IS)*DDIR ! units m^2/(radHz)
END DO
! Calculate K(sigma,theta)
DO IS = 1, MSC
DMAX=-1.0
DO ID = 1, MDC
IF(KTHETA(IS,ID).GT.DMAX)DMAX=KTHETA(IS,ID)
END DO
IF(DMAX.EQ.0.0)THEN
! a fix for freq bins (usually first two or so) that are empty
DMAX=1.0
DO ID = 1, MDC
KTHETA(IS,ID)=1.0
END DO
ENDIF
! Optional: Calculate circular RMS directional spread (for comparison with
! AINV) :
! rmsdir(is)=kuik(ktheta(is,:),ddir,spcdir(:,2),spcdir(:,3),.true./.false)
DO ID = 1, MDC
KTHETA(IS,ID)=KTHETA(IS,ID)/DMAX
END DO
END DO
! Calculate SQRT(Bn(IS))
DO IS=1,MSC
AINV=0.0
DO ID = 1, MDC
AINV=AINV+KTHETA(IS,ID)*DDIR
END DO
! Optional: use RMSDIR in place of AINV: to see how to do this, read notes in
! rev 622 and prior
ANAR(IS)=1.0/AINV
BN=ANAR(IS)*SIGDENS(IS)*(KWAVE(IS,1)**3)*CG(IS,1)
SQRTBN(IS)=SQRT(BN)
END DO
! Calculate S_in for *ALL* directions (not just current sweep). Though this
! is somewhat wasteful (S_in will end up being calculated 4 times for
! each bin), it is nececessary in order to have the correct stress going
! into CALC_LFACTOR
! Update May 2017: Above issue was addressed by Marcel Z. via "memsin".
! Note, if this setting for TEMP2 is modified, it should also be modified in
! subroutine "CALC_TAU_TOTAL"
IF(TRUE_U10)THEN
TEMP2 = WIND10 ! "TEMP2=WIND10" is basically saying "use DBYB as is"
ELSE
TEMP2 = WNDSCL * UFRIC
ENDIF
S_IN=0.0
DO IS = 1, MSC
SIGMA = SPCSIG(IS)
CINV(IS) = KWAVE(IS,1) / SIGMA
UOVERC = TEMP2 * CINV(IS)
DO ID=1,MDC
IF ( ANYWND(ID) ) THEN
COSDIF = SPCDIR(ID,2)*CTW + SPCDIR(ID,3)*STW
!-------------------------------------------------------------------------------
! Yalin Fan, 08/06/13: Below equations calculate W(f,theta) in eq6 of Rogers et al (2012)
! This equation limits the grow rate to be positve, while Tolman and Chalikov (1996)
! gives negative growth rate when wind and wave have large angle
! If the 0 constraint is removed and the sign of TEMP4 is maintained, then W(f,theta)
! is very similar to Beta in Tolman and Chalikov (1996)
!-------------------------------------------------------------------------------
TEMP4=( UOVERC * COSDIF - 1.0 )
TEMP4=MAX(0.0,TEMP4)
WPSI=TEMP4**2
TEMP5=10.0*SQRTBN(IS)*WPSI-11.0
TEMP6=1.0+TANH(TEMP5)
GDONEL=2.8-TEMP6
GAMMAD=GDONEL*SQRTBN(IS)*WPSI
SWINEB = GAMMAD * SIGMA * PWIND(9) ! NEW SWINEB
! In Donelan notation, SWINEB is BETA, so I use BETA=GAMMA*sigma*rhoa/rhow
! (Donelan eq 3)
! Calculate actual wind input term Sin=Beta*Edens
S_IN(ID,IS)=SWINEB*AC2(ID,IS,KCGRD(1))*SPCSIG(IS)
! Note that IMATRA calculation will be done after reduction operation
END IF
ENDDO
ENDDO
! calculate stress for the linear wind input part
TAUX_linear=0.0
TAUY_linear=0.0
DO IS = 1, MSC
CINV2=KWAVE(IS,1)/SPCSIG(IS)
DO ID = 1, MDC
CTH = SPCDIR(ID,2) ! cos(theta)
STH = SPCDIR(ID,3) ! sin(theta)
TAUX_linear=TAUX_linear+CTH*MEMSINA(ID,IS,KCGRD(1))*CINV2* SPCSIG(IS)**2*FRINTF*DDIR
TAUY_linear=TAUY_linear+STH*MEMSINA(ID,IS,KCGRD(1))*CINV2* SPCSIG(IS)**2*FRINTF*DDIR
END DO
END DO
TAUX_linear=TAUX_linear*PWIND(17)*GRAV
TAUY_linear=TAUY_linear*PWIND(17)*GRAV
! Note: test output that was here has been deleted. See rev 622 or earlier.
CALL CALC_LFACTOR(TAUX_linear,TAUY_linear,RDFSIN,S_IN,UFRIC,PWIND,DDIR,SPCSIG,FRINTF,CINV,GRAV, &
WIND10,TESTFL,SPCDIR,VECTOR_TAU,TRUE_U10,CTW,STW,ZE)
ZELEN(KCGRD(1)) = ZE
!-------------------------------------------------------------------------------
! Begin negative wind input
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
! Yalin, Aug30, 2013. Add negative input when wind and wave have large angle
! RDCOEF is the user adjustable reduction coefficient for the negative energy.
! RDCOEF = 0.4 is the suggested value. It is specified in INPUT as "NEGATINP 0.4"
! If no keyword "NEGATINP" specified, RDCOEF is default to 0.0
! ------------------------------------------------------------------------------
! WER, June 8, 2017. "0.4" is from Donelan (1999, Book Chapter, "Wind-induced growth
! and attenuation of laboratory waves") (See also discussion in Donelan et al. (2006).)
! Zieger et al. (2015) use reduction_coef = 0.04 in conjunction with the WW3 v4.18
! swell dissipation (see their table 1). (At time of writing, the WW3 v5.08 swell
! dissipation mentioned in the same table is not implemented in SWAN)
! ------------------------------------------------------------------------------
! WER, June 9, 2017. For efficiency, only go into this loop if ZIEGER. (Changed)
! ------------------------------------------------------------------------------
! WER, Jun 12 2017: I confirmed that use of ANYWND here is a bug, since it disables
! negative wind input for case of component in the +/- 90 deg arc opposing the wind.
! I have removed the ANYWND operation.
IF(ZIEGER)THEN
DO IS = 1, MSC
SIGMA = SPCSIG(IS)
CINV(IS) = KWAVE(IS,1) / SIGMA
UoverC = TEMP2 * CINV(IS)
DO ID=1,MDC
COSDIF = SPCDIR(ID,2)*CTW + SPCDIR(ID,3)*STW
TEMP4 = ( UoverC * COSDIF - 1.0 ) ! modify Rogers_etal_2012 eq 6
TEMP4 = MIN(0.0,TEMP4) ! modify Rogers_etal_2012 eq 6
WPSI = TEMP4**2 ! modify Rogers_etal_2012 eq 6
TEMP5 = 10.0*SQRTBN(IS)*WPSI-11.0
TEMP6 = 1.0+TANH(TEMP5)
GDONEL = 2.8-TEMP6
GAMMAD = GDONEL*SQRTBN(IS)*WPSI
SWINEB = GAMMAD * SIGMA * PWIND(9)
S_IN(ID,IS)= S_IN(ID,IS) - SWINEB*AC2(ID,IS,KCGRD(1))*SPCSIG(IS)*RDCOEF
ENDDO
ENDDO
END IF
!-------------------------------------------------------------------------------
! End negative wind input
!-------------------------------------------------------------------------------
! We want to add B*N to RHS. This is SWINEB*AC2, see SWIND3.
! Above, we had S_IN(ID,IS)=SWINEB*AC2(ID,IS,KCGRD(1))*SPCSIG(IS)
! Thus, we just need to divide out SPCSIG(IS).
! Store this wind input in MEMSINB array for every grid point,
! which can be filled in IMATRA in next four sweeps
! (see subroutine FILSIN).
! Also note that this part, B*N, is explicit and stored in PLWNDS,
! which represents total wind input (=A+BE)
! (see subroutine FILSIN)
DO IS = 1, MSC
DO ID = 1, MDC
S_IN(ID,IS)= S_IN(ID,IS) * RDFSIN(IS)
MEMSINB(ID,IS,KCGRD(1)) = S_IN(ID,IS) / SPCSIG(IS)
!ESMF IF (SAVE_SINBAC) SINBAC(ID,IS,KCGRD(1)) = &
!ESMF SINBAC(ID,IS,KCGRD(1)) + S_IN(ID,IS)/SPCSIG(IS)
ENDDO
ENDDO
! calculate stress (test point only)
!NRL IF(TESTFL)THEN
!NRL TAUX_tmp=0.0
!NRL TAUY_tmp=0.0
!NRL ENCHECK=0.0
!NRL DO IS = 1, MSC
!NRL CINV2=KWAVE(IS,1)/SPCSIG(IS)
!NRL DO ID = 1, MDC
!NRL CTH = SPCDIR(ID,2) ! cos(theta)
!NRL STH = SPCDIR(ID,3) ! sin(theta)
!NRL TAUX_tmp=TAUX_tmp+CTH*CINV2*S_IN(ID,IS)*DDIR*FRINTF*SPCSIG(IS)
!NRL TAUY_tmp=TAUY_tmp+STH*CINV2*S_IN(ID,IS)*DDIR*FRINTF*SPCSIG(IS)
!NRL ENCHECK=ENCHECK+AC2(ID,IS,KCGRD(1))*DDIR*FRINTF*SPCSIG(IS)**2
!NRL END DO
!NRL END DO
!NRL TAUX_tmp=TAUX_tmp*PWIND(17)*GRAV
!NRL TAUY_tmp=TAUY_tmp*PWIND(17)*GRAV
!NRL WRITE(*,*)'SWIND: TAUX,TAUY,HM0 = ',TAUX_tmp,TAUY_tmp,(4*SQRT(ENCHECK))
!NRL ENDIF
! *** test output ***
IF (ITEST.GE. 80.AND.TESTFL) THEN
WRITE(PRTEST,6000) KCGRD(1), THETAW*180./PI
6000 FORMAT(' SWIND_DBYB: POINT THETAW :',I5,E12.4)
WRITE(PRTEST,6100) TEMP2, UFRIC
6100 FORMAT(' SWIND_DBYB: TEMP2 UFRC :',3E12.4, /, &
' IS ID1 ID2 Wind source term')
DO IS = 1, MSC
WRITE(PRTEST,6200) IS, 1, MDC,(MEMSINB(ID,IS,KCGRD(1)), &
ID=1,MDC)
6200 FORMAT(3I4, 600e12.4)
ENDDO
WRITE(PRTEST,*)
END IF
RETURN
END SUBROUTINE SWIND_DBYB
SUBROUTINE CALC_LFACTOR(TAUX_linear,TAUY_linear,LFACTOR_L,S_IN,UFRIC,PWIND,DDIR_RAD,SIGMA_S,FRINTF, &
CINV_S,GRAV,WIND10,TESTFL,SPCDIR,VECTOR_TAU,TRUE_U10,CTHETA_WIND, &
STHETA_WIND, ZE)
!
USE SWCOMM1, ONLY: CHTIME
!
IMPLICIT NONE
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2024 Delft University of Technology
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program. If not, see <http://www.gnu.org/licenses/>.
!
!
! Objective : provide Lfactor (1-d array) by solving for optimal REDUC
! (scalar). It uses a solver method that should work for all monotonic
! relations. For non-monotonic relations, it may get confused. Here,
! tau_normal definitely has a monotonic dependence on REDUC, so this
! potential limitation on the solver is OK. (Monotonicity of S_in, or lack
! thereof, is irrelevant.)
! Points made here:
! 1) The normal stress plus the tangential/viscous stress cannot be greater
! than the total stress.
! 2) We do not know what the contribution is beyond fmax, but we know that it
! is not negative, so again, we say that the normal stress in our
! prognostic range cannot be greater than the total stress...so
! tau_total=tau_normal_prognostic+tau_normal_diagnostic+tau_tangential...
! all are >=0, so maximum allowed value of tau_normal_prognostic is
! tau_total. Update: in this code, we partially address this issue by
! extending the integration (fmax) to 10 Hz for the purposes of this
! integration.
! 3) Similar to Kakha's argument, we use the exponential adjustment that is
! stronger at the higher frequencies (more reduction) since this is where
! our formula is less well-informed by observations (more wiggle room),
! since typical observations are for the dominant waves, f=fp
! 4) Kakha uses f/fp to scale the reduction, which makes sense in the context
! of item (3), but use of fp has disadvantages, so we use U/C instead.
! This means that for young wind sea, S_in for the entire spectrum is
! reduced...will this be a problem? It is probably worthwhile to try using
! an fp value defined in a manner similar to Tolman and Chalikov, which
! does not have the traditional problems with using fp.
! 5) S_in=S_in_linear+S_in_exponential....same argument as (1). Linear wind
! input is assumed small and not included in the stress calculation. This
! was in fact a major source of error in some cases, but I have limited
! the linear wind input source function so that it cannot contribute more
! than 1% of the total stress now. Thus, the error associated with
! excluding it from the stress calculation here will be less than 1% now.
! Having said that, it would be useful to include it, as a future code
! improvement. (Update: Yalin has added linear part to stress calculation.
! This was rev 669, Feb 26 2014.)
! Notes re: "vector_tau":
! In earlier versions of this code, to simplify calcs, I used S_in1D(f)
! S_in1D_S=sum(S_in,1)*DDIR_rad*(2.0*PI) ! in units m2/Hz
! this simplication implies that we will follow the Tsagareli et al.
! (2010 JPO) method for calculating stress:
! tau=integral (Sin(f)/C(f) df
! This calculation is not so good, since it assumes all stresses are in
! the same direction and so never counteract each other.
! That is the method reported in Rogers et al. (JTECH 2012)
! Updated August 2012:
! We replace the non-directional
! tau=integral of Sin(f)/C(f) df
! with
! tau_x=integral of Sin(f,theta)*k_x/(k*C(f)) df dtheta
! tau_y=integral of Sin(f,theta)*k_y/(k*C(f)) df dtheta
! essentially using the unit vector <k> = k_vector/k_scalar
! and it can be written also
! tau_vector=integral of (k_vector/omega) Sin(f,theta) df dtheta
! Now, if vector_tau is true, we use a vector calculation.
! If vector_tau is false, we do not need the directional Sin(f,theta).
! However, to keep the code simple I pass Sin(f,theta) to calc_tau_total
! routine regardless.
! Note: "_S" denotes "short" (fewer freqs)
! "_L" denotes "long" (freqs out to 10 Hz)
! (1d version: was used for calculations, but now used only
! for diagnostic output)
! (2d version: added and now used for calculations)
REAL , INTENT(OUT) :: LFACTOR_L(:)
REAL , INTENT(OUT) :: ZE
REAL , INTENT(IN) :: S_IN(:,:)
REAL , INTENT(IN) :: SPCDIR(:,:)
REAL , INTENT(IN) :: UFRIC
REAL , INTENT(IN) :: DDIR_RAD
REAL , INTENT(IN) :: FRINTF
REAL , INTENT(IN) :: PWIND(:)
REAL , INTENT(IN) :: CINV_S(:)
REAL , INTENT(IN) :: SIGMA_S(:)
REAL , INTENT(IN) :: GRAV
REAL , INTENT(IN) :: TAUX_linear,TAUY_linear,WIND10,CTHETA_WIND,STHETA_WIND
LOGICAL , INTENT(IN) :: TESTFL
LOGICAL , INTENT(IN) :: VECTOR_TAU,TRUE_U10
REAL , ALLOCATABLE :: S_IN1D_L(:)
REAL , ALLOCATABLE :: S_IN_L(:,:)
REAL , ALLOCATABLE :: DF(:)
REAL , ALLOCATABLE :: CINV_L(:)
REAL , ALLOCATABLE :: SIGMA_L(:)
REAL :: DSIGMA
REAL :: REDUC
REAL :: SIGN_NEW,SIGN_OLD
REAL :: FRAT
REAL :: ERR
REAL :: RCHANGE
REAL :: Z0
REAL :: AFA
! tau_total_max is the maximum allowed value for total wind stress
! tau_total is the value for the total wind stress computed from tau_wave and
! tau_visc (i.e. tau_normal and tau_tangential)
! tau_sin_x,tau_sin_y are tau contribution from wind input: in this routine,
! it is only used for test output
REAL :: TAU_TOTAL, TAU_wav
REAL :: TAU_TOTAL_MAX
REAL :: TAU_SIN_X,TAU_SIN_Y
REAL :: DDIR_DEG
REAL :: RHOA,RHOW,PI
REAL :: FREQ_TMP
REAL :: U10PROXY,TAU_VISC,TAU_VISC_X,TAU_VISC_Y
INTEGER :: NF_OLD,NF_NEW,NDIR
INTEGER :: DIR
INTEGER :: ITER
INTEGER :: SLOW_DOWN
INTEGER :: ISOL
INTEGER :: IS,IDIR
!...S_in(theta,freq)
NDIR=SIZE(S_IN,1)
NF_OLD=SIZE(S_IN,2)
PI=2.0*ACOS(0.0)
AFA = 0.01
!...find nf_new
FRAT=FRINTF+1.0 ! =freq(nf_old)/freq(nf_old-1)
FREQ_TMP=SIGMA_S(NF_OLD)/(2.0*PI)
NF_NEW=-99
DO IS=(NF_OLD+1),(NF_OLD+100)
FREQ_TMP=FREQ_TMP*FRAT
IF(FREQ_TMP > 10.0)THEN ! 10.0 here is f=10 Hz
NF_NEW=IS
EXIT
ENDIF
ENDDO
! shortcoming of this method: the maximum frequency changes slightly from
! one simulation to the next, depending on the frequency bands the user
! selects. To correct this, we should stop at 10 Hz and change df of last bin
! df for the last bin would have to use finite differencing, instead of
! FRINTF (as used below), e.g. DF=FREQ(nf_new)-FREQ(nf_new-1)
!... allocate arrays on nf_new
ALLOCATE(S_IN_L(NDIR,NF_NEW))
ALLOCATE(S_IN1D_L(NF_NEW)) ! for diagnostic output only
ALLOCATE(DF(NF_NEW))
ALLOCATE(CINV_L(NF_NEW))
ALLOCATE(SIGMA_L(NF_NEW))
!... extend freqs to nf_new
SIGMA_L(1:NF_OLD)=SIGMA_S(1:NF_OLD)
CINV_L(1:NF_OLD)=CINV_S(1:NF_OLD)
S_IN_L(:,1:NF_OLD)=S_IN(:,1:NF_OLD)
DO IS=(NF_OLD+1),NF_NEW
SIGMA_L(IS)=SIGMA_L(IS-1)*FRAT
CINV_L(IS)=SIGMA_L(IS)*0.102 ! deep water assumption for hf tail is OK
!... For Sin, we use an f-2 approximation
!... ...Sin=Beta*E~sigma*(U/C)^2*E = f * f^2 * f^-5 = f^-2
DO IDIR=1,NDIR
S_IN_L(IDIR,IS)=S_IN_L(IDIR,NF_OLD)*(SIGMA_L(NF_OLD)/SIGMA_L(IS))**2
ENDDO
ENDDO
S_IN_L=S_IN_L*(2.0*PI) ! was in units of m2/(rad-Hz)/rad, so put in
! units m2/Hz/rad
S_IN1D_L=SUM(S_IN_L,1)*DDIR_RAD ! note that variable is S_in(ID,IS)
!NRL IF(TESTFL)THEN
!NRL WRITE(416,*)CHTIME,' % CHTIME'
!NRL WRITE(416,*)' % (freq , old S_in)'
!NRL DO IS=1,NF_NEW
!NRL WRITE(416,211)(SIGMA_L(IS)/(2.0*PI)),S_IN1D_L(IS)
!NRL END DO
!NRL ENDIF
S_IN1D_L=0.0
RHOA=PWIND(16)
RHOW=PWIND(17)
TAU_TOTAL_MAX=(UFRIC**2)*RHOA
DDIR_DEG=DDIR_RAD*180.0/PI
! Method: Tsagareli (thesis, 2009); Tsagareli et al. (JPO 2010, eq 3.7);
! based on fitting to data of Banner and Peirson (JFM 1998).
! New Aug 30 2013: Previously, tau_visc would go to zero for high wind speeds
! (U10>22 m/s), which looks strange when plotted. This is addressed now by
! reformulating such that tau_visc remains constant beyond the wind speed
! at which tau_visc is maximum, i.e. U10 at which dtau/dU=0. This is 14.667
! m/s, which also happens to be close to the maximum wind speed in the plot
! of Banner and Peirson, which Tsagareli used. So we are also holding
! tau_visc constant for wind speeds for which we have no guidance from
! the observations.
U10proxy=min(WIND10,14.66666666)
tau_visc=(1.408e-3)*(U10proxy**2) - (6.4e-5)*(U10proxy**3)
! Notes, Aug 15 2012: Stefan noticed that at low wind speeds, |Cv| is
! slightly larger than |Cd|. This probably is not reasonable. To address
! this, we say "tau_visc cannot be more than 90% of tau_total_max,
! regardless". The 90% is semi-arbitrary, selected because it implies that
! this "fix" has a very minor impact on computations, relative to, say,
! 60%.