SUBROUTINE petsc_solve(nvar,ncells,rtol,atol,dtol,iter,algo,prec,BJ,AIJ,delta,res)
      IMPLICIT NONE
      INTEGER*4,    INTENT(IN)    :: nvar, ncells, iter
      REAL(float),  INTENT(IN)    :: rtol, atol, dtol
      CHARACTER*30, INTENT(IN)    :: algo, prec
      Vec,          INTENT(INOUT) :: BJ
      Mat,          INTENT(INOUT) :: AIJ
      REAL(float),  INTENT(INOUT) :: delta(:,:)
      REAL(float),  INTENT(IN)    :: res(1:nvar,0:ncells)
      REAL(float),  POINTER       :: sol(:)
      INTEGER*4 :: i, ic, ierr, pos(1)
      Vec       :: XJ
      PC        :: pc
      KSP       :: ksp

      !Copy in the RHS vector
      DO ic = 1,ncells
         pos(1) = ic - 1
         CALL VecSetValuesBlockedLocal(BJ,1,pos,-res(:,ic),INSERT_VALUES,ierr)
      ENDDO

      !Solution vector
      CALL VecDuplicate(BJ,XJ,ierr)

      !MATRIX and VECTOR assembly
      CALL MatAssemblyBegin(AIJ,MAT_FINAL_ASSEMBLY,ierr)
      CALL MatAssemblyEnd(AIJ,MAT_FINAL_ASSEMBLY,ierr)
      CALL VecAssemblyBegin(BJ,ierr)
      CALL VecAssemblyEnd(BJ,ierr)

      !Create SOLVER
      CALL KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
      CALL KSPSetOperators(ksp,AIJ,AIJ,ierr)

      !Set the Algorithm to solve the system
      CALL KSPSetType(ksp,algo,ierr)

      !Create and set the Preconditioner
      CALL KSPGetPC(ksp,pc,ierr)
      CALL PCSetType(pc,prec,ierr)

      !Set Tolerances and other options from cmd line
      CALL KSPSetTolerances(ksp,rtol,atol,dtol,iter,ierr)
      CALL KSPSetFromOptions(ksp,ierr)
      CALL KSPSetUp(ksp,ierr)

      !Solve
      CALL KSPSolve(ksp,BJ,XJ,ierr)

      !Copy back the solution
      CALL VecGetArrayF90(XJ,sol,ierr)
      DO ic = 1, ncells
         DO i = 1, nvar
            delta(i,ic) = sol(nvar*(ic-1) + i)
         ENDDO
      ENDDO

      CALL VecRestoreArrayF90(XJ,sol,ierr)
      CALL VecDestroy(XJ,ierr)
      CALL KSPDestroy(ksp,ierr)
   ENDSUBROUTINE petsc_solve