diff --git a/DATA/PROJECT/Montue/DATA/crop_Montue.db b/DATA/PROJECT/Montue/DATA/crop_Montue.db
index f44f7f542..4e8aff6e4 100644
Binary files a/DATA/PROJECT/Montue/DATA/crop_Montue.db and b/DATA/PROJECT/Montue/DATA/crop_Montue.db differ
diff --git a/bin/CRITERIA3D/criteria3DProject.cpp b/bin/CRITERIA3D/criteria3DProject.cpp
index 235bc0e4c..97dd3e8f3 100644
--- a/bin/CRITERIA3D/criteria3DProject.cpp
+++ b/bin/CRITERIA3D/criteria3DProject.cpp
@@ -50,6 +50,7 @@ void Crit3DProcesses::initialize()
     computeMeteo = false;
     computeRadiation = false;
     computeWater = false;
+    computeEvaporation = false;
     computeCrop = false;
     computeSnow = false;
     computeSolutes = false;
@@ -227,6 +228,45 @@ void Crit3DProject::dailyUpdateCrop()
 }
 
 
+/*!
+ * \brief computeRealET
+ *  assign soil evaporation and crop transpiration for the whole domain
+ */
+void Crit3DProject::computeRealET()
+{
+    totalEvaporation = 0;
+    totalTranspiration = 0;
+
+    double area = DEM.header->cellSize * DEM.header->cellSize;
+
+    for (int row = 0; row < indexMap.at(0).header->nrRows; row++)
+    {
+        for (int col = 0; col < indexMap.at(0).header->nrCols; col++)
+        {
+            int surfaceIndex = indexMap.at(0).value[row][col];
+            if (surfaceIndex != indexMap.at(0).header->flag)
+            {
+                float lai = laiMap.value[row][col];
+                if (isEqual(lai, NODATA))
+                {
+                    lai = 0;
+                }
+
+                // assign real evaporation
+                double realEvap = assignEvaporation(row, col, lai);             // [mm]
+                double flow = area * (realEvap / 1000.);                        // [m3 h-1]
+                totalEvaporation += flow;
+
+                // assign real transpiration
+                double realTransp = assignTranspiration(row, col, lai);             // [mm]
+                flow = area * (realTransp / 1000.);                        // [m3 h-1]
+                totalTranspiration += flow;
+            }
+        }
+    }
+}
+
+
 bool Crit3DProject::runModels(QDateTime firstTime, QDateTime lastTime)
 {
     // initialize meteo
@@ -947,36 +987,55 @@ bool Crit3DProject::modelHourlyCycle(QDateTime myTime, const QString& hourlyOutp
         qApp->processEvents();
     }
 
-    if (processes.computeCrop)
+    if (processes.computeSnow)
     {
-        updateDailyTemperatures();
-
-        if (! hourlyMeteoMaps->computeET0PMMap(DEM, radiationMaps))
+        // check conflitti con evaporation
+        // check snowmelt -> surface H0
+        if (! computeSnowModel())
         {
             return false;
         }
+        qApp->processEvents();
+    }
+
+    // initalize sink / source
+    for (unsigned long i = 0; i < nrNodes; i++)
+    {
+        waterSinkSource.at(size_t(i)) = 0.;
+    }
+
+
+    if (processes.computeEvaporation || processes.computeCrop)
+    {
+        if (! hourlyMeteoMaps->computeET0PMMap(DEM, radiationMaps))
+            return false;
 
         if (isSaveOutputRaster())
         {
             saveHourlyMeteoOutput(referenceEvapotranspiration, hourlyOutputPath, myTime);
         }
+    }
 
-        qApp->processEvents();
-
-        // TODO compute evap/transp
+    if (processes.computeEvaporation && (! processes.computeCrop))
+    {
+        computeBareSoilEvaporation();
     }
 
-    if (processes.computeSnow)
+    if (processes.computeCrop)
     {
-        if (! computeSnowModel())
-            return false;
+        updateDailyTemperatures();
+
+        computeRealET();
+
         qApp->processEvents();
     }
 
     // soil water balance
     if (processes.computeWater)
     {
-        if (! computeWaterSinkSource())
+        assignPrecipitation();
+
+        if (! setSinkSource())
         {
             logError();
             return false;
diff --git a/bin/CRITERIA3D/criteria3DProject.h b/bin/CRITERIA3D/criteria3DProject.h
index 5b4a7ba9c..6bafbf728 100644
--- a/bin/CRITERIA3D/criteria3DProject.h
+++ b/bin/CRITERIA3D/criteria3DProject.h
@@ -28,7 +28,7 @@
     public:
 
         bool computeMeteo, computeRadiation, computeWater;
-        bool computeCrop, computeSnow, computeSolutes;
+        bool computeEvaporation, computeCrop, computeSnow, computeSolutes;
         bool computeHeat, computeAdvectiveHeat, computeLatentHeat;
 
         Crit3DProcesses();
@@ -71,6 +71,7 @@
         bool initializeCriteria3DModel();
         void initializeCrop();
         void dailyUpdateCrop();
+        void computeRealET();
 
         bool runModels(QDateTime firstTime, QDateTime lastTime);
 
diff --git a/bin/CRITERIA3D/mainwindow.cpp b/bin/CRITERIA3D/mainwindow.cpp
index 57a781ce8..87d0e4afb 100644
--- a/bin/CRITERIA3D/mainwindow.cpp
+++ b/bin/CRITERIA3D/mainwindow.cpp
@@ -1922,6 +1922,11 @@ bool MainWindow::startModels(QDateTime firstTime, QDateTime lastTime)
     if (myProject.processes.computeCrop)
     {
         // TODO: check on crop
+        if (myProject.landUnitList.size() == 0)
+        {
+            myProject.logError("load land units map before.");
+            return false;
+        }
     }
 
     // Load meteo data
@@ -2213,6 +2218,7 @@ void MainWindow::on_actionCriteria3D_compute_current_hour_triggered()
     myProject.processes.computeMeteo = true;
     myProject.processes.computeRadiation = true;
     myProject.processes.computeWater = true;
+    myProject.processes.computeEvaporation = true;
     myProject.processes.computeCrop = true;
 
     startModels(currentTime, currentTime);
@@ -2235,6 +2241,7 @@ void MainWindow::on_actionCriteria3D_run_models_triggered()
     myProject.processes.computeMeteo = true;
     myProject.processes.computeRadiation = true;
     myProject.processes.computeWater = false;
+    myProject.processes.computeEvaporation = true;
     myProject.processes.computeCrop = true;
 
     startModels(firstTime, lastTime);
diff --git a/bin/CRITERIA3D/mainwindow.ui b/bin/CRITERIA3D/mainwindow.ui
index 4d884ba40..ff9dd946a 100644
--- a/bin/CRITERIA3D/mainwindow.ui
+++ b/bin/CRITERIA3D/mainwindow.ui
@@ -1293,7 +1293,7 @@
        <second>0</second>
        <year>1752</year>
        <month>10</month>
-       <day>2</day>
+       <day>1</day>
       </datetime>
      </property>
      <property name="maximumDate">
@@ -1848,7 +1848,7 @@
      <x>0</x>
      <y>0</y>
      <width>1280</width>
-     <height>29</height>
+     <height>23</height>
     </rect>
    </property>
    <property name="palette">
@@ -2221,8 +2221,8 @@
    <addaction name="menuMeteo_points_2"/>
    <addaction name="menuDataSpatialization"/>
    <addaction name="menuSolar_radiation_model_2"/>
-   <addaction name="menuSnow_model"/>
    <addaction name="menuCriteria3D"/>
+   <addaction name="menuSnow_model"/>
    <addaction name="menuOutput_points"/>
    <addaction name="menuGIS_settings"/>
   </widget>
diff --git a/bin/CRITERIA3D/shared/project3D.cpp b/bin/CRITERIA3D/shared/project3D.cpp
index ba94694bd..163ab2380 100644
--- a/bin/CRITERIA3D/shared/project3D.cpp
+++ b/bin/CRITERIA3D/shared/project3D.cpp
@@ -983,7 +983,8 @@ bool Project3D::aggregateAndSaveDailyMap(meteoVariable myVar, aggregationMethod
 
 
 // compute evaporation [mm]
-double Project3D::computeEvaporation(int row, int col, double lai)
+// TODO check pond
+double Project3D::assignEvaporation(int row, int col, double lai)
 {
     double depthCoeff, thickCoeff, layerCoeff;
     double availableWater;                              // [mm]
@@ -1044,6 +1045,138 @@ double Project3D::computeEvaporation(int row, int col, double lai)
 }
 
 
+// assign real crop transpiration
+// return sum of crop transpiration over the soil column
+double Project3D::assignTranspiration(int row, int col, double lai)
+{
+    double soilColumnTranspirationSum = 0;
+
+    // check
+    /*if (idCrop == "" || ! isLiving) return 0;
+    if (roots.rootDepth <= roots.rootDepthMin) return 0;
+    if (roots.firstRootLayer == NODATA) return 0;
+    if (maxTranspiration < EPSILON) return 0;
+
+    double thetaWP;                                 // [m3 m-3] volumetric water content at Wilting Point
+    double cropWP;                                  // [mm] wilting point specific for crop
+    double waterSurplusThreshold;                        // [mm] water surplus stress threshold
+    double waterScarcityThreshold;                  // [mm] water scarcity stress threshold
+    double WSS;                                     // [] water surplus stress
+
+    double TRs=0.0;                                 // [mm] actual transpiration with only water scarsity stress
+    double TRe=0.0;                                 // [mm] actual transpiration with only water surplus stress
+    double totRootDensityWithoutStress = 0.0;       // [-]
+    double redistribution = 0.0;                    // [mm]
+
+    // initialize
+    unsigned int nrLayers = unsigned(soilLayers.size());
+    bool* isLayerStressed = new bool[nrLayers];
+    for (unsigned int i = 0; i < nrLayers; i++)
+    {
+        isLayerStressed[i] = false;
+        layerTranspiration[i] = 0;
+    }
+
+    // water surplus
+    if (isWaterSurplusResistant())
+        WSS = 0.0;
+    else
+        WSS = 0.5;
+
+    for (unsigned int i = unsigned(roots.firstRootLayer); i <= unsigned(roots.lastRootLayer); i++)
+    {
+        // [mm]
+        waterSurplusThreshold = soilLayers[i].SAT - (WSS * (soilLayers[i].SAT - soilLayers[i].FC));
+
+        thetaWP = soil::thetaFromSignPsi(-soil::cmTokPa(psiLeaf), *(soilLayers[i].horizonPtr));
+        // [mm]
+        cropWP = thetaWP * soilLayers[i].thickness * soilLayers[i].soilFraction * 1000.0;
+
+        // [mm]
+        waterScarcityThreshold = soilLayers[i].FC - fRAW * (soilLayers[i].FC - cropWP);
+
+        if ((soilLayers[i].waterContent - waterSurplusThreshold) > EPSILON)
+        {
+            // WATER SURPLUS
+            layerTranspiration[i] = maxTranspiration * roots.rootDensity[i] *
+                                    ((soilLayers[i].SAT - soilLayers[i].waterContent)
+                                     / (soilLayers[i].SAT - waterSurplusThreshold));
+
+            TRe += layerTranspiration[i];
+            TRs += maxTranspiration * roots.rootDensity[i];
+            isLayerStressed[i] = true;
+        }
+        else if (soilLayers[i].waterContent < waterScarcityThreshold)
+        {
+            // WATER SCARSITY
+            if (soilLayers[i].waterContent <= cropWP)
+            {
+                layerTranspiration[i] = 0;
+            }
+            else
+            {
+                layerTranspiration[i] = maxTranspiration * roots.rootDensity[i] *
+                                        ((soilLayers[i].waterContent - cropWP) / (waterScarcityThreshold - cropWP));
+            }
+
+            TRs += layerTranspiration[i];
+            TRe += maxTranspiration * roots.rootDensity[i];
+            isLayerStressed[i] = true;
+        }
+        else
+        {
+            // normal conditions
+            layerTranspiration[i] = maxTranspiration * roots.rootDensity[i];
+
+            TRs += layerTranspiration[i];
+            TRe += layerTranspiration[i];
+
+            if ((soilLayers[i].waterContent - layerTranspiration[i]) > waterScarcityThreshold)
+            {
+                isLayerStressed[i] = false;
+                totRootDensityWithoutStress +=  roots.rootDensity[i];
+            }
+            else
+            {
+                isLayerStressed[i] = true;
+            }
+        }
+    }
+
+    // WATER STRESS [-]
+    double firstWaterStress = 1 - (TRs / maxTranspiration);
+
+    // Hydraulic redistribution
+    // the movement of water from moist to dry soil through plant roots
+    // TODO add numerical process
+    if (firstWaterStress > EPSILON && totRootDensityWithoutStress > EPSILON)
+    {
+        // redistribution acts on not stressed roots
+        redistribution = MINVALUE(firstWaterStress, totRootDensityWithoutStress) * maxTranspiration;
+
+        for (int i = roots.firstRootLayer; i <= roots.lastRootLayer; i++)
+        {
+            if (! isLayerStressed[i])
+            {
+                double addLayerTransp = redistribution * (roots.rootDensity[unsigned(i)] / totRootDensityWithoutStress);
+                layerTranspiration[unsigned(i)] += addLayerTransp;
+                TRs += addLayerTransp;
+            }
+        }
+    }
+
+    waterStress = 1 - (TRs / maxTranspiration);
+
+    for (int i = roots.firstRootLayer; i <= roots.lastRootLayer; i++)
+    {
+        totalTranspiration += layerTranspiration[unsigned(i)];
+    }
+    */
+
+    return soilColumnTranspirationSum;
+}
+
+
 // input: timeStep [s]
 void Project3D::computeWaterBalance3D(double totalTimeStep)
 {
@@ -1155,17 +1288,10 @@ bool Project3D::interpolateAndSaveHourlyMeteo(meteoVariable myVar, const QDateTi
 }
 
 
-bool Project3D::computeWaterSinkSource()
+void Project3D::assignPrecipitation()
 {
     // initialize
     totalPrecipitation = 0;
-    totalEvaporation = 0;
-    totalTranspiration = 0;
-
-    for (unsigned long i = 0; i < nrNodes; i++)
-    {
-        waterSinkSource.at(size_t(i)) = 0.;
-    }
 
     double area = DEM.header->cellSize * DEM.header->cellSize;
 
@@ -1191,8 +1317,32 @@ bool Project3D::computeWaterSinkSource()
             }
         }
     }
+}
+
+
+bool Project3D::setSinkSource()
+{
+    for (unsigned long i = 0; i < nrNodes; i++)
+    {
+        int myResult = soilFluxes3D::setWaterSinkSource(signed(i), waterSinkSource.at(i));
+        if (isCrit3dError(myResult, errorString))
+        {
+            errorString = "waterBalanceSinkSource:" + errorString;
+            return false;
+        }
+    }
+
+    return true;
+}
+
+
+void Project3D::computeBareSoilEvaporation()
+{
+    // initialize
+    totalEvaporation = 0;
+
+    double area = DEM.header->cellSize * DEM.header->cellSize;
 
-    // evaporation
     for (int row = 0; row < indexMap.at(0).header->nrRows; row++)
     {
         for (int col = 0; col < indexMap.at(0).header->nrCols; col++)
@@ -1200,19 +1350,21 @@ bool Project3D::computeWaterSinkSource()
             int surfaceIndex = indexMap.at(0).value[row][col];
             if (surfaceIndex != indexMap.at(0).header->flag)
             {
-                // TODO crop
                 double lai = 0;
 
-                double realEvap = computeEvaporation(row, col, lai);            // [mm]
+                double realEvap = assignEvaporation(row, col, lai);             // [mm]
                 double flow = area * (realEvap / 1000.);                        // [m3 h-1]
                 totalEvaporation += flow;
             }
         }
     }
+}
 
-    // crop transpiration
-    // TODO
-    /*
+
+// crop transpiration
+// TODO
+/*
+        totalTranspiration = 0;
     for (unsigned int layerIndex=1; layerIndex < nrLayers; layerIndex++)
     {
         for (long row = 0; row < indexMap.at(size_t(layerIndex)).header->nrRows; row++)
@@ -1235,19 +1387,10 @@ bool Project3D::computeWaterSinkSource()
     }
     */
 
-    // assign sink/source
-    for (unsigned long i = 0; i < nrNodes; i++)
-    {
-        int myResult = soilFluxes3D::setWaterSinkSource(signed(i), waterSinkSource.at(i));
-        if (isCrit3dError(myResult, errorString))
-        {
-            errorString = "waterBalanceSinkSource:" + errorString;
-            return false;
-        }
-    }
 
-    return true;
-}
+
+
+
 
 
 bool Project3D::setCriteria3DMap(criteria3DVariable var, int layerIndex)
diff --git a/bin/CRITERIA3D/shared/project3D.h b/bin/CRITERIA3D/shared/project3D.h
index 2a3427faa..abc594589 100644
--- a/bin/CRITERIA3D/shared/project3D.h
+++ b/bin/CRITERIA3D/shared/project3D.h
@@ -142,8 +142,11 @@
                                       const QString& dailyPath, const QString& hourlyPath);
 
         bool interpolateHourlyMeteoVar(meteoVariable myVar, const QDateTime& myTime);
-        double computeEvaporation(int row, int col, double lai);
-        bool computeWaterSinkSource();
+        double assignEvaporation(int row, int col, double lai);
+        double assignTranspiration(int row, int col, double lai);
+        void computeBareSoilEvaporation();
+        void assignPrecipitation();
+        bool setSinkSource();
         void computeWaterBalance3D(double totalTimeStep);
         bool updateCrop(QDateTime myTime);