update set and get statevector tchem. add simple test

CMakeLists.txt

@@ -192,6 +192,9 @@ if(ENABLE_MPI)
 endif()
 add_test(test_rand ${CMAKE_BINARY_DIR}/test_run/run_test_directory.sh rand)
 add_test(test_sedi ${CMAKE_BINARY_DIR}/test_run/run_test_directory.sh sedi)
+if (ENABLE_TCHEM)
+  add_test(test_tchem ${CMAKE_BINARY_DIR}/test_run/run_test_directory.sh tchem)
+endif()
 
 ######################################################################
 # partmc library

src/gas_state.F90

@@ -188,6 +188,20 @@ subroutine gas_state_mole_dens_to_ppb(gas_state, env_state)
 
   end subroutine gas_state_mole_dens_to_ppb
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  !> Convert (ppb) to (mol m^{-3}).
+  subroutine gas_state_ppb_to_mole_dens(gas_state, env_state)
+
+    !> Gas state.
+    type(gas_state_t), intent(inout) :: gas_state
+    !> Environment state.
+    type(env_state_t), intent(in) :: env_state
+
+    call gas_state_scale(gas_state, env_state_air_molar_den(env_state) / 1e9)
+
+  end subroutine gas_state_ppb_to_mole_dens
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 #ifdef PMC_USE_CAMP

src/partmc.F90

@@ -511,6 +511,11 @@ subroutine partmc_part(file)
 #endif
        end if
 
+       if (run_part_opt%do_tchem) then
+          call pmc_tchem_initialize(tchem_config_filename, gas_data, &
+               aero_data)
+       end if
+
        if (do_restart) then
           call input_state(restart_filename, dummy_index, dummy_time, &
                dummy_del_t, dummy_i_repeat, run_part_opt%uuid, aero_data, &
@@ -524,11 +529,7 @@ subroutine partmc_part(file)
             call gas_data_initialize(gas_data, camp_core)
 #endif
           else if (run_part_opt%do_tchem) then
-            ! FIXME: Replacde with something else
-            call spec_file_read_string(file, 'gas_data', sub_filename)
-            call spec_file_open(sub_filename, sub_file)
-            call spec_file_read_gas_data(sub_file, gas_data)
-            call spec_file_close(sub_file)
+            print*, 'do nothing!'
           else
             call spec_file_read_string(file, 'gas_data', sub_filename)
             call spec_file_open(sub_filename, sub_file)
@@ -752,10 +753,6 @@ subroutine partmc_part(file)
 #endif
     end if
 
-    if (run_part_opt%do_tchem) then
-       call pmc_tchem_initialize(tchem_config_filename, gas_data, gas_state_init, aero_data)
-    end if
-
     ! re-initialize RNG with the given seed
     call pmc_rand_finalize()
     call pmc_srand(rand_init, pmc_mpi_rank())

src/tchem_interface.F90

@@ -11,15 +11,49 @@ module pmc_tchem_interface
   use pmc_aero_data
   use pmc_aero_particle
   use pmc_aero_state
-  use pmc_constants,                  only : dp
+  use pmc_constants
   use pmc_gas_data
   use pmc_gas_state
 #ifdef PMC_USE_TCHEM
   use iso_c_binding
-  use tchemdriver
 #endif
   use pmc_util, only : die_msg, warn_assert_msg, assert_msg
 
+interface
+  subroutine initialize(arg_chemfile) bind(c, name="initialize")
+    use iso_c_binding
+    character(kind=c_char) :: arg_chemfile(*)
+  end subroutine initialize
+  subroutine finalize() bind(c, name="finalize")
+  end subroutine finalize
+  function TChem_getNumberOfSpecies() bind(c, name="TChem_getNumberOfSpecies")
+    use iso_c_binding
+    integer(kind=c_int) :: TChem_getNumberOfSpecies
+  end function
+  function TChem_getLengthOfStateVector() bind(c, &
+      name="TChem_getLengthOfStateVector")
+    use iso_c_binding
+    integer(kind=c_int) :: TChem_getLengthOfStateVector
+  end function
+  subroutine TChem_getStateVector(array) bind(c, name="TChem_getStateVector")
+    use iso_c_binding
+    real(kind=c_double) :: array(*)
+  end subroutine
+  subroutine TChem_setStateVector(array) bind(c, name="TChem_setStateVector")
+    use iso_c_binding
+    real(kind=c_double) :: array(*)
+  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
+  subroutine TChem_doTimestep() bind(C, name="TChem_doTimestep")
+  end subroutine
+end interface
+
 contains
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -28,63 +62,67 @@ module pmc_tchem_interface
   subroutine pmc_tchem_interface_solve(env_state, aero_data, aero_state, &
        gas_data, gas_state)
 
+    !> Environment data.
     type(env_state_t), intent(in) :: env_state
+    !> Aerosol data.
     type(aero_data_t), intent(in) :: aero_data
+    !> Aerosol state.
     type(aero_state_t), intent(inout) :: aero_state
+    !> Gas data.
     type(gas_data_t), intent(in) :: gas_data
+    !> Gas state.
     type(gas_state_t), intent(inout) :: gas_state
 
-    call tchem_from_partmc(gas_data, gas_state)
-    call timestep()
-    call tchem_to_partmc(gas_data, gas_state)
+    call tchem_from_partmc(gas_data, gas_state, env_state)
+    call tchem_timestep()
+    call tchem_to_partmc(gas_data, gas_state, env_state)
 
   end subroutine pmc_tchem_interface_solve
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-  subroutine pmc_tchem_initialize(config_filename, gas_data, gas_state, aero_data)
-
+  !> Initialize TChem and PartMC gas and aerosol data.
+  subroutine pmc_tchem_initialize(config_filename, gas_data, aero_data)
     use iso_c_binding
 
+    !>
     character(len=*), intent(in) :: config_filename
+    !> Gas data.
     type(gas_data_t), intent(inout) :: gas_data
-    type(gas_state_t), intent(inout) :: gas_state
+    !> Aerosol data.
     type(aero_data_t), intent(inout) :: aero_data
-    integer(c_int) :: nSpec
-    integer :: i, j
 
-    character(len=5) :: name
-    character(kind=c_char), dimension(5) :: name_data
+    integer(kind=c_int) :: nSpec
+    integer :: i
     real(kind=c_double), dimension(:), allocatable :: array 
+    character(:), allocatable ::  val
 
     ! initialize the model
     call TChem_initialize(trim(config_filename))
 
-    ! feedback to PartMC
     ! what we need:
-    !  - number of gas species
+    !  - number of gas species for gas_data
     !  - name of gas species to set gas_data
-    !  - set the size of gas_state
     nSpec = TChem_getNumberOfSpecies() 
     call ensure_string_array_size(gas_data%name, nSpec)
 
-    call gas_state_set_size(gas_state, nSpec)
-
-
-    call TChem_getSpeciesNames()
     ! name of gas species
-    !do i = 1,nSpec
-    !   gas_data%name(i) = "H2O"
-    !end do   
+    do i = 1,nSpec
+       val = TChem_species_name(i-1) 
+       gas_data%name(i) = trim(val)
+    end do   
+
+    ! We need this
+    allocate(gas_data%mosaic_index(gas_data_n_spec(gas_data)))
+    gas_data%mosaic_index(:) = 0
 
     print*, 'in partmc', nSpec, trim(config_filename), gas_data_n_spec(gas_data)
-
-    call tchem_to_partmc(gas_data, gas_state)
-    print*, 'in partmc'
-    print*, gas_state%mix_rat 
 
-end subroutine
+  end subroutine
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  !> Clean up TChem.
   subroutine pmc_tchem_cleanup()
 
     call finalize()
@@ -93,41 +131,112 @@ end subroutine pmc_tchem_cleanup
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-  !>
-  subroutine tchem_to_partmc(gas_data, gas_state)
+  !> Map all data TChem -> PartMC.
+  subroutine tchem_to_partmc(gas_data, gas_state, env_state)
 
-    !>
+    !> Gas data.
     type(gas_data_t), intent(in) :: gas_data
-    !>
+    !> Gas state.
     type(gas_state_t), intent(inout) :: gas_state
+    !> Environment state.
+    type(env_state_t), intent(in) :: env_state
 
-    integer(c_int) :: nSpec
-    real(kind=c_double), dimension(:), allocatable :: array
+    integer(c_int) :: nSpec, stateVecDim
+    real(kind=c_double), dimension(:), allocatable :: stateVector 
 
     ! Get gas array
+    stateVecDim = TChem_getLengthOfStateVector()
     nSpec = TChem_getNumberOfSpecies()
-    allocate(array(nSpec))
+    allocate(stateVector(stateVecDim))
     array = 0.0d0
-    call TChem_getStateVector(array)
-    gas_state%mix_rat = array
+    call TChem_getStateVector(stateVector)
+
+    ! FIXME: adjust this range later
+    gas_state%mix_rat = stateVector(4:nSpec+3)
+    ! Convert gas_state from mol m^-3 to ppb.
+    call gas_state_mole_dens_to_ppb(gas_state, env_state)
 
   end subroutine tchem_to_partmc
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-  !>
-  subroutine tchem_from_partmc(gas_data, gas_state)
+  !> Map all data PartMC -> TChem.
+  subroutine tchem_from_partmc(gas_data, gas_state, env_state)
 
-    !>
+    !> Gas data.
     type(gas_data_t), intent(in) :: gas_data
-    ! FIXME: Fix intent later
+    !> Gas State.
     type(gas_state_t), intent(inout) :: gas_state
+    !> Environment state.
+    type(env_state_t), intent(in) :: env_state
+
+    real(kind=dp), allocatable :: stateVector(:)
+    integer :: stateVecDim
+
+    ! Make a state vector
+    ! Density
+    ! Pressure
+    ! Temperature
+    ! Concentrations
+    ! get size of stateVector
+    stateVecDim = TChem_getLengthOfStateVector()
+    allocate(stateVector(stateVecDim))
+    stateVector(1) = env_state_air_den(env_state)
+    stateVector(2) = env_state%pressure
+    stateVector(3) = env_state%temp
+
+    call gas_state_ppb_to_mole_dens(gas_state, env_state)
+    stateVector(4:gas_data_n_spec(gas_data)+3) = gas_state%mix_rat
 
-    call TChem_setStateVector(gas_state%mix_rat)
+    call TChem_setStateVector(stateVector)
 
   end subroutine tchem_from_partmc
-#endif
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  !> Do a single timestep of TChem chemistry.
+  subroutine TChem_timestep()
+    use iso_c_binding
+
+    call TChem_doTimestep()
+
+  end subroutine TChem_timestep
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  !> Initialize TChem.
+  subroutine TChem_initialize(chemFile)
+    use iso_c_binding
+
+    !> Chemistry configuration file.
+    character(kind=c_char,len=*), intent(in) :: chemFile
+
+    print*, 'in TChemDriver: ', chemFile
+
+    call initialize(chemFile//c_null_char)
+
+  end subroutine TChem_initialize
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  !> Get i_spec species name from TChem.
+  function TChem_species_name(i_spec) result(species_name)
+    use iso_c_binding
+    ! 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)
+    allocate(character(N) :: species_name)
+    species_name = cbuf(:N)
+
+  end function TChem_species_name
+#endif
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 end module pmc_tchem_interface

test/tchem/README

@@ -0,0 +1,3 @@
+
+TChemn Test-case
+===========================================

test/tchem/aero_back.dat

@@ -0,0 +1,6 @@
+# time (s)
+# rate (s^{-1})
+# aerosol distribution filename
+time    0
+rate    3e-5
+dist    aero_back_dist.dat

test/tchem/aero_back_comp.dat

@@ -0,0 +1,2 @@
+#    mass fractions
+SB                1

test/tchem/aero_back_dist.dat

@@ -0,0 +1,6 @@
+mode_name back1                 # name of mode source
+mass_frac aero_back_comp.dat    # species mass fractions
+diam_type geometric             # type of diameter specified
+mode_type mono                  # type of distribution
+num_conc 1e9                    # particle number concentration (#/m^3)
+diam 1e-4                       # diameter (m)

test/tchem/aero_data.dat

@@ -0,0 +1,4 @@
+#     dens (kg/m^3)   ions in soln (1)    molec wght (kg/mole)   kappa (1)
+SI            10000                0                   18d-3      0
+SB              100                0                   18d-3      0
+SE             1500                0                   18d-3      0

test/tchem/aero_emit.dat

@@ -0,0 +1,6 @@
+# time (s)
+# rate (s^{-1})
+# aerosol distribution filename
+time    0
+rate    1.5e-5
+dist    aero_emit_dist.dat

test/tchem/aero_emit_comp.dat

@@ -0,0 +1,2 @@
+#    mass fractions
+SE                1

test/tchem/aero_emit_dist.dat

@@ -0,0 +1,6 @@
+mode_name emit1                 # name of mode source
+mass_frac aero_emit_comp.dat    # species mass fractions
+diam_type geometric             # type of diameter specified
+mode_type mono                  # type of distribution
+num_conc 1e12                   # particle number concentration (#/m^2)
+diam 6.03e-5                       # diameter (m)

test/tchem/aero_init_comp.dat

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

test/tchem/aero_init_dist.dat

@@ -0,0 +1,6 @@
+mode_name init1                 # name of mode source
+mass_frac aero_init_comp.dat    # species mass fractions
+diam_type geometric             # type of diameter specified
+mode_type mono                  # type of distribution
+num_conc 1e9                    # particle number concentration (#/m^3)
+diam 3.03e-5                       # diameter (m)