From 06597a0420bc631be45b6dd9e4c6e67e838e9992 Mon Sep 17 00:00:00 2001
From: nataraj2 <mn556@cornell.edu>
Date: Mon, 17 Aug 2026 12:38:10 -0700
Subject: [PATCH 01/18] Minor changes and adding readme for preprocessir

---
 .../DataAssimilation/PreProcessing/README.md  | 22 +++++++++++++++++++
 .../DataAssimilation/SquallLine2D/main.cpp    |  7 +++++-
 2 files changed, 28 insertions(+), 1 deletion(-)
 create mode 100644 .Exec_dev/DataAssimilation/PreProcessing/README.md

diff --git a/.Exec_dev/DataAssimilation/PreProcessing/README.md b/.Exec_dev/DataAssimilation/PreProcessing/README.md
new file mode 100644
index 0000000000..cce4a24a00
--- /dev/null
+++ b/.Exec_dev/DataAssimilation/PreProcessing/README.md
@@ -0,0 +1,22 @@
+# Pre-processing for Data Assimilation (2D Squall Line)
+
+This directory contains the pre-processing code that generates the background initial conditions for ensemble-based data assimilation.
+
+## Overview
+
+The code takes a fine-mesh "ground truth" plotfile and a coarse-mesh plotfile, then interpolates the fine-mesh data onto the coarse mesh. 
+
+**Concept:** Satellite observations capture coarse representations of flow features—the features are present, but at a lower resolution. To mimic this, this tool generates a coarsened version of the ground truth and exports it to a custom binary data file used as input for ensemble simulations.
+
+## Build and Run
+
+1. **Compile:**
+   ```bash
+   make -j8
+   ```
+
+2. **Execute:**
+   ```bash
+   ./out -- <fine-mesh-plt-file> <coarse-mesh-plt-file>
+   ```
+Note: The coarse-mesh plotfile **must contain only 1 box**. Ensure it is run on a single MPI rank with a sufficiently large `max_grid_size`.
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
index f7f7bd48f5..932a1ed68a 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
@@ -89,6 +89,11 @@ int main (int argc, char* argv[])
     // ------------------------------------------------------------
     // Ensemble control
     // ------------------------------------------------------------
+
+    ParmParse pp_erf("erf");
+    bool is_init_for_ensemble = false;
+    pp_erf.query("is_init_for_ensemble", is_init_for_ensemble);
+
     ParmParse pp_ens("ensemble");
 
     int n_ens = 1;
@@ -120,7 +125,7 @@ int main (int argc, char* argv[])
     // --------------------------------------------------------
     ParallelDescriptor::Barrier();
 
-    if (ParallelDescriptor::IOProcessor())
+    if (is_init_for_ensemble && ParallelDescriptor::IOProcessor())
     {
         // Create zero-padded member directory
         std::stringstream ss;

From fbb6c9fb6e0879c5e4fe708f05641826cfe959a9 Mon Sep 17 00:00:00 2001
From: nataraj2 <mn556@cornell.edu>
Date: Mon, 17 Aug 2026 12:41:48 -0700
Subject: [PATCH 02/18] Updating readme

---
 .Exec_dev/DataAssimilation/PreProcessing/README.md | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/.Exec_dev/DataAssimilation/PreProcessing/README.md b/.Exec_dev/DataAssimilation/PreProcessing/README.md
index cce4a24a00..46fddb5709 100644
--- a/.Exec_dev/DataAssimilation/PreProcessing/README.md
+++ b/.Exec_dev/DataAssimilation/PreProcessing/README.md
@@ -1,12 +1,12 @@
-# Pre-processing for Data Assimilation (2D Squall Line)
+# Pre-processing for Data Assimilation
 
 This directory contains the pre-processing code that generates the background initial conditions for ensemble-based data assimilation.
 
 ## Overview
 
-The code takes a fine-mesh "ground truth" plotfile and a coarse-mesh plotfile, then interpolates the fine-mesh data onto the coarse mesh. 
+The code takes a fine-mesh "ground truth" plotfile and a coarse-mesh plotfile, then interpolates the fine-mesh data onto the coarse mesh and writes out a custom-written binary file with the coarse data, that can be read-in by the data assimilation code. 
 
-**Concept:** Satellite observations capture coarse representations of flow features—the features are present, but at a lower resolution. To mimic this, this tool generates a coarsened version of the ground truth and exports it to a custom binary data file used as input for ensemble simulations.
+**Concept:** Satellite observations capture coarse representations of flow features ie. the features are present, but at a lower resolution. To mimic this, this tool generates a coarsened version of the ground truth and exports it to a custom binary data file used as input for ensemble simulations.
 
 ## Build and Run
 
@@ -14,8 +14,9 @@ The code takes a fine-mesh "ground truth" plotfile and a coarse-mesh plotfile, t
    ```bash
    make -j8
    ```
+2. `vars.txt` should contain the list of variables to be included in the coarsened data.
 
-2. **Execute:**
+3. **Execute:**
    ```bash
    ./out -- <fine-mesh-plt-file> <coarse-mesh-plt-file>
    ```

From aaad5729de8e02b96edc791953dc4b927d628fb1 Mon Sep 17 00:00:00 2001
From: nataraj2 <mn556@cornell.edu>
Date: Tue, 18 Aug 2026 15:34:37 -0700
Subject: [PATCH 03/18] Minor error correction

---
 Source/Initialization/ERF_InitForEnsemble.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Source/Initialization/ERF_InitForEnsemble.cpp b/Source/Initialization/ERF_InitForEnsemble.cpp
index c99499375c..bb9d3a72ab 100644
--- a/Source/Initialization/ERF_InitForEnsemble.cpp
+++ b/Source/Initialization/ERF_InitForEnsemble.cpp
@@ -639,11 +639,11 @@ MakeFinalMultiFabs (const MultiFab& mf_cc_fine,
     {
         const Box& bx = mfi.tilebox();
         auto const& wface = zvel_pert.array(mfi);
-        //auto const& cc    = mf_cc_fine.const_array(mfi);
+        auto const& cc    = mf_cc_fine.const_array(mfi);
 
         amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
         {
-            wface(i,j,k) = zero; //myhalf * (cc(i,j,k-1,4) + cc(i,j,k,4));
+            wface(i,j,k) = myhalf * (cc(i,j,k-1,4) + cc(i,j,k,4));
         });
     }
 }

From 1b10a791bf1072dbc9fa5256625932f58370bc5d Mon Sep 17 00:00:00 2001
From: nataraj2 <mn556@cornell.edu>
Date: Tue, 18 Aug 2026 15:58:52 -0700
Subject: [PATCH 04/18] Adding functions for ensemble updates

---
 .../SquallLine2D/ERF_DA_EnKFSRF.cpp           | 196 +++++++++++++++---
 1 file changed, 170 insertions(+), 26 deletions(-)

diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
index d4e46fce16..1c3d221f5a 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
@@ -14,6 +14,131 @@
 using namespace amrex;
 namespace fs = std::filesystem;
 
+void ApplyNeumannBCsToEnsembles(const Geometry& geom,
+                               MultiFab& mf_cc)
+{
+
+     // -------------------------------------------------
+    // 2. Fill interior + periodic ghost cells
+    // -------------------------------------------------
+    mf_cc.FillBoundary(geom.periodicity());
+    // -------------------------------------------------
+    // 3. Apply FOExtrap (Neumann) at domain boundaries
+    // -------------------------------------------------
+    const Box& domain = geom.Domain();
+
+    for (MFIter mfi(mf_cc, TilingIfNotGPU()); mfi.isValid(); ++mfi)
+    {
+        const Box& gbx = mfi.growntilebox();   // includes ghost cells
+        const Box& vbx = mfi.validbox();
+
+        auto const& arr = mf_cc.array(mfi);
+        int ncomp = mf_cc.nComp();
+
+        ParallelFor(gbx, ncomp,
+        [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
+        {
+            if (vbx.contains(i,j,k)) return;
+
+            int ii = i;
+            int jj = j;
+            int kk = k;
+
+            // Clamp to domain interior (FOExtrap)
+            ii = amrex::max(domain.smallEnd(0),
+                 amrex::min(i, domain.bigEnd(0)));
+
+            jj = amrex::max(domain.smallEnd(1),
+                 amrex::min(j, domain.bigEnd(1)));
+
+            kk = amrex::max(domain.smallEnd(2),
+                 amrex::min(k, domain.bigEnd(2)));
+
+            arr(i,j,k,n) = arr(ii,jj,kk,n);
+        });
+    }
+}
+
+/**
+ * Split cell-centered interpolated background data into ERF state and face velocities.
+ *
+ * @param mf_cc_fine Fine-grid cell-centered source data
+ * @param cons_pert Conserved-state perturbation MultiFab to fill
+ * @param xvel_pert x-face velocity perturbation MultiFab to fill
+ * @param yvel_pert y-face velocity perturbation MultiFab to fill
+ * @param zvel_pert z-face velocity perturbation MultiFab to fill
+ */
+void
+WriteUpdatedEnsembleToERFClassData (const MultiFab& mf_cc_fine,
+                                    MultiFab& cons_pert,
+                                    MultiFab& xvel_pert,
+                                    MultiFab& yvel_pert,
+                                    MultiFab& zvel_pert,
+                                    const int n_qstate_moist)
+{
+
+    for (MFIter mfi(cons_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi)
+    {
+        const Box& bx = mfi.tilebox();
+
+        auto const& mf_cc_fine_arr  = mf_cc_fine.const_array(mfi);
+        auto const& cons_pert_arr = cons_pert.array(mfi);
+
+        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+        {
+            Real tmp_rho   = mf_cc_fine_arr(i,j,k,0);
+            Real tmp_theta = mf_cc_fine_arr(i,j,k,1);
+            Real tmp_qv    = mf_cc_fine_arr(i,j,k,5);
+            Real tmp_qc    = mf_cc_fine_arr(i,j,k,6);
+            Real tmp_qrain = mf_cc_fine_arr(i,j,k,7);
+            cons_pert_arr(i,j,k,Rho_comp)      = tmp_rho;
+            cons_pert_arr(i,j,k,RhoTheta_comp) = tmp_rho*tmp_theta;
+            if (n_qstate_moist > 0) cons_pert_arr(i,j,k,RhoQ1_comp)    = tmp_rho*tmp_qv;
+            if (n_qstate_moist > 1) cons_pert_arr(i,j,k,RhoQ2_comp)    = tmp_rho*tmp_qc;
+            if (n_qstate_moist > 2) cons_pert_arr(i,j,k,RhoQ3_comp)    = tmp_rho*tmp_qrain;
+        });
+    }
+
+    // --- X-faces (component 2) ---
+    for (MFIter mfi(xvel_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi)
+    {
+        const Box& bx = mfi.tilebox();
+        auto const& uface = xvel_pert.array(mfi);
+        auto const& cc    = mf_cc_fine.const_array(mfi);
+
+        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+        {
+            uface(i,j,k) = myhalf * (cc(i-1,j,k,2) + cc(i,j,k,2));
+        });
+    }
+
+    // --- Y-faces (component 3) ---
+    for (MFIter mfi(yvel_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi)
+    {
+        const Box& bx = mfi.tilebox();
+        auto const& vface = yvel_pert.array(mfi);
+        auto const& cc    = mf_cc_fine.const_array(mfi);
+
+        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+        {
+            vface(i,j,k) = myhalf * (cc(i,j-1,k,3) + cc(i,j,k,3));
+        });
+    }
+
+    // --- Z-faces (component 4) ---
+    for (MFIter mfi(zvel_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi)
+    {
+        const Box& bx = mfi.tilebox();
+        auto const& wface = zvel_pert.array(mfi);
+        auto const& cc    = mf_cc_fine.const_array(mfi);
+
+        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+        {
+            wface(i,j,k) = myhalf * (cc(i,j,k-1,4) + cc(i,j,k,4));
+        });
+    }
+}
+
 void
 ERF::ComputeAndWriteEnsemblePerturbations()
 {
@@ -24,7 +149,7 @@ ERF::ComputeAndWriteEnsemblePerturbations()
     // ie.find the ensemble mean at iteration 100, loop over the plt00100 file in each of the
     // member*/plotfiles/plt00100
     int Nens = solverChoice.n_ensemble;
-    Vector<std::string> varnames = {"density","theta", "x_velocity","y_velocity","z_velocity"};
+    Vector<std::string> varnames = {"density","theta", "rhoQ1", "rhoQ2", "rhoQ3", "x_velocity","y_velocity","z_velocity"};
     const std::string member_prefix = "member_";
     for (const auto& pf_name : pltfiles)
     {
@@ -75,44 +200,41 @@ ERF::PerformDataAssimilation(int da_iter)
     // ie.find the ensemble mean at iteration 100, loop over the plt00100 file in each of the
     // member*-plotfiles-plt00100
     int Nens = solverChoice.n_ensemble;
-    Vector<std::string> varnames = {"density","theta", "x_velocity","y_velocity","z_velocity"};
+    Vector<std::string> varnames = {"density","theta", "x_velocity","y_velocity","z_velocity", "qv", "qc", "qrain"};
 
     // Compute the ensemble mean
     MultiFab xf_bar = compute_ensemble_mean(Nens, last_pf_name, varnames);
 
-   // Construct perturbation plotfile name
-   std::string pltname = "plt_ens_mean";
-   WriteSingleLevelPlotfile(pltname,
-                            xf_bar,
-                            varnames,
-                            geom[0],
-                            0.0,   // time
-                            0);    // level
+    // Construct perturbation plotfile name
+    std::string pltname = "plt_ens_mean";
+    WriteSingleLevelPlotfile(pltname,
+                             xf_bar,
+                             varnames,
+                             geom[0],
+                             0.0,   // time
+                             0);    // level
 
     // Compute the mean of forecast observations yf_bar = Hx_f
     MultiFab mean_H_xf;
     compute_mean_H_xf(mean_H_xf, Nens, last_pf_name, varnames);
 
-     Vector<std::string> varnames1 = {"x_velocity", "y_velocity"};
-   // Construct perturbation plotfile name
-   std::string pltname1 = "plt_mean_H_xf";
-   WriteSingleLevelPlotfile(pltname1,
-                            mean_H_xf,
-                            varnames1,
-                            geom[0],
-                            0.0,   // time
-                            0);    // level
-
-
-
+    Vector<std::string> varnames1 = {"x_velocity", "y_velocity"};
+    // Construct perturbation plotfile name
+    std::string pltname1 = "plt_mean_H_xf";
+    WriteSingleLevelPlotfile(pltname1,
+                             mean_H_xf,
+                             varnames1,
+                             geom[0],
+                             0.0,   // time
+                             0);    // level
 
     // Read in the observation file
     MultiFab y_obs;
     read_in_observations(da_iter, varnames, y_obs);
 
     std::string pltname2 = "plt_y_obs";
-     Vector<std::string> varnames2 = {"x_velocity", "y_velocity"};
-   WriteSingleLevelPlotfile(pltname2,
+    Vector<std::string> varnames2 = {"x_velocity", "y_velocity"};
+    WriteSingleLevelPlotfile(pltname2,
                             y_obs,
                             varnames2,
                             geom[0],
@@ -124,8 +246,8 @@ ERF::PerformDataAssimilation(int da_iter)
     compute_d_vec(y_obs, mean_H_xf, d_vec);
 
     std::string pltname_diff = "plt_rhs_diff";
-     Vector<std::string> varnames_diff = {"x_velocity", "y_velocity"};
-   WriteSingleLevelPlotfile(pltname_diff,
+    Vector<std::string> varnames_diff = {"x_velocity", "y_velocity"};
+    WriteSingleLevelPlotfile(pltname_diff,
                             d_vec,
                             varnames_diff,
                             geom[0],
@@ -164,6 +286,7 @@ ERF::PerformDataAssimilation(int da_iter)
                             0.0,   // time
                             0);    // level
 
+    // Update the ensemble mean
     MultiFab xf_bar_updated;
     add_multifabs(xf_bar, Xf_prime_alpha, xf_bar_updated);
 
@@ -179,12 +302,33 @@ ERF::PerformDataAssimilation(int da_iter)
     Matrix T_mat(Nens);
     compute_T_matrix(S_mat, T_mat);
 
+    // Update all the ensembles and write checkpoint file
     for(int n=0; n< Nens; n++) {
         Print() << "Updating for ensemble " << n << std::endl;
         MultiFab mf_ens_pert;
         update_ensemble(Nens, last_pf_name, varnames, xf_bar, T_mat, n, mf_ens_pert);
+
+        // Update the ensemble
         MultiFab mf_ens_updated;
         add_multifabs(mf_ens_pert, xf_bar_updated, mf_ens_updated);
+
+        // Apply Neumann boundary condition to the updated ensembles
+        ApplyNeumannBCsToEnsembles(geom[0], mf_ens_updated);
+        bool use_moisture = (solverChoice.moisture_type != MoistureType::None);
+        int n_qstate_moist = 0;
+        if (use_moisture) {
+            n_qstate_moist = micro->Get_Qstate_Moist_Size();
+        }
+
+        auto& lev_new = vars_new[0];
+
+        // Copy the cell centered ensemble multifab to the ERF class data structures
+        WriteUpdatedEnsembleToERFClassData(mf_ens_updated,
+                                           lev_new[Vars::cons],
+                                           lev_new[Vars::xvel],
+                                           lev_new[Vars::yvel],
+                                           lev_new[Vars::zvel],
+                                           n_qstate_moist);
         check_file = "chk";
         InitData();
         check_file = MakeEnsembleCheckpointName(da_iter, n);

From 4775b5867c937b0d0a2a87ef5f2c8ca006b2183e Mon Sep 17 00:00:00 2001
From: nataraj2 <mn556@cornell.edu>
Date: Tue, 18 Aug 2026 17:06:37 -0700
Subject: [PATCH 05/18] Updating data assimilation for 2d squall line

---
 .../SquallLine2D/ERF_DA_EnKFSRF.cpp           | 48 +++++++++++++++----
 .../ERF_DA_EnKFSRF_SquallLine2D.cpp           | 16 ++++---
 .../SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp     |  6 +--
 .../DataAssimilation/SquallLine2D/main.cpp    |  4 +-
 4 files changed, 53 insertions(+), 21 deletions(-)

diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
index 1c3d221f5a..1cba0ab13e 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
@@ -196,6 +196,8 @@ ERF::PerformDataAssimilation(int da_iter)
         amrex::Abort("No plotfiles found.");
     }
 
+    Print() << "Reaching here. Beginning data assimilation" << std::endl;
+
     // Step 2: loop over all plotfiles (timestamps at which the plotfiles are written)
     // ie.find the ensemble mean at iteration 100, loop over the plt00100 file in each of the
     // member*-plotfiles-plt00100
@@ -205,6 +207,8 @@ ERF::PerformDataAssimilation(int da_iter)
     // Compute the ensemble mean
     MultiFab xf_bar = compute_ensemble_mean(Nens, last_pf_name, varnames);
 
+    Print() << "Computing ensemble mean complete" << std::endl;
+
     // Construct perturbation plotfile name
     std::string pltname = "plt_ens_mean";
     WriteSingleLevelPlotfile(pltname,
@@ -218,7 +222,9 @@ ERF::PerformDataAssimilation(int da_iter)
     MultiFab mean_H_xf;
     compute_mean_H_xf(mean_H_xf, Nens, last_pf_name, varnames);
 
-    Vector<std::string> varnames1 = {"x_velocity", "y_velocity"};
+    Print() << "compute_mean_H_xf complete" << std::endl;
+
+    Vector<std::string> varnames1 = {"density","theta", "x_velocity","y_velocity","z_velocity", "qv", "qc", "qrain"};
     // Construct perturbation plotfile name
     std::string pltname1 = "plt_mean_H_xf";
     WriteSingleLevelPlotfile(pltname1,
@@ -232,76 +238,96 @@ ERF::PerformDataAssimilation(int da_iter)
     MultiFab y_obs;
     read_in_observations(da_iter, varnames, y_obs);
 
-    std::string pltname2 = "plt_y_obs";
+    Print() << "Observation reading complete" << std::endl;
+
+    /*std::string pltname2 = "plt_y_obs";
     Vector<std::string> varnames2 = {"x_velocity", "y_velocity"};
     WriteSingleLevelPlotfile(pltname2,
                             y_obs,
                             varnames2,
                             geom[0],
                             0.0,   // time
-                            0);    // level
+                            0);    // level*/
 
     // Compute y_obs - yf_bar
     MultiFab d_vec;
     compute_d_vec(y_obs, mean_H_xf, d_vec);
 
-    std::string pltname_diff = "plt_rhs_diff";
+    Print() << "Computing dvec complete" << std::endl;
+
+    /*std::string pltname_diff = "plt_rhs_diff";
     Vector<std::string> varnames_diff = {"x_velocity", "y_velocity"};
     WriteSingleLevelPlotfile(pltname_diff,
                             d_vec,
                             varnames_diff,
                             geom[0],
                             0.0,   // time
-                            0);    // level
+                            0);    // level*/
 
 
     // Assign values for the observation error covarinace matrix
     Vector<Real> R_diag;
     compute_R_diag_vals(R_diag);
 
+    Print() << "compute_R_diag_vals complete" << std::endl;
+
     // Compute d'=R_inv*d
     MultiFab d_prime_vec;
     compute_d_prime_vec(d_prime_vec, d_vec, R_diag);
 
+    Print() << "compute_d_prime_vec complete" << std::endl;
+
     // Compute r = Y'^Td'
     Vector<Real> r_vec;
     compute_r_vec(Nens, last_pf_name, varnames, mean_H_xf, d_prime_vec, r_vec);
 
+    Print() << "compute_r_vec complete" << std::endl;
+
     // Compute the S matrix
     Matrix S_mat(Nens);
     compute_S_matrix(S_mat, Nens, mean_H_xf, R_diag, last_pf_name, varnames);
 
+    Print() << "compute_S_matrix complete" << std::endl;
+
     Vector<Real> alpha_vec;
     compute_alpha_vec(Nens, S_mat, r_vec, alpha_vec);
 
+    Print() << "compute_alpha_vec complete" << std::endl;
+
     MultiFab Xf_prime_alpha;
     compute_Xf_prime_times_vector(Nens, last_pf_name, varnames, xf_bar, alpha_vec, Xf_prime_alpha);
 
-    std::string pltname3 = "plt_Xf_prime_alpha";
+    Print() << "compute_Xf_prime_times_vector complete" << std::endl;
+
+    /*std::string pltname3 = "plt_Xf_prime_alpha";
     Vector<std::string> varnames3 = {"density", "theta", "x_velocity", "y_velocity", "z_velocity"};
     WriteSingleLevelPlotfile(pltname3,
                             Xf_prime_alpha,
                             varnames3,
                             geom[0],
                             0.0,   // time
-                            0);    // level
+                            0);    // level*/
 
     // Update the ensemble mean
     MultiFab xf_bar_updated;
     add_multifabs(xf_bar, Xf_prime_alpha, xf_bar_updated);
 
-    std::string pltname4 = "plt_xf_bar_updated";
+    Print() << "add_multifabs(xf_bar, Xf_prime_alpha, xf_bar_updated) complete" << std::endl;
+
+    /*std::string pltname4 = "plt_xf_bar_updated";
     Vector<std::string> varnames4 = {"density", "theta", "x_velocity", "y_velocity", "z_velocity"};
     WriteSingleLevelPlotfile(pltname4,
                             xf_bar_updated,
                             varnames4,
                             geom[0],
                             0.0,   // time
-                            0);    // level
+                            0);    // level*/
 
     Matrix T_mat(Nens);
     compute_T_matrix(S_mat, T_mat);
 
+    Print() << "compute_T_matrix complete" << std::endl;
+
     // Update all the ensembles and write checkpoint file
     for(int n=0; n< Nens; n++) {
         Print() << "Updating for ensemble " << n << std::endl;
@@ -320,6 +346,9 @@ ERF::PerformDataAssimilation(int da_iter)
             n_qstate_moist = micro->Get_Qstate_Moist_Size();
         }
 
+        m_plot3d_int_1 = -1;
+        m_check_int = -1;
+        InitData();
         auto& lev_new = vars_new[0];
 
         // Copy the cell centered ensemble multifab to the ERF class data structures
@@ -330,7 +359,6 @@ ERF::PerformDataAssimilation(int da_iter)
                                            lev_new[Vars::zvel],
                                            n_qstate_moist);
         check_file = "chk";
-        InitData();
         check_file = MakeEnsembleCheckpointName(da_iter, n);
         WriteCheckpointFile();
     }
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_SquallLine2D.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_SquallLine2D.cpp
index 1428973f9d..d621b013b5 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_SquallLine2D.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_SquallLine2D.cpp
@@ -22,7 +22,7 @@ Apply_H(const amrex::MultiFab& x_mf,
     // Define y_mf with 2 velocity components
     y_mf.define(x_mf.boxArray(),
                 x_mf.DistributionMap(),
-                2,
+                8,
                 x_mf.nGrowVect());
 
     for (amrex::MFIter mfi(y_mf); mfi.isValid(); ++mfi)
@@ -35,10 +35,14 @@ Apply_H(const amrex::MultiFab& x_mf,
         amrex::ParallelFor(bx,
         [=] AMREX_GPU_DEVICE (int i, int j, int k)
         {
-            amrex::Real rho = x(i,j,k,0);
-
-            y(i,j,k,0) = x(i,j,k,2); // u velocity
-            y(i,j,k,1) = x(i,j,k,3); // v velocity
+            y(i,j,k,0) = x(i,j,k,0); // density
+            y(i,j,k,1) = x(i,j,k,1); // theta
+            y(i,j,k,2) = x(i,j,k,2); // x velocity
+            y(i,j,k,3) = x(i,j,k,3); // y velocity
+            y(i,j,k,4) = x(i,j,k,4); // z velocity
+            y(i,j,k,5) = x(i,j,k,5); // qv
+            y(i,j,k,6) = x(i,j,k,6); // qc
+            y(i,j,k,7) = x(i,j,k,7); // qrain
         });
     }
 }
@@ -90,7 +94,7 @@ read_in_observations(const int& da_iter,
 void compute_R_diag_vals(Vector<Real>& R_diag)
 {
     // Set the size to 2
-    R_diag.resize(2);
+    R_diag.resize(8,0.01);
 
     R_diag[0] = 0.01;
     R_diag[1] = 0.01;
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
index c2f031a93e..6b86485207 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
@@ -190,6 +190,7 @@ compute_ensemble_mean(int Nens,
 
     for (int n = 0; n < Nens; ++n)
     {
+        std::cout << "Reading ensembles for compute_ensemble_mean = " << n << std::endl;
         MultiFab mf_tmp = read_member_multifab(n, pf_name, varnames);
 
         if (!initialized) {
@@ -281,7 +282,7 @@ compute_mean_H_xf(MultiFab& mean_H_xf,
         if (n==0) {
             mean_H_xf.define(xf_i.boxArray(),
                              xf_i.DistributionMap(),
-                             2,
+                             8,
                              xf_i.nGrow());
             mean_H_xf.setVal(0.0);
         }
@@ -533,8 +534,7 @@ compute_Xf_prime_times_vector (const int Nens,
 
     for (int n = 0; n < Nens; ++n)
     {
-        MultiFab xf_n =
-            read_member_multifab(n, last_pf_name, varnames);
+        MultiFab xf_n = read_member_multifab(n, last_pf_name, varnames);
 
         AMREX_ALWAYS_ASSERT(xf_n.boxArray() == result.boxArray());
         AMREX_ALWAYS_ASSERT(xf_n.DistributionMap() == result.DistributionMap());
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
index 932a1ed68a..519656f306 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
@@ -140,7 +140,7 @@ int main (int argc, char* argv[])
         // Move plotfiles (plt*) from current directory into member_dir/plotfiles
         for (auto &p : fs::directory_iterator(".")) {
             std::string fname = p.path().filename().string();
-            if (fname.find("plt") == 0) {  // starts with "plt"
+            if (fname.find("plt0") == 0) {  // starts with "plt"
                 fs::rename(p.path(), fs::path(member_dir) / "plotfiles" / fname);
             }
         }
@@ -159,7 +159,7 @@ int main (int argc, char* argv[])
 
    ERF tmp_erf;
    // This is only a post-processing step for visualization
-   tmp_erf.ComputeAndWriteEnsemblePerturbations();
+   //tmp_erf.ComputeAndWriteEnsemblePerturbations();
 
    // Perform data assimilation
    int da_iter = 0;

From f19deb2b986991b666b9f8940ce24d41cb1ea8e1 Mon Sep 17 00:00:00 2001
From: nataraj2 <mn556@cornell.edu>
Date: Tue, 18 Aug 2026 17:43:56 -0700
Subject: [PATCH 06/18] Avoiding moving of plt and chk files for ensembles

---
 .../SquallLine2D/ERF_DA_EnKFSRF.cpp           | 28 ++++++++++++
 .../DataAssimilation/SquallLine2D/main.cpp    | 43 ++++---------------
 Source/ERF.H                                  |  6 ++-
 3 files changed, 40 insertions(+), 37 deletions(-)

diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
index 1cba0ab13e..a5934f8655 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
@@ -363,3 +363,31 @@ ERF::PerformDataAssimilation(int da_iter)
         WriteCheckpointFile();
     }
 }
+
+
+void
+ERF::SetDirsForPlotfilesAndCheckpointsForDA (const int ens_no)
+{
+    // --------------------------------------------------------
+    // Set up member-specific output directories
+    // --------------------------------------------------------
+    std::string member_dir;
+
+    std::stringstream ss;
+    ss << "member_" << std::setw(2) << std::setfill('0') << ens_no;
+    member_dir = ss.str();
+
+    if (ParallelDescriptor::IOProcessor())
+    {
+        fs::create_directories(member_dir + "/plotfiles");
+        fs::create_directories(member_dir + "/chkfiles");
+        fs::create_directories(member_dir + "/pertfiles");
+    }
+
+    ParallelDescriptor::Barrier();
+
+    // Only change the basename/prefix.
+    // ERF will handle the step numbering as before.
+    check_file = member_dir + "/chkfiles/chk";
+    plot3d_file_1  = member_dir + "/plotfiles/plt";
+}
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
index 519656f306..818779f33f 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
@@ -100,13 +100,16 @@ int main (int argc, char* argv[])
     pp_ens.query("n_members", n_ens);
 
     // Ensemble run loop
-    for (int ie = 0; ie < n_ens; ++ie)
+    for (int ens_no = 0; ens_no < n_ens; ++ens_no)
     {
         // --------------------------------------------------------
         // Fresh ERF instance per ensemble
         // --------------------------------------------------------
         ERF erf;
 
+        if(is_init_for_ensemble) {
+            erf.SetDirsForPlotfilesAndCheckpointsForDA(ens_no);
+        }
         erf.InitData();
         erf.Evolve();
 
@@ -115,46 +118,16 @@ int main (int argc, char* argv[])
             end_total, ParallelDescriptor::IOProcessorNumber());
 
         if (erf.Verbose()) {
-            amrex::Print() << "Ensemble " << ie
+            amrex::Print() << "Ensemble " << ens_no
                            << " wallclock time: "
                            << end_total << '\n';
         }
 
         // --------------------------------------------------------
-    // MPI barrier to ensure all ranks finish Evolve
-    // --------------------------------------------------------
-    ParallelDescriptor::Barrier();
-
-    if (is_init_for_ensemble && ParallelDescriptor::IOProcessor())
-    {
-        // Create zero-padded member directory
-        std::stringstream ss;
-        ss << "member_" << std::setw(2) << std::setfill('0') << (ie);
-        std::string member_dir = ss.str();
-
-        fs::create_directory(member_dir);
-        fs::create_directory(member_dir + "/plotfiles");
-        fs::create_directory(member_dir + "/chkfiles");
-        fs::create_directory(member_dir + "/pertfiles");
-
-        // Move plotfiles (plt*) from current directory into member_dir/plotfiles
-        for (auto &p : fs::directory_iterator(".")) {
-            std::string fname = p.path().filename().string();
-            if (fname.find("plt0") == 0) {  // starts with "plt"
-                fs::rename(p.path(), fs::path(member_dir) / "plotfiles" / fname);
-            }
-        }
+        // MPI barrier to ensure all ranks finish Evolve
+        // --------------------------------------------------------
+        ParallelDescriptor::Barrier();
 
-        // Move checkpoint files (chk*) from current directory into member_dir/chkfiles
-        for (auto &c : fs::directory_iterator(".")) {
-            std::string fname = c.path().filename().string();
-            if (fname.find("chk") == 0) {  // starts with "chk"
-                fs::rename(c.path(), fs::path(member_dir) / "chkfiles" / fname);
-            }
-        }
-    }
-    // Optional: barrier after move to ensure rank 0 is done
-    ParallelDescriptor::Barrier();
    } // Ensemble run loop complete
 
    ERF tmp_erf;
diff --git a/Source/ERF.H b/Source/ERF.H
index 33b8987f33..c22a1cf779 100644
--- a/Source/ERF.H
+++ b/Source/ERF.H
@@ -660,9 +660,11 @@ public:
     void apply_gaussian_smoothing_to_perturbations (const int lev,
                                                     amrex::MultiFab& mf_cc_pert);
 
-    void ComputeAndWriteEnsemblePerturbations();
+    void ComputeAndWriteEnsemblePerturbations ();
 
-    void PerformDataAssimilation(int da_iter);
+    void PerformDataAssimilation (int da_iter);
+
+    void SetDirsForPlotfilesAndCheckpointsForDA (const int ens_no);
 
     void init_custom (int lev);
 

From 1d736ca7a238636be5bac09fd948b6a9146be25a Mon Sep 17 00:00:00 2001
From: nataraj2 <mn556@cornell.edu>
Date: Sun, 30 Aug 2026 14:15:53 -0700
Subject: [PATCH 07/18] Getting ready for DA cycles

---
 .../SquallLine2D/ERF_DA_EnKFSRF.cpp           | 61 ++++++------
 .../DataAssimilation/SquallLine2D/main.cpp    | 95 +++++++++++--------
 2 files changed, 88 insertions(+), 68 deletions(-)

diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
index a5934f8655..75e43d4d83 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
@@ -182,6 +182,33 @@ ERF::ComputeAndWriteEnsemblePerturbations()
     }
 }
 
+void
+ERF::SetDirsForPlotfilesAndCheckpointsForDA (const int ens_no)
+{
+    // --------------------------------------------------------
+    // Set up member-specific output directories
+    // --------------------------------------------------------
+    std::string member_dir;
+
+    std::stringstream ss;
+    ss << "member_" << std::setw(2) << std::setfill('0') << ens_no;
+    member_dir = ss.str();
+
+    if (ParallelDescriptor::IOProcessor())
+    {
+        fs::create_directories(member_dir + "/plotfiles");
+        fs::create_directories(member_dir + "/chkfiles");
+        fs::create_directories(member_dir + "/pertfiles");
+    }
+
+    ParallelDescriptor::Barrier();
+
+    // Only change the basename/prefix.
+    // ERF will handle the step numbering as before.
+    check_file = member_dir + "/chkfiles/chk";
+    plot3d_file_1  = member_dir + "/plotfiles/plt";
+}
+
 void
 ERF::PerformDataAssimilation(int da_iter)
 {
@@ -348,6 +375,12 @@ ERF::PerformDataAssimilation(int da_iter)
 
         m_plot3d_int_1 = -1;
         m_check_int = -1;
+        // This handles all the initialization including when doing a
+        // restart. But in the inputs file, we do not specify a
+        // erf.restart. So, the ERF class does not know it is a restart.
+        // So, whenever InitData() is called, it just does the from scratch
+        // initialization. All the restart has to be explicitly handled by
+        // specifying the restart_chkfile variable for the corresponding ensemble
         InitData();
         auto& lev_new = vars_new[0];
 
@@ -363,31 +396,3 @@ ERF::PerformDataAssimilation(int da_iter)
         WriteCheckpointFile();
     }
 }
-
-
-void
-ERF::SetDirsForPlotfilesAndCheckpointsForDA (const int ens_no)
-{
-    // --------------------------------------------------------
-    // Set up member-specific output directories
-    // --------------------------------------------------------
-    std::string member_dir;
-
-    std::stringstream ss;
-    ss << "member_" << std::setw(2) << std::setfill('0') << ens_no;
-    member_dir = ss.str();
-
-    if (ParallelDescriptor::IOProcessor())
-    {
-        fs::create_directories(member_dir + "/plotfiles");
-        fs::create_directories(member_dir + "/chkfiles");
-        fs::create_directories(member_dir + "/pertfiles");
-    }
-
-    ParallelDescriptor::Barrier();
-
-    // Only change the basename/prefix.
-    // ERF will handle the step numbering as before.
-    check_file = member_dir + "/chkfiles/chk";
-    plot3d_file_1  = member_dir + "/plotfiles/plt";
-}
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
index 818779f33f..a0fda86931 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
@@ -96,51 +96,66 @@ int main (int argc, char* argv[])
 
     ParmParse pp_ens("ensemble");
 
-    int n_ens = 1;
-    pp_ens.query("n_members", n_ens);
+    int num_da_cycles = -1;
+    pp_erf.query("num_da_cycles", num_da_cycles);
 
-    // Ensemble run loop
-    for (int ens_no = 0; ens_no < n_ens; ++ens_no)
-    {
-        // --------------------------------------------------------
-        // Fresh ERF instance per ensemble
-        // --------------------------------------------------------
-        ERF erf;
-
-        if(is_init_for_ensemble) {
-            erf.SetDirsForPlotfilesAndCheckpointsForDA(ens_no);
-        }
-        erf.InitData();
-        erf.Evolve();
-
-        Real end_total = amrex::second() - strt_total;
-        ParallelDescriptor::ReduceRealMax(
-            end_total, ParallelDescriptor::IOProcessorNumber());
-
-        if (erf.Verbose()) {
-            amrex::Print() << "Ensemble " << ens_no
-                           << " wallclock time: "
-                           << end_total << '\n';
-        }
-
-        // --------------------------------------------------------
-        // MPI barrier to ensure all ranks finish Evolve
-        // --------------------------------------------------------
-        ParallelDescriptor::Barrier();
-
-   } // Ensemble run loop complete
+    if(num_da_cycles != -1 and !is_init_for_ensemble) {
+        Abort("You are trying to run data assimilation by using the erf.num_da_cycles option. For this, you also "
+              "have to specify erf.is_init_for_ensemble=true, as ensembles have to be run for data assimilation");
+    }
 
-   ERF tmp_erf;
-   // This is only a post-processing step for visualization
-   //tmp_erf.ComputeAndWriteEnsemblePerturbations();
+    int n_ens = 1;
+    pp_ens.query("n_members", n_ens);
 
-   // Perform data assimilation
-   int da_iter = 0;
-   tmp_erf.PerformDataAssimilation(da_iter);
+    for (int da_cycle_no = 0; da_cycle_no < num_da_cycles; ++da_cycle_no) {
+        // Ensemble run loop
+        for (int ens_no = 0; ens_no < n_ens; ++ens_no)
+        {
+            // --------------------------------------------------------
+            // Fresh ERF instance per ensemble
+            // --------------------------------------------------------
+            ERF erf;
+
+            if(is_init_for_ensemble) {
+                erf.SetDirsForPlotfilesAndCheckpointsForDA(ens_no);
+            }
+            // InitData() handles all the initialization including when doing a
+            // restart. But in the inputs file, for DA runs, we do not specify a
+            // erf.restart. So, the ERF class does not know it is a restart.
+            // So, whenever InitData() is called, it just does the from-scratch
+            // initialization. All the restart has to be explicitly handled by
+            // specifying the restart_chkfile variable for the corresponding ensemble
+            erf.InitData();
+            erf.Evolve();
+
+            Real end_total = amrex::second() - strt_total;
+            ParallelDescriptor::ReduceRealMax(
+                end_total, ParallelDescriptor::IOProcessorNumber());
+
+            if (erf.Verbose()) {
+                amrex::Print() << "Ensemble " << ens_no
+                              << " wallclock time: "
+                               << end_total << '\n';
+            }
+
+            // --------------------------------------------------------
+            // MPI barrier to ensure all ranks finish Evolve
+            // --------------------------------------------------------
+            ParallelDescriptor::Barrier();
+        } // Ensemble run loop complete
+
+        ERF tmp_erf;
+        // This is only a post-processing step for visualization
+        //tmp_erf.ComputeAndWriteEnsemblePerturbations();
+
+        // After all ensemble runs are complete, perform data assimilation
+        int da_iter = 0;
+        tmp_erf.PerformDataAssimilation(da_iter);
+    }
 
-   BL_PROFILE_VAR_STOP(pmain);
+        BL_PROFILE_VAR_STOP(pmain);
 
-   amrex::Finalize();
+        amrex::Finalize();
 
 #ifdef AMREX_USE_MPI
     MPI_Finalize();

From 3d8906da995ecf9acba9d0d399408ae206db4996 Mon Sep 17 00:00:00 2001
From: nataraj2 <mn556@cornell.edu>
Date: Sun, 30 Aug 2026 14:23:58 -0700
Subject: [PATCH 08/18] Updating file

---
 .../SquallLine2D/ERF_DA_EnKFSRF.cpp           | 55 +++++++++----------
 1 file changed, 27 insertions(+), 28 deletions(-)

diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
index d20eedf77b..6aead7c3cc 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
@@ -209,6 +209,33 @@ ERF::SetDirsForPlotfilesAndCheckpointsForDA (const int ens_no)
     plot3d_file_1  = member_dir + "/plotfiles/plt";
 }
 
+void
+ERF::SetDirsForPlotfilesAndCheckpointsForDA (const int ens_no)
+{
+    // --------------------------------------------------------
+    // Set up member-specific output directories
+    // --------------------------------------------------------
+    std::string member_dir;
+
+    std::stringstream ss;
+    ss << "member_" << std::setw(2) << std::setfill('0') << ens_no;
+    member_dir = ss.str();
+
+    if (ParallelDescriptor::IOProcessor())
+    {
+        fs::create_directories(member_dir + "/plotfiles");
+        fs::create_directories(member_dir + "/chkfiles");
+        fs::create_directories(member_dir + "/pertfiles");
+    }
+
+    ParallelDescriptor::Barrier();
+
+    // Only change the basename/prefix.
+    // ERF will handle the step numbering as before.
+    check_file = member_dir + "/chkfiles/chk";
+    plot3d_file_1  = member_dir + "/plotfiles/plt";
+}
+
 void
 ERF::PerformDataAssimilation(int da_iter)
 {
@@ -396,31 +423,3 @@ ERF::PerformDataAssimilation(int da_iter)
         WriteCheckpointFile();
     }
 }
-
-
-void
-ERF::SetDirsForPlotfilesAndCheckpointsForDA (const int ens_no)
-{
-    // --------------------------------------------------------
-    // Set up member-specific output directories
-    // --------------------------------------------------------
-    std::string member_dir;
-
-    std::stringstream ss;
-    ss << "member_" << std::setw(2) << std::setfill('0') << ens_no;
-    member_dir = ss.str();
-
-    if (ParallelDescriptor::IOProcessor())
-    {
-        fs::create_directories(member_dir + "/plotfiles");
-        fs::create_directories(member_dir + "/chkfiles");
-        fs::create_directories(member_dir + "/pertfiles");
-    }
-
-    ParallelDescriptor::Barrier();
-
-    // Only change the basename/prefix.
-    // ERF will handle the step numbering as before.
-    check_file = member_dir + "/chkfiles/chk";
-    plot3d_file_1  = member_dir + "/plotfiles/plt";
-}

From 5826c7bb6fe3524d7e69cd6e3b6b114291d180c1 Mon Sep 17 00:00:00 2001
From: nataraj2 <mn556@cornell.edu>
Date: Sun, 30 Aug 2026 14:25:43 -0700
Subject: [PATCH 09/18] Delete duplicate function

---
 .../SquallLine2D/ERF_DA_EnKFSRF.cpp           | 27 -------------------
 1 file changed, 27 deletions(-)

diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
index 6aead7c3cc..75e43d4d83 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
@@ -209,33 +209,6 @@ ERF::SetDirsForPlotfilesAndCheckpointsForDA (const int ens_no)
     plot3d_file_1  = member_dir + "/plotfiles/plt";
 }
 
-void
-ERF::SetDirsForPlotfilesAndCheckpointsForDA (const int ens_no)
-{
-    // --------------------------------------------------------
-    // Set up member-specific output directories
-    // --------------------------------------------------------
-    std::string member_dir;
-
-    std::stringstream ss;
-    ss << "member_" << std::setw(2) << std::setfill('0') << ens_no;
-    member_dir = ss.str();
-
-    if (ParallelDescriptor::IOProcessor())
-    {
-        fs::create_directories(member_dir + "/plotfiles");
-        fs::create_directories(member_dir + "/chkfiles");
-        fs::create_directories(member_dir + "/pertfiles");
-    }
-
-    ParallelDescriptor::Barrier();
-
-    // Only change the basename/prefix.
-    // ERF will handle the step numbering as before.
-    check_file = member_dir + "/chkfiles/chk";
-    plot3d_file_1  = member_dir + "/plotfiles/plt";
-}
-
 void
 ERF::PerformDataAssimilation(int da_iter)
 {

From 22584892eefc0972ead9e4d4cfad2578b93a41c4 Mon Sep 17 00:00:00 2001
From: nataraj2 <mn556@cornell.edu>
Date: Sun, 30 Aug 2026 14:46:28 -0700
Subject: [PATCH 10/18] Some minor changes

---
 .../SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp     | 21 +++++++++++++++++++
 .../DataAssimilation/SquallLine2D/main.cpp    | 15 +++++++++++--
 Source/ERF.H                                  |  2 ++
 3 files changed, 36 insertions(+), 2 deletions(-)

diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
index 6b86485207..a7955ab033 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
@@ -639,3 +639,24 @@ MakeEnsembleCheckpointName (int da_iter, int ens_no)
 
     return chk_name.str();
 }
+
+void
+ERF::GetEnsembleCheckpointName (int da_iter, int ens_no)
+{
+    // Create CheckpointAfterDA if needed
+    const std::string base_dir = "CheckpointAfterDA";
+    fs::create_directories(base_dir);
+
+    // Create CheckpointAfterDA/DACycle_<da_iter>
+    const std::string da_dir =
+        base_dir + "/DACycle_" + std::to_string(da_iter);
+    fs::create_directories(da_dir);
+
+    std::ostringstream chk_name;
+    chk_name << da_dir
+             << "/chk_"
+             << std::setw(2) << std::setfill('0') << ens_no
+             << "_";
+
+    restart_chkfile = chk_name.str();
+}
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
index f78fd79425..118a493c2d 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
@@ -107,7 +107,8 @@ int main (int argc, char* argv[])
     int n_ens = 1;
     pp_ens.query("n_members", n_ens);
 
-    for (int da_cycle_no = 0; da_cycle_no < num_da_cycles; ++da_cycle_no) {
+    // Data assimilation cycle loop
+    for (int da_iter = 0; da_iter < num_da_cycles; ++da_iter) {
         // Ensemble run loop
         for (int ens_no = 0; ens_no < n_ens; ++ens_no)
         {
@@ -125,6 +126,17 @@ int main (int argc, char* argv[])
             // So, whenever InitData() is called, it just does the from-scratch
             // initialization. All the restart has to be explicitly handled by
             // specifying the restart_chkfile variable for the corresponding ensemble
+
+
+            // So, for da_iter=0, it starts from scratch, and for iteration da_iter>=1
+            // there has to be a erf.restart_chkfile specified, which is the checkpoint file
+            // which contains the latest updated ensemble after data assimilation for that ensemble.
+            // The InitData reads that restart file
+
+            // This call will fill the restart_chkfile string
+            if(da_iter > 0){
+                erf.GetEnsembleCheckpointName(da_iter-1, ens_no);
+            }
             erf.InitData();
             erf.Evolve();
 
@@ -149,7 +161,6 @@ int main (int argc, char* argv[])
         //tmp_erf.ComputeAndWriteEnsemblePerturbations();
 
         // After all ensemble runs are complete, perform data assimilation
-        int da_iter = 0;
         tmp_erf.PerformDataAssimilation(da_iter);
     }
 
diff --git a/Source/ERF.H b/Source/ERF.H
index d98972d489..490ae159d0 100644
--- a/Source/ERF.H
+++ b/Source/ERF.H
@@ -714,6 +714,8 @@ public:
 
     void SetDirsForPlotfilesAndCheckpointsForDA (const int ens_no);
 
+    void GetEnsembleCheckpointName (int da_iter, int ens_no);
+
     void init_custom (int lev);
 
 #ifdef ERF_USE_MULTIBLOCK

From 6f5532c5af1b909b90ab0477018a48ce0e190033 Mon Sep 17 00:00:00 2001
From: nataraj2 <mn556@cornell.edu>
Date: Mon, 31 Aug 2026 15:46:35 -0700
Subject: [PATCH 11/18] Some restructuring

---
 .../SquallLine2D/ERF_DA_EnKFSRF.cpp           | 27 -----------------
 .../SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp     | 29 ++++++++++++++++++-
 2 files changed, 28 insertions(+), 28 deletions(-)

diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
index 75e43d4d83..6e4d4950ea 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
@@ -182,33 +182,6 @@ ERF::ComputeAndWriteEnsemblePerturbations()
     }
 }
 
-void
-ERF::SetDirsForPlotfilesAndCheckpointsForDA (const int ens_no)
-{
-    // --------------------------------------------------------
-    // Set up member-specific output directories
-    // --------------------------------------------------------
-    std::string member_dir;
-
-    std::stringstream ss;
-    ss << "member_" << std::setw(2) << std::setfill('0') << ens_no;
-    member_dir = ss.str();
-
-    if (ParallelDescriptor::IOProcessor())
-    {
-        fs::create_directories(member_dir + "/plotfiles");
-        fs::create_directories(member_dir + "/chkfiles");
-        fs::create_directories(member_dir + "/pertfiles");
-    }
-
-    ParallelDescriptor::Barrier();
-
-    // Only change the basename/prefix.
-    // ERF will handle the step numbering as before.
-    check_file = member_dir + "/chkfiles/chk";
-    plot3d_file_1  = member_dir + "/plotfiles/plt";
-}
-
 void
 ERF::PerformDataAssimilation(int da_iter)
 {
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
index a7955ab033..cfc905ebd1 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
@@ -616,6 +616,33 @@ update_ensemble (const int Nens,
     compute_Xf_prime_times_vector(Nens, last_pf_name, varnames, xf_bar, T_colvec, result);
 }
 
+void
+ERF::SetDirsForPlotfilesAndCheckpointsForDA (const int ens_no)
+{
+    // --------------------------------------------------------
+    // Set up member-specific output directories
+    // --------------------------------------------------------
+    std::string member_dir;
+
+    std::stringstream ss;
+    ss << "member_" << std::setw(2) << std::setfill('0') << ens_no;
+    member_dir = ss.str();
+
+    if (ParallelDescriptor::IOProcessor())
+    {
+        fs::create_directories(member_dir + "/plotfiles");
+        fs::create_directories(member_dir + "/chkfiles");
+        fs::create_directories(member_dir + "/pertfiles");
+    }
+
+    ParallelDescriptor::Barrier();
+
+    // Only change the basename/prefix.
+    // ERF will handle the step numbering as before.
+    check_file = member_dir + "/chkfiles/chk";
+    plot3d_file_1  = member_dir + "/plotfiles/plt";
+}
+
 // da_iter : DA iteration number
 // ens_no  : Ensemble number
 
@@ -656,7 +683,7 @@ ERF::GetEnsembleCheckpointName (int da_iter, int ens_no)
     chk_name << da_dir
              << "/chk_"
              << std::setw(2) << std::setfill('0') << ens_no
-             << "_";
+             << "_" <<  "00000";
 
     restart_chkfile = chk_name.str();
 }

From a115ecffddad16d393447e7e1967e18be826fff0 Mon Sep 17 00:00:00 2001
From: nataraj2 <mn556@cornell.edu>
Date: Mon, 31 Aug 2026 16:04:33 -0700
Subject: [PATCH 12/18] Restructuring directories to be more intiuitive

---
 .../SquallLine2D/ERF_DA_EnKFSRF.H             |  3 ++
 .../SquallLine2D/ERF_DA_EnKFSRF.cpp           | 10 ++---
 .../SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp     | 44 +++++++++++--------
 .../DataAssimilation/SquallLine2D/main.cpp    |  2 +-
 Source/ERF.H                                  |  2 +-
 5 files changed, 36 insertions(+), 25 deletions(-)

diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.H b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.H
index 40bc5cec2d..050db1f373 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.H
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.H
@@ -208,4 +208,7 @@ void compute_S_matrix(Matrix& S,
                       const amrex::Vector<std::string>& varnames);
 std::string
 MakeEnsembleCheckpointName (int da_iter, int ens_no);
+
+std::string
+MakeFullPath (const std::string& filename, const int da_iter);
 #endif
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
index 6e4d4950ea..349aa68722 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
@@ -210,7 +210,7 @@ ERF::PerformDataAssimilation(int da_iter)
     Print() << "Computing ensemble mean complete" << std::endl;
 
     // Construct perturbation plotfile name
-    std::string pltname = "plt_ens_mean";
+    std::string pltname = MakeFullPath("plt_ens_mean", da_iter);
     WriteSingleLevelPlotfile(pltname,
                              xf_bar,
                              varnames,
@@ -226,7 +226,7 @@ ERF::PerformDataAssimilation(int da_iter)
 
     Vector<std::string> varnames1 = {"density","theta", "x_velocity","y_velocity","z_velocity", "qv", "qc", "qrain"};
     // Construct perturbation plotfile name
-    std::string pltname1 = "plt_mean_H_xf";
+    std::string pltname1 = MakeFullPath("plt_mean_H_xf", da_iter);
     WriteSingleLevelPlotfile(pltname1,
                              mean_H_xf,
                              varnames1,
@@ -240,7 +240,7 @@ ERF::PerformDataAssimilation(int da_iter)
 
     Print() << "Observation reading complete" << std::endl;
 
-    /*std::string pltname2 = "plt_y_obs";
+    /*std::string pltname2 = WriteSingleLevelPlotfile("plt_y_obs", da_iter);
     Vector<std::string> varnames2 = {"x_velocity", "y_velocity"};
     WriteSingleLevelPlotfile(pltname2,
                             y_obs,
@@ -299,7 +299,7 @@ ERF::PerformDataAssimilation(int da_iter)
 
     Print() << "compute_Xf_prime_times_vector complete" << std::endl;
 
-    /*std::string pltname3 = "plt_Xf_prime_alpha";
+    /*std::string pltname3 = WriteSingleLevelPlotfile("plt_Xf_prime_alpha", da_iter);
     Vector<std::string> varnames3 = {"density", "theta", "x_velocity", "y_velocity", "z_velocity"};
     WriteSingleLevelPlotfile(pltname3,
                             Xf_prime_alpha,
@@ -314,7 +314,7 @@ ERF::PerformDataAssimilation(int da_iter)
 
     Print() << "add_multifabs(xf_bar, Xf_prime_alpha, xf_bar_updated) complete" << std::endl;
 
-    /*std::string pltname4 = "plt_xf_bar_updated";
+    /*std::string pltname4 = WriteSingleLevelPlotfile("plt_xf_bar_updated", da_iter);
     Vector<std::string> varnames4 = {"density", "theta", "x_velocity", "y_velocity", "z_velocity"};
     WriteSingleLevelPlotfile(pltname4,
                             xf_bar_updated,
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
index cfc905ebd1..ccf20ae659 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
@@ -617,16 +617,19 @@ update_ensemble (const int Nens,
 }
 
 void
-ERF::SetDirsForPlotfilesAndCheckpointsForDA (const int ens_no)
+ERF::SetDirsForPlotfilesAndCheckpointsForDA (const int da_iter, 
+                                             const int ens_no)
 {
     // --------------------------------------------------------
-    // Set up member-specific output directories
+    // Set up cycle- and member-specific output directories
     // --------------------------------------------------------
-    std::string member_dir;
+    std::stringstream cycle_ss;
+    cycle_ss << "DA_Cycle_" << da_iter;
+    std::string cycle_dir = cycle_ss.str();
 
-    std::stringstream ss;
-    ss << "member_" << std::setw(2) << std::setfill('0') << ens_no;
-    member_dir = ss.str();
+    std::stringstream member_ss;
+    member_ss << "member_" << std::setw(2) << std::setfill('0') << ens_no;
+    std::string member_dir = cycle_dir + "/" + member_ss.str();
 
     if (ParallelDescriptor::IOProcessor())
     {
@@ -637,12 +640,10 @@ ERF::SetDirsForPlotfilesAndCheckpointsForDA (const int ens_no)
 
     ParallelDescriptor::Barrier();
 
-    // Only change the basename/prefix.
-    // ERF will handle the step numbering as before.
-    check_file = member_dir + "/chkfiles/chk";
-    plot3d_file_1  = member_dir + "/plotfiles/plt";
+    // Set file prefixes under the specific cycle and member path
+    check_file   = member_dir + "/chkfiles/chk";
+    plot3d_file_1 = member_dir + "/plotfiles/plt";
 }
-
 // da_iter : DA iteration number
 // ens_no  : Ensemble number
 
@@ -650,12 +651,12 @@ std::string
 MakeEnsembleCheckpointName (int da_iter, int ens_no)
 {
     // Create CheckpointAfterDA if needed
-    const std::string base_dir = "CheckpointAfterDA";
+    const std::string base_dir = "/DACycle_" + std::to_string(da_iter);
     fs::create_directories(base_dir);
 
-    // Create CheckpointAfterDA/DACycle_<da_iter>
-    const std::string da_dir =
-        base_dir + "/DACycle_" + std::to_string(da_iter);
+    const std::string da_dir = base_dir + "CheckpointAfterDA";
+
+    // Create DACycle_<da_iter>/CheckpointAfterDA
     fs::create_directories(da_dir);
 
     std::ostringstream chk_name;
@@ -671,12 +672,11 @@ void
 ERF::GetEnsembleCheckpointName (int da_iter, int ens_no)
 {
     // Create CheckpointAfterDA if needed
-    const std::string base_dir = "CheckpointAfterDA";
+    const std::string base_dir = "/DACycle_" + std::to_string(da_iter);
     fs::create_directories(base_dir);
 
     // Create CheckpointAfterDA/DACycle_<da_iter>
-    const std::string da_dir =
-        base_dir + "/DACycle_" + std::to_string(da_iter);
+    const std::string da_dir = base_dir + "CheckpointAfterDA";
     fs::create_directories(da_dir);
 
     std::ostringstream chk_name;
@@ -687,3 +687,11 @@ ERF::GetEnsembleCheckpointName (int da_iter, int ens_no)
 
     restart_chkfile = chk_name.str();
 }
+
+std::string
+MakeFullPath (const std::string& filename, const int da_iter)
+{
+    std::stringstream ss;
+    ss << "DA_Cycle_" << da_iter << "/" << filename;
+    return ss.str();
+}
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
index 118a493c2d..747be72307 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
@@ -118,7 +118,7 @@ int main (int argc, char* argv[])
             ERF erf;
 
             if(is_init_for_ensemble) {
-                erf.SetDirsForPlotfilesAndCheckpointsForDA(ens_no);
+                erf.SetDirsForPlotfilesAndCheckpointsForDA(da_iter, ens_no);
             }
             // InitData() handles all the initialization including when doing a
             // restart. But in the inputs file, for DA runs, we do not specify a
diff --git a/Source/ERF.H b/Source/ERF.H
index 490ae159d0..d458086133 100644
--- a/Source/ERF.H
+++ b/Source/ERF.H
@@ -712,7 +712,7 @@ public:
 
     void PerformDataAssimilation (int da_iter);
 
-    void SetDirsForPlotfilesAndCheckpointsForDA (const int ens_no);
+    void SetDirsForPlotfilesAndCheckpointsForDA (const int da_iter, const int ens_no);
 
     void GetEnsembleCheckpointName (int da_iter, int ens_no);
 

From c9dbe30409196d57d9d1f07e4ba981a7c3af0be7 Mon Sep 17 00:00:00 2001
From: nataraj2 <mn556@cornell.edu>
Date: Mon, 31 Aug 2026 16:43:21 -0700
Subject: [PATCH 13/18] Restructuring code for new directory hierarchy

---
 .../SquallLine2D/ERF_DA_EnKFSRF.H             | 33 ++++---
 .../SquallLine2D/ERF_DA_EnKFSRF.cpp           | 22 ++---
 .../SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp     | 97 +++++++++++--------
 .../DataAssimilation/SquallLine2D/main.cpp    |  2 +-
 Source/ERF.H                                  |  2 +-
 5 files changed, 88 insertions(+), 68 deletions(-)

diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.H b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.H
index 050db1f373..408dfd38f4 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.H
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.H
@@ -129,12 +129,14 @@ read_plot_file(amrex::PlotFileData& pf,
                amrex::MultiFab& mf);
 
 amrex::MultiFab
-read_member_multifab(int n,
+read_member_multifab(const int da_iter, 
+                     const int n,
                      const std::string& pf_name,
                      const amrex::Vector<std::string>& varnames);
 
 amrex::MultiFab
-compute_ensemble_mean(int Nens,
+compute_ensemble_mean(const int da_iter, 
+                      const int Nens,
                       const std::string& pf_name,
                       const amrex::Vector<std::string>& varnames);
 
@@ -142,7 +144,7 @@ void Apply_H(const amrex::MultiFab& x_mf, amrex::MultiFab& y_mf);
 
 void read_in_observations(const int& da_iter, amrex::MultiFab& y_obs);
 
-std::vector<std::string> get_plotfile_list();
+std::vector<std::string> get_plotfile_list(const int da_iter);
 void read_in_observations(const int& da_iter,
                           const amrex::Vector<std::string>& varnames,
                           amrex::MultiFab& y_obs);
@@ -156,12 +158,14 @@ void compute_d_prime_vec(amrex::MultiFab& d_prime_vec,
                     const amrex::MultiFab& d_vec,
                     const amrex::Vector<amrex::Real>& R_diag);
 
-void compute_mean_H_xf(amrex::MultiFab& mean_H_xf,
+void compute_mean_H_xf(const int da_iter,
+                       amrex::MultiFab& mean_H_xf,
                        const int Nens,
                        const std::string& last_pf_name,
                        const amrex::Vector<std::string>& varnames);
 
-void compute_r_vec (int Nens,
+void compute_r_vec (const int da_iter, 
+                    int Nens,
                     const std::string& last_pf_name,
                     const amrex::Vector<std::string>& varnames,
                     const amrex::MultiFab& mean_H_xf,
@@ -173,14 +177,16 @@ void compute_alpha_vec (const int& Nens,
                         const amrex::Vector<amrex::Real>& r_vec,
                         amrex::Vector<amrex::Real>& alpha_vec);
 
-void compute_Xf_prime_times_vector (const int Nens,
-                            const std::string& last_pf_name,
-                            const amrex::Vector<std::string>& varnames,
-                            const amrex::MultiFab& xf_bar,
-                            const amrex::Vector<amrex::Real>& vec_in,
-                            amrex::MultiFab& result);
+void compute_Xf_prime_times_vector (const int da_iter,
+                                    const int Nens,
+                                    const std::string& last_pf_name,
+                                    const amrex::Vector<std::string>& varnames,
+                                    const amrex::MultiFab& xf_bar,
+                                    const amrex::Vector<amrex::Real>& vec_in,
+                                    amrex::MultiFab& result);
 
-void update_ensemble (const int Nens,
+void update_ensemble (const int da_iter, 
+                      const int Nens,
                       const std::string& last_pf_name,
                       const amrex::Vector<std::string>& varnames,
                       const amrex::MultiFab& xf_bar,
@@ -200,7 +206,8 @@ Matrix matrix_multiply(const Matrix& A, const Matrix& B);
 void matrix_print(const Matrix& A);
 void lapack_testing();
 
-void compute_S_matrix(Matrix& S,
+void compute_S_matrix(const int da_iter, 
+                      Matrix& S,
                       const int& Nens,
                       const amrex::MultiFab& mean_H_xf,
                       const amrex::Vector<amrex::Real>& R_diag,
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
index 349aa68722..dd7d0cf385 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
@@ -140,10 +140,10 @@ WriteUpdatedEnsembleToERFClassData (const MultiFab& mf_cc_fine,
 }
 
 void
-ERF::ComputeAndWriteEnsemblePerturbations()
+ERF::ComputeAndWriteEnsemblePerturbations(const int da_iter)
 {
 
-    auto pltfiles = get_plotfile_list();
+    auto pltfiles = get_plotfile_list(da_iter);
 
     // Step 2: loop over all plotfiles (timestamps at which the plotfiles are written)
     // ie.find the ensemble mean at iteration 100, loop over the plt00100 file in each of the
@@ -153,13 +153,13 @@ ERF::ComputeAndWriteEnsemblePerturbations()
     const std::string member_prefix = "member_";
     for (const auto& pf_name : pltfiles)
     {
-        MultiFab mf_mean = compute_ensemble_mean(Nens, pf_name, varnames);
+        MultiFab mf_mean = compute_ensemble_mean(da_iter, Nens, pf_name, varnames);
         // -------------------------------
         // Step 3 & 4: Compute perturbations and write plotfiles
         // -------------------------------
         for (int n = 0; n < Nens; ++n)
         {
-            MultiFab mf_pert = read_member_multifab(n, pf_name, varnames);
+            MultiFab mf_pert = read_member_multifab(da_iter, n, pf_name, varnames);
             MultiFab::Subtract(mf_pert, mf_mean, 0, 0, mf_mean.nComp(), mf_mean.nGrow());
 
             // Create output directory
@@ -187,7 +187,7 @@ ERF::PerformDataAssimilation(int da_iter)
 {
     //lapack_testing();
 
-    auto pltfiles = get_plotfile_list();
+    auto pltfiles = get_plotfile_list(da_iter);
     std::string last_pf_name;
     if (!pltfiles.empty()) {
         last_pf_name = pltfiles.back();
@@ -205,7 +205,7 @@ ERF::PerformDataAssimilation(int da_iter)
     Vector<std::string> varnames = {"density","theta", "x_velocity","y_velocity","z_velocity", "qv", "qc", "qrain"};
 
     // Compute the ensemble mean
-    MultiFab xf_bar = compute_ensemble_mean(Nens, last_pf_name, varnames);
+    MultiFab xf_bar = compute_ensemble_mean(da_iter, Nens, last_pf_name, varnames);
 
     Print() << "Computing ensemble mean complete" << std::endl;
 
@@ -220,7 +220,7 @@ ERF::PerformDataAssimilation(int da_iter)
 
     // Compute the mean of forecast observations yf_bar = Hx_f
     MultiFab mean_H_xf;
-    compute_mean_H_xf(mean_H_xf, Nens, last_pf_name, varnames);
+    compute_mean_H_xf(da_iter, mean_H_xf, Nens, last_pf_name, varnames);
 
     Print() << "compute_mean_H_xf complete" << std::endl;
 
@@ -279,13 +279,13 @@ ERF::PerformDataAssimilation(int da_iter)
 
     // Compute r = Y'^Td'
     Vector<Real> r_vec;
-    compute_r_vec(Nens, last_pf_name, varnames, mean_H_xf, d_prime_vec, r_vec);
+    compute_r_vec(da_iter, Nens, last_pf_name, varnames, mean_H_xf, d_prime_vec, r_vec);
 
     Print() << "compute_r_vec complete" << std::endl;
 
     // Compute the S matrix
     Matrix S_mat(Nens);
-    compute_S_matrix(S_mat, Nens, mean_H_xf, R_diag, last_pf_name, varnames);
+    compute_S_matrix(da_iter, S_mat, Nens, mean_H_xf, R_diag, last_pf_name, varnames);
 
     Print() << "compute_S_matrix complete" << std::endl;
 
@@ -295,7 +295,7 @@ ERF::PerformDataAssimilation(int da_iter)
     Print() << "compute_alpha_vec complete" << std::endl;
 
     MultiFab Xf_prime_alpha;
-    compute_Xf_prime_times_vector(Nens, last_pf_name, varnames, xf_bar, alpha_vec, Xf_prime_alpha);
+    compute_Xf_prime_times_vector(da_iter, Nens, last_pf_name, varnames, xf_bar, alpha_vec, Xf_prime_alpha);
 
     Print() << "compute_Xf_prime_times_vector complete" << std::endl;
 
@@ -332,7 +332,7 @@ ERF::PerformDataAssimilation(int da_iter)
     for(int n=0; n< Nens; n++) {
         Print() << "Updating for ensemble " << n << std::endl;
         MultiFab mf_ens_pert;
-        update_ensemble(Nens, last_pf_name, varnames, xf_bar, T_mat, n, mf_ens_pert);
+        update_ensemble(da_iter, Nens, last_pf_name, varnames, xf_bar, T_mat, n, mf_ens_pert);
 
         // Update the ensemble
         MultiFab mf_ens_updated;
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
index ccf20ae659..9bb1afb789 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
@@ -82,12 +82,13 @@ lapack_testing()
     }
 }
 
-
 std::vector<std::string>
-get_plotfile_list()
+get_plotfile_list(const int da_iter)
 {
+    const std::string da_dir = "DA_Cycle_" + std::to_string(da_iter) + "/";
+    
     std::vector<std::string> pltfiles;
-    const std::string member_prefix = "member_";
+    const std::string member_prefix = da_dir + "member_";
 
     std::string pf_dir = member_prefix + "00/plotfiles";
 
@@ -159,11 +160,15 @@ read_plot_file(PlotFileData& pf,
 }
 
 MultiFab
-read_member_multifab(int n,
+read_member_multifab(const int da_iter, 
+                     const int n,
                      const std::string& pf_name,
                      const Vector<std::string>& varnames)
 {
-    const std::string member_prefix = "member_";
+    
+    const std::string da_dir = "DA_Cycle_" + std::to_string(da_iter) + "/";    
+
+    const std::string member_prefix = da_dir + "member_";
     std::string member_dir = member_prefix + amrex::Concatenate("", n, 2);
     std::string pf_path = member_dir + "/plotfiles/" + pf_name;
 
@@ -181,7 +186,8 @@ read_member_multifab(int n,
 }
 
 MultiFab
-compute_ensemble_mean(int Nens,
+compute_ensemble_mean(const int da_iter, 
+                      const int Nens,
                       const std::string& pf_name,
                       const Vector<std::string>& varnames)
 {
@@ -191,7 +197,7 @@ compute_ensemble_mean(int Nens,
     for (int n = 0; n < Nens; ++n)
     {
         std::cout << "Reading ensembles for compute_ensemble_mean = " << n << std::endl;
-        MultiFab mf_tmp = read_member_multifab(n, pf_name, varnames);
+        MultiFab mf_tmp = read_member_multifab(da_iter, n, pf_name, varnames);
 
         if (!initialized) {
             mf_mean.define(mf_tmp.boxArray(),
@@ -270,7 +276,8 @@ compute_d_prime_vec(MultiFab& d_prime_vec,
 }
 
 void
-compute_mean_H_xf(MultiFab& mean_H_xf,
+compute_mean_H_xf(const int da_iter,
+                  MultiFab& mean_H_xf,
                   const int Nens,
                   const std::string& last_pf_name,
                   const Vector<std::string>& varnames)
@@ -278,7 +285,7 @@ compute_mean_H_xf(MultiFab& mean_H_xf,
 
     for (int n = 0; n < Nens; ++n)
     {
-        MultiFab xf_i = read_member_multifab(n, last_pf_name, varnames);
+        MultiFab xf_i = read_member_multifab(da_iter, n, last_pf_name, varnames);
         if (n==0) {
             mean_H_xf.define(xf_i.boxArray(),
                              xf_i.DistributionMap(),
@@ -346,14 +353,15 @@ Compute_yf_prime_i_T_Rinv_yf_prime_j(const MultiFab& Yi,
 }
 
 void
-compute_yf_prime (int i,
-                const std::string& last_pf_name,
-                const Vector<std::string>& varnames,
-                const MultiFab& mean_H_xf,
-                MultiFab& yf_prime_i)
+compute_yf_prime (const int da_iter,
+                  const int i,
+                  const std::string& last_pf_name,
+                  const Vector<std::string>& varnames,
+                  const MultiFab& mean_H_xf,
+                  MultiFab& yf_prime_i)
 {
     // Read ensemble member
-    MultiFab xf_i = read_member_multifab(i, last_pf_name, varnames);
+    MultiFab xf_i = read_member_multifab(da_iter, i, last_pf_name, varnames);
 
     // Apply observation operator
     MultiFab H_xf_i;
@@ -380,7 +388,8 @@ compute_yf_prime (int i,
 }
 
 void
-compute_S_matrix(Matrix& S,
+compute_S_matrix(const int da_iter, 
+                 Matrix& S,
                  const int& Nens,
                  const MultiFab& mean_H_xf,
                  const Vector<Real>& R_diag,
@@ -389,11 +398,11 @@ compute_S_matrix(Matrix& S,
 {
     for (int i = 0; i < Nens; ++i) {
         MultiFab yf_prime_i;
-        compute_yf_prime(i, last_pf_name, varnames, mean_H_xf, yf_prime_i);
+        compute_yf_prime(da_iter, i, last_pf_name, varnames, mean_H_xf, yf_prime_i);
 
         for (int j = 0; j < Nens; ++j) {
             MultiFab yf_prime_j;
-            compute_yf_prime(j, last_pf_name, varnames, mean_H_xf, yf_prime_j);
+            compute_yf_prime(da_iter, j, last_pf_name, varnames, mean_H_xf, yf_prime_j);
             Real val = Compute_yf_prime_i_T_Rinv_yf_prime_j(yf_prime_i, yf_prime_j, R_diag);
             S(i,j) = val/(Nens-1);
             if(i==j) {
@@ -459,7 +468,8 @@ compute_yf_prime_T_d_prime_vec (const MultiFab& yf_prime,
 }
 
 void
-compute_r_vec (int Nens,
+compute_r_vec (const int da_iter, 
+               int Nens,
                const std::string& last_pf_name,
                const Vector<std::string>& varnames,
                const MultiFab& mean_H_xf,
@@ -472,7 +482,7 @@ compute_r_vec (int Nens,
     {
         MultiFab yf_prime_i;
 
-        compute_yf_prime(i, last_pf_name, varnames, mean_H_xf, yf_prime_i);
+        compute_yf_prime(da_iter, i, last_pf_name, varnames, mean_H_xf, yf_prime_i);
 
         r_vec[i] = compute_yf_prime_T_d_prime_vec(yf_prime_i, d_prime_vec);
     }
@@ -510,18 +520,20 @@ compute_alpha_vec (const int& Nens,
 // column-vector element multiply
 
 void
-compute_Xf_prime_times_vector (const int Nens,
-                            const std::string& last_pf_name,
-                            const Vector<std::string>& varnames,
-                            const MultiFab& xf_bar,
-                            const Vector<Real>& vec_in,
-                            MultiFab& result)
+compute_Xf_prime_times_vector (const int da_iter, 
+                               const int Nens,
+                               const std::string& last_pf_name,
+                               const Vector<std::string>& varnames,
+                               const MultiFab& xf_bar,
+                               const Vector<Real>& vec_in,
+                               MultiFab& result)
 {
     AMREX_ALWAYS_ASSERT(vec_in.size() == Nens);
 
     // Read first member to define result
     MultiFab xf_0 =
-        read_member_multifab(0,
+        read_member_multifab(da_iter, 
+                             0,
                              last_pf_name,
                              varnames);
 
@@ -534,7 +546,7 @@ compute_Xf_prime_times_vector (const int Nens,
 
     for (int n = 0; n < Nens; ++n)
     {
-        MultiFab xf_n = read_member_multifab(n, last_pf_name, varnames);
+        MultiFab xf_n = read_member_multifab(da_iter, n, last_pf_name, varnames);
 
         AMREX_ALWAYS_ASSERT(xf_n.boxArray() == result.boxArray());
         AMREX_ALWAYS_ASSERT(xf_n.DistributionMap() == result.DistributionMap());
@@ -601,7 +613,8 @@ compute_T_matrix (const Matrix& S_mat,
 
 
 void
-update_ensemble (const int Nens,
+update_ensemble (const int da_iter, 
+                 const int Nens,
                  const std::string& last_pf_name,
                  const Vector<std::string>& varnames,
                  const MultiFab& xf_bar,
@@ -613,7 +626,7 @@ update_ensemble (const int Nens,
     for (int i = 0; i < Nens; ++i) {
         T_colvec[i] = T(i,n);
     }
-    compute_Xf_prime_times_vector(Nens, last_pf_name, varnames, xf_bar, T_colvec, result);
+    compute_Xf_prime_times_vector(da_iter, Nens, last_pf_name, varnames, xf_bar, T_colvec, result);
 }
 
 void
@@ -651,16 +664,16 @@ std::string
 MakeEnsembleCheckpointName (int da_iter, int ens_no)
 {
     // Create CheckpointAfterDA if needed
-    const std::string base_dir = "/DACycle_" + std::to_string(da_iter);
-    fs::create_directories(base_dir);
+    const std::string da_dir = "DA_Cycle_" + std::to_string(da_iter) + "/";
+    fs::create_directories(da_dir);
 
-    const std::string da_dir = base_dir + "CheckpointAfterDA";
+    const std::string chk_after_da_dir = da_dir + "CheckpointAfterDA";
 
-    // Create DACycle_<da_iter>/CheckpointAfterDA
-    fs::create_directories(da_dir);
+    // Create DA_Cycle_<da_iter>/CheckpointAfterDA
+    fs::create_directories(chk_after_da_dir);
 
     std::ostringstream chk_name;
-    chk_name << da_dir
+    chk_name << chk_after_da_dir
              << "/chk_"
              << std::setw(2) << std::setfill('0') << ens_no
              << "_";
@@ -672,15 +685,15 @@ void
 ERF::GetEnsembleCheckpointName (int da_iter, int ens_no)
 {
     // Create CheckpointAfterDA if needed
-    const std::string base_dir = "/DACycle_" + std::to_string(da_iter);
-    fs::create_directories(base_dir);
-
-    // Create CheckpointAfterDA/DACycle_<da_iter>
-    const std::string da_dir = base_dir + "CheckpointAfterDA";
+    const std::string da_dir = "DA_Cycle_" + std::to_string(da_iter) + "/";
     fs::create_directories(da_dir);
 
+    // Create CheckpointAfterDA/DA_Cycle_<da_iter>
+    const std::string chk_after_da_dir = da_dir + "CheckpointAfterDA";
+    fs::create_directories(chk_after_da_dir);
+
     std::ostringstream chk_name;
-    chk_name << da_dir
+    chk_name << chk_after_da_dir
              << "/chk_"
              << std::setw(2) << std::setfill('0') << ens_no
              << "_" <<  "00000";
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
index 747be72307..625f94447b 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
@@ -158,7 +158,7 @@ int main (int argc, char* argv[])
 
         ERF tmp_erf;
         // This is only a post-processing step for visualization
-        //tmp_erf.ComputeAndWriteEnsemblePerturbations();
+        //tmp_erf.ComputeAndWriteEnsemblePerturbations(da_iter);
 
         // After all ensemble runs are complete, perform data assimilation
         tmp_erf.PerformDataAssimilation(da_iter);
diff --git a/Source/ERF.H b/Source/ERF.H
index d458086133..fc35cd4853 100644
--- a/Source/ERF.H
+++ b/Source/ERF.H
@@ -708,7 +708,7 @@ public:
     void apply_gaussian_smoothing_to_perturbations (const int lev,
                                                     amrex::MultiFab& mf_cc_pert);
 
-    void ComputeAndWriteEnsemblePerturbations ();
+    void ComputeAndWriteEnsemblePerturbations (const int da_iter);
 
     void PerformDataAssimilation (int da_iter);
 

From 5998ab32918e356c19c7855a5f63e2f6d920ca7a Mon Sep 17 00:00:00 2001
From: nataraj2 <mn556@cornell.edu>
Date: Mon, 31 Aug 2026 16:44:41 -0700
Subject: [PATCH 14/18] Correcting tabs and trailing whitespaces

---
 .../SquallLine2D/ERF_DA_EnKFSRF.H             | 10 ++++-----
 .../SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp     | 22 +++++++++----------
 2 files changed, 16 insertions(+), 16 deletions(-)

diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.H b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.H
index 408dfd38f4..27a4550194 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.H
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.H
@@ -129,13 +129,13 @@ read_plot_file(amrex::PlotFileData& pf,
                amrex::MultiFab& mf);
 
 amrex::MultiFab
-read_member_multifab(const int da_iter, 
+read_member_multifab(const int da_iter,
                      const int n,
                      const std::string& pf_name,
                      const amrex::Vector<std::string>& varnames);
 
 amrex::MultiFab
-compute_ensemble_mean(const int da_iter, 
+compute_ensemble_mean(const int da_iter,
                       const int Nens,
                       const std::string& pf_name,
                       const amrex::Vector<std::string>& varnames);
@@ -164,7 +164,7 @@ void compute_mean_H_xf(const int da_iter,
                        const std::string& last_pf_name,
                        const amrex::Vector<std::string>& varnames);
 
-void compute_r_vec (const int da_iter, 
+void compute_r_vec (const int da_iter,
                     int Nens,
                     const std::string& last_pf_name,
                     const amrex::Vector<std::string>& varnames,
@@ -185,7 +185,7 @@ void compute_Xf_prime_times_vector (const int da_iter,
                                     const amrex::Vector<amrex::Real>& vec_in,
                                     amrex::MultiFab& result);
 
-void update_ensemble (const int da_iter, 
+void update_ensemble (const int da_iter,
                       const int Nens,
                       const std::string& last_pf_name,
                       const amrex::Vector<std::string>& varnames,
@@ -206,7 +206,7 @@ Matrix matrix_multiply(const Matrix& A, const Matrix& B);
 void matrix_print(const Matrix& A);
 void lapack_testing();
 
-void compute_S_matrix(const int da_iter, 
+void compute_S_matrix(const int da_iter,
                       Matrix& S,
                       const int& Nens,
                       const amrex::MultiFab& mean_H_xf,
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
index 9bb1afb789..944896fa5c 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_Utils.cpp
@@ -86,7 +86,7 @@ std::vector<std::string>
 get_plotfile_list(const int da_iter)
 {
     const std::string da_dir = "DA_Cycle_" + std::to_string(da_iter) + "/";
-    
+
     std::vector<std::string> pltfiles;
     const std::string member_prefix = da_dir + "member_";
 
@@ -160,13 +160,13 @@ read_plot_file(PlotFileData& pf,
 }
 
 MultiFab
-read_member_multifab(const int da_iter, 
+read_member_multifab(const int da_iter,
                      const int n,
                      const std::string& pf_name,
                      const Vector<std::string>& varnames)
 {
-    
-    const std::string da_dir = "DA_Cycle_" + std::to_string(da_iter) + "/";    
+
+    const std::string da_dir = "DA_Cycle_" + std::to_string(da_iter) + "/";
 
     const std::string member_prefix = da_dir + "member_";
     std::string member_dir = member_prefix + amrex::Concatenate("", n, 2);
@@ -186,7 +186,7 @@ read_member_multifab(const int da_iter,
 }
 
 MultiFab
-compute_ensemble_mean(const int da_iter, 
+compute_ensemble_mean(const int da_iter,
                       const int Nens,
                       const std::string& pf_name,
                       const Vector<std::string>& varnames)
@@ -388,7 +388,7 @@ compute_yf_prime (const int da_iter,
 }
 
 void
-compute_S_matrix(const int da_iter, 
+compute_S_matrix(const int da_iter,
                  Matrix& S,
                  const int& Nens,
                  const MultiFab& mean_H_xf,
@@ -468,7 +468,7 @@ compute_yf_prime_T_d_prime_vec (const MultiFab& yf_prime,
 }
 
 void
-compute_r_vec (const int da_iter, 
+compute_r_vec (const int da_iter,
                int Nens,
                const std::string& last_pf_name,
                const Vector<std::string>& varnames,
@@ -520,7 +520,7 @@ compute_alpha_vec (const int& Nens,
 // column-vector element multiply
 
 void
-compute_Xf_prime_times_vector (const int da_iter, 
+compute_Xf_prime_times_vector (const int da_iter,
                                const int Nens,
                                const std::string& last_pf_name,
                                const Vector<std::string>& varnames,
@@ -532,7 +532,7 @@ compute_Xf_prime_times_vector (const int da_iter,
 
     // Read first member to define result
     MultiFab xf_0 =
-        read_member_multifab(da_iter, 
+        read_member_multifab(da_iter,
                              0,
                              last_pf_name,
                              varnames);
@@ -613,7 +613,7 @@ compute_T_matrix (const Matrix& S_mat,
 
 
 void
-update_ensemble (const int da_iter, 
+update_ensemble (const int da_iter,
                  const int Nens,
                  const std::string& last_pf_name,
                  const Vector<std::string>& varnames,
@@ -630,7 +630,7 @@ update_ensemble (const int da_iter,
 }
 
 void
-ERF::SetDirsForPlotfilesAndCheckpointsForDA (const int da_iter, 
+ERF::SetDirsForPlotfilesAndCheckpointsForDA (const int da_iter,
                                              const int ens_no)
 {
     // --------------------------------------------------------

From e06bfd1ab23f2186bdc85241d28405e4ffa9f301 Mon Sep 17 00:00:00 2001
From: nataraj2 <mn556@cornell.edu>
Date: Wed, 2 Sep 2026 07:37:55 -0700
Subject: [PATCH 15/18] Correcting trailing whitespaces and tabs

---
 .../SquallLine2D/ERF_DA_EnKFSRF.cpp           | 14 +--
 .../ERF_DA_EnKFSRF_SquallLine2D.cpp           | 10 ++-
 .../DataAssimilation/SquallLine2D/main.cpp    | 19 +++--
 Source/Initialization/ERF_InitForEnsemble.cpp | 85 +++++++++++++------
 4 files changed, 86 insertions(+), 42 deletions(-)

diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
index dd7d0cf385..a9f29032a1 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
@@ -348,12 +348,16 @@ ERF::PerformDataAssimilation(int da_iter)
 
         m_plot3d_int_1 = -1;
         m_check_int = -1;
-        // This handles all the initialization including when doing a
-        // restart. But in the inputs file, we do not specify a
+        // InitData() handles all the initialization including when doing a
+        // restart. But when doing data assimilation, in the inputs file, we do not specify a
         // erf.restart. So, the ERF class does not know it is a restart.
-        // So, whenever InitData() is called, it just does the from scratch
-        // initialization. All the restart has to be explicitly handled by
-        // specifying the restart_chkfile variable for the corresponding ensemble
+        // So, whenever InitData() is called without assigning a restart_chkfile,
+        // it just does the from scratch initialization. This function is called by the
+        // temporararily created erf class - tmp_erf, from main.cpp. We simply want to fill all the data
+        // structures here by calling InitData(), and then the WriteUpdatedEnsembleToERFClassData will write
+        // all the updated ensemble data into the lev_new data structure, which will then be
+        // written into the checkpoint files by MakeEnsembleCheckpointName.
+
         InitData();
         auto& lev_new = vars_new[0];
 
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_SquallLine2D.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_SquallLine2D.cpp
index d621b013b5..4e2bac6699 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_SquallLine2D.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF_SquallLine2D.cpp
@@ -96,6 +96,12 @@ void compute_R_diag_vals(Vector<Real>& R_diag)
     // Set the size to 2
     R_diag.resize(8,0.01);
 
-    R_diag[0] = 0.01;
-    R_diag[1] = 0.01;
+    R_diag[0] = 4e-4; //density
+    R_diag[1] = 1.0; // theta
+    R_diag[2] = 1.0; // x_velocity
+    R_diag[3] = 1.0; // y_velocity
+    R_diag[4] = 0.25; // z_velocity
+    R_diag[5] = 2.5e-7; // qv
+    R_diag[6] = 1e-8; // qc
+    R_diag[7] = 1e-8; // qrain
 }
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
index 625f94447b..7df4de80fd 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
@@ -124,14 +124,17 @@ int main (int argc, char* argv[])
             // restart. But in the inputs file, for DA runs, we do not specify a
             // erf.restart. So, the ERF class does not know it is a restart.
             // So, whenever InitData() is called, it just does the from-scratch
-            // initialization. All the restart has to be explicitly handled by
-            // specifying the restart_chkfile variable for the corresponding ensemble
-
-
-            // So, for da_iter=0, it starts from scratch, and for iteration da_iter>=1
-            // there has to be a erf.restart_chkfile specified, which is the checkpoint file
-            // which contains the latest updated ensemble after data assimilation for that ensemble.
-            // The InitData reads that restart file
+            // initialization, if restart_chkfile is not speficied.
+            // So, for InitData () to do the restart pathway, the restart has
+            // to be explicitly handled by specifying the restart_chkfile variable
+            // for the corresponding ensemble.
+
+            // So, for da_iter=0, ie. the first DA cycle, the simulation starts from scratch
+            // for all ensembles, and for iteration da_iter>=1,
+            // there has to be a erf.restart_chkfile specified (for each ensemble),
+            // which is the checkpoint file which contains the latest updated ensemble after
+            // data assimilation for that ensemble. The InitData then reads that checkpoint file
+            // and continues the simulation from that.
 
             // This call will fill the restart_chkfile string
             if(da_iter > 0){
diff --git a/Source/Initialization/ERF_InitForEnsemble.cpp b/Source/Initialization/ERF_InitForEnsemble.cpp
index bb9d3a72ab..b9d9f9e85f 100644
--- a/Source/Initialization/ERF_InitForEnsemble.cpp
+++ b/Source/Initialization/ERF_InitForEnsemble.cpp
@@ -10,6 +10,18 @@
 using namespace amrex;
 namespace fs = std::filesystem;
 
+perturb_scale[8] =
+{
+    0.01,       // density [kg/m3]
+    1.0,        // theta [K]
+    1.0,        // u [m/s]
+    1.0,        // v [m/s]
+    0.5,        // w [m/s]
+    5.0e-4,     // qv [kg/kg]
+    1.0e-4,      // qc [kg/kg]
+    1.0e-4      // qrain [kg/kg]
+};
+
 /**
  * Create a cell-centered MultiFab of random perturbations for one AMR level.
  *
@@ -94,52 +106,66 @@ void NormalizeMultiFabRMS_PerComponent(MultiFab& mf_cc_pert)
 }
 
 /**
- * Apply horizontal Gaussian smoothing to cell-centered perturbations.
+ * Apply 3D isotropic/anisotropic Gaussian smoothing to cell-centered perturbations.
  *
  * @param lev Integer specifying the current level
  * @param mf_cc_pert MultiFab containing perturbations to smooth in place
  */
 void
 ERF::apply_gaussian_smoothing_to_perturbations(const int lev,
-                                               MultiFab& mf_cc_pert)
+                                                MultiFab& mf_cc_pert)
 {
     const Geometry& gm = geom[lev];
     const Real dx = gm.CellSize(0);
     const Real dy = gm.CellSize(1);
+    const Real dz = gm.CellSize(2);
 
-    const Real dmesh = std::min(dx, dy);
+    // Dynamic scale selection based on minimum grid spacing
+    const Real dmesh = std::min({dx, dy, dz});
 
-    // ---- User choice ----
+    // Correlation radius (sigma) from solverChoice
     const Real sigma = solverChoice.ens_pert_correlated_radius;
-    const int  r     = static_cast<int>(3.0 * sigma / dmesh);
+
+    // Truncate stencil at 3*sigma
+    const int rx = static_cast<int>(std::ceil(3.0 * sigma / dx));
+    const int ry = static_cast<int>(std::ceil(3.0 * sigma / dy));
+    const int rz = static_cast<int>(std::ceil(3.0 * sigma / dz));
 
     const int ncomp = mf_cc_pert.nComp();
 
-    // ---- Precompute Gaussian weights ----
-    const int wsize = 2*r + 1;
-    Vector<Real> w_host(wsize * wsize);
+    // ---- Precompute 3D Gaussian weights on Host ----
+    const int wx_size = 2 * rx + 1;
+    const int wy_size = 2 * ry + 1;
+    const int wz_size = 2 * rz + 1;
+
+    Vector<Real> w_host(wx_size * wy_size * wz_size);
 
     Real Z = zero;
-    for (int m = -r; m <= r; ++m) {
-        for (int n = -r; n <= r; ++n) {
-            Real val = std::exp(-(m*m*dx*dx + n*n*dy*dy)
-                                 /(two*sigma*sigma));
-            w_host[(m+r)*wsize + (n+r)] = val;
-            Z += val;
+    for (int m = -rx; m <= rx; ++m) {
+        for (int n = -ry; n <= ry; ++n) {
+            for (int p = -rz; p <= rz; ++p) {
+                Real r_sq = (m * m * dx * dx) + (n * n * dy * dy) + (p * p * dz * dz);
+                Real val  = std::exp(-r_sq / (two * sigma * sigma));
+
+                int idx = (m + rx) * (wy_size * wz_size) + (n + ry) * wz_size + (p + rz);
+                w_host[idx] = val;
+                Z += val;
+            }
         }
     }
 
+    // Normalize weights
     for (auto& v : w_host) {
         v /= Z;
     }
 
+    // Copy weights to Device memory
     Gpu::DeviceVector<Real> w_dev(w_host.size());
     Gpu::copy(Gpu::hostToDevice, w_host.begin(), w_host.end(), w_dev.begin());
-
     Real const* w = w_dev.data();
 
-    // ---- Create a grown copy (for stencil access) ----
-    IntVect ngrow_big(AMREX_D_DECL(r, r, 0));
+    // ---- Create a grown copy in 3D for stencil ghost access ----
+    IntVect ngrow_big(AMREX_D_DECL(rx, ry, rz));
 
     MultiFab mf_copy(mf_cc_pert.boxArray(),
                      mf_cc_pert.DistributionMap(),
@@ -150,7 +176,7 @@ ERF::apply_gaussian_smoothing_to_perturbations(const int lev,
                          IntVect(0), ngrow_big,
                          gm.periodicity());
 
-    // ---- Apply smoothing ----
+    // ---- Apply 3D Smoothing Kernel ----
     for (MFIter mfi(mf_cc_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi)
     {
         const Box& bx = mfi.tilebox();
@@ -163,16 +189,21 @@ ERF::apply_gaussian_smoothing_to_perturbations(const int lev,
         {
             Real sum = zero;
 
-            for (int m = -r; m <= r; ++m) {
-                for (int nn = -r; nn <= r; ++nn) {
-                    Real wij = w[(m+r)*wsize + (nn+r)];
-                    sum += wij * in(i+m, j+nn, k, n);
+            for (int m = -rx; m <= rx; ++m) {
+                for (int nn = -ry; nn <= ry; ++nn) {
+                    for (int p = -rz; p <= rz; ++p) {
+                        int w_idx = (m + rx) * (wy_size * wz_size) + (nn + ry) * wz_size + (p + rz);
+                        Real wij  = w[w_idx];
+                        sum += wij * in(i + m, j + nn, k + p, n);
+                    }
                 }
             }
 
-            out(i,j,k,n) = sum;
+            out(i, j, k, n) = sum;
         });
     }
+
+    // Preserve component variance scaling post-smoothing
     NormalizeMultiFabRMS_PerComponent(mf_cc_pert);
 }
 
@@ -602,9 +633,9 @@ MakeFinalMultiFabs (const MultiFab& mf_cc_fine,
             Real tmp_qrain = mf_cc_fine_arr(i,j,k,7);
             cons_pert_arr(i,j,k,Rho_comp)      = tmp_rho;
             cons_pert_arr(i,j,k,RhoTheta_comp) = tmp_rho*tmp_theta;
-            if (n_qstate_moist > 0) cons_pert_arr(i,j,k,RhoQ1_comp)    = tmp_rho*tmp_qv;
-            if (n_qstate_moist > 1) cons_pert_arr(i,j,k,RhoQ2_comp)    = tmp_rho*tmp_qc;
-            if (n_qstate_moist > 2) cons_pert_arr(i,j,k,RhoQ3_comp)    = tmp_rho*tmp_qrain;
+            if (n_qstate_moist > 0) cons_pert_arr(i,j,k,RhoQ1_comp)    = std::max(tmp_rho*tmp_qv,0.0);
+            if (n_qstate_moist > 1) cons_pert_arr(i,j,k,RhoQ2_comp)    = std::max(tmp_rho*tmp_qc,0.0);
+            if (n_qstate_moist > 2) cons_pert_arr(i,j,k,RhoQ3_comp)    = std::max(tmp_rho*tmp_qrain,0.0);
         });
     }
 
@@ -675,7 +706,7 @@ AddPertToBckgnd(MultiFab& mf_cc_fine,
         amrex::ParallelFor(bx, ncomp,
         [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
         {
-            Real ens_amp = ens_pert_amplitude*std::abs(bg(i,j,k,n));
+            Real ens_amp = ens_pert_amplitude*perturb_scale*std::abs(bg(i,j,k,n));
             bg(i,j,k,n) += ens_amp*pert(i,j,k,n);
         });
     }

From b567d2bb8fb15dd2d912668271b1e58b771787b6 Mon Sep 17 00:00:00 2001
From: Mahesh Natarajan <mahesh@dogora.dhcp.lbl.gov>
Date: Tue, 8 Sep 2026 09:12:04 -0700
Subject: [PATCH 16/18] Better version of Gaussian smoothing for perturbations

---
 .../SquallLine2D/ERF_DA_EnKFSRF.cpp           |   5 +-
 .../DataAssimilation/SquallLine2D/main.cpp    |   1 +
 .../ERF_InitCustomPertState.cpp               |   3 +
 Source/Initialization/ERF_InitForEnsemble.cpp | 136 +++++++++---------
 4 files changed, 74 insertions(+), 71 deletions(-)

diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
index a9f29032a1..394fa0ebe3 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/ERF_DA_EnKFSRF.cpp
@@ -314,11 +314,10 @@ ERF::PerformDataAssimilation(int da_iter)
 
     Print() << "add_multifabs(xf_bar, Xf_prime_alpha, xf_bar_updated) complete" << std::endl;
 
-    /*std::string pltname4 = WriteSingleLevelPlotfile("plt_xf_bar_updated", da_iter);
-    Vector<std::string> varnames4 = {"density", "theta", "x_velocity", "y_velocity", "z_velocity"};
+    std::string pltname4 = MakeFullPath("plt_xf_bar_updated", da_iter);
     WriteSingleLevelPlotfile(pltname4,
                             xf_bar_updated,
-                            varnames4,
+                            varnames,
                             geom[0],
                             0.0,   // time
                             0);    // level*/
diff --git a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
index 7df4de80fd..1f37f0c2ad 100644
--- a/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
+++ b/.Exec_dev/DataAssimilation/SquallLine2D/main.cpp
@@ -109,6 +109,7 @@ int main (int argc, char* argv[])
 
     // Data assimilation cycle loop
     for (int da_iter = 0; da_iter < num_da_cycles; ++da_iter) {
+        Print() << "Doing data assimilation cycle " << da_iter <<std::endl;
         // Ensemble run loop
         for (int ens_no = 0; ens_no < n_ens; ++ens_no)
         {
diff --git a/Source/Initialization/ERF_InitCustomPertState.cpp b/Source/Initialization/ERF_InitCustomPertState.cpp
index 25bd5ebc73..3601e34a27 100644
--- a/Source/Initialization/ERF_InitCustomPertState.cpp
+++ b/Source/Initialization/ERF_InitCustomPertState.cpp
@@ -156,9 +156,12 @@ ERF::init_custom (int lev)
     //    add the perturbations to the background state and then populate the "pert variables
 
     if(solverChoice.is_init_for_ensemble) {
+        Vector<std::string> varnames = {"density","theta", "x_velocity","y_velocity","z_velocity", "qv", "qc", "qrain"};
         MultiFab mf_cc_pert;
         create_random_perturbations(lev, mf_cc_pert);
+        WriteSingleLevelPlotfile("plt_pert_correlated_gaussian", mf_cc_pert, varnames, geom[lev], zero, 0);
         apply_gaussian_smoothing_to_perturbations(lev, mf_cc_pert);
+        //WriteSingleLevelPlotfile("plt_pert_correlated_gaussian", mf_cc_pert, varnames, geom[lev], zero, 0);
         create_background_state_for_ensemble(lev, mf_cc_pert, lev_new[Vars::cons], lev_new[Vars::xvel], lev_new[Vars::yvel], lev_new[Vars::zvel]);
     }
 }
diff --git a/Source/Initialization/ERF_InitForEnsemble.cpp b/Source/Initialization/ERF_InitForEnsemble.cpp
index b9d9f9e85f..6fe8cdca9e 100644
--- a/Source/Initialization/ERF_InitForEnsemble.cpp
+++ b/Source/Initialization/ERF_InitForEnsemble.cpp
@@ -10,7 +10,8 @@
 using namespace amrex;
 namespace fs = std::filesystem;
 
-perturb_scale[8] =
+Vector<Real>
+perturb_scale =
 {
     0.01,       // density [kg/m3]
     1.0,        // theta [K]
@@ -111,99 +112,98 @@ void NormalizeMultiFabRMS_PerComponent(MultiFab& mf_cc_pert)
  * @param lev Integer specifying the current level
  * @param mf_cc_pert MultiFab containing perturbations to smooth in place
  */
+
 void
-ERF::apply_gaussian_smoothing_to_perturbations(const int lev,
-                                                MultiFab& mf_cc_pert)
+ERF::apply_gaussian_smoothing_to_perturbations(const int lev, MultiFab& mf_cc_pert)
 {
     const Geometry& gm = geom[lev];
     const Real dx = gm.CellSize(0);
     const Real dy = gm.CellSize(1);
     const Real dz = gm.CellSize(2);
 
-    // Dynamic scale selection based on minimum grid spacing
-    const Real dmesh = std::min({dx, dy, dz});
-
-    // Correlation radius (sigma) from solverChoice
     const Real sigma = solverChoice.ens_pert_correlated_radius;
 
-    // Truncate stencil at 3*sigma
     const int rx = static_cast<int>(std::ceil(3.0 * sigma / dx));
     const int ry = static_cast<int>(std::ceil(3.0 * sigma / dy));
     const int rz = static_cast<int>(std::ceil(3.0 * sigma / dz));
 
     const int ncomp = mf_cc_pert.nComp();
 
-    // ---- Precompute 3D Gaussian weights on Host ----
-    const int wx_size = 2 * rx + 1;
-    const int wy_size = 2 * ry + 1;
-    const int wz_size = 2 * rz + 1;
-
-    Vector<Real> w_host(wx_size * wy_size * wz_size);
-
-    Real Z = zero;
-    for (int m = -rx; m <= rx; ++m) {
-        for (int n = -ry; n <= ry; ++n) {
-            for (int p = -rz; p <= rz; ++p) {
-                Real r_sq = (m * m * dx * dx) + (n * n * dy * dy) + (p * p * dz * dz);
-                Real val  = std::exp(-r_sq / (two * sigma * sigma));
-
-                int idx = (m + rx) * (wy_size * wz_size) + (n + ry) * wz_size + (p + rz);
-                w_host[idx] = val;
-                Z += val;
-            }
+    // 1. Build 1D Gaussian weight vectors on host
+    auto compute_1d_weights = [](int r, Real d, Real sig) {
+        Vector<Real> w(2*r + 1);
+        Real sum = 0.0;
+        for (int m = -r; m <= r; ++m) {
+            Real val = std::exp(-(m * m * d * d) / (2.0 * sig * sig));
+            w[m + r] = val;
+            sum += val;
         }
+        for (auto& v : w) v /= sum;
+        return w;
+    };
+
+    Vector<Real> wx_h = compute_1d_weights(rx, dx, sigma);
+    Vector<Real> wy_h = compute_1d_weights(ry, dy, sigma);
+    Vector<Real> wz_h = compute_1d_weights(rz, dz, sigma);
+
+    Gpu::DeviceVector<Real> wx_d(wx_h.size()), wy_d(wy_h.size()), wz_d(wz_h.size());
+    Gpu::copy(Gpu::hostToDevice, wx_h.begin(), wx_h.end(), wx_d.begin());
+    Gpu::copy(Gpu::hostToDevice, wy_h.begin(), wy_h.end(), wy_d.begin());
+    Gpu::copy(Gpu::hostToDevice, wz_h.begin(), wz_h.end(), wz_d.begin());
+
+    Real const* wx = wx_d.data();
+    Real const* wy = wy_d.data();
+    Real const* wz = wz_d.data();
+
+    // Temp MultiFab for separable intermediate states
+    IntVect ngrow(AMREX_D_DECL(rx, ry, rz));
+    MultiFab mf_tmp(mf_cc_pert.boxArray(), mf_cc_pert.DistributionMap(), ncomp, ngrow);
+
+    // --- Pass 1: X Direction ---
+    mf_tmp.ParallelCopy(mf_cc_pert, 0, 0, ncomp, IntVect(0), ngrow, gm.periodicity());
+    for (MFIter mfi(mf_cc_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+        const Box& bx = mfi.tilebox();
+        auto const& in  = mf_tmp.const_array(mfi);
+        auto const& out = mf_cc_pert.array(mfi);
+        ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept {
+            Real sum = 0.0;
+            for (int m = -rx; m <= rx; ++m) {
+                sum += wx[m + rx] * in(i + m, j, k, n);
+            }
+            out(i, j, k, n) = sum;
+        });
     }
 
-    // Normalize weights
-    for (auto& v : w_host) {
-        v /= Z;
+    // --- Pass 2: Y Direction ---
+    mf_tmp.ParallelCopy(mf_cc_pert, 0, 0, ncomp, IntVect(0), ngrow, gm.periodicity());
+    for (MFIter mfi(mf_cc_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+        const Box& bx = mfi.tilebox();
+        auto const& in  = mf_tmp.const_array(mfi);
+        auto const& out = mf_cc_pert.array(mfi);
+        ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept {
+            Real sum = 0.0;
+            for (int nn = -ry; nn <= ry; ++nn) {
+                sum += wy[nn + ry] * in(i, j + nn, k, n);
+            }
+            out(i, j, k, n) = sum;
+        });
     }
 
-    // Copy weights to Device memory
-    Gpu::DeviceVector<Real> w_dev(w_host.size());
-    Gpu::copy(Gpu::hostToDevice, w_host.begin(), w_host.end(), w_dev.begin());
-    Real const* w = w_dev.data();
-
-    // ---- Create a grown copy in 3D for stencil ghost access ----
-    IntVect ngrow_big(AMREX_D_DECL(rx, ry, rz));
-
-    MultiFab mf_copy(mf_cc_pert.boxArray(),
-                     mf_cc_pert.DistributionMap(),
-                     ncomp, ngrow_big);
-
-    mf_copy.ParallelCopy(mf_cc_pert,
-                         0, 0, ncomp,
-                         IntVect(0), ngrow_big,
-                         gm.periodicity());
-
-    // ---- Apply 3D Smoothing Kernel ----
-    for (MFIter mfi(mf_cc_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi)
-    {
+    // --- Pass 3: Z Direction ---
+    mf_tmp.ParallelCopy(mf_cc_pert, 0, 0, ncomp, IntVect(0), ngrow, gm.periodicity());
+    for (MFIter mfi(mf_cc_pert, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
         const Box& bx = mfi.tilebox();
-
-        auto const& in  = mf_copy.const_array(mfi);
+        auto const& in  = mf_tmp.const_array(mfi);
         auto const& out = mf_cc_pert.array(mfi);
-
-        ParallelFor(bx, ncomp,
-        [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
-        {
-            Real sum = zero;
-
-            for (int m = -rx; m <= rx; ++m) {
-                for (int nn = -ry; nn <= ry; ++nn) {
-                    for (int p = -rz; p <= rz; ++p) {
-                        int w_idx = (m + rx) * (wy_size * wz_size) + (nn + ry) * wz_size + (p + rz);
-                        Real wij  = w[w_idx];
-                        sum += wij * in(i + m, j + nn, k + p, n);
-                    }
-                }
+        ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept {
+            Real sum = 0.0;
+            for (int p = -rz; p <= rz; ++p) {
+                sum += wz[p + rz] * in(i, j, k + p, n);
             }
-
             out(i, j, k, n) = sum;
         });
     }
 
-    // Preserve component variance scaling post-smoothing
     NormalizeMultiFabRMS_PerComponent(mf_cc_pert);
 }
 
@@ -706,7 +706,7 @@ AddPertToBckgnd(MultiFab& mf_cc_fine,
         amrex::ParallelFor(bx, ncomp,
         [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
         {
-            Real ens_amp = ens_pert_amplitude*perturb_scale*std::abs(bg(i,j,k,n));
+            Real ens_amp = ens_pert_amplitude*perturb_scale[n]*std::abs(bg(i,j,k,n));
             bg(i,j,k,n) += ens_amp*pert(i,j,k,n);
         });
     }

From b246b8ad213c54c9ba115436dbcf4ebca171698e Mon Sep 17 00:00:00 2001
From: nataraj2 <mn556@cornell.edu>
Date: Tue, 8 Sep 2026 12:40:10 -0700
Subject: [PATCH 17/18] Making it work on gpus

---
 Source/Initialization/ERF_InitForEnsemble.cpp | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/Source/Initialization/ERF_InitForEnsemble.cpp b/Source/Initialization/ERF_InitForEnsemble.cpp
index 6fe8cdca9e..145a43b698 100644
--- a/Source/Initialization/ERF_InitForEnsemble.cpp
+++ b/Source/Initialization/ERF_InitForEnsemble.cpp
@@ -696,6 +696,12 @@ AddPertToBckgnd(MultiFab& mf_cc_fine,
     // Optional safety check (recommended)
     AMREX_ALWAYS_ASSERT(mf_cc_pert.nComp() == ncomp);
 
+    // Create a GPU-accessible array and copy the host vector into it
+    GpuArray<Real, 8> scale_gpu; // Use the maximum possible ncomp, e.g., 8
+    for (int n = 0; n < ncomp; ++n) {
+        scale_gpu[n] = perturb_scale[n];
+    }
+
     for (MFIter mfi(mf_cc_fine, TilingIfNotGPU()); mfi.isValid(); ++mfi)
     {
         const Box& bx = mfi.tilebox();
@@ -706,7 +712,7 @@ AddPertToBckgnd(MultiFab& mf_cc_fine,
         amrex::ParallelFor(bx, ncomp,
         [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
         {
-            Real ens_amp = ens_pert_amplitude*perturb_scale[n]*std::abs(bg(i,j,k,n));
+            Real ens_amp = ens_pert_amplitude*scale_gpu[n]*std::abs(bg(i,j,k,n));
             bg(i,j,k,n) += ens_amp*pert(i,j,k,n);
         });
     }

From ada30fdff89427be3e31cd9ca726d20731ae51c5 Mon Sep 17 00:00:00 2001
From: nataraj2 <mn556@cornell.edu>
Date: Tue, 8 Sep 2026 19:35:58 -0700
Subject: [PATCH 18/18] Correcting compilation errors

---
 Source/Initialization/ERF_InitForEnsemble.cpp | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/Source/Initialization/ERF_InitForEnsemble.cpp b/Source/Initialization/ERF_InitForEnsemble.cpp
index 145a43b698..6ea72012c8 100644
--- a/Source/Initialization/ERF_InitForEnsemble.cpp
+++ b/Source/Initialization/ERF_InitForEnsemble.cpp
@@ -633,9 +633,9 @@ MakeFinalMultiFabs (const MultiFab& mf_cc_fine,
             Real tmp_qrain = mf_cc_fine_arr(i,j,k,7);
             cons_pert_arr(i,j,k,Rho_comp)      = tmp_rho;
             cons_pert_arr(i,j,k,RhoTheta_comp) = tmp_rho*tmp_theta;
-            if (n_qstate_moist > 0) cons_pert_arr(i,j,k,RhoQ1_comp)    = std::max(tmp_rho*tmp_qv,0.0);
-            if (n_qstate_moist > 1) cons_pert_arr(i,j,k,RhoQ2_comp)    = std::max(tmp_rho*tmp_qc,0.0);
-            if (n_qstate_moist > 2) cons_pert_arr(i,j,k,RhoQ3_comp)    = std::max(tmp_rho*tmp_qrain,0.0);
+            if (n_qstate_moist > 0) cons_pert_arr(i,j,k,RhoQ1_comp)  = std::max(tmp_rho*tmp_qv,0.0_rt);
+            if (n_qstate_moist > 1) cons_pert_arr(i,j,k,RhoQ2_comp)  = std::max(tmp_rho*tmp_qc,0.0_rt);
+            if (n_qstate_moist > 2) cons_pert_arr(i,j,k,RhoQ3_comp)  = std::max(tmp_rho*tmp_qrain,0.0_rt);
         });
     }