diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 5d7bf4a..238e0c1 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -12,4 +12,14 @@ jobs: - name: build Docker image run: docker build -t musica -f test/docker/Dockerfile.musica . - name: run tests in container - run: docker run --name test-container -t musica bash -c 'make test ARGS="--rerun-failed --output-on-failure -j8"' \ No newline at end of file + run: docker run --name test-container -t musica bash -c 'make test ARGS="--rerun-failed --output-on-failure -j8"' + test_musica_api_no_install: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + with: + submodules: recursive + - name: build Docker image + run: docker build -t musica-no-install -f test/docker/Dockerfile.musica.no_install . + - name: run tests in container + run: docker run --name test-container -t musica-no-install bash -c 'make test ARGS="--rerun-failed --output-on-failure -j8"' \ No newline at end of file diff --git a/schemes/musica/micm/musica_ccpp_micm.F90 b/schemes/musica/micm/musica_ccpp_micm.F90 index d11e86d..37f6673 100644 --- a/schemes/musica/micm/musica_ccpp_micm.F90 +++ b/schemes/musica/micm/musica_ccpp_micm.F90 @@ -30,10 +30,11 @@ subroutine micm_register(constituents, solver_type, num_grid_cells, errmsg, errc integer, intent(out) :: errcode ! local variables - type(error_t) :: error - real(kind=kind_phys) :: molar_mass - logical :: is_advected - integer :: i + type(error_t) :: error + real(kind=kind_phys) :: molar_mass + character(len=:), allocatable :: species_name + logical :: is_advected + integer :: i, species_index errcode = 0 errmsg = '' @@ -41,26 +42,29 @@ subroutine micm_register(constituents, solver_type, num_grid_cells, errmsg, errc micm => micm_t(filename_of_micm_configuration, solver_type, num_grid_cells, error) if (has_error_occurred(error, errmsg, errcode)) return - allocate(constituents(size(micm%species_ordering)), stat=errcode) + allocate(constituents(micm%species_ordering%size()), stat=errcode) if (errcode /= 0) then errmsg = "[MUSICA Error] Failed to allocate memory for constituents." return end if - do i = 1, size(micm%species_ordering) - associate( map => micm%species_ordering(i) ) - molar_mass = micm%get_species_property_double(map%name(), & + do i = 1, micm%species_ordering%size() + associate( map => micm%species_ordering ) + species_name = map%name(i) + species_index = map%index(i) + + molar_mass = micm%get_species_property_double(species_name, & "molecular weight [kg mol-1]", & error) if (has_error_occurred(error, errmsg, errcode)) return - is_advected = micm%get_species_property_bool(map%name(), & - "__is advected", & - error) + is_advected = micm%get_species_property_bool(species_name, & + "__is advected", & + error) if (has_error_occurred(error, errmsg, errcode)) return - call constituents(map%index())%instantiate( & - std_name = map%name(), & - long_name = map%name(), & + call constituents(species_index)%instantiate( & + std_name = species_name, & + long_name = species_name, & units = 'kg kg-1', & vertical_dim = 'vertical_layer_dimension', & default_value = 0.0_kind_phys, & @@ -87,7 +91,7 @@ end subroutine micm_init !> Solve chemistry at the current time step subroutine micm_run(time_step, temperature, pressure, dry_air_density, constituents, & - rate_params, errmsg, errcode) + user_defined_rate_parameters, errmsg, errcode) use musica_micm, only: solver_stats_t use musica_util, only: string_t, error_t @@ -96,7 +100,7 @@ subroutine micm_run(time_step, temperature, pressure, dry_air_density, constitue real(c_double), target, intent(in) :: pressure(:) ! Pa real(c_double), target, intent(in) :: dry_air_density(:) ! kg m-3 real(c_double), target, intent(inout) :: constituents(:) ! mol m-3 - real(c_double), target, intent(inout) :: rate_params(:) + real(c_double), target, intent(in) :: user_defined_rate_parameters(:) ! various units character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode @@ -111,14 +115,14 @@ subroutine micm_run(time_step, temperature, pressure, dry_air_density, constitue errmsg = '' c_time_step = real(time_step, c_double) - call micm%solve(c_time_step, & - temperature, & - pressure, & - dry_air_density, & - constituents, & - rate_params, & - solver_state, & - solver_stats, & + call micm%solve(c_time_step, & + temperature, & + pressure, & + dry_air_density, & + constituents, & + user_defined_rate_parameters, & + solver_state, & + solver_stats, & error) if (has_error_occurred(error, errmsg, errcode)) return diff --git a/schemes/musica/micm/micm_util.F90 b/schemes/musica/micm/musica_ccpp_micm_util.F90 similarity index 70% rename from schemes/musica/micm/micm_util.F90 rename to schemes/musica/micm/musica_ccpp_micm_util.F90 index 3272ee9..08ed0f0 100644 --- a/schemes/musica/micm/micm_util.F90 +++ b/schemes/musica/micm/musica_ccpp_micm_util.F90 @@ -1,4 +1,4 @@ -module micm_util +module musica_ccpp_micm_util implicit none private @@ -8,7 +8,7 @@ module micm_util contains subroutine reshape_into_micm_arr(temperature, pressure, dry_air_density, constituents, & - rate_params, m_temperature, m_pressure, m_dry_air_density, m_constituents, m_rate_params) + m_temperature, m_pressure, m_dry_air_density, m_constituents) use iso_c_binding, only: c_double use ccpp_kinds, only: kind_phys @@ -16,28 +16,24 @@ subroutine reshape_into_micm_arr(temperature, pressure, dry_air_density, constit real(kind_phys), target, intent(in) :: pressure(:,:) ! Pa real(kind_phys), target, intent(in) :: dry_air_density(:,:) ! kg m-3 real(kind_phys), target, intent(in) :: constituents(:,:,:) ! kg kg-1 - real(kind_phys), target, intent(in) :: rate_params(:,:,:) real(c_double), target, intent(out) :: m_temperature(:) ! K real(c_double), target, intent(out) :: m_pressure(:) ! Pa real(c_double), target, intent(out) :: m_dry_air_density(:) ! kg m-3 real(c_double), target, intent(out) :: m_constituents(:) ! kg kg-1 - real(c_double), target, intent(out) :: m_rate_params(:) ! local variables integer :: num_columns, num_layers - integer :: num_constituents, num_rate_params - integer :: i_column, i_layer, i_elem, i_constituents, i_rate_params + integer :: num_constituents + integer :: i_column, i_layer, i_elem, i_constituents num_columns = size(constituents, dim=1) num_layers = size(constituents, dim=2) num_constituents = size(constituents, dim=3) - num_rate_params = size(rate_params, dim=3) ! Reshape into 1-D arry in species-column first order ! refers to: state.variables_[i_cell][i_species] = concentrations[i_species_elem++] i_elem = 1 i_constituents = 1 - i_rate_params = 1 do i_layer = 1, num_layers do i_column = 1, num_columns m_temperature(i_elem) = real(temperature(i_column, i_layer), c_double) @@ -45,56 +41,34 @@ subroutine reshape_into_micm_arr(temperature, pressure, dry_air_density, constit m_dry_air_density(i_elem) = real(dry_air_density(i_column, i_layer), c_double) m_constituents(i_constituents : i_constituents + num_constituents - 1) & = real(constituents(i_column, i_layer, :), c_double) - m_rate_params(i_rate_params : i_rate_params + num_rate_params - 1) & - = real(rate_params(i_column, i_layer, :), c_double) i_elem = i_elem + 1 i_constituents = i_constituents + num_constituents - i_rate_params = i_rate_params + num_rate_params end do end do end subroutine reshape_into_micm_arr - subroutine reshape_into_ccpp_arr(temperature, pressure, dry_air_density, constituents, & - rate_params, m_temperature, m_pressure, m_dry_air_density, m_constituents, m_rate_params) + subroutine reshape_into_ccpp_arr(constituents, m_constituents) use iso_c_binding, only: c_double use ccpp_kinds, only: kind_phys - real(kind_phys), intent(out) :: temperature(:,:) ! K - real(kind_phys), intent(out) :: pressure(:,:) ! Pa - real(kind_phys), intent(out) :: dry_air_density(:,:) ! kg m-3 real(kind_phys), intent(out) :: constituents(:,:,:) ! kg kg-1 - real(kind_phys), intent(out) :: rate_params(:,:,:) - real(c_double), intent(in) :: m_temperature(:) ! K - real(c_double), intent(in) :: m_pressure(:) ! Pa - real(c_double), intent(in) :: m_dry_air_density(:) ! kg m-3 real(c_double), intent(in) :: m_constituents(:) ! kg kg-1 - real(c_double), intent(in) :: m_rate_params(:) ! local variables integer :: num_columns, num_layers - integer :: num_constituents, num_rate_params - integer :: i_column, i_layer, i_elem, i_constituents, i_rate_params + integer :: num_constituents + integer :: i_column, i_layer, i_constituents num_columns = size(constituents, dim=1) num_layers = size(constituents, dim=2) num_constituents = size(constituents, dim=3) - num_rate_params = size(rate_params, dim=3) - i_elem = 1 i_constituents = 1 - i_rate_params = 1 do i_layer = 1, num_layers do i_column = 1, num_columns - temperature(i_column, i_layer) = real(m_temperature(i_elem), kind_phys) - pressure(i_column, i_layer) = real(m_pressure(i_elem), kind_phys) - dry_air_density(i_column, i_layer) = real(m_dry_air_density(i_elem), kind_phys) constituents(i_column, i_layer, :) & = real(m_constituents(i_constituents : i_constituents + num_constituents - 1), kind_phys) - rate_params(i_column, i_layer, :) & - = real(m_rate_params(i_rate_params : i_rate_params + num_rate_params - 1), kind_phys) - i_elem = i_elem + 1 i_constituents = i_constituents + num_constituents - i_rate_params = i_rate_params + num_rate_params end do end do @@ -156,4 +130,4 @@ subroutine convert_to_mass_mixing_ratio(dry_air_density, molar_mass_arr, constit end subroutine convert_to_mass_mixing_ratio -end module micm_util \ No newline at end of file +end module musica_ccpp_micm_util \ No newline at end of file diff --git a/schemes/musica/musica_ccpp.F90 b/schemes/musica/musica_ccpp.F90 index bf4fae3..84b3b44 100644 --- a/schemes/musica/musica_ccpp.F90 +++ b/schemes/musica/musica_ccpp.F90 @@ -25,32 +25,43 @@ end subroutine musica_ccpp_register !> \section arg_table_musica_ccpp_init Argument Table !! \htmlinclude musica_ccpp_init.html - subroutine musica_ccpp_init(errmsg, errcode) + subroutine musica_ccpp_init(vertical_layer_dimension, vertical_interface_dimension, & + errmsg, errcode) + integer, intent(in) :: vertical_layer_dimension ! (count) + integer, intent(in) :: vertical_interface_dimension ! (count) character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode - call tuvx_init(errmsg, errcode) + call tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & + errmsg, errcode) + if (errcode /= 0) return call micm_init(errmsg, errcode) + if (errcode /= 0) return 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, rate_params, height, errmsg, errcode) - use micm_util, only: reshape_into_micm_arr, reshape_into_ccpp_arr - use micm_util, only: convert_to_mol_per_cubic_meter, convert_to_mass_mixing_ratio + 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 use ccpp_kinds, only: kind_phys use iso_c_binding, only: c_double - real(kind_phys), intent(inout) :: time_step ! s - real(kind_phys), target, intent(inout) :: temperature(:,:) ! K - real(kind_phys), target, intent(inout) :: pressure(:,:) ! Pa - real(kind_phys), target, intent(inout) :: dry_air_density(:,:) ! kg m-3 + real(kind_phys), intent(in) :: time_step ! s + real(kind_phys), target, intent(in) :: temperature(:,:) ! K + real(kind_phys), target, intent(in) :: pressure(:,:) ! Pa + 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(inout) :: rate_params(:,:,:) - real(kind_phys), target, intent(in) :: height(:,:) ! km + real(kind_phys), target, intent(inout) :: constituents(:,:,:) ! kg kg-1 + 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 @@ -64,13 +75,18 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co real(c_double), target, dimension(size(constituents, dim=1) & * size(constituents, dim=2) & * size(constituents, dim=3)) :: m_constituents ! mol m-3 - real(c_double), target, dimension(size(rate_params, dim=1) & - * size(rate_params, dim=2) & - * size(rate_params, dim=3)) :: m_rate_params real(kind_phys), target, dimension(size(constituents, dim=3)) :: molar_mass_arr ! kg mol-1 + + ! temporarily dimensioned to Chapman mechanism until mapping between MICM and TUV-x is implemented + real(c_double), target, dimension(size(constituents, dim=1) & + * size(constituents, dim=2) & + * 3) :: photolysis_rate_constants ! s-1 integer :: i_elem - call tuvx_run(height, temperature, dry_air_density, errmsg, errcode) + 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() do i_elem = 1, size(molar_mass_arr) @@ -95,15 +111,15 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co call convert_to_mol_per_cubic_meter(dry_air_density, molar_mass_arr, constituents) ! Reshape array (3D -> 1D) and convert type (kind_phys -> c_double) - call reshape_into_micm_arr(temperature, pressure, dry_air_density, constituents, rate_params, & - m_temperature, m_pressure, m_dry_air_density, m_constituents, m_rate_params) + call reshape_into_micm_arr(temperature, pressure, dry_air_density, constituents, & + m_temperature, m_pressure, m_dry_air_density, m_constituents) + ! temporarily pass in unmapped photolysis rate constants until mapping between MICM and TUV-x is implemented call micm_run(time_step, m_temperature, m_pressure, m_dry_air_density, m_constituents, & - m_rate_params, errmsg, errcode) + photolysis_rate_constants, errmsg, errcode) ! Reshape array (1D -> 3D) and convert type (c_double -> kind_phys) - call reshape_into_ccpp_arr(temperature, pressure, dry_air_density, constituents, rate_params, & - m_temperature, m_pressure, m_dry_air_density, m_constituents, m_rate_params) + call reshape_into_ccpp_arr(constituents, m_constituents) ! Convert MICM unit back to CAM-SIMA unit (mol m-3 -> kg kg-1) call convert_to_mass_mixing_ratio(dry_air_density, molar_mass_arr, constituents) diff --git a/schemes/musica/musica_ccpp.meta b/schemes/musica/musica_ccpp.meta index 30b7c7c..8cca3ab 100644 --- a/schemes/musica/musica_ccpp.meta +++ b/schemes/musica/musica_ccpp.meta @@ -1,12 +1,24 @@ [ccpp-table-properties] name = musica_ccpp type = scheme - dependencies = micm/musica_ccpp_micm.F90,micm/musica_ccpp_tuvx.F90,musica_ccpp_util.F90 + dependencies = micm/musica_ccpp_micm.F90,micm/musica_ccpp_micm_util.F90,tuvx/musica_ccpp_tuvx.F90,tuvx/musica_ccpp_tuvx_height_grid.F90,musica_ccpp_util.F90 dynamic_constituent_routine = musica_ccpp_register [ccpp-arg-table] name = musica_ccpp_init type = scheme +[ vertical_layer_dimension ] + standard_name = vertical_layer_dimension + units = none + type = integer + dimensions = () + intent = in +[ vertical_interface_dimension ] + standard_name = vertical_interface_dimension + units = none + type = integer + dimensions = () + intent = in [ errmsg ] standard_name = ccpp_error_message units = none @@ -59,6 +71,30 @@ type = real | kind = kind_phys dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_ccpp_constituents) intent = inout +[ 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 +[ 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 2413a47..8ebdc2d 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 @@ -3,7 +3,7 @@ module musica_ccpp_tuvx ! Note: "tuvx_t" is included in an external pre-built tuvx library that the host ! model is responsible for linking to during compilation - use musica_tuvx, only: tuvx_t + use musica_tuvx, only: tuvx_t, grid_t use musica_ccpp_util, only: has_error_occurred use ccpp_kinds, only: kind_phys use musica_ccpp_namelist, only: filename_of_tuvx_configuration @@ -14,14 +14,20 @@ module musica_ccpp_tuvx public :: tuvx_init, tuvx_run, tuvx_final type(tuvx_t), pointer :: tuvx => null( ) + type(grid_t), pointer :: height_grid => null( ) contains !> Intitialize TUVX - subroutine tuvx_init(errmsg, errcode) - use musica_tuvx, only: grid_map_t, profile_map_t, radiator_map_t + subroutine tuvx_init(vertical_layer_dimension, & + vertical_interface_dimension, errmsg, errcode) + use musica_tuvx, only: grid_map_t, grid_t, profile_map_t, radiator_map_t use musica_util, only: error_t, mapping_t + use musica_ccpp_tuvx_height_grid, only: create_height_grid, & + height_grid_label, height_grid_units + integer, intent(in) :: vertical_layer_dimension ! (count) + integer, intent(in) :: vertical_interface_dimension ! (count) character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode @@ -35,40 +41,103 @@ subroutine tuvx_init(errmsg, errcode) errmsg = '' grids => grid_map_t( error ) - if (has_error_occurred(error, errmsg, errcode)) return + if (has_error_occurred( error, errmsg, errcode )) return - profiles => profile_map_t( error ) - if (has_error_occurred(error, errmsg, errcode)) return - - radiators =>radiator_map_t( error ) - if (has_error_occurred(error, errmsg, errcode)) return - - tuvx => tuvx_t( filename_of_tuvx_configuration, error ) - if (has_error_occurred(error, errmsg, errcode)) return + height_grid => create_height_grid( vertical_layer_dimension, & + vertical_interface_dimension, errmsg, errcode ) + if (errcode /= 0) return + call grids%add( height_grid, error ) + if (has_error_occurred( error, errmsg, errcode )) return + profiles => profile_map_t( error ) + if (has_error_occurred( error, errmsg, errcode )) then + deallocate( grids ) + return + end if + + radiators => radiator_map_t( error ) + if (has_error_occurred( error, errmsg, errcode )) then + deallocate( grids ) + deallocate( profiles ) + return + end if + + tuvx => tuvx_t( filename_of_tuvx_configuration, grids, profiles, & + radiators, error ) + if (has_error_occurred( error, errmsg, errcode )) then + deallocate( grids ) + deallocate( profiles ) + deallocate( radiators ) + return + end if + + deallocate( height_grid ) deallocate( grids ) deallocate( profiles ) deallocate( radiators ) - deallocate( tuvx ) + + grids => tuvx%get_grids( error ) + if (has_error_occurred( error, errmsg, errcode )) then + deallocate( tuvx ) + tuvx => null() + return + end if + + height_grid => grids%get( height_grid_label, height_grid_units, error ) + if (has_error_occurred( error, errmsg, errcode )) then + deallocate( tuvx ) + tuvx => null() + deallocate( grids ) + return + end if + + deallocate( grids ) end subroutine tuvx_init !> Calculates photolysis rate constants for the current model conditions - subroutine tuvx_run( height, temperature, dry_air_density, 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 - - real(kind_phys), intent(in) :: height(:,:) ! km (layer, column) - real(kind_phys), intent(in) :: temperature(:,:) ! K (layer, column) - real(kind_phys), intent(in) :: dry_air_density(:,:) ! molecule cm-3 (layer, column) + 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) :: 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) character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode ! 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 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 + + ! stand-in until actual photolysis rate constants are calculated + photolysis_rate_constants(:) = 1.0e-6_kind_phys + end subroutine tuvx_run !> Finalize tuvx @@ -78,6 +147,7 @@ subroutine tuvx_final(errmsg, errcode) errcode = 0 errmsg = '' + deallocate( height_grid ) end subroutine tuvx_final diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_height_grid.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_height_grid.F90 new file mode 100644 index 0000000..f61d2d1 --- /dev/null +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_height_grid.F90 @@ -0,0 +1,173 @@ +module musica_ccpp_tuvx_height_grid + + implicit none + + private + public :: create_height_grid, set_height_grid_values, calculate_heights + + ! Conversions between the CAM-SIMA height grid and the TUVX height grid + ! + !----------------------------------------------------------------------- + ! Notes on the conversion between the host-model height grid and the TUVX + ! + ! TUV-x heights are "bottom-up" and require atmospheric constituent + ! concentrations at interfaces. Therefore, CAM-SIMA mid-points are used + ! as TUV-x grid interfaces, with an additional layer introduced between + ! the surface and the lowest CAM-SIMA mid-point, and a layer at the + ! top of the TUV-x grid to hold species densities above the top CAM-SIMA + ! mid-point. + ! + ! Here, + ! - i_int is the index of an interface + ! - i_mid is the index of a mid-point + ! - pver is the CCPP vertical_layer_dimension + ! - pver+1 is the CCPP vertical_interface_dimension + + ! ---- (interface) ===== (mid-point) + ! + ! CAM TUV-x + ! ************************ (exo values) ***************************** + ! ------(top)------ i_int = 1 -------(top)------ i_int = pver + 2 + ! ================== i_mid = pver + 1 + ! ================= i_mid = 1 ------------------ i_int = pver + 1 + ! ----------------- i_int = 2 ================== i_mid = pver + ! ------------------ i_int = pver + ! || + ! || || + ! || + ! ----------------- i_int = pver + ! ================= i_imd = pver ------------------ i_int = 2 + ! ================== i_mid = 1 + ! -----(ground)---- i_int = pver+1 -----(ground)----- i_int = 1 + ! + !----------------------------------------------------------------------- + + !> Label for height grid in TUV-x + character(len=*), parameter, public :: height_grid_label = "height" + !> Units for height grid in TUV-x + character(len=*), parameter, public :: height_grid_units = "km" + +contains + + !> Creates a TUVX height grid from the host-model height grid + function create_height_grid( vertical_layer_dimension, & + vertical_interface_dimension, errmsg, errcode ) result( height_grid ) + + use musica_ccpp_util, only: has_error_occurred + use musica_tuvx_grid, only: grid_t + use musica_util, only: error_t + + integer, intent(in) :: vertical_layer_dimension ! (count) + integer, intent(in) :: vertical_interface_dimension ! (count) + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errcode + type(grid_t), pointer :: height_grid + + type(error_t) :: error + + height_grid => null() + if ( vertical_layer_dimension < 1 ) then + errmsg = "[MUSICA Error] Invalid vertical_layer_dimension." + errcode = 1 + return + end if + if ( vertical_interface_dimension - vertical_layer_dimension /= 1 ) then + errmsg = "[MUSICA Error] Invalid vertical_interface_dimension." + errcode = 1 + return + end if + height_grid => grid_t( height_grid_label, height_grid_units, & + vertical_interface_dimension, error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + + end function create_height_grid + + !> Sets TUVX height grid values from the host-model height grid + subroutine set_height_grid_values( height_grid, host_midpoints, & + host_edges, errmsg, errcode ) + + use ccpp_kinds, only: kind_phys + use musica_ccpp_util, only: has_error_occurred + use musica_tuvx_grid, only: grid_t + use musica_util, only: error_t + + type(grid_t), intent(inout) :: height_grid + real(kind_phys), intent(in) :: host_midpoints(:) ! km + real(kind_phys), intent(in) :: host_edges(:) ! km + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errcode + + type(error_t) :: error + real(kind_phys) :: midpoints(size(host_midpoints)+1) + real(kind_phys) :: edges(size(host_edges)+1) + integer :: n_host_midpoints, n_host_edges + + if ( size(midpoints) /= height_grid%number_of_sections( error ) ) then + errmsg = "[MUSICA Error] Invalid size of TUV-x mid-point heights." + errcode = 1 + return + end if + if ( has_error_occurred( error, errmsg, errcode ) ) return + if ( size(edges) /= height_grid%number_of_sections( error ) + 1 ) then + errmsg = "[MUSICA Error] Invalid size of TUV-x interface heights." + errcode = 1 + return + end if + if ( has_error_occurred( error, errmsg, errcode ) ) return + + n_host_midpoints = size(host_midpoints) + n_host_edges = size(host_edges) + + edges(1) = host_edges(n_host_edges) + edges(2:n_host_edges) = host_midpoints(n_host_midpoints:1:-1) + edges(n_host_edges+1) = host_edges(1) + + midpoints(1) = 0.5_kind_phys * & + ( host_midpoints(n_host_midpoints) + host_edges(n_host_edges) ) + midpoints(2:n_host_midpoints) = host_edges(n_host_midpoints:2:-1) + midpoints(n_host_midpoints+1) = 0.5_kind_phys * & + ( edges(n_host_edges) + edges(n_host_edges+1) ) + + call height_grid%set_edges( edges, error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + call height_grid%set_midpoints( midpoints, error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + + 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/docker/Dockerfile.musica b/test/docker/Dockerfile.musica index 2ad24b7..d323ac7 100644 --- a/test/docker/Dockerfile.musica +++ b/test/docker/Dockerfile.musica @@ -1,6 +1,6 @@ FROM ubuntu:22.04 -ARG MUSICA_GIT_TAG=aa0854ecee54bd7a5aeb7ea1ba0eebb2cf656146 +ARG MUSICA_GIT_TAG=dbbdb22f5f2807e27c2695db85291951ba178634 RUN apt update \ && apt install -y sudo \ @@ -63,15 +63,18 @@ RUN cd musica \ COPY . atmospheric_physics RUN sudo chown -R test_user:test_user atmospheric_physics +# Must make ccpp-framework available before building test RUN cd atmospheric_physics/test \ && mkdir lib \ && cd lib \ && git clone -b add_const_interface --depth 1 https://github.com/peverwhee/ccpp-framework.git ENV CCPP_SRC_PATH="lib/ccpp-framework/src" +ENV CCPP_FORTRAN_TOOLS_PATH="lib/ccpp-framework/scripts/fortran_tools" RUN cd atmospheric_physics/test \ && cmake -S . \ -B build \ + -D CMAKE_BUILD_TYPE=Debug \ -D CCPP_ENABLE_MUSICA_TESTS=ON \ -D CCPP_ENABLE_MEMCHECK=ON \ && cmake --build ./build diff --git a/test/docker/Dockerfile.musica.no_install b/test/docker/Dockerfile.musica.no_install new file mode 100644 index 0000000..efae273 --- /dev/null +++ b/test/docker/Dockerfile.musica.no_install @@ -0,0 +1,72 @@ +FROM ubuntu:22.04 + +ARG MUSICA_GIT_TAG=dbbdb22f5f2807e27c2695db85291951ba178634 + +RUN apt update \ + && apt install -y sudo \ + && adduser test_user \ + && echo "test_user ALL=(root) NOPASSWD: ALL" >> /etc/sudoers.d/test_user \ + && chmod 0440 /etc/sudoers.d/test_user + +USER test_user +WORKDIR /home/test_user + +RUN sudo apt update \ + && sudo apt -y install \ + cmake \ + cmake-curses-gui \ + curl \ + g++ \ + gcc \ + gfortran \ + git \ + libblas-dev \ + liblapack-dev \ + lcov \ + libcurl4-openssl-dev \ + libhdf5-dev \ + libnetcdff-dev \ + libopenmpi-dev \ + m4 \ + make \ + nlohmann-json3-dev \ + openmpi-bin \ + python3 \ + tree \ + valgrind \ + vim \ + zlib1g-dev \ + && sudo apt clean + +ENV PATH="${PATH}:/usr/lib/openmpi/bin" + +ENV FC=mpif90 +ENV FFLAGS="-I/usr/include/" +ENV MUSICA_GIT_TAG=${MUSICA_GIT_TAG} + +COPY . atmospheric_physics +RUN sudo chown -R test_user:test_user atmospheric_physics + +# Must make ccpp-framework available before building test +RUN cd atmospheric_physics/test \ + && mkdir lib \ + && cd lib \ + && git clone -b add_const_interface --depth 1 https://github.com/peverwhee/ccpp-framework.git +ENV CCPP_SRC_PATH="lib/ccpp-framework/src" +ENV CCPP_FORTRAN_TOOLS_PATH="lib/ccpp-framework/scripts/fortran_tools" + +RUN cd atmospheric_physics/test \ + && cmake -S . \ + -B build \ + -D CMAKE_BUILD_TYPE=Debug \ + -D CCPP_ENABLE_MUSICA_TESTS=ON \ + -D CCPP_ENABLE_MEMCHECK=ON \ + && cmake --build ./build + +RUN cd atmospheric_physics/test \ + && mkdir third_party \ + && cd third_party \ + && git clone --depth 1 https://github.com/NCAR/tuv-x.git \ + && cp -r tuv-x/data ../build/data + +WORKDIR /home/test_user/atmospheric_physics/test/build \ No newline at end of file diff --git a/test/musica/CMakeLists.txt b/test/musica/CMakeLists.txt index 3fcdd84..487d181 100644 --- a/test/musica/CMakeLists.txt +++ b/test/musica/CMakeLists.txt @@ -5,8 +5,9 @@ FetchContent_Declare(musica GIT_REPOSITORY https://github.com/NCAR/musica.git GIT_TAG $ENV{MUSICA_GIT_TAG} # Set by docker ) +message(STATUS "Using MUSICA commit: $ENV{MUSICA_GIT_TAG}") -set(MUSICA_BUILD_C_CXX_INTERFACE OFF) +set(MUSICA_BUILD_C_CXX_INTERFACE ON) set(MUSICA_BUILD_FORTRAN_INTERFACE ON) set(MUSICA_ENABLE_TESTS OFF) set(MUSICA_ENABLE_INSTALL OFF) @@ -22,8 +23,9 @@ add_executable(test_musica_api test_musica_api.F90 musica_ccpp_namelist.F90) target_sources(test_musica_api PUBLIC ${MUSICA_SRC_PATH}/micm/musica_ccpp_micm.F90 - ${MUSICA_SRC_PATH}/micm/micm_util.F90 + ${MUSICA_SRC_PATH}/micm/musica_ccpp_micm_util.F90 ${MUSICA_SRC_PATH}/tuvx/musica_ccpp_tuvx.F90 + ${MUSICA_SRC_PATH}/tuvx/musica_ccpp_tuvx_height_grid.F90 ${MUSICA_SRC_PATH}/musica_ccpp.F90 ${MUSICA_SRC_PATH}/musica_ccpp_util.F90 ${CCPP_SRC_PATH}/ccpp_constituent_prop_mod.F90 @@ -63,4 +65,22 @@ add_custom_target( ${CMAKE_CURRENT_SOURCE_DIR}/tuvx/configs ${CMAKE_BINARY_DIR}/configs ) -add_subdirectory(micm) \ No newline at end of file +add_subdirectory(micm) +add_subdirectory(tuvx) + +# Test metdadata +find_package(Python REQUIRED) + +file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/metadata_test) +add_custom_target( + copy_metadata_test_files ALL ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/../../schemes/musica/musica_ccpp.meta + ${CMAKE_CURRENT_SOURCE_DIR}/../../schemes/musica/musica_ccpp.F90 + ${CMAKE_BINARY_DIR}/metadata_test +) + +add_test( + NAME test_metadata + COMMAND ${Python_EXECUTABLE} ${CMAKE_BINARY_DIR}/../$ENV{CCPP_FORTRAN_TOOLS_PATH}/offline_check_fortran_vs_metadata.py + --directory ${CMAKE_BINARY_DIR}/metadata_test +) \ No newline at end of file diff --git a/test/musica/micm/CMakeLists.txt b/test/musica/micm/CMakeLists.txt index d239da2..ed2cf63 100644 --- a/test/musica/micm/CMakeLists.txt +++ b/test/musica/micm/CMakeLists.txt @@ -3,7 +3,7 @@ add_executable(test_micm_util test_micm_util.F90) target_sources(test_micm_util PUBLIC - ${MUSICA_SRC_PATH}/micm/micm_util.F90 + ${MUSICA_SRC_PATH}/micm/musica_ccpp_micm_util.F90 ${CCPP_TEST_SRC_PATH}/ccpp_kinds.F90 ) diff --git a/test/musica/micm/test_micm_util.F90 b/test/musica/micm/test_micm_util.F90 index 76f2c13..19318dc 100644 --- a/test/musica/micm/test_micm_util.F90 +++ b/test/musica/micm/test_micm_util.F90 @@ -1,7 +1,7 @@ program test_micm_util use iso_c_binding - use micm_util + use musica_ccpp_micm_util use ccpp_kinds, only: kind_phys implicit none @@ -28,17 +28,14 @@ subroutine test_reshape() real(kind_phys), target :: pressure(NUM_COLUMNS,NUM_LAYERS) real(kind_phys), target :: dry_air_density(NUM_COLUMNS,NUM_LAYERS) real(kind_phys), target :: constituents(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) - real(kind_phys), target :: rate_params(NUM_COLUMNS,NUM_LAYERS,NUM_RATES) real(c_double), target :: m_temperature(NUM_GRID_CELLS) real(c_double), target :: m_pressure(NUM_GRID_CELLS) real(c_double), target :: m_dry_air_density(NUM_GRID_CELLS) real(c_double), target :: m_constituents(NUM_GRID_CELLS*NUM_SPECIES) - real(c_double), target :: m_rate_params(NUM_GRID_CELLS*NUM_RATES) ! local variables real(c_double), dimension(NUM_GRID_CELLS) :: arr_conditions real(c_double), dimension(NUM_GRID_CELLS*NUM_SPECIES) :: arr_constituents - real(c_double), dimension(NUM_GRID_CELLS*NUM_RATES) :: arr_rates integer :: i_column, i_layer, i_elem, i_arr real :: abs_error = 1e-7 @@ -52,17 +49,12 @@ subroutine test_reshape() constituents(1,2,:) = (/ 0.41_kind_phys, 0.42_kind_phys, 0.43_kind_phys, 0.44_kind_phys /) constituents(2,1,:) = (/ 0.21_kind_phys, 0.22_kind_phys, 0.23_kind_phys, 0.24_kind_phys /) constituents(2,2,:) = (/ 0.31_kind_phys, 0.32_kind_phys, 0.33_kind_phys, 0.34_kind_phys /) - rate_params(1,1,:) = (/ 100._kind_phys, 200._kind_phys, 300._kind_phys /) - rate_params(1,2,:) = (/ 400._kind_phys, 500._kind_phys, 600._kind_phys /) - rate_params(2,1,:) = (/ 700._kind_phys, 800._kind_phys, 900._kind_phys /) - rate_params(2,2,:) = (/ 1000._kind_phys, 1100._kind_phys, 1200._kind_phys /) arr_conditions = (/ 100.0, 200.0, 300.0, 400.0 /) arr_constituents = (/ 0.1, 0.2, 0.3, 0.4, 0.21, 0.22, 0.23, 0.24, 0.41, 0.42, 0.43, 0.44, 0.31, 0.32, 0.33, 0.34 /) - arr_rates = (/ 100., 200., 300., 700., 800., 900., 400., 500., 600., 1000., 1100., 1200. /) call reshape_into_micm_arr(temperature, pressure, dry_air_density, constituents, & - rate_params, m_temperature, m_pressure, m_dry_air_density, m_constituents, m_rate_params) + m_temperature, m_pressure, m_dry_air_density, m_constituents) do i_elem = 1, NUM_GRID_CELLS ASSERT(m_temperature(i_elem) == arr_conditions(i_elem)) @@ -74,22 +66,7 @@ subroutine test_reshape() ASSERT_NEAR(m_constituents(i_elem), arr_constituents(i_elem), abs_error) end do - do i_elem = 1, size(m_rate_params) - ASSERT_NEAR(m_rate_params(i_elem), arr_rates(i_elem), abs_error) - end do - - call reshape_into_ccpp_arr(temperature, pressure, dry_air_density, constituents, & - rate_params, m_temperature, m_pressure, m_dry_air_density, m_constituents, m_rate_params) - - i_elem = 1 - do i_layer = 1, NUM_LAYERS - do i_column = 1, NUM_COLUMNS - ASSERT(temperature(i_column, i_layer) == arr_conditions(i_elem)) - ASSERT(pressure(i_column, i_layer) == arr_conditions(i_elem)) - ASSERT(dry_air_density(i_column, i_layer) == arr_conditions(i_elem)) - i_elem = i_elem + 1 - end do - end do + call reshape_into_ccpp_arr(constituents, m_constituents) i_arr = 1 do i_layer = 1, NUM_LAYERS @@ -101,16 +78,6 @@ subroutine test_reshape() end do end do - i_arr = 1 - do i_layer = 1, NUM_LAYERS - do i_column = 1, NUM_COLUMNS - do i_elem = 1, NUM_RATES - ASSERT_NEAR(rate_params(i_column, i_layer, i_elem), arr_rates(i_arr), abs_error) - i_arr = i_arr + 1 - end do - end do - end do - end subroutine test_reshape subroutine test_unit_conversion() diff --git a/test/musica/test_musica_api.F90 b/test/musica/test_musica_api.F90 index e69f514..3fc25dc 100644 --- a/test/musica/test_musica_api.F90 +++ b/test/musica/test_musica_api.F90 @@ -18,13 +18,15 @@ subroutine test_musica_ccpp_api() integer :: solver_type integer :: errcode character(len=512) :: errmsg - real(kind_phys) :: time_step ! s - real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: height ! km - 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 - real(kind_phys), target, dimension(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) :: constituents ! kg kg-1 - real(kind_phys), target, dimension(NUM_COLUMNS,NUM_LAYERS,NUM_RATES) :: user_defined_reaction_rates + real(kind_phys) :: time_step ! s + 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 + real(kind_phys), target, dimension(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) :: constituents ! kg kg-1 type(ccpp_constituent_prop_ptr_t), allocatable :: constituent_props_ptr(:) ! local variables @@ -39,6 +41,12 @@ subroutine test_musica_ccpp_api() solver_type = Rosenbrock num_grid_cells = NUM_COLUMNS * NUM_LAYERS time_step = 60._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 /) @@ -49,11 +57,7 @@ subroutine test_musica_ccpp_api() constituents(1,2,:) = (/ 0.41_kind_phys, 0.42_kind_phys, 0.43_kind_phys, 0.44_kind_phys /) constituents(2,1,:) = (/ 0.21_kind_phys, 0.22_kind_phys, 0.23_kind_phys, 0.24_kind_phys /) constituents(2,2,:) = (/ 0.31_kind_phys, 0.32_kind_phys, 0.33_kind_phys, 0.34_kind_phys /) - user_defined_reaction_rates(1,1,:) = (/2.7e-19_kind_phys, 1.13e-9_kind_phys, 5.8e-8_kind_phys/) - user_defined_reaction_rates(1,2,:) = (/2.7e-19_kind_phys, 1.13e-9_kind_phys, 5.8e-8_kind_phys/) - user_defined_reaction_rates(2,1,:) = (/2.7e-19_kind_phys, 1.13e-9_kind_phys, 5.8e-8_kind_phys/) - user_defined_reaction_rates(2,2,:) = (/2.7e-19_kind_phys, 1.13e-9_kind_phys, 5.8e-8_kind_phys/) - + call musica_ccpp_register(constituent_props, solver_type, num_grid_cells, errmsg, errcode) ASSERT(allocated(constituent_props)) ASSERT(size(constituent_props) == NUM_SPECIES) @@ -66,7 +70,7 @@ subroutine test_musica_ccpp_api() ASSERT(errcode == 0) call constituent_props(i)%is_advected(is_advected, errcode, errmsg) ASSERT(errcode == 0) - tmp_booL = (trim(species_name) == "O2" .and. molar_mass == 0.032_kind_phys .and. .not. is_advected) .or. & + tmp_bool = (trim(species_name) == "O2" .and. molar_mass == 0.032_kind_phys .and. .not. is_advected) .or. & (trim(species_name) == "O" .and. molar_mass == 0.016_kind_phys .and. .not. is_advected) .or. & (trim(species_name) == "O1D" .and. molar_mass == 0.016_kind_phys .and. .not. is_advected) .or. & (trim(species_name) == "O3" .and. molar_mass == 0.048_kind_phys .and. is_advected) @@ -86,7 +90,7 @@ subroutine test_musica_ccpp_api() call constituent_props_ptr(i)%set(const_prop, errcode, errmsg) end do - call musica_ccpp_init(errmsg, errcode) + call musica_ccpp_init(NUM_LAYERS, NUM_LAYERS+1, errmsg, errcode) if (errcode /= 0) then write(*,*) trim(errmsg) stop 3 @@ -100,11 +104,12 @@ subroutine test_musica_ccpp_api() write(*,fmt="(4(1x,f10.4))") pressure write(*,*) "[MUSICA INFO] Initial Concentrations" write(*,fmt="(4(3x,e13.6))") constituents - write(*,*) "[MUSICA INFO] Initial User-defined Reaction Rates" - write(*,fmt="(3(3x,e13.6))") user_defined_reaction_rates call musica_ccpp_run(time_step, temperature, pressure, dry_air_density, constituent_props_ptr, & - constituents, user_defined_reaction_rates, height, 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 @@ -112,8 +117,6 @@ subroutine test_musica_ccpp_api() write(*,*) "[MUSICA INFO] Solved Concentrations" write(*,fmt="(4(3x,e13.6))") constituents - write(*,*) "[MUSICA INFO] Solved User-defined Reaction Rates" - write(*,fmt="(3(3x,e13.6))") user_defined_reaction_rates call musica_ccpp_final(errmsg, errcode) diff --git a/test/musica/tuvx/CMakeLists.txt b/test/musica/tuvx/CMakeLists.txt new file mode 100644 index 0000000..de5e4c9 --- /dev/null +++ b/test/musica/tuvx/CMakeLists.txt @@ -0,0 +1,27 @@ + +add_executable(test_tuvx_height_grid test_tuvx_height_grid.F90) + +target_sources(test_tuvx_height_grid + PUBLIC + ${MUSICA_SRC_PATH}/tuvx/musica_ccpp_tuvx_height_grid.F90 + ${MUSICA_SRC_PATH}/musica_ccpp_util.F90 + ${CCPP_TEST_SRC_PATH}/ccpp_kinds.F90 +) + +target_link_libraries(test_tuvx_height_grid + PRIVATE + musica::musica-fortran +) + +set_target_properties(test_tuvx_height_grid + PROPERTIES + LINKER_LANGUAGE Fortran +) + +add_test( + NAME test_tuvx_height_grid + COMMAND $ + WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} +) + +add_memory_check_test(test_tuvx_height_grid $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) diff --git a/test/musica/tuvx/configs/ts1_tsmlt.json b/test/musica/tuvx/configs/ts1_tsmlt.json index 7ea1806..817ae96 100644 --- a/test/musica/tuvx/configs/ts1_tsmlt.json +++ b/test/musica/tuvx/configs/ts1_tsmlt.json @@ -4,14 +4,6 @@ "cross section parameters file": "data/cross_sections/O2_parameters.txt" }, "grids": [ - { - "name": "height", - "type": "equal interval", - "units": "km", - "begins at" : 0.0, - "ends at" : 120.0, - "cell delta" : 1.0 - }, { "name": "wavelength", "type": "from csv file", diff --git a/test/musica/tuvx/test_tuvx_height_grid.F90 b/test/musica/tuvx/test_tuvx_height_grid.F90 new file mode 100644 index 0000000..19339c0 --- /dev/null +++ b/test/musica/tuvx/test_tuvx_height_grid.F90 @@ -0,0 +1,107 @@ +program test_tuvx_height_grid + + use iso_c_binding + use musica_ccpp_tuvx_height_grid + use ccpp_kinds, only: kind_phys + + 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 + + call test_create_height_grid() + call test_calculate_height_grid_values() + +contains + + subroutine test_create_height_grid() + + use musica_util, only: error_t + use musica_tuvx_grid, only: grid_t + use iso_c_binding, only: c_double + use ccpp_kinds, only: kind_phys + + integer, parameter :: NUM_HOST_MIDPOINTS = 2 + integer, parameter :: NUM_HOST_EDGES = 3 + real(kind_phys), target :: host_midpoints(NUM_HOST_MIDPOINTS) = [200.8_kind_phys, 100.6_kind_phys] + real(kind_phys), target :: host_edges(NUM_HOST_EDGES) = [250.3_kind_phys, 150.2_kind_phys, 50.1_kind_phys] + type(grid_t), pointer :: height_grid + character(len=512) :: errmsg + integer :: errcode + real(kind_phys) :: abs_error = 1e-5 + integer :: i + + ! local variables + real(kind_phys), dimension(NUM_HOST_MIDPOINTS+1) :: midpoints + real(kind_phys), dimension(NUM_HOST_EDGES+1) :: edges + type(error_t) :: error + + height_grid => create_height_grid(-1, 0, errmsg, errcode) + ASSERT(errcode == 1) + ASSERT(.not. associated(height_grid)) + + height_grid => create_height_grid(3, 3, errmsg, errcode) + ASSERT(errcode == 1) + ASSERT(.not. associated(height_grid)) + + height_grid => create_height_grid(NUM_HOST_MIDPOINTS, NUM_HOST_EDGES, & + errmsg, errcode) + ASSERT(errcode == 0) + ASSERT(associated(height_grid)) + + call set_height_grid_values(height_grid, host_midpoints, host_edges, & + errmsg, errcode) + ASSERT(errcode == 0) + + ASSERT(height_grid%number_of_sections(error) == NUM_HOST_MIDPOINTS + 1) + ASSERT(error%is_success()) + + call height_grid%get_midpoints(midpoints, error) + ASSERT(error%is_success()) + + call height_grid%get_edges(edges, error) + ASSERT(error%is_success()) + + ASSERT_NEAR(midpoints(1), (100.6 + 50.1) * 0.5, abs_error) + ASSERT_NEAR(midpoints(2), 150.2, abs_error) + ASSERT_NEAR(midpoints(3), (250.3 + 200.8) * 0.5, abs_error) + ASSERT_NEAR(edges(1), 50.1, abs_error) + ASSERT_NEAR(edges(2), 100.6, abs_error) + ASSERT_NEAR(edges(3), 200.8, abs_error) + ASSERT_NEAR(edges(4), 250.3, abs_error) + + deallocate( 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