diff --git a/CMakeLists.txt b/CMakeLists.txt index 9bdf6c4e..557fc673 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -54,5 +54,6 @@ elseif( ) endif() +add_subdirectory(util) add_subdirectory(src) add_subdirectory(tests) diff --git a/src/api/catchem.F90 b/src/api/catchem.F90 index db953d4c..6392be00 100644 --- a/src/api/catchem.F90 +++ b/src/api/catchem.F90 @@ -82,6 +82,11 @@ module CATChem use CCPr_DryDep_mod, only: cc_drydep_init => CCPr_DryDep_Init !< DryDep Process Initialization Routine use CCPr_DryDep_mod, only: cc_drydep_run => CCPr_DryDep_Run !< DryDep Process Run Routine use CCPr_DryDep_mod, only: cc_drydep_finalize => CCPr_DryDep_Finalize !< DryDep Process Finalization Routine + ! SUVolcanicEmissions + use CCPr_SUVolcanicEmissions_mod, only: SUVolcanicEmissionsStateType !< SUVolcanicEmissions State + use CCPr_SUVolcanicEmissions_mod, only: cc_suvolcanic_init => CCPr_SUVolcanicEmissions_Init !< SUVolcanicEmissions Process Initialization Routine + use CCPr_SUVolcanicEmissions_mod, only: cc_suvolcanic_run => CCPr_SUVolcanicEmissions_Run !< SUVolcanicEmissions Process Run Routine + use CCPr_SUVolcanicEmissions_mod, only: cc_suvolcanic_finalize => CCPr_SUVolcanicEmissions_Finalize !< SUVolcanicEmissions Process Finalization Routine implicit none diff --git a/src/core/chemstate_mod.F90 b/src/core/chemstate_mod.F90 index 00ac1318..df219bd7 100644 --- a/src/core/chemstate_mod.F90 +++ b/src/core/chemstate_mod.F90 @@ -40,7 +40,7 @@ module ChemState_Mod !! \param nSpeciesGas: The number of gas species. !! \param nSpeciesAero: The number of aerosol species. !! \param nSpeciesDust: The number of dust species. - !! \param nSpeicesSeaSalt: The number of sea salt species. + !! \param nSpeciesSeaSalt: The number of sea salt species. !! \param SpeciesIndex: An array containing the total species index. !! \param AeroIndex: An array containing the aerosol species index. !! \param GasIndex: An array containing the gas species index. @@ -66,6 +66,7 @@ module ChemState_Mod INTEGER :: nSpeciesTracer !< Number of Tracer Species INTEGER :: nSpeciesDust !< Number of Dust Species INTEGER :: nSpeciesSeaSalt !< Number of SeaSalt Species + INTEGER :: nSpeciesSUVolcanic !< Number of Volcanic Sulfur Species INTEGER, ALLOCATABLE :: SpeciesIndex(:) !< Total Species Index INTEGER, ALLOCATABLE :: TracerIndex(:) !< Tracer Species Index INTEGER, ALLOCATABLE :: AeroIndex(:) !< Aerosol Species Index @@ -73,6 +74,7 @@ module ChemState_Mod INTEGER, ALLOCATABLE :: DustIndex(:) !< Dust Species Index INTEGER, ALLOCATABLE :: SeaSaltIndex(:) !< SeaSalt Species Index INTEGER, ALLOCATABLE :: DryDepIndex(:) !< SeaSalt Species Index + INTEGER, ALLOCATABLE :: SUVolcanicIndex(:) !< Volcanic sulfur Species Index CHARACTER(len=50), ALLOCATABLE :: SpeciesNames(:) !< Species Names !--------------------------------------------------------------------- @@ -183,6 +185,7 @@ subroutine Find_Number_of_Species(ChemState, RC) ChemState%nSpeciesGas = 0 ChemState%nSpeciesSeaSalt = 0 ChemState%nSpeciesTracer = 0 + ChemState%nSpeciesSUVolcanic = 0 ! Count number of species do i = 1, ChemState%nSpecies @@ -204,6 +207,9 @@ subroutine Find_Number_of_Species(ChemState, RC) if (ChemState%ChemSpecies(i)%is_drydep .eqv. .true.) then ChemState%nSpeciesAeroDryDep = ChemState%nSpeciesAeroDryDep + 1 endif + if (ChemState%ChemSpecies(i)%is_suvolcanic .eqv. .true.) then + ChemState%nSpeciesSUVolcanic = ChemState%nSpeciesSUVolcanic + 1 + endif enddo end subroutine Find_Number_of_Species @@ -238,6 +244,7 @@ subroutine Find_Index_of_Species(ChemState, RC) integer :: seasalt_index ! Current Seas Salt Index integer :: tracer_index ! Current Tracer Index integer :: drydep_index ! Current DryDep Index + integer :: suvolcanic_index ! Current SUVolcanic Index ! Initialize @@ -253,6 +260,7 @@ subroutine Find_Index_of_Species(ChemState, RC) seasalt_index = 1 tracer_index = 1 drydep_index = 1 + suvolcanic_index = 1 ! Allocate index arrays ALLOCATE(Chemstate%AeroIndex(ChemState%nSpeciesAero), STAT=RC) @@ -297,6 +305,13 @@ subroutine Find_Index_of_Species(ChemState, RC) RETURN ENDIF + ALLOCATE(Chemstate%SUVolcanicIndex(ChemState%nSpeciesSUVolcanic), STAT=RC) + IF ( RC /= CC_SUCCESS ) THEN + errMsg = 'Error allocating Chemstate%SUVolcanicIndex' + call CC_Error(errMsg, RC, thisLoc) + RETURN + ENDIF + ! Find indices for species groups do n = 1, ChemState%nSpecies if (ChemState%ChemSpecies(n)%is_aerosol .eqv. .true.) then @@ -323,6 +338,10 @@ subroutine Find_Index_of_Species(ChemState, RC) Chemstate%DryDepIndex(drydep_index) = n drydep_index = drydep_index + 1 endif + if (ChemState%ChemSpecies(n)%is_suvolcanic .eqv. .true.) then + Chemstate%SUVolcanicIndex(suvolcanic_index) = n + suvolcanic_index = suvolcanic_index + 1 + endif enddo end subroutine Find_index_of_Species diff --git a/src/core/config_mod.F90 b/src/core/config_mod.F90 index d66b1968..fd9dd718 100644 --- a/src/core/config_mod.F90 +++ b/src/core/config_mod.F90 @@ -173,6 +173,15 @@ SUBROUTINE Read_Input_File( Config , GridState, EmisState, ChemState, RC, Config RETURN ENDIF + call Config_Process_SUVolcanicEmissions(ConfigInput, Config, RC) + IF ( RC /= CC_SUCCESS ) THEN + errMsg = 'Error in "Config_Process_SUVolcanicEmissions"!' + CALL CC_Error( errMsg, RC, thisLoc ) + CALL QFYAML_CleanUp( ConfigInput ) + CALL QFYAML_CleanUp( ConfigAnchored ) + RETURN + ENDIF + !======================================================================== ! Config ChemState @@ -569,6 +578,8 @@ SUBROUTINE Config_Chem_State( filename, GridState, ChemState, RC ) write(*,*) '| number of tracers: ', ChemState%nSpeciesTracer write(*,*) '| number of dust: ', ChemState%nSpeciesDust write(*,*) '| number of seasalt: ', ChemState%nSpeciesSeaSalt + write(*,*) '| number of volcanic: ', ChemState%nSpeciesSUVolcanic + write(*,*) '=========================================================' END SUBROUTINE Config_Chem_State @@ -1373,4 +1384,78 @@ SUBROUTINE Config_Process_DryDep( ConfigInput, Config, RC ) END SUBROUTINE Config_Process_DryDep + !> \brief Process SUVolcanicEmissions configuration + !! + !! This function processes the SUVolcanic configuration and performs the necessary actions based on the configuration. + !! + !! \param[in] ConfigInput The YAML configuration object + !! \param[inout] Config The configuration object + !! \param[out] RC The return code + !! + !! \ingroup core_modules + !!!> + SUBROUTINE Config_Process_SUVolcanicEmissions( ConfigInput, Config, RC ) + USE CharPak_Mod, ONLY : StrSplit + USE Error_Mod + USE Config_Opt_Mod, ONLY : ConfigType + + TYPE(QFYAML_t), INTENT(INOUT) :: ConfigInput ! YAML Config object + TYPE(ConfigType), INTENT(INOUT) :: Config ! Input options + + ! + ! !OUTPUT PARAMETERS: + ! + INTEGER, INTENT(OUT) :: RC ! Success or failure + + ! !LOCAL VARIABLES: + ! + ! Scalars + LOGICAL :: v_bool + INTEGER :: v_int + + ! Strings + CHARACTER(LEN=255) :: thisLoc + CHARACTER(LEN=512) :: errMsg + CHARACTER(LEN=QFYAML_StrLen) :: key + + !======================================================================== + ! Config_Process_SUVolcanicEmissions begins here! + !======================================================================== + + ! Initialize + RC = CC_SUCCESS + thisLoc = ' -> at Config_Process_SUVolcanicEmissions (in CATChem/src/core/config_mod.F90)' + errMsg = '' + + ! TODO #105 Fix reading of config file + key = "process%SUVolcanicEmissions%activate" + v_bool = MISSING_BOOL + CALL QFYAML_Add_Get( ConfigInput, TRIM( key ), v_bool, "", RC ) + IF ( RC /= CC_SUCCESS ) THEN + errMsg = 'Error parsing ' // TRIM( key ) // '!' + CALL CC_Error( errMsg, RC, thisLoc ) + RETURN + ENDIF + Config%SUVolcanicEmissions_activate = v_bool + + + key = "process%SUVolcanicEmissions%scheme_opt" + v_int = MISSING_INT + CALL QFYAML_Add_Get( ConfigInput, TRIM( key ), v_int, "", RC ) + IF ( RC /= CC_SUCCESS ) THEN + errMsg = TRIM( key ) // 'Not Found, Setting Default to 1' + v_int = 1 ! default is one + RETURN + ENDIF + Config%SUVolcanicEmissions_scheme = v_int + + + write(*,*) "SUVolcanicEmissions Configuration" + write(*,*) '------------------------------------' + write(*,*) 'Config%SUVolcanicEmissions_activate = ', Config%SUVolcanicEmissions_activate + write(*,*) 'Config%SUVolcanicEmissions_scheme = ', Config%SUVolcanicEmissions_scheme + write(*,*) '------------------------------------' + + END SUBROUTINE Config_Process_SUVolcanicEmissions + END MODULE config_mod diff --git a/src/core/config_opt_mod.F90 b/src/core/config_opt_mod.F90 index cc62f873..90501a52 100644 --- a/src/core/config_opt_mod.F90 +++ b/src/core/config_opt_mod.F90 @@ -45,6 +45,8 @@ MODULE Config_Opt_Mod !! - `drydep_activate` : Activate drydep process !! - `drydep_scheme` : Scheme option for drydep process !! - `drydep_resuspension` : Activate resuspension + !! - `SUVolcanicEmissions_activate` : Activate SU Volcanic process + !! - `SUVolcanicEmissions_scheme` : Scheme option for SU Volcanic process !! !! \ingroup core_modules !!!> @@ -102,6 +104,10 @@ MODULE Config_Opt_Mod INTEGER :: drydep_scheme LOGICAL :: drydep_resuspension !< Turn on resuspension + ! SUVolcanicEmissions Process + LOGICAL :: SUVolcanicEmissions_activate + INTEGER :: SUVolcanicEmissions_scheme + END TYPE ConfigType CONTAINS @@ -176,6 +182,10 @@ SUBROUTINE Set_Config( am_I_Root, Config, RC ) Config%drydep_scheme = 1 Config%drydep_resuspension = .FALSE. + ! SU Volcanic Process + Config%SUVolcanicEmissions_activate = .FALSE. + Config%SUVolcanicEmissions_scheme = 1 + END SUBROUTINE Set_Config !> \brief Cleanup the Config options !! diff --git a/src/core/metstate_mod.F90 b/src/core/metstate_mod.F90 index 164ccdb1..e42a33f9 100644 --- a/src/core/metstate_mod.F90 +++ b/src/core/metstate_mod.F90 @@ -45,6 +45,8 @@ MODULE MetState_Mod ! TIMESTEP !--------- REAL(fp) :: TSTEP !< Time step [s] + INTEGER :: YMD ! Year, month, day + INTEGER :: HMS ! Hour, minute, second ! Logicals !--------- diff --git a/src/core/species_mod.F90 b/src/core/species_mod.F90 index 338d18bc..bf7cf24c 100644 --- a/src/core/species_mod.F90 +++ b/src/core/species_mod.F90 @@ -41,6 +41,7 @@ module species_mod logical :: is_gocart_aero !< if true, species is a GOCART aerosol species logical :: is_dust !< if true, species is a dust logical :: is_seasalt !< if true, species is a seasalt + logical :: is_suvolcanic ! Numerical properties real(kind=fp) :: mw_g !< gaseous molecular weight @@ -59,6 +60,7 @@ module species_mod integer :: drydep_index !< drydep index in drydep array integer :: photolysis_index !< photolysis index in photolysis array integer :: gocart_aero_index !< gocart_aero index in gocart_aero array + integer :: suvolcanic_index !< drydep index in drydep array ! Concentration real(kind=fp), ALLOCATABLE :: conc(:) !< species concentration [v/v] or kg/kg diff --git a/src/process/CMakeLists.txt b/src/process/CMakeLists.txt index d149b457..b7b1459a 100644 --- a/src/process/CMakeLists.txt +++ b/src/process/CMakeLists.txt @@ -4,3 +4,4 @@ add_subdirectory(dust) add_subdirectory(seasalt) add_subdirectory(plumerise) add_subdirectory(drydep) +add_subdirectory(SUVolcanicEmissions) diff --git a/src/process/SUVolcanicEmissions/CCPr_SUVolcanicEmissions_Mod.F90 b/src/process/SUVolcanicEmissions/CCPr_SUVolcanicEmissions_Mod.F90 new file mode 100644 index 00000000..8f9068f7 --- /dev/null +++ b/src/process/SUVolcanicEmissions/CCPr_SUVolcanicEmissions_Mod.F90 @@ -0,0 +1,281 @@ +!> \brief CCPR suvolcanicemissions state types +!! +!! \defgroup catchem_suvolcanicemissions_process +!! +!! \author Lacey Holland +!! \date 10/2024 +!!!> +MODULE CCPR_SUVolcanicEmissions_mod + USE Precision_mod + USE Error_Mod + USE DiagState_Mod, Only : DiagStateType + USE MetState_Mod, Only : MetStateType + USE ChemState_Mod, Only : ChemStateType + USE Config_Opt_Mod, Only : ConfigType + + IMPLICIT NONE + + PRIVATE + + PUBLIC :: CCPR_SUVolcanicEmissions_Init + PUBLIC :: CCPR_SUVolcanicEmissions_Run + PUBLIC :: CCPR_SUVolcanicEmissions_Finalize + PUBLIC :: SUVolcanicEmissionsStateType + + + !> \brief SUVolcanicEmissionsStateType + !! + !! SUVolcanicEmissionsStateType is the process-specific derived type. + !! + !! \param Activate Activate Process (True/False) + !! \param Scheme Scheme Option + !! \param SUVolcanicSpeciesIndex Effected Chemical Species from Volcanic + !! \param nSpc # of species + !! \param SpcIDs CATChem species IDs + !! \param ScaleFactor Scale Factor + !! \param Resuspension Activate resuspension (True/False) + !! + !! \ingroup core_modules + !!!> + TYPE :: SUVolcanicEmissionsStateType + + ! Process Specific Parameters + + ! Namelist parameters for specific Volcanic goes here as well + !================================================================= + ! Module specific variables/arrays/data pointers come below + !================================================================= + + LOGICAL :: Activate ! Activate Process (True/False) + INTEGER :: SchemeOpt ! Scheme Option (if there is only one SchemeOpt always = 1) + + + END TYPE SUVolcanicEmissionsStateType + + +CONTAINS + + !> + !! \brief Initialize the CATChem Volcanic module + !! + !! \param Config CATCHem configuration options + !! \param SUVolcanicEmissionsState CATCHem PROCESS state + !! \param ChemState CATCHem chemical state + !! \param RC Error return code + !! + !! \ingroup catchem_suvolcanicemissions_process + !! + !!!> + SUBROUTINE CCPR_SUVolcanicEmissions_Init( Config, SUVolcanicEmissionsState, ChemState, RC ) + ! USE + + + IMPLICIT NONE + ! INPUT PARAMETERS + !----------------- + TYPE(ConfigType) :: Config ! Module options + TYPE(ChemStateType) :: ChemState ! Chemical state + + ! INPUT/OUTPUT PARAMETERS + !------------------------ + TYPE(SUVolcanicEmissionsStateType) :: SUVolcanicEmissionsState ! Volcanic state + INTEGER, INTENT(INOUT) :: RC ! Success or failure + + ! Error handling + !--------------- + CHARACTER(LEN=255) :: ErrMsg + CHARACTER(LEN=255) :: ThisLoc + + ! LOCAL VARIABLES + !---------------- + + ! Put any local variables here + + + + !================================================================= + ! CCPR_DryDep_Init begins here! + !================================================================= + ErrMsg = '' + ThisLoc = ' -> at CCPR_SUVolcanicEmissions_INIT (in process/SUVolcanicEmissions/CCPr_SUVolcanicEmissions_mod.F90)' + + ! First check if process is activated in config | if not don't allocate arrays or pointers + if (Config%suvolcanicemissions_activate) then + + ! Activate Process + !------------------ + SUVolcanicEmissionsState%Activate = .true. + + ! Set scheme option + !------------------ + ! For now, the only option is SchemeOpt = 1 + SUVolcanicEmissionsState%SchemeOpt = 1 + + else + SUVolcanicEmissionsState%Activate = .false. + end if + + end subroutine CCPR_SUVolcanicEmissions_Init + + !> + !! \brief Run the SUVolcanicEmissions + !! + !! \param [IN] MetState - The MetState object + !! \param [INOUT] DiagState - The DiagState object + !! \param [INOUT] SUVolcanicEmissionsState - The SUVolcanicEmissionsState object + !! \param [INOUT] ChemState - The ChemState object + !! \param [OUT] RC Return code + !! + !! \ingroup catchem_suvolcanicemissions_process + !!!> + SUBROUTINE CCPr_SUVolcanicEmissions_Run( MetState, DiagState, SUVolcanicEmissionsState, ChemState, RC ) + + ! USE + USE constants, only : g0 + use CCPr_Scheme_GOCART_SUVolcanicEmissions_Mod, only : CCPr_Scheme_GOCART_SUVolcanicEmissions + + IMPLICIT NONE + ! INPUT PARAMETERS + TYPE(MetStateType), INTENT(IN) :: MetState !< MetState Instance + + ! INPUT/OUTPUT PARAMETERS + TYPE(DiagStateType), INTENT(INOUT) :: DiagState !< DiagState Instance + TYPE(SUVolcanicEmissionsStateType), INTENT(INOUT) :: SUVolcanicEmissionsState !< SUVolcanicEmissionsState Instance + TYPE(ChemStateType), INTENT(INOUT) :: ChemState !< ChemState Instance + + ! OUTPUT PARAMETERS + INTEGER, INTENT(OUT) :: RC ! Return Code + + ! LOCAL VARIABLES + CHARACTER(LEN=255) :: ErrMsg, thisLoc + INTEGER :: i ! loop index + INTEGER :: hms ! Model time [secs] + + INTEGER, dimension(:), allocatable :: vStart ! Emissions Start time [sec] + INTEGER, dimension(:), allocatable :: vEnd ! Emissions end time [sec] + INTEGER :: nVolc ! number of volcanic sources + INTEGER, dimension(:), allocatable :: iPoint, jPoint ! sub-domain - we only run this at the place/time of eruption?? + INTEGER :: nSO2 ! index of SO2 relative to other sulfate tracers + + REAL, dimension(:), allocatable :: vSO2 ! volcanic emissions [kg] + REAL, dimension(:,:,:),pointer :: SO2 ! SO2 [kg kg-1] + REAL, dimension(:,:,:),pointer :: SU_emis ! SU emissions, kg/m2/s + REAL, dimension(:), allocatable :: vCloud ! top elevation of emissions [m] + REAL, dimension(:), allocatable :: vElev ! bottom elevation of emissions [m] + REAL, dimension(:), allocatable :: vLat ! latitude specified in file [degree] + REAL, dimension(:), allocatable :: VLon ! longitude specified in file [degree] + REAL, dimension(:,:), allocatable :: area ! longitude specified in file [degree] + + + REAL(fp) :: SpecConc ! Temporary Species concentration + + ! Initialize + RC = CC_SUCCESS + errMsg = '' + thisLoc = ' -> at CCPr_SUVolcanicEmissions_Run & + & (in process/SUVolcanicEmissions/ccpr_SUVolcanicEmissions_mod.F90)' + + +! Get pointwise SO2 and altitude of volcanoes from a daily file data base +! if(index(self%volcano_srcfilen,'volcanic_') /= 0) then +! call ReadASCIIPointEmissions (YMD, fname, nVolc, vLat, vLon, & +! vElev, vCloud, vSO2, vStart, vEnd, label, RC) +! call ReadASCIIPointEmissions (fname, label, SUVolcanicEmissionsState, RC) + + nSO2 = ChemState%nSpeciesSUVolcanic !HC'd sulfate tracer number in GOCART??? + + ! Run SUVolcanic + !------------------------- + if (SUVolcanicEmissionsState%Activate) then + ! Run the GOCART SUVolcanic Scheme + !------------------------- + if (SUVolcanicEmissionsState%SchemeOpt == 1) then + ! Run the SU Volcanic Scheme - Only Applicable to AEROSOL species + !------------------------- + + if (ChemState%nSpeciesSUVolcanic > 0) then + + ! loop through aerosol species + do i = 1, ChemState%nSpeciesSUVolcanic + + ! Right now, GOCART has SO2 as the only volcanic species + ! Need to look up which is the index for SO2 concentrations + call CCPr_Scheme_GOCART_SUVolcanicEmissions( MetState%NLEVS, & + MetState%TSTEP, & + VStart, & + VEnd, & + nVolc, & + iPoint, & + jPoint, & + HMS, & + g0, & + MetState%BXHEIGHT, & + MetState%DELP, & + area, & + vSO2, & !volcanic contribution to so2 emissions + nSO2, & ! tracer number for so2 within sulfur trace + SO2, & !total so2 concentration intent(inout) + SU_emis, & !total emission rate for each sulfur species, !SU_emis(:,:,nSO2), nSO2=2, nDMS=1, nSO4=3, nMSA=4 + vCloud, & + vElev, & + vLat, & + VLon, & + RC) + ! nso2 is used to define SU_emis: SU_emis(:,:,nSO2) + + end do ! do i = 1, ChemState%nSpeciesSUVolcanic + + endif ! if (ChemState%nSpeciesSUVolcanic > 0) + + endif ! if (SUVolcanicEmissionsState%SchemeOpt == 1) + + write(*,*) 'TODO: Need to figure out how to add back to the chemical species state ' + + endif ! if (SUVolcanicEmissionsState%Activate) + + end subroutine CCPr_SUVolcanicEmissions_Run + + !> + !! \brief Finalize the DryDep + !! + !! \param [INOUT] SUVolcanicEmissionsState + !! \param [OUT] RC Return code + !!!> + SUBROUTINE CCPr_SUVolcanicEmissions_Finalize( SUVolcanicEmissionsState, RC ) + + ! USE + !---- + + IMPLICIT NONE + + ! INPUT/OUTPUT PARAMETERS + TYPE(SUVolcanicEmissionsStateType), INTENT(INOUT) :: SUVolcanicEmissionsState ! SUVolcanicEmissionsState Instance + + ! OUTPUT PARAMETERS + INTEGER, INTENT(OUT) :: RC ! Return Code + + ! LOCAL VARIABLES + CHARACTER(LEN=255) :: ErrMsg, thisLoc + + ! Initialize + RC = CC_SUCCESS + errMsg = '' + thisLoc = ' -> at CCPr_SUVolcanicEmissions_Finalize & + &(in process/SUVolcanicEmissions/ccpr_SUVolcanicEmissions_mod.F90)' + +! DEALLOCATE(SUVolcanicEmissionsStateType%SU_emiss, STAT=RC) +! IF ( RC /= CC_SUCCESS ) THEN +! ErrMsg = 'Could not Deallocate SUVolcanicEmissionsStateType%SU_emiss(:,:,nSO2)' +! CALL CC_Error( ErrMsg, RC, ThisLoc ) +! ENDIF +! DEALLOCATE(SUVolcanicEmissionsStateType%vSO2, STAT=RC) +! IF ( RC /= CC_SUCCESS ) THEN +! ErrMsg = 'Could not Deallocate SUVolcanicEmissionsStateType%vSO2' +! CALL CC_Error( ErrMsg, RC, ThisLoc ) +! ENDIF + + + end subroutine CCPr_SUVolcanicEmissions_Finalize + + +END MODULE CCPR_SUVolcanicEmissions_Mod diff --git a/src/process/SUVolcanicEmissions/CMakeLists.txt b/src/process/SUVolcanicEmissions/CMakeLists.txt new file mode 100644 index 00000000..70625e50 --- /dev/null +++ b/src/process/SUVolcanicEmissions/CMakeLists.txt @@ -0,0 +1,16 @@ +set( + _srcs + CCPr_SUVolcanicEmissions_Mod.F90 + ccpr_scheme_GOCART_SUVolcanicEmissions_mod.F90 +) + +set(_lib CATChem_process_suvolcanic) + +add_library(${_lib} ${_srcs}) +target_link_libraries(${_lib} PUBLIC CATChem_core) +target_link_libraries(${_lib} PUBLIC GOCART2G) + +set_target_properties( + ${_lib} + PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR}/include +) diff --git a/src/process/SUVolcanicEmissions/ccpr_scheme_GOCART_SUVolcanicEmissions_mod.F90 b/src/process/SUVolcanicEmissions/ccpr_scheme_GOCART_SUVolcanicEmissions_mod.F90 new file mode 100644 index 00000000..20f932ca --- /dev/null +++ b/src/process/SUVolcanicEmissions/ccpr_scheme_GOCART_SUVolcanicEmissions_mod.F90 @@ -0,0 +1,167 @@ +!> +!! \file +!! \brief CCPr Scheme for Volcanic Emissions +!! +!! +!! Reference: Benchmarking GOCART-2G in the Goddard Earth Observing System (GEOS) +!! Allison B. Collow, Peter R. Colarco, Arlindo M. da Silva, Virginie Buchard, +!! Huisheng Bian, M Chin, Sampa Das, Ravi Govindaraju, Dongchul Kim, and Valentina Aquila, +!! Geosci. Model Development, 17, 14431468, 2024 +!! https://doi.org/10.5194/gmd-17-1443-2024 +!! +!! \author Lacey Holland +!! \date 07/2024 +!!!> +module CCPr_Scheme_GOCART_SUVolcanicEmissions_Mod + + implicit none + + private + + public :: CCPr_Scheme_GOCART_SUVolcanicEmissions + +contains + + !> \brief Brief description of the subroutine + !! + !! \param MetState Meteorological Variables + !! \param DiagState Diagnostic Variables + !! \param EmisState Emission State variables + !! \param SUVolcanicEmissions SUVolcanicEmissions Variables + !! \param RC Success or Failure + !! + !! Note that other state types may be required, e.g. one specific to the process group. + !!!> + + ! Need to change this so that it includes only what is necessary + ! May need to do something to connect iPoint, jPoint to vlat, vlon + ! should only run this where we know there is a volcano???? + !subroutine CCPr_Scheme_GOCART_SUVolcanicEmissions( MetState, DiagState, EmisState, & + ! SUVolcanicEmissionsState, & + ! RC) + subroutine CCPr_Scheme_GOCART_SUVolcanicEmissions(km, & + cdt, & + VStart, & + VEnd, & + nVolc, & + iPoint, & + jPoint, & + !YMD, & + HMS, & + g0, & + zbox, & + delp, & + area, & + vSO2, & + nSO2, & + SO2, & + SU_emis, & + vCloud, & + vElev, & + vLat, & + VLon, & + RC ) + + USE GOCART2G_Process, only: SUVolcanicEmissions + USE ReadEmissions, only: ReadASCIIPointEmissions + + IMPLICIT NONE + + ! Arguments + INTEGER, intent(in) :: km ! number of vertical levels + + INTEGER, intent(inout),dimension(1) :: vStart ! Emissions Start time [sec] + INTEGER, intent(inout),dimension(1) :: vEnd ! Emissions end time [sec] + INTEGER, intent(inout) :: nVolc ! number of volcanic sources + INTEGER, intent(inout) :: rc ! error code - is this inout or out??? + INTEGER, intent(inout),dimension(1) :: iPoint, jPoint ! sub-domain - we only run this at the place/time of eruption?? + !INTEGER, intent(in) :: YMD + INTEGER, intent(in) :: HMS ! current model time [sec] + INTEGER, intent(inout) :: nSO2 ! index of SO2 relative to other sulfate tracers + + + REAL, intent(in) :: g0 + REAL, intent(in) :: cdt ! model timestep [sec] + + REAL, intent(inout),dimension(:,:) :: area ! area of grid cell [m^2] + REAL, intent(inout),dimension(:) :: vSO2 ! volcanic emissions [kg] + !!!! Below can be figured out from ChemSpeciesState%nSpeciesSUVolcanicIndex??? + REAL, intent(inout),dimension(:,:,:),pointer :: SO2 ! SO2 [kg kg-1] + REAL, intent(inout),dimension(:,:,:),pointer :: SU_emis ! SU emissions, kg/m2/s + REAL, intent(inout),dimension(:) :: vCloud ! top elevation of emissions [m] + REAL, intent(inout),dimension(:) :: vElev ! bottom elevation of emissions [m] + REAL, intent(inout),dimension(:) :: vLat ! latitude specified in file [degree] + REAL, intent(inout),dimension(:) :: VLon ! longitude specified in file [degree] + + !CHARACTER :: fname + REAL, DIMENSION(:,:),pointer :: SO2EMVN ! non-explosive volcanic emissions [kg m-2 s-1] + REAL, DIMENSION(:,:),pointer :: SO2EMVE ! explosive volcanic emissions [kg m-2 s-1] + REAL, allocatable, DIMENSION(:) :: delp ! Pressure Thickness for layer [Pa] + REAL, allocatable, DIMENSION(:) :: zbox ! geopotential Height difference [m] for layer + + !TYPE(MetStateType), INTENT(IN) :: MetState ! MetState Instance + !TYPE(DiagStateType), INTENT(IN) :: DiagState ! DiagState Instance + !TYPE(SUVolcanicEmissionsStateType), INTENT(IN) :: SUVolcanicEmissionsState ! SUVolcanicEmissionsState Instance + + ! should these be pulled from ChemSpeciesState and the species YAML? + REAL, parameter :: fMassSulfur = 32. ! gram molecular weights of species + REAL, parameter :: fMassSO2 = 64. ! gram molecular weights of species + !REAL, parameter :: fMassSO4 = 96. ! gram molecular weights of species + !CHARACTER(len=7), parameter :: label='volcano' + + + ! Local Variables + real, pointer :: GOCART_ZBOX(:,:,:) + real, pointer :: GOCART_DELP(:,:,:) + + character(len=256) :: errMsg + character(len=256) :: thisLoc + + ! Initialize + errMsg = '' + thisLoc = ' -> at CCPr_Scheme_GOCART_SUVolcanicEmissions & + & (in CCPr_Scheme_GOCART_SUVolcanicEmissions_mod.F90)' + RC = 0 + + ! transform data for GOCART SUVolcanicEmissions call + + ! Need to re-write for individual variables + ! put into the src/core directory as a module??? + call PrepMetVarsForGOCARTSUV(km, & + delp, & + zbox, & + GOCART_DELP, & + GOCART_ZBOX) + + !------------------ + ! Begin Scheme Code + !------------------ + + !!! I don't know if we need the line below, if we are just reading in?? + vSO2 = vSO2 * fMassSO2 / fMassSulfur + +!!!!!!!!!!!!!!!!!!!!! NEED TO EDIT ABOVE !!!!!!!!!!!!!!!!!!!!! +!!!! Most of the volcanic stuff, is going directly into the call below + +!! Need to replace below. Not sure if ipoint, jpoint coordinates are necessary + if (nVolc > 0) then + + iPoint(1) = 0 + jPoint(1) = 0 + + call SUvolcanicEmissions (nVolc, vStart, vEnd, vSO2, & + vElev, vCloud, & + iPoint, jPoint, & + hms, SO2EMVN, SO2EMVE, SO2, nSO2, & + SU_emis, km, cdt, g0, gocart_ZBOX, gocart_DELP, area, & + vLat, vLon, rc) + + end if + + if (associated(GOCART_DELP)) nullify(GOCART_DELP) + if (associated(GOCART_zbox)) nullify(GOCART_zbox) + + + end subroutine CCPr_Scheme_GOCART_SUVolcanicEmissions + +end module CCPr_Scheme_GOCART_SUVolcanicEmissions_Mod diff --git a/tests/CATChem_config.yml b/tests/CATChem_config.yml index a8e4c6ee..228ddb46 100644 --- a/tests/CATChem_config.yml +++ b/tests/CATChem_config.yml @@ -45,3 +45,6 @@ process: activate: true scheme_opt: 1 resuspension: false + SUVolcanicEmissions: + activate: true + scheme_opt: 1 diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 1dcb34a1..cb63b9de 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -81,6 +81,7 @@ add_test( add_executable(test_drydep test_drydep.f90) target_link_libraries(test_drydep PRIVATE CATChem_core) target_link_libraries(test_drydep PRIVATE CATChem) +target_link_libraries(test_drydep PRIVATE CATChem_process_drydep) target_link_libraries(test_drydep PRIVATE testing) set_target_properties( test_drydep @@ -92,3 +93,19 @@ add_test( COMMAND test_drydep WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} ) + +add_executable(test_suvolcanic test_suvolcanic.F90) +target_link_libraries(test_suvolcanic PRIVATE CATChem_core) +target_link_libraries(test_suvolcanic PRIVATE CATChem) +target_link_libraries(test_suvolcanic PRIVATE CATChem_process_suvolcanic) +target_link_libraries(test_suvolcanic PRIVATE testing) +set_target_properties( + test_suvolcanic + PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR}/include +) + +add_test( + NAME test_suvolcanic + COMMAND test_suvolcanic + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} +) diff --git a/tests/Emissions/so2_volcanic_emissions_Carns.20220101.rc b/tests/Emissions/so2_volcanic_emissions_Carns.20220101.rc new file mode 100644 index 00000000..0448bc8e --- /dev/null +++ b/tests/Emissions/so2_volcanic_emissions_Carns.20220101.rc @@ -0,0 +1,93 @@ +### LAT (-90,90), LON (-180,180), SULFUR [kg S/s], ELEVATION [m], CLOUD_COLUMN_HEIGHT [m] +### If elevation=cloud_column_height, emit in layer of elevation +### else, emit in top 1/3 of cloud_column_height +volcano:: +10.030 -83.770 4.347124e+00 3340 3340 +-37.860 -71.160 1.974655e+00 2800 2800 +19.420 -155.290 2.904473e+01 1222 1222 +3.170 98.390 1.889615e+00 2460 2460 +-15.800 -71.860 5.030326e-01 5967 5967 +-15.400 167.830 1.661017e+01 1395 1395 +30.790 130.310 3.388624e+00 614 614 +13.420 -88.470 5.102394e-01 1449 1449 +-19.150 -68.830 4.540266e-01 5550 5550 +31.590 130.660 6.109900e+00 639 639 +-8.820 121.180 4.929431e-01 2124 2124 +56.640 161.340 3.065760e+00 3283 3283 +1.680 127.880 9.992909e+00 1170 1170 +-16.250 168.120 4.256463e+01 1334 1334 +1.360 124.790 1.181910e+00 1580 1580 +32.880 131.110 3.633654e+00 1181 1181 +-8.320 121.710 3.459250e-01 875 875 +-37.520 177.180 1.468740e+00 321 321 +-8.210 119.070 4.107859e-01 1949 1949 +2.930 -76.030 3.629330e+00 5364 5364 +12.700 -87.000 3.591855e+00 1745 1745 +-6.110 105.420 1.754128e+00 813 813 +54.760 -163.970 2.010689e+00 2857 2857 +29.640 129.720 4.991410e+00 560 560 +47.340 152.480 8.013929e-01 870 870 +54.050 159.450 5.272474e+00 1536 1536 +-14.270 167.500 2.512280e+00 797 797 +55.120 160.360 4.116508e+00 2376 2376 +18.140 145.790 3.372769e+00 570 570 +-7.240 109.210 1.190559e+00 3329 3329 +-5.050 151.330 3.648067e+00 2300 2300 +-1.410 29.200 2.044705e+01 2950 2950 +-25.170 -68.500 1.434147e+00 5697 5697 +-8.420 116.470 4.309649e-01 3726 3726 +56.060 160.640 3.358355e+00 4835 4835 +-16.340 -70.900 1.290012e+00 5400 5400 +38.790 15.210 1.044982e+00 870 870 +53.250 158.830 4.090563e+00 2741 2741 +60.490 -152.750 2.133204e+00 3108 3108 +-19.750 -175.070 1.646026e+00 515 515 +-5.530 148.420 3.637978e+00 991 991 +50.330 155.460 2.868295e+00 1816 1816 +-0.080 -77.660 1.192000e+00 3562 3562 +-8.510 124.130 2.156266e+00 862 862 +36.400 138.530 2.600203e+00 2126 2126 +15.550 41.830 5.981620e-01 244 244 +4.900 -75.320 6.213678e+00 5321 5321 +13.600 40.670 3.675453e-01 613 613 +51.790 -178.790 2.998017e-01 1573 1573 +48.080 153.210 1.501891e+00 1200 1200 +52.380 -174.150 1.145877e+00 1533 1533 +-7.940 112.950 4.481170e+00 2258 2258 +13.260 123.690 2.620382e+00 2462 2462 +-1.700 101.260 1.697915e+00 3760 3760 +56.170 -159.380 1.480271e+00 2507 2507 +52.450 158.200 4.361538e+00 2322 2322 +14.470 -90.880 1.461533e+00 3763 3763 +45.390 148.840 1.081016e+00 1125 1125 +10.410 123.130 4.021378e-01 2435 2435 +2.780 125.400 1.808899e+00 1780 1780 +-7.560 110.440 1.859347e-01 2968 2968 +37.730 15.000 1.176001e+01 2711 2711 +-10.380 165.800 1.484595e+00 851 851 +52.830 -169.770 8.821088e-01 1170 1170 +16.350 145.670 7.727100e+00 320 320 +-57.800 -26.490 1.523511e+00 990 990 +-8.270 123.510 3.659598e+00 1339 1339 +12.770 124.050 1.193441e+00 1500 1500 +11.980 -86.160 5.018795e+00 635 635 +59.350 -153.450 4.208754e-01 1252 1252 +14.760 -91.550 1.424058e+00 3772 3772 +-39.420 -71.930 1.625848e+00 2847 2847 +-21.230 55.710 7.783313e-01 2460 2460 +43.420 142.690 7.826553e-01 2077 2077 +-6.090 155.230 2.186967e+01 1106 1106 +61.300 -152.250 6.140169e-01 3374 3374 +-4.240 152.210 1.000876e+01 200 200 +12.280 93.860 1.405320e+00 230 230 +13.850 -89.630 5.578041e-01 2381 2381 +-0.390 100.460 1.960242e-01 2686 2686 +-8.060 114.240 3.650950e+00 2799 2799 +1.200 -77.390 1.261185e+00 4276 4276 +-58.420 -26.330 8.215719e-01 1370 1370 +-4.080 145.040 8.590471e+00 1730 1730 +57.140 -156.990 7.985102e-01 2000 2000 +-77.530 167.170 2.998017e-01 3794 3794 +-1.470 -78.440 1.978979e+00 5023 5023 +16.720 -62.180 7.499366e+00 870 870 +:: diff --git a/tests/test_suvolcanic.F90 b/tests/test_suvolcanic.F90 new file mode 100644 index 00000000..3f2e1683 --- /dev/null +++ b/tests/test_suvolcanic.F90 @@ -0,0 +1,142 @@ +program test_suvolcanic + use CATChem, fp => cc_rk + use testing_mod, only: assert + use precision_mod, only: rae + implicit none + + type(ConfigType) :: Config + type(ChemStateType) :: ChemState + type(MetStateType) :: MetState + type(DiagStateType) :: DiagState + type(SUVolcanicEmissionsStateType) :: SUVolcanicEmissionsState + type(GridStateType) :: GridState + type(EmisStateType) :: EmisState + + ! Integers + INTEGER:: rc ! Success or failure + + character(len=:), allocatable :: title + + ! Error handling + CHARACTER(LEN=512) :: errMsg + CHARACTER(LEN=255) :: thisLoc + CHARACTER(LEN=255), PARAMETER :: configFile ='Configs/Default/CATChem_config.yml' + + + thisLoc = 'test_suvolcanic -> at read CATChem_Config.yml' + errMsg = '' + rc = CC_SUCCESS + + write(*,*) ' CCCCC A TTTTTTT CCCCC H' + write(*,*) ' C A A T C H EEEE M M' + write(*,*) ' C AAAAA T C HHHHH E E M M M M' + write(*,*) ' C A A T C H H E EE M M M' + write(*,*) ' CCCCC A A T CCCCC H H EEEEE M M' + write(*,*) '' + write(*,*) '' + + !---------------------------- + ! Test 1 + !---------------------------- + + ! Read input file and initialize grid + call cc_read_config(Config, GridState, EmisState, ChemState, rc, configFile) + if (rc /= CC_success) then + errMsg = 'Error reading configuration file: ' // TRIM( configFile ) + call cc_emit_error(errMsg, rc, thisLoc) + stop 1 + endif + + + title = 'Volcanic Test 1 | Read Config' + SUVolcanicEmissionsState%Activate = .false. + call print_info(Config, SUVolcanicEmissionsState, MetState, title) + write (*,*) '-- ' + write (*,*) 'Completed ', title + write (*,*) '--' + + !---------------------------- + ! Test 2 + !---------------------------- + ! Set number of Volcanic species + + ChemState%nSpeciesSUVolcanic = 2 + SUVolcanicEmissionsState%Activate = .true. + + ! Meteorological State + MetState%NLEVS = 2 + allocate(MetState%BXHEIGHT(MetState%NLEVS)) + allocate(MetState%DELP(MetState%NLEVS)) + + MetState%DELP(1:MetState%NLEVS)= 10000 ! Need to change to something more reasonable and check units. + MetState%BXHEIGHT(1:MetState%NLEVS) = 100 ! temporary, change to something more reasonable and check units + + SUVolcanicEmissionsState%SchemeOpt = 1 + + ! Allocate DiagState + call cc_allocate_diagstate(Config, DiagState, ChemState, RC) + if (rc /= CC_SUCCESS) then + errMsg = 'Error in cc_allocate_diagstate' + stop 1 + endif + + title = "SUVolcanic Test 2 | Test GOCART SUVolcanic defaults" + + call cc_suvolcanic_init(Config, SUVolcanicEmissionsState, ChemState, rc) + if (rc /= CC_SUCCESS) then + errMsg = 'Error in cc_suvolcanic_init' + call cc_emit_error(errMsg, rc, thisLoc) + stop 1 + end if + + call cc_suvolcanic_run(MetState, DiagState, & + SUVolcanicEmissionsState, ChemState, rc) + if (rc /= CC_SUCCESS) then + errMsg = 'Error in _suvolcanicemissions_run' + call cc_emit_error(errMsg, rc, thisLoc) + stop 1 + end if + +! TODO: Pull from Volcanic Emissions sample ASCII file +! TODO: Put SO2 into SUVolcanicEmissionsState (ccpr...) + call assert(SUVolcanicEmissionsState%SO2 > 0.0_fp, "Test Sulfur Volcanic Emissions") + call print_info(Config, SUVolcanicEmissionsState, MetState, title) + call cc_suvolcanic_finalize( SUVolcanicEmissionsState, rc) + if (rc /= CC_SUCCESS) then + errMsg = 'Error in _suvolcanic_finalize' + call cc_emit_error(errMsg, rc, thisLoc) + stop 1 + end if + +contains + + subroutine print_info(Config_, SUVolcanicEmissionsState_, MetState_, title_) + + type(ConfigType), intent(in) :: Config_ + type(MetStateType), intent(in) :: MetState_ + type(SUVolcanicEmissionsStateType), intent(in) :: SUVolcanicEmissionsState_ + character(len=*), intent(in) :: title_ + + write(*,*) '=======================================' + write(*,*) title_ + write(*,*) '=======================================' + write(*,*) '*************' + write(*,*) 'Configuration ' + write(*,*) '*************' + write(*,*) 'Config%suvolcanicemissions_activate = ', Config_%suvolcanicemissions_activate + write(*,*) 'Config%suvolcanicemissions_scheme = ', Config_%suvolcanicemissions_scheme + + + if (SUVolcanicEmissionsState_%Activate) then + + write(*,*) 'SUVolcanicEmissionsState%Activate = ', SUVolcanicEmissionsState_%Activate + write(*,*) 'SUVolcanicEmissionsState%SchemeOpt = ', SUVolcanicEmissionsState_%SchemeOpt + write(*,*) 'MetState%DELP =', MetState_%DELP + write(*,*) 'MetState%BXHEIGHT = ', MetState_%BXHEIGHT + + end if + + end subroutine print_info + + +end program test_suvolcanic diff --git a/util/CMakeLists.txt b/util/CMakeLists.txt new file mode 100644 index 00000000..2e904543 --- /dev/null +++ b/util/CMakeLists.txt @@ -0,0 +1,12 @@ +set(_srcs reademissions.F90 prepmetvars.F90) + +set(_lib reademissions) + +add_library(${_lib} ${_srcs}) +target_link_libraries(${_lib} PUBLIC CATChem_core) +target_link_libraries(${_lib} PUBLIC GOCART2G) + +set_target_properties( + ${_lib} + PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR}/include +) diff --git a/util/prepmetvars.F90 b/util/prepmetvars.F90 new file mode 100644 index 00000000..5d49b1c2 --- /dev/null +++ b/util/prepmetvars.F90 @@ -0,0 +1,90 @@ +module PrepMetVars + + implicit none + + private + + public :: PrepMetVarsForGOCARTSUV + public :: PrepAnyMetVarForGOCART + +CONTAINS + + + !> + !! \brief PrepMetVarsForGOCARTSUV - Prep the meteorological variables for GOCART SUVolcanicEmissions scheme + !! + !! \param [INOUT] zmid + !! \param [INOUT] pmid + !! + !! \ingroup core_modules + !!!> + + ! Need to fix below subroutine to convert one variable at a time. + + subroutine PrepMetVarsForGOCARTSUV(km, & + delp, & + zbox, & + GOCART_DELP, & + GOCART_ZBOX) + + + IMPLICIT NONE + + ! INPUTS + INTEGER, intent(in) :: km ! number of vertical levels + REAL, intent(in), DIMENSION(:), target :: delp ! Temperature [K] + REAL, intent(in), DIMENSION(:), target :: zbox ! Height [m] + + ! INPUT/OUTPUTS + REAL, intent(inout), pointer :: GOCART_DELP(:,:,:) !< temperature [K] + REAL, intent(inout), pointer, DIMENSION(:,:,:) :: GOCART_ZBOX !< geometric height [m] + + ! OUTPUTS - Add error handling back in late + !INTEGER :: rc !< Return code + + ! Error handling + !character(len=255) :: thisloc + + allocate(GOCART_DELP(1, 1, km)) + allocate(GOCART_ZBOX(1, 1, km)) + + GOCART_DELP(1,1,:) = delp ! pressure in middle of layer + GOCART_ZBOX(1,1,:) = zbox ! mid layer geopotential height [m] + + + end subroutine PrepMetVarsForGOCARTSUV + + + subroutine PrepAnyMetVarForGOCART(km, & + var, & + GOCART_VAR) + + IMPLICIT NONE + + ! INPUTS + INTEGER, intent(in) :: km ! number of vertical levels + REAL, intent(in), DIMENSION(:), target :: var ! any variable + REAL, pointer, DIMENSION(:,:,:) :: GOCART_3D_VAR + REAL, pointer, DIMENSION(:,:) :: GOCART_2D_VAR + REAL, pointer, DIMENSION(:,:,:) :: GOCART_VAR + + + if (km == 0 ) then + ! point + allocate(GOCART_2D_VAR(1, 1)) + Allocate(GOCART_VAR(1,1,1)) + GOCART_2D_VAR(1,1) = var(1) + GOCART_2D_VAR => GOCART_VAR(:,:,1) + else + ! column + allocate(GOCART_3D_VAR(1, 1, km)) + Allocate(GOCART_VAR(1,1,km)) + GOCART_3D_VAR(1,1,:) = var + GOCART_3D_VAR => GOCART_VAR + end if + + + end subroutine PrepAnyMetVarForGOCART + + +end module PrepMetVars diff --git a/util/reademissions.F90 b/util/reademissions.F90 new file mode 100644 index 00000000..e1f5d50e --- /dev/null +++ b/util/reademissions.F90 @@ -0,0 +1,131 @@ +MODULE ReadEmissions + +! +! Author: LDH (13 Jan 2025) +! + + + implicit none + + private + + public :: ReadASCIIPointEmissions + + + ! Define a derived type to hold the data for each emission entry + ! This may need to go into the Emissions State type??? + type :: VolcanicEmissionData + real :: vlat + real :: vlon + real :: VEmis + integer :: vbase + integer :: vtop + integer :: nPts + character(len=255) :: emissfile + character(len=255) :: label + end type VolcanicEmissionData + + integer :: rc + character(len=1055) :: filename + character(len=1055) :: label + + type(VolcanicEmissionData), allocatable :: VolcanicEmiss(:) + + !filename="./so2_volcanic_emissions_Carns.20220101.rc" + !label="volcano" + !call ReadASCIIPointEmissions(filename, label, VolcanicEmiss, rc ) + +contains + + + + subroutine ReadASCIIPointEmissions (filename, label, VolcanicEmissions, rc ) + + implicit none + + integer, intent(out) :: rc + integer :: num_emiss_sources=0 + integer :: num_lines=0 + integer :: num_skip=0 + integer :: i + character(1056) :: line + character(len=1055), intent(in) :: filename + character(len=1055), intent(in) :: label + character(len=255) :: errmsg + + type(VolcanicEmissionData), allocatable :: VolcanicEmissions(:) + type(VolcanicEmissionData), allocatable :: temp_emissions(:) + allocate(VolcanicEmiss(9)) + + + ! Open the file + open(unit=10, file=filename, status='old', action='read', iostat=rc) + + if (rc /= 0) then + print *, "Error opening file:", filename, " RC=", rc + return + end if + + ! Count the number of lines in the file + readloop: do while (rc >= 0) + + read(10, '(A)', iostat=rc) line + num_lines = num_lines+1 + line = trim(line) + + if (rc /= 0) then + print *, "Error reading file:", filename, " RC=", rc + return + end if + + if (line(1:1)=="#") then + num_skip = num_skip + 1 + continue + else if (line==trim(label)//"::") then + num_skip = num_skip + 1 + continue + else if (line(1:2)=="::") then + exit + else + num_emiss_sources = num_emiss_sources + 1 + end if + + end do readloop + + rewind(10) + + ! Allocate the array to hold all entries + allocate(temp_emissions(num_emiss_sources)) + + do i = 1, num_skip + read(10, '(A)', iostat=rc) line + if (rc /= 0) return + end do + + do i = 1, num_emiss_sources + read(10, *, iostat=rc, iomsg=errmsg) temp_emissions(i)%vlat, & + temp_emissions(i)%vlon, & + temp_emissions(i)%vemis, & + temp_emissions(i)%vbase, & + temp_emissions(i)%vtop + if (rc /= 0) then + print *, "Error reading file:", trim(filename), " RC=", rc + print *, "Error message:", trim(errmsg) + return + end if + end do + + ! Close the file and transfer data to output array + close(10) + + temp_emissions%nPts = num_emiss_sources + temp_emissions%emissfile = trim(filename) + temp_emissions%label = trim(label) + + VolcanicEmissions = temp_emissions + + + end subroutine ReadASCIIPointEmissions + + +end MODULE ReadEmissions