SUBROUTINE press3
      INCLUDE 'include3.inc'
C     Calculate pressure and stress distributions
!     Calculate boundary values of pressure
      if(icase.eq.1)then					 ! Cartesian coordiantes
!     x boundaries,left 
      p(1,1)=0.d0
      do j=2,JMAX
	  p(1,j)=p(1,j-1)+(0.5d0/(alf*Re*del))*
     *	  (v(1,j)-2.d0*v(2,j)+v(3,j)
     *      +v(1,j-1)-2.d0*v(2,j-1)+v(3,j-1))
      enddo
!     y boundary, top 
      do i=2,JMAX
	  p(i,JMAX)=p(i-1,JMAX)+(0.5d0*alf2/(Re*del))*
     *      (u(i,JMAX)-2.d0*u(i,JM1)+u(i,JM2)
     *      +u(i-1,JMAX)-2.d0*u(i-1,JM1)+u(i-1,JM2))
      enddo
!     x boundary, right
      do j=JM1,1,-1
	  p(JMAX,j)=p(JMAX,j+1)-(0.5d0/(alf*Re*del))*
     *	  (v(JMAX,j)-2.d0*v(JM1,j)+v(JM2,j)
     *      +v(JMAX,j+1)-2.d0*v(JM1,j+1)+v(JM2,j+1))
      enddo
!     y boundary, bottom
      do i=JM1,1,-1
	  p(i,1)=p(i+1,1)-(0.5d0*alf2/(Re*del))*
     *      (u(i,1)-2.d0*u(i,2)+u(i,3)
     *      +u(i+1,1)-2.d0*u(i+1,2)+u(i+1,3))
      enddo
	else							  ! Cylindrical polar coordinates
      p(1,1)=0.d0
!     r or z boundary, r=1 or z=0
      do j=2,JMAX
	  p(1,j)=p(1,j-1)+(1.d0*gam2/(alf*Re*del*r(1)))*
     *	  ((r(3)*v(3,j)-r(2)*v(2,j))/(r(3)+r(2))
     *      -(r(2)*v(2,j)-r(1)*v(1,j))/(r(2)+r(1))
     *	  +(r(3)*v(3,j-1)-r(2)*v(2,j-1))/(r(3)+r(2))
     *      -(r(2)*v(2,j-1)-r(1)*v(1,j-1))/(r(2)+r(1)))
      enddo
!     theta boundary, theta=1 
      do i=2,JMAX
	  p(i,JMAX)=p(i-1,JMAX)+(0.5d0*alf2/(gamma*Re*del))*
     *      ((u(i,JMAX)  -2.d0*u(i,JM1)  +u(i,JM2))/r(i)
     *      +(u(i-1,JMAX)-2.d0*u(i-1,JM1)+u(i-1,JM2))/r(i-1))
	*      +(0.5d0*del/gamma)*(v(i,JMAX)**2+v(i-1,JMAX)**2)
     *      -0.5d0*(u(i,JMAX)**2-u(i-1,JMAX)**2)
      enddo
!     r or z boundary, r=beta or z=1
      do j=JM1,1,-1
	  p(JMAX,j)=p(JMAX,j+1)-(1.d0*gam2/(alf*Re*del*r(JMAX)))*
     *	  ((r(JMAX)*v(JMAX,j)-r(JM1)*v(JM1,j))/(r(JMAX)+r(JM1))
     *      -(r(JM1)*v(JM1,j)-r(JM2)*v(JM2,j))/(r(JM1)+r(JM2))
     *	  +(r(JMAX)*v(JMAX,j-1)-r(JM1)*v(JM1,j-1))/(r(JMAX)+r(JM1))
     *      -(r(JM1)*v(JM1,j-1)-r(JM2)*v(JM2,j-1))/(r(JM1)+r(JM2)))
      enddo
!     theta boundary, theta=0
      do i=JM1,1,-1
	  p(i,1)=p(i+1,1)-(0.5d0*alf2/(gamma*Re*del))*
     *      ((u(i,1)  -2.d0*u(i,2)  +u(i,3))/r(i)
     *      +(u(i+1,1)-2.d0*u(i+1,2)+u(i+1,3))/r(i+1))
	*      -(0.5d0*del/gamma)*(v(i,1)**2+v(i+1,1)**2) 
     *      +0.5d0*(u(i+1,1)**2-u(i,1)**2)
      enddo
	endif
!  
!     Calculate coefficients
      if(icase.eq.1)then            ! Cartesian coordinates
      sumf=0.
      do i=2,JM1
        do j=2,JM1
         a(i,j,1,1)=1.d0
         b(i,j,1,1)=1.d0
         c(i,j,1,1)=alf2
         d(i,j,1,1)=alf2
         e(i,j,1,1)=-2.d0*(1.d0+alf2)-eps
         f(i,j,1)=(2.d0*alf2/del2)*((s(i+1,j)-2.d0*s(i,j)+s(i-1,j))
     *              *(s(i,j+1)-2.d0*s(i,j)+s(i,j-1))
	*         -((s(i+1,j+1)-s(i-1,j+1)-s(i+1,j-1)+s(i-1,j-1))/4.d0)**2)
         sumf=sumf+f(i,j,1) 
	   enddo
	enddo
	else                         ! Cylindrical coordinates
      sumf=0.
      do i=2,JM1
        do j=2,JM1
         a(i,j,1,1)=gam2/(r(i)**2)
         b(i,j,1,1)=gam2/(r(i)**2)
         c(i,j,1,1)=alf2/(r(i)**2)
         d(i,j,1,1)=alf2/(r(i)**2)
         e(i,j,1,1)=-a(i,j,1,1)-b(i,j,1,1)-c(i,j,1,1)-d(i,j,1,1)-eps
         f(i,j,1)=(2.d0*alf*gamma/(r(i)**2))*
     *            (0.25d0*(u(i+1,j)-u(i-1,j))*(v(i,j+1)-v(i,j-1))
	*            -0.25d0*(v(i+1,j)-v(i-1,j))*(u(i,j+1)-u(i,j-1)))
	*        +(gamma*del/(r(i)**2))
     *        *0.5d0*(u(i+1,j)**2+v(i+1,j)**2-u(i-1,j)**2-v(i-1,j)**2)
         sumf=sumf+f(i,j,1) 
	   enddo
	enddo
	endif
!     x or r boundaries,i.e. left and right
      do i=2,JM1,JM3
        do j=2,JM1
          if(i.eq.2)then
            f(i,j,1)=f(i,j,1)-b(i,j,1,1)*p(1,j)
            b(i,j,1,1)=0.d0
          endif
          if(i.eq.JM1)then
            f(i,j,1)=f(i,j,1)-a(i,j,1,1)*p(JMAX,j)
            a(i,j,1,1)=0.d0
          endif
        enddo
      enddo
!     y or theta boundaries, i.e., top and bottom
      do j=2,JM1,JM3
        do i=2,JM1
          if(j.eq.2)then
            f(i,j,1)=f(i,j,1)-d(i,j,1,1)*p(i,1)
            d(i,j,1,1)=0.d0
         endif
         if(j.eq.JM1)then
           f(i,j,1)=f(i,j,1)-c(i,j,1,1)*p(i,JMAX)
           c(i,j,1,1)=0.d0
         endif
        enddo
      enddo
!     Invert matrix. There is only one dependent variable here
      NN=JM2*JM2
      j=2
      i=1
      do IA=1,NN,1
        i=i+1
        if(i.gt.JM1)then
          i=2
          j=j+1
        endif
        sa(IA)=e(i,j,1,1)     !  Store diagonal elements
      enddo
      ija(1)=NN+2            !  Store size of matrix + 2
      k=NN+1
      j=2
      i=1
      do IA=1,NN,1        !  IA is sequence number 
                            !  i.e., It is row counter for coefficient matrix
        i=i+1
        if(i.gt.JM1)then
          i=2
          j=j+1
        endif
        if(j.gt.2)then
            k=k+1
            sa(k)=d(i,j,1,1)  !  Store j-1 coefficient
            ija(k)=IA-JM2
        endif
        if(i.gt.2)then
            k=k+1
            sa(k)=b(i,j,1,1)  !  Store i-1 coefficient
            ija(k)=IA-1
        endif
        if(i.lt.JM1)then
            k=k+1
            sa(k)=a(i,j,1,1)  !  Store i+1 coefficient
            ija(k)=IA+1
            endif
        if(j.lt.JM1)then
            k=k+1
            sa(k)=c(i,j,1,1)  !  Store j+1 coefficient
            ija(k)=IA+JM2
        endif
            ija(IA+1)=k+1  !  As each row is completed, store index to next
      enddo                !  End of loop on rows
!     Initialize for PBCG
      k =0
      do j=2,JM1
        do i=2,JM1
            k=k+1
			m=1
            xA(k)=p(i,j)
            bA(k)=f(i,j,m)
        enddo
      enddo
      call linbcg(NP1,bA,xA,ITOL,TOL,ITMAX,iter,err)
      write(*,'(/1x,a,e15.6)') 'Estimated error:',err
      write(*,'(/1x,a,i6)') 'Iterations needed:',iter
      call dsprsax(sa,ija,x,bcmp,NP1)
C     this is a double precision version of sprsax
      j=2
      i=1
      do IA=1,NN,1
        i=i+1
        if(i.gt.JM1)then
          i=2
          j=j+1
        endif
        p(i,j)=xA(IA)
      enddo

!     Calculate shear stress
	if(icase.eq.1)then		! Cartesian coordinates 
      do i=2,JM1
	  do j=2,JM1
	tau(i,j)=
	* -(alf*(u(i,j+1)-u(i,j-1))+(v(i+1,j)-v(i-1,j)))/(2.d0*del*Re)
	  enddo
	enddo
!       Store boundary values of shear stress into tau(i,j) array
      do j=2,JM1
        tau(1,j)=-(alf*(u(1,j+1)-u(1,j-1))/(2.d0*del*Re)
     *             +(v(2,j)-v(1,j))/(del*Re))
        tau(JMAX,j)=-(alf*(u(JMAX,j+1)-u(JMAX,j-1))/(2.d0*del*Re)
     *             +(v(JMAX,j)-v(JM1,j))/(del*Re))
      enddo
      do i=2,JM1
        tau(i,1)=-(alf*(u(i,2)-u(i,1))/(del*Re)
     *             +(v(i+1,1)-v(i-1,1))/(2.d0*del*Re))
        tau(i,JMAX)=-(alf*(u(i,JMAX)-u(i,JM1))/(del*Re)
     *             +(v(i+1,JMAX)-v(i-1,JMAX))/(2.d0*del*Re))
      enddo
	tau(1,1)=-(alf*(u(1,2)-u(1,1))+v(2,1)-v(1,1))/(del*Re)
	tau(JMAX,1)=
	* -(alf*(u(JMAX,2)-u(JMAX,1))+v(JMAX,1)-v(JM1,1))/(del*Re)
	tau(JMAX,JMAX)=
	* -(alf*(u(JMAX,JMAX)-u(JMAX,JM1))+v(JMAX,JMAX)-v(JM1,JMAX))
     *  /(del*Re)
	tau(1,JMAX)=
	* -(alf*(u(1,JMAX)-u(1,JM1))+v(2,JMAX)-v(1,JMAX))/(del*Re)
	else		! *******************************            Cylindrical polar coordinates
      do i=2,JM1
	  do j=2,JM1
	tau(i,j)= -(alf*(u(i,j+1)-u(i,j-1))/r(i)
     *         +gamma*(v(i+1,j)/r(i+1)-v(i-1,j)/r(i-1)))/(2.d0*del*Re)
	  enddo
	enddo
!       Store boundary values of shear stress into tau(i,j) array
      do j=2,JM1		   ! r boundaries
       tau(1,j)=-(alf*(u(1,j+1)-u(1,j-1))/(2.d0*del*Re*r(1))
     *             +gamma*(v(2,j)/r(2)-v(1,j)/r(1))/(del*Re))
       tau(JMAX,j)=-(alf*(u(JMAX,j+1)-u(JMAX,j-1))/(2.d0*del*Re*r(JMAX))
     *             +gamma*(v(JMAX,j)/r(JMAX)-v(JM1,j)/r(JM1))/(del*Re))
      enddo
      do i=2,JM1		   ! theta boundaries
        tau(i,1)=-(alf*(u(i,2)-u(i,1))/(del*Re*r(i))
     *           +gamma*(v(i+1,1)/r(i+1)-v(i-1,1)/r(i-1))/(2.d0*del*Re))
        tau(i,JMAX)=-(alf*(u(i,JMAX)-u(i,JM1))/(del*Re*r(i))
     *     +gamma*(v(i+1,JMAX)/r(i+1)-v(i-1,JMAX)/r(i-1))/(2.d0*del*Re))
      enddo
	tau(1,1)=-(alf*(u(1,2)-u(1,1))/r(1)
	*         +gamma*(v(2,1)/r(2)-v(1,1)/r(1)))/(del*Re)
	tau(JMAX,1)=
	* -(alf*(u(JMAX,2)-u(JMAX,1))/r(JMAX)
     * +gamma*(v(JMAX,1)/r(JMAX)-v(JM1,1)/r(JM1)))/(del*Re)
	tau(JMAX,JMAX)=
	* -(alf*(u(JMAX,JMAX)-u(JMAX,JM1))/r(JMAX)
     * +gamma*(v(JMAX,JMAX)/r(JMAX)-v(JM1,JMAX)/r(JM1)))
     *  /(del*Re)
	tau(1,JMAX)=
	* -(alf*(u(1,JMAX)-u(1,JM1))/r(1)
     * +gamma*(v(2,JMAX)/r(2)-v(1,JMAX)/r(1)))/(del*Re)
 	endif
      RETURN
	END