From 4ae1b578fc050cff165251f6da634c3ae0ec0b2f Mon Sep 17 00:00:00 2001 From: vicadmin Date: Fri, 7 May 2004 00:32:27 +0000 Subject: [PATCH] This commit was manufactured by cvs2git to create tag 'Rel_4_0_3'. Sprout from master 2000-08-16 00:35:28 UTC vicadmin 'Converted char flags to int' Cherrypick from master 2004-05-07 00:32:26 UTC tbohn '': README.txt Delete: aurad.c calc_long_shortwave.c calc_netshort.c calc_snow_ground_flux.c calc_trans.c calc_water_density.c func_snow_ground_flux.c initialize_energy.c long_shortwave.c make_out_data.c makefile polint.c priestley.c rad_and_vpd.c rad_and_vpd.h read_PILPS2c.c read_atmosdata.c read_initial_snow.c read_initial_soil_thermal.c read_rosemount.c read_sawd.c read_sawd_binary.c read_snowmodel.c --- README.txt | 247 ++++++++++++++++ aurad.c | 447 ---------------------------- calc_long_shortwave.c | 204 ------------- calc_netshort.c | 36 --- calc_snow_ground_flux.c | 367 ----------------------- calc_trans.c | 22 -- calc_water_density.c | 48 --- func_snow_ground_flux.c | 173 ----------- initialize_energy.c | 561 ------------------------------------ long_shortwave.c | 88 ------ make_out_data.c | 18 -- makefile | 238 --------------- polint.c | 38 --- priestley.c | 38 --- rad_and_vpd.c | 50 ---- rad_and_vpd.h | 61 ---- read_PILPS2c.c | 99 ------- read_atmosdata.c | 90 ------ read_initial_snow.c | 87 ------ read_initial_soil_thermal.c | 162 ----------- read_rosemount.c | 92 ------ read_sawd.c | 120 -------- read_sawd_binary.c | 105 ------- read_snowmodel.c | 94 ------ 24 files changed, 247 insertions(+), 3238 deletions(-) create mode 100644 README.txt delete mode 100644 aurad.c delete mode 100644 calc_long_shortwave.c delete mode 100644 calc_netshort.c delete mode 100644 calc_snow_ground_flux.c delete mode 100644 calc_trans.c delete mode 100644 calc_water_density.c delete mode 100644 func_snow_ground_flux.c delete mode 100644 initialize_energy.c delete mode 100644 long_shortwave.c delete mode 100644 make_out_data.c delete mode 100644 makefile delete mode 100644 polint.c delete mode 100644 priestley.c delete mode 100644 rad_and_vpd.c delete mode 100644 rad_and_vpd.h delete mode 100644 read_PILPS2c.c delete mode 100644 read_atmosdata.c delete mode 100644 read_initial_snow.c delete mode 100644 read_initial_soil_thermal.c delete mode 100644 read_rosemount.c delete mode 100644 read_sawd.c delete mode 100644 read_sawd_binary.c delete mode 100644 read_snowmodel.c diff --git a/README.txt b/README.txt new file mode 100644 index 000000000..e68a7cfc3 --- /dev/null +++ b/README.txt @@ -0,0 +1,247 @@ +# README.txt - Release Notes +# +# $Id$ +#------------------------------------------------------------------------ + +August 15, 2000: VIC release 4.0.3 + +This release fixes a problem with the implementation of va_arg that +causes run time errors on some systems. Previous releases of the code +worked correctly on the LINUX and freeBSD systems where it was tested. +However on some systems (including Sun Ultra-2s) character variables +passed with va_arg are changed into integers so reading a character from +the argument list does not produce the value sent to the routine. The +character flags used by VIC to indicate if there is snow present and if +the frozen soil algorithm has been activated have now been converted to +integers, which should make the va_arg call work on all systems. + +Also fixed in this release was a check in dist_prec.c to see if it is +still raining which actually used the memory address of the +precipitation variable rather than the daily value in the check. + +MODIFIED FILES: + read_atmos_data.c - Fixed input file time step check + write_forcing_files.c - Added free statements for pointers + calc_surf_energy_bal.c - Converted char flags to int + dist_prec.c - Fixed logical statement error + frozen_soil.c - Converted char flags to int + func_surf_energy_bal.c - Converted char flags to int + initialize_atmos.c - Added flag for output forcing + vicNl.h - Converted char flags to int + vicNl_def.h - Converted char flags to int + +July 19, 2000: VIC release 4.0.2 + +Two new pre-processor options have been added to VIC as well as minor +modifications to two subroutines. + +If set to TRUE the NO_REWIND pre-processor option stops the VIC model from +rewinding the soil and vegetation parameter input files for each new grid +cell. This reduces run times but requires that all input files are in the +same order as the soil parameter file. + +If set TRUE the OUTPUT_FORCE pre-processor option blocks off the main +model and only reads the provided forcing files. Once VIC has estimated +the missing forcing parameters the full forcing data set for the defined +simulation period is output to a series of gridded forcing files. The +gridded forcing files are written to the RESULTS directory defined in the +global control file with the prefix "full_data_". The new files are in +Binary or ASCII depending on the setting of BINARY_OUTPUT. + +The error messages in get_global_param.c have been modified so that the +correct file is referenced when telling the user to change values found in +the model source code. + +In read_soilparam.c, the soil parameters are defined only if the current +grid cell is run, otherwise the line in the file is skipped and soil_con +is retruned without new data values. + +May 30, 2000: VIC release 4.0.1 + +Increased use of the released VIC model code has lead to the +discovery of a couple of minor bugs. This release fixes those bugs as +well as introducing a improved precipitation correction algorithm based on +Yang et al. 1998. Unless you have encountered these problems or are +trying to correct precipitation undercatch due to wind, in the VIC +model, your results will not be impacted by these fixes. + +MODIFIED FILES: + correct_precip.c - changed to WMO correction equation for + NWS 8" standard gauge. + full_energy.c - modified to handle WMO correction. + initialize_atmos.c - modified to handle WMO correction. + fixed error in estimating minimum daily + temperature from sub-daily temperatures. + make_in_and_outfiles.c - removed line that opened the state file + again for each new grid cell. + open_state_file.c - modified comments. + put_data.c - modified to handle WMO correction. + snow_utility.c - cleaned source code. + solve_snow.c - modified to handle WMO correction. + surface_fluxes.c - modified to handle WMO correction. + vicNl.h - modified to handle WMO correction. + vicNl_def.h - modified to handle WMO correction. + +REFERENCE: + + Yang, D., B. E. Goodison, J. R. Metcalfe, V. S. Golubev, R. + Bates, T. Pangburn, and C. L. Hanson, Accuracy of NWS 8" Standard + Nonrecording Precipitation Gauge: Results and Application of WMO + Intercomparison, J. Atmos. Oceanic Tech., 15, 54-68, 1998. + + +Date: May 16, 2000 + +From: Keith Cherkauer + +Topic: Release of VIC 4.0.0 + +The code for VIC release 4.0.0 has undergone several months of tests (as +VIC release 3.3.0 Beta) and has now been deemed ready for release to the +general public. This document is designed to provide information +concerning changes in the model between the last release version (3.2.1) +and the current version. + +There is no formal users manual, information about how to use the current +version can be found at +http://www.hydro.washington.edu/Lettenmaier/Models/VIC/VIChome.html. +Information about the basic model design can be found in Liang, et al. +(1994), while the rewrite of the source code as well as the addition of +cold season processes is described in Cherkauer and Lettenmaier (1999). + +The VIC macroscale hydrologic model was developed as a research tool. As +such it is the users responsibility to verify that it is working correctly +when applied to new situations. When using the VIC model reference should +be made to Liang, et al. (1994) and Cherkauer and Lettenmaier (1999) as +well as an acknowledgment that the code was received from the University +of Washington. Other important papers relating to the development of the +VIC model are included on the home page and in the source code. + +Possible bugs in the code should be reported to +vicadmin@hydro.washington.edu. ALWAYS CHECK YOUR INPUT FILES! Most +"bugs" are actually caused by trying to run the model with bad parameters +or forcing data. The VIC model will run limited checks to find common +major errors but in most cases it will attempt to run with the bad values. +If after checking all of your input data you still believe you have found +a bug in the model, send an e-mail including the complete output from the +model as well as a description of the problem and the files necessary to +run the model to recreate the code (if files are large please put a +compressed tar file in +ftp://ftp.ce.washington.edu/pub/HYDRO/vicadmin/TEMP). Outdated and +modified versions of the code are the responsibility of the user to debug. +Modifications made to the code, which may improve the general model +performance, may be submitted to vicadmin@hydro.washington.edu for +possible inclusion in future versions of the code. + +VIC release 4.0.0 represents a major change to the source code from +version 3.2.1. It is strongly recommended that if you were using version +3.2.1 or earlier versions that you update with a complete copy of the new +code. + +Major changes from release version 3.2.1 to 4.0.0: + +- Radiation Estimation Update: The routines to estimate shortwave and +longwave radiation as well as vapor pressure from daily minimum and +maximum temperature and precipitation have been updated to correspond to +the algorithm described by Thornton and Running (1999). These routines +provide significantly improved radiation estimates especially in regions +outside the continental United States. + +- Model Core Update: The core of the VIC model was rewritten so that all +modes of the model (water balance, energy balance, etc.) make use of the +same model code. This makes it easier to modify the model and have +modifications apply to sll modes, it also allows the model to run with new +combinations of algorithms (i.e. full energy balance mode with the finite +difference ground heat flux solution). + +- Soil Moisture Transport Update: The frozen and thawed sub-layers added +to the model for the original frozen soil algorithm have been removed. +This makes the soil drainage routine cleaner and faster. Frozen soils now +estimate full layer ice contents from the ice content at each soil thermal +node. Without being confined by sub-layers, the frozen soil algorithm can +now be applied to regions of permafrost. + +- Forcing File Control Added: Version 4.0.0 moves controls of the forcing +file format and data types into the global control file. The model can +now handle most ASCII column and short int Binary files without writing +new subroutines and recompiling the source code. + +- Pre-processor Options Added: There are now more option flags in the +source code headers to control which parts of the model are in fact +compiled. This allows the model functionality to be adjusted without the +addition of computationally intensive conditional switching statements. + +- Model State File: With the release of version 4.0.0 separate snow and +soil initialization files have been combined into a single model state +file. The state file can be created outside the model for starting +simulations with prescribed initial conditions, or the model state can be +saved by VIC at a specified date. Note that currently there will be small +differences between a full and a warm started simulation because radiation +and vapor pressure are estimated using forcing data from the simulation +period, not from the full dataset included in the forcing file. It also +does not store wet and dry fraction information, when running with +distributed precipitation the model is restarted using average soil +conditions. + +References: + +Liang, X., D. P. Lettenmaier, E. F. Wood, and S. J. Burges, A simple +hydrologically based model of land surface water and energy fluxes for +GSMs, J. Geophys. Res., 99(D7), 14,415-14,428, 1994. + +Cherkauer, K. A., and D. P. Lettenmaier, Hydrologic effects of frozen +soils in the upper Mississippi River basin, J. Geophys. Res., 104(D16), +19,599-19,610, 1999. + +Thornton, P.E. and S.W. Running, An improved algorithm for estimating +incident daily solar radiation from measurements of temperature, humidity, +and precipitation, Ag. For. Met., 93, 211-228, 1999. + +File List: + +CalcAerodynamic.c modify_Ksat.c +Makefile mtclim42_vic.c +SnowPackEnergyBalance.c mtclim42_vic.h +StabilityCorrection.c mtclim42_wrapper.c +alloc_atmos.c nrerror.c +arno_evap.c open_debug.c +calc_air_temperature.c open_file.c +calc_cloud_cover_fraction.c open_state_file.c +calc_longwave.c penman.c +calc_rainonly.c prepare_full_energy.c +calc_root_fraction.c put_data.c +calc_surf_energy_bal.c read_arcinfo_ascii.c +calc_veg_params.c read_atmos_data.c +canopy_evap.c read_forcing_data.c +check_files.c read_initial_model_state.c +check_state_file.c read_snowband.c +close_files.c read_soilparam.c +cmd_proc.c read_soilparam_arc.c +compress_files.c read_veglib.c +compute_dz.c read_vegparam.c +correct_precip.c redistribute_during_storm.c +dist_prec.c root_brent.c +estimate_T1.c runoff.c +free_dist_prcp.c snow.h +free_vegcon.c snow_intercept.c +frozen_soil.c snow_melt.c +full_energy.c snow_utility.c +func_surf_energy_bal.c soil_conduction.c +get_force_type.c soil_thermal_eqn.c +get_global_param.c solve_snow.c +global.h store_moisture_for_debug.c +initialize_atmos.c surface_fluxes.c +initialize_global.c svp.c +initialize_model_state.c user_def.h +initialize_new_storm.c vicNl.c +initialize_snow.c vicNl.h +initialize_soil.c vicNl_def.h +initialize_veg.c vicerror.c +make_cell_data.c write_atmosdata.c +make_dist_prcp.c write_data.c +make_dmy.c write_debug.c +make_energy_bal.c write_layer.c +make_in_and_outfiles.c write_model_state.c +make_snow_data.c write_soilparam.c +make_veg_var.c write_vegparam.c +massrelease.c write_vegvar.c diff --git a/aurad.c b/aurad.c deleted file mode 100644 index 1c3d1e36b..000000000 --- a/aurad.c +++ /dev/null @@ -1,447 +0,0 @@ -/******************************************************************************** - filename : aurad.c - purpose : Calculate daily net radiation - interface : see below - programmer: see below - date : August 4, 1995 - changes : changed include file from bgc.h to rad_and_vpd.h (08/04/1995) - notes : original version received August 4, 1995 from John Kimball, - NTSG, School of Forestry, University of Montana - - unexplained (in list below) variables: - tranday - daily total transmittance - - references: - Buffo, J., L. Fritschen, and J. Murphy, Direct solar radiation on - various slopes from 0 to 60 north latitude. USDA Forest Service - Paper Research Paper PNW-142, 1972. - Swift, L. W. Jr., Algortihm for solar radiation on mountain slopes. - Water Resources Research, v.12, p.108-112, 1976. - - for daily total transmittance: - Bristow, K. L., and G. S. Campbell, On the relationship between - incoming solar radiation and daily maximum and minimum temperature. - Agricultural and Forest Meteorology, v.31, p.159-166, 1984. -********************************************************************************/ - - - -/* ****************************************************** -c * AURAD.F * -c ****************************************************** - -CM CONTAINS: fltrad, shrad -c -c ****************************************************** -c * subroutine fltrad * -c * * -C ****************************************************** -CF FILE : aurad.f -CM CONTAINS: fltrad -CMA AUTHOR : jc coughlan -CMS SYSTEM : autosys -CMV VERSION : #2 -CMX STATUS : -CVF VFILE : -CVI VINPUT : -CVO VOUTPUT : -CR RETURNS : -CD DESCRIP : given below -CH HISTORY : CREATED 5/2/1989 -ch modified 5/ /1990 r nemani jp patterson -ch modified 11/18/98 k cherkauer - fixed comments about UNITS -cd flat surface radiation computation. -cd slope = 0 aspect = 0 -cd this subroutine computes incident shortwave radiation and -cd net shortwave radiation for any given day based on -cd sun-earth geometry, transmissivity; -cd derived from: buffo et al 1972; swift 1976. -cd -cd variable list: -cd sol=solar constant derived from solcon array -cd am=optical air mass for angles > 21 degrees -cd a=optical air mass array for angles between -cd 0 and 21 degrees above horizon -cd decl=declination -cd jday=day of year -cd h=angle of sun from solar noon -cd tran=transmissivity constant -cd tram=transmissivity corrected for air mass -cd nnh=calculation interval in seconds -cd nc=seconds in one day (24 hours) -cd n=number of intervals of length nnh in one day -cd dt=direct solar perpendicular to sun on the -cd outside of atmosphere for interval (kj/m**2) -cd etf=direct solar on outside of atmosphere -cd parallel to earths surface for interval -cd grad=total daily radiation at given location (kj/m**2) -cd hrad=direct solar on earths surface (flat) " -cd tdif=total daily diffuse radiation " -cd difrad=diffuse on slope for interval " -cd drad=direct on slope for interval " -cd cza=cosine zenith angle -cd cbsa=cosine beam slope angle -cd globf=global radiation, determining diffuse -cd diffl=diffuse on flat for interval -cd dayl=daylength (hours) -cd albdo=albedo never used -cd arad=absorbed radiation, kJ/day/m2 removed -cd radn=incident radiation kJ/day/m2 -cd (* . . . . . . . . . . . . . . . . . . . . . . . . . . . */ - -#include -#include -#include - - -static char vcid[] = "$Id$"; - -/* Since K&R C doesn't allow multiple assignments of arrays within functions - this must be done as an external variable, but that's ok since this way - both functions can use the data */ - - double aa[22] = - {0.0, 2.90,3.05,3.21,3.39, 3.69, 3.82, 4.07, 4.37, 4.72, 5.12, - 5.60,6.18,6.88,7.77,8.90,10.39,12.44,15.36,19.79,26.96,30.00}; - - double dec[47] = - { 0.0,-23.0,-22.0,-21.0,-19.0,-17.0,-15.0,-12.0,-9.0, -6.0, -3.0, - 0.0, 3.0, 6.0, 9.0, 12.0, 14.0, 17.0, 19.0, 21.0, 22.0, 23.0, - 23.5, 23.5, 23.0, 21.5, 20.0, 18.0, 16.0, 14.0, 12.0, 9.0, 6.0, - 3.0, 0.0, -3.0, -6.0, -9.0,-12.0,-15.0,-17.0,-19.0, -21.0,-22.0, - -23.0,-23.5,-23.5}; - - double solcon[13] = - { 0.0, 1.445,1.431,1.410,1.389,1.368,1.354, - 1.354,1.375,1.403,1.424,1.438,1.445}; - - -double fltrad(double sehorz, - double swhorz, - double slat, - int jday, - double tranday, - double *dayl) -{ - int i, nnh, nc, n, nh, ml, mo, idec; - - double decl,h, grad, radn; - double xtran, sol, tdif, tdrad, dayl2; - double cosdec, cosdla, sindla, sindec, coshh; - double cza, am, tram, dt, etf, hrad, globf, diffl, cbsa, drad, se; - double difrad; - - xtran=tranday; - nnh=3600; - nc=86400; - n = nc/nnh + 1; - dayl[0]=0.0; - - - mo = jday/30 + 1; - mo = (mo<12) ? mo:12; - - /* solcon array is in units of kw/m**2 */ - - sol=solcon[mo]; - idec= (int) (1.0 + (double) jday/8.0); - decl = dec[idec]; - grad = 0.0; - tdif=0.0; - tdrad=0.0; - nh=0; - dayl2=0.0; - - cosdec = cos(decl*DtoR); - cosdla = cos(slat*DtoR); - sindla = sin(slat*DtoR); - sindec = sin(decl*DtoR); - - for (i=1; i<=n; i++) { - nh = nh+nnh; - - /* determine angle from solar noon */ - - h=(double) (nh-43200)*0.0041667; - coshh = cos(h*DtoR); - cza=cosdec*cosdla*coshh+sindec*sindla; - - if (cza > 0.0) { - - /* (* daylength based on solar elevation above a flat horizon *) */ - - dayl2=dayl2+((double) (nnh)/3600.0); - - /* (* determine optical air mass *) */ - - am=1.0/(cza+0.0000001); - - if (am > 2.9) { - ml= ((int) (acos(cza)/0.0174533)) - 69; - - if(ml < 1) - ml=1; - - if(ml > 21) - ml=21; - - am = aa[ml]; - - } - - tram = pow(xtran,am); - dt = sol* (double) nnh; - etf=cza*dt; - hrad=etf*tram; - dt=dt*tram; - globf=sqrt(hrad*etf); - diffl=globf*(1.0-globf/etf); - - /* (* compute flat surface angle */ - - cbsa = cosdla*cosdec*coshh + sindla*sindec; - - /* (* compute direct incident radiation if solar angle is positive */ - - if(cbsa >= 0.0) { - drad=cbsa*dt; - - /* (* the following three lines computes a topographic - c reduction of direct radiation - c ehe = east horizon elevation (degrees) - c whe = west horizon elevation (degrees) */ - - se = 1.5715-(acos(cza)); - - if (h < 0.0 && se < sehorz *DtoR) - drad=0.0; - - if(h > 0.0 && se < swhorz *DtoR) - drad=0.0; - } - else - drad=0; - - /* the code below is a simplification of the - c commented line assuming a flat surface. */ - - difrad = diffl; - - /* difrad=diffl*(cos(dslop*.5*DtoR))**2.) */ - - dayl[0] =dayl[0]+((double) (nnh)/3600.0); - grad =grad+drad+difrad; - tdif =tdif+difrad; - tdrad =tdrad+drad; - } - } - - /* averad=(grad/(dayl*3600.0)); */ - - radn = grad; - /* arad = averad; */ - - return((double) radn); -} - - - - -/* ****************************************************** -c * subroutine shrad * -c * * -C ****************************************************** -CF FILE : aurad.f -CM CONTAINS: shrad -CMA AUTHOR : jc coughlan -CMS SYSTEM : autosys -CMV VERSION : #1 -CMX STATUS : -CVF VFILE : -CVI VINPUT : -CVO VOUTPUT : -CR RETURNS : -CD DESCRIP : given below -CH HISTORY : CREATED 5/2/1989 -ch modified 5/ /1990 r nemani and jp patterson -c - -cd this subroutine computes incident shortwave radiation and -cd net shortwave radiation for any given day based on surface -cd characteristics, sun-earth geometry, transmissivity; based on -cd buffo et al 1972; swift 1976. -cd -cd variable list: -cd sol=solar constant derived from solcon array -cd am=optical air mass for angles > 21 degrees -cd a=optical air mass array for angles between -cd 0 and 21 degrees above horizon -cd decl=declination -cd jday=day of year -cd asp=aspect in degrees -cd dslop=slope in degrees !!! -cd h=angle of sun from solar noon -cd tran=transmissivity constant -cd tram=transmissivity corrected for air mass -cd nnh=calculation interval in seconds -cd nc=seconds in one day (24 hours) -cd n=number of intervals of length nnh in one day -cd dt=direct solar perpendicular to sun on the -cd outside of atmosphere for interval (kj/m**2) -cd etf=direct solar on outside of atmosphere -cd parallel to earths surface for interval -cd grad=total daily radiation at given location (kj/m**2) -cd hrad=direct solar on earths surface (flat) " -cd tdif=total daily diffuse radiation " -cd difrad=diffuse on slope for interval " -cd drad=direct on slope for interval " -cd cza=cosine zenith angle -cd cbsa=cosine beam slope angle -cd globf=global radiation, determining diffuse -cd diffl=diffuse on flat for interval -cd dayl=daylength (hours) -cd albdo=albedo never used -cd arad=absorbed radiation, kJ/day/m2 removed -cd radn=incident radiation kJ/day/m2 -cd */ - -double shrad (double asp, double dslop, double slat, double sehorz, double swhorz, int jday, double tranday) -{ - int i, n, nh, nc, nnh, ml, mo, idec; - - double decl, h, grad, radn, tempf; - double xtran, dayl, sol, tdif, tdrad, dayl2; - double cosdec, sindsl, sinasp, cosdsl, cosdla, cosasp; - double sindla, sindec, sinhh, coshh, cza, am, tram, dt, etf; - double hrad, globf, diffl, cbsa, drad, se, difrad; - - /* slope in degrees - no conversion! */ - - xtran=tranday; - nnh = 3600; - nc = 86400; - n = nc/nnh + 1; - dayl=0.0; - mo = jday/30+1; - mo = (mo<12) ? mo:12; - - /* solcon array is in units of kw/m**2 */ - - sol = solcon[mo]; - idec = (int) (1.0 + (double) jday/8.0); - decl = dec[idec]; - grad = 0.0; - tdif=0.0; - tdrad=0.0; - nh = 0; - dayl2=0.0; - - tempf= (double) (decl*DtoR); - cosdec = cos((double) tempf); - /* cosdec = cos((double) decl); */ - sindsl = sin((double) (dslop*DtoR)); - sinasp = sin((double) (asp*DtoR)); - cosdsl = cos((double) (dslop*DtoR)); - cosdla = cos((double) (slat*DtoR)); - cosasp = cos((double) (asp*DtoR)); - sindla = sin((double) (slat*DtoR)); - sindec = sin((double) (decl*DtoR)); - - for (i=1;i<=n; i++) { - nh=nh+nnh; - - /* determine angle from solar noon */ - - h = (double) (nh-43200)*0.0041667; - sinhh = sin(h*DtoR); - coshh = cos(h*DtoR); - cza = cosdec*cosdla*coshh+sindec*sindla; - - if(cza > 0.0) { - - /* daylength based on solar elevation above a flat horizon */ - - dayl2=dayl2+((double) (nnh)/3600.0); - - /* next 6 lines, determine optical air mass */ - - am = 1.0/(cza+0.0000001); - - if(am > 2.9) { - ml= ((int) (acos(cza)/0.0174533)) - 69; - - if(ml < 1) - ml=1; - - if(ml > 21) - ml=21; - - am = aa[ml]; - } - - tram = pow(xtran,am); - dt=sol*(double) nnh; - etf=cza*dt; - hrad=etf*tram; - dt=dt*tram; - globf=sqrt(hrad*etf); - diffl=globf*(1.0-globf/etf); - -/* the following line is a simplification of the -c commented computation. It is kept a documentation for -c the derivation of the simpler form. */ - - cbsa=(-sindsl)*sinasp*sinhh*cosdec - + (-cosasp*sindsl*sindla+cosdsl*cosdla)*cosdec*coshh - +(cosasp*sindsl*cosdla+cosdsl*sindla)*sindec; - -/* cbsa=-sin(dslop*DtoR)*sin(asp*DtoR)* -c & sin(h*DtoR)*cos(conv(decl))+(-cos(conv(asp)) -c & *sin(conv(dslop))*sin(conv(dlat))+cos(conv(dslop))* -c & cos(conv(dlat)))*cos(conv(decl)) -c & *cos(conv(h))+(cos(conv(asp))*sin(conv(dslop))* -c & cos(conv(dlat))+cos(conv(dslop))* -c & sin(conv(dlat)))*sin(conv(decl)) */ - - if (cbsa >= 0.0) { - drad=cbsa*dt; - -/* the following three lines computes a topographic reduction of - -c direct radiation -c ehe = east horizon elevation (degrees) -c whe = west horizon elevation (degrees) */ - - se = 1.5715-(acos(cza)); - - if(h < 0.0 && se < sehorz*DtoR) - drad=0.0; - - if(h > 0.0 && se < swhorz*DtoR) - drad=0.0; - } - else - drad=0; - - difrad = diffl * pow(cos(dslop*0.5*DtoR),2.0); - -/* radt = (drad+difrad)/(double) nnh; */ - - dayl=dayl+((double) (nnh)/3600.0); - grad=grad+drad+difrad; - tdif=tdif+difrad; - tdrad=tdrad+drad; - } - } - -/* averad=(grad/(dayl*3600.0)); */ - radn = grad; -/* arad = averad; */ - - return((double) radn); -} - -/*cf ***************************************************** -cf * end of AURAD.F * -cf *****************************************************/ - diff --git a/calc_long_shortwave.c b/calc_long_shortwave.c deleted file mode 100644 index 3e1e63f82..000000000 --- a/calc_long_shortwave.c +++ /dev/null @@ -1,204 +0,0 @@ -#include -#include -#include -#include - -static char vcid[] = "$Id$"; - -void calc_long_shortwave(double *shortwave, - double *longwave, - double *tskc, - double air_temp, - double vp, - double theta_l, - double theta_s, - double phi, - double jdate, - double hour, - char HAVE_SHORTWAVE, - char HAVE_LONGWAVE, - char HAVE_TSKC) { -/********************************************************************** - calc_long_shortwave.c Keith Cherkauer March 7, 1998 - - This routine computes long and short wave radiation based on - geographic location, time and cloud cover. Equations are based on - Bras 21-47. - - shortwave incoming shortwave radiation (W/m^2) - longwave incoming longwave radiation (W/m^2) - tskc fraction of sky covered by clouds (fract) - air_temp air temperature (C) - vp vapor pressure (kPa) - theta_l defined longitude of time zone (degree) - theta_s longitude of gird cell (degree) - phi latitude of grid cell (degree) - jdate day in year of current time step (day) - hour hour of current time step (hour) - HAVE_SHORTWAVE if TRUE do not calculate shortwave - HAVE_LONGWAVE if TRUE do not calculate longwave - HAVE_TSKC if TRUE do not calculate cloud coverage - - Modifications: - 11-30-98 TSKC is a fractional value from 0 to 1, but was being - used as if it was a percent (0 to 100). Based on Gregs - testing have fixed the problem. KAC - -**********************************************************************/ - - static double last_tskc; - - double declination; - double i_var; - double tau; - double sin_alpha; - double radius; - double I0; - double m; - double Ic; - int sum_exceed=0; - int sum_zero=0; - - hour = hour - SOLARTIMEOFFSET; /* assume shortwave - measurements made during previous hour */ - declination = 23.45 * PI / 180.0 * cos(2.0 * PI / 365.0 - * (172.0 - jdate)); - if(fabs(theta_l)==theta_l) i_var = 1.0; - else i_var = -1.0; - - /** Check if sun is east or west of cell longitude **/ - if(((float)hour > (12.+(theta_l-theta_s)*24./360.) - && (float)hour < (24.+(theta_l-theta_s)*24./360.)) - || ((float)hour < (0. +(theta_l-theta_s)*24./360.))) - tau = (hour - 12.0 - (i_var/15.0 * (fabs(theta_s) - - fabs(theta_l)))) * 15.0; - else tau = (hour + 12.0 - (i_var/15.0 * (fabs(theta_s) - - fabs(theta_l)))) * 15.0; - sin_alpha = (sin(declination) * sin(phi*PI/180.0) + cos(declination) - * cos(phi*PI/180.0) * cos(tau*PI/180.0)); - radius = 1.0+0.017*cos((double)(2*PI/365*(186-jdate))); - I0 = 1353.0 * sin_alpha / radius / radius; - - if(!HAVE_SHORTWAVE || (HAVE_SHORTWAVE && *shortwave<0.0)) *shortwave=0.; - if(I0>0.0) { - m = pow((sin_alpha + 0.15*pow((asin(sin_alpha) + 3.885),-1.253)),-1.); - Ic = I0 * exp(-2.0 * (0.128 - 0.054 * log10(m)) * m); - if(HAVE_SHORTWAVE) Ic = I0; /** so far the above eqns appear - to over correct **/ - if(!HAVE_SHORTWAVE && HAVE_TSKC) { - /** Need to Calculate Shortwave Radiation **/ - *shortwave = (1.0 - 0.65 * ((*tskc) * (*tskc)) / 1.00) * Ic; - if(*shortwave < 0.0) { - sum_zero++; - *shortwave = 0.0; - } - if(!HAVE_LONGWAVE) - *longwave = (1.0 + 0.17 * (*tskc) - * (*tskc)) * (0.740+0.0049* (vp)*10.0) - * STEFAN_B * (air_temp+KELVIN) * (air_temp+KELVIN) - * (air_temp+KELVIN) * (air_temp+KELVIN) / LWAVE_COR; - } - - else if(HAVE_SHORTWAVE && !HAVE_LONGWAVE) { - - /** Shortwave Measured, Cloud Cover Needed **/ - if(*shortwave < Ic) { - *tskc = sqrt((1.0 - *shortwave / Ic) / 0.65); - last_tskc = *tskc; - *longwave = (1.0 + 0.17 * (*tskc) * (*tskc)) - * (0.740+0.0049*vp*10.0) * STEFAN_B - * (air_temp+KELVIN) * (air_temp+KELVIN) - * (air_temp+KELVIN) * (air_temp+KELVIN) / LWAVE_COR; - } - else { - /** Measured shortwave exceeds estimated **/ - sum_exceed++; - *tskc = last_tskc; - *longwave = (1.0 + 0.17*(*tskc)*(*tskc)) - * (0.740+0.0049*vp*10.0) * STEFAN_B - * (air_temp+KELVIN) * (air_temp+KELVIN) - * (air_temp+KELVIN) * (air_temp+KELVIN) / LWAVE_COR; - } - } - else { - nrerror("ERROR: To compute long and shortwave radiation, need TSKC (cloud cover fraction), or measured shortwave"); - } - } - else { - if(!HAVE_SHORTWAVE) *shortwave = 0.0; - else if(HAVE_SHORTWAVE && !HAVE_TSKC) *tskc = last_tskc; - if(!HAVE_LONGWAVE) - *longwave = (1.0 + 0.17*(*tskc)*(*tskc)) - * (0.740+0.0049*vp*10.0) * STEFAN_B - * (air_temp+KELVIN) * (air_temp+KELVIN) - * (air_temp+KELVIN) * (air_temp+KELVIN) / LWAVE_COR; - } - -} - -void get_rise_and_set_hours(int *rise_hour, - int *set_hour, - double theta_l, - double theta_s, - double phi, - double jdate) { -/********************************************************************** - get_rise_and_set_hours.c Keith Cherkauer September 8, 1998 - - This routine computes the hours at which the sun rises and sets for - the current day. Equations are based on Bras 21-47. - - rise_hour hour the sun rises above the horizon - set_hour hour the sun sets below the horizon - theta_l defined longitude of time zone (degree) - theta_s longitude of gird cell (degree) - phi latitude of grid cell (degree) - jdate day in year of current time step (day) - -**********************************************************************/ - - double declination; - double i_var; - double cos_tmp; - double sin_tmp; - double tau; - double sin_alpha; - double radius; - double radius2; - double I0; - int hour; - - rise_hour[0] = -1; - set_hour[0] = -1; - - if(fabs(theta_l)==theta_l) - i_var = 1.0; - else - i_var = -1.0; - - declination = 23.45 * PI / 180.0 * cos(2.0 * PI / 365.0 * (172.0 - jdate)); - sin_tmp = sin(declination) * sin(phi*PI/180.0); - cos_tmp = cos(declination) * cos(phi*PI/180.0); - radius = 1.0+0.017*cos((double)(2*PI/365*(186-jdate))); - radius2 = radius*radius; - - for(hour=0;hour<24;hour++) { - hour = hour - SOLARTIMEOFFSET; /* assume shortwave measurements made - during previous hour */ - - /** Check if sun is east or west of cell longitude **/ - if(((float)hour > (12.+(theta_l-theta_s)*24./360.) - && (float)hour < (24.+(theta_l-theta_s)*24./360.)) - || ((float)hour < (0. +(theta_l-theta_s)*24./360.))) - tau = (hour - 12.0 - (i_var/15.0 * (fabs(theta_s) - - fabs(theta_l)))) * 15.0; - else tau = (hour + 12.0 - (i_var/15.0 * (fabs(theta_s) - - fabs(theta_l)))) * 15.0; - sin_alpha = sin_tmp + cos_tmp * cos(tau*PI/180.0); - I0 = 1353.0 * sin_alpha / radius2; - - if(rise_hour[0] < 0 && I0 > 0) rise_hour[0] = hour; - if(set_hour[0] < 0 && rise_hour[0] > 0 && I0 < 0) set_hour[0] = hour; - - } -} diff --git a/calc_netshort.c b/calc_netshort.c deleted file mode 100644 index aa0f81d5b..000000000 --- a/calc_netshort.c +++ /dev/null @@ -1,36 +0,0 @@ -#include -#include -#include - -static char vcid[] = "$Id$"; - -double calc_netshort(double trans, - int day, - double lat, - double *day_len) -/********************************************************************** - calc_netshort Dag Lohmann January 1996 - - This routine calculates the net short wave radiation per day, using - the fltrad routine in the file aurad.c. Fltrad returns net radiation - in kJ/day/m^2, so calc_netshort converts the value into W/m^2 before - returning it. - - Routine returns net shortwave radiation at the surface in W/m^2. - - 09/04/98 Modified to get the length of day in hours from fltrad. KAC - 11/18/98 Comments amended to include references to units. KAC - -**********************************************************************/ -{ - double netshort; - double albedo; - - netshort = 1000.0/SEC_PER_DAY * fltrad((double) 0.0, (double) 0.0, - lat, day, trans, day_len); - albedo = 0.2; - netshort *= (1.0 - albedo); - - return (netshort); - -} diff --git a/calc_snow_ground_flux.c b/calc_snow_ground_flux.c deleted file mode 100644 index 9073f0064..000000000 --- a/calc_snow_ground_flux.c +++ /dev/null @@ -1,367 +0,0 @@ -#include -#include -#include - -static char vcid[] = "$Id$"; - -double calc_snow_ground_flux(int dt, - int Nnodes, - int rec, - int iveg, - double mu, - double dp, - double moist, - double ice0, - double Tsnow_surf, - double *T1, - double *grnd_flux, - double *deltaH, - double *snow_flux, - energy_bal_struct *energy, - snow_data_struct *snow, - layer_data_struct *layer_wet, - layer_data_struct *layer_dry, - soil_con_struct soil_con, - char *FIRST_SOLN) { -/********************************************************************** - calc_snow_ground_flux.c Keith Cherkauer March 31, 1997 - - This subroutine computes the surface temperature under the snowpack - using the thermal flux through the snowpack. - -**********************************************************************/ - - extern option_struct options; - - int Twidth; - - double T2; - double Ts_old; - double T1_old; - double kappa1; - double kappa2; - double Cs1; - double Cs2; - double delta_t; - double snow_density; - double snow_depth; - double surf_temp; - double D1; - double D2; - double max_moist; - double bubble; - double expt; - - double error; - double *T_node; - double Tnew_node[MAX_NODES]; - double *dz_node; - double kappa_node[MAX_NODES]; - double Cs_node[MAX_NODES]; - double moist_node[MAX_NODES]; - double *expt_node; - double *max_moist_node; - double *ice_node; - double *alpha; - double *beta; - double *gamma; - layer_data_struct layer[MAX_LAYERS]; - - /************************************************** - Initialize Snow Flux Variables - **************************************************/ - - T2 = energy->T[Nnodes-1]; - Ts_old = energy->T[0]; - T1_old = energy->T[1]; - kappa1 = energy->kappa[0]; - kappa2 = energy->kappa[1]; - Cs1 = energy->Cs[0]; - Cs2 = energy->Cs[1]; - delta_t = (double)(dt)*3600.; - snow_density = snow->density; - snow_depth = snow->depth; - D1 = soil_con.depth[0]; - D2 = soil_con.depth[1]; - max_moist = soil_con.max_moist[0]/(soil_con.depth[0]*1000.); - bubble = soil_con.bubble; - expt = soil_con.expt[0]; - - /*********************************************************** - Prepare Data Sets for Solving Frozen Soil Thermal Fluxes - ***********************************************************/ - if(options.FROZEN_SOIL) { - - 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; - - setup_frozen_soil(soil_con,layer_wet,layer_dry,layer,energy[0], - rec,iveg,Nnodes,mu,kappa_node,Cs_node, - moist_node); - - T_node = energy->T; - dz_node = energy->dz; - ice_node = energy->ice; - - } - - /************************************************** - Find Surface Temperature using Root Brent Method - **************************************************/ - Twidth = 1; - -#if QUICK_FS - surf_temp = root_brent(energy->T[0]-SURF_DT, 0., func_snow_ground_flux, - T2, Ts_old, T1_old, kappa1, kappa2, Cs1, Cs2, - delta_t, snow_density, snow_depth, Tsnow_surf, - D1, D2, dp,moist, ice0, max_moist, bubble, expt, - grnd_flux, deltaH, - snow_flux, T1, - T_node,Tnew_node,dz_node,kappa_node,Cs_node, - moist_node,expt_node,max_moist_node,ice_node, - alpha,beta,gamma,soil_con.ufwc_table, - Nnodes,FIRST_SOLN); -#else - surf_temp = root_brent(energy->T[0]-SURF_DT, 0., func_snow_ground_flux, - T2, Ts_old, T1_old, kappa1, kappa2, Cs1, Cs2, - delta_t, snow_density, snow_depth, Tsnow_surf, - D1, D2, dp,moist, ice0, max_moist, bubble, expt, - grnd_flux, deltaH, - snow_flux, T1, - T_node,Tnew_node,dz_node,kappa_node,Cs_node, - moist_node,expt_node,max_moist_node,ice_node, - alpha,beta,gamma,Nnodes,FIRST_SOLN); -#endif - - if(surf_temp <= -9998) -#if QUICK_FS - error_calc_snow_ground_flux(surf_temp, T2, Ts_old, T1_old, - kappa1, kappa2, Cs1, Cs2, - delta_t, snow_density, snow_depth, - Tsnow_surf, D1, D2, dp, moist, ice0, - max_moist, bubble, expt, - grnd_flux, deltaH, - snow_flux, T1, - T_node,Tnew_node,dz_node,kappa_node,Cs_node, - moist_node,expt_node,max_moist_node,ice_node, - alpha,beta,gamma,soil_con.ufwc_table, - Nnodes,FIRST_SOLN); -#else - error_calc_snow_ground_flux(surf_temp, T2, Ts_old, T1_old, - kappa1, kappa2, Cs1, Cs2, - delta_t, snow_density, snow_depth, - Tsnow_surf, D1, D2, dp, moist, ice0, - max_moist, bubble, expt, - grnd_flux, deltaH, - snow_flux, T1, - T_node,Tnew_node,dz_node,kappa_node,Cs_node, - moist_node,expt_node,max_moist_node,ice_node, - alpha,beta,gamma,Nnodes,FIRST_SOLN); -#endif - - /************************************************** - Recalculate Energy Fluxes Based on Final Temperature - **************************************************/ -#if QUICK_FS - error = solve_snow_ground_flux(surf_temp, T2, Ts_old, T1_old, - kappa1, kappa2, Cs1, Cs2, - delta_t, snow_density, snow_depth, - Tsnow_surf, D1, D2, dp, moist, ice0, - max_moist, bubble, expt, - grnd_flux, deltaH, - snow_flux, T1, - T_node,Tnew_node,dz_node,kappa_node,Cs_node, - moist_node,expt_node,max_moist_node,ice_node, - alpha,beta,gamma,soil_con.ufwc_table, - Nnodes,FIRST_SOLN); -#else - error = solve_snow_ground_flux(surf_temp, T2, Ts_old, T1_old, - kappa1, kappa2, Cs1, Cs2, - delta_t, snow_density, snow_depth, - Tsnow_surf, D1, D2, dp, moist, ice0, - max_moist, bubble, expt, - grnd_flux, deltaH, - snow_flux, T1, - T_node,Tnew_node,dz_node,kappa_node,Cs_node, - moist_node,expt_node,max_moist_node,ice_node, - alpha,beta,gamma,Nnodes,FIRST_SOLN); -#endif - - energy->error += error; - - /*************************************************** - Recalculate Soil Moisture and Thermal Properties - ***************************************************/ - if(options.FROZEN_SOIL) { - -#if QUICK_FS - finish_frozen_soil_calcs(energy,layer_wet,layer_dry,layer,soil_con, - Nnodes,iveg,mu,Tnew_node,dz_node,kappa_node, - Cs_node,moist_node,expt_node,max_moist_node, - soil_con.ufwc_table); -#else - finish_frozen_soil_calcs(energy,layer_wet,layer_dry,layer,soil_con, - Nnodes,iveg,mu,Tnew_node,dz_node,kappa_node, - Cs_node,moist_node,expt_node,max_moist_node); -#endif - - } - - return (surf_temp); - -} - -double solve_snow_ground_flux(double Tsurf, ...) { - - va_list ap; - - double error; - - va_start(ap,Tsurf); - error = func_snow_ground_flux(Tsurf, ap); - va_end(ap); - - return error; - -} - -double error_calc_snow_ground_flux(double Tsurf, ...) { - - va_list ap; - - double error; - - va_start(ap,Tsurf); - error = error_print_snow_ground_flux(Tsurf, ap); - va_end(ap); - - return error; - -} - -double error_print_snow_ground_flux(double Ts, va_list ap) { - - double T2; - double Ts_old; - double T1_old; - double kappa1; - double kappa2; - double Cs1; - double Cs2; - double delta_t; - double snow_density; - double snow_depth; - double surf_temp; - double D1; - double D2; - double dp; - double moist; - double ice0; - double max_moist; - double bubble; - double expt; - double *grnd_flux; - double *deltaH; - double *snow_flux; - double *TMean; - double *T1; - double *T_node; - double *Tnew_node; - double *dz_node; - double *kappa_node; - double *Cs_node; - double *moist_node; - double *expt_node; - double *max_moist_node; - double *ice_node; - double *alpha; - double *beta; - double *gamma; -#if QUICK_FS - double ***ufwc_table; -#endif - int Nnodes; - char *FIRST_SOLN; - - int i; - - /** Initialize Variables **/ - T2 = (double) va_arg(ap, double); - Ts_old = (double) va_arg(ap, double); - T1_old = (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); - delta_t = (double) va_arg(ap, double); - snow_density = (double) va_arg(ap, double); - snow_depth =(double) va_arg(ap, double); - surf_temp = (double) va_arg(ap, double); - D1 = (double) va_arg(ap, double); - D2 = (double) va_arg(ap, double); - dp = (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); - grnd_flux = (double *) va_arg(ap, double *); - deltaH = (double *) va_arg(ap, double *); - snow_flux = (double *) va_arg(ap, double *); - TMean = (double *) va_arg(ap, double *); - T1 = (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 *); - 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 = (double ***) va_arg(ap, double ***); -#endif - Nnodes = (int) va_arg(ap, int); - FIRST_SOLN = (char *) va_arg(ap, char *); - - /* Print Variables */ - fprintf(stderr,"T2 = %lf\n",T2); - fprintf(stderr,"Ts_old = %lf\n",Ts_old); - fprintf(stderr,"T1_old = %lf\n",T1_old); - fprintf(stderr,"kappa1 = %lf\n",kappa1); - fprintf(stderr,"kappa2 = %lf\n",kappa2); - fprintf(stderr,"Cs1 = %lf\n",Cs1); - fprintf(stderr,"Cs2 = %lf\n",Cs2); - fprintf(stderr,"delta_t = %lf\n",delta_t); - fprintf(stderr,"snow_density = %lf\n",snow_density); - fprintf(stderr,"snow_depth = %lf\n",snow_depth); - fprintf(stderr,"surf_temp = %lf\n",surf_temp); - fprintf(stderr,"D1 = %lf\n",D1); - fprintf(stderr,"D2 = %lf\n",D2); - fprintf(stderr,"dp = %lf\n",dp); - fprintf(stderr,"moist = %lf\n",moist); - fprintf(stderr,"ice0 = %lf\n",ice0); - fprintf(stderr,"max_moist = %lf\n",max_moist); - fprintf(stderr,"bubble = %lf\n",bubble); - fprintf(stderr,"expt = %lf\n",expt); - fprintf(stderr,"grnd_flux = %lf\n",grnd_flux[0]); - fprintf(stderr,"deltaH = %lf\n",deltaH[0]); - fprintf(stderr,"snow_flux = %lf\n",snow_flux[0]); - fprintf(stderr,"TMean = %lf\n",TMean[0]); - fprintf(stderr,"T1 = %lf\n",T1[0]); - for(i=0;i -#include -#include - -static char vcid[] = "$Id$"; - -double calc_trans(double deltat, double elevation) -/********************************************************************** - calc_trans Dag Lohmann January 1996 - - This routine computes the transmissity of the atmosphere for the - current time step. - -**********************************************************************/ -{ - double trans, trans_clear; - - trans_clear = (A1_TRANS + A2_TRANS * elevation); - trans = trans_clear * (1 - exp(B_TRANS * pow(deltat, C_TRANS))); - - return trans; -} diff --git a/calc_water_density.c b/calc_water_density.c deleted file mode 100644 index e54a5360f..000000000 --- a/calc_water_density.c +++ /dev/null @@ -1,48 +0,0 @@ -/* - * SUMMARY: WaterDensity.c - Calculate the density of water - * USAGE: Part of DHSVM - * - * AUTHOR: Bart Nijssen - * ORG: University of Washington, Department of Civil Engineering - * E-MAIL: nijssen@u.washington.edu - * ORIG-DATE: Apr-1996 - * LAST-MOD: Tue Jul 6 16:22:32 1999 by VIC Administrator - * DESCRIPTION: Calculate the density of water as a function of temperature - * DESCRIP-END. - * FUNCTIONS: WaterDensity() - * COMMENTS: - */ - -#include -#include -#include - -static char vcid[] = "$Id$"; - -/***************************************************************************** - Function name: WaterDensity() - - Purpose : Calculate the density of water as a function of temperature - using the Thiesen-Scheel-Diesselhorst equation - - Required : - double T - Temperature in C - - Returns : rho - water density in kg/m3 - - Modifies : NA - - Comments : Reference: McCutcheon, S. C., J. L. Martin, and - T. O. Barnwell, Jr, Water quality, - In: Maidment, D. R. (ed.), Handbook of hydrology, - 1993, McGraw-Hill, New York, etc.. (p. 11.6, Fig. 11.1.1) -*****************************************************************************/ -double WaterDensity(double T) -{ - double rho; - - rho = 1000 * (1 - (T + 288.9414)/(508929.2 * (T + 68.12963)) * - (T - 3.9863) * (T - 3.9863)); - - return rho; -} diff --git a/func_snow_ground_flux.c b/func_snow_ground_flux.c deleted file mode 100644 index bc8901df8..000000000 --- a/func_snow_ground_flux.c +++ /dev/null @@ -1,173 +0,0 @@ -#include -#include -#include - -static char vcid[] = "$Id$"; - -double func_snow_ground_flux(double Ts, va_list ap) { -/********************************************************************** - snow_ground_flux Keith Cherkauer January 24, 1997 - - This subroutine computes the surface temeperature under a snow pack - by balancing the ground heat flux, with the flux coming from the - snow pack. - - Reference: - Gel'fan, A. N., "Comparison of Two Methods of Calculating - Soil Freezing Depth," Meteorologiya i Gidrologiya, No. 2, - pp 98-104, 1989. - - UNITS: mks - -**********************************************************************/ - extern option_struct options; - extern debug_struct debug; - - double T2; - double Ts_old; - double T1_old; - double kappa1; /** thermal conductivity of 1st layer */ - double kappa2; /** thermal conductivity of 2nd layer */ - double Cs1; /** volumetric heat capacity of 1st layer **/ - double Cs2; /** volumetric heat capacity of 2nd layer **/ - double delta_t; /** Time Step in Seconds **/ - double snow_density; /** density of snow pack **/ - double snow_depth; /** depth of snow pack **/ - double surf_temp; /** surface temperature of snow pack **/ - double D1; /** thickness of top layer **/ - double D2; /** thickness of second layer **/ - double dp; /** thermal damping depth of soil column **/ - double moist; /** moisture content of top layer **/ - double ice0; /** ice content of top layer **/ - double max_moist; - double bubble; - double expt; - double *grnd_flux; - double *deltaH; - double *snow_flux; - double TMean; - double *T1; - double *T_node; - double *Tnew_node; - double *dz_node; - double *kappa_node; - double *Cs_node; - double *moist_node; - double *expt_node; - double *max_moist_node; - double *ice_node; - double *alpha; - double *beta; - double *gamma; -#if QUICK_FS - double **ufwc_table; -#endif - int Nnodes; - char *FIRST_SOLN; - - double ice; - double kappa_snow; /* thermal conductivity of snow (W/s/K) */ - double error; - - /** Initialize Variables **/ - T2 = (double) va_arg(ap, double); - Ts_old = (double) va_arg(ap, double); - T1_old = (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); - delta_t = (double) va_arg(ap, double); - snow_density = (double) va_arg(ap, double); - snow_depth = (double) va_arg(ap, double); - surf_temp = (double) va_arg(ap, double); - D1 = (double) va_arg(ap, double); - D2 = (double) va_arg(ap, double); - dp = (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); - grnd_flux = (double *) va_arg(ap, double *); - deltaH = (double *) va_arg(ap, double *); - snow_flux = (double *) va_arg(ap, double *); - T1 = (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 *); - 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 = (double **) va_arg(ap, double **); -#endif - Nnodes = (int) va_arg(ap, int); - FIRST_SOLN = (char *) va_arg(ap, char *); - - TMean = Ts; - - kappa_snow = 2.9302e-6 * (snow_density) * (snow_density); - - *snow_flux = kappa_snow * (TMean - surf_temp) / snow_depth; - - if(!options.FROZEN_SOIL) { - /************************************************* - Use Xu's Equations to Calculate Thermal Fluxes - *************************************************/ - *T1 = estimate_T1(TMean, T1_old, T2, D1, D2, kappa1, kappa2, Cs1, - Cs2, dp, delta_t); - - } - else { - /************************************************************* - Explicitly Solve Thermal Fluxes For all Soil Thermal Nodes - *************************************************************/ - T_node[0] = TMean; -#if QUICK_FS - solve_T_profile(Tnew_node,T_node,dz_node,kappa_node,Cs_node, - moist_node,delta_t,max_moist_node, - bubble,expt_node,ice_node,alpha,beta,gamma,ufwc_table, - Nnodes,FIRST_SOLN,FALSE); -#else - solve_T_profile(Tnew_node,T_node,dz_node,kappa_node,Cs_node, - moist_node,delta_t,max_moist_node, - bubble,expt_node,ice_node,alpha,beta,gamma,Nnodes, - FIRST_SOLN,FALSE); -#endif - *T1 = Tnew_node[1]; - } - - /************************************************************ - Compute the Ground Heat Flux Through the Upper Soil Layer - ************************************************************/ - *grnd_flux = kappa1/D1*((*T1) - (TMean)); - - /************************************************************** - Compute the Change in Heat Storage in the Upper Soil Layer - **************************************************************/ - if(options.FROZEN_SOIL && (TMean+ *T1)/2.<0.) { - ice = moist - maximum_unfrozen_water((TMean+ *T1)/2.,max_moist, - bubble,expt); - if(ice<0.) ice=0.; - if(ice>max_moist) ice=max_moist; - } - else ice=0.; - - *deltaH = Cs1 * ((Ts_old + T1_old)/2. - (TMean + *T1)/2.) * D1 / delta_t; - *deltaH -= ice_density*Lf*(ice0-ice)*D1/delta_t; - - /******************************* - Compute Energy Balance Error - *******************************/ - error = *deltaH + *grnd_flux - *snow_flux; - - return error; - -} diff --git a/initialize_energy.c b/initialize_energy.c deleted file mode 100644 index 76264e399..000000000 --- a/initialize_energy.c +++ /dev/null @@ -1,561 +0,0 @@ -#include -#include -#include -#include - -static char vcid[] = "$Id$"; - -void initialize_energy_bal (energy_bal_struct **energy, - cell_data_struct ***cell, - soil_con_struct *soil_con, - double surf_temp, - double *mu, - int cellnum, - int veg_num, - int Nnodes, - int Ndist, - FILE *fsoil) -/********************************************************************** - initialize_energy_bal Keith Cherkauer July 30, 1996 - - This routine initializes the energy balance data structure for use - with the frozen soils model. - - Soil temperatures are initialized using a linear interpolation - between the surface temperature (assumed = air temperature) and - the average annual air temperature (assumed that bottom of 2nd soil - layer is at const T = average annual air temp). - - UNITS: (m, s, kg, C, moisture in mm) unless otherwise specified - - variable modified: - energy[veg].dz - energy[veg].T - energy[veg].fdepth - cell[veg].layer[index].T - cell[veg].layer[index].T_thaw - cell[veg].layer[index].T_froz - cell[veg].layer[index].moist - cell[veg].layer[index].moist_thaw - cell[veg].layer[index].moist_froz - cell[veg].layer[index].ice - cell[veg].layer[index].tdepth - cell[veg].layer[index].fdepth - cell[veg].layer[index].kappa - cell[veg].layer[index].Cs - - Modifications: - 10-15-96 modified to reflect fixes in the frozen soils code KAC - 07-10-97 modified to read water content in (m/m) instead of (mm) KAC - 09-18-97 modified to initialize with no files, and to initialize - FULL_ENERGY case without FROZEN_SOIL KAC - 06-24-98 modified to use distributed soil moisture only KAC - -**********************************************************************/ -{ - extern option_struct options; - extern debug_struct debug; -#if QUICK_FS - extern double temps[]; -#endif - - char tmpstr[MAXSTRING]; - char ErrStr[MAXSTRING]; - int i, j, ii, veg, index; - int tmpint; - int dry; - int band; - int zindex; - double sum, Lsum, Zsum, dp, Ltotal; - double tmpdp, tmpadj; - double *kappa, *Cs, *M; - double moist[MAXlayer], ice[MAXlayer]; - double unfrozen, frozen; - double **layer_ice; - double **layer_tmp; - double *EMPTY; -#if QUICK_FS - double Aufwc, Bufwc; -#endif - char *EMPTY_C; - - dp = soil_con[0].dp; - Ltotal = 0; - for(index=0;index (Lsum - soil_con[0].depth[index])) { - if(energy[veg][band].fdepth[1] < Lsum) { - cell[dry][veg][band].layer[index].tdepth - = energy[veg][band].fdepth[1] - - (Lsum - soil_con[0].depth[index]); - cell[dry][veg][band].layer[index].moist_thaw = moist[index] - * soil_con[0].depth[index] * 1000.; - if(energy[veg][band].fdepth[0] < Lsum) { - cell[dry][veg][band].layer[index].fdepth - = energy[veg][band].fdepth[0] - - (Lsum - soil_con[0].depth[index]); - cell[dry][veg][band].layer[index].moist_froz - = soil_con[0].depth[index] - * moist[index] * 1000.; - cell[dry][veg][band].layer[index].ice - = soil_con[0].depth[index] - * ice[index] * 1000.; - cell[dry][veg][band].layer[index].moist - = soil_con[0].depth[index] - * moist[index] * 1000.; - } - else { - cell[dry][veg][band].layer[index].fdepth - = soil_con[0].depth[index]; - cell[dry][veg][band].layer[index].moist_froz - = soil_con[0].depth[index] - * moist[index] * 1000.; - cell[dry][veg][band].layer[index].ice - = soil_con[0].depth[index] - * ice[index] * 1000.; - cell[dry][veg][band].layer[index].moist = 0.; - } - } - else { - cell[dry][veg][band].layer[index].tdepth - = soil_con[0].depth[index]; - cell[dry][veg][band].layer[index].fdepth - = soil_con[0].depth[index]; - cell[dry][veg][band].layer[index].moist_thaw - = soil_con[0].depth[index] * moist[index] * 1000.; - cell[dry][veg][band].layer[index].moist_froz = 0.; - cell[dry][veg][band].layer[index].ice = 0.; - cell[dry][veg][band].layer[index].moist = 0.; - } - } - else if(energy[veg][band].fdepth[0] - > (Lsum - soil_con[0].depth[index])) { - if(energy[veg][band].fdepth[0] < Lsum) { - cell[dry][veg][band].layer[index].tdepth = 0.; - cell[dry][veg][band].layer[index].fdepth - = energy[veg][band].fdepth[0] - - (Lsum - soil_con[0].depth[index]); - cell[dry][veg][band].layer[index].moist_thaw = 0.; - cell[dry][veg][band].layer[index].moist_froz - = soil_con[0].depth[index] * moist[index] * 1000.; - cell[dry][veg][band].layer[index].ice - = soil_con[0].depth[index] - * ice[index] * 1000.; - cell[dry][veg][band].layer[index].moist - = soil_con[0].depth[index] - * moist[index] * 1000.; - } - else { - cell[dry][veg][band].layer[index].tdepth = 0.; - cell[dry][veg][band].layer[index].fdepth - = soil_con[0].depth[index]; - cell[dry][veg][band].layer[index].moist_thaw = 0.; - cell[dry][veg][band].layer[index].moist_froz - = soil_con[0].depth[index] * moist[index] * 1000.; - cell[dry][veg][band].layer[index].ice - = soil_con[0].depth[index] - * ice[index] * 1000.; - cell[dry][veg][band].layer[index].moist = 0.; - } - } - else { - cell[dry][veg][band].layer[index].tdepth = 0.; - cell[dry][veg][band].layer[index].fdepth = 0.; - cell[dry][veg][band].layer[index].moist_thaw = 0.; - cell[dry][veg][band].layer[index].moist_froz = 0.; - cell[dry][veg][band].layer[index].ice = 0.; - cell[dry][veg][band].layer[index].moist - = soil_con[0].depth[index] - * moist[index] * 1000.; - } - } - layer_ice = (double **)malloc(options.Nlayer*sizeof(double *)); - for(i=0;i 0.) - cell[dry][veg][band].layer[index].moist_thaw - = cell[dry][veg][band].layer[index].moist; - if(cell[dry][veg][band].layer[index].tdepth - == soil_con[0].depth[index]) - cell[dry][veg][band].layer[index].moist = 0.; - else { - if(cell[dry][veg][band].layer[index].fdepth > 0.) - cell[dry][veg][band].layer[index].moist_froz - = cell[dry][veg][band].layer[index].moist; - if(cell[dry][veg][band].layer[index].fdepth - == soil_con[0].depth[index]) - cell[dry][veg][band].layer[index].moist = 0.; - } - if(cell[dry][veg][band].layer[index].T_froz<0. - && cell[dry][veg][band].layer[index].T_froz != -999.) { - unfrozen - = maximum_unfrozen_water(cell[dry][veg][band].layer[index].T_froz, - soil_con[0].max_moist[index], - soil_con[0].bubble, - soil_con[0].expt[index]); - if(unfrozen>soil_con[0].max_moist[index] || unfrozen<0.) - unfrozen = soil_con[0].max_moist[index]; - cell[dry][veg][band].layer[index].unfrozen = unfrozen; - - frozen = cell[dry][veg][band].layer[index].moist_froz - unfrozen; - if(frozen < 0.0) { - frozen = 0.0; - unfrozen = cell[dry][veg][band].layer[index].moist_froz; - } - cell[dry][veg][band].layer[index].ice = frozen; - cell[dry][veg][band].layer[index].moist_froz = unfrozen; - } - else if(cell[dry][veg][band].layer[index].T_froz == 0.) { - cell[dry][veg][band].layer[index].unfrozen - = soil_con[0].max_moist[index]; - cell[dry][veg][band].layer[index].ice = 0.; - } - else if(cell[dry][veg][band].layer[index].T_froz != -999.) - nrerror("ERROR: Frozen Layer Temperature > 0C"); - - } - - layer_ice = (double **)malloc(options.Nlayer*sizeof(double *)); - for(i=0;i0) { - index = 0; - Zsum = energy[veg][band].dz[index]/2.; - while(index=Nnodes - || (int)((Zsum+energy[veg][band].dz[index+1]/2.)*1000.+0.5) - > (int)((soil_con[0].depth[0])*1000.+0.5)) { - nrerror("One soil temperature depth must equal the depth of the top soil layer"); - } - energy[veg][band].T1_index = index+1; - - for(dry=0;dry0.) energy[veg][band].frozen=TRUE; - } - } /** End loop through elevation bands **/ - } /** End loop through vegetation types **/ - - if(options.FROZEN_SOIL) { - - /*********************************************************** - Prepare soil constants for use in thermal flux solutions - ***********************************************************/ - - for(zindex=0;zindex -#include -#include - -static char vcid[] = "$Id$"; - -#define DAYLENGTH 0.0 /* daylength coefficient: selects a */ - /* definition for sunrise/sunset : */ - /* 0.0 : center of sun even with horizon */ - /* 0.26667 : top of sun even w/ horizon */ - /* 0.8333 : top of sun apparently even */ - -/******************************************************************************* - Function to calculate the daily incoming shortwave in W/m2 (Shuttleworth 1993) - - Modifications: - 08-04-98 Code was modified to use Reiner's corrections so that it would - function correctly in high latitudes during global simulations KAC - - References: - Forsythe, W.C., et al., A Model Comparison for Daylength as a Function of - Latitude and Day of Year, Ecological Modelling 80(1995), pp. 87-95. - - -*******************************************************************************/ -double in_shortwave(float lat, int day, double trans) -{ - double omega_s; /* sunset hour angle in radians */ - double declination; /* solar declination in radians */ - double extra_solar; /* extraterrestrial solar radiation (W/m^2) */ - double dist; /* relative distance bewteen earth and sun */ - double lat_rad; /* latitude in radians */ - - lat_rad = lat * DtoR; - - /** equations established for global simulations (see Forsythe, 1995) **/ - declination = asin(0.39795 * cos(0.2163108 + 2. * - atan(0.9671396 * tan(0.00860*(day-186))))); - - omega_s = (sin(DAYLENGTH*PI/180.)+sin(lat_rad)*sin(declination)) / - (cos(lat_rad)*cos(declination)); - if(omega_s < -1.) - omega_s = -1.; - else if(omega_s > 1.) - omega_s = 1.; - omega_s = PI - acos(omega_s); - - dist = 1.0 + 0.033 * cos(2.0 * PI/DAYS_PER_YEAR * day); - - extra_solar = SOLAR_CONSTANT / PI * dist * - (omega_s * sin(declination) * sin(lat_rad) + - cos(declination) * cos(lat_rad) * sin(omega_s)); - - return (trans*extra_solar); -} - -/******************************************************************************* -* - Function to calculate the daily net outgoing longwave radiation in W/m2 - (Bras 1990) -*******************************************************************************/ -double net_out_longwave(double trans, - double trans_clear, - double tair, - double vapor, - double *cloudiness) -{ - double emissivity; /* emissivity of the atmosphere */ - double cloudfactor; /* cloudiness correction factor */ - double longwave; /* net ourgoing longwave */ - double t_kelvin; /* temperature in Kelvin */ - double v_mbar; /* vapor pressure in mbar */ - - t_kelvin = tair + 273.15; - v_mbar = vapor/100.0; - - *cloudiness = 1.0/0.65 * sqrt(1.0 - trans/trans_clear); - *cloudiness = (*cloudiness > 1.0) ? 1.0 : *cloudiness; - emissivity = 0.70 + 0.0000595 * v_mbar * exp(1500.0/t_kelvin); - cloudfactor = 1.0 + 0.17 * (*cloudiness) * (*cloudiness); - - longwave = (1.0 - emissivity * cloudfactor) * STEFAN * (t_kelvin) - * (t_kelvin) * (t_kelvin) * (t_kelvin); - - return longwave; -} - -#undef DAYLENGTH diff --git a/make_out_data.c b/make_out_data.c deleted file mode 100644 index 30d3227c2..000000000 --- a/make_out_data.c +++ /dev/null @@ -1,18 +0,0 @@ -#include -#include -#include - -out_data_struct *make_out_data(int nrecs) -/********************************************************************** - make_out_data Dag Lohmann January 1996 - - This routine creates an array of out data structures, one for each - time step. - -**********************************************************************/ -{ - out_data_struct *temp; - - temp = (out_data_struct*) calloc(nrecs, sizeof(out_data_struct)); - return temp; -} diff --git a/makefile b/makefile deleted file mode 100644 index 4c78dcae1..000000000 --- a/makefile +++ /dev/null @@ -1,238 +0,0 @@ -SHELL = /bin/csh - -CC = gcc -I. -O -LIBRARY = -lm -UTILDIR = ./ -DISTDIR = ./ -FULLDIR = ./ -FROZDIR = ./ - -INCFILE = vicNl.h vicNl_def.h global.h snow.h user_def.h rad_and_vpd.h - -OBJS = vicNl.o cmd_proc.o check_files.o read_soilparam.o open_file.o \ - make_in_and_outfiles.o read_atmosdata.o read_vegparam.o make_dmy.o \ - get_global_param.o initialize_soil.o close_files.o \ - calc_trans.o calc_netshort.o aurad.o priestley.o svp.o \ - long_shortwave.o nrerror.o vicerror.o rad_and_vpd.o \ - make_veg_var.o make_cell_data.o initialize_veg.o penman.o arno_evap.o \ - put_data.o polint.o initialize_global.o read_snowband.o \ - read_snowmodel.o read_sawd.o read_sawd_binary.o \ - initialize_atmos.o correct_precip.o compress_files.o \ - write_soilparam.o write_vegparam.o write_atmosdata.o write_data.o \ - write_vegvar.o write_layer.o make_dist_prcp.o free_dist_prcp.o \ - initialize_energy.o frozen_soil.o \ - soil_conduction.o make_energy_bal.o canopy_evap.o runoff.o \ - modify_Ksat.o full_energy.o func_surf_energy_bal.o read_rosemount.o \ - initialize_snow.o snow_melt.o root_brent.o calc_water_density.o \ - SnowPackEnergyBalance.o make_snow_data.o StabilityCorrection.o \ - snow_utility.o func_snow_ground_flux.o soil_thermal_eqn.o \ - calc_surf_energy_bal.o write_debug.o dist_prec.o open_debug.o \ - initialize_new_storm.o redistribute_during_storm.o \ - snow_intercept.o massrelease.o CalcAerodynamic.o \ - calc_veg_params.o read_veglib.o calc_rainonly.o \ - calc_long_shortwave.o calc_air_temperature.o read_PILPS2c.o \ - calc_snow_ground_flux.o read_forcing_data.o read_soilparam_arc.o \ - read_arcinfo_ascii.o solve_snow.o calc_root_fraction.o free_vegcon.o - - -all: - make clean - make model - -default: - make model - -clean: - /bin/rm -f *.o core log - -model: $(OBJS) - $(CC) -o vicNl $(OBJS) $(LIBRARY) - -vicNl.o: vicNl.c $(INCFILE) - $(CC) -c vicNl.c -cmd_proc.o: cmd_proc.c $(INCFILE) - $(CC) -c cmd_proc.c -runoff.o: runoff.c $(INCFILE) - $(CC) -c runoff.c -initialize_global.o: initialize_global.c $(INCFILE) - $(CC) -c initialize_global.c - -# DATA HANDLING CODE -make_veg_var.o: make_veg_var.c $(INCFILE) - $(CC) -c make_veg_var.c -make_cell_data.o: make_cell_data.c $(INCFILE) - $(CC) -c make_cell_data.c -make_dmy.o: make_dmy.c $(INCFILE) - $(CC) -c make_dmy.c -initialize_soil.o: initialize_soil.c $(INCFILE) - $(CC) -c initialize_soil.c -initialize_veg.o: initialize_veg.c $(INCFILE) - $(CC) -c initialize_veg.c -initialize_atmos.o: initialize_atmos.c $(INCFILE) - $(CC) -c initialize_atmos.c -write_vegvar.o: write_vegvar.c $(INCFILE) - $(CC) -c write_vegvar.c -write_layer.o: write_layer.c $(INCFILE) - $(CC) -c write_layer.c - -# INPUT/OUTPUT CODE -read_forcing_data.o: read_forcing_data.c $(INCFILE) - $(CC) -c read_forcing_data.c -read_soilparam.o: read_soilparam.c $(INCFILE) - $(CC) -c read_soilparam.c -read_soilparam_arc.o: read_soilparam_arc.c $(INCFILE) - $(CC) -c read_soilparam_arc.c -read_arcinfo_ascii.o: read_arcinfo_ascii.c $(INCFILE) - $(CC) -c read_arcinfo_ascii.c -write_soilparam.o: write_soilparam.c $(INCFILE) - $(CC) -c write_soilparam.c -open_file.o: open_file.c $(INCFILE) - $(CC) -c open_file.c -open_debug.o: open_debug.c $(INCFILE) - $(CC) -c open_debug.c -read_vegparam.o: read_vegparam.c $(INCFILE) - $(CC) -c read_vegparam.c -read_veglib.o: read_veglib.c $(INCFILE) - $(CC) -c read_veglib.c -write_vegparam.o: write_vegparam.c $(INCFILE) - $(CC) -c write_vegparam.c -make_in_and_outfiles.o: make_in_and_outfiles.c $(INCFILE) - $(CC) -c make_in_and_outfiles.c -read_atmosdata.o: read_atmosdata.c $(INCFILE) - $(CC) -c read_atmosdata.c -read_snowmodel.o: read_snowmodel.c $(INCFILE) - $(CC) -c read_snowmodel.c -read_snowband.o: read_snowband.c $(INCFILE) - $(CC) -c read_snowband.c -read_sawd.o: read_sawd.c $(INCFILE) - $(CC) -c read_sawd.c -read_sawd_binary.o: read_sawd_binary.c $(INCFILE) - $(CC) -c read_sawd_binary.c -read_PILPS2c.o: read_PILPS2c.c $(INCFILE) - $(CC) -c read_PILPS2c.c -write_atmosdata.o: write_atmosdata.c $(INCFILE) - $(CC) -c write_atmosdata.c -get_global_param.o: get_global_param.c $(INCFILE) - $(CC) -c get_global_param.c -close_files.o: close_files.c $(INCFILE) - $(CC) -c close_files.c -put_data.o: put_data.c $(INCFILE) - $(CC) -c put_data.c -write_data.o: write_data.c $(INCFILE) - $(CC) -c write_data.c -write_debug.o: write_debug.c $(INCFILE) - $(CC) -c write_debug.c -check_files.o: check_files.c $(INCFILE) - $(CC) -c check_files.c -calc_root_fraction.o: calc_root_fraction.c $(INCFILE) - $(CC) -c calc_root_fraction.c -free_vegcon.o: free_vegcon.c $(INCFILE) - $(CC) -c free_vegcon.c - -# RADIATION AND EVAPORATION CODE -rad_and_vpd.o: rad_and_vpd.c $(INCFILE) - $(CC) -c rad_and_vpd.c -calc_trans.o: calc_trans.c $(INCFILE) - $(CC) -c calc_trans.c -calc_netshort.o: calc_netshort.c $(INCFILE) - $(CC) -c calc_netshort.c -aurad.o: aurad.c $(INCFILE) - $(CC) -c aurad.c -priestley.o: priestley.c $(INCFILE) - $(CC) -c priestley.c -svp.o: svp.c $(INCFILE) - $(CC) -c svp.c -long_shortwave.o: long_shortwave.c $(INCFILE) - $(CC) -c long_shortwave.c -canopy_evap.o: canopy_evap.c $(INCFILE) - $(CC) -c canopy_evap.c -penman.o: penman.c $(INCFILE) - $(CC) -c penman.c -arno_evap.o: arno_evap.c $(INCFILE) - $(CC) -c arno_evap.c - -# DITRIBUTED PRECIPITATION CODE -polint.o: polint.c $(INCFILE) - $(CC) -c polint.c -nrerror.o: nrerror.c $(INCFILE) - $(CC) -c nrerror.c -vicerror.o: vicerror.c $(INCFILE) - $(CC) -c vicerror.c -make_prcp_var.o: make_prcp_var.c $(INCFILE) - $(CC) -c make_prcp_var.c -make_dist_prcp.o: make_dist_prcp.c $(INCFILE) - $(CC) -c make_dist_prcp.c -free_dist_prcp.o: free_dist_prcp.c $(INCFILE) - $(CC) -c free_dist_prcp.c -initialize_new_storm.o: initialize_new_storm.c $(INCFILE) - $(CC) -c initialize_new_storm.c -redistribute_during_storm.o: redistribute_during_storm.c $(INCFILE) - $(CC) -c redistribute_during_storm.c - -# SNOW MELT CODE -make_snow_data.o: make_snow_data.c $(INCFILE) - $(CC) -c make_snow_data.c -initialize_snow.o: initialize_snow.c $(INCFILE) - $(CC) -c initialize_snow.c -solve_snow.o: solve_snow.c $(INCFILE) - $(CC) -c solve_snow.c -snow_melt.o: snow_melt.c $(INCFILE) - $(CC) -c snow_melt.c -root_brent.o: root_brent.c $(INCFILE) - $(CC) -c root_brent.c -calc_water_density.o: calc_water_density.c $(INCFILE) - $(CC) -c calc_water_density.c -SnowPackEnergyBalance.o: SnowPackEnergyBalance.c $(INCFILE) - $(CC) -c SnowPackEnergyBalance.c -StabilityCorrection.o: StabilityCorrection.c $(INCFILE) - $(CC) -c StabilityCorrection.c -func_snow_ground_flux.o: func_snow_ground_flux.c $(INCFILE) - $(CC) -c func_snow_ground_flux.c -calc_snow_ground_flux.o: calc_snow_ground_flux.c $(INCFILE) - $(CC) -c calc_snow_ground_flux.c -snow_utility.o: snow_utility.c $(INCFILE) - $(CC) -c snow_utility.c -soil_thermal_eqn.o: soil_thermal_eqn.c $(INCFILE) - $(CC) -c soil_thermal_eqn.c -snow_intercept.o: snow_intercept.c $(INCFILE) - $(CC) -c snow_intercept.c -massrelease.o: massrelease.c $(INCFILE) - $(CC) -c massrelease.c - -# FROZEN SOIL CODE -frozen_soil.o: frozen_soil.c $(INCFILE) - $(CC) -c frozen_soil.c -initialize_energy.o: initialize_energy.c $(INCFILE) - $(CC) -c initialize_energy.c -soil_conduction.o: soil_conduction.c $(INCFILE) - $(CC) -c soil_conduction.c -make_energy_bal.o: make_energy_bal.c $(INCFILE) - $(CC) -c make_energy_bal.c -modify_Ksat.o: modify_Ksat.c $(INCFILE) - $(CC) -c modify_Ksat.c - -# FULL ENERGY BALANCE CODE -dist_prec.o: dist_prec.c $(INCFILE) - $(CC) -c dist_prec.c -full_energy.o: full_energy.c $(INCFILE) - $(CC) -c full_energy.c -func_surf_energy_bal.o: func_surf_energy_bal.c $(INCFILE) - $(CC) -c func_surf_energy_bal.c -calc_surf_energy_bal.o: calc_surf_energy_bal.c $(INCFILE) - $(CC) -c calc_surf_energy_bal.c -read_rosemount.o: read_rosemount.c $(INCFILE) - $(CC) -c read_rosemount.c -correct_precip.o: correct_precip.c $(INCFILE) - $(CC) -c correct_precip.c -compress_files.o: compress_files.c $(INCFILE) - $(CC) -c compress_files.c -CalcAerodynamic.o: CalcAerodynamic.c $(INCFILE) - $(CC) -c CalcAerodynamic.c -calc_veg_params.o: calc_veg_params.c $(INCFILE) - $(CC) -c calc_veg_params.c -calc_rainonly.o: calc_rainonly.c $(INCFILE) - $(CC) -c calc_rainonly.c -calc_air_temperature.o: calc_air_temperature.c $(INCFILE) - $(CC) -c calc_air_temperature.c -calc_long_shortwave.o: calc_long_shortwave.c $(INCFILE) - $(CC) -c calc_long_shortwave.c diff --git a/polint.c b/polint.c deleted file mode 100644 index a1d31db36..000000000 --- a/polint.c +++ /dev/null @@ -1,38 +0,0 @@ -#include -#include -#include -#include - -void polint(double *xa,double *ya,int n,double x,double *y,double *dy) -{ - int i,m,ns=0; - double den,dif,dift,ho,hp,w; - double *c,*d; - - dif=fabs(x-xa[0]); - c=(double *)malloc(n*sizeof(double)); - d=(double *)malloc(n*sizeof(double)); - for (i=0;i -#include -#include - -static char vcid[] = "$Id$"; - -double priestley(double tair, - double rad) -/********************************************************************** - priestley Dag Lohmann January 1996 - - This routine computes potential evaporation using the Priestly-Taylor - method. - - REFERENCE: Kimball, J. S., et. al., "An Improved method for estimating - surface humidity from daily minimum temperature", - Agriculture and Froest Meteorology, 85(1997) 87-98. - - double tair C Current air temperature - double rad kJ/day/m^2 Net radiation at the surface - - returns potential evapotranspiration in kg-m^2/s - -**********************************************************************/ -{ - - double slope, pot; - - slope = svp_slope(tair); /* gradient in kPa/C */ - - pot = ALPHA_PT * slope/(slope + GAMMA_PT) * 0.9 * rad / LV_PT; - /** The 0.9 is because ground heat flux (G) is assumed to be 10% of the - net radiation (Rn), so the radiation term (Rn - G) can be - approximated as (0.9 * Rn) **/ - - return pot; - -} diff --git a/rad_and_vpd.c b/rad_and_vpd.c deleted file mode 100644 index b0731c985..000000000 --- a/rad_and_vpd.c +++ /dev/null @@ -1,50 +0,0 @@ -#include -#include -#include - -static char vcid[] = "$Id$"; - -double estimate_dew_point(double annual_prec, - double tmin, - double tmax, - double priest, - double day_len, - int day) { -/******************************************************************* - estimate_dew_point Keith Cherkauer November 11, 1998 - - This routine was written to replace the pre-publication version - of dew point temperature estimation routine from Kimball et. al., - with the version that was published. The set-up code for this - routine which was orginally in rad_and_vpd.c has been relocated - into the initialize_atmos.c routines which now handle initializing - all of the atmospheric variables. - - As described in Kimball, J. S., S. W. Running, and R. Nemani, "An - improved method for estimating surface humidity from daily minimum - temperature", Agricultural and Forest Engineering, 85(1997) 87-98. - - Variables: - dmy_struct *dmy N/A Structure containing date information - double *yearly_prec mm/year Precipitation per year - double tmin C Minimum daily temperature - double tmax C Maximum daily temperature - double priest kg-m^2/s Priestly potential evaporation for the day - double day_len s Length of the current day in seconds - int day N/A Day in simulation - - Returns Dew Point Temperature in Degrees C. - -*******************************************************************/ - - double Tdew; - double EF; - - EF = ((priest) * day_len) / (annual_prec); - Tdew = (tmin + KELVIN) * ( -0.127 + 1.121 * (1.003 - 1.444 * EF - + 12.312 * EF * EF - -32.766 * EF * EF * EF) - + 0.0006 * (tmax - tmin) ); - - return(Tdew-KELVIN); -} diff --git a/rad_and_vpd.h b/rad_and_vpd.h deleted file mode 100644 index c6da03963..000000000 --- a/rad_and_vpd.h +++ /dev/null @@ -1,61 +0,0 @@ -/* - filename : rad_and_vpd.h - purpose : include file for rad_and_vpd.h, part of the water balance version - of the VIC-2L model. - programmer: Bart Nijssen - date : August 3, 1995 - changes : -*/ - -/* define statements */ - -/* miscellaneous constants */ -#ifndef MAXSTRING -#define MAXSTRING 512 -#endif - -#ifndef MINSTRING -#define MINSTRING 20 -#endif - -#define DAYS_PER_YEAR 365. -#define MAX_YEARS 100 -#define DtoR 0.017453293 /* degrees to radians */ -#define PI 3.1415927 -#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 */ - -/* constants for daily transmittance, values from Kimball et al., 1995 */ -#define A1_TRANS 0.60 /* maximum clear sky transmittance at MSL */ -#define A2_TRANS 0.0000295 /* lapse rate per meter of maximum clear sky - transmittance */ -#define B_TRANS -0.0030 /* empirical constant */ -#define C_TRANS 2.4 /* empirical constant */ - -/* constants for Priestley-Taylor potential evaporation, values from Kimball - et al., 1995 */ -#define ALPHA_PT 1.26 /* Priestley-Taylor parameter */ -#define GAMMA_PT 0.066 /* psychrometric constant (kPa/K) */ -#define LV_PT 2.5E6 /* latent heat of vaporization (J/kg) */ - -/* 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 */ - -/* define constants for regression equation to correct vapor pressure estimate - based on minimum temperature */ - -#define HUMID_RATIO 2.25 /* ratio of total yearly potential evaporation - to total yearly precipitation*/ - - - - diff --git a/read_PILPS2c.c b/read_PILPS2c.c deleted file mode 100644 index 0ece4241a..000000000 --- a/read_PILPS2c.c +++ /dev/null @@ -1,99 +0,0 @@ -#include -#include -#include - -static char vcid[] = "$Id$"; - -void read_PILPS2c(atmos_data_struct *temp, - FILE *PILPS2c, - int *nrecs, - int dt, - int file_dt, - int fileskip) -/********************************************************************** - read_PILPS2c Dag Lohmann Feb. 12, 1998 - - This routine reads in atmospheric data values from the PILPS2c - data files. - Input Output - Units Units - shortwave W/m2 W/m2 - longwave W/m2 W/m2 - prec - mm mm - air_temp - C C - wind - m/s m/s - pressure - mbar kPa - spec_humid - kg/kg kg/kg - -**********************************************************************/ -{ - extern debug_struct debug; - extern param_set_struct param_set; - - int i, n, rec, maxline = 210; - int year, month, day, hour; - int store_rec; - int skip_bytes; - char str[210]; - - /** locate starting record **/ - skip_bytes = (int)((float)(dt * fileskip)) / (float)file_dt - 1; - if((dt * fileskip) % (24 / file_dt) > 0) - nrerror("Currently unable to handle a model starting date that does not correspond to a line in the forcing file."); - for(i=0;i dt) { - /** Time step used by model finer than that used in forcing file: - Repeat Data Into Extra Columns **/ - store_rec = rec; - for(i=1;i
-#include -#include -#include - -static char vcid[] = "$Id$"; - -void read_atmosdata(atmos_data_struct *temp, - FILE *atmosf, - int *nrecs, - int dt, - int file_dt, - int fileskip) -/********************************************************************** - read_atmosdata Dag Lohmann - - This routine reads in atmospheric data values from gridded daily - precipitation station data records. - - Modifications: - 7-25-96 Routine modified to read data into structure created - for time steps of less than 1 day. KAC - -**********************************************************************/ -{ - extern param_set_struct param_set; - - int i, n, rec, maxline = 210; - int store_rec; - int skip_bytes; - char str[210]; - double junk; - - /** locate starting record **/ - skip_bytes = (int)((float)(dt * fileskip)) / (float)file_dt - 1; - if((dt * fileskip) % (24 / file_dt) > 0) - nrerror("Currently unable to handle a model starting date that does not correspond to a line in the forcing file."); - for(i=0;i0) { - if(temp[rec-1].tmaxtemp[rec].tmax) temp[rec].tmax=temp[rec-1].tmin; - } - rec++; - - if(file_dt < dt) { - /** Time Step in Forcing File Finer than Used by Model: - Skip Records **/ - for(i=0;i
dt) { - /** Time step used by model finer than that used in forcing file: - Repeat Data Into Extra Columns **/ - store_rec = rec - 1; - for(i=1;i -#include -#include -#include - -static char vcid[] = "$Id$"; - -snow_data_struct read_initial_snow(FILE *initsnow, - int cellnum, - int veg, - int band) -/****************************************************************************** - read_initial_snow Keith Cherkauer January 11, 1999 - - This routine reads parameters for an initial snowpack - if present. - - Parameters Read from File: - TYPE NAME UNITS DESCRIPTION - int cellnum N/A cell number for current data - char snow N/A 1 = snow is present, 0 = no snow present - int last_snow (dt) number of time steps since last snowfall - double swq mm snow water equivalent - double surf_temp C surface temperature of snowpack - double density kg.m^3 density of snowpack - - Parameters Computed from those in the File: - TYPE NAME UNITS DESCRIPTION - double coverage fract fraction of grid cell covered in snow - - Modifications: - -**********************************************************************/ -{ - extern option_struct options; - extern debug_struct debug; - - char ErrStr[MAXSTRING]; - char tmpstr[MAXSTRING]; - int index, tmpveg, tmpband; - snow_data_struct temp; - - rewind(initsnow); - - fscanf(initsnow, "%s", tmpstr); - if(tmpstr[0] != '#') { - index = atoi(tmpstr); - fscanf(initsnow, "%i %i", &tmpveg, &tmpband); - } - else index = tmpveg = tmpband = -999; - while(index != cellnum && veg != tmpveg && band != tmpband) { - fgets(tmpstr,MAXSTRING,initsnow); - fscanf(initsnow, "%s", tmpstr); - if(tmpstr[0] != '#') { - index = atoi(tmpstr); - fscanf(initsnow, "%i %i", &tmpveg, &tmpband); - } - else index = tmpveg = tmpband = -999; - } - - if(!feof(initsnow)) { - fscanf(initsnow, "%s", &tmpstr); - if(tmpstr[0]=='1') temp.snow = TRUE; - else temp.snow = FALSE; - fscanf(initsnow, "%i", &temp.last_snow); - fscanf(initsnow, "%lf", &temp.swq); - fscanf(initsnow, "%lf", &temp.surf_temp); - fscanf(initsnow, "%lf", &temp.density); - - if(temp.swq < MAX_FULL_COVERAGE_SWQ) - temp.coverage = 1. / MAX_FULL_COVERAGE_SWQ * temp.swq; - else if(temp.swq > 0.) temp.coverage = 1.0; - else temp.coverage = 0.0; - } - else { - /** Assume no snow cover **/ - fprintf(stderr,"No snow cover defined for cell number %i.\nModel assuming no snow cover.\n",cellnum); - temp.snow = 0; - temp.last_snow = 0; - temp.swq = 0.0; - temp.surf_temp = 0.0; - temp.density = 0.0; - temp.coverage = 0.0; - } - - return temp; - -} diff --git a/read_initial_soil_thermal.c b/read_initial_soil_thermal.c deleted file mode 100644 index 1a24d8aec..000000000 --- a/read_initial_soil_thermal.c +++ /dev/null @@ -1,162 +0,0 @@ -#include -#include -#include -#include - -static char vcid[] = "$Id$"; - -void read_initial_soil_thermal(FILE *initsoil, - int cellnum, - int veg, - int band, - int Nnodes, - double dp, - double *fdepth, - double *moist, - double *ice, - double *dz, - double *T) -/****************************************************************************** - read_initial_snow_thermal Keith Cherkauer January 12, 1999 - - This routine reads parameters for initial soil thermal conditions - if defined - - Parameters Read from File: - TYPE NAME UNITS DESCRIPTION - int cellnum N/A number of current cell - int veg N/A vegetation number within cell - int band N/A snow band number within cell - double *fdepth m depth of freezing and thawing fronts - double *moist mm/mm moisture content of each soil layer - double *ice mm/mm ice content of each soil layer - double *thermdepths m depth from surface of each thermal soln node - double *T C temperature at each thermal soln node depth - - Parameters Computed from those in the File: - TYPE NAME UNITS DESCRIPTION - double *dz m thickness of soil layer between thermal nodes - - Modifications: - -**********************************************************************/ -{ - extern option_struct options; - extern debug_struct debug; - - char ErrStr[MAXSTRING]; - char tmpstr[MAXSTRING]; - int i, tmpint, tmpveg, tmpband, index, tmpnodes; - double *thermdepths; - - rewind(initsoil); - - fscanf(initsoil, "%s", tmpstr); - if(tmpstr[0] != '#') { - index = atoi(tmpstr); - fscanf(initsoil, "%i %i", &tmpveg, &tmpband); - } - else index = tmpveg = tmpband = -999; - while(index != cellnum && veg != tmpveg && band != tmpband) { - fgets(tmpstr,MAXSTRING,initsoil); - fscanf(initsoil, "%s", tmpstr); - if(tmpstr[0] != '#') { - index = atoi(tmpstr); - fscanf(initsoil, "%i %i", &tmpveg, &tmpband); - } - else index = tmpveg = tmpband = -999; - } - - if(!feof(initsoil) && options.FROZEN_SOIL) { - - /** Read Soil Thermal Initialization File for Activated Frozen Soil Model **/ - - thermdepths = (double*)calloc(Nnodes,sizeof(double)); - - fscanf(initsoil, "%i", &tmpnodes); - if(Nnodes!=tmpnodes) { - sprintf(ErrStr,"Number of nodes defined in soil thermal initialization file (%i), not equal to number of nodes defined in model (%i).",tmpnodes,Nnodes); - } - for(i=0;i<2;i++) fscanf(initsoil, "%lf", &fdepth[i]); - for(i=0;i0;j--) { - thermdepths[j] -= thermdepths[j-1]; - thermdepths[j] = rint(thermdepths[j] * 10000.) / 10000.; - } - - sum = 0; - dz[0] = rint(thermdepths[1] * 10000.) / 10000.; - for(j=1;j -#include -#include - -static char vcid[] = "$Id$"; - -void read_rosemount(atmos_data_struct *temp, - FILE *snowf, - int *nrecs, - int starthour, - int dt, - int file_dt, - int fileskip) -/********************************************************************** - read_rosemount Keith Cherkauer January 7, 1997 - - This subroutine reads in hourly meteorolgical data from the - Rosemount weather file (SHAW model hourly met data format). - - NOTE: Use snow file flag to identify this file. - -**********************************************************************/ -{ - extern debug_struct debug; - extern param_set_struct param_set; - - int i, n, rec, maxline = 210; - int FIRST = TRUE; - int day, year, hour; - int store_rec; - int skip_bytes; - char str[210]; - double junk; - - /** locate starting record **/ - skip_bytes = (int)((float)(dt * fileskip)) / (float)file_dt - 1; - if((dt * fileskip) % (24 / file_dt) > 0) - nrerror("Currently unable to handle a model starting date that does not correspond to a line in the forcing file."); - for(i=0;i100) temp[rec].rel_humid=100; - else if(temp[rec].rel_humid<0) temp[rec].rel_humid=0; - if(temp[rec].air_temp>150 || temp[rec].air_temp<-150) { - fprintf(stderr,"ERROR: Invalid Air Temperature %lf\n", - temp[rec].air_temp); - exit(0); - } - if(temp[rec].wind<0) temp[rec].wind=0; - rec++; - - if(file_dt < dt) { - /** Time Step in Forcing File Finer than Used by Model: - Skip Records **/ - for(i=0;i
dt) { - /** Time step used by model finer than that used in forcing file: - Repeat Data Into Extra Columns **/ - store_rec = rec - 1; - for(i=1;i -#include -#include - -static char vcid[] = "$Id$"; - -void read_sawd(atmos_data_struct *temp, - FILE *sawdf, - int *nrecs, - int dt, - int file_dt, - int fileskip) -/********************************************************************** - read_sawd Keith Cherkauer July 25, 1996 - - This routine reads in atmospheric data values from surface - airways gridded hourly data files. If time step is less than - hourly, extra data is skipped. - - Input Output - Units Units - air_temp - C C - rel_humid - % % - tskc - % fraction - wind - m/s m/s - pressure - Pa kPa - -**********************************************************************/ -{ - extern option_struct options; - extern debug_struct debug; - extern param_set_struct param_set; - - int i, n, rec, maxline = 210; - int fixcnt, headcnt; - int store_rec; - int skip_bytes; - char str[210],jnkstr[MAXSTRING]; - - /** Count Records **/ - n = 0; - while (fgets(str,maxline,sawdf) != '\0') n++; - printf("nrecs = %d\n",n); - - rewind(sawdf); - - /** Check for Header, and Skip **/ - fgets(jnkstr,MAXSTRING,sawdf); - if((jnkstr[0]<48 || jnkstr[0]>57) && jnkstr[0]!=46) { - fgets(jnkstr,MAXSTRING,sawdf); - fprintf(stderr,"SAWD ... skipping header\n"); - } - - /** find simulation start record **/ - skip_bytes = (int)((float)(dt * fileskip)) / (float)file_dt - 1; - if((dt * fileskip) % (24 / file_dt) > 0) - nrerror("Currently unable to handle a model starting date that does not correspond to a line in the forcing file."); - for(i=0;i0) { - temp[rec].pressure = temp[rec-1].pressure; - fixcnt++; - } - fscanf(sawdf,"%s",str); - temp[rec].wind = atof(str); - fscanf(sawdf,"%s",str); - temp[rec].rel_humid = atof(str); - if(temp[rec].rel_humid == 0. && rec>0) { - temp[rec].rel_humid = temp[rec-1].rel_humid; - fixcnt++; - } - fscanf(sawdf,"%s",str); - temp[rec].tskc = atof(str) / 100.; - rec++; - - if(file_dt < dt) { - /** Time Step in Forcing File Finer than Used by Model: - Skip Records **/ - for(i=0;i
dt) { - /** Time step used by model finer than that used in forcing file: - Repeat Data Into Extra Columns **/ - store_rec = rec - 1; - for(i=1;i0) { - fprintf(stderr,"WARNING: Had to fix %i values in sawd file.\n",fixcnt); - } - - if(rec < *nrecs) { - fprintf(stderr,"WARNING: Not enough records in the SAWD ASCII forcing file to run the number of records defined in the global file. Check forcing file time step (%i), and global file. Number of records being modified to stop model when available data has run out.\n",file_dt); - *nrecs = rec; - } - - param_set.WIND = param_set.AIR_TEMP = param_set.PRESSURE = TRUE; - param_set.TSKC = TRUE; - param_set.REL_HUMID = TRUE; - -} diff --git a/read_sawd_binary.c b/read_sawd_binary.c deleted file mode 100644 index 6dd627704..000000000 --- a/read_sawd_binary.c +++ /dev/null @@ -1,105 +0,0 @@ -#include -#include -#include - -static char vcid[] = "$Id$"; - -void read_sawd_binary(atmos_data_struct *temp, - FILE *sawdf, - int *nrecs, - int dt, - int file_dt, - int fileskip) -/********************************************************************** - read_sawd_binary Keith Cherkauer April 26, 1998 - - This routine reads in atmospheric data values from surface - airways binary gridded hourly data files. - - Input Output - Units Units - air_temp - C*100 C - pressure - Pa*100 kPa - wind - m/s*1000 m/s - rel_humid - %*100 % - tskc - %*10 fract - -**********************************************************************/ -{ - extern debug_struct debug; - extern param_set_struct param_set; - - int i, n, rec; - int fixcnt; - int store_rec; - int skip_bytes; - char tmpmem[20]; - short int values[5]; - - /** locate starting record **/ - skip_bytes = (int)((float)(dt * fileskip)) / (float)file_dt - 1; - if((dt * fileskip) % (24 / file_dt) > 0) - nrerror("Currently unable to handle a model starting date that does not correspond to a line in the forcing file."); - fseek(sawdf,skip_bytes*12,SEEK_SET); - if(feof(sawdf)) - nrerror("No data for the specified time period in the Daily Precipitation forcing file. Model stopping..."); - - /** Read forcing data **/ - fixcnt = 0; - rec = 0; - while ( !feof(sawdf) && (rec < *nrecs) ) { - fread(values,sizeof(short int),5,sawdf); - temp[rec].air_temp = (double)values[0]/100.; - temp[rec].pressure = (double)values[1]/100.; - if(temp[rec].pressure == 0. && rec>0) { - temp[rec].pressure = temp[rec-1].pressure; - fixcnt++; - } - temp[rec].wind = (double)values[2]/1000.; - temp[rec].rel_humid = (double)values[3]/100.; - if(temp[rec].rel_humid == 0. && rec>0) { - temp[rec].rel_humid = temp[rec-1].rel_humid; - fixcnt++; - } - if(temp[rec].rel_humid > 100.) { - temp[rec].rel_humid = 100.; - fixcnt++; - } - temp[rec].tskc = (double)values[4] / 1000.; - rec++; - - if(file_dt < dt) { - /** Time Step in Forcing File Finer than Used by Model: - Skip Records **/ - for(i=0;i
dt) { - /** Time step used by model finer than that used in forcing file: - Repeat Data Into Extra Columns **/ - store_rec = rec - 1; - for(i=1;i0) { - fprintf(stderr,"WARNING: Had to fix %i values in sawd file.\n",fixcnt); - } - - if(rec < *nrecs) { - fprintf(stderr,"WARNING: Not enough records in the SAWD Binary forcing file to run the number of records defined in the global file. Check forcing file time step (%i), and global file. Number of records being modified to stop model when available data has run out.\n",file_dt); - *nrecs = rec; - } - - param_set.WIND = param_set.AIR_TEMP = param_set.PRESSURE = TRUE; - param_set.TSKC = TRUE; - param_set.REL_HUMID = TRUE; - -} diff --git a/read_snowmodel.c b/read_snowmodel.c deleted file mode 100644 index 5ab401144..000000000 --- a/read_snowmodel.c +++ /dev/null @@ -1,94 +0,0 @@ -#include -#include -#include - -static char vcid[] = "$Id$"; - -void read_snowmodel(atmos_data_struct *temp, - FILE *snowf, - int nrecs, - int dt, - int file_dt, - int fileskip) -/********************************************************************** - read_snowmodel Keith Cherkauer July 25, 1996 - - This routine reads in atmospheric data values from the output - files for the NWS snow melt model. - - Modifications: - 7-25-96 Routine modified to read data into structure created - for time steps of less than 1 day. KAC - -**********************************************************************/ -{ - extern param_set_struct param_set; - - int i, n, rec, maxline = 210; - int store_rec; - int skip_bytes; - char str[210]; - - /** locate starting record **/ - skip_bytes = (int)((float)(dt * fileskip)) / (float)file_dt - 1; - if((dt * fileskip) % (24 / file_dt) > 0) - nrerror("Currently unable to handle a model starting date that does not correspond to a line in the forcing file."); - for(i=0;i dt) { - /** Time step used by model finer than that used in forcing file: - Repeat Data Into Extra Columns **/ - store_rec = rec-1; - for(i=1;i