PSURF and TMPU is not defined

[patricksferraz]

Hi Barron,

I'm trying to run aqmbc with the definition gc12_to_ae7.expr. When I used pncglobal2cmaq it was possible to define AIRMOLDEN with default values ​​for TMPU and PSURF, the definition gc_airmolden.expr has this rule but I was unable to identify how I can define PSURF and TMPU.

I tried adding the content below in the first two lines of the gc_airmolden.expr file:

PSURF = (np.ones_like(O3[:,0][:, None])*1013.25)
TMPU = (np.ones_like(O3[:])*np.array([287.7,286.9,286.0,285.2,284.3,283.4,282.6,281.7,280.7,279.8,278.9,277.9,276.8,275.3,273.6,271.8,270.0,268.2,265.8,262.8,259.7,256.4,252.9,249.2,245.2,241.0,236.4,230.5,223.6,216.8,216.6,216.6,216.6,216.6,216.6,216.6,216.6,217.5,219.8,222.0,225.3,233.6,250.5,269.2,260.0,237.7,214.3])[None, :47, None, None])

but there is a problem with the shape:

Traceback (most recent call last):
  File "bcon.py", line 431, in <module>
    history=str(args))
  File "bcon.py", line 196, in bc
    outf = func(outf)
  File "bcon.py", line 124, in translate
    outf = infile.eval(exprstr, inplace=False)
  File "/home/patrick.ferraz/.conda/envs/pseudonetcdf/lib/python3.6/site-packages/PseudoNetCDF/core/_files.py", line 995, in eval
    exec(comp, None, vardict)
  File "none", line 3, in <module>
ValueError: operands could not be broadcast together with shapes (24,35,164) (1,47,1,1)
make[1]: *** [ts20150117.bpch.BCON.nc] Error 1
make[1]: Leaving directory `/scratch/patrick.ferraz/numerical_models/aqmbc_2/BCON'
make: *** [bcon] Error 2

[barronh]

It looks like you only saved out 35 layers. When you do this, the vertical coordinate can only have 35 layers. By default, the lines you used have 47. Changing the slice :47 to :35 should help. See below.

PSURF = (np.ones_like(O3[:,0][:, None])*1013.25)
TMPU = (np.ones_like(O3[:])*np.array([287.7,286.9,286.0,285.2,284.3,283.4,282.6,281.7,280.7,279.8,278.9,277.9,276.8,275.3,273.6,271.8,270.0,268.2,265.8,262.8,259.7,256.4,252.9,249.2,245.2,241.0,236.4,230.5,223.6,216.8,216.6,216.6,216.6,216.6,216.6,216.6,216.6,217.5,219.8,222.0,225.3,233.6,250.5,269.2,260.0,237.7,214.3])[None, :35, None, None])

Please confirm if this works for you or ask if you need further help.

[patricksferraz]

Barron, did not work, presented the following error:

Traceback (most recent call last):
  File "bcon.py", line 431, in <module>
    history=str(args))
  File "bcon.py", line 196, in bc
    outf = func(outf)
  File "bcon.py", line 124, in translate
    outf = infile.eval(exprstr, inplace=False)
  File "/home/patrick.ferraz/.conda/envs/pseudonetcdf/lib/python3.6/site-packages/PseudoNetCDF/core/_files.py", line 995, in eval
    exec(comp, None, vardict)
  File "none", line 3, in <module>
ValueError: operands could not be broadcast together with shapes (24,35,164) (1,35,1,1)
make[1]: *** [ts20150117.bpch.BCON.nc] Error 1
make[1]: Leaving directory `/scratch/patrick.ferraz/numerical_models/aqmbc_2/BCON'
make: *** [bcon] Error 2

I enabled ND49 considering 47 layers, as shown below:

IMIN, IMAX of region    :   1  72
JMIN, JMAX of region    :   1  46
LMIN, LMAX of region    :   1  47

My GeosChem output files have 4 dimensions for all tracers, eg:

dimensions:
        longitude = 72 ;
        latitude = 46 ;
        layer = 47 ;
        nv = 2 ;
        tnv = 2 ;
        time = 24 ;
        layer47 = 47 ;
        layer_bounds = 48 ;
[...]
float IJ-AVG-$_O3(time, layer47, latitude, longitude);
                IJ-AVG-$_O3:category = "IJ-AVG-$" ;
                IJ-AVG-$_O3:tracerid = 2 ;
                IJ-AVG-$_O3:base_units = b"                                        " ;
                IJ-AVG-$_O3:catoffset = 0 ;
                IJ-AVG-$_O3:noscale = False ;
                IJ-AVG-$_O3:STARTK = 0 ;
                IJ-AVG-$_O3:STARTJ = 0 ;
                IJ-AVG-$_O3:STARTI = 0 ;
                IJ-AVG-$_O3:cattracerid = 2 ;
                IJ-AVG-$_O3:shortname = "O3" ;
                IJ-AVG-$_O3:fullname = "O3 tracer" ;
                IJ-AVG-$_O3:kgpermole = 0.048 ;
                IJ-AVG-$_O3:carbon = 1 ;
                IJ-AVG-$_O3:scale = 1000000000.0 ;
                IJ-AVG-$_O3:units = "ppbv" ;
[...]

I believe that the error is occurring due to considering only three dimensions (24,35,164) but that they also differ from the defined values. I performed the test with the definition gc12_to_cb6r3.expr and there is no error.

[barronh]

I forgot that the expression is being applied after the ij extraction and vertical interpolation. So, the code will work with the update below.

PSURF = (np.ones_like(O3[:,0][:, None])*1013.25)
TMPU = (np.ones_like(O3[:])*np.array([287.7,286.9,286.0,285.2,284.3,283.4,282.6,281.7,280.7,279.8,278.9,277.9,276.8,275.3,273.6,271.8,270.0,268.2,265.8,262.8,259.7,256.4,252.9,249.2,245.2,241.0,236.4,230.5,223.6,216.8,216.6,216.6,216.6,216.6,216.6,216.6,216.6,217.5,219.8,222.0,225.3,233.6,250.5,269.2,260.0,237.7,214.3])[None, :35, None])

Note, however, that the temperature is a gross approximation. Ideally, the ND49 would have been configured to output TMPU and PSURF so that you don't need these approximations.

[patricksferraz]

I had made an attempt using the informed configuration ([None,: 35, None]) and a new error occurred:

Traceback (most recent call last):
  File "bcon.py", line 431, in <module>
    history=str(args))
  File "bcon.py", line 213, in bc
    return saveioapi(wndwf, outf, outpath, metaf, dimkeys)
  File "bcon.py", line 308, in saveioapi
    out = outf.save(outpath, format=outformat, verbose=0, outmode='w')
  File "/home/patrick.ferraz/.conda/envs/pseudonetcdf/lib/python3.6/site-packages/PseudoNetCDF/core/_files.py", line 2325, in save
    return pncwrite(self, *args, **kwds)
  File "/home/patrick.ferraz/.conda/envs/pseudonetcdf/lib/python3.6/site-packages/PseudoNetCDF/_getwriter.py", line 42, in pncwrite
    return pncgen(*args, **kwds)
  File "/home/patrick.ferraz/.conda/envs/pseudonetcdf/lib/python3.6/site-packages/PseudoNetCDF/pncgen.py", line 245, in pncgen
    format=format)
  File "/home/patrick.ferraz/.conda/envs/pseudonetcdf/lib/python3.6/site-packages/PseudoNetCDF/pncgen.py", line 52, in convert
    self.addVariables(pfile, nfile)
  File "/home/patrick.ferraz/.conda/envs/pseudonetcdf/lib/python3.6/site-packages/PseudoNetCDF/pncgen.py", line 174, in addVariables
    self.addVariableData(pfile, nfile, k)
  File "/home/patrick.ferraz/.conda/envs/pseudonetcdf/lib/python3.6/site-packages/PseudoNetCDF/pncgen.py", line 159, in addVariableData
    nvar[:] = pvar[...]
  File "netCDF4/_netCDF4.pyx", line 4893, in netCDF4._netCDF4.Variable.__setitem__
  File "netCDF4/_netCDF4.pyx", line 5151, in netCDF4._netCDF4.Variable._put
IndexError: size of data array does not conform to slice
make[1]: *** [ts20150117.bpch.BCON.nc] Error 1
make[1]: Leaving directory `/scratch/patrick.ferraz/numerical_models/aqmbc_2/BCON'
make: *** [bcon] Error 2

I will follow your recommendation and configure ND49 to generate PSURF and TMPU

[barronh]

patrick - Can you post an example input data file and configuration files somewhere that I can download and test them?

[patricksferraz]

Follow the link: GDrive

I am using GeosChem 12.3.2.

At the root of the folder contains the files referring to GeosChem and in the aqmbc folder, the settings used to execute it.

[barronh]

Four part response:

1. Fix for simplistic TMPU/PRES
2. Updates for combine which is used in some visualizations
3. Description of how to add TMPU/PSURF to ND49 output
4. Request for confirmation

1. Fix for simplistic TMPU/PRES

Ah... The problem is that the BCON is 3D (TSTEP, LAY, PERIM) and the ICON is 4D (TSTEP, LAY, ROW, COL). Thus, the layer has to dynamically reshape... I had only tested one or the other with that quick fix. The code below works with your files.

# Source: GEOS-Chem v8,v9,v10,v11,v12 bpch
newshape = O3.shape
reshape = [1] * O3.ndim
reshape[1] = newshape[1]
PRES = O3[:] * 0 + np.broadcast_to(
    101325. * hybm[:].reshape(reshape) + hyam[:].reshape(reshape),
    newshape
)
PRES.units = 'Pa'
TMPU = (O3[:] * 0 +
    np.array([287.7,286.9, 286.0, 285.2, 284.3, 283.4, 282.6, 281.7, 280.7, 279.8, 278.9, 277.9, 276.8, 275.3, 273.6, 271.8, 270.0, 268.2, 265.8, 262.8, 259.7, 256.4, 252.9, 249.2, 245.2, 241.0, 236.4, 230.5, 223.6, 216.8, 216.6, 216.6, 216.6, 216.6, 216.6, 216.6, 216.6, 217.5, 219.8, 222.0, 225.3, 233.6, 250.5, 269.2, 260.0, 237.7, 214.3])[:35].reshape(reshape)
)
TMPU.units = 'K'
AIRMOLDEN = PRES / 8.3144621 / TMPU
AIRMOLDEN.units ="mole/m3"
del reshape, newshape

2. Updates for combine which is used in some visualizations

Note that BCON/ICON will be created and combine will fail. combine requires you to specify definitions that you want to make quick check files. The GEOS-Chem files don't have all the species in cb6. You would have to update it below. Notice that I removed O3 and NOx because you did not use the gas-phase definitions.

PMIJ = AALJ + ACAJ + ACLJ + AECI + AECJ + AFEJ + AGLYJ + AISO2J + AISO3J + AKJ + AMGJ + AMNJ + ANAJ + ANH4I + ANH4J + ANO3I + ANO3J + AOTHRJ + APNCOMI + APNCOMJ + APOCI + APOCJ + ASIJ + ASO4I + ASO4J + ASQTJ + ATIJ
PMIJ.long_name = 'PMIJ'
PMK = ACLK + ANO3K + ASO4K + ASOIL
PMK.long_name = 'PMK'
ASO4IJ = ASO4I + ASO4J
ASO4IJ.long_name = 'ASO4IJ'
ANO3IJ = ANO3I + ANO3J
ANO3IJ.long_name = 'ANO3IJ'
ANAIJ = ANAJ
ANAIJ.long_name = 'ANAIJ'

3. Description of how to add TMPU/PSURF to ND49 output

To include TMPU and PSURF, you need to configure the input.geos ND49 section to include the associated numbers for pressure and temperature.[1] The associated numbers change by version, but the most recent numbers are in the ND49 list[2] as 512 and 516.

Remember that the goal is to get AIRMOLDEN to convert mixing ratios to concentrations for aerosols. So, you could also use Air density (currently 507), which would have to be converted from molec/cm3 to mole/m3.

4. Request for confirmation

Please confirm that this works for you. I was able to run this with your file on colab.research.google.com.

References:

[1] http://wiki.seas.harvard.edu/geos-chem/index.php/The_input.geos_file#ND49_diagnostic [2] http://wiki.seas.harvard.edu/geos-chem/index.php/List_of_diagnostics_archived_to_bpch_format#ND49:_Instantaneous_timeseries

[patricksferraz]

Thank you very much Barron, it's running now

[patricksferraz]

Hello again Barron, all right?

I hope you can help us. We were able to finish the execution of GeosChem by enabling TMPU and PSURF, we ran aqmbc and everything went without errors. However, when we used the outputs (BCON and ICON) in CMAQ, the following error occurred:

     *** ERROR ABORT in subroutine AQCHEM on PE 000          
     Maximum AQCHEM total iterations exceeded
 PM3EXIT:  DTBUF 1:00:00   Jan. 1, 2015  
     Date and time 1:00:00   Jan. 1, 2015   (2015001:010000)

We believe it is related to scale. When we used the mappings from pncglobal2cmaq, although it was for cb6, it did not contain a scale change:

        [...]
        "AALJ": {
            "expression": "0.05695 * DST1",
            "outunit": "micrograms/m**3"
        },
        "AECI": {
            "expression": "0.1 * BCPI + 0.1 * BCPO",
            "outunit": "micrograms/m**3"
        },
        "AECJ": {
            "expression": "0.9 * BCPI + 0.9 * BCPO",
            "outunit": "micrograms/m**3"
        },
        "APOCI": {
            "expression": "0.1 * OCPI + 0.1 * OCPO",
            "outunit": "micrograms/m**3"
        },
        "APOCJ": {
            "expression": "0.9 * OCPI + 0.9 * OCPO + 0.01075 * DST1",
            "outunit": "micrograms/m**3"
        },
        [...]

In aqmbc, in the ae6 mapping, there is also no change:

ACAJ = (0.0118 * SALA[:] * 0.0314 + 0.0794 * DST1[:] * 0.029) * AIRMOLDEN
ACAJ.units ="micrograms/m**3"
ACLJ = (0.00945 * DST1[:] * 0.029 + 0.5538 * SALA[:] * 0.0314) * AIRMOLDEN
ACLJ.units ="micrograms/m**3"
ACLK = (0.0119 * (DST2[:] * 0.029 + DST3[:] * 0.029 + DST4[:] * 0.029) + 0.5538 * SALC[:] * 0.0314) * AIRMOLDEN
ACLK.units ="micrograms/m**3"
AECI = (0.001 * (BCPI[:] * 0.012 + BCPO[:] * 0.012)) * AIRMOLDEN
AECI.units ="micrograms/m**3"
AECJ = (0.999 * (BCPI[:] * 0.012 + BCPO[:] * 0.012)) * AIRMOLDEN
AECJ.units ="micrograms/m**3"

While in ae7 the change is made in the multiplication by 1e9, could you explain the reason?

AALJ = (0.05695 * SpeciesConc_DST1[:] * 0.029) * AIRMOLDEN * 1e9
AALJ.units ="micrograms/m**3"
ACAJ = (0.0118 * SpeciesConc_SALA[:] * 0.0314 + 0.0794 * SpeciesConc_DST1[:] * 0.029) * AIRMOLDEN * 1e9
ACAJ.units ="micrograms/m**3"
ACLJ = (0.00945 * SpeciesConc_DST1[:] * 0.029 + 0.5538 * SpeciesConc_SALA[:] * 0.0314) * AIRMOLDEN * 1e9
ACLJ.units ="micrograms/m**3"
ACLK = (0.0119 * (SpeciesConc_DST2[:] * 0.029 + SpeciesConc_DST3[:] * 0.029 + SpeciesConc_DST4[:] * 0.029) + 0.5538 * SpeciesConc_SALC[:] * 0.0314) * AIRMOLDEN * 1e9
ACLK.units ="micrograms/m**3"
AECI = (0.001 * (SpeciesConc_BCPI[:] * 0.012 + SpeciesConc_BCPO[:] * 0.012)) * AIRMOLDEN * 1e9
AECI.units ="micrograms/m**3"
AECJ = (0.999 * (SpeciesConc_BCPI[:] * 0.012 + SpeciesConc_BCPO[:] * 0.012)) * AIRMOLDEN * 1e9
AECJ.units ="micrograms/m**3"

[barronh]

I’m on vacation, so will only be able to reply briefly for now. The conversions are taking units of mol/mol from GEOS-Chem and converting to micrograms/m**3. This molecular weights are in kg/mol, so 1e-9 converts from kg to micrograms.

[patricksferraz]

Thanks Barron.