From 654acf8a87108e00314b90846b49dbed546e9f41 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Mon, 17 Oct 2022 07:17:19 -0700 Subject: [PATCH 1/3] use constant number of grid cells in performance driver --- test/performance/driver.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/test/performance/driver.F90 b/test/performance/driver.F90 index 38d7a9ad5..a0194b609 100644 --- a/test/performance/driver.F90 +++ b/test/performance/driver.F90 @@ -21,6 +21,7 @@ program performance_test !> Runs MICM under prescribed conditions and collect performance metrics subroutine run_test( ) + use constants, only : kNumberOfGridCells use micm_environment, only : environment_t use micm_kinetics, only : kinetics_t use micm_ODE_solver, only : ODE_solver_t From 637bcd3d52f17c41a22dce6f2d06e969a8f991e5 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Mon, 17 Oct 2022 07:40:17 -0700 Subject: [PATCH 2/3] use updated preprocessor output --- .../factor_solve_utilities.F90 | 58 ++++++------ .../rate_constants_utility.F90 | 91 ++++++++++--------- .../factor_solve_utilities.F90 | 58 ++++++------ .../rate_constants_utility.F90 | 91 ++++++++++--------- 4 files changed, 152 insertions(+), 146 deletions(-) diff --git a/src/preprocessor_output/factor_solve_utilities.F90 b/src/preprocessor_output/factor_solve_utilities.F90 index 8e82651ff..702902508 100644 --- a/src/preprocessor_output/factor_solve_utilities.F90 +++ b/src/preprocessor_output/factor_solve_utilities.F90 @@ -94,38 +94,40 @@ subroutine factor(LU) real(r8), intent(inout) :: LU(:,:) - integer :: i - - do i = 1, ncell - LU(i,1) = 1./LU(i,1) - LU(i,2) = LU(i,2) * LU(i,1) - LU(i,3) = LU(i,3) * LU(i,1) - LU(i,4) = LU(i,4) * LU(i,1) - LU(i,5) = 1./LU(i,5) - LU(i,6) = 1./LU(i,6) - LU(i,7) = 1./LU(i,7) - LU(i,8) = 1./LU(i,8) - LU(i,9) = LU(i,9) * LU(i,8) - LU(i,10) = LU(i,10) * LU(i,8) - LU(i,11) = 1./LU(i,11) - LU(i,12) = LU(i,12) * LU(i,11) - LU(i,17) = LU(i,17) - LU(i,12)*LU(i,16) - LU(i,21) = LU(i,21) - LU(i,12)*LU(i,20) - LU(i,13) = 1./LU(i,13) - LU(i,14) = LU(i,14) * LU(i,13) - LU(i,15) = LU(i,15) * LU(i,13) - LU(i,18) = LU(i,18) - LU(i,14)*LU(i,17) - LU(i,19) = LU(i,19) - LU(i,15)*LU(i,17) - LU(i,22) = LU(i,22) - LU(i,14)*LU(i,21) - LU(i,23) = LU(i,23) - LU(i,15)*LU(i,21) - LU(i,18) = 1./LU(i,18) - LU(i,19) = LU(i,19) * LU(i,18) - LU(i,23) = LU(i,23) - LU(i,19)*LU(i,22) - LU(i,23) = 1./LU(i,23) +integer :: i + +do i = 1, ncell + LU(i,1) = 1./LU(i,1) + LU(i,2) = LU(i,2) * LU(i,1) + LU(i,3) = LU(i,3) * LU(i,1) + LU(i,4) = LU(i,4) * LU(i,1) + LU(i,5) = 1./LU(i,5) + LU(i,6) = 1./LU(i,6) + LU(i,7) = 1./LU(i,7) + LU(i,8) = 1./LU(i,8) + LU(i,9) = LU(i,9) * LU(i,8) + LU(i,10) = LU(i,10) * LU(i,8) + LU(i,11) = 1./LU(i,11) + LU(i,12) = LU(i,12) * LU(i,11) + LU(i,17) = LU(i,17) - LU(i,12)*LU(i,16) + LU(i,21) = LU(i,21) - LU(i,12)*LU(i,20) + LU(i,13) = 1./LU(i,13) + LU(i,14) = LU(i,14) * LU(i,13) + LU(i,15) = LU(i,15) * LU(i,13) + LU(i,18) = LU(i,18) - LU(i,14)*LU(i,17) + LU(i,19) = LU(i,19) - LU(i,15)*LU(i,17) + LU(i,22) = LU(i,22) - LU(i,14)*LU(i,21) + LU(i,23) = LU(i,23) - LU(i,15)*LU(i,21) + LU(i,18) = 1./LU(i,18) + LU(i,19) = LU(i,19) * LU(i,18) + LU(i,23) = LU(i,23) - LU(i,19)*LU(i,22) + LU(i,23) = 1./LU(i,23) end do end subroutine factor + + subroutine solve(LU,x,b) real(r8), intent(in) :: LU(:,:), b(:,:) ! solve LU * x = b diff --git a/src/preprocessor_output/rate_constants_utility.F90 b/src/preprocessor_output/rate_constants_utility.F90 index bf82d4855..89abb0c39 100644 --- a/src/preprocessor_output/rate_constants_utility.F90 +++ b/src/preprocessor_output/rate_constants_utility.F90 @@ -57,51 +57,52 @@ subroutine calculate_rate_constants( rate_constants, environment ) integer :: i do i = 1, ncell - !O2_1 - !k_O2_1: O2 -> 2*O - photolysis = rate_constant_photolysis_t( & - photolysis_rate_constant_index = 1 ) - rate_constants( i, 1 ) = photolysis%calculate( environment(i) ) - - !O3_1 - !k_O3_1: O3 -> 1*O1D + 1*O2 - photolysis = rate_constant_photolysis_t( & - photolysis_rate_constant_index = 2 ) - rate_constants( i, 2 ) = photolysis%calculate( environment(i) ) - - !O3_2 - !k_O3_2: O3 -> 1*O + 1*O2 - photolysis = rate_constant_photolysis_t( & - photolysis_rate_constant_index = 3 ) - rate_constants( i, 3 ) = photolysis%calculate( environment(i) ) - - !N2_O1D_1 - !k_N2_O1D_1: N2 + O1D -> 1*O + 1*N2 - arrhenius = rate_constant_arrhenius_t( & - A = real( 2.15e-11, kind=musica_dk ), & - C = real( 110, kind=musica_dk ) ) - rate_constants( i, 4 ) = arrhenius%calculate( environment(i) ) - - !O1D_O2_1 - !k_O1D_O2_1: O1D + O2 -> 1*O + 1*O2 - arrhenius = rate_constant_arrhenius_t( & - A = real( 3.3e-11, kind=musica_dk ), & - C = real( 55, kind=musica_dk ) ) - rate_constants( i, 5 ) = arrhenius%calculate( environment(i) ) - - !O_O3_1 - !k_O_O3_1: O + O3 -> 2*O2 - arrhenius = rate_constant_arrhenius_t( & - A = real( 8e-12, kind=musica_dk ), & - C = real( -2060, kind=musica_dk ) ) - rate_constants( i, 6 ) = arrhenius%calculate( environment(i) ) - - !M_O_O2_1 - !k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M - arrhenius = rate_constant_arrhenius_t( & - A = real( 6e-34, kind=musica_dk ), & - B = real( 2.4, kind=musica_dk ) ) - rate_constants( i, 7 ) = arrhenius%calculate( environment(i) ) + !O2_1 + !k_O2_1: O2 -> 2*O + photolysis = rate_constant_photolysis_t( & + photolysis_rate_constant_index = 1 ) + rate_constants( i, 1 ) = photolysis%calculate( environment( i ) ) + + !O3_1 + !k_O3_1: O3 -> 1*O1D + 1*O2 + photolysis = rate_constant_photolysis_t( & + photolysis_rate_constant_index = 2 ) + rate_constants( i, 2 ) = photolysis%calculate( environment( i ) ) + + !O3_2 + !k_O3_2: O3 -> 1*O + 1*O2 + photolysis = rate_constant_photolysis_t( & + photolysis_rate_constant_index = 3 ) + rate_constants( i, 3 ) = photolysis%calculate( environment( i ) ) + + !N2_O1D_1 + !k_N2_O1D_1: N2 + O1D -> 1*O + 1*N2 + arrhenius = rate_constant_arrhenius_t( & + A = real( 2.15e-11, kind=musica_dk ), & + C = real( 110, kind=musica_dk ) ) + rate_constants( i, 4 ) = arrhenius%calculate( environment( i ) ) + + !O1D_O2_1 + !k_O1D_O2_1: O1D + O2 -> 1*O + 1*O2 + arrhenius = rate_constant_arrhenius_t( & + A = real( 3.3e-11, kind=musica_dk ), & + C = real( 55, kind=musica_dk ) ) + rate_constants( i, 5 ) = arrhenius%calculate( environment( i ) ) + + !O_O3_1 + !k_O_O3_1: O + O3 -> 2*O2 + arrhenius = rate_constant_arrhenius_t( & + A = real( 8e-12, kind=musica_dk ), & + C = real( -2060, kind=musica_dk ) ) + rate_constants( i, 6 ) = arrhenius%calculate( environment( i ) ) + + !M_O_O2_1 + !k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M + arrhenius = rate_constant_arrhenius_t( & + A = real( 6e-34, kind=musica_dk ), & + B = real( 2.4, kind=musica_dk ) ) + rate_constants( i, 7 ) = arrhenius%calculate( environment( i ) ) + end do end subroutine calculate_rate_constants diff --git a/test/preprocessor_output/factor_solve_utilities.F90 b/test/preprocessor_output/factor_solve_utilities.F90 index 8e82651ff..702902508 100644 --- a/test/preprocessor_output/factor_solve_utilities.F90 +++ b/test/preprocessor_output/factor_solve_utilities.F90 @@ -94,38 +94,40 @@ subroutine factor(LU) real(r8), intent(inout) :: LU(:,:) - integer :: i - - do i = 1, ncell - LU(i,1) = 1./LU(i,1) - LU(i,2) = LU(i,2) * LU(i,1) - LU(i,3) = LU(i,3) * LU(i,1) - LU(i,4) = LU(i,4) * LU(i,1) - LU(i,5) = 1./LU(i,5) - LU(i,6) = 1./LU(i,6) - LU(i,7) = 1./LU(i,7) - LU(i,8) = 1./LU(i,8) - LU(i,9) = LU(i,9) * LU(i,8) - LU(i,10) = LU(i,10) * LU(i,8) - LU(i,11) = 1./LU(i,11) - LU(i,12) = LU(i,12) * LU(i,11) - LU(i,17) = LU(i,17) - LU(i,12)*LU(i,16) - LU(i,21) = LU(i,21) - LU(i,12)*LU(i,20) - LU(i,13) = 1./LU(i,13) - LU(i,14) = LU(i,14) * LU(i,13) - LU(i,15) = LU(i,15) * LU(i,13) - LU(i,18) = LU(i,18) - LU(i,14)*LU(i,17) - LU(i,19) = LU(i,19) - LU(i,15)*LU(i,17) - LU(i,22) = LU(i,22) - LU(i,14)*LU(i,21) - LU(i,23) = LU(i,23) - LU(i,15)*LU(i,21) - LU(i,18) = 1./LU(i,18) - LU(i,19) = LU(i,19) * LU(i,18) - LU(i,23) = LU(i,23) - LU(i,19)*LU(i,22) - LU(i,23) = 1./LU(i,23) +integer :: i + +do i = 1, ncell + LU(i,1) = 1./LU(i,1) + LU(i,2) = LU(i,2) * LU(i,1) + LU(i,3) = LU(i,3) * LU(i,1) + LU(i,4) = LU(i,4) * LU(i,1) + LU(i,5) = 1./LU(i,5) + LU(i,6) = 1./LU(i,6) + LU(i,7) = 1./LU(i,7) + LU(i,8) = 1./LU(i,8) + LU(i,9) = LU(i,9) * LU(i,8) + LU(i,10) = LU(i,10) * LU(i,8) + LU(i,11) = 1./LU(i,11) + LU(i,12) = LU(i,12) * LU(i,11) + LU(i,17) = LU(i,17) - LU(i,12)*LU(i,16) + LU(i,21) = LU(i,21) - LU(i,12)*LU(i,20) + LU(i,13) = 1./LU(i,13) + LU(i,14) = LU(i,14) * LU(i,13) + LU(i,15) = LU(i,15) * LU(i,13) + LU(i,18) = LU(i,18) - LU(i,14)*LU(i,17) + LU(i,19) = LU(i,19) - LU(i,15)*LU(i,17) + LU(i,22) = LU(i,22) - LU(i,14)*LU(i,21) + LU(i,23) = LU(i,23) - LU(i,15)*LU(i,21) + LU(i,18) = 1./LU(i,18) + LU(i,19) = LU(i,19) * LU(i,18) + LU(i,23) = LU(i,23) - LU(i,19)*LU(i,22) + LU(i,23) = 1./LU(i,23) end do end subroutine factor + + subroutine solve(LU,x,b) real(r8), intent(in) :: LU(:,:), b(:,:) ! solve LU * x = b diff --git a/test/preprocessor_output/rate_constants_utility.F90 b/test/preprocessor_output/rate_constants_utility.F90 index bf82d4855..89abb0c39 100644 --- a/test/preprocessor_output/rate_constants_utility.F90 +++ b/test/preprocessor_output/rate_constants_utility.F90 @@ -57,51 +57,52 @@ subroutine calculate_rate_constants( rate_constants, environment ) integer :: i do i = 1, ncell - !O2_1 - !k_O2_1: O2 -> 2*O - photolysis = rate_constant_photolysis_t( & - photolysis_rate_constant_index = 1 ) - rate_constants( i, 1 ) = photolysis%calculate( environment(i) ) - - !O3_1 - !k_O3_1: O3 -> 1*O1D + 1*O2 - photolysis = rate_constant_photolysis_t( & - photolysis_rate_constant_index = 2 ) - rate_constants( i, 2 ) = photolysis%calculate( environment(i) ) - - !O3_2 - !k_O3_2: O3 -> 1*O + 1*O2 - photolysis = rate_constant_photolysis_t( & - photolysis_rate_constant_index = 3 ) - rate_constants( i, 3 ) = photolysis%calculate( environment(i) ) - - !N2_O1D_1 - !k_N2_O1D_1: N2 + O1D -> 1*O + 1*N2 - arrhenius = rate_constant_arrhenius_t( & - A = real( 2.15e-11, kind=musica_dk ), & - C = real( 110, kind=musica_dk ) ) - rate_constants( i, 4 ) = arrhenius%calculate( environment(i) ) - - !O1D_O2_1 - !k_O1D_O2_1: O1D + O2 -> 1*O + 1*O2 - arrhenius = rate_constant_arrhenius_t( & - A = real( 3.3e-11, kind=musica_dk ), & - C = real( 55, kind=musica_dk ) ) - rate_constants( i, 5 ) = arrhenius%calculate( environment(i) ) - - !O_O3_1 - !k_O_O3_1: O + O3 -> 2*O2 - arrhenius = rate_constant_arrhenius_t( & - A = real( 8e-12, kind=musica_dk ), & - C = real( -2060, kind=musica_dk ) ) - rate_constants( i, 6 ) = arrhenius%calculate( environment(i) ) - - !M_O_O2_1 - !k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M - arrhenius = rate_constant_arrhenius_t( & - A = real( 6e-34, kind=musica_dk ), & - B = real( 2.4, kind=musica_dk ) ) - rate_constants( i, 7 ) = arrhenius%calculate( environment(i) ) + !O2_1 + !k_O2_1: O2 -> 2*O + photolysis = rate_constant_photolysis_t( & + photolysis_rate_constant_index = 1 ) + rate_constants( i, 1 ) = photolysis%calculate( environment( i ) ) + + !O3_1 + !k_O3_1: O3 -> 1*O1D + 1*O2 + photolysis = rate_constant_photolysis_t( & + photolysis_rate_constant_index = 2 ) + rate_constants( i, 2 ) = photolysis%calculate( environment( i ) ) + + !O3_2 + !k_O3_2: O3 -> 1*O + 1*O2 + photolysis = rate_constant_photolysis_t( & + photolysis_rate_constant_index = 3 ) + rate_constants( i, 3 ) = photolysis%calculate( environment( i ) ) + + !N2_O1D_1 + !k_N2_O1D_1: N2 + O1D -> 1*O + 1*N2 + arrhenius = rate_constant_arrhenius_t( & + A = real( 2.15e-11, kind=musica_dk ), & + C = real( 110, kind=musica_dk ) ) + rate_constants( i, 4 ) = arrhenius%calculate( environment( i ) ) + + !O1D_O2_1 + !k_O1D_O2_1: O1D + O2 -> 1*O + 1*O2 + arrhenius = rate_constant_arrhenius_t( & + A = real( 3.3e-11, kind=musica_dk ), & + C = real( 55, kind=musica_dk ) ) + rate_constants( i, 5 ) = arrhenius%calculate( environment( i ) ) + + !O_O3_1 + !k_O_O3_1: O + O3 -> 2*O2 + arrhenius = rate_constant_arrhenius_t( & + A = real( 8e-12, kind=musica_dk ), & + C = real( -2060, kind=musica_dk ) ) + rate_constants( i, 6 ) = arrhenius%calculate( environment( i ) ) + + !M_O_O2_1 + !k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M + arrhenius = rate_constant_arrhenius_t( & + A = real( 6e-34, kind=musica_dk ), & + B = real( 2.4, kind=musica_dk ) ) + rate_constants( i, 7 ) = arrhenius%calculate( environment( i ) ) + end do end subroutine calculate_rate_constants From 2c43a51f2800084ff96a0dbec6f9d781d488e1bf Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Mon, 17 Oct 2022 11:58:41 -0700 Subject: [PATCH 3/3] revert preprocessor output files in src --- .../factor_solve_utilities.F90 | 118 +----- .../kinetics_utilities.F90 | 375 ++---------------- .../rate_constants_utility.F90 | 66 +-- 3 files changed, 50 insertions(+), 509 deletions(-) diff --git a/src/preprocessor_output/factor_solve_utilities.F90 b/src/preprocessor_output/factor_solve_utilities.F90 index 702902508..785707af2 100644 --- a/src/preprocessor_output/factor_solve_utilities.F90 +++ b/src/preprocessor_output/factor_solve_utilities.F90 @@ -1,17 +1,15 @@ module factor_solve_utilities -use musica_constants, only : r8 =>musica_dk -use constants, only : ncell=>kNumberOfGridCells +use musica_constants, only: r8 => musica_dk ! This code was generated by Preprocessor revision NA ! Preprocessor source NA -! This code is generated from tag undefined of the mechanism, undefined. It is named undefined -! This tag was created on undefined by undefined and is marked as not buggy +! This is a stub of the preprocessor-generated module implicit none - integer, parameter :: number_sparse_factor_elements = 23 + integer, parameter :: number_sparse_factor_elements = 0 public :: factor, solve @@ -20,32 +18,11 @@ module factor_solve_utilities subroutine backsolve_L_y_eq_b(LU,b,y) - real(r8), intent(in) :: LU(:,:) - real(r8), intent(in) :: b(:,:) - real(r8), intent(out) :: y(:,:) - - integer :: i - - do i = 1, ncell - y(i,1) = b(i,1) - y(i,2) = b(i,2) - y(i,3) = b(i,3) - y(i,4) = b(i,4) - y(i,5) = b(i,5) - y(i,6) = b(i,6) - y(i,6) = y(i,6) - LU(i,9) * y(i,5) - y(i,7) = b(i,7) - y(i,7) = y(i,7) - LU(i,2) * y(i,1) - y(i,7) = y(i,7) - LU(i,10) * y(i,5) - y(i,7) = y(i,7) - LU(i,12) * y(i,6) - y(i,8) = b(i,8) - y(i,8) = y(i,8) - LU(i,3) * y(i,1) - y(i,8) = y(i,8) - LU(i,14) * y(i,7) - y(i,9) = b(i,9) - y(i,9) = y(i,9) - LU(i,4) * y(i,1) - y(i,9) = y(i,9) - LU(i,15) * y(i,7) - y(i,9) = y(i,9) - LU(i,19) * y(i,8) - end do + + real(r8), intent(in) :: LU(:) + real(r8), intent(in) :: b(:) + real(r8), intent(out) :: y(:) + end subroutine backsolve_L_y_eq_b @@ -53,87 +30,28 @@ end subroutine backsolve_L_y_eq_b subroutine backsolve_U_x_eq_y(LU,y,x) - real(r8), intent(in) :: LU(:,:) - real(r8), intent(in) :: y(:,:) - real(r8), intent(out) :: x(:,:) + + real(r8), intent(in) :: LU(:) + real(r8), intent(in) :: y(:) + real(r8), intent(out) :: x(:) real(r8) :: temporary - integer :: i - - do i = 1, ncell - temporary = y(i,9) - x(i,9) = LU(i,23) * temporary - temporary = y(i,8) - temporary = temporary - LU(i,22) * x(i,9) - x(i,8) = LU(i,18) * temporary - temporary = y(i,7) - temporary = temporary - LU(i,17) * x(i,8) - temporary = temporary - LU(i,21) * x(i,9) - x(i,7) = LU(i,13) * temporary - temporary = y(i,6) - temporary = temporary - LU(i,16) * x(i,8) - temporary = temporary - LU(i,20) * x(i,9) - x(i,6) = LU(i,11) * temporary - temporary = y(i,5) - x(i,5) = LU(i,8) * temporary - temporary = y(i,4) - x(i,4) = LU(i,7) * temporary - temporary = y(i,3) - x(i,3) = LU(i,6) * temporary - temporary = y(i,2) - x(i,2) = LU(i,5) * temporary - temporary = y(i,1) - x(i,1) = LU(i,1) * temporary - end do end subroutine backsolve_U_x_eq_y - - subroutine factor(LU) - real(r8), intent(inout) :: LU(:,:) - -integer :: i - -do i = 1, ncell - LU(i,1) = 1./LU(i,1) - LU(i,2) = LU(i,2) * LU(i,1) - LU(i,3) = LU(i,3) * LU(i,1) - LU(i,4) = LU(i,4) * LU(i,1) - LU(i,5) = 1./LU(i,5) - LU(i,6) = 1./LU(i,6) - LU(i,7) = 1./LU(i,7) - LU(i,8) = 1./LU(i,8) - LU(i,9) = LU(i,9) * LU(i,8) - LU(i,10) = LU(i,10) * LU(i,8) - LU(i,11) = 1./LU(i,11) - LU(i,12) = LU(i,12) * LU(i,11) - LU(i,17) = LU(i,17) - LU(i,12)*LU(i,16) - LU(i,21) = LU(i,21) - LU(i,12)*LU(i,20) - LU(i,13) = 1./LU(i,13) - LU(i,14) = LU(i,14) * LU(i,13) - LU(i,15) = LU(i,15) * LU(i,13) - LU(i,18) = LU(i,18) - LU(i,14)*LU(i,17) - LU(i,19) = LU(i,19) - LU(i,15)*LU(i,17) - LU(i,22) = LU(i,22) - LU(i,14)*LU(i,21) - LU(i,23) = LU(i,23) - LU(i,15)*LU(i,21) - LU(i,18) = 1./LU(i,18) - LU(i,19) = LU(i,19) * LU(i,18) - LU(i,23) = LU(i,23) - LU(i,19)*LU(i,22) - LU(i,23) = 1./LU(i,23) - end do - -end subroutine factor + real(r8), intent(inout) :: LU(:) -subroutine solve(LU,x,b) +end subroutine factor - real(r8), intent(in) :: LU(:,:), b(:,:) ! solve LU * x = b - real(r8), intent(out) :: x(:,:) +subroutine solve(LU, x, b) - real(r8) :: y(ncell,size(b,2)) + real(r8), intent(in) :: LU(:), b(:) ! solve LU * x = b + real(r8), intent(out) :: x(:) + real(r8) :: y(size(b)) call backsolve_L_y_eq_b(LU, b, y) call backsolve_U_x_eq_y(LU, y, x) diff --git a/src/preprocessor_output/kinetics_utilities.F90 b/src/preprocessor_output/kinetics_utilities.F90 index c03b6218b..99816ee90 100644 --- a/src/preprocessor_output/kinetics_utilities.F90 +++ b/src/preprocessor_output/kinetics_utilities.F90 @@ -4,11 +4,9 @@ module kinetics_utilities ! This code was generated by Preprocessor revision NA ! Preprocessor source NA -! This code is generated from tag undefined of the mechanism, undefined. It is named undefined -! This tag was created on undefined by undefined and is marked as not buggy +! This is a stub of the preprocessor generated module - use factor_solve_utilities, only : factor - use constants, only : ncell=>kNumberOfGridCells + use factor_solve_utilities, only: factor implicit none @@ -17,164 +15,35 @@ module kinetics_utilities photolysis_names, dforce_dy, species_names ! Total number of reactions - integer, parameter, public :: number_of_reactions = 7 - integer, parameter, public :: number_of_photolysis_reactions = 3 - integer, parameter, public :: number_of_species = 9 + integer, parameter, public :: number_of_reactions = 0 + integer, parameter, public :: number_of_photolysis_reactions = 0 + integer, parameter, public :: number_of_species = 0 contains subroutine dforce_dy(LU, rate_constant, number_density, number_density_air) + ! Compute the derivative of the Forcing w.r.t. each chemical ! Also known as the Jacobian - real(r8), intent(out) :: LU(:,:) - real(r8), intent(in) :: rate_constant(:,:) - real(r8), intent(in) :: number_density(:,:) - real(r8), intent(in) :: number_density_air(:) - - integer :: i - - LU(:,:) = 0 - - do i = 1, ncell - ! df_O/d(M) - ! k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M - LU(i,2) = LU(i,2) - rate_constant(i,7) * number_density(i,7) * number_density(i,8) - - ! df_O2/d(M) - ! k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M - LU(i,3) = LU(i,3) - rate_constant(i,7) * number_density(i,7) * number_density(i,8) - - ! df_O3/d(M) - ! k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M - LU(i,4) = LU(i,4) + rate_constant(i,7) * number_density(i,7) * number_density(i,8) - - ! df_O1D/d(N2) - ! k_N2_O1D_1: N2 + O1D -> 1*O + 1*N2 - LU(i,9) = LU(i,9) - rate_constant(i,4) * number_density(i,6) - - ! df_O/d(N2) - ! k_N2_O1D_1: N2 + O1D -> 1*O + 1*N2 - LU(i,10) = LU(i,10) + rate_constant(i,4) * number_density(i,6) - - ! df_O1D/d(O1D) - ! k_N2_O1D_1: N2 + O1D -> 1*O + 1*N2 - LU(i,11) = LU(i,11) - rate_constant(i,4) * number_density(i,5) - - ! k_O1D_O2_1: O1D + O2 -> 1*O + 1*O2 - LU(i,11) = LU(i,11) - rate_constant(i,5) * number_density(i,8) - - ! df_O/d(O1D) - ! k_N2_O1D_1: N2 + O1D -> 1*O + 1*N2 - LU(i,12) = LU(i,12) + rate_constant(i,4) * number_density(i,5) - - ! k_O1D_O2_1: O1D + O2 -> 1*O + 1*O2 - LU(i,12) = LU(i,12) + rate_constant(i,5) * number_density(i,8) - - ! df_O/d(O) - ! k_O_O3_1: O + O3 -> 2*O2 - LU(i,13) = LU(i,13) - rate_constant(i,6) * number_density(i,9) - - ! k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M - LU(i,13) = LU(i,13) - rate_constant(i,7) * number_density(i,1) * number_density(i,8) - - ! df_O2/d(O) - ! k_O_O3_1: O + O3 -> 2*O2 - LU(i,14) = LU(i,14) + 2*rate_constant(i,6) * number_density(i,9) - - ! k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M - LU(i,14) = LU(i,14) - rate_constant(i,7) * number_density(i,1) * number_density(i,8) - - ! df_O3/d(O) - ! k_O_O3_1: O + O3 -> 2*O2 - LU(i,15) = LU(i,15) - rate_constant(i,6) * number_density(i,9) - - ! k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M - LU(i,15) = LU(i,15) + rate_constant(i,7) * number_density(i,1) * number_density(i,8) - - ! df_O1D/d(O2) - ! k_O1D_O2_1: O1D + O2 -> 1*O + 1*O2 - LU(i,16) = LU(i,16) - rate_constant(i,5) * number_density(i,6) - - ! df_O/d(O2) - ! k_O2_1: O2 -> 2*O - LU(i,17) = LU(i,17) + 2*rate_constant(i,1) - - ! k_O1D_O2_1: O1D + O2 -> 1*O + 1*O2 - LU(i,17) = LU(i,17) + rate_constant(i,5) * number_density(i,6) - - ! k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M - LU(i,17) = LU(i,17) - rate_constant(i,7) * number_density(i,1) * number_density(i,7) - - ! df_O2/d(O2) - ! k_O2_1: O2 -> 2*O - LU(i,18) = LU(i,18) - rate_constant(i,1) - - ! k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M - LU(i,18) = LU(i,18) - rate_constant(i,7) * number_density(i,1) * number_density(i,7) + real(r8), intent(out) :: LU(:) + real(r8), intent(in) :: rate_constant(:) + real(r8), intent(in) :: number_density(:) + real(r8), intent(in) :: number_density_air - ! df_O3/d(O2) - ! k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M - LU(i,19) = LU(i,19) + rate_constant(i,7) * number_density(i,1) * number_density(i,7) + LU(:) = 0 - ! df_O1D/d(O3) - ! k_O3_1: O3 -> 1*O1D + 1*O2 - LU(i,20) = LU(i,20) + rate_constant(i,2) - - ! df_O/d(O3) - ! k_O3_2: O3 -> 1*O + 1*O2 - LU(i,21) = LU(i,21) + rate_constant(i,3) - - ! k_O_O3_1: O + O3 -> 2*O2 - LU(i,21) = LU(i,21) - rate_constant(i,6) * number_density(i,7) - - ! df_O2/d(O3) - ! k_O3_1: O3 -> 1*O1D + 1*O2 - LU(i,22) = LU(i,22) + rate_constant(i,2) - - ! k_O3_2: O3 -> 1*O + 1*O2 - LU(i,22) = LU(i,22) + rate_constant(i,3) - - ! k_O_O3_1: O + O3 -> 2*O2 - LU(i,22) = LU(i,22) + 2*rate_constant(i,6) * number_density(i,7) - - ! df_O3/d(O3) - ! k_O3_1: O3 -> 1*O1D + 1*O2 - LU(i,23) = LU(i,23) - rate_constant(i,2) - - ! k_O3_2: O3 -> 1*O + 1*O2 - LU(i,23) = LU(i,23) - rate_constant(i,3) - - ! k_O_O3_1: O + O3 -> 2*O2 - LU(i,23) = LU(i,23) - rate_constant(i,6) * number_density(i,7) - - end do end subroutine dforce_dy subroutine factored_alpha_minus_jac(LU, alpha, dforce_dy) - ! Compute LU decomposition of [alpha * I - dforce_dy] + !compute LU decomposition of [alpha * I - dforce_dy] - real(r8), intent(in) :: dforce_dy(:,:) + real(r8), intent(in) :: dforce_dy(:) real(r8), intent(in) :: alpha - real(r8), intent(out) :: LU(:,:) + real(r8), intent(out) :: LU(:) - integer :: i - - LU(:,:) = -dforce_dy(:,:) - - ! add alpha to diagonal elements - do i = 1, ncell - LU(i,1) = -dforce_dy(i,1) + alpha - LU(i,5) = -dforce_dy(i,5) + alpha - LU(i,6) = -dforce_dy(i,6) + alpha - LU(i,7) = -dforce_dy(i,7) + alpha - LU(i,8) = -dforce_dy(i,8) + alpha - LU(i,11) = -dforce_dy(i,11) + alpha - LU(i,13) = -dforce_dy(i,13) + alpha - LU(i,18) = -dforce_dy(i,18) + alpha - LU(i,23) = -dforce_dy(i,23) + alpha - end do + LU(:) = -dforce_dy(:) call factor(LU) @@ -183,133 +52,22 @@ end subroutine factored_alpha_minus_jac subroutine p_force(rate_constant, number_density, number_density_air, force) ! Compute force function for all molecules - real(r8), intent(in) :: rate_constant(:,:) - real(r8), intent(in) :: number_density(:,:) - real(r8), intent(in) :: number_density_air(:) - real(r8), intent(out) :: force(:,:) - - integer :: i - - do i = 1, ncell - - ! M - force(i,1) = 0 - - ! Ar - force(i,2) = 0 - - ! CO2 - force(i,3) = 0 - - ! H2O - force(i,4) = 0 - - ! N2 - force(i,5) = 0 - - ! O1D - force(i,6) = 0 - - ! k_O3_1: O3 -> 1*O1D + 1*O2 - force(i,6) = force(i,6) + rate_constant(i,2) * number_density(i,9) - - ! k_N2_O1D_1: N2 + O1D -> 1*O + 1*N2 - force(i,6) = force(i,6) - rate_constant(i,4) * number_density(i,5) * number_density(i,6) - - ! k_O1D_O2_1: O1D + O2 -> 1*O + 1*O2 - force(i,6) = force(i,6) - rate_constant(i,5) * number_density(i,6) * number_density(i,8) - - ! O - force(i,7) = 0 + real(r8), intent(in) :: rate_constant(:) + real(r8), intent(in) :: number_density(:) + real(r8), intent(in) :: number_density_air + real(r8), intent(out) :: force(:) - ! k_O2_1: O2 -> 2*O - force(i,7) = force(i,7) + 2*rate_constant(i,1) * number_density(i,8) - - ! k_O3_2: O3 -> 1*O + 1*O2 - force(i,7) = force(i,7) + rate_constant(i,3) * number_density(i,9) - - ! k_N2_O1D_1: N2 + O1D -> 1*O + 1*N2 - force(i,7) = force(i,7) + rate_constant(i,4) * number_density(i,5) * number_density(i,6) - - ! k_O1D_O2_1: O1D + O2 -> 1*O + 1*O2 - force(i,7) = force(i,7) + rate_constant(i,5) * number_density(i,6) * number_density(i,8) - - ! k_O_O3_1: O + O3 -> 2*O2 - force(i,7) = force(i,7) - rate_constant(i,6) * number_density(i,7) * number_density(i,9) - - ! k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M - force(i,7) = force(i,7) - rate_constant(i,7) * number_density(i,1) * number_density(i,7) * number_density(i,8) - - ! O2 - force(i,8) = 0 - - ! k_O2_1: O2 -> 2*O - force(i,8) = force(i,8) - rate_constant(i,1) * number_density(i,8) - - ! k_O3_1: O3 -> 1*O1D + 1*O2 - force(i,8) = force(i,8) + rate_constant(i,2) * number_density(i,9) - - ! k_O3_2: O3 -> 1*O + 1*O2 - force(i,8) = force(i,8) + rate_constant(i,3) * number_density(i,9) - - ! k_O_O3_1: O + O3 -> 2*O2 - force(i,8) = force(i,8) + 2*rate_constant(i,6) * number_density(i,7) * number_density(i,9) - - ! k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M - force(i,8) = force(i,8) - rate_constant(i,7) * number_density(i,1) * number_density(i,7) * number_density(i,8) - - ! O3 - force(i,9) = 0 - - ! k_O3_1: O3 -> 1*O1D + 1*O2 - force(i,9) = force(i,9) - rate_constant(i,2) * number_density(i,9) - - ! k_O3_2: O3 -> 1*O + 1*O2 - force(i,9) = force(i,9) - rate_constant(i,3) * number_density(i,9) - - ! k_O_O3_1: O + O3 -> 2*O2 - force(i,9) = force(i,9) - rate_constant(i,6) * number_density(i,7) * number_density(i,9) - - ! k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M - force(i,9) = force(i,9) + rate_constant(i,7) * number_density(i,1) * number_density(i,7) * number_density(i,8) - end do end subroutine p_force function reaction_rates(rate_constant, number_density, number_density_air) ! Compute reaction rates - real(r8), intent(in) :: rate_constant(:,:) - real(r8), intent(in) :: number_density(:,:) - real(r8), intent(in) :: number_density_air(:) - - integer :: i - real(r8) :: reaction_rates(ncell,number_of_reactions) - - do i = 1, ncell - - ! k_O2_1: O2 -> 2*O - reaction_rates(i,1) = rate_constant(i,1) * number_density(i,8) - - ! k_O3_1: O3 -> 1*O1D + 1*O2 - reaction_rates(i,2) = rate_constant(i,2) * number_density(i,9) - - ! k_O3_2: O3 -> 1*O + 1*O2 - reaction_rates(i,3) = rate_constant(i,3) * number_density(i,9) - - ! k_N2_O1D_1: N2 + O1D -> 1*O + 1*N2 - reaction_rates(i,4) = rate_constant(i,4) * number_density(i,5) * number_density(i,6) - - ! k_O1D_O2_1: O1D + O2 -> 1*O + 1*O2 - reaction_rates(i,5) = rate_constant(i,5) * number_density(i,6) * number_density(i,8) + real(r8) :: reaction_rates(number_of_reactions) + real(r8), intent(in) :: rate_constant(:) + real(r8), intent(in) :: number_density(:) + real(r8), intent(in) :: number_density_air - ! k_O_O3_1: O + O3 -> 2*O2 - reaction_rates(i,6) = rate_constant(i,6) * number_density(i,7) * number_density(i,9) - - ! k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M - reaction_rates(i,7) = rate_constant(i,7) * number_density(i,1) * number_density(i,7) * number_density(i,8) - - end do end function reaction_rates @@ -319,13 +77,6 @@ function reaction_names() character(len=128) :: reaction_names(number_of_reactions) - reaction_names(1) = 'O2_1' - reaction_names(2) = 'O3_1' - reaction_names(3) = 'O3_2' - reaction_names(4) = 'N2_O1D_1' - reaction_names(5) = 'O1D_O2_1' - reaction_names(6) = 'O_O3_1' - reaction_names(7) = 'M_O_O2_1' end function reaction_names @@ -335,9 +86,6 @@ function photolysis_names() character(len=28) :: photolysis_names(number_of_photolysis_reactions) - photolysis_names(1) = 'O2_1' - photolysis_names(2) = 'O3_1' - photolysis_names(3) = 'O3_2' end function photolysis_names @@ -347,88 +95,21 @@ function species_names() character(len=128) :: species_names(number_of_species) - species_names(1) = 'M' - species_names(2) = 'Ar' - species_names(3) = 'CO2' - species_names(4) = 'H2O' - species_names(5) = 'N2' - species_names(6) = 'O1D' - species_names(7) = 'O' - species_names(8) = 'O2' - species_names(9) = 'O3' end function species_names pure subroutine dforce_dy_times_vector(dforce_dy, vector, cummulative_product) + ! Compute product of [ dforce_dy * vector ] ! Commonly used to compute time-truncation errors [dforce_dy * force ] - real(r8), intent(in) :: dforce_dy(:,:) ! Jacobian of forcing - real(r8), intent(in) :: vector(:,:) ! Vector ordered as the order of number density in dy - real(r8), intent(out) :: cummulative_product(:,:) ! Product of jacobian with vector - - integer :: i - - cummulative_product(:,:) = 0 - - do i = 1, ncell - - ! df_O/d(M) * M_temporary - cummulative_product(i,7) = cummulative_product(i,7) + dforce_dy(i,2) * vector(i,1) - - ! df_O2/d(M) * M_temporary - cummulative_product(i,8) = cummulative_product(i,8) + dforce_dy(i,3) * vector(i,1) - - ! df_O3/d(M) * M_temporary - cummulative_product(i,9) = cummulative_product(i,9) + dforce_dy(i,4) * vector(i,1) - - ! df_O1D/d(N2) * N2_temporary - cummulative_product(i,6) = cummulative_product(i,6) + dforce_dy(i,9) * vector(i,5) - - ! df_O/d(N2) * N2_temporary - cummulative_product(i,7) = cummulative_product(i,7) + dforce_dy(i,10) * vector(i,5) - - ! df_O1D/d(O1D) * O1D_temporary - cummulative_product(i,6) = cummulative_product(i,6) + dforce_dy(i,11) * vector(i,6) - - ! df_O/d(O1D) * O1D_temporary - cummulative_product(i,7) = cummulative_product(i,7) + dforce_dy(i,12) * vector(i,6) - - ! df_O/d(O) * O_temporary - cummulative_product(i,7) = cummulative_product(i,7) + dforce_dy(i,13) * vector(i,7) - - ! df_O2/d(O) * O_temporary - cummulative_product(i,8) = cummulative_product(i,8) + dforce_dy(i,14) * vector(i,7) - - ! df_O3/d(O) * O_temporary - cummulative_product(i,9) = cummulative_product(i,9) + dforce_dy(i,15) * vector(i,7) - - ! df_O1D/d(O2) * O2_temporary - cummulative_product(i,6) = cummulative_product(i,6) + dforce_dy(i,16) * vector(i,8) - - ! df_O/d(O2) * O2_temporary - cummulative_product(i,7) = cummulative_product(i,7) + dforce_dy(i,17) * vector(i,8) - - ! df_O2/d(O2) * O2_temporary - cummulative_product(i,8) = cummulative_product(i,8) + dforce_dy(i,18) * vector(i,8) - - ! df_O3/d(O2) * O2_temporary - cummulative_product(i,9) = cummulative_product(i,9) + dforce_dy(i,19) * vector(i,8) - - ! df_O1D/d(O3) * O3_temporary - cummulative_product(i,6) = cummulative_product(i,6) + dforce_dy(i,20) * vector(i,9) - - ! df_O/d(O3) * O3_temporary - cummulative_product(i,7) = cummulative_product(i,7) + dforce_dy(i,21) * vector(i,9) - - ! df_O2/d(O3) * O3_temporary - cummulative_product(i,8) = cummulative_product(i,8) + dforce_dy(i,22) * vector(i,9) + real(r8), intent(in) :: dforce_dy(:) ! Jacobian of forcing + real(r8), intent(in) :: vector(:) ! Vector ordered as the order of number density in dy + real(r8), intent(out) :: cummulative_product(:) ! Product of jacobian with vector - ! df_O3/d(O3) * O3_temporary - cummulative_product(i,9) = cummulative_product(i,9) + dforce_dy(i,23) * vector(i,9) + cummulative_product(:) = 0 - end do end subroutine dforce_dy_times_vector diff --git a/src/preprocessor_output/rate_constants_utility.F90 b/src/preprocessor_output/rate_constants_utility.F90 index 89abb0c39..e30974648 100644 --- a/src/preprocessor_output/rate_constants_utility.F90 +++ b/src/preprocessor_output/rate_constants_utility.F90 @@ -1,14 +1,8 @@ ! Copyright (C) 2020 National Center for Atmospheric Research ! SPDX-License-Identifier: Apache-2.0 ! -! This code was generated by Preprocessor revision NA -! Preprocessor source NA -! -! This code is generated from tag undefined of the mechanism, undefined. It is named undefined -! This tag was created on undefined by undefined and is marked as not buggy -! !> \file -!> The mechanism-specific rate constant module +!> Stub functions for the mechanism rate constants !> Mechanism-specific rate constant functions module rate_constants_utility @@ -27,7 +21,6 @@ module rate_constants_utility use micm_rate_constant_wennberg_tunneling, & only : rate_constant_wennberg_tunneling_t use musica_constants, only : musica_dk - use constants, only : ncell=>kNumberOfGridCells implicit none private @@ -42,9 +35,9 @@ module rate_constants_utility subroutine calculate_rate_constants( rate_constants, environment ) !> Rate constant for each reaction [(molec cm-3)^(n-1) s-1] - real(kind=musica_dk), intent(out) :: rate_constants(:,:) - !> Environmental states for each grid cell - type(environment_t), intent(in) :: environment(:) + real(kind=musica_dk), intent(out) :: rate_constants(:) + !> Environmental state + type(environment_t), intent(in) :: environment type( rate_constant_arrhenius_t ) :: arrhenius type( rate_constant_photolysis_t ) :: photolysis @@ -54,57 +47,6 @@ subroutine calculate_rate_constants( rate_constants, environment ) type( rate_constant_wennberg_nitrate_t ) :: wennberg_nitrate type( rate_constant_wennberg_tunneling_t ) :: wennberg_tunneling - integer :: i - - do i = 1, ncell - !O2_1 - !k_O2_1: O2 -> 2*O - photolysis = rate_constant_photolysis_t( & - photolysis_rate_constant_index = 1 ) - rate_constants( i, 1 ) = photolysis%calculate( environment( i ) ) - - !O3_1 - !k_O3_1: O3 -> 1*O1D + 1*O2 - photolysis = rate_constant_photolysis_t( & - photolysis_rate_constant_index = 2 ) - rate_constants( i, 2 ) = photolysis%calculate( environment( i ) ) - - !O3_2 - !k_O3_2: O3 -> 1*O + 1*O2 - photolysis = rate_constant_photolysis_t( & - photolysis_rate_constant_index = 3 ) - rate_constants( i, 3 ) = photolysis%calculate( environment( i ) ) - - !N2_O1D_1 - !k_N2_O1D_1: N2 + O1D -> 1*O + 1*N2 - arrhenius = rate_constant_arrhenius_t( & - A = real( 2.15e-11, kind=musica_dk ), & - C = real( 110, kind=musica_dk ) ) - rate_constants( i, 4 ) = arrhenius%calculate( environment( i ) ) - - !O1D_O2_1 - !k_O1D_O2_1: O1D + O2 -> 1*O + 1*O2 - arrhenius = rate_constant_arrhenius_t( & - A = real( 3.3e-11, kind=musica_dk ), & - C = real( 55, kind=musica_dk ) ) - rate_constants( i, 5 ) = arrhenius%calculate( environment( i ) ) - - !O_O3_1 - !k_O_O3_1: O + O3 -> 2*O2 - arrhenius = rate_constant_arrhenius_t( & - A = real( 8e-12, kind=musica_dk ), & - C = real( -2060, kind=musica_dk ) ) - rate_constants( i, 6 ) = arrhenius%calculate( environment( i ) ) - - !M_O_O2_1 - !k_M_O_O2_1: M + O + O2 -> 1*O3 + 1*M - arrhenius = rate_constant_arrhenius_t( & - A = real( 6e-34, kind=musica_dk ), & - B = real( 2.4, kind=musica_dk ) ) - rate_constants( i, 7 ) = arrhenius%calculate( environment( i ) ) - - end do - end subroutine calculate_rate_constants !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!