diff --git a/README.md b/README.md
index 6d8301d2..5f46d1bd 100644
--- a/README.md
+++ b/README.md
@@ -51,7 +51,10 @@ To cite BiRD, please use these articles on [CO2 interphase mass transfer](https:
     title={Bayesian calibration of bubble size dynamics applied to \ce{CO2} gas fermenters},
     author={Hassanaly, Malik and Parra-Alvarez, John M. and Rahimi, Mohammad J., Municchi, Federico and Sitaraman, Hariswaran},
     journal={Chemical Engineering Research and Design},
+    volume={215},
+    pages={312--328},
     year={2025},
+    publisher={Elsevier}
   }
 
 @article{rahimi2018computational,
diff --git a/bird/__init__.py b/bird/__init__.py
index 1c16af5e..0d1d96eb 100644
--- a/bird/__init__.py
+++ b/bird/__init__.py
@@ -8,6 +8,7 @@
 BIRD_MESH_DIR = os.path.join(BIRD_DIR, "meshing")
 BIRD_POST_DIR = os.path.join(BIRD_DIR, "postprocess")
 BIRD_PRE_DIR = os.path.join(BIRD_DIR, "preprocess")
+BIRD_CONST_DIR = os.path.join(BIRD_DIR, "constants")
 BIRD_BLOCK_CYL_MESH_TEMP_DIR = os.path.join(
     BIRD_MESH_DIR, "block_cyl_mesh_templates"
 )
diff --git a/bird/constants/ch4.yaml b/bird/constants/ch4.yaml
new file mode 100644
index 00000000..82fab403
--- /dev/null
+++ b/bird/constants/ch4.yaml
@@ -0,0 +1,33 @@
+name: CH4
+
+specie:
+    molWeight: 16.04303
+    elements:
+        C: 1
+        H: 4 
+
+liquid:
+    LeLiq: '#calc "$kThermLiq / $rho0MixLiq / $D_CH4 / $CpMixLiq"'
+    thermodynamics:
+        Hf: -1.5879e+07
+
+gas:
+    in-H2O:
+        WC_V: "35e-3" # m3/kmol molar volume at normal boiling temperature (Treybal 1968)
+        D: '#calc "1.173e-16 * pow($WC_psi * $WC_M,0.5) * $T0 / $muMixLiq / pow($WC_V_CH4,0.6)' 
+    H_298: "0.032"
+    DH: "1900"
+    He: '#calc "$H_CH4_298 * exp($DH_CH4 *(1. / $T0 - 1./298.15))"'
+    k: '#calc "$D_CH4*$rho0MixLiq*$CpMixLiq*$LeLiqMix"'
+    Pr: '#calc "$muMixLiq*$CpMixLiq / $kCH4"'
+  
+    thermodynamics:
+        Tlow: 200
+        Thigh: 3500
+        Tcommon: 1000
+        highCpCoeffs: 0.074851495 0.0133909467 -5.73285809e-06 1.22292535e-09 -1.0181523e-13 -9468.34459 18.437318
+        lowCpCoeffs: 5.14987613 -0.0136709788 4.91800599e-05 -4.84743026e-08 1.66693956e-11 -10246.6476 -4.64130376
+    transport:
+        As: 1.512e-06
+        Ts: 120
+   
diff --git a/bird/constants/co.yaml b/bird/constants/co.yaml
new file mode 100644
index 00000000..88d1d2bc
--- /dev/null
+++ b/bird/constants/co.yaml
@@ -0,0 +1,34 @@
+name: CO
+
+specie:
+    molWeight: 28.01055
+    elements:
+        C: 1
+        O: 1
+
+liquid:
+    LeLiq: '#calc "$kThermLiq / $rho0MixLiq / $D_CO / $CpMixLiq"'
+    thermodynamics:
+        Hf: -1.5879e+07
+
+gas:
+    in-H2O:
+        WC_V: "30.7e-3" # m3/kmol molar volume at normal boiling temperature (Treybal 1968)
+        D: '#calc "1.173e-16 * pow($WC_psi * $WC_M,0.5) * $T0 / $muMixLiq / pow($WC_V_CO,0.6)' 
+    H_298: "0.023"
+    DH: "1300"
+    He: '#calc "$H_CO_298 * exp($DH_CO *(1. / $T0 - 1./298.15))"'
+    k: '#calc "$D_CO*$rho0MixLiq*$CpMixLiq*$LeLiqMix"'
+    Pr: '#calc "$muMixLiq*$CpMixLiq / $kCO"'
+   
+    thermodynamics:
+        Tlow: 200
+        Thigh: 3500
+        Tcommon: 1000
+        highCpCoeffs: 2.71518561 0.00206252743 -9.98825771e-07 2.30053008e-10 -2.03647716e-14 -14151.8724 7.81868772
+        lowCpCoeffs: 3.57953347 -0.00061035368 1.01681433e-06 9.07005884e-10 -9.04424499e-13 -14344.086 3.50840928
+    transport:
+        As:  1.512e-06
+        Ts: 120
+
+   
diff --git a/bird/constants/co2.yaml b/bird/constants/co2.yaml
new file mode 100644
index 00000000..aed9085c
--- /dev/null
+++ b/bird/constants/co2.yaml
@@ -0,0 +1,34 @@
+name: CO2
+
+specie:
+    molWeight: 44.00995
+    elements:
+        C: 1
+        O: 2 
+
+liquid:
+    LeLiq: '#calc "$kThermLiq / $rho0MixLiq / $D_CO2 / $CpMixLiq"'
+    thermodynamics:
+        Hf: -1.5879e+07
+
+gas:
+    in-H2O:
+        WC_V: "34e-3" # m3/kmol molar volume at normal boiling temperature (Treybal 1968)
+        D: '#calc "1.173e-16 * pow($WC_psi * $WC_M,0.5) * $T0 / $muMixLiq / pow($WC_V_CO2,0.6)' 
+    H_298: "0.83"
+    DH: "2400"
+    He: '#calc "$H_CO2_298 * exp($DH_CO2 *(1. / $T0 - 1./298.15))"'
+    k: '#calc "$D_CO2*$rho0MixLiq*$CpMixLiq*$LeLiqMix"'
+    Pr: '#calc "$muMixLiq*$CpMixLiq / $kCO2"'
+   
+    thermodynamics:
+        Tlow: 200
+        Thigh: 3500
+        Tcommon: 1000
+        highCpCoeffs: 3.85746029 0.00441437026 -2.21481404e-06 5.23490188e-10 -4.72084164e-14 -48759.166 2.27163806
+        lowCpCoeffs: 2.35677352 0.00898459677 -7.12356269e-06 2.45919022e-09 -1.43699548e-13 -48371.9697 9.90105222
+    transport:
+        As: 1.652e-06
+        Ts: 273
+    
+     
diff --git a/bird/constants/h2.yaml b/bird/constants/h2.yaml
new file mode 100644
index 00000000..f134129a
--- /dev/null
+++ b/bird/constants/h2.yaml
@@ -0,0 +1,32 @@
+name: H2
+
+specie:
+    molWeight: 2.01594
+    elements:
+        H: 2
+
+liquid:
+    LeLiq: '#calc "$kThermLiq / $rho0MixLiq / $D_H2 / $CpMixLiq"'
+    thermodynamics:
+        Hf: -1.5879e+07
+
+gas:
+    in-H2O:
+        WC_V: "14.3e-3" # m3/kmol molar volume at normal boiling temperature (Treybal 1968)
+        D: '#calc "1.173e-16 * pow($WC_psi * $WC_M,0.5) * $T0 / $muMixLiq / pow($WC_V_H2,0.6)' 
+    H_298: "0.019"
+    DH: "500"
+    He: '#calc "$H_H2_298 * exp($DH_H2 *(1. / $T0 - 1./298.15))"'
+    k: '#calc "$D_H2*$rho0MixLiq*$CpMixLiq*$LeLiqMix"'
+    Pr: '#calc "$muMixLiq*$CpMixLiq / $kH2"'
+  
+    thermodynamics:
+        Tlow: 200
+        Thigh: 3500
+        Tcommon: 1000
+        highCpCoeffs: 3.3372792 -4.94024731e-05 4.99456778e-07 -1.79566394e-10 2.00255376e-14 -950.158922 -3.20502331
+        lowCpCoeffs: 2.34433112 0.00798052075 -1.9478151e-05 2.01572094e-08 -7.37611761e-12 -917.935173 0.683010238
+    transport:
+        As: 6.362e-07
+        Ts: 72
+   
diff --git a/bird/constants/h2o.yaml b/bird/constants/h2o.yaml
new file mode 100644
index 00000000..6a501feb
--- /dev/null
+++ b/bird/constants/h2o.yaml
@@ -0,0 +1,30 @@
+name: H2O
+
+specie:
+     molWeight: 18.01534
+     elements:
+         H: 2
+         O: 1
+
+liquid:
+    CpMixLiq: "4181"
+    muMixLiq: '#calc "2.414e-5 * pow(10,247.8/($T0 - 140.0))"' # Viscosity (Pa.s)
+    kThermLiq: "0.62" # thermal conductivity (W/mK)
+    rho0MixLiq: "1000" # density (kg/m3)
+    sigmaLiq: "0.07" # surface tension (N/m)
+    WC_psi: "2.6"
+    thermodynamics:
+         Hf:  -1.5879e+07
+
+
+gas:
+    thermodynamics:
+        Tlow: 200
+        Thigh: 3500
+        Tcommon: 1000
+        highCpCoeffs: 3.03399249 0.00217691804 -1.64072518e-07 -9.7041987e-11 1.68200992e-14 -30004.2971 4.9667701
+        lowCpCoeffs: 4.19864056 -0.0020364341 6.52040211e-06 -5.48797062e-09 1.77197817e-12 -30293.7267 -0.849032208 
+        transport:
+            As: 1.456e-06
+            Ts: 273
+        
diff --git a/bird/constants/n2.yaml b/bird/constants/n2.yaml
new file mode 100644
index 00000000..4ec7d79c
--- /dev/null
+++ b/bird/constants/n2.yaml
@@ -0,0 +1,28 @@
+name: N2
+
+specie:
+    molWeight: 28.0134
+    elements:
+        N: 2 
+
+liquid:
+    LeLiq: '#calc "$kThermLiq / $rho0MixLiq / $D_N2 / $CpMixLiq"'
+
+gas:
+    in-H2O:
+        WC_V: "31.2e-3" # m3/kmol molar volume at normal boiling temperature (Treybal 1968)
+        D: '#calc "1.173e-16 * pow($WC_psi * $WC_M,0.5) * $T0 / $muMixLiq / pow($WC_V_N2,0.6)' 
+    H_298: "0.015"
+    DH: "1300"
+    He: '#calc "$H_N2_298 * exp($DH_N2 *(1. / $T0 - 1./298.15))"'
+    k: '#calc "$D_N2*$rho0MixLiq*$CpMixLiq*$LeLiqMix"'
+    Pr: '#calc "$muMixLiq*$CpMixLiq / $kN2"'
+   
+    thermodynamics:
+        Tlow: 250
+        Thigh: 5000
+        highCpCoeffs: 2.92664 0.0014879768 -5.68476e-07 1.0097038e-10 -6.753351e-15 -922.7977 5.980528
+        lowCpCoeffs:  3.298677 0.0014082404 -3.963222e-06 5.641515e-09 -2.444854e-12 -1020.8999 3.950372
+    transport:
+        As: 1.688e-06
+        Ts: 273
diff --git a/bird/constants/o2.yaml b/bird/constants/o2.yaml
new file mode 100644
index 00000000..8d9704ec
--- /dev/null
+++ b/bird/constants/o2.yaml
@@ -0,0 +1,32 @@
+name: O2
+
+specie:
+    molWeight: 31.999
+    elements:
+        O: 2 
+
+
+liquid:
+    LeLiq: '#calc "$kThermLiq / $rho0MixLiq / $D_O2 / $CpMixLiq"'
+    thermodynamics:
+        Hf: -365639.0
+
+gas:
+    in-H2O:
+        WC_V: "25.6e-3" # m3/kmol molar volume at normal boiling temperature (Treybal 1968)
+        D: '#calc "1.173e-16 * pow($WC_psi * $WC_M,0.5) * $T0 / $muMixLiq / pow($WC_V_O2,0.6)' 
+    H_298: "0.032"
+    DH: "1700"
+    He: '#calc "$H_O2_298 * exp($DH_O2 *(1. / $T0 - 1./298.15))"'
+    k: '#calc "$D_O2*$rho0MixLiq*$CpMixLiq*$LeLiqMix"'
+    Pr: '#calc "$muMixLiq*$CpMixLiq / $kO2"'
+   
+    thermodynamics:
+        Tlow: 200
+        Thigh: 3500
+        Tcommon: 1000
+        highCpCoeffs: 3.28254 0.00148309 -7.57967e-07 2.09471e-10 -2.16718e-14 -1038.31 5.56079
+        lowCpCoeffs: 3.78246 -0.00299673 9.8473e-06 -9.6813e-09 3.24373e-12 -1063.94 3.65768
+    transport: 
+        As: 1.948e-06
+        Ts: 273
diff --git a/bird/preprocess/sim_setup/global_vars.py b/bird/preprocess/sim_setup/global_vars.py
new file mode 100644
index 00000000..33401a3b
--- /dev/null
+++ b/bird/preprocess/sim_setup/global_vars.py
@@ -0,0 +1,592 @@
+import os
+import re
+from collections import defaultdict
+
+from ruamel.yaml import YAML
+
+from bird import BIRD_CONST_DIR
+from bird.utilities.ofio import read_properties
+
+
+def parse_yaml(filename: str):
+    yaml = YAML()
+    with open(filename, "r+") as f:
+        spec = yaml.load(f)
+
+    return spec
+
+
+def global_vars_to_dict(filename: str):
+    with open(filename, "r+") as f:
+        lines = f.readlines()
+    global_vars_dict = {}
+    for line in lines:
+        if not line.startswith("//") and not line.startswith("\n"):
+            line_l = line.split()
+            name = line_l[0]
+            val = line[len(name) + 1 :]
+            ind_end = val.index(";")
+            val = val[:ind_end]
+            global_vars_dict[name] = val
+
+    return global_vars_dict
+
+
+def add_liquid_properties(case_dir, global_vars_dict):
+    print("WARNING: assuming liquid phase is water")
+    h2o_prop = parse_yaml(os.path.join(BIRD_CONST_DIR, "h2o.yaml"))
+    CpMixLiq = h2o_prop["liquid"]["CpMixLiq"]
+    muMixLiq = h2o_prop["liquid"]["muMixLiq"]
+    kThermLiq = h2o_prop["liquid"]["kThermLiq"]
+    rho0MixLiq = h2o_prop["liquid"]["rho0MixLiq"]
+    sigmaLiq = h2o_prop["liquid"]["sigmaLiq"]
+    WC_psi = h2o_prop["liquid"]["WC_psi"]
+    WC_M = h2o_prop["specie"]["molWeight"]
+
+    if "CpMixLiq" not in global_vars_dict:
+        global_vars_dict["CpMixLiq"] = CpMixLiq
+    else:
+        print(
+            f"WARNING: CpMixLiq stays at {global_vars_dict['CpMixLiq']} instead of {CpMixLiq}"
+        )
+    if "muMixLiq" not in global_vars_dict:
+        global_vars_dict["muMixLiq"] = muMixLiq
+    else:
+        print(
+            f"WARNING: muMixLiq stays at {global_vars_dict['muMixLiq']} instead of {muMixLiq}"
+        )
+    if "kThermLiq" not in global_vars_dict:
+        global_vars_dict["kThermLiq"] = kThermLiq
+    else:
+        print(
+            f"WARNING: kThermLiq stays at {global_vars_dict['kThermLiq']} instead of {kThermLiq}"
+        )
+    if "rho0MixLiq" not in global_vars_dict:
+        global_vars_dict["rho0MixLiq"] = rho0MixLiq
+    else:
+        print(
+            f"WARNING: rho0MixLiq stays at {global_vars_dict['rho0MixLiq']} instead of {rho0MixLiq}"
+        )
+    if "sigmaLiq" not in global_vars_dict:
+        global_vars_dict["sigmaLiq"] = sigmaLiq
+    else:
+        print(
+            f"WARNING: sigmaLiq stays at {global_vars_dict['sigmaLiq']} instead of {sigmaLiq}"
+        )
+    if "WC_psi" not in global_vars_dict:
+        global_vars_dict["WC_psi"] = WC_psi
+    else:
+        print(
+            f"WARNING: WC_psi stays at {global_vars_dict['WC_psi']} instead of {WC_psi}"
+        )
+    if "WC_M" not in global_vars_dict:
+        global_vars_dict["WC_M"] = WC_M
+    else:
+        print(
+            f"WARNING: WC_M stays at {global_vars_dict['WC_M']} instead of {WC_M}"
+        )
+
+    return global_vars_dict
+
+
+def add_gas_properties(case_dir, global_vars_dict):
+    gas_properties = read_properties(
+        os.path.join(case_dir, "constant", "thermophysicalProperties.gas")
+    )
+    list_gas_spec = list(
+        set(gas_properties["species"] + [gas_properties["defaultSpecie"]])
+    )
+    try:
+        list_gas_spec.remove("water")
+    except ValueError:
+        pass
+    try:
+        list_gas_spec.remove("H2O")
+    except ValueError:
+        pass
+
+    print(f"INFO: gas species = {list_gas_spec}")
+    for spec in list_gas_spec:
+        if spec == "water":
+            spec_prop = parse_yaml(os.path.join(BIRD_CONST_DIR, "h2o.yaml"))
+        else:
+            spec_prop = parse_yaml(
+                os.path.join(BIRD_CONST_DIR, f"{spec.lower()}.yaml")
+            )
+        name = spec_prop["name"]
+        WC_V = spec_prop["gas"]["in-H2O"]["WC_V"]
+        D = spec_prop["gas"]["in-H2O"]["D"]
+        H_298 = spec_prop["gas"]["H_298"]
+        DH = spec_prop["gas"]["DH"]
+        He = spec_prop["gas"]["He"]
+        LeLiq = spec_prop["liquid"]["LeLiq"]
+        k = spec_prop["gas"]["k"]
+        Pr = spec_prop["gas"]["Pr"]
+        if not f"WC_V_{name}" in global_vars_dict:
+            global_vars_dict[f"WC_V_{name}"] = WC_V
+        else:
+            print(
+                f"WARNING: WC_V_{name} stays at {global_vars_dict[f'WC_V_{name}']} instead of {WC_V}"
+            )
+        if not f"D_{name}" in global_vars_dict:
+            global_vars_dict[f"D_{name}"] = D
+        else:
+            print(
+                f"WARNING: D_{name} stays at {global_vars_dict[f'D_{name}']} instead of {D}"
+            )
+        if not f"H_{name}_298" in global_vars_dict:
+            global_vars_dict[f"H_{name}_298"] = H_298
+        else:
+            print(
+                f"WARNING: H_{name}_298 stays at {global_vars_dict[f'H_{name}_298']} instead of {H_298}"
+            )
+        if not f"DH_{name}" in global_vars_dict:
+            global_vars_dict[f"DH_{name}"] = DH
+        else:
+            print(
+                f"WARNING: DH_{name} stays at {global_vars_dict[f'DH_{name}']} instead of {DH}"
+            )
+        if not f"He_{name}" in global_vars_dict:
+            global_vars_dict[f"He_{name}"] = He
+        else:
+            print(
+                f"WARNING: He_{name} stays at {global_vars_dict[f'He_{name}']} instead of {He}"
+            )
+        if not f"LeLiq{name}" in global_vars_dict:
+            global_vars_dict[f"LeLiq{name}"] = LeLiq
+        else:
+            print(
+                f"WARNING: LeLiq{name} stays at {global_vars_dict[f'LeLiq{name}']} instead of {LeLiq}"
+            )
+        if not f"k{name}" in global_vars_dict:
+            global_vars_dict[f"k{name}"] = k
+        else:
+            print(
+                f"WARNING: k{name} stays at {global_vars_dict[f'k{name}']} instead of {k}"
+            )
+        if not f"Pr{name}" in global_vars_dict:
+            global_vars_dict[f"Pr{name}"] = Pr
+        else:
+            print(
+                f"WARNING: Pr{name} stays at {global_vars_dict[f'Pr{name}']} instead of {Pr}"
+            )
+
+    return global_vars_dict, list_gas_spec
+
+
+def check_globalVars(case_dir, global_vars_dict, list_gas_spec):
+    globalVar_filename = os.path.join(case_dir, "constant", "globalVars")
+    list_supported_key = ["T0", "VVM"]
+    list_supported_key += [
+        "CpMixLiq",
+        "muMixLiq",
+        "kThermLiq",
+        "rho0MixLiq",
+        "sigmaLiq",
+    ]
+    list_supported_key += ["WC_psi", "WC_M"]
+    for spec in list_gas_spec:
+        list_supported_key += [f"WC_V_{spec}"]
+    for spec in list_gas_spec:
+        list_supported_key += [f"D_{spec}"]
+    for spec in list_gas_spec:
+        list_supported_key += [f"H_{spec}_298"]
+    for spec in list_gas_spec:
+        list_supported_key += [f"DH_{spec}"]
+    for spec in list_gas_spec:
+        list_supported_key += [f"He_{spec}"]
+    for spec in list_gas_spec:
+        list_supported_key += [f"f_{spec}"]
+        list_supported_key += [f"x_{spec}"]
+    list_supported_key += ["f_H2O", "f_water"]
+    list_supported_key += [
+        "inletA",
+        "liqVol",
+        "alphaGas",
+        "alphaLiq",
+        "uGasPhase",
+        "P0",
+        "Pbot",
+        "Pmid",
+        "ArbyAs",
+        "uSupVel",
+        "uGas",
+        "rho_gas",
+    ]
+    list_supported_key += [
+        "uGasPhase_Sup",
+        "uLiqPhase_Sup",
+        "uLiqPhase",
+        "mflowRateLiq",
+        "mflowRateGas",
+        "mflowRate",
+    ]
+    for spec in list_gas_spec:
+        list_supported_key += [f"LeLiq{spec}"]
+    list_supported_key += ["PrMixLiq", "LeLiqMix"]
+    for spec in list_gas_spec:
+        list_supported_key += [f"k{spec}", f"Pr{spec}"]
+    list_supported_key += [
+        "l_scale",
+        "intensity",
+        "k_inlet_gas",
+        "k_inlet_liq",
+        "eps_inlet_gas",
+        "eps_inlet_liq",
+        "omega_inlet_gas",
+        "omega_inlet_liq",
+    ]
+    list_supported_key += [
+        "HtBcol",
+        "DiaBcol",
+        "LiqHt",
+        "NPS",
+        "NPD",
+        "NPY",
+        "A_Bcol",
+        "V_flowRate",
+    ]
+
+    keys_unknown = []
+    for key in global_vars_dict:
+        if key not in list_supported_key:
+            print(
+                f"WARNING: no specific treatment for key {key}, will be written at the beginning"
+            )
+            keys_unknown.append(key)
+    return keys_unknown
+
+
+def write_macro_prop(case_dir, global_vars_dict, list_gas_spec):
+    globalVar_filename = os.path.join(case_dir, "constant", "globalVars")
+    try:
+        os.remove(globalVar_filename)
+    except FileNotFoundError:
+        pass
+
+    with open(globalVar_filename, "w+") as f:
+        if "T0" in global_vars_dict:
+            f.write(
+                f"T0      {global_vars_dict['T0']};//initial T(K) which stays constant\n"
+            )
+            print(f"WARNING: Assuming 300K for T0")
+            f.write(f"T0      300;//initial T(K) which stays constant\n")
+        if "VVM" in global_vars_dict:
+            f.write(f"VVM      {global_vars_dict['VVM']};\n")
+
+
+def write_liq_prop(case_dir, global_vars_dict, list_gas_spec):
+    globalVar_filename = os.path.join(case_dir, "constant", "globalVars")
+    with open(globalVar_filename, "a+") as f:
+        f.write("//****water Liquid properties**************\n")
+        if "CpMixLiq" in global_vars_dict:
+            f.write(f"CpMixLiq      {global_vars_dict['CpMixLiq']};\n")
+        if "muMixLiq" in global_vars_dict:
+            f.write(
+                f"muMixLiq      {global_vars_dict['muMixLiq']};//viscosity (Pa.s) of water as a function of T(K) \n"
+            )
+        if "kThermLiq" in global_vars_dict:
+            f.write(
+                f"kThermLiq     {global_vars_dict['kThermLiq']};// W/m-K\n"
+            )
+        if "rho0MixLiq" in global_vars_dict:
+            f.write(
+                f"rho0MixLiq    {global_vars_dict['rho0MixLiq']};// kg/m^3\n"
+            )
+        if "sigmaLiq" in global_vars_dict:
+            f.write(
+                f"sigmaLiq      {global_vars_dict['sigmaLiq']};//surface tension N/m\n"
+            )
+
+
+def write_wc_prop(case_dir, global_vars_dict, list_gas_spec):
+    globalVar_filename = os.path.join(case_dir, "constant", "globalVars")
+    with open(globalVar_filename, "a+") as f:
+        f.write(
+            "//****Wilke-Chang params for diffusion coefficient of a given solute in water (solvent)\n"
+        )
+        if "WC_psi" in global_vars_dict:
+            f.write(f"WC_psi      {global_vars_dict['WC_psi']};\n")
+        if "WC_M" in global_vars_dict:
+            f.write(f"WC_M      {global_vars_dict['WC_M']};// kg/kmol\n")
+        for spec in list_gas_spec:
+            key = f"WC_V_{spec}"
+            if key in global_vars_dict:
+                f.write(f"{key}     {global_vars_dict[key]};// m3/kmol\n")
+
+
+def write_diff_prop(case_dir, global_vars_dict, list_gas_spec):
+    globalVar_filename = os.path.join(case_dir, "constant", "globalVars")
+    with open(globalVar_filename, "a+") as f:
+        f.write("//****** diffusion coeff ***********\n")
+        for spec in list_gas_spec:
+            key = f"D_{spec}"
+            if key in global_vars_dict:
+                f.write(f"{key}     {global_vars_dict[key]};// m3/kmol\n")
+
+
+def write_henry_prop(case_dir, global_vars_dict, list_gas_spec):
+    globalVar_filename = os.path.join(case_dir, "constant", "globalVars")
+    with open(globalVar_filename, "a+") as f:
+        f.write("//****** Henry coeff ***********\n")
+        for spec in list_gas_spec:
+            key = f"H_{spec}_298"
+            if key in global_vars_dict:
+                f.write(f"{key}     {global_vars_dict[key]};\n")
+            key = f"DH_{spec}"
+            if key in global_vars_dict:
+                f.write(f"{key}     {global_vars_dict[key]};\n")
+        for spec in list_gas_spec:
+            key = f"He_{spec}"
+            if key in global_vars_dict:
+                f.write(f"{key}     {global_vars_dict[key]};\n")
+
+
+def write_inlet_prop(case_dir, global_vars_dict, list_gas_spec):
+    globalVar_filename = os.path.join(case_dir, "constant", "globalVars")
+    with open(globalVar_filename, "a+") as f:
+        f.write("//****** inlet gas frac ***********\n")
+        for spec in list_gas_spec:
+            key = f"x_{spec}"
+            if key in global_vars_dict:
+                f.write(f"{key}     {global_vars_dict[key]};\n")
+        for spec in list_gas_spec:
+            key = f"f_{spec}"
+            if key in global_vars_dict:
+                f.write(f"{key}     {global_vars_dict[key]};\n")
+        key = f"f_H2O"
+        if key in global_vars_dict:
+            f.write(f"{key}     {global_vars_dict[key]};\n")
+        key = f"f_water"
+        if key in global_vars_dict:
+            f.write(f"{key}     {global_vars_dict[key]};\n")
+        key_list = [
+            "uGasPhase_Sup",
+            "uLiqPhase_Sup",
+            "alphaGas",
+            "alphaLiq",
+            "uGasPhase",
+            "uLiqPhase",
+            "mflowRateLiq",
+            "mflowRateGas",
+            "mflowRate",
+        ]
+        key_list += ["inletA", "liqVol", "ArbyAs", "uSupVel", "uGas"]
+        for key in key_list:
+            if key in global_vars_dict:
+                val = global_vars_dict[key]
+                f.write(f"{key}     {val};\n")
+
+
+def write_cond_prop(case_dir, global_vars_dict, list_gas_spec):
+    globalVar_filename = os.path.join(case_dir, "constant", "globalVars")
+    with open(globalVar_filename, "a+") as f:
+        f.write("//*****************\n")
+        for ispec, spec in enumerate(list_gas_spec):
+            key = f"k{spec}"
+            if key in global_vars_dict:
+                f.write(f"{key}     {global_vars_dict[key]};\n")
+            key = f"Pr{spec}"
+            if key in global_vars_dict:
+                f.write(f"{key}     {global_vars_dict[key]};\n")
+            if ispec < len(list_gas_spec) - 1:
+                f.write("\n")
+
+
+def write_Le_prop(case_dir, global_vars_dict, list_gas_spec):
+    globalVar_filename = os.path.join(case_dir, "constant", "globalVars")
+    with open(globalVar_filename, "a+") as f:
+        f.write("//*****************\n")
+        for ispec, spec in enumerate(list_gas_spec):
+            key = f"LeLiq{spec}"
+            if key in global_vars_dict:
+                f.write(f"{key}     {global_vars_dict[key]};\n")
+        if "LeLiqMix" in global_vars_dict:
+            f.write(f"LeLiqMix     {global_vars_dict['LeLiqMix']};\n")
+        else:
+            LeLiqMix_val = "#calc "
+            for ispec, spec in enumerate(list_gas_spec):
+                if ispec == 0:
+                    LeLiqMix_val += f'"$f_{spec}*$LeLiq{spec}'
+                else:
+                    LeLiqMix_val += f"+$f_{spec}*$LeLiq{spec}"
+                if ispec == len(list_gas_spec) - 1:
+                    LeLiqMix_val += '"'
+            f.write(f"LeLiqMix     {LeLiqMix_val};\n")
+        if "PrMixLiq" in global_vars_dict:
+            f.write(f"PrMixLiq     {global_vars_dict['PrMixLiq']};\n")
+        else:
+            f.write(
+                f'PrMixLiq     #calc "$CpMixLiq * $muMixLiq / $kThermLiq";\n'
+            )
+
+
+def write_turb_prop(case_dir, global_vars_dict, list_gas_spec):
+    globalVar_filename = os.path.join(case_dir, "constant", "globalVars")
+    with open(globalVar_filename, "a+") as f:
+        f.write("//*****************\n")
+        if "l_scale" in global_vars_dict:
+            f.write(f"l_scale     {global_vars_dict['l_scale']};\n")
+        if "intensity" in global_vars_dict:
+            f.write(f"intensity   {global_vars_dict['intensity']};\n")
+        if "k_inlet_gas" in global_vars_dict:
+            f.write(f"k_inlet_gas {global_vars_dict['k_inlet_gas']};\n")
+        else:
+            if (
+                "uGasPhase" in global_vars_dict
+                and "intensity" in global_vars_dict
+            ):
+                f.write(
+                    f'k_inlet_gas #calc "1.5 * Foam::pow(($uGasPhase), 2) * Foam::pow($intensity, 2)";\n'
+                )
+                global_vars_dict["k_inlet_gas"] = (
+                    '#calc "1.5 * Foam::pow(($uGasPhase), 2) * Foam::pow($intensity, 2)"'
+                )
+            else:
+                print("WARNING: no k_inlet_gas could be printed")
+        if "k_inlet_liq" in global_vars_dict:
+            f.write(f"k_inlet_liq   {global_vars_dict['k_inlet_liq']};\n")
+        else:
+            if (
+                "uGasPhase" in global_vars_dict
+                and "intensity" in global_vars_dict
+            ):
+                f.write(
+                    f'k_inlet_liq #calc "1.5 * Foam::pow(($uGasPhase), 2) * Foam::pow($intensity, 2)";\n'
+                )
+                global_vars_dict["k_inlet_liq"] = (
+                    '#calc "1.5 * Foam::pow(($uGasPhase), 2) * Foam::pow($intensity, 2)"'
+                )
+            else:
+                print("WARNING: no k_inlet_liq could be printed")
+        if "eps_inlet_gas" in global_vars_dict:
+            f.write(f"eps_inlet_gas {global_vars_dict['eps_inlet_gas']};\n")
+        else:
+            if (
+                "k_inlet_gas" in global_vars_dict
+                and "l_scale" in global_vars_dict
+            ):
+                f.write(
+                    f'eps_inlet_gas #calc "pow(0.09,0.75) * Foam::pow($k_inlet_gas, 1.5) / ($l_scale * 0.07)";\n'
+                )
+                global_vars_dict[f"eps_inlet_gas"] = (
+                    '#calc "pow(0.09,0.75) * Foam::pow($k_inlet_gas, 1.5) / ($l_scale * 0.07)"'
+                )
+            else:
+                print("WARNING: no eps_inlet_gas could be printed")
+        if "eps_inlet_liq" in global_vars_dict:
+            f.write(f"eps_inlet_liq   {global_vars_dict['eps_inlet_liq']};\n")
+        else:
+            if (
+                "k_inlet_liq" in global_vars_dict
+                and "l_scale" in global_vars_dict
+            ):
+                f.write(
+                    f'eps_inlet_liq #calc "pow(0.09,0.75) * Foam::pow($k_inlet_liq, 1.5) / ($l_scale * 0.07)";\n'
+                )
+                global_vars_dict[f"eps_inlet_liq"] = (
+                    '#calc "pow(0.09,0.75) * Foam::pow($k_inlet_liq, 1.5) / ($l_scale * 0.07)"'
+                )
+            else:
+                print("WARNING: no eps_inlet_liq could be printed")
+        if "omega_inlet_gas" in global_vars_dict:
+            f.write(
+                f"omega_inlet_gas {global_vars_dict['omega_inlet_gas']};\n"
+            )
+        else:
+            if (
+                "k_inlet_gas" in global_vars_dict
+                and "l_scale" in global_vars_dict
+            ):
+                f.write(
+                    f'omega_inlet_gas #calc "pow(0.09,-0.25) * pow($k_inlet_gas,0.5) / ($l_scale * 0.07)";\n'
+                )
+                global_vars_dict[f"omega_inlet_gas"] = (
+                    '#calc "pow(0.09,-0.25) * pow($k_inlet_gas,0.5) / ($l_scale * 0.07)"'
+                )
+            else:
+                print("WARNING: no omega_inlet_gas could be printed")
+
+        if "omega_inlet_liq" in global_vars_dict:
+            f.write(
+                f"omega_inlet_liq   {global_vars_dict['omega_inlet_liq']};\n"
+            )
+        else:
+            if (
+                "k_inlet_liq" in global_vars_dict
+                and "l_scale" in global_vars_dict
+            ):
+                f.write(
+                    f'omega_inlet_liq #calc "pow(0.09,-0.25) * pow($k_inlet_liq,0.5) / ($l_scale * 0.07)";\n'
+                )
+                global_vars_dict[f"omega_inlet_liq"] = (
+                    '#calc "pow(0.09,-0.25) * pow($k_inlet_liq,0.5) / ($l_scale * 0.07)"'
+                )
+            else:
+                print("WARNING: no omega_inlet_liq could be printed")
+
+
+def write_unknown_prop(
+    case_dir, global_vars_dict, list_gas_spec, unknown_keys
+):
+    globalVar_filename = os.path.join(case_dir, "constant", "globalVars")
+    with open(globalVar_filename, "a+") as f:
+        f.write("//*****************\n")
+        for key in unknown_keys:
+            if key in global_vars_dict:
+                val = global_vars_dict[key]
+                f.write(f"{key}     {val};\n")
+
+
+def write_geom_prop(case_dir, global_vars_dict, list_gas_spec):
+    globalVar_filename = os.path.join(case_dir, "constant", "globalVars")
+    with open(globalVar_filename, "a+") as f:
+        f.write("//*****************\n")
+        for key in [
+            "HtBcol",
+            "DiaBcol",
+            "LiqHt",
+            "NPS",
+            "NPD",
+            "NPY",
+            "P0",
+            "Pbot",
+            "Pmid",
+            "A_Bcol",
+            "V_flowRate",
+            "rho_gas",
+        ]:
+            if key in global_vars_dict:
+                val = global_vars_dict[key]
+                f.write(f"{key}     {val};\n")
+
+
+def write_globalVars(case_dir, global_vars_dict, list_gas_spec):
+    unknown_keys = check_globalVars(case_dir, global_vars_dict, list_gas_spec)
+    write_unknown_prop(case_dir, global_vars_dict, list_gas_spec, unknown_keys)
+    write_macro_prop(case_dir, global_vars_dict, list_gas_spec)
+    write_liq_prop(case_dir, global_vars_dict, list_gas_spec)
+    write_wc_prop(case_dir, global_vars_dict, list_gas_spec)
+    write_diff_prop(case_dir, global_vars_dict, list_gas_spec)
+    write_henry_prop(case_dir, global_vars_dict, list_gas_spec)
+    write_geom_prop(case_dir, global_vars_dict, list_gas_spec)
+    write_inlet_prop(case_dir, global_vars_dict, list_gas_spec)
+    write_Le_prop(case_dir, global_vars_dict, list_gas_spec)
+    write_cond_prop(case_dir, global_vars_dict, list_gas_spec)
+    write_turb_prop(case_dir, global_vars_dict, list_gas_spec)
+
+
+def fill_global_prop(case_dir: str):
+    global_vars_dict = global_vars_to_dict(
+        os.path.join(case_dir, "constant", "globalVars_temp")
+    )
+    global_vars_dict = add_liquid_properties(case_dir, global_vars_dict)
+    global_vars_dict, list_gas_spec = add_gas_properties(
+        case_dir, global_vars_dict
+    )
+    write_globalVars(case_dir, global_vars_dict, list_gas_spec)
+
+
+if __name__ == "__main__":
+    # fill_global_prop("../../../experimental_cases_new/disengagement/bubble_column_pbe_20L/")
+    fill_global_prop("../../../experimental_cases_new/deckwer17")
diff --git a/bird/utilities/ofio.py b/bird/utilities/ofio.py
index a8cd223f..8263a7ce 100644
--- a/bird/utilities/ofio.py
+++ b/bird/utilities/ofio.py
@@ -1,4 +1,5 @@
 import os
+import re
 import sys
 
 import numpy as np
@@ -157,3 +158,98 @@ def getMeshTime(casePath):
     for entry in files_tmp:
         if entry.startswith("meshFaceCentres"):
             return entry[16:-4]
+
+
+def remove_comments(text):
+    text = re.sub(
+        r"/\*.*?\*/", "", text, flags=re.DOTALL
+    )  # Remove /* */ comments
+    text = re.sub(r"//.*", "", text)  # Remove // comments
+    return text
+
+
+def tokenize(text):
+    # Add spaces around braces and semicolons to make them separate tokens
+    text = re.sub(r"([{}();])", r" \1 ", text)
+    return text.split()
+
+
+def parse_tokens(tokens):
+    def parse_block(index):
+        result = {}
+        while index < len(tokens):
+            token = tokens[index]
+            if token == "}":
+                return result, index + 1
+            elif token == "{":
+                raise SyntaxError("Unexpected '{'")
+            else:
+                key = token
+                index += 1
+                if tokens[index] == "{":
+                    index += 1
+                    value, index = parse_block(index)
+                    result[key] = value
+                elif tokens[index] == "(":
+                    # Parse list
+                    index += 1
+                    lst = []
+                    while tokens[index] != ")":
+                        lst.append(tokens[index])
+                        index += 1
+                    index += 1  # Skip ')'
+                    result[key] = lst
+                    if tokens[index] == ";":
+                        index += 1
+                else:
+                    # Parse scalar value
+                    value = tokens[index]
+                    index += 1
+                    if tokens[index] == ";":
+                        index += 1
+                    result[key] = value
+        return result, index
+
+    parsed, _ = parse_block(0)
+    return parsed
+
+
+def parse_openfoam_dict(text):
+    text = remove_comments(text)
+    tokens = tokenize(text)
+    return parse_tokens(tokens)
+
+
+def read_properties(filename: str):
+    with open(filename, "r+") as f:
+        text = f.read()
+    foam_dict = parse_openfoam_dict(text)
+    return foam_dict
+
+
+def write_openfoam_dict(d, filename, indent=0):
+    lines = []
+
+    indent_str = " " * indent
+
+    for key, value in d.items():
+        if isinstance(value, dict):
+            lines.append(f"{indent_str}{key}")
+            lines.append(f"{indent_str}{{")
+            lines.extend(write_openfoam_dict(value, indent + 4))
+            lines.append(f"{indent_str}}}")
+        elif isinstance(value, list):
+            lines.append(f"{indent_str}{key}")
+            lines.append(f"{indent_str}(")
+            for item in value:
+                lines.append(f"{indent_str}    {item}")
+            lines.append(f"{indent_str});")
+        else:
+            lines.append(f"{indent_str}{key}    {value};")
+
+    with open(filename, "w") as f:
+        lines = write_openfoam_dict(foam_dict)
+        f.write("\n".join(lines))
+        f.write(
+            "\n\n// ************************************************************************* //\n"
+        )