diff --git a/src/api/CMakeLists.txt b/src/api/CMakeLists.txt index 2a338f40..8aff5b95 100644 --- a/src/api/CMakeLists.txt +++ b/src/api/CMakeLists.txt @@ -8,6 +8,8 @@ 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_DMS) 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..7ea3bb2e 100644 --- a/src/api/catchem.F90 +++ b/src/api/catchem.F90 @@ -82,6 +82,16 @@ 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 + ! DMS + use CCPr_DMS_mod, only: DMSStateType !< DMS State + use CCPr_DMS_mod, only: cc_dms_init => CCPr_DMS_Init !< DMS Process Initialization Routine + use CCPr_DMS_mod, only: cc_dms_run => CCPr_DMS_Run !< DMS Process Run Routine + use CCPr_DMS_mod, only: cc_dms_finalize => CCPr_DMS_Finalize !< DMS 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..7bc3725c 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 !!!> @@ -362,6 +363,11 @@ subroutine FindSpecByName(ChemState, name, index, RC) exit endif enddo + if (index == 0) then + RC = CC_FAILURE + ErrMsg = 'Species not found: ' // TRIM(name) + call CC_Warning(ErrMsg, RC, thisLoc) + endif end subroutine FindSpecByName diff --git a/src/core/config_mod.F90 b/src/core/config_mod.F90 index d66b1968..36db3653 100644 --- a/src/core/config_mod.F90 +++ b/src/core/config_mod.F90 @@ -173,6 +173,24 @@ SUBROUTINE Read_Input_File( Config , GridState, EmisState, ChemState, RC, Config RETURN ENDIF + call Config_Process_DMS(ConfigInput, Config, RC) + IF ( RC /= CC_SUCCESS ) THEN + errMsg = 'Error in "Config_Process_DMS"!' + CALL CC_Error( errMsg, RC, thisLoc ) + CALL QFYAML_CleanUp( ConfigInput ) + CALL QFYAML_CleanUp( ConfigAnchored ) + 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 +1391,165 @@ 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 + + + !> \brief Process DMS configuration + !! + !! This function processes the DMS 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_DMS( 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_DMS begins here! + !======================================================================== + + ! Initialize + RC = CC_SUCCESS + thisLoc = ' -> at Config_Process_DMS (in CATChem/src/core/config_mod.F90)' + errMsg = '' + + ! TODO #105 Fix reading of config file + key = "process%DMS%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%DMS_activate = v_bool + + + key = "process%DMS%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%DMS_scheme = v_int + + write(*,*) "DMS Configuration" + write(*,*) '------------------------------------' + write(*,*) 'Config%DMS_activate = ', Config%DMS_activate + write(*,*) 'Config%DMS_scheme = ', Config%DMS_scheme + write(*,*) '------------------------------------' + + + END SUBROUTINE Config_Process_DMS + + END MODULE config_mod diff --git a/src/core/config_opt_mod.F90 b/src/core/config_opt_mod.F90 index cc62f873..a6d59096 100644 --- a/src/core/config_opt_mod.F90 +++ b/src/core/config_opt_mod.F90 @@ -45,6 +45,11 @@ 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 + !! - `DMS_activate` : Activate DMS emissions + !! - `DMS_scheme` : Scheme option DMS process !! !! \ingroup core_modules !!!> @@ -102,6 +107,15 @@ 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 + + ! DMS Process + LOGICAL :: DMS_activate + INTEGER :: DMS_scheme + END TYPE ConfigType CONTAINS @@ -176,6 +190,15 @@ 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 = ' ' + + ! DMS Process + Config%DMS_activate = .FALSE. + Config%DMS_scheme = 1 + END SUBROUTINE Set_Config !> \brief Cleanup the Config options !! diff --git a/src/core/emisstate_mod.F90 b/src/core/emisstate_mod.F90 index ec5bcdd7..9041d7ad 100644 --- a/src/core/emisstate_mod.F90 +++ b/src/core/emisstate_mod.F90 @@ -358,7 +358,6 @@ subroutine Emis_Find_Chem_Map_Index(EmisState, ChemState, RC) return endif EmisState%Cats(c)%Species(s)%EmisMapIndex(n) = index - EmisState%Cats(c)%Species(s)%EmisMapIndex(n) = index enddo enddo enddo diff --git a/src/core/metstate_mod.F90 b/src/core/metstate_mod.F90 index 164ccdb1..ee8621e8 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 !--------- @@ -192,6 +194,9 @@ MODULE MetState_Mod REAL(fp), ALLOCATABLE :: PMID(:) !< Average wet air pressure [hPa] defined as arithmetic average of edge pressures REAL(fp), ALLOCATABLE :: PMID_DRY(:) !< Dry air partial pressure [hPa] defined as arithmetic avg of edge pressures + !TODO: temporary put DMS concentration here for DMS emissions from gocart; may read from ChemState in the future + REAL(fp) :: DMSO_CONC !< DMS concentration [mol/L] + END TYPE MetStateType CONTAINS diff --git a/src/core/qfyaml_mod.F90 b/src/core/qfyaml_mod.F90 index a4b50cd5..6a6653b1 100644 --- a/src/core/qfyaml_mod.F90 +++ b/src/core/qfyaml_mod.F90 @@ -1159,15 +1159,9 @@ SUBROUTINE QFYAML_Read_Emis_File( yml, fileName, yml_anchored, EmisState, RC ) cspecies = TRIM(yml%vars(n)%category(ix+1:)) tmpString = TRIM(yml%vars(n)%category(ix+1:)) names = [names, tmpString] + j = j + 1 ! increment species index endif endif - - if ((n == yml%num_vars) .and. (TRIM(yml%vars(n)%category(1:ix-1)) == TRIM(current))) then - cspecies = TRIM(yml%vars(n)%category(ix+1:)) - tmpString = TRIM(yml%vars(n)%category(ix+1:)) - names = [names, tmpString] - j = j + 1 - endif enddo ! Fill last category species names diff --git a/src/process/CMakeLists.txt b/src/process/CMakeLists.txt index a563f0b7..86a8357e 100644 --- a/src/process/CMakeLists.txt +++ b/src/process/CMakeLists.txt @@ -4,4 +4,6 @@ add_subdirectory(dust) add_subdirectory(seasalt) add_subdirectory(plumerise) add_subdirectory(drydep) +add_subdirectory(dms) +add_subdirectory(volcanic) add_subdirectory(chem) diff --git a/src/process/dms/CCPr_DMS_Mod.F90 b/src/process/dms/CCPr_DMS_Mod.F90 new file mode 100644 index 00000000..0957ee75 --- /dev/null +++ b/src/process/dms/CCPr_DMS_Mod.F90 @@ -0,0 +1,307 @@ +!> \brief Driver for CATChem DMS process +!! +!!\defgroup catchem_dms_process +!! The CATChem DMS Process group holds all the CATCHem DMS processes. +!! +!! \author Lacey Holland and Wei Li +!! \date 01/2025 +!!!> +MODULE CCPR_DMS_mod + USE Precision_mod + USE Error_Mod + USE DiagState_Mod, Only : DiagStateType + USE MetState_Mod, Only : MetStateType + USE ChemState_Mod, Only : ChemStateType + USE EmisState_Mod, Only : EmisStateType + USE GridState_Mod, Only : GridStateType + USE Config_Opt_Mod, Only : ConfigType + + IMPLICIT NONE + + PRIVATE + + PUBLIC :: CCPR_DMS_Init + PUBLIC :: CCPR_DMS_Run + PUBLIC :: CCPR_DMS_Finalize + PUBLIC :: DMSStateType + + !> \brief DMSStateType + !! + !! \details Contains all the information needed to run DMS Process + !! + !! This type contains the following variables: + !! \param Activate Activate Process (True/False) + !! \param SchemeOpt Scheme Option + !! \param nDMSSpecies Number of DMS species + !! \param DMSSpeciesIndex Index of DMS species + !! \param DMSSpeciesName Name of DMS species + !! \param SpcIDs CATChem species IDs + !! \param CatIndex Index of emission category in EmisState + !! \param TotalEmission Total emission [kg/m^2/s] + !! \param EmissionPerSpecies Emission per species [kg/m^2/s] + !! + !! \ingroup catchem_dms_process + !!!> + + TYPE :: DMSStateType + LOGICAL :: Activate ! Activate Process (True/False) + INTEGER :: SchemeOpt ! SchemeOption + integer :: nDMSSpecies !< Number of DMS species + integer, allocatable :: DMSSpeciesIndex(:) !< Index of DMS species + character(len=31), allocatable :: DMSSpeciesName(:) !< name of DMS species + integer, allocatable :: SpcIDs(:) !< CATChem species IDs + integer :: CatIndex !< Index of emission category in EmisState + + ! Process Specific Parameters + real(fp) :: TotalEmission !< Total emission [kg/m^2/s] + real(fp), allocatable :: EmissionPerSpecies(:) !< Emission per species [kg/m^2/s] + + END TYPE DMSStateType + +CONTAINS + + !> + !! \brief Initialize the CATChem DMS process + !! + !! \param Config CATCHem configuration options + !! \param DMSState CATCHem DMS state + !! \param ChemState CATCHem chemical state + !! \param RC Error return code + !! + !!!> + SUBROUTINE CCPR_DMS_Init( Config, DMSState, EmisState, RC ) + ! USE + + IMPLICIT NONE + ! INPUT PARAMETERS + !----------------- + TYPE(ConfigType) :: Config ! Module options + TYPE(EmisStateType) :: EmisState ! Chemical state + + ! INPUT/OUTPUT PARAMETERS + !------------------------ + TYPE(DMSStateType) :: DMSState ! DMS state + INTEGER, INTENT(INOUT) :: RC ! Success or failure + + ! Error handling + !--------------- + CHARACTER(LEN=255) :: ErrMsg + CHARACTER(LEN=255) :: ThisLoc + + ! LOCAL VARIABLES + !---------------- + INTEGER :: c + + !================================================================= + ! CCPR_DMS_Init begins here! + !================================================================= + ErrMsg = '' + ThisLoc = ' -> at CCPR_DMS_INIT (in process/DMS/ccpr_dms_mod.F90)' + + ! First check if process is activated in config + if (Config%DMS_activate) then + + ! Activate Process + !------------------ + DMSState%Activate = .true. + + ! Set scheme option + !------------------ + DMSState%SchemeOpt = config%DMS_Scheme + + !Find DMS caterory index in EmisState for future use + !-------------------------------------------- + do c = 1, EmisState%nCats + if (EmisState%Cats(c)%name == 'DMSO') then + DMSState%CatIndex = c + exit + endif + end do + + ! Set number of species from EmisState + !---------------------- + DMSState%nDMSSpecies = EmisState%Cats(DMSState%CatIndex)%nSpecies + + !------------------------------------ + ! Allocate emission species index + ALLOCATE( DMSState%DMSSpeciesIndex(DMSState%nDMSSpecies), STAT=RC ) + CALL CC_CheckVar('DMSState%DMSSpeciesIndex', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + DMSState%DMSSpeciesIndex = -1 + + ! Allocate emission speceis names + ALLOCATE( DMSState%DMSSpeciesName(DMSState%nDMSSpecies), STAT=RC ) + CALL CC_CheckVar('DMSState%DMSSpeciesName', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + DMSState%DMSSpeciesName = '' + + ! Allocate CatChem species index + ALLOCATE( DMSState%SpcIDs(DMSState%nDMSSpecies), STAT=RC ) + CALL CC_CheckVar('DMSState%SpcIDs', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + DMSState%SpcIDs = -1 + + ! Allocate emission flux + ALLOCATE( DMSState%EmissionPerSpecies(DMSState%nDMSSpecies), STAT=RC ) + CALL CC_CheckVar('DMSState%EmissionPerSpecies', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + DMSState%EmissionPerSpecies = ZERO + + !initialize total emissions + DMSState%TotalEmission = ZERO + + else + + DMSState%Activate = .false. + + endif + + end subroutine CCPR_DMS_Init + + !> + !! \brief Run the DMS emission scheme + !! + !! \param [IN] MetState - The MetState object + !! \param [INOUT] DMSState - The DMSState object + !! \param [INOUT] EmisState - The EmisState object + !! \param [OUT] RC Return code + !!!> + SUBROUTINE CCPr_DMS_Run( MetState, DMSState, EmisState, ChemState, RC ) + + ! USE + USE constants, only : g0 + USE CCPr_Scheme_GOCART_DMS_Mod, only : CCPR_Scheme_GOCART_DMS + + IMPLICIT NONE + + ! INPUT PARAMETERS + TYPE(MetStateType), INTENT(IN) :: MetState ! MetState Instance + + ! INPUT/OUTPUT PARAMETERS + TYPE(DMSStateType), INTENT(INOUT) :: DMSState ! DMSState Instance + TYPE(EmisStateType), INTENT(INOUT) :: EmisState ! EmisState Instance + TYPE(ChemStateType), INTENT(INOUT) :: ChemState ! ChemState Instance + + ! OUTPUT PARAMETERS + INTEGER, INTENT(OUT) :: RC ! Return Code + + ! LOCAL VARIABLES + INTEGER :: s, chem_ind + REAL(fp) :: fMassDMS ! Molecular weight of DMS in g/mol + CHARACTER(LEN=255) :: ErrMsg, thisLoc + REAL, dimension(:,:,:),pointer :: SU_emis ! DMS emissions in kg/m2/s + + ! Initialize + RC = CC_SUCCESS + errMsg = '' + thisLoc = ' -> at CCPr_DMS_Run (in process/DMS/ccpr_DMS_mod.F90)' + + ! If DMS is activated + !------------------------- + if (DMSState%Activate) then + ! Run the DMS Scheme + !------------------------- + do s = 1, DMSState%nDMSSpecies !only one species for now + + if (DMSState%SchemeOpt == 1) then + + !Dimensions (lon,lat,nSpecies) to be consistent with GOCART; + !we all have one since it is a column model with only DMS emission + allocate(SU_emis(1,1,1)); SU_emis = ZERO + + !get the index of DMS species in the ChemState array (TODO: this assumes one-to-one mapping only) + chem_ind = EmisState%Cats(DMSState%CatIndex)%Species(s)%EmisMapIndex(1) + fMassDMS = ChemState%ChemSpecies(chem_ind)%mw_g !get the molecular weight of DMS + + call CCPr_Scheme_GOCART_DMS(MetState%NLEVS, & + MetState%TSTEP, & + g0, & + fMassDMS, & + MetState%T, & + MetState%U10M, & + MetState%V10M, & + MetState%LWI, & + MetState%DELP, & + MetState%DMSO_CONC, & + SU_emis, & + RC) + + if (RC /= CC_SUCCESS) then + errMsg = 'Error in CCPr_Scheme_GOCART_DMS' + CALL CC_Error( errMsg, RC, thisLoc ) + endif + + !put it back to DMSState + DMSState%DMSSpeciesIndex(s) = s + DMSState%DMSSpeciesName(s) = EmisState%Cats(DMSState%CatIndex)%Species(s)%name + DMSState%EmissionPerSpecies(s) = SU_emis(1,1,1) + DMSState%TotalEmission = DMSState%TotalEmission + DMSState%EmissionPerSpecies(s) + if (associated(SU_emis)) nullify(SU_emis) + + else + errMsg = 'ERROR: Unknown DMS scheme option' + RC = CC_FAILURE + CALL CC_Error( errMsg, RC, thisLoc ) + return + + endif !DMS scheme option + end do !for each species + + endif + + end subroutine CCPr_DMS_Run + + + !> + !! \brief Finalize the DMS emission scheme + !! + !! \param [INOUT] DMSState + !! \param [OUT] RC Return code + !!!> + SUBROUTINE CCPr_DMS_Finalize( DMSState, RC ) + + ! USE + !---- + + IMPLICIT NONE + + ! INPUT/OUTPUT PARAMETERS + TYPE(DMSStateType), INTENT(INOUT) :: DMSState ! DMSState Instance + + ! OUTPUT PARAMETERS + INTEGER, INTENT(OUT) :: RC ! Return Code + + ! LOCAL VARIABLES + CHARACTER(LEN=255) :: ErrMsg, thisLoc + + ! Initialize + RC = CC_SUCCESS + errMsg = '' + thisLoc = ' -> at CCPr_DMS_Finalize (in process/DMS/ccpr_DMS_mod.F90)' + + !Deallocate DMSState + IF (ALLOCATED(DMSState%DMSSpeciesIndex)) THEN + DEALLOCATE(DMSState%DMSSpeciesIndex, STAT=RC) + CALL CC_CheckVar('DMSState%DMSSpeciesIndex', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + ENDIF + IF (ALLOCATED(DMSState%DMSSpeciesName)) THEN + DEALLOCATE(DMSState%DMSSpeciesName, STAT=RC) + CALL CC_CheckVar('DMSState%DMSSpeciesName', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + ENDIF + IF (ALLOCATED(DMSState%SpcIDs)) THEN + DEALLOCATE(DMSState%SpcIDs, STAT=RC) + CALL CC_CheckVar('DMSState%SpcIDs', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + ENDIF + IF (ALLOCATED(DMSState%EmissionPerSpecies)) THEN + DEALLOCATE(DMSState%EmissionPerSpecies, STAT=RC) + CALL CC_CheckVar('DMSState%EmissionPerSpecies', 0, RC) + IF (RC /= CC_SUCCESS) RETURN + ENDIF + + end subroutine CCPr_DMS_Finalize + + +END MODULE CCPR_DMS_Mod diff --git a/src/process/dms/CMakeLists.txt b/src/process/dms/CMakeLists.txt new file mode 100644 index 00000000..1862fbdc --- /dev/null +++ b/src/process/dms/CMakeLists.txt @@ -0,0 +1,12 @@ +set(_srcs CCPr_DMS_Mod.F90 ccpr_scheme_gocart_dms_mod.F90) + +set(_lib CATChem_process_DMS) + +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/dms/ccpr_scheme_gocart_dms_mod.F90 b/src/process/dms/ccpr_scheme_gocart_dms_mod.F90 new file mode 100644 index 00000000..261ddb5b --- /dev/null +++ b/src/process/dms/ccpr_scheme_gocart_dms_mod.F90 @@ -0,0 +1,156 @@ +!> +!! \file +!! \brief CCPr Scheme for DMS +!! +!! 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 01/2025 +!! \ingroup catchem_dms_process +!!!> +module CCPr_Scheme_GOCART_DMS_Mod + USE Precision_mod, only : ZERO + + implicit none + + private + + public :: CCPr_Scheme_GOCART_DMS + +contains + + !> \brief GOCART DMS emission scheme + !! + !! \param km number of vertical levels + !! \param cdt model timestep [sec] + !! \param g0 gravity [m/s2] + !! \param tmpu Temperature [K] + !! \param u10m 10-m u-wind component [m/sec] + !! \param v10m 10-m v-wind component [m/sec] + !! \param lwi orography flag; Land, ocean, ice mask + !! \param delp Pressure Thickness for layer [Pa] + !! \param dmso_conc DMS source concentration [nmol/L] + !! \param SU_emis SU emissions, kg/m2/s + !! \param RC Success or Failure + !! + !! Note that other state types may be required, e.g. one specific to the process group. + !!!> + + subroutine CCPr_Scheme_GOCART_DMS(km, cdt, g0, fMassDMS, tmpu, u10m, v10m, lwi, delp, & + dmso_conc, SU_emis, RC) + + ! Uses + USE GOCART2G_process, only: DMSemission + + IMPLICIT NONE + + ! Arguments + INTEGER, intent(in) :: km ! number of vertical levels (only surface layer is used) + REAL, intent(in) :: cdt ! model timestep [sec] + REAL, intent(in) :: g0 ! gravity [m/s2] + REAL, intent(in) :: fMassDMS ! DMS molecular weight [g/mol] + REAL, dimension(:), intent(in) :: tmpu ! Temperature [K] + REAL, intent(in) :: u10m ! 10-m u-wind component [m/sec] + REAL, intent(in) :: v10m ! 10-m v-wind component [m/sec] + INTEGER, intent(in) :: lwi ! orography flag; Land, ocean, ice mask + REAL, dimension(:), intent(in) :: delp ! Pressure Thickness for layer [Pa] + REAL, intent(in) :: dmso_conc ! DMS source concentration [nmol/L] + REAL, intent(inout),dimension(:,:,:),pointer :: SU_emis ! DMS emissions in kg/m2/s + integer, intent(out) :: RC ! Success or Failure + + ! Local Variables + INTEGER, parameter :: NDMS = 1 ! index of DMS relative to other sulfate tracers in GOCART + !REAL, parameter :: fMassDMS=62. ! DMS molecular weight [g/mol] + character(len=256) :: errMsg + character(len=256) :: thisLoc + + real, pointer :: GOCART_TMPU(:,:,:) + real, pointer :: GOCART_DELP(:,:,:) + real, pointer :: GOCART_LWI(:,:) + real, pointer :: GOCART_U10(:,:) + real, pointer :: GOCART_V10(:,:) + real, pointer :: GOCART_DMSO_CONC(:,:) + REAL, allocatable, dimension(:,:,:) :: DMS ! DMS [kg kg-1] + + ! Initialize + errMsg = '' + thisLoc = ' -> at CCPr_Scheme_GOCART_DMS (in CCPr_Scheme_GOCART_DMS_mod.F90)' + + call INCR_REAL_RANK3(delp, GOCART_DELP) + call INCR_REAL_RANK3(tmpu, GOCART_TMPU) + call INCR_REAL_RANK2(u10m, GOCART_U10) + call INCR_REAL_RANK2(v10m, GOCART_V10) + call INCR_REAL_RANK2(real(LWI), GOCART_LWI) + call INCR_REAL_RANK2(DMSO_CONC, GOCART_DMSO_CONC) + allocate(DMS(1,1,km)); DMS(:,:,:)= ZERO + + call DMSemission (km, cdt, g0, & + GOCART_TMPU, GOCART_U10, & + GOCART_V10, GOCART_LWI, & + GOCART_DELP, fMassDMS, GOCART_DMSO_CONC, & + dms, SU_emis, ndms, rc) + + if (associated(GOCART_TMPU)) nullify(GOCART_TMPU) + if (associated(GOCART_DELP)) nullify(GOCART_DELP) + if (associated(GOCART_U10)) nullify(GOCART_U10) + if (associated(GOCART_V10)) nullify(GOCART_V10) + if (associated(GOCART_LWI)) nullify(GOCART_LWI) + if (associated(GOCART_DMSO_CONC)) nullify(GOCART_DMSO_CONC) + if (allocated(DMS)) deallocate(DMS) + + end subroutine CCPr_Scheme_GOCART_DMS + + + !> \brief some subroutines to convert MET data to GOCART format + !! TODO: these may be used by other GOCART processes too, + !! so they could be moved to a more general location + !! + !! \param ARR input met data + !! \param RESULT output met data for GOCART format + !! + !!!> + SUBROUTINE INCR_REAL_RANK2(ARR, RESULT) + REAL, INTENT(IN), TARGET :: ARR + REAL, INTENT(INOUT), POINTER :: RESULT(:,:) + + ALLOCATE(RESULT(1, 1)) + RESULT(1,1)=ARR + + END SUBROUTINE INCR_REAL_RANK2 + + SUBROUTINE INCR_REAL_RANK3(ARR, RESULT) + REAL, INTENT(IN), TARGET :: ARR(:) + REAL, INTENT(INOUT), POINTER :: RESULT(:,:,:) + + ALLOCATE(RESULT(1, 1, SIZE(ARR, 1))) + RESULT(1,1,:)=ARR + + END SUBROUTINE INCR_REAL_RANK3 + + !Not used for now but may be needed in the future; + !comment out for now to avoid compiler warnings + + ! SUBROUTINE INCR_INT_RANK2(ARR, RESULT) + ! INTEGER, INTENT(IN), TARGET :: ARR + ! INTEGER, INTENT(INOUT), POINTER :: RESULT(:,:) + + ! ALLOCATE(RESULT(1, 1)) + ! RESULT(1,1)=ARR + + ! END SUBROUTINE INCR_INT_RANK2 + + ! SUBROUTINE INCR_INT_RANK3(ARR, RESULT) + ! INTEGER, INTENT(IN), TARGET :: ARR(:) + ! INTEGER, INTENT(INOUT), POINTER :: RESULT(:,:,:) + + ! ALLOCATE(RESULT(1, 1, SIZE(ARR, 1))) + ! RESULT(1,1,:)=ARR + + ! END SUBROUTINE INCR_INT_RANK3 + + +end module CCPr_Scheme_GOCART_DMS_Mod 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..6829ef67 100644 --- a/tests/CATChem_config.yml +++ b/tests/CATChem_config.yml @@ -45,3 +45,9 @@ process: activate: true scheme_opt: 1 resuspension: false + SUVolcanicEmissions: + activate: true + scheme_opt: 1 + DMS: + activate: true + scheme_opt: 1 diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 487163ec..f37a7a9b 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,40 @@ 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_dms test_dms.F90) +target_link_libraries(test_dms PRIVATE CATChem_core) +target_link_libraries(test_dms PRIVATE CATChem) +target_link_libraries(test_dms PRIVATE CATChem_process_DMS) +target_link_libraries(test_dms PRIVATE testing) +set_target_properties( + test_dms + PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR}/include +) + +add_test( + NAME test_dms + COMMAND test_dms + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} +) + +message(STATUS "Fortran module directory: ${CMAKE_BINARY_DIR}/include") + 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..d7f3aa1b 100644 --- a/tests/Configs/Default/CATChem_config.yml +++ b/tests/Configs/Default/CATChem_config.yml @@ -47,3 +47,10 @@ process: activate: true scheme_opt: 1 resuspension: false + volcanic: + activate: true + scheme_opt: 1 + filedir: ./Emissions/so2_volcanic_emissions_Carns + DMS: + activate: true + scheme_opt: 1 diff --git a/tests/Configs/Default/CATChem_emission.yml b/tests/Configs/Default/CATChem_emission.yml index aee12410..fa543045 100644 --- a/tests/Configs/Default/CATChem_emission.yml +++ b/tests/Configs/Default/CATChem_emission.yml @@ -22,3 +22,13 @@ dust: DU10: long_name: "Dust Species Coarse Mode" map: ["dust1", "dust2"] +VOLCANIC: + SO2: + long_name: "Sulfur Dioxide" + scale: [1.0] + map: ["SO2"] +DMSO: + DMS: + long_name: "Dimethyl Sulfide" + scale: [1.0] + map: ["dms"] diff --git a/tests/Configs/Default/CATChem_species.yml b/tests/Configs/Default/CATChem_species.yml index 77172c76..4c1fed2f 100644 --- a/tests/Configs/Default/CATChem_species.yml +++ b/tests/Configs/Default/CATChem_species.yml @@ -25,3 +25,46 @@ dust2: is_dust: true is_aerosol: true is_drydep: true +dms: + name: dms + long_name: 'Dimethyl Sulfide' + description: dms + is_gas: true + is_drydep: true + mw_g: 62.0 +NO: + name: NO + long_name: 'Nitrogen Oxide' + description: NO + is_gas: true + is_drydep: true + mw_g: 30.0 +CO: + name: CO + long_name: 'Carbon Monoxide' + description: CO + is_gas: true + is_drydep: true + mw_g: 28.0 +ISOP: + name: ISOP + long_name: 'Isoprene' + description: ISOP + is_gas: true + is_drydep: true + mw_g: 68.0 +#dummy species for testing +TEST: + name: TEST + long_name: 'Test Species' + description: TEST + is_gas: true + is_drydep: true + mw_g: 100.0 +SO2: + name: SO2 + long_name: 'Sulfur Dioxide' + description: SO2 + is_gas: true + is_drydep: true + mw_g: 64.0 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_dms.F90 b/tests/test_dms.F90 new file mode 100644 index 00000000..22dc4ae5 --- /dev/null +++ b/tests/test_dms.F90 @@ -0,0 +1,179 @@ +program test_DMS + use CATChem, fp => cc_rk + use testing_mod, only: assert + use precision_mod, only: rae + implicit none + + type(ConfigType) :: Config + type(ChemStateType) :: ChemState + type(DMSStateType) :: DMSState + type(MetStateType) :: MetState + type(DiagStateType) :: DiagState + 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_DMS -> 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 = 'DMS Test 1 | Read Config' + DMSState%Activate = .false. + call print_info(Config, DMSState, 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 + + !use Emis_Allocate subroutines directly + call cc_allocate_emisstate(GridState, EmisState, rc) + if (rc /= CC_SUCCESS) then + errMsg = 'Error in cc_allocate_emisstate' + call cc_emit_error(errMsg, rc, thisLoc) + stop 1 + end if + + ! map emission to chemistry state (This is needed because we need DMS MW from ChemState) + call cc_emis_to_chem_map(EmisState, ChemState, rc) + if (rc /= CC_SUCCESS) then + errMsg = 'Error in cc_map_emis_to_chem' + call cc_emit_error(errMsg, rc, thisLoc) + stop 1 + end if + + !---------------------------- + ! Test 2 + !---------------------------- + + ! Meteorological State + MetState%TSTEP = 300 + MetState%NLEVS = 1 + allocate(MetState%T(MetState%NLEVS)) + allocate(MetState%DELP(MetState%NLEVS)) + MetState%DELP(1:MetState%NLEVS)= 5000 ! Need to change to something more reasonable and check units. + MetState%T(1:MetState%NLEVS) = 300 ! temporary, change to something more reasonable and check units + MetState%U10M = 5.0_fp + MetState%V10M = 5.0_fp + MetState%LWI = 0 !gocart OCEAN=0.0, LAND = 1.0, SEA_ICE = 2.0 + MetState%DMSO_CONC = 3.25_fp !DMS ocean concentration [nmol/L];TODO: may read from ChemState in the future + + ! 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 = "DMS Test 2 | Test GOCART DMS defaults" + !--------------------------------------------- + DMSState%Activate = .true. + DMSState%SchemeOpt = 1 + + call cc_dms_init(Config, DMSState, EmisState, rc) + if (rc /= CC_SUCCESS) then + errMsg = 'Error in cc_dms_init' + call cc_emit_error(errMsg, rc, thisLoc) + stop 1 + end if + + call cc_dms_run(MetState, DMSState, EmisState, ChemState, rc) + if (rc /= CC_SUCCESS) then + errMsg = 'Error in _dms_run' + call cc_emit_error(errMsg, rc, thisLoc) + stop 1 + end if + + call assert( DMSState%TotalEmission > 0.0_fp, "Test DMS Emissions") + call print_info(Config, DMSState, MetState, title) + + call cc_dms_finalize( DMSState, rc) + if (rc /= CC_SUCCESS) then + errMsg = 'Error in _dms_finalize' + call cc_emit_error(errMsg, rc, thisLoc) + stop 1 + end if + +contains + + subroutine print_info(Config_, DMSState_, MetState_, title_) + + type(ConfigType), intent(in) :: Config_ + type(MetStateType), intent(in) :: MetState_ + type(DMSStateType), intent(in) :: DMSState_ + character(len=*), intent(in) :: title_ + + write(*,*) '=======================================' + write(*,*) title_ + write(*,*) '=======================================' + write(*,*) '*************' + write(*,*) 'Configuration ' + write(*,*) '*************' + write(*,*) 'Config%dms_activate = ', Config_%dms_activate + write(*,*) 'Config%dms_scheme = ', Config_%dms_scheme + + + if (DMSState_%Activate) then + + write(*,*) 'DMSState%Activate = ', DMSState_%Activate + write(*,*) 'DMSState%SchemeOpt = ', DMSState_%SchemeOpt + write(*,*) 'MetState%DELP =', MetState_%DELP + write(*,*) 'MetState%T = ', MetState_%T + write(*,*) 'MetState%U10M = ', MetState_%U10M + write(*,*) 'MetState%V10M = ', MetState_%V10M + write(*,*) 'MetState%LWI = ', MetState_%LWI + write(*,*) 'DMSState%CatIndex = ', DMSState_%CatIndex + write(*,*) 'DMSState%nDMSSpecies = ', DMSState_%nDMSSpecies + write(*,*) 'DMSState%DMSSpeciesName = ', DMSState_%DMSSpeciesName + write(*,*) 'DMSState%EmissionPerSpecies = ', DMSState_%EmissionPerSpecies + write(*,*) 'DMSState%TotalEmission = ', DMSState_%TotalEmission + + end if + + end subroutine print_info + + +end program test_dms 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