cccccccccccccccccccccccccccccccccccccccccccc

c   Pressure calculation

cccccccccccccccccccccccccccccccccccccccccccc

    subroutine pressure()

    use mvrs

    implicit none


    F = 0.0

    do i = 2,N-1
    do j = 2,N-1

    F(i,j) = (u(i+1,j+1)-u(i+1,j-1))/(4.0*dL**2.0) -
     c  (u(i-1,j+1)-u(i-1,j-1))/(4.0*dL**2.0) +
     c  (psi(i+1,j)+psi(i-1,j)-2.0*psi(i,j))*
     c  (psi(i,j+1)+psi(i,j-1)-2.0*psi(i,j))/dL**4.0

    enddo
    enddo

    !Calculate dw/dx and dw/dy
    do i = 2, N-1
    do j = 2, N-1
     dwx(i,j) = (w(i+1,j)-w(i-1,j))/(2.0*dL)
     dwy(i,j) = (w(i,j+1)-w(i,j-1))/(2.0*dL)
    enddo
    enddo
    do j = 2, N-1
     dwy(1,j) = (w(1,j+1)-w(1,j-1))/(2.0*dL)
     dwx(j,1) = (w(j+1,N)-w(j-1,N))/(2.0*dL)
    enddo

    dwy(1,1) = (w(1,2)-w(1,1))/dL
    dwy(1,N) = (w(1,N)-w(1,N-1))/dL
    dwy(N,1) = (w(N,2)-w(N,1))/dL
    dwy(N,N) = (w(N,N)-w(N,N-1))/dL

    dwx(1,1) = (w(2,1)-w(1,1))/dL
    dwx(1,N) = (w(2,N)-w(1,N))/dL
    dwx(N,1) = (w(N,1)-w(N-1,1))/dL
    dwx(N,N) = (w(N,N)-w(N-1,N))/dL

    pc = 0

    rp = 1.0


    do while (rp > maxres)

    pc = pc+1
    pn = p



    !----------------------------- Start of matrix solution for current iteration

    do i = 1, N

c   -------------------------

c   Build tridiagonal vectors

c   -------------------------

    ud = 1.0

    ld = 1.0

    ud(1) = 2.0

    ld(N) = 2.0

    md = -4.0



c   Arrange dependent vectors (RHS of EQ)
    if (i == 1) then

        pk = -2.0*(pn(i+1,:)+mu*dL*dwy(i,:))

    else if (i == N) then

        pk = -2.0*(pn(i-1,:)-mu*dL*dwy(i,:))

    else

        pk = 2.0*rho*dL**2.0*F(i,j) -

     c  (pn(i+1,:)+pn(i-1,:))

    end if


    pk(1) = pk(1) - 2.0*mu*dL*dwx(i,1)
    pk(N) = pk(N) + 2.0*mu*dL*dwx(i,N)



c   ---------------------------------------

c   Solve tridagonal system with ld,md,ud,p

c   ---------------------------------------

    do z = 2, N

        ld(z) = ld(z)/md(z-1)

        md(z) = md(z) - ld(z)*ud(z-1)

        pk(z) = pk(z) - ld(z)*psk(z-1)

    enddo



    pni(N) = pk(N)/md(N)



    do z = N-1,1,-1

        pni(z) = (pk(z) - ud(z)*pni(z+1))/md(z)

    enddo



    pn(i,:) = pni(:)


    enddo



    rp = sqrt(sum((pn - p)**2.0))



    if (pc == 1) then

        rpo = rp

    end if



    rp = rp/rpo

    p = pn



c   write(*,'(A, I5, A, e10.3)'), 'Iteration #:', cc,

c     c '     Residual:  ', rs



    enddo





    end subroutine