Skip to content

File ProcessDryDepInterface_Mod.F90

File List > drydep > ProcessDryDepInterface_Mod.F90

Go to the documentation of this file

module processdrydepinterface_mod

   ! Core CATChem infrastructure
   use precision_mod, only: fp
   use processinterface_mod, only: processinterface, columnprocessinterface
   use statemanager_mod, only: statemanagertype
   use gridmanager_mod, only: gridmanagertype
   use error_mod, only: cc_success, cc_failure, cc_error, cc_warning, errormanagertype
   use diagnosticmanager_mod, only: diagnosticmanagertype
   use diagnosticinterface_mod, only: diagnosticregistrytype, diagnosticfieldtype, diagnosticdatatype
   use virtualcolumn_mod, only: virtualcolumntype, virtualmettype
   use constants, only: g0, airmw  ! Gravitational acceleration for tendency physics and air molecular weight for unit conversion

   ! Core utilities (leverage existing infrastructure)
   use configmanager_mod, only: configmanagertype
   use chemstate_mod, only: chemstatetype
   use metstate_mod, only: metstatetype

   ! Common utilities - unified configuration
   use drydepcommon_mod, only: drydepprocessconfig

   ! Scheme modules
   use drydepscheme_wesely_mod, only: compute_wesely
   use drydepscheme_gocart_mod, only: compute_gocart
   use drydepscheme_zhang_mod, only: compute_zhang

   implicit none
   private

   public :: processdrydepinterface

   type, extends(columnprocessinterface) :: processdrydepinterface
      private

      ! Unified process configuration (bridges ConfigManager to process-specific config)
      type(DryDepProcessConfig), public :: process_config

      ! Process utilities (leverage core infrastructure)
      type(ChemStateType), pointer :: chem_state => null()
      type(MetStateType), pointer :: met_state => null()
      ! Note: state_manager pointer removed as it was never used

      ! Process-specific diagnostic indices (base class handles storage)
      integer :: diag_drydep_con_per_species_idx = -1
      integer :: diag_drydep_velocity_per_species_idx = -1

      ! Column-level diagnostic storage for interfacing with DiagManager
      ! These are allocated per-column during processing and can be used to
      ! accumulate data for DiagManager updates
      real(fp), allocatable :: column_drydep_con_per_species(:)           ! 1D: species - per column
      real(fp), allocatable :: column_drydep_velocity_per_species(:)           ! 1D: species - per column

      ! Scheme-specific diagnostic storage (shared across all schemes that use them)

   contains
      ! Required ProcessInterface implementations
      procedure :: init => process_init
      procedure :: run => process_run
      procedure :: finalize => process_finalize
      procedure :: parse_process_config => parse_drydep_config

      ! Required ColumnProcessInterface implementations
      procedure :: init_column_processing => init_column_processing
      procedure :: run_column => run_column
      procedure :: finalize_column_processing => finalize_column_processing

      ! ProcessInterface capability registration
      procedure :: get_required_met_fields => get_required_met_fields
      procedure :: get_required_diagnostic_fields => get_required_diagnostic_fields

      ! Public testing interface for scheme manipulation
      procedure :: set_scheme => set_drydep_scheme
      procedure :: get_scheme => get_drydep_scheme

      ! Process-specific implementations (column virtualization)
      procedure, private :: run_active_scheme_column
      procedure, private :: run_wesely_scheme_column
      procedure, private :: run_gocart_scheme_column
      procedure, private :: run_zhang_scheme_column

      ! Diagnostic procedures (override base class method)
      procedure :: register_diagnostics => register_and_allocate_diagnostics
      procedure, private :: register_and_allocate_diagnostics
      procedure, private :: calculate_and_update_diagnostics

   end type processdrydepinterface

contains

   subroutine process_init(this, container, rc)
      class(ProcessDryDepInterface), intent(inout) :: this
      type(StateManagerType), intent(inout) :: container
      integer, intent(out) :: rc

      type(ErrorManagerType), pointer :: error_manager

      rc = cc_success

      ! Get error manager
      error_manager => container%get_error_manager()

      ! Initialize column processing capabilities
      call this%init_column_processing(container, rc)
      if (rc /= cc_success) return

      ! Set process-specific name and info
      this%name = 'drydep'
      this%version = '1.0.0'
      this%description = 'Process for computing dry deposition of gas and aerosol species'

      ! Parse process-specific configuration using unified approach
      call this%parse_process_config(container, error_manager, rc)
      if (rc /= cc_success) then
         return
      end if

      ! Get state pointers from container (needed for species loading)
      this%chem_state => container%get_chem_state_ptr()
      this%met_state => container%get_met_state_ptr()

      ! Load species from ChemState based on is_drydep property
      call this%process_config%load_species_from_chem_state(this%chem_state, error_manager)
      ! Note: Error handling managed by error_manager internally

      ! Map diagnostic species names to indices in the species array
      call this%process_config%map_diagnostic_species_indices(error_manager)

      ! Validate the configuration with StateManager
      call this%process_config%validate(container, error_manager)
      ! Note: validate doesn't return rc, but error_manager tracks errors

      ! Register diagnostics for this process (only if diagnostics enabled)
      call this%register_diagnostics(container, rc)
      if (rc /= cc_success) then
         return
      end if

      ! Mark process as initialized and active
      call this%activate()

   end subroutine process_init

   subroutine process_run(this, container, rc)
      class(ProcessDryDepInterface), intent(inout) :: this
      type(StateManagerType), intent(inout) :: container
      integer, intent(out) :: rc

      rc = cc_success

      ! Check if process is active
      if (.not. this%process_config%is_active) then
         return
      end if

      ! For ColumnProcessInterface processes, the ProcessManager handles column iteration
      ! and calls run_column() for each virtual column. This method is mainly a placeholder
      ! for any global 3D operations that need to happen before/after column processing.

      ! Currently no global 3D operations needed for drydep process
      ! All processing happens in run_column() method

   end subroutine process_run

   subroutine process_finalize(this, rc)
      class(ProcessDryDepInterface), intent(inout) :: this
      integer, intent(out) :: rc

      rc = cc_success

      ! Finalize column processing
      call this%finalize_column_processing(rc)
      if (rc /= cc_success) return

      ! Deallocate diagnostic class members
      if (allocated(this%column_drydep_con_per_species)) deallocate(this%column_drydep_con_per_species)
      if (allocated(this%column_drydep_velocity_per_species)) deallocate(this%column_drydep_velocity_per_species)
      ! Deallocate scheme-specific diagnostic fields (only deallocate unique fields once)

      ! Finalize unified configuration
      call this%process_config%finalize()

   end subroutine process_finalize

   subroutine parse_drydep_config(this, state_manager, error_manager, rc)
      class(ProcessDryDepInterface), intent(inout) :: this
      type(StateManagerType), intent(inout) :: state_manager  ! Changed to inout for function call
      type(ErrorManagerType), intent(inout) :: error_manager
      integer, intent(out) :: rc

      type(ConfigManagerType), pointer :: config_manager

      rc = cc_success

      ! Get configuration manager from state manager
      config_manager => state_manager%get_config_ptr()
      if (.not. associated(config_manager)) then
         call error_manager%report_error(1003, &
            'ConfigManager not available from StateManager', rc)
         return
      end if

      ! Use the unified configuration loader from DryDepCommon_Mod
      ! This handles the complexity of parsing hierarchical YAML into process-specific types
      call this%process_config%load_from_config(config_manager, error_manager)
      ! Note: Error handling managed by error_manager internally

      ! Process is now configured - the unified config contains all scheme-specific settings

   end subroutine parse_drydep_config

   !========================================================================
   ! Column Processing Interface Implementation
   !========================================================================

   subroutine init_column_processing(this, container, rc)
      class(ProcessDryDepInterface), intent(inout) :: this
      type(StateManagerType), intent(inout) :: container
      integer, intent(out) :: rc

      rc = cc_success

      ! Enable column processing and set batch size for optimal performance
      call this%enable_column_processing()
      call this%set_column_batch_size(50)  ! Process 50 columns at a time

      ! Any process-specific column processing setup would go here
      ! For drydep, no additional setup is needed

   end subroutine init_column_processing

   subroutine run_column(this, column, container, rc)
      class(ProcessDryDepInterface), intent(inout) :: this
      type(VirtualColumnType), intent(inout) :: column
      type(StateManagerType), intent(inout) :: container
      integer, intent(out) :: rc

      rc = cc_success

      ! Check if process is active
      if (.not. this%process_config%is_active) return

      ! Delegate to the active scheme for column processing
      call this%run_active_scheme_column(column, rc)

      ! Calculate and update diagnostics if enabled
      if (this%process_config%drydep_config%diagnostics .and. rc == cc_success) then
         call this%calculate_and_update_diagnostics(column, container, rc)
      end if

   end subroutine run_column

   subroutine finalize_column_processing(this, rc)
      class(ProcessDryDepInterface), intent(inout) :: this
      integer, intent(out) :: rc

      rc = cc_success

      ! Clean up column processing
      call this%disable_column_processing()

      ! Any process-specific cleanup would go here

   end subroutine finalize_column_processing

   subroutine run_active_scheme_column(this, column, rc)
      class(ProcessDryDepInterface), intent(inout) :: this
      type(VirtualColumnType), intent(inout) :: column
      integer, intent(out) :: rc

      rc = cc_success

      ! For gas/aerosol differentiated processing, run schemes based on species type
      ! Run gas scheme for gas species
      select case (trim(this%process_config%drydep_config%gas_scheme))
       case ('wesely')
         call this%run_wesely_scheme_column(column, rc)
         if (rc /= cc_success) return
       case default
         rc = cc_failure
         return
      end select

      ! Run aerosol scheme for aerosol species
      select case (trim(this%process_config%drydep_config%aero_scheme))
       case ('gocart')
         call this%run_gocart_scheme_column(column, rc)
       case ('zhang')
         call this%run_zhang_scheme_column(column, rc)
       case default
         rc = cc_failure
      end select

   end subroutine run_active_scheme_column

   subroutine run_wesely_scheme_column(this, column, rc)
      class(ProcessDryDepInterface), intent(inout) :: this
      type(VirtualColumnType), intent(inout) :: column
      integer, intent(out) :: rc

      ! Local variables for scheme calculation
      type(VirtualMetType), pointer :: met => null()  ! Pointer to meteorological data
      ! Meteorological fields
      real(fp), allocatable :: bxheight(:)
      real(fp), allocatable :: cldfrc(:)
      real(fp), allocatable :: frlai(:)
      real(fp), allocatable :: frlanduse(:)
      integer, allocatable :: iland(:)
      logical, allocatable :: isice(:)
      logical, allocatable :: island(:)
      logical, allocatable :: issnow(:)
      real(fp), allocatable :: lat(:)
      real(fp), allocatable :: lon(:)
      character(len=255), allocatable :: lucname(:)
      real(fp), allocatable :: obk(:)
      real(fp), allocatable :: ps(:)
      real(fp), allocatable :: salinity(:)
      real(fp), allocatable :: suncosmid(:)
      real(fp), allocatable :: swgdn(:)
      real(fp), allocatable :: ts(:)
      real(fp), allocatable :: tskin(:)
      real(fp), allocatable :: tstep(:)
      real(fp), allocatable :: ustar(:)
      real(fp), allocatable :: z0(:)
      ! Species properties
      real(fp), allocatable :: species_mw_g(:)
      real(fp), allocatable :: species_dd_f0(:)
      real(fp), allocatable :: species_dd_hstar(:)
      real(fp), allocatable :: species_dd_DvzAerSnow(:)
      real(fp), allocatable :: species_dd_DvzMinVal_snow(:)
      real(fp), allocatable :: species_dd_DvzMinVal_land(:)
      real(fp), allocatable :: species_conc(:,:)
      real(fp), allocatable :: species_tendencies(:,:)
      integer :: n_species, n_levels, i
      integer, allocatable :: species_indices(:)
      real(fp) :: loss_fraction  ! Loss fraction for multiplicative tendencies

      rc = cc_success

      ! Get dimensions from virtual column
      n_levels = 1  ! Surface-only processing

      ! Get drydep species information from process configuration
      n_species = this%process_config%drydep_config%n_species
      if (n_species <= 0) then
         return
      end if

      ! Get species indices directly from configuration (pre-computed)
      allocate(species_indices(n_species))
      species_indices(1:n_species) = this%process_config%drydep_config%species_indices(1:n_species)

      ! Allocate arrays
      allocate(species_conc(1, n_species))
      allocate(species_tendencies(1, n_species))
      ! Allocate meteorological field arrays based on field type and process configuration
      allocate(bxheight(n_levels))  ! Atmospheric field - always n_levels
      allocate(cldfrc(1))  ! Surface field - always scalar



      allocate(isice(1))  ! Surface field - always scalar
      allocate(island(1))  ! Surface field - always scalar
      allocate(issnow(1))  ! Surface field - always scalar
      allocate(lat(1))  ! Surface field - always scalar
      allocate(lon(1))  ! Surface field - always scalar
      allocate(lucname(1))  ! Surface field - always scalar
      allocate(obk(1))  ! Surface field - always scalar
      allocate(ps(1))  ! Surface field - always scalar
      allocate(salinity(1))  ! Surface field - always scalar
      allocate(suncosmid(1))  ! Surface field - always scalar
      allocate(swgdn(1))  ! Surface field - always scalar
      allocate(ts(1))  ! Surface field - always scalar
      allocate(tskin(1))  ! Surface field - always scalar
      allocate(tstep(1))  ! Special timestep field - scalar
      allocate(ustar(1))  ! Surface field - always scalar
      allocate(z0(1))  ! Surface field - always scalar
      allocate(species_mw_g(n_species))
      allocate(species_dd_f0(n_species))
      allocate(species_dd_hstar(n_species))
      allocate(species_dd_dvzaersnow(n_species))
      allocate(species_dd_dvzminval_snow(n_species))
      allocate(species_dd_dvzminval_land(n_species))
      species_tendencies = 0.0_fp

      ! Get meteorological data pointer from virtual column (VirtualMet pattern)
      met => column%get_met()

      ! Now allocate categorical fields using the met pointer dimensions
      allocate(frlai(size(met%FRLAI)))  ! Categorical field - get size from met pointer
      allocate(frlanduse(size(met%FRLANDUSE)))  ! Categorical field - get size from met pointer
      allocate(iland(size(met%ILAND)))  ! Categorical field - get size from met pointer

      ! Extract required fields from met pointer based on field type and processing mode
      bxheight(1:n_levels) = met%BXHEIGHT(1:n_levels)  ! Atmospheric field - always n_levels
      cldfrc(1) = met%CLDFRC  ! Surface field - scalar access
      frlai(:) = met%FRLAI(:)  ! Categorical field - full dimension
      frlanduse(:) = met%FRLANDUSE(:)  ! Categorical field - full dimension
      iland(:) = met%ILAND(:)  ! Categorical field - full dimension
      isice(1) = met%IsIce  ! Surface field - scalar access
      island(1) = met%IsLand  ! Surface field - scalar access
      issnow(1) = met%IsSnow  ! Surface field - scalar access
      lat(1) = met%LAT  ! Surface field - scalar access
      lon(1) = met%LON  ! Surface field - scalar access
      lucname(1) = met%LUCNAME  ! Surface field - scalar access
      obk(1) = met%OBK  ! Surface field - scalar access
      ps(1) = met%PS  ! Surface field - scalar access
      salinity(1) = met%SALINITY  ! Surface field - scalar access
      suncosmid(1) = met%SUNCOSmid  ! Surface field - scalar access
      swgdn(1) = met%SWGDN  ! Surface field - scalar access
      ts(1) = met%TS  ! Surface field - scalar access
      tskin(1) = met%TSKIN  ! Surface field - scalar access
      tstep(1) = this%get_timestep()  ! Special timestep field - retrieved from ProcessInterface
      ustar(1) = met%USTAR  ! Surface field - scalar access
      z0(1) = met%Z0  ! Surface field - scalar access

      ! Get species concentrations from virtual column
      ! Surface-only processing - get surface level concentrations
      do i = 1, n_species
         species_conc(1, i) = column%get_chem_field(species_indices(i), 1)
      end do

      ! Get species properties from configuration (pre-loaded during initialization)
      ! Use species properties from process configuration
      species_mw_g(1:n_species) = this%process_config%drydep_config%species_mw_g(1:n_species)
      ! Use species properties from process configuration
      species_dd_f0(1:n_species) = this%process_config%drydep_config%species_dd_f0(1:n_species)
      ! Use species properties from process configuration
      species_dd_hstar(1:n_species) = this%process_config%drydep_config%species_dd_hstar(1:n_species)
      ! Use species properties from process configuration
      species_dd_dvzaersnow(1:n_species) = this%process_config%drydep_config%species_dd_DvzAerSnow(1:n_species)
      ! Use species properties from process configuration
      species_dd_dvzminval_snow(1:n_species) = this%process_config%drydep_config%species_dd_DvzMinVal_snow(1:n_species)
      ! Use species properties from process configuration
      species_dd_dvzminval_land(1:n_species) = this%process_config%drydep_config%species_dd_DvzMinVal_land(1:n_species)

      ! Call the science scheme with optional diagnostic parameters
      ! Note: wesely uses the following diagnostic fields (if diagnostics enabled):
      ! - drydep_con_per_species (Dry deposition concentration per species)
      ! - drydep_velocity_per_species (Dry deposition velocity)
      if (this%process_config%drydep_config%diagnostics) then
         ! Call with diagnostic outputs enabled
         call compute_wesely( &
            n_levels, &
            n_species, &
            this%process_config%wesely_config, &
            bxheight, &
            cldfrc(1), &
            frlai, &
            frlanduse, &
            iland, &
            isice(1), &
            island(1), &
            issnow(1), &
            lat(1), &
            lon(1), &
            lucname(1), &
            obk(1), &
            ps(1), &
            salinity(1), &
            suncosmid(1), &
            swgdn(1), &
            ts(1), &
            tskin(1), &
            tstep(1), &
            ustar(1), &
            z0(1)            , &
            species_mw_g, &
            species_dd_f0, &
            this%process_config%drydep_config%species_names, &
            species_dd_hstar, &
            species_dd_dvzaersnow, &
            species_dd_dvzminval_snow, &
            species_dd_dvzminval_land, &
            species_conc, &
            species_tendencies, &
            this%process_config%drydep_config%is_gas, &
            this%column_drydep_con_per_species, &
            this%column_drydep_velocity_per_species, &
            this%process_config%drydep_config%diagnostic_species_id         )
      else
         ! Call without diagnostic outputs (optional parameters not passed)
         call compute_wesely( &
            n_levels, &
            n_species, &
            this%process_config%wesely_config, &
            bxheight, &
            cldfrc(1), &
            frlai, &
            frlanduse, &
            iland, &
            isice(1), &
            island(1), &
            issnow(1), &
            lat(1), &
            lon(1), &
            lucname(1), &
            obk(1), &
            ps(1), &
            salinity(1), &
            suncosmid(1), &
            swgdn(1), &
            ts(1), &
            tskin(1), &
            tstep(1), &
            ustar(1), &
            z0(1)            , &
            species_mw_g, &
            species_dd_f0, &
            this%process_config%drydep_config%species_names, &
            species_dd_hstar, &
            species_dd_dvzaersnow, &
            species_dd_dvzminval_snow, &
            species_dd_dvzminval_land, &
            species_conc, &
            species_tendencies, &
            this%process_config%drydep_config%is_gas &
            )
      end if

      ! Apply tendencies back to virtual column based on tendency_mode
      ! Surface-only processing - apply tendencies to surface level only
      do i = 1, n_species
         ! Skip species that don't match scheme type (gas vs aerosol)
         if (.not. this%process_config%drydep_config%is_gas(i)) cycle
         ! Multiplicative tendency: new_conc = conc * (1 - loss_fraction)
         ! where loss_fraction = 1 - exp(-tendency * dt)
         loss_fraction = max(1.0_fp - exp(-species_tendencies(1, i) * this%get_timestep()), 0.0_fp)
         call column%set_chem_field(1, species_indices(i), &
            species_conc(1, i) * (1.0_fp - loss_fraction))
      end do

   end subroutine run_wesely_scheme_column

   subroutine run_gocart_scheme_column(this, column, rc)
      class(ProcessDryDepInterface), intent(inout) :: this
      type(VirtualColumnType), intent(inout) :: column
      integer, intent(out) :: rc

      ! Local variables for scheme calculation
      type(VirtualMetType), pointer :: met => null()  ! Pointer to meteorological data
      ! Meteorological fields
      real(fp), allocatable :: airden(:)
      real(fp), allocatable :: frlake(:)
      real(fp), allocatable :: gwettop(:)
      real(fp), allocatable :: hflux(:)
      integer, allocatable :: lwi(:)
      real(fp), allocatable :: pblh(:)
      real(fp), allocatable :: t(:)
      real(fp), allocatable :: tstep(:)
      real(fp), allocatable :: u10m(:)
      real(fp), allocatable :: ustar(:)
      real(fp), allocatable :: v10m(:)
      real(fp), allocatable :: z(:)
      real(fp), allocatable :: z0h(:)
      ! Species properties
      real(fp), allocatable :: species_density(:)
      real(fp), allocatable :: species_radius(:)
      logical, allocatable :: species_is_seasalt(:)
      real(fp), allocatable :: species_conc(:,:)
      real(fp), allocatable :: species_tendencies(:,:)
      integer :: n_species, n_levels, i
      integer, allocatable :: species_indices(:)
      real(fp) :: loss_fraction  ! Loss fraction for multiplicative tendencies

      rc = cc_success

      ! Get dimensions from virtual column
      n_levels = 1  ! Surface-only processing

      ! Get drydep species information from process configuration
      n_species = this%process_config%drydep_config%n_species
      if (n_species <= 0) then
         return
      end if

      ! Get species indices directly from configuration (pre-computed)
      allocate(species_indices(n_species))
      species_indices(1:n_species) = this%process_config%drydep_config%species_indices(1:n_species)

      ! Allocate arrays
      allocate(species_conc(1, n_species))
      allocate(species_tendencies(1, n_species))
      ! Allocate meteorological field arrays based on field type and process configuration
      allocate(airden(n_levels))  ! Atmospheric field - always n_levels
      allocate(frlake(1))  ! Surface field - always scalar
      allocate(gwettop(1))  ! Surface field - always scalar
      allocate(hflux(1))  ! Surface field - always scalar
      allocate(lwi(1))  ! Surface field - always scalar
      allocate(pblh(1))  ! Surface field - always scalar
      allocate(t(n_levels))  ! Atmospheric field - always n_levels
      allocate(tstep(1))  ! Special timestep field - scalar
      allocate(u10m(1))  ! Surface field - always scalar
      allocate(ustar(1))  ! Surface field - always scalar
      allocate(v10m(1))  ! Surface field - always scalar
      allocate(z(n_levels+1))  ! Edge field - always n_levels+1
      allocate(z0h(1))  ! Surface field - always scalar
      allocate(species_density(n_species))
      allocate(species_radius(n_species))
      allocate(species_is_seasalt(n_species))
      species_tendencies = 0.0_fp

      ! Get meteorological data pointer from virtual column (VirtualMet pattern)
      met => column%get_met()

      ! Now allocate categorical fields using the met pointer dimensions

      ! Extract required fields from met pointer based on field type and processing mode
      airden(1:n_levels) = met%AIRDEN(1:n_levels)  ! Atmospheric field - always n_levels
      frlake(1) = met%FRLAKE  ! Surface field - scalar access
      gwettop(1) = met%GWETTOP  ! Surface field - scalar access
      hflux(1) = met%HFLUX  ! Surface field - scalar access
      lwi(1) = met%LWI  ! Surface field - scalar access
      pblh(1) = met%PBLH  ! Surface field - scalar access
      t(1:n_levels) = met%T(1:n_levels)  ! Atmospheric field - always n_levels
      tstep(1) = this%get_timestep()  ! Special timestep field - retrieved from ProcessInterface
      u10m(1) = met%U10M  ! Surface field - scalar access
      ustar(1) = met%USTAR  ! Surface field - scalar access
      v10m(1) = met%V10M  ! Surface field - scalar access
      z(1:n_levels+1) = met%Z(1:n_levels+1)  ! Edge field - always n_levels+1
      z0h(1) = met%Z0H  ! Surface field - scalar access

      ! Get species concentrations from virtual column
      ! Surface-only processing - get surface level concentrations
      do i = 1, n_species
         species_conc(1, i) = column%get_chem_field(species_indices(i), 1)
      end do

      ! Get species properties from configuration (pre-loaded during initialization)
      ! Use species properties from process configuration
      species_density(1:n_species) = this%process_config%drydep_config%species_density(1:n_species)
      ! Use species properties from process configuration
      species_radius(1:n_species) = this%process_config%drydep_config%species_radius(1:n_species)
      ! Use species properties from process configuration
      species_is_seasalt(1:n_species) = this%process_config%drydep_config%species_is_seasalt(1:n_species)

      ! Call the science scheme with optional diagnostic parameters
      ! Note: gocart uses the following diagnostic fields (if diagnostics enabled):
      ! - drydep_con_per_species (Dry deposition concentration per species)
      ! - drydep_velocity_per_species (Dry deposition velocity)
      if (this%process_config%drydep_config%diagnostics) then
         ! Call with diagnostic outputs enabled
         call compute_gocart( &
            n_levels, &
            n_species, &
            this%process_config%gocart_config, &
            airden, &
            frlake(1), &
            gwettop(1), &
            hflux(1), &
            lwi(1), &
            pblh(1), &
            t, &
            tstep(1), &
            u10m(1), &
            ustar(1), &
            v10m(1), &
            z, &
            z0h(1)            , &
            species_density, &
            species_radius, &
            species_is_seasalt, &
            species_conc, &
            species_tendencies, &
            this%process_config%drydep_config%is_gas, &
            this%column_drydep_con_per_species, &
            this%column_drydep_velocity_per_species, &
            this%process_config%drydep_config%diagnostic_species_id         )
      else
         ! Call without diagnostic outputs (optional parameters not passed)
         call compute_gocart( &
            n_levels, &
            n_species, &
            this%process_config%gocart_config, &
            airden, &
            frlake(1), &
            gwettop(1), &
            hflux(1), &
            lwi(1), &
            pblh(1), &
            t, &
            tstep(1), &
            u10m(1), &
            ustar(1), &
            v10m(1), &
            z, &
            z0h(1)            , &
            species_density, &
            species_radius, &
            species_is_seasalt, &
            species_conc, &
            species_tendencies, &
            this%process_config%drydep_config%is_gas &
            )
      end if

      ! Apply tendencies back to virtual column based on tendency_mode
      ! Surface-only processing - apply tendencies to surface level only
      do i = 1, n_species
         ! Skip species that don't match scheme type (gas vs aerosol)
         if (this%process_config%drydep_config%is_gas(i)) cycle
         ! Multiplicative tendency: new_conc = conc * (1 - loss_fraction)
         ! where loss_fraction = 1 - exp(-tendency * dt)
         loss_fraction = max(1.0_fp - exp(-species_tendencies(1, i) * this%get_timestep()), 0.0_fp)
         call column%set_chem_field(1, species_indices(i), &
            species_conc(1, i) * (1.0_fp - loss_fraction))
      end do

   end subroutine run_gocart_scheme_column

   subroutine run_zhang_scheme_column(this, column, rc)
      class(ProcessDryDepInterface), intent(inout) :: this
      type(VirtualColumnType), intent(inout) :: column
      integer, intent(out) :: rc

      ! Local variables for scheme calculation
      type(VirtualMetType), pointer :: met => null()  ! Pointer to meteorological data
      ! Meteorological fields
      real(fp), allocatable :: bxheight(:)
      real(fp), allocatable :: frlanduse(:)
      integer, allocatable :: iland(:)
      logical, allocatable :: isice(:)
      logical, allocatable :: issnow(:)
      character(len=255), allocatable :: lucname(:)
      real(fp), allocatable :: obk(:)
      real(fp), allocatable :: ps(:)
      real(fp), allocatable :: rh(:)
      real(fp), allocatable :: ts(:)
      real(fp), allocatable :: tstep(:)
      real(fp), allocatable :: u10m(:)
      real(fp), allocatable :: ustar(:)
      real(fp), allocatable :: v10m(:)
      real(fp), allocatable :: z0(:)
      ! Species properties
      real(fp), allocatable :: species_mw_g(:)
      real(fp), allocatable :: species_radius(:)
      real(fp), allocatable :: species_density(:)
      real(fp), allocatable :: species_dd_hstar(:)
      real(fp), allocatable :: species_dd_DvzAerSnow(:)
      real(fp), allocatable :: species_dd_DvzMinVal_snow(:)
      real(fp), allocatable :: species_dd_DvzMinVal_land(:)
      real(fp), allocatable :: species_lower_radius(:)
      real(fp), allocatable :: species_upper_radius(:)
      logical, allocatable :: species_is_dust(:)
      logical, allocatable :: species_is_seasalt(:)
      real(fp), allocatable :: species_conc(:,:)
      real(fp), allocatable :: species_tendencies(:,:)
      integer :: n_species, n_levels, i
      integer, allocatable :: species_indices(:)
      real(fp) :: loss_fraction  ! Loss fraction for multiplicative tendencies

      rc = cc_success

      ! Get dimensions from virtual column
      n_levels = 1  ! Surface-only processing

      ! Get drydep species information from process configuration
      n_species = this%process_config%drydep_config%n_species
      if (n_species <= 0) then
         return
      end if

      ! Get species indices directly from configuration (pre-computed)
      allocate(species_indices(n_species))
      species_indices(1:n_species) = this%process_config%drydep_config%species_indices(1:n_species)

      ! Allocate arrays
      allocate(species_conc(1, n_species))
      allocate(species_tendencies(1, n_species))
      ! Allocate meteorological field arrays based on field type and process configuration
      allocate(bxheight(n_levels))  ! Atmospheric field - always n_levels


      allocate(isice(1))  ! Surface field - always scalar
      allocate(issnow(1))  ! Surface field - always scalar
      allocate(lucname(1))  ! Surface field - always scalar
      allocate(obk(1))  ! Surface field - always scalar
      allocate(ps(1))  ! Surface field - always scalar
      allocate(rh(n_levels))  ! Atmospheric field - always n_levels
      allocate(ts(1))  ! Surface field - always scalar
      allocate(tstep(1))  ! Special timestep field - scalar
      allocate(u10m(1))  ! Surface field - always scalar
      allocate(ustar(1))  ! Surface field - always scalar
      allocate(v10m(1))  ! Surface field - always scalar
      allocate(z0(1))  ! Surface field - always scalar
      allocate(species_mw_g(n_species))
      allocate(species_radius(n_species))
      allocate(species_density(n_species))
      allocate(species_dd_hstar(n_species))
      allocate(species_dd_dvzaersnow(n_species))
      allocate(species_dd_dvzminval_snow(n_species))
      allocate(species_dd_dvzminval_land(n_species))
      allocate(species_lower_radius(n_species))
      allocate(species_upper_radius(n_species))
      allocate(species_is_dust(n_species))
      allocate(species_is_seasalt(n_species))
      species_tendencies = 0.0_fp

      ! Get meteorological data pointer from virtual column (VirtualMet pattern)
      met => column%get_met()

      ! Now allocate categorical fields using the met pointer dimensions
      allocate(frlanduse(size(met%FRLANDUSE)))  ! Categorical field - get size from met pointer
      allocate(iland(size(met%ILAND)))  ! Categorical field - get size from met pointer

      ! Extract required fields from met pointer based on field type and processing mode
      bxheight(1:n_levels) = met%BXHEIGHT(1:n_levels)  ! Atmospheric field - always n_levels
      frlanduse(:) = met%FRLANDUSE(:)  ! Categorical field - full dimension
      iland(:) = met%ILAND(:)  ! Categorical field - full dimension
      isice(1) = met%IsIce  ! Surface field - scalar access
      issnow(1) = met%IsSnow  ! Surface field - scalar access
      lucname(1) = met%LUCNAME  ! Surface field - scalar access
      obk(1) = met%OBK  ! Surface field - scalar access
      ps(1) = met%PS  ! Surface field - scalar access
      rh(1:n_levels) = met%RH(1:n_levels)  ! Atmospheric field - always n_levels
      ts(1) = met%TS  ! Surface field - scalar access
      tstep(1) = this%get_timestep()  ! Special timestep field - retrieved from ProcessInterface
      u10m(1) = met%U10M  ! Surface field - scalar access
      ustar(1) = met%USTAR  ! Surface field - scalar access
      v10m(1) = met%V10M  ! Surface field - scalar access
      z0(1) = met%Z0  ! Surface field - scalar access

      ! Get species concentrations from virtual column
      ! Surface-only processing - get surface level concentrations
      do i = 1, n_species
         species_conc(1, i) = column%get_chem_field(species_indices(i), 1)
      end do

      ! Get species properties from configuration (pre-loaded during initialization)
      ! Use species properties from process configuration
      species_mw_g(1:n_species) = this%process_config%drydep_config%species_mw_g(1:n_species)
      ! Use species properties from process configuration
      species_radius(1:n_species) = this%process_config%drydep_config%species_radius(1:n_species)
      ! Use species properties from process configuration
      species_density(1:n_species) = this%process_config%drydep_config%species_density(1:n_species)
      ! Use species properties from process configuration
      species_dd_hstar(1:n_species) = this%process_config%drydep_config%species_dd_hstar(1:n_species)
      ! Use species properties from process configuration
      species_dd_dvzaersnow(1:n_species) = this%process_config%drydep_config%species_dd_DvzAerSnow(1:n_species)
      ! Use species properties from process configuration
      species_dd_dvzminval_snow(1:n_species) = this%process_config%drydep_config%species_dd_DvzMinVal_snow(1:n_species)
      ! Use species properties from process configuration
      species_dd_dvzminval_land(1:n_species) = this%process_config%drydep_config%species_dd_DvzMinVal_land(1:n_species)
      ! Use species properties from process configuration
      species_lower_radius(1:n_species) = this%process_config%drydep_config%species_lower_radius(1:n_species)
      ! Use species properties from process configuration
      species_upper_radius(1:n_species) = this%process_config%drydep_config%species_upper_radius(1:n_species)
      ! Use species properties from process configuration
      species_is_dust(1:n_species) = this%process_config%drydep_config%species_is_dust(1:n_species)
      ! Use species properties from process configuration
      species_is_seasalt(1:n_species) = this%process_config%drydep_config%species_is_seasalt(1:n_species)

      ! Call the science scheme with optional diagnostic parameters
      ! Note: zhang uses the following diagnostic fields (if diagnostics enabled):
      ! - drydep_con_per_species (Dry deposition concentration per species)
      ! - drydep_velocity_per_species (Dry deposition velocity)
      if (this%process_config%drydep_config%diagnostics) then
         ! Call with diagnostic outputs enabled
         call compute_zhang( &
            n_levels, &
            n_species, &
            this%process_config%zhang_config, &
            bxheight, &
            frlanduse, &
            iland, &
            isice(1), &
            issnow(1), &
            lucname(1), &
            obk(1), &
            ps(1), &
            rh, &
            ts(1), &
            tstep(1), &
            u10m(1), &
            ustar(1), &
            v10m(1), &
            z0(1)            , &
            species_mw_g, &
            species_radius, &
            species_density, &
            this%process_config%drydep_config%species_names, &
            species_dd_hstar, &
            species_dd_dvzaersnow, &
            species_dd_dvzminval_snow, &
            species_dd_dvzminval_land, &
            species_lower_radius, &
            species_upper_radius, &
            species_is_dust, &
            species_is_seasalt, &
            species_conc, &
            species_tendencies, &
            this%process_config%drydep_config%is_gas, &
            this%column_drydep_con_per_species, &
            this%column_drydep_velocity_per_species, &
            this%process_config%drydep_config%diagnostic_species_id         )
      else
         ! Call without diagnostic outputs (optional parameters not passed)
         call compute_zhang( &
            n_levels, &
            n_species, &
            this%process_config%zhang_config, &
            bxheight, &
            frlanduse, &
            iland, &
            isice(1), &
            issnow(1), &
            lucname(1), &
            obk(1), &
            ps(1), &
            rh, &
            ts(1), &
            tstep(1), &
            u10m(1), &
            ustar(1), &
            v10m(1), &
            z0(1)            , &
            species_mw_g, &
            species_radius, &
            species_density, &
            this%process_config%drydep_config%species_names, &
            species_dd_hstar, &
            species_dd_dvzaersnow, &
            species_dd_dvzminval_snow, &
            species_dd_dvzminval_land, &
            species_lower_radius, &
            species_upper_radius, &
            species_is_dust, &
            species_is_seasalt, &
            species_conc, &
            species_tendencies, &
            this%process_config%drydep_config%is_gas &
            )
      end if

      ! Apply tendencies back to virtual column based on tendency_mode
      ! Surface-only processing - apply tendencies to surface level only
      do i = 1, n_species
         ! Skip species that don't match scheme type (gas vs aerosol)
         if (this%process_config%drydep_config%is_gas(i)) cycle
         ! Multiplicative tendency: new_conc = conc * (1 - loss_fraction)
         ! where loss_fraction = 1 - exp(-tendency * dt)
         loss_fraction = max(1.0_fp - exp(-species_tendencies(1, i) * this%get_timestep()), 0.0_fp)
         call column%set_chem_field(1, species_indices(i), &
            species_conc(1, i) * (1.0_fp - loss_fraction))
      end do

   end subroutine run_zhang_scheme_column



   function get_required_met_fields(this) result(field_names)
      class(ProcessDryDepInterface), intent(in) :: this
      character(len=32), allocatable :: field_names(:)
      character(len=32), allocatable :: process_fields(:)
      character(len=32), allocatable :: gas_scheme_fields(:), aero_scheme_fields(:)
      integer :: gas_scheme_count, aero_scheme_count
      character(len=32), allocatable :: unique_fields(:)
      integer :: total_fields, process_count, i, j, unique_count
      logical :: is_duplicate

      ! No process-level required fields
      process_count = 0
      allocate(process_fields(0))

      ! For gas/aero differentiated processes, get fields from both schemes
      ! Get gas scheme fields
      select case (trim(this%process_config%drydep_config%gas_scheme))
       case ('wesely')
         gas_scheme_count = 21
         allocate(gas_scheme_fields(gas_scheme_count))
         gas_scheme_fields(1) = 'USTAR'
         gas_scheme_fields(2) = 'TSTEP'
         gas_scheme_fields(3) = 'TS'
         gas_scheme_fields(4) = 'SWGDN'
         gas_scheme_fields(5) = 'SUNCOSmid'
         gas_scheme_fields(6) = 'OBK'
         gas_scheme_fields(7) = 'CLDFRC'
         gas_scheme_fields(8) = 'BXHEIGHT'
         gas_scheme_fields(9) = 'Z0'
         gas_scheme_fields(10) = 'PS'
         gas_scheme_fields(11) = 'FRLAI'
         gas_scheme_fields(12) = 'ILAND'
         gas_scheme_fields(13) = 'SALINITY'
         gas_scheme_fields(14) = 'FRLANDUSE'
         gas_scheme_fields(15) = 'TSKIN'
         gas_scheme_fields(16) = 'LON'
         gas_scheme_fields(17) = 'LAT'
         gas_scheme_fields(18) = 'LUCNAME'
         gas_scheme_fields(19) = 'IsSnow'
         gas_scheme_fields(20) = 'IsIce'
         gas_scheme_fields(21) = 'IsLand'
       case default
         gas_scheme_count = 0
         allocate(gas_scheme_fields(0))
      end select

      ! Get aerosol scheme fields
      select case (trim(this%process_config%drydep_config%aero_scheme))
       case ('gocart')
         aero_scheme_count = 13
         allocate(aero_scheme_fields(aero_scheme_count))
         aero_scheme_fields(1) = 'USTAR'
         aero_scheme_fields(2) = 'TSTEP'
         aero_scheme_fields(3) = 'T'
         aero_scheme_fields(4) = 'AIRDEN'
         aero_scheme_fields(5) = 'Z'
         aero_scheme_fields(6) = 'LWI'
         aero_scheme_fields(7) = 'PBLH'
         aero_scheme_fields(8) = 'HFLUX'
         aero_scheme_fields(9) = 'Z0H'
         aero_scheme_fields(10) = 'U10M'
         aero_scheme_fields(11) = 'V10M'
         aero_scheme_fields(12) = 'FRLAKE'
         aero_scheme_fields(13) = 'GWETTOP'
       case ('zhang')
         aero_scheme_count = 15
         allocate(aero_scheme_fields(aero_scheme_count))
         aero_scheme_fields(1) = 'USTAR'
         aero_scheme_fields(2) = 'TSTEP'
         aero_scheme_fields(3) = 'TS'
         aero_scheme_fields(4) = 'OBK'
         aero_scheme_fields(5) = 'BXHEIGHT'
         aero_scheme_fields(6) = 'Z0'
         aero_scheme_fields(7) = 'RH'
         aero_scheme_fields(8) = 'PS'
         aero_scheme_fields(9) = 'U10M'
         aero_scheme_fields(10) = 'V10M'
         aero_scheme_fields(11) = 'FRLANDUSE'
         aero_scheme_fields(12) = 'ILAND'
         aero_scheme_fields(13) = 'LUCNAME'
         aero_scheme_fields(14) = 'IsSnow'
         aero_scheme_fields(15) = 'IsIce'
       case default
         aero_scheme_count = 0
         allocate(aero_scheme_fields(0))
      end select

      ! Combine all fields (process + gas_scheme + aero_scheme) and remove duplicates
      ! First estimate maximum possible size (without duplicates)
      total_fields = process_count + gas_scheme_count + aero_scheme_count
      allocate(unique_fields(total_fields))
      unique_count = 0

      ! Add process-level fields first
      do i = 1, process_count
         unique_count = unique_count + 1
         unique_fields(unique_count) = process_fields(i)
      end do

      ! Add gas scheme fields (check for duplicates)
      do i = 1, gas_scheme_count
         is_duplicate = .false.
         do j = 1, unique_count
            if (trim(gas_scheme_fields(i)) == trim(unique_fields(j))) then
               is_duplicate = .true.
               exit
            end if
         end do
         if (.not. is_duplicate) then
            unique_count = unique_count + 1
            unique_fields(unique_count) = gas_scheme_fields(i)
         end if
      end do

      ! Add aerosol scheme fields (check for duplicates)
      do i = 1, aero_scheme_count
         is_duplicate = .false.
         do j = 1, unique_count
            if (trim(aero_scheme_fields(i)) == trim(unique_fields(j))) then
               is_duplicate = .true.
               exit
            end if
         end do
         if (.not. is_duplicate) then
            unique_count = unique_count + 1
            unique_fields(unique_count) = aero_scheme_fields(i)
         end if
      end do

      ! Allocate final result array with exact size
      allocate(field_names(unique_count))
      field_names(1:unique_count) = unique_fields(1:unique_count)

      ! Clean up temporary arrays
      if (allocated(unique_fields)) deallocate(unique_fields)
      if (allocated(process_fields)) deallocate(process_fields)
      if (allocated(gas_scheme_fields)) deallocate(gas_scheme_fields)
      if (allocated(aero_scheme_fields)) deallocate(aero_scheme_fields)

   end function get_required_met_fields

   function get_required_diagnostic_fields(this) result(field_names)
      class(ProcessDryDepInterface), intent(in) :: this
      character(len=64), allocatable :: field_names(:)

      allocate(field_names(2))
      field_names(1) = 'drydep_con_per_species'
      field_names(2) = 'drydep_velocity_per_species'

   end function get_required_diagnostic_fields


   subroutine register_and_allocate_diagnostics(this, container, rc)
      use diagnosticinterface_mod, only: diagnosticregistrytype, diag_real_2d

      class(ProcessDryDepInterface), intent(inout) :: this
      type(StateManagerType), intent(inout) :: container
      integer, intent(out) :: rc

      type(DiagnosticManagerType), pointer :: diag_mgr
      type(DiagnosticRegistryType), pointer :: registry
      type(GridManagerType), pointer :: grid_mgr
      character(len=256) :: field_name  ! For constructing species-specific field names
      integer :: i  ! Loop variable for diagnostic species
      integer :: nx, ny, nz
      integer :: n_species
      integer :: dims_2d(2)
      integer :: dims_3d_species(3)

      rc = cc_success

      ! Only register diagnostics if enabled in config
      if (.not. this%process_config%drydep_config%diagnostics) then
         return
      endif

      ! Get managers
      diag_mgr => container%get_diagnostic_manager()
      grid_mgr => container%get_grid_manager()

      ! Register this process with diagnostic manager (only once per process)
      call diag_mgr%register_process('drydep', rc)
      if (rc /= cc_success) return

      ! Get the process registry for registering individual diagnostics
      call diag_mgr%get_process_registry('drydep', registry, rc)
      if (rc /= cc_success) return

      ! Get grid dimensions
      call grid_mgr%get_shape(nx, ny, nz)
      dims_2d = [nx, ny]

      ! Get species count for 3D diagnostics
      n_species = this%process_config%drydep_config%n_species
      dims_3d_species = [nx, ny, n_species]

      ! Register drydep_con_per_species
      ! Register individual 2D fields for each diagnostic species (species-only diagnostics)
      if (this%process_config%drydep_config%n_diagnostic_species > 0) then
         do i = 1, this%process_config%drydep_config%n_diagnostic_species
            write(field_name, '(A,A,A)') 'drydep_con_', &
               trim(this%process_config%drydep_config%diagnostic_species(i))
            call this%register_diagnostic_field(registry, trim(field_name), &
               'Dry deposition concentration per species', &
               'ug/kg or ppm', diag_real_2d, &
               'drydep', dims_2d, rc=rc)
            if (rc /= cc_success) return
         end do
      end if
      if (rc /= cc_success) return

      ! Register drydep_velocity_per_species
      ! Register individual 2D fields for each diagnostic species (species-only diagnostics)
      if (this%process_config%drydep_config%n_diagnostic_species > 0) then
         do i = 1, this%process_config%drydep_config%n_diagnostic_species
            write(field_name, '(A,A,A)') 'drydep_velocity_', &
               trim(this%process_config%drydep_config%diagnostic_species(i))
            call this%register_diagnostic_field(registry, trim(field_name), &
               'Dry deposition velocity', &
               'm/s', diag_real_2d, &
               'drydep', dims_2d, rc=rc)
            if (rc /= cc_success) return
         end do
      end if
      if (rc /= cc_success) return

      ! Get selected scheme(s)
      ! For gas/aero differentiated processes, register diagnostics from both schemes
      ! Track registered diagnostics to avoid duplicates
      ! Register gas scheme diagnostics
      select case (trim(this%process_config%drydep_config%gas_scheme))
       case ('wesely')
         ! Register wesely-specific diagnostics (gas)
       case default
         ! Unknown gas scheme
      end select

      ! Register aerosol scheme diagnostics (only if not already registered)
      select case (trim(this%process_config%drydep_config%aero_scheme))
       case ('gocart')
         ! Register gocart-specific diagnostics (aerosol)
       case ('zhang')
         ! Register zhang-specific diagnostics (aerosol)
       case default
         ! Unknown aerosol scheme
      end select

      ! Now allocate diagnostic class members after successful registration
      ! First, deallocate if already allocated (for scheme switching)
      if (allocated(this%column_drydep_con_per_species)) deallocate(this%column_drydep_con_per_species)
      if (allocated(this%column_drydep_velocity_per_species)) deallocate(this%column_drydep_velocity_per_species)

      ! Allocate and initialize scheme-specific diagnostic fields based on selected scheme
      ! For gas/aero differentiated process, allocate diagnostics from both gas and aero schemes
      ! Track allocated diagnostics to avoid duplicates

      ! Allocate common diagnostic fields (used by all schemes)
      ! 1D diagnostic: diagnostic species only - allocated based on n_diagnostic_species
      if (this%process_config%drydep_config%n_diagnostic_species > 0) then
         allocate(this%column_drydep_con_per_species(this%process_config%drydep_config%n_diagnostic_species))
      end if
      if (allocated(this%column_drydep_con_per_species)) this%column_drydep_con_per_species = 0.0_fp
      ! 1D diagnostic: diagnostic species only - allocated based on n_diagnostic_species
      if (this%process_config%drydep_config%n_diagnostic_species > 0) then
         allocate(this%column_drydep_velocity_per_species(this%process_config%drydep_config%n_diagnostic_species))
      end if
      if (allocated(this%column_drydep_velocity_per_species)) this%column_drydep_velocity_per_species = 0.0_fp

      ! Gas scheme diagnostics
      select case (trim(this%process_config%drydep_config%gas_scheme))
       case ("wesely")
         ! Gas scheme-specific diagnostics for wesely
       case default
         ! No gas scheme-specific diagnostics for unknown schemes
      end select

      ! Aerosol scheme diagnostics (only allocate if not already allocated by gas scheme)
      select case (trim(this%process_config%drydep_config%aero_scheme))
       case ("gocart")
         ! Aerosol scheme-specific diagnostics for gocart
       case ("zhang")
         ! Aerosol scheme-specific diagnostics for zhang
       case default
         ! No aerosol scheme-specific diagnostics for unknown schemes
      end select

   end subroutine register_and_allocate_diagnostics

   subroutine calculate_and_update_diagnostics(this, column, container, rc)
      class(ProcessDryDepInterface), intent(inout) :: this
      type(VirtualColumnType), intent(in) :: column
      type(StateManagerType), intent(inout) :: container
      integer, intent(out) :: rc

      integer :: i_col, j_col  ! Column grid position
      integer :: i  ! Loop variable for diagnostic species
      character(len=256) :: field_name  ! For constructing species-specific field names

      rc = cc_success

      ! Skip if diagnostics not enabled
      if (.not. this%process_config%drydep_config%diagnostics) return

      ! Get column grid position (x, y indices)
      call column%get_position(i_col, j_col)

      ! Update common diagnostic fields (used by all schemes)
      ! Update individual species diagnostic fields (species-only diagnostics)
      if (this%process_config%drydep_config%n_diagnostic_species > 0) then
         do i = 1, this%process_config%drydep_config%n_diagnostic_species
            write(field_name, '(A,A,A)') 'drydep_con_', &
               trim(this%process_config%drydep_config%diagnostic_species(i))
            call this%update_scalar_diagnostic_column(trim(field_name), &
               this%column_drydep_con_per_species(i), &
               i_col, j_col, container, rc)
            if (rc /= cc_success) return
         end do
      end if
      ! Update individual species diagnostic fields (species-only diagnostics)
      if (this%process_config%drydep_config%n_diagnostic_species > 0) then
         do i = 1, this%process_config%drydep_config%n_diagnostic_species
            write(field_name, '(A,A,A)') 'drydep_velocity_', &
               trim(this%process_config%drydep_config%diagnostic_species(i))
            call this%update_scalar_diagnostic_column(trim(field_name), &
               this%column_drydep_velocity_per_species(i), &
               i_col, j_col, container, rc)
            if (rc /= cc_success) return
         end do
      end if
      ! Update scheme-specific diagnostic fields based on active scheme
      ! For gas/aero differentiated process, update diagnostics from both gas and aero schemes
      ! Track updated diagnostics to avoid duplicates
      ! Update gas scheme diagnostics
      select case (trim(this%process_config%drydep_config%gas_scheme))
       case ("wesely")
         ! Gas scheme-specific diagnostics for wesely
       case default
         ! No gas scheme diagnostics for unknown schemes
      end select

      ! Update aerosol scheme diagnostics (only if not already updated)
      select case (trim(this%process_config%drydep_config%aero_scheme))
       case ("gocart")
         ! Aerosol scheme-specific diagnostics for gocart
       case ("zhang")
         ! Aerosol scheme-specific diagnostics for zhang
       case default
         ! No aerosol scheme diagnostics for unknown schemes
      end select

   end subroutine calculate_and_update_diagnostics


   subroutine set_drydep_scheme(this, scheme_name, gas_scheme)
      class(ProcessDryDepInterface), intent(inout) :: this
      character(len=*), intent(in) :: scheme_name
      logical, intent(in), optional :: gas_scheme

      logical :: is_gas

      ! Default to setting the single scheme for backward compatibility
      is_gas = .true.
      if (present(gas_scheme)) is_gas = gas_scheme

      if (is_gas) then
         this%process_config%drydep_config%gas_scheme = trim(scheme_name)
      else
         this%process_config%drydep_config%aero_scheme = trim(scheme_name)
      end if

   end subroutine set_drydep_scheme

   function get_drydep_scheme(this, gas_scheme) result(scheme_name)
      class(ProcessDryDepInterface), intent(in) :: this
      logical, intent(in), optional :: gas_scheme
      character(len=64) :: scheme_name

      logical :: is_gas

      ! Default to getting the gas scheme for backward compatibility
      is_gas = .true.
      if (present(gas_scheme)) is_gas = gas_scheme

      if (is_gas) then
         scheme_name = trim(this%process_config%drydep_config%gas_scheme)
      else
         scheme_name = trim(this%process_config%drydep_config%aero_scheme)
      end if

   end function get_drydep_scheme

end module processdrydepinterface_mod