fix(time_integration): update RK4 function

Author: rouson

src/formal/time_integration_t.F90

@@ -2,7 +2,7 @@ module time_integration_m
   implicit none
 
 #ifndef HAVE_TEMPLATE_SUPPORT
-  integer dummy
+  logical prevent_empty_module
 #else
 
   requirement unary_op_r(op, T, U)
@@ -38,21 +38,20 @@ function op(x,y) result(z)
   end requirement
 
   template runge_kutta_4th_order_t(    &
-     R, TFA, TFB, TFC, TFD     & ! types
-    ,rmult, rtmult, tplus, rhs & ! functions
+     TFA, TFB, TFC, TFD     & ! types
+    ,tmult, rtmult, taplus, tcplus, rhs & ! functions
   )
-     type, deferred :: R
      type, deferred :: TFA
      type, deferred :: TFB
      type, deferred :: TFC
      type, deferred :: TFD
 
-     require :: binary_op_r(integer,     TFC, TFC, tmult)
-     require :: binary_op_r(      R, integer,   R, rdiv)
-     require :: binary_op_r(      R,     TFC, TFB, rtmult)
+     require :: binary_op_r(integer,     TFC, TFC, tmult) ! report lfortran bug on intrinsic types
+     require :: binary_op_r(    TFC,       R, TFB, rtmult)
      require :: binary_op_r(    TFA,     TFB, TFA, taplus)
      require :: binary_op_r(    TFC,     TFC, TFC, tcplus)
      require :: binary_op_r(    TFA,     TFD, TFC, rhs)
+     require :: binary_op_r(      R, integer,   R, rdiv)
 
      interface operator(+)
        module procedure taplus, tcplus
@@ -71,19 +70,17 @@ function op(x,y) result(z)
     pure function step(s, v, dt) result(s_next)
       type(TFA), intent(in) :: s
       type(TFD), intent(in) :: v
-      type(R),   intent(in) :: dt
+      type(R)  , intent(in) :: dt
       type(TFA) s_next
 
-      associate(k1 => rhs(s, v)) 
-        associate(k2 => rhs(s + (dt/2)*k1, v)) 
-          associate(k3 => rhs(s + (dt/2)*k2, v)) 
-            associate(k4 => rhs(s + dt*k3, v)) 
-              s_next = s + (dt/6)*(k1 + 2*k2 + 2*k3 + k4) 
-            end associate
-          end associate
-        end associate
-      end associate
-  end function
+      associate(k1 => rhs(s, v))
+      associate(k2 => rhs(s + k1*(dt/2), v))
+      associate(k3 => rhs(s + k2*(dt/2), v))
+      associate(k4 => rhs(s + k3* dt, v))
+        s_next = s + (k1 + 2*k2 + 2*k3 + k4) * (dt/6)
+       !s_next = s + k1*(dt/6) + k2*(dt/3) + k3*(dt/3) + k4*(dt/6)
+      end associate; end associate; end associate; end associate
+    end function
 
   end template