diff --git a/Source/Derive.H b/Source/Derive.H
index f40bdc432..ef974f73b 100644
--- a/Source/Derive.H
+++ b/Source/Derive.H
@@ -45,6 +45,17 @@ void erf_dertemp (
   const int* bcrec,
   const int level);
 
+void erf_dermoisttemp (
+  const amrex::Box& bx,
+  amrex::FArrayBox& derfab,
+  int dcomp,
+  int ncomp,
+  const amrex::FArrayBox& datfab,
+  const amrex::Geometry& geomdata,
+  amrex::Real time,
+  const int* bcrec,
+  const int level);
+
 void erf_dertheta (
   const amrex::Box& bx,
   amrex::FArrayBox& derfab,
diff --git a/Source/Derive.cpp b/Source/Derive.cpp
index 6df9c5b33..5d3ed57d9 100644
--- a/Source/Derive.cpp
+++ b/Source/Derive.cpp
@@ -109,6 +109,29 @@ erf_dertemp (const Box& bx,
         tfab(i,j,k) = getTgivenRandRTh(rho,rhotheta);
     });
 }
+void
+erf_dermoisttemp (const Box& bx,
+             FArrayBox& derfab,
+             int /*dcomp*/,
+             int /*ncomp*/,
+             const FArrayBox& datfab,
+             const Geometry& /*geomdata*/,
+             Real /*time*/,
+             const int* /*bcrec*/,
+             const int /*level*/)
+{
+    auto const dat = datfab.array();
+    auto tfab      = derfab.array();
+
+    ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+    {
+        const Real rho = dat(i, j, k, Rho_comp);
+        const Real rhotheta = dat(i, j, k, RhoTheta_comp);
+        AMREX_ALWAYS_ASSERT(rhotheta > 0.);
+        const Real qv = dat(i, j, k, RhoQ1_comp) / rho;
+        tfab(i,j,k) = getTgivenRandRTh(rho,rhotheta,qv);
+    });
+}
 
 /**
  * Function to define the potential temperature by calling an EOS routine
diff --git a/Source/IO/ERF_WriteBndryPlanes.cpp b/Source/IO/ERF_WriteBndryPlanes.cpp
index eda3353c9..002ce5f2d 100644
--- a/Source/IO/ERF_WriteBndryPlanes.cpp
+++ b/Source/IO/ERF_WriteBndryPlanes.cpp
@@ -143,6 +143,9 @@ void WriteBndryPlanes::write_planes (const int t_step, const Real time,
     Box target_box_shifted(IntVect(0,0,0),new_hi);
     BoxArray ba_shifted(target_box_shifted);
 
+    int n_moist_var = NMOIST_max - (S.nComp() - NVAR_max);
+    bool ismoist = (n_moist_var >= 1);
+
     for (int i = 0; i < m_var_names.size(); i++)
     {
         std::string var_name = m_var_names[i];
@@ -169,7 +172,11 @@ void WriteBndryPlanes::write_planes (const int t_step, const Real time,
             for (MFIter mfi(Temp, TilingIfNotGPU()); mfi.isValid(); ++mfi)
             {
                 const Box& bx = mfi.tilebox();
-                derived::erf_dertemp(bx, Temp[mfi], 0, 1, S[mfi], m_geom[bndry_lev], time, nullptr, bndry_lev);
+                if (ismoist) {
+                    derived::erf_dermoisttemp(bx, Temp[mfi], 0, 1, S[mfi], m_geom[bndry_lev], time, nullptr, bndry_lev);
+                } else {
+                    derived::erf_dertemp(bx, Temp[mfi], 0, 1, S[mfi], m_geom[bndry_lev], time, nullptr, bndry_lev);
+                }
             }
             bndry.copyFrom(Temp, nghost, 0, 0, ncomp, m_geom[bndry_lev].periodicity());
         } else if (var_name == "scalar") {
diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp
index 09be7891a..dfca13bda 100644
--- a/Source/IO/Plotfile.cpp
+++ b/Source/IO/Plotfile.cpp
@@ -300,9 +300,15 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
             }
         };
 
+        bool ismoist = (solverChoice.moisture_type != MoistureType::None);
+
         // Note: All derived variables must be computed in order of "derived_names" defined in ERF.H
         calculate_derived("soundspeed",  vars_new[lev][Vars::cons], derived::erf_dersoundspeed);
-        calculate_derived("temp",        vars_new[lev][Vars::cons], derived::erf_dertemp);
+        if (ismoist) {
+            calculate_derived("temp",        vars_new[lev][Vars::cons], derived::erf_dermoisttemp);
+        } else {
+            calculate_derived("temp",        vars_new[lev][Vars::cons], derived::erf_dertemp);
+        }
         calculate_derived("theta",       vars_new[lev][Vars::cons], derived::erf_dertheta);
         calculate_derived("KE",          vars_new[lev][Vars::cons], derived::erf_derKE);
         calculate_derived("QKE",         vars_new[lev][Vars::cons], derived::erf_derQKE);