Skip to content

Instantly share code, notes, and snippets.

@jhamman

jhamman/0.Summary (RACM.33-fixed)

Last active August 29, 2015 14:23
Show Gist options
  • Save jhamman/97b69ec6084270723077 to your computer and use it in GitHub Desktop.
Save jhamman/97b69ec6084270723077 to your computer and use it in GitHub Desktop.
Download ZIP
Fixed Atmospheric Decoupling VIC Changes in RACM.33
Raw
0.Summary (RACM.33-fixed)
Fixing the land-atmosphere coupling in RACM 33
========
To "fix" the surface layer in RACM. We need to add an upper stability limit to the aerodynamic resistance calculation in both `Fixed calc_surf_energy_bal.c` and `SnowPackEnergyBalance`.
These changes, along with the addition of the constant limit in `vicNl_def.h` are shown here.
Raw
calc_surf_energy_bal.c
#include <stdio.h>
#include <stdlib.h>
#include "vicNl.h"
static char vcid[] =
"$Id: calc_surf_energy_bal.c,v 4.2 2000/08/16 00:35:28 vicadmin Exp vicadmin $";
/*******************************************************************************
calc_surf_energy_bal.c Greg O'Donnell and Keith Cherkauer Sept 9 1997
This function calculates the surface temperature, in the
case of no snow cover. Evaporation is computed using the
previous ground heat flux, and then used to comput latent heat
in the energy balance routine. Surface temperature is found
using the Root Brent method (Numerical Recipies).
modifications:
02-29-00 Included variables needed to compute energy flux
through the snow pack. The ground surface energy
balance will now be a mixture of snow covered
and bare ground, controlled by the snow cover
fraction set in solve_snow.c KAC
modifications:
12/2011 - code reformatting and cleanup
BN
10/2012 - ported changes for aero_resist_used from VIC 4.0.6 to store
the aerodynamic resistance actually used in flux
calculations.
BN
*******************************************************************************/
double
calc_surf_energy_bal(int rec,
int iveg,
int nlayer,
int Nveg,
float dt,
int Nnodes,
int veg_class,
int band,
int hour,
double ice0,
double moist,
double dp,
double surf_atten,
double T0,
double shortwave,
double longwave,
double air_temp,
double Le,
double Ls,
double mu,
double displacement,
double roughness,
double ref_height,
double snow_energy,
double vapor_flux,
double bare_albedo,
double *aero_resist,
double *aero_resist_used,
double *wind,
double *rainfall,
double *ppt,
float *root,
atmos_data_struct *atmos,
veg_var_struct *veg_var_wet,
veg_var_struct *veg_var_dry,
energy_bal_struct *energy,
snow_data_struct *snow,
layer_data_struct *layer_wet,
layer_data_struct *layer_dry,
soil_con_struct *soil_con,
dmy_struct *dmy)
{
extern veg_lib_struct *veg_lib;
extern option_struct options;
char FIRST_SOLN[1];
double T2;
double Ts_old;
double T1_old;
double Tair;
double atmos_density;
double pz; /* kpa, air pressure, M.Huang*/
double albedo;
double emissivity;
double kappa1;
double kappa2;
double Cs1;
double Cs2;
double D1;
double D2;
double delta_t;
double Vapor;
double max_moist;
double bubble;
double expt;
double T1;
double snow_depth;
double snow_density;
double Tsnow_surf;
double snow_cover_fraction;
double snow_flux;
int VEG;
double Tsurf;
double U;
double error;
double Wdew[2];
double *T_node;
double Tnew_node[MAX_NODES];
double *dz_node;
double *kappa_node;
double *Cs_node;
double *moist_node;
double *bubble_node;
double *expt_node;
double *max_moist_node;
double *ice_node;
double *alpha;
double *beta;
double *gamma;
layer_data_struct layer[MAX_LAYERS];
double T_lower, T_upper, Terr;
/**************************************************
Correct Aerodynamic Resistance for Stability
**************************************************/
U = wind[0];
/* RACMNOTE: Uncommented. Note that the HUGE_RESIST line should not be
triggered since there is already a minimum wind speed - BN - 03/2012 */
if (U > 0.0)
*aero_resist_used = aero_resist[0]/
StabilityCorrection(ref_height, displacement, T0,
air_temp, U, roughness);
/* RASMnote - added to avoid decoupling under stable conditions
- BN - 01/2013 */
if (air_temp > T0 && *aero_resist_used > RA_STABLE_UPPER_LIMIT) {
*aero_resist_used = RA_STABLE_UPPER_LIMIT;
}
}
else
*aero_resist_used = HUGE_RESIST;
/**************************************************
Compute Evaporation and Transpiration
**************************************************/
if (iveg != Nveg) {
if (veg_lib[veg_class].LAI[dmy->month - 1] > 0.0) {
VEG = TRUE;
}
else {
VEG = FALSE;
}
}
else {
VEG = FALSE;
}
Vapor = vapor_flux / (double)dt / 3600.;
/**************************************************
Set All Variables For Use
**************************************************/
T2 = energy->T[Nnodes - 1];
Ts_old = energy->T[0];
T1_old = energy->T[1];
Tair = air_temp;
atmos_density = atmos->density[hour];
pz = atmos->pressure[hour]; /* M.Huang*/
albedo = bare_albedo;
emissivity = 1.;
kappa1 = energy->kappa[0];
kappa2 = energy->kappa[1];
Cs1 = energy->Cs[0];
Cs2 = energy->Cs[1];
D1 = soil_con->depth[0];
D2 = soil_con->depth[0];
delta_t = (double)dt * 3600.;
max_moist = soil_con->max_moist[0] / (soil_con->depth[0] * 1000.);
bubble = soil_con->bubble[0];
expt = soil_con->expt[0];
snow_depth = snow->depth;
snow_density = snow->density;
Tsnow_surf = snow->surf_temp;
snow_cover_fraction = snow->coverage;
Wdew[WET] = veg_var_wet->Wdew;
if (options.DIST_PRCP) {
Wdew[DRY] = veg_var_dry->Wdew;
}
FIRST_SOLN[0] = TRUE;
/*************************************************************
Prepare soil node variables for finite difference solution
*************************************************************/
if (!options.QUICK_FLUX) {
bubble_node = soil_con->bubble_node;
expt_node = soil_con->expt_node;
max_moist_node = soil_con->max_moist_node;
alpha = soil_con->alpha;
beta = soil_con->beta;
gamma = soil_con->gamma;
moist_node = energy->moist;
kappa_node = energy->kappa_node;
Cs_node = energy->Cs_node;
T_node = energy->T;
dz_node = soil_con->dz_node;
ice_node = energy->ice;
}
else {
bubble_node = NULL;
expt_node = NULL;
max_moist_node = NULL;
alpha = NULL;
beta = NULL;
gamma = NULL;
moist_node = NULL;
kappa_node = NULL;
Cs_node = NULL;
T_node = NULL;
dz_node = NULL;
ice_node = NULL;
}
/**************************************************
Find Surface Temperature Using Root Brent Method
**************************************************/
if (options.FULL_ENERGY) {
/** Added for temporary backwards compatability **/
if (snow->swq > 0) {
T_lower = 0.;
T_upper = energy->T[0] - SURF_DT;
}
else {
T_lower = 0.5 * (energy->T[0] + air_temp) - SURF_DT;
T_upper = 0.5 * (energy->T[0] + air_temp) + SURF_DT;
}
if (T_lower > T_upper) {
Terr = T_lower;
T_lower = T_upper;
T_upper = Terr;
}
#if QUICK_FS
Tsurf = root_brent(T_lower, T_upper,
func_surf_energy_bal, T2, Ts_old, T1_old, Tair,
*aero_resist_used, atmos_density, pz, shortwave,
longwave, albedo, emissivity, kappa1, kappa2,
Cs1, Cs2, D1, D2, dp, delta_t, Le, Ls, Vapor, moist,
ice0, max_moist, bubble, expt, snow_depth,
snow_density, Tsnow_surf, snow_cover_fraction,
surf_atten, U, displacement, roughness, ref_height,
(double)soil_con->elevation, soil_con->b_infilt,
soil_con->max_infil, (double)dt, atmos->vpd[hour],
snow_energy, mu, rainfall, Wdew, &energy->grnd_flux,
&T1, &energy->latent, &energy->sensible,
&energy->deltaH, &energy->snow_flux, &energy->error,
soil_con->depth, soil_con->Wcr, soil_con->Wpwp,
soil_con->resid_moist, T_node, Tnew_node, dz_node,
kappa_node, Cs_node, moist_node, bubble_node,
expt_node, max_moist_node, ice_node, alpha, beta,
gamma, soil_con->ufwc_table_layer[0],
soil_con->ufwc_table_node, root, layer_wet,
layer_dry, veg_var_wet, veg_var_dry, VEG,
veg_class, dmy->month, Nnodes,
FIRST_SOLN, snow->snow, soil_con->FS_ACTIVE);
#else
Tsurf = root_brent(T_lower, T_upper,
func_surf_energy_bal, T2, Ts_old, T1_old, Tair,
*aero_resist_used, atmos_density, pz, shortwave,
longwave, albedo, emissivity, kappa1, kappa2,
Cs1, Cs2, D1, D2, dp, delta_t, Le, Ls, Vapor, moist,
ice0, max_moist, bubble, expt, snow_depth,
snow_density, Tsnow_surf, snow_cover_fraction,
surf_atten, U, displacement, roughness, ref_height,
(double)soil_con->elevation, soil_con->b_infilt,
soil_con->max_infil, (double)dt, atmos->vpd[hour],
snow_energy, mu, rainfall, Wdew, &energy->grnd_flux,
&T1, &energy->latent, &energy->sensible,
&energy->deltaH, &energy->snow_flux, &energy->error,
soil_con->depth, soil_con->Wcr, soil_con->Wpwp,
soil_con->resid_moist, T_node, Tnew_node, dz_node,
kappa_node, Cs_node, moist_node, bubble_node,
expt_node, max_moist_node, ice_node, alpha, beta,
gamma, root, layer_wet, layer_dry, veg_var_wet,
veg_var_dry, VEG, veg_class,
dmy->month, Nnodes, FIRST_SOLN, snow->snow,
soil_con->FS_ACTIVE);
#endif
if (Tsurf <= -9998) {
#if QUICK_FS
error = error_calc_surf_energy_bal(Tsurf, T2, Ts_old, T1_old, Tair,
*aero_resist_used, atmos_density,
pz, shortwave,
longwave, albedo, emissivity,
kappa1, kappa2, Cs1, Cs2, D1, D2,
dp, delta_t, Le, Ls, Vapor,
moist, ice0, max_moist, bubble,
expt, snow_depth, snow_density,
Tsnow_surf, snow_cover_fraction,
surf_atten, U, displacement,
roughness, ref_height,
(double)soil_con->elevation,
soil_con->b_infilt,
soil_con->max_infil, (double)dt,
atmos->vpd[hour], snow_energy,
mu, rainfall, Wdew,
&energy->grnd_flux,
&T1, &energy->latent,
&energy->sensible,
&energy->deltaH,
&energy->snow_flux,
&energy->error, soil_con->depth,
soil_con->Wcr, soil_con->Wpwp,
soil_con->resid_moist, T_node,
Tnew_node, dz_node, kappa_node,
Cs_node, moist_node, bubble_node,
expt_node, max_moist_node,
ice_node, alpha, beta, gamma,
soil_con->ufwc_table_layer[0],
soil_con->ufwc_table_node,
root, layer_wet, layer_dry,
veg_var_wet, veg_var_dry, VEG,
veg_class, dmy->month, iveg,
Nnodes, FIRST_SOLN, snow->snow,
soil_con->FS_ACTIVE);
#else
error = error_calc_surf_energy_bal(Tsurf, T2, Ts_old, T1_old, Tair,
*aero_resist_used, atmos_density,
pz, shortwave,
longwave, albedo, emissivity,
kappa1, kappa2, Cs1, Cs2, D1, D2,
dp, delta_t, Le, Ls, Vapor,
moist, ice0, max_moist, bubble,
expt, snow_depth, snow_density,
Tsnow_surf, snow_cover_fraction,
surf_atten, U, displacement,
roughness, ref_height,
(double)soil_con->elevation,
soil_con->b_infilt,
soil_con->max_infil, (double)dt,
atmos->vpd[hour], snow_energy,
mu, rainfall, Wdew,
&energy->grnd_flux,
&T1, &energy->latent,
&energy->sensible,
&energy->deltaH,
&energy->snow_flux,
&energy->error, soil_con->depth,
soil_con->Wcr, soil_con->Wpwp,
soil_con->resid_moist, T_node,
Tnew_node, dz_node, kappa_node,
Cs_node, moist_node, bubble_node,
expt_node, max_moist_node,
ice_node, alpha, beta, gamma,
root, layer_wet, layer_dry,
veg_var_wet, veg_var_dry, VEG,
veg_class, dmy->month, iveg,
Nnodes, FIRST_SOLN, snow->snow,
soil_con->FS_ACTIVE);
#endif
}
}
else {
/** Frozen soil model run with no surface energy balance **/
Tsurf = Tair;
}
/**************************************************
Recalculate Energy Balance Terms for Final Surface Temperature
**************************************************/
#if QUICK_FS
error = solve_surf_energy_bal(Tsurf, T2, Ts_old, T1_old, Tair,
*aero_resist_used,
atmos_density, pz, shortwave, longwave,
albedo, emissivity, kappa1, kappa2, Cs1, Cs2,
D1, D2, dp, delta_t, Le, Ls, Vapor, moist,
ice0, max_moist, bubble, expt, snow_depth,
snow_density, Tsnow_surf,
snow_cover_fraction, surf_atten, U,
displacement, roughness, ref_height,
(double)soil_con->elevation,
soil_con->b_infilt, soil_con->max_infil,
(double)dt, atmos->vpd[hour], snow_energy, mu,
rainfall, Wdew, &energy->grnd_flux,
&T1, &energy->latent, &energy->sensible,
&energy->deltaH, &energy->snow_flux,
&energy->error, soil_con->depth,
soil_con->Wcr, soil_con->Wpwp,
soil_con->resid_moist, T_node,
Tnew_node, dz_node, kappa_node, Cs_node,
moist_node, bubble_node, expt_node,
max_moist_node, ice_node, alpha, beta, gamma,
soil_con->ufwc_table_layer[0],
soil_con->ufwc_table_node, root, layer_wet,
layer_dry, veg_var_wet, veg_var_dry, VEG,
veg_class, dmy->month, Nnodes,
FIRST_SOLN, snow->snow,
soil_con->FS_ACTIVE);
#else
error = solve_surf_energy_bal(Tsurf, T2, Ts_old, T1_old, Tair,
*aero_resist_used,
atmos_density, pz, shortwave, longwave,
albedo, emissivity, kappa1, kappa2, Cs1, Cs2,
D1, D2, dp, delta_t, Le, Ls, Vapor, moist,
ice0, max_moist, bubble, expt, snow_depth,
snow_density, Tsnow_surf,
snow_cover_fraction, surf_atten, U,
displacement, roughness, ref_height,
(double)soil_con->elevation,
soil_con->b_infilt, soil_con->max_infil,
(double)dt, atmos->vpd[hour], snow_energy, mu,
rainfall, Wdew, &energy->grnd_flux,
&T1, &energy->latent, &energy->sensible,
&energy->deltaH, &energy->snow_flux,
&energy->error, soil_con->depth,
soil_con->Wcr, soil_con->Wpwp,
soil_con->resid_moist, T_node,
Tnew_node, dz_node, kappa_node, Cs_node,
moist_node, bubble_node, expt_node,
max_moist_node, ice_node, alpha, beta, gamma,
root, layer_wet, layer_dry, veg_var_wet,
veg_var_dry, VEG,
veg_class, dmy->month, Nnodes, FIRST_SOLN,
snow->snow,
soil_con->FS_ACTIVE);
#endif
energy->error = error;
/***************************************************
Recalculate Soil Moisture and Thermal Properties
***************************************************/
if (options.GRND_FLUX) {
if (options.QUICK_FLUX) {
energy->T[0] = Tsurf;
energy->T[1] = T1;
}
else {
finish_frozen_soil_calcs(energy, layer_wet, layer_dry, layer,
soil_con,
Nnodes, iveg, mu, Tnew_node, kappa_node,
Cs_node,
moist_node);
}
}
else {
energy->T[0] = Tsurf;
}
/** Store precipitation that reaches the surface */
if (!snow->snow) {
if (iveg != Nveg) {
if (veg_lib[veg_class].LAI[dmy->month - 1] <= 0.0) {
veg_var_wet->throughfall = rainfall[WET];
if (options.DIST_PRCP) {
veg_var_dry->throughfall = rainfall[DRY];
}
ppt[WET] = veg_var_wet->throughfall;
if (options.DIST_PRCP) {
ppt[DRY] = veg_var_dry->throughfall;
}
}
else {
ppt[WET] = veg_var_wet->throughfall;
if (options.DIST_PRCP) {
ppt[DRY] = veg_var_dry->throughfall;
}
}
}
else {
ppt[WET] = rainfall[WET];
if (options.DIST_PRCP) {
ppt[DRY] = rainfall[DRY];
}
}
}
/** Store net longwave radiation **/
/* RACMnote: hour == 0 in RASM, so the if {} statement is always executed and
the value from which the outgoing is substracted is the value that is
stored as bare_energy->longwave in surface_fluxes () -- BN -- 06/2012 */
if (hour < 24) {
energy->longwave = longwave -
STEFAN_B *
(Tsurf + KELVIN) *
(Tsurf + KELVIN) *
(Tsurf + KELVIN) *
(Tsurf + KELVIN);
}
else {
energy->longwave = longwave;
}
/** Store surface albedo **/
energy->albedo = bare_albedo;
/** Store net shortwave radiation **/
energy->shortwave = (1. - energy->albedo) * shortwave;
/** Return soil surface temperature **/
return (Tsurf);
}
double
solve_surf_energy_bal(double Tsurf,
...)
{
va_list ap;
double error;
va_start(ap, Tsurf);
error = func_surf_energy_bal(Tsurf, ap);
va_end(ap);
return error;
}
double
error_calc_surf_energy_bal(double Tsurf,
...)
{
va_list ap;
double error;
va_start(ap, Tsurf);
error = error_print_surf_energy_bal(Tsurf, ap);
va_end(ap);
return error;
}
double
error_print_surf_energy_bal(double Ts,
va_list ap)
{
extern option_struct options;
/** Thermal Properties **/
double T2;
double Ts_old;
double T1_old;
double Tair;
double ra;
double atmos_density;
double pz; /* kpa, air pressure, M.Huang*/
double shortwave;
double longwave;
double albedo;
double emissivity;
double kappa1;
double kappa2;
double Cs1;
double Cs2;
double D1;
double D2;
double dp;
double delta_t;
double Le;
double Ls;
double Vapor;
double moist;
double ice0;
double max_moist;
double bubble;
double expt;
double snow_depth;
double snow_density;
double Tsnow_surf;
double snow_cover_fraction;
double surf_atten;
double wind;
double displacement;
double roughness;
double ref_height;
double elevation;
double b_infilt;
double max_infil;
double dt;
double vpd;
double snow_energy;
double mu;
double *rainfall;
double *Wdew;
double *grnd_flux;
double *T1;
double *latent_heat;
double *sensible_heat;
double *deltaH;
double *snow_flux;
double *store_error;
double *depth;
double *Wcr;
double *Wpwp;
double *resid_moist;
double *T_node;
double *Tnew_node;
double *dz_node;
double *kappa_node;
double *Cs_node;
double *moist_node;
double *bubble_node;
double *expt_node;
double *max_moist_node;
double *ice_node;
double *alpha;
double *beta;
double *gamma;
#if QUICK_FS
double **ufwc_table_layer;
double ***ufwc_table_node;
#endif
float *root;
layer_data_struct *layer_wet;
layer_data_struct *layer_dry;
veg_var_struct *veg_var_wet;
veg_var_struct *veg_var_dry;
int VEG;
int veg_class;
int month;
int iveg;
int Nnodes;
char *FIRST_SOLN;
int SNOWING;
int FS_ACTIVE;
int i;
/* Initialize Variables */
T2 = (double) va_arg(ap, double);
Ts_old = (double) va_arg(ap, double);
T1_old = (double) va_arg(ap, double);
Tair = (double) va_arg(ap, double);
ra = (double) va_arg(ap, double);
atmos_density = (double) va_arg(ap, double);
pz = (double) va_arg(ap, double); /* M.Huang*/
shortwave = (double) va_arg(ap, double);
longwave = (double) va_arg(ap, double);
albedo = (double) va_arg(ap, double);
emissivity = (double) va_arg(ap, double);
kappa1 = (double) va_arg(ap, double);
kappa2 = (double) va_arg(ap, double);
Cs1 = (double) va_arg(ap, double);
Cs2 = (double) va_arg(ap, double);
D1 = (double) va_arg(ap, double);
D2 = (double) va_arg(ap, double);
dp = (double) va_arg(ap, double);
delta_t = (double) va_arg(ap, double);
Le = (double) va_arg(ap, double);
Ls = (double) va_arg(ap, double);
Vapor = (double) va_arg(ap, double);
moist = (double) va_arg(ap, double);
ice0 = (double) va_arg(ap, double);
max_moist = (double) va_arg(ap, double);
bubble = (double) va_arg(ap, double);
expt = (double) va_arg(ap, double);
snow_depth = (double) va_arg(ap, double);
snow_density = (double) va_arg(ap, double);
Tsnow_surf = (double) va_arg(ap, double);
snow_cover_fraction = (double) va_arg(ap, double);
surf_atten = (double) va_arg(ap, double);
wind = (double) va_arg(ap, double);
displacement = (double) va_arg(ap, double);
roughness = (double) va_arg(ap, double);
ref_height = (double) va_arg(ap, double);
elevation = (double) va_arg(ap, double);
b_infilt = (double) va_arg(ap, double);
max_infil = (double) va_arg(ap, double);
dt = (double) va_arg(ap, double);
vpd = (double) va_arg(ap, double);
snow_energy = (double) va_arg(ap, double);
mu = (double) va_arg(ap, double);
rainfall = (double *) va_arg(ap, double *);
Wdew = (double *) va_arg(ap, double *);
grnd_flux = (double *) va_arg(ap, double *);
T1 = (double *) va_arg(ap, double *);
latent_heat = (double *) va_arg(ap, double *);
sensible_heat = (double *) va_arg(ap, double *);
deltaH = (double *) va_arg(ap, double *);
snow_flux = (double *) va_arg(ap, double *);
store_error = (double *) va_arg(ap, double *);
depth = (double *) va_arg(ap, double *);
Wcr = (double *) va_arg(ap, double *);
Wpwp = (double *) va_arg(ap, double *);
resid_moist = (double *) va_arg(ap, double *);
T_node = (double *) va_arg(ap, double *);
Tnew_node = (double *) va_arg(ap, double *);
dz_node = (double *) va_arg(ap, double *);
kappa_node = (double *) va_arg(ap, double *);
Cs_node = (double *) va_arg(ap, double *);
moist_node = (double *) va_arg(ap, double *);
bubble_node = (double *) va_arg(ap, double *);
expt_node = (double *) va_arg(ap, double *);
max_moist_node = (double *) va_arg(ap, double *);
ice_node = (double *) va_arg(ap, double *);
alpha = (double *) va_arg(ap, double *);
beta = (double *) va_arg(ap, double *);
gamma = (double *) va_arg(ap, double *);
#if QUICK_FS
ufwc_table_layer = (double **) va_arg(ap, double **);
ufwc_table_node = (double ***) va_arg(ap, double ***);
#endif
root = (float *) va_arg(ap, float *);
layer_wet = (layer_data_struct *) va_arg(ap, layer_data_struct *);
layer_dry = (layer_data_struct *) va_arg(ap, layer_data_struct *);
veg_var_wet = (veg_var_struct *) va_arg(ap, veg_var_struct *);
veg_var_dry = (veg_var_struct *) va_arg(ap, veg_var_struct *);
VEG = (int) va_arg(ap, int);
veg_class = (int) va_arg(ap, int);
month = (int) va_arg(ap, int);
iveg = (int) va_arg(ap, int);
Nnodes = (int) va_arg(ap, int);
FIRST_SOLN = (char *) va_arg(ap, char *);
SNOWING = (int) va_arg(ap, int);
FS_ACTIVE = (int) va_arg(ap, int);
/* Print Variables */
fprintf(stderr, "T2 = %f\n", T2);
fprintf(stderr, "Ts_old = %f\n", Ts_old);
fprintf(stderr, "T1_old = %f\n", T1_old);
fprintf(stderr, "Tair = %f\n", Tair);
fprintf(stderr, "ra = %f\n", ra);
fprintf(stderr, "atmos_density = %f\n", atmos_density);
fprintf(stderr, "air pressure = %f\n", pz); /* M.Huang*/
fprintf(stderr, "shortwave = %f\n", shortwave);
fprintf(stderr, "longwave = %f\n", longwave);
fprintf(stderr, "albedo = %f\n", albedo);
fprintf(stderr, "emissivity = %f\n", emissivity);
fprintf(stderr, "kappa1 = %f\n", kappa1);
fprintf(stderr, "kappa2 = %f\n", kappa2);
fprintf(stderr, "Cs1 = %f\n", Cs1);
fprintf(stderr, "Cs2 = %f\n", Cs2);
fprintf(stderr, "D1 = %f\n", D1);
fprintf(stderr, "D2 = %f\n", D2);
fprintf(stderr, "dp = %f\n", dp);
fprintf(stderr, "delta_t = %f\n", delta_t);
fprintf(stderr, "Le = %f\n", Le);
fprintf(stderr, "Ls = %f\n", Ls);
fprintf(stderr, "Vapor = %f\n", Vapor);
fprintf(stderr, "moist = %f\n", moist);
fprintf(stderr, "ice0 = %f\n", ice0);
fprintf(stderr, "max_moist = %f\n", max_moist);
fprintf(stderr, "bubble = %f\n", bubble);
fprintf(stderr, "expt = %f\n", expt);
fprintf(stderr, "surf_atten = %f\n", surf_atten);
fprintf(stderr, "wind = %f\n", wind);
fprintf(stderr, "displacement = %f\n", displacement);
fprintf(stderr, "roughness = %f\n", roughness);
fprintf(stderr, "ref_height = %f\n", ref_height);
fprintf(stderr, "elevation = %f\n", elevation);
fprintf(stderr, "b_infilt = %f\n", b_infilt);
fprintf(stderr, "max_infil = %f\n", max_infil);
fprintf(stderr, "dt = %f\n", dt);
fprintf(stderr, "vpd = %f\n", vpd);
fprintf(stderr, "snow_energy = %f\n", snow_energy);
fprintf(stderr, "mu = %f\n", mu);
fprintf(stderr, "rainfall = %f\n", rainfall[0]);
fprintf(stderr, "Wdew = %f\n", Wdew[0]);
fprintf(stderr, "grnd_flux = %f\n", grnd_flux[0]);
fprintf(stderr, "T1 = %f\n", T1[0]);
fprintf(stderr, "latent_heat = %f\n", latent_heat[0]);
fprintf(stderr, "sensible_heat = %f\n", sensible_heat[0]);
fprintf(stderr, "deltaH = %f\n", deltaH[0]);
fprintf(stderr, "store_error = %f\n", store_error[0]);
fprintf(stderr, "depth = %f\n", depth[0]);
fprintf(stderr, "Wcr = %f\n", Wcr[0]);
fprintf(stderr, "Wpwp = %f\n", Wpwp[0]);
fprintf(stderr, "Residual Moisture = %f\n", resid_moist[0]);
fprintf(stderr, "VEG = %i\n", VEG);
fprintf(stderr, "veg_class = %i\n", veg_class);
fprintf(stderr, "month = %i\n", month);
write_layer(layer_wet, iveg, options.Nlayer, depth);
if (options.DIST_PRCP) {
write_layer(layer_dry, iveg, options.Nlayer, depth);
}
write_vegvar(&(veg_var_wet[0]), iveg);
if (options.DIST_PRCP) {
write_vegvar(&(veg_var_dry[0]), iveg);
}
if (!options.QUICK_FLUX) {
fprintf(
stderr,
"Node\tT\tTnew\tdz\tkappa\tCs\tmoist\tbubble\texpt\tmax_moist\tice\n");
for (i = 0; i < Nnodes; i++) {
fprintf(
stderr,
"%i\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\n",
i, T_node[i], Tnew_node[i], dz_node[i], kappa_node[i],
Cs_node[i], moist_node[i], bubble_node[i], expt_node[i],
max_moist_node[i], ice_node[i]);
}
}
vicerror("Finished writing calc_surf_energy_bal variables.\nTry increasing SURF_DT to get model to complete cell.\nThen check output for instabilities.");
return(0.0);
}
Raw
SnowPackEnergyBalance.c
/*
* SUMMARY: SnowPackEnergyBalance.c - Calculate snow pack energy balance
* USAGE: Part of DHSVM
*
* AUTHOR: Bart Nijssen
* ORG: University of Washington, Department of Civil Engineering
* E-MAIL: nijssen@u.washington.edu
* ORIG-DATE: 8-Oct-1996 at 09:09:29
* LAST-MOD: Thu Apr 6 12:43:30 2000 by Keith Cherkauer <cherkaue@u.washington.edu>
* DESCRIPTION: Calculate snow pack energy balance
* DESCRIP-END.
* FUNCTIONS: SnowPackEnergyBalance()
* COMMENTS:
*/
#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <vicNl.h>
/* RASMnote: Collect these #defines globally to ensure consistency - BN -
06/2012 */
#define GRAMSPKG 1000.
/* RASMnote: Should this be 4.1868 or is this already taking into account the
rest of the units in the equation in which it is used? This is the same in
VIC 4.1.2, so I'll leave it be - BN - 06/2012 */
#define CH_WATER 4186.8e3
#define JOULESPCAL 4.1868 /* Joules per calorie */
#define EPS 0.622 /* ratio of molecular weight of water vapor to
that for dry air */
#define CP 1013.0 /* Specific heat of moist air at constant
pressure (J/(kg*C)) */
/*******************************************************************************
Function name: SnowPackEnergyBalance()
Purpose : Calculate the surface energy balance for the snow pack
Required :
double TSurf - new estimate of effective surface temperature
va_list ap - Argument list initialized by va_start(). For
elements of list and order, see beginning of
routine
Returns :
double RestTerm - Rest term in the energy balance
Modifies :
double *RefreezeEnergy - Refreeze energy (W/m2)
double *VaporMassFlux - Mass flux of water vapor to or from the
intercepted snow (m/s)
Comments :
Reference: Bras, R. A., Hydrology, an introduction to hydrologic
science, Addisson Wesley, Inc., Reading, etc., 1990.
modifications:
06/2012 - code reformatting and cleanup
BN
10/2012 - ported aero_resist_used (Ra_used) from VIC 4.0.6 to store
the aerodynamic resistance actually used in flux
calculations.
BN
*******************************************************************************/
double
SnowPackEnergyBalance(double TSurf,
va_list ap)
{
extern option_struct options;
const char *Routine = "SnowPackEnergyBalance";
/* start of list of arguments in variable argument list */
double Dt; /* Model time step (hours) */
double Ra; /* Aerodynamic resistance (s/m) */
double *Ra_used; /* Aerodynamic resistance (s/m) after stability
correction */
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;
/* end of list of arguments in variable argument list */
double Density; /* Density of water/ice at TMean (kg/m3) */
double EsSnow; /* saturated vapor pressure in the snow pack
(Pa) */
double *LatentHeat; /* Latent heat exchange at surface (W/m2) */
double LongRadOut; /* long wave radiation emitted by surface
(W/m2) */
double Ls; /* Latent heat of sublimation (J/kg) */
double NetRad; /* Net radiation exchange at surface (W/m2) */
double RestTerm; /* Rest term in surface energy balance
(W/m2) */
double *SensibleHeat; /* Sensible heat exchange at surface (W/m2) */
double TMean; /* Average temperature for time step (C) */
/* Assign the elements of the array to the appropriate variables. The list
is traversed as if the elements are doubles, because:
In the variable-length part of variable-length argument lists, the old
``default argument promotions'' apply: arguments of type double are
always promoted (widened) to type double, and types char and short int
are promoted to int. Therefore, it is never correct to invoke
va_arg(argp, double); instead you should always use va_arg(argp,
double).
(quoted from the comp.lang.c FAQ list)
*/
Dt = (double) va_arg(ap, double);
Ra = (double) va_arg(ap, double);
Ra_used = (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 *);
/* Calculate active temp for energy balance as average of old and new */
/* TMean = 0.5 * (OldTSurf + TSurf); */
TMean = TSurf;
Density = RHO_W;
/* Correct aerodynamic conductance for stable conditions
Note: If air temp >> snow temp then aero_cond -> 0 (i.e. very stable)
velocity (vel_2m) is expected to be in m/sec */
/* Apply the stability correction to the aerodynamic resistance
NOTE: In the old code 2m was passed instead of Z-Displacement. I (bart)
think that it is more correct to calculate ALL fluxes at the same
reference level */
if (Wind > 0.0) {
*Ra_used = Ra / StabilityCorrection(Z, 0.f, TMean, Tair, Wind, Z0);
/* RASMnote - added to avoid decoupling under stable conditions
- BN - 01/2013 */
if (Tair > TMean && *Ra_used > RA_STABLE_UPPER_LIMIT) {
*Ra_used = RA_STABLE_UPPER_LIMIT;
}
}
else {
*Ra_used = HUGE_RESIST;
}
/* Calculate longwave exchange and net radiation */
LongRadOut = STEFAN * (TMean + 273.15) * (TMean + 273.15) *
(TMean + 273.15) * (TMean + 273.15);
NetRad = SurfAttenuation * ShortRad + LongRadIn - LongRadOut;
/* Calculate the sensible heat flux */
*SensibleHeat = AirDens * CP * (Tair - TMean) / *Ra_used;
/* Calculate the mass flux of ice to or from the surface layer */
/* Calculate the saturated vapor pressure in the snow pack */
EsSnow = svp(TMean) * 1000.;
*VaporMassFlux = AirDens * (EPS / Press) * (EactAir - EsSnow) / *Ra_used;
*VaporMassFlux /= Density;
if (Vpd == 0.0 && *VaporMassFlux < 0.0) {
*VaporMassFlux = 0.0;
}
/* Calculate latent heat flux */
if (TMean >= 0.0) {
/* Melt conditions: use latent heat of vaporization */
*LatentHeat = Lv * *VaporMassFlux * Density;
}
else {
/* Accumulation: use latent heat of sublimation (Eq. 3.19, Bras 1990 */
Ls = (677. - 0.07 * TMean) * JOULESPCAL * GRAMSPKG;
*LatentHeat = Ls * *VaporMassFlux * Density;
}
/* Calculate advected heat flux from rain
WORK IN PROGRESS: Should the following read (Tair - Tsurf) ?? */
*AdvectedEnergy = (CH_WATER * Tair * Rain) / (Dt * SECPHOUR);
/* Calculate change in cold content */
*DeltaColdContent = CH_ICE * SweSurfaceLayer * (TSurf - OldTSurf) /
(Dt * SECPHOUR);
/* Calculate Ground Heat Flux */
if (SnowDepth > 0.) {
*GroundFlux = 2.9302e-6 * SnowDensity * SnowDensity *
(TGrnd - TMean) / SnowDepth;
}
else {
*GroundFlux = 0;
}
*DeltaColdContent -= *GroundFlux;
/* Calculate net energy exchange at the snow surface */
RestTerm = NetRad + *SensibleHeat + *LatentHeat + *AdvectedEnergy -
*DeltaColdContent;
*RefreezeEnergy = (SurfaceLiquidWater * Lf * Density) / (Dt * SECPHOUR);
if (TSurf == 0.0 && RestTerm > -(*RefreezeEnergy)) {
*RefreezeEnergy = -RestTerm; /* available energy input over cold content
used to melt, i.e. Qrf is negative value
(energy out of pack)*/
RestTerm = 0.0;
}
else {
RestTerm += *RefreezeEnergy; /* add this positive value to the pack */
}
return RestTerm;
}
#undef GRAMSPKG
#undef CH_WATER
#undef JOULESPCAL
#undef EPS
#undef CP
Raw
vicNl_def.h
/*******************************************************************************
modifications:
12/2011 - updated this comment
- code reformatting and cleanup
- removed LINK_DEBUG code
BN
06/2012 - removed SAVE_STATE code
BN
08/2012 - removed ouput_format type
JCR
*******************************************************************************/
#ifndef VICN1_DEF_H
#define VICN1_DEF_H
/***** Model Constants *****/
#define MAXSTRING 512
#define MINSTRING 20
#define HUGE_RESIST 1.e20 /* largest allowable double number */
#define SMALL 1.e-12 /* smallest allowable double number */
#define MISSING -99999. /* missing value for multipliers in
BINARY format */
#define LITTLE 1 /* little-endian flag */
#define BIG 2 /* big-endian flag */
#define BROOKS 1 /* Brooks-Corey parameters for unsaturated flow */
#define ASCII 1 /* met file format flag */
#define BINARY 2 /* met file format flag */
/***** Forcing Variable Types *****/
#define N_FORCING_TYPES 14
#define AIR_TEMP 0 /* air temperature per time step (C) */
#define ALBEDO 1 /* surface albedo (fraction) */
#define DENSITY 2 /* atmospheric density (kg/m^3) */
#define LONGWAVE 3 /* incoming longwave radiation (W/m^2) */
#define PREC 4 /* precipitation (mm) */
#define PRESSURE 5 /* atmospheric pressure (kPa) */
#define SHORTWAVE 6 /* incoming shortwave (W/m^2) */
#define TMAX 7 /* maximum daily temperature (C) */
#define TMIN 8 /* minimum daily temperature (C) */
#define TSKC 9 /* cloud cover (fraction) */
#define VP 10 /* vapor pressure (kPa) */
#define WIND 11 /* wind speed (m/s) */
#define RA 12 /* aerodynaimc resistance (s/m)*/
#define SKIP 13 /* place holder for unused data columns */
/***** Physical Constants *****/
#define BARE_SOIL_ALBEDO 0.2 /* albedo for bare soil */
#define RESID_MOIST 0.0 /* define residual moisture content
of soil column */
#define ice_density 917. /* density of ice (kg/m^3) */
#define T_lapse 6.5 /* tempreature lapse rate of US Std
Atmos in C/km */
#define von_K 0.40 /* Von Karman constant for evapotranspiration */
#define KELVIN 273.15 /* conversion factor C to K */
#define STEFAN_B 5.6696e-8 /* stefan-boltzmann const in unit W/m^2/K^4 */
#define Lf 3.337e5 /* Latent heat of freezing (J/kg) at 0C */
#define RHO_W 1000.0 /* Density of water (kg/m^3) at 0C */
#define Cp 1004.0 /* Specific heat at constant pressure of air
(J/deg/K) */
#define CH_ICE 2100.0e3 /* Volumetric heat capacity (J/(m3*C)) of ice */
#define SECPHOUR 3600 /* seconds per hour */
#define SNOW_DT 5.0 /* Used to bracket snow surface temperatures
while computing the snow surface energy
balance (C) */
#define SURF_DT 1.0 /* Used to bracket soil surface temperatures
while computing energy balance (C) */
#define SOIL_DT 0.25 /* Used to bracket soil temperatures while
solving the soil thermal flux (C) */
#define HOURSPERDAY 24 /* number of hours per day */
#define HOURSPERYEAR 24 * 365 /* number of hours per year */
/***** Physical Constraints *****/
#define MINSOILDEPTH 0.001 /* minimum layer depth with which model can
work (m) */
#define STORM_THRES 0.001 /* threshold at which a new storm is
declared */
#define RA_STABLE_UPPER_LIMIT 100. /* upper limit for the aerodynamic resistance
/***** Define Boolean Values *****/
#ifndef FALSE
#define FALSE 0
#define TRUE !FALSE
#endif
#ifndef WET
#define WET 0
#define DRY 1
#endif
#ifndef SNOW
#define RAIN 0
#define SNOW 1
#endif
#ifndef MSC_PLUSPLUS
#define min(a, b) (a < b) ? a : b
#define max(a, b) (a > b) ? a : b
#endif
#include <user_def.h>
#include <snow.h>
#define DAYS_PER_YEAR 365.
#define DtoR 0.017453293 /* degrees to radians */
#ifndef PI
#define PI 3.1415927
#endif
#define STEFAN 5.6696e-8 /* Stefan boltzmann constant */
#define SOLAR_CONSTANT 1400.0 /* Solar constant in W/m^2 */
#define SEC_PER_DAY 86400. /* seconds per day */
/* define constants for saturated vapor pressure curve (kPa) */
#define A_SVP 0.61078
#define B_SVP 17.269
#define C_SVP 237.3
/* define constants for penman evaporation */
#define CP_PM 1013 /* specific heat of moist air J/kg/C
(Handbook of Hydrology) */
#define PS_PM 101300 /* sea level air pressure in Pa */
#define LAPSE_PM -0.006 /* environmental lapse rate in C/m */
/* global variables */
extern int NR; /* array index for atmos struct that indicates
the model step avarage or sum */
extern int NF; /* array index loop counter limit for atmos
struct that indicates the SNOW_STEP values */
/* number of restart variables in netcdf restart files */
#define RESTART_DIMS 5
/***** Data Structures *****/
/* netcdf restart data structure */
typedef struct {
int NDIMS;
int NNODE;
int NBANDS;
int NLAYER;
int NVEG_SIZE;
int NVEG;
int NDIST;
int IXDIM;
int *vegetat_type_num;
char *STILL_STORM;
float *DRY_TIME;
double *dz_node;
double *mu;
double *moist;
double *ice;
double *Wdew;
int *last_snow;
char *MELTING;
double *coverage;
double *swq;
double *surf_temp;
double *surf_water;
double *pack_temp;
double *pack_water;
double *density;
double *coldcontent;
double *snow_canopy;
double *T;
} netcdf_restart_struct;
/** file structures **/
typedef struct {
FILE *forcing[2]; /* atmospheric forcing data files */
FILE *globalparam; /* global parameters file */
FILE *init_snow; /* snowpack initialization file */
FILE *init_soil; /* soil temp and mosit initialization file */
FILE *snowband; /* snow elevation band data file */
FILE *soilparam; /* soil parameters for all grid cells */
FILE *veglib; /* vegetation parameters for all vege types */
FILE *vegparam; /* fractional coverage info for grid cell */
FILE *statefile; /* initial model state file */
} infiles_struct;
typedef struct {
FILE *fdepth;
FILE *fluxes;
FILE *snow;
FILE *snowband;
} outfiles_struct;
typedef struct {
char fdepth[MAXSTRING]; /* frozen soils depth (output) */
char fluxes[MAXSTRING]; /* grid cell surface fluxes (output) */
char forcing[2][MAXSTRING]; /* atmospheric forcing data file names */
char global[MAXSTRING]; /* global control file name */
char init_state[MAXSTRING]; /* initial model state file name */
char snow[MAXSTRING]; /* snow pack depth and swq (output) */
char snow_band[MAXSTRING]; /* snow band parameter file name */
char snowband[MAXSTRING]; /* snow band pack depth and swq (output) */
char soil[MAXSTRING]; /* soil parameter file name, or name of
file that has a list of all aoil
ARC/INFO files */
char soil_dir[MAXSTRING]; /* directory from which to read ARC/INFO
soil files */
char veg[MAXSTRING]; /* vegetation grid coverage file */
char veglib[MAXSTRING]; /* vegetation parameter library file */
} filenames_struct;
typedef struct {
char ARC_SOIL; /* TRUE = use ARC/INFO gridded ASCII files for soil
parameters*/
char COMPRESS; /* TRUE = Compress all output files */
char CORRPREC; /* TRUE = correct precipitation for gage undercatch */
char DIST_PRCP; /* TRUE = Use distributed precipitation model */
char FROZEN_SOIL; /* TRUE = Use frozen soils code */
char FULL_ENERGY; /* TRUE = Use full energy code */
char GLOBAL_LAI; /* TRUE = read LAI values for each vegetation type
from the veg param file */
char GRND_FLUX; /* TRUE = compute ground heat flux and energy
balance */
char INIT_STATE; /* TRUE = initialize model state from file */
char MOISTFRACT; /* TRUE = output soil moisture as moisture content */
char NOFLUX; /* TRUE = Use no flux lower bondary when computing
soil thermal fluxes */
char PRT_SNOW_BAND; /* TRUE = print snow parameters for each snow band */
char QUICK_FLUX; /* TRUE = Use Liang et al., 1999 formulation for
ground heat flux, if FALSE use explicit finite
difference method */
float MIN_WIND_SPEED; /* Minimum wind speed in m/s that can be used by
the model. **/
float PREC_EXPT; /* Exponential that controls the fraction of a
grid cell that receives rain during a storm
of given intensity */
int GRID_DECIMAL; /* Number of decimal places in grid file extensions */
int Nlayer; /* Number of layers in model */
int Nnode; /* Number of soil thermal nodes in the model */
int ROOT_ZONES; /* Number of root zones used in simulation */
int SNOW_BAND; /* Number of elevation bands over which to solve the
snow model */
float SNOW_STEP; /* Time step in hours to use when solving the
snow model */
} option_struct;
/*******************************************************************************
Stores forcing file input information.
*******************************************************************************/
typedef struct {
char SIGNED;
int SUPPLIED;
double multiplier;
} force_type_struct;
/*******************************************************************************
This structure records the parameters set by the forcing file input routines.
Those filled, are used to estimate the paramters needed for the model run in
initialize_atmos.c.
*******************************************************************************/
typedef struct {
force_type_struct TYPE[N_FORCING_TYPES];
int FORCE_DT[2]; /* forcing file time step */
int FORCE_ENDIAN[2]; /* endian-ness of input file, used for
DAILY_BINARY format */
int FORCE_FORMAT[2]; /* ASCII or BINARY */
int FORCE_INDEX[2][N_FORCING_TYPES];
int N_TYPES[2];
} param_set_struct;
/*******************************************************
This structure stores all model run global parameters.
*******************************************************/
typedef struct {
double MAX_SNOW_TEMP; /* maximum temperature at which snow can fall (C) */
double MIN_RAIN_TEMP; /* minimum temperature at which rain can fall (C) */
double measure_h; /* height of measurements (m) */
double wind_h; /* height of wind measurements (m) */
float resolution; /* Model resolution (degrees) */
float dt; /* Time step in hours (24/dt must be an integer) */
int endday; /* Last day of model simulation */
int endmonth; /* Last month of model simulation */
int endyear; /* Last year of model simulation */
int forceday[2]; /* day forcing files starts */
int forcehour[2]; /* hour forcing files starts */
int forcemonth[2]; /* month forcing files starts */
int forceskip[2]; /* number of model time steps to skip at the start of
the forcing file */
int forceyear[2]; /* year forcing files start */
int nrecs; /* Number of time steps simulated */
int skipyear; /* Number of years to skip before writing output data */
int startday; /* Starting day of the simulation */
int starthour; /* Starting hour of the simulation */
int startmonth; /* Starting month of the simulation */
int startyear; /* Starting year of the simulation */
} global_param_struct;
/***********************************************************
This structure stores the soil parameters for a grid cell.
***********************************************************/
typedef struct {
int FS_ACTIVE; /* if TRUE frozen soil algorithm is
active in current grid cell */
double Ds; /* fraction of maximum subsurface flow
rate */
double Dsmax; /* maximum subsurface flow rate
(mm/day) */
double Ksat[MAX_LAYERS]; /* saturated hydraulic conductivity
(mm/day) */
double Wcr[MAX_LAYERS]; /* critical moisture level for soil
layer, evaporation is no longer
affected moisture stress in the
soil (mm) */
double Wpwp[MAX_LAYERS]; /* soil moisture content at permanent
wilting point (mm) */
double Ws; /* fraction of maximum soil moisture */
double alpha[MAX_NODES]; /* thermal solution constant */
double annual_prec; /* annual average precipitation (mm) */
double avg_temp; /* average soil temperature (C) */
double b_infilt; /* infiltration parameter */
double beta[MAX_NODES]; /* thermal solution constant */
double bubble[MAX_LAYERS]; /* bubbling pressure, HBH 5.15 (cm) */
double bubble_node[MAX_NODES]; /* bubbling pressure (cm) */
double bulk_density[MAX_LAYERS]; /* soil bulk density (kg/m^3) */
double c; /* exponent */
double depth[MAX_LAYERS]; /* thickness of each soil moisture
layer (m) */
double dp; /* soil thermal damping depth (m) */
double dz_node[MAX_NODES]; /* thermal node thickness (m) */
double expt[MAX_LAYERS]; /* pore-size distribution per layer,
HBH 5.15 */
double expt_node[MAX_NODES]; /* pore-size distribution per node */
double gamma[MAX_NODES]; /* thermal solution constant */
double init_moist[MAX_LAYERS]; /* initial layer moisture level (mm) */
double max_infil; /* maximum infiltration rate */
double max_moist[MAX_LAYERS]; /* maximum moisture content (mm) per
layer */
double max_moist_node[MAX_NODES]; /* maximum moisture content (mm/mm) per
node */
double phi_s[MAX_LAYERS]; /* soil moisture diffusion parameter
(mm/mm) */
double porosity[MAX_LAYERS]; /* porosity (fraction) */
double quartz[MAX_LAYERS]; /* quartz content of soil (fraction) */
double resid_moist[MAX_LAYERS]; /* residual moisture content of soil
layer */
double rough; /* soil surface roughness (m) */
double snow_rough; /* snow surface roughness (m) */
double soil_density[MAX_LAYERS]; /* soil partical density (kg/m^3) */
double *AreaFract; /* Fraction of grid cell included in
each elevation band */
double *Pfactor; /* Change in Precipitation due to
elevation (fract) */
double *Tfactor; /* Change in temperature due to
elevation (C) */
char *AboveTreeLine; /* Flag to indicate if band is above
the treeline*/
#if QUICK_FS
double **ufwc_table_layer[MAX_LAYERS];
double **ufwc_table_node[MAX_NODES];
#endif
float elevation; /* grid cell elevation (m) */
float lat; /* grid cell central latitude */
float lng; /* grid cell central longitude */
float time_zone_lng; /* central meridian of the time zone */
float **layer_node_fract; /* fraction of all nodes within each
layer */
int gridcel; /* grid cell number */
} soil_con_struct;
/*******************************************************************************
This structure stores information about the vegetation coverage of the
current grid cell.
*******************************************************************************/
typedef struct {
double Cv; /* fraction of vegetation coverage */
double Cv_sum; /* total fraction of vegetation coverage */
float root[MAX_LAYERS]; /* percent of roots in each soil layer
(fraction) */
float *zone_depth; /* depth of root zone */
float *zone_fract; /* fraction of roots within root zone */
int veg_class; /* vegetation class reference number */
int vegetat_type_num; /* number of vegetation types in the grid cell */
} veg_con_struct;
/*******************************************************************************
This structure stores parameters for individual vegetation types.
*******************************************************************************/
typedef struct {
char overstory; /* TRUE = overstory present, important for snow
accumulation in canopy */
double LAI[12]; /* monthly leaf area index */
double Wdmax[12]; /* maximum monthly dew holding capacity (mm) */
double albedo[12]; /* vegetation albedo (added for full energy)
(fraction) */
double displacement[12]; /* vegetation displacement height (m) */
double emissivity[12]; /* vegetation emissivity (fraction) */
double rad_atten; /* radiation attenuation due to canopy,
default = 0.5 (N/A) */
double rarc; /* architectural resistance (s/m) */
double rmin; /* minimum stomatal resistance (s/m) */
double roughness[12]; /* vegetation roughness length (m) */
double trunk_ratio; /* ratio of trunk height to tree height,
default = 0.2 (fraction) */
double wind_atten; /* wind attenuation through canopy,
default = 0.5 (N/A) */
double wind_h; /* height at which wind is measured (m) */
float RGL; /* Value of solar radiation below which there
will be no transpiration (ranges from
~30 W/m^2 for trees to ~100 W/m^2 for crops) */
int veg_class; /* vegetation class reference number */
double baseLAI[12]; /* monthly leaf area index, added in 04072004,
Mark's e-mail-03-23-04 */
} veg_lib_struct;
/*******************************************************************************
This structure stores the atmospheric forcing data for each model time step
for a single grid cell. Each array stores the values for the SNOW_STEPs
during the current model step and the value for the entire model step. The
latter is referred to by array[NR]. Looping over the SNOW_STEPs is done by
for (i = 0; i < NF; i++)
*******************************************************************************/
typedef struct {
char snowflag[25]; /* TRUE if there is snowfall in any of the snow
bands during the timestep, FALSE otherwise*/
double air_temp[25]; /* air temperature (C) */
double density[25]; /* atmospheric density (kg/m^3) */
double longwave[25]; /* incoming longwave radiation (W/m^2) (net incoming
longwave for water balance model) */
double out_prec; /* Total precipitation for time step - accounts
for corrected precipitation totals */
double prec[25]; /* average precipitation in grid cell (mm) */
double pressure[25]; /* atmospheric pressure (kPa) */
double shortwave[25]; /* incoming shortwave radiation (W/m^2) */
double vp[25]; /* atmospheric vapor pressure (kPa) */
double vpd[25]; /* atmospheric vapor pressure deficit (kPa) */
double wind[25]; /* wind speed (m/s) */
double ra[25]; /* aerodynamic resistance (s/m)*/
} atmos_data_struct;
/*******************************************************************************
This structure stores information about the time and date of the current time
step.
*******************************************************************************/
typedef struct {
int day; /* current day */
int day_in_year; /* julian day in year */
int hour; /* beginning of current hour */
int month; /* current month */
int year; /* current year */
} dmy_struct; /* array of length nrec created */
/*******************************************************************************
This structure stores all soil variables for each layer in the soil column.
*******************************************************************************/
typedef struct {
double Cs; /* average volumetric heat capacity of the
current layer (J/m^3/K) */
double T; /* temperature of the unfrozen sublayer (C) */
double evap; /* evapotranspiration from soil layer (mm) */
double ice; /* ice content of the frozen sublayer (mm) */
double kappa; /* average thermal conductivity of the current
layer (W/m/K) */
double moist; /* moisture content of the unfrozen sublayer
(mm) */
double phi; /* moisture diffusion parameter */
} layer_data_struct;
/*******************************************************************************
This structure stores soil variables for the complete soil column for each
grid cell.
*******************************************************************************/
typedef struct {
double aero_resist[3]; /* aerodynamic resistane (s/m)
[0] = over bare vegetation or soil
[2] = over snow */
double aero_resist_short_ref_crop[3]; /* same for short ref crop*/
double aero_resist_tall_ref_crop[3]; /* same for short ref crop*/
double aero_resist_used; /* The (stability-corrected)
aerodynamic resistance (s/m) that
was actually used in flux
calculations. For cases in which a
cell uses 2 different resistances
for flux computations in the same
time step (i.e. cell contains
overstory and snow is present on
the ground), aero_resist_used will
contain the snow pack's
resistance. */
double baseflow; /* baseflow from current cell
(mm/TS) */
double inflow; /* moisture that reaches the top of
the soil column (mm) */
double runoff; /* runoff from current cell (mm/TS) */
layer_data_struct layer[MAX_LAYERS]; /* structure containing soil variables
for each layer (see above) */
} cell_data_struct;
/*******************************************************************************
This structure stores energy balance components, and variables used to solve
the thermal fluxes through the soil column.
*******************************************************************************/
typedef struct {
char frozen; /* TRUE = frozen soil present */
double Cs[2]; /* heat capacity for top two layers
(J/m^3/K) */
double Cs_node[MAX_NODES]; /* heat capacity of the soil thermal nodes
(J/m^3/K) */
double T[MAX_NODES]; /* thermal node temperatures (C) */
double Trad[2]; /* surface temperature of energy balance
(C) */
double advection; /* advective flux (Wm-2) */
double albedo; /* surface albedo (fraction) */
double deltaCC; /* change in snow heat storage (Wm-2) */
double deltaH; /* change in soil heat storage (Wm-2) */
double error; /* energy balance error (W/m^2) */
double fdepth[MAX_FRONTS]; /* all simulated freezing front depths */
double grnd_flux; /* ground heat flux (Wm-2) */
double ice[MAX_NODES]; /* thermal node ice content */
double kappa[2]; /* soil thermal conductivity for top two
layers (W/m/K) */
double kappa_node[MAX_NODES]; /* thermal conductivity of the soil thermal
nodes (W/m/K) */
double latent; /* net latent heat flux (Wm-2) */
double longwave; /* net longwave flux (Wm-2) */
double moist[MAX_NODES]; /* thermal node moisture content */
double refreeze_energy; /* energy used to refreeze the snowpack
(Wm-2) */
double sensible; /* net sensible heat flux (Wm-2) */
double shortwave; /* incoming shortwave heat (Wm-2) */
double snow_flux; /* thermal flux through the snow pack
(Wm-2) */
double tdepth[MAX_FRONTS]; /* all simulated thawing front depths */
double unfrozen; /* frozen layer water content that is
unfrozen */
int Nfrost; /* number of simulated freezing fronts */
int Nthaw; /* number of simulated thawing fronts */
int T1_index; /* soil node at the bottom of the top layer */
double surf_atten; /* surface attenuation factor */
} energy_bal_struct;
/*******************************************************************************
This structure stores vegetation variables for each vegetation type in a grid
cell.
*******************************************************************************/
typedef struct {
double canopyevap; /* evaporation from canopy (mm/TS) */
double throughfall; /* water that reaches the ground through the
canopy (mm/TS) */
double Wdew; /* dew trapped on vegetation (mm) */
} veg_var_struct;
/*******************************************************************************
This structure stores snow pack variables needed to run the snow model.
*******************************************************************************/
typedef struct {
char MELTING; /* flag indicating that snowpack melted
previously */
int snow; /* TRUE = snow, FALSE = no snow */
double Qnet; /* New energy at snowpack surface */
double albedo; /* snow surface albedo (fraction) */
double canopy_vapor_flux; /* depth of water evaporation, sublimation, or
condensation from intercepted snow (m) */
double coldcontent; /* cold content of snow pack */
double coverage; /* fraction of snow band that is covered with
snow */
double density; /* snow density (kg/m^3) */
double depth; /* snow depth (m) */
double mass_error; /* snow mass balance error */
double pack_temp; /* depth averaged temperature of the snowpack
(C) */
double pack_water; /* liquid water content of the snow pack (m) */
double snow_canopy; /* amount of snow on canopy (m) */
double surf_temp; /* depth averaged temperature of the snow pack
surface layer (C) */
double surf_water; /* liquid water content of the surface layer
(m) */
double swq; /* snow water equivalent of the entire pack (m) */
double tmp_int_storage; /* temporary canopy storage, used in
snow_canopy */
double vapor_flux; /* depth of water evaporation, sublimation, or
condensation from snow pack (m) */
int last_snow; /* time steps since last snowfall */
} snow_data_struct; /* an array of size Nrec */
/*******************************************************************************
This structure stores all variables needed to solve, or save solututions for
all versions of this model. Vegetation and soil variables are created for
both wet and dry fractions of the grid cell (for use with the distributed
precipitation model).
*******************************************************************************/
typedef struct {
cell_data_struct **cell[2]; /* Stores soil layer variables (wet and
dry) */
double *mu; /* fraction of grid cell that receives
precipitation */
energy_bal_struct **energy; /* Stores energy balance variables */
snow_data_struct **snow; /* Stores snow variables */
veg_var_struct **veg_var[2]; /* Stores vegetation variables (wet and
dry) */
double frac_area; /* cell area as fraction of total area */
soil_con_struct *soil; /* soil parameters carried by dist_prcp_st*/
veg_con_struct *veg; /* veg parameters carried by dist_prcp_st*/
char STILL_STORM;
float DRY_TIME;
} dist_prcp_struct;
/*******************************************************************************
This structure stores all variables needed for output.
*******************************************************************************/
typedef struct {
int id; /* gridcell id */
double lat;
double lon;
double Wdew; /* canopy interception of moisture */
double advection[MAX_BANDS + 1]; /* grid cell advection (snow only)
(Wm-2) */
double aero_cond; /* grid cell mean aerodynamic
conductance [m/s] */
double aero_resist; /* grid cell mean aerodynamic
resistence [s/m] */
double air_temp; /* grid cell air temperature */
double albedo; /* grid cell mean albedo */
double baseflow; /* baseflow out of the bottom layer */
double bot_energy_error[2];
double coverage[MAX_BANDS + 1]; /* fractional coverage of grid cell
with snow */
double deltaCC[MAX_BANDS + 1]; /* change of cold content in the
snowpack [Wm-2] */
double deltaH; /* grid cell change in heat storage
(snow only) */
double energy_error; /* energy balance error */
double evap; /* grid cell evaporation */
double evap_bare; /* grid cell net evaporation from
bare ground */
double evap_canop; /* grid cell net evaporation from
canopy interception */
double evap_veg; /* grid cell net evapotraspiration
from vegetation */
double fdepth[MAX_FRONTS]; /* depth of all freezing fronts */
double grnd_flux; /* grid cell ground flux */
double ice[MAX_LAYERS]; /* frozen layer ice content */
double in_long; /* grid cell net incoming longwave
flux */
double inflow; /* moisture that reaches the top of
the soil column */
double latent; /* grid cell net latent heat flux */
double moist[MAX_LAYERS]; /* current moisture in each layer */
double net_long; /* grid cell net longwave flux */
double net_short; /* grid cell net shortwave flux */
double prec; /* incoming precipitation */
double r_net; /* grid cell net radiation W/m^2 */
double rad_temp; /* grid cell average radiative surface
temperature */
double refreeze_energy[MAX_BANDS + 1]; /* energy used to refreeze snowpack
[Wm-2] */
double rel_humid;
double runoff; /* runoff from the surface */
double sensible; /* grid cell net sensible heat flux */
double shortwave; /* grid cell incoming shortwave flux */
double snow_canopy[MAX_BANDS + 1]; /* snow captured by canopy (mm) */
double snow_depth[MAX_BANDS + 1]; /* snow depth (cm) */
double snow_flux[MAX_BANDS + 1]; /* energy flux through the snowpack
[Wm-2] */
double sub_canop; /* grid cell net sublimation from
canopy interception */
double sub_snow; /* grid cell net sublimation from
bare ground from snow pack */
double surf_cond; /* grid cell mean surface conductance
[m/s] */
double surf_temp; /* grid cell average daily surface
temperature */
double swq[MAX_BANDS + 1]; /* snow water equivalent (mm) */
double tdepth[MAX_FRONTS]; /* depth of all thawing fronts */
double wind; /* grid cell wind speed */
double swband[MAX_BANDS + 1]; /* store shortwave by snow band*/
double lwband[MAX_BANDS + 1]; /* store longwave by snow band*/
double albedoband[MAX_BANDS + 1]; /* store snow albedo by snow band*/
double latentband[MAX_BANDS + 1]; /* store latent heat by snow band*/
double sensibleband[MAX_BANDS + 1]; /* store sensible heat by snow band*/
double grndband[MAX_BANDS + 1]; /* store ground heat*/
/* added for rcm interface */
double surf_humid; /* grid cell average surface humidity
kg/kg */
/* RASMnote: Confusing - comment does not appear to match what is done - BN
- 06/2012 */
double Tlayer[3]; /* changed from MAX_LAYERS to 5 by Zhu
08/05/05 */
} out_data_struct;
/*******************************************************************************
This structure holds all variables needed for the error handling routines.
*******************************************************************************/
typedef struct {
atmos_data_struct *atmos;
double dt;
energy_bal_struct *energy;
infiles_struct infp;
int rec;
out_data_struct *out_data;
outfiles_struct outfp;
snow_data_struct *snow;
soil_con_struct soil_con;
veg_con_struct *veg_con;
veg_var_struct *veg_var;
} Error_struct;
/* used in rcm coupling */
typedef struct {
double storage;
double error_water;
double error_energy;
} Error_struct_2;
typedef struct {
double last_storage;
double cum_error_water;
double max_error_water;
double cum_error_energy;
double max_error_energy;
} Error_struct_3;
typedef struct {
double precip;
double evap;
double soil_evap;
double canopy_evap;
double intercept_evap;
double runoff;
double baseflow;
double swq;
double air_temp;
double latent;
double sensible;
double net_rad;
double start_store;
double end_store;
int count;
} basin_struct;
#endif
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment