c***********************  problem name: domains ************************
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine a1xy(x,y,u,ux,uy,rl,itag,values)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            real(kind=rknd), dimension(*) :: values
            character(len=80) :: su
            common /val0/k0,ku,kx,ky,kl
            common /atest2/iu(100),a1x,a1y,a1u,a2x,a2y,a2u,
     +          bux,buy,cu0,cu1,ru(90),su(100)
cy
        values(k0)=a1x*ux+a1y*uy+a1u*u
        values(ku)=a1u
        values(kx)=a1x
        values(ky)=a1y
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine a2xy(x,y,u,ux,uy,rl,itag,values)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            real(kind=rknd), dimension(*) :: values
            character(len=80) :: su
            common /val0/k0,ku,kx,ky,kl
            common /atest2/iu(100),a1x,a1y,a1u,a2x,a2y,a2u,
     +          bux,buy,cu0,cu1,ru(90),su(100)
cy
        values(k0)=a2x*ux+a2y*uy+a2u*u
        values(ku)=a2u
        values(kx)=a2x
        values(ky)=a2y
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine fxy(x,y,u,ux,uy,rl,itag,values)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            real(kind=rknd), dimension(*) :: values
            character(len=80) :: su
            common /val0/k0,ku,kx,ky,kl
            common /atest2/iu(100),a1x,a1y,a1u,a2x,a2y,a2u,
     +          bux,buy,cu0,cu1,ru(90),su(100)
cy
        values(k0)= - bux*ux - buy*uy - cu0 - cu1*u
        values(ku)= - cu1
        values(kx) = - bux
        values(ky) = - buy
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gnxy(x,y,u,rl,itag,values)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            real(kind=rknd), dimension(*) :: values
            common /val1/k0,ku,kl
cy
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gdxy(x,y,rl,itag,values)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            real(kind=rknd), dimension(*) :: values
            common /val2/k0,kl,klb,kub,kic,kim,kil
cy
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine p1xy(x,y,u,ux,uy,rl,itag,values)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            real(kind=rknd), dimension(*) :: values
            common /val0/k0,ku,kx,ky,kl
cy
        values(k0)=1.0e0_rknd
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine p2xy(x,y,dx,dy,u,ux,uy,rl,itag,jtag,values)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            real(kind=rknd), dimension(*) :: values
            common /val0/k0,ku,kx,ky,kl
cy
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine qxy(x,y,u,ux,uy,rl,itag,values)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            real(kind=rknd), dimension(*) :: values
            common /val3/kf,kf1,kf2,kad
cy
        s=x**2+y**2
        values(kf)=real(itag,rknd)
        values(kad)=s
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine sxy(rl,s,itag,values)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            real(kind=rknd), dimension(2,*) :: values
            common /val4/j0,js,jl
cy
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine usrcmd(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            character(len=80), dimension(100) :: sp,su
cy
        iu1=iu(1)
c
        call usrset(iu,ru,su)
c
        ip(41)=0
        if(iu1/=iu(1)) ip(41)=-1
c
        ispd=ip(8)
        a1x=ru(1)
        a1y=ru(2)
        a1u=ru(3)
        a2x=ru(4)
        a2y=ru(5)
        a2u=ru(6)
        bux=ru(7)
        buy=ru(8)
        cu0=ru(9)
        cu1=ru(10)
        if(ispd==1) then
            if((a1y/=a2x).or.(a1u/=0.0e0_rknd).or.
     +      (a2u/=0.0e0_rknd).or.(bux/=0.0e0_rknd).or.
     1      (buy/=0.0e0_rknd)) ip(8)=0
        else
            if((a1y==a2x).and.(a1u==0.0e0_rknd).and.
     +      (a2u==0.0e0_rknd).and.(bux==0.0e0_rknd).and.
     1      (buy==0.0e0_rknd)) ip(8)=1
        endif
        if(ispd/=ip(8)) ip(5)=max(ip(5),1)
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine usrfil(file,len)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), save :: len0
            character(len=80), save, dimension(100) :: file0
            character(len=80), dimension(*) :: file
cy
            data (file0(i),i=  1, 10)/
     +  'n i= 1,n=domain,a= d,t=i',
     1  'n i= 1,n=a1x,   a=x1,t=r',
     2  'n i= 2,n=a1y,   a=y1,t=r',
     3  'n i= 3,n=a1u,   a=u1,t=r',
     4  'n i= 4,n=a2x,   a=x2,t=r',
     5  'n i= 5,n=a2y,   a=y2,t=r',
     6  'n i= 6,n=a2u,   a=u2,t=r',
     7  'n i= 7,n=bux,   a=bx,t=r',
     8  'n i= 8,n=buy,   a=by,t=r',
     9  'n i= 9,n=cu0,   a=c0,t=r'/
            data (file0(i),i= 11, 20)/
     +  'n i=10,n=cu1,   a=c1,t=r',
     1  's n=domain,v=1,l="texas"',
     2  's n=domain,v=2,l="doughnut"',
     3  's n=domain,v=3,l="cmos device"',
     4  's n=domain,v=4,l="lake superior"',
     5  's n=domain,v=5,l="dd20 logo"',
     6  's n=domain,v=6,l="at&t logo"',
     7  's n=domain,v=7,l="north sea"',
     8  's n=domain,v=8,l="airfoil"',
     9  's n=domain,v=9,l="planter"'/
            data (file0(i),i= 21, 30)/
     +  's n=domain,v=10,l="fan"',
     1  's n=domain,v=11,l="arc"',
     2  's n=domain,v=12,l="spiral"',
     3  's n=domain,v=13,l="ucsd logo"',
     4  's n=domain,v=14,l="nozzle"',
     5  's n=domain,v=15,l="pltmg"',
     6  's n=domain,v=16,l="crack"',
     7  's n=domain,v=17,l="monterey bay"',
     9  's n=domain,v=18,l="siam logo"',
     +  's n=domain,v=19,l="mexico"'/
            data (file0(i),i= 31, 34)/
     +  's n=domain,v=20,l="ellipse"',
     1  's n=domain,v=21,l="hole"',
     2  's n=domain,v=22,l="british columbia"',
     3  's n=domain,v=23,l="israel"'/
c
            data len0/34/
c
        len=len0
        do i=1,len
            file(i)=file0(i)
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine box(nvf,vx,vy,nbf,ibndry,sf)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(7,*) :: ibndry
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
cy
c       rescale domain to lie inside the unit square
c
        xmin=vx(1)
        xmax=xmin
        ymin=vy(1)
        ymax=ymin
        do i=1,nvf
            xmin=min(xmin,vx(i))
            xmax=max(xmax,vx(i))
            ymin=min(ymin,vy(i))
            ymax=max(ymax,vy(i))
        enddo
c
c       compute scaled coordinates
c
        xmid=(xmin+xmax)/2.0e0_rknd
        ymid=(ymin+ymax)/2.0e0_rknd
        scale=1.0e0_rknd/max(xmax-xmin,ymax-ymin)
        do i=1,nvf
            vx(i)=0.5e0_rknd+scale*(vx(i)-xmid)
            vy(i)=0.5e0_rknd+scale*(vy(i)-ymid)
        enddo
        do i=1,nbf
            if(ibndry(3,i)==0) cycle
            sf(1,i)=0.5e0_rknd+scale*(sf(1,i)-xmid)
            sf(2,i)=0.5e0_rknd+scale*(sf(2,i)-ymid)
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gdata(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), save :: iprob,ifirst
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            character(len=80), dimension(100) :: sp,su
cy
            external sxy
            data iprob,ifirst/1,1/
c
        if(ip(41)==1) then
            sp(4)='domains'
            sp(5)='domains_figxxx.ext'
            sp(6)='domains_mpixxx.rw'
            sp(7)='domains.jnl'
            sp(9)='domains_mpixxx.out'
c
            iu(1)=15
c
c       initialize as laplacian with dirichlet b.c.
c
            do i=1,10
                ru(i)=0.0e0_rknd
            enddo
            ru(1)=1.0e0_rknd
            ru(5)=1.0e0_rknd
            ru(9)=1.0e0_rknd
            ip(8)=1
        endif
c
        ip(5)=max(ip(5),ifirst)
        ip(6)=iprob
        ip(20)=5
c
        if(iu(1)==1) call gd1(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==2) call gd2(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==3) call gd3(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==4) call gd4(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==5) call gd5(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==6) call gd6(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==7) call gd7(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==8) call gd8(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==9) call gd9(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==10) call gd10(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==11) call gd11(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==12) call gd12(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==13) call gd13(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==14) call gd14(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==15) call gd15(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==16) call gd16(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==17) call gd17(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==18) call gd18(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==19) call gd19(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==20) call gd20(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==21) call gd21(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==22) call gd22(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
        if(iu(1)==23) call gd23(vx,vy,sf,itnode,
     +      ibndry,ip,rp,sp,iu,ru,su,sxy)
c
        nvf=ip(2)
        nbf=ip(3)
        call box(nvf,vx,vy,nbf,ibndry,sf)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd1(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), dimension(12) :: ix,iy
            integer(kind=iknd), dimension(5) :: ixc,iyc
            integer(kind=iknd), save :: ntf,nvf,nbf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), save :: x,y,xc,yc,hmax,grade
            character(len=80), dimension(100) :: sp,su
cy
c       domain in the shape of texas.
c       dirichlet boundary conditions on all boundary edges
c
            external sxy
            data ix/2,50,50,90,90,158,159,160,130,90,62,40/
            data iy/97,97,160,160,125,121,102,62,0,32,60,60/
            data ixc/51,169,106,133,82/
            data iyc/54,82,12,42,50/
            data ntf,nvf,nbf/1,12,12/
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
c
        sp(2)='texas'
        sp(1)='texas'
        sp(3)='texas'
c
        itnode(1,1)=1
        itnode(2,1)=12
        itnode(3,1)=0
        itnode(4,1)=0
        itnode(5,1)=1
        do i=1,nvf
            vx(i)=real(ix(i),rknd)
            vy(i)=real(iy(i),rknd)
        enddo
        do i=1,nbf
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
            ibndry(1,i)=i
            ibndry(2,i)=i+1
            ibndry(3,i)=0
            ibndry(4,i)=2
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=0
        enddo
        ibndry(2,nbf)=1
        ibndry(3,7)=2
        ibndry(3,8)=4
        ibndry(3,9)=3
        ibndry(3,10)=5
        ibndry(3,11)=1
        do i=1,nbf
            if(ibndry(3,i)>0) then
                i1=ibndry(1,i)
                i2=ibndry(2,i)
                i3=ibndry(3,i)
                xc=real(ixc(i3),rknd)
                yc=real(iyc(i3),rknd)
                call centre(vx(i1),vy(i1),vx(i2),vy(i2),
     +              xc,yc,sf(1,i),sf(2,i))
           endif
        enddo
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
c
c       refine long curved edges
c
        call sklutl(1_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd2(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), save :: ntf,nvf,nbf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), save :: hmax,grade
            character(len=80), dimension(100) :: sp,su
cy
c       annular region between circles of radius one and six
c       dirichlet boundary conditions on inner circle
c       neumann boundary conditions on outer circle
c       there are two irregular octagons concentric with the
c       circles. there are 24 regions, and three groups of
c       eight regions are similar.
c
            external sxy
            data ntf,nvf,nbf/24,32,56/
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
c
        sp(2)='doughnut'
        sp(1)='doughnut'
        sp(3)='doughnut'
c
        pi=3.141592653589793e0_rknd
        pi4=pi/4.0e0_rknd
c
        do i=1,8
            r2=2.0e0_rknd
            if((i/2)*2==i) r2=4.0e0_rknd
            r3=3.0e0_rknd
            if((i/2)*2==i) r3=5.0e0_rknd
            r4=6.0e0_rknd
            ang=real(i-1,rknd)*pi4
            k=(i-1)*4+1
            vx(k)=cos(ang)
            vy(k)=sin(ang)
            vx(k+1)=r2*vx(k)
            vy(k+1)=r2*vy(k)
            vx(k+2)=r3*vx(k)
            vy(k+2)=r3*vy(k)
            vx(k+3)=r4*vx(k)
            vy(k+3)=r4*vy(k)
        enddo
        do i=1,nbf
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
c
        do i=1,32
            ibndry(1,i)=i
            ibndry(2,i)=i+4
            if(i+4>32) ibndry(2,i)=i+4-32
            ibndry(3,i)=0
            ibndry(4,i)=0
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=0
        enddo
        do i=1,8
            k=(i-1)*4+1
            ibndry(3,k)=1
            ibndry(4,k)=2
            ibndry(3,k+3)=1
            ibndry(4,k+3)=1
            do j=1,3
                k=32+(i-1)*3+j
                l=(i-1)*4+j
                ibndry(1,k)=l
                ibndry(2,k)=l+1
                ibndry(3,k)=0
                ibndry(4,k)=0
                ibndry(5,k)=0
                ibndry(6,k)=0
                ibndry(7,k)=0
            enddo
        enddo
c
c       comput itnode, find symmetry
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        call sklutl(2_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
c
c       label regions
c
        itnode(5,1)=3
        itnode(5,5)=2
        itnode(5,9)=1
        do i=1,ntf
            j=abs(itnode(3,i))
            if(j>0) itnode(5,i)=itnode(5,j)
        enddo
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd3(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), save, dimension(35) :: b1,b2,b4
            integer(kind=iknd), save, dimension(30) :: ix,iy
            integer(kind=iknd), save :: ntf,nvf,nbf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), save :: hmax,grade
            character(len=80), dimension(100) :: sp,su
cy
c       an irregular region composed of six subregions, two quite small,
c       which define a typical cmos device
c       the data was provided by wolfgang fichtner and donald j. rose
c
            external sxy
            data ntf,nvf,nbf/6,30,35/
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
            data (ix(i),i=   1,  30)/
     +          0,  180,  240,  255,  315,  405,  735,  780,
     1        825, 1185, 1260, 1275, 1320, 1500, 1500, 1260,
     2       1185,  825,  780,  735,  315,  255,  735, 1000,
     3          0,    0, 1000, 1500,    0, 1500/
            data (iy(i),i=   1,  30)/
     +        -11,  -12,  -31,   25,   37,   41,   40,   21,
     1         19,   14,    1,  -55,  -39,  -35,  -63,  -63,
     2        -81,  -83,  -80,  -58,  -56,  -39, -150, -150,
     3        -39, -750, -750, -750,-1000,-1000/
c
            data b1/ 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15,16,16,
     +           17,18,19,20,21,22,22,25,20,23,24,25,26,27,15,28,29,26/
            data b2/ 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15,16,12,17,
     +           18,19,20,21,22, 3,25, 1,23,24,27,26,27,28,28,30,30,29/
            data b4/ 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 0, 0, 0,
     +            0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 2, 1/
c
        sp(2)='cmos device'
        sp(1)='cmos device'
        sp(3)='cmos device'
c
        do i=1,nvf
            vx(i) =real(ix(i),rknd)*2.9563e-4
            vy(i) =real(iy(i),rknd)*2.9563e-4
        enddo
        do i=1,nbf
            ibndry(1,i)=b1(i)
            ibndry(2,i)=b2(i)
            ibndry(3,i)=0
            ibndry(4,i)=b4(i)
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=0
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
c
c       make itnode
c
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd4(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), save, dimension(303) :: ix,iy
            integer(kind=iknd), save :: ntf,nvf,nbf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), save :: hmax,grade
            character(len=80), dimension(100) :: sp,su
cy
c       region in the shape of lake superior
c       there are five subregions and six islands
c       the data was provided by r. bruce simpson
c
            external sxy
            data ntf,nvf,nbf/5,303,313/
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
            data (ix(i),i=   1,  50)/
     +     -891541481,-888470650,-872705746,-867803288,-827638721,
     1     -775894356,-733998108,-720648766,-698716116,-691908741,
     2     -665943623,-647050762,-639824581,-617760324,-596836376,
     3     -605006313,-627510357,-620834494,-635938740,-643797827,
     4     -643896866,-633085680,-592166805,-588535309,-553069878,
     5     -545711184,-527769947,-468826580,-443255520,-398677874,
     6     -362327838,-351568484,-311511421,-260318851,-242016149,
     7     -205767846,-194694924,-165644574,-143858552,-118571806,
     8      -89657599, -46563160,  -7156760,  32215631,  60890228,
     9       60928792,  39424559,  28667879,   7170380,   8968640/
            data (ix(i),i=  51, 100)/
     +       10765760, -14367931, -32358369, -46769229, -93919420,
     1     -104951930,-108823633,-101584518, -94241017, -61523801,
     2      -50619459, -28907421, -50666702, -61600220, -39797199,
     3      -18084030,  28943431,  54311031,  90588319,  94328558,
     4      105342555, 116419089, 131092036, 131172514, 142146897,
     5      193172145, 222260547, 233405209, 240735531, 244494963,
     6      259090638, 262498355, 273434877, 284545517, 293664694,
     7      292019439, 299319220, 331354308, 338635898, 399603295,
     8      407055664, 464906311, 490399790, 544955301, 598504639,
     9      616538048, 634761524, 663559580, 653090858, 646236992/
            data (ix(i),i= 101, 150)/
     +      648030853, 640949774, 641000891, 648826694, 670686483,
     1      685469055, 692969561, 736259174, 741950417, 748102570,
     2      766219997, 773394203, 769277287, 761990309, 757881975,
     3      746499538, 746040964, 766680956, 777094603, 772988129,
     4      767073679, 763567352, 762461424, 755321169, 752808619,
     5      753453922, 771314287, 786032724, 792287445, 773488998,
     6      737138224, 726983166, 718857956, 711526585, 704416895,
     7      693369579, 696761417, 706480742, 708971071, 729173613,
     8      732301521, 720742178, 709879112, 705625916, 683496714,
     9      654260111, 635381174, 664933491, 653851652, 670529938/
            data (ix(i),i= 151, 200)/
     +      670097637, 574411726, 556869459, 478860712, 446368170,
     1      420757627, 360757804, 338736439, 328043842, 317628479,
     2      303819609, 282868814, 285980701, 268455315, 254680014,
     3      244376707, 233904195, 219831228, 216162419, 209189844,
     4      202351665, 195406818, 157050800, 150021160, 125474477,
     5       73096502,  45220131,  17386571,  -6945330, -17374951,
     6      -20846450, -34761500, -45182380, -55581152, -52159631,
     7      -52211881, -48796061, -20926540, -20975170, -52576512,
     8      -56044650, -77162850, -94730818, -98336190,-115933979,
     9     -122879457,-115743721, -69940150, -59331119, -69731730/
            data (ix(i),i= 201, 250)/
     +      -97591662,-111755574,-111940694,-129516613,-154830980,
     1     -193663633,-186377645,-175828385,-172029650,-154273140,
     2     -154171550,-157623363,-199851763,-235065794,-256451440,
     3     -253518987,-275629210,-379806066,-433735085,-498349476,
     4     -592565298,-703822851,-745170450,-820339012,-872353649,
     5     -608656883,-608186007,-579297161,-568086433,-559666157,
     6     -566355705,-583187723,-583458996,-579201603,-587263441,
     7     -591520166,-602223587,-576770115,-572977018,-549145699,
     8     -544981623,-553187752,-553359890,-556836748,-573244572,
     9     -534449291,-531934547,-531668663,-529240274,-519116497/
            data (ix(i),i= 251, 300)/
     +     -521965599,-522571039,-527064848,-259048891,-237567210,
     1     -209000087,-184055698,-173409450,-162766552,-155639541,
     2     -130752015, -95321161, -88289100, -91850460,-120228684,
     3     -116768241,-152397847,-170117927,-191505325,-209371305,
     4     -205889106,-188140857,-216642499,-245249653,-252317262,
     5     -262850904, -17432991,   3484280,  13937110,  48804179,
     6       62737727,  73254800,  82016379,  66322219,  55831802,
     7       38390779,  22677940,   3490070,  -3493550,  -1747930,
     8      -10487590, -17459050, 438131142, 437850761, 455353165,
     9      483649635, 498189831, 508861971, 516141987, 512747717/
            data (ix(i),i= 301, 303)/
     +      498668671, 488142157, 459645033/
            data (iy(i),i=   1,  50)/
     +      166159201, 155386531, 143963063, 133102524, 136420846,
     1      155212080, 158664119, 174001110, 178389883, 167510176,
     2      182378805, 202851081, 202577615, 212360430, 190416360,
     3      164209175, 143823195, 127681434, 112335944,  96729130,
     4       94082898,  91031909, 116069555, 115944648,  88274997,
     5       90689111,  82183212, 106961095, 111578405, 152849865,
     6      157335603, 151815784, 161627734, 192534435, 208145809,
     7      218255091, 239298439, 260162497, 281135678, 296818018,
     8      328396988, 349393415, 365208673, 359945369, 349439454,
     9      338853741, 338789296, 341413832, 333451533, 322866678/
            data (iy(i),i=  51, 100)/
     +      317574883, 301698995, 285845375, 275293469, 206681108,
     1      174987423, 135313761, 132622755, 148459792, 174775326,
     2      190618157, 201153159, 174739587, 153603923, 180003631,
     3      185259604, 179981732, 166811454, 153733039, 132581246,
     4      111473274,  85079402,  69305450,  58719772,  53513449,
     5       54008722,  59652221,  43917981,  41371560,  33483499,
     6       33695781,  49629170,  49798739,  39389670,  39542049,
     7       28925619,  29050219,  71991080,  72132570, 113154662,
     8      105386281, 117447603, 115508866, 114505672, 142784142,
     9      146074021, 144090629, 150483525, 139480269, 128620267/
            data (iy(i),i= 101, 150)/
     +       80988193,  75423539,  70125639,  59818262,  60653287,
     1       55932897,  50931108,  63324392,  58265841,  42624259,
     2       46068209,  49038869,  59460950,  59138030,  69562429,
     3       79670852,  90255553, 117676401, 128750885, 139171314,
     4      149511147, 146702886, 130745316, 127777922, 143574929,
     5      170116031, 176216900, 171585441, 193089068, 208134031,
     6      211820865, 195477509, 216339254, 218681216, 215734005,
     7      220580411, 226020551, 252929902, 279541850, 312215924,
     8      322956347, 341012549, 343199420, 358924270, 373907781,
     9      388636851, 414405775, 479170609, 489328432, 516503143/
            data (iy(i),i= 151, 200)/
     +      527087212, 523458147, 517555094, 515094805, 535365105,
     1      566463518, 718588781, 765763617, 776130056, 773271894,
     2      765061188, 764676332, 785909224, 790898466, 780084085,
     3      769334936, 769178009, 776918602, 790102768, 790009642,
     4      779334068, 776600647, 773530149, 778755760, 794423151,
     5      815268230, 825750160, 830988121, 852151489, 841573811,
     6      844224453, 836309528, 838982296, 846954727, 831064415,
     7      815185833, 794003391, 783356619, 746306610, 704035950,
     8      714633512, 693549252, 688350058, 672492552, 667316341,
     9      677953959, 693780470, 741152763, 772867489, 788788366/
            data (iy(i),i= 201, 250)/
     +      794227743, 762562370, 736098242, 725644398, 641176081,
     1      631007004, 652095509, 651977348, 678403473, 699396944,
     2      709982538, 715309381, 699898672, 689787054, 668937254,
     3      631832743, 573957253, 501914406, 476767540, 457404757,
     4      412823153, 316436315, 275800443, 212452579, 186391926,
     5      143142223, 156373394, 173914695, 184134316, 175908267,
     6      170830834, 166099274, 158160555, 155366385, 152992165,
     7      155788350, 140265095, 195023239, 200193572, 210003328,
     8      204571271, 196889234, 191596723, 195948660, 192254853,
     9      217481208, 228527403, 236995435, 245395660, 246141505/
            data (iy(i),i= 251, 300)/
     +      235634065, 227706409, 222547889, 504862738, 517772865,
     1      533265162, 546199608, 548729992, 551267290, 556491137,
     2      572152567, 587786102, 582453489, 577180147, 561485386,
     3      550874090, 524700975, 524875212, 514524031, 504154873,
     4      498817253, 498604679, 491014862, 478181410, 480935097,
     5      489040852, 788645267, 799221802, 799227524, 791356945,
     6      794052219, 780864763, 772967958, 770248508, 775501299,
     7      772803307, 778066111, 772757530, 756878948, 746292973,
     8      746296310, 764827394, 445106268, 455691290, 466757298,
     9      472859049, 462697220, 463022566, 457950878, 452547693/
            data (iy(i),i= 301, 303)/
     +      446819878, 441208935, 440393257/
c
        sp(2)='lake superior'
        sp(1)='lake superior'
        sp(3)='lake superior'
c
        do i=1,303
            ibndry(1,i)=i-1
            ibndry(2,i)=i
            ibndry(3,i)=0
            ibndry(4,i)=2
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=0
            vx(i)=real(ix(i),rknd)*1.0e-7_rknd
            vy(i)=real(iy(i),rknd)*1.0e-7_rknd
        enddo
        ibndry(1,1)=225
        ibndry(1,226)=237
        ibndry(1,238)=245
        ibndry(1,246)=253
        ibndry(1,254)=276
        ibndry(1,277)=292
        ibndry(1,293)=303
        do i=304,313
            ibndry(3,i)=0
            ibndry(4,i)=0
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=0
        enddo
        ibndry(1,304)=153
        ibndry(2,304)=299
        ibndry(1,305)=154
        ibndry(2,305)=296
        ibndry(1,306)=189
        ibndry(2,306)=291
        ibndry(1,307)=188
        ibndry(2,307)=277
        ibndry(1,308)=254
        ibndry(2,308)=217
        ibndry(1,309)=216
        ibndry(2,309)=256
        ibndry(1,310)=221
        ibndry(2,310)=249
        ibndry(1,311)=246
        ibndry(2,311)=241
        ibndry(1,312)=229
        ibndry(2,312)=245
        ibndry(1,313)=237
        ibndry(2,313)=18
c
        do i=1,nbf
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
c
c       make itnode
c
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd5(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), dimension(200) :: ix,iy,ib1,ib2,ib3,ib4,
     +          ib6,isf1,isf2
            integer(kind=iknd), dimension(5) :: ixc,iyc
            integer(kind=iknd), save :: ntf,nvf,nbf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), save :: x,y,xc,yc,hmax,grade
            character(len=80), dimension(100) :: sp,su
cy
c       domain in the shape of the  dd20 logo
c       dirichlet boundary conditions on all boundary edges
c
            external sxy
            data (ix(i),i=1,100)/
     +    385, 3192, 3516, 5999,10000,10000, 5999, 3516, 3192,  385,
     1   1108, 2083, 4301, 5276, 5276, 4301, 2083, 1108,    0,    0,
     2   6384, 6384, 6598, 6736, 6736, 6874, 6912, 6930, 6948, 6967,
     3   6983, 6996, 7003, 7008, 7010, 7008, 7003, 6996, 6983, 6967,
     4   6948, 6930, 6912, 6874, 6754, 6754, 6885, 6939, 6985, 7021,
     5   7048, 7065, 7087, 7103, 7116, 7127, 7136, 7143, 7148, 7143,
     6   7136, 7127, 7116, 7103, 7087, 7065, 7048, 7021, 6985, 6939,
     7   6885, 6598, 7332, 7470, 7470, 7608, 7647, 7665, 7683, 7701,
     8   7718, 7730, 7738, 7743, 7745, 7743, 7738, 7730, 7718, 7701,
     9   7683, 7665, 7647, 7608, 7488, 7488, 7619, 7674, 7719, 7756/
            data (ix(i),i=101,146)/
     +   7783, 7799, 7821, 7838, 7850, 7861, 7870, 7878, 7883, 7878,
     1   7870, 7861, 7850, 7838, 7821, 7799, 7783, 7756, 7719, 7674,
     2   7619, 7332, 8354, 8476, 8641, 8430, 8630, 8470, 8354, 8238,
     3   8078, 8278, 8067, 8232, 9112, 9234, 9399, 9189, 9389, 9229,
     4   9112, 8996, 8836, 9036, 8825, 8990/
            data (iy(i),i=1,100)/
     +   3379, 1758, 3379, 3379, 3379, 6620, 6620, 6620, 8241, 6620,
     1   2516, 1953, 1953, 2516, 7483, 8046, 8046, 7483, 5562, 4437,
     2   4437, 5562, 4937, 4937, 5488, 5488, 5477, 5470, 5457, 5437,
     3   5413, 5386, 5359, 5322, 5268, 5213, 5177, 5150, 5122, 5099,
     4   5079, 5066, 5059, 5048, 5048, 4937, 4937, 4944, 4959, 4975,
     5   4995, 5013, 5041, 5068, 5095, 5122, 5159, 5195, 5268, 5341,
     6   5377, 5413, 5441, 5468, 5495, 5522, 5541, 5561, 5577, 5591,
     7   5599, 5599, 4937, 4937, 5488, 5488, 5477, 5470, 5457, 5437,
     8   5413, 5386, 5359, 5322, 5268, 5213, 5177, 5150, 5122, 5099,
     9   5079, 5066, 5059, 5048, 5048, 4937, 4937, 4944, 4959, 4975/
            data (iy(i),i=101,146)/
     +   4995, 5013, 5041, 5068, 5095, 5122, 5159, 5195, 5268, 5341,
     1   5377, 5413, 5441, 5468, 5495, 5522, 5541, 5561, 5577, 5591,
     2   5599, 5599, 5170, 4937, 4937, 5275, 5601, 5601, 5382, 5601,
     3   5601, 5275, 4937, 4937, 5170, 4937, 4937, 5275, 5601, 5601,
     4   5382, 5601, 5601, 5275, 4937, 4937/
            data (isf1(i),i=1,24)/
     +   3192, 3192,    0,    0,    0,    0,    0, 3192, 3192, 3192,
     1   3192,    0, 3192, 3192, 3192, 3192, 3192, 3192, 3192, 3192,
     2   3192, 3192, 3192, 3192/
            data (isf2(i),i=1,24)/
     +   5000, 5000,    0,    0,    0,    0,    0, 5000, 5000, 5000,
     1   5000,    0, 5000, 5000, 5000, 5000, 5000, 5000, 5000, 5000,
     2   5000, 5000, 5000, 5000/
            data (ib1(i),i=1,100)/
     +      1,    2,    3,    4,    5,    6,    7,    7,    9,   10,
     1      4,    8,   11,   12,   13,   14,   15,   16,   17,   18,
     2     19,   20,   21,   22,   23,   24,   25,   26,   27,   28,
     3     29,   30,   31,   32,   33,   34,   35,   36,   37,   38,
     4     39,   40,   41,   42,   43,   44,   45,   46,   47,   48,
     5     49,   50,   51,   52,   53,   54,   55,   56,   57,   58,
     6     59,   60,   61,   62,   63,   64,   65,   66,   67,   68,
     7     69,   70,   71,   72,   21,   73,   74,   75,   76,   77,
     8     78,   79,   80,   81,   82,   83,   84,   85,   86,   87,
     9     88,   89,   90,   91,   92,   93,   94,   95,   96,   97/
            data (ib1(i),i=101,153)/
     +     98,   99,  100,  101,  102,  103,  104,  105,  106,  107,
     1    108,  109,  110,  111,  112,  113,  114,  115,  116,  117,
     2    118,  119,  120,  121,  122,   54,  123,  124,  125,  126,
     3    127,  128,  129,  130,  131,  132,  133,  134,  104,  135,
     4    136,  137,  138,  139,  140,  141,  142,  143,  144,  145,
     5    146,  125,    5/
            data (ib2(i),i=1,100)/
     +     11,   13,    4,    5,    6,    7,    8,   15,   17,   19,
     1     21,    3,   12,    2,   14,    4,   16,    9,   18,   10,
     2     20,    1,   22,    7,   24,   25,   26,   27,   28,   29,
     3     30,   31,   32,   33,   34,   35,   36,   37,   38,   39,
     4     40,   41,   42,   43,   44,   45,   46,   47,   48,   49,
     5     50,   51,   52,   53,   54,   55,   56,   57,   58,   59,
     6     60,   61,   62,   63,   64,   65,   66,   67,   68,   69,
     7     70,   71,   72,   23,   23,   74,   75,   76,   77,   78,
     8     79,   80,   81,   82,   83,   84,   85,   86,   87,   88,
     9     89,   90,   91,   92,   93,   94,   95,   96,   97,   98/
            data (ib2(i),i=101,153)/
     +     99,  100,  101,  102,  103,  104,  105,  106,  107,  108,
     1    109,  110,  111,  112,  113,  114,  115,  116,  117,  118,
     2    119,  120,  121,  122,   73,   73,  124,  125,  126,  127,
     3    128,  129,  130,  131,  132,  133,  134,  123,  133,  136,
     4    137,  138,  139,  140,  141,  142,  143,  144,  145,  146,
     5    135,  145,  137/
            data (ib3(i),i=1,24)/
     +  1,1,0,0,0,0,0,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1/
            data (ib4(i),i=1,153)/
     +  2,2,0,2,2,2,0,2,2,2,0,0,2,2,2,2,2,2,2,2,2,2,0,0,2,2,2,2,2,2,
     1  2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
     2  2,2,2,2,2,2,2,2,2,2,2,2,2,2,0,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
     3  2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
     4  2,2,2,2,2,0,2,2,2,2,2,2,2,2,2,2,2,2,0,2,2,2,2,2,2,2,2,2,2,2,
     5  2,0,0/
            data (ib6(i),i=1,153)/
     +   1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12, 1, 1, 2, 2, 8, 8, 9, 9,
     1  10,10,11,11, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     2   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     3   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,
     4   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     5   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     6   1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1,
     7   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0/
            data ntf,nvf,nbf/3,146,153/
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
c
        sp(2)='dd20'
        sp(1)='dd20'
        sp(3)='dd20'
c
        do i=1,nvf
            vx(i)=real(ix(i),rknd)
            vy(i)=real(iy(i),rknd)
        enddo
        do i=1,nbf
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
            ibndry(1,i)=ib1(i)
            ibndry(2,i)=ib2(i)
            ibndry(3,i)=0
            ibndry(4,i)=ib4(i)
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=ib6(i)
        enddo
        do i=1,24
            sf(1,i)=real(isf1(i),rknd)
            sf(2,i)=real(isf2(i),rknd)
            ibndry(3,i)=ib3(i)
        enddo
c
c       uncomment to delete the ddxx symbols
c
cc      nbf=24
cc      nvf=22
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
c
c       refine long curved edges
c
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        itnode(5,1)=3
        itnode(5,2)=2
        itnode(5,3)=1
        itnode(5,4)=1
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd6(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), save, dimension(159) :: b1,b2,b3,b4,b5
            integer(kind=iknd), save, dimension(25) :: c1,c2,c4
            integer(kind=iknd), save, dimension(135) :: ix,iy
            integer(kind=iknd), dimension(93) :: ixp,iyp
            integer(kind=iknd), save :: ntf,ncf,nbf,nvf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), save, dimension(93) :: xp,yp
            real(kind=rknd), save :: hmax,grade,x,y
            character(len=80), dimension(100) :: sp,su
cy
c       an irregular region composed of 25 regions, in the shape of,
c       the at&t logo.
c       the data was provided by r. kent smith
c
c
            external sxy
            data ntf,nvf,nbf,ncf/25,135,159,93/
c
            data (ix(i),i=   1,  50)/
     +      -29088581,  -2115240,  30000000,         0, -40919900,
     1      -31212461, -24240000, -11470000,  -5369100,  20187540,
     2       27160001,  43160391,  33551569,  21120000,  10151761,
     3       -5109150, -13640000, -21440001, -29284019, -36108899,
     4      -49742250, -41352282, -32990000,  -4100000,   1129340,
     5       24291310,  34749999,  50882888,  45282569,  28659999,
     6       22388289,   1792850,  -8660000, -31520000, -38492460,
     7      -46628761, -55230289, -44090900, -37990000,   -510000,
     8        5761710,  24727149,  35179999,  56014829,  52477179,
     9       34890001,  24437151,   5216990,  -2100000, -36550000/
            data (ix(i),i=  51, 100)/
     +      -45265570, -53460388, -58688860, -48302279, -39940000,
     1       -1400000,   6444020,  27805979,  35650001,  58897629,
     2       56731110,  36570001,  27162440,   5342280,  -3020000,
     3      -39630001, -46602459, -57378292, -59963450, -45582852,
     4      -35130000,  -4560000,   3284020,  59990859,  59261298,
     5        5425570,  -3290000, -37790000, -46505570, -59552770,
     6      -59416080, -41115570, -32400000,  -7110000,   -137540,
     7       21172869,  30760000,  59416080,  59990859,  35500000,
     8       26784430,  -1595980,  -9440000, -32820001, -39795649,
     9      -59963450, -57063389, -30675910, -11144090,  13666790/
            data (ix(i),i= 101, 135)/
     +       21340001,  57063389,  58897629,  28210001,  22629480,
     1      -10966530, -34503360, -58897629, -52477179, -26954920,
     2      -16485080,   7010050,  10500000,  52477179,  56014829,
     3       19870000,  14635080, -15375080, -25844920, -56014829,
     4      -45282569,  45282569,  50882888, -50882888, -33551569,
     5       33551569,  43160391, -43160391, -30000000,  30000000,
     6              0,         0, 400000000,         0,-400000000/
            data (iy(i),i=   1,  50)/
     +       52477179,  52805071,  51961522,  60000000,  43881221,
     1       44595580,  44899998,  44899998,  44633632,  42354422,
     2       42049999,  41679502,  49742250,  49850001,  49179149,
     3       47371769,  46849999,  46849999,  47192478,  47918129,
     4       33551569,  34211750,  34649999,  34649999,  34421680,
     5       32256639,  31800001,  31795161,  39363539,  39600000,
     6       39271309,  37097809,  36550000,  36550000,  36854420,
     7       37759221,  23443871,  24183631,  24449999,  24449999,
     8       24121311,  22247810,  21700001,  21502080,  29088581,
     9       29500000,  28952191,  26783471,  26400001,  26400001/
            data (iy(i),i=  51, 100)/
     +       26780529,  27239430,  12474700,  13811750,  14250000,
     1       14250000,  13907520,  11892480,  11550000,  11448539,
     2       19534090,  19349999,  18856970,  16688250,  16250000,
     3       16250000,  16554420,  17542300,   2093970,   3702190,
     4        4250000,   4250000,   3907520,   1047140,   9386070,
     5        6630530,   6250000,   6250000,   6630530,   7312160,
     6       -8350390,  -6480530,  -6100000,  -6100000,  -6404420,
     7       -8331420,  -8750000,  -8350390,  -1047140,   -950000,
     8       -1330530,  -3657520,  -4000000,  -4000000,  -3756410,
     9       -2093970, -18541020, -16861030, -16861030, -18732049/
            data (iy(i),i= 101, 135)/
     +      -19000000, -18541020, -11448539, -11200000, -11394880,
     1      -13924609, -13719200, -11448539, -29088581, -28201380,
     2      -28201380, -29089079, -29150000, -29088581, -21502080,
     3      -21350000, -21441381, -22548621, -22548621, -21502080,
     4      -39363539, -39363539, -31795161, -31795161, -49742250,
     5      -49742250, -41679502, -41679502, -51961522, -51961522,
     6      -60000000, 400000000,         0,-400000000,         0/
c
            data (ixp(i),i=   1,  50)/
     +      -15328120,  15529140, -15022800, -27729549,  -8416640,
     1       23670449,  38567259,  15625629,  -9366600, -25365739,
     2      -38567259, -37176881,  -1482840,  29515669,  48231411,
     3       25519841,  -3426400, -35009551, -48231411, -41043358,
     4        2630160,  29946401,  54378471,  29656401,   1563520,
     5      -40911942, -54378471, -44126878,   2525740,  31724260,
     6       57955551,  31859760,   1166880, -43119550, -58088861,
     7      -40363598,   -634260,  59771681,   1071940, -42151942,
     8      -59815040, -36761940,  -3620450,  25961871,  59815040,
     9       31138060,  -5514260, -36309950, -59771681, -20910001/
            data (ixp(i),i=  51,  93)/
     +       17501060,  58088861,  25418041, -22738979, -58088861,
     1      -21719999,   8754760,  54378471,  17252140, -20610001,
     2      -54378471,  48231411, -48231411,  38567259, -38567259,
     3       15529140, -15529140,  31795161, -32678339,  44236641,
     4      -43881221,  51697750, -51697750,  56381559, -56381559,
     5       59088469, -59177141,  60000000, -60000000,  59177141,
     6      -59177141,  56558490, -56558490,  51697750, -51697750,
     7       44236641, -44236641,  31795161, -31795161, 282842712,
     8      282842712,-282842712,-282842712/
            data (iyp(i),i=   1,  50)/
     +       53786950,  57955551,  58088861,  44823861,  44833379,
     1       42126141,  45962672,  49682131,  46980562,  46935658,
     2       45962672,  34540360,  34592891,  31914210,  35689371,
     3       39517770,  36687050,  36626141,  35689371,  24383380,
     4       24367771,  21837051,  25357101,  29362950,  26495931,
     5       26495180,  25357101,  14140360,  14164340,  11635660,
     6       15529140,  19226660,  16359640,  16326140,  15022800,
     7        4112950,   4164340,   5229340,   6345180,   6345180,
     8        4707550,  -6195180,  -6176140,  -8645300,  -4707550,
     9       -1045180,  -3914340,  -3939080,  -5229340, -16520000/
            data (iyp(i),i=  51,  93)/
     +      -18932990, -15022800, -11248730, -14284290, -15022800,
     1      -28110001, -29134769, -25357101, -21372850, -22639999,
     2      -25357101, -35689371, -35689371, -45962672, -45962672,
     3      -57955551, -57955551,  50882888,  50320230,  40535412,
     4       40919900,  30452299,  30452299,  20521209,  20521209,
     5       10418890,   9902860,         0,         0,  -9902860,
     6       -9902860, -20028410, -20028410, -30452299, -30452299,
     7      -40535412, -40535412, -50882888, -50882888, 282842712,
     8     -282842712,-282842712, 282842712/
c
            data (b1(i),i=1,100)/
     +        1,  2,  3,  4,  5,  6,  7,  8,  9, 10,
     1       11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
     2       21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
     3       31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
     4       41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
     5       51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
     6       61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
     7       71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
     8       81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
     9       91, 92, 93, 94, 95, 96, 97, 98, 99,100/
            data (b1(i),i=101,159)/
     +      101,102,103,104,105,106,107,108,109,110,
     1      111,112,113,114,115,116,117,118,119,120,
     2      121,122,123,124,125,126,127,128,129,130,
     3      131,  3,  1, 12,  5, 28, 21, 44, 37, 60,
     4       53, 74, 69, 88, 81,102, 97,114,109,122,
     5      121,126,125,132,133,134,135,132,134/
            data (b2(i),i=1,100)/
     +        2,  3,  4,  1,  6,  7,  8,  9, 10, 11,
     1       12, 13, 14, 15, 16, 17, 18, 19, 20,  5,
     2       22, 23, 24, 25, 26, 27, 28, 29, 30, 31,
     3       32, 33, 34, 35, 36, 21, 38, 39, 40, 41,
     4       42, 43, 44, 45, 46, 47, 48, 49, 50, 51,
     5       52, 37, 54, 55, 56, 57, 58, 59, 60, 61,
     6       62, 63, 64, 65, 66, 67, 68, 53, 70, 71,
     7       72, 73, 74, 75, 76, 77, 78, 79, 80, 69,
     8       82, 83, 84, 85, 86, 87, 88, 89, 90, 91,
     9       92, 93, 94, 95, 96, 81, 98, 99,100,101/
            data (b2(i),i=101,159)/
     +      102,103,104,105,106,107,108, 97,110,111,
     1      112,113,114,115,116,117,118,119,120,109,
     2      122,123,124,121,126,127,128,125,130,131,
     3      129, 13, 20, 29, 36, 45, 52, 61, 68, 75,
     4       80, 89, 96,103,108,115,120,123,124,127,
     5      128,130,129,133,134,135,132,  4,131/
            data (b3(i),i=1,100)/
     +        1,  0,  2,  3,  0,  4,  0,  5,  0,  6,
     1        0,  7,  0,  8,  0,  9,  0, 10,  0, 11,
     2        0, 12,  0, 13,  0, 14,  0, 15,  0, 16,
     3        0, 17,  0, 18,  0, 19,  0, 20,  0, 21,
     4        0, 22,  0, 23,  0, 24,  0, 25,  0, 26,
     5        0, 27,  0, 28,  0, 29,  0, 30,  0, 31,
     6        0, 32,  0, 33,  0, 34,  0, 35,  0, 36,
     7        0, 37,  0, 38,  0, 39,  0, 40,  0, 41,
     8        0, 42,  0, 43,  0, 44,  0, 45,  0, 46,
     9        0, 47,  0, 48,  0, 49,  0, 50,  0, 51/
            data (b3(i),i=101,159)/
     +        0, 52,  0, 53,  0, 54,  0, 55,  0, 56,
     1        0, 57,  0, 58,  0, 59,  0, 60,  0, 61,
     2        0, 62,  0, 63,  0, 64,  0, 65,  0, 66,
     3       67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
     4       77, 78, 79, 80, 81, 82, 83, 84, 85, 86,
     5       87, 88, 89, 90, 91, 92, 93,  0,  0/
            data (b4(i),i=1,100)/
     +        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     1        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     2        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     3        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     4        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     5        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     6        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     7        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     8        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     9        0,  0,  0,  0,  0,  0,  0,  0,  0,  0/
            data (b4(i),i=101,159)/
     +        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     1        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     2        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     3        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     4        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     5        0,  0,  0,  2,  2,  2,  2,  0,  0/
            data (b5(i),i=1,100)/
     +        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     1        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     2        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     3        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     4        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     5        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     6        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     7        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     8        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     9        0,  0,  0,  0,  0,  0,  0,  0,  0,  0/
            data (b5(i),i=101,159)/
     +        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     1        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     2        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     3        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     4        0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     5        0,  0,  0,  1,  1,  1,  1,  0,  0/
c
            data (c1(i),i=1,25)/
     +        4,  1, 20,  5, 36, 21, 52, 37, 68, 53,
     1       80, 69, 96, 81,108, 97,120,109,124,121,
     2      128,125,131,134,134/
            data (c2(i),i=1,25)/
     +        4,133, 20,135, 36,137, 52,139, 68,141,
     1       80,143, 96,145,108,147,120,149,124,151,
     2      128,153,130,155,159/
            data (c4(i),i=1,25)/
     +        1,  2,  1,  2,  1,  2,  1,  2,  1,  2,
     1        1,  2,  1,  2,  1,  2,  1,  2,  1,  2,
     2        1,  2,  1,  3,  3/
c
            data hmax,grade/0.05e0_rknd,1.5e0_rknd/
c
        sp(2)='at&t logo'
        sp(1)='at&t logo'
        sp(3)='at&t logo'
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
c
        do i=1,nvf
            vx(i)=real(ix(i),rknd)*1.0e-7_rknd
            vy(i)=real(iy(i),rknd)*1.0e-7_rknd
        enddo
        do i=1,ncf
            xp(i)=real(ixp(i),rknd)*1.0e-7_rknd
            yp(i)=real(iyp(i),rknd)*1.0e-7_rknd
        enddo
c
        do i=1,nbf
            ibndry(1,i)=b1(i)
            ibndry(2,i)=b2(i)
            ibndry(3,i)=b3(i)
            ibndry(4,i)=b4(i)
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=b5(i)
c
            if(ibndry(3,i)>0) then
                i1=ibndry(1,i)
                i2=ibndry(2,i)
                i3=ibndry(3,i)
                call centre(vx(i1),vy(i1),vx(i2),vy(i2),
     +              xp(i3),yp(i3),sf(1,i),sf(2,i))
            else
                sf(1,i)=0.0e0_rknd
                sf(2,i)=0.0e0_rknd
           endif
        enddo
c
        do i=1,ntf
            itnode(1,i)=c1(i)
            itnode(2,i)=c2(i)
            itnode(3,i)=0
            itnode(5,i)=c4(i)
        enddo
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd7(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), save, dimension(499) :: ix,iy
            integer(kind=iknd), save :: ntf,nvf,nbf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), save :: hmax,grade
            character(len=80), dimension(100) :: sp,su
cy
c       domain in the shape of the north sea.
c       the domain contains 29 islands.
c       the skeleton has four regions.
c
c       data for this domain was provided by gabriel wittum and
c       wolfgang hoffman from the university of stuttgart.
c
            external sxy
            data ntf,nvf,nbf/4,499,531/
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
c
            data (ix(i),i=   1,  80)/
     +       8820,  8030,     0,     0,  1260,  3510,  3530,  3810,
     1       3840,  3780,  3640,  3700,  3860,  3950,  4140,  4270,
     2       4440,  4600,  4680,  4710,  4960,  5040,  5210,  5320,
     3       5420,  5530,  5890,  6000,  6030,  6190,  6250,  6520,
     4       6550,  6680,  6820,  6770,  6770,  6820,  6880,  6990,
     5       7040,  7120,  7180,  7230,  7260,  7210,  7120,  7040,
     6       7810,  7810,  7700,  7560,  7480,  7480,  7420,  7530,
     7       7560,  7640,  7750,  7840,  7890,  7970,  8050,  8190,
     8       8300,  8440,  8520,  8490,  8410,  8600,  8630,  8630,
     9       8580,  8520,  8410,  8330,  8300,  8300,  8300,  8360/
            data (ix(i),i=  81, 160)/
     +       8410,  8470,  8520,  8520,  8520,  8630,  8710,  8770,
     1       8880,  8960,  9070,  9150,  9260,  9320,  9480,  9590,
     2       9730,  9840,  9890, 10000, 10000,  9890,  9840,  9780,
     3       9730,  9700,  9700,  9700,  9640,  9640,  9590,  9530,
     4       9510,  9480,  9530,  9670,  9730,  9780,  9840,  9860,
     5       9920,  9950,  9920,  9860,  9860,  9840,  9810,  9780,
     6       9780,  9700,  9640,  9560,  9510,  9480,  9450,  9420,
     7       9420,  9400,  9400,  9420,  9480,  9480,  9480,  9480,
     8       9420,  9370,  9320,  9260,  9230,  9150,  9100,  9010,
     9       8990,  8930,  8880,  8790,  8710,  8680,  8660,  8660/
            data (ix(i),i= 161, 240)/
     +       8770,  8820,  8930,  8960,  8930,  8880,  8820,  8790,
     1       8880,  8960,  9070,  9210,  9290,  9370,  9450,  9530,
     2       9620,  9670,  9780,  9840,  9890,  9920, 10000, 10000,
     3       9860,  9810,  9730,  9700,  9640,  9640,  9560,  9450,
     4       9510,  9480,  9420,  9400,  9320,  9210,  9180,  8990,
     5       8990,  9040,  9010,  9010,  8960,  8880,  8880,  8740,
     6       8710,  8740,  8740,  8820,  8850,  8850,  8900,  8880,
     7       8930,  8930,  8770,  8790,  8820,  9010,  9010,  8930,
     8       8850,  8740,  8660,  8660,  8770,  8820,  9010,  8960,
     9       8960,  9010,  8900,  8190,  8080,  8050,  8140,  8410/
            data (ix(i),i= 241, 320)/
     +       8410,  8300,  8330,  7890,  7970,  7950,  7950,  7840,
     1       7890,  7260,  7290,  7210,  7180,  7100,  7100,  7210,
     2       6360,  6550,  6330,  6220,  6160,  6250,  5920,  6050,
     3       5840,  5750,  5730,  5640,  5620,  5590,  5590,  5670,
     4       4900,  4960,  5370,  5420,  5370,  5290,  5180,  5180,
     5       5040,  5040,  4900,  4770,  4820,  4680,  4600,  4600,
     6       4470,  4440,  4300,  4030,  3970,  3890,  3810,  3950,
     7       3150,  3370,  3640,  3510,  3320,  3150,  2960,  2930,
     8       3040,  3100,  3040,  2960,  2930,  2900,  2300,  2490,
     9       2600,  2660,  2710,  2710,  2490,  2470,  2410,  2490/
            data (ix(i),i= 321, 400)/
     +       2330,  2220,  2220,  2030,  2000,  1950,  1920,  1920,
     1       1970,  2030,  1530,  1560,  1700,  1750,  1560,  1530,
     2       1450,  1450,  1480,  1370,  1370,  1340,  1340,  1260,
     3       1230,  6380,  6280,  6270,  6270,  6250,  6330,  6330,
     4       9700,  9750,  9860,  9840,  9810,  9750,  9670,  9480,
     5       9420,  9400,  9420,  9400,  9420,  9510,  9560,  9640,
     6       8580,  8520,  8520,  8580,  8250,  8250,  8300,  8330,
     7       8330,  8300,  8250,  8140,  8190,  8270,  8300,  8300,
     8       8250,  8190,  8220,  8550,  8580,  8580,  8490,  8470,
     9       8380,  8220,  8140,  8190,  8880,  8880,  8820,  8770/
            data (ix(i),i= 401, 480)/
     +       8660,  8580,  8520,  8470,  8440,  8490,  8580,  8660,
     1       8710,  8790,  8900,  9070,  8930,  9070,  7840,  7810,
     2       7860,  7950,  7950,  8000,  8000,  7920,  7890,  7810,
     3       7750,  7700,  7670,  7700,  7670,  7890,  8360,  8410,
     4       8440,  8520,  8550,  8600,  8630,  8600,  8630,  8660,
     5       8600,  8580,  8580,  8220,  8160,  8110,  8030,  8030,
     6       8030,  8080,  8160,  8360,  8220,  8080,  8050,  8140,
     7       8140,  8080,  8000,  7920,  7890,  7860,  7860,  7890,
     8       7950,  8030,  8140,  8220,  8300,  8190,  8110,  8030,
     9       7950,  7920,  7840,  7810,  7810,  7730,  7700,  7700/
            data (ix(i),i= 481, 499)/
     +       7700,  7670,  7670,  7700,  7670,  7640,  7590,  7620,
     1       7620,  7640,  7670,  7730,  7810,  7950,  8030,  8140,
     2       8140,  8220,  8250/
            data (iy(i),i=   1,  80)/
     +      10000, 10000, 10000,     0,     0,     0,    90,   120,
     1        240,   300,   420,   520,   700,   820,   970,  1060,
     2       1090,  1120,  1090,  1120,  1120,  1060,  1060,  1090,
     3       1210,  1240,  1240,  1270,  1300,  1270,  1360,  1360,
     4       1300,  1270,  1240,  1150,  1000,   880,   850,   820,
     5        700,   550,   480,   330,   210,   150,   120,     0,
     6          0,    90,    90,   120,   120,   240,   360,   450,
     7        580,   670,   700,   640,   580,   480,   390,   330,
     8        300,   240,   210,   150,     0,     0,   120,   180,
     9        270,   450,   670,   790,   910,  1030,  1180,  1330/
            data (iy(i),i=  81, 160)/
     +       1450,  1580,  1760,  1940,  2060,  2150,  2270,  2330,
     1       2390,  2390,  2300,  2210,  2120,  1970,  1970,  2000,
     2       2030,  2060,  2060,  2060,  2420,  2480,  2580,  2580,
     3       2670,  2700,  2760,  2820,  2880,  2970,  3030,  3120,
     4       3210,  3270,  3330,  3360,  3330,  3300,  3270,  3300,
     5       3330,  3360,  3580,  3640,  3700,  3760,  3790,  3880,
     6       3940,  3940,  3910,  3910,  3940,  3970,  4030,  4090,
     7       4150,  4180,  4300,  4390,  4450,  4610,  4700,  4880,
     8       4910,  4910,  4880,  4880,  4910,  4940,  4940,  4910,
     9       4880,  4790,  4790,  4850,  4940,  5060,  5180,  5240/
            data (iy(i),i= 161, 240)/
     +       5240,  5270,  5300,  5330,  5390,  5420,  5420,  5520,
     1       5580,  5640,  5640,  5670,  5700,  5730,  5670,  5670,
     2       5670,  5730,  5880,  5910,  5970,  5970,  5940,  6420,
     3       6520,  6610,  6640,  6700,  6850,  6910,  6910,  6910,
     4       6970,  7090,  7210,  7240,  7300,  7390,  7580,  7670,
     5       7730,  7790,  7880,  7970,  8120,  8150,  8240,  8420,
     6       8480,  8520,  8700,  8730,  8790,  8850,  8850,  9150,
     7       9150,  9300,  9610,  9700,  9700,  6550,  6420,  6210,
     8       6180,  6180,  6300,  6480,  6480,  6580,  6610,  3580,
     9       3450,  3330,  3520,  2850,  2760,  2820,  2910,  2640/
            data (iy(i),i= 241, 320)/
     +       2580,  2550,  2670,  2090,  2090,  2060,  1970,  2000,
     1       2120,  1420,  1360,  1360,  1270,  1270,  1360,  1450,
     2       1790,  1730,  1760,  1700,  1790,  1820,  1730,  1700,
     3       1580,  1580,  1610,  1580,  1550,  1550,  1610,  1730,
     4       1450,  1550,  1580,  1480,  1480,  1520,  1520,  1480,
     5       1520,  1390,  1360,  1420,  1330,  1330,  1360,  1450,
     6       1390,  1330,  1300,  1270,  1210,  1240,  1240,  1390,
     7       1120,  1150,  1150,  1090,  1030,  1000,  1000,  1090,
     8        940,   850,   760,   760,   820,   910,   610,   730,
     9        730,   640,   610,   580,   550,   480,   480,   390/
            data (iy(i),i= 321, 400)/
     +        390,   520,   610,   240,   210,   120,   120,   300,
     1        300,   300,   300,   270,   270,   180,   180,   210,
     2        180,   240,   300,   150,    90,    90,    60,    60,
     3        150,  4270,  4220,  4210,  4180,  4180,  4330,  4270,
     4       6390,  6390,  6360,  6210,  6180,  6060,  6000,  6000,
     5       6030,  6060,  6150,  6240,  6270,  6300,  6390,  6390,
     6       6000,  5970,  6060,  6060,  6000,  5940,  5940,  5880,
     7       5820,  5730,  5760,  5880,  6000,  6030,  6450,  6300,
     8       6240,  6330,  6480,  6700,  6670,  6640,  6610,  6670,
     9       6730,  6640,  6640,  6790,  7210,  7150,  7120,  7090/
            data (iy(i),i= 401, 480)/
     +       7090,  7060,  7000,  7030,  7090,  7150,  7150,  7180,
     1       7180,  7210,  7240,  7330,  7360,  7420,  7480,  7390,
     2       7300,  7150,  7060,  7060,  7000,  7000,  6970,  7030,
     3       7090,  7180,  7270,  7300,  7390,  7610,  7820,  7820,
     4       7880,  7880,  7850,  7850,  7790,  7730,  7670,  7640,
     5       7610,  7520,  7390,  7390,  7480,  7480,  7520,  7640,
     6       7730,  7820,  7850,  7850,  9730,  9700,  9640,  9610,
     7       9550,  9520,  9450,  9360,  9270,  9210,  9030,  8850,
     8       8760,  8700,  8640,  8640,  8640,  8550,  8520,  8450,
     9       8520,  8580,  8670,  8640,  8550,  8520,  8480,  8330/
            data (iy(i),i= 481, 499)/
     +       8180,  8090,  8030,  8030,  7880,  7760,  7850,  8060,
     1       8120,  8700,  8850,  9060,  9330,  9640,  9790,  9790,
     2       9760,  9760,  9730/
c
        sp(2)='north sea'
        sp(1)='north sea'
        sp(3)='north sea'
c
        do i=1,nvf
            vx(i)=real(ix(i),rknd)*1.0e-4_rknd
            vy(i)=real(iy(i),rknd)*1.0e-4_rknd
        enddo
        do i=1,499
            ibndry(1,i)=i
            ibndry(2,i)=i+1
            ibndry(3,i)=0
            ibndry(4,i)=2
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=0
        enddo
        ibndry(2,221)=1
        ibndry(2,231)=222
        ibndry(2,235)=232
        ibndry(2,239)=236
        ibndry(2,243)=240
        ibndry(2,249)=244
        ibndry(2,256)=250
        ibndry(2,262)=257
        ibndry(2,272)=263
        ibndry(2,283)=273
        ibndry(2,288)=284
        ibndry(2,296)=289
        ibndry(2,304)=297
        ibndry(2,310)=305
        ibndry(2,323)=311
        ibndry(2,330)=324
        ibndry(2,339)=331
        ibndry(2,345)=340
        ibndry(2,352)=346
        ibndry(2,368)=353
        ibndry(2,372)=369
        ibndry(2,382)=373
        ibndry(2,387)=383
        ibndry(2,393)=388
        ibndry(2,396)=394
        ibndry(2,411)=397
        ibndry(2,414)=412
        ibndry(2,430)=415
        ibndry(2,452)=431
        ibndry(2,499)=453
        do i=500,531
            ibndry(3,i)=0
            ibndry(4,i)=0
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=0
        enddo
        ibndry(1,500)=5
        ibndry(2,500)=344
        ibndry(1,501)=340
        ibndry(2,501)=337
        ibndry(1,502)=334
        ibndry(2,502)=327
        ibndry(1,503)=330
        ibndry(2,503)=321
        ibndry(1,504)=313
        ibndry(2,504)=308
        ibndry(1,505)=305
        ibndry(2,505)=303
        ibndry(1,506)=299
        ibndry(2,506)=295
        ibndry(1,507)=289
        ibndry(2,507)=287
        ibndry(1,508)=285
        ibndry(2,508)=283
        ibndry(1,509)=276
        ibndry(2,509)=270
        ibndry(1,510)=264
        ibndry(2,510)=261
        ibndry(1,511)=258
        ibndry(2,511)=255
        ibndry(1,512)=250
        ibndry(2,512)=248
        ibndry(1,513)=245
        ibndry(2,513)=242
        ibndry(1,514)=243
        ibndry(2,514)=236
        ibndry(1,515)=238
        ibndry(2,515)=349
        ibndry(1,516)=239
        ibndry(2,516)=235
        ibndry(1,517)=234
        ibndry(2,517)=114
        ibndry(1,518)=355
        ibndry(2,518)=184
        ibndry(1,519)=362
        ibndry(2,519)=224
        ibndry(1,520)=226
        ibndry(2,520)=372
        ibndry(1,521)=370
        ibndry(2,521)=376
        ibndry(1,522)=382
        ibndry(2,522)=385
        ibndry(1,523)=387
        ibndry(2,523)=395
        ibndry(1,524)=396
        ibndry(2,524)=393
        ibndry(1,525)=388
        ibndry(2,525)=403
        ibndry(1,526)=411
        ibndry(2,526)=412
        ibndry(1,527)=413
        ibndry(2,527)=443
        ibndry(1,528)=448
        ibndry(2,528)=430
        ibndry(1,529)=429
        ibndry(2,529)=486
        ibndry(1,530)=2
        ibndry(2,530)=495
        ibndry(1,531)=351
        ibndry(2,531)=379
c
        do i=1,nbf
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
c
c       make itnode
c
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd8(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), save, dimension(242) :: ix,iy
            integer(kind=iknd), save :: nvf,ntf,nbf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            character(len=80), dimension(100) :: sp,su
cy
c       domain in the shape of a three piece airfoil
c       the outer boundary is an irregular polygon
c       data for this domain was provided by tony chan from ucla
c
            external sxy
            data nvf,ntf,nbf/242,2,246/
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
c
            data (ix(i),i=   1,  50)/
     +       1423869,  1423991,  1424665,  1425178,  1425859,
     1       1426755,  1427902,  1429365,  1431152,  1433195,
     2       1435414,  1437821,  1440377,  1443041,  1445778,
     3       1448556,  1451382,  1459973,  1837655,  1898567,
     4       1967503,  2440690,  3106079,  3833284,  4532788,
     5       5195898,  5774643,  5842507,  5900000,  5951276,
     6       5988782,  6013725,  6032087,  6035385,  6029606,
     7       6030057,  6035735,  6040665,  6046612,  6053578,
     8       6061420,  6070220,  6080025,  6090760,  6102237,
     9       6114294,  6126850,  6139864,  6153344,  6167255/
            data (ix(i),i=  51, 100)/
     +       6181578,  6196272,  6211289,  6242031,  6257581,
     1       6273259,  6289045,  6304912,  6321770,  6392878,
     2       6423755,  6417585,  6407490,  6377786,  6144193,
     3       5309548,  4401639,  3504295,  2471459,  2220550,
     4       2145875,  2075516,  2009864,  1947819,  1889747,
     5       1835234,  1784572,  1692365,  1612390,  1577555,
     6       1547325,  1521174,  1497263,  1475941,  1457835,
     7       1443479,  1433836,  1428603,  1425751,  1424517,
     8       1423998,  6435090,  6435684,  6436864,  6438792,
     9       6441638,  6445023,  6448381,  6451503,  6453998/
            data (ix(i),i= 101, 150)/
     +       6456095,  6500034,  6519509,  6535504,  6544567,
     1       6554438,  6565221,  6589578,  6651070,  6689324,
     2       6709135,  6727578,  6800035,  6810905,  6844863,
     3       6851923,  6858552,  6860805,  6860022,  6858930,
     4       6857467,  6855525,  6842807,  6832409,  6818113,
     5       6779353,  6725953,  6697470,  6668447,  6636241,
     6       6599383,  6562072,  6528875,  6500614,  6477822,
     7       6461686,  6452649,  6448563,  6446756,  6445690,
     8       6444807,  6438994,  6436979,  6435461,  6613610,
     9       6614642,  6619005,  6625937,  6635061,  6690099/
            data (ix(i),i= 151, 200)/
     +       6901211,  7152194,  7376001,  7768795,  7782484,
     1       7790374,  7788553,  7783294,  7771088,  7137460,
     2       7087985,  7039049,  7014894,  6990902,  6967044,
     3       6943551,  6920739,  6899039,  6878715,  6859989,
     4       6845807,  6832375,  6818686,  6804809,  6790810,
     5       6776770,  6762792,  6748974,  6735394,  6722137,
     6       6709277,  6696883,  6685029,  6673808,  6663302,
     7       6653659,  6644960,  6637259,  6630630,  6625094,
     8       6620691,  6617758,  6615033,        0,   439662,
     9        679643,   137097,   451695,  1296440,  2369967/
            data (ix(i),i= 201, 242)/
     +       3584674,  4993196,  5491886,  6259193,  7002156,
     1       8059564,  9048663,  9714154,  9858797,  9717189,
     2       9646991,  9421564,  9776057,  9837757,  9778552,
     3       9474904,  9602283, 10000000,  9838499,  9353390,
     4       8462743,  7620040,  7015043,  6060112,  5380730,
     5       4499139,  3494099,  2738492,  1768388,   859205,
     6        233243,   251335,   160179,   272935,    94844,
     7        248004,   178170,   291858,   115478,   247026,
     8         53696,     1598/
            data (iy(i),i=   1,  50)/
     +       5254964,  5248606,  5238330,  5233882,  5229771,
     1       5225935,  5222397,  5219099,  5216099,  5213480,
     2       5211248,  5209293,  5207607,  5206144,  5204861,
     3       5203726,  5202719,  5200189,  5137385,  5125430,
     4       5115393,  5066633,  5035633,  5050275,  5103050,
     5       5192077,  5298739,  5308269,  5311492,  5316925,
     6       5327238,  5338626,  5351327,  5377162,  5413147,
     7       5443720,  5469391,  5485436,  5501084,  5516263,
     8       5530949,  5545001,  5558247,  5570642,  5582239,
     9       5593142,  5603388,  5612962,  5621775,  5629776/
            data (iy(i),i=  51, 100)/
     +       5636910,  5643118,  5648364,  5656289,  5659461,
     1       5661924,  5663409,  5663775,  5663349,  5659469,
     2       5659043,  5661356,  5663991,  5670498,  5708333,
     3       5810024,  5862396,  5861083,  5800256,  5769966,
     4       5754943,  5736496,  5716340,  5694184,  5670140,
     5       5644723,  5618731,  5565919,  5513785,  5488443,
     6       5464458,  5441265,  5417161,  5391871,  5365872,
     7       5339320,  5315813,  5297097,  5282008,  5271090,
     8       5262330,  5560329,  5553057,  5545997,  5539125,
     9       5533072,  5528495,  5525126,  5522541,  5520937/
            data (iy(i),i= 101, 150)/
     +       5519998,  5501572,  5494669,  5489811,  5487370,
     1       5485138,  5483336,  5480422,  5469531,  5460928,
     2       5455621,  5449801,  5423853,  5418400,  5395982,
     3       5391642,  5387869,  5389659,  5392076,  5394584,
     4       5397235,  5400257,  5417443,  5430332,  5446495,
     5       5486280,  5535828,  5559779,  5581622,  5601534,
     6       5618159,  5628927,  5632761,  5630887,  5625092,
     7       5617436,  5610627,  5606162,  5603417,  5601282,
     8       5599234,  5582132,  5575334,  5567821,  5167153,
     9       5156108,  5145026,  5134288,  5124168,  5077804/
            data (iy(i),i= 151, 200)/
     +       4898398,  4659643,  4419134,  3942617,  3929639,
     1       3923042,  3933236,  3951430,  3975541,  4989411,
     2       5054514,  5112266,  5138064,  5161310,  5181898,
     3       5199488,  5214413,  5226979,  5237375,  5245711,
     4       5251225,  5255710,  5259456,  5262411,  5264493,
     5       5265598,  5265665,  5264710,  5262919,  5260344,
     6       5257038,  5253061,  5248457,  5243203,  5237352,
     7       5230837,  5223771,  5216286,  5208517,  5200622,
     8       5192757,  5186040,  5177227,  2881950,  2592188,
     9       1793022,  1106241,   258165,   210040,   321775/
            data (iy(i),i= 201, 242)/
     +             0,   413863,   443966,   656610,   223114,
     1        139355,   589955,  1421292,  2288002,  3260701,
     2       3797993,  4254527,  4765088,  5251086,  5694451,
     3       6190597,  6807999,  7503138,  8281126,  9070889,
     4       9623615,  9899876,  9228654,  9173282,  9671695,
     5      10000000,  9896622,  9445850,  9858504,  9651270,
     6       9304286,  8609588,  7904434,  7103745,  6408539,
     7       5841843,  5419985,  5018809,  4774224,  4486439,
     8       4057561,  3433861/
c
        sp(2)='airfoil'
        sp(1)='airfoil'
        sp(3)='airfoil'
c
        do i=1,nvf
            vx(i)=real(ix(i),rknd)*1.0e-7_rknd
            vy(i)=real(iy(i),rknd)*1.0e-7_rknd
        enddo
        do i=1,nbf
            ibndry(1,i)=i
            ibndry(2,i)=i+1
            ibndry(3,i)=0
            ibndry(4,i)=2
            if(i>nvf) ibndry(4,i)=0
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=1
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
        ibndry(2,91)=1
        ibndry(2,144)=92
        ibndry(2,193)=145
        ibndry(2,242)=194
c
        ibndry(1,243)=237
        ibndry(2,243)=2
        ibndry(1,244)=61
        ibndry(2,244)=136
        ibndry(1,245)=117
        ibndry(2,245)=172
        ibndry(1,246)=156
        ibndry(2,246)=212
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
c
c       make itnode
c
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd9(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), save, dimension(39) :: ib1,ib2,ib3,ib4
            integer(kind=iknd), save, dimension(30) :: ix,iy
            integer(kind=iknd), save :: nvf,ntf,nbf
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(100) :: rp,ru
            character(len=80), dimension(100) :: sp,su
cy
c       a domain composed of nine subregions
c       data for this problem was supplied by hans mittelmann
c
            external sxy
            data nvf,ntf,nbf/30,9,39/
            data (ix(i),i=   1,  30)/
     +     -3000,-3000,-3000,-3000, 3600, 7600, 4100, 7100, 4700,
     1      7100, 3600, 4100, 4700, 7100, 9000, 9000, 9000, 7600,
     2       575,  775,    0,    0, -575, -775,    0,    0,  500,
     3         0, -500,    0/
            data (iy(i),i=   1,  30)/
     +         0,-4400, 3200, 3600,    0,    0,  500,  500, 1100,
     1      1100, 3200, 3600, 3600, 3600, 3600, 3200,-4400, 3200,
     2         0,    0,  575,  775,    0,    0, -575, -775,    0,
     3       500,    0, -500/
            data ib1/ 4,12,13,14, 3,18, 9, 7, 5, 4, 3,11,12,
     +               13,14,18,15,16,10, 2, 1, 1,24,24,22,20,
     1               20,26,23,21,19,25,19,23,29,27,27,28,19/
            data ib2/12,13,14,15,11,16,10, 8, 6, 3, 1, 5, 7,
     +                9,10, 6,16,17, 8,17, 2,24,23,22,20, 5,
     1               26,24,21,19,25,23,20,29,30,30,28,29,27/
            data ib3/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     +               1,1,0,1,1,1,1,1,1,0,0,1,1,1,1,0/
            data ib4/1,1,1,1,0,0,0,0,0,1,1,0,0,0,0,0,1,1,0,1,1,0,0,0,0,
     +               0,0,0,0,0,0,0,0,0,2,2,2,2,0/
c
        sp(2)='plant'
        sp(1)='plant'
        sp(3)='plant'
c
        do i=1,nbf
            ibndry(1,i)=ib1(i)
            ibndry(2,i)=ib2(i)
            ibndry(3,i)=ib3(i)
            ibndry(4,i)=ib4(i)
            ibndry(5,i)=0
            ibndry(6,i)=0
            if(i<=4) then
                ibndry(7,i)=1
            else
                ibndry(7,i)=0
            endif
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
        do i=1,nvf
            vx(i)=real(ix(i),rknd)*1.0e-2_rknd
            vy(i)=real(iy(i),rknd)*1.0e-2_rknd
        enddo
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
c
c       make itnode, divide long curved edges
c
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        call sklutl(1_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd10(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(100) :: rp,ru,xm,ym
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), save :: hmax,grade
            character(len=80), dimension(100) :: sp,su
cy
c       circular region with many curved interior edges
c       dirichlet boundary conditions on all boundary edges
c
            external sxy
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
c
        sp(2)='fan'
        sp(1)='fan'
        sp(3)='fan'
c
        pi=3.141592653589793e0_rknd
        r=1.0e0_rknd
        r2=r/sqrt(2.0e0_rknd)
        ntf=7
        nvf=ntf+1
        nbf=2*ntf
        vx(1)=0.0e0_rknd
        vy(1)=0.0e0_rknd
        dt=2.0e0_rknd*pi/real(ntf,rknd)
        xm(1)=0.0e0_rknd
        ym(1)=0.0e0_rknd
        do i=1,ntf
            ang=real(i-1,rknd)*dt
            angm=ang-pi/4.0e0_rknd
            vx(i+1)=r*cos(ang)
            vy(i+1)=r*sin(ang)
            xm(i+1)=r2*cos(angm)
            ym(i+1)=r2*sin(angm)
            k=2*i-1
            ibndry(1,k)=1
            ibndry(2,k)=i+1
            ibndry(3,k)=i+1
            ibndry(4,k)=0
            ibndry(5,k)=0
            ibndry(6,k)=0
            ibndry(7,k)=k
            ibndry(1,k+1)=i+1
            ibndry(2,k+1)=i+2
            if(i==ntf) ibndry(2,k+1)=2
            ibndry(3,k+1)=1
            ibndry(4,k+1)=2
            ibndry(5,k+1)=0
            ibndry(6,k+1)=0
            ibndry(7,k+1)=k+1
        enddo
c
        do i=1,nbf
            if(ibndry(3,i)>0) then
                sf(1,i)=xm(ibndry(3,i))
                sf(2,i)=ym(ibndry(3,i))
            else
                sf(1,i)=0.0e0_rknd
                sf(2,i)=0.0e0_rknd
            endif
        enddo
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
c
c       make itnode, refine long edges, find symmetries
c
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        call sklutl(1_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        call sklutl(2_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd11(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), save :: hmax,grade
            character(len=80), dimension(100) :: sp,su
cy
c       domain which is decomposed into several distinct pieces,
c       and making extensive use of linked edges.
c
            external sxy
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
c
        sp(2)='arc'
        sp(1)='arc'
        sp(3)='arc'
c
        pi=3.141592653589793e0_rknd
        r1=1.0e0_rknd
        r2=2.0e0_rknd
        ntf=6
        nvf=4*ntf
        nbf=nvf
        dt=pi/real(ntf,rknd)
        eps=dt/20.0e0_rknd
c
        do i=1,nbf
            ibndry(1,i)=i
            ibndry(2,i)=i+1
            ibndry(3,i)=0
            ibndry(4,i)=2
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=i
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
        do i=1,ntf
            k=(i-1)*4
            t1=real(i-1,rknd)*dt+eps
            t2=real(i,rknd)*dt-eps
            c1=cos(t1)
            c2=cos(t2)
            s1=sin(t1)
            s2=sin(t2)
            vx(k+1)=r1*c1
            vy(k+1)=r1*s1
            vx(k+2)=r2*c1
            vy(k+2)=r2*s1
            vx(k+3)=r2*c2
            vy(k+3)=r2*s2
            vx(k+4)=r1*c2
            vy(k+4)=r1*s2
            ibndry(2,k+4)=k+1
            ibndry(3,k+2)=1
            ibndry(3,k+4)=1
            if(i>1) ibndry(4,k+1)=-(k-1)
            if(i<ntf) ibndry(4,k+3)=-(k+5)
        enddo
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
c
c       make itnode, refine long edges, find symmetries
c
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        call sklutl(1_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        call sklutl(2_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd12(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), save :: ntf,nvf,nbf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru,xm,ym
            real(kind=rknd), save :: hmax,grade
            character(len=80), dimension(100) :: sp,su
cy
c       a spiral shaped domain for testing tgen
c
            external sxy
            data ntf,nvf,nbf/1,0,0/
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
c
        sp(2)='spiral'
        sp(1)='spiral'
        sp(3)='spiral'
c
        nloops=6
        nvf=8*nloops+2
        nbf=nvf
c
        xm(1)=0.0e0_rknd
        ym(1)=0.0e0_rknd
        xm(2)=1.0e0_rknd
        ym(2)=0.0e0_rknd
c
        do k=1,nbf
            ibndry(1,k)=k
            ibndry(2,k)=k+2
            ibndry(3,k)=1
            ibndry(4,k)=2
            ibndry(5,k)=0
            ibndry(6,k)=0
            ibndry(7,k)=k
        enddo
c
        do i=1,nloops
            k=(i-1)*8
            m=2*i+1
            vx(k+1)=-real(m-1,rknd)
            vy(k+1)=0.0e0_rknd
            vx(k+2)=-real(m,rknd)
            vy(k+2)=0.0e0_rknd
            vx(k+3)=1.0e0_rknd
            vy(k+3)=real(m,rknd)
            vx(k+4)=1.0e0_rknd
            vy(k+4)=real(m+1,rknd)
            vx(k+5)=real(m+1,rknd)
            vy(k+5)=0.0e0_rknd
            vx(k+6)=real(m+2,rknd)
            vy(k+6)=0.0e0_rknd
            vx(k+7)=0.0e0_rknd
            vy(k+7)=-real(m+1,rknd)
            vx(k+8)=0.0e0_rknd
            vy(k+8)=-real(m+2,rknd)
            do j=k+1,k+4
                ibndry(3,j)=2
            enddo
        enddo
        m=2*nloops+3
        k=nloops*8
        vx(k+1)=-real(m-1,rknd)
        vy(k+1)=0.0e0_rknd
        vx(k+2)=-real(m,rknd)
        vy(k+2)=0.0e0_rknd
        ibndry(1,k+1)=1
        ibndry(2,k+1)=2
        ibndry(3,k+1)=0
        ibndry(1,k+2)=k+1
        ibndry(2,k+2)=k+2
        ibndry(3,k+2)=0
c
        itnode(1,1)=1
        itnode(2,1)=1
        itnode(3,1)=0
        itnode(5,1)=1
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
c
        do i=1,nbf
            if(ibndry(3,i)>0) then
                sf(1,i)=xm(ibndry(3,i))
                sf(2,i)=ym(ibndry(3,i))
            else
                sf(1,i)=0.0e0_rknd
                sf(2,i)=0.0e0_rknd
            endif
        enddo
c
c       refine long edges
c
        call sklutl(1_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd13(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), save, dimension(135) :: ix,iy
            integer(kind=iknd), save :: ntf,nvf,nbf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), save :: hmax,grade
            character(len=80), dimension(100) :: sp,su
cy
c       this domain is a facsimile of the newly adopted ucsd logo.
c       dirichlet boundary conditions are imposed.
c
            external sxy
            data (ix(i),i=   1,  90)/
     +       320,  520,  520,  450,  460,  320,  250,  585,  585,
     1       460,  470,  250,  165,  658,  658,  470,  480,  165,
     2       100,  680,  680,  600,  550,  520,  490,  470,  460,
     3       100,  460,  500,  530,  550,  563,  570,  560,  540,
     4       510,  280,  310,  316,  320,  319,  316,  310,  100,
     5        40,   68,   68,   75,  100,  200,  220,  235,  235,
     6       190,  215,  215,  205,  190,  120,  105,  100,  390,
     7       390,  370,  310,  280,  255,  245,  245,  255,  270,
     8       300,  380,  390,  330,  300,  285,  270,  270,  280,
     9       295,  310,  350,  380,  545,  545,  525,  445,  420/
            data (ix(i),i=  91, 135)/
     +       410,  410,  420,  500,  520,  500,  420,  420,  510,
     1       540,  555,  555,  545,  500,  445,  435,  450,  505,
     2       535,  610,  550,  578,  578,  610,  610,  670,  700,
     3       715,  715,  700,  670,  620,  620,  690,  725,  737,
     4       747,  747,  737,  725,  690, -332, 1147, 1147, -332/
            data (iy(i),i=   1,  90)/
     +       520,  520,  490,  490,  500,  500,  482,  482,  450,
     1       450,  460,  460,  440,  440,  410,  410,  420,  420,
     2       400,  400,  370,  360,  350,  345,  330,  310,  290,
     3       380,  270,  260,  250,  240,  230,  219,  210,  207,
     4       210,  280,  281,  283,  288,  293,  296,  300,  150,
     5       150,  140,   30,   10,    0,    0,   13,   25,  150,
     6       150,  140,   30,   20,   10,   10,   20,   30,  110,
     7       145,  150,  150,  145,  130,  100,   50,   30,   15,
     8         0,    0,   10,   10,   20,   30,   50,  100,  120,
     9       130,  140,  140,  115,  110,  145,  150,  150,  140/
            data (iy(i),i=  91, 135)/
     +       130,  100,   85,   50,   30,   10,   10,    0,    0,
     1        10,   25,   50,   60,   80,  100,  120,  140,  140,
     2       115,  150,  150,  140,    0,    0,  140,  140,  130,
     3       105,   45,   20,   10,   10,    0,    0,   15,   30,
     4        50,  100,  120,  135,  150, -400, -400,  920,  920/
c
            data ntf,nvf,nbf/9,135,145/
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
c
        sp(2)='ucsd'
        sp(1)='ucsd'
        sp(3)='ucsd'
c
        do i=1,nvf
            vx(i)=real(ix(i),rknd)
            vy(i)=real(iy(i),rknd)
            ibndry(1,i)=i
            ibndry(2,i)=i+1
            ibndry(3,i)=0
            ibndry(4,i)=2
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=i
        enddo
        ibndry(2,6)=1
        ibndry(2,12)=7
        ibndry(2,18)=13
        ibndry(2,27)=19
        ibndry(2,44)=28
        ibndry(2,62)=45
        ibndry(2,85)=63
        ibndry(2,109)=86
        ibndry(2,131)=110
        ibndry(2,135)=132
c
        do i=nvf+1,nbf
            ibndry(3,i)=0
            ibndry(4,i)=0
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=i
        enddo
        ibndry(1,136)=132
        ibndry(2,136)=49
        ibndry(1,137)=54
        ibndry(2,137)=68
        ibndry(1,138)=75
        ibndry(2,138)=97
        ibndry(1,139)=100
        ibndry(2,139)=113
        ibndry(1,140)=110
        ibndry(2,140)=35
        ibndry(1,141)=29
        ibndry(2,141)=27
        ibndry(1,142)=19
        ibndry(2,142)=18
        ibndry(1,143)=13
        ibndry(2,143)=12
        ibndry(1,144)=6
        ibndry(2,144)=7
        ibndry(1,145)=1
        ibndry(2,145)=135
c
        do i=1,nbf
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
c
c       make itnode
c
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd14(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), save :: ntf,nvf,nbf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(16) :: xx,yy
            character(len=80), dimension(100) :: sp,su
cy
            external sxy
            data ntf,nvf,nbf/4,29,32/
c
c        a domain in the shape of a nozzle.
c        four similar regions, several curved edges.
c
        sp(2)='nozzle'
        sp(1)='nozzle'
        sp(3)='nozzle'
c
        rp(1)=ru(8)
c
        vx(1)=0.0e0_rknd
        vy(1)=0.0e0_rknd
        vx(2)=4.0e0_rknd
        vy(2)=0.0e0_rknd
        vx(3)=4.0e0_rknd
        vy(3)=1.0e0_rknd
        vx(4)=3.0e0_rknd
        vy(4)=1.0e0_rknd
        vx(5)=2.5e0_rknd
        vy(5)=f(2.5e0_rknd)
        vx(6)=2.0e0_rknd
        vy(6)=f(2.0e0_rknd)
        vx(7)=1.5e0_rknd
        vy(7)=f(1.5e0_rknd)
        vx(8)=1.0e0_rknd
        vy(8)=0.5e0_rknd
        vx(9)=0.0e0_rknd
        vy(9)=0.5e0_rknd
        do i=10,16
            vx(i)=-vx(18-i)
            vy(i)=vy(18-i)
        enddo
        do i=17,29
            vx(i)=vx(32-i)
            vy(i)=-vy(32-i)
        enddo
c
        xx(1)=2.75e0_rknd
        yy(1)=f(2.75e0_rknd)
        xx(2)=2.25e0_rknd
        yy(2)=f(2.25e0_rknd)
        xx(3)=1.75e0_rknd
        yy(3)=f(1.75e0_rknd)
        xx(4)=1.25e0_rknd
        yy(4)=f(1.25e0_rknd)
        do i=1,4
            xx(4+i)=-xx(5-i)
            yy(4+i)=yy(5-i)
            xx(8+i)=-xx(i)
            yy(8+i)=-yy(i)
            xx(12+i)=xx(5-i)
            yy(12+i)=-yy(5-i)
        enddo
c
        do i=1,29
            ibndry(1,i)=i
            ibndry(2,i)=i+1
            ibndry(3,i)=0
            ibndry(4,i)=1
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=i
        enddo
        do i=1,4
            ibndry(3,3+i)=i
            ibndry(3,9+i)=i+4
            ibndry(3,17+i)=i+8
            ibndry(3,23+i)=i+12
        enddo
        ibndry(4,1)=0
        ibndry(4,2)=2
        ibndry(4,15)=2
        ibndry(4,16)=2
        ibndry(4,29)=2
        ibndry(2,29)=2
        do i=30,32
            ibndry(1,i)=1
            ibndry(2,i)=9
            ibndry(3,i)=0
            ibndry(4,i)=0
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=i
        enddo
        ibndry(2,30)=16
        ibndry(2,32)=23
c
        do i=1,nbf
            if(ibndry(3,i)>0) then
                i1=ibndry(1,i)
                i2=ibndry(2,i)
                i3=ibndry(3,i)
                call centre(vx(i1),vy(i1),vx(i2),vy(i2),
     +              xx(i3),yy(i3),sf(1,i),sf(2,i))
            else
                sf(1,i)=0.0e0_rknd
                sf(2,i)=0.0e0_rknd
            endif
        enddo
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
c
c       make itnode, find symmetries
c
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        call sklutl(2_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        function f(x)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            real(kind=rknd) :: f
cy
        f=-0.25e0_rknd*x+1.5e0_rknd
        f=f*x-2.25e0_rknd
        f=(f*x+2.0e0_rknd)/2.0e0_rknd
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd15(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), save, dimension(170) :: ix,iy
            integer(kind=iknd), save :: ntf,nvf,nbf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), save :: hmax,grade
            character(len=80), dimension(100) :: sp,su
cy
            external sxy
            data (ix(i),i=   1,  90)/
     +         0,   40,   80,  110,  130,  140,   90,   60,  100,
     1       180,  230,  210,  260,  300,  350,  380,  370,  330,
     2       400,  430,  425,  420,  430,  455,  490,  525,  560,
     3       580,  590,  585,  570,  540,  500,  530,  570,  610,
     4       630,  580,  610,  650,  690,  730,  715,  775,  755,
     5       700,  700,  710,  740,  760,  790,  815,  830,  890,
     6       865,  900,  950,  975, 1000, 1040, 1075, 1085, 1100,
     7      1130, 1115, 1115, 1130, 1170, 1200, 1230, 1250, 1250,
     8      1300, 1310, 1300, 1290, 1420, 1360, 1270, 1240, 1170,
     9      1100, 1000,  900,  800,  700,  600,  500,  400,  350/
            data (ix(i),i=  91, 170)/
     +       320,  310,  320,  350,  400,  500,  600,  700,  800,
     1       900, 1000, 1100, 1200, 1160, 1100, 1050, 1025, 1020,
     2      1010, 1000,  980,  950,  895,  905,  915,  905,  875,
     3       840,  780,  790,  750,  700,  660,  620,  570,  520,
     4       470,  440,  400,  350,  300,  240,  180,  150,  110,
     5        60,  190,  220,  280,  310,  280,  220,  480,  485,
     6       520,  550,  550,  520, 1170, 1200, 1240, 1260, 1220,
     7      1175, 1180, 1100, 1000,  900,  800,  700,  600,  500,
     8       460,  500,  600,  700,  800,  900, 1000, 1100/
            data (iy(i),i=   1,  90)/
     +        70,  100,  160,  200,  260,  300,  270,  310,  340,
     1       380,  380,  340,  355,  370,  360,  330,  280,  240,
     2       270,  300,  340,  370,  400,  450,  500,  545,  575,
     3       550,  505,  450,  400,  350,  300,  265,  280,  320,
     4       350,  350,  390,  400,  480,  500,  410,  420,  360,
     5       350,  300,  260,  250,  270,  300,  340,  390,  400,
     6       360,  385,  400,  380,  390,  390,  370,  300,  270,
     7       260,  300,  320,  360,  380,  385,  380,  365,  380,
     8       390,  300,  240,  200,  220,  150,  130,  100,   50,
     9        10,   -5,  -15,  -10,   -5,    0,   15,   45,   70/
            data (iy(i),i=  91, 170)/
     +       100,  120,  140,  150,  160,  161,  162,  163,  165,
     1       166,  167,  168,  170,  215,  200,  200,  220,  270,
     2       315,  330,  300,  210,  200,  240,  290,  310,  290,
     3       210,  200,  220,  210,  200,  220,  245,  210,  200,
     4       220,  250,  225,  210,  200,  190,  200,  150,  100,
     5        60,  260,  290,  320,  310,  280,  250,  350,  400,
     6       460,  500,  460,  400,  330,  340,  310,  270,  260,
     7       300,  130,   70,   50,   45,   50,   60,   70,   80,
     8        95,  110,  120,  122,  124,  126,  127,  129/
c
            data ntf,nvf,nbf/5,170,178/
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
c
        sp(2)='pltmg'
        sp(1)='pltmg'
        sp(3)='pltmg'
c
        do i=1,nvf
            vx(i)=real(ix(i),rknd)
            vy(i)=real(iy(i),rknd)
            ibndry(1,i)=i
            ibndry(2,i)=i+1
            ibndry(3,i)=0
            ibndry(4,i)=2
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=i
        enddo
        ibndry(2,136)=1
        ibndry(2,142)=137
        ibndry(2,148)=143
        ibndry(2,154)=149
        ibndry(2,170)=155
        do i=nvf+1,nbf
            ibndry(3,i)=0
            ibndry(4,i)=0
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=i
        enddo
        ibndry(1,171)=5
        ibndry(2,171)=137
        ibndry(1,172)=16
        ibndry(2,172)=140
        ibndry(1,173)=20
        ibndry(2,173)=143
        ibndry(1,174)=27
        ibndry(2,174)=146
        ibndry(1,175)=67
        ibndry(2,175)=149
        ibndry(1,176)=75
        ibndry(2,176)=152
        ibndry(1,177)=102
        ibndry(2,177)=170
        ibndry(1,178)=82
        ibndry(2,178)=156
c
        do i=1,nbf
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
c
c       make itnode
c
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd16(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), save :: hmax,grade
            character(len=80), dimension(100) :: sp,su
cy
c       domain which is has many cracks
c
            external sxy
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
c
        sp(2)='crack'
        sp(1)='crack'
        sp(3)='crack'
c
        pi=3.141592653589793e0_rknd
        r1=1.0e0_rknd
        r2=2.0e0_rknd
        ncrack=8
        nvf=1+3*ncrack
        nbf=4*ncrack
        ntf=ncrack
c
        vx(1)=0.0e0_rknd
        vy(1)=0.0e0_rknd
        do i=1,ncrack
            ang=2.0e0_rknd*pi*real(i-1,rknd)/real(ncrack,rknd)
            c=cos(ang)
            s=sin(ang)
            ii=3*i
            vx(ii-1)=r2*c
            vy(ii-1)=r2*s
            vx(ii)=r1*c
            vy(ii)=r1*s
            vx(ii+1)=r2*c
            vy(ii+1)=r2*s
c
            jj=4*i-3
c
            ibndry(1,jj)=1
            ibndry(2,jj)=ii
            ibndry(3,jj)=0
            ibndry(4,jj)=0
            ibndry(5,jj)=0
            ibndry(6,jj)=0
            ibndry(7,jj)=0
c
            ibndry(1,jj+1)=ii
            ibndry(2,jj+1)=ii-1
            ibndry(3,jj+1)=0
            ibndry(4,jj+1)=1
            ibndry(5,jj+1)=0
            ibndry(6,jj+1)=0
            ibndry(7,jj+1)=1
c
            ibndry(1,jj+2)=ii
            ibndry(2,jj+2)=ii+1
            ibndry(3,jj+2)=0
            ibndry(4,jj+2)=2
            ibndry(5,jj+2)=0
            ibndry(6,jj+2)=0
            ibndry(7,jj+2)=2
c
            ibndry(1,jj+3)=ii+1
            ibndry(2,jj+3)=ii+2
            ibndry(3,jj+3)=1
            ibndry(4,jj+3)=2
            ibndry(5,jj+3)=0
            ibndry(6,jj+3)=0
            ibndry(7,jj+3)=3
        enddo
        ibndry(2,nbf)=2
c
        do i=1,nbf
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
c
c       make itnode, refine long arcs, find symmetries
c
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        call sklutl(1_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        call sklutl(2_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd17(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), dimension(1250) :: list
            integer(kind=iknd), save, dimension(300) :: ix,iy
            integer(kind=iknd), save :: ntf,nvf,nbf
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), save :: hmax,grade
            character(len=80), dimension(100) :: sp,su
cy
c       region in the shape of monterey bay
c       the data was provided by francios lekien
c
            external sxy
            data ntf,nvf,nbf/1,300,300/
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
c
            data (ix(i),i=   1,  50)/
     +      -122324501, -122324501, -122135002, -122135002, -122134003,
     1      -122133003, -122128998, -122127998, -122125999, -122125000,
     2      -122122002, -122120003, -122120003, -122116997, -122115997,
     3      -122113998, -122113998, -122111000, -122109001, -122107002,
     4      -122105003, -122102997, -122100998, -122098999, -122097000,
     5      -122094002, -122089996, -122088997, -122086998, -122084999,
     6      -122083000, -122081001, -122080002, -122078003, -122075996,
     7      -122075996, -122073997, -122070999, -122069000, -122065002,
     8      -122060997, -122059998, -122057999, -122056000, -122052002,
     9      -122050003, -122046997, -122044998, -122042000, -122040001/
            data (ix(i),i=  51, 100)/
     +      -122039001, -122035004, -122033997, -122030998, -122028000,
     1      -122024002, -122024002, -122026001, -122026001, -122025002,
     2      -122022003, -122016998, -122014000, -122012001, -122008003,
     3      -122004997, -122001999, -122001999, -122000000, -121997002,
     4      -121992996, -121987999, -121985001, -121980003, -121977997,
     5      -121975998, -121973999, -121973000, -121970001, -121968002,
     6      -121964996, -121961998, -121959999, -121958000, -121955002,
     7      -121953003, -121947998, -121945000, -121941002, -121939003,
     8      -121936996, -121933998, -121931999, -121929001, -121925003,
     9      -121919998, -121914001, -121911003, -121910004, -121906998/
            data (ix(i),i= 101, 150)/
     +      -121903000, -121900002, -121897003, -121892998, -121891998,
     1      -121889000, -121887001, -121884003, -121883003, -121882004,
     2      -121879997, -121876999, -121876999, -121875000, -121872002,
     3      -121869003, -121865997, -121862999, -121859001, -121857002,
     4      -121855003, -121852997, -121849998, -121848000, -121847000,
     5      -121843002, -121838997, -121835999, -121834000, -121830002,
     6      -121828003, -121825996, -121822998, -121820999, -121819000,
     7      -121816002, -121814003, -121811996, -121808998, -121807999,
     8      -121807999, -121806000, -121805000, -121804001, -121804001,
     9      -121803001, -121800003, -121799004, -121799004, -121796997/
            data (ix(i),i= 151, 200)/
     +      -121793999, -121792000, -121790001, -121790001, -121789001,
     1      -121789001, -121788002, -121788002, -121789001, -121791000,
     2      -121792999, -121793999, -121794998, -121794998, -121795998,
     3      -121796997, -121799004, -121800003, -121800003, -121801003,
     4      -121802002, -121804001, -121805000, -121805000, -121804001,
     5      -121805000, -121806999, -121806999, -121806999, -121806999,
     6      -121806999, -121807999, -121808998, -121810997, -121813004,
     7      -121815002, -121817001, -121819000, -121821999, -121824997,
     8      -121829002, -121832001, -121834999, -121838997, -121842003,
     9      -121847000, -121851997, -121859001, -121862999, -121867996/
            data (ix(i),i= 201, 250)/
     +      -121873001, -121876999, -121882004, -121886002, -121889000,
     1      -121892998, -121893997, -121893997, -121896004, -121897003,
     2      -121898003, -121900002, -121902000, -121903999, -121903999,
     3      -121904999, -121906998, -121908997, -121910004, -121910004,
     4      -121911003, -121913002, -121914001, -121915001, -121916000,
     5      -121915001, -121915001, -121916000, -121917999, -121919998,
     6      -121919998, -121921997, -121921997, -121922997, -121924004,
     7      -121926003, -121928001, -121930000, -121931000, -121931999,
     8      -121934998, -121936996, -121936996, -121935997, -121936996,
     9      -121936996, -121939003, -121939003, -121939003, -121941002/
            data (ix(i),i= 251, 300)/
     +      -121941002, -121942001, -121941002, -121941002, -121942001,
     1      -121944000, -121947998, -121950996, -121950996, -121953003,
     2      -121954002, -121956001, -121958000, -121958000, -121959000,
     3      -121959999, -121959999, -121960999, -121962997, -121960999,
     4      -121961998, -121961998, -121963997, -121964996, -121963997,
     5      -121963997, -121966003, -121968002, -121970001, -121973999,
     6      -121974998, -121974998, -121974998, -121975998, -121975998,
     7      -121977997, -121977997, -121975998, -121974998, -121973000,
     8      -121972000, -121971001, -121970001, -121969002, -121967003,
     9      -121963997, -121962997, -121959999, -121958000, -121957001/
            data (iy(i),i=   1,  50)/
     +        36565899,   36973900,   36973900,   36973202,   36971298,
     1        36972500,   36972099,   36971100,   36970001,   36970699,
     2        36970699,   36969799,   36967499,   36967499,   36967899,
     3        36967201,   36965900,   36965500,   36963902,   36962502,
     4        36961399,   36961399,   36960701,   36960499,   36960499,
     5        36959801,   36959099,   36959099,   36957901,   36958401,
     6        36957298,   36957901,   36958599,   36957901,   36957500,
     7        36956402,   36956402,   36954700,   36953800,   36953800,
     8        36953899,   36954700,   36954498,   36953400,   36954498,
     9        36954498,   36954498,   36955200,   36955200,   36957100/
            data (iy(i),i=  51, 100)/
     +        36957901,   36957298,   36956402,   36956799,   36956402,
     1        36955700,   36958801,   36961300,   36963902,   36966599,
     2        36967899,   36967899,   36968601,   36968601,   36967999,
     3        36967999,   36966801,   36965199,   36966099,   36966099,
     4        36964500,   36962502,   36962502,   36962502,   36960701,
     5        36958801,   36958801,   36958801,   36960499,   36962700,
     6        36964600,   36965500,   36967999,   36970501,   36973202,
     7        36975201,   36976601,   36978001,   36980202,   36981998,
     8        36981998,   36981998,   36982700,   36982201,   36980499,
     9        36979099,   36977299,   36975399,   36973900,   36973202/
            data (iy(i),i= 101, 150)/
     +        36972099,   36970699,   36969101,   36967300,   36965500,
     1        36963200,   36961102,   36958801,   36957500,   36957298,
     2        36955700,   36953400,   36952000,   36949100,   36947300,
     3        36944099,   36941101,   36936100,   36932201,   36928799,
     4        36924801,   36921398,   36917999,   36914799,   36911800,
     5        36907700,   36902699,   36897999,   36893799,   36889099,
     6        36884499,   36880600,   36875900,   36872299,   36868599,
     7        36863899,   36860001,   36856998,   36855400,   36852299,
     8        36850700,   36848400,   36845001,   36842499,   36839901,
     9        36837299,   36835400,   36833099,   36829800,   36826500/
            data (iy(i),i= 151, 200)/
     +        36823200,   36820400,   36817299,   36814701,   36811298,
     1        36809898,   36807400,   36804501,   36800900,   36797699,
     2        36793598,   36790401,   36786701,   36785000,   36781101,
     3        36777199,   36772400,   36769001,   36765202,   36761799,
     4        36760201,   36757401,   36753601,   36749699,   36748402,
     5        36745899,   36741600,   36737202,   36732899,   36727200,
     6        36718300,   36709499,   36701302,   36695900,   36689201,
     7        36683601,   36676998,   36670601,   36665600,   36658798,
     8        36651901,   36648300,   36643600,   36637199,   36632198,
     9        36626099,   36620800,   36614300,   36611500,   36609001/
            data (iy(i),i= 201, 250)/
     +        36607201,   36606098,   36604900,   36604900,   36605400,
     1        36605801,   36606998,   36611000,   36612900,   36615799,
     2        36617699,   36619499,   36621498,   36623798,   36624500,
     3        36624199,   36624500,   36623600,   36623600,   36624901,
     4        36626099,   36626099,   36626999,   36627800,   36628502,
     5        36629200,   36629700,   36629700,   36630402,   36631500,
     6        36632900,   36633801,   36635799,   36636799,   36637901,
     7        36638802,   36639900,   36640400,   36640400,   36640400,
     8        36640400,   36637600,   36635101,   36632198,   36631001,
     9        36629902,   36629902,   36628101,   36626900,   36626900/
            data (iy(i),i= 251, 300)/
     +        36624500,   36623100,   36621700,   36621101,   36620399,
     1        36617699,   36614300,   36612900,   36611900,   36611900,
     2        36612598,   36614201,   36612598,   36611099,   36609501,
     3        36609200,   36607601,   36605099,   36603100,   36600800,
     4        36599201,   36597198,   36596100,   36595600,   36593498,
     5        36590099,   36588100,   36585800,   36585800,   36585800,
     6        36584900,   36583599,   36582401,   36582401,   36583099,
     7        36583099,   36581902,   36580799,   36578999,   36578098,
     8        36576099,   36574402,   36573101,   36571899,   36571201,
     9        36570301,   36568501,   36567799,   36567402,   36565899/
c
c
        sp(2)='monterey bay'
        sp(1)='monterey bay'
        sp(3)='monterey bay'
c
        do i=1,nvf
            ibndry(1,i)=i-1
            ibndry(2,i)=i
            ibndry(3,i)=0
            ibndry(4,i)=2
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=0
            vx(i)=real(ix(i),rknd)*1.0e-6_rknd
            vy(i)=real(iy(i),rknd)*1.0e-6_rknd
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
        ibndry(1,1)=nvf
c
c       make itnode
c
        itnode(1,1)=1
        itnode(2,1)=1
        itnode(3,1)=0
        itnode(4,1)=0
        itnode(5,1)=1
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd18(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), save, dimension(89) :: ix,iy
            integer(kind=iknd), save :: ntf,nvf,nbf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), save :: hmax,grade
            character(len=80), dimension(100) :: sp,su
cy
c       siam logo
c
            external sxy
            data ntf,nvf,nbf/2,89,92/
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
c
            data (ix(i),i=   1,  89)/
     +         0,   30,   32,   34,   57,   59,   61,   89,   89,
     1        91,  107,  109,  109,  114,  116,  116,  125,  125,
     2       130,  132,  132,  134,  158,  158,  141,  139,  139,
     3       137,  127,  125,  123,  111,  109,  109,  102,  102,
     4       100,   91,   89,   89,   87,   61,   61,   82,   82,
     5        81,   67,   66,   66,   67,   80,   80,   60,   58,
     6        58,   56,   50,   48,   48,   41,   41,   39,   35,
     7        33,   33,   31,    9,    8,    8,    9,   33,   33,
     8         3,    1,    1,    3,   25,   26,   26,   25,    0,
     9        48,   50,   39,   41,  -60,  220,  220,  -60/
            data (iy(i),i=   1,  89)/
     +         0,    0,    2,    0,    0,    2,    0,    0,    2,
     1         0,    0,    2,   26,   26,   24,    0,    0,   26,
     2        26,   24,    2,    0,    0,    7,    7,    9,   31,
     3        33,   33,   31,   33,   33,   31,   33,   33,    9,
     4         7,    7,    9,   31,   33,   33,   27,   27,    8,
     5         7,    7,    8,   13,   14,   14,   21,   21,   19,
     6         9,    7,    7,    9,   33,   33,    9,    7,    7,
     7         9,   19,   21,   21,   22,   26,   27,   27,   33,
     8        33,   31,   16,   14,   14,   13,    8,    7,    7,
     9        40,   51,   51,   40,  -90,  -90,  140,  140/
c
        sp(2)='siam'
        sp(1)='siam'
        sp(3)='siam'
c
        do i=1,nvf
            vx(i)=real(ix(i),rknd)
            vy(i)=real(iy(i),rknd)
        enddo
        do i=1,nbf
            ibndry(1,i)=i
            ibndry(2,i)=i+1
            ibndry(3,i)=0
            ibndry(4,i)=2
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=0
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
        ibndry(2,81)=1
        ibndry(2,85)=82
        ibndry(2,89)=86
        ibndry(1,90)=86
        ibndry(2,90)=1
        ibndry(1,91)=60
        ibndry(2,91)=85
        ibndry(1,92)=84
        ibndry(2,92)=89
        do i=nvf+1,nbf
            ibndry(4,i)=0
        enddo
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
c
c       make itnode
c
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd19(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), dimension(1000) :: b1,b2
            integer(kind=iknd), dimension(968) :: ix,iy
            integer(kind=iknd), save :: ntf,nvf,nbf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), save :: hmax,grade
            character(len=80), dimension(100) :: sp,su
cy
c       domain in the shape of mexico
c
            external sxy
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
c
            data ntf,nvf,nbf/1,968,999/
            data nbound/  427/
c
            data (ix(i),i=   1,  80)/
     +      10555, 10555, 10699, 11121, 11236, 11562, 11677, 11677,
     1      11792, 11869, 12032, 12032, 11955, 12003, 12358, 12378,
     2      12493, 12445, 12608, 12646, 12531, 12838, 13088, 13155,
     3      13472, 13587, 13481, 13596, 14057, 14172, 14191, 14373,
     4      14776, 14959, 15026, 15170, 15573, 15410, 15007, 14796,
     5      14661, 14383, 14162, 13855, 13529, 12579, 11860, 11764,
     6      11716, 11600, 11437, 11303, 11140, 11351, 10948, 10699,
     7       9355,  9211,  8962,  8616,  8405,  8146,  7686,  7590,
     8       6870,  6937,  7446,  7858,  8041,  8252,  8530,  8578,
     9       8501,  8300,  8070,  8271,  8041,  7954,  8261, 13951/
            data (ix(i),i=  81, 160)/
     +      17396, 16168, 16129, 15832, 15582, 15285, 15131, 14767,
     1      14508, 14182, 14066, 14134, 13731, 13644, 13596, 13337,
     2      13040, 12790, 12224, 12042, 12042, 11668, 11620, 11389,
     3      11399, 10862, 10842, 10670, 10411, 10497, 10497, 10200,
     4       9989,  9451,  9221,  8943,  8396, 26300, 18307, 18394,
     5      18413, 21138, 21570, 21944, 22548, 22961, 23642, 23825,
     6      24016, 24045, 24419, 24698, 25062, 25341, 25696, 25907,
     7      32057, 31635, 31280, 31213, 30599, 29975, 29226, 28651,
     8      27960, 27749, 27519, 27231, 27020, 26549,  4299,  8146,
     9       8002,  7839,  7619,  7667,  7782,  7887,  8146,  8386/
            data (ix(i),i= 161, 240)/
     +       9048,  9451,  9480,  9825, 10036, 10411, 10488, 10555,
     1       8319,  8377,  8540,  8070,  7724,  6976,  6659,  6045,
     2       6045,  5911,  5920,  5604,  5546,  5268,  5239,  4961,
     3       4846,  4826,  4970,  4855,  4654,  4529,  4270,  4318,
     4       4289, 32162, 36452, 36576, 36288, 36010, 35751, 35703,
     5      35924, 36144, 36039, 35962, 35761, 35598, 35665, 35598,
     6      35780, 35655, 35674, 35540, 35559, 35722, 35636, 35626,
     7      35684, 35684, 35463, 35588, 32441, 32565, 32498, 32681,
     8      32662, 33055, 33276, 33755, 34034, 34542, 34907, 35741,
     9      36116, 42477, 41988, 41614, 41220, 40942, 40549, 39464/
            data (ix(i),i= 241, 320)/
     +      39186, 38649, 38227, 37622, 36730, 36336, 36260, 35991,
     1      35828, 36269, 36308, 36317, 36423, 36471, 36308, 35790,
     2      22107, 21877, 21570, 21378, 21426, 20361, 19843, 19660,
     3      19171, 18893, 18586, 18240, 18211, 18077, 18077, 17617,
     4      17281, 16954, 16801, 16686, 16081, 15995, 15995, 15822,
     5      15938, 16187, 22663, 22318, 22184, 22270, 22596, 22702,
     6      22769, 22376, 22021, 22059, 22280, 22347, 22078, 23738,
     7      22568, 22232, 21915, 21982, 22193, 23153, 28075, 34600,
     8      34139, 33794, 33122, 32345, 31942, 31529, 29965, 28939,
     9      28468, 28219, 25197, 25705, 27912, 24026, 24986, 35166/
            data (ix(i),i= 321, 400)/
     +      35626, 36010, 37210, 37305, 37651, 38236, 38601, 40232,
     1      41028, 41230, 41441, 41412, 41719, 42065, 42151, 42708,
     2      43293, 44022, 44166, 44463, 45020, 45336, 48157, 48186,
     3      47524, 48724, 48618, 48762, 48378, 48378, 48244, 47467,
     4      47217, 47505, 47476, 47064, 46488, 46421, 46190, 46286,
     5      46171, 50921, 48954, 49846, 50028, 51708, 52072, 53713,
     6      47476, 48129, 48752, 48915, 49290, 49318, 46968, 46152,
     7      46430, 46363, 46325, 46325, 45317, 45116, 43801, 43111,
     8      54126, 54241, 54289, 54193, 54442, 54701, 54999, 54932,
     9      54999, 54845, 54164, 53991, 53991, 53780, 53665, 53627/
            data (ix(i),i= 401, 480)/
     +      53780, 54087, 54173, 54116, 53838, 53751, 53790, 54068,
     1      54097, 53982, 53608, 53579, 53511, 53521, 53224, 53271,
     2      53204, 53157, 53089, 52792, 52408, 52360, 51938, 51554,
     3      51381, 51362, 51055,  7993,  7580,  7532,  7667, 17281,
     4      17434, 17492, 17665, 17559, 17396, 17434, 17511, 17415,
     5      16820, 16609, 16830, 17175, 17290, 17300, 17549, 16964,
     6      25168, 25600, 24726, 24112, 23853, 23134, 22798, 22213,
     7      22078, 21618, 21455, 21253, 21052, 20936, 20975, 20821,
     8      20754, 20457, 19833, 19344, 19065, 18653, 18576, 18365,
     9      18346, 18240, 17847, 26233, 26089, 26099, 25868, 25907/
            data (ix(i),i= 481, 560)/
     +      26252, 26252, 26952, 27240, 27605, 28344, 28977, 29025,
     1      29409, 29716, 30196, 30378, 30445, 30531, 30397, 30435,
     2      30992, 31088, 30934, 30426, 30359, 30071, 29917, 29754,
     3      30368, 30503, 30752, 30848, 30733, 30551, 30426, 30474,
     4      31002, 31126, 31405, 31702, 31894, 31942, 32172, 32201,
     5      32211, 32441, 32479, 32786, 32815, 33093, 33228, 33391,
     6      33688, 34034, 34149, 34158, 34369, 34197, 33621, 33458,
     7      33295, 33208, 32978, 32450, 32335, 32518, 32489, 32690,
     8      32383, 32105, 31990, 31990, 31760, 31654, 31424, 31194,
     9      31050, 31069, 30829, 30848, 34859, 34676, 34312, 33784/
            data (ix(i),i= 561, 640)/
     +      33160, 32930, 32748, 32537, 31740, 31952, 31952, 34629,
     1      34619, 34369, 34484, 34283, 34417, 34446, 34609, 34724,
     2      34878, 35079, 34964, 34638, 34475, 34542, 34465, 34513,
     3      34580, 35060, 35214, 35329, 35291, 35454, 35578, 35722,
     4      35838, 35857, 36192, 36125, 35847, 35847, 36125, 36471,
     5      36835, 36902, 36663, 36663, 36432, 36691, 37066, 37104,
     6      36787, 36931, 36701, 36701, 37075, 37114, 37459, 37555,
     7      37737, 37881, 37900, 38198, 38783, 39311, 39455, 39320,
     8      39368, 39628, 40501, 40587, 40568, 40894, 40971, 41633,
     9      43072, 43379, 43475, 43379, 43379, 43082, 42813, 42583/
            data (ix(i),i= 641, 720)/
     +      42612, 34264, 33919, 33755, 33477, 33343, 33132, 32738,
     1      32374, 32191, 31107, 30762, 30551, 30081, 29831, 29601,
     2      29198, 29457, 29447, 29543, 29447, 29198, 28871, 28785,
     3      28018, 27989, 27874, 28075, 28209, 28881, 29629, 29879,
     4      27499, 27451, 26147, 25235, 25197, 24803, 24650, 24563,
     5      24381, 24362, 24986, 24870, 24995, 25187, 25072, 25168,
     6      25398, 25494, 25475, 25590, 25955, 26089, 26415, 26434,
     7      25686, 25734, 25648, 25369, 25350, 25513, 26329, 26463,
     8      26501, 26895, 27192, 27240, 27077, 26924, 27116, 27269,
     9      27317, 27231, 27260, 27538, 27768, 28229, 28286, 28526/
            data (ix(i),i= 721, 800)/
     +      28564, 29188, 24199, 24016, 23738, 23738, 23489, 23278,
     1      23191, 23249, 23335, 23067, 22443, 21973, 21925, 21685,
     2      21618, 21080, 20783, 20572, 20428, 19919, 19689, 22347,
     3      22376, 22165, 22337, 22289, 24237, 24371, 24621, 24544,
     4      24890, 24544, 24448, 24343, 24323, 23834, 23489, 22875,
     5      24026, 24208, 24640, 24832, 25091, 25388, 25638, 25705,
     6      25830, 26051, 26147, 26425, 26655, 27231, 27317, 27250,
     7      27048, 26876, 27010, 26962, 26665, 26501, 26482, 26607,
     8      26933, 27164, 27557, 28200, 28363, 28209, 28622, 28891,
     9      29083, 29006, 29121, 29092, 27950, 27595, 27317, 32326/
            data (ix(i),i= 801, 880)/
     +      32230, 31884, 31808, 31597, 31069, 30819, 30982, 31309,
     1      31155, 30954, 30675, 30397, 30215, 29927, 29543, 29418,
     2      29303, 29121, 28910, 28612, 28430, 28411, 31654, 31932,
     3      31952, 32162, 32211, 32211, 32565, 32633, 32633, 33007,
     4      33084, 33525, 35127, 35223, 35022, 35099, 34820, 34753,
     5      34302, 34312, 34197, 33871, 33717, 33717, 33602, 33208,
     6      33064, 37632, 37200, 36807, 36480, 36317, 36087, 35972,
     7      35905, 35933, 35770, 35540, 34916, 34705, 34178, 33669,
     8      33650, 33871, 34264, 34197, 34283, 34235, 34293, 34302,
     9      34321, 34571, 34897, 35147, 35406, 35588, 36116, 36029/
            data (ix(i),i= 881, 960)/
     +      35703, 35694, 35559, 34273, 33256, 33103, 33218, 33496,
     1      33746, 33679, 33554, 32738, 32930, 33113, 33276, 33391,
     2      33909, 33055, 33064, 32441, 32162, 31865, 31443, 31299,
     3      30973, 30877, 31462, 31625, 31664, 31760, 28190, 28651,
     4      28574, 28785, 29313, 29908, 30215, 30301, 30723, 30531,
     5      30579, 30627, 30714, 34859, 35070, 35166, 35463, 35358,
     6      35079, 35012, 35041, 25513, 42775, 42612, 42717, 42928,
     7      43014, 45519, 45864, 45893, 46008, 46085, 46546, 46891,
     8      47073, 47467, 48090, 47179, 47035, 46910, 46757, 46728,
     9      46450, 46267, 46037, 45605, 44991, 44531, 44569, 44060/
            data (ix(i),i= 961, 968)/
     +      43840, 50787, 50700, 49558, 49376, 48896, 53723, 53281/
            data (iy(i),i=   1,  80)/
     +      46567, 46452, 45905, 45454, 45080, 44763, 44715, 44437,
     1      44188, 43842, 43679, 43957, 44159, 44236, 43871, 43622,
     2      43506, 43363, 43017, 42672, 42441, 41799, 41626, 41184,
     3      40570, 40263, 39688, 39304, 38997, 39102, 39333, 39285,
     4      38757, 38690, 38344, 38009, 37500, 37020, 36752, 36550,
     5      36521, 36819, 37769, 38037, 38239, 39198, 39899, 39899,
     6      40139, 40062, 40340, 40388, 40849, 42115, 43113, 43228,
     7      44370, 44303, 44293, 44735, 44802, 45032, 45272, 45464,
     8      46375, 46414, 46289, 46404, 46183, 46135, 46020, 46164,
     9      46231, 46279, 46529, 46711, 53850, 55299, 55145, 52583/
            data (iy(i),i=  81, 160)/
     +      52420, 42950, 43363, 43737, 43670, 43785, 44341, 44562,
     1      44648, 44994, 45608, 46011, 46116, 46068, 46020, 46183,
     2      46471, 47009, 47488, 47690, 47920, 48419, 48841, 48995,
     3      49225, 50233, 50693, 51134, 51710, 52257, 52487, 52698,
     4      52718, 53312, 53428, 53274, 53668, 47690, 52382, 52564,
     5      53169, 53111, 52622, 52276, 51662, 51269, 50770, 50223,
     6      49945, 49369, 48774, 48640, 48314, 48285, 48006, 47959,
     7      45368, 46039, 46279, 46491, 47843, 48525, 49139, 49158,
     8      49340, 49148, 49129, 48937, 48381, 47719, 55798, 55692,
     9      53879, 53907, 52718, 52439, 52228, 50683, 50194, 49964/
            data (iy(i),i= 161, 240)/
     +      49235, 48544, 48333, 48218, 47469, 47239, 46845, 46682,
     1      46922, 47268, 47613, 48237, 48640, 49494, 49734, 50290,
     2      50472, 50664, 51298, 51614, 52420, 52651, 53111, 53572,
     3      53792, 54023, 54310, 54531, 54646, 55107, 55567, 55750,
     4      55807, 45205, 42058, 41942, 40935, 40388, 39189, 39208,
     5      40388, 40858, 40935, 40781, 40849, 40686, 40522, 40206,
     6      40043, 39832, 39467, 39237, 39169, 39074, 38661, 37529,
     7      36560, 36147, 35667, 35188, 45099, 44898, 44639, 44293,
     8      44063, 43353, 42825, 42720, 42489, 42403, 42153, 42086,
     9      41885, 27646, 27406, 27550, 28135, 28202, 28452, 28509/
            data (iy(i),i= 241, 320)/
     +      28740, 29383, 30323, 31148, 32386, 33336, 33432, 33729,
     1      34276, 33854, 33499, 33422, 33432, 33681, 33892, 34420,
     2      35552, 35600, 36061, 36243, 36358, 37481, 38296, 38623,
     3      38862, 39419, 39419, 39784, 40043, 40302, 40446, 40743,
     4      40993, 41309, 41472, 41309, 41645, 41664, 41827, 42173,
     5      42499, 42806, 32060, 32290, 32328, 32510, 32779, 33096,
     6      33556, 33863, 34698, 34890, 34689, 34756, 35159, 29200,
     7      31829, 31733, 31618, 31407, 30812, 29488, 26869, 23846,
     8      24249, 24249, 24422, 24624, 24988, 25036, 25641, 26389,
     9      26869, 26926, 28222, 27598, 26811, 29056, 28490, 23654/
            data (iy(i),i= 321, 400)/
     +      23367, 23232, 23194, 23031, 22877, 22647, 22705, 23338,
     1      24134, 24220, 24134, 24057, 23885, 23865, 24000, 23827,
     2      28039, 28193, 28116, 28260, 28471, 28682, 27329, 26312,
     3      26216, 33096, 32280, 31349, 30918, 30342, 30045, 29411,
     4      29162, 28864, 28768, 28452, 28346, 28481, 28557, 28759,
     5      28826, 27540, 33489, 33979, 34007, 34449, 34717, 34650,
     6      26149, 25583, 25218, 24873, 24691, 24067, 23885, 22340,
     7      21946, 21783, 21275, 20910, 21754, 22023, 23011, 23597,
     8      34660, 34756, 34804, 34890, 34986, 34852, 34487, 34286,
     9      34075, 33566, 32664, 31925, 31532, 31349, 31340, 31225/
            data (iy(i),i= 401, 480)/
     +      31071, 31167, 31119, 30956, 30745, 30553, 30496, 30659,
     1      30544, 30169, 28740, 28298, 28308, 28730, 29056, 29498,
     2      29613, 29565, 29402, 29104, 28816, 28615, 27809, 27953,
     3      27799, 27521, 27492, 53965, 54301, 54531, 55040, 52094,
     4      52027, 51058, 50069, 48851, 48669, 47469, 46615, 46519,
     5      46721, 46491, 45867, 45301, 45099, 44408, 44101, 43516,
     6      45627, 43353, 43526, 42739, 42873, 43075, 42998, 43391,
     7      43334, 43679, 43564, 42988, 42864, 42748, 42336, 42221,
     8      41875, 41578, 41472, 42441, 42499, 42662, 43334, 43679,
     9      43929, 44072, 44245, 42595, 42086, 41530, 41002, 40590/
            data (iy(i),i= 481, 560)/
     +      40225, 40033, 39678, 40273, 40426, 40321, 39880, 39745,
     1      39755, 39726, 39294, 39323, 39812, 40043, 40177, 40359,
     2      40570, 40734, 41002, 41511, 41952, 42307, 42403, 42720,
     3      43286, 43219, 43497, 43842, 43996, 43996, 44207, 44389,
     4      44687, 45109, 45320, 45090, 44859, 44677, 44581, 44197,
     5      43737, 43602, 43324, 43209, 42777, 42633, 42288, 42259,
     6      42067, 42163, 42048, 41453, 41060, 40781, 40302, 40369,
     7      40330, 39784, 39736, 39438, 39265, 39093, 38536, 38143,
     8      37817, 37769, 37605, 37174, 37174, 36828, 36848, 36742,
     9      36838, 37740, 38057, 38584, 35399, 35284, 35303, 35188/
            data (iy(i),i= 561, 640)/
     +      35427, 35773, 35821, 35686, 36061, 36368, 36502, 34967,
     1      34823, 34497, 34180, 33873, 33710, 33393, 33029, 32962,
     2      33010, 32846, 32367, 32367, 32204, 31992, 31714, 31561,
     3      31484, 31858, 32021, 31906, 31647, 31743, 32261, 32328,
     4      32213, 32069, 31868, 31685, 31637, 31311, 31110, 31292,
     5      31119, 31043, 30697, 30285, 30045, 29661, 29430, 29363,
     6      29200, 28701, 28653, 28423, 28241, 28049, 28097, 28135,
     7      27972, 28135, 28270, 28174, 27473, 27358, 27032, 26667,
     8      26485, 26235, 26610, 26523, 26264, 26043, 25784, 25708,
     9      25717, 26063, 26197, 26341, 26379, 26773, 26917, 27109/
            data (iy(i),i= 641, 720)/
     +      27301, 33019, 32971, 33115, 33230, 33566, 33796, 33470,
     1      33624, 33460, 33844, 33547, 33633, 34017, 34017, 34199,
     2      34238, 34573, 34919, 35399, 35629, 35648, 35370, 35370,
     3      36166, 36531, 36761, 37059, 37088, 37605, 38277, 38651,
     4      39553, 39160, 39074, 38392, 37605, 37193, 36684, 36550,
     5      35255, 35120, 35351, 35610, 35926, 35792, 35168, 35053,
     6      35591, 35504, 34737, 34669, 35130, 35360, 35053, 34890,
     7      34094, 33892, 33662, 33566, 33058, 33010, 32770, 32770,
     8      33134, 33163, 33412, 33624, 33806, 34314, 34497, 34516,
     9      34660, 34890, 35025, 35130, 35360, 34986, 34804, 34487/
            data (iy(i),i= 721, 800)/
     +      34535, 34238, 35370, 35293, 35571, 35667, 35715, 35667,
     1      35706, 35907, 36176, 36435, 36598, 37116, 37538, 38037,
     2      38613, 39093, 38910, 39045, 39553, 40052, 40781, 36483,
     3      36281, 36061, 35773, 35629, 34698, 34410, 34113, 33652,
     4      33307, 32981, 32962, 32194, 32156, 32530, 32722, 32501,
     5      29498, 29498, 29853, 29709, 29642, 29814, 29555, 28749,
     6      28797, 28980, 28912, 28874, 29200, 29363, 29661, 29757,
     7      29776, 30169, 30486, 30582, 30601, 30678, 30812, 30995,
     8      30995, 31158, 31302, 31513, 31762, 32252, 32885, 33067,
     9      33249, 33547, 33710, 34103, 33835, 33931, 34017, 33163/
            data (iy(i),i= 801, 880)/
     +      33000, 32952, 32818, 32539, 32482, 32050, 31839, 31033,
     1      30735, 30572, 30688, 30591, 30774, 30908, 30793, 30860,
     2      31225, 31359, 31177, 31129, 31196, 31359, 30640, 30995,
     3      31139, 31206, 31417, 31685, 31887, 31973, 32232, 32597,
     4      32942, 33000, 31369, 31004, 30659, 30217, 30179, 30083,
     5      30054, 30189, 30438, 30438, 30486, 30764, 30831, 30735,
     6      30419, 27646, 27339, 27358, 27080, 27147, 27608, 27579,
     7      27281, 27032, 26965, 27080, 26840, 26993, 27013, 27473,
     8      27655, 27867, 27934, 28250, 28394, 28577, 28797, 28855,
     9      29728, 29670, 29277, 29095, 29277, 29306, 29536, 29632/
            data (iy(i),i= 881, 960)/
     +      29766, 29958, 30093, 29795, 29997, 29392, 29210, 29056,
     1      29056, 29613, 29862, 28318, 28010, 27895, 28039, 28010,
     2      28922, 29008, 28749, 28385, 28145, 28145, 27770, 27790,
     3      28548, 28596, 29527, 29728, 30381, 30515, 27320, 27511,
     4      27857, 28001, 27790, 27943, 27780, 27828, 27761, 28366,
     5      28596, 28548, 28481, 26072, 25641, 25458, 25247, 24576,
     6      24441, 24192, 24077, 28682, 24009, 24393, 24883, 25190,
     7      25650, 28471, 28462, 28462, 28116, 27790, 27473, 27502,
     8      27732, 27665, 27579, 26581, 26581, 26878, 26993, 27157,
     9      27291, 27291, 27080, 26888, 26264, 26763, 27176, 27358/
            data (iy(i),i= 961, 968)/
     +      26715, 29268, 31014, 32309, 32357, 32520, 33710, 33000/
c
            data (b1(i),i=    1,  100)/
     +      1,    2,    3,    4,    5,    6,    7,    8,    9,   10,
     1     11,   12,   13,   14,   15,   16,   17,   18,   19,   20,
     2     21,   22,   23,   24,   25,   26,   27,   28,   29,   30,
     3     31,   32,   33,   34,   35,   36,   37,   38,   39,   40,
     4     41,   42,   43,   44,   45,   46,   47,   48,   49,   50,
     5     51,   52,   53,   54,   55,   56,   57,   58,   59,   60,
     6     61,   62,   63,   64,   65,   66,   67,   68,   69,   70,
     7     71,   72,   73,   74,   75,   78,   79,   80,   82,   83,
     8     84,   85,   86,   87,   88,   89,   90,   91,   92,   93,
     9     94,   95,   96,   97,   98,   99,  100,  101,  102,  103/
            data (b1(i),i=  101,  200)/
     +    104,  105,  106,  107,  108,  109,  110,  111,  112,  113,
     1    114,  115,  116,   77,   81,  119,  120,  121,  122,  123,
     2    124,  125,  126,  127,  128,  129,  130,  131,  132,  133,
     3    134,  135,  118,  137,  138,  139,  140,  141,  142,  143,
     4    144,  145,  146,  147,  148,  149,  118,  151,   78,   77,
     5    153,  154,  155,  156,  157,  158,  159,  160,  161,  162,
     6    163,  164,  165,  166,  167,    1,   76,  169,  170,  171,
     7    172,  173,  174,  175,  176,  177,  178,  179,  180,  181,
     8    182,  183,  184,  185,  186,  187,  188,  189,  190,  191,
     9    192,  151,  137,  195,  196,  197,  198,  199,  200,  201/
            data (b1(i),i=  201,  300)/
     +    202,  203,  204,  205,  206,  207,  208,  209,  210,  211,
     1    212,  213,  214,  215,  216,  217,  218,  219,  194,  221,
     2    222,  223,  224,  225,  226,  227,  228,  229,  230,  231,
     3    232,  195,  234,  235,  236,  237,  238,  239,  240,  241,
     4    242,  243,  244,  245,  246,  247,  248,  249,  250,  251,
     5    252,  253,  254,  255,  220,  257,  258,  259,  260,  261,
     6    262,  263,  264,  265,  266,  267,  268,  269,  270,  271,
     7    272,  273,  274,  275,  276,  277,  278,  279,  280,  281,
     8     82,  283,  284,  285,  286,  287,  288,  289,  290,  291,
     9    292,  293,  294,  257,  283,  297,  298,  299,  300,  301/
            data (b1(i),i=  301,  400)/
     +    296,  304,  305,  306,  307,  308,  309,  310,  311,  312,
     1    313,  303,  315,  316,  303,  296,  318,  315,  304,  320,
     2    321,  322,  323,  324,  325,  326,  327,  328,  329,  330,
     3    331,  332,  333,  334,  335,  234,  337,  338,  339,  340,
     4    341,  343,  344,  346,  347,  348,  349,  350,  351,  352,
     5    353,  354,  355,  356,  357,  358,  359,  360,  342,  343,
     6    346,  363,  364,  365,  366,  367,  345,  369,  370,  371,
     7    372,  373,  374,  375,  376,  377,  378,  379,  380,  381,
     8    382,  383,  336,  368,  385,  386,  387,  388,  389,  390,
     9    391,  392,  393,  394,  395,  396,  397,  398,  399,  400/
            data (b1(i),i=  401,  500)/
     +    401,  402,  403,  404,  405,  406,  407,  408,  409,  410,
     1    411,  412,  413,  414,  415,  416,  417,  418,  419,  420,
     2    421,  422,  423,  424,  425,  426,  362,    1,   77,  428,
     3    429,  430,  431,   81,  432,  433,  434,  435,  436,  437,
     4    438,  439,  440,  441,  442,  443,  444,  445,  446,  447,
     5    448,  118,  449,  450,  451,  452,  453,  454,  455,  456,
     6    457,  458,  459,  460,  461,  462,  463,  464,  465,  466,
     7    467,  468,  469,  470,  471,  472,  473,  474,  447,  450,
     8    476,  477,  478,  479,  480,  481,  482,  483,  484,  485,
     9    486,  487,  488,  489,  490,  491,  492,  493,  494,  495/
            data (b1(i),i=  501,  600)/
     +    496,  497,  498,  499,  500,  501,  502,  503,  504,  505,
     1    506,  507,  508,  509,  510,  511,  512,  513,  514,  515,
     2    137,  194,  517,  518,  519,  520,  521,  522,  523,  524,
     3    525,  526,  527,  528,  529,  530,  531,  532,  533,  534,
     4    535,  536,  537,  538,  539,  540,  541,  542,  543,  544,
     5    545,  546,  547,  548,  549,  550,  551,  552,  553,  554,
     6    555,  492,  220,  557,  558,  559,  560,  561,  562,  563,
     7    564,  565,  566,  550,  559,  568,  569,  570,  571,  572,
     8    573,  574,  575,  576,  577,  578,  579,  580,  581,  582,
     9    583,  584,  585,  586,  587,  588,  589,  590,  591,  592/
            data (b1(i),i=  601,  700)/
     +    593,  594,  595,  596,  597,  598,  599,  600,  601,  602,
     1    603,  604,  605,  606,  607,  608,  609,  610,  611,  612,
     2    613,  614,  615,  616,  617,  618,  619,  620,  621,  622,
     3    623,  624,  625,  626,  627,  628,  629,  630,  631,  632,
     4    633,  634,  635,  636,  637,  638,  639,  640,  234,  574,
     5    642,  643,  644,  645,  646,  647,  648,  649,  650,  651,
     6    652,  653,  654,  655,  656,  657,  658,  659,  660,  661,
     7    662,  663,  664,  665,  666,  667,  668,  669,  670,  671,
     8    492,  484,  673,  674,  675,  676,  677,  678,  679,  680,
     9    681,  682,  683,  684,  685,  686,  687,  688,  689,  690/
            data (b1(i),i=  701,  800)/
     +    691,  692,  693,  694,  695,  696,  697,  698,  699,  700,
     1    701,  702,  703,  704,  705,  706,  707,  708,  709,  710,
     2    711,  712,  713,  714,  715,  716,  717,  718,  719,  720,
     3    721,  657,  681,  723,  724,  725,  726,  727,  728,  729,
     4    730,  731,  732,  733,  734,  735,  736,  737,  738,  739,
     5    740,  741,  742,  467,  733,  744,  745,  746,  747,  257,
     6    682,  749,  750,  751,  752,  753,  754,  755,  756,  757,
     7    758,  759,  760,  296,  761,  762,  763,  764,  765,  766,
     8    767,  768,  769,  770,  771,  772,  773,  774,  775,  776,
     9    777,  778,  779,  780,  781,  782,  783,  784,  785,  786/
            data (b1(i),i=  801,  900)/
     +    787,  788,  789,  790,  791,  792,  793,  794,  795,  657,
     1    720,  797,  798,  712,  650,  800,  801,  802,  803,  804,
     2    805,  806,  807,  808,  809,  810,  811,  812,  813,  814,
     3    815,  816,  817,  818,  819,  820,  821,  788,  823,  824,
     4    825,  826,  827,  828,  829,  830,  831,  832,  833,  645,
     5    808,  589,  835,  836,  837,  838,  839,  840,  841,  842,
     6    843,  844,  845,  846,  847,  848,  826,  617,  850,  851,
     7    852,  853,  854,  855,  856,  857,  858,  859,  860,  861,
     8    862,  863,  864,  865,  866,  867,  868,  869,  870,  871,
     9    872,  873,  874,  875,  876,  877,  878,  879,  880,  881/
            data (b1(i),i=  901, 1000)/
     +    882,  838,  873,  841,  885,  886,  887,  888,  889,  890,
     1    885,  892,  893,  894,  895,  865,  871,  889,  887,  898,
     2    892,  892,  900,  901,  902,  903,  904,  905,  906,  907,
     3    908,  909,  823,  303,  911,  912,  913,  914,  915,  916,
     4    917,  918,  919,  920,  921,  922,  906,  861,  924,  925,
     5    926,  927,  928,  929,  930,  931,  768,  932,  336,  933,
     6    934,  935,  936,  633,  342,  938,  939,  940,  941,  942,
     7    943,  944,  945,  946,  947,  345,  948,  949,  950,  951,
     8    952,  953,  954,  955,  956,  957,  958,  959,  960,  635,
     9    362,  962,  963,  964,  965,  346,  368,  967,  963,    0/
            data (b2(i),i=    1,  100)/
     +      2,    3,    4,    5,    6,    7,    8,    9,   10,   11,
     1     12,   13,   14,   15,   16,   17,   18,   19,   20,   21,
     2     22,   23,   24,   25,   26,   27,   28,   29,   30,   31,
     3     32,   33,   34,   35,   36,   37,   38,   39,   40,   41,
     4     42,   43,   44,   45,   46,   47,   48,   49,   50,   51,
     5     52,   53,   54,   55,   56,   57,   58,   59,   60,   61,
     6     62,   63,   64,   65,   66,   67,   68,   69,   70,   71,
     7     72,   73,   74,   75,   76,   79,   80,   81,   83,   84,
     8     85,   86,   87,   88,   89,   90,   91,   92,   93,   94,
     9     95,   96,   97,   98,   99,  100,  101,  102,  103,  104/
            data (b2(i),i=  101,  200)/
     +    105,  106,  107,  108,  109,  110,  111,  112,  113,  114,
     1    115,  116,  117,  117,  119,  120,  121,  122,  123,  124,
     2    125,  126,  127,  128,  129,  130,  131,  132,  133,  134,
     3    135,  136,  136,  138,  139,  140,  141,  142,  143,  144,
     4    145,  146,  147,  148,  149,  150,  150,  152,  152,  153,
     5    154,  155,  156,  157,  158,  159,  160,  161,  162,  163,
     6    164,  165,  166,  167,  168,  168,  169,  170,  171,  172,
     7    173,  174,  175,  176,  177,  178,  179,  180,  181,  182,
     8    183,  184,  185,  186,  187,  188,  189,  190,  191,  192,
     9    193,  193,  194,  196,  197,  198,  199,  200,  201,  202/
            data (b2(i),i=  201,  300)/
     +    203,  204,  205,  206,  207,  208,  209,  210,  211,  212,
     1    213,  214,  215,  216,  217,  218,  219,  220,  221,  222,
     2    223,  224,  225,  226,  227,  228,  229,  230,  231,  232,
     3    233,  233,  235,  236,  237,  238,  239,  240,  241,  242,
     4    243,  244,  245,  246,  247,  248,  249,  250,  251,  252,
     5    253,  254,  255,  256,  256,  258,  259,  260,  261,  262,
     6    263,  264,  265,  266,  267,  268,  269,  270,  271,  272,
     7    273,  274,  275,  276,  277,  278,  279,  280,  281,  282,
     8    282,  284,  285,  286,  287,  288,  289,  290,  291,  292,
     9    293,  294,  295,  295,  297,  298,  299,  300,  301,  302/
            data (b2(i),i=  301,  400)/
     +    302,  305,  306,  307,  308,  309,  310,  311,  312,  313,
     1    314,  314,  316,  317,  317,  318,  319,  319,  320,  321,
     2    322,  323,  324,  325,  326,  327,  328,  329,  330,  331,
     3    332,  333,  334,  335,  336,  337,  338,  339,  340,  341,
     4    342,  344,  345,  347,  348,  349,  350,  351,  352,  353,
     5    354,  355,  356,  357,  358,  359,  360,  361,  361,  362,
     6    363,  364,  365,  366,  367,  368,  369,  370,  371,  372,
     7    373,  374,  375,  376,  377,  378,  379,  380,  381,  382,
     8    383,  384,  384,  385,  386,  387,  388,  389,  390,  391,
     9    392,  393,  394,  395,  396,  397,  398,  399,  400,  401/
            data (b2(i),i=  401,  500)/
     +    402,  403,  404,  405,  406,  407,  408,  409,  410,  411,
     1    412,  413,  414,  415,  416,  417,  418,  419,  420,  421,
     2    422,  423,  424,  425,  426,  427,  427,   76,  428,  429,
     3    430,  431,   78,  432,  433,  434,  435,  436,  437,  438,
     4    439,  440,  441,  442,  443,  444,  445,  446,  447,  448,
     5     82,  449,  450,  451,  452,  453,  454,  455,  456,  457,
     6    458,  459,  460,  461,  462,  463,  464,  465,  466,  467,
     7    468,  469,  470,  471,  472,  473,  474,  475,  475,  476,
     8    477,  478,  479,  480,  481,  482,  483,  484,  485,  486,
     9    487,  488,  489,  490,  491,  492,  493,  494,  495,  496/
            data (b2(i),i=  501,  600)/
     +    497,  498,  499,  500,  501,  502,  503,  504,  505,  506,
     1    507,  508,  509,  510,  511,  512,  513,  514,  515,  516,
     2    516,  517,  518,  519,  520,  521,  522,  523,  524,  525,
     3    526,  527,  528,  529,  530,  531,  532,  533,  534,  535,
     4    536,  537,  538,  539,  540,  541,  542,  543,  544,  545,
     5    546,  547,  548,  549,  550,  551,  552,  553,  554,  555,
     6    556,  556,  557,  558,  559,  560,  561,  562,  563,  564,
     7    565,  566,  567,  567,  568,  569,  570,  571,  572,  573,
     8    574,  575,  576,  577,  578,  579,  580,  581,  582,  583,
     9    584,  585,  586,  587,  588,  589,  590,  591,  592,  593/
            data (b2(i),i=  601,  700)/
     +    594,  595,  596,  597,  598,  599,  600,  601,  602,  603,
     1    604,  605,  606,  607,  608,  609,  610,  611,  612,  613,
     2    614,  615,  616,  617,  618,  619,  620,  621,  622,  623,
     3    624,  625,  626,  627,  628,  629,  630,  631,  632,  633,
     4    634,  635,  636,  637,  638,  639,  640,  641,  641,  642,
     5    643,  644,  645,  646,  647,  648,  649,  650,  651,  652,
     6    653,  654,  655,  656,  657,  658,  659,  660,  661,  662,
     7    663,  664,  665,  666,  667,  668,  669,  670,  671,  672,
     8    672,  673,  674,  675,  676,  677,  678,  679,  680,  681,
     9    682,  683,  684,  685,  686,  687,  688,  689,  690,  691/
            data (b2(i),i=  701,  800)/
     +    692,  693,  694,  695,  696,  697,  698,  699,  700,  701,
     1    702,  703,  704,  705,  706,  707,  708,  709,  710,  711,
     2    712,  713,  714,  715,  716,  717,  718,  719,  720,  721,
     3    722,  722,  723,  724,  725,  726,  727,  728,  729,  730,
     4    731,  732,  733,  734,  735,  736,  737,  738,  739,  740,
     5    741,  742,  743,  743,  744,  745,  746,  747,  748,  748,
     6    749,  750,  751,  752,  753,  754,  755,  756,  757,  758,
     7    759,  760,  283,  761,  762,  763,  764,  765,  766,  767,
     8    768,  769,  770,  771,  772,  773,  774,  775,  776,  777,
     9    778,  779,  780,  781,  782,  783,  784,  785,  786,  787/
            data (b2(i),i=  801,  900)/
     +    788,  789,  790,  791,  792,  793,  794,  795,  796,  796,
     1    797,  798,  799,  799,  800,  801,  802,  803,  804,  805,
     2    806,  807,  808,  809,  810,  811,  812,  813,  814,  815,
     3    816,  817,  818,  819,  820,  821,  822,  822,  824,  825,
     4    826,  827,  828,  829,  830,  831,  832,  833,  834,  834,
     5    823,  835,  836,  837,  838,  839,  840,  841,  842,  843,
     6    844,  845,  846,  847,  848,  849,  849,  850,  851,  852,
     7    853,  854,  855,  856,  857,  858,  859,  860,  861,  862,
     8    863,  864,  865,  866,  867,  868,  869,  870,  871,  872,
     9    873,  874,  875,  876,  877,  878,  879,  880,  881,  882/
            data (b2(i),i=  901, 1000)/
     +    883,  883,  884,  884,  886,  887,  888,  889,  890,  891,
     1    891,  893,  894,  895,  896,  896,  897,  897,  898,  899,
     2    899,  900,  901,  902,  903,  904,  905,  906,  907,  908,
     3    909,  910,  910,  911,  912,  913,  914,  915,  916,  917,
     4    918,  919,  920,  921,  922,  923,  923,  924,  925,  926,
     5    927,  928,  929,  930,  931,  304,  932,  315,  933,  934,
     6    935,  936,  937,  937,  938,  939,  940,  941,  942,  943,
     7    944,  945,  946,  947,  343,  948,  949,  950,  951,  952,
     8    953,  954,  955,  956,  957,  958,  959,  960,  961,  961,
     9    962,  963,  964,  965,  966,  966,  967,  968,  968,    0/
c
        sp(2)='mexico'
        sp(1)='mexico'
        sp(3)='mexico'
c
c       istate=0/1  dont use/use the internal edges/vertices that
c       define state borders
c
        istate=1
c
        do i=1,nvf
            vx(i)=real(ix(i),rknd)*1.0e-2_rknd
            vy(i)=real(iy(i),rknd)*1.0e-2_rknd
        enddo
        do i=1,nbf
            ibndry(1,i)=b1(i)
            ibndry(2,i)=b2(i)
            ibndry(3,i)=0
            ibndry(4,i)=2
            if(i>nbound) ibndry(4,i)=0
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=0
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
        if(istate==0) then
            ip(2)=nbound
            ip(3)=nbound
        endif
c
c       make itnode
c
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd20(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), save :: ntf,nvf,nbf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), dimension(200) :: xm,ym
            real(kind=rknd), save, dimension(5) :: ratio,diam,xcen,ycen
            real(kind=rknd), save, dimension(5) :: angle
            real(kind=rknd), save :: hmax,grade
            character(len=80), dimension(100) :: sp,su
cy
            external sxy
            data hmax,grade/0.1e0_rknd,1.75e0_rknd/
            data ratio/4.0e0_rknd,1.0e0_rknd,
     +          0.2e0_rknd,2.0e0_rknd,0.7e0_rknd/
            data diam/4.0e0_rknd,3.0e0_rknd,
     +          4.0e0_rknd,2.0e0_rknd,4.0e0_rknd/
            data angle/0.0e0_rknd,0.0e0_rknd,
     +          45.0e0_rknd,30.0e0_rknd,60.0e0_rknd/
            data xcen/0.0e0_rknd,2.0e0_rknd,
     +          5.0e0_rknd,7.0e0_rknd,9.0e0_rknd/
            data ycen/2.0e0_rknd,-1.0e0_rknd,
     +          3.0e0_rknd,-2.0e0_rknd,1.5e0_rknd/
c
        sp(2)='ellipse'
        sp(1)='ellipse'
        sp(3)='ellipse'
c
        ip(1)=0
        ip(2)=0
        ip(3)=0
c
        nel=5
        nn=32
c
        do i=1,nel-1
            ii=i
            do j=i+1,nel
                if(xcen(j)<xcen(ii)) ii=j
                if(xcen(j)==xcen(ii).and.ycen(j)<ycen(ii)) ii=j
            enddo
            if(ii/=i) then
                t=xcen(i)
                xcen(i)=xcen(ii)
                xcen(ii)=t
                t=ycen(i)
                ycen(i)=ycen(ii)
                ycen(ii)=t
                t=diam(i)
                diam(i)=diam(ii)
                diam(ii)=t
                t=ratio(i)
                ratio(i)=ratio(ii)
                ratio(ii)=t
                t=angle(i)
                angle(i)=angle(ii)
                angle(ii)=t
            endif
        enddo
c
        xmin=xcen(1)
        xmax=xmin
        ymin=ycen(1)
        ymax=ymin
        rmax=0.0e0_rknd
        do i=1,nel
            r=diam(i)/2.0e0_rknd
            xmin=min(xcen(i)-r,xmin)
            xmax=max(xcen(i)+r,xmax)
            ymin=min(ycen(i)-r,ymin)
            ymax=max(ycen(i)+r,ymax)
            rmax=max(r,rmax)
        enddo
c
        nbf=0
        nvf=0
        ncf=0
c
        do i=1,nel
            call ellpse(vx(nvf+1),vy(nvf+1),xm(ncf+1),ym(ncf+1),
     +          ratio(i),angle(i),diam(i),xcen(i),ycen(i),nn)
            do k=1,nn
                ibndry(1,nbf+k)=nvf+k
                ibndry(2,nbf+k)=nvf+k+1
                if(k==nn) ibndry(2,nbf+k)=nvf+1
                ibndry(3,nbf+k)=ncf+k
                ibndry(4,nbf+k)=1
                ibndry(5,nbf+k)=0
                ibndry(6,nbf+k)=0
                ibndry(7,nbf+k)=-(nbf+k)
            enddo
            nvf=nvf+nn
            nbf=nbf+nn
            ncf=ncf+nn
        enddo
c
        vx(nvf+1)=xmin-rmax
        vy(nvf+1)=ymin-rmax
        vx(nvf+2)=xmax+rmax
        vy(nvf+2)=ymin-rmax
        vx(nvf+3)=xmax+rmax
        vy(nvf+3)=ymax+rmax
        vx(nvf+4)=xmin-rmax
        vy(nvf+4)=ymax+rmax
c
        do i=1,4
            ibndry(1,nbf+i)=nvf+i
            ibndry(2,nbf+i)=nvf+i+1
            ibndry(3,nbf+i)=0
            ibndry(4,nbf+i)=2
            ibndry(5,nbf+i)=0
            ibndry(6,nbf+i)=0
            ibndry(7,nbf+i)=nbf+i
        enddo
        ibndry(2,nbf+4)=nvf+1
        nvf=nvf+4
        nbf=nbf+4
c
        do i=1,nel+1
            do j=1,7
                ibndry(j,nbf+i)=0
            enddo
        enddo
c
        ii=nvf-3
        i=1
        k1=i
        k2=k1+nn-1
        d0=(vx(ii)-vx(i))**2+(vy(ii)-vy(i))**2
        do k=k1,k2
            dd=(vx(ii)-vx(k))**2+(vy(ii)-vy(k))**2
            if(dd<d0) then
                d0=dd
                i=k
            endif
        enddo
        ibndry(1,nbf+1)=ii
        ibndry(2,nbf+1)=i
        nbf=nbf+1
c
        do m=1,nel-1
            ii=1+(m-1)*nn
            i=1+m*nn
            d0=(vx(ii)-vx(i))**2+(vy(ii)-vy(i))**2
            do mm=1,2
                k1=1+m*nn
                k2=k1+nn-1
                do k=k1,k2
                    dd=(vx(ii)-vx(k))**2+(vy(ii)-vy(k))**2
                    if(dd<d0) then
                        d0=dd
                        i=k
                    endif
                enddo
                k1=1+(m-1)*nn
                k2=k1+nn-1
                do k=k1,k2
                    dd=(vx(i)-vx(k))**2+(vy(i)-vy(k))**2
                    if(dd<d0) then
                        d0=dd
                        ii=k
                    endif
                enddo
            enddo
            ibndry(1,nbf+1)=ii
            ibndry(2,nbf+1)=i
            nbf=nbf+1
        enddo
c
        ii=nvf-1
        i=1+(nel-1)*nn
        k1=i
        k2=k1+nn-1
        d0=(vx(ii)-vx(i))**2+(vy(ii)-vy(i))**2
        do k=k1,k2
            dd=(vx(ii)-vx(k))**2+(vy(ii)-vy(k))**2
            if(dd<d0) then
                d0=dd
                i=k
            endif
        enddo
        ibndry(1,nbf+1)=ii
        ibndry(2,nbf+1)=i
        nbf=nbf+1
c
        ip(1)=2
        ip(2)=nvf
        ip(3)=nbf
c
        do i=1,nbf
            if(ibndry(3,i)>0) then
                sf(1,i)=xm(ibndry(3,i))
                sf(2,i)=ym(ibndry(3,i))
            else
                sf(1,i)=0.0e0_rknd
                sf(2,i)=0.0e0_rknd
            endif
        enddo
c
c       make itnode, find symmetries
c
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine ellpse(vx,vy,xm,ym,ratio,angle,diam,xcen,ycen,nn)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            real(kind=rknd), dimension(*) :: vx,vy,xm,ym
            real(kind=rknd), dimension(1025) :: x,y
cy
        if(nn>1024) stop 5567
        pi=3.141592653589793e0_rknd
        if(ratio<1.0e0_rknd) then
           rat=1.0e0_rknd/ratio
           ang=(angle+90.0e0_rknd)/180.0e0_rknd
        else
           rat=ratio
           ang=angle/180.0e0_rknd
        endif
c
        c=cos(ang*pi)
        s=sin(ang*pi)
        do i=1,2*nn+1
            ang=real(i-1,rknd)*pi/real(nn,rknd)
            xx=cos(ang)/2.0e0_rknd
            yy=sin(ang)/(2.0e0_rknd*rat)
            x(i)=(c*xx+s*yy)*diam+xcen
            y(i)=(c*yy-s*xx)*diam+ycen
        enddo
        do i=1,nn
            ii=2*i
            vx(i)=x(ii-1)
            vy(i)=y(ii-1)
            call centre(x(ii-1),y(ii-1),x(ii+1),y(ii+1),
     +              x(ii),y(ii),xm(i),ym(i))
        enddo
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine fixit(n,x,y)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(20000) :: ix,iy
            real(kind=rknd), dimension(*) :: x,y
            character(len=1), dimension(10) :: mark
cy
            data mark/'+','1','2','3','4','5','6','7','8','9'/
        scale=1.0e1_rknd
        do i=1,n
            if(x(i)>0.0e0_rknd) then
                ix(i)=int(x(i)*scale+0.5e0_rknd)
            else if(x(i)<0) then
                ix(i)=int(x(i)*scale-0.5e0_rknd)
            else
                ix(i)=0
            endif
            if(y(i)>0) then
                iy(i)=int(y(i)*scale+0.5e0_rknd)
            else if(y(i)<0) then
                iy(i)=int(y(i)*scale-0.5e0_rknd)
            else
                iy(i)=0
            endif
        enddo
        do i=n+1,n+100
            ix(i)=0
            iy(i)=0
        enddo
  101   format(12x,'data (ix(i),i=',i4,',',i4,')/')
  104   format(12x,'data (iy(i),i=',i4,',',i4,')/')
  102   format(5x,a1,5x,12(i4,','))
  103   format(5x,a1,5x,11(i4,','),i4,'/')
        nlin=12
        do i=1,n,10*nlin
            write(8,101) i,min(i+10*nlin-1,n)
            do j=1,10
            istrt=i+(j-1)*nlin
            if(j/=10) then
                write(8,102) mark(j),(ix(k),k=istrt,istrt+nlin-1)
            else
                write(8,103) mark(j),(ix(k),k=istrt,istrt+nlin-1)
            endif
            enddo
        enddo
        do i=1,n,10*nlin
            write(8,104) i,min(i+10*nlin-1,n)
            do j=1,10
            istrt=i+(j-1)*nlin
            if(j/=10) then
                write(8,102) mark(j),(iy(k),k=istrt,istrt+nlin-1)
            else
                write(8,103) mark(j),(iy(k),k=istrt,istrt+nlin-1)
            endif
            enddo
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd21(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), save :: ntf,nbf,nvf
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), save :: hmax
            character(len=80), dimension(100) :: sp,su
cy
            external sxy
            data ntf,nbf,nvf/2,12,10/
            data hmax/0.5e0_rknd/
c
        sp(2)='hole'
        sp(1)='hole'
        sp(3)='hole'
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
c
        vx(1)=-2.0e0_rknd
        vy(1)=2.0e0_rknd
        vx(2)=-2.0e0_rknd
        vy(2)=-2.0e0_rknd
        vx(3)=0.0e0_rknd
        vy(3)=-2.0e0_rknd
        vx(4)=2.0e0_rknd
        vy(4)=-2.0e0_rknd
        vx(5)=2.0e0_rknd
        vy(5)=2.0e0_rknd
        vx(6)=0.0e0_rknd
        vy(6)=2.0e0_rknd
        vx(7)=0.0e0_rknd
        vy(7)=1.0e0_rknd
        vx(8)=-1.0e0_rknd
        vy(8)=0.0e0_rknd
        vx(9)=0.0e0_rknd
        vy(9)=-1.0e0_rknd
        vx(10)=1.0e0_rknd
        vy(10)=0.0e0_rknd
c
        do i=1,nbf
            ibndry(1,i)=i-1
            ibndry(2,i)=i
            ibndry(3,i)=0
            if(i>=8.and.i<=11) ibndry(3,i)=1
            ibndry(4,i)=1
            if(i>=2.and.i<=5) ibndry(4,i)=2
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=0
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
c
        ibndry(1,1)=6
        ibndry(1,12)=3
        ibndry(2,11)=7
        ibndry(2,12)=9
        ibndry(4,7)=0
        ibndry(4,12)=0
c
        itnode(1,1)=6
        itnode(2,1)=1
        itnode(3,1)=0
        itnode(4,1)=0
        itnode(5,1)=1
        itnode(1,2)=6
        itnode(2,2)=7
        itnode(3,2)=-1
        itnode(4,2)=0
        itnode(5,2)=2
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd22(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), save, dimension(2000) :: ix,iy
            integer(kind=iknd), save :: ispd=1,iprob=1,itask=0,ifirst=1
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            character(len=80), dimension(100) :: sp,su
cy
            external sxy
c
            data (ix(i),i=    1,  100)/
     +  55110,55856,56437,56519,56418,56905,57294,57852,58624,59675,
     1  60375,60644,60972,61319,61603,61757,61792,61986,62091,62074,
     2  61965,61997,62148,62286,62458,62587,62698,62806,62882,63203,
     3  63123,63000,62778,62628,62325,62092,61833,61745,61287,61222,
     4  61075,60832,60598,59610,58479,58285,57914,57466,56332,55617,
     5  55536,55298,55281,55771,55955,56214,56792,56935,56907,56797,
     6  56656,56612,56612,56716,56721,56511,56335,56231,56030,55670,
     7  55281,55152,55135,54977,54864,54762,54549,54379,54200,53893,
     8  53555,53363,53165,52514,52369,52385,52514,52981,53136,53364,
     9  53380,53293,53135,53006,52931,52951,52745,52765,52680,52534/
            data (ix(i),i=  101,  200)/
     +  52067,51992,51984,52373,52339,52197,52101,51893,51795,51817,
     1  51769,51405,51277,51025,50881,50812,50683,50556,50313,50237,
     2  50133,50074,49982,49842,49659,49664,49801,49934,50284,50811,
     3  51392,51589,51580,51413,50687,50430,50097,50009,49935,49976,
     4  50279,50195,50079,49975,49839,49709,49594,49364,49370,49381,
     5  49421,49439,49440,49381,49324,49272,49130,49047,48743,48589,
     6  48504,48628,48595,48499,48386,48311,48315,48355,48355,48283,
     7  48225,48173,48140,48145,48116,48060,47968,47864,47514,47206,
     8  47080,47095,47214,47191,47405,47420,47304,47120,46977,46887,
     9  46731,46678,46742,46688,46596,46387,46232,46090,45975,46029/
            data (ix(i),i=  201,  300)/
     +  45965,45638,45534,45417,44674,44336,44350,44532,44703,44742,
     1  44573,44430,44273,44116,44038,43972,44129,44286,44627,45309,
     2  45531,45739,45909,45949,45690,45691,45743,45985,46094,46111,
     3  46076,45991,45892,45822,45730,45581,45417,45217,45050,44841,
     4  44697,44642,44679,44795,44946,45078,45209,45327,45373,45336,
     5  45266,45192,44903,44483,44169,44037,43906,43945,43866,43696,
     6  43618,43434,42880,42735,42562,42389,42214,42278,42410,43004,
     7  43689,44033,46013,46617,47061,47858,49646,49999,50637,51485,
     8  52141,52546,54087,54289,54796,28114,28346,28587,28903,29184,
     9  29356,29315,29149,28986,28656,28155,28046,28228,28345,28564/
            data (ix(i),i=  301,  400)/
     +  28845,28923,28744,28928,29349,29759,29958,30035,30720,30891,
     1  30997,31028,31241,31325,31454,31441,31584,31726,31831,31892,
     2  31882,31826,31730,31642,31541,31471,31372,31333,31320,31269,
     3  31186,31102,30913,30836,30915,30881,30816,30733,29899,29697,
     4  29560,29624,29257,28980,28688,28173,27708,27630,27599,27620,
     5  27671,27822,27978,28070,28009,27905,27758,27629,27520,27441,
     6  27290,27208,27065,26816,26620,26614,26683,26818,26950,27105,
     7  27265,27276,27157,26933,26727,26601,26516,26377,26172,26054,
     8  25991,25901,25763,25722,25883,26269,26510,26902,27777,27950,
     9  27986,20188,24195,28203,32210,36217,40225,44232,48240,52247/
            data (ix(i),i=  401,  500)/
     +  56254,60262,60937,61075,64269,64347,64913,65248,65411,65920,
     1  65962,66228,66861,67349,68276,72284,76291,76291,76291,76291,
     2  76291,76291,76291,76291,76291,76291,76291,76291,76291,76291,
     3  76291,76291,76291,76291,76291,76291,76291,76291,76291,76291,
     4  76387,76407,76460,76566,76715,76728,76717,76653,76678,76892,
     5  76926,76983,77042,77073,77153,77334,77339,77394,77415,77411,
     6  77278,77199,77090,76884,76721,76660,76617,76686,76803,76827,
     7  76704,76701,76827,76899,77022,77130,77175,77165,77282,77408,
     8  77474,77591,77633,77708,77834,77868,77847,77863,77906,78037,
     9  78247,78331,78526,78670,78729,78837,78881,78859,78886,78981/
            data (ix(i),i=  501,  600)/
     +  79217,79281,79350,79361,79380,79564,79625,79744,79787,79967,
     1  80085,80076,80132,80226,80237,80206,80282,80283,80288,80370,
     2  80513,80607,80761,80950,81018,81149,81153,81195,81282,81331,
     3  81349,81247,81232,81151,81177,81238,81331,81396,81459,81544,
     4  81627,81661,81699,81681,81695,81633,81619,81702,81792,81827,
     5  81883,81862,81894,81985,82032,82083,82190,82293,82328,82478,
     6  82452,82518,82665,82722,82715,82641,82616,82634,82804,82943,
     7  82951,82935,83075,83099,83104,82901,82873,82882,82956,82943,
     8  82985,83138,83177,83189,83118,83122,83231,83296,83347,83495,
     9  83496,83372,83309,83310,83346,83398,83562,83605,83705,83761/
            data (ix(i),i=  601,  700)/
     +  83890,83918,84039,84070,84089,84074,84170,84236,84295,84373,
     1  84385,84538,84695,84733,84893,85149,85188,85227,85327,85336,
     2  85287,85247,85083,84990,84931,84909,84953,85054,85070,85415,
     3  85492,85584,85710,85833,85877,85915,86089,86387,86603,86683,
     4  86822,86832,86820,86868,86890,87018,87006,86998,87020,87217,
     5  87280,87424,87444,87400,87427,87600,87636,87774,87801,87948,
     6  88155,88162,88294,88330,88413,88431,88455,88444,88581,88606,
     7  88751,88802,88938,89019,89100,89069,89072,89282,89350,89445,
     8  89583,89631,89674,89755,89735,89738,89820,89843,89978,89990,
     9  89979,90016,90169,90342,90453,90483,90530,90549,90693,90744/
            data (ix(i),i=  701,  800)/
     +  90768,90756,90791,90844,90916,91089,91135,91176,91179,91142,
     1  91128,91209,91189,91251,91301,91425,91481,91637,91714,91700,
     2  91730,91827,91856,91982,92053,92171,92220,92314,92314,92226,
     3  92224,92313,92328,92401,92686,92722,92755,93091,93270,93281,
     4  93414,93556,93689,93717,93805,93874,93953,94015,94044,93944,
     5  93637,93623,93652,93780,93984,94038,94093,94093,94227,94355,
     6  94382,94494,94543,94667,94681,94784,94832,94939,95050,95101,
     7  95135,95125,95191,95182,95079,95110,95231,95394,95368,95393,
     8  95481,95528,95607,95615,95722,95880,96021,96134,96243,96320,
     9  96448,96501,96661,96806,96832,96864,96885,96936,97048,97079/
            data (ix(i),i=  801,  900)/
     +  97210,97250,97243,97143,97133,97194,97245,97271,97244,97389,
     1  97374,97365,97390,97492,97641,97620,97672,97511,97534,97531,
     2  97605,97686,97703,97697,97664,97717,97640,97594,97601,97599,
     3  97411,97315,97307,97393,97528,97597,97735,97927,97958,97979,
     4  97882,97855,97863,97871,97934,97994,98114,98126,98158,98242,
     5  98281,98310,98341,98458,98468,98595,98689,98694,98656,98673,
     6  98731,98936,99098,99148,99298,99492,99582,99609,99605,99750,
     7  99894,99974,100000,99995,99950,99890,98675,97417,96319,94463,
     8  92086,88190,86574,83508,80930,79027,78616,77777,77236,76681,
     9  74523,72862,71897,71823,70209,68481,67254,67180,66895,65288/
            data (ix(i),i=  901, 1000)/
     +  65208,64904,64831,64628,64438,64024,63994,64044,64008,63915,
     1  63791,63602,63702,63305,63215,63256,63361,63536,64052,64499,
     2  64670,64816,64709,64603,64493,64274,63684,63187,63075,63024,
     3  63084,63000,63179,63224,63215,63247,63310,63460,63500,63489,
     4  63415,63361,63246,63212,63167,63094,62953,62748,62566,62288,
     5  62186,62124,61906,61581,61236,60470,60372,60150,59919,59839,
     6  59879,60025,60131,60244,60668,60872,60984,61033,61263,61579,
     7  61636,61538,61082,60841,60473,60563,60577,60325,60344,60430,
     8  60817,60995,61017,60902,60482,60367,60246,60292,60731,60812,
     9  60728,60167,60025,60029,60274,60305,60256,60101,59952,59876/
            data (ix(i),i= 1001, 1100)/
     +  59808,59840,59730,59611,59301,58669,58395,57840,57228,56740,
     1  56643,56648,56792,57057,57225,57254,57083,57039,57167,57482,
     2  57561,57585,57303,57376,57474,57767,58027,58321,58635,58683,
     3  58513,58119,57879,57553,57248,56944,56756,56471,56277,56149,
     4  56002,55929,55835,55724,55828,56383,56685,56652,56521,56597,
     5  56751,56758,56585,56475,56231,56374,56250,55646,55424,55265,
     6  55083,54986,54727,54632,54554,54394,54345,54249,54116,53961,
     7  53888,53922,54008,54268,54209,54131,53983,53774,53603,53387,
     8  53228,53022,52948,53033,52930,52188,51783,51519,51284,51167,
     9  51105,51386,52050,52092,51883,51327,51183,51134,51928,52757/
            data (ix(i),i= 1101, 1200)/
     +  53277,53629,53694,53644,53701,53944,54030,53967,53780,53701,
     1  53678,53811,53725,53411,53354,53382,53124,52804,52646,51406,
     2  51355,51533,51933,51987,51871,51476,51310,51071,50690,50532,
     3  50092,49947,49871,49889,50155,51027,51212,51267,51229,51032,
     4  50556,50320,50033,49846,49761,49362,49203,48966,48838,48759,
     5  48413,48196,48116,47917,47469,47153,46441,46229,45979,45966,
     6  46073,46550,47279,48407,48685,49120,49266,49334,49181,49372,
     7  49333,49160,48930,48676,46883,46352,46007,45941,45689,45609,
     8  45622,45569,45263,44958,44691,44823,46116,46997,47465,47492,
     9  47372,46677,46197,46077,46023,46278,46238,46051,45636,45076/
            data (ix(i),i= 1201, 1300)/
     +  44902,44862,44928,45222,45302,45931,45944,46078,46562,46858,
     1  46954,46912,46764,46468,46051,45716,45461,45019,44898,44804,
     2  44817,45044,45246,45286,45218,44949,44586,44423,44433,44658,
     3  44792,45265,45399,45492,45762,45965,46047,46102,46238,46359,
     4  46318,46399,46943,47311,47530,48129,48520,49018,49311,49406,
     5  49476,49450,48968,48549,48321,48376,48567,48963,48950,48787,
     6  48459,47980,47748,47489,47312,47271,47325,47531,48039,48314,
     7  48356,48268,48068,47944,47806,48081,47888,47558,46928,46751,
     8  46573,46176,46039,45971,46150,46042,45280,45072,44936,44777,
     9  44655,44501,44662,44536,44562,44699,45151,45178,45053,44655/
            data (ix(i),i= 1301, 1400)/
     +  44273,44194,44202,44228,44215,44145,44031,43944,43915,43867,
     1  43895,43874,43822,43915,43845,43721,43476,43355,43174,42960,
     2  42757,42578,42628,42977,43107,43061,43564,43769,44322,44446,
     3  44403,43711,43611,43663,43965,44019,43962,43768,43575,43163,
     4  43045,42936,42702,42271,42216,42142,42131,41991,41702,41501,
     5  41245,40963,40734,40453,40419,40467,40393,40485,40411,40143,
     6  39981,39947,40114,40397,40635,40796,41136,41469,41808,42284,
     7  42672,43164,43349,43494,43455,43527,43862,44064,44232,44270,
     8  44183,43325,43156,42933,42668,42371,42177,41869,41572,41227,
     9  40912,40615,40731,40727,40608,40634,41024,41174,41433,41775/
            data (ix(i),i= 1401, 1500)/
     +  41890,41891,41793,41568,41451,41320,41193,41089,41083,41256,
     1  41168,41028,40351,40139,39547,39355,38978,38646,38623,38773,
     2  38856,38684,38459,38322,38165,37831,37577,37313,36951,36328,
     3  35525,35469,35650,35510,35035,34802,34444,34078,33923,34003,
     4  34302,34175,34187,34039,34063,34175,34476,34743,34986,35562,
     5  35803,35974,36020,35942,35677,35050,34604,34596,34717,35224,
     6  35548,35754,35759,35582,35534,35664,35947,36201,36419,36560,
     7  36445,36151,36070,36104,36512,36980,37345,37615,37778,37967,
     8  38005,37732,37294,37013,36831,36795,36730,36554,36553,36691,
     9  36794,36784,36150,35646,35500,35460,35595,35750,36077,36127/
            data (ix(i),i= 1501, 1600)/
     +  36131,36140,36150,36160,36168,36181,36209,36226,36242,36256,
     1  36259,36255,36247,36234,36223,36213,36204,36198,36156,36058,
     2  36057,36055,36053,36052,36046,36042,36039,36035,36034,36031,
     3  36028,36024,36016,36003,35989,35991,35995,35999,35998,35994,
     4  35998,36001,36014,36023,36026,36017,36009,36008,36008,36004,
     5  36009,36008,36003,35996,35994,35983,35966,35964,35955,35951,
     6  35951,35949,35941,35926,35915,35909,35897,35883,35872,35842,
     7  35831,35823,35820,35813,35813,35807,35798,35789,35788,35790,
     8  35788,35781,35780,35783,35780,35777,35772,35764,35760,35769,
     9  35777,35776,35765,35760,35757,35756,35748,35743,35734,35723/
            data (ix(i),i= 1601, 1700)/
     +  35716,35715,35709,35705,35704,35696,35695,35689,35690,35688,
     1  35673,35664,35660,35659,35665,35669,35671,35672,35673,35675,
     2  35674,35671,35665,35654,35653,35650,35650,35646,35647,35641,
     3  35636,35635,35638,35639,35627,35616,35613,35612,35616,35624,
     4  35633,35638,35646,35645,35645,35655,35672,35675,35676,35675,
     5  35675,35671,35663,35661,35671,35686,35688,35695,35706,35713,
     6  35717,35720,35712,35711,35712,35716,35723,35727,35721,35718,
     7  35715,35716,35721,35717,35712,35707,35689,35672,35659,35643,
     8  35620,35613,35586,35569,35534,35525,35522,35534,35541,35554,
     9  35567,35578,35582,35583,35584,35583,35576,35570,35559,35555/
            data (ix(i),i= 1701, 1800)/
     +  35537,35528,35525,35527,35524,35523,35519,35521,35528,35531,
     1  35533,35534,35532,35553,35765,35776,35781,35786,35788,35788,
     2  35812,35828,35836,35842,35851,35854,35867,35872,35892,35908,
     3  35916,35927,35933,35937,35947,35945,35949,35960,35968,35983,
     4  35993,35991,35992,35993,35998,36012,36026,36041,36048,36059,
     5  36079,36092,36127,36164,36132,36113,36121,36124,36136,36156,
     6  36200,36150,36091,36073,35972,35806,35233,34840,34772,34512,
     7  34414,34346,34049,33856,33723,33635,33249,33234,33083,33033,
     8  32969,32885,32685,32145,31940,31917,31877,31405,30342,30003,
     9  29871,28847,28743,28560,28689,27694,27703,27698,27918,27919/
            data (ix(i),i= 1801, 1900)/
     +  27980,26695,26723,26715,27170,27174,26747,26692,25945,25583,
     1  25127,24711,24136,23889,23873,23841,23488,23382,23099,22719,
     2  22357,22377,22350,22346,22655,21403,21023,20769,19055,18817,
     3  18885,18888,18535,18642,18384,18272,17949,17451,17382,16354,
     4  16083,16093,16063,16167,16179,15792,15820,15788,15888,15895,
     5  15983,16075,16075,16070,15721,15638,15603,15567,15528,15512,
     6  15322,15254,15241,15211,15204,15162,14267,13285,12757,12738,
     7  12390,12152,11410,11144,10768,11222,11225,11234,11234,11071,
     8  10963,10738,10707,10272,10302,10193,10223,10207,10211, 9851,
     9   8910, 8413, 8047, 8045, 7803, 7805, 7386, 7071, 6373, 6088/
            data (ix(i),i= 1901, 1920)/
     +   6018, 6021, 6024, 6023, 6135, 6137, 5756, 2686, 1643, 1330,
     1     63,    0,  151, 4158, 7251, 8166,12173,16180,16515,16697/
            data (iy(i),i=    1,  100)/
     +  14254,12799,12005,11759,11024, 9838, 9605, 9429, 8873, 8310,
     1   7667, 6816, 6460, 5152, 4729, 4231, 3624, 3191, 2750, 2405,
     2   1966, 1860, 1863, 2109, 3044, 3118, 3055, 2699, 1949, 1087,
     3    750,  554,  386,  383,  585,  451,    9,    0,  110,  324,
     4    406,  306,  379, 1014, 1604, 2119, 1940, 2093, 3136, 3630,
     5   3876, 3942, 4184, 4482, 4677, 5212, 5812, 6613, 7270, 7794,
     6   7751, 7599, 7220, 6903, 6696, 6106, 5790, 5728, 5775, 5642,
     7   5787, 5708, 5571, 5391, 5449, 5851, 5932, 5787, 5297, 5193,
     8   5366, 5754, 5971, 6329, 6663, 6800, 6880, 6703, 6799, 7235,
     9   7372, 7445, 7263, 7184, 7222, 7531, 7953, 8262, 8421, 8204/
            data (iy(i),i=  101,  200)/
     +   7915, 7953, 8160, 8951, 9142, 9097, 8826, 8697, 8890, 9285,
     1   9425, 9199, 9204, 9386, 9254, 8998, 8917, 9008, 9704, 9741,
     2   9676, 9214, 9114, 9153, 9675,10018,10530,10783,10961,10907,
     3  11057,11221,11428,11452,11170,11179,11311,11382,11591,11745,
     4  12147,12477,12532,12467,12041,11959,12015,12694,12841,13060,
     5  13228,13349,13443,13485,13444,13370,13296,13354,14033,14105,
     6  14013,13667,13537,13421,13355,13405,13509,13660,13831,13933,
     7  14009,13953,13841,13614,13452,13259,13072,13006,12996,13226,
     8  13572,13795,14085,14411,15146,15283,15423,15220,15257,15413,
     9  15399,15279,15106,14901,14714,14580,14651,14859,15169,15546/
            data (iy(i),i=  201,  300)/
     +  15633,15380,15398,15537,15097,15219,15442,15647,16058,16384,
     1  16574,17038,17107,17005,17040,17212,17743,18017,18238,18250,
     2  18093,17713,17694,17848,18538,18709,18829,18980,19051,19178,
     3  19310,19362,19375,19334,19294,19322,19392,19455,19516,19522,
     4  19458,19374,19272,19196,19242,19248,19159,19077,19010,18878,
     5  18807,18782,18631,18531,18326,18327,18584,18739,18859,18705,
     6  18139,18191,19117,19236,19768,20041,20744,21001,21088,21175,
     7  21536,21587,20629,20157,19669,19105,18304,18311,18067,17951,
     8  17693,17409,16856,16089,14616,41722,41841,41790,41127,40838,
     9  40473,40179,40151,40345,40511,40202,39903,39368,38743,38279/
            data (iy(i),i=  301,  400)/
     +  38213,38047,37744,37123,36142,35946,35514,35039,33849,33705,
     1  33780,34021,34170,34141,34011,33651,33381,33360,33338,33235,
     2  33130,33063,33047,33072,33097,33119,33025,32931,32847,32784,
     3  32752,32844,33193,33051,32799,32643,32554,32583,34107,34316,
     4  34615,35013,35195,35792,36114,36455,37140,37306,37563,37667,
     5  37779,37864,37882,37973,38112,38183,38202,38123,38035,37944,
     6  37851,37851,37975,38535,38974,39260,39352,39449,39484,39535,
     7  39672,39785,39886,39919,40008,40097,40093,40046,40037,40079,
     8  40150,40294,40332,40459,40575,40675,40848,40861,41218,41385,
     9  41459,99737,99737,99737,99737,99737,99737,99737,99737,99737/
            data (iy(i),i=  401,  500)/
     +  99737,99737,99742,99742,99737,99738,99741,99743,99744,99745,
     1  99745,99745,99745,99743,99737,99737,99737,98948,96158,95674,
     2  91982,91179,88901,85941,82620,82078,76769,74062,73553,66780,
     3  65504,64576,62635,59930,56946,54731,51323,49193,48387,46734,
     4  46723,46580,46534,46537,46465,46220,46159,45970,45881,45926,
     5  45812,45821,45868,45865,45698,45480,45385,45213,45164,45031,
     6  44906,44926,44953,45080,45138,45146,44978,44829,44658,44569,
     7  44325,44286,44161,44215,44130,44007,43868,43744,43416,43134,
     8  43118,43040,42949,42971,42924,42966,43102,43139,43142,43015,
     9  42949,42931,42906,42942,42910,42732,42592,42470,42427,42164/
            data (iy(i),i=  501,  600)/
     +  41670,41332,41354,41493,41506,41435,41434,41324,41247,41215,
     1  41010,40839,40847,40944,41146,41213,41399,41697,41821,41849,
     2  41806,41690,41559,41336,41279,41103,41001,40932,40872,40810,
     3  40737,40469,40440,40209,40236,40157,40253,40291,40250,40105,
     4  40147,40127,40003,39935,39823,39668,39584,39415,39331,39114,
     5  38996,38897,38845,38839,38887,38989,38958,39028,39007,38789,
     6  38699,38563,38572,38540,38462,38176,38101,37902,37711,37572,
     7  37532,37424,37168,37101,36896,36724,36618,36515,36402,36263,
     8  36193,36154,36124,36052,35875,35780,35591,35455,35416,35392,
     9  35337,35169,35014,34959,34765,34733,34700,34708,34810,34833/
            data (iy(i),i=  601,  700)/
     +  34804,34784,34858,34987,35188,35292,35410,35540,35586,35518,
     1  35391,35177,34932,34886,34805,34723,34748,34718,34458,34425,
     2  34174,34126,34044,33950,33834,33610,33391,33131,33113,33039,
     3  33089,33058,32762,32490,32459,32476,32665,32606,32615,32632,
     4  32772,32810,32985,33039,33019,32463,32183,32098,31906,31682,
     5  31600,31551,31516,31375,31291,31082,31091,30956,30860,30603,
     6  30426,30212,29932,29886,29659,29554,29330,29230,29016,28917,
     7  28788,28762,28760,28745,28894,29041,29127,29424,29615,29558,
     8  29633,29640,29600,29419,29290,29211,29093,28979,28811,28761,
     9  28662,28458,28240,28049,27831,27621,27572,27584,27493,27427/
            data (iy(i),i=  701,  800)/
     +  27219,27111,26955,26850,26765,26675,26619,26406,26335,26107,
     1  25985,25717,25643,25473,25392,25296,25277,25261,25213,25044,
     2  24998,24923,24925,24844,24806,24663,24629,24382,24380,24246,
     3  24170,23916,23750,23751,23396,23403,23451,23399,23102,23056,
     4  22888,22859,22711,22684,22547,22297,21989,21921,21828,21724,
     5  21580,21512,21420,21321,21331,21296,21126,21032,20893,20809,
     6  20740,20613,20587,20465,20344,20242,20255,20372,20284,20186,
     7  20037,19961,19789,19720,19587,19509,19277,19146,18987,18863,
     8  18686,18691,18869,18930,19047,19065,19134,19068,19106,18938,
     9  18539,18370,18041,17873,17804,17561,17501,17497,17423,17266/
            data (iy(i),i=  801,  900)/
     +  17211,17170,17109,17008,16924,16736,16685,16482,16363,15755,
     1  15641,15393,15104,14828,14611,14452,14015,13561,13430,13235,
     2  13013,12765,12542,12190,12112,11817,11617,11385,11313,11078,
     3  10929,10796,10602,10432,10297,10265,10246,10339,10316,10256,
     4   9724, 9472, 9181, 9116, 8905, 8804, 8781, 8700, 8645, 8600,
     5   8551, 8301, 8239, 8060, 7917, 7822, 7728, 7679, 7556, 7387,
     6   7318, 7304, 7130, 7125, 7164, 7020, 6903, 6810, 6693, 6232,
     7   6061, 5947, 5846, 5730, 5687, 5553, 5557, 5562, 5561, 5560,
     8   5537, 5552, 5594, 5588, 5580, 5583, 5584, 5583, 5577, 5581,
     9   5584, 5580, 5565, 5564, 5558, 5564, 5537, 5535, 5583, 5591/
            data (iy(i),i=  901, 1000)/
     +   5592, 5601, 5679, 6150, 6169, 6021, 5920, 5655, 5614, 5610,
     1   5681, 6201, 6675, 6869, 7362, 7929, 7988, 7919, 8041, 8290,
     2   8636, 9469, 9618, 9560, 8914, 8624, 8337, 8350, 8412, 8677,
     3   9380, 9750,10413,10483,10671,10849,11024,11263,11415,11531,
     4  11553,11505,11406,11287,11102,10882,10719,10542,10417,10144,
     5   9877, 9381, 9090, 9240, 9546, 9903,10015,10516,10845,11093,
     6  11452,11872,11933,11872,10940,10302,10240,10772,11029,11209,
     7  11325,11437,11303,11375,12423,12727,13071,13851,13987,14205,
     8  14467,14694,14917,15272,15789,15850,15738,15597,15217,14968,
     9  14837,14433,14099,13891,13233,13041,12803,12798,13259,13299/
            data (iy(i),i= 1001, 1100)/
     +  13218,13026,12879,12474,12463,12096,12357,12706,13695,14209,
     1  14405,14577,14620,14396,14368,14469,14791,15018,15389,15868,
     2  16294,16603,17155,17409,17593,17763,17746,17915,18308,18546,
     3  18575,18136,17997,18018,17883,17368,17173,17259,17736,17830,
     4  17700,17446,17348,17579,18020,19108,20003,20660,21133,21473,
     5  21774,22032,21974,21825,21133,20160,19495,18548,18078,17983,
     6  18045,18325,18425,18792,18830,18735,18324,18140,18060,18136,
     7  18347,19257,19563,20015,20276,20315,20184,19678,18515,18215,
     8  18119,18249,18459,18765,18873,18422,18473,18397,18422,18478,
     9  18652,19038,19443,19596,19638,19384,19424,19649,19963,19999/
            data (iy(i),i= 1101, 1200)/
     +  20354,20820,21368,21508,21712,21942,22248,22974,23499,23537,
     1  23056,22499,22193,21708,21504,20969,20602,20409,20399,20155,
     2  20294,20701,20944,21062,21204,21219,21586,21439,21590,21578,
     3  21265,21303,21512,21907,22071,22008,22037,22155,22259,22352,
     4  22402,22513,22848,22733,22254,22007,21994,22103,22536,22572,
     5  22357,21709,21573,21491,21637,21866,22186,22394,22964,23170,
     6  23237,23111,23100,23250,23210,22924,22972,23142,23731,24276,
     7  24380,24315,23874,23674,23450,23543,23839,24011,24151,24101,
     8  23809,23689,23709,23985,24502,24758,25023,25320,25656,25759,
     9  25898,25548,25502,25555,25693,25982,26086,26139,25920,25478/
            data (iy(i),i= 1201, 1300)/
     +  25496,25599,25770,25957,26266,26707,27531,27959,28263,28362,
     1  28575,28653,28689,28504,28371,28117,27569,27332,27384,27625,
     2  27762,27950,27948,28017,28189,28448,28432,28844,29497,30543,
     3  30714,30970,31140,31758,32443,32613,32578,32183,32181,32472,
     4  32747,32970,33393,33835,33969,33736,32733,32052,31289,31304,
     5  31560,31802,32792,33351,34111,34316,34433,34459,34666,34738,
     6  34710,34460,34447,34588,35123,35570,35776,36116,36641,37152,
     7  37564,38433,38720,38773,38690,37912,37108,36563,35918,35250,
     8  35183,35410,35325,35154,34706,34381,33905,34421,34335,33546,
     9  33512,34062,34594,34989,35091,35177,35073,35176,35399,35538/
            data (iy(i),i= 1301, 1400)/
     +  35298,34989,34795,34705,34503,34307,34227,34301,34416,34556,
     1  34654,34740,35228,35503,35675,35726,35571,35364,34385,33990,
     2  33817,33919,34296,34865,35466,35912,37613,37785,37615,37649,
     3  37838,38094,38334,38626,38816,38936,39073,39124,39089,38829,
     4  38365,38210,38191,38394,38528,38701,39387,40295,40978,41284,
     5  41555,41722,42010,42742,42982,43360,43530,43806,43976,44022,
     6  44311,44550,44656,44576,44598,44309,44177,43651,43518,43558,
     7  43802,43720,43533,43190,42950,42779,42814,42335,42266,42593,
     8  42884,44168,44235,44096,44059,44177,44054,44034,44150,44539,
     9  44774,45455,45920,46092,46364,46467,46715,46940,47039,46975/
            data (iy(i),i= 1401, 1500)/
     +  47012,47126,47196,47206,47209,47245,47304,47539,47796,48297,
     1  48415,48327,47028,46870,46428,46217,45539,45325,45135,44797,
     2  44370,43474,43502,43671,43925,44618,44813,45114,45839,46656,
     3  48013,49074,49741,49847,49929,50348,51343,51497,51747,52196,
     4  52570,52839,53199,53586,53690,53764,53521,53036,52756,51374,
     5  51093,51117,51410,51578,52065,52930,53870,54041,54252,54583,
     6  55454,55720,55943,56004,56105,56453,56636,56492,56501,56789,
     7  56913,56987,57156,57397,58099,59146,60257,60472,60443,60518,
     8  60673,60853,60752,60794,61250,61712,61847,61823,61514,61073,
     9  60391,59945,57999,57190,57184,57422,57994,58446,59008,59218/
            data (iy(i),i= 1501, 1600)/
     +  59222,59225,59235,59248,59255,59266,59286,59307,59341,59371,
     1  59387,59428,59442,59472,59492,59512,59531,59541,59852,60055,
     2  60072,60101,60133,60158,60186,60196,60208,60235,60272,60300,
     3  60325,60351,60369,60368,60367,60374,60403,60415,60440,60464,
     4  60483,60489,60496,60498,60505,60535,60557,60590,60623,60659,
     5  60705,60723,60741,60753,60771,60795,60811,60825,60847,60867,
     6  60886,60896,60910,60933,60953,60972,60987,61003,61017,61056,
     7  61072,61084,61090,61106,61122,61140,61152,61157,61170,61198,
     8  61217,61255,61274,61301,61313,61329,61357,61373,61390,61408,
     9  61405,61423,61441,61443,61451,61465,61502,61510,61511,61515/
            data (iy(i),i= 1601, 1700)/
     +  61526,61535,61553,61567,61573,61579,61599,61625,61632,61646,
     1  61659,61671,61677,61687,61714,61723,61741,61774,61799,61825,
     2  61834,61846,61854,61865,61875,61881,61894,61910,61928,61940,
     3  61954,61969,61985,61990,62007,62027,62035,62047,62054,62064,
     4  62073,62094,62121,62148,62160,62165,62170,62186,62199,62207,
     5  62221,62247,62300,62319,62340,62357,62367,62382,62407,62422,
     6  62430,62447,62471,62491,62518,62541,62564,62585,62613,62648,
     7  62668,62697,62744,62768,62789,62809,62841,62869,62897,62920,
     8  62950,62972,63003,63019,63046,63057,63067,63076,63089,63102,
     9  63136,63171,63192,63220,63253,63276,63304,63328,63342,63350/
            data (iy(i),i= 1701, 1800)/
     +  63359,63367,63381,63394,63408,63424,63438,63447,63459,63468,
     1  63476,63490,63500,63566,63948,63946,63951,63957,63972,63989,
     2  64032,64027,64025,64024,64028,64039,64058,64073,64119,64159,
     3  64175,64203,64224,64236,64259,64288,64298,64309,64328,64349,
     4  64364,64401,64426,64432,64440,64458,64483,64500,64521,64548,
     5  64572,64589,64687,64788,64782,64785,64867,64910,65199,65372,
     6  65447,65652,65819,65847,66090,66502,66333,66592,66615,66708,
     7  67211,67556,67626,67732,67797,67909,68439,68457,68649,68664,
     8  68679,68708,68758,68910,68960,68897,68909,69385,70177,70541,
     9  70684,70590,71493,71937,72347,72949,72987,72990,74116,74121/
            data (iy(i),i= 1801, 1900)/
     +  74415,74820,74875,74877,75859,75869,76960,77054,78320,79240,
     1  80039,81231,82219,82618,82672,82726,83839,83994,84424,85049,
     2  85824,85848,85903,85920,86263,87720,88459,88829,90009,90442,
     3  90771,90786,90946,91446,91876,92225,92220,92792,93240,93497,
     4  94039,94071,94134,94404,94430,94754,94814,94840,95043,95058,
     5  95148,95241,95764,96003,96514,96520,96524,96529,96576,96589,
     6  96864,97145,97152,97152,97152,97181,98022,97423,97100,97084,
     7  96860,96811,96655,96516,96308,95960,95785,95673,95669,95356,
     8  95151,95028,95013,95151,93611,93515,93413,93379,93321,92541,
     9  92479,92049,91794,91798,91628,91625,91332,91153,90362,90344/
            data (iy(i),i= 1901, 1920)/
     +  90340,90359,90371,90371,90940,91135,93166,96692,97808,98936,
     1  99625,99703,99737,99737,99743,99737,99737,99737,99739,99740/
c
c       map of british columbia
c       this data was provided by jeff ovall
c
        nvf=1920
        nbf=nvf
        ntf=3
        rp(12)=0.1e0_rknd
        rp(13)=1.5e0_rknd
c
        sp(1)='British Columbia'
        sp(2)='British Columbia'
        sp(3)='British Columbia'
c
        do i=1,nvf-1
            if(ix(i) /= ix(i+1)) cycle
            if(iy(i) /= iy(i+1)) cycle
            write(6,*) 'error', i,ix(i),iy(i)
        enddo
        do i=1,nvf
            vx(i)=real(ix(i),rknd)
            vy(i)=real(iy(i),rknd)
        enddo
c
c       Vancouver Island
c
        do i=1,285
            ibndry(1,i)=i
            ibndry(2,i)=i+1
            ibndry(3,i)=0
            ibndry(4,i)=2
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=1
        enddo
        ibndry(2,285)=1
c
c       Graham Island
c
        do i=286,391
            ibndry(1,i)=i
            ibndry(2,i)=i+1
            ibndry(3,i)=0
            ibndry(4,i)=2
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=2
        enddo
        ibndry(2,391)=286
c
c       Mainland
c
        do i=392,nbf
            ibndry(1,i)=i
            ibndry(2,i)=i+1
            ibndry(3,i)=0
            ibndry(4,i)=2
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=3
        enddo
        ibndry(2,nbf)=392
c
        eps=1.0e-3_rknd
        call fixdg2(nvf,nbf,vx,vy,sf,ibndry,eps)
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
c
c       make itnode
c
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine gd23(vx,vy,sf,itnode,ibndry,ip,rp,sp,iu,ru,su,sxy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), dimension(600), save :: ix,iy
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            character(len=80), dimension(100) :: sp,su
cy
            external sxy
            data (ix(i),i=    1,  100)/
     +  34609,35071,38054,38231,38260,38405,38494,39439,40660,40845,
     1  40951,41434,42059,42748,42911,43039,43075,43728,43827,43991,
     2  44239,44232,44012,44076,44140,44282,44438,44488,44516,44566,
     3  44829,44758,45163,45191,45219,45305,45397,45468,45503,45511,
     4  45511,45553,45788,45986,45993,46050,46057,46079,46086,46093,
     5  46100,46128,46135,46185,46192,46199,46221,46228,46235,46228,
     6  46235,46952,47172,47513,47641,47464,47527,47712,47883,47975,
     7  48223,49495,49502,49509,49502,49097,49211,49225,49331,49353,
     8  49360,49445,49480,49630,49878,49992,50070,50105,50091,50091,
     9  50355,51225,51897,52031,52936,53025,53172,53468,54414,55071/
            data (ix(i),i=  101,  200)/
     +  55398,55603,55905,56335,56576,56961,57109,57759,58461,58645,
     1  59081,59207,59417,59528,59978,60411,60471,60489,60505,60537,
     2  60563,61050,61178,61229,61390,61553,61680,61659,61624,61567,
     3  61579,61797,61946,62186,62431,62491,62671,63276,63375,63469,
     4  63565,63685,63774,64426,64814,65209,65472,65852,66065,66224,
     5  66865,67254,67767,67954,68205,68393,68472,68591,68810,69099,
     6  69199,70606,70739,70774,70682,70414,69906,69688,69609,69669,
     7  69693,69732,69694,69401,69590,69968,70133,70218,70389,70422,
     8  70390,70371,70275,70046,69965,70015,70226,70576,70771,71310,
     9  71250,70841,70793,70824,70851,71140,71169,71134,70997,70709/
            data (ix(i),i=  201,  300)/
     +  70199,70091,70030,70160,70085,70173,70218,70437,70882,70969,
     1  70971,70854,70787,70756,70293,69918,69503,68876,68077,68054,
     2  67892,67432,67334,67069,66941,66843,66738,66156,65787,65382,
     3  65127,64945,64312,64006,63286,63120,63166,62878,62890,62827,
     4  62601,62232,62198,62262,62240,62299,62488,62529,62538,62417,
     5  62400,62663,62638,62381,62337,62219,62070,62016,62032,62204,
     6  62148,61890,61753,61874,62039,62202,62393,62443,62411,62259,
     7  62096,62074,62114,62325,62377,62314,62189,62102,62166,61865,
     8  61939,61895,61704,61921,61751,61120,60367,60045,59464,59081,
     9  58761,58589,58597,58415,58267,58299,58218,57803,57390,57173/
            data (ix(i),i=  301,  400)/
     +  56681,56617,55431,55134,54931,54528,54136,53680,53326,53233,
     1  53079,52943,52801,52474,52384,52078,51575,50938,50174,49910,
     2  49715,49563,49478,49365,49246,49023,48971,48918,48830,48739,
     3  48074,47904,47633,47605,47725,47640,47586,47424,47073,46907,
     4  46610,46566,46578,46926,47065,47258,47289,47256,47208,47000,
     5  46986,46971,47051,47318,47353,47403,47612,47697,47637,47432,
     6  47438,47487,47403,47403,47486,47403,47136,47199,47403,47590,
     7  48174,48298,48071,47403,46795,46480,46693,47187,47403,47894,
     8  48545,49025,49374,49512,50075,50589,51693,52880,53345,53268,
     9  53001,52966,51678,50588,49483,49252,48567,48053,47629,47403/
            data (ix(i),i=  401,  500)/
     +  46849,46616,46501,46220,45996,46106,46011,45500,44933,44436,
     1  44505,44793,45203,45838,46999,47403,48071,49115,50590,51130,
     2  51485,52038,52516,53131,53507,54252,54889,55717,56565,56813,
     3  57470,57873,59562,59559,59554,59551,59549,59458,59037,58070,
     4  57775,57500,57369,57661,57912,58246,58802,59005,59133,59158,
     5  59099,58986,58827,58164,57849,57746,57805,57939,58015,58013,
     6  57873,57570,57495,57070,56629,56282,56085,55926,55918,56046,
     7  55916,55427,55331,55045,54820,54819,54818,54650,54457,53830,
     8  53694,53326,53118,53047,52876,52757,52690,52560,52520,52498,
     9  51920,51744,51442,51391,51290,51363,51580,52175,52309,52228/
            data (ix(i),i=  501,  572)/
     +  51659,51497,51584,51406,51097,51519,51898,49920,49753,49387,
     1  49468,49176,48714,48687,48600,48657,48674,48600,48164,48164,
     2  48231,48238,47903,47446,47403,47057,46782,46610,46542,46377,
     3  46256,44475,44138,44016,43914,44155,44274,44159,43701,43718,
     4  41103,41140,41156,40327,38126,37605,35821,35756,36284,35175,
     5  35249,34509,32257,30798,28689,29839,30033,30977,31260,31254,
     6  31190,31313,31413,32265,32469,33405,34831,35576,35803,35799,
     7  35452,34351/
            data (iy(i),i=    1,  100)/
     +  53333,54036,58922,58709,59029,59000,59299,61252,63943,64462,
     1  64675,65783,67522,69390,69852,70306,70505,72863,73239,73978,
     2  74745,75349,75349,75384,75377,75384,75384,75903,75974,76250,
     3  77046,77174,78772,78885,79148,79439,79951,80306,80540,80703,
     4  81094,81257,81414,82003,82067,82251,82308,82429,82500,82592,
     5  82678,82990,83047,83359,83430,83509,83693,83750,83807,84275,
     6  84332,84730,84595,84460,84453,84439,84346,84183,84112,84133,
     7  84034,86413,86498,86598,86669,86882,87194,87322,87833,87983,
     8  88082,88402,88558,89133,89978,90347,90752,90937,91058,91186,
     9  91204,91214,91262,91272,91139,91172,91306,91386,91413,91490/
            data (iy(i),i=  101,  200)/
     +  91529,91488,91429,90740,90508,90315,90279,90523,90570,90652,
     1  91007,91152,91281,91291,91200,91337,91398,91565,91721,91780,
     2  91829,92003,92259,92817,93438,93678,93963,94273,94407,94622,
     3  94821,95239,95343,95609,96062,96117,95993,95268,95242,95318,
     4  95520,95653,95688,95945,96140,96339,96552,97062,97224,97275,
     5  97230,97286,97465,97598,97765,97961,98163,98814,99108,99285,
     6  99345,100000,99914,99731,99406,99034,98548,98073,97610,97108,
     7  96900,96573,96399,95977,95786,95516,95322,95198,94678,94377,
     8  94249,94174,93966,93757,93620,93343,93279,93248,93171,92516,
     9  91846,91541,91434,91359,91295,91019,90909,90846,90778,90749/
            data (iy(i),i=  201,  300)/
     +  90863,90790,90656,90208,89403,89162,89038,87336,86815,86636,
     1  86496,86223,86125,86080,85398,84701,84326,83376,82592,82677,
     2  82755,82589,82560,82534,82442,82251,81903,81837,81522,81278,
     3  81213,81112,81074,80908,80642,80482,80122,80025,79869,79807,
     4  79754,79517,79254,79191,78795,78677,78589,78482,78430,78191,
     5  77951,77490,77437,77440,77199,77201,77284,77259,77107,76871,
     6  76790,76744,76581,76504,76508,76192,76126,76057,75927,75679,
     7  75445,75315,75200,75113,75026,74901,74897,74830,74492,74412,
     8  73835,73740,73693,73301,73315,73463,73761,73844,73866,73771,
     9  73813,73943,74072,74271,74641,75152,75700,76295,76557,76602/
            data (iy(i),i=  301,  400)/
     +  76609,76599,76414,76416,76462,76629,76874,77274,77438,77425,
     1  77405,77348,77211,76898,76742,76491,76180,75906,75661,75576,
     2  75452,75209,75074,74718,74063,73540,73221,72891,72819,72744,
     3  72491,72309,71876,71187,70455,70265,70233,70134,69702,69385,
     4  68813,68515,68218,67775,67528,66799,66471,66361,66203,65869,
     5  65767,65665,65495,65253,65195,65114,64776,64482,64105,63724,
     6  63552,63414,63356,63270,63084,62932,62445,61866,61631,61416,
     7  61019,60609,60012,59910,59817,59419,58734,58554,58527,58467,
     8  58666,59111,59287,59357,58913,58702,58322,58328,57912,56987,
     9  56630,56583,56554,56174,55239,55043,54768,54551,54234,54025/
            data (iy(i),i=  401,  500)/
     +  53513,53181,53015,52675,51715,50673,50227,49395,48767,47912,
     1  47402,47005,46713,46567,46698,46784,46927,46852,46771,47011,
     2  46941,47023,47058,47157,47350,48006,48620,48935,49467,49654,
     3  50151,50273,50252,49809,49215,48834,48501,47986,46956,45465,
     4  44950,44535,44019,43474,43183,42797,42154,41779,41246,40847,
     5  40620,40191,39723,38828,38424,37903,37487,37170,36989,36461,
     6  36350,36111,36068,35828,35578,34944,34568,34267,33804,33382,
     7  32867,32715,32615,32317,31784,31294,30741,30403,30015,29025,
     8  28456,28161,27824,27706,27419,27221,27109,26892,26616,26470,
     9  24149,23691,23315,23252,22906,22791,22449,21464,21065,20760/
            data (iy(i),i=  501,  572)/
     +  20215,19775,19265,18923,18327,16705,15249,12780,12203,10939,
     1   8939, 8430, 7625, 7508, 7144, 6853, 6762, 6246, 4662, 3800,
     2   3260, 2762, 2404, 2110, 2067, 1717, 1235,  807,  637,  631,
     3    688,    0,  556,  984, 1729, 2427, 3126, 3606, 4303, 5535,
     4  11933,12203,12320,14311,19721,21419,22432,23217,24509,25635,
     5  27031,29694,33779,37462,43486,44228,44371,45072,45516,46205,
     6  46942,47296,47583,48523,48692,49465,50653,50805,51151,51642,
     7  52083,52958/
c
c       map of israel
c       this data was provided by jeff ovall
c
        sp(1)='Israel'
        sp(2)='Israel'
        sp(3)='Israel'
c
        nvf=572
        nbf=nvf
        ntf=1
        rp(12)=0.1e0_rknd
        rp(13)=1.5e0_rknd
c
        do i=1,nvf
            vx(i)=real(ix(i),rknd)
            vy(i)=real(iy(i),rknd)
        enddo
c
        do i=1,nbf
            ibndry(1,i)=i
            ibndry(2,i)=i+1
            ibndry(3,i)=0
            ibndry(4,i)=2
            ibndry(5,i)=0
            ibndry(6,i)=1
        enddo
        ibndry(2,nbf)=1
c
        eps=1.0e-3_rknd
        call fixdg2(nvf,nbf,vx,vy,sf,ibndry,eps)
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
c
c       make itnode
c
        call sklutl(0_iknd,vx,vy,sf,itnode,ibndry,ip,rp,iflag,sxy)
        return
        end
