#include "wwm_functions.h"
!**********************************************************************
!* Yasser Eldeberky, Nonlinear Transformation of Wave Spectra in *
!* the Nearshore Zone, PhD thesis, TU Delft *
!**********************************************************************
subroutine triad_eldeberky(ip, hs, smespc, acloc, imatra, imatda, ssnl3, dssnl3)
use datapool
implicit none
integer, intent(in) :: ip
real(rkind), intent(in) :: hs, smespc
real(rkind), intent(out) :: ssnl3(msc,mdc), dssnl3(msc,mdc)
real(rkind), intent(in) :: acloc(msc,mdc)
real(rkind), intent(inout) :: imatra(msc,mdc), imatda(msc,mdc)
integer i1, i2, id, is, ismax, ij1, ij2
real(rkind) aux1, aux2, biph, c0, cm, dep_2, dep_3, e0
real(rkind) em,ft, rint, sigpi, sinbph, stri
real(rkind) w0, wm, wn0, wnm, ursell, c1, c2, c3
real(rkind) eCont
real(rkind) :: E(MSC)
real(rkind) :: SA(1:MSC+TRI_ISP1,1:MDC)
ssnl3 = ZERO
dssnl3 = ZERO
IF (HS .LT. SMALL) RETURN
CALL URSELL_NUMBER(HS,SMESPC,DEP(IP),URSELL)
IF ( URSELL .le. TRI_ARR(5) ) RETURN
E = 0.
SA = 0.
ISMAX = TRI_ISP1
DO IS = TRI_ISP1, MSC
IF ( SPSIG(IS) .LT. ( TRI_ARR(2) * SMESPC) ) THEN
ISMAX = IS
ENDIF
ENDDO
c1 = G9*(DEP(IP)**2)
c2 = (TWO/15._rkind)*G9*(DEP(IP)**4)
c3 = (TWO/5._rkind)*(DEP(IP)**3)
BIPH = PIHALF*(MyTANH(TRI_ARR(4)/URSELL)-1.)
SINBPH = ABS( SIN(BIPH) )
DO ID = 1, MDC
E = ACLOC(:,ID) * PI2 * SPSIG
DO IS=TRI_ISBEGIN, ISMAX
E0 = E(IS)
W0 = SPSIG(IS)
WN0 = WK(IS,IP)
C0 = W0 / WN0
EM = TRI_WISM * E(IS+TRI_ISM1) + TRI_WISM1 * E(IS+TRI_ISM)
WM = TRI_WISM * SPSIG(IS+TRI_ISM1) + TRI_WISM1 * SPSIG(IS+TRI_ISM)
WNM = TRI_WISM * WK(IS+TRI_ISM1,IP) + TRI_WISM1 * WK(IS+TRI_ISM,IP)
CM = WM / WNM
AUX1 = WNM**2 * ( G9 * DEP(IP) + TWO*CM**2 )
AUX2 = WN0 * ( c1 + c2 * WN0**2 - c3 * W0**2) ! (m/s² * m + m/s² * m³*1/m² - 1/s² * m²)
RINT = AUX1 / AUX2
FT = TRI_ARR(1) * C0 * CG(IS,IP) * RINT**2 * SINBPH
SA(IS,ID) = MAX(ZERO, FT*(EM*(EM - 2*E0)))
END DO
END DO
DO ID = 1, MDC
DO IS = 1, MSC
SIGPI = SPSIG(IS) * PI2
IF (ACLOC(IS,ID) .LT. THR) CYCLE
STRI = SA(IS,ID) - TWO*(TRI_WISP*SA(IS+TRI_ISP1,ID) + TRI_WISP1*SA(IS+TRI_ISP,ID))
IF (ABS(STRI) .LT. THR) CYCLE
eCont = STRI / SIGPI
IF (ICOMP .GE. 2) THEN
IF (STRI .GT. 0.) THEN
IMATRA(IS,ID) = IMATRA(IS,ID) + eCont
SSNL3(IS,ID) = eCont
ELSE
IMATDA(IS,ID) = IMATDA(IS,ID) - eCont / ACLOC(IS,ID)
DSSNL3(IS,ID) = - eCont/ACLOC(IS,ID)
END IF
ELSE
IMATRA(IS,ID) = IMATRA(IS,ID) + eCont
IMATDA(IS,ID) = ZERO !IMATDA(IS,ID) + STRI / (ACLOC(IS,ID)*SIGPI)
SSNL3(IS,ID) = eCont
DSSNL3(IS,ID) = eCont / ACLOC(IS,ID)
END IF
END DO
END DO
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE triad_dingemans (ip, acloc, imatra, imatda, ssnl3)
use datapool
implicit none
integer, intent(in) :: ip
real(rkind), intent(in) :: acloc(msc,mdc)
real(rkind), intent(inout) :: imatra(msc,mdc), imatda(msc,mdc)
real(rkind), intent(out) :: ssnl3(msc,mdc)
integer :: is, is2, id
real(rkind) :: ecloc(msc,mdc), e2(msc,mdc), d20
real(rkind) :: df, domega, omega, omega1, fac, z1a, z1b
integer :: j1, j2, jmin, j2abs
do is = 1, msc
do id = 1, mdc
ecloc(is,id) = acloc(is,id) * spsig(is) * ddir
end do
end do
SSNL3 = 0.
do id = 1, mdc
do is = 1,msc-1
df = ((spsig(is+1) - spsig(is)))/pi2
domega = pi2 * df
omega = spsig(is) * pi2
jmin = nint(0.5*MyREAL(is))
if (2*jmin .eq. is) then
fac = 0.5
else
fac = 1.
endif
z1a = dep(ip) * (jmin*domega)**2 / g9
z1b = dep(ip) * ((jmin-is-1)*domega)**2 / g9
e2(is,id) = 0.
do is2 = jmin, msc
j1 = is2
j2 = is2 - is
j2abs = iabs (j2)
! zero frequencies are skipped
if (j2 .eq. 0) cycle
omega1 = is2 * domega
e2(is,id) = e2(is,id) + fac * d20(omega,omega1,dep(ip),z1a,z1b)* ecloc(j1,id)*ecloc(j2abs,id)
fac = 1.
end do
e2(is,id) = e2(is,id) * df
end do
end do
end subroutine
!**********************************************************************
!* *
!**********************************************************************
function d20(omega, omega1, h, z1a, z1b)
use datapool, only : g9, rkind
implicit none
real(rkind), intent(in) :: omega, omega1, h
real(rkind), intent(inout) :: z1a, z1b
real(rkind) :: k, k1, k2
real(rkind) :: omega2
real(rkind) :: aome2
real(rkind) :: cthkh, cthk1h, cthk2h
real(rkind) :: d2, d20
omega2 = omega - omega1
aome2 = abs (omega2)
z1a = abs(z1a)
z1b = abs(z1b)
call dispu2 (omega1, h, z1a, k1)
call dispu2 (aome2, h, z1b, k2)
if (omega2 .lt. 0.) k2 = - k2
k = k1 + k2
cthk1h = 1. / MyTANH (k1*h)
cthk2h = 1. / MyTANH (k2*h)
cthkh = 1. / MyTANH (k*h)
d2 = 0.5 * (omega1**2 + omega2**2 + omega1*omega2 - &
& omega1 * omega2 * cthk1h * cthk2h - omega * &
& omega1 * cthkh * cthk1h - omega * omega2 * &
& cthkh * cthk2h)
d20 = d2 / (g9 * (1. - omega**2 * cthkh / (g9*k)))
d20 = 2. * d20**2
end function
!**********************************************************************
!* *
!**********************************************************************
subroutine dispu2(omega,d,z1,k)
use datapool, only : g9, rkind
implicit none
real(rkind), intent(in) :: omega, d
real(rkind), intent(inout) :: z1
real(rkind), intent(out) :: k
real(rkind), parameter :: eps = 0.0001
real(rkind) :: z0, z2, fak1, fak2, sig
z0 = d*omega*omega/g9
10 sig = MyTANH(z1)
fak1 = z1*sig
fak2 = z1 + sig*(1.-fak1)
z2 = z1 + (z0-fak1)/fak2
if (abs((z2-z1)/z2).gt.eps) goto 40
goto 60
40 z1 = z2
goto 10
60 k = z2/d
z1 = z2
end subroutine
!**********************************************************************
!* *
!**********************************************************************
function delta(ip, is, is1, is2) result(res)
use datapool, only : rkind, wk
implicit none
integer, intent(in) :: ip, is, is1, is2
real(rkind) :: res !return
res = wk(ip,is) - wk(ip,is1) - wk(ip,is2)
end function delta
!**********************************************************************
!* *
!**********************************************************************
function ddelta_dx(ip, is, is1, is2, id) result(res)
use datapool, only : rkind, sinth, costh, dwkdx, dwkdy
implicit none
integer, intent(in) :: ip, is, is1, is2, id
real(rkind) :: res
real(rkind) :: ka_1abl, kb_1abl, kc_1abl
ka_1abl = costh(id) * dwkdx(ip,is) + sinth(id) * dwkdy(ip,is)
kb_1abl = costh(id) * dwkdx(ip,is1) + sinth(id) * dwkdy(ip,is1)
kc_1abl = costh(id) * dwkdx(ip,is2) + sinth(id) * dwkdy(ip,is2)
res = ka_1abl - kb_1abl - kc_1abl
end function ddelta_dx
!**********************************************************************
!* *
!**********************************************************************
function w(ip, is, is1, is2, n1, emf) result(res)
use datapool, only : rkind, wk, cg, spsig, g9, one
implicit none
integer, intent(in) :: ip, is, is1, is2, n1, emf
real(rkind) :: res
real(rkind) :: omegaa
real(rkind) :: omegab
real(rkind) :: omegac
real(rkind) :: cga
real(rkind) :: cgb
real(rkind) :: cgc
real(rkind) :: ka
real(rkind) :: kb
real(rkind) :: kc
real(rkind) :: tau
real(rkind) :: a,b,c,d,e
! n1=0 means +
! n1=1 means -
! em=0 without eldeberky & madsen
! em=1 with eldeberky & madsen
omegaa = spsig(is)
omegab = spsig(is)
omegac = spsig(is)
cga = cg(ip,is)
cgc = cg(ip,is1)
cgb = cg(ip,is2)
ka = wk(ip,is)
kb = wk(ip,is1)
kc = wk(ip,is2)
!AR: reduce the +- stuff ...
!
res = one / (8 * sqrt(cga * cgb *cgc)) * &
& ( (-1)**n1 * (2 - emf * tau(ip, is, is1, is2, n1)) * kb * kc + &
& ( 1 - emf * tau(ip, is, is1, is2, n1)) * ((omegab * omegac)**2 / g9**2) + &
& ( kb**2 * omegac / omegaa) + ( (-1)**n1 * kc**2 * omegab / omegaa) + ( (-1)**(n1+1) * &
& ( (-1)**(n1+1) * ( 1- emf * tau(ip, is, is1, is2,n1)) * ((omegaa**2 * omegab * omegac) / g9**2) ) ) )
a = one / (8 * sqrt(cga * cgb *cgc))
b = (-1)**n1 * (2 - emf * tau(ip, is, is1, is2, n1)) * kb * kc
c = ( 1 - emf * tau(ip, is, is1, is2, n1)) * ((omegab * omegac)**2 / g9**2)
d = ( kb**2 * omegac / omegaa) + ( (-1)**n1 * kc**2 * omegab / omegaa)
e = ( (-1)**(n1+1) * ( (-1)**(n1+1) * ( 1- emf * tau(ip, is, is1, is2,n1)) * ((omegaa**2 * omegab * omegac) / g9**2) ))
end function w
!**********************************************************************
!* *
!**********************************************************************
function tau(ip, is, is1, is2, n1) result(res)
use datapool, only : rkind, wk, cg, spsig
implicit none
integer, intent(in) :: ip, is, is1, is2, n1
real(rkind) :: res
res = 2 * cg(ip,is) * ( wk(ip,is) + (-1)**(n1+1) * wk(ip,is1) - wk(ip,is2) ) / spsig(is)
end function tau
!**********************************************************************
!* *
!**********************************************************************
function dwdx(ip, is, is1, is2, id, n1, switch) result(res)
use datapool, only : rkind, wk, cg, spsig, g9, one, costh, sinth, dcgdx, dcgdy, dwkdx, dwkdy, one, two, three
implicit none
integer, intent(in) :: ip, is, is1, is2, id, n1
real(rkind) :: res
integer :: switch
real(rkind) :: omegaa, omegab, omegac
real(rkind) :: cga, cgb, cgc, cga_, cgb_, cgc_
real(rkind) :: ka, kb, kc, ka_, kb_, kc_
omegaa = spsig(is)
omegab = spsig(is)
omegac = spsig(is)
cga = cg(ip,is)
cgc = cg(ip,is1)
cgb = cg(ip,is2)
ka = wk(ip,is)
kb = wk(ip,is1)
kc = wk(ip,is2)
cga_ = costh(id) * dcgdx(ip,is) + sinth(id) * dcgdy(ip,is)
cgb_ = costh(id) * dcgdx(ip,is1) + sinth(id) * dcgdy(ip,is1)
cgc_ = costh(id) * dcgdx(ip,is2) + sinth(id) * dcgdy(ip,is2)
ka_ = costh(id) * dwkdx(ip,is) + sinth(id) * dwkdy(ip,is)
kb_ = costh(id) * dwkdx(ip,is1) + sinth(id) * dwkdy(ip,is1)
kc_ = costh(id) * dwkdx(ip,is2) + sinth(id) * dwkdy(ip,is2)
res = (((4*Cga*Cgb*Cgc*((-1)**n1*((kb_*kc + kb*kc_)*omegaa +kc*kc_*omegab) +kb*kb_*omegac))/&
omegaa -(Cga*Cgb_*Cgc +Cgb*(Cga_*Cgc + Cga*Cgc_))*(2*(-1)**n1*kb*kc -((-1)**n1*omegaa**2*omegab*omegac)/g9**2+&
(omegab**2*omegac**2)/g9**2 +((-1)**n1*kc**2*omegab +kb**2*omegac)/omegaa)))/(16.*(Cga*Cgb*Cgc)**1.5)
end function dwdx
!**********************************************************************
!* *
!**********************************************************************
function k(ip, is, is1, is2, id, is3, is4, is5, n1, emf) result(res)
use datapool, only : rkind, g9, one, stat
implicit none
integer, intent(in) :: ip, is, is1, is2, is3, is4, is5, id, n1, emf
real(rkind) :: res
real(rkind) :: delta, dwdx, w, ddelta_dx
res = one / ((delta(ip, is3, is4, is5))**2) * dwdx(ip,is,is1,is2,id,n1,emf) - (w(ip,is,is1,is2,n1,emf) / ((delta(ip, is, is1, is2))**3)) * ddelta_dx(ip, is3, is4, is5, id)
end function k
!**********************************************************************
!* *
!**********************************************************************
function j(ip, is, is1, is2, id) result(res)
use datapool, only : rkind, one, stat
implicit none
integer, intent(in) :: ip, is, is1, is2, id
real(rkind) :: res
!!! function declaration
real(rkind) :: delta, ddelta_dx
res = - ( one / (delta(ip, is, is1, is2)**3) ) * ddelta_dx(ip, is, is1, is2, id)
end function j
!**********************************************************************
!* *
!**********************************************************************
subroutine snl3ta(ip,snl3,dsnl3)
use datapool, only : rkind, msc, mdc, ac2, ZERO, spsig, cg, frintf, ddir, fr, stat
implicit none
real(rkind), intent(out) :: snl3(msc,mdc), dsnl3(msc,mdc)
integer, intent(in) :: ip
integer :: is, is1, is2, id
real(rkind) :: SUPER, SUB, f, f1, f2, SUPERD, SUBD, k, w, JAC
do id = 1, mdc
snl3(:,id) = zero
dsnl3(:,id) = zero
#ifdef SNL3_DOUBLE
do is = 1, msc
SUPER = ZERO; SUPERD = ZERO
SUB = ZERO; SUBD = ZERO
f = (ac2(ip,is,id) * spsig(is)**2. * frintf * ddir)* cg(ip,is)
do is1 = 1, msc
f1 = (ac2(ip,is1,id) * spsig(is1)**2. * frintf * ddir) * cg(ip,is1)
do is2 = 1, msc
f2 = (ac2(ip,is2,id) * spsig(is2)**2. * frintf * ddir) * cg(ip,is2)
SUPER = SUPER + ( k(ip, is, is1, is2, id, is, is1, is2, 0, 0) * f1 * f2 &
& + k(ip, is1, is2, is, id, is, is1, is2, 1, 0) * f1 * f &
& + k(ip, is2, is1, is, id, is, is1, is2, 1, 0) * f2 * f ) * w(ip, is, is1, is2, 0, 0) * kron_delta(is,is1+is2)
SUPERD = SUPERD + ( k(ip, is1, is2, is, id, is, is1, is2, 1, 0) * f1 + &
& k(ip, is2, is1, is, id, is, is1, is2, 1, 0) * f2 ) * w(ip, is, is1, is2, 0, 0) * kron_delta(is,is1+is2)
enddo
enddo
SUPER = 4 * SUPER
do is1 = 1, msc
f1 = ac2(ip,is1,id)
do is2 = 1, msc
f1 = ac2(ip,is2,id)
SUB = SUB + ( k(ip, is, is1, is2, id, is2, is, is1, 1, 0) * f1 * f2 + k(ip, is1, is, is2, id, is2, is, is1, 1, 0) * f2 * f + &
& k(ip, is2, is1, is, id, is2, is, is1, 0, 0) * f1 * f ) * w(ip, is, is1, is2, 1, 0) * kron_delta(is2,is+is1)
SUBD = SUBD + ( k(ip, is, is1, is2, id, is2, is, is1, 1, 0) * f1 * f2 + k(ip, is1, is, is2, id, is2, is, is1, 1, 0) * f2 * f + &
& k(ip, is2, is1, is, id, is2, is, is1, 0, 0) * f1 * f ) * w(ip, is, is1, is2, 1, 0) * kron_delta(is2,is+is1)
enddo
enddo
#else
do is = 2, msc
SUPER = ZERO; SUPERD = ZERO
f = (ac2(ip,is,id) * spsig(is)**2. * frintf * ddir) * cg(ip,is)
do is1 = 1, is - 1
is2 = is - is1
f1 = (ac2(ip,is1,id) * spsig(is1)**2. * frintf * ddir) * cg(ip,is1)
f2 = (ac2(ip,is2,id) * spsig(is2)**2. * frintf * ddir) * cg(ip,is2)
SUPER = SUPER + ( k(ip, is , is1, is2, id, is, is1, is2, 0, 0) * f1 * f2 &
& + k(ip, is1, is2, is , id, is, is1, is2, 1, 0) * f1 * f &
& + k(ip, is2, is1, is , id, is, is1, is2, 1, 0) * f2 * f ) * w(ip, is, is1, is2, 0, 0)
SUPERD = SUPERD + ( k(ip, is1, is2, is , id, is, is1, is2, 1, 0) * f1 &
& + k(ip, is2, is1, is , id, is, is1, is2, 1, 0) * f2 ) * w(ip, is, is1, is2, 0, 0)
enddo
SUPER = 4 * SUPER
SUPERD = 4 * SUPERD
JAC = 1./(SPSIG(IS)*DDIR*FRINTF)
snl3(is,id) = SUPER * JAC
dsnl3(is,id) = SUPERD
end do
do is = 1, msc-1
SUB = ZERO; SUBD = ZERO
f = (ac2(ip,is,id) * spsig(is)**2. * frintf * ddir) * cg(ip,is)
do is1 = 1, msc - is
is2 = is + is1
f1 = (ac2(ip,is1,id) * spsig(is1)**2. * frintf * ddir) * cg(ip,is1)
f2 = (ac2(ip,is2,id) * spsig(is2)**2. * frintf * ddir) * cg(ip,is2)
SUB = SUB + ( k(ip, is , is1, is2, id, is2, is, is1, 1, 0) * f1 * f2 &
& + k(ip, is1, is , is2, id, is2, is, is1, 1, 0) * f2 * f + &
& k(ip, is2, is1, is , id, is2, is, is1, 0, 0) * f1 * f ) * w(ip, is, is1, is2, 1, 0)
SUBD = SUBD + ( k(ip, is , is1, is2, id, is2, is, is1, 1, 0) * f1 * f2 &
& + k(ip, is1, is , is2, id, is2, is, is1, 1, 0) * f2 * f &
& + k(ip, is2, is1, is , id, is2, is, is1, 0, 0) * f1 * f ) * w(ip, is, is1, is2, 1, 0)
enddo ! is1
SUB = 8 * SUB
SUBD = 8 * SUBD
JAC = 1./(SPSIG(IS)*DDIR*FRINTF)
snl3(is,id) = snl3(is,id) + SUB * JAC
dsnl3(is,id) = dsnl3(is,id) + SUBD
#endif
enddo
enddo
end subroutine snl3ta
!**********************************************************************
!* *
!**********************************************************************
function kron_delta(i, j) result(res)
use datapool, only : stat
implicit none
integer, intent(in) :: i,j
integer :: res
res = int((float(i+j)-abs(i-j)))/(float((i+j)+abs(i-j)))
end function kron_delta
!**********************************************************************
!* *
!**********************************************************************
subroutine snl31 (ip, hs, smespc, acloc, imatra, imatda, ssnl3)
use datapool
implicit none
integer, intent(in) :: ip
real(rkind), intent(in) :: hs, smespc
real(rkind), intent(out) :: ssnl3(msc,mdc)
real(rkind), intent(in) :: acloc(msc,mdc)
real(rkind), intent(inout) :: imatra(msc,mdc), imatda(msc,mdc)
INTEGER :: IS, ID, I1, I2
REAL(rkind) :: BIPH, ASINB, XISTRI
REAL(rkind) :: ED(MSC)
INTEGER :: IS1, IS2, ISMAX
INTEGER :: IRES
REAL(rkind) :: AUX1, AUX2, FAC
REAL(rkind) :: JAC, DNL3IS1, DNL3IS2
REAL(rkind) :: NL3IS1, NL3IS2
REAL(rkind) :: cgl(msc),cl(msc),wkl(msc), AA
REAL(rkind) :: URSELL, BB, DEP1, DEP2, DEP3
PTRIAD(1) = 1.
PTRIAD(2) = 2.5
PTRIAD(3) = 10.
PTRIAD(4) = 0.2
PTRIAD(5) = 0.01
IF (TRICO .GT. 0.) PTRIAD(1) = TRICO
IF (TRIRA .GT. 0.) PTRIAD(2) = TRIRA
IF (TRIURS .GT. 0.) PTRIAD(5) = TRIURS
ssnl3 = 0
CALL URSELL_NUMBER(HS,SMESPC,DEP(IP),URSELL)
IF (URSELL .le. PTRIAD(5)) RETURN
BB = ONE/15._rkind
AA = TWO/FIVE
DEP1 = DEP(IP)
DEP2 = DEP1**2
DEP3 = DEP2*DEP1
wkl = wk(ip,:)
cgl = cg(ip,:)
cl = wc(ip,:)
IRES = NINT(LOG(TWO)/LOG(XIS))
ISMAX = 1
DO IS = 1, MSC
IF (SPSIG(IS) < (PTRIAD(2)*SMESPC)) ISMAX = IS
END DO
ISMAX = MAX(ISMAX,IRES+1)
BIPH = PI/TWO*(MyTANH(0.2/URSELL)-ONE)
ASINB = ABS(SIN(BIPH))
DO ID = 1, MDC
ED = ACLOC(:,ID)*SPSIG
DO IS = 1, ISMAX-IRES
IS1 = IS+IRES
IS2 = IS
AUX1 = WKL(IS2)**2*(G9*DEP1+TWO*CL(IS2)**2)
AUX2 = WKL(IS1)*DEP1*(G9*DEP1+TWO*BB*G9*DEP3*WKL(IS1)**2-AA*SIGPOW(IS,2)*DEP2)
JAC = AUX1/AUX2
FAC = PTRIAD(1)*PI2*CL(IS1)*CGL(IS1)*JAC**2*ASINB
NL3IS1 = MAX(ZERO,FAC*(ED(IS2)*ED(IS2)-TWO*ED(IS2)*ED(IS1)))
DNL3IS1 = MAX(ZERO,FAC*(ED(IS2)-TWO*ED(IS1)))
NL3IS2 = TWO*NL3IS1
DNL3IS2 = TWO*DNL3IS1
IF (ICOMP .LT. 2) THEN
IMATRA(IS1,ID) = IMATRA(IS1,ID) + NL3IS1 / SPSIG(IS1)
IMATRA(IS2,ID) = IMATRA(IS2,ID) - NL3IS2 / SPSIG(IS2)
IMATDA(IS1,ID) = IMATDA(IS1,ID) + DNL3IS1
IMATDA(IS2,ID) = IMATDA(IS2,ID) - DNL3IS2
ELSE
IMATRA(IS1,ID) = IMATRA(IS1,ID) + NL3IS1 / SPSIG(IS1)
IMATRA(IS2,ID) = IMATRA(IS2,ID) - NL3IS2 / SPSIG(IS2)
IMATDA(IS1,ID) = IMATDA(IS1,ID) - DNL3IS1
IMATDA(IS2,ID) = IMATDA(IS2,ID) + DNL3IS2
ENDIF
SSNL3(IS1,ID) = SSNL3(IS1,ID) + NL3IS1
SSNL3(IS2,ID) = SSNL3(IS2,ID) - NL3IS2
END DO
END DO
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
subroutine snl32 (ip, hs, smespc, acloc, imatra, imatda, ssnl3)
use datapool
implicit none
integer, intent(in) :: ip
real(rkind), intent(in) :: hs, smespc
real(rkind), intent(out) :: ssnl3(msc,mdc)
real(rkind), intent(in) :: acloc(msc,mdc)
real(rkind), intent(inout) :: imatra(msc,mdc), imatda(msc,mdc)
INTEGER :: I, J, I1, I2
INTEGER :: IRES, ISMAX
INTEGER :: IS, ID
REAL(rkind) :: PTTRIAD(5)
REAL(rkind) :: DEP_2, DEP_3, BB, BIPH
REAL(rkind) :: TMN, URS
REAL(rkind) :: WiT,WjT,WNi,WNj,CGi,CGj,JACi,JACj,XISTRI
REAL(rkind) :: ALPH, BETA, SINBPH, FT, RHV_i, RHV_j, DIA_i, DIA_j
REAL(rkind) :: E(MSC), URSELL
LOGICAL :: TRIEXP
PTTRIAD(1) = 0.1
PTTRIAD(2) = 2.2
PTTRIAD(3) = 10.
PTTRIAD(4) = 0.2
PTTRIAD(5) = 0.01
IF (TRICO .GT. 0.) PTTRIAD(1) = TRICO
IF (TRIRA .GT. 0.) PTTRIAD(2) = TRIRA
IF (TRIURS .GT. 0.) PTTRIAD(5) = TRIURS
ssnl3=0
CALL URSELL_NUMBER(HS,SMESPC,DEP(IP),URSELL)
URS = MIN ( URSELL , TEN )
IF ( URS .lt. PTTRIAD(5) ) RETURN
DEP_2 = DEP(IP)**2
DEP_3 = DEP(IP)**3
BB = 1. / 15.
TRIEXP = .FALSE.
IRES = NINT ( LOG(TWO) / XISLN )
ISMAX = 1
DO IS = 1, MSC
IF ( SPSIG(IS) .LT. ( PTTRIAD(2) * SMESPC) ) THEN
ISMAX = IS
ENDIF
ENDDO
ISMAX = MAX ( ISMAX , IRES + 1 )
TMN = PI2 / SMESPC
IF (URS .LT. PTTRIAD(5)) THEN
BIPH = 0.0
ELSE
BIPH = - PI2 / 8. * ( LOG10(URS) + 1.)
ENDIF
SINBPH = ABS( SIN(BIPH) )
DO ID = 1, MDC
DO IS = 1, MSC
E(IS) = ACLOC(IS,ID) * PI2 * SPSIG(IS)
ENDDO
DO I = 1, ISMAX-IRES
J = I + IRES
WiT = SPSIG(I)
WjT = SPSIG(J)
WNi = WK(I,IP)
WNj = WK(J,IP)
CGi = CG(I,IP)
CGj = CG(J,IP)
JACi = PI2 * WiT
JACj = PI2 * WjT
ALPH = 4. * WNi**2 * ( 0.5 + ( WiT**2 / ( WNi**2 * G9 * DEP(IP) ) ) )
BETA = -2. * WNj * ( G9 * DEP(IP) + 2. * BB * G9 * DEP_3 * WNj**2 - ( BB + 1./3. ) * WjT**2 * DEP_2 )
FT = PTTRIAD(1) * CGj * SINBPH * ( G9 * ALPH / BETA )**2
RHV_i = 0.
RHV_j = 0.
DIA_i = 0.
DIA_j = 0.
IF ( TRIEXP ) THEN
RHV_j = FT * ( (WjT/WNj) * E(i) * E(i) - 2. * (WiT/WNi) * E(j) * E(i) )
RHV_i = RHV_j
IF ( RHV_j .LE. 0. ) THEN
RHV_j = 0.
RHV_i = 0.
ENDIF
IMATRA(i,ID) = IMATRA(i,ID) - RHV_i / JACi
IMATRA(j,ID) = IMATRA(j,ID) + 0.5 * RHV_j / JACj
ssnl3(i,ID) = 0.5 * RHV_j / JACj
ELSE
RHV_j = FT * ((WjT/WNj) * E(i) * E(i) - 2. * (WiT/WNi) * E(j) * E(i))
DIA_i = FT * ((WjT/WNj) * E(i) - 2. * (WiT/WNi) * E(j))
IF ( RHV_j .LE. 0. ) THEN
RHV_j = 0.
RHV_i = 0.
DIA_j = 0.
DIA_i = 0.
ENDIF
IMATRA(i,ID) = IMATRA(i,ID) - 0.5 * RHV_j / JACj
IMATDA(j,ID) = IMATDA(j,ID) - DIA_i
ssnl3(i,ID) = 0.5 * RHV_j / JACj
ENDIF
ENDDO
ENDDO
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
subroutine snl33 (ip, hs, smespc, acloc, imatra, imatda, ssnl3)
use datapool
implicit none
integer, intent(in) :: ip
real(rkind), intent(in) :: hs, smespc
real(rkind), intent(out) :: ssnl3(msc,mdc)
real(rkind), intent(in) :: acloc(msc,mdc)
real(rkind), intent(inout) :: imatra(msc,mdc), imatda(msc,mdc)
INTEGER :: I, J
INTEGER :: IRES, ISMAX
INTEGER :: IS, ID
REAL(rkind) :: PTTRIAD(5)
REAL(rkind) :: DEP_2, DEP_3, BIPH, RINT
REAL(rkind) :: AUX1, AUX2
REAL(rkind) :: WiT,WjT,WNi,WNj,CGi,CGj,JACi,JACj,XISTRI,Ci,Cj
REAL(rkind) :: SINBPH, FT, RHV_i, RHV_j, DIA_i, DIA_j
REAL(rkind) :: E(MSC), URSELL
LOGICAL :: TRIEXP
PTTRIAD(1) = 0.25
PTTRIAD(2) = 2.5
PTTRIAD(3) = 10.
PTTRIAD(4) = 0.2
PTTRIAD(5) = 0.01
IF (TRICO .GT. 0.) PTTRIAD(1) = TRICO
IF (TRIRA .GT. 0.) PTTRIAD(2) = TRIRA
IF (TRIURS .GT. 0.) PTTRIAD(5) = TRIURS
DEP_2 = DEP(IP)
DEP_3 = DEP(IP)
TRIEXP = .FALSE.
CALL URSELL_NUMBER(HS,SMESPC,DEP(IP),URSELL)
IF ( URSELL .lt. PTTRIAD(5) ) RETURN
E = 0.
IRES = NINT ( LOG(TWO) / XISLN )
ISMAX = 1
DO IS = 1, MSC
IF ( SPSIG(IS) .LT. ( PTTRIAD(2) * SMESPC) ) THEN
ISMAX = IS
ENDIF
ENDDO
ISMAX = MAX ( ISMAX , IRES + 1 )
BIPH = (0.5*PI)*(MyTANH(PTTRIAD(4)/URSELL)-1)
SINBPH = ABS( SIN(BIPH) )
DO ID = 1, MDC
DO IS = 1, MSC
E(IS) = ACLOC(IS,ID) * PI2 * SPSIG(IS)
ENDDO
DO I = 1, ISMAX-IRES
J = I + IRES
WiT = SPSIG(I)
WjT = SPSIG(J)
WNi = WK(I,IP)
WNj = WK(J,IP)
Ci = WiT / WNi
Cj = WjT / WNj
CGi = CG(i,IP)
CGj = CG(j,IP)
JACi = PI2 * WiT
JACj = PI2 * WjT
AUX1 = WNi**2 * ( G9 * DEP(IP) + 2. * Ci**2 )
AUX2 = WNj * DEP(IP) * ( G9 * DEP(IP) + (2./15.) * G9 * DEP_3 * WNj**2 - (2./ 5.) * WjT**2 * DEP_2 )
RINT = AUX1 / AUX2
FT = PTTRIAD(1) * Cj * CGj * RINT**2 * SINBPH
RHV_i = 0.
RHV_j = 0.
DIA_i = 0.
DIA_j = 0.
IF ( TRIEXP ) THEN
RHV_j = FT * ( E(i) * E(i) - 2. * E(j) * E(i) )
RHV_i = RHV_j
IF ( RHV_j .LE. 0. ) THEN
RHV_j = 0.
RHV_i = 0.
ENDIF
IMATRA(i,ID) = IMATRA(i,ID) - 2. * RHV_i / JACi
IMATRA(j,ID) = IMATRA(j,ID) + RHV_j / JACj
ssnl3(i,id) = - 2. * RHV_i / JACi
ELSE
RHV_j = FT * ( E(i) * E(i) - 2. * E(j) * E(i) )
DIA_i = FT * ( E(i) - 2. * E(j) )
IF ( RHV_j .LE. 0. ) THEN
RHV_j = 0.
RHV_i = 0.
DIA_j = 0.
DIA_i = 0.
ENDIF
IMATRA(j,ID) = IMATRA(j,ID) + RHV_j / JACj
IMATDA(i,ID) = IMATDA(i,ID) - 2. * DIA_i
ssnl3(j,id) = RHV_j / JACj
ENDIF
ENDDO
ENDDO
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
subroutine triad_polnikov (ip, hs, smespc, acloc, imatra, imatda, ssnl3)
use datapool
implicit none
integer, intent(in) :: ip
real(rkind), intent(in) :: hs, smespc
real(rkind), intent(out) :: ssnl3(msc,mdc)
real(rkind), intent(in) :: acloc(msc,mdc)
real(rkind), intent(inout) :: imatra(msc,mdc), imatda(msc,mdc)
INTEGER :: J, IS, ID, IT, ITER
REAL(rkind) :: KI, KJ , ETOT
REAL(rkind) :: PROPFAK, DBETA, DSIGMA
REAL(rkind) :: BETA_0(MSC), BETA_1(MSC), SNL3(MSC,MDC)
REAL(rkind) :: DELTAK(MSC), DK, H, TMP1, TMP2, TMP3
REAL(rkind) :: PART1(MSC), PART2(MSC), EPS(MSC)
REAL(rkind) :: SUMAC, KJKIKJ, KJKJKI
ITER = 10
ETOT = hs**2/FOUR
DELTAK(1) = WK(IP,1)
DO IS = 2, MSC
DELTAK(IS) = WK(IS,IP) - WK(IS-1,IP)
END DO
H = DEP(IP)
IF (ETOT .GT. THR) THEN
DO ID = 1, MDC
DO IS = 1, MSC
BETA_0(IS) = 0.1 * SPSIG(IS)/PI2
KI = WK(IS,IP)
DK = DELTAK(IS)
PROPFAK = 9. * KI * SQRT(G9) * DK / (32. * PI * SQRT(H))
DO IT = 1, ITER
BETA_1(IS) = 0.
PART1(IS) = 0.
PART2(IS) = 0.
DO J = 1, IS - 1
SUMAC = MyREAL(( ACLOC(J,ID) * CG(J,IP) - ACLOC(IS-J,ID) * CG(IS-J,IP) ))
KJ = WK(J,IP)
KJKIKJ = KJ*(KI-KJ)
KJKJKI = KJ*(KJ-KI)
DSIGMA = 0.5 * SQRT(G9*H**5.) * ABS(KI*KJKIKJ)
DBETA = BETA_0(IS) / ( PI * DSIGMA**2. + BETA_0(IS)**2. )
PART1(IS) = PART1(IS) + KJKIKJ * DBETA * SUMAC
IF (PART1(IS) .NE. PART1(IS)) WRITE (*,*) 'PART1', PART1(IS), KJKIKJ, DBETA, DSIGMA, BETA_0(IS)
END DO
DO J = IS+1, MSC
SUMAC = MyREAL(( ACLOC(J-IS,ID) * CG(J-IS,IP) - ACLOC(J,ID) * CG(J,IP) ))
KJ = WK(J,IP)
KJKIKJ = KJ*(KI-KJ)
KJKJKI = KJ*(KJ-KI)
DSIGMA = 0.5 * SQRT(G9*H**5.) * ABS(KI*KJKIKJ)
DBETA = BETA_0(IS) / ( PI * DSIGMA**2. + BETA_0(IS)**2. )
PART2(IS) = PART2(IS) + KJKJKI * DBETA * SUMAC
IF (PART2(IS) .NE. PART2(IS)) WRITE (*,*)'PART2', PART2(IS), KJKIKJ, DBETA, DSIGMA, BETA_0(IS)
END DO
BETA_1(IS) = PROPFAK * ( PART1(IS) - 2. * PART2(IS) )
IF (ABS(BETA_1(IS)) .GT. THR) THEN
EPS(IS) = ( ABS( BETA_0(IS) - BETA_1(IS) ) ) / ABS(BETA_0(IS))
ELSE
EPS(IS) = 0.
END IF
IF (ABS(EPS(IS)) .LT. 10E-3) THEN
EXIT
END IF
BETA_0(IS) = BETA_1(IS)
END DO
PART1(IS) = 0.
PART2(IS) = 0.
DO J = 1, IS - 1
IF (BETA_0(IS) .LT. THR) CYCLE
KJ = WK(J,IP)
KJKIKJ = KJ*(KI-KJ)
KJKJKI = KJ*(KJ-KI)
DSIGMA = 0.5 * SQRT(G9*H**5.) * ABS(KI*KJKIKJ)
DBETA = BETA_0(IS) / ( PI * DSIGMA**2. + BETA_0(IS)**2. )
TMP1 = MyREAL(ACLOC(J,ID)*CG(J,IP)*ACLOC(IS-J,ID)*CG(IS-J,IP))
TMP2 = MyREAL(ACLOC(J,ID)*CG(J,IP)+ACLOC(IS-J,ID)*CG(IS-J,IP))
TMP3 = (TMP1 - ACLOC(IS,ID)*CG(IS,IP)*TMP2)
PART1(IS) = PART1(IS) + KJKIKJ * DBETA * TMP3
IF (PART1(IS) .NE. PART1(IS)) WRITE (*,*) 'PART1', PART1(IS), KJKIKJ, DBETA, DSIGMA, BETA_0(IS)
END DO
DO J = IS+1, MSC
IF (BETA_0(IS) .LT. THR) CYCLE
KJ = WK(J,IP)
KJKIKJ = KJ*(KI-KJ)
KJKJKI = KJ*(KJ-KI)
DSIGMA = 0.5 * SQRT(G9*H**5.) * ABS(KI*KJKIKJ)
DBETA = BETA_0(IS) / ( PI * DSIGMA**2. + BETA_0(IS)**2. )
TMP1 = MyREAL(ACLOC(IS,ID)*CG(IS,IP)*ACLOC(J-IS,ID)*CG(J-IS,IP))
TMP2 = MyREAL(ACLOC(IS,ID)*CG(IS,IP)+ACLOC(J-IS,ID)*CG(J-IS,IP))
TMP3 = (TMP1 - ACLOC(J,ID)*CG(J,IP)*TMP2)
PART2(IS) = PART2(IS) + KJKIKJ * DBETA * TMP3
IF (PART2(IS) .NE. PART2(IS)) WRITE (*,*)'PART2', PART2(IS), KJKIKJ, DBETA, DSIGMA, BETA_0(IS)
END DO
SNL3(IS,ID) = PROPFAK * ( PART1(IS) - 2. * PART2(IS) )
END DO
END DO
DO IS = 1, MSC
DO ID = 1, MDC
IMATRA(IS,ID) = IMATRA(IS,ID) + SNL3(IS,ID)
ssnl3(is,id) = SNL3(IS,ID)
IF (ACLOC(IS,ID) .GT. 0.) THEN
IMATDA(IS,ID) = IMATDA(IS,ID) + SNL3(IS,ID)/ACLOC(IS,ID)
ELSE
IMATDA(IS,ID) = 0.
END IF
END DO
END DO
END IF
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************