diff --git a/schemes/musica/musica_ccpp.F90 b/schemes/musica/musica_ccpp.F90 index 4db0f18c..84b3b441 100644 --- a/schemes/musica/musica_ccpp.F90 +++ b/schemes/musica/musica_ccpp.F90 @@ -43,7 +43,10 @@ end subroutine musica_ccpp_init !> \section arg_table_musica_ccpp_run Argument Table !! \htmlinclude musica_ccpp_run.html subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, constituent_props, & - constituents, height_midpoints, height_interfaces, errmsg, errcode) + constituents, geopotential_height_wrt_surface_at_midpoint, & + geopotential_height_wrt_surface_at_interface, & + surface_geopotential, reciprocal_of_gravitational_acceleration, & + errmsg, errcode) use musica_ccpp_micm_util, only: reshape_into_micm_arr, reshape_into_ccpp_arr use musica_ccpp_micm_util, only: convert_to_mol_per_cubic_meter, convert_to_mass_mixing_ratio use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t @@ -55,8 +58,10 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co real(kind_phys), target, intent(in) :: dry_air_density(:,:) ! kg m-3 type(ccpp_constituent_prop_ptr_t), intent(in) :: constituent_props(:) real(kind_phys), target, intent(inout) :: constituents(:,:,:) ! kg kg-1 - real(kind_phys), target, intent(in) :: height_midpoints(:,:) ! km - real(kind_phys), target, intent(in) :: height_interfaces(:,:) ! km + real(kind_phys), target, intent(in) :: geopotential_height_wrt_surface_at_midpoint(:,:) ! m + real(kind_phys), target, intent(in) :: geopotential_height_wrt_surface_at_interface(:,:) ! m + real(kind_phys), target, intent(in) :: surface_geopotential(:) ! m2 s-2 + real(kind_phys), target, intent(in) :: reciprocal_of_gravitational_acceleration ! s2 m-1 character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode @@ -78,7 +83,9 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co * 3) :: photolysis_rate_constants ! s-1 integer :: i_elem - call tuvx_run(temperature, dry_air_density, height_midpoints, height_interfaces, & + call tuvx_run(temperature, dry_air_density, geopotential_height_wrt_surface_at_midpoint, & + geopotential_height_wrt_surface_at_interface, & + surface_geopotential, reciprocal_of_gravitational_acceleration, & photolysis_rate_constants, errmsg, errcode) ! Get the molar_mass that is set in the call to instantiate() diff --git a/schemes/musica/musica_ccpp.meta b/schemes/musica/musica_ccpp.meta index 331f8500..8cca3abb 100644 --- a/schemes/musica/musica_ccpp.meta +++ b/schemes/musica/musica_ccpp.meta @@ -71,18 +71,30 @@ type = real | kind = kind_phys dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_ccpp_constituents) intent = inout -[ height_midpoints ] - standard_name = ccpp_height_midpoints +[ geopotential_height_wrt_surface_at_midpoint ] + standard_name = geopotential_height_wrt_surface units = km type = real | kind = kind_phys dimensions = (horizontal_loop_extent,vertical_layer_dimension) intent = in -[ height_interfaces ] - standard_name = ccpp_height_interfaces +[ geopotential_height_wrt_surface_at_interface ] + standard_name = geopotential_height_wrt_surface_at_interface units = km type = real | kind = kind_phys dimensions = (horizontal_loop_extent,vertical_layer_dimension) intent = in +[ surface_geopotential ] + standard_name = surface_geopotential + type = real | kind = kind_phys + units = m2 s-2 + dimensions = (horizontal_loop_extent) + intent = in +[ reciprocal_of_gravitational_acceleration ] + standard_name = reciprocal_of_gravitational_acceleration + units = s2 m-1 + type = real | kind = kind_phys + dimensions = () + intent = in [ errmsg ] standard_name = ccpp_error_message units = none diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 index 94799f0c..8ebdc2d0 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 @@ -96,15 +96,20 @@ subroutine tuvx_init(vertical_layer_dimension, & end subroutine tuvx_init !> Calculates photolysis rate constants for the current model conditions - subroutine tuvx_run( temperature, dry_air_density, height_midpoints, & - height_interfaces, photolysis_rate_constants, errmsg, errcode ) + subroutine tuvx_run( temperature, dry_air_density, & + geopotential_height_wrt_surface_at_midpoint, & + geopotential_height_wrt_surface_at_interface, & + surface_geopotential, reciprocal_of_gravitational_acceleration, & + photolysis_rate_constants, errmsg, errcode ) use musica_util, only: error_t - use musica_ccpp_tuvx_height_grid, only: set_height_grid_values + use musica_ccpp_tuvx_height_grid, only: set_height_grid_values, calculate_heights real(kind_phys), intent(in) :: temperature(:,:) ! K (column, layer) real(kind_phys), intent(in) :: dry_air_density(:,:) ! molecule cm-3 (column, layer) - real(kind_phys), intent(in) :: height_midpoints(:,:) ! km (column, layer) - real(kind_phys), intent(in) :: height_interfaces(:,:) ! km (column, interface) + 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 + real(kind_phys), intent(in) :: reciprocal_of_gravitational_acceleration ! s2 m-1 ! temporarily set to Chapman mechanism and 1 dimension ! until mapping between MICM and TUV-x is implemented real(kind_phys), intent(out) :: photolysis_rate_constants(:) ! s-1 (column, reaction) @@ -113,14 +118,20 @@ subroutine tuvx_run( temperature, dry_air_density, height_midpoints, & ! local variables type(error_t) :: error + 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 integer :: i_col errcode = 0 errmsg = '' do i_col = 1, size(temperature, dim=1) - call set_height_grid_values( height_grid, height_midpoints(i_col,:), & - height_interfaces(i_col,:), errmsg, errcode ) + 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 ) end do if (errcode /= 0) return diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_height_grid.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_height_grid.F90 index ce085cb5..f61d2d19 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_height_grid.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_height_grid.F90 @@ -3,7 +3,7 @@ module musica_ccpp_tuvx_height_grid implicit none private - public :: create_height_grid, set_height_grid_values + public :: create_height_grid, set_height_grid_values, calculate_heights ! Conversions between the CAM-SIMA height grid and the TUVX height grid ! @@ -135,4 +135,39 @@ subroutine set_height_grid_values( height_grid, host_midpoints, & end subroutine set_height_grid_values + !> Calculates the heights needed for the TUV-x grid based on available data + !! + !! Uses the reciprocal of gravitational acceleration, the surface geopotential, + !! and the geopotential height wrt the surface and midpoints/interfaces to calculate + !! the midpoint/interface height above sea level in km needed for the TUV-x grid. + !! + !! The equation used is taked from CAMChem + !! (see https://github.com/ESCOMP/CAM/blob/f0e489e9708ce7b91635f6d4997fbf1e390b0dbb/src/chemistry/mozart/mo_gas_phase_chemdr.F90#L514-L526) + subroutine calculate_heights( geopotential_height_wrt_surface_at_midpoint, & + geopotential_height_wrt_surface_at_interface, & + surface_geopotential, reciprocal_of_gravitational_acceleration, & + height_midpoints, height_interfaces ) + + use ccpp_kinds, only: kind_phys + use musica_ccpp_util, only: has_error_occurred + use musica_util, only: error_t + + real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_midpoint(:) ! m + real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_interface(:) ! m + real(kind_phys), intent(in) :: surface_geopotential ! m2 s-2 + real(kind_phys), intent(in) :: reciprocal_of_gravitational_acceleration ! s2 m-1 + real(kind_phys), intent(out) :: height_midpoints(:) ! Pa + real(kind_phys), intent(out) :: height_interfaces(:) ! Pa + real(kind_phys) :: surface_height ! km + + surface_height = surface_geopotential * reciprocal_of_gravitational_acceleration + height_midpoints(:) = 0.001_kind_phys * & + ( geopotential_height_wrt_surface_at_midpoint(:) & + + surface_height ) + height_interfaces(:) = 0.001_kind_phys * & + ( geopotential_height_wrt_surface_at_interface(:) & + + surface_height ) + + end subroutine calculate_heights + end module musica_ccpp_tuvx_height_grid diff --git a/test/musica/test_musica_api.F90 b/test/musica/test_musica_api.F90 index ad6d8b40..3fc25dc7 100644 --- a/test/musica/test_musica_api.F90 +++ b/test/musica/test_musica_api.F90 @@ -19,8 +19,10 @@ subroutine test_musica_ccpp_api() integer :: errcode character(len=512) :: errmsg real(kind_phys) :: time_step ! s - real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: height_midpoints ! km - real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS+1) :: height_interfaces ! km + real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: geopotential_height_wrt_surface_at_midpoint ! m + 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) :: reciprocal_of_gravitational_acceleration ! s2 m-1 real(kind_phys), target, dimension(NUM_COLUMNS,NUM_LAYERS) :: temperature ! K real(kind_phys), target, dimension(NUM_COLUMNS,NUM_LAYERS) :: pressure ! Pa real(kind_phys), target, dimension(NUM_COLUMNS,NUM_LAYERS) :: dry_air_density ! kg m-3 @@ -39,10 +41,12 @@ subroutine test_musica_ccpp_api() solver_type = Rosenbrock num_grid_cells = NUM_COLUMNS * NUM_LAYERS time_step = 60._kind_phys - height_midpoints(1,:) = (/ 3.0_kind_phys, 1.0_kind_phys /) - height_midpoints(2,:) = (/ 4.0_kind_phys, 1.5_kind_phys /) - height_interfaces(1,:) = (/ 4.0_kind_phys, 2.0_kind_phys, 1.0_kind_phys /) - height_interfaces(2,:) = (/ 5.0_kind_phys, 2.5_kind_phys, 0.5_kind_phys /) + geopotential_height_wrt_surface_at_midpoint(1,:) = (/ 2000.0_kind_phys, 500.0_kind_phys /) + geopotential_height_wrt_surface_at_midpoint(2,:) = (/ 2000.0_kind_phys, -500.0_kind_phys /) + geopotential_height_wrt_surface_at_interface(1,:) = (/ 3000.0_kind_phys, 1000.0_kind_phys, 0.0_kind_phys /) + geopotential_height_wrt_surface_at_interface(2,:) = (/ 3000.0_kind_phys, 500.0_kind_phys, -1500.0_kind_phys /) + surface_geopotential = (/ 100.0_kind_phys, 200.0_kind_phys /) + reciprocal_of_gravitational_acceleration = 10.0_kind_phys temperature(:,1) = (/ 100._kind_phys, 200._kind_phys /) temperature(:,2) = (/ 300._kind_phys, 400._kind_phys /) pressure(:,1) = (/ 6000.04_kind_phys, 7000.04_kind_phys /) @@ -102,7 +106,10 @@ subroutine test_musica_ccpp_api() write(*,fmt="(4(3x,e13.6))") constituents call musica_ccpp_run(time_step, temperature, pressure, dry_air_density, constituent_props_ptr, & - constituents, height_midpoints, height_interfaces, errmsg, errcode) + constituents, geopotential_height_wrt_surface_at_midpoint, & + geopotential_height_wrt_surface_at_interface, & + surface_geopotential, reciprocal_of_gravitational_acceleration, & + errmsg, errcode) if (errcode /= 0) then write(*,*) trim(errmsg) stop 3 diff --git a/test/musica/tuvx/test_tuvx_height_grid.F90 b/test/musica/tuvx/test_tuvx_height_grid.F90 index 9f0fed15..19339c0a 100644 --- a/test/musica/tuvx/test_tuvx_height_grid.F90 +++ b/test/musica/tuvx/test_tuvx_height_grid.F90 @@ -10,6 +10,7 @@ program test_tuvx_height_grid #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 call test_create_height_grid() + call test_calculate_height_grid_values() contains @@ -73,4 +74,34 @@ subroutine test_create_height_grid() end subroutine test_create_height_grid + subroutine test_calculate_height_grid_values() + + use ccpp_kinds, only: kind_phys + + integer, parameter :: NUM_LAYERS = 2 + real(kind_phys), dimension(NUM_LAYERS) :: geopotential_height_wrt_surface_at_midpoint ! m + real(kind_phys), dimension(NUM_LAYERS+1) :: geopotential_height_wrt_surface_at_interface ! m + real(kind_phys) :: surface_geopotential ! m2 s-2 + real(kind_phys) :: reciprocal_of_gravitational_acceleration ! s2 m-1 + real(kind_phys), dimension(NUM_LAYERS) :: height_midpoints + real(kind_phys), dimension(NUM_LAYERS+1) :: height_interfaces + + geopotential_height_wrt_surface_at_midpoint(:) = (/ 2000.0_kind_phys, 500.0_kind_phys /) + geopotential_height_wrt_surface_at_interface(:) = (/ 3000.0_kind_phys, 1000.0_kind_phys, 0.0_kind_phys /) + surface_geopotential = 100.0_kind_phys + reciprocal_of_gravitational_acceleration = 10.0_kind_phys + + call calculate_heights(geopotential_height_wrt_surface_at_midpoint, & + geopotential_height_wrt_surface_at_interface, & + surface_geopotential, reciprocal_of_gravitational_acceleration, & + height_midpoints, height_interfaces) + + ASSERT_NEAR(height_midpoints(1), 3.0, 1e-5) + ASSERT_NEAR(height_midpoints(2), 1.5, 1e-5) + ASSERT_NEAR(height_interfaces(1), 4.0, 1e-5) + ASSERT_NEAR(height_interfaces(2), 2.0, 1e-5) + ASSERT_NEAR(height_interfaces(3), 1.0, 1e-5) + + end subroutine test_calculate_height_grid_values + end program test_tuvx_height_grid