From 8072cd9d03872fa553bbb4e815112dd28527ccbb Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Fri, 22 Nov 2024 14:39:34 -0800 Subject: [PATCH 01/14] draft cloud optics calculations --- schemes/musica/musica_ccpp.F90 | 38 ++-- schemes/musica/musica_ccpp.meta | 18 ++ schemes/musica/tuvx/musica_ccpp_tuvx.F90 | 176 ++++++++++++++---- .../tuvx/musica_ccpp_tuvx_cloud_optics.F90 | 129 +++++++++++++ schemes/musica/tuvx/musica_ccpp_tuvx_util.F90 | 11 +- test/CMakeLists.txt | 1 + test/musica/CMakeLists.txt | 1 + test/musica/test_musica_api.F90 | 62 ++++-- test/musica/tuvx/CMakeLists.txt | 32 +++- test/musica/tuvx/test_tuvx_cloud_optics.F90 | 115 ++++++++++++ to_be_ccppized/ccpp_const_utils.F90 | 10 +- 11 files changed, 524 insertions(+), 69 deletions(-) create mode 100644 schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 create mode 100644 test/musica/tuvx/test_tuvx_cloud_optics.F90 diff --git a/schemes/musica/musica_ccpp.F90 b/schemes/musica/musica_ccpp.F90 index d3512f4a..ed676bdc 100644 --- a/schemes/musica/musica_ccpp.F90 +++ b/schemes/musica/musica_ccpp.F90 @@ -2,7 +2,7 @@ module musica_ccpp use musica_ccpp_micm, only: micm_register, micm_init, micm_run, micm_final use musica_ccpp_namelist, only: filename_of_tuvx_micm_mapping_configuration - use musica_ccpp_tuvx, only: tuvx_init, tuvx_run, tuvx_final + use musica_ccpp_tuvx, only: tuvx_register, tuvx_init, tuvx_run, tuvx_final use musica_util, only: index_mappings_t implicit none @@ -23,28 +23,40 @@ subroutine musica_ccpp_register(solver_type, num_grid_cells, constituent_props, character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode - call micm_register(solver_type, num_grid_cells, constituent_props, errmsg, errcode) + type(ccpp_constituent_properties_t), allocatable :: constituent_set(:) + + call micm_register(solver_type, num_grid_cells, constituent_set, errmsg, errcode) + if (errcode /= 0) return + constituent_props = constituent_set + + call tuvx_register(constituent_set, errmsg, errcode) + if (errcode /= 0) return + constituent_props = [ constituent_props, constituent_set ] end subroutine musica_ccpp_register !> \section arg_table_musica_ccpp_init Argument Table !! \htmlinclude musica_ccpp_init.html subroutine musica_ccpp_init(vertical_layer_dimension, vertical_interface_dimension, & - photolysis_wavelength_grid_interfaces, errmsg, errcode) + photolysis_wavelength_grid_interfaces, & + constituent_props, errmsg, errcode) + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t use ccpp_kinds, only : kind_phys use musica_ccpp_micm, only: micm use musica_ccpp_util, only: has_error_occurred - integer, intent(in) :: vertical_layer_dimension ! (count) - integer, intent(in) :: vertical_interface_dimension ! (count) - real(kind_phys), intent(in) :: photolysis_wavelength_grid_interfaces(:) ! m - character(len=512), intent(out) :: errmsg - integer, intent(out) :: errcode + integer, intent(in) :: vertical_layer_dimension ! (count) + integer, intent(in) :: vertical_interface_dimension ! (count) + real(kind_phys), intent(in) :: photolysis_wavelength_grid_interfaces(:) ! m + type(ccpp_constituent_prop_ptr_t), intent(in) :: constituent_props(:) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode call micm_init(errmsg, errcode) if (errcode /= 0) return call tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & photolysis_wavelength_grid_interfaces, & - micm%user_defined_reaction_rates, errmsg, errcode) + micm%user_defined_reaction_rates, & + constituent_props, errmsg, errcode) if (errcode /= 0) return end subroutine musica_ccpp_init @@ -59,7 +71,8 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co constituents, geopotential_height_wrt_surface_at_midpoint, & geopotential_height_wrt_surface_at_interface, surface_temperature, & surface_geopotential, surface_albedo, & - standard_gravitational_acceleration, errmsg, errcode) + standard_gravitational_acceleration, cloud_area_fraction, & + air_pressure_thickness, 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 @@ -78,6 +91,8 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co real(kind_phys), intent(in) :: surface_geopotential(:) ! m2 s-2 real(kind_phys), intent(in) :: surface_albedo ! unitless 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) character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode @@ -95,7 +110,8 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co surface_temperature, surface_geopotential, & surface_albedo, & standard_gravitational_acceleration, & - rate_parameters, & + cloud_area_fraction, constituents, & + air_pressure_thickness, rate_parameters, & 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 85e01f3e..5fcfcb32 100644 --- a/schemes/musica/musica_ccpp.meta +++ b/schemes/musica/musica_ccpp.meta @@ -59,6 +59,12 @@ type = real | kind = kind_phys dimensions = (photolysis_wavelength_grid_interface_dimension) intent = in +[ constituent_props ] + standard_name = ccpp_constituent_properties + units = None + type = ccpp_constituent_prop_ptr_t + dimensions = (number_of_ccpp_constituents) + intent = in [ errmsg ] standard_name = ccpp_error_message units = none @@ -147,6 +153,18 @@ type = real | kind = kind_phys dimensions = () intent = in +[ cloud_area_fraction ] + standard_name = cloud_area_fraction + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + intent = in +[ air_pressure_thickness ] + standard_name = air_pressure_thickness + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + 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 5a12a3ee..e6384e30 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 @@ -5,29 +5,73 @@ module musica_ccpp_tuvx use ccpp_kinds, only: kind_phys use musica_ccpp_namelist, only: filename_of_tuvx_configuration use musica_ccpp_util, only: has_error_occurred - use musica_tuvx, only: tuvx_t, grid_t, profile_t + use musica_tuvx, only: tuvx_t, grid_t, profile_t, radiator_t use musica_util, only: mappings_t, index_mappings_t implicit none private - public :: tuvx_init, tuvx_run, tuvx_final + public :: tuvx_register, tuvx_init, tuvx_run, tuvx_final - type(tuvx_t), pointer :: tuvx => null() - type(grid_t), pointer :: height_grid => null() - type(grid_t), pointer :: wavelength_grid => null() - type(profile_t), pointer :: temperature_profile => null() - type(profile_t), pointer :: surface_albedo_profile => null() + type(tuvx_t), pointer :: tuvx => null() + type(grid_t), pointer :: height_grid => null() + type(grid_t), pointer :: wavelength_grid => null() + type(profile_t), pointer :: temperature_profile => null() + type(profile_t), pointer :: surface_albedo_profile => null() + type(radiator_t), pointer :: cloud_optics => null() type(index_mappings_t), pointer :: photolysis_rate_constants_mapping => null( ) - integer :: number_of_photolysis_rate_constants = 0 + integer, parameter :: DEFAULT_NUM_PHOTOLYSIS_RATE_CONSTANTS = 0 + integer :: number_of_photolysis_rate_constants = DEFAULT_NUM_PHOTOLYSIS_RATE_CONSTANTS + integer, parameter :: DEFAULT_INDEX_NOT_FOUND = -1 + character(len=*), parameter :: CLOUD_LIQUID_WATER_CONTENT_LABEL = & + 'cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water' + character(len=*), parameter :: CLOUD_LIQUID_WATER_CONTENT_LONG_NAME = & + 'Cloud water mass mixing ratio with respect to moist air plus all airborne condensates' + character(len=*), parameter :: CLOUD_LIQUID_WATER_CONTENT_UNITS = 'kg kg-1' + real(kind_phys), parameter :: CLOUD_LIQUID_WATER_CONTENT_MOLAR_MASS = 0.018_kind_phys ! kg mol-1 + integer :: index_cloud_liquid_water_content = DEFAULT_INDEX_NOT_FOUND contains + !> Registers constituent properties with the CCPP needed by TUV-x + subroutine tuvx_register(constituent_props, errmsg, errcode) + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + use musica_util, only: error_t + + type(ccpp_constituent_properties_t), allocatable, intent(out) :: constituent_props(:) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode + + allocate(constituent_props(1), stat=errcode) + if (errcode /= 0) then + errmsg = "[MUSICA Error] Failed to allocate memory for constituent properties." + return + end if + + ! Register cloud liquid water content needed for cloud optics calculations + call constituent_props(1)%instantiate( & + std_name = CLOUD_LIQUID_WATER_CONTENT_LABEL, & + long_name = CLOUD_LIQUID_WATER_CONTENT_LONG_NAME, & + units = CLOUD_LIQUID_WATER_CONTENT_UNITS, & + vertical_dim = "vertical_layer_dimension", & + default_value = 0.0_kind_phys, & + min_value = 0.0_kind_phys, & + molar_mass = CLOUD_LIQUID_WATER_CONTENT_MOLAR_MASS, & + advected = .true., & + errcode = errcode, & + errmsg = errmsg & + ) + if (errcode /= 0) return + + end subroutine tuvx_register + !> Initializes TUV-x subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & wavelength_grid_interfaces, micm_rate_parameter_ordering, & - errmsg, errcode) + constituent_props, errmsg, errcode) + use ccpp_const_utils, only: ccpp_const_get_idx + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t 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 @@ -40,13 +84,16 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & only: create_temperature_profile, temperature_label, temperature_unit use musica_ccpp_tuvx_surface_albedo, & only: create_surface_albedo_profile, surface_albedo_label, surface_albedo_unit + use musica_ccpp_tuvx_cloud_optics, & + only: create_cloud_optics_radiator, cloud_optics_label - integer, intent(in) :: vertical_layer_dimension ! (count) - integer, intent(in) :: vertical_interface_dimension ! (count) - real(kind_phys), intent(in) :: wavelength_grid_interfaces(:) ! m - type(mappings_t), intent(in) :: micm_rate_parameter_ordering ! index mappings for MICM rate parameters - character(len=512), intent(out) :: errmsg - integer, intent(out) :: errcode + integer, intent(in) :: vertical_layer_dimension ! (count) + integer, intent(in) :: vertical_interface_dimension ! (count) + real(kind_phys), intent(in) :: wavelength_grid_interfaces(:) ! m + type(mappings_t), intent(in) :: micm_rate_parameter_ordering ! index mappings for MICM rate parameters + type(ccpp_constituent_prop_ptr_t), intent(in) :: constituent_props(:) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode ! local variables type(grid_map_t), pointer :: grids @@ -56,6 +103,11 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & type(mappings_t), pointer :: photolysis_rate_constants_ordering type(error_t) :: error + ! Get needed indices in constituents array + call ccpp_const_get_idx(constituent_props, CLOUD_LIQUID_WATER_CONTENT_LABEL, & + index_cloud_liquid_water_content, errmsg, errcode) + if (errcode /= 0) return + grids => grid_map_t( error ) if (has_error_occurred( error, errmsg, errcode )) return @@ -69,7 +121,7 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & call grids%add( height_grid, error ) if (has_error_occurred( error, errmsg, errcode )) then call tuvx_deallocate( grids, null(), null(), null(), height_grid, null(), & - null(), null() ) + null(), null(), null() ) return end if @@ -77,35 +129,35 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & errmsg, errcode ) if (errcode /= 0) then call tuvx_deallocate( grids, null(), null(), null(), height_grid, null(), & - null(), null() ) + null(), null(), null() ) return endif call grids%add( wavelength_grid, error ) if (has_error_occurred( error, errmsg, errcode )) then call tuvx_deallocate( grids, null(), null(), null(), height_grid, & - wavelength_grid, null(), null() ) + wavelength_grid, null(), null(), null() ) return end if profiles => profile_map_t( error ) if (has_error_occurred( error, errmsg, errcode )) then call tuvx_deallocate( grids, null(), null(), null(), height_grid, & - wavelength_grid, null(), null() ) + wavelength_grid, null(), null(), null() ) return end if temperature_profile => create_temperature_profile( height_grid, errmsg, errcode ) if (errcode /= 0) then call tuvx_deallocate( grids, profiles, null(), null(), height_grid, & - wavelength_grid, null(), null() ) + wavelength_grid, null(), null(), null() ) return endif call profiles%add( temperature_profile, error ) if (has_error_occurred( error, errmsg, errcode )) then call tuvx_deallocate( grids, profiles, null(), null(), height_grid, & - wavelength_grid, temperature_profile, null() ) + wavelength_grid, temperature_profile, null(), null() ) return end if @@ -113,21 +165,37 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & errmsg, errcode ) if (errcode /= 0) then call tuvx_deallocate( grids, profiles, null(), null(), height_grid, & - wavelength_grid, temperature_profile, null() ) + wavelength_grid, temperature_profile, null(), null() ) return endif call profiles%add( surface_albedo_profile, error ) if (has_error_occurred( error, errmsg, errcode )) then call tuvx_deallocate( grids, profiles, null(), null(), height_grid, & - wavelength_grid, temperature_profile, surface_albedo_profile ) + wavelength_grid, temperature_profile, surface_albedo_profile, null() ) return end if radiators => radiator_map_t( error ) if (has_error_occurred( error, errmsg, errcode )) then call tuvx_deallocate( grids, profiles, null(), null(), height_grid, & - wavelength_grid, temperature_profile, surface_albedo_profile ) + wavelength_grid, temperature_profile, surface_albedo_profile, null() ) + return + end if + + cloud_optics => create_cloud_optics_radiator( height_grid, wavelength_grid, & + errmsg, errcode ) + if (errcode /= 0) then + call tuvx_deallocate( grids, profiles, radiators, null(), height_grid, & + wavelength_grid, temperature_profile, surface_albedo_profile, null() ) + return + endif + + call radiators%add( cloud_optics, error ) + if (has_error_occurred( error, errmsg, errcode )) then + call tuvx_deallocate( grids, profiles, radiators, null(), height_grid, & + wavelength_grid, temperature_profile, surface_albedo_profile, & + cloud_optics ) return end if @@ -135,12 +203,12 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & radiators, error ) if (has_error_occurred( error, errmsg, errcode )) then call tuvx_deallocate( grids, profiles, radiators, null(), height_grid, & - wavelength_grid, temperature_profile, surface_albedo_profile ) + wavelength_grid, temperature_profile, surface_albedo_profile, null() ) return end if call tuvx_deallocate( grids, profiles, radiators, null(), height_grid, & - wavelength_grid, temperature_profile, surface_albedo_profile ) + wavelength_grid, temperature_profile, surface_albedo_profile, cloud_optics ) grids => tuvx%get_grids( error ) if (has_error_occurred( error, errmsg, errcode )) return @@ -148,7 +216,7 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & height_grid => grids%get( height_grid_label, height_grid_unit, error ) if (has_error_occurred( error, errmsg, errcode )) then call tuvx_deallocate( grids, null(), null(), tuvx, null(), null(), & - null(), null() ) + null(), null(), null() ) return end if @@ -156,33 +224,47 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & error ) if (has_error_occurred( error, errmsg, errcode )) then call tuvx_deallocate( grids, null(), null(), tuvx, height_grid, null(), & - null(), null() ) + null(), null(), null() ) return end if profiles => tuvx%get_profiles( error ) if (has_error_occurred( error, errmsg, errcode )) then call tuvx_deallocate( grids, null(), null(), tuvx, height_grid, & - wavelength_grid, null(), null() ) + wavelength_grid, null(), null(), null() ) return end if temperature_profile => profiles%get( temperature_label, temperature_unit, error ) if (has_error_occurred( error, errmsg, errcode )) then call tuvx_deallocate( grids, profiles, null(), tuvx, height_grid, & - wavelength_grid, null(), null() ) + wavelength_grid, null(), null(), null() ) return end if surface_albedo_profile => profiles%get( surface_albedo_label, surface_albedo_unit, error ) if (has_error_occurred( error, errmsg, errcode )) then call tuvx_deallocate( grids, profiles, null(), tuvx, height_grid, & - wavelength_grid, temperature_profile, null() ) + wavelength_grid, temperature_profile, null(), null() ) + return + end if + + radiators => tuvx%get_radiators( error ) + if (has_error_occurred( error, errmsg, errcode )) then + call tuvx_deallocate( grids, profiles, null(), tuvx, height_grid, & + wavelength_grid, temperature_profile, surface_albedo_profile, null() ) + return + end if + + cloud_optics => radiators%get( cloud_optics_label, error ) + if (has_error_occurred( error, errmsg, errcode )) then + call tuvx_deallocate( grids, profiles, radiators, tuvx, height_grid, & + wavelength_grid, temperature_profile, surface_albedo_profile, null() ) return end if - call tuvx_deallocate( grids, profiles, null(), null(), null(), null(), & - null(), null() ) + call tuvx_deallocate( grids, profiles, radiators, null(), null(), null(), & + null(), null(), null() ) photolysis_rate_constants_ordering => & tuvx%get_photolysis_rate_constants_ordering( error ) @@ -212,12 +294,15 @@ subroutine tuvx_run(temperature, dry_air_density, & surface_temperature, surface_geopotential, & surface_albedo, & standard_gravitational_acceleration, & - 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 + 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_tuvx_surface_albedo, only: set_surface_albedo_values + use musica_ccpp_tuvx_cloud_optics, only: set_cloud_optics_values real(kind_phys), intent(in) :: temperature(:,:) ! K (column, layer) real(kind_phys), intent(in) :: dry_air_density(:,:) ! kg m-3 (column, layer) @@ -227,6 +312,9 @@ subroutine tuvx_run(temperature, dry_air_density, & real(kind_phys), intent(in) :: surface_geopotential(:) ! m2 s-2 real(kind_phys), intent(in) :: surface_albedo ! unitless real(kind_phys), intent(in) :: standard_gravitational_acceleration ! m s-2 + real(kind_phys), intent(in) :: cloud_area_fraction(:,:) ! unitless (column, layer) + real(kind_phys), intent(in) :: constituents(:,:,:) ! various (column, layer, constituent) + real(kind_phys), intent(in) :: air_pressure_thickness(:,:) ! Pa (column, layer) real(kind_phys), intent(inout) :: rate_parameters(:,:,:) ! various units (column, layer, reaction) character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode @@ -263,6 +351,13 @@ subroutine tuvx_run(temperature, dry_air_density, & surface_temperature(i_col), errmsg, errcode ) if (errcode /= 0) return + call set_cloud_optics_values( cloud_optics, cloud_area_fraction(i_col,:), & + air_pressure_thickness(i_col,:), & + constituents(i_col,:,index_cloud_liquid_water_content), & + reciprocal_of_gravitational_acceleration, & + errmsg, errcode ) + if (errcode /= 0) return + ! temporary values until these are available from the host model solar_zenith_angle = 0.0_kind_phys earth_sun_distance = 1.0_kind_phys @@ -315,6 +410,11 @@ subroutine tuvx_final(errmsg, errcode) surface_albedo_profile => null() end if + if (associated( cloud_optics )) then + deallocate( cloud_optics ) + cloud_optics => null() + end if + if (associated( tuvx )) then deallocate( tuvx ) tuvx => null() diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 new file mode 100644 index 00000000..a6459e3d --- /dev/null +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 @@ -0,0 +1,129 @@ +! Copyright (C) 2024 National Science Foundation-National Center for Atmospheric Research +! SPDX-License-Identifier: Apache-2.0 +module musica_ccpp_tuvx_cloud_optics + implicit none + + private + public :: create_cloud_optics_radiator, set_cloud_optics_values, & + cloud_optics_label + + ! This module is used to set the optical properties of clouds in TUV-x. + ! Optical properties are defined as a function of wavelength and height, + ! and include the cloud optical depth, single scattering albedo, + ! and asymmetry parameter. + ! + ! See module_ccpp_tuvx_height_grid for the definition of the height grid + ! and its mapping to the CAM-SIMA vertical grid. + + !> Label for cloud optical properties in TUV-x + character(len=*), parameter :: cloud_optics_label = "clouds" + !> Default value of number of vertical levels + integer, parameter :: DEFAULT_NUM_VERTICAL_LEVELS = 0 + !> Number of vertical levels + integer, protected :: num_vertical_levels = DEFAULT_NUM_VERTICAL_LEVELS + !> Default value of number of wavelength bins + integer, parameter :: DEFAULT_NUM_WAVELENGTH_BINS = 0 + !> Number of wavelength bins + integer, protected :: num_wavelength_bins = DEFAULT_NUM_WAVELENGTH_BINS + +contains + + !> Creates a TUV-x cloud optics radiator from the host-model wavelength grid + function create_cloud_optics_radiator( height_grid, wavelength_grid, & + errmsg, errcode ) result( radiator ) + use musica_ccpp_util, only: has_error_occurred + use musica_tuvx_grid, only: grid_t + use musica_tuvx_radiator, only: radiator_t + use musica_util, only: error_t + + type(grid_t), intent(inout) :: height_grid + type(grid_t), intent(inout) :: wavelength_grid + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errcode + type(radiator_t), pointer :: radiator + + ! local variables + type(error_t) :: error + + num_vertical_levels = height_grid%number_of_sections( error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + + num_wavelength_bins = wavelength_grid%number_of_sections( error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + + radiator => radiator_t( cloud_optics_label, height_grid, wavelength_grid, & + error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + + end function create_cloud_optics_radiator + + !> Sets TUV-x cloud optics values + subroutine set_cloud_optics_values( radiator, cloud_fraction, delta_pressure, & + cloud_liquid_water_content, & + reciprocal_of_gravitational_acceleration, & + errmsg, errcode ) + use ccpp_kinds, only: kind_phys + use musica_ccpp_util, only: has_error_occurred + use musica_tuvx_radiator, only: radiator_t + use musica_util, only: error_t + + type(radiator_t), intent(inout) :: radiator + real(kind_phys), intent(in) :: delta_pressure(:) ! pressure delta about vertical level midpoints (Pa) + real(kind_phys), intent(in) :: cloud_fraction(:) ! (unitless) + real(kind_phys), intent(in) :: cloud_liquid_water_content(:) ! (kg/kg) + real(kind_phys), intent(in) :: reciprocal_of_gravitational_acceleration ! (s^2/m) + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errcode + + ! local variables + type(error_t) :: error + real(kind_phys) :: od(num_vertical_levels) + real(kind_phys) :: cloud_optical_depth(num_vertical_levels, num_wavelength_bins) + integer :: n_host_midpoints + integer :: i_level + + n_host_midpoints = size(cloud_fraction) + if ( size(delta_pressure) /= n_host_midpoints ) then + errmsg = "[MUSICA Error] Invalid size of cloud pressure delta for TUV-x." + errcode = 1 + return + end if + if ( size(cloud_liquid_water_content) /= n_host_midpoints ) then + errmsg = "[MUSICA Error] Invalid size of cloud liquid water content for TUV-x." + errcode = 1 + return + end if + if ( n_host_midpoints + 1 /= num_vertical_levels ) then + errmsg = "[MUSICA Error] Invalid size of cloud fraction for TUV-x." + errcode = 1 + return + end if + + ! Estimate cloud optical depth (od) [unitless] from cloud fraction (cf) + ! [unitless] and liquid water content (lwc) [kg kg-1] by first calculating + ! the cloud liquid water path (lwp) [kg m-2]: + ! lwp = 1/g * lwc * dP / cf + ! where g is the gravitational acceleration [m s-2] and dP is the change in + ! pressure across the vertical level [Pa]. + ! The cloud optical depth is then estimated as: + ! od = lwp * 155 * cf^1.5 + ! A constant cloud optical depth is used for all wavelengths. + do i_level = 1, n_host_midpoints + if ( cloud_fraction(i_level) > 0.0_kind_phys ) then + od(i_level) = ( reciprocal_of_gravitational_acceleration & + * cloud_liquid_water_content(i_level) * delta_pressure(i_level) & + / cloud_fraction(i_level) ) * 155.0_kind_phys * cloud_fraction(i_level)**1.5_kind_phys + else + od(i_level) = 0.0_kind_phys + end if + end do + do i_level = 1, n_host_midpoints + cloud_optical_depth(i_level, :) = od(n_host_midpoints-i_level+1) + end do + cloud_optical_depth(num_vertical_levels, :) = 0.0_kind_phys + call radiator%set_optical_depths( cloud_optical_depth, error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + + end subroutine set_cloud_optics_values + +end module musica_ccpp_tuvx_cloud_optics diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_util.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_util.F90 index 25ac1a28..4bcf814f 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_util.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_util.F90 @@ -8,9 +8,10 @@ module musica_ccpp_tuvx_util !> This is a helper subroutine created to deallocate objects associated with TUV-x subroutine tuvx_deallocate(grids, profiles, radiators, tuvx, height_grid, & - wavelength_grid, temperature_profile, surface_albedo_profile) + wavelength_grid, temperature_profile, surface_albedo_profile, & + cloud_optics) use musica_tuvx, only: tuvx_t, grid_map_t, profile_map_t, radiator_map_t, & - grid_t, profile_t + grid_t, profile_t, radiator_t type(grid_map_t), pointer :: grids type(profile_map_t), pointer :: profiles @@ -20,6 +21,7 @@ subroutine tuvx_deallocate(grids, profiles, radiators, tuvx, height_grid, & type(grid_t), pointer :: wavelength_grid type(profile_t), pointer :: temperature_profile type(profile_t), pointer :: surface_albedo_profile + type(radiator_t), pointer :: cloud_optics if (associated( grids )) deallocate( grids ) if (associated( profiles )) deallocate( profiles ) @@ -50,6 +52,11 @@ subroutine tuvx_deallocate(grids, profiles, radiators, tuvx, height_grid, & surface_albedo_profile => null() end if + if (associated( cloud_optics )) then + deallocate( cloud_optics ) + cloud_optics => null() + end if + end subroutine tuvx_deallocate end module musica_ccpp_tuvx_util \ No newline at end of file diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 6a7d1e30..3d541d63 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -26,6 +26,7 @@ set(BUILD_GMOCK OFF CACHE BOOL "" FORCE) if (CCPP_ENABLE_MUSICA_TESTS) set(MUSICA_VERSION $ENV{MUSICA_VERSION}) set(MUSICA_SRC_PATH ${CMAKE_SOURCE_DIR}/../schemes/musica) + set(TO_BE_CCPPIZED_PATH ${CMAKE_SOURCE_DIR}/../to_be_ccppized) set(CCPP_SRC_PATH ${CMAKE_SOURCE_DIR}/$ENV{CCPP_SRC_PATH}) set(CCPP_TEST_SRC_PATH ${CMAKE_SOURCE_DIR}/include) diff --git a/test/musica/CMakeLists.txt b/test/musica/CMakeLists.txt index 74b29693..ceac05dc 100644 --- a/test/musica/CMakeLists.txt +++ b/test/musica/CMakeLists.txt @@ -33,6 +33,7 @@ target_sources(test_musica_api ${CCPP_SRC_PATH}/ccpp_hash_table.F90 ${CCPP_SRC_PATH}/ccpp_hashable.F90 ${CCPP_SRC_PATH}/ccpp_types.F90 + ${TO_BE_CCPPIZED_PATH}/ccpp_const_utils.F90 ${CCPP_TEST_SRC_PATH}/ccpp_kinds.F90 ) diff --git a/test/musica/test_musica_api.F90 b/test/musica/test_musica_api.F90 index 61cca2f7..68be9237 100644 --- a/test/musica/test_musica_api.F90 +++ b/test/musica/test_musica_api.F90 @@ -145,6 +145,7 @@ subroutine test_chapman() implicit none integer, parameter :: NUM_SPECIES = 5 + integer, parameter :: NUM_TUVX_CONSTITUENTS = 1 ! This test requires that the number of grid cells = 4, which is the default ! vector dimension for MICM. This restriction will be removed once ! https://github.com/NCAR/musica/issues/217 is finished. @@ -166,8 +167,12 @@ subroutine test_chapman() real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: temperature ! K real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: pressure ! Pa real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: dry_air_density ! kg m-3 - real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) :: constituents ! kg kg-1 - real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) :: initial_constituents ! kg kg-1 + real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: cloud_area_fraction ! unitless + real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: air_pressure_thickness ! Pa + real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS, & + 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 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 @@ -196,6 +201,10 @@ subroutine test_chapman() pressure(:,2) = (/ 8000.04_kind_phys, 9000.04_kind_phys /) dry_air_density(:,1) = (/ 3.5_kind_phys, 4.5_kind_phys /) dry_air_density(:,2) = (/ 5.5_kind_phys, 6.5_kind_phys /) + cloud_area_fraction(:,1) = (/ 0.1_kind_phys, 0.2_kind_phys /) + 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 /) filename_of_micm_configuration = 'musica_configurations/chapman/micm/config.json' filename_of_tuvx_configuration = 'musica_configurations/chapman/tuvx/config.json' @@ -207,7 +216,7 @@ subroutine test_chapman() stop 3 endif ASSERT(allocated(constituent_props)) - ASSERT(size(constituent_props) == NUM_SPECIES) + ASSERT(size(constituent_props) == NUM_SPECIES+NUM_TUVX_CONSTITUENTS) do i = 1, size(constituent_props) ASSERT(constituent_props(i)%is_instantiated(errcode, errmsg)) ASSERT(errcode == 0) @@ -221,7 +230,9 @@ subroutine test_chapman() (trim(species_name) == "O" .and. molar_mass == 0.0159994_kind_phys .and. .not. is_advected) .or. & (trim(species_name) == "O1D" .and. molar_mass == 0.0159994_kind_phys .and. .not. is_advected) .or. & (trim(species_name) == "O3" .and. molar_mass == 0.0479982_kind_phys .and. is_advected) .or. & - (trim(species_name) == "N2" .and. molar_mass == 0.0280134_kind_phys .and. is_advected) + (trim(species_name) == "N2" .and. molar_mass == 0.0280134_kind_phys .and. is_advected) .or. & + (trim(species_name) == "cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water" .and. & + molar_mass == 0.018_kind_phys .and. is_advected) ASSERT(tmp_bool) call constituent_props(i)%units(units, errcode, errmsg) if (errcode /= 0) then @@ -242,7 +253,7 @@ subroutine test_chapman() end do call musica_ccpp_init(NUM_LAYERS, NUM_LAYERS+1, photolysis_wavelength_grid_interfaces, & - errmsg, errcode) + constituent_props_ptr, errmsg, errcode) if (errcode /= 0) then write(*,*) trim(errmsg) stop 3 @@ -275,6 +286,12 @@ subroutine test_chapman() end do end do end do + ! set initial cloud liquid water mixing ratio to ~1e-3 kg kg-1 + do j = 1, NUM_COLUMNS + do k = 1, NUM_LAYERS + constituents(j,k,NUM_SPECIES+1) = 1.0e-3_kind_phys * (1.0 + 0.1 * (j-1) + 0.01 * (k-1)) + end do + end do initial_constituents(:,:,:) = constituents(:,:,:) write(*,*) "[MUSICA INFO] Initial Time Step" @@ -290,7 +307,7 @@ subroutine test_chapman() constituents, geopotential_height_wrt_surface_at_midpoint, & geopotential_height_wrt_surface_at_interface, surface_temperature, & surface_geopotential, surface_albedo, standard_gravitational_acceleration, & - errmsg, errcode) + cloud_area_fraction, air_pressure_thickness, errmsg, errcode) if (errcode /= 0) then write(*,*) trim(errmsg) stop 3 @@ -322,6 +339,8 @@ subroutine test_chapman() constituents(i,j,O2_index) + constituents(i,j,O3_index) total_O_init = initial_constituents(i,j,O_index) + initial_constituents(i,j,O1D_index) + & initial_constituents(i,j,O2_index) + initial_constituents(i,j,O3_index) + ! cloud liquid water mixing ratio should be unchanged + ASSERT_NEAR(constituents(i,j,NUM_SPECIES+1), initial_constituents(i,j,NUM_SPECIES+1), 1.0e-13) ASSERT_NEAR(total_O, total_O_init, 1.0e-13) end do end do @@ -344,6 +363,7 @@ subroutine test_terminator() implicit none integer, parameter :: NUM_SPECIES = 2 + integer, parameter :: NUM_TUVX_CONSTITUENTS = 1 ! This test requires that the number of grid cells = 4, which is the default ! vector dimension for MICM. This restriction will be removed once ! https://github.com/NCAR/musica/issues/217 is finished. @@ -365,8 +385,12 @@ subroutine test_terminator() real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: temperature ! K real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: pressure ! Pa real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: dry_air_density ! kg m-3 - real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) :: constituents ! kg kg-1 - real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) :: initial_constituents ! kg kg-1 + real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: cloud_area_fraction ! unitless + real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: air_pressure_thickness ! Pa + real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS, & + 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 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 @@ -395,6 +419,10 @@ subroutine test_terminator() pressure(:,2) = (/ 8000.04_kind_phys, 9000.04_kind_phys /) dry_air_density(:,1) = (/ 3.5_kind_phys, 4.5_kind_phys /) dry_air_density(:,2) = (/ 5.5_kind_phys, 6.5_kind_phys /) + cloud_area_fraction(:,1) = (/ 0.1_kind_phys, 0.2_kind_phys /) + 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 /) filename_of_micm_configuration = 'musica_configurations/terminator/micm/config.json' filename_of_tuvx_configuration = 'musica_configurations/terminator/tuvx/config.json' @@ -406,7 +434,7 @@ subroutine test_terminator() stop 3 endif ASSERT(allocated(constituent_props)) - ASSERT(size(constituent_props) == NUM_SPECIES) + ASSERT(size(constituent_props) == NUM_SPECIES+NUM_TUVX_CONSTITUENTS) do i = 1, size(constituent_props) ASSERT(constituent_props(i)%is_instantiated(errcode, errmsg)) ASSERT(errcode == 0) @@ -417,7 +445,9 @@ subroutine test_terminator() call constituent_props(i)%is_advected(is_advected, errcode, errmsg) ASSERT(errcode == 0) tmp_bool = (trim(species_name) == "Cl" .and. molar_mass == 0.035453_kind_phys .and. is_advected) .or. & - (trim(species_name) == "Cl2" .and. molar_mass == 0.070906_kind_phys .and. is_advected) + (trim(species_name) == "Cl2" .and. molar_mass == 0.070906_kind_phys .and. is_advected) .or. & + (trim(species_name) == "cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water" & + .and. molar_mass == 0.018_kind_phys .and. is_advected) ASSERT(tmp_bool) call constituent_props(i)%units(units, errcode, errmsg) if (errcode /= 0) then @@ -438,7 +468,7 @@ subroutine test_terminator() end do call musica_ccpp_init(NUM_LAYERS, NUM_LAYERS+1, photolysis_wavelength_grid_interfaces, & - errmsg, errcode) + constituent_props_ptr, errmsg, errcode) if (errcode /= 0) then write(*,*) trim(errmsg) stop 3 @@ -462,6 +492,12 @@ subroutine test_terminator() end do end do end do + ! set initial cloud liquid water mixing ratio to ~1e-3 kg kg-1 + do j = 1, NUM_COLUMNS + do k = 1, NUM_LAYERS + constituents(j,k,NUM_SPECIES+1) = 1.0e-3_kind_phys * (1.0 + 0.1 * (j-1) + 0.01 * (k-1)) + end do + end do initial_constituents(:,:,:) = constituents(:,:,:) write(*,*) "[MUSICA INFO] Initial Time Step" @@ -477,7 +513,7 @@ subroutine test_terminator() constituents, geopotential_height_wrt_surface_at_midpoint, & geopotential_height_wrt_surface_at_interface, surface_temperature, & surface_geopotential, surface_albedo, standard_gravitational_acceleration, & - errmsg, errcode) + cloud_area_fraction, air_pressure_thickness, errmsg, errcode) if (errcode /= 0) then write(*,*) trim(errmsg) stop 3 @@ -504,6 +540,8 @@ subroutine test_terminator() total_Cl = constituents(i,j,Cl_index) + constituents(i,j,Cl2_index) total_Cl_init = initial_constituents(i,j,Cl_index) + initial_constituents(i,j,Cl2_index) ASSERT_NEAR(total_Cl, total_Cl_init, 1.0e-13) + ! cloud liquid water should be unchanged + ASSERT_NEAR(constituents(i,j,NUM_SPECIES+1), initial_constituents(i,j,NUM_SPECIES+1), 1.0e-13) end do end do diff --git a/test/musica/tuvx/CMakeLists.txt b/test/musica/tuvx/CMakeLists.txt index 2c8623fd..bb780942 100644 --- a/test/musica/tuvx/CMakeLists.txt +++ b/test/musica/tuvx/CMakeLists.txt @@ -110,4 +110,34 @@ add_test( WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} ) -add_memory_check_test(test_tuvx_surface_albedo $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) \ No newline at end of file +add_memory_check_test(test_tuvx_surface_albedo $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) + +# Cloud optics +add_executable(test_tuvx_cloud_optics test_tuvx_cloud_optics.F90) + +target_sources(test_tuvx_cloud_optics + PUBLIC + ${MUSICA_SRC_PATH}/tuvx/musica_ccpp_tuvx_height_grid.F90 + ${MUSICA_SRC_PATH}/tuvx/musica_ccpp_tuvx_wavelength_grid.F90 + ${MUSICA_SRC_PATH}/tuvx/musica_ccpp_tuvx_cloud_optics.F90 + ${MUSICA_SRC_PATH}/musica_ccpp_util.F90 + ${CCPP_TEST_SRC_PATH}/ccpp_kinds.F90 +) + +target_link_libraries(test_tuvx_cloud_optics + PRIVATE + musica::musica-fortran +) + +set_target_properties(test_tuvx_cloud_optics + PROPERTIES + LINKER_LANGUAGE Fortran +) + +add_test( + NAME test_tuvx_cloud_optics + COMMAND $ + WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} +) + +add_memory_check_test(test_tuvx_cloud_optics $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) diff --git a/test/musica/tuvx/test_tuvx_cloud_optics.F90 b/test/musica/tuvx/test_tuvx_cloud_optics.F90 new file mode 100644 index 00000000..b01e75a8 --- /dev/null +++ b/test/musica/tuvx/test_tuvx_cloud_optics.F90 @@ -0,0 +1,115 @@ +! Copyright (C) 2024 National Science Foundation-National Center for Atmospheric Research +! SPDX-License-Identifier: Apache-2.0 +program test_tuvx_cloud_optics + + use musica_ccpp_tuvx_cloud_optics + +#define ASSERT(x) if (.not.(x)) then; write(*,*) "Assertion failed[", __FILE__, ":", __LINE__, "]: x"; stop 1; endif +#define ASSERT_NEAR( a, b, abs_error ) if( (abs(a - b) >= abs_error) .and. (abs(a - b) /= 0.0) ) then; write(*,*) "Assertion failed[", __FILE__, ":", __LINE__, "]: a, b"; stop 1; endif + + real, parameter :: ABS_ERROR = 1e-5 + + call test_create_cloud_optics_radiator() + +contains + + subroutine test_create_cloud_optics_radiator() + + use musica_util, only: error_t + use musica_ccpp_tuvx_height_grid, only: create_height_grid + use musica_ccpp_tuvx_wavelength_grid, only: create_wavelength_grid + use musica_tuvx_grid, only: grid_t + use musica_tuvx_radiator, only: radiator_t + use ccpp_kinds, only: kind_phys + + integer, parameter :: NUM_HOST_HEIGHT_MIDPOINTS = 2 + integer, parameter :: NUM_HOST_HEIGHT_INTERFACES = 3 + integer, parameter :: NUM_WAVELENGTH_MIDPOINTS = 3 + integer, parameter :: NUM_WAVELENGTH_INTERFACES = 4 + real(kind_phys) :: host_wavelength_interfaces(NUM_WAVELENGTH_INTERFACES) = [180.0e-9_kind_phys, 200.0e-9_kind_phys, 240.0e-9_kind_phys, 300.0e-9_kind_phys] + real(kind_phys) :: delta_pressure(NUM_HOST_HEIGHT_MIDPOINTS) = [100.0_kind_phys, 200.0_kind_phys] + real(kind_phys) :: cloud_fraction(NUM_HOST_HEIGHT_MIDPOINTS) = [0.1_kind_phys, 0.0_kind_phys] + real(kind_phys) :: liquid_water_content(NUM_HOST_HEIGHT_MIDPOINTS) = [0.0003_kind_phys, 0.0004_kind_phys] + real(kind_phys) :: reciprocal_of_gravitational_acceleration = 0.1_kind_phys + real(kind_phys) :: cloud_optical_depth(NUM_HOST_HEIGHT_MIDPOINTS+1, NUM_WAVELENGTH_MIDPOINTS) + real(kind_phys) :: expected_cloud_optical_depth(NUM_HOST_HEIGHT_MIDPOINTS+1, NUM_WAVELENGTH_MIDPOINTS) = & + reshape([ 0.0_kind_phys, 0.14704591_kind_phys, 0.0_kind_phys, 0.0_kind_phys, 0.14704591_kind_phys, 0.0_kind_phys, 0.0_kind_phys, 0.14704591_kind_phys, 0.0_kind_phys, 0.0_kind_phys, 0.14704591_kind_phys, 0.0_kind_phys ], & + [ NUM_HOST_HEIGHT_MIDPOINTS+1, NUM_WAVELENGTH_MIDPOINTS ]) + real(kind_phys) :: single_scattering_albedo(NUM_HOST_HEIGHT_MIDPOINTS+1, NUM_WAVELENGTH_MIDPOINTS) + real(kind_phys) :: asymmetry_parameter(NUM_HOST_HEIGHT_MIDPOINTS+1, NUM_WAVELENGTH_MIDPOINTS,1) + type(grid_t), pointer :: height_grid => null() + type(grid_t), pointer :: wavelength_grid => null() + type(radiator_t), pointer :: clouds => null() + type(error_t) :: error + character(len=512) :: errmsg + integer :: errcode + integer :: i + + height_grid => create_height_grid(NUM_HOST_HEIGHT_MIDPOINTS, NUM_HOST_HEIGHT_INTERFACES, & + errmsg, errcode) + ASSERT(errcode == 0) + ASSERT(associated(height_grid)) + + wavelength_grid => create_wavelength_grid(host_wavelength_interfaces, errmsg, errcode) + ASSERT(errcode == 0) + ASSERT(associated(wavelength_grid)) + + clouds => create_cloud_optics_radiator(height_grid, wavelength_grid, errmsg, errcode) + ASSERT(errcode == 0) + ASSERT(associated(clouds)) + + call set_cloud_optics_values(clouds, cloud_fraction, [ 0.0_kind_phys ], & + liquid_water_content, & + reciprocal_of_gravitational_acceleration, & + errmsg, errcode) + ASSERT(errcode == 1) + + call set_cloud_optics_values(clouds, cloud_fraction, delta_pressure, & + [ 1.0_kind_phys ], & + reciprocal_of_gravitational_acceleration, & + errmsg, errcode) + ASSERT(errcode == 1) + + call set_cloud_optics_values(clouds, [ 1.0_kind_phys ], delta_pressure, & + liquid_water_content, & + reciprocal_of_gravitational_acceleration, & + errmsg, errcode) + ASSERT(errcode == 1) + + call set_cloud_optics_values(clouds, cloud_fraction, delta_pressure, & + liquid_water_content, & + reciprocal_of_gravitational_acceleration, & + errmsg, errcode) + ASSERT(errcode == 0) + + call clouds%get_optical_depths(cloud_optical_depth, error) + ASSERT(error%is_success()) + do i = 1, size(cloud_optical_depth, dim=1) + do j = 1, size(cloud_optical_depth, dim=2) + ASSERT_NEAR(cloud_optical_depth(i,j), expected_cloud_optical_depth(i,j), ABS_ERROR) + end do + end do + + call clouds%get_single_scattering_albedos(single_scattering_albedo, error) + ASSERT(error%is_success()) + do i = 1, size(single_scattering_albedo, dim=1) + do j = 1, size(single_scattering_albedo, dim=2) + ASSERT_NEAR(single_scattering_albedo(i,j), 0.0_kind_phys, ABS_ERROR) + end do + end do + + call clouds%get_asymmetry_factors(asymmetry_parameter, error) + ASSERT(error%is_success()) + do i = 1, size(asymmetry_parameter, dim=1) + do j = 1, size(asymmetry_parameter, dim=2) + ASSERT_NEAR(asymmetry_parameter(i,j,1), 0.0_kind_phys, ABS_ERROR) + end do + end do + + deallocate( height_grid ) + deallocate( wavelength_grid ) + deallocate( clouds ) + + end subroutine test_create_cloud_optics_radiator + +end program test_tuvx_cloud_optics diff --git a/to_be_ccppized/ccpp_const_utils.F90 b/to_be_ccppized/ccpp_const_utils.F90 index 902d6059..e20aed31 100644 --- a/to_be_ccppized/ccpp_const_utils.F90 +++ b/to_be_ccppized/ccpp_const_utils.F90 @@ -13,13 +13,13 @@ subroutine ccpp_const_get_idx(constituent_props, name, cindex, errmsg, errflg) use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t ! Input arguments - type(ccpp_constituent_prop_ptr_t), pointer, intent(in) :: constituent_props(:) - character(len=*), intent(in) :: name ! constituent name + type(ccpp_constituent_prop_ptr_t), intent(in) :: constituent_props(:) + character(len=*), intent(in) :: name ! constituent name ! Output arguments - integer, intent(out) :: cindex ! global constituent index - character(len=512), intent(out) :: errmsg ! error message - integer, intent(out) :: errflg ! error flag + integer, intent(out) :: cindex ! global constituent index + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag ! Local variables integer :: t_cindex From 2a9d9ab205cca8366ce2a2b6d838323063af2505 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Fri, 22 Nov 2024 16:10:26 -0800 Subject: [PATCH 02/14] add valgrind suppressions --- test/valgrind.supp | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/test/valgrind.supp b/test/valgrind.supp index ed4ee497..ee1ba85e 100644 --- a/test/valgrind.supp +++ b/test/valgrind.supp @@ -35,4 +35,47 @@ fun:__tuvx_core_MOD_constructor fun:InternalCreateTuvx ... +} +{ + Suppress_MUSICA_TUV-x_CreateRadiator + Memcheck:Leak + match-leak-kinds: definite + fun:malloc + fun:__musica_string_MOD_string_assign_char + fun:__tuvx_radiator_from_host_MOD_constructor_char + fun:__tuvx_radiator_from_host_MOD_constructor_string + fun:InternalCreateRadiator + ... +} +{ + Suppress_MUSICA_TUV-x_AddRadiator + Memcheck:Leak + match-leak-kinds: definite + fun:malloc + fun:__tuvx_radiator_from_host_MOD___copy_tuvx_radiator_from_host_Radiator_from_host_t + fun:__tuvx_radiator_warehouse_MOD_add_radiator + fun:InternalAddRadiator + ... +} +{ + Suppress_MUSICA_TUV-x_GetRadiator + Memcheck:Leak + match-leak-kinds: definite + fun:malloc + fun:__tuvx_radiator_from_host_MOD___copy_tuvx_radiator_from_host_Radiator_from_host_t + fun:InternalGetRadiator + ... +} +{ + Suppress_MUSICA_TUV-x_CreateTuvx-RadiatorFromHost + Memcheck:Leak + match-leak-kinds: definite + fun:malloc + fun:__tuvx_radiator_from_host_MOD___copy_tuvx_radiator_from_host_Radiator_from_host_t + fun:__tuvx_radiator_warehouse_MOD_add_radiator + fun:__tuvx_radiator_warehouse_MOD_add_radiators + fun:__tuvx_radiative_transfer_MOD_constructor + fun:__tuvx_core_MOD_constructor + fun:InternalCreateTuvx + ... } \ No newline at end of file From 645ec882b28cff822c5cde73b2ecbb93e99ca9d0 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Fri, 22 Nov 2024 16:49:41 -0800 Subject: [PATCH 03/14] fix typo --- schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 index 09ae02e8..95b695ed 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 @@ -11,7 +11,7 @@ module musica_ccpp_tuvx_cloud_optics ! and include the cloud optical depth, single scattering albedo, ! and asymmetry parameter. ! - ! See module_ccpp_tuvx_height_grid for the definition of the height grid + ! See musica_ccpp_tuvx_height_grid for the definition of the height grid ! and its mapping to the CAM-SIMA vertical grid. !> Label for cloud optical properties in TUV-x From 4216e13c86832d32f061bb3be1562688d2e95b72 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Fri, 22 Nov 2024 16:53:59 -0800 Subject: [PATCH 04/14] fix formatting --- schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 index 95b695ed..72f63f09 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 @@ -35,11 +35,11 @@ function create_cloud_optics_radiator( height_grid, wavelength_grid, & use musica_tuvx_radiator, only: radiator_t use musica_util, only: error_t - type(grid_t), intent(inout) :: height_grid - type(grid_t), intent(inout) :: wavelength_grid - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errcode - type(radiator_t), pointer :: radiator + type(grid_t), intent(inout) :: height_grid + type(grid_t), intent(inout) :: wavelength_grid + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errcode + type(radiator_t), pointer :: radiator ! local variables type(error_t) :: error From 4c11684344e434ac5f277d1dc5e8d25a443bee55 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Mon, 25 Nov 2024 09:09:29 -0800 Subject: [PATCH 05/14] address reviewer comments --- schemes/musica/musica_ccpp.F90 | 4 ++-- .../tuvx/musica_ccpp_tuvx_cloud_optics.F90 | 22 +++++++++---------- test/musica/CMakeLists.txt | 2 +- 3 files changed, 13 insertions(+), 15 deletions(-) diff --git a/schemes/musica/musica_ccpp.F90 b/schemes/musica/musica_ccpp.F90 index 321356f2..f2aaedb2 100644 --- a/schemes/musica/musica_ccpp.F90 +++ b/schemes/musica/musica_ccpp.F90 @@ -72,8 +72,8 @@ 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, & - standard_gravitational_acceleration, cloud_area_fraction, & + photolysis_wavelength_grid_interfaces, extraterrestrial_flux, & + standard_gravitational_acceleration, cloud_area_fraction, & air_pressure_thickness, errmsg, errcode) use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t use ccpp_kinds, only: kind_phys diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 index 72f63f09..dad013f5 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 @@ -67,8 +67,8 @@ subroutine set_cloud_optics_values( radiator, cloud_fraction, delta_pressure, & use musica_util, only: error_t type(radiator_t), intent(inout) :: radiator - real(kind_phys), intent(in) :: delta_pressure(:) ! pressure delta about vertical level midpoints (Pa) real(kind_phys), intent(in) :: cloud_fraction(:) ! (unitless) + real(kind_phys), intent(in) :: delta_pressure(:) ! pressure delta about vertical level midpoints (Pa) real(kind_phys), intent(in) :: cloud_liquid_water_content(:) ! (kg/kg) real(kind_phys), intent(in) :: reciprocal_of_gravitational_acceleration ! (s^2/m) character(len=*), intent(out) :: errmsg @@ -76,23 +76,21 @@ subroutine set_cloud_optics_values( radiator, cloud_fraction, delta_pressure, & ! local variables type(error_t) :: error - real(kind_phys) :: od(num_vertical_levels) + real(kind_phys) :: optical_depth(num_vertical_levels) ! working array for cloud optical depth real(kind_phys) :: cloud_optical_depth(num_vertical_levels, num_wavelength_bins) - integer :: n_host_midpoints integer :: i_level - n_host_midpoints = size(cloud_fraction) - if ( size(delta_pressure) /= n_host_midpoints ) then + if ( size(delta_pressure) /= size(cloud_fraction) ) then errmsg = "[MUSICA Error] Invalid size of cloud pressure delta for TUV-x." errcode = 1 return end if - if ( size(cloud_liquid_water_content) /= n_host_midpoints ) then + if ( size(cloud_liquid_water_content) /= size(cloud_fraction) ) then errmsg = "[MUSICA Error] Invalid size of cloud liquid water content for TUV-x." errcode = 1 return end if - if ( n_host_midpoints + 1 /= num_vertical_levels ) then + if ( size(cloud_fraction) + 1 /= num_vertical_levels ) then errmsg = "[MUSICA Error] Invalid size of cloud fraction for TUV-x." errcode = 1 return @@ -107,17 +105,17 @@ subroutine set_cloud_optics_values( radiator, cloud_fraction, delta_pressure, & ! The cloud optical depth is then estimated as: ! od = lwp * 155 * cf^1.5 ! A constant cloud optical depth is used for all wavelengths. - do i_level = 1, n_host_midpoints + do i_level = 1, size(cloud_fraction) if ( cloud_fraction(i_level) > 0.0_kind_phys ) then - od(i_level) = ( reciprocal_of_gravitational_acceleration & + optical_depth(i_level) = ( reciprocal_of_gravitational_acceleration & * cloud_liquid_water_content(i_level) * delta_pressure(i_level) & / cloud_fraction(i_level) ) * 155.0_kind_phys * cloud_fraction(i_level)**1.5_kind_phys else - od(i_level) = 0.0_kind_phys + optical_depth(i_level) = 0.0_kind_phys end if end do - do i_level = 1, n_host_midpoints - cloud_optical_depth(i_level, :) = od(n_host_midpoints-i_level+1) + do i_level = 1, size(cloud_fraction) + cloud_optical_depth(i_level, :) = optical_depth(size(cloud_fraction)-i_level+1) end do cloud_optical_depth(num_vertical_levels, :) = 0.0_kind_phys call radiator%set_optical_depths( cloud_optical_depth, error ) diff --git a/test/musica/CMakeLists.txt b/test/musica/CMakeLists.txt index 0d36d9f0..951b0b1e 100644 --- a/test/musica/CMakeLists.txt +++ b/test/musica/CMakeLists.txt @@ -29,11 +29,11 @@ target_sources(test_musica_api PUBLIC ${MUSICA_CCPP_SOURCES} ${TO_BE_CCPPIZED_SRC_PATH}/ccpp_tuvx_utils.F90 + ${TO_BE_CCPPIZED_SRC_PATH}/ccpp_const_utils.F90 ${CCPP_SRC_PATH}/ccpp_constituent_prop_mod.F90 ${CCPP_SRC_PATH}/ccpp_hash_table.F90 ${CCPP_SRC_PATH}/ccpp_hashable.F90 ${CCPP_SRC_PATH}/ccpp_types.F90 - ${TO_BE_CCPPIZED_SRC_PATH}/ccpp_const_utils.F90 ${CCPP_TEST_SRC_PATH}/ccpp_kinds.F90 ) From 2f8afe6a26c9eefac81b0a27b72007949dcb42f4 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Mon, 25 Nov 2024 09:55:54 -0800 Subject: [PATCH 06/14] address reviewer comments --- .../tuvx/musica_ccpp_tuvx_cloud_optics.F90 | 21 ++++++++++--------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 index dad013f5..8632176f 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_cloud_optics.F90 @@ -78,20 +78,21 @@ subroutine set_cloud_optics_values( radiator, cloud_fraction, delta_pressure, & type(error_t) :: error real(kind_phys) :: optical_depth(num_vertical_levels) ! working array for cloud optical depth real(kind_phys) :: cloud_optical_depth(num_vertical_levels, num_wavelength_bins) - integer :: i_level + integer :: i_level, size_cloud_fraction - if ( size(delta_pressure) /= size(cloud_fraction) ) then - errmsg = "[MUSICA Error] Invalid size of cloud pressure delta for TUV-x." + size_cloud_fraction = size(cloud_fraction) + if ( size_cloud_fraction + 1 /= num_vertical_levels ) then + errmsg = "[MUSICA Error] Invalid size of cloud fraction for TUV-x." errcode = 1 return end if - if ( size(cloud_liquid_water_content) /= size(cloud_fraction) ) then - errmsg = "[MUSICA Error] Invalid size of cloud liquid water content for TUV-x." + if ( size(delta_pressure) /= size_cloud_fraction ) then + errmsg = "[MUSICA Error] Invalid size of cloud pressure delta for TUV-x." errcode = 1 return end if - if ( size(cloud_fraction) + 1 /= num_vertical_levels ) then - errmsg = "[MUSICA Error] Invalid size of cloud fraction for TUV-x." + if ( size(cloud_liquid_water_content) /= size_cloud_fraction ) then + errmsg = "[MUSICA Error] Invalid size of cloud liquid water content for TUV-x." errcode = 1 return end if @@ -105,7 +106,7 @@ subroutine set_cloud_optics_values( radiator, cloud_fraction, delta_pressure, & ! The cloud optical depth is then estimated as: ! od = lwp * 155 * cf^1.5 ! A constant cloud optical depth is used for all wavelengths. - do i_level = 1, size(cloud_fraction) + do i_level = 1, size_cloud_fraction if ( cloud_fraction(i_level) > 0.0_kind_phys ) then optical_depth(i_level) = ( reciprocal_of_gravitational_acceleration & * cloud_liquid_water_content(i_level) * delta_pressure(i_level) & @@ -114,8 +115,8 @@ subroutine set_cloud_optics_values( radiator, cloud_fraction, delta_pressure, & optical_depth(i_level) = 0.0_kind_phys end if end do - do i_level = 1, size(cloud_fraction) - cloud_optical_depth(i_level, :) = optical_depth(size(cloud_fraction)-i_level+1) + do i_level = 1, size_cloud_fraction + cloud_optical_depth(i_level, :) = optical_depth(size_cloud_fraction-i_level+1) end do cloud_optical_depth(num_vertical_levels, :) = 0.0_kind_phys call radiator%set_optical_depths( cloud_optical_depth, error ) From 4454a161859f65ec2467b580ca16ad824c5d377c Mon Sep 17 00:00:00 2001 From: Jiwon Gim Date: Mon, 2 Dec 2024 11:17:48 -0700 Subject: [PATCH 07/14] add species --- .gitignore | 3 +- schemes/musica/musica_ccpp.F90 | 1 + schemes/musica/tuvx/musica_ccpp_tuvx.F90 | 34 +- .../musica_ccpp_tuvx_gas_species_profile.F90 | 455 ++++++++++++++++++ test/musica/test_musica_api.F90 | 26 +- 5 files changed, 502 insertions(+), 17 deletions(-) create mode 100644 schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 diff --git a/.gitignore b/.gitignore index 83c5d6be..8639122f 100644 --- a/.gitignore +++ b/.gitignore @@ -4,4 +4,5 @@ **/CMakeCache.txt **/CMakeFiles/ .vscode -xcode/ \ No newline at end of file +xcode/ +*.mod \ No newline at end of file diff --git a/schemes/musica/musica_ccpp.F90 b/schemes/musica/musica_ccpp.F90 index f2aaedb2..9d922226 100644 --- a/schemes/musica/musica_ccpp.F90 +++ b/schemes/musica/musica_ccpp.F90 @@ -122,6 +122,7 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co air_pressure_thickness, rate_parameters, & errmsg, errcode) + ! TODO(jiwon) this might not be correct because it doesn't know the index ! Get the molar mass that is set in the call to instantiate() do i_elem = 1, size(molar_mass_arr) call constituent_props(i_elem)%molar_mass(molar_mass_arr(i_elem), errcode, errmsg) diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 index df20b664..2805d4fa 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 @@ -48,7 +48,7 @@ subroutine reset_tuvx_map_state( grids, profiles, radiators ) end subroutine reset_tuvx_map_state - !> This is a helper subroutine created to deallocate objects associated with TUV-x + !> Deallocates objects associated with TUV-x subroutine cleanup_tuvx_resources() if (associated( height_grid )) then @@ -142,6 +142,8 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & extraterrestrial_flux_unit use musica_ccpp_tuvx_cloud_optics, & only: create_cloud_optics_radiator, cloud_optics_label + use musica_ccpp_tuvx_gas_species_profilev, & + only: gas_species_t, configure_gas_species, create_gas_species_profile_group integer, intent(in) :: vertical_layer_dimension ! (count) integer, intent(in) :: vertical_interface_dimension ! (count) @@ -152,12 +154,14 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & integer, intent(out) :: errcode ! local variables - type(grid_map_t), pointer :: grids - type(profile_map_t), pointer :: profiles - type(radiator_map_t), pointer :: radiators - type(configuration_t) :: config - type(mappings_t), pointer :: photolysis_rate_constants_ordering - type(error_t) :: error + type(error_t) :: error + type(grid_map_t), pointer :: grids + type(profile_map_t), pointer :: profiles + type(radiator_map_t), pointer :: radiators + type(configuration_t) :: config + type(mappings_t), pointer :: photolysis_rate_constants_ordering + type(gas_species_t), allocatable :: gas_species_group(:) + type(profile_t), pointer :: profile_gas_species_group(:) ! Get needed indices in constituents array call ccpp_const_get_idx(constituent_props, CLOUD_LIQUID_WATER_CONTENT_LABEL, & @@ -247,6 +251,22 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & return end if + call configure_gas_species( constituent_props, gas_species_group, & + errmsg, errcode ) + if (errcode /= 0) then + call reset_tuvx_map_state( grids, profiles, null() ) + call cleanup_tuvx_resources() + return + end if + + call create_gas_species_profile_group( height_grid, & + gas_species_group, profile_gas_species_group, errmsg, errcode) + if (errcode /= 0) then + call reset_tuvx_map_state( grids, profiles, null() ) + call cleanup_tuvx_resources() + return + endif + radiators => radiator_map_t( error ) if (has_error_occurred( error, errmsg, errcode )) then call reset_tuvx_map_state( grids, profiles, null() ) diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 new file mode 100644 index 00000000..a5e99220 --- /dev/null +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 @@ -0,0 +1,455 @@ +module musica_ccpp_tuvx_gas_species_profile + use ccpp_kinds, only: kind_phys + + implicit none + + private + public :: create_gas_species, configure_gas_species, create_gas_species_profile_group + + !> Conversion factor from km to cm + real(kind_phys), parameter, public :: km_to_cm = 1.0e5 + !> Conversion factor from m3 to cm3 + real(kind_phys), parameter, public :: m_3_to_cm_3 = 1.0e6 + !> Default value of number of wavelength grid bins + integer, parameter :: DEFAULT_NUM_GRID_SECTIONS = 0 + !> Number of wavelength grid bins + integer, protected :: num_midpoint_heights_ = DEFAULT_NUM_GRID_SECTIONS + integer, protected :: num_interfaces_ = DEFAULT_NUM_GRID_SECTIONS + + type, public :: gas_species_t + character(len=50) :: name + character(len=20) :: unit + real(kind_phys) :: molar_mass + real(kind_phys) :: scale_height + integer :: index_constituent_props + end type gas_species_t + +contains + + function create_gas_species(name, unit, molar_mass, scale_height, & + index_constituent_props) result( gas_species ) + use ccpp_kinds, only: kind_phys + + character(len=*), intent(in) :: name + character(len=*), intent(in) :: unit + real(kind_phys), intent(in) :: molar_mass + real(kind_phys), intent(in) :: scale_height + integer, intent(in) :: index_constituent_props + type(gas_species_t) :: gas_species + + gas_species%name = name + gas_species%unit = unit + gas_species%molar_mass = molar_mass + gas_species%scale_height = scale_height + gas_species%index_constituent_props = index_constituent_props + + end function create_gas_species + + !> TODO(jiwon) - write function description, error message + subroutine configure_gas_species(constituent_props, gas_species_group, errmsg, errcode) + use ccpp_kinds, only: kind_phys + use ccpp_const_utils, only: ccpp_const_get_idx + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t + + type(ccpp_constituent_prop_ptr_t), intent(in) :: constituent_props(:) + type(gas_species_t), allocatable, intent(out) :: gas_species_group(:) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode + + ! local variable + integer :: number_of_gas_species = 3 + ! The value of molar mass of dry air is from 'CAM-SIMA/src/utils/std_atm_profile.F90' + real(kind_phys), parameter :: SCALE_HEIGHT_DRY_AIR = 8.01_kind_phys ! km + real(kind_phys), parameter :: SCALE_HEIGHT_O2 = 7.0_kind_phys ! km + real(kind_phys), parameter :: SCALE_HEIGHT_O3 = 7.0_kind_phys ! km + real(kind_phys), parameter :: MOLAR_MASS_DRY_AIR = 0.0289644_kind_phys ! kg mol-1 + real(kind_phys) :: molar_mass_O2 ! kg mol-1 + real(kind_phys) :: molar_mass_O3 ! kg mol-1 + + integer :: index_air, index_O2, index_O3 + allocate(gas_species_group(number_of_gas_species)) ! TODO(jiwon) - when deallocate this? + + ! Air + ! Find the index of the species and create a 'gas_species_t' for it + ! call ccpp_const_get_idx(constituent_props, "air", index_air, errmsg, errcode) + ! write(*,*) "air index: ", index_air + ! if (errcode /= 0) return + + ! TODO(jiwon) air may not be registered + ! call constituent_props(index_air)%molar_mass(molar_mass_air, errcode, errmsg) + ! if (errcode /= 0) return + + gas_species_group(1) = create_gas_species("air", "molecule cm-3", & + MOLAR_MASS_DRY_AIR, SCALE_HEIGHT_DRY_AIR, index_air) + ! O2 + call ccpp_const_get_idx(constituent_props, "O2", index_O2, errmsg, errcode) + write(*,*) "air index: ", index_O2 + if (errcode /= 0) return + + call constituent_props(index_O2)%molar_mass(molar_mass_O2, errcode, errmsg) + if (errcode /= 0) return + + gas_species_group(2) = create_gas_species("O2", "molecule cm-3", & + molar_mass_O2, SCALE_HEIGHT_O2, index_O2) + + ! O3 + call ccpp_const_get_idx(constituent_props, "O3", index_O3, errmsg, errcode) + if (errcode /= 0) return + + call constituent_props(index_O3)%molar_mass(molar_mass_O3, errcode, errmsg) + if (errcode /= 0) return + + gas_species_group(3) = create_gas_species("O3", "molecule cm-3", & + molar_mass_O3, SCALE_HEIGHT_O3, index_O3) + + end subroutine configure_gas_species + + !> Creates + subroutine create_gas_species_profile_group(height_grid, gas_species_group, & + profile_group, errmsg, errcode) + use musica_ccpp_util, only: has_error_occurred + use musica_tuvx_grid, only: grid_t + use musica_tuvx_profile, only: profile_t + use musica_util, only: error_t + + type(grid_t), intent(in) :: height_grid + type(gas_species_t), intent(in) :: gas_species_group(:) + type(profile_t), pointer, intent(out) :: profile_group(:) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode + + ! local variables + type(error_t) :: error + integer :: number_of_gas_species + integer :: i_species + + number_of_gas_species = size(gas_species_group) + allocate( profile_group( number_of_gas_species ) ) + + do i_species = 1, number_of_gas_species + allocate(profile_group(i_species)) + write(*,*) "gas_species_group(i_species)%name: ", gas_species_group(i_species)%name + profile_group(i_species) = profile_t( gas_species_group(i_species)%name, & + gas_species_group(i_species)%unit, height_grid, error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + ! TODO(jiwon) deallocate the previously allocated one? + end do + + end subroutine create_gas_species_profile_group + + ! function get_height_delta(height_grid, num_interfaces, errmsg, errcode) result( height_delta ) + ! 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(in) :: height_grid + ! integer, intent(in) :: num_interfaces + ! character(len=*), intent(out) :: errmsg + ! integer, intent(out) :: errcode + + ! ! local variables + ! type(error_t) :: error + ! real(kind_phys) :: interfaces(num_interfaces) + ! real(kind_phys) :: height_delta(num_interfaces-1) + ! integer :: i_elem + + ! call height_grid%get_edges(interfaces, error) + ! if ( has_error_occurred( error, errmsg, errcode ) ) return + + ! do i_elem = 1, size(height_delta) + ! height_delta[i] = interfaces[i+1] - interfaces[i] + ! end do + + ! end function get_height_delta + + ! !> Creates a TUV-x surface albedo profile from the host-model wavelength grid + ! function create_air_profile( height_grid, errmsg, errcode ) & + ! result( profile ) + ! use musica_ccpp_util, only: has_error_occurred + ! use musica_tuvx_grid, only: grid_t + ! use musica_tuvx_profile, only: profile_t + ! use musica_util, only: error_t + + ! type(grid_t), intent(in) :: height_grid + ! character(len=*), intent(out) :: errmsg + ! integer, intent(out) :: errcode + ! type(profile_t), pointer :: profile + + ! ! local variables + ! type(error_t) :: error + + ! profile => profile_t( air_label, air_unit, height_grid, error ) + ! if ( has_error_occurred( error, errmsg, errcode ) ) return + + ! end function create_air_profile + + ! !> Creates a TUV-x surface albedo profile from the host-model wavelength grid + ! function create_O2_profile( height_grid, errmsg, errcode ) & + ! result( profile ) + ! use musica_ccpp_util, only: has_error_occurred + ! use musica_tuvx_grid, only: grid_t + ! use musica_tuvx_profile, only: profile_t + ! use musica_util, only: error_t + + ! type(grid_t), intent(in) :: height_grid + ! character(len=*), intent(out) :: errmsg + ! integer, intent(out) :: errcode + ! type(profile_t), pointer :: profile + + ! ! local variables + ! type(error_t) :: error + + ! profile => profile_t( O2_label, O2_unit, height_grid, error ) + ! if ( has_error_occurred( error, errmsg, errcode ) ) return + + ! end function create_O2_profile + + ! !> Creates a TUV-x surface albedo profile from the host-model wavelength grid + ! function create_O3_profile( height_grid, errmsg, errcode ) & + ! result( profile ) + ! use musica_ccpp_util, only: has_error_occurred + ! use musica_tuvx_grid, only: grid_t + ! use musica_tuvx_profile, only: profile_t + ! use musica_util, only: error_t + + ! type(grid_t), intent(in) :: height_grid + ! character(len=*), intent(out) :: errmsg + ! integer, intent(out) :: errcode + ! type(profile_t), pointer :: profile + + ! ! local variables + ! type(error_t) :: error + + ! profile => profile_t( O3_label, O3_unit, height_grid, error ) + ! if ( has_error_occurred( error, errmsg, errcode ) ) return + + ! end function create_O3_profile + + ! subroutine create_gas_species_profiles( air_profile, O2_profile, O3_profile, & + ! height_grid, errmsg, errcode ) + ! use musica_ccpp_util, only: has_error_occurred + ! use musica_tuvx_grid, only: grid_t + ! use musica_tuvx_profile, only: profile_t + ! use musica_util, only: error_t + + ! type(profile_t), intent(inout) :: air_profile + ! type(profile_t), intent(inout) :: O2_profile + ! type(profile_t), intent(inout) :: O3_profile + ! type(grid_t), intent(in) :: height_grid + ! character(len=*), intent(out) :: errmsg + ! integer, intent(out) :: errcode + + ! ! local variables + ! type(error_t) :: error + + ! air_profile = create_air_profile( height_grid, errmsg, errcode ) + ! if (errcode /= 0) return + + ! O2_profile = create_O2_profile( height_grid, errmsg, errcode ) + ! if (errcode /= 0) return + + ! O3_profile = create_O3_profile( height_grid, errmsg, errcode ) + ! if (errcode /= 0) return + + ! num_midpoint_heights_ = height_grid%number_of_sections( error ) + ! if ( has_error_occurred( error, errmsg, errcode ) ) return + + ! num_interfaces_ = num_midpoint_heights_ + 1 + + ! ! Set height delta values + ! height_delta_ = get_height_delta( height_grid, num_interfaces_, errmsg, errcode ) + ! if ( has_error_occurred( error, errmsg, errcode ) ) return + + ! end subroutine create_gas_species_profiles + + ! !> Sets TUV-x surface albedo values + ! subroutine set_dry_air_density_values( profile, dry_air_density, errmsg, errcode ) + ! use ccpp_kinds, only: kind_phys + ! use musica_ccpp_util, only: has_error_occurred + ! use musica_tuvx_profile, only: profile_t + ! use musica_util, only: error_t + + ! type(profile_t), intent(inout) :: profile + ! real(kind_phys), intent(in) :: dry_air_density(:) ! TODO(jiwon) : molecule cm-3 for CAM, kg m-3 as input + ! integer, intent(in) :: index_air + ! character(len=*), intent(out) :: errmsg + ! integer, intent(out) :: errcode + + ! ! local variables + ! type(error_t) :: error + ! real(kind_phys) :: scale_height = 8.01_kind_phys ! km + ! integer :: num_midpoints = size(dry_air_density) + ! real(kind_phys) :: interface(num_midpoints + 2) + ! real(kind_phys) :: densities(num_midpoints + 1) + + ! interfaces(1) = dry_air_density(num_midpoints) + ! interfaces(2:num_midpoints+1) = dry_air_density(num_midpoints:1:-1) + ! interfaces(num_midpoints+2) = dry_air_density(1) ! use upper mid-point value for top edge + + ! densities(1:num_midpoints+1) = & + ! height_delta_(1:num_midpoints+1) * km_to_cm & ! TODO(jiwon) : molecule cm-3 for CAM, kg m-3 as input + ! * sqrt(interfaces(1:num_midpoints+1)) & + ! * sqrt(interfaces(2:num_midpoints+2)) + + ! call profile%set_edge_values( interfaces, error ) + ! if ( has_error_occurred( error, errmsg, errcode ) ) return + + ! call profile%set_layer_densities( densities, error) + ! if ( has_error_occurred( error, errmsg, errcode ) ) return + + ! call profile%calculate_exo_layer_density( scale_height, error ) + ! if ( has_error_occurred( error, errmsg, errcode ) ) return + + ! end subroutine set_air_density_values + + ! !> Sets TUV-x surface albedo values + ! subroutine set_O2_values(profile, O2_constituent, height_delta, errmsg, errcode) + ! use ccpp_kinds, only: kind_phys + ! use musica_ccpp_util, only: has_error_occurred + ! use musica_tuvx_profile, only: profile_t + ! use musica_util, only: error_t + + ! type(profile_t), intent(inout) :: profile + ! real(kind_phys), intent(in) :: O2_constituent(:) ! mol m-3 + ! real(kind_phys), intent(in) :: height_delta(:) ! km + ! character(len=*), intent(out) :: errmsg + ! integer, intent(out) :: errcode + + ! ! local variables + ! type(error_t) :: error + ! real(kind_phys), parameter :: SCALE_HEIGHT = 7.0_kind_phys ! km + ! real(kind_phys) :: O2_constituent_mol_cm-3(size(O2_constituent)) ! mol cm-3 + ! integer :: num_midpoints = size(O2_constituent) + ! real(kind_phys) :: interface(size(O2_constituent) + 2) + ! real(kind_phys) :: densities(size(O2_constituent) + 1) + + ! O2_constituent_mol_cm-3(:) = O2_constituent(:) * m_3_to_cm_3 + + ! ! TODO(jiwon) if ( is_parameterized ) then + + ! interfaces(1) = O2_constituent_mol_cm-3(num_midpoints) + ! interfaces(2:num_midpoints+1) = O2_constituent_mol_cm-3(num_midpoints:1:-1) + ! interfaces(num_midpoints+2) = O2_constituent_mol_cm-3(1) + + ! densities(1:num_midpoints+1) = & + ! height_delta(1:num_midpoints+1) * km_to_cm & + ! * sqrt(interfaces(1:num_midpoints+1)) & + ! * sqrt(interfaces(2:num_midpoints+2)) + + ! call profile%set_edge_values( interfaces, error ) + ! if ( has_error_occurred( error, errmsg, errcode ) ) return + + ! call profile%set_layer_densities( densities, error) + ! if ( has_error_occurred( error, errmsg, errcode ) ) return + + ! call profile%calculate_exo_layer_density( SCALE_HEIGHT, error ) + ! if ( has_error_occurred( error, errmsg, errcode ) ) return + + ! end subroutine set_O2_values + + ! !> Sets TUV-x surface albedo values + ! subroutine set_O3_values(profile, O3_constituent, height_delta, above_column_density, & + ! errmsg, errcode) + ! use ccpp_kinds, only: kind_phys + ! use musica_ccpp_util, only: has_error_occurred + ! use musica_tuvx_profile, only: profile_t + ! use musica_util, only: error_t + + ! type(profile_t), intent(inout) :: profile + ! real(kind_phys), intent(in) :: O2_constituent(:) ! mol m-3 + ! real(kind_phys), intent(in) :: height_delta(:) ! km + ! real(kind_phys), intent(in) :: above_column_density(:,:) ! molecule cm-2 + ! character(len=*), intent(out) :: errmsg + ! integer, intent(out) :: errcode + + ! ! local variables + ! type(error_t) :: error + ! real(kind_phys), parameter :: SCALE_HEIGHT = 7.0_kind_phys ! km + ! real(kind_phys) :: O2_constituent_mol_cm-3(size(O3_constituent)) ! mol cm-3 + ! integer :: num_midpoints = size(O3_constituent) + ! real(kind_phys) :: interface(size(O3_constituent) + 2) + ! real(kind_phys) :: densities(size(O3_constituent) + 1) + + ! O3_constituent_mol_cm-3(:) = O3_constituent(:) * m_3_to_cm_3 !TODO(devided maybe) + + ! interfaces(1) = O3_constituent_mol_cm-3(num_midpoints) + ! interfaces(2:num_midpoints+1) = O3_constituent_mol_cm-3(num_midpoints:1:-1) + ! interfaces(num_midpoints+2) = O3_constituent_mol_cm-3(1) + + ! call profile%set_edge_values( interfaces, error ) + ! if ( has_error_occurred( error, errmsg, errcode ) ) return + + ! densities(1:num_midpoints+1) = & + ! height_delta(1:num_midpoints+1) * km_to_cm & + ! * sqrt(interfaces(1:num_midpoints+1)) & + ! * sqrt(interfaces(2:num_midpoints+2)) + + ! if ( num_absorbing_species >= 1 ) then ! TODO - nabscol = 2 in chem_mods + ! densities(1) = 0.5_kind_phys * above_column_density(num_midpoints, 1) + ! densities(2:num_midpoints) = 0.5_kind_phys * ( above_column_density(num_midpoints-1:1:-1,1) & + ! + above_column_density(num_midpoints:2:-1,1) ) + ! densities(num_midpoints+1) = above_column_density(0,1) & + ! + 0.5_kind_phys * above_column_density(1,1) + + ! call profile%set_layer_densities( densities, error) + ! if ( has_error_occurred( error, errmsg, errcode ) ) return + + ! call profile%set_exo_layer_density( above_column_density(0, 1), error) ! TODO: should be incorrect? exo_density = above_column_density(0,1) ) + ! if ( has_error_occurred( error, errmsg, errcode ) ) return + ! else + ! densities(1:num_midpoints+1) = height_delta(1:num_midpoints+1) * km_to_cm & + ! * ( interfaces(1:num_midpoints+1) + interfaces(2:num_midpoints+2) ) * 0.5_kind_phys + + ! call profile%set_layer_densities( densities, error) + ! if ( has_error_occurred( error, errmsg, errcode ) ) return + + ! call profile%calculate_exo_layer_density( SCALE_HEIGHT, error ) + ! if ( has_error_occurred( error, errmsg, errcode ) ) return + ! end if + + ! end subroutine set_O3_values + + ! subroutine set_gas_species_values( air_profile, O2_profile, O3_profile & + ! fixed_species_density, species_volume_mixing_ratio, above_column_density, & + ! index_air, index_O2, index_O3, errmsg, errcode) + ! use ccpp_kinds, only: kind_phys + ! use musica_tuvx_profile, only: profile_t + ! use musica_util, only: error_t + + ! type(profile_t), intent(inout) :: air_profile + ! type(profile_t), intent(inout) :: O2_profile + ! type(profile_t), intent(inout) :: O3_profile + ! real(kind_phys), intent(in) :: fixed_species_density(:,:) ! molecule cm-3 + ! real(kind_phys), intent(in) :: species_volume_mixing_ratio(:,:) ! mol mol-1 + ! real(kind_phys), intent(in) :: above_column_density(:,:) ! molecule cm-2 + ! integer, intent(in) :: index_air + ! integer, intent(in) :: index_O2 + ! integer, intent(in) :: index_O3 + ! character(len=*), intent(out) :: errmsg + ! integer, intent(out) :: errcode + + ! ! local variables + ! real(kind_phys) :: air_density = fixed_species_density(:, index_air) + + + ! ! do the conversion to mol / cm3 + ! call set_air_density_values( air_profile, & + ! fixed_species_density(:, index_air), errmsg, errcode ) + ! if (errcode /= 0) return + + ! call set_O2_values( O2_profile, fixed_species_density(:, index_air), & + ! fixed_species_density(:, index_O2), & + ! species_volume_mixing_ratio(:, index_O2), & + ! errmsg, errcode ) + ! if (errcode /= 0) return + + ! call set_O3_values( O3_profile, fixed_species_density(:, index_air), & + ! fixed_species_density(:, index_O3), & + ! species_volume_mixing_ratio(:, index_O3), & + ! above_column_density(:, 1), errmsg, errcode )! TODO (jiwon) what is this '1'? + ! if (errcode /= 0) return + + ! ! TODO (jiwon) deallocate height delta + ! end subroutine set_gas_species_values + +end module musica_ccpp_tuvx_gas_species_profile \ No newline at end of file diff --git a/test/musica/test_musica_api.F90 b/test/musica/test_musica_api.F90 index eb185fbc..95d73c80 100644 --- a/test/musica/test_musica_api.F90 +++ b/test/musica/test_musica_api.F90 @@ -144,8 +144,8 @@ subroutine test_chapman() implicit none - integer, parameter :: NUM_SPECIES = 5 - integer, parameter :: NUM_TUVX_CONSTITUENTS = 1 + integer, parameter :: NUM_SPECIES = 5 + integer, parameter :: NUM_TUVX_CONSTITUENTS = 1 ! This test requires that the number of grid cells = 4, which is the default ! vector dimension for MICM. This restriction will be removed once ! https://github.com/NCAR/musica/issues/217 is finished. @@ -235,11 +235,11 @@ subroutine test_chapman() 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.0319988_kind_phys .and. is_advected) .or. & + tmp_bool = (trim(species_name) == "O2" .and. molar_mass == 0.0319988_kind_phys .and. .not. is_advected) .or. & (trim(species_name) == "O" .and. molar_mass == 0.0159994_kind_phys .and. .not. is_advected) .or. & (trim(species_name) == "O1D" .and. molar_mass == 0.0159994_kind_phys .and. .not. is_advected) .or. & (trim(species_name) == "O3" .and. molar_mass == 0.0479982_kind_phys .and. is_advected) .or. & - (trim(species_name) == "N2" .and. molar_mass == 0.0280134_kind_phys .and. is_advected) .or. & + (trim(species_name) == "N2" .and. molar_mass == 0.0280134_kind_phys .and. .not. is_advected) .or. & (trim(species_name) == "cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water" .and. & molar_mass == 0.018_kind_phys .and. is_advected) ASSERT(tmp_bool) @@ -372,8 +372,8 @@ subroutine test_terminator() implicit none - integer, parameter :: NUM_SPECIES = 2 - integer, parameter :: NUM_TUVX_CONSTITUENTS = 1 + integer, parameter :: NUM_SPECIES = 4 + integer, parameter :: NUM_TUVX_CONSTITUENTS = 1 ! This test requires that the number of grid cells = 4, which is the default ! vector dimension for MICM. This restriction will be removed once ! https://github.com/NCAR/musica/issues/217 is finished. @@ -412,7 +412,7 @@ subroutine test_terminator() character(len=:), allocatable :: micm_species_name logical :: tmp_bool, is_advected integer :: i, j, k - integer :: Cl_index, Cl2_index + integer :: Cl_index, Cl2_index, O2_index, O3_index real(kind_phys) :: total_Cl, total_Cl_init call get_wavelength_edges(photolysis_wavelength_grid_interfaces) @@ -465,8 +465,10 @@ subroutine test_terminator() ASSERT(errcode == 0) tmp_bool = (trim(species_name) == "Cl" .and. molar_mass == 0.035453_kind_phys .and. is_advected) .or. & (trim(species_name) == "Cl2" .and. molar_mass == 0.070906_kind_phys .and. is_advected) .or. & - (trim(species_name) == "cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water" & - .and. molar_mass == 0.018_kind_phys .and. is_advected) + (trim(species_name) == "O2" .and. molar_mass == 0.0319988_kind_phys .and. .not. is_advected) .or. & + (trim(species_name) == "O3" .and. molar_mass == 0.0479982_kind_phys .and. .not. is_advected) .or. & + (trim(species_name) == "cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water" & + .and. molar_mass == 0.018_kind_phys .and. is_advected) ASSERT(tmp_bool) call constituent_props(i)%units(units, errcode, errmsg) if (errcode /= 0) then @@ -501,6 +503,12 @@ subroutine test_terminator() else if (micm_species_name == "Cl2") then Cl2_index = i base_conc = 1.0e-6_kind_phys + else if (micm_species_name == "O2") then + O2_index = i + base_conc = 0.21_kind_phys + else if (micm_species_name == "O3") then + O3_index = i + base_conc = 1.0e-4_kind_phys else write(*,*) "Unknown species: ", micm_species_name stop 3 From dbf019c24756e9717c3e3e60d9acf5818750febe Mon Sep 17 00:00:00 2001 From: Jiwon Gim Date: Wed, 4 Dec 2024 11:35:57 -0700 Subject: [PATCH 08/14] stash for profile group --- schemes/musica/musica_ccpp.F90 | 28 +- schemes/musica/tuvx/musica_ccpp_tuvx.F90 | 51 +- .../musica_ccpp_tuvx_gas_species_profile.F90 | 460 +++++++----------- 3 files changed, 226 insertions(+), 313 deletions(-) diff --git a/schemes/musica/musica_ccpp.F90 b/schemes/musica/musica_ccpp.F90 index 9d922226..a93c202e 100644 --- a/schemes/musica/musica_ccpp.F90 +++ b/schemes/musica/musica_ccpp.F90 @@ -108,20 +108,6 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co number_of_rate_parameters) :: rate_parameters ! various units integer :: i_elem - ! Calculate photolysis rate constants using TUV-x - call tuvx_run(temperature, dry_air_density, & - geopotential_height_wrt_surface_at_midpoint, & - geopotential_height_wrt_surface_at_interface, & - surface_geopotential, surface_temperature, & - surface_albedo, & - number_of_photolysis_wavelength_grid_sections, & - photolysis_wavelength_grid_interfaces, & - extraterrestrial_flux, & - standard_gravitational_acceleration, & - cloud_area_fraction, constituents, & - air_pressure_thickness, rate_parameters, & - errmsg, errcode) - ! TODO(jiwon) this might not be correct because it doesn't know the index ! Get the molar mass that is set in the call to instantiate() do i_elem = 1, size(molar_mass_arr) @@ -145,6 +131,20 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co ! Convert CAM-SIMA unit to MICM unit (kg kg-1 -> mol m-3) call convert_to_mol_per_cubic_meter(dry_air_density, molar_mass_arr, constituents) + ! Calculate photolysis rate constants using TUV-x + call tuvx_run(temperature, dry_air_density, & + geopotential_height_wrt_surface_at_midpoint, & + geopotential_height_wrt_surface_at_interface, & + surface_geopotential, surface_temperature, & + surface_albedo, & + number_of_photolysis_wavelength_grid_sections, & + photolysis_wavelength_grid_interfaces, & + extraterrestrial_flux, & + standard_gravitational_acceleration, & + cloud_area_fraction, constituents, & + air_pressure_thickness, rate_parameters, & + errmsg, errcode) + ! Solve chemistry at the current time step call micm_run(time_step, temperature, pressure, dry_air_density, rate_parameters, & constituents, errmsg, errcode) diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 index 2805d4fa..688165bf 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 @@ -7,6 +7,7 @@ module musica_ccpp_tuvx use musica_ccpp_util, only: has_error_occurred use musica_tuvx, only: tuvx_t, grid_t, profile_t, radiator_t use musica_util, only: mappings_t, index_mappings_t + use musica_ccpp_tuvx_gas_species_profile, only: gas_species_t, profile_group_t implicit none private @@ -21,6 +22,9 @@ module musica_ccpp_tuvx type(profile_t), pointer :: extraterrestrial_flux_profile => null() type(radiator_t), pointer :: cloud_optics => null() type(index_mappings_t), pointer :: photolysis_rate_constants_mapping => null( ) + type(gas_species_t), allocatable :: gas_species_group(:) + type(profile_t), pointer :: profile_gas_species_group(:) + integer, parameter :: DEFAULT_NUM_PHOTOLYSIS_RATE_CONSTANTS = 0 integer :: number_of_photolysis_rate_constants = DEFAULT_NUM_PHOTOLYSIS_RATE_CONSTANTS integer, parameter :: DEFAULT_INDEX_NOT_FOUND = -1 @@ -50,6 +54,8 @@ end subroutine reset_tuvx_map_state !> Deallocates objects associated with TUV-x subroutine cleanup_tuvx_resources() + use musica_ccpp_tuvx_gas_species_profile, & + only: deallocate_gas_species_profile_group if (associated( height_grid )) then deallocate( height_grid ) @@ -86,6 +92,8 @@ subroutine cleanup_tuvx_resources() photolysis_rate_constants_mapping => null() end if + call deallocate_gas_species_profile_group( profile_gas_species_group ) + end subroutine cleanup_tuvx_resources !> Registers constituent properties with the CCPP needed by TUV-x @@ -142,8 +150,8 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & extraterrestrial_flux_unit use musica_ccpp_tuvx_cloud_optics, & only: create_cloud_optics_radiator, cloud_optics_label - use musica_ccpp_tuvx_gas_species_profilev, & - only: gas_species_t, configure_gas_species, create_gas_species_profile_group + use musica_ccpp_tuvx_gas_species_profile, & + only: configure_gas_species, create_gas_species_profile_group integer, intent(in) :: vertical_layer_dimension ! (count) integer, intent(in) :: vertical_interface_dimension ! (count) @@ -160,8 +168,7 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & type(radiator_map_t), pointer :: radiators type(configuration_t) :: config type(mappings_t), pointer :: photolysis_rate_constants_ordering - type(gas_species_t), allocatable :: gas_species_group(:) - type(profile_t), pointer :: profile_gas_species_group(:) + integer :: i_species ! Get needed indices in constituents array call ccpp_const_get_idx(constituent_props, CLOUD_LIQUID_WATER_CONTENT_LABEL, & @@ -365,6 +372,20 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & return end if + allocate( profile_gas_species_group( number_of_gas_species ) ) + do i_species = 1, size(gas_species_group) + write(*,*) " [init] profiles%get( gas_species_group(i_species)%name: ", gas_species_group(i_species)%name + profile_gas_species_group(i_species) = & + profiles%get( gas_species_group(i_species)%name, gas_species_group(i_species)%unit, error ) + if (has_error_occurred( error, errmsg, errcode )) then + deallocate( tuvx ) + tuvx => null() + call reset_tuvx_map_state( grids, profiles, null() ) + call cleanup_tuvx_resources() + return + end if + end do + radiators => tuvx%get_radiators( error ) if (has_error_occurred( error, errmsg, errcode )) then deallocate( tuvx ) @@ -440,7 +461,7 @@ subroutine tuvx_run(temperature, dry_air_density, & 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 - + use musica_ccpp_tuvx_gas_species_profile, only: set_gas_species_values real(kind_phys), intent(in) :: temperature(:,:) ! K (column, layer) real(kind_phys), intent(in) :: dry_air_density(:,:) ! kg m-3 (column, layer) real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_midpoint(:,:) ! m (column, layer) @@ -469,7 +490,7 @@ subroutine tuvx_run(temperature, dry_air_density, & real(kind_phys) :: solar_zenith_angle ! degrees real(kind_phys) :: earth_sun_distance ! AU type(error_t) :: error - integer :: i_col, i_level + integer :: i_col, i_level, i_gas_species reciprocal_of_gravitational_acceleration = 1.0_kind_phys / standard_gravitational_acceleration @@ -484,6 +505,7 @@ 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), & @@ -504,6 +526,20 @@ subroutine tuvx_run(temperature, dry_air_density, & errmsg, errcode ) if (errcode /= 0) return + !!! THIS DOESN'T WORK......... + do i_gas_species = 1, size(gas_species_group) + write(*,*) "gas_species_group(i_gas_species)%index_constituent_props)", gas_species_group(i_gas_species)%index_constituent_props + write(*,*) "gas_species_group(i_gas_species)%index_constituent_props)", gas_species_group(i_gas_species)%name + write(*,*) "size(gas_species_group)", gas_species_group(1)%name + + call set_gas_species_values(profile_gas_species_group(i_gas_species), & + gas_species_group(i_gas_species), & + constituents(i_col,:,gas_species_group(i_gas_species)%index_constituent_props), & ! mol m-3 + errmsg, errcode) + if (errcode /= 0) return + end do + !!! + ! temporary values until these are available from the host model solar_zenith_angle = 0.0_kind_phys earth_sun_distance = 1.0_kind_phys @@ -530,12 +566,15 @@ end subroutine tuvx_run !> Finalizes TUV-x subroutine tuvx_final(errmsg, errcode) + use musica_ccpp_tuvx_gas_species_profile, only: deallocate_gas_species + character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode errmsg = '' errcode = 0 + call deallocate_gas_species( gas_species_group ) call cleanup_tuvx_resources() if (associated( tuvx )) then diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 index a5e99220..b037b2a7 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 @@ -4,7 +4,8 @@ module musica_ccpp_tuvx_gas_species_profile implicit none private - public :: create_gas_species, configure_gas_species, create_gas_species_profile_group + public :: create_gas_species, configure_gas_species, create_gas_species_profile_group, & + deallocate_gas_species, deallocate_gas_species_profile_group, set_gas_species_values !> Conversion factor from km to cm real(kind_phys), parameter, public :: km_to_cm = 1.0e5 @@ -15,15 +16,20 @@ module musica_ccpp_tuvx_gas_species_profile !> Number of wavelength grid bins integer, protected :: num_midpoint_heights_ = DEFAULT_NUM_GRID_SECTIONS integer, protected :: num_interfaces_ = DEFAULT_NUM_GRID_SECTIONS + real(kind_phys), allocatable :: height_delta_(:) type, public :: gas_species_t character(len=50) :: name character(len=20) :: unit real(kind_phys) :: molar_mass - real(kind_phys) :: scale_height + real(kind_phys) :: scale_height ! km integer :: index_constituent_props end type gas_species_t + type, public :: profile_group_t + type(profile_t), pointer :: ptr + end type gas_species_t + contains function create_gas_species(name, unit, molar_mass, scale_height, & @@ -67,7 +73,7 @@ subroutine configure_gas_species(constituent_props, gas_species_group, errmsg, e real(kind_phys) :: molar_mass_O3 ! kg mol-1 integer :: index_air, index_O2, index_O3 - allocate(gas_species_group(number_of_gas_species)) ! TODO(jiwon) - when deallocate this? + allocate(gas_species_group(number_of_gas_species)) ! Air ! Find the index of the species and create a 'gas_species_t' for it @@ -83,7 +89,6 @@ subroutine configure_gas_species(constituent_props, gas_species_group, errmsg, e MOLAR_MASS_DRY_AIR, SCALE_HEIGHT_DRY_AIR, index_air) ! O2 call ccpp_const_get_idx(constituent_props, "O2", index_O2, errmsg, errcode) - write(*,*) "air index: ", index_O2 if (errcode /= 0) return call constituent_props(index_O2)%molar_mass(molar_mass_O2, errcode, errmsg) @@ -112,11 +117,11 @@ subroutine create_gas_species_profile_group(height_grid, gas_species_group, & use musica_tuvx_profile, only: profile_t use musica_util, only: error_t - type(grid_t), intent(in) :: height_grid - type(gas_species_t), intent(in) :: gas_species_group(:) - type(profile_t), pointer, intent(out) :: profile_group(:) - character(len=512), intent(out) :: errmsg - integer, intent(out) :: errcode + type(grid_t), intent(in) :: height_grid + type(gas_species_t), intent(in) :: gas_species_group(:) + type(profile_t), pointer, intent(out) :: profile_group(:) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode ! local variables type(error_t) :: error @@ -127,329 +132,198 @@ subroutine create_gas_species_profile_group(height_grid, gas_species_group, & allocate( profile_group( number_of_gas_species ) ) do i_species = 1, number_of_gas_species - allocate(profile_group(i_species)) write(*,*) "gas_species_group(i_species)%name: ", gas_species_group(i_species)%name profile_group(i_species) = profile_t( gas_species_group(i_species)%name, & gas_species_group(i_species)%unit, height_grid, error ) if ( has_error_occurred( error, errmsg, errcode ) ) return + ! profile_group(i_species) = ptr ! TODO(jiwon) deallocate the previously allocated one? end do - end subroutine create_gas_species_profile_group - - ! function get_height_delta(height_grid, num_interfaces, errmsg, errcode) result( height_delta ) - ! 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(in) :: height_grid - ! integer, intent(in) :: num_interfaces - ! character(len=*), intent(out) :: errmsg - ! integer, intent(out) :: errcode - - ! ! local variables - ! type(error_t) :: error - ! real(kind_phys) :: interfaces(num_interfaces) - ! real(kind_phys) :: height_delta(num_interfaces-1) - ! integer :: i_elem - - ! call height_grid%get_edges(interfaces, error) - ! if ( has_error_occurred( error, errmsg, errcode ) ) return - - ! do i_elem = 1, size(height_delta) - ! height_delta[i] = interfaces[i+1] - interfaces[i] - ! end do - - ! end function get_height_delta - - ! !> Creates a TUV-x surface albedo profile from the host-model wavelength grid - ! function create_air_profile( height_grid, errmsg, errcode ) & - ! result( profile ) - ! use musica_ccpp_util, only: has_error_occurred - ! use musica_tuvx_grid, only: grid_t - ! use musica_tuvx_profile, only: profile_t - ! use musica_util, only: error_t - - ! type(grid_t), intent(in) :: height_grid - ! character(len=*), intent(out) :: errmsg - ! integer, intent(out) :: errcode - ! type(profile_t), pointer :: profile - - ! ! local variables - ! type(error_t) :: error - - ! profile => profile_t( air_label, air_unit, height_grid, error ) - ! if ( has_error_occurred( error, errmsg, errcode ) ) return - - ! end function create_air_profile - - ! !> Creates a TUV-x surface albedo profile from the host-model wavelength grid - ! function create_O2_profile( height_grid, errmsg, errcode ) & - ! result( profile ) - ! use musica_ccpp_util, only: has_error_occurred - ! use musica_tuvx_grid, only: grid_t - ! use musica_tuvx_profile, only: profile_t - ! use musica_util, only: error_t - - ! type(grid_t), intent(in) :: height_grid - ! character(len=*), intent(out) :: errmsg - ! integer, intent(out) :: errcode - ! type(profile_t), pointer :: profile - - ! ! local variables - ! type(error_t) :: error - - ! profile => profile_t( O2_label, O2_unit, height_grid, error ) - ! if ( has_error_occurred( error, errmsg, errcode ) ) return - - ! end function create_O2_profile - - ! !> Creates a TUV-x surface albedo profile from the host-model wavelength grid - ! function create_O3_profile( height_grid, errmsg, errcode ) & - ! result( profile ) - ! use musica_ccpp_util, only: has_error_occurred - ! use musica_tuvx_grid, only: grid_t - ! use musica_tuvx_profile, only: profile_t - ! use musica_util, only: error_t - - ! type(grid_t), intent(in) :: height_grid - ! character(len=*), intent(out) :: errmsg - ! integer, intent(out) :: errcode - ! type(profile_t), pointer :: profile - - ! ! local variables - ! type(error_t) :: error - - ! profile => profile_t( O3_label, O3_unit, height_grid, error ) - ! if ( has_error_occurred( error, errmsg, errcode ) ) return - - ! end function create_O3_profile - - ! subroutine create_gas_species_profiles( air_profile, O2_profile, O3_profile, & - ! height_grid, errmsg, errcode ) - ! use musica_ccpp_util, only: has_error_occurred - ! use musica_tuvx_grid, only: grid_t - ! use musica_tuvx_profile, only: profile_t - ! use musica_util, only: error_t - - ! type(profile_t), intent(inout) :: air_profile - ! type(profile_t), intent(inout) :: O2_profile - ! type(profile_t), intent(inout) :: O3_profile - ! type(grid_t), intent(in) :: height_grid - ! character(len=*), intent(out) :: errmsg - ! integer, intent(out) :: errcode - - ! ! local variables - ! type(error_t) :: error - - ! air_profile = create_air_profile( height_grid, errmsg, errcode ) - ! if (errcode /= 0) return - - ! O2_profile = create_O2_profile( height_grid, errmsg, errcode ) - ! if (errcode /= 0) return + num_midpoint_heights_ = height_grid%number_of_sections( error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return - ! O3_profile = create_O3_profile( height_grid, errmsg, errcode ) - ! if (errcode /= 0) return + num_interfaces_ = num_midpoint_heights_ + 1 - ! num_midpoint_heights_ = height_grid%number_of_sections( error ) - ! if ( has_error_occurred( error, errmsg, errcode ) ) return + ! Set height delta values + height_delta_ = get_height_delta( height_grid, num_interfaces_, errmsg, errcode ) + if ( has_error_occurred( error, errmsg, errcode ) ) return - ! num_interfaces_ = num_midpoint_heights_ + 1 + ! do i_species = 1, number_of_gas_species + ! deallocate( profile_group(i_species)) + ! end do - ! ! Set height delta values - ! height_delta_ = get_height_delta( height_grid, num_interfaces_, errmsg, errcode ) - ! if ( has_error_occurred( error, errmsg, errcode ) ) return + end subroutine create_gas_species_profile_group - ! end subroutine create_gas_species_profiles + subroutine deallocate_gas_species(gas_species_group) + type(gas_species_t), allocatable, intent(inout) :: gas_species_group(:) - ! !> Sets TUV-x surface albedo values - ! subroutine set_dry_air_density_values( profile, dry_air_density, errmsg, errcode ) - ! use ccpp_kinds, only: kind_phys - ! use musica_ccpp_util, only: has_error_occurred - ! use musica_tuvx_profile, only: profile_t - ! use musica_util, only: error_t + if (allocated( gas_species_group )) then + write(*,*) "deallocate( gas_species_group )" + deallocate( gas_species_group ) + end if - ! type(profile_t), intent(inout) :: profile - ! real(kind_phys), intent(in) :: dry_air_density(:) ! TODO(jiwon) : molecule cm-3 for CAM, kg m-3 as input - ! integer, intent(in) :: index_air - ! character(len=*), intent(out) :: errmsg - ! integer, intent(out) :: errcode + end subroutine - ! ! local variables - ! type(error_t) :: error - ! real(kind_phys) :: scale_height = 8.01_kind_phys ! km - ! integer :: num_midpoints = size(dry_air_density) - ! real(kind_phys) :: interface(num_midpoints + 2) - ! real(kind_phys) :: densities(num_midpoints + 1) - - ! interfaces(1) = dry_air_density(num_midpoints) - ! interfaces(2:num_midpoints+1) = dry_air_density(num_midpoints:1:-1) - ! interfaces(num_midpoints+2) = dry_air_density(1) ! use upper mid-point value for top edge - - ! densities(1:num_midpoints+1) = & - ! height_delta_(1:num_midpoints+1) * km_to_cm & ! TODO(jiwon) : molecule cm-3 for CAM, kg m-3 as input - ! * sqrt(interfaces(1:num_midpoints+1)) & - ! * sqrt(interfaces(2:num_midpoints+2)) - - ! call profile%set_edge_values( interfaces, error ) - ! if ( has_error_occurred( error, errmsg, errcode ) ) return - - ! call profile%set_layer_densities( densities, error) - ! if ( has_error_occurred( error, errmsg, errcode ) ) return - - ! call profile%calculate_exo_layer_density( scale_height, error ) - ! if ( has_error_occurred( error, errmsg, errcode ) ) return - - ! end subroutine set_air_density_values - - ! !> Sets TUV-x surface albedo values - ! subroutine set_O2_values(profile, O2_constituent, height_delta, errmsg, errcode) - ! use ccpp_kinds, only: kind_phys - ! use musica_ccpp_util, only: has_error_occurred + ! subroutine deallocate_gas_species_profile_group(profile_group) ! use musica_tuvx_profile, only: profile_t - ! use musica_util, only: error_t - ! type(profile_t), intent(inout) :: profile - ! real(kind_phys), intent(in) :: O2_constituent(:) ! mol m-3 - ! real(kind_phys), intent(in) :: height_delta(:) ! km - ! character(len=*), intent(out) :: errmsg - ! integer, intent(out) :: errcode + ! type(profile_t), allocatable, intent(inout) :: profile_group(:) ! ! local variables - ! type(error_t) :: error - ! real(kind_phys), parameter :: SCALE_HEIGHT = 7.0_kind_phys ! km - ! real(kind_phys) :: O2_constituent_mol_cm-3(size(O2_constituent)) ! mol cm-3 - ! integer :: num_midpoints = size(O2_constituent) - ! real(kind_phys) :: interface(size(O2_constituent) + 2) - ! real(kind_phys) :: densities(size(O2_constituent) + 1) + ! integer :: i_species + ! type(profile_t), pointer :: ptr + + ! if ( allocated(profile_group) ) then + ! write(*,*) "allocated(profile_group)" + ! ! Deallocate each element of the pointer array + ! do i_species = 1, size(profile_group) + + ! ! ptr = profile_group(i_species) + ! if ( allocated( profile_group(i_species)) ) then + ! deallocate( profile_group(i_species) ) + ! write(*,*) "deallocate( profile_group(i_species) ) ", i_species + ! ! if ( associated( ptr) ) then + ! ! write(*,*) "deallocate( profile_group(i_species) ) ", i_species + ! ! deallocate( ptr) + ! end if + ! end do + + ! ! Deallocate the entire profile_group array + ! ! if (allocated(profile_group)) then + ! ! nullify(profile_group) + ! ! end if - ! O2_constituent_mol_cm-3(:) = O2_constituent(:) * m_3_to_cm_3 + ! else + ! write(*,*) "allocated(profile_group) is false" + ! return + ! end if - ! ! TODO(jiwon) if ( is_parameterized ) then + ! end subroutine - ! interfaces(1) = O2_constituent_mol_cm-3(num_midpoints) - ! interfaces(2:num_midpoints+1) = O2_constituent_mol_cm-3(num_midpoints:1:-1) - ! interfaces(num_midpoints+2) = O2_constituent_mol_cm-3(1) - ! densities(1:num_midpoints+1) = & - ! height_delta(1:num_midpoints+1) * km_to_cm & - ! * sqrt(interfaces(1:num_midpoints+1)) & - ! * sqrt(interfaces(2:num_midpoints+2)) + subroutine deallocate_gas_species_profile_group(profile_group) + use musica_util, only: error_t + use musica_tuvx_profile, only: profile_t + + type(profile_t), pointer, intent(inout) :: profile_group(:) + character(len=512) :: errmsg + integer :: errcode + + ! type(error_t) :: error + integer :: i_species + logical :: deallocation_error + type(profile_t), pointer :: ptr + + write(*,*) "Called deallocate_gas_species_profile_group(profile_group)" + deallocation_error = .false. + + ! Check if the profile_group is allocated + ! if (.not. allocated(profile_group)) then + ! write(*,*) 'profile_group is not allocated.' + ! errcode = 1 + ! return + ! end if + + ! Deallocate each element of profile_group (array of pointers) + ! do i_species = 1, size(profile_group) + ! if (allocated( profile_group(i_species)) ) then + ! ! Deallocate each individual pointer + ! deallocate(profile_group(i_species), stat=errcode) + ! if (errcode /= 0) then + ! ! errmsg = 'Error deallocating profile_group(' // trim(adjustl(i_species)) // ')' + ! deallocation_error = .true. + ! exit + ! end if + ! end if + ! end do + deallocate(profile_group, stat=errcode) + !do i_species = 1, size(profile_group) + !ptr => profile_group(i_species) + ! deallocate(profile_group(i_species), stat=errcode) + !end do + ! write(*,*) "size:L profile_group ", size(profile_group) + ! ! Finally, deallocate the array of pointers itself + ! if (.not. deallocation_error) then + ! deallocate(profile_group, stat=errcode) + ! if (errcode /= 0) then + ! write(*,*) 'Error deallocating the profile_group array.' + ! else + ! write(*,*) 'profile_group successfully deallocated.' + ! end if + !end if + + end subroutine deallocate_gas_species_profile_group - ! call profile%set_edge_values( interfaces, error ) - ! if ( has_error_occurred( error, errmsg, errcode ) ) return + function get_height_delta(height_grid, num_interfaces, errmsg, errcode) result( height_delta ) + use musica_ccpp_util, only: has_error_occurred + use musica_tuvx_grid, only: grid_t + use musica_util, only: error_t - ! call profile%set_layer_densities( densities, error) - ! if ( has_error_occurred( error, errmsg, errcode ) ) return + type(grid_t), intent(in) :: height_grid + integer, intent(in) :: num_interfaces + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errcode - ! call profile%calculate_exo_layer_density( SCALE_HEIGHT, error ) - ! if ( has_error_occurred( error, errmsg, errcode ) ) return + ! local variables + type(error_t) :: error + real(kind_phys) :: interfaces(num_interfaces) + real(kind_phys) :: height_delta(num_interfaces-1) + integer :: i_elem - ! end subroutine set_O2_values + call height_grid%get_edges(interfaces, error) + if ( has_error_occurred( error, errmsg, errcode ) ) return - ! !> Sets TUV-x surface albedo values - ! subroutine set_O3_values(profile, O3_constituent, height_delta, above_column_density, & - ! errmsg, errcode) - ! use ccpp_kinds, only: kind_phys - ! use musica_ccpp_util, only: has_error_occurred - ! use musica_tuvx_profile, only: profile_t - ! use musica_util, only: error_t + do i_elem = 1, size(height_delta) + height_delta(i_elem) = interfaces(i_elem+1) - interfaces(i_elem) + end do - ! type(profile_t), intent(inout) :: profile - ! real(kind_phys), intent(in) :: O2_constituent(:) ! mol m-3 - ! real(kind_phys), intent(in) :: height_delta(:) ! km - ! real(kind_phys), intent(in) :: above_column_density(:,:) ! molecule cm-2 - ! character(len=*), intent(out) :: errmsg - ! integer, intent(out) :: errcode + end function get_height_delta - ! ! local variables - ! type(error_t) :: error - ! real(kind_phys), parameter :: SCALE_HEIGHT = 7.0_kind_phys ! km - ! real(kind_phys) :: O2_constituent_mol_cm-3(size(O3_constituent)) ! mol cm-3 - ! integer :: num_midpoints = size(O3_constituent) - ! real(kind_phys) :: interface(size(O3_constituent) + 2) - ! real(kind_phys) :: densities(size(O3_constituent) + 1) - - ! O3_constituent_mol_cm-3(:) = O3_constituent(:) * m_3_to_cm_3 !TODO(devided maybe) - - ! interfaces(1) = O3_constituent_mol_cm-3(num_midpoints) - ! interfaces(2:num_midpoints+1) = O3_constituent_mol_cm-3(num_midpoints:1:-1) - ! interfaces(num_midpoints+2) = O3_constituent_mol_cm-3(1) - - ! call profile%set_edge_values( interfaces, error ) - ! if ( has_error_occurred( error, errmsg, errcode ) ) return - - ! densities(1:num_midpoints+1) = & - ! height_delta(1:num_midpoints+1) * km_to_cm & - ! * sqrt(interfaces(1:num_midpoints+1)) & - ! * sqrt(interfaces(2:num_midpoints+2)) - - ! if ( num_absorbing_species >= 1 ) then ! TODO - nabscol = 2 in chem_mods - ! densities(1) = 0.5_kind_phys * above_column_density(num_midpoints, 1) - ! densities(2:num_midpoints) = 0.5_kind_phys * ( above_column_density(num_midpoints-1:1:-1,1) & - ! + above_column_density(num_midpoints:2:-1,1) ) - ! densities(num_midpoints+1) = above_column_density(0,1) & - ! + 0.5_kind_phys * above_column_density(1,1) - - ! call profile%set_layer_densities( densities, error) - ! if ( has_error_occurred( error, errmsg, errcode ) ) return - - ! call profile%set_exo_layer_density( above_column_density(0, 1), error) ! TODO: should be incorrect? exo_density = above_column_density(0,1) ) - ! if ( has_error_occurred( error, errmsg, errcode ) ) return - ! else - ! densities(1:num_midpoints+1) = height_delta(1:num_midpoints+1) * km_to_cm & - ! * ( interfaces(1:num_midpoints+1) + interfaces(2:num_midpoints+2) ) * 0.5_kind_phys + !> Sets TUV-x surface albedo values + subroutine set_gas_species_values(profile, gas_species, constituent, errmsg, errcode) + use ccpp_kinds, only: kind_phys + use musica_ccpp_util, only: has_error_occurred + use musica_tuvx_profile, only: profile_t + use musica_util, only: error_t - ! call profile%set_layer_densities( densities, error) - ! if ( has_error_occurred( error, errmsg, errcode ) ) return + type(profile_t), intent(inout) :: profile + type(gas_species_t), intent(in) :: gas_species + real(kind_phys), intent(in) :: constituent(:) ! mol m-3 + ! real(kind_phys), intent(in) :: height_delta(:) ! km + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errcode - ! call profile%calculate_exo_layer_density( SCALE_HEIGHT, error ) - ! if ( has_error_occurred( error, errmsg, errcode ) ) return - ! end if + ! local variables + type(error_t) :: error + integer :: num_constituents + real(kind_phys) :: constituent_mol_per_cm_3(size(constituent)) ! mol cm-3 + real(kind_phys) :: interfaces(size(constituent) + 2) + real(kind_phys) :: densities(size(constituent) + 1) - ! end subroutine set_O3_values + num_constituents = size(constituent) + constituent_mol_per_cm_3(:) = constituent(:) / m_3_to_cm_3 - ! subroutine set_gas_species_values( air_profile, O2_profile, O3_profile & - ! fixed_species_density, species_volume_mixing_ratio, above_column_density, & - ! index_air, index_O2, index_O3, errmsg, errcode) - ! use ccpp_kinds, only: kind_phys - ! use musica_tuvx_profile, only: profile_t - ! use musica_util, only: error_t - - ! type(profile_t), intent(inout) :: air_profile - ! type(profile_t), intent(inout) :: O2_profile - ! type(profile_t), intent(inout) :: O3_profile - ! real(kind_phys), intent(in) :: fixed_species_density(:,:) ! molecule cm-3 - ! real(kind_phys), intent(in) :: species_volume_mixing_ratio(:,:) ! mol mol-1 - ! real(kind_phys), intent(in) :: above_column_density(:,:) ! molecule cm-2 - ! integer, intent(in) :: index_air - ! integer, intent(in) :: index_O2 - ! integer, intent(in) :: index_O3 - ! character(len=*), intent(out) :: errmsg - ! integer, intent(out) :: errcode + ! TODO(jiwon) if ( is_parameterized ) then - ! ! local variables - ! real(kind_phys) :: air_density = fixed_species_density(:, index_air) + interfaces(1) = constituent_mol_per_cm_3(num_constituents) + interfaces(2:num_constituents+1) = constituent_mol_per_cm_3(num_constituents:1:-1) + interfaces(num_constituents+2) = constituent_mol_per_cm_3(1) + densities(1:num_constituents+1) = & + height_delta_(1:num_constituents+1) * km_to_cm & + * sqrt(interfaces(1:num_constituents+1)) & + * sqrt(interfaces(2:num_constituents+2)) - ! ! do the conversion to mol / cm3 - ! call set_air_density_values( air_profile, & - ! fixed_species_density(:, index_air), errmsg, errcode ) - ! if (errcode /= 0) return + call profile%set_edge_values( interfaces, error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return - ! call set_O2_values( O2_profile, fixed_species_density(:, index_air), & - ! fixed_species_density(:, index_O2), & - ! species_volume_mixing_ratio(:, index_O2), & - ! errmsg, errcode ) - ! if (errcode /= 0) return + call profile%set_layer_densities( densities, error) + if ( has_error_occurred( error, errmsg, errcode ) ) return - ! call set_O3_values( O3_profile, fixed_species_density(:, index_air), & - ! fixed_species_density(:, index_O3), & - ! species_volume_mixing_ratio(:, index_O3), & - ! above_column_density(:, 1), errmsg, errcode )! TODO (jiwon) what is this '1'? - ! if (errcode /= 0) return + call profile%calculate_exo_layer_density( gas_species%scale_height, error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return - ! ! TODO (jiwon) deallocate height delta - ! end subroutine set_gas_species_values + end subroutine set_gas_species_values end module musica_ccpp_tuvx_gas_species_profile \ No newline at end of file From 0ed49958c86a545b62b20af07ba7ca0c50d0114f Mon Sep 17 00:00:00 2001 From: Jiwon Gim Date: Wed, 4 Dec 2024 18:54:13 -0700 Subject: [PATCH 09/14] fix bugs - test passed --- schemes/musica/tuvx/musica_ccpp_tuvx.F90 | 81 +++-- .../musica_ccpp_tuvx_gas_species_profile.F90 | 278 ++++++------------ test/musica/test_musica_api.F90 | 5 + 3 files changed, 154 insertions(+), 210 deletions(-) diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 index 688165bf..8077c33f 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 @@ -22,8 +22,10 @@ module musica_ccpp_tuvx type(profile_t), pointer :: extraterrestrial_flux_profile => null() type(radiator_t), pointer :: cloud_optics => null() type(index_mappings_t), pointer :: photolysis_rate_constants_mapping => null( ) - type(gas_species_t), allocatable :: gas_species_group(:) - type(profile_t), pointer :: profile_gas_species_group(:) + + real(kind_phys), allocatable :: height_delta(:) ! tuvx_run needs this + type(gas_species_t), allocatable :: gas_species_group(:) + type(profile_group_t), allocatable :: profile_gas_species_group(:) integer, parameter :: DEFAULT_NUM_PHOTOLYSIS_RATE_CONSTANTS = 0 integer :: number_of_photolysis_rate_constants = DEFAULT_NUM_PHOTOLYSIS_RATE_CONSTANTS @@ -92,8 +94,11 @@ subroutine cleanup_tuvx_resources() photolysis_rate_constants_mapping => null() end if - call deallocate_gas_species_profile_group( profile_gas_species_group ) + if (allocated( profile_gas_species_group )) then + deallocate( profile_gas_species_group ) + end if + call deallocate_gas_species_profile_group( profile_gas_species_group ) end subroutine cleanup_tuvx_resources !> Registers constituent properties with the CCPP needed by TUV-x @@ -151,7 +156,7 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & use musica_ccpp_tuvx_cloud_optics, & only: create_cloud_optics_radiator, cloud_optics_label use musica_ccpp_tuvx_gas_species_profile, & - only: configure_gas_species, create_gas_species_profile_group + only: configure_gas_species, create_gas_species_profile_group, set_height_delta integer, intent(in) :: vertical_layer_dimension ! (count) integer, intent(in) :: vertical_interface_dimension ! (count) @@ -162,13 +167,15 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & integer, intent(out) :: errcode ! local variables - type(error_t) :: error - type(grid_map_t), pointer :: grids - type(profile_map_t), pointer :: profiles - type(radiator_map_t), pointer :: radiators - type(configuration_t) :: config - type(mappings_t), pointer :: photolysis_rate_constants_ordering - integer :: i_species + type(error_t) :: error + type(grid_map_t), pointer :: grids + type(profile_map_t), pointer :: profiles + type(radiator_map_t), pointer :: radiators + type(configuration_t) :: config + type(mappings_t), pointer :: photolysis_rate_constants_ordering + real(kind_phys) :: interfaces_height(vertical_interface_dimension) + integer :: number_of_gas_species + integer :: i_species ! Get needed indices in constituents array call ccpp_const_get_idx(constituent_props, CLOUD_LIQUID_WATER_CONTENT_LABEL, & @@ -266,6 +273,17 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & return end if + ! Gets interfaces of height to set height delta values + ! allocate( interfaces_height( vertical_interface_dimension) ) + call height_grid%get_edges(interfaces_height, error) + + allocate( height_delta( vertical_interface_dimension - 1) ) + call set_height_delta(interfaces_height, height_delta) + ! deallocate( interfaces_height ) + + number_of_gas_species = size(gas_species_group) + allocate( profile_gas_species_group( number_of_gas_species ) ) + call create_gas_species_profile_group( height_grid, & gas_species_group, profile_gas_species_group, errmsg, errcode) if (errcode /= 0) then @@ -274,6 +292,15 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & return endif + do i_species = 1, number_of_gas_species + call profiles%add( profile_gas_species_group(i_species)%profile, error ) + if (has_error_occurred( error, errmsg, errcode )) then + call reset_tuvx_map_state( grids, profiles, null() ) + call cleanup_tuvx_resources() + return + end if + end do + radiators => radiator_map_t( error ) if (has_error_occurred( error, errmsg, errcode )) then call reset_tuvx_map_state( grids, profiles, null() ) @@ -372,11 +399,14 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & return end if - allocate( profile_gas_species_group( number_of_gas_species ) ) - do i_species = 1, size(gas_species_group) - write(*,*) " [init] profiles%get( gas_species_group(i_species)%name: ", gas_species_group(i_species)%name - profile_gas_species_group(i_species) = & - profiles%get( gas_species_group(i_species)%name, gas_species_group(i_species)%unit, error ) + if (.not. allocated( profile_gas_species_group)) then + allocate( profile_gas_species_group( number_of_gas_species ) ) + end if + + do i_species = 1, number_of_gas_species + write(*,*) "gas_species_group(i_species)%label",gas_species_group(i_species)%label + profile_gas_species_group(i_species)%profile => & + profiles%get( gas_species_group(i_species)%label, gas_species_group(i_species)%unit, error ) if (has_error_occurred( error, errmsg, errcode )) then deallocate( tuvx ) tuvx => null() @@ -526,19 +556,13 @@ subroutine tuvx_run(temperature, dry_air_density, & errmsg, errcode ) if (errcode /= 0) return - !!! THIS DOESN'T WORK......... do i_gas_species = 1, size(gas_species_group) - write(*,*) "gas_species_group(i_gas_species)%index_constituent_props)", gas_species_group(i_gas_species)%index_constituent_props - write(*,*) "gas_species_group(i_gas_species)%index_constituent_props)", gas_species_group(i_gas_species)%name - write(*,*) "size(gas_species_group)", gas_species_group(1)%name - - call set_gas_species_values(profile_gas_species_group(i_gas_species), & - gas_species_group(i_gas_species), & - constituents(i_col,:,gas_species_group(i_gas_species)%index_constituent_props), & ! mol m-3 - errmsg, errcode) + call set_gas_species_values(profile_gas_species_group(i_gas_species)%profile, & + gas_species_group(i_gas_species), & + constituents(i_col,:,gas_species_group(i_gas_species)%index_constituent_props), & ! mol m-3 + height_delta, errmsg, errcode) if (errcode /= 0) return end do - !!! ! temporary values until these are available from the host model solar_zenith_angle = 0.0_kind_phys @@ -577,6 +601,11 @@ subroutine tuvx_final(errmsg, errcode) call deallocate_gas_species( gas_species_group ) call cleanup_tuvx_resources() + ! TODO(jiwon) is it okay to deallocate this in the final call? + if (allocated( height_delta )) then + deallocate( height_delta ) + end if + if (associated( tuvx )) then deallocate( tuvx ) tuvx => null() diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 index b037b2a7..273041f8 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 @@ -1,49 +1,46 @@ module musica_ccpp_tuvx_gas_species_profile use ccpp_kinds, only: kind_phys + use musica_tuvx_profile, only: profile_t implicit none private public :: create_gas_species, configure_gas_species, create_gas_species_profile_group, & - deallocate_gas_species, deallocate_gas_species_profile_group, set_gas_species_values + set_height_delta, set_gas_species_values, deallocate_gas_species, & + deallocate_gas_species_profile_group !> Conversion factor from km to cm real(kind_phys), parameter, public :: km_to_cm = 1.0e5 !> Conversion factor from m3 to cm3 real(kind_phys), parameter, public :: m_3_to_cm_3 = 1.0e6 - !> Default value of number of wavelength grid bins - integer, parameter :: DEFAULT_NUM_GRID_SECTIONS = 0 - !> Number of wavelength grid bins - integer, protected :: num_midpoint_heights_ = DEFAULT_NUM_GRID_SECTIONS - integer, protected :: num_interfaces_ = DEFAULT_NUM_GRID_SECTIONS - real(kind_phys), allocatable :: height_delta_(:) + !> Definition of the gas species type type, public :: gas_species_t - character(len=50) :: name + character(len=50) :: label character(len=20) :: unit - real(kind_phys) :: molar_mass - real(kind_phys) :: scale_height ! km - integer :: index_constituent_props + real(kind_phys) :: molar_mass ! kg mol-1 + real(kind_phys) :: scale_height ! km + integer :: index_constituent_props ! index of the constituents end type gas_species_t + !> Defines a type to contaion a pointer to 'profile_t' type, public :: profile_group_t - type(profile_t), pointer :: ptr - end type gas_species_t + type(profile_t), pointer :: profile + end type profile_group_t contains - - function create_gas_species(name, unit, molar_mass, scale_height, & + !> Creates 'gas_species_t' object + function create_gas_species(label, unit, molar_mass, scale_height, & index_constituent_props) result( gas_species ) - use ccpp_kinds, only: kind_phys - character(len=*), intent(in) :: name + character(len=*), intent(in) :: label character(len=*), intent(in) :: unit - real(kind_phys), intent(in) :: molar_mass - real(kind_phys), intent(in) :: scale_height + real(kind_phys), intent(in) :: molar_mass ! kg mol-1 + real(kind_phys), intent(in) :: scale_height ! km integer, intent(in) :: index_constituent_props type(gas_species_t) :: gas_species - gas_species%name = name + gas_species%label = label gas_species%unit = unit gas_species%molar_mass = molar_mass gas_species%scale_height = scale_height @@ -51,9 +48,11 @@ function create_gas_species(name, unit, molar_mass, scale_height, & end function create_gas_species - !> TODO(jiwon) - write function description, error message + !> Configures gas species. + ! Note: The user can customize this function to configure gas species in the system + ! This function allocates memory for the gas species that is needed in tuvx_init and tuvx_run. + ! The user is responsible for deallocating the memory. subroutine configure_gas_species(constituent_props, gas_species_group, errmsg, errcode) - use ccpp_kinds, only: kind_phys use ccpp_const_utils, only: ccpp_const_get_idx use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t @@ -63,17 +62,22 @@ subroutine configure_gas_species(constituent_props, gas_species_group, errmsg, e integer, intent(out) :: errcode ! local variable - integer :: number_of_gas_species = 3 - ! The value of molar mass of dry air is from 'CAM-SIMA/src/utils/std_atm_profile.F90' + integer :: num_gas_species = 3 real(kind_phys), parameter :: SCALE_HEIGHT_DRY_AIR = 8.01_kind_phys ! km real(kind_phys), parameter :: SCALE_HEIGHT_O2 = 7.0_kind_phys ! km real(kind_phys), parameter :: SCALE_HEIGHT_O3 = 7.0_kind_phys ! km + ! molar mass of dry air is from 'CAM-SIMA/src/utils/std_atm_profile.F90' real(kind_phys), parameter :: MOLAR_MASS_DRY_AIR = 0.0289644_kind_phys ! kg mol-1 real(kind_phys) :: molar_mass_O2 ! kg mol-1 real(kind_phys) :: molar_mass_O3 ! kg mol-1 integer :: index_air, index_O2, index_O3 - allocate(gas_species_group(number_of_gas_species)) + allocate(gas_species_group(num_gas_species)) ! need to deallocate + + ! TODO(jiwon) - I commented out this block of code that searches for the molar mass + ! of air and instead hard-coded it using the value imported from CAM-SIMA. + ! Should we register air in the register phase like cloud liquid water content, + ! or should we include it in the configuration file, or is it good as it is? ! Air ! Find the index of the species and create a 'gas_species_t' for it @@ -81,7 +85,6 @@ subroutine configure_gas_species(constituent_props, gas_species_group, errmsg, e ! write(*,*) "air index: ", index_air ! if (errcode /= 0) return - ! TODO(jiwon) air may not be registered ! call constituent_props(index_air)%molar_mass(molar_mass_air, errcode, errmsg) ! if (errcode /= 0) return @@ -109,7 +112,7 @@ subroutine configure_gas_species(constituent_props, gas_species_group, errmsg, e end subroutine configure_gas_species - !> Creates + !> Creates an array of profile group that contains a 'profile_t' for each gas species subroutine create_gas_species_profile_group(height_grid, gas_species_group, & profile_group, errmsg, errcode) use musica_ccpp_util, only: has_error_occurred @@ -117,182 +120,51 @@ subroutine create_gas_species_profile_group(height_grid, gas_species_group, & use musica_tuvx_profile, only: profile_t use musica_util, only: error_t - type(grid_t), intent(in) :: height_grid - type(gas_species_t), intent(in) :: gas_species_group(:) - type(profile_t), pointer, intent(out) :: profile_group(:) - character(len=512), intent(out) :: errmsg - integer, intent(out) :: errcode + type(grid_t), intent(in) :: height_grid + type(gas_species_t), intent(in) :: gas_species_group(:) + type(profile_group_t), intent(inout) :: profile_group(:) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode ! local variables type(error_t) :: error - integer :: number_of_gas_species integer :: i_species - number_of_gas_species = size(gas_species_group) - allocate( profile_group( number_of_gas_species ) ) - - do i_species = 1, number_of_gas_species - write(*,*) "gas_species_group(i_species)%name: ", gas_species_group(i_species)%name - profile_group(i_species) = profile_t( gas_species_group(i_species)%name, & + do i_species = 1, size(gas_species_group) + profile_group(i_species)%profile => profile_t( gas_species_group(i_species)%label, & gas_species_group(i_species)%unit, height_grid, error ) if ( has_error_occurred( error, errmsg, errcode ) ) return - ! profile_group(i_species) = ptr - ! TODO(jiwon) deallocate the previously allocated one? end do - num_midpoint_heights_ = height_grid%number_of_sections( error ) - if ( has_error_occurred( error, errmsg, errcode ) ) return - - num_interfaces_ = num_midpoint_heights_ + 1 - - ! Set height delta values - height_delta_ = get_height_delta( height_grid, num_interfaces_, errmsg, errcode ) - if ( has_error_occurred( error, errmsg, errcode ) ) return - - ! do i_species = 1, number_of_gas_species - ! deallocate( profile_group(i_species)) - ! end do - end subroutine create_gas_species_profile_group - subroutine deallocate_gas_species(gas_species_group) - type(gas_species_t), allocatable, intent(inout) :: gas_species_group(:) - - if (allocated( gas_species_group )) then - write(*,*) "deallocate( gas_species_group )" - deallocate( gas_species_group ) - end if - - end subroutine - - ! subroutine deallocate_gas_species_profile_group(profile_group) - ! use musica_tuvx_profile, only: profile_t - - ! type(profile_t), allocatable, intent(inout) :: profile_group(:) - - ! ! local variables - ! integer :: i_species - ! type(profile_t), pointer :: ptr - - ! if ( allocated(profile_group) ) then - ! write(*,*) "allocated(profile_group)" - ! ! Deallocate each element of the pointer array - ! do i_species = 1, size(profile_group) - - ! ! ptr = profile_group(i_species) - ! if ( allocated( profile_group(i_species)) ) then - ! deallocate( profile_group(i_species) ) - ! write(*,*) "deallocate( profile_group(i_species) ) ", i_species - ! ! if ( associated( ptr) ) then - ! ! write(*,*) "deallocate( profile_group(i_species) ) ", i_species - ! ! deallocate( ptr) - ! end if - ! end do - - ! ! Deallocate the entire profile_group array - ! ! if (allocated(profile_group)) then - ! ! nullify(profile_group) - ! ! end if - - ! else - ! write(*,*) "allocated(profile_group) is false" - ! return - ! end if - - ! end subroutine - - - subroutine deallocate_gas_species_profile_group(profile_group) - use musica_util, only: error_t - use musica_tuvx_profile, only: profile_t - - type(profile_t), pointer, intent(inout) :: profile_group(:) - character(len=512) :: errmsg - integer :: errcode - - ! type(error_t) :: error - integer :: i_species - logical :: deallocation_error - type(profile_t), pointer :: ptr - - write(*,*) "Called deallocate_gas_species_profile_group(profile_group)" - deallocation_error = .false. - - ! Check if the profile_group is allocated - ! if (.not. allocated(profile_group)) then - ! write(*,*) 'profile_group is not allocated.' - ! errcode = 1 - ! return - ! end if - - ! Deallocate each element of profile_group (array of pointers) - ! do i_species = 1, size(profile_group) - ! if (allocated( profile_group(i_species)) ) then - ! ! Deallocate each individual pointer - ! deallocate(profile_group(i_species), stat=errcode) - ! if (errcode /= 0) then - ! ! errmsg = 'Error deallocating profile_group(' // trim(adjustl(i_species)) // ')' - ! deallocation_error = .true. - ! exit - ! end if - ! end if - ! end do - deallocate(profile_group, stat=errcode) - !do i_species = 1, size(profile_group) - !ptr => profile_group(i_species) - ! deallocate(profile_group(i_species), stat=errcode) - !end do - ! write(*,*) "size:L profile_group ", size(profile_group) - ! ! Finally, deallocate the array of pointers itself - ! if (.not. deallocation_error) then - ! deallocate(profile_group, stat=errcode) - ! if (errcode /= 0) then - ! write(*,*) 'Error deallocating the profile_group array.' - ! else - ! write(*,*) 'profile_group successfully deallocated.' - ! end if - !end if - - end subroutine deallocate_gas_species_profile_group - - function get_height_delta(height_grid, num_interfaces, errmsg, errcode) result( height_delta ) - 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(in) :: height_grid - integer, intent(in) :: num_interfaces - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errcode + !> Sets the height delta, which is the difference in height between interfaces + subroutine set_height_delta(interfaces_height, height_delta) + real(kind_phys), intent(in) :: interfaces_height(:) + real(kind_phys), intent(inout) :: height_delta(:) ! local variables - type(error_t) :: error - real(kind_phys) :: interfaces(num_interfaces) - real(kind_phys) :: height_delta(num_interfaces-1) - integer :: i_elem - - call height_grid%get_edges(interfaces, error) - if ( has_error_occurred( error, errmsg, errcode ) ) return + integer :: i_elem do i_elem = 1, size(height_delta) - height_delta(i_elem) = interfaces(i_elem+1) - interfaces(i_elem) + height_delta(i_elem) = interfaces_height(i_elem + 1) - interfaces_height(i_elem) end do - end function get_height_delta + end subroutine set_height_delta - !> Sets TUV-x surface albedo values - subroutine set_gas_species_values(profile, gas_species, constituent, errmsg, errcode) - use ccpp_kinds, only: kind_phys + !> Sets interface, density, and above-column density values for gas species + subroutine set_gas_species_values(profile, gas_species, constituent, & + height_delta, errmsg, errcode) use musica_ccpp_util, only: has_error_occurred use musica_tuvx_profile, only: profile_t use musica_util, only: error_t - type(profile_t), intent(inout) :: profile - type(gas_species_t), intent(in) :: gas_species - real(kind_phys), intent(in) :: constituent(:) ! mol m-3 - ! real(kind_phys), intent(in) :: height_delta(:) ! km - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errcode + type(profile_t), intent(inout) :: profile + type(gas_species_t), intent(in) :: gas_species + real(kind_phys), intent(in) :: constituent(:) ! mol m-3 + real(kind_phys), intent(in) :: height_delta(:) ! km + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errcode ! local variables type(error_t) :: error @@ -304,14 +176,12 @@ subroutine set_gas_species_values(profile, gas_species, constituent, errmsg, err num_constituents = size(constituent) constituent_mol_per_cm_3(:) = constituent(:) / m_3_to_cm_3 - ! TODO(jiwon) if ( is_parameterized ) then - interfaces(1) = constituent_mol_per_cm_3(num_constituents) interfaces(2:num_constituents+1) = constituent_mol_per_cm_3(num_constituents:1:-1) interfaces(num_constituents+2) = constituent_mol_per_cm_3(1) densities(1:num_constituents+1) = & - height_delta_(1:num_constituents+1) * km_to_cm & + height_delta(1:num_constituents+1) * km_to_cm & * sqrt(interfaces(1:num_constituents+1)) & * sqrt(interfaces(2:num_constituents+2)) @@ -326,4 +196,44 @@ subroutine set_gas_species_values(profile, gas_species, constituent, errmsg, err end subroutine set_gas_species_values + !> Deallocates gas species object + subroutine deallocate_gas_species(gas_species_group) + type(gas_species_t), allocatable, intent(inout) :: gas_species_group(:) + + if (allocated( gas_species_group )) deallocate( gas_species_group ) + + end subroutine + + !> TODO(jiwon) in progress + subroutine deallocate_gas_species_profile_group(profile_group) + use musica_util, only: error_t + use musica_tuvx_profile, only: profile_t + + type(profile_group_t), allocatable, intent(inout) :: profile_group(:) + ! character(len=512) :: errmsg + integer :: errcode + + ! type(error_t) :: error + integer :: i_species + logical :: deallocation_error + type(profile_t), pointer :: profile + + write(*,*) "Called deallocate_gas_species_profile_group(profile_group)" + deallocation_error = .false. + + ! do i_species = 1, size(profile_group) + ! if (allocated( profile_group(i_species))%profile ) then + ! deallocate(profile_group(i_species)%profile, stat=errcode) + ! if (errcode /= 0) then + ! ! errmsg = 'Error deallocating profile_group(' // trim(adjustl(i_species)) // ')' + ! deallocation_error = .true. + ! exit + ! end if + ! end if + ! end do + deallocate(profile_group, stat=errcode) + + end subroutine deallocate_gas_species_profile_group + + end module musica_ccpp_tuvx_gas_species_profile \ No newline at end of file diff --git a/test/musica/test_musica_api.F90 b/test/musica/test_musica_api.F90 index 95d73c80..2976c42a 100644 --- a/test/musica/test_musica_api.F90 +++ b/test/musica/test_musica_api.F90 @@ -7,8 +7,13 @@ program run_test_musica_ccpp #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 + write(*,*) "[MUSICA Test] Running the Chapman test" call test_chapman() + write(*,*) "[MUSICA Test] Ends the Chapman test" + + write(*,*) "[MUSICA Test] Running the Terminator test" call test_terminator() + write(*,*) "[MUSICA Test] Ends the Terminator test" contains From b24f5ee550fefba4a545a9a10a9c3d02e906cb1c Mon Sep 17 00:00:00 2001 From: Jiwon Gim Date: Thu, 5 Dec 2024 17:56:14 -0700 Subject: [PATCH 10/14] test passed --- schemes/musica/tuvx/musica_ccpp_tuvx.F90 | 65 +++++++------- .../musica_ccpp_tuvx_gas_species_profile.F90 | 89 +++++++------------ .../tuvx/musica_ccpp_tuvx_height_grid.F90 | 15 +++- test/musica/tuvx/test_tuvx_height_grid.F90 | 3 +- 4 files changed, 76 insertions(+), 96 deletions(-) diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 index 8077c33f..73b60fff 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 @@ -23,7 +23,6 @@ module musica_ccpp_tuvx type(radiator_t), pointer :: cloud_optics => null() type(index_mappings_t), pointer :: photolysis_rate_constants_mapping => null( ) - real(kind_phys), allocatable :: height_delta(:) ! tuvx_run needs this type(gas_species_t), allocatable :: gas_species_group(:) type(profile_group_t), allocatable :: profile_gas_species_group(:) @@ -59,6 +58,9 @@ subroutine cleanup_tuvx_resources() use musica_ccpp_tuvx_gas_species_profile, & only: deallocate_gas_species_profile_group + ! local variable + integer :: i_species + if (associated( height_grid )) then deallocate( height_grid ) height_grid => null() @@ -94,11 +96,8 @@ subroutine cleanup_tuvx_resources() photolysis_rate_constants_mapping => null() end if - if (allocated( profile_gas_species_group )) then - deallocate( profile_gas_species_group ) - end if - call deallocate_gas_species_profile_group( profile_gas_species_group ) + end subroutine cleanup_tuvx_resources !> Registers constituent properties with the CCPP needed by TUV-x @@ -156,7 +155,7 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & use musica_ccpp_tuvx_cloud_optics, & only: create_cloud_optics_radiator, cloud_optics_label use musica_ccpp_tuvx_gas_species_profile, & - only: configure_gas_species, create_gas_species_profile_group, set_height_delta + only: configure_gas_species, create_gas_species_profile_group integer, intent(in) :: vertical_layer_dimension ! (count) integer, intent(in) :: vertical_interface_dimension ! (count) @@ -173,7 +172,6 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & type(radiator_map_t), pointer :: radiators type(configuration_t) :: config type(mappings_t), pointer :: photolysis_rate_constants_ordering - real(kind_phys) :: interfaces_height(vertical_interface_dimension) integer :: number_of_gas_species integer :: i_species @@ -273,14 +271,6 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & return end if - ! Gets interfaces of height to set height delta values - ! allocate( interfaces_height( vertical_interface_dimension) ) - call height_grid%get_edges(interfaces_height, error) - - allocate( height_delta( vertical_interface_dimension - 1) ) - call set_height_delta(interfaces_height, height_delta) - ! deallocate( interfaces_height ) - number_of_gas_species = size(gas_species_group) allocate( profile_gas_species_group( number_of_gas_species ) ) @@ -404,7 +394,6 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & end if do i_species = 1, number_of_gas_species - write(*,*) "gas_species_group(i_species)%label",gas_species_group(i_species)%label profile_gas_species_group(i_species)%profile => & profiles%get( gas_species_group(i_species)%label, gas_species_group(i_species)%unit, error ) if (has_error_occurred( error, errmsg, errcode )) then @@ -484,14 +473,15 @@ subroutine tuvx_run(temperature, dry_air_density, & 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_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 - use musica_ccpp_tuvx_gas_species_profile, only: set_gas_species_values + 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_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 + use musica_ccpp_tuvx_gas_species_profile, only: set_gas_species_values, MOLAR_MASS_DRY_AIR + real(kind_phys), intent(in) :: temperature(:,:) ! K (column, layer) real(kind_phys), intent(in) :: dry_air_density(:,:) ! kg m-3 (column, layer) real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_midpoint(:,:) ! m (column, layer) @@ -513,6 +503,7 @@ subroutine tuvx_run(temperature, dry_air_density, & ! local variables 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 + real(kind_phys), dimension(size(height_interfaces)) :: height_deltas ! km real(kind_phys), dimension(size(rate_parameters, dim=2)+2, & number_of_photolysis_rate_constants) :: photolysis_rate_constants, & ! s-1 heating_rates ! K s-1 (TODO: check units) @@ -542,7 +533,7 @@ subroutine tuvx_run(temperature, dry_air_density, & reciprocal_of_gravitational_acceleration, & height_midpoints, height_interfaces ) call set_height_grid_values( height_grid, height_midpoints, height_interfaces, & - errmsg, errcode ) + height_deltas, errmsg, errcode ) if (errcode /= 0) return call set_temperature_values( temperature_profile, temperature(i_col,:), & @@ -557,11 +548,20 @@ subroutine tuvx_run(temperature, dry_air_density, & if (errcode /= 0) return do i_gas_species = 1, size(gas_species_group) - call set_gas_species_values(profile_gas_species_group(i_gas_species)%profile, & - gas_species_group(i_gas_species), & - constituents(i_col,:,gas_species_group(i_gas_species)%index_constituent_props), & ! mol m-3 - height_delta, errmsg, errcode) - if (errcode /= 0) return + ! TODO(jiwon) + if (gas_species_group(i_gas_species)%label == "air") then + call set_gas_species_values(profile_gas_species_group(i_gas_species)%profile, & + gas_species_group(i_gas_species), & + dry_air_density(i_col,:) * MOLAR_MASS_DRY_AIR, & ! mol m-3 + height_deltas, errmsg, errcode) + if (errcode /= 0) return + else + call set_gas_species_values(profile_gas_species_group(i_gas_species)%profile, & + gas_species_group(i_gas_species), & + constituents(i_col,:,gas_species_group(i_gas_species)%index_constituent_props), & ! mol m-3 + height_deltas, errmsg, errcode) + if (errcode /= 0) return + end if end do ! temporary values until these are available from the host model @@ -601,11 +601,6 @@ subroutine tuvx_final(errmsg, errcode) call deallocate_gas_species( gas_species_group ) call cleanup_tuvx_resources() - ! TODO(jiwon) is it okay to deallocate this in the final call? - if (allocated( height_delta )) then - deallocate( height_delta ) - end if - if (associated( tuvx )) then deallocate( tuvx ) tuvx => null() diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 index 273041f8..4e0020f6 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 @@ -6,13 +6,15 @@ module musica_ccpp_tuvx_gas_species_profile private public :: create_gas_species, configure_gas_species, create_gas_species_profile_group, & - set_height_delta, set_gas_species_values, deallocate_gas_species, & - deallocate_gas_species_profile_group + set_gas_species_values, deallocate_gas_species, deallocate_gas_species_profile_group !> Conversion factor from km to cm real(kind_phys), parameter, public :: km_to_cm = 1.0e5 !> Conversion factor from m3 to cm3 real(kind_phys), parameter, public :: m_3_to_cm_3 = 1.0e6 + !> Molar mass value of dry air is obtained from 'CAM-SIMA/src/utils/std_atm_profile.F90' + ! TODO(jiwon) + real(kind_phys), parameter, public :: MOLAR_MASS_DRY_AIR = 0.0289644_kind_phys ! kg mol-1 !> Definition of the gas species type type, public :: gas_species_t @@ -66,12 +68,11 @@ subroutine configure_gas_species(constituent_props, gas_species_group, errmsg, e real(kind_phys), parameter :: SCALE_HEIGHT_DRY_AIR = 8.01_kind_phys ! km real(kind_phys), parameter :: SCALE_HEIGHT_O2 = 7.0_kind_phys ! km real(kind_phys), parameter :: SCALE_HEIGHT_O3 = 7.0_kind_phys ! km - ! molar mass of dry air is from 'CAM-SIMA/src/utils/std_atm_profile.F90' - real(kind_phys), parameter :: MOLAR_MASS_DRY_AIR = 0.0289644_kind_phys ! kg mol-1 real(kind_phys) :: molar_mass_O2 ! kg mol-1 real(kind_phys) :: molar_mass_O3 ! kg mol-1 + integer, parameter :: INDEX_NOT_KNOWN = -9999 + integer :: index_air, index_O2, index_O3 - integer :: index_air, index_O2, index_O3 allocate(gas_species_group(num_gas_species)) ! need to deallocate ! TODO(jiwon) - I commented out this block of code that searches for the molar mass @@ -88,6 +89,7 @@ subroutine configure_gas_species(constituent_props, gas_species_group, errmsg, e ! call constituent_props(index_air)%molar_mass(molar_mass_air, errcode, errmsg) ! if (errcode /= 0) return + index_air = INDEX_NOT_KNOWN gas_species_group(1) = create_gas_species("air", "molecule cm-3", & MOLAR_MASS_DRY_AIR, SCALE_HEIGHT_DRY_AIR, index_air) ! O2 @@ -138,23 +140,9 @@ subroutine create_gas_species_profile_group(height_grid, gas_species_group, & end subroutine create_gas_species_profile_group - !> Sets the height delta, which is the difference in height between interfaces - subroutine set_height_delta(interfaces_height, height_delta) - real(kind_phys), intent(in) :: interfaces_height(:) - real(kind_phys), intent(inout) :: height_delta(:) - - ! local variables - integer :: i_elem - - do i_elem = 1, size(height_delta) - height_delta(i_elem) = interfaces_height(i_elem + 1) - interfaces_height(i_elem) - end do - - end subroutine set_height_delta - !> Sets interface, density, and above-column density values for gas species subroutine set_gas_species_values(profile, gas_species, constituent, & - height_delta, errmsg, errcode) + height_deltas, errmsg, errcode) use musica_ccpp_util, only: has_error_occurred use musica_tuvx_profile, only: profile_t use musica_util, only: error_t @@ -162,28 +150,28 @@ subroutine set_gas_species_values(profile, gas_species, constituent, & type(profile_t), intent(inout) :: profile type(gas_species_t), intent(in) :: gas_species real(kind_phys), intent(in) :: constituent(:) ! mol m-3 - real(kind_phys), intent(in) :: height_delta(:) ! km + real(kind_phys), intent(in) :: height_deltas(:) ! km, change in height in each vertical layer character(len=*), intent(out) :: errmsg integer, intent(out) :: errcode ! local variables - type(error_t) :: error - integer :: num_constituents - real(kind_phys) :: constituent_mol_per_cm_3(size(constituent)) ! mol cm-3 - real(kind_phys) :: interfaces(size(constituent) + 2) - real(kind_phys) :: densities(size(constituent) + 1) + type(error_t) :: error + integer :: num_constituents + real(kind_phys) :: constituent_mol_per_cm_3(size(constituent)) ! mol cm-3 + real(kind_phys) :: interfaces(size(constituent) + 2) + real(kind_phys) :: densities(size(constituent) + 1) num_constituents = size(constituent) constituent_mol_per_cm_3(:) = constituent(:) / m_3_to_cm_3 + interfaces(1) = constituent_mol_per_cm_3(num_constituents) interfaces(2:num_constituents+1) = constituent_mol_per_cm_3(num_constituents:1:-1) interfaces(num_constituents+2) = constituent_mol_per_cm_3(1) - densities(1:num_constituents+1) = & - height_delta(1:num_constituents+1) * km_to_cm & - * sqrt(interfaces(1:num_constituents+1)) & - * sqrt(interfaces(2:num_constituents+2)) + densities(:) = height_deltas(:) * km_to_cm & + * sqrt(interfaces(1:num_constituents+1)) & + * sqrt(interfaces(2:num_constituents+2)) call profile%set_edge_values( interfaces, error ) if ( has_error_occurred( error, errmsg, errcode ) ) return @@ -204,36 +192,23 @@ subroutine deallocate_gas_species(gas_species_group) end subroutine - !> TODO(jiwon) in progress + !> Deallocates profiles of each gas species subroutine deallocate_gas_species_profile_group(profile_group) - use musica_util, only: error_t - use musica_tuvx_profile, only: profile_t - type(profile_group_t), allocatable, intent(inout) :: profile_group(:) - ! character(len=512) :: errmsg - integer :: errcode - - ! type(error_t) :: error + + ! local variables integer :: i_species - logical :: deallocation_error - type(profile_t), pointer :: profile - - write(*,*) "Called deallocate_gas_species_profile_group(profile_group)" - deallocation_error = .false. - - ! do i_species = 1, size(profile_group) - ! if (allocated( profile_group(i_species))%profile ) then - ! deallocate(profile_group(i_species)%profile, stat=errcode) - ! if (errcode /= 0) then - ! ! errmsg = 'Error deallocating profile_group(' // trim(adjustl(i_species)) // ')' - ! deallocation_error = .true. - ! exit - ! end if - ! end if - ! end do - deallocate(profile_group, stat=errcode) - + + if (allocated( profile_group )) then + do i_species = 1, size(profile_group) + if (associated( profile_group(i_species)%profile) ) then + deallocate( profile_group(i_species)%profile) + end if + end do + + deallocate( profile_group ) + end if + end subroutine deallocate_gas_species_profile_group - end module musica_ccpp_tuvx_gas_species_profile \ No newline at end of file diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_height_grid.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_height_grid.F90 index cff0d2a2..ec341196 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_height_grid.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_height_grid.F90 @@ -84,7 +84,7 @@ end function create_height_grid !> Sets TUV-x height grid values from the host-model height grid subroutine set_height_grid_values(height_grid, host_midpoints, & - host_interfaces, errmsg, errcode) + host_interfaces, height_deltas, errmsg, errcode) use ccpp_kinds, only: kind_phys use musica_ccpp_util, only: has_error_occurred use musica_tuvx_grid, only: grid_t @@ -93,6 +93,7 @@ subroutine set_height_grid_values(height_grid, host_midpoints, & type(grid_t), intent(inout) :: height_grid real(kind_phys), intent(in) :: host_midpoints(:) ! km real(kind_phys), intent(in) :: host_interfaces(:) ! km + real(kind_phys), intent(inout) :: height_deltas(:) ! km character(len=*), intent(out) :: errmsg integer, intent(out) :: errcode @@ -102,6 +103,9 @@ subroutine set_height_grid_values(height_grid, host_midpoints, & real(kind_phys) :: interfaces(size(host_interfaces)+1) integer :: n_host_midpoints, n_host_interfaces + n_host_midpoints = size(host_midpoints) + n_host_interfaces = size(host_interfaces) + if ( size(midpoints) /= height_grid%number_of_sections( error ) ) then errmsg = "[MUSICA Error] Invalid size of TUV-x mid-point heights." errcode = 1 @@ -114,8 +118,11 @@ subroutine set_height_grid_values(height_grid, host_midpoints, & return end if - n_host_midpoints = size(host_midpoints) - n_host_interfaces = size(host_interfaces) + if ( size(height_deltas) /= n_host_interfaces ) then + errmsg = "[MUSICA Error] Invalid size of TUV-x height deltas." + errcode = 1 + return + end if interfaces(1) = host_interfaces(n_host_interfaces) interfaces(2:n_host_interfaces) = host_midpoints(n_host_midpoints:1:-1) @@ -127,6 +134,8 @@ subroutine set_height_grid_values(height_grid, host_midpoints, & midpoints(n_host_midpoints+1) = 0.5_kind_phys * & ( interfaces(n_host_interfaces) + interfaces(n_host_interfaces+1) ) + height_deltas(:) = interfaces(2:size(interfaces)) - interfaces(1:size(interfaces)-1) + call height_grid%set_edges( interfaces, error ) if ( has_error_occurred( error, errmsg, errcode ) ) return diff --git a/test/musica/tuvx/test_tuvx_height_grid.F90 b/test/musica/tuvx/test_tuvx_height_grid.F90 index 6da26960..ef774ce9 100644 --- a/test/musica/tuvx/test_tuvx_height_grid.F90 +++ b/test/musica/tuvx/test_tuvx_height_grid.F90 @@ -27,6 +27,7 @@ subroutine test_create_height_grid() real(kind_phys) :: expected_interfaces(NUM_HOST_INTERFACES+1) = [50.1_kind_phys, 100.6_kind_phys, 200.8_kind_phys, 250.3_kind_phys] real(kind_phys) :: midpoints(NUM_HOST_MIDPOINTS+1) real(kind_phys) :: interfaces(NUM_HOST_INTERFACES+1) + real(kind_phys) :: height_deltas(NUM_HOST_INTERFACES) type(grid_t), pointer :: height_grid => null() type(error_t) :: error character(len=512) :: errmsg @@ -47,7 +48,7 @@ subroutine test_create_height_grid() ASSERT(associated(height_grid)) call set_height_grid_values(height_grid, host_midpoints, host_interfaces, & - errmsg, errcode) + height_deltas, errmsg, errcode) ASSERT(errcode == 0) ASSERT(height_grid%number_of_sections(error) == NUM_HOST_MIDPOINTS + 1) From ce5190fbff1bb081b1f1ab91efaf89e1755b96c0 Mon Sep 17 00:00:00 2001 From: Jiwon Gim Date: Thu, 5 Dec 2024 17:58:09 -0700 Subject: [PATCH 11/14] update cam sima data commit --- test/docker/Dockerfile.musica | 2 +- test/docker/Dockerfile.musica.no_install | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/test/docker/Dockerfile.musica b/test/docker/Dockerfile.musica index e98255f7..15e46b98 100644 --- a/test/docker/Dockerfile.musica +++ b/test/docker/Dockerfile.musica @@ -6,7 +6,7 @@ FROM ubuntu:22.04 ARG MUSICA_GIT_TAG=326b5119768d5be9654baf96ae3bd6a1b757fdc8 -ARG CAM_SIMA_CHEMISTRY_DATA_TAG=abc7cacbec3d33d5c0ed5bb79a157e93b42c45c0 +ARG CAM_SIMA_CHEMISTRY_DATA_TAG=9e84370d917fa4357cb3e7d9a0f3b57869fb342c ARG BUILD_TYPE=Debug RUN apt update \ diff --git a/test/docker/Dockerfile.musica.no_install b/test/docker/Dockerfile.musica.no_install index cb4edf56..940f252a 100644 --- a/test/docker/Dockerfile.musica.no_install +++ b/test/docker/Dockerfile.musica.no_install @@ -9,7 +9,7 @@ FROM ubuntu:22.04 ARG MUSICA_GIT_TAG=326b5119768d5be9654baf96ae3bd6a1b757fdc8 -ARG CAM_SIMA_CHEMISTRY_DATA_TAG=abc7cacbec3d33d5c0ed5bb79a157e93b42c45c0 +ARG CAM_SIMA_CHEMISTRY_DATA_TAG=9e84370d917fa4357cb3e7d9a0f3b57869fb342c ARG BUILD_TYPE=Debug RUN apt update \ From 8afeebd40a070cf9bcc5fa844357ca28b20e3465 Mon Sep 17 00:00:00 2001 From: Jiwon Gim Date: Fri, 6 Dec 2024 16:11:20 -0700 Subject: [PATCH 12/14] add a unit test for gas species profile --- schemes/musica/tuvx/musica_ccpp_tuvx.F90 | 19 ++- .../musica_ccpp_tuvx_gas_species_profile.F90 | 13 +- test/musica/tuvx/CMakeLists.txt | 34 +++++ .../tuvx/test_tuvx_gas_species_profile.F90 | 125 ++++++++++++++++++ test/musica/tuvx/test_tuvx_height_grid.F90 | 8 ++ to_be_ccppized/ccpp_tuvx_utils.F90 | 2 +- 6 files changed, 185 insertions(+), 16 deletions(-) create mode 100644 test/musica/tuvx/test_tuvx_gas_species_profile.F90 diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 index a5bb7d75..c47c1766 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 @@ -509,6 +509,7 @@ subroutine tuvx_run(temperature, dry_air_density, & 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 real(kind_phys), dimension(size(height_interfaces)) :: height_deltas ! km + real(kind_phys), dimension(size(constituents, dim = 2)) :: gas_species_constituents real(kind_phys), dimension(size(rate_parameters, dim=2)+2, & number_of_photolysis_rate_constants) :: photolysis_rate_constants, & ! s-1 heating_rates ! K s-1 (TODO: check units) @@ -553,20 +554,16 @@ subroutine tuvx_run(temperature, dry_air_density, & if (errcode /= 0) return do i_gas_species = 1, size(gas_species_group) - ! TODO(jiwon) if (gas_species_group(i_gas_species)%label == "air") then - call set_gas_species_values(profile_gas_species_group(i_gas_species)%profile, & - gas_species_group(i_gas_species), & - dry_air_density(i_col,:) * MOLAR_MASS_DRY_AIR, & ! mol m-3 - height_deltas, errmsg, errcode) - if (errcode /= 0) return + gas_species_constituents = dry_air_density(i_col,:) * MOLAR_MASS_DRY_AIR else - call set_gas_species_values(profile_gas_species_group(i_gas_species)%profile, & - gas_species_group(i_gas_species), & - constituents(i_col,:,gas_species_group(i_gas_species)%index_constituent_props), & ! mol m-3 - height_deltas, errmsg, errcode) - if (errcode /= 0) return + gas_species_constituents = constituents(i_col,:,gas_species_group(i_gas_species)%index_constituent_props) end if + call set_gas_species_values(profile_gas_species_group(i_gas_species)%profile, & + gas_species_group(i_gas_species), & + gas_species_constituents, & ! mol m-3 + height_deltas, errmsg, errcode) + if (errcode /= 0) return end do ! temporary values until these are available from the host model diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 index 4e0020f6..5be6bc81 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 @@ -169,10 +169,15 @@ subroutine set_gas_species_values(profile, gas_species, constituent, & interfaces(2:num_constituents+1) = constituent_mol_per_cm_3(num_constituents:1:-1) interfaces(num_constituents+2) = constituent_mol_per_cm_3(1) - densities(:) = height_deltas(:) * km_to_cm & - * sqrt(interfaces(1:num_constituents+1)) & - * sqrt(interfaces(2:num_constituents+2)) - + if (gas_species%label == "O3") then + densities(:) = height_deltas(:) * km_to_cm & + * ( interfaces(1:num_constituents+1) & + + interfaces(2:num_constituents+2) ) * 0.5_kind_phys + else + densities(:) = height_deltas(:) * km_to_cm & + * sqrt(interfaces(1:num_constituents+1)) & + * sqrt(interfaces(2:num_constituents+2)) + end if call profile%set_edge_values( interfaces, error ) if ( has_error_occurred( error, errmsg, errcode ) ) return diff --git a/test/musica/tuvx/CMakeLists.txt b/test/musica/tuvx/CMakeLists.txt index 21301b0b..a2980e63 100644 --- a/test/musica/tuvx/CMakeLists.txt +++ b/test/musica/tuvx/CMakeLists.txt @@ -171,3 +171,37 @@ add_test( ) add_memory_check_test(test_tuvx_cloud_optics $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) + +# Gas species profile + +add_executable(test_tuvx_gas_species_profile test_tuvx_gas_species_profile.F90) + +target_sources(test_tuvx_gas_species_profile + PUBLIC + ${MUSICA_CCPP_SOURCES} + ${TO_BE_CCPPIZED_SRC_PATH}/ccpp_const_utils.F90 + ${TO_BE_CCPPIZED_SRC_PATH}/ccpp_tuvx_utils.F90 + ${CCPP_SRC_PATH}/ccpp_constituent_prop_mod.F90 + ${CCPP_SRC_PATH}/ccpp_hash_table.F90 + ${CCPP_SRC_PATH}/ccpp_hashable.F90 + ${CCPP_TEST_SRC_PATH}/ccpp_kinds.F90 + ${CMAKE_CURRENT_SOURCE_DIR}/../musica_ccpp_namelist.F90 +) + +target_link_libraries(test_tuvx_gas_species_profile + PRIVATE + musica::musica-fortran +) + +set_target_properties(test_tuvx_gas_species_profile + PROPERTIES + LINKER_LANGUAGE Fortran +) + +add_test( + NAME test_tuvx_gas_species_profile + COMMAND $ + WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} +) + +add_memory_check_test(test_tuvx_gas_species_profile $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) \ No newline at end of file diff --git a/test/musica/tuvx/test_tuvx_gas_species_profile.F90 b/test/musica/tuvx/test_tuvx_gas_species_profile.F90 new file mode 100644 index 00000000..14ce341e --- /dev/null +++ b/test/musica/tuvx/test_tuvx_gas_species_profile.F90 @@ -0,0 +1,125 @@ +program test_tuvx_gas_species_profile + + use musica_ccpp_tuvx_gas_species_profile + + 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_gas_species_profile() + +contains + + subroutine test_create_gas_species_profile() + use musica_micm, only: Rosenbrock, RosenbrockStandardOrder + use musica_tuvx_grid, only: grid_t + use ccpp_kinds, only: kind_phys + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + use musica_ccpp, only: musica_ccpp_register + use musica_ccpp_tuvx_height_grid, only: create_height_grid + use musica_ccpp_namelist, only: filename_of_micm_configuration, & + filename_of_tuvx_configuration, & + filename_of_tuvx_micm_mapping_configuration + + integer, parameter :: NUM_COLUMNS = 2 + integer, parameter :: NUM_LAYERS = 2 + integer :: NUM_GRID_CELLS = NUM_COLUMNS * NUM_LAYERS + integer :: solver_type = Rosenbrock + ! Molar mass for O2, O3 is read from 'cam-sima-chemistry-data/mechanisms/chapman/micm/species.json' + real(kind_phys), parameter :: MOLAR_MASS_O2 = 0.0319988_kind_phys + real(kind_phys), parameter :: MOLAR_MASS_O3 = 0.0479982_kind_phys + integer, parameter :: NUMBER_GAS_SPECIES = 3 + integer, parameter :: INDEX_NOT_KNOWN = -9999 + 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 + type(gas_species_t), allocatable :: gas_species_group(:) + type(profile_group_t), allocatable :: profile_gas_species_group(:) + type(grid_t), pointer :: height_grid => null() + real(kind_phys) :: dry_air_density(NUM_COLUMNS,NUM_LAYERS) ! kg m-3 + real(kind_phys) :: gas_species_constituents(NUM_COLUMNS,NUM_LAYERS) + real(kind_phys) :: height_deltas(NUM_LAYERS) ! km + integer :: errcode + character(len=512) :: errmsg + character(len=50) :: label + character(len=50) :: unit + real(kind_phys) :: molar_mass ! kg mol-1 + real(kind_phys) :: scale_height ! km + integer :: index + logical :: tmp_bool + integer :: i, i_gas_species, i_col + + filename_of_micm_configuration = 'musica_configurations/chapman/micm/config.json' + filename_of_tuvx_configuration = 'musica_configurations/chapman/tuvx/config.json' + filename_of_tuvx_micm_mapping_configuration = 'musica_configurations/chapman/tuvx_micm_mapping.json' + + dry_air_density(:,1) = (/ 3.5_kind_phys, 4.5_kind_phys /) + dry_air_density(:,2) = (/ 5.5_kind_phys, 6.5_kind_phys /) + gas_species_constituents(:,1) = (/ 0.35_kind_phys, 0.45_kind_phys /) + gas_species_constituents(:,2) = (/ 0.55_kind_phys, 0.65_kind_phys /) + height_deltas(:) = (/ 0.5_kind_phys, 1.5_kind_phys /) + + call musica_ccpp_register(solver_type, NUM_GRID_CELLS, constituent_props, errmsg, errcode) + if (errcode /= 0) then + write(*,*) trim(errmsg) + stop 3 + endif + + allocate(constituent_props_ptr(size(constituent_props))) + do i = 1, size(constituent_props) + const_prop => constituent_props(i) + call constituent_props_ptr(i)%set( const_prop, errcode, errmsg ) + end do + + height_grid => create_height_grid( NUM_LAYERS, NUM_LAYERS + 1 , errmsg, errcode ) + ASSERT(errcode == 0) + ASSERT(associated(height_grid)) + + call configure_gas_species( constituent_props_ptr, gas_species_group, errmsg, errcode ) + ASSERT(errcode == 0) + ASSERT(size(gas_species_group) == NUMBER_GAS_SPECIES) + + allocate( profile_gas_species_group( size(gas_species_group) ) ) + + call create_gas_species_profile_group( height_grid, gas_species_group, & + profile_gas_species_group, errmsg, errcode ) + ASSERT(errcode == 0) + ASSERT(size(profile_gas_species_group) == NUMBER_GAS_SPECIES) + + do i_gas_species = 1, size(gas_species_group) + label = gas_species_group(i_gas_species)%label + unit = gas_species_group(i_gas_species)%unit + molar_mass = gas_species_group(i_gas_species)%molar_mass + scale_height = gas_species_group(i_gas_species)%scale_height + index = gas_species_group(i_gas_species)%index_constituent_props + + tmp_bool = (trim(label) == "air" .and. trim(unit) == "molecule cm-3" .and. molar_mass == MOLAR_MASS_DRY_AIR .and. & + scale_height == 8.01_kind_phys .and. index == INDEX_NOT_KNOWN) .or. & + (trim(label) == "O2" .and. trim(unit) == "molecule cm-3" .and. molar_mass == MOLAR_MASS_O2 .and. & + scale_height == 7.0_kind_phys .and. index /= INDEX_NOT_KNOWN .and. index /= 0) .or. & + (trim(label) == "O3" .and. trim(unit) == "molecule cm-3" .and. molar_mass == MOLAR_MASS_O3 .and. & + scale_height == 7.0_kind_phys .and. index /= INDEX_NOT_KNOWN .and. index /= 0) + ASSERT(tmp_bool) + + ASSERT(associated(profile_gas_species_group(i_gas_species)%profile)) + end do + + do i_col = 1, NUM_COLUMNS + do i_gas_species = 1, size(gas_species_group) + call set_gas_species_values(profile_gas_species_group(i_gas_species)%profile, & + gas_species_group(i_gas_species), & + gas_species_constituents(i_col, :), & ! mol m-3 + height_deltas, errmsg, errcode) + if (errcode /= 0) return + end do + end do + + call deallocate_gas_species(gas_species_group) + call deallocate_gas_species_profile_group(profile_gas_species_group) + + end subroutine test_create_gas_species_profile + +end program test_tuvx_gas_species_profile \ No newline at end of file diff --git a/test/musica/tuvx/test_tuvx_height_grid.F90 b/test/musica/tuvx/test_tuvx_height_grid.F90 index ef774ce9..c18ce338 100644 --- a/test/musica/tuvx/test_tuvx_height_grid.F90 +++ b/test/musica/tuvx/test_tuvx_height_grid.F90 @@ -25,6 +25,7 @@ subroutine test_create_height_grid() real(kind_phys) :: host_interfaces(NUM_HOST_INTERFACES) = [250.3_kind_phys, 150.2_kind_phys, 50.1_kind_phys] real(kind_phys) :: expected_midpoints(NUM_HOST_MIDPOINTS+1) = [(100.6 + 50.1) * 0.5, 150.2, (250.3 + 200.8) * 0.5] real(kind_phys) :: expected_interfaces(NUM_HOST_INTERFACES+1) = [50.1_kind_phys, 100.6_kind_phys, 200.8_kind_phys, 250.3_kind_phys] + real(kind_phys) :: expected_height_deltas(NUM_HOST_INTERFACES) real(kind_phys) :: midpoints(NUM_HOST_MIDPOINTS+1) real(kind_phys) :: interfaces(NUM_HOST_INTERFACES+1) real(kind_phys) :: height_deltas(NUM_HOST_INTERFACES) @@ -34,6 +35,9 @@ subroutine test_create_height_grid() integer :: errcode integer :: i + expected_height_deltas(:) = expected_interfaces(2:size(expected_interfaces)) & + - expected_interfaces(1:size(expected_interfaces)-1) + height_grid => create_height_grid(-1, 0, errmsg, errcode) ASSERT(errcode == 1) ASSERT(.not. associated(height_grid)) @@ -54,6 +58,10 @@ subroutine test_create_height_grid() ASSERT(height_grid%number_of_sections(error) == NUM_HOST_MIDPOINTS + 1) ASSERT(error%is_success()) + do i = 1, size(height_deltas) + ASSERT_NEAR(height_deltas(i), expected_height_deltas(i), ABS_ERROR) + end do + call height_grid%get_midpoints(midpoints, error) ASSERT(error%is_success()) do i = 1, size(midpoints) diff --git a/to_be_ccppized/ccpp_tuvx_utils.F90 b/to_be_ccppized/ccpp_tuvx_utils.F90 index 837801bc..e7329d56 100644 --- a/to_be_ccppized/ccpp_tuvx_utils.F90 +++ b/to_be_ccppized/ccpp_tuvx_utils.F90 @@ -10,7 +10,7 @@ module ccpp_tuvx_utils contains !> Regrids normalized flux data to match a specified wavelength grid - !! This function is copied from CAM/src/chemistry/utils/mo_util.F90 + ! This function is copied from CAM/src/chemistry/utils/mo_util.F90 subroutine rebin( source_dimension, target_dimension, source_coordinates, & target_coordinates, source, target ) use ccpp_kinds, only: kind_phys From 274f446451a4346a30c3365a5235302f30a5ebeb Mon Sep 17 00:00:00 2001 From: Jiwon Gim Date: Fri, 6 Dec 2024 16:46:09 -0700 Subject: [PATCH 13/14] minor fix --- schemes/musica/tuvx/musica_ccpp_tuvx.F90 | 32 +++++++++---------- .../musica_ccpp_tuvx_gas_species_profile.F90 | 6 ++-- test/musica/tuvx/CMakeLists.txt | 1 - 3 files changed, 19 insertions(+), 20 deletions(-) diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 index c47c1766..7045cf64 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 @@ -14,28 +14,28 @@ module musica_ccpp_tuvx public :: tuvx_register, tuvx_init, tuvx_run, tuvx_final - type(tuvx_t), pointer :: tuvx => null() - type(grid_t), pointer :: height_grid => null() - type(grid_t), pointer :: wavelength_grid => null() - type(profile_t), pointer :: temperature_profile => null() - type(profile_t), pointer :: surface_albedo_profile => null() - type(profile_t), pointer :: extraterrestrial_flux_profile => null() - type(radiator_t), pointer :: cloud_optics => null() - type(index_mappings_t), pointer :: photolysis_rate_constants_mapping => null( ) + type(tuvx_t), pointer :: tuvx => null() + type(grid_t), pointer :: height_grid => null() + type(grid_t), pointer :: wavelength_grid => null() + type(profile_t), pointer :: temperature_profile => null() + type(profile_t), pointer :: surface_albedo_profile => null() + type(profile_t), pointer :: extraterrestrial_flux_profile => null() + type(radiator_t), pointer :: cloud_optics => null() + type(index_mappings_t), pointer :: photolysis_rate_constants_mapping => null( ) type(gas_species_t), allocatable :: gas_species_group(:) type(profile_group_t), allocatable :: profile_gas_species_group(:) - integer, parameter :: DEFAULT_NUM_PHOTOLYSIS_RATE_CONSTANTS = 0 - integer :: number_of_photolysis_rate_constants = DEFAULT_NUM_PHOTOLYSIS_RATE_CONSTANTS - integer, parameter :: DEFAULT_INDEX_NOT_FOUND = -1 - character(len=*), parameter :: CLOUD_LIQUID_WATER_CONTENT_LABEL = & + integer, parameter :: DEFAULT_NUM_PHOTOLYSIS_RATE_CONSTANTS = 0 + integer :: number_of_photolysis_rate_constants = DEFAULT_NUM_PHOTOLYSIS_RATE_CONSTANTS + integer, parameter :: DEFAULT_INDEX_NOT_FOUND = -1 + character(len=*), parameter :: CLOUD_LIQUID_WATER_CONTENT_LABEL = & 'cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water' - character(len=*), parameter :: CLOUD_LIQUID_WATER_CONTENT_LONG_NAME = & + character(len=*), parameter :: CLOUD_LIQUID_WATER_CONTENT_LONG_NAME = & 'Cloud water mass mixing ratio with respect to moist air plus all airborne condensates' - character(len=*), parameter :: CLOUD_LIQUID_WATER_CONTENT_UNITS = 'kg kg-1' - real(kind_phys), parameter :: CLOUD_LIQUID_WATER_CONTENT_MOLAR_MASS = 0.018_kind_phys ! kg mol-1 - integer :: index_cloud_liquid_water_content = DEFAULT_INDEX_NOT_FOUND + character(len=*), parameter :: CLOUD_LIQUID_WATER_CONTENT_UNITS = 'kg kg-1' + real(kind_phys), parameter :: CLOUD_LIQUID_WATER_CONTENT_MOLAR_MASS = 0.018_kind_phys ! kg mol-1 + integer :: index_cloud_liquid_water_content = DEFAULT_INDEX_NOT_FOUND contains diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 index 5be6bc81..cc58e500 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 @@ -77,9 +77,9 @@ subroutine configure_gas_species(constituent_props, gas_species_group, errmsg, e ! TODO(jiwon) - I commented out this block of code that searches for the molar mass ! of air and instead hard-coded it using the value imported from CAM-SIMA. - ! Should we register air in the register phase like cloud liquid water content, - ! or should we include it in the configuration file, or is it good as it is? - + ! Should we register air during the registration phase, similar to cloud liquid water content? + ! Or should we include it in the configuration file, or is it fine as it is? + ! Air ! Find the index of the species and create a 'gas_species_t' for it ! call ccpp_const_get_idx(constituent_props, "air", index_air, errmsg, errcode) diff --git a/test/musica/tuvx/CMakeLists.txt b/test/musica/tuvx/CMakeLists.txt index a2980e63..14f26b46 100644 --- a/test/musica/tuvx/CMakeLists.txt +++ b/test/musica/tuvx/CMakeLists.txt @@ -173,7 +173,6 @@ add_test( add_memory_check_test(test_tuvx_cloud_optics $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) # Gas species profile - add_executable(test_tuvx_gas_species_profile test_tuvx_gas_species_profile.F90) target_sources(test_tuvx_gas_species_profile From 3f55081e2512108d0a4a2b9fa2569452ecf4d28c Mon Sep 17 00:00:00 2001 From: Jiwon Gim Date: Fri, 13 Dec 2024 17:28:31 -0700 Subject: [PATCH 14/14] git stash --- schemes/musica/micm/musica_ccpp_micm.F90 | 60 ++++++++++- schemes/musica/musica_ccpp.F90 | 100 +++++++++++++----- .../musica_ccpp_tuvx_gas_species_profile.F90 | 44 ++------ 3 files changed, 138 insertions(+), 66 deletions(-) diff --git a/schemes/musica/micm/musica_ccpp_micm.F90 b/schemes/musica/micm/musica_ccpp_micm.F90 index d2c6044d..5ba9eaf8 100644 --- a/schemes/musica/micm/musica_ccpp_micm.F90 +++ b/schemes/musica/micm/musica_ccpp_micm.F90 @@ -24,10 +24,14 @@ subroutine micm_register(solver_type, number_of_grid_cells, constituent_props, & use musica_micm, only: Rosenbrock, RosenbrockStandardOrder use musica_util, only: error_t use iso_c_binding, only: c_int + use musica_ccpp_gas_species, only: species_t, species_group_t, & + species_profiled_t, species_profiled_groupt_t,INDEX_UNINITIALIZED integer(c_int), intent(in) :: solver_type integer(c_int), intent(in) :: number_of_grid_cells type(ccpp_constituent_properties_t), allocatable, intent(out) :: constituent_props(:) + type(species_groupt_t), allocatable, intent(out) :: species_group(:) + type(species_profiled_group_t), allocatable, intent(out) :: species_profiled_group(:) character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode @@ -35,20 +39,36 @@ subroutine micm_register(solver_type, number_of_grid_cells, constituent_props, & type(error_t) :: error real(kind=kind_phys) :: molar_mass character(len=:), allocatable :: species_name - logical :: is_advected - integer :: i, species_index + logical :: is_advected, is_profiled + integer :: i, i_species_profiled + integer :: species_index + integer :: number_of_species + real(kind=kind_phys) :: species_unit, species_scale_height + + type(species_profiled_group_t), allocatable, :: temp_species_profiled_group(:) micm => micm_t(trim(filename_of_micm_configuration), solver_type, & number_of_grid_cells, error) if (has_error_occurred(error, errmsg, errcode)) return - allocate(constituent_props(micm%species_ordering%size()), stat=errcode) + number_of_species = micm%species_ordering%size() + + allocate(constituent_props(number_of_species), stat=errcode) if (errcode /= 0) then errmsg = "[MUSICA Error] Failed to allocate memory for constituent properties." return end if - do i = 1, micm%species_ordering%size() + allocate(temp_species_profiled_group(number_of_species), stat=errcode) + if (errcode /= 0) then + errmsg = "[MUSICA Error] Failed to allocate memory for musica species group." + return + end if + + i_species_profiled = 0 + + ! Sets constituents, creates species_t and species_profiled_t + do i = 1, number_of_species associate( map => micm%species_ordering ) species_name = map%name(i) species_index = map%index(i) @@ -74,8 +94,40 @@ subroutine micm_register(solver_type, number_of_grid_cells, constituent_props, & errcode = errcode, & errmsg = errmsg) if (errcode /= 0) return + + ! set musica species + musica_species_group(i) = musica_species_t(species_name, molar_mass, & + species_index, INDEX_UNINITIALIZED) + + is_profiled = micm%get_species_property_bool(species_name, & + "__is profiled", error) + if (has_error_occurred(error, errmsg, errcode)) return + + if (is_profiled) then + i_species_profiled = i_species_profiled + 1 + + species_unit = micm%get_species_property_double(species_name, & + "__unit", error) + if (has_error_occurred(error, errmsg, errcode)) return + + species_scale_height = micm%get_species_property_double(species_name, & + "__scale height [km]", error) + if (has_error_occurred(error, errmsg, errcode)) return + + musica_species_profiled_group(i_species_profiled) = & + musica_species_profiled_t( species_name, species_unit, molar_mass, & + species_scale_height, species_index, INDEX_UNINITIALIZED) + end if + end associate ! map end do + + if (i_species_profiled /= 0) then + allocate( species_profiled_group(i_species_profiled) ) + species_profiled_group(:) = temp_species_profiled_group(1:i_species_profiled) + end if + deallocate(temp_species_profiled_group) + number_of_rate_parameters = micm%user_defined_reaction_rates%size() end subroutine micm_register diff --git a/schemes/musica/musica_ccpp.F90 b/schemes/musica/musica_ccpp.F90 index d9172703..55ec4806 100644 --- a/schemes/musica/musica_ccpp.F90 +++ b/schemes/musica/musica_ccpp.F90 @@ -10,24 +10,28 @@ module musica_ccpp public :: musica_ccpp_register, musica_ccpp_init, musica_ccpp_run, musica_ccpp_final + type(musica_sequence_t) :: musica_sequence contains !> \section arg_table_musica_ccpp_register Argument Table !! \htmlinclude musica_ccpp_register.html subroutine musica_ccpp_register(micm_solver_type, number_of_grid_cells, & - constituent_props, errmsg, errcode) + constituent_props, species_group, species_profiled_group, errmsg, errcode) use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t - + use musica_ccpp_gas_species, only: species_groupt_t, species_profiled_group_t + integer, intent(in) :: micm_solver_type integer, intent(in) :: number_of_grid_cells type(ccpp_constituent_properties_t), allocatable, intent(out) :: constituent_props(:) + type(species_groupt_t), allocatable, intent(out) :: species_group(:) + type(species_profiled_group_t), allocatable, intent(out) :: species_profiled_group(:) character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode type(ccpp_constituent_properties_t), allocatable :: constituent_props_subset(:) call micm_register(micm_solver_type, number_of_grid_cells, constituent_props_subset, & - errmsg, errcode) + species_group, species_profiled_group, errmsg, errcode) if (errcode /= 0) return constituent_props = constituent_props_subset deallocate(constituent_props_subset) @@ -44,9 +48,11 @@ subroutine musica_ccpp_init(vertical_layer_dimension, vertical_interface_dimensi photolysis_wavelength_grid_interfaces, & constituent_props, errmsg, errcode) use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t - use ccpp_kinds, only : kind_phys - use musica_ccpp_micm, only: micm - use musica_ccpp_util, only: has_error_occurred + use ccpp_kinds, only: kind_phys + use ccpp_const_utils, only: ccpp_const_get_idx + use musica_ccpp_micm, only: micm + use musica_ccpp_util, only: has_error_occurred + integer, intent(in) :: vertical_layer_dimension ! (count) integer, intent(in) :: vertical_interface_dimension ! (count) real(kind_phys), intent(in) :: photolysis_wavelength_grid_interfaces(:) ! m @@ -54,6 +60,31 @@ subroutine musica_ccpp_init(vertical_layer_dimension, vertical_interface_dimensi character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode + !!!!!!!!! JIWON + ! Search for the index based on the species name + do i_species = 1, size(species) + name = species%name + call ccpp_const_get_idx(constituent_props, name, index_constituents, errmsg, errcode) + if (errcode /= 0) return + species%index_constituents = index_constituents + array(i_species) = index_constituents + end do + + ! Get the molar mass that is set in the call to instantiate() + do i_species = 1, size(species) + call constituent_props(array(i_species))%molar_mass(molar_mass_arr(i_elem), errcode, errmsg) + if (errcode /= 0) then + errmsg = "[MUSICA Error] Unable to get molar mass." + return + end if + end do + !!!!!!!!! + !!!!!!!!! + + ! set molar mass, set species from tuvx + call set_musica_species() + + call micm_init(errmsg, errcode) if (errcode /= 0) return call tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & @@ -111,25 +142,35 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co number_of_rate_parameters) :: rate_parameters ! various units integer :: i_elem - ! TODO(jiwon) this might not be correct because it doesn't know the index - ! Get the molar mass that is set in the call to instantiate() - do i_elem = 1, size(molar_mass_arr) - call constituent_props(i_elem)%molar_mass(molar_mass_arr(i_elem), errcode, errmsg) - if (errcode /= 0) then - errmsg = "[MUSICA Error] Unable to get molar mass." - return - end if - end do - - ! TODO(jiwon) Check molar mass is non zero as it becomes a denominator for unit converison - ! this code will be deleted when the framework does the check - do i_elem = 1, size(molar_mass_arr) - if (molar_mass_arr(i_elem) <= 0) then - errcode = 1 - errmsg = "[MUSICA Error] Molar mass must be greater than zero." - return - end if - end do + real(kind_phys), dimension(size(constituents, dim=1), & + size(constituents, dim=2), & + number_of_musica_speices) :: constituents_musica + + ! ! TODO(jiwon) this might not be correct because it doesn't know the index + ! ! Get the molar mass that is set in the call to instantiate() + ! do i_elem = 1, size(molar_mass_arr) + ! call constituent_props(i_elem)%molar_mass(molar_mass_arr(i_elem), errcode, errmsg) + ! if (errcode /= 0) then + ! errmsg = "[MUSICA Error] Unable to get molar mass." + ! return + ! end if + ! end do + + ! ! TODO(jiwon) Check molar mass is non zero as it becomes a denominator for unit converison + ! ! this code will be deleted when the framework does the check + ! do i_elem = 1, size(molar_mass_arr) + ! if (molar_mass_arr(i_elem) <= 0) then + ! errcode = 1 + ! errmsg = "[MUSICA Error] Molar mass must be greater than zero." + ! return + ! end if + ! end do + + ! ! Convert CAM-SIMA unit to MICM unit (kg kg-1 -> mol m-3) + ! call convert_to_mol_per_cubic_meter(dry_air_density, molar_mass_arr, constituents) + + ! find musica constituents array + call get_musica_constituents(constituents, constituents_musica, sequence_musica_species) ! Convert CAM-SIMA unit to MICM unit (kg kg-1 -> mol m-3) call convert_to_mol_per_cubic_meter(dry_air_density, molar_mass_arr, constituents) @@ -144,17 +185,20 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co photolysis_wavelength_grid_interfaces, & extraterrestrial_flux, & standard_gravitational_acceleration, & - cloud_area_fraction, constituents, & + cloud_area_fraction, constituents_musica, & ! need to reset the arrangement air_pressure_thickness, rate_parameters, & errmsg, errcode) ! Solve chemistry at the current time step call micm_run(time_step, temperature, pressure, dry_air_density, rate_parameters, & - constituents, errmsg, errcode) + constituents_musica, errmsg, errcode) ! 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) + ! return as CAM-SIMA constituents array + call convert_to_mass_mixing_ratio(dry_air_density, molar_mass_arr, constituents_musica) + call update_musica_constituents(constituents_musica, constituents, sequence_musica_species) + end subroutine musica_ccpp_run !> \section arg_table_musica_ccpp_final Argument Table diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 index cc58e500..8a13d230 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species_profile.F90 @@ -1,11 +1,11 @@ module musica_ccpp_tuvx_gas_species_profile use ccpp_kinds, only: kind_phys use musica_tuvx_profile, only: profile_t - + use musica_ccpp_gas_species, only: gas_species_t implicit none private - public :: create_gas_species, configure_gas_species, create_gas_species_profile_group, & + public :: configure_gas_species, create_gas_species_profile_group, & set_gas_species_values, deallocate_gas_species, deallocate_gas_species_profile_group !> Conversion factor from km to cm @@ -16,39 +16,12 @@ module musica_ccpp_tuvx_gas_species_profile ! TODO(jiwon) real(kind_phys), parameter, public :: MOLAR_MASS_DRY_AIR = 0.0289644_kind_phys ! kg mol-1 - !> Definition of the gas species type - type, public :: gas_species_t - character(len=50) :: label - character(len=20) :: unit - real(kind_phys) :: molar_mass ! kg mol-1 - real(kind_phys) :: scale_height ! km - integer :: index_constituent_props ! index of the constituents - end type gas_species_t - !> Defines a type to contaion a pointer to 'profile_t' type, public :: profile_group_t type(profile_t), pointer :: profile end type profile_group_t contains - !> Creates 'gas_species_t' object - function create_gas_species(label, unit, molar_mass, scale_height, & - index_constituent_props) result( gas_species ) - - character(len=*), intent(in) :: label - character(len=*), intent(in) :: unit - real(kind_phys), intent(in) :: molar_mass ! kg mol-1 - real(kind_phys), intent(in) :: scale_height ! km - integer, intent(in) :: index_constituent_props - type(gas_species_t) :: gas_species - - gas_species%label = label - gas_species%unit = unit - gas_species%molar_mass = molar_mass - gas_species%scale_height = scale_height - gas_species%index_constituent_props = index_constituent_props - - end function create_gas_species !> Configures gas species. ! Note: The user can customize this function to configure gas species in the system @@ -89,8 +62,13 @@ subroutine configure_gas_species(constituent_props, gas_species_group, errmsg, e ! call constituent_props(index_air)%molar_mass(molar_mass_air, errcode, errmsg) ! if (errcode /= 0) return + is_advected = micm%get_species_property_bool(species_name, & + "__is advected", & + error) + + index_air = INDEX_NOT_KNOWN - gas_species_group(1) = create_gas_species("air", "molecule cm-3", & + gas_species_group(1) = gas_species_t("air", "molecule cm-3", & MOLAR_MASS_DRY_AIR, SCALE_HEIGHT_DRY_AIR, index_air) ! O2 call ccpp_const_get_idx(constituent_props, "O2", index_O2, errmsg, errcode) @@ -99,7 +77,7 @@ subroutine configure_gas_species(constituent_props, gas_species_group, errmsg, e call constituent_props(index_O2)%molar_mass(molar_mass_O2, errcode, errmsg) if (errcode /= 0) return - gas_species_group(2) = create_gas_species("O2", "molecule cm-3", & + gas_species_group(2) = gas_species_t("O2", "molecule cm-3", & molar_mass_O2, SCALE_HEIGHT_O2, index_O2) ! O3 @@ -109,7 +87,7 @@ subroutine configure_gas_species(constituent_props, gas_species_group, errmsg, e call constituent_props(index_O3)%molar_mass(molar_mass_O3, errcode, errmsg) if (errcode /= 0) return - gas_species_group(3) = create_gas_species("O3", "molecule cm-3", & + gas_species_group(3) = gas_species_t("O3", "molecule cm-3", & molar_mass_O3, SCALE_HEIGHT_O3, index_O3) end subroutine configure_gas_species @@ -163,8 +141,6 @@ subroutine set_gas_species_values(profile, gas_species, constituent, & num_constituents = size(constituent) constituent_mol_per_cm_3(:) = constituent(:) / m_3_to_cm_3 - - interfaces(1) = constituent_mol_per_cm_3(num_constituents) interfaces(2:num_constituents+1) = constituent_mol_per_cm_3(num_constituents:1:-1) interfaces(num_constituents+2) = constituent_mol_per_cm_3(1)