diff --git a/src/time_integration_mod.F90 b/src/time_integration_mod.F90
index 12b829109a1306f72778851be780902fe3bd667a..02ba5bf98e36cb10854da9f64256cd62f17cd685 100644
--- a/src/time_integration_mod.F90
+++ b/src/time_integration_mod.F90
@@ -60,56 +60,104 @@ CONTAINS
IMPLICIT NONE
SELECT CASE (numerical_scheme)
- CASE ('RK4')
- CALL RK4
- CASE ('SSPx_RK3')
- CALL SSPx_RK3
+ ! Order 2 methods
+ CASE ('RK2')
+ CALL RK2
CASE ('SSPx_RK2')
CALL SSPx_RK2
+ ! Order 3 methods
+ CASE ('RK3')
+ CALL RK3
+ CASE ('SSP_RK3')
+ CALL SSP_RK3
+ CASE ('SSPx_RK3')
+ CALL SSPx_RK3
+ CASE ('IMEX_SSP2')
+ CALL IMEX_SSP2
+ CASE ('ARK2')
+ CALL ARK2
+ ! Order 4 methods
+ CASE ('RK4')
+ CALL RK4
+ ! Order 5 methods
CASE ('DOPRI5')
CALL DOPRI5
CASE DEFAULT
IF (my_id .EQ. 0) WRITE(*,*) 'Cannot initialize time integration scheme. Name invalid.'
END SELECT
-
IF (my_id .EQ. 0) WRITE(*,*) " Time integration with ", numerical_scheme
END SUBROUTINE set_numerical_scheme
- SUBROUTINE RK4
- ! Butcher coeff for RK4 (default)
-
+ !!! second order time schemes
+ SUBROUTINE RK2
+ ! Butcher coeff for clasical RK2 (Heun's)
USE basic
USE prec_const
IMPLICIT NONE
- INTEGER,PARAMETER :: nbstep = 4
-
+ INTEGER,PARAMETER :: nbstep = 2
CALL allocate_array_dp1(c_E,1,nbstep)
CALL allocate_array_dp1(b_E,1,nbstep)
CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
-
- ntimelevel = 4
-
+ ntimelevel = 2
c_E(1) = 0.0_dp
- c_E(2) = 1.0_dp/2.0_dp
- c_E(3) = 1.0_dp/2.0_dp
- c_E(4) = 1.0_dp
+ c_E(2) = 1.0_dp
+ b_E(1) = 1._dp/2._dp
+ b_E(2) = 1._dp/2._dp
+ A_E(2,1) = 1._dp
+ END SUBROUTINE RK2
- b_E(1) = 1.0_dp/6.0_dp
- b_E(2) = 1.0_dp/3.0_dp
- b_E(3) = 1.0_dp/3.0_dp
- b_E(4) = 1.0_dp/6.0_dp
+ SUBROUTINE SSPx_RK2
+ ! DOESNT WORK
+ ! Butcher coeff for modified strong stability preserving RK2
+ ! used in GX (Hammett 2022, Mandell et al. 2022)
+ USE basic
+ USE prec_const
+ IMPLICIT NONE
+ INTEGER,PARAMETER :: nbstep = 2
+ REAL(dp) :: alpha, beta
+ alpha = 1._dp/SQRT(2._dp)
+ beta = SQRT(2._dp) - 1._dp
+ CALL allocate_array_dp1(c_E,1,nbstep)
+ CALL allocate_array_dp1(b_E,1,nbstep)
+ CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
+ ntimelevel = 2
+ c_E(1) = 0.0_dp
+ c_E(2) = 1.0_dp/2.0_dp
+ b_E(1) = alpha*beta/2._dp
+ b_E(2) = alpha/2._dp
+ A_E(2,1) = alpha
+ ! b_E(1) = 1._dp
+ ! b_E(2) = 1._dp/SQRT(2._dp)
+ ! A_E(2,1) = 1._dp/SQRT(2._dp)
+ END SUBROUTINE SSPx_RK2
+ !!! third order time schemes
+ SUBROUTINE RK3
+ ! Butcher coeff for classical RK3
+ USE basic
+ USE prec_const
+ IMPLICIT NONE
+ INTEGER,PARAMETER :: nbstep = 3
+ CALL allocate_array_dp1(c_E,1,nbstep)
+ CALL allocate_array_dp1(b_E,1,nbstep)
+ CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
+ ntimelevel = 3
+ c_E(1) = 0.0_dp
+ c_E(2) = 1.0_dp/2.0_dp
+ c_E(3) = 1.0_dp
+ b_E(1) = 1._dp/6._dp
+ b_E(2) = 2._dp/3._dp
+ b_E(3) = 1._dp/6._dp
A_E(2,1) = 1.0_dp/2.0_dp
- A_E(3,2) = 1.0_dp/2.0_dp
- A_E(4,3) = 1.0_dp
-
- END SUBROUTINE RK4
+ A_E(3,1) = -1._dp
+ A_E(3,2) = 2._dp
+ END SUBROUTINE RK3
SUBROUTINE SSPx_RK3
+ ! DOESNT WORK
! Butcher coeff for modified strong stability preserving RK3
! used in GX (Hammett 2022, Mandell et al. 2022)
-
USE basic
USE prec_const
IMPLICIT NONE
@@ -125,125 +173,122 @@ CONTAINS
! w2 = 0.3769892220587804931852815570891834795475_dp ! (6^(2/3)-1-sqrt(9-2*6^(2/3)))/2
w3 = 1._dp/a1 - w2 * (1._dp + w1)
! w3 = 1.3368459739528868457369981115334667265415_dp
-
CALL allocate_array_dp1(c_E,1,nbstep)
CALL allocate_array_dp1(b_E,1,nbstep)
CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
-
ntimelevel = 3
-
c_E(1) = 0.0_dp
c_E(2) = 1.0_dp/2.0_dp
c_E(3) = 1.0_dp/2.0_dp
-
b_E(1) = a1 * (w1*w2 + w3)
b_E(2) = a2 * w2
b_E(3) = a3
-
A_E(2,1) = a1
A_E(3,1) = a1 * w1
A_E(3,2) = a2
-
END SUBROUTINE SSPx_RK3
- SUBROUTINE SSPx_RK2
- ! Butcher coeff for modified strong stability preserving RK2
- ! used in GX (Hammett 2022, Mandell et al. 2022)
-
+ SUBROUTINE IMEX_SSP2
+ !! Version of Rokhzadi 2017 (An Optimally Stable and Accurate Second-Order
+ ! SSP Runge-Kutta IMEX Scheme for Atmospheric Applications)
USE basic
USE prec_const
IMPLICIT NONE
- INTEGER,PARAMETER :: nbstep = 2
- REAL(dp) :: alpha, beta
- alpha = 1._dp/SQRT(2._dp)
- beta = SQRT(2._dp) - 1._dp
-
+ INTEGER,PARAMETER :: nbstep = 3
CALL allocate_array_dp1(c_E,1,nbstep)
CALL allocate_array_dp1(b_E,1,nbstep)
CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
-
- ntimelevel = 2
-
- c_E(1) = 0.0_dp
- c_E(2) = 1.0_dp/2.0_dp
-
- b_E(1) = alpha*beta/2._dp
- b_E(2) = alpha/2._dp
-
- A_E(2,1) = alpha
-
- END SUBROUTINE SSPx_RK2
-
- SUBROUTINE DOPRI5
- ! Butcher coeff for DOPRI5 --> Stiffness detection
- ! DOPRI5 used for stiffness detection.
- ! 5 order method/7 stages
-
+ ntimelevel = 3
+ c_E(1) = 0._dp
+ c_E(2) = 0.711664700366941_dp
+ c_E(3) = 0.994611536833690_dp
+ b_E(1) = 0.398930808264688_dp
+ b_E(2) = 0.345755244189623_dp
+ b_E(3) = 0.255313947545689_dp
+ A_E(2,1) = 0.711664700366941_dp
+ A_E(3,1) = 0.077338168947683_dp
+ A_E(3,2) = 0.917273367886007_dp
+ END SUBROUTINE IMEX_SSP2
+
+ SUBROUTINE ARK2
+ !! Version of Rokhzadi 2017 (An Optimally Stable and Accurate Second-Order
+ ! SSP Runge-Kutta IMEX Scheme for Atmospheric Applications)
USE basic
+ USE prec_const
IMPLICIT NONE
- INTEGER,PARAMETER :: nbstep =7
-
+ INTEGER,PARAMETER :: nbstep = 3
CALL allocate_array_dp1(c_E,1,nbstep)
CALL allocate_array_dp1(b_E,1,nbstep)
CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
+ ntimelevel = 3
+ c_E(1) = 0._dp
+ c_E(2) = 2._dp*(1._dp - 1._dp/SQRT2)
+ c_E(3) = 1._dp
+ b_E(1) = 1._dp/(2._dp*SQRT2)
+ b_E(2) = 1._dp/(2._dp*SQRT2)
+ b_E(3) = 1._dp - 1._dp/SQRT2
+ A_E(2,1) = 2._dp*(1._dp - 1._dp/SQRT2)
+ A_E(3,1) = 1._dp - (3._dp + 2._dp*SQRT2)/6._dp
+ A_E(3,2) = (3._dp + 2._dp*SQRT2)/6._dp
+ END SUBROUTINE ARK2
+
+ SUBROUTINE SSP_RK3
+ ! Butcher coeff for strong stability preserving RK3
+ USE basic
+ USE prec_const
+ IMPLICIT NONE
+ INTEGER,PARAMETER :: nbstep = 3
+ CALL allocate_array_dp1(c_E,1,nbstep)
+ CALL allocate_array_dp1(b_E,1,nbstep)
+ CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
+ ntimelevel = 3
+ c_E(1) = 0.0_dp
+ c_E(2) = 1.0_dp
+ c_E(3) = 1.0_dp/2.0_dp
+ b_E(1) = 1._dp/6._dp
+ b_E(2) = 1._dp/6._dp
+ b_E(3) = 2._dp/3._dp
+ A_E(2,1) = 1._dp
+ A_E(3,1) = 1._dp/4._dp
+ A_E(3,2) = 1._dp/4._dp
+ END SUBROUTINE SSP_RK3
+
+ !!! fourth order time schemes
+ SUBROUTINE RK4
+ ! Butcher coeff for RK4 (default)
+ USE basic
+ USE prec_const
+ IMPLICIT NONE
+ INTEGER,PARAMETER :: nbstep = 4
+ CALL allocate_array_dp1(c_E,1,nbstep)
+ CALL allocate_array_dp1(b_E,1,nbstep)
+ CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
+ ntimelevel = 4
+ c_E(1) = 0.0_dp
+ c_E(2) = 1.0_dp/2.0_dp
+ c_E(3) = 1.0_dp/2.0_dp
+ c_E(4) = 1.0_dp
+ b_E(1) = 1.0_dp/6.0_dp
+ b_E(2) = 1.0_dp/3.0_dp
+ b_E(3) = 1.0_dp/3.0_dp
+ b_E(4) = 1.0_dp/6.0_dp
+ A_E(2,1) = 1.0_dp/2.0_dp
+ A_E(3,2) = 1.0_dp/2.0_dp
+ A_E(4,3) = 1.0_dp
+ END SUBROUTINE RK4
- ntimelevel = 7
- !
- c_E(1) = 0._dp
- c_E(2) = 1.0_dp/5.0_dp
- c_E(3) = 3.0_dp /10.0_dp
- c_E(4) = 4.0_dp/5.0_dp
- c_E(5) = 8.0_dp/9.0_dp
- c_E(6) = 1.0_dp
- c_E(7) = 1.0_dp
- !
- A_E(2,1) = 1.0_dp/5.0_dp
- A_E(3,1) = 3.0_dp/40.0_dp
- A_E(3,2) = 9.0_dp/40.0_dp
- A_E(4,1) = 44.0_dp/45.0_dp
- A_E(4,2) = -56.0_dp/15.0_dp
- A_E(4,3) = 32.0_dp/9.0_dp
- A_E(5,1 ) = 19372.0_dp/6561.0_dp
- A_E(5,2) = -25360.0_dp/2187.0_dp
- A_E(5,3) = 64448.0_dp/6561.0_dp
- A_E(5,4) = -212.0_dp/729.0_dp
- A_E(6,1) = 9017.0_dp/3168.0_dp
- A_E(6,2)= -355.0_dp/33.0_dp
- A_E(6,3) = 46732.0_dp/5247.0_dp
- A_E(6,4) = 49.0_dp/176.0_dp
- A_E(6,5) = -5103.0_dp/18656.0_dp
- A_E(7,1) = 35.0_dp/384.0_dp
- A_E(7,3) = 500.0_dp/1113.0_dp
- A_E(7,4) = 125.0_dp/192.0_dp
- A_E(7,5) = -2187.0_dp/6784.0_dp
- A_E(7,6) = 11.0_dp/84.0_dp
- !
- b_E(1) = 35.0_dp/384.0_dp
- b_E(2) = 0._dp
- b_E(3) = 500.0_dp/1113.0_dp
- b_E(4) = 125.0_dp/192.0_dp
- b_E(5) = -2187.0_dp/6784.0_dp
- b_E(6) = 11.0_dp/84.0_dp
- b_E(7) = 0._dp
- !
- END SUBROUTINE DOPRI5
-
- SUBROUTINE DOPRI5_ADAPT
+ !!! fifth order time schemes
+ SUBROUTINE DOPRI5
! Butcher coeff for DOPRI5 --> Stiffness detection
! DOPRI5 used for stiffness detection.
! 5 order method/7 stages
-
USE basic
IMPLICIT NONE
INTEGER,PARAMETER :: nbstep =7
-
CALL allocate_array_dp1(c_E,1,nbstep)
CALL allocate_array_dp1(b_E,1,nbstep)
- CALL allocate_array_dp1(b_Es,1,nbstep)
CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
-
ntimelevel = 7
- !
c_E(1) = 0._dp
c_E(2) = 1.0_dp/5.0_dp
c_E(3) = 3.0_dp /10.0_dp
@@ -251,7 +296,6 @@ CONTAINS
c_E(5) = 8.0_dp/9.0_dp
c_E(6) = 1.0_dp
c_E(7) = 1.0_dp
- !
A_E(2,1) = 1.0_dp/5.0_dp
A_E(3,1) = 3.0_dp/40.0_dp
A_E(3,2) = 9.0_dp/40.0_dp
@@ -272,7 +316,6 @@ CONTAINS
A_E(7,4) = 125.0_dp/192.0_dp
A_E(7,5) = -2187.0_dp/6784.0_dp
A_E(7,6) = 11.0_dp/84.0_dp
- !
b_E(1) = 35.0_dp/384.0_dp
b_E(2) = 0._dp
b_E(3) = 500.0_dp/1113.0_dp
@@ -280,15 +323,6 @@ CONTAINS
b_E(5) = -2187.0_dp/6784.0_dp
b_E(6) = 11.0_dp/84.0_dp
b_E(7) = 0._dp
- !
- b_Es(1) = 5179.0_dp/57600.0_dp
- b_Es(2) = 0._dp
- b_Es(3) = 7571.0_dp/16695.0_dp
- b_Es(4) = 393.0_dp/640.0_dp
- b_Es(5) = -92097.0_dp/339200.0_dp
- b_Es(6) = 187.0_dp/2100.0_dp
- b_Es(7) = 1._dp/40._dp
- !
- END SUBROUTINE DOPRI5_ADAPT
+ END SUBROUTINE DOPRI5
END MODULE time_integration