subroutine GW_AMatrix_RHS(t_s,N,RHS_vector, &
                          temp_solution,A_Matrix)
!subroutine GW_AMatrix_RHS(t_s,N,RECH,DEEPRECH,RHS_vector,old_vector,residual_vector, &
!                          temp_solution,diff_vector,norm_vector,A_Matrix,ierr)
#include <petsc/finclude/petsc.h>
#include <petsc/finclude/petscvec.h>
#include <petsc/finclude/petscmat.h>

    use petsc
    use petscvec
    use petscmat
    use GW_constants
    use GW_param_by_user
    use GW_variables

    implicit none

    integer,intent(in)             :: N
    real,intent(in)                :: t_s
   ! real,allocatable,dimension(:,:),intent(in)      :: RECH
   ! real,allocatable,dimension(:,:),intent(in)      :: DEEPRECH

    Vec, intent(inout) ::RHS_vector
    Vec, intent(inout) ::temp_solution
    Mat, intent(inout) ::A_Matrix


    real,allocatable,dimension(:,:,:) :: HCOF ! Head Coefficent
    real,allocatable,dimension(:,:,:) :: B !The constant cofficient from all exteral sources
    real,allocatable,dimension(:,:,:) :: P !Head dependent cofficient from all external sources

    ! Declare PETSc objects
    PetscErrorCode   ierr
    PetscInt   i, j, k, indexGW
    PetscScalar A_value,h,RHS

   if (.not. allocated (HCOF))then
        allocate(HCOF(NCOL,NROW,NLAY))
    end if
    if (.not. allocated(B))then
        allocate (B(NCOL,NROW,NLAY))
    end if
    if (.not. allocated(P))then
        allocate(P(NCOL,NROW,NLAY))
    end if
           ! Initialize PETSc
    !call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

    ! Create and set up the matrix A
  ! Create the matrix
 !       call MatCreate(PETSC_COMM_WORLD, A_Matrix, ierr)
 !       call MatSetSizes(A_Matrix, PETSC_DECIDE, PETSC_DECIDE, N, N, ierr)
 !       call MatSetFromOptions(A_Matrix,ierr)
 !       call MatSetType(A_Matrix, MATMPIAIJ, ierr)
  ! Create and set up the vectors for x and b
 ! call VecCreate(PETSC_COMM_WORLD, temp_solution, ierr)
 ! call VecCreate(PETSC_COMM_WORLD, RHS_vector, ierr)
 ! call VecSetSizes(temp_solution, PETSC_DECIDE, N, ierr)
 ! call VecSetSizes(RHS_vector, PETSC_DECIDE, N, ierr)
 ! call VecSetUp(temp_solution,ierr)
 ! call VecSetUp(RHS_vector,ierr)




        ! Update the matrix A and vector RHS
    do k = 1, NLAY
        do i = 1, NROW
            do j = 1, NCOL

                ! assign the old HNEW to HOLD
               ! HOLD(j,i,k)=HNEW(j,i,k) ! it is already in elevation
                ! Calculate the global index of the current cell
                !indexGW = ((k - 1) * NROW + i - 1) * NCOL + j
                indexGW= (k-1)*NCOL*NROW+(i-1)*NCOL+j
                ! Calculate P(j,i,k) and B(j,i,k)
                ! No need for P_Recharge because it does not depend on
                ! head
                P(j, i, k) = P_RIVER(j, i, k)
                B(j, i, k) = B_RIVER(j, i, k) + B_Rech(j, i) + B_pump(j,i, k)

                ! Calculate HCOF(j,i,k)
                if (k == 1) then
                    HCOF(j, i, k) = P(j, i, k) - (SY(j, i) * DELTAC *DELTAR * DELTAV(j, i, k) / DELTAT)
                else
                    HCOF(j, i, k) = P(j, i, k) - (SS(j, i) * DELTAC *DELTAR * DELTAV(j, i, k) / DELTAT)
                end if

                ! update diagonal element of A_Matrix
                     if(indexGW>=1 .and. indexGW<N)then
                    A_value = -(CC(j, i, k) + CR(j, i,k) + CV(j, i, k)) + HCOF(j, i, k)
                     call MatSetValue(A_Matrix,indexGW-1,indexGW-1,A_value,INSERT_VALUES,ierr)
                !     end if
                ! Update the neighboring elements of A_Matrix
                if (j < NCOL) then
                    A_value  = CC(j, i, k)
                   call MatSetValue(A_Matrix,indexGW-1,indexGW+1-1,A_value,INSERT_VALUES,ierr)
                end if

                if (j > 1) then
                    A_value  = CC(j - 1, i, k)
                call MatSetValue(A_Matrix,indexGW-1,indexGW-1-1,A_value,INSERT_VALUES,ierr)
                end if

                if (i < NROW) then
                    A_value = CR(j,i,k)
                call MatSetValue(A_Matrix,indexGW-1,indexGW+NCOL-1,A_value,INSERT_VALUES,ierr)
                end if

                if (i > 1) then
                    A_value  = CR(j, i - 1, k)
                call MatSetValue(A_Matrix,indexGW-1,indexGW-NCOL-1,A_value,INSERT_VALUES,ierr)

                end if

                if (k < NLAY) then
                    A_value = CV(j, i, k)
                call MatSetValue(A_Matrix,indexGW-1,indexGW+NCOL*NROW-1,A_value,INSERT_VALUES,ierr)
                end if

                if (k > 1) then
                    A_value = CV(j, i, k - 1)
                call MatSetValue(A_Matrix,indexGW-1,indexGW-NCOL*NROW-1,A_value,INSERT_VALUES,ierr)
                end if

                ! Calculate RHS based on conditions
                ! If the confined aquifer is fully saturated (the
                ! piezometric head higher than the top of the confined
                ! aquifer)
                if (HOLD(j, i, 2) > (GWBOT(j, i, (k - 1)) - BED_THK))then
                    if (k == 1) then
                        RHS = -B(j, i, k) - ((SY(j, i) *DELTAC * DELTAR * DELTAV(j, i, k) * HOLD(j, i, k) / DELTAT))
                   call VecSetValue(RHS_vector, indexGW-1, RHS, INSERT_VALUES, ierr)
                    else

                        RHS = -B(j,i,k) -((SS(j,i)*DELTAC*DELTAR*DELTAV(j,i,k)*HOLD(j,i,k)/DELTAT))
                        call VecSetValue(RHS_vector, indexGW-1, RHS, INSERT_VALUES, ierr)

                   end if
               end if

               !! if the confined aquifer has dewatered (the peizomteric
               !head is below the top of the confined aquifer, start
               !behave as unconfined)
               if (HOLD (j,i,2) .LT. (GWBOT (j,i,(k-1))-BED_THK)) then
                   if ( k .EQ. 1 ) then

                     RHS = -B(j,i,k) -((SY(j,i)*DELTAC*DELTAR*DELTAV(j,i,k)*HOLD(i,j,k)/DELTAT))+ CV(j,i,k)*(HOLD (j,i,k+1)-GWTOP(j,i,k+1))
                       call VecSetValue(RHS_vector, indexGW-1, RHS,INSERT_VALUES, ierr)
                   else

                     RHS = -B(j,i,k) -((SS(j,i)*DELTAC*DELTAR*DELTAV(j,i,k)*HOLD(j,i,k)/DELTAT))+ CV(j,i,k-1)*(GWTOP(j,i,k)-HOLD(j,i,k))
                    call VecSetValue(RHS_vector, indexGW-1, RHS,INSERT_VALUES, ierr)
                   end if
               end if
              ! this the temporary h vector that we want to find
              h =HOLD(j,i,k)
                call VecSetValue(temp_solution, indexGW-1, h,INSERT_VALUES, ierr)
                end if
           end do
       end do
   end do

   ! Set values in A_Martrix (A_Matrix)
   call MatAssemblyBegin(A_Matrix,MAT_FINAL_ASSEMBLY, ierr)
   call MatAssemblyEnd(A_Matrix,MAT_FINAL_ASSEMBLY ,ierr)
  ! Set values in the solution vector (temp_solution)
   call VecAssemblyBegin(temp_solution, ierr)
   call VecAssemblyEnd(temp_solution, ierr)

   ! Set values in the right-hand side vector (RHS_vector)
   call VecAssemblyBegin(RHS_vector, ierr)
   call VecAssemblyEnd(RHS_vector, ierr)


end subroutine GW_AMatrix_RHS