Add KPP Standalone into GEOS-Chem 14.5.1

.gitignore

@@ -24,4 +24,5 @@ build/
 build_*/
 *___.h
 *___.rc
-core.*
+core.*
+*.exe

.gitmodules

@@ -0,0 +1,4 @@
+[submodule "KPP/standalone"]
+	path = KPP/standalone
+	url = https://github.com/geoschem/KPP-Standalone
+	tag = "GC_Interface"

CHANGELOG.md

@@ -11,11 +11,15 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 - Added `State_Diag%Archive_SatDiagnEdgeCount` field
 - Added `State_Diag%Archive_SatDiagnEdge` field
 - Added routine `SatDiagn_or_SatDiagnEdge` in `History/history_utils_mod.F90`
+- Added KPP standalone interface (archives model state to selected locations)
+- Added `https://github/geoschem/KPP-Standalone` as a Git submodule (`GC_Interface` tag)
+- Added comments in `./run/sharedcleanRunDir.sh` describing the `--force` option (i.e. remove files w/o user confirmation)
 
 ### Changed
 - Renamed `Emiss_Carbon_Gases` to `CO2_Production` in `carbon_gases_mod.F90`
 - Updated start date and restart file for CO2 and tagCO simulations for consistency with carbon simulations
 - Allocated `State_Diag%SatDiagnPEDGE` ffield with vertical dimension `State_Grid%NZ+1`
+- Modified `run/GCClassic/cleanRunDir.sh` to skip removing bpch files, as well as now removing `fort.*` and `OutputDir/*.txt` files
 
 ### Fixed
 - Added a fix to skip the call to KPP when only CO2 is defined in the carbon simulation

GeosCore/CMakeLists.txt

@@ -19,6 +19,7 @@ add_library(GeosCore
   STATIC EXCLUDE_FROM_ALL
   aero_drydep.F90
   aerosol_mod.F90
+  aerosol_thermodynamics_mod.F90
   airs_ch4_mod.F90
   calc_met_mod.F90
   carbon_mod.F90
@@ -47,7 +48,7 @@ add_library(GeosCore
   hco_interface_gc_mod.F90
   hco_utilities_gc_mod.F90
   input_mod.F90
-  aerosol_thermodynamics_mod.F90
+  kppsa_interface_mod.F90
   land_mercury_mod.F90
   linear_chem_mod.F90
   linoz_mod.F90

GeosCore/fullchem_mod.F90

@@ -120,6 +120,7 @@ SUBROUTINE Do_FullChem( Input_Opt,  State_Chm, State_Diag,                 &
     USE GcKpp_Rates,              ONLY : UPDATE_RCONST, RCONST
     USE GcKpp_Util,               ONLY : Get_OHreactivity
     USE Input_Opt_Mod,            ONLY : OptInput
+    USE KppSa_Interface_Mod
     USE Photolysis_Mod,           ONLY : Do_Photolysis, PhotRate_Adj
     USE PhysConstants,            ONLY : AVO, AIRMW
     USE PRESSURE_MOD
@@ -178,7 +179,7 @@ SUBROUTINE Do_FullChem( Input_Opt,  State_Chm, State_Diag,                 &
     INTEGER                :: errorCount, previous_units
     REAL(fp)               :: SO4_FRAC,   T,         TIN
     REAL(fp)               :: TOUT,       SR,        LWC
-
+    REAL(dp)               :: KPPH_before_integrate
     ! Strings
     CHARACTER(LEN=255)     :: errMsg,     thisLoc
 
@@ -200,6 +201,7 @@ SUBROUTINE Do_FullChem( Input_Opt,  State_Chm, State_Diag,                 &
     REAL(dp)               :: RCNTRL (20)
     REAL(dp)               :: RSTATE (20)
     REAL(dp)               :: C_before_integrate(NSPEC)
+    REAL(dp)               :: local_RCONST(NREACT)
 
     ! For tagged CO saving
     REAL(fp)               :: LCH4, PCO_TOT, PCO_CH4, PCO_NMVOC
@@ -235,6 +237,8 @@ SUBROUTINE Do_FullChem( Input_Opt,  State_Chm, State_Diag,                 &
     ! (assuming Rosenbrock solver).  Define this locally in order to break
     ! a compile-time dependency.  -- Bob Yantosca (05 May 2022)
     INTEGER,     PARAMETER :: Nhnew = 3
+    ! Add Nhexit, the last timestep length -- Obin Sturm (30 April 2024)
+    INTEGER,     PARAMETER :: Nhexit = 2
 
     ! Suppress printing out KPP error messages after this many errors occur
     INTEGER,     PARAMETER :: INTEGRATE_FAIL_TOGGLE = 20
@@ -438,6 +442,30 @@ SUBROUTINE Do_FullChem( Input_Opt,  State_Chm, State_Diag,                 &
        mapData => NULL()
     ENDIF
 
+    !=======================================================================
+    ! Setup for the KPP standalone interface (Obin Sturm, Bob Yantosca)
+    !
+    ! NOTE: These routines return immediately if the KPP standalone
+    ! interface has been disabled (or if the *.yml file is missing.)
+    !=======================================================================
+
+    ! Get the (I,J) grid box indices for active cells that are on this CPU
+    ! so that we can print the full chemical state to text files.
+    !
+    ! For computational efficency, only do this on the first call, as
+    ! this information does not change with time.
+    IF ( FirstChem ) THEN
+       CALL KppSa_Check_Domain( RC )
+       IF ( RC /= GC_SUCCESS ) THEN
+          ErrMsg = 'Error encountered in "Check_Domain"!'
+          CALL GC_Error( ErrMsg, RC, ThisLoc )
+          RETURN
+       ENDIF
+    ENDIF
+
+    ! Are we within the time window for archiving model state?
+    CALL KppSa_Check_Time( RC )
+
     !========================================================================
     ! Set up integration convergence conditions and timesteps
     ! (cf. M. J. Evans)
@@ -514,6 +542,7 @@ SUBROUTINE Do_FullChem( Input_Opt,  State_Chm, State_Diag,                 &
     !$OMP DEFAULT( SHARED                                                   )&
     !$OMP PRIVATE( I,        J,        L,       N                           )&
     !$OMP PRIVATE( ICNTRL,   C_before_integrate                             )&
+    !$OMP PRIVATE( KPPH_before_integrate,       local_RCONST                )&
     !$OMP PRIVATE( SO4_FRAC, IERR,     RCNTRL,  ISTATUS,   RSTATE           )&
     !$OMP PRIVATE( SpcID,    KppID,    F,       P,         Vloc             )&
     !$OMP PRIVATE( Aout,     Thread,   RC,      S,         LCH4             )&
@@ -569,6 +598,11 @@ SUBROUTINE Do_FullChem( Input_Opt,  State_Chm, State_Diag,                 &
        ! atmosphere if keepActive option is enabled. (hplin, 2/9/22)
        CALL fullchem_AR_SetKeepActive( option=.TRUE. )
 
+       ! Check if the current grid cell in this loop should have its
+       ! full chemical state printed (concentrations, rates, constants)
+       ! for use with the KPP Standalone (psturm, 03/22/24)
+       CALL KppSa_Check_ActiveCell( I, J, L )
+
        ! Start measuring KPP-related routine timing for this grid box
        IF ( State_Diag%Archive_KppTime ) THEN
           call cpu_time(TimeStart)
@@ -990,6 +1024,11 @@ SUBROUTINE Do_FullChem( Input_Opt,  State_Chm, State_Diag,                 &
        ! let us reset concentrations before calling "Integrate" a 2nd time.
        C_before_integrate = C
 
+       ! Do the same for the KPP initial timestep
+       ! Save local rate constants too
+       KPPH_before_integrate = State_Chm%KPPHvalue(I,J,L)
+       local_RCONST          = RCONST
+
        ! Call the Rosenbrock integrator
        ! (with optional auto-reduce functionality)
        CALL Integrate( TIN,    TOUT,    ICNTRL,                              &
@@ -1260,6 +1299,26 @@ SUBROUTINE Do_FullChem( Input_Opt,  State_Chm, State_Diag,                 &
           State_Diag%KppTime(I,J,L) = TimeEnd - TimeStart
        ENDIF
 
+       ! Write chemical state to file for the kpp standalone interface
+       ! No external logic needed, this subroutine exits early if the
+       ! chemical state should not be printed (psturm, 03/23/24)
+       CALL KppSa_Write_Samples(                                             &
+            I            = I,                                                &
+            J            = J,                                                &
+            L            = L,                                                &
+            initC        = C_before_integrate,                               &
+            localRCONST  = local_RCONST,                                     &
+            initHvalue   = KPPH_before_integrate,                            &
+            exitHvalue   = RSTATE(Nhexit),                                   &
+            ICNTRL       = ICNTRL,                                           &
+            RCNTRL       = RCNTRL,                                           &
+            State_Grid   = State_Grid,                                       &
+            State_Chm    = State_Chm,                                        &
+            State_Met    = State_Met,                                        &
+            Input_Opt    = Input_Opt,                                        &
+            KPP_TotSteps = ISTATUS(3),                                       &
+            RC           = RC                                               )
+
        !=====================================================================
        ! Check we have no negative values and copy the concentrations
        ! calculated from the C array back into State_Chm%Species%Conc
@@ -2664,6 +2723,7 @@ SUBROUTINE Init_FullChem( Input_Opt, State_Chm, State_Diag, RC )
     USE Gckpp_Parameters,         ONLY : nFam, nReact
     USE Gckpp_Global,             ONLY : Henry_K0, Henry_CR, MW, SR_MW
     USE Input_Opt_Mod,            ONLY : OptInput
+    USE KppSa_Interface_Mod,      ONLY : KppSa_Config
     USE State_Chm_Mod,            ONLY : ChmState
     USE State_Chm_Mod,            ONLY : Ind_
     USE State_Diag_Mod,           ONLY : DgnState
@@ -2693,9 +2753,9 @@ SUBROUTINE Init_FullChem( Input_Opt, State_Chm, State_Diag, RC )
     ! Strings
     CHARACTER(LEN=255) :: ErrMsg,   ThisLoc
 
-    !=======================================================================
+    !========================================================================
     ! Init_FullChem begins here!
-    !=======================================================================
+    !========================================================================
 
     ! Assume success
     RC       = GC_SUCCESS
@@ -2705,9 +2765,9 @@ SUBROUTINE Init_FullChem( Input_Opt, State_Chm, State_Diag, RC )
     ! modify the IF statement accordingly to allow initialization
     IF ( .not. Input_Opt%ITS_A_FULLCHEM_SIM ) RETURN
 
-    !=======================================================================
+    !========================================================================
     ! Initialize variables
-    !=======================================================================
+    !========================================================================
     ErrMsg   = ''
     ThisLoc  = ' -> at Init_FullChem (in module GeosCore/FullChem_mod.F90)'
 
@@ -2823,10 +2883,10 @@ SUBROUTINE Init_FullChem( Input_Opt, State_Chm, State_Diag, RC )
                                State_Diag%Archive_O3PconcAfterChem          )
 
 
-    !=======================================================================
+    !========================================================================
     ! Assign default values for KPP absolute and relative tolerances
     ! for species where these have not been explicitly defined.
-    !=======================================================================
+    !========================================================================
     WHERE( State_Chm%KPP_AbsTol == MISSING_DBLE )
        State_Chm%KPP_AbsTol = 1.0e-2_f8
     ENDWHERE
@@ -2835,10 +2895,10 @@ SUBROUTINE Init_FullChem( Input_Opt, State_Chm, State_Diag, RC )
        State_Chm%KPP_RelTol = 0.5e-2_f8
     ENDWHERE
 
-    !=======================================================================
+    !========================================================================
     ! Save physical parameters from the species database into KPP arrays
     ! in gckpp_Global.F90.  These are for the hetchem routines.
-    !=======================================================================
+    !========================================================================
     DO KppId = 1, State_Chm%nKppSpc + State_Chm%nOmitted
        N                  = State_Chm%Map_KppSpc(KppId)
        IF ( N > 0 ) THEN
@@ -2848,18 +2908,18 @@ SUBROUTINE Init_FullChem( Input_Opt, State_Chm, State_Diag, RC )
           HENRY_CR(KppId) = State_Chm%SpcData(N)%Info%Henry_CR
        ENDIF
     ENDDO
-    !=======================================================================
+    !========================================================================
     ! Allocate arrays
-    !=======================================================================
+    !========================================================================
 
     ! Initialize
     id_PSO4 = -1
     id_PCO  = -1
     id_LCH4 = -1
 
-    !--------------------------------------------------------------------
+    !------------------------------------------------------------------------
     ! Pre-store the KPP indices for each KPP prod/loss species or family
-    !--------------------------------------------------------------------
+    !------------------------------------------------------------------------
 
     IF ( nFam > 0 ) THEN
 
@@ -2900,11 +2960,11 @@ SUBROUTINE Init_FullChem( Input_Opt, State_Chm, State_Diag, RC )
     ENDIF
 
 #ifdef MODEL_CESM
-    !--------------------------------------------------------------------
+    !------------------------------------------------------------------------
     ! If we are finding H2SO4_RATE from a fullchem
     ! simulation for the CESM, throw an error if we cannot find
     ! the PSO4 prod family in this KPP mechanism.
-    !--------------------------------------------------------------------
+    !------------------------------------------------------------------------
     IF ( id_PSO4 < 1 ) THEN
        ErrMsg = 'Could not find PSO4 in list of KPP families!  This   ' // &
                 'is needed for State_Chm%H2SO4_PRDR and coupling to CESM!'
@@ -2913,11 +2973,11 @@ SUBROUTINE Init_FullChem( Input_Opt, State_Chm, State_Diag, RC )
     ENDIF
 #endif
 
-    !--------------------------------------------------------------------
+    !------------------------------------------------------------------------
     ! If we are archiving the P(CO) from CH4 and from NMVOC from a fullchem
     ! simulation for the tagCO simulation, throw an error if we cannot find
     ! the PCO or LCH4 prod/loss families in this KPP mechanism.
-    !--------------------------------------------------------------------
+    !------------------------------------------------------------------------
     IF ( State_Diag%Archive_ProdCOfromCH4    .or.                            &
          State_Diag%Archive_ProdCOfromNMVOC ) THEN
 
@@ -2937,19 +2997,19 @@ SUBROUTINE Init_FullChem( Input_Opt, State_Chm, State_Diag, RC )
 
     ENDIF
 
-    !--------------------------------------------------------------------
+    !------------------------------------------------------------------------
     ! Initialize sulfate chemistry code (cf Mike Long)
-    !--------------------------------------------------------------------
+    !------------------------------------------------------------------------
     CALL fullchem_InitSulfurChem( RC )
     IF ( RC /= GC_SUCCESS ) THEN
        ErrMsg = 'Error encountered in "fullchem_InitSulfurCldChem"!'
        CALL GC_Error( ErrMsg, RC, ThisLoc )
        RETURN
     ENDIF
 
-    !--------------------------------------------------------------------
+    !------------------------------------------------------------------------
     ! Initialize dust acid uptake code (Mike Long, Bob Yantosca)
-    !--------------------------------------------------------------------
+    !------------------------------------------------------------------------
     IF ( Input_Opt%LDSTUP ) THEN
        CALL aciduptake_InitDustChem( RC )
        IF ( RC /= GC_SUCCESS ) THEN
@@ -2959,6 +3019,18 @@ SUBROUTINE Init_FullChem( Input_Opt, State_Chm, State_Diag, RC )
        ENDIF
     ENDIF
 
+    !------------------------------------------------------------------------
+    ! Initialize the KPP standalone interface, which will save model state
+    ! for the grid cells specified in kpp_standalone_interface.yml.
+    ! This is needed for input to the KPP standalone box model.
+    !------------------------------------------------------------------------
+    CALL KppSa_Config( Input_Opt, RC )
+    IF ( RC /= GC_SUCCESS ) THEN
+       ErrMsg = 'Error encountered in "KPP_Standalone"!'
+       CALL GC_Error( ErrMsg, RC, ThisLoc )
+       RETURN
+    ENDIF
+
   END SUBROUTINE Init_FullChem
 !EOC
 !------------------------------------------------------------------------------
@@ -2978,6 +3050,7 @@ SUBROUTINE Cleanup_FullChem( RC )
 ! !USES:
 !
     USE ErrCode_Mod
+    USE KppSa_Interface_Mod, ONLY : KppSa_Cleanup
 !
 ! !OUTPUT PARAMETERS:
 !
@@ -3027,6 +3100,10 @@ SUBROUTINE Cleanup_FullChem( RC )
        IF ( RC /= GC_SUCCESS ) RETURN
     ENDIF
 
+    ! Deallocate variables from kpp standalone module
+    ! psturm, 03/22/2024
+    CALL KppSa_Cleanup( RC )
+
   END SUBROUTINE Cleanup_FullChem
 !EOC
 END MODULE FullChem_Mod