solve_snow.c

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

	    * (canopy_temp+KELVIN);

	}
	else {
	  /** Compute Canopy Evaporation, if Canopy and No Snow **/
	  tmp_Wdew[WET] = veg_var_wet->Wdew;
	  if(options.DIST_PRCP) tmp_Wdew[DRY] = 0.;
	  Evap[0] = canopy_evap(layer_wet, layer_dry, veg_var_wet, 
				veg_var_dry, FALSE, veg_class, month, mu, 
				tmp_Wdew, (double)dt, rad[0], vpd, 
				net_short[0], air_temp, aero_resist[0], 
				displacement, roughness, ref_height, 
				(double)soil_con->elevation, rainfall, 
				soil_con->depth, soil_con->Wcr, 
				soil_con->Wpwp, root);

	  /* Store throughfall from canopy */
	  rainfall[WET] = veg_var_wet->throughfall;
	  if(options.DIST_PRCP) 
	    rainfall[DRY] = veg_var_dry->throughfall;

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

	}

      }

      else if(snowfall[0] > 0. && veg_var_wet->Wdew > 0.) {

	/** If No Overstory, Empty Vegetation of Stored Water **/

	tmp_rain = calc_rainonly(air_temp, veg_var_wet->Wdew,
				 MAX_SNOW_TEMP, MIN_RAIN_TEMP, mu);
	rainfall[WET]            += tmp_rain;
	snowfall[WET]            += veg_var_wet->Wdew - tmp_rain;
	veg_var_wet->throughfall  = rainfall[WET] + snowfall[WET];
	veg_var_wet->Wdew         = 0.;
	energy->longwave          = longwave;

      }
      else {

	/** Precipitation "Passes Through" Vegetation which 
	    is Under Snow (used only for accounting purposes)**/

	veg_var_wet->throughfall = rainfall[WET] + snowfall[WET];
	veg_var_dry->throughfall = rainfall[DRY] + snowfall[DRY];
	energy->longwave         = longwave;

      }

    }
    else energy->longwave = longwave;
    
    old_swq = snow->swq;
    
    if(snow->swq>0.0 || snowfall[0] > 0) {
      
      /******************************
	Snow Pack Present on Ground
      ******************************/

      // store snowfall reaching the ground for determining the albedo
      store_snowfall            = snowfall[WET];

      /** Age Snowpack **/
      if( snowfall[WET] > 0 ) curr_snow = 1; // new snow - reset pack age
      else curr_snow = snow->last_snow + 1; // age pack by one time step

      if(overstory && snowfall[0] > 0.) {

	/** recompute surface properties if overstory drops snow **/

	snow->albedo    = NEW_SNOW_ALB;
	energy->albedo  = snow->albedo;
	snow->last_snow = 1;
	(*out_short)    = energy->albedo * shortwave;
	(*net_short)    = (1.0 - energy->albedo) * shortwave;
	(*rad)          = (*net_short) + longwave - STEFAN_B 
	  * (snow->surf_temp+KELVIN) * (snow->surf_temp+KELVIN) 
	  * (snow->surf_temp+KELVIN) * (snow->surf_temp+KELVIN);
      }
    
      (*snow_inflow) += rainfall[WET] + snowfall[WET];
      if(options.FULL_ENERGY || options.FROZEN_SOIL) grnd_temp = energy->T[0];
      else grnd_temp = air_temp;

      /** Call snow pack accumulation and ablation algorithm **/

      snow_melt(soil_con, rec, iveg, wind_h+soil_con->snow_rough, 
		aero_resist[2], Le[0], snow, (double)dt, 0.00, 
		soil_con->snow_rough, surf_atten, rainfall[WET], 
		snowfall[WET], tmp_wind[2], grnd_temp, air_temp, net_short[0], 
		energy->longwave, density, pressure, vpd, vp, &melt, 
		&energy->advection, &energy->deltaCC, &energy->grnd_flux, 
		&energy->latent, &energy->sensible, &energy->error, 
		&energy->refreeze_energy);

      ppt[WET] += melt;
      energy->albedo   = snow->albedo;
      tmp_snow_energy[0] = energy->advection - energy->deltaCC 
	+ energy->refreeze_energy;
      (*Ls) = (677. - 0.07 * snow->surf_temp) * 4.1868 * 1000;
      
      /** Compute Snow Parameters **/
      if(snow->swq > 0.) {

	/** Calculate Snow Density **/
	if ( snow->surf_temp <= 0 )
	  // snowpack present, compress and age density
	  snow->density = snow_density(day_in_year, snowfall[WET], air_temp, 
				       old_swq, snow->depth, snow->coldcontent, 
				       (double)dt, snow->surf_temp);
	else 
	  // no snowpack present, start with new snow density
	  if ( curr_snow == 1 ) 
	    snow->density = new_snow_density(air_temp);

	/** Calculate Snow Depth (H.B.H. 7.2.1) **/
	snow->depth = 1000. * snow->swq / snow->density; 
	
	/** Record if snowpack is melting this time step **/
	if ( snow->coldcontent >= 0 && day_in_year > 60 // ~ March 1
	     && day_in_year < 273 // ~ October 1
	     ) snow->MELTING = TRUE;
	else if ( snow->MELTING && snowfall[WET] > TraceSnow ) 
	  snow->MELTING = FALSE;
	
	/** Check for Thin Snowpack which only Partially Covers Grid Cell **/
	if(snow->swq < MAX_FULL_COVERAGE_SWQ) {
	  snow->coverage = 1.; 
/* 	  snow->coverage = 1. / MAX_FULL_COVERAGE_SWQ   */
/* 	    * snow->swq;  */
/* 	  if(snow_coverage <= snow->coverage) { */
/* 	    energy->albedo = bare_albedo; */
/* 	  } */
/* 	  else { */
/* 	    energy->albedo = (snow_coverage - snow->coverage)  */
/* 	      / (1. - snow->coverage) * snow->albedo; */
/* 	    energy->albedo += (1. - snow_coverage)  */
/* 	      / (1. - snow->coverage) * bare_albedo; */
/* 	  } */
	}
	else {
	  snow->coverage = 1.;
	}

	/** Estimate net longwave at the snow surface */
	energy->longwave -= (STEFAN_B * (snow->surf_temp+KELVIN) 
			     * (snow->surf_temp+KELVIN) 
			     * (snow->surf_temp+KELVIN) 
			     * (snow->surf_temp+KELVIN));

      }
      else {

	/** Reset Snow Pack Variables after Complete Melt **/

	snow->density    = 0.;
	snow->depth      = 0.;
	snow->surf_water = 0;
	snow->pack_water = 0;
	snow->surf_temp  = 0;
	snow->pack_temp  = 0;
	snow->coverage   = 0;
	snow->MELTING    = FALSE;

 	energy->albedo = (snow_coverage - snow->coverage)  
 	  / (1. - snow->coverage) * snow->albedo; 
 	energy->albedo += (1. - snow_coverage)  
 	  / (1. - snow->coverage) * bare_albedo; 
	
      }
    }

    else {

      /** Ground Snow not Present, and Falling Snow Does not Reach Ground **/

      ppt[WET] += rainfall[WET];
      tmp_snow_energy[0] = 0.;
      curr_snow               = 0;
      snow->MELTING           = FALSE;

    }
    
    if ( store_snowfall > TraceSnow || store_snowfall == 0 ) {
      // reset snow albedo ago if new snow is sufficiently deep
      snow->last_snow = curr_snow;
    }
    else {
      snow->last_snow++;
    }

  }
  else {
    
    /*****************************
      No Snow Present or Falling
    *****************************/

    /** Initialize variables **/
    snow->snow         = FALSE;
    tmp_snow_energy[0] = 0;
    energy->albedo     = bare_albedo;
    energy->longwave   = longwave;

    /** Compute Radiation Balance for Bare Surface **/ 
    out_short[0] = energy->albedo * shortwave;
    net_short[0] = (1.0 - energy->albedo) * shortwave;
    if(options.FULL_ENERGY || options.FROZEN_SOIL) {
      rad[0]       = net_short[0] + longwave 
	- STEFAN_B * (energy->T[0]+KELVIN) * (energy->T[0]+KELVIN) 
	* (energy->T[0]+KELVIN) * (energy->T[0]+KELVIN);
    }
    else {
      rad[0]       = net_short[0] + longwave 
	- STEFAN_B * (air_temp+KELVIN) * (air_temp+KELVIN) 
	* (air_temp+KELVIN) * (air_temp+KELVIN);
    }
    snow->MELTING        = FALSE;
  }

  energy->shortwave = (*net_short);

  return(melt);

}