.gitignore

@@ -15,5 +15,13 @@
 *.dylib
 
 vicNl
+DEPEND
+.depend
+VIC.sublime-project
+VIC.sublime-workspace
+*.sublime-project
+*.sublime-workspace
 
 .DS_Store
+
+*.uncrustify

samples/global.param.clean.sample

@@ -0,0 +1,154 @@
+#######################################################################
+# VIC Model Parameters - 4.1.x
+#######################################################################
+
+#######################################################################
+# Simulation Parameters
+#######################################################################
+NLAYER              3
+NODES               10
+TIME_STEP           3
+SNOW_STEP           3
+STARTYEAR           2000
+STARTMONTH          01
+STARTDAY            01 
+STARTHOUR           00
+ENDYEAR             2000
+ENDMONTH            12
+ENDDAY              31
+
+#######################################################################
+# Energy Balance Parameters
+#######################################################################
+FULL_ENERGY         FALSE
+#CLOSE_ENERGY       FALSE
+
+#######################################################################
+# Soil Temperature Parameters
+#######################################################################
+FROZEN_SOIL         FALSE
+#QUICK_FLUX         FALSE
+#IMPLICIT           TRUE
+#QUICK_SOLVE        FALSE
+#NO_FLUX            FALSE
+#EXP_TRANS          TRUE
+#GRND_FLUX_TYPE     GF_410
+
+#TFALLBACK          TRUE
+#SPATIAL_FROST      FALSE   (Nfrost)
+
+#######################################################################
+# Precip (Rain and Snow) Parameters
+# Generally these default values do not need to be overridden
+#######################################################################
+#SNOW_DENSITY       DENS_BRAS
+#BLOWING            FALSE
+#GLACIER            FALSE  
+#GLACIER_OVERFLOW   TRUE
+#COMPUTE_TREELINE   FALSE
+#DIST_PRCP          FALSE
+#PREC_EXPT          0.6
+#CORRPREC           FALSE
+#MAX_SNOW_TEMP      0.5
+#MIN_RAIN_TEMP      -0.5
+
+#######################################################################
+# Turbulent Flux Parameters
+# Generally these default values do not need to be overridden
+#######################################################################
+#MIN_WIND_SPEED     0.1
+#AERO_RESIST_CANSNOW    AR_406_FULL
+
+#######################################################################
+# Meteorological Forcing Disaggregation Parameters
+# Generally these default values do not need to be overridden
+#######################################################################
+#OUTPUT_FORCE       FALSE
+#PLAPSE             TRUE
+#SW_PREC_THRESH     0
+#MTCLIM_SWE_CORR    TRUE
+#VP_ITER            VP_ITER_ANNUAL      
+#VP_INTERP          TRUE
+#LW_TYPE            LW_PRATA    
+#LW_CLOUD           LW_CLOUD_DEARDORFF
+
+#######################################################################
+# Miscellaneous Simulation Parameters
+# Generally these default values do not need to be overridden
+#######################################################################
+#CONTINUEONERROR    TRUE
+#LOG_MATRIC         FALSE
+#SPATIAL_SNOW       FALSE
+
+#######################################################################
+# State Files and Parameters
+#######################################################################
+#INIT_STATE         (put the initial state path/filename here)
+#STATENAME          (put the path/prefix of output state file here)
+#STATEYEAR          2000
+#STATEMONTH         12
+#STATEDAY           31
+#BINARY_STATE_FILE  FALSE
+
+#######################################################################
+# Forcing Files and Parameters
+#######################################################################
+FORCING1            (put the forcing path/prefix here)
+FORCE_FORMAT        BINARY
+FORCE_ENDIAN        LITTLE
+N_TYPES             4
+FORCE_TYPE          PREC    UNSIGNED    40
+FORCE_TYPE          TMAX    SIGNED      100
+FORCE_TYPE          TMIN    SIGNED      100
+FORCE_TYPE          WIND    SIGNED      100
+FORCE_DT            24
+FORCEYEAR           2000
+FORCEMONTH          01
+FORCEDAY            01
+FORCEHOUR           00
+GRID_DECIMAL        4
+WIND_H              10.0   
+MEASURE_H           2.0    
+ALMA_INPUT          FALSE
+
+#######################################################################
+# Land Surface Files and Parameters
+#######################################################################
+SOIL                (put the soil parameter path/file here)
+ARC_SOIL            FALSE  
+#SOIL_DIR           (soil param directory)  
+BASEFLOW            ARNO
+JULY_TAVG_SUPPLIED  FALSE
+ORGANIC_FRACT       FALSE
+VEGLIB              (put the veg library path/file here)
+VEGPARAM            (put the veg parameter path/file here)
+ROOT_ZONES          3
+VEGPARAM_LAI        TRUE   
+LAI_SRC             LAI_FROM_VEGLIB   
+SNOW_BAND           1
+#EQUAL_AREA         FALSE
+#RESOLUTION         0.125
+
+#######################################################################
+# Lake Simulation Parameters
+# These need to be un-commented and set to correct values only when running lake model (LAKES is not FALSE)
+#######################################################################
+#LAKES              (put lake parameter path/file here)
+#LAKE_PROFILE       FALSE
+
+#######################################################################
+# Output Files and Parameters
+#######################################################################
+RESULT_DIR          (put the result directory path here)
+OUT_STEP            0      
+SKIPYEAR            0
+COMPRESS            FALSE
+BINARY_OUTPUT       FALSE
+ALMA_OUTPUT         FALSE
+MOISTFRACT          FALSE
+PRT_HEADER          FALSE  
+PRT_SNOW_BAND       FALSE  
+
+#######################################################################
+# Output File Contents
+#######################################################################

samples/global.param.sample

@@ -48,6 +48,8 @@ FROZEN_SOIL	FALSE	# TRUE = calculate frozen soils.  Default = FALSE.
 #######################################################################
 #SNOW_DENSITY	DENS_BRAS	# DENS_BRAS = use traditional VIC algorithm taken from Bras, 1990; DENS_SNTHRM = use algorithm taken from SNTHRM model.
 #BLOWING		FALSE	# TRUE = compute evaporative fluxes due to blowing snow
+#GLACIER        FALSE   # DYNAMIC = use dynamic glacier model, SCALING = use V-A scaling glacier model
+#GLACIER_OVERFLOW   TRUE
 #COMPUTE_TREELINE	FALSE	# Can be either FALSE or the id number of an understory veg class; FALSE = turn treeline computation off; VEG_CLASS_ID = replace any overstory veg types with the this understory veg type in all snow bands for which the average July Temperature <= 10 C (e.g. "COMPUTE_TREELINE 10" replaces any overstory veg cover with class 10)
 #CORRPREC	FALSE	# TRUE = correct precipitation for gauge undercatch
 #MAX_SNOW_TEMP	0.5	# maximum temperature (C) at which snow can fall
@@ -170,15 +172,15 @@ ROOT_ZONES      3	# Number of root zones (must match format of veg param file)
 #ALB_SRC 	FROM_VEGLIB    # FROM_VEGPARAM = read albedo from veg param file; FROM_VEGLIB = read albedo from veg library file
 #VEGCOVER_SRC 	FROM_VEGLIB    # FROM_VEGPARAM = read veg_cover from veg param file; FROM_VEGLIB = read veg_cover from veg library file
 SNOW_BAND	1	# Number of snow bands; if number of snow bands > 1, you must insert the snow band path/file after the number of bands (e.g. SNOW_BAND 5 my_path/my_snow_band_file)
+#EQUAL_AREA FALSE   # TRUE = grid cells are from an equal-area projection; FALSE = grid cells are on a regular lat-lon grid.  Only needed when either LAKES or GLACIER is being used.  
+#RESOLUTION 0.125   # Grid cell resolution (degrees if EQUAL_AREA is FALSE, km^2 if EQUAL_AREA is TRUE); ignored if LAKES is FALSE.  Only needed when either LAKES or GLACIER is being used. 
 
 #######################################################################
 # Lake Simulation Parameters
 # These need to be un-commented and set to correct values only when running lake model (LAKES is not FALSE)
 #######################################################################
 #LAKES		(put lake parameter path/file here)	# Lake parameter path/file
 #LAKE_PROFILE	FALSE	# TRUE = User-specified depth-area parameters in lake parameter file; FALSE = VIC computes a parabolic depth-area profile
-#EQUAL_AREA	FALSE	# TRUE = grid cells are from an equal-area projection; FALSE = grid cells are on a regular lat-lon grid
-#RESOLUTION	0.125	# Grid cell resolution (degrees if EQUAL_AREA is FALSE, km^2 if EQUAL_AREA is TRUE); ignored if LAKES is FALSE
 
 #######################################################################
 # Output Files and Parameters

src/CalcBlowingSnow.c

@@ -32,13 +32,7 @@
 
 static char vcid[] = "$Id$";
 
-#define GRAMSPKG 1000.
-#define CH_WATER 4186.8e3
-#define JOULESPCAL     4.1868   /* Joules per calorie */
-#define Ka  .0245187            /* thermal conductivity of air (W/mK) */
-#define CSALT 0.68              /* saltation constant m/s */
 #define UTHRESH 0.25            /* threshold shear velocity m/s */
-#define KIN_VIS 1.3e-5          /* Kinemativ viscosity of air (m2/s) */
 #define MAX_ITER 100             /* Max. iterations for numerical integration */
 #define K 5
 #define MACHEPS 1.0e-6          /* Accuracy tolerance for numerical integration */
@@ -146,7 +140,7 @@ double CalcBlowingSnow( double Dt,
   double Transport;
   int count=0;
 
-  Lv = (2.501e6 - 0.002361e6 * Tsnow);
+  Lv = calc_latent_heat_of_vaporization(Tsnow);
   /*******************************************************************/
   /* Calculate some general variables, that don't depend on wind speed. */
 
@@ -156,7 +150,7 @@ double CalcBlowingSnow( double Dt,
   /* Saturation density of water vapor, Liston A-8 */
   es = svp(Tair);
 
-  Tk = Tair  + KELVIN;
+  Tk = Tair + KELVIN;
 
   Ros = 0.622*es/(287*Tk);
 

src/IceEnergyBalance.c

@@ -201,9 +201,7 @@ double IceEnergyBalance(double TSurf, va_list ap)
     *Ra_used = HUGE_RESIST;
 
   /* Calculate longwave exchange and net radiation */
-
-  *LongRadOut = LongRadIn - STEFAN_B * (TMean+273.15) * (TMean+273.15) 
-    * (TMean+273.15) * (TMean+273.15);
+  *LongRadOut = LongRadIn - calc_outgoing_longwave(TMean+KELVIN, EMISS); 
   NetRad = ShortRad + *LongRadOut;
 
   /* Calculate the sensible heat flux */

src/LAKE.h

@@ -35,11 +35,6 @@
 
 #define LAKE_SET
 #define TMELT 0.0
-#define EMICE 0.97      /* Ice emissivity */
-#define EMH2O .98
-#define RHOSNOW   250.  /* densities of water and snow */
-#define RHOICE   917.   /* ice density*/
-#define rhosurf 1.275   /* surface air density */
 #define MAX_SURFACE_LAKE   .6  /* max. surface layer thickness for E-B (m) */
 #define BETA 0.001       /* Curve shape parameter for lake profile. */
 #define FRACMIN  0.10   /* min ice thickness in meters */
@@ -65,7 +60,6 @@
 #define  a1 0.7        /* Percent of radiation in visible band. */
 #define  a2 0.3        /* Percent of radiation in infrared band. */
 #define QWTAU 86400./2.   /* D. Pollard sub-ice time constant. */
-#define RADIUS 6371.228 /* Earth radius in km. */
 
 //#endif // LAKE_SET
 

src/Makefile

@@ -105,15 +105,15 @@ OBJS =  CalcAerodynamic.o CalcBlowingSnow.o SnowPackEnergyBalance.o \
 	display_current_settings.o estimate_T1.o faparl.o free_all_vars.o \
 	free_vegcon.o frozen_soil.o full_energy.o func_atmos_energy_bal.o \
 	func_atmos_moist_bal.o func_canopy_energy_bal.o \
-	func_surf_energy_bal.o get_dist.o get_force_type.o get_global_param.o \
+	func_surf_energy_bal.o get_dist.o get_force_type.o get_global_param.o glacier.o \
 	initialize_atmos.o initialize_model_state.o \
 	initialize_global.o initialize_snow.o \
 	initialize_soil.o initialize_veg.o latent_heat_from_snow.o \
 	make_cell_data.o make_all_vars.o make_dmy.o make_energy_bal.o \
 	make_in_and_outfiles.o make_snow_data.o make_veg_var.o massrelease.o \
 	modify_Ksat.o mtclim_vic.o mtclim_wrapper.o newt_raph_func_fast.o \
 	nrerror.o open_file.o open_state_file.o \
-	output_list_utils.o parse_output_info.o penman.o photosynth.o \
+	output_list_utils.o parse_output_info.o penman.o photosynth.o physics.o\
 	prepare_full_energy.o print_library.o put_data.o \
 	read_atmos_data.o read_forcing_data.o read_initial_model_state.o \
 	read_snowband.o read_soilparam.o read_veglib.o \
@@ -147,7 +147,7 @@ full:
 	make model
 
 clean::
-	/bin/rm -f *.o core log *~
+	/bin/rm -f *.o core log vicNl DEPEND *~
 
 model: $(OBJS)
 	$(CC) -o vicNl$(EXT) $(OBJS) $(CFLAGS) $(LIBRARY)

src/SnowPackEnergyBalance.c

@@ -247,25 +247,25 @@ double SnowPackEnergyBalance(double TSurf, va_list ap)
   /* Calculate longwave exchange and net radiation */
 
   Tmp = TMean + KELVIN;
-  (*NetLongUnder) = LongSnowIn - STEFAN_B * Tmp * Tmp * Tmp * Tmp;
+  (*NetLongUnder) = LongSnowIn - calc_outgoing_longwave(Tmp, EMISS);
   NetRad = NetShortUnder + (*NetLongUnder);
 
   /* Calculate the sensible heat flux */
 
-  *SensibleHeat = AirDens * Cp * (Tair - TMean) / Ra_used[0];
+  *SensibleHeat = calc_sensible_heat(AirDens, Tair, TMean, Ra_used[0]);
 
   if (options.SPATIAL_SNOW) {
     /* Add in Sensible heat flux turbulent exchange from surrounding 
        snow free patches - if present */
-    if ( SnowCoverFract > 0 ) {
+    if ( SnowCoverFract > 0. ) {
       *(AdvectedSensibleHeat) = advected_sensible_heat(SnowCoverFract, 
 						       AirDens, Tair, TGrnd, 
 						       Ra_used[0]);
     }
-    else (*AdvectedSensibleHeat) = 0;
+    else (*AdvectedSensibleHeat) = 0.;
   }
   else {
-    (*AdvectedSensibleHeat) = 0;
+    (*AdvectedSensibleHeat) = 0.;
   }
 
   /* Convert sublimation terms from m/timestep to kg/m2s */
@@ -290,7 +290,7 @@ double SnowPackEnergyBalance(double TSurf, va_list ap)
   /* Calculate advected heat flux from rain 
      Equation 7.3.12 from H.B.H. for rain falling on melting snowpack */
 
-  if ( TMean == 0 )
+  if ( TMean == 0. )
     *AdvectedEnergy = (CH_WATER * (Tair) * Rain) / (Dt);
   else
     *AdvectedEnergy = 0.;
@@ -304,7 +304,7 @@ double SnowPackEnergyBalance(double TSurf, va_list ap)
     *GroundFlux = 2.9302e-6 * SnowDensity * SnowDensity
         * (TGrnd - TMean) / SnowDepth / (Dt);
   }
-  else *GroundFlux=0;
+  else *GroundFlux = 0.;
   *DeltaColdContent -= *GroundFlux;
 
   /* Calculate energy balance error at the snowpack surface */

src/advected_sensible_heat.c

@@ -27,7 +27,7 @@ double advected_sensible_heat(double SnowCoverFract,
   double val;
 
   /* Compute sensible heat flux from bare patches */
-  Qbare = AirDens * Cp * (Tair - TGrnd) / Ra;
+  Qbare = calc_sensible_heat(AirDens, Tair, TGrnd, Ra);
 
   /* Compute fraction of sensible heat that contributes to snowmelt */
   if ( SnowCoverFract > 0.6 ) Fs = 1.0;

src/arno_evap.c

@@ -183,10 +183,10 @@ double arno_evap(layer_data_struct *layer,
 
   /* only consider positive evaporation; we won't put limits on condensation */
   if (evap > 0.0) {
-    if (moist > moist_resid * depth1 * 1000.) {
+    if (moist > moist_resid * depth1 * MMPERMETER) {
       /* there is liquid moisture available; cap evap at available liquid moisture */
-      if (evap > moist -  moist_resid * depth1 * 1000.) {
-        evap = moist -  moist_resid * depth1 * 1000.;
+      if (evap > moist -  moist_resid * depth1 * MMPERMETER) {
+        evap = moist -  moist_resid * depth1 * MMPERMETER;
       }
     }
     else {
@@ -196,7 +196,7 @@ double arno_evap(layer_data_struct *layer,
   }
 
   layer[0].evap = evap;
-  Evap += evap / 1000. / delta_t;
+  Evap += evap / MMPERMETER / delta_t;
 
   return(Evap);
 

src/calc_longwave.c

@@ -70,7 +70,7 @@ void calc_longwave(double *longwave,
     emissivity = cloudfactor*emissivity_clear;
   }
 
-  *longwave = emissivity * STEFAN_B * air_temp * air_temp * air_temp * air_temp; 
+  *longwave = calc_outgoing_longwave(air_temp, emissivity); 
 
 }