camp_solver_data.F90

Go to the documentation of this file.

! Copyright (C) 2021 Barcelona Supercomputing Center and University of
! Illinois at Urbana-Champaign
! SPDX-License-Identifier: MIT
 
!> \file
!> The camp_camp_solver_data module.
 
!> The camp_solver_data_t structure and associated subroutines.
module camp_camp_solver_data
 
  use camp_aero_phase_data
  use camp_aero_rep_data
  use camp_aero_rep_factory
  use camp_constants,                   only : i_kind, dp
  use camp_mechanism_data
  use camp_camp_state
  use camp_rxn_data
  use camp_rxn_factory
  use camp_solver_stats
  use camp_sub_model_data
  use camp_sub_model_factory
  use camp_util,                        only : assert_msg, to_string, &
                                              warn_assert_msg, die_msg
 
  use iso_c_binding
 
  implicit none
  private
 
  public :: camp_solver_data_t
 
  !> Default relative tolerance for integration
  real(kind=dp), parameter :: camp_solver_default_rel_tol = 1.0d-8
  !> Default max number of integration steps
  integer(kind=i_kind), parameter :: camp_solver_default_max_steps = 10000
  !> Default maximum number of integration convergence failures
  integer(kind=i_kind), parameter :: camp_solver_default_max_conv_fails = 1000
 
  !> Result code indicating successful completion
  integer, parameter :: camp_solver_success = 0
 
  !> Interface to c ODE solver functions
  interface
    !> Get a new solver
    type(c_ptr) function solver_new(n_state_var, n_cells, var_type, &
                    n_rxn, n_rxn_int_param, n_rxn_float_param, &
                    n_rxn_env_param, n_aero_phase, n_aero_phase_int_param, &
                    n_aero_phase_float_param, n_aero_rep, &
                    n_aero_rep_int_param, n_aero_rep_float_param, &
                    n_aero_rep_env_param, n_sub_model, n_sub_model_int_param,&
                    n_sub_model_float_param, n_sub_model_env_param) bind (c)
      use iso_c_binding
      !> Number of variables on the state array per grid cell
      !! (including const, PSSA, etc.)
      integer(kind=c_int), value :: n_state_var
      !> Number of cells to compute
      integer(kind=c_int), value :: n_cells
      !> Pointer to array of state variable types (solver, constant, PSSA)
      type(c_ptr), value :: var_type
      !> Number of reactions to solve
      integer(kind=c_int), value :: n_rxn
      !> Total number of integer parameters for all reactions
      integer(kind=c_int), value :: n_rxn_int_param
      !> Total number of floating-point parameters for all reactions
      integer(kind=c_int), value :: n_rxn_float_param
      !> Total number of environment-dependent parameters for all reactions
      integer(kind=c_int), value :: n_rxn_env_param
      !> Number of aerosol phases
      integer(kind=c_int), value :: n_aero_phase
      !> Total number of integer parameters for all aerosol phases
      integer(kind=c_int), value :: n_aero_phase_int_param
      !> Total number of floating-point parameters for all aerosol phases
      integer(kind=c_int), value :: n_aero_phase_float_param
      !> Number of aerosol representations
      integer(kind=c_int), value :: n_aero_rep
      !> Total number of integer parameters for all aerosol representations
      integer(kind=c_int), value :: n_aero_rep_int_param
      !> Total number of floating-point parameters for all aerosol
      !! representations
      integer(kind=c_int), value :: n_aero_rep_float_param
      !> Total number of environment-dependent parameters for all aerosol
      !> representations
      integer(kind=c_int), value :: n_aero_rep_env_param
      !> Number of sub models
      integer(kind=c_int), value :: n_sub_model
      !> Total number of integer parameters for all sub models
      integer(kind=c_int), value :: n_sub_model_int_param
      !> Total number of floating-point parameters for all sub models
      integer(kind=c_int), value :: n_sub_model_float_param
      !> Total number of environment-dependent parameters for all sub models
      integer(kind=c_int), value :: n_sub_model_env_param
    end function solver_new
 
    !> Solver initialization
    subroutine solver_initialize(solver_data, abs_tol, rel_tol, max_steps, &
                    max_conv_fails) bind (c)
      use iso_c_binding
      !> Pointer to a SolverData object
      type(c_ptr), value :: solver_data
      !> Pointer to array of absolute toleracnes
      type(c_ptr), value :: abs_tol
      !> Relative integration tolerance
      real(kind=c_double), value :: rel_tol
      !> Maximum number of internal integration steps
      integer(kind=c_int), value :: max_steps
      !> Maximum number of convergence failures
      integer(kind=c_int), value :: max_conv_fails
    end subroutine solver_initialize
 
#ifdef CAMP_DEBUG
    !> Set the debug output flag for the solver
    integer(kind=c_int) function solver_set_debug_out(solver_data, &
                    do_output) bind (c)
      use iso_c_binding
      !> Pointer to a SolverData object
      type(c_ptr), value :: solver_data
      !> Flag indicating whether to output debugging information
      integer(kind=c_int), value :: do_output
    end function solver_set_debug_out
 
    !> Set the Jacobian evaluation flag for the solver
    integer(kind=c_int) function solver_set_eval_jac(solver_data, &
                    eval_Jac) bind (c)
      use iso_c_binding
      !> Pointer to a SolverData object
      type(c_ptr), value :: solver_data
      !> Flag indicating whether to evaluate the Jacobian during solving
      integer(kind=c_int), value :: eval_Jac
    end function solver_set_eval_jac
#endif
 
    !> Run the solver
    integer(kind=c_int) function solver_run(solver_data, state, env, &
                    t_initial, t_final) bind (c)
      use iso_c_binding
      !> Pointer to the initialized solver data
      type(c_ptr), value :: solver_data
      !> Pointer to the state array
      type(c_ptr), value :: state
      !> Pointer to the environmental state array
      type(c_ptr), value :: env
      !> Initial time (s)
      real(kind=c_double), value :: t_initial
      !> Final time (s)
      real(kind=c_double), value :: t_final
    end function solver_run
 
    !> Reset the solver function timers
    subroutine solver_reset_timers( solver_data ) bind(c)
      use iso_c_binding
      !> Pointer to the solver data
      type(c_ptr), value :: solver_data
    end subroutine solver_reset_timers
 
    !> Get the solver statistics
    subroutine solver_get_statistics( solver_data, solver_flag, num_steps, &
                    RHS_evals, LS_setups, error_test_fails, NLS_iters, &
                    NLS_convergence_fails, DLS_Jac_evals, DLS_RHS_evals, &
                    last_time_step__s, next_time_step__s, Jac_eval_fails, &
                    RHS_evals_total, Jac_evals_total, RHS_time__s, &
                    Jac_time__s, max_loss_precision) bind (c)
      use iso_c_binding
      !> Pointer to the solver data
      type(c_ptr), value :: solver_data
      !> Last flag returned by the solver
      type(c_ptr), value :: solver_flag
      !> Number of steps
      type(c_ptr), value :: num_steps
      !> Right-hand side evaluations
      type(c_ptr), value :: RHS_evals
      !> Linear solver setups
      type(c_ptr), value :: LS_setups
      !> Error test failures
      type(c_ptr), value :: error_test_fails
      !> Non-Linear solver iterations
      type(c_ptr), value :: NLS_iters
      !> Non-Linear solver failures
      type(c_ptr), value :: NLS_convergence_fails
      !> Direct Linear Solver Jacobian evaluations
      type(c_ptr), value :: DLS_Jac_evals
      !> Direct Linear Solver right-hand side evaluations
      type(c_ptr), value :: DLS_RHS_evals
      !> Last time step [s]
      type(c_ptr), value :: last_time_step__s
      !> Next time step [s]
      type(c_ptr), value :: next_time_step__s
      !> Number of Jacobian evaluation failures
      type(c_ptr), value :: Jac_eval_fails
      !> Total number of calls to `f()`
      type(c_ptr), value :: RHS_evals_total
      !> Total number of calls to `Jac()`
      type(c_ptr), value :: Jac_evals_total
      !> Compute time for calls to `f()`
      type(c_ptr), value :: RHS_time__s
      !> Compute time for calls to `Jac()`
      type(c_ptr), value :: Jac_time__s
      !> Maximum loss of precision on last call the f()
      type(c_ptr), value :: max_loss_precision
    end subroutine solver_get_statistics
 
    !> Add condensed reaction data to the solver data block
    subroutine rxn_add_condensed_data(rxn_type, n_int_param, &
                    n_float_param, n_env_param, int_param, float_param, &
                    solver_data)  bind (c)
      use iso_c_binding
      !> Reaction type
      integer(kind=c_int), value :: rxn_type
      !> Number of integer parameters to add
      integer(kind=c_int), value :: n_int_param
      !> Number of floating-point parameters to add
      integer(kind=c_int), value :: n_float_param
      !> Number of environment-dependent parameters to add
      integer(kind=c_int), value :: n_env_param
      !> Pointer to the integer parameter array
      type(c_ptr), value :: int_param
      !> Pointer to the floating-point parameter array
      type(c_ptr), value :: float_param
      !> Pointer to the solver data
      type(c_ptr), value :: solver_data
    end subroutine rxn_add_condensed_data
 
    !> Update reaction data
    subroutine rxn_update_data(cell_id, rxn_id, rxn_type, update_data, &
        solver_data) bind(c)
      use iso_c_binding
      !> Grid cell to update
      integer(kind=c_int), value :: cell_id
      !> Reaction id
      integer(kind=c_int) :: rxn_id
      !> Reaction type to updateto update
      integer(kind=c_int), value :: rxn_type
      !> Data required by reaction for updates
      type(c_ptr), value :: update_data
      !> Solver data
      type(c_ptr), value :: solver_data
    end subroutine rxn_update_data
 
    !> Print the solver data
    subroutine rxn_print_data(solver_data) bind(c)
      use iso_c_binding
      !> Solver data
      type(c_ptr), value :: solver_data
    end subroutine rxn_print_data
 
    !> Add condensed sub model data to the solver data block
    subroutine sub_model_add_condensed_data(sub_model_type, n_int_param, &
                  n_float_param, n_env_param, int_param, float_param, &
                  solver_data) bind(c)
      use iso_c_binding
      !> Sub model  type
      integer(kind=c_int), value :: sub_model_type
      !> Number of integer parameters to add
      integer(kind=c_int), value :: n_int_param
      !> Number of floating-point parameters to add
      integer(kind=c_int), value :: n_float_param
      !> Number of environment-dependent parameters to add
      integer(kind=c_int), value :: n_env_param
      !> Pointer to the integer parameter array
      type(c_ptr), value :: int_param
      !> Pointer to the floating-point parameter array
      type(c_ptr), value :: float_param
      !> Pointer to the solver data
      type(c_ptr), value :: solver_data
    end subroutine sub_model_add_condensed_data
 
    !> Update reaction data
    subroutine sub_model_update_data(cell_id, sub_model_id, sub_model_type, &
        update_data, solver_data) bind(c)
      use iso_c_binding
      !> Grid cell to update
      integer(kind=c_int), value :: cell_id
      !> Sub model solver id
      integer(kind=c_int) :: sub_model_id
      !> Reaction type to updateto update
      integer(kind=c_int), value :: sub_model_type
      !> Data required by reaction for updates
      type(c_ptr), value :: update_data
      !> Solver data
      type(c_ptr), value :: solver_data
    end subroutine sub_model_update_data
 
    !> Print the solver data
    subroutine sub_model_print_data(solver_data) bind(c)
      use iso_c_binding
      !> Solver data
      type(c_ptr), value :: solver_data
    end subroutine sub_model_print_data
 
    !> Add condensed aerosol phase data to the solver data block
    subroutine aero_phase_add_condensed_data(n_int_param, n_float_param, &
                  int_param, float_param, solver_data) bind(c)
      use iso_c_binding
      !> Number of integer parameters to add
      integer(kind=c_int), value :: n_int_param
      !> Number of floating-point parameters to add
      integer(kind=c_int), value :: n_float_param
      !> Pointer to the integer parameter array
      type(c_ptr), value :: int_param
      !> Pointer to the floating-point parameter array
      type(c_ptr), value :: float_param
      !> Pointer to the solver data
      type(c_ptr), value :: solver_data
    end subroutine aero_phase_add_condensed_data
 
    !> Print the solver data
    subroutine aero_phase_print_data(solver_data) bind(c)
      use iso_c_binding
      !> Solver data
      type(c_ptr), value :: solver_data
    end subroutine aero_phase_print_data
 
    !> Add condensed aerosol representation data to the solver data block
    subroutine aero_rep_add_condensed_data(aero_rep_type, n_int_param, &
                  n_float_param, n_env_param, int_param, float_param, &
                  solver_data) bind(c)
      use iso_c_binding
      !> Aerosol representation type
      integer(kind=c_int), value :: aero_rep_type
      !> Number of integer parameters to add
      integer(kind=c_int), value :: n_int_param
      !> Number of floating-point parameters to add
      integer(kind=c_int), value :: n_float_param
      !> Number of environment-dependent parameters to add
      integer(kind=c_int), value :: n_env_param
      !> Pointer to the integer parameter array
      type(c_ptr), value :: int_param
      !> Pointer to the floating-point parameter array
      type(c_ptr), value :: float_param
      !> Pointer to the solver data
      type(c_ptr), value :: solver_data
    end subroutine aero_rep_add_condensed_data
 
    !> Update aerosol representation data
    subroutine aero_rep_update_data(cell_id, aero_rep_id, aero_rep_type, &
        update_data, solver_data)  bind(c)
      use iso_c_binding
      !> Grid cell to update
      integer(kind=c_int), value :: cell_id
      !> Aerosol representation solver id
      integer(kind=c_int) :: aero_rep_id
      !> Aerosol representation type to updateto update
      integer(kind=c_int), value :: aero_rep_type
      !> Data required by aerosol representation for updates
      type(c_ptr), value :: update_data
      !> Solver data
      type(c_ptr), value :: solver_data
    end subroutine aero_rep_update_data
 
    !> Print the aerosol representation data
    subroutine aero_rep_print_data(solver_data) bind(c)
      use iso_c_binding
      !> Pointer to the solver data
      type(c_ptr), value :: solver_data
    end subroutine aero_rep_print_data
 
    !> Free the memory associated with a solver
    pure subroutine solver_free(solver_data) bind(c)
      use iso_c_binding
      !> Pointer to the solver data
      type(c_ptr), value, intent(in) :: solver_data
    end subroutine solver_free
 
  end interface
 
  !> Solver data
  !!
  !! Acts as the interface between the camp-chem module and the solver.
  !! Instances of the type hold a pointer to a solver c object and provide
  !! functions to initialize the solver with model data and run the solver
  !! over a specified time step.
  type :: camp_solver_data_t
    private
    !> C Solver object
    type(c_ptr), public :: solver_c_ptr
    !> Relative tolerance for the integration
    real(kind=dp), public :: rel_tol = camp_solver_default_rel_tol
    !> Maximum number of timesteps
    integer(kind=i_kind), public :: max_steps = camp_solver_default_max_steps
    !> Maximum number of convergence failures
    integer(kind=i_kind), public :: max_conv_fails = &
            camp_solver_default_max_conv_fails
    !> Flag indicating whether the solver was intialized
    logical :: initialized = .false.
  contains
    !> Initialize the solver
    procedure :: initialize
    !> Update sub-model data
    procedure :: update_sub_model_data
    !> Update reactions data
    procedure :: update_rxn_data
    !> Update aerosol representation data
    procedure :: update_aero_rep_data
    !> Integrate over a given time step
    procedure :: solve
    !> Reset the solver function timers
    procedure, private :: reset_timers
    !> Get the solver statistics from the last run
    procedure, private :: get_solver_stats
    !> Checks whether a solver is available
    procedure :: is_solver_available
    !> Print the solver data
    procedure :: print => do_print
    !> Finalize the solver data
    final :: finalize
  end type camp_solver_data_t
 
  ! Constructor for camp_solver_data_t
  interface camp_solver_data_t
    procedure :: constructor
  end interface camp_solver_data_t
 
contains
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  !> Constructor for camp_solver_data_t
  function constructor() result (new_obj)
 
    !> New solver variable
    type(camp_solver_data_t), pointer :: new_obj
 
    !> Allocate space for the new object
    allocate(new_obj)
 
  end function constructor
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  !> Initialize the solver
  subroutine initialize(this, var_type, abs_tol, mechanisms, aero_phases, &
                  aero_reps, sub_models, rxn_phase, n_cells)
 
    !> Solver data
    class(camp_solver_data_t), intent(inout) :: this
    !> Variable type for each species in the state array. This array must be
    !! of the same length as the state array.
    integer(kind=i_kind), allocatable, intent(in) :: var_type(:)
    !> Absolute tolerance for each species in the state array. This array must
    !! be of the same length as the state array. Values for CONST and PSSA
    !! species will be ignored by the solver.
    real(kind=dp), allocatable, intent(in) :: abs_tol(:)
    !> Mechanisms to include in solver
    type(mechanism_data_ptr), pointer, intent(in) :: mechanisms(:)
    !> Aerosol phases to include
    type(aero_phase_data_ptr), pointer, intent(in) :: aero_phases(:)
    !> Aerosol representations to include
    type(aero_rep_data_ptr), pointer, intent(in) :: aero_reps(:)
    !> Sub models to include
    type(sub_model_data_ptr), pointer, intent(in) :: sub_models(:)
    !> Reactions phase to solve -- gas, aerosol, or both (default)
    !! Use parameters in camp_rxn_data to specify phase:
    !! GAS_RXN, AERO_RXN, GAS_AERO_RXN
    integer(kind=i_kind), intent(in) :: rxn_phase
    !> Number of cells to compute
    integer(kind=i_kind), optional :: n_cells
 
    ! Variable types
    integer(kind=c_int), pointer :: var_type_c(:)
    ! Absolute tolerances
    real(kind=c_double), pointer :: abs_tol_c(:)
    ! Indices for iteration
    integer(kind=i_kind) :: i_mech, i_rxn, i_aero_phase, i_aero_rep, &
            i_sub_model
    ! Reaction pointer
    class(rxn_data_t), pointer :: rxn
    ! Reaction factory object for getting reaction type
    type(rxn_factory_t) :: rxn_factory
    ! Aerosol phase pointer
    type(aero_phase_data_t), pointer :: aero_phase
    ! Aerosol representation pointer
    class(aero_rep_data_t), pointer :: aero_rep
    ! Aerosol representation factory object for getting aerosol
    ! representation type
    type(aero_rep_factory_t) :: aero_rep_factory
    ! Sub model pointer
    class(sub_model_data_t), pointer :: sub_model
    ! Sub model factory for getting sub model type
    type(sub_model_factory_t) :: sub_model_factory
    ! Integer parameters being transfered
    integer(kind=c_int), pointer :: int_param(:)
    ! Floating point parameters being transfered
    real(kind=c_double), pointer :: float_param(:)
    ! Number of reactions
    integer(kind=c_int) :: n_rxn
    ! Number of integer reaction parameters
    integer(kind=c_int) :: n_rxn_int_param
    ! Number of floating-point reaction parameters
    integer(kind=c_int) :: n_rxn_float_param
    ! Number of environment-dependent reaction parameters
    integer(kind=c_int) :: n_rxn_env_param
    ! Number of aerosol phases
    integer(kind=c_int) :: n_aero_phase
    ! Number of integer aerosol phase  parameters
    integer(kind=c_int) :: n_aero_phase_int_param
    ! Number of floating-point aerosol phase parameters
    integer(kind=c_int) :: n_aero_phase_float_param
    ! Number of aerosol representations
    integer(kind=c_int) :: n_aero_rep
    ! Number of integer aerosol representation parameters
    integer(kind=c_int) :: n_aero_rep_int_param
    ! Number of floating-point aerosol representation parameters
    integer(kind=c_int) :: n_aero_rep_float_param
    ! Number of environment-dependent aerosol representation parameters
    integer(kind=c_int) :: n_aero_rep_env_param
    ! Number of sub models
    integer(kind=c_int) :: n_sub_model
    ! Number of integer sub model parameters
    integer(kind=c_int) :: n_sub_model_int_param
    ! Number of floating-point sub model parameters
    integer(kind=c_int) :: n_sub_model_float_param
    ! Number of environment-dependent sub model parameters
    integer(kind=c_int) :: n_sub_model_env_param
    ! Number of cells to compute
    integer(kind=c_int) :: l_n_cells
 
    if (present(n_cells)) then
      l_n_cells = n_cells
    else
      l_n_cells = 1
    end if
 
    ! Make sure the variable type and absolute tolerance arrays are of
    ! equal length
    call assert_msg(825843466, size(abs_tol).eq.size(var_type), &
            "Mismatched absolute tolerance and variable type arrays: "// &
            "abs_tol size: "//trim(to_string(size(abs_tol)))// &
            "; var_type: "//trim(to_string(size(var_type))))
 
    ! Set the absolute tolerance and variable type arrays
    allocate(var_type_c(size(var_type)))
    allocate(abs_tol_c(size(abs_tol)))
    var_type_c(:) = int(var_type(:), kind=c_int)
    abs_tol_c(:) = real(abs_tol(:), kind=c_double)
 
    ! Initialize the counters
    n_rxn = 0
    n_rxn_int_param = 0
    n_rxn_float_param = 0
    n_rxn_env_param = 0
 
    ! Calculate the number of reactions and the size of the condensed data
    do i_mech=1, size(mechanisms)
      do i_rxn=1, mechanisms(i_mech)%val%size()
        rxn => mechanisms(i_mech)%val%get_rxn(i_rxn)
        select case (rxn%rxn_phase)
          case (gas_rxn)
            if (rxn_phase.eq.aero_rxn) cycle
          case (aero_rxn)
            if (rxn_phase.eq.gas_rxn) cycle
          case (gas_aero_rxn)
            if (rxn_phase.eq.gas_rxn) cycle
        end select
        n_rxn = n_rxn + 1
        n_rxn_int_param = n_rxn_int_param + size(rxn%condensed_data_int)
        n_rxn_float_param = n_rxn_float_param + size(rxn%condensed_data_real)
        n_rxn_env_param = n_rxn_env_param + rxn%num_env_params
      end do
    end do
    rxn => null()
 
    ! Initialize the aerosol phase counters
    n_aero_phase = size(aero_phases)
    n_aero_phase_int_param = 0
    n_aero_phase_float_param = 0
 
    ! Calculate the size of the aerosol phases condensed data
    do i_aero_phase=1, n_aero_phase
      aero_phase => aero_phases(i_aero_phase)%val
      n_aero_phase_int_param = n_aero_phase_int_param + &
              size(aero_phase%condensed_data_int)
      n_aero_phase_float_param = n_aero_phase_float_param + &
              size(aero_phase%condensed_data_real)
    end do
    aero_phase => null()
 
    ! Initialize the aerosol representation counters
    n_aero_rep = size(aero_reps)
    n_aero_rep_int_param = 0
    n_aero_rep_float_param = 0
    n_aero_rep_env_param = 0
 
    ! Calculate the size of the aerosol representations condensed data
    do i_aero_rep=1, n_aero_rep
      aero_rep => aero_reps(i_aero_rep)%val
      n_aero_rep_int_param = n_aero_rep_int_param + &
              size(aero_rep%condensed_data_int)
      n_aero_rep_float_param = n_aero_rep_float_param + &
              size(aero_rep%condensed_data_real)
      n_aero_rep_env_param = n_aero_rep_env_param + &
              aero_rep%num_env_params
    end do
    aero_rep => null()
 
    ! Initialize the sub model counters
    n_sub_model = size(sub_models)
    n_sub_model_int_param = 0
    n_sub_model_float_param = 0
    n_sub_model_env_param = 0
 
    ! Calculate the size of the sub model condensed data
    do i_sub_model=1, n_sub_model
      sub_model => sub_models(i_sub_model)%val
      n_sub_model_int_param = n_sub_model_int_param + &
              size(sub_model%condensed_data_int)
      n_sub_model_float_param = n_sub_model_float_param + &
              size(sub_model%condensed_data_real)
      n_sub_model_env_param = n_sub_model_env_param + sub_model%num_env_params
    end do
    sub_model => null()
 
    ! Get a new solver object
    this%solver_c_ptr = solver_new( &
            int(size(var_type_c), kind=c_int), & ! Size of the state variable
            l_n_cells,                         & ! # of cells computed at once
            c_loc(var_type_c),                 & ! Variable types
            n_rxn,                             & ! # of reactions
            n_rxn_int_param,                   & ! # of rxn data int params
            n_rxn_float_param,                 & ! # of rxn data real params
            n_rxn_env_param,                   & ! # of rxn env-dependent params
            n_aero_phase,                      & ! # of aero phases
            n_aero_phase_int_param,            & ! # of aero phase int params
            n_aero_phase_float_param,          & ! # of aero phase real params
            n_aero_rep,                        & ! # of aero reps
            n_aero_rep_int_param,              & ! # of aero rep int params
            n_aero_rep_float_param,            & ! # of aero rep real params
            n_aero_rep_env_param,              & ! # of aero rep env params
            n_sub_model,                       & ! # of sub models
            n_sub_model_int_param,             & ! # of sub model int params
            n_sub_model_float_param,           & ! # of sub model real params
            n_sub_model_env_param              & ! # of sub model env params
            )
 
    ! Add all the condensed reaction data to the solver data block for
    ! reactions of the specified phase
    do i_mech=1, size(mechanisms)
      do i_rxn=1, mechanisms(i_mech)%val%size()
 
        ! FIXME Put ZSR aerosol water first, so water is available for other
        ! reactions - and find a better way to account for inter-dependence
        ! of reactions/sub-models
 
        ! Assign rxn to the current reaction
        rxn => mechanisms(i_mech)%val%get_rxn(i_rxn)
 
        ! Check reaction phase
        select case (rxn%rxn_phase)
          case (gas_rxn)
            if (rxn_phase.eq.aero_rxn) cycle
          case (aero_rxn)
            if (rxn_phase.eq.gas_rxn) cycle
          case (gas_aero_rxn)
            if (rxn_phase.eq.gas_rxn) cycle
        end select
 
        ! Load temporary data arrays
        allocate(int_param(size(rxn%condensed_data_int)))
        allocate(float_param(size(rxn%condensed_data_real)))
        int_param(:) = int(rxn%condensed_data_int(:), kind=c_int)
        float_param(:) = real(rxn%condensed_data_real(:), kind=c_double)
 
        ! Send the condensed data to the solver
        call rxn_add_condensed_data ( &
                int(rxn_factory%get_type(rxn), kind=c_int),& ! Rxn type
                int(size(int_param), kind=c_int),          & ! Int array size
                int(size(float_param), kind=c_int),        & ! Real array size
                rxn%num_env_params,                        & ! Env-dep array size
                c_loc(int_param),                          & ! Int array ptr
                c_loc(float_param),                        & ! Real array ptr
                this%solver_c_ptr                          & ! Solver data ptr
                )
 
        ! Deallocate temporary arrays
        deallocate(int_param)
        deallocate(float_param)
 
      end do
    end do
    rxn => null()
 
    ! Add all the condensed aerosol phase data to the solver data block
    do i_aero_phase=1, size(aero_phases)
 
      ! Assign aero_phase to the current aerosol phase
      aero_phase => aero_phases(i_aero_phase)%val
 
      ! Load temporary data arrays
      allocate(int_param(size(aero_phase%condensed_data_int)))
      allocate(float_param(size(aero_phase%condensed_data_real)))
      int_param(:) = int(aero_phase%condensed_data_int(:), kind=c_int)
      float_param(:) = real(aero_phase%condensed_data_real(:), kind=c_double)
 
      ! Send the condensed data to the solver
      call aero_phase_add_condensed_data ( &
              int(size(int_param), kind=c_int),   & ! Int array size
              int(size(float_param), kind=c_int), & ! Real array size
              c_loc(int_param),                   & ! Int array ptr
              c_loc(float_param),                 & ! Real array ptr
              this%solver_c_ptr                   & ! Solver data ptr
              )
 
      ! Deallocate temporary arrays
      deallocate(int_param)
      deallocate(float_param)
 
    end do
    aero_phase => null()
 
    ! Add all the condensed aerosol representation data to the solver data
    ! block
    do i_aero_rep=1, size(aero_reps)
 
      ! Assign aero_rep to the current aerosol representation
      aero_rep => aero_reps(i_aero_rep)%val
 
      ! Load temporary data arrays
      allocate(int_param(size(aero_rep%condensed_data_int)))
      allocate(float_param(size(aero_rep%condensed_data_real)))
      int_param(:) = int(aero_rep%condensed_data_int(:), kind=c_int)
      float_param(:) = real(aero_rep%condensed_data_real(:), kind=c_double)
 
      ! Send the condensed data to the solver
      call aero_rep_add_condensed_data ( &
              int(aero_rep_factory%get_type(aero_rep), kind=c_int), &
                                                    ! Aero rep type
              int(size(int_param), kind=c_int),   & ! Int array size
              int(size(float_param), kind=c_int), & ! Real array size
              aero_rep%num_env_params,            & ! Env-dep array size
              c_loc(int_param),                   & ! Int array ptr
              c_loc(float_param),                 & ! Real array ptr
              this%solver_c_ptr                   & ! Solver data ptr
              )
 
      ! Deallocate temporary arrays
      deallocate(int_param)
      deallocate(float_param)
 
    end do
    aero_rep => null()
 
    ! Add all the condensed sub model data to the solver data block
    do i_sub_model=1, size(sub_models)
 
      ! Assign sub_model to the current sub model
      sub_model => sub_models(i_sub_model)%val
 
      ! Load temporary data arrays
      allocate(int_param(size(sub_model%condensed_data_int)))
      allocate(float_param(size(sub_model%condensed_data_real)))
      int_param(:) = int(sub_model%condensed_data_int(:), kind=c_int)
      float_param(:) = real(sub_model%condensed_data_real(:), kind=c_double)
 
      ! Send the condensed data to the solver
      call sub_model_add_condensed_data ( &
              int(sub_model_factory%get_type(sub_model), kind=c_int), &
                                                    ! Sub model type
              int(size(int_param), kind=c_int),   & ! Int array size
              int(size(float_param), kind=c_int), & ! Real array size
              sub_model%num_env_params,           & ! Env-dep array size
              c_loc(int_param),                   & ! Int array ptr
              c_loc(float_param),                 & ! Real array ptr
              this%solver_c_ptr                   & ! Solver data ptr
              )
 
      ! Deallocate temporary arrays
      deallocate(int_param)
      deallocate(float_param)
 
    end do
    sub_model => null()
 
    ! Initialize the solver
    call solver_initialize( &
            this%solver_c_ptr,                  & ! Pointer to solver data
            c_loc(abs_tol_c),                   & ! Absolute tolerances
            real(this%rel_tol, kind=c_double),  & ! Relative tolerance
            int(this%max_steps, kind=c_int),    & ! Max # of integration steps
            int(this%max_conv_fails, kind=c_int)& ! Max # of convergence fails
            )
 
    ! Flag the solver as initialized
    this%initialized = .true.
 
    ! Free memory allocated for solver initialization
    deallocate(abs_tol_c)
    deallocate(var_type_c)
 
  end subroutine initialize
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  !> Update sub-model data
  subroutine update_sub_model_data(this, update_data)
 
    !> Solver data
    class(camp_solver_data_t), intent(inout) :: this
    !> Update data
    class(sub_model_update_data_t), intent(in) :: update_data
 
    call sub_model_update_data( &
            update_data%get_cell_id()-1,     & ! Grid cell to update
            update_data%sub_model_solver_id, & ! Solver's sub model id
            update_data%get_type(),          & ! Sub-model type to update
            update_data%get_data(),          & ! Data needed to perform update
            this%solver_c_ptr                & ! Pointer to solver data
            )
 
  end subroutine update_sub_model_data
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  !> Update reaction data
  subroutine update_rxn_data(this, update_data)
 
    !> Solver data
    class(camp_solver_data_t), intent(inout) :: this
    !> Update data
    class(rxn_update_data_t), intent(in) :: update_data
 
    call rxn_update_data( &
            update_data%get_cell_id()-1,     & ! Grid cell to update
            update_data%rxn_solver_id,       & ! Solver's reaction id
            update_data%get_type(),          & ! Reaction type to update
            update_data%get_data(),          & ! Data needed to perform update
            this%solver_c_ptr                & ! Pointer to solver data
            )
 
  end subroutine update_rxn_data
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  !> Update aerosol representation data based on data passed from the host
  !! model related to aerosol properties
  subroutine update_aero_rep_data(this, update_data)
 
    !> Solver data
    class(camp_solver_data_t), intent(inout) :: this
    !> Update data
    class(aero_rep_update_data_t), intent(in) :: update_data
 
    call aero_rep_update_data( &
            update_data%get_cell_id()-1,     & ! Grid cell to update
            update_data%aero_rep_solver_id,  & ! Solver's aero rep id
            update_data%get_type(),          & ! Aerosol representation type
            update_data%get_data(),          & ! Data needed to perform update
            this%solver_c_ptr                & ! Pointer to solver data
            )
 
  end subroutine update_aero_rep_data
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  !> Solve the mechanism(s) for a specified timestep
  subroutine solve(this, camp_state, t_initial, t_final, solver_stats)
 
    !> Solver data
    class(camp_solver_data_t), intent(inout) :: this
    !> Model state
    type(camp_state_t), target, intent(inout) :: camp_state
    !> Start time (s)
    real(kind=dp), intent(in) :: t_initial
    !> End time (s)
    real(kind=dp), intent(in) :: t_final
    !> Solver statistics
    type(solver_stats_t), intent(inout), optional, target :: solver_stats
 
    integer(kind=c_int) :: solver_status
 
#ifdef CAMP_DEBUG
    if (present(solver_stats)) then
      ! Update the debugging output flag in the solver data
      if (solver_stats%debug_out) then
        solver_status = solver_set_debug_out( &
            this%solver_c_ptr,              & ! Pointer to the solver data
            int(1, kind=c_int)              & ! Debug flag
            )
      else
        solver_status = solver_set_debug_out( &
            this%solver_c_ptr,              & ! Pointer to the solver data
            int(0, kind=c_int)              & ! Debug flag
            )
      end if
      ! Update Jacobian evaluation flag in the solver data
      if (solver_stats%eval_Jac) then
        solver_stats = solver_set_eval_jac( &
            this%solver_c_ptr,              & ! Pointer to the solver data
            int(1, kind=c_int)              & ! Jac eval flag
            )
      else
        solver_stats = solver_set_eval_jac( &
            this%solver_c_ptr,              & ! Pointer to the solver data
            int(0, kind=c_int)              & ! Jac eval flag
            )
      end if
    end if
 
    ! Reset the solver function timers
    call this%reset_timers( )
#endif
 
    ! Run the solver
    solver_status = solver_run( &
            this%solver_c_ptr,              & ! Pointer to intialized solver
            c_loc(camp_state%state_var),    & ! Pointer to state array
            c_loc(camp_state%env_var),      & ! Pointer to environmental vars
            real(t_initial, kind=c_double), & ! Start time (s)
            real(t_final, kind=c_double)    & ! Final time (s)
            )
 
    ! Get the solver statistics
    if (present(solver_stats)) then
      call this%get_solver_stats( solver_stats )
      solver_stats%status_code   = solver_status
      solver_stats%start_time__s = t_initial
      solver_stats%end_time__s   = t_final
    else
      call warn_assert_msg(997420005, solver_status.eq.0, "Solver failed")
    end if
 
  end subroutine solve
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  !> Reset the solver function timers
  subroutine reset_timers( this )
 
    !> Solver data
    class(camp_solver_data_t), intent(inout) :: this
 
    call solver_reset_timers( this%solver_c_ptr )
 
  end subroutine reset_timers
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  !> Get solver statistics
  subroutine get_solver_stats( this, solver_stats )
 
    !> Solver data
    class(camp_solver_data_t), intent(inout) :: this
    !> Solver statistics
    type(solver_stats_t), intent(inout), target :: solver_stats
 
    call solver_get_statistics( &
            this%solver_c_ptr,                             & ! Solver data
            c_loc( solver_stats%solver_flag           ),   & ! Last flag returned CVode
            c_loc( solver_stats%num_steps             ),   & ! Number of steps
            c_loc( solver_stats%RHS_evals             ),   & ! Right-hand side evals
            c_loc( solver_stats%LS_setups             ),   & ! Linear solver setups
            c_loc( solver_stats%error_test_fails      ),   & ! Error test failures
            c_loc( solver_stats%NLS_iters             ),   & ! Non-Linear solver interations
            c_loc( solver_stats%NLS_convergence_fails ),   & ! Non-Linear solver fails
            c_loc( solver_stats%DLS_Jac_evals         ),   & ! Jacobian evals
            c_loc( solver_stats%DLS_RHS_evals         ),   & ! DLS Right-hand side evals
            c_loc( solver_stats%last_time_step__s     ),   & ! Last time step [s]
            c_loc( solver_stats%next_time_step__s     ),   & ! Next time step [s]
            c_loc( solver_stats%Jac_eval_fails        ),   & ! Number of Jac eval fails
            c_loc( solver_stats%RHS_evals_total       ),   & ! total f() calls
            c_loc( solver_stats%Jac_evals_total       ),   & ! total Jac() calls
            c_loc( solver_stats%RHS_time__s           ),   & ! Compute time f() [s]
            c_loc( solver_stats%Jac_time__s           ),   & ! Compute time Jac() [s]
            c_loc( solver_stats%max_loss_precision    ) )    ! Maximum loss of precision
 
  end subroutine get_solver_stats
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  !> Check whether a solver is available for the integration
  logical function is_solver_available(this)
 
    !> Solver data
    class(camp_solver_data_t), intent(in) :: this
 
#ifdef CAMP_USE_SUNDIALS
    is_solver_available = .true.
#else
    is_solver_available = .false.
#endif
 
  end function is_solver_available
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  !> Print the solver data
  subroutine do_print(this)
 
    !> Solver data
    class(camp_solver_data_t), intent(in) :: this
 
    call rxn_print_data(this%solver_c_ptr)
    call aero_phase_print_data(this%solver_c_ptr)
    call aero_rep_print_data(this%solver_c_ptr)
    call sub_model_print_data(this%solver_c_ptr)
 
  end subroutine do_print
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  !> Finalize the solver data
  elemental subroutine finalize(this)
 
    !> Solver data
    type(camp_solver_data_t), intent(inout) :: this
 
    if (this%initialized) call solver_free(this%solver_c_ptr)
 
  end subroutine finalize
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
end module camp_camp_solver_data