diff --git a/src/api/CMakeLists.txt b/src/api/CMakeLists.txt index 2a338f40..853b5664 100644 --- a/src/api/CMakeLists.txt +++ b/src/api/CMakeLists.txt @@ -8,6 +8,7 @@ target_link_libraries(${_lib} PUBLIC CATChem_process_dust) target_link_libraries(${_lib} PUBLIC CATChem_process_seasalt) target_link_libraries(${_lib} PUBLIC CATChem_process_plumerise) target_link_libraries(${_lib} PUBLIC CATChem_process_drydep) +target_link_libraries(${_lib} PUBLIC CATChem_process_volcanic) target_link_libraries(${_lib} PUBLIC CATChem_process_chem) set_target_properties( ${_lib} diff --git a/src/api/catchem.F90 b/src/api/catchem.F90 index 47324865..3673f824 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_Volcanic_mod, only: VolcanicStateType !< SUVolcanic State + use CCPr_Volcanic_mod, only: cc_volcanic_init => CCPr_Volcanic_Init !< SUVolcanicEmissions Process Initialization Routine + use CCPr_Volcanic_mod, only: cc_volcanic_run => CCPr_Volcanic_Run !< SUVolcanicEmissions Process Run Routine + use CCPr_Volcanic_mod, only: cc_volcanic_finalize => CCPr_Volcanic_Finalize !< SUVolcanicEmissions Process Finalization Routine ! Chemical mechanism solver use CCPr_Chem_mod, only: cc_get_micm_version => get_micm_version diff --git a/src/core/chemstate_mod.F90 b/src/core/chemstate_mod.F90 index 00ac1318..dbc81d7f 100644 --- a/src/core/chemstate_mod.F90 +++ b/src/core/chemstate_mod.F90 @@ -40,13 +40,14 @@ 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. !! \param DustIndex: An array containing the dust species index. !! \param SeaSaltIndex: An array containing the sea salt species index. - !! \param chemSpecies: A 2-D array containing the concentration of each species. + !! \param DryDepIndex: An array containing the dry deposition species index. + !! \param SpeciesNames: A character array containing the names of the species. !! !! \ingroup core_modules !!!> diff --git a/src/core/config_mod.F90 b/src/core/config_mod.F90 index d66b1968..df920969 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_Volcanic(ConfigInput, Config, RC) + IF ( RC /= CC_SUCCESS ) THEN + errMsg = 'Error in "Config_Process_Volcanic"!' + CALL CC_Error( errMsg, RC, thisLoc ) + CALL QFYAML_CleanUp( ConfigInput ) + CALL QFYAML_CleanUp( ConfigAnchored ) + RETURN + ENDIF + !======================================================================== ! Config ChemState @@ -1373,4 +1382,89 @@ SUBROUTINE Config_Process_DryDep( ConfigInput, Config, RC ) END SUBROUTINE Config_Process_DryDep + !> \brief Process Volcanic configuration + !! + !! This function processes the Volcanic 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_Volcanic( 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 + CHARACTER(LEN=1055) :: v_str + + ! Strings + CHARACTER(LEN=255) :: thisLoc + CHARACTER(LEN=512) :: errMsg + CHARACTER(LEN=QFYAML_StrLen) :: key + + !======================================================================== + ! Config_Process_Volcanic begins here! + !======================================================================== + + ! Initialize + RC = CC_SUCCESS + thisLoc = ' -> at Config_Process_Volcanic (in CATChem/src/core/config_mod.F90)' + errMsg = '' + + ! TODO #105 Fix reading of config file + key = "process%volcanic%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%volcanic_activate = v_bool + + + key = "process%volcanic%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%volcanic_scheme = v_int + + key = "process%volcanic%filedir" + v_str = MISSING_STR + CALL QFYAML_Add_Get( ConfigInput, TRIM( key ), v_str, "", RC ) + IF ( RC /= CC_SUCCESS ) THEN + errMsg = 'Error parsing ' // TRIM( key ) // '!' + CALL CC_Warning( errMsg, RC, thisLoc ) + ENDIF + Config%volcanic_filedir = TRIM( v_str ) + + + write(*,*) "Volcanic Configuration" + write(*,*) '------------------------------------' + write(*,*) 'Config%volcanic_activate = ', Config%volcanic_activate + write(*,*) 'Config%volcanic_scheme = ', Config%volcanic_scheme + write(*,*) 'Config%volcanic_filedir = ', Config%volcanic_filedir + write(*,*) '------------------------------------' + + END SUBROUTINE Config_Process_Volcanic + END MODULE config_mod diff --git a/src/core/config_opt_mod.F90 b/src/core/config_opt_mod.F90 index cc62f873..dda9bb3a 100644 --- a/src/core/config_opt_mod.F90 +++ b/src/core/config_opt_mod.F90 @@ -45,6 +45,9 @@ MODULE Config_Opt_Mod !! - `drydep_activate` : Activate drydep process !! - `drydep_scheme` : Scheme option for drydep process !! - `drydep_resuspension` : Activate resuspension + !! - `volcanic_activate` : Activate Volcanic process + !! - `volcanic_scheme` : Scheme option for Volcanic process + !! - `volcanic_filedir` : File directory for Volcanic process !! !! \ingroup core_modules !!!> @@ -102,6 +105,11 @@ MODULE Config_Opt_Mod INTEGER :: drydep_scheme LOGICAL :: drydep_resuspension !< Turn on resuspension + ! VolcanicEmissions Process + LOGICAL :: volcanic_activate + INTEGER :: volcanic_scheme + character(len=1055) :: volcanic_filedir + END TYPE ConfigType CONTAINS @@ -176,6 +184,11 @@ SUBROUTINE Set_Config( am_I_Root, Config, RC ) Config%drydep_scheme = 1 Config%drydep_resuspension = .FALSE. + ! SU Volcanic Process + Config%volcanic_activate = .FALSE. + Config%volcanic_scheme = 1 + Config%volcanic_filedir = ' ' + 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 a563f0b7..ff09ffb8 100644 --- a/src/process/CMakeLists.txt +++ b/src/process/CMakeLists.txt @@ -4,4 +4,5 @@ add_subdirectory(dust) add_subdirectory(seasalt) add_subdirectory(plumerise) add_subdirectory(drydep) +add_subdirectory(volcanic) add_subdirectory(chem) diff --git a/src/process/volcanic/CCPr_Volcanic_Mod.F90 b/src/process/volcanic/CCPr_Volcanic_Mod.F90 new file mode 100644 index 00000000..68b51511 --- /dev/null +++ b/src/process/volcanic/CCPr_Volcanic_Mod.F90 @@ -0,0 +1,379 @@ +!> \brief CCPR Volcanice state types +!! +!! \defgroup catchem_Volcanic_process +!! +!! \author Lacey Holland and Wei Li +!! \date 10/2024 +!!!> +MODULE CCPR_Volcanic_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 + USE EmisState_Mod, Only: EmisStateType + + IMPLICIT NONE + + PRIVATE + + PUBLIC :: CCPR_Volcanic_Init + PUBLIC :: CCPR_Volcanic_Run + PUBLIC :: CCPR_Volcanic_Finalize + PUBLIC :: VolcanicStateType + + + !> \brief VolcanicStateType + !! + !! VolcanicStateType is the process-specific derived type. + !! + !! \param Activate Activate Process (True/False) + !! \param SchemeOpt Scheme Option + !! \param nVolcanicSpecies # of Volcanic species + !! \param VolcanicSpeciesIndex Index of Volcanic species + !! \param VolcanicSpeciesName Name of Volcanic species + !! \param SpcIDs CATChem species IDs + !! \param CatIndex Index of emission category in EmisState + !! \param TotalEmission Total emission of all species at each level [kg/m^2/s] + !! \param EmissionPerSpecies Emission per species at each level [kg/m^2/s] + !! \param FileDir Input file directory for reading in emissions + !! + !! \ingroup catchem_Volcanic_process + !!!> + TYPE :: VolcanicStateType + + ! Generic Variables for Every Process + LOGICAL :: Activate ! Activate Process (True/False) + INTEGER :: SchemeOpt ! Scheme Option (if there is only one SchemeOpt always = 1) + integer :: nVolcanicSpecies !< Number of Volcanic species + integer, allocatable :: VolcanicSpeciesIndex(:) !< Index of Volcanic species + character(len=31), allocatable :: VolcanicSpeciesName(:) !< name of Volcanic species + integer, allocatable :: SpcIDs(:) !< CATChem species IDs + integer :: CatIndex !< Index of emission category in EmisState + + ! Process Specific Parameters + real(fp), allocatable :: TotalEmission(:) !< Total emission of all species at each level [kg/m^2/s] + real(fp), allocatable :: EmissionPerSpecies(:,:) !< Emission per species at each level [kg/m^2/s] + character(len=1055) :: FileDir !< Input file directory for reading in emissions + + + END TYPE VolcanicStateType + + +CONTAINS + + !> + !! \brief Initialize the CATChem Volcanic module + !! + !! \param Config CATCHem configuration options + !! \param VolcanicState CATCHem PROCESS state + !! \param EmisState CATCHem emission state + !! \param RC Error return code + !! + !! \ingroup catchem_Volcanicemissions_process + !! + !!!> + SUBROUTINE CCPR_Volcanic_Init( Config, VolcanicState, EmisState, RC ) + ! USE + + + IMPLICIT NONE + ! INPUT PARAMETERS + !----------------- + TYPE(ConfigType) :: Config ! Module options + TYPE(EmisStateType) :: EmisState ! Emission state + + ! INPUT/OUTPUT PARAMETERS + !------------------------ + TYPE(VolcanicStateType) :: VolcanicState ! Volcanic state + INTEGER, INTENT(INOUT) :: RC ! Success or failure + + ! Error handling + !--------------- + CHARACTER(LEN=255) :: ErrMsg + CHARACTER(LEN=255) :: ThisLoc + + ! LOCAL VARIABLES + !---------------- + INTEGER :: c + + !================================================================= + ! CCPR_DryDep_Init begins here! + !================================================================= + ErrMsg = '' + ThisLoc = ' -> at CCPR_Volcanic_INIT (in process/Volcanic/CCPr_Volcanic_mod.F90)' + + ! First check if process is activated in config | if not don't allocate arrays or pointers + if (Config%volcanic_activate) then + + ! Activate Process + !------------------ + VolcanicState%Activate = .true. + + ! Set scheme option + !------------------ + ! For now, the only option is SchemeOpt = 1 + VolcanicState%SchemeOpt = Config%volcanic_scheme + + !Find VOLCANIC caterory index in EmisState for future use + !-------------------------------------------- + do c = 1, EmisState%nCats + if (EmisState%Cats(c)%name == 'VOLCANIC') then + VolcanicState%CatIndex = c + exit + endif + end do + + ! Set number of species from EmisState + !---------------------- + VolcanicState%nVolcanicSpecies = EmisState%Cats(VolcanicState%CatIndex)%nSpecies + + !------------------------------------ + ! Allocate emission species index + ALLOCATE( VolcanicState%VolcanicSpeciesIndex(VolcanicState%nVolcanicSpecies), STAT=RC ) + CALL CC_CheckVar('VolcanicState%VolcanicSpeciesIndex', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + VolcanicState%VolcanicSpeciesIndex = -1 + + ! Allocate emission speceis names + ALLOCATE( VolcanicState%VolcanicSpeciesName(VolcanicState%nVolcanicSpecies), STAT=RC ) + CALL CC_CheckVar('VolcanicState%VolcanicSpeciesName', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + VolcanicState%VolcanicSpeciesName = '' + + ! Allocate CatChem species index + ALLOCATE( VolcanicState%SpcIDs(VolcanicState%nVolcanicSpecies), STAT=RC ) + CALL CC_CheckVar('VolcanicState%SpcIDs', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + VolcanicState%SpcIDs = -1 + + ! Allocate emission flux + ALLOCATE( VolcanicState%EmissionPerSpecies(VolcanicState%nVolcanicSpecies, & + SIZE(EmisState%Cats(VolcanicState%CatIndex)%Species(1)%Flux)), STAT=RC ) + CALL CC_CheckVar('VolcanicState%EmissionPerSpecies', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + VolcanicState%EmissionPerSpecies = ZERO + + ! Allocate total emissions + ALLOCATE( VolcanicState%TotalEmission(SIZE(EmisState%Cats(VolcanicState%CatIndex)%Species(1)%Flux)), STAT=RC ) + CALL CC_CheckVar('VolcanicState%TotalEmission', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + VolcanicState%TotalEmission = ZERO + + ! Set the file directory + VolcanicState%FileDir = TRIM(Config%Volcanic_filedir) + + else + VolcanicState%Activate = .false. + end if + + end subroutine CCPR_Volcanic_Init + + !> + !! \brief Run the VolcanicEmissions + !! + !! \param [IN] MetState - The MetState object + !! \param [INOUT] VolcanicState - The VolcanicState object + !! \param [INOUT] EmisState - The EmisState object + !! \param [INOUT] RC Return code + !! + !! \ingroup catchem_Volcanicemissions_process + !!!> + SUBROUTINE CCPr_Volcanic_Run( MetState, VolcanicState, EmisState, RC ) + + ! USE + USE constants, only : g0 + use CCPr_Scheme_Volcanic_GOCART_Mod, only : CCPr_Scheme_Volcanic_GOCART, VolcanicEmisData, ReadASCIIPointEmis + + IMPLICIT NONE + ! INPUT PARAMETERS + TYPE(MetStateType), INTENT(IN) :: MetState !< MetState Instance + + ! INPUT/OUTPUT PARAMETERS + TYPE(VolcanicStateType), INTENT(INOUT) :: VolcanicState !< VolcanicState Instance + TYPE(EmisStateType), INTENT(INOUT) :: EmisState !< ChemState Instance + + ! OUTPUT PARAMETERS + INTEGER, INTENT(INOUT) :: RC ! Return Code + + ! LOCAL VARIABLES + CHARACTER(LEN=255) :: ErrMsg, thisLoc + CHARACTER(LEN=1055) :: fname + type(VolcanicEmisData), allocatable :: VolcanicEmis(:) + CHARACTER(len=7), parameter :: label='volcano' + INTEGER :: i ! loop index + INTEGER :: hms ! Model time [secs] TODO: format is right? + INTEGER :: ymd ! Model date [YYYYMMDD] TODO: format is right? + CHARACTER(LEN=18) :: ymd_str ! Model date in string + + 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, parameter :: nSO2 =1 ! 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(:), 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 ! area of current grid cell [m^2] + + ! Initialize + RC = CC_SUCCESS + errMsg = '' + thisLoc = ' -> at CCPr_VolcanicEmissions_Run & + & (in process/VolcanicEmissions/ccpr_VolcanicEmissions_mod.F90)' + + ! Run Volcanic + !------------------------- + if (VolcanicState%Activate) then + ! Run the GOCART Volcanic Scheme + !------------------------- + if (VolcanicState%SchemeOpt == 1) then + ! Run the SU Volcanic GOCART Scheme + !------------------------- + if (VolcanicState%nVolcanicSpecies > 0) then + + !TODO: read file name from config file based on ymd? + ymd = MetState%YMD; hms = MetState%HMS + write(ymd_str, "(i0)") ymd ! converting integer to string + fname = TRIM(VolcanicState%FileDir) // '.' // TRIM(ymd_str) // '.rc' + call ReadASCIIPointEmis (fname, label, VolcanicEmis, RC) + nVolc = VolcanicEmis(1)%nPts + allocate(vSO2(nVolc), vCloud(nVolc), vElev(nVolc), vLat(nVolc), VLon(nVolc)) + vSO2 = VolcanicEmis(:)%VEmis + vCloud = VolcanicEmis(:)%Vtop + vElev = VolcanicEmis(:)%Vbase + vLon = VolcanicEmis(:)%Vlon + vLat = VolcanicEmis(:)%Vlat + if (allocated(VolcanicEmis)) deallocate(VolcanicEmis) + + !TODO: iPoint and jPoint needs to be determined; currently gives all one. In real run, + !we could pre-select sources for the current grid cell so giving all ones can work fine. + allocate(iPoint(nVolc), jPoint(nVolc)) + iPoint = 1; jPoint = 1 + !TODO: Not sure how to get Vstart and VEnd format and values. + !Here I just start from 00:00:00 and end at 23:59:59 + allocate(vStart(nVolc), vEnd(nVolc)) + vStart = ymd + 000000 + vEnd = ymd + 235959 + hms = ymd + hms + !set area and SO2 + allocate(area(1,1), SO2(1,1, MetState%NLEVS)) + area(1,1) = MetState%AREA_M2 + SO2 = ZERO + + ! loop through all species. Right now, GOCART only has SO2 + do i = 1, VolcanicState%nVolcanicSpecies + + !Need to look up which is the index for SO2 concentrations + !TODO: is level index reversed in the GOCART??? + call CCPr_Scheme_Volcanic_GOCART( 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) + vCloud, & + vElev, & + vLat, & + VLon, & + RC) + ! nso2 is used to define SU_emis: SU_emis(:,:,nSO2). We only use SO2 for now and assign nSO2=1 + + !put it back to VolcanicState + VolcanicState%VolcanicSpeciesIndex(i) = i + VolcanicState%VolcanicSpeciesName(i) = EmisState%Cats(VolcanicState%CatIndex)%Species(i)%name + !TODO: convert unit from kg kg-1 to kg m-2 s-1; + !TODO: The test only has 8 levels and we give EmissionPerSpecies 28 levels from GridState. In real run, + ! "1:8" should be changed to ":" for all the vertical levels. + VolcanicState%EmissionPerSpecies(i,1:8) = SO2(1, 1, :) * MetState%DELP / g0 / MetState%TSTEP + VolcanicState%TotalEmission(:) = VolcanicState%TotalEmission(:) + VolcanicState%EmissionPerSpecies(i,:) + + end do ! do i = 1, VolcanicState%nVolcanicSpecies + + endif ! if (VolcanicState%nVolcanicSpecies > 0) + + endif ! if (VolcanicState%SchemeOpt == 1) + + write(*,*) 'TODO: Need to figure out how to add back to the chemical species state ' + + endif ! if (VolcanicState%Activate) + + end subroutine CCPr_Volcanic_Run + + !> + !! \brief Finalize the DryDep + !! + !! \param [INOUT] VolcanicState + !! \param [INOUT] RC Return code + !!!> + SUBROUTINE CCPr_Volcanic_Finalize( VolcanicState, RC ) + + ! USE + !---- + + IMPLICIT NONE + + ! INPUT/OUTPUT PARAMETERS + TYPE(VolcanicStateType), INTENT(INOUT) :: VolcanicState ! VolcanicState Instance + + ! OUTPUT PARAMETERS + INTEGER, INTENT(INOUT) :: RC ! Return Code + + ! LOCAL VARIABLES + CHARACTER(LEN=255) :: ErrMsg, thisLoc + + ! Initialize + RC = CC_SUCCESS + errMsg = '' + thisLoc = ' -> at CCPr_VolcanicEmissions_Finalize & + &(in process/VolcanicEmissions/ccpr_VolcanicEmissions_mod.F90)' + + ! Deallocate any arrays here + IF (ALLOCATED(VolcanicState%SpcIDs)) THEN + DEALLOCATE( VolcanicState%SpcIDs, STAT=RC ) + CALL CC_CheckVar('VolcanicState%SpcIDs', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + END IF + + IF (ALLOCATED(VolcanicState%VolcanicSpeciesIndex)) THEN + DEALLOCATE( VolcanicState%VolcanicSpeciesIndex, STAT=RC ) + CALL CC_CheckVar('VolcanicState%VolcanicSpeciesIndex', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + END IF + + IF (ALLOCATED(VolcanicState%VolcanicSpeciesName)) THEN + DEALLOCATE( VolcanicState%VolcanicSpeciesName, STAT=RC ) + CALL CC_CheckVar('VolcanicState%VolcanicSpeciesName', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + END IF + + IF (ALLOCATED(VolcanicState%EmissionPerSpecies)) THEN + DEALLOCATE( VolcanicState%EmissionPerSpecies, STAT=RC ) + CALL CC_CheckVar('VolcanicState%EmissionPerSpecies', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + END IF + + IF (ALLOCATED(VolcanicState%TotalEmission)) THEN + DEALLOCATE( VolcanicState%TotalEmission, STAT=RC ) + CALL CC_CheckVar('VolcanicState%TotalEmission', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + END IF + + end subroutine CCPr_Volcanic_Finalize + + +END MODULE CCPR_Volcanic_Mod diff --git a/src/process/volcanic/CMakeLists.txt b/src/process/volcanic/CMakeLists.txt new file mode 100644 index 00000000..c650c7f2 --- /dev/null +++ b/src/process/volcanic/CMakeLists.txt @@ -0,0 +1,12 @@ +set(_srcs CCPr_Volcanic_Mod.F90 ccpr_scheme_volcanic_gocart_mod.F90) + +set(_lib CATChem_process_volcanic) + +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/volcanic/ccpr_scheme_volcanic_gocart_mod.F90 b/src/process/volcanic/ccpr_scheme_volcanic_gocart_mod.F90 new file mode 100644 index 00000000..e800beab --- /dev/null +++ b/src/process/volcanic/ccpr_scheme_volcanic_gocart_mod.F90 @@ -0,0 +1,299 @@ +!> +!! \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 and Wei Li +!! \date 07/2024 +!!!> +module CCPr_Scheme_Volcanic_GOCART_Mod + + implicit none + + private + + public :: CCPr_Scheme_Volcanic_GOCART + public :: VolcanicEmisData + public :: ReadASCIIPointEmis + + !> \brief VolcanicEmissionData + !! + !! VolcanicEmissionData is to hold volcanic emission data. + !! + !! \param vlat Volcano latitude + !! \param vlon Volcano longitude + !! \param VEmis Volcanic emissions [kg S/s] + !! \param vbase Bottom elevation of emissions [m] + !! \param vtop Top elevation of emissions [m] + !! \param nPts Number of volcanic sources in the current file + !! \param emissfile Emissions file name + !! \param label Label for emissions + !! + !! \ingroup CCPr_Scheme_Volcanic_GOCART_Mod + !!!> + type :: VolcanicEmisData + real :: vlat !volcano latitude + real :: vlon !volcano longitude + real :: VEmis !volcanic emissions [kg S/s] + integer :: vbase !bottom elevation of emissions [m] + integer :: vtop !top elevation of emissions [m] + integer :: nPts !number of volcanic sources in the current file + character(len=255) :: emissfile !emissions file name + character(len=255) :: label !label for emissions + end type VolcanicEmisData + +contains + + !> \brief Brief description of the subroutine + !! + !! \param km Number of vertical levels + !! \param cdt Model timestep [sec] + !! \param VStart Emissions Start time [sec] + !! \param VEnd Emissions end time [sec] + !! \param nVolc Number of volcanic sources + !! \param iPoint Grid cell index i of each volcanic source + !! \param jPoint Grid cell index j of each volcanic source + !! \param hms Current model time [sec] + !! \param g0 Gravity [m/s^2] + !! \param zbox Geopotential Height difference [m] for layer + !! \param delp Pressure Thickness for layer [Pa] + !! \param area Area of grid cell [m^2] + !! \param vSO2 Volcanic emissions [kg S/s] + !! \param nSO2 Index of SO2 relative to other sulfate tracers + !! \param SO2 SO2 emissions [kg kg-1] + !! \param vCloud Top elevation of emissions [m] + !! \param vElev Bottom elevation of emissions [m] + !! \param vLat Latitude specified in file [degree] + !! \param VLon Longitude specified in file [degree] + !! \param RC Success or Failure + !! + !!!> + subroutine CCPr_Scheme_Volcanic_GOCART(km, & + cdt, & + VStart, & + VEnd, & + nVolc, & + iPoint, & + jPoint, & + hms, & + g0, & + zbox, & + delp, & + area, & + vSO2, & + nSO2, & + SO2, & + vCloud, & + vElev, & + vLat, & + VLon, & + RC ) + + USE GOCART2G_Process, only: SUVolcanicEmissions + + IMPLICIT NONE + + ! Arguments + INTEGER, intent(in) :: km ! number of vertical levels + REAL, intent(in) :: cdt ! model timestep [sec] + INTEGER, intent(inout),dimension(:) :: vStart ! Emissions Start time [sec] + INTEGER, intent(inout),dimension(:) :: vEnd ! Emissions end time [sec] + INTEGER, intent(inout) :: nVolc ! number of volcanic sources + INTEGER, intent(inout),dimension(:) :: iPoint, jPoint ! grid cell index of each volcanic source + INTEGER, intent(in) :: hms ! current model time [sec] + REAL, intent(in) :: g0 + REAL, intent(in), dimension(:) :: zbox ! geopotential Height difference [m] for layer + REAL, intent(in), dimension(:) :: delp ! Pressure Thickness for layer [Pa] + REAL, intent(inout),dimension(:,:) :: area ! area of grid cell [m^2] + REAL, intent(inout),dimension(:) :: vSO2 ! volcanic emissions [kg S/s] + INTEGER, intent(in) :: nSO2 ! index of SO2 relative to other sulfate tracers + REAL, intent(inout),dimension(:,:,:),pointer :: SO2 ! SO2 emissions [kg kg-1] + 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] + INTEGER, intent(inout) :: rc ! error code + + !local variables + REAL, dimension(:,:,:),pointer :: SU_emis ! SU emissions [kg/m2/s] + REAL, dimension(:,:),pointer :: SO2EMVol ! volcanic emissions [kg m-2 s-1] + REAL, parameter :: fMassSulfur = 32. ! gram molecular weights of species + REAL, parameter :: fMassSO2 = 64. ! gram molecular weights of species + 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 + call PrepMetVarsForGOCARTSUV(km, & + delp, & + zbox, & + GOCART_DELP, & + GOCART_ZBOX) + + !convert SO2 unit from Kg S to Kg SO2 + vSO2 = vSO2 * fMassSO2 / fMassSulfur + + !call gocart emission function + if (nVolc > 0) then + + !iPoint(1) = 0 + !jPoint(1) = 0 + + allocate(SU_emis(1,1,nSO2)) !TODO: nSO2 =1 for now + allocate(SO2EMVol, mold=area) + !allocate(SO2EMVN, SO2EMVE, mold=area) + + call SUvolcanicEmissions (nVolc, vStart, vEnd, vSO2, & + vElev, vCloud, & + iPoint, jPoint, & + hms, SO2EMVol, 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) + if (associated(SU_emis)) nullify(SU_emis) + if (associated(SO2EMVol)) nullify(SO2EMVol) + + + end subroutine CCPr_Scheme_Volcanic_GOCART + + + !> \brief Brief description of the subroutine + !! + !! \param filename Emissions file name + !! \param label Label for emissions + !! \param VolcanicEmissions Volcanic emissions data + !! \param rc Success or Failure + !! + !!!> + subroutine ReadASCIIPointEmis (filename, label, VolcanicEmissions, rc ) + + implicit none + + character(len=1055), intent(in) :: filename + character(len=7), intent(in) :: label + type(VolcanicEmisData), intent(inout), allocatable :: VolcanicEmissions(:) + integer, intent(inout) :: rc + !local variables + integer :: num_emiss_sources=0 + integer :: num_lines=0 + integer :: num_skip=0 + integer :: i + character(1056) :: line + character(len=255) :: errmsg + + ! 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 (trim(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( VolcanicEmissions(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) VolcanicEmissions(i)%vlat, & + VolcanicEmissions(i)%vlon, & + VolcanicEmissions(i)%vemis, & + VolcanicEmissions(i)%vbase, & + VolcanicEmissions(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) + + VolcanicEmissions%nPts = num_emiss_sources + VolcanicEmissions%emissfile = trim(filename) + VolcanicEmissions%label = trim(label) + + end subroutine ReadASCIIPointEmis + + + !> \brief Brief description of the subroutine + !! + !! \param km Number of vertical levels + !! \param delp Pressure Thickness for layer [Pa] + !! \param zbox Geopotential Height difference [m] for layer + !! \param GOCART_DELP Pressure Thickness for layer in GOCART format [Pa] + !! \param GOCART_ZBOX Geopotential Height difference in GOCART format [m] for layer + !! + !!!> + 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 ! Pressure Thickness for layer [Pa] + REAL, intent(in), DIMENSION(:), target :: zbox ! Geopotential Height difference [m] for layer + + ! INPUT/OUTPUTS + REAL, intent(inout), pointer :: GOCART_DELP(:,:,:) !< pressure thickness for layer in GOCART format [Pa] + REAL, intent(inout), pointer, DIMENSION(:,:,:) :: GOCART_ZBOX !< Geopotential Height difference in GOCART format [m] for layer + + 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 + + +end module CCPr_Scheme_Volcanic_GOCART_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 487163ec..ac393673 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -76,6 +76,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 @@ -87,6 +88,22 @@ add_test( WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} ) +add_executable(test_volcanic test_volcanic.F90) +target_link_libraries(test_volcanic PRIVATE CATChem_core) +target_link_libraries(test_volcanic PRIVATE CATChem) +target_link_libraries(test_volcanic PRIVATE CATChem_process_volcanic) +target_link_libraries(test_volcanic PRIVATE testing) +set_target_properties( + test_volcanic + PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR}/include +) + +add_test( + NAME test_volcanic + COMMAND test_volcanic + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} +) + add_executable(test_chem test_chem.f90) target_link_libraries(test_chem PRIVATE CATChem) target_link_libraries(test_chem PRIVATE testing) diff --git a/tests/Configs/Default/CATChem_config.yml b/tests/Configs/Default/CATChem_config.yml index ca532a07..4c9b2c22 100644 --- a/tests/Configs/Default/CATChem_config.yml +++ b/tests/Configs/Default/CATChem_config.yml @@ -47,3 +47,7 @@ process: activate: true scheme_opt: 1 resuspension: false + volcanic: + activate: true + scheme_opt: 1 + filedir: ./Emissions/so2_volcanic_emissions_Carns diff --git a/tests/Configs/Default/CATChem_emission.yml b/tests/Configs/Default/CATChem_emission.yml index aee12410..78135377 100644 --- a/tests/Configs/Default/CATChem_emission.yml +++ b/tests/Configs/Default/CATChem_emission.yml @@ -22,3 +22,8 @@ dust: DU10: long_name: "Dust Species Coarse Mode" map: ["dust1", "dust2"] +VOLCANIC: + SO2: + long_name: "Sulfur Dioxide" + scale: [1.0] + map: ["SO2"] 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_volcanic.F90 b/tests/test_volcanic.F90 new file mode 100644 index 00000000..1d9387f0 --- /dev/null +++ b/tests/test_volcanic.F90 @@ -0,0 +1,170 @@ +program test_volcanic + 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(VolcanicStateType) :: VolcanicState + type(GridStateType) :: GridState + type(EmisStateType) :: EmisState + + ! Integers + INTEGER:: rc ! Success or failure + + character(len=:), allocatable :: title + integer :: c ,s ! Loop counter for emission state + + ! Error handling + CHARACTER(LEN=512) :: errMsg + CHARACTER(LEN=255) :: thisLoc + CHARACTER(LEN=255), PARAMETER :: configFile ='Configs/Default/CATChem_config.yml' + + + thisLoc = 'test_volcanic -> 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' + VolcanicState%Activate = .false. + call print_info(Config, VolcanicState, MetState, title) + write (*,*) '-- ' + write (*,*) 'Completed ', title + write (*,*) '--' + + !allocate emission state + if (EmisState%nCats > 0) then + do c = 1, EmisState%nCats + do s = 1, EmisState%Cats(c)%nSpecies + ALLOCATE(EmisState%Cats(c)%Species(s)%Flux(GridState%number_of_levels), STAT=RC) + if (RC /= CC_SUCCESS) then + ErrMsg = 'Error allocating "EmisState%Cats%Species%Flux"!' + call cc_emit_error(ErrMsg, RC, ThisLoc) + stop 1 !!Note here is not 'return' + endif + end do + end do + end if + + !---------------------------- + ! Test 2 + !---------------------------- + ! Set number of Volcanic species + + !ChemState%nSpeciesVolcanic = 2 + !VolcanicState%Activate = .true. + + ! Meteorological State + MetState%YMD = 20220101 + MetState%HMS = 120000 + MetState%TSTEP = 300 + MetState%AREA_M2 = 10000.0_fp + MetState%NLEVS = 8 + allocate(MetState%BXHEIGHT(MetState%NLEVS)) + + !TODO: is the layer index reversed in the GOCART? + allocate(MetState%DELP(MetState%NLEVS)) + MetState%DELP(1:MetState%NLEVS)= (/25000, 20000, 15000, 110000, 10000, 9000,5000,4000/)! Need to change to something more reasonable and check units. + MetState%BXHEIGHT(1:MetState%NLEVS) =(/10000, 8000, 7000, 5000, 3000, 1000, 100, 50/) ! temporary, change to something more reasonable and check units + + !VolcanicState%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 = "Volcanic Test 2 | Test GOCART Volcanic defaults" + Config%volcanic_activate = .TRUE. + Config%volcanic_scheme = 1 + + call cc_volcanic_init(Config, VolcanicState, EmisState, rc) + if (rc /= CC_SUCCESS) then + errMsg = 'Error in cc_volcanic_init' + call cc_emit_error(errMsg, rc, thisLoc) + stop 1 + end if + + call cc_volcanic_run(MetState, VolcanicState, EmisState, rc) + if (rc /= CC_SUCCESS) then + errMsg = 'Error in _volcanicemissions_run' + call cc_emit_error(errMsg, rc, thisLoc) + stop 1 + end if + + !TODO: This is specific to the inputs in this test only. Change it when inputs are changed. + call assert( sum(VolcanicState%TotalEmission) > 1.0e-2_fp, "Test non-zero Sulfur Volcanic Emissions") + call assert( rae(sum(VolcanicState%TotalEmission(1:3)), 0.0_fp) , "Test1 zero Sulfur Volcanic Emissions") + call assert( rae(sum(VolcanicState%TotalEmission(8:28)), 0.0_fp), "Test2 zero Sulfur Volcanic Emissions") + call print_info(Config, VolcanicState, MetState, title) + + call cc_volcanic_finalize( VolcanicState, rc) + if (rc /= CC_SUCCESS) then + errMsg = 'Error in _volcanic_finalize' + call cc_emit_error(errMsg, rc, thisLoc) + stop 1 + end if + +contains + + subroutine print_info(Config_, VolcanicState_, MetState_, title_) + + type(ConfigType), intent(in) :: Config_ + type(MetStateType), intent(in) :: MetState_ + type(VolcanicStateType), intent(in) :: VolcanicState_ + character(len=*), intent(in) :: title_ + + write(*,*) '=======================================' + write(*,*) title_ + write(*,*) '=======================================' + write(*,*) '*************' + write(*,*) 'Configuration ' + write(*,*) '*************' + write(*,*) 'Config%volcanic_activate = ', Config_%volcanic_activate + write(*,*) 'Config%volcanic_scheme = ', Config_%volcanic_scheme + + + if (VolcanicState_%Activate) then + + write(*,*) 'VolcanicState%Activate = ', VolcanicState_%Activate + write(*,*) 'VolcanicState%SchemeOpt = ', VolcanicState_%SchemeOpt + write(*,*) 'MetState%DELP =', MetState_%DELP + write(*,*) 'MetState%BXHEIGHT = ', MetState_%BXHEIGHT + write(*,*) 'VolcanicState%CatIndex = ', VolcanicState_%CatIndex + write(*,*) 'VolcanicState%nVolcanicSpecies = ', VolcanicState_%nVolcanicSpecies + write(*,*) 'VolcanicState%VolcanicSpeciesName = ', VolcanicState_%VolcanicSpeciesName + write(*,*) 'VolcanicState%EmissionPerSpecies = ', VolcanicState_%EmissionPerSpecies + write(*,*) 'VolcanicState%TotalEmission = ', VolcanicState_%TotalEmission + + end if + + end subroutine print_info + + +end program test_volcanic diff --git a/util/metutils_mod.F90 b/util/metutils_mod.F90 deleted file mode 100644 index 2c1584bd..00000000 --- a/util/metutils_mod.F90 +++ /dev/null @@ -1,186 +0,0 @@ -@ -0,0 +1,188 @@ -MODULE MetUtils_Mod - USE Precision_Mod - USE Constants - USE Error_Mod, ONLY : CC_SUCCESS - - IMPLICIT NONE - - PRIVATE - -CONTAINS - - !> - !! \brief Calculate the weibull distribution for 10m wind speed (u10, v10) - !! - !! The Weibull distribution correction ends up being a multiplicative constant - !! (g) times our present source function (see Eq. 12 in Fan & Toon, 2011 and notes for - !! (9/22/11). This constant is derived from the incomplete and complete forms of the gamma - !! function, hence the utilities pasted below. The Weibull function and shape - !! parameters (k, c) assumed are from Justus 1978. - !! - !! \param[inout] gweibull - !! \param[in] weibullFlag - !! \param[in] u10 - !! \param[in] v10 - !! \param[out] RC - !!!> - subroutine weibullDistribution(gweibull, weibullFlag, u10, v10, RC) - - implicit none - - ! Input/Output - !------------- - real(fp), intent(inout) :: gweibull - - ! Input - !------ - integer, intent(in) :: weibullFlag - real(fp), intent(in) :: u10 - real(fp), intent(in) :: v10 - errMsg = '' - thisLoc = ' -> at weibullDistribution (in util/metutils_mod.F90)' - - ! Output - !------- - integer, intent(out) :: RC - - ! Local Variables - real(kind=f8) :: a, c, k, wt, x - real(kind=f8) :: wm - integer :: i, j - - integer :: status - - ! Initialize - RC = CC_SUCCESS - gweibull = 1.0_fp - - wm = sqrt(u10 ** 2 + v10 ** 2) ! Mean wind speed - wt = 4.d0 - - if (weibullFlag) then - gweibull = 0.0_fp - - if (wm > 0.01) then - k = 0.94d0 * sqrt(wm) - c = wm / gamma(1.d0 + 1.d0 / k) - x = (wt / c) ** k - a = 3.41d0 / k + 1.d0 - gweibull = (c / wm) ** 3.41d0 * igamma(a,x, RC) - endif - endif - - if (RC /= CC_SUCCESS) then - errMsg = 'Error Calculating Weibull Distribution' - call CC_Error(errMsg, RC, thisLoc) - return - endif - - end subroutine weibullDistribution - - !> - !! \brief Calculate the incomplete Gamma function - !! - !! The incomplete Gamma function is defined as - !! \int_x^\infty t^{A-1}\exp(-t) dt - !! - !! \param[in] A - !! \param[in] X - !! \param[out] RC - DOUBLE PRECISION function igamma(A, X, rc) -!----------------------------------------------------------------------- -! incomplete (upper) Gamma function -! \int_x^\infty t^{A-1}\exp(-t) dt -!----------------------------------------------------------------------- - IMPLICIT NONE - double precision, intent(in) :: A - DOUBLE PRECISION, INTENT(IN) :: X - - integer, intent(out) :: rc -! LOCAL VARIABLE - DOUBLE PRECISION :: XAM, GIN, S, R, T0 - INTEGER K - rc = __SUCCESS__ - - XAM=-X+A*LOG(X) - IF (XAM.GT.700.0.OR.A.GT.170.0) THEN - WRITE(*,*)'IGAMMA: a and/or x too large, X = ', X - WRITE(*,*) 'A = ', A - rc = __FAIL__ - return - ENDIF - - IF (X.EQ.0.0) THEN - IGAMMA=GAMMA(A) - - ELSE IF (X.LE.1.0+A) THEN - S=1.0/A - R=S - DO K=1,60 - R=R*X/(A+K) - S=S+R - IF (ABS(R/S).LT.1.0e-15) EXIT - END DO - GIN=EXP(XAM)*S - IGAMMA=GAMMA(A)-GIN - ELSE IF (X.GT.1.0+A) THEN - T0=0.0 - DO K=60,1,-1 - T0=(K-A)/(1.0+K/(X+T0)) - end do - - IGAMMA=EXP(XAM)/(X+T0) - - ENDIF - - end function igamma - - !> - !! \brief Calculate the friction velocity - !! - !! \param[in] U10 - !! \param[in] V10 - !! \param[in] z - !! \param[in] z0 - !! \param[out] USTAR - !! \param[out] RC - !!!> - SUBROUTINE Calc_USTAR(U10, V10, z, z0, USTAR, RC) - - IMPLICIT NONE - - ! Input - !------ - real(fp), intent(in) :: U10 ! Wind Speed (m/s) - real(fp), intent(in) :: V10 ! Wind Speed (m/s) - real(fp), intent(in) :: z ! Reference Height for wind speed [m] (typically 10m) - real(fp), intent(in) :: z0 ! Surface roughness [m] - - ! Output - !------- - real(fp), intent(out) :: USTAR - integer, intent(out) :: RC - - ! Local Variables - !---------------- - CHARACTER(LEN=512) :: errMsg, thisLoc - real(fp) :: wind_speed ! wind speed - errMsg = '' - thisLoc = ' -> at Calc_USTAR (in util/metutils_mod.F90)' - - ! Initialize - RC = CC_SUCCESS - USTAR = 0.0_fp - - ! Calculate Mean wind speed from u10 and v10 - !------------------------------------------- - wind_speed = sqrt(U10**2 + V10**2) - - ! Calculate USTAR - !---------------- - USTAR = VON_KARMAN / log (z / z0) * wind_speed - - RETURN - - END SUBROUTINE Calc_USTAR