From e8da3510ec2874a189fa488e5f6fdc94eb34fefa Mon Sep 17 00:00:00 2001
From: Magne Sjaastad <magne.sjaastad@ceetron.com>
Date: Fri, 12 Oct 2012 14:27:43 +0200
Subject: [PATCH] Compute SOIL if one or none of SGAS or SWAT is present p4#:
 19173

---
 .../RigReservoirCellResults.cpp               | 70 ++++++++++++++-----
 1 file changed, 51 insertions(+), 19 deletions(-)

diff --git a/ApplicationCode/ReservoirDataModel/RigReservoirCellResults.cpp b/ApplicationCode/ReservoirDataModel/RigReservoirCellResults.cpp
index 6aa4f5c3d5..c82a69e5b4 100644
--- a/ApplicationCode/ReservoirDataModel/RigReservoirCellResults.cpp
+++ b/ApplicationCode/ReservoirDataModel/RigReservoirCellResults.cpp
@@ -335,35 +335,67 @@ size_t RigReservoirCellResults::findScalarResultIndex(const QString& resultName)
 //--------------------------------------------------------------------------------------------------
 void RigReservoirCellResults::loadOrComputeSOIL()
 {
-    size_t resultGridIndex = findOrLoadScalarResult(RimDefines::DYNAMIC_NATIVE, "SOIL");
+    size_t soilResultGridIndex = findOrLoadScalarResult(RimDefines::DYNAMIC_NATIVE, "SOIL");
 
-    if (resultGridIndex == cvf::UNDEFINED_SIZE_T)
+    if (soilResultGridIndex == cvf::UNDEFINED_SIZE_T)
     {
+        const std::vector< std::vector<double> >* swat = NULL;
+        const std::vector< std::vector<double> >* sgas = NULL;
+
         size_t scalarIndexSWAT = findOrLoadScalarResult(RimDefines::DYNAMIC_NATIVE, "SWAT");
-        size_t scalarIndexSGAS = findOrLoadScalarResult(RimDefines::DYNAMIC_NATIVE, "SGAS");
-
-        if (scalarIndexSGAS != cvf::UNDEFINED_SIZE_T && scalarIndexSWAT != cvf::UNDEFINED_SIZE_T)
+        if (scalarIndexSWAT != cvf::UNDEFINED_SIZE_T)
         {
-            size_t timeStepCount = m_resultInfos[scalarIndexSWAT].m_timeStepDates.size();
-            resultGridIndex = addEmptyScalarResult(RimDefines::DYNAMIC_NATIVE, "SOIL");
-            m_cellScalarResults[resultGridIndex].resize(timeStepCount);
+            swat = &(cellScalarResults(scalarIndexSWAT));
+        }
 
-            const std::vector< std::vector<double> >& sgas = cellScalarResults(scalarIndexSGAS);
-            const std::vector< std::vector<double> >& swat = cellScalarResults(scalarIndexSWAT);
-            std::vector< std::vector<double> >& soil = cellScalarResults(resultGridIndex);
+        size_t scalarIndexSGAS = findOrLoadScalarResult(RimDefines::DYNAMIC_NATIVE, "SGAS");
+        if (scalarIndexSGAS != cvf::UNDEFINED_SIZE_T)
+        {
+            sgas = &(cellScalarResults(scalarIndexSGAS));
+        }
 
-            size_t resultValueCount = sgas[0].size();
+        size_t soilResultValueCount = 0;
+        size_t soilTimeStepCount = 0;
+        if (swat)
+        {
+            soilResultValueCount = swat->at(0).size();
+            soilTimeStepCount = m_resultInfos[scalarIndexSWAT].m_timeStepDates.size();
+        }
 
-            int timeStepIdx = 0;
-            for (timeStepIdx = 0; timeStepIdx < static_cast<int>(timeStepCount); timeStepIdx++)
-            {
-                soil[timeStepIdx].resize(resultValueCount);
-                int idx = 0;
+        if (sgas)
+        {
+            soilResultValueCount = qMax(soilResultValueCount, sgas->at(0).size());
+            
+            size_t sgasTimeStepCount = m_resultInfos[scalarIndexSGAS].m_timeStepDates.size();
+            soilTimeStepCount = qMax(soilTimeStepCount, sgasTimeStepCount);
+        }
+
+        soilResultGridIndex = addEmptyScalarResult(RimDefines::DYNAMIC_NATIVE, "SOIL");
+        m_cellScalarResults[soilResultGridIndex].resize(soilTimeStepCount);
+
+        std::vector< std::vector<double> >& soil = cellScalarResults(soilResultGridIndex);
+
+        int timeStepIdx = 0;
+        for (timeStepIdx = 0; timeStepIdx < static_cast<int>(soilTimeStepCount); timeStepIdx++)
+        {
+            soil[timeStepIdx].resize(soilResultValueCount);
+
+            int idx = 0;
 #pragma omp parallel for
-                for (idx = 0; idx < static_cast<int>(resultValueCount); idx++)
+            for (idx = 0; idx < static_cast<int>(soilResultValueCount); idx++)
+            {
+                double soilValue = 1.0;
+                if (sgas)
                 {
-                    soil[timeStepIdx][idx] = 1.0 - sgas[timeStepIdx][idx] - swat[timeStepIdx][idx];
+                    soilValue -= sgas->at(timeStepIdx)[idx];
                 }
+
+                if (swat)
+                {
+                    soilValue -= swat->at(timeStepIdx)[idx];
+                }
+
+                soil[timeStepIdx][idx] = soilValue;
             }
         }
     }