#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,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,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
!**********************************************************************
!*                                                                    *
!**********************************************************************