snow_intercept.c

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

  /* physical depth */
  MaxWaterInt = LIQUID_WATER_CAPACITY * (*IntSnow) + MaxInt;
  
  if ((*IntRain + *RainFall) <= MaxWaterInt) {
    /* physical depth */
    *IntRain += *RainFall;
    /* pixel depth */
    RainThroughFall = *RainFall * (1 - F);      
  }
  else {
    /* pixel depth */
    RainThroughFall = (*IntRain + *RainFall - MaxWaterInt) * F + 
      (*RainFall * (1 - F));
    /* physical depth */
    *IntRain = MaxWaterInt;
  }

  /* at this point we have calculated the amount of snowfall intercepted and
     the amount of rainfall intercepted.  These values have been 
     appropriately subtracted from SnowFall and RainFall to determine 
     SnowThroughfall and RainThroughfall.  However, we can end up with the 
     condition that the total intercepted rain plus intercepted snow is 
     greater than the maximum bearing capacity of the tree regardless of air 
     temp (Imax1).  The following routine will adjust *IntRain and *IntSnow 
     by triggering mass release due to overloading.  Of course since *IntRain
     and *IntSnow are mixed, we need to slough them of as fixed fractions  */

  if (*IntRain + *IntSnow > Imax1) { /*then trigger structural unloading*/
    Overload = (*IntSnow + *IntRain) - Imax1;
    IntRainFract= *IntRain/(*IntRain + *IntSnow);
    IntSnowFract = *IntSnow/(*IntRain + *IntSnow);
    *IntRain = *IntRain - Overload*IntRainFract;
    *IntSnow = *IntSnow - Overload*IntSnowFract;
    RainThroughFall = RainThroughFall + (Overload*IntRainFract)*F;
    SnowThroughFall = SnowThroughFall + (Overload*IntSnowFract)*F;
  }
  
  /* The canopy temperature is assumed to be equal to the air temperature if 
     the air temperature is below 0C, otherwise the canopy temperature is 
     equal to 0C */
  
  if (Tair > 0.)
    *Tcanopy = 0.;
  else
    *Tcanopy = Tair;

  /* Calculate the net radiation at the canopy surface, using the canopy 
     temperature.  The outgoing longwave is subtracted twice, because the 
     canopy radiates in two directions */

  Tmp = *Tcanopy + 273.15;
  LongOut = STEFAN_B * (Tmp * Tmp * Tmp * Tmp);
  NetRadiation = (1.-NEW_SNOW_ALB)*Shortwave + Longwave - 2 * F * LongOut;
  NetRadiation /= F;

  /* Calculate the vapor mass flux between the canopy and the surrounding 
     air mass */
  
  EsSnow = svp(*Tcanopy); 

  /** Added division by 10 to incorporate change in canopy resistance due
      to smoothing by intercepted snow **/
  *VaporMassFlux = AirDens * (0.622/Press) * (EactAir - EsSnow) / Ra / 10.; 
  *VaporMassFlux /= RHO_W; 

  if (Vpd == 0.0 && *VaporMassFlux < 0.0)
    *VaporMassFlux = 0.0;
  
  /* Calculate the latent heat flux */

  Ls = (677. - 0.07 * *Tcanopy) * 4.1868 * 1000.;
  LatentHeat = Ls * *VaporMassFlux * RHO_W;

  /* Calculate the sensible heat flux */

  SensibleHeat = AirDens * Cp * (Tair - *Tcanopy)/Ra;
  
  /* Calculate the advected energy */

  AdvectedEnergy = (4186.8 * Tair * *RainFall)/(Dt * SECPHOUR);

  /* Calculate the amount of energy available for refreezing */

  RefreezeEnergy = SensibleHeat + LatentHeat + NetRadiation + AdvectedEnergy;

  RefreezeEnergy *= Dt * SECPHOUR;

  /* if RefreezeEnergy is positive it means energy is available to melt the
     intercepted snow in the canopy.  If it is negative, it means that 
     intercepted water will be refrozen */
  
  /* Update maximum water interception storage */
  
  MaxWaterInt = LIQUID_WATER_CAPACITY * (*IntSnow) + MaxInt;

  /* Convert the vapor mass flux from a flux to a depth per interval */
  *VaporMassFlux *= Dt * SECPHOUR;
  
  if (RefreezeEnergy > 0.0) {

    if (-(*VaporMassFlux) > *IntRain) {
      *VaporMassFlux = -(*IntRain);
      *IntRain = 0.;
    }
    else
      *IntRain += *VaporMassFlux;

    PotSnowMelt = min((RefreezeEnergy/Lf/RHO_W), *IntSnow);

    *MeltEnergy -= (Lf * PotSnowMelt * RHO_W) / (Dt *SECPHOUR);
    
    if ((*IntRain + PotSnowMelt) <= MaxWaterInt) {

      *IntSnow -= PotSnowMelt;
      *IntRain += PotSnowMelt;
      PotSnowMelt = 0.0;

    }
    
    else {

      ExcessSnowMelt = PotSnowMelt + *IntRain - MaxWaterInt;
      
      *IntSnow -= MaxWaterInt - (*IntRain);
      *IntRain = MaxWaterInt;
      if (*IntSnow < 0.0) 
        *IntSnow = 0.0;
      
      if (SnowThroughFall > 0.0 && 
	  InitialSnowInt <= MIN_INTERCEPTION_STORAGE) {
        /* Water in excess of MaxWaterInt has been generated.  If it is 
           snowing and there was little intercepted snow at the beginning 
	   of the time step ( <= MIN_INTERCEPTION_STORAGE), then allow the 
	   snow to melt as it is intercepted */
        Drip += ExcessSnowMelt; 
        *IntSnow -= ExcessSnowMelt;
        if (*IntSnow < 0.0) 
          *IntSnow = 0.0;
      }
      else 
      /* Else, SnowThroughFall = 0.0 or SnowThroughFall > 0.0 and there is a 
         substantial amount of intercepted snow at the beginning of the time 
         step ( > MIN_INTERCEPTION_STORAGE).  Snow melt may generate mass 
         release. */
        *TempIntStorage += ExcessSnowMelt;
      
      MassRelease(IntSnow, TempIntStorage, &ReleasedMass, &Drip);
    }
    
    /* If intercepted snow has melted, add the water it held to drip */
    
    MaxWaterInt = LIQUID_WATER_CAPACITY * (*IntSnow) + MaxInt;
    if (*IntRain > MaxWaterInt) {
      Drip += *IntRain - MaxWaterInt;
      *IntRain = MaxWaterInt;
    }
  }  
  
  else /* else (RefreezeEnergy <= 0.0) */ {
    
    /* Reset *TempIntStorage to 0.0 when energy balance is negative */
    
    *TempIntStorage = 0.0;
    
    /* Refreeze as much surface water as you can */
    
    if (RefreezeEnergy > - (*IntRain) * Lf) {
      *IntSnow += fabs(RefreezeEnergy) / Lf;
      *IntRain -= fabs(RefreezeEnergy) / Lf;

      *MeltEnergy += (fabs(RefreezeEnergy) * RHO_W) / (Dt *SECPHOUR);

      RefreezeEnergy = 0.0;
    }

    else {
      
      /* All of the water in the surface layer has been frozen. */
      
      *IntSnow += *IntRain;
     
      /* Added on April 8 as a test */
      /*       RefreezeEnergy += *IntRain*Lf; */
      /*       *VaporMassFlux = MAX(*VaporMassFlux,  */
      /*                            RefreezeEnergy/(Ls * RHO_W)); */
      
      /* Energy released by freezing of intercepted water is added to the 
         MeltEnergy */

      *MeltEnergy += (Lf * *IntRain * RHO_W) / (Dt *SECPHOUR);
      *IntRain = 0.0;
      
    } 
    
    if (-(*VaporMassFlux) > *IntSnow) {
      *VaporMassFlux = -(*IntSnow);
      *IntSnow = 0.0;
    }
    else
      *IntSnow += *VaporMassFlux;
  } 
  
  *IntSnow *= F;
  *IntRain *= F;
  *MeltEnergy *= F;
  *VaporMassFlux *= F;
  Drip           *= F;
  ReleasedMass   *= F;
  
  /* Calculate intercepted H2O balance. */  
  
  MassBalanceError = (InitialWaterInt - (*IntSnow + *IntRain))
                   + (*SnowFall + *RainFall) - (SnowThroughFall
                   + RainThroughFall + Drip + ReleasedMass)
                   + *VaporMassFlux;

  *RainFall = RainThroughFall + Drip;
  *SnowFall = SnowThroughFall + ReleasedMass;

  /* Convert Units to VIC (m -> mm) */
  *VaporMassFlux *= -1.;
  *RainFall *= 1000.;
  *SnowFall *= 1000.;
  *IntRain  *= 1000.;

}