diff --git a/DATA/PROJECT/VINE3D_test/SETTINGS/parameters.ini b/DATA/PROJECT/VINE3D_test/SETTINGS/parameters.ini
index 7e18454d8..9a9d71277 100644
--- a/DATA/PROJECT/VINE3D_test/SETTINGS/parameters.ini
+++ b/DATA/PROJECT/VINE3D_test/SETTINGS/parameters.ini
@@ -37,7 +37,7 @@ relhum_tolerance=102
 water_table_maximum_depth=300
 
 [meteo]
-min_percentage=60
+min_percentage=80
 prec_threshold=0.2
 samani_coefficient=0.17
 thom_threshold=24
diff --git a/agrolib/crop/landUnit.cpp b/agrolib/crop/landUnit.cpp
index be2794188..d39f754d3 100644
--- a/agrolib/crop/landUnit.cpp
+++ b/agrolib/crop/landUnit.cpp
@@ -59,13 +59,3 @@ bool loadLandUnitList(const QSqlDatabase &dbCrop, std::vector<Crit3DLandUnit> &l
 
     return true;
 }
-
-
-int getLandUnitIndex(const std::vector<Crit3DLandUnit> &landUnitList, int idLandUnit)
-{
-    for (int index = 0; index < int(landUnitList.size()); index++)
-        if (landUnitList[index].id == idLandUnit)
-            return index;
-
-    return NODATA;
-}
diff --git a/agrolib/crop/landUnit.h b/agrolib/crop/landUnit.h
index be563b3df..fb11f46c3 100644
--- a/agrolib/crop/landUnit.h
+++ b/agrolib/crop/landUnit.h
@@ -23,7 +23,5 @@
 
     bool loadLandUnitList(const QSqlDatabase &dbCrop, std::vector<Crit3DLandUnit> &landUnitList, QString &errorStr);
 
-    int getLandUnitIndex(const std::vector<Crit3DLandUnit> &landUnitList, int idLandUnit);
-
 
 #endif // LANDUNIT_H
diff --git a/agrolib/grapevine/grapevine.cpp b/agrolib/grapevine/grapevine.cpp
index 22125020c..83e04c5df 100644
--- a/agrolib/grapevine/grapevine.cpp
+++ b/agrolib/grapevine/grapevine.cpp
@@ -21,9 +21,7 @@ const float LAIMIN = 0.01f;
 
 
 Vine3D_Grapevine::Vine3D_Grapevine()
-{
-}
-
+{}
 
 bool Vine3D_Grapevine::compute(bool computeDaily, int secondsPerStep, Crit3DModelCase* modelCase, double chlorophyll)
 {
@@ -106,9 +104,9 @@ bool Vine3D_Grapevine::compute(bool computeDaily, int secondsPerStep, Crit3DMode
     grassTranspiration(modelCase);
 
     return(true);
-
 }
 
+
 void Vine3D_Grapevine::resetLayers()
 {
     for (int i=0 ; i < nrMaxLayers ; i++)
@@ -121,6 +119,7 @@ void Vine3D_Grapevine::resetLayers()
     }
 }
 
+
 void Vine3D_Grapevine::initializeLayers(int myMaxLayers)
 {
     nrMaxLayers = myMaxLayers;
@@ -156,22 +155,25 @@ bool Vine3D_Grapevine::setWeather(double meanDailyTemp, double temp, double irra
     myMeanDailyTemperature = meanDailyTemp;
     double deltaRelHum = MAXVALUE(100.0 - myRelativeHumidity, 0.01);
     myVaporPressureDeficit = 0.01 * deltaRelHum * 613.75 * exp(17.502 * myInstantTemp / (240.97 + myInstantTemp));
-    //globalRadiation = globRad;
-    if ((int(prec) != NODATA) && (int(temp) != NODATA) && (int(windSpeed) != NODATA) && (int(irradiance) != NODATA) && (int(relativeHumidity) != NODATA) && (int(atmosphericPressure) != NODATA)) isReadingOK = true ;
-    return isReadingOK ;
+
+    if ((int(prec) != NODATA) && (int(temp) != NODATA) && (int(windSpeed) != NODATA)
+        && (int(irradiance) != NODATA) && (int(relativeHumidity) != NODATA) && (int(atmosphericPressure) != NODATA))
+        isReadingOK = true;
+
+    return isReadingOK;
 }
 
+
 bool Vine3D_Grapevine::setDerivedVariables(double diffuseIrradiance, double directIrradiance, double cloudIndex, double sunElevation)
 {
     bool isReadingOK = false;
     myDiffuseIrradiance = diffuseIrradiance ;
     myDirectIrradiance = directIrradiance ;
-    //myLongWaveIrradiance = longWaveIrradiance ;
     myCloudiness = MINVALUE(1,MAXVALUE(0,cloudIndex)) ;
-    //myAirVapourPressure = airVapourPressure ;
     mySunElevation =  sunElevation;
-    if (int(sunElevation) != NODATA && int(diffuseIrradiance) != NODATA && int(directIrradiance) != NODATA
-            && int(cloudIndex) != NODATA) isReadingOK = true ;
+    if (int(sunElevation) != NODATA && int(diffuseIrradiance) != NODATA && int(directIrradiance) != NODATA && int(cloudIndex) != NODATA)
+        isReadingOK = true ;
+
     return isReadingOK ;
 }
 
@@ -185,32 +187,45 @@ void Vine3D_Grapevine::initializeWaterStress(Crit3DModelCase* modelCase)
 }
 
 
-bool Vine3D_Grapevine::setSoilProfile(Crit3DModelCase* modelCase, double* myWiltingPoint, double* myFieldCapacity,
-                            double* myPsiSoilProfile, double* mySoilWaterContentProfile,
-                            double* mySoilWaterContentFC, double* mySoilWaterContentWP)
+bool Vine3D_Grapevine::setSoilProfile(Crit3DModelCase* modelCase, std::vector<double>& myWiltingPoint, std::vector<double>& myFieldCapacity,
+                            std::vector<double>& myPsiSoilProfile, std::vector<double>& mySoilWaterContentProfile,
+                            std::vector<double>& mySoilWaterContentFC, std::vector<double>& mySoilWaterContentWP)
 {
-    double psiSoilProfile;
-    double logPsiSoilAverage = 0.;
-    double logPsiFCAverage = 0.;
-    double soilFieldCapacity;
+    // check last Soil layer
+    int lastLayer = modelCase->soilLayersNr - 1;
+    for (int i = 1; i < modelCase->soilLayersNr; i++)
+    {
+        if (isEqual(myPsiSoilProfile[i], NODATA) || isEqual(myFieldCapacity[i], NODATA) || isEqual(mySoilWaterContentProfile[i], NODATA)
+            || isEqual(mySoilWaterContentProfile[i], NODATA) || isEqual(mySoilWaterContentWP[i], NODATA))
+        {
+            lastLayer = i-1;
+            break;
+        }
+    }
+    if (lastLayer == 0)
+        return false;
 
     psiSoilAverage = 0.;
     psiFieldCapacityAverage = 0.;
 
-    if (int(myWiltingPoint[int(modelCase->soilLayersNr / 2)]) == NODATA)
-        return false;
+    // initialize vector
+    for (int i = 1; i < modelCase->soilLayersNr; i++)
+    {
+        transpirationLayer[i] = 0.;
+        transpirationCumulatedGrass[i] = 0.;
+        fractionTranspirableSoilWaterProfile[i] = 0.;
+    }
 
-    wiltingPoint = myWiltingPoint[int(modelCase->soilLayersNr / 2)] / 101.97;     // conversion from mH2O to MPa
+    int midLayer = (lastLayer + 1) / 2;
+    wiltingPoint = myWiltingPoint[midLayer] / 101.97;                   // conversion from mH2O to MPa
 
     //layer 0: surface, no soil
-    for (int i = 1; i < modelCase->soilLayersNr; i++)
+    double logPsiSoilAverage = 0.;
+    double logPsiFCAverage = 0.;
+    for (int i = 1; i <= lastLayer; i++)
     {
-        if (isEqual(myPsiSoilProfile[i], NODATA) || isEqual(myFieldCapacity[i], NODATA) || isEqual(mySoilWaterContentProfile[i], NODATA)
-                || isEqual(mySoilWaterContentProfile[i], NODATA) || isEqual(mySoilWaterContentWP[i], NODATA))
-            return false;
-
-        soilFieldCapacity = myFieldCapacity[i]/101.97; // conversion from mH2O to MPa
-        psiSoilProfile = MINVALUE(myPsiSoilProfile[i],-1.)/101.97 ; // conversion from mH2O to MPa
+        double soilFieldCapacity = myFieldCapacity[i] / 101.97;                 // conversion from mH2O to MPa
+        double psiSoilProfile = MINVALUE(myPsiSoilProfile[i], -1.) / 101.97 ;   // conversion from mH2O to MPa
         logPsiSoilAverage += log(-psiSoilProfile) * modelCase->rootDensity[i];
         logPsiFCAverage += log(-soilFieldCapacity) * modelCase->rootDensity[i];
     }
@@ -219,28 +234,26 @@ bool Vine3D_Grapevine::setSoilProfile(Crit3DModelCase* modelCase, double* myWilt
     psiFieldCapacityAverage = -exp(logPsiFCAverage);
     fractionTranspirableSoilWaterAverage = 0;
 
-    double waterContent, waterContentFC, waterContentWP;
-
-    for (int i = 0; i < modelCase->soilLayersNr; i++)
+    for (int i = 0; i <= lastLayer; i++)
     {
-        waterContent = mySoilWaterContentProfile[i];
-        waterContentFC = mySoilWaterContentFC[i];
-        waterContentWP = mySoilWaterContentWP[i];
+        double waterContent = mySoilWaterContentProfile[i];
+        double waterContentFC = mySoilWaterContentFC[i];
+        double waterContentWP = mySoilWaterContentWP[i];
 
         fractionTranspirableSoilWaterProfile[i] = MAXVALUE(0, MINVALUE(1, (waterContent - waterContentWP) / (waterContentFC - waterContentWP)));
         fractionTranspirableSoilWaterAverage += fractionTranspirableSoilWaterProfile[i] * modelCase->rootDensity[i];
-        transpirationLayer[i] = 0.;
-        transpirationCumulatedGrass[i] = 0. ;
     }
+
     return true ;
 }
 
-bool Vine3D_Grapevine::setStatePlant(TstatePlant myStatePlant, bool isVineyard)
+
+bool Vine3D_Grapevine::setStatePlant(TstatePlant myStatePlant, bool isVineyard_)
 {
     statePlant = myStatePlant;
     this->statePlant.outputPlant.transpirationNoStress = 0.;
 
-    if (! isVineyard)
+    if (! isVineyard_)
     {
         statePlant.outputPlant.brixBerry = NODATA;
         statePlant.statePheno.stage = NODATA;
@@ -248,17 +261,8 @@ bool Vine3D_Grapevine::setStatePlant(TstatePlant myStatePlant, bool isVineyard)
         statePlant.stateGrowth.leafAreaIndex = NODATA;
         statePlant.stateGrowth.fruitBiomass = NODATA;
     }
-    return true;
-}
-
-TstatePlant Vine3D_Grapevine::getStatePlant()
-{
-    return(this->statePlant);
-}
 
-ToutputPlant Vine3D_Grapevine::getOutputPlant()
-{
-    return (this->statePlant.outputPlant);
+    return true;
 }
 
 
@@ -268,13 +272,13 @@ void Vine3D_Grapevine::getFixSimulationParameters()
     parameterBindiMigliettaFix.a = -0.28;
     parameterBindiMigliettaFix.b = 0.04;
     parameterBindiMigliettaFix.c = -0.015;
-    //parameterBindiMigliettaFix.baseTemperature = 10;
-    //parameterBindiMigliettaFix.tempMaxThreshold = 35;
     parameterBindiMigliettaFix.extinctionCoefficient = 0.5;
     parameterBindiMigliettaFix.shadedSurface = 0.8;
+
     // Wang Leuning parameters
     parameterWangLeuningFix.stomatalConductanceMin = 0.008;
     parameterWangLeuningFix.optimalTemperatureForPhotosynthesis = 298.15;
+
     // fenovitis parameters
     parameterPhenoVitisFix.a= 0.005;
     parameterPhenoVitisFix.optimalChillingTemp = 2.8;
@@ -525,6 +529,8 @@ void Vine3D_Grapevine::radiationAbsorption()
     {
          sunlit.leafAreaIndex =  0.0 ;
          sunlit.absorbedPAR = 0.0 ;
+
+         // TODO: non servono?
          sunlitAbsorbedNIR = 0.0 ;
          sunlitAbsorbedLW = 0.0 ;
          sunlit.isothermalNetRadiation =  0.0 ;
@@ -537,12 +543,13 @@ void Vine3D_Grapevine::radiationAbsorption()
          shadedAbsorbedLW= dum[16] * (UPSCALINGFUNC(diffuseLightK,statePlant.stateGrowth.leafAreaIndex) - UPSCALINGFUNC(directLightK+diffuseLightK,statePlant.stateGrowth.leafAreaIndex)) ;
          shaded.isothermalNetRadiation = shaded.absorbedPAR + shadedAbsorbedNIR + shadedAbsorbedLW ;
     }
+
      // Convert absorbed PAR into units of mol m-2 s-1
      sunlit.absorbedPAR *= 4.57E-6 ;
      shaded.absorbedPAR *= 4.57E-6 ;
-
 }
 
+
 void Vine3D_Grapevine::leafTemperature()
 {
     // taken from Hydrall Model, Magnani UNIBO
@@ -695,7 +702,8 @@ void Vine3D_Grapevine::aerodynamicalCoupling()
             {
                 sunlitDeltaTemp = 0.0;
             }
-            sunlitDeltaTemp = 0.0;
+            sunlitDeltaTemp = 0.0;  // TODO: check
+
             sunlit.leafTemperature = myInstantTemp + sunlitDeltaTemp	+ ZEROCELSIUS ; //sunlit big-leaf
             stomatalConductanceWater = 10.0/shaded.leafAreaIndex ; //dummy stom res for shaded big-leaf
             //if (shaded.isothermalNetRadiation > 100) stomatalConductanceWater *= pow(100/shaded.isothermalNetRadiation,0.5);
@@ -802,10 +810,9 @@ void Vine3D_Grapevine::upscale(TVineCultivar *cultivar)
         sunlit.darkRespiration = 0.0;
         sunlit.maximalCarboxylationRate = 0.0;
     }
-
-
 }
 
+
 void Vine3D_Grapevine::photosynthesisKernel(TVineCultivar* cultivar, double COMP,double GAC,double GHR,double GSCD,double J,double KC,double KO
                      ,double RD,double RNI,double STOMWL,double VCmax,double *ASS,double *GSC,double *TR)
 {
diff --git a/agrolib/grapevine/grapevine.h b/agrolib/grapevine/grapevine.h
index a0d18faa1..1ff7442c9 100644
--- a/agrolib/grapevine/grapevine.h
+++ b/agrolib/grapevine/grapevine.h
@@ -41,9 +41,9 @@
     enum phenoStage {endoDormancy, ecoDormancy , budBurst , flowering , fruitSet, veraison, physiologicalMaturity, vineSenescence};
     enum TfieldOperation {irrigationOperation, grassSowing, grassRemoving, trimming, leafRemoval,
                           clusterThinning, harvesting, tartaricAnalysis};
-    enum Crit3DLanduse {landuse_nodata, landuse_bare, landuse_vineyard};
+    enum GrapevineLanduse {landuse_nodata, landuse_bare, landuse_vineyard};
 
-    const std::map<std::string, Crit3DLanduse> landuseNames = {
+    const std::map<std::string, GrapevineLanduse> landuseNames = {
         { "UNDEFINED", landuse_nodata },
         { "BARESOIL", landuse_bare },
         { "VINEYARD", landuse_vineyard}
@@ -191,7 +191,7 @@
 
     struct Crit3DModelCase {
         int id;
-        Crit3DLanduse landuse;
+        GrapevineLanduse landuse;
         int soilIndex;
 
         float shootsPerPlant;
@@ -300,8 +300,7 @@
         double psiSoilAverage;
         double psiFieldCapacityAverage;
 
-        //double* layerRootDensity;
-        double totalStomatalConductance, totalStomatalConductanceNoStress ;
+        double totalStomatalConductance, totalStomatalConductanceNoStress;
         double transpirationInstant;
         double* currentProfile;
         double* transpirationInstantLayer;          //molH2O m^-2 s^-1
@@ -327,8 +326,6 @@
         bool isAmphystomatic ;
         double specificLeafArea ;
         double alphaLeuning ;
-        //double leafNitrogen ;
-        //double entropicFactorCarboxyliation,entropicFactorElectronTransporRate ;
 
         Vine3D_SunShade shaded ;
         Vine3D_SunShade sunlit ;
@@ -374,7 +371,7 @@
 
         double getWaterStressByPsiSoil(double myPsiSoil,double psiSoilStressParameter,double exponentialFactorForPsiRatio);
         double getWaterStressSawFunction(int index, TVineCultivar *cultivar);
-        //bool getExtractedWaterFromGrassTranspirationandEvaporation(double* myWaterExtractionProfile);
+
         double getWaterStressSawFunctionAverage(TVineCultivar* cultivar);
         double getGrassTranspiration(double stress, double laiGrassMax, double sensitivityToVPD, double fieldCoverByPlant);
         double getFallowTranspiration(double stress, double laiGrassMax, double sensitivityToVPD);
@@ -392,6 +389,11 @@
     public:
         Vine3D_Grapevine();
 
+        TstatePlant getStatePlant(){ return statePlant; }
+        ToutputPlant getOutputPlant() { return statePlant.outputPlant; }
+
+        bool setStatePlant(TstatePlant myStatePlant, bool isVineyard_);
+
         //void initializeGrapevineModel(TVineCultivar* cultivar, double secondsPerStep);
         void initializeLayers(int myMaxLayers);
         bool initializeStatePlant(int doy, Crit3DModelCase *vineField);
@@ -409,15 +411,12 @@
                 double prec , double relativeHumidity , double windSpeed, double atmosphericPressure);
         bool setDerivedVariables (double diffuseIrradiance, double directIrradiance,
                 double cloudIndex, double sunElevation);
-        bool setSoilProfile(Crit3DModelCase *modelCase, double* myWiltingPoint, double *myFieldCapacity,
-                            double *myPsiSoilProfile , double *mySoilWaterContentProfile,
-                            double* mySoilWaterContentFC, double* mySoilWaterContentWP);
-        bool setStatePlant(TstatePlant myStatePlant, bool isVineyard);
 
-        TstatePlant getStatePlant();
-        ToutputPlant getOutputPlant();
+        bool setSoilProfile(Crit3DModelCase* modelCase, std::vector<double>& myWiltingPoint, std::vector<double>& myFieldCapacity,
+                            std::vector<double>& myPsiSoilProfile, std::vector<double>& mySoilWaterContentProfile,
+                            std::vector<double>& mySoilWaterContentFC, std::vector<double>& mySoilWaterContentWP);
+
         double* getExtractedWater(Crit3DModelCase* modelCase);
-        //bool getOutputPlant(int hour, ToutputPlant *outputPlant);
         double getStressCoefficient();
         double getRealTranspirationGrapevine(Crit3DModelCase *modelCase);
         double getRealTranspirationGrass(Crit3DModelCase *modelCase);
diff --git a/agrolib/project/project.cpp b/agrolib/project/project.cpp
index dc86f2746..73410dd64 100644
--- a/agrolib/project/project.cpp
+++ b/agrolib/project/project.cpp
@@ -934,8 +934,13 @@ Crit3DTime Project::getCrit3DCurrentTime()
 QDateTime Project::getCurrentTime()
 {
     QDateTime myDateTime;
-    myDateTime.setDate(this->currentDate);
-    return myDateTime.addSecs(this->currentHour * HOUR_SECONDS);
+    if (gisSettings.isUTC)
+    {
+        myDateTime.setTimeSpec(Qt::UTC);
+    }
+
+    myDateTime.setDate(currentDate);
+    return myDateTime.addSecs(currentHour * HOUR_SECONDS);
 }
 
 
diff --git a/bin/CRITERIA3D/criteria3DProject.cpp b/bin/CRITERIA3D/criteria3DProject.cpp
index 9f6752114..6e58092ed 100644
--- a/bin/CRITERIA3D/criteria3DProject.cpp
+++ b/bin/CRITERIA3D/criteria3DProject.cpp
@@ -422,13 +422,14 @@ bool Crit3DProject::runModels(QDateTime firstTime, QDateTime lastTime, bool isRe
         for (int hour = firstHour; hour <= lastHour; hour++)
         {
             setCurrentHour(hour);
+            if (currentSeconds == 3600)
+                isRestart = false;
 
-            if (! runModelHour(getCurrentTime(), currentOutputPath, isRestart))
+            if (! runModelHour(currentOutputPath, isRestart))
             {
                 logError();
                 return false;
             }
-            isRestart = false;
 
             // output points
             if (isSaveOutputPoints() && currentSeconds == 3600)
@@ -933,10 +934,13 @@ bool Crit3DProject::updateDailyTemperatures()
 }
 
 
-bool Crit3DProject::runModelHour(const QDateTime &myDateTime, const QString& hourlyOutputPath, bool isRestart)
+bool Crit3DProject::runModelHour(const QString& hourlyOutputPath, bool isRestart)
 {
     if (! isRestart)
     {
+        QDateTime myDateTime = getCurrentTime();
+        currentSeconds = 0;
+
         hourlyMeteoMaps->setComputed(false);
         radiationMaps->setComputed(false);
 
@@ -955,14 +959,14 @@ bool Crit3DProject::runModelHour(const QDateTime &myDateTime, const QString& hou
                 return false;
 
             hourlyMeteoMaps->setComputed(true);
-            emit updateOutputSignal();
+            qApp->processEvents();
         }
 
         if (processes.computeRadiation)
         {
             if (! interpolateAndSaveHourlyMeteo(globalIrradiance, myDateTime, hourlyOutputPath, isSaveOutputRaster()))
                 return false;
-            emit updateOutputSignal();
+            qApp->processEvents();
         }
 
         if (processes.computeSnow)
@@ -973,7 +977,7 @@ bool Crit3DProject::runModelHour(const QDateTime &myDateTime, const QString& hou
             {
                 return false;
             }
-            emit updateOutputSignal();
+            qApp->processEvents();
         }
 
         // initalize sink / source
@@ -1002,6 +1006,7 @@ bool Crit3DProject::runModelHour(const QDateTime &myDateTime, const QString& hou
         if (processes.computeCrop)
         {
             updateDailyTemperatures();
+            qApp->processEvents();
         }
 
         if (processes.computeWater)
@@ -1011,6 +1016,8 @@ bool Crit3DProject::runModelHour(const QDateTime &myDateTime, const QString& hou
             if (! setSinkSource())
                 return false;
         }
+
+        emit updateOutputSignal();
     }
 
     // soil fluxes
@@ -1018,12 +1025,10 @@ bool Crit3DProject::runModelHour(const QDateTime &myDateTime, const QString& hou
     {
         if (! isRestart)
         {
-            logInfo("\nCompute soil fluxes: " + myDateTime.toString());
+            logInfo("\nCompute soil fluxes: " + getCurrentTime().toString());
         }
 
-        runModel(3600, isRestart);
-
-        qApp->processEvents();
+        runWaterFluxes3DModel(3600, isRestart);
     }
 
     // soil heat
diff --git a/bin/CRITERIA3D/criteria3DProject.h b/bin/CRITERIA3D/criteria3DProject.h
index 4eae14520..a2fc43cc9 100644
--- a/bin/CRITERIA3D/criteria3DProject.h
+++ b/bin/CRITERIA3D/criteria3DProject.h
@@ -90,7 +90,7 @@
         bool computeSnowModel();
         void computeSnowPoint(int row, int col);
 
-        bool runModelHour(const QDateTime &myDateTime, const QString& hourlyOutputPath, bool isRestart = false);
+        bool runModelHour(const QString& hourlyOutputPath, bool isRestart = false);
 
         void setAllHourlyMeteoMapsComputed(bool value);
 
diff --git a/bin/CRITERIA3D/mainwindow.cpp b/bin/CRITERIA3D/mainwindow.cpp
index 5d18a6462..139bda912 100644
--- a/bin/CRITERIA3D/mainwindow.cpp
+++ b/bin/CRITERIA3D/mainwindow.cpp
@@ -433,7 +433,7 @@ bool MainWindow::contextMenuRequested(QPoint localPos)
             int id = myProject.getLandUnitIdGeo(geoPos.latitude(), geoPos.longitude());
             if (id != NODATA)
             {
-                int index = getLandUnitIndex(myProject.landUnitList, id);
+                int index = myProject.getLandUnitListIndex(id);
                 if (index != NODATA)
                 {
                     Crit3DLandUnit landUnit = myProject.landUnitList[index];
@@ -461,7 +461,7 @@ bool MainWindow::contextMenuRequested(QPoint localPos)
             int id = myProject.getLandUnitIdGeo(geoPos.latitude(), geoPos.longitude());
             if (id != NODATA)
             {
-                int index = getLandUnitIndex(myProject.landUnitList, id);
+                int index = myProject.getLandUnitListIndex(id);
                 if (index != NODATA)
                 {
                     QString infoStr;
@@ -762,6 +762,7 @@ void MainWindow::updateModelTime()
         }
         minutes = 0;
     }
+
     QDateTime currentDateTime;
     currentDateTime.setDate(date);
     currentDateTime.setTime(QTime(hour, minutes, seconds));
diff --git a/bin/CRITERIA3D/shared/project3D.cpp b/bin/CRITERIA3D/shared/project3D.cpp
index 27af3ce6a..3f20230a5 100644
--- a/bin/CRITERIA3D/shared/project3D.cpp
+++ b/bin/CRITERIA3D/shared/project3D.cpp
@@ -404,12 +404,15 @@ bool Project3D::loadLandUseMap(const QString &fileName)
     landUseMapFileName = getCompleteFileName(fileName, PATH_GEO);
 
     std::string errorStr;
-    if (! gis::openRaster(landUseMapFileName.toStdString(), &landUseMap, gisSettings.utmZone, errorStr))
+    gis::Crit3DRasterGrid raster;
+    if (! gis::openRaster(landUseMapFileName.toStdString(), &raster, gisSettings.utmZone, errorStr))
     {
         logError("Load land use map failed: " + landUseMapFileName + "\n" + QString::fromStdString(errorStr));
         return false;
     }
 
+    gis::resampleGrid(raster, &landUseMap, DEM.header, aggrPrevailing, 0);
+
     logInfo("Land use map = " + landUseMapFileName);
     return true;
 }
@@ -426,13 +429,14 @@ bool Project3D::loadSoilMap(const QString &fileName)
     soilMapFileName = getCompleteFileName(fileName, PATH_GEO);
 
     std::string errorStr;
-    if (! gis::openRaster(soilMapFileName.toStdString(), &soilMap, gisSettings.utmZone, errorStr))
+    gis::Crit3DRasterGrid raster;
+    if (! gis::openRaster(soilMapFileName.toStdString(), &raster, gisSettings.utmZone, errorStr))
     {
         logError("Loading soil map failed: " + soilMapFileName + "\n" + QString::fromStdString(errorStr));
         return false;
     }
 
-    gis::updateMinMaxRasterGrid(&(soilMap));
+    gis::resampleGrid(raster, &soilMap, DEM.header, aggrPrevailing, 0);
 
     logInfo("Soil map = " + soilMapFileName);
 
@@ -529,6 +533,11 @@ void Project3D::setIndexMaps()
                         indexMap.at(layer).value[row][col] = currentIndex;
                         currentIndex++;
                     }
+                    else
+                    {
+                        // test
+                        indexMap.at(layer).value[row][col] = indexMap.at(layer).header->flag;
+                    }
                 }
             }
         }
@@ -930,14 +939,14 @@ bool Project3D::initializeSoilMoisture(int month)
 }
 
 
-/*! \brief runModel
+/*! \brief runWaterFluxes3DModel
  *  \param totalTimeStep [s]
  */
-void Project3D::runModel(double totalTimeStep, bool isRestart)
+void Project3D::runWaterFluxes3DModel(double totalTimeStep, bool isRestart)
 {
     if (! isRestart)
     {
-        currentSeconds = 0;                 // [s]
+        currentSeconds = 0;                                 // [s]
         soilFluxes3D::initializeBalance();
 
         previousTotalWaterContent = soilFluxes3D::getTotalWaterContent();       // [m3]
@@ -955,7 +964,7 @@ void Project3D::runModel(double totalTimeStep, bool isRestart)
     {
         currentSeconds += soilFluxes3D::computeStep(totalTimeStep - currentSeconds);
 
-        if (modelPause)
+        if (modelPause && currentSeconds < totalTimeStep)
         {
             emit updateOutputSignal();
             return;
@@ -1039,7 +1048,7 @@ bool Project3D::loadCropDatabase(const QString &fileName)
 }
 
 
-int Project3D::getLandUnitIdUTM(double x, double y)
+int Project3D::getLandUnitFromUtm(double x, double y)
 {
     if (! landUseMap.isLoaded)
         return NODATA;
@@ -1062,7 +1071,7 @@ int Project3D::getLandUnitIdGeo(double lat, double lon)
     double x, y;
     gis::latLonToUtmForceZone(gisSettings.utmZone, lat, lon, &x, &y);
 
-    return getLandUnitIdUTM(x, y);
+    return getLandUnitFromUtm(x, y);
 }
 
 
@@ -1075,15 +1084,25 @@ int Project3D::getLandUnitIndexRowCol(int row, int col)
     }
 
     double x, y;
-    DEM.getXY(row, col, x, y);
+    gis::getUtmXYFromRowCol(DEM.header, row, col, &x, &y);
 
-    int id = getLandUnitIdUTM(x, y);
+    int id = getLandUnitFromUtm(x, y);
     if (id == NODATA)
-    {
         return NODATA;
+
+    return getLandUnitListIndex(id);
+}
+
+
+int Project3D::getLandUnitListIndex(int id)
+{
+    for (int index = 0; index < int(landUnitList.size()); index++)
+    {
+        if (landUnitList[index].id == id)
+            return index;
     }
 
-    return getLandUnitIndex(landUnitList, id);
+    return NODATA;
 }
 
 
diff --git a/bin/CRITERIA3D/shared/project3D.h b/bin/CRITERIA3D/shared/project3D.h
index 87c869201..ae7be8a50 100644
--- a/bin/CRITERIA3D/shared/project3D.h
+++ b/bin/CRITERIA3D/shared/project3D.h
@@ -158,7 +158,7 @@
         double getSoilLayerTop(unsigned int i);
         double getSoilLayerBottom(unsigned int i);
         int getSoilLayerIndex(double depth);
-        int getLandUnitIdUTM(double x, double y);
+        int getLandUnitFromUtm(double x, double y);
         int getLandUnitIdGeo(double lat, double lon);
         int getLandUnitIndexRowCol(int row, int col);
 
@@ -168,6 +168,8 @@
         int getSoilListIndex(double x, double y);
         QString getSoilCode(double x, double y);
 
+        int getLandUnitListIndex(int id);
+
         int getSoilIndex(long row, long col);
         bool isWithinSoil(int soilIndex, double depth);
 
@@ -184,7 +186,7 @@
         double assignTranspiration(int row, int col, double currentLai, double currentDegreeDays);
 
         bool setSinkSource();
-        void runModel(double totalTimeStep, bool isRestart = false);
+        void runWaterFluxes3DModel(double totalTimeStep, bool isRestart = false);
         bool updateCrop(QDateTime myTime);
 
         bool computeCriteria3DMap(gis::Crit3DRasterGrid &outputRaster, criteria3DVariable var, int layerIndex);
diff --git a/bin/VINE3D/modelCore.cpp b/bin/VINE3D/modelCore.cpp
index e6bf9b9c9..9036ccfe2 100644
--- a/bin/VINE3D/modelCore.cpp
+++ b/bin/VINE3D/modelCore.cpp
@@ -3,6 +3,7 @@
 #include <qstring.h>
 #include <QDate>
 
+#include "commonConstants.h"
 #include "vine3DProject.h"
 #include "waterBalance.h"
 #include "crit3dDate.h"
@@ -20,31 +21,21 @@ extern Vine3DProject myProject;
 
 bool setSoilProfileCrop(Vine3DProject* myProject, int row, int col, Crit3DModelCase* modelCase)
 {
-    double* soilWPProfile = getSoilVarProfile(myProject, row, col, soil::soilWaterPotentialWP);
-    double* soilFCProfile = getSoilVarProfile(myProject, row, col, soil::soilWaterPotentialFC) ;
-    double* matricPotentialProfile = getCriteria3DVarProfile(myProject, row, col, waterMatricPotential);
-    double* waterContentProfile = getCriteria3DVarProfile(myProject, row, col, volumetricWaterContent);
-    double* waterContentProfileWP = getSoilVarProfile(myProject, row, col, soil::soilWaterContentWP);
-    double* waterContentProfileFC = getSoilVarProfile(myProject, row, col, soil::soilWaterContentFC);
-
-    if (! myProject->grapevine.setSoilProfile(modelCase, soilWPProfile, soilFCProfile,
-            matricPotentialProfile, waterContentProfile, waterContentProfileFC,
-            waterContentProfileWP)) return false;
-
-    free(soilWPProfile);
-    free(soilFCProfile);
-    free(matricPotentialProfile);
-    free(waterContentProfile);
-    free(waterContentProfileWP);
-    free(waterContentProfileFC);
+    std::vector<double> soilWPProfile = getSoilVarProfile(myProject, row, col, soil::soilWaterPotentialWP);
+    std::vector<double> soilFCProfile = getSoilVarProfile(myProject, row, col, soil::soilWaterPotentialFC) ;
+    std::vector<double> waterContentProfileWP = getSoilVarProfile(myProject, row, col, soil::soilWaterContentWP);
+    std::vector<double> waterContentProfileFC = getSoilVarProfile(myProject, row, col, soil::soilWaterContentFC);
 
-    return true;
+    std::vector<double> matricPotentialProfile = getCriteria3DVarProfile(myProject, row, col, waterMatricPotential);
+    std::vector<double> waterContentProfile = getCriteria3DVarProfile(myProject, row, col, volumetricWaterContent);
+
+    return myProject->grapevine.setSoilProfile(modelCase, soilWPProfile, soilFCProfile,
+                                               matricPotentialProfile, waterContentProfile, waterContentProfileFC, waterContentProfileWP);
 }
 
 
 bool assignIrrigation(Vine3DProject* myProject, Crit3DTime myTime)
 {
-    int fieldIndex;
     float nrHours;
     float irrigationRate, rate;
     int hour = myTime.getHour();
@@ -61,10 +52,11 @@ bool assignIrrigation(Vine3DProject* myProject, Crit3DTime myTime)
                 //initialize
                 myProject->vine3DMapsH->mapHourlyIrrigation->value[row][col] = 0.0;
 
-                fieldIndex = myProject->getModelCaseIndex(unsigned(row), unsigned(col));
-                if (fieldIndex > 0)
+                int caseIndex = myProject->getModelCaseIndex(row, col);
+                if (caseIndex != NODATA)
                 {
                     idBook = 0;
+                    int fieldIndex = myProject->modelCases[caseIndex].id;
                     while (myProject->getFieldBookIndex(idBook, myDate, fieldIndex, &idBook))
                     {
                         if (myProject->fieldBook[idBook].operation == irrigationOperation)
@@ -152,9 +144,9 @@ bool modelDailyCycle(bool isInitialState, Crit3DDate myDate, int nrHours,
         {
             for (long col = 0; col < myProject->DEM.header->nrCols; col++)
             {
-                if (int(myProject->DEM.value[row][col]) != int(myProject->DEM.header->flag))
+                modelCaseIndex = myProject->getModelCaseIndex(row, col);
+                if (modelCaseIndex != NODATA)
                 {
-                    modelCaseIndex = myProject->getModelCaseIndex(unsigned(row), unsigned(col));
                     isNewModelCase = (int(myProject->statePlantMaps->fruitBiomassMap->value[row][col])
                                         == int(myProject->statePlantMaps->fruitBiomassMap->header->flag));
 
@@ -165,12 +157,13 @@ bool modelDailyCycle(bool isInitialState, Crit3DDate myDate, int nrHours,
                                 double(myProject->hourlyMeteoMaps->mapHourlyPrec->value[row][col]),
                                 double(myProject->hourlyMeteoMaps->mapHourlyRelHum->value[row][col]),
                                 double(myProject->hourlyMeteoMaps->mapHourlyWindScalarInt->value[row][col]),
-                                PRESS)) {
+                                PRESS))
+                    {
                         myProject->errorString = grapevineError(myCurrentTime, row, col, "Weather data missing");
                         return(false);
                     }
 
-                    if (!myProject->grapevine.setDerivedVariables(
+                    if (! myProject->grapevine.setDerivedVariables(
                                 double(myProject->radiationMaps->diffuseRadiationMap->value[row][col]),
                                 double(myProject->radiationMaps->beamRadiationMap->value[row][col]),
                                 double(myProject->radiationMaps->transmissivityMap->value[row][col] / CLEAR_SKY_TRANSMISSIVITY_DEFAULT),
@@ -182,21 +175,23 @@ bool modelDailyCycle(bool isInitialState, Crit3DDate myDate, int nrHours,
 
                     myProject->grapevine.resetLayers();
 
-                    if (! setSoilProfileCrop(myProject, row, col, &(myProject->modelCases[modelCaseIndex]))) {
-                        myProject->errorString = grapevineError(myCurrentTime, row, col, "Error setting soil profile");
+                    if (! setSoilProfileCrop(myProject, row, col, &(myProject->modelCases[modelCaseIndex])))
+                    {
+                        myProject->errorString = grapevineError(myCurrentTime, row, col, "Error in soil profile setting");
                         return false;
                     }
 
                     if ((isInitialState) || (isNewModelCase))
                     {
-                        if(!myProject->grapevine.initializeStatePlant(getDoyFromDate(myDate), &(myProject->modelCases[modelCaseIndex])))
+                        if(! myProject->grapevine.initializeStatePlant(getDoyFromDate(myDate), &(myProject->modelCases[modelCaseIndex])))
                         {
                             myProject->logInfo("Could not initialize grapevine in the present growing season.\nIt will be replaced by a complete grass cover.");
                         }
                     }
                     else
                     {
-                        if (!setStatePlantfromMap(row, col, myProject)) return(false);
+                        if (! setStatePlantfromMap(row, col, myProject))
+                            return false;
                         myProject->grapevine.setStatePlant(myProject->statePlant, true);
                     }
                     double chlorophyll = NODATA;
@@ -204,7 +199,7 @@ bool modelDailyCycle(bool isInitialState, Crit3DDate myDate, int nrHours,
                     if (! myProject->grapevine.compute((myCurrentTime == myFirstTime), myTimeStep, &(myProject->modelCases[modelCaseIndex]), chlorophyll))
                     {
                         myProject->errorString = grapevineError(myCurrentTime, row, col, "Error in grapevine computation");
-                        return(false);
+                        return false;
                     }
 
                     // check field book (first hour)
@@ -255,7 +250,6 @@ bool modelDailyCycle(bool isInitialState, Crit3DDate myDate, int nrHours,
         // Irrigation
         assignIrrigation(myProject, myCurrentTime);
 
-
         if (! myProject->computeVine3DWaterSinkSource())
             return false;
 
diff --git a/bin/VINE3D/vine3DProject.cpp b/bin/VINE3D/vine3DProject.cpp
index 396c56dce..0b836299c 100644
--- a/bin/VINE3D/vine3DProject.cpp
+++ b/bin/VINE3D/vine3DProject.cpp
@@ -458,6 +458,7 @@ bool Vine3DProject::setModelCasesMap()
     // set model cases
     int nrModelCases = nrInputCases * nrSoils;
     modelCases.resize(nrModelCases);
+    landUnitList.resize(nrModelCases);
 
     for (int i = 0; i < nrInputCases; i++)
     {
@@ -466,6 +467,7 @@ bool Vine3DProject::setModelCasesMap()
             int index = nrSoils * i + j;
             modelCases[index] = inputModelCases[i];
             modelCases[index].soilIndex = soilIndexList[j];
+            landUnitList[index].id = index;
         }
     }
 
@@ -499,7 +501,7 @@ bool Vine3DProject::setModelCasesMap()
 }
 
 
-bool Vine3DProject::readFieldQuery(QSqlQuery &myQuery, int &idField, Crit3DLanduse &landuse, int &vineIndex, int &trainingIndex,
+bool Vine3DProject::readFieldQuery(QSqlQuery &myQuery, int &idField, GrapevineLanduse &landuse, int &vineIndex, int &trainingIndex,
                                    float &maxLaiGrass, float &maxIrrigationRate)
 {
     idField = myQuery.value("id_field").toInt();
@@ -555,7 +557,7 @@ bool Vine3DProject::loadFieldsProperties()
 
     int idField, vineIndex, trainingIndex;
     float maxLaiGrass, maxIrrigationRate;
-    Crit3DLanduse landuse;
+    GrapevineLanduse landuse;
 
     QSqlQuery myQuery(dbVine3D);
 
@@ -1352,24 +1354,24 @@ bool Vine3DProject::saveStateAndOutput(QDate myDate)
 }
 
 
-int Vine3DProject::getModelCaseIndex(unsigned row, unsigned col)
+int Vine3DProject::getModelCaseIndex(int row, int col)
 {
-    if (gis::isOutOfGridRowCol(int(row), int(col), landUseMap)) return NODATA;
+    if (gis::isOutOfGridRowCol(row, col, landUseMap))
+        return NODATA;
 
-    int caseIndex = int(landUseMap.value[row][col]);
-    if (caseIndex == int(landUseMap.header->flag))
-    {
-        //DEFAULT
-        caseIndex = 0;
-    }
+    if (isEqual(landUseMap.value[row][col], landUseMap.header->flag))
+        return NODATA;
 
-    return caseIndex;
+    return int(landUseMap.value[row][col]);
 }
 
 
 bool Vine3DProject::isVineyard(unsigned row, unsigned col)
 {
     int caseIndex = getModelCaseIndex(row, col);
+    if (caseIndex == NODATA)
+        return false;
+
     return (modelCases[caseIndex].landuse == landuse_vineyard);
 }
 
@@ -1458,17 +1460,20 @@ bool Vine3DProject::computeVine3DWaterSinkSource()
             surfaceIndex = long(indexMap.at(0).value[row][col]);
             if (surfaceIndex != long(indexMap.at(0).header->flag))
             {
-                // LAI
-                int idField = getModelCaseIndex(row, col);
-                float laiGrass = modelCases[idField].maxLAIGrass;
-                float laiVine = statePlantMaps->leafAreaIndexMap->value[row][col];
-                double laiTot = double(laiVine + laiGrass);
+                int caseIndex = getModelCaseIndex(row, col);
+                if (caseIndex != NODATA)
+                {
+                    // LAI
+                    float laiGrass = modelCases[caseIndex].maxLAIGrass;
+                    float laiVine = statePlantMaps->leafAreaIndexMap->value[row][col];
+                    double laiTot = double(laiVine + laiGrass);
 
-                int soilIndex = getSoilIndex(row, col);
+                    int soilIndex = getSoilIndex(row, col);
 
-                double realEvap = assignEvaporation(row, col, laiTot, soilIndex);       // [mm]
-                flow = area * (realEvap / 1000.0);                                      // [m3/h]
-                totalEvaporation += flow;
+                    double realEvap = assignEvaporation(row, col, laiTot, soilIndex);       // [mm]
+                    flow = area * (realEvap / 1000.0);                                      // [m3/h]
+                    totalEvaporation += flow;
+                }
             }
         }
     }
diff --git a/bin/VINE3D/vine3DProject.h b/bin/VINE3D/vine3DProject.h
index 65c5eea4f..9406e21ce 100644
--- a/bin/VINE3D/vine3DProject.h
+++ b/bin/VINE3D/vine3DProject.h
@@ -47,7 +47,6 @@
     {
 
     public:
-
         QString dbVine3DFileName;
         QSqlDatabase dbVine3D;
 
@@ -111,8 +110,8 @@
 
         int queryFieldPoint(double x, double y);
 
-        bool readFieldQuery(QSqlQuery &myQuery, int &idField, Crit3DLanduse &landuse, int &vineIndex, int &trainingIndex, float &maxLaiGrass,  float &maxIrrigationRate);
-        bool setField(int fieldIndex, int fieldId, Crit3DLanduse landuse, int soilIndex, int vineIndex, int trainingIndex,
+        bool readFieldQuery(QSqlQuery &myQuery, int &idField, GrapevineLanduse &landuse, int &vineIndex, int &trainingIndex, float &maxLaiGrass,  float &maxIrrigationRate);
+        bool setField(int fieldIndex, int fieldId, GrapevineLanduse landuse, int soilIndex, int vineIndex, int trainingIndex,
                             float maxLaiGrass,  float maxIrrigationRate);
         bool getFieldBookIndex(int firstIndex, QDate myQDate, int fieldIndex, int* outputIndex);
 
@@ -134,7 +133,7 @@
 
         float getTimeStep();
 
-        int getModelCaseIndex(unsigned row, unsigned col);
+        int getModelCaseIndex(int row, int col);
 
         bool isVineyard(unsigned row, unsigned col);
 
diff --git a/bin/VINE3D/waterBalance.cpp b/bin/VINE3D/waterBalance.cpp
index 1622101f6..e7e851296 100644
--- a/bin/VINE3D/waterBalance.cpp
+++ b/bin/VINE3D/waterBalance.cpp
@@ -1,6 +1,7 @@
 #include <math.h>
 #include <vector>
 #include "commonConstants.h"
+#include "basicMath.h"
 #include "gis.h"
 #include "dataHandler.h"
 #include "vine3DProject.h"
@@ -47,8 +48,8 @@ gis::Crit3DRasterGrid* Crit3DWaterBalanceMaps::getMapFromVar(criteria3DVariable
 double getSoilVar(Vine3DProject* myProject, int soilIndex, int layerIndex, soil::soilVariable myVar)
 {
     int horizonIndex = soil::getHorizonIndex(myProject->soilList[soilIndex], myProject->layerDepth[unsigned(layerIndex)]);
-
-    if (int(horizonIndex) == int(NODATA)) return NODATA;
+    if (horizonIndex == NODATA)
+        return NODATA;
 
     if (myVar == soil::soilWaterPotentialWP)
     {
@@ -88,41 +89,39 @@ double getSoilVar(Vine3DProject* myProject, int soilIndex, int layerIndex, soil:
 }
 
 
-double* getSoilVarProfile(Vine3DProject* myProject, int row, int col, soil::soilVariable myVar)
+std::vector<double> getSoilVarProfile(Vine3DProject* myProject, int row, int col, soil::soilVariable myVar)
 {
-    double* myProfile = static_cast<double*> (calloc(size_t(myProject->nrLayers), sizeof(double)));
+    std::vector<double> varProfile;
+    varProfile.resize(myProject->nrLayers);
 
     for (unsigned int layerIndex = 0; layerIndex < myProject->nrLayers; layerIndex++)
-        myProfile[layerIndex] = NODATA;
+        varProfile[layerIndex] = NODATA;
 
-    int nodeIndex;
     int soilIndex = myProject->getSoilIndex(row, col);
-
-    for (unsigned int layerIndex = 0; layerIndex < myProject->nrLayers; layerIndex++)
+    if (soilIndex != NODATA)
     {
-        nodeIndex = int(myProject->indexMap.at(size_t(layerIndex)).value[row][col]);
-
-        if (nodeIndex != int(myProject->indexMap.at(size_t(layerIndex)).header->flag))
+        for (unsigned int layerIndex = 0; layerIndex < myProject->nrLayers; layerIndex++)
         {
             if (myVar == soil::soilWaterPotentialFC || myVar == soil::soilWaterPotentialWP
                 || myVar == soil::soilWaterContentSat  || myVar == soil::soilWaterContentFC
                 || myVar == soil::soilWaterContentWP)
             {
-                myProfile[layerIndex] = getSoilVar(myProject, soilIndex, layerIndex, myVar);
+                varProfile[layerIndex] = getSoilVar(myProject, soilIndex, layerIndex, myVar);
             }
         }
     }
 
-    return myProfile;
+    return varProfile;
 }
 
 
-double* getCriteria3DVarProfile(Vine3DProject* myProject, int row, int col, criteria3DVariable myVar)
+std::vector<double> getCriteria3DVarProfile(Vine3DProject* myProject, int row, int col, criteria3DVariable myVar)
 {
-    double* myProfile = (double *) calloc(myProject->nrLayers, sizeof(double));
+    std::vector<double> varProfile;
+    varProfile.resize(myProject->nrLayers);
 
     for (unsigned int layerIndex = 0; layerIndex < myProject->nrLayers; layerIndex++)
-        myProfile[layerIndex] = NODATA;
+        varProfile[layerIndex] = NODATA;
 
     for (unsigned int layerIndex = 0; layerIndex < myProject->nrLayers; layerIndex++)
     {
@@ -130,11 +129,11 @@ double* getCriteria3DVarProfile(Vine3DProject* myProject, int row, int col, crit
 
         if (nodeIndex != myProject->indexMap.at(layerIndex).header->flag)
         {
-             myProfile[layerIndex] = getCriteria3DVar(myVar, nodeIndex);
+             varProfile[layerIndex] = getCriteria3DVar(myVar, nodeIndex);
         }
     }
 
-    return myProfile;
+    return varProfile;
 }
 
 
@@ -241,7 +240,6 @@ bool getRootZoneAWCmap(Vine3DProject* myProject, gis::Crit3DRasterGrid* outputMa
     long nodeIndex;
     double awc, thickness, sumAWC;
     int soilIndex, horizonIndex;
-    int modelCaseIndex;
 
     for (int row = 0; row < outputMap->header->nrRows; row++)
         for (int col = 0; col < outputMap->header->nrCols; col++)
@@ -253,9 +251,9 @@ bool getRootZoneAWCmap(Vine3DProject* myProject, gis::Crit3DRasterGrid* outputMa
             {
                 sumAWC = 0.0;
                 soilIndex = myProject->getSoilIndex(row, col);
-                modelCaseIndex = myProject->getModelCaseIndex(row,col);
+                int caseIndex = myProject->getModelCaseIndex(row,col);
 
-                if (soilIndex != NODATA)
+                if (soilIndex != NODATA && caseIndex != NODATA)
                 {
                     for (unsigned int layer = 1; layer < myProject->nrLayers; layer++)
                     {
@@ -263,7 +261,7 @@ bool getRootZoneAWCmap(Vine3DProject* myProject, gis::Crit3DRasterGrid* outputMa
 
                         if (nodeIndex != myProject->indexMap.at(layer).header->flag)
                         {
-                            if (myProject->grapevine.getRootDensity(&(myProject->modelCases[modelCaseIndex]), layer) > 0.0)
+                            if (myProject->grapevine.getRootDensity(&(myProject->modelCases[caseIndex]), layer) > 0.0)
                             {
                                 awc = soilFluxes3D::getAvailableWaterContent(nodeIndex);  //[m3 m-3]
                                 if (awc != NODATA)
diff --git a/bin/VINE3D/waterBalance.h b/bin/VINE3D/waterBalance.h
index 9a72193ef..2ed141f6a 100644
--- a/bin/VINE3D/waterBalance.h
+++ b/bin/VINE3D/waterBalance.h
@@ -48,8 +48,8 @@
     bool saveWaterBalanceCumulatedOutput(Vine3DProject* myProject, QDate myDate, criteria3DVariable myVar,
                                 QString varName, QString notes, QString outputPath);
 
-    double* getCriteria3DVarProfile(Vine3DProject* myProject, int myRow, int myCol, criteria3DVariable myVar);
-    double* getSoilVarProfile(Vine3DProject* myProject, int myRow, int myCol, soil::soilVariable myVar);
+    std::vector<double> getCriteria3DVarProfile(Vine3DProject* myProject, int myRow, int myCol, criteria3DVariable myVar);
+    std::vector<double> getSoilVarProfile(Vine3DProject* myProject, int myRow, int myCol, soil::soilVariable myVar);
 
 
 #endif // WATERBALANCE_H