func_surf_energy_bal.c

能量分析的完美工具

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <vicNl.h>

static char vcid[] = "$Id: func_surf_energy_bal.c,v 4.2.2.2 2004/06/23 18:37:22 tbohn Exp $";

double func_surf_energy_bal(double Ts, va_list ap)
/**********************************************************************
	func_surf_energy_bal	Keith Cherkauer		January 3, 1996

  This subroutine computes the surface energy balance for bare soil
  and vegetation uncovered by snow.  It computes outgoing longwave,
  sensible heat flux, ground heat flux, and storage of heat in the thin
  upper layer, based on the given surface temperature.

  The Energy Balance Equation used comes from Xu Liang's Paper 
  "Insights of the Ground Heat Flux in Land Surface Parameterization
  Schemes."

  Modifications:
  04-14-98 modified to compute evapotranspiration within this routine
           in the hopes of reducing the number of iteration 
	  needed to find a solution surface temperature.       KAC
  07-13-98 modified to include elevation bands for vegetation 
           and snow                                             KAC
  01-20-00 modified to work with the updated radiation estimation
           routines, as well as the simplified frozen soil moisture
           storage                                              KAC
  07-May-04 Added check that both FS_ACTIVE and FROZEN_SOIL are true
	    before adjusting *deltaH for ice.  This is just a safety
	    measure; ice and ice0 should both be 0 if FS_ACTIVE is
	    FALSE.						TJB

**********************************************************************/
{
  extern option_struct options;
#if LINK_DEBUG
  extern debug_struct  debug;
#endif

  /** Thermal Properties **/
  double             T2;       	/** average soil temperature (C) **/
  double             Ts_old;	/** last temperature (C) **/
  double             T1_old;	/** last layer 1 soil temperature (C) **/
  double             Tair;     	/** Air Temperature **/
  double             ra;       	/** aerodynamic reisistance (s/m) **/
  double             atmos_density;	/** atmospheric density (kg/m^3) **/
  double             shortwave;
  double             longwave;
  double             albedo;
  double             emissivity;
  double             kappa1;	/** thermal conductivity of 1st layer */
  double             kappa2;	/** thermal conductivity of 2nd layer */
  double             Cs1;      	/** volumetric heat capacity of 1st layer **/
  double             Cs2;      	/** volumetric heat capacity of 2nd layer **/
  double             D1;       	/** thickness of 1st layer (m) **/
  double             D2;       	/** thickness of 2nd layer (m) **/
  double             dp;       	/** depth to constant temperature (m) */
  double             delta_t;	/** Time Step in Seconds **/
  double             Le;       	/** Latent heat of vapoization **/
  double             Ls;       	/** Latent heat of sublimation **/
  double             Vapor;    	/** Total vapor flux from snow in m/s **/
  double             moist;    	/** layer moisture content in m/m **/
  double             ice0;     	/** layer ice content in m/m **/
  double             max_moist;	/** layer maximum moisture content in m/m **/
  double             bubble;	/** bubbling pressure in cm **/
  double             expt;
  double             snow_depth;
  double             snow_density;
  double             Tsnow_surf;
  double             snow_cover_fraction;
  double             surf_atten;
  double             wind;
  double             displacement;
  double             roughness;
  double             ref_height;
  double             elevation;
  double             b_infilt;
  double             max_infil;
  double             dt;
  double             vpd;
  double             snow_energy;
  double             mu;
  double            *rainfall;
  double            *Wdew;
  double            *grnd_flux;
  double            *T1;
  double            *latent_heat;
  double            *sensible_heat;
  double            *deltaH;
  double            *snow_flux;
  double            *store_error;
  double             TMean;
  double             rad;
  double            *depth;
  double            *Wcr;
  double            *Wpwp;
  double            *resid_moist;
  double            *T_node;
  double            *Tnew_node;
  double            *dz_node;
  double            *kappa_node;
  double            *Cs_node;
  double            *moist_node;
  double            *bubble_node;
  double            *expt_node;
  double            *max_moist_node;
  double            *ice_node;
  double            *alpha;
  double            *beta;
  double            *gamma;
#if QUICK_FS
  double           **ufwc_table_layer;
  double          ***ufwc_table_node;
#endif
  float             *root;
  layer_data_struct *layer_wet;
  layer_data_struct *layer_dry;
  veg_var_struct    *veg_var_wet;
  veg_var_struct    *veg_var_dry;
  int                VEG;
  int                veg_class;
  int                month;
  int                Nnodes;
  char              *FIRST_SOLN;
  int                SNOWING;
  int                FS_ACTIVE;

  double             error;
  double             ice;
  double             Evap;		/** Total evap in m/s **/
  double             kappa_snow;

  T2                  = (double) va_arg(ap, double);
  Ts_old              = (double) va_arg(ap, double);
  T1_old              = (double) va_arg(ap, double);
  Tair                = (double) va_arg(ap, double);
  ra                  = (double) va_arg(ap, double);
  atmos_density       = (double) va_arg(ap, double);
  shortwave           = (double) va_arg(ap, double);
  longwave            = (double) va_arg(ap, double);
  albedo              = (double) va_arg(ap, double);
  emissivity          = (double) va_arg(ap, double);
  kappa1              = (double) va_arg(ap, double);
  kappa2              = (double) va_arg(ap, double);
  Cs1                 = (double) va_arg(ap, double);
  Cs2                 = (double) va_arg(ap, double);
  D1                  = (double) va_arg(ap, double);
  D2                  = (double) va_arg(ap, double);
  dp                  = (double) va_arg(ap, double);
  delta_t             = (double) va_arg(ap, double);
  Le                  = (double) va_arg(ap, double);
  Ls                  = (double) va_arg(ap, double);
  Vapor               = (double) va_arg(ap, double);
  moist               = (double) va_arg(ap, double);
  ice0                = (double) va_arg(ap, double);
  max_moist           = (double) va_arg(ap, double);
  bubble              = (double) va_arg(ap, double);
  expt                = (double) va_arg(ap, double);
  snow_depth          = (double) va_arg(ap, double);
  snow_density        = (double) va_arg(ap, double);
  Tsnow_surf          = (double) va_arg(ap, double);
  snow_cover_fraction = (double) va_arg(ap, double);
  surf_atten          = (double) va_arg(ap, double);
  wind                = (double) va_arg(ap, double);
  displacement        = (double) va_arg(ap, double);
  roughness           = (double) va_arg(ap, double);
  ref_height          = (double) va_arg(ap, double);
  elevation           = (double) va_arg(ap, double);
  b_infilt            = (double) va_arg(ap, double);
  max_infil           = (double) va_arg(ap, double);
  dt                  = (double) va_arg(ap, double);
  vpd                 = (double) va_arg(ap, double);
  snow_energy         = (double) va_arg(ap, double);
  mu                  = (double) va_arg(ap, double);
  rainfall            = (double *) va_arg(ap, double *);
  Wdew                = (double *) va_arg(ap, double *);
  grnd_flux           = (double *) va_arg(ap, double *);
  T1                  = (double *) va_arg(ap, double *);
  latent_heat         = (double *) va_arg(ap, double *);