#include "wwm_functions.h"
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE PARAMETER4SNL
USE DATAPOOL
IMPLICIT NONE
INTEGER :: IS
REAL(rkind) :: AUX1, FREQ
REAL(rkind) :: LAMBDA, LAMM2, LAMP2
REAL(rkind) :: DELTH3, DELTH4
INTEGER :: IDP, IDP1, IDM, IDM1
REAL(rkind) :: CIDP, WIDP, WIDP1, CIDM, WIDM, WIDM1
INTEGER :: ISP, ISP1, ISM, ISM1
REAL(rkind) :: WISP, WISP1, WISM, WISM1
REAL(rkind) :: AWG1, AWG2, AWG3, AWG4, AWG5, AWG6, AWG7, AWG8
REAL(rkind) :: SWG1, SWG2, SWG3, SWG4, SWG5, SWG6, SWG7, SWG8
INTEGER :: ISLOW, ISHGH, ISCLW, ISCHG, IDLOW, IDHGH
!
! *** set values for the nonlinear four-wave interactions ***
!
LAMBDA = PQUAD(1)
LAMM2 = (1.0-LAMBDA)**2
LAMP2 = (1.0+LAMBDA)**2
DELTH3 = MyACOS((LAMM2**2+4.0-LAMP2**2)/(4.0*LAMM2))
AUX1 = MySIN(DELTH3)
DELTH4 = MyASIN(-AUX1*LAMM2/LAMP2)
!
! *** Compute directional indices in sigma and theta space ***
!
CIDP = ABS(DELTH4/DDIR)
IDP = INT(CIDP)
IDP1 = IDP + 1
WIDP = CIDP - MyREAL(IDP)
WIDP1 = 1.0 - WIDP
CIDM = ABS(DELTH3/DDIR)
IDM = INT(CIDM)
IDM1 = IDM + 1
WIDM = CIDM - MyREAL(IDM)
WIDM1 = 1.0 - WIDM
ISP = INT(LOG(ONE+LAMBDA)/XISLN)
ISP1 = ISP + 1
WISP = (1.0+LAMBDA-XIS**ISP)/(XIS**ISP1-XIS**ISP)
WISP1 = 1.0 - WISP
ISM = INT(LOG(ONE-LAMBDA)/XISLN)
ISM1 = ISM - 1
WISM = (XIS**ISM-(1.0-LAMBDA))/(XIS**ISM-XIS**ISM1)
WISM1 = 1.0 - WISM
!
! *** Range of calculations ***
!
ISLOW = 1 + ISM1
ISHGH = MSC + ISP1 - ISM1
ISCLW = 1
ISCHG = MSC - ISM1
IDLOW = 1 - MAX(IDM1,IDP1)
IDHGH = MDC + MAX(IDM1,IDP1)
MSC4MI = ISLOW
MSC4MA = ISHGH
MDC4MI = IDLOW
MDC4MA = IDHGH
MSCMAX = MSC4MA - MSC4MI + 1
MDCMAX = MDC4MA - MDC4MI + 1
!
! *** Interpolation weights ***
!
AWG1 = WIDP *WISP
AWG2 = WIDP1*WISP
AWG3 = WIDP *WISP1
AWG4 = WIDP1*WISP1
AWG5 = WIDM *WISM
AWG6 = WIDM1*WISM
AWG7 = WIDM *WISM1
AWG8 = WIDM1*WISM1
!
! *** quadratic interpolation
!
SWG1 = AWG1**2
SWG2 = AWG2**2
SWG3 = AWG3**2
SWG4 = AWG4**2
SWG5 = AWG5**2
SWG6 = AWG6**2
SWG7 = AWG7**2
SWG8 = AWG8**2
!
! *** fill the arrays
!
WWINT(1) = IDP
WWINT(2) = IDP1
WWINT(3) = IDM
WWINT(4) = IDM1
WWINT(5) = ISP
WWINT(6) = ISP1
WWINT(7) = ISM
WWINT(8) = ISM1
WWINT(9) = ISLOW
WWINT(10) = ISHGH
WWINT(11) = ISCLW
WWINT(12) = ISCHG
WWINT(13) = IDLOW
WWINT(14) = IDHGH
WWINT(15) = MSC4MI
WWINT(16) = MSC4MA
WWINT(17) = MDC4MI
WWINT(18) = MDC4MA
WWINT(19) = MSCMAX
WWINT(20) = MDCMAX
WWAWG(1) = AWG1
WWAWG(2) = AWG2
WWAWG(3) = AWG3
WWAWG(4) = AWG4
WWAWG(5) = AWG5
WWAWG(6) = AWG6
WWAWG(7) = AWG7
WWAWG(8) = AWG8
WWSWG(1) = SWG1
WWSWG(2) = SWG2
WWSWG(3) = SWG3
WWSWG(4) = SWG4
WWSWG(5) = SWG5
WWSWG(6) = SWG6
WWSWG(7) = SWG7
WWSWG(8) = SWG8
!
! *** Fill scaling array (f**11) ***
! *** compute the radian frequency**11 for IS=ISHGH, ISLOW ***
!
IF (LTEST) THEN
WRITE(STAT%FHNDL,*) 'PARAMETER 4 SNL'
WRITE(STAT%FHNDL,*) IDP, IDP1, IDM, IDM1
WRITE(STAT%FHNDL,*) ISP, ISP1, ISM, ISM1
WRITE(STAT%FHNDL,*) ISLOW, ISHGH, ISCLW, ISCHG, IDLOW, IDHGH
WRITE(STAT%FHNDL,*) XIS
WRITE(STAT%FHNDL,*) AWG1, AWG2, AWG3, AWG4
WRITE(STAT%FHNDL,*) AWG5, AWG6, AWG7, AWG8
WRITE(STAT%FHNDL,*) '---------------------------------------'
END IF
!
IF (ALLOCATED (AF11)) DEALLOCATE (AF11)
ALLOCATE( AF11(MSC4MI:MSC4MA), stat=istat)
IF (istat/=0) CALL WWM_ABORT('wwm_snl4, allocate error 1')
DO IS = 1, MSC
AF11(IS) = (SPSIG(IS)/PI2)**11.0
END DO
FREQ = SPSIG(MSC)/PI2
DO IS = MSC+1, ISHGH
FREQ = FREQ*XIS
AF11(IS) = FREQ**11.0
END DO
FREQ = SPSIG(1)/PI2
DO IS = 0, ISLOW, -1
FREQ = FREQ/XIS
AF11(IS) = FREQ**11.0
END DO
RETURN
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE SNL4(IP, KMESPC, ACLOC, IMATRA, IMATDA)
USE DATAPOOL
IMPLICIT NONE
INTEGER :: IP
REAL(rkind) :: KMESPC
REAL(rkind), INTENT(INOUT) :: IMATRA(MSC,MDC), IMATDA(MSC,MDC)
REAL(rkind), INTENT(IN) :: ACLOC(MSC,MDC)
REAL(rkind) :: PWTAIL, FACHFR
REAL(rkind) :: LAMBDA
REAL(rkind) :: SNLC1, SNLCS1, SNLCS2, SNLCS3
REAL(rkind) :: CONS, FACTOR
INTEGER :: K3SF, K3SB, K4SF, K4SB
INTEGER :: K3D1, K3D2, K3D3, K3D4
INTEGER :: K4D1, K4D2, K4D3, K4D4
INTEGER :: ISHGH, ISCLO, ISCHG, IDLOW, IDHGH
INTEGER :: IDP, IDP1, IDM, IDM1
INTEGER :: ISP, ISP1, ISM, ISM1
INTEGER :: IS, ID, ID0
REAL(rkind) :: AWG1, AWG2, AWG3, AWG4, AWG5, AWG6, AWG7, AWG8
REAL(rkind) :: E00, EP1, EM1, EP2, EM2, SIGPI2
REAL(rkind) :: SA11, SA22, AUX, AUX2
REAL(rkind) :: SA1A, SA1B, SA2A, SA2B, SFNL(MSC,MDC)
REAL(rkind) :: UE(MSC4MI:MSC4MA, MDC4MI:MDC4MA)
IDP = WWINT(1)
IDP1 = WWINT(2)
IDM = WWINT(3)
IDM1 = WWINT(4)
ISP = WWINT(5)
ISP1 = WWINT(6)
ISM = WWINT(7)
ISM1 = WWINT(8)
ISHGH = WWINT(10)
ISCLO = WWINT(11)
ISCHG = WWINT(12)
IDLOW = WWINT(13)
IDHGH = WWINT(14)
AWG1 = WWAWG(1)
AWG2 = WWAWG(2)
AWG3 = WWAWG(3)
AWG4 = WWAWG(4)
AWG5 = WWAWG(5)
AWG6 = WWAWG(6)
AWG7 = WWAWG(7)
AWG8 = WWAWG(8)
UE(:,:) = 0.0
SFNL(:,:) = 0.0
!
! *** Calculate factor R(X) to calculate the NL wave-wave ***
! *** interaction for shallow water ***
! *** SNLC1 = CONSTANT * GRAV**-4 (CONSTANT = 3.E7) ***
!
PWTAIL = PTAIL(1)
LAMBDA = PQUAD(1)
DAL1 = 1.0/(1.0+LAMBDA)**4.0
DAL2 = 1.0/(1.0-LAMBDA)**4.0
DAL3 = 2.0*DAL1*DAL2
SNLC1 = PQUAD(2)/G9**4.0
SNLCS1 = PQUAD(3)
SNLCS2 = PQUAD(4)
SNLCS3 = PQUAD(5)
!AUX = MAX(DEP(IP)*KMESPC,1.363_rkind) ! ... check this kind of stuff with peters paper ...
AUX = MAX(0.75*DEP(IP)*KMESPC, 0.5_rkind)
AUX2 = MAX ( -1.E15_rkind, SNLCS3*AUX)
CONS = SNLC1 * ( 1. + SNLCS1/AUX * (1.-SNLCS2*AUX) * EXP(AUX2))
FACHFR = 1.0/XIS**PWTAIL
DO ID = 1, MDC
DO IS = 1, MSC
UE(IS,ID) = ACLOC(IS,ID)*SPSIG(IS)*PI2
END DO
END DO
DO ID = 1, IDHGH - MDC
ID0 = 1 - ID
DO IS = 1, MSC
UE(IS,MDC+ID) = UE(IS,ID )
UE(IS,ID0 ) = UE(IS,MDC+ID0)
END DO
END DO
DO IS = MSC + 1, ISHGH
DO ID = IDLOW, IDHGH
UE(IS,ID) = UE(IS-1,ID)*FACHFR
END DO
END DO
DO IS = ISCLO, ISCHG
DO ID = 1, MDC
K3SF = IS+ISP1
K3SB = IS+ISP
K4SF = IS+ISM1
K4SB = IS+ISM
K3D1 = ID+IDP1
K3D2 = ID+IDP
K4D1 = ID-IDM1
K4D2 = ID-IDM
K3D3 = ID-IDP1
K3D4 = ID-IDP
K4D3 = ID+IDM1
K4D4 = ID+IDM
E00 = UE(IS,ID)
EP1 = AWG1 * UE(K3SF ,K3D1 ) &
& +AWG2 * UE(K3SF ,K3D2 ) &
& +AWG3 * UE(K3SB ,K3D1 ) &
& +AWG4 * UE(K3SB ,K3D2 )
EM1 = AWG5 * UE(K4SF ,K4D1 ) &
& +AWG6 * UE(K4SF ,K4D2 ) &
& +AWG7 * UE(K4SB ,K4D1 ) &
& +AWG8 * UE(K4SB ,K4D2 )
EP2 = AWG1 * UE(K3SF ,K3D3 ) &
& +AWG2 * UE(K3SF ,K3D4 ) &
& +AWG3 * UE(K3SB ,K3D3 ) &
& +AWG4 * UE(K3SB ,K3D4 )
EM2 = AWG5 * UE(K4SF ,K4D3 ) &
& +AWG6 * UE(K4SF ,K4D4 ) &
& +AWG7 * UE(K4SB ,K4D3 ) &
& +AWG8 * UE(K4SB ,K4D4 )
FACTOR = CONS*AF11(IS)*E00
SA1A = E00*(EP1*DAL1+EM1*DAL2)
SA1B = SA1A-EP1*EM1*DAL3
SA2A = E00*(EP2*DAL1+EM2*DAL2)
SA2B = SA2A-EP2*EM2*DAL3
SA11 = FACTOR*SA1B
SA22 = FACTOR*SA2B
IF ((IS >= 1) .AND. (IS <= MSC)) THEN
SFNL(IS,ID) = SFNL(IS,ID) - 2.0 * ( SA11 + SA22 )
END IF
IF (K3D1 > MDC) K3D1 = K3D1 - MDC
IF (K3D2 > MDC) K3D2 = K3D2 - MDC
IF (K4D1 < 1 ) K4D1 = MDC + K4D1
IF (K4D2 < 1 ) K4D2 = MDC + K4D2
IF (K3D3 < 1 ) K3D3 = MDC + K3D3
IF (K3D4 < 1 ) K3D4 = MDC + K3D4
IF (K4D3 > MDC) K4D3 = K4D3 - MDC
IF (K4D4 > MDC) K4D4 = K4D4 - MDC
IF ((K3SF >= 1) .AND. (K3SF <= MSC) .AND. (K3SB >= 1) .AND. (K3SB <= MSC)) THEN
SFNL(K3SF,K3D1) = SFNL(K3SF,K3D1) + AWG1 * SA11
SFNL(K3SF,K3D2) = SFNL(K3SF,K3D2) + AWG2 * SA11
SFNL(K3SB,K3D1) = SFNL(K3SB,K3D1) + AWG3 * SA11
SFNL(K3SB,K3D2) = SFNL(K3SB,K3D2) + AWG4 * SA11
SFNL(K3SF,K3D3) = SFNL(K3SF,K3D3) + AWG1 * SA22
SFNL(K3SF,K3D4) = SFNL(K3SF,K3D4) + AWG2 * SA22
SFNL(K3SB,K3D3) = SFNL(K3SB,K3D3) + AWG3 * SA22
SFNL(K3SB,K3D4) = SFNL(K3SB,K3D4) + AWG4 * SA22
END IF
IF ((K4SF >= 1) .AND. (K4SF <= MSC) .AND. (K4SB >= 1) .AND. (K4SB <= MSC)) THEN
SFNL(K4SF,K4D1) = SFNL(K4SF,K4D1) + AWG5 * SA11
SFNL(K4SF,K4D2) = SFNL(K4SF,K4D2) + AWG6 * SA11
SFNL(K4SB,K4D1) = SFNL(K4SB,K4D1) + AWG7 * SA11
SFNL(K4SB,K4D2) = SFNL(K4SB,K4D2) + AWG8 * SA11
SFNL(K4SF,K4D3) = SFNL(K4SF,K4D3) + AWG5 * SA22
SFNL(K4SF,K4D4) = SFNL(K4SF,K4D4) + AWG6 * SA22
SFNL(K4SB,K4D3) = SFNL(K4SB,K4D3) + AWG7 * SA22
SFNL(K4SB,K4D4) = SFNL(K4SB,K4D4) + AWG8 * SA22
END IF
END DO
END DO
!
! *** SFNL -> energy density, IMATRA -> action density
!
DO IS = 1, MSC
SIGPI2 = SPSIG(IS)*PI2
DO ID = 1, MDC
IF (ICOMP >= 2) THEN
IF (SFNL(IS,ID).GT.0.) THEN
IMATRA(IS,ID) = IMATRA(IS,ID) + SFNL(IS,ID) / SIGPI2
ELSE
IF (ACLOC(IS,ID) .GT. 10.E-10) IMATDA(IS,ID) = IMATDA(IS,ID) - SFNL(IS,ID) / ACLOC(IS,ID)*SIGPI2
END IF
ELSE
IMATRA(IS,ID) = IMATRA(IS,ID) + SFNL(IS,ID) / SIGPI2
IF (ACLOC(IS,ID) .GT. 10.E-10) IMATDA(IS,ID) = IMATDA(IS,ID) + SFNL(IS,ID) / ACLOC(IS,ID)*SIGPI2
END IF
END DO
END DO
RETURN
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************