Skip to content

Commit

Permalink
checking if there is an no photolysis rate reaction in micm
Browse files Browse the repository at this point in the history
  • Loading branch information
@K20shores
K20shores committed Feb 5, 2025
1 parent 5399a16 commit 8cb899d
Show file tree
Hide file tree
Showing 5 changed files with 47 additions and 69 deletions.
15 changes: 14 additions & 1 deletion schemes/musica/micm/musica_ccpp_micm.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,16 +19,18 @@ module musica_ccpp_micm

!> Registers MICM constituent properties with the CCPP
subroutine micm_register(solver_type, number_of_grid_cells, constituent_props, &
micm_species, errmsg, errcode)
micm_species, includes_no_photolysis, 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

Expand Down Expand Up @@ -98,6 +100,17 @@ subroutine micm_register(solver_type, number_of_grid_cells, constituent_props, &
end do
number_of_rate_parameters = micm%user_defined_reaction_rates%size()

! 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
Expand Down
15 changes: 12 additions & 3 deletions schemes/musica/musica_ccpp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ 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

! Temporary fix until the number of grid cells is only needed to create a MICM state
Expand All @@ -36,10 +37,17 @@ subroutine musica_ccpp_register(constituent_props, errmsg, errcode)
! 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, errmsg, errcode)
micm_species, includes_no_photolysis,errmsg, errcode)
if (errcode /= 0) return
constituent_props = constituent_props_subset
deallocate(constituent_props_subset)

if (includes_no_photolysis) then
print *, "Do NO photolysis? yes"
else
print *, "Do NO photolysis? no"
end if


call tuvx_register(micm_species, tuvx_species, constituent_props_subset, &
errmsg, errcode)
Expand All @@ -50,7 +58,7 @@ 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)
call check_NO_exist(micm_species, includes_no_photolysis)

end subroutine musica_ccpp_register

Expand Down Expand Up @@ -81,14 +89,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, errmsg, errcode)
micm_species, includes_no_photolysis, errmsg, errcode)
call micm_init(errmsg, errcode)
if (errcode /= 0) return
call tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
Expand Down
6 changes: 5 additions & 1 deletion schemes/musica/tuvx/musica_ccpp_tuvx.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -528,7 +528,8 @@ subroutine tuvx_run(temperature, dry_air_density, &
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
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

real(kind_phys), intent(in) :: temperature(:,:) ! K (column, layer)
Expand Down Expand Up @@ -647,6 +648,9 @@ subroutine tuvx_run(temperature, dry_air_density, &
photolysis_rate_constants(size(rate_parameters, dim=2)-i_level+2,:), &
rate_parameters(i_col,i_level,:) )
end do
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 subroutine tuvx_run
Expand Down
15 changes: 12 additions & 3 deletions schemes/musica/tuvx/musica_ccpp_tuvx_no_photolysis_rate.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,12 +23,13 @@ module musica_ccpp_tuvx_no_photolysis_rate

private
public :: check_NO_exist, set_NO_index_constituent_props, check_NO_initialization, &
calculate_NO_photolysis_rate
calculate_NO_photolysis_rate, set_NO_photolysis_index

character(len=*), parameter, public :: N2_label = 'N2'
character(len=*), parameter, public :: NO_label = 'NO'
integer, parameter, public :: num_NO_photolysis_constituents = 2
logical, protected, public :: is_NO_photolysis_active = .false.
integer, protected, public :: index_NO_photolysis_rate = MUSICA_INT_UNASSIGNED
integer, protected, public :: index_N2 = MUSICA_INT_UNASSIGNED
integer, protected, public :: index_NO = MUSICA_INT_UNASSIGNED
integer, protected, public :: constituent_props_index_N2 = MUSICA_INT_UNASSIGNED
Expand All @@ -42,10 +43,11 @@ module musica_ccpp_tuvx_no_photolysis_rate

!> Checks for the presence of NO and N2, and sets the flag to indicate
!! whether the NO photolysis calculation is enabled or disabled
subroutine check_NO_exist(micm_species)
subroutine check_NO_exist(micm_species, includes_no_photolysis)
use musica_ccpp_species, only: musica_species_t

type(musica_species_t), intent(in) :: micm_species(:)
logical, intent(in) :: includes_no_photolysis

! local variables
logical :: is_N2_registered = .false.
Expand All @@ -72,12 +74,19 @@ subroutine check_NO_exist(micm_species)
end if
end do

if (is_N2_registered .and. is_NO_registered) then
if (is_N2_registered .and. is_NO_registered .and. includes_no_photolysis) then
is_NO_photolysis_active = .true.
end if

end subroutine check_NO_exist

subroutine set_NO_photolysis_index(no_index)
integer, intent(in) :: no_index

index_NO_photolysis_rate = no_index

end subroutine set_NO_photolysis_index

!> Gets the indices of NO constituents from the CCPP constituent properties,
!! and stores them in the indices array along with their relative positions
!! in the NO constituents array
Expand Down
65 changes: 4 additions & 61 deletions test/musica/tuvx/test_tuvx_no_photolysis.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ subroutine test_NO_and_N2_absence()
index_musica_species = i_species )
end do

call check_NO_exist( micm_species )
call check_NO_exist( micm_species, .true. )
ASSERT(.not. is_NO_photolysis_active)

end subroutine test_NO_and_N2_absence
Expand Down Expand Up @@ -67,7 +67,7 @@ subroutine test_only_NO_exist()
index_musica_species = i_species )
end do

call check_NO_exist( micm_species )
call check_NO_exist( micm_species, .true. )
ASSERT(.not. is_NO_photolysis_active)

end subroutine test_only_NO_exist
Expand Down Expand Up @@ -95,7 +95,7 @@ subroutine test_only_N2_exist()
index_musica_species = i_species )
end do

call check_NO_exist( micm_species )
call check_NO_exist( micm_species, .true. )
ASSERT(.not. is_NO_photolysis_active)

end subroutine test_only_N2_exist
Expand Down Expand Up @@ -123,7 +123,7 @@ subroutine test_NO_and_N2_exist()
index_musica_species = i_species )
end do

call check_NO_exist( micm_species )
call check_NO_exist( micm_species, .true. )
ASSERT(is_NO_photolysis_active)
ASSERT(MOLAR_MASS_N2 == 0.0280134_kind_phys)
ASSERT(MOLAR_MASS_NO == 0.0300100_kind_phys)
Expand Down Expand Up @@ -340,61 +340,4 @@ subroutine test_calculate_NO_photolysis_rate()

end subroutine test_calculate_NO_photolysis_rate

! subroutine test_calculate_NO_photolysis_rate()
! use ccpp_kinds, only: kind_phys
! implicit none

! integer, parameter :: number_of_vertical_layers = 2

! real(kind_phys) :: solar_zenith_angle
! real(kind_phys), dimension(4) :: extraterrestrial_flux
! real(kind_phys), dimension(1,number_of_vertical_layers,4) :: constituents ! (column, layers, constituents) kg kg-1
! real(kind_phys), dimension(1,number_of_vertical_layers,2) :: constituents_NO_photolysis
! real(kind_phys), dimension(3) :: height_at_interfaces ! (layers + 1) km
! real(kind_phys), dimension(1,number_of_vertical_layers) :: dry_air_density ! (column, layers) kg m-3
! integer :: N2_index, O2_index, O3_index, NO_index
! real(kind_phys) :: molar_mass_N2, molar_mass_O2, molar_mass_O3, molar_mass_NO ! kg mol-1
! real(kind_phys), dimension(number_of_vertical_layers) :: jNO
! ! Some of this initialization data corresponds to values listed in the Minschwaner and Siskind (1993) paper
! ! see the musica_ccpp_tuvx_no_photolysis_rate.F90 file for the full citation

! ! Initialize test data
! solar_zenith_angle = 60.0_kind_phys
! extraterrestrial_flux = (/ 1.5e13_kind_phys, 1.4e13_kind_phys, 1.3e13_kind_phys, 1.2e13_kind_phys /)
! height_at_interfaces = reshape([40.0_kind_phys, 30.0_kind_phys, 20.0_kind_phys], shape(height_at_interfaces))
! dry_air_density = reshape([1.2_kind_phys, 1.1_kind_phys], shape(dry_air_density))
! N2_index = 1
! O2_index = 2
! O3_index = 3
! NO_index = 4
! molar_mass_N2 = 0.0280134_kind_phys
! molar_mass_O2 = 0.0319988_kind_phys
! molar_mass_O3 = 0.0479982_kind_phys
! molar_mass_NO = 0.3000610_kind_phys

! constituents(1,:,N2_index) = 0.78_kind_phys
! constituents(1,:,O2_index) = 0.21_kind_phys
! constituents(1,:,O3_index) = 20.0e-9_kind_phys
! constituents(1,:,NO_index) = 0.3e-9_kind_phys
! constituents_NO_photolysis(1,:,1) = constituents(1,:,N2_index)
! constituents_NO_photolysis(1,:,2) = constituents(1,:,NO_index)

! ! Call the function to test
! ! jNO = calculate_NO_photolysis_rate(size(constituents, dim=2), solar_zenith_angle, extraterrestrial_flux, constituents(1,:,:), height_at_interfaces, &
! ! dry_air_density(1,:), N2_index, O2_index, O3_index, NO_index, molar_mass_N2, molar_mass_O2, molar_mass_O3, molar_mass_NO)

! jNO = calculate_NO_photolysis_rate(size(constituents, dim=2), &
! solar_zenith_angle, extraterrestrial_flux, height_at_interfaces, &
! dry_air_density(1,:), constituents(1,:,O2_index), &
! constituents(1,:,O3_index), constituents_NO_photolysis(1,:,:))

! ! Validate the results
! print *, jNO
! ASSERT(jNO(1) .ne. 0.0_kind_phys)
! ASSERT(jNO(2) .ne. 0.0_kind_phys)
! ASSERT(jNO(1) .lt. 1.0_kind_phys)
! ASSERT(jNO(2) .lt. 1.0_kind_phys)

! end subroutine test_calculate_NO_photolysis_rate

end program test_tuvx_no_photolysis

0 comments on commit 8cb899d

Please sign in to comment.