! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! Numerical Integrator (Time-Stepping) File
!
! Generated by KPP-2.1 symbolic chemistry Kinetics PreProcessor
!       (http://www.cs.vt.edu/~asandu/Software/KPP)
! KPP is distributed under GPL, the general public licence
!       (http://www.gnu.org/copyleft/gpl.html)
! (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa
! (C) 1997-2005, A. Sandu, Michigan Tech, Virginia Tech
!     With important contributions from:
!        M. Damian, Villanova University, USA
!        R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany
!
! File                 : radm2sorg_Integrator.f90
! Time                 : Mon Feb 27 11:08:59 2017
! Working directory    : /data/ksetigui/setigui/WRFV3_171015/chem/KPP/mechanisms/radm2sorg
! Equation file        : radm2sorg.kpp
! Output root filename : radm2sorg
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! INTEGRATE - Integrator routine
!   Arguments :
!      TIN       - Start Time for Integration
!      TOUT      - End Time for Integration
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

MODULE radm2sorg_Integrator

 USE radm2sorg_Parameters
 USE radm2sorg_Precision
 USE radm2sorg_JacobianSP

  IMPLICIT NONE

  INTEGER, PARAMETER :: ifun=1, ijac=2, istp=3, iacc=4, &
    irej=5, idec=6, isol=7, isng=8, itexit=1, ihexit=2

  ! description of the error numbers IERR
  CHARACTER(LEN=50), PARAMETER, DIMENSION(-8:1) :: IERR_NAMES = (/ &
    'Matrix is repeatedly singular                     ', & ! -8
    'Step size too small                               ', & ! -7
    'No of steps exceeds maximum bound                 ', & ! -6
    'Improper tolerance values                         ', & ! -5
    'FacMin/FacMax/FacRej must be positive             ', & ! -4
    'Hmin/Hmax/Hstart must be positive                 ', & ! -3
    'Selected Rosenbrock method not implemented        ', & ! -2
    'Improper value for maximal no of steps            ', & ! -1
    '                                                  ', & !  0 (not used)
    'Success                                           ' /) !  1

CONTAINS

SUBROUTINE  radm2sorg_INTEGRATE( TIN, TOUT, &
  FIX, VAR,  RCONST, ATOL, RTOL,           &
  ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U  )

   USE radm2sorg_Parameters
!!   USE radm2sorg_Global
   IMPLICIT NONE
   REAL(kind=dp), INTENT(INOUT), DIMENSION(NFIX) :: FIX
   REAL(kind=dp), INTENT(INOUT), DIMENSION(NVAR) :: VAR
   REAL(kind=dp), INTENT(IN), DIMENSION(NSPEC) :: ATOL, RTOL
   REAL(kind=dp), INTENT(IN), DIMENSION(NREACT) :: RCONST
   REAL(kind=dp), INTENT(IN) :: TIN  ! Start Time
   REAL(kind=dp), INTENT(IN) :: TOUT ! End Time
   ! Optional input parameters and statistics
   INTEGER,  INTENT(IN),  OPTIONAL :: ICNTRL_U(20)
   REAL(kind=dp), INTENT(IN),  OPTIONAL :: RCNTRL_U(20)
   INTEGER,  INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
   REAL(kind=dp), INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
   INTEGER,  INTENT(OUT), OPTIONAL :: IERR_U

   REAL(kind=dp)  :: STEPMIN

   INTEGER :: N_stp, N_acc, N_rej, N_sng
   SAVE N_stp, N_acc, N_rej, N_sng
   INTEGER :: i, IERR
   REAL(kind=dp) :: RCNTRL(20), RSTATUS(20)
   INTEGER :: ICNTRL(20), ISTATUS(20)

   ICNTRL(:)  = 0
   RCNTRL(:)  = 0.0_dp
   ISTATUS(:) = 0
   RSTATUS(:) = 0.0_dp

   ! If optional parameters are given, and if they are >0,
   ! then they overwrite default settings.
   IF (PRESENT(ICNTRL_U)) THEN
     WHERE(ICNTRL_U(:) > 0) ICNTRL(:) = ICNTRL_U(:)
   END IF
   IF (PRESENT(RCNTRL_U)) THEN
     WHERE(RCNTRL_U(:) > 0) RCNTRL(:) = RCNTRL_U(:)
   END IF

   CALL radm2sorg_Rosenbrock(VAR, FIX, RCONST, TIN,TOUT,   &
         ATOL,RTOL,               &
         RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)

   STEPMIN = RCNTRL(ihexit)
   ! if optional parameters are given for output they to return information
   IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(:)
   IF (PRESENT(RSTATUS_U)) RSTATUS_U(:) = RSTATUS(:)
   IF (PRESENT(IERR_U))    IERR_U       = IERR

END SUBROUTINE  radm2sorg_INTEGRATE

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE  radm2sorg_Rosenbrock(Y, FIX, RCONST, 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
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE radm2sorg_Parameters
!!  USE radm2sorg_LinearAlgebra
  IMPLICIT NONE

!~~~>  Arguments
   REAL(kind=dp), INTENT(INOUT) :: Y(NVAR)
   REAL(kind=dp), INTENT(IN), DIMENSION(NFIX) :: FIX
   REAL(kind=dp), INTENT(IN), DIMENSION(NREACT) :: RCONST
   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

!~~~>  Statistics on the work performed by the Rosenbrock method
  INTEGER :: Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng

!~~~>  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 radm2sorg_ros_ErrorMsg(-2,Tstart,ZERO,IERR)
      RETURN
   END IF

!~~~>   The maximum number of steps admitted
   IF (ICNTRL(4) == 0) THEN
      Max_no_steps = 100000
   ELSEIF (ICNTRL(4) > 0) THEN
      Max_no_steps=ICNTRL(4)
   ELSE
      PRINT * ,'User-selected max no. of steps: ICNTRL(4)=',ICNTRL(4)
      CALL radm2sorg_ros_ErrorMsg(-1,Tstart,ZERO,IERR)
      RETURN
   END IF

!~~~>  Unit roundoff (1+Roundoff>1)
   Roundoff = radm2sorg_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 radm2sorg_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 radm2sorg_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 radm2sorg_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 radm2sorg_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 radm2sorg_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 radm2sorg_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 radm2sorg_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 radm2sorg_ros_ErrorMsg(-5,Tstart,ZERO,IERR)
        RETURN
      END IF
    END DO

!~~~>   Initialize the particular Rosenbrock method
   SELECT CASE (Method)
     CASE (1)
       CALL radm2sorg_Ros2(ros_S, ros_A, ros_C, ros_M, ros_E,   &
          ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
     CASE (2)
       CALL radm2sorg_Ros3(ros_S, ros_A, ros_C, ros_M, ros_E,   &
          ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
     CASE (3)
       CALL radm2sorg_Ros4(ros_S, ros_A, ros_C, ros_M, ros_E,   &
          ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
     CASE (4)
       CALL radm2sorg_Rodas3(ros_S, ros_A, ros_C, ros_M, ros_E, &
          ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
     CASE (5)
       CALL radm2sorg_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(4)=', Method
       CALL radm2sorg_ros_ErrorMsg(-2,Tstart,ZERO,IERR)
       RETURN
   END SELECT

!~~~>  CALL Rosenbrock method
   CALL radm2sorg_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,                                    &
!  Statistics on the work performed by the Rosenbrock method
         Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng,&
!~~~>
         RCONST, FIX &
)

!~~~>  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  radm2sorg_ros_ErrorMsg(Code,T,H,IERR)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Handles all error messages
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   USE radm2sorg_Precision

   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  radm2sorg_ros_ErrorMsg

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 SUBROUTINE  radm2sorg_ros_Integrator (Y, Tstart, Tend, T,     &
        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,                                    &
!~~~>   Statistics on the work performed by the Rosenbrock method
        Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng, &
!~~~>
        RCONST, FIX &
 )
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!   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}
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  IMPLICIT NONE

!~~~> Input: the initial condition at Tstart; Output: the solution at T
   REAL(kind=dp), INTENT(INOUT) :: Y(NVAR)
!~~~> 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) ::  T
!~~~> 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
!~~~> Input
   REAL(kind=dp), INTENT(IN), DIMENSION(NFIX) :: FIX
!~~~> Input
   REAL(kind=dp), INTENT(IN), DIMENSION(NREACT) :: RCONST

!~~~>  Statistics on the work performed by the Rosenbrock method
  INTEGER, INTENT(INOUT)  :: Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng

! ~~~~ Local variables
   REAL(kind=dp) :: Ynew(NVAR), Fcn0(NVAR), Fcn(NVAR)
   REAL(kind=dp) :: K(NVAR*ros_S), dFdT(NVAR)
#ifdef FULL_ALGEBRA
   REAL(kind=dp) :: Jac0(NVAR,NVAR), Ghimj(NVAR,NVAR)
#else
   REAL(kind=dp) :: Jac0(LU_NONZERO), Ghimj(LU_NONZERO)
#endif
   REAL(kind=dp) :: H, Hnew, HC, HG, Fac, Tau
   REAL(kind=dp) :: Err, Yerr(NVAR)
   INTEGER :: Pivot(NVAR), Direction, ioffset, j, istage
   LOGICAL :: RejectLastH, RejectMoreH, Singular
!~~~>  Local parameters
   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
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

!~~~>  Initial preparations
   T = Tstart
   Hexit = 0.0_dp
   H = MIN(Hstart,Hmax)
   IF (ABS(H) <= 10.0_dp*Roundoff) H = DeltaMin

   IF (Tend  >=  Tstart) THEN
     Direction = +1
   ELSE
     Direction = -1
   END IF

   RejectLastH=.FALSE.
   RejectMoreH=.FALSE.

!~~~> Time loop begins below

TimeLoop: DO WHILE ( (Direction > 0).AND.((T-Tend)+Roundoff <= ZERO) &
       .OR. (Direction < 0).AND.((Tend-T)+Roundoff <= ZERO) )

   IF ( Nstp > Max_no_steps ) THEN  ! Too many steps
      CALL radm2sorg_ros_ErrorMsg(-6,T,H,IERR)
      RETURN
   END IF
   IF ( ((T+0.1_dp*H) == T).OR.(H <= Roundoff) ) THEN  ! Step size too small
      CALL radm2sorg_ros_ErrorMsg(-7,T,H,IERR)
      RETURN
   END IF

!~~~>  Limit H if necessary to avoid going beyond Tend
   Hexit = H
   H = MIN(H,ABS(Tend-T))

!~~~>   Compute the function at current time
   CALL radm2sorg_FunTemplate(T,Y,Fcn0, RCONST, FIX, Nfun)

!~~~>  Compute the function derivative with respect to T
   IF (.NOT.Autonomous) THEN
      CALL radm2sorg_ros_FunTimeDeriv ( T, Roundoff, Y, &
                Fcn0, dFdT, RCONST, FIX, Nfun )
   END IF

!~~~>   Compute the Jacobian at current time
   CALL radm2sorg_JacTemplate(T,Y,Jac0, FIX, Njac, RCONST)

!~~~>  Repeat step calculation until current step accepted
UntilAccepted: DO

   CALL radm2sorg_ros_PrepareMatrix(H,Direction,ros_Gamma(1), &
          Jac0,Ghimj,Pivot,Singular, Ndec,  Nsng )
   IF (Singular) THEN ! More than 5 consecutive failed decompositions
       CALL radm2sorg_ros_ErrorMsg(-8,T,H,IERR)
       RETURN
   END IF

!~~~>   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
         CALL radm2sorg_WCOPY(NVAR,Fcn0,1,Fcn,1)
      ! istage>1 and a new function evaluation is needed at the current istage
       ELSEIF ( ros_NewF(istage) ) THEN
         CALL radm2sorg_WCOPY(NVAR,Y,1,Ynew,1)
         DO j = 1, istage-1
           CALL radm2sorg_WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), &
            K(NVAR*(j-1)+1),1,Ynew,1)
         END DO
         Tau = T + ros_Alpha(istage)*Direction*H
         CALL radm2sorg_FunTemplate(Tau,Ynew,Fcn, RCONST, FIX, Nfun)
       END IF ! if istage == 1 elseif ros_NewF(istage)
       CALL radm2sorg_WCOPY(NVAR,Fcn,1,K(ioffset+1),1)
       DO j = 1, istage-1
         HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H)
         CALL radm2sorg_WAXPY(NVAR,HC,K(NVAR*(j-1)+1),1,K(ioffset+1),1)
       END DO
       IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
         HG = Direction*H*ros_Gamma(istage)
         CALL radm2sorg_WAXPY(NVAR,HG,dFdT,1,K(ioffset+1),1)
       END IF
       CALL radm2sorg_ros_Solve(Ghimj, Pivot, K(ioffset+1), Nsol)

   END DO Stage

!~~~>  Compute the new solution
   CALL radm2sorg_WCOPY(NVAR,Y,1,Ynew,1)
   DO j=1,ros_S
         CALL radm2sorg_WAXPY(NVAR,ros_M(j),K(NVAR*(j-1)+1),1,Ynew,1)
   END DO

!~~~>  Compute the error estimation
   CALL radm2sorg_WSCAL(NVAR,ZERO,Yerr,1)
   DO j=1,ros_S
        CALL radm2sorg_WAXPY(NVAR,ros_E(j),K(NVAR*(j-1)+1),1,Yerr,1)
   END DO
   Err = radm2sorg_ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol )

!~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
   Fac  = MIN(FacMax,MAX(FacMin,FacSafe/Err**(ONE/ros_ELO)))
   Hnew = H*Fac

!~~~>  Check the error magnitude and adjust step size
   Nstp = Nstp+1
   IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN  !~~~> Accept step
      Nacc = Nacc+1
      CALL radm2sorg_WCOPY(NVAR,Ynew,1,Y,1)
      T = T + Direction*H
      Hnew = MAX(Hmin,MIN(Hnew,Hmax))
      IF (RejectLastH) THEN  ! No step size increase after a rejected step
         Hnew = MIN(Hnew,H)
      END IF
      RejectLastH = .FALSE.
      RejectMoreH = .FALSE.
      H = Hnew
      EXIT UntilAccepted ! EXIT THE LOOP: WHILE STEP NOT ACCEPTED
   ELSE           !~~~> Reject step
      IF (RejectMoreH) THEN
         Hnew = H*FacRej
      END IF
      RejectMoreH = RejectLastH
      RejectLastH = .TRUE.
      H = Hnew
      IF (Nacc >= 1) THEN
         Nrej = Nrej+1
      END IF
   END IF ! Err <= 1

   END DO UntilAccepted

   END DO TimeLoop

!~~~> Succesful exit
   IERR = 1  !~~~> The integration was successful

  END SUBROUTINE  radm2sorg_ros_Integrator

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  REAL(kind=dp) FUNCTION  radm2sorg_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(NVAR), Ynew(NVAR), &
          Yerr(NVAR), AbsTol(NVAR), RelTol(NVAR)
   LOGICAL, INTENT(IN) ::  VectorTol
! Local variables
   REAL(kind=dp) :: Err, Scale, Ymax
   INTEGER  :: i
   REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp

   Err = ZERO
   DO i=1,NVAR
     Ymax = MAX(ABS(Y(i)),ABS(Ynew(i)))
     IF (VectorTol) THEN
       Scale = AbsTol(i)+RelTol(i)*Ymax
     ELSE
       Scale = AbsTol(1)+RelTol(1)*Ymax
     END IF
     Err = Err+(Yerr(i)/Scale)**2
   END DO
   Err  = SQRT(Err/NVAR)

    radm2sorg_ros_ErrorNorm = Err

  END FUNCTION  radm2sorg_ros_ErrorNorm

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  SUBROUTINE radm2sorg_ros_FunTimeDeriv ( T, Roundoff, Y, &
                Fcn0, dFdT, RCONST, FIX, Nfun )
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~> The time partial derivative of the function by finite differences
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   IMPLICIT NONE

!~~~> Input arguments
   REAL(kind=dp), INTENT(IN) :: T, Roundoff, Y(NVAR), Fcn0(NVAR)
   REAL(kind=dp), INTENT(IN) :: RCONST(NREACT), FIX(NFIX)
!~~~> Output arguments
   REAL(kind=dp), INTENT(OUT) :: dFdT(NVAR)
!~~~> InOut args
   INTEGER, INTENT(INOUT) ::Nfun
!~~~> Local variables
   REAL(kind=dp) :: Delta
   REAL(kind=dp), PARAMETER :: ONE = 1.0_dp, DeltaMin = 1.0E-6_dp

   Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T))
   CALL radm2sorg_FunTemplate(T+Delta,Y,dFdT, RCONST, FIX, Nfun)
   CALL radm2sorg_WAXPY(NVAR,(-ONE),Fcn0,1,dFdT,1)
   CALL radm2sorg_WSCAL(NVAR,(ONE/Delta),dFdT,1)

  END SUBROUTINE  radm2sorg_ros_FunTimeDeriv

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  SUBROUTINE  radm2sorg_ros_PrepareMatrix ( H, Direction, gam, &
             Jac0, Ghimj, Pivot, Singular, Ndec,  Nsng  )
! --- --- --- --- --- --- --- --- --- --- --- --- ---
!  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
#ifdef FULL_ALGEBRA
   REAL(kind=dp), INTENT(IN) ::  Jac0(NVAR,NVAR)
#else
   REAL(kind=dp), INTENT(IN) ::  Jac0(LU_NONZERO)
#endif
   REAL(kind=dp), INTENT(IN) ::  gam
   INTEGER, INTENT(IN) ::  Direction
!~~~> Output arguments
#ifdef FULL_ALGEBRA
   REAL(kind=dp), INTENT(OUT) :: Ghimj(NVAR,NVAR)
#else
   REAL(kind=dp), INTENT(OUT) :: Ghimj(LU_NONZERO)
#endif
   LOGICAL, INTENT(OUT) ::  Singular
   INTEGER, INTENT(OUT) ::  Pivot(NVAR)
!~~~> Inout arguments
   REAL(kind=dp), INTENT(INOUT) :: H   ! step size is decreased when LU fails
   INTEGER, INTENT(INOUT) ::  Ndec, Nsng
!~~~> Local variables
   INTEGER  :: i, ising, Nconsecutive
   REAL(kind=dp) :: ghinv
   REAL(kind=dp), PARAMETER :: ONE  = 1.0_dp, HALF = 0.5_dp

   Nconsecutive = 0
   Singular = .TRUE.

   DO WHILE (Singular)

!~~~>    Construct Ghimj = 1/(H*gam) - Jac0
#ifdef FULL_ALGEBRA
     CALL radm2sorg_WCOPY(NVAR*NVAR,Jac0,1,Ghimj,1)
     CALL radm2sorg_WSCAL(NVAR*NVAR,(-ONE),Ghimj,1)
     ghinv = ONE/(Direction*H*gam)
     DO i=1,NVAR
       Ghimj(i,i) = Ghimj(i,i)+ghinv
     END DO
#else
     CALL radm2sorg_WCOPY(LU_NONZERO,Jac0,1,Ghimj,1)
     CALL radm2sorg_WSCAL(LU_NONZERO,(-ONE),Ghimj,1)
     ghinv = ONE/(Direction*H*gam)
     DO i=1,NVAR
       Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv
     END DO
#endif
!~~~>    Compute LU decomposition
     CALL radm2sorg_ros_Decomp( Ghimj, Pivot, ising, Ndec )
     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
        Nsng = 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

  END SUBROUTINE  radm2sorg_ros_PrepareMatrix

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  SUBROUTINE  radm2sorg_ros_Decomp( A, Pivot, ising, Ndec )
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  Template for the LU decomposition
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   IMPLICIT NONE
!~~~> Inout variables
   REAL(kind=dp), INTENT(INOUT) :: A(LU_NONZERO)
!~~~> Output variables
   INTEGER, INTENT(OUT) :: Pivot(NVAR), ising
   INTEGER, INTENT(INOUT) :: Ndec

#ifdef FULL_ALGEBRA
   CALL  radm2sorg_DGETRF( NVAR, NVAR, A, NVAR, Pivot, ising )
#else
   CALL radm2sorg_KppDecomp ( A, ising )
   Pivot(1) = 1
#endif
   Ndec = Ndec + 1

  END SUBROUTINE  radm2sorg_ros_Decomp

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  SUBROUTINE  radm2sorg_ros_Solve( A, Pivot, b, Nsol )
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  Template for the forward/backward substitution (using pre-computed LU decomposition)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   IMPLICIT NONE
!~~~> Input variables
   REAL(kind=dp), INTENT(IN) :: A(LU_NONZERO)
   INTEGER, INTENT(IN) :: Pivot(NVAR)
!~~~~>  InOut args
   INTEGER, INTENT(INOUT) :: nsol
!~~~> InOut variables
   REAL(kind=dp), INTENT(INOUT) :: b(NVAR)

#ifdef FULL_ALGEBRA
   CALL  radm2sorg_DGETRS( 'N', NVAR , 1, A, NVAR, Pivot, b, NVAR, 0 )
#else
   CALL radm2sorg_KppSolve( A, b )
#endif

   Nsol = Nsol+1

  END SUBROUTINE  radm2sorg_ros_Solve

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  SUBROUTINE  radm2sorg_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
!cms   DOUBLE PRECISION g
    REAL(kind=dp) :: 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  radm2sorg_Ros2

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  SUBROUTINE  radm2sorg_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  radm2sorg_Ros3

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  SUBROUTINE  radm2sorg_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
!cms   DOUBLE PRECISION g
   REAL(kind=dp) :: 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  radm2sorg_Ros4

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  SUBROUTINE  radm2sorg_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
!cms   DOUBLE PRECISION g
   REAL(kind=dp) :: 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  radm2sorg_Rodas3

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  SUBROUTINE  radm2sorg_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
!cms   DOUBLE PRECISION g
    REAL(kind=dp) :: 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  radm2sorg_Rodas4

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!   End of the set of internal Rosenbrock subroutines
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
END SUBROUTINE  radm2sorg_Rosenbrock
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE  radm2sorg_FunTemplate( T, Y, Ydot, RCONST, FIX, Nfun )
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  Template for the ODE function call.
!  Updates the rate coefficients (and possibly the fixed species) at each call
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   USE radm2sorg_Parameters
!! USE radm2sorg_Global
!! USE radm2sorg_Function
!! USE radm2sorg_Rates
!~~~> Input variables
   REAL(kind=dp) :: T, Y(NVAR)
   REAL(kind=dp) :: RCONST(NREACT)
   REAL(kind=dp) :: FIX(NFIX)
!~~~> Output variables
   REAL(kind=dp) :: Ydot(NVAR)
   INTEGER :: Nfun

!~~~> Local variables
!!   REAL(kind=dp) :: Told

!!   Told = TIME
!!   TIME = T
!!   CALL Update_SUN()
!!   CALL Update_RCONST()
   CALL radm2sorg_Fun( Y, FIX, RCONST, Ydot )
!!   TIME = Told

   Nfun = Nfun+1

END SUBROUTINE  radm2sorg_FunTemplate

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE  radm2sorg_JacTemplate( T, Y, Jcb, FIX, Njac, RCONST )
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  Template for the ODE Jacobian call.
!  Updates the rate coefficients (and possibly the fixed species) at each call
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 USE radm2sorg_Parameters
 !!USE radm2sorg_Global
 USE radm2sorg_Jacobian
!! USE radm2sorg_LinearAlgebra
!! USE radm2sorg_Rates
!~~~> Input variables
    REAL(kind=dp) :: T, Y(NVAR)
    REAL(kind=dp) :: FIX(NFIX)
    REAL(kind=dp) :: RCONST(NREACT)

    INTEGER :: Njac

!~~~> Output variables
#ifdef FULL_ALGEBRA
    REAL(kind=dp) :: JV(LU_NONZERO), Jcb(NVAR,NVAR)
#else
    REAL(kind=dp) :: Jcb(LU_NONZERO)
#endif
!~~~> Local variables
    REAL(kind=dp) :: Told
#ifdef FULL_ALGEBRA
    INTEGER :: i, j
#endif

!!    Told = TIME
!!    TIME = T
!!    CALL Update_SUN()
!!    CALL Update_RCONST()
#ifdef FULL_ALGEBRA
    CALL radm2sorg_Jac_SP(Y, FIX, RCONST, JV)
    DO j=1,NVAR
      DO i=1,NVAR
         Jcb(i,j) = 0.0d0
      END DO
    END DO
    DO i=1,LU_NONZERO
       Jcb(LU_IROW(i),LU_ICOL(i)) = JV(i)
    END DO
#else
    CALL radm2sorg_Jac_SP( Y, FIX, RCONST, Jcb )
#endif
!!    TIME = Told

    Njac = Njac+1

END SUBROUTINE  radm2sorg_JacTemplate

!!!END MODULE radm2sorg_Integrator
! End of INTEGRATE function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! radm2sorg_Fun - time derivatives of variables - Agregate form
!   Arguments :
!      V         - Concentrations of variable species (local)
!      F         - Concentrations of fixed species (local)
!      RCT       - Rate constants (local)
!      Vdot      - Time derivative of variable species concentrations
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE radm2sorg_Fun ( V, F, RCT, Vdot )

! V - Concentrations of variable species (local)
  REAL(kind=dp) :: V(NVAR)
! F - Concentrations of fixed species (local)
  REAL(kind=dp) :: F(NFIX)
! RCT - Rate constants (local)
  REAL(kind=dp) :: RCT(NREACT)
! Vdot - Time derivative of variable species concentrations
  REAL(kind=dp) :: Vdot(NVAR)

! Local variables
! A - Rate for each equation
  REAL(kind=dp) :: A(NREACT)

! Computation of equation rates
  A(1) = RCT(1)*V(59)
  A(2) = RCT(2)*V(46)
  A(3) = RCT(3)*V(46)
  A(4) = RCT(4)*V(13)
  A(5) = RCT(5)*V(25)
  A(6) = RCT(6)*V(20)
  A(7) = RCT(7)*V(54)
  A(8) = RCT(8)*V(54)
  A(9) = RCT(9)*V(14)
  A(10) = RCT(10)*V(43)
  A(11) = RCT(11)*V(43)
  A(12) = RCT(12)*V(48)
  A(13) = RCT(13)*V(21)
  A(14) = RCT(14)*V(42)
  A(15) = RCT(15)*V(18)
  A(16) = RCT(16)*V(33)
  A(17) = RCT(17)*V(31)
  A(18) = RCT(18)*V(31)
  A(19) = RCT(19)*V(34)
  A(20) = RCT(20)*V(30)
  A(21) = RCT(21)*V(47)
  A(22) = RCT(22)*V(19)*F(2)
  A(23) = RCT(23)*V(19)*V(59)
  A(24) = RCT(24)*V(7)*F(2)
  A(25) = RCT(25)*V(7)*F(1)
  A(26) = RCT(26)*V(46)*V(58)
  A(27) = RCT(27)*V(46)*V(53)
  A(28) = RCT(28)*V(46)*V(51)
  A(29) = RCT(29)*V(51)*V(58)
  A(30) = RCT(30)*V(51)*V(59)
  A(31) = RCT(31)*V(20)
  A(32) = RCT(32)*V(51)*V(51)
  A(33) = RCT(33)*V(51)*V(51)*F(1)
  A(34) = RCT(34)*V(14)*V(53)
  A(35) = RCT(35)*V(53)*V(58)
  A(36) = RCT(36)*V(58)*V(58)*F(2)
  A(37) = RCT(37)*V(46)*V(59)
  A(38) = RCT(38)*V(54)*V(58)
  A(39) = RCT(39)*V(54)*V(59)
  A(40) = RCT(40)*V(51)*V(54)
  A(41) = RCT(41)*V(54)*V(59)
  A(42) = RCT(42)*V(15)
  A(43) = RCT(43)*V(15)
  A(44) = RCT(44)*V(53)*V(59)
  A(45) = RCT(45)*V(25)*V(53)
  A(46) = RCT(46)*V(20)*V(53)
  A(47) = RCT(47)*V(51)*V(53)
  A(48) = RCT(48)*V(5)*V(53)
  A(49) = RCT(49)*V(26)*V(53)
  A(50) = RCT(50)*V(17)*V(53)
  A(51) = RCT(51)*V(6)*V(53)
  A(52) = RCT(52)*V(16)*V(53)
  A(53) = RCT(53)*V(8)*V(53)
  A(54) = RCT(54)*V(9)*V(53)
  A(55) = RCT(55)*V(24)*V(53)
  A(56) = RCT(56)*V(28)*V(53)
  A(57) = RCT(57)*V(29)*V(53)
  A(58) = RCT(58)*V(10)*V(53)
  A(59) = RCT(59)*V(11)*V(53)
  A(60) = RCT(60)*V(22)*V(53)
  A(61) = RCT(61)*V(43)*V(53)
  A(62) = RCT(62)*V(48)*V(53)
  A(63) = RCT(63)*V(33)*V(53)
  A(64) = RCT(64)*V(31)*V(53)
  A(65) = RCT(65)*V(34)*V(53)
  A(66) = RCT(66)*V(30)*V(53)
  A(67) = RCT(67)*V(21)*V(53)
  A(68) = RCT(68)*V(42)*V(53)
  A(69) = RCT(69)*V(18)*V(53)
  A(70) = RCT(70)*V(23)*V(53)
  A(71) = RCT(71)*V(47)*V(53)
  A(72) = RCT(72)*V(27)*V(53)
  A(73) = RCT(73)*V(55)*V(59)
  A(74) = RCT(74)*V(23)
  A(75) = RCT(75)*V(45)*V(59)
  A(76) = RCT(76)*V(12)
  A(77) = RCT(77)*V(52)*V(58)
  A(78) = RCT(78)*V(56)*V(58)
  A(79) = RCT(79)*V(41)*V(58)
  A(80) = RCT(80)*V(44)*V(58)
  A(81) = RCT(81)*V(40)*V(58)
  A(82) = RCT(82)*V(37)*V(58)
  A(83) = RCT(83)*V(49)*V(58)
  A(84) = RCT(84)*V(55)*V(58)
  A(85) = RCT(85)*V(45)*V(58)
  A(86) = RCT(86)*V(35)*V(58)
  A(87) = RCT(87)*V(36)*V(58)
  A(88) = RCT(88)*V(57)*V(58)
  A(89) = RCT(89)*V(50)*V(58)
  A(90) = RCT(90)*V(38)*V(58)
  A(91) = RCT(91)*V(43)*V(54)
  A(92) = RCT(92)*V(48)*V(54)
  A(93) = RCT(93)*V(31)*V(54)
  A(94) = RCT(94)*V(34)*V(54)
  A(95) = RCT(95)*V(30)*V(54)
  A(96) = RCT(96)*V(22)*V(54)
  A(97) = RCT(97)*V(24)*V(54)
  A(98) = RCT(98)*V(28)*V(54)
  A(99) = RCT(99)*V(29)*V(54)
  A(100) = RCT(100)*V(27)*V(54)
  A(101) = RCT(101)*V(24)*V(46)
  A(102) = RCT(102)*V(28)*V(46)
  A(103) = RCT(103)*V(29)*V(46)
  A(104) = RCT(104)*V(27)*V(46)
  A(105) = RCT(105)*V(51)*V(52)
  A(106) = RCT(106)*V(51)*V(57)
  A(107) = RCT(107)*V(51)*V(56)
  A(108) = RCT(108)*V(41)*V(51)
  A(109) = RCT(109)*V(44)*V(51)
  A(110) = RCT(110)*V(40)*V(51)
  A(111) = RCT(111)*V(37)*V(51)
  A(112) = RCT(112)*V(49)*V(51)
  A(113) = RCT(113)*V(50)*V(51)
  A(114) = RCT(114)*V(51)*V(55)
  A(115) = RCT(115)*V(35)*V(51)
  A(116) = RCT(116)*V(36)*V(51)
  A(117) = RCT(117)*V(45)*V(51)
  A(118) = RCT(118)*V(38)*V(51)
  A(119) = RCT(119)*V(52)*V(52)
  A(120) = RCT(120)*V(52)*V(57)
  A(121) = RCT(121)*V(52)*V(56)
  A(122) = RCT(122)*V(41)*V(52)
  A(123) = RCT(123)*V(44)*V(52)
  A(124) = RCT(124)*V(40)*V(52)
  A(125) = RCT(125)*V(37)*V(52)
  A(126) = RCT(126)*V(49)*V(52)
  A(127) = RCT(127)*V(50)*V(52)
  A(128) = RCT(128)*V(52)*V(55)
  A(129) = RCT(129)*V(35)*V(52)
  A(130) = RCT(130)*V(36)*V(52)
  A(131) = RCT(131)*V(45)*V(52)
  A(132) = RCT(132)*V(55)*V(57)
  A(133) = RCT(133)*V(55)*V(56)
  A(134) = RCT(134)*V(41)*V(55)
  A(135) = RCT(135)*V(44)*V(55)
  A(136) = RCT(136)*V(40)*V(55)
  A(137) = RCT(137)*V(37)*V(55)
  A(138) = RCT(138)*V(49)*V(55)
  A(139) = RCT(139)*V(50)*V(55)
  A(140) = RCT(140)*V(55)*V(55)
  A(141) = RCT(141)*V(35)*V(55)
  A(142) = RCT(142)*V(36)*V(55)
  A(143) = RCT(143)*V(45)*V(55)
  A(144) = RCT(144)*V(39)*V(51)
  A(145) = RCT(145)*V(39)*V(52)
  A(146) = RCT(146)*V(39)*V(55)
  A(147) = RCT(147)*V(39)*V(39)
  A(148) = RCT(148)*V(39)*V(58)
  A(149) = RCT(149)*V(32)*V(59)
  A(150) = RCT(150)*V(32)*V(51)
  A(151) = RCT(151)*V(32)*V(52)
  A(152) = RCT(152)*V(32)*V(55)
  A(153) = RCT(153)*V(32)*V(32)
  A(154) = RCT(154)*V(38)*V(52)
  A(155) = RCT(155)*V(38)*V(55)
  A(156) = RCT(156)*V(38)*V(38)

! Aggregate function
  Vdot(1) = A(48)
  Vdot(2) = 0.4*A(101)+0.2*A(102)+0.06*A(103)+0.2*A(104)
  Vdot(3) = 0.2*A(102)+0.29*A(103)+0.2*A(104)+0.5*A(128)+0.5*A(131)+0.5*A(132)+0.5*A(133)+0.5*A(134)+0.5*A(135)+0.5&
              &*A(136)+0.5*A(137)+0.5*A(138)+0.5*A(139)+0.5*A(155)
  Vdot(4) = A(49)
  Vdot(5) = -A(48)
  Vdot(6) = -A(51)
  Vdot(7) = A(2)-A(24)-A(25)
  Vdot(8) = -A(53)
  Vdot(9) = -A(54)
  Vdot(10) = -A(58)
  Vdot(11) = -A(59)
  Vdot(12) = A(75)-A(76)
  Vdot(13) = -A(4)+A(35)
  Vdot(14) = -A(9)+A(32)+A(33)-A(34)
  Vdot(15) = A(41)-A(42)-A(43)
  Vdot(16) = -A(52)
  Vdot(17) = -A(50)+0.06*A(102)+0.09*A(103)
  Vdot(18) = -A(15)-A(69)+A(114)
  Vdot(19) = A(1)+A(3)+A(8)-A(22)-A(23)+A(24)
  Vdot(20) = -A(6)+A(30)-A(31)-A(46)
  Vdot(21) = -A(13)-A(67)+A(105)
  Vdot(22) = 0.25*A(58)+0.17*A(59)-A(60)-0.5*A(96)
  Vdot(23) = -A(70)+A(73)-A(74)
  Vdot(24) = -A(55)-A(97)-A(101)
  Vdot(25) = -A(5)+A(40)+2*A(43)+A(44)-A(45)+A(91)+A(92)+A(93)+A(94)+A(95)+A(96)
  Vdot(26) = A(10)+A(11)+A(12)+1.87*A(17)+1.55*A(18)+A(19)-A(49)+A(61)+2*A(64)+A(65)+0.95*A(85)+A(91)+2*A(93)+A(94)+0.42&
               &*A(101)+0.33*A(102)+0.23*A(103)+0.33*A(104)+0.475*A(131)+0.95*A(143)
  Vdot(27) = -A(72)-A(100)-A(104)
  Vdot(28) = -A(56)-A(98)-A(102)
  Vdot(29) = -A(57)-A(99)-A(103)
  Vdot(30) = -A(20)-A(66)+0.7*A(86)+0.806*A(87)-A(95)+0.7*A(129)+0.806*A(130)+A(141)+A(142)
  Vdot(31) = -A(17)-A(18)-A(64)+0.89*A(85)+0.16*A(86)-A(93)+0.16*A(129)+0.445*A(131)+0.2*A(141)+0.89*A(143)
  Vdot(32) = A(96)-A(149)-A(150)-A(151)-A(152)-2*A(153)
  Vdot(33) = -A(16)+0.8*A(21)+0.025*A(52)-A(63)+0.25*A(78)+0.69*A(79)+1.06*A(80)+0.1*A(83)+0.1*A(103)+0.6*A(121)+0.75&
               &*A(122)+1.39*A(123)+0.55*A(126)+0.8*A(133)+0.86*A(134)+0.9*A(135)+0.55*A(138)
  Vdot(34) = -A(19)-A(65)+0.11*A(85)+0.17*A(86)+0.45*A(87)+A(89)-A(94)+0.75*A(127)+0.17*A(129)+0.45*A(130)+0.055*A(131)&
               &+A(139)+0.8*A(141)+A(142)+0.11*A(143)
  Vdot(35) = 0.75*A(58)-A(86)-A(115)-A(129)-A(141)
  Vdot(36) = 0.83*A(59)-A(87)-A(116)-A(130)-A(142)
  Vdot(37) = A(56)+A(72)-A(82)-A(111)-A(125)-A(137)
  Vdot(38) = -A(90)+A(97)+A(98)+A(99)+A(100)-A(118)-A(154)-A(155)-2*A(156)
  Vdot(39) = 0.25*A(53)+0.75*A(54)+0.9*A(60)+A(70)+2*A(85)+A(131)+2*A(143)-A(144)-A(145)-A(146)-2*A(147)-A(148)
  Vdot(40) = A(55)-A(81)-A(110)-A(124)-A(136)
  Vdot(41) = A(53)-A(79)-A(108)-A(122)-A(134)
  Vdot(42) = -A(14)-A(68)+A(106)+A(107)+A(108)+A(109)+A(110)+A(111)+A(112)+A(113)+A(115)+A(116)+A(117)+A(144)+A(150)
  Vdot(43) = -A(10)-A(11)+A(13)+0.13*A(17)+0.45*A(18)+0.009*A(52)-A(61)+0.5*A(67)+A(70)+A(77)+0.09*A(78)+0.04*A(80)+1.6&
               &*A(81)+A(82)+0.28*A(83)+A(90)-A(91)+A(101)+0.53*A(102)+0.18*A(103)+0.53*A(104)+1.5*A(119)+0.75*A(120)+0.75&
               &*A(121)+0.77*A(122)+0.8*A(123)+1.55*A(124)+1.25*A(125)+0.89*A(126)+0.75*A(127)+A(128)+A(129)+A(130)+0.5&
               &*A(131)+0.8*A(136)+0.5*A(137)+0.14*A(138)+A(145)+A(151)+1.75*A(154)+A(155)+2*A(156)
  Vdot(44) = A(54)-A(80)-A(109)-A(123)-A(135)
  Vdot(45) = A(20)+0.9*A(60)+A(66)-A(75)+A(76)-A(85)+A(95)-A(117)-A(131)-A(143)
  Vdot(46) = -A(2)-A(3)+A(22)-A(26)-A(27)-A(28)-A(37)-A(101)-A(102)-A(103)-A(104)
  Vdot(47) = -A(21)-A(71)+0.036*A(78)+0.08*A(79)+0.24*A(80)+A(118)+A(149)
  Vdot(48) = -A(12)+A(14)+0.2*A(21)+0.075*A(52)-A(62)+0.5*A(68)+0.75*A(78)+0.38*A(79)+0.35*A(80)+0.2*A(81)+A(82)+1.45&
               &*A(83)+A(88)+A(90)-A(92)+0.5*A(102)+0.72*A(103)+0.5*A(104)+0.75*A(120)+0.15*A(121)+0.41*A(122)+0.46*A(123)&
               &+0.35*A(124)+0.75*A(125)+0.725*A(126)+A(132)+0.2*A(133)+0.14*A(134)+0.1*A(135)+0.6*A(136)+A(137)+0.725&
               &*A(138)+A(154)+A(155)+2*A(156)
  Vdot(49) = A(57)-A(83)-A(112)-A(126)-A(138)
  Vdot(50) = A(63)-A(89)-A(113)-A(127)-A(139)
  Vdot(51) = 0.65*A(6)+2*A(11)+A(12)+A(13)+A(14)+0.8*A(18)+A(19)+A(20)+A(21)+A(27)-A(28)-A(29)-A(30)+A(31)-2*A(32)-2&
               &*A(33)+A(34)-A(40)-A(47)+A(48)+A(49)+0.17*A(52)+0.25*A(58)+0.17*A(59)+0.1*A(60)+A(61)+A(64)+A(77)+0.964&
               &*A(78)+0.92*A(79)+0.76*A(80)+A(81)+A(82)+A(83)+0.92*A(85)+A(86)+A(87)+A(88)+A(89)+A(91)+A(93)+0.12*A(101)&
               &+0.23*A(102)+0.26*A(103)+0.23*A(104)-A(105)-A(106)-A(107)-A(108)-A(109)-A(110)-A(111)-A(112)-A(113)-A(114)&
               &-A(115)-A(116)-A(117)-A(118)+A(119)+A(120)+A(121)+A(122)+A(123)+A(124)+A(125)+A(126)+A(127)+0.5*A(128)+2&
               &*A(129)+2*A(130)+0.46*A(131)+0.5*A(132)+0.5*A(133)+0.5*A(134)+0.5*A(135)+0.5*A(136)+0.5*A(137)+0.5*A(138)&
               &+0.5*A(139)+A(141)+A(142)+0.92*A(143)-A(144)+A(145)-A(150)+A(151)+0.5*A(154)
  Vdot(52) = A(12)+A(15)+A(50)+0.5*A(67)-A(77)+A(84)+0.22*A(102)+0.31*A(103)+0.22*A(104)-A(105)-2*A(119)-A(120)-A(121)&
               &-A(122)-A(123)-A(124)-A(125)-A(126)-A(127)-0.5*A(128)-A(129)-A(130)-A(131)+0.5*A(132)+0.5*A(133)+0.5*A(134)&
               &+0.5*A(135)+0.5*A(136)+0.5*A(137)+0.5*A(138)+0.5*A(139)+2*A(140)+A(141)+A(142)+A(143)-A(145)+A(146)-A(151)&
               &+A(152)-A(154)+0.5*A(155)
  Vdot(53) = A(4)+A(5)+0.35*A(6)+2*A(9)+A(13)+A(14)+A(15)+2*A(25)-A(27)+A(28)+A(29)-A(34)-A(35)-A(44)-A(45)-A(46)-A(47)&
               &-A(48)-A(49)-A(50)-A(51)-A(52)-A(53)-A(54)-A(55)-A(56)-A(57)-A(58)-A(59)-1.9*A(60)-A(61)-A(62)-A(63)-A(64)&
               &-A(65)-A(66)-0.5*A(67)-0.5*A(68)-A(69)-A(70)-A(71)-A(72)+0.1*A(102)+0.14*A(103)+0.1*A(104)
  Vdot(54) = 0.35*A(6)-A(7)-A(8)+A(37)-A(38)-A(39)-A(40)-A(41)+A(42)+A(45)+A(70)-A(91)-A(92)-A(93)-A(94)-A(95)-A(96)&
               &-A(97)-A(98)-A(99)-A(100)
  Vdot(55) = A(16)+A(19)+A(62)+A(65)+A(69)-A(73)+A(74)-A(84)+0.05*A(85)+A(92)+A(94)-A(114)-A(128)+0.025*A(131)-A(132)&
               &-A(133)-A(134)-A(135)-A(136)-A(137)-A(138)-A(139)-2*A(140)-A(141)-A(142)-0.95*A(143)-A(146)-A(152)-A(155)
  Vdot(56) = 0.83*A(52)+0.5*A(68)+A(71)-A(78)-A(107)-A(121)-A(133)
  Vdot(57) = A(16)+A(51)-A(88)-A(106)-A(120)-A(132)
  Vdot(58) = A(1)+A(4)+A(7)+A(23)-A(26)-A(29)-A(35)-2*A(36)-A(38)+A(39)-A(77)-A(78)-A(79)-A(80)-A(81)-A(82)-A(83)-A(84)&
               &-A(85)-A(86)-A(87)-A(88)-A(89)-A(90)-A(148)
  Vdot(59) = -A(1)+A(5)+0.65*A(6)+A(8)+A(21)-A(23)+A(26)+A(29)-A(30)+A(31)+2*A(36)-A(37)+2*A(38)-A(41)+A(42)-A(44)+A(46)&
               &+A(71)-A(73)+A(74)-A(75)+A(76)+A(77)+0.964*A(78)+0.92*A(79)+0.76*A(80)+A(81)+A(82)+A(83)+A(84)+A(85)+A(86)&
               &+A(87)+A(88)+A(89)+2*A(90)+A(148)-A(149)+A(154)+A(155)+2*A(156)

END SUBROUTINE radm2sorg_Fun

! End of radm2sorg_Fun function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! radm2sorg_Jac_SP - the Jacobian of Variables in sparse matrix representation
!   Arguments :
!      V         - Concentrations of variable species (local)
!      F         - Concentrations of fixed species (local)
!      RCT       - Rate constants (local)
!      JVS       - sparse Jacobian of variables
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE radm2sorg_Jac_SP ( V, F, RCT, JVS )

! V - Concentrations of variable species (local)
  REAL(kind=dp) :: V(NVAR)
! F - Concentrations of fixed species (local)
  REAL(kind=dp) :: F(NFIX)
! RCT - Rate constants (local)
  REAL(kind=dp) :: RCT(NREACT)
! JVS - sparse Jacobian of variables
  REAL(kind=dp) :: JVS(LU_NONZERO)

! Local variables
! B - Temporary array
  REAL(kind=dp) :: B(280)

! B(1) = dA(1)/dV(59)
  B(1) = RCT(1)
! B(2) = dA(2)/dV(46)
  B(2) = RCT(2)
! B(3) = dA(3)/dV(46)
  B(3) = RCT(3)
! B(4) = dA(4)/dV(13)
  B(4) = RCT(4)
! B(5) = dA(5)/dV(25)
  B(5) = RCT(5)
! B(6) = dA(6)/dV(20)
  B(6) = RCT(6)
! B(7) = dA(7)/dV(54)
  B(7) = RCT(7)
! B(8) = dA(8)/dV(54)
  B(8) = RCT(8)
! B(9) = dA(9)/dV(14)
  B(9) = RCT(9)
! B(10) = dA(10)/dV(43)
  B(10) = RCT(10)
! B(11) = dA(11)/dV(43)
  B(11) = RCT(11)
! B(12) = dA(12)/dV(48)
  B(12) = RCT(12)
! B(13) = dA(13)/dV(21)
  B(13) = RCT(13)
! B(14) = dA(14)/dV(42)
  B(14) = RCT(14)
! B(15) = dA(15)/dV(18)
  B(15) = RCT(15)
! B(16) = dA(16)/dV(33)
  B(16) = RCT(16)
! B(17) = dA(17)/dV(31)
  B(17) = RCT(17)
! B(18) = dA(18)/dV(31)
  B(18) = RCT(18)
! B(19) = dA(19)/dV(34)
  B(19) = RCT(19)
! B(20) = dA(20)/dV(30)
  B(20) = RCT(20)
! B(21) = dA(21)/dV(47)
  B(21) = RCT(21)
! B(22) = dA(22)/dV(19)
  B(22) = RCT(22)*F(2)
! B(24) = dA(23)/dV(19)
  B(24) = RCT(23)*V(59)
! B(25) = dA(23)/dV(59)
  B(25) = RCT(23)*V(19)
! B(26) = dA(24)/dV(7)
  B(26) = RCT(24)*F(2)
! B(28) = dA(25)/dV(7)
  B(28) = RCT(25)*F(1)
! B(30) = dA(26)/dV(46)
  B(30) = RCT(26)*V(58)
! B(31) = dA(26)/dV(58)
  B(31) = RCT(26)*V(46)
! B(32) = dA(27)/dV(46)
  B(32) = RCT(27)*V(53)
! B(33) = dA(27)/dV(53)
  B(33) = RCT(27)*V(46)
! B(34) = dA(28)/dV(46)
  B(34) = RCT(28)*V(51)
! B(35) = dA(28)/dV(51)
  B(35) = RCT(28)*V(46)
! B(36) = dA(29)/dV(51)
  B(36) = RCT(29)*V(58)
! B(37) = dA(29)/dV(58)
  B(37) = RCT(29)*V(51)
! B(38) = dA(30)/dV(51)
  B(38) = RCT(30)*V(59)
! B(39) = dA(30)/dV(59)
  B(39) = RCT(30)*V(51)
! B(40) = dA(31)/dV(20)
  B(40) = RCT(31)
! B(41) = dA(32)/dV(51)
  B(41) = RCT(32)*2*V(51)
! B(42) = dA(33)/dV(51)
  B(42) = RCT(33)*2*V(51)*F(1)
! B(44) = dA(34)/dV(14)
  B(44) = RCT(34)*V(53)
! B(45) = dA(34)/dV(53)
  B(45) = RCT(34)*V(14)
! B(46) = dA(35)/dV(53)
  B(46) = RCT(35)*V(58)
! B(47) = dA(35)/dV(58)
  B(47) = RCT(35)*V(53)
! B(48) = dA(36)/dV(58)
  B(48) = RCT(36)*2*V(58)*F(2)
! B(50) = dA(37)/dV(46)
  B(50) = RCT(37)*V(59)
! B(51) = dA(37)/dV(59)
  B(51) = RCT(37)*V(46)
! B(52) = dA(38)/dV(54)
  B(52) = RCT(38)*V(58)
! B(53) = dA(38)/dV(58)
  B(53) = RCT(38)*V(54)
! B(54) = dA(39)/dV(54)
  B(54) = RCT(39)*V(59)
! B(55) = dA(39)/dV(59)
  B(55) = RCT(39)*V(54)
! B(56) = dA(40)/dV(51)
  B(56) = RCT(40)*V(54)
! B(57) = dA(40)/dV(54)
  B(57) = RCT(40)*V(51)
! B(58) = dA(41)/dV(54)
  B(58) = RCT(41)*V(59)
! B(59) = dA(41)/dV(59)
  B(59) = RCT(41)*V(54)
! B(60) = dA(42)/dV(15)
  B(60) = RCT(42)
! B(61) = dA(43)/dV(15)
  B(61) = RCT(43)
! B(62) = dA(44)/dV(53)
  B(62) = RCT(44)*V(59)
! B(63) = dA(44)/dV(59)
  B(63) = RCT(44)*V(53)
! B(64) = dA(45)/dV(25)
  B(64) = RCT(45)*V(53)
! B(65) = dA(45)/dV(53)
  B(65) = RCT(45)*V(25)
! B(66) = dA(46)/dV(20)
  B(66) = RCT(46)*V(53)
! B(67) = dA(46)/dV(53)
  B(67) = RCT(46)*V(20)
! B(68) = dA(47)/dV(51)
  B(68) = RCT(47)*V(53)
! B(69) = dA(47)/dV(53)
  B(69) = RCT(47)*V(51)
! B(70) = dA(48)/dV(5)
  B(70) = RCT(48)*V(53)
! B(71) = dA(48)/dV(53)
  B(71) = RCT(48)*V(5)
! B(72) = dA(49)/dV(26)
  B(72) = RCT(49)*V(53)
! B(73) = dA(49)/dV(53)
  B(73) = RCT(49)*V(26)
! B(74) = dA(50)/dV(17)
  B(74) = RCT(50)*V(53)
! B(75) = dA(50)/dV(53)
  B(75) = RCT(50)*V(17)
! B(76) = dA(51)/dV(6)
  B(76) = RCT(51)*V(53)
! B(77) = dA(51)/dV(53)
  B(77) = RCT(51)*V(6)
! B(78) = dA(52)/dV(16)
  B(78) = RCT(52)*V(53)
! B(79) = dA(52)/dV(53)
  B(79) = RCT(52)*V(16)
! B(80) = dA(53)/dV(8)
  B(80) = RCT(53)*V(53)
! B(81) = dA(53)/dV(53)
  B(81) = RCT(53)*V(8)
! B(82) = dA(54)/dV(9)
  B(82) = RCT(54)*V(53)
! B(83) = dA(54)/dV(53)
  B(83) = RCT(54)*V(9)
! B(84) = dA(55)/dV(24)
  B(84) = RCT(55)*V(53)
! B(85) = dA(55)/dV(53)
  B(85) = RCT(55)*V(24)
! B(86) = dA(56)/dV(28)
  B(86) = RCT(56)*V(53)
! B(87) = dA(56)/dV(53)
  B(87) = RCT(56)*V(28)
! B(88) = dA(57)/dV(29)
  B(88) = RCT(57)*V(53)
! B(89) = dA(57)/dV(53)
  B(89) = RCT(57)*V(29)
! B(90) = dA(58)/dV(10)
  B(90) = RCT(58)*V(53)
! B(91) = dA(58)/dV(53)
  B(91) = RCT(58)*V(10)
! B(92) = dA(59)/dV(11)
  B(92) = RCT(59)*V(53)
! B(93) = dA(59)/dV(53)
  B(93) = RCT(59)*V(11)
! B(94) = dA(60)/dV(22)
  B(94) = RCT(60)*V(53)
! B(95) = dA(60)/dV(53)
  B(95) = RCT(60)*V(22)
! B(96) = dA(61)/dV(43)
  B(96) = RCT(61)*V(53)
! B(97) = dA(61)/dV(53)
  B(97) = RCT(61)*V(43)
! B(98) = dA(62)/dV(48)
  B(98) = RCT(62)*V(53)
! B(99) = dA(62)/dV(53)
  B(99) = RCT(62)*V(48)
! B(100) = dA(63)/dV(33)
  B(100) = RCT(63)*V(53)
! B(101) = dA(63)/dV(53)
  B(101) = RCT(63)*V(33)
! B(102) = dA(64)/dV(31)
  B(102) = RCT(64)*V(53)
! B(103) = dA(64)/dV(53)
  B(103) = RCT(64)*V(31)
! B(104) = dA(65)/dV(34)
  B(104) = RCT(65)*V(53)
! B(105) = dA(65)/dV(53)
  B(105) = RCT(65)*V(34)
! B(106) = dA(66)/dV(30)
  B(106) = RCT(66)*V(53)
! B(107) = dA(66)/dV(53)
  B(107) = RCT(66)*V(30)
! B(108) = dA(67)/dV(21)
  B(108) = RCT(67)*V(53)
! B(109) = dA(67)/dV(53)
  B(109) = RCT(67)*V(21)
! B(110) = dA(68)/dV(42)
  B(110) = RCT(68)*V(53)
! B(111) = dA(68)/dV(53)
  B(111) = RCT(68)*V(42)
! B(112) = dA(69)/dV(18)
  B(112) = RCT(69)*V(53)
! B(113) = dA(69)/dV(53)
  B(113) = RCT(69)*V(18)
! B(114) = dA(70)/dV(23)
  B(114) = RCT(70)*V(53)
! B(115) = dA(70)/dV(53)
  B(115) = RCT(70)*V(23)
! B(116) = dA(71)/dV(47)
  B(116) = RCT(71)*V(53)
! B(117) = dA(71)/dV(53)
  B(117) = RCT(71)*V(47)
! B(118) = dA(72)/dV(27)
  B(118) = RCT(72)*V(53)
! B(119) = dA(72)/dV(53)
  B(119) = RCT(72)*V(27)
! B(120) = dA(73)/dV(55)
  B(120) = RCT(73)*V(59)
! B(121) = dA(73)/dV(59)
  B(121) = RCT(73)*V(55)
! B(122) = dA(74)/dV(23)
  B(122) = RCT(74)
! B(123) = dA(75)/dV(45)
  B(123) = RCT(75)*V(59)
! B(124) = dA(75)/dV(59)
  B(124) = RCT(75)*V(45)
! B(125) = dA(76)/dV(12)
  B(125) = RCT(76)
! B(126) = dA(77)/dV(52)
  B(126) = RCT(77)*V(58)
! B(127) = dA(77)/dV(58)
  B(127) = RCT(77)*V(52)
! B(128) = dA(78)/dV(56)
  B(128) = RCT(78)*V(58)
! B(129) = dA(78)/dV(58)
  B(129) = RCT(78)*V(56)
! B(130) = dA(79)/dV(41)
  B(130) = RCT(79)*V(58)
! B(131) = dA(79)/dV(58)
  B(131) = RCT(79)*V(41)
! B(132) = dA(80)/dV(44)
  B(132) = RCT(80)*V(58)
! B(133) = dA(80)/dV(58)
  B(133) = RCT(80)*V(44)
! B(134) = dA(81)/dV(40)
  B(134) = RCT(81)*V(58)
! B(135) = dA(81)/dV(58)
  B(135) = RCT(81)*V(40)
! B(136) = dA(82)/dV(37)
  B(136) = RCT(82)*V(58)
! B(137) = dA(82)/dV(58)
  B(137) = RCT(82)*V(37)
! B(138) = dA(83)/dV(49)
  B(138) = RCT(83)*V(58)
! B(139) = dA(83)/dV(58)
  B(139) = RCT(83)*V(49)
! B(140) = dA(84)/dV(55)
  B(140) = RCT(84)*V(58)
! B(141) = dA(84)/dV(58)
  B(141) = RCT(84)*V(55)
! B(142) = dA(85)/dV(45)
  B(142) = RCT(85)*V(58)
! B(143) = dA(85)/dV(58)
  B(143) = RCT(85)*V(45)
! B(144) = dA(86)/dV(35)
  B(144) = RCT(86)*V(58)
! B(145) = dA(86)/dV(58)
  B(145) = RCT(86)*V(35)
! B(146) = dA(87)/dV(36)
  B(146) = RCT(87)*V(58)
! B(147) = dA(87)/dV(58)
  B(147) = RCT(87)*V(36)
! B(148) = dA(88)/dV(57)
  B(148) = RCT(88)*V(58)
! B(149) = dA(88)/dV(58)
  B(149) = RCT(88)*V(57)
! B(150) = dA(89)/dV(50)
  B(150) = RCT(89)*V(58)
! B(151) = dA(89)/dV(58)
  B(151) = RCT(89)*V(50)
! B(152) = dA(90)/dV(38)
  B(152) = RCT(90)*V(58)
! B(153) = dA(90)/dV(58)
  B(153) = RCT(90)*V(38)
! B(154) = dA(91)/dV(43)
  B(154) = RCT(91)*V(54)
! B(155) = dA(91)/dV(54)
  B(155) = RCT(91)*V(43)
! B(156) = dA(92)/dV(48)
  B(156) = RCT(92)*V(54)
! B(157) = dA(92)/dV(54)
  B(157) = RCT(92)*V(48)
! B(158) = dA(93)/dV(31)
  B(158) = RCT(93)*V(54)
! B(159) = dA(93)/dV(54)
  B(159) = RCT(93)*V(31)
! B(160) = dA(94)/dV(34)
  B(160) = RCT(94)*V(54)
! B(161) = dA(94)/dV(54)
  B(161) = RCT(94)*V(34)
! B(162) = dA(95)/dV(30)
  B(162) = RCT(95)*V(54)
! B(163) = dA(95)/dV(54)
  B(163) = RCT(95)*V(30)
! B(164) = dA(96)/dV(22)
  B(164) = RCT(96)*V(54)
! B(165) = dA(96)/dV(54)
  B(165) = RCT(96)*V(22)
! B(166) = dA(97)/dV(24)
  B(166) = RCT(97)*V(54)
! B(167) = dA(97)/dV(54)
  B(167) = RCT(97)*V(24)
! B(168) = dA(98)/dV(28)
  B(168) = RCT(98)*V(54)
! B(169) = dA(98)/dV(54)
  B(169) = RCT(98)*V(28)
! B(170) = dA(99)/dV(29)
  B(170) = RCT(99)*V(54)
! B(171) = dA(99)/dV(54)
  B(171) = RCT(99)*V(29)
! B(172) = dA(100)/dV(27)
  B(172) = RCT(100)*V(54)
! B(173) = dA(100)/dV(54)
  B(173) = RCT(100)*V(27)
! B(174) = dA(101)/dV(24)
  B(174) = RCT(101)*V(46)
! B(175) = dA(101)/dV(46)
  B(175) = RCT(101)*V(24)
! B(176) = dA(102)/dV(28)
  B(176) = RCT(102)*V(46)
! B(177) = dA(102)/dV(46)
  B(177) = RCT(102)*V(28)
! B(178) = dA(103)/dV(29)
  B(178) = RCT(103)*V(46)
! B(179) = dA(103)/dV(46)
  B(179) = RCT(103)*V(29)
! B(180) = dA(104)/dV(27)
  B(180) = RCT(104)*V(46)
! B(181) = dA(104)/dV(46)
  B(181) = RCT(104)*V(27)
! B(182) = dA(105)/dV(51)
  B(182) = RCT(105)*V(52)
! B(183) = dA(105)/dV(52)
  B(183) = RCT(105)*V(51)
! B(184) = dA(106)/dV(51)
  B(184) = RCT(106)*V(57)
! B(185) = dA(106)/dV(57)
  B(185) = RCT(106)*V(51)
! B(186) = dA(107)/dV(51)
  B(186) = RCT(107)*V(56)
! B(187) = dA(107)/dV(56)
  B(187) = RCT(107)*V(51)
! B(188) = dA(108)/dV(41)
  B(188) = RCT(108)*V(51)
! B(189) = dA(108)/dV(51)
  B(189) = RCT(108)*V(41)
! B(190) = dA(109)/dV(44)
  B(190) = RCT(109)*V(51)
! B(191) = dA(109)/dV(51)
  B(191) = RCT(109)*V(44)
! B(192) = dA(110)/dV(40)
  B(192) = RCT(110)*V(51)
! B(193) = dA(110)/dV(51)
  B(193) = RCT(110)*V(40)
! B(194) = dA(111)/dV(37)
  B(194) = RCT(111)*V(51)
! B(195) = dA(111)/dV(51)
  B(195) = RCT(111)*V(37)
! B(196) = dA(112)/dV(49)
  B(196) = RCT(112)*V(51)
! B(197) = dA(112)/dV(51)
  B(197) = RCT(112)*V(49)
! B(198) = dA(113)/dV(50)
  B(198) = RCT(113)*V(51)
! B(199) = dA(113)/dV(51)
  B(199) = RCT(113)*V(50)
! B(200) = dA(114)/dV(51)
  B(200) = RCT(114)*V(55)
! B(201) = dA(114)/dV(55)
  B(201) = RCT(114)*V(51)
! B(202) = dA(115)/dV(35)
  B(202) = RCT(115)*V(51)
! B(203) = dA(115)/dV(51)
  B(203) = RCT(115)*V(35)
! B(204) = dA(116)/dV(36)
  B(204) = RCT(116)*V(51)
! B(205) = dA(116)/dV(51)
  B(205) = RCT(116)*V(36)
! B(206) = dA(117)/dV(45)
  B(206) = RCT(117)*V(51)
! B(207) = dA(117)/dV(51)
  B(207) = RCT(117)*V(45)
! B(208) = dA(118)/dV(38)
  B(208) = RCT(118)*V(51)
! B(209) = dA(118)/dV(51)
  B(209) = RCT(118)*V(38)
! B(210) = dA(119)/dV(52)
  B(210) = RCT(119)*2*V(52)
! B(211) = dA(120)/dV(52)
  B(211) = RCT(120)*V(57)
! B(212) = dA(120)/dV(57)
  B(212) = RCT(120)*V(52)
! B(213) = dA(121)/dV(52)
  B(213) = RCT(121)*V(56)
! B(214) = dA(121)/dV(56)
  B(214) = RCT(121)*V(52)
! B(215) = dA(122)/dV(41)
  B(215) = RCT(122)*V(52)
! B(216) = dA(122)/dV(52)
  B(216) = RCT(122)*V(41)
! B(217) = dA(123)/dV(44)
  B(217) = RCT(123)*V(52)
! B(218) = dA(123)/dV(52)
  B(218) = RCT(123)*V(44)
! B(219) = dA(124)/dV(40)
  B(219) = RCT(124)*V(52)
! B(220) = dA(124)/dV(52)
  B(220) = RCT(124)*V(40)
! B(221) = dA(125)/dV(37)
  B(221) = RCT(125)*V(52)
! B(222) = dA(125)/dV(52)
  B(222) = RCT(125)*V(37)
! B(223) = dA(126)/dV(49)
  B(223) = RCT(126)*V(52)
! B(224) = dA(126)/dV(52)
  B(224) = RCT(126)*V(49)
! B(225) = dA(127)/dV(50)
  B(225) = RCT(127)*V(52)
! B(226) = dA(127)/dV(52)
  B(226) = RCT(127)*V(50)
! B(227) = dA(128)/dV(52)
  B(227) = RCT(128)*V(55)
! B(228) = dA(128)/dV(55)
  B(228) = RCT(128)*V(52)
! B(229) = dA(129)/dV(35)
  B(229) = RCT(129)*V(52)
! B(230) = dA(129)/dV(52)
  B(230) = RCT(129)*V(35)
! B(231) = dA(130)/dV(36)
  B(231) = RCT(130)*V(52)
! B(232) = dA(130)/dV(52)
  B(232) = RCT(130)*V(36)
! B(233) = dA(131)/dV(45)
  B(233) = RCT(131)*V(52)
! B(234) = dA(131)/dV(52)
  B(234) = RCT(131)*V(45)
! B(235) = dA(132)/dV(55)
  B(235) = RCT(132)*V(57)
! B(236) = dA(132)/dV(57)
  B(236) = RCT(132)*V(55)
! B(237) = dA(133)/dV(55)
  B(237) = RCT(133)*V(56)
! B(238) = dA(133)/dV(56)
  B(238) = RCT(133)*V(55)
! B(239) = dA(134)/dV(41)
  B(239) = RCT(134)*V(55)
! B(240) = dA(134)/dV(55)
  B(240) = RCT(134)*V(41)
! B(241) = dA(135)/dV(44)
  B(241) = RCT(135)*V(55)
! B(242) = dA(135)/dV(55)
  B(242) = RCT(135)*V(44)
! B(243) = dA(136)/dV(40)
  B(243) = RCT(136)*V(55)
! B(244) = dA(136)/dV(55)
  B(244) = RCT(136)*V(40)
! B(245) = dA(137)/dV(37)
  B(245) = RCT(137)*V(55)
! B(246) = dA(137)/dV(55)
  B(246) = RCT(137)*V(37)
! B(247) = dA(138)/dV(49)
  B(247) = RCT(138)*V(55)
! B(248) = dA(138)/dV(55)
  B(248) = RCT(138)*V(49)
! B(249) = dA(139)/dV(50)
  B(249) = RCT(139)*V(55)
! B(250) = dA(139)/dV(55)
  B(250) = RCT(139)*V(50)
! B(251) = dA(140)/dV(55)
  B(251) = RCT(140)*2*V(55)
! B(252) = dA(141)/dV(35)
  B(252) = RCT(141)*V(55)
! B(253) = dA(141)/dV(55)
  B(253) = RCT(141)*V(35)
! B(254) = dA(142)/dV(36)
  B(254) = RCT(142)*V(55)
! B(255) = dA(142)/dV(55)
  B(255) = RCT(142)*V(36)
! B(256) = dA(143)/dV(45)
  B(256) = RCT(143)*V(55)
! B(257) = dA(143)/dV(55)
  B(257) = RCT(143)*V(45)
! B(258) = dA(144)/dV(39)
  B(258) = RCT(144)*V(51)
! B(259) = dA(144)/dV(51)
  B(259) = RCT(144)*V(39)
! B(260) = dA(145)/dV(39)
  B(260) = RCT(145)*V(52)
! B(261) = dA(145)/dV(52)
  B(261) = RCT(145)*V(39)
! B(262) = dA(146)/dV(39)
  B(262) = RCT(146)*V(55)
! B(263) = dA(146)/dV(55)
  B(263) = RCT(146)*V(39)
! B(264) = dA(147)/dV(39)
  B(264) = RCT(147)*2*V(39)
! B(265) = dA(148)/dV(39)
  B(265) = RCT(148)*V(58)
! B(266) = dA(148)/dV(58)
  B(266) = RCT(148)*V(39)
! B(267) = dA(149)/dV(32)
  B(267) = RCT(149)*V(59)
! B(268) = dA(149)/dV(59)
  B(268) = RCT(149)*V(32)
! B(269) = dA(150)/dV(32)
  B(269) = RCT(150)*V(51)
! B(270) = dA(150)/dV(51)
  B(270) = RCT(150)*V(32)
! B(271) = dA(151)/dV(32)
  B(271) = RCT(151)*V(52)
! B(272) = dA(151)/dV(52)
  B(272) = RCT(151)*V(32)
! B(273) = dA(152)/dV(32)
  B(273) = RCT(152)*V(55)
! B(274) = dA(152)/dV(55)
  B(274) = RCT(152)*V(32)
! B(275) = dA(153)/dV(32)
  B(275) = RCT(153)*2*V(32)
! B(276) = dA(154)/dV(38)
  B(276) = RCT(154)*V(52)
! B(277) = dA(154)/dV(52)
  B(277) = RCT(154)*V(38)
! B(278) = dA(155)/dV(38)
  B(278) = RCT(155)*V(55)
! B(279) = dA(155)/dV(55)
  B(279) = RCT(155)*V(38)
! B(280) = dA(156)/dV(38)
  B(280) = RCT(156)*2*V(38)

! Construct the Jacobian terms from B's
! JVS(1) = Jac_FULL(1,1)
  JVS(1) = 0
! JVS(2) = Jac_FULL(1,5)
  JVS(2) = B(70)
! JVS(3) = Jac_FULL(1,53)
  JVS(3) = B(71)
! JVS(4) = Jac_FULL(2,2)
  JVS(4) = 0
! JVS(5) = Jac_FULL(2,24)
  JVS(5) = 0.4*B(174)
! JVS(6) = Jac_FULL(2,27)
  JVS(6) = 0.2*B(180)
! JVS(7) = Jac_FULL(2,28)
  JVS(7) = 0.2*B(176)
! JVS(8) = Jac_FULL(2,29)
  JVS(8) = 0.06*B(178)
! JVS(9) = Jac_FULL(2,46)
  JVS(9) = 0.4*B(175)+0.2*B(177)+0.06*B(179)+0.2*B(181)
! JVS(10) = Jac_FULL(3,3)
  JVS(10) = 0
! JVS(11) = Jac_FULL(3,27)
  JVS(11) = 0.2*B(180)
! JVS(12) = Jac_FULL(3,28)
  JVS(12) = 0.2*B(176)
! JVS(13) = Jac_FULL(3,29)
  JVS(13) = 0.29*B(178)
! JVS(14) = Jac_FULL(3,37)
  JVS(14) = 0.5*B(245)
! JVS(15) = Jac_FULL(3,38)
  JVS(15) = 0.5*B(278)
! JVS(16) = Jac_FULL(3,40)
  JVS(16) = 0.5*B(243)
! JVS(17) = Jac_FULL(3,41)
  JVS(17) = 0.5*B(239)
! JVS(18) = Jac_FULL(3,44)
  JVS(18) = 0.5*B(241)
! JVS(19) = Jac_FULL(3,45)
  JVS(19) = 0.5*B(233)
! JVS(20) = Jac_FULL(3,46)
  JVS(20) = 0.2*B(177)+0.29*B(179)+0.2*B(181)
! JVS(21) = Jac_FULL(3,49)
  JVS(21) = 0.5*B(247)
! JVS(22) = Jac_FULL(3,50)
  JVS(22) = 0.5*B(249)
! JVS(23) = Jac_FULL(3,52)
  JVS(23) = 0.5*B(227)+0.5*B(234)
! JVS(24) = Jac_FULL(3,55)
  JVS(24) = 0.5*B(228)+0.5*B(235)+0.5*B(237)+0.5*B(240)+0.5*B(242)+0.5*B(244)+0.5*B(246)+0.5*B(248)+0.5*B(250)+0.5&
              &*B(279)
! JVS(25) = Jac_FULL(3,56)
  JVS(25) = 0.5*B(238)
! JVS(26) = Jac_FULL(3,57)
  JVS(26) = 0.5*B(236)
! JVS(27) = Jac_FULL(4,4)
  JVS(27) = 0
! JVS(28) = Jac_FULL(4,26)
  JVS(28) = B(72)
! JVS(29) = Jac_FULL(4,53)
  JVS(29) = B(73)
! JVS(30) = Jac_FULL(5,5)
  JVS(30) = -B(70)
! JVS(31) = Jac_FULL(5,53)
  JVS(31) = -B(71)
! JVS(32) = Jac_FULL(6,6)
  JVS(32) = -B(76)
! JVS(33) = Jac_FULL(6,53)
  JVS(33) = -B(77)
! JVS(34) = Jac_FULL(7,7)
  JVS(34) = -B(26)-B(28)
! JVS(35) = Jac_FULL(7,46)
  JVS(35) = B(2)
! JVS(36) = Jac_FULL(8,8)
  JVS(36) = -B(80)
! JVS(37) = Jac_FULL(8,53)
  JVS(37) = -B(81)
! JVS(38) = Jac_FULL(9,9)
  JVS(38) = -B(82)
! JVS(39) = Jac_FULL(9,53)
  JVS(39) = -B(83)
! JVS(40) = Jac_FULL(10,10)
  JVS(40) = -B(90)
! JVS(41) = Jac_FULL(10,53)
  JVS(41) = -B(91)
! JVS(42) = Jac_FULL(11,11)
  JVS(42) = -B(92)
! JVS(43) = Jac_FULL(11,53)
  JVS(43) = -B(93)
! JVS(44) = Jac_FULL(12,12)
  JVS(44) = -B(125)
! JVS(45) = Jac_FULL(12,45)
  JVS(45) = B(123)
! JVS(46) = Jac_FULL(12,59)
  JVS(46) = B(124)
! JVS(47) = Jac_FULL(13,13)
  JVS(47) = -B(4)
! JVS(48) = Jac_FULL(13,53)
  JVS(48) = B(46)
! JVS(49) = Jac_FULL(13,58)
  JVS(49) = B(47)
! JVS(50) = Jac_FULL(14,14)
  JVS(50) = -B(9)-B(44)
! JVS(51) = Jac_FULL(14,51)
  JVS(51) = B(41)+B(42)
! JVS(52) = Jac_FULL(14,53)
  JVS(52) = -B(45)
! JVS(53) = Jac_FULL(15,15)
  JVS(53) = -B(60)-B(61)
! JVS(54) = Jac_FULL(15,54)
  JVS(54) = B(58)
! JVS(55) = Jac_FULL(15,59)
  JVS(55) = B(59)
! JVS(56) = Jac_FULL(16,16)
  JVS(56) = -B(78)
! JVS(57) = Jac_FULL(16,53)
  JVS(57) = -B(79)
! JVS(58) = Jac_FULL(17,17)
  JVS(58) = -B(74)
! JVS(59) = Jac_FULL(17,28)
  JVS(59) = 0.06*B(176)
! JVS(60) = Jac_FULL(17,29)
  JVS(60) = 0.09*B(178)
! JVS(61) = Jac_FULL(17,46)
  JVS(61) = 0.06*B(177)+0.09*B(179)
! JVS(62) = Jac_FULL(17,53)
  JVS(62) = -B(75)
! JVS(63) = Jac_FULL(18,18)
  JVS(63) = -B(15)-B(112)
! JVS(64) = Jac_FULL(18,51)
  JVS(64) = B(200)
! JVS(65) = Jac_FULL(18,53)
  JVS(65) = -B(113)
! JVS(66) = Jac_FULL(18,55)
  JVS(66) = B(201)
! JVS(67) = Jac_FULL(19,7)
  JVS(67) = B(26)
! JVS(68) = Jac_FULL(19,19)
  JVS(68) = -B(22)-B(24)
! JVS(69) = Jac_FULL(19,46)
  JVS(69) = B(3)
! JVS(70) = Jac_FULL(19,54)
  JVS(70) = B(8)
! JVS(71) = Jac_FULL(19,59)
  JVS(71) = B(1)-B(25)
! JVS(72) = Jac_FULL(20,20)
  JVS(72) = -B(6)-B(40)-B(66)
! JVS(73) = Jac_FULL(20,51)
  JVS(73) = B(38)
! JVS(74) = Jac_FULL(20,53)
  JVS(74) = -B(67)
! JVS(75) = Jac_FULL(20,59)
  JVS(75) = B(39)
! JVS(76) = Jac_FULL(21,21)
  JVS(76) = -B(13)-B(108)
! JVS(77) = Jac_FULL(21,51)
  JVS(77) = B(182)
! JVS(78) = Jac_FULL(21,52)
  JVS(78) = B(183)
! JVS(79) = Jac_FULL(21,53)
  JVS(79) = -B(109)
! JVS(80) = Jac_FULL(22,10)
  JVS(80) = 0.25*B(90)
! JVS(81) = Jac_FULL(22,11)
  JVS(81) = 0.17*B(92)
! JVS(82) = Jac_FULL(22,22)
  JVS(82) = -B(94)-0.5*B(164)
! JVS(83) = Jac_FULL(22,53)
  JVS(83) = 0.25*B(91)+0.17*B(93)-B(95)
! JVS(84) = Jac_FULL(22,54)
  JVS(84) = -0.5*B(165)
! JVS(85) = Jac_FULL(23,23)
  JVS(85) = -B(114)-B(122)
! JVS(86) = Jac_FULL(23,53)
  JVS(86) = -B(115)
! JVS(87) = Jac_FULL(23,55)
  JVS(87) = B(120)
! JVS(88) = Jac_FULL(23,59)
  JVS(88) = B(121)
! JVS(89) = Jac_FULL(24,24)
  JVS(89) = -B(84)-B(166)-B(174)
! JVS(90) = Jac_FULL(24,46)
  JVS(90) = -B(175)
! JVS(91) = Jac_FULL(24,53)
  JVS(91) = -B(85)
! JVS(92) = Jac_FULL(24,54)
  JVS(92) = -B(167)
! JVS(93) = Jac_FULL(25,15)
  JVS(93) = 2*B(61)
! JVS(94) = Jac_FULL(25,22)
  JVS(94) = B(164)
! JVS(95) = Jac_FULL(25,25)
  JVS(95) = -B(5)-B(64)
! JVS(96) = Jac_FULL(25,30)
  JVS(96) = B(162)
! JVS(97) = Jac_FULL(25,31)
  JVS(97) = B(158)
! JVS(98) = Jac_FULL(25,34)
  JVS(98) = B(160)
! JVS(99) = Jac_FULL(25,43)
  JVS(99) = B(154)
! JVS(100) = Jac_FULL(25,48)
  JVS(100) = B(156)
! JVS(101) = Jac_FULL(25,51)
  JVS(101) = B(56)
! JVS(102) = Jac_FULL(25,53)
  JVS(102) = B(62)-B(65)
! JVS(103) = Jac_FULL(25,54)
  JVS(103) = B(57)+B(155)+B(157)+B(159)+B(161)+B(163)+B(165)
! JVS(104) = Jac_FULL(25,59)
  JVS(104) = B(63)
! JVS(105) = Jac_FULL(26,24)
  JVS(105) = 0.42*B(174)
! JVS(106) = Jac_FULL(26,26)
  JVS(106) = -B(72)
! JVS(107) = Jac_FULL(26,27)
  JVS(107) = 0.33*B(180)
! JVS(108) = Jac_FULL(26,28)
  JVS(108) = 0.33*B(176)
! JVS(109) = Jac_FULL(26,29)
  JVS(109) = 0.23*B(178)
! JVS(110) = Jac_FULL(26,31)
  JVS(110) = 1.87*B(17)+1.55*B(18)+2*B(102)+2*B(158)
! JVS(111) = Jac_FULL(26,34)
  JVS(111) = B(19)+B(104)+B(160)
! JVS(112) = Jac_FULL(26,43)
  JVS(112) = B(10)+B(11)+B(96)+B(154)
! JVS(113) = Jac_FULL(26,45)
  JVS(113) = 0.95*B(142)+0.475*B(233)+0.95*B(256)
! JVS(114) = Jac_FULL(26,46)
  JVS(114) = 0.42*B(175)+0.33*B(177)+0.23*B(179)+0.33*B(181)
! JVS(115) = Jac_FULL(26,48)
  JVS(115) = B(12)
! JVS(116) = Jac_FULL(26,52)
  JVS(116) = 0.475*B(234)
! JVS(117) = Jac_FULL(26,53)
  JVS(117) = -B(73)+B(97)+2*B(103)+B(105)
! JVS(118) = Jac_FULL(26,54)
  JVS(118) = B(155)+2*B(159)+B(161)
! JVS(119) = Jac_FULL(26,55)
  JVS(119) = 0.95*B(257)
! JVS(120) = Jac_FULL(26,58)
  JVS(120) = 0.95*B(143)
! JVS(121) = Jac_FULL(27,27)
  JVS(121) = -B(118)-B(172)-B(180)
! JVS(122) = Jac_FULL(27,46)
  JVS(122) = -B(181)
! JVS(123) = Jac_FULL(27,53)
  JVS(123) = -B(119)
! JVS(124) = Jac_FULL(27,54)
  JVS(124) = -B(173)
! JVS(125) = Jac_FULL(28,28)
  JVS(125) = -B(86)-B(168)-B(176)
! JVS(126) = Jac_FULL(28,46)
  JVS(126) = -B(177)
! JVS(127) = Jac_FULL(28,53)
  JVS(127) = -B(87)
! JVS(128) = Jac_FULL(28,54)
  JVS(128) = -B(169)
! JVS(129) = Jac_FULL(29,29)
  JVS(129) = -B(88)-B(170)-B(178)
! JVS(130) = Jac_FULL(29,46)
  JVS(130) = -B(179)
! JVS(131) = Jac_FULL(29,53)
  JVS(131) = -B(89)
! JVS(132) = Jac_FULL(29,54)
  JVS(132) = -B(171)
! JVS(133) = Jac_FULL(30,30)
  JVS(133) = -B(20)-B(106)-B(162)
! JVS(134) = Jac_FULL(30,35)
  JVS(134) = 0.7*B(144)+0.7*B(229)+B(252)
! JVS(135) = Jac_FULL(30,36)
  JVS(135) = 0.806*B(146)+0.806*B(231)+B(254)
! JVS(136) = Jac_FULL(30,52)
  JVS(136) = 0.7*B(230)+0.806*B(232)
! JVS(137) = Jac_FULL(30,53)
  JVS(137) = -B(107)
! JVS(138) = Jac_FULL(30,54)
  JVS(138) = -B(163)
! JVS(139) = Jac_FULL(30,55)
  JVS(139) = B(253)+B(255)
! JVS(140) = Jac_FULL(30,58)
  JVS(140) = 0.7*B(145)+0.806*B(147)
! JVS(141) = Jac_FULL(31,31)
  JVS(141) = -B(17)-B(18)-B(102)-B(158)
! JVS(142) = Jac_FULL(31,35)
  JVS(142) = 0.16*B(144)+0.16*B(229)+0.2*B(252)
! JVS(143) = Jac_FULL(31,45)
  JVS(143) = 0.89*B(142)+0.445*B(233)+0.89*B(256)
! JVS(144) = Jac_FULL(31,52)
  JVS(144) = 0.16*B(230)+0.445*B(234)
! JVS(145) = Jac_FULL(31,53)
  JVS(145) = -B(103)
! JVS(146) = Jac_FULL(31,54)
  JVS(146) = -B(159)
! JVS(147) = Jac_FULL(31,55)
  JVS(147) = 0.2*B(253)+0.89*B(257)
! JVS(148) = Jac_FULL(31,58)
  JVS(148) = 0.89*B(143)+0.16*B(145)
! JVS(149) = Jac_FULL(32,22)
  JVS(149) = B(164)
! JVS(150) = Jac_FULL(32,32)
  JVS(150) = -B(267)-B(269)-B(271)-B(273)-2*B(275)
! JVS(151) = Jac_FULL(32,51)
  JVS(151) = -B(270)
! JVS(152) = Jac_FULL(32,52)
  JVS(152) = -B(272)
! JVS(153) = Jac_FULL(32,53)
  JVS(153) = 0
! JVS(154) = Jac_FULL(32,54)
  JVS(154) = B(165)
! JVS(155) = Jac_FULL(32,55)
  JVS(155) = -B(274)
! JVS(156) = Jac_FULL(32,59)
  JVS(156) = -B(268)
! JVS(157) = Jac_FULL(33,16)
  JVS(157) = 0.025*B(78)
! JVS(158) = Jac_FULL(33,29)
  JVS(158) = 0.1*B(178)
! JVS(159) = Jac_FULL(33,33)
  JVS(159) = -B(16)-B(100)
! JVS(160) = Jac_FULL(33,41)
  JVS(160) = 0.69*B(130)+0.75*B(215)+0.86*B(239)
! JVS(161) = Jac_FULL(33,44)
  JVS(161) = 1.06*B(132)+1.39*B(217)+0.9*B(241)
! JVS(162) = Jac_FULL(33,46)
  JVS(162) = 0.1*B(179)
! JVS(163) = Jac_FULL(33,47)
  JVS(163) = 0.8*B(21)
! JVS(164) = Jac_FULL(33,49)
  JVS(164) = 0.1*B(138)+0.55*B(223)+0.55*B(247)
! JVS(165) = Jac_FULL(33,52)
  JVS(165) = 0.6*B(213)+0.75*B(216)+1.39*B(218)+0.55*B(224)
! JVS(166) = Jac_FULL(33,53)
  JVS(166) = 0.025*B(79)-B(101)
! JVS(167) = Jac_FULL(33,54)
  JVS(167) = 0
! JVS(168) = Jac_FULL(33,55)
  JVS(168) = 0.8*B(237)+0.86*B(240)+0.9*B(242)+0.55*B(248)
! JVS(169) = Jac_FULL(33,56)
  JVS(169) = 0.25*B(128)+0.6*B(214)+0.8*B(238)
! JVS(170) = Jac_FULL(33,58)
  JVS(170) = 0.25*B(129)+0.69*B(131)+1.06*B(133)+0.1*B(139)
! JVS(171) = Jac_FULL(34,34)
  JVS(171) = -B(19)-B(104)-B(160)
! JVS(172) = Jac_FULL(34,35)
  JVS(172) = 0.17*B(144)+0.17*B(229)+0.8*B(252)
! JVS(173) = Jac_FULL(34,36)
  JVS(173) = 0.45*B(146)+0.45*B(231)+B(254)
! JVS(174) = Jac_FULL(34,45)
  JVS(174) = 0.11*B(142)+0.055*B(233)+0.11*B(256)
! JVS(175) = Jac_FULL(34,50)
  JVS(175) = B(150)+0.75*B(225)+B(249)
! JVS(176) = Jac_FULL(34,52)
  JVS(176) = 0.75*B(226)+0.17*B(230)+0.45*B(232)+0.055*B(234)
! JVS(177) = Jac_FULL(34,53)
  JVS(177) = -B(105)
! JVS(178) = Jac_FULL(34,54)
  JVS(178) = -B(161)
! JVS(179) = Jac_FULL(34,55)
  JVS(179) = B(250)+0.8*B(253)+B(255)+0.11*B(257)
! JVS(180) = Jac_FULL(34,58)
  JVS(180) = 0.11*B(143)+0.17*B(145)+0.45*B(147)+B(151)
! JVS(181) = Jac_FULL(35,10)
  JVS(181) = 0.75*B(90)
! JVS(182) = Jac_FULL(35,35)
  JVS(182) = -B(144)-B(202)-B(229)-B(252)
! JVS(183) = Jac_FULL(35,51)
  JVS(183) = -B(203)
! JVS(184) = Jac_FULL(35,52)
  JVS(184) = -B(230)
! JVS(185) = Jac_FULL(35,53)
  JVS(185) = 0.75*B(91)
! JVS(186) = Jac_FULL(35,55)
  JVS(186) = -B(253)
! JVS(187) = Jac_FULL(35,58)
  JVS(187) = -B(145)
! JVS(188) = Jac_FULL(36,11)
  JVS(188) = 0.83*B(92)
! JVS(189) = Jac_FULL(36,36)
  JVS(189) = -B(146)-B(204)-B(231)-B(254)
! JVS(190) = Jac_FULL(36,51)
  JVS(190) = -B(205)
! JVS(191) = Jac_FULL(36,52)
  JVS(191) = -B(232)
! JVS(192) = Jac_FULL(36,53)
  JVS(192) = 0.83*B(93)
! JVS(193) = Jac_FULL(36,55)
  JVS(193) = -B(255)
! JVS(194) = Jac_FULL(36,58)
  JVS(194) = -B(147)
! JVS(195) = Jac_FULL(37,27)
  JVS(195) = B(118)
! JVS(196) = Jac_FULL(37,28)
  JVS(196) = B(86)
! JVS(197) = Jac_FULL(37,37)
  JVS(197) = -B(136)-B(194)-B(221)-B(245)
! JVS(198) = Jac_FULL(37,46)
  JVS(198) = 0
! JVS(199) = Jac_FULL(37,51)
  JVS(199) = -B(195)
! JVS(200) = Jac_FULL(37,52)
  JVS(200) = -B(222)
! JVS(201) = Jac_FULL(37,53)
  JVS(201) = B(87)+B(119)
! JVS(202) = Jac_FULL(37,54)
  JVS(202) = 0
! JVS(203) = Jac_FULL(37,55)
  JVS(203) = -B(246)
! JVS(204) = Jac_FULL(37,58)
  JVS(204) = -B(137)
! JVS(205) = Jac_FULL(38,24)
  JVS(205) = B(166)
! JVS(206) = Jac_FULL(38,27)
  JVS(206) = B(172)
! JVS(207) = Jac_FULL(38,28)
  JVS(207) = B(168)
! JVS(208) = Jac_FULL(38,29)
  JVS(208) = B(170)
! JVS(209) = Jac_FULL(38,38)
  JVS(209) = -B(152)-B(208)-B(276)-B(278)-2*B(280)
! JVS(210) = Jac_FULL(38,46)
  JVS(210) = 0
! JVS(211) = Jac_FULL(38,51)
  JVS(211) = -B(209)
! JVS(212) = Jac_FULL(38,52)
  JVS(212) = -B(277)
! JVS(213) = Jac_FULL(38,53)
  JVS(213) = 0
! JVS(214) = Jac_FULL(38,54)
  JVS(214) = B(167)+B(169)+B(171)+B(173)
! JVS(215) = Jac_FULL(38,55)
  JVS(215) = -B(279)
! JVS(216) = Jac_FULL(38,58)
  JVS(216) = -B(153)
! JVS(217) = Jac_FULL(39,8)
  JVS(217) = 0.25*B(80)
! JVS(218) = Jac_FULL(39,9)
  JVS(218) = 0.75*B(82)
! JVS(219) = Jac_FULL(39,22)
  JVS(219) = 0.9*B(94)
! JVS(220) = Jac_FULL(39,23)
  JVS(220) = B(114)
! JVS(221) = Jac_FULL(39,39)
  JVS(221) = -B(258)-B(260)-B(262)-2*B(264)-B(265)
! JVS(222) = Jac_FULL(39,45)
  JVS(222) = 2*B(142)+B(233)+2*B(256)
! JVS(223) = Jac_FULL(39,51)
  JVS(223) = -B(259)
! JVS(224) = Jac_FULL(39,52)
  JVS(224) = B(234)-B(261)
! JVS(225) = Jac_FULL(39,53)
  JVS(225) = 0.25*B(81)+0.75*B(83)+0.9*B(95)+B(115)
! JVS(226) = Jac_FULL(39,54)
  JVS(226) = 0
! JVS(227) = Jac_FULL(39,55)
  JVS(227) = 2*B(257)-B(263)
! JVS(228) = Jac_FULL(39,58)
  JVS(228) = 2*B(143)-B(266)
! JVS(229) = Jac_FULL(39,59)
  JVS(229) = 0
! JVS(230) = Jac_FULL(40,24)
  JVS(230) = B(84)
! JVS(231) = Jac_FULL(40,40)
  JVS(231) = -B(134)-B(192)-B(219)-B(243)
! JVS(232) = Jac_FULL(40,46)
  JVS(232) = 0
! JVS(233) = Jac_FULL(40,51)
  JVS(233) = -B(193)
! JVS(234) = Jac_FULL(40,52)
  JVS(234) = -B(220)
! JVS(235) = Jac_FULL(40,53)
  JVS(235) = B(85)
! JVS(236) = Jac_FULL(40,54)
  JVS(236) = 0
! JVS(237) = Jac_FULL(40,55)
  JVS(237) = -B(244)
! JVS(238) = Jac_FULL(40,58)
  JVS(238) = -B(135)
! JVS(239) = Jac_FULL(41,8)
  JVS(239) = B(80)
! JVS(240) = Jac_FULL(41,41)
  JVS(240) = -B(130)-B(188)-B(215)-B(239)
! JVS(241) = Jac_FULL(41,51)
  JVS(241) = -B(189)
! JVS(242) = Jac_FULL(41,52)
  JVS(242) = -B(216)
! JVS(243) = Jac_FULL(41,53)
  JVS(243) = B(81)
! JVS(244) = Jac_FULL(41,55)
  JVS(244) = -B(240)
! JVS(245) = Jac_FULL(41,58)
  JVS(245) = -B(131)
! JVS(246) = Jac_FULL(42,32)
  JVS(246) = B(269)
! JVS(247) = Jac_FULL(42,35)
  JVS(247) = B(202)
! JVS(248) = Jac_FULL(42,36)
  JVS(248) = B(204)
! JVS(249) = Jac_FULL(42,37)
  JVS(249) = B(194)
! JVS(250) = Jac_FULL(42,39)
  JVS(250) = B(258)
! JVS(251) = Jac_FULL(42,40)
  JVS(251) = B(192)
! JVS(252) = Jac_FULL(42,41)
  JVS(252) = B(188)
! JVS(253) = Jac_FULL(42,42)
  JVS(253) = -B(14)-B(110)
! JVS(254) = Jac_FULL(42,44)
  JVS(254) = B(190)
! JVS(255) = Jac_FULL(42,45)
  JVS(255) = B(206)
! JVS(256) = Jac_FULL(42,46)
  JVS(256) = 0
! JVS(257) = Jac_FULL(42,49)
  JVS(257) = B(196)
! JVS(258) = Jac_FULL(42,50)
  JVS(258) = B(198)
! JVS(259) = Jac_FULL(42,51)
  JVS(259) = B(184)+B(186)+B(189)+B(191)+B(193)+B(195)+B(197)+B(199)+B(203)+B(205)+B(207)+B(259)+B(270)
! JVS(260) = Jac_FULL(42,52)
  JVS(260) = 0
! JVS(261) = Jac_FULL(42,53)
  JVS(261) = -B(111)
! JVS(262) = Jac_FULL(42,54)
  JVS(262) = 0
! JVS(263) = Jac_FULL(42,55)
  JVS(263) = 0
! JVS(264) = Jac_FULL(42,56)
  JVS(264) = B(187)
! JVS(265) = Jac_FULL(42,57)
  JVS(265) = B(185)
! JVS(266) = Jac_FULL(42,58)
  JVS(266) = 0
! JVS(267) = Jac_FULL(42,59)
  JVS(267) = 0
! JVS(268) = Jac_FULL(43,16)
  JVS(268) = 0.009*B(78)
! JVS(269) = Jac_FULL(43,21)
  JVS(269) = B(13)+0.5*B(108)
! JVS(270) = Jac_FULL(43,23)
  JVS(270) = B(114)
! JVS(271) = Jac_FULL(43,24)
  JVS(271) = B(174)
! JVS(272) = Jac_FULL(43,27)
  JVS(272) = 0.53*B(180)
! JVS(273) = Jac_FULL(43,28)
  JVS(273) = 0.53*B(176)
! JVS(274) = Jac_FULL(43,29)
  JVS(274) = 0.18*B(178)
! JVS(275) = Jac_FULL(43,31)
  JVS(275) = 0.13*B(17)+0.45*B(18)
! JVS(276) = Jac_FULL(43,32)
  JVS(276) = B(271)
! JVS(277) = Jac_FULL(43,35)
  JVS(277) = B(229)
! JVS(278) = Jac_FULL(43,36)
  JVS(278) = B(231)
! JVS(279) = Jac_FULL(43,37)
  JVS(279) = B(136)+1.25*B(221)+0.5*B(245)
! JVS(280) = Jac_FULL(43,38)
  JVS(280) = B(152)+1.75*B(276)+B(278)+2*B(280)
! JVS(281) = Jac_FULL(43,39)
  JVS(281) = B(260)
! JVS(282) = Jac_FULL(43,40)
  JVS(282) = 1.6*B(134)+1.55*B(219)+0.8*B(243)
! JVS(283) = Jac_FULL(43,41)
  JVS(283) = 0.77*B(215)
! JVS(284) = Jac_FULL(43,43)
  JVS(284) = -B(10)-B(11)-B(96)-B(154)
! JVS(285) = Jac_FULL(43,44)
  JVS(285) = 0.04*B(132)+0.8*B(217)
! JVS(286) = Jac_FULL(43,45)
  JVS(286) = 0.5*B(233)
! JVS(287) = Jac_FULL(43,46)
  JVS(287) = B(175)+0.53*B(177)+0.18*B(179)+0.53*B(181)
! JVS(288) = Jac_FULL(43,49)
  JVS(288) = 0.28*B(138)+0.89*B(223)+0.14*B(247)
! JVS(289) = Jac_FULL(43,50)
  JVS(289) = 0.75*B(225)
! JVS(290) = Jac_FULL(43,51)
  JVS(290) = 0
! JVS(291) = Jac_FULL(43,52)
  JVS(291) = B(126)+1.5*B(210)+0.75*B(211)+0.75*B(213)+0.77*B(216)+0.8*B(218)+1.55*B(220)+1.25*B(222)+0.89*B(224)+0.75&
               &*B(226)+B(227)+B(230)+B(232)+0.5*B(234)+B(261)+B(272)+1.75*B(277)
! JVS(292) = Jac_FULL(43,53)
  JVS(292) = 0.009*B(79)-B(97)+0.5*B(109)+B(115)
! JVS(293) = Jac_FULL(43,54)
  JVS(293) = -B(155)
! JVS(294) = Jac_FULL(43,55)
  JVS(294) = B(228)+0.8*B(244)+0.5*B(246)+0.14*B(248)+B(279)
! JVS(295) = Jac_FULL(43,56)
  JVS(295) = 0.09*B(128)+0.75*B(214)
! JVS(296) = Jac_FULL(43,57)
  JVS(296) = 0.75*B(212)
! JVS(297) = Jac_FULL(43,58)
  JVS(297) = B(127)+0.09*B(129)+0.04*B(133)+1.6*B(135)+B(137)+0.28*B(139)+B(153)
! JVS(298) = Jac_FULL(43,59)
  JVS(298) = 0
! JVS(299) = Jac_FULL(44,9)
  JVS(299) = B(82)
! JVS(300) = Jac_FULL(44,44)
  JVS(300) = -B(132)-B(190)-B(217)-B(241)
! JVS(301) = Jac_FULL(44,51)
  JVS(301) = -B(191)
! JVS(302) = Jac_FULL(44,52)
  JVS(302) = -B(218)
! JVS(303) = Jac_FULL(44,53)
  JVS(303) = B(83)
! JVS(304) = Jac_FULL(44,55)
  JVS(304) = -B(242)
! JVS(305) = Jac_FULL(44,58)
  JVS(305) = -B(133)
! JVS(306) = Jac_FULL(45,12)
  JVS(306) = B(125)
! JVS(307) = Jac_FULL(45,22)
  JVS(307) = 0.9*B(94)
! JVS(308) = Jac_FULL(45,30)
  JVS(308) = B(20)+B(106)+B(162)
! JVS(309) = Jac_FULL(45,35)
  JVS(309) = 0
! JVS(310) = Jac_FULL(45,36)
  JVS(310) = 0
! JVS(311) = Jac_FULL(45,45)
  JVS(311) = -B(123)-B(142)-B(206)-B(233)-B(256)
! JVS(312) = Jac_FULL(45,51)
  JVS(312) = -B(207)
! JVS(313) = Jac_FULL(45,52)
  JVS(313) = -B(234)
! JVS(314) = Jac_FULL(45,53)
  JVS(314) = 0.9*B(95)+B(107)
! JVS(315) = Jac_FULL(45,54)
  JVS(315) = B(163)
! JVS(316) = Jac_FULL(45,55)
  JVS(316) = -B(257)
! JVS(317) = Jac_FULL(45,58)
  JVS(317) = -B(143)
! JVS(318) = Jac_FULL(45,59)
  JVS(318) = -B(124)
! JVS(319) = Jac_FULL(46,19)
  JVS(319) = B(22)
! JVS(320) = Jac_FULL(46,24)
  JVS(320) = -B(174)
! JVS(321) = Jac_FULL(46,27)
  JVS(321) = -B(180)
! JVS(322) = Jac_FULL(46,28)
  JVS(322) = -B(176)
! JVS(323) = Jac_FULL(46,29)
  JVS(323) = -B(178)
! JVS(324) = Jac_FULL(46,46)
  JVS(324) = -B(2)-B(3)-B(30)-B(32)-B(34)-B(50)-B(175)-B(177)-B(179)-B(181)
! JVS(325) = Jac_FULL(46,51)
  JVS(325) = -B(35)
! JVS(326) = Jac_FULL(46,53)
  JVS(326) = -B(33)
! JVS(327) = Jac_FULL(46,54)
  JVS(327) = 0
! JVS(328) = Jac_FULL(46,58)
  JVS(328) = -B(31)
! JVS(329) = Jac_FULL(46,59)
  JVS(329) = -B(51)
! JVS(330) = Jac_FULL(47,32)
  JVS(330) = B(267)
! JVS(331) = Jac_FULL(47,38)
  JVS(331) = B(208)
! JVS(332) = Jac_FULL(47,41)
  JVS(332) = 0.08*B(130)
! JVS(333) = Jac_FULL(47,44)
  JVS(333) = 0.24*B(132)
! JVS(334) = Jac_FULL(47,46)
  JVS(334) = 0
! JVS(335) = Jac_FULL(47,47)
  JVS(335) = -B(21)-B(116)
! JVS(336) = Jac_FULL(47,51)
  JVS(336) = B(209)
! JVS(337) = Jac_FULL(47,52)
  JVS(337) = 0
! JVS(338) = Jac_FULL(47,53)
  JVS(338) = -B(117)
! JVS(339) = Jac_FULL(47,54)
  JVS(339) = 0
! JVS(340) = Jac_FULL(47,55)
  JVS(340) = 0
! JVS(341) = Jac_FULL(47,56)
  JVS(341) = 0.036*B(128)
! JVS(342) = Jac_FULL(47,58)
  JVS(342) = 0.036*B(129)+0.08*B(131)+0.24*B(133)
! JVS(343) = Jac_FULL(47,59)
  JVS(343) = B(268)
! JVS(344) = Jac_FULL(48,16)
  JVS(344) = 0.075*B(78)
! JVS(345) = Jac_FULL(48,27)
  JVS(345) = 0.5*B(180)
! JVS(346) = Jac_FULL(48,28)
  JVS(346) = 0.5*B(176)
! JVS(347) = Jac_FULL(48,29)
  JVS(347) = 0.72*B(178)
! JVS(348) = Jac_FULL(48,37)
  JVS(348) = B(136)+0.75*B(221)+B(245)
! JVS(349) = Jac_FULL(48,38)
  JVS(349) = B(152)+B(276)+B(278)+2*B(280)
! JVS(350) = Jac_FULL(48,40)
  JVS(350) = 0.2*B(134)+0.35*B(219)+0.6*B(243)
! JVS(351) = Jac_FULL(48,41)
  JVS(351) = 0.38*B(130)+0.41*B(215)+0.14*B(239)
! JVS(352) = Jac_FULL(48,42)
  JVS(352) = B(14)+0.5*B(110)
! JVS(353) = Jac_FULL(48,44)
  JVS(353) = 0.35*B(132)+0.46*B(217)+0.1*B(241)
! JVS(354) = Jac_FULL(48,45)
  JVS(354) = 0
! JVS(355) = Jac_FULL(48,46)
  JVS(355) = 0.5*B(177)+0.72*B(179)+0.5*B(181)
! JVS(356) = Jac_FULL(48,47)
  JVS(356) = 0.2*B(21)
! JVS(357) = Jac_FULL(48,48)
  JVS(357) = -B(12)-B(98)-B(156)
! JVS(358) = Jac_FULL(48,49)
  JVS(358) = 1.45*B(138)+0.725*B(223)+0.725*B(247)
! JVS(359) = Jac_FULL(48,50)
  JVS(359) = 0
! JVS(360) = Jac_FULL(48,51)
  JVS(360) = 0
! JVS(361) = Jac_FULL(48,52)
  JVS(361) = 0.75*B(211)+0.15*B(213)+0.41*B(216)+0.46*B(218)+0.35*B(220)+0.75*B(222)+0.725*B(224)+B(277)
! JVS(362) = Jac_FULL(48,53)
  JVS(362) = 0.075*B(79)-B(99)+0.5*B(111)
! JVS(363) = Jac_FULL(48,54)
  JVS(363) = -B(157)
! JVS(364) = Jac_FULL(48,55)
  JVS(364) = B(235)+0.2*B(237)+0.14*B(240)+0.1*B(242)+0.6*B(244)+B(246)+0.725*B(248)+B(279)
! JVS(365) = Jac_FULL(48,56)
  JVS(365) = 0.75*B(128)+0.15*B(214)+0.2*B(238)
! JVS(366) = Jac_FULL(48,57)
  JVS(366) = B(148)+0.75*B(212)+B(236)
! JVS(367) = Jac_FULL(48,58)
  JVS(367) = 0.75*B(129)+0.38*B(131)+0.35*B(133)+0.2*B(135)+B(137)+1.45*B(139)+B(149)+B(153)
! JVS(368) = Jac_FULL(48,59)
  JVS(368) = 0
! JVS(369) = Jac_FULL(49,29)
  JVS(369) = B(88)
! JVS(370) = Jac_FULL(49,46)
  JVS(370) = 0
! JVS(371) = Jac_FULL(49,49)
  JVS(371) = -B(138)-B(196)-B(223)-B(247)
! JVS(372) = Jac_FULL(49,51)
  JVS(372) = -B(197)
! JVS(373) = Jac_FULL(49,52)
  JVS(373) = -B(224)
! JVS(374) = Jac_FULL(49,53)
  JVS(374) = B(89)
! JVS(375) = Jac_FULL(49,54)
  JVS(375) = 0
! JVS(376) = Jac_FULL(49,55)
  JVS(376) = -B(248)
! JVS(377) = Jac_FULL(49,58)
  JVS(377) = -B(139)
! JVS(378) = Jac_FULL(49,59)
  JVS(378) = 0
! JVS(379) = Jac_FULL(50,33)
  JVS(379) = B(100)
! JVS(380) = Jac_FULL(50,41)
  JVS(380) = 0
! JVS(381) = Jac_FULL(50,44)
  JVS(381) = 0
! JVS(382) = Jac_FULL(50,46)
  JVS(382) = 0
! JVS(383) = Jac_FULL(50,47)
  JVS(383) = 0
! JVS(384) = Jac_FULL(50,49)
  JVS(384) = 0
! JVS(385) = Jac_FULL(50,50)
  JVS(385) = -B(150)-B(198)-B(225)-B(249)
! JVS(386) = Jac_FULL(50,51)
  JVS(386) = -B(199)
! JVS(387) = Jac_FULL(50,52)
  JVS(387) = -B(226)
! JVS(388) = Jac_FULL(50,53)
  JVS(388) = B(101)
! JVS(389) = Jac_FULL(50,54)
  JVS(389) = 0
! JVS(390) = Jac_FULL(50,55)
  JVS(390) = -B(250)
! JVS(391) = Jac_FULL(50,56)
  JVS(391) = 0
! JVS(392) = Jac_FULL(50,58)
  JVS(392) = -B(151)
! JVS(393) = Jac_FULL(50,59)
  JVS(393) = 0
! JVS(394) = Jac_FULL(51,5)
  JVS(394) = B(70)
! JVS(395) = Jac_FULL(51,10)
  JVS(395) = 0.25*B(90)
! JVS(396) = Jac_FULL(51,11)
  JVS(396) = 0.17*B(92)
! JVS(397) = Jac_FULL(51,14)
  JVS(397) = B(44)
! JVS(398) = Jac_FULL(51,16)
  JVS(398) = 0.17*B(78)
! JVS(399) = Jac_FULL(51,20)
  JVS(399) = 0.65*B(6)+B(40)
! JVS(400) = Jac_FULL(51,21)
  JVS(400) = B(13)
! JVS(401) = Jac_FULL(51,22)
  JVS(401) = 0.1*B(94)
! JVS(402) = Jac_FULL(51,24)
  JVS(402) = 0.12*B(174)
! JVS(403) = Jac_FULL(51,26)
  JVS(403) = B(72)
! JVS(404) = Jac_FULL(51,27)
  JVS(404) = 0.23*B(180)
! JVS(405) = Jac_FULL(51,28)
  JVS(405) = 0.23*B(176)
! JVS(406) = Jac_FULL(51,29)
  JVS(406) = 0.26*B(178)
! JVS(407) = Jac_FULL(51,30)
  JVS(407) = B(20)
! JVS(408) = Jac_FULL(51,31)
  JVS(408) = 0.8*B(18)+B(102)+B(158)
! JVS(409) = Jac_FULL(51,32)
  JVS(409) = -B(269)+B(271)
! JVS(410) = Jac_FULL(51,34)
  JVS(410) = B(19)
! JVS(411) = Jac_FULL(51,35)
  JVS(411) = B(144)-B(202)+2*B(229)+B(252)
! JVS(412) = Jac_FULL(51,36)
  JVS(412) = B(146)-B(204)+2*B(231)+B(254)
! JVS(413) = Jac_FULL(51,37)
  JVS(413) = B(136)-B(194)+B(221)+0.5*B(245)
! JVS(414) = Jac_FULL(51,38)
  JVS(414) = -B(208)+0.5*B(276)
! JVS(415) = Jac_FULL(51,39)
  JVS(415) = -B(258)+B(260)
! JVS(416) = Jac_FULL(51,40)
  JVS(416) = B(134)-B(192)+B(219)+0.5*B(243)
! JVS(417) = Jac_FULL(51,41)
  JVS(417) = 0.92*B(130)-B(188)+B(215)+0.5*B(239)
! JVS(418) = Jac_FULL(51,42)
  JVS(418) = B(14)
! JVS(419) = Jac_FULL(51,43)
  JVS(419) = 2*B(11)+B(96)+B(154)
! JVS(420) = Jac_FULL(51,44)
  JVS(420) = 0.76*B(132)-B(190)+B(217)+0.5*B(241)
! JVS(421) = Jac_FULL(51,45)
  JVS(421) = 0.92*B(142)-B(206)+0.46*B(233)+0.92*B(256)
! JVS(422) = Jac_FULL(51,46)
  JVS(422) = B(32)-B(34)+0.12*B(175)+0.23*B(177)+0.26*B(179)+0.23*B(181)
! JVS(423) = Jac_FULL(51,47)
  JVS(423) = B(21)
! JVS(424) = Jac_FULL(51,48)
  JVS(424) = B(12)
! JVS(425) = Jac_FULL(51,49)
  JVS(425) = B(138)-B(196)+B(223)+0.5*B(247)
! JVS(426) = Jac_FULL(51,50)
  JVS(426) = B(150)-B(198)+B(225)+0.5*B(249)
! JVS(427) = Jac_FULL(51,51)
  JVS(427) = -B(35)-B(36)-B(38)-2*B(41)-2*B(42)-B(56)-B(68)-B(182)-B(184)-B(186)-B(189)-B(191)-B(193)-B(195)-B(197)&
               &-B(199)-B(200)-B(203)-B(205)-B(207)-B(209)-B(259)-B(270)
! JVS(428) = Jac_FULL(51,52)
  JVS(428) = B(126)-B(183)+B(210)+B(211)+B(213)+B(216)+B(218)+B(220)+B(222)+B(224)+B(226)+0.5*B(227)+2*B(230)+2*B(232)&
               &+0.46*B(234)+B(261)+B(272)+0.5*B(277)
! JVS(429) = Jac_FULL(51,53)
  JVS(429) = B(33)+B(45)-B(69)+B(71)+B(73)+0.17*B(79)+0.25*B(91)+0.17*B(93)+0.1*B(95)+B(97)+B(103)
! JVS(430) = Jac_FULL(51,54)
  JVS(430) = -B(57)+B(155)+B(159)
! JVS(431) = Jac_FULL(51,55)
  JVS(431) = -B(201)+0.5*B(228)+0.5*B(235)+0.5*B(237)+0.5*B(240)+0.5*B(242)+0.5*B(244)+0.5*B(246)+0.5*B(248)+0.5*B(250)&
               &+B(253)+B(255)+0.92*B(257)
! JVS(432) = Jac_FULL(51,56)
  JVS(432) = 0.964*B(128)-B(187)+B(214)+0.5*B(238)
! JVS(433) = Jac_FULL(51,57)
  JVS(433) = B(148)-B(185)+B(212)+0.5*B(236)
! JVS(434) = Jac_FULL(51,58)
  JVS(434) = -B(37)+B(127)+0.964*B(129)+0.92*B(131)+0.76*B(133)+B(135)+B(137)+B(139)+0.92*B(143)+B(145)+B(147)+B(149)&
               &+B(151)
! JVS(435) = Jac_FULL(51,59)
  JVS(435) = -B(39)
! JVS(436) = Jac_FULL(52,17)
  JVS(436) = B(74)
! JVS(437) = Jac_FULL(52,18)
  JVS(437) = B(15)
! JVS(438) = Jac_FULL(52,21)
  JVS(438) = 0.5*B(108)
! JVS(439) = Jac_FULL(52,27)
  JVS(439) = 0.22*B(180)
! JVS(440) = Jac_FULL(52,28)
  JVS(440) = 0.22*B(176)
! JVS(441) = Jac_FULL(52,29)
  JVS(441) = 0.31*B(178)
! JVS(442) = Jac_FULL(52,32)
  JVS(442) = -B(271)+B(273)
! JVS(443) = Jac_FULL(52,35)
  JVS(443) = -B(229)+B(252)
! JVS(444) = Jac_FULL(52,36)
  JVS(444) = -B(231)+B(254)
! JVS(445) = Jac_FULL(52,37)
  JVS(445) = -B(221)+0.5*B(245)
! JVS(446) = Jac_FULL(52,38)
  JVS(446) = -B(276)+0.5*B(278)
! JVS(447) = Jac_FULL(52,39)
  JVS(447) = -B(260)+B(262)
! JVS(448) = Jac_FULL(52,40)
  JVS(448) = -B(219)+0.5*B(243)
! JVS(449) = Jac_FULL(52,41)
  JVS(449) = -B(215)+0.5*B(239)
! JVS(450) = Jac_FULL(52,44)
  JVS(450) = -B(217)+0.5*B(241)
! JVS(451) = Jac_FULL(52,45)
  JVS(451) = -B(233)+B(256)
! JVS(452) = Jac_FULL(52,46)
  JVS(452) = 0.22*B(177)+0.31*B(179)+0.22*B(181)
! JVS(453) = Jac_FULL(52,48)
  JVS(453) = B(12)
! JVS(454) = Jac_FULL(52,49)
  JVS(454) = -B(223)+0.5*B(247)
! JVS(455) = Jac_FULL(52,50)
  JVS(455) = -B(225)+0.5*B(249)
! JVS(456) = Jac_FULL(52,51)
  JVS(456) = -B(182)
! JVS(457) = Jac_FULL(52,52)
  JVS(457) = -B(126)-B(183)-2*B(210)-B(211)-B(213)-B(216)-B(218)-B(220)-B(222)-B(224)-B(226)-0.5*B(227)-B(230)-B(232)&
               &-B(234)-B(261)-B(272)-B(277)
! JVS(458) = Jac_FULL(52,53)
  JVS(458) = B(75)+0.5*B(109)
! JVS(459) = Jac_FULL(52,54)
  JVS(459) = 0
! JVS(460) = Jac_FULL(52,55)
  JVS(460) = B(140)-0.5*B(228)+0.5*B(235)+0.5*B(237)+0.5*B(240)+0.5*B(242)+0.5*B(244)+0.5*B(246)+0.5*B(248)+0.5*B(250)+2&
               &*B(251)+B(253)+B(255)+B(257)+B(263)+B(274)+0.5*B(279)
! JVS(461) = Jac_FULL(52,56)
  JVS(461) = -B(214)+0.5*B(238)
! JVS(462) = Jac_FULL(52,57)
  JVS(462) = -B(212)+0.5*B(236)
! JVS(463) = Jac_FULL(52,58)
  JVS(463) = -B(127)+B(141)
! JVS(464) = Jac_FULL(52,59)
  JVS(464) = 0
! JVS(465) = Jac_FULL(53,5)
  JVS(465) = -B(70)
! JVS(466) = Jac_FULL(53,6)
  JVS(466) = -B(76)
! JVS(467) = Jac_FULL(53,7)
  JVS(467) = 2*B(28)
! JVS(468) = Jac_FULL(53,8)
  JVS(468) = -B(80)
! JVS(469) = Jac_FULL(53,9)
  JVS(469) = -B(82)
! JVS(470) = Jac_FULL(53,10)
  JVS(470) = -B(90)
! JVS(471) = Jac_FULL(53,11)
  JVS(471) = -B(92)
! JVS(472) = Jac_FULL(53,13)
  JVS(472) = B(4)
! JVS(473) = Jac_FULL(53,14)
  JVS(473) = 2*B(9)-B(44)
! JVS(474) = Jac_FULL(53,16)
  JVS(474) = -B(78)
! JVS(475) = Jac_FULL(53,17)
  JVS(475) = -B(74)
! JVS(476) = Jac_FULL(53,18)
  JVS(476) = B(15)-B(112)
! JVS(477) = Jac_FULL(53,20)
  JVS(477) = 0.35*B(6)-B(66)
! JVS(478) = Jac_FULL(53,21)
  JVS(478) = B(13)-0.5*B(108)
! JVS(479) = Jac_FULL(53,22)
  JVS(479) = -1.9*B(94)
! JVS(480) = Jac_FULL(53,23)
  JVS(480) = -B(114)
! JVS(481) = Jac_FULL(53,24)
  JVS(481) = -B(84)
! JVS(482) = Jac_FULL(53,25)
  JVS(482) = B(5)-B(64)
! JVS(483) = Jac_FULL(53,26)
  JVS(483) = -B(72)
! JVS(484) = Jac_FULL(53,27)
  JVS(484) = -B(118)+0.1*B(180)
! JVS(485) = Jac_FULL(53,28)
  JVS(485) = -B(86)+0.1*B(176)
! JVS(486) = Jac_FULL(53,29)
  JVS(486) = -B(88)+0.14*B(178)
! JVS(487) = Jac_FULL(53,30)
  JVS(487) = -B(106)
! JVS(488) = Jac_FULL(53,31)
  JVS(488) = -B(102)
! JVS(489) = Jac_FULL(53,33)
  JVS(489) = -B(100)
! JVS(490) = Jac_FULL(53,34)
  JVS(490) = -B(104)
! JVS(491) = Jac_FULL(53,35)
  JVS(491) = 0
! JVS(492) = Jac_FULL(53,36)
  JVS(492) = 0
! JVS(493) = Jac_FULL(53,41)
  JVS(493) = 0
! JVS(494) = Jac_FULL(53,42)
  JVS(494) = B(14)-0.5*B(110)
! JVS(495) = Jac_FULL(53,43)
  JVS(495) = -B(96)
! JVS(496) = Jac_FULL(53,44)
  JVS(496) = 0
! JVS(497) = Jac_FULL(53,45)
  JVS(497) = 0
! JVS(498) = Jac_FULL(53,46)
  JVS(498) = -B(32)+B(34)+0.1*B(177)+0.14*B(179)+0.1*B(181)
! JVS(499) = Jac_FULL(53,47)
  JVS(499) = -B(116)
! JVS(500) = Jac_FULL(53,48)
  JVS(500) = -B(98)
! JVS(501) = Jac_FULL(53,49)
  JVS(501) = 0
! JVS(502) = Jac_FULL(53,50)
  JVS(502) = 0
! JVS(503) = Jac_FULL(53,51)
  JVS(503) = B(35)+B(36)-B(68)
! JVS(504) = Jac_FULL(53,52)
  JVS(504) = 0
! JVS(505) = Jac_FULL(53,53)
  JVS(505) = -B(33)-B(45)-B(46)-B(62)-B(65)-B(67)-B(69)-B(71)-B(73)-B(75)-B(77)-B(79)-B(81)-B(83)-B(85)-B(87)-B(89)&
               &-B(91)-B(93)-1.9*B(95)-B(97)-B(99)-B(101)-B(103)-B(105)-B(107)-0.5*B(109)-0.5*B(111)-B(113)-B(115)-B(117)&
               &-B(119)
! JVS(506) = Jac_FULL(53,54)
  JVS(506) = 0
! JVS(507) = Jac_FULL(53,55)
  JVS(507) = 0
! JVS(508) = Jac_FULL(53,56)
  JVS(508) = 0
! JVS(509) = Jac_FULL(53,57)
  JVS(509) = 0
! JVS(510) = Jac_FULL(53,58)
  JVS(510) = B(37)-B(47)
! JVS(511) = Jac_FULL(53,59)
  JVS(511) = -B(63)
! JVS(512) = Jac_FULL(54,15)
  JVS(512) = B(60)
! JVS(513) = Jac_FULL(54,20)
  JVS(513) = 0.35*B(6)
! JVS(514) = Jac_FULL(54,22)
  JVS(514) = -B(164)
! JVS(515) = Jac_FULL(54,23)
  JVS(515) = B(114)
! JVS(516) = Jac_FULL(54,24)
  JVS(516) = -B(166)
! JVS(517) = Jac_FULL(54,25)
  JVS(517) = B(64)
! JVS(518) = Jac_FULL(54,27)
  JVS(518) = -B(172)
! JVS(519) = Jac_FULL(54,28)
  JVS(519) = -B(168)
! JVS(520) = Jac_FULL(54,29)
  JVS(520) = -B(170)
! JVS(521) = Jac_FULL(54,30)
  JVS(521) = -B(162)
! JVS(522) = Jac_FULL(54,31)
  JVS(522) = -B(158)
! JVS(523) = Jac_FULL(54,34)
  JVS(523) = -B(160)
! JVS(524) = Jac_FULL(54,35)
  JVS(524) = 0
! JVS(525) = Jac_FULL(54,36)
  JVS(525) = 0
! JVS(526) = Jac_FULL(54,43)
  JVS(526) = -B(154)
! JVS(527) = Jac_FULL(54,44)
  JVS(527) = 0
! JVS(528) = Jac_FULL(54,45)
  JVS(528) = 0
! JVS(529) = Jac_FULL(54,46)
  JVS(529) = B(50)
! JVS(530) = Jac_FULL(54,48)
  JVS(530) = -B(156)
! JVS(531) = Jac_FULL(54,49)
  JVS(531) = 0
! JVS(532) = Jac_FULL(54,50)
  JVS(532) = 0
! JVS(533) = Jac_FULL(54,51)
  JVS(533) = -B(56)
! JVS(534) = Jac_FULL(54,52)
  JVS(534) = 0
! JVS(535) = Jac_FULL(54,53)
  JVS(535) = B(65)+B(115)
! JVS(536) = Jac_FULL(54,54)
  JVS(536) = -B(7)-B(8)-B(52)-B(54)-B(57)-B(58)-B(155)-B(157)-B(159)-B(161)-B(163)-B(165)-B(167)-B(169)-B(171)-B(173)
! JVS(537) = Jac_FULL(54,55)
  JVS(537) = 0
! JVS(538) = Jac_FULL(54,56)
  JVS(538) = 0
! JVS(539) = Jac_FULL(54,57)
  JVS(539) = 0
! JVS(540) = Jac_FULL(54,58)
  JVS(540) = -B(53)
! JVS(541) = Jac_FULL(54,59)
  JVS(541) = B(51)-B(55)-B(59)
! JVS(542) = Jac_FULL(55,18)
  JVS(542) = B(112)
! JVS(543) = Jac_FULL(55,23)
  JVS(543) = B(122)
! JVS(544) = Jac_FULL(55,32)
  JVS(544) = -B(273)
! JVS(545) = Jac_FULL(55,33)
  JVS(545) = B(16)
! JVS(546) = Jac_FULL(55,34)
  JVS(546) = B(19)+B(104)+B(160)
! JVS(547) = Jac_FULL(55,35)
  JVS(547) = -B(252)
! JVS(548) = Jac_FULL(55,36)
  JVS(548) = -B(254)
! JVS(549) = Jac_FULL(55,37)
  JVS(549) = -B(245)
! JVS(550) = Jac_FULL(55,38)
  JVS(550) = -B(278)
! JVS(551) = Jac_FULL(55,39)
  JVS(551) = -B(262)
! JVS(552) = Jac_FULL(55,40)
  JVS(552) = -B(243)
! JVS(553) = Jac_FULL(55,41)
  JVS(553) = -B(239)
! JVS(554) = Jac_FULL(55,44)
  JVS(554) = -B(241)
! JVS(555) = Jac_FULL(55,45)
  JVS(555) = 0.05*B(142)+0.025*B(233)-0.95*B(256)
! JVS(556) = Jac_FULL(55,46)
  JVS(556) = 0
! JVS(557) = Jac_FULL(55,47)
  JVS(557) = 0
! JVS(558) = Jac_FULL(55,48)
  JVS(558) = B(98)+B(156)
! JVS(559) = Jac_FULL(55,49)
  JVS(559) = -B(247)
! JVS(560) = Jac_FULL(55,50)
  JVS(560) = -B(249)
! JVS(561) = Jac_FULL(55,51)
  JVS(561) = -B(200)
! JVS(562) = Jac_FULL(55,52)
  JVS(562) = -B(227)+0.025*B(234)
! JVS(563) = Jac_FULL(55,53)
  JVS(563) = B(99)+B(105)+B(113)
! JVS(564) = Jac_FULL(55,54)
  JVS(564) = B(157)+B(161)
! JVS(565) = Jac_FULL(55,55)
  JVS(565) = -B(120)-B(140)-B(201)-B(228)-B(235)-B(237)-B(240)-B(242)-B(244)-B(246)-B(248)-B(250)-2*B(251)-B(253)-B(255)&
               &-0.95*B(257)-B(263)-B(274)-B(279)
! JVS(566) = Jac_FULL(55,56)
  JVS(566) = -B(238)
! JVS(567) = Jac_FULL(55,57)
  JVS(567) = -B(236)
! JVS(568) = Jac_FULL(55,58)
  JVS(568) = -B(141)+0.05*B(143)
! JVS(569) = Jac_FULL(55,59)
  JVS(569) = -B(121)
! JVS(570) = Jac_FULL(56,16)
  JVS(570) = 0.83*B(78)
! JVS(571) = Jac_FULL(56,42)
  JVS(571) = 0.5*B(110)
! JVS(572) = Jac_FULL(56,44)
  JVS(572) = 0
! JVS(573) = Jac_FULL(56,45)
  JVS(573) = 0
! JVS(574) = Jac_FULL(56,46)
  JVS(574) = 0
! JVS(575) = Jac_FULL(56,47)
  JVS(575) = B(116)
! JVS(576) = Jac_FULL(56,49)
  JVS(576) = 0
! JVS(577) = Jac_FULL(56,50)
  JVS(577) = 0
! JVS(578) = Jac_FULL(56,51)
  JVS(578) = -B(186)
! JVS(579) = Jac_FULL(56,52)
  JVS(579) = -B(213)
! JVS(580) = Jac_FULL(56,53)
  JVS(580) = 0.83*B(79)+0.5*B(111)+B(117)
! JVS(581) = Jac_FULL(56,54)
  JVS(581) = 0
! JVS(582) = Jac_FULL(56,55)
  JVS(582) = -B(237)
! JVS(583) = Jac_FULL(56,56)
  JVS(583) = -B(128)-B(187)-B(214)-B(238)
! JVS(584) = Jac_FULL(56,57)
  JVS(584) = 0
! JVS(585) = Jac_FULL(56,58)
  JVS(585) = -B(129)
! JVS(586) = Jac_FULL(56,59)
  JVS(586) = 0
! JVS(587) = Jac_FULL(57,6)
  JVS(587) = B(76)
! JVS(588) = Jac_FULL(57,33)
  JVS(588) = B(16)
! JVS(589) = Jac_FULL(57,41)
  JVS(589) = 0
! JVS(590) = Jac_FULL(57,44)
  JVS(590) = 0
! JVS(591) = Jac_FULL(57,46)
  JVS(591) = 0
! JVS(592) = Jac_FULL(57,47)
  JVS(592) = 0
! JVS(593) = Jac_FULL(57,49)
  JVS(593) = 0
! JVS(594) = Jac_FULL(57,51)
  JVS(594) = -B(184)
! JVS(595) = Jac_FULL(57,52)
  JVS(595) = -B(211)
! JVS(596) = Jac_FULL(57,53)
  JVS(596) = B(77)
! JVS(597) = Jac_FULL(57,54)
  JVS(597) = 0
! JVS(598) = Jac_FULL(57,55)
  JVS(598) = -B(235)
! JVS(599) = Jac_FULL(57,56)
  JVS(599) = 0
! JVS(600) = Jac_FULL(57,57)
  JVS(600) = -B(148)-B(185)-B(212)-B(236)
! JVS(601) = Jac_FULL(57,58)
  JVS(601) = -B(149)
! JVS(602) = Jac_FULL(57,59)
  JVS(602) = 0
! JVS(603) = Jac_FULL(58,13)
  JVS(603) = B(4)
! JVS(604) = Jac_FULL(58,19)
  JVS(604) = B(24)
! JVS(605) = Jac_FULL(58,35)
  JVS(605) = -B(144)
! JVS(606) = Jac_FULL(58,36)
  JVS(606) = -B(146)
! JVS(607) = Jac_FULL(58,37)
  JVS(607) = -B(136)
! JVS(608) = Jac_FULL(58,38)
  JVS(608) = -B(152)
! JVS(609) = Jac_FULL(58,39)
  JVS(609) = -B(265)
! JVS(610) = Jac_FULL(58,40)
  JVS(610) = -B(134)
! JVS(611) = Jac_FULL(58,41)
  JVS(611) = -B(130)
! JVS(612) = Jac_FULL(58,44)
  JVS(612) = -B(132)
! JVS(613) = Jac_FULL(58,45)
  JVS(613) = -B(142)
! JVS(614) = Jac_FULL(58,46)
  JVS(614) = -B(30)
! JVS(615) = Jac_FULL(58,49)
  JVS(615) = -B(138)
! JVS(616) = Jac_FULL(58,50)
  JVS(616) = -B(150)
! JVS(617) = Jac_FULL(58,51)
  JVS(617) = -B(36)
! JVS(618) = Jac_FULL(58,52)
  JVS(618) = -B(126)
! JVS(619) = Jac_FULL(58,53)
  JVS(619) = -B(46)
! JVS(620) = Jac_FULL(58,54)
  JVS(620) = B(7)-B(52)+B(54)
! JVS(621) = Jac_FULL(58,55)
  JVS(621) = -B(140)
! JVS(622) = Jac_FULL(58,56)
  JVS(622) = -B(128)
! JVS(623) = Jac_FULL(58,57)
  JVS(623) = -B(148)
! JVS(624) = Jac_FULL(58,58)
  JVS(624) = -B(31)-B(37)-B(47)-2*B(48)-B(53)-B(127)-B(129)-B(131)-B(133)-B(135)-B(137)-B(139)-B(141)-B(143)-B(145)&
               &-B(147)-B(149)-B(151)-B(153)-B(266)
! JVS(625) = Jac_FULL(58,59)
  JVS(625) = B(1)+B(25)+B(55)
! JVS(626) = Jac_FULL(59,12)
  JVS(626) = B(125)
! JVS(627) = Jac_FULL(59,15)
  JVS(627) = B(60)
! JVS(628) = Jac_FULL(59,19)
  JVS(628) = -B(24)
! JVS(629) = Jac_FULL(59,20)
  JVS(629) = 0.65*B(6)+B(40)+B(66)
! JVS(630) = Jac_FULL(59,23)
  JVS(630) = B(122)
! JVS(631) = Jac_FULL(59,25)
  JVS(631) = B(5)
! JVS(632) = Jac_FULL(59,30)
  JVS(632) = 0
! JVS(633) = Jac_FULL(59,31)
  JVS(633) = 0
! JVS(634) = Jac_FULL(59,32)
  JVS(634) = -B(267)
! JVS(635) = Jac_FULL(59,34)
  JVS(635) = 0
! JVS(636) = Jac_FULL(59,35)
  JVS(636) = B(144)
! JVS(637) = Jac_FULL(59,36)
  JVS(637) = B(146)
! JVS(638) = Jac_FULL(59,37)
  JVS(638) = B(136)
! JVS(639) = Jac_FULL(59,38)
  JVS(639) = 2*B(152)+B(276)+B(278)+2*B(280)
! JVS(640) = Jac_FULL(59,39)
  JVS(640) = B(265)
! JVS(641) = Jac_FULL(59,40)
  JVS(641) = B(134)
! JVS(642) = Jac_FULL(59,41)
  JVS(642) = 0.92*B(130)
! JVS(643) = Jac_FULL(59,43)
  JVS(643) = 0
! JVS(644) = Jac_FULL(59,44)
  JVS(644) = 0.76*B(132)
! JVS(645) = Jac_FULL(59,45)
  JVS(645) = -B(123)+B(142)
! JVS(646) = Jac_FULL(59,46)
  JVS(646) = B(30)-B(50)
! JVS(647) = Jac_FULL(59,47)
  JVS(647) = B(21)+B(116)
! JVS(648) = Jac_FULL(59,48)
  JVS(648) = 0
! JVS(649) = Jac_FULL(59,49)
  JVS(649) = B(138)
! JVS(650) = Jac_FULL(59,50)
  JVS(650) = B(150)
! JVS(651) = Jac_FULL(59,51)
  JVS(651) = B(36)-B(38)
! JVS(652) = Jac_FULL(59,52)
  JVS(652) = B(126)+B(277)
! JVS(653) = Jac_FULL(59,53)
  JVS(653) = -B(62)+B(67)+B(117)
! JVS(654) = Jac_FULL(59,54)
  JVS(654) = B(8)+2*B(52)-B(58)
! JVS(655) = Jac_FULL(59,55)
  JVS(655) = -B(120)+B(140)+B(279)
! JVS(656) = Jac_FULL(59,56)
  JVS(656) = 0.964*B(128)
! JVS(657) = Jac_FULL(59,57)
  JVS(657) = B(148)
! JVS(658) = Jac_FULL(59,58)
  JVS(658) = B(31)+B(37)+2*B(48)+2*B(53)+B(127)+0.964*B(129)+0.92*B(131)+0.76*B(133)+B(135)+B(137)+B(139)+B(141)+B(143)&
               &+B(145)+B(147)+B(149)+B(151)+2*B(153)+B(266)
! JVS(659) = Jac_FULL(59,59)
  JVS(659) = -B(1)-B(25)-B(39)-B(51)-B(59)-B(63)-B(121)-B(124)-B(268)

END SUBROUTINE radm2sorg_Jac_SP

! End of radm2sorg_Jac_SP function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! SPARSE_UTIL - SPARSE utility functions
!   Arguments :
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE radm2sorg_KppDecomp( JVS, IER )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Sparse LU factorization
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

!  USE radm2sorg_Parameters
!  USE radm2sorg_JacobianSP

      INTEGER  :: IER
      REAL(kind=dp) :: JVS(659), W(59), a
      INTEGER  :: k, kk, j, jj

      a = 0. ! mz_rs_20050606
      IER = 0
      DO k=1,NVAR
        ! mz_rs_20050606: don't check if real value == 0
        ! IF ( JVS( LU_DIAG(k) ) .EQ. 0. ) THEN
        IF ( ABS(JVS(LU_DIAG(k))) < TINY(a) ) THEN
            IER = k
            RETURN
        END IF
        DO kk = LU_CROW(k), LU_CROW(k+1)-1
              W( LU_ICOL(kk) ) = JVS(kk)
        END DO
        DO kk = LU_CROW(k), LU_DIAG(k)-1
            j = LU_ICOL(kk)
            a = -W(j) / JVS( LU_DIAG(j) )
            W(j) = -a
            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
               W( LU_ICOL(jj) ) = W( LU_ICOL(jj) ) + a*JVS(jj)
            END DO
         END DO
         DO kk = LU_CROW(k), LU_CROW(k+1)-1
            JVS(kk) = W( LU_ICOL(kk) )
         END DO
      END DO

END SUBROUTINE radm2sorg_KppDecomp

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE radm2sorg_KppDecompCmplx( JVS, IER )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Sparse LU factorization, complex
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

!  USE radm2sorg_Parameters
!  USE radm2sorg_JacobianSP

      INTEGER  :: IER
      DOUBLE COMPLEX :: JVS(659), W(59), a
      INTEGER  :: k, kk, j, jj

      IER = 0
      DO k=1,NVAR
        IF ( JVS( LU_DIAG(k) ) .EQ. 0. ) THEN
            IER = k
            RETURN
        END IF
        DO kk = LU_CROW(k), LU_CROW(k+1)-1
              W( LU_ICOL(kk) ) = JVS(kk)
        END DO
        DO kk = LU_CROW(k), LU_DIAG(k)-1
            j = LU_ICOL(kk)
            a = -W(j) / JVS( LU_DIAG(j) )
            W(j) = -a
            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
               W( LU_ICOL(jj) ) = W( LU_ICOL(jj) ) + a*JVS(jj)
            END DO
         END DO
         DO kk = LU_CROW(k), LU_CROW(k+1)-1
            JVS(kk) = W( LU_ICOL(kk) )
         END DO
      END DO

END SUBROUTINE radm2sorg_KppDecompCmplx

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE radm2sorg_KppSolveIndirect( JVS, X )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Sparse solve subroutine using indirect addressing
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

!  USE radm2sorg_Parameters
!  USE radm2sorg_JacobianSP

      INTEGER i, j
      REAL(kind=dp) JVS(659), X(59), sum

      DO i=1,NVAR
         DO j = LU_CROW(i), LU_DIAG(i)-1
             X(i) = X(i) - JVS(j)*X(LU_ICOL(j));
         END DO
      END DO

      DO i=NVAR,1,-1
        sum = X(i);
        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
          sum = sum - JVS(j)*X(LU_ICOL(j));
        END DO
        X(i) = sum/JVS(LU_DIAG(i));
      END DO

END SUBROUTINE radm2sorg_KppSolveIndirect

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
SUBROUTINE radm2sorg_KppSolveCmplx( JVS, X )
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!        Complex sparse solve subroutine using indirect addressing
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

!  USE radm2sorg_Parameters
!  USE radm2sorg_JacobianSP

      INTEGER i, j
      DOUBLE COMPLEX JVS(659), X(59), sum

      DO i=1,NVAR
         DO j = LU_CROW(i), LU_DIAG(i)-1
             X(i) = X(i) - JVS(j)*X(LU_ICOL(j));
         END DO
      END DO

      DO i=NVAR,1,-1
        sum = X(i);
        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
          sum = sum - JVS(j)*X(LU_ICOL(j));
        END DO
        X(i) = sum/JVS(LU_DIAG(i));
      END DO

END SUBROUTINE radm2sorg_KppSolveCmplx
! End of SPARSE_UTIL function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! radm2sorg_KppSolve - sparse back substitution
!   Arguments :
!      JVS       - sparse Jacobian of variables
!      X         - Vector for variables
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE radm2sorg_KppSolve ( JVS, X )

! JVS - sparse Jacobian of variables
  REAL(kind=dp) :: JVS(LU_NONZERO)
! X - Vector for variables
  REAL(kind=dp) :: X(NVAR)

  X(19) = X(19)-JVS(67)*X(7)
  X(22) = X(22)-JVS(80)*X(10)-JVS(81)*X(11)
  X(25) = X(25)-JVS(93)*X(15)-JVS(94)*X(22)
  X(26) = X(26)-JVS(105)*X(24)
  X(32) = X(32)-JVS(149)*X(22)
  X(33) = X(33)-JVS(157)*X(16)-JVS(158)*X(29)
  X(35) = X(35)-JVS(181)*X(10)
  X(36) = X(36)-JVS(188)*X(11)
  X(37) = X(37)-JVS(195)*X(27)-JVS(196)*X(28)
  X(38) = X(38)-JVS(205)*X(24)-JVS(206)*X(27)-JVS(207)*X(28)-JVS(208)*X(29)
  X(39) = X(39)-JVS(217)*X(8)-JVS(218)*X(9)-JVS(219)*X(22)-JVS(220)*X(23)
  X(40) = X(40)-JVS(230)*X(24)
  X(41) = X(41)-JVS(239)*X(8)
  X(42) = X(42)-JVS(246)*X(32)-JVS(247)*X(35)-JVS(248)*X(36)-JVS(249)*X(37)-JVS(250)*X(39)-JVS(251)*X(40)-JVS(252)*X(41)
  X(43) = X(43)-JVS(268)*X(16)-JVS(269)*X(21)-JVS(270)*X(23)-JVS(271)*X(24)-JVS(272)*X(27)-JVS(273)*X(28)-JVS(274)*X(29)&
            &-JVS(275)*X(31)-JVS(276)*X(32)-JVS(277)*X(35)-JVS(278)*X(36)-JVS(279)*X(37)-JVS(280)*X(38)-JVS(281)*X(39)&
            &-JVS(282)*X(40)-JVS(283)*X(41)
  X(44) = X(44)-JVS(299)*X(9)
  X(45) = X(45)-JVS(306)*X(12)-JVS(307)*X(22)-JVS(308)*X(30)-JVS(309)*X(35)-JVS(310)*X(36)
  X(46) = X(46)-JVS(319)*X(19)-JVS(320)*X(24)-JVS(321)*X(27)-JVS(322)*X(28)-JVS(323)*X(29)
  X(47) = X(47)-JVS(330)*X(32)-JVS(331)*X(38)-JVS(332)*X(41)-JVS(333)*X(44)-JVS(334)*X(46)
  X(48) = X(48)-JVS(344)*X(16)-JVS(345)*X(27)-JVS(346)*X(28)-JVS(347)*X(29)-JVS(348)*X(37)-JVS(349)*X(38)-JVS(350)*X(40)&
            &-JVS(351)*X(41)-JVS(352)*X(42)-JVS(353)*X(44)-JVS(354)*X(45)-JVS(355)*X(46)-JVS(356)*X(47)
  X(49) = X(49)-JVS(369)*X(29)-JVS(370)*X(46)
  X(50) = X(50)-JVS(379)*X(33)-JVS(380)*X(41)-JVS(381)*X(44)-JVS(382)*X(46)-JVS(383)*X(47)-JVS(384)*X(49)
  X(51) = X(51)-JVS(394)*X(5)-JVS(395)*X(10)-JVS(396)*X(11)-JVS(397)*X(14)-JVS(398)*X(16)-JVS(399)*X(20)-JVS(400)*X(21)&
            &-JVS(401)*X(22)-JVS(402)*X(24)-JVS(403)*X(26)-JVS(404)*X(27)-JVS(405)*X(28)-JVS(406)*X(29)-JVS(407)*X(30)&
            &-JVS(408)*X(31)-JVS(409)*X(32)-JVS(410)*X(34)-JVS(411)*X(35)-JVS(412)*X(36)-JVS(413)*X(37)-JVS(414)*X(38)&
            &-JVS(415)*X(39)-JVS(416)*X(40)-JVS(417)*X(41)-JVS(418)*X(42)-JVS(419)*X(43)-JVS(420)*X(44)-JVS(421)*X(45)&
            &-JVS(422)*X(46)-JVS(423)*X(47)-JVS(424)*X(48)-JVS(425)*X(49)-JVS(426)*X(50)
  X(52) = X(52)-JVS(436)*X(17)-JVS(437)*X(18)-JVS(438)*X(21)-JVS(439)*X(27)-JVS(440)*X(28)-JVS(441)*X(29)-JVS(442)*X(32)&
            &-JVS(443)*X(35)-JVS(444)*X(36)-JVS(445)*X(37)-JVS(446)*X(38)-JVS(447)*X(39)-JVS(448)*X(40)-JVS(449)*X(41)&
            &-JVS(450)*X(44)-JVS(451)*X(45)-JVS(452)*X(46)-JVS(453)*X(48)-JVS(454)*X(49)-JVS(455)*X(50)-JVS(456)*X(51)
  X(53) = X(53)-JVS(465)*X(5)-JVS(466)*X(6)-JVS(467)*X(7)-JVS(468)*X(8)-JVS(469)*X(9)-JVS(470)*X(10)-JVS(471)*X(11)&
            &-JVS(472)*X(13)-JVS(473)*X(14)-JVS(474)*X(16)-JVS(475)*X(17)-JVS(476)*X(18)-JVS(477)*X(20)-JVS(478)*X(21)&
            &-JVS(479)*X(22)-JVS(480)*X(23)-JVS(481)*X(24)-JVS(482)*X(25)-JVS(483)*X(26)-JVS(484)*X(27)-JVS(485)*X(28)&
            &-JVS(486)*X(29)-JVS(487)*X(30)-JVS(488)*X(31)-JVS(489)*X(33)-JVS(490)*X(34)-JVS(491)*X(35)-JVS(492)*X(36)&
            &-JVS(493)*X(41)-JVS(494)*X(42)-JVS(495)*X(43)-JVS(496)*X(44)-JVS(497)*X(45)-JVS(498)*X(46)-JVS(499)*X(47)&
            &-JVS(500)*X(48)-JVS(501)*X(49)-JVS(502)*X(50)-JVS(503)*X(51)-JVS(504)*X(52)
  X(54) = X(54)-JVS(512)*X(15)-JVS(513)*X(20)-JVS(514)*X(22)-JVS(515)*X(23)-JVS(516)*X(24)-JVS(517)*X(25)-JVS(518)*X(27)&
            &-JVS(519)*X(28)-JVS(520)*X(29)-JVS(521)*X(30)-JVS(522)*X(31)-JVS(523)*X(34)-JVS(524)*X(35)-JVS(525)*X(36)&
            &-JVS(526)*X(43)-JVS(527)*X(44)-JVS(528)*X(45)-JVS(529)*X(46)-JVS(530)*X(48)-JVS(531)*X(49)-JVS(532)*X(50)&
            &-JVS(533)*X(51)-JVS(534)*X(52)-JVS(535)*X(53)
  X(55) = X(55)-JVS(542)*X(18)-JVS(543)*X(23)-JVS(544)*X(32)-JVS(545)*X(33)-JVS(546)*X(34)-JVS(547)*X(35)-JVS(548)*X(36)&
            &-JVS(549)*X(37)-JVS(550)*X(38)-JVS(551)*X(39)-JVS(552)*X(40)-JVS(553)*X(41)-JVS(554)*X(44)-JVS(555)*X(45)&
            &-JVS(556)*X(46)-JVS(557)*X(47)-JVS(558)*X(48)-JVS(559)*X(49)-JVS(560)*X(50)-JVS(561)*X(51)-JVS(562)*X(52)&
            &-JVS(563)*X(53)-JVS(564)*X(54)
  X(56) = X(56)-JVS(570)*X(16)-JVS(571)*X(42)-JVS(572)*X(44)-JVS(573)*X(45)-JVS(574)*X(46)-JVS(575)*X(47)-JVS(576)*X(49)&
            &-JVS(577)*X(50)-JVS(578)*X(51)-JVS(579)*X(52)-JVS(580)*X(53)-JVS(581)*X(54)-JVS(582)*X(55)
  X(57) = X(57)-JVS(587)*X(6)-JVS(588)*X(33)-JVS(589)*X(41)-JVS(590)*X(44)-JVS(591)*X(46)-JVS(592)*X(47)-JVS(593)*X(49)&
            &-JVS(594)*X(51)-JVS(595)*X(52)-JVS(596)*X(53)-JVS(597)*X(54)-JVS(598)*X(55)-JVS(599)*X(56)
  X(58) = X(58)-JVS(603)*X(13)-JVS(604)*X(19)-JVS(605)*X(35)-JVS(606)*X(36)-JVS(607)*X(37)-JVS(608)*X(38)-JVS(609)*X(39)&
            &-JVS(610)*X(40)-JVS(611)*X(41)-JVS(612)*X(44)-JVS(613)*X(45)-JVS(614)*X(46)-JVS(615)*X(49)-JVS(616)*X(50)&
            &-JVS(617)*X(51)-JVS(618)*X(52)-JVS(619)*X(53)-JVS(620)*X(54)-JVS(621)*X(55)-JVS(622)*X(56)-JVS(623)*X(57)
  X(59) = X(59)-JVS(626)*X(12)-JVS(627)*X(15)-JVS(628)*X(19)-JVS(629)*X(20)-JVS(630)*X(23)-JVS(631)*X(25)-JVS(632)*X(30)&
            &-JVS(633)*X(31)-JVS(634)*X(32)-JVS(635)*X(34)-JVS(636)*X(35)-JVS(637)*X(36)-JVS(638)*X(37)-JVS(639)*X(38)&
            &-JVS(640)*X(39)-JVS(641)*X(40)-JVS(642)*X(41)-JVS(643)*X(43)-JVS(644)*X(44)-JVS(645)*X(45)-JVS(646)*X(46)&
            &-JVS(647)*X(47)-JVS(648)*X(48)-JVS(649)*X(49)-JVS(650)*X(50)-JVS(651)*X(51)-JVS(652)*X(52)-JVS(653)*X(53)&
            &-JVS(654)*X(54)-JVS(655)*X(55)-JVS(656)*X(56)-JVS(657)*X(57)-JVS(658)*X(58)
  X(59) = X(59)/JVS(659)
  X(58) = (X(58)-JVS(625)*X(59))/(JVS(624))
  X(57) = (X(57)-JVS(601)*X(58)-JVS(602)*X(59))/(JVS(600))
  X(56) = (X(56)-JVS(584)*X(57)-JVS(585)*X(58)-JVS(586)*X(59))/(JVS(583))
  X(55) = (X(55)-JVS(566)*X(56)-JVS(567)*X(57)-JVS(568)*X(58)-JVS(569)*X(59))/(JVS(565))
  X(54) = (X(54)-JVS(537)*X(55)-JVS(538)*X(56)-JVS(539)*X(57)-JVS(540)*X(58)-JVS(541)*X(59))/(JVS(536))
  X(53) = (X(53)-JVS(506)*X(54)-JVS(507)*X(55)-JVS(508)*X(56)-JVS(509)*X(57)-JVS(510)*X(58)-JVS(511)*X(59))/(JVS(505))
  X(52) = (X(52)-JVS(458)*X(53)-JVS(459)*X(54)-JVS(460)*X(55)-JVS(461)*X(56)-JVS(462)*X(57)-JVS(463)*X(58)-JVS(464)&
            &*X(59))/(JVS(457))
  X(51) = (X(51)-JVS(428)*X(52)-JVS(429)*X(53)-JVS(430)*X(54)-JVS(431)*X(55)-JVS(432)*X(56)-JVS(433)*X(57)-JVS(434)&
            &*X(58)-JVS(435)*X(59))/(JVS(427))
  X(50) = (X(50)-JVS(386)*X(51)-JVS(387)*X(52)-JVS(388)*X(53)-JVS(389)*X(54)-JVS(390)*X(55)-JVS(391)*X(56)-JVS(392)&
            &*X(58)-JVS(393)*X(59))/(JVS(385))
  X(49) = (X(49)-JVS(372)*X(51)-JVS(373)*X(52)-JVS(374)*X(53)-JVS(375)*X(54)-JVS(376)*X(55)-JVS(377)*X(58)-JVS(378)&
            &*X(59))/(JVS(371))
  X(48) = (X(48)-JVS(358)*X(49)-JVS(359)*X(50)-JVS(360)*X(51)-JVS(361)*X(52)-JVS(362)*X(53)-JVS(363)*X(54)-JVS(364)&
            &*X(55)-JVS(365)*X(56)-JVS(366)*X(57)-JVS(367)*X(58)-JVS(368)*X(59))/(JVS(357))
  X(47) = (X(47)-JVS(336)*X(51)-JVS(337)*X(52)-JVS(338)*X(53)-JVS(339)*X(54)-JVS(340)*X(55)-JVS(341)*X(56)-JVS(342)&
            &*X(58)-JVS(343)*X(59))/(JVS(335))
  X(46) = (X(46)-JVS(325)*X(51)-JVS(326)*X(53)-JVS(327)*X(54)-JVS(328)*X(58)-JVS(329)*X(59))/(JVS(324))
  X(45) = (X(45)-JVS(312)*X(51)-JVS(313)*X(52)-JVS(314)*X(53)-JVS(315)*X(54)-JVS(316)*X(55)-JVS(317)*X(58)-JVS(318)&
            &*X(59))/(JVS(311))
  X(44) = (X(44)-JVS(301)*X(51)-JVS(302)*X(52)-JVS(303)*X(53)-JVS(304)*X(55)-JVS(305)*X(58))/(JVS(300))
  X(43) = (X(43)-JVS(285)*X(44)-JVS(286)*X(45)-JVS(287)*X(46)-JVS(288)*X(49)-JVS(289)*X(50)-JVS(290)*X(51)-JVS(291)&
            &*X(52)-JVS(292)*X(53)-JVS(293)*X(54)-JVS(294)*X(55)-JVS(295)*X(56)-JVS(296)*X(57)-JVS(297)*X(58)-JVS(298)&
            &*X(59))/(JVS(284))
  X(42) = (X(42)-JVS(254)*X(44)-JVS(255)*X(45)-JVS(256)*X(46)-JVS(257)*X(49)-JVS(258)*X(50)-JVS(259)*X(51)-JVS(260)&
            &*X(52)-JVS(261)*X(53)-JVS(262)*X(54)-JVS(263)*X(55)-JVS(264)*X(56)-JVS(265)*X(57)-JVS(266)*X(58)-JVS(267)&
            &*X(59))/(JVS(253))
  X(41) = (X(41)-JVS(241)*X(51)-JVS(242)*X(52)-JVS(243)*X(53)-JVS(244)*X(55)-JVS(245)*X(58))/(JVS(240))
  X(40) = (X(40)-JVS(232)*X(46)-JVS(233)*X(51)-JVS(234)*X(52)-JVS(235)*X(53)-JVS(236)*X(54)-JVS(237)*X(55)-JVS(238)&
            &*X(58))/(JVS(231))
  X(39) = (X(39)-JVS(222)*X(45)-JVS(223)*X(51)-JVS(224)*X(52)-JVS(225)*X(53)-JVS(226)*X(54)-JVS(227)*X(55)-JVS(228)&
            &*X(58)-JVS(229)*X(59))/(JVS(221))
  X(38) = (X(38)-JVS(210)*X(46)-JVS(211)*X(51)-JVS(212)*X(52)-JVS(213)*X(53)-JVS(214)*X(54)-JVS(215)*X(55)-JVS(216)&
            &*X(58))/(JVS(209))
  X(37) = (X(37)-JVS(198)*X(46)-JVS(199)*X(51)-JVS(200)*X(52)-JVS(201)*X(53)-JVS(202)*X(54)-JVS(203)*X(55)-JVS(204)&
            &*X(58))/(JVS(197))
  X(36) = (X(36)-JVS(190)*X(51)-JVS(191)*X(52)-JVS(192)*X(53)-JVS(193)*X(55)-JVS(194)*X(58))/(JVS(189))
  X(35) = (X(35)-JVS(183)*X(51)-JVS(184)*X(52)-JVS(185)*X(53)-JVS(186)*X(55)-JVS(187)*X(58))/(JVS(182))
  X(34) = (X(34)-JVS(172)*X(35)-JVS(173)*X(36)-JVS(174)*X(45)-JVS(175)*X(50)-JVS(176)*X(52)-JVS(177)*X(53)-JVS(178)&
            &*X(54)-JVS(179)*X(55)-JVS(180)*X(58))/(JVS(171))
  X(33) = (X(33)-JVS(160)*X(41)-JVS(161)*X(44)-JVS(162)*X(46)-JVS(163)*X(47)-JVS(164)*X(49)-JVS(165)*X(52)-JVS(166)&
            &*X(53)-JVS(167)*X(54)-JVS(168)*X(55)-JVS(169)*X(56)-JVS(170)*X(58))/(JVS(159))
  X(32) = (X(32)-JVS(151)*X(51)-JVS(152)*X(52)-JVS(153)*X(53)-JVS(154)*X(54)-JVS(155)*X(55)-JVS(156)*X(59))/(JVS(150))
  X(31) = (X(31)-JVS(142)*X(35)-JVS(143)*X(45)-JVS(144)*X(52)-JVS(145)*X(53)-JVS(146)*X(54)-JVS(147)*X(55)-JVS(148)&
            &*X(58))/(JVS(141))
  X(30) = (X(30)-JVS(134)*X(35)-JVS(135)*X(36)-JVS(136)*X(52)-JVS(137)*X(53)-JVS(138)*X(54)-JVS(139)*X(55)-JVS(140)&
            &*X(58))/(JVS(133))
  X(29) = (X(29)-JVS(130)*X(46)-JVS(131)*X(53)-JVS(132)*X(54))/(JVS(129))
  X(28) = (X(28)-JVS(126)*X(46)-JVS(127)*X(53)-JVS(128)*X(54))/(JVS(125))
  X(27) = (X(27)-JVS(122)*X(46)-JVS(123)*X(53)-JVS(124)*X(54))/(JVS(121))
  X(26) = (X(26)-JVS(107)*X(27)-JVS(108)*X(28)-JVS(109)*X(29)-JVS(110)*X(31)-JVS(111)*X(34)-JVS(112)*X(43)-JVS(113)&
            &*X(45)-JVS(114)*X(46)-JVS(115)*X(48)-JVS(116)*X(52)-JVS(117)*X(53)-JVS(118)*X(54)-JVS(119)*X(55)-JVS(120)&
            &*X(58))/(JVS(106))
  X(25) = (X(25)-JVS(96)*X(30)-JVS(97)*X(31)-JVS(98)*X(34)-JVS(99)*X(43)-JVS(100)*X(48)-JVS(101)*X(51)-JVS(102)*X(53)&
            &-JVS(103)*X(54)-JVS(104)*X(59))/(JVS(95))
  X(24) = (X(24)-JVS(90)*X(46)-JVS(91)*X(53)-JVS(92)*X(54))/(JVS(89))
  X(23) = (X(23)-JVS(86)*X(53)-JVS(87)*X(55)-JVS(88)*X(59))/(JVS(85))
  X(22) = (X(22)-JVS(83)*X(53)-JVS(84)*X(54))/(JVS(82))
  X(21) = (X(21)-JVS(77)*X(51)-JVS(78)*X(52)-JVS(79)*X(53))/(JVS(76))
  X(20) = (X(20)-JVS(73)*X(51)-JVS(74)*X(53)-JVS(75)*X(59))/(JVS(72))
  X(19) = (X(19)-JVS(69)*X(46)-JVS(70)*X(54)-JVS(71)*X(59))/(JVS(68))
  X(18) = (X(18)-JVS(64)*X(51)-JVS(65)*X(53)-JVS(66)*X(55))/(JVS(63))
  X(17) = (X(17)-JVS(59)*X(28)-JVS(60)*X(29)-JVS(61)*X(46)-JVS(62)*X(53))/(JVS(58))
  X(16) = (X(16)-JVS(57)*X(53))/(JVS(56))
  X(15) = (X(15)-JVS(54)*X(54)-JVS(55)*X(59))/(JVS(53))
  X(14) = (X(14)-JVS(51)*X(51)-JVS(52)*X(53))/(JVS(50))
  X(13) = (X(13)-JVS(48)*X(53)-JVS(49)*X(58))/(JVS(47))
  X(12) = (X(12)-JVS(45)*X(45)-JVS(46)*X(59))/(JVS(44))
  X(11) = (X(11)-JVS(43)*X(53))/(JVS(42))
  X(10) = (X(10)-JVS(41)*X(53))/(JVS(40))
  X(9) = (X(9)-JVS(39)*X(53))/(JVS(38))
  X(8) = (X(8)-JVS(37)*X(53))/(JVS(36))
  X(7) = (X(7)-JVS(35)*X(46))/(JVS(34))
  X(6) = (X(6)-JVS(33)*X(53))/(JVS(32))
  X(5) = (X(5)-JVS(31)*X(53))/(JVS(30))
  X(4) = (X(4)-JVS(28)*X(26)-JVS(29)*X(53))/(JVS(27))
  X(3) = (X(3)-JVS(11)*X(27)-JVS(12)*X(28)-JVS(13)*X(29)-JVS(14)*X(37)-JVS(15)*X(38)-JVS(16)*X(40)-JVS(17)*X(41)-JVS(18)&
           &*X(44)-JVS(19)*X(45)-JVS(20)*X(46)-JVS(21)*X(49)-JVS(22)*X(50)-JVS(23)*X(52)-JVS(24)*X(55)-JVS(25)*X(56)-JVS(26)&
           &*X(57))/(JVS(10))
  X(2) = (X(2)-JVS(5)*X(24)-JVS(6)*X(27)-JVS(7)*X(28)-JVS(8)*X(29)-JVS(9)*X(46))/(JVS(4))
  X(1) = (X(1)-JVS(2)*X(5)-JVS(3)*X(53))/(JVS(1))

END SUBROUTINE radm2sorg_KppSolve

! End of radm2sorg_KppSolve function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! BLAS_UTIL - BLAS-LIKE utility functions
!   Arguments :
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

!--------------------------------------------------------------
!
! BLAS/LAPACK-like subroutines used by the integration algorithms
! It is recommended to replace them by calls to the optimized
!      BLAS/LAPACK library for your machine
!
!  (C) Adrian Sandu, Aug. 2004
!      Virginia Polytechnic Institute and State University
!--------------------------------------------------------------

!--------------------------------------------------------------
      SUBROUTINE radm2sorg_WCOPY(N,X,incX,Y,incY)
!--------------------------------------------------------------
!     copies a vector, x, to a vector, y:  y <- x
!     only for incX=incY=1
!     after BLAS
!     replace this by the function from the optimized BLAS implementation:
!         CALL  SCOPY(N,X,1,Y,1)   or   CALL  DCOPY(N,X,1,Y,1)
!--------------------------------------------------------------
!     USE radm2sorg_Precision

      INTEGER i,incX,incY,M,MP1,N
      REAL(kind=dp) X(N),Y(N)

      IF (N.LE.0) RETURN

      M = MOD(N,8)
      IF( M .NE. 0 ) THEN
        DO i = 1,M
          Y(i) = X(i)
        END DO
        IF( N .LT. 8 ) RETURN
      END IF
      MP1 = M+1
      DO i = MP1,N,8
        Y(i) = X(i)
        Y(i + 1) = X(i + 1)
        Y(i + 2) = X(i + 2)
        Y(i + 3) = X(i + 3)
        Y(i + 4) = X(i + 4)
        Y(i + 5) = X(i + 5)
        Y(i + 6) = X(i + 6)
        Y(i + 7) = X(i + 7)
      END DO

      END SUBROUTINE radm2sorg_WCOPY

!--------------------------------------------------------------
      SUBROUTINE radm2sorg_WAXPY(N,Alpha,X,incX,Y,incY)
!--------------------------------------------------------------
!     constant times a vector plus a vector: y <- y + Alpha*x
!     only for incX=incY=1
!     after BLAS
!     replace this by the function from the optimized BLAS implementation:
!         CALL SAXPY(N,Alpha,X,1,Y,1) or  CALL DAXPY(N,Alpha,X,1,Y,1)
!--------------------------------------------------------------
!      USE radm2sorg_Precision

      INTEGER i,incX,incY,M,MP1,N
      REAL(kind=dp) X(N),Y(N),Alpha
      REAL(kind=dp) ZERO
      PARAMETER( ZERO = 0.0_dp )

      IF (Alpha .EQ. ZERO) RETURN
      IF (N .LE. 0) RETURN

      M = MOD(N,4)
      IF( M .NE. 0 ) THEN
        DO i = 1,M
          Y(i) = Y(i) + Alpha*X(i)
        END DO
        IF( N .LT. 4 ) RETURN
      END IF
      MP1 = M + 1
      DO i = MP1,N,4
        Y(i) = Y(i) + Alpha*X(i)
        Y(i + 1) = Y(i + 1) + Alpha*X(i + 1)
        Y(i + 2) = Y(i + 2) + Alpha*X(i + 2)
        Y(i + 3) = Y(i + 3) + Alpha*X(i + 3)
      END DO

      END SUBROUTINE radm2sorg_WAXPY

!--------------------------------------------------------------
      SUBROUTINE radm2sorg_WSCAL(N,Alpha,X,incX)
!--------------------------------------------------------------
!     constant times a vector: x(1:N) <- Alpha*x(1:N)
!     only for incX=incY=1
!     after BLAS
!     replace this by the function from the optimized BLAS implementation:
!         CALL SSCAL(N,Alpha,X,1) or  CALL DSCAL(N,Alpha,X,1)
!--------------------------------------------------------------
!      USE radm2sorg_Precision

      INTEGER i,incX,M,MP1,N
      REAL(kind=dp) X(N),Alpha
      REAL(kind=dp) ZERO, ONE
      PARAMETER( ZERO = 0.0_dp )
      PARAMETER( ONE  = 1.0_dp )

      IF (Alpha .EQ. ONE) RETURN
      IF (N .LE. 0) RETURN

      M = MOD(N,5)
      IF( M .NE. 0 ) THEN
        IF (Alpha .EQ. (-ONE)) THEN
          DO i = 1,M
            X(i) = -X(i)
          END DO
        ELSEIF (Alpha .EQ. ZERO) THEN
          DO i = 1,M
            X(i) = ZERO
          END DO
        ELSE
          DO i = 1,M
            X(i) = Alpha*X(i)
          END DO
        END IF
        IF( N .LT. 5 ) RETURN
      END IF
      MP1 = M + 1
      IF (Alpha .EQ. (-ONE)) THEN
        DO i = MP1,N,5
          X(i)     = -X(i)
          X(i + 1) = -X(i + 1)
          X(i + 2) = -X(i + 2)
          X(i + 3) = -X(i + 3)
          X(i + 4) = -X(i + 4)
        END DO
      ELSEIF (Alpha .EQ. ZERO) THEN
        DO i = MP1,N,5
          X(i)     = ZERO
          X(i + 1) = ZERO
          X(i + 2) = ZERO
          X(i + 3) = ZERO
          X(i + 4) = ZERO
        END DO
      ELSE
        DO i = MP1,N,5
          X(i)     = Alpha*X(i)
          X(i + 1) = Alpha*X(i + 1)
          X(i + 2) = Alpha*X(i + 2)
          X(i + 3) = Alpha*X(i + 3)
          X(i + 4) = Alpha*X(i + 4)
        END DO
      END IF

      END SUBROUTINE radm2sorg_WSCAL

!--------------------------------------------------------------
      REAL(kind=dp) FUNCTION radm2sorg_WLAMCH( C )
!--------------------------------------------------------------
!     returns epsilon machine
!     after LAPACK
!     replace this by the function from the optimized LAPACK implementation:
!          CALL SLAMCH('E') or CALL DLAMCH('E')
!--------------------------------------------------------------
!      USE radm2sorg_Precision

      CHARACTER C
      INTEGER   i
      REAL(kind=dp)  ONE, HALF, Eps, Sum
      PARAMETER (ONE  = 1.0_dp)
      PARAMETER (HALF = 0.5_dp)
      LOGICAL   First
      SAVE     First, Eps
      DATA     First /.TRUE./

      IF (First) THEN
        First = .FALSE.
        Eps = HALF**(16)
        DO i = 17, 80
          Eps = Eps*HALF
          CALL radm2sorg_WLAMCH_ADD(ONE,Eps,Sum)
          IF (Sum.LE.ONE) GOTO 10
        END DO
        PRINT*,'ERROR IN WLAMCH. EPS < ',Eps
        RETURN
10      Eps = Eps*2
        i = i-1
      END IF

      radm2sorg_WLAMCH = Eps

      END FUNCTION radm2sorg_WLAMCH

      SUBROUTINE radm2sorg_WLAMCH_ADD( A, B, Sum )
!      USE radm2sorg_Precision

      REAL(kind=dp) A, B, Sum
      Sum = A + B

      END SUBROUTINE radm2sorg_WLAMCH_ADD
!--------------------------------------------------------------

!--------------------------------------------------------------
      SUBROUTINE radm2sorg_SET2ZERO(N,Y)
!--------------------------------------------------------------
!     copies zeros into the vector y:  y <- 0
!     after BLAS
!--------------------------------------------------------------

      INTEGER ::  i,M,MP1,N
      REAL(kind=dp) ::  Y(N)
      REAL(kind=dp), PARAMETER :: ZERO = 0.0d0

      IF (N.LE.0) RETURN

      M = MOD(N,8)
      IF( M .NE. 0 ) THEN
        DO i = 1,M
          Y(i) = ZERO
        END DO
        IF( N .LT. 8 ) RETURN
      END IF
      MP1 = M+1
      DO i = MP1,N,8
        Y(i)     = ZERO
        Y(i + 1) = ZERO
        Y(i + 2) = ZERO
        Y(i + 3) = ZERO
        Y(i + 4) = ZERO
        Y(i + 5) = ZERO
        Y(i + 6) = ZERO
        Y(i + 7) = ZERO
      END DO

      END SUBROUTINE radm2sorg_SET2ZERO

!--------------------------------------------------------------
      REAL(kind=dp) FUNCTION radm2sorg_WDOT (N, DX, incX, DY, incY)
!--------------------------------------------------------------
!     dot produce: wdot = x(1:N)*y(1:N)
!     only for incX=incY=1
!     after BLAS
!     replace this by the function from the optimized BLAS implementation:
!         CALL SDOT(N,X,1,Y,1) or  CALL DDOT(N,X,1,Y,1)
!--------------------------------------------------------------
!      USE messy_mecca_kpp_Precision
!--------------------------------------------------------------
      IMPLICIT NONE
      INTEGER :: N, incX, incY
      REAL(kind=dp) :: DX(N), DY(N)

      INTEGER :: i, IX, IY, M, MP1, NS

      radm2sorg_WDOT = 0.0D0
      IF (N .LE. 0) RETURN
      IF (incX .EQ. incY) IF (incX-1) 5,20,60
!
!     Code for unequal or nonpositive increments.
!
    5 IX = 1
      IY = 1
      IF (incX .LT. 0) IX = (-N+1)*incX + 1
      IF (incY .LT. 0) IY = (-N+1)*incY + 1
      DO i = 1,N
        radm2sorg_WDOT = radm2sorg_WDOT + DX(IX)*DY(IY)
        IX = IX + incX
        IY = IY + incY
      END DO
      RETURN
!
!     Code for both increments equal to 1.
!
!     Clean-up loop so remaining vector length is a multiple of 5.
!
   20 M = MOD(N,5)
      IF (M .EQ. 0) GO TO 40
      DO i = 1,M
         radm2sorg_WDOT = radm2sorg_WDOT + DX(i)*DY(i)
      END DO
      IF (N .LT. 5) RETURN
   40 MP1 = M + 1
      DO i = MP1,N,5
          radm2sorg_WDOT = radm2sorg_WDOT + DX(i)*DY(i) + DX(i+1)*DY(i+1) +&
                   DX(i+2)*DY(i+2) +  &
                   DX(i+3)*DY(i+3) + DX(i+4)*DY(i+4)
      END DO
      RETURN
!
!     Code for equal, positive, non-unit increments.
!
   60 NS = N*incX
      DO i = 1,NS,incX
        radm2sorg_WDOT = radm2sorg_WDOT + DX(i)*DY(i)
      END DO

      END FUNCTION radm2sorg_WDOT
! End of BLAS_UTIL function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

END MODULE radm2sorg_Integrator