geoschem · yantosca · Oct 6, 2023 · Sep 27, 2023 · Oct 4, 2023 · Oct 6, 2023
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -16,6 +16,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 - Added warning in GCHP HISTORY.rc about outputting area-dependent variables on custom grids
 - Added option to use a single advected species in the carbon simulation
 - Added option to perturb CH4 boundary conditions in CH4 simulation
+- Added option to perturb OH in CH4 simulation using scale factor in HEMCO_Config.rc
 
 ### Changed
 - Update `DiagnFreq` in GCClassic integration tests to ensure HEMCO diagnostic output

diff --git a/GeosCore/global_ch4_mod.F90 b/GeosCore/global_ch4_mod.F90
@@ -149,7 +149,7 @@ SUBROUTINE EMISSCH4( Input_Opt, State_Chm, State_Grid, State_Met, RC )
     ! =================================================================
     ! Get fields for CH4 analytical inversions if needed
     ! =================================================================
-    IF ( Input_Opt%AnalyticalInv ) THEN
+    IF ( Input_Opt%DoAnalyticalInv ) THEN
 
        ! Evaluate the state vector field from HEMCO
        CALL HCO_GC_EvalFld( Input_Opt, State_Grid, 'CH4_STATE_VECTOR', &
@@ -506,23 +506,22 @@ SUBROUTINE EMISSCH4( Input_Opt, State_Chm, State_Grid, State_Met, RC )
     ! =================================================================
     ! Do scaling for analytical inversion
     ! =================================================================
-    IF ( Input_Opt%AnalyticalInv  .or. &
+    IF ( Input_Opt%DoAnalyticalInv  .or. &
          Input_Opt%UseEmisSF      .or. &
          Input_Opt%UseOHSF        ) THEN
 
        ! Don't optimize for soil absorption so remove from the total
        ! emissions array
        State_Chm%CH4_EMIS(:,:,1) = State_Chm%CH4_EMIS(:,:,1) + State_Chm%CH4_EMIS(:,:,15)
 
-       ! Rescale emissions
        !$OMP PARALLEL DO       &
        !$OMP DEFAULT( SHARED ) &
        !$OMP PRIVATE( I, J)	 
        DO J = 1, State_Grid%NY
        DO I = 1, State_Grid%NX
 
           !------------------------------------------------------------
-          ! For applying scale factors from analytical inversion
+          ! Apply emission scale factors from a previous inversion
           !------------------------------------------------------------
           IF ( Input_Opt%UseEmisSF ) THEN
              ! Scale total emissions
@@ -532,16 +531,20 @@ SUBROUTINE EMISSCH4( Input_Opt, State_Chm, State_Grid, State_Met, RC )
           !------------------------------------------------------------
           ! Perturb emissions for analytical inversion
           !------------------------------------------------------------
-          IF ( Input_Opt%AnalyticalInv ) THEN
+          IF ( Input_Opt%DoAnalyticalInv ) THEN
 
              ! Only apply emission perturbation to current state vector
              ! element number
              IF ( Input_Opt%StateVectorElement .GT. 0 ) THEN
                 IF ( STATE_VECTOR(I,J) == Input_Opt%StateVectorElement ) THEN
-                   State_Chm%CH4_EMIS(I,J,1) = State_Chm%CH4_EMIS(I,J,1) * Input_Opt%PerturbEmis
-                   !Print*, 'Analytical Inversion: Scaled state vector element ', &
-                   !        Input_Opt%StateVectorElement, ' by ', &
-                   !        Input_Opt%PerturbEmis
+                   State_Chm%CH4_EMIS(I,J,1) = State_Chm%CH4_EMIS(I,J,1) * &
+                                               Input_Opt%EmisPerturbFactor
+                   IF ( Input_Opt%Verbose ) THEN
+                      Print*, 'Analytical Inversion: ',              &
+                              'Scaled state vector element ',       &
+                              Input_Opt%StateVectorElement, ' by ', &
+                              Input_Opt%EmisPerturbFactor
+                   ENDIF
                 ENDIF
              ENDIF
           ENDIF
@@ -930,10 +933,13 @@ SUBROUTINE CH4_DECAY( Input_Opt,  State_Chm, State_Diag,  &
           ! BOH from HEMCO in units of kg/m3, convert to molec/cm3
           C_OH = State_Chm%BOH(I,J,L) * XNUMOL_OH / CM3PERM3
 
-          ! Apply scale factors from analytical inversion
+          ! Apply OH scale factors from a previous inversion
           IF ( Input_Opt%UseOHSF ) THEN
              C_OH = C_OH * OH_SF(I,J)
-             !Print*, 'Applying scale factor to OH: ', OH_SF(I,J)
+             IF ( Input_Opt%Verbose ) THEN
+                !This will print over every grid box; comment out for now
+                !Print*, 'Applying scale factor to OH: ', OH_SF(I,J)
+             ENDIF
           ENDIF
 
           ! Cl in [molec/cm3]
@@ -973,9 +979,11 @@ SUBROUTINE CH4_DECAY( Input_Opt,  State_Chm, State_Diag,  &
     ENDDO
     !$OMP END PARALLEL DO
 
-    !print*,'% --- CHEMCH4: CH4_DECAY: TROP DECAY (Tg): ',TROPOCH4/1e9
-    !print*,'Trop decay should be over 1Tg per day globally'
-    !print*,'    ~ 500Tg/365d ~ 1.37/d'
+    IF ( Input_Opt%Verbose ) THEN
+       print*,'% --- CHEMCH4: CH4_DECAY: TROP DECAY (Tg): ',TROPOCH4/1e9
+       print*,'Trop decay should be over 1Tg per day globally'
+       print*,'    ~ 500Tg/365d ~ 1.37/d'
+    ENDIF
 
     ! Free pointers
     Spc => NULL()

diff --git a/GeosCore/hco_interface_gc_mod.F90 b/GeosCore/hco_interface_gc_mod.F90
@@ -4171,7 +4171,7 @@ SUBROUTINE CheckSettings( HcoConfig, Input_Opt, State_Met, State_Chm, RC )
     !-----------------------------------------------------------------------
     IF ( Input_Opt%ITS_A_CH4_SIM .or. Input_Opt%ITS_A_TAGCH4_SIM) THEN
 
-       IF ( Input_Opt%AnalyticalInv ) THEN
+       IF ( Input_Opt%DoAnalyticalInv ) THEN
           CALL GetExtOpt( HcoConfig, -999, 'AnalyticalInv', &
                           OptValBool=LTMP, FOUND=FOUND,  RC=HMRC )
 

diff --git a/GeosCore/hco_utilities_gc_mod.F90 b/GeosCore/hco_utilities_gc_mod.F90
@@ -2576,7 +2576,7 @@ SUBROUTINE Get_Boundary_Conditions( Input_Opt, State_Chm, State_Grid, &
       ! Each time the boundary conditions are read, they have no longer been perturbed
       ! This value is sometimes set to True in set_boundary_conditions_mod.F90
       IF ( ( Input_Opt%ITS_A_CH4_SIM .OR. Input_Opt%ITS_A_CARBON_SIM ) .AND. & 
-           Input_Opt%PerturbCH4BoundaryConditions ) THEN
+           Input_Opt%DoPerturbCH4BoundaryConditions ) THEN
         State_Chm%IsCH4BCPerturbed = .FALSE.
       ENDIF
 

diff --git a/GeosCore/input_mod.F90 b/GeosCore/input_mod.F90
@@ -4899,20 +4899,7 @@ SUBROUTINE Config_CH4( Config, Input_Opt, RC )
        CALL GC_Error( errMsg, RC, thisLoc )
        RETURN
     ENDIF
-    Input_Opt%AnalyticalInv = v_bool
-
-    !------------------------------------------------------------------------
-    ! Emission perturbation
-    !------------------------------------------------------------------------
-    key   = "CH4_simulation_options%analytical_inversion%emission_perturbation"
-    v_str = MISSING_STR
-    CALL QFYAML_Add_Get( Config, TRIM( key ), v_str, "", RC )
-    IF ( RC /= GC_SUCCESS ) THEN
-       errMsg = 'Error parsing ' // TRIM( key ) // '!'
-       CALL GC_Error( errMsg, RC, thisLoc )
-       RETURN
-    ENDIF
-    Input_Opt%PerturbEmis = Cast_and_RoundOff( v_str, places=4 )
+    Input_Opt%DoAnalyticalInv = v_bool
 
     !------------------------------------------------------------------------
     ! Current state vector element number
@@ -4929,31 +4916,17 @@ SUBROUTINE Config_CH4( Config, Input_Opt, RC )
     Input_Opt%StateVectorElement = v_int
 
     !------------------------------------------------------------------------
-    ! Use emission scale factor?
+    ! Emission perturbation factor
     !------------------------------------------------------------------------
-    key = &
-     "CH4_simulation_options%analytical_inversion%use_emission_scale_factor"
-    v_bool = MISSING_BOOL
-    CALL QFYAML_Add_Get( Config, TRIM( key ), v_bool, "", RC )
-    IF ( RC /= GC_SUCCESS ) THEN
-       errMsg = 'Error parsing ' // TRIM( key ) // '!'
-       CALL GC_Error( errMsg, RC, thisLoc )
-       RETURN
-    ENDIF
-    Input_Opt%UseEmisSF = v_bool
-
-    !------------------------------------------------------------------------
-    ! Use OH scale factors?
-    !------------------------------------------------------------------------
-    key    = "CH4_simulation_options%analytical_inversion%use_OH_scale_factors"
-    v_bool = MISSING_BOOL
-    CALL QFYAML_Add_Get( Config, TRIM( key ), v_bool, "", RC )
+    key   = "CH4_simulation_options%analytical_inversion%emission_perturbation_factor"
+    v_str = MISSING_STR
+    CALL QFYAML_Add_Get( Config, TRIM( key ), v_str, "", RC )
     IF ( RC /= GC_SUCCESS ) THEN
        errMsg = 'Error parsing ' // TRIM( key ) // '!'
        CALL GC_Error( errMsg, RC, thisLoc )
        RETURN
     ENDIF
-    Input_Opt%UseOHSF = v_bool
+    Input_Opt%EmisPerturbFactor = Cast_and_RoundOff( v_str, places=4 )
 
     !------------------------------------------------------------------------
     ! Perturb CH4 boundary conditions?
@@ -4966,7 +4939,7 @@ SUBROUTINE Config_CH4( Config, Input_Opt, RC )
        CALL GC_Error( errMsg, RC, thisLoc )
        RETURN
     ENDIF
-    Input_Opt%PerturbCH4BoundaryConditions = v_bool
+    Input_Opt%DoPerturbCH4BoundaryConditions = v_bool
 
     !------------------------------------------------------------------------
     ! How much to perturb CH4 boundary conditions by?
@@ -4984,25 +4957,52 @@ SUBROUTINE Config_CH4( Config, Input_Opt, RC )
     Input_Opt%CH4BoundaryConditionIncreaseEast  = Cast_and_RoundOff( a_str(3), places=4 )
     Input_Opt%CH4BoundaryConditionIncreaseWest  = Cast_and_RoundOff( a_str(4), places=4 )
 
+    !------------------------------------------------------------------------
+    ! Use emission scale factors from a previous inversion?
+    !------------------------------------------------------------------------
+    key = &
+     "CH4_simulation_options%analytical_inversion%use_emission_scale_factor"
+    v_bool = MISSING_BOOL
+    CALL QFYAML_Add_Get( Config, TRIM( key ), v_bool, "", RC )
+    IF ( RC /= GC_SUCCESS ) THEN
+       errMsg = 'Error parsing ' // TRIM( key ) // '!'
+       CALL GC_Error( errMsg, RC, thisLoc )
+       RETURN
+    ENDIF
+    Input_Opt%UseEmisSF = v_bool
+
+    !------------------------------------------------------------------------
+    ! Use OH scale factors from a previous inversion?
+    !------------------------------------------------------------------------
+    key    = "CH4_simulation_options%analytical_inversion%use_OH_scale_factors"
+    v_bool = MISSING_BOOL
+    CALL QFYAML_Add_Get( Config, TRIM( key ), v_bool, "", RC )
+    IF ( RC /= GC_SUCCESS ) THEN
+       errMsg = 'Error parsing ' // TRIM( key ) // '!'
+       CALL GC_Error( errMsg, RC, thisLoc )
+       RETURN
+    ENDIF
+    Input_Opt%UseOHSF = v_bool
+
     !========================================================================
     ! Print to screen
     !========================================================================
     IF ( Input_Opt%amIRoot ) THEN
        WRITE(6,90 ) 'CH4 SIMULATION SETTINGS'
        WRITE(6,95 ) '-----------------------'
-       WRITE(6,100) 'Use AIRS obs operator    : ', Input_Opt%AIRS_CH4_OBS
-       WRITE(6,100) 'Use GOSAT obs operator   : ', Input_Opt%GOSAT_CH4_OBS
-       WRITE(6,100) 'Use TCCON obs operator   : ', Input_Opt%TCCON_CH4_OBS
-       WRITE(6,100) 'Do analytical inversion  : ', Input_Opt%AnalyticalInv
-       WRITE(6,110) 'Emission perturbation    : ', Input_Opt%PerturbEmis
+       WRITE(6,100) 'Use AIRS obs operator?   : ', Input_Opt%AIRS_CH4_OBS
+       WRITE(6,100) 'Use GOSAT obs operator?  : ', Input_Opt%GOSAT_CH4_OBS
+       WRITE(6,100) 'Use TCCON obs operator?  : ', Input_Opt%TCCON_CH4_OBS
+       WRITE(6,100) 'Do analytical inversion? : ', Input_Opt%DoAnalyticalInv
        WRITE(6,120) 'Current state vector elem: ', Input_Opt%StateVectorElement
-       WRITE(6,100) 'Use emis scale factors   : ', Input_Opt%UseEmisSF
-       WRITE(6,100) 'Use OH scale factors     : ', Input_Opt%UseOHSF
-       WRITE(6,100) 'Perturb CH4 BCs          : ', Input_Opt%PerturbCH4BoundaryConditions
+       WRITE(6,110) 'Emiss perturbation factor: ', Input_Opt%EmisPerturbFactor
+       WRITE(6,100) 'Perturb CH4 BCs?         : ', Input_Opt%DoPerturbCH4BoundaryConditions
        WRITE(6,130) 'CH4 BC ppb increase NSEW : ', Input_Opt%CH4BoundaryConditionIncreaseNorth,&
                                                    Input_Opt%CH4BoundaryConditionIncreaseSouth,&
                                                    Input_Opt%CH4BoundaryConditionIncreaseEast,&
                                                    Input_Opt%CH4BoundaryConditionIncreaseWest
+       WRITE(6,100) 'Use emis scale factors?  : ', Input_Opt%UseEmisSF
+       WRITE(6,100) 'Use OH scale factors?    : ', Input_Opt%UseOHSF
     ENDIF
 
     ! FORMAT statements

diff --git a/GeosCore/set_boundary_conditions_mod.F90 b/GeosCore/set_boundary_conditions_mod.F90
@@ -141,7 +141,7 @@ SUBROUTINE Set_Boundary_Conditions( Input_Opt, State_Chm, State_Grid, RC )
       ! Convert to [kg/kg dry] (nbalasus, 8/31/2023)
       Perturb_CH4_BC = ( State_Chm%SpcData(N)%Info%Name == "CH4" .AND. &
                          ( Input_Opt%ITS_A_CH4_SIM .OR. Input_Opt%ITS_A_CARBON_SIM ) .AND. &
-                         Input_Opt%PerturbCH4BoundaryConditions .AND. &
+                         Input_Opt%DoPerturbCH4BoundaryConditions .AND. &
                          ( .NOT. State_Chm%IsCH4BCPerturbed ) )
       MW_g_CH4       =   State_Chm%SpcData(N)%Info%MW_g
 

diff --git a/Headers/input_opt_mod.F90 b/Headers/input_opt_mod.F90
@@ -365,16 +365,16 @@ MODULE Input_Opt_Mod
      LOGICAL                     :: GOSAT_CH4_OBS
      LOGICAL                     :: AIRS_CH4_OBS
      LOGICAL                     :: TCCON_CH4_OBS
-     LOGICAL                     :: AnalyticalInv
-     REAL(fp)                    :: PerturbEmis
+     LOGICAL                     :: DoAnalyticalInv
      INTEGER                     :: StateVectorElement
-     LOGICAL                     :: UseEmisSF
-     LOGICAL                     :: UseOHSF
-     LOGICAL                     :: PerturbCH4BoundaryConditions
+     REAL(fp)                    :: EmisPerturbFactor
+     LOGICAL                     :: DoPerturbCH4BoundaryConditions
      REAL(fp)                    :: CH4BoundaryConditionIncreaseNorth
      REAL(fp)                    :: CH4BoundaryConditionIncreaseSouth
      REAL(fp)                    :: CH4BoundaryConditionIncreaseEast
      REAL(fp)                    :: CH4BoundaryConditionIncreaseWest
+     LOGICAL                     :: UseEmisSF
+     LOGICAL                     :: UseOHSF
 
      !----------------------------------------
      ! POPS MENU fields
@@ -906,16 +906,16 @@ SUBROUTINE Set_Input_Opt( am_I_Root, Input_Opt, RC )
     Input_Opt%GOSAT_CH4_OBS                     = .FALSE.
     Input_Opt%AIRS_CH4_OBS                      = .FALSE.
     Input_Opt%TCCON_CH4_OBS                     = .FALSE.
-    Input_Opt%AnalyticalInv                     = .FALSE.
-    Input_Opt%PerturbEmis                       = 1.0
+    Input_Opt%DoAnalyticalInv                   = .FALSE.
     Input_Opt%StateVectorElement                = 0
-    Input_Opt%UseEmisSF                         = .FALSE.
-    Input_Opt%UseOHSF                           = .FALSE.
-    Input_Opt%PerturbCH4BoundaryConditions      = .FALSE.
+    Input_Opt%EmisPerturbFactor                 = 1.0
+    Input_Opt%DoPerturbCH4BoundaryConditions    = .FALSE.
     Input_Opt%CH4BoundaryConditionIncreaseNorth = 0.0_fp
     Input_Opt%CH4BoundaryConditionIncreaseSouth = 0.0_fp
     Input_Opt%CH4BoundaryConditionIncreaseEast  = 0.0_fp
     Input_Opt%CH4BoundaryConditionIncreaseWest  = 0.0_fp
+    Input_Opt%UseEmisSF                         = .FALSE.
+    Input_Opt%UseOHSF                           = .FALSE.
 
     !----------------------------------------
     ! POPS MENU fields

diff --git a/run/GCClassic/HEMCO_Config.rc.templates/HEMCO_Config.rc.CH4 b/run/GCClassic/HEMCO_Config.rc.templates/HEMCO_Config.rc.CH4
@@ -656,7 +656,7 @@ ${RUNDIR_CH4_LOSS}
 
 # --- Global OH from GEOS-Chem v5-07 [kg/m3] ---
 (((GLOBAL_OH
-* GLOBAL_OH  $ROOT/OH/v2014-09/v5-07-08/OH_3Dglobal.geos5.47L.4x5.nc OH           1985/1-12/1/0 C xyz kg/m3 * - 1 1
+* GLOBAL_OH  $ROOT/OH/v2014-09/v5-07-08/OH_3Dglobal.geos5.47L.4x5.nc OH           1985/1-12/1/0 C xyz kg/m3 * 2 1 1
 )))GLOBAL_OH
 
 # --- Global Cl [mol/mol dry air] ---
@@ -884,6 +884,14 @@ ${RUNDIR_GLOBAL_Cl}
 #==============================================================================
 1 NEGATIVE       -1.0 - - - xy 1 1
 
+#==============================================================================
+# --- Perturbation factors ---
+#
+# Add factors to perturb OH, emissions, and other fields here for
+# analytical inversions.
+#==============================================================================
+2 OH_pert_factor  1.0 - - - xy 1 1
+
 #==============================================================================
 # --- Seasonal scaling factors ----
 #==============================================================================

diff --git a/run/GCClassic/HEMCO_Config.rc.templates/HEMCO_Config.rc.carbon b/run/GCClassic/HEMCO_Config.rc.templates/HEMCO_Config.rc.carbon
@@ -1296,7 +1296,7 @@ Mask fractions:              false
 #------------------------------------------------------------------------------
 # --- OH from GEOS-Chem v5-07 [kg/m3], needed for CH4/IMI ---
 (((GLOBAL_OH_GCv5
-* GLOBAL_OH  $ROOT/OH/v2022-11/v5-07-08/OH_3Dglobal.geos5.47L.4x5.nc OH  1985/1-12/1/0 C xyz kg/m3 * - 1 1
+* GLOBAL_OH  $ROOT/OH/v2022-11/v5-07-08/OH_3Dglobal.geos5.47L.4x5.nc OH  1985/1-12/1/0 C xyz kg/m3 * 2 1 1
 )))GLOBAL_OH_GCv5
 
 # --- OH from the last 10-yr benchmark [mol/mol dry air] ---
@@ -1376,6 +1376,14 @@ ${RUNDIR_CO2_COPROD}
 #==============================================================================
 1 NEGATIVE -1.0 - - - xy 1 1
 
+#==============================================================================
+# --- Perturbation factors ---
+#
+# Add factors to perturb OH, emissions, and other fields here for
+# analytical inversions.
+#==============================================================================
+2 OH_pert_factor  1.0 - - - xy 1 1
+
 #==============================================================================
 # --- Seasonal scaling factors ----
 #==============================================================================

diff --git a/run/GCClassic/HEMCO_Config.rc.templates/HEMCO_Config.rc.tagCH4 b/run/GCClassic/HEMCO_Config.rc.templates/HEMCO_Config.rc.tagCH4
@@ -718,7 +718,7 @@ ${RUNDIR_CH4_LOSS}
 
 # --- Global OH from GEOS-Chem v5-07 [kg/m3] ---
 (((GLOBAL_OH
-* GLOBAL_OH  $ROOT/OH/v2014-09/v5-07-08/OH_3Dglobal.geos5.47L.4x5.nc OH           1985/1-12/1/0 C xyz kg/m3 * - 1 1
+* GLOBAL_OH  $ROOT/OH/v2014-09/v5-07-08/OH_3Dglobal.geos5.47L.4x5.nc OH           1985/1-12/1/0 C xyz kg/m3 * 2 1 1
 )))GLOBAL_OH
 
 # --- Global Cl [mol/mol dry air] ---
@@ -946,6 +946,14 @@ ${RUNDIR_GLOBAL_Cl}
 #==============================================================================
 1 NEGATIVE       -1.0 - - - xy 1 1
 
+#==============================================================================
+# --- Perturbation factors ---
+#
+# Add factors to perturb OH, emissions, and other fields here for
+# analytical inversions.
+#==============================================================================
+2 OH_pert_factor  1.0 - - - xy 1 1
+
 #==============================================================================
 # --- Seasonal scaling factors ----
 #==============================================================================