Subroutine poiss

Use comm
Implicit None
Real(Kind=8) :: normalizer
Integer :: i,j,k,l,ii

pc = -1d8
Call fill
      
Do j=1,ny
  Do k=1,nz
    ap(ijk(j,k)+nn,8) = ap(ijk(j,k)+nn,8) + orne(j,k,1)
    af(ijk(j,k)+nn) = af(ijk(j,k)+nn) + orne(j,k,2)
    ap(ijk(j,k)+2*nn,13) = ap(ijk(j,k)+2*nn,13) + ornp(j,k,1)
    af(ijk(j,k)+2*nn) = af(ijk(j,k)+2*nn) + ornp(j,k,2)
    ap(ijk(j,k)+3*nn,18) = ap(ijk(j,k)+3*nn,18) + ornn(j,k,1)
    af(ijk(j,k)+3*nn) = af(ijk(j,k)+3*nn) + ornn(j,k,2)
  End Do
End Do
Do i=1,nnn
  normalizer = Sqrt(Sum(ap(i,1:20)**2))
  ap(i,1:20) = ap(i,1:20)/normalizer
  af(i) = af(i)/normalizer
End Do

Call coup_dsolve(ny,nz,nv,temp_sol,ap,af,connect,11,directCGS)

Do k=1,nz
  Do j=1,ny
    phi(j,k)=temp_sol(ijk(j,k))
    rnp(j,k)=temp_sol(ijk(j,k)+2*nn)
    rnn(j,k)=temp_sol(ijk(j,k)+3*nn)
    rne(j,k)=temp_sol(ijk(j,k)+nn)
  End Do
End Do

Call fill
ForAll(j=1:ny)
  rne(j,0) = rne(j,1)
  rnn(j,0) = rnn(j,1)
  rnp(j,0) = rnp(j,1)
  rne(j,nz+1) = rne(j,nz)
  rnn(j,nz+1) = rnn(j,nz)
  rnp(j,nz+1) = rnp(j,nz)
End ForAll
ForAll(k=1:nz)
  rne(0,k) = rne(1,k)
  rnn(0,k) = rnn(1,k)
  rnp(0,k) = rnp(1,k)
  phi(0,k) = phi(1,k)
  rne(ny+1,k) = rne(ny,k)
  rnn(ny+1,k) = rnn(ny,k)
  rnp(ny+1,k) = rnp(ny,k)
  phi(ny+1,k) = phi(ny,k)
End ForAll
Do j=1,ny
  Do k=1,nz
    ap(ijk(j,k)+nn,8) = ap(ijk(j,k)+nn,8) + orne(j,k,1)
    af(ijk(j,k)+nn) = af(ijk(j,k)+nn) + orne(j,k,2)
    ap(ijk(j,k)+2*nn,13) = ap(ijk(j,k)+2*nn,13) + ornp(j,k,1)
    af(ijk(j,k)+2*nn) = af(ijk(j,k)+2*nn) + ornp(j,k,2)
    ap(ijk(j,k)+3*nn,18) = ap(ijk(j,k)+3*nn,18) + ornn(j,k,1)
    af(ijk(j,k)+3*nn) = af(ijk(j,k)+3*nn) + ornn(j,k,2)
  End Do
End Do
Do j=0,ny
  Do k=0,nz
    cdy(j,k)=-(2d0*sigy(j,k)/(dy(j+1)+dy(j)))*(phi(j+1,k)-phi(j,k))
    cdz(j,k)=-(2d0*sigz(j,k)/(dz(k+1)+dz(k)))*(phi(j,k+1)-phi(j,k))
  End Do
End Do
Do i=1,nnn
  normalizer = Sqrt(Sum(ap(i,1:20)**2))
  ap(i,1:20) = ap(i,1:20)/normalizer
  af(i) = af(i)/normalizer
End Do
ForAll(k=0:nz)
cdy(ny+1,k) = cdy(ny,k)
cdz(ny+1,k) = cdz(ny,k)
ey(ny+1,k) = ey(ny,k)
ez(ny+1,k) = ez(ny,k)
End ForAll
ForAll(j=0:ny+1)
cdy(j,nz+1) = cdy(j,nz)
cdz(j,nz+1) = cdz(j,nz)
ey(j,nz+1) = ey(j,nz)
ez(j,nz+1) = ez(j,nz)
End ForAll
If (MinVal(rne,Mask=ian==0)<0d0.OR.MinVal(rnn)<0d0.OR.MinVal(rnp)<0d0) Then
  Stop
End If

res = 0d0
rel = 0d0
res_calc=0d0
Do k=1,nz
  Do j=1,ny
    Do l=1,nv
      res_calc(j,k,l) = af(ijk(j,k)+(l-1)*nn)
      Do ii=1,20
        If (connect(ijk(j,k)+(l-1)*nn,ii)>0) res_calc(j,k,l) = res_calc(j,k,l) - ap(ijk(j,k)+(l-1)*nn,ii)*temp_sol(connect(ijk(j,k)+(l-1)*nn,ii))
      End Do
      res_calc(j,k,l) = res_calc(j,k,l)**2
    End Do
  End Do      
End Do
res=Sum(res_calc)
rel=res/Sum(af**2)
Print *, "poiss: res and rel: ", res, rel

End Subroutine poiss

