c***********************  problem name: naca    ************************
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),rminf,angle,
     +          uinf,gamma,theta,eps,ru(94),su(100)
cy
        call rho(ux,uy,r,rx,ry,rxx,rxy,ryy)
c
        values(k0)=ux*(r*theta+(1.0e0_rknd-theta))
        values(kx)=(r+ux*rx)*theta+(1.0e0_rknd-theta) 
        values(ky)=ux*ry*theta
        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),rminf,angle,
     +          uinf,gamma,theta,eps,ru(94),su(100)
cy
        call rho(ux,uy,r,rx,ry,rxx,rxy,ryy)
c
        values(k0)=uy*(r*theta+(1.0e0_rknd-theta))
        values(ky)=(r+uy*ry)*theta+(1.0e0_rknd-theta) 
        values(kx)=uy*rx*theta
        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
            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 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
            character(len=80) :: su
            common /val1/k0,ku,kl
            common /atest2/iu(100),rminf,angle,
     +          uinf,gamma,theta,eps,ru(94),su(100)
cy
        if(rl/=rminf) then
            rminf=rl
            call cuinf
        end if
c
        if(itag<=0) return
        pi=3.141592653589793e0_rknd
        ang=(1.0e0_rknd/8.0e0_rknd+real(itag-1,rknd)/4.0e0_rknd
     +      -angle/180.0e0_rknd)*pi
        cc=cos(ang)
        call rho(uinf,0.0e0_rknd,r,rx,ry,rxx,rxy,ryy)
        values(k0)=r*uinf*cc
        values(kl)=(rx*uinf+r)*cc
        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
            character(len=80) :: su
            common /val2/k0,kl,klb,kub,kic,kim,kil
            common /atest2/iu(100),rminf,angle,
     +          uinf,gamma,theta,eps,ru(94),su(100)
cy
        call rho(uinf,0.0e0_rknd,r,rx,ry,rxx,rxy,ryy)
        values(kic)=r*uinf*x
        values(kil)=rminf
        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
            character(len=80) :: su,sp
            common /val0/k0,ku,kx,ky,kl
cy
c
        values(k0)=u**2
        values(ku)=2.0e0_rknd*u
        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
            character(len=80) :: su,sp
            common /val0/k0,ku,kx,ky,kl
cy
        if(itag>0) return
        values(k0)=ux**2+uy**2
        values(kx)=ux*2.0e0_rknd
        values(ky)=uy*2.0e0_rknd
        return
c
        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
            character(len=80) :: su
            common /val3/kf,kf1,kf2,kad
            common /atest2/iu(100),rminf,angle,
     +          uinf,gamma,theta,eps,ru(94),su(100)
cy
        ss=ux**2+uy**2
        uu=1.0e0_rknd-ss
        if(uu>0.0e0_rknd) then
            c=ss/uu
        else
            c=1.0e20_rknd
        endif
        rmach=sqrt(2.0e0_rknd*c/(gamma-1.0e0_rknd))
        values(kf)=rmach
        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 rho(ux,uy,r,rx,ry,rxx,rxy,ryy)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            character(len=80) :: su
            common /atest2/iu(100),rminf,angle,
     +          uinf,gamma,theta,eps,ru(94),su(100)
cy
        uu=1.0e0_rknd-ux**2-uy**2
        uux=-2.0e0_rknd*ux
        uuy=-2.0e0_rknd*uy
        uuxx=-2.0e0_rknd
        uuyy=-2.0e0_rknd
        uuxy=0.0e0_rknd
        al=1.0e0_rknd/(gamma-1.0e0_rknd)
c
        if(uu>eps) then
            r=uu**al
            r1=al*r/uu
            r2=(al-1.0e0_rknd)*r1/uu
            rx=r1*uux
            ry=r1*uuy
            rxx=r2*uux**2+r1*uuxx
            ryy=r2*uuy**2+r1*uuyy
            rxy=r2*uux*uuy+r1*uuxy
        else
            ff=eps*exp((uu-eps)/eps)
            fx=ff*uux/eps
            fy=ff*uuy/eps
            fxx=ff*((uux/eps)**2+uuxx/eps)
            fyy=ff*((uuy/eps)**2+uuyy/eps)
            fxy=ff*((uux*uuy)/eps**2+uuxy/eps)
c
            r=ff**al
            r1=al*ff**(al-1.0e0_rknd)
            r2=(al-1.0e0_rknd)*al*ff**(al-2.0e0_rknd)
            rx=r1*fx
            ry=r1*fy
            rxx=r2*fx**2+r1*fxx
            ryy=r2*fy**2+r1*fyy
            rxy=r2*fx*fy+r1*fxy
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        function wing(x)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            real(kind=rknd) :: wing
cy
        z=x+0.5e0_rknd
        t=0.12e0_rknd
        q=(((.1015e0_rknd*z-.2843e0_rknd)*z+.3516e0_rknd)*z+.126e0_rknd)*z
        wing=5.0e0_rknd*t*(0.2969e0_rknd*sqrt(z)-q)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine cuinf
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            character(len=80) :: su
            common /atest2/iu(100),rminf,angle,
     +          uinf,gamma,theta,eps,ru(94),su(100)
cy
        a=2.0e0_rknd/(gamma-1.0e0_rknd)
        c2=rminf**2+a
        c=sqrt(c2)
        uinf=rminf/c
c
        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
        ii=iu(1)
        s=ru(7)
c
        call usrset(iu,ru,su)
c
        ru(5)=max(1.0e0_rknd,ru(5))
        ru(5)=min(100e0_rknd,ru(5))
c
        ic=0
        if(ii/=iu(1)) ic=1
        if(s/=ru(7)) ic=1
        if(ic==0) then
            call cuinf
        else
            iu(1)=max(1,iu(1))
            iu(1)=min(4,iu(1))
            ip(41)=-1
        endif
        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=minf  ,a= m,t=r',
     2  'n i= 2,n=angle ,a= a,t=r',
     3  'n i= 5,n=theta,a=t,t=r',
     4  'n i= 6,n=eps   ,a= e,t=r',
     5  'n i= 7,n=size  ,a= s,t=r',
     6  's n=domain,v=1,l="naca 0012"',
     7  's n=domain,v=2,l="bi naca 0012"',
     8  's n=domain,v=3,l="three element airfoil"',
     9  's n=domain,v=4,l="three element airfoil"'/
c
            data len0/10/
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 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 :: ispd,iprob,itask,ising,ifirst
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), save :: size,ratio,theta,eps,gamma,angle
            real(kind=rknd), save :: rminf
            character(len=80), dimension(100) :: sp,su
cy
            external sxy
            data ispd,iprob,itask,ising,ifirst/1,1,0,1,1/
            data gamma,angle,rminf/1.4e0_rknd,0.0e0_rknd,0.72e0_rknd/
            data theta,eps/0.5e0_rknd,1.0e-2_rknd/
            data size,ratio/8.0e0_rknd,1.0e0_rknd/
c
        if(ip(41)==1) then
            sp(4)='naca'
            sp(5)='naca_figxxx.ext'
            sp(6)='naca_mpixxx.rw'
            sp(7)='naca.jnl'
            iu(1)=1
            ru(1)=rminf
            ru(2)=angle
            ru(4)=gamma
            ru(5)=theta
            ru(6)=eps
            ru(7)=size
        endif
c
        call cuinf
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)
c
        ip(5)=max(ip(5),ifirst)
        ip(6)=iprob
        ip(12)=ising
        ip(7)=itask
        ip(8)=ispd
        ip(20)=5
        sp(6)='naca_mpixxx.rw'
        sp(7)='naca.jnl'
        sp(9)='naca_mpixxx.out'
c
c
        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), save, dimension(16) :: 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), dimension(16) :: xw
            real(kind=rknd), save, dimension(16) :: yw
            real(kind=rknd), save :: hmax,grade,xy
            character(len=80), dimension(100) :: sp,su
cy
            external sxy
            data ntf,nvf,nbf/2,42,44/
            data hmax,grade/0.1e0_rknd,1.75e0_rknd/
c
c       definition of naca0012
c
            data (ix(i),i=   1,  16)/
     +         351381,   1346558,   3066379,   5517089,   8673215,
     1       12514616,  17041744,  22262234,  28171958,  34769244,
     2       42050736,  50013380,  58655040,  67974792,  77972808,
     3       88648616/
            data (iy(i),i=   1,  16)/
     +         919038,   1914215,   2881318,   3721619,   4423213,
     1        5026549,   5524749,   5855709,   6009519,   5966640,
     2        5719906,   5277705,   4655302,   3855282,   2845609,
     3        1568145/
c
c
        sp(2)='naca 0012'
        sp(1)='naca 0012'
        sp(3)='naca 0012'
c
        rp(1)=ru(1)
        pi=3.141592653589793e0_rknd
c
        do i=1,16
            xw(i)=real(ix(i),rknd)*1.0e-8_rknd
            yw(i)=real(iy(i),rknd)*1.0e-8_rknd
        enddo
c
        size=ru(7)
        vx(1)=-0.5e0_rknd
        vy(1)=0.0e0_rknd
        vx(2)=0.5e0_rknd
        vy(2)=0.0e0_rknd
        do i=3,10
            arg=pi*real(i-3,rknd)/4.0e0_rknd
            vx(i)=size*cos(arg)
            vy(i)=size*sin(arg)
        enddo
        do i=1,16
            vx(10+i)=xw(i)-.5e0_rknd
            vy(10+i)=yw(i)
            vx(26+i)=xw(i)-.5e0_rknd
            vy(26+i)=-yw(i)
        enddo
c
c
        do i=1,8
            ibndry(1,i)=i+2
            ibndry(2,i)=i+3
            ibndry(3,i)=0
            ibndry(4,i)=1
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=i
        enddo
        ibndry(2,8)=3
        ibndry(1,9)=7
        ibndry(2,9)=1
        ibndry(3,9)=0
        ibndry(4,9)=0
        ibndry(7,9)=0
        ibndry(1,10)=2
        ibndry(2,10)=3
        ibndry(3,10)=0
        ibndry(4,10)=0
        ibndry(7,10)=0
        do i=11,27
            ibndry(1,i)=i-1
            ibndry(2,i)=i
            ibndry(3,i)=0
            ibndry(4,i)=1
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=-1
        enddo
        ibndry(1,11)=1
        ibndry(2,27)=2
        do i=28,44
            ibndry(1,i)=i-2
            ibndry(2,i)=i-1
            ibndry(3,i)=0
            ibndry(4,i)=1
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=-2
        enddo
        ibndry(1,28)=1
        ibndry(2,44)=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, 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-----------------------------------------------------------------------
        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, dimension(16) :: 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, dimension(16) :: xw,yw
            real(kind=rknd), save :: hmax,grade
            character(len=80), dimension(100) :: sp,su
cy
            external sxy
            data ntf,nvf,nbf/4,78,83/
            data hmax,grade/0.1e0_rknd,1.75e0_rknd/
c
c       definition of naca0012
c
            data (ix(i),i=   1,  16)/
     +         351381,   1346558,   3066379,   5517089,   8673215,
     1       12514616,  17041744,  22262234,  28171958,  34769244,
     2       42050736,  50013380,  58655040,  67974792,  77972808,
     3       88648616/
            data (iy(i),i=   1,  16)/
     +         919038,   1914215,   2881318,   3721619,   4423213,
     1        5026549,   5524749,   5855709,   6009519,   5966640,
     2        5719906,   5277705,   4655302,   3855282,   2845609,
     3        1568145/
c
        sp(2)='bi naca 0012'
        sp(1)='bi naca 0012'
        sp(3)='bi naca 0012'
c
        rp(1)=ru(1)
        size=ru(7)
c
        pi=3.141592653589793e0_rknd
c
        do i=1,16
            xw(i)=real(ix(i),rknd)*1.0e-8_rknd
            yw(i)=real(iy(i),rknd)*1.0e-8_rknd
        enddo
c
        vx(1)=-0.5e0_rknd
        vy(1)=0.0e0_rknd
        vx(2)=0.5e0_rknd
        vy(2)=0.0e0_rknd
        vx(75)=vx(1)
        vy(75)=0.25e0_rknd
        vx(76)=vx(2)
        vy(76)=0.25e0_rknd
        vx(77)=vx(1)
        vy(77)=-0.25e0_rknd
        vx(78)=vx(2)
        vy(78)=-0.25e0_rknd
        do i=3,10
            arg=pi*real(i-3,rknd)/4.0e0_rknd
            vx(i)=size*cos(arg)
            vy(i)=size*sin(arg)
        enddo
        do i=1,16
            vx(10+i)=xw(i)-.5e0_rknd
            vy(10+i)=yw(i)+.25e0_rknd
            vx(26+i)=xw(i)-.5e0_rknd
            vy(26+i)=-yw(i)+.25e0_rknd
            vx(42+i)=xw(i)-.5e0_rknd
            vy(42+i)=yw(i)-.25e0_rknd
            vx(58+i)=xw(i)-.5e0_rknd
            vy(58+i)=-yw(i)-.25e0_rknd
        enddo
c
c
        do i=1,8
            ibndry(1,i)=i+2
            ibndry(2,i)=i+3
            ibndry(3,i)=0
            ibndry(4,i)=1
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=i
        enddo
        ibndry(2,8)=3
        ibndry(1,9)=7
        ibndry(2,9)=1
        ibndry(3,9)=0
        ibndry(4,9)=0
        ibndry(7,9)=0
        ibndry(1,10)=2
        ibndry(2,10)=3
        ibndry(3,10)=0
        ibndry(4,10)=0
        ibndry(7,10)=0
        do i=11,27
            ibndry(1,i)=i-1
            ibndry(2,i)=i
            ibndry(3,i)=0
            ibndry(4,i)=1
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=-1
        enddo
        ibndry(1,11)=75
        ibndry(2,27)=76
        do i=28,44
            ibndry(1,i)=i-2
            ibndry(2,i)=i-1
            ibndry(3,i)=0
            ibndry(4,i)=1
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=-2
        enddo
        ibndry(1,28)=75
        ibndry(2,44)=76
        do i=45,61
            ibndry(1,i)=i-3
            ibndry(2,i)=i-2
            ibndry(3,i)=0
            ibndry(4,i)=1
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=-3
        enddo
        ibndry(1,45)=77
        ibndry(2,61)=78
        do i=62,78
            ibndry(1,i)=i-4
            ibndry(2,i)=i-3
            ibndry(3,i)=0
            ibndry(4,i)=1
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=-4
        enddo
        ibndry(1,62)=77
        ibndry(2,78)=78
c
        do i=79,83
            ibndry(3,i)=0
            ibndry(4,i)=0
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=0
        enddo
        ibndry(1,79)=1
        ibndry(2,79)=75
        ibndry(1,80)=1
        ibndry(2,80)=77
        ibndry(1,81)=2
        ibndry(2,81)=76
        ibndry(1,82)=2
        ibndry(2,82)=78
        ibndry(1,83)=2
        ibndry(2,83)=1
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, 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-----------------------------------------------------------------------
        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(193) :: 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 ntf,nvf,nbf/2,201,205/
            data hmax,grade/0.1e0_rknd,1.75e0_rknd/
c
            data (ix(i),i=   1,  50)/
     +       14238695,  14239908,  14246652,  14251777,  14258588,
     1       14267555,  14279019,  14293652,  14311521,  14331954,
     2       14354139,  14378213,  14403770,  14430406,  14457782,
     3       14485565,  14513819,  14599728,  18376550,  18985666,
     4       19675030,  24406898,  31060792,  38332840,  45327880,
     5       51958980,  57746436,  58425076,  59000004,  59512764,
     6       59887820,  60137256,  60320872,  60353844,  60296060,
     7       60300572,  60357352,  60406648,  60466124,  60535776,
     8       60614204,  60702200,  60800248,  60907596,  61022372,
     9       61142940,  61268496,  61398644,  61533440,  61672552/
            data (ix(i),i=  51, 100)/
     +       61815780,  61962720,  62112888,  62420316,  62575812,
     1       62732588,  62890448,  63049120,  63217700,  63928776,
     2       64237548,  64175852,  64074900,  63777864,  61441936,
     3       53095484,  44016384,  35042948,  24714592,  22205500,
     4       21458752,  20755162,  20098640,  19478194,  18897466,
     5       18352342,  17845722,  16923652,  16123903,  15775547,
     6       15473247,  15211744,  14972629,  14759408,  14578352,
     7       14434788,  14338360,  14286031,  14257508,  14245168,
     8       14239976,  64350896,  64356840,  64368636,  64387924,
     9       64416380,  64450228,  64483808,  64515032,  64539980/
            data (ix(i),i= 101, 150)/
     +       64560952,  65000344,  65195092,  65355040,  65445668,
     1       65544384,  65652212,  65895780,  66510696,  66893240,
     2       67091352,  67275784,  68000344,  68109048,  68448632,
     3       68519232,  68585520,  68608048,  68600216,  68589304,
     4       68574664,  68555248,  68428072,  68324088,  68181136,
     5       67793528,  67259536,  66974700,  66684468,  66362412,
     6       65993828,  65620720,  65288752,  65006144,  64778220,
     7       64616860,  64526492,  64485632,  64467560,  64456904,
     8       64448072,  64389944,  64369788,  64354612,  66136104,
     9       66146420,  66190052,  66259372,  66350608,  66900992/
            data (ix(i),i= 151, 193)/
     +       69012104,  71521936,  73760016,  77687952,  77824840,
     1       77903744,  77885536,  77832936,  77710888,  71374592,
     2       70879848,  70390488,  70148944,  69909016,  69670448,
     3       69435512,  69207392,  68990384,  68787144,  68599888,
     4       68458072,  68323752,  68186864,  68048088,  67908096,
     5       67767704,  67627920,  67489744,  67353944,  67221368,
     6       67092772,  66968828,  66850288,  66738080,  66633016,
     7       66536588,  66449600,  66372592,  66306304,  66250944,
     8       66206916,  66177576,  66150332/
            data (iy(i),i=   1,  50)/
     +       52549636,  52486056,  52383300,  52338820,  52297704,
     1       52259352,  52223976,  52190988,  52160992,  52134800,
     2       52112484,  52092936,  52076072,  52061444,  52048608,
     3       52037264,  52027192,  52001896,  51373844,  51254296,
     4       51153928,  50666328,  50356336,  50502748,  51030500,
     5       51920772,  52987392,  53082688,  53114920,  53169252,
     6       53272384,  53386260,  53513272,  53771620,  54131468,
     7       54437204,  54693912,  54854356,  55010844,  55162632,
     8       55309492,  55450012,  55582472,  55706424,  55822392,
     9       55931420,  56033880,  56129624,  56217756,  56297756/
            data (iy(i),i=  51, 100)/
     +       56369096,  56431180,  56483644,  56562896,  56594612,
     1       56619240,  56634088,  56637748,  56633492,  56594688,
     2       56590432,  56613564,  56639912,  56704980,  57083332,
     3       58100240,  58623964,  58610828,  58002560,  57699656,
     4       57549436,  57364964,  57163400,  56941844,  56701400,
     5       56447228,  56187308,  55659188,  55137856,  54884428,
     6       54644584,  54412652,  54171612,  53918708,  53658720,
     7       53393196,  53158128,  52970972,  52820080,  52710904,
     8       52623296,  55603292,  55530572,  55459972,  55391252,
     9       55330720,  55284948,  55251264,  55225416,  55209376/
            data (iy(i),i= 101, 150)/
     +       55199980,  55015724,  54946692,  54898108,  54873704,
     1       54851376,  54833364,  54804220,  54695308,  54609280,
     2       54556216,  54498012,  54238536,  54183996,  53959824,
     3       53916424,  53878688,  53896588,  53920760,  53945840,
     4       53972352,  54002572,  54174428,  54303320,  54464948,
     5       54862796,  55358280,  55597788,  55816220,  56015336,
     6       56181592,  56289276,  56327604,  56308872,  56250920,
     7       56174364,  56106268,  56061620,  56034172,  56012820,
     8       55992336,  55821324,  55753340,  55678208,  51671528,
     9       51561080,  51450264,  51342880,  51241684,  50778044/
            data (iy(i),i= 151, 193)/
     +       48983984,  46596432,  44191344,  39426168,  39296392,
     1       39230424,  39332364,  39514296,  39755412,  49894108,
     2       50545140,  51122660,  51380640,  51613100,  51818984,
     3       51994880,  52144128,  52269792,  52373748,  52457108,
     4       52512252,  52557100,  52594560,  52624112,  52644932,
     5       52655976,  52656652,  52647096,  52629192,  52603440,
     6       52570384,  52530612,  52484568,  52432032,  52373524,
     7       52308376,  52237708,  52162856,  52085172,  52006220,
     8       51927568,  51860404,  51772272/
c
        sp(2)='wing'
        sp(1)='wing'
        sp(3)='wing'
c
        pi=3.141592653589793e0_rknd
        nw=9
        nf1=100
        nf2=153
c
        size=ru(7)
        do i=1,8
            arg=pi*real(i-1,rknd)/4.0e0_rknd
            vx(i)=size*cos(arg)
            vy(i)=size*sin(arg)
        enddo
        xmin=real(ix(1),rknd)
        xmax=real(ix(1),rknd)
        ymin=real(iy(1),rknd)
        ymax=real(iy(1),rknd)
        do i=1,nvf-8
           xmax=max(xmax,real(ix(i),rknd))
           xmin=min(xmin,real(ix(i),rknd))
           ymax=max(ymax,real(iy(i),rknd))
           ymin=min(ymin,real(iy(i),rknd))
        enddo
        xx=(xmin+xmax)/2.0e0_rknd
        yy=(ymin+ymax)/2.0e0_rknd
        ss=max(xmax-xmin,ymax-ymin)
        do i=9,nvf
            vx(i)=(real(ix(i-8),rknd)-xx)/ss
            vy(i)=(real(iy(i-8),rknd)-yy)/ss
        enddo
c
        do i=1,nbf
            ibndry(1,i)=i
            ibndry(2,i)=i+1
            ibndry(3,i)=0
            if(i<=nvf) then
                ibndry(4,i)=1
            else
                ibndry(4,i)=0
            endif
            ibndry(5,i)=0
            ibndry(6,i)=0
            if(i<nw) then
                ibndry(7,i)=i
            else
                ibndry(7,i)=-1
            endif
        enddo
        ibndry(2,nw-1)=1
        ibndry(2,nf1-1)=nw
        ibndry(2,nf2-1)=nf1
        ibndry(2,nvf)=nf2
c
        ibndry(1,nvf+1)=5
        ibndry(2,nvf+1)=10
        ibndry(1,nvf+2)=69
        ibndry(2,nvf+2)=144
        ibndry(1,nvf+3)=125
        ibndry(2,nvf+3)=180
        ibndry(1,nvf+4)=164
        ibndry(2,nvf+4)=1
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        rp(12)=hmax
        rp(13)=grade
c
        do i=1,nbf
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
c
c       make itnode
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 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(232) :: 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 ntf,nvf,nbf/4,240,244/
            data hmax,grade/0.1e0_rknd,1.75e0_rknd/
c
c       wing
c
            data (ix(i),i=   1,  50)/
     +       -7341201,  -7525801,  -7734500,  -7897900,  -7964801,
     1       -7751900,  -7210900,  -5903800,  -4446800,  -2864700,
     2       -1615700,  -1197700,  -1272700,  -2115900,  -2953400,
     3       -3785300,  -4610400,  -5428500,  -6237800,  -7035100,
     4       -7814601,  -8565500,  -9265900,  -9584401,  -9867100,
     5       -9988501, -10090600, -10122300, -10147800, -10157701,
     6      -10161901, -10162600, -10150400, -10126700, -10098301,
     7      -10000000,  -9803300,  -9714200,  -9577101,  -9466301,
     8       -9365001,  -9272201,  -9185200,  -9023101,  -8873101,
     9       -8519600,  -8184500,  -7552900,  -6960800,  48000000/
            data (ix(i),i=  51, 100)/
     +       46999500,  45999900,  45500100,  44499900,  43499800,
     1       42499800,  41500000,  40500300,  39500000,  39000100,
     2       38499800,  38000100,  37500000,  37000000,  36500000,
     3       36000000,  35000000,  32999900,  30000000,  28000000,
     4       26000000,  24000000,  20000000,  16000000,  12000000,
     5       11000000,  10000000,   9000000,   8000000,   7000000,
     6        6000100,   5000100,   4000000,   2999800,   1999600,
     7        1500000,   1000400,    750500,    500900,    399800,
     8         299500,    250400,    199600,    150700,     98600,
     9          49900,     24300,     -3900,     22900,     49700/
            data (ix(i),i= 101, 150)/
     +         101000,    149700,    200000,    250100,    300100,
     1         400100,    498500,    748600,   1001300,   1504300,
     2        2003500,   2998700,   3998200,   4999400,   6000100,
     3        7000100,   8000000,   9000000,  10000000,  11000000,
     4       12000000,  16000000,  19999700,  23998900,  25999000,
     5       28000700,  29001402,  30001300,  30998100,  31498500,
     6       32000300,  32250000,  32500202,  32749600,  32999600,
     7       33250000,  33500102,  33625300,  33749400,  33875500,
     8       33999800,  34125200,  34248800,  34499100,  34751200,
     9       35000700,  35250600,  35497400,  35749900,  36001200/
            data (ix(i),i= 151, 200)/
     +       36250700,  36501700,  36995400,  37503300,  37997700,
     1       38500800,  39001800,  39501200,  40499800,  41500200,
     2       42500200,  43500000,  44499100,  45499800,  46000400,
     3       47000200,  48000000,  58745600,  58064504,  57368904,
     4       56657300,  55920400,  55160000,  54770204,  54374300,
     5       53972900,  53565400,  53152600,  52734400,  52311500,
     6       51883900,  51450000,  51013400,  50576300,  50137200,
     7       49695300,  49252000,  48809100,  48366200,  48014000,
     8       47750700,  47579500,  47412100,  47338000,  47288400,
     9       47258000,  47267200,  47352200,  47498500,  47593400/
            data (ix(i),i= 201, 232)/
     +       47754800,  47902100,  48040400,  48299000,  48542200,
     1       48891600,  49335700,  49847100,  50335800,  50799700,
     2       51244800,  51676500,  52094704,  52501500,  52897100,
     3       53282500,  53434004,  53657900,  53879104,  54024404,
     4       54240400,  54383300,  54595696,  54736604,  54945204,
     5       55082704,  55286904,  55421004,  55621300,  55752804,
     6       55949700,  56079900/
            data (iy(i),i=   1,  50)/
     +       -8805301,  -8565400,  -8224000,  -7741001,  -6625000,
     1       -5693800,  -4630900,  -2894700,  -1418400,   -158600,
     2         678000,    954000,   1083900,    544400,     -5000,
     3        -564200,  -1135100,  -1718100,  -2316300,  -2935300,
     4       -3585300,  -4284600,  -5071500,  -5519900,  -6030100,
     5       -6319900,  -6643100,  -6788100,  -6944000,  -7026800,
     6       -7119600,  -7218400,  -7339400,  -7480500,  -7579700,
     7       -7800000,  -8090600,  -8195000,  -8332400,  -8424300,
     8       -8499801,  -8560501,  -8611300,  -8692100,  -8751900,
     9       -8864200,  -8944600,  -9038601,  -9064000,   -562900/
            data (iy(i),i=  51, 100)/
     +        -340200,   -114800,     -1000,    228100,    456300,
     1         679600,    896000,   1104900,   1305600,   1402000,
     2        1495500,   1585700,   1673200,   1758300,   1841000,
     3        1921100,   2073700,   2349000,   2693900,   2883300,
     4        3043000,   3174400,   3353000,   3411500,   3329900,
     5        3284300,   3227400,   3158700,   3077100,   2981000,
     6        2867300,   2730000,   2557600,   2328000,   1996800,
     7        1766900,   1465700,   1275500,   1046400,    937700,
     8         815600,    748200,    670400,    584600,    474700,
     9         340000,    240100,         0,   -350300,   -485100/
            data (iy(i),i= 101, 150)/
     +        -672200,   -807600,   -923500,  -1022500,  -1109900,
     1       -1260900,  -1389100,  -1663300,  -1897500,  -2290100,
     2       -2608100,  -3082600,  -3412700,  -3656400,  -3849300,
     3       -4008600,  -4141900,  -4253000,  -4345100,  -4420300,
     4       -4480100,  -4580000,  -4462600,  -4128600,  -3890500,
     5       -3614000,  -3457800,  -3282100,  -3087700,  -2992100,
     6       -2906800,  -2867400,  -2826300,  -2778200,  -2718300,
     7       -2645100,  -2560200,  -2513400,  -2463500,  -2408500,
     8       -2349500,  -2285300,  -2217600,  -2068500,  -1906800,
     9       -1741000,  -1570800,  -1396400,  -1209900,  -1021500/
            data (iy(i),i= 151, 200)/
     +        -841000,   -672700,   -385100,   -161400,    -20300,
     1          87700,    182200,    257800,    336500,    329800,
     2         261500,    145000,    -14900,   -211200,   -318200,
     3        -539300,   -763000, -10754700, -10010300,  -9283500,
     4       -8576001,  -7898300,  -7249000,  -6935400,  -6629600,
     5       -6330300,  -6037600,  -5752100,  -5472700,  -5198700,
     6       -4930600,  -4669400,  -4412300,  -4155500,  -3900900,
     7       -3649800,  -3400200,  -3150300,  -2900300,  -2697700,
     8       -2544900,  -2437800,  -2326200,  -2258900,  -2162500,
     9       -2043100,  -1954200,  -1775000,  -1678300,  -1643000/
            data (iy(i),i= 201, 232)/
     +       -1606100,  -1586100,  -1576700,  -1579500,  -1600700,
     1       -1650900,  -1743900,  -1912100,  -2107400,  -2332300,
     2       -2579700,  -2843200,  -3122600,  -3415700,  -3722100,
     3       -4040700,  -4171200,  -4370500,  -4574300,  -4711700,
     4       -4921200,  -5062100,  -5275400,  -5419100,  -5637100,
     5       -5784200,  -6008400,  -6159200,  -6388300,  -6541600,
     6       -6774300,  -6929800/
c
        sp(2)='wing'
        sp(1)='wing'
        sp(3)='wing'
c
        pi=3.141592653589793e0_rknd
c
        size=ru(7)
        do i=1,8
            arg=pi*real(i-1,rknd)/4.0e0_rknd
            vx(i)=size*cos(arg)
            vy(i)=size*sin(arg)
        enddo
        xmin=real(ix(1),rknd)
        xmax=real(ix(1),rknd)
        ymin=real(iy(1),rknd)
        ymax=real(iy(1),rknd)
        do i=1,nvf-8
           xmax=max(xmax,real(ix(i),rknd))
           xmin=min(xmin,real(ix(i),rknd))
           ymax=max(ymax,real(iy(i),rknd))
           ymin=min(ymin,real(iy(i),rknd))
        enddo
        xx=(xmin+xmax)/2.0e0_rknd
        yy=(ymin+ymax)/2.0e0_rknd
        ss=max(xmax-xmin,ymax-ymin)
        do i=9,nvf
            vx(i)=(real(ix(i-8),rknd)-xx)/ss
            vy(i)=(real(iy(i-8),rknd)-yy)/ss
        enddo
c
        do i=1,nbf
            ibndry(1,i)=i
            ibndry(2,i)=i+1
            ibndry(3,i)=0
            if(i<=nvf) then
                ibndry(4,i)=1
            else
                ibndry(4,i)=0
            endif
            ibndry(5,i)=0
            ibndry(6,i)=0
            if(i<=8) then
                ibndry(7,i)=i
            else
                ibndry(7,i)=-1
            endif
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
        ns=9
        nw=ns+49
        nf=nw+118
        nn=nf+73
        ibndry(2,ns-1)=1
        ibndry(2,nw-1)=ns
        ibndry(2,nf-1)=nw
        ibndry(2,nvf)=nf
        imn=nw
        imx=nw
        do i=nw,nf-1
            if(vx(i)<=vx(imn)) imn=i
            if(vx(i)>=vx(imx)) imx=i
        enddo
        jmn=ns
        jmx=ns
        do i=ns,nw-1
           qj=(vx(5)-vx(jmn))**2+(vy(5)-vy(jmn))**2
           qi=(vx(5)-vx(i))**2+(vy(5)-vy(i))**2
           if(qi<=qj) jmn=i
           qj=(vx(imn)-vx(jmx))**2+(vy(imn)-vy(jmx))**2
           qi=(vx(imn)-vx(i))**2+(vy(imn)-vy(i))**2
           if(qi<=qj) jmx=i
        enddo
        kmn=nf
        kmx=nf
        do i=nf,nn-1
           qj=(vx(imx)-vx(kmn))**2+(vy(imx)-vy(kmn))**2
           qi=(vx(imx)-vx(i))**2+(vy(imx)-vy(i))**2
           if(qi<=qj) kmn=i
           qj=(vx(1)-vx(kmx))**2+(vy(1)-vy(kmx))**2
           qi=(vx(1)-vx(i))**2+(vy(1)-vy(i))**2
           if(qi<=qj) kmx=i
        enddo
        imn=nw
        imx=nw
        do i=nw,nf-1
           qj=(vx(jmx)-vx(imn))**2+(vy(jmx)-vy(imn))**2
           qi=(vx(jmx)-vx(i))**2+(vy(jmx)-vy(i))**2
           if(qi<=qj) imn=i
           qj=(vx(kmn)-vx(imx))**2+(vy(kmn)-vy(imx))**2
           qi=(vx(kmn)-vx(i))**2+(vy(kmn)-vy(i))**2
           if(qi<=qj) imx=i
        enddo
c
        ibndry(1,nvf+1)=5
        ibndry(2,nvf+1)=jmn
        ibndry(1,nvf+2)=jmx
        ibndry(2,nvf+2)=imn
        ibndry(1,nvf+3)=imx
        ibndry(2,nvf+3)=kmn
        ibndry(1,nvf+4)=kmx
        ibndry(2,nvf+4)=1
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)
c
        return
        end