Merge branch 'master' into reconnection

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
 !

sources/forcing.F90

@@ -35,8 +35,8 @@ module forcing
 ! interfaces for procedure pointers
 !
   abstract interface
-    subroutine update_forcing_iface(dt)
+    subroutine update_forcing_iface(t, dt)
-      real(kind=8), intent(in)  :: 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()

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 :: a1 = 5.01241562476900742404893381352899d-01
-    real(kind=8), parameter :: a2 = 2.5452055302783787784904423d-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 :: a1 = 5.00765434358710549355767793738926d-01
-    real(kind=8), parameter :: a2 = 3.0649545768561641522180678d-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 :: a1 = 5.00415305771646324666320357024935d-01
-    real(kind=8), parameter :: a2 = 3.4091832657841234002695365d-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 :: a1 = 6.70862177813208397304414380169977d-01
-    real(kind=8), parameter :: a2 = 3.4097048071278782662422942d-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 :: a1 = 6.69938726201951413971659371944416d-01
-    real(kind=8), parameter :: a2 = 4.1002518667894543684869155d-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,                                    &