calc_surf_energy_bal.c

能量分析的完美工具

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

static char vcid[] = "$Id: calc_surf_energy_bal.c,v 4.2.2.2 2004/09/22 00:40:32 vicadmin Exp $";

double calc_surf_energy_bal(int                rec,
			    int                iveg,
			    int                nlayer,
			    int                Nveg,
			    int                dt,
			    int                Nnodes,
			    int                veg_class,
			    int                band,
			    int                hour,
			    double             ice0,
			    double             moist,
			    double             dp,
			    double             surf_atten,
			    double             T0,
			    double             shortwave,
			    double             longwave,
			    double             air_temp,
			    double             Le,
			    double             Ls,
			    double             mu,
			    double             displacement,
			    double             roughness,
			    double             ref_height,
			    double             snow_energy,
			    double             vapor_flux,
			    double             bare_albedo,
			    double            *aero_resist,
			    double            *wind,
			    double            *rainfall,
			    double            *ppt,
			    float             *root,
			    atmos_data_struct *atmos,
			    veg_var_struct    *veg_var_wet,
			    veg_var_struct    *veg_var_dry,
			    energy_bal_struct *energy,
			    snow_data_struct  *snow,
			    layer_data_struct *layer_wet,
			    layer_data_struct *layer_dry,
			    soil_con_struct   *soil_con,
			    dmy_struct        *dmy)
/**************************************************************
  calc_surf_energy_bal.c  Greg O'Donnell and Keith Cherkauer  Sept 9 1997
  
  This function calculates the surface temperature, in the
  case of no snow cover.  Evaporation is computed using the
  previous ground heat flux, and then used to comput latent heat
  in the energy balance routine.  Surface temperature is found
  using the Root Brent method (Numerical Recipies).
  

  modifications:
    02-29-00  Included variables needed to compute energy flux
              through the snow pack.  The ground surface energy
              balance will now be a mixture of snow covered
	      and bare ground, controlled by the snow cover 
	      fraction set in solve_snow.c                 KAC
    04-Jun-04 Placed "ERROR" at beginning of screen dump in
	      error_print_surf_energy_bal.		TJB
    21-Sep-04 Added ErrorString to store error messages from
	      root_brent.				TJB

***************************************************************/
{
  extern veg_lib_struct *veg_lib;
  extern option_struct   options;

  char     FIRST_SOLN[1];
  double   T2;
  double   Ts_old;
  double   T1_old;
  double   Tair;
  double   ra;
  double   atmos_density;
  double   albedo;
  double   emissivity;
  double   kappa1;
  double   kappa2;
  double   Cs1;
  double   Cs2;
  double   D1;
  double   D2;
  double   delta_t;
  double   Vapor;
  double   max_moist;
  double   bubble;
  double   expt;
  double   T1;
  double   snow_depth;
  double   snow_density;
  double   Tsnow_surf;
  double   snow_cover_fraction;
  double   snow_flux;

  int      VEG;
  double   Tsurf;
  double   U;
  double   error;
  double   Wdew[2];
  double  *T_node;
  double  Tnew_node[MAX_NODES];
  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;
  layer_data_struct layer[MAX_LAYERS];

  double   T_lower, T_upper;
  char    ErrorString[MAXSTRING];

  /**************************************************
    Correct Aerodynamic Resistance for Stability
  **************************************************/
  ra = aero_resist[0];
  U = wind[0];
  if (U > 0.0)
    ra /= StabilityCorrection(ref_height, displacement, T0,
          air_temp, U, roughness);
  else
    ra = HUGE_RESIST;
  
  /**************************************************
    Compute Evaporation and Transpiration 
  **************************************************/
  if(iveg!=Nveg) {
    if(veg_lib[veg_class].LAI[dmy->month-1] > 0.0) VEG = TRUE;
    else VEG = FALSE;
  }
  else VEG = FALSE;
  Vapor = vapor_flux / (double)dt / 3600.;

  /**************************************************
    Set All Variables For Use
  **************************************************/

  T2                  = energy->T[Nnodes-1];
  Ts_old              = energy->T[0];
  T1_old              = energy->T[1];
  Tair                = air_temp;
  atmos_density       = atmos->density[hour];
  albedo              = bare_albedo;
  emissivity          = 1.;
  kappa1              = energy->kappa[0];
  kappa2              = energy->kappa[1];
  Cs1                 = energy->Cs[0];
  Cs2                 = energy->Cs[1];
  D1                  = soil_con->depth[0];
/*   D2                  = soil_con->depth[1]; */
  D2                  = soil_con->depth[0];
  delta_t             = (double)dt*3600.;
  max_moist           = soil_con->max_moist[0]/(soil_con->depth[0]*1000.);
  bubble              = soil_con->bubble[0];
  expt                = soil_con->expt[0];
  snow_depth          = snow->depth;
  snow_density        = snow->density;
  Tsnow_surf          = snow->surf_temp;
  snow_cover_fraction = snow->coverage;
  Wdew[WET]           = veg_var_wet->Wdew;
  if(options.DIST_PRCP) Wdew[DRY] = veg_var_dry->Wdew;
  FIRST_SOLN[0] = TRUE;

  /*************************************************************
    Prepare soil node variables for finite difference solution
  *************************************************************/

  if(!options.QUICK_FLUX) {

    bubble_node    = soil_con->bubble_node; 
    expt_node      = soil_con->expt_node; 
    max_moist_node = soil_con->max_moist_node;  
    alpha          = soil_con->alpha; 
    beta           = soil_con->beta; 
    gamma          = soil_con->gamma; 
    moist_node     = energy->moist;
    kappa_node     = energy->kappa_node;
    Cs_node        = energy->Cs_node;
    T_node         = energy->T;
    dz_node        = soil_con->dz_node;
    ice_node       = energy->ice;

  }
  else {

    bubble_node    = NULL; 
    expt_node      = NULL; 
    max_moist_node = NULL;  
    alpha          = NULL; 
    beta           = NULL; 
    gamma          = NULL; 
    moist_node     = NULL;
    kappa_node     = NULL;
    Cs_node        = NULL;
    T_node         = NULL;
    dz_node        = NULL;
    ice_node       = NULL;

  }

  /**************************************************
    Find Surface Temperature Using Root Brent Method
  **************************************************/
  if(options.FULL_ENERGY) {

    /** Added for temporary backwards compatability **/
    if(snow->swq > 0) {
      T_lower = 0.;
      T_upper = energy->T[0]-SURF_DT;
    }
    else {
      T_lower = 0.5*(energy->T[0]+air_temp)-SURF_DT;
      T_upper = 0.5*(energy->T[0]+air_temp)+SURF_DT;
    }

#if QUICK_FS
    Tsurf = root_brent(T_upper, T_lower, ErrorString,
		       func_surf_energy_bal, T2, Ts_old, T1_old, Tair, 
		       ra, atmos_density, shortwave, longwave, albedo, 
		       emissivity, kappa1, kappa2, Cs1, Cs2, D1, D2, dp, 
		       delta_t, Le, Ls, Vapor, moist, ice0, max_moist, 
		       bubble, expt, snow_depth, snow_density, Tsnow_surf, 
		       snow_cover_fraction, surf_atten, U, displacement, 
		       roughness, ref_height, (double)soil_con->elevation, 
		       soil_con->b_infilt, soil_con->max_infil, 
		       (double)dt, atmos->vpd[hour], snow_energy, mu, 
		       rainfall, Wdew, &energy->grnd_flux, &T1, 
		       &energy->latent, &energy->sensible, 
		       &energy->deltaH, &energy->snow_flux,
		       &energy->error, soil_con->depth, 
		       soil_con->Wcr, soil_con->Wpwp, 
		       soil_con->resid_moist, T_node, Tnew_node, dz_node, 
		       kappa_node, Cs_node, moist_node, bubble_node, 
		       expt_node, max_moist_node, ice_node, alpha, beta, 
		       gamma, soil_con->ufwc_table_layer[0], 
		       soil_con->ufwc_table_node, root, layer_wet, 
		       layer_dry, veg_var_wet, veg_var_dry, VEG, 
		       veg_class, dmy->month, Nnodes, 
		       FIRST_SOLN, snow->snow, soil_con->FS_ACTIVE);
#else
    Tsurf = root_brent(T_upper, T_lower, ErrorString,
		       func_surf_energy_bal, T2, Ts_old, T1_old, Tair, 
		       ra, atmos_density, shortwave, longwave, albedo, 
		       emissivity, kappa1, kappa2, Cs1, Cs2, D1, D2, dp, 
		       delta_t, Le, Ls, Vapor, moist, ice0, max_moist, 
		       bubble, expt, snow_depth, snow_density, Tsnow_surf, 
		       snow_cover_fraction, surf_atten, U, displacement, 
		       roughness, ref_height, (double)soil_con->elevation, 
		       soil_con->b_infilt, soil_con->max_infil, 
		       (double)dt, atmos->vpd[hour], snow_energy, mu, 
		       rainfall, Wdew, &energy->grnd_flux, &T1, 
		       &energy->latent, &energy->sensible, &energy->deltaH, 
		       &energy->snow_flux, &energy->error, soil_con->depth, 
		       soil_con->Wcr, soil_con->Wpwp, 
		       soil_con->resid_moist, T_node, Tnew_node, dz_node, 
		       kappa_node, Cs_node, moist_node, bubble_node, 
		       expt_node, max_moist_node, ice_node, alpha, beta, 
		       gamma, root, layer_wet, layer_dry, veg_var_wet, 
		       veg_var_dry, VEG, veg_class, 
		       dmy->month, Nnodes, FIRST_SOLN, snow->snow, 
		       soil_con->FS_ACTIVE);
#endif

    if(Tsurf <= -9998) {  
#if QUICK_FS
      error = error_calc_surf_energy_bal(Tsurf, T2, Ts_old, T1_old, Tair, 
					 ra, atmos_density, shortwave, 
					 longwave, albedo, emissivity, kappa1, 
					 kappa2, Cs1, Cs2, D1, D2, dp, 
					 delta_t, Le, Ls, Vapor, moist, ice0, 
					 max_moist, bubble, expt, snow_depth, 
					 snow_density, Tsnow_surf, 
					 snow_cover_fraction, surf_atten, 
					 U, displacement, roughness, 
					 ref_height, 
					 (double)soil_con->elevation, 
					 soil_con->b_infilt, 
					 soil_con->max_infil, (double)dt, 
					 atmos->vpd[hour], snow_energy, mu, 
					 rainfall, Wdew, &energy->grnd_flux, 
					 &T1, &energy->latent, 
					 &energy->sensible, &energy->deltaH, 
					 &energy->snow_flux,
					 &energy->error, soil_con->depth, 
					 soil_con->Wcr, soil_con->Wpwp, 
					 soil_con->resid_moist, T_node, 
					 Tnew_node, dz_node, kappa_node, 
					 Cs_node, moist_node, bubble_node, 
					 expt_node, max_moist_node, 
					 ice_node, alpha, beta, gamma, 
					 soil_con->ufwc_table_layer[0], 
					 soil_con->ufwc_table_node, 
					 root, layer_wet, layer_dry, 
					 veg_var_wet, veg_var_dry, VEG, 
					 veg_class, 
					 dmy->month, iveg, Nnodes, FIRST_SOLN,
					 snow->snow, soil_con->FS_ACTIVE, ErrorString);
#else
      error = error_calc_surf_energy_bal(Tsurf, T2, Ts_old, T1_old, Tair, 
					 ra, atmos_density, shortwave, 
					 longwave, albedo, emissivity, kappa1, 
					 kappa2, Cs1, Cs2, D1, D2, dp, 
					 delta_t, Le, Ls, Vapor, moist, ice0, 
					 max_moist, bubble, expt, snow_depth, 
					 snow_density, Tsnow_surf, 
					 snow_cover_fraction, surf_atten, 
					 U, displacement, roughness, 
					 ref_height, 
					 (double)soil_con->elevation, 
					 soil_con->b_infilt, 
					 soil_con->max_infil, (double)dt, 
					 atmos->vpd[hour], snow_energy, mu, 
					 rainfall, Wdew, &energy->grnd_flux, 
					 &T1, &energy->latent, 
					 &energy->sensible, &energy->deltaH, 
					 &energy->snow_flux,
					 &energy->error, soil_con->depth, 
					 soil_con->Wcr, soil_con->Wpwp, 
					 soil_con->resid_moist, T_node, 
					 Tnew_node, dz_node, kappa_node, 
					 Cs_node, moist_node, bubble_node, 
					 expt_node, max_moist_node, ice_node, 
					 alpha, beta, gamma, root, layer_wet, 
					 layer_dry, veg_var_wet, veg_var_dry, 
					 VEG, veg_class, 
					 dmy->month, iveg, Nnodes, FIRST_SOLN, 
					 snow->snow, soil_con->FS_ACTIVE, ErrorString);
#endif
    }
  }
  else {

    /** Frozen soil model run with no surface energy balance **/
    Tsurf = Tair;

  }

  /**************************************************
    Recalculate Energy Balance Terms for Final Surface Temperature
  **************************************************/
#if QUICK_FS
  error = solve_surf_energy_bal(Tsurf, T2, Ts_old, T1_old, Tair, ra, 
				atmos_density, shortwave, longwave, albedo, 
				emissivity, kappa1, kappa2, Cs1, Cs2, D1, D2, 
				dp, delta_t, Le, Ls, Vapor, moist, ice0, 
				max_moist, bubble, expt, snow_depth, 
				snow_density, Tsnow_surf, 
				snow_cover_fraction, surf_atten, U, 
				displacement, roughness, ref_height, 
				(double)soil_con->elevation, 
				soil_con->b_infilt, soil_con->max_infil, 
				(double)dt, atmos->vpd[hour], snow_energy, mu, 
				rainfall, Wdew, &energy->grnd_flux, 
				&T1, &energy->latent, &energy->sensible, 
				&energy->deltaH, &energy->snow_flux, 
				&energy->error, soil_con->depth, 
				soil_con->Wcr, soil_con->Wpwp, 
				soil_con->resid_moist, T_node, 
				Tnew_node, dz_node, kappa_node, Cs_node, 
				moist_node, bubble_node, expt_node, 
				max_moist_node, ice_node, alpha, beta, gamma,  
				soil_con->ufwc_table_layer[0],  
				soil_con->ufwc_table_node, root, layer_wet, 
				layer_dry, veg_var_wet, veg_var_dry, VEG, 
				veg_class, dmy->month, Nnodes, 
				FIRST_SOLN, snow->snow, soil_con->FS_ACTIVE);
#else
  error = solve_surf_energy_bal(Tsurf, T2, Ts_old, T1_old, Tair, ra, 
				atmos_density, shortwave, longwave, albedo,