c***********************  problem name: circle  ************************
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 11.0 - - - june, 2012
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 /atest2/iu(100),ru(100),su(100)
            common /val0/k0,ku,kx,ky,kl,kuu,kux,kuy,kul,
     +          kxu,kxx,kxy,kxl,kyu,kyx,kyy,kyl,klu,klx,kly,kll
cy
        values(k0)=ux*ru(itag)
        values(kx)=ru(itag)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 11.0 - - - june, 2012
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 /atest2/iu(100),ru(100),su(100)
            common /val0/k0,ku,kx,ky,kl,kuu,kux,kuy,kul,
     +          kxu,kxx,kxy,kxl,kyu,kyx,kyy,kyl,klu,klx,kly,kll
cy
        values(k0)=uy*ru(itag)
        values(ky)=ru(itag)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 11.0 - - - june, 2012
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,kuu,kux,kuy,kul,
     +          kxu,kxx,kxy,kxl,kyu,kyx,kyy,kyl,klu,klx,kly,kll
cy
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 11.0 - - - june, 2012
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,kuu,kul,klu,kll
            common /atest2/iu(100),ru(100),su(100)
cy
        if(itag>iu(1)) return
c
c       (x,y) is unit outward normal
c
        if(itag<1) return
        call uexact(x,y,itag,r,rx,ry,rxx,ryy,rxy)
        if(itag<=iu(1)) then
            values(k0)=(x*rx+y*ry)*ru(itag)
        else
            values(k0)=ry*ru(itag-1)
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 11.0 - - - june, 2012
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,kll,klb,kub,kic,kim,kil
            common /atest2/iu(100),ru(100),su(100)
cy
        if(itag<1.or.itag>iu(1)) return
        call uexact(x,y,itag,r,rx,ry,rxx,ryy,rxy)
        values(k0)=r
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 11.0 - - - june, 2012
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,kuu,kux,kuy,kul,
     +          kxu,kxx,kxy,kxl,kyu,kyx,kyy,kyl,klu,klx,kly,kll
cy
        xx=0.5e0_rknd
        yy=0.5e0_rknd
        e=1.0e4_rknd*exp(-5.0e0_rknd*((x-xx)**2+(y-yy)**2))
        values(k0)=e*u**2
        values(ku)=e*2.0e0_rknd*u
        values(kuu)=e*2.0e0_rknd
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 11.0 - - - june, 2012
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,kuu,kux,kuy,kul,
     +          kxu,kxx,kxy,kxl,kyu,kyx,kyy,kyl,klu,klx,kly,kll
cy
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 11.0 - - - june, 2012
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),ru(100),su(100)
cy
        call uexact(x,y,itag,r,rx,ry,rxx,ryy,rxy)
        values(kf)=r-u
        values(kf1)=rx-ux
        values(kf2)=ry-uy
        values(kad)=r
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 11.0 - - - june, 2012
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(*) :: values
            common /val4/jx,jy,jxs,jys,jxl,jyl,jxss,jyss,jxll,jyll,
     +          jxsl,jysl,jxls,jyls
cy
        values(jx)=cos(s)
        values(jy)=sin(s)
        values(jxs)=-values(jy)
        values(jys)=values(jx)
        values(jxss)=-values(jx)
        values(jyss)=-values(jy)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 11.0 - - - june, 2012
c
c-----------------------------------------------------------------------
        subroutine uexact(x,y,itag,u,ux,uy,uxx,uyy,uxy)
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),ru(100),su(100)
cy
        u=0.0e0_rknd
        ux=0.0e0_rknd
        uy=0.0e0_rknd
        uxx=0.0e0_rknd
        uxy=0.0e0_rknd
        uyy=0.0e0_rknd
        r=sqrt(x**2+y**2)
        if(r<=0.0e0_rknd) return
c
        pi=3.141592653589793e0_rknd
        al=ru(25)
        arg=min(x/r,1.0e0_rknd)
        arg=max(-1.0e0_rknd,arg)
        theta=acos(arg)
        if(itag>=5) theta=2.0e0_rknd*pi-theta
c
        s=r**al*sin(theta*al)
        al1=al-1.0e0_rknd
        al2=al-2.0e0_rknd
        sx=al*r**al1*sin(theta*al1)
        sy=al*r**al1*cos(theta*al1)
        sxx=al*al1*r**al2*sin(theta*al2)
        syy=-sxx
        sxy=al*al1*r**al2*cos(theta*al2)
c
        c=r**al*cos(theta*al)
        cx=sy
        cy=-sx
        cxx=sxy
        cyy=-sxy
        cxy=syy
c
        cf=ru(itag+8)
        sf=ru(itag+16)
        u=cf*c+sf*s
        ux=cf*cx+sf*sx
        uy=cf*cy+sf*sy
        uxx=cf*cxx+sf*sxx
        uyy=cf*cyy+sf*syy
        uxy=cf*cxy+sf*sxy
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 11.0 - - - june, 2012
c
c-----------------------------------------------------------------------
        subroutine cnu(ntri,af,cf,sf,vv)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            real(kind=rknd), dimension(8) :: cf,sf,af,cd,sd
            real(kind=rknd), dimension(1000) :: val
cy
        eps=1.0e1_rknd*epsilon(1.0e0_rknd)
        cf(1)=0.0e0_rknd
        sf(1)=1.0e0_rknd
        cd(1)=0.0e0_rknd
        sd(1)=0.0e0_rknd
        vmax=2.0e0_rknd
        n=1000
        do i=1,n
            vv=real(i,rknd)*vmax/real(n,rknd)
            call veval(ntri,cf,sf,af,vv,val(i),dp)
        enddo
        do i=2,n
           if(val(i-1)>0.0e0_rknd.and.val(i)<=0.0e0_rknd) then
                vmin=real(i-1,rknd)*vmax/real(n,rknd)
                vmax=real(i,rknd)*vmax/real(n,rknd)
                go to 10
            endif
        enddo
        vv=vmax
        return
c
   10   itmax=100
        vv=(vmin+vmax)/2.0e0_rknd
        do k=1,itmax
            call veval(ntri,cf,sf,af,vv,dd,dp)
            if(dd>0.0e0_rknd) then
                vmin=vv
            else
                vmax=vv
            endif
            vv=vv-dd/dp
            if(vv<=vmin.or.vv>=vmax)
     +          vv=(vmin+vmax)/2.0e0_rknd
            if(abs(dd)<=eps) return
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 11.0 - - - june, 2012
c
c-----------------------------------------------------------------------
        subroutine veval(ntri,cf,sf,af,vv,dd,dp)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            real(kind=rknd), dimension(8) :: cf,sf,af,cd,sd
cy
        pi=3.141592653589793e0_rknd
        cf(1)=0.0e0_rknd
        sf(1)=1.0e0_rknd
        cd(1)=0.0e0_rknd
        sd(1)=0.0e0_rknd
c
        do i=1,ntri-1
            if(af(i+1)==af(i)) then
                cf(i+1)=cf(i)
                sf(i+1)=sf(i)
                cd(i+1)=cd(i)
                sd(i+1)=sd(i)
            else
                theta=real(i,rknd)*pi/4.0e0_rknd
                c=cos(theta*vv)
                s=sin(theta*vv)
                dc=-theta*s
                ds=theta*c
c
                ff=c*cf(i)+s*sf(i)
                dd=-s*cf(i)+c*sf(i)
                df=c*cd(i)+s*sd(i)+dc*cf(i)+ds*sf(i)
                dp=-s*cd(i)+c*sd(i)-ds*cf(i)+dc*sf(i)
c
                dd=dd*af(i)/af(i+1)
                dp=dp*af(i)/af(i+1)
c
                cf(i+1)=c*ff-s*dd
                sf(i+1)=s*ff+c*dd
                cd(i+1)=dc*ff-ds*dd+c*df-s*dp
                sd(i+1)=ds*ff+dc*dd+s*df+c*dp
            endif
        enddo
        theta=real(ntri,rknd)*pi/4.0e0_rknd
        c=cos(theta*vv)
        s=sin(theta*vv)
        dc=-theta*s
        ds=theta*c
c
        dd=-s*cf(ntri)+c*sf(ntri)
        dp=-ds*cf(ntri)+dc*sf(ntri)-s*cd(ntri)+c*sd(ntri)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 11.0 - - - june, 2012
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
            integer(kind=iknd), save :: len
            real(kind=rknd), dimension(*) :: vx,vy
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(100) :: rp,ru
            real(kind=rknd), dimension(8) :: psave
            character(len=80), dimension(100) :: sp,su
            character(len=80), save, dimension(12) :: file
cy
            data len/12/
            data (file(i),i=1,10)/
     +      'n i=1,n=a1,a=a1,t=r',
     1      'n i=2,n=a2,a=a2,t=r',
     2      'n i=3,n=a3,a=a3,t=r',
     3      'n i=4,n=a4,a=a4,t=r',
     4      'n i=5,n=a5,a=a5,t=r',
     5      'n i=6,n=a6,a=a6,t=r',
     6      'n i=7,n=a7,a=a7,t=r',
     7      'n i=8,n=a8,a=a8,t=r',
     8      'n i=1,n=ntri,a=nt,t=i',
     9      'n i=2,n=ibc,a=bc,t=i'/
            data (file(i),i=11,12)/
     +      's n=ibc,v=1,l="neumann"',
     1      's n=ibc,v=2,l="dirichlet"'/
c
        ntf=iu(1)
        ibc=iu(2)
        do i=1,8
            psave(i)=ru(i)
        enddo
c
        call usrset(file,len,iu,ru,su)
c
        iu(1)=min(8,iu(1))
        iu(1)=max(1,iu(1))
        if(iu(2)/=1) iu(2)=2
        do i=1,8
            if(ru(i)<=0.0e0_rknd) ru(i)=1.0e0_rknd
        enddo
c
        if(ntf/=iu(1).or.ibc/=iu(2)) ip(41)=-1
        do i=1,8
            if(psave(i)/=ru(i)) ip(41)=-1
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 11.0 - - - june, 2012
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=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
c       unit circle divided into 8 equal parts.
c       this problem can be run with ntf=1,2,...8,
c       nvf=ntf+2,nbf=ntf+2
c
        if(ip(41)==1) then
            sp(1)='circle'
            sp(2)='circle'
            sp(3)='circle'
            sp(6)='circle_mpixxx.rw'
            sp(7)='circle.jnl'
            sp(9)='circle_mpixxx.out'
c
            iu(1)=8
            iu(2)=2
            do i=1,8
                ru(i)=1.0e0_rknd
            enddo
        endif
c
        ntf=iu(1)
        ibc=iu(2)
        nvf=ntf+2
        nbf=ntf+2
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        ip(5)=max(ip(5),ifirst)
        ip(22)=100000
        ip(6)=iprob
        ip(8)=ispd
        ip(7)=itask
c
        pi=3.141592653589793e0_rknd
        do i=1,ntf
            itnode(1,i)=1
            itnode(2,i)=i+1
            itnode(3,i)=i+2
            itnode(4,i)=0
            itnode(5,i)=i
        enddo
        vx(1)=0.0e0_rknd
        vy(1)=0.0e0_rknd
        do i=2,nvf
            arg=pi*real(i-2,rknd)/4.0e0_rknd
            vx(i)=cos(arg)
            vy(i)=sin(arg)
        enddo
        do i=1,nbf
            ibndry(1,i)=i
            ibndry(2,i)=i+1
            ibndry(3,i)=1
            ibndry(4,i)=ibc
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=i-1
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
        ibndry(2,nbf)=1
        ibndry(3,1)=0
        ibndry(3,nbf)=0
        ibndry(4,1)=2
        ibndry(4,nbf)=1
c
        pi4=pi/4.0e0_rknd
        do i=2,nbf-1
            sf(1,i)=real(i-2,rknd)*pi4
            sf(2,i)=real(i-1,rknd)*pi4
            ibndry(3,i)=-i
        enddo
c
c
        call cnu(ntf,ru(1),ru(9),ru(17),ru(25))
c
        sp(1)='alpha = '
        call sreal(sp(1)(9:9),nn,ru(25),3,1)
        sp(2)=sp(1)
        sp(3)=sp(1)
        return
        end