Index: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/Rosenbrock_vec.f90
===================================================================
--- palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/Rosenbrock_vec.f90	(revision 3298)
+++ palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/Rosenbrock_vec.f90	(revision 3298)
@@ -0,0 +1,1258 @@
+! NOTE: This routine is a vectorised version of Rosenbrock
+!       and can only be used with KP4.
+!
+SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, &
+           AbsTol,RelTol,              &
+           RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!
+!    Solves the system y'=F(t,y) using a Rosenbrock method defined by:
+!
+!     G = 1/(H*gamma(1)) - Jac(t0,Y0)
+!     T_i = t0 + Alpha(i)*H
+!     Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
+!     G * K_i = Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j +
+!         gamma(i)*dF/dT(t0, Y0)
+!     Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
+!
+!    For details on Rosenbrock methods and their implementation consult:
+!      E. Hairer and G. Wanner
+!      "Solving ODEs II. Stiff and differential-algebraic problems".
+!      Springer series in computational mathematics, Springer-Verlag, 1996.
+!    The codes contained in the book inspired this implementation.
+!
+!    (C)  Adrian Sandu, August 2004
+!    Virginia Polytechnic Institute and State University
+!    Contact: sandu@cs.vt.edu
+!    Revised by Philipp Miehe and Adrian Sandu, May 2006                  
+!    This implementation is part of KPP - the Kinetic PreProcessor
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!
+!~~~>   INPUT ARGUMENTS:
+!
+!-     Y(N)    = vector of initial conditions (at T=Tstart)
+!-    [Tstart,Tend]  = time range of integration
+!     (if Tstart>Tend the integration is performed backwards in time)
+!-    RelTol, AbsTol = user precribed accuracy
+!- SUBROUTINE  Fun( T, Y, Ydot ) = ODE function,
+!                       returns Ydot = Y' = F(T,Y)
+!- SUBROUTINE  Jac( T, Y, Jcb ) = Jacobian of the ODE function,
+!                       returns Jcb = dFun/dY
+!-    ICNTRL(1:20)    = integer inputs parameters
+!-    RCNTRL(1:20)    = real inputs parameters
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!
+!~~~>     OUTPUT ARGUMENTS:
+!
+!-    Y(N)    -> vector of final states (at T->Tend)
+!-    ISTATUS(1:20)   -> integer output parameters
+!-    RSTATUS(1:20)   -> real output parameters
+!-    IERR            -> job status upon return
+!                        success (positive value) or
+!                        failure (negative value)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!
+!~~~>     INPUT PARAMETERS:
+!
+!    Note: For input parameters equal to zero the default values of the
+!       corresponding variables are used.
+!
+!    ICNTRL(1) = 1: F = F(y)   Independent of T (AUTONOMOUS)
+!              = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
+!
+!    ICNTRL(2) = 0: AbsTol, RelTol are N-dimensional vectors
+!              = 1: AbsTol, RelTol are scalars
+!
+!    ICNTRL(3)  -> selection of a particular Rosenbrock method
+!        = 0 :    Rodas3 (default)
+!        = 1 :    Ros2
+!        = 2 :    Ros3
+!        = 3 :    Ros4F
+!        = 4 :    Rodas3
+!        = 5 :    Rodas4
+!
+!    ICNTRL(4)  -> maximum number of integration steps
+!        For ICNTRL(4)=0) the default value of 100000 is used
+!
+!    RCNTRL(1)  -> Hmin, lower bound for the integration step size
+!          It is strongly recommended to keep Hmin = ZERO
+!    RCNTRL(2)  -> Hmax, upper bound for the integration step size
+!    RCNTRL(3)  -> Hstart, starting value for the integration step size
+!
+!    RCNTRL(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
+!    RCNTRL(5)  -> FacMax, upper bound on step increase factor (default=6)
+!    RCNTRL(6)  -> FacRej, step decrease factor after multiple rejections
+!                          (default=0.1)
+!    RCNTRL(7)  -> FacSafe, by which the new step is slightly smaller
+!         than the predicted value  (default=0.9)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!
+!
+!    OUTPUT ARGUMENTS:
+!    -----------------
+!
+!    T           -> T value for which the solution has been computed
+!                     (after successful return T=Tend).
+!
+!    Y(N)        -> Numerical solution at T
+!
+!    IDID        -> Reports on successfulness upon return:
+!                    = 1 for success
+!                    < 0 for error (value equals error code)
+!
+!    ISTATUS(1)  -> No. of function calls
+!    ISTATUS(2)  -> No. of jacobian calls
+!    ISTATUS(3)  -> No. of steps
+!    ISTATUS(4)  -> No. of accepted steps
+!    ISTATUS(5)  -> No. of rejected steps (except at very beginning)
+!    ISTATUS(6)  -> No. of LU decompositions
+!    ISTATUS(7)  -> No. of forward/backward substitutions
+!    ISTATUS(8)  -> No. of singular matrix decompositions
+!
+!    RSTATUS(1)  -> Texit, the time corresponding to the
+!                     computed Y upon return
+!    RSTATUS(2)  -> Hexit, last accepted step before exit
+!    RSTATUS(3)  -> Hnew, last predicted step (not yet taken)
+!                   For multiple restarts, use Hnew as Hstart 
+!                     in the subsequent run
+!
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  USE cpulog,             ONLY: cpu_log, log_point_s
+
+  IMPLICIT NONE
+
+!~~~>  Arguments
+   INTEGER,       INTENT(IN)    :: N
+   REAL(kind=dp), INTENT(INOUT) :: Y(:,:)
+   REAL(kind=dp), INTENT(IN)    :: Tstart,Tend
+   REAL(kind=dp), INTENT(IN)    :: AbsTol(N),RelTol(N)
+   INTEGER,       INTENT(IN)    :: ICNTRL(20)
+   REAL(kind=dp), INTENT(IN)    :: RCNTRL(20)
+   INTEGER,       INTENT(INOUT) :: ISTATUS(20)
+   REAL(kind=dp), INTENT(INOUT) :: RSTATUS(20)
+   INTEGER, INTENT(OUT)   :: IERR
+!~~~>  Parameters of the Rosenbrock method, up to 6 stages
+   INTEGER ::  ros_S, rosMethod
+   INTEGER, PARAMETER :: RS2=1, RS3=2, RS4=3, RD3=4, RD4=5, RG3=6
+   REAL(kind=dp) :: ros_A(15), ros_C(15), ros_M(6), ros_E(6), &
+                    ros_Alpha(6), ros_Gamma(6), ros_ELO
+   LOGICAL :: ros_NewF(6)
+   CHARACTER(LEN=12) :: ros_Name
+!~~~>  Local variables
+   REAL(kind=dp) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
+   REAL(kind=dp) :: Hmin, Hmax, Hstart
+   REAL(kind=dp) :: Texit
+   INTEGER       :: i, UplimTol, Max_no_steps
+   LOGICAL       :: Autonomous, VectorTol
+!~~~>   Parameters
+   REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp, ONE  = 1.0_dp
+   REAL(kind=dp), PARAMETER :: DeltaMin = 1.0E-5_dp
+   INTEGER        :: solver_method
+
+
+!~~~>  Autonomous or time dependent ODE. Default is time dependent.
+   Autonomous = .NOT.(ICNTRL(1) == 0)
+
+!~~~>  For Scalar tolerances (ICNTRL(2).NE.0)  the code uses AbsTol(1) and RelTol(1)
+!   For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:N) and RelTol(1:N)
+   IF (ICNTRL(2) == 0) THEN
+      VectorTol = .TRUE.
+      UplimTol  = N
+   ELSE
+      VectorTol = .FALSE.
+      UplimTol  = 1
+   END IF
+
+   solver_method = ICNTRL(3)
+
+!~~~>   Initialize the particular Rosenbrock method selected
+   SELECT CASE (solver_method)
+     CASE (1)
+       CALL Ros2
+     CASE (2)
+       CALL Ros3
+     CASE (3)
+       CALL Ros4
+     CASE (0,4)
+       CALL Rodas3
+     CASE (5)
+       CALL Rodas4
+     CASE (6)
+       CALL Rang3
+     CASE DEFAULT
+       PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=',ICNTRL(3) 
+       CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
+       IERRV(:) = -2
+       RETURN
+   END SELECT
+
+!~~~>   The maximum number of steps admitted
+   IF (ICNTRL(4) == 0) THEN
+      Max_no_steps = 200000
+   ELSEIF (ICNTRL(4) > 0) THEN
+      Max_no_steps=ICNTRL(4)
+   ELSE
+      PRINT * ,'User-selected max no. of steps: ICNTRL(4)=',ICNTRL(4)
+      CALL ros_ErrorMsg(-1,Tstart,ZERO,IERR)
+      IERRV(:) = -1
+      RETURN
+   END IF
+
+!~~~>  Unit roundoff (1+Roundoff>1)
+   Roundoff = WLAMCH('E')
+
+!~~~>  Lower bound on the step size: (positive value)
+   IF (RCNTRL(1) == ZERO) THEN
+      Hmin = ZERO
+   ELSEIF (RCNTRL(1) > ZERO) THEN
+      Hmin = RCNTRL(1)
+   ELSE
+      PRINT * , 'User-selected Hmin: RCNTRL(1)=', RCNTRL(1)
+      CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
+      IERRV(:) = -3
+      RETURN
+   END IF
+!~~~>  Upper bound on the step size: (positive value)
+   IF (RCNTRL(2) == ZERO) THEN
+      Hmax = ABS(Tend-Tstart)
+   ELSEIF (RCNTRL(2) > ZERO) THEN
+      Hmax = MIN(ABS(RCNTRL(2)),ABS(Tend-Tstart))
+   ELSE
+      PRINT * , 'User-selected Hmax: RCNTRL(2)=', RCNTRL(2)
+      CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
+      IERRV(:) = -3
+      RETURN
+   END IF
+!~~~>  Starting step size: (positive value)
+   IF (RCNTRL(3) == ZERO) THEN
+      Hstart = MAX(Hmin,DeltaMin)
+   ELSEIF (RCNTRL(3) > ZERO) THEN
+      Hstart = MIN(ABS(RCNTRL(3)),ABS(Tend-Tstart))
+   ELSE
+      PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
+      CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
+      RETURN
+   END IF
+!~~~>  Step size can be changed s.t.  FacMin < Hnew/Hold < FacMax
+   IF (RCNTRL(4) == ZERO) THEN
+      FacMin = 0.2_dp
+   ELSEIF (RCNTRL(4) > ZERO) THEN
+      FacMin = RCNTRL(4)
+   ELSE
+      PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
+      CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
+      IERRV(:) = -4
+      RETURN
+   END IF
+   IF (RCNTRL(5) == ZERO) THEN
+      FacMax = 6.0_dp
+   ELSEIF (RCNTRL(5) > ZERO) THEN
+      FacMax = RCNTRL(5)
+   ELSE
+      PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
+      CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
+      IERRV(:) = -4
+      RETURN
+   END IF
+!~~~>   FacRej: Factor to decrease step after 2 succesive rejections
+   IF (RCNTRL(6) == ZERO) THEN
+      FacRej = 0.1_dp
+   ELSEIF (RCNTRL(6) > ZERO) THEN
+      FacRej = RCNTRL(6)
+   ELSE
+      PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
+      CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
+      IERRV(:) = -4
+      RETURN
+   END IF
+!~~~>   FacSafe: Safety Factor in the computation of new step size
+   IF (RCNTRL(7) == ZERO) THEN
+      FacSafe = 0.9_dp
+   ELSEIF (RCNTRL(7) > ZERO) THEN
+      FacSafe = RCNTRL(7)
+   ELSE
+      PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
+      CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
+      IERRV(:) = -4
+      RETURN
+   END IF
+!~~~>  Check if tolerances are reasonable
+    DO i=1,UplimTol
+      IF ( (AbsTol(i) <= ZERO) .OR. (RelTol(i) <= 10.0_dp*Roundoff) &
+         .OR. (RelTol(i) >= 1.0_dp) ) THEN
+        PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
+        PRINT * , ' RelTol(',i,') = ',RelTol(i)
+        CALL ros_ErrorMsg(-5,Tstart,ZERO,IERR)
+        RETURN
+      END IF
+    END DO
+
+
+!~~~>  CALL Rosenbrock method
+   CALL ros_Integrator(Y, Tstart, Tend, Texit,   &
+        AbsTol, RelTol,                          &
+!  Integration parameters
+        Autonomous, VectorTol, Max_no_steps,     &
+        Roundoff, Hmin, Hmax, Hstart,            &
+        FacMin, FacMax, FacRej, FacSafe,         &
+!  Error indicator
+        IERR)
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+CONTAINS !  SUBROUTINES internal to Rosenbrock
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+   
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
+ SUBROUTINE ros_ErrorMsg(Code,T,H,IERR)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
+!    Handles all error messages
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
+  
+   REAL(kind=dp), INTENT(IN) :: T, H
+   INTEGER, INTENT(IN)  :: Code
+   INTEGER, INTENT(OUT) :: IERR
+   
+   IERR = Code
+   PRINT * , &
+     'Forced exit from Rosenbrock due to the following error:' 
+     
+   SELECT CASE (Code)
+    CASE (-1)    
+      PRINT * , '--> Improper value for maximal no of steps'
+    CASE (-2)    
+      PRINT * , '--> Selected Rosenbrock method not implemented'
+    CASE (-3)    
+      PRINT * , '--> Hmin/Hmax/Hstart must be positive'
+    CASE (-4)    
+      PRINT * , '--> FacMin/FacMax/FacRej must be positive'
+    CASE (-5) 
+      PRINT * , '--> Improper tolerance values'
+    CASE (-6) 
+      PRINT * , '--> No of steps exceeds maximum bound'
+    CASE (-7) 
+      PRINT * , '--> Step size too small: T + 10*H = T', &
+            ' or H < Roundoff'
+    CASE (-8)    
+      PRINT * , '--> Matrix is repeatedly singular'
+    CASE DEFAULT
+      PRINT *, 'Unknown Error code: ', Code
+   END SELECT
+   
+!   PRINT *, "T=", T, "and H=", H
+     
+ END SUBROUTINE ros_ErrorMsg
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ SUBROUTINE ros_Integrator (Y, Tstart, Tend, Tout,  &
+        AbsTol, RelTol,                          &
+!~~~> Integration parameters
+        Autonomous, VectorTol, Max_no_steps,     &
+        Roundoff, Hmin, Hmax, Hstart,            &
+        FacMin, FacMax, FacRej, FacSafe,         &
+!~~~> Error indicator
+        IERR_out )
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!   Template for the implementation of a generic Rosenbrock method
+!      defined by ros_S (no of stages)
+!      and its coefficients ros_{A,C,M,E,Alpha,Gamma}
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+  USE kp4_compress,    ONLY: kco_initialize, kco_compress, &
+                             kco_finalize, cell_done
+  IMPLICIT NONE
+
+!~~~> Input: the initial condition at Tstart; Output: the solution at T
+   REAL(kind=dp), INTENT(INOUT) :: Y(:,:)
+!~~~> Input: integration interval
+   REAL(kind=dp), INTENT(IN) :: Tstart,Tend
+!~~~> Output: time at which the solution is returned (T=Tend if success)
+   REAL(kind=dp), INTENT(OUT) ::  Tout
+!~~~> Input: tolerances
+   REAL(kind=dp), INTENT(IN) ::  AbsTol(N), RelTol(N)
+!~~~> Input: integration parameters
+   LOGICAL, INTENT(IN) :: Autonomous, VectorTol
+   REAL(kind=dp), INTENT(IN) :: Hstart, Hmin, Hmax
+   INTEGER, INTENT(IN) :: Max_no_steps
+   REAL(kind=dp), INTENT(IN) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
+!~~~> Output: Error indicator
+   INTEGER, INTENT(OUT) :: IERR_out
+! ~~~~ Local variables
+   REAL(kind=dp) :: Ynew(VL_DIM,N), Fcn0(VL_DIM,N), Fcn(VL_DIM,N)
+   REAL(kind=dp) :: K(VL_DIM,N*ros_S), dFdT(VL_DIM,N)
+   REAL(kind=dp) :: Jac0(VL_DIM,LU_NONZERO), Ghimj(VL_DIM,LU_NONZERO)
+!kk   REAL(kind=dp) :: H, Hnew, HC, HG, Fac, Tau
+   REAL(kind=dp) :: Yerr(VL_DIM,N)
+   REAL(kind=dp) :: Y_save(VL_DIM,N),Ynew_save(VL_DIM,N)
+   REAL(kind=dp) :: Hnew_save(VL_DIM)
+   INTEGER :: Pivot(N), Direction, ioffset, j, istage, i
+   LOGICAL :: Singular
+!~~~>  Local parameters
+   REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp, ONE  = 1.0_dp
+   REAL(kind=dp), PARAMETER :: DeltaMin = 1.0E-5_dp
+
+   LOGICAL,SAVE             :: lfirst = .TRUE.
+!~~~>  Locally called functions
+!    REAL(kind=dp) WLAMCH
+!    EXTERNAL WLAMCH
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+!  vector compression algorithm
+
+   REAL(kind=dp), dimension(VL_dim)      :: T
+   REAL(kind=dp), dimension(VL_dim)      :: H, Hnew, HC, HC2, HG, Fac, Tau, Err,ros_v
+   LOGICAL, dimension(VL_dim)            :: RejectLastH, RejectMoreH
+   LOGICAL, dimension(VL_dim)            :: Accept_step,y_copied
+   integer                               :: vl_save
+   integer                               :: loop_count
+
+   if (lfirst)  then
+      write(9,'(a,i4,/)') 'kpp4palm vector version   ==> Vector length = ',vl_dim
+   end if
+
+!~~~>  Initial preparations
+   T(1:VL) = Tstart
+   RSTATUS(Nhexit) = ZERO
+   Tout    = Tend
+   H(1:VL) = MIN( MAX(ABS(Hmin),ABS(Hstart)) , ABS(Hmax) )
+   where (ABS(H(1:VL)) <= 10.0_dp*Roundoff) H(1:VL) = DeltaMin
+
+
+   IF (Tend  >=  Tstart) THEN
+     Direction = +1
+   ELSE
+     Direction = -1
+   END IF
+   H = Direction*H
+
+   RejectLastH(1:VL)=.FALSE.
+   RejectMoreH(1:VL)=.FALSE.
+
+   loop_count = 0
+!~~~> Time loop begins below
+
+   IERRV(1:VL) = 1  ! successful integration
+   vl_save          = VL
+   Kacc       = 0
+   Krej       = 0
+
+!kk   IF (.not. l_fixed_step) then
+      call kco_initialize (VL, y, IERRV, Kacc, Krej)
+!kk   end if
+
+TimeLoop: DO
+!kk Eventuell als Where einbauen
+   where ( (Direction > 0).AND.((T(1:VL)-Tend)+Roundoff <= ZERO) &
+        .OR. (Direction < 0).AND.((Tend-T(1:VL))+Roundoff <= ZERO) )
+      cell_done(1:VL) = .FALSE.
+   elsewhere
+      cell_done(1:VL) = .TRUE.
+   end where
+!   do i=1,vl
+!      if ( (Direction > 0).AND.((T(i)-Tend)+Roundoff <= ZERO) &
+!        .OR. (Direction < 0).AND.((Tend-T(i))+Roundoff <= ZERO) )  then
+!         cell_done(i) = .FALSE.
+!      else
+!         cell_done(i) = .TRUE.
+!       end if
+!   end do
+
+   if(count(cell_done(1:VL)) == vl_save)   then                      ! EXIT 1, all cells converge at the same time
+      EXIT                                                           ! Iteration terminates without without compress
+   end if
+
+   loop_count = loop_count+1
+
+   call kco_compress (VL, T, H, Hnew, IERRV, y, RCONST, FIX, RejectLastH, RejectMoreH, Kacc, Krej)
+
+   if(all(cell_done(1:VL)))  then                                    ! EXIT 2, all cells converge have been converged
+     EXIT
+   end if
+
+   IF ( ISTATUS(Nstp) > Max_no_steps ) THEN  ! Too many steps
+      CALL ros_ErrorMsg(-6,T(1),H(1),IERR)
+      IERRV(1:VL) = -6
+      RETURN
+   END IF
+
+   IERRV = 0
+   WHERE ( ((T(1:VL)+0.1_dp*H(1:VL)) == T(1)).OR.(H(1:VL) <= Roundoff) )   ! Step size too small
+      IERRV(1:VL) = -7
+   END WHERE
+
+   IF(MINVAL(IERRV) == -7)   then
+      CALL ros_ErrorMsg(-7,T(1),H(1),IERR)
+      RETURN
+    END IF
+
+!~~~>  Limit H if necessary to avoid going beyond Tend
+   H(1:VL) = MIN(H(1:VL),ABS(Tend-T(1:VL)))
+
+!~~~>   Compute the function at current time
+   CALL FunTemplate(T,Y,Fcn0)
+   ISTATUS(Nfun) = ISTATUS(Nfun) + 1
+
+!~~~>  Compute the function derivative with respect to T
+   IF (.NOT.Autonomous) THEN
+      CALL ros_FunTimeDerivative ( T, Roundoff, Y, &
+                Fcn0, dFdT )
+   END IF
+
+!~~~>   Compute the Jacobian at current time
+   CALL JacTemplate(T,Y,Jac0)
+   ISTATUS(Njac) = ISTATUS(Njac) + 1                            ! This is not correct in vector mode, please use xnacc and xnrej
+
+!~~~>  Repeat step calculation until current step accepted
+    Accept_step(1:VL) = .FALSE.
+    y_copied(1:VL)    = .FALSE.
+
+UntilAccepted: DO
+
+   CALL ros_PrepareMatrix(H,Direction,ros_Gamma(1), &
+          Jac0,Ghimj,Pivot,Singular)
+   IF (Singular) THEN ! More than 5 consecutive failed decompositions
+       CALL ros_ErrorMsg(-8,T(1),H(1),IERR)
+       RETURN
+   END IF
+
+!~~~>   Compute the stages
+Stage: DO istage = 1, ros_S
+
+!      Current istage offset. Current istage vector is K(ioffset+1:ioffset+N)
+       ioffset = N*(istage-1)
+
+!      For the 1st istage the function has been computed previously
+       IF ( istage == 1 ) THEN
+         Fcn(1:VL,1:N) = Fcn0(1:VL,1:N) 
+      ! istage>1 and a new function evaluation is needed at the current istage
+       ELSEIF ( ros_NewF(istage) ) THEN
+	   Ynew(1:VL,1:N) = Y(1:VL,1:N)
+        DO j = 1, istage-1
+!           ros_v = ros_A((istage-1)*(istage-2)/2+j)
+!           CALL WAXPY(N,ros_v,K(:,N*(j-1)+1:N*(j-1)+N),1,Ynew,1)
+            Ynew(1:vl,1:N) = Ynew(1:vl,1:N)+ros_A((istage-1)*(istage-2)/2+j)*K(1:VL,N*(j-1)+1:N*(j-1)+N)
+         END DO
+         Tau = T + ros_Alpha(istage)*Direction*H
+         CALL FunTemplate(Tau,Ynew,Fcn)
+         ISTATUS(Nfun) = ISTATUS(Nfun) + 1
+       END IF ! if istage == 1 elseif ros_NewF(istage)
+       K(1:VL,ioffset+1:ioffset+N) = Fcn(1:VL,1:N)
+       DO j = 1, istage-1
+         HC(1:vl) = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H(1:vl))
+!         CALL WAXPY(N,HC,K(:,N*(j-1)+1:N*(j-1)+N),1,K(:,ioffset+1:ioffset+n),1)
+         DO i=1,N
+            K(1:vl,ioffset+i) = K(1:vl,ioffset+i) + HC(1:vl)*K(1:vl,N*(j-1)+i)
+         END DO
+
+       END DO
+       IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
+         HG(1:vl) = Direction*H(1:vl)*ros_Gamma(istage)
+         CALL WAXPY(N,HG,dFdT,1,K(:,ioffset+1:ioffset+n),1)
+       END IF
+
+       CALL ros_Solve(Ghimj, K(:,ioffset+1:ioffset+n))
+
+   END DO Stage
+
+
+!~~~>  Compute the new solution
+   Ynew(1:vl,1:N) = Y(1:vl,1:N)
+   DO j=1,ros_S
+         ros_v = ros_M(j)
+         CALL WAXPY(N,ros_v,K(:,N*(j-1)+1:),1,Ynew,1)
+   END DO
+
+!~~~>  Compute the error estimation
+   Yerr(1:vl,1:N) = ZERO
+   DO j=1,ros_S
+        ros_v = ros_E(j)
+        CALL WAXPY(N,ros_v,K(:,N*(j-1)+1:N*(j-1)+N),1,Yerr,1)
+   END DO
+   Err(1:VL) = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol )
+
+!~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
+   Fac(1:vl)  = MIN(FacMax,MAX(FacMin,FacSafe/Err(1:vl)**(ONE/ros_ELO)))
+   Hnew(1:vl) = H(1:vl)*Fac(1:vl)
+
+!~~~>  Check the error magnitude and adjust step size
+   ISTATUS(Nstp) = ISTATUS(Nstp) + 1
+!   write(9,*) 'Hmin ',Hmin,Hmax
+ do i=1,VL
+    IF(Accept_step(i))   then
+       Ynew(i,:) = Ynew_save(i,:)
+       Hnew(i)   = Hnew_save(i)
+    ELSE IF ( (Err(i) <= ONE).OR.(H(i) <= Hmin) ) THEN  !~~~> Accept step
+      ISTATUS(Nacc) = ISTATUS(Nacc) + 1
+      Y(i,1:N) = Ynew(i,1:N)
+      T(i) = T(i) + Direction*H(i)
+      Hnew(i) = MAX(Hmin,MIN(Hnew(i),Hmax))
+      IF (RejectLastH(i)) THEN  ! No step size increase after a rejected step
+         Hnew(i) = MIN(Hnew(i),H(i))
+      END IF
+      RSTATUS(Nhexit) = H(i)                    ! ueberdenken
+      RSTATUS(Nhnew)  = Hnew(i)
+      RSTATUS(Ntexit) = T(i)                    ! ueberdenken ende
+      RejectLastH(i) = .FALSE.
+      RejectMoreH(i) = .FALSE.
+      H(i) = Hnew(i)
+
+      Accept_step(i) = .TRUE.                   ! EXIT THE LOOP: WHILE STEP NOT ACCEPTED
+      Kacc(i) = Kacc(i) +1
+   ELSE           !~~~> Reject step
+      IF (RejectMoreH(i)) THEN
+         Hnew (i)= H(i)*FacRej
+      END IF
+      RejectMoreH(i) = RejectLastH(i)
+      RejectLastH(i) = .TRUE.
+      H(i) = Hnew(i)
+
+      IF (ISTATUS(Nacc) >= 1)  ISTATUS(Nrej) = ISTATUS(Nrej) + 1
+      Krej(i) = Krej(i)+1
+
+   END IF ! Err <= 1
+  END DO
+   if(all(Accept_step(1:VL)))   EXIT
+
+   do i=1,VL
+      IF(Accept_step(i))   then
+         if(.not. y_copied(i))  then
+            y_save(i,:)    = y(i,:)
+            Ynew_save(i,:) = Ynew(i,:)
+            Hnew_save(i)   = Hnew(i)
+            y_copied(i) = .TRUE.
+         end if
+      end if
+   end do
+
+
+   END DO UntilAccepted
+
+    do i=1,VL
+       if(y_copied(i))  then
+           y(i,:)  = y_save(i,:)
+       end if
+    end do
+
+   END DO TimeLoop
+
+   call kco_finalize (y, IERRV, Kacc, Krej)
+
+   VL = VL_save
+
+!~~~> Succesful exit
+   IERR_out = 1  !~~~> The integration was successful
+
+   lfirst = .FALSE.
+
+  END SUBROUTINE ros_Integrator
+
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  FUNCTION ros_ErrorNorm ( Y, Ynew, Yerr, &
+               AbsTol, RelTol, VectorTol )
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~> Computes the "scaled norm" of the error vector Yerr
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+   IMPLICIT NONE
+
+! Input arguments
+   REAL(kind=dp), INTENT(IN) :: Y(:,:), Ynew(:,:), &
+          Yerr(:,:), AbsTol(N), RelTol(N)
+   LOGICAL, INTENT(IN) ::  VectorTol
+! Local variables
+   REAL(kind=dp), dimension(VL)       :: Err, Scale, Ymax
+   INTEGER  :: i,j
+   REAL(kind=dp), PARAMETER           :: ZERO = 0.0_dp
+   REAL(kind=dp), dimension(VL)       :: ros_ErrorNorm
+
+   Err = ZERO
+   DO i=1,NVAR
+       do j=1,VL
+           Ymax(j) = MAX( ABS(Y(j,i)),ABS(Ynew(j,i)))
+           IF (VectorTol) THEN
+               Scale(j) = AbsTol(i)+ RelTol(i)*Ymax(j)
+           ELSE
+               Scale(j) = AbsTol(1)+ RelTol(1)*Ymax(1)
+           END IF
+           Err(j) = Err(j)+ (Yerr(j,i)/Scale(j))**2
+       end do
+   END DO
+   Err  = SQRT(Err/N)
+
+   ros_ErrorNorm = max(err, 1.0d-10)
+
+  END FUNCTION ros_ErrorNorm
+
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, &
+                Fcn0, dFdT )
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~> The time partial derivative of the function by finite differences
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+   IMPLICIT NONE
+
+!~~~> Input arguments
+   REAL(kind=dp), INTENT(IN) :: T(:), Roundoff, Y(:,:), Fcn0(:,:)
+!~~~> Output arguments
+   REAL(kind=dp), INTENT(OUT) :: dFdT(:,:)
+!~~~> Local variables
+   REAL(kind=dp) :: Delta(1:vl_dim),one_v(1:vl_dim)
+   REAL(kind=dp), PARAMETER :: ONE = 1.0_dp, DeltaMin = 1.0E-6_dp
+   INTEGER   :: i
+
+   Delta(1:vl) = SQRT(Roundoff)*MAX(DeltaMin,ABS(T(1:vl)))
+   CALL FunTemplate(T(1:vl)+Delta(1:vl),Y,dFdT)
+   ISTATUS(Nfun) = ISTATUS(Nfun) + 1
+   one_v = -ONE
+   CALL WAXPY(N,one_v,Fcn0,1,dFdT,1)
+!   CALL WSCAL(N,(ONE/Delta),dFdT,1)
+   do i=1,size(dFdT,2)
+      dFdT(1:vl,i) = 1.0_dp/Delta(1:vl)*dFdT(1:vl,i)
+   end do
+
+  END SUBROUTINE ros_FunTimeDerivative
+
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, &
+             Jac0, Ghimj, Pivot, Singular )
+! --- --- --- --- --- --- --- --- --- --- --- --- ---
+!  Prepares the LHS matrix for stage calculations
+!  1.  Construct Ghimj = 1/(H*ham) - Jac0
+!      "(Gamma H) Inverse Minus Jacobian"
+!  2.  Repeat LU decomposition of Ghimj until successful.
+!       -half the step size if LU decomposition fails and retry
+!       -exit after 5 consecutive fails
+! --- --- --- --- --- --- --- --- --- --- --- --- ---
+   IMPLICIT NONE
+
+!~~~> Input arguments
+   REAL(kind=dp), INTENT(IN) ::  Jac0(:,:)
+   REAL(kind=dp), INTENT(IN) ::  gam
+   INTEGER, INTENT(IN) ::  Direction
+!~~~> Output arguments
+   REAL(kind=dp), INTENT(OUT) :: Ghimj(:,:)
+
+   LOGICAL, INTENT(OUT) ::  Singular
+   INTEGER, INTENT(OUT) ::  Pivot(N)
+!~~~> Inout arguments
+   REAL(kind=dp), INTENT(INOUT) :: H(:)   ! step size is decreased when LU fails
+!~~~> Local variables
+   INTEGER  :: i, ISING, Nconsecutive
+   REAL(kind=dp) :: ghinv(1:VL)
+   REAL(kind=dp), PARAMETER :: ONE  = 1.0_dp, HALF = 0.5_dp
+
+   Nconsecutive = 0
+   Singular = .FALSE.
+
+!kk   DO WHILE (Singular)                        !test of Singular currently disabled
+
+!~~~>    Construct Ghimj = 1/(H*gam) - Jac0
+   Ghimj(1:VL,1:LU_NONZERO) = -JAC0(1:VL,1:LU_NONZERO)
+   ghinv(1:VL) = 1.0_dp/(H(1:VL)*gam)
+   DO i=1,N
+     Ghimj(1:VL,LU_DIAG(i)) = Ghimj(1:VL,LU_DIAG(i))+ghinv(1:VL)
+   END DO
+
+!~~~>    Compute LU decomposition
+   CALL KppDecomp( Ghimj, ising )
+
+!kk Original code, implement singular check if necessary
+!!~~~>    Compute LU decomposition
+!     CALL ros_Decomp( Ghimj, Pivot, ISING )
+!     IF (ISING == 0) THEN
+!!~~~>    If successful done
+!        Singular = .FALSE.
+!     ELSE ! ISING .ne. 0
+!!~~~>    If unsuccessful half the step size; if 5 consecutive fails then return
+!        ISTATUS(Nsng) = ISTATUS(Nsng) + 1
+!        Nconsecutive = Nconsecutive+1
+!        Singular = .TRUE.
+!        PRINT*,'Warning: LU Decomposition returned ISING = ',ISING
+!        IF (Nconsecutive <= 5) THEN ! Less than 5 consecutive failed decompositions
+!           H = H*HALF
+!        ELSE  ! More than 5 consecutive failed decompositions
+!           RETURN
+!        END IF  ! Nconsecutive
+!      END IF    ! ISING
+
+!    END DO ! WHILE Singular
+
+     Pivot = 0                                !Set to avoid compiler warnings
+
+  END SUBROUTINE ros_PrepareMatrix
+
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!  SUBROUTINE ros_Decomp( A, Pivot, ISING )
+!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!!  Template for the LU decomposition
+!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!   IMPLICIT NONE
+!!~~~> Inout variables
+!#ifdef FULL_ALGEBRA
+!   REAL(kind=dp), INTENT(INOUT) :: A(N,N)
+!#else
+!   REAL(kind=dp), INTENT(INOUT) :: A(LU_NONZERO)
+!#endif
+!!~~~> Output variables
+!   INTEGER, INTENT(OUT) :: Pivot(N), ISING
+!
+!#ifdef FULL_ALGEBRA
+!   CALL  DGETRF( N, N, A, N, Pivot, ISING )
+!#else
+!   CALL KppDecomp ( A, ISING )
+!   Pivot(1) = 1
+!#endif
+!   ISTATUS(Ndec) = ISTATUS(Ndec) + 1
+!
+!  END SUBROUTINE ros_Decomp
+
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE ros_Solve( A, b )
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!  Template for the forward/backward substitution (using pre-computed LU decomposition)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+   IMPLICIT NONE
+!~~~> Input variables
+   REAL(kind=dp), INTENT(IN)         :: A(:,:)
+!~~~> InOut variables
+   REAL(kind=dp), INTENT(INOUT)      :: b(:,:)
+
+   CALL KppSolve( A, b )
+
+   ISTATUS(Nsol) = ISTATUS(Nsol) + VL
+
+  END SUBROUTINE ros_Solve
+
+
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE Ros2
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! --- AN L-STABLE METHOD, 2 stages, order 2
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+   IMPLICIT NONE
+   DOUBLE PRECISION g
+
+    g = 1.0_dp + 1.0_dp/SQRT(2.0_dp)
+    rosMethod = RS2
+!~~~> Name of the method
+    ros_Name = 'ROS-2'
+!~~~> Number of stages
+    ros_S = 2
+
+!~~~> The coefficient matrices A and C are strictly lower triangular.
+!   The lower triangular (subdiagonal) elements are stored in row-wise order:
+!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
+!   The general mapping formula is:
+!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
+!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
+
+    ros_A(1) = (1.0_dp)/g
+    ros_C(1) = (-2.0_dp)/g
+!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
+!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
+    ros_NewF(1) = .TRUE.
+    ros_NewF(2) = .TRUE.
+!~~~> M_i = Coefficients for new step solution
+    ros_M(1)= (3.0_dp)/(2.0_dp*g)
+    ros_M(2)= (1.0_dp)/(2.0_dp*g)
+! E_i = Coefficients for error estimator
+    ros_E(1) = 1.0_dp/(2.0_dp*g)
+    ros_E(2) = 1.0_dp/(2.0_dp*g)
+!~~~> ros_ELO = estimator of local order - the minimum between the
+!    main and the embedded scheme orders plus one
+    ros_ELO = 2.0_dp
+!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
+    ros_Alpha(1) = 0.0_dp
+    ros_Alpha(2) = 1.0_dp
+!~~~> Gamma_i = \sum_j  gamma_{i,j}
+    ros_Gamma(1) = g
+    ros_Gamma(2) =-g
+
+ END SUBROUTINE Ros2
+
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE Ros3
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+   IMPLICIT NONE
+   rosMethod = RS3
+!~~~> Name of the method
+   ros_Name = 'ROS-3'
+!~~~> Number of stages
+   ros_S = 3
+
+!~~~> The coefficient matrices A and C are strictly lower triangular.
+!   The lower triangular (subdiagonal) elements are stored in row-wise order:
+!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
+!   The general mapping formula is:
+!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
+!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
+
+   ros_A(1)= 1.0_dp
+   ros_A(2)= 1.0_dp
+   ros_A(3)= 0.0_dp
+
+   ros_C(1) = -0.10156171083877702091975600115545E+01_dp
+   ros_C(2) =  0.40759956452537699824805835358067E+01_dp
+   ros_C(3) =  0.92076794298330791242156818474003E+01_dp
+!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
+!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
+   ros_NewF(1) = .TRUE.
+   ros_NewF(2) = .TRUE.
+   ros_NewF(3) = .FALSE.
+!~~~> M_i = Coefficients for new step solution
+   ros_M(1) =  0.1E+01_dp
+   ros_M(2) =  0.61697947043828245592553615689730E+01_dp
+   ros_M(3) = -0.42772256543218573326238373806514_dp
+! E_i = Coefficients for error estimator
+   ros_E(1) =  0.5_dp
+   ros_E(2) = -0.29079558716805469821718236208017E+01_dp
+   ros_E(3) =  0.22354069897811569627360909276199_dp
+!~~~> ros_ELO = estimator of local order - the minimum between the
+!    main and the embedded scheme orders plus 1
+   ros_ELO = 3.0_dp
+!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
+   ros_Alpha(1)= 0.0_dp
+   ros_Alpha(2)= 0.43586652150845899941601945119356_dp
+   ros_Alpha(3)= 0.43586652150845899941601945119356_dp
+!~~~> Gamma_i = \sum_j  gamma_{i,j}
+   ros_Gamma(1)= 0.43586652150845899941601945119356_dp
+   ros_Gamma(2)= 0.24291996454816804366592249683314_dp
+   ros_Gamma(3)= 0.21851380027664058511513169485832E+01_dp
+
+  END SUBROUTINE Ros3
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE Ros4
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!     L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
+!     L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
+!
+!      E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
+!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
+!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
+!      SPRINGER-VERLAG (1990)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+   IMPLICIT NONE
+
+   rosMethod = RS4
+!~~~> Name of the method
+   ros_Name = 'ROS-4'
+!~~~> Number of stages
+   ros_S = 4
+
+!~~~> The coefficient matrices A and C are strictly lower triangular.
+!   The lower triangular (subdiagonal) elements are stored in row-wise order:
+!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
+!   The general mapping formula is:
+!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
+!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
+
+   ros_A(1) = 0.2000000000000000E+01_dp
+   ros_A(2) = 0.1867943637803922E+01_dp
+   ros_A(3) = 0.2344449711399156_dp
+   ros_A(4) = ros_A(2)
+   ros_A(5) = ros_A(3)
+   ros_A(6) = 0.0_dp
+
+   ros_C(1) =-0.7137615036412310E+01_dp
+   ros_C(2) = 0.2580708087951457E+01_dp
+   ros_C(3) = 0.6515950076447975_dp
+   ros_C(4) =-0.2137148994382534E+01_dp
+   ros_C(5) =-0.3214669691237626_dp
+   ros_C(6) =-0.6949742501781779_dp
+!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
+!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
+   ros_NewF(1)  = .TRUE.
+   ros_NewF(2)  = .TRUE.
+   ros_NewF(3)  = .TRUE.
+   ros_NewF(4)  = .FALSE.
+!~~~> M_i = Coefficients for new step solution
+   ros_M(1) = 0.2255570073418735E+01_dp
+   ros_M(2) = 0.2870493262186792_dp
+   ros_M(3) = 0.4353179431840180_dp
+   ros_M(4) = 0.1093502252409163E+01_dp
+!~~~> E_i  = Coefficients for error estimator
+   ros_E(1) =-0.2815431932141155_dp
+   ros_E(2) =-0.7276199124938920E-01_dp
+   ros_E(3) =-0.1082196201495311_dp
+   ros_E(4) =-0.1093502252409163E+01_dp
+!~~~> ros_ELO  = estimator of local order - the minimum between the
+!    main and the embedded scheme orders plus 1
+   ros_ELO  = 4.0_dp
+!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
+   ros_Alpha(1) = 0.0_dp
+   ros_Alpha(2) = 0.1145640000000000E+01_dp
+   ros_Alpha(3) = 0.6552168638155900_dp
+   ros_Alpha(4) = ros_Alpha(3)
+!~~~> Gamma_i = \sum_j  gamma_{i,j}
+   ros_Gamma(1) = 0.5728200000000000_dp
+   ros_Gamma(2) =-0.1769193891319233E+01_dp
+   ros_Gamma(3) = 0.7592633437920482_dp
+   ros_Gamma(4) =-0.1049021087100450_dp
+
+  END SUBROUTINE Ros4
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE Rodas3
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! --- A STIFFLY-STABLE METHOD, 4 stages, order 3
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+   IMPLICIT NONE
+
+   rosMethod = RD3
+!~~~> Name of the method
+   ros_Name = 'RODAS-3'
+!~~~> Number of stages
+   ros_S = 4
+
+!~~~> The coefficient matrices A and C are strictly lower triangular.
+!   The lower triangular (subdiagonal) elements are stored in row-wise order:
+!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
+!   The general mapping formula is:
+!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
+!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
+
+   ros_A(1) = 0.0_dp
+   ros_A(2) = 2.0_dp
+   ros_A(3) = 0.0_dp
+   ros_A(4) = 2.0_dp
+   ros_A(5) = 0.0_dp
+   ros_A(6) = 1.0_dp
+
+   ros_C(1) = 4.0_dp
+   ros_C(2) = 1.0_dp
+   ros_C(3) =-1.0_dp
+   ros_C(4) = 1.0_dp
+   ros_C(5) =-1.0_dp
+   ros_C(6) =-(8.0_dp/3.0_dp)
+
+!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
+!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
+   ros_NewF(1)  = .TRUE.
+   ros_NewF(2)  = .FALSE.
+   ros_NewF(3)  = .TRUE.
+   ros_NewF(4)  = .TRUE.
+!~~~> M_i = Coefficients for new step solution
+   ros_M(1) = 2.0_dp
+   ros_M(2) = 0.0_dp
+   ros_M(3) = 1.0_dp
+   ros_M(4) = 1.0_dp
+!~~~> E_i  = Coefficients for error estimator
+   ros_E(1) = 0.0_dp
+   ros_E(2) = 0.0_dp
+   ros_E(3) = 0.0_dp
+   ros_E(4) = 1.0_dp
+!~~~> ros_ELO  = estimator of local order - the minimum between the
+!    main and the embedded scheme orders plus 1
+   ros_ELO  = 3.0_dp
+!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
+   ros_Alpha(1) = 0.0_dp
+   ros_Alpha(2) = 0.0_dp
+   ros_Alpha(3) = 1.0_dp
+   ros_Alpha(4) = 1.0_dp
+!~~~> Gamma_i = \sum_j  gamma_{i,j}
+   ros_Gamma(1) = 0.5_dp
+   ros_Gamma(2) = 1.5_dp
+   ros_Gamma(3) = 0.0_dp
+   ros_Gamma(4) = 0.0_dp
+
+  END SUBROUTINE Rodas3
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE Rodas4
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!     STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
+!
+!      E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
+!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
+!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
+!      SPRINGER-VERLAG (1996)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+   IMPLICIT NONE
+
+    rosMethod = RD4
+!~~~> Name of the method
+    ros_Name = 'RODAS-4'
+!~~~> Number of stages
+    ros_S = 6
+
+!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
+    ros_Alpha(1) = 0.000_dp
+    ros_Alpha(2) = 0.386_dp
+    ros_Alpha(3) = 0.210_dp
+    ros_Alpha(4) = 0.630_dp
+    ros_Alpha(5) = 1.000_dp
+    ros_Alpha(6) = 1.000_dp
+
+!~~~> Gamma_i = \sum_j  gamma_{i,j}
+    ros_Gamma(1) = 0.2500000000000000_dp
+    ros_Gamma(2) =-0.1043000000000000_dp
+    ros_Gamma(3) = 0.1035000000000000_dp
+    ros_Gamma(4) =-0.3620000000000023E-01_dp
+    ros_Gamma(5) = 0.0_dp
+    ros_Gamma(6) = 0.0_dp
+
+!~~~> The coefficient matrices A and C are strictly lower triangular.
+!   The lower triangular (subdiagonal) elements are stored in row-wise order:
+!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
+!   The general mapping formula is:  A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
+!                  C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
+
+    ros_A(1) = 0.1544000000000000E+01_dp
+    ros_A(2) = 0.9466785280815826_dp
+    ros_A(3) = 0.2557011698983284_dp
+    ros_A(4) = 0.3314825187068521E+01_dp
+    ros_A(5) = 0.2896124015972201E+01_dp
+    ros_A(6) = 0.9986419139977817_dp
+    ros_A(7) = 0.1221224509226641E+01_dp
+    ros_A(8) = 0.6019134481288629E+01_dp
+    ros_A(9) = 0.1253708332932087E+02_dp
+    ros_A(10) =-0.6878860361058950_dp
+    ros_A(11) = ros_A(7)
+    ros_A(12) = ros_A(8)
+    ros_A(13) = ros_A(9)
+    ros_A(14) = ros_A(10)
+    ros_A(15) = 1.0_dp
+
+    ros_C(1) =-0.5668800000000000E+01_dp
+    ros_C(2) =-0.2430093356833875E+01_dp
+    ros_C(3) =-0.2063599157091915_dp
+    ros_C(4) =-0.1073529058151375_dp
+    ros_C(5) =-0.9594562251023355E+01_dp
+    ros_C(6) =-0.2047028614809616E+02_dp
+    ros_C(7) = 0.7496443313967647E+01_dp
+    ros_C(8) =-0.1024680431464352E+02_dp
+    ros_C(9) =-0.3399990352819905E+02_dp
+    ros_C(10) = 0.1170890893206160E+02_dp
+    ros_C(11) = 0.8083246795921522E+01_dp
+    ros_C(12) =-0.7981132988064893E+01_dp
+    ros_C(13) =-0.3152159432874371E+02_dp
+    ros_C(14) = 0.1631930543123136E+02_dp
+    ros_C(15) =-0.6058818238834054E+01_dp
+
+!~~~> M_i = Coefficients for new step solution
+    ros_M(1) = ros_A(7)
+    ros_M(2) = ros_A(8)
+    ros_M(3) = ros_A(9)
+    ros_M(4) = ros_A(10)
+    ros_M(5) = 1.0_dp
+    ros_M(6) = 1.0_dp
+
+!~~~> E_i  = Coefficients for error estimator
+    ros_E(1) = 0.0_dp
+    ros_E(2) = 0.0_dp
+    ros_E(3) = 0.0_dp
+    ros_E(4) = 0.0_dp
+    ros_E(5) = 0.0_dp
+    ros_E(6) = 1.0_dp
+
+!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
+!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
+    ros_NewF(1) = .TRUE.
+    ros_NewF(2) = .TRUE.
+    ros_NewF(3) = .TRUE.
+    ros_NewF(4) = .TRUE.
+    ros_NewF(5) = .TRUE.
+    ros_NewF(6) = .TRUE.
+
+!~~~> ros_ELO  = estimator of local order - the minimum between the
+!        main and the embedded scheme orders plus 1
+    ros_ELO = 4.0_dp
+
+  END SUBROUTINE Rodas4
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE Rang3
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! STIFFLY-STABLE W METHOD OF ORDER 3, WITH 4 STAGES
+!
+! J. RANG and L. ANGERMANN
+! NEW ROSENBROCK W-METHODS OF ORDER 3
+! FOR PARTIAL DIFFERENTIAL ALGEBRAIC
+!        EQUATIONS OF INDEX 1                                             
+! BIT Numerical Mathematics (2005) 45: 761-787
+!  DOI: 10.1007/s10543-005-0035-y
+! Table 4.1-4.2
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+   IMPLICIT NONE
+
+    rosMethod = RG3
+!~~~> Name of the method
+    ros_Name = 'RANG-3'
+!~~~> Number of stages
+    ros_S = 4
+
+    ros_A(1) = 5.09052051067020d+00;
+    ros_A(2) = 5.09052051067020d+00;
+    ros_A(3) = 0.0d0;
+    ros_A(4) = 4.97628111010787d+00;
+    ros_A(5) = 2.77268164715849d-02;
+    ros_A(6) = 2.29428036027904d-01;
+
+    ros_C(1) = -1.16790812312283d+01;
+    ros_C(2) = -1.64057326467367d+01;
+    ros_C(3) = -2.77268164715850d-01;
+    ros_C(4) = -8.38103960500476d+00;
+    ros_C(5) = -8.48328409199343d-01;
+    ros_C(6) =  2.87009860433106d-01;
+
+    ros_M(1) =  5.22582761233094d+00;
+    ros_M(2) = -5.56971148154165d-01;
+    ros_M(3) =  3.57979469353645d-01;
+    ros_M(4) =  1.72337398521064d+00;
+
+    ros_E(1) = -5.16845212784040d+00;
+    ros_E(2) = -1.26351942603842d+00;
+    ros_E(3) = -1.11022302462516d-16;
+    ros_E(4) =  2.22044604925031d-16;
+
+    ros_Alpha(1) = 0.0d00;
+    ros_Alpha(2) = 2.21878746765329d+00;
+    ros_Alpha(3) = 2.21878746765329d+00;
+    ros_Alpha(4) = 1.55392337535788d+00;
+
+    ros_Gamma(1) =  4.35866521508459d-01;
+    ros_Gamma(2) = -1.78292094614483d+00;
+    ros_Gamma(3) = -2.46541900496934d+00;
+    ros_Gamma(4) = -8.05529997906370d-01;
+
+
+!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
+!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
+    ros_NewF(1) = .TRUE.
+    ros_NewF(2) = .TRUE.
+    ros_NewF(3) = .TRUE.
+    ros_NewF(4) = .TRUE.
+
+!~~~> ros_ELO  = estimator of local order - the minimum between the
+!        main and the embedded scheme orders plus 1
+    ros_ELO = 3.0_dp
+
+  END SUBROUTINE Rang3
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!   End of the set of internal Rosenbrock subroutines
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+END SUBROUTINE Rosenbrock
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+
+
+
Index: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/initialize_kpp_ctrl_template.f90
===================================================================
--- palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/initialize_kpp_ctrl_template.f90	(revision 2718)
+++ palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/initialize_kpp_ctrl_template.f90	(revision 3298)
@@ -32,5 +32,5 @@
 CONTAINS
 
-SUBROUTINE initialize_kpp_ctrl(status, iou, modstr)
+SUBROUTINE initialize_kpp_ctrl(status)
 
   IMPLICIT NONE
@@ -38,6 +38,4 @@
   ! I/O
   INTEGER,          INTENT(OUT) :: status
-  INTEGER,          INTENT(IN)  :: iou     ! logical I/O unit
-  CHARACTER(LEN=*), INTENT(IN)  :: modstr  ! read <modstr>.nml
 
   ! LOCAL
Index: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/module_header
===================================================================
--- palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/module_header	(revision 2718)
+++ palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/module_header	(revision 3298)
@@ -3,5 +3,22 @@
 ! ****** Module chem_gasphase_mod is automatically generated by kpp4palm ******
 !
-!           ********* Please do NOT change this Code *********
+!   ********* Please do NOT change this Code, it will be ovewritten *********
+!
+!------------------------------------------------------------------------------!
+! This file was created by KPP (http://people.cs.vt.edu/asandu/Software/Kpp/)
+! and kpp4palm (created by Klaus Ketelsen). kpp4palm is an adapted version 
+! of KP4 (JÃ¶ckel, P., Kerkweg, A., Pozzer, A., Sander, R., Tost, H., Riede, 
+! H., Baumgaertner, A., Gromov, S., and Kern, B., 2010: Development cycle 2 of
+! the Modular Earth Submodel System (MESSy2), Geosci. Model Dev., 3, 717-752, 
+! https://doi.org/10.5194/gmd-3-717-2010). KP4 is part of the Modular Earth 
+! Submodel System (MESSy), which is is available under the  GNU General Public 
+! License (GPL).
+!
+! KPP is free software; you can redistribute it and/or modify it under the terms 
+! of the General Public Licence as published by the Free Software Foundation; 
+! either version 2 of the License, or (at your option) any later version. 
+! KPP is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; 
+! without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
+! PURPOSE. See the GNU General Public Licence for more details. 
 !
 !------------------------------------------------------------------------------!
@@ -30,10 +47,11 @@
 ! -----------------
 ! $Id: module_header 2460 2017-09-13 14:47:48Z forkel $
+! forkel June 2018: qvap, fakt added
+! forkel June 2018: reset case in  Initialize, Integrate, Update_rconst
+! 
+! 
+! 2460 2017-09-13 14:47:48Z forkel
 !
-! 
-! Variables for photolyis added
-!
-!
-!
+! forkel Sept. 2017: Variables for photolyis added
 !
 !
@@ -64,10 +82,13 @@
   PUBLIC :: nspec, nreact
   PUBLIC :: temp
+  PUBLIC :: qvap
+  PUBLIC :: fakt
   PUBLIC :: phot 
   PUBLIC :: rconst
   PUBLIC :: nvar
   PUBLIC :: nphot
+  PUBLIC :: vl_dim                     ! Public to ebable other modules to distiguish between scalar and vec
  
-  PUBLIC :: initialize, integrate, update_rconst
+  PUBLIC :: Initialize, Integrate, Update_rconst
   PUBLIC :: chem_gasphase_integrate
   PUBLIC :: initialize_kpp_ctrl
Index: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/rosenbrock_vec.f90
===================================================================
--- palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/rosenbrock_vec.f90	(revision 2718)
+++ 	(revision )
@@ -1,1161 +1,0 @@
-! NOTE: This routine is a vectorised version of the routine in
-!       rosenbrock_mz.f90 and can only be used with KP4.
-!
-SUBROUTINE Rosenbrock(Y,Tstart,Tend, &
-           AbsTol,RelTol,            &
-           RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!
-!    Solves the system y'=F(t,y) using a Rosenbrock method defined by:
-!
-!     G = 1/(H*gamma(1)) - Jac(t0,Y0)
-!     T_i = t0 + Alpha(i)*H
-!     Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
-!     G * K_i = Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j +
-!         gamma(i)*dF/dT(t0, Y0)
-!     Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
-!
-!    For details on Rosenbrock methods and their implementation consult:
-!      E. Hairer and G. Wanner
-!      "Solving ODEs II. Stiff and differential-algebraic problems".
-!      Springer series in computational mathematics, Springer-Verlag, 1996.
-!    The codes contained in the book inspired this implementation.
-!
-!    (C)  Adrian Sandu, August 2004
-!    Virginia Polytechnic Institute and State University
-!    Contact: sandu@cs.vt.edu
-!    This implementation is part of KPP - the Kinetic PreProcessor
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!
-!~~~>   INPUT ARGUMENTS:
-!
-!-     Y(NVAR)    = vector of initial conditions (at T=Tstart)
-!-    [Tstart,Tend]  = time range of integration
-!     (if Tstart>Tend the integration is performed backwards in time)
-!-    RelTol, AbsTol = user precribed accuracy
-!- SUBROUTINE  Fun( T, Y, Ydot ) = ODE function,
-!                       returns Ydot = Y' = F(T,Y)
-!- SUBROUTINE  Jac( T, Y, Jcb ) = Jacobian of the ODE function,
-!                       returns Jcb = dFun/dY
-!-    ICNTRL(1:20)    = integer inputs parameters
-!-    RCNTRL(1:20)    = real inputs parameters
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!
-!~~~>     OUTPUT ARGUMENTS:
-!
-!-    Y(NVAR)    -> vector of final states (at T->Tend)
-!-    ISTATUS(1:20)   -> integer output parameters
-!-    RSTATUS(1:20)   -> real output parameters
-!-    IERR            -> job status upon return
-!                        success (positive value) or
-!                        failure (negative value)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!
-!~~~>     INPUT PARAMETERS:
-!
-!    Note: For input parameters equal to zero the default values of the
-!       corresponding variables are used.
-!
-!    ICNTRL(1) = 1: F = F(y)   Independent of T (AUTONOMOUS)
-!              = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
-!
-!    ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
-!              = 1: AbsTol, RelTol are scalars
-!
-!    ICNTRL(3)  -> selection of a particular Rosenbrock method
-!        = 0 :  default method is Rodas3
-!        = 1 :  method is  Ros2
-!        = 2 :  method is  Ros3
-!        = 3 :  method is  Ros4
-!        = 4 :  method is  Rodas3
-!        = 5:   method is  Rodas4
-!
-!    ICNTRL(4)  -> maximum number of integration steps
-!        For ICNTRL(4)=0) the default value of 100000 is used
-!
-!    RCNTRL(1)  -> Hmin, lower bound for the integration step size
-!          It is strongly recommended to keep Hmin = ZERO
-!    RCNTRL(2)  -> Hmax, upper bound for the integration step size
-!    RCNTRL(3)  -> Hstart, starting value for the integration step size
-!
-!    RCNTRL(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
-!    RCNTRL(5)  -> FacMax, upper bound on step increase factor (default=6)
-!    RCNTRL(6)  -> FacRej, step decrease factor after multiple rejections
-!            (default=0.1)
-!    RCNTRL(7)  -> FacSafe, by which the new step is slightly smaller
-!         than the predicted value  (default=0.9)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!
-!~~~>     OUTPUT PARAMETERS:
-!
-!    Note: each call to Rosenbrock adds the current no. of fcn calls
-!      to previous value of ISTATUS(1), and similar for the other params.
-!      Set ISTATUS(1:20) = 0 before call to avoid this accumulation.
-!
-!    ISTATUS(1) = No. of function calls
-!    ISTATUS(2) = No. of jacobian calls
-!    ISTATUS(3) = No. of steps
-!    ISTATUS(4) = No. of accepted steps
-!    ISTATUS(5) = No. of rejected steps (except at the beginning)
-!    ISTATUS(6) = No. of LU decompositions
-!    ISTATUS(7) = No. of forward/backward substitutions
-!    ISTATUS(8) = No. of singular matrix decompositions
-!
-!    RSTATUS(1)  -> Texit, the time corresponding to the
-!                   computed Y upon return
-!    RSTATUS(2)  -> Hexit, last accepted step before exit
-!    For multiple restarts, use Hexit as Hstart in the following run
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-  IMPLICIT NONE
-
-!~~~>  Arguments
-   REAL(kind=dp), INTENT(INOUT) :: Y(:,:)
-   REAL(kind=dp), INTENT(IN)   :: Tstart,Tend
-   REAL(kind=dp), INTENT(IN)   :: AbsTol(NVAR),RelTol(NVAR)
-   INTEGER, INTENT(IN)    :: ICNTRL(20)
-   REAL(kind=dp), INTENT(IN)   :: RCNTRL(20)
-   INTEGER, INTENT(INOUT) :: ISTATUS(20)
-   REAL(kind=dp), INTENT(INOUT) :: RSTATUS(20)
-   INTEGER, INTENT(OUT)   :: IERR
-!~~~>  The method parameters
-   INTEGER, PARAMETER :: Smax = 6
-   INTEGER  :: Method, ros_S
-   REAL(kind=dp), DIMENSION(Smax) :: ros_M, ros_E, ros_Alpha, ros_Gamma
-   REAL(kind=dp), DIMENSION(Smax*(Smax-1)/2) :: ros_A, ros_C
-   REAL(kind=dp) :: ros_ELO
-   LOGICAL, DIMENSION(Smax) :: ros_NewF
-   CHARACTER(LEN=12) :: ros_Name
-!~~~>  Local variables
-   REAL(kind=dp) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
-   REAL(kind=dp) :: Hmin, Hmax, Hstart, Hexit
-   REAL(kind=dp) :: Texit
-   INTEGER :: i, UplimTol, Max_no_steps
-   LOGICAL :: Autonomous, VectorTol
-!~~~>   Parameters
-   REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp, ONE  = 1.0_dp
-   REAL(kind=dp), PARAMETER :: DeltaMin = 1.0E-5_dp
-
-!~~~>  Initialize statistics
-   Nfun = ISTATUS(ifun)
-   Njac = ISTATUS(ijac)
-   Nstp = ISTATUS(istp)
-   Nacc = ISTATUS(iacc)
-   Nrej = ISTATUS(irej)
-   Ndec = ISTATUS(idec)
-   Nsol = ISTATUS(isol)
-   Nsng = ISTATUS(isng)
-
-!~~~>  Autonomous or time dependent ODE. Default is time dependent.
-   Autonomous = .NOT.(ICNTRL(1) == 0)
-
-!~~~>  For Scalar tolerances (ICNTRL(2).NE.0)  the code uses AbsTol(1) and RelTol(1)
-!   For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR)
-   IF (ICNTRL(2) == 0) THEN
-      VectorTol = .TRUE.
-         UplimTol  = NVAR
-   ELSE
-      VectorTol = .FALSE.
-         UplimTol  = 1
-   END IF
-
-!~~~>  The particular Rosenbrock method chosen
-   IF (ICNTRL(3) == 0) THEN
-      Method = 4
-   ELSEIF ( (ICNTRL(3) >= 1).AND.(ICNTRL(3) <= 5) ) THEN
-      Method = ICNTRL(3)
-   ELSE
-      PRINT * , 'User-selected Rosenbrock method: ICNTRL(3)=', Method
-      !CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
-      IERRV(:) = -2
-      RETURN
-   END IF
-
-!~~~>   The maximum number of steps admitted
-   IF (ICNTRL(4) == 0) THEN
-      Max_no_steps = 100000 ! 200
-   ELSEIF (ICNTRL(4) > 0) THEN
-      Max_no_steps=ICNTRL(4)
-   ELSE
-      PRINT * ,'User-selected max no. of steps: ICNTRL(4)=',ICNTRL(4)
-      !CALL ros_ErrorMsg(-1,Tstart,ZERO,IERR)
-      IERRV(:) = -1
-      RETURN
-   END IF
-
-!~~~>  Unit roundoff (1+Roundoff>1)
-   Roundoff = WLAMCH('E')
-
-!~~~>  Lower bound on the step size: (positive value)
-   IF (RCNTRL(1) == ZERO) THEN
-      Hmin = ZERO
-   ELSEIF (RCNTRL(1) > ZERO) THEN
-      Hmin = RCNTRL(1)
-   ELSE
-      PRINT * , 'User-selected Hmin: RCNTRL(1)=', RCNTRL(1)
-      !CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
-      IERRV(:) = -3
-      RETURN
-   END IF
-!~~~>  Upper bound on the step size: (positive value)
-   IF (RCNTRL(2) == ZERO) THEN
-      Hmax = ABS(Tend-Tstart)
-   ELSEIF (RCNTRL(2) > ZERO) THEN
-      Hmax = MIN(ABS(RCNTRL(2)),ABS(Tend-Tstart))
-   ELSE
-      PRINT * , 'User-selected Hmax: RCNTRL(2)=', RCNTRL(2)
-      !CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
-      IERRV(:) = -3
-      RETURN
-   END IF
-!~~~>  Starting step size: (positive value)
-   IF (RCNTRL(3) == ZERO) THEN
-      Hstart = MAX(Hmin,DeltaMin)
-   ELSEIF (RCNTRL(3) > ZERO) THEN
-      Hstart = MIN(ABS(RCNTRL(3)),ABS(Tend-Tstart))
-   ELSE
-      PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
-      !CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
-      IERRV(:) = -3
-      RETURN
-   END IF
-!~~~>  Step size can be changed s.t.  FacMin < Hnew/Hexit < FacMax
-   IF (RCNTRL(4) == ZERO) THEN
-      FacMin = 0.2_dp
-   ELSEIF (RCNTRL(4) > ZERO) THEN
-      FacMin = RCNTRL(4)
-   ELSE
-      PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
-      !CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
-      IERRV(:) = -4
-      RETURN
-   END IF
-   IF (RCNTRL(5) == ZERO) THEN
-      FacMax = 6.0_dp
-   ELSEIF (RCNTRL(5) > ZERO) THEN
-      FacMax = RCNTRL(5)
-   ELSE
-      PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
-      !CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
-      IERRV(:) = -4
-      RETURN
-   END IF
-!~~~>   FacRej: Factor to decrease step after 2 succesive rejections
-   IF (RCNTRL(6) == ZERO) THEN
-      FacRej = 0.1_dp
-   ELSEIF (RCNTRL(6) > ZERO) THEN
-      FacRej = RCNTRL(6)
-   ELSE
-      PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
-      !CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
-      IERRV(:) = -4
-      RETURN
-   END IF
-!~~~>   FacSafe: Safety Factor in the computation of new step size
-   IF (RCNTRL(7) == ZERO) THEN
-      FacSafe = 0.9_dp
-   ELSEIF (RCNTRL(7) > ZERO) THEN
-      FacSafe = RCNTRL(7)
-   ELSE
-      PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
-      !CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
-      IERRV(:) = -4
-      RETURN
-   END IF
-!~~~>  Check if tolerances are reasonable
-    DO i=1,UplimTol
-      IF ( (AbsTol(i) <= ZERO) .OR. (RelTol(i) <= 10.0_dp*Roundoff) &
-         .OR. (RelTol(i) >= 1.0_dp) ) THEN
-        PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
-        PRINT * , ' RelTol(',i,') = ',RelTol(i)
-        !CALL ros_ErrorMsg(-5,Tstart,ZERO,IERR)
-        IERRV(:) = -5
-        RETURN
-      END IF
-    END DO
-
-
-!~~~>   Initialize the particular Rosenbrock method
-   SELECT CASE (Method)
-     CASE (1)
-       CALL Ros2(ros_S, ros_A, ros_C, ros_M, ros_E,   &
-          ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
-     CASE (2)
-       CALL Ros3(ros_S, ros_A, ros_C, ros_M, ros_E,   &
-          ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
-     CASE (3)
-       CALL Ros4(ros_S, ros_A, ros_C, ros_M, ros_E,   &
-          ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
-     CASE (4)
-       CALL Rodas3(ros_S, ros_A, ros_C, ros_M, ros_E, &
-          ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
-     CASE (5)
-       CALL Rodas4(ros_S, ros_A, ros_C, ros_M, ros_E, &
-          ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
-     CASE DEFAULT
-       PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=', Method
-       !CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
-       IERRV(:) = -2
-       RETURN
-   END SELECT
-
-!~~~>  CALL Rosenbrock method
-   CALL ros_Integrator(Y,Tstart,Tend,Texit,      &
-        AbsTol, RelTol,                          &
-!  Rosenbrock method coefficients
-        ros_S, ros_M, ros_E, ros_A, ros_C,       &
-        ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
-!  Integration parameters
-        Autonomous, VectorTol, Max_no_steps,     &
-        Roundoff, Hmin, Hmax, Hstart, Hexit,     &
-        FacMin, FacMax, FacRej, FacSafe,         &
-!  Error indicator
-        IERR)
-
-
-!~~~>  Collect run statistics
-   ISTATUS(ifun) = Nfun
-   ISTATUS(ijac) = Njac
-   ISTATUS(istp) = Nstp
-   ISTATUS(iacc) = Nacc
-   ISTATUS(irej) = Nrej
-   ISTATUS(idec) = Ndec
-   ISTATUS(isol) = Nsol
-   ISTATUS(isng) = Nsng
-!~~~> Last T and H
-   RSTATUS(itexit) = Texit
-   RSTATUS(ihexit) = Hexit
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-CONTAINS !  SUBROUTINES internal to Rosenbrock
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- SUBROUTINE ros_ErrorMsg(Code,T,H,IERR)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!    Handles all error messages
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-   REAL(kind=dp), INTENT(IN) :: T, H
-   INTEGER, INTENT(IN)  :: Code
-   INTEGER, INTENT(INOUT) :: IERR
-
-   IERR = Code
-   PRINT * , &
-     'Forced exit from Rosenbrock due to the following error:'
-   IF ((Code>=-8).AND.(Code<=-1)) THEN
-     PRINT *, IERR_NAMES(Code)
-   ELSE
-     PRINT *, 'Unknown Error code: ', Code
-   ENDIF
-
-   PRINT *, "T=", T, "and H=", H
-
- END SUBROUTINE ros_ErrorMsg
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- SUBROUTINE ros_Integrator (Y, Tstart, Tend, Tout,     &
-        AbsTol, RelTol,                          &
-!~~~> Rosenbrock method coefficients
-        ros_S, ros_M, ros_E, ros_A, ros_C,       &
-        ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
-!~~~> Integration parameters
-        Autonomous, VectorTol, Max_no_steps,     &
-        Roundoff, Hmin, Hmax, Hstart, Hexit,     &
-        FacMin, FacMax, FacRej, FacSafe,         &
-!~~~> Error indicator
-        IERR )
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!   Template for the implementation of a generic Rosenbrock method
-!      defined by ros_S (no of stages)
-!      and its coefficients ros_{A,C,M,E,Alpha,Gamma}
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-  USE messy_main_tools_kp4_compress, ONLY: kco_initialize, kco_compress, &
-                                           kco_finalize, cell_done
-  IMPLICIT NONE
-
-!~~~> Input: the initial condition at Tstart; Output: the solution at T
-   REAL(kind=dp), INTENT(INOUT) :: Y(:,:)
-!~~~> Input: integration interval
-   REAL(kind=dp), INTENT(IN) :: Tstart,Tend
-!~~~> Output: time at which the solution is returned (T=Tend if success)
-   REAL(kind=dp), INTENT(OUT) ::  Tout
-!~~~> Input: tolerances
-   REAL(kind=dp), INTENT(IN) ::  AbsTol(NVAR), RelTol(NVAR)
-!~~~> Input: The Rosenbrock method parameters
-   INTEGER, INTENT(IN) ::  ros_S
-   REAL(kind=dp), INTENT(IN) :: ros_M(ros_S), ros_E(ros_S),  &
-       ros_Alpha(ros_S), ros_A(ros_S*(ros_S-1)/2), &
-       ros_Gamma(ros_S), ros_C(ros_S*(ros_S-1)/2), ros_ELO
-   LOGICAL, INTENT(IN) :: ros_NewF(ros_S)
-!~~~> Input: integration parameters
-   LOGICAL, INTENT(IN) :: Autonomous, VectorTol
-   REAL(kind=dp), INTENT(IN) :: Hstart, Hmin, Hmax
-   INTEGER, INTENT(IN) :: Max_no_steps
-   REAL(kind=dp), INTENT(IN) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
-!~~~> Output: last accepted step
-   REAL(kind=dp), INTENT(OUT) :: Hexit
-!~~~> Output: Error indicator
-   INTEGER, INTENT(OUT) :: IERR
-! ~~~~ Local variables
-   REAL(kind=dp) :: Ynew(VL_DIM,NVAR), Fcn0(VL_DIM,NVAR), Fcn(VL_DIM,NVAR)
-   REAL(kind=dp) :: K(VL_DIM,NVAR*ros_S)
-   REAL(kind=dp) :: Jac0(VL_DIM,LU_NONZERO), Ghimj(VL_DIM,LU_NONZERO)
-   REAL(kind=dp) :: Yerr(VL_DIM,NVAR)
-
-!~~~>  Local parameters
-   INTEGER :: ioffset, i, j, istage
-   REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp, ONE  = 1.0_dp
-   REAL(kind=dp), PARAMETER :: DeltaMin = 1.0E-5_dp
-!~~~>  Locally called functions
-!    REAL(kind=dp) WLAMCH
-!    EXTERNAL WLAMCH
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-   real(kind=dp)                  :: deltat
-!kk   logical,parameter              :: fixed_step=.true.
-   
-   integer,save             :: icount=0
-
-!  vector compression algorithm
-
-   REAL(kind=dp), dimension(VL_dim)      :: T
-   REAL(kind=dp), dimension(VL_dim)      :: H, Hnew, HC, HC2, HG, Fac, Tau, Err
-   LOGICAL, dimension(VL_dim)            :: RejectLastH, RejectMoreH
-   LOGICAL, dimension(VL_dim)            :: Accept_step
-   integer                               :: vl_save
-
-   REAL(kind=dp), dimension(size(FIX,1),size(FIX,2))      :: FIX_lo
-
-!~~~>  Initial preparations
-   T(1:VL)  = Tstart
-   Hexit    = 0.0_dp
-   H(1:VL)  = MIN(Hstart,Hmax)
-   Tout     = Tend
-   icount   = icount+1
-   where (ABS(H(1:VL)) <= 10.0_dp*Roundoff) H(1:VL) = DeltaMin
-
-   deltat=Tend-Tstart
-
-   RejectLastH(1:VL)=.FALSE.
-   RejectMoreH(1:VL)=.FALSE.
-
-!~~~> Time loop begins below
-
-   Nstp       = 0
-   IERRV(1:VL) = 1  ! successful integration
-   vl_save          = VL
-   FIX_lo(1:VL,:)   = FIX(1:VL,:)
-   Kacc(1:VL)       = 0
-   Krej(1:VL)       = 0
-
-!  intialize compression algorithm
-
-   IF (.not. l_fixed_step) then
-      call kco_initialize (VL, y, IERRV, Kacc, Krej)
-   end if
-
-TimeLoop: DO 
-
-   IF ( Nstp > Max_no_steps ) THEN  ! Too many steps
-      !CALL ros_ErrorMsg(-6,T(1),H(1),IERR)
-      IERRV(1:VL) = -6
-!      call abort()                                   !Vorlaeufig
-      RETURN
-   END IF
-   
-   Where ( (T(1:VL)+0.1_dp*H(1:VL)) == T(1:VL))      ! Step size too small
-      IERRV(1:VL) = -7
-   endwhere
-
-!!$   if( any(error_flag))  then
-!!$      CALL ros_ErrorMsg(-7,T(1),H(1),IERR)
-!!$      call abort()                                   !Vorlaeufig
-!!$      RETURN
-!!$   end if
-
-!~~~>  Limit H if necessary to avoid going beyond Tend
-!kk   Hexit = H
-   H(1:VL) = MIN(H(1:VL),ABS(Tend-T(1:VL)))
-
-!~~~>   Compute the function at current time
-   CALL Fun( Y, FIX_lo, RCONST, Fcn0)
-
-
-!~~~>   Compute the Jacobian at current time
-   CALL Jac_SP( Y, FIX_lo, RCONST, Jac0)
-
-
-!~~~>  Repeat step calculation until current step accepted
-
-   CALL ros_PrepareMatrix(H,ros_Gamma(1), Jac0,Ghimj)
-
-!~~~>   Compute the stages
-Stage: DO istage = 1, ros_S
-
-      ! Current istage offset. Current istage vector is K(ioffset+1:ioffset+NVAR)
-       ioffset = NVAR*(istage-1)
-
-      ! For the 1st istage the function has been computed previously
-       IF ( istage == 1 ) THEN
-         Fcn(1:VL,1:NVAR) = Fcn0(1:VL,1:NVAR) 
-      ! istage>1 and a new function evaluation is needed at the current istage
-       ELSEIF ( ros_NewF(istage) ) THEN
-         select case (istage-1)
-          case (1)
-           Ynew(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
-           Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) +  K(1:VL,NVAR*(1- 1)+ 1:NVAR*1) &
-                        *ros_A((istage-1)*(istage-2)/2+ 1)
-          case (2)
-           Ynew(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
-           Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) +  K(1:VL,NVAR*(2- 1)+ 1:NVAR*2) &
-                        *ros_A((istage-1)*(istage-2)/2+ 2)
-           Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) +  K(1:VL,NVAR*(2- 1)+ 1:NVAR*2) &
-                        *ros_A((istage-1)*(istage-2)/2+ 2)
-          case DEFAULT
-           Ynew(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
-           DO j = 1, istage-1
-             Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) +  K(1:VL,NVAR*(j- 1)+ 1:NVAR*j) &
-                          *ros_A((istage-1)*(istage-2)/2+ j)
-           END DO
-         end select
-
-         CALL Fun( Ynew, FIX_lo, RCONST, Fcn)
-       END IF ! if istage == 1 elseif ros_NewF(istage)
-
-       select case (istage-1)
-        case (1)
-         K(1:VL,ioffset+ 1:ioffset+ NVAR) = Fcn(1:VL,1:NVAR)
-         j = 1
-         HC(1:VL) = ros_C((istage-1)*(istage-2)/2+ j)/(H(1:VL))
-         do i=1,NVAR
-           K(1:VL,ioffset+ i) = K(1:VL,ioffset+ i) +  K(1:VL,NVAR*(j- 1)+ i) *HC(1:VL)
-         end do
-        case (2)
-         K(1:VL,ioffset+ 1:ioffset+ NVAR) = Fcn(1:VL,1:NVAR)
-         HC(1:VL)  = ros_C((istage-1)*(istage-2)/2+ 1)/(H(1:VL))
-         HC2(1:VL) = ros_C((istage-1)*(istage-2)/2+ 2)/(H(1:VL))
-         do i=1,NVAR
-           K(1:VL,ioffset+ i) = K(1:VL,ioffset+ i) +  K(1:VL,NVAR*(1- 1)+ i) *HC(1:VL)
-           K(1:VL,ioffset+ i) = K(1:VL,ioffset+ i) +  K(1:VL,NVAR*(2- 1)+ i) *HC2(1:VL)
-         end do
-        case DEFAULT
-         K(1:VL,ioffset+ 1:ioffset+ NVAR) = Fcn(1:VL,1:NVAR)
-         DO j = 1, istage-1
-           HC(1:VL) = ros_C((istage-1)*(istage-2)/2+ j)/(H(1:VL))
-           do i=1,NVAR
-             K(1:VL,ioffset+ i) = K(1:VL,ioffset+ i) +  K(1:VL,NVAR*(j- 1)+ i) *HC(1:VL)
-           end do
-         END DO
-       end select
-
-       CALL ros_Solve(Ghimj, K(:,ioffset+1:ioffset+NVAR))
-
-   END DO Stage
-
-
-   if(ros_S == 3)   then
-     YNEW(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
-     Yerr(1:VL,1:NVAR) = ZERO
-
-!~~~>  Compute the new solution
-     Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR*(1-1)+1:NVAR*1) * ros_m(1)
-     Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR*(2-1)+1:NVAR*2) * ros_m(2)
-     Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR*(3-1)+1:NVAR*3) * ros_m(3)
-
-!~~~>  Compute the error estimation
-
-     Yerr(1:VL,1:NVAR) = Yerr(1:VL,1:NVAR) + K(1:VL,NVAR*(1-1)+1:NVAR*1) * ros_E(1)
-     Yerr(1:VL,1:NVAR) = Yerr(1:VL,1:NVAR) + K(1:VL,NVAR*(2-1)+1:NVAR*2) * ros_E(2)
-     Yerr(1:VL,1:NVAR) = Yerr(1:VL,1:NVAR) + K(1:VL,NVAR*(3-1)+1:NVAR*3) * ros_E(3)
-   else
-     YNEW(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
-     Yerr(1:VL,1:NVAR) = ZERO
-     DO j=1,ros_S
-!~~~>  Compute the new solution
-         Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + &
-                        K(1:VL,NVAR*(j-1)+1:NVAR*j) * ros_m(j)
-
-!~~~>  Compute the error estimation
-
-       Yerr(1:VL,1:NVAR) = Yerr(1:VL,1:NVAR) + &
-                    K(1:VL,NVAR*(j-1)+1:NVAR*j) * ros_E(j)
-     END DO
-   end if
-
-
-!~~~>  Check the error magnitude and adjust step size
-   Nstp = Nstp+1
-
-   if(l_fixed_step)   then
-     Err  = 0.0
-     H    = t_steps(Nstp)*deltat
-     Hnew = t_steps(Nstp)*deltat
-   else
-     Err(1:VL) = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol )
-!~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
-     where(Err(1:VL) >= 10.0_dp*Roundoff)
-       Fac(1:VL)  = MIN(FacMax,MAX(FacMin,FacSafe/Err(1:VL)**(ONE/ros_ELO)))
-     elsewhere
-       Fac(1:VL)  = 1000.                        !Force large time step
-     endwhere
-     Hnew(1:VL) = H(1:VL)*Fac(1:VL)
-   end if
-
-   Accept_step(1:VL) = (Err(1:VL) <= ONE).OR.(H(1:VL) <= Hmin)
-   where (Accept_step(1:VL))
-      Kacc(1:VL) = Kacc(1:VL)+1
-      T(1:VL) = T(1:VL) + H(1:VL)
-      Hnew(1:VL) = MAX(Hmin,MIN(Hnew(1:VL),Hmax))
-   endwhere
-   do i=1,NVAR
-     where (Accept_step(1:VL))
-       Y(1:VL,i) = MAX(Ynew(1:VL,i),ZERO)
-     endwhere
-   end do
-   where (Accept_step(1:VL) .and. RejectLastH(1:VL))
-      Hnew(1:VL) = MIN(Hnew(1:VL),H(1:VL))
-   endwhere
-   where (Accept_step(1:VL))
-      RejectLastH(1:VL) = .FALSE.
-      RejectMoreH(1:VL) = .FALSE.
-      H(1:VL) = Hnew(1:VL)
-   endwhere
-!kk   ELSE           !~~~> Reject step
-   where (.NOT. Accept_step(1:VL) .and. RejectMoreH(1:VL))
-      Hnew(1:VL) = H(1:VL)*FacRej
-   endwhere
-   where (.NOT. Accept_step(1:VL))
-      RejectMoreH(1:VL) = RejectLastH(1:VL)
-      RejectLastH(1:VL) = .TRUE.
-      H(1:VL) = Hnew(1:VL)
-   endwhere
-   where ( .NOT. Accept_step(1:VL) .and. Kacc(1:VL) >= 1)
-      Krej(1:VL) = Krej(1:VL)+1
-   endwhere
-
-
-   if(l_fixed_step)   then
-      if(Nstp >= Nfsteps) EXIT
-   else
-     cell_done(1:VL) = ( .NOT. (T(1:VL)-Tend)+Roundoff <= ZERO  )  .OR. &
-          (IERRV(1:VL) < 0) ! mz_pj_20070516
-     call kco_compress (VL, T, H, Hnew, IERRV, y, RCONST, FIX_lo,         &
-                                     RejectLastH, RejectMoreH, Kacc, Krej)
-     if(VL == 0)   EXIT
-   end if
-
-!kk   write(0,*) '  aaa ',icount,Nstp,count(cell_done),count(Accept_step),VL,minval(T)
-   END DO TimeLoop
-
-!~~~> Succesful exit
-   IERR = 1  !~~~> The integration was successful
-
-   IF (.not. l_fixed_step) then
-     call kco_finalize (y, IERRV, Kacc, Krej)
-   end if
-
-   VL = VL_save
-
-!kk   write(0,*) 'bbb ',icount,Nstp,maxval(Kacc),maxval(Krej)
-!kk   if(icount >= 10)  call abort()
-
-  END SUBROUTINE ros_Integrator
-
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  FUNCTION ros_ErrorNorm ( Y, Ynew, Yerr, &
-               AbsTol, RelTol, VectorTol )
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!~~~> Computes the "scaled norm" of the error vector Yerr
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-   IMPLICIT NONE
-
-! Input arguments
-   REAL(kind=dp), INTENT(IN) :: Y(:,:), Ynew(:,:), &
-          Yerr(:,:), AbsTol(NVAR), RelTol(NVAR)
-   LOGICAL, INTENT(IN) ::  VectorTol
-! Local variables
-   REAL(kind=dp), dimension(VL)       :: Err, Scale, Ymax
-   INTEGER  :: i
-   REAL(kind=dp), PARAMETER           :: ZERO = 0.0_dp
-   REAL(kind=dp), dimension(VL)       :: ros_ErrorNorm
-
-   Err = ZERO
-   DO i=1,NVAR
-     Ymax = MAX( ABS(Y(1:VL,i)),ABS(Ynew(1:VL,i)))
-     IF (VectorTol) THEN
-       Scale = AbsTol(i)+ RelTol(i)*Ymax
-     ELSE
-       Scale = AbsTol(1)+ RelTol(1)*Ymax
-     END IF
-     Err = Err+ (Yerr(1:VL,i)/Scale)**2
-   END DO
-   Err  = SQRT(Err/NVAR)
-
-   ros_ErrorNorm = Err
-
-  END FUNCTION ros_ErrorNorm
-
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-  SUBROUTINE ros_PrepareMatrix ( H, gam, Jac0, Ghimj )
-
-! --- --- --- --- --- --- --- --- --- --- --- --- ---
-!  Prepares the LHS matrix for stage calculations
-!  1.  Construct Ghimj = 1/(H*ham) - Jac0
-!      "(Gamma H) Inverse Minus Jacobian"
-!  2.  Repeat LU decomposition of Ghimj until successful.
-!       -half the step size if LU decomposition fails and retry
-!       -exit after 5 consecutive fails
-! --- --- --- --- --- --- --- --- --- --- --- --- ---
-   IMPLICIT NONE
-
-!~~~> Input arguments
-   REAL(kind=dp), INTENT(IN) ::  Jac0(:,:)
-   REAL(kind=dp), INTENT(IN) ::  gam
-!~~~> Output arguments
-   REAL(kind=dp), INTENT(OUT) :: Ghimj(:,:)
-!~~~> Inout arguments
-   REAL(kind=dp), INTENT(INOUT) :: H(:)   ! step size is decreased when LU fails
-!~~~> Local variables
-   INTEGER  :: i, ising, Nconsecutive
-   REAL(kind=dp), dimension(size(H)) :: ghinv
-   REAL(kind=dp), PARAMETER :: ONE  = 1.0_dp, HALF = 0.5_dp
-
-   Nconsecutive = 0
-
-
-!~~~>    Construct Ghimj = 1/(H*gam) - Jac0
-   Ghimj(1:VL,1:LU_NONZERO) = -ONE * JAC0(1:VL,1:LU_NONZERO)
-   ghinv(1:VL) = ONE/(H(1:VL)*gam)
-   DO i=1,NVAR
-     Ghimj(1:VL,LU_DIAG(i)) = Ghimj(1:VL,LU_DIAG(i))+ghinv(1:VL)
-   END DO
-
-!~~~>    Compute LU decomposition
-!CDIR noiexpand
-   CALL KppDecomp( Ghimj, ising )
-
-   Ndec = Ndec + VL
-
-  END SUBROUTINE ros_PrepareMatrix
-
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  SUBROUTINE ros_Solve( A, b )
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!  Template for the forward/backward substitution (using pre-computed LU decomposition)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-   IMPLICIT NONE
-!~~~> Input variables
-   REAL(kind=dp), INTENT(IN)         :: A(:,:)
-!~~~> InOut variables
-   REAL(kind=dp), INTENT(INOUT)      :: b(:,:)
-
-   CALL KppSolve( A, b )
-
-   Nsol = Nsol+VL
-
-  END SUBROUTINE ros_Solve
-
-
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  SUBROUTINE Ros2 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
-            ros_Gamma,ros_NewF,ros_ELO,ros_Name)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! --- AN L-STABLE METHOD, 2 stages, order 2
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-  IMPLICIT NONE
-
-   INTEGER, PARAMETER :: S=2
-   INTEGER, INTENT(OUT) ::  ros_S
-   REAL(kind=dp), DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
-   REAL(kind=dp), DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
-   REAL(kind=dp), INTENT(OUT) :: ros_ELO
-   LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
-   CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
-   DOUBLE PRECISION g
-
-    g = 1.0_dp + 1.0_dp/SQRT(2.0_dp)
-
-!~~~> Name of the method
-    ros_Name = 'ROS-2'
-!~~~> Number of stages
-    ros_S = S
-
-!~~~> The coefficient matrices A and C are strictly lower triangular.
-!   The lower triangular (subdiagonal) elements are stored in row-wise order:
-!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
-!   The general mapping formula is:
-!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
-!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
-
-    ros_A(1) = (1.0_dp)/g
-    ros_C(1) = (-2.0_dp)/g
-!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
-!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
-    ros_NewF(1) = .TRUE.
-    ros_NewF(2) = .TRUE.
-!~~~> M_i = Coefficients for new step solution
-    ros_M(1)= (3.0_dp)/(2.0_dp*g)
-    ros_M(2)= (1.0_dp)/(2.0_dp*g)
-! E_i = Coefficients for error estimator
-    ros_E(1) = 1.0_dp/(2.0_dp*g)
-    ros_E(2) = 1.0_dp/(2.0_dp*g)
-!~~~> ros_ELO = estimator of local order - the minimum between the
-!    main and the embedded scheme orders plus one
-    ros_ELO = 2.0_dp
-!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
-    ros_Alpha(1) = 0.0_dp
-    ros_Alpha(2) = 1.0_dp
-!~~~> Gamma_i = \sum_j  gamma_{i,j}
-    ros_Gamma(1) = g
-    ros_Gamma(2) =-g
-
- END SUBROUTINE Ros2
-
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  SUBROUTINE Ros3 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
-           ros_Gamma,ros_NewF,ros_ELO,ros_Name)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-  IMPLICIT NONE
-
-   INTEGER, PARAMETER :: S=3
-   INTEGER, INTENT(OUT) ::  ros_S
-   REAL(kind=dp), DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
-   REAL(kind=dp), DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
-   REAL(kind=dp), INTENT(OUT) :: ros_ELO
-   LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
-   CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
-
-!~~~> Name of the method
-   ros_Name = 'ROS-3'
-!~~~> Number of stages
-   ros_S = S
-
-!~~~> The coefficient matrices A and C are strictly lower triangular.
-!   The lower triangular (subdiagonal) elements are stored in row-wise order:
-!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
-!   The general mapping formula is:
-!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
-!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
-
-   ros_A(1)= 1.0_dp
-   ros_A(2)= 1.0_dp
-   ros_A(3)= 0.0_dp
-
-   ros_C(1) = -0.10156171083877702091975600115545E+01_dp
-   ros_C(2) =  0.40759956452537699824805835358067E+01_dp
-   ros_C(3) =  0.92076794298330791242156818474003E+01_dp
-!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
-!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
-   ros_NewF(1) = .TRUE.
-   ros_NewF(2) = .TRUE.
-   ros_NewF(3) = .FALSE.
-!~~~> M_i = Coefficients for new step solution
-   ros_M(1) =  0.1E+01_dp
-   ros_M(2) =  0.61697947043828245592553615689730E+01_dp
-   ros_M(3) = -0.42772256543218573326238373806514E+00_dp
-! E_i = Coefficients for error estimator
-   ros_E(1) =  0.5E+00_dp
-   ros_E(2) = -0.29079558716805469821718236208017E+01_dp
-   ros_E(3) =  0.22354069897811569627360909276199E+00_dp
-!~~~> ros_ELO = estimator of local order - the minimum between the
-!    main and the embedded scheme orders plus 1
-   ros_ELO = 3.0_dp
-!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
-   ros_Alpha(1)= 0.0E+00_dp
-   ros_Alpha(2)= 0.43586652150845899941601945119356E+00_dp
-   ros_Alpha(3)= 0.43586652150845899941601945119356E+00_dp
-!~~~> Gamma_i = \sum_j  gamma_{i,j}
-   ros_Gamma(1)= 0.43586652150845899941601945119356E+00_dp
-   ros_Gamma(2)= 0.24291996454816804366592249683314E+00_dp
-   ros_Gamma(3)= 0.21851380027664058511513169485832E+01_dp
-
-  END SUBROUTINE Ros3
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  SUBROUTINE Ros4 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
-           ros_Gamma,ros_NewF,ros_ELO,ros_Name)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!     L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
-!     L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
-!
-!      E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
-!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
-!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
-!      SPRINGER-VERLAG (1990)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-  IMPLICIT NONE
-
-   INTEGER, PARAMETER :: S=4
-   INTEGER, INTENT(OUT) ::  ros_S
-   REAL(kind=dp), DIMENSION(4), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
-   REAL(kind=dp), DIMENSION(6), INTENT(OUT) :: ros_A, ros_C
-   REAL(kind=dp), INTENT(OUT) :: ros_ELO
-   LOGICAL, DIMENSION(4), INTENT(OUT) :: ros_NewF
-   CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
-   DOUBLE PRECISION g
-
-!~~~> Name of the method
-   ros_Name = 'ROS-4'
-!~~~> Number of stages
-   ros_S = S
-
-!~~~> The coefficient matrices A and C are strictly lower triangular.
-!   The lower triangular (subdiagonal) elements are stored in row-wise order:
-!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
-!   The general mapping formula is:
-!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
-!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
-
-   ros_A(1) = 0.2000000000000000E+01_dp
-   ros_A(2) = 0.1867943637803922E+01_dp
-   ros_A(3) = 0.2344449711399156E+00_dp
-   ros_A(4) = ros_A(2)
-   ros_A(5) = ros_A(3)
-   ros_A(6) = 0.0_dp
-
-   ros_C(1) =-0.7137615036412310E+01_dp
-   ros_C(2) = 0.2580708087951457E+01_dp
-   ros_C(3) = 0.6515950076447975E+00_dp
-   ros_C(4) =-0.2137148994382534E+01_dp
-   ros_C(5) =-0.3214669691237626E+00_dp
-   ros_C(6) =-0.6949742501781779E+00_dp
-!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
-!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
-   ros_NewF(1)  = .TRUE.
-   ros_NewF(2)  = .TRUE.
-   ros_NewF(3)  = .TRUE.
-   ros_NewF(4)  = .FALSE.
-!~~~> M_i = Coefficients for new step solution
-   ros_M(1) = 0.2255570073418735E+01_dp
-   ros_M(2) = 0.2870493262186792E+00_dp
-   ros_M(3) = 0.4353179431840180E+00_dp
-   ros_M(4) = 0.1093502252409163E+01_dp
-!~~~> E_i  = Coefficients for error estimator
-   ros_E(1) =-0.2815431932141155E+00_dp
-   ros_E(2) =-0.7276199124938920E-01_dp
-   ros_E(3) =-0.1082196201495311E+00_dp
-   ros_E(4) =-0.1093502252409163E+01_dp
-!~~~> ros_ELO  = estimator of local order - the minimum between the
-!    main and the embedded scheme orders plus 1
-   ros_ELO  = 4.0_dp
-!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
-   ros_Alpha(1) = 0.0_dp
-   ros_Alpha(2) = 0.1145640000000000E+01_dp
-   ros_Alpha(3) = 0.6552168638155900E+00_dp
-   ros_Alpha(4) = ros_Alpha(3)
-!~~~> Gamma_i = \sum_j  gamma_{i,j}
-   ros_Gamma(1) = 0.5728200000000000E+00_dp
-   ros_Gamma(2) =-0.1769193891319233E+01_dp
-   ros_Gamma(3) = 0.7592633437920482E+00_dp
-   ros_Gamma(4) =-0.1049021087100450E+00_dp
-
-  END SUBROUTINE Ros4
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  SUBROUTINE Rodas3 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
-            ros_Gamma,ros_NewF,ros_ELO,ros_Name)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! --- A STIFFLY-STABLE METHOD, 4 stages, order 3
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-  IMPLICIT NONE
-
-   INTEGER, PARAMETER :: S=4
-   INTEGER, INTENT(OUT) ::  ros_S
-   REAL(kind=dp), DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
-   REAL(kind=dp), DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
-   REAL(kind=dp), INTENT(OUT) :: ros_ELO
-   LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
-   CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
-   DOUBLE PRECISION g
-
-!~~~> Name of the method
-   ros_Name = 'RODAS-3'
-!~~~> Number of stages
-   ros_S = S
-
-!~~~> The coefficient matrices A and C are strictly lower triangular.
-!   The lower triangular (subdiagonal) elements are stored in row-wise order:
-!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
-!   The general mapping formula is:
-!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
-!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
-
-   ros_A(1) = 0.0E+00_dp
-   ros_A(2) = 2.0E+00_dp
-   ros_A(3) = 0.0E+00_dp
-   ros_A(4) = 2.0E+00_dp
-   ros_A(5) = 0.0E+00_dp
-   ros_A(6) = 1.0E+00_dp
-
-   ros_C(1) = 4.0E+00_dp
-   ros_C(2) = 1.0E+00_dp
-   ros_C(3) =-1.0E+00_dp
-   ros_C(4) = 1.0E+00_dp
-   ros_C(5) =-1.0E+00_dp
-   ros_C(6) =-(8.0E+00_dp/3.0E+00_dp)
-
-!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
-!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
-   ros_NewF(1)  = .TRUE.
-   ros_NewF(2)  = .FALSE.
-   ros_NewF(3)  = .TRUE.
-   ros_NewF(4)  = .TRUE.
-!~~~> M_i = Coefficients for new step solution
-   ros_M(1) = 2.0E+00_dp
-   ros_M(2) = 0.0E+00_dp
-   ros_M(3) = 1.0E+00_dp
-   ros_M(4) = 1.0E+00_dp
-!~~~> E_i  = Coefficients for error estimator
-   ros_E(1) = 0.0E+00_dp
-   ros_E(2) = 0.0E+00_dp
-   ros_E(3) = 0.0E+00_dp
-   ros_E(4) = 1.0E+00_dp
-!~~~> ros_ELO  = estimator of local order - the minimum between the
-!    main and the embedded scheme orders plus 1
-   ros_ELO  = 3.0E+00_dp
-!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
-   ros_Alpha(1) = 0.0E+00_dp
-   ros_Alpha(2) = 0.0E+00_dp
-   ros_Alpha(3) = 1.0E+00_dp
-   ros_Alpha(4) = 1.0E+00_dp
-!~~~> Gamma_i = \sum_j  gamma_{i,j}
-   ros_Gamma(1) = 0.5E+00_dp
-   ros_Gamma(2) = 1.5E+00_dp
-   ros_Gamma(3) = 0.0E+00_dp
-   ros_Gamma(4) = 0.0E+00_dp
-
-  END SUBROUTINE Rodas3
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  SUBROUTINE Rodas4 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
-             ros_Gamma,ros_NewF,ros_ELO,ros_Name)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!     STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
-!
-!      E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
-!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
-!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
-!      SPRINGER-VERLAG (1996)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-  IMPLICIT NONE
-
-   INTEGER, PARAMETER :: S=6
-   INTEGER, INTENT(OUT) ::  ros_S
-   REAL(kind=dp), DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
-   REAL(kind=dp), DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
-   REAL(kind=dp), INTENT(OUT) :: ros_ELO
-   LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
-   CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
-   DOUBLE PRECISION g
-
-!~~~> Name of the method
-    ros_Name = 'RODAS-4'
-!~~~> Number of stages
-    ros_S = 6
-
-!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
-    ros_Alpha(1) = 0.000_dp
-    ros_Alpha(2) = 0.386_dp
-    ros_Alpha(3) = 0.210_dp
-    ros_Alpha(4) = 0.630_dp
-    ros_Alpha(5) = 1.000_dp
-    ros_Alpha(6) = 1.000_dp
-
-!~~~> Gamma_i = \sum_j  gamma_{i,j}
-    ros_Gamma(1) = 0.2500000000000000E+00_dp
-    ros_Gamma(2) =-0.1043000000000000E+00_dp
-    ros_Gamma(3) = 0.1035000000000000E+00_dp
-    ros_Gamma(4) =-0.3620000000000023E-01_dp
-    ros_Gamma(5) = 0.0_dp
-    ros_Gamma(6) = 0.0_dp
-
-!~~~> The coefficient matrices A and C are strictly lower triangular.
-!   The lower triangular (subdiagonal) elements are stored in row-wise order:
-!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
-!   The general mapping formula is:  A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
-!                  C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
-
-    ros_A(1) = 0.1544000000000000E+01_dp
-    ros_A(2) = 0.9466785280815826E+00_dp
-    ros_A(3) = 0.2557011698983284E+00_dp
-    ros_A(4) = 0.3314825187068521E+01_dp
-    ros_A(5) = 0.2896124015972201E+01_dp
-    ros_A(6) = 0.9986419139977817E+00_dp
-    ros_A(7) = 0.1221224509226641E+01_dp
-    ros_A(8) = 0.6019134481288629E+01_dp
-    ros_A(9) = 0.1253708332932087E+02_dp
-    ros_A(10) =-0.6878860361058950E+00_dp
-    ros_A(11) = ros_A(7)
-    ros_A(12) = ros_A(8)
-    ros_A(13) = ros_A(9)
-    ros_A(14) = ros_A(10)
-    ros_A(15) = 1.0E+00_dp
-
-    ros_C(1) =-0.5668800000000000E+01_dp
-    ros_C(2) =-0.2430093356833875E+01_dp
-    ros_C(3) =-0.2063599157091915E+00_dp
-    ros_C(4) =-0.1073529058151375E+00_dp
-    ros_C(5) =-0.9594562251023355E+01_dp
-    ros_C(6) =-0.2047028614809616E+02_dp
-    ros_C(7) = 0.7496443313967647E+01_dp
-    ros_C(8) =-0.1024680431464352E+02_dp
-    ros_C(9) =-0.3399990352819905E+02_dp
-    ros_C(10) = 0.1170890893206160E+02_dp
-    ros_C(11) = 0.8083246795921522E+01_dp
-    ros_C(12) =-0.7981132988064893E+01_dp
-    ros_C(13) =-0.3152159432874371E+02_dp
-    ros_C(14) = 0.1631930543123136E+02_dp
-    ros_C(15) =-0.6058818238834054E+01_dp
-
-!~~~> M_i = Coefficients for new step solution
-    ros_M(1) = ros_A(7)
-    ros_M(2) = ros_A(8)
-    ros_M(3) = ros_A(9)
-    ros_M(4) = ros_A(10)
-    ros_M(5) = 1.0E+00_dp
-    ros_M(6) = 1.0E+00_dp
-
-!~~~> E_i  = Coefficients for error estimator
-    ros_E(1) = 0.0E+00_dp
-    ros_E(2) = 0.0E+00_dp
-    ros_E(3) = 0.0E+00_dp
-    ros_E(4) = 0.0E+00_dp
-    ros_E(5) = 0.0E+00_dp
-    ros_E(6) = 1.0E+00_dp
-
-!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
-!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
-    ros_NewF(1) = .TRUE.
-    ros_NewF(2) = .TRUE.
-    ros_NewF(3) = .TRUE.
-    ros_NewF(4) = .TRUE.
-    ros_NewF(5) = .TRUE.
-    ros_NewF(6) = .TRUE.
-
-!~~~> ros_ELO  = estimator of local order - the minimum between the
-!        main and the embedded scheme orders plus 1
-    ros_ELO = 4.0_dp
-
-  END SUBROUTINE Rodas4
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!   End of the set of internal Rosenbrock subroutines
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-END SUBROUTINE Rosenbrock
Index: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/rosenbrock_vec.f90.save
===================================================================
--- palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/rosenbrock_vec.f90.save	(revision 2718)
+++ 	(revision )
@@ -1,1152 +1,0 @@
-SUBROUTINE Rosenbrock(Y,Tstart,Tend, &
-           AbsTol,RelTol,            &
-           RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!
-!    Solves the system y'=F(t,y) using a Rosenbrock method defined by:
-!
-!     G = 1/(H*gamma(1)) - Jac(t0,Y0)
-!     T_i = t0 + Alpha(i)*H
-!     Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
-!     G * K_i = Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j +
-!         gamma(i)*dF/dT(t0, Y0)
-!     Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
-!
-!    For details on Rosenbrock methods and their implementation consult:
-!      E. Hairer and G. Wanner
-!      "Solving ODEs II. Stiff and differential-algebraic problems".
-!      Springer series in computational mathematics, Springer-Verlag, 1996.
-!    The codes contained in the book inspired this implementation.
-!
-!    (C)  Adrian Sandu, August 2004
-!    Virginia Polytechnic Institute and State University
-!    Contact: sandu@cs.vt.edu
-!    This implementation is part of KPP - the Kinetic PreProcessor
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!
-!~~~>   INPUT ARGUMENTS:
-!
-!-     Y(NVAR)    = vector of initial conditions (at T=Tstart)
-!-    [Tstart,Tend]  = time range of integration
-!     (if Tstart>Tend the integration is performed backwards in time)
-!-    RelTol, AbsTol = user precribed accuracy
-!- SUBROUTINE  Fun( T, Y, Ydot ) = ODE function,
-!                       returns Ydot = Y' = F(T,Y)
-!- SUBROUTINE  Jac( T, Y, Jcb ) = Jacobian of the ODE function,
-!                       returns Jcb = dFun/dY
-!-    ICNTRL(1:20)    = integer inputs parameters
-!-    RCNTRL(1:20)    = real inputs parameters
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!
-!~~~>     OUTPUT ARGUMENTS:
-!
-!-    Y(NVAR)    -> vector of final states (at T->Tend)
-!-    ISTATUS(1:20)   -> integer output parameters
-!-    RSTATUS(1:20)   -> real output parameters
-!-    IERR            -> job status upon return
-!                        success (positive value) or
-!                        failure (negative value)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!
-!~~~>     INPUT PARAMETERS:
-!
-!    Note: For input parameters equal to zero the default values of the
-!       corresponding variables are used.
-!
-!    ICNTRL(1) = 1: F = F(y)   Independent of T (AUTONOMOUS)
-!              = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
-!
-!    ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
-!              = 1: AbsTol, RelTol are scalars
-!
-!    ICNTRL(3)  -> selection of a particular Rosenbrock method
-!        = 0 :  default method is Rodas3
-!        = 1 :  method is  Ros2
-!        = 2 :  method is  Ros3
-!        = 3 :  method is  Ros4
-!        = 4 :  method is  Rodas3
-!        = 5:   method is  Rodas4
-!
-!    ICNTRL(4)  -> maximum number of integration steps
-!        For ICNTRL(4)=0) the default value of 100000 is used
-!
-!    RCNTRL(1)  -> Hmin, lower bound for the integration step size
-!          It is strongly recommended to keep Hmin = ZERO
-!    RCNTRL(2)  -> Hmax, upper bound for the integration step size
-!    RCNTRL(3)  -> Hstart, starting value for the integration step size
-!
-!    RCNTRL(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
-!    RCNTRL(5)  -> FacMax, upper bound on step increase factor (default=6)
-!    RCNTRL(6)  -> FacRej, step decrease factor after multiple rejections
-!            (default=0.1)
-!    RCNTRL(7)  -> FacSafe, by which the new step is slightly smaller
-!         than the predicted value  (default=0.9)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!
-!~~~>     OUTPUT PARAMETERS:
-!
-!    Note: each call to Rosenbrock adds the current no. of fcn calls
-!      to previous value of ISTATUS(1), and similar for the other params.
-!      Set ISTATUS(1:20) = 0 before call to avoid this accumulation.
-!
-!    ISTATUS(1) = No. of function calls
-!    ISTATUS(2) = No. of jacobian calls
-!    ISTATUS(3) = No. of steps
-!    ISTATUS(4) = No. of accepted steps
-!    ISTATUS(5) = No. of rejected steps (except at the beginning)
-!    ISTATUS(6) = No. of LU decompositions
-!    ISTATUS(7) = No. of forward/backward substitutions
-!    ISTATUS(8) = No. of singular matrix decompositions
-!
-!    RSTATUS(1)  -> Texit, the time corresponding to the
-!                   computed Y upon return
-!    RSTATUS(2)  -> Hexit, last accepted step before exit
-!    For multiple restarts, use Hexit as Hstart in the following run
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-  IMPLICIT NONE
-
-!~~~>  Arguments
-   REAL(kind=dp), INTENT(INOUT) :: Y(:,:)
-   REAL(kind=dp), INTENT(IN)   :: Tstart,Tend
-   REAL(kind=dp), INTENT(IN)   :: AbsTol(NVAR),RelTol(NVAR)
-   INTEGER, INTENT(IN)    :: ICNTRL(20)
-   REAL(kind=dp), INTENT(IN)   :: RCNTRL(20)
-   INTEGER, INTENT(INOUT) :: ISTATUS(20)
-   REAL(kind=dp), INTENT(INOUT) :: RSTATUS(20)
-   INTEGER, INTENT(OUT)   :: IERR
-!~~~>  The method parameters
-   INTEGER, PARAMETER :: Smax = 6
-   INTEGER  :: Method, ros_S
-   REAL(kind=dp), DIMENSION(Smax) :: ros_M, ros_E, ros_Alpha, ros_Gamma
-   REAL(kind=dp), DIMENSION(Smax*(Smax-1)/2) :: ros_A, ros_C
-   REAL(kind=dp) :: ros_ELO
-   LOGICAL, DIMENSION(Smax) :: ros_NewF
-   CHARACTER(LEN=12) :: ros_Name
-!~~~>  Local variables
-   REAL(kind=dp) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
-   REAL(kind=dp) :: Hmin, Hmax, Hstart, Hexit
-   REAL(kind=dp) :: Texit
-   INTEGER :: i, UplimTol, Max_no_steps
-   LOGICAL :: Autonomous, VectorTol
-!~~~>   Parameters
-   REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp, ONE  = 1.0_dp
-   REAL(kind=dp), PARAMETER :: DeltaMin = 1.0E-5_dp
-
-!~~~>  Initialize statistics
-   Nfun = ISTATUS(ifun)
-   Njac = ISTATUS(ijac)
-   Nstp = ISTATUS(istp)
-   Nacc = ISTATUS(iacc)
-   Nrej = ISTATUS(irej)
-   Ndec = ISTATUS(idec)
-   Nsol = ISTATUS(isol)
-   Nsng = ISTATUS(isng)
-
-!~~~>  Autonomous or time dependent ODE. Default is time dependent.
-   Autonomous = .NOT.(ICNTRL(1) == 0)
-
-!~~~>  For Scalar tolerances (ICNTRL(2).NE.0)  the code uses AbsTol(1) and RelTol(1)
-!   For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR)
-   IF (ICNTRL(2) == 0) THEN
-      VectorTol = .TRUE.
-         UplimTol  = NVAR
-   ELSE
-      VectorTol = .FALSE.
-         UplimTol  = 1
-   END IF
-
-!~~~>  The particular Rosenbrock method chosen
-   IF (ICNTRL(3) == 0) THEN
-      Method = 4
-   ELSEIF ( (ICNTRL(3) >= 1).AND.(ICNTRL(3) <= 5) ) THEN
-      Method = ICNTRL(3)
-   ELSE
-      PRINT * , 'User-selected Rosenbrock method: ICNTRL(3)=', Method
-      CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
-      RETURN
-   END IF
-
-!~~~>   The maximum number of steps admitted
-   IF (ICNTRL(4) == 0) THEN
-      Max_no_steps = 100000 ! 200
-   ELSEIF (ICNTRL(4) > 0) THEN
-      Max_no_steps=ICNTRL(4)
-   ELSE
-      PRINT * ,'User-selected max no. of steps: ICNTRL(4)=',ICNTRL(4)
-      CALL ros_ErrorMsg(-1,Tstart,ZERO,IERR)
-      RETURN
-   END IF
-
-!~~~>  Unit roundoff (1+Roundoff>1)
-   Roundoff = WLAMCH('E')
-
-!~~~>  Lower bound on the step size: (positive value)
-   IF (RCNTRL(1) == ZERO) THEN
-      Hmin = ZERO
-   ELSEIF (RCNTRL(1) > ZERO) THEN
-      Hmin = RCNTRL(1)
-   ELSE
-      PRINT * , 'User-selected Hmin: RCNTRL(1)=', RCNTRL(1)
-      CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
-      RETURN
-   END IF
-!~~~>  Upper bound on the step size: (positive value)
-   IF (RCNTRL(2) == ZERO) THEN
-      Hmax = ABS(Tend-Tstart)
-   ELSEIF (RCNTRL(2) > ZERO) THEN
-      Hmax = MIN(ABS(RCNTRL(2)),ABS(Tend-Tstart))
-   ELSE
-      PRINT * , 'User-selected Hmax: RCNTRL(2)=', RCNTRL(2)
-      CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
-      RETURN
-   END IF
-!~~~>  Starting step size: (positive value)
-   IF (RCNTRL(3) == ZERO) THEN
-      Hstart = MAX(Hmin,DeltaMin)
-   ELSEIF (RCNTRL(3) > ZERO) THEN
-      Hstart = MIN(ABS(RCNTRL(3)),ABS(Tend-Tstart))
-   ELSE
-      PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
-      CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
-      RETURN
-   END IF
-!~~~>  Step size can be changed s.t.  FacMin < Hnew/Hexit < FacMax
-   IF (RCNTRL(4) == ZERO) THEN
-      FacMin = 0.2_dp
-   ELSEIF (RCNTRL(4) > ZERO) THEN
-      FacMin = RCNTRL(4)
-   ELSE
-      PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
-      CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
-      RETURN
-   END IF
-   IF (RCNTRL(5) == ZERO) THEN
-      FacMax = 6.0_dp
-   ELSEIF (RCNTRL(5) > ZERO) THEN
-      FacMax = RCNTRL(5)
-   ELSE
-      PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
-      CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
-      RETURN
-   END IF
-!~~~>   FacRej: Factor to decrease step after 2 succesive rejections
-   IF (RCNTRL(6) == ZERO) THEN
-      FacRej = 0.1_dp
-   ELSEIF (RCNTRL(6) > ZERO) THEN
-      FacRej = RCNTRL(6)
-   ELSE
-      PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
-      CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
-      RETURN
-   END IF
-!~~~>   FacSafe: Safety Factor in the computation of new step size
-   IF (RCNTRL(7) == ZERO) THEN
-      FacSafe = 0.9_dp
-   ELSEIF (RCNTRL(7) > ZERO) THEN
-      FacSafe = RCNTRL(7)
-   ELSE
-      PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
-      CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
-      RETURN
-   END IF
-!~~~>  Check if tolerances are reasonable
-    DO i=1,UplimTol
-      IF ( (AbsTol(i) <= ZERO) .OR. (RelTol(i) <= 10.0_dp*Roundoff) &
-         .OR. (RelTol(i) >= 1.0_dp) ) THEN
-        PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
-        PRINT * , ' RelTol(',i,') = ',RelTol(i)
-        CALL ros_ErrorMsg(-5,Tstart,ZERO,IERR)
-        RETURN
-      END IF
-    END DO
-
-
-!~~~>   Initialize the particular Rosenbrock method
-   SELECT CASE (Method)
-     CASE (1)
-       CALL Ros2(ros_S, ros_A, ros_C, ros_M, ros_E,   &
-          ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
-     CASE (2)
-       CALL Ros3(ros_S, ros_A, ros_C, ros_M, ros_E,   &
-          ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
-     CASE (3)
-       CALL Ros4(ros_S, ros_A, ros_C, ros_M, ros_E,   &
-          ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
-     CASE (4)
-       CALL Rodas3(ros_S, ros_A, ros_C, ros_M, ros_E, &
-          ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
-     CASE (5)
-       CALL Rodas4(ros_S, ros_A, ros_C, ros_M, ros_E, &
-          ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
-     CASE DEFAULT
-       PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=', Method
-       CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
-       RETURN
-   END SELECT
-
-!~~~>  CALL Rosenbrock method
-   CALL ros_Integrator(Y,Tstart,Tend,Texit,      &
-        AbsTol, RelTol,                          &
-!  Rosenbrock method coefficients
-        ros_S, ros_M, ros_E, ros_A, ros_C,       &
-        ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
-!  Integration parameters
-        Autonomous, VectorTol, Max_no_steps,     &
-        Roundoff, Hmin, Hmax, Hstart, Hexit,     &
-        FacMin, FacMax, FacRej, FacSafe,         &
-!  Error indicator
-        IERR)
-
-
-!~~~>  Collect run statistics
-   ISTATUS(ifun) = Nfun
-   ISTATUS(ijac) = Njac
-   ISTATUS(istp) = Nstp
-   ISTATUS(iacc) = Nacc
-   ISTATUS(irej) = Nrej
-   ISTATUS(idec) = Ndec
-   ISTATUS(isol) = Nsol
-   ISTATUS(isng) = Nsng
-!~~~> Last T and H
-   RSTATUS(itexit) = Texit
-   RSTATUS(ihexit) = Hexit
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-CONTAINS !  SUBROUTINES internal to Rosenbrock
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- SUBROUTINE ros_ErrorMsg(Code,T,H,IERR)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!    Handles all error messages
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-   REAL(kind=dp), INTENT(IN) :: T, H
-   INTEGER, INTENT(IN)  :: Code
-   INTEGER, INTENT(OUT) :: IERR
-
-   IERR = Code
-   PRINT * , &
-     'Forced exit from Rosenbrock due to the following error:'
-   IF ((Code>=-8).AND.(Code<=-1)) THEN
-     PRINT *, IERR_NAMES(Code)
-   ELSE
-     PRINT *, 'Unknown Error code: ', Code
-   ENDIF
-
-   PRINT *, "T=", T, "and H=", H
-
- END SUBROUTINE ros_ErrorMsg
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- SUBROUTINE ros_Integrator (Y, Tstart, Tend, Tout,     &
-        AbsTol, RelTol,                          &
-!~~~> Rosenbrock method coefficients
-        ros_S, ros_M, ros_E, ros_A, ros_C,       &
-        ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
-!~~~> Integration parameters
-        Autonomous, VectorTol, Max_no_steps,     &
-        Roundoff, Hmin, Hmax, Hstart, Hexit,     &
-        FacMin, FacMax, FacRej, FacSafe,         &
-!~~~> Error indicator
-        IERR )
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!   Template for the implementation of a generic Rosenbrock method
-!      defined by ros_S (no of stages)
-!      and its coefficients ros_{A,C,M,E,Alpha,Gamma}
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-  USE      messy_kpp_compress,         ONLY: kco_initialize, kco_compress, &
-                                             kco_finalize, cell_done
-  IMPLICIT NONE
-
-!~~~> Input: the initial condition at Tstart; Output: the solution at T
-   REAL(kind=dp), INTENT(INOUT) :: Y(:,:)
-!~~~> Input: integration interval
-   REAL(kind=dp), INTENT(IN) :: Tstart,Tend
-!~~~> Output: time at which the solution is returned (T=Tend if success)
-   REAL(kind=dp), INTENT(OUT) ::  Tout
-!~~~> Input: tolerances
-   REAL(kind=dp), INTENT(IN) ::  AbsTol(NVAR), RelTol(NVAR)
-!~~~> Input: The Rosenbrock method parameters
-   INTEGER, INTENT(IN) ::  ros_S
-   REAL(kind=dp), INTENT(IN) :: ros_M(ros_S), ros_E(ros_S),  &
-       ros_Alpha(ros_S), ros_A(ros_S*(ros_S-1)/2), &
-       ros_Gamma(ros_S), ros_C(ros_S*(ros_S-1)/2), ros_ELO
-   LOGICAL, INTENT(IN) :: ros_NewF(ros_S)
-!~~~> Input: integration parameters
-   LOGICAL, INTENT(IN) :: Autonomous, VectorTol
-   REAL(kind=dp), INTENT(IN) :: Hstart, Hmin, Hmax
-   INTEGER, INTENT(IN) :: Max_no_steps
-   REAL(kind=dp), INTENT(IN) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
-!~~~> Output: last accepted step
-   REAL(kind=dp), INTENT(OUT) :: Hexit
-!~~~> Output: Error indicator
-   INTEGER, INTENT(OUT) :: IERR
-! ~~~~ Local variables
-   REAL(kind=dp) :: Ynew(VL_DIM,NVAR), Fcn0(VL_DIM,NVAR), Fcn(VL_DIM,NVAR)
-   REAL(kind=dp) :: K(VL_DIM,NVAR*ros_S)
-   REAL(kind=dp) :: Jac0(VL_DIM,LU_NONZERO), Ghimj(VL_DIM,LU_NONZERO)
-   REAL(kind=dp) :: Yerr(VL_DIM,NVAR)
-
-!~~~>  Local parameters
-   INTEGER :: ioffset, i, j, istage
-   REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp, ONE  = 1.0_dp
-   REAL(kind=dp), PARAMETER :: DeltaMin = 1.0E-5_dp
-!~~~>  Locally called functions
-!    REAL(kind=dp) WLAMCH
-!    EXTERNAL WLAMCH
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-   real(kind=dp)                  :: deltat
-!kk   logical,parameter              :: fixed_step=.true.
-   
-   integer,save             :: icount=0
-
-!  vector compression algorithm
-
-   REAL(kind=dp), dimension(VL_dim)      :: T
-   REAL(kind=dp), dimension(VL_dim)      :: H, Hnew, HC, HC2, HG, Fac, Tau, Err
-   LOGICAL, dimension(VL_dim)            :: RejectLastH, RejectMoreH
-   LOGICAL, dimension(VL_dim)            :: Accept_step, error_flag
-   integer                               :: vl_save
-
-   REAL(kind=dp), dimension(size(FIX,1),size(FIX,2))      :: FIX_lo
-
-!~~~>  Initial preparations
-   T(1:VL)  = Tstart
-   Hexit    = 0.0_dp
-   H(1:VL)  = MIN(Hstart,Hmax)
-   Tout     = Tend
-   icount   = icount+1
-   where (ABS(H(1:VL)) <= 10.0_dp*Roundoff) H(1:VL) = DeltaMin
-
-   deltat=Tend-Tstart
-
-   RejectLastH(1:VL)=.FALSE.
-   RejectMoreH(1:VL)=.FALSE.
-
-!~~~> Time loop begins below
-
-   Nstp       = 0
-   error_flag(1:VL) = .false.
-   vl_save          = VL
-   FIX_lo(1:VL,:)   = FIX(1:VL,:)
-   Kacc(1:VL)       = 0
-   Krej(1:VL)       = 0
-
-!  intialize compression algorithm
-
-   IF (.not. l_fixed_step) then
-     if( return_A ) then
-       call kco_initialize (VL, y, Kacc, Krej, a=a)
-     else
-       call kco_initialize (VL, y, Kacc, Krej)
-     end if 
-   end if
-
-TimeLoop: DO 
-
-   IF ( Nstp > Max_no_steps ) THEN  ! Too many steps
-      CALL ros_ErrorMsg(-6,T(1),H(1),IERR)
-! mz_pj_20070514: replace 'call abort()' by
-!     error_flag(1:VL) = .TRUE.
-! NOTE: Nstp is the maximum number of steps of the remaining current vector
-      call abort()                                   !Vorlaeufig
-      RETURN
-   END IF
-   
-   Where ( (T(1:VL)+0.1_dp*H(1:VL)) == T(1:VL))      ! Step size too small
-      error_flag(1:VL) = .TRUE.  
-   endwhere
-
-   if( any(error_flag))  then
-      CALL ros_ErrorMsg(-7,T(1),H(1),IERR)
-      call abort()                                   !Vorlaeufig
-      RETURN
-   end if
-
-!~~~>  Limit H if necessary to avoid going beyond Tend
-!kk   Hexit = H
-   H(1:VL) = MIN(H(1:VL),ABS(Tend-T(1:VL)))
-
-!~~~>   Compute the function at current time
-   CALL Fun( Y, FIX_lo, RCONST, Fcn0)
-
-
-!~~~>   Compute the Jacobian at current time
-   CALL Jac_SP( Y, FIX_lo, RCONST, Jac0)
-
-
-!~~~>  Repeat step calculation until current step accepted
-
-   CALL ros_PrepareMatrix(H,ros_Gamma(1), Jac0,Ghimj)
-
-!~~~>   Compute the stages
-Stage: DO istage = 1, ros_S
-
-      ! Current istage offset. Current istage vector is K(ioffset+1:ioffset+NVAR)
-       ioffset = NVAR*(istage-1)
-
-      ! For the 1st istage the function has been computed previously
-       IF ( istage == 1 ) THEN
-         Fcn(1:VL,1:NVAR) = Fcn0(1:VL,1:NVAR) 
-      ! istage>1 and a new function evaluation is needed at the current istage
-       ELSEIF ( ros_NewF(istage) ) THEN
-         select case (istage-1)
-          case (1)
-           Ynew(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
-           Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) +  K(1:VL,NVAR*(1- 1)+ 1:NVAR*1) &
-                        *ros_A((istage-1)*(istage-2)/2+ 1)
-          case (2)
-           Ynew(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
-           Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) +  K(1:VL,NVAR*(2- 1)+ 1:NVAR*2) &
-                        *ros_A((istage-1)*(istage-2)/2+ 2)
-           Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) +  K(1:VL,NVAR*(2- 1)+ 1:NVAR*2) &
-                        *ros_A((istage-1)*(istage-2)/2+ 2)
-          case DEFAULT
-           Ynew(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
-           DO j = 1, istage-1
-             Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) +  K(1:VL,NVAR*(j- 1)+ 1:NVAR*j) &
-                          *ros_A((istage-1)*(istage-2)/2+ j)
-           END DO
-         end select
-
-         CALL Fun( Ynew, FIX_lo, RCONST, Fcn)
-       END IF ! if istage == 1 elseif ros_NewF(istage)
-
-       select case (istage-1)
-        case (1)
-         K(1:VL,ioffset+ 1:ioffset+ NVAR) = Fcn(1:VL,1:NVAR)
-         j = 1
-         HC(1:VL) = ros_C((istage-1)*(istage-2)/2+ j)/(H(1:VL))
-         do i=1,NVAR
-           K(1:VL,ioffset+ i) = K(1:VL,ioffset+ i) +  K(1:VL,NVAR*(j- 1)+ i) *HC(1:VL)
-         end do
-        case (2)
-         K(1:VL,ioffset+ 1:ioffset+ NVAR) = Fcn(1:VL,1:NVAR)
-         HC(1:VL)  = ros_C((istage-1)*(istage-2)/2+ 1)/(H(1:VL))
-         HC2(1:VL) = ros_C((istage-1)*(istage-2)/2+ 2)/(H(1:VL))
-         do i=1,NVAR
-           K(1:VL,ioffset+ i) = K(1:VL,ioffset+ i) +  K(1:VL,NVAR*(1- 1)+ i) *HC(1:VL)
-           K(1:VL,ioffset+ i) = K(1:VL,ioffset+ i) +  K(1:VL,NVAR*(2- 1)+ i) *HC2(1:VL)
-         end do
-        case DEFAULT
-         K(1:VL,ioffset+ 1:ioffset+ NVAR) = Fcn(1:VL,1:NVAR)
-         DO j = 1, istage-1
-           HC(1:VL) = ros_C((istage-1)*(istage-2)/2+ j)/(H(1:VL))
-           do i=1,NVAR
-             K(1:VL,ioffset+ i) = K(1:VL,ioffset+ i) +  K(1:VL,NVAR*(j- 1)+ i) *HC(1:VL)
-           end do
-         END DO
-       end select
-
-       CALL ros_Solve(Ghimj, K(:,ioffset+1:ioffset+NVAR))
-
-   END DO Stage
-
-
-   if(ros_S == 3)   then
-     YNEW(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
-     Yerr(1:VL,1:NVAR) = ZERO
-
-!~~~>  Compute the new solution
-     Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR*(1-1)+1:NVAR*1) * ros_m(1)
-     Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR*(2-1)+1:NVAR*2) * ros_m(2)
-     Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR*(3-1)+1:NVAR*3) * ros_m(3)
-
-!~~~>  Compute the error estimation
-
-     Yerr(1:VL,1:NVAR) = Yerr(1:VL,1:NVAR) + K(1:VL,NVAR*(1-1)+1:NVAR*1) * ros_E(1)
-     Yerr(1:VL,1:NVAR) = Yerr(1:VL,1:NVAR) + K(1:VL,NVAR*(2-1)+1:NVAR*2) * ros_E(2)
-     Yerr(1:VL,1:NVAR) = Yerr(1:VL,1:NVAR) + K(1:VL,NVAR*(3-1)+1:NVAR*3) * ros_E(3)
-   else
-     YNEW(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
-     Yerr(1:VL,1:NVAR) = ZERO
-     DO j=1,ros_S
-!~~~>  Compute the new solution
-         Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + &
-                        K(1:VL,NVAR*(j-1)+1:NVAR*j) * ros_m(j)
-
-!~~~>  Compute the error estimation
-
-       Yerr(1:VL,1:NVAR) = Yerr(1:VL,1:NVAR) + &
-                    K(1:VL,NVAR*(j-1)+1:NVAR*j) * ros_E(j)
-     END DO
-   end if
-
-
-!~~~>  Check the error magnitude and adjust step size
-   Nstp = Nstp+1
-
-   if(l_fixed_step)   then
-     Err  = 0.0
-     H    = t_steps(Nstp)*deltat
-     Hnew = t_steps(Nstp)*deltat
-   else
-     Err(1:VL) = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol )
-!~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
-     where(Err(1:VL) >= 10.0_dp*Roundoff)
-       Fac(1:VL)  = MIN(FacMax,MAX(FacMin,FacSafe/Err(1:VL)**(ONE/ros_ELO)))
-     elsewhere
-       Fac(1:VL)  = 1000.                        !Force large time step
-     endwhere
-     Hnew(1:VL) = H(1:VL)*Fac(1:VL)
-   end if
-
-   Accept_step(1:VL) = (Err(1:VL) <= ONE).OR.(H(1:VL) <= Hmin)
-   where (Accept_step(1:VL))
-      Kacc(1:VL) = Kacc(1:VL)+1
-      T(1:VL) = T(1:VL) + H(1:VL)
-      Hnew(1:VL) = MAX(Hmin,MIN(Hnew(1:VL),Hmax))
-   endwhere
-   do i=1,NVAR
-     where (Accept_step(1:VL))
-       Y(1:VL,i) = MAX(Ynew(1:VL,i),ZERO)
-     endwhere
-   end do
-   where (Accept_step(1:VL) .and. RejectLastH(1:VL))
-      Hnew(1:VL) = MIN(Hnew(1:VL),H(1:VL))
-   endwhere
-   where (Accept_step(1:VL))
-      RejectLastH(1:VL) = .FALSE.
-      RejectMoreH(1:VL) = .FALSE.
-      H(1:VL) = Hnew(1:VL)
-   endwhere
-!kk   ELSE           !~~~> Reject step
-   where (.NOT. Accept_step(1:VL) .and. RejectMoreH(1:VL))
-      Hnew(1:VL) = H(1:VL)*FacRej
-   endwhere
-   where (.NOT. Accept_step(1:VL))
-      RejectMoreH(1:VL) = RejectLastH(1:VL)
-      RejectLastH(1:VL) = .TRUE.
-      H(1:VL) = Hnew(1:VL)
-   endwhere
-   where ( .NOT. Accept_step(1:VL) .and. Kacc(1:VL) >= 1)
-      Krej(1:VL) = Krej(1:VL)+1
-   endwhere
-
-
-   if(l_fixed_step)   then
-      if(Nstp >= Nfsteps) EXIT
-   else
-     cell_done(1:VL) = ( .NOT. (T(1:VL)-Tend)+Roundoff <= ZERO  )  
-     call kco_compress (VL, T, H, Hnew, error_flag, y, RCONST, FIX_lo,         &
-                                     RejectLastH, RejectMoreH, Kacc, Krej, a=a)
-     if(VL == 0)   EXIT
-   end if
-
-!kk   write(0,*) '  aaa ',icount,Nstp,count(cell_done),count(Accept_step),VL,minval(T)
-   END DO TimeLoop
-
-!~~~> Succesful exit
-   IERR = 1  !~~~> The integration was successful
-
-   IF (.not. l_fixed_step) then
-     call kco_finalize (y, Kacc, Krej, a=a)
-   end if
-
-   VL = VL_save
-
-!kk   write(0,*) 'bbb ',icount,Nstp,maxval(Kacc),maxval(Krej)
-!kk   if(icount >= 10)  call abort()
-
-  END SUBROUTINE ros_Integrator
-
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  FUNCTION ros_ErrorNorm ( Y, Ynew, Yerr, &
-               AbsTol, RelTol, VectorTol )
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!~~~> Computes the "scaled norm" of the error vector Yerr
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-   IMPLICIT NONE
-
-! Input arguments
-   REAL(kind=dp), INTENT(IN) :: Y(:,:), Ynew(:,:), &
-          Yerr(:,:), AbsTol(NVAR), RelTol(NVAR)
-   LOGICAL, INTENT(IN) ::  VectorTol
-! Local variables
-   REAL(kind=dp), dimension(VL)       :: Err, Scale, Ymax
-   INTEGER  :: i
-   REAL(kind=dp), PARAMETER           :: ZERO = 0.0_dp
-   REAL(kind=dp), dimension(VL)       :: ros_ErrorNorm
-
-   Err = ZERO
-   DO i=1,NVAR
-     Ymax = MAX( ABS(Y(1:VL,i)),ABS(Ynew(1:VL,i)))
-     IF (VectorTol) THEN
-       Scale = AbsTol(i)+ RelTol(i)*Ymax
-     ELSE
-       Scale = AbsTol(1)+ RelTol(1)*Ymax
-     END IF
-     Err = Err+ (Yerr(1:VL,i)/Scale)**2
-   END DO
-   Err  = SQRT(Err/NVAR)
-
-   ros_ErrorNorm = Err
-
-  END FUNCTION ros_ErrorNorm
-
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-  SUBROUTINE ros_PrepareMatrix ( H, gam, Jac0, Ghimj )
-
-! --- --- --- --- --- --- --- --- --- --- --- --- ---
-!  Prepares the LHS matrix for stage calculations
-!  1.  Construct Ghimj = 1/(H*ham) - Jac0
-!      "(Gamma H) Inverse Minus Jacobian"
-!  2.  Repeat LU decomposition of Ghimj until successful.
-!       -half the step size if LU decomposition fails and retry
-!       -exit after 5 consecutive fails
-! --- --- --- --- --- --- --- --- --- --- --- --- ---
-   IMPLICIT NONE
-
-!~~~> Input arguments
-   REAL(kind=dp), INTENT(IN) ::  Jac0(:,:)
-   REAL(kind=dp), INTENT(IN) ::  gam
-!~~~> Output arguments
-   REAL(kind=dp), INTENT(OUT) :: Ghimj(:,:)
-!~~~> Inout arguments
-   REAL(kind=dp), INTENT(INOUT) :: H(:)   ! step size is decreased when LU fails
-!~~~> Local variables
-   INTEGER  :: i, ising, Nconsecutive
-   REAL(kind=dp), dimension(size(H)) :: ghinv
-   REAL(kind=dp), PARAMETER :: ONE  = 1.0_dp, HALF = 0.5_dp
-
-   Nconsecutive = 0
-
-
-!~~~>    Construct Ghimj = 1/(H*gam) - Jac0
-   Ghimj(1:VL,1:LU_NONZERO) = -ONE * JAC0(1:VL,1:LU_NONZERO)
-   ghinv(1:VL) = ONE/(H(1:VL)*gam)
-   DO i=1,NVAR
-     Ghimj(1:VL,LU_DIAG(i)) = Ghimj(1:VL,LU_DIAG(i))+ghinv(1:VL)
-   END DO
-
-!~~~>    Compute LU decomposition
-!CDIR noiexpand
-   CALL KppDecomp( Ghimj, ising )
-
-   Ndec = Ndec + VL
-
-  END SUBROUTINE ros_PrepareMatrix
-
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  SUBROUTINE ros_Solve( A, b )
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!  Template for the forward/backward substitution (using pre-computed LU decomposition)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-   IMPLICIT NONE
-!~~~> Input variables
-   REAL(kind=dp), INTENT(IN)         :: A(:,:)
-!~~~> InOut variables
-   REAL(kind=dp), INTENT(INOUT)      :: b(:,:)
-
-   CALL KppSolve( A, b )
-
-   Nsol = Nsol+VL
-
-  END SUBROUTINE ros_Solve
-
-
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  SUBROUTINE Ros2 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
-            ros_Gamma,ros_NewF,ros_ELO,ros_Name)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! --- AN L-STABLE METHOD, 2 stages, order 2
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-  IMPLICIT NONE
-
-   INTEGER, PARAMETER :: S=2
-   INTEGER, INTENT(OUT) ::  ros_S
-   REAL(kind=dp), DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
-   REAL(kind=dp), DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
-   REAL(kind=dp), INTENT(OUT) :: ros_ELO
-   LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
-   CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
-   DOUBLE PRECISION g
-
-    g = 1.0_dp + 1.0_dp/SQRT(2.0_dp)
-
-!~~~> Name of the method
-    ros_Name = 'ROS-2'
-!~~~> Number of stages
-    ros_S = S
-
-!~~~> The coefficient matrices A and C are strictly lower triangular.
-!   The lower triangular (subdiagonal) elements are stored in row-wise order:
-!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
-!   The general mapping formula is:
-!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
-!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
-
-    ros_A(1) = (1.0_dp)/g
-    ros_C(1) = (-2.0_dp)/g
-!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
-!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
-    ros_NewF(1) = .TRUE.
-    ros_NewF(2) = .TRUE.
-!~~~> M_i = Coefficients for new step solution
-    ros_M(1)= (3.0_dp)/(2.0_dp*g)
-    ros_M(2)= (1.0_dp)/(2.0_dp*g)
-! E_i = Coefficients for error estimator
-    ros_E(1) = 1.0_dp/(2.0_dp*g)
-    ros_E(2) = 1.0_dp/(2.0_dp*g)
-!~~~> ros_ELO = estimator of local order - the minimum between the
-!    main and the embedded scheme orders plus one
-    ros_ELO = 2.0_dp
-!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
-    ros_Alpha(1) = 0.0_dp
-    ros_Alpha(2) = 1.0_dp
-!~~~> Gamma_i = \sum_j  gamma_{i,j}
-    ros_Gamma(1) = g
-    ros_Gamma(2) =-g
-
- END SUBROUTINE Ros2
-
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  SUBROUTINE Ros3 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
-           ros_Gamma,ros_NewF,ros_ELO,ros_Name)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-  IMPLICIT NONE
-
-   INTEGER, PARAMETER :: S=3
-   INTEGER, INTENT(OUT) ::  ros_S
-   REAL(kind=dp), DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
-   REAL(kind=dp), DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
-   REAL(kind=dp), INTENT(OUT) :: ros_ELO
-   LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
-   CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
-
-!~~~> Name of the method
-   ros_Name = 'ROS-3'
-!~~~> Number of stages
-   ros_S = S
-
-!~~~> The coefficient matrices A and C are strictly lower triangular.
-!   The lower triangular (subdiagonal) elements are stored in row-wise order:
-!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
-!   The general mapping formula is:
-!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
-!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
-
-   ros_A(1)= 1.0_dp
-   ros_A(2)= 1.0_dp
-   ros_A(3)= 0.0_dp
-
-   ros_C(1) = -0.10156171083877702091975600115545E+01_dp
-   ros_C(2) =  0.40759956452537699824805835358067E+01_dp
-   ros_C(3) =  0.92076794298330791242156818474003E+01_dp
-!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
-!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
-   ros_NewF(1) = .TRUE.
-   ros_NewF(2) = .TRUE.
-   ros_NewF(3) = .FALSE.
-!~~~> M_i = Coefficients for new step solution
-   ros_M(1) =  0.1E+01_dp
-   ros_M(2) =  0.61697947043828245592553615689730E+01_dp
-   ros_M(3) = -0.42772256543218573326238373806514E+00_dp
-! E_i = Coefficients for error estimator
-   ros_E(1) =  0.5E+00_dp
-   ros_E(2) = -0.29079558716805469821718236208017E+01_dp
-   ros_E(3) =  0.22354069897811569627360909276199E+00_dp
-!~~~> ros_ELO = estimator of local order - the minimum between the
-!    main and the embedded scheme orders plus 1
-   ros_ELO = 3.0_dp
-!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
-   ros_Alpha(1)= 0.0E+00_dp
-   ros_Alpha(2)= 0.43586652150845899941601945119356E+00_dp
-   ros_Alpha(3)= 0.43586652150845899941601945119356E+00_dp
-!~~~> Gamma_i = \sum_j  gamma_{i,j}
-   ros_Gamma(1)= 0.43586652150845899941601945119356E+00_dp
-   ros_Gamma(2)= 0.24291996454816804366592249683314E+00_dp
-   ros_Gamma(3)= 0.21851380027664058511513169485832E+01_dp
-
-  END SUBROUTINE Ros3
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  SUBROUTINE Ros4 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
-           ros_Gamma,ros_NewF,ros_ELO,ros_Name)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!     L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
-!     L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
-!
-!      E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
-!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
-!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
-!      SPRINGER-VERLAG (1990)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-  IMPLICIT NONE
-
-   INTEGER, PARAMETER :: S=4
-   INTEGER, INTENT(OUT) ::  ros_S
-   REAL(kind=dp), DIMENSION(4), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
-   REAL(kind=dp), DIMENSION(6), INTENT(OUT) :: ros_A, ros_C
-   REAL(kind=dp), INTENT(OUT) :: ros_ELO
-   LOGICAL, DIMENSION(4), INTENT(OUT) :: ros_NewF
-   CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
-   DOUBLE PRECISION g
-
-!~~~> Name of the method
-   ros_Name = 'ROS-4'
-!~~~> Number of stages
-   ros_S = S
-
-!~~~> The coefficient matrices A and C are strictly lower triangular.
-!   The lower triangular (subdiagonal) elements are stored in row-wise order:
-!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
-!   The general mapping formula is:
-!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
-!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
-
-   ros_A(1) = 0.2000000000000000E+01_dp
-   ros_A(2) = 0.1867943637803922E+01_dp
-   ros_A(3) = 0.2344449711399156E+00_dp
-   ros_A(4) = ros_A(2)
-   ros_A(5) = ros_A(3)
-   ros_A(6) = 0.0_dp
-
-   ros_C(1) =-0.7137615036412310E+01_dp
-   ros_C(2) = 0.2580708087951457E+01_dp
-   ros_C(3) = 0.6515950076447975E+00_dp
-   ros_C(4) =-0.2137148994382534E+01_dp
-   ros_C(5) =-0.3214669691237626E+00_dp
-   ros_C(6) =-0.6949742501781779E+00_dp
-!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
-!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
-   ros_NewF(1)  = .TRUE.
-   ros_NewF(2)  = .TRUE.
-   ros_NewF(3)  = .TRUE.
-   ros_NewF(4)  = .FALSE.
-!~~~> M_i = Coefficients for new step solution
-   ros_M(1) = 0.2255570073418735E+01_dp
-   ros_M(2) = 0.2870493262186792E+00_dp
-   ros_M(3) = 0.4353179431840180E+00_dp
-   ros_M(4) = 0.1093502252409163E+01_dp
-!~~~> E_i  = Coefficients for error estimator
-   ros_E(1) =-0.2815431932141155E+00_dp
-   ros_E(2) =-0.7276199124938920E-01_dp
-   ros_E(3) =-0.1082196201495311E+00_dp
-   ros_E(4) =-0.1093502252409163E+01_dp
-!~~~> ros_ELO  = estimator of local order - the minimum between the
-!    main and the embedded scheme orders plus 1
-   ros_ELO  = 4.0_dp
-!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
-   ros_Alpha(1) = 0.0_dp
-   ros_Alpha(2) = 0.1145640000000000E+01_dp
-   ros_Alpha(3) = 0.6552168638155900E+00_dp
-   ros_Alpha(4) = ros_Alpha(3)
-!~~~> Gamma_i = \sum_j  gamma_{i,j}
-   ros_Gamma(1) = 0.5728200000000000E+00_dp
-   ros_Gamma(2) =-0.1769193891319233E+01_dp
-   ros_Gamma(3) = 0.7592633437920482E+00_dp
-   ros_Gamma(4) =-0.1049021087100450E+00_dp
-
-  END SUBROUTINE Ros4
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  SUBROUTINE Rodas3 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
-            ros_Gamma,ros_NewF,ros_ELO,ros_Name)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! --- A STIFFLY-STABLE METHOD, 4 stages, order 3
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-  IMPLICIT NONE
-
-   INTEGER, PARAMETER :: S=4
-   INTEGER, INTENT(OUT) ::  ros_S
-   REAL(kind=dp), DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
-   REAL(kind=dp), DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
-   REAL(kind=dp), INTENT(OUT) :: ros_ELO
-   LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
-   CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
-   DOUBLE PRECISION g
-
-!~~~> Name of the method
-   ros_Name = 'RODAS-3'
-!~~~> Number of stages
-   ros_S = S
-
-!~~~> The coefficient matrices A and C are strictly lower triangular.
-!   The lower triangular (subdiagonal) elements are stored in row-wise order:
-!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
-!   The general mapping formula is:
-!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
-!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
-
-   ros_A(1) = 0.0E+00_dp
-   ros_A(2) = 2.0E+00_dp
-   ros_A(3) = 0.0E+00_dp
-   ros_A(4) = 2.0E+00_dp
-   ros_A(5) = 0.0E+00_dp
-   ros_A(6) = 1.0E+00_dp
-
-   ros_C(1) = 4.0E+00_dp
-   ros_C(2) = 1.0E+00_dp
-   ros_C(3) =-1.0E+00_dp
-   ros_C(4) = 1.0E+00_dp
-   ros_C(5) =-1.0E+00_dp
-   ros_C(6) =-(8.0E+00_dp/3.0E+00_dp)
-
-!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
-!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
-   ros_NewF(1)  = .TRUE.
-   ros_NewF(2)  = .FALSE.
-   ros_NewF(3)  = .TRUE.
-   ros_NewF(4)  = .TRUE.
-!~~~> M_i = Coefficients for new step solution
-   ros_M(1) = 2.0E+00_dp
-   ros_M(2) = 0.0E+00_dp
-   ros_M(3) = 1.0E+00_dp
-   ros_M(4) = 1.0E+00_dp
-!~~~> E_i  = Coefficients for error estimator
-   ros_E(1) = 0.0E+00_dp
-   ros_E(2) = 0.0E+00_dp
-   ros_E(3) = 0.0E+00_dp
-   ros_E(4) = 1.0E+00_dp
-!~~~> ros_ELO  = estimator of local order - the minimum between the
-!    main and the embedded scheme orders plus 1
-   ros_ELO  = 3.0E+00_dp
-!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
-   ros_Alpha(1) = 0.0E+00_dp
-   ros_Alpha(2) = 0.0E+00_dp
-   ros_Alpha(3) = 1.0E+00_dp
-   ros_Alpha(4) = 1.0E+00_dp
-!~~~> Gamma_i = \sum_j  gamma_{i,j}
-   ros_Gamma(1) = 0.5E+00_dp
-   ros_Gamma(2) = 1.5E+00_dp
-   ros_Gamma(3) = 0.0E+00_dp
-   ros_Gamma(4) = 0.0E+00_dp
-
-  END SUBROUTINE Rodas3
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  SUBROUTINE Rodas4 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
-             ros_Gamma,ros_NewF,ros_ELO,ros_Name)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!     STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
-!
-!      E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
-!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
-!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
-!      SPRINGER-VERLAG (1996)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-  IMPLICIT NONE
-
-   INTEGER, PARAMETER :: S=6
-   INTEGER, INTENT(OUT) ::  ros_S
-   REAL(kind=dp), DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
-   REAL(kind=dp), DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
-   REAL(kind=dp), INTENT(OUT) :: ros_ELO
-   LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
-   CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
-   DOUBLE PRECISION g
-
-!~~~> Name of the method
-    ros_Name = 'RODAS-4'
-!~~~> Number of stages
-    ros_S = 6
-
-!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
-    ros_Alpha(1) = 0.000_dp
-    ros_Alpha(2) = 0.386_dp
-    ros_Alpha(3) = 0.210_dp
-    ros_Alpha(4) = 0.630_dp
-    ros_Alpha(5) = 1.000_dp
-    ros_Alpha(6) = 1.000_dp
-
-!~~~> Gamma_i = \sum_j  gamma_{i,j}
-    ros_Gamma(1) = 0.2500000000000000E+00_dp
-    ros_Gamma(2) =-0.1043000000000000E+00_dp
-    ros_Gamma(3) = 0.1035000000000000E+00_dp
-    ros_Gamma(4) =-0.3620000000000023E-01_dp
-    ros_Gamma(5) = 0.0_dp
-    ros_Gamma(6) = 0.0_dp
-
-!~~~> The coefficient matrices A and C are strictly lower triangular.
-!   The lower triangular (subdiagonal) elements are stored in row-wise order:
-!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
-!   The general mapping formula is:  A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
-!                  C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
-
-    ros_A(1) = 0.1544000000000000E+01_dp
-    ros_A(2) = 0.9466785280815826E+00_dp
-    ros_A(3) = 0.2557011698983284E+00_dp
-    ros_A(4) = 0.3314825187068521E+01_dp
-    ros_A(5) = 0.2896124015972201E+01_dp
-    ros_A(6) = 0.9986419139977817E+00_dp
-    ros_A(7) = 0.1221224509226641E+01_dp
-    ros_A(8) = 0.6019134481288629E+01_dp
-    ros_A(9) = 0.1253708332932087E+02_dp
-    ros_A(10) =-0.6878860361058950E+00_dp
-    ros_A(11) = ros_A(7)
-    ros_A(12) = ros_A(8)
-    ros_A(13) = ros_A(9)
-    ros_A(14) = ros_A(10)
-    ros_A(15) = 1.0E+00_dp
-
-    ros_C(1) =-0.5668800000000000E+01_dp
-    ros_C(2) =-0.2430093356833875E+01_dp
-    ros_C(3) =-0.2063599157091915E+00_dp
-    ros_C(4) =-0.1073529058151375E+00_dp
-    ros_C(5) =-0.9594562251023355E+01_dp
-    ros_C(6) =-0.2047028614809616E+02_dp
-    ros_C(7) = 0.7496443313967647E+01_dp
-    ros_C(8) =-0.1024680431464352E+02_dp
-    ros_C(9) =-0.3399990352819905E+02_dp
-    ros_C(10) = 0.1170890893206160E+02_dp
-    ros_C(11) = 0.8083246795921522E+01_dp
-    ros_C(12) =-0.7981132988064893E+01_dp
-    ros_C(13) =-0.3152159432874371E+02_dp
-    ros_C(14) = 0.1631930543123136E+02_dp
-    ros_C(15) =-0.6058818238834054E+01_dp
-
-!~~~> M_i = Coefficients for new step solution
-    ros_M(1) = ros_A(7)
-    ros_M(2) = ros_A(8)
-    ros_M(3) = ros_A(9)
-    ros_M(4) = ros_A(10)
-    ros_M(5) = 1.0E+00_dp
-    ros_M(6) = 1.0E+00_dp
-
-!~~~> E_i  = Coefficients for error estimator
-    ros_E(1) = 0.0E+00_dp
-    ros_E(2) = 0.0E+00_dp
-    ros_E(3) = 0.0E+00_dp
-    ros_E(4) = 0.0E+00_dp
-    ros_E(5) = 0.0E+00_dp
-    ros_E(6) = 1.0E+00_dp
-
-!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
-!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
-    ros_NewF(1) = .TRUE.
-    ros_NewF(2) = .TRUE.
-    ros_NewF(3) = .TRUE.
-    ros_NewF(4) = .TRUE.
-    ros_NewF(5) = .TRUE.
-    ros_NewF(6) = .TRUE.
-
-!~~~> ros_ELO  = estimator of local order - the minimum between the
-!        main and the embedded scheme orders plus 1
-    ros_ELO = 4.0_dp
-
-  END SUBROUTINE Rodas4
-
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-!   End of the set of internal Rosenbrock subroutines
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-END SUBROUTINE Rosenbrock
