PROGRAM inf_tied1
!
!      INFINITE SLOPE STABILITY ANALYSIS
!-------------------------------------------------------------------------
 USE main 
 USE geom   
 IMPLICIT NONE    
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 INTEGER::nels,nye,neq,nn,nip=4,nodof=2,nod=4,nst=4,ndof=8,i,j,k,iel,     &
   nosurf,iters,limit,iy,ndim=2,nprop=6,np_types,min_pos(1),nsrf=20,nlen
 LOGICAL::converged
 CHARACTER(LEN=80)::element='quadrilateral',trash,argv
 REAL(iwp)::e,v,det,phi,psi,gamma,f,dsbar,dq1,dq2,dq3,lode_theta,   &
   sigm,pi,tnph,phif,tnps,psif,snph,cf,tol=0.0001_iwp,h1,w1,s1,area,w,fmax, &
   gam_w,depth,k_h,srf,srf_lo,srf_hi,fos_tol,jump,srf_tmp,r_u, &
   zero=0.0_iwp,one=1.0_iwp,two=2.0_iwp,d4=4.0_iwp,d10=10.0_iwp,rtp,sigy,  &
   d180=180.0_iwp,pt999=0.999_iwp,pt5=0.5_iwp,oldsrf,beta,w_d
!---------------------------- dynamic arrays------------------------------
 REAL(iwp),ALLOCATABLE::kv(:),loads(:),points(:,:),bdylds(:),srf_h(:),    &
   evpt(:,:,:),oldis(:),gravlo(:),dee(:,:),coord(:,:),fun(:),jac(:,:),    &
   weights(:),der(:,:),deriv(:,:),bee(:,:),km(:,:),eld(:),eps(:),sigma(:),&
   bload(:),eload(:),erate(:),g_coord(:,:),evp(:),devp(:),m1(:,:),m2(:,:),&
   m3(:,:),flow(:,:),prop(:,:),surf(:,:),pore(:),gc(:),maxd_h(:),       &
   elastic(:),oldlds(:),srfm(:),ddt(:),xprop(:,:)
 INTEGER,ALLOCATABLE::nf(:,:),g(:),num(:),g_num(:,:),g_g(:,:),etype(:),   &
   kdiag(:),iters_h(:),gg(:)
!-------------------------input and initialisation------------------------
 CALL getname(argv,nlen)
 OPEN(10,FILE=argv(1:nlen)//'.dat')
 OPEN(11,FILE=argv(1:nlen)//'.res')
!
 READ(10,*)trash
 WRITE(11,'(a)')trash
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)h1
 WRITE(11,'(a)')"Height of soil column"
 WRITE(11,'(f8.2)')h1
 READ(10,*)trash
 READ(10,*)beta
 WRITE(11,'(/a)')"Slope angle"
 WRITE(11,'(f8.2)')beta
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)nye
 WRITE(11,'(a)')"Number of finite elements in column"
 WRITE(11,'(i5)')nye
 READ(10,*)trash
 READ(10,*)rtp
 READ(10,*)trash
 IF(rtp==1)THEN
   READ(10,*)np_types
 ELSE
   np_types=nye
 END IF
 WRITE(11,*)
 pi=d4*ATAN(one)
 w1=h1/10._iwp
 s1=w1*TAN(beta*pi/d180)
 nels=nye
 nn=2*(nye+1)
 neq=nye*4
 area= w1*h1/nye
 !-------------------------------------------------------------------------
 ALLOCATE(nf(nodof,nn),points(nip,ndim),weights(nip),g_coord(ndim,nn),    &
   num(nod),dee(nst,nst),evpt(nst,nip,nels),coord(nod,ndim),fun(nod),     &
   g_g(ndof,nels),jac(ndim,ndim),der(ndim,nod),etype(nels),               &
   deriv(ndim,nod),g_num(nod,nels),bee(nst,ndof),km(ndof,ndof),eld(ndof), &
   eps(nst),sigma(nst),bload(ndof),eload(ndof),erate(nst),evp(nst),       &
   devp(nst),g(ndof),m1(nst,nst),m2(nst,nst),m3(nst,nst),flow(nst,nst),   &
   prop(nprop,np_types),pore(nels),gc(ndim),gg(ndof),ddt(nels),       &
   surf(2,2),xprop(2,nprop))
!
! prop(1)=phi, prop(2)=c, prop(3)=psi, prop(4)=gamma, prop(5)=e, prop(6)=v
!
  READ(10,*)trash
  READ(10,*)trash
  IF(rtp==1)THEN
    DO i=1,np_types
      READ(10,*)prop(1:2,i),prop(4:6,i)
      prop(3,i)=prop(1,i)
    END DO
    WRITE(11,'(a)')"Group  phi      c      psi    gamma      e          v"
    DO i=1,np_types
     WRITE(11,'(i3,4f8.2,e12.4,f8.2)')i,prop(:,i)
    END DO
  ELSE
    READ(10,*)xprop(:,1:2),xprop(:,4:6)
    xprop(:,3)=xprop(:,1)
    WRITE(11,'(a)')"Prop     xprop(1,:)    xprop(2,:)"
    DO i=1,nprop
      WRITE(11,'(i3,2f13.2)')i,xprop(:,i)
    END DO
  ENDIF
 etype=1
 READ(10,*)trash
 IF(np_types>1.and.rtp==1)THEN
   READ(10,*)etype
 END IF
 WRITE(11,*)
 IF(np_types>1.AND.rtp==1)THEN
   WRITE(11,'(a)')"Property group assigned to each element"
   WRITE(11,'(20i3)')etype
 WRITE(11,*)
 END IF
 READ(10,*)trash
 READ(10,*)k_h
 WRITE(11,'(a)')"Pseudo horizontal acceleration"
 WRITE(11,'(f5.2)')k_h
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)nosurf
 !   0 means no water pressures
 !   1 means read r_u
 !   2 means read the vertical depth below ground level of a parallel free surface
 IF(nosurf==1)THEN
   READ(10,*)r_u
   WRITE(11,'(a)')"r_u parameter"
   WRITE(11,'(f5.2)')r_u
 END IF
 IF(nosurf==2)THEN
   READ(10,*)w_d
   surf(1,1)=zero
   surf(2,1)=-w_d
   surf(1,2)=w1
   surf(2,2)=-s1-w_d
   WRITE(11,'(a)')"Free-surface depth"
   WRITE(11,'(f5.2)')w_d
 END IF
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)gam_w
 READ(10,*)trash
 READ(10,*)limit
 WRITE(11,'(a)')"Iteration ceiling"
 WRITE(11,'(i5)')limit
 WRITE(11,*)
 READ(10,*)trash 
 READ(10,*)fos_tol
 IF(fos_tol>pt999)THEN
   nsrf=NINT(fos_tol)
   ALLOCATE(srfm(nsrf))
   READ(10,*)srfm
 ELSE
   WRITE(11,'(a)')"Factor of Safety tolerance"
   WRITE(11,'(f6.3)')fos_tol
 END IF
 ALLOCATE(loads(0:neq),bdylds(0:neq),oldis(0:neq),gravlo(0:neq),          &
   kdiag(neq),elastic(0:neq),oldlds(0:neq))
!---------- loop the element to set up global arrays --------------------
 kdiag=0
 nf=0
 IF(rtp/=1)THEN
 WRITE(11,*)
 WRITE(11,'(a)')"Element   phi      c      psi   gamma       e          v"
 ENDIF
 elements_1: DO iel=1,nels
   CALL inf_geometry_tied(iel,nye,h1,w1,s1,coord,num,g)
   g_num(:,iel)=num
   g_coord(:,num)=TRANSPOSE(coord)
   g_g(:,iel)=g
   IF(rtp/=1)THEN
    etype(iel)=iel
    prop(:,iel)=-xprop(1,:)*g_coord(2,g_num(1,iel))+xprop(2,:)
    WRITE(11,'(i5,4f8.2,e12.4,f8.2)')iel,prop(:,iel)
   ENDIF
   CALL fkdiag(kdiag,g)
   DO i=1,4
     nf(1,num(i))=g(2*i-1)
     nf(2,num(i))=g(2*i)
   END DO
 END DO elements_1
 DO i=2,neq
   kdiag(i)=kdiag(i)+kdiag(i-1)
 END DO
 ALLOCATE(kv(kdiag(neq)))
 kv=zero
 oldis=zero
 gravlo=zero
 pore=zero
 CALL sample(element,points,weights)
 pi=ACOS(-one)
 sigy=zero
!-gravity loads generation and element stiffness integration and assembly-
 elements_2: DO iel=1,nels
   gamma=prop(4,etype(iel))
   e=prop(5,etype(iel))
   v=prop(6,etype(iel))
   CALL deemat(dee,e,v)
   num=g_num(:,iel)
   coord=TRANSPOSE(g_coord(:,num)) 
   g=g_g(:,iel)
   km=zero
   w=area*gamma
   IF(iel==1)THEN
     gravlo(4*iel-2)=-w*pt5
     gravlo(4*iel)=  -w*pt5
     gravlo(4*iel-3)= w*pt5*k_h
     gravlo(4*iel-1)= w*pt5*k_h
   ELSE
     gravlo(4*iel-2)= gravlo(4*(iel-1)-2)-w*pt5
     gravlo(4*iel)=   gravlo(4*(iel-1))  -w*pt5
     gravlo(4*iel-3)= gravlo(4*(iel-1)-3)+w*pt5*k_h
     gravlo(4*iel-1)= gravlo(4*(iel-1)-1)+w*pt5*k_h
   END IF 
   gauss_pts_1: DO i=1,nip
     CALL shape_fun(fun,points,i)
     CALL shape_der(der,points,i)
     jac=MATMUL(der,coord) 
     det=determinant(jac)
     CALL invert(jac)
     deriv=MATMUL(jac,der)
     CALL beemat(bee,deriv)           
     km=km+MATMUL(MATMUL(TRANSPOSE(bee),dee),bee)*det*weights(i)
   END DO gauss_pts_1   
   IF(nosurf==1)THEN
     sigy=sigy+prop(4,etype(iel))*w1
     pore(iel)=-sigy*r_u
   ELSE IF(nosurf>1)THEN
     IF(w_d<zero)w_d=zero
     depth=w1*iel-w_d
     pore(iel)=-depth*gam_w*COS(ATAN(s1/w1))**2
   END IF
   CALL fsparv(kv,km,g,kdiag)
 END DO elements_2
 CALL sparin(kv,kdiag)                                                  
!------------------------trial strength reduction factor loop-------------
 oldsrf=d180
 jump=d10*fos_tol
 srf_lo=-jump
 srf_hi=zero
 srf=zero
 iters=0
 iy=0
!
 WRITE(11,*)
 WRITE(11,'(a)')"trial factor     max displ          iterations"
 ALLOCATE(srf_h(20),maxd_h(20),iters_h(20))
 srf_trials: do
   iy=iy+1
   IF(fos_tol>pt999)THEN
     srf=srfm(iy)
   ELSE
     IF(iters.NE.limit)THEN
       IF((srf-srf_lo)>pt999*jump)THEN
         srf=srf_hi+jump
         srf_lo=srf_hi
         srf_hi=srf
       ELSE
         srf_lo=srf
         srf=pt5*(srf_lo+srf_hi)
       END IF     
     ELSE 
       IF((srf-srf_lo)>pt999*jump)THEN
         srf=pt5*(srf_lo+srf_hi)
       ELSE
         srf_hi=srf
         srf=pt5*(srf_lo+srf_hi)
       END IF
     END IF
   END IF
   DO iel=1,nels
     phi=prop(1,etype(iel))
     tnph=TAN(phi*pi/d180)
     phif=ATAN(tnph/srf)
     snph=SIN(phif)
     e=prop(5,etype(iel))
     v=prop(6,etype(iel))
     ddt(iel)=d4*(one+v)*(one-two*v)/(e*(one-two*v+snph**2))
   END DO
   iters=0
   bdylds=zero
   evpt=zero
 nip=1
 DEALLOCATE(points,weights)
 ALLOCATE(points(nip,ndim),weights(nip))
 CALL sample(element,points,weights)
 !----------------------------iteration loop  -----------------------------
   its: DO    
     fmax=zero
     iters=iters+1
     loads=gravlo+bdylds
     CALL spabac(kv,loads,kdiag)
     loads(0)=zero
     IF(iy==1.and.iters==1)elastic=loads
!-------------------------check convergence ------------------------------
     CALL checon1(loads,oldis,tol,converged)
     IF(iters==1)converged=.FALSE.
     IF(converged.OR.iters==limit)bdylds=zero
!----------------------- go round the Gauss Points  ----------------------
     elements_3: DO iel=1,nels
       phi=prop(1,etype(iel))
       tnph=TAN(phi*pi/d180)
       phif=ATAN(tnph/srf)*d180/pi
       psi=prop(3,etype(iel))
       tnps=TAN(psi*pi/d180)
       psif=ATAN(tnps/srf)*d180/pi
       cf=prop(2,etype(iel))/srf
       e=prop(5,etype(iel))
       v=prop(6,etype(iel))
       CALL deemat(dee,e,v)
       num=g_num(:,iel)
       coord=TRANSPOSE(g_coord(:,num))
       g=g_g(:,iel)
       eld=loads(g)
       bload=zero
       gauss_points_2: DO i=1,nip
         CALL shape_der(der,points,i)
         jac=MATMUL(der,coord)
         det=determinant(jac)
         CALL invert(jac)
         deriv=MATMUL(jac,der) 
         CALL beemat(bee,deriv)
         eps=MATMUL(bee,eld)
         eps=eps-evpt(:,i,iel)
         sigma=MATMUL(dee,eps)
         sigma(1)=sigma(1)-pore(iel)
         sigma(2)=sigma(2)-pore(iel)
         sigma(4)=sigma(4)-pore(iel)
         CALL invar(sigma,sigm,dsbar,lode_theta)                             
!------------------  check whether yield is violated  --------------------
         CALL mocouf(phif,cf,sigm,dsbar,lode_theta,f)
         IF(converged.OR.iters==limit)THEN
           sigma(1)=sigma(1)+pore(iel)
           sigma(2)=sigma(2)+pore(iel)
           sigma(4)=sigma(4)+pore(iel)
           devp=sigma 
           IF(f>fmax)fmax=f
         else
           IF(f>=zero)THEN
!            write(11,'(f10.4,i5,f10.4)')srf,iel,f
            IF(f>fmax)fmax=f
             CALL mocouq(psif,dsbar,lode_theta,dq1,dq2,dq3)
             CALL formm(sigma,m1,m2,m3)
             flow=f*(m1*dq1+m2*dq2+m3*dq3) 
             erate=MATMUL(flow,sigma)
             evp=erate*ddt(iel)
             evpt(:,i,iel)=evpt(:,i,iel)+evp
             devp=MATMUL(dee,evp) 
           END IF
         END IF
         IF(f>=zero.OR.(converged.OR.iters==limit))THEN
           eload=MATMUL(devp,bee)
           bload=bload+eload*det*weights(i)
         END IF
       END DO gauss_points_2
!-------------------  compute the total bodyloads vector -----------------
       bdylds(g)=bdylds(g)+bload
       bdylds(0)=zero
     END DO elements_3             
     WRITE(*,'(a,f9.6,a,i4,a,f10.2)')                                      &
       "   SRF=",srf,"   ITERS=",iters,"   FMAX=",fmax    
     IF(iters==limit)THEN
       IF(srf<oldsrf)THEN
         oldlds=loads
         oldsrf=srf
       END IF
       EXIT
     END IF
     IF(converged)EXIT
   END DO its
   WRITE(*,*)
   srf_h(iy)=srf
   DO iel=nels,2,-1
    g=g_g(:,iel)
    gg=g_g(:,iel-1)
    loads(gg)=loads(gg)+loads(g)
   ENDDO                                                            
   maxd_h(iy)=MAXVAL(abs(loads))
   iters_h(iy)=iters
   IF(fos_tol>pt999)THEN
     WRITE(11,'(f10.4,e16.5,i15)')srf_h(iy),maxd_h(iy),iters_h(iy)
     IF(iy==nsrf)EXIT
   ELSE
     IF((srf_hi-srf_lo<fos_tol).OR.iy==nsrf)EXIT
   END IF
 END DO srf_trials
!
 IF(fos_tol<pt999)THEN
   DO i=1,iy
     min_pos=MINLOC(srf_h(i:iy))
     srf_tmp=srf_h(i+min_pos(1)-1)
     srf_h(i+min_pos(1)-1)=srf_h(i)
     srf_h(i)=srf_tmp
     srf_tmp=maxd_h(i+min_pos(1)-1)
     maxd_h(i+min_pos(1)-1)=maxd_h(i)
     maxd_h(i)=srf_tmp
     j=iters_h(i+min_pos(1)-1)
     iters_h(i+min_pos(1)-1)=iters_h(i)
     iters_h(i)=j
     IF(iters_h(i)<limit)                                                   &
       WRITE(11,'(f10.4,e16.5,i15)')srf_h(i),maxd_h(i),iters_h(i)
     IF(i>1.AND.iters_h(i)==limit.AND.iters_h(i-1)<limit)k=i
   END DO
   WRITE(11,'(f10.4,e16.5,i15)')srf_h(k),maxd_h(k),iters_h(k)
   WRITE(11,*)
   WRITE(11,'(a)')"Estimated Factor of Safety"
   WRITE(11,'(f8.2)')srf_h(k)
 END IF
!
 oldlds=oldlds-elastic 
 DO iel=nels,2,-1
 g=g_g(:,iel)
 gg=g_g(:,iel-1)
 oldlds(gg)=oldlds(gg)+oldlds(g)
 ENDDO                                                            
 CALL mesh(g_coord,g_num,argv,nlen,12)
 CALL dismsh(loads,nf,0.2_iwp,g_coord,g_num,argv,nlen,13)
 CALL vecmsh(loads,nf,0.05_iwp,0.1_iwp,g_coord,g_num,argv,nlen,14)
!
STOP
CONTAINS
!
   SUBROUTINE inf_geometry_tied(iel,nels,h1,w1,s1,coord,num,g)
!
!      this subroutine forms the coordinates and steering vector
!      for 4-node quadrilaterals in a 'infintite slope' geometry
!      (elements and nodes numbered in the y-direction)
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::h1,w1,s1
 INTEGER,INTENT(IN)::iel,nels
 REAL(iwp),INTENT(OUT)::coord(:,:) 
 REAL(iwp)::zero=0.0_iwp
 INTEGER,INTENT(OUT)::num(:),g(:)
 INTEGER::iq
 iq=iel 
 num(1)=iq*2+1
 num(2)=num(1)-2
 num(3)=num(2)+1
 num(4)=num(1)+1
 g(1)=0
 g(2)=0
 g(3)=4*(iq-1)+1
 g(4)=g(3)+1
 g(5)=g(4)+1
 g(6)=g(5)+1
 g(7)=0
 g(8)=0
 coord(1,1)=zero
 coord(2,1)=zero
 coord(3,1)=w1
 coord(4,1)=w1
 coord(1,2)=    -iq*h1/nels
 coord(2,2)=-(iq-1)*h1/nels
 coord(3,2)=coord(2,2)-s1
 coord(4,2)=coord(1,2)-s1
RETURN
END SUBROUTINE inf_geometry_tied
!
SUBROUTINE vecs(loads,nf,ratmax,cutoff,g_coord,g_num,ips)
!
! This subroutine produces a PostScript output file "*.vec" displaying
! the nodal displacement vectors.
!
IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::g_coord(:,:),loads(0:),ratmax,cutoff
 INTEGER,INTENT(IN)::g_num(:,:),ips,nf(:,:)
 REAL(iwp)::width,height,scale=72,sxy,xo,yo,x1,y1,x2,y2,dismag,           &
   zero=0.0_iwp,pt5=0.5_iwp,opt5=1.5_iwp,fpt5=5.5_iwp,d8=8.0_iwp,         &
   ept5=8.5_iwp,d11=11.0_iwp,xmin,xmax,ymin,ymax,dmax,vlen,vmax
 INTEGER::i,j,k,l,nn,nels,nod,ns,i1,i2,j1,j2
 INTEGER,ALLOCATABLE::corner(:,:)
 LOGICAL::draw
!                       formats
 OPEN(ips,FILE='fe95.vec')
!                       open output file and compute scale factors
 nn=UBOUND(nf,2)
!
 xmin=g_coord(1,1)
 xmax=g_coord(1,1)
 ymin=g_coord(2,1)
 ymax=g_coord(2,1)
 DO i=2,nn
   IF(g_coord(1,i)<xmin)xmin=g_coord(1,i)      
   IF(g_coord(1,i)>xmax)xmax=g_coord(1,i)      
   IF(g_coord(2,i)<ymin)ymin=g_coord(2,i)      
   IF(g_coord(2,i)>ymax)ymax=g_coord(2,i)      
 END DO
 width =xmax-xmin
 height=ymax-ymin
 dmax=ratmax*width
 IF(height>width)dmax=ratmax*height
!
 vmax=zero
 DO i=1,nn
   DO j=1,2
     IF(ABS(loads(nf(j,i)))>vmax)vmax=ABS(loads(nf(j,i)))
   END DO
 END DO
 dismag=dmax/vmax
!
 xmin=g_coord(1,1)
 xmax=g_coord(1,1)
 ymin=g_coord(2,1)
 ymax=g_coord(2,1)
!
 DO i=1,nn
   IF(g_coord(1,i)+dismag*loads(nf(1,i)) < xmin)                          &
     xmin=g_coord(1,i)+dismag*loads(nf(1,i))      
   IF(g_coord(1,i)+dismag*loads(nf(1,i)) > xmax)                          &
     xmax=g_coord(1,i)+dismag*loads(nf(1,i))      
   IF(g_coord(2,i)+dismag*loads(nf(2,i)) < ymin)                          &
     ymin=g_coord(2,i)+dismag*loads(nf(2,i))      
   IF(g_coord(2,i)+dismag*loads(nf(2,i)) > ymax)                          &
     ymax=g_coord(2,i)+dismag*loads(nf(2,i))      
!
   IF(g_coord(1,i)<xmin)xmin=g_coord(1,i)
   IF(g_coord(1,i)>xmax)xmax=g_coord(1,i)
   IF(g_coord(2,i)<ymin)ymin=g_coord(2,i)
   IF(g_coord(2,i)>ymax)ymax=g_coord(2,i)
 END DO
!
 width=xmax-xmin
 height=ymax-ymin
!
!                       allow 1.5" margin minimum on each side of figure
!
!                       Portrait mode 
!
 IF(height.GE.d11/ept5*width)THEN
!
!                       height governs the scale
!
   sxy=scale*d8/height
   xo=scale*pt5*(ept5-d8*width/height)
   yo=scale*opt5
 ELSE
!
!                       width governs the scale
!
   sxy=scale*fpt5/width
   xo=scale*opt5
   yo=scale*pt5*(d11-fpt5*height/width)
 END IF
!
 WRITE(ips,'(a)')'%!PS-Adobe-1.0'
 WRITE(ips,'(a)')'%%DocumentFonts: none'
 WRITE(ips,'(a)')'%%Pages: 1'
 WRITE(ips,'(a)')'%%EndComments'
 WRITE(ips,'(a)')'/m {moveto} def'
 WRITE(ips,'(a)')'/l {lineto} def'
 WRITE(ips,'(a)')'/s {stroke} def'
 WRITE(ips,'(a)')'/c {closepath} def'
 WRITE(ips,'(a)')'/edef {exch def} bind def'
 WRITE(ips,'(a)')                                                         &
   '/arrow {/@to_y edef /@to_x edef /@from_y edef /@from_x edef'
 WRITE(ips,'(a)')'/@dx @to_x @from_x sub def /@dy @to_y @from_y sub def'
 WRITE(ips,'(a)')'/@length @dx @dx mul @dy @dy mul add sqrt def'
 WRITE(ips,'(a)')'/@angle @dy @dx atan def'
 WRITE(ips,'(a)')'gsave @from_x @from_y translate @angle rotate'
 WRITE(ips,'(a)')                                                         &
   '0 0 moveto @length 0 lineto currentpoint stroke newpath moveto'
 WRITE(ips,'(a)')'-4 -2 rlineto @length 0 moveto'
 WRITE(ips,'(a)')'-4  2 rlineto stroke grestore'
 WRITE(ips,'(a)')'} def'
 WRITE(ips,'(a)')'/*sf {'
 WRITE(ips,'(a)')'exch findfont exch'
 WRITE(ips,'(a)')                                                         &
   'dup type /arraytype eq {makefont}{scalefont} ifelse setfont'
 WRITE(ips,'(a)')'} bind def'
 WRITE(ips,'(a)')'/languagelevel where'
 WRITE(ips,'(a)')'{pop languagelevel} {1} ifelse'
 WRITE(ips,'(a)')'2 lt { % ifelse'
 WRITE(ips,'(a)')'/sf /*sf load def'
 WRITE(ips,'(a)')'} { % else'
 WRITE(ips,'(a)')'/sf /selectfont load def'
 WRITE(ips,'(a)')'} ifelse'
 WRITE(ips,'(a)')'%%EndProlog'
 WRITE(ips,'(a)')'%%Page: 0 1'
 WRITE(ips,'(a)')'gsave'
!
 WRITE(ips,'(2f9.2,a)')xo,yo, ' translate'
 WRITE(ips,'(f9.2,a)')0.5,' setlinewidth'
!
!                       draw the displacement vectors
!
 vmax=zero
 DO i=1,nn
   vlen=loads(nf(1,i))**2+loads(nf(2,i))**2
   IF(vlen>vmax)vmax=vlen
 END DO
 vmax=SQRT(vmax)*cutoff      
 DO i=1,nn
   vlen=SQRT(loads(nf(1,i))**2+loads(nf(2,i))**2)
   x1=sxy*(g_coord(1,i)-xmin)
   y1=sxy*(g_coord(2,i)-ymin)
   x2=sxy*(g_coord(1,i)+dismag*loads(nf(1,i))-xmin)
   y2=sxy*(g_coord(2,i)+dismag*loads(nf(2,i))-ymin)
   IF(vlen>vmax)THEN
     WRITE(ips,'(2f9.2,a,2f9.2,a)') x1, y1,' ', x2, y2, ' arrow'
     WRITE(ips,'(a)') 's'
   END IF
 END DO
!
!                       draw the mesh border                          
!
 nels=UBOUND(g_num,2)
 nod=UBOUND(g_num,1)
 IF(nod==3.OR.nod==6.OR.nod==10.OR.nod==15)ns=3
 IF(nod==4.OR.nod==8.OR.nod==9)ns=4
 ALLOCATE(corner(ns,2))
 IF(nod== 3)corner=RESHAPE((/1,2,3,2,3,1/),(/3,2/))
 IF(nod== 6)corner=RESHAPE((/1,3,5,3,5,1/),(/3,2/))
 IF(nod==10)corner=RESHAPE((/1,4,7,4,7,1/),(/3,2/))
 IF(nod==15)corner=RESHAPE((/1,5,9,5,9,1/),(/3,2/))
 IF(nod== 4)corner=RESHAPE((/1,2,3,4,2,3,4,1/),(/4,2/))
 IF(nod== 8)corner=RESHAPE((/1,3,5,7,3,5,7,1/),(/4,2/))
 IF(nod== 9)corner=RESHAPE((/1,3,5,7,3,5,7,1/),(/4,2/))
 DO i=1,nels
   DO j=1,ns
     draw=.TRUE.
     i1=g_num(corner(j,1),i)
     i2=g_num(corner(j,2),i)
     DO k=1,nels
       DO l=1,ns
         j1=g_num(corner(l,1),k)
         j2=g_num(corner(l,2),k)
         IF((i1==j2).AND.(i2==j1))THEN
           draw=.FALSE.
           EXIT
         END IF
       END DO
       IF(.NOT.draw)EXIT
     END DO
     IF(draw)THEN
       x1=sxy*(g_coord(1,i1)-xmin)
       y1=sxy*(g_coord(2,i1)-ymin)
       WRITE(ips,'(2f9.2,a)')x1, y1,' m'
       x1=sxy*(g_coord(1,i2)-xmin)
       y1=sxy*(g_coord(2,i2)-ymin)
       WRITE(ips,'(2f9.2,a)')x1, y1,' l'
       WRITE(ips,'(a)')' s'
     END IF
   END DO
 END DO
!                       close output file?
 WRITE(ips,'(a)')'grestore'
 WRITE(ips,'(a)')'showpage'
 CLOSE(ips)
RETURN
END SUBROUTINE vecs
SUBROUTINE checon1(loads,oldlds,tol,converged)
!
! This subroutine sets converged to .FALSE. if relative change in loads
! and oldlds is greater than tol.
!
 IMPLICIT NONE
 INTEGER::i,neq
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::loads(0:),tol
 REAL(iwp),INTENT(IN OUT)::oldlds(0:)
 LOGICAL,INTENT(OUT)::converged
  converged=.TRUE.
 neq=UBOUND(loads,1)
 DO i=1,neq
  IF(ABS(loads(i)-oldlds(i))/ABS(oldlds(i))>=tol)THEN
  converged=.FALSE.
  EXIT
  ENDIF
 ENDDO
 oldlds=loads
RETURN
END SUBROUTINE checon1
END PROGRAM inf_tied1





















