change testcases to do SIMPOL. finalize aero_data from tchem

src/tchem_interface.F90

@@ -35,6 +35,29 @@ function TChem_getNumberOfSpecies() bind(c, name="TChem_getNumberOfSpecies")
     use iso_c_binding
     integer(kind=c_int) :: TChem_getNumberOfSpecies
   end function
+  function TChem_getNumberOfAeroSpecies() bind(c, &
+       name="TChem_getNumberOfAeroSpecies")
+    use iso_c_binding
+    integer(kind=c_int) :: TChem_getNumberOfAeroSpecies
+  end function
+  function TChem_getAerosolSpeciesDensity(i_spec) bind(c, &
+       name="TChem_getAerosolSpeciesDensity")
+    use iso_c_binding
+    integer(kind=c_int) :: i_spec
+    real(kind=c_double) :: TChem_getAerosolSpeciesDensity
+  end function
+  function TChem_getAerosolSpeciesMW(i_spec) bind(c, &
+       name="TChem_getAerosolSpeciesMW")
+    use iso_c_binding
+    integer(kind=c_int) :: i_spec
+    real(kind=c_double) :: TChem_getAerosolSpeciesMW
+  end function
+  function TChem_getAerosolSpeciesKappa(i_spec) bind(c, &
+       name="TChem_getAerosolSpeciesKappa")
+    use iso_c_binding
+    integer(kind=c_int) :: i_spec
+    real(kind=c_double) :: TChem_getAerosolSpeciesKappa
+  end function
   function TChem_getLengthOfStateVector() bind(c, &
       name="TChem_getLengthOfStateVector")
     use iso_c_binding
@@ -52,13 +75,26 @@ subroutine TChem_setStateVector(array, i_batch) bind(c, &
     real(kind=c_double) :: array(*)
     integer(c_int), value :: i_batch
   end subroutine
+  subroutine TChem_setNumberConcentrationVector(array, i_batch) bind(c, &
+       name="TChem_setNumberConcentrationVector")
+    use iso_c_binding
+    real(kind=c_double) :: array(*)
+    integer(c_int), value :: i_batch
+  end subroutine
   integer(kind=c_size_t) function TChem_getSpeciesName(index, result, &
        buffer_size) bind(C, name="TChem_getSpeciesName")
      use iso_c_binding
      integer(kind=c_int), intent(in) :: index
      character(kind=c_char), intent(out) :: result(*)
      integer(kind=c_size_t), intent(in), value :: buffer_size
   end function
+  integer(kind=c_size_t) function TChem_getAerosolSpeciesName(index, result, &
+       buffer_size) bind(C, name="TChem_getAerosolSpeciesName")
+     use iso_c_binding
+     integer(kind=c_int), intent(in) :: index
+     character(kind=c_char), intent(out) :: result(*)
+     integer(kind=c_size_t), intent(in), value :: buffer_size
+  end function
   subroutine TChem_doTimestep(del_t) bind(C, name="TChem_doTimestep")
      use iso_c_binding
      real(kind=c_double), intent(in) :: del_t
@@ -116,7 +152,7 @@ subroutine pmc_tchem_initialize(gas_config_filename, aero_config_filename, &
     integer, intent(in) :: n_grid_cells
 
     integer(kind=c_int) :: nSpec, nAeroSpec
-    integer :: n_species
+    integer :: n_aero_spec
     integer :: i
     real(kind=c_double), dimension(:), allocatable :: array 
     character(:), allocatable ::  val
@@ -127,33 +163,38 @@ subroutine pmc_tchem_initialize(gas_config_filename, aero_config_filename, &
          n_grid_cells)
 
     ! Get size that gas_data should be
-    nSpec = TChem_getNumberOfSpecies() 
-    call ensure_string_array_size(gas_data%name, nSpec)
+    n_gas_spec = TChem_getNumberOfSpecies()
+    call ensure_string_array_size(gas_data%name, n_gas_spec)
 
     ! Populate gas_data with gas species from TChem
-    do i = 1,nSpec
-       val = TChem_species_name(i-1) 
+    do i = 1,n_gas_spec
+       val = TChem_species_name(i-1, .true.)
        gas_data%name(i) = trim(val)
     end do   
 
     ! For output and MPI, this needs to be allocated (for now)
-    allocate(gas_data%mosaic_index(gas_data_n_spec(gas_data)))
+    allocate(gas_data%mosaic_index(n_gas_spec))
     gas_data%mosaic_index(:) = 0
 
     ! TODO: Create aero_data based on TChem input.
     ! From TChem we need:
     !   Species names
     !   Species properties - density, kappa, molecular weight
-    ! n_species = 10
-    ! call ensure_string_array_size(aero_data%name, n_species)
-    ! call ensure_integer_array_size(aero_data%mosaic_index, n_species)
-    ! call ensure_real_array_size(aero_data%wavelengths, n_swbands)
-    ! call ensure_real_array_size(aero_data%density, n_species)
-    ! call ensure_integer_array_size(aero_data%num_ions, n_species)
-    ! call ensure_real_array_size(aero_data%molec_weight, n_species)
-    ! call ensure_real_array_size(aero_data%kappa, n_species)
-    !do i = 1,n_species
-    !end do 
+    n_aero_spec = TChem_getNumberOfAeroSpecies()
+    call ensure_string_array_size(aero_data%name, n_aero_spec)
+    call ensure_integer_array_size(aero_data%mosaic_index, n_aero_spec)
+    call ensure_real_array_size(aero_data%wavelengths, n_swbands)
+    call ensure_real_array_size(aero_data%density, n_aero_spec)
+    call ensure_integer_array_size(aero_data%num_ions, n_aero_spec)
+    call ensure_real_array_size(aero_data%molec_weight, n_aero_spec)
+    call ensure_real_array_size(aero_data%kappa, n_aero_spec)
+    do i = 1,n_aero_spec
+       val = TChem_species_name(i-1, .false.)
+       aero_data%name(i) =  trim(val)
+       aero_data%density(i) = TChem_getAerosolSpeciesDensity(i-1)
+       aero_data%molec_weight(i) = TChem_getAerosolSpeciesMW(i-1)
+       aero_data%kappa(i) = TChem_getAerosolSpeciesKappa(i-1)
+    end do
 
   end subroutine pmc_tchem_initialize
 
@@ -186,20 +227,26 @@ subroutine tchem_to_partmc(aero_data, aero_state, gas_data, gas_state, &
     integer(c_int) :: nSpec, stateVecDim
     integer :: i_part
     real(kind=c_double), dimension(:), allocatable :: stateVector 
+    integer :: n_gas_spec, n_aero_spec
+
+    n_gas_spec = gas_data_n_spec(gas_data)
+    n_aero_spec = aero_data_n_spec(aero_data)
 
     ! Get gas array
     stateVecDim = TChem_getLengthOfStateVector()
-    nSpec = TChem_getNumberOfSpecies()
     allocate(stateVector(stateVecDim))
     call TChem_getStateVector(stateVector, 0)
 
     gas_state%mix_rat = 0.0
     ! Convert from ppm to ppb.
-    gas_state%mix_rat = stateVector(4:nSpec+3) * 1000.d0
+    gas_state%mix_rat = stateVector(4:n_gas_spec+3) * 1000.d0
 
-    ! Map aerosols
     do i_part = 1,aero_state_n_part(aero_state)
-
+       do i_spec = 1,n_aero_spec
+          aero_state%apa%particle(i_part)%vol(i_spec) = stateVector( &
+               n_gas_spec + 3 + i_spec + (i_part -1) * n_aero_spec) &
+               / aero_data%density(i_spec)
+       end do
     end do
 
   end subroutine tchem_to_partmc
@@ -221,33 +268,58 @@ subroutine tchem_from_partmc(aero_data, aero_state, gas_data, gas_state, &
     !> Environment state.
     type(env_state_t), intent(in) :: env_state
 
-    real(kind=dp), allocatable :: stateVector(:)
-    integer :: stateVecDim
+    real(kind=dp), allocatable :: stateVector(:), number_concentration(:)
+    integer :: stateVecDim, tchem_n_part
     integer :: i_part
     integer :: i_water
 
+    integer :: n_gas_spec, n_aero_spec
+
+    n_gas_spec = gas_data_n_spec(gas_data)
+    n_aero_spec = aero_data_n_spec(aero_data)
+
     ! Get size of stateVector
     stateVecDim = TChem_getLengthOfStateVector()
     allocate(stateVector(stateVecDim))
+
+    ! Get size of number concentration
+    tchem_n_part = 20 
+    allocate(number_concentration(tchem_n_part))
+
     ! First three elements are density, pressure and temperature
     stateVector(1) = env_state_air_den(env_state)
     stateVector(2) = env_state%pressure
     stateVector(3) = env_state%temp
 
     ! PartMC uses relative humidity and not H2O mixing ratio.
-    ! Equation 1.10 from Seinfeld and Pandis - Second Edition.
     i_water = gas_data_spec_by_name(gas_data, "H2O")
     gas_state%mix_rat(i_water) = env_state_rel_humid_to_mix_rat(env_state)
     ! Add gas species to state vector. Convert from ppb to ppm.
-    stateVector(4:gas_data_n_spec(gas_data)+3) = gas_state%mix_rat / 1000.d0
+    stateVector(4:n_gas_spec + 3) = gas_state%mix_rat / 1000.d0
 
-    ! TODO: Map aerosols
     do i_part = 1,aero_state_n_part(aero_state)
+      do i_spec = 1,n_aero_spec
+         stateVector(n_gas_spec + 3 + i_spec + (i_part - 1) * n_aero_spec) = &
+              aero_particle_species_mass(aero_state%apa%particle(i_part), &
+              i_spec, aero_data)
+      end do
+      number_concentration(i_part) = aero_state_particle_num_conc( &
+            aero_state, aero_state%apa%particle(i_part), aero_data)
+    end do
 
+    ! FIXME: What do we have to do here for this to work well?
+    do i_part = aero_state_n_part(aero_state)+1,tchem_n_part
+       do i_spec = 1,n_aero_spec
+          stateVector(n_gas_spec + 3 + i_spec + (i_part-1) * n_aero_spec) = &
+1d-10
+       end do
+       number_concentration(i_part) = 0.0d0
     end do
 
     call TChem_setStateVector(stateVector, 0)
 
+    call TChem_setNumberConcentrationVector(number_concentration, 0)
+
   end subroutine tchem_from_partmc
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -286,19 +358,26 @@ end subroutine TChem_initialize
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
   !> Get species name from TChem for a given index.
-  function TChem_species_name(i_spec) result(species_name)
+  function TChem_species_name(i_spec, is_gas) result(species_name)
     use iso_c_binding
 
-    ! Species name.
+    !> Index of species.
+    integer(kind=c_int), intent(in) :: i_spec
+    !> Logical for if the species is a gas
+    logical :: is_gas
+    !> Species name.
     character(:), allocatable :: species_name
 
-    integer(kind=c_int), intent(in) :: i_spec
     character(kind=c_char, len=:), allocatable :: cbuf
     integer(kind=c_size_t) :: N
 
     allocate(character(256) :: cbuf)
     N = len(cbuf)
-    N = TChem_getSpeciesName(i_spec, cbuf, N)
+    if (is_gas) then
+       N = TChem_getSpeciesName(i_spec, cbuf, N)
+    else
+       N = TChem_getAerosolSpeciesName(i_spec, cbuf, N)
+    end if
     allocate(character(N) :: species_name)
     species_name = cbuf(:N)
 

test/tchem/aero_init_comp.dat

@@ -1,2 +1,2 @@
 #    mass fractions
-SI                1
+POA                1

test/tchem/config_aero_cb05cl_ae5_with_SIMPOL.yaml

@@ -0,0 +1,147 @@
+NCAR-version:
+- absolute integration tolerance: 0.001
+  description: gas-phase product from isoprene oxidation by OH
+  molecular weight [kg mol-1]: 0.08806
+  diffusion coeff [m2 s-1]: 9.5e-06
+  name: ISOP-P1
+  notes:
+  - using diffusion coefficient from dry deposition
+  - see notes on ISOP-P1_aero species
+  type: CHEM_SPEC
+- absolute integration tolerance: 0.001
+  description: gas-phase product from isoprene oxidation by O3
+  molecular weight [kg mol-1]: 0.09003
+  name: ISOP-P2
+  diffusion coeff [m2 s-1]: 9.5e-06
+  notes:
+  - using diffusion coefficient from dry deposition
+  - see notes on ISOP-P2_aero species
+  type: CHEM_SPEC
+- absolute integration tolerance: 0.001
+  description: gas-phase product from monoterpene oxidation by OH
+  molecular weight [kg mol-1]: 0.17025
+  diffusion coeff [m2 s-1]: 9.5e-06
+  name: TERP-P1
+  notes:
+  - using diffusion coefficient from dry deposition
+  - see notes on TERP-P1_aero species
+  type: CHEM_SPEC
+- absolute integration tolerance: 0.001
+  description: gas-phase product from monoterpene oxidation by O3
+  molecular weight [kg mol-1]: 0.202162
+  diffusion coeff [m2 s-1]: 9.5e-06
+  name: TERP-P2
+  diffusion coeff [m2 s-1]: 9.5e-06
+  notes:
+  - using diffusion coefficient from dry deposition
+  - see notes on TERP-P2_aero species
+  type: CHEM_SPEC
+- absolute integration tolerance: 1.0e-05
+  density [kg m-3]: 1800.0
+  description: lumped hydrophobic particulate matter
+  kappa: 0.0
+  molecular weight [kg mol-1]: 0.41475
+  name: POA
+  note: Using C30H54 for molecular weight. TODO find best surrogate
+  num_ions: 0
+  phase: AEROSOL
+  type: CHEM_SPEC
+- absolute integration tolerance: 1.0e-05
+  density [kg m-3]: 1800.0
+  description: First isoprene SOA species in 2-product scheme
+  kappa: 0.1
+  molecular weight [kg mol-1]: 0.08806
+  name: ISOP-P1_aero
+  diffusion coeff [m2 s-1]: 9.5e-06
+  note:
+  - Using ketopropanoic acid for molecular weight. TODO find best surrogate
+  - TODO update SIMPOL parameters based on MW of new surrogate
+  num_ions: 0
+  phase: AEROSOL
+  type: CHEM_SPEC
+- absolute integration tolerance: 1.0e-05
+  density [kg m-3]: 1800.0
+  description: Second isoprene SOA species in 2-product scheme
+  kappa: 0.1
+  molecular weight [kg mol-1]: 0.09003
+  name: ISOP-P2_aero
+  diffusion coeff [m2 s-1]: 9.5e-06
+  note:
+  - Using oxalic acid for molecular weight. TODO find best surrogate
+  - TODO update SIMPOL parameters based on MW of new surrogate
+  num_ions: 0
+  phase: AEROSOL
+  type: CHEM_SPEC
+- absolute integration tolerance: 1.0e-05
+  density [kg m-3]: 1800.0
+  description: First monoterpene SOA species in 2-product scheme
+  kappa: 0.1
+  diffusion coeff [m2 s-1]: 9.5e-06
+  molecular weight [kg mol-1]: 0.17025
+  name: TERP-P1_aero
+  note:
+  - Using 2-hydroxy-3-isopropyl-6-methyl-cyclohexanone for molecular weight
+  - TODO find best surrogate
+  - TODO update SIMPOL parameters based on MW of new surrogate
+  num_ions: 0
+  phase: AEROSOL
+  type: CHEM_SPEC
+- absolute integration tolerance: 1.0e-05
+  density [kg m-3]: 1800.0
+  description: Second monoterpene SOA species in 2-product scheme
+  diffusion coeff [m2 s-1]: 9.5e-06
+  kappa: 0.1
+  molecular weight [kg mol-1]: 0.202162
+  name: TERP-P2_aero
+  note:
+  - Using 2-methyl-5-carboxy-2,4-hexodienoic acid. TODO find best surrogate
+  - TODO update SIMPOL parameters based on MW of new surrogate
+  num_ions: 0
+  phase: AEROSOL
+  type: CHEM_SPEC
+- maximum computational particles: 20 
+  name: my aero rep 2
+  type: AERO_REP_SINGLE_PARTICLE
+- name: MONARCH mod37
+  reactions:
+  - B: [3810.0, -21.3, 0.0, 0.0]
+    aerosol phase: organic_matter
+    aerosol-phase species: ISOP-P1_aero
+    gas-phase species: ISOP-P1
+    type: SIMPOL_PHASE_TRANSFER
+  - B:
+    - 3810.0
+    - -20.9
+    - 0.0
+    - 0.0
+    aerosol phase: organic_matter
+    aerosol-phase species: ISOP-P2_aero
+    gas-phase species: ISOP-P2
+    type: SIMPOL_PHASE_TRANSFER
+  - B:
+    - 2190.0
+    - -17.5
+    - 0.0
+    - 0.0
+    aerosol phase: organic_matter
+    aerosol-phase species: TERP-P1_aero
+    gas-phase species: TERP-P1
+    type: SIMPOL_PHASE_TRANSFER
+  - B:
+    - 2190.0
+    - -15.3
+    - 0.0
+    - 0.0
+    aerosol phase: organic_matter
+    aerosol-phase species: TERP-P2_aero
+    gas-phase species: TERP-P2
+    type: SIMPOL_PHASE_TRANSFER
+  type: MECHANISM
+notes:
+- 2-product SOA scheme from MONARCH model. Based on Tsigaridis and Kanakidou (2007).
+- Details in Spada et al. (2013) in prep for Geosci. Model Dev.
+- Gas-phase rate constants taken from CB05 reactions with same reactants
+- TODO the CB05 reactions should be updated/removed to avoid competition with SOA
+  scheme
+- Clausius clapyron parameters (C* and -dH/R) converted to SIMPOL.1 paramaters
+