Merge pull request #4 from BramDr/develop

Merge develop into master

vic/plugins/crops/include/crop.h

@@ -13,13 +13,17 @@ typedef struct {
 } crop_con_map_struct;
 
 typedef struct {
+    double tsum1;
+    double tsum2;
     double *DVS_point;
     double *N_amount;
     double *P_amount;
     double *K_amount;
 } crop_con_struct;
 
 typedef struct {
+    double *tsum1;
+    double *tsum2;
     double **DVS_point;
     double **N_amount;
     double **P_amount;

vic/plugins/crops/src/crop_alloc_free.c

@@ -51,6 +51,16 @@ crop_alloc(void)
         crop_con[i] = malloc(crop_con_map[i].nc_active * sizeof(*crop_con[i]));
         check_alloc_status(crop_con[i], "Memory allocation error");
 
+        crop_force[i].tsum1 =
+            malloc(crop_con_map[i].nc_active *
+                   sizeof(*crop_force[i].tsum1));
+        check_alloc_status(crop_force[i].tsum1, "Memory allocation error");
+
+        crop_force[i].tsum2 =
+            malloc(crop_con_map[i].nc_active *
+                   sizeof(*crop_force[i].tsum2));
+        check_alloc_status(crop_force[i].tsum2, "Memory allocation error");
+
         crop_force[i].DVS_point =
             malloc(crop_con_map[i].nc_active *
                    sizeof(*crop_force[i].DVS_point));
@@ -153,6 +163,8 @@ crop_finalize(void)
         free(crop_con_map[i].veg_class);
         free(crop_con_map[i].Cc);
         free(crop_con[i]);
+        free(crop_force[i].tsum1);
+        free(crop_force[i].tsum2);
         free(crop_force[i].DVS_point);
         free(crop_force[i].N_amount);
         free(crop_force[i].P_amount);

vic/plugins/crops/src/crop_forcing.c

@@ -14,6 +14,8 @@ crop_forcing(void)
 
     double                           *dvar;
 
+    size_t                            d4count[4];
+    size_t                            d4start[4];
     size_t                            d5count[5];
     size_t                            d5start[5];
 
@@ -25,6 +27,15 @@ crop_forcing(void)
     dvar = malloc(local_domain.ncells_active * sizeof(*dvar));
     check_alloc_status(dvar, "Memory allocation error.");
 
+    d4start[0] = 0;
+    d4start[1] = 0;
+    d4start[2] = 0;
+    d4start[3] = 0;
+    d4count[0] = 1;
+    d4count[1] = 1;
+    d4count[2] = global_domain.n_ny;
+    d4count[3] = global_domain.n_nx;
+
     d5start[0] = 0;
     d5start[1] = 0;
     d5start[2] = 0;
@@ -36,6 +47,51 @@ crop_forcing(void)
     d5count[3] = global_domain.n_ny;
     d5count[4] = global_domain.n_nx;
 
+    // Get tsum forcing data
+    if (plugin_options.WOFOST_FORCE_TSUM) {
+        // Get tsum1
+        if (plugin_global_param.forcerun[FORCING_TSUM_1]) {
+            d4start[0] = plugin_global_param.forceskip[FORCING_TSUM_1] +
+                         plugin_global_param.forceoffset[FORCING_TSUM_1];
+            for (j = 0; j < plugin_options.NCROPTYPES; j++) {
+                d4start[1] = j;
+                get_scatter_nc_field_double(&(plugin_filenames.forcing[
+                                                  FORCING_TSUM_1]),
+                                            plugin_filenames.f_varname[
+                                                FORCING_TSUM_1], d4start,
+                                            d4count, dvar);
+
+                for (i = 0; i < local_domain.ncells_active; i++) {
+                    iCrop = crop_con_map[i].cidx[j];
+                    if (iCrop != NODATA_VEG) {
+                        crop_force[i].tsum1[iCrop] = dvar[i];
+                    }
+                }
+            }
+        }
+
+        // Get tsum2
+        if (plugin_global_param.forcerun[FORCING_TSUM_2]) {
+            d4start[0] = plugin_global_param.forceskip[FORCING_TSUM_2] +
+                         plugin_global_param.forceoffset[FORCING_TSUM_2];
+            for (j = 0; j < plugin_options.NCROPTYPES; j++) {
+                d4start[1] = j;
+                get_scatter_nc_field_double(&(plugin_filenames.forcing[
+                                                  FORCING_TSUM_2]),
+                                            plugin_filenames.f_varname[
+                                                FORCING_TSUM_2], d4start,
+                                            d4count, dvar);
+
+                for (i = 0; i < local_domain.ncells_active; i++) {
+                    iCrop = crop_con_map[i].cidx[j];
+                    if (iCrop != NODATA_VEG) {
+                        crop_force[i].tsum2[iCrop] = dvar[i];
+                    }
+                }
+            }
+        }
+    }
+
     // Get fertilizer forcing data
     if (plugin_options.WOFOST_FORCE_FERT) {
         // Get DVS_point

vic/plugins/crops/src/crop_get_global_param.c

@@ -33,6 +33,10 @@ crop_get_global_param(char *cmdstr)
         sscanf(cmdstr, "%*s %s", flgstr);
         plugin_options.WOFOST_PFERT = str_to_bool(flgstr);
     }
+    else if (strcasecmp("WOFOST_POTENTIAL_TEMPERATURE", optstr) == 0) {
+        sscanf(cmdstr, "%*s %s", flgstr);
+        plugin_options.WOFOST_PTEMP = str_to_bool(flgstr);
+    }
     else if (strcasecmp("WOFOST_DIST_SEASON", optstr) == 0) {
         sscanf(cmdstr, "%*s %s", flgstr);
         plugin_options.WOFOST_DIST_SEASON = str_to_bool(flgstr);
@@ -41,6 +45,10 @@ crop_get_global_param(char *cmdstr)
         sscanf(cmdstr, "%*s %s", flgstr);
         plugin_options.WOFOST_DIST_TSUM = str_to_bool(flgstr);
     }
+    else if (strcasecmp("WOFOST_FORCE_TSUM", optstr) == 0) {
+        sscanf(cmdstr, "%*s %s", flgstr);
+        plugin_options.WOFOST_FORCE_TSUM = str_to_bool(flgstr);
+    }
     else if (strcasecmp("WOFOST_DIST_FERTILIZER", optstr) == 0) {
         sscanf(cmdstr, "%*s %s", flgstr);
         plugin_options.WOFOST_DIST_FERT = str_to_bool(flgstr);
@@ -57,6 +65,10 @@ crop_get_global_param(char *cmdstr)
         sscanf(cmdstr, "%*s %s", flgstr);
         plugin_options.WOFOST_CONTINUE = str_to_bool(flgstr);
     }
+    else if (strcasecmp("WOFOST_FORCE_TSUM", optstr) == 0) {
+        sscanf(cmdstr, "%*s %s", flgstr);
+        plugin_options.WOFOST_FORCE_TSUM = str_to_bool(flgstr);
+    }
     else if (strcasecmp("WOFOST_FORCE_FERTILIZER", optstr) == 0) {
         sscanf(cmdstr, "%*s %s", flgstr);
         plugin_options.WOFOST_FORCE_FERT = str_to_bool(flgstr);
@@ -103,6 +115,19 @@ crop_validate_global_param(void)
     if (!plugin_options.FORCE_CO2) {
         log_err("WOFOST = TRUE but FORCE_CO2 = FALSE");
     }
+    if (plugin_options.WOFOST_FORCE_TSUM) {
+        if (!plugin_options.WOFOST_DIST_TSUM) {
+            log_err("WOFOST_FORCE_TSUM = TRUE but WOFOST_DIST_TSUM = FALSE");
+        }
+        if (strcasecmp(plugin_filenames.f_path_pfx[FORCING_TSUM_1],
+                       MISSING_S) == 0) {
+            log_err("WOFOST_FORCE_TSUM = TRUE but tsum1 forcing file is missing");
+        }
+        if (strcasecmp(plugin_filenames.f_path_pfx[FORCING_TSUM_2],
+                       MISSING_S) == 0) {
+            log_err("WOFOST_FORCE_TSUM = TRUE but tsum2 forcing file is missing");
+        }
+    }
     if (plugin_options.WOFOST_FORCE_FERT) {
         if (!plugin_options.WOFOST_DIST_FERT) {
             log_err("WOFOST_FORCE_FERT = TRUE but WOFOST_DIST_FERT = FALSE");

vic/plugins/crops/src/crop_init_library.c

@@ -37,6 +37,8 @@ initialize_crop_con(crop_con_struct *crop_con)
 
     size_t                      i;
 
+    crop_con->tsum1 = 0.;
+    crop_con->tsum2 = 0.;
     for (i = 0; i < plugin_options.NFERTTIMES; i++) {
         crop_con->DVS_point[i] = 0.;
         crop_con->N_amount[i] = 0.;
@@ -58,6 +60,8 @@ initialize_crop_force(crop_force_struct *crop_force,
     size_t                      iFert;
 
     for (iCrop = 0; iCrop < ncrops; iCrop++) {
+        crop_force->tsum1[iCrop] = 0.;
+        crop_force->tsum2[iCrop] = 0.;
         for (iFert = 0; iFert < plugin_options.NFERTTIMES; iFert++) {
             crop_force->DVS_point[iCrop][iFert] = 0.;
             crop_force->N_amount[iCrop][iFert] = 0.;

vic/plugins/crops/src/crop_output.c

@@ -944,7 +944,7 @@ crop_put_state_data(size_t iCell)
                     cgrid->growing;
                 out_data[iCell][N_OUTVAR_TYPES +
                                 OUT_CROP_TMAXSTRESS][crop_class] +=
-                    cgrid->met->TmaxStress *
+                    cgrid->met->TempStress *
                     area_fract *
                     cgrid->growing;
                 out_data[iCell][N_OUTVAR_TYPES +
@@ -1255,7 +1255,7 @@ crop_put_rate_data(size_t iCell)
                     cgrid->growing;
                 out_data[iCell][N_OUTVAR_TYPES +
                                 OUT_CROP_TMAXSTRESS][crop_class] +=
-                    cgrid->met->TmaxStress *
+                    cgrid->met->TempStress *
                     area_fract *
                     cgrid->growing;
                 out_data[iCell][N_OUTVAR_TYPES +

vic/plugins/crops/src/crop_run.c

@@ -234,6 +234,11 @@ crop_run(size_t iCell)
                     /* Determine if the sowing already has occurred */
                     IfSowing(cgrid, &cgrid->start);
 
+                    if (plugin_options.WOFOST_DIST_TSUM) {
+                        cgrid->crp->prm.TempSum1 = crop_con[iCell][iCrop].tsum1;
+                        cgrid->crp->prm.TempSum2 = crop_con[iCell][iCrop].tsum2;
+                    }
+
                     if (plugin_options.WOFOST_DIST_FERT) {
                         for (iTime = 0;
                              iTime < plugin_options.NFERTTIMES;

vic/plugins/crops/src/crop_start.c

@@ -37,6 +37,9 @@ crop_start(void)
     extern domain_struct           global_domain;
     extern plugin_filenames_struct plugin_filenames;
 
+    size_t                         tsum_idx[2] = {
+        FORCING_TSUM_1, FORCING_TSUM_2
+    };
     size_t                         fert_idx[4] = {
         FORCING_FERT_DVS, FORCING_FERT_N, FORCING_FERT_P, FORCING_FERT_K
     };
@@ -61,6 +64,33 @@ crop_start(void)
         plugin_options.NFERTTIMES = get_nc_dimension(&(plugin_filenames.crop),
                                                      "fertilizer_times");
     }
+    if (plugin_options.WOFOST_FORCE_TSUM) {
+        for (i = 0; i < 2; i++) {
+            // Get information from the forcing file(s)
+            // Open first-year forcing files and get info
+            snprintf(plugin_filenames.forcing[tsum_idx[i]].nc_filename,
+                     MAXSTRING, "%s%4d.nc",
+                     plugin_filenames.f_path_pfx[tsum_idx[i]],
+                     global_param.startyear);
+            status = nc_open(plugin_filenames.forcing[tsum_idx[i]].nc_filename,
+                             NC_NOWRITE,
+                             &(plugin_filenames.forcing[tsum_idx[i]].nc_id));
+            check_nc_status(status, "Error opening %s",
+                            plugin_filenames.forcing[tsum_idx[i]].nc_filename);
+
+            dim_len = get_nc_dimension(&(plugin_filenames.forcing[tsum_idx[i]]),
+                                       "crop_class");
+            if (dim_len != plugin_options.NCROPTYPES) {
+                log_err(
+                    "Tsum forcing crop_class length is not equal to NCROPTYPES");
+            }
+
+            // Close first-year forcing files
+            status = nc_close(plugin_filenames.forcing[tsum_idx[i]].nc_id);
+            check_nc_status(status, "Error closing %s",
+                            plugin_filenames.forcing[tsum_idx[i]].nc_filename);
+        }
+    }
     if (plugin_options.WOFOST_FORCE_FERT) {
         for (i = 0; i < 4; i++) {
             // Get information from the forcing file(s)

vic/plugins/crops/src/wofost_copy_library.c

@@ -307,6 +307,7 @@ copy_wofost_parameters(Parameters *from,
     to->TempEffMax = from->TempEffMax;
     to->TSumEmergence = from->TSumEmergence;
     to->IdentifyAnthesis = from->IdentifyAnthesis;
+    to->MaxOptimumTemp = from->MaxOptimumTemp;
     to->OptimumDaylength = from->OptimumDaylength;
     to->CriticalDaylength = from->CriticalDaylength;
     to->SatVernRequirement = from->SatVernRequirement;
@@ -493,7 +494,7 @@ copy_wofost_vic(Meteo *from,
     to->CosLD = from->CosLD;
     to->DiffRadPP = from->DiffRadPP;
     to->DSinBE = from->DSinBE;
-    to->TmaxStress = from->TmaxStress;
+    to->TempStress = from->TempStress;
     to->TminStress = from->TminStress;
 }
 

vic/plugins/crops/src/wofost_init.c

@@ -205,14 +205,16 @@ wofost_set_data(void)
             for (i = 0; i < local_domain.ncells_active; i++) {
                 iCrop = crop_con_map[i].cidx[j];
                 if (iCrop != NODATA_VEG) {
-                    for (k = 0; k < options.SNOW_BAND; k++) {
-                        iGrid = Grid[i][k];
-                        for (l = 0; l < (size_t)iCrop; l++) {
-                            iGrid = iGrid->next;
-                        }
-
-                        iGrid->crp->prm.TempSum1 = dvar[i];
-                    }
+                    crop_con[i][iCrop].tsum1 = dvar[i];
+
+//                    for (k = 0; k < options.SNOW_BAND; k++) {
+//                        iGrid = Grid[i][k];
+//                        for (l = 0; l < (size_t)iCrop; l++) {
+//                            iGrid = iGrid->next;
+//                        }
+//
+//                        iGrid->crp->prm.TempSum1 = dvar[i];
+//                    }
                 }
             }
         }
@@ -227,14 +229,16 @@ wofost_set_data(void)
             for (i = 0; i < local_domain.ncells_active; i++) {
                 iCrop = crop_con_map[i].cidx[j];
                 if (iCrop != NODATA_VEG) {
-                    for (k = 0; k < options.SNOW_BAND; k++) {
-                        iGrid = Grid[i][k];
-                        for (l = 0; l < (size_t)iCrop; l++) {
-                            iGrid = iGrid->next;
-                        }
-
-                        iGrid->crp->prm.TempSum2 = dvar[i];
-                    }
+                    crop_con[i][iCrop].tsum2 = dvar[i];
+
+//                    for (k = 0; k < options.SNOW_BAND; k++) {
+//                        iGrid = Grid[i][k];
+//                        for (l = 0; l < (size_t)iCrop; l++) {
+//                            iGrid = iGrid->next;
+//                        }
+//
+//                        iGrid->crp->prm.TempSum2 = dvar[i];
+//                    }
                 }
             }
         }

vic/plugins/crops/src/wofost_init_library.c

@@ -202,6 +202,7 @@ initialize_wofost_parameters(Parameters *prm)
     prm->TempEffMax = 0.;
     prm->TSumEmergence = 0.;
     prm->IdentifyAnthesis = 0;
+    prm->MaxOptimumTemp = 0;
     prm->OptimumDaylength = 0.;
     prm->CriticalDaylength = 0.;
     prm->SatVernRequirement = 0.;
@@ -382,7 +383,7 @@ initialize_wofost_vic(Meteo *met)
     met->CosLD = 0.;
     met->DiffRadPP = 0.;
     met->DSinBE = 0.;
-    met->TmaxStress = 0.;
+    met->TempStress = 0.;
     met->TminStress = 0.;
 }