From 513861dc73feb91698a62eb0e40c7167e9a12cd8 Mon Sep 17 00:00:00 2001
From: Mahesh Natarajan <mn556@cornell.edu>
Date: Sat, 23 Mar 2024 16:51:48 -0700
Subject: [PATCH] Make rain accumulation work with restart (#1526)

Co-authored-by: Mahesh Natarajan <mahesh@dogora.dhcp.lbl.gov>
---
 Source/IO/Checkpoint.cpp | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)

diff --git a/Source/IO/Checkpoint.cpp b/Source/IO/Checkpoint.cpp
index 0830d6ade..020972485 100644
--- a/Source/IO/Checkpoint.cpp
+++ b/Source/IO/Checkpoint.cpp
@@ -142,6 +142,16 @@ ERF::WriteCheckpointFile () const
             VisMF::Write(z_height, MultiFabFileFullPrefix(lev, checkpointname, "Level_", "Z_Phys_nd"));
         }
 
+         // We must read and write qmoist with ghost cells because we don't directly impose BCs on these vars
+         // Write the precipitation accumulation component only
+        if (solverChoice.moisture_type == MoistureType::Kessler) {
+            ng = qmoist[lev][0]->nGrowVect();
+            int nvar = 1;
+             MultiFab moist_vars(grids[lev],dmap[lev],nvar,ng);
+             MultiFab::Copy(moist_vars,*(qmoist[lev][0]),0,0,nvar,ng);
+             VisMF::Write(moist_vars, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MoistVars"));
+        }
+
         if (solverChoice.lsm_type != LandSurfaceType::None) {
             for (int mvar(0); mvar<lsm_data[lev].size(); ++mvar) {
                 BoxArray ba = lsm_data[lev][mvar]->boxArray();
@@ -356,6 +366,15 @@ ERF::ReadCheckpointFile ()
            update_terrain_arrays(lev, t_new[lev]);
         }
 
+        // Read in the precipitation accumulation component
+        if (solverChoice.moisture_type == MoistureType::Kessler) {
+            ng = qmoist[lev][0]->nGrowVect();
+            int nvar = 1;
+            MultiFab moist_vars(grids[lev],dmap[lev],nvar,ng);
+            VisMF::Read(moist_vars, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MoistVars"));
+            MultiFab::Copy(*(qmoist[lev][0]),moist_vars,0,0,nvar,ng);
+        }
+
         if (solverChoice.lsm_type != LandSurfaceType::None) {
             for (int mvar(0); mvar<lsm_data[lev].size(); ++mvar) {
                 BoxArray ba = lsm_data[lev][mvar]->boxArray();