利用VIC模型将气象数据日值插为小时值

搭建水文模型,最让人头疼的就是数据不全了,国内能够下载的气象站数据只有温度、相对湿度、风速和降雨,而且还全都是日值。而更为不幸的是,我最近在用的DHSVM模型需要的是小时值,而且还需要短波辐射和长波辐射数据。如何有效利用现有的气象站数据为水文模型做输入就成为我迫切要解决的问题。方法可能有很多,我这里实现了一个解决方案:利用VIC模型来进行插值转换。
题外话:按理说,VIC模型和DHSVM模型与华盛顿大学都渊源很深,可是VIC模型中有相应的插值模块,但是DHSVM中却没有。
VIC的插值算法大体是(可能我的表述不够准确):
(1)首先VIC利用每日平均气温和平均相对湿度估算出每日平均水汽压,在与每日最高温度和最低温度数据以及每日降水数据一起作为输入,使用蒙大拿大学开发的MT-CLIM模型,算出每小时短波辐射 。同时,会算出逐小时云量(后面会用到)
(2)根据上一步算出的每小时短波辐射判断出每日最高温度和最低温度时刻,建立Hermite多项式,保证在最低和最高时刻的一阶导数为0(从而保证了最低温度和最高温度与输入一致),插值得到逐小时的温度。
(3)VIC自身会通过经验公式,利用逐小时的温度和逐日的相对湿度计算逐小时的水汽压,然后利用逐小时的云量、温度和水汽压直接算出逐小时的长波辐射。
以上,是在我详细查看了VIC相关源码的基础上做的一点总结,可能不够准确。
在了解了VIC的工作原理之后,需要做的就是编写相应的接口程序来调用VIC的相应模块生成我们所需要的结果。而且,最好不要修改VIC的源码。所以,我就花了些时间,编写了一个 调用程序(源码贴在后面),没有改动VIC的源文件,只需要将我写的代码文件与VIC的源码文件一起编译即可。
它的使用方法是:
(1)有六个可选项,分别是,”temp”, “prec”, “srad”, “lrad” , “rh”, “wind” ,表示,你希望程序输出温度、降水、短波辐射、长波辐射、相对湿度或者风速。可以多选。如果要计算长波辐射、逐小时温度或者短波辐射,那么需要准备的输入数据是:逐日的降水、相对湿度、最高温度、最低温度和平均温度。如果要计算逐小时的降水,则只需要逐日的降水数据即可。实际上,此程序只会计算短波、长波和温度,其余变量只是简单的平均或填充一天。
(2)此外,用户还需要设置一些必要的参数,程序会在当前目录下寻找一个名为“parameter.txt”的文件,该文件内容如下:

annual_prec    383.96
elevation    4432.4
slope    0
aspect    0
ehoriz    0
whoriz    0
StartYear    2000
StartMonth    1
StartDay    1
latitude    28.58
sw_prec_thresh    0
mtclim_swe_corr     0
have_vp     1
其各项的含义如下,
annual_prec :年平均降水量(mm)
elevation : Site elevation, meters
slope : Site slope, degrees
aspect : Site aspect, degrees (0=N,90=E,180=S,270=W)
ehoriz :Site east horizon, degrees
whoriz :Site west horizon, degrees
StartYear :起始年
StartMonth:起始月
    StartDay:起始日
    latitude :Site latitude, degrees (- for south)
    sw_prec_thresh:Minimum daily precipitation [mm] that can cause “dimming” of incoming shortwave radiation
    mtclim_swe_corr :TRUE = correct MTCLIM’s downward shortwave  radiation estimate in presence of snow
    have_vp :1表示用户提供水汽压数据,0表示使用本程序计算水汽压注意: 通常,slope、aspect、ehoriz以及whoriz这四项是未知的,我查阅了VIC的源码,其在计算时,都是简单地设定为0,如下
1

因此,我在案例中也将这四项全部设为0,至于这四项具体影响有多大。进过简单测试,slope的影响很大,aspect几乎没有影响。slope为0应是假定没有地形影响,以这样的数据再在DHSVM中进行全流域distibute可能会更加合理。

我编写的源码如下:
注意:一个可行的编译运行方法是,下载VIC源码,并解压,将下面的代码复制到名为“vic_invoke.c”的文本文件中,并移动到解压后的目录。Linux下,输入命令:gcc vic_invoke.c Svp.c mtclim_wrapper.c mtclim_vic.c Vicerror.c Close_files.c Nrerror.c Compress_files.c Calc_air_temperature.c Calc_longwave.c -o vic_invoke -Wall -g -I . -lm即可顺利编译,然后./vic_invoke [option]运行。
如果要在windows下运行该程序,可下载安装MinGW,然后将下面源码中的LINUX宏设置为0,然后切换目录到VIC文件夹下,输入同样的编译命令即可生成win32可执行程序。
windows下直接运行,会弹出如下界面
2

所以,需要启动cmd,切换到数据、程序所在目录,以命令行的形式运行。

/*
 * =====================================================================================
 *       Filename:  zhb_invoke.c
 *    Description:  invoke VIC_routines to compute hourly parameters
 *        Version:  1.0
 *        Created:  2013年10月12日 16时58分17秒
 *       Revision:  2013年10月17日 17时37分29秒
 *       Compiler:  gcc
 *         Author:  zhanghongbo (zhb), zhanghongbo@itpcas.ac.cn
 *   Organization:  ITPCAS
 * =====================================================================================
 */
/* Detailed Description:
 * (1) Temperature is calculated based on Tmax, Tmin and hourly shortwave
 * (2) Short wave radiation is computed by MT-CLIM
 * (3) Wind is just constant duting the day
 * (4) Prec is averaged to every hour of a day
 * (5) The format of input “temp” file should be like:
 * max min mean
 * 2.1-10.1 -4.4
 * 3.2-10.3 -3.9
 * …… …
 * and its name should be “temp.txt” if you didn’t modify source codes.
 * (6) The format of input “humidity”, “wind” and “prec” file should be like:
 * humidity // any name is OK, just skipped, but this line is needed
 * 95
 * 60
 * 36
 * …
 * and theirs names must be “rh.txt”, “wind.txt” and “prec.txt”
 * respectively, if you didn’t modify source codes.
 * */
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>
#if LINUX
#include <termios.h>
#endif
#include <unistd.h>
#include <vicNl.h>
#define StrSize 255// Buffer size of a string
#define HeadLines 1// Header lines count of input data file
#define DBG 1// set this to print debug info
#define KELVIN 273.15
#define LINUX 0 // 0 for windows
char* LowerCase(char *);
int CheckOptions(char *);
void PrintOptions(void);
int CheckInputs(void);
int GetRecords(char *, int);
int Lines(char *);
int MaxLineLen(char *);
void err(char *format, …);
int SkipLines(FILE *f, int num, int linebuf);
int InitInputs();
int ComputeVP();
int YearDays(int);
int LeapYear(int);
int MonthDays(int mon, int year);
int InitParams();
int DayInYear(int, int, int);
int InitDates();
int ComputeSrad();
int ComputeTemp();
int ComputeHourlyVP();
int ComputeLrad();
int ComputePrec();
int ComputeWind();
int ComputeRH();
void Output(char *, double *);
void PrintUsage();
int getchar_ontime();
enum DataType {character=0, ucharacter, shortint, ushortint,
longint, ulongint, longlongint, ulonglongint, floatpoint,
doublefloat};
typedef struct {
int Srad;// short wave radiation
int Temp;// temperature
int Lrad;// long wave radiation
int Wind;// wind spedd
int Prec;// precipitation
int RH; // Relative Humidity
char TempFile[StrSize+1]; // input file for Temp
char RhFile[StrSize+1]; // input file for Relative Humidity
char WindFile[StrSize+1]; // input file for Wind
char PrecFile[StrSize+1]; // input file for Prec
char ParamFile[StrSize+1]; // input file of parameters for MTCLIM
char VpFile[StrSize+1];
} OPTIONS;
// fake options to suit for vic
option_struct options;
Error_struct Error;
OPTIONS Options = {
0,0,0,0,0,0,
“temp.txt”,
“rh.txt”,
“wind.txt”,
“prec.txt”,
“parameter.txt”,
“vp.txt”
};
struct ParamsStruct {
char name[StrSize+1];
int id;
int supplied;
};
enum ParamsSeq {APREC = 0, ELEV, SLOPE, ASPECT,
EHORIZ, WHORIZ, SYEAR, SMONTH, SDAY, LAT,
SW_PREC_THRESH, MTCLIM_SWE_CORR, HAVE_VP
};
struct ParamsStruct Params[] = {
{“annual_prec”, APREC, 0},
{“elevation”, ELEV, 0},
{“slope”, SLOPE, 0},
{“aspect”, ASPECT, 0},
{“ehoriz”, EHORIZ, 0},
{“whoriz”, WHORIZ, 0},
{“startyear”, SYEAR, 0},
{“startday”, SDAY, 0},
{“startmonth”, SMONTH, 0},
{“latitude”, LAT, 0},
{“sw_prec_thresh”, SW_PREC_THRESH, 0},
{“mtclim_swe_corr”, MTCLIM_SWE_CORR, 0},
{“have_vp”, HAVE_VP, 0}
};
double *Tmax;// max temperature(c degree)
double *Tmin;// min temperature(c degree)
double *Tair;// mean temperature(c degree)
float *RH;// relative humidity(fraction)
double *Wind;// wind speed(m/s)
double *Prec;// precipitation(mm)
// prepare vars for VIC to use
double *vhourlyrad;
double *vlrad;
double *vtair; // hourly temperature
double *vprec;
double *vtskc;// cloudness
double *vdaily_vp;
double *vhourly_vp;
int *vtmaxhour;
int *vtminhour;
int Ndays;
int have_dewpt;
int have_shortwave;
int hour_offset;
int have_vp; // 1 means we have site vapor pressure, 0 means the opposite
double elevation;// Site elevation, meters
double slope;// Site slope, degrees
double aspect;// Site aspect, degrees (0=N,90=E,180=S,270=W)
double ehoriz;// Site east horizon, degrees
double whoriz;// Site west horizon, degrees
double annual_prec;// mm
double latitude;// latitude
char mtclim_swe_corr;// TRUE = correct MTCLIM’s downward shortwave radiation estimate in presence of snow
float sw_prec_thresh;// Minimum daily precipitation [mm] that can cause “dimming” of incoming shortwave radiation
dmy_struct *vdmy;// vic-type date struct
int startyear;
int startmon;
int startday;
int main(int argn, char *args[]){
int i = 0;
++args;
int count = 0;
for (i=1; i<argn; i++) {
count += CheckOptions(LowerCase(*args++));
}
if (count == 0) {
err(“ERROR: No options are detected. You must input at least one option.\n\n”);
PrintUsage();
printf(“Press any key to continue…\n”);
getchar_ontime();
exit(1);
}
PrintOptions();
if (CheckInputs() != 0) {
fprintf(stderr, “Inputs check falied\nProgram exits\n”);
exit(1);
};
if (!InitInputs()) {
err(“ERROR: Initialisation of input data failed.\n”);
exit(1);
}
if (Options.Lrad) {
ComputeSrad();
ComputeTemp();
ComputeLrad();
Output(“hourly_lrad.txt”, vlrad);
}
else {
if (Options.Temp) {
ComputeSrad();
ComputeTemp();
} else {
if (Options.Srad) {
ComputeSrad();
}
}
}
if (Options.Temp) {
Output(“hourly_temp.txt”, vtair);
}
if (Options.Srad) {
Output(“hourly_srad.txt”, vhourlyrad);
}
if (Options.Prec)
ComputePrec();
if (Options.Wind)
ComputeWind();
if (Options.RH) {
ComputeRH();
}
return 0;
}
int CheckOptions(char *op){
if (strcmp(op, “temp”) == 0) {
Options.Temp = 1;
return 1;
}
else if (strcmp(op, “wind”) == 0) {
Options.Wind = 1;
return 1;
}
else if (strcmp(op, “srad”) == 0) {
Options.Srad = 1;
return 1;
}
else if (strcmp(op, “lrad”) == 0) {
Options.Lrad = 1;
return 1;
}
else if (strcmp(op, “prec”) == 0) {
Options.Prec = 1;
return 1;
}
else if (strcmp(op, “rh”) == 0) {
Options.RH = 1;
return 1;
}
return 0;
}
int CheckInputs(void){
if (Options.Temp || Options.Srad || Options.Lrad) {
if (access(Options.TempFile, F_OK) != 0) {
fprintf(stderr, “Error: input file \”%s\” is not found!\n”, Options.TempFile);
return -1;
}
if (access(Options.RhFile, F_OK) != 0) {
fprintf(stderr, “Error: input file \”%s\” is not found!\n”, Options.RhFile);
return -1;
}
if (access(Options.PrecFile, F_OK) != 0) {
fprintf(stderr, “Error: input file \”%s\” is not found!\n”, Options.PrecFile);
return -1;
}
if (access(Options.ParamFile, F_OK) != 0) {
fprintf(stderr, “Error: input file \”%s\” is not found!\n”, Options.ParamFile);
return -1;
}
}
if (Options.Wind) {
if (access(Options.WindFile, F_OK) != 0) {
fprintf(stderr, “\”wind\” option is set, but %s is not found!\n”, Options.WindFile);
return -1;
}
}
if (Options.Prec) {
if (access(Options.PrecFile, F_OK) != 0) {
fprintf(stderr, “\”prec\” option is set, but %s is not found!\n”, Options.PrecFile);
return -1;
}
}
if (Options.RH) {
if (access(Options.RhFile, F_OK) != 0) {
fprintf(stderr, “\”rh\” option is set, but %s is not found!\n”, Options.RhFile);
return -1;
}
}
return 0;
}
void PrintOptions(void){
printf(“Options are listed below.\n”);
printf(“Srad: %d\n”, Options.Srad);
printf(“Temp: %d\n”, Options.Temp);
printf(“Lrad: %d\n”, Options.Lrad);
printf(“Wind: %d\n”, Options.Wind);
printf(“Prec: %d\n”, Options.Prec);
printf(“RH: %d\n”, Options.RH);
}
void PrintUsage(){
printf(“Usage: vic_invoke [option]\n”);
printf(“Valid options are listed below:\n”);
printf(“temp\t\t:compute hourly temperaure\n\
srad\t\t:compute hourly short wave radiation\n\
lrad\t\t:compute hourly long wave radiation\n\
prec\t\t:compute hourly precipitation\n\n”);
printf(“Note:\n”);
printf(“If you set \”temp\”, \”srad\” or \”lrad\”, you should prepare four input files, they are:\n”);
printf(“(1) temp.txt\t\t:max min mean of temperature\n\
(2) prec.txt\t\t:prec\n\
(3) rh.txt\t\t:Relative Humidity\n\
(4) parameter.txt\t:twelve \”key value\” paired parameters\n”);
printf(“\n”);
printf(“Author:  zhb\nContact:  zhanghongbo@itpcas.ac.cn\n”);
printf(“\n\n”);
}
char* LowerCase(char *str) {
char *sav = str;
int n = strlen(str);
int i = 0;
for (; i<n; i++) {
if (*str >= ‘A’ && *str <= ‘Z’)
 *str += 32;
++str;
}
return sav;
}
// read necessary input files and initialize needed arrays
int InitInputs(){
if (Options.Temp || Options.Srad || Options.Lrad) {
int lines = Lines(Options.TempFile);
if (Lines(Options.RhFile) != lines) {
err(“Error: RH file and temp file have different lines’ count\n”);
return 0;
}
if (Lines(Options.PrecFile) != lines) {
err(“Error: Prec file and temp file have different lines’ count\n”);
return 0;
}
int BufSize = lines – HeadLines;
Tmax = calloc(BufSize, sizeof(double));
Tmin = calloc(BufSize, sizeof(double));
Tair = calloc(BufSize, sizeof(double));
RH = calloc(BufSize, sizeof(float));
Prec = calloc(BufSize, sizeof(double));
BufSize = MaxLineLen(Options.TempFile)+1;
char line[BufSize];
FILE *f_temp = fopen(Options.TempFile, “r”);
if (!SkipLines(f_temp, HeadLines, BufSize)) {
fclose(f_temp);
err(“Current file is %s\n”, Options.TempFile);
return 0;
}
FILE *f_rh = fopen(Options.RhFile, “r”);
if (!SkipLines(f_rh, HeadLines, BufSize)) {
fclose(f_rh);
err(“Current file is %s\n”, Options.RhFile);
return 0;
}
FILE *f_prec = fopen(Options.PrecFile, “r”);
if (!SkipLines(f_prec, HeadLines, BufSize)) {
fclose(f_prec);
err(“Current file is %s\n”, Options.PrecFile);
return 0;
}
lines = lines – HeadLines;
Ndays = lines;
int BufRh = MaxLineLen(Options.RhFile)+1;
char lineRh[BufRh];
int BufPrec = MaxLineLen(Options.PrecFile)+1;
char linePrec[BufPrec];
while (lines–) {
fgets(line, BufSize, f_temp);
sscanf(line, “%lf %lf %lf”, Tmax++, Tmin++, Tair++);
#if DBG
printf(“Temp %d: %lf %lf %lf\n”, lines, *(Tmax-1), *(Tmin-1), *(Tair-1));
#endif
fgets(lineRh, BufRh, f_rh);
sscanf(lineRh, “%f”, RH++);
#if DBG
printf(“RH %d: %f\n”, lines, *(RH-1));
#endif
fgets(linePrec, BufPrec, f_prec);
sscanf(linePrec, “%lf”, Prec++);
#if DBG
printf(“Prec %d: %lf\n”, lines, *(Prec-1));
#endif
}
fclose(f_temp);
fclose(f_rh);
fclose(f_prec);
// recovery
Tmax -= Ndays;
Tmin -= Ndays;
Tair -= Ndays;
RH -= Ndays;
Prec -= Ndays;
// covert c to kelvin, of course, here is a joke
#if TO_KEVIN
int i = 0;
for (; i<Ndays; i++) {
Tmax[i] += KELVIN;
Tmin[i] += KELVIN;
Tair[i] += KELVIN;
}
#endif
// prepare for VIC
vtskc = calloc(Ndays*24, sizeof(double));
vhourlyrad = calloc(Ndays*24, sizeof(double));
if (!InitParams()) {
err(“Parameters check failed, program exits.\n”);
return 0;
}
}
if (Options.Prec) {
if (Prec == NULL) {
int lines = Lines(Options.PrecFile) – HeadLines;
if (Ndays != 0) {
if (Ndays != lines – HeadLines) {
err(“ERROR: num of recs in %s doesn’t match Ndays(%d)\n”, Options.WindFile, Ndays);
return 0;
}
}
Ndays = lines;
int linebuf = MaxLineLen(Options.PrecFile) + 1;
FILE *f_prec = fopen(Options.PrecFile, “r”);
if (!SkipLines(f_prec, HeadLines, linebuf)) {
fclose(f_prec);
err(“Current file is %s\n”, Options.PrecFile);
return 0;
}
Prec = calloc(Ndays, sizeof(double));
char linePrec[linebuf];
while (lines–) {
fgets(linePrec, linebuf, f_prec);
sscanf(linePrec, “%lf”, Prec++);
#if DBG
printf(“Prec %d: %lf\n”, lines, *(Prec-1));
#endif
}
Prec -= Ndays;
fclose(f_prec);
}
}
if (Options.Wind) {
int lines = Lines(Options.WindFile) – HeadLines;
if (Ndays != 0) {
if (Ndays != lines) {
err(“ERROR: num of recs in %s doesn’t match Ndays(%d)\n”, Options.WindFile, Ndays);
return 0;
}
}
Ndays = lines;
int linebuf = MaxLineLen(Options.WindFile) + 1;
char linewind[linebuf];
FILE *f_wind = fopen(Options.WindFile, “r”);
if (!SkipLines(f_wind, HeadLines, linebuf)) {
fclose(f_wind);
err(“Current file is %s\n”, Options.WindFile);
return 0;
}
Wind = calloc(Ndays, sizeof(double));
while (lines–) {
fgets(linewind, linebuf, f_wind);
sscanf(linewind, “%lf”, Wind++);
#if DBG
printf(“WindSpeed %d: %lf\n”, lines, *(Wind-1));
#endif
}
Wind -= Ndays;
fclose(f_wind);
}
if (Options.RH) {
if (RH == NULL) {
int lines = Lines(Options.RhFile) – HeadLines;
if (Ndays != 0) {
if (lines != Ndays) {
err(“Error: num of recs in %s doesn’t match Ndays(%d)\n”, Options.RhFile, Ndays);
return 0;
}
}
Ndays = lines;
int linebuf = MaxLineLen(Options.RhFile) + 1;
char linerh[linebuf];
FILE *f_rh = fopen(Options.RhFile, “r”);
if (!SkipLines(f_rh, HeadLines, linebuf)) {
fclose(f_rh);
err(“Current file is %s\n”, Options.RhFile);
return 0;
}
RH = calloc(Ndays, sizeof(double));
while (lines–) {
fgets(linerh, linebuf, f_rh);
sscanf(linerh, “%f”, RH++);
#if DBG
printf(“RH %d: %f\n”, lines, *(RH-1));
#endif
}
RH -= Ndays;
fclose(f_rh);
}
}
return 1;
}
int InitDates(){
vdmy = (dmy_struct*)calloc(Ndays*24, sizeof(dmy_struct));
int i = 0;
int year = startyear;
int mon = startmon;
int day_in_year = DayInYear(startday, mon, year);
int day = startday;
int hour = 0;
for (; i<Ndays*24; i++) {
vdmy[i].day = day;
vdmy[i].day_in_year = day_in_year;
vdmy[i].hour = hour;
vdmy[i].month = mon;
vdmy[i].year = year;
++hour;
if (hour == 24) {
++day;
++day_in_year;
hour = 0;
if (day > MonthDays(mon, year)) {
day = 1;
++mon;
if (mon > 12) {
mon = 1;
day_in_year = 1;
++year;
}
}
}
#if DBG_DATE
printf(“Time %d: %d – %d – %d (%d), %d\n”, i, vdmy[i].year,
vdmy[i].month, vdmy[i].day, vdmy[i].day_in_year, vdmy[i].hour);
#endif
}
return 1;
}
int InitParams(){
int lines = Lines(Options.ParamFile);
int linebuf = MaxLineLen(Options.ParamFile)+1;
FILE *f_pm = fopen(Options.ParamFile, “r”);
char line[linebuf];
char tmp[linebuf];
int i = 0;
int npm = sizeof(Params) / sizeof(Params[0]);
char sep[] = ” \t”;
for (; i<lines; i++) {
fgets(line, linebuf, f_pm);
strcpy(tmp, line);
char *first = strtok(tmp, sep);
int j = 0;
for(; j<npm; j++) {
if(strcmp(Params[j].name, LowerCase(first)) == 0) {
char *sec = strtok(NULL, sep);
Params[j].supplied = 1;
switch(Params[j].id) {
case APREC:
annual_prec = atof(sec);
break;
case ELEV:
elevation = atof(sec);
break;
case SLOPE:
slope = atof(sec);
break;
case ASPECT:
aspect = atof(sec);
break;
case EHORIZ:
ehoriz = atof(sec);
break;
case WHORIZ:
whoriz = atof(sec);
break;
case SYEAR:
startyear = atoi(sec);
break;
case SMONTH:
startmon = atoi(sec);
break;
case SDAY:
startday = atoi(sec);
break;
case LAT:
latitude = atof(sec);
break;
case SW_PREC_THRESH:
sw_prec_thresh = atof(sec);
break;
case MTCLIM_SWE_CORR:
mtclim_swe_corr = atof(sec);
break;
case HAVE_VP:
have_vp = atoi(sec);
break;
default:
err(“ERROR: unidentified parameter ID\n”);
}
#if DBG
printf(“%s: %lf\n”, Params[j].name, atof(sec));
#endif
break;
}
}
}
// check which param is missing
for (i=0; i<npm; i++) {
if (Params[i].supplied == 0) {
err(“Error, missing parameter: %s, you must specify it.\n”, Params[i].name);
return 0;
}
}
if (have_vp) {
if (access(Options.VpFile, F_OK) != 0) {
fprintf(stderr, “Error: input file \”%s\” is not found!\n”, Options.VpFile);
return 0;
}
}
InitDates();
return 1;
}
// compute daily Vapor Pressure
int ComputeVP(){
int i = 0;
double tmp;
for (; i<Ndays; i++) {
tmp = svp(Tair[i]);
vdaily_vp[i] = RH[i] * tmp / 100;
#if DBG
printf(“vp %d: %10.6lf, with RH: %5.2f, Tair: %5.2lf, svp: %10.6lf\n”, i, vdaily_vp[i], RH[i], Tair[i], tmp);
#endif
}
// then we have vapor pressure
return 1;
}
int ComputeSrad() {
if (have_vp) {
int lines = Lines(Options.VpFile);
lines -= HeadLines;
if (lines != Ndays) {
err(“ERROR: num of recs in %s dosen’t match Ndays”,
Options.VpFile);
exit(1);
}
int linebuf = MaxLineLen(Options.VpFile) + 1;
FILE *f_vp = fopen(Options.VpFile, “r”);
if (!SkipLines(f_vp, HeadLines, linebuf)) {
fclose(f_vp);
err(“Current file is %s\n”, Options.VpFile);
exit(1);
}
vdaily_vp = calloc(Ndays, sizeof(double));
char linevdaily_vp[linebuf];
while (lines–) {
fgets(linevdaily_vp, linebuf, f_vp);
sscanf(linevdaily_vp, “%lf”, vdaily_vp++);
#if DBG
printf(“vdaily_vp %d: %lf\n”, lines, *(vdaily_vp-1));
#endif
}
vdaily_vp -= Ndays;
fclose(f_vp);
}
else
ComputeVP();
// use MTCLIM
hour_offset = 0; // maybe not so exact
options.SW_PREC_THRESH = sw_prec_thresh;
options.MTCLIM_SWE_CORR = mtclim_swe_corr;
// the two params below are temporarily set internally
options.VP_ITER = VP_ITER_ANNUAL;
options.COMPRESS = 0;
have_dewpt = 2;// because we have RH, so here must be 2
have_shortwave = 0;// here must be 0, for we don’t have site srad
mtclim_wrapper(have_dewpt, have_shortwave, hour_offset,
elevation, slope, aspect, ehoriz, whoriz, annual_prec, latitude,
Ndays, vdmy, Prec, Tmax, Tmin, vtskc, vdaily_vp, vhourlyrad);
#if DBG
int dbg = 0;
for (; dbg < Ndays*24; dbg++) {
printf(“srad %d: %lf\n”, dbg, vhourlyrad[dbg]);
}
#endif
return 1;
}
int ComputeTemp(){
vtmaxhour = calloc(Ndays, sizeof(int));
vtminhour = calloc(Ndays, sizeof(int));
vtair = calloc(Ndays*24, sizeof(double));
set_max_min_hour(vhourlyrad, Ndays, vtmaxhour, vtminhour);
#if DBG_MAX_HOUR
int dbg1 = 0;
for (; dbg1 < Ndays; dbg1++) {
printf(“MaxTempHour:  %d, MinTempHour %d\n”, vtmaxhour[dbg1], vtminhour[dbg1]);
}
#endif
HourlyT(1, Ndays, vtmaxhour, Tmax, vtminhour, Tmin, vtair);
#if DBG
int dbg = 0;
for (; dbg < Ndays*24; dbg++) {
printf(“temp %d: %lf\n”, dbg, vtair[dbg]);
}
#endif
return 1;
}
int ComputeHourlyVP(){
vhourly_vp = calloc(Ndays*24, sizeof(double));
int i = 0;
int day = 0;
for (; i<Ndays*24; i++) {
day = i / 24;
vhourly_vp[i] = RH[day] * svp(vtair[i]) / 100;
}
#if DBG
int dbg = 0;
for (; dbg < Ndays*24; dbg++) {
printf(“Hourly_VP %d: %lf\n”, dbg, vhourly_vp[dbg]);
}
#endif
return 1;
}
int ComputeLrad(){
ComputeHourlyVP();
vlrad = calloc(Ndays*24, sizeof(double));
int i = 0;
options.LW_CLOUD = LW_CLOUD_DEARDORFF;
options.LW_TYPE = LW_TVA;
for (; i<Ndays*24; i++) {
calc_longwave(&vlrad[i], vtskc[i], vtair[i], vhourly_vp[i]);
}
#if DBG
int dbg = 0;
for (; dbg < Ndays*24; dbg++) {
printf(“Lrad %d: %lf\n”, dbg, vlrad[dbg]);
}
#endif
return 1;
}
int ComputePrec(){
int i = 0;
int day = 0;
double tmp = 0;
FILE *out = fopen(“hourly_prec.txt”, “w”);
for (; i<Ndays*24; i++) {
day = i / 24;
tmp = Prec[day]/24;
fprintf(out, “%lf\n”, tmp);
#if DBG
printf(“Prec %d: %lf\n”, i, tmp);
#endif
}
fclose(out);
return 1;
}
int ComputeWind(){
int i = 0;
int day = 0;
double tmp = 0;
FILE *out = fopen(“hourly_wind.txt”, “w”);
for (; i<Ndays*24; i++) {
day = i / 24;
tmp = Wind[day];
fprintf(out, “%lf\n”, tmp);
#if DBG
printf(“WindSpeed %d: %lf\n”, i, tmp);
#endif
}
fclose(out);
return 1;
}
int ComputeRH(){
int i = 0;
int day = 0;
double tmp = 0;
FILE *out = fopen(“hourly_rh.txt”, “w”);
for (; i<Ndays*24; i++) {
day = i / 24;
tmp = RH[day];
fprintf(out, “%lf\n”, tmp);
#if DBG1
printf(“RH %d: %lf\n”, i, tmp);
#endif
}
fclose(out);
return 1;
}
void Output(char *fpath, double *data) {
FILE *out = fopen(fpath, “w”);
fprintf(out,”DIY\tYear\tMonth\tDay\tHour\tData\n”);
int i = 0;
char format[] = “%3d\t%4d\t%5d\t%3d\t%4d\t%lf\n”;
for (; i<Ndays*24; i++) {
fprintf(out, format,
vdmy[i].day_in_year, vdmy[i].year, vdmy[i].month, vdmy[i].day, vdmy[i].hour, data[i]);
}
fclose(out);
}
int SkipLines(FILE *f, int num, int linebuf){
char line[linebuf];
int n0 = num;
while (num–) {
if (!fgets(line, linebuf, f)) {
printf(“\n”);
err(“Error ocurrs when skipping lines remaining %d rounds\n”, num+1);
return 0;
}
}
#if DBG
printf(“In file no.%d: %d lines are skipped.\n”, fileno(f), n0);
#endif
return 1;
}
int Lines(char *fpath) {
FILE *f = fopen(fpath, “r”);
int recs = 0;
int buf;
int reserve = 0;
while ((buf = fgetc(f)) != EOF) {
reserve = buf;
if (buf == ‘\n’) {
++recs;
continue;
}
}
recs = ((reserve != ‘\r’ && reserve != ‘\n’) ? recs + 1 : recs);
fclose(f);
#if DBG
printf(“Total count of lines in file \”%s\” is %d\n”, fpath, recs);
#endif
return recs;
}
int MaxLineLen(char *fpath) {
FILE *f = fopen(fpath, “r”);
int len = 0;
int count = 0;
int c;
while ((c = fgetc(f)) != EOF) {
++count;
if (c == ‘\n’){
if (count > len)
 len = count;
count = 0;
continue;
}
}
len = count > len ? count : len;
fclose(f);
#if DBG
printf(“Max line lenth of file \”%s\” is %d\n”, fpath, len);
#endif
return len;
}
int DayInYear(int date, int mon, int year){
int count = 0;
int i = 1;
for (; i<=mon-1; i++) {
count += MonthDays(i, year);
}
count += date;
return count;
}
int YearDays(int year){
if (LeapYear(year))
 return 366;
return 365;
}
int MonthDays(int mon, int year){
if (mon > 12 || mon < 1)
return -1;
switch (mon) {
case 4:
case 6:
case 9:
case 11: {
return 30;
}
case 2: {
if (LeapYear(year))
 return 29;
return 28;
}
default:
return 31;
}
return -1;
}
int LeapYear(int year){
if ((year % 4 == 0 && year % 100 != 0) || year % 400 == 0)
 return 1;
return 0;
}
int getchar_ontime(){
#if LINUX
struct termios pret;
struct termios curt;
tcgetattr(STDIN_FILENO, &pret);
curt = pret;
curt.c_lflag &= !(ICANON | ECHO);
tcsetattr(STDIN_FILENO, TCSANOW, &curt);
int c = getchar();
tcsetattr(STDIN_FILENO, TCSANOW, &pret);
return c;
#else
system(“pause”);
return 1;
#endif
}
void err(char *format, …){
va_list ap;
va_start(ap, format);
vfprintf(stderr, format, ap);
va_end(ap);
}

Leave a Reply

Your email address will not be published.

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

Post Navigation