      subroutine ldecom (is,i2, p1,p2,p3,p4, d,d1,d2,d3,d4) !, r,z)
c
c     approximate decomposition of the finite element matrix
c     outer codiagonals of decomposed matrix at the same position as for
c     the original f.E.M. Matrix, in diagonal units
c     diagonal array d stores reciprocal values of the diagonal of the
c     decomposed matrix    (detail see sem/86, p.906)
c     ***** iselin's correction made, diagonal array p=1.d0
c

!z      implicit real*8 (a-h,o-z)
      implicit none
      integer*4  is,i2,                        i,ik
      real*8     P1(is),p2(is),p3(is),p4(is),                           &
     &     d(is),d1(is),d2(is),d3(is),d4(is),  sum   ! ,  r,z
c
c     is  - number of unknowns
c     i2  - number of unknowns in r direction
c     F.E.M matrix is symmetrical, on the diagonal are P(i)=1.
c     system is described as:
c     P(i)*V(i)+P1(i)*V(i+1)  +P2(i)*V(i+i2-1)+    P3(i)*V(i+i2)   +P4(i+i2+1)*V(i+i2+1)+
c               P1(i-1)*V(i-1)+P2(i-i2+1)V(i-i2+1)+P3(i-i2)*V(i-i2)+P4(i-i2-1)*V(i-i2-1)=B(i)
c
c  B(i) = P4(i-i2-1) * V(i-i2-1)
c       + P3(i-i2)   * V(i-i2)
c       + P2(i-i2+1) * V(i-i2+1)
c       + P1(i-1)    * V(i-1)
c       + P (i)      * V(i)
c       + P1(i)      * V(i+1)
c       + P2(i)      * V(i+i2-1)
c       + P3(i)      * V(i+i2)
c       + P4(i+i2+1) * V(i+i2+1)
c
c   input
c     1.,P1,... - The coefficients of the f.E.M. Matrix (diagonals)
c   output
c     d1,..,D4 - the coefficients of the approx.Decomposed upper triang. Matrix
c     d - reciprocal of diagonal terms from d
c
c     first i2 lines
c          (i2 - is the number of unknowns in R direction)
c
!        if ( i == 1 ) then
            d(1)  = 1.d0
            d1(1) = p1(1)
            d4(1) = p4(1)
            d3(1) = p3(1)
            d2(1) = p2(1)
!
      do i = 2,is

        if ( i < i2 ) then             ! 2:i2-1
            sum   = d(i-1) * d1(i-1)**2
            d(i)  = (1.d0  - sum )
            d1(i) =  p1(i)                              / d(i)

        else if ( i == i2 ) then       ! i2 , ik=0
            sum   = d(i-1)  * d1(i-1)**2                                &
     &            + d(1)    * d2(1)**2                                  &
            d(i)  = (1.d0  - sum )
            d1(i) = (p1(i) - d(1)  * d3(1) * d2(1)   ) / d(i)

        else if ( i == i2+1 ) then     ! i2+1 , ik=1
            sum   = d(i-1) * d1(i-1)**2                                 &
     &            + d(2)   * d2(2)**2                                   &
     &            + d(1)   * d3(1)**2
            d(i)  = (1.d0  - sum )
            d1(i) = (p1(i) - d(1)  * d3(1) * d4(1)                      &
     &                     - d(2)  * d3(2) * d2(2)   ) / d(i)
     
        else                           ! i2+2:is
            ik    = i-i2
            sum   = d(i-1)  * d1(i-1)**2                                &
     &            + d(ik+1) * d2(ik+1)**2                               &
     &            + d(ik)   * d3(ik)**2                                 &
     &            + d(ik-1) * d4(ik-1)**2
            d(i)  = (1.d0  - sum )
            d1(i) = (p1(i) - d(ik)   * d3(ik)   * d4(ik)                &
     &                     - d(ik+1) * d3(ik+1) * d2(ik+1)  )/d(i)
        end if
            d4(i) = (p4(i)                                ) / d(i)
            d3(i) = (p3(i) - d(i-1)  * d1(i-1)  * d4(i-1) ) / d(i)
            d2(i) = (p2(i) - d(i-1)  * d1(i-1)  * d3(i-1) ) / d(i)
      end do

      end subroutine ldecom
      