#include "wwm_functions.h"
!**********************************************************************
!*                                                                    *
!**********************************************************************
!2do add mean quantities for 
      SUBROUTINE SOURCETERMS (IP, ACLOC, IMATRA, IMATDA, LRECALC, ISELECT, CALLFROM)
      
      USE DATAPOOL, ONLY : MSC, MDC, MNP, WK, LINID, THR, UFRIC, RADFLAG
      USE DATAPOOL, ONLY : CD, TAUTOT, TAUWX, TAUWY, DEP, AC1, AC2
      USE DATAPOOL, ONLY : PI2, CG, G9, ZERO, ALPHA_CH, QBLOCAL
      USE DATAPOOL, ONLY : USTDIR, Z0, SMALL, VERYSMALL, MSC_HF
      USE DATAPOOL, ONLY : DDIR, SPSIG, SPDIR, FRINTF, LMAXETOT, LADVTEST
      USE DATAPOOL, ONLY : USTDIR, Z0, SMALL, VERYSMALL, MSC_HF, DDIR
      USE DATAPOOL, ONLY : SPSIG, SPDIR, FRINTF, ONE, RHOA, RHOAW, TAUHF
      USE DATAPOOL, ONLY : TAUW, FR, MESNL, MESIN, MESDS, MESBF, MESBR, IBREAK
      USE DATAPOOL, ONLY : MESTR, ISHALLOW, DS_INCR, IOBP, IOBPD, MEVEG
      USE DATAPOOL, ONLY : LNANINFCHK, DBG, WRITEDBGFLAG, IFRIC, RTIME, DISSIPATION
      USE DATAPOOL, ONLY : AIRMOMENTUM, ONEHALF, NSPEC, RKIND, ICOMP, IROLLER
      USE DATAPOOL, ONLY : WAVE_SBRTOT, WAVE_SBFTOT, WAVE_SINTOT, WAVE_SDSTOT
      USE DATAPOOL, ONLY : EPS_W, EPS_BR
      USE DATAPOOL, ONLY : RHOW
#ifdef WWM_MPI
      USE DATAPOOL, ONLY : myrank
#endif
         USE SdsBabanin
#ifdef SNL4_TSA
         USE W3SNLXMD
#endif

#ifdef ST41
         USE W3SRC4MD_OLD
#elif ST42
         USE W3SRC4MD
#endif
         IMPLICIT NONE

         INTEGER, INTENT(IN) :: IP, ISELECT

         REAL(rkind), INTENT(OUT) :: IMATRA(MSC,MDC), IMATDA(MSC,MDC)
         REAL(rkind), INTENT(IN)  :: ACLOC(MSC,MDC)

         LOGICAL, INTENT(IN) :: LRECALC
         CHARACTER(LEN=*), INTENT(IN) :: CALLFROM

         INTEGER        :: IS, ID, IS0, IK, ITH, IDISP, JU, NZZ
         REAL(rkind)    :: WIND10, WINDTH
         REAL(rkind)    :: FPM
         REAL(rkind)    :: SME01, SME10, SMECP, KME01, KMWAM, KMWAM2
         REAL(rkind)    :: SME01WS, SME10WS
         REAL(rkind)    :: HS, ETOT, FPMH,FPM4 
         REAL(rkind)    :: LPOW, MPOW, a1, a2, ETOTWS
         REAL(rkind)    :: XRR, XPP, XFLT, XREL, FACP, XFILT
         REAL(rkind)    :: TEMP2(MSC), SSBRL(MSC,MDC)
         REAL(rkind)    :: AWW3(NSPEC), IMATRAWW3(MSC,MDC), IMATDAWW3(MSC,MDC)
         REAL(rkind)    :: EDENS(MSC), KDS(MSC), ABAB(MSC)
         REAL(rkind)    :: WHITECAP(1:4),AKMEAN,XKMEAN,F1MEAN,TMPAC(MDC,MSC),TEMP(MSC), FCONST(MSC)
         REAL(rkind)    :: ACLOC1(MSC,MDC)


         REAL(rkind)    :: SSNL3(MSC,MDC), SSNL4(MSC,MDC), SSINL(MSC,MDC), SSDS(MSC,MDC), DSSNL4(MSC,MDC)
         REAL(rkind)    :: SSVEG(MSC,MDC),DSSVEG(MSC,MDC), DSSBF(MSC,MDC)
         REAL(rkind)    :: SSBF(MSC,MDC), SSBR(MSC,MDC), SSINE(MSC,MDC), DSSNL3(MSC,MDC), DSSBR(MSC,MDC)
         REAL(rkind)    :: TMP_IN(MSC), TMP_DS(MSC), WN2(MSC*MDC),SPRDD(MDC),AK2VGM1,AKM1
         REAL(rkind)    :: SINTOT,SDSTOT,SBRTOT,SBFTOT

         REAL(rkind)    :: EMEAN, FMEAN, FMEAN1, WNMEAN, AMAX, AS
         REAL(rkind)    :: FMEANWS, TAUWAX, TAUWAY, XJ, FLLOWEST, GADIAG
         REAL(rkind)    :: IMATDA1D(NSPEC), IMATRA1D(NSPEC), SUMACLOC, IMATRAT(MSC,MDC), BRLAMBDA(NSPEC)
         REAL(rkind)    :: IMATRA_WAM(MDC,MSC), IMATDA_WAM(MDC,MSC), TAILFACTOR, FLOGSPRDM1, SNL3(MSC,MDC), DSNL3(MSC,MDC)

#ifdef SNL4_TSA
         REAL           :: IMATRA_TSA(MDC,MSC), IMATDA_TSA(MDC,MSC), TMPAC_TSA(MDC,MSC), CG_TSA(MSC), WK_TSA(MSC), DEP_TSA
#endif

         REAL           :: XNL(MSC,MDC), DDIAG(MSC,MDC), ACLOC_WRT(MSC,MDC), DEP_WRT, SPSIG_WRT(MSC), SPDIR_WRT(MDC)
         INTEGER        :: IERR_WRT

#ifdef TIMINGS 
         REAL(rkind)        :: T1, T2
         REAL(rkind), SAVE  :: TIME1, TIME2, TIME3, TIME4, TIME5, TIME6, TIME7, TIME8, TIME9
#endif

         LOGICAL        :: LWINDSEA(MSC,MDC)
         REAL(rkind)    :: XLCKS(MDC,MSC)

#ifdef TIMINGS
         call WAV_MY_WTIME(TIME1)
#endif 

!         IF (LMAXETOT .AND. .NOT. LADVTEST .AND. ISHALLOW(IP) .EQ. 1 .AND. .NOT. LRECALC) THEN
!           CALL BREAK_LIMIT(IP,ACLOC,SSBRL) ! Miche to reduce stiffness of source terms ...
!         END IF

         IDISP = 999

         IF (.NOT. LRECALC) THEN
           ACLOC1 = AC1(:,:,IP)
           CALL MEAN_WAVE_PARAMETER(IP,ACLOC1,HS,ETOT,SME01,SME10,SMECP,KME01,KMWAM,KMWAM2) ! 1st guess ... 
         END IF

         SUMACLOC    = SUM(ACLOC)
!
         WIND10      = ZERO
         BRLAMBDA    = ZERO
         SSINE       = zero
         SSINL       = zero
         SSNL3       = zero
         SSNL4       = zero
         DSSNL4      = zero
         SSBR        = zero
         SSBF        = zero
         SSDS        = zero
         DSSBF       = zero
         SSVEG       = zero
         DSSVEG      = zero
         IMATRA_WAM  = zero
         IMATDA_WAM  = zero
         TMPAC       = zero
         IMATRA      = zero 
         IMATDA      = zero 
         IMATRAWW3   = zero 
         IMATDAWW3   = zero 
         QBLOCAL(IP) = zero 

#ifdef ST_DEF
         IF (MESDS == 1 .OR. MESIN .EQ. 1) THEN
           DO IS = 1, MSC
             DO ID = 1, MDC
               AWW3(ID + (IS-1) * MDC) = ACLOC(IS,ID) * CG(IS,IP)
             END DO
           END DO
         END IF
         XPP     = 0.15
         XRR     = 0.10
         XFILT  = 0.05
         XPP     = MAX ( 1.E-6_rkind , XPP )
         XRR     = MAX ( 1.E-6_rkind , XRR )
         XREL   = XRR
         XFILT  = MAX ( ZERO , XFILT )
         XFLT   = XFILT
         FACP   = 2*XPP / PI2 * 0.62E-3 * PI2**4 / G9**2
         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
#endif

#ifdef DEBUG
         WRITE(12,*) '0', SUM(IMATRA), SUM(IMATDA), ISELECT, MESIN, CALLFROM
#endif

#ifdef TIMINGS
         call WAV_MY_WTIME(TIME2)
#endif 
         IF ((ISELECT .EQ. 1 .OR. ISELECT .EQ. 10 .OR. ISELECT .EQ. 20) .AND. .NOT. LRECALC) THEN
           IF (IOBP(IP) .EQ. 0) THEN
             IF (MESIN == 1) THEN ! Ardhuin et al. 2010
               CALL SET_WIND( IP, WIND10, WINDTH )
               TAUWX(IP) = ZERO
               TAUWY(IP) = ZERO               
               IF (SUMACLOC .LT. THR .AND. WIND10 .GT. THR) THEN
                 CALL SET_FRICTION( IP, ACLOC, WIND10, WINDTH, FPM )
                 CALL SIN_LIN_CAV(IP,WIND10,WINDTH,FPM,IMATRA,SSINL)
                 !WRITE(*,'(I10,4F20.10)') IP, WIND10, FPM, SUM(IMATRA)
               ELSE
                 MSC_HF(IP) = MSC
                 AS      = 0. 
#ifdef ST41
                 CALL W3SPR4_OLD ( AWW3, CG(:,IP), WK(:,IP), EMEAN, FMEAN, WNMEAN, AMAX, WIND10, WINDTH, UFRIC(IP), USTDIR(IP), TAUWX(IP), TAUWY(IP), CD(IP), Z0(IP), ALPHA_CH(IP), LLWS, FMEANWS)
                 CALL W3SIN4_OLD ( AWW3, CG(:,IP), WK(:,IP), WIND10, UFRIC(IP), RHOAW, AS, WINDTH, Z0(IP), CD(IP), TAUWX(IP), TAUWY(IP), TAUWAX, TAUWAY, IMATRA1D, IMATDA1D, LLWS)
                 CALL W3SPR4_OLD ( AWW3, CG(:,IP), WK(:,IP), EMEAN, FMEAN, WNMEAN, AMAX, WIND10, WINDTH, UFRIC(IP), USTDIR(IP), TAUWX(IP), TAUWY(IP), CD(IP), Z0(IP), ALPHA_CH(IP), LLWS, FMEANWS)
                 CALL W3SIN4_OLD ( AWW3, CG(:,IP), WK(:,IP), WIND10, UFRIC(IP), RHOAW, AS, WINDTH, Z0(IP), CD(IP), TAUWX(IP), TAUWY(IP), TAUWAX, TAUWAY, IMATRA1D, IMATDA1D, LLWS)
#elif ST42
                 CALL W3SPR4 ( AWW3, CG(:,IP), WK(:,IP), EMEAN, FMEAN, FMEAN1, WNMEAN, AMAX, WIND10, WINDTH, UFRIC(IP), USTDIR(IP), TAUWX(IP), TAUWY(IP), CD(IP), Z0(IP), ALPHA_CH(IP), LLWS, FMEANWS)
                 CALL W3SIN4 ( IP, AWW3, CG(:,IP), WN2,  WIND10, UFRIC(IP), RHOAW, AS, WINDTH, Z0(IP), CD(IP), TAUWX(IP), TAUWY(IP), TAUWAX, TAUWAY, IMATRA1D, IMATDA1D, LLWS, BRLAMBDA)
                 CALL W3SPR4 ( AWW3, CG(:,IP), WK(:,IP), EMEAN, FMEAN, FMEAN1, WNMEAN, AMAX, WIND10, WINDTH, UFRIC(IP), USTDIR(IP), TAUWX(IP), TAUWY(IP), CD(IP), Z0(IP), ALPHA_CH(IP), LLWS, FMEANWS)  
                 CALL W3SIN4 ( IP, AWW3, CG(:,IP), WN2, WIND10, UFRIC(IP), RHOAW, AS, WINDTH, Z0(IP), CD(IP), TAUWX(IP), TAUWY(IP), TAUWAX, TAUWAY, IMATRA1D, IMATDA1D, LLWS, BRLAMBDA) 
#else
                 IF (WRITEDBGFLAG == 1) WRITE(DBG%FHNDL,*) 'NO ST42 or ST41 chosen but MESIN == 1'
                 CALL WWM_ABORT('stop wwm_sourceterms l.176')
#endif
                 CALL ONED2TWOD(IMATRA1D,IMATRAWW3)
                 SSINE = IMATRAWW3
                 CALL ONED2TWOD(IMATDA1D,IMATDAWW3)
                 DO ID = 1, MDC
                   IMATRA(:,ID) = IMATRAWW3(:,ID) / CG(:,IP) + IMATRA(:,ID)
                   IMATDA(:,ID) = ZERO!IMATDAWW3(:,ID) !/ CG(:,IP) 
                 END DO
               END IF
               IF (LNANINFCHK) THEN
                 IF (SUM(IMATRA) .NE. SUM(IMATRA) .OR. SUM(IMATDA) .NE. SUM(IMATDA)) THEN
                   IF (WRITEDBGFLAG == 1) WRITE(DBG%FHNDL,*) 'NAN AT GRIDPOINT NORMAL', IP, 'DUE TO SIN', SUM(IMATRA), SUM(IMATDA)
                   CALL WWM_ABORT('NAN AT wwm_sourceterms.F90 l.189')
                 END IF
               ENDIF
             ELSE IF (MESIN == 2) THEN ! Cycle 4, Bidlot et al. ...
               CALL WWM_ABORT('PLEASE USE "LSOURCEWAM = T" FOR ECWAM SOURCE TERM FORMULATION') 
             ELSE IF (MESIN == 3) THEN ! Makin & Stam
               CALL SET_WIND( IP, WIND10, WINDTH )
               IFRIC = 4
               CALL SET_FRICTION( IP, ACLOC, WIND10, WINDTH, FPM )
               CALL SIN_LIN_CAV( IP, WIND10, WINDTH, FPM, IMATRA, SSINL)
               CALL SIN_MAKIN( IP, WIND10, WINDTH, KME01,ETOT,ACLOC,IMATRA,IMATDA,SSINE)
             ELSE IF (MESIN == 4) THEN ! Donealan et al.
               CALL SET_WIND( IP, WIND10, WINDTH )
               IFRIC = 1
               CALL SET_FRICTION( IP, ACLOC, WIND10, WINDTH, FPM )
               IF (.NOT. LINID) CALL SIN_LIN_CAV( IP, WIND10, WINDTH,FPM,IMATRA,SSINL)
               CALL SWIND_DBYB (IP,WIND10,WINDTH,IMATRA,SSINE)
             ELSE IF (MESIN == 5) THEN ! Cycle 3
               CALL SET_WIND( IP, WIND10, WINDTH )
               IFRIC = 1
               CALL SET_FRICTION( IP, ACLOC, WIND10, WINDTH, FPM )
               IF (.NOT. LINID) CALL SIN_LIN_CAV( IP, WIND10, WINDTH,FPM,IMATRA,SSINL)
               CALL SIN_EXP_KOMEN( IP, WINDTH, ACLOC, IMATRA, IMATDA, SSINE )
             END IF ! MESIN
           END IF ! IOBP
           !IF (IDISP == IP) THEN
             !CALL PLOT_SHADED_CONTOUR_POLAR(SPSIG/PI2,SPDIR*RADDEG,MSC,MDC,IMATRA/MAXVAL(IMATRA),50,MSC,MDC,'BEFORE ANY CALL')
             !PAUSE
           !END IF
         END IF ! ISELECT 

#ifdef DEBUG
         WRITE(12,*) '1', SUM(IMATRA), SUM(IMATDA)
#endif

#ifdef TIMINGS
         call WAV_MY_WTIME(TIME3)
#endif 
         IF (LNANINFCHK) THEN
           IF (SUM(IMATRA) .NE. SUM(IMATRA) .OR. SUM(IMATDA) .NE. SUM(IMATDA)) THEN
             IF (WRITEDBGFLAG == 1) WRITE(DBG%FHNDL,*) 'NAN AT GRIDPOINT', IP, '   DUE TO SIN', SUM(IMATRA), SUM(IMATDA)
             CALL WWM_ABORT('NAN in wwm_sourceterm.F90 l.311')
           END IF
         ENDIF

         IF ((ISELECT.EQ.2 .OR. ISELECT.EQ.10 .OR. ISELECT.EQ.20) .AND. .NOT. LRECALC) THEN
           IF (IOBP(IP) .EQ. 0) THEN
             IF (MESNL .EQ. 1) THEN
               CALL SNL41 (IP, KMWAM, ACLOC, IMATRA, IMATDA, SSNL4, DSSNL4)
             ELSE IF (MESNL .EQ. 2) THEN
               CALL SNL4(IP, KMWAM, ACLOC, IMATRA, IMATDA)
             ELSE IF (MESNL .EQ. 3) THEN
               CALL SNL42(IP, KMWAM, ACLOC, IMATRA, IMATDA)
             ELSE IF (MESNL .EQ. 4) THEN
               CALL SNL43(IP, KMWAM, ACLOC, IMATRA, IMATDA)
             ELSE IF (MESNL .EQ. 5) THEN
               ACLOC_WRT = REAL(ACLOC)
               SPSIG_WRT = REAL(SPSIG)
               SPDIR_WRT = REAL(SPDIR)
               DEP_WRT   = DEP(IP)
               CALL WWMQUAD_WRT (ACLOC_WRT,SPSIG_WRT,SPDIR_WRT,MDC,MSC,DEP_WRT,1,XNL,DDIAG,IERR_WRT)
               IF (IERR_WRT .GT. 0) THEN
                 WRITE (DBG%FHNDL,*) 'XNL_WRT ERROR', IERR_WRT
                 CALL WWM_ABORT('XNL_WRT ERROR')
               ELSE
                 IMATRA(:,:) = IMATRA(:,:) + XNL (:,:)
                 IMATDA(:,:) = IMATDA(:,:) + DDIAG(:,:)
               END IF
             ELSE IF (MESNL .EQ. 6) THEN
#ifdef SNL4_TSA
               DO IS = 1, MSC
                 DO ID = 1, MDC
                   TMPAC_TSA(ID,IS) = ACLOC(IS,ID) * CG(IS,IP)
                 END DO
               END DO
               CG_TSA = CG(:,IP)
               WK_TSA = WK(:,IP)
               DEP_TSA = DEP(IP)
               NZZ = (MSC*(MSC+1))/2
               CALL W3SNLX ( TMPAC_TSA, CG_TSA, WK_TSA, DEP_TSA, NZZ, IMATRA_TSA, IMATDA_TSA)
               DO IS = 1, MSC
                 DO ID = 1, MDC
                   IMATRA(IS,ID) = IMATRA(IS,ID) + IMATRA_TSA(ID,IS) / CG(IS,IP)
                   IMATDA(IS,ID) = IMATDA(IS,ID) + IMATDA_TSA(ID,IS)
                 END DO
               END DO
#endif
             END IF
           END IF
           IF (IDISP == IP) THEN
             !CALL PLOT_SHADED_CONTOUR_POLAR(SPSIG/PI2,SPDIR*RADDEG,MSC,MDC,IMATRA/MAXVAL(IMATRA),50,MSC,MDC,'BEFORE ANY CALL')
             !PAUSE
           END IF
         END IF ! ISELECT

         IF (LNANINFCHK) THEN
           IF (SUM(IMATRA) .NE. SUM(IMATRA) .OR. SUM(IMATDA) .NE. SUM(IMATDA)) THEN
             IF (WRITEDBGFLAG == 1) WRITE(DBG%FHNDL,*) 'NAN AT GRIDPOINT', IP, 'DUE TO SNL4'
             CALL WWM_ABORT('NAN at wwm_sourceterms.F90 l.368')
           END IF
         ENDIF

#ifdef DEBUG
         WRITE(12,*) '2', SUM(IMATRA), SUM(IMATDA)
#endif


#ifdef TIMINGS
         call WAV_MY_WTIME(TIME4)
#endif 
         IF ((ISELECT.EQ.3 .OR. ISELECT.EQ.10 .OR. ISELECT.EQ.20) .AND. .NOT. LRECALC) THEN
           
           IMATRAT = IMATRA

           IF (IOBP(IP) .EQ. 0 .OR. IOBP(IP) .EQ. 4) THEN

             IF (MESDS == 1) THEN
#ifdef ST41
               CALL W3SDS4_OLD(AWW3,WK(:,IP),CG(:,IP),UFRIC(IP),USTDIR(IP),DEP(IP),IMATRA1D, IMATDA1D) 
#elif ST42
               CALL W3SDS4(AWW3,WK(:,IP),CG(:,IP),UFRIC(IP),USTDIR(IP),DEP(IP),IMATRA1D,IMATDA1D,BRLAMBDA,WHITECAP)
#endif
               CALL ONED2TWOD(IMATRA1D,IMATRAWW3)
               CALL ONED2TWOD(IMATDA1D,IMATDAWW3)
               DO ID = 1, MDC
                 SSDS(:,ID)   = IMATRAWW3(:,ID) / CG(:,IP)
                 IF (ISELECT.EQ.3 .OR. ISELECT.EQ.20) THEN
                   IMATRA(:,ID) = IMATRA(:,ID)+IMATRAWW3(:,ID) / CG(:,IP)
                   !IMATDA(:,ID) = IMATDA(:,ID) + IMATDAWW3(:,ID) !/ CG(:,IP)
                 ELSE IF (ISELECT.EQ.10) THEN
                   IMATDA(:,ID) = IMATDA(:,ID) - IMATDAWW3(:,ID)
                 END IF
               END DO
               CALL COMPUTE_WAVE_SDSTOT(SSDS,SDSTOT)
               WAVE_SDSTOT(IP) = SDSTOT
             ELSE IF (MESDS == 2) THEN
               CALL WWM_ABORT('PLEASE USE LSOURCEWAM FOR ECWAM SOURCE TERM FORMULATION')
             ELSE IF (MESDS == 3) THEN
               CALL SDS_NEDWAM_CYCLE4( IP, KMWAM, SME01, ETOT, ACLOC, IMATRA, IMATDA, SSDS  ) ! NEDWAM 
             ELSE IF (MESDS == 4) THEN
               ABAB = 1.
               LPOW = 4.
               MPOW = 4.
               a1  = 0.00000045
               a2  = 0.0000095
!              LPOW = 2.
!              MPOW = 2.
!              a1  = 0.0002
!              a2  = 0.003
!  0.0002 0.003 2.0 2.0 KOM
!  0.00000045 0.0000095 4.0 4.0 BD
               DO IS = 1, MSC
                 EDENS(IS) = 0.
                 DO ID = 1, MDC
                   EDENS(IS) = EDENS(IS) + ACLOC(IS,ID) *  SPSIG(IS) * PI2 * DDIR
                 END DO
               END DO
               CALL CALC_SDS(IP,MSC,EDENS,FR,Kds,ABAB,LPOW,MPOW,a1,a2,ACLOC,IMATRA,IMATDA,SSDS)
!              CALL SSWELL(IP,ETOT,ACLOC,IMATRA,IMATDA,URMSTOP,CG0)
             ELSE IF (MESDS == 5) THEN
               CALL SDS_CYCLE3 ( IP, KMWAM, SME10, ETOT, ACLOC, IMATRA, IMATDA, SSDS ) 
             END IF
             IF ( MESDS/=0 .AND. RADFLAG .EQ. 'VOR' ) THEN
                 CALL COMPUTE_SDS(IP,SSDS)
             END IF
           END IF
#ifdef VDISLIN
           IF (IDISP == IP) THEN
             CALL PLOT_SHADED_CONTOUR_POLAR(SPSIG/PI2,SPDIR*RADDEG,MSC,MDC,IMATRAT-IMATRA,50,MSC,MDC,'BEFORE ANY CALL')
           END IF
#endif
         END IF

#ifdef DEBUG
         WRITE(12,*) '3', SUM(IMATRA), SUM(IMATDA)
#endif

#ifdef TIMINGS
         call WAV_MY_WTIME(TIME5)
#endif 
         IF (((ISELECT.EQ.4 .OR. ISELECT.EQ.10) .AND.ISHALLOW(IP).EQ.1) .AND. .NOT. LRECALC) THEN
           IF (SUMACLOC .GT. VERYSMALL) THEN
             IF ((MESTR .EQ. 1).or.(MESTR .eq. 5) ) THEN
               CALL TRIAD_ELDEBERKY(IP,HS,SME01,ACLOC,IMATRA,IMATDA,SSNL3,DSSNL3)
             ELSE IF (MESTR .EQ. 2) THEN
               CALL SNL31 (IP,HS,SME01,ACLOC,IMATRA,IMATDA,SSNL3)
             ELSE IF (MESTR .EQ. 3) THEN
               CALL SNL32 (IP,HS,SME01,ACLOC,IMATRA,IMATDA,SSNL3)
             ELSE IF (MESTR .EQ. 4) THEN
               CALL SNL33 (IP,HS,SME01,ACLOC,IMATRA,IMATDA,SSNL3)
             ELSE IF (MESTR .EQ. 6) THEN
               CALL WWM_ABORT('NOT READ YET')
               CALL TRIAD_DINGEMANS (IP,ACLOC,IMATRA,IMATDA,SSNL3)
             ELSE IF (MESTR .EQ. 7) THEN
               CALL WWM_ABORT('NOT READ YET')
               CALL snl3ta(ip,snl3,dsnl3)
               SSNL3 = SNL3
               IMATRA = IMATRA + SNL3
               IMATDA = IMATDA + DSNL3
             END IF
           END IF
         END IF

#ifdef DEBUG
        WRITE(12,*) '4', SUM(IMATRA), SUM(IMATDA)
#endif


#ifdef TIMINGS
         call WAV_MY_WTIME(TIME6)
#endif 
         IF (MESBR .EQ. 1) THEN
           IF (((ISELECT.EQ.5 .OR. ISELECT.EQ.10 .OR. ISELECT.EQ.30) .AND. ISHALLOW(IP) .EQ. 1) .AND. .NOT. LRECALC) THEN
             IF (IOBP(IP) == 0 .OR. IOBP(IP) == 4 .OR. IOBP(IP) == 3) THEN
               ! Compute source term
               IF (IBREAK .EQ. 1 .OR. IBREAK .EQ. 4) THEN
                 CALL SDS_SWB(IP,SME01,KMWAM,ETOT,HS,ACLOC,IMATRA,IMATDA,SSBR,DSSBR)
               ELSE
                 CALL SDS_SWB(IP,SMECP,KMWAM,ETOT,HS,ACLOC,IMATRA,IMATDA,SSBR,DSSBR)
               ENDIF
               ! In case of VOR + implicit mode, we deal with the SBR term here
               ! In implicit, ISELECT.EQ.10 => call from SOURCE_INT_IMP_WWM
               ! In exlicit, we should use SMETHOD = 6, and SBR is computed in SUBROUTINE INT_SHALLOW_SOURCETERMS
               IF ( ISELECT.EQ.10 ) THEN 
                 CALL COMPUTE_WAVE_SBRTOT(-(DSSBR*ACLOC-SSBR),SBRTOT)
                 WAVE_SBRTOT(IP) = SBRTOT
                 EPS_W(IP) = - WAVE_SBRTOT(IP)/RHOW
                 IF (IROLLER == 0) EPS_BR(IP) = EPS_W(IP)
               ENDIF
               IF ( ISELECT.EQ.10 .AND. RADFLAG .EQ. 'VOR' ) CALL COMPUTE_SBR(IP, -(DSSBR*ACLOC-SSBR))
             END IF
           ENDIF
         END IF
#ifdef TIMINGS
         call WAV_MY_WTIME(TIME7)
#endif 

#ifdef DEBUG
         WRITE(12,*) '5', SUM(IMATRA), SUM(IMATDA)
#endif

         IF (MESBF .GE. 1) THEN
           IF (((ISELECT.EQ.6 .OR. ISELECT.EQ.10 .OR. ISELECT.EQ.30) .AND. ISHALLOW(IP) .EQ. 1) .AND. .NOT. LRECALC) THEN
             IF (IOBP(IP) == 0 .OR. IOBP(IP) == 4 .OR. IOBP(IP) == 3) THEN
              CALL SDS_BOTF(IP,ACLOC,IMATRA,IMATDA,SSBF,DSSBF)
              ! In case of VOR + implicit mode, we deal with the SBF term here
              ! In implicit, ISELECT.EQ.10 => call from SOURCE_INT_IMP_WWM
              ! In exlicit, we should use SMETHOD = 6, and SBF is computed in SUBROUTINE INT_SHALLOW_SOURCETERMS
              ! In contrast to the breaking source term, the friction source term is not linearized (hence SSBF) 
              IF ( ISELECT.EQ.10 ) THEN 
                CALL COMPUTE_WAVE_SBFTOT(SSBF,SBFTOT)
                WAVE_SBFTOT(IP) = SBFTOT
              ENDIF
              IF ( ISELECT.EQ.10 .AND. RADFLAG .EQ. 'VOR' ) THEN
                 CALL COMPUTE_SBF(IP,SSBF)
              END IF
             END IF
           END IF
         ENDIF

         IF (MEVEG .GE. 1) THEN
           IF (((ISELECT.EQ.6 .OR. ISELECT.EQ.10 .OR. ISELECT.EQ.30) .AND. ISHALLOW(IP) .EQ. 1) .AND. .NOT. LRECALC) THEN
             IF (IOBP(IP) == 0 .OR. IOBP(IP) == 4 .OR. IOBP(IP) == 3) THEN
              CALL VEGDISSIP (IP,IMATRA,IMATDA,SSVEG,DSSVEG,ACLOC,DEP(IP),ETOT,SME01,KME01)
             END IF
           END IF
         ENDIF

         IF (LNANINFCHK) THEN
           IF (SUM(IMATRA) .NE. SUM(IMATRA) .OR. SUM(IMATDA) .NE. SUM(IMATDA)) THEN
             IF (WRITEDBGFLAG == 1) WRITE(DBG%FHNDL,*) 'NAN AT GRIDPOINT', IP, '   DUE TO SBF' 
             CALL WWM_ABORT('NAN in wwm_sourceterms.F90 at l.419')
           END IF
         ENDIF

#ifdef DEBUG
         WRITE(12,*) '6', SUM(IMATRA), SUM(IMATDA)
#endif

#ifdef TIMINGS
        call WAV_MY_WTIME(TIME8)
#endif 
!------------------------------------------------------------------------------------------------------------------------!
!-------------------------------- RECALCULATE ALL SOURCE TERMS BASED ON THE NEW SPECTRA ---------------------------------! 
!------------------------------------------------------------------------------------------------------------------------!
         IF (LRECALC .and. IOBP(IP) .EQ. 0) THEN

           IF (MESIN == 1) THEN
             AS      = 0.
#ifdef ST41
             CALL W3SPR4_OLD ( AWW3, CG(:,IP), WK(:,IP), EMEAN, FMEAN, WNMEAN, AMAX, WIND10, WINDTH, UFRIC(IP), USTDIR(IP), TAUWX(IP), TAUWY(IP), CD(IP), Z0(IP), ALPHA_CH(IP), LLWS, FMEANWS)
             CALL W3SIN4_OLD ( AWW3, CG(:,IP), WK(:,IP), WIND10, UFRIC(IP), RHOAW, AS, WINDTH, Z0(IP), CD(IP), TAUWX(IP), TAUWY(IP), TAUWAX, TAUWAY, IMATRA1D, IMATDA1D, LLWS)
             CALL W3SPR4_OLD ( AWW3, CG(:,IP), WK(:,IP), EMEAN, FMEAN, WNMEAN, AMAX, WIND10, WINDTH, UFRIC(IP), USTDIR(IP), TAUWX(IP), TAUWY(IP), CD(IP), Z0(IP), ALPHA_CH(IP), LLWS, FMEANWS)
#elif ST42
             CALL W3SPR4 ( AWW3, CG(:,IP), WK(:,IP), EMEAN, FMEAN, FMEAN1, WNMEAN, AMAX, WIND10, WINDTH, UFRIC(IP), USTDIR(IP), TAUWX(IP), TAUWY(IP), CD(IP), Z0(IP), ALPHA_CH(IP), LLWS, FMEANWS)
             CALL W3SIN4 ( IP, AWW3, CG(:,IP), WN2,  WIND10, UFRIC(IP), RHOAW, AS, WINDTH, Z0(IP), CD(IP), TAUWX(IP), TAUWY(IP), TAUWAX, TAUWAY, IMATRA1D, IMATDA1D, LLWS, BRLAMBDA)
             CALL W3SPR4 ( AWW3, CG(:,IP), WK(:,IP), EMEAN, FMEAN, FMEAN1, WNMEAN, AMAX, WIND10, WINDTH, UFRIC(IP), USTDIR(IP), TAUWX(IP), TAUWY(IP), CD(IP), Z0(IP), ALPHA_CH(IP), LLWS, FMEANWS)
#else
             IF (WRITEDBGFLAG == 1) WRITE(DBG%FHNDL,*) 'NO ST42 or ST41 chosen but MESIN == 1'
             CALL WWM_ABORT('stop wwm_sourceterms l.186')
#endif
           ELSEIF (MESIN == 2) THEN
           ELSEIF (MESIN == 3) THEN
           ELSEIF (MESIN == 4) THEN
           ELSEIF (MESIN == 5) THEN
           ENDIF

           ! Only used for nesting, will be removed in the future
           !DISSIPATION(IP) = 0.
           !AIRMOMENTUM(IP) = 0.
           !DO ID = 1, MDC
           !  TMP_DS = ( SSBR(:,ID) + SSBF(:,ID) + SSDS(:,ID) ) * SPSIG * DDIR
           !  TMP_IN = ( SSINE(:,ID) + SSINL(:,ID) ) * SPSIG * DDIR
           !  DO IS = 2, MSC
           !    DISSIPATION(IP) = DISSIPATION(IP) + ONEHALF * ( TMP_DS(IS) + TMP_DS(IS-1) ) * DS_INCR(IS)
           !    AIRMOMENTUM(IP) = AIRMOMENTUM(IP) + ONEHALF * ( TMP_IN(IS) + TMP_IN(IS-1) ) * DS_INCR(IS)
           !  END DO
           !END DO

         ENDIF

#ifdef TIMINGS 
         call WAV_MY_WTIME(TIME9)
#ifdef WWM_MPI 
         IF (IP == MNP .AND. myrank == 0 ) THEN
#else 
         IF (IP == MNP) THEN
#endif
           IF (WRITESTATFLAG == 1) THEN
             WRITE(STAT%FHNDL,'("+TRACE...",A,F15.6)') '-----SOURCE TIMINGS-----'
             WRITE(STAT%FHNDL,'("+TRACE...",A,F15.6)') 'PREPARATION        ', TIME2-TIME1
             WRITE(STAT%FHNDL,'("+TRACE...",A,F15.6)') 'SIN                ', TIME3-TIME2
             WRITE(STAT%FHNDL,'("+TRACE...",A,F15.6)') 'SDS                ', TIME4-TIME3
             WRITE(STAT%FHNDL,'("+TRACE...",A,F15.6)') 'SNL4               ', TIME5-TIME4
             WRITE(STAT%FHNDL,'("+TRACE...",A,F15.6)') 'SNL3               ', TIME6-TIME5
             WRITE(STAT%FHNDL,'("+TRACE...",A,F15.6)') 'SBR                ', TIME7-TIME6
             WRITE(STAT%FHNDL,'("+TRACE...",A,F15.6)') 'SBF                ', TIME8-TIME7
             WRITE(STAT%FHNDL,'("+TRACE...",A,F15.6)') 'RECALC             ', TIME9-TIME8
             WRITE(STAT%FHNDL,'("+TRACE...",A,F15.6)') 'TOTAL              ', TIME9-TIME1
             WRITE(STAT%FHNDL,'("+TRACE...",A,F15.6)') '------END-TIMINGS-  ---'
             FLUSH(STAT%FHNDL)
           END IF
         END IF
#endif 

      END SUBROUTINE
!**********************************************************************
!*                                                                    *
!**********************************************************************
#ifdef SCHISM
    SUBROUTINE COMPUTE_WAVE_SBFTOT(SSBF,SBFTOT)
      USE DATAPOOL
      IMPLICIT NONE
      INTEGER                  :: IS, ID
      REAL(rkind), INTENT(IN)  :: SSBF(MSC,MDC)
      REAL(rkind), INTENT(OUT) :: SBFTOT

      ! Initialization
      SBFTOT = ZERO
       
      ! Loop over frequencies and directions
      DO IS = 1, MSC
        DO ID = 1, MDC
          SBFTOT = SBFTOT + SPSIG(IS)*SSBF(IS,ID)*DS_INCR(IS)*DDIR*G9*RHOW
        END DO
      END DO
    END SUBROUTINE
#endif
!**********************************************************************
!*                                                                    *
!**********************************************************************
#ifdef SCHISM
    SUBROUTINE COMPUTE_WAVE_SBRTOT(SSBR, SBRTOT)
      !This serves mostly for implicit solvers. For the explicit
      !solver, a more accurate computation is performed during 
      !the call to routine INT_SHALLOW_SOURCETERMS.
      !As the source term is linearized with implicit solver,
      !take as input : DSSBR*ACLOC - SSBR
      USE DATAPOOL
      IMPLICIT NONE
      INTEGER                  :: IS, ID
      REAL(rkind), INTENT(IN)  :: SSBR(MSC,MDC)
      REAL(rkind), INTENT(OUT) :: SBRTOT
	  
      ! Initialization
      SBRTOT = ZERO

      ! Loop over frequencies and directions
      DO IS = 1, MSC
        DO ID = 1, MDC
          SBRTOT = SBRTOT + SPSIG(IS)*SSBR(IS,ID)*DS_INCR(IS)*DDIR*G9*RHOW
        END DO
      END DO
    END SUBROUTINE
#endif
!**********************************************************************
!*                                                                    *
!**********************************************************************
#ifdef SCHISM
    SUBROUTINE COMPUTE_WAVE_SINTOT(SSIN,SINTOT)
      USE DATAPOOL
      IMPLICIT NONE
      INTEGER                  :: IS, ID
      REAL(rkind), INTENT(IN)  :: SSIN(MSC,MDC)
      REAL(rkind), INTENT(OUT) :: SINTOT
	  
      ! Initialization
      SINTOT = ZERO

      ! Loop over frequencies and directions
      DO IS = 1, MSC
        DO ID = 1, MDC
          SINTOT = SINTOT + SPSIG(IS)*SSIN(IS,ID)*DS_INCR(IS)*DDIR*G9*RHOW
        END DO
      END DO
    END SUBROUTINE
#endif
!**********************************************************************
!*                                                                    *
!**********************************************************************
#ifdef SCHISM
    SUBROUTINE COMPUTE_WAVE_SDSTOT(SSDS,SDSTOT)
      USE DATAPOOL
      IMPLICIT NONE
      INTEGER                 :: IS, ID
      REAL(rkind), INTENT(IN) :: SSDS(MSC,MDC)
      REAL(rkind), INTENT(OUT):: SDSTOT
	  
      ! Initialization
      SDSTOT = ZERO

      ! Loop over frequencies and directions
      DO IS = 1, MSC
        DO ID = 1, MDC
          SDSTOT = SDSTOT + SPSIG(IS)*SSDS(IS,ID)*DS_INCR(IS)*DDIR*G9*RHOW
        END DO
      END DO
    END SUBROUTINE
#endif
!**********************************************************************
!*                                                                    *
!**********************************************************************
#ifdef SCHISM
    SUBROUTINE COMPUTE_SDS(IP,SSDS)
      USE DATAPOOL
      IMPLICIT NONE
      INTEGER                 :: IS, ID
      INTEGER, INTENT(IN)     :: IP
      REAL(rkind), INTENT(IN) :: SSDS(MSC,MDC)

      ! Initialization
      SDS(:,IP) = ZERO

      ! Loop over frequencies and directions
      DO IS = 1, MSC
        DO ID = 1, MDC
          SDS(1,IP) = SDS(1,IP) + G9*COSTH(ID)*WK(IS,IP)*SSDS(IS,ID)*DS_INCR(IS)*DDIR
          SDS(2,IP) = SDS(2,IP) + G9*SINTH(ID)*WK(IS,IP)*SSDS(IS,ID)*DS_INCR(IS)*DDIR
        END DO
      END DO
    END SUBROUTINE
#endif