!------------------------ BEGIN saprc99_main.f90 BEGIN -----------------------
! @file saprc99_main.f90                                                      
! @author jlinford                                                            
! @date 2017-02-27 19:38:52.825324                                            
! @brief Fortran90 driver with plplot                                         
!                                                                             
! Fortran90 driver with plplot                                                
!                                                                             
! This file was generated by Kppa: http://www.paratools.com/Kppa              
!-----------------------------------------------------------------------------



PROGRAM main

    USE saprc99_parameters
    USE saprc99_integrate
    USE saprc99_initialize
    USE saprc99_driver_parameters
    USE saprc99_monitor

    IMPLICIT NONE

    !---------------------------------------------------------------------------
    ! Lookat data                                                               
    !---------------------------------------------------------------------------
    
    INTEGER, PARAMETER :: lookatFile = 100
    CHARACTER(80), PARAMETER :: fname = "saprc99.dat"
    

    !---------------------------------------------------------------------------
    ! Program data                                                               
    !---------------------------------------------------------------------------
    
  ! Integration status code 
  INTEGER :: retval
  ! Integration status code per timestep iteration 
  INTEGER :: tsretval
  ! Integration time 
  REAL(8) :: time

  ! Absolute integration tolerances for variable species 
  REAL(8) :: abstol(74)
  ! Relative integration tolerances for variable species 
  REAL(8) :: reltol(74)

  ! Integer integration in/out parameters 
  INTEGER :: idata(20)
  ! Real value integration in/out parameters 
  REAL(8) :: rdata(20)

    ! Species concentrations for all grid cells
    REAL(8), ALLOCATABLE :: conc(:,:)

    ! Last stepsize in each grid cell
    REAL(8), ALLOCATABLE :: lastH(:)
    
    ! Integrator statistics
    INTEGER :: nFun = 0       ! Number of function evaluations
    INTEGER :: nJac = 0       ! Number of Jacobian evaluations
    INTEGER :: nStp = 0       ! Number of solver steps
    INTEGER :: nAcc = 0       ! Number of accepted steps
    INTEGER :: nRej = 0       ! Number of rejected steps
    INTEGER :: nDec = 0       ! Number of matrix decompositions
    INTEGER :: nSol = 0       ! Number of Ax=b solves
    INTEGER :: nSng = 0       ! Number of singular decomposition results
    
    ! Iterator
    INTEGER :: i
   
    !---------------------------------------------------------------------------
    ! Program text                                                               
    !---------------------------------------------------------------------------

    WRITE(*,101) NCELLS,TSTART,TEND
101 FORMAT("Kppa: Integrating ",I8," grid cells, time interval [",E8.2,",",E8.2,"]")

    ! Initialize example grid data
    ALLOCATE(conc(NSPEC,NCELLS))
    DO i=1,ncells
        CALL Initialize(conc(1:NVAR,i), conc(NVAR+1:NSPEC,i))
    END DO
    
    ! Allocate storage for last step size (useful in reentrant solvers)
    ALLOCATE(lastH(ncells))

    ! Set time interval
    time = TSTART

    ! Initialize vector tolerances
    abstol = ATOLS
    reltol = RTOLS

    ! Initialize integration parameters to 0 = use defaults
    idata = 0
    rdata = 0

    ! Rosenbrock default parameters
!    idata(1) = 0       ! System is non-autonomous: F = F(t,y)
!    idata(2) = 0       ! Use vector tolerances
    idata(3) = 200000  ! Maximum number of integration steps
    idata(4) = 4       ! Rodas4 Rosenbrock method
!    idata(5) = 0       ! Tolerance vectors will not be checked
!
!    rdata(1) = 0       ! Integration step size lower bound: 0 recommended
!    rdata(2) = 0       ! Integration step size upper bound: 0 recommended
!    rdata(3) = TDEL    ! Starting integration step size
!    rdata(4) = 0.2     ! Lower bound on step decrease factor
!    rdata(5) = 6       ! Upper bound on step increase factor
!    rdata(6) = 0.1     ! Step decrease factor after step rejection
!    rdata(7) = 0.9     ! Safety factor in the computation of new step size

    ! Initialize LookAt file
    OPEN(unit=lookatFile,file=fname)

    CALL LookAt(rdata(12), rdata(13))

    ! Initialize species monitor
  CALL InitMonitor(rdata(3))



    ! Perform time integration
    retval = 0
    DO WHILE (time < TEND)
        tsretval = GridIntegrate(NCELLS, conc, time, time+TDEL, &
                                 abstol, reltol, idata, rdata, lastH)

        ! Accumulate statistics
        nStp = nStp + idata(13)
        nAcc = nAcc + idata(14)
        nRej = nRej + idata(15)
        nSng = nSng + idata(18) 
        nFun = nFun + idata(11)
        nJac = nJac + idata(12)
        nDec = nDec + idata(16)
        nSol = nSol + idata(17)

        ! Update time
        time = rdata(11)

        ! Write species concentrations to file
      CALL LookAt(rdata(12), rdata(13))

        
        ! Report species concentrations to stdout
      CALL Monitor(rdata(12), rdata(13))

        
        ! Process return code
        IF (tsretval < 0) THEN
            ! Error occured: save code and abort
            retval = MIN(retval, tsretval)
            EXIT
        ELSE IF (tsretval > 0) THEN
            ! Warning occured: save code and continue
            IF (retval >= 0 .AND. retval < tsretval) THEN
                retval = tsretval;
            END IF
        END IF

        ! Use last integration step size as the starting step size
        rdata(3) = rdata(12)
    END DO

    ! Report statistics                       
    WRITE(*,102) nStp,nAcc,nRej,nSng,nFun,nJac,nDec,nSol
102 FORMAT("Kppa: Stp=",I12,", Acc=",I12,", Rej=",I12,", Sng=",I12, &
               ", Fun=",I12,", Jac=",I12,", Dec=",I12,", Sol=",I12)

    ! Clean up and exit
    CLOSE(lookatFile)
    DEALLOCATE(conc)
    DEALLOCATE(lastH)
    CALL EXIT(retval)
    STOP


  CONTAINS

!---------------------------------- Monitor ----------------------------------
! Reports on species concentrations as they are calculated                    
!                                                                             
! @param[in]     step Last accepted integrator time step (H)                  
! @param[in]     err  Error in solution vector                                
!-----------------------------------------------------------------------------
  SUBROUTINE Monitor(step, err)
    IMPLICIT NONE

    REAL(8), INTENT(IN) :: step
    REAL(8), INTENT(IN) :: err

        REAL(8) :: perc
        REAL(8) :: cspc(NMONITOR)

        cspc(1) = conc(1,1) / CFACTOR
        cspc(2) = conc(2,1) / CFACTOR
        cspc(3) = conc(3,1) / CFACTOR
        cspc(4) = conc(4,1) / CFACTOR
        cspc(5) = conc(5,1) / CFACTOR
        cspc(6) = conc(6,1) / CFACTOR
        cspc(7) = conc(7,1) / CFACTOR
        cspc(8) = conc(8,1) / CFACTOR
        cspc(9) = conc(9,1) / CFACTOR
        cspc(10) = conc(10,1) / CFACTOR
        cspc(11) = conc(11,1) / CFACTOR
        cspc(12) = conc(12,1) / CFACTOR
        cspc(13) = conc(13,1) / CFACTOR
        cspc(14) = conc(14,1) / CFACTOR
        cspc(15) = conc(15,1) / CFACTOR
        cspc(16) = conc(16,1) / CFACTOR
        cspc(17) = conc(17,1) / CFACTOR
        cspc(18) = conc(18,1) / CFACTOR
        cspc(19) = conc(19,1) / CFACTOR
        cspc(20) = conc(20,1) / CFACTOR
        cspc(21) = conc(21,1) / CFACTOR
        cspc(22) = conc(22,1) / CFACTOR
        cspc(23) = conc(23,1) / CFACTOR
        cspc(24) = conc(24,1) / CFACTOR
        cspc(25) = conc(25,1) / CFACTOR
        cspc(26) = conc(26,1) / CFACTOR
        cspc(27) = conc(27,1) / CFACTOR
        cspc(28) = conc(28,1) / CFACTOR
        cspc(29) = conc(29,1) / CFACTOR
        cspc(30) = conc(30,1) / CFACTOR
        cspc(31) = conc(31,1) / CFACTOR
        cspc(32) = conc(32,1) / CFACTOR
        cspc(33) = conc(33,1) / CFACTOR
        cspc(34) = conc(34,1) / CFACTOR
        cspc(35) = conc(35,1) / CFACTOR
        cspc(36) = conc(36,1) / CFACTOR
        cspc(37) = conc(37,1) / CFACTOR
        cspc(38) = conc(38,1) / CFACTOR
        cspc(39) = conc(39,1) / CFACTOR
        cspc(40) = conc(40,1) / CFACTOR
        cspc(41) = conc(41,1) / CFACTOR
        cspc(42) = conc(42,1) / CFACTOR
        cspc(43) = conc(43,1) / CFACTOR
        cspc(44) = conc(44,1) / CFACTOR
        cspc(45) = conc(45,1) / CFACTOR
        cspc(46) = conc(46,1) / CFACTOR
        cspc(47) = conc(47,1) / CFACTOR
        cspc(48) = conc(48,1) / CFACTOR
        cspc(49) = conc(49,1) / CFACTOR
        cspc(50) = conc(50,1) / CFACTOR
        cspc(51) = conc(51,1) / CFACTOR
        cspc(52) = conc(52,1) / CFACTOR
        cspc(53) = conc(53,1) / CFACTOR
        cspc(54) = conc(54,1) / CFACTOR
        cspc(55) = conc(55,1) / CFACTOR
        cspc(56) = conc(56,1) / CFACTOR
        cspc(57) = conc(57,1) / CFACTOR
        cspc(58) = conc(58,1) / CFACTOR
        cspc(59) = conc(59,1) / CFACTOR
        cspc(60) = conc(60,1) / CFACTOR
        cspc(61) = conc(61,1) / CFACTOR
        cspc(62) = conc(62,1) / CFACTOR
        cspc(63) = conc(63,1) / CFACTOR
        cspc(64) = conc(64,1) / CFACTOR
        cspc(65) = conc(65,1) / CFACTOR
        cspc(66) = conc(66,1) / CFACTOR
        cspc(67) = conc(67,1) / CFACTOR
        cspc(68) = conc(68,1) / CFACTOR
        cspc(69) = conc(69,1) / CFACTOR
        cspc(70) = conc(70,1) / CFACTOR
        cspc(71) = conc(71,1) / CFACTOR
        cspc(72) = conc(72,1) / CFACTOR
        cspc(73) = conc(73,1) / CFACTOR
        cspc(74) = conc(74,1) / CFACTOR
        cspc(75) = conc(75,1) / CFACTOR
        cspc(76) = conc(76,1) / CFACTOR
        cspc(77) = conc(77,1) / CFACTOR
        cspc(78) = conc(78,1) / CFACTOR
        cspc(79) = conc(79,1) / CFACTOR

        perc = ((time-TSTART)/ABS(TEND-TSTART))*100.0
        WRITE(*,301) perc,time,step,err,cspc
301     FORMAT(F5.1,"% | ",F9.2," | ",E9.3," | ",E9.3," | ",6(E11.4," |"))
  END SUBROUTINE Monitor



!-------------------------------- InitMonitor --------------------------------
! Initializes the species monitor                                             
!                                                                             
! @param[in]     step Last accepted integrator time step (H)                  
!-----------------------------------------------------------------------------
  SUBROUTINE InitMonitor(step)
    IMPLICIT NONE

    REAL(8), INTENT(IN) :: step


        WRITE(*,401)
        WRITE(*,402) "H2SO4","HCOOH","CCO_OH","RCO_OH","CCO_OOH","RCO_OOH","XN","XC","SO2","O1D","ALK1","BACL",&
                     "PAN","PAN2","PBZN","MA_PAN","H2O2","N2O5","HONO","ALK2","ALK3","TBU_O","ALK5","ARO2","HNO4",&
                     "COOH","HOCOO","BZNO2_O","MEOH","ALK4","ARO1","DCB2","DCB3","CRES","DCB1","NPHE","ROOH","BALD",&
                     "PHEN","CO","MGLY","ACET","HNO3","ETHENE","GLY","BZ_O","ISOPRENE","R2O2","TERP","METHACRO",&
                     "OLE1","ISOPROD","OLE2","MVK","CCHO","HCHO","RNO3","O3P","RCHO","MEK","PROD2","O3","HO2",&
                     "RO2_N","MA_RCO3","C_O2","BZCO_O2","RO2_R","NO","NO2","NO3","CCO_O2","RCO_O2","OH","AIR","O2",&
                     "H2O","H2","CH4"
        WRITE(*,401)
401     FORMAT(1072("-"))
402     FORMAT("       |     T     |     H     | Err. Norm | ",568(A11," |"))


    CALL Monitor(step, ZERO)

  END SUBROUTINE InitMonitor



!----------------------------------- LookAt ----------------------------------
! Writes species concentrations to file                                       
!                                                                             
! @param[in]     step Last accepted integrator time step (H)                  
! @param[in]     err  Error in solution vector                                
!-----------------------------------------------------------------------------
  SUBROUTINE LookAt(step, err)
    IMPLICIT NONE

    REAL(8), INTENT(IN) :: step
    REAL(8), INTENT(IN) :: err

    REAL(8) :: data(3 + NLOOKAT)

    write(lookatFile,"(E24.18,',',100(E24.18,','))") time, (conc(ilookat(i), 1)/CFACTOR, i=1,NLOOKAT)

!    data(1) = time
!    data(2) = step
!    data(3) = err
!    data(4) = conc(1, 1) / CFACTOR
!    data(5) = conc(2, 1) / CFACTOR
!    data(6) = conc(3, 1) / CFACTOR
!    data(7) = conc(4, 1) / CFACTOR
!    data(8) = conc(5, 1) / CFACTOR
!    data(9) = conc(6, 1) / CFACTOR
!    data(10) = conc(7, 1) / CFACTOR
!    data(11) = conc(8, 1) / CFACTOR
!    data(12) = conc(9, 1) / CFACTOR
!    data(13) = conc(10, 1) / CFACTOR
!    data(14) = conc(11, 1) / CFACTOR
!    data(15) = conc(12, 1) / CFACTOR
!    data(16) = conc(13, 1) / CFACTOR
!    data(17) = conc(14, 1) / CFACTOR
!    data(18) = conc(15, 1) / CFACTOR
!    data(19) = conc(16, 1) / CFACTOR
!    data(20) = conc(17, 1) / CFACTOR
!    data(21) = conc(18, 1) / CFACTOR
!    data(22) = conc(19, 1) / CFACTOR
!    data(23) = conc(20, 1) / CFACTOR
!    data(24) = conc(21, 1) / CFACTOR
!    data(25) = conc(22, 1) / CFACTOR
!    data(26) = conc(23, 1) / CFACTOR
!    data(27) = conc(24, 1) / CFACTOR
!    data(28) = conc(25, 1) / CFACTOR
!    data(29) = conc(26, 1) / CFACTOR
!    data(30) = conc(27, 1) / CFACTOR
!    data(31) = conc(28, 1) / CFACTOR
!    data(32) = conc(29, 1) / CFACTOR
!    data(33) = conc(30, 1) / CFACTOR
!    data(34) = conc(31, 1) / CFACTOR
!    data(35) = conc(32, 1) / CFACTOR
!    data(36) = conc(33, 1) / CFACTOR
!    data(37) = conc(34, 1) / CFACTOR
!    data(38) = conc(35, 1) / CFACTOR
!    data(39) = conc(36, 1) / CFACTOR
!    data(40) = conc(37, 1) / CFACTOR
!    data(41) = conc(38, 1) / CFACTOR
!    data(42) = conc(39, 1) / CFACTOR
!    data(43) = conc(40, 1) / CFACTOR
!    data(44) = conc(41, 1) / CFACTOR
!    data(45) = conc(42, 1) / CFACTOR
!    data(46) = conc(43, 1) / CFACTOR
!    data(47) = conc(44, 1) / CFACTOR
!    data(48) = conc(45, 1) / CFACTOR
!    data(49) = conc(46, 1) / CFACTOR
!    data(50) = conc(47, 1) / CFACTOR
!    data(51) = conc(48, 1) / CFACTOR
!    data(52) = conc(49, 1) / CFACTOR
!    data(53) = conc(50, 1) / CFACTOR
!    data(54) = conc(51, 1) / CFACTOR
!    data(55) = conc(52, 1) / CFACTOR
!    data(56) = conc(53, 1) / CFACTOR
!    data(57) = conc(54, 1) / CFACTOR
!    data(58) = conc(55, 1) / CFACTOR
!    data(59) = conc(56, 1) / CFACTOR
!    data(60) = conc(57, 1) / CFACTOR
!    data(61) = conc(58, 1) / CFACTOR
!    data(62) = conc(59, 1) / CFACTOR
!    data(63) = conc(60, 1) / CFACTOR
!    data(64) = conc(61, 1) / CFACTOR
!    data(65) = conc(62, 1) / CFACTOR
!    data(66) = conc(63, 1) / CFACTOR
!    data(67) = conc(64, 1) / CFACTOR
!    data(68) = conc(65, 1) / CFACTOR
!    data(69) = conc(66, 1) / CFACTOR
!    data(70) = conc(67, 1) / CFACTOR
!    data(71) = conc(68, 1) / CFACTOR
!    data(72) = conc(69, 1) / CFACTOR
!    data(73) = conc(70, 1) / CFACTOR
!    data(74) = conc(71, 1) / CFACTOR
!    data(75) = conc(72, 1) / CFACTOR
!    data(76) = conc(73, 1) / CFACTOR
!    data(77) = conc(74, 1) / CFACTOR
!    data(78) = conc(75, 1) / CFACTOR
!    data(79) = conc(76, 1) / CFACTOR
!    data(80) = conc(77, 1) / CFACTOR
!    data(81) = conc(78, 1) / CFACTOR
!    data(82) = conc(79, 1) / CFACTOR
!    WRITE(lookatFile,*) data
  END SUBROUTINE LookAt



END PROGRAM main


!-------------------------- END saprc99_main.f90 END -------------------------
