modified to import plant_uptake and export wn,sw,psi_m,aet

Author: jaideep777

src/EVAP.cpp

@@ -94,6 +94,9 @@ void EVAP::calculate_daily_fluxes(double sw, int n, int y, double sw_in,
     rn_d = d_sr.rn_d;
     rnn_d = d_sr.rnn_d;
     double ts = d_sr.ts;
+    
+    printf("Rnl in EVAP function: %0.6f W/m^2\n", rnl);
+
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // 1.1. Assume water temperature 0.0 if air temperature < 0.0
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

src/SOLAR.cpp

@@ -1,6 +1,7 @@
 #include <cmath>
 #include <stdio.h>
 #include <vector>
+#include <iostream>
 
 #include "global.h"
 #include "SOLAR.h"
@@ -84,6 +85,8 @@ void SOLAR::calculate_daily_fluxes(int n, int y, double sw_in, double tc, double
     Output:
     Features: Calculates the daily solar radiation fluxes
     *********************************************************************** */
+    cout << "SOLAR calc inputs: " << n << " " << y << " " << sw_in << " " << tc << " " << slop << " " << asp << " " << snow << " " << nd << " " << sw << '\n';
+
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // 0. Save day of year
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -263,6 +266,7 @@ void SOLAR::calculate_daily_fluxes(int n, int y, double sw_in, double tc, double
     //max air temeperature
     //ts = tc/dcos(hs/2);
 
+    cout << "End of SOLAR calc: rnl = " << rnl << endl;
 
 }
 
@@ -388,6 +392,7 @@ srad SOLAR::get_vals(){
     dsr.rn_d = rn_d;    // daytime net radiation, J/m^2
     dsr.rnn_d = rnn_d;  // nighttime net radiation, J/m^2
     dsr.ts = ts;        // surface temperature, C
+    dsr.ppfd_d = ppfd_d;  // daily mean ppfd
 
     return dsr;
 }

src/SPLASH.cpp

@@ -925,9 +925,9 @@ void SPLASH::quick_run(int n, int y, double wn, double sw_in, double tc,
 
 }
 
-void SPLASH::run_one_day(int n, int y, double wn, double sw_in, double tc,
+etr SPLASH::run_one_day(int n, int y, double wn, double sw_in, double tc,
                          double pn,smr &dsoil, double slop, double asp,double snow, double snowfall, vector <double> &soil_info,
-                         double qin, double td, double nd) {
+                         double qin, double td, double nd, double plant_uptake) {
     /* ***********************************************************************
     Name:     SPLASH.run_one_day
     Input:    - int, day of year (n)
@@ -1277,7 +1277,7 @@ void SPLASH::run_one_day(int n, int y, double wn, double sw_in, double tc,
     // 5.1.5. calculate Hortonian (infiltration excess) runoff
     double ro_h = max(inflow-infi,0.0);
     // 5.1.6. calculate recharge
-    double R = infi - dvap.aet;
+    double R = infi - dvap.aet - plant_uptake;
     // 5.1.7. define the hydraulic gradient assuming steady state: z: I=Ks(dh/dz +1), x: tan(slop), if theta<fc no vertical flow outside the column
 
     double Kunsat = Ksat_visc * pow((theta_m/theta_s),(3.0+(2.0/lambda)));
@@ -1589,10 +1589,21 @@ void SPLASH::run_one_day(int n, int y, double wn, double sw_in, double tc,
     dsoil.sqout = qin_nday; 
     dsoil.bflow = T;
     dsoil.nd = nd;
-    
-   
 
-       
+    dsoil.stress_factor = (sm-RES)/(Wmax-RES);   
+    cout << "sw_end(calc at start) / stress_end (calc at end) = " << sw << " / " << dsoil.stress_factor << '\n';
+
+    theta_i = (sm)/(depth*1000.0);
+    // correct theta_i for NA error reaching boundary conditions
+    if (theta_i>=theta_s){
+        theta_i = theta_s - 0.001;
+    } else if (theta_i<=theta_r){
+        theta_i = theta_r + 0.001;
+    }
+    dsoil.psi_m = bub_press/pow((((theta_i-theta_r)/(theta_s-theta_r))),(1/lambda));
+    dsoil.psi_m *= 0.00000980665; // convert from mmH2o --> MPa
+
+    return dvap;
 }
 
 vector<vector<double>> SPLASH::spin_up_cpp(int n, int y, vector<double> &sw_in, vector <double> &tair, vector <double> &pn, 

src/SPLASH.h

@@ -89,8 +89,8 @@ class SPLASH {
         void quick_run(int n, int y, double wn, double sw_in, double tc,
                        double pn, smr &dsm, double slop, double asp, double snow, double snowfall, vector <double> &soil_info,
                        double qin, double td, double nd);
-        void run_one_day(int n, int y, double wn, double sw_in, double tc,
-                         double pn, smr &dsoil, double slop, double asp,double snow, double snowfall, vector <double> &soil_info,double qin, double td, double nd);
+        etr run_one_day(int n, int y, double wn, double sw_in, double tc,
+                         double pn, smr &dsoil, double slop, double asp,double snow, double snowfall, vector <double> &soil_info,double qin, double td, double nd, double plant_uptake = 0);
 
         vector<vector<double>> spin_up_cpp(int n, int y, vector<double> &sw_in, vector <double> &tair, vector <double> &pn, double slop, double asp, vector <double> &snowfall, vector <double> &soil_info);        
 

src/smr.h

@@ -57,5 +57,7 @@ struct smr {
     double tdr;     // day left to drain totally water from upslope
     double nd;     // snow age (days)
     double pet;     // potential evapotrasnpiration mm
+    double stress_factor; // soil moisture stress factor (wn-RES)/(Wmax-RES), unitless
+    double psi_m;    // soil matrix potential, MPa
 };
 #endif

src/srad.h

@@ -50,5 +50,6 @@ struct srad {
     double rn_d;    // daytime net radiation, J/m^2
     double rnn_d;   // nighttime net radiation, J/m^2
     double ts;     // surface temperature
+    double ppfd_d;   // photosynthetic photon flux density, mol/m^2/day
 };
 #endif