WW3 2dm File Format

model/src/CMakeLists.txt

@@ -47,16 +47,15 @@ target_include_directories(ww3_lib
 
 # Set compiler flags.
 if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$")
-  set(compile_flags -no-fma -ip -g -traceback -i4 -real-size 32 -fp-model precise
-    -assume byterecl -fno-alias -fno-fnalias)
+  set(compile_flags)
   # -sox only works on Linux
   if(LINUX)
     list(APPEND compile_flags -sox)
   endif()
-  set(compile_flags_release -O3)
+  set(compile_flags_release)
   # SHELL: prefix fixes CMake attempting to de-duplicate the repeated uses of 'all' in -warn, -debug, -check
   # See https://cmake.org/cmake/help/latest/command/target_compile_options.html#option-de-duplication
-  set(compile_flags_debug -O0 "SHELL:-debug all" "SHELL:-warn all" "SHELL:-check all" -check noarg_temp_created -fp-stack-check -heap-arrays -traceback -fpe0)
+  set(compile_flags_debug)
 
   if(APPLE)
     # The linker on macOS does not include `common symbols` (usually module variables without a default value) by default

model/src/w3gdatmd.F90

@@ -1043,6 +1043,7 @@ MODULE W3GDATMD
     LOGICAL :: FSFREQSHIFT  = .FALSE.
     LOGICAL :: FSSOURCE     = .FALSE.
     LOGICAL :: FSBCCFL      = .FALSE.
+    LOGICAL :: L2DM         = .TRUE. 
     LOGICAL :: DO_CHANGE_WLV
     REAL(8) :: SOLVERTHR_STP
     REAL(8) :: CRIT_DEP_STP
@@ -1403,6 +1404,7 @@ MODULE W3GDATMD
   LOGICAL, POINTER  :: FSN,FSPSI,FSFCT,FSNIMP,FSTOTALIMP,FSTOTALEXP
   LOGICAL, POINTER  :: FSREFRACTION, FSFREQSHIFT, FSSOURCE, FSBCCFL
   LOGICAL, POINTER  :: DO_CHANGE_WLV
+  LOGICAL, POINTER  :: L2DM 
   REAL(8), POINTER  :: SOLVERTHR_STP
   REAL(8), POINTER  :: CRIT_DEP_STP
   LOGICAL, POINTER :: B_JGS_TERMINATE_MAXITER
@@ -2840,6 +2842,7 @@ SUBROUTINE W3SETG ( IMOD, NDSE, NDST )
     FSREFRACTION => MPARS(IMOD)%SCHMS%FSREFRACTION
     FSFREQSHIFT => MPARS(IMOD)%SCHMS%FSFREQSHIFT
     FSSOURCE => MPARS(IMOD)%SCHMS%FSSOURCE
+    L2DM => MPARS(IMOD)%SCHMS%L2DM
     DO_CHANGE_WLV => MPARS(IMOD)%SCHMS%DO_CHANGE_WLV
     SOLVERTHR_STP => MPARS(IMOD)%SCHMS%SOLVERTHR_STP
     CRIT_DEP_STP => MPARS(IMOD)%SCHMS%CRIT_DEP_STP

model/src/w3gridmd.F90

@@ -926,6 +926,7 @@ MODULE W3GRIDMD
   LOGICAL :: IMPREFRACTION
   LOGICAL :: IMPFREQSHIFT
   LOGICAL :: IMPSOURCE
+  LOGICAL :: LREAD2DM
   LOGICAL :: SETUP_APPLY_WLV
   INTEGER :: JGS_MAXITER
   INTEGER :: nbSel
@@ -1087,7 +1088,7 @@ MODULE W3GRIDMD
        UGBCCFL, EXPFSN, EXPFSPSI, EXPFSFCT,       &
        IMPFSN, IMPTOTAL, EXPTOTAL,                &
        IMPREFRACTION, IMPFREQSHIFT,               &
-       IMPSOURCE,                                 &
+       IMPSOURCE, LREAD2DM,                       &
        JGS_TERMINATE_MAXITER,                     &
        JGS_TERMINATE_DIFFERENCE,                  &
        JGS_TERMINATE_NORM,                        &
@@ -2441,6 +2442,7 @@ SUBROUTINE W3GRID()
     IMPREFRACTION = .FALSE.
     IMPFREQSHIFT = .FALSE.
     IMPSOURCE = .FALSE.
+    LREAD2DM   = .FALSE. 
     SETUP_APPLY_WLV = .TRUE.
     SOLVERTHR_SETUP=1E-6
     CRIT_DEP_SETUP=0.1
@@ -2460,6 +2462,7 @@ SUBROUTINE W3GRID()
     ! read data from the unstructured devoted namelist
     CALL READNL ( NDSS, 'UNST', STATUS )
 
+    L2DM             = LREAD2DM
     B_JGS_USE_JACOBI = JGS_USE_JACOBI
     B_JGS_TERMINATE_MAXITER = JGS_TERMINATE_MAXITER
     B_JGS_TERMINATE_DIFFERENCE = JGS_TERMINATE_DIFFERENCE
@@ -3318,10 +3321,10 @@ SUBROUTINE W3GRID()
       WRITE (NDSO,2953) CFLTM, WDTHCG, WDTHTH
 #endif
       !
-      WRITE (NDSO,2956) UGBCCFL, UGOBCAUTO, UGOBCDEPTH,TRIM(UGOBCFILE), &
+      WRITE (NDSO,*) UGBCCFL, UGOBCAUTO, UGOBCDEPTH,TRIM(UGOBCFILE), &
            EXPFSN, EXPFSPSI, EXPFSFCT, IMPFSN, EXPTOTAL,&
            IMPTOTAL, IMPREFRACTION, IMPFREQSHIFT,      &
-           IMPSOURCE, SETUP_APPLY_WLV,                 &
+           IMPSOURCE, LREAD2DM, SETUP_APPLY_WLV,        &
            JGS_TERMINATE_MAXITER,                      &
            JGS_TERMINATE_DIFFERENCE,                   &
            JGS_TERMINATE_NORM,                         &
@@ -4000,7 +4003,11 @@ SUBROUTINE W3GRID()
       !
       ! Reading depths on unstructured grid (this also sets number of mesh points, NX)
       !
-      CALL READMSH(NDSG,FNAME)
+      IF (L2DM) THEN
+        CALL READ2DM(NDSG,FNAME)
+      ELSE
+        CALL READMSH(NDSG,FNAME)
+      ENDIF
       ALLOCATE(ZBIN(NX, NY),OBSX(NX,NY),OBSY(NX,NY))
       ZBIN(:,1) = VSC * ZB(:)
       !
@@ -4440,6 +4447,9 @@ SUBROUTINE W3GRID()
          CALL READMSHOBC(NDSG,UGOBCFILE,TMPSTA,UGOBCOK)
     IF ((GTYPE.EQ.UNGTYPE).AND.UGOBCAUTO.AND.(.NOT.UGOBCOK))  &
          CALL UG_GETOPENBOUNDARY(TMPSTA,ZBIN,UGOBCDEPTH)
+    IF ((GTYPE.EQ.UNGTYPE).AND.L2DM.AND.(.NOT.UGOBCOK))  &
+         CALL READ2DM_TMPSTA(23956,trim("meshbnd.2dm"),TMPSTA)     
+
     !
     ! 8.b Determine where to get the data
     !

model/src/w3initmd.F90

@@ -1,3 +1,4 @@
+
 !> @file
 !> @brief Contains module W3INITMD.
 !>

model/src/w3profsmd_pdlib.F90

@@ -2854,12 +2854,24 @@ SUBROUTINE PDLIB_W3XYPUG_BLOCK_EXPLICIT(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC)
     !
     USE W3ODATMD, only: IAPROC
     USE W3GDATMD, only: B_JGS_USE_JACOBI
+    USE W3TIMEMD, only: DSEC21
+    USE W3ODATMD, only: TBPI0, TBPIN, FLBPI
+    USE W3WDATMD, only: TIME
 
     LOGICAL, INTENT(IN) :: LCALC
     INTEGER, INTENT(IN) :: IMOD
     REAL, INTENT(IN) :: FACX, FACY, DTG, VGX, VGY
+    REAL             :: RD1, RD2
 
-    CALL PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC)
+    IF ( FLBPI ) THEN
+      RD1  = DSEC21 ( TBPI0, TIME )
+      RD2  = DSEC21 ( TBPI0, TBPIN )
+    ELSE
+      RD1=1.
+      RD2=0.
+    END IF
+
+    CALL PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, RD1, RD2, DTG, VGX, VGY, LCALC)
     !/
     !/ End of W3XYPFSN ----------------------------------------------------- /
     !/
@@ -5546,26 +5558,26 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL
     INTEGER :: nbIter, ISPnextDir, ISPprevDir
     INTEGER :: ISPp1, ISPm1, JP, ICOUNT1, ICOUNT2
     ! for the exchange
-    REAL*8  :: CCOS, CSIN, CCURX, CCURY
-    REAL*8  :: eSum(NSPEC), FRLOCAL
-    REAL*8  :: eA_THE, eC_THE, eA_SIG, eC_SIG, eSI
-    REAL*8  :: CAD(NSPEC), CAS(NSPEC), ACLOC(NSPEC)
-    REAL*8  :: CP_SIG(NSPEC), CM_SIG(NSPEC)
-    REAL*8  :: eFactM1, eFactP1
-    REAL*8  :: Sum_Prev, Sum_New, p_is_converged, DiffNew, prop_conv
-    REAL*8  :: Sum_L2, Sum_L2_GL
+    REAL  :: CCOS, CSIN, CCURX, CCURY
+    REAL  :: eSum(NSPEC), FRLOCAL
+    REAL  :: eA_THE, eC_THE, eA_SIG, eC_SIG, eSI
+    REAL  :: CAD(NSPEC), CAS(NSPEC), ACLOC(NSPEC)
+    REAL  :: CP_SIG(NSPEC), CM_SIG(NSPEC)
+    REAL  :: eFactM1, eFactP1
+    REAL  :: Sum_Prev, Sum_New, p_is_converged, DiffNew, prop_conv
+    REAL  :: Sum_L2, Sum_L2_GL
     REAL  :: DMM(0:NK2), DAM(NSPEC), DAM2(NSPEC), SPEC(NSPEC)
-    REAL*8  :: eDiff(NSPEC), eProd(NSPEC), eDiffB(NSPEC)
-    REAL*8  :: DWNI_M2(NK), CWNB_M2(1-NTH:NSPEC)
+    REAL  :: eDiff(NSPEC), eProd(NSPEC), eDiffB(NSPEC)
+    REAL  :: DWNI_M2(NK), CWNB_M2(1-NTH:NSPEC)
     REAL  :: VAnew(NSPEC), VFLWN(1-NTH:NSPEC), JAC, JAC2
     REAL  :: VAAnew(1-NTH:NSPEC+NTH), VAAacloc(1-NTH:NSPEC+NTH)
     REAL  :: VAinput(NSPEC), VAacloc(NSPEC), ASPAR_DIAG(NSPEC)
     REAL  :: aspar_diag_local(nspec), aspar_off_diag_local(nspec), b_jac_local(nspec)
-    REAL*8 :: eDiffSing, eSumPart
+    REAL  :: eDiffSing, eSumPart
     REAL  :: EMEAN, FMEAN, FMEAN1, WNMEAN, AMAX, U10ABS, U10DIR, TAUA, TAUADIR
     REAL  :: USTAR, USTDIR, TAUWX, TAUWY, CD, Z0, CHARN, FMEANWS, DLWMEAN
-    REAL*8  :: eVal1, eVal2
-    REAL*8  :: eVA, eVO, CG2, NEWDAC, NEWAC, OLDAC, MAXDAC
+    REAL  :: eVal1, eVal2
+    REAL  :: eVA, eVO, CG2, NEWDAC, NEWAC, OLDAC, MAXDAC
     REAL  :: CG1(0:NK+1), WN1(0:NK+1)
     LOGICAL :: LCONVERGED(NSEAL), lexist, LLWS(NSPEC)
 #ifdef WEIGHTS
@@ -6328,7 +6340,7 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL
 #endif
   END SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK
   !/ ------------------------------------------------------------------- /
-  SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC)
+  SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, RD10, RD20, DTG, VGX, VGY, LCALC)
     !/
     !/                  +-----------------------------------+
     !/                  | WAVEWATCH III           NOAA/NCEP |
@@ -6402,7 +6414,7 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC)
 
     INTEGER, INTENT(IN) :: IMOD
 
-    REAL, INTENT(IN)    :: FACX, FACY, DTG, VGX, VGY
+    REAL, INTENT(IN)    :: FACX, FACY, DTG, VGX, VGY, RD10, RD20
 
     REAL              :: KTMP(3), UTILDE(NTH), ST(NTH,NPA)
     REAL              :: FL11(NTH), FL12(NTH), FL21(NTH), FL22(NTH), FL31(NTH), FL32(NTH), KKSUM(NTH,NPA)
@@ -6411,7 +6423,7 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC)
     REAL              :: KSIG(NPA), CGSIG(NPA), CXX(NTH,NPA), CYY(NTH,NPA)
     REAL              :: LAMBDAX(NTH), LAMBDAY(NTH)
     REAL              :: DTMAX(NTH), DTMAXEXP(NTH), DTMAXOUT, DTMAXGL
-    REAL              :: FIN(1), FOUT(1), REST, CFLXY, RD1, RD2, RD10, RD20
+    REAL              :: FIN(1), FOUT(1), REST, CFLXY, RD1, RD2
     REAL              :: UOLD(NTH,NPA), U(NTH,NPA)
 
     REAL, PARAMETER   :: ONESIXTH = 1.0/6.0
@@ -6570,8 +6582,8 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC)
         IF ( FLBPI ) THEN
           DO ITH = 1, NTH
             ISP = ITH + (IK-1) * NTH
-            RD1 = RD10 - DTG * REAL(ITER(IK)-IT)/REAL(ITER(IK))
-            RD2 = RD20
+            RD1=RD10 - DTMAXGL * REAL(ITER(IK)-IT)/REAL(ITER(IK))
+            RD2=RD20
             IF ( RD2 .GT. 0.001 ) THEN
               RD2    = MIN(1.,MAX(0.,RD1/RD2))
               RD1    = 1. - RD2

model/src/w3ref1md.F90

@@ -251,6 +251,7 @@ SUBROUTINE W3SREF(A, CG, WN, EMEAN, FMEAN, DEPTH, CX1, CY1, REFLC, REFLD,     &
     IGBCOVERWRITE =(MOD( NINT(IGPARS(4)),2).EQ.1)
     IGSWELLMAX =( NINT(IGPARS(4)).GE.2)
     ! This following line is a quick fix before the bug is understood ....
+    ! AR: which bug? 
     IF (GTYPE.EQ.UNGTYPE) IGSWELLMAX =.FALSE.
     IGFAC1 = 0.25
     IGFAC2 = 0.25
@@ -321,8 +322,8 @@ SUBROUTINE W3SREF(A, CG, WN, EMEAN, FMEAN, DEPTH, CX1, CY1, REFLC, REFLD,     &
 #endif
     HS=4.*SQRT(EMEANA)
 #ifdef W3_IG1
-    ATMP(:)=A(:)      ! the IG energy will be added to this ATMP
-    ATMP2(:)=A(:)     ! this is really to keep in memory the original spectrum
+    ATMP(:) = A(:)      ! the IG energy will be added to this ATMP
+    ATMP2(:) = A(:)     ! this is really to keep in memory the original spectrum
     IF (IGBCOVERWRITE.AND.REFLC(1).GT.0) THEN
       IGFAC1 = 1.
       ATMP2(1:NSPECIGSTART) = 0.
@@ -349,12 +350,11 @@ SUBROUTINE W3SREF(A, CG, WN, EMEAN, FMEAN, DEPTH, CX1, CY1, REFLC, REFLD,     &
         ENDIF
         FREQIG=SIG(NINT(IGPARS(5)))*TPIINV
         !
-        NSPECIG=NKIG*NTH
-        ATMP(1:NSPECIGSTART)=0.               ! flat bottom approximation (Hasselmann 1962)
-        !                 is not valid for long waves
+        NSPECIG = NKIG * NTH
+        ATMP(1:NSPECIGSTART) = 0.             ! flat bottom approximation (Hasselmann 1962), is not valid for long waves (long?)
         IF (NINT(IGPARS(3)).EQ.1) THEN        ! IGPARS(3) = IGSOURCE
           IF (NINT(IGPARS(8)).EQ.1) THEN      ! in this case, uses depth at break point
-            DEPTHIG=MAX(1.,HS/0.3)            ! to be modified later with a proper gamma
+            DEPTHIG=MAX(1.,HS/0.3)            ! to be modified later with a proper gamma. AR: Strange stuff ...  0.3? 
           ELSE
             DEPTHIG=DEPTH
           END IF

model/src/w3sdb1md.F90

@@ -187,6 +187,7 @@ SUBROUTINE W3SDB1 (IX, A, DEPTH, EMEAN, FMEAN, WNMEAN, CG, LBREAK, S, D )
     USE W3ODATMD, ONLY: NDST
     USE W3GDATMD, ONLY: SIG
     USE W3ODATMD, only : IAPROC
+    USE W3PARALL, only : THR 
 #ifdef W3_S
     USE W3SERVMD, ONLY: STRACE
 #endif
@@ -218,7 +219,7 @@ SUBROUTINE W3SDB1 (IX, A, DEPTH, EMEAN, FMEAN, WNMEAN, CG, LBREAK, S, D )
     INTEGER, SAVE           :: IENT = 0
 #endif
     REAL*8                    :: HM, BB, ARG, Q0, QB, B, CBJ, HRMS, EB(NK)
-    REAL*8                    :: AUX, CBJ2, RATIO, S0, S1, THR, BR1, BR2, FAK
+    REAL*8                    :: AUX, CBJ2, RATIO, S0, S1, BR1, BR2, FAK
     REAL                      :: ETOT, FMEAN2
 #ifdef W3_T0
     REAL                    :: DOUT(NK,NTH)
@@ -235,8 +236,11 @@ SUBROUTINE W3SDB1 (IX, A, DEPTH, EMEAN, FMEAN, WNMEAN, CG, LBREAK, S, D )
     S = 0.
     D = 0.
 
-    THR = DBLE(1.E-15)
-    IF (SUM(A) .LT. THR) RETURN
+    IF (EMEAN .LT. TINY(1.d0)) THEN
+      S = 0 
+      D = 0 
+      RETURN
+    ENDIF
 
     IWB = 1
     !

model/src/w3srcemd.F90

@@ -1244,7 +1244,7 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD,          &
       IF (.NOT. FSSOURCE .or. LSLOC) THEN
 #endif
 #ifdef W3_TR1
-        CALL W3STR1 ( SPEC, SPECOLD, CG1, WN1, DEPTH, IX,        VSTR, VDTR )
+        CALL W3STR1 ( SPEC, CG1, WN1, DEPTH, IX, VSTR, VDTR )
 #endif
 #ifdef W3_PDLIB
       ENDIF
@@ -1534,8 +1534,13 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD,          &
                   DVS  = SIGN(MIN(MAXDAC,ABS(DVS)),DVS)
                 ENDIF
                 PreVS  = DVS / FAKS
-                eVS    = PreVS / CG1(IK) * CLATSL
-                eVD    = MIN(0.,VD(ISP))
+                IF (IOBP_LOC(JSEA) .EQ. 3) THEN
+                  eVS = 0
+                  eVD = 0
+                ELSE
+                  eVS    = PreVS / CG1(IK) * CLATSL
+                  eVD    = MIN(0.,VD(ISP))
+                ENDIF
                 B_JAC(ISP,JSEA)                   = B_JAC(ISP,JSEA) + SIDT * (eVS - eVD*SPEC(ISP)*JAC)
                 ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) = ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) - SIDT * eVD
 #ifdef W3_DB1
@@ -1548,9 +1553,9 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD,          &
                   evS = -evS
                   evD = 2*evD
                 ENDIF
-#endif
                 B_JAC(ISP,JSEA)                   = B_JAC(ISP,JSEA) + SIDT * eVS
                 ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) = ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) - SIDT * eVD
+#endif
 
 #ifdef W3_TR1
                 eVS = VSTR(ISP) * JAC
@@ -1562,9 +1567,9 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD,          &
                   evS = -evS
                   evD = 2*evD
                 ENDIF
-#endif
                 B_JAC(ISP,JSEA)                   = B_JAC(ISP,JSEA) + SIDT * eVS
                 ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) = ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) - SIDT * eVD
+#endif
               END DO
             END DO
 

model/src/w3str1md.F90

@@ -180,7 +180,7 @@ MODULE W3STR1MD
   !>
   !> @author A. J. van der Westhuysen  @date 13-Jan-2013
   !>
-  SUBROUTINE W3STR1 (A, AOLD, CG, WN, DEPTH, IX, S, D)
+  SUBROUTINE W3STR1 (A, CG, WN, DEPTH, IX, S, D)
     !/
     !/                  +-----------------------------------+
     !/                  | WAVEWATCH III           NOAA/NCEP |
@@ -259,7 +259,6 @@ SUBROUTINE W3STR1 (A, AOLD, CG, WN, DEPTH, IX, S, D)
     !       CG      R.A.  I   Group velocities.
     !       WN      R.A.  I   Wavenumbers.
     !       DEPTH   Real  I   Mean water depth.
-    !       EMEAN   Real  I   Mean wave energy.
     !       FMEAN   Real  I   Mean wave frequency.
     !       S       R.A.  O   Source term (1-D version).
     !       D       R.A.  O   Diagonal term of derivative (1-D version).
@@ -320,7 +319,7 @@ SUBROUTINE W3STR1 (A, AOLD, CG, WN, DEPTH, IX, S, D)
     !/ ------------------------------------------------------------------- /
     !/ Parameter list
     !/
-    REAL, INTENT(IN)        :: CG(NK), WN(NK), DEPTH, A(NSPEC), AOLD(NSPEC) 
+    REAL, INTENT(IN)        :: CG(NK), WN(NK), DEPTH, A(NSPEC)
     INTEGER, INTENT(IN)     :: IX
     REAL, INTENT(OUT)       :: S(NSPEC), D(NSPEC)
     !/
@@ -391,11 +390,15 @@ SUBROUTINE W3STR1 (A, AOLD, CG, WN, DEPTH, IX, S, D)
 #ifdef W3_S
     CALL STRACE (IENT, 'W3STR1')
 #endif
-
-!AR: todo: check all PRX routines for differences, check original thesis of elderberky. 
 !
 ! 1.  Integral over directions
 !
+    IF (MAXVAL(A) .LT. TINY(1.)) THEN
+      S = 0 
+      D = 0 
+      RETURN 
+    ENDIF 
+
     SIGM01 = 0.
     EMEAN  = 0.
     JACEPS = 1E-12