cc****************************************************************************
cc****************************************************************************
cc                                                                          cc
cc                                                                          cc
cc                               R E T P L A                                cc
cc                                                                          cc
cc            Programa de calcul d'estructures reticulades planes           cc
cc                                                                          cc
cc                                                                          cc
c*****************************************************************************
c*****************************************************************************
      implicit real*8(a-h,o-z)
c
c***  Dimensions maximes
c
      parameter  ( maxpoin= 100, maxelem=  50, maxmats=   4,
     .             maxvfix=  20, maxband=  26, maxtotv=3*maxpoin )
c
c***  Dimensionament de les variables principals i significat
c
      integer*4  ifpre(3*maxvfix),        ! codis de vinculacio als suports
     .           lnods(2*maxelem),        ! connectivitats de les barres
     .           matno(  maxelem),        ! tipus de barra
     .           nofix(  maxvfix)         ! numeracio global dels suports
      real*8     coord(2*maxpoin),        ! coordenades nodal 
     .           gload(3*maxpoin),        ! carregues exteriors
     .           presc(3*maxvfix),        ! valor de les condicions de contorn
     .           props(3*maxmats),        ! propietats dels materials (E,A,I)
     .           rigid(maxband*maxtotv),  ! matriu de rigidesa global
     .           tdisp(3*maxpoin),        ! vector de moviments global
     .           treac(3*maxvfix),        ! reaccions als suports
     .                                    ! per a cada barra:
     .           estif(36),               ! matriu de rigidesa global (Kg)
     .           estil(36),               ! matriu de rigidesa local (Kl)
     .           etmat(36),               ! matriu per calcular esforcos (Kl.T)
     .           tmatr(9)                 ! matriu de canvi de base (T)
c
c***  Blocs COMMON que empra el programa
c
      character c_coor*12,c_prop*14,c_dofn*9,c_load*13,c_stre*10
      common /arxius/ n_lect,n_escr,n_kltt,n_tran
      common /contrl/ npoin,nelem,nmats,nvfix,nload,nband,ntotv
      common /titols/ c_coor(2),c_prop(3),c_dofn(3),c_load(3),c_stre(3)
c
c***  Programa principal
c
      call obrarx                              ! Obre arxius
      call lecdad (coord,props,presc,gload,    ! Llegeix les dades
     .             lnods,matno,nofix,ifpre)
      call prevop (lnods,rigid,treac,maxband)  ! Algunes operacions previes
      call outdat (coord,props,presc,gload,    ! Escriu les dades
     .             lnods,matno,nofix,ifpre)
      call genrig (coord,props,rigid,estif,    ! Genera la matriu de rigidesa
     .             etmat,estil,tmatr,lnods,    ! de l'estructura
     .             matno,nofix)
      call vincul (rigid,gload,treac,nofix,    ! Imposa les condicions de
     .             ifpre,presc)                ! vinculacio
      call solsis (ntotv,nband,rigid,rigid,    ! Calcula els moviments
     .                 1,gload,tdisp)
      call outdes (tdisp)                      ! Escriu els moviments
      call calesf (tdisp,treac,etmat,tmatr,    ! Calcula i escriu reaccions i
     .             lnods,matno,nofix,ifpre)    ! esforcos d'extrem de barra
      call tanarx                              ! Tanca arxius

      stop
      end
      subroutine obrarx
c*****************************************************************************
c
c   Aquesta subrutina obre els arxius d'entrada i sortida de dades
c
c*****************************************************************************
      implicit real*8(a-h,o-z)
      common /arxius/ n_lect,n_escr,n_kltt,n_tran
      character*40 arx_lect,arx_escr,titol
c
c***  Obre arxius d'entrada i sortida de dades
c
      n_lect=20
 2000 write(6,'(a,$)') ' >>> Nom de l''arxiu de dades :'
      read(5,'(a)') arx_lect
      open(unit=n_lect,err=2000,file=arx_lect,status='old')

      n_escr=21
 2010 write(6,'(a,$)') ' >>> Nom de l''arxiu de resultats :'
      read(5,'(a)') arx_escr
      open(unit=n_escr,err=2010,file=arx_escr,status='unknown')
c
c***  Escriu el titol del problema
c
      write(6,'(a,$)') ' >>> Titol del problema (a40) :'
      read(5,'(a)') titol
      write(n_escr,10) titol
  10  format(///,10x,60('*'),//,20x,a40,//,10x,60('*'),//)
c
c***  Obre arxius temporals d'emmagatzematge de les matrius Kl.Tt i T
c
      n_kltt=10
      n_tran=11
      open(unit=n_kltt,form='unformatted',status='scratch')
      open(unit=n_tran,form='unformatted',status='scratch')

      return
      end
      subroutine lecdad(coord,props,presc,gload,lnods,matno,nofix,ifpre)
c*****************************************************************************
c
c   Lectura de les dades del problema
c
c*****************************************************************************
      implicit real*8(a-h,o-z)
      common /contrl/ npoin,nelem,nmats,nvfix,nload,nband,ntotv
      common /arxius/ n_lect,n_escr,n_kltt,n_tran
      real*8     coord(*), props(*), presc(*), gload(*)
      integer*4  lnods(*), matno(*), nofix(*), ifpre(*)
c
c***  Lectura de les variables de control
c
      read(n_lect,10) npoin,nelem,nmats,nvfix,nload
  10  format(/,5i5)
      ntotv=3*npoin
c
c***  Lectura de la resta de dades
c
      call auxdad (coord,props,presc,gload,lnods,matno,nofix,ifpre)

      return
      end
      subroutine auxdad(coord,props,presc,gload,lnods,matno,nofix,ifpre)
c*****************************************************************************
c
c   Llegeix la geometria, propietats de les barres, condicions de contorn
c   i carregues nodals
c
c*****************************************************************************
      implicit real*8(a-h,o-z)
      common /contrl/ npoin,nelem,nmats,nvfix,nload,nband,ntotv
      common /arxius/ n_lect,n_escr,n_kltt,n_tran
      real*8     coord(2,npoin), props(3,nmats), presc(3,nvfix),
     .           gload(ntotv)
      integer*4  lnods(2,nelem), matno(nelem), nofix(nvfix),
     .           ifpre(3,nvfix)
c
c***  Lectura de la matriu de connectivitats i tipus de barra
c
      do ielem=1,nelem
        read(n_lect,10) jelem,matno(jelem),(lnods(i,jelem),i=1,2)
      end do
c
c***  Lectura de les coordenades nodals
c
      do ipoin=1,npoin
        read(n_lect,11) jpoin,(coord(i,jpoin),i=1,2)
      end do
c
c***  Lectura de les propietats de les barres
c
      do imats=1,nmats
        read(n_lect,12) jmats,(props(i,jmats),i=1,3)
      end do
c
c***  Lectura de les condicions de vinculacio
c
      do ivfix=1,nvfix
        read(n_lect,13) jvfix,nofix(jvfix),(ifpre(i,jvfix),i=1,3),
     .                                     (presc(i,jvfix),i=1,3)
      end do
c
c***  Inicialitzacio del vector de carrega nodal 
c
      do itotv=1,ntotv
        gload(itotv)=0.0
      end do
c
c***  Lectura de les carregues nodals
c
      do iload=1,nload
        read(n_lect,14) jpoin, (gload(3*(jpoin-1)+i),i=1,3)
      end do
c
c***  Formats de lectura
c
  10  format(4i5)                              ! connectivitats
  11  format(i5,2f10.5)                        ! coordenades
  12  format(i5,3e15.6)                        ! propietats de les barres
  13  format(2i5,3i2,3f10.5)                   ! condicions de vinculacio
  14  format(i5,3e15.6)                        ! carregues nodals

      return
      end
      subroutine prevop(lnods,rigid,treac,maxband)
c*****************************************************************************
c
c   Calcula el semiample de banda, inicialitza vectors i defineix variables
c
c*****************************************************************************
      implicit real*8(a-h,o-z)
      character c_coor*12,c_prop*14,c_dofn*9,c_load*13,c_stre*10
      common /contrl/ npoin,nelem,nmats,nvfix,nload,nband,ntotv
      common /arxius/ n_lect,n_escr,n_kltt,n_tran
      common /titols/ c_coor(2),c_prop(3),c_dofn(3),c_load(3),c_stre(3)
      real*8     rigid(nband+1,ntotv), treac(3,nvfix)
      integer*4  lnods(2,nelem)
c
c***  Calcula el semiample de banda
c
      nband=0
      do ielem=1,nelem
        iband=abs(lnods(1,ielem)-lnods(2,ielem))
        if (iband.gt.nband) nband=iband
      end do
      nband=(nband+1)*3-1
      if(nband.ge.maxband) then
        write(6,'(a,i5,a)') 
     .  ' >>> Cal augmentar MAXBAND a ',nband+1,' !!!'
        stop
      end if
c
c***  Inicialitzacio de les matrius de treball
c
      do iband=1,nband+1
        do itotv=1,ntotv
          rigid(iband,itotv)=0.0
        end do
      end do

      do i=1,3
        do ivfix=1,nvfix
          treac(i,ivfix)=0.0
        end do
      end do
c
c***  Defineix algunes variables caracter
c
      c_coor(1)='Coordenada-X'
      c_coor(2)='Coordenada-Y'
      c_prop(1)='Mod. Young (E)'
      c_prop(2)='      Area (A)'
      c_prop(3)='   Inercia (I)'
      c_dofn(1)='Despl. -X'
      c_dofn(2)='Despl. -Y'
      c_dofn(3)='      Gir'
      c_load(1)='Forca Ext. -X'
      c_load(2)='Forca Ext. -Y'
      c_load(3)='  Moment Ext.'
      c_stre(1)='  E. axial'
      c_stre(2)='E. tallant'
      c_stre(3)='    Moment'

      return
      end
      subroutine outdat(coord,props,presc,gload,lnods,matno,nofix,ifpre)
c*****************************************************************************
c
c   Sortida de les dades del problema a l'arxiu de resultats
c
c*****************************************************************************
      implicit real*8(a-h,o-z)
      character c_coor*12,c_prop*14,c_dofn*9,c_load*13,c_stre*10
      common /contrl/ npoin,nelem,nmats,nvfix,nload,nband,ntotv
      common /arxius/ n_lect,n_escr,n_kltt,n_tran
      common /titols/ c_coor(2),c_prop(3),c_dofn(3),c_load(3),c_stre(3)
      real*8     coord(2,npoin), props(3,nmats), presc(3,nvfix),
     .           gload(ntotv)
      integer*4  lnods(2,nelem), matno(nelem), nofix(nvfix),
     .           ifpre(3,nvfix)
c
c***  Escriptura de les variables de control 
c
      write(n_escr,100) npoin,nelem,nmats,nvfix,nload,nband,ntotv
  100 format(//,' VARIABLES DE CONTROL DEL PROGRAMA',/,
     .          ' --------- -- ------- --- --------',//,
     .    5x,' Num. de nodes             :',i5,/,
     .    5x,' Num. de barres            :',i5,/,
     .    5x,' Num. de tipus barra       :',i5,/,
     .    5x,' Num. de nodes coaccionats :',i5,/,
     .    5x,' Num. de nodes carregats   :',i5,/,
     .    5x,' Semiample de banda de K   :',i5,/,
     .    5x,' Num. de g.d.l total       :',i5)
c
c***  Escriptura de la matriu de connectivitats
c
      write(n_escr,101)
  101 format(//' MATRIU DE CONNECTIVITAT',/,
     .         ' ------ -- -------------',//,
     .         '          Barra',
     .         '          Tipus',
     .         '    Node dorsal',
     .         '   Node frontal')
      do ielem=1,nelem
        write(n_escr,10) ielem,matno(ielem),(lnods(i,ielem),i=1,2)
      end do
c
c***  Escriptura de les coordenades nodals
c
      write(n_escr,102) (c_coor(i),i=1,2)
  102 format(//,' COORDENADES NODALS',/,
     .          ' ----------- ------',//,
     .          '        Num. de node',2(8x,a12))
      do ipoin=1,npoin
        write(n_escr,11) ipoin,(coord(i,ipoin),i=1,2)
      end do
c
c***  Escriptura de les propietats de les barres
c
      write(n_escr,103) (c_prop(i),i=1,3)
  103 format(//,' PROPIETATS DE LES BARRES',/,
     .          ' ---------- -- --- ------',//,
     .          ' Tipus de barra num.',3(2x,a14))
      do imats=1,nmats
        write(n_escr,12) imats,(props(i,imats),i=1,3)
      end do
c
c***  Escriptura de les condicions de vinculacio
c
      write(n_escr,104) (c_dofn(idofn),idofn=1,3)
  104 format(//,' CONDICIONS DE VINCULACIO',/,
     .          ' ---------- -- ----------',//,
     .          ' Node coaccionat','  Codis',3(5x,a9))
      do ivfix=1,nvfix
        write(n_escr,13) nofix(ivfix),(ifpre(i,ivfix),i=1,3),
     .                                (presc(i,ivfix),i=1,3)
      end do
c
c***  Escriptura de les carregues nodals
c
      write(n_escr,105) (c_load(i),i=1,3)
  105 format(//,' CARREGUES NODALS',/,
     .          ' --------- ------',//,
     .          ' Node carregat',3(2x,a13)) 
      do ipoin=1,npoin
        itotv=3*(ipoin-1)
        fx=gload(itotv+1) 
        fy=gload(itotv+2) 
        zm=gload(itotv+3) 
        if(fx.ne.0.0 .or. fy.ne.0.0 .or. zm.ne.0.0)
     .  write(n_escr,14) ipoin,fx,fy,zm
      end do
c
c***  Formats
c
  10  format(4i15)     
  11  format(i20,2f20.6)
  12  format(i20,3e16.6)
  13  format(i16,4x,3i1,3(2x,f12.6))
  14  format(i14,3(e15.6))

      return
      end
      subroutine genrig(coord,props,rigid,estif,etmat,estil,tmatr,
     .                  lnods,matno,nofix)
c*****************************************************************************
c
c   Aquesta rutina genera la matriu de rigidesa per a tota l'estructura
c
c*****************************************************************************
      implicit real*8(a-h,o-z)
      common /contrl/ npoin,nelem,nmats,nvfix,nload,nband,ntotv
      common /arxius/ n_lect,n_escr,n_kltt,n_tran
      real*8     coord(2,npoin), props(3,nmats), rigid(nband+1,ntotv),
     .           estif(6,6), etmat(6,6), estil(6,6), tmatr(3,3)
      integer*4  lnods(2,nelem), matno(nelem), nofix(nvfix)
c
c***  Llac sobre el nombre de barres
c
      do ielem=1,nelem
c
c***  Calcul de la matriu de rigidesa de cada barra
c
        node1=lnods(1,ielem)
        node2=lnods(2,ielem)
        call rigbar (coord,props(1,matno(ielem)),tmatr,estil,etmat,
     .               estif,node1,node2)
c
c***  Emmagatzematge de la matriu Kl.Tt per determinar esforcos
c
        write(n_kltt) ((etmat(i,j),i=1,6),j=1,6)
c
c***  Si la barra te algun node prescrit, es guarda la matriu de canvi i
c***  es canvia el signe de MATNO per identificar-la
c
        ivfix=0
        do while(ivfix.le.nvfix .and. matno(ielem).gt.0)
          ivfix=ivfix+1
          ipoin=nofix(ivfix)
          if(node1.eq.ipoin .or. node2.eq.ipoin) then
            matno(ielem)=-matno(ielem)
            write(n_tran) ((tmatr(i,j),i=1,3),j=1,3)
          end if
        end do
c
c***  Assembly de la matriu de la barra en la matriu de l'estructura
c
        call assemb (ielem,lnods(1,ielem),estif,rigid)
      end do

      return
      end
      subroutine rigbar(coord,props,tmatr,estil,etmat,estif,node1,node2)
c*****************************************************************************
c
c   Calcula la matriu de rigidesa de cada barra en coordenades locals (estil),
c   locals-globals (etmat) i globals (estif)
c
c*****************************************************************************
      implicit real*8(a-h,o-z)
      common /contrl/ npoin,nelem,nmats,nvfix,nload,nband,ntotv
      real*8    coord(2,npoin), props(3),   tmatr(3,3), 
     .          estil(6,6),     etmat(6,6), estif(6,6)
c
c***  Inicialitzacio de les matrius de rigidesa i de la matriu de canvi
c
      do i=1,6
        do j=1,6
          estil(i,j)=0.0
          etmat(i,j)=0.0
          estif(i,j)=0.0
        end do
      end do
      do i=1,3
        do j=1,3
          tmatr(i,j)=0.0
        end do
      end do
c
c***  Triangle inferior de la matriu de rigidesa local
c
      x1=coord(1,node1)
      y1=coord(2,node1)
      x2=coord(1,node2)
      y2=coord(2,node2)
      blong=dsqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
      eal =props(1)*props(2)/blong
      eil1=props(1)*props(3)/blong
      eil2=eil1/blong
      eil3=eil2/blong
      estil(1,1)=       eal
      estil(4,1)= -     eal
      estil(2,2)=  12.0*eil3
      estil(3,2)=   6.0*eil2
      estil(5,2)= -12.0*eil3
      estil(6,2)=   6.0*eil2
      estil(3,3)=   4.0*eil1
      estil(5,3)= - 6.0*eil2
      estil(6,3)=   2.0*eil1
      estil(4,4)=       eal
      estil(5,5)=  12.0*eil3
      estil(6,5)= - 6.0*eil2
      estil(6,6)=   4.0*eil1
c
c***  Triangle superior de la matriu de rigidesa local
c
      do i=2,6
        do j=1,i-1
          estil(j,i)=estil(i,j)
        end do
      end do
c
c***  Matriu de canvi de coordenades locals a globals
c
      calfa=(x2-x1)/blong
      salfa=(y2-y1)/blong
      tmatr(1,1)= calfa
      tmatr(2,1)= salfa
      tmatr(1,2)=-salfa
      tmatr(2,2)= calfa
      tmatr(3,3)=   1.0
c
c***  Calcul de la matriu Kl.Tt 
c
      do i=1,3
        do j=1,3
          do k=1,3
            etmat(i  ,j  )=etmat(i  ,j  )+estil(i  ,k  )*tmatr(j,k)
            etmat(i+3,j  )=etmat(i+3,j  )+estil(i+3,k  )*tmatr(j,k)
            etmat(i  ,j+3)=etmat(i  ,j+3)+estil(i  ,k+3)*tmatr(j,k)
            etmat(i+3,j+3)=etmat(i+3,j+3)+estil(i+3,k+3)*tmatr(j,k)
          end do
        end do
      end do
c
c***  Calcul de la matriu de rigidesa global, K = T.Kl.Tt
c
      do i=1,3
        do j=1,3
          do k=1,3
            estif(i  ,j  )=estif(i  ,j  )+tmatr(i,k)*etmat(k  ,j  )
            estif(i+3,j  )=estif(i+3,j  )+tmatr(i,k)*etmat(k+3,j  )
            estif(i  ,j+3)=estif(i  ,j+3)+tmatr(i,k)*etmat(k  ,j+3)
            estif(i+3,j+3)=estif(i+3,j+3)+tmatr(i,k)*etmat(k+3,j+3)
          end do
        end do
      end do

      return
      end
      subroutine assemb(ielem,lnods,estif,rigid)
c*****************************************************************************
c
c   Fa l'assembly de les matrius de rigidesa de cada barra en la matriu global
c
c*****************************************************************************
      implicit real*8(a-h,o-z)
      common /contrl/ npoin,nelem,nmats,nvfix,nload,nband,ntotv
      real*8    estif(6,6), rigid(nband+1,ntotv)
      integer*4 lnods(2)

      do inode=1,2
        ipoin=lnods(inode)
        do idofn=1,3
          ievab=(inode-1)*3+idofn
          itotv=(ipoin-1)*3+idofn
          do jnode=1,2
            jpoin=lnods(jnode)
            do jdofn=1,3
              jevab=(jnode-1)*3+jdofn
              jtotv=(jpoin-1)*3+jdofn
              jband=jtotv-itotv+1
              if (jband.ge.1)
     .        rigid(jband,itotv)=rigid(jband,itotv)+estif(ievab,jevab)
            end do
          end do
        end do
      end do

      return
      end
      subroutine vincul(rigid,gload,treac,nofix,ifpre,presc)
c*****************************************************************************
c
c   Imposa les condicions de contorn
c
c*****************************************************************************
      implicit real*8(a-h,o-z)
      common /contrl/ npoin,nelem,nmats,nvfix,nload,nband,ntotv
      real*8     rigid(nband+1,ntotv), gload(ntotv), 
     .           treac(3,nvfix),       presc(3,nvfix)
      integer*4  nofix(nvfix), ifpre(3,nvfix)
c
c***  Inicialitzacio de les reaccions per als graus de llibertat prescrits
c
      do ivfix=1,nvfix
        do idofn=1,3
          if (ifpre(idofn,ivfix).eq.1) then
            itotv=(nofix(ivfix)-1)*3+idofn
            treac(idofn,ivfix)=-gload(itotv)
          end if
        end do
      end do
c
c***  Llac sobre els nodes coaccionats
c
      do ivfix=1,nvfix
        jpoin=nofix(ivfix)
        do jdofn=1,3
          if (ifpre(jdofn,ivfix).eq.1) then
c
c***  Cas de grau de llibertat prescrit
c
            jtotv=(jpoin-1)*3+jdofn
c
c***  Modificacio de les equacions anteriors a la del g.d.ll. prescrit
c
            if (jtotv.gt.1) then 
              itot1=jtotv-nband
              if (itot1.lt.1) itot1=1
              do itotv=itot1,jtotv-1
                jband=jtotv-itotv+1
                gload(itotv)=gload(itotv)
     .                      -rigid(jband,itotv)*presc(jdofn,ivfix)
                rigid(jband,itotv)=0.0
              end do
            end if
c
c***  Modificacio de les equacions posteriors a la del g.d.ll. prescrit
c
            if (jtotv.lt.ntotv) then 
              itot2=jtotv+nband
              if (itot2.gt.ntotv) itot2=ntotv
              do itotv=jtotv+1,itot2
                jband=itotv-jtotv+1
                gload(itotv)=gload(itotv)
     .                      -rigid(jband,jtotv)*presc(jdofn,ivfix)
                rigid(jband,jtotv)=0.0
              end do
            end if
c
c***  Equacio trivial per al grau de llibertat prescrit
c
            rigid(1,jtotv)=1.0
            gload(jtotv)=presc(jdofn,ivfix)
          end if
        end do
      end do 

      return
      end
      subroutine solsis(nequa,nsban,amatr,cmatr,nsist,bvect,xvect)
c*****************************************************************************
c
c     Resolucio de sistemes lineals amb matriu en banda, simetriques i  
c     definides positives emprant la descomposicio de Cholesky.
c
c             A . x = b                  L . z = b
c                             -----> 
c             L . Lt= A                  Lt. x = z
c
c     nequa = num. d'equacions = num.d'incognites
c     nsban = semiample de banda
c     amatr = triangle superior de la matriu de coeficients emmagatzemat
c             en banda per files. Si A es la matriu original,
c             amatr(i,j) = A(j,i+j-1)
c     cmatr = matriu L transposta de la descomposicio de Cholesky, 
c             emmagatzemada en banda per files
c     nsist = num. de sistemes que cal resoldre
c     bvect = vector de termes independents
c     xvect = vector d'incognites
c
c*****************************************************************************
      implicit real*8(a-h,o-z)
      real*8  amatr(nsban+1,nequa),bvect(nequa,nsist),
     .        cmatr(nsban+1,nequa),xvect(nequa,nsist)

      call cholesky (nequa,nsban,amatr,cmatr)
      call sub1 (nequa,nsban,cmatr,nsist,bvect,xvect)
      call sub2 (nequa,nsban,cmatr,nsist,xvect,xvect)
c
      return
      end
c
c***  Descomposicio de Cholesky 
c
      subroutine cholesky(n,l,a,c)
      implicit real*8(a-h,o-z)
      real*8  a(l+1,n),c(l+1,n)

      do i = 1,n
        c(1,i) = a(1,i)
        k1 = i-l
        if (k1.lt.1) k1 = 1
        if (k1.le.i-1) then 
          do k = k1,i-1
            c(1,i) = c(1,i) - c(i-k+1,k)*c(i-k+1,k)
          end do
        end if
        if (c(1,i).le.0.) call error(' Element no positiu a la fila ',i)
        c(1,i) = dsqrt(c(1,i))
        j2 = i + l
        if (j2.gt.n) j2 = n
        if (i+1.le.j2) then 
          do j = i+1,j2
            c(j-i+1,i) = a(j-i+1,i)
            k1 = j - l
            if (k1.lt.1) k1 = 1
            if (k1.le.i-1) then 
              do k = k1,i-1
                c(j-i+1,i) = c(j-i+1,i) - c(i-k+1,k)*c(j-k+1,k)
              end do
            end if
            c(j-i+1,i) = c(j-i+1,i)/c(1,i)
          end do
        end if
      end do

      return
      end
c
c***  Substitucio endavant
c
      subroutine sub1(n,l,c,m,b,z)
      implicit real*8(a-h,o-z)
      real*8 c(l+1,n),b(n,m),z(n,m)

      do is = 1,m
        z(1,is) = b(1,is)/c(1,1)
        do i = 2,n
          z(i,is) = b(i,is)
          k1 = i-l
          if (k1.lt.1) k1 = 1
          if (k1.le.i-1) then 
            do k = k1,i-1
              z(i,is) = z(i,is)-c(i-k+1,k)*z(k,is)
            end do
          end if
          if (c(1,i).le.0.) call error(' Pivot no positiu a la fila ',i)
          z(i,is) = z(i,is)/c(1,i)
        end do
      end do
 
      return
      end
c
c***  Substitucio enrere
c
      subroutine sub2(n,l,c,m,z,x)
      implicit real*8(a-h,o-z)
      real*8  c(l+1,n),z(n,m),x(n,m)

      do is = 1,m
        x(n,is) = z(n,is)/c(1,n)
        do i = n-1,1,-1
          x(i,is) = z(i,is)
          k2 = i+l
          if (k2.gt.n) k2 = n
          if (i+1.le.k2) then 
            do k = i+1,k2
              x(i,is) = x(i,is)-c(k-i+1,i)*x(k,is)
            end do
          end if
          if (c(1,i).le.0.) call error(' Pivot no positiu a la fila ',i)
          x(i,is) = x(i,is)/c(1,i)
        end do
      end do
 
      return
      end
c
c***  Sortida de possibles errors
c
      subroutine error(alfa,i)
      implicit real*8(a-h,o-z)
      character alfa*(*)
      write(6,'(a,i5)') alfa,i
      stop
      end
      subroutine outdes(tdisp)
c*****************************************************************************
c
c   Escriptura dels moviments de l'estructura a l'arxiu de resultats
c
c*****************************************************************************
      implicit real*8(a-h,o-z)
      character c_coor*12,c_prop*14,c_dofn*9,c_load*13,c_stre*10
      common /contrl/ npoin,nelem,nmats,nvfix,nload,nband,ntotv
      common /arxius/ n_lect,n_escr,n_kltt,n_tran
      common /titols/ c_coor(2),c_prop(3),c_dofn(3),c_load(3),c_stre(3)
      real*8 tdisp(ntotv)

      write(n_escr,100) (c_dofn(i),i=1,3)
  100 format(//' MOVIMENTS DELS NODES DE L''ESTRUCTURA',/,
     .         ' --------- ---- ----- -- ------------',//,
     .         '        Num. de node',3(5x,a9))
      do ipoin=1,npoin
        itotv=3*(ipoin-1)
        write(n_escr,10) ipoin,(tdisp(itotv+i),i=1,3)
      end do
  10  format(i20,3(f14.8))

      return
      end
      subroutine calesf(tdisp,treac,etmat,tmatr,lnods,matno,
     .                  nofix,ifpre)
c*****************************************************************************
c
c   Calcula i escriu els esforcos d'extrem de  barra i les reaccions
c
c*****************************************************************************
      implicit real*8(a-h,o-z)
      character c_coor*12,c_prop*14,c_dofn*9,c_load*13,c_stre*10
      common /contrl/ npoin,nelem,nmats,nvfix,nload,nband,ntotv
      common /arxius/ n_lect,n_escr,n_kltt,n_tran
      common /titols/ c_coor(2),c_prop(3),c_dofn(3),c_load(3),c_stre(3)
      real*8    tdisp(ntotv), treac(3,nvfix),
     .          etmat(6,6), tmatr(3,3), sextr(6), sglob(6) 
      integer*4 lnods(2,nelem),matno(nelem),nofix(nvfix),
     .          ifpre(3,nvfix)
c
c***  Rebobina els arxius de Kl.Tt i T
c
      rewind(n_kltt)
      rewind(n_tran)
c
c***  Titols per escriure els esforcos d'extrem de barra
c
      write(n_escr,100) (c_stre(i),i=1,3)
  100 format(//,' ESFORCOS EN EXTREMS DE BARRA',/,
     .          ' -------- -- ------- -- -----',//,
     .          '  Barra','            ',3(5x,a10))
c
c***  Llac sobre les barres
c
      do ielem=1,nelem
c
c***  Calcul dels esforcos d'extrem de barra (P_l = K_l . T^t . d_g): SEXTR
c
        read(n_kltt) ((etmat(i,j),i=1,6),j=1,6)
        do inode=1,2
          do idofn=1,3
            stres=0.0
            ievab=(inode-1)*3+idofn
            do jnode=1,2
              jpoin=lnods(jnode,ielem)
              do jdofn=1,3
                jevab=(jnode-1)*3+jdofn
                jtotv=(jpoin-1)*3+jdofn
                stres=stres+etmat(ievab,jevab)*tdisp(jtotv)
              end do
            end do
            sextr(ievab)=stres
          end do
        end do
c
c***  Escriptura dels esforcos d'extrem de barra en coordenades locals
c
        ipoi1=lnods(1,ielem)
        ipoi2=lnods(2,ielem)
        write(n_escr,10) ielem,'Node 1',ipoi1,(sextr(i),i=1,3)
        write(n_escr,11)       'Node 2',ipoi2,(sextr(i),i=4,6)
c
c***  Calcul de reaccions en el cas de barres coaccionades
c
        if (matno(ielem).lt.0) then
          read(n_tran) ((tmatr(i,j),i=1,3),j=1,3)
c
c***  Esforcos d'extrem de barra en coordenades globals: SGLOB
c
          do i=1,3
            sglob(i  )=0.0
            sglob(i+3)=0.0
            do j=1,3
              sglob(i  )=sglob(i  )+tmatr(i,j)*sextr(j  )
              sglob(i+3)=sglob(i+3)+tmatr(i,j)*sextr(j+3)
            end do
          end do 
c
c***  Identificacio del g.d.ll. prescrit i obtencio de la reaccio
c
          do ivfix=1,nvfix
            if(nofix(ivfix).eq.ipoi1) then
              do i=1,3
                if(ifpre(i,ivfix).eq.1)
     .          treac(i,ivfix)=treac(i,ivfix)+sglob(i  )
              end do
            else if(nofix(ivfix).eq.ipoi2) then
              do i=1,3
                if(ifpre(i,ivfix).eq.1)
     .          treac(i,ivfix)=treac(i,ivfix)+sglob(i+3)
              end do
            end if
          end do
        end if
      end do
c
c***  Escriptura de les reaccions
c
      write(n_escr,101) (c_load(i),i=1,3)
  101 format(//' REACCIONS EN NODES COACCIONATS',/,
     .         ' --------- -- ----- -----------',//,
     .         '    Node prescrit',3(7x,a13))
      do ivfix=1,nvfix
        ipoin=nofix(ivfix)
        write(n_escr,12) ipoin,(treac(i,ivfix),i=1,3)
      end do
c
c***  Formats
c
  10  format(i7,2x,a6,i4,3f15.8)
  11  format(7x,2x,a6,i4,3f15.8)
  12  format(i17,3f20.8)

      return
      end
      subroutine tanarx
c*****************************************************************************
c
c   Tancament dels arxius emprats en el programa
c
c*****************************************************************************
      implicit real*8(a-h,o-z)
      common /arxius/ n_lect,n_escr,n_kltt,n_tran

      close(unit=n_lect,status='keep')
      close(unit=n_escr,status='keep')
      close(unit=n_kltt,status='delete')
      close(unit=n_tran,status='delete')

      return
      end
