#include "wwm_functions.h"
!**********************************************************************
!*                                                                    *
!**********************************************************************
      SUBROUTINE SOURCE_INT_EXP
         USE DATAPOOL
#ifdef ST41
         USE W3SRC4MD_OLD, ONLY : LFIRSTSOURCE
#elif ST42
         USE W3SRC4MD, ONLY : LFIRSTSOURCE
#endif
         IMPLICIT NONE

         INTEGER        :: IP, IS, ID, I, K, M
         INTEGER        :: NIT_SIN, NIT_SDS, NIT_SNL4, NIT_SNL3, NIT_SBR, NIT_SBF, NIT_ALL
         INTEGER, SAVE  :: IFIRST, ISELECT
         REAL(rkind)    :: ACLOC(MSC,MDC), IMATRA(MSC,MDC), IMATDA(MSC,MDC), SSBRL2(MSC,MDC)
         REAL(rkind)    :: DT4S_T, DT4S_E, DT4S_Q, DT4S_H, DT4S_TQ, DT4S_TS, VEC2RAD

!$OMP PARALLEL DEFAULT(NONE)  &
!$OMP&         SHARED(MNP,MSC,MDC,DEP,DMIN,IOBP,SMETHOD, &
!$OMP&         DT4S,DTMIN_SDS,DTMIN_SIN,DTMIN_SBR,DTMIN_DYN, &
!$OMP&         DTMIN_SNL3, DTMIN_SNL4, DTMIN_SBF, WINDXY, WINDFAC,&
!$OMP&         LSOUBOUND,ISHALLOW,LADVTEST,LMAXETOT, &
!$OMP&         NDYNITER_SIN,NDYNITER_SNL4, NDYNITER_SDS, &
!$OMP&         NDYNITER_SBR, NDYNITER_SNL3, NDYNITER_SBF, &
!$OMP&         NDYNITER,LSOURCESWAM,LLIMT,MESTR,MESBR,MESBF,AC2,AC1,&
!$OMP&         FL,SL,FL3,THWNEW,THWOLD,U10NEW,U10OLD,Z0NEW,Z0OLD) &
!$OMP&         PRIVATE(IP,IS,ID,ACLOC,NIT_SIN,NIT_SDS,&
!$OMP&         NIT_SNL4,NIT_SNL3,NIT_SBR,NIT_SBF,NIT_ALL,SSBRL2,&
!$OMP&         IMATRA,IMATDA,VEC2RAD)
!$OMP DO SCHEDULE(DYNAMIC,1)
         DO IP = 1, MNP
!           IF (IP_IS_STEADY(IP) .EQ. 1) CYCLE
           IF (DEP(IP) .LT. DMIN) CYCLE
           IF (IOBP(IP) .EQ. 0) THEN
               ACLOC  = AC2(:,:,IP)
               IF (SMETHOD == 1) THEN
                 IF (MESBR .GT. 0 .OR. MESBF .GT. 0) CALL RKS_SP3(IP,DT4S,.FALSE.,ACLOC,30)
                 IF (LSOURCESWAM) THEN
                   CALL SOURCE_INT_EXP_WAM(IP, ACLOC)  
                 ELSE
                   CALL INT_IP_STAT(IP,DT4S,LLIMT,ACLOC,20)
                 ENDIF
                 IF (MESTR .GT. 0) CALL INT_IP_DYN(IP, DT4S, LLIMT, DTMIN_DYN, NDYNITER, ACLOC, NIT_ALL, 4)
               ELSE IF (SMETHOD == 2) THEN
                 CALL INT_IP_STAT(IP,DT4S, LLIMT,ACLOC, 10)
               ELSE IF (SMETHOD == 3) THEN
                 CALL RKS_SP3(IP,DT4S,LLIMT,ACLOC, 10)
               ELSE IF (SMETHOD == 4) THEN
                 CALL INT_IP_DYN(IP, DT4S, LLIMT, DTMIN_DYN, NDYNITER, ACLOC, NIT_ALL, 10)
               ELSE IF (SMETHOD == 5) THEN ! Full splitting of all source embedded within a dynamic RK-3 Integration ... 
                 CALL INT_IP_DYN(IP, DT4S, LLIMT, DTMIN_SIN,  NDYNITER_SIN  , ACLOC, NIT_SIN, 1) ! Sin
                 CALL INT_IP_DYN(IP, DT4S, LLIMT, DTMIN_SNL4,  NDYNITER_SNL4 , ACLOC, NIT_SNL4, 2)! Snl4b
                 CALL INT_IP_DYN(IP, DT4S, LLIMT, DTMIN_SDS,  NDYNITER_SDS  , ACLOC, NIT_SDS, 3) ! Sds
                 CALL INT_IP_DYN(IP, DT4S, LLIMT, DTMIN_SNL3,  NDYNITER_SNL3 , ACLOC, NIT_SNL3, 4)! Snl3
                 CALL INT_IP_DYN(IP, DT4S, LLIMT, DTMIN_SBR,  NDYNITER_SBR  , ACLOC, NIT_SBR, 5) ! Sbr
                 CALL INT_IP_DYN(IP, DT4S, LLIMT, DTMIN_SBF,  NDYNITER_SBF  , ACLOC, NIT_SBF, 6) ! Sbf
               END IF
               CALL SOURCETERMS(IP, ACLOC, IMATRA, IMATDA, .TRUE.,1,'RECALC DOMAIN') ! Update everything based on the new spectrum ...
               IF (LMAXETOT .AND. .NOT. LADVTEST .AND. ISHALLOW(IP) .EQ. 1) THEN
                 CALL BREAK_LIMIT(IP,ACLOC,SSBRL2)
               ENDIF
               AC2(:,:,IP) = ACLOC
           ELSE !Boundary node ... 
             IF (LSOUBOUND) THEN ! Source terms on boundary ...
               IF (IOBP(IP) .NE. 2) THEN
                 ACLOC  = AC2(:,:,IP)
                 IF (SMETHOD == 1) THEN
                   CALL RKS_SP3(IP,DT4S,.FALSE.,ACLOC,30)
                   IF (LSOURCESWAM) THEN
                     CALL SOURCE_INT_EXP_WAM(IP, ACLOC)
                   ELSE
                     CALL INT_IP_STAT(IP,DT4S,LLIMT,ACLOC,20)
                   ENDIF
                   CALL INT_IP_DYN(IP, DT4S, LLIMT, DTMIN_DYN, NDYNITER, ACLOC, NIT_ALL,4)
                 ELSE IF (SMETHOD == 2) THEN
                   CALL INT_IP_STAT(IP,DT4S, LLIMT,ACLOC,10)
                 ELSE IF (SMETHOD == 3) THEN
                   CALL RKS_SP3(IP,DT4S,LLIMT,ACLOC,10)
                 ELSE IF (SMETHOD == 4) THEN
                   CALL INT_IP_DYN(IP, DT4S, LLIMT, DTMIN_DYN, NDYNITER, ACLOC, NIT_ALL,10)
                 ELSE IF (SMETHOD == 5) THEN ! Full splitting of all source embedded within a dynamic RK-3 Integration ... 
                   CALL INT_IP_DYN(IP, DT4S, LLIMT, DTMIN_SIN,  NDYNITER_SIN  , ACLOC, NIT_SIN,1) ! Sin
                   CALL INT_IP_DYN(IP, DT4S, LLIMT, DTMIN_SNL4,  NDYNITER_SNL4 , ACLOC, NIT_SNL4,2)! Snl4b
                   CALL INT_IP_DYN(IP, DT4S, LLIMT, DTMIN_SDS,  NDYNITER_SDS  , ACLOC, NIT_SDS,3) ! Sds
                   CALL INT_IP_DYN(IP, DT4S, LLIMT, DTMIN_SNL3,  NDYNITER_SNL3 , ACLOC, NIT_SNL3,4)! Snl3
                   CALL INT_IP_DYN(IP, DT4S, LLIMT, DTMIN_SBR,  NDYNITER_SBR  , ACLOC, NIT_SBR,5) ! Sbr
                   CALL INT_IP_DYN(IP, DT4S, LLIMT, DTMIN_SBF,  NDYNITER_SBF  , ACLOC, NIT_SBF,6) ! Sbf
                 END IF
                 CALL SOURCETERMS(IP, ACLOC, IMATRA, IMATDA, .TRUE., 1, 'RECALC BOUNDARY') ! Update everything based on the new spectrum ...
                 IF (LMAXETOT .AND. .NOT. LADVTEST .AND. ISHALLOW(IP) .EQ. 1) THEN
                   CALL BREAK_LIMIT(IP,ACLOC,SSBRL2)
                 ENDIF
                 AC2(:,:,IP) = ACLOC
               ENDIF
             ELSE
               ACLOC = AC2(:,:,IP)
               IF (LMAXETOT .AND. .NOT. LADVTEST .AND. ISHALLOW(IP) .EQ. 1) THEN ! limit wave height on the boundary ...
                 CALL BREAK_LIMIT(IP,ACLOC,SSBRL2)
                 AC2(:,:,IP) = ACLOC
               ENDIF
             ENDIF
           ENDIF
           AC1(:,:,IP) = AC2(:,:,IP)
         ENDDO
!$OMP END PARALLEL 

#if defined ST41 || defined ST42
         LFIRSTSOURCE = .FALSE.
#endif
      END SUBROUTINE
!**********************************************************************
!*                                                                    *
!**********************************************************************
      SUBROUTINE SOURCE_INT_IMP_WWM
         USE DATAPOOL
#ifdef ST41
         USE W3SRC4MD_OLD, ONLY : LFIRSTSOURCE
#endif
#ifdef ST42
         USE W3SRC4MD, ONLY : LFIRSTSOURCE
#endif
         IMPLICIT NONE

         INTEGER :: IP

         REAL(rkind)    :: ACLOC(MSC,MDC)
         REAL(rkind)    :: IMATRA(MSC,MDC), IMATDA(MSC,MDC)

!$OMP WORKSHARE
         IMATDAA = 0.
         IMATRAA = 0.
!$OMP END WORKSHARE

!$OMP PARALLEL DEFAULT(NONE)  &
!$OMP&         SHARED(MNP,MSC,MDC,DEP,DMIN,IOBP,SMETHOD, &
!$OMP&         LSOUBOUND) &
!$OMP&         PRIVATE(IP,ACLOC,AC2,IMATDA,IMATRA,IMATRAA,IMATDAA)
!$OMP DO SCHEDULE(DYNAMIC,1)
         DO IP = 1, MNP
           IF (DEP(IP) .LT. DMIN) CYCLE
           IF ((ABS(IOBP(IP)) .NE. 1 .AND. IOBP(IP) .NE. 3)) THEN
             IF ( DEP(IP) .GT. DMIN .AND. IOBP(IP) .NE. 2) THEN
               ACLOC = AC2(:,:,IP)
               CALL SOURCETERMS(IP, ACLOC, IMATRA, IMATDA, .FALSE., 10, 'SOURCE_INT_IMP_WWM DOMAIN') 
               !CALL CYCLE3(IP, ACLOC, IMATRA, IMATDA)
               IMATDAA(:,:,IP) = IMATDA 
               IMATRAA(:,:,IP) = IMATRA 
             ENDIF
           ELSE
             IF (LSOUBOUND) THEN ! Source terms on boundary ...
               IF ( DEP(IP) .GT. DMIN .AND. IOBP(IP) .NE. 2) THEN
                 ACLOC = AC2(:,:,IP)
                 CALL SOURCETERMS(IP, ACLOC, IMATRA, IMATDA, .FALSE.,10, 'SOURCE_INT_IMP_WWM BOUNDARY')
                 !CALL CYCLE3(IP, ACLOC, IMATRA, IMATDA)
                 IMATDAA(:,:,IP) = IMATDA
                 IMATRAA(:,:,IP) = IMATRA 
               ENDIF
             ENDIF
           ENDIF
         END DO
!$OMP END PARALLEL 
#if defined ST41 || defined ST42
         LFIRSTSOURCE = .FALSE.
#endif

!         PAUSE
      END SUBROUTINE
!**********************************************************************
!*                                                                    *
!**********************************************************************      
      SUBROUTINE SOURCE_INT_EXP_WAM(IP, ACLOC)
         USE DATAPOOL
         IMPLICIT NONE
         
         INTEGER, INTENT(IN) :: IP
         REAL(rkind), INTENT(INOUT) :: ACLOC(MSC,MDC)
         
         INTEGER      :: IS, ID, IMETHOD, ICODE
         REAL(rkind)  :: VEC2RAD
         REAL(rkind)  :: IMATRA(MSC,MDC), IMATDA(MSC,MDC)
         REAL(rkind)  :: ETOT,SME01,SME10,KME01,KMWAM,KMWAM2,HS,WIND10
         REAL(rkind)  :: ETAIL,EFTAIL,EMAX,LIMAC,NEWDAC,FPM,WINDTH,TEMP,GTEMP1
         REAL(rkind)  :: RATIO,LIMFAC,LIMDAC,GTEMP2,FLHAB,DELFL,USFM, NEWDACDT
         REAL(rkind)  :: MAXDAC, MAXDACDT, MAXDACDTDA, SC, SP, DNEWDACDTDA, JAC, FF

         REAL(rkind),DIMENSION(MDC,MSC)  :: SSDS,DSSDS,SSNL4,DSSNL4,SSIN,DSSIN
     
         ICODE = 1 
         
         IF (MESIN .GT. 0 .OR. MESDS .GT. 0 .OR. MESNL .GT. 0) THEN
           DO IS = 1, MSC
             DO ID = 1, MDC
               FL3(IP,ID,IS) = ACLOC(IS,ID) * PI2 * SPSIG(IS)
               FL(IP,ID,IS)  = FL3(IP,ID,IS)
               SL(IP,ID,IS)  = FL(IP,ID,IS)
             END DO
           END DO
           THWOLD(IP,1) = THWNEW(IP)
           U10NEW(IP) = MAX(TWO,SQRT(WINDXY(IP,1)**2+WINDXY(IP,2)**2))*WINDFAC
           Z0NEW(IP) = Z0OLD(IP,1)
           THWNEW(IP) = VEC2RAD(WINDXY(IP,1),WINDXY(IP,2))
           IF (ICODE == 1) THEN
             CALL IMPLSCH (FL3(IP,:,:), FL(IP,:,:), IP, IP, 1, &
     &                     THWOLD(IP,1), USOLD(IP,1), &
     &                     TAUW(IP), Z0OLD(IP,1), &
     &                     ROAIRO(IP,1), ZIDLOLD(IP,1), &
     &                     U10NEW(IP), THWNEW(IP), USNEW(IP), &
     &                     Z0NEW(IP), ROAIRN(IP), ZIDLNEW(IP), &
     &                     SL(IP,:,:), FCONST(IP,:))
           ELSE IF (ICODE == 2) THEN
             CALL PREINTRHS (FL3(IP,:,:), FL(IP,:,:), IP, IP, 1, &
     &                       THWOLD(IP,1), USOLD(IP,1), &
     &                       TAUW(IP), Z0OLD(IP,1), &
     &                       ROAIRO(IP,1), ZIDLOLD(IP,1), &
     &                       U10NEW(IP), THWNEW(IP), USNEW(IP), &
     &                       Z0NEW(IP), ROAIRN(IP), ZIDLNEW(IP), &
     &                       SL(IP,:,:), FCONST(IP,:), FMEANWS(IP), MIJ(IP), &
     &                       SSDS, DSSDS, SSIN, DSSIN, &
     &                       SSNL4, DSSNL4)
             CALL INTSPECWAM (FL3(IP,:,:), FL(IP,:,:), IP, IP, 1, &
     &                      THWOLD(IP,1), USOLD(IP,1), &
     &                      TAUW(IP), Z0OLD(IP,1), &
     &                      ROAIRO(IP,1), ZIDLOLD(IP,1), &
     &                      U10NEW(IP), THWNEW(IP), USNEW(IP), &
     &                      Z0NEW(IP), ROAIRN(IP), ZIDLNEW(IP), &
     &                      SL(IP,:,:), FCONST(IP,:), FMEANWS(IP), MIJ(IP))
             CALL POSTINTRHS (FL3(IP,:,:), FL(IP,:,:), IP, IP, 1, &
     &                      THWOLD(IP,1), USOLD(IP,1), &
     &                      TAUW(IP), Z0OLD(IP,1), &
     &                      ROAIRO(IP,1), ZIDLOLD(IP,1), &
     &                      U10NEW(IP), THWNEW(IP), USNEW(IP), &
     &                      Z0NEW(IP), ROAIRN(IP), ZIDLNEW(IP), &
     &                      SL(IP,:,:), FCONST(IP,:), FMEANWS(IP), MIJ(IP))
           ELSE IF (ICODE == 3) THEN
             CALL IMPLSCH_LOCAL (IP, FL3(IP,:,:), FL(IP,:,:), 1, &
     &                     SL(IP,:,:))
           ENDIF ! ICODE
           DO IS = 1, MSC
             DO ID = 1, MDC
               ACLOC(IS,ID) =  MAX(ZERO, FL3(IP,ID,IS) / PI2 / SPSIG(IS))
             END DO
           END DO
         ENDIF
      END SUBROUTINE
!**********************************************************************
!*                                                                    *
!**********************************************************************
      SUBROUTINE SOURCE_INT_IMP_WAM_PRE
         USE DATAPOOL
         IMPLICIT NONE
         INTEGER      :: IP, IS, ID, IMETHOD
         REAL(rkind)  :: ACLOC(MSC,MDC), VEC2RAD
         REAL(rkind)  :: IMATRA(MSC,MDC), IMATDA(MSC,MDC)
         REAL(rkind)  :: SSNL3(MSC,MDC),DSSNL3(MSC,MDC), DSSVEG(MSC,MDC)
         REAL(rkind)  :: SSBR(MSC,MDC),DSSBR(MSC,MDC), SSVEG(MSC,MDC)
         REAL(rkind)  :: SSBF(MSC,MDC),DSSBF(MSC,MDC), SSINL(MSC,MDC)
         REAL(rkind)  :: ETOT,SME01,SME10,KME01,KMWAM,KMWAM2,HS,WIND10
         REAL(rkind)  :: ETAIL,EFTAIL,EMAX,LIMAC,NEWDAC,FPM,WINDTH,TEMP,GTEMP1
         REAL(rkind)  :: RATIO,LIMFAC,LIMDAC,GTEMP2,FLHAB,DELFL,USFM, NEWDACDT
         REAL(rkind)  :: MAXDAC, MAXDACDT, MAXDACDTDA, SC, SP, DNEWDACDTDA, JAC, FF

         REAL(rkind),DIMENSION(MDC,MSC)  :: SSDS,DSSDS,SSNL4,DSSNL4,SSIN,DSSIN

         IMETHOD = 4 

!$OMP WORKSHARE
         IMATDAA = ZERO
         IMATRAA = ZERO
!$OMP END WORKSHARE

!$OMP PARALLEL DEFAULT(NONE)  &
!$OMP&         SHARED(MNP,MSC,MDC,DEP,DMIN,IOBP,SMETHOD, &
!$OMP&         LSOUBOUND,LINID,MESTR,MESBR,MESBF,LIMFAK,&
!$OMP&         COFRM4,WINDFAC,MESIN,MESDS,MESNL,SPSIG) &
!$OMP&         PRIVATE(IP,ACLOC,AC2,IMATDA,IMATRA,IMATRAA,IMATDAA, &
!$OMP&         SSNL3,DSSNL3,SSBR,DSSBR,SSBF,DSSBF,FL3,FL,Z0NEW,&
!$OMP&         SSDS,DSSDS,JAC,USNEW,ZIDLNEW,FF,FMEAN,USFM,LIMFAC,&
!$OMP&         FCONST,ZIDLOLD,TAUW,USOLD,WINDXY,WIND10,WINDTH,&
!$OMP&         TEMP,HS,ETOT,FPM,ISHALLOW,KME01,SME10,KMWAM, &
!$OMP&         NEWDAC,KMWAM2,MAXDAC,WK,SC,CG,SME01,GTEMP2,&
!$OMP&         DELFL,FLHAB,GTEMP1,IMETHOD,DT4A,SSNL4,DSSNL4,&
!$OMP&         SSIN,DSSIN,FMEANWS,ROAIRN,MIJ,Z0OLD,ROAIRO, &
!$OMP&         THWOLD,THWNEW,U10NEW,SL,SSINL)
!$OMP DO SCHEDULE(DYNAMIC,1)
         DO IP = 1, MNP
           SSNL3 = ZERO; DSSNL3 = ZERO
           SSBR  = ZERO; DSSBR  = ZERO
           SSBF  = ZERO; DSSBF  = ZERO
           SSINL = ZERO
           IF (DEP(IP) .LT. DMIN) THEN
             IMATRAA(:,:,IP) = ZERO
             IMATDAA(:,:,IP) = ZERO
             CYCLE
           ENDIF
           IF (MESIN .GT. 0 .OR. MESDS .GT. 0 .OR. MESNL .GT. 0) THEN
             DO IS = 1, MSC
               DO ID = 1, MDC
                 FL3(IP,ID,IS) = AC2(IS,ID,IP) * PI2 * SPSIG(IS)
                 FL(IP,ID,IS)  = FL3(IP,ID,IS)
                 SL(IP,ID,IS)  = FL(IP,ID,IS)
               END DO
             END DO
             THWOLD(IP,1) = THWNEW(IP)
             U10NEW(IP) = MAX(TWO,SQRT(WINDXY(IP,1)**2+WINDXY(IP,2)**2))*WINDFAC
             Z0NEW(IP) = Z0OLD(IP,1)
             THWNEW(IP) = VEC2RAD(WINDXY(IP,1),WINDXY(IP,2))
             CALL PREINTRHS (FL3(IP,:,:), FL(IP,:,:), IP, IP, 1, &
     &                       THWOLD(IP,1), USOLD(IP,1), &
     &                       TAUW(IP), Z0OLD(IP,1), &
     &                       ROAIRO(IP,1), ZIDLOLD(IP,1), &
     &                       U10NEW(IP), THWNEW(IP), USNEW(IP), &
     &                       Z0NEW(IP), ROAIRN(IP), ZIDLNEW(IP), &
     &                       SL(IP,:,:), FCONST(IP,:), FMEANWS(IP), MIJ(IP), &
     &                       SSDS, DSSDS, SSIN, DSSIN, &
     &                       SSNL4, DSSNL4)
           ENDIF ! MESIN .EQ. 0 .AND. MESDS .AND. 0 .AND. MESNL = 0
           IF (IOBP(IP) .EQ. 0) THEN
             DO ID = 1, MDC
               DO IS = 1, MSC 
                 IF (AC2(IS,ID,IP) .LT. THR) CYCLE
                 JAC = ONE/PI2/SPSIG(IS)
                 !IMATDAA(IS,ID,IP) =  FL(IP,ID,IS) !... this is not working right, reason is unknown, there signchanges that should not be
                 !IMATRAA(IS,ID,IP) =  SL(IP,ID,IS)/PI2/SPSIG(IS) 
                 FF = FL3(IP,ID,IS)
                 !WRITE(11140,'(2I10,7E20.10)') IS, ID, FF, SSDS(ID,IS)/FF, DSSDS(ID,IS), SSIN(ID,IS)/FF, DSSNL4(ID,IS), SSNL4(ID,IS)/FF, DSSNL4(ID,IS) 
                 IF (IMETHOD == 0) THEN
                   IMATRAA(IS,ID,IP) =  SL(IP,ID,IS)/PI2/SPSIG(IS)
                   IMATDAA(IS,ID,IP) = ZERO 
                 ELSE IF (IMETHOD == 1) THEN 
                   IMATRAA(IS,ID,IP) = (SSIN(ID,IS)+SSDS(ID,IS)+SSNL4(ID,IS))*JAC
                 ELSE IF (IMETHOD == 2) THEN
                   IMATRAA(IS,ID,IP) = (SSIN(ID,IS)+SSNL4(ID,IS))*JAC
                   IMATDAA(IS,ID,IP) = -TWO*DSSDS(ID,IS)
                 ELSE IF (IMETHOD == 3) THEN
                   IF (SSIN(ID,IS) .GT. ZERO) THEN
                     IMATRAA(IS,ID,IP) = SSIN(ID,IS)*JAC
                   ELSE
                     IMATDAA(IS,ID,IP) = -TWO*DSSIN(ID,IS)
                   ENDIF
                 ELSE IF (IMETHOD == 4) THEN
                   GTEMP1 = MAX((1.-DT4A*FL(IP,ID,IS)),1.)
                   GTEMP2 = DT4A*SL(IP,ID,IS)/GTEMP1
                   FLHAB  = ABS(GTEMP2)
                   DELFL  = COFRM4(IS)*DT4A
                   USFM   = USNEW(IP)*MAX(FMEANWS(IP),FMEAN(IP))
                   TEMP   = USFM*DELFL
                   FLHAB  = MIN(FLHAB,TEMP)
                   IMATRAA(IS,ID,IP) = SIGN(FLHAB,GTEMP2)/DT4A*JAC
                   LIMFAC            = MIN(ONE,ABS(SIGN(FLHAB,GTEMP2/DT4A))/MAX(SMALL,ABS(IMATRAA(IS,ID,IP))))
                   IMATDAA(IS,ID,IP) = ZERO ! -FL(IP,ID,IS)
                 ENDIF 
               ENDDO ! ID
             ENDDO ! IS
             IF (LNANINFCHK) THEN
               IF (SUM(IMATRAA(:,:,IP)) .NE. SUM(IMATRAA(:,:,IP))) THEN
                 WRITE(12,*) 'NAN AT NODE', IP, 'AT DEPTH', DEP(IP)
                 CALL WWM_ABORT('NAN IN SOURCES 1a') 
               ENDIF
               IF (SUM(IMATDAA(:,:,IP)) .NE. SUM(IMATDAA(:,:,IP))) THEN
                 WRITE(12,*) 'NAN AT NODE', IP, 'AT DEPTH', DEP(IP)
                 CALL WWM_ABORT('NAN IN SOURCES 1b') 
               ENDIF
             ENDIF
             ACLOC = AC2(:,:,IP)
             IF (.NOT. LINID) THEN
               CALL SET_WIND( IP, WIND10, WINDTH )
               CALL SET_FRICTION( IP, ACLOC, WIND10, WINDTH, FPM )
               CALL SIN_LIN_CAV(IP,WIND10,WINDTH,FPM,IMATRA,SSINL)
               IMATRAA(:,:,IP) = IMATRAA(:,:,IP) + SSINL
             ENDIF
             IF (ISHALLOW(IP) .EQ. 1) THEN
               CALL MEAN_WAVE_PARAMETER(IP,ACLOC,HS,ETOT,SME01,SME10,KME01,KMWAM,KMWAM2)
               IF (MESTR .GT. 0) THEN
                 CALL TRIAD_ELDEBERKY(IP, HS, SME01, ACLOC, IMATRA, IMATDA, SSNL3, DSSNL3)
                 DO IS = 1, MSC
                   DO ID = 1, MDC
                     NEWDAC = SSNL3(IS,ID)*DT4A/MAX((1.-DT4A*DSSNL3(IS,ID)),1.)
                     MAXDAC = 0.0081*LIMFAK/(TWO*SPSIG(IS)*WK(IS,IP)**3*CG(IS,IP))*100
                     LIMFAC = ONE/MAX(ONE,NEWDAC/MAXDAC)
                     SC = SIGN(MIN(ABS(NEWDAC),MAXDAC),NEWDAC)/DT4A
                     !SSNL3(IS,ID)  = SC
                     !DSSNL3(IS,ID) = DSSNL3(IS,ID)*LIMFAC
                     !IF (ABS(SC) .GT. THR) WRITE(*,'(2I10,5F20.8)') IS, ID, NEWDAC, MAXDAC, DSSNL3(IS,ID), LIMFAC
                   END DO
                 END DO
               ENDIF ! MESTR
               IF (MESBR .GT. 0) CALL SDS_SWB(IP, SME01, KMWAM, ETOT, HS, ACLOC, IMATRA, IMATDA, SSBR, DSSBR)
               IF (MESBF .GT. 0) CALL SDS_BOTF(IP,ACLOC,IMATRA,IMATDA,SSBF,DSSBF)
               IMATDAA(:,:,IP) = IMATDAA(:,:,IP) + DSSBR  + DSSNL3 + DSSBF + DSSVEG 
               IMATRAA(:,:,IP) = IMATRAA(:,:,IP) + SSBR + SSNL3 + SSVEG 
             ENDIF ! ISHALLOW(IP) .EQ. 1
           ELSE ! IOBP(IP) .NE. 0
             IF (LSOUBOUND) THEN ! Source terms on boundary ...
               IF (IOBP(IP) .NE. 2) THEN
                 DO ID = 1, MDC
                   DO IS = 1, MSC
                     IF (AC2(IS,ID,IP) .LT. THR) CYCLE
                     JAC = ONE/PI2/SPSIG(IS)
                     !IMATDAA(IS,ID,IP) =  FL(IP,ID,IS) !... this is not working right, reason is unknown, there signchanges that should not be
                     !IMATRAA(IS,ID,IP) =  SL(IP,ID,IS)/PI2/SPSIG(IS) 
                     FF = FL3(IP,ID,IS)
                     !WRITE(11140,'(2I10,7E20.10)') IS, ID, FF, SSDS(ID,IS)/FF, DSSDS(ID,IS), SSIN(ID,IS)/FF, DSSNL4(ID,IS), SSNL4(ID,IS)/FF, DSSNL4(ID,IS) 
                     IF (IMETHOD == 0) THEN
                       IMATRAA(IS,ID,IP) =  SL(IP,ID,IS)/PI2/SPSIG(IS)
                       IMATDAA(IS,ID,IP) = ZERO
                     ELSE IF (IMETHOD == 1) THEN
                       IMATRAA(IS,ID,IP) = (SSIN(ID,IS)+SSDS(ID,IS)+SSNL4(ID,IS))*JAC
                     ELSE IF (IMETHOD == 2) THEN
                       IMATRAA(IS,ID,IP) = (SSIN(ID,IS)+SSNL4(ID,IS))*JAC
                       IMATDAA(IS,ID,IP) = -TWO*DSSDS(ID,IS)
                     ELSE IF (IMETHOD == 3) THEN
                       IF (SSIN(ID,IS) .GT. ZERO) THEN
                         IMATRAA(IS,ID,IP) = SSIN(ID,IS)*JAC
                       ELSE
                         IMATDAA(IS,ID,IP) = -TWO*DSSIN(ID,IS)
                       ENDIF
                     ELSE IF (IMETHOD == 4) THEN
                       GTEMP1 = MAX((1.-DT4A*FL(IP,ID,IS)),1.)
                       GTEMP2 = DT4A*SL(IP,ID,IS)/GTEMP1
                       FLHAB  = ABS(GTEMP2)
                       DELFL  = COFRM4(IS)*DT4A
                       USFM   = USNEW(IP)*MAX(FMEANWS(IP),FMEAN(IP))
                       TEMP   = USFM*DELFL
                       FLHAB  = MIN(FLHAB,TEMP)
                       IMATRAA(IS,ID,IP) = SIGN(FLHAB,GTEMP2)/DT4A*JAC
                       LIMFAC            = MIN(ONE,ABS(SIGN(FLHAB,GTEMP2/DT4A))/MAX(THR,ABS(IMATRAA(IS,ID,IP))))
                       IMATDAA(IS,ID,IP) = ZERO ! -FL(IP,ID,IS)
                     ENDIF
                   ENDDO ! ID
                 ENDDO ! IS
                 ACLOC = AC2(:,:,IP)
                 IF (.NOT. LINID) THEN
                   CALL SET_WIND( IP, WIND10, WINDTH )
                   CALL SET_FRICTION( IP, ACLOC, WIND10, WINDTH, FPM )
                   CALL SIN_LIN_CAV(IP,WIND10,WINDTH,FPM,IMATRA,SSINL)
                   IMATRAA(:,:,IP) = IMATRAA(:,:,IP) + SSINL
                 ENDIF
                 IF (ISHALLOW(IP) .EQ. 1) THEN
                   CALL MEAN_WAVE_PARAMETER(IP,ACLOC,HS,ETOT,SME01,SME10,KME01,KMWAM,KMWAM2)
                   IF (MESTR .GT. 0) THEN
                     CALL TRIAD_ELDEBERKY(IP, HS, SME01, ACLOC, IMATRA, IMATDA, SSNL3, DSSNL3)
                     DO IS = 1, MSC
                       DO ID = 1, MDC
                         NEWDAC = SSNL3(IS,ID)*DT4A/MAX((1.-DT4A*DSSNL3(IS,ID)),1.)
                         MAXDAC = 0.0081*LIMFAK/(TWO*SPSIG(IS)*WK(IS,IP)**3*CG(IS,IP))*100
                         LIMFAC = ONE/MAX(ONE,NEWDAC/MAXDAC)
                         SC = SIGN(MIN(ABS(NEWDAC),MAXDAC),NEWDAC)/DT4A
                         !SSNL3(IS,ID)  = SC
                         !DSSNL3(IS,ID) = DSSNL3(IS,ID)*LIMFAC
                         !IF (ABS(SC) .GT. THR) WRITE(*,'(2I10,5F20.8)') IS, ID, NEWDAC, MAXDAC, DSSNL3(IS,ID), LIMFAC
                       END DO
                     END DO
                   ENDIF ! MESTR
                   IF (MESBR .GT. 0) CALL SDS_SWB(IP, SME01, KMWAM, ETOT, HS, ACLOC, IMATRA, IMATDA, SSBR, DSSBR)
                   IF (MESBF .GT. 0) CALL SDS_BOTF(IP,ACLOC,IMATRA,IMATDA,SSBF,DSSBF)
                   IMATDAA(:,:,IP) = IMATDAA(:,:,IP) + DSSBR  + DSSNL3 + DSSBF + DSSVEG
                   IMATRAA(:,:,IP) = IMATRAA(:,:,IP) + SSBR + SSNL3 + SSVEG
                 ENDIF ! ISHALLOW(IP) .EQ. 1
               ENDIF
             ENDIF
           ENDIF
           !IF (IP==TESTNODE) WRITE(*,*) IP, SUM(IMATRAA(:,:,IP)), DEP(IP), IOBP(IP), SUM(SSBR), SUM(SSBF), SUM(SSINL), SUM(SSNL3)
         ENDDO
!$OMP END PARALLEL
         !WRITE(*,'(A20,3F15.10)') 'FROM SPECINT', SUM(IMATRAA), SUM(IMATDAA), SUM(AC2)
      END SUBROUTINE
!**********************************************************************
!*                                                                    *
!**********************************************************************
      SUBROUTINE SOURCE_INT_IMP_WAM_POST
         USE DATAPOOL
         IMPLICIT NONE
         INTEGER        :: IP, IS, ID
         REAL(rkind)    :: ACLOC(MSC,MDC), VEC2RAD
         REAL(rkind)    :: IMATRA(MSC,MDC), IMATDA(MSC,MDC)

         IMATDAA = 0.
         IMATRAA = 0.

         DO IP = 1, MNP
           IF ((ABS(IOBP(IP)) .NE. 1 .AND. IOBP(IP) .NE. 3)) THEN
             IF ( DEP(IP) .GT. DMIN .AND. IOBP(IP) .NE. 2) THEN
               THWOLD(IP,1) = THWNEW(IP)
               THWNEW(IP) = VEC2RAD(WINDXY(IP,1),WINDXY(IP,2))
               U10NEW(IP) = MAX(TWO,SQRT(WINDXY(IP,1)**2+WINDXY(IP,2)**2)) * WINDFAC
               Z0NEW(IP) = Z0OLD(IP,1)
               DO IS = 1, MSC
                 DO ID = 1, MDC
                   FL3(IP,ID,IS) =  AC2(IS,ID,IP) * PI2 * SPSIG(IS)
                   FL(IP,ID,IS)  =  IMATDAA(IS,ID,IP)
                   SL(IP,ID,IS)  =  IMATRAA(IS,ID,IP) * PI2 * SPSIG(IS)
                 END DO
               END DO
               CALL POSTINTRHS (FL3(IP,:,:), FL(IP,:,:), IP, IP, 1, &
     &                          THWOLD(IP,1), USOLD(IP,1), &
     &                          TAUW(IP), Z0OLD(IP,1), &
     &                          ROAIRO(IP,1), ZIDLOLD(IP,1), &
     &                          U10NEW(IP), THWNEW(IP), USNEW(IP), &
     &                          Z0NEW(IP), ROAIRN(IP), ZIDLNEW(IP), &
     &                          SL(IP,:,:), FCONST(IP,:), FMEANWS(IP), MIJ(IP))
               DO IS = 1, MSC
                 DO ID = 1, MDC
                   AC2(IS,ID,IP) = FL3(IP,ID,IS) / PI2 / SPSIG(IS)
                 END DO
               END DO
             END IF !
           ELSE
             IF (LSOUBOUND) THEN ! Source terms on boundary ...
               IF ( DEP(IP) .GT. DMIN .AND. IOBP(IP) .NE. 2) THEN
                 FL = FL3
                 THWOLD(:,1) = THWNEW
                 U10NEW = MAX(TWO,SQRT(WINDXY(:,1)**2+WINDXY(:,2)**2)) * WINDFAC
                 Z0NEW(IP) = Z0OLD(IP,1)
                 THWNEW(IP) = VEC2RAD(WINDXY(IP,1),WINDXY(IP,2))
                 DO IS = 1, MSC
                   DO ID = 1, MDC
                     FL3(IP,ID,IS) =  AC2(IS,ID,IP) * PI2 * SPSIG(IS)
                     FL(IP,ID,IS)  =  IMATDAA(IS,ID,IP)
                     SL(IP,ID,IS)  =  IMATRAA(IS,ID,IP) * PI2 * SPSIG(IS)
                   END DO
                 END DO
                 CALL POSTINTRHS (FL3(IP,:,:), FL(IP,:,:), IP, IP, 1, &
     &                            THWOLD(IP,1), USOLD(IP,1), &
     &                            TAUW(IP), Z0OLD(IP,1), &
     &                            ROAIRO(IP,1), ZIDLOLD(IP,1), &
     &                            U10NEW(IP), THWNEW(IP), USNEW(IP), &
     &                            Z0NEW(IP), ROAIRN(IP), ZIDLNEW(IP), &
     &                            SL(IP,:,:), FCONST(IP,:), FMEANWS(IP), MIJ(IP))
                 DO IS = 1, MSC
                   DO ID = 1, MDC
                     AC2(IS,ID,IP) = FL3(IP,ID,IS) / PI2 / SPSIG(IS)
                   END DO
                 END DO
               ENDIF
             ENDIF
           ENDIF
         END DO
      END SUBROUTINE
!**********************************************************************
!*                                                                    *
!**********************************************************************
      SUBROUTINE INT_IP_STAT(IP,DT,LIMITER,ACLOC,ISELECT)

         USE DATAPOOL
         IMPLICIT NONE

         INTEGER, INTENT(IN) :: IP, ISELECT
         LOGICAL, INTENT(IN) :: LIMITER
         REAL(rkind), INTENT(IN) :: DT
         REAL(rkind), INTENT(INOUT) :: ACLOC(MSC,MDC)
         REAL(rkind)                :: IMATRA(MSC,MDC), IMATDA(MSC,MDC)

         INTEGER :: IS, ID

         REAL(rkind)    :: NPF(MSC)
         REAL(rkind)    :: OLDAC
         REAL(rkind)    :: NEWDAC
         REAL(rkind)    :: MAXDAC, CONST, SND, USTAR

         CALL SOURCETERMS(IP, ACLOC, IMATRA, IMATDA, .FALSE., ISELECT, 'INT_IP_STAT')  ! 1. CALL

         CONST = PI2**2*3.0*1.0E-7*DT*SPSIG(MSC_HF(IP))
         SND   = PI2*5.6*1.0E-3

         IF (LIMITER) THEN
           IF (MELIM == 1) THEN
             NPF = 0.0081_rkind*LIMFAK/(TWO*SPSIG*WK(:,IP)**3*CG(:,IP))
           ELSE IF (MELIM == 2) THEN
             DO IS = 1, MSC
               USTAR = MAX(UFRIC(IP), G9*SND/SPSIG(IS))
               NPF(IS) = ABS((CONST*USTAR)/(SPSIG(IS)**3*WK(IS,IP)))
             END DO
           END IF
         END IF

!         if (SUM(ACLOC) .NE. SUM(ACLOC)) STOP 'NAN l. 174 wwm_specint.F90'

         DO IS = 1, MSC
           MAXDAC = NPF(IS)
           DO ID = 1, MDC
             OLDAC  = ACLOC(IS,ID)
             NEWDAC = IMATRA(IS,ID) * DT / ( 1.0 - DT * MIN(ZERO,IMATDA(IS,ID))) 
             IF (LIMITER) NEWDAC = SIGN(MIN(MAXDAC,ABS(NEWDAC)),NEWDAC)
             ACLOC(IS,ID) = MAX( 0.0_rkind, OLDAC + NEWDAC ) 
           END DO
         END DO

!         write(*,'(A10,2I10,L10,I10,2F15.6)') 'after', ip, iobp(ip), limiter, iselect, sum(acloc), sum(imatra)

         RETURN
      END SUBROUTINE
!**********************************************************************
!*                                                                    *
!**********************************************************************
      SUBROUTINE RKS_SP3(IP,DTSII,LIMITER,ACLOC,ISELECT)

         USE DATAPOOL
         IMPLICIT NONE

         INTEGER, INTENT(IN) :: IP, ISELECT
         LOGICAL, INTENT(IN) :: LIMITER
         REAL(rkind), INTENT(IN)    :: DTSII
         REAL(rkind), INTENT(INOUT) :: ACLOC(MSC,MDC)

         INTEGER :: IS, ID
         REAL(rkind)    :: ACOLD(MSC,MDC)
         REAL(rkind)    :: NPF(MSC)
         REAL(rkind)    :: NEWDAC, MAXDAC, CONST, SND, USTAR
         REAL(rkind)    :: IMATRA(MSC,MDC), IMATDA(MSC,MDC)

         CONST = PI2**2*3.0*1.0E-7*DTSII*SPSIG(MSC_HF(IP))
         SND   = PI2*5.6*1.0E-3

         IF (LIMITER) THEN
           IF (MELIM == 1) THEN
             NPF = 0.0081_rkind*LIMFAK/(TWO*SPSIG*WK(:,IP)**3*CG(:,IP))
           ELSE IF (MELIM == 2) THEN
             DO IS = 1, MSC
               USTAR = MAX(UFRIC(IP), G9*SND/SPSIG(IS))
               NPF(IS) = ABS((CONST*USTAR)/(SPSIG(IS)**3*WK(IS,IP)))
             END DO
           END IF
         END IF
       
         CALL SOURCETERMS(IP, ACLOC, IMATRA, IMATDA, .FALSE., ISELECT, 'RKS_SP3 1')  ! 1. CALL

         ACOLD = ACLOC

         DO IS = 1, MSC
           MAXDAC = NPF(IS)
           DO ID = 1, MDC
             NEWDAC = IMATRA(IS,ID) * DTSII / ( ONE - DTSII * MIN(ZERO,IMATDA(IS,ID)) )
             IF (LIMITER) NEWDAC = SIGN(MIN(MAXDAC,ABS(NEWDAC)),NEWDAC)
             ACLOC(IS,ID) = MAX( ZERO, ACLOC(IS,ID) + NEWDAC )
           END DO
         END DO

         !WRITE(*,*) '1 RK-TVD', SUM(ACOLD), SUM(ACLOC)

         CALL SOURCETERMS(IP, ACLOC, IMATRA, IMATDA, .FALSE., ISELECT, 'RKS_SP3 2') ! 2. CALL

         DO IS = 1, MSC
           MAXDAC = NPF(IS)
           DO ID = 1, MDC
             NEWDAC = IMATRA(IS,ID) * DTSII / ( ONE - DTSII * MIN(ZERO,IMATDA(IS,ID)) )
             IF (LIMITER) NEWDAC = SIGN(MIN(MAXDAC,ABS(NEWDAC)),NEWDAC) 
             ACLOC(IS,ID) = MAX( ZERO, 0.75_rkind * ACOLD(IS,ID) +  0.25_rkind * ACLOC(IS,ID) + 0.25_rkind * NEWDAC)
           END DO
         END DO

         !WRITE(*,*) '2 RK-TVD', SUM(ACOLD), SUM(ACLOC)

         CALL SOURCETERMS(IP, ACLOC, IMATRA, IMATDA, .FALSE., ISELECT, 'RKS_SP3 3') ! 3. CALL

         DO IS = 1, MSC
           MAXDAC = NPF(IS)
           DO ID = 1, MDC
             NEWDAC = IMATRA(IS,ID) * DTSII / ( ONE - DTSII * MIN(ZERO,IMATDA(IS,ID)) )
             IF (LIMITER) NEWDAC = SIGN(MIN(MAXDAC,ABS(NEWDAC)),NEWDAC)
             ACLOC(IS,ID) = MAX( ZERO, ONE/THREE * ACOLD(IS,ID) + TWO/THREE * ACLOC(IS,ID) + TWO/THREE * NEWDAC)
           END DO
         END DO

         !IF (IP == 1701) WRITE(*,*) SUM(IMATRA), SUM(IMATDA), SUM(ACLOC)

         !WRITE(*,*) '3 RK-TVD', SUM(ACOLD), SUM(ACLOC)

      END SUBROUTINE
!**********************************************************************
!*                                                                    *
!**********************************************************************
      SUBROUTINE INT_IP_DYN(IP, DT, LIMIT, DTMIN, ITRMX, ACLOC, ITER, ISELECT)

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

         INTEGER, INTENT(IN)        :: IP, ITRMX, ISELECT
         INTEGER, INTENT(OUT)       :: ITER
         LOGICAL,INTENT(IN)         :: LIMIT
         REAL(rkind), INTENT(IN)    :: DTMIN, DT
         REAL(rkind), INTENT(INOUT) :: ACLOC(MSC,MDC)
         REAL(rkind)                :: IMATRA(MSC,MDC), IMATDA(MSC,MDC)

         INTEGER :: IS, ID

         REAL(rkind)    :: ACOLD(MSC,MDC)
         REAL(rkind)    :: NPF(MSC)
         REAL(rkind)    :: AFILT
         REAL(rkind)    :: NEWDAC, CONST, SND, DAMAX, AFAC
         REAL(rkind)    :: MAXDAC, DTMAX, DTTOT, DTLEFT, DT4SI, USTAR

         REAL(rkind), PARAMETER :: MAXDTFAC = VERYLARGE 

         LOGICAL :: LSTABLE

         REAL(rkind) :: XPP, XRR, XFILT, XREL, XFLT, FACP, DAM(MSC)

#ifdef ST_DEF
         XPP    = 0.15
         XRR    = 0.10
         XFILT  = 0.0001 ! AR: check why it must be so small ..
         XPP    = MAX ( 1.E-6_rkind , XPP )
         XRR    = MAX ( 1.E-6_rkind , XRR )
         XREL   = XRR
         XFLT   = MAX ( ZERO , XFILT ) 
         FACP   = 2*XPP / PI2 * 0.62E-3 * PI2**4 / G9**2  ! s4/m4
         DAM    = FACP / ( SPSIG * WK(:,IP)**3 ) / CG(:,IP) ! s * m³ * s4/m4 = 
         AFILT  = MAX ( DAM(MSC) , XFLT*MAXVAL(ACLOC))!*PI2*SPSIG(MSC_HF(IP)) )
#endif

         CONST = PI2**2*3.0*1.0E-7*DT*SPSIG(MSC_HF(IP))
         SND   = PI2*5.6*1.0E-3
         IMATRA = ZERO 
         IMATDA = ZERO

         IF (LIMIT) THEN
           IF (MELIM == 1) THEN
             NPF = 0.0081_rkind*LIMFAK/(TWO*SPSIG*WK(:,IP)**3*CG(:,IP))
           ELSE IF (MELIM == 2) THEN
             DO IS = 1, MSC
               USTAR = MAX(UFRIC(IP), G9*SND/SPSIG(IS))
               NPF(IS) = ABS((CONST*USTAR)/(SPSIG(IS)**3*WK(IS,IP)))
             END DO
           END IF
         END IF

         DTTOT = 0.
         ITER  = 0

         DO WHILE ( DTTOT < DT )

           ACOLD = ACLOC
           IF (ITER == 0) DT4SI = DT
           IF (ITER == 0) DTLEFT = DT
           ITER  = ITER + 1
           DTMAX =  DT

           CALL SOURCETERMS(IP, ACLOC, IMATRA, IMATDA, .FALSE.,ISELECT, 'INT_IP_DYN')  ! 1. CALL
        
           DO ID = 1, MDC
             DO IS = 1, MSC_HF(IP)
               IF (ABS(IMATRA(IS,ID)) .GT. VERYSMALL) THEN                
                 DAMAX  = MIN ( DAM(IS) , MAX ( XREL*ACLOC(IS,ID), AFILT ) )
                 AFAC  = ONE / MAX( 1.E-10_rkind , ABS(IMATRA(IS,ID)/DAMAX) )
                 DTMAX = MIN (DTMAX ,AFAC/(MAX(1.E-10_rkind, ONE + AFAC*MIN(ZERO,IMATDA(IS,ID))))) 
               END IF
             END DO
           END DO

           DT4SI  = DTMAX
           DT4SI  = MAX(DTMIN,DT4SI) ! This is not entirely stable !!!
           DTLEFT = DT - DTTOT

           IF ( DTLEFT > THR .AND. DTLEFT < DT4SI) THEN
             DT4SI = (DT - DTTOT)
           ELSE IF ( DTLEFT .GE. DT4SI .AND. ITER .EQ. ITRMX) THEN
             DT4SI = DTLEFT 
             LSTABLE = .FALSE. 
           END IF 

           DTTOT = DTTOT + DT4SI

           IF ( DT4SI .LT. DTMAX + SMALL) THEN
             LSTABLE = .TRUE. 
           ELSE
             LSTABLE = .FALSE.
           END IF

           IF (LSTABLE) THEN
             DO IS = 1, MSC_HF(IP)
               DO ID = 1, MDC
                 NEWDAC = IMATRA(IS,ID) * DT4SI / ( 1.0 - DT4SI * MIN(ZERO,IMATDA(IS,ID)) )
                 ACLOC(IS,ID) = MAX( ZERO, ACLOC(IS,ID) + NEWDAC )
               END DO
             END DO
             CALL SOURCETERMS(IP,ACLOC,IMATRA,IMATDA,.FALSE.,ISELECT, 'INT_IP_DYN - STABLE1') ! 2. CALL
             DO IS = 1, MSC_HF(IP)
               DO ID = 1, MDC
                 NEWDAC = IMATRA(IS,ID) * DT4SI / ( 1.0 - DT4SI * MIN(ZERO,IMATDA(IS,ID)) )
                 ACLOC(IS,ID) = MAX( ZERO, 3._rkind/4._rkind * ACOLD(IS,ID) + 1._rkind/4._rkind * ACLOC(IS,ID) + 1._rkind/4._rkind * NEWDAC)
               END DO
             END DO
             CALL SOURCETERMS(IP,ACLOC,IMATRA,IMATDA, .FALSE.,ISELECT, 'INT_IP_DYN STABLE2') ! 3. CALL
             DO IS = 1, MSC_HF(IP)
               MAXDAC = NPF(IS)
               DO ID = 1, MDC
                 NEWDAC = IMATRA(IS,ID) * DT4SI / ( 1.0 - DT4SI * MIN(ZERO,IMATDA(IS,ID)) )
                 ACLOC(IS,ID) = MAX( ZERO,  1._rkind/3._rkind * ACOLD(IS,ID) + 2._rkind/3._rkind * ACLOC(IS,ID) + 2._rkind/3._rkind * NEWDAC)
               END DO
             END DO
           ELSE ! .NOT. LSTABLE
             DO IS = 1, MSC_HF(IP)
               MAXDAC = NPF(IS)
               DO ID = 1, MDC
                 NEWDAC = IMATRA(IS,ID) * DT4SI / ( 1.0 - DT4SI * MIN(ZERO,IMATDA(IS,ID)) )
                 IF (LIMIT) NEWDAC = SIGN(MIN(MAXDAC,ABS(NEWDAC)),NEWDAC)
                 ACLOC(IS,ID) = MAX( ZERO, ACLOC(IS,ID) + NEWDAC )
               END DO
             END DO
             CALL SOURCETERMS(IP, ACLOC, IMATRA, IMATDA, .FALSE., ISELECT, 'INT_IP_DYN - LIM1 ') ! 2. CALL
             DO IS = 1, MSC_HF(IP)
               MAXDAC = NPF(IS)
               DO ID = 1, MDC
                 NEWDAC = IMATRA(IS,ID) * DT4SI / ( 1.0 - DT4SI * MIN(ZERO,IMATDA(IS,ID)) )
                 IF (LIMIT) NEWDAC = SIGN(MIN(MAXDAC,ABS(NEWDAC)),NEWDAC)
                 ACLOC(IS,ID) = MAX( ZERO, 3./4. * ACOLD(IS,ID) +  1./4. * ACLOC(IS,ID) + 1./4. * NEWDAC)
               END DO
             END DO
             CALL SOURCETERMS(IP, ACLOC, IMATRA, IMATDA, .FALSE., ISELECT, 'INT_IP_DYN - LIM2 ') ! 3. CALL
             DO IS = 1, MSC_HF(IP)
               MAXDAC = NPF(IS)
               DO ID = 1, MDC
                 NEWDAC = IMATRA(IS,ID) * DT4SI / ( 1.0 - DT4SI * MIN(ZERO,IMATDA(IS,ID)) )
                 IF (LIMIT) NEWDAC = SIGN(MIN(MAXDAC,ABS(NEWDAC)),NEWDAC)
                 ACLOC(IS,ID) = MAX( ZERO, 1./3. * ACOLD(IS,ID) +  2./3. * ACLOC(IS,ID) + 2./3. * NEWDAC)
               END DO
             END DO
           END IF
         END DO
         END SUBROUTINE
!**********************************************************************
!*                                                                    *
!**********************************************************************
         SUBROUTINE ACTION_LIMITER
         USE DATAPOOL
         IMPLICIT NONE

         INTEGER                 :: IP, IS, ID
         REAL(rkind)             :: NEWDAC, OLDAC, NEWAC, DELT, XIMP, DELFL(MSC)
         REAL(rkind)             :: MAXDAC, CONST, SND, UFR_LIM, DELT5, USFM
         REAL(rkind)             :: ACLOC(MSC,MDC), ACOLD(MSC,MDC), MAXDACOLD


         CONST = PI2**2*3.0*1.0E-7*DT4S*SPSIG(MSC)
         SND   = PI2*5.6*1.0E-3

         DELT = DT4S
         XIMP = 1._rkind
         DELT5 = XIMP*DELT
         DELFL= COFRM4*DELT
         MAXDAC = ZERO

!$OMP PARALLEL DO SCHEDULE(DYNAMIC,1) PRIVATE(IP,ACOLD,ACLOC, MAXDAC, UFR_LIM, NEWAC, OLDAC, NEWDAC)
         DO IP = 1, MNP
           IF (DEP(IP) .LT. DMIN .OR. IOBP(IP) .EQ. 2) CYCLE
           ACOLD = AC1(:,:,IP)
           ACLOC = AC2(:,:,IP)
           DO IS = 1, MSC
             IF (MELIM .EQ. 1) THEN
               MAXDAC = 0.0081*LIMFAK/(TWO*SPSIG(IS)*WK(IS,IP)**3*CG(IS,IP))
             ELSE IF (MELIM .EQ. 2) THEN
               UFR_LIM = MAX(UFRIC(IP),G9*SND/SPSIG(IS))
               MAXDAC  = LIMFAK*ABS((CONST*UFR_LIM)/(SPSIG(IS)**3*WK(IS,IP)))
             ELSE IF (MELIM .EQ. 3) THEN
               IF (USNEW(IP) .GT. SMALL) THEN
                 MAXDAC = COFRM4(IS)*USNEW(IP)*MAX(FMEANWS(IP),FMEAN(IP))/PI2/SPSIG(IS)*DT4A
               ELSE
                 MAXDAC = 0.0081*LIMFAK/(TWO*SPSIG(IS)*WK(IS,IP)**3*CG(IS,IP))
               ENDIF
             END IF
             DO ID = 1, MDC
               NEWAC  = ACLOC(IS,ID)
               OLDAC  = ACOLD(IS,ID)
               NEWDAC = NEWAC - OLDAC
!               IF (NEWDAC .GT. 0.) THEN
                 NEWDAC = SIGN(MIN(MAXDAC,ABS(NEWDAC)), NEWDAC)
!               ELSE
!                 IF (QBLOCAL(IP) .LT. THR) NEWDAC = SIGN(MIN(MAXDAC,ABS(NEWDAC)), NEWDAC)
!               END IF
               AC2(IS,ID,IP) = MAX( zero, OLDAC + NEWDAC )
               !IF (MELIM .EQ. 3) FL3(IP,ID,IS) = AC2(IP,IS,ID) * PI2 * SPSIG(IS)
             END DO
           END DO
         END DO
         END SUBROUTINE
!**********************************************************************
!*                                                                    *
!**********************************************************************
         SUBROUTINE ACTION_LIMITER_LOCAL(IP,ACLOC,ACOLD)
         USE DATAPOOL
         IMPLICIT NONE

         INTEGER, INTENT(IN)        :: IP
         INTEGER                    :: IS, ID
         REAL(rkind), INTENT(INOUT) :: ACLOC(MSC,MDC)
         REAL(rkind), INTENT(IN)    :: ACOLD(MSC,MDC)
         REAL(rkind)                :: NEWDAC, OLDAC, NEWAC, DELT, XIMP, DELFL(MSC)
         REAL(rkind)                :: MAXDAC, CONST, SND, UFR_LIM, DELT5, USFM
         REAL(rkind)                :: MAXDACOLD

         CONST = PI2**2*3.0*1.0E-7*DT4S*SPSIG(MSC)
         SND   = PI2*5.6*1.0E-3

         DELT = DT4S
         XIMP = 1._rkind
         DELT5 = XIMP*DELT
         DELFL= COFRM4*DELT
         MAXDAC = ZERO

         DO IS = 1, MSC

           IF (MELIM .EQ. 1) THEN
             MAXDAC = 0.0081*LIMFAK/(TWO*SPSIG(IS)*WK(IS,IP)**3*CG(IS,IP))
           ELSE IF (MELIM .EQ. 2) THEN
             UFR_LIM = MAX(UFRIC(IP),G9*SND/SPSIG(IS))
             MAXDAC  = LIMFAK*ABS((CONST*UFR_LIM)/(SPSIG(IS)**3*WK(IS,IP)))
           ELSE IF (MELIM .EQ. 3) THEN
             IF (USNEW(IP) .GT. SMALL) THEN
               MAXDAC = COFRM4(IS)*USNEW(IP)*MAX(FMEANWS(IP),FMEAN(IP))/PI2/SPSIG(IS)*DT4A
             ELSE
               MAXDAC = 0.0081*LIMFAK/(TWO*SPSIG(IS)*WK(IS,IP)**3*CG(IS,IP))
             ENDIF
           END IF

           DO ID = 1, MDC
             NEWAC  = ACLOC(IS,ID)
             OLDAC  = ACOLD(IS,ID)
             NEWDAC = NEWAC - OLDAC
             NEWDAC = SIGN(MIN(MAXDAC,ABS(NEWDAC)), NEWDAC)
             ACLOC(IS,ID) = MAX( zero, OLDAC + NEWDAC )
           END DO

         END DO

         END SUBROUTINE
!**********************************************************************
!*                                                                    *
!**********************************************************************
         SUBROUTINE BREAK_LIMIT(IP,ACLOC,SSBRL)
         USE DATAPOOL
         IMPLICIT NONE

         INTEGER, INTENT(IN)  :: IP

         REAL(rkind), INTENT(INOUT)  :: ACLOC(MSC,MDC)
         REAL(rkind), INTENT(OUT)    :: SSBRL(MSC,MDC)

         REAL(rkind)                 :: EFTAIL, HS
         REAL(rkind)                 :: EMAX, RATIO, ETOT
         REAL(rkind)                 :: DINTSPEC

         ETOT   = 0.0
         EFTAIL = 1.0 / (PTAIL(1)-1.0)

         ETOT = DINTSPEC(IP,ACLOC)

         HS = 4.*SQRT(ETOT)

         !IF (LMONO_IN) HMAX(IP) = HMAX(IP) * SQRT(2.)

         EMAX = 1./16. * (HMAX(IP))**2

         IF (ETOT .GT. EMAX .AND. ETOT .GT. THR) THEN
           RATIO = EMAX/ETOT
           SSBRL = ACLOC - RATIO * ACLOC
           ACLOC = RATIO * ACLOC
         END IF
         END SUBROUTINE
!**********************************************************************
!*                                                                    *
!**********************************************************************
         SUBROUTINE BREAK_LIMIT_ALL
         USE DATAPOOL
         IMPLICIT NONE

         INTEGER              :: IP
         REAL(rkind)          :: HS
         REAL(rkind)          :: EMAX, RATIO, ETOT
         REAL(rkind)          :: DINTSPEC
         REAL(rkind)          :: ACLOC(MSC, MDC)
!      Print *, 'Passing BREAK_LIMIT_ALL'
         DO IP = 1, MNP
           ACLOC = AC2(:,:,IP)
           IF (ISHALLOW(IP) .EQ. 0) CYCLE
           ETOT = DINTSPEC(IP,ACLOC)
           HS = 4.*SQRT(ETOT)
           EMAX = 1./16. * (HMAX(IP))**2
!        WRITE(300,*) 'IP=', IP, ' HMAX=', HMAX(IP), ' DEP=', DEP(IP)
!        WRITE(300,*) '   ', IP, ' EMAX=', EMAX, ' ETOT=', ETOT
!        WRITE(300,*) '   ', IP, ' HS=', HS, ' BRHD=', BRHD

           IF (ETOT .GT. EMAX) THEN
             if(myrank==0) WRITE(300,*) '   break XP=', XP(IP)
             RATIO = EMAX/ETOT
             AC2(:,:,IP) = RATIO * ACLOC(:,:)
             AC1(:,:,IP) = RATIO * ACLOC(:,:)
           END IF
         END DO
         END SUBROUTINE
!**********************************************************************
!*                                                                    *
!**********************************************************************
         SUBROUTINE RESCALE_SPECTRUM()
         USE DATAPOOL
         IMPLICIT NONE
         INTEGER :: IP, IS, ID
         REAL(rkind)    :: ETOTAL, EPOSIT
         REAL(rkind)    :: FACTOR
         LOGICAL :: ENEG

         DO IP = 1, MNP
            DO IS = 1, MSC
               ETOTAL = ZERO 
               EPOSIT = ZERO 
               ENEG   = .FALSE.
               DO ID = 1, MDC
                  ETOTAL = ETOTAL + AC2(IS,ID,IP)
                  IF (AC2(IS,ID,IP) > ZERO) THEN
                     EPOSIT = EPOSIT + AC2(IS,ID,IP)
                  ELSE
                     ENEG = .TRUE.
                  END IF
               END DO
               IF (ENEG) THEN
                  IF (EPOSIT .GT. VERYSMALL) THEN
                    FACTOR = ETOTAL/EPOSIT
                  ELSE 
                    FACTOR = 0.
                  END IF
                  DO ID = 1, MDC
                     IF (AC2(IS,ID,IP) < ZERO) AC2(IS,ID,IP) = ZERO 
                     IF (FACTOR >= ZERO)  AC2(IS,ID,IP) = AC2(IS,ID,IP)*FACTOR
                     AC2(IS,ID,IP) = MAX(zero,AC2(IS,ID,IP))
                  END DO
               END IF
            END DO
         END DO
         END SUBROUTINE
!**********************************************************************
!*                                                                    *
!**********************************************************************
         SUBROUTINE SETSHALLOW
         USE DATAPOOL
         IMPLICIT NONE
         INTEGER :: IP

         DO IP = 1, MNP
           IF (WK(1,IP)*DEP(IP) .LT. PI) THEN
             ISHALLOW(IP) = 1
           ELSE
             ISHALLOW(IP) = 0
           END IF
         END DO
         END SUBROUTINE
!**********************************************************************
!*                                                                    *
!**********************************************************************