UW-Hydro · bartnijssen · Sep 13, 2018 · Sep 10, 2018 · Sep 10, 2018 · Sep 10, 2018
diff --git a/vic/drivers/cesm/src/display_current_settings.c b/vic/drivers/cesm/src/display_current_settings.c
@@ -283,6 +283,12 @@ display_current_settings(int mode)
     else {
         fprintf(LOG_DEST, "ORGANIC_FRACT\t\tFALSE\n");
     }
+    if (options.BULK_DENSITY_COMB) {
+        fprintf(LOG_DEST, "BULK_DENSITY_COMB\t\tTRUE\n");
+    }
+    else {
+        fprintf(LOG_DEST, "BULK_DENSITY_COMB\t\tFALSE\n");
+    }
 
     fprintf(LOG_DEST, "\n");
     if (options.VEGLIB_PHOTO) {

diff --git a/vic/drivers/cesm/src/get_global_param.c b/vic/drivers/cesm/src/get_global_param.c
@@ -338,6 +338,10 @@ get_global_param(FILE *gp)
                 sscanf(cmdstr, "%*s %s", flgstr);
                 options.ORGANIC_FRACT = str_to_bool(flgstr);
             }
+            else if (strcasecmp("BULK_DENSITY_COMB", optstr) == 0) {
+                sscanf(cmdstr, "%*s %s", flgstr);
+                options.BULK_DENSITY_COMB = str_to_bool(flgstr);
+            }
             else if (strcasecmp("VEGLIB_PHOTO", optstr) == 0) {
                 sscanf(cmdstr, "%*s %s", flgstr);
                 options.VEGLIB_PHOTO = str_to_bool(flgstr);

diff --git a/vic/drivers/classic/src/display_current_settings.c b/vic/drivers/classic/src/display_current_settings.c
@@ -295,6 +295,13 @@ display_current_settings(int mode)
     else {
         fprintf(LOG_DEST, "ORGANIC_FRACT\t\tFALSE\n");
     }
+    if (options.BULK_DENSITY_COMB) {
+        fprintf(LOG_DEST, "BULK_DENSITY_COMB\t\tTRUE\n");
+    }
+    else {
+        fprintf(LOG_DEST, "BULK_DENSITY_COMB\t\tFALSE\n");
+    }
+
 
     fprintf(LOG_DEST, "\n");
     fprintf(LOG_DEST, "Input Veg Data:\n");

diff --git a/vic/drivers/classic/src/get_global_param.c b/vic/drivers/classic/src/get_global_param.c
@@ -446,6 +446,10 @@ get_global_param(FILE *gp)
                 sscanf(cmdstr, "%*s %s", flgstr);
                 options.ORGANIC_FRACT = str_to_bool(flgstr);
             }
+            else if (strcasecmp("BULK_DENSITY_COMB", optstr) == 0) {
+                sscanf(cmdstr, "%*s %s", flgstr);
+                options.BULK_DENSITY_COMB = str_to_bool(flgstr);
+            }
             else if (strcasecmp("VEGLIB", optstr) == 0) {
                 sscanf(cmdstr, "%*s %s", filenames.veglib);
             }

diff --git a/vic/drivers/classic/src/read_soilparam.c b/vic/drivers/classic/src/read_soilparam.c
@@ -438,6 +438,29 @@ read_soilparam(FILE            *soilparam,
             }
         }
 
+        if (options.BULK_DENSITY_COMB) {
+            /* read soil bulk density */
+            for (layer = 0; layer < options.Nlayer; layer++) {
+                token = strtok(NULL, delimiters);
+                while (token != NULL && (length = strlen(token)) == 0) { 
+                    token = strtok(NULL, delimiters);
+                }    
+                if (token == NULL) {
+                    log_err("Can't find values for SOIL BULK DENSITY for "
+                            "layer %zu in soil file", layer);
+                }    
+                sscanf(token, "%lf", &(temp->bulk_density)[layer]);
+            }    
+        }    
+        else {
+            for (layer = 0; layer < options.Nlayer; layer++) {
+                temp->bulk_density[layer] =
+                    (1 - 
+                     temp->organic[layer]) * temp->bulk_dens_min[layer] +
+                    temp->organic[layer] * temp->bulk_dens_org[layer];
+            }    
+        }
+
         /* read cell gmt offset */
         token = strtok(NULL, delimiters);
         while (token != NULL && (length = strlen(token)) == 0) {
@@ -593,10 +616,6 @@ read_soilparam(FILE            *soilparam,
            Compute Soil Layer Properties
         *******************************************/
         for (layer = 0; layer < options.Nlayer; layer++) {
-            temp->bulk_density[layer] =
-                (1 -
-                 temp->organic[layer]) * temp->bulk_dens_min[layer] +
-                temp->organic[layer] * temp->bulk_dens_org[layer];
             temp->soil_density[layer] =
                 (1 -
                  temp->organic[layer]) * temp->soil_dens_min[layer] +

diff --git a/vic/drivers/image/src/display_current_settings.c b/vic/drivers/image/src/display_current_settings.c
@@ -292,6 +292,12 @@ display_current_settings(int mode)
     else {
         fprintf(LOG_DEST, "ORGANIC_FRACT\t\tFALSE\n");
     }
+    if (options.BULK_DENSITY_COMB) {
+        fprintf(LOG_DEST, "BULK_DENSITY_COMB\t\tTRUE\n");
+    }
+    else {
+        fprintf(LOG_DEST, "BULK_DENSITY_COMB\t\tFALSE\n");
+    }
 
     fprintf(LOG_DEST, "\n");
     if (options.VEGLIB_PHOTO) {

diff --git a/vic/drivers/image/src/get_global_param.c b/vic/drivers/image/src/get_global_param.c
@@ -429,6 +429,10 @@ get_global_param(FILE *gp)
                 sscanf(cmdstr, "%*s %s", flgstr);
                 options.ORGANIC_FRACT = str_to_bool(flgstr);
             }
+            else if (strcasecmp("BULK_DENSITY_COMB", optstr) == 0) {
+                sscanf(cmdstr, "%*s %s", flgstr);
+                options.BULK_DENSITY_COMB = str_to_bool(flgstr);
+            }
             else if (strcasecmp("VEGLIB_PHOTO", optstr) == 0) {
                 sscanf(cmdstr, "%*s %s", flgstr);
                 options.VEGLIB_PHOTO = str_to_bool(flgstr);

diff --git a/vic/drivers/python/src/display_current_settings.c b/vic/drivers/python/src/display_current_settings.c
@@ -257,6 +257,13 @@ display_current_settings(int mode)
     else {
         fprintf(LOG_DEST, "ORGANIC_FRACT\t\tFALSE\n");
     }
+    if (options.BULK_DENSITY_COMB) {
+        fprintf(LOG_DEST, "BULK_DENSITY_COMB\t\tTRUE\n");
+    }
+    else {
+        fprintf(LOG_DEST, "BULK_DENSITY_COMB\t\tFALSE\n");
+    }
+
 
     fprintf(LOG_DEST, "\n");
     if (options.VEGLIB_PHOTO) {

diff --git a/vic/drivers/shared_all/src/initialize_options.c b/vic/drivers/shared_all/src/initialize_options.c
@@ -84,6 +84,7 @@ initialize_options()
     options.JULY_TAVG_SUPPLIED = false;
     options.LAI_SRC = FROM_VEGLIB;
     options.ORGANIC_FRACT = false;
+    options.BULK_DENSITY_COMB = false;
     options.VEGLIB_FCAN = false;
     options.VEGLIB_PHOTO = false;
     options.VEGPARAM_ALB = false;

diff --git a/vic/drivers/shared_all/src/print_library_shared.c b/vic/drivers/shared_all/src/print_library_shared.c
@@ -516,6 +516,8 @@ print_option(option_struct *option)
             option->LAKE_PROFILE ? "true" : "false");
     fprintf(LOG_DEST, "\tORGANIC_FRACT        : %s\n",
             option->ORGANIC_FRACT ? "true" : "false");
+    fprintf(LOG_DEST, "\tBULK_DENSITY_COMB        : %s\n",
+            option->BULK_DENSITY_COMB ? "true" : "false");
     fprintf(LOG_DEST, "\tSTATE_FORMAT         : %d\n", option->STATE_FORMAT);
     fprintf(LOG_DEST, "\tINIT_STATE           : %s\n",
             option->INIT_STATE ? "true" : "false");

diff --git a/vic/drivers/shared_image/src/get_global_domain.c b/vic/drivers/shared_image/src/get_global_domain.c
@@ -36,7 +36,7 @@ get_global_domain(nameid_struct *domain_nc_nameid,
 {
     int    *run = NULL;
     int    *mask = NULL;
-    int typeid;
+    int     typeid;
     double *var = NULL;
     size_t  i;
     size_t  j;

diff --git a/vic/drivers/shared_image/src/vic_init.c b/vic/drivers/shared_image/src/vic_init.c
@@ -609,6 +609,29 @@ vic_init(void)
             }
         }
     }
+    // bulk density of mineral and organic soil
+    if (options.BULK_DENSITY_COMB) {
+        for (j = 0; j < options.Nlayer; j++) {
+            d3start[0] = j;
+            get_scatter_nc_field_double(&(filenames.params),
+                                        "bulk_density_comb",
+                                        d3start, d3count, dvar);
+            for (i = 0; i < local_domain.ncells_active; i++) {
+                soil_con[i].bulk_density[j] = (double) dvar[i];
+            }
+        }
+    }
+    else {
+        for (i = 0; i < local_domain.ncells_active; i++) {
+            for (j = 0; j < options.Nlayer; j++) {
+                soil_con[i].bulk_density[j] =
+                    (1 - soil_con[i].organic[j]) *
+                    soil_con[i].bulk_dens_min[j] +
+                    soil_con[i].organic[j] * soil_con[i].bulk_dens_org[j];
+            }
+        }
+    }
+
 
     // Wcr: critical point for each layer
     // Note this value is  multiplied with the maximum moisture in each layer
@@ -714,9 +737,6 @@ vic_init(void)
     for (i = 0; i < local_domain.ncells_active; i++) {
         for (j = 0; j < options.Nlayer; j++) {
             // compute layer properties
-            soil_con[i].bulk_density[j] =
-                (1 - soil_con[i].organic[j]) * soil_con[i].bulk_dens_min[j] +
-                soil_con[i].organic[j] * soil_con[i].bulk_dens_org[j];
             soil_con[i].soil_density[j] =
                 (1 - soil_con[i].organic[j]) * soil_con[i].soil_dens_min[j] +
                 soil_con[i].organic[j] * soil_con[i].soil_dens_org[j];

diff --git a/vic/drivers/shared_image/src/vic_mpi_support.c b/vic/drivers/shared_image/src/vic_mpi_support.c
@@ -488,7 +488,7 @@ create_MPI_option_struct_type(MPI_Datatype *mpi_type)
     MPI_Datatype   *mpi_types;
 
     // nitems has to equal the number of elements in option_struct
-    nitems = 53;
+    nitems = 54;
     blocklengths = malloc(nitems * sizeof(*blocklengths));
     check_alloc_status(blocklengths, "Memory allocation error.");
 
@@ -706,6 +706,10 @@ create_MPI_option_struct_type(MPI_Datatype *mpi_type)
     offsets[i] = offsetof(option_struct, ORGANIC_FRACT);
     mpi_types[i++] = MPI_C_BOOL;
 
+    // book BULK_DENSITY_COMB;
+    offsets[i] = offsetof(option_struct, BULK_DENSITY_COMB);
+    mpi_types[i++] = MPI_C_BOOL;
+
     // unsigned short STATE_FORMAT;
     offsets[i] = offsetof(option_struct, STATE_FORMAT);
     mpi_types[i++] = MPI_UNSIGNED_SHORT;

diff --git a/vic/vic_run/include/vic_def.h b/vic/vic_run/include/vic_def.h
@@ -279,6 +279,7 @@ typedef struct {
                                           FROM_VEGPARAM = use LAI values from the veg param file */
     bool LAKE_PROFILE;   /**< TRUE = user-specified lake/area profile */
     bool ORGANIC_FRACT;  /**< TRUE = organic matter fraction of each layer is read from the soil parameter file; otherwise set to 0.0. */
+    bool BULK_DENSITY_COMB; /**< TRUE = soil bulk density (combined mineral and organic matter) read from soil parameter file; otherwise set to 0.0 */
 
     // state options
     unsigned short int STATE_FORMAT;  /**< TRUE = model state file is binary (default) */