#include "wwm_functions.h"
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE ST6_PRE (IP, WALOC, PHI, DPHIDN, SSINE, DSSINE, SSDS, DSSDS, SSNL4, DSSNL4, SSINL)
USE DATAPOOL
USE W3SRC6MD
IMPLICIT NONE
INTEGER, INTENT(IN) :: IP
REAL(rkind), INTENT(IN) :: WALOC(NUMSIG,NUMDIR)
REAL(rkind), INTENT(OUT) :: PHI(NUMSIG,NUMDIR), DPHIDN(NUMSIG,NUMDIR)
REAL(rkind), INTENT(OUT) :: SSINE(NUMSIG,NUMDIR), DSSINE(NUMSIG,NUMDIR)
REAL(rkind), INTENT(OUT) :: SSDS(NUMSIG,NUMDIR),DSSDS(NUMSIG,NUMDIR)
REAL(rkind), INTENT(OUT) :: SSNL4(NUMSIG,NUMDIR),DSSNL4(NUMSIG,NUMDIR)
REAL(rkind), INTENT(OUT) :: SSINL(NUMSIG,NUMDIR)
INTEGER :: IS, ID, ITH, IK, IS0
REAL(rkind) :: AWW3(NSPEC)
REAL(rkind) :: VDDS(NSPEC), VSDS(NSPEC), BRLAMBDA(NSPEC)
REAL(rkind) :: WN2(NUMSIG*NUMDIR), WHITECAP(1:4)
REAL(rkind) :: VSIN(NSPEC), VDIN(NSPEC)
REAL(rkind) :: ETOT, FAVG, FMEAN1, AS, SUMWALOC, FAVGWS
REAL(rkind) :: TAUWAX, TAUWAY, AMAX, FPM, WIND10, WINDTH, FP
REAL(rkind) :: HS,SME01,SME10,SMECP,KME01,KMWAM,KMWAM2,WNMEAN
DO IS = 1, NUMSIG
DO ID = 1, NUMDIR
AWW3(ID + (IS-1) * NUMDIR) = WALOC(IS,ID) * CG(IS,IP)
END DO
END DO
DO IK=1, NK
WN2(1+(IK-1)*NTH) = WK(IK,IP)
END DO
DO IK=1, NK
IS0 = (IK-1)*NTH
DO ITH=2, NTH
WN2(ITH+IS0) = WN2(1+IS0)
END DO
END DO
!
! wind input
!
TAUWX(IP) = ZERO
TAUWY(IP) = ZERO
SSINL = ZERO
NUMSIG_HF(IP) = NUMSIG
AS = 0.
BRLAMBDA = ZERO
IF (MESIN .GT. 0) THEN
CALL SET_WIND( IP, WIND10, WINDTH )
CALL SET_FRICTION( IP, WALOC, WIND10, WINDTH, FPM )
#ifdef DEBUGSRC
WRITE(740+myrank,*) '1: input value USTAR=', UFRIC(IP), ' USTDIR=', USTDIR(IP)
#endif
CALL W3SPR6(AWW3, CG(:,IP), WK(:,IP), EMEAN(IP), FMEAN(IP), WNMEAN, AMAX, FP)
#ifdef DEBUGSRC
WRITE(740+myrank,*) '1: out value USTAR=', UFRIC(IP), ' USTDIR=', USTDIR(IP)
WRITE(740+myrank,*) '1: out value EMEAN=', EMEAN(IP), ' FMEAN=', FMEAN(IP)
WRITE(740+myrank,*) '1: out value FMEAN1=', FMEAN1, ' WNMEAN=', WNMEAN
WRITE(740+myrank,*) '1: out value CD=', CD(IP), ' Z0=', Z0(IP)
WRITE(740+myrank,*) '1: out value ALPHA=', ALPHA_CH(IP), ' FMEANWS=', FMEANWS(IP)
#endif
IF (EMEAN(IP) .LT. THR .AND. WIND10 .GT. THR) CALL SIN_LIN_CAV(IP,WIND10,WINDTH,FPM,SSINL)
CALL W3SIN6 (AWW3, CG(:,IP), WN2, UFRIC(IP), WINDTH, CD(IP), TAUWX(IP), TAUWY(IP), TAUWAX, TAUWAY, VSIN, VDIN )
#ifdef DEBUGSRC
WRITE(740+myrank,*) '1: WINDTH=', WINDTH, ' Z0=', Z0(IP), ' CD=', CD(IP)
WRITE(740+myrank,*) '1: UFRIC=', UFRIC(IP), 'WIND10=', WIND10, ' RHOAW=', RHOAW
WRITE(740+myrank,*) '1: TAUWX=', TAUWX(IP), ' TAUWY=', TAUWY(IP)
WRITE(740+myrank,*) '1: TAUWAX=', TAUWAX, ' TAUWAY=', TAUWAY
WRITE(740+myrank,*) '1: W3SIN4min/max/sum(VSIN)=', minval(VSIN), maxval(VSIN), sum(VSIN)
WRITE(740+myrank,*) '1: W3SIN4min/max/sum(VDIN)=', minval(VDIN), maxval(VDIN), sum(VDIN)
#endif
CALL W3SPR6(AWW3, CG(:,IP), WK(:,IP), EMEAN(IP), FMEAN(IP), WNMEAN, AMAX, FP)
CALL W3SIN6 (AWW3, CG(:,IP), WN2, UFRIC(IP), WINDTH, CD(IP), TAUWX(IP), TAUWY(IP), TAUWAX, TAUWAY, VSIN, VDIN )
#ifdef DEBUGSRC
WRITE(740+myrank,*) '2: W3SIN4min/max/sum(VSIN)=', minval(VSIN), maxval(VSIN), sum(VSIN)
WRITE(740+myrank,*) '2: W3SIN4min/max/sum(VDIN)=', minval(VDIN), maxval(VDIN), sum(VDIN)
#endif
CALL CONVERT_VS_VD_WWM(IP, VSIN, VDIN, SSINE, DSSINE)
ENDIF
IF (MESNL .GT. 0) THEN
CALL MEAN_WAVE_PARAMETER(IP,WALOC,HS,ETOT,SME01,SME10,SMECP,KME01,KMWAM,KMWAM2)
CALL DIASNL4WW3(IP, KMWAM, WALOC, SSNL4, DSSNL4)
END IF
IF (MESDS .GT. 0) THEN
CALL W3SDS6 (AWW3, CG(:,IP), WK(:,IP), VSDS, VDDS)
#ifdef DEBUGSRC
WRITE(740+myrank,*) '2: W3SDS4min/max/sum(VSDS)=', minval(VSDS), maxval(VSDS), sum(VSDS)
WRITE(740+myrank,*) '2: W3SDS4min/max/sum(VDDS)=', minval(VDDS), maxval(VDDS), sum(VDDS)
#endif
CALL CONVERT_VS_VD_WWM(IP, VSDS, VDDS, SSDS, DSSDS)
ENDIF
!
PHI = SSINL + SSINE + SSNL4 + SSDS
DPHIDN = DSSINE + DSSNL4 + DSSDS
!
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************
SUBROUTINE ST6_POST (IP, WALOCOLD, WALOC)
!Now here quite some memory and indexing overhead ... will consolidate this after validation study ... now it fits our work in WW3
USE DATAPOOL
USE W3SRC6MD
IMPLICIT NONE
INTEGER, INTENT(IN) :: IP
REAL(rkind), INTENT(IN) :: WALOCOLD(NUMSIG,NUMDIR)
REAL(rkind), INTENT(OUT) :: WALOC(NUMSIG,NUMDIR)
REAL(rkind) :: SSINE(NUMSIG,NUMDIR),DSSINE(NUMSIG,NUMDIR)
REAL(rkind) :: SSDS(NUMSIG,NUMDIR),DSSDS(NUMSIG,NUMDIR)
REAL(rkind) :: SSINL(NUMSIG,NUMDIR)
INTEGER :: IS, ID, IK, ITH, ITH2, IS0, NKH, NKH1
REAL(rkind) :: PHI(NUMSIG,NUMDIR), DPHIDN(NUMSIG,NUMDIR)
REAL(rkind) :: AWW3(NSPEC), AWW3OLD(NSPEC), WN2(NUMSIG*NUMDIR), BRLAMBDA(NSPEC)
REAL(rkind) :: SSINE_WW3(NSPEC), DSSINE_WW3(NSPEC), TMP_DS(NUMSIG)
REAL(rkind) :: SSDS_WW3(NSPEC), DSSDS_WW3(NSPEC), SSNL4_WW3(NSPEC), DSSNL4_WW3(NSPEC)
REAL(rkind) :: SSBF_WW3(NSPEC), DSSBF_WW3(NSPEC), SSNL4(NUMSIG,NUMDIR), DSSNL4(NUMSIG,NUMDIR)
REAL(rkind) :: SSBF(NUMSIG,NUMDIR), DSSBF(NUMSIG,NUMDIR)
REAL(rkind) :: ETOT, FAVG, FMEAN1, WNMEAN, AS, FAVGWS
REAL(rkind) :: TAUWAX, TAUWAY, AMAX, WIND10, WINDTH, EBAND, B1BAND
REAL(rkind) :: WHITECAP(1:4), SUMWALOC, FPM, FH1, FH2, TAUBBL(2)
REAL(rkind) :: PHIAW, CHARN, PHINL, PHIBBL, TAUWIX, TAUWIY, TAUWNX, TAUWNY, TAUOX, TAUOY
REAL(rkind) :: FACTOR, FACTOR2, MWXFINISH, MWYFINISH, EFINISH, DIFF, CONST2, TEMP2
REAL(rkind) :: A1BAND, PHIOC, HSTOT, PIBBL, FAGE, FHIGH, FACHFA, FP
DO IS = 1, NUMSIG
DO ID = 1, NUMDIR
AWW3(ID + (IS-1) * NUMDIR) = WALOC(IS,ID) * CG(IS,IP)
AWW3OLD(ID + (IS-1) * NUMDIR) = WALOCOLD(IS,ID) * CG(IS,IP)
END DO
END DO
DO IK=1, NK
WN2(1+(IK-1)*NTH) = WK(IK,IP)
END DO
DO IK=1, NK
IS0 = (IK-1)*NTH
DO ITH=2, NTH
WN2(ITH+IS0) = WN2(1+IS0)
END DO
END DO
!
! wind input
!
AS = ZERO
NUMSIG_HF(IP) = NUMSIG
CALL SET_WIND( IP, WIND10, WINDTH )
CALL SET_FRICTION( IP, WALOC, WIND10, WINDTH, FPM )
CALL W3SPR6(AWW3, CG(:,IP), WK(:,IP), EMEAN(IP), FMEAN(IP), WNMEAN, AMAX, FP)
IF (EMEAN(IP) .LT. THR .AND. WIND10 .GT. THR) CALL SIN_LIN_CAV(IP,WIND10, WINDTH,FPM,SSINL)
CALL W3SIN6 (AWW3, CG(:,IP), WN2, UFRIC(IP), WINDTH, CD(IP), TAUWX(IP), TAUWY(IP), TAUWAX, TAUWAY, SSINE_WW3, DSSINE_WW3 )
!
! dissipation
!
CALL W3SDS6 (AWW3, CG(:,IP), WK(:,IP),SSDS_WW3,DSSDS_WW3)
!
! snl4
!
IF (MESNL .GT. 0) CALL DIASNL4WW3(IP, WNMEAN, WALOC, SSNL4, DSSNL4)
!
! and last but not least the bottom friction dissipation ...
!
IF (MESBF .GT. 0) CALL SDS_BOTF(IP,WALOC,SSBF,DSSBF)
!
! Resio & Perrie ... scaling ...
!
FAGE = FFXFA*TANH(0.3*WIND10*FMEANWS(IP)*PI2/G9)
FH1 = (FFXFM+FAGE) * FMEAN(IP)
FH2 = FFXPM / UFRIC(IP)
FHIGH = MIN ( SIG(NK) , MAX ( FH1 , FH2 ) )
NKH = MIN ( NK , INT(FACTI2+FACTI1*LOG(MAX(1.E-7,FHIGH))) )
NKH1 = MIN ( NK , NKH+1 )
NKH = MAX ( 2 , MIN ( NKH1 , INT ( FACTI2 + FACTI1*LOG(MAX(1.E-7_rkind,FHIGH)) ) ) )
!
! Add tail
!
FACHF = 5.
FACHFA = XFR**(-FACHF-2)
DO IK=NKH+1, NK
DO ITH=1, NTH
AWW3(ITH+(IK-1)*NTH) = AWW3(ITH+(IK-2)*NTH) * FACHFA + 0. ! Adding a magic zero without a comment ...
END DO
END DO
!
! convert to wwm convetion ...
!
DO IS = 1, NUMSIG
DO ID = 1, NUMDIR
WALOC(IS,ID) = AWW3(ID + (IS-1) * NUMDIR) / CG(IS,IP)
END DO
END DO
CALL CONVERT_VS_VD_WWM(IP, SSINE_WW3, DSSINE_WW3, SSINE, DSSINE)
CALL CONVERT_VS_VD_WWM(IP, SSDS_WW3, DSSDS_WW3, SSDS, DSSDS)
IF (MESNL .GT. 0) CALL CONVERT_WWM_VS_VD(IP, SSNL4_WW3, DSSNL4_WW3, SSNL4, DSSNL4)
IF (MESBF .GT. 0) CALL CONVERT_WWM_VS_VD(IP, SSBF_WW3, DSSBF_WW3, SSBF, DSSBF)
!
! compute momentum to BBL
!
DO IK=1, NK
CONST2=DDEN(IK)/CG(IK,IP)*G9/(SIG(IK)/WK(IK,IP))
DO ITH=1,NTH
IS=ITH+(IK-1)*NTH
TEMP2=CONST2*DSSBF_WW3(IS)*AWW3(IS)
TAUBBL(1) = TAUBBL(1) - TEMP2*ECOS(IS)
TAUBBL(2) = TAUBBL(2) - TEMP2*ESIN(IS)
END DO
END DO
!
! adding the fluxes from waves to ocean ...
!
WHITECAP(3) = ZERO
PHIAW = ZERO
CHARN = ZERO
PHINL = ZERO
PHIBBL = ZERO
TAUWIX = ZERO
TAUWIY = ZERO
TAUWNX = ZERO
TAUWNY = ZERO
HSTOT = ZERO
DO IK = 1, NK
FACTOR = DDEN(IK)/CG(IK,IP)
FACTOR2 = FACTOR*G9*WK(IK,IP)/SIG(IK)
DO ITH = 1, NTH
IS = (IK-1) * NTH + ITH
PHIAW = PHIAW + SSINE_WW3(IS) * DT4S * FACTOR / MAX ( ONE , (ONE-DT4S*DSSINE_WW3(IS)))
PIBBL = PHIBBL - SSBF_WW3(IS) * DT4S * FACTOR / MAX ( ONE , (ONE-DT4S*DSSBF_WW3(IS)))
PHINL = PHINL + SSNL4_WW3(IS) * DT4S * FACTOR / MAX ( ONE , (ONE-DT4S*DSSNL4_WW3(IS)))
IF (SSINE_WW3(IS) .GT. ZERO) THEN
WHITECAP(3) = WHITECAP(3) + AWW3(IS) * FACTOR
ELSE
! computes the upward energy flux (counted > 0 upward)
CHARN = CHARN - (SSINE_WW3(IS))* DT4S * FACTOR / MAX ( ONE , (ONE-DT4S*DSSINE_WW3(IS)))
END IF
HSTOT = HSTOT + AWW3(IS) * FACTOR
END DO
END DO
WHITECAP(3) = 4.*SQRT(WHITECAP(3))
HSTOT = 4.*SQRT(HSTOT)
TAUWIX = TAUWIX + TAUWX(IP) * RHOAW * DT4S
TAUWIY = TAUWIY + TAUWY(IP) * RHOAW * DT4S
TAUWNX = TAUWNX + TAUWAX * RHOAW * DT4S
TAUWNY = TAUWNY + TAUWAY * RHOAW * DT4S
!
! The wave to ocean flux is the difference between initial energy
! and final energy, plus wind input plus the SNL flux to high freq.,
! minus the energy lost to the bottom boundary layer (BBL)
!
EFINISH = 0.
MWXFINISH = 0.
MWYFINISH = 0.
DO IK = 1, NK
EBAND = 0.
A1BAND = 0.
B1BAND = 0.
DO ITH = 1, NTH
DIFF = AWW3OLD(ITH+(IK-1)*NTH)-AWW3(ITH+(IK-1)*NTH) ! Check that the difference is taken in the right direction!
EBAND = EBAND + DIFF
A1BAND = A1BAND + DIFF*ECOS(ITH)
B1BAND = B1BAND + DIFF*ESIN(ITH)
END DO
EFINISH = EFINISH + EBAND * DDEN(IK) / CG(IK,IP)
MWXFINISH = MWXFINISH + A1BAND * DDEN(IK) / CG(IK,IP) * WK(IK,IP)/SIG(IK)
MWYFINISH = MWYFINISH + B1BAND * DDEN(IK) / CG(IK,IP) * WK(IK,IP)/SIG(IK)
END DO
!
! Transfoe rmation in momentum flux in m^2 / s^2
!
TAUOX = (G9*MWXFINISH+TAUWIX-TAUBBL(1))/DT4S
TAUOY = (G9*MWYFINISH+TAUWIY-TAUBBL(2))/DT4S
TAUWIX = TAUWIX/DT4S
TAUWIY = TAUWIY/DT4S
TAUWNX = TAUWNX/DT4S
TAUWNY = TAUWNY/DT4S
!
! Transformation in wave energy flux in W/m^2=kg / s^3
!
PHIOC = RHOW*G9*(EFINISH+PHIAW-PHIBBL)/DT4S
PHIAW = RHOW*G9*PHIAW /DT4S
PHINL = RHOW*G9*PHINL /DT4S
PHIBBL= RHOW*G9*PHIBBL/DT4S
#ifdef SCHISM
OUTT_INTPAR(IP,32) = TAUOX ! x-component of the wave-ocean momentum flux (tauox in m2.s-2)
OUTT_INTPAR(IP,33) = TAUOX ! y-component of the wave-ocean momentum flux (tauoy in m2.s-2)
OUTT_INTPAR(IP,34) = PHIOC ! Wave-to-ocean TKE flux (phioc in W.m-2)
OUTT_INTPAR(IP,35) = PHIBBL ! phibbl / wave_phibbl : Energy flux due to bottom friction (phioc in W.m-2)
#endif
END SUBROUTINE
!**********************************************************************
!* *
!**********************************************************************