INTEGER JMAX,JM1,JM2,JMAX2,NSTEP
      PARAMETER(JMAX=20, JM1=JMAX-1, JM2=JMAX-2, JMAX2=2*JM2*JM2
     *          ,JM3=JMAX-3, NSTEP=4)
      INTEGER i,j,k,IA,NN
      DOUBLE PRECISION
     *      a(JMAX,JMAX,2,2),b(JMAX,JMAX,2,2),c(JMAX,JMAX,2,2)
     *     ,d(JMAX,JMAX,2,2),e(JMAX,JMAX,2,2),f(JMAX,JMAX,2)
     *     ,s(JMAX,JMAX),w(JMAX,JMAX),u(JMAX,JMAX),v(JMAX,JMAX)
     *     ,sbx(2,JMAX),sby(JMAX,2)
     *     ,vbx(2,JMAX),vby(JMAX,2)
     *     ,del,alf,del2,alf2
C     Specification statements for routine linbcg below
C     bA & xA are the right side & unknown of linear system
C     If you change NMAX here, you must also change it in the subroutines
      INTEGER NP,NMAX,ITOL,ITMAX
      DOUBLE PRECISION TOL
      PARAMETER (NP= JMAX2,NMAX=50000,ITOL=3,TOL=1.d-9,ITMAX=75)
      INTEGER iter,ija
      DOUBLE PRECISION bA(NP),xA(NP),bcmp(NP),sa,err
      COMMON /mat/ sa(NMAX),ija(NMAX)
      COMMON i,j,k,IA,NN,iter
     * ,a,b,c,d,e,f,s,w,u,v,sbx,sby,vbx,vby,del,alf,del2,alf2,
     * ,bA,xA,bcmp,err