Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add function for no photolysis calculation #191

Draft
wants to merge 1 commit into
base: 101-calculate-no-photolysis-rate-constants
Choose a base branch
from
+322 −124
Draft

Add function for no photolysis calculation #191

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 38 additions & 4 deletions schemes/musica/musica_ccpp.F90
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading