Using data from GEOS-FP

This example demonstrates how to use the GEOS-FP data loader in the EarthSciML ecosystem. The GEOS-FP data loader is used to load data from the GEOS-FP dataset.

First, let's initialize some packages and set up the GEOS-FP equation system.

using EarthSciData, EarthSciMLBase
using DomainSets, ModelingToolkit, MethodOfLines, DifferentialEquations
using Dates, Plots, DataFrames

# Set up system
@parameters t lev lon lat
geosfp = GEOSFP("4x5", t)

\[ \begin{align} A3dyn_{+}DTRAIN\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}DTRAIN_{interp} \right), t, lon, lat, lev \right) \\ A3dyn_{+}OMEGA\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}OMEGA_{interp} \right), t, lon, lat, lev \right) \\ A3dyn_{+}RH\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}RH_{interp} \right), t, lon, lat, lev \right) \\ A3dyn_{+}U\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}U_{interp} \right), t, lon, lat, lev \right) \\ A3dyn_{+}V\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}V_{interp} \right), t, lon, lat, lev \right) \\ A3cld_{+}CLOUD\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}CLOUD_{interp} \right), t, lon, lat, lev \right) \\ A3cld_{+}OPTDEPTH\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}OPTDEPTH_{interp} \right), t, lon, lat, lev \right) \\ A3cld_{+}QCCU\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}QCCU_{interp} \right), t, lon, lat, lev \right) \\ A3cld_{+}QI\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}QI_{interp} \right), t, lon, lat, lev \right) \\ A3cld_{+}QL\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}QL_{interp} \right), t, lon, lat, lev \right) \\ A3cld_{+}TAUCLI\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}TAUCLI_{interp} \right), t, lon, lat, lev \right) \\ A3cld_{+}TAUCLW\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}TAUCLW_{interp} \right), t, lon, lat, lev \right) \\ I3_{+}PS\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PS_{interp} \right), t, lon, lat \right) \\ I3_{+}PV\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PV_{interp} \right), t, lon, lat, lev \right) \\ I3_{+}QV\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}QV_{interp} \right), t, lon, lat, lev \right) \\ I3_{+}T\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}T_{interp} \right), t, lon, lat, lev \right) \\ A3mstE_{+}CMFMC\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}CMFMC_{interp} \right), t, lon, lat, lev \right) \\ A3mstE_{+}PFICU\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PFICU_{interp} \right), t, lon, lat, lev \right) \\ A3mstE_{+}PFILSAN\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PFILSAN_{interp} \right), t, lon, lat, lev \right) \\ A3mstE_{+}PFLCU\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PFLCU_{interp} \right), t, lon, lat, lev \right) \\ A3mstE_{+}PFLLSAN\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PFLLSAN_{interp} \right), t, lon, lat, lev \right) \\ A1_{+}ALBEDO\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}ALBEDO_{interp} \right), t, lon, lat \right) \\ A1_{+}CLDTOT\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}CLDTOT_{interp} \right), t, lon, lat \right) \\ A1_{+}EFLUX\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}EFLUX_{interp} \right), t, lon, lat \right) \\ A1_{+}EVAP\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}EVAP_{interp} \right), t, lon, lat \right) \\ A1_{+}FRSEAICE\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}FRSEAICE_{interp} \right), t, lon, lat \right) \\ A1_{+}FRSNO\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}FRSNO_{interp} \right), t, lon, lat \right) \\ A1_{+}GRN\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}GRN_{interp} \right), t, lon, lat \right) \\ A1_{+}GWETROOT\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}GWETROOT_{interp} \right), t, lon, lat \right) \\ A1_{+}GWETTOP\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}GWETTOP_{interp} \right), t, lon, lat \right) \\ A1_{+}HFLUX\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}HFLUX_{interp} \right), t, lon, lat \right) \\ A1_{+}LAI\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}LAI_{interp} \right), t, lon, lat \right) \\ A1_{+}LWI\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}LWI_{interp} \right), t, lon, lat \right) \\ A1_{+}LWGNT\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}LWGNT_{interp} \right), t, lon, lat \right) \\ A1_{+}LWTUP\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}LWTUP_{interp} \right), t, lon, lat \right) \\ A1_{+}PARDF\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PARDF_{interp} \right), t, lon, lat \right) \\ A1_{+}PARDR\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PARDR_{interp} \right), t, lon, lat \right) \\ A1_{+}PBLH\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PBLH_{interp} \right), t, lon, lat \right) \\ A1_{+}PRECANV\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PRECANV_{interp} \right), t, lon, lat \right) \\ A1_{+}PRECCON\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PRECCON_{interp} \right), t, lon, lat \right) \\ A1_{+}PRECLSC\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PRECLSC_{interp} \right), t, lon, lat \right) \\ A1_{+}PRECSNO\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PRECSNO_{interp} \right), t, lon, lat \right) \\ A1_{+}PRECTOT\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PRECTOT_{interp} \right), t, lon, lat \right) \\ A1_{+}QV2M\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}QV2M_{interp} \right), t, lon, lat \right) \\ A1_{+}SEAICE00\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE00_{interp} \right), t, lon, lat \right) \\ A1_{+}SEAICE10\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE10_{interp} \right), t, lon, lat \right) \\ A1_{+}SEAICE20\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE20_{interp} \right), t, lon, lat \right) \\ A1_{+}SEAICE30\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE30_{interp} \right), t, lon, lat \right) \\ A1_{+}SEAICE40\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE40_{interp} \right), t, lon, lat \right) \\ A1_{+}SEAICE50\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE50_{interp} \right), t, lon, lat \right) \\ A1_{+}SEAICE60\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE60_{interp} \right), t, lon, lat \right) \\ A1_{+}SEAICE70\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE70_{interp} \right), t, lon, lat \right) \\ A1_{+}SEAICE80\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE80_{interp} \right), t, lon, lat \right) \\ A1_{+}SEAICE90\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE90_{interp} \right), t, lon, lat \right) \\ A1_{+}SLP\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SLP_{interp} \right), t, lon, lat \right) \\ A1_{+}SNODP\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SNODP_{interp} \right), t, lon, lat \right) \\ A1_{+}SNOMAS\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SNOMAS_{interp} \right), t, lon, lat \right) \\ A1_{+}SWGDN\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SWGDN_{interp} \right), t, lon, lat \right) \\ A1_{+}TO3\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}TO3_{interp} \right), t, lon, lat \right) \\ A1_{+}TROPPT\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}TROPPT_{interp} \right), t, lon, lat \right) \\ A1_{+}TS\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}TS_{interp} \right), t, lon, lat \right) \\ A1_{+}T2M\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}T2M_{interp} \right), t, lon, lat \right) \\ A1_{+}U10M\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}U10M_{interp} \right), t, lon, lat \right) \\ A1_{+}USTAR\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}USTAR_{interp} \right), t, lon, lat \right) \\ A1_{+}V10M\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}V10M_{interp} \right), t, lon, lat \right) \\ A1_{+}Z0M\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}Z0M_{interp} \right), t, lon, lat \right) \\ A1_{+}T10M\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}T10M_{interp} \right), t, lon, lat \right) \\ A1_{+}Q850\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}Q850_{interp} \right), t, lon, lat \right) \\ A3mstC_{+}DQRCU\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}DQRCU_{interp} \right), t, lon, lat, lev \right) \\ A3mstC_{+}DQRLSAN\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}DQRLSAN_{interp} \right), t, lon, lat, lev \right) \\ A3mstC_{+}REEVAPCN\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}REEVAPCN_{interp} \right), t, lon, lat, lev \right) \\ A3mstC_{+}REEVAPLS\left( t \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}REEVAPLS_{interp} \right), t, lon, lat, lev \right) \\ P\left( t \right) =& P_{unit} \mathrm{DataInterpolations.LinearInterpolation{Vector{Float64}, UnitRange{Int64}, Float64}([0.0, 0.04804826, 6.593752, 13.1348, 19.61311, 26.09201, 32.57081, 38.98201, 45.33901, 51.69611, 58.05321, 64.36264, 70.62198, 78.83422, 89.09992, 99.36521, 109.1817, 118.9586, 128.6959, 142.91, 156.26, 169.609, 181.619, 193.097, 203.259, 212.15, 218.776, 223.898, 224.363, 216.865, 201.192, 176.93, 150.393, 127.837, 108.663, 92.36572, 78.51231, 66.60341, 56.38791, 47.64391, 40.17541, 33.81001, 28.36781, 23.73041, 19.7916, 16.4571, 13.6434, 11.2769, 9.292942, 7.619842, 6.216801, 5.046801, 4.076571, 3.276431, 2.620211, 2.08497, 1.65079, 1.30051, 1.01944, 0.7951341, 0.6167791, 0.4758061, 0.3650411, 0.2785261, 0.211349, 0.159495, 0.119703, 0.08934502, 0.06600001, 0.04758501, 0.0327, 0.02, 0.01], 1:73, false)}\left( lev \right) + I3_{+}PS\left( t \right) \mathrm{DataInterpolations.LinearInterpolation{Vector{Float64}, UnitRange{Int64}, Float64}([1.0, 0.984952, 0.963406, 0.941865, 0.920387, 0.898908, 0.877429, 0.856018, 0.8346609, 0.8133039, 0.7919469, 0.7706375, 0.7493782, 0.721166, 0.6858999, 0.6506349, 0.6158184, 0.5810415, 0.5463042, 0.4945902, 0.4437402, 0.3928911, 0.3433811, 0.2944031, 0.2467411, 0.2003501, 0.1562241, 0.1136021, 0.06372006, 0.02801004, 0.006960025, 8.175413e-9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 1:73, false)}\left( lev \right) \end{align} \]

We can see above the different variables that are available in the GEOS-FP dataset. But also, here they are in table form:

vars = states(geosfp)
DataFrame(
        :Name => [string(Symbolics.tosymbol(v, escape=false)) for v ∈ vars],
        :Units => [ModelingToolkit.get_unit(v) for v ∈ vars],
        :Description => [ModelingToolkit.getdescription(v) for v ∈ vars],
)

The GEOS-FP equation system isn't an ordinary differential equation (ODE) system, so we can't run it by itself. To fix this, we create another equation system that is an ODE. (We don't actually end up using this system for anything, it's just necessary to get the system to compile.)

function Example(t)
    @variables c(t) = 5.0
    D = Differential(t)
    ODESystem([D(c) ~ sin(lat * π / 180.0 * 6) + sin(lon * π / 180 * 6)], t, name=:Docs₊Example)
end
examplesys = Example(t)

\[ \begin{align} \frac{\mathrm{d} c\left( t \right)}{\mathrm{d}t} =& \sin\left( 0.10472 lon \right) + \sin\left( 0.10472 lat \right) \end{align} \]

Now, let's couple these two systems together, and also add in advection and some information about the domain:

domain = DomainInfo(
    partialderivatives_δxyδlonlat,
    constIC(0.0, t ∈ Interval(Dates.datetime2unix(DateTime(2022, 1, 1)), Dates.datetime2unix(DateTime(2022, 1, 3)))),
    zerogradBC(lat ∈ Interval(-80.0f0, 80.0f0)),
    periodicBC(lon ∈ Interval(-180.0f0, 180.0f0)),
    zerogradBC(lev ∈ Interval(1.0f0, 11.0f0)),
)

composed_sys = couple(examplesys, domain, geosfp)
pde_sys = get_mtk(composed_sys)

\[ \begin{align} \frac{\mathrm{d}}{\mathrm{d}t} Docs_{+}Example_{+}c\left( t, lat, lon, lev \right) =& \sin\left( 0.10472 lon \right) + \sin\left( 0.10472 lat \right) \\ EarthSciData_{+}GEOSFP_{+}A3dyn_{+}DTRAIN\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}DTRAIN_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3dyn_{+}OMEGA\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}OMEGA_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3dyn_{+}RH\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}RH_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3dyn_{+}U\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}U_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3dyn_{+}V\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}V_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3cld_{+}CLOUD\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}CLOUD_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3cld_{+}OPTDEPTH\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}OPTDEPTH_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3cld_{+}QCCU\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}QCCU_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3cld_{+}QI\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}QI_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3cld_{+}QL\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}QL_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3cld_{+}TAUCLI\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}TAUCLI_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3cld_{+}TAUCLW\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}TAUCLW_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}I3_{+}PS\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PS_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}I3_{+}PV\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PV_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}I3_{+}QV\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}QV_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}I3_{+}T\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}T_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3mstE_{+}CMFMC\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}CMFMC_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3mstE_{+}PFICU\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PFICU_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3mstE_{+}PFILSAN\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PFILSAN_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3mstE_{+}PFLCU\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PFLCU_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3mstE_{+}PFLLSAN\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PFLLSAN_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}ALBEDO\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}ALBEDO_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}CLDTOT\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}CLDTOT_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}EFLUX\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}EFLUX_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}EVAP\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}EVAP_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}FRSEAICE\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}FRSEAICE_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}FRSNO\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}FRSNO_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}GRN\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}GRN_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}GWETROOT\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}GWETROOT_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}GWETTOP\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}GWETTOP_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}HFLUX\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}HFLUX_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}LAI\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}LAI_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}LWI\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}LWI_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}LWGNT\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}LWGNT_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}LWTUP\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}LWTUP_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}PARDF\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PARDF_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}PARDR\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PARDR_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}PBLH\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PBLH_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}PRECANV\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PRECANV_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}PRECCON\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PRECCON_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}PRECLSC\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PRECLSC_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}PRECSNO\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PRECSNO_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}PRECTOT\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}PRECTOT_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}QV2M\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}QV2M_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}SEAICE00\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE00_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}SEAICE10\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE10_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}SEAICE20\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE20_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}SEAICE30\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE30_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}SEAICE40\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE40_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}SEAICE50\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE50_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}SEAICE60\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE60_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}SEAICE70\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE70_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}SEAICE80\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE80_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}SEAICE90\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SEAICE90_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}SLP\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SLP_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}SNODP\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SNODP_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}SNOMAS\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SNOMAS_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}SWGDN\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}SWGDN_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}TO3\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}TO3_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}TROPPT\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}TROPPT_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}TS\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}TS_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}T2M\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}T2M_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}U10M\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}U10M_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}USTAR\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}USTAR_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}V10M\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}V10M_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}Z0M\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}Z0M_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}T10M\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}T10M_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A1_{+}Q850\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}Q850_{interp} \right), t, lon, lat \right) \\ EarthSciData_{+}GEOSFP_{+}A3mstC_{+}DQRCU\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}DQRCU_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3mstC_{+}DQRLSAN\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}DQRLSAN_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3mstC_{+}REEVAPCN\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}REEVAPCN_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}A3mstC_{+}REEVAPLS\left( t, lat, lon, lev \right) =& \mathrm{interp!}\left( \mathrm{EarthSciData}\left( GEOSFPFileSet_{x}REEVAPLS_{interp} \right), t, lon, lat, lev \right) \\ EarthSciData_{+}GEOSFP_{+}P\left( t, lat, lon, lev \right) =& EarthSciData_{+}GEOSFP_{+}P_{unit} \mathrm{DataInterpolations.LinearInterpolation{Vector{Float64}, UnitRange{Int64}, Float64}([0.0, 0.04804826, 6.593752, 13.1348, 19.61311, 26.09201, 32.57081, 38.98201, 45.33901, 51.69611, 58.05321, 64.36264, 70.62198, 78.83422, 89.09992, 99.36521, 109.1817, 118.9586, 128.6959, 142.91, 156.26, 169.609, 181.619, 193.097, 203.259, 212.15, 218.776, 223.898, 224.363, 216.865, 201.192, 176.93, 150.393, 127.837, 108.663, 92.36572, 78.51231, 66.60341, 56.38791, 47.64391, 40.17541, 33.81001, 28.36781, 23.73041, 19.7916, 16.4571, 13.6434, 11.2769, 9.292942, 7.619842, 6.216801, 5.046801, 4.076571, 3.276431, 2.620211, 2.08497, 1.65079, 1.30051, 1.01944, 0.7951341, 0.6167791, 0.4758061, 0.3650411, 0.2785261, 0.211349, 0.159495, 0.119703, 0.08934502, 0.06600001, 0.04758501, 0.0327, 0.02, 0.01], 1:73, false)}\left( lev \right) + EarthSciData_{+}GEOSFP_{+}I3_{+}PS\left( t, lat, lon, lev \right) \mathrm{DataInterpolations.LinearInterpolation{Vector{Float64}, UnitRange{Int64}, Float64}([1.0, 0.984952, 0.963406, 0.941865, 0.920387, 0.898908, 0.877429, 0.856018, 0.8346609, 0.8133039, 0.7919469, 0.7706375, 0.7493782, 0.721166, 0.6858999, 0.6506349, 0.6158184, 0.5810415, 0.5463042, 0.4945902, 0.4437402, 0.3928911, 0.3433811, 0.2944031, 0.2467411, 0.2003501, 0.1562241, 0.1136021, 0.06372006, 0.02801004, 0.006960025, 8.175413e-9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 1:73, false)}\left( lev \right) \end{align} \]

Now, finally, we can run the simulation and plot the GEOS-FP wind fields in the result:

(The code below is commented out because it is very slow right now. A faster solution is coming soon!)

# discretization = MOLFiniteDifference([lat => 10, lon => 10, lev => 10], t, approx_order=2)
# @time pdeprob = discretize(pde_sys, discretization)

# pdesol = solve(pdeprob, Tsit5(), saveat=3600.0)

# discrete_lon = pdesol[lon]
# discrete_lat = pdesol[lat]
# discrete_lev = pdesol[lev]
# discrete_t = pdesol[t]

# @variables meanwind₊u(..) meanwind₊v(..) examplesys₊c(..)
# sol_u = pdesol[meanwind₊u(t, lat, lon, lev)]
# sol_v = pdesol[meanwind₊v(t, lat, lon, lev)]
# sol_c = pdesol[examplesys₊c(t, lat, lon, lev)]

# anim = @animate for k in 1:length(discrete_t)
#     p1 = heatmap(discrete_lon, discrete_lat, sol_c[k, 1:end, 1:end, 2], clim=(minimum(sol_c[:, :, :, 2]), maximum(sol_c[:, :, :, 2])),
#             xlabel="Longitude", ylabel="Latitude", title="examplesys.c: $(Dates.unix2datetime(discrete_t[k]))")
#     p2 = heatmap(discrete_lon, discrete_lat, sol_u[k, 1:end, 1:end, 2], clim=(minimum(sol_u[:, :, :, 2]), maximum(sol_u[:, :, :, 2])), 
#             title="U")
#     p3 = heatmap(discrete_lon, discrete_lat, sol_v[k, 1:end, 1:end, 2], clim=(minimum(sol_v[:, :, :, 2]), maximum(sol_v[:, :, :, 2])),
#             title="V")
#     plot(p1, p2, p3, size=(800, 500))
# end
# gif(anim, "animation.gif", fps = 8)