diff --git a/bird/meshing/block_cyl_mesh.py b/bird/meshing/block_cyl_mesh.py
index 518ba23d..6b2aa7da 100644
--- a/bird/meshing/block_cyl_mesh.py
+++ b/bird/meshing/block_cyl_mesh.py
@@ -549,4 +549,3 @@ def writeBlockMeshDict(out_folder, geom_dict, mesh_dict):
     fw.write("defaultPatch\n{type wall;}\n\n")
 
     fw.close()
-
diff --git a/bird/postprocess/early_pred.py b/bird/postprocess/early_pred.py
index 2e491048..84958620 100644
--- a/bird/postprocess/early_pred.py
+++ b/bird/postprocess/early_pred.py
@@ -265,4 +265,3 @@ def bayes_fit(data_dict, num_warmup=1000, num_samples=500):
         data_dict[datf]["std16_real"] = std16_real
 
     return data_dict
-
diff --git a/bird/postprocess/kla_utils.py b/bird/postprocess/kla_utils.py
index c7a5af4d..eca199ac 100644
--- a/bird/postprocess/kla_utils.py
+++ b/bird/postprocess/kla_utils.py
@@ -370,4 +370,3 @@ def print_res_dict(res_dict: dict) -> None:
         logger.info(f"Without data bootstrap")
     logger.info(f"\tkla = {kla_nb*3600:.4g} +/- {kla_err_nb*3600:.4g} [h-1]")
     logger.info(f"\tcstar = {cs_nb:.4g} +/- {cs_err_nb:.4g} [mol/m3]")
-
diff --git a/bird/preprocess/dynamic_mixer/mixing_fvModels.py b/bird/preprocess/dynamic_mixer/mixing_fvModels.py
index 69cc48f4..738fe7a5 100644
--- a/bird/preprocess/dynamic_mixer/mixing_fvModels.py
+++ b/bird/preprocess/dynamic_mixer/mixing_fvModels.py
@@ -54,4 +54,3 @@ def write_fvModel(input_dict, output_folder=".", force_sign=False):
                     write_mixer(mixer, output_folder)
 
     write_end(output_folder)
-
diff --git a/bird/preprocess/species_gen/setup_thermo_prop.py b/bird/preprocess/species_gen/setup_thermo_prop.py
index 31c68dc6..e329c8e6 100644
--- a/bird/preprocess/species_gen/setup_thermo_prop.py
+++ b/bird/preprocess/species_gen/setup_thermo_prop.py
@@ -419,4 +419,3 @@ def write_species_properties(case_folder: str, phase: str = "gas") -> None:
         case_folder, "constant", f"thermophysicalProperties.{phase}"
     )
     write_openfoam_dict(thermo_properties_update, filename=filename)
-
diff --git a/bird/utilities/ofio.py b/bird/utilities/ofio.py
index a1e51b3c..478a3882 100644
--- a/bird/utilities/ofio.py
+++ b/bird/utilities/ofio.py
@@ -1325,19 +1325,25 @@ def repl(match):
                 expr = new_expr
             # Replace 'exp' with 'math.exp'
             expr = expr.replace("exp", "math.exp")
-
+            # Replace 'sin' with 'math.sin'
+            expr = expr.replace("sin", "math.sin")
+            # Replace 'cos' with 'math.cos'
+            expr = expr.replace("cos", "math.cos")
+            # Replace 'tan' with 'math.tan'
+            expr = expr.replace("tan", "math.tan")
+
+            degToRad = lambda x: x * np.pi / 180.0
             try:
-                result = eval(expr, {"math": math})
-            except Exception as e:
-                raise RuntimeError(
-                    f"Error evaluating globalVars expression for {key}: {expr}"
-                ) from e
-
-            # Convert to int if whole number
-            if isinstance(result, float) and result.is_integer():
-                result = int(result)
+                result = eval(expr, {"math": math, "degToRad": degToRad})
+                # Convert to int if whole number
+                if isinstance(result, float) and result.is_integer():
+                    result = int(result)
+                cross_referenced_globalVars_dict[key] = result
 
-            cross_referenced_globalVars_dict[key] = result
+            except Exception as e:
+                logger.warning(
+                    f"Could not evaluate globalVars expression for {key}: {expr}"
+                )
 
     return cross_referenced_globalVars_dict
 
diff --git a/bird/version.py b/bird/version.py
index 6825d034..47726bba 100644
--- a/bird/version.py
+++ b/bird/version.py
@@ -1,3 +1,3 @@
 """Bio reactor design version"""
 
-__version__ = "0.0.50"
+__version__ = "0.0.51"
diff --git a/tests/io/putida_test/constant/globalVars b/tests/io/putida_test/constant/globalVars
new file mode 100644
index 00000000..7f44f721
--- /dev/null
+++ b/tests/io/putida_test/constant/globalVars
@@ -0,0 +1,152 @@
+T0    303.15;  // initial T(K) which stays constant 
+// water thermophysical properties
+muMixLiq #calc "2.414e-5 * pow(10,247.8/($T0 - 140.0))"; //viscosity (Pa.s) of water as a function of T(K)
+CpMixLiq        4181; // J/kg-K
+kThermLiq       0.62; // W/m-K
+rho0MixLiq      1000; // kg/m^3
+sigmaLiq        0.07; //surface tension N/m
+//Wilke-Chang params for diffusion coefficient of a given solute in water (solvent)
+WC_psi          2.6;
+WC_M            18; // kg/kmol
+WC_V_O2         25.6e-3; // m3/kmol molar volume at normal boiling temperature (Treybal 1968) 
+WC_V_H2         14.3e-3; 
+WC_V_CO2        34e-3;
+WC_V_CO         30.7e-3;
+WC_V_CH4        35e-3; // V_b[cm3/mol]=0.285*V_critical^1.048 (Tyn and Calus; ESTIMATING LIQUID MOLAL VOLUME; Processing, Volume 21, Issue 4, Pages 16 - 17)
+WC_V_N2		31.2e-3;
+//****** diffusion coeff ***********
+D_H2 #calc "1.173e-16 * pow($WC_psi * $WC_M,0.5) * $T0 / $muMixLiq / pow($WC_V_H2,0.6)"; 
+D_CO2 #calc "1.173e-16 * pow($WC_psi * $WC_M,0.5) * $T0 / $muMixLiq / pow($WC_V_CO2,0.6)"; 
+D_O2 #calc "1.173e-16 * pow($WC_psi * $WC_M,0.5) * $T0 / $muMixLiq / pow($WC_V_O2,0.6)"; 
+D_CO #calc "1.173e-16 * pow($WC_psi * $WC_M,0.5) * $T0 / $muMixLiq / pow($WC_V_CO,0.6)";
+D_CH4 #calc "1.173e-16 * pow($WC_psi * $WC_M,0.5) * $T0 / $muMixLiq / pow($WC_V_CH4,0.6)";
+D_N2 #calc "1.173e-16 * pow($WC_psi * $WC_M,0.5) * $T0 / $muMixLiq / pow($WC_V_N2,0.6)";
+// D_H2=3.4e-9, D_CO2=2e-9, D_CO=2.16e-9 at 25C
+// Looks like the H2 diffusion is less than litr reported values ~ 4.5e-9 
+//****** Henry coeff ***************
+H_O2_298 0.032;
+DH_O2   1700;
+H_CO2_298 0.83;
+DH_CO2  2400;
+H_CO_298 0.023;
+DH_CO  1300;
+H_H2_298 0.019;
+DH_H2   500;
+H_CH4_298 0.032;
+DH_CH4   1900;
+H_N2_298 0.015;
+DH_N2   1300;
+
+He_H2  #calc "$H_H2_298 * exp($DH_H2 *(1. / $T0 - 1./298.15))";
+He_CO  #calc "$H_CO_298 * exp($DH_CO *(1. / $T0 - 1./298.15))";
+He_CO2 #calc "$H_CO2_298 * exp($DH_CO2 *(1. / $T0 - 1./298.15))";
+He_O2 #calc "$H_O2_298 * exp($DH_O2 *(1. / $T0 - 1./298.15))";
+He_CH4 #calc "$H_CH4_298 * exp($DH_CH4 *(1. / $T0 - 1./298.15))";
+He_N2 #calc "$H_N2_298 * exp($DH_N2 *(1. / $T0 - 1./298.15))";
+//*******inlet gas species mass frac*************
+f_H2    0.0;
+f_CO2   0.0;
+f_CO    0.0;
+f_O2    1.0;
+f_N2    0.0;
+//*********************************
+LeLiqH2 #calc "$kThermLiq / $rho0MixLiq / $D_H2 / $CpMixLiq";
+LeLiqCO #calc "$kThermLiq / $rho0MixLiq / $D_CO / $CpMixLiq";
+LeLiqCO2 #calc "$kThermLiq / $rho0MixLiq / $D_CO2 / $CpMixLiq"; // = 74
+LeLiqO2 #calc "$kThermLiq / $rho0MixLiq / $D_O2 / $CpMixLiq"; 
+LeLiqCH4 #calc "$kThermLiq / $rho0MixLiq / $D_CH4 / $CpMixLiq";
+LeLiqN2 #calc "$kThermLiq / $rho0MixLiq / $D_N2 / $CpMixLiq";
+
+LeLiqMix #calc "$f_N2*$LeLiqN2+$f_CO2*$LeLiqCO2+$f_O2*$LeLiqO2";
+PrMixLiq #calc "$CpMixLiq * $muMixLiq / $kThermLiq";
+// Pr number is ~ 7 for water and ~ 0.7 for air
+//*********************************
+kH2 #calc "$D_H2*$rho0MixLiq*$CpMixLiq*$LeLiqMix";
+PrH2 #calc "$muMixLiq*$CpMixLiq / $kH2";
+
+kCO #calc "$D_CO*$rho0MixLiq*$CpMixLiq*$LeLiqMix";
+PrCO #calc "$muMixLiq*$CpMixLiq / $kCO";
+
+kCO2 #calc "$D_CO2*$rho0MixLiq*$CpMixLiq*$LeLiqMix";
+PrCO2 #calc "$muMixLiq*$CpMixLiq / $kCO2";
+
+kO2 #calc "$D_O2*$rho0MixLiq*$CpMixLiq*$LeLiqMix";
+PrO2 #calc "$muMixLiq*$CpMixLiq / $kO2";
+
+kCH4 #calc "$D_CH4*$rho0MixLiq*$CpMixLiq*$LeLiqMix";
+PrCH4 #calc "$muMixLiq*$CpMixLiq / $kCH4";
+
+kN2 #calc "$D_N2*$rho0MixLiq*$CpMixLiq*$LeLiqMix";
+PrN2 #calc "$muMixLiq*$CpMixLiq / $kN2";
+//*****Gas transport*******************
+muMixGas        1.88e-05;
+PrMixGas        0.7;
+LeGas           1.0;
+//dbubGas         0.001;
+//*********************************
+HtBcol  0.135;
+DiaBcol 0.085;
+LiqHt   0.051;
+LiqHt_limit #calc "0.9*$HtBcol";
+liqVol	0.000337201;
+inletA	0;
+//*********************************
+P0 #calc "101325*0.84974"; //86100; conditions in Golden. CO
+//P0 101325.0;
+//P0 85000.0;
+Pbot #calc "$P0+1000.0*$LiqHt*9.8";
+Pmid #calc "$P0+1000.0*0.5*$LiqHt*9.8";
+
+//bubble column cross-section divided by sparger area
+ArbyAs 2.0; 
+uSupVel 0.02;
+uGas #calc "$Pmid / $Pbot * $ArbyAs * $uSupVel";     
+A_Bcol #calc "3.14159 * $DiaBcol * $DiaBcol / 4.0"; // 0.00567450173
+rho_gas #calc "$Pmid / 287.0 / $T0";
+mflowRate #calc "$uSupVel * $A_Bcol * $rho_gas"; //kg/s
+//at .250 L working fluid, from t=6 to t=22.92h, 15 mL base added
+// -> 4000x scale up for 1000 L working fluid? ignore for now
+//Vb_flowRate #calc "0.015/(22.92 - 6)/3600/1000"; //m^3/s
+Vb_flowRate 8.865e-06;
+
+VVM 1.0; //gas feed per vessel volume, m3 gas/m3 vessel min 
+V_flowRate #calc "$VVM * $liqVol / 60"; // 5.000217e-6 m^3/s => ~300 ccm
+M_flowRate #calc "$V_flowRate * $rho_gas";
+base_flowRate #calc "$Vb_flowRate * $rho0MixLiq"; // 2.4635e-7 kg/s of base solution
+
+// Sparge holes
+ang1 #calc "degToRad(22.5)";
+ang2 #calc "degToRad(45)";
+ang3 #calc "degToRad(67.5)";
+spargeRad  0.0219;
+
+x1 #calc "$spargeRad*cos($ang1)";
+x1neg #calc "-1.0*$x1";
+z1 #calc "$spargeRad*sin($ang1)";
+z1neg #calc "-1.0*$z1";
+x2 #calc "$spargeRad*cos($ang2)";
+x2neg #calc "-1.0*$x2";
+z2 #calc "$spargeRad*sin($ang2)";
+z2neg #calc "-1.0*$z2";
+x3 #calc "$spargeRad*cos($ang3)";
+x3neg #calc "-1.0*$x3";
+z3 #calc "$spargeRad*sin($ang3)";
+z3neg #calc "-1.0*$z3";
+x4 0.0;
+z4 0.0219;
+z4neg -0.0219;
+
+
+
+//*********************************
+//intensity 0.01;
+//k_inlet #calc "1.5 * Foam::pow(($uSupVel * Foam::pow($ArbyAs, 2)), 2) * Foam::pow($intensity, 2)";
+intensity 0.02;
+k_inlet #calc "1.5 * pow($uSupVel * $intensity, 2)";
+//eps_gas_inlet #calc "0.09 * Foam::pow($k_inlet, 1.5) / (($DiaBcol / $ArbyAs) * 0.07)";
+eps_gas_inlet #calc "pow(0.09,0.75) * pow($k_inlet,1.5) / ($DiaBcol / $ArbyAs * 0.07)";
+//eps_liq_inlet #calc "0.09 * 1000 * Foam::pow($k_inlet, 2) * 0.05 / 0.00089";
+eps_liq_inlet #calc "pow(0.09,0.75) * pow($k_inlet,1.5) / ($DiaBcol / $ArbyAs * 0.07)";
+//omega_liq_inlet #calc "1000 * $k_inlet * 0.05 / 0.00089";
+omega_liq_inlet #calc "pow(0.09,-0.25) * pow($k_inlet,0.5) / ($DiaBcol / $ArbyAs * 0.07)";
+//************************************
diff --git a/tests/io/test_case.py b/tests/io/test_case.py
index b0b9c636..e665dfea 100644
--- a/tests/io/test_case.py
+++ b/tests/io/test_case.py
@@ -79,4 +79,3 @@ def test_mesh_vol():
     volumes2, _ = read_cell_volumes(case_folder, time_folder=vol_time)
 
     assert np.linalg.norm(volumes - volumes2) < 1e-12
-
diff --git a/tests/io/test_read_foam_dict.py b/tests/io/test_read_foam_dict.py
index b93f7e23..28d2d58c 100644
--- a/tests/io/test_read_foam_dict.py
+++ b/tests/io/test_read_foam_dict.py
@@ -141,4 +141,3 @@ def test_read_controlDict():
 
     assert foam_dict["writeControl"] == "adjustableRunTime"
     assert foam_dict["maxCo"] == "0.5"
-
diff --git a/tests/io/test_read_foam_fields.py b/tests/io/test_read_foam_fields.py
index e75f1ce4..0eda5014 100644
--- a/tests/io/test_read_foam_fields.py
+++ b/tests/io/test_read_foam_fields.py
@@ -212,4 +212,3 @@ def test_read_mu_liquid():
         os.path.join(case_folder, "thermo:mu.liquid_tmp"),
         os.path.join(case_folder, "80", "thermo:mu.liquid"),
     )
-
diff --git a/tests/io/test_read_global_vars.py b/tests/io/test_read_global_vars.py
index 560e7bd4..f19f8eee 100644
--- a/tests/io/test_read_global_vars.py
+++ b/tests/io/test_read_global_vars.py
@@ -12,8 +12,8 @@ def test_read_global_vars():
     """
     case_folder = os.path.join(
         Path(__file__).parent,
-        "..",
-        "..",
+        " .. ".strip(),
+        " .. ".strip(),
         "tutorial_cases",
         "loop_reactor_mixing",
     )
@@ -105,3 +105,20 @@ def test_read_global_vars():
         < 1e-3
     )
 
+
+def test_read_global_vars_putida():
+    """
+    Test for reading content of `constant/globalVars` for the Putida case
+    """
+    case_folder = os.path.join(Path(__file__).parent, "putida_test")
+    # Read globalVars from case_folder path
+    globalVars_dict = read_global_vars(
+        case_folder=case_folder, cross_ref=False
+    )
+    globalVars_dict = read_global_vars(case_folder=case_folder, cross_ref=True)
+
+    assert abs(globalVars_dict["ang1"] - 22.5 * np.pi / 180) < 1e-12
+    assert (
+        abs(globalVars_dict["z3"] - 0.0219 * np.sin(67.5 * np.pi / 180))
+        < 1e-12
+    )
diff --git a/tests/postprocess/test_cond_mean.py b/tests/postprocess/test_cond_mean.py
index f89e6c6f..bc3218d9 100644
--- a/tests/postprocess/test_cond_mean.py
+++ b/tests/postprocess/test_cond_mean.py
@@ -43,4 +43,3 @@ def test_compute_cond():
         plot_name = sequencePlot(cond, [case_folder], field_name)
         pretty_labels(plot_name, "y [m]", 14)
         plt.close()
-
diff --git a/tests/postprocess/test_post_quantities.py b/tests/postprocess/test_post_quantities.py
index faf8e205..97bafb28 100644
--- a/tests/postprocess/test_post_quantities.py
+++ b/tests/postprocess/test_post_quantities.py
@@ -366,4 +366,3 @@ def test_fitted_kla():
     )
     for time_folder in [str(entry) for entry in range(81, 91)]:
         shutil.rmtree(os.path.join(case_folder, time_folder))
-
diff --git a/tests/preprocess/test_stl_patch.py b/tests/preprocess/test_stl_patch.py
index bb5f9e90..51dbb6b2 100644
--- a/tests/preprocess/test_stl_patch.py
+++ b/tests/preprocess/test_stl_patch.py
@@ -75,4 +75,3 @@ def test_loop_reactor_branch():
         write_boundaries(input_dict, output_folder=tmpdirname)
         axes = plot_stl(os.path.join(tmpdirname, "inlets.stl"))
         pretty_labels("x", "y", zlabel="z", fontsize=14, ax=axes)
-