subroutine GW_solver(t_s,N)
!subroutine GW_solver(t_s,N,RECH,DEEPRECH)
!

#include <petsc/finclude/petsc.h>
#include <petsc/finclude/petscksp.h>

#include <petsc/finclude/petscvec.h>
#include <petsc/finclude/petscmat.h>
    use GW_constants
    use petsc
    use petscvec
    use petscmat
    use GW_param_by_user
    use GW_variables


    implicit none

    ! Declare arguments for the solver subroutine

    !Integer, intent(in) :: N,NCOL,NROW,max_iter
    real, intent(in) ::t_s
    integer, intent(in) ::N
    !real, allocatable,dimension(:,:),intent(inout) ::RECH
    !real, allocatable,dimension(:,:),intent(inout) ::DEEPRECH

    real, allocatable,dimension(:) :: H_vector

    Mat                            :: A_Matrix
    Vec                            :: RHS_vector
    Vec                            :: old_solution
    Vec                            :: temp_solution
    Vec                            :: residual_vector
    Vec                            :: diff_vector
    Vec                            :: norm_vector

    PetscInt          iter,i,j,k, indexGW
    KSP               ksp
    PC                pc
    PetscReal         tol,norm
    PetscErrorCode    ierr
    PetscReal         norm_picard
    PetscReal         diff_norm
    PetscScalar       neg_one
    PetscScalar       zero
    Vec               abs_diff_vector


    !Set tolerance
    tol = 1.0e-6
    ierr = 0.0
    zero = 0.0
    neg_one = -1.0
    !do iter = 1, max_iter

    ! Initialize PETSc
    call PetscInitialize(ierr)

    ! create AMatrix and RHS_vector and temp_solution vectors

    call MatCreate(PETSC_COMM_WORLD, A_Matrix, ierr)
    call MatSetSizes(A_Matrix, PETSC_DECIDE, PETSC_DECIDE, N, N, ierr)
    call MatSetType(A_Matrix, MATAIJ, ierr) ! as sparse matrix
    call MatSetFromOptions(A_Matrix,ierr)
    call MatSetUp(A_Matrix,ierr)

    call VecCreate(PETSC_COMM_WORLD, RHS_vector, ierr)
    call VecSetSizes(RHS_vector, PETSC_DECIDE, N, ierr)
    call VecSetFromOptions(RHS_vector,ierr)
    call VecSet(RHS_vector, zero, ierr) ! initialize RHS_vector with zero

    call VecCreate(PETSC_COMM_WORLD, temp_solution, ierr)
    call VecSetSizes(temp_solution, PETSC_DECIDE, N, ierr)
    call VecSetFromOptions(temp_solution,ierr)
    call VecSet(temp_solution, zero, ierr) ! initialize temp_solution with zero

     ! Create vectors for old solution and residual
    call VecDuplicate(temp_solution, old_solution, ierr)
    call VecDuplicate(RHS_vector, residual_vector, ierr)


    ! Create the KSP solver cotext
    ! Create the KSP solver cotext
    call KSPCreate(PETSC_COMM_WORLD, ksp, ierr)
    call KSPSetOperators(ksp, A_Matrix, A_Matrix, ierr)

    ! Set the solver type (CG for Conjugate Gradient)
    call KSPSetType(ksp, KSPCG, ierr)
    !call KSPSetType(ksp, KSP_CG_SYMMETRIC, ierr)
    ! Create and set up a preconditioner object
    call PCCreate(PETSC_COMM_WORLD, pc, ierr)
    call PCSetType(pc, PCILU,ierr) ! Choose a preconditioner type (ILU)
    call KSPSetPC(ksp, pc,ierr)  ! Associate the preconditioner with the KSP solver


    ! Set solver options (e.g tolerance, max_iteration for the solver (default))
    call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)
    call KSPSetUp(ksp,ierr)


     call MatGetSize(A_matrix,N,N,ierr)


    ! Nonlinear Picard Iteration
        do iter = 1, MAX_ITER

        call VecCopy(temp_solution, old_solution, ierr)

        ! Update the RHS_vector and A-Matrix based on the current
        ! temp solution_vector

        ! Call GW_AMatrix_RHS to updqate A_matrix and RHS_vector

    call GW_AMatrix_RHS(t_s,N,RHS_vector, &
                          temp_solution,A_Matrix)
    !call
    !GW_AMatrix_RHS(t_s,N,RECH,DEEPRECH,RHS_vector, &
    !                      temp_solution,A_Matrix)



        ! Solve the linear system using the PCG solver
        call KSPSolve(ksp, RHS_vector, temp_solution, ierr)
        !_______________________ Check for convergence_________________________!
!____________________________________________________________________________________!

        ! Compute the difference between the residual and the RHS vector
            ! Calculate the residual vector
        call MatMult(A_Matrix, temp_solution, residual_vector,ierr)

        ! Compute the difference between the residual and the RHS vector
        !call VecAXPY(residual_vector, neg_one, RHS_vector,ierr)

        ! Compute the L2 norm of the residual vector

        !call VecNorm(residual_vector, NORM_2,norm_picard, ierr)

        ! Calculate the absolute difference between temp_solution and
        ! old_solution
        call VecWAXPY(abs_diff_vector, neg_one, temp_solution, old_solution,ierr)

        ! Compute the L-infinity (maximum absolute) norm of the absolute
        ! difference vector
        call VecNorm(abs_diff_vector, NORM_INFINITY, norm_picard, ierr)

        ! Now, norm_picard contains the maximum absolute difference
            if (norm_picard .LT. 1.0E-4) exit !then

             ! exit
            !else
            !  continue
            !endif




    end do
    call KSPDestroy(ksp,ierr)
    call PCDestroy(pc, ierr) ! Destroy the preconditioner object
  ! Destroy the difference_vector when it is no longer needed
    call VecDestroy(abs_diff_vector,ierr)

    call VecDestroy(old_solution, ierr)
    call VecDestroy(residual_vector, ierr)
    ! store the values from temp_solution into H_vector as array

     call VecGetArray(temp_solution, H_vector, ierr)
     call VecDestroy(temp_solution, ierr)
     call VecDestroy(RHS_vector, ierr)
     call MatDestroy(A_Matrix,ierr)
   ! Now, you can use the 'H_vector' array in your main program or another
   ! subroutine

   ! Now, you can use the 'H_vector' array in your main program or another
   ! subroutine
        do k=1, NLAY
                do i=1, NROW
                        do j= 1, NCOL
                                !indexGW = ((k - 1) * NROW + i - 1) * NCOL + j
                                indexGW= (k-1)*NCOL*NROW+(i-1)*NCOL+j
                                HNEW (j,i,k) = H_vector (indexGW)
                                if (k==1) then
                                   UNCONF_H(j,i) = HNEW(j,i,k)
                                   WTD(j,i) = SurElev (j,i)
                                end if
                       end do
                end do
        end do
        GW_HEAD(:,:,:,t_s)=HNEW(:,:,:)
call VecRestoreArray (temp_solution, H_vector, ierr)
call PetscFinalize(ierr)
end subroutine GW_solver