[PULL REQUEST] HEMCO grid feature and interface restructuring

[jimmielin]

This is the successor to #364, which has been reposted under a new branch due to the rebase to a newer version of GEOS-Chem (13.0.0 release on the main branch)

Please find attached the original description of the feature below:

Hi GCST,

This update includes the HEMCO 3.0 planned feature, the "HEMCO Intermediate Grid", which allows HEMCO to operate on a different grid than the GC-Classic simulation, at runtime. The update is rather hefty and here goes a quick summary of what's happening under the hood:

Changes summary

• The intermediate grid feature is controlled by three new input options, Input_Opt%IMGRID_XSCALE, Input_Opt%IMGRID_YSCALE, which defines the refinement scale for the HEMCO grid (setting them to 2, 2 means that HEMCO works on 1x1.25 on a 2x2.5 simulation). This number is hardcoded as 1, 1 (= do nothing) under hco_interface_gc_mod.F90 right now because the feature requires testing and validation and its best kept hidden for now. If any of the scales are not equal to 1, then Input_Opt%LIMGRID is enabled, which enables the intermediate grid feature. In the future these should be switches in input.geos but they are not implemented yet as the feature has not passed scientific validation.

• The HEMCO grid is described in State_Grid_HCO stored in hco_state_gc_mod.F90. It reuses the GrdState nicely. Generation of a "refined"/scaled grid is using the new subroutine gc_grid_mod.F90 : Compute_Scaled_Grid.

• The main idea: Data is stored in HEMCO resolution in HcoState as usual, and only regridded on-demand. The regridded data lies in a single buffer (actually, two.) (TMP_MDL in hco_utilities_gc_mod.F90) which is updated when a new field is requested via the new HCO_GC_* subroutines below...

• Everything intermediate-grid related is wrapped in MODEL_CLASSIC and IF ( Input_Opt%LIMGRID ) checks. There are minor changes to the code that affect all simulations, though:

  • All calls to HCO_GetPtr should now be HCO_GC_GetPtr, implemented in hco_utilities_gc_mod.F90.
  • All calls to HCO_EvalFld should now be HCO_GC_EvalFld.
  • All calls to GetHcoDiagn should now be HCO_GC_GetDiagn.
  • The above calls redirect to their in-HEMCO counterparts when intermediate grid is disabled. Nothing new happens.
  • Compute_Sflx_For_Vdiff has been moved to hco_interface_gc_mod.F90 due to circular dependency.
  • Loop orders in mixing_mod.F90 : DO_TEND and hco_interface_gc_mod.F90 : Compute_Sflx_For_Vdiff have been switched slightly: the species loop must be outermost, as HEMCO intermediate grid only stores up to one species in the "regridded" memory buffer, so species must be read contiguously before moving on to the next.
  • The changes above have been validated to minimal differences, attributed to switch of some precisions from f4 to fp, and optimization changes due to change of loop orders. See attached log below.

• When intermediate grid is enabled, the above HCO_GC_* calls retrieves data from HEMCO, compares it against the buffer, and regrids as necessary, discarding old data in the buffer. Why a buffer? Because some calls to HCO_GetPtr expect a POINTER and do their own data copying. To maximize compatibility, they're given the temporary buffer to point to. Care should be taken to "update" these pointers at every time step now, as they won't be pointing to the right things if the data wasn't copied from them before the next HCO_GC_* call. I've checked all the places in G-C that do this and changed everywhere, except strat_chem_mod.F90, pending decision on tropchem.

• Detailed list of changes below. I've also fixed some debug messages as I went probing through the entire GC source.

List of changes made to the code in all commits

Structural changes (validated to "zero" diff, see below, when IMGrid is not enabled):

• Implements a new Input_Opt field LIMGRID to indicate intermediate grid functionality.
• Reads HEMCO extension met fields from HEMCO state when possible
• Remove FRAC_OF_PBL duplicate and read through State_Met.
• Add IMGRID_XSCALE, IMGRID_YSCALE HEMCO grid scaling (refinement) factors to Input_Opt.
• Move SUMCOSZA to State_Met%SUNCOSsum.
• Move Calc_SUMCOSZA to Calc_Met_Mod.
• Move Read_Met_Fields to HCO_Interface_GC_Mod to avoid circular dependencies.
• Add State_Grid_HCO to store HEMCO grid state in HCO_State_GC_Mod.
• Modified method signatures in a few subroutines to accept Input_Opt, State_Grid.
• Add InquireHco subroutine in HCO_GC_Utilities_Mod for checking if given species has Deposition/Emissions, instead of calling with I,J,L=1 to retrieve data for probing.
• Now use EmisSpec, DepSpec in Compute_Sflx_For_Vdiff which retrieves a species is emis/dep enabled or not from the InquireHco subroutine, to save CPU.
• Moved Compute_Sflx_For_Vdiff from HCO_Utilities_GC_Mod to HCO_Interface_GC_Mod to avoid a circular dependence issue (the subroutine depends on Global_CH4 and Global_Br)
• Now read PNOXLOSS_O3 and PNOXLOSS_HNO3 diagnostics in MIXING_MOD and COMPUTE_SFLX_FOR_VDIFF at every time step, in preparation of HEMCO changes.
• Add HCO_GC_GetDiagn interface in HCO_Utilities_GC_Mod. It is a shim for GetHcoDiagn in HCO_Interface_Common, performing on-demand regrid when IMGrid is enabled.
• Add a secondary R4 buffer for model regridded data.
• Fixed wrong output type and added extra out variables for N_/L2_L2-cutoff error in FAST_JX_Mod.
• Add a debug message for Compute_Sflx_for_Vdiff in MAIN.F90

Functionality changes, all MODEL_CLASSIC:

• Add intermediate grid assignment in HCOI_GC_Init.
• Add intermediate grid variables and allocation in HCO_Interface_GC_Mod.
• Add Compute_Scaled_Grid to compute refined grid based on State_Grid in GC_Grid_Mod.
• Add on-demand regridding functionality for GetHcoValDep, GetHcoValEmis in HCO_Utilities_GC_Mod.
• Add on-demand regridding shim HCO_GC_GetPtr in HCO_Utilities_GC_Mod.
• Add on-demand regridding shim HCO_GC_EvalFld in HCO_Utilities_GC_Mod.
• Use on-demand regridding for reading met fields in FlexGrid and HCO modules.
• Use on-demand regridding in a variety of modules that used HCO_GetPtr.
• Use on-demand regridding for various modules (TOMS_MOD, TAGGED_O3, TAGGED_CO, Sulfate_Mod, SfcVmr_Mod)
• Use on-demand regridding for various modules (RPMARES_Mod, ISORROPIAII_Mod, Global_CH4_Mod, Global_Br_Mod)
• Use on-demand regridding for various modules (Aerosol_Mod, Carbon_Mod)
• Use HCO_GC_EvalFld for OLSON_LANDMAP_MOD and MODIS_LAI_MOD.
• Use on-demand regridding for Global_CH4_Mod::EmissCH4, Land_Mercury_Mod::BIOMASSHg, Mercury_Mod
• Use on-demand regridding for Sulfate_Mod HEMCO diagnostics

Difference test output

Only nonzero species attached.
Comparison was made against GEOS-Chem Classic "Wrapper", with GEOS-Chem Classic commit 28aa0358 in dev/13.0.0.

Variable               Ref=geosfp_2x25_standard Dev=dev_geosfp_2x25_standard Dev-Ref
SpeciesConc_Br2      : 1.3926416e-07          | 1.3926413e-07          | -2.842170943040401e-14
SpeciesConc_HI       : 5.0519042e-09          | 5.0519047e-09          | 4.440892098500626e-16
SpeciesConc_HNO3     : 0.0009782148           | 0.0009782154           | 5.820766091346741e-10
SpeciesConc_IHN4     : 1.115168e-09           | 1.1151681e-09          | 1.1102230246251565e-16
SpeciesConc_INPB     : 4.9958273e-09          | 4.995828e-09           | 8.881784197001252e-16
SpeciesConc_LVOC     : 2.276081e-10           | 2.2760828e-10          | 1.8041124150158794e-16
SpeciesConc_NH3      : 2.5541443e-05          | 2.5541423e-05          | -2.000888343900442e-11
SpeciesConc_NIT      : 9.795141e-06           | 9.7945585e-06          | -5.820766091346741e-10

Please feel free to squash the commits and let me know of any questions regarding this implementation... Apologies in advance for a 6,000 line commit in one PR. Hopefully there aren't too many merge conflicts. All changes are under the hood if the functionality isn't enabled at code level in hco_interface_gc_mod.F90 to set the XSCALE/YSCALE to non-1 value, so afterwards GEOS-Chem should look and work like GEOS-Chem as usual.

Have a nice day 😄
Haipeng

[jimmielin]

Hi @lizziel this should replace the original pull request #364. I've checked that it compiles and runs as-is. However I'm not sure how to do a full difference test like we did before using the GC Unit Tester.

As intended, the code should yield minimal differences (there are a few minor precision switches around from r4 to r8) when the following hardcoded parameters for scaling the HEMCO Grid are set, in hco_interface_gc_mod.F90:

    ! Initialize the intermediate grid descriptor.
    ! To disable the HEMCO intermediate grid feature, simply set this DY, DX to
    ! equal to State_Grid%DY, State_Grid%DX (e.g. 2x2.5, 4x5, ...)
    !
    ! TODO: Read in the grid parameters via input.geos. For now, hardcode the scale factor.
    Input_Opt%IMGRID_XSCALE = 1
    Input_Opt%IMGRID_YSCALE = 1

If both scale parameters are 1, then the entire code goes through a path exactly like the reference GEOS-Chem version (Input_Opt%LIMGRID = .false.) - if set to 2, HEMCO grid is 2x finer in that dimension, and so on.

Additionally, strat_chem_mod.F90 has been left untouched as in 13.0.0 there was a pending decision on what to do with that particular module. I see it is still there. May I know what are the plans for it? I will be happy to adapt the intermediate grid changes to strat_chem_mod.F90. Right now, if modules call HCO_GetPtr/HCO_EvalFld directly instead of their _GC_ counterparts, the code will act normally if HEMCO Grid is disabled (xscale, yscale = 1, 1) but the pointer sizes will be invalid if it is enabled, so code has to be adapted.

All modules that "save" the pointer from within HEMCO (using HCO_GetPtr and SAVE or in-module variables) and directly read from it at every time step need to be remodeled to copy the data instead (HCO_EvalFld is good). I've changed this in most of the modules in the standard simulation.

Please let me know what I can do to further help with this!

[lizziel]

Thanks @jimmielin. I will do the difference testing and some other tests. Thanks so much for dealing with the conflicts. Much appreciated!

[lizziel]

I am now getting zero diffs for GC-Classic benchmark simulation with this update. The only change I had to make was to revert to getting pointer rather than evaluating field when retrieving Olson land map data from HEMCO during land type initialization. This avoided triggering an error trap for PBLHEIGHT in HEMCO not yet associated.

I am testing with GCHP now. I had to make several minor changes for compilation, mostly regarding placement of MODEL_CLASSIC ifdefs. It now runs but I am seeing before and after differences in mixing. For example, first day average ozone zonal mean:

The differences are small but clearly not numerical noise. I am looking through the changes to see if anything pops out. @jimmielin, does anything come to mind on this?

[jimmielin]

Hi @lizziel thanks so much for doing this! This indeed looks weird. Just to confirm, does this only run PBL mixing or there's also emissions/deposition/chemistry? Nothing immediately comes to mind however it might be useful to look at HNO3 to see if this also happens there. If so it might be related to the way that diagnostics are now copied from the pointers and not directly used inside the loop:

    ! If using HEMCO Intermediate grid feature, then the call needs to be
    ! refreshed at every time step for regridding. (hplin, 6/21/20)
    IF ( .NOT. Input_Opt%LNLPBL ) THEN
      CALL HCO_GC_GetDiagn( Input_Opt, State_Grid, 'PARANOX_O3_DEPOSITION_FLUX', &
                            .FALSE.,   RC, Ptr2D = Ptr2D          )
      IF( ASSOCIATED( Ptr2D )) THEN
        ALLOCATE ( PNOxLoss_O3( State_Grid%NX, State_Grid%NY ), STAT=RC )
        PNOxLoss_O3(:,:) = Ptr2D(:,:)
      ENDIF
      Ptr2D => NULL()

      CALL HCO_GC_GetDiagn( Input_Opt, State_Grid, 'PARANOX_HNO3_DEPOSITION_FLUX',&
                            .FALSE.,   RC, Ptr2D = Ptr2D        )
      IF( ASSOCIATED( Ptr2D )) THEN
        ALLOCATE ( PNOxLoss_HNO3( State_Grid%NX, State_Grid%NY ), STAT=RC )
        PNOxLoss_HNO3(:,:) = Ptr2D(:,:)
      ENDIF
      Ptr2D => NULL()
    ENDIF

It may also be useful to look at some other species that are emitted (e.g. NO, CO) to see if it also happens there. If so it might be something related to the way the loops are rearranged, and maybe during the merge conflict resolution something was shifted around. I'll also look into this more closely to see if I messed something up.

Thanks again for looking into this!

[lizziel]

Diffs occur for many species, including NO and CO, with a similar signature. It only occurs if Turn on PBL mixing? in input.geos is F (Input_Opt%LTURB false). My tests so far keep LNLPBL true to match benchmark and help hone in on this issue. If everything except LTURB is true there are no differences.

Using LTURB as the condition narrows it down to subroutines vdiff or compute_sflx_for_vdiff. These diffs are not in GC-Classic which makes me think it is has to do with MODEL_CLASSIC ifdefs. I'll keep looking.

[lizziel]

I traced the problem to the new interface HCO_GC_GetDiagn. Optional arguments iCOL and iAF are not passed in the call within compute_sflx_for_vdiff but the case of not present is only handled if running GEOS-Chem Classic. This prevented getting the pointer to the deposition flux for O3 and HNO3 in GCHP. After making a fix to the placement of the MODEL_CLASSIC ifdef I am getting zero diffs in GCHP.

I will merge this PR as is into dev and make the various fixes (change evalFld to getPtr for Olson land type, compilation fixes for GCHP, and the HCO_GC_GetDiagn fix) in a subsequent commit so we can get it into 13.1.

[jimmielin]

Hi Lizzie, thank you for the update! This sounds great. Please let me know if I can help with any other fixes.