diff --git a/C/projects/CARDAMOM_GENERAL/CARDAMOM_NETCDF_DATA_STRUCTURE.c b/C/projects/CARDAMOM_GENERAL/CARDAMOM_NETCDF_DATA_STRUCTURE.c
index 3112cf79..c81792d7 100644
--- a/C/projects/CARDAMOM_GENERAL/CARDAMOM_NETCDF_DATA_STRUCTURE.c
+++ b/C/projects/CARDAMOM_GENERAL/CARDAMOM_NETCDF_DATA_STRUCTURE.c
@@ -114,6 +114,9 @@ TIMESERIES_DRIVER_STRUCT TIME_INDEX;
 TIMESERIES_DRIVER_STRUCT TOTAL_PREC;
 TIMESERIES_DRIVER_STRUCT VPD;
 TIMESERIES_DRIVER_STRUCT YIELD;
+TIMESERIES_DRIVER_STRUCT ET_REF;//For precip-recycling experiments only
+TWO_DIM_VAR_STRUCT ET_REF2;//For precip-recycling experiments only
+
 
 //Summary variables
 
diff --git a/C/projects/CARDAMOM_GENERAL/CARDAMOM_READ_NETCDF_DATA.c b/C/projects/CARDAMOM_GENERAL/CARDAMOM_READ_NETCDF_DATA.c
index e6e93703..7cc30e79 100644
--- a/C/projects/CARDAMOM_GENERAL/CARDAMOM_READ_NETCDF_DATA.c
+++ b/C/projects/CARDAMOM_GENERAL/CARDAMOM_READ_NETCDF_DATA.c
@@ -325,8 +325,46 @@ printf("Done reading all other edc ");
                                 DEFAULT_REFERENCE_MEAN(&DATA->YIELD);
 
         
+
+	DATA->ET_REF.values=ncdf_read_double_var(ncid, "ET_REF", &(DATA->ET_REF.length));
+		DATA->ET_REF.reference_mean=ncdf_read_double_attr(ncid, "ET_REF","reference_mean");
+                                DEFAULT_REFERENCE_MEAN(&DATA->ET_REF);
         
-        
+      printf("Made it here 1\n");
+     DATA->ET_REF2.values=ncdf_read_double_2D(ncid, "ET_REF2", DATA->ET_REF2.dimensions);
+               // DATA->ET_REF2.reference_mean=ncdf_read_double_attr(ncid, "ET_REF2","reference_mean");
+                                //DEFAULT_REFERENCE_MEAN(&DATA->ET_REF2);
+      //printf("Made it here 2\n");
+
+        //           printf("%9.4f\t", DATA->ET_REF.values[0]);
+
+
+      printf("Made it here 3\n");
+                   
+		   printf("%9i\n", DATA->ET_REF2.dimensions[0]);
+		   printf("%9i\n", DATA->ET_REF2.dimensions[1]);
+
+		   printf("%9.4f\n", DATA->ET_REF2.values[0][0]);
+                   printf("%9.4f\n", DATA->ET_REF2.values[0][1]);
+                   printf("%9.4f\n", DATA->ET_REF2.values[0][2]);
+		   printf("%9.4f\n", DATA->ET_REF2.values[1][0]);
+                   printf("%9.4f\n", DATA->ET_REF2.values[1][1]);
+                   printf("%9.4f\n", DATA->ET_REF2.values[1][2]);
+		   printf("%9.4f\n", DATA->ET_REF2.values[2][0]);
+                   printf("%9.4f\n", DATA->ET_REF2.values[2][1]);
+                   printf("%9.4f\n", DATA->ET_REF2.values[2][2]);
+ 
+
+     printf("Made it here 4\n");
+
+
+//   printf("%9.4f\t", DATA->ET_REF2.values[0][2]);
+                   // printf("%9.4f\t", DATA->ET_REF2.values[2][2]);
+
+
+
+
+
         //Summary & derived variables
 // double EDC;
 // double EDCDIAG;
diff --git a/C/projects/CARDAMOM_GENERAL/CARDAMOM_RUN_MODEL.c b/C/projects/CARDAMOM_GENERAL/CARDAMOM_RUN_MODEL.c
index 6564e15e..29745f2e 100644
--- a/C/projects/CARDAMOM_GENERAL/CARDAMOM_RUN_MODEL.c
+++ b/C/projects/CARDAMOM_GENERAL/CARDAMOM_RUN_MODEL.c
@@ -359,6 +359,20 @@ for (n=0;n<N;n++){
 FAILONERROR(nc_get_vara_double(cbrId,cbrParsVarID,(const size_t[]){n,0}, (const size_t[]){1,CARDADATA.nopars},pars));
 
 
+//Aditional step: get ref ET for recycling
+CARDADATA.ncdf_data.ET_REF.values=CARDADATA.ncdf_data.ET_REF2.values[n]; 
+printf("CARDADATA.ncdf_data.ET_REF.values[0]=%2.2f\n",CARDADATA.ncdf_data.ET_REF.values[0]);
+printf("CARDADATA.ncdf_data.ET_REF.values[1]=%2.2f\n",CARDADATA.ncdf_data.ET_REF.values[1]);
+printf("CARDADATA.ncdf_data.ET_REF.values[2]=%2.2f\n",CARDADATA.ncdf_data.ET_REF.values[2]);
+
+printf("CARDADATA.ncdf_data.ET_REF2.values[n][0]=%2.2f\n",CARDADATA.ncdf_data.ET_REF2.values[n][0]);
+printf("CARDADATA.ncdf_data.ET_REF2.values[n][1]=%2.2f\n",CARDADATA.ncdf_data.ET_REF2.values[n][1]);
+printf("CARDADATA.ncdf_data.ET_REF2.values[n][2]=%2.2f\n",CARDADATA.ncdf_data.ET_REF2.values[n][2]);
+
+
+
+printf("Made it here, n = %i\n",n);
+
 /*step 4.2 - running MODEL_LIKELIHOOD_FUNCTION (e.g. DALEC_CDEA_MLF) with parameter vector from file*/
 /*This is stored in CARDADATA.MLF*/
 /*Setting EDCDIAG = 1 to ensure full model run*/
diff --git a/C/projects/CARDAMOM_GENERAL/NETCDF_AUXILLIARY_FUNCTIONS.c b/C/projects/CARDAMOM_GENERAL/NETCDF_AUXILLIARY_FUNCTIONS.c
index 8ca942e6..eb035c5b 100644
--- a/C/projects/CARDAMOM_GENERAL/NETCDF_AUXILLIARY_FUNCTIONS.c
+++ b/C/projects/CARDAMOM_GENERAL/NETCDF_AUXILLIARY_FUNCTIONS.c
@@ -287,6 +287,10 @@ double ** ncdf_read_double_2D(int ncid, const char * varName, size_t * dimLen ){
 	for(int idx=0; idx<dimLen[0]; idx++){
 		//Allocate and populate the child arrays that hold the data
 		resultArray[idx] = (double*)calloc(dimLen[1], sizeof(double));
+
+       // *** THE FIX IS HERE ***
+        // Update the starting index for the first dimension (the row).
+        // On the next loop, it will read from {1,0}, then {2,0}, and so on.
   // THIS IS THE FIX: Update which row you are about to read
         startIndexes[0] = idx;
 
diff --git a/C/projects/CARDAMOM_MODELS/DALEC/DALEC_1100/DALEC_1100.c b/C/projects/CARDAMOM_MODELS/DALEC/DALEC_1100/DALEC_1100.c
index d53ecb5c..1279ea42 100644
--- a/C/projects/CARDAMOM_MODELS/DALEC/DALEC_1100/DALEC_1100.c
+++ b/C/projects/CARDAMOM_MODELS/DALEC/DALEC_1100/DALEC_1100.c
@@ -269,6 +269,8 @@ double *SKT=DATA.ncdf_data.SKT.values;
 double *STRD=DATA.ncdf_data.STRD.values;
 double *DIST=DATA.ncdf_data.DISTURBANCE_FLUX.values;
 double *YIELD=DATA.ncdf_data.YIELD.values;
+double *ET_REF=DATA.ncdf_data.ET_REF.values;
+
 
 /*C-pools, fluxes, meteorology indices*/
 int p=0,f,m,nxp, i;
@@ -456,6 +458,7 @@ int nofluxes=((DALEC *)DATA.MODEL)->nofluxes;
 
 
     /*repeating loop for each timestep*/
+
 for (n=0; n < N_timesteps; n++){
     /*pool index*/
 p=nopools*n;
@@ -671,7 +674,10 @@ FLUXES[f+F.gh_in] =FLUXES[f+F.ground_heat] *DGCM_SEC_DAY;
 FLUXES[f+F.sensible_heat] = Rn - FLUXES[f+F.ground_heat] - FLUXES[f+F.latent_heat];
 
     // Infiltration (mm/day)
-double liquid_in = (PREC[n] - SNOWFALL[n] + FLUXES[f+F.melt]);
+double et_recycled=(FLUXES[f+F.ets]-ET_REF[n])*0.5;
+/*printf ET_rec*/
+//printf("et_recycled = %2.2f\n",et_recycled);
+double liquid_in = (PREC[n] - SNOWFALL[n] + FLUXES[f+F.melt]+et_recycled);
 FLUXES[f+F.infil] = pars[P.max_infil]*(1 - exp(-liquid_in/pars[P.max_infil]));
 
     // Surface runoff (mm/day)
@@ -1135,6 +1141,7 @@ POOLS[nxp+S.D_PSI_LY3]=fmax(HYDROFUN_MOI2PSI(  POOLS[nxp+S.D_SM_LY3],psi_porosit
 
 
 
+
 return 0;
 }
 
@@ -1562,4 +1569,4 @@ DALECmodel->OBSOPE=OBSOPE;
 
 
 
-return 0;}
\ No newline at end of file
+return 0;}
diff --git a/MATLAB/auxi_fun/plotfun_freeze_xyticks.m b/MATLAB/auxi_fun/plotfun_freeze_xyticks.m
index 35b97160..3d73b2a3 100644
--- a/MATLAB/auxi_fun/plotfun_freeze_xyticks.m
+++ b/MATLAB/auxi_fun/plotfun_freeze_xyticks.m
@@ -1,2 +1,3 @@
+function plotfun_freeze_xyticks
 set(gca,'xtick',get(gca,'xtick'),'ytick',get(gca,'ytick'));
-
+end
diff --git a/MATLAB/io_fun/CARDAMOM_WRITE_NC_CBF_FILE.m b/MATLAB/io_fun/CARDAMOM_WRITE_NC_CBF_FILE.m
index 91ae5cb7..4d7b6d2e 100644
--- a/MATLAB/io_fun/CARDAMOM_WRITE_NC_CBF_FILE.m
+++ b/MATLAB/io_fun/CARDAMOM_WRITE_NC_CBF_FILE.m
@@ -19,6 +19,9 @@ function CARDAMOM_WRITE_NC_CBF_FILE(CBF,fname)
 for f=1:numel(fieldnames)
     
 
+s1=size(CBF.(fieldnames{f}).values,1);
+s2=size(CBF.(fieldnames{f}).values,2);
+
        if numel(CBF.(fieldnames{f}).values)==notimesteps
 
        %Enforcing "double" type
@@ -33,6 +36,17 @@ function CARDAMOM_WRITE_NC_CBF_FILE(CBF,fname)
         nccreate(fname,fieldnames{f}); 
        ncwrite(fname,fieldnames{f},CBF.(fieldnames{f}).values);
 
+
+        elseif s1>1 &  s2>1
+
+      %Enforcing "double" type
+       nccreate(fname,fieldnames{f},'Dimensions',{'dim1',s1,'dim2',s2},'FillValue',-9999,'Datatype','double' ); 
+
+       ncwrite(fname,fieldnames{f},CBF.(fieldnames{f}).values);
+
+       
+    disp('Wrote 2D array')
+
            
        
    end