#include "wwm_functions.h"
MODULE SdsBabanin
CONTAINS
subroutine calc_Sds(ip,nfreq,EDENS,f,Kds,ANAR_IN,LPOW,MPOW,a1,a2,ACLOC,IMATRA,IMATDA,SSDS)
USE DATAPOOL
IMPLICIT NONE
integer, intent(in) :: ip
real(rkind) , INTENT(IN) :: EDENS(:) ! E(f) m2/Hz
real(rkind) , INTENT(IN) :: f(:) ! f in Hz
real(rkind) , INTENT(IN) :: ANAR_IN(:) ! directional narrowness as defined in Babanin publications (input value)
real(rkind) , INTENT(IN) :: a1 ! coefficient on T1
real(rkind) , INTENT(IN) :: a2 ! coefficient on T2
! power on threshold exceedence in T1 (precise definition depends
! on whether using "U" vs "D" method)
real(rkind) , INTENT(IN) :: LPOW
real(rkind) , INTENT(IN) :: MPOW
! power on threshold exceedence in T2 (precise definition depends
! on whether using "U" vs "D" method)
integer , INTENT(IN) :: nfreq ! # freqs
real(rkind) , INTENT(OUT) :: SSDS(MSC,MDC) !
real(rkind) , INTENT(IN) :: ACLOC(MSC,MDC)
real(rkind) , INTENT(INOUT) :: IMATRA(MSC,MDC), IMATDA(MSC,MDC)
real(rkind) , INTENT(OUT) :: Kds(:) ! Kds(f)=Sds(f)/E(f)
real(rkind) :: Sds(nfreq) ! Sds(f), the source term
real(rkind) :: ndedens(nfreq) ! NDEDENS(f)=DEDENS(f)/EDENST(f)
real(rkind) :: dedens(nfreq) ! DEDENS(f)=EDENS(f)-EDENST(f)
real(rkind) :: edenst(nfreq) ! EDENST(f)=(2*g.^2)*((2*pi).^-4)*(f.^-5).*(A.^-1)*Bnt
real(rkind) :: T1(nfreq) ! inherent dissipation (divided by E(f))
real(rkind) :: T2(nfreq) ! induced dissipation (divided by E(f))
real(rkind) :: ST1(nfreq) ! inherent dissipation
real(rkind) :: ST2(nfreq) ! induced dissipation
real(rkind) :: ANAR(nfreq) ! directional narrowness as defined in Babanin publications (value used)
real(rkind) :: xff(nfreq),ADF(nfreq) ! temporary arrays for integration
REAL(rkind) :: ASUM ! temporary variable for integration
REAL(rkind) :: BNT ! is an empirical coefficient related to the spectral density in the saturation range (Banner et al 2007)
real(rkind) :: ST1_INT,ST2_INT,Sds_int ! integrated values (for test output only)
INTEGER :: II,IS,ID ! counters
! real :: ELIM ! needed for 42D (can comment but don't delete)
Bnt=(0.035**2) ! Use the Bnt given by Banner et al 2007
! ----------------------------------------------------------------------------------------------------------------
! #1 a1 and a2 get a1 and a2 as a function of U/cp
! ------------------------------------------------------------------------------------------------------------------
! find the fp and determine the cp
! imax=maxloc(EDENS,1)
! fpcheck=f(imax)
! imax=-999
! fp=fpinput
! needed: A(f), see Young and Babanin eq 19.
! Normally, it is read in and used, but per recommendation by Alex that it 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(.false.)then
write(DBG%FHNDL,*)'a1,a2 = ',a1,a2
endif
! two matrix operations follow:
EDENST= (2.*g9**2) * ((2.0*pi)**(-4)) * (f**(-5)) * (ANAR**(-1)) * &
& Bnt ! see eq 9 of "Implementation of new..." Bnt is sigma_thr
DEDENS=EDENS-EDENST
! if DEDENS < 0 set to zero
DEDENS=MAX(0.0_rkind,DEDENS)
! DO IS = 1, nfreq
! WRITE(DBG%FHNDL,*) IS, EDENS(IS), EDENST(IS), DEDENS(IS)
! END DO
! ELIM=maxval(Edens)*1.0e-5 ! needed for v42D (can comment but don't delete)
do is=1,nfreq
!v42U normalize by EDENST, a la Fabrice, this means that, with L=2
! and large exceedence, the dissipation will get very strong
!v42U NDEDENS(is)=DEDENS(is)/EDENST(is)
!v42D: normalize by EDENS, a la Kakha, 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
! if(Edens(is).gt.ELIM)then
! NDEDENS(is)=DEDENS(is)/EDENS(is) ! v42D
! else
! NDEDENS(is)=0.0
! end if
NDEDENS(is)=DEDENS(is)/EDENST(is) ! v42U
enddo
NDEDENS=MAX(0.0_rkind,NDEDENS)
! -------------------------------------------------------------------------------
! calculate T1
! -------------------------------------------------------------------------------
! T1=a1*f*A*DEDENS ! original matrix operation
do is=1,nfreq
T1(is)=a1*f(is)*ANAR(is)*NDEDENS(is)**LPOW
end do
! ------------------------------------------------------------------------------------
! calculate T2 an integration from fp to f
! ------------------------------------------------------------------------------------
! OLD version: uses fp that falls on a node
! New version: provided peak does not necessarily fall on a node
! New New version: don't worry about fp, since stuff below fp is not breaking, so it isn't part of the calculation anyway.
xff=0.
do is=1,nfreq
ASUM=0.
do ii=1,is
xff(ii)=f(ii)
ADF(ii)=ANAR(ii)*NDEDENS(ii)**MPOW
enddo
CALL integrate(ASUM,xff,ADF,is)
T2(is)=a2*ASUM
enddo
T1=max(0._rkind,T1) ! WER a matrix operation
T2=max(0._rkind,T2) ! WER a matrix operation
Kds=T1+T2 ! WER a matrix operation
! (point output write location 2)
if(.false.)then
do is=1,nfreq
write(DBG%FHNDL,205)f(is),EDENS(is),EDENST(is),T1(is),T2(is)
end do
endif
205 format('threshold: ',5(1x,D12.5))
! figure out integrated T1 and T2 (for output purposes only)
! 3/9/2009 Since we have already non-dimensionalized by using NDEDENS instead of DEDENS, we do it backwards: Sds=Kds*Edens
do is=1,nfreq
Sds(is)=Kds(is)*Edens(is) ! note that in most cases, the calling routine will not use Sds
ST1(is)=T1(is)*Edens(is)
ST2(is)=T2(is)*Edens(is)
end do
DO IS = 1, MSC
DO ID = 1, MDC
SSDS(IS,ID) = KDS(IS)
IF (ICOMP .GE. 2) THEN
IMATDA(IS,ID) = IMATDA(IS,ID) + Kds(IS)
ELSE
IMATDA(IS,ID) = IMATDA(IS,ID) - kds(IS)
IMATRA(IS,ID) = IMATRA(IS,ID) - kds(IS) * ACLOC(IS,ID)
END IF
END DO
END DO
CALL integrate(ST1_INT,f,ST1,nfreq)
CALL integrate(ST2_INT,f,ST2,nfreq)
CALL integrate(Sds_int,f,Sds,nfreq)
! (point output write location 3)
if(.false.)then
WRITE(DBG%FHNDL,*)'integral of T1,T2,Sds = ',ST1_INT,ST2_INT,Sds_INT
endif
END subroutine calc_Sds
! ************************************************************************************************
!
! ************************************************************************************************
SUBROUTINE integrate(ansNum, x,y,np)
! ====================================================
USE DATAPOOL, ONLY : RHOA, RHOW, RKIND
IMPLICIT NONE
!
! 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(rkind) , INTENT(IN) :: X(:)
REAL(rkind) , INTENT(IN) :: Y(:)
integer i
! real sum2, esterr
REAL(rkind) , INTENT(out) :: ANSNUM
INTEGER , INTENT(IN) :: NP
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
!
! The following only works if np is an odd integer
!
! sum2=0.0
! do i=1,np-1,2
! sum2=sum2+.5*(y(i)+y(i+2))*(x(i+2)-x(i))
! end do
return
end SUBROUTINE integrate
subroutine calc_Lfactor(ip,Lfactor_S,S_in,DDIR_RAD,SIGMA_S,FRINTF,CINV_S,grav,WIND10,TESTFL)
use datapool, only : rhoa, rhow, pi, ufric, rkind, DBG
use datapool, only : TAUTOT, CD, Z0, ALPHA_CH, G9, ac2
implicit none
! 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. (Monoticity of S_in, or lack thereof, is irrelevant.)
! points made here:
! 1) we have no idea what the tangential stress is, but we no that
! it isn't less
! that zero, so we say that the normal stress cannot be
! greater than the total stress.
! 2) we don't know what the contribution is beyone fmax, but we know
! that it isn't negative, so again,
! we 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
! 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? (the original thinking
! was to use an fp value in a manner similar to
! Tolman and Chalikov, which is based on U/C somehow.
! 5) S_in=S_in_linear+S_in_exponential....same argument as (1)
Integer, Intent(In) :: ip
REAL(rkind) , INTENT(OUT) :: Lfactor_S(:)
REAL(rkind) , INTENT(IN) :: S_in(:,:)
REAL(rkind) , INTENT(IN) :: DDIR_RAD
REAL(rkind) , INTENT(IN) :: FRINTF
REAL(rkind) , INTENT(IN) :: CINV_S(:)
REAL(rkind) , INTENT(IN) :: SIGMA_S(:)
REAL(rkind) , INTENT(IN) :: grav
REAL(rkind) , INTENT(IN) :: WIND10
LOGICAL , INTENT(INOUT) :: TESTFL
REAL(rkind) , ALLOCATABLE :: S_in1D_S(:)
REAL(rkind) , ALLOCATABLE :: S_in1D_L(:)
REAL(rkind) , ALLOCATABLE :: DF(:)
REAL(rkind) , ALLOCATABLE :: CINV_L(:)
REAL(rkind) , ALLOCATABLE :: SIGMA_L(:)
REAL(rkind) , ALLOCATABLE :: LFACTOR_L(:)
REAL(rkind) :: DSIGMA
REAL(rkind) :: REDUC
REAL(rkind) :: tau_normal
REAL(rkind) :: sign_new,sign_old
REAL(rkind) :: frat
REAL(rkind) :: err
REAL(rkind) :: RCHANGE
REAL(rkind) :: tau_total
REAL(rkind) :: TAUWLIM
REAL(rkind) :: U10MOD
REAL(rkind) :: freq_tmp
REAL(rkind) :: Cv,tauv
REAL(rkind), PARAMETER :: KAPPA = 0.4_rkind
INTEGER :: nf_old,nf_new
INTEGER :: iter
INTEGER :: slow_down
INTEGER :: isol
INTEGER :: IS
integer istat
TESTFL = .false.
nf_old=size(S_in,1)
! Initialize LFACTOR_S (for error catching)
LFACTOR_S = HUGE(LFACTOR_S)
ALLOCATE(S_in1D_S(nf_old), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_babanin, allocate error 4')
! To simplify calcs, use S_in1D(f)
S_in1D_S = sum(S_in,1)*DDIR_rad ! note that variable is S_in(ID,IS)
S_in1D_S = S_in1D_S*(2.0*PI) ! was in units of m2/(rad-Hz), so put in units m2/Hz
! 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
! allocate arrays on nf_new
ALLOCATE(S_in1D_L(nf_new), DF(nf_new), CINV_L(nf_new), SIGMA_L(nf_new), LFACTOR_L(nf_new), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_babanin, allocate error 5')
! extend freqs to nf_new
SIGMA_L=SIGMA_S
CINV_L=CINV_S
S_in1D_L=S_in1D_S
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
! For Sin, we use an f-2 approximation...Sin=Beta*E~sigma*(U/C)^2*E=f1*f2*f-5=f-2
S_in1D_L(IS)=S_in1D_L(nf_old)*(SIGMA_L(nf_old)/SIGMA_L(IS))**2
enddo
tau_total=(UFRIC(IP)**2)*rhoa
TAUTOT(IP)=tau_total
U10MOD=UFRIC(IP)*28.0_rkind
Cv=(-5.0e-5_rkind)*WIND10+1.1e-3_rkind
Cv=max(Cv,0.0_rkind) ! otherwise Cv<0 for U10>22 m/s
tauv=rhoa*Cv*(WIND10**2)
CD(IP)=(UFRIC(IP)/WIND10)**2
Z0(IP) = 10.0_rkind/EXP(KAPPA*WIND10 /UFRIC(IP))
ALPHA_CH(IP)=G9 * Z0(IP) /(UFRIC(IP)**2)
TAUWLIM=tau_total-tauv
DO IS=1,NF_new
DSIGMA=SIGMA_L(IS)*FRINTF
DF(IS)=DSIGMA/(2.0*PI)
END DO
! calculate without reduction
REDUC=0.0
call calc_tau_normal(tau_normal,Lfactor_L,REDUC,S_in1D_L,df,CINV_L,U10MOD,grav,rhow)
if(TESTFL)then
206 format('tau_total = ',f9.6,' ; tauv = ',f9.6, &
& ' ; tau_normal = ',f9.6,' ; TAUWLIM = ',f9.6)
write(DBG%FHNDL,206)tau_total,tauv,tau_normal,TAUWLIM
endif
if(tau_normal < TAUWLIM)then
DO IS=1,NF_OLD
Lfactor_S(IS)=1.0_rkind
END DO
RETURN
else
if(TESTFL)then
write(DBG%FHNDL,*)'reducing Sin to make tau_normal match TAUWLIM'
endif
endif
! if(tau_normal > 10.0) then
! write(DBG%FHNDL,*) ip, wind10, ufric, tau_normal
! end if
! calculate with reduction, start value arbitrary
REDUC=0.05
call calc_tau_normal(tau_normal,Lfactor_L,REDUC,S_in1D_L,df,CINV_L,U10MOD,grav,rhow)
err=tau_normal-TAUWLIM
sign_new=sign(1.0_rkind,err)
slow_down=0
RCHANGE=2.0
isol=0
do iter=1,500
if(tau_normal > TAUWLIM)then ! increase REDUC
REDUC=REDUC*RCHANGE
else ! then reduce REDUC
REDUC=REDUC/RCHANGE
endif
call calc_tau_normal(tau_normal,Lfactor_L,REDUC,S_in1D_L,df, &
& CINV_L,U10MOD,grav,rhow)
err=tau_normal-TAUWLIM
sign_old=sign_new
sign_new=sign(1.0_rkind,err)
if(TESTFL)then
207 format('iter = ',i3,' ; REDUC = ',f9.6,' ; RCHANGE = ', &
& f9.6,' ; tau_normal = ',f9.6,' ; err = ',f9.6)
write(DBG%FHNDL,207)iter,REDUC,RCHANGE,tau_normal,err
endif
if(abs(sign_new-sign_old) .gt. tiny(1.))then
RCHANGE=0.5_rkind*(1.0_rkind+RCHANGE)
endif
if((abs(err)/TAUWLIM) < 1.656e-06_rkind)then
isol=1
exit
endif
end do
if(isol==0)then
write(DBG%FHNDL,*)'no solution found at gridpoint', IP, SUM(AC2(:,:,IP))
write(DBG%FHNDL,'(A20,F15.8,A20,F15.8,A10,F15.8,A20,F15.8)')'tau_normal = ',tau_normal,' TAUWLIM =', TAUWLIM,'err =',err,'(abs(err)/TAUWLIM) = ',abs(err)/TAUWLIM
write(DBG%FHNDL,*) wind10, ufric(ip)
! if((abs(err)/TAUWLIM).ge.1.0_rkind) then
open(401,file='debug.dat',form='unformatted')
write(401)nf_new ! scalar
write(401)Lfactor_L ! LFACTOR_L(nf_new)
write(401)S_in1D_L ! S_in1D_L(nf_new)
write(401)df ! DF(nf_new)
write(401)CINV_L ! CINV_L(nf_new)
write(401)TAUWLIM ! scalar
write(401)tau_normal ! scalar
write(401)REDUC ! scalar
write(401)U10MOD ! scalar
write(401)rhow ! scalar
close(401)
stop 'wwm_babanin l.466'
! end if
end if
DO IS=1,NF_OLD
Lfactor_S(IS)=Lfactor_L(IS)
END DO
DEALLOCATE(S_in1D_S, S_in1D_L, DF, CINV_L, SIGMA_L, LFACTOR_L)
end subroutine calc_Lfactor
subroutine calc_tau_normal(tau_normal,Lfactor,REDUC,S_in1D,df,CINV, &
& U10MOD,grav,rhow)
! 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
use datapool, only : rkind
implicit none
REAL(rkind) , INTENT(OUT) :: tau_normal
REAL(rkind) , INTENT(OUT) :: Lfactor(:)
REAL(rkind) , INTENT(IN) :: REDUC
REAL(rkind) , INTENT(IN) :: S_in1D(:)
REAL(rkind) , INTENT(IN) :: df(:)
REAL(rkind) , INTENT(IN) :: CINV(:)
REAL(rkind) , INTENT(IN) :: U10MOD
REAL(rkind) , INTENT(IN) :: grav
REAL(rkind) , INTENT(IN) :: rhow
REAL(rkind) , ALLOCATABLE :: A(:)
REAL(rkind) , ALLOCATABLE :: UoverC(:)
REAL(rkind) , ALLOCATABLE :: S_in1D_red(:)
INTEGER :: nf
INTEGER :: IS
integer istat
nf=size(df)
ALLOCATE(UoverC(nf), A(nf), S_in1D_red(nf), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_babanin, allocate error 6')
UoverC=U10MOD * CINV ! matrix op
A=exp((1-UoverC) * REDUC) ! matrix op
do IS=1,nf
Lfactor(IS)=min(1.0_rkind,A(IS)) ! cannot be done as matrix op?
enddo
S_in1D_red=S_in1D * Lfactor ! matrix op
tau_normal=0.0_rkind
do IS=1,nf
tau_normal=tau_normal+rhow*grav*S_in1D_red(IS)*CINV(IS)*df(IS)
enddo
DEALLOCATE(UoverC, A, S_in1D_red)
end subroutine calc_tau_normal
!****************************************************************
SUBROUTINE SSWELL (IP,ETOT,ACLOC,IMATRA, IMATDA,URMSTOP,CGo)
! note that CGo is for diagnostic purposes only, may be excluded
! excluded : SPCDIR AC2 DEP2 IMATRA
!****************************************************************
USE DATAPOOL, ONLY : RKIND, RHOAW, RHOW, SPSIG, G9, WK, THR, MSC, SIGPOW
IMPLICIT NONE
INTEGER, INTENT(IN) :: IP
REAL(rkind) , INTENT(IN) :: ETOT
REAL(rkind) , INTENT(IN) :: CGo(:,:) ! CGo(MSC,MICMAX) ! group velocity without currents, but includes depth effects
REAL(rkind) , INTENT(IN) :: URMSTOP
REAL(rkind) , INTENT(INOUT) :: IMATDA(:,:) ! IMATDA(MDC,MSC)
REAL(rkind) , INTENT(INOUT) :: IMATRA(:,:) ! IMATDA(MDC,MSC)
REAL(rkind) , INTENT(IN) :: ACLOC(:,:) ! ACLOC(MDC,MSC)
! local constants:
REAL(rkind), parameter :: feswell=0.0035_rkind ! "fe is of the order 0.002 to 0.008" (Fabrice proposes 0.0035)
REAL(rkind), parameter :: Re_crit=1.0e+5_rkind ! Re_c or SWELLF4
REAL(rkind), parameter :: Cdsv=1.2_rkind ! C_dsv or SWELLF5
REAL(rkind), parameter :: nu_air=1.4E-5_rkind
! local variables:
INTEGER :: IS
INTEGER :: ID
REAL(rkind) :: aorb,USIGTOP
REAL(rkind) :: Re
REAL(rkind) :: SWDIS(MSC) ! this is S_SWELL / E(f,theta)
! REAL(rkind) :: myu,SWDIS_CHECK ! error checking, may be omitted
! Method:
! In general....
! @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)
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, MSC
SWDIS(IS) = RHOAW * 16.0 * feswell * SIGPOW(IS,2) & ! note that "-1" omitted since LHS
& * URMSTOP / g9 ! units of radian^2/s (radians/sec is conventional)
END DO
else
DO IS=1, MSC
SWDIS(IS) = RHOAW * Cdsv * 2.0 * WK(IS,IP) & ! note that "-1" omitted since LHS
& * SQRT(2.0 * NU_AIR * SPSIG(IS)) ! units of radian^1.5/s (radians/sec is conventional)
! start error check block (may be omitted)
! myu=2.0*RHOAW*(1.0/g9)*(1.0/CGo(IS,1))
! & *(SPCSIG(IS)**2.5)*sqrt(2.0*NU_AIR)
! write(DBG%FHNDL,*)'myu(',IS,') = ',myu
! SWDIS_CHECK=Cdsv*myu*CGo(IS,1)
! write(DBG%FHNDL,*)'SWDIS(IS),SWDIS_CHECK = ',
! & SWDIS(IS),SWDIS_CHECK
! end error check block (may be omitted)
END DO
endif
DO IS=1, MSC
DO ID = 1, MSC
IMATDA(IS,ID) = IMATDA(IS,ID) + SWDIS(IS)
IF (ACLOC(IS,ID) .GT. THR) IMATRA(IS,ID) = SWDIS(IS) * ACLOC(IS,ID)
END DO
END DO
RETURN
END SUBROUTINE SSWELL
SUBROUTINE SWIND_DBYB (IP,WIND10,THETAW,IMATRA,SSINE) ! 30.21
use datapool
implicit none
INTEGER, INTENT(IN) :: IP
INTEGER :: ID, IS
REAL(rkind) :: THETAW, SIGMA, SWINEB, CTW, STW, COSDIF, TEMP2, UoverC, WIND10
!
REAL(rkind) :: S_IN(MSC,MDC), IMATRA(MSC,MDC)
REAL(rkind), INTENT(OUT) :: SSINE(MSC,MDC)
!
REAL(rkind) :: DMAX,AINV
REAL(rkind) :: KTHETA(MSC,MDC) ! like D(theta), except max value at each freq is unity
REAL(rkind) :: ASPR(MSC), SIGDENS(MSC), SQRTBN(MSC), CINV(MSC), LFACTOR(MSC)
REAL(rkind) :: GAMMAD, GDONEL, TEMP4, TEMP42, TEMP5, TEMP6, BN
LOGICAL :: TESTFL
!WER REAL(rkind) HM0,ETOT,EDENS2D ! FOR test calcs (remove later)
CTW = COS(THETAW) ! 40.41
STW = SIN(THETAW) ! 40.41
IF (WIND10 .LT. VERYSMALL) RETURN
DO IS = 1, MSC
SIGDENS(IS) = 0.
DO ID = 1, MDC
SIGDENS(IS) = SIGDENS(IS) + SPSIG(IS) * AC2(IS,ID,IP)
KTHETA(IS,ID) = AC2(IS,ID,IP)
END DO
SIGDENS(IS)=SIGDENS(IS)*DDIR
END DO
! 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 ! 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
! DO ID = 1, MDC
! KTHETA(IS,ID)=KTHETA(IS,ID)/DMAX
! END DO
! END DO
DO IS = 1, MSC
DMAX = 0.
DMAX = MAXVAL(KTHETA(IS,:))
IF (DMAX .LT. 10.E-10) THEN
KTHETA(IS,:) = 1.
ELSE
KTHETA(IS,:)=KTHETA(IS,:)/DMAX
END IF
END DO
DO IS = 1, MSC
AINV=0.0_rkind
DO ID = 1, MDC
AINV=AINV+KTHETA(IS,ID)*DDIR
END DO
ASPR(IS)=1./AINV
BN=ASPR(IS)*SIGPOW(IS,5)*SIGDENS(IS)/(2.*G9**2)
SQRTBN(IS)=SQRT(BN)
END DO
TEMP2 = 28.0_rkind * UFRIC(IP)
S_IN=0.0_rkind
DO IS = 1, MSC
SIGMA = SPSIG(IS)
CINV(IS) = WK(IS,IP) / SIGMA
UoverC = TEMP2 * CINV(IS)
DO ID=1,MDC
IF ( WIND10 .GT. THR ) THEN
COSDIF = COSTH(ID)*CTW+SINTH(ID)*STW
TEMP4=( UoverC * COSDIF - 1.0_rkind )
TEMP4=MAX(0.0_rkind,TEMP4)
TEMP42=TEMP4**2
TEMP5=10.0_rkind*SQRTBN(IS)*TEMP42-11.0_rkind
TEMP6=1.0_rkind+MyTANH(TEMP5)
GDONEL=2.8_rkind-TEMP6
GAMMAD=GDONEL*SQRTBN(IS)*TEMP42
SWINEB = GAMMAD * SIGMA * RHOAW ! NEW SWINEB
S_IN(IS,ID)=SWINEB*AC2(IS,ID,IP)*SPSIG(IS)
END IF
ENDDO
ENDDO
call calc_Lfactor(ip,Lfactor,S_in,DDIR,SPSIG,FRINTF,CINV,G9,WIND10,TESTFL)
SSINE = S_IN
DO IS = 1, MSC
DO ID = 1, MDC
IF ( WIND10 .GT. VERYSMALL ) THEN
S_IN(IS,ID)= S_IN(IS,ID) * Lfactor(IS)
IMATRA(IS,ID) = IMATRA(IS,ID) + S_IN(IS,ID)/SPSIG(IS)
SWINEB=S_IN(IS,ID)/(AC2(IS,ID,IP)*SPSIG(IS))
END IF
ENDDO
ENDDO
END SUBROUTINE SWIND_DBYB
END MODULE SdsBabanin