module Fortran_CVODE
  IMPLICIT NONE
  contains
    
    SUBROUTINE INITALL(NEU, NNY, NX, NY, NZ, Y, YDOT, Euro, A1, A2, W1, W1A, W2, W2EXT,X)
    use CONSTANTS
    USE MathSub, ONLY: MATRIZW, RSAMPLE 
    IMPLICIT NONE

    REAL(DBL) :: Y(2*NEU), YDOT(2*NEU), W2(NEU), W1(NEU), W1A(NEU), W2EXT(NEU,NEU), X(6)
    REAL(DBL) :: A1, A2
    INTEGER :: euro(NEU)
    INTEGER :: NX, NY, NZ, NNY, NEU
    
    call  RSAMPLE(NNy, NEU, euro)
    Y=0d0
    YDOT=0
    y(neu+1:2*neu)=1d-4
    CALL MatrizW(NEU,NX,NY,NZ,euro,X,W1,W1A,W2,W2EXT)
    A1=X(1)
    A2=X(2)

    return
    END SUBROUTINE INITALL
    
    SUBROUTINE INITALLO_LIN(NEU, NO, NNY, NX, NY, NZ, Y, YDOT, Euro, A1, A2, W1, W1A, W1O, W1AO, W2, X)
    use CONSTANTS
    use GLOBALS, only: neq, EUROPIOS
    USE MathSub, ONLY: MATRIZEUO_LIN, RSAMPLE 
    IMPLICIT NONE

    REAL(DBL) :: Y(NEq), YDOT(NEq), X(6)
    REAL(DBL) :: W2(*), W1O(*), W1AO(*), W1(*), W1A(*)
    REAL(DBL) :: A1, A2
    INTEGER :: euro(NEU)
    INTEGER :: NX, NY, NZ, NNY, NEU, NO
    
    call  RSAMPLE(NNy, NEU, euro)
    Y=0d0
    YDOT=0
    y(neu+1:NEU+NO)=1d-4
    CALL MatrizEUO_lin(NEU, NO, NX, NY, NZ, EUROPIOS, X, W1, W1A, W1O, W1AO, W2)
    A1=X(1)
    A2=X(2)

    return
    END SUBROUTINE INITALLO_LIN
    
    SUBROUTINE INITALLO_NOLIN(NEU, NO, NNY, NX, NY, NZ, Y, YDOT, Euro, A1, A2, W1_ext, W1A_ext, W1O_ext, W1AO_ext, W2ext, X)
    use CONSTANTS
    use GLOBALS, only: neq, EUROPIOS
    USE MathSub, ONLY: MATRIZEUO_NOLIN, RSAMPLE 
    IMPLICIT NONE

    REAL(DBL) :: Y(NEq), YDOT(NEq), X(6)
    REAL(DBL) :: W1O_EXT(*), W1AO_EXT(*), W2EXT(*), W1_EXT(*), W1A_EXT(*)

    REAL(DBL) :: A1, A2
    INTEGER :: euro(NEU)
    INTEGER :: NX, NY, NZ, NNY, NEU, NO
    
    call  RSAMPLE(NNy, NEU, euro)
    Y=0d0
    YDOT=0
    y(neu+1:NEU+NO)=1d-4
    CALL MatrizEUO_NOlin(NEU, NO, NX, NY, NZ, EUROPIOS, X, W1_EXT, W1A_EXT, W1O_EXT, W1AO_EXT, W2EXT)
    A1=X(1)
    A2=X(2)

    return
    END SUBROUTINE INITALLO_NOLIN
    
    SUBROUTINE Output_Summary_Files(NEU, T, y)
    use CONSTANTS
    use globals
    IMPLICIT NONE
    real(dbl) :: t, y(2*neu)
    INTEGER :: NEU
    
    write(1112,*) t, sum(y(1:neu))
    
    end SUBROUTINE Output_Summary_Files
    
    SUBROUTINE LDE_Solution(number, NX, NY, NZ, concentracion, repeticiones, THETA_OUT, x, VERBOSE)
    use CONSTANTS
    use GLOBALS
    USE ExpAdj, ONLY : AJUSTE_EXPONENCIAL
    USE GLOBALS_LDE_SOLUTION, ONLY: NLINEAR, val_Exp, NLINEAR, use_jac
    IMPLICIT NONE
    LOGICAL, INTENT(IN) :: VERBOSE
    INTEGER, INTENT(IN) :: NX, NY, NZ, CONCENTRACION, REPETICIONES, number
    REAL(DBL), INTENT(in) :: x(6)
    REAL(DBL), INTENT(OUT) :: THETA_OUT
    !INTEGER, PARAMETER :: MX=10, MY=5
    !REAL(DBL), PARAMETER :: MXMY=MX*MY
    !REAL(DBL), PARAMETER :: XMAX=2.0D0, YMAX=1.0D0
    REAL(DBL) :: THETA(3)
    INTEGER :: METH, ITMETH, ITOL, ITASK, NNY
    !C The following declaration specification should match C type long int.
    INTEGER(c_long) :: IOUT(21), IPAR, bdf_order
    rEAL(DBL) :: RTOL, ATOL, T, TOUT
    rEAL(DBL) :: ROUT(6)
    REAL(DBL), ALLOCATABLE, dimension(:) :: Y, YDOT, y_acc, times_out
    REAL(DBL) :: RPAR, there
    INTEGER :: USE_PRECOND, USE_JACOBIAN, III, JJJ, IERR
    

    ierr=0

    NNY=NX*NY*NZ*32
    NNEU=NNY*CONCENTRACION/100
    NNO=NX*NY*NZ*48
    NEQ=NNO+NNEU
    
    t=0
    tout=1d-1
    !TOUT_TIME = 1d-6
    !delta_t=TOUT_TIME
    !max_time_step=TOUT_TIME
    !N_ITERATIONS=CEILING(TOUT/delta_t)
    bdf_order=2    
    USE_PRECOND=0
    USE_JACOBIAN=0
    IF (USE_JAC==1) USE_JACOBIAN= 1
    

    METH = 2 !this should be equal to 2 for stiff problems
    ITMETH = 2
    ITOL = 1
    !ALLOCATE(times_out(N_ITERATIONS), STAT=IERR)
    if (ierr/=0) PAUSE "error -1"
    times_out=0d0
    ALLOCATE(Y_acc(500), STAT=IERR)
    if (ierr/=0) PAUSE "error 0"
    y_acc=0d0
    ALLOCATE(Y(NEQ), STAT=IERR)
    if (ierr/=0) PAUSE "error 1"
    ALLOCATE(YDOT(NEQ), STAT=IERR)
    if (ierr/=0) PAUSE "error 2"
    ALLOCATE(europios(NNEU), STAT=IERR)
    if (ierr/=0) PAUSE "error 3"
    
    IF (NLINEAR==.FALSE.) THEN
        ALLOCATE (W2(NNEU),stat=ierr)
        if (ierr/=0) PAUSE "error 4"
        ALLOCATE (W1(NNEU,NNO),stat=ierr)
        if (ierr/=0) PAUSE "error 5"
        ALLOCATE (W1A(NNEU),stat=ierr)
        if (ierr/=0) PAUSE "error 6"
        ALLOCATE (W1O(NNO),stat=ierr)
        if (ierr/=0) PAUSE "error 5O"
        ALLOCATE (W1AO(NNO,NNEU),stat=ierr)
        if (ierr/=0) PAUSE "error 6O"
    ELSE
        ALLOCATE (W2eXT(NNEU,NNEU),stat=ierr)
        if (ierr/=0) PAUSE "error 4"
        ALLOCATE (W1_EXT(NNEU,NNO),stat=ierr)
        if (ierr/=0) PAUSE "error 5"
        ALLOCATE (W1A_EXT(NNEU,NNO),stat=ierr)
        if (ierr/=0) PAUSE "error 6"
        ALLOCATE (W1O_EXT(NNO,NNEU),stat=ierr)
        if (ierr/=0) PAUSE "error 5O"
        ALLOCATE (J11(NNEU),stat=ierr)
        if (ierr/=0) PAUSE "error 5O"
        ALLOCATE (J12(NNEU),stat=ierr)
        if (ierr/=0) PAUSE "error 5O"
        ALLOCATE (J21(NNO),stat=ierr)
        if (ierr/=0) PAUSE "error 5O"
        ALLOCATE (J22(NNO),stat=ierr)
        if (ierr/=0) PAUSE "error 5O"
        ALLOCATE (W1AO_eXT(NNO,NNEU),stat=ierr)
        if (ierr/=0) PAUSE "error 6O"
    END IF
        
    
    RTOL = 1D-10
    ATOL = 1D-10
    ! (unit = 1112, file = "./files/run_info.txt")  
    
    do  jjj=1,repeticiones
     !   MAXVAL_HERE=0D0
        if (verbose) write(13,*) JJJ
        T=0
        ITASK = 1
        iout=0
        rout=0
        ipar=0
        rpar=0
        ierr=0

    IF (NLINEAR==.FALSE.) THEN
        call INITALLO_LIN(nNEU, nNO, NNY, NX, NY, NZ, Y, YDOT, Europios, A1, A2, W1, W1A, W1O, W1AO, W2, X)
    ELSE
        call INITALLO_NOLIN(nNEU, nNO, NNY, NX, NY, NZ, Y, YDOT, Europios, A1, A2, W1_ext, W1A_ext, W1O_ext, W1AO_ext, W2ext, X)
    END IF    
   
        if (verbose) WRITE(13,*) "All values initialized"
        CALL FNVINITS_(1, NEQ, ierr)
        IF (ierr .NE. 0) THEN
           WRITE(6,20) ierr
    20     FORMAT(///' SUNDIALS_ERROR: FNVINITS_ returned IER = ', I5)
           READ(*,*)
           STOP
        ENDIF
        CALL FCVMALLOC_(T, y, METH, ITMETH, ITOL, RTOL, ATOL, IOUT, ROUT, IPAR, RPAR, ierr)
        IF (ierr .NE. 0) THEN
           WRITE(6,30) ierr
    30     FORMAT(///' SUNDIALS_ERROR: FCVMALLOC returned IER = ', I5)
           STOP
        ENDIF
        CALL FCVSETIIN_("MAX_NSTEPS",100000,ierr)
        CALL FCVSETIIN_("MAX_ORD",bdf_order,ierr)
        !CALL FCVSETRIN_("MAX_STEP",max_time_step,IER)
      !  CALL FCVSETIIN_("STAB_LIM",1,IER)
        CALL FCVSETRIN_("INIT_STEP", 1d-17,ierr)
        CALL FCVSETRIN_("MIN_STEP", 1d-18,ierr)

        CALL FCVSPGMR_(use_precond,1,200,0,ierr)
        !if (verbose) WRITE(*,*) "Using GMRES for linear solver"
    
        IF (ierr/=0) STOP "ERROR in the set up of the linear solver"

        !if(use_jacobian==1) then !
        !   CALL FCVSPILSSETJAC(1,IER) !COMMENT IF YOU ARE NOT PROVIDING JACOBIAN############
        !   IF (IER/=0) STOP "ERROR SETTING UP THE JACOBIAN"
        !end if

        !If (use_precond/=0) THEN
        !   CALL FCVSPILSSETPREC(1,IER) !COMMENT IF YOU ARE NOT PROVIDING THE PRECON FOR THE I- GAMMA*JACOBIAN
        !   IF (IER/=0) STOP "ERROR SETTING UP THE THE PRECON FOR THE I- GAMMA*JACOBIAN"
       !    call P_STATIC(Y, Y, Y, 1d0, 0, Y	)
       ! END IF

        IF (ierr .NE. 0) THEN
           WRITE(6,40) ierr
    40     FORMAT(///' SUNDIALS_ERROR: FCVS returned IER = ', I5)
           CALL FCVFREE_
           STOP
       ! ELSE
       !    WRITE(*,*) " "
        ENDIF
 !
    
        CALL FCVSPILSSETJAC_(USE_JACOBIAN,ierr) !COMMENT IF YOU ARE NOT PROVIDING JACOBIAN############
    
        DO iii=1,500 !N_ITERATIONS
              there=mytime(iii)!MIN(DBLE(iii-1)*DELTA_t+1d-15,tout)
              !times_out(III)=there
             ! write(*,*) "here", iii, there, y(1)
              CALL FCVSETRIN_("STOP_TIME",there,ierr)
      !        IF(verbose_mode) WRITE(*,*) "Solving until", MIN(DBLE(iii)*inp_output_time,tout), tout
              CALL FCVODE_(there, T, y, ITASK, ierr)
     !         call Output_Summary_Files(there,y)
    !          MAXVAL_HERE=MAXVAL(Y(1:NNEU),1)
    !          Y=Y/MAXVAL_HERE
              y_acc(iii)=y_acc(iii)+sum(abs(y(1:nneu)))
        END DO

      !        pause

        
    
        if (verbose) then
            NST = IOUT(LNST)
            NFE = IOUT(LNFE)
            NPSET = IOUT(LNSETUP)
            NPE = IOUT(LNPE)
            NPS = IOUT(LNPS)
            NNI = IOUT(LNNI)
            NLI = IOUT(LNLI)
            NCFN = IOUT(LNCF)
            NCFL = IOUT(LNCFL)
            NETF = IOUT(LNETF)
            NJTV = IOUT(LNJTV)

            WRITE(13,*) 'Final statistics (LLG):'
            WRITE(13,*)    ' number of steps = ', NST
            WRITE(13,*)    ' number of f evals. =', NFE
            WRITE(13,*)    ' number of prec. setups = ', NPSET
            WRITE(13,*)    ' number of prec. evals. = ', NPE
            WRITE(13,*)    ' number of prec. solves = ', NPS
            WRITE(13,*)    ' number of nonl. iters. = ', NNI
            WRITE(13,*)    ' number of lin. iters. = ', NLI
            WRITE(13,*)    ' number of JtV. = ', NJTV
            WRITE(13,*)    ' number of conv. failures.. nonlinear = ', NCFN
            WRITE(13,*)    ' linear = ', NCFL
            WRITE(13,*)    ' number of error test failures = ', NETF
        end if

        CALL FCVFREE_   
    end do
    y_acc=dabs(y_acc)!/dble(repeticiones)
    
    !if (times_here==1) tHEN
    !    output_Array(:,1)=times_out
    !    output_array(:,2)=y_acc
    !    times_here=times_here+1
    !else
    !    output_array(:,1+times_here)=y_acc
    !    times_here=times_here+1
    !end if
    !
    !call AJUSTE_EXPONENCIAL(N_ITERATIONS, times_out, y_acc, Theta, VERBOSE)
    call AJUSTE_EXPONENCIAL(500, mytime, y_acc, Theta, VERBOSE)
  !  write(*,*) Theta
    
    THETA_OUT=THETA(2)
    WRITE(*,*) number, concentracion, (theta(2)-val_exp(number))/val_exp(number)
    
    
   

    !CLOSE(1112)
    !PAUSE
    IF (NLINEAR==.FALSE.) THEN
        deallocate(y,ydot, europios,W2,W1, W1A,y_acc,W1O, W1AO)
    else
      deallocate(y,ydot, europios, W1_Ext, W1A_Ext, y_acc, W2EXT, W1O_ext, W1AO_ext,J11,J12,J21,J22)
    end if
    End SUBROUTINE  LDE_Solution
    
    
    SUBROUTINE PRE_LDE(Icon, concentracion, REPETICIONES, NX, NY, NZ, R_SUM)
    use CONSTANTS
    use MathSub, only: rsample, sumr
    use GLOBALS, only: europios, NNEU
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: NX, NY, NZ, REPETICIONES, icon, concentracion
    REAL(DBL), INTENT(OUT) :: R_SUM
    integer :: nny, ierr, iii
    real(Dbl) :: rone
       
    NNY=NX*NY*NZ*32
    NNEU=NNY*CONCENTRACION/100
    ALLOCATE(europios(NNEU), STAT=IERR)
    r_sum=0d0
    DO III=1,REPETICIONES
        call  RSAMPLE(NNy, NNEU, europios)
        CALL SUMR(nneu, NX, NY, NZ, europios, Rone)
        r_sum=r_sum+rone
    end do
    r_sum=r_sum/dble(repeticiones)
    if (ierr/=0) PAUSE "error 3"
    
    DEALLOCATE(europios)
    
    End SUBROUTINE  PRE_LDE

    SUBROUTINE findthetas(M,N,X,FVEC)
    USE GLOBALS_LDE_SOLUTION
    IMPLICIT NONE
    iNTEGER, INTENT(IN) :: M, N
    double precision, INTENT(OUT) :: FVEC(M)
    double precision, INTENT(IN) :: X(N)
    INTEGER :: III
    
    
    
    theta_b=0d0
    DO III=1,6
        CALL LDE_Solution(iii, nx, ny, nz, concentracionES(III), repeticiones, THETA_B(III), x, VERBOSE)
        FVEC(iii)=THETA_B(iii)
    END DO
    
    END SUBROUTINE findthetas
    
    
 SUBROUTINE findthetas_jacobian(M,N,X,FVEC)
    USE GLOBALS_LDE_SOLUTION
    IMPLICIT NONE
    iNTEGER, INTENT(IN) :: M, N
    double precision, INTENT(OUT) :: FVEC(M)
    double precision, INTENT(IN) :: X(N)
    INTEGER :: III
    
    
    

    DO III=1,6
        CALL LDE_Solution(iii, nx, ny, nz, concentracionES(III), repeticiones/2, FVEC(III), x, VERBOSE)
    END DO
    
END SUBROUTINE findthetas_jacobian

    
end module Fortran_CVODE

    
    