      subroutine elmt05(d,ul,xl,ix,tl,s,p,ndf,ndm,nst,isw)

c     Purpose: Geometrical nonlinear axisymmetric shell: (c) w.wagner

c     Inputs:
c        d(*)      - Material set parameters
c        ul(ndf,*) - Solution variables for element
c        xl(ndf,*) - Nodal coordinates for element
c        ix(*)     - Nodal connections for element
c        tl(*)     - Nodal temperatures for element
c        ndf       - Number dof/node
c        ndm       - Spatial dimension of FEM mesh
c        nst       - Size of element matrix
c        isw       - Switch parameter to control action taken by routine

c     Outputs:
c        s(nst,*)  - Element matrix
c        p(*)      - Element vector

      implicit real*8(a-h,o-z)

      character*4 head,yyy*80
      logical     bs
      integer     ix(*)
      real*8      xl(ndm,*),tl(*),dm
      real*8      d(7),s(nst,*),shp(2),dd(5,5),btd(5,3),
     1            bm(5,3),sig(5),ul(ndf,*),vl(3,2),p(nst),dot
      common /bdata/  head(20)
      common /cdata/  numnp,numel,nummat,nen,neq
      common /eldata/ dm,n,ma,mct,iel,nel
      common /iofile/ ior,iow
      common /ydata/  yyy
c.... input material properties
      if(isw.eq.1) then
        call matl05(d)
        call pconsd(dd,25,0.0d0)
        call pconsd(bm,15,0.0d0)
      else
c.... length, angle, radius, jacobian
        sn = xl(2,2)-xl(2,1)
        cs = xl(1,2)-xl(1,1)
        sl = sqrt(cs*cs + sn*sn)
        bs = d(13).eq.0.0d0
        if(isw.eq.2) then
c.... check mesh for error
          if(sl.le.0.0d0) then
            write(iow,3001) n
            if(ior.lt.0) write(*,3001) n
          endif
        elseif(isw.eq.5) then
c.... form a lumped mass matrix; length, volume*density
          dv    = d(11)*sl/8.0
          if(bs) then
            p(1)     = 4.0*dv
            p(ndf+1) = p(1)
          else
            p(1)     = dv*(3.*xl(1,1) +    xl(1,2))
            p(ndf+1) = dv*(   xl(1,1) + 3.*xl(1,2))
          endif
          p(2)     = p(1)
          p(3)     = p(1)*d(12)
          p(ndf+2) = p(ndf+1)
          p(ndf+3) = p(ndf+1)*d(12)
        else
c.... form stiffness/residual
          sn   = sn/sl
          cs   = cs/sl
          if(bs) then
            r    = 1.0
            recr = 0.0
          else
            r    = 0.5*(xl(1,1) + xl(1,2))
            recr = 1.0/r
          endif
          dv   = sl*r
c.... shape function derivatives
          shp(2) =  1.0/sl
          shp(1) = -shp(2)
c.... local  displacements
          do 300 k = 1,2
            vl(1,k) = cs*ul(1,k) + sn*ul(2,k)
            vl(2,k) =-sn*ul(1,k) + cs*ul(2,k)
            vl(3,k) = ul(3,k)
300       continue
c.... derivative w,s
          lin = d(5)
          if(lin.eq.0) then
            wks = 0.0
          else
            wks = (vl(2,2) - vl(2,1))/sl
          endif
c.... stresses, strains, d-matrix
          call modl05(sig,vl,dd,d,sn,cs,sl,recr,wks)
          if(mod(isw,3).eq.0) then
            i    = ndf + 2
c.... Load vector in local coordinates (reference system)
            f = d(4)*sl/8.0*dm
            if(bs) then
              p(2) = f*4.0
              p(i) = p(2)
            else
              p(2) = f*(3.*xl(1,1) +    xl(1,2))
              p(i) = f*(   xl(1,1) + 3.*xl(1,2))
            endif
c.... K-sigma tangent matrix
            if(lin.ne.0) then
              s(2,2) =   r*sig(1)/sl
              s(i,i) =   s(2,2)
              s(2,i) = - s(2,2)
            endif
c.... multiply stress and moduli by jacobian
            do 302 k = 1,5
              sig(k) = sig(k)*dv
              do 301 j = 1,5
                dd(j,k) = dd(j,k)*dv
301           continue
302         continue
            i1=0
            do 309 ii = 1,2
c.... residual G = P - Bt*S 
              call bmat05(bm,shp(ii),sn,cs,recr,wks)
              do 303 i = 1,3
                p(i1+i) = p(i1+i) - dot(bm(1,i),sig,5)
303           continue
c.... tangent stiffness matrix
              if(isw.eq.3) then
                do  305 i = 1,3
                  do 304 k = 1,5
                    btd(k,i) = dot(bm(1,i),dd(1,k),5)
304               continue
305             continue
                j1 = i1
                do 308 jj = ii,2
                  call bmat05(bm,shp(jj),sn,cs,recr,wks)
                  do 307  i = 1,3
                    do 306  j = 1,3
                      s(i1+i,j1+j)=s(i1+i,j1+j)+dot(btd(1,i),bm(1,j),5)
306                 continue
307               continue
                  j1 = j1 + ndf
308             continue
              endif
              i1 = i1 + ndf
309         continue
c.... lower part of stiffness matrix and transformation to global frame
            if(isw.eq.3) then
              do 310 i = 1,3
              do 310 j = 1,3
                s(i+ndf,j) = s(j,i+ndf)
310           continue
              call tran05(s,cs,sn,nst,ndf,1)
            endif
            call tran05(p,cs,sn,nst,ndf,2)
c.... output stresses (N, M, Q)
          elseif(isw.eq.4) then
            mct = mct - 1
            if(mct.le.0) then
               write(iow,2001) head
               if(ior.lt.0) write(*,2001) head
               mct = 50
            endif
            r = 0.5*(xl(1,1) + xl(1,2))
            z = 0.5*(xl(2,1) + xl(2,2))
            write(iow,2002) n,ma,r,z,sig
            if(ior.lt.0) write(*,2002) n,ma,r,z,sig
          endif
        endif
      endif
c.... format statements
2001  format(1x,20a4//2x,'E L E M E N T   S T R E S S E S'//
     1  '  El  Mat  1-Coor  2-Coor   **NS**    **NPHI**'
     2  '    **MS**    **MPHI**    **QS**'/1x)
2002  format(i5,i4,0p2f8.3,5(1x,1p1e10.3))
3001  format(' *ERROR* Element',i5,' has zero length.')
      end
c
      subroutine bmat05(bm,shp,sn,cs,recr,wks)
      implicit real*8(a-h,o-z)
      real*8 bm(5,3)
c
c.... nonlinear B-matrix for axisymmetric shells
      bm(1,1) = shp
      bm(1,2) = shp*wks
      bm(2,1) = 0.5*cs*recr
      bm(2,2) =-0.5*sn*recr
      bm(3,3) = shp
      bm(4,1) = bm(2,1)*sn*recr
      bm(4,2) = bm(2,2)*sn*recr
      bm(4,3) = bm(2,1)
      bm(5,2) = shp
      bm(5,3) =-0.5
c
      end
c
      subroutine matl05(d)
      character yyy*80,typ*5,type*5
      real*8    d(6)
      common /iofile/ ior,iow
      common /ydata/  yyy
      data   type/'beam '/
c
 1    if(ior.lt.0) then
        write(*,3000)
        read(*,1000) typ
      else
        read(ior,1000) typ
      endif
      call pintio(yyy,10)
      read(yyy,1001,err=2) d
      if(typ.eq.type) then
        d(2)  = 0.0
        d(13) = 0.0
      else
        d(13) = 1.0
      endif
      write(iow,2000) typ,d
      if(ior.lt.0) write(*,2000) typ,d
c.... set beam/shell in-plane and bending stiffness values
      d(7)  = d(1)*d(3)/(1.0 - d(2)*d(2))
      d(8)  = d(2)*d(7)
      d(9)  = d(7)*d(3)*d(3)/12.
      d(10) = d(2)*d(9)
      d(11) = d(6)*d(3)
      d(12) = d(3)*d(3)/12.
      return
 2    call perror('PCELM5',yyy)
      go to 1
c
1000  format(a5)
1001  format(7f10.0)
2000  format(5x,'Rectangular Beam/Axisymmetric Shell Model'/
     x 5x,'Type: ',a5/5x,'Elastic Modulus',g15.4/
     1 5x,'Poisson Ratio',g17.4/5x,'Thickness  ',g19.4/
     2 5x,'Normal Load  ',g17.4/
     3 5x,'Linear(0=l 1=nl)',f9.1/5x,'Mass Density',g18.4/1x)
3000  format(' Input: 1. Type (beam or shell)'/
     1       '        2. E nu h press lin rho'/3x,'>',$)
      end
c
      subroutine modl05(sig,vl,dd,d,sn,cs,sl,recr,wks)
      implicit real*8(a-h,o-z)
      real*8 sig(5),eps(5),vl(3,2),dd(5,5),d(1)
c
c.... elasticity matrix for axisymmetric shells
      dd(1,1) = d(7)
      dd(2,2) = d(7)
      dd(1,2) = d(8)
      dd(2,1) = d(8)
      dd(3,3) = d(9)
      dd(4,4) = d(9)
      dd(3,4) = d(10)
      dd(4,3) = d(10)
      dd(5,5) = 5.*d(7)*(1.0 - d(2))/12.0
c...  strains
      v1     = 0.5*((vl(1,1) + vl(1,2))*cs - (vl(2,1) + vl(2,2))*sn)
      eps(1) = (vl(1,2) - vl(1,1))/sl + 0.5*wks*wks
      eps(2) =  v1*recr
      eps(3) = (vl(3,2) - vl(3,1))/sl
      eps(4) = (v1*sn*recr +0.5*(vl(3,1) + vl(3,2))*cs)*recr
      eps(5) = (vl(2,2) - vl(2,1))/sl - 0.5*(vl(3,1) + vl(3,2))
c.... stresses
      sig(1) = dd(1,1)*eps(1) + dd(1,2)*eps(2)
      sig(2) = dd(2,1)*eps(1) + dd(2,2)*eps(2)
      sig(3) = dd(3,3)*eps(3) + dd(3,4)*eps(4)
      sig(4) = dd(4,3)*eps(3) + dd(4,4)*eps(4)
      sig(5) = dd(5,5)*eps(5)
c
      end
c
      subroutine tran05(s,cs,sn,nst,ndf,itype)
      implicit real*8 (a-h,o-z)
      real*8 s(nst,1)
c
c.... itype: 1  Transform matrix s(nst,nst)
c            2  Transform vector s(nst,1)
c
      if(itype.eq.1) then
        do 12 i = 1,nst,ndf
          do 11 j = 1,nst
            t        = s(j,i)*cs - s(j,i+1)*sn
            s(j,i+1) = s(j,i)*sn + s(j,i+1)*cs
            s(j,i  ) = t
11        continue
12      continue
        nn = nst
      else
        nn = 1
      endif
c
      do 14 i = 1,nst,ndf
        do 13 j = 1,nn
          t        = s(i,j)*cs - s(i+1,j)*sn
          s(i+1,j) = s(i,j)*sn + s(i+1,j)*cs
          s(i  ,j) = t
13      continue
14    continue
c
      end
