grid dependent calculation fine

ESMValGroup · Dec 21, 2023 · f5e0377 · f5e0377
1 parent 194b08a
commit f5e0377
Showing 1 changed file with 28 additions and 5 deletions.
diff --git a/esmvaltool/diag_scripts/lifetime/lifetime.py b/esmvaltool/diag_scripts/lifetime/lifetime.py
@@ -632,6 +632,7 @@ def _calculate_rho(self, variables):
         Calculate number density with respect to the given input
         variables.
         """
+        logger.info("Calculate number density (rho)")
         # model levels
         if 'grmassdry' in variables and 'grvol' in variables:
             #if self.z_coord == 'atmosphere_hybrid_sigma_pressure_coordinate'
@@ -693,15 +694,36 @@ def _number_density_dryair_by_press(self, ta, hus, press=None):
         - m_air   Molarmass of Air
         - m_h2o   Molarmass of watervapor
         """
+        logger.info('Calculate number density of dry air by pressure')
 
         if not press:
-            press = broadcast_to_shape(
-                ta.coord('air_pressure').points,
-                ta.shape,
-                ta.coord_dims('air_pressure')
-            )
+            logger.info('Pressure not given')
+
+            if ta.coord('air_pressure').points.shape == ta.shape:
+                press = ta.coord('air_pressure').points
+            else:
+                press = broadcast_to_shape(
+                    ta.coord('air_pressure').points,
+                    ta.shape,
+                    ta.coord_dims('air_pressure')
+                )
+        print(press)
+        print('nominator')
+        nominator = (press * N_A * 10**(-6))
+        print(nominator)
+        print('denominator')
+        denominator = (R * ta * (1. + (self.m_air / self.m_h2o - 1.) * hus))
+        print(denominator)
+        print('rho')
+        rho = nominator / denominator
+        print(rho)
+        sys.exit(1)
+
         rho = ((press * N_A * 10**(-6)) /
                (R * ta * (1. + (self.m_air / self.m_h2o - 1.) * hus)))
+
+        print(rho)
+
         # correct metadata
         rho.var_name = 'rho'
         rho.units = 'cm-3'
@@ -724,6 +746,7 @@ def _number_density_dryair_by_grid(self, grmassdry, grvol):
         - N_A    Avogrado constant from scipy
         - m_air   Molarmass of Air
         """
+        logger.info('Calculate number density of dry air by grid information')
         rho = ((grmassdry / grvol)
                * (N_A / self.m_air) * 10**(-3))  # [ 1 / cm^3 ]
         # correct metadata