Set gas-species profiles in TUV-x and map indices between constituents and MICM

schemes/musica/micm/musica_ccpp_micm.F90

@@ -19,15 +19,16 @@ module musica_ccpp_micm
 
   !> Registers MICM constituent properties with the CCPP
   subroutine micm_register(solver_type, number_of_grid_cells, constituent_props, &
-                           errmsg, errcode)
+                           micm_species, errmsg, errcode)
     use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t
-    use musica_micm,               only: Rosenbrock, RosenbrockStandardOrder
+    use musica_ccpp_species,       only: musica_species_t
     use musica_util,               only: error_t
     use iso_c_binding,             only: c_int
 
     integer(c_int),                                   intent(in)  :: solver_type
     integer(c_int),                                   intent(in)  :: number_of_grid_cells
     type(ccpp_constituent_properties_t), allocatable, intent(out) :: constituent_props(:)
+    type(musica_species_t),              allocatable, intent(out) :: micm_species(:)
     character(len=512),                               intent(out) :: errmsg
     integer,                                          intent(out) :: errcode
 
@@ -36,6 +37,7 @@ subroutine micm_register(solver_type, number_of_grid_cells, constituent_props, &
     real(kind=kind_phys)          :: molar_mass
     character(len=:), allocatable :: species_name
     logical                       :: is_advected
+    integer                       :: number_of_species
     integer                       :: i, species_index
 
     if (associated( micm )) then
@@ -46,13 +48,20 @@ subroutine micm_register(solver_type, number_of_grid_cells, constituent_props, &
                    number_of_grid_cells, error)
     if (has_error_occurred(error, errmsg, errcode)) return
 
-    allocate(constituent_props(micm%species_ordering%size()), stat=errcode)
+    number_of_species = micm%species_ordering%size()
+    allocate(constituent_props(number_of_species), stat=errcode)
     if (errcode /= 0) then
       errmsg = "[MUSICA Error] Failed to allocate memory for constituent properties."
       return
     end if
 
-    do i = 1, micm%species_ordering%size()
+    allocate(micm_species(number_of_species), stat=errcode)
+    if (errcode /= 0) then
+      errmsg = "[MUSICA Error] Failed to allocate memory for micm species."
+      return
+    end if
+
+    do i = 1, number_of_species
     associate( map => micm%species_ordering )
       species_name = map%name(i)
       species_index = map%index(i)
@@ -78,6 +87,13 @@ subroutine micm_register(solver_type, number_of_grid_cells, constituent_props, &
         errcode = errcode, &
         errmsg = errmsg)
       if (errcode /= 0) return
+
+      ! Species are ordered to match the sequence of the MICM state array
+      micm_species(species_index) = musica_species_t( &
+        name = species_name, &
+        unit = 'kg kg-1', &
+        molar_mass = molar_mass, &
+        index_musica_species = species_index )
     end associate ! map
     end do
     number_of_rate_parameters = micm%user_defined_reaction_rates%size()

schemes/musica/musica_ccpp.F90

@@ -14,16 +14,20 @@ module musica_ccpp
 
   !> \section arg_table_musica_ccpp_register Argument Table
   !! \htmlinclude musica_ccpp_register.html
-  subroutine musica_ccpp_register(constituent_props, errmsg, &
-                                  errcode)
-    use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t
-    use musica_ccpp_namelist,      only: micm_solver_type
+  subroutine musica_ccpp_register(constituent_props, errmsg, errcode)
+    use ccpp_constituent_prop_mod,     only: ccpp_constituent_properties_t
+    use musica_ccpp_namelist,          only: micm_solver_type
+    use musica_ccpp_species,           only: musica_species_t, register_musica_species
+    use musica_ccpp_tuvx_load_species, only: check_tuvx_species_initialization
 
     type(ccpp_constituent_properties_t), allocatable, intent(out) :: constituent_props(:)
     character(len=512),                               intent(out) :: errmsg
     integer,                                          intent(out) :: errcode
 
+    ! local variables
     type(ccpp_constituent_properties_t), allocatable :: constituent_props_subset(:)
+    type(musica_species_t),              allocatable :: micm_species(:)
+    type(musica_species_t),              allocatable :: tuvx_species(:)
     integer :: number_of_grid_cells
 
     ! Temporary fix until the number of grid cells is only needed to create a MICM state
@@ -32,50 +36,67 @@ subroutine musica_ccpp_register(constituent_props, errmsg, &
     ! the solver when the number of grid cells is known at the init stage.
     number_of_grid_cells = 1
     call micm_register(micm_solver_type, number_of_grid_cells, constituent_props_subset, &
-                       errmsg, errcode)
+                       micm_species, errmsg, errcode)
     if (errcode /= 0) return
     constituent_props = constituent_props_subset
     deallocate(constituent_props_subset)
 
-    call tuvx_register(constituent_props_subset, errmsg, errcode)
+    call tuvx_register(micm_species, tuvx_species, constituent_props_subset, &
+                       errmsg, errcode)
     if (errcode /= 0) return
     constituent_props = [ constituent_props, constituent_props_subset ]
 
+    call register_musica_species(micm_species, tuvx_species)
+    call check_tuvx_species_initialization(errmsg, errcode)
+    if (errcode /= 0) return
+
   end subroutine musica_ccpp_register
 
   !> \section arg_table_musica_ccpp_init Argument Table
   !! \htmlinclude musica_ccpp_init.html
   subroutine musica_ccpp_init(horizontal_dimension, vertical_layer_dimension, &
                               vertical_interface_dimension, &
                               photolysis_wavelength_grid_interfaces, &
-                              constituent_props, errmsg, errcode)
-    use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t, ccpp_constituent_prop_ptr_t
-    use ccpp_kinds, only : kind_phys
-    use musica_ccpp_micm, only: micm
-    use musica_ccpp_namelist, only: micm_solver_type
-    use musica_ccpp_util, only: has_error_occurred
+                              constituent_props_ptr, errmsg, errcode)
+    use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t, ccpp_constituent_properties_t
+    use ccpp_kinds,                only: kind_phys
+    use musica_ccpp_micm,          only: micm
+    use musica_ccpp_namelist,      only: micm_solver_type
+    use musica_ccpp_util,          only: has_error_occurred
+    use musica_ccpp_species,       only: initialize_musica_species_indices, initialize_molar_mass_array, &
+                                         check_initialization, musica_species_t
+
     integer,                           intent(in)  :: horizontal_dimension                     ! (count)
     integer,                           intent(in)  :: vertical_layer_dimension                 ! (count)
     integer,                           intent(in)  :: vertical_interface_dimension             ! (count)
     real(kind_phys),                   intent(in)  :: photolysis_wavelength_grid_interfaces(:) ! m
-    type(ccpp_constituent_prop_ptr_t), intent(in)  :: constituent_props(:)
+    type(ccpp_constituent_prop_ptr_t), intent(in)  :: constituent_props_ptr(:)
     character(len=512),                intent(out) :: errmsg
     integer,                           intent(out) :: errcode
 
-    integer :: number_of_grid_cells
-    type(ccpp_constituent_properties_t), allocatable :: micm_species_props(:)
+    ! local variables
+    type(ccpp_constituent_properties_t), allocatable :: constituent_props(:)
+    type(musica_species_t),              allocatable :: micm_species(:)
+    integer                                          :: number_of_grid_cells
 
     ! Temporary fix until the number of grid cells is only needed to create a MICM state
     ! instead of when the solver is created.
     ! Re-create the MICM solver with the correct number of grid cells
     number_of_grid_cells = horizontal_dimension * vertical_layer_dimension
-    call micm_register(micm_solver_type, number_of_grid_cells, micm_species_props, errmsg, errcode)
+    call micm_register(micm_solver_type, number_of_grid_cells, constituent_props, &
+                       micm_species, errmsg, errcode)
     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, &
-                   constituent_props, errmsg, errcode)
+                   photolysis_wavelength_grid_interfaces,                  &
+                   micm%user_defined_reaction_rates, errmsg, errcode)
+    if (errcode /= 0) return
+
+    call initialize_musica_species_indices(constituent_props_ptr, errmsg, errcode)
+    if (errcode /= 0) return
+    call initialize_molar_mass_array(constituent_props_ptr, errmsg, errcode)
+    if (errcode /= 0) return
+    call check_initialization(errmsg, errcode)
     if (errcode /= 0) return
 
   end subroutine musica_ccpp_init
@@ -98,6 +119,9 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
     use ccpp_kinds,                only: kind_phys
     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
+    use musica_ccpp_species,       only: number_of_micm_species, number_of_tuvx_species, &
+      micm_indices_constituent_props, tuvx_indices_constituent_props, micm_molar_mass_array, &
+      extract_subset_constituents, update_constituents
 
     real(kind_phys),         intent(in)    :: time_step                                         ! s
     real(kind_phys), target, intent(in)    :: temperature(:,:)                                  ! K
@@ -122,69 +146,71 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
     integer,                 intent(out)   :: errcode
 
     ! local variables
-    real(kind_phys), dimension(size(constituents, dim=3))    :: molar_mass_arr    ! kg mol-1
     real(kind_phys), dimension(size(constituents, dim=1), &
                                size(constituents, dim=2), &
-                               number_of_rate_parameters)    :: rate_parameters ! various units
-    integer :: i_elem
+                               number_of_rate_parameters) :: rate_parameters ! various units
+    real(kind_phys), dimension(size(constituents, dim=1), &
+                               size(constituents, dim=2), &
+                               number_of_micm_species)    :: constituents_micm_species ! kg kg-1
+    real(kind_phys), dimension(size(constituents, dim=1), &
+                               size(constituents, dim=2), &
+                               number_of_tuvx_species)    :: constituents_tuvx_species ! kg kg-1
+
+    call extract_subset_constituents(tuvx_indices_constituent_props, constituents, &
+                                     constituents_tuvx_species, errmsg, errcode)
+    if (errcode /= 0) return
 
     ! Calculate photolysis rate constants using TUV-x
-    call tuvx_run(temperature, dry_air_density,                  &
-                  constituents,                                  &
-                  geopotential_height_wrt_surface_at_midpoint,   &
-                  geopotential_height_wrt_surface_at_interface,  &
-                  surface_geopotential, surface_temperature,     &
-                  surface_albedo,                                &
-                  photolysis_wavelength_grid_interfaces,         &
-                  extraterrestrial_flux,                         &
-                  standard_gravitational_acceleration,           &
-                  cloud_area_fraction,                           &
-                  air_pressure_thickness,                        &
-                  solar_zenith_angle,                            &
-                  earth_sun_distance,                            &
-                  rate_parameters,                               &
+    call tuvx_run(temperature, dry_air_density,                 &
+                  constituents_tuvx_species,                    &
+                  geopotential_height_wrt_surface_at_midpoint,  &
+                  geopotential_height_wrt_surface_at_interface, &
+                  surface_geopotential, surface_temperature,    &
+                  surface_albedo,                               &
+                  photolysis_wavelength_grid_interfaces,        &
+                  extraterrestrial_flux,                        &
+                  standard_gravitational_acceleration,          &
+                  cloud_area_fraction,                          &
+                  air_pressure_thickness,                       &
+                  solar_zenith_angle,                           &
+                  earth_sun_distance,                           &
+                  rate_parameters,                              &
                   errmsg, errcode)
 
-    ! Get the molar mass that is set in the call to instantiate()
-    do i_elem = 1, size(molar_mass_arr)
-      call constituent_props(i_elem)%molar_mass(molar_mass_arr(i_elem), errcode, errmsg)
-      if (errcode /= 0) then
-        errmsg = "[MUSICA Error] Unable to get molar mass."
-        return
-      end if
-    end do
-
-    ! TODO(jiwon) Check molar mass is non zero as it becomes a denominator for unit converison
-    ! this code will be deleted when the framework does the check
-    do i_elem = 1, size(molar_mass_arr)
-      if (molar_mass_arr(i_elem) <= 0) then
-        errcode = 1
-        errmsg = "[MUSICA Error] Molar mass must be greater than zero."
-        return
-      end if
-    end do
+    call update_constituents(tuvx_indices_constituent_props, constituents_tuvx_species, &
+                             constituents, errmsg, errcode)
+    if (errcode /= 0) return
+    call extract_subset_constituents(micm_indices_constituent_props, constituents, &
+                                     constituents_micm_species, errmsg, errcode)
+    if (errcode /= 0) return
 
     ! 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)
+    call convert_to_mol_per_cubic_meter(dry_air_density, micm_molar_mass_array, constituents_micm_species)
 
     ! Solve chemistry at the current time step
     call micm_run(time_step, temperature, pressure, dry_air_density, rate_parameters, &
-                  constituents, errmsg, errcode)
+                  constituents_micm_species, 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)
+    call convert_to_mass_mixing_ratio(dry_air_density, micm_molar_mass_array, constituents_micm_species)
 
+    call update_constituents(micm_indices_constituent_props, constituents_micm_species, &
+                             constituents, errmsg, errcode)
+    if (errcode /= 0) return
+
   end subroutine musica_ccpp_run
 
   !> \section arg_table_musica_ccpp_final Argument Table
   !! \htmlinclude musica_ccpp_final.html
   subroutine musica_ccpp_final(errmsg, errcode)
+    use musica_ccpp_species, only: cleanup_musica_species
     character(len=512), intent(out) :: errmsg
     integer,            intent(out) :: errcode
 
+    call cleanup_musica_species()
     call tuvx_final(errmsg, errcode)
     call micm_final(errmsg, errcode)
 
   end subroutine musica_ccpp_final
 
-end module musica_ccpp
+end module musica_ccpp

schemes/musica/musica_ccpp.meta

@@ -53,7 +53,7 @@
   type = real | kind = kind_phys
   dimensions = (photolysis_wavelength_grid_interface_dimension)
   intent = in
-[ constituent_props ]
+[ constituent_props_ptr ]
   standard_name = ccpp_constituent_properties
   units = None
   type = ccpp_constituent_prop_ptr_t