func_surf_energy_bal.c

能量分析的完美工具

  sensible_heat       = (double *) va_arg(ap, double *);
  deltaH              = (double *) va_arg(ap, double *);
  snow_flux           = (double *) va_arg(ap, double *);
  store_error         = (double *) va_arg(ap, double *);
  depth               = (double *) va_arg(ap, double *);
  Wcr                 = (double *) va_arg(ap, double *);
  Wpwp                = (double *) va_arg(ap, double *);
  resid_moist         = (double *) va_arg(ap, double *);
  T_node              = (double *) va_arg(ap, double *);
  Tnew_node           = (double *) va_arg(ap, double *);
  dz_node             = (double *) va_arg(ap, double *);
  kappa_node          = (double *) va_arg(ap, double *);
  Cs_node             = (double *) va_arg(ap, double *);
  moist_node          = (double *) va_arg(ap, double *);
  bubble_node         = (double *) va_arg(ap, double *);
  expt_node           = (double *) va_arg(ap, double *);
  max_moist_node      = (double *) va_arg(ap, double *);
  ice_node            = (double *) va_arg(ap, double *);
  alpha               = (double *) va_arg(ap, double *);
  beta                = (double *) va_arg(ap, double *);
  gamma               = (double *) va_arg(ap, double *);
#if QUICK_FS
  ufwc_table_layer    = (double **) va_arg(ap, double **);
  ufwc_table_node     = (double ***) va_arg(ap, double ***);
#endif
  root                = (float  *) va_arg(ap, float  *);
  layer_wet           = (layer_data_struct *) va_arg(ap, layer_data_struct *);
  layer_dry           = (layer_data_struct *) va_arg(ap, layer_data_struct *);
  veg_var_wet         = (veg_var_struct *) va_arg(ap, veg_var_struct *);
  veg_var_dry         = (veg_var_struct *) va_arg(ap, veg_var_struct *);
  VEG                 = (int) va_arg(ap, int);
  veg_class           = (int) va_arg(ap, int);
  month               = (int) va_arg(ap, int);
  Nnodes              = (int) va_arg(ap, int);
  FIRST_SOLN          = (char *)va_arg(ap, char *);
  SNOWING             = (int) va_arg(ap, int);
  FS_ACTIVE           = (int) va_arg(ap, int);

  TMean = Ts;

  if(options.GRND_FLUX) {
  
    /**********************************************
      Compute Surface Temperature at Half Time Step
    **********************************************/
    if(snow_cover_fraction > 0) {

      /****************************************
        Compute energy flux through snow pack
      ****************************************/

      kappa_snow = 2.9302e-6 * (snow_density) * (snow_density);
 
      *snow_flux = kappa_snow * (TMean - Tsnow_surf) / snow_depth;

    }

    if(options.QUICK_FLUX) {
      /**************************************************************
        Use Liang et al. 1999 Equations to Calculate Ground Heat 
	Flux
      **************************************************************/
      *T1 = estimate_T1(TMean, T1_old, T2, D1, D2, kappa1, kappa2, Cs1, 
			Cs2, dp, delta_t);
    
    }
    else {
    /*************************************************************
      Use Finite Difference Method to Explicitly Solve Ground Heat
      Flux at Soil Thermal Nodes (Cherkauer and Lettenmaier, 1999)
    *************************************************************/
      T_node[0] = TMean;
#if QUICK_FS
      solve_T_profile(Tnew_node, T_node, dz_node, kappa_node, Cs_node, 
		      moist_node, delta_t, max_moist_node, bubble_node, 
		      expt_node, ice_node, alpha, beta, gamma, 
		      ufwc_table_node, Nnodes, FIRST_SOLN, FALSE, FS_ACTIVE);
#else
      solve_T_profile(Tnew_node, T_node, dz_node, kappa_node, Cs_node, 
		      moist_node, delta_t, max_moist_node, bubble_node, 
		      expt_node, ice_node, alpha, beta, gamma, Nnodes, 
		      FIRST_SOLN, FALSE, FS_ACTIVE);
#endif
      *T1 = Tnew_node[1];
    }

    /*****************************************************
      Compute the Ground Heat Flux from the Top Soil Layer
    *****************************************************/
/*     *grnd_flux = surf_atten * (kappa1/D1*((*T1) - (TMean))); */
/*     *grnd_flux = (kappa1/D1*((*T1) - (TMean))); */
    *grnd_flux = (snow_cover_fraction + (1. - snow_cover_fraction) 
		  * surf_atten) * (kappa1 / D1 * ((*T1) - TMean));

    /******************************************************
      Compute the Current Ice Content of the Top Soil Layer
    ******************************************************/
    if((FS_ACTIVE && options.FROZEN_SOIL) && (TMean+ *T1)/2.<0.) {
      ice = moist - maximum_unfrozen_water((TMean+ *T1)/2.,
					   max_moist,bubble,expt);
      if(ice<0.) ice=0.;
    }
    else ice=0.;
 
    *deltaH = Cs1 * ((Ts_old + T1_old)/2. - (TMean + *T1)/2.) * D1 / delta_t;
    /* Only adjust for ice if both FS_ACTIVE and FROZEN_SOIL are true */
    if(FS_ACTIVE && options.FROZEN_SOIL)
      *deltaH -= ice_density*Lf*(ice0-ice)*D1/delta_t;

    /** Compute net surface radiation for evaporation estimates **/
    rad = (1.0 - albedo) * shortwave + longwave 
      - STEFAN_B * (TMean+KELVIN) * (TMean+KELVIN) * (TMean+KELVIN) 
      * (TMean+KELVIN) + *grnd_flux + *deltaH;

  } /* End computation for ground heat flux */
  else { /* ground heat flux not estimated */

    if(dt < 24)

      /** Compute net surface radiation for evaporation estimates **/
      rad = (1.0 - albedo) * shortwave + longwave 
	- STEFAN_B * (TMean+KELVIN) * (TMean+KELVIN) * (TMean+KELVIN) 
	* (TMean+KELVIN);

    else

      /** Daily water balance model provides average shortwave and 
	  net longwave radiation **/
      rad = (1.0 - albedo) * shortwave + longwave;

  }

  /*************************************************
    Compute Evapotranspiration if not snow covered
  *************************************************/
  if(VEG && !SNOWING) 
    Evap = canopy_evap(layer_wet, layer_dry, veg_var_wet, veg_var_dry, TRUE, 
		       veg_class, month, mu, Wdew, dt, rad, vpd, 
		       (1.0 - albedo) * shortwave, Tair, ra, displacement, 
		       roughness, ref_height, elevation, rainfall, depth, 
		       Wcr, Wpwp, root);
  else if(!SNOWING)
    Evap = arno_evap(layer_wet, layer_dry, rad, Tair, vpd, 
		     (1.0 - albedo) * shortwave, D1, 
		     max_moist * depth[0] * 1000., elevation, b_infilt,  
		     Tair, displacement, roughness, ref_height, ra, dt, mu,
		     resid_moist[0]);
  else Evap = 0.;
  
  /**********************************************************************
    Compute the Latent Heat Flux from the Surface and Covering Vegetation
  **********************************************************************/
  *latent_heat  = -RHO_W*Le*Evap;
  *latent_heat += -atmos_density*Ls*Vapor;

  if(options.GRND_FLUX) {
  
    /************************************************
      Compute the Sensible Heat Flux from the Surface
    ************************************************/
    *sensible_heat = atmos_density*Cp*(Tair - (TMean))/ra;

    /*************************************
      Compute Surface Energy Balance Error
    *************************************/
    error = (1. - snow_cover_fraction) 
      * ((1.-albedo)*shortwave 
	 + emissivity*(longwave-STEFAN_B*(TMean + KELVIN)*(TMean + KELVIN)
		       *(TMean + KELVIN)*(TMean + KELVIN))
	 + *sensible_heat + *latent_heat + snow_energy) 
      - snow_cover_fraction * *snow_flux + *grnd_flux + *deltaH;
    
    *store_error = error;
  }
  else error = MISSING;

  return error;

}