Index: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/bin/kpp4palm.ksh
===================================================================
--- palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/bin/kpp4palm.ksh	(revision 2768)
+++ palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/bin/kpp4palm.ksh	(revision 3298)
@@ -21,4 +21,15 @@
 # Copyright 2017-2018  Leibniz Universitaet Hannover
 #------------------------------------------------------------------------------#
+# Nov. 2016: Initial Version of KPP chemistry convertor adapted for PALM 
+# by Klaus Ketelsen
+#
+# This code is a modified 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).
+#
+#------------------------------------------------------------------------------#
 #
 # Current revisions:
@@ -29,7 +40,28 @@
 # -----------------
 # $Id$ 
+# forkel   25. September 2018: Added cat for $MECH to pass mechanism name to kpp4palm
+# ketelsen 18. September 2018: Added cat for '#INLINE F90_GLOBAL' 
+# (moved here from mechanisms/def_MECH/chem_gasphase.kpp
+#
+# forkel: 14. September  2018: WCOPY removed
+# ketelsen: July 2018: Adaptations for vektor mode
+# forkel June 2018: re-established original case of subroutine names
+# forkel May 2018: additional copying of chem_gasphase_mod.f90 into $DEFDIR
+# forkel 20.04.2018: removed  wlamch and wlamch_add from $KPP_SUBROUTINE_LIST
+#                    (epsilon(one) is used now)
+# forkel March 2017
+# Re-introduced relative path for KPP_HOME
+# Subroutine list adapted to lowercase subroutine names
+# Added arr2, removed update_sun and k_3rd from subroutine list
+# Renamed output file to chem_gasphase_mod
+# Renamed this file from kp4/ksh to kpp4kpp.ksh
+# changed location of def_mechanism directories to gasphase_preproc/mechanisms
+# 
+# 
+# 2718 2018-01-02 08:49:38Z maronga
 # Initial revision
 # 
 # 
+##########################################################################
 #
 #
@@ -63,4 +95,5 @@
 # Default
 
+MECH=smog
 OUTDIR=`pwd`/../../../SOURCE
 OUTFILE=chem_gasphase_mod
@@ -70,6 +103,7 @@
 VLEN=1
 KEEP="NO"
-DE_INDEX="NO"
-DE_INDEX_FAST="NO"
+UPDT="NO"
+DE_INDEX=0
+DE_INDEX_FAST="YES"
 
 export KPP_SOLVER=Rosenbrock
@@ -77,10 +111,9 @@
 # get Command line option
 
-echo xxxxxxxxxx
-while  getopts :d:ifkp:o:s:v:w:  c     # get options
+while  getopts :m:i:fkup:o:s:vl:w:  c     # get options
 do case $c in
-      d)   DEFDIR=$OPTARG;;          # directory of definition files
-
-      i)   DE_INDEX="YES";;          # if set, deindexing
+      m)   MECH=$OPTARG;;            # mechanism
+
+      i)   DE_INDEX=$OPTARG;;        # if set, deindexing
 
       f)   DE_INDEX_FAST="YES";;     # if set, fast deindexing
@@ -92,13 +125,16 @@
       p)   PREFIX=$OPTARG;;          # Name Prefix
 
-      s)   KPP_SOLVER=$OPTARG;;      # Name Prefix
-
-      v)   MODE="vector"
-           VLEN=$OPTARG;;            # Set to vector Mode
+      s)   KPP_SOLVER=$OPTARG;;      # Update sample f90 code in the def_MECH directory
+
+      u)   UPDT="YES";;              # keep Working directory
+
+      v)   MODE="vector";;           # Set to vector Mode
+
+      l)   VLEN=$OPTARG;;            # Set vector length
 
       w)   WORK=$OPTARG;;            # Working directory
 
       \?)  print ${0##*/} "unknown option:" $OPTARG
-           print "USAGE: ${0##*/} [ -d dir -e -k -o dir -p name -s solver -v length -w dir ] "
+           print "USAGE: ${0##*/} [ -m dir -e -k -u -o dir -p name -s solver -v -l length -w dir ] "
            exit 1;;
    esac
@@ -106,7 +142,14 @@
 shift OPTIND-1
 
-echo $DEFDIR
+echo MECHANISM = $MECH
+echo DE_INDEX = $DE_INDEX
+echo KEEP = $KEEP
+echo UPDT = $UPDT
+echo MODE = $MODE
+echo VLEN = $VLEN
 
 DEF_PREFIX=${PREFIX}.kpp
+DEFDIR=`pwd`/mechanisms/def_$MECH
+echo DEFDIR = $DEFDIR
 
 # Create or clean working directory
@@ -122,9 +165,9 @@
 
 
-KPP_SUBROUTINE_LIST="initialize"
-KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST integrate fun"
-KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST kppsolve kppdecomp wlamch wlamch_add"
-KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST jac_sp k_arr "
-KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST update_rconst arr2"
+KPP_SUBROUTINE_LIST="Initialize"
+KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST INTEGRATE Fun"
+KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST KppSolve KppDecomp"
+KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST Jac_SP k_arr "
+KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST Update_RCONST ARR2"
 KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST initialize_kpp_ctrl error_output"
 
@@ -149,5 +192,9 @@
 
 # Interface ignore list
-KPP_INTERFACE_IGNORE="waxpy wcopy"
+KPP_INTERFACE_IGNORE=" "
+
+echo " "
+echo KPP_SOLVER $KPP_SOLVER
+echo " "
 
 case $KPP_SOLVER in
@@ -155,9 +202,12 @@
 
     Rosenbrock)   
-      KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST wcopy wscal waxpy"
-    if [[ $MODE != "vector" ]]
-    then
-      KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST rosenbrock  funtemplate jactemplate"
-    fi;;
+      KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST WSCAL WAXPY"
+      if [[ $MODE != "vector" ]]
+      then
+         KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST Rosenbrock  FunTemplate JacTemplate"
+        KPP_INTERFACE_IGNORE="WAXPY"
+      else 
+        KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST FunTemplate JacTemplate"
+      fi;;
 
     rosenbrock_mz)
@@ -167,5 +217,5 @@
     rosenbrock)
       KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST WCOPY WSCAL WAXPY"
-      KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST rosenbrock  funtemplate jactemplate";;
+      KPP_SUBROUTINE_LIST="$KPP_SUBROUTINE_LIST Rosenbrock  FunTemplate JacTemplate";;
 
     kpp_lsode)
@@ -203,4 +253,23 @@
 cp $DEFDIR/${PREFIX}.kpp     .
 
+# Global variable are defined here 
+# This has the advantage that it is not necessary to include these variables in all .kpp definition files
+
+cat  >> ${PREFIX}.kpp  <<  EOF
+#INLINE F90_GLOBAL
+! QVAP - Water vapor
+  REAL(kind=dp) :: QVAP
+! FAKT - Conversion factor
+  REAL(kind=dp) :: FAKT
+
+#ENDINLINE
+EOF
+
+# Store mechanism name in file mech_list
+cat  >> mech_list  <<  EOF
+!   Mechanism: $MECH
+!
+EOF
+
 # Run kpp
 
@@ -240,6 +309,14 @@
 done
 
+echo start kp4.exe with arguments
+echo $PREFIX $MODE $VLEN $DE_INDEX $DE_INDEX_FAST
+
 $BASE/bin/kpp4palm.exe $PREFIX $MODE $VLEN $DE_INDEX $DE_INDEX_FAST
 
+#Prelimanary, substitution has to be moved into kpp4palm.exe
+if [[ $MODE = "vector" ]]
+then
+  sed -i -e 's/phot(nphot/phot(vl_dim,nphot/g' kk_kpp.f90
+fi
 
 if [[ -e $OUTDIR/${OUTFILE}.f90 ]] 
@@ -248,8 +325,13 @@
 fi
 cp -p kk_kpp.f90    $OUTDIR/${OUTFILE}.f90
-#cp -p kk_kpp.f90    $MY_PWD/../SOURCE/${OUTFILE}.f90
-
 echo " "
 echo "Write kpp module -- > " $OUTDIR/${OUTFILE}.f90
+
+if [[ $UPDT = "YES" ]]
+then
+cp -p kk_kpp.f90    $DEFDIR/${OUTFILE}.f90
+echo " "
+echo "Write kpp module -- > " $DEFDIR/${OUTFILE}.f90
+fi
 
 if [[ $KEEP = "NO" ]]
Index: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/create_kpp_module.C
===================================================================
--- palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/create_kpp_module.C	(revision 2768)
+++ palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/create_kpp_module.C	(revision 3298)
@@ -17,7 +17,24 @@
 //-----------------
 //$Id: create_kpp_module.C 2470 2017-09-14 13:56:42Z forkel $
+// forkel 25.09.2018:  added  automatic line with meechanism name (read mech_list)
+// forkel 20.09.2018:  added  vl_glo = size(tempi,1)
+// ketelsen 18.09.2018: Removed creation of fill_ Subroutine and creation of calls thereof
+//
+// forkel 12.09.2018: Fix in order not to loose the values of qvap and fakt
+//
+// forkel 03.09.2018: Bug fix: moved kppi.add_line("    CALL initialize after fakt = fakti(is)
+// forkel June 2018: added qvap and fakt, re-established original uppercase/lowercase 
+//                    Deleted definition of qvap,fakt in create_kpp_integrate again
+//
+// ketelsen July 2018: Changes for vector mode (edit_WAXPY, edit_FunTemplate, edit_JacTemplate, 
+//                            some cleanup of comments, various changes in create_kpp_integrate)
+//
+// forkel June 2018: Added qvap and fakt, re-established original uppercase/lowercase 
+//
+//
+/ 2017-09-14 13:56:42Z forkel $
 // Removed preprocessor directive __chem again 
 //
-/ 2017-09-14 13:56:42Z forkel $
+// 2017-09-14 13:56:42Z forkel $
 //
 // 
@@ -25,5 +42,5 @@
 //change of some output to lowercase with uppercase Fortran
 //
-// Nov 2016: Intial version (Klaus Ketelsen)
+// Nov 2016: Intial version of KP4 adapted to PALM (Klaus Ketelsen)
 //
 
@@ -49,15 +66,13 @@
    cout << "Create " << module_name << " from kpp Fortran sources" <<endl;
    cout << "Vector mode " << kpp_switches.is_vector() <<endl;
+   cout << "De_indexing " << kpp_switches.de_indexing() <<endl;
 
    create_fortran_files_and_read();
-   cout << "## after create_fortran_files_and_read " <<endl;
 
 // Generate first module lines 
 
      string first_line="MODULE " + module_name;
-     mz_kpp.add_line(first_line);
-     mz_kpp.add_line("!");
-//   mz_kpp.add_line("#if defined( __chem )");
-//   mz_kpp.add_line(" ");
+   mz_kpp.add_line(first_line);
+   mz_kpp.add_line(" ");
 
 //    string e5_line = first_line +"_e5";
@@ -93,5 +108,4 @@
      it->edit_fortran ();
    }
-   cout << "## after edit FORTRAN files            " <<endl;
 
 // Generate a list of single subroutines from kpp-files
@@ -109,7 +123,7 @@
    }
 
-
-   if(kpp_switches.is_vector()) {
-
+   if(kpp_switches.is_vector()) {
+
+   cout << "##### Hier kpp_switches.is_vector          " <<endl;
 //   Change header section
      for(it=kpp_includes.begin();it!=kpp_includes.end();it++) {
@@ -138,6 +152,15 @@
          it->edit_Fun();
        }
-     }
-   cout << "## after Edit individual subroutines    " <<endl;
+       if(it->get_name() == "WAXPY" ) {
+         it->edit_WAXPY();
+       }
+       if(it->get_name() == "FunTemplate" ) {
+         it->edit_FunTemplate();
+       }
+       if(it->get_name() == "JacTemplate" ) {
+         it->edit_JacTemplate();
+       }
+     }
+// cout << "## after Edit individual subroutines    " <<endl;
 
    }
@@ -146,7 +169,12 @@
 
    for(it=kpp_subroutines.begin();it!=kpp_subroutines.end();it++) {
-     if(it->get_name() == "update_rconst") {
+     if(it->get_name() == "Update_RCONST") {
        it->edit_Update_RCONST(Vvar_list);
      }
+
+     if(it->get_name() == "Initialize") {
+       it->edit_Initialize(Vvar_list);
+     }
+
    }
 
@@ -184,32 +212,15 @@
      }
 
-     buf = prefix + "INTERFACE            " + it->get_name() ;
+     buf = prefix + "interface            " + it->get_name() ;
      mz_kpp.add_line(buf);
-     buf = prefix + "  MODULE PROCEDURE   " + it->get_name();
+     buf = prefix + "  module procedure   " + it->get_name();
      mz_kpp.add_line(buf);
-     buf = prefix + "END INTERFACE        " + it->get_name();
+     buf = prefix + "end interface        " + it->get_name();
      mz_kpp.add_line(buf);
      mz_kpp.add_line(" ");
    }
-   for(iv=Vvar_list.begin();iv!=Vvar_list.end();iv++) {
-     string          buf;
-
-     string sub_name = "fill_" + iv->name;
-     buf = "  INTERFACE            " + sub_name;
-     mz_kpp.add_line(buf);
-     buf = "    MODULE PROCEDURE   " + sub_name;
-     mz_kpp.add_line(buf);
-     buf = "  END INTERFACE        " + sub_name;
-     mz_kpp.add_line(buf);
-     buf = "  PUBLIC               " + sub_name;
-     mz_kpp.add_line(buf);
-     mz_kpp.add_line(" ");
-   }
 
    mz_kpp.add_line(" ");
 
-   for(iv=Vvar_list.begin();iv!=Vvar_list.end();iv++) {
-     create_fill_routine(kpp_subroutines, *iv);
-   }
 
 // Copy FORTRAN subroutines to mz_kpp
@@ -224,19 +235,7 @@
 // Finish module
 
-//   mz_kpp.add_line("#endif");
-     string last_line="END MODULE " + module_name;
-     mz_kpp.add_line("");
-     mz_kpp.add_line(last_line);
-
-// Handle e5 module
-
-//    for(it=e5_subroutines.begin();it!=e5_subroutines.end();it++) {
-//      e5_kpp.add_line(" ");
-//      it->copy_to_MZ_KPP(e5_kpp);
-//    }
-
-//    last_line = last_line + "_e5";
-//    e5_kpp.add_line(" ");
-//    e5_kpp.add_line(last_line);
+   string last_line="end module " + module_name;
+   mz_kpp.add_line("");
+   mz_kpp.add_line(last_line);
 
 // Write the complete module to file: mz_kpp.f
@@ -376,4 +375,21 @@
    string                          diline;
 
+// Read mechanism from mech_list
+
+   in.open("mech_list");
+   if( !in ) {
+      cout << "cannot open " << endl; my_abort("mech_list");
+   }
+
+   while ( 1 ) {
+     getline (in, buf);
+     if( in.eof() ) break;
+     if( in.bad() ) my_abort("ERROR_READ_4");
+     line.set_line(buf);
+     mz_kpp.add_line(line);
+   }
+   in.close();
+
+
 // Read Modul Header from file $MZ_KPP_HOME/templates/module_header
 
@@ -394,40 +410,27 @@
    mz_kpp.add_line("                                                                 "); 
    if(kpp_switches.is_vector()) {
-       mz_kpp.add_line("  LOGICAL,  PARAMETER          :: l_vector = .TRUE.             ");
-   } else {
-       mz_kpp.add_line("  LOGICAL,  PARAMETER          :: l_vector = .FALSE.            ");
+       mz_kpp.add_line("  logical,parameter          :: L_VECTOR = .TRUE.             ");
+   } else {
+       mz_kpp.add_line("  logical,parameter          :: L_VECTOR = .FALSE.            ");
    }
 //  mz_pj_20070531+
    sprintf(distr,"%i",kpp_switches.de_indexing());
    diline = distr ;
-   mz_kpp.add_line("  INTEGER,  PARAMETER          :: i_lu_di = " + diline );
+   mz_kpp.add_line("  integer,parameter          :: I_LU_DI = " + diline );
 //  mz_pj_20070531-
 
-   mz_kpp.add_line("  INTEGER,  PARAMETER          :: vl_dim = " 
-                   + kpp_switches.get_vector_length() ); 
-   mz_kpp.add_line("  INTEGER                     :: vl                              "); 
+   mz_kpp.add_line("  integer,parameter          :: VL_DIM = " 
+                 + kpp_switches.get_vector_length() ); 
+   mz_kpp.add_line("  integer                     :: vl                              "); 
    mz_kpp.add_line("                                                                 "); 
-   mz_kpp.add_line("  INTEGER                     :: vl_glo                          "); 
-   mz_kpp.add_line("  INTEGER                     :: is,ie                           "); 
+   mz_kpp.add_line("  integer                     :: VL_glo                          "); 
+   mz_kpp.add_line("  integer                     :: is,ie                           "); 
    mz_kpp.add_line("                                                                 "); 
-   mz_kpp.add_line("  INTEGER,  DIMENSION(vl_dim)  :: kacc, krej                      "); 
-   mz_kpp.add_line("  INTEGER,  DIMENSION(vl_dim)  :: ierrv                           "); 
-   mz_kpp.add_line("  LOGICAL                     :: data_loaded = .false.             "); 
-   
+   mz_kpp.add_line("                                                                 "); 
+   mz_kpp.add_line("  integer, dimension(VL_dim)   :: Kacc,Krej                       "); 
+   mz_kpp.add_line("  integer, dimension(VL_dim)   :: IERRV                           "); 
+   mz_kpp.add_line("  logical                     :: data_loaded = .false.             "); 
+
    in.close();
-
-//    in_e5.open("module_header_e5");
-//    if( !in_e5 ) {
-//       cout << "cannot open " << endl; my_abort("module_header_e5");
-//    }
-
-//    while ( 1 ) {
-//      getline (in_e5, buf);
-//      if( in_e5.eof() ) break;
-//      if( in_e5.bad() ) my_abort("ERROR_READ_4");
-//      line.set_line(buf);
-//      e5_kpp.add_line(line);
-//    }
-//    in_e5.close();
 
    return;
@@ -448,13 +451,4 @@
    out.close();
    
-//    out_file  = "kk_mecca_kpp_e5.f90";
-//    out_e5.open(out_file.c_str(), ios::out);
-//    if( !out_e5 ) {
-//       cout << "cannot open " << endl; my_abort(out_file);
-//    }
-
-//    e5_kpp.write_file (out_e5);
-
-//    out_e5.close();
 
    return;
@@ -468,48 +462,73 @@
    kppi.set_name("chem_gasphase_integrate");
 
-   kppi.add_line("SUBROUTINE chem_gasphase_integrate (time_step_len, conc, tempk, photo, ierrf, xnacc, xnrej, istatus, l_debug, pe )  ");
+   kppi.add_line("SUBROUTINE chem_gasphase_integrate (time_step_len, conc, tempi, qvapi, fakti, photo, ierrf, xnacc, xnrej, istatus, l_debug, pe, icntrl_i, rcntrl_i )  ");
    kppi.add_line("                                                                    ");
    kppi.add_line("  IMPLICIT NONE                                                     ");
    kppi.add_line("                                                                    ");
 
-   kppi.add_line("  REAL(dp),  INTENT(IN)                    :: time_step_len           ");
+   kppi.add_line("  REAL(dp), INTENT(IN)                   :: time_step_len           ");
    kppi.add_line("  REAL(dp),  DIMENSION(:,:),  INTENT(INOUT) :: conc                    ");
-   kppi.add_line("  REAL(dp),  DIMENSION(:,:),  INTENT(INOUT) :: photo                   ");
-   kppi.add_line("  REAL(dp),  DIMENSION(:),  INTENT(IN)      :: tempk                   ");
-   kppi.add_line("  INTEGER,  INTENT(OUT),  OPTIONAL          :: ierrf(:)                ");
-   kppi.add_line("  INTEGER,  INTENT(OUT),  OPTIONAL          :: xnacc(:)                ");
-   kppi.add_line("  INTEGER,  INTENT(OUT),  OPTIONAL          :: xnrej(:)                ");
-   kppi.add_line("  INTEGER,  INTENT(INOUT),  OPTIONAL        :: istatus(:)              ");
-   kppi.add_line("  INTEGER,  INTENT(IN),  OPTIONAL           :: pe                      ");
-   kppi.add_line("  LOGICAL,  INTENT(IN),  OPTIONAL           :: l_debug                 ");
+   kppi.add_line("  REAL(dp),  DIMENSION(:,:),  INTENT(IN)    :: photo                   ");
+   kppi.add_line("  REAL(dp),  DIMENSION(:),  INTENT(IN)      :: tempi                   ");
+   kppi.add_line("  REAL(dp),  DIMENSION(:),  INTENT(IN)      :: qvapi                   ");
+   kppi.add_line("  REAL(dp),  DIMENSION(:),  INTENT(IN)      :: fakti                   ");
+   kppi.add_line("  INTEGER,  INTENT(OUT), OPTIONAL        :: ierrf(:)                ");
+   kppi.add_line("  INTEGER,  INTENT(OUT), OPTIONAL        :: xNacc(:)                ");
+   kppi.add_line("  INTEGER,  INTENT(OUT), OPTIONAL        :: xNrej(:)                ");
+   kppi.add_line("  INTEGER,  INTENT(INOUT), OPTIONAL      :: istatus(:)              ");
+   kppi.add_line("  INTEGER,  INTENT(IN), OPTIONAL         :: PE                      ");
+   kppi.add_line("  LOGICAL,  INTENT(IN), OPTIONAL         :: l_debug                 ");
+   kppi.add_line("  INTEGER,  DIMENSION(nkppctrl),INTENT(IN), OPTIONAL  :: icntrl_i         ");
+   kppi.add_line("  REAL(dp), DIMENSION(nkppctrl),INTENT(IN), OPTIONAL  :: rcntrl_i         ");
    kppi.add_line("                                                                    ");
    kppi.add_line("  INTEGER                                 :: k   ! loop variable     ");
    kppi.add_line("  REAL(dp)                                :: dt                      ");
-   kppi.add_line("  INTEGER,  DIMENSION(20)                  :: istatus_u               ");
-   kppi.add_line("  INTEGER                                 :: ierr_u                  ");
-   kppi.add_line("                                                                    ");
-   kppi.add_line("                                                                    ");
-   kppi.add_line("  if (present (istatus) )  istatus = 0                              ");
-   kppi.add_line("                                                                    ");
-// kppi.add_line("  vk_glo = size(tempk,1)                                            ");
-// kppi.add_line("                                                                    ");
-
-   kppi.add_line("  DO k=1,vl_glo,vl_dim                                              ");
+   kppi.add_line("  integer, dimension(20)                 :: istatus_u               ");
+   kppi.add_line("  integer                                :: ierr_u                  ");
+   kppi.add_line("  integer                                :: istatf                  ");
+   kppi.add_line("  integer                                :: vl_dim_lo               ");
+   kppi.add_line("                                                                    ");
+   kppi.add_line("                                                                    ");
+   kppi.add_line("  if (present (istatus) )   istatus = 0                             ");
+   kppi.add_line("  if (present (icntrl_i) )  icntrl  = icntrl_i                      ");
+   kppi.add_line("  if (present (rcntrl_i) )  rcntrl  = rcntrl_i                      ");
+   kppi.add_line("                                                                    ");
+   kppi.add_line("  vl_glo = size(tempi,1)                                            ");
+   kppi.add_line("                                                                    ");
+   kppi.add_line("  vl_dim_lo = VL_DIM                                                ");
+   kppi.add_line("  DO k=1,VL_glo,vl_dim_lo                                           ");
    kppi.add_line("    is = k                                                          ");
-   kppi.add_line("    ie = min(k+vl_dim-1,vl_glo)                                     ");
+   kppi.add_line("    ie = min(k+vl_dim_lo-1,VL_glo)                                  ");
    kppi.add_line("    vl = ie-is+1                                                    ");
 
    kppi.add_line("                                                                    ");
    if(kpp_switches.is_vector()) {
-     kppi.add_line("    c(1:vl,:) = conc(is:ie,:)                                     ");
-   } else {
-     kppi.add_line("    c(:) = conc(is,:)                                             ");
-   }
-   kppi.add_line("                                                                    ");
-   if(kpp_switches.is_vector()) {
-     kppi.add_line("    temp(1:vl) = tempk(is:ie)                                     ");
-   } else {
-     kppi.add_line("    temp = tempk(is)                                              ");
-   }
+     kppi.add_line("    C(1:vl,:) = Conc(is:ie,:)                                     ");
+   } else {
+     kppi.add_line("    C(:) = Conc(is,:)                                             ");
+   }
+
+   kppi.add_line("                                                                    ");
+   if(kpp_switches.is_vector()) {
+     kppi.add_line("    temp(1:vl) = tempi(is:ie)                                     ");
+   } else {
+     kppi.add_line("    temp = tempi(is)                                              ");
+   }
+   kppi.add_line("                                                                    ");
+   if(kpp_switches.is_vector()) {
+     kppi.add_line("    qvap(1:vl) = qvapi(is:ie)                                     ");
+   } else {
+     kppi.add_line("    qvap = qvapi(is)                                              ");
+   }
+   kppi.add_line("                                                                    ");
+   if(kpp_switches.is_vector()) {
+     kppi.add_line("    fakt(1:vl) = fakti(is:ie)                                     ");
+   } else {
+     kppi.add_line("    fakt = fakti(is)                                              ");
+   }
+
+   kppi.add_line("                                                                    ");
+   kppi.add_line("    CALL initialize                                                 ");
+
    kppi.add_line("                                                                    ");
    if(kpp_switches.is_vector()) {
@@ -528,10 +547,11 @@
    kppi.add_line("                                                                    ");
    if(kpp_switches.is_vector()) {
-     kppi.add_line("    conc(is:ie,:) = c(1:VL,:)                                     ");
-   } else {
-     kppi.add_line("    IF (PRESENT(l_debug) .AND. PRESENT(pe)) THEN                       ");
-     kppi.add_line("       IF (l_debug) CALL error_output(conc(is,:),ierr_u, pe)           ");
-     kppi.add_line("    ENDIF                                                              ");
-     kppi.add_line("    conc(is,:) = c(:)                                                  ");
+     kppi.add_line("    Conc(is:ie,:) = C(1:VL,:)                                     ");
+   } else {
+     kppi.add_line("   IF (PRESENT(l_debug) .AND. PRESENT(PE)) THEN                       ");
+     kppi.add_line("      IF (l_debug) CALL error_output(Conc(is,:),ierr_u, PE)           ");
+     kppi.add_line("   ENDIF                                                              ");
+     kppi.add_line("                                                                      ");
+     kppi.add_line("    Conc(is,:) = C(:)                                                 ");
    }
 
@@ -540,17 +560,17 @@
    kppi.add_line("                                                                    ");
    if(kpp_switches.is_vector()) {
-     kppi.add_line("    if(PRESENT(ierrf))    ierrf(is:ie) = ierrv(1:vl)              ");
-     kppi.add_line("    if(PRESENT(xnacc))    xnacc(is:ie) = kacc(1:vl)               ");
-     kppi.add_line("    if(PRESENT(xnrej))    xnrej(is:ie) = krej(1:vl)               ");
-   } else {
-     kppi.add_line("    if(PRESENT(ierrf))    ierrf(is) = ierr_u                      ");
-     kppi.add_line("    if(PRESENT(xnacc))    xnacc(is) = istatus_u(4)                ");
-     kppi.add_line("    if(PRESENT(xnrej))    xnrej(is) = istatus_u(5)                ");
-   }
-   kppi.add_line("                                                                    ");
-   kppi.add_line("    if (PRESENT (istatus) )  then                                   ");
-   if(kpp_switches.is_vector()) {
-     kppi.add_line("      istatus(4) =   istatus(4) + sum(kacc(1:vl))                  ");
-     kppi.add_line("      istatus(5) =   istatus(5) + sum(krej(1:vl))                  ");
+     kppi.add_line("    if(Present(ierrf))    ierrf(is:ie) = IERRV(1:VL)              ");
+     kppi.add_line("    if(Present(xNacc))    xNacc(is:ie) = Kacc(1:VL)               ");
+     kppi.add_line("    if(Present(xNrej))    xNrej(is:ie) = Krej(1:VL)               ");
+   } else {
+     kppi.add_line("    if(Present(ierrf))    ierrf(is) = IERR_U                      ");
+     kppi.add_line("    if(Present(xNacc))    xNacc(is) = istatus_u(4)                ");
+     kppi.add_line("    if(Present(xNrej))    xNrej(is) = istatus_u(5)                ");
+   }
+   kppi.add_line("                                                                    ");
+   kppi.add_line("    if (present (istatus) )  then                                   ");
+   if(kpp_switches.is_vector()) {
+     kppi.add_line("      istatus(4) =   istatus(4) + sum(Kacc(1:VL))                  ");
+     kppi.add_line("      istatus(5) =   istatus(5) + sum(Krej(1:VL))                  ");
      kppi.add_line("      istatus(3) =   istatus(4) + istatus(5)                       ");
      kppi.add_line("      istatus(6) =   istatus(6) + istatus_u(6)                     ");
@@ -559,5 +579,5 @@
      kppi.add_line("      istatus(1:8) = istatus(1:8) + istatus_u(1:8)                 ");
    }
-   kppi.add_line("    ENDIF                                                          ");
+   kppi.add_line("    end if                                                          ");
    kppi.add_line("                                                                    ");
    kppi.add_line("  END DO                                                            ");
@@ -568,5 +588,5 @@
    kppi.add_line("                                                                    ");
    for(iv=Vvar_list.begin();iv!=Vvar_list.end();iv++) {
-     kppi.add_line("  if (ALLOCATED("+ iv->name +"))   DEALLOCATE("+ iv->name +" )    ");
+//     kppi.add_line("  if (allocated("+ iv->name +"))   deallocate("+ iv->name +" )    ");
    }
 
@@ -574,5 +594,5 @@
    kppi.add_line("  data_loaded = .false.                                             ");
    kppi.add_line("                                                                    ");
-   kppi.add_line("  RETURN                                                            ");
+   kppi.add_line("  return                                                            ");
    kppi.add_line("END SUBROUTINE chem_gasphase_integrate                                        ");
 
@@ -583,60 +603,2 @@
 }
 
-void create_kpp_module::create_fill_routine(vector<fortran_file> &fi_list, Vvar & var) {
-   fortran_file                         fi;
-   vector<string>::iterator             is;
-   string                               line;
-
-   cout << "Generate fill subroutine for " << var.name << endl;
-
-   fi.set_name(var.name);
-   line = "  SUBROUTINE fill_" + var.name;
-   fi.add_line(line + "(status, array) ");
-     fi.add_line(" ");
-     fi.add_line("    INTEGER,  INTENT(OUT)               :: status ");
-   if(var.nr_dim() == 0) {
-     fi.add_line("    REAL(dp),  INTENT(IN),  DIMENSION(:) :: array ");
-     fi.add_line(" ");
-     fi.add_line("    status = 0");
-     fi.add_line("    IF (.not. ALLOCATED("+var.name+")) & ");
-     fi.add_line("       ALLOCATE("+var.name+"(size(array))) ");
-   } else if(var.nr_dim() == 1) {
-     fi.add_line("    REAL (dp), INTENT(IN), DIMENSION(:,:) :: array ");
-     fi.add_line(" ");
-     fi.add_line("    status = 0 ");
-     fi.add_line("    if (.not. ALLOCATED("+var.name+")) & ");
-     fi.add_line("        ALLOCATE("+var.name+"(size(array,1),"+var.dim_var[0]+")) ");
-   } else if(var.nr_dim() == 2) {
-     fi.add_line(" ");
-     fi.add_line("    REAL (dp), INTENT(IN), DIMENSION(:,:,:) :: array ");
-     fi.add_line(" ");
-     fi.add_line("    status = 0 ");
-     fi.add_line("    if (.not. ALLOCATED("+var.name+")) & ");
-     fi.add_line("        ALLOCATE("+var.name+"(size(array,1),"+var.dim_var[0]+var.dim_var[1]+")) ");
-   } else {
-     fi.add_line("    REAL (dp), INTENT(IN), DIMENSION(:,:,:,:) :: array ");
-     fi.add_line(" ");
-     fi.add_line("    status = 0 ");
-     fi.add_line("    IF (.not. ALLOCATED("+var.name+")) & ");
-     fi.add_line("        ALLOCATE("+var.name+"(size(array,1),"+var.dim_var[0]
-                               +var.dim_var[1]+var.dim_var[3]+")) ");
-   }
-
-   fi.add_line(" ");
-   fi.add_line("    IF (data_loaded .AND. (vl_glo /= size(array,1)) )  THEN ");
-   fi.add_line("       status = 1 ");
-   fi.add_line("       RETURN ");
-   fi.add_line("    END IF ");
-   fi.add_line(" ");
-   fi.add_line("    vl_glo = size(array,1) ");
-   fi.add_line("    "+var.name+ " = array ");
-   fi.add_line("    data_loaded = .TRUE. ");
-   fi.add_line(" ");
-   fi.add_line("    RETURN");
-   fi.add_line(" ");
-   fi.add_line("  END " + line);
-
-   fi_list.push_back(fi);
-
-   return;
-}
Index: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/create_kpp_module.h
===================================================================
--- palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/create_kpp_module.h	(revision 2768)
+++ palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/create_kpp_module.h	(revision 3298)
@@ -1,14 +1,20 @@
 #ifndef mecca
-#define mecca 1
+// ketelsen 18.09.2018: Removed create_fill_routine
 
-// ############################################################################
-//
-//     create_mz_kpp_module 
-//
-//     create scalar code from .f90 sources created by KPP to be used in MESSy
-//
-//     COPYRIGHT Klaus Ketelsen and MPI-CH   April 2007
-//
-// ############################################################################
+#include <iostream>
+#include <fstream>
+
+#include <string>
+#include <list>
+#include <vector>
+
+#include "fortran_file.h"
+#include "expand_decomp.h"
+
+// Class to create module which contains code generated by kpp.
+
+/ 
+/ 
+// ketelsen 18.09.2018: Removed create_fill_routine
 
 #include <iostream>
@@ -48,6 +54,6 @@
 
   void create_kpp_integrate();
-  void create_fill_routine(vector<fortran_file> &e5_list, Vvar &var );
- public:
+
+public:
 
   void do_work (string s) ;
Index: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/expand_decomp.C
===================================================================
--- palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/expand_decomp.C	(revision 2768)
+++ palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/expand_decomp.C	(revision 3298)
@@ -165,5 +165,5 @@
   vector<fortran_file>::iterator  it;
   fortran_file                    de;
-  int                             k,kk,j,jj,i,ii;
+  int                             k,kk,j,jj,i;
   string                          line;
   char                            cline[80];
Index: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file.C
===================================================================
--- palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file.C	(revision 2768)
+++ palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file.C	(revision 3298)
@@ -13,9 +13,19 @@
 //Current revisions:
 //------------------
-// 
-// 
+/ 
+/ 
 //Former revisions:
 //-----------------
 //$Id: fortran_file.C 2470 2017-09-14 13:56:42Z forkel $
+//  ketelsen 18.09.2018: (Rev 3260) added vector switch and creation of dimension statement
+//  forkel Sept.2018:  'fill_' routine applied to qvap and fakt
+//                      removal of unnecessary variables (LOOKAT, monitor etc. from Fortran code)
+//  forkel June 2018:  Moved adaption to PALM conventions to the end of the process
+//                     in order to make future use of vector code adaptations possible
+//  forkel 20.04.2018: Replace Roundoff = WLAMCH('E') by Roundoff = epsilon(one)
+// 
+// 
+// 
+/ 2017-09-14 13:56:42Z forkel $
 //
 // 
@@ -26,5 +36,5 @@
 //
 //
-// Nov 2016: Intial version (Klaus Ketelsen)
+// Nov 2016: Intial version of KP4 adapted to PALM (Klaus Ketelsen)
 //
 
@@ -62,7 +72,24 @@
      }
 
-// Now do some cosmetics to adapt the KPP generated output a bit o the looks of PALM, 
-// i.e. add some blanks, convert all to lowercase except Fortran language elements, etc.
-       if(line.find("'",0) == string::npos)  {     // No substitution in line with strings
+     if(line.find("'",0) == string::npos)  {     // No substitution in line with strings
+//     if(line.substr(0,4) !="!KPP") {
+//      global_substitute(line,"!"," ! ");
+//     if(line.substr(0,2) =="! ")   continue;           // No Substitute of Comments
+//     global_substitute(line,"("," ( ");
+//     global_substitute(line,")"," ) ");
+//     global_substitute(line,",",", ");
+//     global_substitute(line,"*","* ");
+//     global_substitute(line,"* *","**");
+//     global_substitute(line,"/JVS","/ JVS");
+//     global_substitute(line,"-","- ");
+//     global_substitute(line,"e- ","e-");
+//     global_substitute(line,"+","+ ");
+//     global_substitute(line,"d- ","d-");
+//     global_substitute(line,"D- ","D-");
+//     global_substitute(line,"e+ ","e+");
+//     global_substitute(line,"E+ ","E+");
+//     global_substitute(line,"E- ","E-");
+//      }
+
         if(line.substr(0,2) =="! ") {
          global_substitute(line,"#","#");      // just a dummy so comments do not get lost
@@ -77,48 +104,7 @@
          global_substitute(line,"/JVS","/ JVS");
          global_substitute(line,"-","- ");
-         global_substitute(line,"e- ","e-");
          global_substitute(line,"+","+ ");
          global_substitute(line,"+  ","+ ");
-         global_substitute(line,"d- ","d-");
-         global_substitute(line,"D- ","D-");
-         global_substitute(line,"e+ ","e+");
-         global_substitute(line,"E+ ","E+");
-         global_substitute(line,"E- ","E-");
-         global_substitute(line,")=",") =");
-         global_subtolower(line);
-         global_substitute(line,"call ","CALL ");
-         global_substitute(line,"case","CASE");
-         global_substitute(line,"character","CHARACTER");
-         global_substitute(line,"contains","CONTAINS");
-         global_substitute(line,"dimension","DIMENSION");
-         global_substitute(line,"do ","DO ");
-         global_substitute(line,"elseif","ELSEIF");
-         global_substitute(line,"else","ELSE");
-         global_substitute(line,"#ELSE","#else");
-         global_substitute(line,"end do","ENDDO");
-         global_substitute(line,"end if","ENDIF");
-         global_substitute(line,"end ","END ");    // Modify "end" after all other strings containing "end..." are done!
-         global_substitute(line,"#ENDIF","#endif");
-         global_substitute(line,"function","FUNCTION");
-         global_substitute(line,"if(","IF (");
-         global_substitute(line," if "," IF ");
-         global_substitute(line,"implicit","IMPLICIT");
-         global_substitute(line,"include","INCLUDE");
-         global_substitute(line,"intent","INTENT");
-         global_substitute(line,"integer","INTEGER");
-         global_substitute(line,"logical","LOGICAL");
-         global_substitute(line,"module","MODULE");
-         global_substitute(line,"none","NONE");
-         global_substitute(line,"optional","OPTIONAL");
-         global_substitute(line,"parameter","PARAMETER");
-         global_substitute(line,"public","PUBLIC");
-         global_substitute(line,"real","REAL");
-         global_substitute(line,"return","RETURN");
-         global_substitute(line,"use ","USE ");
-         global_substitute(line,"save","SAVE");
-         global_substitute(line,"subroutine","SUBROUTINE");
-         global_substitute(line,"then","THEN");
-         global_substitute(line,"while","WHILE");
-       }
+         }
      }
 
@@ -167,6 +153,5 @@
 //  Update_RCONST has only to be called once per outer timeloop in KPP_FOR_PALM
 
-//  if(ip->get_token(0) == "CALL" && ip->get_token(1) == "Update_RCONST" ) {
-    if(ip->get_token(0) == "CALL" && ip->get_token(1) == "update_rconst" ) {
+    if(ip->get_token(0) == "CALL" && ip->get_token(1) == "Update_RCONST" ) {
       lo_line.insert(0,"!DELETE ");
     cout << lo_line << endl;
@@ -175,6 +160,12 @@
 //  Update_SUN must not be called within in KPP_FOR_PALM
 
-//  if(ip->get_token(0) == "CALL" && ip->get_token(1) == "Update_SUN" ) {
-    if(ip->get_token(0) == "CALL" && ip->get_token(1) == "update_sun" ) {
+    if(ip->get_token(0) == "CALL" && ip->get_token(1) == "Update_SUN" ) {
+      lo_line.insert(0,"!DELETE ");
+    cout << lo_line << endl;
+    }
+
+//  Remove "    var(i) = x"  (Make sure  that var does not occur as first token anywhere else)
+
+    if(ip->get_token(0) == "VAR") {
       lo_line.insert(0,"!DELETE ");
     cout << lo_line << endl;
@@ -338,4 +329,93 @@
     }
 
+//  Remove  REAL(kind=dp):: sun
+
+    if(ip->get_token_number_from_string("SUN") > 0) {
+      lo_line.insert(0,"!DELETE ");
+    cout << lo_line << endl;
+    }
+
+//  Remove    INTEGER :: ddmtype
+
+    if(ip->get_token_number_from_string("DDMTYPE") > 0) {
+    cout << lo_line << endl;
+      lo_line.insert(0,"!DELETE ");
+    cout << lo_line << endl;
+    }
+
+//  Remove    REAL(kind=dp) :: dt
+
+    if(ip->get_token_number_from_string("DT") > 0) {
+    cout << lo_line << endl;
+      lo_line.insert(0,"!DELETE ");
+    cout << lo_line << endl;
+    }
+
+//  Remove    LOOKAT
+
+    if(ip->get_token_number_from_string("LOOKAT") > 0) {
+    cout << lo_line << endl;
+      lo_line.insert(0,"!DELETE ");
+    cout << lo_line << endl;
+    }
+
+//  Remove    NLOOKAT
+
+    if(ip->get_token_number_from_string("NLOOKAT") > 0) {
+    cout << lo_line << endl;
+      lo_line.insert(0,"!DELETE ");
+    cout << lo_line << endl;
+    }
+
+//  Remove    MONITOR
+
+    if(ip->get_token_number_from_string("MONITOR") > 0) {
+    cout << lo_line << endl;
+      lo_line.insert(0,"!DELETE ");
+    cout << lo_line << endl;
+    }
+
+//  Remove    NMONITOR
+
+    if(ip->get_token_number_from_string("NMONITOR") > 0) {
+    cout << lo_line << endl;
+      lo_line.insert(0,"!DELETE ");
+    cout << lo_line << endl;
+    }
+
+//  Remove    SMASS
+
+    if(ip->get_token_number_from_string("SMASS") > 0) {
+    cout << lo_line << endl;
+      lo_line.insert(0,"!DELETE ");
+    cout << lo_line << endl;
+    }
+
+//  Remove    RTOLS
+
+    if(ip->get_token_number_from_string("RTOLS") > 0) {
+    cout << lo_line << endl;
+      lo_line.insert(0,"!DELETE ");
+    cout << lo_line << endl;
+    }
+
+//  Remove    TEND
+
+    if(ip->get_token_number_from_string("TEND") > 0) {
+    cout << lo_line << endl;
+      lo_line.insert(0,"!DELETE ");
+    cout << lo_line << endl;
+    }
+
+//  Remove    STEPMAX
+
+    if(ip->get_token_number_from_string("STEPMAX") > 0) {
+    cout << lo_line << endl;
+      lo_line.insert(0,"!DELETE ");
+    cout << lo_line << endl;
+    }
+
+
+
 //  Make ind_ variables public
     if(ip->get_token(3).substr(0,4) == "ind_") {
@@ -464,8 +544,16 @@
        to_do = false;
     }
-    if(ip->get_token(0).substr(0,1) != "!" && ip->get_token_number_from_string("temp") > 2) {
-//  if(ip->get_token(0).substr(0,1) != "!" && ip->get_token_number_from_string("XXXX") > 2) {
-      todo_1 = true;
-    }
+    if(kpp_switches.is_vector()) {
+       if(ip->get_token(0).substr(0,1) != "!" && ip->get_token_number_from_string("TEMP") > 2) {
+         todo_1 = true;
+       }
+       if(ip->get_token(0).substr(0,1) != "!" && ip->get_token_number_from_string("QVAP") > 2) {
+         todo_1 = true;
+       }
+       if(ip->get_token(0).substr(0,1) != "!" && ip->get_token_number_from_string("FAKT") > 2) {
+         todo_1 = true;
+       }
+    }
+
     if(to_do || todo_1) {
       if(ip->get_token(0).substr(0,1) == "!") continue;       // skip comment limes
@@ -499,19 +587,24 @@
       cout << "Vector variable " << ip->get_token(pos+1) <<" " << vari.nr_dim() <<endl;
       lo_line = ip->get_line() ;
-      if(vari.nr_dim() == 0) {
-        lo_line.clear();
-        lo_line = "  REAL(dp),dimension(:),allocatable             :: " + vari.name;
-      }
-      if(vari.nr_dim() == 1) {
-        lo_line.clear();
-        lo_line = "  REAL(dp),dimension(:,:),allocatable           :: " + vari.name;
-      }
-      if(vari.nr_dim() == 2) {
-        lo_line.clear();
-        lo_line = "  REAL(dp),dimension(:,:,:),allocatable         :: " + vari.name;
-      }
-      if(vari.nr_dim() == 3) {
-        lo_line.clear();
-        lo_line = "  REAL(dp),dimension(:,:,:,:),allocatable       :: " + vari.name;
+      if(todo_1)   {
+          lo_line.clear();
+          lo_line = "  REAL(dp),dimension(VL_DIM)                  :: " + vari.name;
+      } else {
+         if(vari.nr_dim() == 0) {
+           lo_line.clear();
+           lo_line = "  REAL(dp),dimension(:),allocatable             :: " + vari.name;
+         }
+         if(vari.nr_dim() == 1) {
+           lo_line.clear();
+           lo_line = "  REAL(dp),dimension(:,:),allocatable           :: " + vari.name;
+         }
+         if(vari.nr_dim() == 2) {
+           lo_line.clear();
+           lo_line = "  REAL(dp),dimension(:,:,:),allocatable         :: " + vari.name;
+         }
+         if(vari.nr_dim() == 3) {
+           lo_line.clear();
+           lo_line = "  REAL(dp),dimension(:,:,:,:),allocatable       :: " + vari.name;
+         }
       }
       ip->set_line(lo_line);
@@ -541,5 +634,7 @@
 
   vector<program_line>::iterator     ip;
+  program_line                    pl;
   string                             lo_line;
+  string                             line;
 
   for (ip=pline.begin(); ip != pline.end(); ip++) {
@@ -555,6 +650,6 @@
     ip->global_substitute("* ","*");
     ip->global_substitute("* ","*");
-    ip->global_substitute("d- ","d-");
-    ip->global_substitute("d+ ","d+");
+//  ip->global_substitute("d- ","d-");
+//  ip->global_substitute("d+ ","d+");
     ip->global_substitute(", ",",");
     ip->global_substitute(")=",") =");
@@ -562,9 +657,108 @@
     ip->global_substitute(", - ",", -");
 
+// Now do some cosmetics to adapt the KPP generated output a bit o the looks of PALM, 
+// i.e. add some blanks, convert all to lowercase except Fortran language elements, etc.
+//  Replace Roundoff = WLAMCH('E') by Roundoff = epsilon(one)
+       ip->global_substitute("Roundoff = WLAMCH('E')","roundoff = epsilon(one)");
+       lo_line = ip->get_line();
+       if(lo_line.find("'",0) == string::npos)  {     // No substitution in line with strings
+        if(lo_line.substr(0,2) =="! ") {
+         ip->global_substitute("#","#");      // just a dummy so comments do not get lost
+         } else {
+//       cout << "HIER0 " << lo_line <<endl;
+//       global_substitute(lo_line,"("," ( ");
+//       global_substitute(lo_line,")"," ) ");
+         global_substitute(lo_line,",",", ");
+         global_substitute(lo_line,",  ",", ");
+         global_substitute(lo_line,")  ",") ");
+         global_substitute(lo_line,"*","* ");
+         global_substitute(lo_line,"* *","**");
+         global_substitute(lo_line,")+",") +");
+         global_substitute(lo_line,")-",") -");
+         global_substitute(lo_line,")/",") /");
+         global_substitute(lo_line,")*",") *");
+         global_substitute(lo_line,")=",") =");
+//       global_substitute(lo_line,"-","- ");
+         global_substitute(lo_line,"-  ","- ");
+         global_substitute(lo_line,"+","+ ");
+         global_substitute(lo_line,"+  ","+ ");
+         global_substitute(lo_line,"=- ","= -");
+         global_substitute(lo_line,"/JVS","/ JVS");
+         global_substitute(lo_line,"d- ","d-");
+         global_substitute(lo_line,"d+ ","d+");
+         global_substitute(lo_line,"D- ","D-");
+         global_substitute(lo_line,"D+ ","D+");
+         global_substitute(lo_line,"e- ","e-");
+         global_substitute(lo_line,"e+ ","e+");
+         global_substitute(lo_line,"E+ ","E+");
+         global_substitute(lo_line,"E- ","E-");
+         ip->global_subtolower(lo_line);
+         global_substitute(lo_line,"allocated","ALLOCATED");
+         global_substitute(lo_line,"allocatable","ALLOCATABLE");
+         global_substitute(lo_line,".and. ",".AND. ");
+         global_substitute(lo_line,"call ","CALL ");
+         global_substitute(lo_line,"case","CASE");
+         global_substitute(lo_line,"character","CHARACTER");
+         global_substitute(lo_line,"contains","CONTAINS");
+         global_substitute(lo_line,"deallocate","DEALLOCATE");
+         global_substitute(lo_line,"allocate","ALLOCATE");
+         global_substitute(lo_line,"dimension","DIMENSION");
+         global_substitute(lo_line,"do ","DO ");
+         global_substitute(lo_line,"elseif","ELSEIF");
+         global_substitute(lo_line,"else","ELSE");
+         global_substitute(lo_line,"#ELSE","#else");
+         global_substitute(lo_line,"end do","ENDDO");
+         global_substitute(lo_line,"end if","ENDIF");
+         global_substitute(lo_line,"endif","ENDIF");
+         global_substitute(lo_line,"endwhere","ENDWHERE");
+         global_substitute(lo_line,"end ","END ");    // Modify "end" after all other strings containing "end..." are done!
+         global_substitute(lo_line,"tEND","tend");
+         global_substitute(lo_line,"#ENDIF","#endif");
+         global_substitute(lo_line,"equivalence","EQUIVALENCE");
+         global_substitute(lo_line,"function","FUNCTION");
+         global_substitute(lo_line,"if(","IF (");
+         global_substitute(lo_line," if "," IF ");
+         global_substitute(lo_line,"implicit","IMPLICIT");
+         global_substitute(lo_line,"include","INCLUDE");
+         global_substitute(lo_line,"intent","INTENT");
+         global_substitute(lo_line,"integer","INTEGER");
+         global_substitute(lo_line,"interface","INTERFACE");
+         global_substitute(lo_line,"logical","LOGICAL");
+         global_substitute(lo_line,"module","MODULE");
+         global_substitute(lo_line,"none","NONE");
+         global_substitute(lo_line,"only","ONLY");
+         global_substitute(lo_line,"optional","OPTIONAL");
+         global_substitute(lo_line,"parameter","PARAMETER");
+         global_substitute(lo_line,"present","PRESENT");
+         global_substitute(lo_line,"private","PRIVATE");
+         global_substitute(lo_line,"procedure","PROCEDURE");
+         global_substitute(lo_line,"public","PUBLIC");
+         global_substitute(lo_line,"real","REAL");
+         global_substitute(lo_line,"return","RETURN");
+         global_substitute(lo_line,"use ","USE ");
+         global_substitute(lo_line,"save","SAVE");
+         global_substitute(lo_line,"subroutine","SUBROUTINE");
+         global_substitute(lo_line,"then","THEN");
+         global_substitute(lo_line,"where","WHERE");
+         global_substitute(lo_line,"while","WHILE");
+         global_substitute(lo_line,".false.",".FALSE.");
+         global_substitute(lo_line,".true.",".TRUE.");
+         global_substitute(lo_line,"(in)","(IN)");
+         global_substitute(lo_line,"(out)","(OUT)");
+         global_substitute(lo_line,"(inout)","(INOUT)");
+         global_substitute(lo_line,"\t","      ");
+//       cout << "HIER1 " << lo_line <<endl;
+       }
+     }
+//   pl.set_line(lo_line);
+//   pline.push_back(pl);
+
+
 //  line break if more than 130 character
 
-    lo_line = ip->get_line();
+//  lo_line = ip->get_line();
     if( lo_line.size() < 130 ) {
-      out << ip->get_line() <<endl;
+//    out << ip->get_line() <<endl;
+      out <<  lo_line    <<endl;
     } else {
       int cp  = lo_line.rfind("!",129);
Index: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file.h
===================================================================
--- palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file.h	(revision 2768)
+++ palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file.h	(revision 3298)
@@ -12,6 +12,16 @@
 // ############################################################################
 //
-//  Modifications:
-//    RFo added global_subtolower
+//Current revisions:
+//------------------
+/ 
+/ 
+//Former revisions:
+//-----------------
+// Nov 2016: Intial version of KP4 adapted to PALM (Klaus Ketelsen)
+//  forkel Sept.2018:  added edit_Initialize
+//  forkel          :  added global_subtolower
+//
+/ Intial version of KP4
+
 
 #include <iostream>
@@ -42,6 +52,10 @@
   void edit_KppSolve () ;
   void edit_Fun () ;
+  void edit_WAXPY () ;
+  void edit_FunTemplate () ;
+  void edit_JacTemplate () ;
   void edit_FUNC () ;
   void edit_Update_RCONST (vector <Vvar> &var_list) ;
+  void edit_Initialize (vector <Vvar> &var_list) ;
   void edit_inc (fortran_file & header_var) ;
   void create_species_list(vector <string> &species_list);
Index: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file_vec.C
===================================================================
--- palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file_vec.C	(revision 2768)
+++ palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/fortran_file_vec.C	(revision 3298)
@@ -9,4 +9,21 @@
 //
 // ############################################################################
+//Current revisions:
+//------------------
+/ 
+/ 
+//Former revisions:
+//-----------------
+//$Id: fortran_file.C 2470 2017-09-14 13:56:42Z forkel $
+//  ketelsen 18.09.2018: Line 112: do k=is,ie bydo k=1,vl; line 156 ff: replaced index k by j
+//
+//  forkel Sept.2018:  added edit_Initialize
+//                     changed loop direction for update_rconst in edit_Update_RCONST
+/ 2017-09-14 13:56:42Z forkel $
+//
+//
+// Nov 2016: Intial version of KP4 adapted to PALM (Klaus Ketelsen)
+//
+
 
 #include <fstream>
@@ -77,8 +94,19 @@
     }
   }
+
+  for (ip=pline.begin(); ip != pline.end(); ip++) {
+      if(kpp_switches.is_vector() ) {
+    	  ip->global_substitute ("phot (","phot(j,");
+      }
+  }
+
   ip = pline.begin()+1 ;
   lo_line = ip->get_line() ;
   lo_line.erase();
-  lo_line = " INTEGER         :: j,k";
+  if(kpp_switches.is_vector() ) {
+    lo_line = " INTEGER         :: j,k";
+  } else {
+    lo_line = " INTEGER         :: k";
+  }
   ip->set_line(lo_line);
   
@@ -87,5 +115,5 @@
     lo_line = ip->get_line() ;
     lo_line.erase();
-    lo_line = " do k=is,ie";
+    lo_line = " do k=1,vl";
     ip->set_line(lo_line);
 
@@ -93,5 +121,5 @@
     lo_line = ip->get_line() ;
     lo_line.erase();
-    lo_line = "  j = k-is+1";
+    lo_line = "  j = k";
     ip->set_line(lo_line);
 
@@ -103,4 +131,68 @@
   } else {
     ip = pline.begin()+3 ;
+    lo_line = ip->get_line() ;
+    lo_line.erase();
+    lo_line = "  k = is";
+    ip->set_line(lo_line);
+  }
+
+  return;
+}
+
+void fortran_file::edit_Initialize (vector <Vvar> &var_list) {
+  vector<program_line>::iterator     ip;
+  vector<Vvar>::iterator             iv;
+  string                             lo_line;
+
+  for (ip=pline.begin(); ip != pline.end(); ip++) {
+    ip->global_substitute ("*"," *");
+    ip->global_substitute ("* *","**");
+    ip->global_substitute (","," , ");
+    ip->global_substitute ("/"," / ");
+    ip->global_substitute ("1:VL","j");
+  }
+  for (ip=pline.begin(); ip != pline.end(); ip++) {
+    for (iv=var_list.begin(); iv != var_list.end(); iv++) {
+      ip->change_variable_to_vector_g (*iv);
+    }
+  }
+
+  for (ip=pline.begin(); ip != pline.end(); ip++) {
+    if(kpp_switches.is_vector() ) {
+          ip->global_substitute ("qvap","qvap(j)");
+    }
+  }
+  for (ip=pline.begin(); ip != pline.end(); ip++) {
+    if(kpp_switches.is_vector() ) {
+          ip->global_substitute ("fakt","fakt(j)");
+    }
+  }
+
+  ip = pline.begin()+3 ;
+  lo_line = ip->get_line() ;
+  lo_line.erase();
+  lo_line = "  INTEGER         :: j,k";
+  ip->set_line(lo_line);
+
+  if(kpp_switches.is_vector() ) {
+    ip = pline.begin()+9 ;
+    lo_line = ip->get_line() ;
+    lo_line.erase();
+    lo_line = " do k = is,ie";
+    ip->set_line(lo_line);
+
+    ip = pline.begin()+10 ;
+    lo_line = ip->get_line() ;
+    lo_line.erase();
+    lo_line = "  j = k - is +1";
+    ip->set_line(lo_line);
+
+    ip = pline.end()-2 ;
+    lo_line = ip->get_line() ;
+//  lo_line.erase();
+    lo_line = " end do";
+    ip->set_line(lo_line);
+  } else {
+    ip = pline.begin()+7 ;
     lo_line = ip->get_line() ;
     lo_line.erase();
@@ -236,2 +328,73 @@
 }
 
+void fortran_file::edit_WAXPY () {
+
+  vector<program_line>::iterator     ip;
+
+  cout << "Handling subroutine: WAXPY" <<endl;
+
+  for (ip=pline.begin(); ip != pline.end(); ip++) {
+
+    if(ip->get_token(0) == "REAL") {
+       ip->substitute("N",":,:");
+       ip->substitute("N",":,:");
+       ip->substitute("Alpha","Alpha(:)");
+    } else {
+        ip->change_variable_to_vector ("Y");
+        ip->change_variable_to_vector ("X");
+        if(ip->get_token(0) != "SUBROUTINE") {
+        	if(ip->get_token(0) == "IF")  {
+        		ip->substitute("Alpha","SUM(alpha(1:VL))");
+        	} else {
+        		ip->substitute("Alpha","alpha(1:VL)");
+        	}
+        }
+    }
+
+
+  }
+
+  return;
+}
+
+void fortran_file::edit_FunTemplate () {
+
+  vector<program_line>::iterator     ip;
+
+  cout << "Handling subroutine: FunTemplate" <<endl;
+
+  for (ip=pline.begin(); ip != pline.end(); ip++) {
+
+    if(ip->get_token(0) == "REAL") {
+       ip->substitute("NVAR",":,:");
+       ip->substitute("T,","T(:),");
+       ip->substitute("Told","Told(size(T)),Time(size(T))");
+    }
+
+
+  }
+
+  return;
+}
+
+void fortran_file::edit_JacTemplate () {
+
+  vector<program_line>::iterator     ip;
+
+  cout << "Handling subroutine: JacTemplate" <<endl;
+
+  for (ip=pline.begin(); ip != pline.end(); ip++) {
+
+    if(ip->get_token(0) == "REAL") {
+        ip->substitute("NVAR",":,:");
+        ip->substitute("T,","T(:),");
+        ip->substitute("Told","Told(size(T)),Time(size(T))");
+        ip->substitute("LU_NONZERO",":,:");
+    }
+
+
+  }
+
+  return;
+}
+
Index: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/program_line.C
===================================================================
--- palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/program_line.C	(revision 2768)
+++ palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/program_line.C	(revision 3298)
@@ -134,3 +134,17 @@
   return;
 }
+void   program_line::global_subtolower(string &line) {
 
+   int start = line.size()-1;
+   char c;
+
+    int i = 0;
+    while (line[i])
+    {
+      c = line[i];
+      line[i] = tolower(c);
+      i++;
+    }
+   return;
+}
+
Index: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/program_line.h
===================================================================
--- palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/program_line.h	(revision 2768)
+++ palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/src/program_line.h	(revision 3298)
@@ -11,4 +11,6 @@
 //
 // ############################################################################
+//
+// forkel June 2018: added    void   global_subtolower(string &line);
 
 #include <iostream>
@@ -42,4 +44,5 @@
    void   change_variable_to_vector    (string var);
    void   change_variable_to_vector_g  (Vvar  &var);
+   void   global_subtolower(string &line);
 };
 
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 2768)
+++ 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 2768)
+++ 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 2768)
+++ 	(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 2768)
+++ 	(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
