c***********************  problem name: jcn     ************************
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
            common /val0/k0,ku,kx,ky,kl
cy
        call cbeta(x,y,betax,betay,itag)
        values(k0)=ux+u*betax
        values(ku)=betax
        values(kx)=1.0e0_rknd
        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
            common /val0/k0,ku,kx,ky,kl
cy
        call cbeta(x,y,betax,betay,itag)
        values(k0)=uy+u*betay
        values(ku)=betay
        values(ky)=1.0e0_rknd
        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
            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
            character(len=80) :: su
            common /val2/k0,kl,klb,kub,kic,kim,kil
            common /atest2/ival,iob,iu(98),top,bottom,du,
     +          angle,ct,st,xt(7),yt(7),ru(80),su(100)
cy
        values(k0)=top
        if(itag==1) values(k0)=bottom
        values(kic)=values(k0)
        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
        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
        za=0.5e0_rknd*abs(u)
        if(za>1.0e5_rknd) then
            as=2.0e0_rknd*za
        else
            as=za+sqrt(1.0e0_rknd+za*za)
        endif
c
        if(u<0.0e0_rknd) then
             v=-log10(as)
        else
             v=log10(as)
        endif
c
        values(kf)=v
        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 cbeta(x,y,betax,betay,itag)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            character(len=80) :: su
            common /atest2/ival,iob,iu(98),top,bottom,du,
     +          angle,ct,st,xt(7),yt(7),ru(80),su(100)
cy
        betax=du*(ct*xt(itag)+st*yt(itag))
        betay=du*(ct*yt(itag)-st*xt(itag))
        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
        aa=ru(4)
c
        call usrset(iu,ru,su)
c
        if(iu(1)/=0) ip(41)=-1
        if(ru(4)/=aa) then
            pi=3.141592653589793e0_rknd
            theta=ru(4)*pi/180.0e0_rknd
            ru(5)=cos(theta)
            ru(6)=sin(theta)
        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,  7)/
     +  'n i= 1,n=obtuse,a= o,t=i',
     1  'n i= 1,n=top   ,a= t,t=r',
     2  'n i= 2,n=bottom,a= b,t=r',
     3  'n i= 3,n=du    ,a= d,t=r',
     4  'n i= 4,n=angle ,a= a,t=r',
     5  's n=obtuse,v=0,l="do not trisect elements"',
     6  's n=obtuse,v=1,l="trisect all elements"'/
c
            data len0/7/
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, dimension(23) :: b1,b2,b4
            integer(kind=iknd), save, dimension(7) :: ib1,ib2
            integer(kind=iknd), save, dimension(11) :: ix,iy
            integer(kind=iknd), save :: ntf,nvf,nbf,ispd,iprob
            integer(kind=iknd), save :: ifirst
            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
c
            data b1/ 1, 2, 3, 4, 1, 6, 7, 6, 8, 9,
     +               8,10,11,12,13,14,15,16,17, 7,
     1               12,13,14/
            data b2/ 2, 3, 4, 5, 6, 7,12, 8, 9,15,
     +              10,11, 5,13,14, 3,16,17, 4, 9,
     1              15,16,17/
            data b4/ 2, 1, 1, 1, 1, 0, 0, 1, 0, 0,
     +               1, 2, 1, 0, 0, 0, 0, 0, 0, 0,
     1               0, 0, 0/
c
            data ib1/ 8,6, 7,12,14,14,3/
            data ib2/11,8,20,21,15,23,2/
c
            data ix/175,200,225,250,275,175,200,175,200,175,275/
            data iy/0,0,0,0,0,-25,-25,-50,-50,-100,-100/
            data ntf,nvf,nbf,ispd,iprob,ifirst/7,17,23,0,1,1/
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
c
c
        if(ip(41)==1) then
            sp(1)='continuity equation'
            sp(2)='continuity equation'
            sp(3)='continuity equation'
            sp(4)='jcn'
            sp(5)='jcn_figxxx.ext'
            sp(6)='jcn_mpixxx.rw'
            sp(7)='jcn.jnl'
            sp(9)='jcn_mpixxx.out'
c
            iu(1)=0
            ru(1)=1.0e10_rknd
            ru(2)=1.0e-5_rknd
            ru(3)=(40.0e0_rknd+15.0e0_rknd*log(10.0e0_rknd))
            ru(4)=0.0e0_rknd
            ru(5)=1.0e0_rknd
            ru(6)=0.0e0_rknd
        endif
c
        if(iu(1)/=0) go to 100
c
c       set up regions
c
        ip(5)=max(ip(5),ifirst)
        ip(6)=iprob
        ip(8)=ispd
        ip(19)=1
        ip(20)=5
        rp(12)=hmax
        rp(13)=grade
c
        ss=2.99563e-2_rknd
        pi=3.141592653589793e0_rknd
c
        s1=ss*1.0e-2_rknd
        do i=1,11
            vx(i)=real(ix(i),rknd)*s1
            vy(i)=real(iy(i),rknd)*s1
        enddo
        do i=1,3
            arg=(12.0e0_rknd+real(i,rknd))*pi/8.0e0_rknd
            c=cos(arg)*ss
            s=sin(arg)*ss
            vx(11+i)=vx(2)+0.25e0_rknd*c
            vy(11+i)=0.25e0_rknd*s
            vx(14+i)=vx(2)+0.5e0_rknd*c
            vy(14+i)=0.5e0_rknd*s
        enddo
c
c       constants for a1xy and a2xy
c
        ru(7)=0.0e0_rknd
        ru(14)=0.0e0_rknd
c
        ru(8)=0.0e0_rknd
        ru(15)=1.0e0_rknd/(vy(9)-vy(7))
c
        xx=(vx(9)-vx(7)+vx(15)-vx(12))/2.0e0_rknd
        yy=(vy(9)-vy(7)+vy(15)-vy(12))/2.0e0_rknd
        dd=1.0e0_rknd/(xx*xx+yy*yy)
        ru(9)=xx*dd
        ru(16)=yy*dd
c
        xx=(vx(15)-vx(12)+vx(16)-vx(13))/2.0e0_rknd
        yy=(vy(15)-vy(12)+vy(16)-vy(13))/2.0e0_rknd
        dd=1.0e0_rknd/(xx*xx+yy*yy)
        ru(10)=xx*dd
        ru(17)=yy*dd
c
        xx=(vx(16)-vx(13)+vx(17)-vx(14))/2.0e0_rknd
        yy=(vy(16)-vy(13)+vy(17)-vy(14))/2.0e0_rknd
        dd=1.0e0_rknd/(xx*xx+yy*yy)
        ru(11)=xx*dd
        ru(18)=yy*dd
c
        xx=(vx(17)-vx(14)+vx(4)-vx(3))/2.0e0_rknd
        yy=(vy(17)-vy(14)+vy(4)-vy(3))/2.0e0_rknd
        dd=1.0e0_rknd/(xx*xx+yy*yy)
        ru(12)=xx*dd
        ru(19)=yy*dd
c
        ru(13)=0.0e0_rknd
        ru(20)=0.0e0_rknd
c
        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
        ibndry(7,12)=1
        ibndry(7,1)=3
c
        do i=1,ntf
            itnode(1,i)=ib1(i)
            itnode(2,i)=ib2(i)
            itnode(3,i)=0
            itnode(4,i)=0
            itnode(5,i)=i
        enddo
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        return
c
c       make obtuse angles
c
  100   iu(1)=0
        ntf=ip(1)
        nt0=ntf
        nvf=ip(2)
        nbf=ip(3)
        maxt=ip(83)
        maxv=ip(84)
        if (ntf*3>maxt) return
        if (nvf+ntf>maxv) return
        do it=1,nt0
            iv1=itnode(1,it)
            iv2=itnode(2,it)
            iv3=itnode(3,it)
            iv4=itnode(4,it)
            iv5=itnode(5,it)
c
            nvf=nvf+1
            vx(nvf)=(vx(iv1)+vx(iv2)+vx(iv3))/3.0e0_rknd
            vy(nvf)=(vy(iv1)+vy(iv2)+vy(iv3))/3.0e0_rknd
c
            itnode(1,ntf+1)=nvf
            itnode(2,ntf+1)=iv2
            itnode(3,ntf+1)=iv3
            itnode(4,ntf+1)=iv4
            itnode(5,ntf+1)=iv5
            itnode(1,ntf+2)=iv1
            itnode(2,ntf+2)=nvf
            itnode(3,ntf+2)=iv3
            itnode(4,ntf+2)=iv4
            itnode(5,ntf+2)=iv5
            itnode(3,it)=nvf
            ntf=ntf+2
        enddo
        ip(1)=ntf
        ip(2)=nvf
        return
        end