diff --git a/fftma_module/gen/include/geostat.h b/fftma_module/gen/include/geostat.h
index 33cf755..f1d2922 100644
--- a/fftma_module/gen/include/geostat.h
+++ b/fftma_module/gen/include/geostat.h
@@ -327,7 +327,7 @@ void coordinates(int maille, int i[3], struct grid_mod grid);
 /*variogram: structure defined above     */
 /*grid: structure defined above          */
 /*n: number of gridblocks along X,Y and Z*/
-void covariance(double* covar, struct vario_mod variogram, struct grid_mod grid, int n[3], int cores);
+void covariance(chunk_array_t* covar, struct vario_mod variogram, struct grid_mod grid, int n[3], int cores);
 
 /*computation of the covariance matrix for the well data*/
 /*well coordinates are given as a number of cells       */
@@ -403,7 +403,7 @@ double exponential(double h);
 /*workr: utility real part vector for storage         */
 /*worki: utility imaginary part vector for storage    */
 /*The transformed data are returned in datar and datai*/
-void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores);
+void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores);
 
 /*calculates F(x) = (1/a)*exp(-x*x/2)*/
 double funtrun1(double x);
diff --git a/fftma_module/gen/include/toolsFFTMA.h b/fftma_module/gen/include/toolsFFTMA.h
index 1731ebb..42af6d4 100644
--- a/fftma_module/gen/include/toolsFFTMA.h
+++ b/fftma_module/gen/include/toolsFFTMA.h
@@ -51,7 +51,7 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r
 /*        must be a Gaussian white noise       */
 /*realization: structure defining a realization*/
 
-void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, double* realization, int solver, int cores, long* seed);
+void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, chunk_array_t* realization, int solver, int cores, long* seed);
 
 /* build_real   */
 /* build a realization in the spectral domain */
@@ -65,8 +65,8 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin
 /*realization: vector defining the real part   */
 /*ireal: vector defining the i-part   */
 
-void build_real(int n[3], int NTOT, double* covar, double* realization, double* ireal, int cores);
+void build_real(int n[3], int NTOT, chunk_array_t* covar, chunk_array_t* realization, double* ireal, int cores);
 
-void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, double* vectorresult, struct realization_mod* realout);
+void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, chunk_array_t* vectorresult, struct realization_mod* realout, int cores);
 
 #endif // define _TOOLSFFTMA_H
diff --git a/fftma_module/gen/lib_src/Py_kgeneration.c b/fftma_module/gen/lib_src/Py_kgeneration.c
index c79ad46..a20469a 100644
--- a/fftma_module/gen/lib_src/Py_kgeneration.c
+++ b/fftma_module/gen/lib_src/Py_kgeneration.c
@@ -27,6 +27,8 @@ void Py_kgeneration(long seed, struct grid_mod grid, struct statistic_mod stat,
     n[1] = 0;
     n[2] = 0;
 
+    printf("aaa");
+
     generate(&seed, N, Z, cores);
 
     /*FFTMA*/
diff --git a/fftma_module/gen/lib_src/build_real.c b/fftma_module/gen/lib_src/build_real.c
index f563f85..6274e3b 100644
--- a/fftma_module/gen/lib_src/build_real.c
+++ b/fftma_module/gen/lib_src/build_real.c
@@ -1,4 +1,5 @@
 #include "geostat.h"
+#include "chunk_array.h"
 #include <math.h>
 #include <stdarg.h>
 #include <stddef.h>
@@ -19,22 +20,26 @@
 /*realization: vector defining the real part   */
 /*ireal: vector defining the i-part   */
 
-void build_real(int n[3], int NTOT, double* covar, double* realization, double* ireal, int cores) {
+void build_real(int n[3], int NTOT, chunk_array_t* covar, chunk_array_t* realization, double* ireal, int cores) {
     int i, j, k, maille1;
     double temp;
 
+    chunk_array_read(realization);
+
     /*decomposition and multiplication in the spectral domain*/
     for (k = 1; k <= n[2]; k++) {
         for (j = 1; j <= n[1]; j++) {
             for (i = 1; i <= n[0]; i++) {
                 maille1 = i + (j - 1 + (k - 1) * n[1]) * n[0];
-                temp = covar[maille1];
+                chunk_aray_get(covar, maille1, &temp);
                 if (temp > 0.) {
                     temp = sqrt(temp) / (double)NTOT;
                 } else if (temp < 0.) {
                     temp = sqrt(-temp) / (double)NTOT;
                 }
-                realization[maille1] *= temp;
+                double valuemaille1;
+                chunk_array_get(realization, maille1, &valuemaille1);
+                chunk_array_save(realization, maille1, temp * valuemaille1);
                 ireal[maille1] *= temp;
             }
         }
diff --git a/fftma_module/gen/lib_src/clean_real.c b/fftma_module/gen/lib_src/clean_real.c
index 9917554..a6de2b7 100644
--- a/fftma_module/gen/lib_src/clean_real.c
+++ b/fftma_module/gen/lib_src/clean_real.c
@@ -1,4 +1,5 @@
 #include "geostat.h"
+#include "chunk_array.h"
 #include <math.h>
 #include <stdarg.h>
 #include <stddef.h>
@@ -7,7 +8,7 @@
 #include <string.h>
 #include <time.h>
 
-void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, double* vectorresult, struct realization_mod* realout, int cores) {
+void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid, chunk_array_t* vectorresult, struct realization_mod* realout, int cores) {
     int i, j, k, maille0, maille1;
     double NTOT;
 
@@ -15,6 +16,8 @@ void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid,
     /*is the output realization already allocated?*/
     /*if not, memory allocation*/
 
+    chunk_array_read(vectorresult);
+
     if (realout->vector == NULL || realout->n != realin->n) {
         realout->vector = (double*)malloc(realin->n * sizeof(double));
         if (realout->vector == NULL) {
@@ -31,7 +34,9 @@ void clean_real(struct realization_mod* realin, int n[3], struct grid_mod grid,
                 maille0 = i - 1 + (j - 1 + (k - 1) * grid.NY) * grid.NX;
                 /* Modif du 18 juin 2003 */
                 /*realout->vector[maille0] = vectorresult[maille1]/(double) NTOT;*/
-                realout->vector[maille0] = vectorresult[maille1];
+                double valuemaille1;
+                chunk_array_get(vectorresult, maille1, &valuemaille1);
+                realout->vector[maille0] = valuemaille1;
             }
         }
     }
diff --git a/fftma_module/gen/lib_src/covariance.c b/fftma_module/gen/lib_src/covariance.c
index 1bc9064..97ef87a 100644
--- a/fftma_module/gen/lib_src/covariance.c
+++ b/fftma_module/gen/lib_src/covariance.c
@@ -1,15 +1,18 @@
 #include "geostat.h"
+#include "chunk_array.h"
 #include <time.h>
 
 /*builds the sampled covariance function*/
 /*dimensions are even*/
-void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh, int n[3], int cores) {
+void covariance(chunk_array_t* covar, struct vario_mod variogram, struct grid_mod mesh, int n[3], int cores) {
     int i, j, k, maille, n2[3], symmetric;
     double di, dj, dk;
 
     for (i = 0; i < 3; i++)
         n2[i] = n[i] / 2;
 
+    chunk_array_read(covar);
+
     for (i = 0; i <= n2[0]; i++) {
         for (j = 0; j <= n2[1]; j++) {
             for (k = 0; k <= n2[2]; k++) {
@@ -19,12 +22,14 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh,
                 di = (double)i * mesh.DX;
                 dj = (double)j * mesh.DY;
                 dk = (double)k * mesh.DZ;
-                covar[maille] = (double)cov_value(variogram, di, dj, dk, cores);
+                chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores));
 
                 if (k > 0 && k < n2[2] && j > 0 && j < n2[1] && i > 0 && i < n2[0]) {
                     /*area 2*/
                     symmetric = 1 + n[0] - i + n[0] * (n[1] - j + n[1] * (n[2] - k));
-                    covar[symmetric] = covar[maille];
+                    double value;
+                    chunk_array_get(covar, maille, &value);
+                    chunk_array_save(covar, symmetric, value);
                 }
 
                 if (i > 0 && i < n2[0]) {
@@ -33,13 +38,15 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh,
                     dj = (double)j * mesh.DY;
                     dk = (double)k * mesh.DZ;
                     maille = 1 + (n[0] - i) + n[0] * (j + n[1] * k);
-                    covar[maille] = (double)cov_value(variogram, di, dj, dk, cores);
+                    chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores));
                 }
 
                 if (k > 0 && k < n2[2] && j > 0 && j < n2[1]) {
                     /*area 8*/
                     symmetric = 1 + i + n[0] * (n[1] - j + n[1] * (n[2] - k));
-                    covar[symmetric] = covar[maille];
+                    double value;
+                    chunk_array_get(covar, maille, &value);
+                    chunk_array_save(covar, symmetric, value);               
                 }
 
                 if (i > 0 && i < n2[0] && j > 0 && j < n2[1]) {
@@ -48,13 +55,15 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh,
                     dj = -(double)j * mesh.DY;
                     dk = (double)k * mesh.DZ;
                     maille = 1 + (n[0] - i) + n[0] * (n[1] - j + n[1] * k);
-                    covar[maille] = (double)cov_value(variogram, di, dj, dk, cores);
+                    chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores));
                 }
 
                 if (k > 0 && k < n2[2]) {
                     /*area 6*/
                     symmetric = 1 + i + n[0] * (j + n[1] * (n[2] - k));
-                    covar[symmetric] = covar[maille];
+                    double value;
+                    chunk_array_get(covar, maille, &value);
+                    chunk_array_save(covar, symmetric, value);                  
                 }
 
                 if (j > 0 && j < n2[1]) {
@@ -63,13 +72,15 @@ void covariance(double* covar, struct vario_mod variogram, struct grid_mod mesh,
                     dj = -(double)j * mesh.DY;
                     dk = (double)k * mesh.DZ;
                     maille = 1 + i + n[0] * (n[1] - j + n[1] * k);
-                    covar[maille] = (double)cov_value(variogram, di, dj, dk, cores);
+                    chunk_array_save(covar, maille, (double)cov_value(variogram, di, dj, dk, cores));
                 }
 
                 if (k > 0 && k < n2[2] && i > 0 && i < n2[0]) {
                     /*area 7*/
                     symmetric = 1 + n[0] - i + n[0] * (j + n[1] * (n[2] - k));
-                    covar[symmetric] = covar[maille];
+                    double value;
+                    chunk_array_get(covar, maille, &value);
+                    chunk_array_save(covar, symmetric, value);                      
                 }
             }
         }
diff --git a/fftma_module/gen/lib_src/fftma2.c b/fftma_module/gen/lib_src/fftma2.c
index d699498..5c8b224 100644
--- a/fftma_module/gen/lib_src/fftma2.c
+++ b/fftma_module/gen/lib_src/fftma2.c
@@ -1,4 +1,5 @@
 #include "geostat.h"
+#include "chunk_array.h"
 #include <math.h>
 #include <stdlib.h>
 #include <string.h>
@@ -26,7 +27,10 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r
     int NTOT, i, j, k, NMAX, NDIM, ntot, nmax, NXYZ, nxyz;
     int solver;
     double temp;
-    double *ireal, *covar, *workr, *worki, *realization;
+    chunk_array_t *covar, *ireal, *realization;
+    double *workr, *worki;
+
+    printf("Ayuda");
 
     /*covariance axis normalization*/
     axes(variogram.ap, variogram.scf, variogram.Nvario);
@@ -50,14 +54,9 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r
     nxyz = NXYZ + 1;
 
     /*array initialization*/
-    covar = (double*)malloc(ntot * sizeof(double));
-    testmemory(covar);
-
-    ireal = (double*)malloc(ntot * sizeof(double));
-    testmemory(ireal);
-
-    realization = (double*)malloc(ntot * sizeof(double));
-    testmemory(realization);
+    covar = chunk_array_create("covar.txt", ntot, 512);
+    ireal = chunk_array_create("ireal.txt", ntot, 512);
+    realization = chunk_array_create("realization.txt", ntot, 512);
 
     workr = (double*)malloc(nmax * sizeof(double));
     testmemory(workr);
@@ -83,12 +82,12 @@ void FFTMA2(struct vario_mod variogram, struct grid_mod grid, int n[3], struct r
     /* build realization in spectral domain   */
     build_real(n, NTOT, covar, realization, ireal, cores);
 
-    free(covar);
+    chunk_array_free(covar);
 
     /*backward fourier transform*/
     fourt(realization, ireal, n, NDIM, 0, 1, workr, worki, cores);
 
-    free(ireal);
+    chunk_array_free(ireal);
     free(workr);
     free(worki);
 
diff --git a/fftma_module/gen/lib_src/fourt.c b/fftma_module/gen/lib_src/fourt.c
index fc4184c..9d4ed0f 100644
--- a/fftma_module/gen/lib_src/fourt.c
+++ b/fftma_module/gen/lib_src/fourt.c
@@ -1,6 +1,7 @@
 #include <math.h>
 #include <stdio.h>
 #include <time.h>
+#include "chunk_array.h"
 
 /*fast fourier transform                                                   */
 /*     THE COOLEY-TUKEY FAST FOURIER TRANSFORM                             */
@@ -90,7 +91,7 @@
 /*     PROGRAM MODIFIED FROM A SUBROUTINE OF BRENNER                       */
 /*     10-06-2000, MLR                                                     */
 
-void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores) {
+void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores) {
     int ifact[21], ntot, idim, np1, n, np2, m, ntwo, iff, idiv, iquot, irem, inon2, non2p, np0, nprev, icase, ifmin, i, j, jmax, np2hf, i2, i1max, i3, j3, i1, ifp1, ifp2, i2max, i1rng, istep, imin, imax, mmax, mmin, mstep, j1, j2max, j2, jmin, j3max, nhalf;
     double theta, wstpr, wstpi, wminr, wmini, wr, wi, wtemp, thetm, wmstr, wmsti, twowr, sr, si, oldsr, oldsi, stmpr, stmpi, tempr, tempi, difi, difr, sumr, sumi, TWOPI = 6.283185307179586476925286766559;
 
@@ -99,6 +100,9 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
         ntot *= nn[idim];
     }
 
+    chunk_array_read(datar);
+    chunk_array_read(datai);
+
     /*main loop for each dimension*/
     np1 = 1;
     for (idim = 1; idim <= ndim; idim++) {
@@ -184,8 +188,11 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
         ntot /= 2;
         i = 1;
         for (j = 1; j <= ntot; j++) {
-            datar[j] = datar[i];
-            datai[j] = datar[i + 1];
+            double valuei, valuei1;
+            chunk_array_get(datar, i, &valuei);
+            chunk_array_get(datar, i+1, &valuei1);
+            chunk_array_save(datar, j, valuei);
+            chunk_array_save(datai, j, valuei1);
             i += 2;
         }
 
@@ -202,12 +209,20 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
             for (i1 = i2; i1 <= i1max; i1++) {
                 for (i3 = i1; i3 <= ntot; i3 += np2) {
                     j3 = j + i3 - i2;
-                    tempr = datar[i3];
-                    tempi = datai[i3];
-                    datar[i3] = datar[j3];
-                    datai[i3] = datai[j3];
-                    datar[j3] = tempr;
-                    datai[j3] = tempi;
+                    double valueri3, valueii3, valuerj3, valueij3;
+
+                    chunk_array_get(datar, i3, &valueri3);
+                    chunk_array_get(datai, i3, &valueii3);
+                    chunk_array_get(datar, j3, &valuerj3);
+                    chunk_array_get(datai, j3, &valueij3);
+
+                    tempr = valueri3;
+                    tempi = valueii3;
+
+                    chunk_array_save(datar, i3, valuerj3);
+                    chunk_array_save(datai, i3, valueij3);
+                    chunk_array_save(datar, j3, tempr);
+                    chunk_array_save(datai, j3, tempi);
                 }
             }
 
@@ -232,11 +247,14 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
             for (i3 = i1; i3 <= ntot; i3 += np2) {
                 j = i3;
                 for (i = 1; i <= n; i++) {
+                    double valuerj, valueij;
+                    chunk_array_get(datar, j, &valuerj);
+                    chunk_array_get(datai, j, &valueij);
                     if (icase != 3) {
-                        workr[i] = datar[j];
-                        worki[i] = datai[j];
+                        workr[i] = valuerj;
+                        worki[i] = valueij;
                     } else {
-                        workr[i] = datar[j];
+                        workr[i] = valuerj;
                         worki[i] = 0.;
                     }
                     ifp2 = np2;
@@ -255,8 +273,8 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
                 i2max = i3 + np2 - np1;
                 i = 1;
                 for (i2 = i3; i2 <= i2max; i2 += np1) {
-                    datar[i2] = workr[i];
-                    datai[i2] = worki[i];
+                    chunk_save(datar, i2, workr[i]);
+                    chunk_save(datai, i2, worki[i]);
                     i++;
                 }
             }
@@ -277,12 +295,16 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
         L310:
             j = i1;
             for (i = imin; i <= ntot; i += istep) {
-                tempr = datar[i];
-                tempi = datai[i];
-                datar[i] = datar[j] - tempr;
-                datai[i] = datai[j] - tempi;
-                datar[j] = datar[j] + tempr;
-                datai[j] = datai[j] + tempi;
+                chunk_array_get(datar, i, &tempr);
+                chunk_array_get(datai, i, &tempi);
+
+                double valuerj, valueij;
+                chunk_array_get(datar, j, &valuerj);
+                chunk_array_get(datai, j, &valueij);
+                chunk_array_save(datar, i, valuerj - tempr);
+                chunk_array_save(datai, i, valueij - tempi);
+                chunk_array_save(datar, j, valuerj + tempr);
+                chunk_array_save(datai, j, valueij + tempi);
                 j += istep;
             }
             imin = 2 * imin - i1;
@@ -301,16 +323,26 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
             j = imin - istep / 2;
             for (i = imin; i <= ntot; i += istep) {
                 if (ifrwd != 0) {
-                    tempr = datai[i];
-                    tempi = -datar[i];
+                    chunk_array_get(datai, i, &tempr);
+                    chunk_array_get(datar, i, &tempi);
+                    tempi = -tempi;
                 } else {
-                    tempr = -datai[i];
-                    tempi = datar[i];
+                    chunk_array_get(datai, i, &tempr);
+                    chunk_array_get(datar, i, &tempi);
                 }
-                datar[i] = datar[j] - tempr;
-                datai[i] = datai[j] - tempi;
-                datar[j] += tempr;
-                datai[j] += tempi;
+                double valuerj, valueij;
+                chunk_array_get(datar, j, &valuerj);
+                chunk_array_get(datai, j, &valueij);
+
+                chunk_array_save(datar, i, valuerj - tempr);
+                chunk_array_save(datai, i, valueij - tempi);
+
+                chunk_array_get(datar, j, &valuerj);
+                chunk_array_get(datai, j, &valueij);
+
+                chunk_array_save(datar, j, valuerj + tempr);
+                chunk_array_save(datai, j, valueij + tempi);
+
                 j += istep;
             }
 
@@ -347,12 +379,18 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
             L510:
                 j = imin - istep / 2;
                 for (i = imin; i <= ntot; i += istep) {
-                    tempr = datar[i] * wr - datai[i] * wi;
-                    tempi = datar[i] * wi + datai[i] * wr;
-                    datar[i] = datar[j] - tempr;
-                    datai[i] = datai[j] - tempi;
-                    datar[j] += tempr;
-                    datai[j] += tempi;
+                    double valueri, valueii, valuerj, valueij;
+                    chunk_array_get(datar, i, &valueri);
+                    chunk_array_get(datar, j, &valuerj);
+                    chunk_array_get(datai, i, &valueii);
+                    chunk_array_get(datai, j, &valueij);
+
+                    tempr = valueri * wr - valueii * wi;
+                    tempi = valueri * wi + valueii * wr;
+                    chunk_data_save(datar, i, valuerj - tempr);
+                    chunk_data_save(datai, i, valueij - tempi);
+                    chunk_data_save(datar, j, valuerj + tempr);
+                    chunk_data_save(datai, j, valueij + tempi);
                     j += istep;
                 }
                 imin = 2 * imin - i1;
@@ -405,23 +443,27 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
                         j3max = j2 + np2 - ifp2;
                         for (j3 = j2; j3 <= j3max; j3 += ifp2) {
                             j = jmin + ifp2 - ifp1;
-                            sr = datar[j];
-                            si = datai[j];
+                            chunk_array_get(datar, j, &sr);
+                            chunk_array_get(datai, j, &si);
                             oldsr = 0.;
                             oldsi = 0.;
                             j -= ifp1;
                         L620:
                             stmpr = sr;
                             stmpi = si;
-                            sr = twowr * sr - oldsr + datar[j];
-                            si = twowr * si - oldsi + datai[j];
+                            double valuerj, valueij;
+                            chunk_array_get(datar, j, &valuerj);
+                            chunk_array_get(datai, j, &valueij);
+
+                            sr = twowr * sr - oldsr + valuerj;
+                            si = twowr * si - oldsi + valueij;
                             oldsr = stmpr;
                             oldsi = stmpi;
                             j -= ifp1;
                             if (j > jmin)
                                 goto L620;
-                            workr[i] = wr * sr - wi * si - oldsr + datar[j];
-                            worki[i] = wi * sr + wr * si - oldsi + datai[j];
+                            workr[i] = wr * sr - wi * si - oldsr + valuerj;
+                            worki[i] = wi * sr + wr * si - oldsi + valueij;
                             jmin += ifp2;
                             i++;
                         }
@@ -433,8 +475,8 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
                     for (j2 = i3; j2 <= j2max; j2 += ifp1) {
                         j3max = j2 + np2 - ifp2;
                         for (j3 = j2; j3 <= j3max; j3 += ifp2) {
-                            datar[j3] = workr[i];
-                            datai[j3] = worki[i];
+                            chunk_array_save(datar, j3, workr[i]);
+                            chunk_array_save(datai, j3, worki[i]);
                             i++;
                         }
                     }
@@ -478,16 +520,22 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
     L710:
         j = jmin;
         for (i = imin; i <= ntot; i += np2) {
-            sumr = (datar[i] + datar[j]) / 2.;
-            sumi = (datai[i] + datai[j]) / 2.;
-            difr = (datar[i] - datar[j]) / 2.;
-            difi = (datai[i] - datai[j]) / 2.;
+            double valueri, valueii, valuerj, valueij;
+            chunk_array_get(datar, i, &valueri);
+            chunk_array_get(datai, i, &valueii);
+            chunk_array_get(datar, j, &valuerj);
+            chunk_array_get(datai, j, &valueij);
+
+            sumr = (valueri + valuerj) / 2.;
+            sumi = (valueii + valueij) / 2.;
+            difr = (valueri - valuerj) / 2.;
+            difi = (valueii - valueij) / 2.;
             tempr = wr * sumi + wi * difr;
             tempi = wi * sumi - wr * difr;
-            datar[i] = sumr + tempr;
-            datai[i] = difi + tempi;
-            datar[j] = sumr - tempr;
-            datai[j] = tempi - difi;
+            chunk_data_save(datar, i, sumr + tempr);
+            chunk_data_save(datai, i, difi + tempi);
+            chunk_data_save(datar, j, sumr - tempr);
+            chunk_data_save(datai, j, tempi - difi);
             j += np2;
         }
         imin++;
@@ -504,7 +552,9 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
         if (ifrwd == 0)
             goto L740;
         for (i = imin; i <= ntot; i += np2) {
-            datai[i] = -datai[i];
+            double valueii;
+            chunk_data_get(datai, i, &valueii);
+            chunk_data_save(datai, i, -valueii);
         }
     L740:
         np2 *= 2;
@@ -515,36 +565,56 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
         imin = imax - nhalf;
         i = imin;
         goto L755;
-    L750:
-        datar[j] = datar[i];
-        datai[j] = -datai[i];
+    L750: ;
+        double valueri, valueii;
+        chunk_data_get(datar, i, &valueri);
+        chunk_data_get(datai, i, &valueii);
+        chunk_data_save(datar, j, valueri);
+        chunk_data_save(datai, j, -valueii);
     L755:
         i++;
         j--;
         if (i < imax)
             goto L750;
-        datar[j] = datar[imin] - datai[imin];
-        datai[j] = 0.;
+        double valuerimin, valueiimin;
+        chunk_data_get(datar, imin, &valuerimin);
+        chunk_data_get(datai, imin, &valueiimin);
+
+        chunk_data_save(datar, j, valuerimin - valueiimin);
+        chunk_data_save(datai, j, 0.);
         if (i >= j) {
             goto L780;
         } else {
             goto L770;
         }
     L765:
-        datar[j] = datar[i];
-        datai[j] = datai[i];
+        chunk_data_get(datar, i, &valueri);
+        chunk_data_get(datai, i, &valueii);
+
+        chunk_data_save(datar, j, valueri);
+        chunk_data_save(datai, j, valueii);
     L770:
         i--;
         j--;
         if (i > imin)
             goto L765;
-        datar[j] = datar[imin] + datai[imin];
-        datai[j] = 0.;
+
+        chunk_data_get(datar, imin, &valuerimin);
+        chunk_data_get(datai, imin, &valueiimin);
+
+        chunk_data_save(datar, j, valuerimin + valueiimin);
+        chunk_data_save(datai, j, 0.);
+
         imax = imin;
         goto L745;
-    L780:
-        datar[1] += datai[1];
-        datai[1] = 0.;
+    L780: ;
+        double valuei1, valuer1;
+        chunk_data_get(datai, 1, &valuei1);
+        chunk_data_get(datar, 1, &valuer1);
+
+
+        chunk_data_save(datar, 1, valuei1 + valuer1);
+        chunk_data_save(datai, 1, 0.);
         goto L900;
 
         /*complete a real transform for the 2nd, 3rd, ... dimension by conjugate symmetries*/
@@ -562,15 +632,24 @@ void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icp
                 if (idim > 2) {
                     j = jmax + np0;
                     for (i = imin; i <= imax; i++) {
-                        datar[i] = datar[j];
-                        datai[i] = -datai[j];
+                        double valuerj, valueij;
+                        chunk_data_get(datar, j, &valuerj);
+                        chunk_data_get(datai, j, &valueij);
+
+                        chunk_data_save(datar, i, valuerj);
+                        chunk_data_save(datai, i, -valueij);
                         j--;
                     }
                 }
                 j = jmax;
                 for (i = imin; i <= imax; i += np0) {
-                    datar[i] = datar[j];
-                    datai[i] = -datai[j];
+                    double valuerj, valueij;
+                    chunk_data_get(datar, j, &valuerj);
+                    chunk_data_get(datai, j, &valueij);
+
+                    chunk_data_save(datar, i, valuerj);
+                    chunk_data_save(datai, i, -valueij);
+
                     j -= np0;
                 }
             }
diff --git a/fftma_module/gen/lib_src/geostat.h b/fftma_module/gen/lib_src/geostat.h
index 35d8bcc..c7244c9 100644
--- a/fftma_module/gen/lib_src/geostat.h
+++ b/fftma_module/gen/lib_src/geostat.h
@@ -328,7 +328,7 @@ void coordinates(int maille, int i[3], struct grid_mod grid);
 /*variogram: structure defined above     */
 /*grid: structure defined above          */
 /*n: number of gridblocks along X,Y and Z*/
-void covariance(double* covar, struct vario_mod variogram, struct grid_mod grid, int n[3], int cores);
+void covariance(chunk_array_t* covar, struct vario_mod variogram, struct grid_mod grid, int n[3], int cores);
 
 /*computation of the covariance matrix for the well data*/
 /*well coordinates are given as a number of cells       */
@@ -404,7 +404,7 @@ double exponential(double h);
 /*workr: utility real part vector for storage         */
 /*worki: utility imaginary part vector for storage    */
 /*The transformed data are returned in datar and datai*/
-void fourt(double* datar, double* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores);
+void fourt(chunk_array_t* datar, chunk_array_t* datai, int nn[3], int ndim, int ifrwd, int icplx, double* workr, double* worki, int cores);
 
 /*calculates F(x) = (1/a)*exp(-x*x/2)*/
 double funtrun1(double x);
diff --git a/fftma_module/gen/lib_src/prebuild_gwn.c b/fftma_module/gen/lib_src/prebuild_gwn.c
index 02d9028..a6382cf 100644
--- a/fftma_module/gen/lib_src/prebuild_gwn.c
+++ b/fftma_module/gen/lib_src/prebuild_gwn.c
@@ -1,4 +1,5 @@
 #include "geostat.h"
+#include "chunk_array.h"
 #include <math.h>
 #include <stdarg.h>
 #include <stddef.h>
@@ -20,7 +21,7 @@
 /*        must be a Gaussian white noise       */
 /*realization: structure defining a realization*/
 
-void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, double* realization, int solver, int cores, long* seed) {
+void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin, chunk_array_t* realization, int solver, int cores, long* seed) {
     int i, j, k, maille0, maille1;
     int ntot;
 
@@ -32,14 +33,14 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin
         *seed = -(*seed);
 
     ntot = n[0] * n[1] * n[2];
-    realization[0] = 0.;
+    chunk_array_save(realization, 0, 0.);
     /*printf("Antes de llamar a chunkarray read\n");
     chunk_array_read((*realin).vector_2);
     printf("Despues de llamar a chunkarray read\n");*/
     if (solver == 1) {
         for (i = 0; i < ntot; i++) {
             double value = gasdev(seed, &idum2, &iy, iv, cores);
-            realization[i+1] = value;
+            chunk_array_save(realization, i+1, value);
             //chunk_array_get((*realin).vector_2, i, &realization[i + 1]);
         }
     } else {
@@ -51,10 +52,10 @@ void prebuild_gwn(struct grid_mod grid, int n[3], struct realization_mod* realin
                         //maille0 = i - 1 + (j - 1 + (k - 1) * grid.NY) * grid.NX;
                         //printf("Maille0 es %d", maille0);
                         double value = gasdev(seed, &idum2, &iy, iv, cores);
-                        realization[maille1] = value;
+                        chunk_array_save(realization, maille1, value);
                         //chunk_array_get((*realin).vector_2, maille0, &realization[maille1]);
                     } else {
-                        realization[maille1] = 0.;
+                        chunk_array_save(realization, maille1, 0.);
                     }
                 }
             }