solve_snow.c

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

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

static char vcid[] = "$Id: solve_snow.c,v 4.2.2.3 2004/05/06 19:57:35 tbohn Exp $";

double solve_snow(snow_data_struct    *snow,
		  layer_data_struct   *layer_wet,
		  layer_data_struct   *layer_dry,
		  veg_var_struct      *veg_var_wet,
		  veg_var_struct      *veg_var_dry,
		  int                  month,
		  int                  day_in_year,
		  energy_bal_struct   *energy,
		  soil_con_struct     *soil_con,
		  char                 overstory,
		  int                  dt,
		  int                  rec,
		  int                  veg_class,
		  int                  iveg,
		  int                  Nveg,
		  int                  band,
		  int                  hour,
		  int                  Nnodes,
		  double               shortwave,
		  double               longwave,
		  double               air_temp,
		  double               prec,
		  double               density,
		  double               vp,
		  double               vpd,
		  double               pressure,
		  double               mu,
		  double               roughness,
		  double               displacement,
		  double               ref_height,
		  double               surf_atten,
		  double               MAX_SNOW_TEMP,
		  double               MIN_RAIN_TEMP,
		  double               wind_h,
		  double               moist,
		  double               ice0,
		  double               dp,
		  double               bare_albedo,
		  double              *rainfall,
		  double              *out_prec,
		  double              *Le,
		  double              *Ls,
		  double              *aero_resist,
		  double              *tmp_wind,
		  double              *net_short,
		  double              *out_short,
		  double              *rad,
		  double              *Evap,
		  double              *tmp_snow_energy,
		  double              *snow_inflow,
		  double              *ppt,
		  double              *gauge_correction,
		  float               *root) {
/*********************************************************************
  solve_snow.c                Keith Cherkauer       July 2, 1998

  This routine was written to handle the various calls and data 
  handling needed to solve the various components of the new VIC
  snow code for both the full_energy and water_balance models.

  Returns snow, veg_var, and energy variables for each elevation
  band.  Variable ppt[] is defined for elevation bands with snow.

  07-13-98 modified to use elevation bands when solving the 
           snow model                                        KAC
  11-30-98 reworked the way the snow/rain fraction is computed
           and added to check to assure that very small amounts 
           of snow do not fall, causing snow sublimation to 
           be calculated.  (found by Greg)                   KAC
  02-29-00 removed ground heat flux computations, will now 
           make those outside of this routine, in the same function
	   that is used to compute the ground heat flux when
	   there is no snow cover.                           KAC
  06-04-03 Counter for number of days since last snowfall was 
           incremented twice in the MELTING update.  This has been
           fixed.                                            KAC
  06-04-03 Added check so that MELTING flag is only TRUE if melt
           occurs in the melt season - currently this is defined
           between March 1 and October 1.  Otherwise the MELTING
           flag can trigger rapid very early season melt

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

  extern option_struct   options;
  extern veg_lib_struct *veg_lib;

  char                ErrStr[MAXSTRING];
  char                FIRST_SOLN[1];
  double              rainonly;
  double              canopy_temp;
  double              tmp_energy_val = 0.;
  double              old_swq;
  double              tmp_Wdew[2];
  double              melt;
  double              snowfall[2];
  double              snow_coverage;
  double              grnd_temp;
  double              tmp_rain;
  double              surf_long;
  double              store_snowfall;
  int                 curr_snow;

  melt     = 0.;
  ppt[WET] = 0.; 
  ppt[DRY] = 0.; 

  /** Calculate Fraction of Precipitation that falls as Rain **/
  rainonly      = calc_rainonly(air_temp, prec, 
				MAX_SNOW_TEMP, MIN_RAIN_TEMP, mu);
  snowfall[WET] = gauge_correction[SNOW] * (prec - rainonly);
  rainfall[WET] = gauge_correction[RAIN] * rainonly;
  snowfall[DRY] = 0.;
  rainfall[DRY] = 0.;
  if(snowfall[WET] < 1e-5) snowfall[WET] = 0.;
  (*out_prec) = snowfall[WET] + rainfall[WET];
  store_snowfall = snowfall[WET];

  /** Compute latent heats **/
  (*Le) = (2.501 - 0.002361 * air_temp) * 1.0e6;
  (*Ls) = 0.;

  /** Error checks **/
  if((snow->swq > 0 || snowfall[WET] > 0.
      || (snow->snow_canopy>0. && overstory))) {
    if(mu!=1 && options.FULL_ENERGY) {
      sprintf(ErrStr,"Snow model cannot be used if mu (%f) is not equal to 1.\n\tsolve_snow.c: record = %i,\t vegetation type = %i",
	      mu, rec, iveg);
      vicerror(ErrStr);
    }
    else if(mu!=1) {
      fprintf(stderr,"WARNING: Snow is falling, but mu not equal to 1 (%f)\n", mu);
      fprintf(stderr,"\trec = %i, veg = %i, hour = %i\n",rec,iveg,hour);
    }
  }

  if((snow->swq > 0 || snowfall[WET] > 0.
      || (snow->snow_canopy>0. && overstory)) && mu==1) {
    
    /************************************************
      Snow is Present or Falling 
    ************************************************/

    /** Snow is present or falling **/
    snow->snow = TRUE;
	
    /** Initialize variables **/
    (*tmp_snow_energy) = 0.;
    if(!overstory) surf_atten = 1.;
    snow_coverage = snow->coverage;
      
    /** Compute Snow Pack Albedo **/
    if(snow->swq > 0 || snowfall[WET] > 0.) {
      snow->albedo   = snow_albedo( snowfall[WET], snow->swq, 
				    snow->coldcontent, dt, snow->last_snow,
				    snow->MELTING );
      energy->albedo = snow->albedo;
    }
    else {
      snow->albedo   = NEW_SNOW_ALB;
      energy->albedo = bare_albedo;
    }
    
    /** Age Snow Pack **/
    //if( snowfall[WET] > 0 )
      //snow->last_snow = 1;
    //else snow->last_snow++;
    
    /** Compute Radiation Balance over Snow **/ 
    (*out_short) = energy->albedo * shortwave;
    (*net_short) = (1.0 - energy->albedo) * shortwave;
    if(snow->swq > 0) {
      (*rad) = (*net_short) + longwave - STEFAN_B 
	* (snow->surf_temp+KELVIN) * (snow->surf_temp+KELVIN) 
	* (snow->surf_temp+KELVIN) * (snow->surf_temp+KELVIN);
    }
    else {
      if(options.FULL_ENERGY || options.FROZEN_SOIL) {
	(*rad) = (*net_short) + longwave - STEFAN_B 
	  * (energy->T[0]+KELVIN) * (energy->T[0]+KELVIN) 
	  * (energy->T[0]+KELVIN) * (energy->T[0]+KELVIN);
      }
      else {
	(*rad) = (*net_short) + longwave - STEFAN_B 
	  * (air_temp+KELVIN) * (air_temp+KELVIN) 
	  * (air_temp+KELVIN) * (air_temp+KELVIN);
      }
    }

    if(iveg!=Nveg) {
      
      /****************************************
	Check Vegetation for Intercepted Snow
      ****************************************/
      
      if(overstory) {
	if( snowfall[WET] > 0. || snow->snow_canopy > 0. ) {
	  
	  /** Compute Canopy Interception, if Snow and Canopy **/
	  if(snow->swq > 0) 
	    surf_long = STEFAN_B * (snow->surf_temp + KELVIN) 
	      * (snow->surf_temp + KELVIN) * (snow->surf_temp + KELVIN) 
	      * (snow->surf_temp + KELVIN);
	  else {
	    if(options.FULL_ENERGY || options.FROZEN_SOIL)
	      surf_long = STEFAN_B * (energy->T[0] + KELVIN) 
		* (energy->T[0] + KELVIN) * (energy->T[0] + KELVIN) 
		* (energy->T[0] + KELVIN);
	    else
	      surf_long = STEFAN_B * (air_temp + KELVIN) 
		* (air_temp + KELVIN) * (air_temp + KELVIN) 
		* (air_temp + KELVIN);
	  }
	  snow_intercept((double)dt, 1., veg_lib[veg_class].LAI[month-1], 
			 veg_lib[veg_class].Wdmax[month-1], 
			 aero_resist[1], density, vp, (*Le), shortwave, 
			 longwave+surf_long, pressure, air_temp, vpd, 
			 tmp_wind[1], rainfall, snowfall, &veg_var_wet->Wdew, 
			 &snow->snow_canopy, &snow->tmp_int_storage, 
			 &snow->canopy_vapor_flux, &canopy_temp, 
			 &tmp_energy_val, month, rec, hour);

	  /* Store throughfall from canopy */
	  veg_var_wet->throughfall = rainfall[0] + snowfall[0];

	  /* Estimate longwave radiation from canopy */
	  energy->longwave = STEFAN_B * (canopy_temp+KELVIN) 
	    * (canopy_temp+KELVIN) * (canopy_temp+KELVIN)