Repository: parflow/parflow (language: python)

Issue:
ParFlow/TCL is not assigning SSAT values to a constant
Here is the TCL section where values are assigned to the Van Genuchten parameters:
```
#-----------------------------------------------------------------------------
# Saturation, VA [1/m] , VN [-] , SSat [-] , & SRes [-] 
#-----------------------------------------------------------------------------
pfset Phase.Saturation.GeomNames                                      "domain"
pfset Phase.Saturation.Type                                     "VanGenuchten"
pfset Phase.Saturation.VanGenuchten.File                                     0
pfset Geom.domain.Saturation.Alpha                                         1.5
pfset Geom.domain.Saturation.N                                             2.0
pfset Geom.domain.Saturation.SRes                                          0.1
pfset Geom.domain.Saturation.SSat                                          1.0
```

Here is the prove of the issue reading ParFlow output files in R:
```
> VA <- PriorityFlow::readpfb(filename = paste0(rd, "KWHM_parkinglot.out.alpha.pfb"), verbose = F)
> summary(VA)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    1.5     1.5     1.5     1.5     1.5     1.5 
> 
> VN <- PriorityFlow::readpfb(filename = paste0(rd, "KWHM_parkinglot.out.n.pfb"), verbose = F)
> summary(VN)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      2       2       2       2       2       2 
> 
> TR <- PriorityFlow::readpfb(filename = paste0(rd, "KWHM_parkinglot.out.sres.pfb"), verbose = F)
> summary(TR)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.1     0.1     0.1     0.1     0.1     0.1 
> 
> TS <- PriorityFlow::readpfb(filename = paste0(rd, "KWHM_parkinglot.out.ssat.pfb"), verbose = F)
> summary(TS) ### Issue
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       0       0       0       0       0
```
As you can see, SSAT values are zero despite being set to one in the TCL file.



Discussion / hints:
Thanks for the bug report.    This appears to an issue with writing ssat the file, the internal calculation looks to be correct.

I'll work on a patch.

Fix patch:
diff --git a/pfsimulator/parflow_lib/problem_saturation.c b/pfsimulator/parflow_lib/problem_saturation.c
index 30231738f..236745008 100644
--- a/pfsimulator/parflow_lib/problem_saturation.c
+++ b/pfsimulator/parflow_lib/problem_saturation.c
@@ -1299,6 +1299,7 @@ void  SaturationOutputStatic(
       alphas = (dummy1->alphas);
       ns = (dummy1->ns);
       s_ress = (dummy1->s_ress);
+      double* s_difs = (dummy1->s_difs);
       data_from_file = (dummy1->data_from_file);
 
       if (data_from_file == 0) /* Soil parameters given by region */
@@ -1338,10 +1339,13 @@ void  SaturationOutputStatic(
 	      double alpha = alphas[ir];
 	      double n = ns[ir];
 	      double s_res = s_ress[ir];
+	      double s_dif = s_difs[ir];
 	      
 	      pd_alpha_dat[ips] = alpha;  //BB
 	      pd_n_dat[ips] = n;  //BB
 	      pd_sres_dat[ips] = s_res;  //BB  // no ssat???
+	      // Storing s_dif in the stucture, convert back to s_sat for output
+	      pd_ssat_dat[ips] = s_dif + s_res;
 	      
 	    });
           }     /* End subgrid loop */


Test patch:
diff --git a/test/correct_output/crater2D_vangtable_linear.out.ssat.pfb b/test/correct_output/crater2D_vangtable_linear.out.ssat.pfb
index d3a5da4d0..6fe7b60ac 100644
Binary files a/test/correct_output/crater2D_vangtable_linear.out.ssat.pfb and b/test/correct_output/crater2D_vangtable_linear.out.ssat.pfb differ