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

Merged
merged 23 commits into from
Dec 14, 2024
Merged

Solar zenith angle and Earth-Sun distance #171

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
23 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
2a87a0c
update musica namelist parameters
mattldawson Dec 11, 2024
cbc590a
recreate micm solver at init
mattldawson Dec 11, 2024
f797cc1
update horizontal loop argument
mattldawson Dec 12, 2024
39da138
update meta data
mattldawson Dec 12, 2024
e3c190e
create musica dependencies module
mattldawson Dec 12, 2024
5122681
modfy musica dependencies module to be a scheme
mattldawson Dec 12, 2024
89587af
fix bug in et flux interpolation
mattldawson Dec 12, 2024
50e0e90
move musica dependencies to cam-sima
mattldawson Dec 12, 2024
b920c7e
update metadata
mattldawson Dec 12, 2024
5582de2
fix meta data
mattldawson Dec 13, 2024
f82d4c0
address reviewer comments
mattldawson Dec 13, 2024
3aa9953
address reviewer comments
mattldawson Dec 13, 2024
5b7c25c
use photons in units
mattldawson Dec 13, 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
25 changes: 23 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 @@ -75,13 +75,16 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
geopotential_height_wrt_surface_at_interface, surface_geopotential, &
surface_temperature, surface_albedo, &
number_of_photolysis_wavelength_grid_sections, &
photolysis_wavelength_grid_interfaces, extraterrestrial_flux, &
photolysis_wavelength_grid_interfaces, extraterrestrial_flux, &
standard_gravitational_acceleration, cloud_area_fraction, &
air_pressure_thickness, errmsg, errcode)
air_pressure_thickness, latitude, longitude, &
earth_eccentricity, earth_obliquity, perihelion_longitude, &
moving_vernal_equinox_longitude, calendar_day, 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
use musica_ccpp_micm_util, only: convert_to_mol_per_cubic_meter, convert_to_mass_mixing_ratio
use musica_ccpp_util, only: calculate_solar_zenith_angle_and_earth_sun_distance

real(kind_phys), intent(in) :: time_step ! s
real(kind_phys), target, intent(in) :: temperature(:,:) ! K
Expand All @@ -101,6 +104,14 @@ 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) :: latitude(:) ! radians (column)
real(kind_phys), intent(in) :: longitude(:) ! radians (column)

real(kind_phys), intent(in) :: earth_eccentricity ! Earth's eccentricity factor (unitless) (typically 0 to 0.1)
real(kind_phys), intent(in) :: earth_obliquity ! Earth's obliquity in radians
real(kind_phys), intent(in) :: perihelion_longitude ! Earth's mean perihelion longitude at the vernal equinox (radians)
real(kind_phys), intent(in) :: moving_vernal_equinox_longitude ! Earth's moving vernal equinox longitude of perihelion plus pi (radians)
real(kind_phys), intent(in) :: calendar_day ! fractional calendar day
character(len=512), intent(out) :: errmsg
integer, intent(out) :: errcode

Expand All @@ -109,8 +120,16 @@ 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_rate_parameters) :: rate_parameters ! various units
real(kind_phys), dimension(size(latitude)) :: solar_zenith_angle ! radians
real(kind_phys) :: earth_sun_distance ! AU
integer :: i_elem

! Calculate solar zenith angle and Earth-Sun distance
call calculate_solar_zenith_angle_and_earth_sun_distance(calendar_day, &
latitude, longitude, earth_eccentricity, earth_obliquity, perihelion_longitude, &
moving_vernal_equinox_longitude, solar_zenith_angle, earth_sun_distance, &
errmsg, errcode)

! Calculate photolysis rate constants using TUV-x
call tuvx_run(temperature, dry_air_density, &
geopotential_height_wrt_surface_at_midpoint, &
Expand All @@ -121,6 +140,8 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
photolysis_wavelength_grid_interfaces, &
extraterrestrial_flux, &
standard_gravitational_acceleration, &
solar_zenith_angle, &
earth_sun_distance, &
cloud_area_fraction, constituents, &
air_pressure_thickness, rate_parameters, &
errmsg, errcode)
Expand Down
42 changes: 42 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,48 @@
type = real | kind = kind_phys
dimensions = (horizontal_loop_extent,vertical_layer_dimension)
intent = in
[ latitude ]
standard_name = latitude
units = radians
type = real | kind = kind_phys
dimensions = (horizontal_loop_extent)
intent = in
[ longitude ]
standard_name = longitude
units = radians
type = real | kind = kind_phys
dimensions = (horizontal_loop_extent)
intent = in
[ earth_eccentricity ]
standard_name = eccentricity_factor
units = 1
type = real | kind = kind_phys
dimensions = ()
intent = in
[ earth_obliquity ]
standard_name = obliquity
units = radians
type = real | kind = kind_phys
dimensions = ()
intent = in
[ perihelion_longitude ]
standard_name = mean_longitude_of_perihelion_at_vernal_equinox
units = radians
type = real | kind = kind_phys
dimensions = ()
intent = in
[ moving_vernal_equinox_longitude ]
standard_name = moving_vernal_equinox_longitude_perihelion_plus_pi
units = radians
type = real | kind = kind_phys
dimensions = ()
intent = in
[ calendar_day ]
standard_name = fractional_calendar_days_on_end_of_current_timestep
units = 1
type = real | kind = kind_phys
dimensions = ()
intent = in
[ errmsg ]
standard_name = ccpp_error_message
units = none
Expand Down
56 changes: 55 additions & 1 deletion schemes/musica/musica_ccpp_util.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,10 +2,15 @@
! SPDX-License-Identifier: Apache-2.0
module musica_ccpp_util

use ccpp_kinds, only: kind_phys

implicit none

private
public :: has_error_occurred
public :: has_error_occurred, calculate_solar_zenith_angle_and_earth_sun_distance

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

Expand Down Expand Up @@ -37,4 +42,53 @@ logical function has_error_occurred(error, error_message, error_code)

end function has_error_occurred

!> Calculate the solar zenith angle and Earth-Sun distance
!> @param[in] calendar_day Calendar day, including fraction
!> @param[in] latitude Latitude in radians
!> @param[in] longitude Longitude in radians
!> @param[in] earth_eccentricity Earth's eccentricity factor (unitless)
!> @param[in] earth_obliquity Earth's obliquity in radians
!> @param[in] perihelion_longitude Earth's mean perihelion longitude at the vernal equinox (radians)
!> @param[in] moving_vernal_equinox_longitude Earth's moving vernal equinox longitude of perihelion plus pi (radians)
!> @param[out] solar_zenith_angle Solar zenith angle in radians
!> @param[out] earth_sun_distance Earth-Sun distance in AU
!> @param[out] errmsg Error message
!> @param[out] errcode Error code
subroutine calculate_solar_zenith_angle_and_earth_sun_distance(calendar_day, &
latitude, longitude, earth_eccentricity, earth_obliquity, perihelion_longitude, &
moving_vernal_equinox_longitude, solar_zenith_angle, earth_sun_distance, &
errmsg, errcode)
use shr_orb_mod, only: shr_orb_decl, shr_orb_cosz
use musica_util, only: error_t

real(kind_phys), intent(in) :: calendar_day
real(kind_phys), intent(in) :: latitude(:)
real(kind_phys), intent(in) :: longitude(:)
real(kind_phys), intent(in) :: earth_eccentricity
real(kind_phys), intent(in) :: earth_obliquity
real(kind_phys), intent(in) :: perihelion_longitude
real(kind_phys), intent(in) :: moving_vernal_equinox_longitude
real(kind_phys), intent(out) :: solar_zenith_angle(:)
real(kind_phys), intent(out) :: earth_sun_distance
character(len=512), intent(out) :: errmsg
integer, intent(out) :: errcode

real(kind_phys) :: delta
integer :: i_sza

errcode = 0
errmsg = ''

! Calculate earth/orbit parameters
call shr_orb_decl(calendar_day, earth_eccentricity, earth_obliquity, &
perihelion_longitude, moving_vernal_equinox_longitude, &
delta, earth_sun_distance)

! Calculate solar zenith angle
do i_sza = 1, size(solar_zenith_angle)
solar_zenith_angle(i_sza) = acos(shr_orb_cosz(calendar_day, latitude(i_sza), longitude(i_sza), delta))
end do

end subroutine calculate_solar_zenith_angle_and_earth_sun_distance

end module musica_ccpp_util
58 changes: 34 additions & 24 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,14 @@ subroutine tuvx_run(temperature, dry_air_density, &
photolysis_wavelength_grid_interfaces, &
extraterrestrial_flux, &
standard_gravitational_acceleration, &
solar_zenith_angle, earth_sun_distance, &
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_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 @@ -437,6 +442,8 @@ subroutine tuvx_run(temperature, dry_air_density, &
real(kind_phys), intent(in) :: photolysis_wavelength_grid_interfaces(:) ! nm
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) :: solar_zenith_angle(:) ! radians
real(kind_phys), intent(in) :: earth_sun_distance ! m
real(kind_phys), intent(in) :: cloud_area_fraction(:,:) ! unitless (column, layer)
real(kind_phys), intent(in) :: constituents(:,:,:) ! various (column, layer, constituent)
real(kind_phys), intent(in) :: air_pressure_thickness(:,:) ! Pa (column, layer)
Expand All @@ -451,8 +458,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,18 +475,25 @@ 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
if (errcode /= 0) return

call set_cloud_optics_values( cloud_optics, cloud_area_fraction(i_col,:), &
air_pressure_thickness(i_col,:), &
Expand All @@ -489,19 +502,16 @@ subroutine tuvx_run(temperature, dry_air_density, &
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
! 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 )
! 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
27 changes: 27 additions & 0 deletions test/musica/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,32 @@ set(MUSICA_ENABLE_INSTALL OFF)

FetchContent_MakeAvailable(musica)

# MUSICA utilities
add_executable(test_musica_ccpp_util test_musica_ccpp_util.F90)

target_sources(test_musica_ccpp_util
PUBLIC
${MUSICA_SRC_PATH}/musica_ccpp_util.F90
${CCPP_TEST_SRC_PATH}/ccpp_kinds.F90
${TO_BE_CCPPIZED_SRC_PATH}/shr_orb_mod.F90
)

target_link_libraries(test_musica_ccpp_util
PRIVATE
musica::musica-fortran
)

set_target_properties(test_musica_ccpp_util
PROPERTIES
LINKER_LANGUAGE Fortran
)

add_test(
NAME test_musica_ccpp_util
COMMAND $<TARGET_FILE:test_musica_ccpp_util>
WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
)

# ---------------------------------------------------------
# Create a test for MUSICA CCPP wrapper
# ---------------------------------------------------------
Expand All @@ -30,6 +56,7 @@ target_sources(test_musica_api
${MUSICA_CCPP_SOURCES}
${TO_BE_CCPPIZED_SRC_PATH}/ccpp_tuvx_utils.F90
${TO_BE_CCPPIZED_SRC_PATH}/ccpp_const_utils.F90
${TO_BE_CCPPIZED_SRC_PATH}/shr_orb_mod.F90
${CCPP_SRC_PATH}/ccpp_constituent_prop_mod.F90
${CCPP_SRC_PATH}/ccpp_hash_table.F90
${CCPP_SRC_PATH}/ccpp_hashable.F90
Expand Down
Loading