Skip to content

Commit

Permalink
update to use solar zenith angle and earth-sun distance from host model
Browse files Browse the repository at this point in the history
  • Loading branch information
@mattldawson
mattldawson committed Dec 11, 2024
1 parent d7220d6 commit 174190e
Show file tree
Hide file tree
Showing 6 changed files with 21 additions and 227 deletions.
23 changes: 4 additions & 19 deletions schemes/musica/musica_ccpp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,14 +77,12 @@ 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, latitude, longitude, &
earth_eccentricity, earth_obliquity, perihelion_longitude, &
moving_vernal_equinox_longitude, calendar_day, 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
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 @@ -104,13 +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) :: 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
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 @@ -119,16 +112,8 @@ 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 Down
38 changes: 4 additions & 34 deletions schemes/musica/musica_ccpp.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -183,46 +183,16 @@
type = real | kind = kind_phys
dimensions = (horizontal_loop_extent,vertical_layer_dimension)
intent = in
[ latitude ]
standard_name = latitude
[ solar_zenith_angle ]
standard_name = solar_zenith_angle
units = radians
type = real | kind = kind_phys
dimensions = (horizontal_loop_extent)
intent = in
[ longitude ]
standard_name = longitude
[ earth_sun_distance ]
standard_name = earth_sun_distance
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 ]
Expand Down
51 changes: 1 addition & 50 deletions schemes/musica/musica_ccpp_util.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ module musica_ccpp_util
implicit none

private
public :: has_error_occurred, calculate_solar_zenith_angle_and_earth_sun_distance
public :: has_error_occurred

real(kind_phys), parameter, public :: PI = 3.14159265358979323846_kind_phys
real(kind_phys), parameter, public :: DEGREE_TO_RADIAN = PI / 180.0_kind_phys
Expand Down Expand Up @@ -42,53 +42,4 @@ 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
26 changes: 0 additions & 26 deletions test/musica/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,32 +13,6 @@ 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 Down
52 changes: 12 additions & 40 deletions test/musica/test_musica_api.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,8 +164,6 @@ subroutine test_chapman()
real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS+1) :: geopotential_height_wrt_surface_at_interface ! m
real(kind_phys), dimension(NUM_COLUMNS) :: surface_geopotential ! m2 s-2
real(kind_phys), dimension(NUM_COLUMNS) :: surface_temperature ! K
real(kind_phys), dimension(NUM_COLUMNS) :: latitude ! radians
real(kind_phys), dimension(NUM_COLUMNS) :: longitude ! radians
real(kind_phys) :: surface_albedo ! unitless
integer, parameter :: num_photolysis_wavelength_grid_sections = 8 ! (count)
real(kind_phys), dimension(num_photolysis_wavelength_grid_sections+1) :: flux_data_photolysis_wavelength_interfaces ! nm
Expand All @@ -180,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 All @@ -190,11 +190,6 @@ subroutine test_chapman()
integer :: i, j, k
integer :: N2_index, O2_index, O_index, O1D_index, O3_index
real(kind_phys) :: total_O, total_O_init
real(kind_phys) :: earth_eccentricity
real(kind_phys) :: earth_obliquity
real(kind_phys) :: perihelion_longitude
real(kind_phys) :: moving_vernal_equinox_longitude
real(kind_phys) :: calendar_day

call get_wavelength_edges(photolysis_wavelength_grid_interfaces)
solver_type = Rosenbrock
Expand Down Expand Up @@ -223,16 +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 /)

! Set conditions for one daytime and one nighttime column
! Greenwich, UK and Wellington, NZ
latitude = (/ 51.5_kind_phys, -41.3_kind_phys /)
longitude = (/ 0.0_kind_phys, 174.8_kind_phys /)
earth_eccentricity = 0.0167_kind_phys
earth_obliquity = 23.5_kind_phys * DEGREE_TO_RADIAN
perihelion_longitude = 102.9_kind_phys * DEGREE_TO_RADIAN
moving_vernal_equinox_longitude = 182.7_kind_phys * DEGREE_TO_RADIAN
calendar_day = 183.5_kind_phys ! noon GMT Jul 1
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 @@ -337,9 +324,8 @@ subroutine test_chapman()
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, latitude, longitude, earth_eccentricity, &
earth_obliquity, perihelion_longitude, moving_vernal_equinox_longitude, &
calendar_day, errmsg, errcode )
air_pressure_thickness, solar_zenith_angle, earth_sun_distance, errmsg, &
errcode )
if (errcode /= 0) then
write(*,*) trim(errmsg)
stop 3
Expand Down Expand Up @@ -416,8 +402,6 @@ subroutine test_terminator()
real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS+1) :: geopotential_height_wrt_surface_at_interface ! m
real(kind_phys), dimension(NUM_COLUMNS) :: surface_geopotential ! m2 s-2
real(kind_phys), dimension(NUM_COLUMNS) :: surface_temperature ! K
real(kind_phys), dimension(NUM_COLUMNS) :: latitude ! radians
real(kind_phys), dimension(NUM_COLUMNS) :: longitude ! radians
real(kind_phys) :: surface_albedo ! unitless
integer, parameter :: num_photolysis_wavelength_grid_sections = 8 ! (count)
real(kind_phys), dimension(num_photolysis_wavelength_grid_sections+1) :: flux_data_photolysis_wavelength_interfaces ! nm
Expand All @@ -432,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 All @@ -442,11 +428,6 @@ subroutine test_terminator()
integer :: i, j, k
integer :: Cl_index, Cl2_index
real(kind_phys) :: total_Cl, total_Cl_init
real(kind_phys) :: earth_eccentricity
real(kind_phys) :: earth_obliquity
real(kind_phys) :: perihelion_longitude
real(kind_phys) :: moving_vernal_equinox_longitude
real(kind_phys) :: calendar_day

call get_wavelength_edges(photolysis_wavelength_grid_interfaces)
solver_type = Rosenbrock
Expand Down Expand Up @@ -475,16 +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 /)

! Set conditions for one daytime and one nighttime column
! Greenwich, UK and Wellington, NZ
latitude = (/ 51.5_kind_phys, -41.3_kind_phys /)
longitude = (/ 0.0_kind_phys, 174.8_kind_phys /)
earth_eccentricity = 0.0167_kind_phys
earth_obliquity = 23.5_kind_phys * DEGREE_TO_RADIAN
perihelion_longitude = 102.9_kind_phys * DEGREE_TO_RADIAN
moving_vernal_equinox_longitude = 182.7_kind_phys * DEGREE_TO_RADIAN
calendar_day = 183.5_kind_phys ! noon GMT Jul 1
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 @@ -577,9 +550,8 @@ subroutine test_terminator()
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, latitude, longitude, earth_eccentricity, &
earth_obliquity, perihelion_longitude, moving_vernal_equinox_longitude, &
calendar_day, errmsg, errcode )
air_pressure_thickness, solar_zenith_angle, earth_sun_distance, errmsg, &
errcode )
if (errcode /= 0) then
write(*,*) trim(errmsg)
stop 3
Expand Down
Loading

0 comments on commit 174190e

Please sign in to comment.