      subroutine pplotf(x,ix,b,lci,ct,ndf,ndm,nen1)

c     Purpose: Plot control subroutine for feap

c     Inputs:
c        x(ndm,*)   - Nodal coordinates for mesh
c        ix(nen1,*) - Element nodal connection lists
c        b(*)       - Current solution vecotr
c        lci        - Character plot option
c        ct(2)      - Parameters for plot
c        ndf        - Number dof/node
c        ndm        - Spatial dimension of mesh
c        nen1       - Dimension if ix array

c     Outputs:
c        none       - Graphics to plot area

      implicit  none

      logical   pcomp,oflg
      character lci*4
      integer*2 status,vslcol,coli
      integer   i,ic,n1, ndf,ndm,nen1
      real*8    c

      integer   ix(*)
      real*8    x(ndm,*),ct(2)
      real*8    b(*)

      real*8         dm
      common         dm(1000)

      real*8         dr
      integer              il
      common /adata/ dr(1),il(63998)

      integer        numnp,numel,nummat,nen,neq
      common /cdata/ numnp,numel,nummat,nen,neq

      integer*2       ixy
      common /pdata2/ ixy(4)

c     Open kernel system and plot mesh or outline of parts

      call pdevop()
      call frame(x,ndm,numnp)

c     Plot mesh or outline of parts

      oflg = .not.pcomp(lci,'outl')
      if(ct(1).ne.0.0.or.pcomp(lci,'eigv')) then
        ic = 2
        if(ct(1).eq.0.0) ct(1) = 1.0
      else
        ic = 1
      endif
      c  = 0.0
      n1 = 2*ndm*numnp
      do 100 i = ic,1,-1
        coli   = 8 - 2*i
        status = vslcol(coli)
        if(pcomp(lci,'eigv')) then
          call pdefm(x,dm,c,ndm,ndf,numnp, dr)
        else
          call pdefm(x, b,c,ndm,ndf,numnp, dr)
        endif
        call pline(dr,ix,il(n1),numnp,numel,ndm,nen1,nen,ct(2),oflg)
        c  = ct(1)
100   continue

c     Close plot

      call pdevcl()

      end

      subroutine pdevcl()

c     Purpose: Close the plotting device

      implicit  none

      integer*2 status,vclrwk,vencur,vclswk,vrqstr,ixy(2)

      character*1 xxx

      status = vrqstr(1,0,ixy,xxx)
      status = vclrwk()
      status = vencur()
      status = vclswk()

      end

      subroutine pdevop()

c     Purpose: Open graphics workstation

      implicit  none

      integer*2 status,vopnwk,vclrwk,vslcol
      integer*2 icl

c     Open kernel system

      status = vopnwk()
      if(status.lt.0) then
        write(*,2000)
        return
      endif
      status = vclrwk()
      icl    = 2
      status = vslcol(icl)
      call dplot(0.0000d0,0.0000d0,1)
      call dplot(1.4545d0,0.0000d0,2)
      call dplot(1.4545d0,1.0090d0,2)
      call dplot(0.0000d0,1.0090d0,2)
      call dplot(0.0000d0,0.0000d0,2)

2000  format('  Graphics Device Driver not installed correctly.')

      end

      subroutine plotl(x1,x2,x3,ipen)

c     Line drawing command

c     Inputs:
c        x1        - X-mesh coordinate to plot
c        x2        - Y-mesh coordinate to plot
c        x3        - Z-mesh coordinate to plot
c        ipen      - Pen action: = 2 draw; = 3 move

c     Outputs:
c        none

      implicit  none

      integer   ipen
      real*8    x1,x2,x3, s1,s2

      real*8          scale,dx   ,sx
      common /pdata1/ scale,dx(2),sx(2)

c     Compute the normal coordinates

      s1 = x3
      s1 = max(0.0,min(1.45,scale*(x1 + x1 - sx(1)) + 0.725))
      s2 = max(0.0,min(1.00,scale*(x2 + x2 - sx(2)) + 0.500))
      call dplot(s1,s2,ipen)

      end

      subroutine frame(x,ndm,numnp)

c     Purpose: Compute scaling for plot area

c     Inputs:
c        x(ndm,*)  - Nodal coordinates for mesh
c        ndm       - Spatial dimension of mesh
c        numnp     - Number of nodes in mesh

c     Outputs:
c        none      - Passed through common blocks

      implicit  none

      logical   iflg
      real*8    fact,xcen
      integer   i,ii,ij, n,ndm,numnp

      real*8    x(ndm,*),xmn(2),xmx(2),xmin(3),xmax(3)

      real*8          scale,dx   ,sx
      common /pdata1/ scale,dx(2),sx(2)

c     Determine window coordinates

      if(ndm.eq.1) then
        dx(2) = 0.
        sx(2) = 0.0
      endif
      ii = min(ndm,3)
      ij = min(ndm,2)

c     Find the minimum and maximum coordinate of input nodes

      iflg = .true.
      do 104 n = 1,numnp
        if(x(1,n).ne. -999.) then
          if(iflg) then
            do 100 i = 1,ii
              xmin(i) = x(i,n)
              xmax(i) = x(i,n)
100         continue
            iflg = .false.
          else
            do 102 i = 1,ii
              xmin(i) = min(xmin(i),x(i,n))
              xmax(i) = max(xmax(i),x(i,n))
102         continue
          endif
        endif
104   continue
      scale  = max(xmax(1)-xmin(1),xmax(2)-xmin(2))

c     Plot region determination

      do 110 i = 1,ij
        xmn(i) = min(xmin(i),xmax(i))
        xmx(i) = max(xmin(i),xmax(i))
        dx(i) = xmx(i) - xmn(i)
        sx(i) = xmx(i) + xmn(i)
110   continue

c     Rescale window

      if(dx(1).gt.1.45*dx(2)) then
        xmn(2) = (sx(2) - dx(1))/2.0
        xmx(2) = (sx(2) + dx(1))/2.0
        fact   = 0.58
      else
        xmn(1) = (sx(1) - dx(2))/2.0
        xmx(1) = (sx(1) + dx(2))/2.0
        fact   = 0.40
      endif
      do 112 i = 1,ij
        xmin(i) = max(xmin(i),xmn(i)) - scale/100.
        xmax(i) = min(xmax(i),xmx(i)) + scale/100.
112   continue

c     Default values

      scale  = max(xmax(1)-xmin(1),xmax(2)-xmin(2))

c     Reset values for deformed plotting

      do 114 i = 1,ij
        xcen = xmax(i)+xmin(i)
        xmax(i) = (xcen + 1.1*scale)/2.
        xmin(i) = (xcen - 1.1*scale)/2.
114   continue
      scale = fact/scale

      end

      subroutine pdefm(x,b,c,ndm,ndf,numnp, dr)

c     Purpose: Compute deformed position of two-dimensional meshes

c     Inputs:
c        x(ndm,*)  - Nodal coordinates of mesh
c        b(ndf,*)  - Solution vector
c        c         - Scale factor for added displacements
c        ndm       - Spatial dimension of mesh
c        ndf       - Number dof/node
c        numnp     - Number of nodes in mesh

c     Outputs:
c        dr(ndf,*) - Deformed nodal coordinates

      implicit  none

      real*8    c
      integer   i, n,ndf,ndm,numnp

      real*8 x(ndm,*),dr(ndm,*)
      real*8 b(ndf,*)

      call pconsd(dr,ndm*numnp,0.0d0)

      do 120 n = 1,numnp
        if(x(1,n).ne. -999.) then
          do 110 i = 1,ndm
            dr(i,n) = x(i,n) + c*b(i,n)
110       continue
        endif
120   continue

      end

      subroutine pline(x,ix,ic,numnp,numel,ndm,nen1,nen,ct,isw)

c     Purpose: Plot mesh or outline

c     Inputs:
c       x(ndm,*)   - Nodal coordinates of mesh
c       ix(nen1,*) - Element connection list
c       ic(numnp,*)- Working array
c       numnp      - Number of nodes in mesh
c       numel      - Number of elements in mesh
c       ndm        - Spatial dimension of mesh
c       nen1       - Dimension of ix array
c       nen        - Number of nodes/element, maximum
c       ct         - Color selection for plots
c       isw        - Flag, plot mesh if true, otherwise outline

c     Outputs:
c       none

      implicit  none

      logical   ifl,iend,isw
      real*8    ct,x3
      integer   i,ii,ij, j,jj, k, n,n1,n2,ni,nen,ndm,nen1,numnp,numel

      integer   ix(nen1,*),iplt(8),ic(numnp,*)
      real*8    x(ndm,*)

      data iplt/5,2,6,3,7,4,8,1/

c     Initialize connection array

      call pconsi(ic,numnp*4,0)

c     Loop through elements to set up list

      do 206 n = 1,numel
        jj = abs(ct)
        if(jj.eq.0) go to 197
        ii = abs(ix(nen1,n))
        if(ii.ne.jj) go to 206
197     i = 1
        ii = abs(ix(i,n))
        do 205 ij = 1,8
          j = iplt(ij)
          if(j.le.nen.and.ix(j,n).ne.0) then
            jj = abs(ix(j,n))
            if(jj.ne.ii) then
              n1 = min(ii,jj)
              n2 = max(ii,jj)
              do 203 k = 1,4
                if(ic(n1,k).eq.0) then
                  ic(n1,k) = n2
                  go to 204
                elseif(ic(n1,k).eq.n2) then
                  ic(n1,k) = -n2
                  go to 204
                endif
203           continue
            endif
204         ii = jj
          endif
205     continue
206   continue

c     change signs to permit mesh plot

      if(isw) then
        do 250 n = 1,numnp
        do 250 i = 1,4
          ic(n,i) = abs(ic(n,i))
250     continue
      endif

c     plot outline of part with continuous lines

      x3 = 0.
      do 304 ni = 1,numnp
        iend = .true.
        do 303 n = 1,numnp
          ifl = .true.
          n1 = n
300       do 301 i = 1,4
            if(ic(n1,i)) 301,303,302
301       continue
          go to 303
302       iend = .false.
          if(ndm.ge.3) x3 = x(3,n1)
          if(ifl) call plotl(x(1,n1),x(2,n1),x3,1)
          ifl = .false.
          n2 = ic(n1,i)
          ic(n1,i) = -n2
          if(ndm.ge.3) x3 = x(3,n2)
          call plotl(x(1,n2),x(2,n2),x3,2)
          n1 = n2
          go to 300
303     continue
        if(iend) go to 305
304   continue
305   return

      end
