Adapt to use GCHP outputs for IMI

[1Dandan]

Name and Institution (Required)

Name: Dandan Zhang
Institution: Harvard University

Describe the update

These updates are primarily for using stretched GCHP simulation for Jacobian runs in IMI, with several additional features/fixes for both IMI-GCC and IMI-GCHP. All modifications are on top of IMI v2.2.0.

Specific modifications for IMI-GCHP Jacobian runs:

As the feature of scaling grid boxes with specific value from a text file in HEMCO does not work with GCHP, as stated in a ticket of geoschem/GCHP#432, I modified the HEMCO code to apply gridded scaling factor based on the value of state vector ID at geoschem/HEMCO#330. This update of HEMCO would work for both GCC and GCHP. Thus, to minimize the differences between IMI-GCC and IMI-GCHP for synchronizing updates later, I modified the code for creating Jacobian runs in IMI-GCC as well. I tested with a Global GCC 2x2.5 and it works well. The changes would be:

• Preparing the GCHP restart file regridded from a satellite restart file with a regridding function.
• Creating gridded perturbation scaling factors for emissions by src/components/jacobian_component/make_perturbation_sf.py and gridded OH perturbation (If OptimizeOH is enabled)
• Selecting a scaling factor value based on state vector ID to apply to that state vector element with the functionality enabled by Modify for GCHP IMI perturbation runs HEMCO#330.
• The ExtData.rc and HEMCO_Config.rc for GCHP carbon simulation for $\rm CH_4$ should be updated by the PR at Fix for emission configurations for carbon simulations geos-chem#2893. In that PR, I also fixed some time-range and unit errors I found for GCC carbon-run and methane-run. The updated config would not interfere with other carbon simulation, and thus should be merged generally for carbon fix too.

Specific modifications for IMI-GCHP to create a grid file first:

• GridSpec and SparSelt were used to create state vector file with grid box area calculation in advance as implemented in src/utilities/make_state_vector_file.py.

Specific modifications for IMI-GCHP inversion (calculating Jacobian matrix):

• Utilities are updated in src/inversion_scripts/classify_TROPOMI_obs_to_CSgrids.py, src/inversion_scripts/operators/TROPOMI_operator.py, and src/inversion_scripts/operators/operator_utilities.py.
• Using SphericalPolygon for calculating overlap area for easier handling across the dateline. However, SphericalPolygon is 3D polygon and thus would be much slower than Shapely Polygon. Currently I am restricting the usage of SphericalPolygon only when the GC grid cross the dateline for GCHP for both polygon_gchp and polygon_tropomi, while Shapely Polygon was used elsewhere.

Fixes/features for both IMI-GCC and IMI-GCHP:

• Use time-averaged outputs for methane concentration and level edge diagnostics instead of instantaneous, as satellite observations are averaged to GC grid hourly to get super observations. Using time-average mode would not be that different or less preferred than the instantaneous mode. Using time-average mode would also eliminate the need to append outputs from spin-up for the first hour diagnostics.
• As a result of eliminating the need of appending outputs from spin-up, I modified the IMI code to disable the output from spinup run as well as we only need a new restart file from spin-up.
• Excluding permanent ice for both original TROPOMI and blended satellite observations using the surface classification from TROPOMI dataset. It is equivalent to Megan's filter as specified in Filter out permanent snow/ice scenes #349 for blended satellite observations, but this filter is applied to both blended and operational TROPOMI dataset in this fix.
• Excluding permanent ice or lakes when constructing state vector file.
• Add an option for optimizing soil when using normal errors. Note this option should not be used with lognormal errors as negative flux is not allowed with lognormal error. This option should not be used for regional IMI with buffer clusters as well, as using median scaling factor could potentially lead to large negative flux is one of the buffer state vector elements is negative.
• Update_prior_emis.py at src/components/posterior_component/update_prior_emis.py was never used. Instead, this function was called from https://github.com/geoschem/integrated_methane_inversion/blob/main/src/inversion_scripts/utils.py#L552. Thus, I deleted this file.
• The perturbation scaling factor starting at 2000 would interfere with the scaling factors used in CEDS 0.1 for carbon simulations, and thus modify it to start at 3000.
• For ReDoJacobian, check if the last-day output file is valid with 24 time entries.
• Format gridded netCDF files to have correct time dimension to be compatible with MAPL ExtData to read, including state vector, gridded perturbation scaling factor, gridded OH scaling factor, and prior emission outputs.
• In imi_preview.py, the inversion area cannot be calculated correctly with the function used there. Instead, an alternate method of estimating number of grid box is used as implemented in src/inversion_scripts/imi_preview.py. This will fix the error stated in Potential bugs for the calculation for averaging kernel sensitivity and cost estimation in preview? #331 for the cost estimation in imi_preview.py.
• Fix the error for calculating averaging kernel density as well as stated in Potential bugs for the calculation for averaging kernel sensitivity and cost estimation in preview? #331.
• The message printed out for TROPOMI file time range would be one day more, and thus confusing. Fix it accordingly.
• Add the extraction of GC_index during the calculation of Jacobian matrix to be used later for super observation plotting.
• Plotting super observations on GC grid, and add plotting for soil absorption when OptimizeSoil is turned on.

Yaml file changes: (An example yaml file)

• For GCHP specific changes, two additional entries of UseGCHP and STRETCH_GRID were added to config.yml. Need to add two lines of: UseGCHP: false and STRETCH_GRID: false when using GCC too, as these are mandatory to know which code should be applied.
• Same applies to OptimizeSoil too. We need to add an entry of OptimizeSoil: false in config.yml if not using OptimizeSoil.
• Use blended TROPOMI dataset by default.
• Need to build a new environment with a new yaml file supplied to use SphericalPolygon package along with other tools needed for GCHP.

Instructions/fixes to use stretched GCHP simulations

Note: Below would be fixed in a new version of GC or MAPL ExtData, and thus you only need to make these modifications on top of the raw GCHP v14.6.2 or GCC v14.6.2.

• After cloning GCHP v14.6.2 or GCC v14.6.2 (if you want to use what I modified for GCC) git repository, we need to update the geos-chem repository with changes for HEMCO_Config.rc and ExtData.rc (checking out this one), and update HEMCO with changes in this one.
• Update the scripts in <your-gchp-repository>/src/MAPL/gridcomps/ExtData/ExtDataGridCompMod.F90 as (removing the conversion of radians to degrees) (Specific to GCHP):

!call MAPL_ConfigSetAttribute(cflocal,value=temp_real*MAPL_RADIANS_TO_DEGREES, label=trim(COMP_Name)//   MAPL_CF_COMPONENT_SEPARATOR//"TARGET_LON:",rc=status)
call MAPL_ConfigSetAttribute(cflocal,value=temp_real, label=trim(COMP_Name)//                            MAPL_CF_COMPONENT_SEPARATOR//"TARGET_LON:",rc=status)
...
!call MAPL_ConfigSetAttribute(cflocal,value=temp_real*MAPL_RADIANS_TO_DEGREES, label=trim(COMP_Name)//   MAPL_CF_COMPONENT_SEPARATOR//"TARGET_LAT:",rc=status)
call MAPL_ConfigSetAttribute(cflocal,value=temp_real, label=trim(COMP_Name)//                            MAPL_CF_COMPONENT_SEPARATOR//"TARGET_LAT:",rc=status)

Further features to be added:

• Currently supplying GCHP emission outputs would need an additional day. This would be fixed with a new version of MAPL ExtData with ExtData2G
• The current stretched GCHP simulation itself was problematic for version 14.6.2. Only use the feature for uniform GCHP now.