! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! Main Program File
! 
! Generated by KPP-2.2.3 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                 : cbm_Main.f90
! Time                 : Fri Mar 15 14:04:58 2013
! Working directory    : /home/sandu/kpp-2.2.3/examples/Cbm_fortran
! Equation file        : cbm.kpp
! Output root filename : cbm
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! MAIN - Main program - driver routine
!   Arguments :
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

PROGRAM cbm_Driver

  USE cbm_Model
  USE cbm_Initialize, ONLY: Initialize

      REAL(kind=dp) :: T, REF(NVAR), threshold, error
      INTEGER :: i, k, itest
      REAL(kind=dp) :: RCNTRL(20), RSTATUS(20)
      INTEGER       :: ICNTRL(20), ISTATUS(20), IERR
      REAL :: time_start, time_finish

      
!~~~> Choose a particular Rosenbrock method      
      ICNTRL(3) = 7 
!        = 0 :    Rodas3 (default)
!        = 1 :    Ros2
!        = 2 :    Ros3
!        = 3 :    Ros4
!        = 4 :    Rodas3
!        = 5 :    Rodas4
!        = 6 :    Rang3
!        = 7 :    Rok4a
!        = 8 :    Rok4b
!        = 9 :    Rok4p

!~~~> Open order results file
      OPEN(10,file='Errors.m')
      WRITE(10,*) '% Steps   CPU-time   Error'

  
test: DO itest = 1,20  
  
!~~~> Initialization 
      STEPMIN = 0.0d0
      STEPMAX = 0.0d0
      RCNTRL(1:20) = 0.0_dp
      ICNTRL(1:20) = 0
      CALL Initialize()

!~~~> Choose tolerances here 
      DO i=1,NVAR
        RTOL(i) = 10_dp**(-1.d0-dble(itest)/4)
        ATOL(i) = 1.0d-1
      END DO     


      PRINT*, ' '
      PRINT*, 'vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv'
      WRITE(6,993) RTOL(1)

      T = TSTART
        TIME = T
        WRITE(6,991)  T/3600_dp,       &
                   ( TRIM(SPC_NAMES(MONITOR(i))),           &
                     C(MONITOR(i))/CFACTOR, i=1,NMONITOR )

        CALL Update_SUN() 
        CALL Update_RCONST()

        CALL cpu_time(time_start)
        CALL Rock(NVAR,VAR,T,TEND, ATOL,RTOL, &
              RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)
        CALL cpu_time(time_finish)

        T = RSTATUS(1)

      WRITE(6,991) T/3600_dp,     &
               ( TRIM(SPC_NAMES(MONITOR(i))),           &
                 C(MONITOR(i))/CFACTOR, i=1,NMONITOR ) 
      
      ! Read reference and compute error
      OPEN(15,file='Reference.m')
      READ(15,996)  ( REF(i), i=1,NVAR )
      CLOSE(15)
      k = 0
      error = 0_dp
      threshold = 1.d+11
      DO i=1,NVAR
         ! The following line prints errors for individual species
         ! print*,SPC_NAMES(i), ABS( VAR(i)-REF(i) )/MAX( REF(i), threshold )
         IF ( REF(i) > threshold ) THEN
            error = error + ( ( VAR(i)-REF(i) )/MAX( REF(i), threshold ) )**2
            k = k+1
         END IF
      END DO
      error = SQRT( error/k )
      
      ! Accepted steps, time, and error
      WRITE(10,999) ISTATUS(Nacc), time_finish - time_start, error

      WRITE(6,997) error
      WRITE(6,995) time_finish - time_start
      WRITE(6,992) ISTATUS(Nfun), ISTATUS(Ndec), ISTATUS(Nstp), ISTATUS(Nacc), ISTATUS(Nrej)
      WRITE(6,994) RTOL(1)
      PRINT*, '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
      PRINT*, ' '

      END DO test

!~~~> Close order results file
      CLOSE(10)

991   FORMAT(' ... T=',F9.3,'h  ',200(A,'=',ES11.4,'; '))
992   FORMAT(' > Func=',I8,'; Decomp=',I8,'; Steps=',I8,' [Accepted=',I8,', Rejected=',I5,']')
993   FORMAT('  Start simulation with RTOL =',ES9.3)
994   FORMAT('  END simulation with RTOL =',ES9.3)
995   FORMAT(' > CPU time =',ES12.5,' seconds')
996   FORMAT(1000(E24.14,2X))
997   FORMAT(' > Solution error =',ES12.5)
999   FORMAT(I8,ES12.5,ES12.5)

END PROGRAM cbm_Driver

! End of MAIN function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~