snow_melt.c

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

/*
 * SUMMARY:      SnowMelt.c - Calculate snow accumulation and melt
 * USAGE:        
 *
 * AUTHOR:       Mark Wigmosta and Pascal Storck
 * ORG:          University of Washington, Department of Civil Engineering
 * E-MAIL:       [e-mail]
 * ORIG-DATE:     8-Oct-1996 at 08:50:06
 * LAST-MOD: Tue Mar  7 17:02:40 2000 by Keith Cherkauer <[e-mail]>
 * DESCRIPTION:  Calculate snow accumulation and melt using an energy balance
 *               approach for a two layer snow model
 * DESCRIP-END.
 * FUNCTIONS:    SnowMelt()
 * COMMENTS:     
 */

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

static char vcid[] = "$Id: snow_melt.c,v 4.1.2.2 2004/09/22 00:40:33 vicadmin Exp $";

/*****************************************************************************
  Function name: SnowMelt()

  Purpose      : Calculate snow accumulation and melt using an energy balance
                 approach for a two layer snow model

  Required     :
    double delta_t               - Model timestep (hours)
    double z2           - Reference height (m) 
    double displacement          - Displacement height (m)
    double aero_resist  - Aerodynamic resistance (uncorrected for
                                   stability) (s/m)
    double atmos->density        - Density of air (kg/m3)
    double atmos->vp             - Actual vapor pressure of air (Pa) 
    double Le           - Latent heat of vaporization (J/kg3)
    double atmos->net_short      - Net exchange of shortwave radiation (W/m2)
    double atmos->longwave       - Incoming long wave radiation (W/m2)
    double atmos->pressure       - Air pressure (Pa)
    double RainFall              - Amount of rain (m)
    double Snowfall              - Amount of snow (m)
    double atmos->air_temp       - Air temperature (C)
    double atmos->vpd            - Vapor pressure deficit (Pa)
    double wind                  - Wind speed (m/s)
    double snow->pack_water      - Liquid water content of snow pack 
    double snow->surf_water	 - Liquid water content of surface layer 
    double snow->swq             - Snow water equivalent at current pixel (m)
    double snow->vapor_flux;     - Mass flux of water vapor to or from the
                                   intercepted snow (m/time step)
    double snow->pack_temp       - Temperature of snow pack (C)
    double snow->surf_temp       - Temperature of snow pack surface layer (C)
    double snow->melt_energy     - Energy used for melting and heating of
                                   snow pack (W/m2)

  Modifies     :
    double atmos->melt           - Amount of snowpack outflow (m)
    double snow->pack_water      - Liquid water content of snow pack 
    double snow->surf_water	 - Liquid water content of surface layer 
    double snow->swq             - Snow water equivalent at current pixel (m)
    double snow->vapor_flux;     - Mass flux of water vapor to or from the
                                   intercepted snow (m)
    double snow->pack_temp       - Temperature of snow pack (C)
    double snow->surf_temp       - Temperature of snow pack surface layer (C)
    double snow->melt_energy     - Energy used for melting and heating of
                                   snow pack (W/m2)

  Comments     :
    04-Jun-04 Added descriptive error message to beginning of screen dump
	      in ErrorPrintSnowPackEnergyBalance.			TJB
    21-Sep-04 Added ErrorString to store error messages from
	      root_brent.						TJB
*****************************************************************************/
void snow_melt(soil_con_struct  *soil_con, 
	       int               rec,
	       int               iveg,
               double            z2,
               double            aero_resist,
               double            Le,
               snow_data_struct *snow, 
               double            delta_t,
               double            displacement,
               double            Z0,
               double            surf_atten,
               double            rainfall,
               double            snowfall,
               double            wind,
               double            grnd_surf_temp,
	       double            air_temp,
	       double            net_short,
	       double            longwave,
	       double            density,
	       double            pressure,
	       double            vpd,
	       double            vp,
	       double           *melt,
               double           *save_advection,
               double           *save_deltaCC,
               double           *save_grnd_flux,
               double           *save_latent,
               double           *save_sensible,
               double           *save_Qnet,
	       double           *save_refreeze_energy)
{
  int    Twidth;

  double DeltaPackCC;            /* Change in cold content of the pack */
  double DeltaPackSwq;           /* Change in snow water equivalent of the
                                   pack (m) */
  double Ice;                    /* Ice content of snow pack (m)*/
  double InitialSwq;             /* Initial snow water equivalent (m) */
  double MassBalanceError;       /* Mass balance error (m) */
  double MaxLiquidWater;         /* Maximum liquid water content of pack (m) */
  double OldTSurf;               /* Old snow surface temperature (C) */
  double PackCC;                 /* Cold content of snow pack (J) */
  double PackSwq;                /* Snow pack snow water equivalent (m) */
  double Qnet;                   /* Net energy exchange at the surface (W/m2) */
  double RefreezeEnergy;         /* refreeze energy (W/m2) */
  double RefrozenWater;          /* Amount of refrozen water (m) */
  double SnowFallCC;             /* Cold content of new snowfall (J) */
  double SnowMelt;               /* Amount of snow melt during time interval
                                   (m water equivalent) */
  double SurfaceCC;              /* Cold content of snow pack (J) */
  double SurfaceSwq;             /* Surface layer snow water equivalent (m) */
  double SnowFall;
  double RainFall;
  double vapor_flux;
  double advection;
  double deltaCC;
  double grnd_flux;		/* thermal flux through snowpack from ground */
  double latent_heat;		/* thermal flux through snowpack from ground */
  double sensible_heat;		/* thermal flux through snowpack from ground */
  double melt_energy = 0.;
  char ErrorString[MAXSTRING];

  SnowFall = snowfall / 1000.; /* convet to m */
  RainFall = rainfall / 1000.; /* convet to m */

  InitialSwq = snow->swq;
  OldTSurf = snow->surf_temp;

  /* Initialize snowpack variables */
  
  Ice  = snow->swq - snow->pack_water - snow->surf_water;
  
  /* Reconstruct snow pack */
  if (Ice > MAX_SURFACE_SWE)
    SurfaceSwq = MAX_SURFACE_SWE;
  else
    SurfaceSwq = Ice;
  PackSwq = Ice - SurfaceSwq;
  
  /* Calculate cold contents */
  SurfaceCC = CH_ICE * SurfaceSwq * snow->surf_temp;
  PackCC = CH_ICE * PackSwq * snow->pack_temp;
  if (air_temp > 0.0)
    SnowFallCC = 0.0;
  else
    SnowFallCC = CH_ICE * SnowFall * air_temp;
  
  /* Distribute fresh snowfall */
  if (SnowFall > (MAX_SURFACE_SWE - SurfaceSwq)) {
    DeltaPackSwq = SurfaceSwq + SnowFall - MAX_SURFACE_SWE;
    if (DeltaPackSwq > SurfaceSwq)
      DeltaPackCC = SurfaceCC + (SnowFall - MAX_SURFACE_SWE)/SnowFall *
        SnowFallCC;
    else
      DeltaPackCC = DeltaPackSwq/SurfaceSwq * SurfaceCC;
    SurfaceSwq = MAX_SURFACE_SWE;
    SurfaceCC += SnowFallCC - DeltaPackCC;
    PackSwq += DeltaPackSwq;
    PackCC += DeltaPackCC;
  }
  else {
    SurfaceSwq += SnowFall;
    SurfaceCC += SnowFallCC;
  }
  if (SurfaceSwq > 0.0)
    snow->surf_temp = SurfaceCC/(CH_ICE * SurfaceSwq);
  else 
    snow->surf_temp = 0.0;
  if (PackSwq > 0.0)    
    snow->pack_temp = PackCC/(CH_ICE * PackSwq);
  else
    snow->pack_temp = 0.0;

  /* Adjust ice and snow->surf_water */
  Ice += SnowFall;
  snow->surf_water += RainFall;
  
  /* Calculate the surface energy balance for snow_temp = 0.0 */
  
  vapor_flux = snow->vapor_flux;

  Qnet = CalcSnowPackEnergyBalance((double)0.0, delta_t, aero_resist,
				   z2, displacement, Z0, wind, net_short, 
				   longwave, density, Le, air_temp,
				   pressure * 1000., vpd * 1000., vp * 1000.,
				   RainFall, SurfaceSwq, snow->surf_water, 
				   OldTSurf, &RefreezeEnergy, &vapor_flux, 
				   &advection, &deltaCC, grnd_surf_temp,
				   snow->depth, snow->density, surf_atten, 
				   &grnd_flux, &latent_heat, &sensible_heat);

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

  /* If Qnet == 0.0, then set the surface temperature to 0.0 */
  if (Qnet == 0.0) {
    snow->surf_temp = 0.0;
    if (RefreezeEnergy >= 0.0) {
      RefrozenWater = RefreezeEnergy/(Lf * RHO_W) 
          * delta_t * SECPHOUR; 
      if (RefrozenWater > snow->surf_water) {
        RefrozenWater = snow->surf_water;
        RefreezeEnergy = RefrozenWater * Lf * RHO_W/
          (delta_t * SECPHOUR);
      } 
      melt_energy  += RefreezeEnergy;
      SurfaceSwq   += RefrozenWater;
      Ice          += RefrozenWater;
      snow->surf_water   -= RefrozenWater;
      assert(snow->surf_water >= 0.0);
      SnowMelt      = 0.0;
    }
    else {
      
      /* Calculate snow melt */      
      SnowMelt = fabs(RefreezeEnergy)/(Lf * RHO_W) * 
        delta_t * SECPHOUR;
      melt_energy += RefreezeEnergy;
    }

    /* Convert vapor mass flux to a depth per timestep and adjust snow->surf_water */
    snow->vapor_flux *= delta_t * SECPHOUR;
    
    if (snow->surf_water < -(snow->vapor_flux)) {
      snow->vapor_flux = -(snow->surf_water);
      snow->surf_water    = 0.0;
    }
    else
      snow->surf_water += snow->vapor_flux;
    
    /* If SnowMelt < Ice, there was incomplete melting of the pack */
    
    if (SnowMelt < Ice) {
      if (SnowMelt <= PackSwq) {
        snow->surf_water += SnowMelt;
        PackSwq    -= SnowMelt;
        Ice        -= SnowMelt;
      }
      else {
        snow->surf_water += SnowMelt + snow->pack_water;
        snow->pack_water  = 0.0;
        PackSwq     = 0.0;
        Ice        -= SnowMelt;
        SurfaceSwq  = Ice;
      }
    }
    
    /* Else, SnowMelt > Ice and there was complete melting of the pack */
    else {
      SnowMelt    = Ice;
      snow->surf_water += Ice;
      SurfaceSwq  = 0.0;
      snow->surf_temp      = 0.0;
      PackSwq     = 0.0;
      snow->pack_temp      = 0.0;
      Ice         = 0.0;
    }
  }
  
  /* Else, SnowPackEnergyBalance(T=0.0) <= 0.0 */
  else  {
    /* Calculate surface layer temperature using "Brent method" */

    vapor_flux = snow->vapor_flux;

    snow->surf_temp = root_brent((double)(snow->surf_temp-SNOW_DT), 
				 (double)0.0, ErrorString,
				 SnowPackEnergyBalance, delta_t, 
				 aero_resist, z2, 
				 displacement, Z0, wind, net_short, longwave,
				 density, Le, air_temp, pressure * 1000.,
				 vpd * 1000., vp * 1000., RainFall, 
				 SurfaceSwq,
				 snow->surf_water, OldTSurf, &RefreezeEnergy, 
				 &vapor_flux, &advection, &deltaCC, 
				 grnd_surf_temp, snow->depth,
				 snow->density,surf_atten,&grnd_flux,
				 &latent_heat, &sensible_heat);

    if(snow->surf_temp <= -9998)
      ErrorSnowPackEnergyBalance(snow->surf_temp, delta_t, aero_resist,
				 z2, displacement, Z0, wind, net_short,
				 longwave, density, Le, air_temp,
				 pressure * 1000., vpd * 1000., vp * 1000.,
				 RainFall, SurfaceSwq, snow->surf_water, 
				 OldTSurf,
				 &RefreezeEnergy, &vapor_flux, &advection, 
				 &deltaCC, grnd_surf_temp,
				 snow->depth, snow->density,surf_atten,
				 &grnd_flux,&latent_heat,
				 &sensible_heat, ErrorString);

    Qnet = CalcSnowPackEnergyBalance(snow->surf_temp, delta_t, aero_resist,
				     z2, displacement, Z0, wind, net_short,
				     longwave, density, Le, air_temp,
				     pressure * 1000., vpd * 1000., 
				     vp * 1000.,
				     RainFall, SurfaceSwq, snow->surf_water, 
				     OldTSurf,
				     &RefreezeEnergy, &vapor_flux, 
				     &advection, &deltaCC, grnd_surf_temp,