! 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-2017 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 2 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.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 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-2017 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 2 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.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 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-2017 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 2 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.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 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%.
tau_visc=min((0.9*tau_total_max),tau_visc)
! Note, Aug 30 2013: We may consider using a formulation based on "Law of the
! Wall", (the subset that deals with hydraulically smooth flow). This is
! used by Makin et al. (BLM 1995), Donelan et al. (JGR 2012), Edson et al
! (JPO 2013), Hayes et al. (Icarus 2013), and George Mellor, among
! others. The method results in tau_visc very similar to that of our
! existing method in cases of low winds (U10<3 m/s) which is the regime
! where tau_visc is expected to be important. However, it is very different
! at high wind speeds, since it grows indefinitely with U10. The result is
! that skin drag (tau_visc) is greater than 28% of total stress when U10
! is greater than 25 m/s. See my matlab script Makin_iter.m. If we do this,
! we may want to reduce tau_visc to account for sheltering. Donelan et al.
! (JPO 2012) do this, though they use Cd, and I think using mss would be
! better. But even if we reduce tau_visc to 15% of total stress in high
! winds, this still implies a 15% reduction of wind stress, which would be
! similar to a 8% reduction in Ustar for our hurricane modeling vs. our
! current approach based on Banner and Peirson.
! Example computations for U10=40 m/s:
! tau_visc as implemented prior to Aug 30 2013 : zero
! tau_visc as implemented after Aug 30 2013: 0.1 Pa
! tau_visc from Makin_iter.m : 1.4 Pa
! tau_total: 4.3 Pa
! We will assume that viscous stress is in wind direction
TAU_VISC_X=TAU_VISC*CTHETA_WIND
TAU_VISC_Y=TAU_VISC*STHETA_WIND
! DF is needed by calc_tau_total
DO IS=1,NF_NEW
DSIGMA=SIGMA_L(IS)*FRINTF ! since DF=FREQ*FRINTF
DF(IS)=DSIGMA/(2.0*PI)
END DO
! calculate without reduction
REDUC=0.0
CALL CALC_TAU_TOTAL(TAU_TOTAL,LFACTOR_L,REDUC,S_IN_L,DF,CINV_L,WIND10,GRAV,&
RHOW,SPCDIR,DDIR_RAD,VECTOR_TAU,TRUE_U10,UFRIC,TAU_VISC_X,TAU_VISC_Y, &
TAU_SIN_X,TAU_SIN_Y,TAUX_linear,TAUY_linear,SIGMA_L,TESTFL)
206 FORMAT('tau_total_max = ',F9.6,' ; tau_total = ',F10.6,' ; tau_visc = ',&
F9.6)
!NRL IF(TESTFL)THEN
!NRL WRITE(*,206)TAU_TOTAL_MAX,TAU_TOTAL,TAU_VISC
!NRL ENDIF
IF(TAU_TOTAL < TAU_TOTAL_MAX .OR. TAU_TOTAL < 1E-10)THEN
DO IS=1,NF_OLD
LFACTOR_L(IS)=1.0
END DO
DO IS=NF_OLD+1,NF_NEW
LFACTOR_L(IS)=0.0
END DO
! Deallocate and return
DEALLOCATE(S_IN_L)
DEALLOCATE(S_IN1D_L)
DEALLOCATE(DF)
DEALLOCATE(CINV_L)
DEALLOCATE(SIGMA_L)
RETURN
ELSE
!NRL IF(TESTFL)THEN
!NRL WRITE(*,*)'REDUCING SIN TO MAKE TAU_TOTAL = TAU_TOTAL_MAX'
!NRL ENDIF
ENDIF
! test output
! e.g. 132+ for NGOM Katrina case, 502+ for GMEX HRD_buoy Ivan case...
! note that if Edens is too large due to weak Sds, this will be large as
! a result
!NRL IF(TESTFL)THEN
!NRL IF(TAU_TOTAL > 500.0)THEN
!NRL WRITE(*,*)'WARNING: tau_total without reduction seems large: ', &
!NRL TAU_TOTAL
!NRL ENDIF
!NRL ENDIF
! calculate *with* reduction, start value arbitrary
REDUC=0.05
CALL CALC_TAU_TOTAL(TAU_TOTAL,LFACTOR_L,REDUC,S_IN_L,DF,CINV_L,WIND10,GRAV,&
RHOW,SPCDIR,DDIR_RAD,VECTOR_TAU,TRUE_U10,UFRIC,TAU_VISC_X,TAU_VISC_Y,&
TAU_SIN_X,TAU_SIN_Y,TAUX_linear,TAUY_linear,SIGMA_L,TESTFL)
ERR=TAU_TOTAL-TAU_TOTAL_MAX
SIGN_NEW=SIGN(1.0,ERR)
!NRL IF(TESTFL)THEN
!NRL WRITE(*,207)0,REDUC,0.0,TAU_TOTAL,ERR
!NRL ENDIF
SLOW_DOWN=0
RCHANGE=2.0
ISOL=0
! Notes: Objective: find value of REDUC that yields desired tau_total.
! Current method is to change REDUC sequentially, and slow down the
! incrementing of REDUC if there is a sign change in the err value.
! [incrementing of REDUC=RCHANGE]. It might be useful to experiment with
! Newton-Raphson as an alternative
DO ITER=1,50
IF(TAU_TOTAL > TAU_TOTAL_MAX)THEN ! increase REDUC
REDUC=REDUC*RCHANGE
ELSE ! then reduce REDUC
REDUC=REDUC/RCHANGE
ENDIF
CALL CALC_TAU_TOTAL(TAU_TOTAL,LFACTOR_L,REDUC,S_IN_L,DF,CINV_L,WIND10, &
GRAV,RHOW,SPCDIR,DDIR_RAD,VECTOR_TAU,TRUE_U10,UFRIC,TAU_VISC_X, &
TAU_VISC_Y,TAU_SIN_X,TAU_SIN_Y,TAUX_linear,TAUY_linear,SIGMA_L,TESTFL)
ERR=TAU_TOTAL-TAU_TOTAL_MAX
SIGN_OLD=SIGN_NEW
SIGN_NEW=SIGN(1.0,ERR)
207 FORMAT('iter = ',I3,' ; REDUC = ',F9.6,' ; RCHANGE = ',F9.6, &
' ; tau_total = ', F9.6,' ; err = ',F7.4,' ; tau_sin_x = ',F9.6, &
' ; tau_sin_y = ',F9.6)
!NRL IF(TESTFL)THEN
!NRL WRITE(*,207)ITER,REDUC,RCHANGE,TAU_TOTAL,ERR,TAU_SIN_X,TAU_SIN_Y
!NRL ENDIF
! once "slow_down" is set for a given iteration, it stays in this state
! for future iterations
IF(SIGN_NEW /= SIGN_OLD)THEN
SLOW_DOWN=1
ENDIF
IF(SLOW_DOWN==1)THEN
RCHANGE=0.5*(1.0+RCHANGE)
ENDIF
! For notes re: alternative method that was experimented with,
! see rev 622 and prior.
IF((ABS(ERR)/TAU_TOTAL_MAX) < 5.E-04)THEN
ISOL=1
EXIT
ENDIF
END DO ! DO ITER=...
!! BEGIN: force stop in case of failure
IF(ISOL==0)THEN
IF((ABS(ERR)/TAU_TOTAL_MAX).GE.0.001)THEN
WRITE(*,*)'warning: no solution found'
WRITE(*,*)' tau_total = ',TAU_TOTAL,' tau_total_max = ',TAU_TOTAL_MAX, &
' err = ',ERR,'(abs(err)/tau_total_max) = ',(ABS(ERR)/TAU_TOTAL_MAX)
ENDIF
!NRL IF((ABS(ERR)/TAU_TOTAL_MAX).GE.1.0)THEN
!NRL OPEN(411,FILE='DEBUG.DAT',FORM='UNFORMATTED')
!NRL WRITE(411)NF_NEW
!NRL WRITE(411)LFACTOR_L ! LFACTOR_L(nf_new)
!NRL WRITE(411)S_IN_L ! S_in_L(ndir,nf_new)
!NRL WRITE(411)DF ! DF(nf_new)
!NRL WRITE(411)CINV_L ! CINV_L(nf_new)
!NRL WRITE(411)TAU_TOTAL_MAX
!NRL WRITE(411)TAU_TOTAL
!NRL WRITE(411)REDUC
!NRL WRITE(411)WIND10
!NRL WRITE(411)RHOW
!NRL CLOSE(411)
!NRL WRITE(*,*)'stopping, see DEBUG.DAT'
!NRL STOP
!NRL END IF
IF(TAU_TOTAL.NE.TAU_TOTAL)THEN ! check for NaN
WRITE(*,*)'stopping due to tau_total=NaN'
STOP
ENDIF
END IF
!! END: force stop in case of failure
!! BEGIN: test output for NRL purposes
!NRL210 FORMAT(F9.6,' ',E12.6,' % (freq , Lfactor)')
!NRL IF(TESTFL)THEN
!NRL DO IS=1,NF_NEW
!NRL WRITE(413,210)(SIGMA_L(IS)/(2.0*PI)),LFACTOR_L(IS)
!NRL END DO
!NRL ENDIF
S_IN1D_L=SUM(S_IN_L,1)*DDIR_RAD
!NRL211 FORMAT(F9.6,1X,E12.6)
!NRL IF(TESTFL)THEN
!NRL WRITE(414,*)CHTIME,' % CHTIME'
!NRL WRITE(414,*)' % (freq , new S_in)'
!NRL DO IS=1,NF_NEW
!NRL WRITE(414,211)(SIGMA_L(IS)/(2.0*PI)),(S_IN1D_L(IS)*LFACTOR_L(IS))
!NRL END DO
!NRL ENDIF
! S_IN1D_L=0.0
!! END: test output
TAU_wav=0.0
DO IS=1,NF_NEW
TAU_wav=TAU_wav+RHOW*GRAV*S_IN1D_L(IS)*LFACTOR_L(IS) &
*CINV_L(IS)*DF(IS)
END DO
Z0 = AFA * UFRIC**2 / GRAV
ZE = Z0/SQRT(1-TAU_wav/TAU_TOTAL)
! Deallocate and return
DEALLOCATE(S_IN_L)
DEALLOCATE(S_IN1D_L)
DEALLOCATE(DF)
DEALLOCATE(CINV_L)
DEALLOCATE(SIGMA_L)
END SUBROUTINE CALC_LFACTOR
SUBROUTINE CALC_TAU_TOTAL(TAU_TOTAL,LFACTOR,REDUC,S_IN,DF,CINV,WIND10,GRAV, &
RHOW,SPCDIR,DDIR_RAD,VECTOR_TAU,TRUE_U10,UFRIC,TAU_VISC_X,TAU_VISC_Y, &
TAU_SIN_X,TAU_SIN_Y,TAUX_l,TAUY_l,SIGMA,TESTFL)
! Objective : apply exponential reduction...more reduction at higher U/C
! with factor smoothly intersecting value=1 at U/C=1. It is not applied
! for U/C<1 since that would produce an increase...so factor 1
! for U/C<=1
! On methods used (Aug 14 2012):
! Old method (reasonable enough if tau_wave, tau_total, and tau_visc are
! all in the same direction) (does not make much sense if tau_wave is
! opposite to wind direction, which can be the case if we have negative
! wind input)
! * tau_total=(UFRIC**2)*rhoa
! * tau_wave_maximum=tau_total-tauv;
! * pass Sin(f,theta) and reduction factor to subroutine
! calc_tau_normal
! * subroutine calc_tau_normal integrates Sin with reduction factor
! and returns tau_normal
! * repeat until tau_normal is less than tau_wave_maximum
! New method:
! * there is no tau_wave_maximum
! * tau_total_maximum=(UFRIC**2)*rhoa (this is a scalar)
! * tau_visc is now a vector, in the wind direction
! * pass Sin(f,theta), reduction factor, x-component of tauv,
! y-component of tauv to subroutine calc_tau_total
! * subroutine calc_tau_total integrates Sin with reduction factor,
! plus tau_visc components and returns tau_total (a scalar)
! * repeat until tau_total (a scalar) is less than tau_total_maximum
! (a scalar)
! On computational efficiency, David Johnson, March 1 2013, writes: "... I
! profiled your SWAN code and found out that fully 50% of the total run
! time is spent in subroutine calc_tau_total, with an effective doubling
! of runtimes with the new physics - much of this is probably the
! repeated subroutine calls and memalloc etc.. By inlining the shear
! stress calculation within calc_Lfactor, I managed to reduce run time
! by about 40%. " <end quote> I have experimented with this and do not
! find any improvement to run time (it actually increases run time by
! 1%).Since calc_tau_total is called in 3 separate places, inlining
! makes the code itself less concise, less readable and more prone to
! bugs. Therefore, I keep calc_tau_total as a separate subroutine for
! now.
! On run time:
! I find that for a square domain, deep water, run in serial mode,
! there is the following increase in computation time (expressed as a
! factor between total run time with new physics and total run time
! with old physics):
! domain size 7x7: factor 1.74
! domain size 101x101: factor 1.27 (with substantial i/o, this drops
! to 1.21)
USE SWCOMM3, ONLY: WNDSCL
IMPLICIT NONE
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 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 2 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.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
REAL , INTENT(OUT) :: TAU_TOTAL
REAL , INTENT(OUT) :: LFACTOR(:)
REAL , INTENT(IN) :: REDUC
REAL , INTENT(IN) :: S_IN(:,:)
REAL , INTENT(IN) :: SIGMA(:)
REAL , INTENT(IN) :: DF(:)
REAL , INTENT(IN) :: CINV(:)
REAL , INTENT(IN) :: WIND10,UFRIC
REAL , INTENT(IN) :: GRAV
REAL , INTENT(IN) :: RHOW
REAL , INTENT(IN) :: DDIR_RAD
REAL , INTENT(IN) :: TAU_VISC_X,TAU_VISC_Y,TAUX_l,TAUY_l
REAL , INTENT(IN) :: SPCDIR(:,:)
REAL , INTENT(OUT) :: TAU_SIN_X,TAU_SIN_Y
! SPCDIR explained:
! SPCDIR: (*,1); spectral directions (radians)
! ECOS : =SPCDIR(*,2); cosine of spectral directions
! ESIN : =SPCDIR(*,3); sine of spectral directions
LOGICAL , INTENT(IN) :: VECTOR_TAU, TRUE_U10, TESTFL
REAL , ALLOCATABLE :: AEXP(:)
REAL , ALLOCATABLE :: UOVERC(:)
REAL , ALLOCATABLE :: S_IN_RED(:,:)
INTEGER :: NF,NDIR,IS,IDIR
LOGICAL :: IS_THIS_A_NAN
REAL :: PI,CTHETA,STHETA,TAU_TOTAL_X,TAU_TOTAL_Y,TAU_SIN_SCALAR,TEMP2,TAU_TMP
NF=SIZE(DF)
NDIR=SIZE(S_IN,1)
PI=2.0*ACOS(0.0)
ALLOCATE(UOVERC(NF))
ALLOCATE(AEXP(NF))
ALLOCATE(S_IN_RED(NDIR,NF))
IF(TRUE_U10)THEN
TEMP2 = WIND10
ELSE
TEMP2 = WNDSCL * UFRIC
ENDIF
UOVERC=TEMP2 * CINV ! matrix op
AEXP=EXP((1-UOVERC) * REDUC) ! MATRIX OP
DO IS=1,NF
LFACTOR(IS)=MIN(1.0,AEXP(IS))
DO IDIR=1,NDIR
S_IN_RED(IDIR,IS)=S_IN(IDIR,IS) * LFACTOR(IS)
ENDDO
ENDDO
! see notes above about this calculation
TAU_SIN_X=0.0
TAU_SIN_Y=0.0
TAU_SIN_SCALAR=0.0
DO IS=1,NF
DO IDIR=1,NDIR
CTHETA=SPCDIR(IDIR,2)
STHETA=SPCDIR(IDIR,3)
TAU_SIN_SCALAR=TAU_SIN_SCALAR+RHOW*GRAV*S_IN_RED(IDIR,IS)*CINV(IS) &
*DF(IS)*DDIR_RAD
TAU_SIN_X=TAU_SIN_X+RHOW*GRAV*S_IN_RED(IDIR,IS)*CINV(IS)*DF(IS) &
*CTHETA*DDIR_RAD
TAU_SIN_Y=TAU_SIN_Y+RHOW*GRAV*S_IN_RED(IDIR,IS)*CINV(IS)*DF(IS) &
*STHETA*DDIR_RAD
END DO
!NRL! IF(TESTFL)THEN
!NRL! output stress in cumulative-with-freq format
!NRL! this generates a large file, so use only when needed
!NRL! TAU_TMP=SQRT(TAU_SIN_X**2+TAU_SIN_Y**2)
!NRL! WRITE(417,*)IS,(SIGMA(IS)/(2.0*PI)),TAU_TMP
!NRL! END IF
ENDDO
IF(VECTOR_TAU)THEN
TAU_TOTAL_X=TAU_SIN_X+TAU_VISC_X+TAUX_l
TAU_TOTAL_Y=TAU_SIN_Y+TAU_VISC_Y+TAUY_l
TAU_TOTAL=SQRT(TAU_TOTAL_X**2+TAU_TOTAL_Y**2)
ELSE ! recover old method, based on scalar addition
TAU_TOTAL=TAU_SIN_SCALAR+SQRT(TAU_VISC_X**2+TAU_VISC_Y**2)
ENDIF
! NaN/Inf check
IS_THIS_A_NAN = .NOT. ( TAU_TOTAL .GE. -HUGE(TAU_TOTAL) .AND. TAU_TOTAL &
.LE. HUGE(TAU_TOTAL) )
IF(IS_THIS_A_NAN)THEN
WRITE(*,*)'UFRIC = ',UFRIC
WRITE(*,*)'stopping ; tau_total = ',TAU_TOTAL
STOP
ENDIF
! Deallocate and return
DEALLOCATE(UOVERC)
DEALLOCATE(AEXP)
DEALLOCATE(S_IN_RED)
END SUBROUTINE CALC_TAU_TOTAL
!****************************************************************
!
SUBROUTINE SWIND0_NRL (SPCSIG,THETAW,ANYWND, &
UFRIC,FPM,MEMSINA,SPCDIR,KWAVE)
!
!****************************************************************
!
USE SWCOMM3
USE SWCOMM4
USE OCPCOMM4
IMPLICIT NONE
!
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmers: The SWAN team |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 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 2 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.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! For background info on SWIND0_NRL, see SWIND0
! Q: How is SWIND0_NRL different from SWIND0 ?
! A: linear wind input term is reduced in some cases
!
! See inline comments re: methods of reducing the linear wind input term
! Method 1: introduce a dependence on freq after f=1 Hz to draw it down
! Method 2: limit the overall contribution to stress to 1% of the stress
! that one would calculate from Ustar
! Note that it is OK to use both together.
! However, at high wind speeds (e.g. 38 m/s), Method 1 by itself is not
! enough, since the contribution can grow to 37% of the total stress.
! Note, Method 1 by itself actually causes a small *increase* in SWH for my
! U=12 m/s case. I'm not sure why that is. Maybe via SNL.
! Without either method, overall contribution to stress can be 1400% of
! the stress calculated from Ufric! (factor 14)
REAL SPCDIR(MDC,6)
REAL SPCSIG(MSC)
INTEGER ID ,IS
REAL FPM ,UFRIC ,THETA ,THETAW
REAL SWINEA(MDC,MSC)
REAL TEMP1 ,TEMP2
REAL CTW ,STW ,COSDIF
REAL TEMP3 ,FILTER
REAL MEMSINA(MDC,MSC,MCGRD)
REAL SIGMASIGMA1,REDUCTION
REAL KWAVE(MSC,MICMAX)
LOGICAL ANYWND(MDC)
REAL TAUX_linear,TAUY_linear,CINV2,CTH,STH,TAU_FROM_UFRIC,TAU_MAGNITUDE
INTEGER IENT
REAL ARGU
SAVE IENT
DATA IENT/0/
IF (LTRACE) CALL STRACE (IENT,'SWIND0_NRL')
!
! *** calculate linear wind input term ***
!
CTW = COS(THETAW)
STW = SIN(THETAW)
FPM = GRAV / ( 28.0 * UFRIC )
TEMP1 = PWIND(31) / ( GRAV**2 * 2. * PI )
SWINEA=HUGE(SWINEA) ! initialize for error catching
DO IS = 1, MSC
ARGU = MIN (2., FPM / SPCSIG(IS))
FILTER = EXP ( - ARGU**4 )
! note that this SIGMA is not in eq A-2 of Ris (1997). Thus, we are calculating
! "A/sigma" here, not "A"
TEMP2 = TEMP1 / SPCSIG(IS)
DO ID = 1, MDC
! Note that we include the entire range for ID, since we need to calc stress
IF ( ANYWND(ID) .AND. SPCSIG(IS) .GE. (0.7 * FPM) ) THEN
THETA = SPCDIR(ID,1)
COSDIF = SPCDIR(ID,2)*CTW + SPCDIR(ID,3)*STW
TEMP3 = ( UFRIC * MAX( 0. , COSDIF))**4
SWINEA(ID,IS) = MAX( 0. , TEMP2 * TEMP3 * FILTER )
! begin method 1 for reduction
REDUCTION=1.0
SIGMASIGMA1=SPCSIG(IS)/(2.0*PI)
! SIGMASIGMA1 is SIGMA/SIGMA1, where SIGMA1 corresponds to f=1 Hz.
IF(SIGMASIGMA1.GT.1.0)THEN
REDUCTION=SIGMASIGMA1**(-3)
ENDIF
SWINEA(ID,IS)=SWINEA(ID,IS)*REDUCTION
! end method 1 for reduction
ELSE
SWINEA(ID,IS) = 0.0
END IF
ENDDO
ENDDO
! calculate total stress
! recall: PWIND(16) 1D Rho air (1.28)
TAU_FROM_UFRIC=PWIND(16)*UFRIC**2
! calculate stress from linear wind input
TAUX_linear=0.0
TAUY_linear=0.0
DO IS = 1, MSC
CINV2=KWAVE(IS,1)/SPCSIG(IS)
DO ID = 1, MDC
! Note that we include the entire range for ID, since we need to calc stress
CTH = SPCDIR(ID,2) ! cos(theta)
STH = SPCDIR(ID,3) ! sin(theta)
TAUX_linear=TAUX_linear+CTH*CINV2*SWINEA(ID,IS)*SPCSIG(IS)*DDIR*FRINTF*SPCSIG(IS)
TAUY_linear=TAUY_linear+STH*CINV2*SWINEA(ID,IS)*SPCSIG(IS)*DDIR*FRINTF*SPCSIG(IS)
END DO
END DO
! recall: PWIND(17) Rho water (1025)
TAU_MAGNITUDE=SQRT(TAUX_linear**2+TAUY_linear**2)*PWIND(17)*GRAV
IF (TAU_MAGNITUDE.GT.(0.01*TAU_FROM_UFRIC))THEN
REDUCTION=(0.01*TAU_FROM_UFRIC)/TAU_MAGNITUDE ! method 2 for reduction
SWINEA=SWINEA*REDUCTION
!NRL IF(TESTFL)THEN
!NRL WRITE(*,*)'swind0 : tau_magnitude,tau_from_ufric,reduction 2 = ' ,&
!NRL TAU_MAGNITUDE,TAU_FROM_UFRIC,REDUCTION
!NRL ENDIF
ENDIF
DO IS = 1, MSC
IF ( SPCSIG(IS) .GE. (0.7 * FPM) ) THEN
DO ID = 1, MDC
MEMSINA(ID,IS,KCGRD(1)) = SWINEA(ID,IS)
! *** test output ***
IF (ITEST .GE. 80 .AND. TESTFL ) &
WRITE (PRTEST, 333) ID, IS, FILTER, SWINEA(ID,IS)
333 FORMAT (' ID IS FILTER WIND SOURCE ',2I4, 1X, 2(1X,E11.4))
ENDDO
ENDIF
ENDDO
!
! calculate stress (test point only)
!NRL IF(TESTFL)THEN
!NRL TAUX_linear=0.0
!NRL TAUY_linear=0.0
!NRL DO IS = 1, MSC
!NRL CINV2=KWAVE(IS,1)/SPCSIG(IS)
!NRL DO ID = 1, MDC
!NRL! Note that we include the entire range for ID, since we need to calc stress
!NRL CTH = SPCDIR(ID,2) ! cos(theta)
!NRL STH = SPCDIR(ID,3) ! sin(theta)
!NRL! SPCSIG(IS) replaces "EN" as used in SWIND_Donelan: Thus, factor AC2 is
!NRL! removed, EN=sig*ac2
!NRL TAUX_linear =TAUX_linear +CTH*CINV2*SWINEA(ID,IS)*SPCSIG(IS) &
!NRL *DDIR*FRINTF*SPCSIG(IS)
!NRL TAUY_linear =TAUY_linear +STH*CINV2*SWINEA(ID,IS)*SPCSIG(IS) &
!NRL *DDIR*FRINTF*SPCSIG(IS)
!NRL END DO
!NRL END DO
!NRL! recall: PWIND(17) 1D Rho water (1025)
!NRL TAUX_linear=TAUX_linear*PWIND(17)*GRAV
!NRL TAUY_linear=TAUY_linear*PWIND(17)*GRAV
!NRL WRITE(*,*)'SWIND0: TAUX,TAUY = ',TAUX_linear,TAUY_linear
!NRL ENDIF
!
! *** test output ***
!
IF (ITEST.GE. 60.AND.TESTFL) THEN
WRITE(PRINTF,400) KCGRD(1), THETAW*180./PI
400 FORMAT(' SWIND0: POINT THETAW :',I5,E12.4)
WRITE(PRINTF,500) TEMP1, FPM, UFRIC
500 FORMAT(' SWIND0: TEMP1 FPM UFRC :',3E12.4)
WRITE(PRINTF,*)
IF (ITEST.GE. 120.AND.TESTFL) THEN
DO IS = 1, MSC
DO ID = 1, MDC
WRITE(PRINTF,100) IS,ID,ANYWND(ID)
100 FORMAT(' IS ID ANYWND : ', 2I5,1X,L2)
ENDDO
ENDDO
ENDIF
END IF
!
RETURN
! end of subroutine SWIND0_NRL
END SUBROUTINE SWIND0_NRL
!
subroutine filsin ( memsin, idcmin, idcmax, imatra, anywnd, plwnds, isstop, genc0 )
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering and Geosciences |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmer: Marcel Zijlema |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 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 2 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.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! Authors
!
! 40.88: Marcel Zijlema
! 40.88: Erick Rogers
!
! Updates
!
! 40.88, May 2017: New subroutine
! 40.88, Jun 2017: Check for ZIEGER to allow case of negative wind input,
! wave direction opposing wind direction
!
! Purpose
!
! Fills IMATRA and GENC0 arrays with wind input term for every gridpoint per sweep direction
!
! Modules used
!
use swcomm3
use swcomm4
use ocpcomm4
!
implicit none
!
! Argument variables
!
integer, intent(in) :: isstop ! maximum frequency that is propagated within a sweep
!
integer, dimension(MSC), intent(in) :: idcmax ! maximum frequency-dependent counter in directional space
integer, dimension(MSC), intent(in) :: idcmin ! minimum frequency-dependent counter in directional space
!
real, dimension(MDC,MSC) , intent(out) :: imatra ! coefficients of right hand side of action balance equation
real, dimension(MDC,MSC,MGENR), intent(out) :: genc0 ! explicit part of generation in present vertex for output purposes
real, dimension(MDC,MSC,MCGRD), intent(in) :: memsin ! wind input stored in all active grid points
real, dimension(MDC,MSC,NPTST), intent(out) :: plwnds ! explicit part of wind input for test output
!
logical, dimension(MDC) , intent(in) :: anywnd ! determine if wind input is active for bin
!
! Local variables
!
integer :: id ! loop counter over direction bins
integer :: iddum ! counter in directional space for considered sweep
integer, save :: ient = 0 ! number of entries in this subroutine
integer :: is ! loop counter over frequency bins
!
! Structure
!
! Description of the pseudo code
!
! Source text
!
if (ltrace) call strace (ient,'filsin')
!
do is = 1, isstop
do iddum = idcmin(is), idcmax(is)
id = mod ( iddum - 1 + MDC , MDC ) + 1
if ( anywnd(id).or.ZIEGER ) then
if(TESTFL) plwnds(id,is,IPTST) = plwnds(id,is,IPTST) + &
memsin(id,is,KCGRD(1))
imatra(id,is) = imatra(id,is) + memsin(id,is,KCGRD(1))
genc0(id,is,1)= genc0(id,is,1) + memsin(id,is,KCGRD(1))
endif
enddo
enddo
!
if ( TESTFL .and. ITEST > 50 ) then
write (PRINTF,101) idcmin(1), idcmax(1), MSC, isstop
if ( ITEST > 100 ) then
do is = 1, isstop
do iddum = idcmin(is), idcmax(is)
id = mod ( iddum - 1 + MDC , MDC ) + 1
write (PRINTF,102) is, id, memsin(id,is,KCGRD(1))
enddo
enddo
endif
endif
!
101 format(' FILSIN: ID_MIN ID_MAX MSC ISTOP :',4i6)
102 format(' FILSIN: IS ID MEMSIN() :',2i6,e12.4)
!
end subroutine filsin
!
!****************************************************************
SUBROUTINE SSWELL_ROGERS (SPCSIG ,KWAVE ,IDCMIN ,IDCMAX,ISSTOP , DISSC1 ,ETOT &
,IMATDA,URMSTOP, GRAV , RHOAW , MDC,TESTFL,IPTST,PLSWEL,CGO,CDSV,FESWELL)
! note that CGo is for diagnostic purposes only
! excluded : SPCDIR AC2 DEP2 IMATRA
! WARNING: This source term may have a large impact in stationary
! computations with SWAN, e.g. on scales of 2 deg x 2 deg
!****************************************************************
IMPLICIT NONE
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 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 2 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.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
LOGICAL, INTENT(IN) :: TESTFL
INTEGER, INTENT(IN) :: ISSTOP,MDC,IPTST
INTEGER, INTENT(IN) :: IDCMIN(:) ! IDCMIN(MSC)
INTEGER, INTENT(IN) :: IDCMAX(:) ! IDCMAX(MSC)
REAL , INTENT(IN) :: ETOT
REAL , INTENT(IN) :: KWAVE(:,:) ! KWAVE(MSC,MICMAX)
REAL , INTENT(IN) :: CGo(:,:) ! CGo(MSC,MICMAX)
! CGo: group velocity without currents, but includes depth effects
REAL , INTENT(IN) :: SPCSIG(:) ! SPCSIG(MSC)
REAL , INTENT(IN) :: URMSTOP
! URMSTOP: this should be the "rms velocity amplitude" *not* the "rms velocity"
REAL , INTENT(IN) :: RHOAW
REAL , INTENT(IN) :: GRAV
REAL , INTENT(IN) :: CDSV
REAL , INTENT(IN) :: FESWELL
REAL , INTENT(INOUT) :: IMATDA(:,:) ! IMATDA(MDC,MSC)
REAL , INTENT(INOUT) :: DISSC1(:,:,:) ! DISSC1(1:MDC,1:MSC,1:MDISP)
REAL , INTENT(OUT) :: PLSWEL(:,:,:) ! PLSWEL(MDC,MSC,NPTST)
! parameters:
REAL, PARAMETER :: RE_CRIT=2.0E+5 ! Re_c or SWELLF4
! REAL, parameter :: Cdsv=1.2 ! C_dsv or SWELLF5
REAL, PARAMETER :: NU_AIR=1.51E-5 ! 1.51E-5 m2/s at 20 deg C
REAL, PARAMETER :: NU_WATER=1.0E-6 ! 1.0E-6 m2/s at 20 deg C
! local variables:
INTEGER :: IS
INTEGER :: IDDUM,ID
REAL :: AORB,USIGTOP
REAL :: RE
REAL :: SWDIS(ISSTOP) ! this is S_SWELL / E(f,theta)
REAL :: KDS_PHILLIPS
! Q: What should fe be?
! 1) 0.004 < fe < 0.013 according to JPO 2010, near final version and
! March 2009 version
! 2) "for a smooth surface, fe is of the order 0.002 to 0.008" (Fabrice
! proposes 0.0035), GRL 2009
! 3) "Most of the values of fe fall in the range 0.005 to 0.01" (2008
! Conf paper)
! 4) Supplemental information for "Ocean swell evolution from distant
! storms" : experiments with 3 constant values ( 0 , 0.0035 , 0.0070)
! 5) In Ardhuin et al. (JPO 2010), it is calculated using a formula (eq. 10)
! 6) Ardhuin et al. (JPO 2010) say: "Assuming a constant fe in (9),
! Ardhuin et al. (2009b) found that swell observations are consistent
! with 0.004 < fe < 0.013.
! Q: What should Re_crit be ?
! 1) my original code (up to and including v54):
! ==>Re_crit=1e+5
! 2) change to this code in v55 :
! ==>Re_crit=2e+5
! 3) JPO 2010, near-final version: "Semi-empirical dissipation source
! functions for ocean waves: Part I, definition, calibration and
! validation."
! ==>Re_crit=2e+5
! 4) JPO 2010, March 2009 draft
! ==>Re_crit=1e+5
! 5) GRL 2009: Observation of swell dissipation across oceans
! ==>Re_crit=1e+5
! 6) conf paper 2008: Spectral wave dissipation based on observations: a
! global validation
! ==>Re_crit=1e+5
! 7) auxiliary material for "Strong decay of steep swells observed across
! oceans"
! ==>Re_crit=1e+5
! 8) Nature submission (?), "Ocean swell evolution from distant storms"
! ==>Re_crit=2e+5
! 9) Supplemental information for "Ocean swell evolution from distant
! storms"
! ==>Re_crit=1e+5
! my notes from phone call w/Fabrice:
! "....10^5 w/bug ....2x10^5 w/out bug.....sqrt orbital displacement aorb"
!
! Method:
! @N/@t=S/sig=C(sig)*E(sig,theta)/sig=C(sig)*N(sig,theta)
! where C(sig) has units of rad/sec
! ...so if we are passing to IMATDA, we just need to pass C(sig)
! Here, SWDIS(IS) is my C(sig)
! Notes about proportionality:
! Both S_{swell,Babanin} and S_{swell,Ardhuin} are proportional
! to a^1 and omega^3.
! The only difference is in
! 1) the proportionality coefficient
! 2) the manner in which amplitude is converted to E(f)
! The dissipation by viscosity in the water is different. It is proportional to
! a^0 and omega^4.
AORB=2.0*SQRT(ETOT)
USIGTOP=URMSTOP*SQRT(2.0)
RE=4.0*USIGTOP*AORB/NU_AIR
IF(RE > RE_CRIT)THEN
DO IS=1, ISSTOP
! note that "-1" omitted since LHS
SWDIS(IS) = RHOAW * 16.0 * FESWELL * SPCSIG(IS)**2 &
* USIGTOP / GRAV
END DO
ELSE
DO IS=1, ISSTOP
SWDIS(IS) = RHOAW * CDSV * 2.0 * KWAVE(IS,1) &
* SQRT(2.0 * NU_AIR * SPCSIG(IS))
END DO
ENDIF
! Add dissipation by viscosity in the water.
! Method: Phillips pg 37 "the energy density of the wave field decreases
! as exp(-4*nu*k^2*t)"
! Thus, Kds_Phillips=4*nuw*(k^2), where Kds is (@E/@t)/E
! For other ideas about S_visc, see Pierson et al. (2013+),
! "Rain-induced attenuation of water waves"
! Alternate reference: Dulov and Kosnnik (Izvestiya, Atmospheric and
! Oceanic Physics, 2009, vol 45, No 3, pp 380-391) :
! gamma=-4*nu*k^2, where nu=1.3e-6 m2/s
! They also say, "This theory suggests that local three-wave
! interactions are a dominant mechanism for forming the spectrum in the
! entire capillary range except for its shortest wavelength portion,
! where viscous dissipation becomes important"
! The implication is that viscosity is only important for the shortest
! waves in the capillary range, for "forming the spectrum", and perhaps
! for dissipation also.
! Problem identified: if we use this for calculation of tau_wave_to_atm,
! then we should not include Sds,visc...that part would go to
! tau_wave_to_ocean
! DO IS=1, ISSTOP
! KDS_PHILLIPS=4.0*NU_WATER*(KWAVE(IS,1)**2) ! UNITS OF RADIAN^2/S
! SWDIS(IS) = SWDIS(IS) + KDS_PHILLIPS
! END DO
DO IS=1, ISSTOP
! Only fill the values for the current sweep
DO IDDUM = IDCMIN(IS), IDCMAX(IS)
ID = MOD(IDDUM - 1 + MDC, MDC) + 1
IMATDA(ID,IS) = IMATDA(ID,IS) + SWDIS(IS)
DISSC1(ID,IS,4) = DISSC1(ID,IS,4) + SWDIS(IS)
IF (TESTFL) PLSWEL(ID,IS,IPTST) = -1.*(SWDIS(IS))
END DO
END DO
RETURN
END SUBROUTINE SSWELL_ROGERS
SUBROUTINE SSWELL_ARDHUIN (SPCSIG, THETAW, KWAVE, IDCMIN, IDCMAX, ISSTOP, DISSC1, ETOT, IMATDA &
,SPCDIR, UFRIC, URMSTOP, GRAV , RHOAW , TESTFL,IPTST,PLSWEL,MDC,CGO,CDSV)
!---------------------------------------------------------------------------------------------
! SUBROUTINE SSWELL_ARDHUIN calculates swell dissipation using equations (8) and (9) in
! Ardhuin et al (2010). These equations are the same as Rogers et al 2012 in Subroutine
! SSWELL_ROGERS, except fe is calculated using wind-wave angle, friction velocity and
! the significant surface orbital velocity.
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 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 2 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.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
! Yalin Fan
!
! 1. Updates
!
! 2. PURPOSE
! Compute swell dissipation / attenuation for a third generation wave model
!
! 3. METHOD
!
! Re > Recr
! Sswell(f,theta) = 16 * (Rho_a/Rho_w) * fe * sigma^2 *uorb * E(f,theta) / g (9)
! fe = s1 * {feGM + [|s3| + s2 * cos(theta - theta_wnd)] * u* /uorb }
!
! Ardhuin calculates using method by Grant and Madsen, which is quite involved,
! so a shortcut taken here, based on what Ardhuin says: "fe,GM [is] of the order
! of 0.003" (pg 1921, right column)
!
! Re <= Recr
! Sswell(f,theta) = 2 * (Rho_a/Rho_w) * Cdsv * k * sqrt(2*nv*sigma) * F(f,theta) (8)
!
! 4. Argument variables
!
! SPCSIG: Relative frequencies in computational domain in sigma-space
! SPCDIR: (*,1); spectral directions (radians) 30.82
! (*,2); cosine of spectral directions 30.82
! (*,3); sine of spectral directions 30.82
! (*,4); cosine^2 of spectral directions 30.82
! (*,5); cosine*sine of spectral directions 30.82
! (*,6); sine^2 of spectral directions 30.82
! IS Counter of relative frequency band
! ID Counter of directional distribution
! MSC Maximum counter of relative frequency
! MDC Maximum counter of directional distribution
!
! one and more dimensional arrays:
! ---------------------------------
! KWAVE 2D Wavenumber
! IMATDA 2D Coefficients of main diagonal of matrix
! IDCMIN 1D Frequency dependent counter
! IDCMAX 1D Frequency dependent counter
!
! 5. SUBROUTINES CALLING
! SOURCE
!-------------------------------------------------------------------------------------
IMPLICIT NONE
!
LOGICAL, INTENT(IN) :: TESTFL
INTEGER, INTENT(IN) :: ISSTOP,MDC,IPTST
INTEGER, INTENT(IN) :: IDCMIN(:) ! IDCMIN(MSC)
INTEGER, INTENT(IN) :: IDCMAX(:) ! IDCMAX(MSC)
REAL , INTENT(IN) :: ETOT
REAL , INTENT(IN) :: KWAVE(:,:) ! KWAVE(MSC,MICMAX)
REAL , INTENT(IN) :: CGo(:,:) ! CGo(MSC,MICMAX)
! CGo: group velocity without currents, but includes depth effects
REAL , INTENT(IN) :: SPCSIG(:) ! SPCSIG(MSC)
REAL , INTENT(IN) :: SPCDIR(:,:) ! SPCDIR(MDC,6)
REAL , INTENT(IN) :: URMSTOP
! URMSTOP: this should be the "rms velocity amplitude" *not* the "rms velocity"
REAL , INTENT(IN) :: RHOAW
REAL , INTENT(IN) :: GRAV
REAL , INTENT(IN) :: CDSV
REAL , INTENT(IN) :: THETAW
REAL , INTENT(IN) :: UFRIC
REAL , INTENT(INOUT) :: IMATDA(:,:) ! IMATDA(MDC,MSC)
REAL , INTENT(INOUT) :: DISSC1(:,:,:) ! DISSC1(1:MDC,1:MSC,1:MDISP)
REAL , INTENT(OUT) :: PLSWEL(:,:,:) ! PLSWEL(MDC,MSC,NPTST)
! parameters:
REAL, PARAMETER :: RE_CRIT=2.0E+5 ! Re_c or SWELLF4
! REAL, parameter :: Cdsv=1.2 ! C_dsv or SWELLF5
REAL, PARAMETER :: NU_AIR=1.51E-5 ! 1.51E-5 m2/s at 20 deg C
REAL, PARAMETER :: NU_WATER=1.0E-6 ! 1.0E-6 m2/s at 20 deg C
! local variables:
INTEGER :: IS
INTEGER :: IDDUM,ID
REAL :: AORB,USIGTOP
REAL :: RE
REAL :: SWDIS(MDC,ISSTOP) ! this is S_SWELL / E(f,theta)
REAL :: KDS_PHILLIPS
REAL :: FESWELL, CTW, STW, COSDIF
REAL, PARAMETER :: s1 = 0.8 ! Ardhuin, table A1
REAL, PARAMETER :: s2 = -0.018 ! fixed by Ardhuin (pg 1921, right column)
REAL, PARAMETER :: s3 = 0.015 ! fixed by Ardhuin (pg 1921, right column)
REAL, PARAMETER :: feGM = 0.003 ! fe Grant and Madsen, see notes above
CTW = COS(THETAW)
STW = SIN(THETAW)
AORB=2.0*SQRT(ETOT)
USIGTOP=URMSTOP*SQRT(2.0)
RE=4.0*USIGTOP*AORB/NU_AIR
IF(RE > RE_CRIT)THEN
DO IS=1, ISSTOP
DO IDDUM = IDCMIN(IS), IDCMAX(IS)
ID = MOD(IDDUM - 1 + MDC, MDC) + 1
COSDIF = SPCDIR(ID,2)*CTW + SPCDIR(ID,3)*STW
FESWELL = s1*(feGM + (abs(s3) + s2*COSDIF)*UFRIC/USIGTOP)
! note that "-1" omitted since LHS
SWDIS(ID,IS) = RHOAW * 16.0 * FESWELL * SPCSIG(IS)**2 &
* USIGTOP / GRAV
END DO
END DO
ELSE
DO IS=1, ISSTOP
DO IDDUM = IDCMIN(IS), IDCMAX(IS)
ID = MOD(IDDUM - 1 + MDC, MDC) + 1
SWDIS(ID,IS) = RHOAW * CDSV * 2.0 * KWAVE(IS,1) &
* SQRT(2.0 * NU_AIR * SPCSIG(IS))
END DO
END DO
ENDIF
DO IS=1, ISSTOP
! Only fill the values for the current sweep
DO IDDUM = IDCMIN(IS), IDCMAX(IS)
ID = MOD(IDDUM - 1 + MDC, MDC) + 1
IMATDA(ID,IS) = IMATDA(ID,IS) + SWDIS(ID,IS)
DISSC1(ID,IS,4) = DISSC1(ID,IS,4) + SWDIS(ID,IS)
IF (TESTFL) PLSWEL(ID,IS,IPTST) = -1.*(SWDIS(ID,IS))
END DO
END DO
RETURN
END SUBROUTINE SSWELL_ARDHUIN
SUBROUTINE SSWELL_ZIEGER (SPCSIG, KWAVE, AC2, CG, ISSTOP ,IDCMIN ,IDCMAX, MDC, DISSC1, IMATDA, TESTFL, IPTST, PLSWEL)
!---------------------------------------------------------------------------------------------
! SUBROUTINE SSWELL_ZIEGER calculates swell dissipation using equations (21) - (23) in
! Zieger et al (2013), which are originally from Bowden (1950) and Babanin (2011)
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 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 2 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.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Authors
! Yalin Fan
!
! 1. Updates
!
! 2. PURPOSE
! Compute swell dissipation / attenuation for a third generation wave model
!
! 3. METHOD
!
! Sswell(sigma,d) = -2.0*b1*sigma*sqrt(Bn)*E(sigma,d)
!
! 4. Argument variables
!
! 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
!
! one and more dimensional arrays:
! ---------------------------------
! KWAVE 2D Wavenumber
! IMATDA 2D Coefficients of main diagonal of matrix
! IDCMIN 1D Frequency dependent counter
! IDCMAX 1D Frequency dependent counter
!
! 5. SUBROUTINES CALLING
! SOURCE
!-------------------------------------------------------------------------------------
USE SWCOMM3, ONLY : KCGRD, DDIR, MSC, B1Z
IMPLICIT NONE
!
! Subroutine arguments:
LOGICAL, INTENT(IN) :: TESTFL
INTEGER, INTENT(IN) :: ISSTOP,MDC,IPTST
INTEGER, INTENT(IN) :: IDCMIN(:) ! IDCMIN(MSC)
INTEGER, INTENT(IN) :: IDCMAX(:) ! IDCMAX(MSC)
REAL , INTENT(IN) :: SPCSIG(:) ! SPCSIG(MSC)
REAL , INTENT(IN) :: CG(:,:) ! CG(MSC,ICMAX)
REAL , INTENT(IN) :: KWAVE(:,:) ! KWAVE(MSC,MICMAX)
REAL , INTENT(IN) :: AC2(:,:,:) ! AC2(MDC,MSC,MCGRD)
REAL , INTENT(INOUT) :: IMATDA(:,:) ! IMATDA(MDC,MSC)
REAL , INTENT(INOUT) :: DISSC1(:,:,:) ! DISSC1(1:MDC,1:MSC,1:MDISP)
REAL , INTENT(OUT) :: PLSWEL(:,:,:) ! PLSWEL(MDC,MSC,NPTST)
! Local variables:
INTEGER :: IDDUM ,ID ,IS
REAL :: DMAX,AINV,BN
! 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 :: SWDIS(ISSTOP) ! this is S_SWELL / E(f,theta)
! 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)
! SQRTBN(IS)=ANAR(IS)*(SIGDENS(IS)*(KWAVE(IS,1)**3)*CG(IS,1))
END DO
! Calculate Swell dissipation Sswl divided by E(sigma,theta) (F(sigma,theta)
! This is from WW3v4 and Zieger et al. (2015) eq. 23.
! This does *not* include the steepness-dependent B1 that is introduced in
! WW3v5 and Zieger et al. (2015) eq. 28.
DO IS=1, ISSTOP
SWDIS(IS) = 2.0*B1Z*SPCSIG(IS)*SQRTBN(IS)/3.0
END DO
! Re: usage of DISCC1(_,_,4): In case of Zieger Sswell,
! this should go to tau_wave_to_ocean, not tau_wave_to_atm
DO IS=1, ISSTOP
! Only fill the values for the current sweep
DO IDDUM = IDCMIN(IS), IDCMAX(IS)
ID = MOD(IDDUM - 1 + MDC, MDC) + 1
IMATDA(ID,IS) = IMATDA(ID,IS) + SWDIS(IS)
DISSC1(ID,IS,4) = DISSC1(ID,IS,4) + SWDIS(IS)
IF (TESTFL) PLSWEL(ID,IS,IPTST) = -1.*(SWDIS(IS))
END DO
END DO
RETURN
END SUBROUTINE SSWELL_ZIEGER
!----------------------------------------------------------------------------
SUBROUTINE integrate(ansNum, x,y,np)
!-------------------------------------------------------------------------
IMPLICIT NONE
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 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 2 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.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! Use Trapezoidal Rule to Integrate the area under the curve
! specified by the data points in x and y
!
! John Mahaffy 4/9/95
!
!
! Local Variables
!
! esterr - estimated error for Trapizoidal integration
! sum2 - Trapizoidal integration using every other available point
!
REAL , INTENT(IN) :: X(:)
REAL , INTENT(IN) :: Y(:)
REAL , INTENT(OUT) :: ANSNUM
INTEGER , INTENT(IN) :: NP
INTEGER :: I
ANSNUM=0.0
! np=size(x) ! we cannot do this because array may be partially filled
DO I=1,NP-1
ANSNUM=ANSNUM+.5*(Y(I)+Y(I+1))*(X(I+1)-X(I))
END DO
RETURN
END SUBROUTINE INTEGRATE
!-------------------------------------------------------------------------
SUBROUTINE SURF_ROUGH_FAN (FPI,U,UST,CD)
!
! Parameter list
! ----------------------------------------------------------------
! FPI Real I Peak-input frequency.
! U Real I True Wind speed.
! UST Real O Friction velocity.
! CD Real O Drag coefficient.
! ----------------------------------------------------------------
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
REAL, INTENT(IN) :: FPI, U
REAL, INTENT(OUT) :: UST, CD
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
INTEGER :: I1, ITT, I
REAL :: SQRTH, SIX, R1, WNP, CP, UNZ, ALPHA, &
RDCH, AFP, Z0, &
G,CAPA,PII,CW,Z0FTN,CDFTN,CPEAK, &
ZCH,X,WAGE,A1,B1,A2,B2,FX,xx1,slope
!/
!/ ------------------------------------------------------------------- /
!/
!
! 0. Definition of constant
G=9.8
CAPA=0.4
PII=3.141517
!
!---- A1=0.57 produce too high Hs ----------
!---- First try, use 0.4521, still too large
!---- Second try, use 0.3931, Hs is still 1.5m higher
!---- Third try, use 0.354
A1=0.354
B1=-1.5
B2=1.5
!
!---- a transition point of wave age
CW=11
!
! 1. Defining CD function ---------------------------------------------- *
IF(U.LE.0.) THEN
Z0FTN=0.000001
CDFTN=0.000001
ELSE IF(U.LE.12.5) THEN
Z0FTN=0.0185/G*(0.001*U*U+0.028*U)**2
CDFTN=0.4*0.4*(LOG(10./Z0FTN)**(-2))
ELSE
Z0FTN=(0.085*U-0.58)*0.001
CDFTN=0.4*0.4*(LOG(10./Z0FTN)**(-2))
END IF
!
!---- UST
IF(UST.LE.0.) UST=SQRT(CDFTN)*U
!
!---- Phase velocity at peak
CPEAK=SQRT(G**2/(2.*PII*FPI)**2)
!
!---- First check----------------------------------
IF(U.LE.0.) THEN
Z0 = 0.000001
ZCH = 0.000001
CD = 0.000001
UST = 0.000001
GO TO 900
ELSE
!
!----- Initial guess of UST using a function
X=SQRT(CDFTN)*U
!
!----- Iteration to get UST------------------------------------------- *
!
DO I=1,5
!
!----- Wave age
WAGE=CPEAK/X
!
!----- Relationship between wave age and Charnock Coeff.
xx1=100.0
slope= 0.15*(U/12.5)
ZCH=0.023/xx1**slope*WAGE**slope
!
!----- Calculation of Charnock Coefficient(Zch-Wage relationship)
!
FX=X/CAPA*LOG(10.*G/(X**2*ZCH))-U
!
IF(ABS(FX).GE. 0.001) THEN
X=X-(X/CAPA*LOG(10*G/(X**2*ZCH))-U)/ &
((1/CAPA)*(LOG(10*G/(X**2*ZCH))-2.))
ELSE
UST=X
END IF
END DO
!------------------ End of Iteration --------------------------------- *
!
UST=X
Z0=UST**2*ZCH/G
IF (Z0.LT.0.00001) Z0 = 0.00001
CD=(1/CAPA*LOG(10/Z0))**(-2)
!
!----- Second Check--------
! IF((WAGE.GT.3.).AND.(WAGE.LT.32.)) THEN
!
!------- Roughness length and Drag coefficient
! Z0=UST**2*ZCH/G
! CD=(1/CAPA*LOG(10/Z0))**(-2)
! ELSE
! Z0=Z0FTN
! CD=CDFTN
! UST=SQRT(CD)*U
! ZCH=Z0*G/UST**2
! END IF
!---------------------------
!
END IF
!--------------------------- First Check Ending
!
!----- Final check Cd
! IF((ABS(CDFTN-CD).GT.0.001).OR.(CD.GT.0.005)) THEN
! Z0=Z0FTN
! CD=CDFTN
! UST=SQRT(CD)*U
! UST=SQRT(CD)*U/ASF(ISEA)
! ZCH=Z0*G/UST**2
! END IF
!
!---- Direction of Stress
!
900 continue
RETURN
END SUBROUTINE SURF_ROUGH_FAN
SUBROUTINE SURF_ROUGH_ECMWF(U, UST, S_IN, SPCSIG, KWAVE, CD)
!------------------------------------------------------------------------------
! variable list
! ----------------------------------------------------------------
! S_IN Real I Wind input source function
! SPCSIG Real I Relative frequencies in computational domain in sigma-space
! KWAVE Real I 2D Wavenumber
! U Real I True Wind speed.
! UST Real O Friction velocity.
! CD Real O Drag coefficient.
! ----------------------------------------------------------------
!/
USE SWCOMM3
USE SWCOMM4
USE OCPCOMM4
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
REAL, INTENT(IN) :: S_IN(MDC,MSC,MGENR)
REAL, INTENT(IN) :: SPCSIG(MSC)
REAL, INTENT(IN) :: KWAVE(MSC,ICMAX)
REAL, INTENT(IN) :: U
REAL, INTENT(INOUT) :: UST
REAL, INTENT(OUT) :: CD
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
INTEGER :: I1, ITT, I, IS, IDIR, NDIR, NF_OLD, NF_NEW, K, II
REAL :: G, CAPA, PII, afa, TAU_wav, RHOW, RHOA, TAU_total
REAL :: FREQ_TMP, FRAT, Z0FTN, CDFTN, z0, ze, a, b, X, CDRAG
REAL, ALLOCATABLE :: S_IN_L(:,:,:)
REAL, ALLOCATABLE :: CINV_L(:), SIGMA_L(:), DF_L(:)
! REAL :: SQRTH, SIX, R1, WNP, CP, UNZ, ALPHA, &
! RDCH, AFP, Z0, afa, &
! G,CAPA,PII,CW,Z0FTN,CDFTN,CPEAK, &
! ZCH,X,WAGE,A1,B1,A2,B2,FX,xx1,slope
!/
!/ ------------------------------------------------------------------- /
!/
!
! 0. Definition of constant
G=9.81
CAPA=0.4
PII=3.141517
AFA = 0.01
RHOW = 1025.0
RHOA = 1.28
!
! 1. Difining CD function ---------------------------------------------- *
IF(U.LE.0.) THEN
Z0FTN=0.000001
CDFTN=0.000001
ELSE IF(U.LE.12.5) THEN
Z0FTN=0.0185/G*(0.001*U*U+0.028*U)**2
CDFTN=0.4*0.4*(LOG(10./Z0FTN)**(-2))
ELSE
Z0FTN=(0.085*U-0.58)*0.001
CDFTN=0.4*0.4*(LOG(10./Z0FTN)**(-2))
END IF
!
!---- UST
IF(UST.LE.1E-10) THEN
CDRAG = CDFTN
UST = SQRT ( CDRAG ) * U
GOTO 110
ENDIF
!...S_in(theta,freq)
NDIR=SIZE(S_IN,1)
NF_OLD=SIZE(S_IN,2)
!...find nf_new
FRAT=FRINTF+1.0 ! =freq(nf_old)/freq(nf_old-1)
FREQ_TMP=SPCSIG(NF_OLD)/(2.0*PII)
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
ALLOCATE(CINV_L(NF_NEW))
ALLOCATE(SIGMA_L(NF_NEW))
ALLOCATE(DF_L(NF_NEW))
ALLOCATE(S_IN_L(NDIR,NF_NEW,MGENR))
CINV_L = 0.0
SIGMA_L = 0.0
DF_L = 0.0
S_IN_L = 0.0
SIGMA_L(1:NF_OLD)=SPCSIG(1:NF_OLD)
DO IS=1,NF_OLD
CINV_L(IS)=KWAVE(IS,1)/SPCSIG(IS)
ENDDO
DO IS = 1, NF_OLD
DO IDIR = 1, NDIR
DO II = 1, MGENR
S_IN_L(IDIR,IS,II)=S_IN(IDIR,IS,II)
ENDDO
ENDDO
ENDDO
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
!.... ...since GENC0 (S_IN here) = Sin / f, we use therefore f^-3 as tail
!... Also, because LFACTOR is applied to the tail part of S_IN before calculate
!... tau_sin_x and tau_sin_y, and the reduction is significant (on the order of 10^-4)
!... we should apply the same reduction to GENC0 as well. This requires to bring array
!... LFACTOR to this subroutine.
! *** Note that LFACTOR_SRC is renamed here as RDFSIN ***
!Yalin: S_IN(NF_OLD) / GENC0(NF_OLD) was already drawn down by LFACTOR(NF_OLD)
! before pass into routine ADDDIS. But, in CALC_LFACTOR, the tail is added using
! the original S_IN(NF_OLD) before add tail then calculate LFACTOR.
! Thus, we need to restore S_IN(NF_OLD) to its original value before add
! tail here
DO IDIR=1,NDIR
DO II = 1, MGENR
S_IN_L(IDIR,IS,II) = (S_IN_L(IDIR,NF_OLD,II)/RDFSIN(NF_OLD)) &
* (SIGMA_L(NF_OLD)/SIGMA_L(IS))**3 * RDFSIN(IS)
ENDDO
ENDDO
ENDDO
TAU_wav=0.0
DO IDIR=1,NDIR
DO IS=1,NF_NEW
DF_L(IS)=SIGMA_L(IS)*FRINTF
TAU_wav=TAU_wav+RHOW*G*S_IN_L(IDIR,IS,1)*CINV_L(IS)*SIGMA_L(IS) &
*DF_L(IS)*DDIR
END DO
ENDDO
! now we got tau_w, need to cal z0 use tau_w, then iteration
Z0 = AFA * UST**2 / G
TAU_total = RHOA * UST**2
ZE = Z0/SQRT(1-TAU_wav/TAU_total)
X = SQRT(G*Z0/AFA)*LOG((10.0+ZE-Z0)/Z0)/CAPA - U
IF (ABS(X) .LT. 1E-6) GO TO 100
IF (X.GT.0.0) THEN
A = 0
B = Z0
ELSE
A = Z0
B = 0.1
ENDIF
K = 0
80 Z0 = (A+B)/2
UST = SQRT(Z0*G/AFA)
TAU_total = RHOA * UST**2
ZE = Z0/SQRT(1-TAU_wav/TAU_total)
X = SQRT(G*Z0/AFA)*LOG((10.0+ZE-Z0)/Z0)/CAPA - U
IF (X .GT. 0.0) THEN
B = (A+B)/2
ELSE
A = (A+B)/2
ENDIF
K = K + 1
IF (ABS(X) > u*1E-6 .or. X < 0) GOTO 80
100 CONTINUE
UST = SQRT(Z0*G/AFA)
DEALLOCATE(CINV_L)
DEALLOCATE(SIGMA_L)
DEALLOCATE(DF_L)
DEALLOCATE(S_IN_L)
110 CONTINUE
RETURN
END SUBROUTINE SURF_ROUGH_ECMWF
!------------------------------------------------------------------------------
REAL FUNCTION KUIK(E_AT_FREQ,DTHETA_RAD,CTH,STH,DIAGNOSTICS)
!------------------------------------------------------------------------------
!
IMPLICIT NONE
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2017 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 2 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.
!
! A copy of the GNU General Public License is available at
! http://www.gnu.org/copyleft/gpl.html#SEC3
! or by writing to the Free Software Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
!
! 0. Update history
!
! Jan 3 2013 Initial version by Erick Rogers
!
! 1. Purpose:
!
! Compute circular rms directional spreading (Kuik et al.)
!
! 2. Method
!
! Kiuk et al. (JPO 1988), see also Rogers and Wang (2007) and Mardia book
! referenced therein.
!
! 3. Parameter list:
!
! Type I/O Name Description
!--------------------------------------------------------------------
REAL, INTENT(IN) :: E_AT_FREQ(:) ! E(f,theta) for one frequency f
REAL, INTENT(IN) :: DTHETA_RAD ! delta_theta
REAL, INTENT(IN) :: CTH(:) ! cos(theta)
REAL, INTENT(IN) :: STH(:) ! sin(theta)
LOGICAL, INTENT(IN) :: DIAGNOSTICS ! Flag for diagnostic output
!
! 4. Error messages
!
! 5. Called by:
!
! 6. Subroutines used:
!
! 7. Remarks
!
! 8. Structure
!
! 9. Switches
!
! 10. Source code:
!------------------------------------------------------------------------------
! local variables
!
INTEGER NDIR,IDIR ! number of directions and counter for directions
REAL EF ! E(f) at this frequency
REAL INTEG3,INTEG4,a1,b1,m1,DSPR ! temporary variables for integration
!---------------------------------------------------------------------
NDIR=SIZE(E_AT_FREQ)
INTEG3=0.0
INTEG4=0.0
EF=0.0
DO IDIR=1,NDIR
INTEG3=INTEG3+E_AT_FREQ(IDIR)*CTH(IDIR)*DTHETA_RAD
INTEG4=INTEG4+E_AT_FREQ(IDIR)*STH(IDIR)*DTHETA_RAD
EF=EF+E_AT_FREQ(IDIR)*DTHETA_RAD
END DO
! note that if dtheta_rad is constant, it cancels out, so it can be omitted from this routine
A1=INTEG3/EF
B1=INTEG4/EF
M1=SQRT(A1**2+B1**2)
DSPR=SQRT(2*(1-M1))
!NRL IF(DIAGNOSTICS)THEN
!NRL WRITE(415,*)'% calc_Kuik: using dtheta_rad = ',dtheta_rad,' radians'
!NRL WRITE(415,*)'% calc_Kuik: DSPR = ',DSPR,' radians'
!NRL END IF
KUIK=DSPR
RETURN
END FUNCTION KUIK
END MODULE SDSBABANIN