c***********************  problem name: square  ************************
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 12.0 - - - june, 2016
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/ibc(4),icont,iu(95),a1x,a1y,a1u,a2x,a2y,a2u,
     +          bux,buy,cu0,cu1,cu2,cahn,cexp,cir,csin,du0,du1,eu0,f0,
     1          ru(81),su(100)
cy
        call setrl(rl)
        values(k0)=a1x*ux+a1y*uy+a1u*u
        values(ku)=a1u
        values(kx)=a1x
        values(ky)=a1y
        if(icont==1) then
            values(kl)=ux
        else if(icont==2) then
            values(kl)=uy
        else if(icont==3) then
            values(kl)=u
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 12.0 - - - june, 2016
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/ibc(4),icont,iu(95),a1x,a1y,a1u,a2x,a2y,a2u,
     +          bux,buy,cu0,cu1,cu2,cahn,cexp,cir,csin,du0,du1,eu0,f0,
     1          ru(81),su(100)
cy
        call setrl(rl)
        values(k0)=a2x*ux+a2y*uy+a2u*u
        values(ku)=a2u
        values(kx)=a2x
        values(ky)=a2y
        if(icont==4) then
            values(kl)=ux
        else if(icont==5) then
            values(kl)=uy
        else if(icont==6) then
            values(kl)=u
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 12.0 - - - june, 2016
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/ibc(4),icont,iu(95),a1x,a1y,a1u,a2x,a2y,a2u,
     +          bux,buy,cu0,cu1,cu2,cahn,cexp,cir,csin,du0,du1,eu0,f0,
     1          ru(81),su(100)
cy
        call setrl(rl)
        values(k0)= - bux*ux - buy*uy
     +      - cir*(ux*(y-0.5e0_rknd)-uy*(x-0.5e0_rknd))
     1      - cu0 - cu1*u - cu2*u**2
     2      - cahn*u*(1.0e0_rknd-u**2) -f0*(y-x)
        values(ku)= - cu1 - cu2*u*2.0e0_rknd
     +      - cahn*(1.0e0_rknd-3.0e0_rknd*u**2)
        values(kx) = - bux - cir*(y-0.5e0_rknd)
        values(ky) = - buy + cir*(x-0.5e0_rknd)
c
        if(cexp/=0.0e0_rknd.or.icont==13) then
            ee=exp(u)
            values(k0)=values(k0) - cexp*ee
            values(ku)=values(ku) - cexp*ee
        endif
        if(csin/=0.0e0_rknd.or.icont==15) then
            ss=sin(u)
            cc=cos(u)
            values(k0)=values(k0) - csin*ss
            values(ku)=values(ku) - csin*cc
        endif
c
        if(icont==7) then
            values(kl) = - ux
        else if(icont==8) then
            values(kl) = - uy
        else if(icont==9) then
            values(kl) = - 1.0e0_rknd
        else if(icont==10) then
            values(kl) = - u
        else if(icont==11) then
            values(kl) = - u**2
        else if(icont==12) then
            values(kl) = - u*(1.0e0_rknd-u**2)
        else if(icont==13) then
            values(kl) = - ee
        else if(icont==14) then
            values(kl) = - (ux*(y-0.5e0_rknd)-uy*(x-0.5e0_rknd))
        else if(icont==15) then
            values(kl) = - ss
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 12.0 - - - june, 2016
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/ibc(4),icont,iu(95),a1x,a1y,a1u,a2x,a2y,a2u,
     +          bux,buy,cu0,cu1,cu2,cahn,cexp,cir,csin,du0,du1,eu0,f0,
     1          ru(81),su(100)
cy
        call setrl(rl)
        values(k0)= - du0 - du1*u
        values(ku)= - du1
        if(icont==16) then
            values(kl)=-1.0e0_rknd
        else if(icont==17) then
            values(kl)=-u
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 12.0 - - - june, 2016
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/ibc(4),icont,iu(95),a1x,a1y,a1u,a2x,a2y,a2u,
     +          bux,buy,cu0,cu1,cu2,cahn,cexp,cir,csin,du0,du1,eu0,f0,
     1          ru(81),su(100)
cy
        call setrl(rl)
        values(k0)=-eu0
        if(icont==18) values(kl)=-1.0e0_rknd
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 12.0 - - - june, 2016
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)=u*u
        values(ku)=2.0e0_rknd*u
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 12.0 - - - june, 2016
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 12.0 - - - june, 2016
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
        r=(x-0.5e0_rknd)**2+(y-0.5e0_rknd)**2
        s=64.0e0_rknd
        d=1.0e0_rknd/16.0e0_rknd
        t=tanh(s*(d-r))
        values(kf)=t
        values(kad)=t
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 12.0 - - - june, 2016
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 12.0 - - - june, 2016
c
c-----------------------------------------------------------------------
        subroutine setrl(rl)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            character(len=80) :: su
            common /atest2/ibc(4),icont,iu(95),ru(100),su(100)
cy
        if(icont>0) ru(icont)=rl
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 12.0 - - - june, 2016
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), dimension(5) :: isv
            integer(kind=iknd), save :: len
            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
            character(len=80), save, dimension(80) :: file
cy
            data len/56/
            data (file(i),i=  1, 10)/
     +      'n i= 1,n=left,  a= l,t=i',
     1      'n i= 2,n=top,   a= t,t=i',
     2      'n i= 3,n=right, a= r,t=i',
     3      'n i= 4,n=bottom,a= b,t=i',
     4      'n i= 5,n=icont, a= i,t=i',
     5      'n i= 1,n=a1x,   a=x1,t=r',
     6      'n i= 2,n=a1y,   a=y1,t=r',
     7      'n i= 3,n=a1u,   a=u1,t=r',
     8      'n i= 4,n=a2x,   a=x2,t=r',
     9      'n i= 5,n=a2y,   a=y2,t=r'/
            data (file(i),i= 11, 20)/
     +      'n i= 6,n=a2u,   a=u2,t=r',
     1      'n i= 7,n=bux,   a=bx,t=r',
     2      'n i= 8,n=buy,   a=by,t=r',
     3      'n i= 9,n=cu0,   a=c0,t=r',
     4      'n i=10,n=cu1,   a=c1,t=r',
     5      'n i=11,n=cu2,   a=c2,t=r',
     6      'n i=12,n=cahn,  a=ca,t=r',
     7      'n i=13,n=cexp,  a=cx,t=r',
     8      'n i=14,n=cir,   a=cr,t=r',
     9      'n i=15,n=csin,  a=cs,t=r'/
            data (file(i),i= 21, 30)/
     +      'n i=16,n=du0,   a=d0,t=r',
     1      'n i=17,n=du1,   a=d1,t=r',
     2      'n i=18,n=eu0,   a=e0,t=r',
     3      'n i=19,n=f0,    a=f0,t=r',
     4      's n=left  ,v=1,l="neumann"',
     5      's n=left  ,v=2,l="dirichlet"',
     6      's n=left  ,v=0,l="periodic"',
     7      's n=top   ,v=1,l="neumann"',
     8      's n=top   ,v=2,l="dirichlet"',
     9      's n=top   ,v=0,l="periodic"'/
            data (file(i),i= 31, 40)/
     +      's n=right ,v=1,l="neumann"',
     1      's n=right ,v=2,l="dirichlet"',
     2      's n=right ,v=0,l="periodic"',
     3      's n=bottom,v=1,l="neumann"',
     4      's n=bottom,v=2,l="dirichlet"',
     5      's n=bottom,v=0,l="periodic"',
     6      's n=icont ,v= 0,l="none"',
     7      's n=icont ,v= 1,l="a1x"',
     8      's n=icont ,v= 2,l="a1y"',
     9      's n=icont ,v= 3,l="a1u"'/
            data (file(i),i= 41, 50)/
     +      's n=icont ,v= 4,l="a2x"',
     1      's n=icont ,v= 5,l="a2y"',
     2      's n=icont ,v= 6,l="a2u"',
     3      's n=icont ,v= 7,l="bux"',
     4      's n=icont ,v= 8,l="buy"',
     5      's n=icont ,v= 9,l="cu0"',
     6      's n=icont ,v=10,l="cu1"',
     7      's n=icont ,v=11,l="cu2"',
     8      's n=icont ,v=12,l="cahn"',
     9      's n=icont ,v=13,l="cexp"'/
            data (file(i),i= 51, 56)/
     +      's n=icont ,v=14,l="cir"',
     1      's n=icont ,v=15,l="csin"',
     2      's n=icont ,v=16,l="du0"',
     3      's n=icont ,v=17,l="du1"',
     4      's n=icont ,v=18,l="eu0"',
     5      's n=icont ,v=19,l="f0"'/
c
c
c       save integer parameters
c
        do i=1,5
            isv(i)=iu(i)
        enddo
c
c       enter input mode
c
        call usrset(file,len,iu,ru,su)
c
c       if any of the integer parameters have changed, call gdata
c
        do i=1,4
            if(iu(i)<0.or.iu(i)>2) iu(i)=2
            if(iu(i)/=isv(i)) ip(41)=-1
        enddo
        if(iu(5)==0) then
            if(isv(5)/=0) ip(41)=-1
        else
            if(isv(5)==0) ip(41)=-1
            rp(1)=ru(iu(5))
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                    edition 12.0 - - - june, 2016
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(9) :: ix,iy
            integer(kind=iknd), save :: ntf,nvf,nbf,ispd,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 ix/0,0,1,2,2,2,1,0,1/
            data iy/1,2,2,2,1,0,0,0,1/
            data ntf,nvf,nbf,ispd,ifirst/8,9,8,0,1/
c
c           common /atest2/ibc(4),icont,iu(95),a1x,a1y,a1u,a2x,a2y,a2u,
c    +          bux,buy,cu0,cu1,cu2,cahn,cexp,cir,csin,du0,du1,eu0,f0,
c    1          ru(81),su(100)
c
        if(ip(41)==1) then
            sp(2)='square'
            sp(1)='square'
            sp(3)='square'
            sp(6)='square_mpixxx.rw'
            sp(7)='square.jnl'
            sp(9)='square_mpixxx.out'
c
c       initialize as laplacian with dirichlet b.c.
c
            do i=1,19
                ru(i)=0.0e0_rknd
            enddo
            ru(1)=1.0e0_rknd
            ru(5)=1.0e0_rknd
            ru(9)=1.0e0_rknd
            do i=1,4
               iu(i)=2
            enddo
            iu(5)=0
        endif
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        ip(5)=max(ip(5),ifirst)
        ip(6)=1
        if(iu(5)/=0) ip(6)=3
        ip(7)=0
        ip(8)=ispd
        ip(22)=100000
        do i=1,ntf
            itnode(1,i)=9
            itnode(2,i)=i
            itnode(3,i)=i-1
            itnode(4,i)=0
            itnode(5,i)=i
ccc         if(i>4) itnode(5,i)=1
            ibndry(1,i)=i
            ibndry(2,i)=i-1
            ibndry(3,i)=0
            k=(i+1)/2
            ibndry(4,i)=iu(k)
            ibndry(5,i)=0
            ibndry(6,i)=0
            ibndry(7,i)=(i+1)/2
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
        itnode(3,1)=8
        ibndry(2,1)=8
c
        do i=1,nvf
            vx(i)=real(ix(i),rknd)/2.0e0_rknd
            vy(i)=real(iy(i),rknd)/2.0e0_rknd
        enddo
c
c       check for periodic boundary conditions
c
        if(iu(1)==0.and.iu(3)==0) then
            ibndry(4,1)=-6
            ibndry(4,6)=-1
            ibndry(4,2)=-5
            ibndry(4,5)=-2
        endif
        if(iu(2)==0.and.iu(4)==0) then
            ibndry(4,3)=-8
            ibndry(4,8)=-3
            ibndry(4,4)=-7
            ibndry(4,7)=-4
        endif
c
        return
        end