Index: trunk/tools/Makefile.llamp
===================================================================
--- trunk/tools/Makefile.llamp	(revision 1608)
+++ trunk/tools/Makefile.llamp	(revision 1608)
@@ -0,0 +1,90 @@
+# L. Fita, LMD. January 2017
+# Compilation of all required Fortran modules in `llamp' (debian 9) with gcc compilers
+# A following python 'library' creation is required
+# Diagnostics compilation _______
+# f2py -m module_ForDiag -h module_ForDiag.pyf module_generic.f90 module_ForDiagnosticsVars.f90 module_ForDiagnostics.f90 >& run_f2py_ForDiag_wrapper.log
+# make
+# f2py -m module_ForDiag -c --f90exec=/usr/bin/gfortran-6 module_ForDiag.pyf module_generic.f90 module_ForDiagnosticsVars.f90 module_ForDiagnostics.f90 >& run_f2py_ForDiag_compile.log
+# Interpolation compilation _______
+# f2py -m module_ForInt -h module_ForInt.pyf module_generic.f90 module_ForInterpolate.f90  >& run_f2py_ForInt_wrapper.log
+# f2py -m module_ForInt -c --f90exec=/usr/bin/gfortran-6 module_ForInt.pyf module_generic.f90 module_ForInterpolate.f90  >& run_f2py_ForInt_compile.log
+
+FC 		= gfortran
+FCF 	= -c
+LIB_INC = 
+RM 		= rm -f
+#DBGFLAGS = -g -Wall -Wextra -Warray-temporaries -Wconversion -fimplicit-none -fbacktrace -ffree-line-length-0 -fcheck=all -ffpe-trap=zero,overflow,underflow -finit-real=nan
+NCLIBFOLD	= /usr/lib/x86_64-linux-gnu
+NCINCFOLD	= /usr
+LIB_NETCDF	= -L$(NCLIBFOLD)/lib -lnetcdff -lnetcdf -I$(NCINCFOLD)/include
+
+FCFLAGS = $(FCF) $(DBGFLAGS)
+
+# Sources for f2py
+distrisrcs = module_definitions.f90 module_generic.f90 module_scientific.f90 module_DistriCorrection.f90
+diagsrcs = module_definitions.f90 module_generic.f90 module_scientific.f90 module_ForDiagnosticsVars.f90 module_ForDiagnostics.f90
+intsrcs = module_definitions.f90 module_generic.f90 module_scientific.f90 module_ForInterpolate.f90
+
+####### ###### ##### #### ### ## #
+
+MODULES = \
+	module_definitions.o \
+	module_scientific.o \
+	module_generic.o \
+	module_ForDiagnosticsVars.o \
+	module_ForDiagnostics.o \
+	module_DistriCorrection.o
+
+all : \
+	module_ForDiagnostics.o \
+	DistriCorrection.o \
+	interpolate.o \
+	pydistrimods.o \
+	pydiagmods.o \
+	pyintmods.o
+
+diags : \
+	module_ForDiagnostics.o
+
+clean :
+	$(RM) *.mod *.o interpolate module_ForDiag.so module_ForInt.so module_ForDistriCorrect.so
+
+########    #######
+##
+    ########
+
+module_definitions.o:
+	$(FC) $(FCFLAGS) module_definitions.f90
+
+module_generic.o: module_definitions.o
+	$(FC) $(FCFLAGS) $(LIB_NETCDF) $(LIB_INC) module_generic.f90
+
+module_scientific.o: module_definitions.o module_generic.o
+	$(FC) $(FCFLAGS) module_scientific.f90
+
+module_ForDiagnosticsVars.o: module_definitions.o module_generic.o
+	$(FC) $(FCFLAGS) $(LIB_INC) module_ForDiagnosticsVars.f90
+
+module_ForDiagnostics.o: module_definitions.o module_generic.o module_ForDiagnosticsVars.o
+	$(FC) $(FCFLAGS) $(LIB_INC) module_ForDiagnostics.f90
+
+module_DistriCorrection.o: module_definitions.o module_generic.o module_scientific.o
+	$(FC) $(FCFLAGS) $(LIB_NETCDF) module_DistriCorrection.f90
+
+FCEXECFLAGS = $(MODULES) $(LIB_NETCDF) $(LIB_INC)
+
+interpolate.o: module_definitions.o module_generic.o module_scientific.o
+	$(FC) $(FCEXECFLAGS) interpolate.f90 -o interpolate
+
+DistriCorrection.o: module_definitions.o module_generic.o module_scientific.o module_DistriCorrection.o
+	$(FC) $(FCEXECFLAGS) DistriCorrection.f90 -o DistriCorrection
+
+pydistrimods.o:
+	f2py -c -m module_ForDistriCorrect $(distrisrcs)
+
+pydiagmods.o:
+	f2py -c -m module_ForDiag $(diagsrcs)
+
+pyintmods.o:
+	f2py -c -m module_ForInt $(intsrcs)
+
Index: trunk/tools/interpolate.f90
===================================================================
--- trunk/tools/interpolate.f90	(revision 1607)
+++ trunk/tools/interpolate.f90	(revision 1608)
@@ -16,4 +16,5 @@
 ! Function to search a given value from a coarser version of the data
 
+  USE module_definitions
   USE module_generic
 
@@ -37,5 +38,4 @@
   REAL(r_k)                                              :: fraclonv, fraclatv
   REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: difflonlat, lon, lat
-  CHARACTER(LEN=50)                                      :: fname
   
 ! Variables
@@ -225,5 +225,4 @@
   REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: fractionlon, fractionlat
   INTEGER                                                :: fracdx, fracdy
-  CHARACTER(LEN=50)                                      :: fname
 
 !!!!!!! Variables
Index: trunk/tools/module_ForDiagnostics.f90
===================================================================
--- trunk/tools/module_ForDiagnostics.f90	(revision 1607)
+++ trunk/tools/module_ForDiagnostics.f90	(revision 1608)
@@ -7,4 +7,5 @@
 MODULE module_ForDiagnostics
 
+  USE module_definitions
   USE module_generic
   USE module_ForDiagnosticsVars
@@ -31,5 +32,4 @@
     IMPLICIT NONE
 
-    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
     INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
     REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: cldfra4D, pres4D
@@ -38,5 +38,4 @@
 ! Local
     INTEGER                                              :: i,j,k, zdim, Ndim
-    CHARACTER(LEN=50)                                    :: fname
 
 !!!!!!! Variables
@@ -71,5 +70,4 @@
     IMPLICIT NONE
 
-    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
     INTEGER, INTENT(in)                                  :: d1, d2, d3
     REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: cldfra3D, pres3D
@@ -78,5 +76,4 @@
 ! Local
     INTEGER                                              :: i,j,k, zdim, Ndim
-    CHARACTER(LEN=50)                                    :: fname
 
 !!!!!!! Variables
@@ -109,5 +106,4 @@
     IMPLICIT NONE
 
-    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
     INTEGER, INTENT(in)                                  :: Ndim, d1, d2, d3, d4, zdim
     REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in)       :: cldfra1D, pres1D
@@ -133,5 +129,4 @@
 ! Local
     INTEGER                                              :: i,j,k
-    CHARACTER(LEN=50)                                    :: fname
 
 !!!!!!! Variables
@@ -253,5 +248,4 @@
     IMPLICIT NONE
 
-    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
     INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
     REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: cldfra4D
@@ -260,5 +254,4 @@
 ! Local
     INTEGER                                              :: i,j,k, zdim, Ndim
-    CHARACTER(LEN=50)                                    :: fname
 
 !!!!!!! Variables
@@ -292,5 +285,4 @@
     IMPLICIT NONE
 
-    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
     INTEGER, INTENT(in)                                  :: d1, d2, d3
     REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: cldfra3D
@@ -299,5 +291,4 @@
 ! Local
     INTEGER                                              :: i,j,k, zdim, Ndim
-    CHARACTER(LEN=50)                                    :: fname
 
 !!!!!!! Variables
@@ -328,5 +319,4 @@
     IMPLICIT NONE
 
-    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
     INTEGER, INTENT(in)                                  :: Ndim, d1, d2, d3, d4, zdim
     REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in)       :: cldfra1D
@@ -348,5 +338,4 @@
 ! Local
     INTEGER                                              :: i,j,k
-    CHARACTER(LEN=50)                                    :: fname
 
 !!!!!!! Variables
Index: trunk/tools/module_ForDiagnosticsVars.f90
===================================================================
--- trunk/tools/module_ForDiagnosticsVars.f90	(revision 1607)
+++ trunk/tools/module_ForDiagnosticsVars.f90	(revision 1608)
@@ -7,9 +7,9 @@
 MODULE module_ForDiagnosticsVars
 
+  USE module_definitions
   USE module_generic
 
   IMPLICIT NONE
 
-!  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
   REAL(r_k), PARAMETER                                   :: ZEPSEC=1.0D-12
 ! Low limit pressure for medium clouds [Pa]
@@ -17,8 +17,4 @@
 ! Low limit pressure for High clouds [Pa]
   REAL(r_k), PARAMETER                                   :: prmlc = 68000.d0
-
-  REAL(r_k), PARAMETER                                   :: zero=0.d0
-  REAL(r_k), PARAMETER                                   :: one=1.d0
-  REAL(r_k), PARAMETER                                   :: two=2.d0
 
   CONTAINS
@@ -40,5 +36,4 @@
     IMPLICIT NONE
 
-!    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
     INTEGER, INTENT(in)                                  :: dz
     REAL(r_k), DIMENSION(dz), INTENT(in)                 :: clfra, p
@@ -47,5 +42,4 @@
 ! Local
     INTEGER                                              :: iz
-    CHARACTER(LEN=50)                                    :: fname
     REAL(r_k)                                            :: zclearl, zcloudl, zclearm, zcloudm,       &
       zclearh, zcloudh
@@ -56,22 +50,22 @@
     fname = 'var_cllmh'
 
-    zclearl = one
-    zcloudl = zero
-    zclearm = one
-    zcloudm = zero
-    zclearh = one
-    zcloudh = zero
+    zclearl = oneRK
+    zcloudl = zeroRK
+    zclearm = oneRK
+    zcloudm = zeroRK
+    zclearh = oneRK
+    zcloudh = zeroRK
 
-    var_cllmh = one
+    var_cllmh = oneRK
 
     DO iz=1, dz
       IF (p(iz) < prmhc) THEN
-        var_cllmh(3) = var_cllmh(3)*(one-MAX(clfra(iz),zcloudh))/(one-MIN(zcloudh,one-ZEPSEC))
+        var_cllmh(3) = var_cllmh(3)*(oneRK-MAX(clfra(iz),zcloudh))/(oneRK-MIN(zcloudh,oneRK-ZEPSEC))
         zcloudh = clfra(iz)
       ELSE IF ( (p(iz) >= prmhc) .AND. (p(iz) < prmlc) ) THEN
-        var_cllmh(2) = var_cllmh(2)*(one-MAX(clfra(iz),zcloudm))/(one-MIN(zcloudm,one-ZEPSEC))
+        var_cllmh(2) = var_cllmh(2)*(oneRK-MAX(clfra(iz),zcloudm))/(oneRK-MIN(zcloudm,oneRK-ZEPSEC))
         zcloudm = clfra(iz)
       ELSE IF (p(iz) >= prmlc) THEN
-        var_cllmh(1) = var_cllmh(1)*(one-MAX(clfra(iz),zcloudl))/(one-MIN(zcloudl,one-ZEPSEC))
+        var_cllmh(1) = var_cllmh(1)*(oneRK-MAX(clfra(iz),zcloudl))/(oneRK-MIN(zcloudl,oneRK-ZEPSEC))
         zcloudl = clfra(iz)
       ELSE
@@ -85,5 +79,5 @@
     END DO
 
-    var_cllmh = one - var_cllmh
+    var_cllmh = oneRK - var_cllmh
 
     RETURN 
@@ -97,5 +91,4 @@
     IMPLICIT NONE
 
-!    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
     INTEGER, INTENT(in)                                  :: dz
     REAL(r_k), DIMENSION(dz), INTENT(in)                 :: clfra
@@ -103,5 +96,5 @@
     INTEGER                                              :: iz
     REAL(r_k)                                            :: zclear, zcloud
-    CHARACTER(LEN=50)                                    :: fname
+
 !!!!!!! Variables
 ! cfra: 1-column cloud fraction values
@@ -109,10 +102,10 @@
     fname = 'var_clt'
 
-    zclear = one
-    zcloud = zero
+    zclear = oneRK
+    zcloud = zeroRK
 
     DO iz=1,dz
-      zclear = zclear*(one-MAX(clfra(iz),zcloud))/(one-MIN(zcloud,1.-ZEPSEC))
-      var_clt = one - zclear
+      zclear = zclear*(oneRK-MAX(clfra(iz),zcloud))/(oneRK-MIN(zcloud,1.-ZEPSEC))
+      var_clt = oneRK - zclear
       zcloud = clfra(iz)
     END DO
Index: trunk/tools/module_ForInterpolate.f90
===================================================================
--- trunk/tools/module_ForInterpolate.f90	(revision 1607)
+++ trunk/tools/module_ForInterpolate.f90	(revision 1608)
@@ -10,4 +10,5 @@
 MODULE module_ForInterpolate
 
+  USE module_definitions
   USE module_generic
 
@@ -37,5 +38,4 @@
   REAL(r_k)                                              :: fraclonv, fraclatv
   REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: difflonlat, lon, lat
-  CHARACTER(LEN=50)                                      :: fname
   
 ! Variables
@@ -227,5 +227,4 @@
   REAL(r_k)                                              :: fraclonv, fraclatv
   REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: difflonlat, lon, lat
-  CHARACTER(LEN=50)                                      :: fname
   
 ! Variables
@@ -446,5 +445,4 @@
 ! Local
   REAL(r_k), DIMENSION(dx,dy)                            :: difflonlat
-  CHARACTER(LEN=50)                                      :: fname
   
 ! Variables
@@ -503,5 +501,4 @@
   REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: fractionlon, fractionlat
   INTEGER                                                :: dfracdx, dfracdy, fracdx, fracdy
-  CHARACTER(LEN=50)                                      :: fname
 
 !!!!!!! Variables
@@ -609,5 +606,4 @@
   REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: fractionlon, fractionlat
   INTEGER                                                :: dfracdx, dfracdy, fracdx, fracdy
-  CHARACTER(LEN=50)                                      :: fname
 
 !!!!!!! Variables
@@ -735,5 +731,4 @@
   REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat, ipos
   INTEGER, DIMENSION(2)                                  :: iLl
-  CHARACTER(LEN=50)                                      :: fname
 
 !!!!!!! Variables
@@ -889,5 +884,4 @@
   REAL(r_k)                                              :: w
   REAL(r_k), DIMENSION(3,16,pdimx,pdimy)                 :: outLlw
-  CHARACTER(LEN=50)                                      :: fname
 
 !!!!!!! Variables
@@ -967,5 +961,4 @@
   REAL(r_k)                                              :: w
   REAL(r_k), DIMENSION(3,16,pdimx,pdimy)                 :: outLlw
-  CHARACTER(LEN=50)                                      :: fname
 
 !!!!!!! Variables
@@ -1048,5 +1041,4 @@
   REAL(r_k)                                              :: w
   REAL(r_k), DIMENSION(3,16,pdimx,pdimy)                 :: outLlw
-  CHARACTER(LEN=50)                                      :: fname
 
 !!!!!!! Variables
@@ -1133,5 +1125,4 @@
   REAL(r_k)                                              :: w
   REAL(r_k), DIMENSION(3,16,pdimx,pdimy)                 :: outLlw
-  CHARACTER(LEN=50)                                      :: fname
 
 !!!!!!! Variables
@@ -1225,5 +1216,4 @@
   REAL(r_k), DIMENSION(dimx,dimy)                        :: difflonlat
   REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
-  CHARACTER(LEN=50)                                      :: fname
 
 !!!!!!! Variables
@@ -1324,5 +1314,4 @@
   REAL(r_k), DIMENSION(dimy)                             :: difflat
   REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
-  CHARACTER(LEN=50)                                      :: fname
 
 !!!!!!! Variables
Index: trunk/tools/module_generic.f90
===================================================================
--- trunk/tools/module_generic.f90	(revision 1607)
+++ trunk/tools/module_generic.f90	(revision 1608)
@@ -1,15 +1,74 @@
 MODULE module_generic
 ! Module with generic functions
+
+!!!!!!! Subroutines/Functions
+! ErrMsg: Subroutine to stop execution and provide an error message
+! ErrWarnMsg: Function to print error/warning message
 ! Index1DArrayR: Function to provide the first index of a given value inside a 1D real array
+! Index1DArrayR_K: Function to provide the first index of a given value inside a 1D real(r_k) array
 ! Index2DArrayR: Function to provide the first index of a given value inside a 2D real array
-! ErrWarnMsg: Function to print error/warning message
-
-  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
-  ! Fill value at 64 bits
-  REAL(r_k)                                              :: fillval64 = 1.e20
-  CHARACTER(len=50), PARAMETER                           :: errormsg = 'ERROR -- error -- ERROR -- error'
-  CHARACTER(len=50), PARAMETER                           :: warnmsg = 'WARNING -- warning -- WARNING -- warning'
+! Nstrings: Function to repeat a number of times a given string
+! RangeI: Function to provide a range of d1 values from 'iniv' to 'endv', of integer values in a vector
+! RangeR: Function to provide a range of d1 values from 'iniv' to 'endv', of real values in a vector
+! RangeR_K: Function to provide a range of d1 from 'iniv' to 'endv', of real(r_k) values in a vector
+! vectorR_KS: Function to transform a vector of reals to a string of characters
+
+  USE module_definitions
 
   CONTAINS
+
+  CHARACTER(len=1000) FUNCTION vectorR_KS(d1, vector)
+  ! Function to transform a vector of reals(r_k) to a string of characters
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: d1
+    REAL(r_k), DIMENSION(d1), INTENT(in)                 :: vector
+
+! Local
+    INTEGER                                              :: iv
+    CHARACTER(len=50)                                    :: RS
+
+!!!!!!! Variables
+! d1: length of the vector
+! vector: values to transform
+
+    fname = 'vectorR_KS'
+
+    vectorR_KS = ''
+    DO iv=1, d1
+      WRITE(RS, '(F50.25)')vector(iv)
+      IF (iv == 1) THEN
+        vectorR_KS = TRIM(RS)
+      ELSE
+        vectorR_KS = TRIM(vectorR_KS) // ', ' // TRIM(RS)
+      END IF
+    END DO    
+
+  END FUNCTION vectorR_KS
+
+CHARACTER(len=1000) FUNCTION Nstrings(Strval, Ntimes)
+! Function to repeat a number of times a given string
+
+  IMPLICIT NONE
+
+  CHARACTER(LEN=50), INTENT(in)                        :: Strval
+  INTEGER, INTENT(in)                                  :: Ntimes
+
+! Local
+  INTEGER                                              :: i
+
+!!!!!!! Variables
+! Strval: String to repeat
+! Ntimes: number of repetitions
+
+  fname = 'Nstrings'
+
+  Nstrings = ''
+  Do i=1, Ntimes
+    Nstrings = TRIM(Nstrings) // TRIM(Strval)
+  END DO
+
+END FUNCTION Nstrings
 
   INTEGER FUNCTION Index1DArrayR(array1D, d1, val)
@@ -25,5 +84,4 @@
 ! Local
     INTEGER                                              :: i
-    CHARACTER(LEN=50)                                    :: fname
 
     fname = 'Index1DArrayR'
@@ -40,4 +98,29 @@
   END FUNCTION Index1DArrayR
 
+  INTEGER FUNCTION Index1DArrayR_K(array1D, d1, val)
+! Function to provide the first index of a given value inside a 1D real(r_k) array
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: d1
+    REAL(r_k), INTENT(in)                                :: val
+    REAL(r_k), DIMENSION(d1), INTENT(in)                 :: array1D
+
+! Local
+    INTEGER                                              :: i
+
+    fname = 'Index1DArrayR_K'
+
+    Index1DArrayR_K = -1
+
+    DO i=1,d1
+      IF (array1d(i) == val) THEN
+        Index1DArrayR_K = i
+        EXIT
+      END IF
+    END DO
+
+  END FUNCTION Index1DArrayR_K
+
   FUNCTION Index2DArrayR(array2D, d1, d2, val)
 ! Function to provide the first index of a given value inside a 2D real array
@@ -53,5 +136,4 @@
 ! Local
     INTEGER                                              :: i, j
-    CHARACTER(LEN=50)                                    :: fname
 
     fname = 'Index2DArrayR'
@@ -71,4 +153,105 @@
   END FUNCTION Index2DArrayR
 
+  FUNCTION RangeI(d1, iniv, endv)
+! Function to provide a range of d1 values from 'iniv' to 'endv', of integer values in a vector
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: d1, iniv, endv
+    INTEGER, DIMENSION(d1)                               :: RangeI
+
+! Local
+    INTEGER                                              :: i, intv
+
+    fname = 'RangeI'
+
+    intv = (endv - iniv) / (d1*1 - 1)
+
+    RangeI(1) = iniv
+    DO i=2,d1
+      RangeI(i) = RangeI(i-1) + intv
+    END DO
+
+  END FUNCTION RangeI 
+
+  FUNCTION RangeR(d1, iniv, endv)
+! Function to provide a range of d1 from 'iniv' to 'endv', of real values in a vector
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: d1
+    REAL, INTENT(in)                                     :: iniv, endv
+    REAL, DIMENSION(d1)                                  :: RangeR
+
+! Local
+    INTEGER                                              :: i
+    REAL                                                 :: intv
+
+    fname = 'RangeR'
+
+    intv = (endv - iniv) / (d1*1. - 1.)
+
+    RangeR(1) = iniv
+    DO i=2,d1
+      RangeR(i) = RangeR(i-1) + intv
+    END DO
+
+  END FUNCTION RangeR
+
+
+  FUNCTION RangeR_K(d1, iniv, endv)
+! Function to provide a range of d1 from 'iniv' to 'endv', of real(r_k) values in a vector
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: d1
+    REAL(r_k), INTENT(in)                                :: iniv, endv
+    REAL(r_k), DIMENSION(d1)                             :: RangeR_K
+
+! Local
+    INTEGER                                              :: i
+    REAL(r_k)                                            :: intv
+
+    fname = 'RangeR_K'
+
+    intv = (endv - iniv) / (d1*oneRK-oneRK)
+
+    RangeR_K(1) = iniv
+    DO i=2,d1
+      RangeR_K(i) = RangeR_K(i-1) + intv
+    END DO
+
+  END FUNCTION RangeR_K
+
+
+SUBROUTINE ErrMsg(msg, funcn, errN)
+! Subroutine to stop execution and provide an error message
+
+  IMPLICIT NONE
+
+  CHARACTER(LEN=*), INTENT(in)                           :: msg, funcn
+  INTEGER, INTENT(in)                                    :: errN
+
+! Local
+  CHARACTER(LEN=50)                                      :: emsg
+
+!!!!!!! Variables
+! msg: message to print with the error
+! funcn: name of the funciton, section to localize the error
+! errN: number of the error provided for a given FORTRAN function
+
+  emsg = 'ERORR -- error -- ERROR -- error'
+
+  IF (errN /= 0) THEN
+    PRINT *,TRiM(emsg)
+    PRINT *,'  ' // TRIM(funcn) // ': ' // TRIM(msg) // ' !!'
+    PRINT *,'    error number:', errN
+    STOP
+  END IF
+
+  RETURN
+
+END SUBROUTINE ErrMsg
+
   CHARACTER(LEN=50) FUNCTION ErrWarnMsg(msg)
 ! Function to print error/warning message
@@ -78,5 +261,4 @@
     CHARACTER(LEN=3), INTENT(in)                         :: msg
 ! Local
-    CHARACTER(LEN=50)                                    :: fname
 
     fname = 'ErrWarnMsg'
