GEOS-ESM / GEOSgcm

GEOS Earth System Model GEOSgcm Fixture
Apache License 2.0
36 stars 13 forks source link

New MBundle_IncrementMod.F90 has unintended consequences and interaction with bundleParser.py #581

Open bena-nasa opened 1 year ago

bena-nasa commented 1 year ago

@sdrabenh @amdasilva @amolod @wmputman @atrayano @mathomp4 @tclune I'm putting this issue here as I don't know where else to put it since it involves source code or input files in at least 3 different repos and I apologize as this might be long and for all the people I tagged but I'm pretty sure every in the list either chimed in on this or is affected. Don't know where to start, this is related to the work Atanas did before he left for PTO.

The problem was that there were some increments, some sort of wet scavenging mumbo jumbo being computed in moist but @wputman removed these in his grand refactoring of moist. Apparently however, @sdrabenh count not release a tag without these increments being output as somebody said so is all I know.

Anyway, these moist increments were something from gocart before and after moist and History needed to write compute 2D column-summed results because that is what our file spec says, but the fields coming from gocart were either 4D (multiple bins) or there were multiple 3D fields from gocart that needed to be aggregated along the bins and vertical, but the problem is that some of these variables need to be aggregated but are not as a single variable with bins! This was previously done directly in moist which is only saw the 3D variables used hard coded logic about what to add together.

Atanas was instructed to make a new routine bundle routine to handle MTRI to aggregate these 3D or 4D fields into the 2D diagnostics needed and output in the MTRI bundle outside of moist which he did and is working (MBundle_IncrementMod.F90). That required adding this to AGCM.rc to tell it to aggregate these separate variables into 1 variable if that was needed which is was for the Nitrates:

# Enable wet scavenging
MTRI_increments::
DU::DU default
SS::SS default
SU::SO4 default
CA.bc::CA.bcphilic default
CA.br::CA.brphilic default
CA.oc::CA.ocphilic default
NI::NO3an1 "NI::NO3an2,NI::NO3an3"
::

and then this was added to History

                               'MTRI%DU::DUIM'             , 'PHYSICS' , 'DUSV' ,
                               'MTRI%SS::SSIM'             , 'PHYSICS' , 'SSSV' ,
                               'MTRI%CA.oc::CA.ocphilicIM' , 'PHYSICS' , 'OCSV' ,
                               'MTRI%CA.bc::CA.bcphilicIM' , 'PHYSICS' , 'BCSV' ,
                               'MTRI%CA.br::CA.brphilicIM' , 'PHYSICS' , 'BRSV' ,
                               'MTRI%SU::SO4IM'            , 'PHYSICS' , 'SUSV' ,
                               'MTRI%NI::NO3an1IM'         , 'PHYSICS' , 'NISV' ,

in the tavg2d_aer collection.

So far so good, in my and Scott's testing it worked, but @pcolarco reported in issue #579 he tried GOCART2G.data with the tag that had this and it failed unsurprisingly since these variables are not export from GOCART2D with those names when GOCART2G is run as the data version, so I looked into it. I said, I'll just remove the MTRI_increment stuff from AGCM.rc but to my surprise, it still fails and not only that I saw in the AGCM.rc in scratch it had this in it how in addition to the above:

MTRI_increments::
DU::DU
SS::SS
CA.oc::CA.ocphilic
CA.bc::CA.bcphilic
CA.br::CA.brphilic
SU::SO4
NI::NO3an1

it took a while but I remembered that there is this bundleParser.py business run in gcm_run.j That looks for anything in History.rc that starts with either TRADVI%, CHEMTRI%,H2ORTRI%, TRI% or MTRI% If it finds those, suppose it finds TRADVI%foo it add this to the AGCM.rc

TRADVI_increment:
foo
::

which then triggers various places in code to go compute the bundle increment. With this new thing now, either bundleParser.py can't add these lines for MTRI as there is no way it can know from History.rc that it needs to add the new thing that was added for the MTRI bundle for nitrates or I have to hard code this in bundleParser.py and have it add the MTRI bundle.

But either way, we can't have this hard coded in AGCM.rc.tmpl and have bundleParser.py add this as it would appear 2x in the AGCM.rc file!

But wait, it gets worse. I was looking in the DAS Histories and saw stuff like this: 'MTRI%PCHEM::OXIM' , 'PHYSICS' , 'DOXDTMST' , in the tavg_3d_odt collection. Now I have to assume this will break the above since now MTRI is handled by MBundle_IncrementMod.F90 which wants to aggregate everything to 2D but this field in MTRI needs to be 3D. So now two different routines are going to do something to the MTRI bundle.

I will check on the consequences of all this.

In any event this is quite a mess. I will do my best to come up with the least kludgy solution to the already kludgy solution.

bena-nasa commented 1 year ago

Note, one option would be to just have physics create a new bundle, not named MTRI that gets passed to the new MBundle_IncrementMod.F90 to store the increments, then the new lines in AGMC.rc.tmpl are kept. Then the other fields the DAS needs that are also in MTRI will not be affected nor will there be any name collision with what is in BundleParser.py with existing things increment calculation enabled via bundleParser.py

Anyone have any opinions? That's what I would recommend as a short term solution until a better solution to this situation can be done in the long run. I did confirm that simply changing the name works...

I've looked at this and can't see a good solution short of this other than a massive overhaul of this bundle increment calculator in chem_shared and hard coded stuff in bundleParser.py as there is simply no way to infer what must be put in these bundles in terms of aggregation from just looking at History.rc

amdasilva commented 1 year ago

For the record, this implementation that specifically names GOCART species is very GOCART specific and violates the abstraction of the aerosol state. I didn't like when it was first implemented, I like it even less now given that this is a refactoring exercise. It is probably OK as a place holder to honor the file spec, but a more general solution is needed. The general principle is that service providers must export the tendencies associated with each service provided, in this case vertical transport and scavenging. However, the variable names in these exports should be derived, abstractly, from the variable names in the bundle. to be serviced. Here an exchange I had with Tom recently:

TOM: Going forward we will of course want to split this into at least two separate calculations. 1st, as we have discussed, any component that provides a service would also provide tendencies related to that service._

ARLINDO: Precisely. Both 3D and vertically averaged – the vertical average is the most common use case, the 3D tendencies are needed only occasionally for debugging.

TOM: Presumably in MOIST this would then be

(delp(t+1) * X(t+1)) – (delp_t X_t)) / (g dt)

rather than a simple raw delta.

ARLINDO: I am not sure if I’d multiply by delp/g here. We must have an option for vertical average. Thinking of it, MAPL should be smart enough to compute vertical averages without any explicit code in the components. We need to design how this would be done but it could have a large payoff in terms of fewer lines of code in component code. Suppose a component export a variable called ‘QV’. The vertical average of ‘QV’, named ‘TQV’, could be added to the export spec with the indication that if is a vertical average of ‘QV’:

     call MAPL_AddExportSpec(GC,                                                 &
                             LONG_NAME             = 'Total Precipitable Water', &
                             UNITS                 = 'kg s-1',                   &
                             SHORT_NAME            = 'TQV',                      &
                             VERTICAL_INTEGRAL_OF  = ‘QV’                        &
                             DIMS                  = MAPL_DimsHorzOnly,          &
                             VLOCATION             = MAPL_VLocationNone,         &
                             __RC__ )

(Would also implement VERTICAL_AVERAGE_OF, VERTICAL_SUM_OF, with INTEGRAL implying mass weighted). In this case, QV must be one of the exports, and it must be automatically allocated if one selects TQV. One may need a call to MAPL at the end of the run method to finish these housekeeping chores. Can you think of something better?

TOM: The component would also choose its own convention for how to decorate the short names to produce short names for these tendencies.

ARLINDO: I don’t think this is necessary. History already provides a device for renaming variables and having a standard convention for denoting tendencies from services would not be a bad thing.

TOM: Then, some other layer would do further reductions on these tendencies to produce the required diagnostics.

ARLINDO: Agreed. The service provider should not know anything about summing over bins or species – this is so component specific. Again, history has an expression parser already, although it requires that only variables appearing in the file can enter the expression. In the very least we need to introduce the notion of “virtual variables” which are “defined” in history.rc but are never written to file – the purpose of virtual variables is to enter the expression parser. A virtual variable could be named with ~ in the beginning: ~DUEXT …

TOM: This could be some other parent GC, or HIST using the expression parser. The current expression parser is not quite powerful enough for this case, but we will improve this.

ARLINDO: Agreed, making the expression parser in history smarter will reduce the number of lines of code in component code. Reusability of code is not a bad thing.

bena-nasa commented 1 year ago

I ended up just making new bundle in physics that these wet scavenging things can be stored in so as to not conflict with MTRI and so Scott can release a tag. That's the only short term solution.