From 99d5196736b8cbc71fb1cc7b40460b78f8d979d7 Mon Sep 17 00:00:00 2001
From: Grzegorz Kowal <grzegorz@amuncode.org>
Date: Mon, 27 Feb 2023 10:04:06 -0300
Subject: [PATCH] FORCING: Allow to specify the period of turbulence injection.

The period is controlled by two parameters:

'injection_time_start' - the beginning of the injection period,
'injection_time_stop'  - the end of the injection period.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
---
 sources/evolution.F90 |  2 +-
 sources/forcing.F90   | 51 +++++++++++++++++++++++++++++++++----------
 2 files changed, 41 insertions(+), 12 deletions(-)

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()