move NO photolysis logic to TUVx modules

schemes/musica/micm/musica_ccpp_micm.F90

@@ -19,18 +19,16 @@ module musica_ccpp_micm
 
   !> Registers MICM constituent properties with the CCPP
   subroutine micm_register(solver_type, number_of_grid_cells, constituent_props, &
-                           micm_species, includes_no_photolysis, errmsg, errcode)
+                           micm_species, errmsg, errcode)
     use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t
     use musica_ccpp_species,       only: musica_species_t
     use musica_util,               only: error_t
     use iso_c_binding,             only: c_int
-    use musica_ccpp_tuvx_no_photolysis_rate, only : set_NO_photolysis_index
 
     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(:)
-    logical,                                          intent(out) :: includes_no_photolysis
     character(len=512),                               intent(out) :: errmsg
     integer,                                          intent(out) :: errcode
 
@@ -100,18 +98,6 @@ subroutine micm_register(solver_type, number_of_grid_cells, constituent_props, &
     end do
     number_of_rate_parameters = micm%user_defined_reaction_rates%size()
 
-    includes_no_photolysis = .false.
-    ! iterate through the user defined reaction rates and print their names and indices
-    do i = 1, number_of_rate_parameters
-      associate( map => micm%user_defined_reaction_rates )
-        ! check if any rate is named PHOTO.jNO
-        if (index(map%name(i), "PHOTO.jNO") > 0) then
-          includes_no_photolysis = .true.
-          call set_NO_photolysis_index(map%index(i))
-        end if
-      end associate ! map
-    end do
-
   end subroutine micm_register
 
   !> Initializes MICM

schemes/musica/musica_ccpp.F90

@@ -18,7 +18,6 @@ subroutine musica_ccpp_register(constituent_props, errmsg, errcode)
     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
-    use musica_ccpp_tuvx_no_photolysis_rate, only: check_NO_exist
 
     type(ccpp_constituent_properties_t), allocatable, intent(out) :: constituent_props(:)
     character(len=512),                               intent(out) :: errmsg
@@ -28,18 +27,15 @@ subroutine musica_ccpp_register(constituent_props, errmsg, errcode)
     type(ccpp_constituent_properties_t), allocatable :: constituent_props_subset(:)
     type(musica_species_t),              allocatable :: micm_species(:)
     type(musica_species_t),              allocatable :: tuvx_species(:)
-    logical                                          :: includes_no_photolysis
     integer :: number_of_grid_cells
 
-    includes_no_photolysis = .false.
-
     ! Temporary fix until the number of grid cells is only needed to create a MICM state
     ! instead of when the solver is created.
     ! The number of grid cells is not known at this point, so we set it to 1 and recreate
     ! 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, &
-                       micm_species, includes_no_photolysis,errmsg, errcode)
+                       micm_species, errmsg, errcode)
     if (errcode /= 0) return
     constituent_props = constituent_props_subset
     deallocate(constituent_props_subset)
@@ -53,8 +49,6 @@ subroutine musica_ccpp_register(constituent_props, errmsg, errcode)
     call check_tuvx_species_initialization(errmsg, errcode)
     if (errcode /= 0) return
 
-    call check_NO_exist(micm_species, includes_no_photolysis)
-
   end subroutine musica_ccpp_register
 
   !> \section arg_table_musica_ccpp_init Argument Table
@@ -70,8 +64,6 @@ subroutine musica_ccpp_init(horizontal_dimension, vertical_layer_dimension, &
     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
-    use musica_ccpp_tuvx_no_photolysis_rate, &
-      only: is_NO_photolysis_active, set_NO_index_constituent_props, check_NO_initialization
 
     integer,                           intent(in)  :: horizontal_dimension                     ! (count)
     integer,                           intent(in)  :: vertical_layer_dimension                 ! (count)
@@ -84,15 +76,15 @@ subroutine musica_ccpp_init(horizontal_dimension, vertical_layer_dimension, &
     ! local variables
     type(ccpp_constituent_properties_t), allocatable :: constituent_props(:)
     type(musica_species_t),              allocatable :: micm_species(:)
-    logical                                          :: includes_no_photolysis
     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, constituent_props, &
-                       micm_species, includes_no_photolysis, errmsg, errcode)
+                       micm_species, errmsg, errcode)
+    if (errcode /= 0) return
     call micm_init(errmsg, errcode)
     if (errcode /= 0) return
     call tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
@@ -107,13 +99,6 @@ subroutine musica_ccpp_init(horizontal_dimension, vertical_layer_dimension, &
     call check_initialization(errmsg, errcode)
     if (errcode /= 0) return
 
-    if (is_NO_photolysis_active) then
-      call set_NO_index_constituent_props(constituent_props_ptr, errmsg, errcode)
-      if (errcode /= 0) return
-      call check_NO_initialization(errmsg, errcode)
-      if (errcode /= 0) return
-    end if
-
   end subroutine musica_ccpp_init
 
   !> \section arg_table_musica_ccpp_run Argument Table
@@ -137,8 +122,6 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
     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
-    use musica_ccpp_tuvx_no_photolysis_rate, &
-      only: is_NO_photolysis_active, num_NO_photolysis_constituents, NO_photolysis_indices_constituent_props
 
     real(kind_phys),         intent(in)    :: time_step                                         ! s
     real(kind_phys), target, intent(in)    :: temperature(:,:)                                  ! K (column, layer)
@@ -172,24 +155,14 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
     real(kind_phys), dimension(size(constituents, dim=1), &
                                size(constituents, dim=2), &
                                number_of_tuvx_species)    :: constituents_tuvx_species ! kg kg-1
-    real(kind_phys), dimension(size(constituents, dim=1), &
-                               size(constituents, dim=2), &
-                               num_NO_photolysis_constituents) :: constituents_NO_photolysis ! kg kg-1
 
     call extract_subset_constituents(tuvx_indices_constituent_props, constituents, &
                                      constituents_tuvx_species, errmsg, errcode)
     if (errcode /= 0) return
 
-    if (is_NO_photolysis_active) then
-      call extract_subset_constituents(NO_photolysis_indices_constituent_props, constituents, &
-                                      constituents_NO_photolysis, errmsg, errcode)
-      if (errcode /= 0) return
-    end if
-
     ! Calculate photolysis rate constants using TUV-x
     call tuvx_run(temperature, dry_air_density,                 &
                   constituents_tuvx_species,                    &
-                  constituents_NO_photolysis,                   &
                   geopotential_height_wrt_surface_at_midpoint,  &
                   geopotential_height_wrt_surface_at_interface, &
                   surface_geopotential, surface_temperature,    &
@@ -203,6 +176,7 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
                   earth_sun_distance,                           &
                   rate_parameters,                              &
                   errmsg, errcode)
+    if (errcode /= 0) return
 
     call extract_subset_constituents(micm_indices_constituent_props, constituents, &
                                      constituents_micm_species, errmsg, errcode)
@@ -214,6 +188,7 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
     ! Solve chemistry at the current time step
     call micm_run(time_step, temperature, pressure, dry_air_density, rate_parameters, &
                   constituents_micm_species, errmsg, errcode)
+    if (errcode /= 0) return
 
     ! Convert MICM unit back to CAM-SIMA unit (mol m-3  ->  kg kg-1)
     call convert_to_mass_mixing_ratio(dry_air_density, micm_molar_mass_array, constituents_micm_species)

schemes/musica/musica_ccpp_species.F90

@@ -1,5 +1,6 @@
 module musica_ccpp_species
   use ccpp_kinds, only: kind_phys
+  use musica_ccpp_util, only: MUSICA_INT_UNASSIGNED
 
   implicit none
   private
@@ -8,8 +9,6 @@ module musica_ccpp_species
             initialize_molar_mass_array, extract_subset_constituents, update_constituents, &
             check_initialization
 
-  integer, parameter, public :: MUSICA_INT_UNASSIGNED = -99999
-
   !> Definition of musica species object
   type, public :: musica_species_t
     character(len=:), allocatable :: name

schemes/musica/musica_ccpp_util.F90

@@ -11,7 +11,9 @@ module musica_ccpp_util
 
   real(kind_phys), parameter, public :: PI = 3.14159265358979323846_kind_phys
   real(kind_phys), parameter, public :: DEGREE_TO_RADIAN = PI / 180.0_kind_phys
-  real(kind_phys), parameter, public :: EARTH_RADIUS_KM = 6378_kind_phys
+  real(kind_phys), parameter, public :: EARTH_RADIUS_M = 6378000.0_kind_phys
+  real(kind_phys), parameter, public :: AVOGADRO = 6.02214076e23_kind_phys
+  integer,         parameter, public :: MUSICA_INT_UNASSIGNED = -99999
 
 contains
 

schemes/musica/tuvx/musica_ccpp_tuvx.F90

@@ -133,7 +133,7 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
                        errmsg, errcode)
     use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
     use musica_tuvx, only: grid_map_t, profile_map_t, radiator_map_t
-    use musica_util, only: error_t, configuration_t
+    use musica_util, only: error_t, configuration_t, MUSICA_INDEX_MAPPINGS_MAP_ANY
     use musica_ccpp_namelist, only: filename_of_tuvx_micm_mapping_configuration
     use musica_ccpp_util, only: PI
     use musica_ccpp_tuvx_height_grid, &
@@ -155,6 +155,7 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
       only: create_aerosol_optics_radiator, aerosol_optics_label
     use musica_ccpp_tuvx_load_species, &
       only: DRY_AIR_LABEL, O2_LABEL, O3_LABEL, TUVX_GAS_SPECIES_UNITS
+    use musica_ccpp_tuvx_no_photolysis, only: NO_photolysis_init
 
     integer,            intent(in)  :: vertical_layer_dimension      ! (count)
     integer,            intent(in)  :: vertical_interface_dimension  ! (count)
@@ -486,24 +487,30 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
     end if
 
     photolysis_rate_constants_mapping => &
-        index_mappings_t( config, photolysis_rate_constants_ordering, &
+        index_mappings_t( config, MUSICA_INDEX_MAPPINGS_MAP_ANY, &
+                          photolysis_rate_constants_ordering, &
                           micm_rate_parameter_ordering, error )
+    deallocate( photolysis_rate_constants_ordering )
     if (has_error_occurred( error, errmsg, errcode )) then
       deallocate( tuvx )
       tuvx => null()
       call cleanup_tuvx_resources()
-      deallocate( photolysis_rate_constants_ordering )
       return
     end if
 
-    deallocate( photolysis_rate_constants_ordering )
+    call NO_photolysis_init(config, micm_rate_parameter_ordering, errmsg, errcode)
+    if (errcode /= 0) then
+      deallocate( tuvx )
+      tuvx => null()
+      call cleanup_tuvx_resources()
+      return
+    endif
 
   end subroutine tuvx_init
 
   !> Calculates photolysis rate constants for the current model conditions
   subroutine tuvx_run(temperature, dry_air_density,                 &
                       constituents,                                 &
-                      constituents_NO_photolysis,                   &
                       geopotential_height_wrt_surface_at_midpoint,  &
                       geopotential_height_wrt_surface_at_interface, &
                       surface_geopotential, surface_temperature,    &
@@ -520,22 +527,19 @@ subroutine tuvx_run(temperature, dry_air_density,                 &
     use musica_util,                            only: error_t
     use musica_ccpp_tuvx_height_grid,           only: set_height_grid_values, calculate_heights
     use musica_ccpp_tuvx_temperature,           only: set_temperature_values
-    use musica_ccpp_util,                       only: has_error_occurred, PI
+    use musica_ccpp_util,                       only: has_error_occurred, PI, MUSICA_INT_UNASSIGNED
     use musica_ccpp_tuvx_surface_albedo,        only: set_surface_albedo_values
     use musica_ccpp_tuvx_extraterrestrial_flux, only: set_extraterrestrial_flux_values
     use musica_ccpp_tuvx_cloud_optics,          only: set_cloud_optics_values
     use musica_ccpp_tuvx_aerosol_optics,        only: set_aerosol_optics_values
     use musica_ccpp_tuvx_load_species,          only: index_cloud_liquid_water_content, &
                                                       index_dry_air, index_O2, index_O3
     use musica_ccpp_tuvx_gas_species,           only: set_gas_species_values
-    use musica_ccpp_tuvx_no_photolysis_rate,    only: is_NO_photolysis_active, calculate_NO_photolysis_rate, &
-                                                      index_NO_photolysis_rate
-    use musica_ccpp_species,                    only: MUSICA_INT_UNASSIGNED
+    use musica_ccpp_tuvx_no_photolysis,         only: calculate_NO_photolysis_rate_constants
 
     real(kind_phys),    intent(in)    :: temperature(:,:)                                  ! K (column, layer)
     real(kind_phys),    intent(in)    :: dry_air_density(:,:)                              ! kg m-3 (column, layer)
     real(kind_phys),    intent(in)    :: constituents(:,:,:)                               ! kg kg-1 (column, layer, constituent)
-    real(kind_phys),    intent(in)    :: constituents_NO_photolysis(:,:,:)                 ! kg kg-1 (column, layer, constituent)
     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_geopotential(:)                           ! m2 s-2 (column)
@@ -548,14 +552,15 @@ subroutine tuvx_run(temperature, dry_air_density,                 &
     real(kind_phys),    intent(in)    :: air_pressure_thickness(:,:)                       ! Pa (column, layer)
     real(kind_phys),    intent(in)    :: solar_zenith_angle(:)                             ! radians (column)
     real(kind_phys),    intent(in)    :: earth_sun_distance                                ! m
-    real(kind_phys),    intent(inout) :: rate_parameters(:,:,:)                            ! various units (column, layer, reaction)
+    real(kind_phys),    intent(inout) :: rate_parameters(:,:,:)                            ! various units (column, layer, parameter)
     character(len=512), intent(out)   :: errmsg
     integer,            intent(out)   :: errcode
 
     ! local variables
-    real(kind_phys), dimension(size(geopotential_height_wrt_surface_at_midpoint, dim = 2))  :: height_midpoints
-    real(kind_phys), dimension(size(geopotential_height_wrt_surface_at_interface, dim = 2)) :: height_interfaces
-    real(kind_phys), dimension(size(height_interfaces))                                     :: height_deltas ! km
+    real(kind_phys), dimension(size(geopotential_height_wrt_surface_at_midpoint, dim = 2))  :: height_midpoints ! km (layer)
+    real(kind_phys), dimension(size(geopotential_height_wrt_surface_at_interface, dim = 1), &
+                               size(geopotential_height_wrt_surface_at_interface, dim = 2)) :: height_interfaces ! km (column, interface)
+    real(kind_phys), dimension(size(height_interfaces, dim=2))                              :: height_deltas ! km (layer)
     real(kind_phys), dimension(size(rate_parameters, dim=2)+2, &
                                number_of_photolysis_rate_constants) :: photolysis_rate_constants, & ! s-1
                                                                        heating_rates                ! K s-1 (TODO: check units)
@@ -584,9 +589,9 @@ subroutine tuvx_run(temperature, dry_air_density,                 &
                                 geopotential_height_wrt_surface_at_interface(i_col,:), &
                                 surface_geopotential(i_col),                           &
                                 reciprocal_of_gravitational_acceleration,              &
-                                height_midpoints, height_interfaces )
+                                height_midpoints, height_interfaces(i_col,:) )
 
-        call set_height_grid_values( height_grid, height_midpoints, height_interfaces, &
+        call set_height_grid_values( height_grid, height_midpoints, height_interfaces(i_col,:), &
                                      height_deltas, errmsg, errcode )
         if (errcode /= 0) return
 
@@ -634,30 +639,24 @@ subroutine tuvx_run(temperature, dry_air_density,                 &
             max( photolysis_rate_constants(:,:), 0.0_kind_phys )
       end if ! solar zenith angle check
 
-      if (is_NO_photolysis_active) then
-        ! TODO: How do I ensure that the extraterrestrial flux matches the wavelength grid needed for NO photolysis?
-        NO_photolysis_rate_constant = calculate_NO_photolysis_rate(              &
-            solar_zenith_angle(i_col), extraterrestrial_flux, height_interfaces, &
-            dry_air_density(i_col,:), constituents(i_col,:,index_O2),            &
-            constituents(i_col,:,index_O3), constituents_NO_photolysis(i_col,:,:) )
-      end if
-
-      ! TODO: Where does jno go into the photolysis_rate_constants array?
       ! map photolysis rate constants to the host model's rate parameters and vertical grid
       do i_level = 1, size(rate_parameters, dim=2)
         call photolysis_rate_constants_mapping%copy_data( &
             photolysis_rate_constants(size(rate_parameters, dim=2)-i_level+2,:), &
             rate_parameters(i_col,i_level,:) )
-        if (is_NO_photolysis_active) then
-          rate_parameters(i_col,i_level,index_NO_photolysis_rate) = NO_photolysis_rate_constant(i_level)
-        end if
       end do
     end do
 
+    ! TODO: How do I ensure that the extraterrestrial flux matches the wavelength grid needed for NO photolysis?
+    call calculate_NO_photolysis_rate_constants( solar_zenith_angle, extraterrestrial_flux, &
+      height_interfaces, dry_air_density, constituents, rate_parameters )
+
   end subroutine tuvx_run
 
   !> Finalizes TUV-x
   subroutine tuvx_final(errmsg, errcode)
+    use musica_ccpp_tuvx_no_photolysis, only: NO_photolysis_final
+
     character(len=512), intent(out) :: errmsg
     integer,            intent(out) :: errcode
 
@@ -671,6 +670,8 @@ subroutine tuvx_final(errmsg, errcode)
       tuvx => null()
     end if
 
+    call NO_photolysis_final()
+
   end subroutine tuvx_final
 
 end module musica_ccpp_tuvx