Code help; loop problem
Hi,
I have the following program. It worked great until I added the For loop at the begining, and modified if((i=0)&&(doy <=1)). Now the output looks good for a couple iterations, but then goes into this endless output of day nan nan nan ~0.035 nan nan nan nan. Any suggestions on how to fix this problem? Many thanks!
The code:
Code:#define SMA 0.00036 /* DIM scaling parameter for SM curve */ #define SMB 105.0 /* DIM scaling parameter for SM curve */ #define FMET 0.59 /* DIM metabolic fraction of NPP; Potter et al. 1993, CHANGE FOR DIFFERENT BIOMES */ #define FSTR 0.70 /* DIM slow SOC fraction of structural SOC pool; Ise and Moorcroft 2006 */ #define KSTR 0.40 /* DIM proportional rate K for structural SOC decomposition; Ise and Moorcroft 2006 */ #define KSLW 0.0093 /* DIM proportional rate K for slow SOC decomposition; Ise and Moorcroft 2006 */ #define FRAUT 0.53 /* DIM proportion of GPP represented by autotrophic respiration. This value is taken from Choudhury (2001). The origional model used Gifford et al. 2005, ** /* value of 0.46. CHANGE FOR DIFFERENT BIOMES */ #define ME 0.45 /* DIM Microbial efficiency; McGill et al. 1981; Potter et al. 1993 */ #define KOPT 0.0301 /* Optimal decomposition rate (1/d) for soil metabolic fraction; Ise and Moorcroft 2006 */ #include <math.h> #include <stdio.h> #include <stdlib.h> #include <malloc.h> #include <unistd.h> int main() { int doy; doy=0; double soil_t; /* surface soil temperature; deg K */ double soil_m; /* surface soil moisture; % saturation */ double gpp; /* GPP; g C m-2 d-1 */ double npp; double rmet; double rstr; double rslw; double npp_ann_init; /* Initial NPP pool; g C m-2 yr-1 */ double soil_t_av; /* long-term avg surface soil T; deg K */ double soil_m_av; /* long-term avg JJA soil moisture; % saturation */ double met_c_init; /* initial soil metabolic C pool; g C m-2; Ise and Moorcroft, 2006 */ double struct_c_init; /* initial structural soil C pool; g C m-2; Ise and Moorcroft, 2006 */ double slow_c_init; /* initial slow soil C pool; g C m-2; Ise and Moorcroft, 2006 */ double r_aut; /* autotrophic respiration; g C m-2 d-1 */ double r_het; /* heterotrophic respiration; g C m-2 d-1 */ double dCmet; double dCstr; double dCslw; double cmet; double cstr; double cslw; double r_tot; /* total respiration; g C m-2 d-1 */ double k_met_opt; /* optimal soil metabolic C decomposition rate; d-1 */ double k_met_av; /* average soil metabolic C decomposition rate; d-1 */ double k_met; /* actual soil metabolic C decomp rate; d-1 */ double k_struct; /* soil structural C decomp rate; d-1 */ double k_slow; /* soil slow C decomp rate; d-1 */ double tmult; /* DIM 0-1, soil temperature rate response curve */ double tmult_av; /* DIM 0-1, soil temperature rate response under avg conditions */ double wmult; /* DIM 0-1, soil moisture rate response curve */ double wmult_av; /* DIM 0-1, soil moisture rate response under avg conditions */ int flag; /* DIM 0.0 or 1.0 for spinup or daily run, respectively */ /* 1.0 open data input file */ FILE *datain; int i; i=1; for (i = 0 ; i < 2; i++) { /* datain =fopen("/home/jennifer.watts/Carbon_Models/TCF_model_modified/Testing_Spinup/", "r"); */ datain =fopen("tcf_spinup.txt", "r"); if (datain == NULL) { fprintf(stderr, "Can't open input file!\n"); exit(-1); } while(fscanf(datain,"%d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf", &flag, &doy, &gpp, &soil_t, &soil_m, &soil_t_av, &soil_m_av, &npp_ann_init, &met_c_init, &struct_c_init, &slow_c_init) != EOF) { doy++; /* address inital values for year 1, day 1, of spin-up*/ if((i=0) && (doy <= 1 )) { cmet = met_c_init; cstr = struct_c_init; cslw = slow_c_init; rmet = 0; rstr = 0; rslw = 0; } /* 3.0 Environmental Constraint Calculations*/ k_met_opt = KOPT; /* calculate soil temperature rate response (Lloyd and Taylor 1994) */ tmult = exp(308.56 * ((1.0 / 66.02) - (1.0 / ((soil_t) - 227.13)))); /* calculate soil moisture response (Monsi 1969, Bunnell and Tait 1974, Oberbauer et al. 1991, 1996) */ wmult = SMA * (SMB * soil_m - pow(soil_m, 2.0)); /* 4.0 Calculate soil CO2 input and output; soil C pool values */ /* calculate decomposition rate constants */ k_met = k_met_opt * tmult * wmult; k_struct = KSTR * k_met; k_slow = KSLW * k_met; rmet = (k_met * cmet * (1.0 - ME)); rstr = (k_struct * cstr * (1.0 - ME)); rslw = (k_slow * cslw * (1.0 - ME)); r_het = (rmet + rstr + rslw); r_aut = (gpp * FRAUT); r_tot = (r_het + r_aut); npp = (gpp - r_aut); dCmet = (FMET * npp - rmet); dCstr = (((1 - FMET) * npp) - rstr - FSTR * cstr); dCslw = ((FSTR * cstr) - rslw); cmet += (cmet + dCmet); cstr += (cstr + dCstr); cslw += (cslw + dCslw); /* print out daily updating of carbon pools for spinup run */ printf("%d %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f\n", doy, cmet, cstr, cslw, npp, r_tot, dCmet, dCstr, dCslw); } fclose(datain); } }
