initial changes for chem potential energy rates

modified: src/core.F90 modified: src/photolysis_rates.F90
NCAR · Feb 1, 2024 · 06cb2bc · 06cb2bc
1 parent e1e7f1d
commit 06cb2bc
Show file tree

Hide file tree

Showing 2 changed files with 15 additions and 3 deletions.
diff --git a/src/core.F90 b/src/core.F90
@@ -191,7 +191,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, diagnostic_label,  bde, cpe_rates )
     ! Performs calculations for specified photolysis and dose rates for a
     ! given set of conditions
 
@@ -207,6 +207,10 @@ subroutine run( this, solar_zenith_angle, earth_sun_distance,               &
     real(dk),         optional, intent(out)   :: dose_rates(:,:)                ! (vertical level, reaction) [s-1]
     character(len=*), optional, intent(in)    :: diagnostic_label               ! label used in diagnostic file names
 
+    !> for photolysis heating rates
+    real(dk), optional,         intent(in)    :: bde(:,:) ! bond disociation energies (joules)
+    real(dk), optional,         intent(inout) :: cpe_rates(:,:) ! chem potential energy rates (joules/sec)
+
     ! Local variables
     character(len=*), parameter         :: Iam = 'Photolysis core run: '
     character(len=2)                    :: number
@@ -245,7 +249,7 @@ subroutine run( this, solar_zenith_angle, earth_sun_distance,               &
                                        this%profile_warehouse_,               &
                                        this%radiation_field_,                 &
                                        photolysis_rate_constants,             &
-                                       diag_label )
+                                       diag_label, bde=bde, cpe_rates=cpe_rates )
     end if
     if( associated( this%dose_rates_ ) .and. present( dose_rates ) ) then
       call this%dose_rates_%get( this%grid_warehouse_,                        &

diff --git a/src/photolysis_rates.F90 b/src/photolysis_rates.F90
@@ -264,7 +264,7 @@ end subroutine add
 
   !> calculate photolysis rate constants
   subroutine get( this, la_srb, spherical_geometry, grid_warehouse,           &
-      profile_warehouse, radiation_field, photolysis_rates, file_tag )
+      profile_warehouse, radiation_field, photolysis_rates, file_tag, bde, cpe_rates )
 
     use musica_assert,                 only : assert_msg, die_msg
     use tuvx_diagnostic_util,          only : diagout
@@ -291,6 +291,10 @@ subroutine get( this, la_srb, spherical_geometry, grid_warehouse,           &
     !> Calculated photolysis rate constants
     real(dk),                   intent(inout) :: photolysis_rates(:,:)
 
+    !> for photolysis heating rates
+    real(dk), optional,         intent(in)    :: bde(:,:) ! bond disociation energies (joules)
+    real(dk), optional,         intent(inout) :: cpe_rates(:,:) ! chem potential energy rates (joules/sec)
+
     !> Local variables
     character(len=*), parameter :: Iam = "photolysis rates calculator"
     integer               :: vertNdx, rateNdx, nRates
@@ -373,6 +377,10 @@ subroutine get( this, la_srb, spherical_geometry, grid_warehouse,           &
       do vertNdx = 1, zGrid%ncells_ + 1
         photolysis_rates( vertNdx, rateNdx ) =                                &
             dot_product( actinicFlux( :, vertNdx ), xsqy( :, vertNdx ) )
+        if (present(bde).and.present(cpe_rates)) then
+           cpe_rates( vertNdx, rateNdx ) =                                &
+                dot_product( actinicFlux(:,vertNdx), bde(:,rateNdx)*xsqy(:,vertNdx)  )
+        end if
       enddo
       if( allocated( cross_section ) ) deallocate( cross_section )
       if( allocated( quantum_yield ) ) deallocate( quantum_yield )