diff --git a/.gitignore b/.gitignore
index 2f0e9c116e63c87469bb54c12a0521ffadb8f722..ed0815650fb51e9d4d826d72dfc99122e0bd3c0e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,4 +2,14 @@
 .deps
 .dirstamp
 *.mod
-*~
\ No newline at end of file
+*~
+dmm
+*.nml
+*.m
+*.PRI
+*.DIS
+*.FST
+*.INN
+*.ML
+*.PAR
+*.UNB
diff --git a/Makefile b/Makefile
index 1404d6b24667958275009685dfca620bb725ef57..1863a68a263a0b7bf6292dc5fe290002cc1653f1 100644
--- a/Makefile
+++ b/Makefile
@@ -106,30 +106,43 @@ RANDLIB_OBJS = \
 	sexpo.o \
 	spofa.o
 
-MATLAB_OBJS = ReadMatLabDesign.o
-
 EXEC = dmm
 
 VPATH := $(VPATH) randlib
-LIBS = -llapack -ldl
 
-all: $(MOD_OBJS) $(OBJS) $(RANDLIB_OBJS)
+ifdef DLL
+UNAME_S := $(shell uname -s)
+ifeq ($(UNAME_S), Linux)
+MATLAB_LIBS = -L$(MATLABROOT)/bin/glnxa64 -Wl,-rpath-link,$(MATLABROOT)/bin/glnxa64 -Wl,-rpath,$(MATLABROOT)/bin/glnxa64
+endif
+ifeq ($(UNAME_S), Darwin)
+MATLAB_LIBS = -L$(MATLABROOT)/bin/maci64
+endif
+LIBS = $(MATLAB_LIBS) -llapack -leng -lmx
+DLL_OBJS += design.o setfilem.o geterrstr.o
+INCLUDE = -I$(MATLABROOT)/extern/include
+DEFINE = -DORIGDLL
+else
+LIBS = -llapack
+endif
+
+all: $(MOD_OBJS) $(OBJS) $(RANDLIB_OBJS) $(DLL_OBJS)
 	$(FC) $(FCFLAGS) $^ $(LIBS) -o $(EXEC)
 
 %.o: %.f90 %.mod
-	$(FC) $(FCFLAGS)  -c $<
+	$(FC) $(FCFLAGS) $(INCLUDE) $(DEFINE) -c $<
 
 %.o: %.f90
-	$(FC) $(FCFLAGS) -c $<
+	$(FC) $(FCFLAGS) $(INCLUDE) $(DEFINE) -c $<
 
 %.mod: %.f90 %.o
 	@true
 
 %.o : %.f
-	$(FC) $(FFFLAGS) -c $<
+	$(FC) $(FFFLAGS) $(INCLUDE) $(DEFINE) -c $<
 
 %.o : %.for
-	$(FC) $(FFFLAGS) -c $<
+	$(FC) $(FFFLAGS) $(INCLUDE) $(DEFINE) -c $<
 
 clean:
-	rm -f *.o $(EXEC)
+	rm -f *.o $(EXEC) *.mod *.PRI *.DIS *.FST *.INN *.ML *.PAR *.UNB
diff --git a/amh.for b/amh.for
index 06b63f71ecda4495c8f9160aed8c9b9ca76c94e9..d4bd92bea2c3043edd84170807d1f1bb4e68db3f 100644
--- a/amh.for
+++ b/amh.for
@@ -58,7 +58,7 @@ C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C -------------------------------------------------------------
 	SUBROUTINE AMH(HFIX,nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np,
-	1               yk,IYK,theta,psi,PTR,PM,INFOS,pdll,Z,S,ACCRATE)
+	1               yk,IYK,theta,psi,PTR,PM,INFOS,Z,S,ACCRATE)
 #if defined(__CYGWIN32__) || defined(_WIN32)
 #ifdef __INTEL_COMPILER
       USE dfwin
@@ -68,22 +68,6 @@ C -------------------------------------------------------------
       USE ISO_C_UTILITIES
       USE DLFCN
 #endif
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  POINTER (pdll,fittizia)
-	  POINTER (pdesign,DESIGN)
-#else
-	  TYPE(C_PTR) :: pdll
-	  TYPE(C_FUNPTR) :: pdesign
-	  PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-#endif
 
 C INPUT
 	INTEGER HFIX,nobs,d(2),ny,nz,nx,nu,nv,ns(6),nstot,nt,np,
@@ -123,16 +107,9 @@ C LOCALS
 	id1   = max(1,d(1))
 	Z0    = Z
 	delta = 1.D-3
-#if defined(__CYGWIN32__) || defined(_WIN32)
-      pdesign = getprocaddress(pdll, "design_"C)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 	  CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-	  pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-      IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-         WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-	  END IF
-	  CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-	  CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
 	CALL DESIGNZ(nv,np,psi,INFOS,P1,P2,P3,P4,P5,P6)
 C PALL(i,j) = Pr[Z(t+1)=i|Z(t)=j], Z = S1 x S2 x ... x Snv
diff --git a/amh2.for b/amh2.for
index 76ac605528b38a1790d7410bee7594c1f579b757..83218e3136bf55e6f9fc04782ba8e29d591f6345 100644
--- a/amh2.for
+++ b/amh2.for
@@ -58,7 +58,7 @@ C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C -------------------------------------------------------------
 	SUBROUTINE AMH2(HFIX,nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np,yk,
-	1                theta,psi,PTR,PM,INFOS,pdll,Z,S,ACCRATE)
+	1                theta,psi,PTR,PM,INFOS,Z,S,ACCRATE)
 #if defined(__CYGWIN32__) || defined(_WIN32)
 #ifdef __INTEL_COMPILER
       USE dfwin
@@ -68,22 +68,6 @@ C -------------------------------------------------------------
       USE ISO_C_UTILITIES
       USE DLFCN
 #endif
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  POINTER (pdll,fittizia)	!  ASSOCIATE  pointer P alla DLL ad una varibile fittizia
-	  POINTER (pdesign,DESIGN)	! IMPORTANT associo il puntatore pdesign alla Interface definita
-#else
-	  TYPE(C_PTR) :: pdll
-      TYPE(C_FUNPTR) :: pdesign=C_NULL_FUNPTR
-	  PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-#endif
 
 C INPUT
 	INTEGER HFIX,nobs,d(2),ny,nz,nx,nu,nv,ns(6),nstot,nt,np,
@@ -123,16 +107,9 @@ C LOCALS
 	id1   = max(d(1),1)
 	Z0    = Z
 	delta = 1.D-3
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  pdesign = getprocaddress(pdll, "design_"C)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 	  CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-	  pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-      IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-         WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-	  END IF
-	  CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-	  CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
 	CALL DESIGNZ(nv,np,psi,INFOS,P1,P2,P3,P4,P5,P6)
 C PMAT(i,j) = Pr[Z(t+1)=i|Z(t)=j], Z = S1 x S2 x ... x Snv
diff --git a/checkdesign.for b/checkdesign.for
index 438b042f8fe1daf1d64717ed51c0e75ff4e96374..2638169adfb45696e2a14edb6fc7a76273b792d1 100644
--- a/checkdesign.for
+++ b/checkdesign.for
@@ -34,7 +34,7 @@ C
 C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C -------------------------------------------------------------
-	SUBROUTINE CHECKDESIGN(ny,nz,nx,nu,ns,nt,d,theta,pdll,PATH,NMLNAME)
+	SUBROUTINE CHECKDESIGN(ny,nz,nx,nu,ns,nt,d,theta,PATH,NMLNAME)
 #if defined(__CYGWIN32__) || defined(_WIN32)
 #ifdef __INTEL_COMPILER
       USE dfwin
@@ -48,23 +48,6 @@ C -------------------------------------------------------------
       USE DLFCN
 #endif
 
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  POINTER (pdll,fittizia)	! ASSOCIATE  pointer pdll alla DLL ad una varibile fittizia
-	  POINTER (pdesign,DESIGN)
-#else
-	  TYPE(C_PTR) :: pdll
-      TYPE(C_FUNPTR) :: pdesign=C_NULL_FUNPTR
-	  PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-#endif
-
 C INPUT
 	INTEGER ny,nz,nx,nu,ns(6),nt,d(2)
 	DOUBLE PRECISION theta(nt)
@@ -87,16 +70,9 @@ C EXTERNAL SUBROUTINES
 	ALLOCATE(c(ny,max(nz,1),ns(1)),H(ny,nx,ns(2)),
 	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))) !,HRG(ny,nu),HRGRH(ny,ny))
 
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  pdesign = getprocaddress(pdll, "design_"C)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 	  CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-	  pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-      IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-         WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-	  END IF
-	  CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-	  CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
 
 	maxnz = max(1,nz)
diff --git a/chi2inv.for b/chi2inv.for
index 842c48701a98aa4ad97b6958eb2d5ce9a4675ac3..40c4b7ed162746602c7d843224b97b48c67f4743 100644
--- a/chi2inv.for
+++ b/chi2inv.for
@@ -430,10 +430,19 @@ C LOCALS
       DATA T(198,6:10)/199.8400, 205.0857, 211.0344, 218.6401, 231.8292/
       DATA T(199,6:10)/200.8463, 206.1051, 212.0684, 219.6922, 232.9118/
       DATA T(200,6:10)/201.8526,207.1244, 213.1022, 220.7441, 233.9943/
-
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+      CHARACTER(len=200) :: MEXPRINT
+      INTEGER*4 mexPrintf
+      INTEGER*4 mpfout
+#endif
       IF (V.GT.200) THEN
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+         WRITE(MEXPRINT,*) 'CHI2INV: Too many degrees of freedom'
+         mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
        WRITE(*,*) 'CHI2INV: Too many degrees of freedom'
        PAUSE
+#endif
        RETURN
       ENDIF
 
diff --git a/ReadMatLabDesign.F b/design.for
similarity index 51%
rename from ReadMatLabDesign.F
rename to design.for
index 577270c21bfa5840b91b582c1823d57dfc0af5c4..76ede64648428e64191e580a743711bedbaddd09 100644
--- a/ReadMatLabDesign.F
+++ b/design.for
@@ -50,18 +50,36 @@ C
 	DOUBLE PRECISION nsd(6)
 	INTEGER ny_ptr,nz_ptr,nx_ptr,nu_ptr,ns_ptr,nt_ptr,theta_ptr
 	INTEGER C_ptr, H_ptr, G_ptr, A_ptr, F_ptr, R_ptr
-
+      CHARACTER*1024 matlaberror
 
 C Try to open MatLab (just the first time)
       IF (ep .eq.0 ) THEN
-        ep = engOpen('matlab ')
-        IF (ep .eq. 0) THEN
-          ny = 0 ! ' Can''t start MatLab engine'
-          RETURN
+        ep = engOpen('/Applications/MATLAB_R2014b.app/bin/matlab ')
+        IF (ep .eq. 0) THEN ! Can''t start Matlab engine
+#ifdef __GFORTRAN__
+		   WRITE(*,*) ' '
+		   WRITE(*,*) ' Can''t start MATLAB engine'
+		   WRITE(*,*) ' Program aborting'
+#else
+		   TYPE *, ' '
+		   TYPE *, ' Can''t start MATLAB engine'
+		   TYPE *, ' Program aborting'
+		   PAUSE
+#endif
+		   STOP
         ENDIF
-        IF (engEvalString(ep,'cd ' // pathmfile).ne. 0) then
-          ny = -7 ! ' Can''t find or open the MatLab funtion'
-          RETURN
+        IF (engEvalString(ep,'cd ' // pathmfile).ne. 0) then ! Can''t find or open the MatLab funtion
+#ifdef __GFORTRAN__
+		   WRITE(*,*) ' '
+		   WRITE(*,*) ' Can''t find or open the MatLab function'
+		   WRITE(*,*) ' Program aborting'
+#else
+		   TYPE *, ' '
+		   TYPE *, ' Can''t find or open the MatLab function'
+		   TYPE *, ' Program aborting'
+		   PAUSE
+#endif
+		   STOP
         ENDIF
       ENDIF
 
@@ -74,20 +92,37 @@ C
 #else
       status = engPutVariable(ep, 'ny'C, ny_ptr)
 #endif
-      IF (status .ne. 0) THEN
-         ny = -1 ! ' Can''t read ny in the MatLab file'
-         RETURN
+      IF (status .ne. 0) THEN ! Can''t read ny in the Matlab file
+#ifdef __GFORTRAN__
+		 WRITE(*,*) ' '
+		 WRITE(*,*) ' Can''t read ny in the MATLAB file'
+		 WRITE(*,*) ' Program aborting'
+#else
+		 TYPE *, ' '
+		 TYPE *, ' Can''t read ny in the MATLAB file'
+		 TYPE *, ' Program aborting'
+		 PAUSE
+#endif
+		 STOP
       ENDIF
-
 	nz_ptr = mxCreateDoubleScalar(nz*1.0d0)
 #ifdef __GFORTRAN__
       status = engPutVariable(ep, 'nz', nz_ptr)
 #else
       status = engPutVariable(ep, 'nz'C, nz_ptr)
 #endif
-      IF (status .ne. 0) THEN
-         ny = -2 ! ' Can''t read nz in the MatLab file'
-         RETURN
+      IF (status .ne. 0) THEN ! Can''t read nz in the Matlab file
+#ifdef __GFORTRAN__
+		 WRITE(*,*) ' '
+		 WRITE(*,*) ' Can''t read nz in the MATLAB file'
+		 WRITE(*,*) ' Program aborting'
+#else
+		 TYPE *, ' '
+		 TYPE *, ' Can''t read nz in the MATLAB file'
+		 TYPE *, ' Program aborting'
+		 PAUSE
+#endif
+		 STOP
       ENDIF
 
 	nx_ptr = mxCreateDoubleScalar(nx*1.0d0)
@@ -96,9 +131,18 @@ C
 #else
       status = engPutVariable(ep, 'nx'C, nx_ptr)
 #endif
-      IF (status .ne. 0) THEN
-         ny = -3 ! ' Can''t read nx in the MatLab file'
-         RETURN
+      IF (status .ne. 0) THEN ! Can''t read nx in the Matlab file
+#ifdef __GFORTRAN__
+		 WRITE(*,*) ' '
+		 WRITE(*,*) ' Can''t read nx in the MATLAB file'
+		 WRITE(*,*) ' Program aborting'
+#else
+		 TYPE *, ' '
+		 TYPE *, ' Can''t read nx in the MATLAB file'
+		 TYPE *, ' Program aborting'
+		 PAUSE
+#endif
+		 STOP
       ENDIF
 
 	nu_ptr = mxCreateDoubleScalar(nu*1.0d0)
@@ -107,9 +151,18 @@ C
 #else
       status = engPutVariable(ep, 'nu'C, nu_ptr)
 #endif
-      IF (status .ne. 0) THEN
-         ny = -4 ! ' Can''t read nu in the MatLab file'
-         RETURN
+      IF (status .ne. 0) THEN ! Can''t read nu in the MatLab file
+#ifdef __GFORTRAN__
+		 WRITE(*,*) ' '
+		 WRITE(*,*) ' Can''t read nu in the MATLAB file'
+		 WRITE(*,*) ' Program aborting'
+#else
+		 TYPE *, ' '
+		 TYPE *, ' Can''t read nu in the MATLAB file'
+		 TYPE *, ' Program aborting'
+		 PAUSE
+#endif
+		 STOP
       ENDIF
 
       ns_ptr = mxCreateDoubleMatrix(1, 6, 0)
@@ -122,9 +175,18 @@ C
 #else
       status = engPutVariable(ep, 'ns'C, ns_ptr)
 #endif
-      IF (status .ne. 0) THEN
-         ny = -5 ! ' Can''t read ns in the MatLab file'
-         RETURN
+      IF (status .ne. 0) THEN ! Can''t read ns in the Matlab file
+#ifdef __GFORTRAN__
+		 WRITE(*,*) ' '
+		 WRITE(*,*) ' Can''t read ns in the MATLAB file'
+		 WRITE(*,*) ' Program aborting'
+#else
+		 TYPE *, ' '
+		 TYPE *, ' Can''t read ns in the MATLAB file'
+		 TYPE *, ' Program aborting'
+		 PAUSE
+#endif
+		 STOP
       ENDIF
 
       theta_ptr = mxCreateDoubleMatrix(1, nt, 0)
@@ -134,9 +196,18 @@ C
 #else
       status = engPutVariable(ep, 'theta'C, theta_ptr)
 #endif
-      IF (status .ne. 0) THEN
-         ny = -6 ! ' Can''t read theta in the MatLab file'
-         RETURN
+      IF (status .ne. 0) THEN ! Can''t read theta in the Matlab file
+#ifdef __GFORTRAN__
+		 WRITE(*,*) ' '
+		 WRITE(*,*) ' Can''t read nt in the MATLAB file'
+		 WRITE(*,*) ' Program aborting'
+#else
+		 TYPE *, ' '
+		 TYPE *, ' Can''t read nt in the MATLAB file'
+		 TYPE *, ' Program aborting'
+		 PAUSE
+#endif
+		 STOP
       ENDIF
 
 C
@@ -146,19 +217,51 @@ C
       status = engOutputBuffer(ep, buffer1)
       IF (engEvalString(ep, 'clear success;'//
      &   '[C,H,G,A,F,R]='//TRIM(mfile)//'( ny,nz,nx,'//
-     &   'nu,ns,theta);'//'success=1;') .ne. 0) then
-         ny = -8   ! engEvalString failed
-         RETURN
+     &   'nu,ns,theta);'//'success=1;') .ne. 0) then ! engEvalString failed
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+		 CALL GETERRSTR(matlaberror)
+#else
+#endif
+
+#ifdef __GFORTRAN__
+		 WRITE(*,*) ' '
+		 WRITE(*,*) ' the MATLAB funtion can not be executed:'
+		 WRITE(*,*) trim(matlaberror)
+		 WRITE(*,*) ' Program aborting'
+#else
+		 TYPE *, ' '
+		 TYPE *, ' the MATLAB funtion can not be executed:'
+		 TYPE *, trim(matlaberror)
+		 TYPE *, ' Program aborting'
+		 PAUSE
+#endif
+		 STOP
       ENDIF
 #ifdef __GFORTRAN__
       C_ptr = engGetVariable(ep, 'success')
 #else
       C_ptr = engGetVariable(ep, 'success'C)
 #endif
-      IF (C_ptr .eq. 0) then
+      IF (C_ptr .eq. 0) then ! engEvalString failed
           buffer=buffer1
-          ny=-8   ! engEvalString failed
-          RETURN
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+		 CALL GETERRSTR(matlaberror)
+#else
+#endif
+
+#ifdef __GFORTRAN__
+		 WRITE(*,*) ' '
+		 WRITE(*,*) ' the MATLAB funtion can not be executed:'
+		 WRITE(*,*) trim(matlaberror)
+		 WRITE(*,*) ' Program aborting'
+#else
+		 TYPE *, ' '
+		 TYPE *, ' the MATLAB funtion can not be executed:'
+		 TYPE *, trim(matlaberror)
+		 TYPE *, ' Program aborting'
+		 PAUSE
+#endif
+		 STOP
       ENDIF
 
 C
@@ -172,7 +275,17 @@ C
       IF(C_ptr.NE.0) THEN
        CALL mxCopyPtrToReal8(mxGetPr(C_ptr), c, ny*max(1,nz)*ns(1))
       ELSE
-       ny = -101
+#ifdef __GFORTRAN__
+		 WRITE(*,*) ' '
+		 WRITE(*,*) ' C could not be assigned during the call'
+		 WRITE(*,*) ' Program aborting'
+#else
+		 TYPE *, ' '
+		 TYPE *, ' C could not be assigned during the call '
+		 TYPE *, ' Program aborting'
+		 PAUSE
+#endif
+		 STOP
       ENDIF
 
 #ifdef __GFORTRAN__
@@ -183,7 +296,17 @@ C
       IF(H_ptr.NE.0) THEN
        CALL mxCopyPtrToReal8(mxGetPr(H_ptr), H, ny*nx*ns(2))
       ELSE
-       ny = -102
+#ifdef __GFORTRAN__
+		 WRITE(*,*) ' '
+		 WRITE(*,*) ' H could not be assigned during the call'
+		 WRITE(*,*) ' Program aborting'
+#else
+		 TYPE *, ' '
+		 TYPE *, ' H could not be assigned during the call '
+		 TYPE *, ' Program aborting'
+		 PAUSE
+#endif
+		 STOP
       ENDIF
 
 #ifdef __GFORTRAN__
@@ -194,7 +317,17 @@ C
       IF(G_ptr.NE.0) THEN
        CALL mxCopyPtrToReal8(mxGetPr(G_ptr), G, ny*nu*ns(3))
       ELSE
-       ny = -103
+#ifdef __GFORTRAN__
+		 WRITE(*,*) ' '
+		 WRITE(*,*) ' G could not be assigned during the call'
+		 WRITE(*,*) ' Program aborting'
+#else
+		 TYPE *, ' '
+		 TYPE *, ' G could not be assigned during the call '
+		 TYPE *, ' Program aborting'
+		 PAUSE
+#endif
+		 STOP
       ENDIF
 
 #ifdef __GFORTRAN__
@@ -205,7 +338,17 @@ C
       IF(A_ptr.NE.0) THEN
        CALL mxCopyPtrToReal8(mxGetPr(A_ptr), a, nx*ns(4))
       ELSE
-       ny = -104
+#ifdef __GFORTRAN__
+		 WRITE(*,*) ' '
+		 WRITE(*,*) ' A could not be assigned during the call'
+		 WRITE(*,*) ' Program aborting'
+#else
+		 TYPE *, ' '
+		 TYPE *, ' A could not be assigned during the call '
+		 TYPE *, ' Program aborting'
+		 PAUSE
+#endif
+		 STOP
       ENDIF
 
 #ifdef __GFORTRAN__
@@ -216,20 +359,39 @@ C
       IF(F_ptr.NE.0) THEN
        CALL mxCopyPtrToReal8(mxGetPr(F_ptr), F, nx*nx*ns(5))
       ELSE
-       ny = -105
+#ifdef __GFORTRAN__
+		 WRITE(*,*) ' '
+		 WRITE(*,*) ' F could not be assigned during the call'
+		 WRITE(*,*) ' Program aborting'
+#else
+		 TYPE *, ' '
+		 TYPE *, ' F could not be assigned during the call '
+		 TYPE *, ' Program aborting'
+		 PAUSE
+#endif
+		 STOP
       ENDIF
 
 #ifdef __GFORTRAN__
-      r_ptr = enggetvariable(ep, 'r')
+      r_ptr = enggetvariable(ep, 'R')
 #else
-      r_ptr = enggetvariable(ep, 'r'c)
+      r_ptr = enggetvariable(ep, 'R'c)
 #endif
       IF(R_ptr.NE.0) THEN
        CALL mxCopyPtrToReal8(mxGetPr(R_ptr), R, nx*nu*ns(6))
       ELSE
-       ny = -106
+#ifdef __GFORTRAN__
+		 WRITE(*,*) ' '
+		 WRITE(*,*) ' R could not be assigned during the call'
+		 WRITE(*,*) ' Program aborting'
+#else
+		 TYPE *, ' '
+		 TYPE *, ' R could not be assigned during the call '
+		 TYPE *, ' Program aborting'
+		 PAUSE
+#endif
+		 STOP
       ENDIF
-
 C
 C Free dynamic memory allocated by MXCREATE function
 C
@@ -248,27 +410,3 @@ C
 
       RETURN
       END
-
-C -----------------------------------------
-C To make dynamic the name of the .m file
-C -----------------------------------------
-	SUBROUTINE SETFILEM(string1,string2)
-!DEC$ ATTRIBUTES DLLEXPORT, ALIAS:'setfilem_' :: SETFILEM
-	CHARACTER*200 string1,string2,mfile,pathmfile
-      COMMON /M/ mfile, pathmfile
-      mfile     = string1
-      pathmfile = string2
-	RETURN
-      END
-
-C --------------------
-C To get MatLab errors
-C --------------------
-      SUBROUTINE GETERRSTR(matlaberror)
-!DEC$ ATTRIBUTES DLLEXPORT, ALIAS:'geterrstr_' :: GETERRSTR
-      CHARACTER*1024 matlaberror
-      CHARACTER*1024 buffer
-      COMMON /ERRBUFFER/ buffer
-	 matlaberror = buffer
-      RETURN
-      END
diff --git a/drawtheta.for b/drawtheta.for
index 1e17c32190560d5068b37005be62c51f71c17476..475fb179d648b6aaadb08165dff164ccb5eaef19 100644
--- a/drawtheta.for
+++ b/drawtheta.for
@@ -105,7 +105,9 @@ C CHEK theta
 	   IF (NN.LE.1000) THEN
 	     GOTO 7777
 	   ELSE
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nReduce skcriterium or use Slice sampling\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
           WRITE(*,*) ' '
           WRITE(*,*) 'Reduce skcriterium or use Slice sampling'
           WRITE(*,*) 'Program aborting'
diff --git a/drawtheta2.for b/drawtheta2.for
index c2813d4a175d7e260d4e5d4451a6bf98f4c7a858..f8c1a30bf2a2dd67058e182a92a918afe1c91785 100644
--- a/drawtheta2.for
+++ b/drawtheta2.for
@@ -105,7 +105,9 @@ C CHEK theta
 	   IF (NN.LE.1000) THEN
 	     GOTO 7777
 	   ELSE
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nReduce skcriterium or use Slice sampling\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
           WRITE(*,*) ' '
           WRITE(*,*) 'Reduce skcriterium or use Slice sampling'
           WRITE(*,*) 'Program aborting'
diff --git a/examples/NILE.NML b/examples/NILE.NML
index c79eb74a6cb652b9594a58cc97499e30d2c212d8..d0cf2cfdaf3aaeeb44e237792fa77e97d8e03fd1 100644
--- a/examples/NILE.NML
+++ b/examples/NILE.NML
@@ -36,7 +36,7 @@
  check   = to debug the system matrices c, H, G, a, F, and R
 
  &ssm
-  nx=1 nu=2 d=1 1 nv=2 check=n estimation=BA dllname = path\NILE1.m
+  nx=1 nu=2 d=1 1 nv=2 check='n' estimation='BA' dllname = 'path\NILE1.m'
  &end
 
  Namelist Sj describes discrete latent variables Sj, j = 1,...,nv:
@@ -47,9 +47,9 @@
             of dimension k ordered by columns
  matSj    = matrices impacted by Sj (one or more of c H G a F R)
 
- &S1 dynS1=I nS1=2 hypS1(1,1)=16 2 matS1=G &end
+ &S1 dynS1='I' nS1=2 hypS1(1,1)=16 2 matS1='G' &end
 
- &S2 dynS2=I nS2=2 hypS2(1,1)=16 2 matS2=R &end
+ &S2 dynS2='I' nS2=2 hypS2(1,1)=16 2 matS2='R' &end
 
  Namelist prior describes priors pdf of model parameters (except transition probs):
  nt       = number of theta parameters
@@ -59,25 +59,25 @@
             value is fixed at lb=ub
 
  &prior nt = 3
-  pdftheta(1) = IG hyptheta(1,1) = 6000 6  0 500000
-  pdftheta(2) = IG hyptheta(1,2) = 6000 6  0 500000
-  pdftheta(3) = BE hyptheta(1,3) = 2 4  1 20
+  pdftheta(1) = 'IG' hyptheta(1,1) = 6000 6  0 500000
+  pdftheta(2) = 'IG' hyptheta(1,2) = 6000 6  0 500000
+  pdftheta(3) = 'BE' hyptheta(1,3) = 2 4  1 20
  &end
 
   ML:
-  pdftheta(1) = IG hyptheta(1,1) = 10000 6  0 50000
-  pdftheta(2) = IG hyptheta(1,2) = 10000 6  0 50000
-  pdftheta(3) = BE hyptheta(1,3) = 6     4  1 20
+  pdftheta(1) = 'IG' hyptheta(1,1) = 10000 6  0 50000
+  pdftheta(2) = 'IG' hyptheta(1,2) = 10000 6  0 50000
+  pdftheta(3) = 'BE' hyptheta(1,3) = 6     4  1 20
 
   Fixed:
-  pdftheta(1) = IG hyptheta(1,1) = 60000 6 12000 12000
-  pdftheta(2) = IG hyptheta(1,2) = 60000 6 9000 9000
-  pdftheta(3) = BE hyptheta(1,3) = 2 4 4 4
+  pdftheta(1) = 'IG' hyptheta(1,1) = 60000 6 12000 12000
+  pdftheta(2) = 'IG' hyptheta(1,2) = 60000 6 9000 9000
+  pdftheta(3) = 'BE' hyptheta(1,3) = 2 4 4 4
 
   Bayes:
-  pdftheta(1) = IG hyptheta(1,1) = 60000 6 0 50000
-  pdftheta(2) = IG hyptheta(1,2) = 60000 6 0 50000
-  pdftheta(3) = BE hyptheta(1,3) = 2     4 1 20
+  pdftheta(1) = 'IG' hyptheta(1,1) = 60000 6 0 50000
+  pdftheta(2) = 'IG' hyptheta(1,2) = 60000 6 0 50000
+  pdftheta(3) = 'BE' hyptheta(1,3) = 2     4 1 20
 
  Namelist mcmc contains the Markov Chain Monte Carlo options:
  seed        = seed of random number generator (0-999; 0 default)
@@ -87,7 +87,7 @@
  hbl         = block length discrete latent variable (1:GCK, >1: AMH)
 
  &mcmc
-  seed=0 thin=1 burnin=100 simulrec=500 hbl=1 MargLik=y
+  seed=0 thin=1 burnin=100 simulrec=500 hbl=1 MargLik='y'
  &end
 
  Namelist dataset provides data:
@@ -102,7 +102,7 @@
 
  Note: use can be made of -99999 to assign missing values to the endogenous variables
 
- &dataset T=100 ny=1 nz=0 nf=10 datasim=n obs=
+ &dataset T=100 ny=1 nz=0 nf=10 datasim='n' obs=
   1120
   1160
   963
diff --git a/findinput.for b/findinput.for
index c3441d1fe63fff89be8a4d472c82c8fc101966f6..adfee8dfb973cdabe727bb6429c25056ead50aba 100644
--- a/findinput.for
+++ b/findinput.for
@@ -48,9 +48,13 @@ C LOCALS
 	ERRMSG = 'INPUT ERROR: '// STR2 //' is not set'
 	IPOS = INDEX(TRIM(STR1),STR2)
 	IF (IPOS.EQ.0) THEN
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt(ERRMSG)
+#else
 	 WRITE(*,*) ERRMSG
 	 PAUSE
 	 STOP
+#endif
       ELSE
 	 IPOS = IPOS+2
 	 DO WHILE ((OSTR.NE.'=').AND.(OSTR.NE.'$'))
@@ -66,9 +70,13 @@ C LOCALS
 	  READ(STR1(IPOS:IPOS2-2),'(F20.10)') NUM
 	  READ(STR1(IPOS:IPOS2-2),*) NUM
 	 ELSE
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt(ERRMSG)
+#else
 	  WRITE(*,*) ERRMSG
 	  PAUSE
 	  STOP
+#endif
 	 ENDIF
 	ENDIF
 
diff --git a/forecast.for b/forecast.for
index 0309e677eddf950fd52e4937a7493b348836ef48..5a34b5dc798bc1914f80e8134821b842358b41bc 100644
--- a/forecast.for
+++ b/forecast.for
@@ -36,7 +36,7 @@ C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 C GNU General Public License for more details.
 C ----------------------------------------------------------------------------
 	SUBROUTINE FORECAST(zk,nf,ny,nz,nx,nu,nv,ns,nstot,nt,np,
-	1                    theta,psi,INFOS,Z,STATE,pdll,FORE)
+	1                    theta,psi,INFOS,Z,STATE,FORE)
 #if defined(__CYGWIN32__) || defined(_WIN32)
 #ifdef __INTEL_COMPILER
       USE dfwin
@@ -46,23 +46,6 @@ C ----------------------------------------------------------------------------
       USE ISO_C_UTILITIES
       USE DLFCN
 #endif
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-	CHARACTER*1 fittizia
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  POINTER (pdll,fittizia)	! ASSOCIATE  pointer pdll alla DLL ad una varibile fittizia
-	  POINTER (pdesign,DESIGN)
-#else
-	  TYPE(C_PTR) :: pdll
-      TYPE(C_FUNPTR) :: pdesign=C_NULL_FUNPTR
-	  PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-#endif
 
 C INPUT
 	INTEGER nf,ny,nz,nx,nu,nv,nt,np(3),ns(6),nstot,Z,INFOS(9,6)
@@ -90,17 +73,9 @@ C EXTERNAL SUBROUTINES
 	ALLOCATE(R(nx,nu,ns(6)),c(ny,max(nz,1),ns(1)),H(ny,nx,ns(2)),
 	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)))
 
-C Call DESIGN
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  pdesign = getprocaddress(pdll, "design_"C)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 	  CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-	  pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-      IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-         WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-	  END IF
-	  CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-	  CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
 
 	IF (nv.GT.0) THEN
diff --git a/gck.for b/gck.for
index ac564e68331856759cd7f1dfae42f485c6744671..a07ef50157c5bc520168af8ea010de3ee59f0226 100644
--- a/gck.for
+++ b/gck.for
@@ -59,7 +59,7 @@ C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C ----------------------------------------------------------------------
 	SUBROUTINE GCK(nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np,yk,IYK,
-	1               theta,psi,INFOS,pdll,Z,S)
+	1               theta,psi,INFOS,Z,S)
 #if defined(__CYGWIN32__) || defined(_WIN32)
 #ifdef __INTEL_COMPILER
       USE dfwin
@@ -69,23 +69,6 @@ C ----------------------------------------------------------------------
       USE ISO_C_UTILITIES
       USE DLFCN
 #endif
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  POINTER (pdll,fittizia)	!  ASSOCIATE  pointer P alla DLL ad una varibile fittizia
-	  POINTER (pdesign,DESIGN)	! IMPORTANT associo il puntatore pdesign alla Interface definita
-#else
-	  TYPE(C_PTR) :: pdll
-      TYPE(C_FUNPTR) :: pdesign=C_NULL_FUNPTR
-	  PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-#endif
-
 
 C INPUT
 	INTEGER nobs,d(2),ny,nz,nx,nu,nv,ns(6),nstot,nt,np,IYK(nobs,ny+1),
@@ -127,16 +110,9 @@ C LOCALS
 	DATA EPS/1.D-14/,ONE/1.0D0/,ZERO/0.0D0/
 	DOUBLE PRECISION genunf,LEMMA4,MARKOVP
 
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  pdesign = getprocaddress(pdll, "design_"C)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 	  CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-	  pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-      IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-         WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-	  END IF
-	  CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-	  CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
 	CALL DESIGNZ(nv,np,psi,INFOS,P1,P2,P3,P4,P5,P6)
 C PALL(i,j) = Pr[Z(t+1)=i|Z(t)=j], Z = S1 x S2 x ... x Snv
diff --git a/gck2.for b/gck2.for
index bfe3c5948169553108321553a52e332fc2075eeb..cddb24b31839e00d1e861cd1ae6cd02b89d20034 100644
--- a/gck2.for
+++ b/gck2.for
@@ -59,7 +59,7 @@ C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C ----------------------------------------------------------------------
 	SUBROUTINE GCK2(nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np,yk,
-	1                theta,psi,INFOS,pdll,Z,S)
+	1                theta,psi,INFOS,Z,S)
 #if defined(__CYGWIN32__) || defined(_WIN32)
 #ifdef __INTEL_COMPILER
       USE dfwin
@@ -69,22 +69,6 @@ C ----------------------------------------------------------------------
       USE ISO_C_UTILITIES
       USE DLFCN
 #endif
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  POINTER (pdll,fittizia)	!  ASSOCIATE  pointer P alla DLL ad una varibile fittizia
-	  POINTER (pdesign,DESIGN)	! IMPORTANT associo il puntatore pdesign alla Interface definita
-#else
-	  TYPE(C_PTR) :: pdll
-      TYPE(C_FUNPTR) :: pdesign=C_NULL_FUNPTR
-	  PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-#endif
 
 C INPUT
 	INTEGER nobs,d(2),ny,nz,nx,nu,nv,ns(6),nstot,nt,np,
@@ -117,16 +101,9 @@ C LOCALS
 	DATA EPS/1.D-14/,ONE/1.0D0/,ZERO/0.0D0/
 	DOUBLE PRECISION genunf,LEMMA4,MARKOVP
 
-#if defined(__CYGWIN32__) || defined(_WIN32)
-      pdesign = getprocaddress(pdll, "design_"C)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 	  CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-	  pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-      IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-         WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-	  END IF
-	  CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-	  CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
 	CALL DESIGNZ(nv,np,psi,INFOS,P1,P2,P3,P4,P5,P6)
 C PALL(i,j) = Pr[Z(t+1)=i|Z(t)=j], Z = S1 x S2 x ... x Snv
diff --git a/geterrstr.for b/geterrstr.for
new file mode 100644
index 0000000000000000000000000000000000000000..5770e0d9434f32d0c340802c38b6485dd65a5fee
--- /dev/null
+++ b/geterrstr.for
@@ -0,0 +1,32 @@
+C --------------------------------------------------------------------
+C This file crates matlabdll.dll whih is ment to read the .m file
+C
+C Copyright (C) 2010-2014 European Commission
+C
+C This file is part of Program DMM
+C
+C DMM is free software developed at the Joint Research Centre of the
+C European Commission: you can redistribute it and/or modify it under
+C the terms of the GNU General Public License as published by
+C the Free Software Foundation, either version 3 of the License, or
+C (at your option) any later version.
+C
+C DMM is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
+C ---------------------------------------------------------------------
+C --------------------
+C To get MatLab errors
+C --------------------
+      SUBROUTINE GETERRSTR(matlaberror)
+!DEC$ ATTRIBUTES DLLEXPORT, ALIAS:'geterrstr_' :: GETERRSTR
+      CHARACTER*1024 matlaberror
+      CHARACTER*1024 buffer
+      COMMON /ERRBUFFER/ buffer
+	 matlaberror = buffer
+      RETURN
+      END
diff --git a/harmonic.for b/harmonic.for
index 2167f6d4c1fdf4f514c5d5bb52d2114e8d30145c..3ae3e252dadf0ee5dea28e6f4f70edbe3cc106ae 100644
--- a/harmonic.for
+++ b/harmonic.for
@@ -27,7 +27,7 @@ C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C --------------------------------------------------------------------------
 	SUBROUTINE HARMONIC(G,nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np,
 	1                    INFOS,yk,IYK,gibpar,gibZ,thetaprior,
-     2                    psiprior,tipo,pdll,MLH)
+     2                    psiprior,tipo,MLH)
 
 #ifdef __GFORTRAN__
       USE gfortran
@@ -39,7 +39,6 @@ C INPUT
 	DOUBLE PRECISION yk(nobs,ny+nz),gibpar(G,nt+np(1)),
 	1 thetaprior(nt,4),psiprior(np(2),np(3))
 	CHARACTER*2 tipo(nt)
-	POINTER (pdll,fittizia)  !  ASSOCIATE  pointer P alla DLL ad una varibile fittizia
 
 C OUTPUT
 	DOUBLE PRECISION MLH(11,2)
@@ -245,7 +244,7 @@ C Mothod of Moments: a0 = m1(1-m1)/V1+1, ai = mi*a0, i=1,2,..,N
 
 	 Fpar = PTHETA(NPOS(1),nobs,d,ny,nz,nx,nu,ns,nt,IS,yk,IYK,
 	1               gibpar(IG,1:nt),thetaprior(NPOS(1),:),
-     2               tipo(NPOS(1)),pdll)
+     2               tipo(NPOS(1)))
 	 Fpar = Fpar + SUM(Ppar(2:NPARTH+K))  ! log[f(y|par,S)f(par)]
 
 	 DEN2(IG) = QTHETA+QPSI+QS-Fpar-PS
diff --git a/harmonic2.for b/harmonic2.for
index be9b3a8b6225398e4ac1b871cfc180fe17240439..72f17f9ce42a4f7dc06d1d60aa2b56e5e394a2b0 100644
--- a/harmonic2.for
+++ b/harmonic2.for
@@ -28,7 +28,7 @@ C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C --------------------------------------------------------------------------
 	SUBROUTINE HARMONIC2(G,nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np,
 	1                     INFOS,yk,gibpar,gibZ,thetaprior,psiprior,
-     2                     tipo,pdll,MLH)
+     2                     tipo,MLH)
 #ifdef __GFORTRAN__
       USE gfortran
 #endif
@@ -38,7 +38,6 @@ C INPUT
 	DOUBLE PRECISION yk(nobs,ny+nz),gibpar(G,nt+np(1)),
 	1 thetaprior(nt,4),psiprior(np(2),np(3))
 	CHARACTER*2 tipo(nt)
-	POINTER (pdll,fittizia)  !  ASSOCIATE  pointer P alla DLL ad una varibile fittizia
 
 C OUTPUT
 	DOUBLE PRECISION MLH(11,2)
@@ -244,7 +243,7 @@ C Mothod of Moments: a0 = m1(1-m1)/V1+1, ai = mi*a0, i=1,2,..,N
 
 	 Fpar = PTHETA2(NPOS(1),nobs,d,ny,nz,nx,nu,ns,nt,IS,yk,
 	1                gibpar(IG,1:nt),thetaprior(NPOS(1),:),
-     2                tipo(NPOS(1)),pdll)
+     2                tipo(NPOS(1)))
 	 Fpar = Fpar + SUM(Ppar(2:NPARTH+K))  ! log[f(y|par,S)f(par)]
 
 	 DEN2(IG) = QTHETA+QPSI+QS-Fpar-PS
diff --git a/innov.for b/innov.for
index d420494ee5a96552e2c29a25913a3693592f6fcb..f0196f6e7b9a0cdfb7948d50f2f904e70c6ba1fb 100644
--- a/innov.for
+++ b/innov.for
@@ -38,7 +38,7 @@ C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C --------------------------------------------------------------------
       SUBROUTINE INNOV(nobs,d,ny,nz,nx,nu,ns,nt,S,yk,IYK,
-	1                 theta,pdll,INN)
+	1                 theta,INN)
 #if defined(__CYGWIN32__) || defined(_WIN32)
 #ifdef __INTEL_COMPILER
       USE dfwin
@@ -48,22 +48,6 @@ C --------------------------------------------------------------------
       USE ISO_C_UTILITIES
       USE DLFCN
 #endif
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-#if defined(__CYGWIN32__) || defined(_WIN32)
-      POINTER (pdll,fittizia)   !  ASSOCIATE  pointer P alla DLL ad una varibile fittizia
-      POINTER (pdesign,DESIGN)
-#else
-	  TYPE(C_PTR) :: pdll
-      TYPE(C_FUNPTR) :: pdesign=C_NULL_FUNPTR
-      PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-#endif
 
 C INPUT
 	INTEGER nobs,d(2),ny,nz,nx,nu,nt
@@ -88,16 +72,9 @@ C LOCALS
      1 Xdd(max(d(1),1),nx),Pdd(max(d(1),1),nx,nx),LIKE(max(d(1),1)),
      1 XT(nx),PT(nx,nx),INN0(ny))
 
-#if defined(__CYGWIN32__) || defined(_WIN32)
-      pdesign = getprocaddress(pdll, "design_"C)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
       CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-      pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-      IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-         WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-	  END IF
-      CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-	  CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
 	INN(:,:) = 0.D0
 	CALL IKF(d,ny,nz,nx,nu,ns,S(1:max(d(1),1),1:6),
diff --git a/innov2.for b/innov2.for
index ea62123065be6c4fba66e2beebb09d2786ce274e..b55116e3b4d1b252f38ee24a1a6c24eebf51c8de 100644
--- a/innov2.for
+++ b/innov2.for
@@ -38,7 +38,7 @@ C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C --------------------------------------------------------------------
       SUBROUTINE INNOV2(nobs,d,ny,nz,nx,nu,ns,nt,S,yk,
-	1                  theta,pdll,INN)
+	1                  theta,INN)
 #if defined(__CYGWIN32__) || defined(_WIN32)
 #ifdef __INTEL_COMPILER
       USE dfwin
@@ -48,22 +48,6 @@ C --------------------------------------------------------------------
       USE ISO_C_UTILITIES
       USE DLFCN
 #endif
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-#if defined(__CYGWIN32__) || defined(_WIN32)
-      POINTER (pdll,fittizia)   !  ASSOCIATE  pointer P alla DLL ad una varibile fittizia
-      POINTER (pdesign,DESIGN)  ! IMPORTANT associo il puntatore pdesign alla Interface definita
-#else
-	  TYPE(C_PTR) :: pdll
-      TYPE(C_FUNPTR) :: pdesign=C_NULL_FUNPTR
-      PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-#endif
 
 C INPUT
 	INTEGER nobs,d(2),ny,nz,nx,nu,nt,ns(6),S(nobs,6)
@@ -87,16 +71,9 @@ C LOCALS
      1 Xdd(max(d(1),1),nx),Pdd(max(d(1),1),nx,nx),LIKE(max(d(1),1)),
      1 XT(nx),PT(nx,nx))
 
-#if defined(__CYGWIN32__) || defined(_WIN32)
-      pdesign = getprocaddress(pdll, "design_"C)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
       CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-      pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-      IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-         WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-	  END IF
-      CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-	  CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
 	INN(:,:) = 0.D0
 	CALL IKF2(d,ny,nz,nx,nu,ns,S(1:max(d(1),1),1:6),
diff --git a/input.for b/input.for
index 2c645bded31e34df1de29124504b5ad290cd77d3..7918dd421edfa9e41dfc481aa496c4f3e174179c 100644
--- a/input.for
+++ b/input.for
@@ -106,7 +106,9 @@ C FIND namelist ssm
 #endif
 	CLOSE(1)
 	IF (IFAIL.EQ.-1) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nNamelist ssm not found\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Namelist ssm not found'
        WRITE(*,*) ' Program aborting'
 #else
@@ -134,7 +136,9 @@ C READ namelist ssm
 	estimation = 'BA'
 	READ(1,NML=ssm,END=5001,ERR=5001)
 	IF (nx.LE.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck nx in namelist ssm\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
 	 WRITE(*,*) ' Check nx in namelist ssm'
 	 WRITE(*,*) ' Program aborting'
 #else
@@ -145,7 +149,9 @@ C READ namelist ssm
 	 STOP
 	ENDIF
 	IF(nu.LE.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck nu in namelist ssm\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Check nu in namelist ssm'
        WRITE(*,*) ' Program aborting'
 #else
@@ -156,7 +162,9 @@ C READ namelist ssm
 	 STOP
 	ENDIF
 	IF((d(1).LT.0).OR.(d(2).LT.0).OR.(d(2).GT.nx)) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck d in namelist ssm\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Check d in namelist ssm'
        WRITE(*,*) ' Program aborting'
 #else
@@ -167,7 +175,9 @@ C READ namelist ssm
 	 STOP
 	ENDIF
 	IF((nv.LT.0).OR.(nv.GT.6)) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck nv in namelist ssm\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Check nv in namelist ssm'
        WRITE(*,*) ' Program aborting'
 #else
@@ -178,7 +188,9 @@ C READ namelist ssm
 	 STOP
       ENDIF
       IF(dllname.EQ.'') THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck dllname in namelist ssm\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
          WRITE(*,*) ' Check dllname in namelist ssm'
          WRITE(*,*) ' Program aborting'
 #else
@@ -214,7 +226,9 @@ C FIND namelist S1
 #endif
 	 CLOSE(1)
 	 IF (IFAIL.EQ.-1) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nNamelist S1 not found\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Namelist S1 not found'
         WRITE(*,*) ' Program aborting'
 #else
@@ -234,7 +248,9 @@ C READ namelist S1
 	 READ(1,NML=S1,END=5002,ERR=5002)
 
 	 IF ((dynS1.NE.'I').AND.(dynS1.NE.'M'))THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck dynS1 in namelist S1\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check dynS1 in namelist S1'
         WRITE(*,*) ' Program aborting'
 #else
@@ -246,7 +262,9 @@ C READ namelist S1
 	 ENDIF
 
 	 IF (ns1.LT.2) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck ns1 in namelist S1\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check ns1 in namelist S1'
         WRITE(*,*) ' Program aborting'
 #else
@@ -260,7 +278,9 @@ C READ namelist S1
 	 IF (dynS1.EQ.'I') THEN
         DO J = 1,ns1
 	   IF (hypS1(J,1).LE.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hypS1 in namelist S1\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
           WRITE(*,*) ' Check hypS1 in namelist S1'
           WRITE(*,*) ' Program aborting'
 #else
@@ -275,7 +295,9 @@ C READ namelist S1
 	  DO J = 1,ns1
 	  DO K = 1,ns1
 	   IF (hypS1(J,K).LE.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hypS1 in namelist S1\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
           WRITE(*,*) ' Check hypS1 in namelist S1'
           WRITE(*,*) ' Program aborting'
 #else
@@ -293,7 +315,9 @@ C READ namelist S1
      #	  .AND.(matS1(1).NE.'H').AND.(matS1(1).NE.'G')
      #	  .AND.(matS1(1).NE.'c').AND.(matS1(1).NE.'F')
      #	  .AND.(matS1(1).NE.'R'))) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck matS1 in namelist S1\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check matS1 in namelist S1'
         WRITE(*,*) ' Program aborting'
 #else
@@ -308,7 +332,9 @@ C READ namelist S1
      #	  .AND.(matS1(I).NE.'H').AND.(matS1(I).NE.'G')
      #	  .AND.(matS1(I).NE.'c').AND.(matS1(I).NE.'F')
      #	  .AND.(matS1(I).NE.'R')) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck matS1 in namelist S1\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
          WRITE(*,*) ' Check matS1 in namelist S1'
          WRITE(*,*) ' Program aborting'
 #else
@@ -341,7 +367,9 @@ C FIND namelist S2
 #endif
 	 CLOSE(1)
 	 IF (IFAIL.EQ.-1) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nNamelist S2 not found\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Namelist S2 not found'
         WRITE(*,*) ' Program aborting'
 #else
@@ -360,7 +388,9 @@ C READ namelist S2
 	 READ(1,NML=S2,END=5003,ERR=5003)
 
 	 IF ((dynS2.NE.'I').AND.(dynS2.NE.'M'))THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck dynS2 in namelist S2\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check dynS2 in namelist S2'
         WRITE(*,*) ' Program aborting'
 #else
@@ -372,7 +402,9 @@ C READ namelist S2
 	 ENDIF
 
 	 IF (ns2.LT.2) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck ns2 in namelist S2\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check ns2 in namelist S2'
         WRITE(*,*) ' Program aborting'
 #else
@@ -386,7 +418,9 @@ C READ namelist S2
 	 IF (dynS2.EQ.'I') THEN
         DO J = 1,ns2
 	   IF (hypS2(J,1).LE.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hypS2 in namelist S2\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
           WRITE(*,*) ' Check hypS2 in namelist S2'
           WRITE(*,*) ' Program aborting'
 #else
@@ -401,7 +435,9 @@ C READ namelist S2
 	  DO J = 1,ns2
 	  DO K = 1,ns2
 	   IF (hypS2(J,K).LE.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hypS2 in namelist S2\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
           WRITE(*,*) ' Check hypS2 in namelist S2'
           WRITE(*,*) ' Program aborting'
 #else
@@ -419,7 +455,9 @@ C READ namelist S2
      #	  .AND.(matS2(I).NE.'H').AND.(matS2(I).NE.'G')
      #	  .AND.(matS2(I).NE.'c').AND.(matS2(I).NE.'F')
      #	  .AND.(matS2(I).NE.'R'))) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck matS2 in namelist S2\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check matS2 in namelist S2'
         WRITE(*,*) ' Program aborting'
 #else
@@ -434,7 +472,9 @@ C READ namelist S2
      #	  .AND.(matS2(I).NE.'H').AND.(matS2(I).NE.'G')
      #	  .AND.(matS2(I).NE.'c').AND.(matS2(I).NE.'F')
      #	  .AND.(matS2(I).NE.'R')) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck matS2 in namelist S2\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Check matS2 in namelist S2'
        WRITE(*,*) ' Program aborting'
 #else
@@ -466,7 +506,9 @@ C FIND namelist S3
 #endif
 	 CLOSE(1)
 	 IF (IFAIL.EQ.-1) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck S3 in namelist S3\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Namelist S3 not found'
         WRITE(*,*) ' Program aborting'
 #else
@@ -486,7 +528,9 @@ C READ namelist S3
 	 READ(1,NML=S3,END=5004,ERR=5004)
 
 	 IF ((dynS3.NE.'I').AND.(dynS3.NE.'M'))THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck dynS3 in namelist S3\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check dynS3 in namelist S3'
         WRITE(*,*) ' Program aborting'
 #else
@@ -498,7 +542,9 @@ C READ namelist S3
 	 ENDIF
 
 	 IF (ns3.LT.2) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck ns3 in namelist S3\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check ns3 in namelist S3'
         WRITE(*,*) ' Program aborting'
 #else
@@ -512,7 +558,9 @@ C READ namelist S3
 	 IF (dynS3.EQ.'I') THEN
         DO J = 1,ns3
 	   IF (hypS3(J,1).LE.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hypS3 in namelist S3\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
           WRITE(*,*) ' Check hypS3 in namelist S3'
           WRITE(*,*) ' Program aborting'
 #else
@@ -527,7 +575,9 @@ C READ namelist S3
 	  DO J = 1,ns3
 	  DO K = 1,ns3
 	   IF (hypS3(J,K).LE.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hypS3 in namelist S3\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
           WRITE(*,*) ' Check hypS3 in namelist S3'
           WRITE(*,*) ' Program aborting'
 #else
@@ -545,7 +595,9 @@ C READ namelist S3
      #	  .AND.(matS3(1).NE.'H').AND.(matS3(1).NE.'G')
      #	  .AND.(matS3(1).NE.'c').AND.(matS3(1).NE.'F')
      #	  .AND.(matS3(1).NE.'R'))) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck matS3 in namelist S3\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check matS3 in namelist S3'
         WRITE(*,*) ' Program aborting'
 #else
@@ -560,7 +612,9 @@ C READ namelist S3
      #	  .AND.(matS3(I).NE.'H').AND.(matS3(I).NE.'G')
      #	  .AND.(matS3(I).NE.'c').AND.(matS3(I).NE.'F')
      #	  .AND.(matS3(I).NE.'R')) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck matS3 in namelist S3\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
          WRITE(*,*) ' Check matS3 in namelist S3'
          WRITE(*,*) ' Program aborting'
 #else
@@ -592,11 +646,13 @@ C FIND namelist S4
 #endif
 	 CLOSE(1)
 	 IF (IFAIL.EQ.-1) THEN
-#ifdef __GFORTRAN__
-        WRITE(*,*) ' Namlist S4 not found'
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nNamelist S4 not found\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
+        WRITE(*,*) ' Namelist S4 not found'
         WRITE(*,*) ' Program aborting'
 #else
-	  TYPE *, ' Namlist S4 not found'
+	  TYPE *, ' Namelist S4 not found'
 	  TYPE *, ' Program aborting'
 	  PAUSE
 #endif
@@ -612,7 +668,9 @@ C READ namelist S4
 	 READ(1,NML=S4,END=5005,ERR=5005)
 
 	 IF ((dynS4.NE.'I').AND.(dynS4.NE.'M'))THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck dynS4 in namelist S4\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check dynS4 in namelist S4'
         WRITE(*,*) ' Program aborting'
 #else
@@ -624,7 +682,9 @@ C READ namelist S4
 	 ENDIF
 
 	 IF (ns4.LT.2) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck ns4 in namelist S4\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check ns4 in namelist S4'
         WRITE(*,*) ' Program aborting'
 #else
@@ -638,7 +698,9 @@ C READ namelist S4
 	 IF (dynS4.EQ.'I') THEN
         DO J = 1,ns4
 	   IF (hypS4(J,1).LE.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hypS4 in namelist S4\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
           WRITE(*,*) ' Check hypS4 in namelist S4'
           WRITE(*,*) ' Program aborting'
 #else
@@ -653,7 +715,9 @@ C READ namelist S4
 	  DO J = 1,ns4
 	  DO K = 1,ns4
 	   IF (hypS4(J,K).LE.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hypS4 in namelist S4\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
           WRITE(*,*) ' Check hypS4 in namelist S4'
           WRITE(*,*) ' Program aborting'
 #else
@@ -671,7 +735,9 @@ C READ namelist S4
      #	  .AND.(matS4(1).NE.'H').AND.(matS4(1).NE.'G')
      #	  .AND.(matS4(1).NE.'c').AND.(matS4(1).NE.'F')
      #	  .AND.(matS4(1).NE.'R'))) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck matS4 in namelist S4\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check matS4 in namelist S4'
         WRITE(*,*) ' Program aborting'
 #else
@@ -686,7 +752,9 @@ C READ namelist S4
      #	  .AND.(matS4(I).NE.'H').AND.(matS4(I).NE.'G')
      #	  .AND.(matS4(I).NE.'c').AND.(matS4(I).NE.'F')
      #	  .AND.(matS4(I).NE.'R')) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck matS4 in namelist S4\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
          WRITE(*,*) ' Check matS4 in namelist S4'
          WRITE(*,*) ' Program aborting'
 #else
@@ -718,11 +786,13 @@ C FIND namelist S5
 #endif
 	 CLOSE(1)
 	 IF (IFAIL.EQ.-1) THEN
-#ifdef __GFORTRAN__
-        WRITE(*,*) ' Namlist S5 not found'
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nNamelist S5 not found\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
+        WRITE(*,*) ' Namelist S5 not found'
         WRITE(*,*) ' Program aborting'
 #else
-	  TYPE *, ' Namlist S5 not found'
+	  TYPE *, ' Namelist S5 not found'
 	  TYPE *, ' Program aborting'
 	  PAUSE
 #endif
@@ -738,7 +808,9 @@ C READ namelist S5
 	 READ(1,NML=S5,END=5006,ERR=5006)
 
 	 IF ((dynS5.NE.'I').AND.(dynS5.NE.'M'))THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck dynS5 in namelist S5\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check dynS5 in namelist S5'
         WRITE(*,*) ' Program aborting'
 #else
@@ -750,7 +822,9 @@ C READ namelist S5
 	 ENDIF
 
 	 IF (ns5.LT.2) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck ns5 in namelist S5\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check ns5 in namelist S5'
         WRITE(*,*) ' Program aborting'
 #else
@@ -764,7 +838,9 @@ C READ namelist S5
 	 IF (dynS5.EQ.'I') THEN
         DO J = 1,ns5
 	   IF (hypS5(J,1).LE.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hypS5 in namelist S5\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
           WRITE(*,*) ' Check hypS5 in namelist S5'
           WRITE(*,*) ' Program aborting'
 #else
@@ -779,7 +855,9 @@ C READ namelist S5
 	  DO J = 1,ns5
 	  DO K = 1,ns5
 	   IF (hypS5(J,K).LE.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hypS5 in namelist S5\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
           WRITE(*,*) ' Check hypS5 in namelist S5'
           WRITE(*,*) ' Program aborting'
 #else
@@ -797,7 +875,9 @@ C READ namelist S5
      #	  .AND.(matS5(1).NE.'H').AND.(matS5(1).NE.'G')
      #	  .AND.(matS5(1).NE.'c').AND.(matS5(1).NE.'F')
      #	  .AND.(matS5(1).NE.'R'))) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck matS5 in namelist S5\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check matS5 in namelist S5'
         WRITE(*,*) ' Program aborting'
 #else
@@ -812,7 +892,9 @@ C READ namelist S5
      #	  .AND.(matS5(I).NE.'H').AND.(matS5(I).NE.'G')
      #	  .AND.(matS5(I).NE.'c').AND.(matS5(I).NE.'F')
      #	  .AND.(matS5(I).NE.'R')) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck matS5 in namelist S5\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
          WRITE(*,*) ' Check matS5 in namelist S5'
          WRITE(*,*) ' Program aborting'
 #else
@@ -844,11 +926,13 @@ C FIND namelist S6
 #endif
 	 CLOSE(1)
 	 IF (IFAIL.EQ.-1) THEN
-#ifdef __GFORTRAN__
-        WRITE(*,*) ' Namlist S6 not found'
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nNamelist S6 not found\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
+        WRITE(*,*) ' Namelist S6 not found'
         WRITE(*,*) ' Program aborting'
 #else
-	  TYPE *, ' Namlist S6 not found'
+	  TYPE *, ' Namelist S6 not found'
 	  TYPE *, ' Program aborting'
 	  PAUSE
 #endif
@@ -864,7 +948,9 @@ C READ namelist S6
 	 READ(1,NML=S6,END=5007,ERR=5007)
 
 	 IF ((dynS6.NE.'I').AND.(dynS6.NE.'M'))THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck dynS6 in namelist S6\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check dynS6 in namelist S6'
         WRITE(*,*) ' Program aborting'
 #else
@@ -876,7 +962,9 @@ C READ namelist S6
 	 ENDIF
 
 	 IF (ns6.LT.2) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck ns6 in namelist S6\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check ns6 in namelist S6'
         WRITE(*,*) ' Program aborting'
 #else
@@ -890,7 +978,9 @@ C READ namelist S6
 	 IF (dynS6.EQ.'I') THEN
         DO J = 1,ns6
 	   IF (hypS6(J,1).LE.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hypS6 in namelist S6\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
           WRITE(*,*) ' Check hypS6 in namelist S6'
           WRITE(*,*) ' Program aborting'
 #else
@@ -905,7 +995,9 @@ C READ namelist S6
 	  DO J = 1,ns6
 	  DO K = 1,ns6
 	   IF (hypS6(J,K).LE.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hypS6 in namelist S6\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
           WRITE(*,*) ' Check hypS6 in namelist S6'
           WRITE(*,*) ' Program aborting'
 #else
@@ -923,7 +1015,9 @@ C READ namelist S6
      #	  .AND.(matS6(1).NE.'H').AND.(matS6(1).NE.'G')
      #	  .AND.(matS6(1).NE.'c').AND.(matS6(1).NE.'F')
      #	  .AND.(matS6(1).NE.'R'))) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck matS6 in namelist S6\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
         WRITE(*,*) ' Check matS6 in namelist S6'
         WRITE(*,*) ' Program aborting'
 #else
@@ -938,7 +1032,9 @@ C READ namelist S6
      #	  .AND.(matS6(I).NE.'H').AND.(matS6(I).NE.'G')
      #	  .AND.(matS6(I).NE.'c').AND.(matS6(I).NE.'F')
      #	  .AND.(matS6(I).NE.'R')) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck matS6 in namelist S6\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
          WRITE(*,*) ' Check matS6 in namelist S6'
          WRITE(*,*) ' Program aborting'
 #else
@@ -969,7 +1065,9 @@ C FIND namelist prior
 #endif
 	CLOSE(1)
 	IF (IFAIL.EQ.-1) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nNamelist prior not found\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Namelist prior not found'
        WRITE(*,*) ' Program aborting'
 #else
@@ -986,7 +1084,9 @@ C READ namelist prior
 	hyptheta(:,:) = -1
 	READ(1,NML = prior,END=5008,ERR=5008)
 	IF (nt.LE.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck nt in namelist prior\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Check nt in namelist prior'
        WRITE(*,*) ' Program aborting'
 #else
@@ -997,7 +1097,9 @@ C READ namelist prior
 	 STOP
       ENDIF
 	IF (nt.GT.200) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nnt is too large\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' nt is too large '
        WRITE(*,*) ' Program aborting'
 #else
@@ -1012,7 +1114,9 @@ C READ namelist prior
 	  WRITE(IC,'(I3)') I
 	  IF ((pdftheta(I).NE.'BE').AND.(pdftheta(I).NE.'NT').AND.
      #     (pdftheta(I).NE.'IG')) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck pdftheta('//IC//') in namelist prior\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
          WRITE(*,*) ' Check pdftheta('//IC//') in namelist prior'
          WRITE(*,*) ' Program aborting'
 #else
@@ -1023,7 +1127,9 @@ C READ namelist prior
 	   STOP
 	  ENDIF
 	  IF (hyptheta(3,I).GT.hyptheta(4,I)) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hyptheta('//IC//') in namelist prior\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
          WRITE(*,*) ' Check hyptheta('//IC//') in namelist prior'
          WRITE(*,*) ' Program aborting'
 #else
@@ -1037,7 +1143,9 @@ C READ namelist prior
 	   IF (hyptheta(3,I).LT.hyptheta(4,I)) THEN
 	    IF ((hyptheta(1,I).LE.0.).OR.(hyptheta(2,I).LE.0.).OR.
      #       (hyptheta(3,I).GT.hyptheta(4,I))) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hyptheta('//IC//') in namelist prior\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
            WRITE(*,*) ' Check hyptheta('//IC//') in namelist prior'
            WRITE(*,*) ' Program aborting'
 #else
@@ -1052,7 +1160,9 @@ C READ namelist prior
 	   IF (hyptheta(3,I).LT.hyptheta(4,I)) THEN
 	    IF ((hyptheta(2,I).LE.0.).OR.(hyptheta(3,I).GT.hyptheta(4,I)))
      #     THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hyptheta('//IC//') in namelist prior\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
            WRITE(*,*) ' Check hyptheta('//IC//') in namelist prior'
            WRITE(*,*) ' Program aborting'
 #else
@@ -1067,7 +1177,9 @@ C READ namelist prior
  	   IF (hyptheta(3,I).LT.hyptheta(4,I)) THEN
 	    IF ((hyptheta(1,I).LE.0.).OR.(hyptheta(2,I).LE.0.).OR.
      #    (hyptheta(3,I).GT.hyptheta(4,I)).OR.(hyptheta(3,I).LT.0.))THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hyptheta('//IC//') in namelist prior\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
            WRITE(*,*) ' Check hyptheta('//IC//') in namelist prior'
            WRITE(*,*) ' Program aborting'
 #else
@@ -1084,7 +1196,9 @@ C READ namelist prior
        DO I = 1,nt  ! ML check
         WRITE(IC,'(I3)') I
         IF (hyptheta(3,I).GT.hyptheta(4,I)) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hyptheta('//IC//') in namelist prior\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
            WRITE(*,*) ' Check hyptheta('//IC//') in namelist prior'
            WRITE(*,*) ' Program aborting'
 #else
@@ -1115,7 +1229,9 @@ C FIND namelist mcmc
 #endif
 	CLOSE(1)
 	IF (IFAIL.EQ.-1) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nNamelist mcmc not found\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Namelist mcmc not found'
        WRITE(*,*) ' Program aborting'
 #else
@@ -1138,7 +1254,9 @@ C READ namelist mcmc
 	MargLik      = 'N'
 	READ(1,NML=mcmc,END=5009,ERR=5009)
 	IF ((seed.LT.0).OR.(seed.GT.999)) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck seed in namelist mcmc\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Check seed in namelist mcmc'
        WRITE(*,*) ' Program aborting'
 #else
@@ -1149,7 +1267,9 @@ C READ namelist mcmc
 	 STOP
 	ENDIF
 	IF (thin.LT.1) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck thin in namelist mcmc\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Check thin in namelist mcmc'
        WRITE(*,*) ' Program aborting'
 #else
@@ -1160,7 +1280,9 @@ C READ namelist mcmc
 	 STOP
 	ENDIF
 	IF (burnin.LE.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck burnin in namelist mcmc\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Check burnin in namelist mcmc'
        WRITE(*,*) ' Program aborting'
 #else
@@ -1171,7 +1293,9 @@ C READ namelist mcmc
 	 STOP
 	ENDIF
 	IF (simulrec.LE.1) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck simulrec in namelist mcmc\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Check simulrec in namelist mcmc'
        WRITE(*,*) ' Program aborting'
 #else
@@ -1182,7 +1306,9 @@ C READ namelist mcmc
 	 STOP
 	ENDIF
 	IF ((sampler.NE.'SL').AND.(sampler.NE.'MH')) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck sampler in namelist mcmc\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Check sampler in namelist mcmc'
        WRITE(*,*) ' Program aborting'
 #else
@@ -1226,7 +1352,9 @@ C FIND namelist dataset
 #endif
 	CLOSE(1)
 	IF (IFAIL.EQ.-1) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nNamelist dataset not found\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Namelist dataset not found'
        WRITE(*,*) ' Program aborting'
 #else
@@ -1244,7 +1372,9 @@ C READ namelist dataset
 	datasim = 'N'
 	READ(1,NML=dataset,END=5010,ERR=5010)
 	IF ((T.LE.0).OR.(T.GT.3000)) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck T in namelist dataset\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Check T in namelist dataset (T<=3000)'
        WRITE(*,*) ' Program aborting'
 #else
@@ -1255,7 +1385,9 @@ C READ namelist dataset
 	 STOP
 	ENDIF
 	IF (ny.LE.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck ny in namelist dataset\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Check ny in namelist dataset'
        WRITE(*,*) ' Program aborting'
 #else
@@ -1266,7 +1398,9 @@ C READ namelist dataset
 	 STOP
 	ENDIF
 	IF (nz.LT.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck nz in namelist dataset\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Check nz in namelist dataset'
        WRITE(*,*) ' Program aborting'
 #else
@@ -1277,7 +1411,9 @@ C READ namelist dataset
 	 STOP
 	ENDIF
 	IF (nf.LT.0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck nf in namelist dataset\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Check nf in namelist dataset'
        WRITE(*,*) ' Program aborting'
 #else
@@ -1288,7 +1424,9 @@ C READ namelist dataset
 	 STOP
 	ENDIF
 	IF (T.LT.hbl) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nCheck hbl in namelist mcmc (hbl > T)\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' Check hbl in namelist mcmc (hbl > T)'
        WRITE(*,*) ' Program aborting'
 #else
@@ -1406,8 +1544,19 @@ C -----------------------------------------------------------------------
 	ENDIF
 
 	GO TO 7777
-
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+ 5000 CALL mexErrMsgTxt('\nInput file not found\nProgram aborting\n')
+ 5001 CALL mexErrMsgTxt('\nInput error in namelist ssm\nProgram aborting\n')
+ 5002 CALL mexErrMsgTxt('\nInput error in namelist S1\nProgram aborting\n')
+ 5003 CALL mexErrMsgTxt('\nInput error in namelist S2\nProgram aborting\n')
+ 5004 CALL mexErrMsgTxt('\nInput error in namelist S3\nProgram aborting\n')
+ 5005 CALL mexErrMsgTxt('\nInput error in namelist S4\nProgram aborting\n')
+ 5006 CALL mexErrMsgTxt('\nInput error in namelist S5\nProgram aborting\n')
+ 5007 CALL mexErrMsgTxt('\nInput error in namelist S6\nProgram aborting\n')
+ 5008 CALL mexErrMsgTxt('\nInput error in namelist prior\nProgram aborting\n')
+ 5009 CALL mexErrMsgTxt('\nInput error in namelist mcmc\nProgram aborting\n')
+ 5010 CALL mexErrMsgTxt('\nInput error in namelist dataset\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
  5000 WRITE(*,*)'Input file not found'
       WRITE(*,*)'Program aborting'
       STOP
diff --git a/kim.for b/kim.for
index 6cd5067f32dbe97127d880efa35f89c4dfb353f2..1b27eed70d2dfbad4ddd167fba8120380c180b8f 100644
--- a/kim.for
+++ b/kim.for
@@ -136,7 +136,10 @@ C X-filter initialization
 
         ENDDO
       ELSE
-        WRITE(*,*) 'WARNING: d(1)=2 not implemeted yet'
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nERROR: d(1)=2 not implemeted yet\n')
+#endif
+        WRITE(*,*) 'ERROR: d(1)=2 not implemeted yet'
         PAUSE
         STOP
       ENDIF
diff --git a/lyapunov.for b/lyapunov.for
index 59a745e344c5d5b0862e2b09db1b0c24fe2c9a7e..f271a9cad9cbb010ec30afebb877553c5cd97c31 100644
--- a/lyapunov.for
+++ b/lyapunov.for
@@ -63,7 +63,9 @@ C	CALL F02EAF('V',nx,T,nx,WR,WI,Z,nx,WORK,LWORK,IFAIL) ! F = ZTZ'
      #           BWORK,IFAIL)
 	DO I = 1,nx
 	IF (WI(I)**2+WR(I)**2.GE.1.D0) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nLYAPUNOV SUBROUTINE: Some parameters out of\nstationary region. Check hyptheta in namelist prior.\nProgram aborting')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' '
        WRITE(*,*) ' LYAPUNOV SUBROUTINE: Some parameters out of '
        WRITE(*,*) ' stationary region. Check hyptheta in namelist '
diff --git a/main.for b/main.for
index 34e0c5d40fec2f471a34433f9f0dd583ae88af9f..586203724fe2d071de2f7cc401b4526748c2e63d 100644
--- a/main.for
+++ b/main.for
@@ -34,7 +34,12 @@ C
 C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C --------------------------------------------------------------------------------------
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+#include "fintrf.h"
+      SUBROUTINE DMMMAIN(FILEIN)
+#else
       PROGRAM DMM
+#endif
 #if defined(__CYGWIN32__) || defined(_WIN32)
 #ifdef __INTEL_COMPILER
       USE dfwin
@@ -43,46 +48,6 @@ C ------------------------------------------------------------------------------
       USE ISO_C_BINDING
       USE ISO_C_UTILITIES
       USE DLFCN
-#endif
-C DECLARE an "interface block" to the .DLL that contains DESIGN
-
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-	CHARACTER*1 fittizia
-C DECLARE an "interface block" to the .DLL that contains SETFILEM
-	INTERFACE
-	 SUBROUTINE SETFILEM(mfile,pathmfile)
-	 CHARACTER*200 mfile,pathmfile
-       END SUBROUTINE
-      END INTERFACE
-C DECLARE an "interface block" to the .DLL that contains GETERRSTR
-	INTERFACE
-	 SUBROUTINE GETERRSTR(matlaberror)
-	 CHARACTER*1024 matlaberror
-       END SUBROUTINE
-      END INTERFACE
-
-#if defined(__CYGWIN32__) || defined(_WIN32)
-      LOGICAL STATUS
-      POINTER (pdll,fittizia)   ! ASSOCIATE  pointer pdll alla DLL ad una varibile fittizia
-      POINTER (pdesign,DESIGN)
-      POINTER (psetfilem,SETFILEM)
-      POINTER (pgeterrstr,GETERRSTR)
-#else
-      INTEGER(C_INT) :: STATUS
-      TYPE(C_PTR) :: pdll=C_NULL_PTR
-      TYPE(C_FUNPTR) :: pdesign=C_NULL_FUNPTR
-      TYPE(C_FUNPTR) :: psetfilem=C_NULL_FUNPTR
-      TYPE(C_FUNPTR) :: pgeterrstr=C_NULL_FUNPTR
-      PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-      PROCEDURE(DESIGN), POINTER :: ptrsetfilem=>NULL()
-      PROCEDURE(DESIGN), POINTER :: ptrgeterrstr=>NULL()
 #endif
 	CHARACTER*200 DLLNAME    ! name of the DLL (defined by the user)
 
@@ -110,14 +75,20 @@ C LOCALS
 	INTEGER ntf,nstot,nmis,indmis,IT,I,J,K,L1,jjj,IND,IFAIL,IMAX(1),
      1 IMIN(1),IMSVAR
 	DOUBLE PRECISION AUX,lastl,lasth
-	CHARACTER*1 DEB
       CHARACTER*3 DLLEXT
       CHARACTER*200 mfile,pathmfile
-	CHARACTER*200 FILEIN,NMLNAME,PATH,FILEOUT,DMMTITLE,CURDIR
-      CHARACTER*1024 matlaberror
+	CHARACTER*200 FILEIN,NMLNAME,PATH,FILEOUT,CURDIR
+#ifndef __GFORTRAN__
+      CHARACTER*200 DMMTITLE
+#endif
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+      CHARACTER(len=200) :: MEXPRINT
+      INTEGER*4 mexPrintf
+      INTEGER*4 mpfout
+#endif
 
 #ifdef __GFORTRAN__
-      CHARACTER*12 fmt
+      CHARACTER*16 fmt
 #endif
 
 #ifdef __INTEL_COMPILER
@@ -134,19 +105,15 @@ C TIME
      1                   DATE_ITIME)
       IT1(1:3) = DATE_ITIME(1:3)
       IT1(4:7) = DATE_ITIME(5:8)
-
+#if !(defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE))
 C GET the namelist specified by FILEIN
-	DEB = 'D'
-	IF (DEB.EQ.'R') THEN
-       CALL GETARG(1,FILEIN)    ! load name of input file
-      ELSE
-      FILEIN = 'H:\AROSSI\DMM\NILE\nile.nml'
-C     FILEIN = 'H:\arossi\dmm\tfpf\tfpf_es.nml'
-	ENDIF
-
+      CALL GETARG(1,FILEIN)     ! load name of input file
+#endif
 C CHECK FILEIN
 	IF (TRIM(FILEIN).EQ.'') THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nNo input file provided\nProgram Aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' '
        WRITE(*,*) ' No input file provided'
        WRITE(*,*) ' Program aborting'
@@ -174,242 +141,32 @@ C CHECK DLL NAME AND FIND FILE EXTENSION (.dll or .m)
       IF ((DLLEXT.EQ.'M  ').OR.(DLLEXT.EQ.'m  ')) THEN
        mfile     = DLLNAME(J+1:I-1)
        pathmfile = DLLNAME(1:J-1)
-       DLLNAME   = 'H:\arossi\dmm64\matlabdll\debug\matlabdll.dll' ! provvisorio
        IND = GETCWD(CURDIR)  ! current directory
-C       DLLNAME = TRIM(CURDIR) // '\matlabdll.dll'                 ! definitivo
-      ENDIF
-
-C FIND the DLL and LOAD it into the memory
-#if defined(__CYGWIN32__) || defined(_WIN32)
-      pdll = loadlibrary(DLLNAME)
-      IF (pdll.EQ.0) THEN
-         TYPE *, ' '
-         TYPE *, TRIM(DLLNAME) // ' cannot be found or opened'
-         TYPE *, ' Program aborting'
-         PAUSE
-#else
-      pdll = DLOpen(TRIM(DLLNAME)//C_NULL_CHAR, RTLD_NOW)
-      IF(.NOT.C_ASSOCIATED(pdll)) THEN
-         WRITE(*,*) ' '
-         WRITE(*,*) ' Error in dlopen: ', C_F_STRING(DLError())
-         WRITE(*,*) TRIM(DLLNAME) // ' cannot be found or opened'
-         WRITE(*,*) ' Program aborting'
-#endif
-         STOP
-      ENDIF
-
-C SET UP the pointer to the DLL function
-#if defined(__CYGWIN32__) || defined(_WIN32)
-      pdesign = getprocaddress(pdll, "design_"C)
-      IF (pdesign.EQ.0) THEN
-         TYPE *, ' '
-         TYPE *, ' Sub DESIGN cannot be found into '// DLLNAME
-         TYPE *, ' Program aborting'
-         PAUSE
-         STOP
-      ENDIF
-#else
-      pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-      IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-         WRITE(*,*) ' '
-         WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-         WRITE(*,*) ' Sub DESIGN cannot be found into '// DLLNAME
-         WRITE(*,*) ' Program aborting'
-         STOP
       ENDIF
-      CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-#endif
 
 C CHECK the MatLab file if needed
       IF ((DLLEXT.EQ.'M  ').OR.(DLLEXT.EQ.'m  ')) THEN
-C SET UP the pointer to the DLL function
-#if defined(__CYGWIN32__) || defined(_WIN32)
-         psetfilem = getprocaddress(pdll, "setfilem_"C)
-         IF (psetfilem.EQ.0) THEN
-            TYPE *, ' '
-            TYPE *, ' Sub SETFILEM cannot be found into '// DLLNAME
-            TYPE *, ' Program aborting'
-            PAUSE
-            STOP
-         ENDIF
-#else
-         psetfilem = DLSym(pdll, 'setfilem_'//C_NULL_CHAR)
-         IF(.NOT.C_ASSOCIATED(psetfilem)) THEN
-            WRITE(*,*) ' '
-            WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-            WRITE(*,*) ' Sub SETFILEM cannot be found into '// DLLNAME
-            WRITE(*,*) ' Program aborting'
-            STOP
-         ENDIF
-         CALL C_F_PROCPOINTER(CPTR=psetfilem, FPTR=ptrsetfilem)
-#endif
-
-C SET UP the pointer to the DLL function
-#if defined(__CYGWIN32__) || defined(_WIN32)
-         pgeterrstr = getprocaddress(pdll, "geterrstr_"C)
-         IF (pgeterrstr.EQ.0) THEN
-            TYPE *, ' '
-            TYPE *, ' Sub GETERRSTR cannot be found into '// DLLNAME
-            TYPE *, ' Program aborting'
-            PAUSE
-            STOP
-         ENDIF
-#else
-         pgeterrstr = DLSym(pdll, 'geterrstr_'//C_NULL_CHAR)
-         IF(.NOT.C_ASSOCIATED(pgeterrstr)) THEN
-            WRITE(*,*) ' '
-            WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-            WRITE(*,*) ' Sub GETERRSTR cannot be found into '// DLLNAME
-            WRITE(*,*) ' Program aborting'
-            STOP
-         ENDIF
-         CALL C_F_PROCPOINTER(CPTR=pgeterrstr, FPTR=ptrgeterrstr)
-#endif
-
 C Assign the name of the matlab file
        ALLOCATE( c(ny,max(nz,1),ns(1)),H(ny,nx,ns(2)),
 	1  G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6)),
      1  theta(nt))
-       CALL SETFILEM(mfile,pathmfile)  ! ONLY THE FIRST TIME
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL SETFILEM(mfile,pathmfile) ! ONLY THE FIRST TIME
+#endif
+
        theta(:) = 1.D0
-#if defined(__CYGWIN32__) || defined(_WIN32)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
        CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-#else
-       CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
        DEALLOCATE(c,H,G,a,F,R,theta)
-       IF (ny.EQ.0) THEN
-#ifdef __GFORTRAN__
-          WRITE(*,*) ' '
-          WRITE(*,*) ' Can''t start MATLAB engine'
-          WRITE(*,*) ' Program aborting'
-#else
-        TYPE *, ' '
-	  TYPE *, ' Can''t start MATLAB engine'
-	  TYPE *, ' Program aborting'
-	  PAUSE
-#endif
-	  STOP
-       ELSEIF (ny.EQ.-1) THEN
-#ifdef __GFORTRAN__
-          WRITE(*,*) ' '
-          WRITE(*,*) ' Can''t read ny in the MATLAB file'
-          WRITE(*,*) ' Program aborting'
-#else
-        TYPE *, ' '
-	  TYPE *, ' Can''t read ny in the MATLAB file'
-	  TYPE *, ' Program aborting'
-	  PAUSE
-#endif
-	  STOP
-       ELSEIF (ny.EQ.-2) THEN
-#ifdef __GFORTRAN__
-          WRITE(*,*) ' '
-          WRITE(*,*) ' Can''t read nz in the MATLAB file'
-          WRITE(*,*) ' Program aborting'
-#else
-        TYPE *, ' '
-	  TYPE *, ' Can''t read nz in the MATLAB file'
-	  TYPE *, ' Program aborting'
-	  PAUSE
-#endif
-	  STOP
-       ELSEIF (ny.EQ.-3) THEN
-#ifdef __GFORTRAN__
-          WRITE(*,*) ' '
-          WRITE(*,*) ' Can''t read nx in the MATLAB file'
-          WRITE(*,*) ' Program aborting'
-#else
-        TYPE *, ' '
-	  TYPE *, ' Can''t read nx in the MATLAB file'
-	  TYPE *, ' Program aborting'
-	  PAUSE
-#endif
-	  STOP
-       ELSEIF (ny.EQ.-4) THEN
-#ifdef __GFORTRAN__
-          WRITE(*,*) ' '
-          WRITE(*,*) ' Can''t read nu in the MATLAB file'
-          WRITE(*,*) ' Program aborting'
-#else
-        TYPE *, ' '
-	  TYPE *, ' Can''t read nu in the MATLAB file'
-	  TYPE *, ' Program aborting'
-	  PAUSE
-#endif
-	  STOP
-       ELSEIF (ny.EQ.-5) THEN
-#ifdef __GFORTRAN__
-          WRITE(*,*) ' '
-          WRITE(*,*) ' Can''t read ns in the MATLAB file'
-          WRITE(*,*) ' Program aborting'
-#else
-        TYPE *, ' '
-	  TYPE *, ' Can''t read ns in the MATLAB file'
-	  TYPE *, ' Program aborting'
-	  PAUSE
-#endif
-	  STOP
-       ELSEIF (ny.EQ.-6) THEN
-#ifdef __GFORTRAN__
-          WRITE(*,*) ' '
-          WRITE(*,*) ' Can''t read nt in the MATLAB file'
-          WRITE(*,*) ' Program aborting'
-#else
-        TYPE *, ' '
-	  TYPE *, ' Can''t read nt in the MATLAB file'
-	  TYPE *, ' Program aborting'
-	  PAUSE
-#endif
-	  STOP
-       ELSEIF (ny.EQ.-7) THEN
-#ifdef __GFORTRAN__
-          WRITE(*,*) ' '
-          WRITE(*,*) ' Can''t find or open the MatLab function'
-          WRITE(*,*) ' Program aborting'
-#else
-        TYPE *, ' '
-	  TYPE *, ' Can''t find or open the MatLab function'
-	  TYPE *, ' Program aborting'
-	  PAUSE
-#endif
-	  STOP
-       ELSEIF (ny.EQ.-8) THEN
-        CALL GETERRSTR(matlaberror)
-#ifdef __GFORTRAN__
-        WRITE(*,*) ' '
-        WRITE(*,*) ' the MATLAB funtion can not be executed:'
-        WRITE(*,*) trim(matlaberror)
-        WRITE(*,*) ' Program aborting'
-#else
-        TYPE *, ' '
-	  TYPE *, ' the MATLAB funtion can not be executed:'
-        TYPE *, trim(matlaberror)
-	  TYPE *, ' Program aborting'
-	  PAUSE
-#endif
-	  STOP
-       ELSEIF (ny.LT.-100) THEN
-#ifdef __GFORTRAN__
-          WRITE(*,*) ' '
-          WRITE(*,*) ' One of the output canot be assigned during the call'
-          WRITE(*,*) ' ' // trim(DLLNAME)
-          WRITE(*,*) ' Program aborting'
-#else
-        TYPE *, ' '
-	  TYPE *, ' One of the output canot be assigned during the call '
-        TYPE *, ' ' // trim(DLLNAME)
-        TYPE *, ' Program aborting'
-	  PAUSE
-#endif
-	  STOP
-       ENDIF
       ENDIF
 
 C SET SHELL title
+#ifndef __GFORTRAN__
 	DMMTITLE = 'title DMM input:' // TRIM(PATH) // TRIM(NMLNAME)
-     #     // '.nml' // ' - ' // TRIM(DLLNAME)
+     #     // '.nml' // ' - '
 	CALL system(DMMTITLE)
+#endif
 
 C INITIALISE THE RANDOM NUMBER GENERATOR
       CALL INITRAND(SEED,DATE_ITIME)
@@ -519,7 +276,7 @@ C WRITE HYPERPARAMTERS for THETA and PSI plus DATA
 	FILEOUT = TRIM(PATH)//TRIM(NMLNAME)//'.PRI'
  	OPEN(10,FILE = FILEOUT, ACCESS='SEQUENTIAL')
 #ifdef __GFORTRAN__
-      WRITE(fmt, '(a,i4,a)') '(', nv+11, '(I6))'
+      WRITE(fmt, '(a,i4,a)') '(', nv+11, 'I6)'
       WRITE(10,fmt) nt,np(1:3),nf,nz,seed,nx,ny,nobs,nv,INFOS(8,1:nv)
 #else
 	 WRITE(10,'(<11+nv>(I6))') nt,np(1:3),nf,nz,seed,nx,ny,nobs,
@@ -528,8 +285,8 @@ C WRITE HYPERPARAMTERS for THETA and PSI plus DATA
        WRITE(10,'(A2)') estimation
        DO I =1,nt
 #ifdef __GFORTRAN__
-          WRITE(fmt, '(a,i4,a)') '(', 4, '(F25.12)), ''  '',A2'
-          WRITE(10,fmt) thetaprior(I,1:4),pdftheta(I)
+      WRITE(10, '(4F25.12)', advance='no') thetaprior(I,1:4)
+      WRITE(10,'(A,A)') '  ', pdftheta(I)
 #else
 	  WRITE(10,1111) thetaprior(I,1:4),pdftheta(I)
 #endif
@@ -538,8 +295,10 @@ C WRITE HYPERPARAMTERS for THETA and PSI plus DATA
 	 DO I = 1,nv
 	  IF (INFOS(9,I).EQ.1) THEN
 #ifdef __GFORTRAN__
-         WRITE(fmt, '(a,i4,a)') 'I10,(', np(3), '(F25.12)), ''  '',I2'
-         WRITE(10,fmt) INFOS(8,I),psiprior(K+1,:),INFOS(9,I)
+         WRITE(10,'(I10)',advance='no') INFOS(8,I)
+         WRITE(fmt, '(a,i4,a)') '(', np(3), 'F25.12)'
+         WRITE(10, fmt,advance='no') psiprior(K+1,:)
+         WRITE(10, '(A, I2)') '  ', INFOS(9,I)
 #else
 	    WRITE(10,1112) INFOS(8,I),psiprior(K+1,:),INFOS(9,I)
 #endif
@@ -547,8 +306,10 @@ C WRITE HYPERPARAMTERS for THETA and PSI plus DATA
 	  ELSEIF (INFOS(9,I).EQ.2) THEN
 	    DO J = 1,INFOS(8,I)
 #ifdef __GFORTRAN__
-         WRITE(fmt, '(a,i4,a)') 'I10,(', np(3), '(F25.12)), ''  '',I2'
-         WRITE(10,fmt) INFOS(8,I),psiprior(K+1,:),INFOS(9,I)
+         WRITE(10,'(I10)',advance='no') INFOS(8,I)
+         WRITE(fmt, '(a,i4,a)') '(', np(3), 'F25.12)'
+         WRITE(10, fmt,advance='no') psiprior(K+1,:)
+         WRITE(10, '(A, I2)') '  ', INFOS(9,I)
 #else
 	     WRITE(10,1112) INFOS(8,I),psiprior(K+1,:),INFOS(9,I)
 #endif
@@ -558,8 +319,8 @@ C WRITE HYPERPARAMTERS for THETA and PSI plus DATA
        END DO
        DO I =1,nobs+nf
 #ifdef __GFORTRAN__
-      WRITE(fmt, '(a,i4,a)') '(', ny+nz, '(F20.10))'
-      WRITE(10,fmt) yk(I,1:ny+nz)
+          WRITE(fmt, '(a,i4,a)') '(', ny+nz, 'F20.10)'
+          WRITE(10, fmt) yk(I,1:ny+nz)
 #else
 	  WRITE(10,'(<ny+nz>(F20.10))') yk(I,1:ny+nz)
 #endif
@@ -570,7 +331,7 @@ C WRITE HYPERPARAMTERS for THETA and PSI plus DATA
 
 C CHECK DESIGN.dll
 	IF ((check.EQ.'Y').OR.(check.EQ.'y')) THEN
-	 CALL CHECKDESIGN(ny,nz,nx,nu,ns,nt,d,theta0,pdll,PATH,NMLNAME)
+	 CALL CHECKDESIGN(ny,nz,nx,nu,ns,nt,d,theta0,PATH,NMLNAME)
 	 GOTO 7777
 	ENDIF
 
@@ -578,23 +339,23 @@ C SIMULATION of DATA and UNOBSERVABLES
 	IF ((datasim.EQ.'Y').OR.(datasim.EQ.'y')) THEN
 	 CALL OPENFILES(ESTIMATION,SEED,NV,0,0,datasim,MARGLIK,
 	1                PATH,NMLNAME)
-	 CALL SIMDATA(nobs,d,ny,nz,nx,nu,ns,nstot,nt,nv,np,INFOS,pdll,
+	 CALL SIMDATA(nobs,d,ny,nz,nx,nu,ns,nstot,nt,nv,np,INFOS,
      2              theta0,psi0,Z,STATE,yk)
 	 IF (nv.EQ.0) THEN
 	  WRITE(9,'((F25.15))') theta0(1:nt)
 	 ELSE
 	  WRITE(9,'((F25.15))') theta0(1:nt),psi0(1:np(1))
 #ifdef __GFORTRAN__
-      WRITE(fmt, '(a,i4,a)') '(', 1, '(I3))'
+      WRITE(fmt, '(a,i4,a)') '(', 1, 'I3)'
 	  WRITE(11,fmt) Z(:)
 #else
 	  WRITE(11,'(<1>(I3))') Z(:)
 #endif
 	 ENDIF
 #ifdef __GFORTRAN__
-      WRITE(fmt, '(a,i4,a)') '(', nx, '(F20.10))'
+      WRITE(fmt, '(a,i4,a)') '(', nx, 'F20.10)'
       WRITE(10,fmt) (STATE(I,1:nx),I=1,nobs)
-      WRITE(fmt, '(a,i4,a)') '(', nx, '(F20.10))'
+      WRITE(fmt, '(a,i4,a)') '(', nx, 'F20.10)'
       WRITE(15,fmt) (yk(I,1:ny),I=1,nobs)
 #else
 	 WRITE(10,'(<nx>(F20.10))') (STATE(I,1:nx),I=1,nobs)
@@ -610,7 +371,9 @@ C SIMULATION of DATA and UNOBSERVABLES
 C MAXIMUM LIKELIHOOD ESTIMATION
 	IF ((estimation.EQ.'ML').OR.(estimation.EQ.'ml').OR.
      &    (estimation.EQ.'Ml').OR.(estimation.EQ.'mL')) THEN
-#ifdef __GFORTRAN__
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       CALL mexErrMsgTxt('\nMaximum Likelihood inference not allowed\nProgram aborting\n')
+#elif defined(__GFORTRAN__)
        WRITE(*,*) ' '
        WRITE(*,*) ' Maximum Likelihood inference not allowed '
        WRITE(*,*) ' Program aborting'
@@ -628,16 +391,16 @@ c	 CALL ML(nobs,d,ny,nz,nx,nu,nt,nv,ns,np(1),INFOS,pdll,INDT,yk,IYK,S,
 c	1         thetaprior,theta0,psi0,IMSVAR,HESS,AUX)
 	 ALLOCATE(THETASE(nt),AKMSE(nobs,nx),INN(nobs,ny))
 	 IF (nv.EQ.0) THEN
-        CALL OPG(nobs,d,ny,nz,nx,nu,nt,ns,pdll,yk,IYK,S,
+        CALL OPG(nobs,d,ny,nz,nx,nu,nt,ns,yk,IYK,S,
 	1           theta0,thetaprior,HESS,thetase,STATE,AKMSE,INN,IFAIL)
 #ifdef __GFORTRAN__
-        WRITE(fmt, '(a,i4,a)') '(', 2, '(F25.15))'
+        WRITE(fmt, '(a,i4,a)') '(', 2, 'F25.15)'
         WRITE(9,fmt) (theta0(I),thetase(I),I=1,nt)
         WRITE(9,fmt) AUX,IFAIL
-        WRITE(fmt, '(a,i4,a)') '(', nx, '(F20.10))'
+        WRITE(fmt, '(a,i4,a)') '(', nx, 'F20.10)'
         WRITE(10,fmt) (STATE(I,1:nx),I=1,nobs)
         WRITE(10,fmt) (AKMSE(I,1:nx),I=1,nobs)
-        WRITE(fmt, '(a,i4,a)') '(', ny, '(F20.10))'
+        WRITE(fmt, '(a,i4,a)') '(', ny, 'F20.10)'
         WRITE(12,fmt) (INN(I,1:ny),I=1,nobs)
 #else
         WRITE(9,'(<2>(F25.15))') (theta0(I),thetase(I),I=1,nt)
@@ -649,25 +412,25 @@ c	1         thetaprior,theta0,psi0,IMSVAR,HESS,AUX)
        ELSE
         ALLOCATE(psise(np(1)),SSMOOTH(nobs,nstot))
         IF(IMSVAR.EQ.1)THEN
-         CALL OPGH(nobs,ny,nz,nx,nu,nt,nv,ns,nstot,np(1),pdll,yk,IYK,
+         CALL OPGH(nobs,ny,nz,nx,nu,nt,nv,ns,nstot,np(1),yk,IYK,
      1             INFOS,theta0,psi0,thetaprior,HESS,thetase,psise,
      1             SSMOOTH,INN,IFAIL)
         ELSE
-         CALL OPGKIM(nobs,d,ny,nz,nx,nu,nt,nv,ns,nstot,np(1),pdll,
+         CALL OPGKIM(nobs,d,ny,nz,nx,nu,nt,nv,ns,nstot,np(1),
      1               yk,IYK,INFOS,theta0,psi0,thetaprior,HESS,
      1               thetase,psise,STATE,AKMSE,SSMOOTH,INN,IFAIL)
 #ifdef __GFORTRAN__
-         WRITE(fmt, '(a,i4,a)') '(', nx, '(F20.10))'
+         WRITE(fmt, '(a,i4,a)') '(', nx, 'F20.10)'
          WRITE(10,fmt) (STATE(I,1:nx),I=1,nobs)
          WRITE(10,fmt) (AKMSE(I,1:nx),I=1,nobs)
       ENDIF
-      WRITE(fmt, '(a,i4,a)') '(', 2, '(F25.15))'
+      WRITE(fmt, '(a,i4,a)') '(', 2, 'F25.15)'
       WRITE(9,fmt) (theta0(I),thetase(I),I=1,nt)
       WRITE(9,fmt) (psi0(I),psise(I),I=1,np(1))
 	  WRITE(9,fmt) AUX,IFAIL
-      WRITE(fmt, '(a,i4,a)') '(', nstot, '(F20.10))'
+      WRITE(fmt, '(a,i4,a)') '(', nstot, 'F20.10)'
       WRITE(11,fmt) (SSMOOTH(I,1:nstot),I=1,nobs)
-      WRITE(fmt, '(a,i4,a)') '(', ny, '(F20.10))'
+      WRITE(fmt, '(a,i4,a)') '(', ny, 'F20.10)'
       WRITE(12,fmt) (INN(I,1:ny),I=1,nobs)
 #else
 	   WRITE(10,'(<nx>(F20.10))') (STATE(I,1:nx),I=1,nobs)
@@ -707,7 +470,7 @@ C MCMC BURN-IN
 	  IF (nv.GT.0) THEN
 	   CALL GCK(nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np(1),
      1            yk(1:nobs,:),IYK(1:nobs,:),theta0,psi0,
-     2            INFOS,pdll,Z,S)
+     2            INFOS,Z,S)
 	   IF (HBL.GT.1) THEN
 		CALL RECPR(jjj,nstot,nobs,Z,ZW,PM,PTR)
 	   ENDIF
@@ -717,7 +480,7 @@ C MCMC BURN-IN
 	   IF (thetaprior(it,3).LT.thetaprior(it,4)) THEN
 	    CALL SLICE(it,nobs,d,ny,nz,nx,nu,ns,nt,S,
 	1               yk(1:nobs,:),IYK(1:nobs,:),theta0,
-	2               thetaprior(it,:),pdftheta(it),pdll,
+	2               thetaprior(it,:),pdftheta(it),
      3               NEVAL(it),theta(it))
           theta0(it) = theta(it)
          ENDIF
@@ -728,16 +491,27 @@ C MCMC BURN-IN
 	   IMIN(1) = CUMN(INDT(IMIN(1))) !CUMN(IMIN(1))
 	   IMAX    = MAXLOC(CUMN(INDT(1:ntf)))
 	   IMAX(1) = CUMN(INDT(IMAX(1))) !CUMN(IMAX(1))
+#if defined(__CYGWIN32__) || defined(_WIN32)
 	   CALL system('cls')
+#endif
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       WRITE(MEXPRINT,11131) jjj
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+       WRITE(MEXPRINT,11132) ntf
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+       WRITE(MEXPRINT,11133) IMIN(1)/dfloat(jjj),IMAX(1)/dfloat(jjj)
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
 	   WRITE(6,1113) jjj,ntf,IMIN(1)/dfloat(jjj),
      #             IMAX(1)/dfloat(jjj)
+#endif
         ENDIF
        ENDDO
 	ELSE  ! NO MISSING
 	 DO jjj = 1,burnin
 	  IF (nv.GT.0) THEN
 	   CALL GCK2(nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np(1),
-	1             yk(1:nobs,:),theta0,psi0,INFOS,pdll,Z,S)
+	1             yk(1:nobs,:),theta0,psi0,INFOS,Z,S)
 	   IF (HBL.GT.1) THEN
 	    CALL RECPR(jjj,nstot,nobs,Z,ZW,PM,PTR)
 	   ENDIF
@@ -746,7 +520,7 @@ C MCMC BURN-IN
         DO it = 1,nt
 	   IF (thetaprior(it,3).LT.thetaprior(it,4)) THEN
 	    CALL SLICE2(it,nobs,d,ny,nz,nx,nu,ns,nt,S,yk(1:nobs,:),
-	1             theta0,thetaprior(it,:),pdftheta(it),pdll,
+	1             theta0,thetaprior(it,:),pdftheta(it),
      2                NEVAL(it),theta(it))
           theta0(it) = theta(it)
          ENDIF
@@ -757,9 +531,20 @@ C MCMC BURN-IN
 	   IMIN(1) = CUMN(INDT(IMIN(1))) !CUMN(IMIN(1))
 	   IMAX    = MAXLOC(CUMN(INDT(1:ntf)))
 	   IMAX(1) = CUMN(INDT(IMAX(1))) !CUMN(IMAX(1))
+#if defined(__CYGWIN32__) || defined(_WIN32)
 	   CALL system('cls')
+#endif
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       WRITE(MEXPRINT,11131) jjj
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+       WRITE(MEXPRINT,11132) ntf
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+       WRITE(MEXPRINT,11133) IMIN(1)/dfloat(jjj),IMAX(1)/dfloat(jjj)
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
 	   WRITE(6,1113) jjj,ntf,IMIN(1)/dfloat(jjj),
      #             IMAX(1)/dfloat(jjj)
+#endif
         ENDIF
        ENDDO
 	ENDIF
@@ -791,11 +576,11 @@ C MCMC RECORDING phase
 	   IF (HBL.EQ.1) THEN
 	    CALL GCK(nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np(1),
 	1             yk(1:nobs,:),IYK(1:nobs,:),theta0,psi0,
-     2             INFOS,pdll,Z,S)
+     2             INFOS,Z,S)
 	   ELSE
 	    CALL AMH(HBL,nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np(1),
 	1             yk(1:nobs,:),IYK(1:nobs,:),theta0,psi0,
-	2             PTR,PM,INFOS,pdll,Z,S,ACCRATE)
+	2             PTR,PM,INFOS,Z,S,ACCRATE)
 		CALL RECPR(jjj+burnin,nstot,nobs,Z,ZW,PM,PTR)
 	   ENDIF
 	   CALL DRAWPSI(nobs,nv,np,INFOS,Z,psiprior,psi0,psi)
@@ -804,19 +589,19 @@ C MCMC RECORDING phase
 	   IF (thetaprior(it,3).LT.thetaprior(it,4)) THEN
 	    CALL SLICE(it,nobs,d,ny,nz,nx,nu,ns,nt,S,yk(1:nobs,:),
 	1               IYK(1:nobs,:),theta0,thetaprior(it,:),
-     2               pdftheta(it),pdll,NEVAL(it),theta(it))
+     2               pdftheta(it),NEVAL(it),theta(it))
           theta0(it) = theta(it)
          ENDIF
 	  END DO
 	  CUMN = CUMN+NEVAL
 	  CALL SIMSTATE(nobs,d,ny,nz,nx,nu,ns,nt,yk(1:nobs,:),
-	1                IYK(1:nobs,:),theta,S,pdll,STATE)
+	1                IYK(1:nobs,:),theta,S,STATE)
 	  CALL INNOV(nobs,d,ny,nz,nx,nu,ns,nt,S,
-	1             yk(1:nobs,:),IYK(1:nobs,:),theta,pdll,INN)
+	1             yk(1:nobs,:),IYK(1:nobs,:),theta,INN)
 	  IF (nf.GT.0) THEN
 	   CALL FORECAST(yk(nobs+1:nobs+nf,ny+1:ny+nz),nf,ny,nz,nx,nu,nv,
 	1                 ns,nstot,nt,np,theta,psi,INFOS,Z(nobs),
-     2                 STATE(nobs,:),pdll,FORE)
+     2                 STATE(nobs,:),FORE)
 	  ENDIF
 	  IF (INDMIS*nmis.GE.1) THEN
 	   J = 1
@@ -824,7 +609,7 @@ C MCMC RECORDING phase
 	    IF (IYK(I,ny+1).LT.ny) THEN
 		  K = ny-IYK(I,ny+1)
 	      CALL MISSING(yk(I,:),ny,nz,nx,nu,ns,nt,K,theta,
-	1                  S(I,1:6),STATE(I,:),pdll,ykmis(J:J+K-1))
+	1                  S(I,1:6),STATE(I,:),ykmis(J:J+K-1))
 	     J = J+K
 	    ENDIF
          ENDDO
@@ -834,22 +619,53 @@ C MCMC RECORDING phase
 	   IMIN(1) = CUMN(INDT(L1+IMIN(1)))
 	   IMAX    = MAXLOC(CUMN(INDT(L1+1:ntf)))
 	   IMAX(1) = CUMN(INDT(L1+IMAX(1)))
+#if defined(__CYGWIN32__) || defined(_WIN32)
 	   CALL system('cls')
+#endif
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       WRITE(MEXPRINT,11131) BURNIN
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+       WRITE(MEXPRINT,11132) ntf
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+       WRITE(MEXPRINT,11133) lastl,lasth
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
 	   WRITE(6,1113) BURNIN,ntf,lastl,lasth
+#endif
 	   IF ((HBL.EQ.1).OR.(nv.EQ.0)) THEN
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+          WRITE(MEXPRINT,11141) jjj
+          mpfout = mexPrintf(MEXPRINT//achar(13))
+          WRITE(MEXPRINT,11142) ntf
+          mpfout = mexPrintf(MEXPRINT//achar(13))
+          WRITE(MEXPRINT,11143) IMIN(1)/dfloat(jjj),IMAX(1)/dfloat(jjj)
+          mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
 	    WRITE(6,1114) jjj,ntf,IMIN(1)/dfloat(jjj),
      #           IMAX(1)/dfloat(jjj)
+#endif
          ELSEIF ((HBL.GT.1).AND.(nv.GT.0)) THEN
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+            WRITE(MEXPRINT,11151) jjj
+            mpfout = mexPrintf(MEXPRINT//achar(13))
+            WRITE(MEXPRINT,11152) ntf
+            mpfout = mexPrintf(MEXPRINT//achar(13))
+            WRITE(MEXPRINT,11153) IMIN(1)/dfloat(jjj),IMAX(1)/dfloat(jjj)
+            mpfout = mexPrintf(MEXPRINT//achar(13))
+            WRITE(MEXPRINT,11154) SUM(1.D0-ACCRATE(1:nobs)/DFLOAT(jjj))/DFLOAT(nobs)
+            mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
 	    WRITE(6,1115) jjj,ntf,IMIN(1)/dfloat(jjj),
      #           IMAX(1)/dfloat(jjj),
      #           SUM(1.D0-ACCRATE(1:nobs)/DFLOAT(jjj))/DFLOAT(nobs)
+#endif
 	   ENDIF
         ENDIF
         IF (jjj/thin*thin.EQ.jjj) THEN
 #ifdef __GFORTRAN__
-           WRITE(fmt, '(a,i4,a)') '(', nobs*ny, '(F20.10))'
+           WRITE(fmt, '(a,i4,a)') '(', nobs*ny, 'F20.10)'
            WRITE(12,fmt) (INN(1:nobs,I),I=1,ny)
-           WRITE(fmt, '(a,i4,a)') '(', nobs*nx, '(F20.10))'
+           WRITE(fmt, '(a,i4,a)') '(', nobs*nx, 'F20.10)'
            WRITE(10,fmt) (STATE(1:nobs,I),I=1,nx)
 #else
 	   WRITE(12,'(<nobs*ny>(F20.10))') (INN(1:nobs,I),I=1,ny)
@@ -860,7 +676,7 @@ C MCMC RECORDING phase
 	   ENDIF
 	   IF (nv.EQ.0) THEN
 #ifdef __GFORTRAN__
-          WRITE(fmt, '(a,i4,a)') '(', nt, '(F25.15))'
+          WRITE(fmt, '(a,i4,a)') '(', nt, 'F25.15)'
           WRITE(9,fmt) theta(1:nt)
 #else
 	    WRITE(9,'(<nt>(F25.15))') theta(1:nt)
@@ -871,9 +687,9 @@ C MCMC RECORDING phase
 	    gibZ(jjj/thin,1:nobs) = Z(1:nobs)
 	   ENDIF
 #ifdef __GFORTRAN__
-       WRITE(fmt, '(a,i4,a)') '(', nt+np(1), '(F25.15))'
+       WRITE(fmt, '(a,i4,a)') '(', nt+np(1), 'F25.15)'
        WRITE(9,fmt) theta(1:nt),psi(1:np(1))
-       WRITE(fmt, '(a,i4,a)') '(', nobs, '(I3))'
+       WRITE(fmt, '(a,i4,a)') '(', nobs, 'I3)'
 	   WRITE(11,fmt) Z(:)
 #else
 	   WRITE(9,'(<nt+np(1)>(F25.15))') theta(1:nt),psi(1:np(1))
@@ -883,7 +699,7 @@ C MCMC RECORDING phase
 	  IF (nf.GT.0) THEN
 	   J = min(nv,1)
 #ifdef __GFORTRAN__
-       WRITE(fmt, '(a,i4,a)') '(', nf*(nx+ny+J), '(F20.10))'
+       WRITE(fmt, '(a,i4,a)') '(', nf*(nx+ny+J), 'F20.10)'
        WRITE(13,fmt) (FORE(1:nf,I),I=1,nx+ny+J)
 #else
 	   WRITE(13,'(<nf*(nx+ny+J)>(F20.10))') (FORE(1:nf,I),I=1,nx+ny+J)
@@ -891,7 +707,7 @@ C MCMC RECORDING phase
 	  ENDIF
 #ifdef __GFORTRAN__
       IF (INDMIS*nmis.GE.1) THEN
-         WRITE(fmt, '(a,i4,a)') '(', nmis, '(F20.10))'
+         WRITE(fmt, '(a,i4,a)') '(', nmis, 'F20.10)'
          WRITE(14,fmt) ykmis(1:nmis)
       END IF
 #else
@@ -904,11 +720,11 @@ C MCMC RECORDING phase
         IF (nv.GT.0) THEN
 	   IF (HBL.EQ.1) THEN
 	    CALL GCK2(nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np(1),
-	1              yk(1:nobs,:),theta0,psi0,INFOS,pdll,Z,S)
+	1              yk(1:nobs,:),theta0,psi0,INFOS,Z,S)
 	   ELSE
 		CALL AMH2(hbl,nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np(1),
 	1              yk(1:nobs,:),theta0,psi0,
-	2              PTR,PM,INFOS,pdll,Z,S,ACCRATE)
+	2              PTR,PM,INFOS,Z,S,ACCRATE)
 		CALL RECPR(jjj+burnin,nstot,nobs,Z,ZW,PM,PTR)
 	   ENDIF
 	   CALL DRAWPSI(nobs,nv,np,INFOS,Z,psiprior,psi0,psi)
@@ -916,35 +732,66 @@ C MCMC RECORDING phase
 	  DO it = 1,nt
 	   IF (thetaprior(it,3).LT.thetaprior(it,4)) THEN
 	    CALL SLICE2(it,nobs,d,ny,nz,nx,nu,ns,nt,S,yk(1:nobs,:),
-	1                 theta0,thetaprior(it,:),pdftheta(it),pdll,
+	1                 theta0,thetaprior(it,:),pdftheta(it),
      2                NEVAL(it),theta(it))
           theta0(it) = theta(it)
          ENDIF
 	  END DO
 	  CUMN = CUMN+NEVAL
 	  CALL SIMSTATE2(nobs,d,ny,nz,nx,nu,ns,nt,yk(1:nobs,:),
-	1                 theta,S,pdll,STATE)
+	1                 theta,S,STATE)
 	  CALL INNOV2(nobs,d,ny,nz,nx,nu,ns,nt,S,
-	1              yk(1:nobs,:),theta,pdll,INN)
+	1              yk(1:nobs,:),theta,INN)
 	  IF (nf.GT.0) THEN
 	   CALL FORECAST(yk(nobs+1:nobs+nf,ny+1:ny+nz),nf,ny,nz,nx,nu,nv,
 	1                 ns,nstot,nt,np,theta,psi,INFOS,Z(nobs),
-     2                 STATE(nobs,:),pdll,FORE)
+     2                 STATE(nobs,:),FORE)
 	  ENDIF
 	  IF (jjj/IND*IND.EQ.jjj) THEN
 	   IMIN    = MINLOC(CUMN(INDT(L1+1:ntf)))
 	   IMIN(1) = CUMN(INDT(L1+IMIN(1)))
 	   IMAX    = MAXLOC(CUMN(INDT(L1+1:ntf)))
 	   IMAX(1) = CUMN(INDT(L1+IMAX(1)))
+#if defined(__CYGWIN32__) || defined(_WIN32)
 	   CALL system('cls')
+#endif
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       WRITE(MEXPRINT,11131) BURNIN
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+       WRITE(MEXPRINT,11132) ntf
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+       WRITE(MEXPRINT,11133) lastl,lasth
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
 	   WRITE(6,1113) BURNIN,ntf,lastl,lasth
+#endif
 	   IF ((HBL.EQ.1).OR.(nv.EQ.0)) THEN
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+          WRITE(MEXPRINT,11141) jjj
+          mpfout = mexPrintf(MEXPRINT//achar(13))
+          WRITE(MEXPRINT,11142) ntf
+          mpfout = mexPrintf(MEXPRINT//achar(13))
+          WRITE(MEXPRINT,11143) IMIN(1)/dfloat(jjj),IMAX(1)/dfloat(jjj)
+          mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
 	    WRITE(6,1114) jjj,ntf,IMIN(1)/dfloat(jjj),
      #           IMAX(1)/dfloat(jjj)
+#endif
          ELSEIF ((HBL.GT.1).AND.(nv.GT.0)) THEN
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+            WRITE(MEXPRINT,11151) jjj
+            mpfout = mexPrintf(MEXPRINT//achar(13))
+            WRITE(MEXPRINT,11152) ntf
+            mpfout = mexPrintf(MEXPRINT//achar(13))
+            WRITE(MEXPRINT,11153) IMIN(1)/dfloat(jjj),IMAX(1)/dfloat(jjj)
+            mpfout = mexPrintf(MEXPRINT//achar(13))
+            WRITE(MEXPRINT,11154) SUM(1.D0-ACCRATE(1:nobs)/DFLOAT(jjj))/DFLOAT(nobs)
+            mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
 	    WRITE(6,1115) jjj,ntf,IMIN(1)/dfloat(jjj),
      #           IMAX(1)/dfloat(jjj),
      #           SUM(1.D0-ACCRATE(1:nobs)/DFLOAT(jjj))/DFLOAT(nobs)
+#endif
 	   ENDIF
         ENDIF
         IF (jjj/thin*thin.EQ.jjj) THEN
@@ -952,18 +799,18 @@ C MCMC RECORDING phase
 	    gibtheta(jjj/thin,1:nt) = theta(1:nt)
 	   ENDIF
 #ifdef __GFORTRAN__
-       WRITE(fmt, '(a,i4,a)') '(', nobs*ny, '(F20.10))'
+       WRITE(fmt, '(a,i4,a)') '(', nobs*ny, 'F20.10)'
        WRITE(12,fmt) (INN(1:nobs,I),I=1,ny)
-       WRITE(fmt, '(a,i4,a)') '(', nobs*nx, '(F20.10))'
-	   WRITE(10,fmt) (STATE(1:nobs,I),I=1,nx)
+       WRITE(fmt, '(a,i4,a)') '(', nobs*nx, 'F20.10)'
+       WRITE(10,fmt) (STATE(1:nobs,I),I=1,nx)
 #else
 	   WRITE(12,'(<nobs*ny>(F20.10))') (INN(1:nobs,I),I=1,ny)
 	   WRITE(10,'(<nobs*nx>(F20.10))') (STATE(1:nobs,I),I=1,nx)
 #endif
 	   IF (nv.EQ.0) THEN
 #ifdef __GFORTRAN__
-       WRITE(fmt, '(a,i4,a)') '(', nt, '(F25.15))'
-       WRITE(9,fmt) theta(1:nt)
+          WRITE(fmt, '(a,i4,a)') '(', nt, 'F25.15)'
+          WRITE(9,fmt) theta(1:nt)
 #else
 	    WRITE(9,'(<nt>(F25.15))') theta(1:nt)
 #endif
@@ -973,10 +820,10 @@ C MCMC RECORDING phase
 	     gibZ(jjj/thin,1:nobs) = Z(1:nobs)
 	    ENDIF
 #ifdef __GFORTRAN__
-        WRITE(fmt, '(a,i4,a)') '(', nt+np(1), '(F25.15))'
-	    WRITE(9,fmt) theta(1:nt),psi(1:np(1))
-        WRITE(fmt, '(a,i4,a)') '(', nobs, '(F25.15))'
-	    WRITE(11,fmt) Z(:)
+        WRITE(fmt, '(a,i4,a)') '(', nt+np(1), 'F25.15)'
+        WRITE(9,fmt) theta(1:nt),psi(1:np(1))
+        WRITE(fmt, '(a,i4,a)') '(', nobs, 'I3)'
+        WRITE(11,fmt) Z(:)
 #else
 	    WRITE(9,'(<nt+np(1)>(F25.15))') theta(1:nt),psi(1:np(1))
 	    WRITE(11,'(<nobs>(I3))') Z(:)
@@ -985,7 +832,7 @@ C MCMC RECORDING phase
 	   IF (nf.GT.0) THEN
 	   J = min(nv,1)
 #ifdef __GFORTRAN__
-       WRITE(fmt, '(a,i4,a)') '(', nf*(nx+ny+J), '(F20.10))'
+       WRITE(fmt, '(a,i4,a)') '(', nf*(nx+ny+J), 'F20.10)'
        WRITE(13,fmt) (FORE(1:nf,I),I=1,nx+ny+J)
 #else
 	   WRITE(13,'(<nf*(nx+ny+J)>(F20.10))') (FORE(1:nf,I),I=1,nx+ny+J)
@@ -1010,35 +857,65 @@ C MCMC RECORDING phase
 
 C MARGINAL LIKELIHOOD
 	IF ((MargLik.EQ.'Y').OR.(MargLik.EQ.'y')) THEN
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       WRITE(MEXPRINT,*) ' '
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+       WRITE(MEXPRINT,*) 'Computing the marginal likelihood. Please wait ...'
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
 	 WRITE(*,*) ' '
        WRITE(*,*) 'Computing the marginal likelihood. Please wait ...'
+#endif
 	 IF (nmis.GT.0) THEN
 	  CALL HARMONIC(GGG,nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np,
 	1                INFOS,yk(1:nobs,:),IYK(1:nobs,:),gibtheta,gibZ,
-     2                thetaprior,psiprior,pdftheta,pdll,MLHM)
+     2                thetaprior,psiprior,pdftheta,MLHM)
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       WRITE(MEXPRINT,*) 'Modified harmonic mean: done!'
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
 	  WRITE(*,*) 'Modified harmonic mean: done!'
+#endif
 	  CALL MENGWONG(GGG,nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np,
 	1                INFOS,yk(1:nobs,:),IYK(1:nobs,:),gibtheta,gibZ,
-     2                thetaprior,psiprior,pdftheta,pdll,MLHM(5,1),MLMW)
+     2                thetaprior,psiprior,pdftheta,MLHM(5,1),MLMW)
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       WRITE(MEXPRINT,*) 'Bridge sampling: done!'
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+       WRITE(MEXPRINT,*) ' '
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
         WRITE(*,*) 'Bridge sampling: done!'
 	  WRITE(*,*) ' '
+#endif
 	 ELSE
 	  CALL HARMONIC2(GGG,nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np,
 	1                 INFOS,yk(1:nobs,:),gibtheta,gibZ,thetaprior,
-     2                 psiprior,pdftheta,pdll,MLHM)
+     2                 psiprior,pdftheta,MLHM)
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       WRITE(MEXPRINT,*) 'Modified harmonic mean: done!'
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
 	  WRITE(*,*) 'Modified harmonic mean: done!'
+#endif
 	  CALL MENGWONG2(GGG,nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np,
 	1                 INFOS,yk(1:nobs,:),gibtheta,gibZ,thetaprior,
-     2                 psiprior,pdftheta,pdll,MLHM(5,1),MLMW)
+     2                 psiprior,pdftheta,MLHM(5,1),MLMW)
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       WRITE(MEXPRINT,*) 'Bridge sampling: done!'
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+       WRITE(MEXPRINT,*) ' '
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
         WRITE(*,*) 'Bridge sampling: done!'
 	  WRITE(*,*) ' '
+#endif
 	 ENDIF
 	 WRITE(15,*) 'Modified Harmonic mean (ML and Var)'
 #ifdef __GFORTRAN__
-      WRITE(fmt, '(a,i4,a)') '(', 2, '(F20.10))'
-      WRITE(15,fmt) (MLHM(I,:),I=1,11)
+      WRITE(15,'(2F20.10)') (MLHM(I,:),I=1,11)
       WRITE(15,*) 'Bridge Sampling'
-      WRITE(15,fmt) (MLMW(I,:),I=1,2)
+      WRITE(15,'(2F20.10)') (MLMW(I,:),I=1,2)
 #else
 	 WRITE(15,'(<2>(F20.10))') (MLHM(I,:),I=1,11)
 	 WRITE(15,*) 'Bridge Sampling'
@@ -1054,15 +931,6 @@ C MARGINAL LIKELIHOOD
 	 DEALLOCATE(psi0,psi,psiprior,ZW)
 	ENDIF
 
-#if defined(__CYGWIN32__) || defined(_WIN32)
-      STATUS = freelibrary(pdll) !libero la DLL dalla memoria alla fine del programma
-#else
-      STATUS = DLClose(pdll)
-      IF(STATUS/=0) THEN
-         WRITE(*,*) ' Error in dlclose: ', C_F_STRING(DLError())
-         STOP
-      END IF
-#endif
 	IF (TRIM(PATH).EQ.'') THEN
 	 STATUS = getcwd(PATH) ! get current directory
       ENDIF
@@ -1073,16 +941,48 @@ C MARGINAL LIKELIHOOD
       IT2(4:7) = DATE_ITIME(5:8)
 	IT=(IT2(4)-IT1(4))*3600+(IT2(5)-IT1(5))*60+(IT2(6)-IT1(6))
       IF ((check.EQ.'Y').OR.(check.EQ.'y')) THEN
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+       WRITE(MEXPRINT,11171)
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+       WRITE(MEXPRINT,11172) TRIM(PATH)
+       mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
         WRITE(6,1117) TRIM(PATH)
+#endif
       ELSE
           IF ((datasim.EQ.'Y').OR.(datasim.EQ.'y')) THEN
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+             WRITE(MEXPRINT,11181)
+             mpfout = mexPrintf(MEXPRINT//achar(13))
+             WRITE(MEXPRINT,11182) TRIM(PATH)
+             mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
             WRITE(6,1118) TRIM(PATH)
+#endif
           ELSE
             IF ((estimation.EQ.'ML').OR.(estimation.EQ.'ml').OR.
      &          (estimation.EQ.'Ml').OR.(estimation.EQ.'mL')) THEN
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+               WRITE(MEXPRINT,11191)
+               mpfout = mexPrintf(MEXPRINT//achar(13))
+               WRITE(MEXPRINT,11192) IT
+               mpfout = mexPrintf(MEXPRINT//achar(13))
+               WRITE(MEXPRINT,11193) TRIM(PATH)
+               mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
               WRITE(6,1119) IT,TRIM(PATH)
+#endif
             ELSE
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+               WRITE(MEXPRINT,11161)
+               mpfout = mexPrintf(MEXPRINT//achar(13))
+               WRITE(MEXPRINT,11162) IT
+               mpfout = mexPrintf(MEXPRINT//achar(13))
+               WRITE(MEXPRINT,11163) TRIM(PATH)
+               mpfout = mexPrintf(MEXPRINT//achar(13))
+#else
               WRITE(6,1116) IT,TRIM(PATH)
+#endif
             ENDIF
           ENDIF
       ENDIF
@@ -1092,6 +992,28 @@ C MARGINAL LIKELIHOOD
 1111  FORMAT((<4>(F25.12)), '  ',A2)
 1112  FORMAT(I10,(<np(3)>(F25.12)), '  ',I2)
 #endif
+#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
+11131 FORMAT(' Burn-in draws = ',I8)
+11132 FORMAT(' Parameters sampled by SLICE ',I5)
+11133 FORMAT(' SLICE likelihood eval. Min/Max = ',F6.2, ' / ',F6.2)
+11141 FORMAT(' Recording draws = ',I8)
+11142 FORMAT(' Parameters sampled by SLICE ',I5)
+11143 FORMAT(' SLICE likelihood eval. Min/Max = ',F6.2, ' / ',F6.2)
+11151 FORMAT(' Recording draws = ',I8)
+11152 FORMAT(' Parameters sampled by SLICE ',I5)
+11153 FORMAT(' SLICE likelihood eval. Min/Max = ',F6.2, ' / ',F6.2)
+11154 FORMAT(' Adaptive MH accettance rate = ',F6.2)
+11161 FORMAT(' MCMC completed')
+11162 FORMAT(' CPU-time (sec)=', I10)
+11163 FORMAT(' Output printed in ',A)
+11171 FORMAT(' Check completed')
+11172 FORMAT(' Output printed in ',A)
+11181 FORMAT(' Data simulation completed')
+11182 FORMAT(' Output printed in 'A)
+11191 FORMAT(' Maximum Likelihood completed')
+11192 FORMAT(' CPU-time (sec)=', I10)
+11193 FORMAT(' Output printed in ',A)
+#else
 1113  FORMAT(/,' Burn-in draws = ',I8,
      #       /,' Parameters sampled by SLICE ',I5,
      #       /,' SLICE likelihood eval. Min/Max = ',F6.2, ' / ',F6.2)
@@ -1112,8 +1034,11 @@ C MARGINAL LIKELIHOOD
 1119  FORMAT(/,' Maximum Likelihood completed',
      #       /,' CPU-time (sec)=', I10,
      #       /,' Output printed in ',A)
+#endif
 #ifdef __INTEL_COMPILER
       PAUSE
 #endif
+#if !(defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE))
 	STOP
+#endif
       END
diff --git a/mengwong.for b/mengwong.for
index a7ac08a03f511e98890ebf99ac13df15d6ae61f3..e68c1e09c70ff28ec7ad49eb43e33a311cbcbb70 100644
--- a/mengwong.for
+++ b/mengwong.for
@@ -34,7 +34,7 @@ C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C -------------------------------------------------------------------
 	SUBROUTINE MENGWONG(G,nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np,
 	1                    INFOS,yk,IYK,gibpar,gibZ,thetaprior,psiprior,
-     2                    tipo,pdll,MLSTART,MLMW)
+     2                    tipo,MLSTART,MLMW)
 #ifdef __GFORTRAN__
       USE gfortran
 #endif
@@ -44,7 +44,6 @@ C INPUT
 	DOUBLE PRECISION yk(nobs,ny+nz),gibpar(G,nt+np(1)),
 	1 thetaprior(nt,4),psiprior(np(2),np(3)),MLSTART
 	CHARACTER*2 tipo(nt)
-	POINTER (pdll,fittizia)  !  ASSOCIATE  pointer P alla DLL ad una varibile fittizia
 
 C OUTPUT
 	DOUBLE PRECISION MLMW(2,2)
@@ -332,7 +331,7 @@ C PRIOR for PSI and Q(PSI)~Dirichlet(a1,a2,...,aN)
 
 	 Fpar = PTHETA(NPOS(1),nobs,d,ny,nz,nx,nu,ns,nt,IS,yk,IYK,
 	1               par(1:nt),thetaprior(NPOS(1),:),
-     2               tipo(NPOS(1)),pdll)
+     2               tipo(NPOS(1)))
        Fpar = Fpar + SUM(Ppar(2:NPARTH+K)) ! log f(y|par,S)f(par,S)
 
 	 VQN(IG,1) = DEXP(QPSI)*QS*C*DEXP(-.5D0*MUC)/TRC
@@ -393,7 +392,7 @@ C PRIOR for PSI and Q(PSI)~Dirichlet(a1,a2,...,aN)
 
 	 Fpar = PTHETA(NPOS(1),nobs,d,ny,nz,nx,nu,ns,nt,IS,yk,IYK,
 	1               gibpar(IG,1:nt),thetaprior(NPOS(1),:),
-     2               tipo(NPOS(1)),pdll)
+     2               tipo(NPOS(1)))
        Fpar = Fpar + SUM(Ppar(2:NPARTH+K)) ! log f(y|par,S)f(par,S)
 
 	 VHD(IG,1) = Fpar + DLOG(PS)
diff --git a/mengwong2.for b/mengwong2.for
index 2778e1dd357bbc07505a597e5c6b6811b65e9ac9..3c333707c47927474736e43f502e981017f7ee1a 100644
--- a/mengwong2.for
+++ b/mengwong2.for
@@ -34,7 +34,7 @@ C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C -------------------------------------------------------------------
 	SUBROUTINE MENGWONG2(G,nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np,
 	1                     INFOS,yk,gibpar,gibZ,thetaprior,psiprior,
-     2                     tipo,pdll,MLSTART,MLMW)
+     2                     tipo,MLSTART,MLMW)
 #ifdef __GFORTRAN__
       USE gfortran
 #endif
@@ -44,7 +44,6 @@ C INPUT
 	DOUBLE PRECISION yk(nobs,ny+nz),gibpar(G,nt+np(1)),
 	1 thetaprior(nt,4),psiprior(np(2),np(3)),MLSTART
 	CHARACTER*2 tipo(nt)
-	POINTER (pdll,fittizia)
 
 C OUTPUT
 	DOUBLE PRECISION MLMW(2,2)
@@ -331,7 +330,7 @@ C PRIOR for PSI and Q(PSI)~Dirichlet(a1,a2,...,aN)
 
 	 Fpar = PTHETA2(NPOS(1),nobs,d,ny,nz,nx,nu,ns,nt,IS,yk,
 	1                par(1:nt),thetaprior(NPOS(1),:),
-     2                tipo(NPOS(1)),pdll)
+     2                tipo(NPOS(1)))
        Fpar = Fpar + SUM(Ppar(2:NPARTH+K)) ! log f(y|par,S)f(par,S)
 
 	 VQN(IG,1) = DEXP(QPSI)*QS*C*DEXP(-.5D0*MUC)/TRC
@@ -392,7 +391,7 @@ C PRIOR for PSI and Q(PSI)~Dirichlet(a1,a2,...,aN)
 
 	 Fpar = PTHETA2(NPOS(1),nobs,d,ny,nz,nx,nu,ns,nt,IS,yk,
 	1                gibpar(IG,1:nt),thetaprior(NPOS(1),:),
-     2                tipo(NPOS(1)),pdll)
+     2                tipo(NPOS(1)))
        Fpar = Fpar + SUM(Ppar(2:NPARTH+K)) ! log f(y|par,S)f(par,S)
 
 	 VHD(IG,1) = Fpar + DLOG(PS)
diff --git a/missing.for b/missing.for
index ba8657235b00b8947baabc45a8064736cef50347..e84805bef76277ab61cc2ce67abcb3a40a458d4a 100644
--- a/missing.for
+++ b/missing.for
@@ -20,7 +20,7 @@ C
 C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C --------------------------------------------------------------------
-	SUBROUTINE MISSING(yk,ny,nz,nx,nu,ns,nt,nmis,theta,S,STATE,pdll,ykmis)
+	SUBROUTINE MISSING(yk,ny,nz,nx,nu,ns,nt,nmis,theta,S,STATE,ykmis)
 #if defined(__CYGWIN32__) || defined(_WIN32)
 #ifdef __INTEL_COMPILER
       USE dfwin
@@ -30,23 +30,6 @@ C --------------------------------------------------------------------
       USE ISO_C_UTILITIES
       USE DLFCN
 #endif
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-	CHARACTER*1 fittizia
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  POINTER (pdll,fittizia)
-	  POINTER (pdesign,DESIGN)
-#else
-	  TYPE(C_PTR) :: pdll
-      TYPE(C_FUNPTR) :: pdesign=C_NULL_FUNPTR
-	  PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-#endif
 
 C INPUT
 	INTEGER ny,nz,nx,nu,nt,ns(6),nmis
@@ -63,16 +46,9 @@ C LOCALS
 
 	ALLOCATE(R(nx,nu,ns(6)),c(ny,max(nz,1),ns(1)),H(ny,nx,ns(2)),
 	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)))
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  pdesign = getprocaddress(pdll, "design_"C)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 	  CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-	  pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-      IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-         WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-	  END IF
-	  CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-	  CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
 
 C NIYK = not(IYK)
diff --git a/opg.for b/opg.for
index 205e87bb70aba7664f68927fab473b5286ed39d5..60c3ae1ac298622e45ab468194fd80c95825ac3b 100644
--- a/opg.for
+++ b/opg.for
@@ -27,7 +27,7 @@ C
 C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C ------------------------------------------------------------
-	SUBROUTINE OPG(nobs,d,ny,nz,nx,nu,nt,ns,pdll,yk,IYK,S,
+	SUBROUTINE OPG(nobs,d,ny,nz,nx,nu,nt,ns,yk,IYK,S,
 	1 theta,thetaprior,HESS,SE,XS,AKMSE,INN,IFAIL)
 #if defined(__CYGWIN32__) || defined(_WIN32)
 #ifdef __INTEL_COMPILER
@@ -38,28 +38,8 @@ C ------------------------------------------------------------
       USE ISO_C_UTILITIES
       USE DLFCN
 #endif
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-	CHARACTER*1 fittizia
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  POINTER (pdll,fittizia)	! ASSOCIATE  pointer pdll alla DLL ad una varibile fittizia
-	  POINTER (pdesign,DESIGN)	! IMPORTANT associo il puntatore pdesign alla Interface definita
-#else
-	  TYPE(C_PTR) :: pdll
-      TYPE(C_FUNPTR) :: pdesign=C_NULL_FUNPTR
-	  PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-#endif
 
 ! Input
-#ifdef __INTEL_COMPILER
-      INTEGER pdll
-#endif
 	INTEGER nobs,d(2),ny,nz,nx,nu,nt,ns(6),IYK(nobs,ny+1),
 	1 S(nobs,6)
       DOUBLE PRECISION yk(nobs,ny+nz),theta(nt),thetaprior(nt,4),
@@ -105,16 +85,9 @@ C Using Hessian from E04UCF
      1 XT(0:nobs,nx),PT(0:nobs,nx,nx),Xdd(max(d(1),1),nx),
      1 Pdd(max(d(1),1),nx,nx))
 
-#if defined(__CYGWIN32__) || defined(_WIN32)
-      pdesign = getprocaddress(pdll, "design_"C)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 	  CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-	  pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-      IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-         WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-	  END IF
-	  CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-	  CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
       CALL IKF(d,ny,nz,nx,nu,ns,S(1:max(d(1),1),1:6),
 	1         yk(1:max(d(1),1),1:ny+nz),IYK(1:max(d(1),1),1:ny+1),
@@ -150,10 +123,9 @@ C Using Hessian from E04UCF
 	DO 1000 I=1,NFREE
 	 THETAV(I) = THETAV(I) + P(I)
        theta(IFREE(I)) = THETAV(I)
-#if defined(__CYGWIN32__) || defined(_WIN32)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 	   CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-	   CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
 	 CALL IKF(d,ny,nz,nx,nu,ns,S(1:max(d(1),1),1:6),
 	1          yk(1:max(d(1),1),1:ny+nz),IYK(1:max(d(1),1),1:ny+1),
diff --git a/opgh.for b/opgh.for
index d99573846232f7fefa1f3d73bb98a3e4e91ff1eb..fa3ca697b2c79da22d4f3a07f3fe87d07e256225 100644
--- a/opgh.for
+++ b/opgh.for
@@ -28,7 +28,7 @@ C
 C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C ------------------------------------------------------------
-      SUBROUTINE OPGH(nobs,ny,nz,nx,nu,nt,nv,ns,nstot,np,pdll,yk,IYK,
+      SUBROUTINE OPGH(nobs,ny,nz,nx,nu,nt,nv,ns,nstot,np,yk,IYK,
 	1                INFOS,theta,psi,thetaprior,HESS,thetase,psise,
      1                SSMOOTH,INN,IFAIL)
 #if defined(__CYGWIN32__) || defined(_WIN32)
@@ -40,28 +40,8 @@ C ------------------------------------------------------------
       USE ISO_C_UTILITIES
       USE DLFCN
 #endif
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-	CHARACTER*1 fittizia
-#if defined(__CYGWIN32__) || defined(_WIN32)
-      POINTER (pdll,fittizia)   ! ASSOCIATE  pointer pdll alla DLL ad una varibile fittizia
-      POINTER (pdesign,DESIGN)  ! IMPORTANT associo il puntatore pdesign alla Interface definita
-#else
-	  TYPE(C_PTR) :: pdll
-      TYPE(C_FUNPTR) :: pdesign=C_NULL_FUNPTR
-      PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-#endif
 
 ! Input
-#ifdef __INTEL_COMPILER
-      INTEGER pdll
-#endif
 	INTEGER nobs,ny,nz,nx,nu,nt,nv,ns(6),nstot,np,IYK(nobs,ny+1),
      1 INFOS(9,6)
       DOUBLE PRECISION yk(nobs,ny+nz),theta(nt),psi(np),
@@ -108,16 +88,9 @@ C ------------------------------------------------------------
       IFAIL = 0
       CALL SYMINV(LTR,NFREE,LTR,W,J,IFAIL,RMAX)
 
-#if defined(__CYGWIN32__) || defined(_WIN32)
-      pdesign = getprocaddress(pdll, "design_"C)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
       CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-	  pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-      IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-         WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-	  END IF
-      CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-	  CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
       CALL HF(nobs,nx,nstot,nz,nu,ns,nv,np,psi,1,yk,IYK,INFOS,
      1        c,a,F,R,SSMOOTH,INN,DLL)
@@ -162,10 +135,9 @@ C -----------
         ELSE
          psi(I-NFT) = PAR(I)
         ENDIF
-#if defined(__CYGWIN32__) || defined(_WIN32)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 	  CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-      CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
         CALL HF(nobs,nx,nstot,nz,nu,ns,nv,np,psi,0,yk,IYK,INFOS,
      1        c,a,F,R,SSMOOTH,INN,DLLM)
diff --git a/opgkim.for b/opgkim.for
index 267662eaa50657aa3aeb2b12e846adb810e5f650..014afe7ad155a418781af85ca8a2ed256829c5d0 100644
--- a/opgkim.for
+++ b/opgkim.for
@@ -27,7 +27,7 @@ C
 C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C ------------------------------------------------------------
-      SUBROUTINE OPGKIM(nobs,d,ny,nz,nx,nu,nt,nv,ns,nstot,np,pdll,yk,
+      SUBROUTINE OPGKIM(nobs,d,ny,nz,nx,nu,nt,nv,ns,nstot,np,yk,
      1                  IYK,INFOS,theta,psi,thetaprior,HESS,thetase,
      1                  psise,XS,XSSE,SSMOOTH,INN,IFAIL)
 #if defined(__CYGWIN32__) || defined(_WIN32)
@@ -39,28 +39,8 @@ C ------------------------------------------------------------
       USE ISO_C_UTILITIES
       USE DLFCN
 #endif
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-	CHARACTER*1 fittizia
-#if defined(__CYGWIN32__) || defined(_WIN32)
-      POINTER (pdll,fittizia)   ! ASSOCIATE  pointer pdll alla DLL ad una varibile fittizia
-      POINTER (pdesign,DESIGN)  ! IMPORTANT associo il puntatore pdesign alla Interface definita
-#else
-	  TYPE(C_PTR) :: pdll
-      TYPE(C_FUNPTR) :: pdesign=C_NULL_FUNPTR
-      PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-#endif
 
 ! Input
-#ifdef __INTEL_COMPILER
-      INTEGER pdll
-#endif
 	INTEGER nobs,d(2),ny,nz,nx,nu,nt,nv,ns(6),nstot,np,IYK(nobs,ny+1),
      1 INFOS(9,6)
       DOUBLE PRECISION yk(nobs,ny+nz),theta(nt),psi(np),
@@ -107,16 +87,9 @@ C ------------------------------------------------------------
       IFAIL = 0
       CALL SYMINV(LTR,NFREE,LTR,W,J,IFAIL,RMAX)
 
-#if defined(__CYGWIN32__) || defined(_WIN32)
-      pdesign = getprocaddress(pdll, "design_"C)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
       CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-      pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-      IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-         WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-	  END IF
-      CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-	  CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
       CALL KIM(nobs,d,ny,nz,nx,nu,ns,nstot,nv,np,INFOS,yk,IYK,
      1         c,H,G,a,F,R,psi,1,XS,XSSE,SSMOOTH,INN,DLL)
@@ -162,10 +135,9 @@ C -----------
         ELSE
          psi(I-NFT) = PAR(I)
         ENDIF
-#if defined(__CYGWIN32__) || defined(_WIN32)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
         CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-        CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
         CALL KIM(nobs,d,ny,nz,nx,nu,ns,nstot,nv,np,INFOS,yk,IYK,
      1           c,H,G,a,F,R,psi,0,XS,XSSE,SSMOOTH,INN,DLLM)
diff --git a/ptheta.for b/ptheta.for
index 51ea229d53494256bef8084be2b7d13f6ed127b8..537b16809afd126aa37887a24a910b4f1f6a32b8 100644
--- a/ptheta.for
+++ b/ptheta.for
@@ -21,7 +21,7 @@ C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C ----------------------------------------------------------------
 	DOUBLE PRECISION FUNCTION PTHETA(it,nobs,d,ny,nz,nx,nu,ns,nt,
-	1 S,yk,IYK,theta,thetaprior,tipo,pdll)
+	1 S,yk,IYK,theta,thetaprior,tipo)
 #if defined(__CYGWIN32__) || defined(_WIN32)
 #ifdef __INTEL_COMPILER
       USE dfwin
@@ -31,23 +31,6 @@ C ----------------------------------------------------------------
       USE ISO_C_UTILITIES
       USE DLFCN
 #endif
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-	CHARACTER*1 fittizia
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  POINTER (pdll,fittizia)	! ASSOCIATE  pointer pdll alla DLL ad una varibile fittizia
-	  POINTER (pdesign,DESIGN)	! IMPORTANT associo il puntatore pdesign alla Interface definita
-#else
-	  TYPE(C_PTR) :: pdll
-      TYPE(C_FUNPTR) :: pdesign=C_NULL_FUNPTR
-	  PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-#endif
 
 C INPUT
 	INTEGER it,nobs,d(2),ny,nz,nx,nu,ns(6),nt,S(nobs,6),
@@ -66,16 +49,9 @@ C LOCALS
      3 Pdd(max(d(1),1),nx,nx))
 
 C computes the log-posterior
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  pdesign = getprocaddress(pdll, "design_"C)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 	  CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-	  pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-      IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-         WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-	  END IF
-	  CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-	  CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
 	CALL IKF(d,ny,nz,nx,nu,ns,S(1:max(d(1),1),1:6),
 	1         yk(1:max(d(1),1),1:ny+nz),IYK(1:max(d(1),1),1:ny+1),
diff --git a/ptheta2.for b/ptheta2.for
index 8b36c62cb5e0cb13427cc25eda991b1f461d7985..c1e44ed5f5ffb52daa8c5b3c5f65ae4d1e64278a 100644
--- a/ptheta2.for
+++ b/ptheta2.for
@@ -21,7 +21,7 @@ C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C ----------------------------------------------------------------
 	DOUBLE PRECISION FUNCTION PTHETA2(it,nobs,d,ny,nz,nx,nu,ns,nt,
-	1 S,yk,theta,thetaprior,tipo,pdll)
+	1 S,yk,theta,thetaprior,tipo)
 #if defined(__CYGWIN32__) || defined(_WIN32)
 #ifdef __INTEL_COMPILER
       USE dfwin
@@ -31,23 +31,6 @@ C ----------------------------------------------------------------
       USE ISO_C_UTILITIES
       USE DLFCN
 #endif
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-	CHARACTER*1 fittizia
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  POINTER (pdll,fittizia)	! ASSOCIATE  pointer pdll alla DLL ad una varibile fittizia
-	  POINTER (pdesign,DESIGN)	! IMPORTANT associo il puntatore pdesign alla Interface definita
-#else
-	  TYPE(C_PTR) :: pdll
-      TYPE(C_FUNPTR) :: pdesign=C_NULL_FUNPTR
-	  PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-#endif
 
 C INPUT
 	INTEGER it,nobs,d(2),ny,nz,nx,nu,ns(6),nt,S(nobs,6)
@@ -66,16 +49,9 @@ C LOCALS
 
 
 C computes the log-posterior
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  pdesign = getprocaddress(pdll, "design_"C)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 	  CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-	  pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-      IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-         WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-	  END IF
-	  CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-	  CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
 	CALL IKF2(d,ny,nz,nx,nu,ns,S(1:max(d(1),1),1:6),
 	1          yk(1:max(d(1),1),1:ny+nz),
diff --git a/setfilem.for b/setfilem.for
new file mode 100644
index 0000000000000000000000000000000000000000..7c9f20b7bb36f24819fb967243ca1b0bde510860
--- /dev/null
+++ b/setfilem.for
@@ -0,0 +1,32 @@
+C --------------------------------------------------------------------
+C This file crates matlabdll.dll whih is ment to read the .m file
+C
+C Copyright (C) 2010-2014 European Commission
+C
+C This file is part of Program DMM
+C
+C DMM is free software developed at the Joint Research Centre of the
+C European Commission: you can redistribute it and/or modify it under
+C the terms of the GNU General Public License as published by
+C the Free Software Foundation, either version 3 of the License, or
+C (at your option) any later version.
+C
+C DMM is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
+C ---------------------------------------------------------------------
+C -----------------------------------------
+C To make dynamic the name of the .m file
+C -----------------------------------------
+	SUBROUTINE SETFILEM(string1,string2)
+!DEC$ ATTRIBUTES DLLEXPORT, ALIAS:'setfilem_' :: SETFILEM
+	CHARACTER*200 string1,string2,mfile,pathmfile
+      COMMON /M/ mfile, pathmfile
+      mfile     = string1
+      pathmfile = string2
+	RETURN
+      END
diff --git a/simdata.for b/simdata.for
index 9e1cbbd8a418557f01271d32efb3e00a7f895a9d..4c88df50b18246b5867ae0d84f6c32d0a7bc847a 100644
--- a/simdata.for
+++ b/simdata.for
@@ -48,7 +48,7 @@ C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C --------------------------------------------------------------------
 	SUBROUTINE SIMDATA(nobs,d,ny,nz,nx,nu,ns,nstot,nt,nv,np,INFOS,
-	1                   pdll,theta,psi,Z,STATE,yk)
+	1                   theta,psi,Z,STATE,yk)
 #if defined(__CYGWIN32__) || defined(_WIN32)
 #ifdef __INTEL_COMPILER
       USE dfwin
@@ -58,22 +58,6 @@ C --------------------------------------------------------------------
       USE ISO_C_UTILITIES
       USE DLFCN
 #endif
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	POINTER (pdll,fittizia)  ! ASSOCIATE  pointer pdll alla DLL ad una varibile fittizia
-	POINTER (pdesign,DESIGN) ! IMPORTANT associo il puntatore pdesign alla Interface definita
-#else
-	  TYPE(C_PTR) :: pdll
-      TYPE(C_FUNPTR) :: pdesign=C_NULL_FUNPTR
-	  PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-#endif
 
 C INPUT
 	INTEGER nobs,d(2),ny,nz,nx,nu,ns(6),nstot,nt,nv,np(3),INFOS(9,6)
@@ -108,16 +92,9 @@ C EXTERNAL SUBROUTINES
      3 Xdd(MAX(d(1),1),nx),Pdd(MAX(d(1),1),nx,nx),
      3 WORK((nx+2)*(nx+1)/2),FP(nx,nx),WORK1(64*nx),UP(nu) )
 
-#if defined(__CYGWIN32__) || defined(_WIN32)
-      pdesign = getprocaddress(pdll, "design_"C)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 	  CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-	  pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-      IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-         WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-	  END IF
-	  CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-	  CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
 
 C DRAW Z ~ Pr(S1 x ... x Snv|psi)
diff --git a/simstate.for b/simstate.for
index 6ee4c354c7d4cb778150a1060989cce8edf1d006..f46a772d456c9f35d74a25ff16af85b8e51df236 100644
--- a/simstate.for
+++ b/simstate.for
@@ -42,7 +42,7 @@ C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C --------------------------------------------------------------------
 	SUBROUTINE SIMSTATE(nobs,d,ny,nz,nx,nu,ns,nt,yk,IYK,
-	1                    theta,S,pdll,STATE)
+	1                    theta,S,STATE)
 #if defined(__CYGWIN32__) || defined(_WIN32)
 #ifdef __INTEL_COMPILER
       USE dfwin
@@ -52,22 +52,6 @@ C --------------------------------------------------------------------
       USE ISO_C_UTILITIES
       USE DLFCN
 #endif
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  POINTER (pdll,fittizia)	! ASSOCIATE  pointer pdll alla DLL ad una varibile fittizia
-	  POINTER (pdesign,DESIGN)	! IMPORTANT associo il puntatore pdesign alla Interface definita
-#else
-	  TYPE(C_PTR) :: pdll
-      TYPE(C_FUNPTR) :: pdesign=C_NULL_FUNPTR
-	  PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-#endif
 
 C INPUT
 	INTEGER nobs,d(2),ny,nz,nx,nu,nt,ns(6),S(nobs,6),
@@ -93,17 +77,9 @@ C EXTERNAL FUNCTIONS
      3 Xdd(MAX(d(1),1),nx),Pdd(MAX(d(1),1),nx,nx),
      3 WORK((nx+2)*(nx+1)/2),FP(nx,nx),WORK1(64*nx),UP(nu) )
 
-C	pdesign = getprocaddress(pdll, "DESIGN"C)
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  pdesign = getprocaddress(pdll, "design_"C)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
       CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-	  pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-      IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-         WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-	  END IF
-	  CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-	  CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
 
 	ykP(:,:) = 0.D0
diff --git a/simstate2.for b/simstate2.for
index 8b580d4dc3c84726867caf30876d467e4d0587ba..1c4e8e9f65a8446bf11d004a599aabb81f197339 100644
--- a/simstate2.for
+++ b/simstate2.for
@@ -42,7 +42,7 @@ C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C --------------------------------------------------------------------
 	SUBROUTINE SIMSTATE2(nobs,d,ny,nz,nx,nu,ns,nt,yk,
-	1                     theta,S,pdll,STATE)
+	1                     theta,S,STATE)
 #if defined(__CYGWIN32__) || defined(_WIN32)
 #ifdef __INTEL_COMPILER
       USE dfwin
@@ -52,22 +52,6 @@ C --------------------------------------------------------------------
       USE ISO_C_UTILITIES
       USE DLFCN
 #endif
-	INTERFACE
-	 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
-	 INTEGER ny,nz,nx,nu,ns(6),nt
-	 DOUBLE PRECISION theta(nt)
-	 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)),
-	1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))
-	 END SUBROUTINE
-	END INTERFACE
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  POINTER (pdll,fittizia)	! ASSOCIATE  pointer pdll alla DLL ad una varibile fittizia
-	  POINTER (pdesign,DESIGN)	! IMPORTANT associo il puntatore pdesign alla Interface definita
-#else
-	  TYPE(C_PTR) :: pdll
-      TYPE(C_FUNPTR) :: pdesign=C_NULL_FUNPTR
-	  PROCEDURE(DESIGN), POINTER :: ptrdesign=>NULL()
-#endif
 
 C INPUT
 	INTEGER nobs,d(2),ny,nz,nx,nu,nt,ns(6),S(nobs,6)
@@ -92,16 +76,9 @@ C EXTERNAL FUNCTIONS
      3 Xdd(MAX(d(1),1),nx),Pdd(MAX(d(1),1),nx,nx),
      3 WORK((nx+2)*(nx+1)/2),FP(nx,nx),WORK1(64*nx),UP(nu) )
 
-#if defined(__CYGWIN32__) || defined(_WIN32)
-	  pdesign = getprocaddress(pdll, "design_"C)
+#if defined(ORIGDLL) || defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
 	  CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #else
-	  pdesign = DLSym(pdll, 'design_'//C_NULL_CHAR)
-	  IF(.NOT.C_ASSOCIATED(pdesign)) THEN
-		 WRITE(*,*) ' Error in dlsym: ', C_F_STRING(DLError())
-	  END IF
-	  CALL C_F_PROCPOINTER(CPTR=pdesign, FPTR=ptrdesign)
-	  CALL ptrdesign(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R)
 #endif
 
 	ykP(:,:) = 0.D0
diff --git a/slice.for b/slice.for
index 616a43354369cd445f9612e3351e8fd682a8c9d4..9d4f571e2979fa9732f371f9d3be3d71585fdfd8 100644
--- a/slice.for
+++ b/slice.for
@@ -22,14 +22,12 @@ C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C -------------------------------------------------------------
 	SUBROUTINE SLICE(it,nobs,d,ny,nz,nx,nu,ns,nt,S,yk,IYK,
-	1                 theta,thetaprior,tipo,pdll,NEVAL,XSIM)
+	1                 theta,thetaprior,tipo,NEVAL,XSIM)
 C INPUT
 	INTEGER it,nobs,d(2),ny,nz,nx,nu,ns(6),nt,S(nobs,6),
 	1 IYK(nobs,ny+1)
 	DOUBLE PRECISION yk(nobs,ny+nz),theta(nt),thetaprior(4)
 	CHARACTER*2 tipo
-	CHARACTER*1 fittizia
-	POINTER (pdll,fittizia)
 
 C OUTPUT
 	INTEGER NEVAL
@@ -51,7 +49,7 @@ C    THIS DEFINES THE SLICE S={x: z < ln(f(x))}
 C -------------------------------------------------------
 	theta(it) = XOLD
 	FXOLD=PTHETA(it,nobs,d,ny,nz,nx,nu,ns,nt,S,yk,IYK,
-	1             theta,thetaprior,tipo,pdll)
+	1             theta,thetaprior,tipo)
 	NEVAL = NEVAL + 1
       U = genunf(0.D0,1.D0)
 	Z = FXOLD + DLOG(U)
@@ -72,14 +70,14 @@ C	U = G05CAF(U)
 	 DO 100 WHILE(L.GT.XLB)
         theta(it) = L
  	  FXL=PTHETA(it,nobs,d,ny,nz,nx,nu,ns,nt,S,yk,IYK,
-	1             theta,thetaprior,tipo,pdll)
+	1             theta,thetaprior,tipo)
 	  NEVAL = NEVAL + 1
 	  IF (FXL.LE.Z) GOTO 110
 100	 L = L - W
 110     DO 200 WHILE(R.LT.XUB)
         theta(it) = R
  	  FXR=PTHETA(it,nobs,d,ny,nz,nx,nu,ns,nt,S,yk,IYK,
-	1             theta,thetaprior,tipo,pdll)
+	1             theta,thetaprior,tipo)
 	  NEVAL = NEVAL + 1
 	  IF (FXR.LE.Z) GOTO 210
 200	 R = R + W
@@ -93,7 +91,7 @@ C	 U = G05CAF(U)
 	  IF (L.LE.XLB) GOTO 310
 	  theta(it) = L
  	  FXL=PTHETA(it,nobs,d,ny,nz,nx,nu,ns,nt,S,yk,IYK,
-	1             theta,thetaprior,tipo,pdll)
+	1             theta,thetaprior,tipo)
 	  NEVAL = NEVAL + 1
 	  IF (FXL.LE.Z) GOTO 310
 	  L = L - W
@@ -103,7 +101,7 @@ C	 U = G05CAF(U)
 	  IF (R.GE.XLB) GOTO 410
 	  theta(it) = R
  	  FXR=PTHETA(it,nobs,d,ny,nz,nx,nu,ns,nt,S,yk,IYK,
-	1             theta,thetaprior,tipo,pdll)
+	1             theta,thetaprior,tipo)
 	  NEVAL = NEVAL + 1
 	  IF (FXR.LE.Z) GOTO 410
 	  R = R + W
@@ -123,7 +121,7 @@ C	 U = G05CAF(U)
 	 XSIM = L + U*(R - L)
 	 theta(it) = XSIM
 	 FXSIM=PTHETA(it,nobs,d,ny,nz,nx,nu,ns,nt,S,yk,IYK,
-	1              theta,thetaprior,tipo,pdll)
+	1              theta,thetaprior,tipo)
 	 NEVAL = NEVAL + 1
 	 IF (FXSIM.GE.Z) OK = 1
 	 IF(XSIM.GT.XOLD) THEN
diff --git a/slice2.for b/slice2.for
index 6810e700ab2387ff42e66865c4d1bc03d0d0f304..80d673ef0ec08a9763a27650d8a9ff2e538a5958 100644
--- a/slice2.for
+++ b/slice2.for
@@ -22,13 +22,11 @@ C You should have received a copy of the GNU General Public License
 C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
 C -------------------------------------------------------------
 	SUBROUTINE SLICE2(it,nobs,d,ny,nz,nx,nu,ns,nt,S,yk,
-	1                  theta,thetaprior,tipo,pdll,NEVAL,XSIM)
+	1                  theta,thetaprior,tipo,NEVAL,XSIM)
 C INPUT
 	INTEGER it,nobs,d(2),ny,nz,nx,nu,ns(6),nt,S(nobs,6)
 	DOUBLE PRECISION yk(nobs,ny+nz),theta(nt),thetaprior(4)
 	CHARACTER*2 tipo
-	CHARACTER*1 fittizia
-	POINTER (pdll,fittizia)
 
 C OUTPUT
 	INTEGER NEVAL
@@ -50,7 +48,7 @@ C    THIS DEFINES THE SLICE S={x: z < ln(f(x))}
 C -------------------------------------------------------
 	theta(it) = XOLD
 	FXOLD=PTHETA2(it,nobs,d,ny,nz,nx,nu,ns,nt,S,yk,
-	1              theta,thetaprior,tipo,pdll)
+	1              theta,thetaprior,tipo)
 	NEVAL = NEVAL + 1
       U = genunf(0.D0,1.D0)
 	Z = FXOLD + DLOG(U)
@@ -71,14 +69,14 @@ C	U = G05CAF(U)
 	 DO 100 WHILE(L.GT.XLB)
         theta(it) = L
  	  FXL=PTHETA2(it,nobs,d,ny,nz,nx,nu,ns,nt,S,yk,
-	1              theta,thetaprior,tipo,pdll)
+	1              theta,thetaprior,tipo)
 	  NEVAL = NEVAL + 1
 	  IF (FXL.LE.Z) GOTO 110
 100	   L = L - W
 110     DO 200 WHILE(R.LT.XUB)
         theta(it) = R
  	  FXR=PTHETA2(it,nobs,d,ny,nz,nx,nu,ns,nt,S,yk,
-	1              theta,thetaprior,tipo,pdll)
+	1              theta,thetaprior,tipo)
 	  NEVAL = NEVAL + 1
 	  IF (FXR.LE.Z) GOTO 210
 200	 R = R + W
@@ -92,7 +90,7 @@ C	 U = G05CAF(U)
 	  IF (L.LE.XLB) GOTO 310
 	  theta(it) = L
  	  FXL=PTHETA2(it,nobs,d,ny,nz,nx,nu,ns,nt,S,yk,
-	1              theta,thetaprior,tipo,pdll)
+	1              theta,thetaprior,tipo)
 	  NEVAL = NEVAL + 1
 	  IF (FXL.LE.Z) GOTO 310
 	  L = L - W
@@ -102,7 +100,7 @@ C	 U = G05CAF(U)
 	  IF (R.GE.XLB) GOTO 410
 	  theta(it) = R
  	  FXR=PTHETA2(it,nobs,d,ny,nz,nx,nu,ns,nt,S,yk,
-	1              theta,thetaprior,tipo,pdll)
+	1              theta,thetaprior,tipo)
 	  NEVAL = NEVAL + 1
 	  IF (FXR.LE.Z) GOTO 410
 	  R = R + W
@@ -122,7 +120,7 @@ C	 U = G05CAF(U)
 	 XSIM = L + U*(R - L)
 	 theta(it) = XSIM
 	 FXSIM=PTHETA2(it,nobs,d,ny,nz,nx,nu,ns,nt,S,yk,
-	1               theta,thetaprior,tipo,pdll)
+	1               theta,thetaprior,tipo)
 	 NEVAL = NEVAL + 1
 	 IF (FXSIM.GE.Z) OK = 1
 	 IF(XSIM.GT.XOLD) THEN