diff --git a/.gitignore b/.gitignore index 9cd76c482..a8d8153ad 100644 --- a/.gitignore +++ b/.gitignore @@ -24,4 +24,5 @@ build/ build_*/ *___.h *___.rc -core.* \ No newline at end of file +core.* +*.exe \ No newline at end of file diff --git a/CHANGELOG.md b/CHANGELOG.md index 05d62b156..bef5b5c1f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,6 +17,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Added computation of water concentration to use in photolysis for application of UV absorption by water in Cloud-J v8 - Added ACO3, ACR, ACRO2, ALK4N{1,2,O}2, ALK4P, ALK7, APAN, APINN, APINO2, APINP, AROCMCHO, AROMCO3, AROMPN, BPINN, BPINO2, BPINON, BPINOO2, BPINOOH, BPINP, BUTN, BUTO2, C4H6, C96N, C96O2, C9602H, EBZ, GCO3, HACTA, LIMAL, LIMKB, LIMKET, LIMKO2, LIMN, LIMNB, LIMO2H, LIMO3, LIMO3H, LIMPAN, MEKCO3, MEKPN, MYRCO, PHAN, PIN, PINAL, PINO3, PINONIC, PINPAN, R7N{1,2}, R7O2, R7P, RNO3, STYR, TLFUO2, TLFUONE, TMB, ZRO2 to `species_database.yml` following Travis et al. 2024. - Added TSOIL1 field to `State_Met` for use in HEMCO soil NOx extension. This should only be read in when the `UseSoilTemperature` option is true in HEMCO config. +- Added KPP standalone interface (archives model state to selected locations) ### Changed - Copy values from `State_Chm%KPP_AbsTol` to `ATOL` and `State_Chm%KPP_RelTol` to `RTOL` for fullchem and Hg simulations @@ -32,6 +33,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), - Read aerosol optical properties files from new data directory specified in geoschem_config.yml rather than directory containing photolysis input files - Call `RD_AOD` and `CALC_AOD` from `Init_Aerosol` rather than `Init_Photolysis` - Moved PINO3H to be in alphabetical order in `species_database.yml` +- Modified `run/GCClassic/cleanRunDir.sh` to skip removing bpch files, as well as now removing `fort.*` and `OutputDir/*.txt` files ### Fixed - Simplified SOA representations and fixed related AOD and TotalOA/OC calculations in benchmark. diff --git a/GeosCore/CMakeLists.txt b/GeosCore/CMakeLists.txt old mode 100755 new mode 100644 index f01fbbd42..ebb7ae2dd --- a/GeosCore/CMakeLists.txt +++ b/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 diff --git a/GeosCore/fullchem_mod.F90 b/GeosCore/fullchem_mod.F90 index ecdda088d..aacee485f 100644 --- a/GeosCore/fullchem_mod.F90 +++ b/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,37 @@ 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 ) + + ! test the force write option on the root node + ! example use case: printing chemical state under conditions + ! without knowing where those conditions will happen + ! IF ( Input_Opt%amIRoot .AND. L == 1 ) & + ! CALL Write_Samples( I, J, L, C_before_integrate, & + ! local_RCONST, KPPH_before_integrate, & + ! RSTATE(Nhexit), & + ! State_Grid, State_Chm, State_Met, & + ! Input_Opt, ISTATUS(3), RC, & + ! FORCE_WRITE = .TRUE., CELL_NAME = 'root') + !===================================================================== ! Check we have no negative values and copy the concentrations ! calculated from the C array back into State_Chm%Species%Conc @@ -2664,6 +2734,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 +2764,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 +2776,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 +2894,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 +2906,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 +2919,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 +2971,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 +2984,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,9 +3008,9 @@ 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"!' @@ -2947,9 +3018,9 @@ SUBROUTINE Init_FullChem( Input_Opt, State_Chm, State_Diag, RC ) 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 +3030,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 +3061,7 @@ SUBROUTINE Cleanup_FullChem( RC ) ! !USES: ! USE ErrCode_Mod + USE KppSa_Interface_Mod, ONLY : KppSa_Cleanup ! ! !OUTPUT PARAMETERS: ! @@ -3027,6 +3111,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 diff --git a/GeosCore/kppsa_interface_mod.F90 b/GeosCore/kppsa_interface_mod.F90 new file mode 100644 index 000000000..226626d85 --- /dev/null +++ b/GeosCore/kppsa_interface_mod.F90 @@ -0,0 +1,1015 @@ +!------------------------------------------------------------------------------ +! GEOS-Chem Global Chemical Transport Model ! +!------------------------------------------------------------------------------ +!BOP +! +! !MODULE: kppsa_interface_mod.F90 +! +! !DESCRIPTION: Contains routines to print the full chemical state +! which can be used as input to the KPP Standalone. +!\\ +!\\ +! !INTERFACE: +! +MODULE KppSa_Interface_Mod +! +! !USES: +! + USE Precision_Mod + USE HCO_Error_Mod, ONLY : hp + + IMPLICIT NONE + PRIVATE +! +! !PUBLIC MEMBERS: +! + PUBLIC :: KppSa_Check_ActiveCell + PUBLIC :: KppSa_Check_Domain + PUBLIC :: KppSa_Check_Time + PUBLIC :: KppSa_Cleanup + PUBLIC :: KppSa_Config + PUBLIC :: KppSa_Write_Samples +! +! !DERIVED TYPES: +! + ! Type to hold information read from the YAML config file + TYPE, PRIVATE :: KppSa_Interface_Type + INTEGER :: NLOC + INTEGER :: Start_Output(2) + INTEGER :: Stop_Output(2) + LOGICAL :: SkipIt + LOGICAL :: SkipWriteAtThisTime + CHARACTER(LEN=255) :: Output_Directory + CHARACTER(LEN=255), ALLOCATABLE :: LocationName(:) + REAL(hp), ALLOCATABLE :: LocationLons(:) + REAL(hp), ALLOCATABLE :: LocationLats(:) + INTEGER, ALLOCATABLE :: IDX(:) + INTEGER, ALLOCATABLE :: JDX(:) + INTEGER, ALLOCATABLE :: Levels(:) + END TYPE KppSa_Interface_Type + + ! Type to denote active cells + TYPE, PRIVATE :: KppSa_ActiveCell_Type + LOGICAL :: Active_Cell + CHARACTER(LEN=255) :: Active_Cell_Name + END TYPE KppSa_ActiveCell_Type +! +! !PRIVATE DATA MEMBERS: +! + TYPE(KppSa_Interface_Type), PRIVATE :: KppSa_State + TYPE(KppSa_ActiveCell_Type), PRIVATE :: KppSa_ActiveCell + !$OMP THREADPRIVATE( KppSa_ActiveCell ) +! +! !AUTHORS: +! P. Obin Sturm (psturm@usc.edu) +! +! !REVISION HISTORY: +! 11 Mar 2024 - P. Obin Sturm - Initial version +! See https://github.com/geoschem/geos-chem for complete history +!EOP +!------------------------------------------------------------------------------ +!BOC +CONTAINS +!EOC +!------------------------------------------------------------------------------ +! GEOS-Chem Global Chemical Transport Model ! +!------------------------------------------------------------------------------ +!BOP +! +! !IROUTINE: kppsa_check_domain +! +! !DESCRIPTION: Subroutine Check_Domain is used to identify if a +! specified latitude and longitude falls within a grid cell on the +! current CPU. Multiple lat/lon pairs can be checked simultaneously. +!\\ +!\\ +! !INTERFACE: +! + SUBROUTINE KppSa_Check_Domain( RC ) +! +! !USES: +! + USE HCO_GeoTools_Mod, ONLY : HCO_GetHorzIJIndex + USE HCO_State_GC_Mod, ONLY : HcoState +! +! !OUTPUT PARAMETERS +! + INTEGER, INTENT(out) :: RC +! +! !REVISION HISTORY: +! 11 Mar 2024 - P. Obin Sturm - Initial version +! See https://github.com/geoschem/geos-chem for complete history +!EOP +!------------------------------------------------------------------------------ +!BOC + + ! Early exit if no locations + IF ( KppSa_State%SkipIt ) RETURN + + ! Compute (I,J) indices of grid boxes + CALL HCO_GetHorzIJIndex( HcoState, & + KppSa_State%NLOC, & + KppSa_State%LocationLons, & + KppSa_State%LocationLats, & + KppSa_State%IDX, & + KppSa_State%JDX, & + RC ) + + END SUBROUTINE KppSa_Check_Domain +!EOC +!------------------------------------------------------------------------------ +! GEOS-Chem Global Chemical Transport Model ! +!------------------------------------------------------------------------------ +!BOP +! +! !IROUTINE: kppsa_check_time +! +! !DESCRIPTION: Subroutine Check_Domain is used to identify if a +! specified latitude and longitude falls within a grid cell on the +! current CPU. Multiple lat/lon pairs can be checked simultaneously. +!\\ +!\\ +! !INTERFACE: +! + SUBROUTINE KppSa_Check_Time( RC ) +! +! !USES: +! + USE ErrCode_Mod + USE Time_Mod, ONLY : Get_Nymd, Get_Nhms +! +! !OUTPUT PARAMETERS +! + INTEGER, INTENT(OUT) :: RC ! Success or failure? +! +! !REVISION HISTORY: +! 11 Mar 2024 - P. Obin Sturm - Initial version +! See https://github.com/geoschem/geos-chem for complete history +!EOP +!------------------------------------------------------------------------------ +!BOC +! +! !LOCAL VARIABLES: +! + INTEGER :: yyyymmdd, hhmmss + + ! Initialize + RC = GC_SUCCESS + + ! Early exit if no locations + IF ( KppSa_State%SkipIt ) RETURN + + ! Assume we will not write to disk at this date/time + KppSa_State%SkipWriteAtThisTime = .TRUE. + + ! Get current date & time + yyyymmdd = Get_Nymd() + hhmmss = Get_Nhms() + + ! Exit if we are outside the window for archiving model state + IF ( yyyymmdd < KppSa_State%Start_Output(1) ) RETURN + IF ( yyyymmdd > KppSa_State%Stop_Output(1) ) RETURN + IF ( hhmmss < KppSa_State%Start_Output(2) ) RETURN + IF ( hhmmss > KppSa_State%Stop_Output(2) ) RETURN + + ! If we get this far, we're in the time window where we + ! archive the chemical state for the KPP standalone + KppSa_State%SkipWriteAtThisTime = .FALSE. + + END SUBROUTINE KppSa_Check_Time +!EOC +!------------------------------------------------------------------------------ +! GEOS-Chem Global Chemical Transport Model ! +!------------------------------------------------------------------------------ +!BOP +! +! !IROUTINE: kppsa_check_activecell +! +! !DESCRIPTION: Identifies if a grid cell is within a specified latitude +! and longitude to print the full chemical state (all concentrations, +! reaction rates, rate constants, and meteo metadata). +!\\ +!\\ +! !INTERFACE: +! + SUBROUTINE KppSa_Check_ActiveCell( I, J, L ) +! +! !INPUT PARAMETERS: +! + INTEGER, INTENT(IN) :: I, J, L ! Grid Indices +! +! !REVISION HISTORY: +! 11 Mar 2024 - P. Obin Sturm - Initial version +! See https://github.com/geoschem/geos-chem for complete history +!EOP +!------------------------------------------------------------------------------ +!BOC +! +! !LOCAL VARIABLES +! + INTEGER :: K + + ! Early exit if KPP standalone interface is disabled + IF ( KppSa_State%SkipIt ) RETURN + + ! Initialize + KppSa_ActiveCell%Active_Cell = .FALSE. + KppSa_ActiveCell%Active_Cell_Name = '' + + ! Skip if we are outside the time interval + IF ( KppSa_State%SkipWriteAtThisTime ) RETURN + + ! Flag active cells + IF ( ANY( L == KppSa_State%Levels ) ) THEN + DO K = 1, KppSa_State%NLOC + IF ( KppSa_State%IDX(K) == I .AND. & + KppSa_State%JDX(K) == J ) THEN + KppSa_ActiveCell%Active_Cell = .TRUE. + KppSa_ActiveCell%Active_Cell_Name = & + KppSa_State%LocationName(K) + ENDIF + ENDDO + ENDIF + + END SUBROUTINE KppSa_Check_ActiveCell +!EOC +!------------------------------------------------------------------------------ +! GEOS-Chem Global Chemical Transport Model ! +!------------------------------------------------------------------------------ +!BOP +! +! !IROUTINE: kppsa_config +! +! !DESCRIPTION: Subroutine Config_KPP_Standalone reads a set of gridcells +! to be sampled and the full chemical state printed. +!\\ +!\\ +! !INTERFACE: +! + SUBROUTINE KppSa_Config( Input_Opt, RC ) +! +! !USES: +! + USE QfYaml_Mod + USE ErrCode_Mod + USE Input_Opt_Mod, ONLY : OptInput + USE RoundOff_Mod, ONLY : Cast_and_RoundOff + USE inquireMod, ONLY : findFreeLUN +! +! !INPUT PARAMETERS: +! + TYPE(OptInput), INTENT(IN) :: Input_Opt ! Input Options object +! +! !OUTPUT PARAMETERS: +! + INTEGER, INTENT(OUT) :: RC ! Success or failure +! +! !REVISION HISTORY: +! 11 Mar 2024 - P. Obin Sturm - Initial version +! See https://github.com/geoschem/geos-chem for complete history +!EOP +!------------------------------------------------------------------------------ +!BOC +! +! !LOCAL VARIABLES: +! + ! Scalars + INTEGER :: I, N + INTEGER :: IU_FILE ! Available unit for writing + INTEGER :: path_exists + LOGICAL :: file_exists + LOGICAL :: v_bool + + ! Strings + CHARACTER(LEN=255) :: thisLoc + CHARACTER(LEN=512) :: errMsg + CHARACTER(LEN=QFYAML_NamLen) :: key + CHARACTER(LEN=QFYAML_StrLen) :: v_str + + ! Objects + TYPE(QFYAML_t) :: Config, ConfigAnchored + + ! Arrays + INTEGER :: a_int(QFYAML_MaxArr) + + ! String arrays + CHARACTER(LEN=QFYAML_NamLen) :: a_str(QFYAML_MaxArr) + + ! YAML configuration file name to be read + CHARACTER(LEN=30), PARAMETER :: configFile = & + './kpp_standalone_interface.yml' + + ! Inquire if YAML interface exists -- if not, skip initializing + KppSa_State%SkipIt = .FALSE. + INQUIRE( FILE=configFile, EXIST=file_exists ) + IF ( .NOT. file_exists ) THEN + KppSa_State%SkipIt = .TRUE. + IF ( Input_Opt%amIRoot ) THEN + WRITE( 6, 100 ) TRIM( configFile ) + 100 FORMAT( "Config file ", a ", not found, ", & + "skipping KPP standalone interface" ) + RETURN + ENDIF + ENDIF + + ! Assume success + RC = GC_SUCCESS + errMsg = '' + thisLoc = ' -> at Config_KPP_Standalone (in module GeosCore/kpp_standalone_interface.F90)' + + !======================================================================== + ! Read the YAML file into the Config object + !======================================================================== + CALL QFYAML_Init( configFile, Config, ConfigAnchored, RC ) + IF ( RC /= GC_SUCCESS ) THEN + errMsg = 'Error reading configuration file: ' // TRIM( configFile ) + CALL GC_Error( errMsg, RC, thisLoc ) + RETURN + ENDIF + + !======================================================================== + ! Read the main on/off switch; Exit if the switch is turned off + !======================================================================== + key = "settings%activate" + 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 + KppSa_State%SkipIt = ( .not. v_bool ) + IF ( KppSa_State%SkipIt ) THEN + WRITE( 6, 110 ) + 110 FORMAT( "KPP standalone interface was manually disabled" ) + RETURN + ENDIF + + !======================================================================== + ! Read the list of active cells + !======================================================================== + key = "active_cells" + a_str = MISSING_STR + CALL QFYAML_Add_Get( Config, key, a_str, "", RC, dynamic_size=.TRUE. ) + IF ( RC /= GC_SUCCESS ) THEN + errMsg = 'Error parsing ' // TRIM( key ) // '!' + CALL GC_Error( errMsg, RC, thisLoc ) + RETURN + ENDIF + + !======================================================================== + ! Get the number of active cells (if 0, return) and the list of names + !======================================================================== + KppSa_State%NLOC = Find_Number_of_Locations( a_str ) + IF ( KppSa_State%NLOC .eq. 0 ) THEN + ! Set SkipIt flag to short circuit other subroutines + KppSa_State%SkipIt = .TRUE. + IF ( Input_Opt%amIRoot ) THEN + WRITE( 6, 120 ) + 120 FORMAT( "No active cells for box modeling ", & + "in kpp_standalone_interface.yml") + RETURN + ENDIF + ENDIF + ALLOCATE( KppSa_State%LocationName( KppSa_State%NLOC ), STAT=RC ) + CALL GC_CheckVar( 'KppSa_State%LocationName', 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + DO I = 1,KppSa_State%NLOC + KppSa_State%LocationName(I) = TRIM( a_str(I) ) + END DO + + !======================================================================== + ! Read latitude and longitude of active cells + !======================================================================== + + ! Allocate number of locations for lats and lons + ALLOCATE( KppSa_State%LocationLons( KppSa_State%NLOC ), STAT=RC ) + CALL GC_CheckVar( 'KppSa_State%LocationLons', 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + + ALLOCATE( KppSa_State%LocationLats( KppSa_State%NLOC ), STAT=RC ) + CALL GC_CheckVar( 'KppSa_State%LocationLats', 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + + ! Read coordinates + DO I = 1,KppSa_State%NLOC + ! Read longitudes + key = "locations%"//TRIM( KppSa_State%LocationName(I) )//"%longitude" + 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 + KppSa_State%LocationLons( I ) = Cast_and_RoundOff( TRIM( v_str ), places=-1 ) + ! Read latitudes + key = "locations%"//TRIM( KppSa_State%LocationName(I) )//"%latitude" + 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 + KppSa_State%LocationLats( I ) = Cast_and_RoundOff( TRIM( v_str ), places=-1 ) + END DO + + ! Allocate IDX and JDX (masks for whether a location is on the CPU) + ALLOCATE( KppSa_State%IDX( KppSa_State%NLOC ), STAT=RC ) + CALL GC_CheckVar( 'KppSa_State%IDX', 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + + ALLOCATE( KppSa_State%JDX( KppSa_State%NLOC ), STAT=RC ) + CALL GC_CheckVar( 'KppSa_State%JDX', 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + + KppSa_State%IDX(:) = -1 + KppSa_State%JDX(:) = -1 + + !======================================================================== + ! Get the list of levels and number of levels + !======================================================================== + ! TODO: could add capability for location specific levels + key = "settings%levels" + a_int = MISSING_INT + CALL QFYAML_Add_Get( Config, key, a_int, "", RC, dynamic_size=.TRUE. ) + IF ( RC /= GC_SUCCESS ) THEN + errMsg = 'Error parsing ' // TRIM( key ) // '!' + CALL GC_Error( errMsg, RC, thisLoc ) + RETURN + ENDIF + N = Find_Number_of_Levels( a_int ) + ! if no specified levels, print the surface + IF ( N .eq. 0 ) THEN + N = 1 + a_int(1) = 1 + END IF + ALLOCATE( KppSa_State%Levels( N ), STAT=RC ) + CALL GC_CheckVar( 'KppSa_State%Levels', 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + DO I = 1,N + KppSa_State%Levels(I) = a_int(I) + END DO + + !======================================================================== + ! Get the start & stop date/time for which output will be printed + !======================================================================== + key = "settings%start_output_at" + a_int = MISSING_INT + CALL QFYAML_Add_Get( Config, key, a_int(1:2), "", RC ) + IF ( RC /= GC_SUCCESS ) THEN + errMsg = 'Error parsing ' // TRIM( key ) // '!' + CALL GC_Error( errMsg, RC, thisLoc ) + RETURN + ENDIF + KppSa_State%Start_Output = a_int(1:2) + + key = "settings%stop_output_at" + a_int = MISSING_INT + CALL QFYAML_Add_Get( Config, key, a_int(1:2), "", RC ) + IF ( RC /= GC_SUCCESS ) THEN + errMsg = 'Error parsing ' // TRIM( key ) // '!' + CALL GC_Error( errMsg, RC, thisLoc ) + RETURN + ENDIF + KppSa_State%Stop_Output = a_int(1:2) + + !======================================================================== + ! Set the output directory + !======================================================================== + ! Get that value + key = "settings%output_directory" + 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 + ! Check to see if the directory exists + ! Do this in a portable way that works across compilers + ! The directory specifier in inquire might be specific to ifort + ! So instead try to open a test file within the output directory + ! Check ./OutputDir (which exists for GEOS-Chem and GCHP) as backup + IU_FILE = findFreeLUN() + OPEN( IU_FILE, FILE = trim(v_str)//'/.test_directory_existence', & + ACTION = "WRITE", & + IOSTAT = path_exists, & + ACCESS = 'SEQUENTIAL' ) + + ! If the specified folder doesn't exist, try OutputDir + IF ( path_exists /= 0 ) THEN + OPEN( IU_FILE, FILE = './OutputDir/.test_directory_existence', & + ACTION = "WRITE", & + IOSTAT = path_exists, & + ACCESS ='SEQUENTIAL' ) + KppSa_State%Output_Directory = "./OutputDir" + IF ( Input_Opt%amIRoot ) THEN + WRITE( 6, '(a)' ) & + "KPP Standalone Interface warning: Specified output directory ",& + trim(v_str), & + " was not found, trying default output path './OutputDir' " + ENDIF + + ! If OutputDir doesn't exist, write to the current directory + IF ( path_exists /= 0 ) THEN + IF ( Input_Opt%amIRoot ) THEN + WRITE( 6, '(a)' ) & + "KPP Standalone Interface warning: Specified output directory ", & + trim(v_str), & + " and default output directory './OutputDir' " // & + "were not found, writing output to the current directory './'" + KppSa_State%Output_Directory = "./" + ENDIF + ENDIF + ELSE + KppSa_State%Output_Directory = trim(v_str) + close(IU_FILE) + END IF + + !======================================================================= + ! Print information about sites that will be archived + !======================================================================= + IF ( Input_Opt%amIRoot ) THEN + WRITE( 6, '(a)' ) REPEAT( "=", 79 ) + WRITE( 6, '(a,/)' ) "KPP STANDALONE INTERFACE" + WRITE( 6, '(a,/)' ) "Model state will be archived at these sites:" + DO I = 1, KppSa_State%NLOC + WRITE( 6, 150 ) KppSa_State%LocationName(I), & + KppSa_State%LocationLons(I), & + KppSa_State%LocationLats(I) + 150 FORMAT( a25, "( ", f9.4, ", ", f9.4, " )") + ENDDO + WRITE( 6, '(/,a)' ) "For GEOS-Chem vertical levels:" + WRITE( 6, '(100i4)' ) KppSa_State%Levels + WRITE( 6, 160 ) KppSa_State%Start_Output + 160 FORMAT( "Starting at ", i8.8, 1x, i6.6 ) + WRITE( 6, 170 ) KppSa_State%Stop_Output + 170 FORMAT( "Ending at ", i8.8, 1x, i6.6 ) + WRITE( 6, '(a)' ) REPEAT( "=", 79 ) + ENDIF + + END SUBROUTINE KppSa_Config +!EOC +!------------------------------------------------------------------------------ +! GEOS-Chem Global Chemical Transport Model ! +!------------------------------------------------------------------------------ +!BOP +! +! !IROUTINE: kppsa_write_samples +! +! !DESCRIPTION: Subroutine Write_Samples writes the full chemical state +! (concentrations, reaction rates and rate constants, meteorological +! conditions). +!\\ +!\\ +! !INTERFACE: +! + SUBROUTINE KppSa_Write_Samples( I, J, L, & + initC, localRCONST, initHvalue, & + exitHvalue, ICNTRL, RCNTRL, & + State_Grid, State_Chm, State_Met, & + Input_Opt, KPP_TotSteps, RC, & + FORCE_WRITE, CELL_NAME ) +! +! !USES: +! + USE ErrCode_Mod + USE State_Grid_Mod, ONLY : GrdState + USE State_Chm_Mod, ONLY : ChmState + USE State_Met_Mod, ONLY : MetState + USE Input_Opt_Mod, ONLY : OptInput + USE GcKpp_Global, ONLY : ATOL + USE GcKpp_Function + USE GcKpp_Parameters, ONLY : NSPEC, NREACT, NVAR + USE TIME_MOD, ONLY : GET_TS_CHEM + USE TIME_MOD, ONLY : TIMESTAMP_STRING + USE TIME_MOD, ONLY : Get_Minute + USE TIME_MOD, ONLY : Get_Hour + USE TIME_MOD, ONLY : Get_Day + USE TIME_MOD, ONLY : Get_Month + USE TIME_MOD, ONLY : Get_Year + USE Pressure_Mod, ONLY : Get_Pcenter + USE inquireMod, ONLY : findFreeLUN +! +! !INPUT PARAMETERS: +! + INTEGER, INTENT(IN) :: I ! Longitude index + INTEGER, INTENT(IN) :: J ! Latitude index + INTEGER, INTENT(IN) :: L ! Vertical level + INTEGER, INTENT(IN) :: KPP_TotSteps ! Total integr. steps + INTEGER, INTENT(IN) :: ICNTRL(20) ! Integrator options + TYPE(GrdState), INTENT(IN) :: State_Grid ! Grid State object + TYPE(ChmState), INTENT(IN) :: State_Chm ! Chem State obj + TYPE(MetState), INTENT(IN) :: State_Met ! Met State obj + TYPE(OptInput), INTENT(IN) :: Input_Opt ! Input Options obj + REAL(dp), INTENT(IN) :: initC(NSPEC) ! Initial conc. + REAL(dp), INTENT(IN) :: localRCONST(NREACT) ! Rate constants + REAL(dp) :: initHvalue ! Initial timestep + REAL(dp) :: exitHvalue ! Final timestep: + ! RSTATE(Nhexit) + REAL(dp), INTENT(IN) :: RCNTRL(20) ! Integrator options + LOGICAL, OPTIONAL :: FORCE_WRITE ! Write even if not + ! in an active cell + CHARACTER(LEN=*), OPTIONAL :: CELL_NAME ! Customize name of + ! this file +! +! !OUTPUT PARAMETERS: +! + INTEGER, INTENT(OUT) :: RC ! Success or failure +! +! !REVISION HISTORY: +! 11 Mar 2024 - P. Obin Sturm - Initial version +! See https://github.com/geoschem/geos-chem for complete history +!EOP +!------------------------------------------------------------------------------ +!BOC +! +! !LOCAL VARIABLES: +! + ! Integers + INTEGER :: N + INTEGER :: IU_FILE + INTEGER :: SpcID + REAL(fp) :: DT + LOGICAL :: FORCE_WRITE_AUX + CHARACTER(LEN=255) :: CELL_NAME_AUX + + ! Strings + CHARACTER(LEN=255) :: YYYYMMDD_hhmmz + CHARACTER(LEN=255) :: level_string + CHARACTER(LEN=512) :: errMsg, filename + + ! Arrays + REAL(dp) :: Aout(NREACT) + REAL(dp) :: Vloc(NVAR) + + !====================================================================== + ! Write_Samples begins here! + !====================================================================== + + ! Did a user want to write the chemical state + ! even if not in an active cell? + FORCE_WRITE_AUX = .FALSE. + IF ( PRESENT( FORCE_WRITE ) ) FORCE_WRITE_AUX = FORCE_WRITE + + ! Quit early if there's no writing to be done + IF ( .not. KppSa_ActiveCell%Active_Cell .AND. & + .not. FORCE_WRITE_AUX ) THEN + RETURN + END IF + + ! Did the call include an optional cell name? + CELL_NAME_AUX = '' + IF ( PRESENT( CELL_NAME ) ) CELL_NAME_AUX = CELL_NAME + + ! Get KPP state + CALL Fun( V = initC(1:NVAR), & + F = initC(NVAR+1:NSPEC), & + RCT = localRCONST, & + Vdot = Vloc, & + Aout = Aout ) + + ! Chemistry timestep (seconds) + DT = GET_TS_CHEM() + + !====================================================================== + ! Write the file. We need to place this into an !$OMP CRITICAL + ! block to ensure that only one thread can open & write to the file + ! at a time. Otherwise we will get corrupted files + !====================================================================== + !$OMP CRITICAL + + ! Find a free file LUN + IU_FILE = findFreeLUN() + WRITE(level_string,'(I0)') L + WRITE( YYYYMMDD_hhmmz,'(I0.4,I0.2,I0.2,a,I0.2,I0.2)' ) & + Get_Year(), Get_Month(), Get_Day(), '_', Get_Hour(), Get_Minute() + + ! Filename for output + filename = TRIM( KppSa_State%Output_Directory ) // & + '/' // & + TRIM( Cell_Name_Aux ) // & + TRIM( KppSa_ActiveCell%Active_Cell_Name ) // & + '_L' // & + trim( level_string ) // & + '_' // & + TRIM( YYYYMMDD_hhmmz ) // & + '.txt' + + ! Open the file + ! NOTE: We cannot exit from within an !$OMP CRITICAL block + OPEN( IU_FILE, FILE=TRIM(filename), ACTION="WRITE", & + IOSTAT=RC, ACCESS='SEQUENTIAL') + + ! Write header to file + WRITE( IU_FILE, '(a)' ) '60' + WRITE( IU_FILE, '(a)' ) REPEAT("=", 76 ) + WRITE( IU_FILE, '(a)' ) '' + WRITE( IU_FILE, '(a)' ) & + ' KPP Standalone Atmospheric Chemical State' + WRITE( IU_FILE, '(a)' ) 'File Description:' + WRITE( IU_FILE, '(a)' ) & + 'This file contains model output of the atmospheric chemical state' + WRITE( IU_FILE, '(a)' ) & + 'as simulated by the GEOS-Chem chemistry module in a 3D setting.' + WRITE( IU_FILE, '(a)' ) & + 'Each grid cell represents the chemical state of an individual location,' + WRITE( IU_FILE, '(a)' ) & + 'suitable for input into a separate KPP Standalone program which will' + WRITE( IU_FILE, '(a)' ) & + 'replicate the chemical evolution of that grid cell for mechanism analysis.' + WRITE( IU_FILE, '(a)' ) & + 'Note that the KPP Standalone will only use concentrations, rate constants,' + WRITE( IU_FILE, '(a)' ) & + 'and KPP-specific fields. All other fields are for reference. The first line' + WRITE( IU_FILE, '(a)' ) & + 'contains the number of lines in this header. If wanting to use this output' + WRITE( IU_FILE, '(a)' ) & + 'for other analysis, a Python class to read these fields is available by' + WRITE( IU_FILE, '(a)' ) & + 'request, contact Obin Sturm (psturm@usc.edu).' + WRITE( IU_FILE, '(a)' ) '' + WRITE( IU_FILE, '(a)' ) 'Generated by the GEOS-Chem Model' + WRITE( IU_FILE, '(a)' ) ' (https://geos-chem.org/)' + WRITE( IU_FILE, '(a)' ) 'Using the KPP Standalone Interface' + WRITE( IU_FILE, '(a)' ) 'github.com/GEOS-ESM/geos-chem/tree/feature/psturm/kpp_standalone_interface' + WRITE( IU_FILE, '(a)' ) ' With contributions from:' + WRITE( IU_FILE, '(a)' ) ' Obin Sturm (psturm@usc.edu)' + WRITE( IU_FILE, '(a)' ) ' Christoph Keller' + WRITE( IU_FILE, '(a)' ) ' Michael Long' + WRITE( IU_FILE, '(a)' ) ' Sam Silva' + WRITE( IU_FILE, '(a)' ) '' + + ! Write the grid cell metadata as part of the header + WRITE( IU_FILE, '(a,/)' ) & + 'Meteorological and general grid cell metadata ' + WRITE( IU_FILE, '(a,a)' ) & + 'Location: ' // & + TRIM( CELL_NAME_AUX ) // & + TRIM( KppSa_ActiveCell%ACTIVE_CELL_NAME ) + WRITE( IU_FILE, '(a,a)' ) & + 'Timestamp: ', & + TIMESTAMP_STRING() + WRITE( IU_FILE, '(a,f11.4)' ) & + 'Longitude (degrees): ', & + State_Grid%XMid(I,J) + WRITE( IU_FILE, '(a,f11.4)' ) & + 'Latitude (degrees): ', & + State_Grid%YMid(I,J) + WRITE( IU_FILE, '(a,i6)' ) & + 'GEOS-Chem Vertical Level: ', & + L + WRITE( IU_FILE, '(a,f11.4)' ) & + 'Pressure (hPa): ', & + Get_Pcenter( I, J, L ) + WRITE( IU_FILE, '(a,f11.2)' ) & + 'Temperature (K): ', & + State_Met%T(I,J,L) + WRITE( IU_FILE, '(a,e11.4)' ) & + 'Dry air density (molec/cm3): ', & + State_Met%AIRNUMDEN(I,J,L) + WRITE( IU_FILE, '(a,e11.4)' ) & + 'Water vapor mixing ratio (vol H2O/vol dry air): ', & + State_Met%AVGW(I,J,L) + WRITE( IU_FILE, '(a,e11.4)' ) & + 'Cloud fraction: ', & + State_Met%CLDF(I,J,L) + WRITE( IU_FILE, '(a,e11.4)' ) & + 'Cosine of solar zenith angle: ', & + State_Met%SUNCOSmid(I,J) + WRITE( IU_FILE, '(/,a,/)' ) & + 'KPP Integrator-specific parameters ' + WRITE( IU_FILE, '(a,f11.4)' ) & + 'Init KPP Timestep (seconds): ', & + initHvalue + WRITE( IU_FILE, '(a,f11.4)' ) & + 'Exit KPP Timestep (seconds): ', & + exitHvalue + WRITE( IU_FILE, '(a,f11.4)' ) & + 'Chemistry operator timestep (seconds): ', & + DT + WRITE( IU_FILE, '(a,i6)' ) & + 'Number of internal timesteps: ', & + KPP_TotSteps + WRITE( IU_FILE, '(a)' ) 'ICNTRL integrator options used:' + WRITE( IU_FILE, '(10i6)' ) ICNTRL( 1:10) + WRITE( IU_FILE, '(10i6)' ) ICNTRL(11:20) + WRITE( IU_FILE, '(a)' ) 'RCNTRL integrator options used:' + WRITE( IU_FILE, '(5F13.6)' ) RCNTRL( 1: 5) + WRITE( IU_FILE, '(5F13.6)' ) RCNTRL( 6:10) + WRITE( IU_FILE, '(5F13.6)' ) RCNTRL(11:15) + WRITE( IU_FILE, '(5F13.6)' ) RCNTRL(16:20) + WRITE( IU_FILE, '(/,a)' ) & + 'CSV data of full chemical state, including species concentrations,' + WRITE( IU_FILE, '(a)' ) & + 'rate constants (R) and instantaneous reaction rates (A).' + WRITE( IU_FILE, '(a)' ) & + 'All concentration units are in molec/cm3 and rates in molec/cm3/s.' + WRITE( IU_FILE, '(a)' ) '' + WRITE( IU_FILE, '(a)' ) REPEAT("=", 76 ) + WRITE( IU_FILE, '(a)' ) 'Name, Value, Absolute Tolerance' + + ! Write species concentrations and absolute tolerances + DO N = 1, NSPEC + SpcID = State_Chm%Map_KppSpc(N) + IF ( SpcID <= 0 ) THEN + WRITE( IU_FILE, 120 ) N, initC(N), ATOL(N) + 120 FORMAT( "C", i0, ",", es25.16e3, ",", es10.2e2 ) + CYCLE + ENDIF + WRITE( IU_FILE, 130 ) & + TRIM(State_Chm%SpcData(SpcID)%Info%Name), initC(N), ATOL(N) + 130 FORMAT( a, ",", es25.16e3, ",", es10.2e2 ) + ENDDO + + ! Write reaction rates + DO N = 1, NREACT + WRITE( IU_FILE, 140 ) N, localRCONST(N) + 140 FORMAT( "R", i0, ",", es25.16e3 ) + ENDDO + + ! Write instantaneous reaction rates + DO N = 1, NREACT + WRITE( IU_FILE, 150 ) N, Aout(N) + 150 FORMAT( "A", i0, ",", es25.16e3 ) + ENDDO + + ! Close file + CLOSE( IU_FILE ) + !$OMP END CRITICAL + + END SUBROUTINE KppSa_Write_Samples +!EOC +!------------------------------------------------------------------------------ +! GEOS-Chem Global Chemical Transport Model ! +!------------------------------------------------------------------------------ +!BOP +! +! !IROUTINE: kppsa_cleanup +! +! !DESCRIPTION: Deallocates module variables that may have been allocated +! at run time and unnecessary files required during the process +!\\ +!\\ +! !INTERFACE: +! + SUBROUTINE KppSa_Cleanup( RC ) +! +! !USES: +! + USE ErrCode_Mod + USE inquireMod, ONLY : findFreeLUN +! +! !OUTPUT PARAMETERS: +! + INTEGER, INTENT(OUT) :: RC ! Success or failure? +! +! !REVISION HISTORY: +! 11 Mar 2024 - P. Obin Sturm - Initial version +!EOP +!------------------------------------------------------------------------------ +!BOC +! +! !LOCAL VARIABLES:: +! + ! Strings + CHARACTER(LEN=255) :: arrayId + + ! Assume success + RC = GC_SUCCESS + + IF ( ALLOCATED( KppSa_State%LocationName ) ) THEN + arrayId = 'kpp_standalone_interface.F90:KppSa_State%LocationName' + DEALLOCATE( KppSa_State%LocationName, STAT=RC ) + CALL GC_CheckVar( arrayId, 2, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + ENDIF + + IF ( ALLOCATED( KppSa_State%LocationLons ) ) THEN + arrayId = 'kpp_standalone_interface.F90:KppSa_State%LocationLons' + DEALLOCATE( KppSa_State%LocationLons, STAT=RC ) + CALL GC_CheckVar( arrayId, 2, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + ENDIF + + IF ( ALLOCATED( KppSa_State%LocationLats ) ) THEN + arrayId = 'kpp_standalone_interface.F90:KppSa_State%LocationLats' + DEALLOCATE( KppSa_State%LocationLats, STAT=RC ) + CALL GC_CheckVar( arrayId, 2, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + ENDIF + + IF ( ALLOCATED( KppSa_State%IDX ) ) THEN + arrayId = 'kpp_standalone_interface.F90:KppSa_State%IDX' + DEALLOCATE( KppSa_State%IDX, STAT=RC ) + CALL GC_CheckVar( arrayId, 2, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + ENDIF + + IF ( ALLOCATED( KppSa_State%JDX ) ) THEN + arrayId = 'kpp_standalone_interface.F90:KppSa_State%JDX' + DEALLOCATE( KppSa_State%JDX, STAT=RC ) + CALL GC_CheckVar( arrayId, 2, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + ENDIF + + IF ( ALLOCATED( KppSa_State%Levels ) ) THEN + arrayId = 'kpp_standalone_interface.F90:KppSa_State%Levels' + DEALLOCATE( KppSa_State%Levels, STAT=RC ) + CALL GC_CheckVar( arrayId, 2, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + ENDIF + + END SUBROUTINE KppSa_Cleanup +!EOC +!------------------------------------------------------------------------------ +! GEOS-Chem Global Chemical Transport Model ! +!------------------------------------------------------------------------------ +!BOP +! +! !IROUTINE: Find_Number_of_Locations +! +! !DESCRIPTION: Searches a string array containing location names and returns +! the number of valid locations (i.e. char that do not match MISSING_STR). +! Assumes all the valid locations will be listed contiguously at the front +! of the array. Taken from Get_Number_of_Species from input_mod.F90 +!\\ +!\\ +! !INTERFACE: + FUNCTION Find_Number_of_Locations( a_str ) RESULT( n_valid ) +! +! !INPUT PARAMETERS: +! + CHARACTER(LEN=*), INTENT(IN) :: a_str(:) +! +! !RETURN VALUE: +! + INTEGER :: n_valid +! +! !REVISION HISTORY: +! 11 Mar 2024 - P. Obin Sturm - Initial version +! See https://github.com/geoschem/geos-chem for complete history +!EOP +!------------------------------------------------------------------------------ +!BOC +! +! !LOCAL VARIABLES: +! + INTEGER :: N + + ! Return the number of valid locations + n_valid = 0 + DO N = 1, SIZE( a_str ) + IF ( TRIM( a_str(N) ) == MISSING_STR ) EXIT + n_valid = n_valid + 1 + ENDDO + + END FUNCTION Find_Number_of_Locations +!EOC +!------------------------------------------------------------------------------ +! GEOS-Chem Global Chemical Transport Model ! +!------------------------------------------------------------------------------ +!BOP +! +! !IROUTINE: Find_Number_of_Levels +! +! !DESCRIPTION: Searches an integer array containing location names and returns +! the number of valid levels (i.e. int that do not match MISSING_INT). +! Assumes all the valid levels will be listed contiguously at the front +! of the array. Taken from Get_Number_of_Species from input_mod.F90 +!\\ +!\\ +! !INTERFACE: + FUNCTION Find_Number_of_Levels( a_int ) RESULT( n_valid ) +! +! !INPUT PARAMETERS: +! + INTEGER, INTENT(IN) :: a_int(:) +! +! !RETURN VALUE: +! + INTEGER :: n_valid +! +! !REVISION HISTORY: +! 11 Mar 2024 - P. Obin Sturm - Initial version +! See https://github.com/geoschem/geos-chem for complete history +!EOP +!------------------------------------------------------------------------------ +!BOC +! +! !LOCAL VARIABLES: +! + INTEGER :: N + + ! Return the number of valid locations + n_valid = 0 + DO N = 1, SIZE( a_int ) + IF ( a_int(N) == MISSING_INT ) EXIT + n_valid = n_valid + 1 + ENDDO + + END FUNCTION Find_Number_of_Levels +!EOC +END MODULE KppSa_Interface_Mod diff --git a/run/GCClassic/createRunDir.sh b/run/GCClassic/createRunDir.sh index 39d34e84b..638eb8825 100755 --- a/run/GCClassic/createRunDir.sh +++ b/run/GCClassic/createRunDir.sh @@ -876,11 +876,18 @@ if [[ ${met} = "ModelE2.1" ]] || [[ ${met} = "ModelE2.2" ]]; then cp ${gcdir}/run/shared/download_data.gcap2.40L.yml ${rundir}/download_data.yml fi +# Copy the OH metrics Python script to the rundir (fullchem/CH4 only) if [[ "x${sim_name}" == "xfullchem" || "x${sim_name}" == "xCH4" ]]; then cp -r ${gcdir}/run/shared/metrics.py ${rundir} chmod 744 ${rundir}/metrics.py fi +# Copy the KPP standalone interface config file to ther rundir (fullchem only) +if [[ "x${sim_name}" == "xfullchem" ]]; then + cp -r ${gcdir}/run/shared/kpp_standalone_interface.yml ${rundir} + chmod 644 ${rundir}/kpp_standalone_interface.yml +fi + # Set permissions chmod 744 ${rundir}/cleanRunDir.sh chmod 744 ${rundir}/archiveRun.sh diff --git a/run/GCHP/createRunDir.sh b/run/GCHP/createRunDir.sh index 095d59385..6a14c52db 100755 --- a/run/GCHP/createRunDir.sh +++ b/run/GCHP/createRunDir.sh @@ -622,6 +622,12 @@ if [[ "x${sim_name}" == "xfullchem" || "x${sim_name}" == "xcarbon" ]]; then chmod 744 ${rundir}/metrics.py fi +# Copy the KPP standalone interface config file to ther rundir (fullchem only) +if [[ "x${sim_name}" == "xfullchem" ]]; then + cp -r ${gcdir}/run/shared/kpp_standalone_interface.yml ${rundir} + chmod 644 ${rundir}/kpp_standalone_interface.yml +fi + # Set permissions chmod 744 ${rundir}/cleanRunDir.sh chmod 744 ${rundir}/archiveRun.sh diff --git a/run/shared/cleanRunDir.sh b/run/shared/cleanRunDir.sh index 2c2ca77e0..e0b216848 100755 --- a/run/shared/cleanRunDir.sh +++ b/run/shared/cleanRunDir.sh @@ -1,8 +1,20 @@ #!/bin/bash -rm -fv trac_avg.* -rm -fv tracerinfo.dat -rm -fv diaginfo.dat +#============================================================================ +# cleanRunDir.sh: Removes files created by GEOS-Chem from a run directory +# +# Usage: +# ------ +# $ ./cleanRunDir.sh # Removes model output files in the run directory. +# # Also prompts the user before removing diagnostic +# # output files in OutputDir/. +# +# $ ./cleanRunDir.sh 1 # Removes model ouptut files in the run directory, +# # but will remove diagnostic output files without +# # prompting first. USE WITH CAUTION! +#============================================================================ + +# Clean model output files in the run directory rm -fv gcchem* rm -fv *.rcx rm -fv *~ @@ -22,14 +34,21 @@ rm -fv EGRESS rm -fv core.* rm -fv PET*.ESMF_LogFile rm -fv allPEs.log +rm -fv fort.* -# Clean data too. If an argument is passed, then prompt user to confirm -# perhaps asking if they want to archive before deletion. -if [[ "x${1}" == "x" ]]; then - rm -Iv ./OutputDir/*.nc* # Get confirmation from user -else - rm -fv ./OutputDir/*.nc* # Skip confirmation from user +#---------------------------------------------------------------------------- +# Clean data files in OutputDir. +# These are netCDF files (*.nc) and KPP standalone interface files (*.txt). +#---------------------------------------------------------------------------- +if [[ "x${1}" == "x" ]]; then # User confirmation required + rm -Iv ./OutputDir/*.nc* + rm -Iv ./OutputDir/*.txt +else # User Confirmation not required + rm -fv ./OutputDir/*.nc* + rm -fv ./OutputDir/*.txt* fi +#--------------------------------------------------------------------------- # Give instruction to reset start date if using GCHP +#--------------------------------------------------------------------------- echo "Reset simulation start date in cap_restart if using GCHP" diff --git a/run/shared/kpp_standalone_interface.yml b/run/shared/kpp_standalone_interface.yml new file mode 100644 index 000000000..634751772 --- /dev/null +++ b/run/shared/kpp_standalone_interface.yml @@ -0,0 +1,99 @@ +--- +# ============================================================================ +# Configuration file for KPP standalone interface +# +# This file specifies at which locations we will archive the model +# state so that we can initialize KPP standalone box model simulations. +# ============================================================================ + +# ------------------------------------ +# General settngs +# ------------------------------------ +settings: + activate: false # Main on/off switch + start_output_at: [19000101, 000000] # Save model state for KPP standalone + stop_output_at: [21000101, 000000] # ... if between these 2 datetimes + output_directory: "./OutputDir/" # This directory should already exist + levels: # Model levels to archive + - 1 + - 2 + - 10 + - 23 + - 35 + - 48 + - 56 + timestep: 15 # Timestep (mins) for KPP standalone + +# ------------------------------------ +# Where to archive model state? +# ------------------------------------ +active_cells: + - LosAngeles + - McMurdo + - Paris + - Beijing + - Kinshasa + - Kennaook + - Graciosa + - Utqiagvik + - Ozarks + - Amazon + - Congo + - Borneo + - IndianOcean + - AtlanticOcean + - PacificOcean + - ElDjouf + +# ------------------------------------ +# Active cell geographic coordinates +# ------------------------------------ +locations: + LosAngeles: + longitude: -118.243 + latitude: 34.0522 + Paris: + longitude: 2.3522 + latitude: 48.8566 + Beijing: + longitude: 116.4074 + latitude: 39.9042 + Kinshasa: + longitude: 15.3105 + latitude: -4.3033 + Kennaook: + longitude: 144.6833 + latitude: -40.6833 + Graciosa: + longitude: -28.0069 + latitude: 39.0525 + McMurdo: + longitude: 166.6698 + latitude: -77.8455 + Utqiagvik: + longitude: -156.7886 + latitude: 71.2906 + Ozarks: + longitude: -91.259 + latitude: 37.502 + Amazon: + longitude: -62.2159 + latitude: -3.4653 + Congo: + longitude: 12.5484 + latitude: -5.9175 + Borneo: + longitude: 114.0 + latitude: 0.0 + IndianOcean: + longitude: 87.2 + latitude: 23.0 + AtlanticOcean: + longitude: -41.574755 + latitude: 34.707874 + PacificOcean: + longitude: -121.964508 + latitude: 0.0 + ElDjouf: + longitude: -6.6661 + latitude: 21.5008