Skip to content

Commit

Permalink
add function for no photolysis calculation
Browse files Browse the repository at this point in the history
  • Loading branch information
@boulderdaze
boulderdaze committed Jan 23, 2025
1 parent ecdb406 commit 4d977b4
Show file tree
Hide file tree
Showing 8 changed files with 322 additions and 124 deletions.
42 changes: 38 additions & 4 deletions schemes/musica/musica_ccpp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,10 +14,11 @@ 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
use musica_ccpp_species, only: musica_species_t, register_musica_species
use musica_ccpp_tuvx_load_species, only: check_tuvx_species_initialization
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
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
Expand Down Expand Up @@ -49,6 +50,8 @@ 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)

end subroutine musica_ccpp_register

!> \section arg_table_musica_ccpp_init Argument Table
Expand All @@ -64,6 +67,8 @@ 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_present, set_NO_index_constituent_props, check_NO_initialization

integer, intent(in) :: horizontal_dimension ! (count)
integer, intent(in) :: vertical_layer_dimension ! (count)
Expand Down Expand Up @@ -98,6 +103,13 @@ subroutine musica_ccpp_init(horizontal_dimension, vertical_layer_dimension, &
call check_initialization(errmsg, errcode)
if (errcode /= 0) return

if (is_NO_present) 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
Expand All @@ -121,6 +133,8 @@ 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_present, 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
Expand Down Expand Up @@ -154,14 +168,24 @@ 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_present) 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, &
Expand All @@ -176,9 +200,19 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
rate_parameters, &
errmsg, errcode)

! TODO (jiwon) - update constituents doens't have an effect since tuvx_run takes it as intent in
! is this correct?
call update_constituents(tuvx_indices_constituent_props, constituents_tuvx_species, &
constituents, errmsg, errcode)
if (errcode /= 0) return

if (is_NO_present) then
! TODO(jiwon) - does it need to update NO constituents?
call update_constituents(NO_photolysis_indices_constituent_props, constituents_NO_photolysis, &
constituents, errmsg, errcode)
if (errcode /= 0) return
end if

call extract_subset_constituents(micm_indices_constituent_props, constituents, &
constituents_micm_species, errmsg, errcode)
if (errcode /= 0) return
Expand Down
20 changes: 12 additions & 8 deletions schemes/musica/tuvx/musica_ccpp_tuvx.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -503,6 +503,7 @@ 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, &
Expand All @@ -525,16 +526,16 @@ subroutine tuvx_run(temperature, dry_air_density, &
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, &
index_N2, index_NO, MOLAR_MASS_N2, &
MOLAR_MASS_O2, MOLAR_MASS_O3, MOLAR_MASS_NO
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: calculate_NO_photolysis_rate
use musica_ccpp_species, only: MUSICA_INT_UNASSIGNED
use musica_ccpp_tuvx_no_photolysis_rate, only: is_NO_present

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
Expand Down Expand Up @@ -562,7 +563,7 @@ subroutine tuvx_run(temperature, dry_air_density, &
real(kind_phys) :: solar_zenith_angle_degrees
type(error_t) :: error
integer :: i_col, i_level
real(kind_phys) :: jno(size(constituents, dim=2)) ! s-1
real(kind_phys) :: jno(size(constituents, dim=2)) ! s-1

reciprocal_of_gravitational_acceleration = 1.0_kind_phys / standard_gravitational_acceleration

Expand Down Expand Up @@ -631,11 +632,14 @@ subroutine tuvx_run(temperature, dry_air_density, &
max( photolysis_rate_constants(:,:), 0.0_kind_phys )
end if ! solar zenith angle check

if (index_N2 .ne. MUSICA_INT_UNASSIGNED .and. index_O2 .ne. MUSICA_INT_UNASSIGNED .and. &
index_O3 .ne. MUSICA_INT_UNASSIGNED .and. index_NO .ne. MUSICA_INT_UNASSIGNED) then
if (is_NO_present) then
! TODO: How do I ensure that the extraterrestrial flux matches the wavelength grid needed for NO photolysis?
jno = calculate_NO_photolysis_rate(size(constituents, dim=2), solar_zenith_angle(i_col), extraterrestrial_flux, constituents(i_col,:,:), height_interfaces, &
dry_air_density(i_col,:), index_N2, index_O2, index_O3, index_NO, MOLAR_MASS_N2, MOLAR_MASS_O2, MOLAR_MASS_O3, MOLAR_MASS_NO)
! jno = calculate_NO_photolysis_rate(size(constituents, dim=2), solar_zenith_angle(i_col), extraterrestrial_flux, constituents(i_col,:,:), height_interfaces, &
! dry_air_density(i_col,:), index_N2, index_O2, index_O3, index_NO, MOLAR_MASS_N2, MOLAR_MASS_O2, MOLAR_MASS_O3, MOLAR_MASS_NO)
jno = calculate_NO_photolysis_rate( size(constituents, dim=2), &
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?
Expand Down
68 changes: 32 additions & 36 deletions schemes/musica/tuvx/musica_ccpp_tuvx_load_species.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,8 @@ module musica_ccpp_tuvx_load_species
integer, protected, public :: index_dry_air = MUSICA_INT_UNASSIGNED
integer, protected, public :: index_O2 = MUSICA_INT_UNASSIGNED
integer, protected, public :: index_O3 = MUSICA_INT_UNASSIGNED
integer, protected, public :: index_N2 = MUSICA_INT_UNASSIGNED
integer, protected, public :: index_NO = MUSICA_INT_UNASSIGNED
! integer, protected, public :: index_N2 = MUSICA_INT_UNASSIGNED
! integer, protected, public :: index_NO = MUSICA_INT_UNASSIGNED

! Constants
! Cloud liquid water
Expand All @@ -34,9 +34,6 @@ module musica_ccpp_tuvx_load_species
real(kind_phys), parameter, public :: MOLAR_MASS_DRY_AIR = 0.0289644_kind_phys ! kg mol-1
real(kind_phys), parameter, public :: MOLAR_MASS_O2 = 0.0319988_kind_phys ! kg mol-1
real(kind_phys), parameter, public :: MOLAR_MASS_O3 = 0.0479982_kind_phys ! kg mol-1
real(kind_phys), parameter, public :: MOLAR_MASS_N2 = 0.0280134_kind_phys ! kg mol-1
real(kind_phys), parameter, public :: MOLAR_MASS_NO = 0.0300100_kind_phys ! kg mol-1


contains

Expand Down Expand Up @@ -65,8 +62,8 @@ subroutine configure_tuvx_species(micm_species, tuvx_species, constituent_props,
logical :: is_dry_air_registered = .false.
logical :: is_O2_registered = .false.
logical :: is_O3_registered = .false.
logical :: is_N2_registered = .false.
logical :: is_NO_registered = .false.
! logical :: is_N2_registered = .false.
! logical :: is_NO_registered = .false.
integer :: i_new, i_species, i_tuvx_species

num_micm_species = size(micm_species)
Expand All @@ -92,8 +89,7 @@ subroutine configure_tuvx_species(micm_species, tuvx_species, constituent_props,
! Iterate through the MICM species to check if any TUV-x gas
! species are included; if present, updates the scale height and profiled status.
do i_species = 1, num_micm_species
if (is_dry_air_registered .and. is_O2_registered .and. is_O3_registered .and. &
is_N2_registered .and. is_NO_registered) exit
if ( is_dry_air_registered .and. is_O2_registered .and. is_O3_registered ) exit

if ( micm_species(i_species)%name == DRY_AIR_LABEL ) then
is_dry_air_registered = .true.
Expand All @@ -107,12 +103,12 @@ subroutine configure_tuvx_species(micm_species, tuvx_species, constituent_props,
is_O3_registered = .true.
micm_species(i_species)%profiled = .true.
micm_species(i_species)%scale_height = SCALE_HEIGHT_O3
else if ( micm_species(i_species)%name == 'N2' ) then
num_new_species = num_new_species + 1
is_N2_registered = .true.
else if ( micm_species(i_species)%name == 'NO' ) then
num_new_species = num_new_species + 1
is_NO_registered = .true.
! else if ( micm_species(i_species)%name == 'N2' ) then
! num_new_species = num_new_species + 1
! is_N2_registered = .true.
! else if ( micm_species(i_species)%name == 'NO' ) then
! num_new_species = num_new_species + 1
! is_NO_registered = .true.
end if
end do

Expand Down Expand Up @@ -206,27 +202,27 @@ subroutine configure_tuvx_species(micm_species, tuvx_species, constituent_props,
profiled = .true., &
scale_height = SCALE_HEIGHT_O3 )

if (is_NO_registered) then
i_tuvx_species = i_tuvx_species + 1
index_NO = i_tuvx_species
tuvx_species(i_tuvx_species) = musica_species_t( &
name = 'NO', &
unit = TUVX_GAS_SPECIES_UNITS, &
molar_mass = MOLAR_MASS_NO, & ! kg mol-1
index_musica_species = i_tuvx_species, &
profiled = .true.)
end if

if (is_N2_registered) then
i_tuvx_species = i_tuvx_species + 1
index_N2 = i_tuvx_species
tuvx_species(i_tuvx_species) = musica_species_t( &
name = 'N2', &
unit = TUVX_GAS_SPECIES_UNITS, &
molar_mass = MOLAR_MASS_N2, & ! kg mol-1
index_musica_species = i_tuvx_species, &
profiled = .true.)
end if
! if (is_NO_registered) then
! i_tuvx_species = i_tuvx_species + 1
! index_NO = i_tuvx_species
! tuvx_species(i_tuvx_species) = musica_species_t( &
! name = 'NO', &
! unit = TUVX_GAS_SPECIES_UNITS, &
! molar_mass = MOLAR_MASS_NO, & ! kg mol-1
! index_musica_species = i_tuvx_species, &
! profiled = .true.)
! end if

! if (is_N2_registered) then
! i_tuvx_species = i_tuvx_species + 1
! index_N2 = i_tuvx_species
! tuvx_species(i_tuvx_species) = musica_species_t( &
! name = 'N2', &
! unit = TUVX_GAS_SPECIES_UNITS, &
! molar_mass = MOLAR_MASS_N2, & ! kg mol-1
! index_musica_species = i_tuvx_species, &
! profiled = .true.)
! end if

end subroutine configure_tuvx_species

Expand Down
Loading

0 comments on commit 4d977b4

Please sign in to comment.