Add ability to calculate heating rates for photolysis reactions (#50)

* draft heating rates class and tests

* finish heating rate tests

* add heating rate output to driver, and include in TS example

examples/ts1_tsmlt.json

@@ -199,6 +199,9 @@
                   "value": 1.0
                 }
               ]
+            },
+            "heating" : {
+               "energy term": 175.05
             }
          },
          {
@@ -228,6 +231,9 @@
                    "value": 0.0
                  }
                ]
+            },
+            "heating" : {
+               "energy term": 242.37
             }
          },
          {
@@ -244,6 +250,9 @@
             },
             "quantum yield": {
                "type": "O3+hv->O2+O(1D)"
+            },
+            "heating" : {
+               "energy term": 310.32
             }
          },
          {
@@ -260,6 +269,9 @@
             },
             "quantum yield": {
                "type": "O3+hv->O2+O(3P)"
+            },
+            "heating" : {
+               "energy term": 1179.87
             }
          },
          {

examples/ts1_tsmlt.yml

@@ -210,6 +210,8 @@ photolysis:
         lower extrapolation:
           type: boundary
       type: base
+    heating:
+      energy term: 175.05
     name: jo2_a
     quantum yield:
       constant value: 0
@@ -229,6 +231,8 @@ photolysis:
         lower extrapolation:
           type: boundary
       type: base
+    heating:
+      energy term: 242.37
     name: jo2_b
     quantum yield:
       constant value: 1.0
@@ -246,6 +250,8 @@ photolysis:
       - file path: data/cross_sections/O3_3.nc
       - file path: data/cross_sections/O3_4.nc
       type: O3
+    heating:
+      energy term: 310.32
     name: jo3_a
     quantum yield:
       type: O3+hv->O2+O(1D)
@@ -257,6 +263,8 @@ photolysis:
       - file path: data/cross_sections/O3_3.nc
       - file path: data/cross_sections/O3_4.nc
       type: O3
+    heating:
+      energy term: 1179.87
     name: jo3_b
     quantum yield:
       type: O3+hv->O2+O(3P)

src/CMakeLists.txt

@@ -46,6 +46,7 @@ target_sources(tuvx_object
     grid.F90
     grid_factory.F90
     grid_warehouse.F90
+    heating_rates.F90
     interpolate.F90
     la_sr_bands.F90
     linear_algebra.F90

src/core.F90

@@ -9,6 +9,7 @@ module tuvx_core
   use musica_constants,                only : dk => musica_dk
   use tuvx_dose_rates,                 only : dose_rates_t
   use tuvx_grid_warehouse,             only : grid_warehouse_t
+  use tuvx_heating_rates,              only : heating_rates_t
   use tuvx_la_sr_bands,                only : la_sr_bands_t
   use tuvx_photolysis_rates,           only : photolysis_rates_t
   use tuvx_profile_warehouse,          only : profile_warehouse_t
@@ -32,10 +33,11 @@ module tuvx_core
     type(radiative_transfer_t),  pointer :: radiative_transfer_ => null()
     type(photolysis_rates_t),    pointer :: photolysis_rates_ => null()
     type(dose_rates_t),          pointer :: dose_rates_ => null()
+    type(heating_rates_t),       pointer :: heating_rates_ => null()
     type(radiation_field_t),     pointer :: radiation_field_ => null()
     logical                              :: enable_diagnostics_ ! determines if diagnostic output is written or not
   contains
-    ! Calculate photolysis rate constants and dose rates
+    ! Calculate photolysis rate constants, dose rates, and heating rates
     procedure :: run
     ! Returns a grid from the warehouse
     procedure :: get_grid
@@ -50,10 +52,14 @@ module tuvx_core
     procedure :: number_of_photolysis_reactions
     ! Returns the number of dose rates
     procedure :: number_of_dose_rates
+    ! Returns the number of heating rates
+    procedure :: number_of_heating_rates
     ! Returns the set of photolysis reaction labels
     procedure :: photolysis_reaction_labels
     ! Returns the set of dose rate labels
     procedure :: dose_rate_labels
+    ! Returns the set of heating rate labels
+    procedure :: heating_rate_labels
     ! Returns the photolysis reaction cross section for the current conditions
     procedure :: get_photolysis_cross_section
     ! Returns the photolysis reaction quantum yield for the current conditions
@@ -165,6 +171,9 @@ function constructor( config, grids, profiles, radiators ) result( new_core )
           photolysis_rates_t( child_config,                                   &
                               new_core%grid_warehouse_,                       &
                               new_core%profile_warehouse_ )
+      new_core%heating_rates_ => heating_rates_t( child_config,             &
+                                                  new_core%grid_warehouse_, &
+                                                  new_core%profile_warehouse_ )
     end if
 
     ! dose rates
@@ -191,7 +200,7 @@ end function constructor
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
   subroutine run( this, solar_zenith_angle, earth_sun_distance,               &
-      photolysis_rate_constants, dose_rates, diagnostic_label )
+      photolysis_rate_constants, dose_rates, heating_rates, diagnostic_label )
     ! Performs calculations for specified photolysis and dose rates for a
     ! given set of conditions
 
@@ -205,6 +214,7 @@ subroutine run( this, solar_zenith_angle, earth_sun_distance,               &
     real(dk),                   intent(in)    :: earth_sun_distance             ! [AU]
     real(dk),         optional, intent(out)   :: photolysis_rate_constants(:,:) ! (vertical level, reaction) [s-1]
     real(dk),         optional, intent(out)   :: dose_rates(:,:)                ! (vertical level, reaction) [s-1]
+    real(dk),         optional, intent(out)   :: heating_rates(:,:)             ! (vertical level, reaction) [J s-1]  
     character(len=*), optional, intent(in)    :: diagnostic_label               ! label used in diagnostic file names
 
     ! Local variables
@@ -247,6 +257,14 @@ subroutine run( this, solar_zenith_angle, earth_sun_distance,               &
                                        photolysis_rate_constants,             &
                                        diag_label )
     end if
+    if( associated( this%heating_rates_ ) .and. present( heating_rates ) ) then
+      call this%heating_rates_%get( this%la_sr_bands_,                        &
+                                    this%spherical_geometry_,                 &
+                                    this%grid_warehouse_,                     &
+                                    this%profile_warehouse_,                  &
+                                    this%radiation_field_,                    &
+                                    heating_rates )
+    end if
     if( associated( this%dose_rates_ ) .and. present( dose_rates ) ) then
       call this%dose_rates_%get( this%grid_warehouse_,                        &
                                  this%profile_warehouse_,                     &
@@ -410,6 +428,20 @@ integer function number_of_dose_rates( this )
 
   end function number_of_dose_rates
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  integer function number_of_heating_rates( this )
+    ! Returns the number of heating rates
+
+    class(core_t), intent(in) :: this
+
+    number_of_heating_rates = 0
+    if( associated( this%heating_rates_ ) ) then
+      number_of_heating_rates = this%heating_rates_%size( )
+    end if
+
+  end function number_of_heating_rates
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
   function photolysis_reaction_labels( this ) result( labels )
@@ -442,6 +474,22 @@ function dose_rate_labels( this ) result( labels )
 
   end function dose_rate_labels
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  function heating_rate_labels( this ) result( labels )
+    ! Returns the set of heating rate labels
+
+    class(core_t),  intent(in)  :: this
+    type(string_t), allocatable :: labels(:)
+
+    if( associated( this%heating_rates_ ) ) then
+      labels = this%heating_rates_%labels( )
+    else
+      allocate( labels( 0 ) )
+    end if
+
+  end function heating_rate_labels
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
   function get_photolysis_cross_section( this, reaction_label, found )        &
@@ -558,6 +606,11 @@ integer function pack_size( this, comm )
     if( associated( this%dose_rates_ ) ) then
       pack_size = pack_size + this%dose_rates_%pack_size( comm )
     end if
+    pack_size = pack_size +                                                   &
+      musica_mpi_pack_size( associated( this%heating_rates_ ), comm )
+    if( associated( this%heating_rates_ ) ) then
+      pack_size = pack_size + this%heating_rates_%pack_size( comm )
+    end if
 #else
     pack_size = 0
 #endif
@@ -617,6 +670,11 @@ subroutine mpi_pack( this, buffer, position, comm )
     if( associated( this%dose_rates_ ) ) then
       call this%dose_rates_%mpi_pack( buffer, position, comm )
     end if
+    call musica_mpi_pack( buffer, position,                                   &
+                          associated( this%heating_rates_ ), comm )
+    if( associated( this%heating_rates_ ) ) then
+      call this%heating_rates_%mpi_pack( buffer, position, comm )
+    end if
     call assert( 332208077, position - prev_pos <= this%pack_size( comm ) )
 #endif
 
@@ -676,6 +734,11 @@ subroutine mpi_unpack( this, buffer, position, comm )
       allocate( this%dose_rates_ )
       call this%dose_rates_%mpi_unpack( buffer, position, comm )
     end if
+    call musica_mpi_unpack( buffer, position, alloced, comm )
+    if( alloced ) then
+      allocate( this%heating_rates_ )
+      call this%heating_rates_%mpi_unpack( buffer, position, comm )
+    end if
     call assert( 332208077, position - prev_pos <= this%pack_size( comm ) )
 #endif
 
@@ -713,6 +776,9 @@ subroutine finalize( this )
     if( associated( this%radiation_field_ ) ) then
       deallocate( this%radiation_field_ )
     end if
+    if( associated( this%heating_rates_ ) ) then
+      deallocate( this%heating_rates_ )
+    end if
 
   end subroutine finalize
 

src/grid_warehouse.F90

@@ -152,10 +152,10 @@ function get_grid_string( this, name, units ) result( a_grid_ptr )
     use musica_string,                 only : string_t
     use tuvx_grid,                     only : grid_t
 
-    class(grid_warehouse_t), intent(inout) :: this ! This :f:type:`~tuvx_grid_warehouse/grid_warehouse_t`
-    type(string_t),          intent(in)    :: name ! The name of a grid, see :ref:`configuration-grids` for grid names
-    type(string_t),          intent(in)    :: units ! The units of the grid
-    class(grid_t), pointer                 :: a_grid_ptr ! The :f:type:`~tuvx_grid/grid_t` which matches the name passed in
+    class(grid_warehouse_t), intent(in) :: this ! This :f:type:`~tuvx_grid_warehouse/grid_warehouse_t`
+    type(string_t),          intent(in) :: name ! The name of a grid, see :ref:`configuration-grids` for grid names
+    type(string_t),          intent(in) :: units ! The units of the grid
+    class(grid_t), pointer              :: a_grid_ptr ! The :f:type:`~tuvx_grid/grid_t` which matches the name passed in
 
     a_grid_ptr => this%get_grid_char( name%to_char( ), units%to_char( ) )
 
@@ -169,9 +169,9 @@ function get_grid_ptr( this, ptr ) result( grid )
     use musica_assert,                 only : assert_msg
     use tuvx_grid,                     only : grid_t
 
-    class(grid_warehouse_t),  intent(inout) :: this ! This grid warehouse
-    type(grid_warehouse_ptr), intent(in)    :: ptr  ! Pointer to a grid in the warehouse
-    class(grid_t),            pointer       :: grid
+    class(grid_warehouse_t),  intent(in) :: this ! This grid warehouse
+    type(grid_warehouse_ptr), intent(in) :: ptr  ! Pointer to a grid in the warehouse
+    class(grid_t),            pointer    :: grid
 
     call assert_msg( 870082797, ptr%index_ > 0, "Invalid grid pointer" )
     allocate( grid, source = this%grids_( ptr%index_ )%val_ )