snow_melt.c

降雪及积雪计算的源代码

				     snow->depth, snow->density,surf_atten,
				     &grnd_flux,&latent_heat,
				     &sensible_heat);

    snow->vapor_flux = vapor_flux;
    save_refreeze_energy[0] = RefreezeEnergy;

    /* since we iterated, the surface layer is below freezing and no snowmelt
     */ 
    
    SnowMelt = 0.0;
    
    /* Since updated snow_temp < 0.0, all of the liquid water in the surface
       layer has been frozen */ 
    
    SurfaceSwq += snow->surf_water;
    Ice        += snow->surf_water;
    snow->surf_water  = 0.0;
    melt_energy += snow->surf_water * Lf 
                       * RHO_W/(delta_t * SECPHOUR);
    
    /* Convert mass flux to a depth per timestep and adjust SurfaceSwq */
    
    snow->vapor_flux *= delta_t * SECPHOUR;

    if (SurfaceSwq < -(snow->vapor_flux)) {
      snow->vapor_flux = -SurfaceSwq;
      SurfaceSwq    = 0.0;
      Ice           = PackSwq;
    }
    else {
      SurfaceSwq += snow->vapor_flux;
      Ice += snow->vapor_flux;
    }

  }
  
  /* Done with iteration etc, now Update the liquid water content of the
     surface layer */ 
  
  MaxLiquidWater = LIQUID_WATER_CAPACITY * SurfaceSwq;
  if  (snow->surf_water > MaxLiquidWater) {
    melt[0]    = snow->surf_water - MaxLiquidWater;
    snow->surf_water = MaxLiquidWater;
  }
  else
    melt[0] = 0.0;
  
  /* Refreeze liquid water in the pack.                                   
     variable 'RefreezeEnergy' is the heat released to the snow pack
     if all liquid water were refrozen.                                   
     if RefreezeEnergy < PackCC then all water IS refrozen           
     PackCC always <=0.0 

     WORK IN PROGRESS: This energy is NOT added to MeltEnergy, since this does 
     not involve energy transported to the pixel.  Instead heat from the snow 
     pack is used to refreeze water */
  
  snow->pack_water += melt[0]; /* add surface layer outflow to pack 
                                      liquid water*/
  RefreezeEnergy = snow->pack_water * Lf * RHO_W;

  /* calculate energy released to freeze*/
  
  if (PackCC < -RefreezeEnergy) { /* cold content not fully depleted*/
    PackSwq   += snow->pack_water;    /* refreeze all water and update*/
    Ice       += snow->pack_water;
    snow->pack_water = 0.0;
    if (PackSwq > 0.0) {
      PackCC = PackSwq * CH_ICE * snow->pack_temp + RefreezeEnergy;
      snow->pack_temp = PackCC / (CH_ICE * PackSwq);
      if(snow->pack_temp > 0.) snow->pack_temp = 0.;
    }
    else 
      snow->pack_temp = 0.0;
  }
  else { 
    /* cold content has been either exactly satisfied or exceeded. If
       PackCC = refreeze then pack is ripe and all pack water is
       refrozen, else if energy released in refreezing exceeds PackCC 
       then exactly the right amount of water is refrozen to satify PackCC.
       The refrozen water is added to PackSwq and Ice */

    snow->pack_temp      = 0.0;    
    DeltaPackSwq = -PackCC/(Lf * RHO_W); 
    snow->pack_water -= DeltaPackSwq;
    PackSwq += DeltaPackSwq;
    Ice += DeltaPackSwq;
  }
  
  /* Update the liquid water content of the pack */
  
  MaxLiquidWater = LIQUID_WATER_CAPACITY * PackSwq;
  if (snow->pack_water > MaxLiquidWater) {
    melt[0]    = snow->pack_water - MaxLiquidWater;
    snow->pack_water = MaxLiquidWater;
  }
  else
    melt[0] = 0.0;

  /* Update snow properties */
  
  Ice  = PackSwq + SurfaceSwq;

  if (Ice > MAX_SURFACE_SWE) {
    SurfaceCC   = CH_ICE * snow->surf_temp * SurfaceSwq;
    PackCC      = CH_ICE * snow->pack_temp * PackSwq;
    if (SurfaceSwq > MAX_SURFACE_SWE) {
      PackCC     += SurfaceCC * (SurfaceSwq - MAX_SURFACE_SWE) / SurfaceSwq; 
      SurfaceCC  -= SurfaceCC * (SurfaceSwq - MAX_SURFACE_SWE) / SurfaceSwq; 
      PackSwq    += SurfaceSwq - MAX_SURFACE_SWE;
      SurfaceSwq -= SurfaceSwq - MAX_SURFACE_SWE;
    }
    else if ( SurfaceSwq < MAX_SURFACE_SWE) {
      PackCC     -= PackCC * (MAX_SURFACE_SWE - SurfaceSwq) / PackSwq; 
      SurfaceCC  += PackCC * (MAX_SURFACE_SWE - SurfaceSwq) / PackSwq; 
      PackSwq    -= MAX_SURFACE_SWE - SurfaceSwq;
      SurfaceSwq += MAX_SURFACE_SWE - SurfaceSwq;
    }
    snow->pack_temp      = PackCC / (CH_ICE * PackSwq);
    snow->surf_temp      = SurfaceCC / (CH_ICE * SurfaceSwq);
  }
  else {
    PackSwq = 0.0;
    PackCC  = 0.0;
    snow->pack_temp  = 0.0;
  }

  snow->swq = Ice + snow->pack_water + snow->surf_water;

  if (snow->swq == 0.0) {
    snow->surf_temp = 0.0;
    snow->pack_temp = 0.0;
  }
    
  /* Mass balance test */
  
  MassBalanceError = (InitialSwq - snow->swq) + (RainFall + SnowFall) 
    - melt[0] + snow->vapor_flux; 
  
/*  printf("%d %d %g\n", y, x, MassBalanceError);*/

  melt[0] *= 1000.; /* converts back to mm */
  snow->mass_error = MassBalanceError;
  snow->coldcontent = SurfaceCC;
  snow->vapor_flux *= -1.;
  *save_advection = advection;
  *save_deltaCC = deltaCC;
  *save_grnd_flux = grnd_flux;
  *save_latent = latent_heat;
  *save_sensible = sensible_heat;
  *save_Qnet = Qnet;

}

/*****************************************************************************
  Function name: CalcSnowPackEnergyBalance()

  Purpose      : Dummy function to make a direct call to
                 SnowEnergyBalance() possible.

  Required     : 
    double TSurf - SnowPack surface temperature (C)
    other arguments required by SnowPackEnergyBalance()

  Returns      :
    double Qnet - Net energy exchange at the SnowPack snow surface (W/m^2)

  Modifies     : none

  Comments     : function is local to this module
*****************************************************************************/
double CalcSnowPackEnergyBalance(double Tsurf, ...)
{

  va_list ap;                   /* Used in traversing variable argument list
                                 */ 
  double Qnet;                   /* Net energy exchange at the SnowPack snow
                                   surface (W/m^2) */

  va_start(ap, Tsurf);
  Qnet = SnowPackEnergyBalance(Tsurf, ap);
  va_end(ap);
  
  return Qnet;
}

double ErrorSnowPackEnergyBalance(double Tsurf, ...)
{

  va_list ap;                   /* Used in traversing variable argument list
                                 */ 
  double Qnet;                   /* Net energy exchange at the SnowPack snow
                                   surface (W/m^2) */

  va_start(ap, Tsurf);
  Qnet = ErrorPrintSnowPackEnergyBalance(Tsurf, ap);
  va_end(ap);
  
  return Qnet;
}

double ErrorPrintSnowPackEnergyBalance(double TSurf, va_list ap)
{


  double Dt;                     /* Model time step (hours) */
  double Ra;                     /* Aerodynamic resistance (s/m) */
  double Z;                      /* Reference height (m) */
  double Displacement;           /* Displacement height (m) */
  double Z0;                     /* surface roughness height (m) */
  double Wind;                   /* Wind speed (m/s) */
  double ShortRad;               /* Net incident shortwave radiation (W/m2) */
  double LongRadIn;              /* Incoming longwave radiation (W/m2) */
  double AirDens;                /* Density of air (kg/m3) */
  double Lv;                     /* Latent heat of vaporization (J/kg3) */
  double Tair;                   /* Air temperature (C) */
  double Press;                  /* Air pressure (Pa) */
  double Vpd;			/* Vapor pressure deficit (Pa) */
  double EactAir;                /* Actual vapor pressure of air (Pa) */
  double Rain;                   /* Rain fall (m/timestep) */
  double SweSurfaceLayer;        /* Snow water equivalent in surface layer (m)
                                 */ 
  double SurfaceLiquidWater;     /* Liquid water in the surface layer (m) */
  double OldTSurf;               /* Surface temperature during previous time
                                   step */ 
  double *GroundFlux;		/* Ground Heat Flux (W/m2) */
  double *RefreezeEnergy;        /* Refreeze energy (W/m2) */
  double *VaporMassFlux;          /* Mass flux of water vapor to or from the
                                   intercepted snow */
  double *AdvectedEnergy;         /* Energy advected by precipitation (W/m2) */
  double *DeltaColdContent;       /* Change in cold content (W/m2) */
  double TGrnd;     
  double SnowDepth;  
  double SnowDensity; 
  double SurfAttenuation; 
  char *ErrorString;
  
  /* end of list of arguments in variable argument list */

  double *LatentHeat;		/* Latent heat exchange at surface (W/m2) */
  double *SensibleHeat;		/* Sensible heat exchange at surface (W/m2) */

  /* initialize variables */
  Dt                 = (double) va_arg(ap, double);
  Ra                 = (double) va_arg(ap, double);
  Z                  = (double) va_arg(ap, double);
  Displacement       = (double) va_arg(ap, double);
  Z0                 = (double) va_arg(ap, double);
  Wind               = (double) va_arg(ap, double);
  ShortRad           = (double) va_arg(ap, double);
  LongRadIn          = (double) va_arg(ap, double);
  AirDens            = (double) va_arg(ap, double);
  Lv                 = (double) va_arg(ap, double);
  Tair               = (double) va_arg(ap, double);
  Press              = (double) va_arg(ap, double);
  Vpd                = (double) va_arg(ap, double);
  EactAir            = (double) va_arg(ap, double);
  Rain               = (double) va_arg(ap, double);
  SweSurfaceLayer    = (double) va_arg(ap, double);
  SurfaceLiquidWater = (double) va_arg(ap, double);
  OldTSurf           = (double) va_arg(ap, double);
  RefreezeEnergy     = (double *) va_arg(ap, double *);
  VaporMassFlux      = (double *) va_arg(ap, double *);
  AdvectedEnergy     = (double *) va_arg(ap, double *);
  DeltaColdContent   = (double *) va_arg(ap, double *);
  TGrnd              = (double) va_arg(ap, double);
  SnowDepth          = (double) va_arg(ap, double);
  SnowDensity        = (double) va_arg(ap, double);
  SurfAttenuation    = (double) va_arg(ap, double);
  GroundFlux         = (double *) va_arg(ap, double *);
  LatentHeat         = (double *) va_arg(ap, double *);
  SensibleHeat       = (double *) va_arg(ap, double *);
  ErrorString        = (char *) va_arg(ap, char *);
  
  /* print variables */
  fprintf(stderr, "%s", ErrorString);
  fprintf(stderr, "ERROR: snow_melt failed to converge to a solution in root_brent.  Variable values will be dumped to the screen, check for invalid values.\n");

  fprintf(stderr,"Dt = %f\n",Dt);
  fprintf(stderr,"Ra = %f\n",Ra);
  fprintf(stderr,"Z = %f\n",Z);
  fprintf(stderr,"Displacement = %f\n",Displacement);
  fprintf(stderr,"Z0 = %f\n",Z0);
  fprintf(stderr,"Wind = %f\n",Wind);
  fprintf(stderr,"ShortRad = %f\n",ShortRad);
  fprintf(stderr,"LongRadIn = %f\n",LongRadIn);
  fprintf(stderr,"AirDens = %f\n",AirDens);
  fprintf(stderr,"Lv = %f\n",Lv);
  fprintf(stderr,"Tair = %f\n",Tair);
  fprintf(stderr,"Press = %f\n",Press);
  fprintf(stderr,"Vpd = %f\n",Vpd);
  fprintf(stderr,"EactAir = %f\n",EactAir);
  fprintf(stderr,"Rain = %f\n",Rain);
  fprintf(stderr,"SweSurfaceLayer = %f\n",SweSurfaceLayer);
  fprintf(stderr,"SurfaceLiquidWater = %f\n",SurfaceLiquidWater);
  fprintf(stderr,"OldTSurf = %f\n",OldTSurf);
  fprintf(stderr,"RefreezeEnergy = %f\n",RefreezeEnergy[0]);
  fprintf(stderr,"VaporMassFlux = %f\n",VaporMassFlux[0]);
  fprintf(stderr,"AdvectedEnergy = %f\n",AdvectedEnergy[0]);
  fprintf(stderr,"DeltaColdContent = %f\n",DeltaColdContent[0]);
  fprintf(stderr,"TGrnd = %f\n",TGrnd);
  fprintf(stderr,"SnowDepth = %f\n",SnowDepth);
  fprintf(stderr,"SnowDensity = %f\n",SnowDensity);
  fprintf(stderr,"SurfAttenuation = %f\n",SurfAttenuation);
  fprintf(stderr,"GroundFlux = %f\n",GroundFlux[0]);
  fprintf(stderr,"LatentHeat = %f\n",LatentHeat[0]);
  fprintf(stderr,"SensibleHeat = %f\n",SensibleHeat[0]);
  
  vicerror("Finished dumping snow_melt variables.\nTry increasing SNOW_DT to get model to complete cell.\nThen check output for instabilities.\n");

  return(0.0);

}