#!/usr/bin/env python3
"""
Test program to compare planeflight vs. History output
"""
import pandas as pd
import xarray as xr

# Read data from Planeflight diagnostic
dframe = pd.read_csv("plane.log.20190701.csv")

# Read data from GEOS-Chem History output
dset = xr.open_mfdataset("OutputDir/GEOSChem.SpeciesConc.20190701*.nc4")

for idx in range(24):

    # Get the indices from the plane.log file
    i_ind = dframe["I-IND"][idx]
    j_ind = dframe["J-IND"][idx]
    p_ind = dframe["P-I"][idx]
    t_ind = dframe["T-IND"][idx]

    # Planeflight output, v/v dry
    tra_077 = dframe["TRA_077"][idx]

    # Planeflight output, molec/cm3
    fura = dframe["FURA"][idx]

    # SpeciesConc output, v/v dry
    vvdry = dset["SpeciesConcVV_FURA"].isel(
        time=t_ind-1,
        lev=p_ind-1,
        lat=j_ind-1,
        lon=i_ind-1,
    ).values

    # SpeciesConc output, molec/cm3
    molec_cm3 = dset["SpeciesConcMND_FURA"].isel(
        time=t_ind-1,
        lev=p_ind-1,
        lat=j_ind-1,
        lon=i_ind-1,
    ).values

    print("-"*70)
    print(f"time, lev, lat, lon indices: {t_ind}, {p_ind}, {j_ind}, {i_ind}")
    print(f"Planefight TRA_077  : {tra_077:>13.6e} v/v dry")
    print(f"SpeciesConcVV_FURA  : {vvdry:>13.6e} v/v dry")
    print(f"Planeflight FURA    : {fura:>13.6f}, molec/cm3")
    print(f"SpeciesConcMND_FURA : {molec_cm3:>13.6f} molec/cm3")