Include 'mkl_pardiso.f90'
Subroutine connectivity_table(n_y,n_z,n_v,connect,iter)
Implicit None
Integer, Intent(In) :: n_y,n_z,n_v
Integer, Intent(Out) :: connect(n_y*n_z*n_v,5*n_v),iter(n_y,n_z)
Integer :: i,j,k,m,ijk,nxyz,nxy

nxy=n_y
nxyz=nxy*n_z
Write(*,*) "connectivity_table: nxy, nxyz = ", nxy, nxyz

Do j=1,n_y
  Do k=1,n_z
    ijk = j+(k-1)*n_y
    iter(j,k)=ijk
    If (k>1) Then
      connect(ijk,1) = ijk-n_y
    Else
      connect(ijk,1) = -1
    End If

    If (j>1) Then
      connect(ijk,2) = ijk-1
    Else
      connect(ijk,2) = -2
    End If
    
    connect(ijk,3) = ijk

    If (j<n_y) Then
      connect(ijk,4) = ijk+1
    Else
      connect(ijk,4) = -3
    End If

    If (k<n_z) Then
      connect(ijk,5) = ijk+n_y
    Else
      connect(ijk,5) = -4
    End If

    If (n_v>1) Then
      Do i=1,n_v-1
        Do m=1,5
          If (connect(ijk,m)>0) Then
            connect(ijk,m+(i*5)) = connect(ijk,m) + i*nxyz
          Else
            connect(ijk,m+(i*5)) = connect(ijk,m)
          End If
        End Do
      End Do
      ForAll(i=1:n_v-1) connect(ijk+i*nxyz,:) = connect(ijk,:)
    End If
    If (Any(connect(ijk,:) > nxyz*n_v)) Then
      Write(*,*) "connectivity_table: Error!"
      Write(*,*) "j,k", j,k
      Write(*,*) "connect(ijk,:)"
      Write(*,*) connect(ijk,:)
    End If
  End Do
End Do
End Subroutine connectivity_table

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Subroutine coup_dsolve(n_y,n_z,n_v,variable,a_p,a_f,connect,mtype,stype)
Use mkl_pardiso
Implicit None
Integer, Intent(In) :: n_y,n_z,n_v                  ! Number of y points, zpoints and variables respectively
Real(Kind=8), Intent(In) :: a_p(n_y*n_z*n_v,5*n_v)  ! Control volume values
Real(Kind=8), Intent(In) :: a_f(n_y*n_z*n_v)        ! source term values
Integer, Intent(In) :: connect(n_y*n_z*n_v,5*n_v)   ! connectivity table
Real(Kind=8), Intent(Out) :: variable(n_y*n_z*n_v)  ! solution
Integer,Intent(In) :: mtype                         ! type of matrix (see pardiso manual for full list)
Logical,Intent(In) :: stype                         ! Direct (true) or CGS (false) solver
Type(MKL_PARDISO_HANDLE), Automatic :: pt(64)
Integer :: maxfct
Integer :: mnum
Integer :: phase
Integer :: ia(n_v*n_y*n_z+1)
Integer, Target :: ja(5*n_y*n_z*n_v**2)
Integer, Pointer :: jb(:)
Integer :: nrhs
Integer :: iparm(64)
Integer :: msglvl
Real(Kind=8), Dimension(n_v*n_y*n_z) :: q,fi
Real(Kind=8), Target, Dimension(5*(n_v**2)*n_y*n_z) :: mat
Real(Kind=8), Pointer :: mat2(:)
Integer :: i,j,l,nn,nn2,error
Integer :: idum(1)
Real(Kind=8) :: ddum(1)

nn=n_y*n_z*n_v
nn2=5*n_v
iparm=0

ia=0
ja=0
q=a_f
fi=0d0
mat=0d0
l=1

Do i=1,nn
  ia(i) = l
  Do j=1,nn2
    If (a_p(i,j)/=0d0.AND.connect(i,j)>0) Then
      ja(l) = connect(i,j)
      mat(l) = a_p(i,j)
      l=l+1
    End If
  End Do
End Do
      
ia(nn+1) = l
jb=>ja(1:l-1)
mat2=>mat(1:l-1)

maxfct=1
mnum=1
msglvl=0
nrhs=1

Call pardisoinit(pt,mtype,iparm)
iparm(1)=1
iparm(2)=3
If (stype) Then
  iparm(4)=0
Else
  iparm(4)=81
End If
iparm(24) = 1     ! Two level factorization (may not speed up unless there are
iparm(25) = 0     ! Parallel back substitution (use 1 for serial)
iparm(27)=1
iparm(60)=0
phase=11

!****** Program USUALLY crashes somewhere in the next line **************!
!****** Program USUALLY crashes somewhere in the next line **************!
!****** Program USUALLY crashes somewhere in the next line **************!
Call pardiso(pt,maxfct,mnum,mtype,phase,nn,mat2,ia,jb,idum,nrhs,iparm,msglvl,ddum,ddum,error)
If (error/=0) Then
  Write(*,*) "coup_dsolve: Error0!", error
  Stop
End If
phase=22
If (stype) Then
  iparm(4)=0
Else
  iparm(4)=81
End If
iparm(24) = 1     ! Two level factorization (may not speed up unless there are
iparm(25) = 0     ! Parallel back substitution (use 1 for serial)

Call pardiso(pt,maxfct,mnum,mtype,phase,nn,mat2,ia,jb,idum,nrhs,iparm,msglvl,ddum,ddum,error)
If (error/=0) Then
  Write(*,*) "coup_dsolve: Error1!", error
  Stop
End If
phase=33
iparm(24) = 1     ! Two level factorization (may not speed up unless there are
iparm(25) = 0     ! Parallel back substitution (use 1 for serial)

Call pardiso(pt,maxfct,mnum,mtype,phase,nn,mat2,ia,jb,idum,nrhs,iparm,msglvl,q,fi,error)
If (error/=0) Then
  Write(*,*) "coup_dsolve: Error2!", error
  Stop
End If
variable = fi
phase=-1
Call pardiso(pt,maxfct,mnum,mtype,phase,nn,ddum,idum,idum,idum,nrhs,iparm,msglvl,ddum,ddum,error)
If (error/=0) Then
  Write(*,*) "coup_dsolve: Error3!", error
  Stop
End If
Nullify(jb)
Nullify(mat2)

End Subroutine coup_dsolve