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

Set gas-species profiles in TUV-x prior to calculating rate constants #174

Closed
wants to merge 16 commits into from
Closed

Set gas-species profiles in TUV-x prior to calculating rate constants #174

Show file tree
Hide file tree
Changes from 14 commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
8072cd9
draft cloud optics calculations
mattldawson Nov 22, 2024
18723a1
Merge branch 'development' into develop-100-cloud-optics
mattldawson Nov 22, 2024
2a9d9ab
add valgrind suppressions
mattldawson Nov 23, 2024
645ec88
fix typo
mattldawson Nov 23, 2024
4216e13
fix formatting
mattldawson Nov 23, 2024
4c11684
address reviewer comments
mattldawson Nov 25, 2024
2f8afe6
address reviewer comments
mattldawson Nov 25, 2024
4454a16
add species
boulderdaze Dec 2, 2024
dbf019c
stash for profile group
boulderdaze Dec 4, 2024
0ed4995
fix bugs - test passed
boulderdaze Dec 5, 2024
b24f5ee
test passed
boulderdaze Dec 6, 2024
ce5190f
update cam sima data commit
boulderdaze Dec 6, 2024
b47819d
merge dev
boulderdaze Dec 6, 2024
8afeebd
add a unit test for gas species profile
boulderdaze Dec 6, 2024
274f446
minor fix
boulderdaze Dec 6, 2024
3f55081
git stash
boulderdaze Dec 14, 2024
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
3 changes: 2 additions & 1 deletion .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,4 +4,5 @@
**/CMakeCache.txt
**/CMakeFiles/
.vscode
xcode/
xcode/
*.mod
29 changes: 15 additions & 14 deletions schemes/musica/musica_ccpp.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,20 +111,7 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
number_of_rate_parameters) :: rate_parameters ! various units
integer :: i_elem

! Calculate photolysis rate constants using TUV-x
call tuvx_run(temperature, dry_air_density, &
geopotential_height_wrt_surface_at_midpoint, &
geopotential_height_wrt_surface_at_interface, &
surface_geopotential, surface_temperature, &
surface_albedo, &
number_of_photolysis_wavelength_grid_sections, &
photolysis_wavelength_grid_interfaces, &
extraterrestrial_flux, &
standard_gravitational_acceleration, &
cloud_area_fraction, constituents, &
air_pressure_thickness, rate_parameters, &
errmsg, errcode)

! TODO(jiwon) this might not be correct because it doesn't know the index
Copy link
Collaborator Author

@boulderdaze boulderdaze Dec 6, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @peverwhee, is this code still valid (line 116 - 122)? We discussed that the sequence of the constituents array is not guaranteed. Do we now need to search for the index of the MICM constituents before attempting to retrieve the molar mass?

    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

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

hi @boulderdaze sorry for the delay!

Yes, you're correct that the order is not guaranteed. So you'll have to grab the index by the standard name. Haipeng created ccpp_const_get_idx in https://github.com/ESCOMP/atmospheric_physics/blob/main/to_be_ccppized/ccpp_const_utils.F90 for this purpose.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great, thank you for confirming!

! 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)
Expand All @@ -147,6 +134,20 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
! 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)

! Calculate photolysis rate constants using TUV-x
call tuvx_run(temperature, dry_air_density, &
geopotential_height_wrt_surface_at_midpoint, &
geopotential_height_wrt_surface_at_interface, &
surface_geopotential, surface_temperature, &
surface_albedo, &
number_of_photolysis_wavelength_grid_sections, &
photolysis_wavelength_grid_interfaces, &
extraterrestrial_flux, &
standard_gravitational_acceleration, &
cloud_area_fraction, constituents, &
air_pressure_thickness, rate_parameters, &
errmsg, errcode)

! Solve chemistry at the current time step
call micm_run(time_step, temperature, pressure, dry_air_density, rate_parameters, &
constituents, errmsg, errcode)
Expand Down
112 changes: 96 additions & 16 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 @@ -7,6 +7,7 @@ module musica_ccpp_tuvx
use musica_ccpp_util, only: has_error_occurred
use musica_tuvx, only: tuvx_t, grid_t, profile_t, radiator_t
use musica_util, only: mappings_t, index_mappings_t
use musica_ccpp_tuvx_gas_species_profile, only: gas_species_t, profile_group_t

implicit none
private
Expand All @@ -21,6 +22,10 @@ module musica_ccpp_tuvx
type(profile_t), pointer :: extraterrestrial_flux_profile => null()
type(radiator_t), pointer :: cloud_optics => null()
type(index_mappings_t), pointer :: photolysis_rate_constants_mapping => null( )

type(gas_species_t), allocatable :: gas_species_group(:)
type(profile_group_t), allocatable :: profile_gas_species_group(:)

integer, parameter :: DEFAULT_NUM_PHOTOLYSIS_RATE_CONSTANTS = 0
integer :: number_of_photolysis_rate_constants = DEFAULT_NUM_PHOTOLYSIS_RATE_CONSTANTS
integer, parameter :: DEFAULT_INDEX_NOT_FOUND = -1
Expand Down Expand Up @@ -48,8 +53,13 @@ subroutine reset_tuvx_map_state( grids, profiles, radiators )

end subroutine reset_tuvx_map_state

!> This is a helper subroutine created to deallocate objects associated with TUV-x
!> Deallocates objects associated with TUV-x
subroutine cleanup_tuvx_resources()
use musica_ccpp_tuvx_gas_species_profile, &
only: deallocate_gas_species_profile_group

! local variable
integer :: i_species

if (associated( height_grid )) then
deallocate( height_grid )
Expand Down Expand Up @@ -86,6 +96,8 @@ subroutine cleanup_tuvx_resources()
photolysis_rate_constants_mapping => null()
end if

call deallocate_gas_species_profile_group( profile_gas_species_group )

end subroutine cleanup_tuvx_resources

!> Registers constituent properties with the CCPP needed by TUV-x
Expand Down Expand Up @@ -142,6 +154,8 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
extraterrestrial_flux_unit
use musica_ccpp_tuvx_cloud_optics, &
only: create_cloud_optics_radiator, cloud_optics_label
use musica_ccpp_tuvx_gas_species_profile, &
only: configure_gas_species, create_gas_species_profile_group

integer, intent(in) :: vertical_layer_dimension ! (count)
integer, intent(in) :: vertical_interface_dimension ! (count)
Expand All @@ -152,12 +166,14 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
integer, intent(out) :: errcode

! local variables
type(grid_map_t), pointer :: grids
type(profile_map_t), pointer :: profiles
type(radiator_map_t), pointer :: radiators
type(configuration_t) :: config
type(mappings_t), pointer :: photolysis_rate_constants_ordering
type(error_t) :: error
type(error_t) :: error
type(grid_map_t), pointer :: grids
type(profile_map_t), pointer :: profiles
type(radiator_map_t), pointer :: radiators
type(configuration_t) :: config
type(mappings_t), pointer :: photolysis_rate_constants_ordering
integer :: number_of_gas_species
integer :: i_species

! Get needed indices in constituents array
call ccpp_const_get_idx(constituent_props, CLOUD_LIQUID_WATER_CONTENT_LABEL, &
Expand Down Expand Up @@ -252,6 +268,34 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
return
end if

call configure_gas_species( constituent_props, gas_species_group, &
errmsg, errcode )
if (errcode /= 0) then
call reset_tuvx_map_state( grids, profiles, null() )
call cleanup_tuvx_resources()
return
end if

number_of_gas_species = size(gas_species_group)
allocate( profile_gas_species_group( number_of_gas_species ) )

call create_gas_species_profile_group( height_grid, &
gas_species_group, profile_gas_species_group, errmsg, errcode)
if (errcode /= 0) then
call reset_tuvx_map_state( grids, profiles, null() )
call cleanup_tuvx_resources()
return
endif

do i_species = 1, number_of_gas_species
call profiles%add( profile_gas_species_group(i_species)%profile, error )
if (has_error_occurred( error, errmsg, errcode )) then
call reset_tuvx_map_state( grids, profiles, null() )
call cleanup_tuvx_resources()
return
end if
end do

radiators => radiator_map_t( error )
if (has_error_occurred( error, errmsg, errcode )) then
call reset_tuvx_map_state( grids, profiles, null() )
Expand Down Expand Up @@ -350,6 +394,22 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
return
end if

if (.not. allocated( profile_gas_species_group)) then
allocate( profile_gas_species_group( number_of_gas_species ) )
end if

do i_species = 1, number_of_gas_species
profile_gas_species_group(i_species)%profile => &
profiles%get( gas_species_group(i_species)%label, gas_species_group(i_species)%unit, error )
if (has_error_occurred( error, errmsg, errcode )) then
deallocate( tuvx )
tuvx => null()
call reset_tuvx_map_state( grids, profiles, null() )
call cleanup_tuvx_resources()
return
end if
end do

radiators => tuvx%get_radiators( error )
if (has_error_occurred( error, errmsg, errcode )) then
deallocate( tuvx )
Expand Down Expand Up @@ -418,13 +478,14 @@ subroutine tuvx_run(temperature, dry_air_density, &
cloud_area_fraction, constituents, &
air_pressure_thickness, rate_parameters, &
errmsg, errcode)
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
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_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
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_gas_species_profile, only: set_gas_species_values, MOLAR_MASS_DRY_AIR

real(kind_phys), intent(in) :: temperature(:,:) ! K (column, layer)
real(kind_phys), intent(in) :: dry_air_density(:,:) ! kg m-3 (column, layer)
Expand All @@ -447,14 +508,16 @@ subroutine tuvx_run(temperature, dry_air_density, &
! 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(constituents, dim = 2)) :: gas_species_constituents
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)
real(kind_phys) :: reciprocal_of_gravitational_acceleration ! s2 m-1
real(kind_phys) :: solar_zenith_angle ! degrees
real(kind_phys) :: earth_sun_distance ! AU
type(error_t) :: error
integer :: i_col, i_level
integer :: i_col, i_level, i_gas_species

reciprocal_of_gravitational_acceleration = 1.0_kind_phys / standard_gravitational_acceleration

Expand All @@ -469,13 +532,14 @@ subroutine tuvx_run(temperature, dry_air_density, &
if (errcode /= 0) return

do i_col = 1, size(temperature, dim=1)

call calculate_heights( geopotential_height_wrt_surface_at_midpoint(i_col,:), &
geopotential_height_wrt_surface_at_interface(i_col,:), &
surface_geopotential(i_col), &
reciprocal_of_gravitational_acceleration, &
height_midpoints, height_interfaces )
call set_height_grid_values( height_grid, height_midpoints, height_interfaces, &
errmsg, errcode )
height_deltas, errmsg, errcode )
if (errcode /= 0) return

call set_temperature_values( temperature_profile, temperature(i_col,:), &
Expand All @@ -489,6 +553,19 @@ subroutine tuvx_run(temperature, dry_air_density, &
errmsg, errcode )
if (errcode /= 0) return

do i_gas_species = 1, size(gas_species_group)
if (gas_species_group(i_gas_species)%label == "air") then
gas_species_constituents = dry_air_density(i_col,:) * MOLAR_MASS_DRY_AIR
else
gas_species_constituents = constituents(i_col,:,gas_species_group(i_gas_species)%index_constituent_props)
end if
call set_gas_species_values(profile_gas_species_group(i_gas_species)%profile, &
gas_species_group(i_gas_species), &
gas_species_constituents, & ! mol m-3
height_deltas, errmsg, errcode)
if (errcode /= 0) return
end do

! temporary values until these are available from the host model
solar_zenith_angle = 0.0_kind_phys
earth_sun_distance = 1.0_kind_phys
Expand All @@ -515,12 +592,15 @@ end subroutine tuvx_run

!> Finalizes TUV-x
subroutine tuvx_final(errmsg, errcode)
use musica_ccpp_tuvx_gas_species_profile, only: deallocate_gas_species

character(len=512), intent(out) :: errmsg
integer, intent(out) :: errcode

errmsg = ''
errcode = 0

call deallocate_gas_species( gas_species_group )
call cleanup_tuvx_resources()

if (associated( tuvx )) then
Expand Down
Loading
Loading