subroutine ushell(thick,xintp,ncrd,m,nn)
c
c CONTACT: ap@marc.de
c      
c
      implicit real*8 (a-h,o-z)                                               dp
      dimension xintp(*),m(*)
c
c* * * * *
c
c Subroutine to specify shell element thickness
c for each integration point.
c
c* * * * *
c
c INPUT
c   xintp     integration point coordinates
c   ncrd      number of coordinates
c   mel(1)    user element number
c   mel(2)    internal element number
c   nn        integration point number
c
c OUTPUT
c   thick     total element thickness
c
c* * * * *
c
c
      return
      end
      subroutine ulthck(ht,mel,nip,lay,matno,ktp,hk,xip,ncrd)
c
      implicit real*8 (a-h,o-z)
      dimension xip(*),mel(*)
c
c* * * * *
c
c Subroutine to redefine the layer thickness as
c a fraction of the total element thickness in
c a composite shell element.
c
c* * * * *
c
c INPUT
c   mel(1)    user element number
c   mel(2)    internal element number
c   nip       integration point number
c   lay       layer number
c   matno     composite group number
c   ktp       total number of layers in this group
c   hk        total element thickness
c   xip       integration point coordinates
c   ncrd      number of coordinates
c   ht        fraction of total thickness for this layer as input
c             on composite input
c
c OUTPUT
c   ht        redefined fraction of total thickness for this layer
c
c                   current layer thickness
c             ht =  -----------------------
c                   total element thickness
c
c
c NOTE
c   The user should only define variable ht and leave all
c   other variables unaltered.
c   The user must also make sure that the sum of all
c   layer thickness fractions will equal to one (i.e.
c   the sum of all layer thicknesses equals to the
c   total element thickness).
c
c* * * * *
c
c
      return
      end
      subroutine ssect2(kcx,ktp,geom,hk,fkc,fkc2,htq,ar,fkcold)
c
      implicit real*8 (a-h,o-z)
c
      include '../common/matdat'
      include '../common/maters'
      include '../common/arrays'
      include '../common/elmcom'
      include '../common/space'
      include '../common/pival'
      include '../common/nzro1'
      include '../common/heat'
      include '../common/lass'
      include '../common/dimen'
      include '../common/array4'
      include '../common/far'
      dimension geom(*),mn(2)
c
c input:
c   kcx   =layer number
c   ktp   =number of layers
c   fkcold=position of previous layer
c   geom  =geometry parameters
c   hk    =total thickness or area
c output:
c   fkc   =absolute position of current layer
c   fkc2  =position of current layer for bending element
c   htq   =thickness of current layer
c   ar    =relative position of current layer
c
      if(ishell.eq.0) then
         htq=1.d0
      else if(lclass.eq.12 .or. lclass.eq.13) then
         htq=1.d0
         hk=geom(1)
         hk1=geom(2)
         call csect(hk,hk1,kcx,m,jtype,xc,yc,tc,phi,wg,xcb,ycb,xjot,
     *       xsh,ysh,rad,r1,r2,inpsxc,npp,nsecs)
         if(lclass.eq.13) phi=0.0d0
      else
         if(icomps.eq.1) then
            if(kcx.eq.1) then
               fkc0=vars(icompd-1 + ncompd*(matno-1) + 3)
               fkcold=hk/2.-fkc0
            endif
            ht = vars(ilycmp-1 + 2*nseqst*(matno-1) + 2*(kcx-1) + 1)
            la1 = icrxpt + (n-1)*nelstr + (nn-1)*ncrd
            mn(1) = m
            mn(2) = n
            call ulthck(ht,mn,nn,kcx,matno,ktp,hk,vars(la1),ncrd)
            htq=hk*ht
            q=1.d0
            fkc=fkcold-htq/2.
            fkc2=fkc
            ar=2.0d0*fkc/hk
            fkcold=fkcold-htq
         else if(lbend.eq.1) then
            kt2=ktp/2
            fko=real(kt2-kcx+1)*hk/real(ktp-1)
            fkc=sin(fko*pi)
            fkc2=cos(fko*pi)
            q=1.0d0
            if(kcx.eq.1.or.kcx.eq.ktp) q=0.5d0
            htq=hk*q/real(ktp-1)
       else if(ktp.eq.1) then
          fac=1.d0
          htq=hk
          ar=0.d0
          fkc=0.d0
          fkc2=0.d0
       else
            fac=2.d0/3.d0
            if(mod(kcx,2).eq.0) fac=4.d0/3.d0
            if(kcx.eq.1.or.kcx.eq.ktp) fac=1.d0/3.d0
            htq=hk*fac/real(ktp-1)
            ar=1.d0-2.d0*real(kcx-1)/real(ktp-1)
            fkc=ar*hk/2.
            fkc2=fkc
         endif
      endif
c
      return
      end