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

Solar zenith angle and Earth-Sun distance #171

Draft
wants to merge 10 commits into
base: development
Choose a base branch
from
+100 −41
Draft

Solar zenith angle and Earth-Sun distance #171

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
8063a01
draft sza and earth-sun distance code
mattldawson Nov 26, 2024
073eb1d
finish sza and earth-sun distance changes
mattldawson Nov 27, 2024
da7b2c6
finish tests
mattldawson Nov 27, 2024
372c4ca
Merge branch 'development' into develop-163-solar-zenith-angle
mattldawson Nov 27, 2024
19a0cc2
fix build script
mattldawson Nov 27, 2024
6538d52
formatting and test conditions
mattldawson Nov 27, 2024
47dc398
address reviewer comments
mattldawson Dec 4, 2024
d7220d6
fix Dockerfile
mattldawson Dec 10, 2024
174190e
update to use solar zenith angle and earth-sun distance from host model
mattldawson Dec 11, 2024
8c1a7ce
remove shr_orb_mod
mattldawson Dec 11, 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
10 changes: 8 additions & 2 deletions schemes/musica/musica_ccpp.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,8 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
number_of_photolysis_wavelength_grid_sections, &
photolysis_wavelength_grid_interfaces, extraterrestrial_flux, &
standard_gravitational_acceleration, cloud_area_fraction, &
air_pressure_thickness, errmsg, errcode)
air_pressure_thickness, solar_zenith_angle, &
earth_sun_distance, errmsg, errcode)
use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
use ccpp_kinds, only: kind_phys
use musica_ccpp_micm, only: number_of_rate_parameters
Expand All @@ -101,6 +102,8 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
real(kind_phys), intent(in) :: standard_gravitational_acceleration ! m s-2
real(kind_phys), intent(in) :: cloud_area_fraction(:,:) ! unitless (column, level)
real(kind_phys), intent(in) :: air_pressure_thickness(:,:) ! Pa (column, level)
real(kind_phys), intent(in) :: solar_zenith_angle(:) ! radians (column)
real(kind_phys), intent(in) :: earth_sun_distance ! AU
character(len=512), intent(out) :: errmsg
integer, intent(out) :: errcode

Expand All @@ -121,7 +124,10 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
photolysis_wavelength_grid_interfaces, &
extraterrestrial_flux, &
standard_gravitational_acceleration, &
cloud_area_fraction, constituents, &
cloud_area_fraction, &
solar_zenith_angle, &
earth_sun_distance, &
boulderdaze marked this conversation as resolved.
Show resolved Hide resolved
constituents, &
air_pressure_thickness, rate_parameters, &
errmsg, errcode)

Expand Down
12 changes: 12 additions & 0 deletions schemes/musica/musica_ccpp.meta
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -183,6 +183,18 @@
type = real | kind = kind_phys
dimensions = (horizontal_loop_extent,vertical_layer_dimension)
intent = in
[ solar_zenith_angle ]
standard_name = solar_zenith_angle
units = radians
type = real | kind = kind_phys
dimensions = (horizontal_loop_extent)
intent = in
[ earth_sun_distance ]
standard_name = earth_sun_distance
units = radians
type = real | kind = kind_phys
dimensions = ()
intent = in
[ errmsg ]
standard_name = ccpp_error_message
units = none
Expand Down
5 changes: 5 additions & 0 deletions schemes/musica/musica_ccpp_util.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,16 @@
! SPDX-License-Identifier: Apache-2.0
module musica_ccpp_util

use ccpp_kinds, only: kind_phys

implicit none

private
public :: has_error_occurred

real(kind_phys), parameter, public :: PI = 3.14159265358979323846_kind_phys
Copy link
Member

Choose a reason for hiding this comment

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

I would suggest avoid hardcoding this in a separate location and use the CCPP standard name pi_constant if possible, or if not, maybe include from shr_const_pi (these are the same numerically). I realize this is used in another module and maybe the CCPP variable can be brought in from there.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I guess my preference would be to use the shr_const_pi instead of adding yet another function argument to a list that is getting pretty long, and then having to pass constants down through all the internal functions, but I'll go with the approach the group prefers.

real(kind_phys), parameter, public :: DEGREE_TO_RADIAN = PI / 180.0_kind_phys

contains

!> @brief Evaluate a MUSICA error for failure and convert to CCPP error data
Expand Down
79 changes: 45 additions & 34 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 @@ -13,6 +13,9 @@ module musica_ccpp_tuvx

public :: tuvx_register, tuvx_init, tuvx_run, tuvx_final

real(kind_phys), parameter :: MAX_SOLAR_ZENITH_ANGLE = 110.0_kind_phys ! degrees
real(kind_phys), parameter :: MIN_SOLAR_ZENITH_ANGLE = 0.0_kind_phys ! degrees

type(tuvx_t), pointer :: tuvx => null()
type(grid_t), pointer :: height_grid => null()
type(grid_t), pointer :: wavelength_grid => null()
Expand Down Expand Up @@ -129,6 +132,7 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
use musica_tuvx, only: grid_map_t, profile_map_t, radiator_map_t
use musica_util, only: error_t, configuration_t
use musica_ccpp_namelist, only: filename_of_tuvx_micm_mapping_configuration
use musica_ccpp_util, only: PI
Copy link
Member

Choose a reason for hiding this comment

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

Maybe use CCPP variable pi_constant as an argument, if possible?

use musica_ccpp_tuvx_height_grid, &
only: create_height_grid, height_grid_label, height_grid_unit
use musica_ccpp_tuvx_wavelength_grid, &
Expand Down Expand Up @@ -415,13 +419,15 @@ subroutine tuvx_run(temperature, dry_air_density, &
photolysis_wavelength_grid_interfaces, &
extraterrestrial_flux, &
standard_gravitational_acceleration, &
cloud_area_fraction, constituents, &
cloud_area_fraction, &
solar_zenith_angle, earth_sun_distance, &
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_util, only: has_error_occurred, PI
Copy link
Member

Choose a reason for hiding this comment

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

Similarly, would it be possible to use the CCPP pi_constant here as well?

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
Expand All @@ -438,6 +444,8 @@ subroutine tuvx_run(temperature, dry_air_density, &
real(kind_phys), intent(in) :: extraterrestrial_flux(:) ! photons cm-2 s-1 nm-1
real(kind_phys), intent(in) :: standard_gravitational_acceleration ! m s-2
real(kind_phys), intent(in) :: cloud_area_fraction(:,:) ! unitless (column, layer)
real(kind_phys), intent(in) :: solar_zenith_angle(:) ! radians
real(kind_phys), intent(in) :: earth_sun_distance ! m
real(kind_phys), intent(in) :: constituents(:,:,:) ! various (column, layer, constituent)
real(kind_phys), intent(in) :: air_pressure_thickness(:,:) ! Pa (column, layer)
real(kind_phys), intent(inout) :: rate_parameters(:,:,:) ! various units (column, layer, reaction)
Expand All @@ -451,8 +459,7 @@ subroutine tuvx_run(temperature, dry_air_density, &
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
real(kind_phys) :: solar_zenith_angle_degrees
type(error_t) :: error
integer :: i_col, i_level

Expand All @@ -469,39 +476,43 @@ 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, &

! check if solar zenith angle is within the range to calculate photolysis rate constants
solar_zenith_angle_degrees = solar_zenith_angle(i_col) * 180.0_kind_phys / PI
if (solar_zenith_angle_degrees > MAX_SOLAR_ZENITH_ANGLE .or. &
solar_zenith_angle_degrees < MIN_SOLAR_ZENITH_ANGLE) then
photolysis_rate_constants(:,:) = 0.0_kind_phys
else
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 )
if (errcode /= 0) return
if (errcode /= 0) return

call set_temperature_values( temperature_profile, temperature(i_col,:), &
call set_temperature_values( temperature_profile, temperature(i_col,:), &
surface_temperature(i_col), errmsg, errcode )
if (errcode /= 0) return

call set_cloud_optics_values( cloud_optics, cloud_area_fraction(i_col,:), &
air_pressure_thickness(i_col,:), &
constituents(i_col,:,index_cloud_liquid_water_content), &
reciprocal_of_gravitational_acceleration, &
errmsg, errcode )
if (errcode /= 0) return

! temporary values until these are available from the host model
solar_zenith_angle = 0.0_kind_phys
earth_sun_distance = 1.0_kind_phys

! calculate photolysis rate constants and heating rates
call tuvx%run( solar_zenith_angle, earth_sun_distance, &
photolysis_rate_constants(:,:), heating_rates(:,:), &
error )
if (has_error_occurred( error, errmsg, errcode )) return

! filter out negative photolysis rate constants
photolysis_rate_constants(:,:) = &
max( photolysis_rate_constants(:,:), 0.0_kind_phys )
if (errcode /= 0) return

call set_cloud_optics_values( cloud_optics, cloud_area_fraction(i_col,:), &
air_pressure_thickness(i_col,:), &
constituents(i_col,:,index_cloud_liquid_water_content), &
reciprocal_of_gravitational_acceleration, &
errmsg, errcode )
if (errcode /= 0) return

! calculate photolysis rate constants and heating rates
call tuvx%run( solar_zenith_angle(i_col), earth_sun_distance, &
photolysis_rate_constants(:,:), heating_rates(:,:), &
error )
if (has_error_occurred( error, errmsg, errcode )) return

! filter out negative photolysis rate constants
photolysis_rate_constants(:,:) = &
max( photolysis_rate_constants(:,:), 0.0_kind_phys )
end if ! solar zenith angle check

! map photolysis rate constants to the host model's rate parameters and vertical grid
do i_level = 1, size(rate_parameters, dim=2)
Expand Down
2 changes: 1 addition & 1 deletion test/docker/Dockerfile.musica
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ ENV CCPP_STD_NAMES_PATH="lib/CCPPStandardNames"

RUN cd atmospheric_physics/test \
&& cmake -S . -B build \
-D CMAKE_BUILD_TYPE={BUILD_TYPE} \
-D CMAKE_BUILD_TYPE=${BUILD_TYPE} \
-D CCPP_ENABLE_MUSICA_TESTS=ON \
-D CCPP_ENABLE_MEMCHECK=ON \
&& cmake --build ./build
Expand Down
32 changes: 28 additions & 4 deletions test/musica/test_musica_api.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,12 +1,15 @@
program run_test_musica_ccpp

use ccpp_kinds, only: kind_phys
use musica_ccpp

implicit none

#define ASSERT(x) if (.not.(x)) then; write(*,*) "Assertion failed[", __FILE__, ":", __LINE__, "]: x"; stop 1; endif
#define ASSERT_NEAR( a, b, abs_error ) if( (abs(a - b) >= abs_error) .and. (abs(a - b) /= 0.0) ) then; write(*,*) "Assertion failed[", __FILE__, ":", __LINE__, "]: a, b"; stop 1; endif

real(kind_phys), parameter :: DEGREE_TO_RADIAN = 3.14159265358979323846_kind_phys / 180.0_kind_phys

call test_chapman()
call test_terminator()

Expand Down Expand Up @@ -134,7 +137,6 @@ end subroutine get_wavelength_edges
!> Tests the Chapman chemistry scheme
subroutine test_chapman()
use musica_micm, only: Rosenbrock, RosenbrockStandardOrder
use ccpp_kinds, only: kind_phys
use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t
use musica_ccpp_micm, only: micm
Expand Down Expand Up @@ -176,6 +178,8 @@ subroutine test_chapman()
NUM_SPECIES+NUM_TUVX_CONSTITUENTS) :: constituents ! kg kg-1
real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS, &
NUM_SPECIES+NUM_TUVX_CONSTITUENTS) :: initial_constituents ! kg kg-1
real(kind_phys), dimension(NUM_COLUMNS) :: solar_zenith_angle ! radians
real(kind_phys) :: earth_sun_distance ! AU
type(ccpp_constituent_prop_ptr_t), allocatable :: constituent_props_ptr(:)
type(ccpp_constituent_properties_t), allocatable, target :: constituent_props(:)
type(ccpp_constituent_properties_t), pointer :: const_prop
Expand Down Expand Up @@ -214,6 +218,8 @@ subroutine test_chapman()
cloud_area_fraction(:,2) = (/ 0.3_kind_phys, 0.4_kind_phys /)
air_pressure_thickness(:,1) = (/ 900.0_kind_phys, 905.0_kind_phys /)
air_pressure_thickness(:,2) = (/ 910.0_kind_phys, 915.0_kind_phys /)
solar_zenith_angle = (/ 0.0_kind_phys, 2.1_kind_phys /)
earth_sun_distance = 1.04_kind_phys

filename_of_micm_configuration = 'musica_configurations/chapman/micm/config.json'
filename_of_tuvx_configuration = 'musica_configurations/chapman/tuvx/config.json'
Expand Down Expand Up @@ -317,7 +323,9 @@ subroutine test_chapman()
geopotential_height_wrt_surface_at_interface, surface_geopotential, &
surface_temperature, surface_albedo, num_photolysis_wavelength_grid_sections, &
flux_data_photolysis_wavelength_interfaces, extraterrestrial_flux, &
standard_gravitational_acceleration, cloud_area_fraction, air_pressure_thickness, errmsg, errcode )
standard_gravitational_acceleration, cloud_area_fraction, &
air_pressure_thickness, solar_zenith_angle, earth_sun_distance, errmsg, &
errcode )
if (errcode /= 0) then
write(*,*) trim(errmsg)
stop 3
Expand Down Expand Up @@ -354,6 +362,11 @@ subroutine test_chapman()
ASSERT_NEAR(total_O, total_O_init, 1.0e-13)
end do
end do
do j = 1, NUM_LAYERS
! O and O1D should be lower in the nighttime column
ASSERT(constituents(2,j,O_index) < constituents(1,j,O_index))
ASSERT(constituents(2,j,O1D_index) < constituents(1,j,O1D_index))
end do

deallocate(constituent_props_ptr)

Expand All @@ -362,7 +375,6 @@ end subroutine test_chapman
!> Tests the simple Terminator chemistry scheme
subroutine test_terminator()
use musica_micm, only: Rosenbrock, RosenbrockStandardOrder
use ccpp_kinds, only: kind_phys
use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t
use musica_ccpp_micm, only: micm
Expand Down Expand Up @@ -404,6 +416,8 @@ subroutine test_terminator()
NUM_SPECIES+NUM_TUVX_CONSTITUENTS) :: constituents ! kg kg-1
real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS, &
NUM_SPECIES+NUM_TUVX_CONSTITUENTS) :: initial_constituents ! kg kg-1
real(kind_phys), dimension(NUM_COLUMNS) :: solar_zenith_angle ! radians
real(kind_phys) :: earth_sun_distance ! AU
type(ccpp_constituent_prop_ptr_t), allocatable :: constituent_props_ptr(:)
type(ccpp_constituent_properties_t), allocatable, target :: constituent_props(:)
type(ccpp_constituent_properties_t), pointer :: const_prop
Expand Down Expand Up @@ -442,6 +456,8 @@ subroutine test_terminator()
cloud_area_fraction(:,2) = (/ 0.3_kind_phys, 0.4_kind_phys /)
air_pressure_thickness(:,1) = (/ 900.0_kind_phys, 905.0_kind_phys /)
air_pressure_thickness(:,2) = (/ 910.0_kind_phys, 915.0_kind_phys /)
solar_zenith_angle = (/ 0.0_kind_phys, 2.1_kind_phys /)
earth_sun_distance = 1.04_kind_phys

filename_of_micm_configuration = 'musica_configurations/terminator/micm/config.json'
filename_of_tuvx_configuration = 'musica_configurations/terminator/tuvx/config.json'
Expand Down Expand Up @@ -533,7 +549,9 @@ subroutine test_terminator()
geopotential_height_wrt_surface_at_interface, surface_geopotential, &
surface_temperature, surface_albedo, num_photolysis_wavelength_grid_sections, &
flux_data_photolysis_wavelength_interfaces, extraterrestrial_flux, &
standard_gravitational_acceleration, cloud_area_fraction, air_pressure_thickness, errmsg, errcode )
standard_gravitational_acceleration, cloud_area_fraction, &
air_pressure_thickness, solar_zenith_angle, earth_sun_distance, errmsg, &
errcode )
if (errcode /= 0) then
write(*,*) trim(errmsg)
stop 3
Expand Down Expand Up @@ -564,6 +582,12 @@ subroutine test_terminator()
ASSERT_NEAR(constituents(i,j,NUM_SPECIES+1), initial_constituents(i,j,NUM_SPECIES+1), 1.0e-13)
end do
end do
do j = 1, NUM_LAYERS
! Cl should be lower in the nighttime column
ASSERT(constituents(2,j,Cl_index) < constituents(1,j,Cl_index))
! Cl2 should be higher in the nighttime column
ASSERT(constituents(2,j,Cl2_index) > constituents(1,j,Cl2_index))
end do

deallocate(constituent_props_ptr)

Expand Down
1 change: 1 addition & 0 deletions test/musica/tuvx/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -151,6 +151,7 @@ target_sources(test_tuvx_cloud_optics
${MUSICA_SRC_PATH}/tuvx/musica_ccpp_tuvx_wavelength_grid.F90
${MUSICA_SRC_PATH}/tuvx/musica_ccpp_tuvx_cloud_optics.F90
${MUSICA_SRC_PATH}/musica_ccpp_util.F90
${TO_BE_CCPPIZED_SRC_PATH}/ccpp_tuvx_utils.F90
${CCPP_TEST_SRC_PATH}/ccpp_kinds.F90
)

Expand Down
Loading