Merge pull request dailyerosion#255 from wils0674/patch-8

Update magic.R

scripts/RT/magic.R

@@ -109,6 +109,15 @@ Rockpct<-Rockpct/100 #don't think I need this, but including here for consistenc
 #compute silt by difference
 Siltpct<-round((1-(Sandpct+Claypct)),digits=2)
 
+#wilting point & saturation calculations using pedotransfer functions from Saxton and Rawls 2006
+wilting_point_i<-(-0.024*Sandpct+0.487*Claypct+0.006*Organicpct+0.0005*(Sandpct*Organicpct)-0.013(Claypct*Organicpct)+0.068*(Sandpct*Claypct)+0.31)
+wilting_point<-(wilting_point_i+(0.14*wilting_point_i-0.02))
+theta_s1<-0.278*Sandpct+0.034*Claypct+0.022*Organicpct-0.018*Sandpct*Claypct-0.027*Claypct*Organicpct-0.584*Sandpct*Claypct+0.078
+theta_s2<-theta_s1+0.636*theta_s1-0.107
+field_capacity1<-(-0.251*Sandpct+0.195*Claypct+0.011*Organicpct+0.006*Sandpct*Organicpct-0.027*Claypct*Organicpct+0.452*Sandpct*Claypct+0.299)
+field_capacity<-field_capacity1+(1.283*field_capacity1^2)-0.374*field_capacity1-0.015
+theta_s<-field_capacity+theta_s2-0.097*Sandpct+0.045
+
 # Compute soil CaCO3 based on pedotransfer function found in:
 # Wang, Q., Y. Li, and W. Klassen. 2005. Determination of Cation Exchange Capacity on Low to Highly
 # Calcareous Soils. Communications in Soil Science and Pland Anlaysis, 36: 1479-1498.
@@ -311,10 +320,11 @@ if (SEag < 0.1){
   SoilMass<-(Depthmm/10)*Soilpb  #soil mass in top layer from SSURGO file (not necessarily the same layer depth as WEPP output)
   GravWaterContent<-round((Volwater/Soilpb),digits=2) #assumes water densit of 1 g/cm3
   #GravWaterContent<-round((SoilWater/SoilMass),digits=2) #old calculation
-  wilting_point<-0.22 #JOhnston soil=0.123 probably volumetric #bearden silt loam =0.22; crete 0.11
-  field_capacity<-0.41 #Johnston=0.216 volumetric; 0.41 bearden silt loam? SSURGO value is 0.34; crete 0.24
+  #7/3/24 turned off wilting_point<-0.22 #JOhnston soil=0.123 probably volumetric #bearden silt loam =0.22; crete 0.11
+  #7/3/24 turned off field_capacity<-0.41 #Johnston=0.216 volumetric; 0.41 bearden silt loam? SSURGO value is 0.34; crete 0.24
   HRwcw<-(wilting_point)/Soilpb
-  HRwcs<-(field_capacity)/Soilpb #needs to use saturated content not FC
+  HRwcs<-(theta_s)/Soilpb
+  #HRwcs<-(field_capacity)/Soilpb #needs to use saturated content not FC
   HRwc<-(GravWaterContent-HRwcw)/(HRwcs-HRwcw) #non-dimensional soil moisture content
   if (HRwc < 0){
     HRwc <- 0