Add TChem interface for aerosols

Dockerfile.tchem

@@ -28,6 +28,7 @@ RUN apt-get update \
 
 RUN git clone https://github.com/PCLAeroParams/TChem-atm.git /tchem_dir/
 WORKDIR /tchem_dir/
+RUN git checkout -b add_aerosol_driver origin/add_aerosol_driver
 RUN git submodule update --init --recursive
 
 RUN cmake -S /tchem_dir/external/Tines/ext/kokkos -B /build/kokkos_build \

src/run_part.F90

@@ -593,18 +593,12 @@ subroutine spec_file_read_run_part(file, run_part_opt, aero_data, &
        call spec_file_read_gas_state(sub_file, gas_data, gas_state_init)
        call spec_file_close(sub_file)
 
-       if (.not. run_part_opt%do_camp_chem) then
-          call spec_file_read_string(file, 'aerosol_data', sub_filename)
-          call spec_file_open(sub_filename, sub_file)
-          call spec_file_read_aero_data(sub_file, aero_data)
-          call spec_file_close(sub_file)
-       ! FIXME: Temporary to run PartMC. Replace with initialization from TChem
-       else if (run_part_opt%do_tchem) then
+       if (.not. (run_part_opt%do_camp_chem .or. run_part_opt%do_tchem)) then
           call spec_file_read_string(file, 'aerosol_data', sub_filename)
           call spec_file_open(sub_filename, sub_file)
           call spec_file_read_aero_data(sub_file, aero_data)
           call spec_file_close(sub_file)
-       else
+       else if (run_part_opt%do_camp_chem) then
 #ifdef PMC_USE_CAMP
           call aero_data_initialize(aero_data, camp_core)
           call aero_state_initialize(aero_state, aero_data, camp_core)

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,31 @@ subroutine TChem_setStateVector(array, i_batch) bind(c, &
     real(kind=c_double) :: array(*)
     integer(c_int), value :: i_batch
   end subroutine
+  function TChem_getNumberConcentrationVectorSize() bind(c, &
+       name="TChem_getNumberConcentrationVectorSize")
+    use iso_c_binding
+    integer(kind=c_int) :: TChem_getNumberConcentrationVectorSize
+  end function
+  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
@@ -115,8 +156,7 @@ subroutine pmc_tchem_initialize(gas_config_filename, aero_config_filename, &
     !> Number of cells to solve.
     integer, intent(in) :: n_grid_cells
 
-    integer(kind=c_int) :: nSpec, nAeroSpec
-    integer :: n_species
+    integer(kind=c_int) :: n_aero_spec, n_gas_spec
     integer :: i
     real(kind=c_double), dimension(:), allocatable :: array 
     character(:), allocatable ::  val
@@ -127,33 +167,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
 
@@ -184,22 +229,28 @@ subroutine tchem_to_partmc(aero_data, aero_state, gas_data, gas_state, &
     type(env_state_t), intent(in) :: env_state
 
     integer(c_int) :: nSpec, stateVecDim
-    integer :: i_part
+    integer :: i_part, i_spec
     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 +272,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, i_spec
     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 = TChem_getNumberConcentrationVectorSize()
+    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 +362,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