runoff.c

超强的大尺度水文模拟工具

	      
	  /** Brooks & Corey relation for hydraulic conductivity **/

	  if((tmp_moist = moist[lindex] - layer[lindex].evap / (double)dt) 
	     < resid_moist[lindex])
	    tmp_moist = resid_moist[lindex];
	  
	  if(moist[lindex] > resid_moist[lindex]) {
#if LOW_RES_MOIST
	    avg_matric = pow( 10, (soil_con->depth[lindex+1] 
				     * log10(fabs(matric[lindex]))
				     + soil_con->depth[lindex]
				     * log10(fabs(matric[lindex+1])))
				    / (soil_con->depth[lindex] 
				       + soil_con->depth[lindex+1]) );
	    tmp_moist = resid_moist[lindex]
	      + ( soil_con->max_moist[lindex] - resid_moist[lindex] )
	      * pow( ( avg_matric / soil_con->bubble[lindex] ), -1/b[lindex] );
#endif
	    Q12[lindex] 
	      = Ksat[lindex] * pow(((tmp_moist - resid_moist[lindex]) 
				    / (soil_con->max_moist[lindex] 
				       - resid_moist[lindex])),
				   soil_con->expt[lindex]); 
	  }
	  else Q12[lindex] = 0.;
	  last_layer[last_cnt] = lindex;
	}

	/**************************************************
          Solve for Current Soil Layer Moisture, and
          Check Versus Maximum and Minimum Moisture
          Contents.  
	**************************************************/
	
	firstlayer = TRUE;
	last_index = 0;
	for(lindex=0;lindex<options.Nlayer-1;lindex++) {

	  if(lindex==0) dt_runoff = *runoff / (double) dt;
	  else dt_runoff = 0;

	  /* Store moistures for water balance debugging */
#if LINK_DEBUG
	  if(debug.PRT_BALANCE) {
	    if(time_step==0) {
	      if(firstlayer)
		debug.inflow[dist][band][lindex+2] = inflow - dt_runoff;
	      else {
		debug.inflow[dist][band][lindex+2] = inflow;
		debug.outflow[dist][band][lindex+1] = inflow;
	      }
	    }
	    else {
	      if(firstlayer)
		debug.inflow[dist][band][lindex+2]  += inflow - dt_runoff;
	      else {
		debug.inflow[dist][band][lindex+2]  += inflow;
		debug.outflow[dist][band][lindex+1] += inflow;
	      }
	    }
	  }
#endif

	  /* transport moisture for all sublayers **/
#if LINK_DEBUG
	  if(debug.DEBUG || debug.PRT_BALANCE) 
	    last_moist = moist[lindex];
#endif
	      
	  tmp_inflow = 0.;

	  /** Update soil layer moisture content **/
	  moist[lindex] = moist[lindex] + (inflow - dt_runoff) 
	    - (Q12[lindex] + layer[lindex].evap/(double)dt);

	  /** Verify that soil layer moisture is less than maximum **/
	  if((moist[lindex]+ice[lindex]) > max_moist[lindex]) {
	    tmp_inflow = (moist[lindex]+ice[lindex]) - max_moist[lindex];
	    moist[lindex] = max_moist[lindex] - ice[lindex];

	    if(lindex==0) {
	      Q12[lindex] += tmp_inflow;
	      tmp_inflow = 0;
	    }
	    else {
	      tmplayer = lindex;
	      while(tmp_inflow > 0) {
		tmplayer--;
#if LINK_DEBUG
		if(debug.PRT_BALANCE) {
		  /** Update debugging storage terms **/
		  debug.inflow[dist][band][lindex+2]  -= tmp_inflow;
		  debug.outflow[dist][band][lindex+1] -= tmp_inflow;
		}
#endif
		if(tmplayer<0) {
		  /** If top layer saturated, add to top layer **/
		  *runoff += tmp_inflow;
		  tmp_inflow = 0;
		}
		else {
		  /** else add excess soil moisture to next higher layer **/
		  moist[tmplayer] += tmp_inflow;
		  if((moist[tmplayer]+ice[tmplayer]) > max_moist[tmplayer]) {
		    tmp_inflow = ((moist[tmplayer] + ice[tmplayer])
				  - max_moist[tmplayer]);
		    moist[tmplayer] = max_moist[tmplayer] - ice[tmplayer];
		  }
		  else tmp_inflow=0;
		}
	      }
	    } /** end trapped excess moisture **/
	  } /** end check if excess moisture in top layer **/
	      
	  firstlayer=FALSE;
	  
	  /** verify that current layer moisture is greater than minimum **/
	  if ((moist[lindex]+ice[lindex]) < resid_moist[lindex]) {
	    /** moisture cannot fall below residual moisture content **/
	    Q12[lindex] += moist[lindex] - resid_moist[lindex];
	    moist[lindex] = resid_moist[lindex];
	  }
	      
	  inflow = (Q12[lindex]+tmp_inflow);
	  Q12[lindex] += tmp_inflow;
	      
	  last_index++;
	      
	} /* end loop through soil layers */
	  
	/**************************************************
	   Compute Baseflow
        **************************************************/
      
	/** ARNO model for the bottom soil layer (based on bottom
	 soil layer moisture from previous time step) **/
      
	lindex = options.Nlayer-1;
	Dsmax = soil_con->Dsmax / 24.;

#if LINK_DEBUG
	if(debug.DEBUG || debug.PRT_BALANCE) {
	  last_moist = moist[lindex];
	  /** Update debugging storage terms **/
	  debug.outflow[dist][band][lindex+1] += Q12[lindex-1];
	  debug.inflow[dist][band][lindex+2]  += Q12[lindex-1];
	}
#endif
      
	frac = soil_con->Ds * Dsmax 
	  / (soil_con->Ws * soil_con->max_moist[lindex]);
	dt_baseflow = frac * ( moist[lindex] - resid_moist[lindex] );
	if (moist[lindex] > soil_con->Ws * soil_con->max_moist[lindex]) {
	  frac = (moist[lindex] - soil_con->Ws * soil_con->max_moist[lindex]) 
	    / (soil_con->max_moist[lindex] - soil_con->Ws 
	       * soil_con->max_moist[lindex]);
	  dt_baseflow += (Dsmax - soil_con->Ds * Dsmax / soil_con->Ws) 
	    * pow(frac,soil_con->c);
	}

        /** Make sure baseflow isn't negative **/
	if(dt_baseflow < 0) dt_baseflow = 0;

	/** Extract baseflow from the bottom soil layer **/ 
	
	moist[lindex] += Q12[lindex-1] - (layer[lindex].evap/(double)dt 
					  + dt_baseflow);
      
	/** Check Lower Sub-Layer Moistures **/
	tmp_moist = 0;
	
        /* If baseflow > 0 and soil moisture has gone below residual,
         * take water out of baseflow and add back to soil to make up the difference
         * Note: if baseflow is small, soil moisture may still be < residual after this */
	if(dt_baseflow > 0 && (moist[lindex]+ice[lindex]) < resid_moist[lindex]) {
          if ( dt_baseflow > resid_moist[lindex] - (moist[lindex]+ice[lindex]) ) {
            dt_baseflow -= resid_moist[lindex] - (moist[lindex]+ice[lindex]);
            moist[lindex] += resid_moist[lindex] - (moist[lindex]+ice[lindex]);
	  }
          else {
            moist[lindex] += dt_baseflow;
            dt_baseflow = 0.0;
          }
	}

	if((moist[lindex]+ice[lindex]) > max_moist[lindex]) {
	  /* soil moisture above maximum */
	  tmp_moist = ((moist[lindex]+ice[lindex]) - max_moist[lindex]);
	  moist[lindex] = max_moist[lindex] - ice[lindex];
	  tmplayer = lindex;
	  while(tmp_moist > 0) {
	    tmplayer--;
#if LINK_DEBUG
	    if(debug.PRT_BALANCE) {
	      /** Update debugging storage terms **/
	      debug.inflow[dist][band][lindex+2]  -= tmp_moist;
	      debug.outflow[dist][band][lindex+1] -= tmp_moist;
	    }
#endif
	    if(tmplayer<0) {
	      /** If top layer saturated, add to top layer **/
	      *runoff += tmp_moist;
	      tmp_moist = 0;
	    }
	    else {
	      /** else if sublayer exists, add excess soil moisture **/
	      moist[tmplayer] += tmp_moist ;
	      if((moist[tmplayer]+ice[tmplayer]) > max_moist[tmplayer]) {
		tmp_moist = ((moist[tmplayer] + ice[tmplayer])
				 - max_moist[tmplayer]);
		moist[tmplayer] = max_moist[tmplayer] - ice[tmplayer];
	      }
	      else tmp_moist=0;
	    }
	  }
	}

	for(lindex=0;lindex<options.Nlayer;lindex++) 
	  layer[lindex].moist      = moist[lindex] + ice[lindex];      
	*baseflow += dt_baseflow;

      } /* end of hourly time step loop */
      
#if LINK_DEBUG
      if(debug.PRT_BALANCE) {
	debug.outflow[dist][band][options.Nlayer+2] = (*runoff + *baseflow);
	debug.outflow[dist][band][options.Nlayer+1] = *baseflow;
      }
#endif
    }
  } /** Loop over wet and dry fractions **/

  /** Recompute Thermal Parameters Based on New Moisture Distribution **/
  if(options.FULL_ENERGY || options.FROZEN_SOIL) {

    for(lindex=0;lindex<options.Nlayer;lindex++) {
      tmp_layer = find_average_layer(&(layer_wet[lindex]), 
				     &(layer_dry[lindex]), 
				     soil_con->depth[lindex], tmp_mu);
      moist[lindex] = tmp_layer.moist;
    }

    distribute_node_moisture_properties(energy->moist, energy->ice,
					energy->kappa_node, energy->Cs_node,
					soil_con->dz_node, energy->T,
					soil_con->max_moist_node,
#if QUICK_FS
					soil_con->ufwc_table_node,
#else
					soil_con->expt_node,
					soil_con->bubble_node, 
#endif
					moist, soil_con->depth, 
					soil_con->soil_density,
					soil_con->bulk_density,
					soil_con->quartz, Nnodes, 
					options.Nlayer, soil_con->FS_ACTIVE);

  }

}