Add TUV-x height grid updates

.github/workflows/test.yaml

@@ -12,4 +12,14 @@ jobs:
     - name: build Docker image
       run: docker build -t musica -f test/docker/Dockerfile.musica .
     - name: run tests in container
-      run: docker run --name test-container -t musica bash -c 'make test ARGS="--rerun-failed --output-on-failure -j8"'
+      run: docker run --name test-container -t musica bash -c 'make test ARGS="--rerun-failed --output-on-failure -j8"'
+  test_musica_api_no_install:
+    runs-on: ubuntu-latest
+    steps:
+    - uses: actions/checkout@v4
+      with:
+        submodules: recursive
+    - name: build Docker image
+      run: docker build -t musica-no-install -f test/docker/Dockerfile.musica.no_install .
+    - name: run tests in container
+      run: docker run --name test-container -t musica-no-install bash -c 'make test ARGS="--rerun-failed --output-on-failure -j8"'

schemes/musica/micm/musica_ccpp_micm.F90

@@ -30,37 +30,41 @@ subroutine micm_register(constituents, solver_type, num_grid_cells, errmsg, errc
     integer,                                          intent(out) :: errcode
 
     ! local variables
-    type(error_t)        :: error
-    real(kind=kind_phys) :: molar_mass
-    logical              :: is_advected
-    integer              :: i
+    type(error_t)                 :: error
+    real(kind=kind_phys)          :: molar_mass
+    character(len=:), allocatable :: species_name
+    logical                       :: is_advected
+    integer                       :: i, species_index
 
     errcode = 0
     errmsg = ''
 
     micm => micm_t(filename_of_micm_configuration, solver_type, num_grid_cells, error)
     if (has_error_occurred(error, errmsg, errcode)) return
 
-    allocate(constituents(size(micm%species_ordering)), stat=errcode)
+    allocate(constituents(micm%species_ordering%size()), stat=errcode)
     if (errcode /= 0) then
       errmsg = "[MUSICA Error] Failed to allocate memory for constituents."
       return
     end if
 
-    do i = 1, size(micm%species_ordering)
-    associate( map => micm%species_ordering(i) )
-      molar_mass = micm%get_species_property_double(map%name(), &
+    do i = 1, micm%species_ordering%size()
+    associate( map => micm%species_ordering )
+      species_name = map%name(i)
+      species_index = map%index(i)
+
+      molar_mass = micm%get_species_property_double(species_name, &
                                                     "molecular weight [kg mol-1]", &
                                                     error)
       if (has_error_occurred(error, errmsg, errcode)) return
-      is_advected = micm%get_species_property_bool(map%name(), &
-                                                      "__is advected", &
-                                                      error)
+      is_advected = micm%get_species_property_bool(species_name, &
+                                                   "__is advected", &
+                                                   error)
       if (has_error_occurred(error, errmsg, errcode)) return
 
-      call constituents(map%index())%instantiate( &
-        std_name = map%name(), &
-        long_name = map%name(), &
+      call constituents(species_index)%instantiate( &
+        std_name = species_name, &
+        long_name = species_name, &
         units = 'kg kg-1', &
         vertical_dim = 'vertical_layer_dimension', &
         default_value = 0.0_kind_phys, &
@@ -87,7 +91,7 @@ end subroutine micm_init
 
   !> Solve chemistry at the current time step
   subroutine micm_run(time_step, temperature, pressure, dry_air_density, constituents, &
-                      rate_params, errmsg, errcode)
+                      user_defined_rate_parameters, errmsg, errcode)
     use musica_micm, only: solver_stats_t
     use musica_util, only: string_t, error_t
 
@@ -96,7 +100,7 @@ subroutine micm_run(time_step, temperature, pressure, dry_air_density, constitue
     real(c_double), target, intent(in)    :: pressure(:)        ! Pa
     real(c_double), target, intent(in)    :: dry_air_density(:) ! kg m-3
     real(c_double), target, intent(inout) :: constituents(:)    ! mol m-3
-    real(c_double), target, intent(inout) :: rate_params(:)
+    real(c_double), target, intent(in)    :: user_defined_rate_parameters(:) ! various units
     character(len=512),     intent(out)   :: errmsg
     integer,                intent(out)   :: errcode
 
@@ -111,14 +115,14 @@ subroutine micm_run(time_step, temperature, pressure, dry_air_density, constitue
     errmsg = ''
     c_time_step = real(time_step, c_double) 
 
-    call micm%solve(c_time_step,     &
-                    temperature,     &
-                    pressure,        &
-                    dry_air_density, &
-                    constituents,    &
-                    rate_params,     &
-                    solver_state,    &
-                    solver_stats,    &
+    call micm%solve(c_time_step,                  &
+                    temperature,                  &
+                    pressure,                     &
+                    dry_air_density,              &
+                    constituents,                 &
+                    user_defined_rate_parameters, &
+                    solver_state,                 &
+                    solver_stats,                 &
                     error)
     if (has_error_occurred(error, errmsg, errcode)) return
 

schemes/musica/micm/micm_util.F90 → schemes/musica/micm/musica_ccpp_micm_util.F90

@@ -1,4 +1,4 @@
-module micm_util
+module musica_ccpp_micm_util
   implicit none
 
   private
@@ -8,93 +8,67 @@ module micm_util
 contains
 
   subroutine reshape_into_micm_arr(temperature, pressure, dry_air_density, constituents, &
-      rate_params, m_temperature, m_pressure, m_dry_air_density, m_constituents, m_rate_params)
+      m_temperature, m_pressure, m_dry_air_density, m_constituents)
     use iso_c_binding, only: c_double
     use ccpp_kinds,    only: kind_phys
 
     real(kind_phys), target, intent(in)  :: temperature(:,:)     ! K
     real(kind_phys), target, intent(in)  :: pressure(:,:)        ! Pa
     real(kind_phys), target, intent(in)  :: dry_air_density(:,:) ! kg m-3
     real(kind_phys), target, intent(in)  :: constituents(:,:,:)  ! kg kg-1
-    real(kind_phys), target, intent(in)  :: rate_params(:,:,:)
     real(c_double),  target, intent(out) :: m_temperature(:)     ! K
     real(c_double),  target, intent(out) :: m_pressure(:)        ! Pa
     real(c_double),  target, intent(out) :: m_dry_air_density(:) ! kg m-3
     real(c_double),  target, intent(out) :: m_constituents(:)    ! kg kg-1
-    real(c_double),  target, intent(out) :: m_rate_params(:)
 
     ! local variables
     integer :: num_columns, num_layers
-    integer :: num_constituents, num_rate_params
-    integer :: i_column, i_layer, i_elem, i_constituents, i_rate_params
+    integer :: num_constituents
+    integer :: i_column, i_layer, i_elem, i_constituents
 
     num_columns = size(constituents, dim=1)
     num_layers = size(constituents, dim=2)
     num_constituents = size(constituents, dim=3)
-    num_rate_params = size(rate_params, dim=3)
 
     ! Reshape into 1-D arry in species-column first order
     ! refers to: state.variables_[i_cell][i_species] = concentrations[i_species_elem++]
     i_elem = 1
     i_constituents = 1
-    i_rate_params = 1
     do i_layer = 1, num_layers
       do i_column = 1, num_columns
         m_temperature(i_elem) = real(temperature(i_column, i_layer), c_double)
         m_pressure(i_elem) = real(pressure(i_column, i_layer), c_double)
         m_dry_air_density(i_elem) = real(dry_air_density(i_column, i_layer), c_double)
         m_constituents(i_constituents : i_constituents + num_constituents - 1) &
                             = real(constituents(i_column, i_layer, :), c_double)
-        m_rate_params(i_rate_params : i_rate_params + num_rate_params - 1) &
-                            = real(rate_params(i_column, i_layer, :), c_double)
         i_elem = i_elem + 1
         i_constituents = i_constituents + num_constituents
-        i_rate_params = i_rate_params + num_rate_params
       end do
     end do
 
   end subroutine reshape_into_micm_arr
 
-  subroutine reshape_into_ccpp_arr(temperature, pressure, dry_air_density, constituents, &
-      rate_params, m_temperature, m_pressure, m_dry_air_density, m_constituents, m_rate_params)
+  subroutine reshape_into_ccpp_arr(constituents, m_constituents)
     use iso_c_binding, only: c_double
     use ccpp_kinds,    only: kind_phys
-    real(kind_phys), intent(out) :: temperature(:,:)     ! K
-    real(kind_phys), intent(out) :: pressure(:,:)        ! Pa
-    real(kind_phys), intent(out) :: dry_air_density(:,:) ! kg m-3
     real(kind_phys), intent(out) :: constituents(:,:,:)  ! kg kg-1
-    real(kind_phys), intent(out) :: rate_params(:,:,:)
-    real(c_double),  intent(in)  :: m_temperature(:)     ! K
-    real(c_double),  intent(in)  :: m_pressure(:)        ! Pa
-    real(c_double),  intent(in)  :: m_dry_air_density(:) ! kg m-3
     real(c_double),  intent(in)  :: m_constituents(:)    ! kg kg-1
-    real(c_double),  intent(in)  :: m_rate_params(:)
 
     ! local variables
     integer :: num_columns, num_layers
-    integer :: num_constituents, num_rate_params
-    integer :: i_column, i_layer, i_elem, i_constituents, i_rate_params
+    integer :: num_constituents
+    integer :: i_column, i_layer, i_constituents
 
     num_columns = size(constituents, dim=1)
     num_layers = size(constituents, dim=2)
     num_constituents = size(constituents, dim=3)
-    num_rate_params = size(rate_params, dim=3)
 
-    i_elem = 1
     i_constituents = 1
-    i_rate_params = 1
     do i_layer = 1, num_layers
       do i_column = 1, num_columns
-        temperature(i_column, i_layer) = real(m_temperature(i_elem), kind_phys)
-        pressure(i_column, i_layer) = real(m_pressure(i_elem), kind_phys)
-        dry_air_density(i_column, i_layer) = real(m_dry_air_density(i_elem), kind_phys)
         constituents(i_column, i_layer, :) &
             = real(m_constituents(i_constituents : i_constituents + num_constituents - 1), kind_phys)
-        rate_params(i_column, i_layer, :) &
-            = real(m_rate_params(i_rate_params : i_rate_params + num_rate_params - 1), kind_phys)
-        i_elem = i_elem + 1
         i_constituents = i_constituents + num_constituents
-        i_rate_params = i_rate_params + num_rate_params
       end do
     end do
 
@@ -156,4 +130,4 @@ subroutine convert_to_mass_mixing_ratio(dry_air_density, molar_mass_arr, constit
 
   end subroutine convert_to_mass_mixing_ratio
 
-end module micm_util
+end module musica_ccpp_micm_util

schemes/musica/musica_ccpp.F90

@@ -25,32 +25,43 @@ end subroutine musica_ccpp_register
 
   !> \section arg_table_musica_ccpp_init Argument Table
   !! \htmlinclude musica_ccpp_init.html
-  subroutine musica_ccpp_init(errmsg, errcode)
+  subroutine musica_ccpp_init(vertical_layer_dimension, vertical_interface_dimension, &
+      errmsg, errcode)
+    integer,            intent(in)  :: vertical_layer_dimension     ! (count)
+    integer,            intent(in)  :: vertical_interface_dimension ! (count)
     character(len=512), intent(out) :: errmsg
     integer,            intent(out) :: errcode
 
-    call tuvx_init(errmsg, errcode)
+    call tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
+                   errmsg, errcode)
+    if (errcode /= 0) return
     call micm_init(errmsg, errcode)
+    if (errcode /= 0) return
 
   end subroutine musica_ccpp_init
 
   !> \section arg_table_musica_ccpp_run Argument Table
   !! \htmlinclude musica_ccpp_run.html
   subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, constituent_props, &
-                      constituents, rate_params, height, errmsg, errcode)
-    use micm_util,                 only: reshape_into_micm_arr, reshape_into_ccpp_arr
-    use micm_util,                 only: convert_to_mol_per_cubic_meter, convert_to_mass_mixing_ratio
+                      constituents, geopotential_height_wrt_surface_at_midpoint, &
+                      geopotential_height_wrt_surface_at_interface, &
+                      surface_geopotential, reciprocal_of_gravitational_acceleration, &
+                      errmsg, errcode)
+    use musica_ccpp_micm_util,     only: reshape_into_micm_arr, reshape_into_ccpp_arr
+    use musica_ccpp_micm_util,     only: convert_to_mol_per_cubic_meter, convert_to_mass_mixing_ratio
     use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
     use ccpp_kinds,                only: kind_phys
     use iso_c_binding,             only: c_double
-    real(kind_phys),                   intent(inout) :: time_step            ! s
-    real(kind_phys), target,           intent(inout) :: temperature(:,:)     ! K
-    real(kind_phys), target,           intent(inout) :: pressure(:,:)        ! Pa
-    real(kind_phys), target,           intent(inout) :: dry_air_density(:,:) ! kg m-3
+    real(kind_phys),                   intent(in)    :: time_step              ! s
+    real(kind_phys), target,           intent(in)    :: temperature(:,:)       ! K
+    real(kind_phys), target,           intent(in)    :: pressure(:,:)          ! Pa
+    real(kind_phys), target,           intent(in)    :: dry_air_density(:,:)   ! kg m-3
     type(ccpp_constituent_prop_ptr_t), intent(in)    :: constituent_props(:)
-    real(kind_phys), target,           intent(inout) :: constituents(:,:,:)  ! kg kg-1
-    real(kind_phys), target,           intent(inout) :: rate_params(:,:,:)
-    real(kind_phys), target,           intent(in)    :: height(:,:)          ! km
+    real(kind_phys), target,           intent(inout) :: constituents(:,:,:)    ! kg kg-1
+    real(kind_phys), target,           intent(in)    :: geopotential_height_wrt_surface_at_midpoint(:,:)  ! m
+    real(kind_phys), target,           intent(in)    :: geopotential_height_wrt_surface_at_interface(:,:) ! m
+    real(kind_phys), target,           intent(in)    :: surface_geopotential(:)    ! m2 s-2
+    real(kind_phys), target,           intent(in)    :: reciprocal_of_gravitational_acceleration ! s2 m-1
     character(len=512),                intent(out)   :: errmsg
     integer,                           intent(out)   :: errcode
 
@@ -64,13 +75,18 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
     real(c_double), target, dimension(size(constituents, dim=1)    &
                                     * size(constituents, dim=2)    & 
                                     * size(constituents, dim=3))     :: m_constituents ! mol m-3
-    real(c_double), target, dimension(size(rate_params, dim=1)     &
-                                    * size(rate_params, dim=2)     & 
-                                    * size(rate_params, dim=3))      :: m_rate_params
     real(kind_phys), target, dimension(size(constituents, dim=3))    :: molar_mass_arr ! kg mol-1
+
+    ! temporarily dimensioned to Chapman mechanism until mapping between MICM and TUV-x is implemented
+    real(c_double), target, dimension(size(constituents, dim=1) &
+                                    * size(constituents, dim=2) &
+                                    * 3) :: photolysis_rate_constants ! s-1
     integer :: i_elem
 
-    call tuvx_run(height, temperature, dry_air_density, errmsg, errcode)
+    call tuvx_run(temperature, dry_air_density, geopotential_height_wrt_surface_at_midpoint, &
+                  geopotential_height_wrt_surface_at_interface, &
+                  surface_geopotential, reciprocal_of_gravitational_acceleration, &
+                  photolysis_rate_constants, errmsg, errcode)
 
     ! Get the molar_mass that is set in the call to instantiate()
     do i_elem = 1, size(molar_mass_arr)
@@ -95,15 +111,15 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
     call convert_to_mol_per_cubic_meter(dry_air_density, molar_mass_arr, constituents)
 
     ! Reshape array (3D -> 1D) and convert type (kind_phys -> c_double)
-    call reshape_into_micm_arr(temperature, pressure, dry_air_density, constituents, rate_params, &
-                      m_temperature, m_pressure, m_dry_air_density, m_constituents, m_rate_params)
+    call reshape_into_micm_arr(temperature, pressure, dry_air_density, constituents, &
+                      m_temperature, m_pressure, m_dry_air_density, m_constituents)
 
+    ! temporarily pass in unmapped photolysis rate constants until mapping between MICM and TUV-x is implemented
     call micm_run(time_step, m_temperature, m_pressure, m_dry_air_density, m_constituents, &
-                      m_rate_params, errmsg, errcode)
+                  photolysis_rate_constants, errmsg, errcode)
 
     ! Reshape array (1D -> 3D) and convert type (c_double -> kind_phys)
-    call reshape_into_ccpp_arr(temperature, pressure, dry_air_density, constituents, rate_params, &
-                      m_temperature, m_pressure, m_dry_air_density, m_constituents, m_rate_params)
+    call reshape_into_ccpp_arr(constituents, m_constituents)
 
     ! Convert MICM unit back to CAM-SIMA unit (mol m-3  ->  kg kg-1)
     call convert_to_mass_mixing_ratio(dry_air_density, molar_mass_arr, constituents)

schemes/musica/musica_ccpp.meta

@@ -1,12 +1,24 @@
 [ccpp-table-properties]
   name = musica_ccpp
   type = scheme
-  dependencies = micm/musica_ccpp_micm.F90,micm/musica_ccpp_tuvx.F90,musica_ccpp_util.F90
+  dependencies = micm/musica_ccpp_micm.F90,micm/musica_ccpp_micm_util.F90,tuvx/musica_ccpp_tuvx.F90,tuvx/musica_ccpp_tuvx_height_grid.F90,musica_ccpp_util.F90
   dynamic_constituent_routine = musica_ccpp_register
 
 [ccpp-arg-table]
   name  = musica_ccpp_init
   type  = scheme
+[ vertical_layer_dimension ]
+  standard_name = vertical_layer_dimension
+  units = none
+  type = integer
+  dimensions = ()
+  intent = in
+[ vertical_interface_dimension ]
+  standard_name = vertical_interface_dimension
+  units = none
+  type = integer
+  dimensions = ()
+  intent = in
 [ errmsg ]
   standard_name = ccpp_error_message
   units = none
@@ -59,6 +71,30 @@
   type = real | kind = kind_phys
   dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_ccpp_constituents)
   intent = inout
+[ geopotential_height_wrt_surface_at_midpoint ]
+  standard_name = geopotential_height_wrt_surface
+  units = km
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  intent = in
+[ geopotential_height_wrt_surface_at_interface ]
+  standard_name = geopotential_height_wrt_surface_at_interface
+  units = km
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  intent = in
+[ surface_geopotential ]
+  standard_name = surface_geopotential
+  type = real | kind = kind_phys
+  units = m2 s-2
+  dimensions = (horizontal_loop_extent)
+  intent = in
+[ reciprocal_of_gravitational_acceleration ]
+  standard_name = reciprocal_of_gravitational_acceleration
+  units = s2 m-1
+  type = real | kind = kind_phys
+  dimensions = ()
+  intent = in
 [ errmsg ]
   standard_name = ccpp_error_message
   units = none