diff --git a/applications/sopflow_main.cpp b/applications/sopflow_main.cpp
index df98dcc7..b3240000 100644
--- a/applications/sopflow_main.cpp
+++ b/applications/sopflow_main.cpp
@@ -10,9 +10,10 @@ int main(int argc, char **argv) {
   SOPFLOW sopflow;
   char file[PETSC_MAX_PATH_LEN];
   char scenfile[PETSC_MAX_PATH_LEN];
+  char loadfile[PETSC_MAX_PATH_LEN];
   char outputdir[PETSC_MAX_PATH_LEN];
   PetscBool outputdir_set;
-  PetscBool flg = PETSC_FALSE, flgscen = PETSC_FALSE;
+  PetscBool flg = PETSC_FALSE, flgscen = PETSC_FALSE, flgload = PETSC_FALSE;
   PetscBool print_output = PETSC_FALSE, save_output = PETSC_FALSE;
   PetscLogStage stages[3];
   char appname[] = "sopflow";
@@ -54,6 +55,11 @@ int main(int argc, char **argv) {
                                PETSC_MAX_PATH_LEN, &flgscen);
   CHKERRQ(ierr);
 
+  /* Get load data file from command line */
+  ierr = PetscOptionsGetString(NULL, NULL, "-loadfile", loadfile,
+                               PETSC_MAX_PATH_LEN, &flgload);
+  CHKERRQ(ierr);
+
   /* Stage 1 - Application creation and reading data */
   ierr = PetscLogStagePush(stages[0]);
   CHKERRQ(ierr);
@@ -78,6 +84,13 @@ int main(int argc, char **argv) {
     CHKERRQ(ierr);
   }
 
+  /* Set Load Data file */
+  if (flgload) {
+    ierr = SOPFLOWSetScenarioData(sopflow, SOPFLOW_NATIVE_SINGLEPERIOD, LOAD,
+                                  loadfile);
+    CHKERRQ(ierr);
+  }
+
   /* Set a subset of scenarios to be selected. Can use the option -sopflow_Ns
    * instead */
   /*   ierr = SOPFLOWSetNumScenarios(sopflow,2);CHKERRQ(ierr); */
diff --git a/datafiles/case9/10_scenarios_9bus_load.txt b/datafiles/case9/10_scenarios_9bus_load.txt
new file mode 100644
index 00000000..a9374ae4
--- /dev/null
+++ b/datafiles/case9/10_scenarios_9bus_load.txt
@@ -0,0 +1,11 @@
+scenario_nr,weight,kv_list
+1,0.101,8_Load_1_185.0_1000,6_Load_1_105.5_2000
+2,0.102,6_Load_1_105.0_10
+3,0.103,8_Load_1_75.0_20
+4,0.104,6_Load_1_95.0_30
+5,0.105,8_Load_1_115.0_40
+6,0.106,6_Load_1_65.0_50
+7,0.107,8_Load_1_33.0_60
+8,0.108,6_Load_1_125.0_70
+9,0.109,8_Load_1_55.0_80
+10,0.1110,6_Load_1_250.0_90
diff --git a/include/common.h b/include/common.h
index 357a2e58..b906419f 100644
--- a/include/common.h
+++ b/include/common.h
@@ -21,7 +21,7 @@ typedef enum {
 typedef enum {
   FORECAST_WIND = 1,
   FORECAST_LOAD_P = 2,
-  FORECAST_LOAD_Q = 3
+  FORECAST_LOAD_PQ = 3
 } ForecastType;
 
 /**
diff --git a/include/opflow.h b/include/opflow.h
index bce416c8..f46bc3a1 100644
--- a/include/opflow.h
+++ b/include/opflow.h
@@ -65,6 +65,7 @@ namespace OPFLOWOptions {
 const auto model = ExaGOStringOption("-opflow_model", "OPFLOW model name",
                                      "POWER_BALANCE_POLAR",
                                      {
+                                         "DCOPF",
 #ifdef EXAGO_ENABLE_HIOP
                                          "POWER_BALANCE_HIOP",
                                          "PBPOLRAJAHIOP",
diff --git a/include/private/psimpl.h b/include/private/psimpl.h
index 225a67cf..fc091a70 100644
--- a/include/private/psimpl.h
+++ b/include/private/psimpl.h
@@ -497,7 +497,8 @@ struct _p_PS {
 
   PetscBool setupcalled; /* Is setup called on PS? */
 
-  PetscLogDouble solve_real_time;
+  PetscLogDouble solve_real_time; /* Solve time */
+  PetscInt numits;                /* Number of solver iterations */
 };
 
 extern PetscErrorCode PSCheckTopology(PS);
diff --git a/include/private/scenariolist.h b/include/private/scenariolist.h
index 4e6bae44..5e07cc81 100644
--- a/include/private/scenariolist.h
+++ b/include/private/scenariolist.h
@@ -18,7 +18,8 @@ struct _p_Forecast {
   PetscInt nele;     /* Number of devices/elements involved in this forecast */
   PetscInt *buses;   /* Bus numbers */
   char **id;         /* Device ids */
-  PetscScalar *val;  /* forecast values */
+  PetscScalar *val1 = nullptr; /* forecast values */
+  PetscScalar *val2 = nullptr; /* 2nd set of forecast values */
 };
 
 typedef struct _p_Forecast Forecast;
diff --git a/src/opflow/interface/opflow.cpp b/src/opflow/interface/opflow.cpp
index 64c7ed60..9c42ef36 100644
--- a/src/opflow/interface/opflow.cpp
+++ b/src/opflow/interface/opflow.cpp
@@ -3196,6 +3196,7 @@ PetscErrorCode OPFLOWSetSummaryStats(OPFLOW opflow) {
   CHKERRQ(ierr);
 
   ps->solve_real_time = opflow->solve_real_time;
+  ps->numits = opflow->numits;
 
   PetscFunctionReturn(0);
 }
diff --git a/src/ps/ps.cpp b/src/ps/ps.cpp
index 15a7c3d7..2053afcf 100644
--- a/src/ps/ps.cpp
+++ b/src/ps/ps.cpp
@@ -1609,6 +1609,7 @@ PetscErrorCode PSApplyScenario(PS ps, Scenario scenario) {
   PetscInt i, j;
   Forecast *forecast;
   PSGEN gen;
+  PSLOAD load;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;
@@ -1621,7 +1622,7 @@ PetscErrorCode PSApplyScenario(PS ps, Scenario scenario) {
         CHKERRQ(ierr);
         if (gen) {
           gen->pg = gen->pt =
-              forecast->val[j] /
+              forecast->val1[j] /
               ps->MVAbase;    /* Set real power generation. Note that
                    Pg upper limit is also set to the forecast value. This allows
                    the    wind generator to be dispatched at its limit */
@@ -1638,6 +1639,34 @@ PetscErrorCode PSApplyScenario(PS ps, Scenario scenario) {
                  forecast->buses[j], forecast->id[j]);
         }
       }
+    } else if (forecast->type == FORECAST_LOAD_P) {
+      /* Load forecast */
+      for (j = 0; j < forecast->nele; j++) {
+        ierr = PSGetLoad(ps, forecast->buses[j], forecast->id[j], &load);
+        CHKERRQ(ierr);
+        if (load) {
+          load->pl = forecast->val1[j] / ps->MVAbase; /* Set real power load. */
+        } else {
+          printf("No load on bus %d with id %s. Cannot apply the "
+                 "requested scenario\n",
+                 forecast->buses[j], forecast->id[j]);
+        }
+      }
+    } else if (forecast->type == FORECAST_LOAD_PQ) {
+      /* Load forecast */
+      for (j = 0; j < forecast->nele; j++) {
+        ierr = PSGetLoad(ps, forecast->buses[j], forecast->id[j], &load);
+        CHKERRQ(ierr);
+        if (load) {
+          load->pl = forecast->val1[j] / ps->MVAbase; /* Set real power load. */
+          load->ql =
+              forecast->val2[j] / ps->MVAbase; /* Set reactive power load. */
+        } else {
+          printf("No load on bus %d with id %s. Cannot apply the "
+                 "requested scenario\n",
+                 forecast->buses[j], forecast->id[j]);
+        }
+      }
     }
   }
 
diff --git a/src/ps/psoutput.cpp b/src/ps/psoutput.cpp
index a7f3f7e7..d2642fba 100644
--- a/src/ps/psoutput.cpp
+++ b/src/ps/psoutput.cpp
@@ -272,6 +272,9 @@ PetscErrorCode PSSaveSolution_MATPOWER(PS ps, const char outfile[]) {
   fprintf(fd, "\n%%%% solve time\n");
   fprintf(fd, "%ssolve_time = %.5g;\n", prefix, ps->solve_real_time);
 
+  fprintf(fd, "\n%%%% number of iterations\n");
+  fprintf(fd, "%snum_iter = %d;\n", prefix, ps->numits);
+
   fclose(fd);
   PetscFunctionReturn(0);
 }
@@ -888,7 +891,8 @@ PetscErrorCode PSSaveSolution_JSON(PS ps, const char outfile[]) {
   PrintJSONArray(fd, "LOAD", 2, &ps->sys_info.total_load[0], true);
   PrintJSONArray(fd, "LOADSHED", 2, &ps->sys_info.total_loadshed[0], true);
 
-  PrintJSONDouble(fd, "SolveRealTime", ps->solve_real_time, false);
+  PrintJSONDouble(fd, "SolveRealTime", ps->solve_real_time, true);
+  PrintJSONInt(fd, "NumIter", ps->numits, false);
 
   PrintJSONObjectEnd(fd, false); // System summary object end
 
@@ -944,6 +948,7 @@ PetscErrorCode PSSaveSolution_MINIMAL(PS ps, const char outfile[]) {
   fprintf(fd, "\tTotal Load Shed P, Q: %9g, %9g\n",
           ps->sys_info.total_loadshed[0], ps->sys_info.total_loadshed[1]);
   fprintf(fd, "\tSolve Time: %5g\n", ps->solve_real_time);
+  fprintf(fd, "\tNumber of iterations: %d\n", ps->numits);
   fprintf(fd, "\tNzones: %d\n", ps->nzones);
   fprintf(fd, "\tNareas: %d\n", ps->nareas);
 
diff --git a/src/scopflow/interface/scopflow.cpp b/src/scopflow/interface/scopflow.cpp
index dcf45fce..35346192 100644
--- a/src/scopflow/interface/scopflow.cpp
+++ b/src/scopflow/interface/scopflow.cpp
@@ -489,6 +489,7 @@ PetscErrorCode SCOPFLOWSetUp(SCOPFLOW scopflow) {
   PetscInt c, i, j;
   PS ps;
   OPFLOW opflow;
+  PetscBool issolver_ipopt;
 
   char ploadprofile[PETSC_MAX_PATH_LEN];
   char qloadprofile[PETSC_MAX_PATH_LEN];
@@ -652,8 +653,8 @@ PetscErrorCode SCOPFLOWSetUp(SCOPFLOW scopflow) {
 
   ExaGOLog(EXAGO_LOG_INFO,
            "SCOPFLOW running with {:d} subproblems (base case + {:d} "
-           "contingencies)",
-           scopflow->nc, scopflow->nc - 1);
+           "contingencies) globally; {:d} subproblems locally",
+           scopflow->Nc, scopflow->Nc - 1, scopflow->nc);
 
   /* Set model */
   if (!scopflow->ismultiperiod) {
@@ -694,11 +695,9 @@ PetscErrorCode SCOPFLOWSetUp(SCOPFLOW scopflow) {
     ierr = OPFLOWIgnoreLineflowConstraints(
         scopflow->opflow0, scopflow->ignore_lineflow_constraints);
     CHKERRQ(ierr);
-    /* Base-case problem model should be POWER_BALANCE_POLAR */
-    ierr = OPFLOWSetModel(scopflow->opflow0, OPFLOWMODEL_PBPOL);
+    ierr = OPFLOWSetModel(scopflow->opflow0, scopflow->subproblem_model);
     CHKERRQ(ierr);
-    /* Base-case problem solver should be IPOPT */
-    ierr = OPFLOWSetSolver(scopflow->opflow0, OPFLOWSOLVER_IPOPT);
+    ierr = OPFLOWSetSolver(scopflow->opflow0, scopflow->subproblem_solver);
     CHKERRQ(ierr);
     ierr = OPFLOWReadMatPowerData(scopflow->opflow0, scopflow->netfile);
     CHKERRQ(ierr);
@@ -800,18 +799,20 @@ PetscErrorCode SCOPFLOWSetUp(SCOPFLOW scopflow) {
       }
 
       if (scopflow->cstart + c == 0) { /* First stage */
-        ierr = OPFLOWSetModel(scopflow->opflows[c], OPFLOWMODEL_PBPOL);
+        ierr = OPFLOWSetModel(scopflow->opflows[c], scopflow->subproblem_model);
         CHKERRQ(ierr);
-        ierr = OPFLOWSetSolver(scopflow->opflows[c], OPFLOWSOLVER_IPOPT);
+        ierr =
+            OPFLOWSetSolver(scopflow->opflows[c], scopflow->subproblem_solver);
         CHKERRQ(ierr);
         ierr = OPFLOWSetObjectiveType(scopflow->opflows[c], MIN_GEN_COST);
         CHKERRQ(ierr);
       } else { /* Second stages */
         ierr = OPFLOWHasGenSetPoint(scopflow->opflows[c], PETSC_TRUE);
         CHKERRQ(ierr); /* Activates ramping variables */
-        ierr = OPFLOWSetModel(scopflow->opflows[c], OPFLOWMODEL_PBPOL);
+        ierr = OPFLOWSetModel(scopflow->opflows[c], scopflow->subproblem_model);
         CHKERRQ(ierr);
-        ierr = OPFLOWSetSolver(scopflow->opflows[c], OPFLOWSOLVER_IPOPT);
+        ierr =
+            OPFLOWSetSolver(scopflow->opflows[c], scopflow->subproblem_solver);
         CHKERRQ(ierr);
         //	ierr = OPFLOWSetObjectiveType(scopflow->opflows[c], NO_OBJ);
         CHKERRQ(ierr);
@@ -935,40 +936,48 @@ PetscErrorCode SCOPFLOWSetUp(SCOPFLOW scopflow) {
     }
   }
 
-  /* Create vector X */
-  ierr = VecCreate(scopflow->comm->type, &scopflow->X);
-  CHKERRQ(ierr);
-  ierr = VecSetSizes(scopflow->X, scopflow->nx, PETSC_DECIDE);
-  CHKERRQ(ierr);
-  ierr = VecSetFromOptions(scopflow->X);
-  CHKERRQ(ierr);
-  ierr = VecGetSize(scopflow->X, &scopflow->Nx);
+  ierr = PetscStrcmp(scopflow->solvername.c_str(), "IPOPT", &issolver_ipopt);
   CHKERRQ(ierr);
+  if (issolver_ipopt) {
+    /* Create vector X */
+    ierr = VecCreate(scopflow->comm->type, &scopflow->X);
+    CHKERRQ(ierr);
+    ierr = VecSetSizes(scopflow->X, scopflow->nx, PETSC_DECIDE);
+    CHKERRQ(ierr);
+    ierr = VecSetFromOptions(scopflow->X);
+    CHKERRQ(ierr);
+    ierr = VecGetSize(scopflow->X, &scopflow->Nx);
+    CHKERRQ(ierr);
 
-  ierr = VecDuplicate(scopflow->X, &scopflow->Xl);
-  CHKERRQ(ierr);
-  ierr = VecDuplicate(scopflow->X, &scopflow->Xu);
-  CHKERRQ(ierr);
-  ierr = VecDuplicate(scopflow->X, &scopflow->gradobj);
-  CHKERRQ(ierr);
+    ierr = VecDuplicate(scopflow->X, &scopflow->Xl);
+    CHKERRQ(ierr);
+    ierr = VecDuplicate(scopflow->X, &scopflow->Xu);
+    CHKERRQ(ierr);
+    ierr = VecDuplicate(scopflow->X, &scopflow->gradobj);
+    CHKERRQ(ierr);
 
-  /* Vector for constraints */
-  ierr = VecCreate(scopflow->comm->type, &scopflow->G);
-  CHKERRQ(ierr);
-  ierr = VecSetSizes(scopflow->G, scopflow->ncon, PETSC_DECIDE);
-  CHKERRQ(ierr);
+    /* Vector for constraints */
+    ierr = VecCreate(scopflow->comm->type, &scopflow->G);
+    CHKERRQ(ierr);
+    ierr = VecSetSizes(scopflow->G, scopflow->ncon, PETSC_DECIDE);
+    CHKERRQ(ierr);
 
-  ierr = VecSetFromOptions(scopflow->G);
-  CHKERRQ(ierr);
+    ierr = VecSetFromOptions(scopflow->G);
+    CHKERRQ(ierr);
 
-  ierr = VecGetSize(scopflow->G, &scopflow->Ncon);
-  CHKERRQ(ierr);
+    ierr = VecGetSize(scopflow->G, &scopflow->Ncon);
+    CHKERRQ(ierr);
 
-  /* Constraint bounds vectors  */
-  ierr = VecDuplicate(scopflow->G, &scopflow->Gl);
-  CHKERRQ(ierr);
-  ierr = VecDuplicate(scopflow->G, &scopflow->Gu);
-  CHKERRQ(ierr);
+    /* Constraint bounds vectors  */
+    ierr = VecDuplicate(scopflow->G, &scopflow->Gl);
+    CHKERRQ(ierr);
+    ierr = VecDuplicate(scopflow->G, &scopflow->Gu);
+    CHKERRQ(ierr);
+
+    /* Lagrangian multipliers */
+    ierr = VecDuplicate(scopflow->G, &scopflow->Lambda);
+    CHKERRQ(ierr);
+  }
 
   /* The matrices are not used in parallel, so we don't need to create them */
   if (scopflow->comm->size == 1) {
@@ -1007,10 +1016,6 @@ PetscErrorCode SCOPFLOWSetUp(SCOPFLOW scopflow) {
     CHKERRQ(ierr);
   }
 
-  /* Lagrangian multipliers */
-  ierr = VecDuplicate(scopflow->G, &scopflow->Lambda);
-  CHKERRQ(ierr);
-
   ierr = (*scopflow->solverops.setup)(scopflow);
   CHKERRQ(ierr);
   ExaGOLog(EXAGO_LOG_INFO, "{}", "SCOPFLOW: Setup completed");
diff --git a/src/sopflow/interface/sopflow.cpp b/src/sopflow/interface/sopflow.cpp
index dfdbc745..7c00f8ad 100644
--- a/src/sopflow/interface/sopflow.cpp
+++ b/src/sopflow/interface/sopflow.cpp
@@ -204,7 +204,9 @@ PetscErrorCode SOPFLOWDestroy(SOPFLOW *sopflow) {
         }
         ierr = PetscFree((*sopflow)->scenlist.scen[s].forecastlist[i].id);
         CHKERRQ(ierr);
-        ierr = PetscFree((*sopflow)->scenlist.scen[s].forecastlist[i].val);
+        ierr = PetscFree((*sopflow)->scenlist.scen[s].forecastlist[i].val1);
+        CHKERRQ(ierr);
+        ierr = PetscFree((*sopflow)->scenlist.scen[s].forecastlist[i].val2);
         CHKERRQ(ierr);
       }
     }
diff --git a/src/sopflow/interface/sopflowscen.cpp b/src/sopflow/interface/sopflowscen.cpp
index df0d642b..52e6d635 100644
--- a/src/sopflow/interface/sopflowscen.cpp
+++ b/src/sopflow/interface/sopflowscen.cpp
@@ -4,6 +4,9 @@
 #include <private/sopflowimpl.h>
 #include <private/tcopflowimpl.h>
 
+#include <map>
+using namespace std;
+
 extern void clean2Char(char *);
 extern char **blankTokenizer(const char *str, int *numtok, int maxtokens,
                              int maxchar);
@@ -35,6 +38,152 @@ PetscErrorCode SOPFLOWSetScenarioData(SOPFLOW sopflow,
   PetscFunctionReturn(0);
 }
 
+/*
+  SOPFLOWReadScenarioData_LoadPQ_SinglePeriod - Reads the load data and
+populates the scenario list Input Parameters
++ sopflow - SOPFLOW object
+. loadprofile - load profile file
+
+  Note: This function reads the load scenario data for "single period" format
+files.
+*/
+PetscErrorCode
+SOPFLOWReadScenarioData_LoadPQ_SinglePeriod(SOPFLOW sopflow,
+                                            const char windgenprofile[]) {
+
+  PetscErrorCode ierr;
+  FILE *fp;
+  char line[MAXLINE];
+  char *out;
+  PetscInt ngenwind, nw = 0;
+  char *tok, *tok2, *tok3, *tok4, *tok5, *tok6;
+  char sep_comma[] = ",", sep2_dash[] = "_";
+  PetscReal pl, ql, weight;
+  int scen_num = 0;
+  int genid;
+  int bus[1000];
+  char genids[1000][3];
+  int i;
+  ScenarioList *scenlist = &sopflow->scenlist;
+  Scenario *scenario;
+  Forecast *forecast;
+
+  PetscFunctionBegin;
+
+  ngenwind = 1000; // This should be increased for larger cases (?)
+  fp = fopen(windgenprofile, "r");
+  if (fp == NULL) {
+    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN,
+            "Cannot open wind generation profile file %s", windgenprofile);
+  }
+
+  /* First line -- has the bus numbers */
+  out = fgets(line, MAXLINE, fp);
+  /* Parse wind generator numbers */
+
+  // printf("Done Parsing Header. nw = %d\n", nw);
+
+  // Body Example
+  // 1,0.1,8_Load_1_185.0_10,6_Load_1_105.5_20
+  while ((out = fgets(line, MAXLINE, fp)) != NULL) {
+    if (strcmp(line, "\r\n") == 0 || strcmp(line, "\n") == 0) {
+      continue; /* Skip blank lines */
+    }
+
+    tok = strtok(line, sep_comma);
+    sscanf(tok, "%d", &scen_num); /* Scenario number */
+    scen_num -= 1; /* Scenario numbers start with 1 in the file, convert to
+                      zero-based start */
+
+    if (scen_num < scenlist->Nscen || scen_num == sopflow->Ns) {
+      fclose(fp);
+      PetscFunctionReturn(0);
+    }
+
+    struct ForecastData {
+      int busnum;
+      char type[16];
+      PetscReal value1;
+      PetscReal value2;
+      char id[3];
+    };
+
+    vector<ForecastData> gendata;
+
+    tok = strtok(NULL, sep_comma);
+    sscanf(tok, "%lf", &weight);
+
+    tok = strtok(NULL, sep_comma);
+    while (tok != NULL) {
+
+      ForecastData gd;
+
+      /* Parse generator info */
+      tok2 = strsep(&tok, sep2_dash); // Bus Number
+      sscanf(tok2, "%d", &gd.busnum);
+
+      tok3 = strsep(&tok, sep2_dash); // Type (Load or Gen)
+      sscanf(tok3, "%s", gd.type);
+
+      tok4 = strsep(&tok, sep2_dash); // CKT
+      sscanf(tok4, "%d", &genid);
+      snprintf(gd.id, 3, "%-2d", genid);
+
+      if (strcmp(gd.type, "Load") == 0) {
+        /* Load forecast */
+        tok5 = strsep(&tok, sep2_dash); // Value1
+        sscanf(tok5, "%lf", &gd.value1);
+        tok6 = strsep(&tok, sep2_dash); // Value2
+        sscanf(tok6, "%lf", &gd.value2);
+      } else if (strcmp(gd.type, "Gen") == 0) {
+        /* Gen forecast */
+        tok5 = strsep(&tok, sep2_dash); // Value1
+        sscanf(tok5, "%lf", &gd.value1);
+      }
+
+      gendata.push_back(gd);
+
+      tok = strtok(NULL, sep_comma);
+    }
+
+    scenario = &scenlist->scen[scen_num];
+
+    printf("scenario %d has %d forecasts\n", scen_num, scenario->nforecast);
+
+    forecast = &scenario->forecastlist[scenario->nforecast];
+    forecast->num = scen_num;
+    forecast->type = FORECAST_LOAD_PQ;
+    forecast->nele = gendata.size();
+    ierr = PetscCalloc1(forecast->nele, &forecast->buses);
+    CHKERRQ(ierr);
+    ierr = PetscCalloc1(forecast->nele, &forecast->id);
+    CHKERRQ(ierr);
+    for (i = 0; i < forecast->nele; i++) {
+      ierr = PetscCalloc1(3, &forecast->id[i]);
+      CHKERRQ(ierr);
+    }
+    ierr = PetscCalloc1(forecast->nele, &forecast->val1);
+    CHKERRQ(ierr);
+    ierr = PetscCalloc1(forecast->nele, &forecast->val2);
+    CHKERRQ(ierr);
+
+    for (i = 0; i < gendata.size(); i++) {
+      forecast->buses[i] = gendata[i].busnum;
+      ierr = PetscStrcpy(forecast->id[i], gendata[i].id);
+      CHKERRQ(ierr);
+      forecast->val1[i] = gendata[i].value1;
+      forecast->val2[i] = gendata[i].value2;
+    }
+
+    // /* Read scenario weight */
+    scenario->prob = weight;
+
+    scenario->nforecast++;
+    scenlist->Nscen++;
+  }
+  PetscFunctionReturn(0);
+}
+
 /*
   SOPFLOWReadScenarioData_Wind_SinglePeriod - Reads the wind data and populates
 the scenario list Input Parameters
@@ -126,7 +275,7 @@ SOPFLOWReadScenarioData_Wind_SinglePeriod(SOPFLOW sopflow,
     for (i = 0; i < forecast->nele; i++) {
       ierr = PetscCalloc1(3, &forecast->id[i]);
     }
-    ierr = PetscCalloc1(forecast->nele, &forecast->val);
+    ierr = PetscCalloc1(forecast->nele, &forecast->val1);
     CHKERRQ(ierr);
 
     tok = strtok(NULL, sep);
@@ -135,7 +284,7 @@ SOPFLOWReadScenarioData_Wind_SinglePeriod(SOPFLOW sopflow,
       ierr = PetscStrcpy(forecast->id[i], windgenid[i]);
       CHKERRQ(ierr);
       sscanf(tok, "%lf", &pg);
-      forecast->val[i] = pg;
+      forecast->val1[i] = pg;
       tok = strtok(NULL, sep);
     }
 
@@ -238,7 +387,7 @@ SOPFLOWReadScenarioData_Wind_MultiPeriod(SOPFLOW sopflow,
     for (i = 0; i < forecast->nele; i++) {
       ierr = PetscCalloc1(3, &forecast->id[i]);
     }
-    ierr = PetscCalloc1(forecast->nele, &forecast->val);
+    ierr = PetscCalloc1(forecast->nele, &forecast->val1);
     CHKERRQ(ierr);
 
     for (i = 0; i < nw; i++) {
@@ -251,7 +400,7 @@ SOPFLOWReadScenarioData_Wind_MultiPeriod(SOPFLOW sopflow,
     nw = 0;
     while (tok != NULL) {
       sscanf(tok, "%lf", &pg);
-      forecast->val[nw] = pg;
+      forecast->val1[nw] = pg;
       nw++;
       tok = strtok(NULL, sep);
     }
@@ -284,6 +433,9 @@ PetscErrorCode SOPFLOWReadScenarioData(SOPFLOW sopflow,
       ierr = SOPFLOWReadScenarioData_Wind_MultiPeriod(sopflow, scenfile);
       CHKERRQ(ierr);
     }
+  } else if (sopflow->scenunctype == LOAD) {
+    ierr = SOPFLOWReadScenarioData_LoadPQ_SinglePeriod(sopflow, scenfile);
+    CHKERRQ(ierr);
   }
 
   PetscFunctionReturn(0);
diff --git a/tests/functionality/sopflow/selfcheck.cpp b/tests/functionality/sopflow/selfcheck.cpp
index 924779b8..8eb1d22a 100644
--- a/tests/functionality/sopflow/selfcheck.cpp
+++ b/tests/functionality/sopflow/selfcheck.cpp
@@ -9,6 +9,7 @@ struct SopflowFunctionalityTestParameters {
   std::string solver = "";
   std::string network = "";
   std::string scenfile = "";
+  std::string loadfile = "";
   std::string contingencies = "";
   std::string pload = "";
   std::string qload = "";
@@ -54,6 +55,7 @@ struct SopflowFunctionalityTestParameters {
     set_if_found(solver, values, "solver");
     set_if_found(network, values, "network");
     set_if_found(scenfile, values, "scenfile");
+    set_if_found(loadfile, values, "loadfile");
     set_if_found(num_scenarios, values, "num_scenarios");
     set_if_found(tolerance, values, "tolerance");
     set_if_found(warning_tolerance, values, "warning_tolerance");
@@ -138,8 +140,7 @@ struct SopflowFunctionalityTests
       }
     };
 
-    for (const auto &opt :
-         {"solver", "network", "scenfile", "num_scenarios", "tolerance"})
+    for (const auto &opt : {"solver", "network", "num_scenarios", "tolerance"})
       ensure_option_available(opt);
 
     bool is_multicontingency = false;
@@ -171,7 +172,10 @@ struct SopflowFunctionalityTests
     testcase["description"] = params.description;
     testcase["solver"] = params.solver;
     testcase["network"] = params.network;
-    testcase["scenfile"] = params.scenfile;
+    if (params.loadfile != "")
+      testcase["loadfile"] = params.loadfile;
+    else if (params.scenfile != "")
+      testcase["scenfile"] = params.scenfile;
     testcase["num_scenarios"] = params.num_scenarios;
     testcase["initialization_type"] = params.initialization_type;
 
@@ -231,10 +235,21 @@ struct SopflowFunctionalityTests
     ExaGOCheckError(ierr);
 
     // Prepend installation directory to scenario data
-    resolve_datafiles_path(params.scenfile);
-    ierr = SOPFLOWSetScenarioData(sopflow, SOPFLOW_NATIVE_SINGLEPERIOD, WIND,
-                                  params.scenfile.c_str());
-    ExaGOCheckError(ierr);
+    if (params.loadfile != "") {
+      std::cout << "Network: " << params.network
+                << " Using Load File: " << params.loadfile << std::endl;
+      resolve_datafiles_path(params.loadfile);
+      ierr = SOPFLOWSetScenarioData(sopflow, SOPFLOW_NATIVE_SINGLEPERIOD, LOAD,
+                                    params.loadfile.c_str());
+      ExaGOCheckError(ierr);
+    } else if (params.scenfile != "") {
+      std::cout << "Network: " << params.network
+                << " Using Scenario File: " << params.scenfile << std::endl;
+      resolve_datafiles_path(params.scenfile);
+      ierr = SOPFLOWSetScenarioData(sopflow, SOPFLOW_NATIVE_SINGLEPERIOD, WIND,
+                                    params.scenfile.c_str());
+      ExaGOCheckError(ierr);
+    }
 
     ierr = SOPFLOWSetInitializationType(sopflow, params.initialization_type);
     ExaGOCheckError(ierr);
diff --git a/tests/functionality/sopflow/sopflow_multiscenario.toml b/tests/functionality/sopflow/sopflow_multiscenario.toml
index 6beb2e6b..39c4ddc5 100644
--- a/tests/functionality/sopflow/sopflow_multiscenario.toml
+++ b/tests/functionality/sopflow/sopflow_multiscenario.toml
@@ -53,3 +53,13 @@ solver = 'EMPAR'
 num_scenarios = 3
 num_iters = 0
 obj_value = 14316.630954933302
+
+[[testcase]]
+description = 'LoadTest'
+tolerance = 1e-3
+network = 'datafiles/case9/case9mod.m'
+loadfile = 'datafiles/case9/10_scenarios_9bus_load.txt'
+solver = 'EMPAR'
+num_scenarios = 3
+num_iters = 0
+obj_value = 722.106958693