c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine mpiutl(isw)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c       initialize
c
        if(isw==1) then
            call MPI_INIT(ierr)
            call MPI_COMM_RANK(MPI_COMM_WORLD,id,ierr)
            call MPI_COMM_SIZE(MPI_COMM_WORLD,np,ierr)
            nproc=np
            myid=id
            mpirgn=0
            if(nproc<=1) then
                mpisw=-1
            else
                mpisw=1
            endif
            if(rknd==rsngl) then
                mpiflt=MPI_REAL
            else if(rknd==rdble) then
                mpiflt=MPI_DOUBLE_PRECISION
            else
                mpiflt=MPI_LONG_DOUBLE
            endif
            if(iknd==isngl) then
                mpiint=MPI_INTEGER
            else
                mpiint=MPI_LONG_LONG_INT
            endif
        else
c
c       quit
c
            call MPI_FINALIZE(ierr)
        endif
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine star0(cmd)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            integer(kind=iknd), save :: mpijnl=0,lowera=97,lowerz=122
            character(len=6) :: cmdtyp
            character(len=80) :: list,cmd,cmd0
            common /atest3/mode,jnlsw,jnlr,jnlw,ibatch
            common /atest4/jcmd,cmdtyp,list
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c       master
c
        if(myid==0) then
            cmd0=cmd
            if(mode==jnlsw) mpijnl=0
            if(mode==1) mpijnl=1
            if(mpijnl==1.and.cmdtyp/='mpicmd') return
            if(mpisw<0) then
                if(cmdtyp=='popup '.or.cmdtyp=='file  '
     +              .or.cmdtyp=='journl') then
                    ii=ichar(cmd0(1:1))
                    if(ii>=lowera.and.ii<=lowerz) then
                        cmd0(1:1)=char(ii-32)
                    endif
                endif
            else
                if(cmdtyp=='journl') mpijnl=1
            endif
            call MPI_BCAST(cmd0,80,MPI_CHARACTER,0,MPI_COMM_WORLD,ierr)
        else
c
c       slave
c
            if(jnlsw>=1.and.cmdtyp/='mpicmd') return
            call MPI_BCAST(cmd,80,MPI_CHARACTER,0,MPI_COMM_WORLD,ierr)
        endif
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine expth(ip,ipath)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            integer(kind=iknd), dimension(6,*) :: ipath
            integer(kind=iknd), dimension(100) :: ip
            integer(kind=iknd), allocatable, dimension(:,:) :: ipath0
            integer(kind=iknd), allocatable, dimension(:) :: bsize
            integer(kind=isngl), allocatable, dimension(:) ::
     +          counts,offset
            integer(kind=iknd), dimension(4) :: id,jd
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c       condense data structure
c
        maxpth=ip(82)
        allocate(ipath0(6,maxpth),counts(nproc),offset(nproc),
     +      bsize(nproc))
        ishift=nproc+1
        ipath(3,1)=ipath(3,nproc+2)
        ipath(4,1)=ipath(4,nproc+2)
        ipath(1,1)=ipath(1,nproc+2)-ishift
        ipath(2,1)=ipath(2,nproc+2)-ishift
        do j=ipath(1,1),ipath(2,1)
            do k=1,6
                ipath(k,j)=ipath(k,j+ishift)
            enddo
            if(ipath(2,j)>0) ipath(2,j)=ipath(2,j)-ishift
        enddo
        len=ipath(2,1)*6
c
        call MPI_ALLGATHER(len,1,mpiint,bsize,1,
     +      mpiint,MPI_COMM_WORLD,ierr)
c
        offset(1)=0
        do irgn=1,nproc-1
            offset(irgn+1)=offset(irgn)+bsize(irgn)
            counts(irgn)=bsize(irgn)
        enddo
        counts(nproc)=bsize(nproc)
        isum=(offset(nproc)+counts(nproc))/6+nproc+2
        if(isum>maxpth) then
            ip(25)=72
            go to 20
        endif
c
        call MPI_ALLGATHERV(ipath,len,mpiint,
     +      ipath0,counts,offset,mpiint,MPI_COMM_WORLD,ierr)
c
        ipath(3,nproc+2)=0
        ipath(4,nproc+2)=0
        ipath(1,nproc+2)=nproc+3
        ipath(2,nproc+2)=nproc+2
        do irgn=1,nproc
            ii=offset(irgn)/6
c
            ipath(3,nproc+2)=max(ipath(3,nproc+2),ipath0(3,1+ii))
c
            n0=ipath(4,nproc+2)
            ipath(3,irgn)=n0+1
            ipath(4,irgn)=n0+ipath0(4,1+ii)
            ipath(4,nproc+2)=ipath(4,irgn)
c
            ipath(1,irgn)=ipath(2,nproc+2)+1
            jb0=ipath(1,irgn)-ipath0(1,1+ii)
            ipath(2,irgn)=ipath0(2,1+ii)+jb0
            ipath(2,nproc+2)=ipath(2,irgn)
c
            do j=ipath(1,irgn),ipath(2,irgn)
                do k=1,6
                    ipath(k,j)=ipath0(k,j-jb0+ii)
                enddo
                if(ipath(2,j)>0) then
                    ipath(2,j)=ipath(2,j)+jb0
                endif
                ipath(3,j)=ipath(3,j)+n0
                ipath(4,j)=ipath(4,j)+n0
                if(ipath(5,j)>0) then
                    ipath(5,j)=ipath(5,j)+n0
                else if(ipath(5,j)<0) then
                    ipath(5,j)=ipath(5,j)-n0
                endif
            enddo
        enddo
c
        ip(74)=ipath(4,nproc+2)
        ip(75)=ipath(2,nproc+2)
c
c       global dimensions
c
        id(1)=ip(27)
        id(2)=ip(28)
        id(3)=ip(29)
        id(4)=ip(30)
        call MPI_ALLREDUCE(id,jd,4,mpiint,MPI_SUM,MPI_COMM_WORLD,ierr)
        ip(37)=jd(1)
        ip(38)=jd(2)
        ip(39)=jd(3)
        ip(40)=jd(4)
c
   20   deallocate(ipath0,counts,offset,bsize)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine exbdy(ipath,ir0,map,gf,nn,num)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            integer(kind=iknd), dimension(6,*) :: ipath
            integer(kind=iknd), dimension(*) :: ir0,map
            integer(kind=isngl), allocatable, dimension(:) ::
     +          counts,offset
            real(kind=rknd), allocatable,  dimension(:) :: vin,vout
            real(kind=rknd), dimension(nn,*) :: gf
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c       communicate interface data
c
        allocate(vin(nn*num),vout(nn*num),counts(nproc),offset(nproc))
        do irgn=1,nproc
            ii=ipath(3,irgn)-1
            nvdd=ipath(4,irgn)-ii
            counts(irgn)=num*nvdd
            offset(irgn)=num*ii
        enddo
c
        irgn=myid+1
        ii=ipath(3,irgn)-1
        nvdd=ipath(4,irgn)-ii
        len=counts(irgn)
        do j=1,num
            jj=(j-1)*nvdd
            do i=1,nvdd
                vin(i+jj)=gf(i+ii,j)
            enddo
        enddo
c
        call MPI_ALLGATHERV(vin,len,mpiflt,
     +      vout,counts,offset,mpiflt,MPI_COMM_WORLD,ierr)
c
        do irgn=1,nproc
            ii=ipath(3,irgn)-1
            nvdd=ipath(4,irgn)-ii
            do j=1,num
                jj=(j-1)*nvdd+offset(irgn)
                 do i=1,nvdd
                     gf(i+ii,j)=vout(i+jj)
                 enddo
             enddo
        enddo
c
c       interpolate to mesh on this processor
c
        lipath=ipath(2,nproc+1)
        call intrpi(lipath,ipath,ir0,map,nn,num,gf,vin)
c
        deallocate(vin,vout,counts,offset)
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine pl2ip(t,num)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            real(kind=rknd), dimension(*) :: t
            real(kind=rknd), dimension(200) :: temp
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c       compute all inner products/norms
c
        do i=1,num
            temp(i)=t(i)
        enddo
c
        call MPI_ALLREDUCE(temp,t,num,mpiflt,
     +      MPI_SUM,MPI_COMM_WORLD,ierr)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine bcast(vx,vy,sf,ibndry,itnode,itdof,
     +      ip,rp,sp,iu,ru,su,gf,e)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(100) :: ip,iu
            integer(kind=iknd), dimension(8,*) :: itdof
            real(kind=rknd), dimension(*) :: vx,vy,gf,sf,e
            real(kind=rknd), dimension(100) :: rp,ru
            character(len=1), allocatable, dimension(:) :: buff
            character(len=80), dimension(100) :: sp,su
            common /pltmg7/time(3,50),hist(22,30)
            common /pltmg6/path(101,6)
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c
c       inital fan-out of coarse grid data
c
c       master node
c
        if(myid==0) then
            ntf=ip(1)
            nvf=ip(2)
            nbf=ip(3)
            ndf=ip(4)
            maxt=ip(83)
            maxd=ip(85)
            ngf=ip(77)
c
            leni=(200+13*ntf+7*nbf)*iknd
            lenr=(200+ngf*ndf+2*nvf+2*nbf+2*ntf+606+660+150)*rknd
            lens=200*80
            len=leni+lenr+lens+16
c
c       integers
c
            allocate(buff(len))
            lenz=len
            len=0
            call MPI_PACK(ip,100,mpiint,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            call MPI_PACK(iu,100,mpiint,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            call MPI_PACK(itnode,5*ntf,mpiint,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            call MPI_PACK(ibndry,7*nbf,mpiint,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            call MPI_PACK(itdof,8*ntf,mpiint,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
c
c       floats
c
            call MPI_PACK(rp,100,mpiflt,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            call MPI_PACK(ru,100,mpiflt,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            call MPI_PACK(sf,2*nbf,mpiflt,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            call MPI_PACK(vx,nvf,mpiflt,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            call MPI_PACK(vy,nvf,mpiflt,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            do k=1,ngf
                iptr=1+(k-1)*maxd
                call MPI_PACK(gf(iptr),ndf,mpiflt,
     +              buff,lenz,len,MPI_COMM_WORLD,ierr)
            enddo
            do k=1,2
                iptr=1+(k-1)*maxt
                call MPI_PACK(e(iptr),ntf,mpiflt,
     +              buff,lenz,len,MPI_COMM_WORLD,ierr)
            enddo
            call MPI_PACK(path,606,mpiflt,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            call MPI_PACK(hist,660,mpiflt,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            call MPI_PACK(time,150,mpiflt,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
c
c       characters
c
            call MPI_PACK(sp,8000,MPI_CHARACTER,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            call MPI_PACK(su,8000,MPI_CHARACTER,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
c
c       the big broadcast
c
            call MPI_BCAST(len,1,mpiint,0,MPI_COMM_WORLD,ierr)
            call MPI_BCAST(buff,len,MPI_PACKED,0,MPI_COMM_WORLD,ierr)
            deallocate(buff)
c
c       slave node
c
        else
            call MPI_BCAST(len,1,mpiint,0,MPI_COMM_WORLD,ierr)
            allocate(buff(len))
            call MPI_BCAST(buff,len,MPI_PACKED,0,MPI_COMM_WORLD,ierr)
c
c       integers
c
            lenz=len
            len=0
c
            call MPI_UNPACK(buff,lenz,len,ip,100,
     +          mpiint,MPI_COMM_WORLD,ierr)
            ip(48)=mpisw
            ip(49)=nproc
            ip(50)=myid+1
            ntf=ip(1)
            nvf=ip(2)
            nbf=ip(3)
            ndf=ip(4)
            maxt=ip(83)
            maxd=ip(85)
            ngf=ip(77)
c
            call MPI_UNPACK(buff,lenz,len,iu,100,
     +          mpiint,MPI_COMM_WORLD,ierr)
            call MPI_UNPACK(buff,lenz,len,itnode,5*ntf,
     +          mpiint,MPI_COMM_WORLD,ierr)
            call MPI_UNPACK(buff,lenz,len,ibndry,7*nbf,
     +          mpiint,MPI_COMM_WORLD,ierr)
            call MPI_UNPACK(buff,lenz,len,itdof,8*ntf,
     +          mpiint,MPI_COMM_WORLD,ierr)
c
c       floats
c
            call MPI_UNPACK(buff,lenz,len,rp,100,
     +          mpiflt,MPI_COMM_WORLD,ierr)
            call MPI_UNPACK(buff,lenz,len,ru,100,
     +          mpiflt,MPI_COMM_WORLD,ierr)
            call MPI_UNPACK(buff,lenz,len,sf,2*nbf,
     +          mpiflt,MPI_COMM_WORLD,ierr)
            call MPI_UNPACK(buff,lenz,len,vx,nvf,
     +          mpiflt,MPI_COMM_WORLD,ierr)
            call MPI_UNPACK(buff,lenz,len,vy,nvf,
     +          mpiflt,MPI_COMM_WORLD,ierr)
            do k=1,ngf
                iptr=1+(k-1)*maxd
                call MPI_UNPACK(buff,lenz,len,gf(iptr),ndf,
     +              mpiflt,MPI_COMM_WORLD,ierr)
            enddo
            do k=1,2
                iptr=1+(k-1)*maxt
                call MPI_UNPACK(buff,lenz,len,e(iptr),ntf,
     +              mpiflt,MPI_COMM_WORLD,ierr)
            enddo
            call MPI_UNPACK(buff,lenz,len,path,606,
     +          mpiflt,MPI_COMM_WORLD,ierr)
            call MPI_UNPACK(buff,lenz,len,hist,660,
     +          mpiflt,MPI_COMM_WORLD,ierr)
            call MPI_UNPACK(buff,lenz,len,time,150,
     +          mpiflt,MPI_COMM_WORLD,ierr)
c
c       characters
c
            call MPI_UNPACK(buff,lenz,len,sp,8000,
     +          MPI_CHARACTER,MPI_COMM_WORLD,ierr)
            call MPI_UNPACK(buff,lenz,len,su,8000,
     +          MPI_CHARACTER,MPI_COMM_WORLD,ierr)
            deallocate(buff)
        endif
c
c       this uses many bcasts but no buffer storage
c
c       call MPI_BCAST(ip,100,mpiint,0,MPI_COMM_WORLD,ierr)
c       ip(48)=mpisw
c       ip(49)=nproc
c       ip(50)=myid+1
c       ntf=ip(1)
c       nvf=ip(2)
c       nbf=ip(3)
c       ndf=ip(4)
c       maxt=ip(83)
c       maxd=ip(85)
c       ngf=ip(77)
c
c       call MPI_BCAST(iu,100,mpiint,0,MPI_COMM_WORLD,ierr)
c       call MPI_BCAST(itnode,5*ntf,mpiint,0,MPI_COMM_WORLD,ierr)
c       call MPI_BCAST(ibndry,7*nbf,mpiint,0,MPI_COMM_WORLD,ierr)
c       call MPI_BCAST(itdof,8*ntf,mpiint,0,MPI_COMM_WORLD,ierr)
c
c       floats
c
c       call MPI_BCAST(rp,100,mpiflt,0,MPI_COMM_WORLD,ierr)
c       call MPI_BCAST(ru,100,mpiflt,0,MPI_COMM_WORLD,ierr)
c       call MPI_BCAST(sf,2*nbf,mpiflt,0,MPI_COMM_WORLD,ierr)
c       call MPI_BCAST(vx,nvf,mpiflt,0,MPI_COMM_WORLD,ierr)
c       call MPI_BCAST(vy,nvf,mpiflt,0,MPI_COMM_WORLD,ierr)
c       do k=1,ngf
c           iptr=1+(k-1)*maxd
c           call MPI_BCAST(gf(iptr),ndf,mpiflt,0,MPI_COMM_WORLD,ierr)
c       enddo
c       do k=1,2
c           iptr=1+(k-1)*maxt
c           call MPI_BCAST(e(iptr),ntf,mpiflt,0,MPI_COMM_WORLD,ierr)
c       enddo
c       call MPI_BCAST(path,606,mpiflt,0,MPI_COMM_WORLD,ierr)
c       call MPI_BCAST(hist,660,mpiflt,0,MPI_COMM_WORLD,ierr)
c       call MPI_BCAST(time,150,mpiflt,0,MPI_COMM_WORLD,ierr)
c
c       characters
c
c       call MPI_BCAST(sp,8000,MPI_CHARACTER,0,MPI_COMM_WORLD,ierr)
c       call MPI_BCAST(su,8000,MPI_CHARACTER,0,MPI_COMM_WORLD,ierr)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine exstat(rp,pstat)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            real(kind=rknd), dimension(100) :: rp
            real(kind=rknd), dimension(4,*) :: pstat
            real(kind=rknd), dimension(4) :: sendb
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c       exchange load balance statistics
c
        sendb(1)=rp(82)
        sendb(2)=rp(83)
        sendb(3)=rp(84)
        sendb(4)=rp(85)
c
        call MPI_ALLGATHER(sendb,4,mpiflt,
     +      pstat,4,mpiflt,MPI_COMM_WORLD,ierr)
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine exbox(bmin,bmax,num)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            real(kind=rknd), dimension(5) :: tin,tout
            real(kind=rknd), dimension(*) :: bmin,bmax
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c       compute the global bounding box
c
        do i=1,num
            tin(i)=bmin(i)
        enddo
        call MPI_ALLREDUCE(tin,tout,num,mpiflt,
     +      MPI_MIN,MPI_COMM_WORLD,ierr)
c
        do i=1,num
            bmin(i)=tout(i)
            tin(i)=bmax(i)
        enddo
c
        call MPI_ALLREDUCE(tin,tout,num,mpiflt,
     +      MPI_MAX,MPI_COMM_WORLD,ierr)
        do i=1,num
            bmax(i)=tout(i)
        enddo
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine exdist(kdist,num)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            integer(kind=iknd), dimension(*) :: kdist
            integer(kind=iknd), dimension(22) :: itemp
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c       compute all inner products/norms
c
        do i=1,num
            itemp(i)=kdist(i)
        enddo
c
        call MPI_ALLREDUCE(itemp,kdist,num,mpiint,
     +      MPI_SUM,MPI_COMM_WORLD,ierr)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine exqual(num,val)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            integer(kind=iknd), dimension(5) :: num,itemp
            real(kind=rknd), dimension(2) :: temp,val
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c       compute global mesh quality statistics (binit)
c
        do i=1,5
            itemp(i)=num(i)
        enddo
        call MPI_ALLREDUCE(itemp,num,5,mpiint,
     +      MPI_SUM,MPI_COMM_WORLD,ierr)
c
        do i=1,2
            temp(i)=val(i)
        enddo
        call MPI_ALLREDUCE(temp,val,1,mpiflt,
     +      MPI_MIN,MPI_COMM_WORLD,ierr)
        call MPI_ALLREDUCE(temp(2),val(2),1,mpiflt,
     +      MPI_SUM,MPI_COMM_WORLD,ierr)
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine exibox(ibmin,ibmax,num)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            integer(kind=iknd), dimension(5) :: itemp
            integer(kind=iknd), dimension(*) :: ibmin,ibmax
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c       compute min and max color (binit)
c
        do i=1,num
            itemp(i)=ibmin(i)
        enddo
        call MPI_ALLREDUCE(itemp,ibmin,num,mpiint,
     +      MPI_MIN,MPI_COMM_WORLD,ierr)
c
        do i=1,num
            itemp(i)=ibmax(i)
        enddo
        call MPI_ALLREDUCE(itemp,ibmax,num,mpiint,
     +      MPI_MAX,MPI_COMM_WORLD,ierr)
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine exflag(iflag)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c       global dimensions
c
        itemp=iflag
        call MPI_ALLREDUCE(itemp,iflag,1,mpiint,
     +      MPI_MAX,MPI_COMM_WORLD,ierr)
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine exstep(step)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c       compute the maximum step (tpickd)
c
        temp=step
        call MPI_ALLREDUCE(temp,step,1,mpiflt,
     +      MPI_MIN,MPI_COMM_WORLD,ierr)
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine extim(atime,ptime)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            integer(kind=iknd), save :: len=50
            real(kind=rknd), dimension(*) :: ptime
            real(kind=rknd), dimension(3,50) :: temp,atemp
            real(kind=rknd), dimension(50) :: tmax,tmin,smax,smin
            real(kind=rknd), dimension(5,*) :: atime
            real(kind=rknd), dimension(2) :: sendb
            common /pltmg7/time(3,50),hist(22,30)
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c     exchange times
c
        sum=0.0e0_rknd
        do i=1,len
            sum=sum+time(2,i)
            temp(1,i)=time(1,i)
            temp(2,i)=time(2,i)
            temp(3,i)=time(3,i)
            tmax(i)=time(2,i)
            tmin(i)=time(2,i)
        enddo
c
        sendb(1)=sum
        call MPI_ALLGATHER(sendb,1,mpiflt,
     +      ptime,1,mpiflt,MPI_COMM_WORLD,ierr)
c
        num=3*len
        call MPI_ALLREDUCE(temp,atemp,num,mpiflt,
     +      MPI_SUM,MPI_COMM_WORLD,ierr)
        call MPI_ALLREDUCE(tmax,smax,len,mpiflt,
     +      MPI_MAX,MPI_COMM_WORLD,ierr)
        call MPI_ALLREDUCE(tmin,smin,len,mpiflt,
     +      MPI_MIN,MPI_COMM_WORLD,ierr)
c
        do i=1,len
            atime(1,i)=atemp(1,i)/real(nproc,rknd)
            atime(2,i)=atemp(2,i)/real(nproc,rknd)
            atime(3,i)=atemp(3,i)/real(nproc,rknd)
            atime(4,i)=smin(i)
            atime(5,i)=smax(i)
        enddo
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine exerr(hist,numhst,e)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            real(kind=rknd), dimension(*):: e
            real(kind=rknd), dimension(22,*) :: hist
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c     collect convergence data from all processes
c
        ee=0.0e0_rknd
        mxhist=20
        do j=1,numhst
            num=min(int(hist(mxhist+1,j)),mxhist)
            if(num>0) then
                e1=abs(hist(mxhist+2,j))
                if(e1>0.0e0_rknd) e1=1.0e0_rknd/e1
                qq=abs(hist(num,j))*e1
                if(qq>0.0e0_rknd) qq=log10(qq)
                if(ee==0.0e0_rknd) ee=qq
                ee=max(ee,qq)
            endif
        enddo
c
        call MPI_GATHER(ee,1,mpiflt,e,1,mpiflt,0,MPI_COMM_WORLD,ierr)
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine exsze(jd,isw)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            integer(kind=iknd), dimension(3) :: id,jd
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c       global dimensions
c
        id(1)=jd(1)
        id(2)=jd(2)
        id(3)=jd(3)
        if(isw==0) then
            call MPI_ALLREDUCE(id,jd,3,mpiint,
     +          MPI_MAX,MPI_COMM_WORLD,ierr)
        else
            call MPI_ALLREDUCE(id,jd,3,mpiint,
     +          MPI_SUM,MPI_COMM_WORLD,ierr)
        endif
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine exint(ia,len)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            integer(kind=isngl), dimension(MPI_STATUS_SIZE) :: status
            integer(kind=iknd), dimension(*) :: ia
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c       2-way integer transfer
c
        if(myid==0) then
            idsrc=mpirgn-1
            call MPI_RECV(ia,len,mpiint,idsrc,
     +          0,MPI_COMM_WORLD,status,ierr)
         else
            call MPI_SEND(ia,len,mpiint,0,
     +          0,MPI_COMM_WORLD,ierr)
        endif
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine exflt(ra,len)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            integer(kind=isngl), dimension(MPI_STATUS_SIZE) :: status
            real(kind=rknd), dimension(*) :: ra
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c       2-way real transfer
c
        if(myid==0) then
            idsrc=mpirgn-1
            call MPI_RECV(ra,len,mpiflt,idsrc,
     +          0,MPI_COMM_WORLD,status,ierr)
         else
            call MPI_SEND(ra,len,mpiflt,0,
     +          0,MPI_COMM_WORLD,ierr)
        endif
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine exchr(ca,ilen)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            integer(kind=isngl), dimension(MPI_STATUS_SIZE) :: status
            character(len=1), dimension(*) :: ca
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c       2-way character transfer
c
        if(myid==0) then
            idsrc=mpirgn-1
            call MPI_RECV(ca,ilen,MPI_CHARACTER,idsrc,
     +          0,MPI_COMM_WORLD,status,ierr)
         else
            call MPI_SEND(ca,ilen,MPI_CHARACTER,0,
     +          0,MPI_COMM_WORLD,ierr)
        endif
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine glbpix(vx,vy,sf,ibndry,itnode,icolor,
     +      iuvptr,ut,vt,jp,isw,jsw)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            integer(kind=iknd), dimension(7,*) :: ibndry
            integer(kind=iknd), dimension(5,*) :: itnode
            integer(kind=iknd), dimension(25) :: jp
            integer(kind=iknd), dimension(5) :: ib
            integer(kind=iknd), dimension(*) :: icolor
            integer(kind=isngl), dimension(MPI_STATUS_SIZE) :: status
            integer(kind=iknd), dimension(2,*) :: iuvptr
            real(kind=rknd), dimension(2,*) :: sf
            real(kind=rknd), dimension(*) :: vx,vy,ut,vt
            character(len=1), allocatable, dimension(:) :: buff
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c       collect global graph on master process
c
        itag=4
c
        if(myid==0) then
            if(jsw==1) then
                num=nproc
            else
                num=2
            endif
c
c       integers
c
            do i=1,num-1
                if(jsw==1) then
                    ntf=jp(1)
                    nvf=jp(2)
                    nbf=jp(3)
                    ndf=jp(24)
                    if(isw/=0) isv=iuvptr(1,ntf+1)-1
                else
                    ntf=0
                    nvf=0
                    nbf=0
                    ndf=0
                    isv=0
                endif
c
                call MPI_RECV(len,1,mpiint,MPI_ANY_SOURCE,
     +              itag,MPI_COMM_WORLD,status,ierr)
                allocate(buff(len))
                idsrc=status(MPI_SOURCE)
                call MPI_RECV(buff,len,MPI_PACKED,
     +              idsrc,itag,MPI_COMM_WORLD,status,ierr)
c
                lenz=len
                len=0
                call MPI_UNPACK(buff,lenz,len,
     +              ib,5,mpiint,MPI_COMM_WORLD,ierr)
c
                ntf0=ib(1)
                nvf0=ib(2)
                nbf0=ib(3)
                ndf0=ib(4)
                nuv0=ib(5)
c
                call MPI_UNPACK(buff,lenz,len,itnode(1,ntf+1),5*ntf0,
     +              mpiint,MPI_COMM_WORLD,ierr)
                call MPI_UNPACK(buff,lenz,len,ibndry(1,nbf+1),7*nbf0,
     +              mpiint,MPI_COMM_WORLD,ierr)
                if(isw==0) then
                    call MPI_UNPACK(buff,lenz,len,icolor(ntf+1),ntf0,
     +                  mpiint,MPI_COMM_WORLD,ierr)
                else
                    call MPI_UNPACK(buff,lenz,len,iuvptr(1,ntf+1),
     +                  2*(ntf0+1),mpiint,MPI_COMM_WORLD,ierr)
                endif
c
c       update itnode and ibndry
c
                if(jsw==1) then
                    do k=ntf+1,ntf+ntf0
                        do j=1,3
                            itnode(j,k)=itnode(j,k)+nvf
                        enddo
                    enddo
                    do k=nbf+1,nbf+nbf0
                        do j=1,2
                            ibndry(j,k)=ibndry(j,k)+nvf
                        enddo
                        if(ibndry(4,k)<0) then
                            ibndry(4,k)=ibndry(4,k)-nbf
                        endif
                    enddo
c
c       fixup iuvptr
c
                    if(isw/=0) then
                        do k=ntf+1,ntf+ntf0+1
                            iuvptr(1,k)=iuvptr(1,k)+isv
                        enddo
                    endif
                endif
c
c       floats
c
                call MPI_UNPACK(buff,lenz,len,
     +              sf(1,nbf+1),2*nbf0,mpiflt,MPI_COMM_WORLD,ierr)
                call MPI_UNPACK(buff,lenz,len,
     +              vx(nvf+1),nvf0,mpiflt,MPI_COMM_WORLD,ierr)
                call MPI_UNPACK(buff,lenz,len,
     +              vy(nvf+1),nvf0,mpiflt,MPI_COMM_WORLD,ierr)
                if(isw/=0) then
                    call MPI_UNPACK(buff,lenz,len,
     +                  ut(isv+1),nuv0,mpiflt,MPI_COMM_WORLD,ierr)
                    call MPI_UNPACK(buff,lenz,len,
     +                  vt(isv+1),nuv0,mpiflt,MPI_COMM_WORLD,ierr)
                endif
c
                jp(1)=ntf+ntf0
                jp(2)=nvf+nvf0
                jp(3)=nbf+nbf0
                jp(24)=ndf+ndf0
c
                deallocate(buff)
c
            enddo
        else
c
c       slave node
c
c       integers
c
            ntf=jp(1)
            nvf=jp(2)
            nbf=jp(3)
            ndf=jp(24)
            nuv=0
            if(isw/=0) nuv=iuvptr(1,ntf+1)-iuvptr(1,1)
c
            lenz=(5*ntf+7*nbf+5)*iknd+2*(nbf+nvf)*rknd
            if(isw==0) lenz=lenz+ntf*iknd
            if(isw/=0) lenz=lenz+(ntf+1)*2*iknd+rknd*2*nuv
c
            allocate(buff(lenz))
c
            ib(1)=ntf
            ib(2)=nvf
            ib(3)=nbf
            ib(4)=ndf
            ib(5)=nuv
            len=0
            call MPI_PACK(ib,5,mpiint,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            call MPI_PACK(itnode,5*ntf,mpiint,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            call MPI_PACK(ibndry,7*nbf,mpiint,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            if(isw==0) then
                call MPI_PACK(icolor,ntf,mpiint,
     +              buff,lenz,len,MPI_COMM_WORLD,ierr)
            else
                call MPI_PACK(iuvptr,2*(ntf+1),mpiint,
     +              buff,lenz,len,MPI_COMM_WORLD,ierr)
            endif
c
c       floats
c
            call MPI_PACK(sf,2*nbf,mpiflt,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            call MPI_PACK(vx,nvf,mpiflt,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            call MPI_PACK(vy,nvf,mpiflt,
     +          buff,lenz,len,MPI_COMM_WORLD,ierr)
            if(isw/=0) then
                call MPI_PACK(ut,nuv,mpiflt,
     +              buff,lenz,len,MPI_COMM_WORLD,ierr)
                call MPI_PACK(vt,nuv,mpiflt,
     +              buff,lenz,len,MPI_COMM_WORLD,ierr)
            endif
c
            call MPI_SEND(len,1,mpiint,
     +          0,itag,MPI_COMM_WORLD,ierr)
            call MPI_SEND(buff,len,MPI_PACKED,
     +          0,itag,MPI_COMM_WORLD,ierr)
c
            deallocate(buff)
c
        endif
c
        return
        end
c-----------------------------------------------------------------------
c
c            piecewise lagrange triangle multi grid package
c
c                  edition 13.0 - - - september, 2018
c
c-----------------------------------------------------------------------
        subroutine exgph(ip,rp,sp,hist,time,path)
cx
            use mthdef
            implicit real(kind=rknd) (a-h,o-z)
            implicit integer(kind=iknd) (i-n)
            include "mpif.h"
            integer(kind=iknd), dimension(100) :: ip
            integer(kind=isngl), dimension(MPI_STATUS_SIZE) :: status
            real(kind=rknd), dimension(100) :: rp
            real(kind=rknd), dimension(3,50) :: time
            real(kind=rknd), dimension(22,30) :: hist
            real(kind=rknd), dimension(101,6) :: path
            character(len=1), allocatable, dimension(:) :: buff
            character(len=80), dimension(100) :: sp
            common /atest6/nproc,myid,mpisw,mpirgn,mpiint,mpiflt
cy
c
c       inital fan-out of coarse grid data
c
        leni=100*iknd
        lenr=(100+606+660+150)*rknd
        lens=100*80
        len0=leni+lenr+lens
        allocate(buff(len0))
c
c       slave  node
c
        if(myid/=0) then
c
c       integers
c
            lenz=len0
            len1=0
            call MPI_PACK(ip,100,mpiint,
     +          buff,lenz,len1,MPI_COMM_WORLD,ierr)
c
c       floats
c
            call MPI_PACK(rp,100,mpiflt,
     +          buff,lenz,len1,MPI_COMM_WORLD,ierr)
            call MPI_PACK(time,150,mpiflt,
     +          buff,lenz,len1,MPI_COMM_WORLD,ierr)
            call MPI_PACK(hist,660,mpiflt,
     +          buff,lenz,len1,MPI_COMM_WORLD,ierr)
            call MPI_PACK(path,606,mpiflt,
     +          buff,lenz,len1,MPI_COMM_WORLD,ierr)
c
c       characters
c
            call MPI_PACK(sp,8000,MPI_CHARACTER,
     +          buff,lenz,len1,MPI_COMM_WORLD,ierr)
c
c       the big send
c
            call MPI_SEND(buff,len0,MPI_PACKED,
     +          0,itag,MPI_COMM_WORLD,ierr)
c
c       master node
c
        else
c
c       the big receive
c
            idsrc=mpirgn-1
            call MPI_RECV(buff,len0,MPI_PACKED,
     +              idsrc,itag,MPI_COMM_WORLD,status,ierr)
c
c       integers
c
            lenz=len0
            len1=0
c
            call MPI_UNPACK(buff,lenz,len1,ip,100,
     +          mpiint,MPI_COMM_WORLD,ierr)
c
c       floats
c
            call MPI_UNPACK(buff,lenz,len1,rp,100,
     +          mpiflt,MPI_COMM_WORLD,ierr)
            call MPI_UNPACK(buff,lenz,len1,time,150,
     +          mpiflt,MPI_COMM_WORLD,ierr)
            call MPI_UNPACK(buff,lenz,len1,hist,660,
     +          mpiflt,MPI_COMM_WORLD,ierr)
            call MPI_UNPACK(buff,lenz,len1,path,606,
     +          mpiflt,MPI_COMM_WORLD,ierr)
c
c       characters
c
            call MPI_UNPACK(buff,lenz,len1,sp,8000,
     +          MPI_CHARACTER,MPI_COMM_WORLD,ierr)
        endif
c
        deallocate(buff)
        return
        end

