diff --git a/sources/evolution.F90 b/sources/evolution.F90 index c953f32..8b7fccc 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -1109,7 +1109,7 @@ module evolution ! add forcing contribution ! - if (forcing_enabled) call update_forcing(dt) + if (forcing_enabled) call update_forcing(time, dt) ! check if we need to perform the refinement step ! diff --git a/sources/forcing.F90 b/sources/forcing.F90 index 5a5858b..af248b5 100644 --- a/sources/forcing.F90 +++ b/sources/forcing.F90 @@ -35,8 +35,8 @@ module forcing ! interfaces for procedure pointers ! abstract interface - subroutine update_forcing_iface(dt) - real(kind=8), intent(in) :: dt + subroutine update_forcing_iface(t, dt) + real(kind=8), intent(in) :: t, dt end subroutine end interface @@ -89,6 +89,8 @@ module forcing real(kind=8), save :: zth = 0.0d+00 real(kind=8), save :: zfc = 1.0d+00 #endif /* NDIMS == 3 */ + real(kind=8) :: tinj_start = 0.00000d+00 + real(kind=8) :: tinj_stop = 9.99999d+99 ! Fourier profile parameters ! @@ -244,6 +246,21 @@ module forcing end if #endif /* NDIMS == 3 */ +! get the injection period +! + call get_parameter("injection_time_start", tinj_start) + call get_parameter("injection_time_stop" , tinj_stop ) + if (tinj_start >= tinj_stop) then + if (verbose) then + write(*,*) + write(*,wfmt) "Injection cannot stop before it starts!" + write(*,wfmt) "Please, modify parameter 'injection_time_start'" // & + " or 'injection_time_stop'." + write(*,wfmt) "Resetting 'injection_time_start' to 0.0!" + tinj_start = 0.0d+00 + end if + end if + ! select the energy injection method ! select case(trim(injection)) @@ -929,6 +946,9 @@ module forcing call print_parameter(verbose, "amplitude threshold", fmin) end if + call print_parameter(verbose, "injection starts at", tinj_start) + call print_parameter(verbose, "injection stops at", tinj_stop) + if (.not. periodic(1)) & call print_parameter(verbose, "profile X transition between", xlo, xup) if (.not. periodic(2)) & @@ -953,18 +973,19 @@ module forcing ! ! Arguments: ! -! dt - time increment; +! t - time; +! dt - time step; ! !=============================================================================== ! - subroutine update_forcing_eddies(dt) + subroutine update_forcing_eddies(t, dt) use coordinates, only : xmin, ymin, zmin, xlen, ylen, zlen use random , only : randuni, randsym implicit none - real(kind=8), intent(in) :: dt + real(kind=8), intent(in) :: t, dt integer :: ni, n #if NDIMS == 3 @@ -974,6 +995,8 @@ module forcing !------------------------------------------------------------------------------- ! + if (t < tinj_start .or. t > tinj_stop) return + ! calculate the number of eddies to be injected during this interval ! dinj = dinj + rate * dt @@ -1029,11 +1052,12 @@ module forcing ! ! Arguments: ! -! dt - time increment; +! t - time; +! dt - time step; ! !=============================================================================== ! - subroutine update_forcing_oh(dt) + subroutine update_forcing_oh(t, dt) #if NDIMS == 3 use constants, only : pi2 @@ -1042,7 +1066,7 @@ module forcing implicit none - real(kind=8), intent(in) :: dt + real(kind=8), intent(in) :: t, dt integer :: l real(kind=8) :: acoeff, dcoeff @@ -1055,6 +1079,8 @@ module forcing !------------------------------------------------------------------------------- ! + if (t < tinj_start .or. t > tinj_stop) return + ! calculate drift and diffusion coefficients ! acoeff = exp(- dt / tscale) @@ -1113,18 +1139,19 @@ module forcing ! ! Arguments: ! -! dt - time increment; +! t - time; +! dt - time step; ! !=============================================================================== ! - subroutine update_forcing_alvelius(dt) + subroutine update_forcing_alvelius(t, dt) use constants, only : pi2 use random , only : randuni implicit none - real(kind=8), intent(in) :: dt + real(kind=8), intent(in) :: t, dt integer :: l real(kind=8) :: th1, dinj, sqdt @@ -1139,6 +1166,8 @@ module forcing !------------------------------------------------------------------------------- ! + if (t < tinj_start .or. t > tinj_stop) return + ! determine velocify coefficients in Fourier space ! call get_vcoefs() diff --git a/sources/interpolations.F90 b/sources/interpolations.F90 index ab77d7a..4393c46 100644 --- a/sources/interpolations.F90 +++ b/sources/interpolations.F90 @@ -4283,7 +4283,7 @@ module interpolations ! The method uses an implicit tri-diagonal solver and its coefficients are ! optimized for the maximum value of imaginary part of the scaled modified ! wave number (Im(κ)) to be 2e-6 and the integrated error of compact -! discretization E to be 6.4e-12. +! discretization E to be 6.017e-12 (at local minimum of E/r). ! ! Arguments are described in subroutine reconstruct(). ! @@ -4316,8 +4316,8 @@ module interpolations real(kind=8), dimension(size(fc)) :: r real(kind=8), dimension(size(fc)) :: u - real(kind=8), parameter :: a1 = 5.0148583220385203994203557d-01 - real(kind=8), parameter :: a2 = 2.5452055302783787784904423d-01 + real(kind=8), parameter :: a1 = 5.01241562476900742404893381352899d-01 + real(kind=8), parameter :: a2 = 2.54639102475156530564823003088612d-01 real(kind=8), dimension(3), parameter :: di5 = [ a1, 1.0d+00, a2 ] real(kind=8), dimension(5), parameter :: & ci5 = [- 3.0d+00 * a1 - 3.0d+00 * a2 + 2.0d+00, & @@ -4415,7 +4415,7 @@ module interpolations ! The method uses an implicit tri-diagonal solver and its coefficients are ! optimized for the maximum value of imaginary part of the scaled modified ! wave number (Im(κ)) to be 2e-6 and the integrated error of compact -! discretization E to be 1.0e-11. +! discretization E to be 9.4859e-12 (at local minimum of E/r). ! ! Arguments are described in subroutine reconstruct(). ! @@ -4439,8 +4439,8 @@ module interpolations real(kind=8), dimension(size(fc)) :: r real(kind=8), dimension(size(fc)) :: u - real(kind=8), parameter :: a1 = 5.0119544102147968726471782d-01 - real(kind=8), parameter :: a2 = 3.0649545768561641522180678d-01 + real(kind=8), parameter :: a1 = 5.00765434358710549355767793738926d-01 + real(kind=8), parameter :: a2 = 3.06739643875254943361647028328465d-01 real(kind=8), dimension(3), parameter :: di7 = [ a1, 1.0d+00, a2 ] real(kind=8), dimension(7), parameter :: & ci7 = [ 4.00d+00 * a1 + 4.00d+00 * a2 - 3.00d+00, & @@ -4544,7 +4544,7 @@ module interpolations ! The method uses an implicit tri-diagonal solver and its coefficients are ! optimized for the maximum value of imaginary part of the scaled modified ! wave number (Im(κ)) to be 2e-6 and the integrated error of compact -! discretization E to be 1.32e-11. +! discretization E to be 1.2513e-11 (at local minimum of E/r²). ! ! Arguments are described in subroutine reconstruct(). ! @@ -4568,8 +4568,8 @@ module interpolations real(kind=8), dimension(size(fc)) :: r real(kind=8), dimension(size(fc)) :: u - real(kind=8), parameter :: a1 = 5.0087697436175480833840357d-01 - real(kind=8), parameter :: a2 = 3.4091832657841234002695365d-01 + real(kind=8), parameter :: a1 = 5.00415305771646324666320357024935d-01 + real(kind=8), parameter :: a2 = 3.41203006456869647299948174549092d-01 real(kind=8), dimension(3), parameter :: di9 = [ a1, 1.0d+00, a2 ] real(kind=8), dimension(9), parameter :: & ci9 = [ - 5.000d+00 * a1 - 5.000d+00 * a2 + 4.000d+00, & @@ -4679,7 +4679,7 @@ module interpolations ! The method uses an implicit penta-diagonal solver and its coefficients are ! optimized for the maximum value of imaginary part of the scaled modified ! wave number (Im(κ)) to be 2e-6 and the integrated error of compact -! discretization E to be 1.38e-11. +! discretization E to be 1.2984e-11 (at local minimum of E/r). ! ! Arguments are described in subroutine reconstruct(). ! @@ -4703,8 +4703,8 @@ module interpolations real(kind=8), dimension(size(fc)) :: r real(kind=8), dimension(size(fc)) :: u - real(kind=8), parameter :: a1 = 6.7100477315625665203995935d-01 - real(kind=8), parameter :: a2 = 3.4097048071278782662422942d-01 + real(kind=8), parameter :: a1 = 6.70862177813208397304414380169977d-01 + real(kind=8), parameter :: a2 = 3.41040434946989987878686778898846d-01 real(kind=8), dimension(5), parameter :: & di7 = [ (3.0d+00 * a1 + 2.0d+00 * a2 - 2.0d+00) / 8.0d+00, & a1, 1.0d+00, a2, & @@ -4815,7 +4815,7 @@ module interpolations ! The method uses an implicit penta-diagonal solver and its coefficients are ! optimized for the maximum value of imaginary part of the scaled modified ! wave number (Im(κ)) to be 2e-6 and the integrated error of compact -! discretization E to be 1.54e-11. +! discretization E to be 1.4776e-11 (at local minimum of E/r). ! ! Arguments are described in subroutine reconstruct(). ! @@ -4839,8 +4839,8 @@ module interpolations real(kind=8), dimension(size(fc)) :: r real(kind=8), dimension(size(fc)) :: u - real(kind=8), parameter :: a1 = 6.7028117788580347100541103d-01 - real(kind=8), parameter :: a2 = 4.1002518667894543684869155d-01 + real(kind=8), parameter :: a1 = 6.69938726201951413971659371944416d-01 + real(kind=8), parameter :: a2 = 4.10221751253870603082971042746238d-01 real(kind=8), dimension(5), parameter :: & di9 = [ (9.0d+00 * a1 + 5.0d+00 * a2 - 6.0d+00) / 2.0d+01, & a1, 1.0d+00, a2, &