Add option for 0.125x0.15625 resolution following updates from Xiaoli…

…n Wang

The grid resolution option is expanded to include 0.125x0.15625. These
modifications were obtained from Xiaolin Wang. These simulations require
using the feature/0.125x0.15625_resolution in the GCClassic/src/GEOS-Chem
repository (see also geoschem/geos-chem#1980).

Signed-off-by: Melissa Sulprizio <mpayer@seas.harvard.edu>

config.yml

@@ -82,7 +82,8 @@ Gamma: 1.0
 PrecomputedJacobian: false
 
 ## Grid
-## Options are 0.25x0.3125 (GEOSFP only), 0.5x0.625, 2.0x2.5, or 4.0x5.0
+## Options are 0.125x0.15625 (GEOS-FP only), 0.25x0.3125 (GEOS-FP only),
+## 0.5x0.625, 2.0x2.5, or 4.0x5.0
 Res: "0.25x0.3125"
 
 ## Meteorology

envs/Harvard-Cannon/config.harvard-cannon.12km.yml

@@ -0,0 +1,217 @@
+## IMI configuration file
+## Documentation @ https://imi.readthedocs.io/en/latest/getting-started/imi-config-file.html
+
+## General
+RunName: "Test_IMI_Philadelphia"
+isAWS: false
+UseSlurm: true
+SafeMode: true
+S3Upload: false
+
+## Period of interest
+StartDate: 20220101
+EndDate: 20230101
+SpinupMonths: 1
+
+## Use blended TROPOMI+GOSAT data (true)? Or use operational TROPOMI data (false)?
+BlendedTROPOMI: false
+
+## Is this a regional inversion? Set to false for global inversion
+isRegional: true
+
+## Select two character region ID (for using pre-cropped meteorological fields)
+## Current options are listed below with ([lat],[lon]) bounds:
+## "AF" : Africa ([-37,40], [-20,53])
+## "AS" : Asia ([-11,55],[60,150]) 
+## "EU" : Europe ([33,61],[-30,70])
+## "ME" : Middle East ([12,50], [-20,70])
+## "NA" : North America ([10,70],[-140,-40])
+## "OC" : Oceania ([-50,5], [110,180])
+## "RU" : Russia ([41,83], [19,180])
+## "SA" : South America ([-59,16], [-88,-31])
+## "" : Use for global global simulation or custom regions
+## For example, if the region of interest is in Europe ([33,61],[-30,70]), select "EU".
+RegionID: "NA"
+
+## Region of interest
+## These lat/lon bounds are only used if CreateAutomaticRectilinearStateVectorFile: true
+## Otherwise lat/lon bounds are determined from StateVectorFile
+LonMin: -105
+LonMax: -103
+LatMin: 31
+LatMax: 33
+
+## Kalman filter options
+KalmanMode: false
+UpdateFreqDays: 7
+NudgeFactor: 0.1
+
+## State vector
+CreateAutomaticRectilinearStateVectorFile: true
+nBufferClusters: 8
+BufferDeg: 5
+LandThreshold: 0.25
+OffshoreEmisThreshold: 0
+OptimizeBCs: false
+OptimizeOH: false
+
+## Point source datasets
+## Used for visualizations and state vector clustering
+PointSourceDatasets: ["SRON"]
+
+## Clustering Options
+ReducedDimensionStateVector: false
+DynamicKFClustering: false
+ClusteringMethod: "kmeans"
+NumberOfElements: 45
+ForcedNativeResolutionElements: 
+ - [31.5, -104]
+
+## Custom state vector
+StateVectorFile: "/path/to/StateVector.nc"
+ShapeFile: "resources/shapefiles/tmp.shp"
+
+## Inversion
+## Note PriorError and PriorErrorOH are relative fractions (e.g. 0.5 = 50%)
+## and PriorErrorBCs is in ppb
+PriorError: 0.5
+PriorErrorBCs: 10.0
+PriorErrorOH: 0.5
+ObsError: 15
+Gamma: 1.0
+PrecomputedJacobian: false
+
+## Grid
+## Options are 0.125x0.15625 (GEOS-FP only), 0.25x0.3125 (GEOS-FP only),
+## 0.5x0.625, 2.0x2.5, or 4.0x5.0
+Res: "0.125x0.15625"
+
+## Meteorology
+## Options are GEOSFP or MERRA2
+Met: "GEOSFP"
+
+## Setup modules
+## Turn on/off different steps in setting up the inversion 
+RunSetup: true
+SetupTemplateRundir: true
+SetupSpinupRun: false
+SetupJacobianRuns: false
+SetupInversion: false
+SetupPosteriorRun: false
+
+## Run modules
+## Turn on/off different steps in performing the inversion
+DoPriorEmis: true
+DoSpinup: false
+DoJacobian: false
+DoInversion: false
+DoPosterior: false
+
+## IMI preview
+## NOTE: RunSetup must be true to run preview
+DoPreview: true
+DOFSThreshold: 0
+
+## Resource allocation settings for slurm jobs
+PriorMemory: 50000
+SimulationCPUs: 16
+SimulationMemory: 200000
+JacobianCPUs: 16
+JacobianMemory: 100000
+RequestedTime: "0-24:00"
+SchedulerPartition: "sapphire,huce_cascade,huce_intel,seas_compute,shared"
+
+## Max number of simultaneous Jacobian runs from the job array (-1: no limit)
+MaxSimultaneousRuns: -1
+
+## Specify number of Jacobians runs to perform
+## The default is -1 which will create and submit a jacobian run
+## for each state vector element. Specifying a value >= 1 will
+## combine state vector elements into fewer runs.
+NumJacobianRuns: -1
+
+##====================================================================
+##
+## Advanced Settings (optional)
+##
+##====================================================================
+
+## These settings are intended for advanced users who wish to:
+## a. modify additional GEOS-Chem options, or
+## b. run the IMI on a local cluster.
+## They can be ignored for any standard cloud application of the IMI.
+
+##--------------------------------------------------------------------
+## Additional settings for GEOS-Chem simulations
+##--------------------------------------------------------------------
+
+## Jacobian settings
+## Note PerturbValue and PerturbValueOH are relative scale factors and 
+## PerturbValueBCs is in ppb
+PerturbValue: 1.5
+PerturbValueOH: 1.5
+PerturbValueBCs: 10.0
+
+## Apply scale factors from a previous inversion?
+UseEmisSF: false
+UseOHSF: false
+
+## Save out hourly diagnostics from GEOS-Chem?
+## For use in satellite operators via post-processing -- required for TROPOMI
+## inversions
+HourlyCH4: true
+
+## Turn on planeflight diagnostic in GEOS-Chem?
+## For use in comparing GEOS-Chem against planeflight data. The path
+## to those data must be specified in input.geos.
+PLANEFLIGHT: false
+
+## Turn on old observation operators in GEOS-Chem?
+## These will save out text files comparing GEOS-Chem to observations, but have
+## to be manually incorporated into the IMI
+GOSAT: false
+TCCON: false
+AIRS: false
+
+##------------------------------------------------------------------
+## Settings for running on local cluster
+##------------------------------------------------------------------
+
+## Path for IMI runs and output
+OutputPath: "/n/holyscratch01/jacob_lab/$USER"
+
+## Path to GEOS-Chem input data
+DataPath: "/n/holyscratch01/external_repos/GEOS-CHEM/gcgrid/gcdata/ExtData"
+
+## Path to TROPOMI Data
+DataPathTROPOMI: "/n/holylfs05/LABS/jacob_lab/imi/ch4/tropomi"
+
+## Conda environment file
+## See envs/README to create the Conda environment specified below
+CondaEnv: "imi_env"
+
+## GEOS-Chem environment file (with fortran compiler, netcdf libraries, etc.)
+## NOTE: Copy your own file in the envs/ directory within the IMI
+GEOSChemEnv: "envs/Harvard-Cannon/gcclassic.rocky+gnu12.minimal.env"
+
+## Download initial restart file from AWS S3?
+## NOTE: Must have AWS CLI enabled
+RestartDownload: false
+
+## Path to initial GEOS-Chem restart file + prefix
+## ("YYYYMMDD_0000z.nc4" will be appended)
+RestartFilePrefix: "/n/holylfs05/LABS/jacob_lab/imi/ch4/tropomi-boundary-conditions/v2023-10/GEOSChem.BoundaryConditions."
+RestartFilePreviewPrefix: "/n/holylfs05/LABS/jacob_lab/imi/ch4/tropomi-boundary-conditions/v2023-10/GEOSChem.BoundaryConditions."
+
+## Path to GEOS-Chem boundary condition files (for regional simulations)
+## BCversion will be appended to the end of this path. ${BCpath}/${BCversion}
+BCpath: "/n/holylfs05/LABS/jacob_lab/imi/ch4/tropomi-boundary-conditions"
+BCversion: "v2023-10"
+
+## Options to download missing GEOS-Chem input data from AWS S3
+## NOTE: Must have AWS CLI enabled
+PreviewDryRun: false
+SpinupDryrun: false
+ProductionDryRun: false
+PosteriorDryRun: false
+BCdryrun: false

envs/Harvard-Cannon/config.harvard-cannon.global_inv.yml

@@ -82,7 +82,8 @@ Gamma: 1.0
 PrecomputedJacobian: false
 
 ## Grid
-## Options are 0.25x0.3125 (GEOSFP only), 0.5x0.625, 2.0x2.5, or 4.0x5.0
+## Options are 0.125x0.15625 (GEOS-FP only), 0.25x0.3125 (GEOS-FP only),
+## 0.5x0.625, 2.0x2.5, or 4.0x5.0
 Res: "2.0x2.5"
 
 ## Meteorology

envs/Harvard-Cannon/config.harvard-cannon.yml

@@ -82,7 +82,8 @@ Gamma: 1.0
 PrecomputedJacobian: false
 
 ## Grid
-## Options are 0.25x0.3125 (GEOSFP only), 0.5x0.625, 2.0x2.5, or 4.0x5.0
+## Options are 0.125x0.15625 (GEOS-FP only), 0.25x0.3125 (GEOS-FP only),
+## 0.5x0.625, 2.0x2.5, or 4.0x5.0
 Res: "0.25x0.3125"
 
 ## Meteorology

src/components/prior_component/prior.sh

@@ -40,10 +40,12 @@ run_prior() {
  resnum="3"
  elif [ "$Res" == "0.25x0.3125" ]; then
  resnum="4"
+ elif [ "$Res" == "0.125x0.15625" ]; then
+ resnum="4"
  else
- printf "\nERROR: Grid resolution ${Res} is not supported by the IMI. "
- printf "\n Options are 0.25x0.3125, 0.5x0.625, 2.0x2.5, or 4.0x5.0.\n"
- exit 1
+  printf "\nERROR: Grid resolution ${Res} is not supported by the IMI."
+  printf "\n Options are 0.125x0.15625, 0.25x0.3125, 0.5x0.625, 2.0x2.5, or 4.0x5.0.\n"
+  exit 1
  fi
  HEMCOconfig=${RunTemplate}/HEMCO_Config.rc
 

src/components/setup_component/setup.sh

@@ -97,7 +97,10 @@ setup_imi() {
  exit 1
  fi
 
- if [ "$Res" == "0.25x0.3125" ]; then
+ if [ "$Res" == "0.125x0.15625" ]; then
+ gridDir="${Res}"
+ gridFile="0125x015625" 
+ elif [ "$Res" == "0.25x0.3125" ]; then
  gridDir="${Res}"
  gridFile="025x03125"
  elif [ "$Res" == "0.5x0.625" ]; then
@@ -111,7 +114,7 @@ setup_imi() {
  gridFile="4x5"
  else
  printf "\nERROR: Grid resolution ${Res} is not supported by the IMI. "
- printf "\n Options are 0.25x0.3125, 0.5x0.625, 2.0x2.5, or 4.0x5.0.\n"
+ printf "\n Options are 0.125x0.15625, 0.25x0.3125, 0.5x0.625, 2.0x2.5, or 4.0x5.0.\n"
  exit 1
  fi
  # Use cropped met for regional simulations instead of using global met

src/components/statevector_component/aggregation.py

@@ -224,7 +224,9 @@ def get_max_aggregation_level(config, sensitivities, desired_element_num):
  desired_element_num int : desired number of state vector elements
  Returns: int : max gridcells per cluster
  """
- if config["Res"] == "0.25x0.3125":
+ if config["Res"] == "0.125x0.15625":
+ max_aggregation_level = 512
+ elif config["Res"] == "0.25x0.3125":
  max_aggregation_level = 256
  elif config["Res"] == "0.5x0.625":
  max_aggregation_level = 64
@@ -328,7 +330,10 @@ def force_native_res_pixels(config, clusters, sensitivities):
  )
  return sensitivities
 
- if config["Res"] == "0.25x0.3125":
+ if config["Res"] == "0.125x0.15625":
+ lat_step = 0.125
+ lon_step = 0.15625
+ elif config["Res"] == "0.25x0.3125":
  lat_step = 0.25
  lon_step = 0.3125
  elif config["Res"] == "0.5x0.625":

src/components/statevector_component/statevector.sh

@@ -12,12 +12,18 @@ create_statevector() {
 
  # Use GEOS-FP or MERRA-2 CN file to determine ocean/land grid boxes
  if "$isRegional"; then
- LandCoverFile="${DataPath}/GEOS_${gridDir}/${metDir}/${constYr}/01/${Met}.${constYr}0101.CN.${gridFile}.${RegionID}.${LandCoverFileExtension}"
+ if [ "$Res" = "0.125x0.15625" ]; then
+ LandCoverFile="/n/home00/xlwang/xlwang/methane_inversion/InputData/land_cover/data/IMERG_land_sea_mask_0125x015625.nc"
+ HemcoDiagFile="/n/home00/xlwang/xlwang/Test_12km/gc_0125x015625_NA_47L_geosfp_CH4_emission/OutputDir/HEMCO_diagnostics.202205010000.nc"
+ else
+ LandCoverFile="${DataPath}/GEOS_${gridDir}/${metDir}/${constYr}/01/${Met}.${constYr}0101.CN.${gridFile}.${RegionID}.${LandCoverFileExtension}"
+ HemcoDiagFile="${DataPath}/HEMCO/CH4/v2023-04/HEMCO_SA_Output/HEMCO_sa_diagnostics.${gridFile}.20190101.nc"
+ fi
  else
  LandCoverFile="${DataPath}/GEOS_${gridDir}/${metDir}/${constYr}/01/${Met}.${constYr}0101.CN.${gridFile}.${LandCoverFileExtension}"
+ HemcoDiagFile="${DataPath}/HEMCO/CH4/v2023-04/HEMCO_SA_Output/HEMCO_sa_diagnostics.${gridFile}.20190101.nc"
  fi
- HemcoDiagFile="${DataPath}/HEMCO/CH4/v2023-04/HEMCO_SA_Output/HEMCO_sa_diagnostics.${gridFile}.20190101.nc"
-
+
  if "$isAWS"; then
  # Download land cover and HEMCO diagnostics files
  s3_lc_path="s3://gcgrid/GEOS_${gridDir}/${metDir}/${constYr}/01/${Met}.${constYr}0101.CN.${gridFile}.${RegionID}.${LandCoverFileExtension}"

src/components/template_component/template.sh

@@ -53,9 +53,11 @@ setup_template() {
  else
  cmd="3\n${metNum}\n4\n1\n2\n${RunDirs}\n${runDir}\nn\n"
  fi
+ elif [ "$Res" == "0.125x0.15625" ]; then
+ cmd="3\n${metNum}\n5\n4\n2\n${RunDirs}\n${runDir}\nn\n" #regional run
  else
  printf "\nERROR: Grid resolution ${Res} is not supported by the IMI. "
- printf "\n Options are 0.25x0.3125, 0.5x0.625, 2.0x2.5, or 4.0x5.0.\n"
+ printf "\n Options are 0.125x0.15625, 0.25x0.3125, 0.5x0.625, 2.0x2.5, or 4.0x5.0.\n"
  exit 1
  fi
 
@@ -115,6 +117,15 @@ setup_template() {
  sed -i "s/$OLD/$NEW/g" geoschem_config.yml
  fi
 
+ # Modify the METDIR for 0.125x0.15625 simulation
+ if [ "$Res" = "0.125x0.15625" ]; then
+ sed -i -e "s:GEOS_0.25x0.3125\/GEOS_FP:GEOS_0.25x0.3125_NA\/GEOS_FP:g" HEMCO_Config.rc.gmao_metfields_0125
+ OLD="/n/holyscratch01/external_repos/GEOS-CHEM/gcgrid/gcdata/ExtData/GEOS_0.125x0.15625/GEOS_FP"
+ NEW="/n/holylfs05/LABS/jacob_lab/Users/xlwang/methane_inversion/InputData/GEOS_0.125x0.15625_NA/GEOS_FP_DerivedWinds"
+ sed -i "s|$OLD|$NEW|g" HEMCO_Config.rc.gmao_metfields_0125
+ sed -i '/METDIR/d' HEMCO_Config.rc
+ fi
+
  # Modify HEMCO_Config.rc based on settings in config.yml
  # Use cropped met fields (add the region to both METDIR and the met files)
  if "$isRegional"; then

src/inversion_scripts/imi_preview.py

@@ -165,7 +165,9 @@ def imi_preview(
  ]
  inversion_area_km = calculate_area_in_km(coords)
 
- if config["Res"] == "0.25x0.3125":
+ if config["Res"] == "0.125x0.15625":
+ res_factor = 2
+ elif config["Res"] == "0.25x0.3125":
  res_factor = 1
  elif config["Res"] == "0.5x0.625":
  res_factor = 0.5
@@ -492,7 +494,11 @@ def estimate_averaging_kernel(
 
  # Set resolution specific variables
  # L_native = Rough length scale of native state vector element [m]
- if config["Res"] == "0.25x0.3125":
+ if config["Res"] == "0.125x0.15625":
+ L_native = 12 * 1000
+ lat_step = 0.125
+ lon_step = 0.15625
+ elif config["Res"] == "0.25x0.3125":
  L_native = 25 * 1000
  lat_step = 0.25
  lon_step = 0.3125

src/inversion_scripts/invert.py

@@ -56,8 +56,11 @@ def do_inversion(
  # Need to ignore data in the GEOS-Chem 3 3 3 3 buffer zone
  # Shave off one or two degrees of latitude/longitude from each side of the domain
  # ~1 degree if 0.25x0.3125 resolution, ~2 degrees if 0.5x0.6125 resolution
- # This assumes 0.25x0.3125 and 0.5x0.625 simulations are always regional
- if "0.25x0.3125" in res:
+ # This assumes 0.125x0.15625, 0.25x0.3125, and 0.5x0.625 simulations are always regional
+ if "0.125x0.15625" in res:
+ degx = 4 * 0.15625
+ degy = 4 * 0.125
+ elif "0.25x0.3125" in res:
  degx = 4 * 0.3125
  degy = 4 * 0.25
  elif "0.5x0.625" in res: