/*--------------------------------------------------------------*/
/* 											*/
/*					compute_subsurface_routing_hourly		*/
/*											*/
/*	compute_subsurface_routing_hourly.c - do subsurface computation during the end of each hour	*/
/*											*/
/*	NAME										*/
/*	compute_subsurface_routing_hourly.c - do subsurface computation during the end of each hour	*/
/*											*/
/*	SYNOPSIS									*/
/*	void compute_subsurface_routing_hourly( 								*/
/*						struct command_line_object command_line, */
/*							struct basin_object *basin)	*/
/*				 			int,			 	*/
/*							struct date *current_date)	*/
/*											*/
/* 											*/
/*											*/
/*	OPTIONS										*/
/*											*/
/*											*/
/*	DESCRIPTION									*/
/*	this function is called in basin_hourly_test at the end of each hour during	*/
/* 	hourly calculation, it is doing the compute_subsurface_routing for each hour	*/
/*											*/
/*											*/
/*											*/
/*	PROGRAMMER NOTES								*/
/*											*/
/*											*/
/*											*/
/*--------------------------------------------------------------*/
#include <stdio.h>
#include "rhessys.h"

void compute_subsurface_routing_hourly(
		struct command_line_object *command_line,
		struct basin_object *basin, 
		int n_timesteps, 
		struct date current_date) {
	/*--------------------------------------------------------------*/
	/*	Local function definition.				*/
	/*--------------------------------------------------------------*/

	void update_drainage_stream(struct patch_object *,
			struct command_line_object *, double, int);

	void update_drainage_road(struct patch_object *,
			struct command_line_object *, double, int);

	void update_drainage_land(struct patch_object *,
			struct command_line_object *, double, int);

	double compute_infiltration(int, double, double, double, double, double,
			double, double, double, double, double);

	double compute_z_final(int, double, double, double, double, double);

    double compute_N_leached(
                             int verbose_flag,
                             double total_nitrate,
                             double Qout,
                             double N_decay_rate,
                             double activedepthz,
                             double N_absorption_rate,
                             int signal,
                             struct patch_object *patch);

	double compute_layer_field_capacity(int, int, double, double, double,
			double, double, double, double, double, double);

	double compute_unsat_zone_drainage(int, int, double, double, double, double,
			double, double, double);

	/*--------------------------------------------------------------*/
	/*	Local variable definition.				*/
	/*--------------------------------------------------------------*/
	int i, d;
	int j, k;
	int grow_flag, verbose_flag;
	double time_int, tmp;
	double theta, m, Ksat, Nout;
	double NO3_out, NH4_out, DON_out, DOC_out;
	double return_flow, excess;
	double water_balance, infiltration;
	double innundation_depth;
	double basin_outflow;
	double basin_rz_storage;
	double basin_unsat_storage;
	double basin_sat_deficit;
	double basin_return_flow;
	double basin_detention_store;
	double basin_area;
	double preday_basin_unsat_storage;
	double preday_basin_rz_storage;
	double preday_basin_sat_deficit;
	double preday_sat_deficit;
	double preday_basin_return_flow;
	double preday_basin_detention_store;	
	double add_field_capacity, rz_drainage, unsat_drainage;
	double streamflow, Qout, Qin_total, Qstr_total;
	struct patch_object *patch;
	struct hillslope_object *hillslope;
	struct patch_object *neigh;
	struct litter_object *litter;
    double totalfc, temp;
	d=0;
	/*--------------------------------------------------------------*/
	/*	initializations						*/
	/*--------------------------------------------------------------*/
	
		grow_flag = command_line[0].grow_flag;
		verbose_flag = command_line[0].verbose_flag;

		time_int = 1.0 / n_timesteps;

	if (current_date.hour==1)
	{
		basin_outflow = 0.0;
		basin_area = 0.0;
		basin_unsat_storage = 0.0;
		basin_rz_storage = 0.0;
		basin_sat_deficit = 0.0;
		basin_return_flow = 0.0;
		basin_detention_store = 0.0;
		preday_basin_rz_storage = 0.0;
		preday_basin_unsat_storage = 0.0;
		preday_basin_sat_deficit = 0.0;
		preday_basin_return_flow = 0.0;
		preday_basin_detention_store = 0.0;

		basin[0].basin_outflow = 0.0;
		basin[0].basin_area = 0.0;
		basin[0].basin_unsat_storage = 0.0;
		basin[0].basin_rz_storage = 0.0;
		basin[0].basin_sat_deficit = 0.0;
		basin[0].basin_return_flow = 0.0;
		basin[0].basin_detention_store = 0.0;
		basin[0].preday_basin_rz_storage = 0.0;
		basin[0].preday_basin_unsat_storage = 0.0;
		basin[0].preday_basin_sat_deficit = 0.0;
		basin[0].preday_basin_return_flow = 0.0;
		basin[0].preday_basin_detention_store = 0.0;
		streamflow = 0.0;
		Qin_total = 0.0;
		Qstr_total = 0.0;
		d = 0;
		// Note: this assumes that the set of patches in the surface routing table is identical to
		//       the set of patches in the subsurface flow table
		for (i = 0; i < basin->route_list->num_patches; i++) {
			patch = basin->route_list->list[i];

			patch[0].streamflow = 0.0;
			patch[0].return_flow = 0.0;
			patch[0].base_flow = 0.0;
			patch[0].infiltration_excess = 0.0;
			basin[0].preday_basin_rz_storage += patch[0].rz_storage * patch[0].area;
			basin[0].preday_basin_unsat_storage += patch[0].unsat_storage * patch[0].area;
			basin[0].preday_basin_sat_deficit += patch[0].sat_deficit * patch[0].area;
			basin[0].preday_basin_return_flow += patch[0].return_flow * patch[0].area;
			basin[0].preday_basin_detention_store += patch[0].detention_store
					* patch[0].area;
			basin[0].basin_area += patch[0].area;
			patch[0].Qin_total = 0.0;
			patch[0].Qout_total = 0.0;
			patch[0].Qin = 0.0;
			patch[0].Qout = 0.0;
			patch[0].surface_Qin = 0.0;
			patch[0].surface_Qout = 0.0;



			patch[0].interim_sat = patch[0].sat_deficit - patch[0].unsat_storage;
//			if ((patch[0].sat_deficit - patch[0].unsat_storage) < ZERO)
//				patch[0].S = 1.0;
//			else
//				patch[0].S = patch[0].unsat_storage / patch[0].sat_deficit;

			if (grow_flag > 0) {
				patch[0].soil_ns.NO3_Qin = 0.0;
				patch[0].soil_ns.NO3_Qout = 0.0;
				patch[0].soil_ns.NH4_Qin = 0.0;
				patch[0].soil_ns.NH4_Qout = 0.0;
				patch[0].soil_ns.NO3_Qin_total = 0.0;
				patch[0].soil_ns.NO3_Qout_total = 0.0;
				patch[0].soil_ns.NH4_Qin_total = 0.0;
				patch[0].soil_ns.NH4_Qout_total = 0.0;
				patch[0].streamflow_DON = 0.0;
				patch[0].streamflow_DOC = 0.0;
				patch[0].streamflow_NO3 = 0.0;
				patch[0].streamflow_NH4 = 0.0;
				patch[0].soil_ns.DON_Qin_total = 0.0;
				patch[0].soil_ns.DON_Qout_total = 0.0;
				patch[0].soil_cs.DOC_Qin_total = 0.0;
				patch[0].soil_cs.DOC_Qout_total = 0.0;
				patch[0].surface_DON_Qin_total = 0.0;
				patch[0].surface_DON_Qout_total = 0.0;
				patch[0].surface_DOC_Qin_total = 0.0;
				patch[0].surface_DOC_Qout_total = 0.0;
				patch[0].soil_ns.leach = 0.0;
				patch[0].surface_ns_leach = 0.0;
				patch[0].soil_ns.DON_Qout = 0.0;
				patch[0].soil_ns.DON_Qin = 0.0;
				patch[0].soil_cs.DOC_Qout = 0.0;
				patch[0].soil_cs.DOC_Qin = 0.0;
				patch[0].surface_DON_Qout = 0.0;
				patch[0].surface_DON_Qin = 0.0;
				patch[0].surface_DOC_Qout = 0.0;
				patch[0].surface_DOC_Qin = 0.0;
			
				patch[0].streamNO3_from_surface	= 0.0;
				patch[0].streamNO3_from_sub = 0.0;
			}
		}
	}
	/*--------------------------------------------------------------*/
	/*	calculate Qout for each patch and add appropriate	*/
	/*	proportion of subsurface outflow to each neighbour	*/
	/*--------------------------------------------------------------*/
		for (i = 0; i < basin->route_list->num_patches; i++) {
			patch = basin->route_list->list[i];
						
//			patch[0].preday_sat_deficit = patch[0].sat_deficit;
//			patch[0].preday_sat_deficit_z = compute_z_final(verbose_flag,
//					patch[0].soil_defaults[0][0].porosity_0,
//					patch[0].soil_defaults[0][0].porosity_decay,
//					patch[0].soil_defaults[0][0].soil_depth, 0.0,
//					-1.0 * patch[0].sat_deficit);
			
            if(patch[0].soil_ns.nitrate!=patch[0].soil_ns.nitrate || patch[0].soil_ns.nitrate<0 ||
               patch[0].soil_ns.sminn!=patch[0].soil_ns.sminn || patch[0].soil_ns.sminn<0 ||
               patch[0].soil_cs.DOC!=patch[0].soil_cs.DOC || patch[0].soil_cs.DOC<0 ||
               patch[0].sat_NO3!=patch[0].sat_NO3 || patch[0].sat_NO3<0 ||
               patch[0].sat_NH4!=patch[0].sat_NH4 || patch[0].sat_NH4<0 ||
               patch[0].sat_DOC!=patch[0].sat_DOC || patch[0].sat_DOC<0 ||
               patch[0].surface_NO3!=patch[0].surface_NO3 || patch[0].surface_NO3<0 ||
               patch[0].surface_NH4!=patch[0].surface_NH4 || patch[0].surface_NH4<0 ||
               patch[0].surface_DOC!=patch[0].surface_DOC || patch[0].surface_DOC<0){
                printf("sub hourly routing1 (%d,%d) [%e %e %e] [%e %e %e] [%e %e %e]\n",
                       patch[0].ID, k,
                       patch[0].soil_ns.nitrate,patch[0].soil_ns.sminn,patch[0].soil_cs.DOC,
                       patch[0].sat_NO3,patch[0].sat_NH4,patch[0].sat_DOC,
                       patch[0].surface_NO3,patch[0].surface_NH4,patch[0].surface_DOC);
            }//debug
            if(patch[0].ID==193917 || patch[0].ID==131724 || patch[0].ID==182853 || patch[0].ID==167273){
                printf("sub hourly routing1 (forced) (%d,%d) [%e %e %e] [%e %e %e] [%e %e %e]\n",
                       patch[0].ID, k,
                       patch[0].soil_ns.nitrate,patch[0].soil_ns.sminn,patch[0].soil_cs.DOC,
                       patch[0].sat_NO3,patch[0].sat_NH4,patch[0].sat_DOC,
                       patch[0].surface_NO3,patch[0].surface_NH4,patch[0].surface_DOC);
            }//debug
            
            if(patch[0].sat_deficit >= 0){
                patch[0].sat_deficit = min(patch[0].sat_deficit,patch[0].soil_defaults[0][0].soil_water_cap);
                patch[0].preday_sat_deficit = patch[0].sat_deficit;
                patch[0].preday_sat_deficit_z = patch[0].sat_def_pct_indexM * patch[0].soil_defaults[0][0].sat_def_z[patch[0].sat_def_pct_index+1] + (1.0-patch[0].sat_def_pct_indexM) * patch[0].soil_defaults[0][0].sat_def_z[patch[0].sat_def_pct_index];
                temp = patch[0].sat_def_pct_indexM * patch[0].soil_defaults[0][0].fc1_0z[patch[0].sat_def_pct_index+1] + (1.0-patch[0].sat_def_pct_indexM) * patch[0].soil_defaults[0][0].fc1_0z[patch[0].sat_def_pct_index]; // total fc before ET water consumption to SAT
            }else{
                // surface
                patch[0].preday_sat_deficit = patch[0].sat_deficit;
                patch[0].preday_sat_deficit_z = patch[0].sat_deficit;
                temp = patch[0].soil_defaults[0][0].fc1_0z[0];
            }
            
            
            patch[0].hourly_subsur2stream_flow = 0;
			patch[0].hourly_sur2stream_flow = 0;
			patch[0].hourly_stream_flow = 0;
			patch[0].hourly[0].streamflow_NO3 = 0;
			patch[0].hourly[0].streamflow_NO3_from_sub = 0;
			patch[0].hourly[0].streamflow_NO3_from_surface = 0;
		}

		for (i = 0; i < basin->route_list->num_patches; i++) {
			patch = basin->route_list->list[i];
			litter=&(patch[0].litter);
			/*--------------------------------------------------------------*/
			/*	for roads, saturated throughflow beneath road cut	*/
			/*	is routed to downslope patches; saturated throughflow	*/
			/*	above the cut and overland flow is routed to the stream	*/
			/*								*/
			/*	for streams, no routing - all exported from basin	*/
			/*								*/
			/*	regular land patches - route to downslope neighbours    */
			/*--------------------------------------------------------------*/
			if ((patch[0].drainage_type == ROAD) && (command_line[0].road_flag == 1)) {
				update_drainage_road(patch, command_line, time_int,
						verbose_flag);
			} else if (patch[0].drainage_type == STREAM) {
				update_drainage_stream(patch, command_line, time_int,
						verbose_flag);
			} else {
				update_drainage_land(patch, command_line, time_int,
						verbose_flag);
			}



		} /* end i */

		/*--------------------------------------------------------------*/
		/*	update soil moisture and nitrogen stores		*/
		/*	check water balance					*/
		/*--------------------------------------------------------------*/
		for (i = 0; i < basin->route_list->num_patches; i++) {
			patch = basin->route_list->list[i];

			/*--------------------------------------------------------------*/
			/*	update subsurface 				*/
			/*-------------------------------------------------------------------------*/

			/*-------------------------------------------------------------------------*/
			/*	Recompute current actual depth to water table				*/
			/*-------------------------------------------------------------------------*/
			patch[0].sat_deficit += (patch[0].Qout - patch[0].Qin); // this part need to put into some where else
            if(patch[0].sat_deficit>patch[0].soil_defaults[0][0].soil_water_cap){
                patch[0].sat_deficit=patch[0].soil_defaults[0][0].soil_water_cap;
            }//if
            
//			patch[0].sat_deficit_z = compute_z_final(verbose_flag,
//					patch[0].soil_defaults[0][0].porosity_0,
//					patch[0].soil_defaults[0][0].porosity_decay,
//					patch[0].soil_defaults[0][0].soil_depth, 0.0,
//					-1.0 * patch[0].sat_deficit);
            
            if(patch[0].sat_deficit >= 0){
                patch[0].sat_deficit = min(patch[0].sat_deficit,patch[0].soil_defaults[0][0].soil_water_cap);
                patch[0].available_soil_water = patch[0].soil_defaults[0][0].soil_water_cap - patch[0].sat_deficit;
                patch[0].sat_def_pct = patch[0].sat_deficit * patch[0].soil_defaults[0][0].max_sat_def_1;
                patch[0].sat_def_pct_index = (int)(patch[0].sat_def_pct*1000);
                patch[0].sat_def_pct_indexM = 1000*(patch[0].sat_def_pct - patch[0].sat_def_pct_index*0.001);
                
                patch[0].sat_deficit_z = patch[0].sat_def_pct_indexM * patch[0].soil_defaults[0][0].sat_def_z[patch[0].sat_def_pct_index+1] + (1.0-patch[0].sat_def_pct_indexM) * patch[0].soil_defaults[0][0].sat_def_z[patch[0].sat_def_pct_index];
            }else{
                // surface
                patch[0].available_soil_water = patch[0].soil_defaults[0][0].soil_water_cap;
                patch[0].sat_deficit_z = patch[0].sat_deficit;
                patch[0].sat_def_pct = 0.0;
                patch[0].sat_def_pct_index = 0;
                patch[0].sat_def_pct_indexM = 0;
            }
            // fc & SatPct
            totalfc = patch[0].sat_def_pct_indexM * patch[0].soil_defaults[0][0].fc1_0z[patch[0].sat_def_pct_index+1] + (1.0-patch[0].sat_def_pct_indexM) * patch[0].soil_defaults[0][0].fc1_0z[patch[0].sat_def_pct_index];
            //totalfc *= (1.0-patch[0].basementFrac); // <---- second thought on this, Oct 8, 2019; basement is 3m at most
            
            if (patch[0].sat_deficit < ZERO) {
                //patch[0].aboveWT_SatPct = 1.0;
                //patch[0].rootzone.SatPct = 1.0;
                patch[0].rootzone.field_capacity = 0.0;
                patch[0].field_capacity = 0.0;
            } else {
                patch[0].rootzone.field_capacity = totalfc * patch[0].zeroRootCoef * (patch[0].rootzone_scale_ref*patch[0].rootzone_end_reffc[patch[0].sat_def_pct_index] + (1.0-patch[0].rootzone_scale_ref)*patch[0].rootzone_start_reffc[patch[0].sat_def_pct_index]);
                patch[0].rootzone.field_capacity = min(patch[0].rootzone.field_capacity,patch[0].rootzone.potential_sat);
                patch[0].field_capacity = max(0.0,min(patch[0].sat_deficit-patch[0].rootzone.potential_sat, totalfc - patch[0].rootzone.field_capacity));
            }//if else
            
			if (grow_flag > 0) {
				patch[0].soil_ns.nitrate += (patch[0].soil_ns.NO3_Qin
						- patch[0].soil_ns.NO3_Qout);
				patch[0].soil_ns.sminn += (patch[0].soil_ns.NH4_Qin
						- patch[0].soil_ns.NH4_Qout);
				patch[0].soil_cs.DOC += (patch[0].soil_cs.DOC_Qin
						- patch[0].soil_cs.DOC_Qout);
				patch[0].soil_ns.DON += (patch[0].soil_ns.DON_Qin
						- patch[0].soil_ns.DON_Qout);
			}

			/*--------------------------------------------------------------*/
			/*      Recompute 	soil moisture storage                   */
			/*--------------------------------------------------------------*/

//			if (patch[0].sat_deficit > patch[0].rootzone.potential_sat) {
//				patch[0].rootzone.S =
//						min(patch[0].rz_storage / patch[0].rootzone.potential_sat, 1.0);
//				patch[0].S = patch[0].unsat_storage
//						/ (patch[0].sat_deficit
//								- patch[0].rootzone.potential_sat);
//			} else {
//				patch[0].rootzone.S =
//						min((patch[0].rz_storage + patch[0].rootzone.potential_sat - patch[0].sat_deficit)
//								/ patch[0].rootzone.potential_sat, 1.0);
//				patch[0].S =
//						min(patch[0].rz_storage / patch[0].sat_deficit, 1.0);
//			}

			/*--------------------------------------------------------------*/
			/*	reset iterative  patch fluxes to zero			*/
			/*--------------------------------------------------------------*/
			patch[0].soil_ns.leach += (patch[0].soil_ns.DON_Qout
					+ patch[0].soil_ns.NH4_Qout + patch[0].soil_ns.NO3_Qout
					- patch[0].soil_ns.NH4_Qin - patch[0].soil_ns.NO3_Qin
					- patch[0].soil_ns.DON_Qin);
			patch[0].surface_ns_leach += ((patch[0].surface_NO3_Qout
					- patch[0].surface_NO3_Qin)
					+ (patch[0].surface_NH4_Qout - patch[0].surface_NH4_Qin)
					+ (patch[0].surface_DON_Qout - patch[0].surface_DON_Qin));
			patch[0].Qin_total += patch[0].Qin + patch[0].surface_Qin;
			patch[0].Qout_total += patch[0].Qout + patch[0].surface_Qout;

			patch[0].surface_Qin = 0.0;
			patch[0].surface_Qout = 0.0;
			patch[0].Qin = 0.0;
			patch[0].Qout = 0.0;
			if (grow_flag > 0) {
				patch[0].soil_cs.DOC_Qin_total += patch[0].soil_cs.DOC_Qin;
				patch[0].soil_cs.DOC_Qout_total += patch[0].soil_cs.DOC_Qout;
				patch[0].soil_ns.NH4_Qin_total += patch[0].soil_ns.NH4_Qin;
				patch[0].soil_ns.NH4_Qout_total += patch[0].soil_ns.NH4_Qout;
				patch[0].soil_ns.NO3_Qin_total += patch[0].soil_ns.NO3_Qin;
				patch[0].soil_ns.NO3_Qout_total += patch[0].soil_ns.NO3_Qout;
				patch[0].soil_ns.DON_Qin_total += patch[0].soil_ns.DON_Qin;
				patch[0].soil_ns.DON_Qout_total += patch[0].soil_ns.DON_Qout;
				patch[0].surface_DON_Qin_total += patch[0].surface_DON_Qin;
				patch[0].surface_DON_Qout_total += patch[0].surface_DON_Qout;
				patch[0].surface_DOC_Qin_total += patch[0].surface_DOC_Qin;
				patch[0].surface_DOC_Qout_total += patch[0].surface_DOC_Qout;

				patch[0].soil_ns.NH4_Qin = 0.0;
				patch[0].soil_ns.NH4_Qout = 0.0;
				patch[0].soil_ns.NO3_Qin = 0.0;
				patch[0].soil_ns.NO3_Qout = 0.0;
				patch[0].soil_ns.DON_Qout = 0.0;
				patch[0].soil_ns.DON_Qin = 0.0;
				patch[0].soil_cs.DOC_Qout = 0.0;
				patch[0].soil_cs.DOC_Qin = 0.0;
				patch[0].surface_NH4_Qout = 0.0;
				patch[0].surface_NH4_Qin = 0.0;
				patch[0].surface_NO3_Qout = 0.0;
				patch[0].surface_NO3_Qin = 0.0;
				patch[0].surface_DON_Qout = 0.0;
				patch[0].surface_DON_Qin = 0.0;
				patch[0].surface_DOC_Qout = 0.0;
				patch[0].surface_DOC_Qin = 0.0;

			}
			/*--------------------------------------------------------------*/
			/*	finalize streamflow and saturation deficits		*/
			/*								*/
			/*	move any saturation excess into detention store		*/
			/*	(i.e this needs to be routed on the next time step)	*/
			/* 	some streamflow may have already been accumulated from 	*/
			/* 	redirected streamflow					*/
			/*	water balance calculations				*/
			/* only on last iteration					*/
			/* **** note that streamflow is updated sequentially		*/
			/*	i.e not at the end; it is similar to Qout, in		*/
			/*	that it accumulates flux in from patches		*/
			/*	(roads) that direct water to the stream			*/
			/*--------------------------------------------------------------*/
            double excess = patch[0].rz_storage+patch[0].unsat_storage - patch[0].sat_deficit;// + patch[0].constraintWaterTableTopDepth_def;
		      	if (excess > 0) {
				  
					patch[0].detention_store += excess;
					patch[0].sat_deficit -= patch[0].unsat_storage + patch[0].rz_storage - excess;
					patch[0].unsat_storage = 0.0;
					patch[0].rz_storage = 0.0;

				
					if (grow_flag > 0) {
						Nout =
								compute_N_leached(verbose_flag,
										patch[0].sat_DOC,//soil_cs.DOC,
                                        excess,
                                        patch[0].soil_defaults[0][0].DOMdecayRate,
                                        patch[0].soil_defaults[0][0].active_zone_z,
                                        patch[0].soil_defaults[0][0].DOC_adsorption_rate,
										26,patch);
						patch[0].surface_DOC += Nout;
						patch[0].soil_cs.DOC -= Nout;
					}
	
					if (grow_flag > 0) {
						Nout =
								compute_N_leached(verbose_flag,
										patch[0].sat_DON,//soil_ns.DON,
                                        excess,
                                        patch[0].soil_defaults[0][0].DOMdecayRate,
                                        patch[0].soil_defaults[0][0].active_zone_z,
                                        patch[0].soil_defaults[0][0].DON_adsorption_rate,
										23,patch);
						patch[0].surface_DON += Nout;
						patch[0].soil_ns.DON -= Nout;
					}
					if (grow_flag > 0) {
						Nout =
								compute_N_leached(verbose_flag,
										patch[0].sat_NO3,//soil_ns.nitrate,
                                        excess,
                                        patch[0].soil_defaults[0][0].NO3decayRate,
                                        patch[0].soil_defaults[0][0].active_zone_z,
                                        patch[0].soil_defaults[0][0].NO3_adsorption_rate,
										17,patch);
						patch[0].surface_NO3 += Nout;
						patch[0].soil_ns.nitrate -= Nout;
					}

					if (grow_flag > 0) {
						Nout =
								compute_N_leached(verbose_flag,
										patch[0].sat_NH4, //soil_ns.sminn,
                                        excess,
                                        patch[0].soil_defaults[0][0].NH4decayRate,
                                        patch[0].soil_defaults[0][0].active_zone_z,
                                        patch[0].soil_defaults[0][0].NH4_adsorption_rate,
										20,patch);
						patch[0].surface_NH4 += Nout;
						patch[0].soil_ns.sminn -= Nout;
					}
				}
				/*--------------------------------------------------------------*/
				/*	final overland flow routing				*/
				/*--------------------------------------------------------------*/
				if (((excess = patch[0].detention_store
						- patch[0].landuse_defaults[0][0].detention_store_size* (1.0 - patch[0].Ksat_vertical)-patch[0].landuse_defaults[0][0].pond_size * patch[0].waterFrac)
						> ZERO) && (patch[0].detention_store > ZERO)) {

					if (patch[0].drainage_type == STREAM) {
						if (grow_flag > 0) {
							patch[0].streamflow_DON += (excess
									/ patch[0].detention_store)
									* patch[0].surface_DON;
							patch[0].streamflow_DOC += (excess
									/ patch[0].detention_store)
									* patch[0].surface_DOC;

							patch[0].streamflow_NO3 += (excess
									/ patch[0].detention_store)
									* patch[0].surface_NO3;
							patch[0].streamNO3_from_surface +=(excess
									/ patch[0].detention_store)
									* patch[0].surface_NO3;
							patch[0].hourly[0].streamflow_NO3 += (excess
									/ patch[0].detention_store)
									* patch[0].surface_NO3;
							patch[0].hourly[0].streamflow_NO3_from_surface +=(excess
									/ patch[0].detention_store)
									* patch[0].surface_NO3;

							patch[0].streamflow_NH4 += (excess
									/ patch[0].detention_store)
									* patch[0].surface_NH4;
							patch[0].surface_DON -= (excess
									/ patch[0].detention_store)
									* patch[0].surface_DON;
							patch[0].surface_DOC -= (excess
									/ patch[0].detention_store)
									* patch[0].surface_DOC;
							patch[0].surface_NO3 -= (excess
									/ patch[0].detention_store)
									* patch[0].surface_NO3;
							patch[0].surface_NH4 -= (excess
									/ patch[0].detention_store)
									* patch[0].surface_NH4;
						}
						patch[0].return_flow += excess;
						patch[0].detention_store -= excess;
						patch[0].Qout_total += excess;
						patch[0].hourly_sur2stream_flow += excess;
					} else {
						/*--------------------------------------------------------------*/
						/* determine which innundation depth to consider		*/
						/*--------------------------------------------------------------*/
						if (patch[0].num_innundation_depths > 0) {
							innundation_depth = patch[0].detention_store;
							d = 0;
							while ((innundation_depth
									> patch[0].innundation_list[d].critical_depth)
									&& (d < patch[0].num_innundation_depths - 1)) {
								d++;
							}
						} else {
							d = 0;
						}

						for (j = 0; j < patch->surface_innundation_list[d].num_neighbours; j++) {
							neigh = patch->surface_innundation_list[d].neighbours[j].patch;
							Qout = excess * patch->surface_innundation_list[d].neighbours[j].gamma;
							if (grow_flag > 0) {
								NO3_out = Qout / patch[0].detention_store
										* patch[0].surface_NO3;
								NH4_out = Qout / patch[0].detention_store
										* patch[0].surface_NH4;
								DON_out = Qout / patch[0].detention_store
										* patch[0].surface_DON;
								DOC_out = Qout / patch[0].detention_store
										* patch[0].surface_DOC;
								Nout = NO3_out + NH4_out + DON_out;
							}
							if (neigh[0].drainage_type == STREAM) {
								neigh[0].Qin_total += Qout * patch[0].area
										/ neigh[0].area;
								neigh[0].return_flow += Qout * patch[0].area
										/ neigh[0].area;
								if (grow_flag > 0) {
									neigh[0].streamflow_DOC += (DOC_out
											* patch[0].area / neigh[0].area);
									neigh[0].streamflow_DON += (DON_out
											* patch[0].area / neigh[0].area);

									neigh[0].streamflow_NO3 += (NO3_out
											* patch[0].area / neigh[0].area);
									neigh[0].streamNO3_from_surface +=(NO3_out
											* patch[0].area / neigh[0].area);
									neigh[0].hourly[0].streamflow_NO3 += (NO3_out
											* patch[0].area / neigh[0].area);
									neigh[0].hourly[0].streamflow_NO3_from_sub +=(NO3_out
											* patch[0].area / neigh[0].area);


									neigh[0].streamflow_NH4 += (NH4_out
											* patch[0].area / neigh[0].area);
									neigh[0].surface_ns_leach += (Nout
											* patch[0].area / neigh[0].area);
								}
							} else {
								neigh[0].Qin_total += Qout * patch[0].area
										/ neigh[0].area;
								neigh[0].detention_store += Qout * patch[0].area
										/ neigh[0].area;
								if (grow_flag > 0) {
									neigh[0].surface_DOC += (DOC_out
											* patch[0].area / neigh[0].area);
									neigh[0].surface_DON += (DON_out
											* patch[0].area / neigh[0].area);
									neigh[0].surface_NO3 += (NO3_out
											* patch[0].area / neigh[0].area);
									neigh[0].surface_ns_leach -= (Nout
											* patch[0].area / neigh[0].area);
									neigh[0].surface_NH4 += (NH4_out
											* patch[0].area / neigh[0].area);

								}

							}
						}
						if (grow_flag > 0) {
							patch[0].surface_DOC -= (excess
									/ patch[0].detention_store)
									* patch[0].surface_DOC;
							patch[0].surface_DON -= (excess
									/ patch[0].detention_store)
									* patch[0].surface_DON;
							patch[0].surface_NO3 -= (excess
									/ patch[0].detention_store)
									* patch[0].surface_NO3;
							patch[0].surface_NH4 -= (excess
									/ patch[0].detention_store)
									* patch[0].surface_NH4;
							patch[0].surface_ns_leach += (excess
									/ patch[0].detention_store)
									* patch[0].surface_NO3;
						}
						patch[0].detention_store -= excess;
						patch[0].Qout_total += excess;
					}
				}
				/*-------------------------------------------------------------------------*/
				/*Recompute current actual depth to water table				*/
				/*-------------------------------------------------------------------------*/
//				patch[0].sat_deficit_z = compute_z_final(verbose_flag,
//						patch[0].soil_defaults[0][0].porosity_0,
//						patch[0].soil_defaults[0][0].porosity_decay,
//						patch[0].soil_defaults[0][0].soil_depth, 0.0,
//						-1.0 * patch[0].sat_deficit);
            
            // need to be careful here: patch[0].sat_deficit could be negative.
            if(patch[0].sat_deficit >= 0){
                patch[0].sat_deficit = min(patch[0].sat_deficit,patch[0].soil_defaults[0][0].soil_water_cap);
                patch[0].available_soil_water = patch[0].soil_defaults[0][0].soil_water_cap - patch[0].sat_deficit;
                patch[0].sat_def_pct = patch[0].sat_deficit * patch[0].soil_defaults[0][0].max_sat_def_1;
                patch[0].sat_def_pct_index = (int)(patch[0].sat_def_pct*1000);
                patch[0].sat_def_pct_indexM = 1000*(patch[0].sat_def_pct - patch[0].sat_def_pct_index*0.001);

                patch[0].sat_deficit_z = patch[0].sat_def_pct_indexM * patch[0].soil_defaults[0][0].sat_def_z[patch[0].sat_def_pct_index+1] + (1.0-patch[0].sat_def_pct_indexM) * patch[0].soil_defaults[0][0].sat_def_z[patch[0].sat_def_pct_index];
            }else{
                // surface
                patch[0].available_soil_water = patch[0].soil_defaults[0][0].soil_water_cap;
                patch[0].sat_deficit_z = patch[0].sat_deficit;
                patch[0].sat_def_pct = 0.0;
                patch[0].sat_def_pct_index = 0;
                patch[0].sat_def_pct_indexM = 0;
            }

            totalfc = patch[0].sat_def_pct_indexM * patch[0].soil_defaults[0][0].fc1_0z[patch[0].sat_def_pct_index+1] + (1.0-patch[0].sat_def_pct_indexM) * patch[0].soil_defaults[0][0].fc1_0z[patch[0].sat_def_pct_index]; // total fc after ET water consumption to SAT
            //totalfc *= (1.0-patch[0].basementFrac); // <---- second thought on this, Oct 8, 2019; basement is 3m at most
            if (patch[0].sat_deficit < ZERO) {
                //patch[0].aboveWT_SatPct = 1.0;
                //patch[0].rootzone.SatPct = 1.0;
                patch[0].rootzone.field_capacity = 0.0;
                patch[0].field_capacity = 0.0;
            } else {
                patch[0].rootzone.field_capacity = totalfc * patch[0].zeroRootCoef * (patch[0].rootzone_scale_ref*patch[0].rootzone_end_reffc[patch[0].sat_def_pct_index] + (1.0-patch[0].rootzone_scale_ref)*patch[0].rootzone_start_reffc[patch[0].sat_def_pct_index]);
                patch[0].rootzone.field_capacity = min(patch[0].rootzone.field_capacity,patch[0].rootzone.potential_sat);
                patch[0].field_capacity = max(0.0,min(patch[0].sat_deficit-patch[0].rootzone.potential_sat, totalfc - patch[0].rootzone.field_capacity));
            }//if else
            
            /*--------------------------------------------------------------*/
				/* 	leave behind field capacity			*/
				/*	if sat deficit has been lowered			*/
				/*	this should be an interactive process, we will use 	*/
				/*	0th order approximation					*/
				/* 	we do not do this once sat def is below 0.9 soil depth	*/
				/*     we use 0.9 to prevent numerical instability		*/
				/*--------------------------------------------------------------*/
				if ((patch[0].sat_deficit_z > patch[0].preday_sat_deficit_z)
						&& (patch[0].sat_deficit_z
								< patch[0].soil_defaults[0][0].soil_depth * 0.9)) {
//					add_field_capacity = compute_layer_field_capacity(
//							command_line[0].verbose_flag,
//							patch[0].soil_defaults[0][0].theta_psi_curve,
//							patch[0].soil_defaults[0][0].psi_air_entry,
//							patch[0].soil_defaults[0][0].pore_size_index,
//							patch[0].soil_defaults[0][0].p3,
//							patch[0].soil_defaults[0][0].p4,
//							patch[0].soil_defaults[0][0].porosity_0,
//							patch[0].soil_defaults[0][0].porosity_decay,
//							patch[0].sat_deficit_z, patch[0].sat_deficit_z,
//							patch[0].preday_sat_deficit_z);
//
//					add_field_capacity = max(add_field_capacity, 0.0);
                    
                    add_field_capacity = max(0.0, totalfc - temp); // only acount for leave-behind water
                    
					patch[0].sat_deficit += add_field_capacity;

					if ((patch[0].sat_deficit_z > patch[0].rootzone.depth)
							&& (patch[0].preday_sat_deficit_z
									> patch[0].rootzone.depth))
						patch[0].unsat_storage += add_field_capacity;
					else
						patch[0].rz_storage += add_field_capacity;
				}				
			
				if (patch[0].rootzone.depth > ZERO) {
					if ((patch[0].sat_deficit > ZERO)
							&& (patch[0].rz_storage == 0.0)) {
//						add_field_capacity = compute_layer_field_capacity(
//								command_line[0].verbose_flag,
//								patch[0].soil_defaults[0][0].theta_psi_curve,
//								patch[0].soil_defaults[0][0].psi_air_entry,
//								patch[0].soil_defaults[0][0].pore_size_index,
//								patch[0].soil_defaults[0][0].p3,
//								patch[0].soil_defaults[0][0].p4,
//								patch[0].soil_defaults[0][0].porosity_0,
//								patch[0].soil_defaults[0][0].porosity_decay,
//								patch[0].sat_deficit_z, patch[0].sat_deficit_z,
//								0.0);
//						add_field_capacity = max(add_field_capacity, 0.0);
						patch[0].sat_deficit += patch[0].rootzone.field_capacity;
                        patch[0].rz_storage += patch[0].rootzone.field_capacity;
                        patch[0].sat_deficit += patch[0].field_capacity;
                        patch[0].unsat_storage += patch[0].field_capacity;
					}
				} else {
					if ((patch[0].sat_deficit > ZERO)
							&& (patch[0].unsat_storage == 0.0)) {
//						add_field_capacity = compute_layer_field_capacity(
//								command_line[0].verbose_flag,
//								patch[0].soil_defaults[0][0].theta_psi_curve,
//								patch[0].soil_defaults[0][0].psi_air_entry,
//								patch[0].soil_defaults[0][0].pore_size_index,
//								patch[0].soil_defaults[0][0].p3,
//								patch[0].soil_defaults[0][0].p4,
//								patch[0].soil_defaults[0][0].porosity_0,
//								patch[0].soil_defaults[0][0].porosity_decay,
//								patch[0].sat_deficit_z, patch[0].sat_deficit_z,
//								0.0);
//						add_field_capacity = max(add_field_capacity, 0.0);
//						patch[0].sat_deficit += add_field_capacity;
//						patch[0].unsat_storage += add_field_capacity;
                        patch[0].sat_deficit += patch[0].field_capacity;
                        patch[0].unsat_storage += patch[0].field_capacity;
					}
				}

				/*--------------------------------------------------------------*/
				/* try to infiltrate this water					*/
				/* use time_int as duration */
				/*--------------------------------------------------------------*/

            // need to be careful here: patch[0].sat_deficit could be negative.
            if(patch[0].sat_deficit >= 0){
                patch[0].sat_deficit = min(patch[0].sat_deficit,patch[0].soil_defaults[0][0].soil_water_cap);
                patch[0].available_soil_water = patch[0].soil_defaults[0][0].soil_water_cap - patch[0].sat_deficit;
                patch[0].sat_def_pct = patch[0].sat_deficit * patch[0].soil_defaults[0][0].max_sat_def_1;
                patch[0].sat_def_pct_index = (int)(patch[0].sat_def_pct*1000);
                patch[0].sat_def_pct_indexM = 1000*(patch[0].sat_def_pct - patch[0].sat_def_pct_index*0.001);
                
                patch[0].sat_deficit_z = patch[0].sat_def_pct_indexM * patch[0].soil_defaults[0][0].sat_def_z[patch[0].sat_def_pct_index+1] + (1.0-patch[0].sat_def_pct_indexM) * patch[0].soil_defaults[0][0].sat_def_z[patch[0].sat_def_pct_index];
            }else{
                // surface
                patch[0].available_soil_water = patch[0].soil_defaults[0][0].soil_water_cap;
                patch[0].sat_deficit_z = patch[0].sat_deficit;
                patch[0].sat_def_pct = 0.0;
                patch[0].sat_def_pct_index = 0;
                patch[0].sat_def_pct_indexM = 0;
            }
            // fc & SatPct
            totalfc = patch[0].sat_def_pct_indexM * patch[0].soil_defaults[0][0].fc1_0z[patch[0].sat_def_pct_index+1] + (1.0-patch[0].sat_def_pct_indexM) * patch[0].soil_defaults[0][0].fc1_0z[patch[0].sat_def_pct_index];
            //totalfc *= (1.0-patch[0].basementFrac); // <---- second thought on this, Oct 8, 2019; basement is 3m at most
            
            if (patch[0].sat_deficit < ZERO) {
                //patch[0].aboveWT_SatPct = 1.0;
                //patch[0].rootzone.SatPct = 1.0;
                patch[0].rootzone.field_capacity = 0.0;
                patch[0].field_capacity = 0.0;
            } else {
                patch[0].rootzone.field_capacity = totalfc * patch[0].zeroRootCoef * (patch[0].rootzone_scale_ref*patch[0].rootzone_end_reffc[patch[0].sat_def_pct_index] + (1.0-patch[0].rootzone_scale_ref)*patch[0].rootzone_start_reffc[patch[0].sat_def_pct_index]);
                patch[0].rootzone.field_capacity = min(patch[0].rootzone.field_capacity,patch[0].rootzone.potential_sat);
                patch[0].field_capacity = max(0.0,min(patch[0].sat_deficit-patch[0].rootzone.potential_sat, totalfc - patch[0].rootzone.field_capacity));
            }//if else
            
            
            if (patch[0].detention_store > ZERO){
                
                infiltration = compute_infiltration(
                    command_line[0].verbose_flag,
                    patch[0].sat_deficit_z,
                    0.0,//patch[0].aboveWT_SatPct, // initiated in daily_I()
                    patch[0].Ksat_vertical, // 1- impervious
                    patch[0].sat_def_pct_indexM * patch[0].soil_defaults[0][0].vksat_0zm[patch[0].sat_def_pct_index+1] + (1.0-patch[0].sat_def_pct_indexM) * patch[0].soil_defaults[0][0].vksat_0zm[patch[0].sat_def_pct_index],
                    patch[0].rz_storage+patch[0].unsat_storage,
                    patch[0].sat_def_pct_indexM * patch[0].soil_defaults[0][0].sat_def_0zm[patch[0].sat_def_pct_index+1] + (1.0-patch[0].sat_def_pct_indexM) * patch[0].soil_defaults[0][0].sat_def_0zm[patch[0].sat_def_pct_index],
                    patch[0].sat_deficit,
                    patch[0].detention_store,
                    time_int,
                    patch[0].soil_defaults[0][0].psi_air_entry);
                
//                if (patch[0].rootzone.depth > ZERO) {
//                    infiltration = compute_infiltration(verbose_flag,
//                            patch[0].sat_deficit_z, patch[0].rootzone.S,
//                            patch[0].Ksat_vertical,
//                            patch[0].soil_defaults[0][0].Ksat_0_v,
//                            patch[0].soil_defaults[0][0].mz_v,
//                            patch[0].soil_defaults[0][0].porosity_0,
//                            patch[0].soil_defaults[0][0].porosity_decay,
//                            (patch[0].detention_store), time_int,
//                            patch[0].soil_defaults[0][0].psi_air_entry);
//                } else {
//                    infiltration = compute_infiltration(verbose_flag,
//                            patch[0].sat_deficit_z, patch[0].S,
//                            patch[0].Ksat_vertical,
//                            patch[0].soil_defaults[0][0].Ksat_0_v,
//                            patch[0].soil_defaults[0][0].mz_v,
//                            patch[0].soil_defaults[0][0].porosity_0,
//                            patch[0].soil_defaults[0][0].porosity_decay,
//                            (patch[0].detention_store), time_int,
//                            patch[0].soil_defaults[0][0].psi_air_entry);
//                }
            }else infiltration = 0.0;
				/*--------------------------------------------------------------*/
				/* added an surface N flux to surface N pool	and		*/
				/* allow infiltration of surface N				*/
				/*--------------------------------------------------------------*/
				if ((grow_flag > 0) && (infiltration > ZERO)) {
					patch[0].soil_ns.DON += ((infiltration
							/ patch[0].detention_store) * patch[0].surface_DON);
					patch[0].soil_cs.DOC += ((infiltration
							/ patch[0].detention_store) * patch[0].surface_DOC);
					patch[0].soil_ns.nitrate += ((infiltration
							/ patch[0].detention_store) * patch[0].surface_NO3);
					patch[0].surface_NO3 -= ((infiltration
							/ patch[0].detention_store) * patch[0].surface_NO3);
					patch[0].soil_ns.sminn += ((infiltration
							/ patch[0].detention_store) * patch[0].surface_NH4);
					patch[0].surface_NH4 -= ((infiltration
							/ patch[0].detention_store) * patch[0].surface_NH4);
					patch[0].surface_DOC -= ((infiltration
							/ patch[0].detention_store) * patch[0].surface_DOC);
					patch[0].surface_DON -= ((infiltration
							/ patch[0].detention_store) * patch[0].surface_DON);
				}
				/*--------------------------------------------------------------*/
				/*	Determine if the infifltration will fill up the unsat	*/
				/*	zone or not.						*/
				/*	We use the strict assumption that sat deficit is the	*/
				/*	amount of water needed to saturate the soil.		*/
				/*--------------------------------------------------------------*/

				if (infiltration
						> patch[0].sat_deficit - patch[0].unsat_storage
								- patch[0].rz_storage) {
					/*--------------------------------------------------------------*/
					/*		Yes the unsat zone will be filled so we may	*/
					/*		as well treat the unsat_storage and infiltration*/
					/*		as water added to the water table.		*/
					/*--------------------------------------------------------------*/
					patch[0].sat_deficit -= (infiltration
							+ patch[0].unsat_storage + patch[0].rz_storage);
					/*--------------------------------------------------------------*/
					/*		There is no unsat_storage left.			*/
					/*--------------------------------------------------------------*/
					patch[0].unsat_storage = 0.0;
					patch[0].rz_storage = 0.0;
					patch[0].field_capacity = 0.0;
					patch[0].rootzone.field_capacity = 0.0;
				} else if ((patch[0].sat_deficit
						> patch[0].rootzone.potential_sat)
						&& (infiltration
								> patch[0].rootzone.potential_sat
										- patch[0].rz_storage)) {
					/*------------------------------------------------------------------------------*/
					/*		Just add the infiltration to the rz_storage and unsat_storage	*/
					/*------------------------------------------------------------------------------*/
					patch[0].unsat_storage += infiltration
							- (patch[0].rootzone.potential_sat
									- patch[0].rz_storage);
					patch[0].rz_storage = patch[0].rootzone.potential_sat;
				}
				/* Only rootzone layer saturated - perched water table case */
				else if ((patch[0].sat_deficit > patch[0].rootzone.potential_sat)
						&& (infiltration
								<= patch[0].rootzone.potential_sat
										- patch[0].rz_storage)) {
					/*--------------------------------------------------------------*/
					/*		Just add the infiltration to the rz_storage	*/
					/*--------------------------------------------------------------*/
					patch[0].rz_storage += infiltration;
				}

				else if ((patch[0].sat_deficit
						<= patch[0].rootzone.potential_sat)
						&& (infiltration
								<= patch[0].sat_deficit - patch[0].rz_storage
										- patch[0].unsat_storage)) {
					patch[0].rz_storage += patch[0].unsat_storage;
					/* transfer left water in unsat storage to rootzone layer */
					patch[0].unsat_storage = 0;
					patch[0].rz_storage += infiltration;
					patch[0].field_capacity = 0;
				}

				if (patch[0].sat_deficit < 0.0) {
					patch[0].detention_store -= (patch[0].sat_deficit
							- patch[0].unsat_storage);
					patch[0].sat_deficit = 0.0;
					patch[0].unsat_storage = 0.0;
				}

				patch[0].detention_store -= infiltration;
				/*--------------------------------------------------------------*/
				/* recompute saturation deficit					*/
				/*--------------------------------------------------------------*/
//				patch[0].sat_deficit_z = compute_z_final(verbose_flag,
//						patch[0].soil_defaults[0][0].porosity_0,
//						patch[0].soil_defaults[0][0].porosity_decay,
//						patch[0].soil_defaults[0][0].soil_depth, 0.0,
//						-1.0 * patch[0].sat_deficit);
//
//
//
//				/*--------------------------------------------------------------*/
//				/*	compute new field capacity				*/
//				/*--------------------------------------------------------------*/
//				if (patch[0].sat_deficit_z < patch[0].rootzone.depth) {
//					patch[0].rootzone.field_capacity =
//							compute_layer_field_capacity(
//									command_line[0].verbose_flag,
//									patch[0].soil_defaults[0][0].theta_psi_curve,
//									patch[0].soil_defaults[0][0].psi_air_entry,
//									patch[0].soil_defaults[0][0].pore_size_index,
//									patch[0].soil_defaults[0][0].p3,
//									patch[0].soil_defaults[0][0].p4,
//									patch[0].soil_defaults[0][0].porosity_0,
//									patch[0].soil_defaults[0][0].porosity_decay,
//									patch[0].sat_deficit_z,
//									patch[0].rootzone.depth, 0.0);
//
//					patch[0].field_capacity = 0.0;
//				} else {
//
//					patch[0].rootzone.field_capacity =
//							compute_layer_field_capacity(
//									command_line[0].verbose_flag,
//									patch[0].soil_defaults[0][0].theta_psi_curve,
//									patch[0].soil_defaults[0][0].psi_air_entry,
//									patch[0].soil_defaults[0][0].pore_size_index,
//									patch[0].soil_defaults[0][0].p3,
//									patch[0].soil_defaults[0][0].p4,
//									patch[0].soil_defaults[0][0].porosity_0,
//									patch[0].soil_defaults[0][0].porosity_decay,
//									patch[0].sat_deficit_z,
//									patch[0].rootzone.depth, 0.0);
//
//					patch[0].field_capacity = compute_layer_field_capacity(
//							command_line[0].verbose_flag,
//							patch[0].soil_defaults[0][0].theta_psi_curve,
//							patch[0].soil_defaults[0][0].psi_air_entry,
//							patch[0].soil_defaults[0][0].pore_size_index,
//							patch[0].soil_defaults[0][0].p3,
//							patch[0].soil_defaults[0][0].p4,
//							patch[0].soil_defaults[0][0].porosity_0,
//							patch[0].soil_defaults[0][0].porosity_decay,
//							patch[0].sat_deficit_z, patch[0].sat_deficit_z, 0)
//							- patch[0].rootzone.field_capacity;
//				}
            
//            /*--------------------------------------------------------------*/
//				/*      Recompute patch soil moisture storage                   */
//				/*--------------------------------------------------------------*/
//				if (patch[0].sat_deficit < ZERO) {
//					patch[0].S = 1.0;
//					patch[0].rootzone.S = 1.0;
//					rz_drainage = 0.0;
//					unsat_drainage = 0.0;
//				} else if (patch[0].sat_deficit_z > patch[0].rootzone.depth) { /* Constant vertical profile of soil porosity */
//
//					/*-------------------------------------------------------*/
//					/*	soil drainage and storage update	     	 */
//					/*-------------------------------------------------------*/
//					patch[0].rootzone.S =
//							min(patch[0].rz_storage / patch[0].rootzone.potential_sat, 1.0);
//					rz_drainage = compute_unsat_zone_drainage(
//							command_line[0].verbose_flag,
//							patch[0].soil_defaults[0][0].theta_psi_curve,
//							patch[0].soil_defaults[0][0].pore_size_index,
//							patch[0].rootzone.S,
//							patch[0].soil_defaults[0][0].mz_v,
//							patch[0].rootzone.depth,
//							patch[0].soil_defaults[0][0].Ksat_0 / n_timesteps
//									/ 2,
//							patch[0].rz_storage
//									- patch[0].rootzone.field_capacity);
//
//					patch[0].rz_storage -= rz_drainage;
//					patch[0].unsat_storage += rz_drainage;
//
//					patch[0].S =
//							min(patch[0].unsat_storage / (patch[0].sat_deficit - patch[0].rootzone.potential_sat), 1.0);
//					unsat_drainage = compute_unsat_zone_drainage(
//							command_line[0].verbose_flag,
//							patch[0].soil_defaults[0][0].theta_psi_curve,
//							patch[0].soil_defaults[0][0].pore_size_index,
//							patch[0].S, patch[0].soil_defaults[0][0].mz_v,
//							patch[0].sat_deficit_z,
//							patch[0].soil_defaults[0][0].Ksat_0 / n_timesteps
//									/ 2,
//							patch[0].unsat_storage - patch[0].field_capacity);
//
//					patch[0].unsat_storage -= unsat_drainage;
//					patch[0].sat_deficit -= unsat_drainage;
//				} else {
//					patch[0].sat_deficit -= patch[0].unsat_storage; /* transfer left water in unsat storage to rootzone layer */
//					patch[0].unsat_storage = 0.0;
//
//					patch[0].S =
//							min(patch[0].rz_storage / patch[0].sat_deficit, 1.0);
//					rz_drainage = compute_unsat_zone_drainage(
//							command_line[0].verbose_flag,
//							patch[0].soil_defaults[0][0].theta_psi_curve,
//							patch[0].soil_defaults[0][0].pore_size_index,
//							patch[0].S, patch[0].soil_defaults[0][0].mz_v,
//							patch[0].sat_deficit_z,
//							patch[0].soil_defaults[0][0].Ksat_0 / n_timesteps
//									/ 2,
//							patch[0].rz_storage
//									- patch[0].rootzone.field_capacity);
//
//					unsat_drainage = 0.0;
//
//					patch[0].rz_storage -= rz_drainage;
//					patch[0].sat_deficit -= rz_drainage;
//				}
            
            // need to be careful here: patch[0].sat_deficit could be negative.
            if(patch[0].sat_deficit >= 0){
                patch[0].sat_deficit = min(patch[0].sat_deficit,patch[0].soil_defaults[0][0].soil_water_cap);
                patch[0].available_soil_water = patch[0].soil_defaults[0][0].soil_water_cap - patch[0].sat_deficit;
                patch[0].sat_def_pct = patch[0].sat_deficit * patch[0].soil_defaults[0][0].max_sat_def_1;
                patch[0].sat_def_pct_index = (int)(patch[0].sat_def_pct*1000);
                patch[0].sat_def_pct_indexM = 1000*(patch[0].sat_def_pct - patch[0].sat_def_pct_index*0.001);
                
                patch[0].sat_deficit_z = patch[0].sat_def_pct_indexM * patch[0].soil_defaults[0][0].sat_def_z[patch[0].sat_def_pct_index+1] + (1.0-patch[0].sat_def_pct_indexM) * patch[0].soil_defaults[0][0].sat_def_z[patch[0].sat_def_pct_index];
            }else{
                // surface
                patch[0].available_soil_water = patch[0].soil_defaults[0][0].soil_water_cap;
                patch[0].sat_deficit_z = patch[0].sat_deficit;
                patch[0].sat_def_pct = 0.0;
                patch[0].sat_def_pct_index = 0;
                patch[0].sat_def_pct_indexM = 0;
            }
            
            // fc & SatPct
            totalfc = patch[0].sat_def_pct_indexM * patch[0].soil_defaults[0][0].fc1_0z[patch[0].sat_def_pct_index+1] + (1.0-patch[0].sat_def_pct_indexM) * patch[0].soil_defaults[0][0].fc1_0z[patch[0].sat_def_pct_index];
            //totalfc *= (1.0-patch[0].basementFrac); // <---- second thought on this, Oct 8, 2019; basement is 3m at most
            
            if (patch[0].sat_deficit < ZERO) {
                 //patch[0].aboveWT_SatPct = 1.0;
                 //patch[0].rootzone.SatPct = 1.0;
                 patch[0].rootzone.field_capacity = 0.0;
                 patch[0].field_capacity = 0.0;
                unsat_drainage = 0.0;
                rz_drainage = 0.0;
             } else if (patch[0].sat_deficit_z > patch[0].rootzone.depth)  {
                 patch[0].rootzone.field_capacity = totalfc * patch[0].zeroRootCoef *  (patch[0].rootzone_scale_ref*patch[0].rootzone_end_reffc[patch[0].sat_def_pct_index] + (1.0-patch[0].rootzone_scale_ref)*patch[0].rootzone_start_reffc[patch[0].sat_def_pct_index]);
                 patch[0].rootzone.field_capacity = min(patch[0].rootzone.field_capacity,patch[0].rootzone.potential_sat);
                 patch[0].field_capacity = max(0.0,min(patch[0].sat_deficit-patch[0].rootzone.potential_sat, totalfc - patch[0].rootzone.field_capacity));
                 
                 //patch[0].rootzone.SatPct = min(patch[0].rz_storage/patch[0].rootzone.potential_sat/(1.0-patch[0].basementFrac), 1.0);
                 //patch[0].aboveWT_SatPct = (patch[0].rz_storage+patch[0].unsat_storage)/patch[0].sat_deficit;
                 
                 /// drainage
                 rz_drainage = compute_unsat_zone_drainage(
                         command_line[0].verbose_flag,
                         patch[0].soil_defaults[0][0].theta_psi_curve,
                         patch[0].soil_defaults[0][0].pore_size_index,
                         patch[0].rootzone.potential_sat, //patch[0].rootzone.SatPct,
                         0.04166667*(patch[0].rootdepth_indexM * patch[0].soil_defaults[0][0].vksat_0zm[patch[0].rootdepth_index+1] + (1.0-patch[0].rootdepth_indexM)* patch[0].soil_defaults[0][0].vksat_0zm[patch[0].rootdepth_index]),
                         0.04166667*(patch[0].rootdepth_indexM * patch[0].soil_defaults[0][0].vksat_z[patch[0].rootdepth_index+1] + (1.0-patch[0].rootdepth_indexM) * patch[0].soil_defaults[0][0].vksat_z[patch[0].rootdepth_index]),
                         patch[0].rz_storage,
                         patch[0].rootzone.field_capacity,
                         patch[0].sat_deficit);
                 patch[0].rz_storage -=  rz_drainage;
                 patch[0].unsat_storage +=  rz_drainage;
                 
                 unsat_drainage = compute_unsat_zone_drainage(
                         command_line[0].verbose_flag,
                         patch[0].soil_defaults[0][0].theta_psi_curve,
                         patch[0].soil_defaults[0][0].pore_size_index,
                         patch[0].sat_deficit - patch[0].rootzone.potential_sat,//patch[0].aboveWT_SatPct,
                         0.04166667*(patch[0].sat_def_pct_indexM * patch[0].soil_defaults[0][0].vksat_0zm[patch[0].sat_def_pct_index+1] + (1.0-patch[0].sat_def_pct_indexM)* patch[0].soil_defaults[0][0].vksat_0zm[patch[0].sat_def_pct_index]),
                         0.04166667*(patch[0].sat_def_pct_indexM * patch[0].soil_defaults[0][0].vksat_z[patch[0].sat_def_pct_index+1] + (1.0-patch[0].sat_def_pct_indexM) * patch[0].soil_defaults[0][0].vksat_z[patch[0].sat_def_pct_index]),
                         patch[0].unsat_storage,
                         patch[0].field_capacity,
                         patch[0].sat_deficit);
                 patch[0].unsat_storage -=  unsat_drainage;
                 patch[0].sat_deficit -=  unsat_drainage;
            
             }//if else
            
            patch[0].unsat_drainage += unsat_drainage;
            patch[0].rz_drainage += rz_drainage;

            // need to be careful here: patch[0].sat_deficit could be negative.
            if(patch[0].sat_deficit >= 0){
                patch[0].sat_deficit = min(patch[0].sat_deficit,patch[0].soil_defaults[0][0].soil_water_cap);
                patch[0].available_soil_water = patch[0].soil_defaults[0][0].soil_water_cap - patch[0].sat_deficit;
                patch[0].sat_def_pct = patch[0].sat_deficit * patch[0].soil_defaults[0][0].max_sat_def_1;
                patch[0].sat_def_pct_index = (int)(patch[0].sat_def_pct*1000);
                patch[0].sat_def_pct_indexM = 1000*(patch[0].sat_def_pct - patch[0].sat_def_pct_index*0.001);
                
                patch[0].sat_deficit_z = patch[0].sat_def_pct_indexM * patch[0].soil_defaults[0][0].sat_def_z[patch[0].sat_def_pct_index+1] + (1.0-patch[0].sat_def_pct_indexM) * patch[0].soil_defaults[0][0].sat_def_z[patch[0].sat_def_pct_index];
            }else{
                // surface
                patch[0].available_soil_water = patch[0].soil_defaults[0][0].soil_water_cap;
                patch[0].sat_deficit_z = patch[0].sat_deficit;
                patch[0].sat_def_pct = 0.0;
                patch[0].sat_def_pct_index = 0;
                patch[0].sat_def_pct_indexM = 0;
            }
            
            // fc & SatPct
            totalfc = patch[0].sat_def_pct_indexM * patch[0].soil_defaults[0][0].fc1_0z[patch[0].sat_def_pct_index+1] + (1.0-patch[0].sat_def_pct_indexM) * patch[0].soil_defaults[0][0].fc1_0z[patch[0].sat_def_pct_index];
            //totalfc *= (1.0-patch[0].basementFrac); // <---- second thought on this, Oct 8, 2019; basement is 3m at most
            
            if (patch[0].sat_deficit < ZERO) {
                //patch[0].aboveWT_SatPct = 1.0;
                //patch[0].rootzone.SatPct = 1.0;
                patch[0].rootzone.field_capacity = 0.0;
                patch[0].field_capacity = 0.0;
            } else {
                patch[0].rootzone.field_capacity = totalfc * patch[0].zeroRootCoef * (patch[0].rootzone_scale_ref*patch[0].rootzone_end_reffc[patch[0].sat_def_pct_index] + (1.0-patch[0].rootzone_scale_ref)*patch[0].rootzone_start_reffc[patch[0].sat_def_pct_index]);
                patch[0].rootzone.field_capacity = min(patch[0].rootzone.field_capacity,patch[0].rootzone.potential_sat);
                patch[0].field_capacity = max(0.0,min(patch[0].sat_deficit-patch[0].rootzone.potential_sat, totalfc - patch[0].rootzone.field_capacity));
            }//if else
            
//				if (patch[0].sat_deficit > patch[0].rootzone.potential_sat)
//					patch[0].rootzone.S =
//							min(patch[0].rz_storage / patch[0].rootzone.potential_sat, 1.0);
//				else
//					patch[0].rootzone.S =
//							min((patch[0].rz_storage + patch[0].rootzone.potential_sat - patch[0].sat_deficit)
//									/ patch[0].rootzone.potential_sat, 1.0);
//
////				/*-------------------c------------------------------------------------------*/
//				/*	Recompute current actual depth to water table				*/
//				/*-------------------------------------------------------------------------*/
//				patch[0].sat_deficit_z = compute_z_final(verbose_flag,
//						patch[0].soil_defaults[0][0].porosity_0,
//						patch[0].soil_defaults[0][0].porosity_decay,
//						patch[0].soil_defaults[0][0].soil_depth, 0.0,
//						-1.0 * patch[0].sat_deficit);



			/*--------------------------------------------------------------*/
			/*--------------------------------------------------------------*/
			/*--------------------------------------------------------------*/
			/*--------------------------------------------------------------*/
			/*--------------------------------------------------------------*/
			/*--------------------------------------------------------------*/			
			patch[0].hourly_stream_flow += patch[0].hourly_subsur2stream_flow
		      				+ patch[0].hourly_sur2stream_flow;

			basin[0].basin_return_flow += (patch[0].return_flow) * patch[0].area;



			/* ******************************** */
			/* accumulate the daily returnflow and baseflow calculated from update_drainage*/
			/* The N calculation has been completed in update_drainage_***.c routing*/
			if (current_date.hour == n_timesteps){
				    if (patch[0].drainage_type == STREAM) {
					patch[0].streamflow += patch[0].return_flow
							+ patch[0].base_flow;
				    }
			}
		    

			
	

		} /* end i */


	basin[0].basin_outflow /= basin[0].basin_area;
	basin[0].preday_basin_rz_storage /= basin[0].basin_area;
	basin[0].preday_basin_unsat_storage /= basin[0].basin_area;
	basin[0].preday_basin_detention_store /= basin[0].basin_area;
	basin[0].preday_basin_sat_deficit /= basin[0].basin_area;
	basin[0].basin_rz_storage /= basin[0].basin_area;
	basin[0].basin_unsat_storage /= basin[0].basin_area;
	basin[0].basin_detention_store /= basin[0].basin_area;
	basin[0].basin_sat_deficit /= basin[0].basin_area;
	water_balance = basin[0].preday_basin_rz_storage + basin[0].preday_basin_unsat_storage
			+ basin[0].preday_basin_detention_store - basin[0].preday_basin_sat_deficit
			- (basin[0].basin_rz_storage + basin[0].basin_unsat_storage + basin[0].basin_detention_store
					- basin[0].basin_sat_deficit) - basin[0].basin_outflow;


	return;

} /*end compute_subsurface_routing_hourly.c*/


