soil_conduction.c

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

	    and the model exits.				TJB

*********************************************************************/

  extern option_struct options;

  char ErrStr[MAXSTRING];
  char PAST_BOTTOM;
  int nidx, lidx, i;
  int tmplidx;
  double Lsum; /* cumulative depth of moisture layer */
  double Zsum; /* cumulative depth of thermal node */
  double tmp_moist;
  double tmp_T;
  double fract;
  double tmpfract;
  double ice_sum;
  double dz, dz_int;
  double T_upper;
  double T_mid;
  double T_lower;
  double T_int;
  double factor;
  double soil_fract;

  lidx = 0;
  Lsum = 0.;
  Zsum = 0.;
  PAST_BOTTOM = FALSE;

  /* node estimates */
  for(nidx=0;nidx<Nnodes;nidx++) {
 
    if(Zsum == Lsum + depth[lidx] && nidx != 0 && lidx != Nlayers-1) {
      /* node on layer boundary */
      moist_node[nidx] = (moist[lidx] / depth[lidx] 
			      + moist[lidx+1] / depth[lidx+1]) / 1000 / 2.;
      soil_fract = (bulk_density[lidx] / soil_density[lidx] 
		    + bulk_density[lidx+1] / soil_density[lidx+1]) / 2.;
    }
    else { 
      /* node completely in layer */
      moist_node[nidx] = moist[lidx] / depth[lidx] / 1000;
      soil_fract = (bulk_density[lidx] / soil_density[lidx]);
    }      

    // Check that node moisture does not exceed maximum node moisture
    if (abs(moist_node[nidx]-max_moist_node[nidx]) > SMALL) {
      sprintf( ErrStr, "Node soil moisture, %f, exceeds maximum node soil moisuttre, %f.", moist_node[nidx], max_moist_node[nidx] );
      vicerror(ErrStr);
    }

    if(T_node[nidx] < 0 && (FS_ACTIVE && options.FROZEN_SOIL)) {
      /* compute moisture and ice contents */
#if QUICK_FS
      ice_node[nidx] 
	= moist_node[nidx] - maximum_unfrozen_water_quick(T_node[nidx],
						   max_moist_node[nidx], 
						   ufwc_table_node[nidx]);
#else
      ice_node[nidx] 
	= moist_node[nidx] - maximum_unfrozen_water(T_node[nidx],
					     max_moist_node[nidx], 
					     bubble_node[nidx],
					     expt_node[nidx]);
#endif
      if(ice_node[nidx]<0) ice_node[nidx]=0;

      /* compute thermal conductivity */
      kappa_node[nidx] 
	= soil_conductivity(moist_node[nidx], moist_node[nidx] 
			    - ice_node[nidx], soil_density[lidx],
			    bulk_density[lidx], quartz[lidx]);

    }
    else {
      /* compute moisture and ice contents */
      ice_node[nidx]   = 0;
      /* compute thermal conductivity */
      kappa_node[nidx] 
	= soil_conductivity(moist_node[nidx], moist_node[nidx], 
			    soil_density[lidx], bulk_density[lidx], 
			    quartz[lidx]);
    }
    /* compute volumetric heat capacity */
    Cs_node[nidx] = volumetric_heat_capacity(bulk_density[lidx] 
					     / soil_density[lidx],
					     moist_node[nidx] - ice_node[nidx],
					     ice_node[nidx]);

    Zsum += (dz_node[nidx] + dz_node[nidx+1]) / 2.;
    if(Zsum > Lsum + depth[lidx] && !PAST_BOTTOM) {
      Lsum += depth[lidx];
      lidx++;
      if( lidx == Nlayers ) {
	PAST_BOTTOM = TRUE;
	lidx = Nlayers-1;
      }
    }
  }

}

#undef N_INTS

void estimate_layer_ice_content(layer_data_struct *layer,
				double            *dz,
				double            *T,
				double            *max_moist_node,
#if QUICK_FS
				double          ***ufwc_table_node,
#else
				double            *expt_node,
				double            *bubble_node,
#endif
				double            *depth,
				double            *max_moist,
#if QUICK_FS
				double          ***ufwc_table_layer,
#else
				double            *expt,
				double            *bubble,
#endif
				double            *bulk_density,
				double            *soil_density,
				double            *quartz,
				float            **layer_node_fract,
				int                Nnodes, 
				int                Nlayers,
				char               FS_ACTIVE) {
/**************************************************************
  This subroutine estimates the ice content of all soil 
  moisture layers based on the distribution of soil thermal
  node temperatures.

  layer_struct *layer           structure with all soil moisture layer info
  double       *dz              soil thermal node thicknesses (m)
  double       *T               soil thermal node temperatures (C)
  double       *max_moist_node  soil thermal node max moisture content (mm/mm)
  double       *expt_node       soil thermal node exponential ()
  double       *bubble_node     soil thermal node bubbling pressure (cm)
  double       *depth           soil moisture layer thickness (m)
  double       *max_moist       soil layer maximum soil moisture (mm)
  double       *expt            soil layer exponential ()
  double       *bubble          soil layer bubling pressure (cm)
  double       *bulk_density    soil layer bulk density (kg m-3)
  double        soil_density    soil layer soil density (kg m-3)
  double        quartz          soil layer quartz content (fract)
  int           Nnodes          number of soil thermal nodes
  int           Nlayer          number of soil moisture layers

**************************************************************/

  extern option_struct options;

  int    nidx;
  int    lidx;
  double Lsum; /* cumulative depth of moisture layer */
  double Zsum; /* cumulative depth of thermal node */
  double deltaz;
  double fract; /* fract of internodal region in layer */
  double boundT; /* soil temperature between layers */
  double tmp_ice;
  double ice_content[2]; /* stores estimated nodal ice content */
  double kappa_layer[2]; /* stores node thermal conductivity */
  double Cs_layer[2]; /* stores node volumetric heat capacity */
  double T_sum; /* summation of nodal temperatures within the layer */
  double ice_sum; /* summation of nodal ice contents within the layer */
  double dz_sum; /* summation of depth used in computing statistics */

  for(lidx=0;lidx<Nlayers;lidx++) {
    layer[lidx].T = 0.;
    layer[lidx].ice = 0.;
    for(nidx=0;nidx<Nnodes;nidx++) {
      if(layer_node_fract[lidx][nidx] > 0) {
	layer[lidx].T += T[nidx] * layer_node_fract[lidx][nidx] * dz[nidx];
	if(T[nidx] < 0 && options.FROZEN_SOIL && FS_ACTIVE) {
	  tmp_ice = layer[lidx].moist 
#if QUICK_FS
		     - maximum_unfrozen_water_quick(T[nidx], 
						    max_moist[lidx], 
						    ufwc_table_layer[lidx]);
#else
	- maximum_unfrozen_water(T[nidx], max_moist[lidx], bubble[lidx], 
				 expt[lidx]);
#endif
	  if(tmp_ice < 0) tmp_ice = 0;
	}
	else tmp_ice = 0;
	layer[lidx].ice += tmp_ice * layer_node_fract[lidx][nidx] * dz[nidx];
      }
    }
    layer[lidx].T /= depth[lidx];
    layer[lidx].ice /= depth[lidx];
  }

}

void compute_soil_layer_thermal_properties(layer_data_struct *layer,
					   double            *depth,
					   double            *bulk_density,
					   double            *soil_density,
					   double            *quartz,
					   int                Nlayers) {
/********************************************************************
  This subroutine computes the thermal conductivity and volumetric
  heat capacity of each soil layer based on its current moisture
  and ice contents.  Ice is only present if the frozen soil
  algorithm is activated.

  layer_data_struct *layer          structure with all soil layer variables
  double            *depth          soil layer depths (m)
  double            *bulk_density   soil layer bulk density (kg/m^3)
  double             soil_density   soil layer soil density (kg/m^3)
  double             quartz         soil layer quartz content (fract)
  int                Nlayers        number of soil layers

********************************************************************/

  int lidx;
  double moist, ice;

  /* compute layer thermal properties */
  for(lidx=0;lidx<Nlayers;lidx++) {
    moist = layer[lidx].moist / depth[lidx] / 1000;
    ice = layer[lidx].ice / depth[lidx] / 1000;
    layer[lidx].kappa 
      = soil_conductivity(moist, moist - ice, soil_density[lidx], 
			  bulk_density[lidx], quartz[lidx]);
    layer[lidx].Cs 
      = volumetric_heat_capacity(bulk_density[lidx] / soil_density[lidx], 
				 moist - ice, ice);
  }
}

void find_0_degree_fronts(energy_bal_struct *energy,
			  double            *dz,
			  double            *T,
			  int                Nnodes) {
/***********************************************************************
  This subroutine reads through the soil thermal nodes and determines 
  the depths of all thawing and freezing fronts that are present.

  energy_bal_struct *energy  energy balance variable structure
  double            *dz      thermal node thicknesses (m)
  double            *T       thermal node temperatures (C)
  int                Nnodes  number of defined thermal nodes

***********************************************************************/

  int    nidx, fidx; 
  int    Nthaw; /* number of thawing fronts found */
  int    Nfrost; /* number of frost fronts found */
  double tdepth[MAX_FRONTS]; /* thawing frost depths */
  double fdepth[MAX_FRONTS]; /* freezing front depths */
  double Zsum;
  double deltaz;
  
  /* Initialize parameters */
  Zsum = 0;
  for(nidx=0;nidx<Nnodes-1;nidx++)
    Zsum += (dz[nidx] + dz[nidx+1]) / 2.;
  Nthaw = Nfrost = 0;
  for(fidx=0;fidx<MAX_FRONTS;fidx++) {
    fdepth[fidx] = MISSING;
    tdepth[fidx] = MISSING;
  }

  /* find 0 degree fronts */
  for(nidx=Nnodes-2;nidx>=0;nidx--) {
    deltaz = (dz[nidx] + dz[nidx+1]) / 2.;
    if(T[nidx] > 0 && T[nidx+1] <= 0 && Nthaw<MAX_FRONTS) {
      tdepth[Nthaw] = linear_interp(0,T[nidx],T[nidx+1],Zsum-deltaz,Zsum);
      Nthaw++;
    }
    else if(T[nidx] < 0 && T[nidx+1] >= 0 && Nfrost<MAX_FRONTS) {
      fdepth[Nfrost] = linear_interp(0,T[nidx],T[nidx+1],Zsum-deltaz,Zsum);
      Nfrost++;
    }
    Zsum -= deltaz;
  }

  /* store thaw depths */
  for(fidx=0;fidx<MAX_FRONTS;fidx++) energy->tdepth[fidx] = tdepth[fidx];
  /* store frost depths */
  for(fidx=0;fidx<MAX_FRONTS;fidx++) energy->fdepth[fidx] = fdepth[fidx];
  energy->Nthaw = Nthaw;
  energy->Nfrost = Nfrost;

}

double maximum_unfrozen_water(double T,
                              double max_moist,
                              double bubble,
                              double expt) {
/**********************************************************************
  This subroutine computes the maximum amount of unfrozen water that
  can exist at the current temperature.
**********************************************************************/

  double unfrozen;

  unfrozen = max_moist * pow((-Lf * T) / (T
      + 273.16) / (9.18 * bubble / 100.), -(2.0 / (expt - 3.0)));
  if(unfrozen > max_moist) unfrozen = max_moist;
  if(unfrozen < 0) unfrozen = 0;

  return (unfrozen);
}

#if QUICK_FS
double maximum_unfrozen_water_quick(double   T,
				    double   max_moist,
				    double **table) {
/**********************************************************************
  This subroutine computes the maximum amount of unfrozen water that
  can exist at the current temperature.
**********************************************************************/

  extern double temps[];

  int i;
  double unfrozen;

  i = 1;
  while(T < temps[i] && i < QUICK_FS_TEMPS) i++;
  unfrozen = max_moist * (table[i-1][0] + table[i-1][1] * T);
  if(unfrozen > max_moist) unfrozen = max_moist;
  if(unfrozen < 0) unfrozen = 0;

  return (unfrozen);
}
#endif

layer_data_struct find_average_layer(layer_data_struct *wet,
				     layer_data_struct *dry,
				     double             depth,
				     double             mu) {
/*************************************************************
  This subroutine computes the average soil layer moistures
  between the wet and dry fraction for use in computing 
  energy balance parameters.  Other layer variables are copied 
  from the wet fraction structure since they are they same for 
  wet and dry fractions.
**************************************************************/

  layer_data_struct layer;

  layer = *wet;

  layer.ice = ((wet->ice * mu) + (dry->ice * (1. - mu)));
  layer.moist = ((wet->moist * mu) + (dry->moist * (1. - mu)));

  return(layer);

}