Add Chapman and Terminator configurations for chemistry

.gitignore

@@ -0,0 +1,7 @@
+**/*.DS_Store
+**/*.pyc
+**/build/
+**/CMakeCache.txt
+**/CMakeFiles/
+.vscode
+xcode/

schemes/musica/micm/musica_ccpp_micm.F90

@@ -2,17 +2,18 @@ module musica_ccpp_micm
 
   ! Note: "micm_t" is included in an external pre-built MICM library that the host
   ! model is responsible for linking to during compilation
-  use musica_micm,          only: micm_t
+  use ccpp_kinds,           only: kind_phys
   use musica_ccpp_util,     only: has_error_occurred
   use musica_ccpp_namelist, only: filename_of_micm_configuration
-  use ccpp_kinds,           only: kind_phys
+  use musica_micm,          only: micm_t
 
   implicit none
   private
 
-  public :: micm_register, micm_init, micm_run, micm_final
+  public :: micm_register, micm_init, micm_run, micm_final, micm, number_of_rate_parameters
 
   type(micm_t), pointer  :: micm => null( )
+  integer :: number_of_rate_parameters = 0
 
 contains
 
@@ -36,7 +37,7 @@ subroutine micm_register(solver_type, num_grid_cells, constituent_props, errmsg,
     logical                       :: is_advected
     integer                       :: i, species_index
 
-    micm => micm_t(filename_of_micm_configuration, solver_type, num_grid_cells, error)
+    micm => micm_t(trim(filename_of_micm_configuration), solver_type, num_grid_cells, error)
     if (has_error_occurred(error, errmsg, errcode)) return
 
     allocate(constituent_props(micm%species_ordering%size()), stat=errcode)
@@ -73,6 +74,7 @@ subroutine micm_register(solver_type, num_grid_cells, constituent_props, errmsg,
       if (errcode /= 0) return
     end associate ! map
     end do
+    number_of_rate_parameters = micm%user_defined_reaction_rates%size()
 
   end subroutine micm_register
 
@@ -91,34 +93,31 @@ subroutine micm_run(time_step, temperature, pressure, dry_air_density, &
                       user_defined_rate_parameters, constituents, errmsg, errcode)
     use musica_micm,   only: solver_stats_t
     use musica_util,   only: string_t, error_t
-    use iso_c_binding, only: c_double
+    use iso_c_binding, only: c_double, c_loc
 
-    real(kind_phys),    intent(in)    :: time_step          ! s
-    real(c_double),     intent(in)    :: temperature(:)     ! K
-    real(c_double),     intent(in)    :: pressure(:)        ! Pa
-    real(c_double),     intent(in)    :: dry_air_density(:) ! kg m-3
-    real(c_double),     intent(in)    :: user_defined_rate_parameters(:) ! various units
-    real(c_double),     intent(inout) :: constituents(:)    ! mol m-3
-    character(len=512), intent(out)   :: errmsg
-    integer,            intent(out)   :: errcode
+    real(kind_phys),         intent(in)    :: time_step            ! s
+    real(kind_phys), target, intent(in)    :: temperature(:,:)     ! K
+    real(kind_phys), target, intent(in)    :: pressure(:,:)        ! Pa
+    real(kind_phys), target, intent(in)    :: dry_air_density(:,:) ! kg m-3
+    real(kind_phys), target, intent(in)    :: user_defined_rate_parameters(:,:,:) ! various units
+    real(kind_phys), target, intent(inout) :: constituents(:,:,:)  ! mol m-3
+    character(len=512),      intent(out)   :: errmsg
+    integer,                 intent(out)   :: errcode
 
     ! local variables
     type(string_t)       :: solver_state
     type(solver_stats_t) :: solver_stats
     type(error_t)        :: error
-    real(c_double)       :: c_time_step
     integer              :: i_elem
 
-    c_time_step = real(time_step, c_double) 
-
-    call micm%solve(c_time_step,                  &
-                    temperature,                  &
-                    pressure,                     &
-                    dry_air_density,              &
-                    constituents,                 &
-                    user_defined_rate_parameters, &
-                    solver_state,                 &
-                    solver_stats,                 &
+    call micm%solve(real(time_step, kind=c_double),      &
+                    c_loc(temperature),                  &
+                    c_loc(pressure),                     &
+                    c_loc(dry_air_density),              &
+                    c_loc(constituents),                 &
+                    c_loc(user_defined_rate_parameters), &
+                    solver_state,                        &
+                    solver_stats,                        &
                     error)
     if (has_error_occurred(error, errmsg, errcode)) return
 

schemes/musica/micm/musica_ccpp_micm_util.F90

@@ -2,80 +2,10 @@ module musica_ccpp_micm_util
   implicit none
 
   private
-  public :: reshape_into_micm_arr, reshape_into_ccpp_arr, convert_to_mol_per_cubic_meter, &
-            convert_to_mass_mixing_ratio
+  public :: convert_to_mol_per_cubic_meter, convert_to_mass_mixing_ratio
 
 contains
 
-  !> Reshape array (2D/3D -> 1D) and convert type (kind_phys -> c_double)
-  subroutine reshape_into_micm_arr(temperature, pressure, dry_air_density, constituents,  &
-                                   micm_temperature, micm_pressure, micm_dry_air_density, &
-                                   micm_constituents)
-    use iso_c_binding, only: c_double
-    use ccpp_kinds,    only: kind_phys
-
-    real(kind_phys), intent(in)  :: temperature(:,:)        ! K
-    real(kind_phys), intent(in)  :: pressure(:,:)           ! Pa
-    real(kind_phys), intent(in)  :: dry_air_density(:,:)    ! kg m-3
-    real(kind_phys), intent(in)  :: constituents(:,:,:)     ! kg kg-1
-    real(c_double),  intent(out) :: micm_temperature(:)     ! K
-    real(c_double),  intent(out) :: micm_pressure(:)        ! Pa
-    real(c_double),  intent(out) :: micm_dry_air_density(:) ! kg m-3
-    real(c_double),  intent(out) :: micm_constituents(:)    ! kg kg-1
-
-    ! local variables
-    integer :: num_columns, num_layers, num_constituents
-    integer :: i_column, i_layer, i_elem, i_constituents
-
-    num_columns = size(constituents, dim=1)
-    num_layers = size(constituents, dim=2)
-    num_constituents = size(constituents, dim=3)
-
-    ! Reshape into 1-D arry in species-column first order, referring to
-    ! state.variables_[i_cell][i_species] = concentrations[i_species_elem++]
-    i_elem = 1
-    i_constituents = 1
-    do i_layer = 1, num_layers
-      do i_column = 1, num_columns
-        micm_temperature(i_elem) = real(temperature(i_column, i_layer), c_double)
-        micm_pressure(i_elem) = real(pressure(i_column, i_layer), c_double)
-        micm_dry_air_density(i_elem) = real(dry_air_density(i_column, i_layer), c_double)
-        micm_constituents(i_constituents : i_constituents + num_constituents - 1) &
-                            = real(constituents(i_column, i_layer, :), c_double)
-        i_elem = i_elem + 1
-        i_constituents = i_constituents + num_constituents
-      end do
-    end do
-
-  end subroutine reshape_into_micm_arr
-
-  !> Reshape array (1D -> 3D) and convert type (c_double -> kind_phys)
-  subroutine reshape_into_ccpp_arr(micm_constituents, constituents)
-    use iso_c_binding, only: c_double
-    use ccpp_kinds,    only: kind_phys
-
-    real(c_double),  intent(in)    :: micm_constituents(:) ! kg kg-1
-    real(kind_phys), intent(inout) :: constituents(:,:,:)  ! kg kg-1
-
-    ! local variables
-    integer :: num_columns, num_layers, num_constituents
-    integer :: i_column, i_layer, i_constituents
-
-    num_columns = size(constituents, dim=1)
-    num_layers = size(constituents, dim=2)
-    num_constituents = size(constituents, dim=3)
-
-    i_constituents = 1
-    do i_layer = 1, num_layers
-      do i_column = 1, num_columns
-        constituents(i_column, i_layer, :) &
-            = real(micm_constituents(i_constituents : i_constituents + num_constituents - 1), kind_phys)
-        i_constituents = i_constituents + num_constituents
-      end do
-    end do
-
-  end subroutine reshape_into_ccpp_arr
-
   ! Convert CAM-SIMA unit to MICM unit (kg kg-1  ->  mol m-3)
   subroutine convert_to_mol_per_cubic_meter(dry_air_density, molar_mass_arr, constituents)
     use ccpp_kinds,    only: kind_phys

schemes/musica/musica_ccpp.F90

@@ -1,7 +1,9 @@
 !> Top-level wrapper for MUSICA chemistry components
 module musica_ccpp
-  use musica_ccpp_micm, only: micm_register, micm_init, micm_run, micm_final
-  use musica_ccpp_tuvx, only: tuvx_init, tuvx_run, tuvx_final
+  use musica_ccpp_micm,     only: micm_register, micm_init, micm_run, micm_final
+  use musica_ccpp_namelist, only: filename_of_tuvx_micm_mapping_configuration
+  use musica_ccpp_tuvx,     only: tuvx_init, tuvx_run, tuvx_final
+  use musica_util,          only: index_mappings_t
 
   implicit none
   private
@@ -30,18 +32,20 @@ end subroutine musica_ccpp_register
   subroutine musica_ccpp_init(vertical_layer_dimension, vertical_interface_dimension, &
                               photolysis_wavelength_grid_interfaces, errmsg, errcode)
     use ccpp_kinds, only : kind_phys
-
+    use musica_ccpp_micm, only: micm
+    use musica_ccpp_util, only: has_error_occurred
     integer,            intent(in)  :: vertical_layer_dimension                 ! (count)
     integer,            intent(in)  :: vertical_interface_dimension             ! (count)
     real(kind_phys),    intent(in)  :: photolysis_wavelength_grid_interfaces(:) ! m
     character(len=512), intent(out) :: errmsg
     integer,            intent(out) :: errcode
 
-    call tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
-                   photolysis_wavelength_grid_interfaces, errmsg, errcode)
-    if (errcode /= 0) return
     call micm_init(errmsg, errcode)
     if (errcode /= 0) return
+    call tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
+                   photolysis_wavelength_grid_interfaces, &
+                   micm%user_defined_reaction_rates, errmsg, errcode)
+    if (errcode /= 0) return
 
   end subroutine musica_ccpp_init
 
@@ -56,53 +60,42 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
                              geopotential_height_wrt_surface_at_interface, surface_temperature,    &
                              surface_geopotential, surface_albedo,                                 &
                              standard_gravitational_acceleration, errmsg, errcode)
-    use musica_ccpp_micm_util,     only: reshape_into_micm_arr, reshape_into_ccpp_arr
-    use musica_ccpp_micm_util,     only: convert_to_mol_per_cubic_meter, convert_to_mass_mixing_ratio
     use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
     use ccpp_kinds,                only: kind_phys
-    use iso_c_binding,             only: c_double
+    use musica_ccpp_micm,          only: number_of_rate_parameters
+    use musica_ccpp_micm_util,     only: convert_to_mol_per_cubic_meter, convert_to_mass_mixing_ratio
 
-    real(kind_phys),    intent(in)    :: time_step                                         ! s
-    real(kind_phys),    intent(in)    :: temperature(:,:)                                  ! K
-    real(kind_phys),    intent(in)    :: pressure(:,:)                                     ! Pa
-    real(kind_phys),    intent(in)    :: dry_air_density(:,:)                              ! kg m-3
+    real(kind_phys),         intent(in)    :: time_step                                         ! s
+    real(kind_phys), target, intent(in)    :: temperature(:,:)                                  ! K
+    real(kind_phys), target, intent(in)    :: pressure(:,:)                                     ! Pa
+    real(kind_phys), target, intent(in)    :: dry_air_density(:,:)                              ! kg m-3
     type(ccpp_constituent_prop_ptr_t), &
-                        intent(in)    :: constituent_props(:)
-    real(kind_phys),    intent(inout) :: constituents(:,:,:)                               ! kg kg-1
-    real(kind_phys),    intent(in)    :: geopotential_height_wrt_surface_at_midpoint(:,:)  ! m (column, layer)
-    real(kind_phys),    intent(in)    :: geopotential_height_wrt_surface_at_interface(:,:) ! m (column, interface)
-    real(kind_phys),    intent(in)    :: surface_temperature(:)                            ! K
-    real(kind_phys),    intent(in)    :: surface_geopotential(:)                           ! m2 s-2
-    real(kind_phys),    intent(in)    :: surface_albedo                                    ! unitless
-    real(kind_phys),    intent(in)    :: standard_gravitational_acceleration               ! m s-2
-    character(len=512), intent(out)   :: errmsg
-    integer,            intent(out)   :: errcode
+                             intent(in)    :: constituent_props(:)
+    real(kind_phys), target, intent(inout) :: constituents(:,:,:)                               ! kg kg-1
+    real(kind_phys),         intent(in)    :: geopotential_height_wrt_surface_at_midpoint(:,:)  ! m (column, layer)
+    real(kind_phys),         intent(in)    :: geopotential_height_wrt_surface_at_interface(:,:) ! m (column, interface)
+    real(kind_phys),         intent(in)    :: surface_temperature(:)                            ! K
+    real(kind_phys),         intent(in)    :: surface_geopotential(:)                           ! m2 s-2
+    real(kind_phys),         intent(in)    :: surface_albedo                                    ! unitless
+    real(kind_phys),         intent(in)    :: standard_gravitational_acceleration               ! m s-2
+    character(len=512),      intent(out)   :: errmsg
+    integer,                 intent(out)   :: errcode
 
     ! local variables
-    real(c_double), dimension(size(temperature, dim=1)     &
-                            * size(temperature, dim=2))      :: micm_temperature
-    real(c_double), dimension(size(pressure, dim=1)        &
-                            * size(pressure, dim=2))         :: micm_pressure
-    real(c_double), dimension(size(dry_air_density, dim=1) &
-                            * size(dry_air_density, dim=2))  :: micm_dry_air_density
-    real(c_double), dimension(size(constituents, dim=1)    &
-                            * size(constituents, dim=2)    &
-                            * size(constituents, dim=3))     :: micm_constituents ! mol m-3
     real(kind_phys), dimension(size(constituents, dim=3))    :: molar_mass_arr    ! kg mol-1
-
-    ! temporarily dimensioned to Chapman mechanism until mapping between MICM and TUV-x is implemented
-    real(c_double), dimension(size(constituents, dim=1) &
-                            * size(constituents, dim=2) &
-                            * 3) :: photolysis_rate_constants ! s-1
+    real(kind_phys), dimension(size(constituents, dim=1), &
+                               size(constituents, dim=2), &
+                               number_of_rate_parameters)    :: rate_parameters ! various units
     integer :: i_elem
 
+    ! Calculate photolysis rate constants using TUV-x
     call tuvx_run(temperature, dry_air_density,                 &
                   geopotential_height_wrt_surface_at_midpoint,  &
                   geopotential_height_wrt_surface_at_interface, &
                   surface_temperature, surface_geopotential,    &
                   surface_albedo,                               &
                   standard_gravitational_acceleration,          &
-                  photolysis_rate_constants,                    &
+                  rate_parameters,                              &
                   errmsg, errcode)
 
     ! Get the molar mass that is set in the call to instantiate()
@@ -127,16 +120,9 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
     ! Convert CAM-SIMA unit to MICM unit (kg kg-1  ->  mol m-3)
     call convert_to_mol_per_cubic_meter(dry_air_density, molar_mass_arr, constituents)
 
-    ! Reshape array (2D/3D -> 1D) and convert type (kind_phys -> c_double)
-    call reshape_into_micm_arr(temperature, pressure, dry_air_density, constituents, &
-                      micm_temperature, micm_pressure, micm_dry_air_density, micm_constituents)
-
-    ! temporarily pass in unmapped photolysis rate constants until mapping between MICM and TUV-x is implemented
-    call micm_run(time_step, micm_temperature, micm_pressure, micm_dry_air_density, &
-                  photolysis_rate_constants, micm_constituents, errmsg, errcode)
-
-    ! Reshape array (1D -> 3D) and convert type (c_double -> kind_phys)
-    call reshape_into_ccpp_arr(micm_constituents, constituents)
+    ! Solve chemistry at the current time step
+    call micm_run(time_step, temperature, pressure, dry_air_density, rate_parameters, &
+                  constituents, errmsg, errcode)
 
     ! Convert MICM unit back to CAM-SIMA unit (mol m-3  ->  kg kg-1)
     call convert_to_mass_mixing_ratio(dry_air_density, molar_mass_arr, constituents)

schemes/musica/musica_ccpp_namelist.xml

@@ -103,5 +103,17 @@
       <value>UNSET_PATH</value>
     </values>
   </entry>
+  <entry id="filename_of_tuvx_micm_mapping_configuration">
+    <type>char*512</type>
+    <category>musica_ccpp</category>
+    <group>musica_ccpp</group>
+    <standard_name>filename_of_tuvx_micm_mapping_configuration</standard_name>
+    <units>none</units>
+    <desc>
+     A configuration file for the mapping of TUV-x photolysis rates to MICM custom rate parameters
+    </desc>
+    <values>
+      <value>UNSET_PATH</value>
+    </values>
 
 </entry_id_pg>