diff --git a/components/eamxx/src/diagnostics/potential_temperature.cpp b/components/eamxx/src/diagnostics/potential_temperature.cpp
index ff197eaedea..988bb555224 100644
--- a/components/eamxx/src/diagnostics/potential_temperature.cpp
+++ b/components/eamxx/src/diagnostics/potential_temperature.cpp
@@ -7,7 +7,26 @@ namespace scream
 PotentialTemperatureDiagnostic::PotentialTemperatureDiagnostic (const ekat::Comm& comm, const ekat::ParameterList& params)
   : AtmosphereDiagnostic(comm,params)
 {
-  // Nothing to do here
+  EKAT_REQUIRE_MSG(params.isParameter("Temperature Kind"),
+      "Error! PotentialTemperatureDiagnostic requires 'Temperature Kind' in its input parameters.\n");
+  
+  auto pt_type = params.get<std::string>("Temperature Kind");
+
+  if (pt_type=="Tot"){
+    m_ptype = "PotentialTemperature";
+  } else if (pt_type=="Liq") {
+    m_ptype = "LiqPotentialTemperature";
+  } else {
+    EKAT_ERROR_MSG (
+        "Error! Invalid choice for 'TemperatureKind' in PotentialTemperatureDiagnostic.\n"
+        "  - input value: " + pt_type + "\n"
+        "  - valid values: Tot, Liq\n");
+  }
+}
+
+std::string PotentialTemperatureDiagnostic::name() const
+{
+  return m_ptype;
 }
 
 // =========================================================================================
@@ -30,6 +49,9 @@ void PotentialTemperatureDiagnostic::set_grids(const std::shared_ptr<const Grids
   // The fields required for this diagnostic to be computed
   add_field<Required>("T_mid",          scalar3d_layout_mid, K,  grid_name, ps);
   add_field<Required>("p_mid",          scalar3d_layout_mid, Pa, grid_name, ps);
+  // Only needed for LiqPotentialTemperature, but put it here for ease
+  // TODO: only request it if it is needed
+  add_field<Required>("qc",             scalar3d_layout_mid, Q,  grid_name, ps);
 
   // Construct and allocate the diagnostic field
   FieldIdentifier fid (name(), scalar3d_layout_mid, K, grid_name);
@@ -41,17 +63,27 @@ void PotentialTemperatureDiagnostic::set_grids(const std::shared_ptr<const Grids
 // =========================================================================================
 void PotentialTemperatureDiagnostic::compute_diagnostic_impl()
 {
+  bool is_liq = (m_ptype=="LiqPotentialTemperature");
 
   const auto npacks  = ekat::npack<Pack>(m_num_levs);
   auto theta = m_diagnostic_output.get_view<Pack**>();
   auto T_mid = get_field_in("T_mid").get_view<const Pack**>();
   auto p_mid = get_field_in("p_mid").get_view<const Pack**>();
+  auto q_mid = get_field_in("qc").get_view<const Pack**>();
+
   Kokkos::parallel_for("PotentialTemperatureDiagnostic",
                        Kokkos::RangePolicy<>(0,m_num_cols*npacks),
                        KOKKOS_LAMBDA (const int& idx) {
       const int icol  = idx / npacks;
       const int jpack = idx % npacks;
-      theta(icol,jpack) = PF::calculate_theta_from_T(T_mid(icol,jpack),p_mid(icol,jpack));
+      auto temp = PF::calculate_theta_from_T(T_mid(icol,jpack),p_mid(icol,jpack));
+      if (is_liq) {
+        // Liquid potential temperature (consistent with how it is calculated in SHOC)
+        theta(icol,jpack) = PF::calculate_thetal_from_theta(temp,T_mid(icol,jpack),q_mid(icol,jpack));
+      } else {
+        // The total potential temperature
+        theta(icol,jpack) = temp;
+      }
   });
   Kokkos::fence();
 }
diff --git a/components/eamxx/src/diagnostics/potential_temperature.hpp b/components/eamxx/src/diagnostics/potential_temperature.hpp
index 933a6935005..37fd1a30806 100644
--- a/components/eamxx/src/diagnostics/potential_temperature.hpp
+++ b/components/eamxx/src/diagnostics/potential_temperature.hpp
@@ -16,12 +16,13 @@ class PotentialTemperatureDiagnostic : public AtmosphereDiagnostic
 public:
   using Pack          = ekat::Pack<Real,SCREAM_PACK_SIZE>;
   using PF            = scream::PhysicsFunctions<DefaultDevice>;
+  using C             = physics::Constants<Real>;
 
   // Constructors
   PotentialTemperatureDiagnostic (const ekat::Comm& comm, const ekat::ParameterList& params);
 
   // The name of the diagnostic
-  std::string name () const { return "PotentialTemperature"; }
+  std::string name () const;
 
   // Set the grid
   void set_grids (const std::shared_ptr<const GridsManager> grids_manager);
@@ -37,6 +38,9 @@ class PotentialTemperatureDiagnostic : public AtmosphereDiagnostic
   Int m_num_cols;
   Int m_num_levs;
 
+  // What type of potential temperature to compute
+  std::string m_ptype;
+
 }; // class PotentialTemperatureDiagnostic
 
 } //namespace scream
diff --git a/components/eamxx/src/diagnostics/tests/potential_temperature_test.cpp b/components/eamxx/src/diagnostics/tests/potential_temperature_test.cpp
index b11d8aa72ef..1f584ec7e71 100644
--- a/components/eamxx/src/diagnostics/tests/potential_temperature_test.cpp
+++ b/components/eamxx/src/diagnostics/tests/potential_temperature_test.cpp
@@ -40,7 +40,7 @@ create_gm (const ekat::Comm& comm, const int ncols, const int nlevs) {
 
 //-----------------------------------------------------------------------------------------------//
 template<typename DeviceT>
-void run(std::mt19937_64& engine)
+void run(std::mt19937_64& engine, int int_ptype)
 {
   using PF         = scream::PhysicsFunctions<DeviceT>;
   using PC         = scream::physics::Constants<Real>;
@@ -67,7 +67,8 @@ void run(std::mt19937_64& engine)
 
   // Input (randomized) views
   view_1d temperature("temperature",num_mid_packs),
-          pressure("pressure",num_mid_packs);
+          pressure("pressure",num_mid_packs),
+          condensate("condensate",num_mid_packs);
 
   auto dview_as_real = [&] (const view_1d& v) -> rview_1d {
     return rview_1d(reinterpret_cast<Real*>(v.data()),v.size()*packsize);
@@ -85,10 +86,11 @@ void run(std::mt19937_64& engine)
   ekat::ParameterList params;
   register_diagnostics();
   auto& diag_factory = AtmosphereDiagnosticFactory::instance();
+  std::string ptype = int_ptype == 0 ? "Tot" : "Liq";
+  params.set("Temperature Kind", ptype);
   auto diag = diag_factory.create("PotentialTemperature",comm,params);
   diag->set_grids(gm);
 
-
   // Set the required fields for the diagnostic.
   std::map<std::string,Field> input_fields;
   for (const auto& req : diag->get_required_field_requests()) {
@@ -114,13 +116,18 @@ void run(std::mt19937_64& engine)
     const auto& T_mid_v = T_mid_f.get_view<Pack**>();
     const auto& p_mid_f = input_fields["p_mid"];
     const auto& p_mid_v = p_mid_f.get_view<Pack**>();
+    const auto& q_mid_f = input_fields["qc"];
+    const auto& q_mid_v = q_mid_f.get_view<Pack**>();
     for (int icol=0;icol<ncols;icol++) {
       const auto& T_sub = ekat::subview(T_mid_v,icol);
       const auto& p_sub = ekat::subview(p_mid_v,icol);
+      const auto& q_sub = ekat::subview(q_mid_v,icol);
       ekat::genRandArray(dview_as_real(temperature), engine, pdf_temp);
       ekat::genRandArray(dview_as_real(pressure),    engine, pdf_pres);
+      ekat::genRandArray(dview_as_real(condensate),  engine, pdf_pres);
       Kokkos::deep_copy(T_sub,temperature);
       Kokkos::deep_copy(p_sub,pressure);
+      Kokkos::deep_copy(q_sub,condensate);
     }
 
     // Run diagnostic and compare with manual calculation
@@ -133,7 +140,10 @@ void run(std::mt19937_64& engine)
     Kokkos::parallel_for("", policy, KOKKOS_LAMBDA(const MemberType& team) {
       const int icol = team.league_rank();
       Kokkos::parallel_for(Kokkos::TeamVectorRange(team,num_mid_packs), [&] (const Int& jpack) {
-        theta_v(icol,jpack) = PF::calculate_theta_from_T(T_mid_v(icol,jpack),p_mid_v(icol,jpack));
+        auto theta = PF::calculate_theta_from_T(T_mid_v(icol,jpack),p_mid_v(icol,jpack));
+        if (int_ptype==1) {
+          theta_v(icol,jpack) = theta - (theta/T_mid_v(icol,jpack)) * (PC::LatVap/PC::Cpair) * q_mid_v(icol,jpack);
+        } else { theta_v(icol,jpack) = theta; }
       });
       team.team_barrier();
     });
@@ -159,7 +169,9 @@ TEST_CASE("potential_temp_test", "potential_temp_test]"){
 
   printf(" -> Testing Pack<Real,%d> scalar type...",SCREAM_PACK_SIZE);
   for (int irun=0; irun<num_runs; ++irun) {
-    run<Device>(engine);
+    for (const int int_ptype : {0,1}) {
+    run<Device>(engine, int_ptype);
+  }
   }
   printf("ok!\n");
 
diff --git a/components/eamxx/src/physics/shoc/eamxx_shoc_process_interface.hpp b/components/eamxx/src/physics/shoc/eamxx_shoc_process_interface.hpp
index 099794adafc..d3f5b05f88d 100644
--- a/components/eamxx/src/physics/shoc/eamxx_shoc_process_interface.hpp
+++ b/components/eamxx/src/physics/shoc/eamxx_shoc_process_interface.hpp
@@ -75,7 +75,6 @@ class SHOCMacrophysics : public scream::AtmosphereProcess
       const int i = team.league_rank();
 
       const Real zvir = C::ZVIR;
-      const Real latvap = C::LatVap;
       const Real cpair = C::Cpair;
       const Real ggr = C::gravit;
       const Real inv_ggr = 1/ggr;
@@ -107,8 +106,9 @@ class SHOCMacrophysics : public scream::AtmosphereProcess
         qw(i,k) = qv(i,k) + qc(i,k);
 
         // Temperature
+        // NOTE: theta_v (thv) is intentionally different from one in HOMME
         const auto theta_zt = PF::calculate_theta_from_T(T_mid(i,k),p_mid(i,k));
-        thlm(i,k) = theta_zt-(theta_zt/T_mid(i,k))*(latvap/cpair)*qc(i,k);
+        thlm(i,k) = PF::calculate_thetal_from_theta(theta_zt,T_mid(i,k),qc(i,k));
         thv(i,k)  = theta_zt*(1 + zvir*qv(i,k) - qc(i,k));
 
         // Vertical layer thickness
diff --git a/components/eamxx/src/share/io/scorpio_output.cpp b/components/eamxx/src/share/io/scorpio_output.cpp
index b1ce7d3f326..f35c7fb0c56 100644
--- a/components/eamxx/src/share/io/scorpio_output.cpp
+++ b/components/eamxx/src/share/io/scorpio_output.cpp
@@ -1369,6 +1369,14 @@ AtmosphereOutput::create_diagnostic (const std::string& diag_field_name) {
     diag_name = "VaporFlux";
     // split will return the list [X, ''], with X being whatever is before 'VapFlux'
     params.set<std::string>("Wind Component",ekat::split(diag_field_name,"VapFlux").front());
+  } else if (diag_field_name=="PotentialTemperature" or
+             diag_field_name=="LiqPotentialTemperature") {
+    diag_name = "PotentialTemperature";
+    if (diag_field_name == "LiqPotentialTemperature") {
+      params.set<std::string>("Temperature Kind", "Liq");
+    } else {
+      params.set<std::string>("Temperature Kind", "Tot");
+    }
   } else {
     diag_name = diag_field_name;
   }
diff --git a/components/eamxx/src/share/util/scream_common_physics_functions.hpp b/components/eamxx/src/share/util/scream_common_physics_functions.hpp
index 92050bccee8..87ae6604d99 100644
--- a/components/eamxx/src/share/util/scream_common_physics_functions.hpp
+++ b/components/eamxx/src/share/util/scream_common_physics_functions.hpp
@@ -94,6 +94,19 @@ struct PhysicsFunctions
   KOKKOS_INLINE_FUNCTION
   static ScalarT calculate_theta_from_T(const ScalarT& temperature, const ScalarT& pressure);
 
+  //-----------------------------------------------------------------------------------------------//
+  // Converts potential temperature to liquid potental temperature:
+  //   theta_l = theta - (theta / temperature) * (LatVap/Cpair) * qc,
+  // where
+  //   theta  is the potential temperature, [K]
+  //   temperature  is the temperature, [K]
+  //   qc is the cloud liquid mixing ratio, [kg/kg]
+  //   and the others are constants 
+  //-----------------------------------------------------------------------------------------------//
+  template<typename ScalarT>
+  KOKKOS_INLINE_FUNCTION
+  static ScalarT calculate_thetal_from_theta(const ScalarT& theta, const ScalarT& temperature, const ScalarT& qc);
+
   //-----------------------------------------------------------------------------------------------//
   // Converts potential temperature to temperature using Exners function:
   //   temperature = theta*exner(pressure),
@@ -121,7 +134,7 @@ struct PhysicsFunctions
 
   //-----------------------------------------------------------------------------------------------//
   // Compute virtual temperature
-  //   T_virtual = temperature * (qv+ep_2)/(qv+1)
+  //   T_virtual = temperature * (1+(-1+1/ep_2)*qv)
   // where
   //   ep_2        is ratio of molecular mass of water to the molecular mass of dry air
   //   temperature is the atmospheric temperature.  Units in [K].
@@ -392,6 +405,13 @@ struct PhysicsFunctions
                                       const InputProviderP& pressure,
                                       const view_1d<ScalarT>& theta);
 
+  template<typename ScalarT, typename InputProviderTheta, typename InputProviderT, typename InputProviderQ>
+  KOKKOS_INLINE_FUNCTION
+  static void calculate_thetal_from_theta(const MemberType& team,
+                                      const InputProviderTheta& theta,
+                                      const InputProviderT& temperature,
+                                      const InputProviderQ& qc,
+                                      const view_1d<ScalarT>& thetal);  
   template<typename ScalarT, typename InputProviderT, typename InputProviderP>
   KOKKOS_INLINE_FUNCTION
   static void calculate_T_from_theta (const MemberType& team,
diff --git a/components/eamxx/src/share/util/scream_common_physics_functions_impl.hpp b/components/eamxx/src/share/util/scream_common_physics_functions_impl.hpp
index 4a6a654056c..f41e389bf67 100644
--- a/components/eamxx/src/share/util/scream_common_physics_functions_impl.hpp
+++ b/components/eamxx/src/share/util/scream_common_physics_functions_impl.hpp
@@ -117,6 +117,16 @@ ScalarT PhysicsFunctions<DeviceT>::calculate_theta_from_T(const ScalarT& tempera
   return temperature/exner_function(pressure);
 }
 
+template<typename DeviceT>
+template<typename ScalarT>
+KOKKOS_INLINE_FUNCTION
+ScalarT PhysicsFunctions<DeviceT>::calculate_thetal_from_theta(const ScalarT& theta, const ScalarT& temperature, const ScalarT& qc)
+{
+  using C = scream::physics::Constants<Real>;
+
+  return theta - (theta / temperature) * (C::LatVap/C::Cpair) * qc;
+}
+
 template<typename DeviceT>
 template<typename ScalarT, typename InputProviderT, typename InputProviderP>
 KOKKOS_INLINE_FUNCTION
@@ -131,6 +141,21 @@ void PhysicsFunctions<DeviceT>::calculate_theta_from_T(const MemberType& team,
   });
 }
 
+template<typename DeviceT>
+template<typename ScalarT, typename InputProviderTheta, typename InputProviderT, typename InputProviderQ>
+KOKKOS_INLINE_FUNCTION
+void PhysicsFunctions<DeviceT>::calculate_thetal_from_theta(const MemberType& team,
+                                                            const InputProviderTheta& theta,
+                                                            const InputProviderT& temperature,
+                                                            const InputProviderQ& qc,
+                                                            const view_1d<ScalarT>& thetal)
+{
+  Kokkos::parallel_for(Kokkos::TeamVectorRange(team,thetal.extent(0)),
+                       [&] (const int k) {
+    thetal(k) = calculate_thetal_from_theta(theta(k),temperature(k),qc(k));
+  });
+}
+
 template<typename DeviceT>
 template<typename ScalarT>
 KOKKOS_INLINE_FUNCTION