Skip to content

Commit d9cab01

Browse files
HishamEldardiryHishamEldardiry
authored
Add files via upload
1 parent 58d1fe9 commit d9cab01

File tree

97 files changed

+19905
-0
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

97 files changed

+19905
-0
lines changed

‎Rainfall_Runoff/IceEnergyBalance.c‎

Lines changed: 274 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,274 @@
1+
/*
2+
* SUMMARY: IceEnergyBalance.c - Calculate lake ice energy balance
3+
* USAGE: Part of the lake algorithm
4+
*
5+
* AUTHOR: Laura Bowling
6+
* ORG: University of Washington, Department of Civil Engineering
7+
* E-MAIL: [email protected]
8+
* ORIG-DATE: March 16, 2001
9+
* LAST-MOD: Mon Apr 21 15:51:12 2003 by Keith Cherkauer <[email protected]>
10+
* DESCRIPTION: Calculate ice energy balance
11+
* DESCRIP-END.
12+
* FUNCTIONS: IceEnergyBalance()
13+
* COMMENTS:
14+
*/
15+
16+
#include <stdarg.h>
17+
#include <stdio.h>
18+
#include <stdlib.h>
19+
#include <vicNl.h>
20+
21+
static char vcid[] = "$Id$";
22+
23+
/*****************************************************************************
24+
Function name: IceEnergyBalance()
25+
26+
Purpose : Calculate the surface energy balance for the snow pack
27+
28+
Required :
29+
double TSurf - new estimate of effective surface temperature
30+
va_list ap - Argument list initialized by va_start(). For
31+
elements of list and order, see beginning of
32+
routine
33+
34+
Returns :
35+
double RestTerm - Rest term in the energy balance
36+
37+
Modifies :
38+
double *RefreezeEnergy - Refreeze energy (W/m2)
39+
double *VaporMassFlux - Mass flux of water vapor to or from the
40+
intercepted snow (m/s)
41+
42+
Comments :
43+
Reference: Bras, R. A., Hydrology, an introduction to hydrologic
44+
science, Addisson Wesley, Inc., Reading, etc., 1990.
45+
46+
Modifications:
47+
16-Jul-04 Renamed VaporMassFlux to vapor_flux to denote fact that its
48+
units are m/timestep rather than kg/m2s. Created new variable
49+
VaporMassFlux with units of kg/m2s. After VaporMassFlux is
50+
computed, vapor_flux is derived from it by unit conversion.
51+
vapor_flux is the variable that is passed in/out of this
52+
function. TJB
53+
24-Aug-04 Added logic to handle blowing_flux and vapor_flux. TJB
54+
28-Sep-04 Added Ra_used to store the aerodynamic resistance used in flux
55+
calculations. TJB
56+
04-Oct-04 Merged with Laura Bowling's updated lake model code. Now
57+
blowing snow sublimation is calculated for lakes. TJB
58+
2006-Sep-23 Replaced redundant STEFAN constant with STEFAN_B. TJB
59+
2006-Nov-07 Removed LAKE_MODEL option. TJB
60+
61+
*****************************************************************************/
62+
double IceEnergyBalance(double TSurf, va_list ap)
63+
{
64+
65+
extern option_struct options;
66+
67+
const char *Routine = "IceEnergyBalance";
68+
69+
/* start of list of arguments in variable argument list */
70+
71+
double Dt; /* Model time step (hours) */
72+
double Ra; /* Aerodynamic resistance (s/m) */
73+
double *Ra_used; /* Aerodynamic resistance (s/m) after stability correction */
74+
double Z; /* Reference height (m) */
75+
double Displacement; /* Displacement height (m) */
76+
double Z0; /* surface roughness height (m) */
77+
double Wind; /* Wind speed (m/s) */
78+
double ShortRad; /* Net incident shortwave radiation (W/m2) */
79+
double LongRadIn; /* Incoming longwave radiation (W/m2) */
80+
double AirDens; /* Density of air (kg/m3) */
81+
double Lv; /* Latent heat of vaporization (J/kg3) */
82+
double Tair; /* Air temperature (C) */
83+
double Press; /* Air pressure (Pa) */
84+
double Vpd; /* Vapor pressure deficit (Pa) */
85+
double EactAir; /* Actual vapor pressure of air (Pa) */
86+
double Rain; /* Rain fall (m/timestep) */
87+
double SweSurfaceLayer; /* Snow water equivalent in surface layer (m)
88+
*/
89+
double SurfaceLiquidWater; /* Liquid water in the surface layer (m) */
90+
double OldTSurf; /* Surface temperature during previous time
91+
step */
92+
double *RefreezeEnergy; /* Refreeze energy (W/m2) */
93+
double *vapor_flux; /* Total mass flux of water vapor to or from
94+
snow (m/timestep) */
95+
double *blowing_flux; /* Mass flux of water vapor to or from
96+
blowing snow (m/timestep) */
97+
double *surface_flux; /* Mass flux of water vapor to or from
98+
snow pack (m/timestep) */
99+
double *AdvectedEnergy; /* Energy advected by precipitation (W/m2) */
100+
double DeltaColdContent; /* Change in cold content (W/m2) */
101+
double Tfreeze;
102+
double AvgCond;
103+
double SWconducted;
104+
double SnowDepth;
105+
double SnowDensity;
106+
double SurfAttenuation;
107+
double *qf; /* Ground Heat Flux (W/m2) */
108+
double *LatentHeat; /* Latent heat exchange at surface (W/m2) */
109+
double *LatentHeatSub; /* Latent heat exchange at surface (W/m2) due to sublimation */
110+
double *SensibleHeat; /* Sensible heat exchange at surface (W/m2) */
111+
double *LongRadOut;
112+
113+
/* end of list of arguments in variable argument list */
114+
115+
double Density; /* Density of water/ice at TMean (kg/m3) */
116+
double EsSnow; /* saturated vapor pressure in the snow pack
117+
(Pa) */
118+
119+
double Ls; /* Latent heat of sublimation (J/kg) */
120+
double NetRad; /* Net radiation exchange at surface (W/m2) */
121+
double RestTerm; /* Rest term in surface energy balance
122+
(W/m2) */
123+
double TMean; /* Average temperature for time step (C) */
124+
double qnull;
125+
double VaporMassFlux; /* Total mass flux of water vapor to or from
126+
snow (kg/m2s) */
127+
double BlowingMassFlux; /* Mass flux of water vapor to or from
128+
blowing snow (kg/m2s) */
129+
double SurfaceMassFlux; /* Mass flux of water vapor to or from
130+
snow pack (kg/m2s) */
131+
132+
/* Assign the elements of the array to the appropriate variables. The list
133+
is traversed as if the elements are doubles, because:
134+
135+
In the variable-length part of variable-length argument lists, the old
136+
``default argument promotions'' apply: arguments of type double are
137+
always promoted (widened) to type double, and types char and short int
138+
are promoted to int. Therefore, it is never correct to invoke
139+
va_arg(argp, double); instead you should always use va_arg(argp,
140+
double).
141+
142+
(quoted from the comp.lang.c FAQ list)
143+
*/
144+
Dt = (double) va_arg(ap, double);
145+
Ra = (double) va_arg(ap, double);
146+
Ra_used = (double *) va_arg(ap, double *);
147+
Z = (double) va_arg(ap, double);
148+
Displacement = (double) va_arg(ap, double);
149+
Z0 = (double) va_arg(ap, double);
150+
Wind = (double) va_arg(ap, double);
151+
ShortRad = (double) va_arg(ap, double);
152+
LongRadIn = (double) va_arg(ap, double);
153+
AirDens = (double) va_arg(ap, double);
154+
Lv = (double) va_arg(ap, double);
155+
Tair = (double) va_arg(ap, double);
156+
Press = (double) va_arg(ap, double);
157+
Vpd = (double) va_arg(ap, double);
158+
EactAir = (double) va_arg(ap, double);
159+
Rain = (double) va_arg(ap, double);
160+
SweSurfaceLayer = (double) va_arg(ap, double);
161+
SurfaceLiquidWater = (double) va_arg(ap, double);
162+
OldTSurf = (double) va_arg(ap, double);
163+
RefreezeEnergy = (double *) va_arg(ap, double *);
164+
vapor_flux = (double *) va_arg(ap, double *);
165+
blowing_flux = (double *) va_arg(ap, double *);
166+
surface_flux = (double *) va_arg(ap, double *);
167+
AdvectedEnergy = (double *) va_arg(ap, double *);
168+
DeltaColdContent = (double) va_arg(ap, double );
169+
Tfreeze = (double) va_arg(ap, double);
170+
AvgCond = (double) va_arg(ap, double);
171+
SWconducted = (double) va_arg(ap, double);
172+
SnowDepth = (double) va_arg(ap, double);
173+
SnowDensity = (double) va_arg(ap, double);
174+
SurfAttenuation = (double) va_arg(ap, double);
175+
qf = (double *) va_arg(ap, double *);
176+
LatentHeat = (double *) va_arg(ap, double *);
177+
LatentHeatSub = (double *) va_arg(ap, double *);
178+
SensibleHeat = (double *) va_arg(ap, double *);
179+
LongRadOut = (double *) va_arg(ap, double *);
180+
181+
/* Calculate active temp for energy balance as average of old and new */
182+
183+
/* TMean = 0.5 * (OldTSurf + TSurf); */
184+
TMean = TSurf;
185+
Density = RHO_W;
186+
187+
/* Correct aerodynamic conductance for stable conditions
188+
Note: If air temp >> snow temp then aero_cond -> 0 (i.e. very stable)
189+
velocity (vel_2m) is expected to be in m/sec */
190+
191+
/* Apply the stability correction to the aerodynamic resistance
192+
NOTE: In the old code 2m was passed instead of Z-Displacement. I (bart)
193+
think that it is more correct to calculate ALL fluxes at the same
194+
reference level */
195+
196+
197+
if (Wind > 0.0)
198+
*Ra_used = Ra / StabilityCorrection(Z, 0.f, TMean, Tair, Wind, Z0);
199+
/* *Ra_used = Ra / StabilityCorrection(2.f, 0.f, TMean, Tair, Wind, Z0);*/
200+
else
201+
*Ra_used = HUGE_RESIST;
202+
203+
/* Calculate longwave exchange and net radiation */
204+
205+
*LongRadOut = LongRadIn - STEFAN_B * (TMean+273.15) * (TMean+273.15)
206+
* (TMean+273.15) * (TMean+273.15);
207+
NetRad = ShortRad + *LongRadOut;
208+
209+
/* Calculate the sensible heat flux */
210+
211+
*SensibleHeat = AirDens * CP_PM * (Tair - TMean) / *Ra_used;
212+
213+
/* Calculate the mass flux of ice to or from the surface layer */
214+
215+
/* Calculate the saturated vapor pressure in the snow pack,
216+
(Equation 3.32, Bras 1990) */
217+
218+
EsSnow = svp(TMean) /* * 1000. */;
219+
220+
/* EsSnow = 610.78 * exp((double)((17.269 * TMean) / (237.3 + TMean))); */
221+
222+
/* if (TMean < 0.0) */
223+
/* EsSnow *= 1.0 + .00972 * TMean + .000042 * pow((double)TMean,(double)2.0); */
224+
225+
/* Calculate sublimation terms and latent heat flux */
226+
227+
/* blowing_flux was calculated outside of the root_brent iteration */
228+
BlowingMassFlux = *blowing_flux * Density / (Dt * SECPHOUR);
229+
230+
latent_heat_from_snow(AirDens, EactAir, Lv, Press, Ra, TMean, Vpd,
231+
LatentHeat, LatentHeatSub, &VaporMassFlux, &BlowingMassFlux,
232+
&SurfaceMassFlux);
233+
234+
/* Convert sublimation terms from kg/m2s to m/timestep */
235+
*vapor_flux = VaporMassFlux * Dt * SECPHOUR / Density;
236+
*surface_flux = SurfaceMassFlux * Dt * SECPHOUR / Density;
237+
238+
/* Calculate advected heat flux from rain
239+
WORK IN PROGRESS: Should the following read (Tair - Tsurf) ?? */
240+
241+
// Temporary fix for lake model.
242+
*AdvectedEnergy = (CH_WATER * Tair * Rain) / (Dt*SECPHOUR);
243+
//*AdvectedEnergy = 0.0;
244+
245+
/* Calculate change in cold content */
246+
/* No change in cold content in lake model */
247+
248+
/* Changes for lake model start here. Actually, equals qo-Io (P&H eq. 7)*/
249+
/* because Io (SWnet) is included in NetRad below. */
250+
qnull = (1/AvgCond)*(Tfreeze - TMean + SWconducted);
251+
*qf = qnull;
252+
253+
/* Changes for lake model end here. */
254+
255+
/* Calculate net energy exchange at the snow surface */
256+
257+
RestTerm = ( NetRad + *SensibleHeat + *LatentHeat + *LatentHeatSub + *AdvectedEnergy
258+
+ qnull );
259+
260+
*RefreezeEnergy = (SurfaceLiquidWater * Lf * Density)/(Dt * SECPHOUR);
261+
262+
/* Melting, or partially refreeze surface water. */
263+
if (TSurf == 0.0 && RestTerm > -(*RefreezeEnergy)) {
264+
*RefreezeEnergy = -RestTerm; /* available energy input over cold content
265+
used to melt, i.e. Qrf is negative value
266+
(energy out of pack)*/
267+
RestTerm = 0.0;
268+
}
269+
else { /* Pack is getting colder. */
270+
RestTerm += *RefreezeEnergy; /* add this positive value to the pack */
271+
}
272+
273+
return RestTerm;
274+
}

‎Rainfall_Runoff/IceEnergyBalance.o‎

15.8 KB
Binary file not shown.

‎Rainfall_Runoff/LAKE.h‎

Lines changed: 116 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,116 @@
1+
/******************************************************************************
2+
// $Id$
3+
Modifications:
4+
2006-Nov-07 Removed LAKE_MODEL option. TJB
5+
2007-Aug-16 Made return value of initialize_prcp an int. JCA
6+
2007-Aug-21 Added features for EXCESS_ICE option. JCA
7+
2007-Oct-24 Changed get_sarea, get_volume, and get_depth to return exit
8+
status so that errors can be trapped and communicated up the
9+
chain of function calls. KAC via TJB
10+
2007-Oct-24 Changed lakeice() to return exit status. KAC via TJB
11+
2007-Nov-06 New lake physics parameters. Modified argument lists for
12+
various functions. Moved get_dist() to vicNl.h. LCB via TJB
13+
2008-Apr-21 Added argument to alblake. LCB via TJB
14+
2008-Sep-10 Updated values of CONDS and lamwlw to match Laura
15+
Bowling's lake work. LCB via TJB
16+
2009-Jul-31 Removed lakemain() and wetland_energy(); initialize_lake
17+
no longer takes a snow structure as input. TJB
18+
2009-Sep-28 Removed initialize_prcp and update_prcp. Modified
19+
argument list of initialize_lake. TJB
20+
2009-Sep-30 Miscellaneous fixes for lake model. TJB
21+
2009-Oct-05 Added functions for updating/rescaling lake and wetland
22+
fluxes and storages when lake area changes. TJB
23+
2010-Nov-11 Added skip_hydro flag to initialize_lake() arg list.
24+
Removed rescale_lake_fluxes(). TJB
25+
2010-Nov-26 Changed the argument list of water_balance(). TJB
26+
2010-Dec-28 Added latitude to alblake() arglist. TJB
27+
2011-Mar-01 Added rescale_snow_storage(). Added terms to argument
28+
list of initialize_lake(). TJB
29+
2013-Jul-25 Added advect_carbon_storage(). TJB
30+
2013-Dec-26 Removed EXCESS_ICE option. TJB
31+
2014-Mar-28 Removed DIST_PRCP option. TJB
32+
******************************************************************************/
33+
34+
//#ifndef LAKE_SET
35+
36+
#define LAKE_SET
37+
#define TMELT 0.0
38+
#define EMICE 0.97 /* Ice emissivity */
39+
#define EMH2O .98
40+
#define RHOSNOW 250. /* densities of water and snow */
41+
#define RHOICE 917. /* ice density*/
42+
#define rhosurf 1.275 /* surface air density */
43+
#define MAX_SURFACE_LAKE .6 /* max. surface layer thickness for E-B (m) */
44+
#define BETA 0.001 /* Curve shape parameter for lake profile. */
45+
#define FRACMIN 0.10 /* min ice thickness in meters */
46+
#define FRACLIM 0.02 /* lower limit on fractional ice cover */
47+
#define CPW_ICE 4200. /* specific heat of ice */
48+
#define DM 1.38889E-07 /* molecular diffusivity of water */
49+
#define SNOWCRIT 0.05 /* for albedo, in m */
50+
//#define G 9.80616
51+
#define ZWATER 0.0045 // 0.004 - original value
52+
#define ZSNOW 0.005
53+
#define CONDI 2.3 /* thermal conductivity of ice */
54+
#define CONDS 0.7 /* thermal conductivity of snow */
55+
56+
// attenuation of short and longwave radiation through ice (1/m)
57+
#define lamisw 1.5 // 1.5 in Patterson & Hamblin
58+
#define lamilw 20 // 20.0 in Patterson & Hamblin
59+
// attenuation of short and longwave radiation through snow (1/m)
60+
#define lamssw 6.0 // 6.0 in Patterson & Hamblin
61+
#define lamslw 20 // 20.0 in Patterson & Hamblin
62+
// attenuation of short and longwave radiation through water (1/m)
63+
#define lamwsw .3 // San Fran Bay data: 0.31 - 29.9 1/m (visible)
64+
#define lamwlw 1.4 // Hostetler and Bartlein assume 0.85 1/m (total)
65+
#define a1 0.7 /* Percent of radiation in visible band. */
66+
#define a2 0.3 /* Percent of radiation in infrared band. */
67+
#define QWTAU 86400./2. /* D. Pollard sub-ice time constant. */
68+
#define RADIUS 6371.228 /* Earth radius in km. */
69+
70+
//#endif // LAKE_SET
71+
72+
/*** Subroutine prototypes ***/
73+
74+
double adjflux(double, double, double ,double, double, double, double,
75+
double, double, double, double *, double *);
76+
void advect_carbon_storage(double, double, lake_var_struct *, cell_data_struct *);
77+
void advect_soil_veg_storage(double, double, double, double *, soil_con_struct *, veg_con_struct *, cell_data_struct *, veg_var_struct *, lake_con_struct);
78+
void advect_snow_storage(double, double, double, snow_data_struct *);
79+
void alblake(double, double, double *, double *, float *, float *, double, double,
80+
int, int *, double, double, char *, int, double);
81+
double calc_density(double);
82+
double CalcIcePackEnergyBalance(double Tsurf, ...);
83+
void colavg (double *, double *, double *, float, double *, int, double, double);
84+
float dragcoeff(float, double, double);
85+
void eddy (int, double, double * , double *, double *, double, int, double, double);
86+
void energycalc(double *, double *, int, double, double,double *, double *, double *);
87+
double ErrorIcePackEnergyBalance(double Tsurf, ...);
88+
double ErrorPrintIcePackEnergyBalance(double, va_list);
89+
int get_depth(lake_con_struct, double, double *);
90+
int get_sarea(lake_con_struct, double, double *);
91+
int get_volume(lake_con_struct, double, double *);
92+
void iceform (double *,double *,double ,double,double *,int, int, double, double, double *, double *, double *, double *, double *, double);
93+
void icerad(double,double ,double,double *, double *,double *);
94+
int ice_melt(double, double, double *, double, snow_data_struct *, lake_var_struct *, int, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double *, double *, double *, double *, double *, double *, double *, double *, double *, double);
95+
double IceEnergyBalance(double, va_list);
96+
int initialize_lake(lake_var_struct *, lake_con_struct, soil_con_struct *, cell_data_struct *, double, int);
97+
int lakeice(double *, double, double, double, double, int,
98+
double, double, double *, double, double, int, dmy_struct, double *, double *, double, double);
99+
void latsens(double,double, double, double, double, double, double, double,
100+
double *, double *, double);
101+
float lkdrag(float, double, double, double, double);
102+
lake_con_struct read_lakeparam(FILE *, soil_con_struct, veg_con_struct *);
103+
void rescale_soil_veg_fluxes(double, double, cell_data_struct *, veg_var_struct *);
104+
void rescale_snow_energy_fluxes(double, double, snow_data_struct *, energy_bal_struct *);
105+
void rescale_snow_storage(double, double, snow_data_struct *);
106+
void rhoinit(double *, double);
107+
int solve_lake(double, double, double, double, double, double, double, double,
108+
double, double, lake_var_struct *, lake_con_struct,
109+
soil_con_struct, int, int, double, dmy_struct, double);
110+
double specheat (double);
111+
void temp_area(double, double, double, double *, double *, double *, double *, int, double *, int, double, double, double*, double *, double *);
112+
void tracer_mixer(double *, int *, int, double*, int, double, double, double *);
113+
void tridia(int, double *, double *, double *, double *, double *);
114+
int water_balance (lake_var_struct *, lake_con_struct, int, all_vars_struct *, int, int, int, double, soil_con_struct, veg_con_struct);
115+
int water_energy_balance(int, double*, double*, int, int, double, double, double, double, double, double, double, double, double, double, double, double, double, double *, double *, double *, double*, double *, double *, double *, double, double *, double *, double *, double *, double *, double);
116+
int water_under_ice(int, double, double, double *, double *, double, int, double, double, double, double *, double *, double *, double *, int, double, double, double, double *);

0 commit comments

Comments
 (0)