c***********************  problem name: mnsurf  ************************
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
            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 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,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 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
            common /val1/k0,ku,kl,kuu,kul,klu,kll
cy
        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
        theta=ru(1)
c
        if(itag==2.or.itag==4) values(k0)=4.0e0_rknd*x*(1.0e0_rknd-x)
        if (itag==2) then
           values(klb)=1.0e0_rknd
        else if(itag==3) then
            s=(1.0e0_rknd-theta)/2.0e0_rknd
            b=(x-s)/(0.25e0_rknd-s)
            values(klb)=2.0e0_rknd*b-1.0e0_rknd
        else if(itag==4) then
            s=(1.0e0_rknd+theta)/2.0e0_rknd
            b=(y-s)/(0.75e0_rknd-s)
            values(klb)=2.0e0_rknd*b-1.0e0_rknd
        else if(itag==5) then
            s=(1.0e0_rknd+theta)/2.0e0_rknd
            b=(x-s)/(0.75e0_rknd-s)
            values(klb)=2.0e0_rknd*b-1.0e0_rknd
        else if(itag==6) then
            s=(1.0e0_rknd-theta)/2.0e0_rknd
            b=(y-s)/(0.25e0_rknd-s)
            values(klb)=2.0e0_rknd*b-1.0e0_rknd
        else
           values(klb)=-1.0e0_rknd
        endif
        values(kub)=1.5e0_rknd
        values(kic)=(values(klb)+values(kub))/2.0e0_rknd
        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
            character(len=80) :: su
            common /val0/k0,ku,kx,ky,kl,kuu,kux,kuy,kul,
     +          kxu,kxx,kxy,kxl,kyu,kyx,kyy,kyl,klu,klx,kly,kll
            common /atest2/iu(100),ru(100),su(100)
cy
        dshift=ru(2)
        ss=sqrt(1.0e0_rknd+ux**2+uy**2)
        values(k0)=ss
        values(kx)=ux/ss
        values(ky)=uy/ss
        values(kxx)=(1.0e0_rknd+uy**2)/ss**3+dshift
        values(kxy)=-ux*uy/ss**3
        values(kyy)=(1.0e0_rknd+ux**2)/ss**3+dshift
        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
            common /val3/kf,kf1,kf2,kad
cy
        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
        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
            character(len=80), dimension(100) :: sp,su
            character(len=80), save, dimension(10) :: file
cy
            data len/2/
            data (file(i),i=  1, 2)/
     +      'n i= 1,n=theta, a=t ,t=r',
     1      'n i= 2,n=dshift,a=d ,t=r'/
c
c       enter input mode
c
        rsv=ru(1)
        call usrset(file,len,iu,ru,su)
        if(ru(1)/=rsv) ip(41)=-1
        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, dimension(9) :: ix,iy
            integer(kind=iknd), dimension(200) :: jb
            integer(kind=iknd), save :: ntf,nvf,nbf,ispd,iprob
            integer(kind=iknd), save :: iadapt,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
            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/24,25,48,1/
            data iprob,iadapt,ifirst/2,5,1/
            data hmax,grade/0.1e0_rknd,1.5e0_rknd/
c
        if(ip(41)==1) then
            sp(1)='minimal surface'
            sp(2)='minimal surface'
            sp(3)='minimal surface'
            sp(6)='mnsurf.rw'
            sp(7)='mnsurf.jnl'
            sp(9)='mnsurf.out'
c
            ru(1)=0.55e0_rknd
            ru(2)=1.0e-2_rknd
        endif
c
        rp(3)=1.0e-1_rknd
c
        do i=1,9
            vx(i)=real(ix(i),rknd)/2.0e0_rknd
            vy(i)=real(iy(i),rknd)/2.0e0_rknd
        enddo
        t0=ru(1)/2.0e0_rknd
        t1=(1.0e0_rknd-ru(1))/2.0e0_rknd
        do i=1,8
            vx(9+i)=(real(ix(i)+ix(9),rknd))/4.0e0_rknd
            vy(9+i)=(real(iy(i)+iy(9),rknd))/4.0e0_rknd
            vx(17+i)=t0*real(ix(i),rknd)+t1*real(ix(9),rknd)
            vy(17+i)=t0*real(iy(i),rknd)+t1*real(iy(9),rknd)
        enddo
        do i=1,nbf
            do j=1,7
                ibndry(j,i)=0
            enddo
            sf(1,i)=0.0e0_rknd
            sf(2,i)=0.0e0_rknd
        enddo
        do i=1,8
            ibndry(1,i)=i
            ibndry(2,i)=i-1
            ibndry(4,i)=2
            ibndry(7,i)=(i+1)/2
c
            ibndry(1,i+8)=i
            ibndry(2,i+8)=i+17
            ibndry(1,i+16)=i+17
            ibndry(2,i+16)=i+9
            ibndry(1,i+24)=i+9
            ibndry(2,i+24)=9
            ibndry(1,i+32)=i+17
            ibndry(2,i+32)=i+16
            ibndry(1,i+40)=i+9
            ibndry(2,i+40)=i+8
        enddo
        ibndry(2,1)=8
        ibndry(2,33)=25
        ibndry(2,41)=17
c
c
        ip(1)=ntf
        ip(2)=nvf
        ip(3)=nbf
        ip(5)=max(ip(5),ifirst)
        ip(6)=iprob
        ip(8)=ispd
        rp(15)=hmax
        rp(16)=grade
        ip(19)=1
        ip(20)=iadapt
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
c       compute correct region labels
c
        rl=0.0e0_rknd
        call makjb(nvf,nbf,ntf,vx,vy,sf,ibndry,itnode,1,jb,
     +      iflag,rl,sxy)
c
        do i=1,ntf
           do j=jb(i),jb(i+1)-1
                if(jb(j)==10) itnode(5,i)=1
                if(jb(j)==12) itnode(5,i)=1
                if(jb(j)==14) itnode(5,i)=1
                if(jb(j)==16) itnode(5,i)=1
                if(jb(j)==26) itnode(5,i)=2
                if(jb(j)==28) itnode(5,i)=2
                if(jb(j)==30) itnode(5,i)=2
                if(jb(j)==32) itnode(5,i)=2
                if(jb(j)==17) itnode(5,i)=3
                if(jb(j)==19) itnode(5,i)=4
                if(jb(j)==21) itnode(5,i)=5
                if(jb(j)==23) itnode(5,i)=6
            enddo
        enddo
c
        return
        end

.