snow_intercept.c

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

/*
 * SUMMARY:      SnowInterception.c - simulates snow interception and release
 * USAGE:        
 *
 * AUTHOR:       Brian Connelly and Pascal Storck
 * ORG:          University of Washington, Department of Civil Engineering
 * E-MAIL:       [e-mail]
 * ORIG-DATE:    29-Aug-1996 at 13:42:17
 * LAST-MOD: Tue Feb 29 14:36:36 2000 by Keith Cherkauer <[e-mail]>
 * DESCRIPTION:  Calculates the interception and subsequent release of
 *               by the forest canopy using an energy balance approach
 * DESCRIP-END.
 * FUNCTIONS:    SnowInterception()
 * COMMENTS:     Modified for use with VIC-NL code by Keith Cherkauer
 *               on 4-9-98   
 */

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

static char vcid[] = "$Id: snow_intercept.c,v 4.1 2000/05/16 21:07:16 vicadmin Exp $";

/*****************************************************************************
  Function name: SnowInterception()

  Purpose      : Calculate snow interception and release by the canopy

  Required     :
    int Dt                 - Model timestep (hours)
    double F                - Fractional coverage
    double LAI              - Leaf Area Index
    double MaxInt           - Maximum rainfall interception storage (m)
    double BaseRa           - Aerodynamic resistance (uncorrected for
                             stability) (s/m)
    double AirDens          - Density of air (kg/m3)
    double EactAir          - Actual vapor pressure of air (Pa) 
    double Lv               - Latent heat of vaporization (J/kg3)
    PIXRAD *LocalRad       - Components of radiation balance for current pixel
                             (W/m2) 
    double Press            - Air pressure (Pa)
    double Tair             - Air temperature (C) 
    double Vpd	           - Vapor pressure deficit (Pa) 
    double Wind             - Wind speed (m/s)
    double *RainFall        - Amount of rain (m)
    double *Snowfall        - Amount of snow (m)
    double *IntRain         - Intercepted rain (m) 
    double *IntSnow         - Snow water equivalent of intercepted snow (m)
    double *TempIntStorage  - Temporary storage for snowmelt and rainfall
                             involved in mass release calculations (m)
    double *VaporMassFlux   - Vapor mass flux to/from intercepted snow
                             (m/timestep)
    double *Tcanopy         - Canopy temperature (C)
    double *MeltEnergy      - Energy used in heating and melting of the snow 
                             (W/m2)

  Returns      : none

  Modifies     :
    double *RainFall        - Amount of rain (m)
    double *Snowfall        - Amount of snow (m)
    double *IntRain         - Intercepted rain (m) 
    double *IntSnow         - Snow water equivalent of intercepted snow (m)
    double *TempIntStorage  - Temporary storage for snowmelt and rainfall
                             involved in mass release calculations (m)
    double *VaporMassFlux   - Vapor mass flux to/from intercepted snow
                             (m/timestep)  
    double *Tcanopy         - Canopy temperature (C)

  Comments     : Only the top canopy layer is taken into account for snow
                 interception.  Snow interception by lower canopy is
                 disregarded.  Rain water CAN be intercepted by lower canopy
                 layers (similar to InterceptionStorage()).
                 Of course:  NO vegetation -> NO interception

  Modifications:
  06-98 included maximum structural loading to prevent the model
        from loading the canopy with more snow than it can handle
        structurally.                                             PXS
  09-98 aerodynamic resistances in the canopy when snow has been  
        intercepted is increased by a factor of 10:  include REF
        Journal of Hydrology, 1998                            KAC, GO'D

*****************************************************************************/
void snow_intercept(double Dt, 
		    double F,  
		    double LAI, 
		    double MaxInt, 
		    double Ra, 
		    double AirDens,
		    double EactAir, 
		    double Lv, 
		    double Shortwave,
		    double Longwave, 
		    double Press, 
		    double Tair, 
		    double Vpd, 
		    double Wind,  
		    double *RainFall,
		    double *SnowFall, 
		    double *IntRain, 
		    double *IntSnow,
		    double *TempIntStorage, 
		    double *VaporMassFlux,
		    double *Tcanopy, 
		    double *MeltEnergy, 
		    int month, 
		    int rec,
		    int hour)
{
  double AdvectedEnergy;         /* Energy advected by the rain (W/m2) */
  double BlownSnow;              /* Depth of snow blown of the canopy (m) */
  double DeltaSnowInt;           /* Change in the physical swe of snow
				    interceped on the branches. (m) */
  double Drip;                   /* Amount of drip from intercepted snow as a
				    result of snowmelt (m) */
  double ExcessSnowMelt;         /* Snowmelt in excess of the water holding
				    capacity of the tree (m) */
  double EsSnow;                 /* saturated vapor pressure in the snow pack
                                   (Pa)  */
  double InitialSnowInt;         /* Initial intercepted snow (m) */ 
  double InitialWaterInt;        /* Initial intercepted water (snow and rain)
                                    (m) */ 
  double LatentHeat;             /* Latent heat flux (W/m2) */
  double LongOut;                /* Longwave radiation emitted by canopy 
                                   (W/m2) */
  double Ls;                     /* Latent heat of sublimation (J/(kg K) */
  double MassBalanceError;       /* Mass blalnce to make sure no water is
				    being destroyed/created (m) */
  double MaxWaterInt;            /* Water interception capacity (m) */  
  double MaxSnowInt;             /* Snow interception capacity (m) */
  double NetRadiation;
  double PotSnowMelt;            /* Potential snow melt (m) */
  double RainThroughFall;        /* Amount of rain reaching to the ground (m)
				  */ 
  double RefreezeEnergy;         /* Energy available for refreezing or melt */
  double ReleasedMass;           /* Amount of mass release of intercepted snow
				    (m) */ 
  double SensibleHeat;           /* Sensible heat flux (W/m2) */
  double SnowThroughFall;        /* Amount of snow reaching to the ground (m)
				  */ 
  double Tmp;                    /* Temporary variable */

  double Imax1;                  /* maxium water intecept regardless of temp */
  double IntRainFract;           /* Fraction of intercpeted water which is 
				    liquid */
  double IntSnowFract;           /* Fraction of intercepted water which is 
				    solid */
  double Overload;               /* temp variable to calculated structural 
				    overloading */

  /* Convert Units from VIC (mm -> m) */
  *RainFall /= 1000.;
  *SnowFall /= 1000.;
  *IntRain  /= 1000.;
  MaxInt    /= 1000.;

  /* Initialize Drip, H2O balance, and mass release variables. */
  
  InitialWaterInt = *IntSnow + *IntRain;
  
  *IntSnow /= F;
  *IntRain /= F;
  
  InitialSnowInt = *IntSnow;
  
  Drip = 0.0;
  ReleasedMass = 0.0;
  
  /* Determine the maximum snow interception water equivalent.           
     Kobayashi, D., 1986, Snow Accumulation on a Narrow Board,           
     Cold Regions Science and Technology, (13), pp. 239-245.           
     Figure 4. */  
  
  Imax1 = 4.0* LAI_SNOW_MULTIPLIER * LAI;

  if (Tair < -1.0 && Tair > -3.0)
    MaxSnowInt = (Tair*3.0/2.0) + (11.0/2.0);
  else if (Tair > -1.0) 
    MaxSnowInt = 4.0;  
  else
    MaxSnowInt = 1.0;
  
  /* therefore LAI_ratio decreases as temp decreases */
  
  MaxSnowInt *= LAI_SNOW_MULTIPLIER * LAI;
  
  /* Calculate snow interception. */  
  
  DeltaSnowInt = (1-*IntSnow/MaxSnowInt) * *SnowFall; 
  if (DeltaSnowInt + *IntSnow > MaxSnowInt) 
    DeltaSnowInt = MaxSnowInt - *IntSnow;
  if (DeltaSnowInt < 0.0)  
    DeltaSnowInt = 0.0;
  
  /* Reduce the amount of intercepted snow if windy and cold.         
     Ringyo Shikenjo Tokyo, #54, 1952.                                
     Bulletin of the Govt. Forest Exp. Station,                       
     Govt. Forest Exp. Station, Meguro, Tokyo, Japan.                 
     FORSTX 634.9072 R475r #54.                                       
     Page 146, Figure 10.                                               
     
     Reduce the amount of intercepted snow if snowing, windy, and     
     cold (< -3 to -5 C).                                             
     Schmidt and Troendle 1992 western snow conference paper. */  
  
  if (Tair < -3.0 && DeltaSnowInt > 0.0 && Wind > 1.0) {
    BlownSnow = (0.2 * Wind - 0.2) * DeltaSnowInt;
    if (BlownSnow >= DeltaSnowInt) 
      BlownSnow = DeltaSnowInt;
    DeltaSnowInt -= BlownSnow;
  }
  
  /* now update snowfall and total accumulated intercepted snow amounts */

  if (*IntSnow +  DeltaSnowInt > Imax1) DeltaSnowInt =0.0; 
  
  /* pixel depth    */ 
  SnowThroughFall = (*SnowFall - DeltaSnowInt) * F + (*SnowFall) * (1 - F);
  
  /* physical depth */
  *IntSnow += DeltaSnowInt;
  
  /* Calculate amount of rain intercepted on branches and stored in
     intercepted snow. */