SUBROUTINE fddisc(nLay,nSubLayR,nCells,
     &                  subMode,l,lambda,rho,cp,C,R)

C     Last change:  LE   1 Nov 2001


C Purpose:    The subroutine calculates the thermal resistances and
C             the heat capacities of a multilayer building component
C             on basis that the layers are manually subdivided.

C Called by:  fdwall.nmf
C Calls:      None

C Limits:     No automatical discretization.

C Date:       991001
C Made by:    Jan Akander
C             Axel Bring
C             
C Input:
C      nLay   = number of material layers in the building component
C   nSubLay   = vector with number of subdivisions (cells) of material layer
C    nCells   = number of nodes (cells) in the RC-chain
C   subMode   = method of dividing layers to cells
C                   0 - uniform;  1 - tighter at the beginning
C                   2 - tighter at the end;  3 - tighter at both sides
C         l   = vector with material layer thickness
C    lambda   = vector with material layer heat conductivity
C       rho   = vector with material layer density
C        cp   = vector with material layer specific heat capacity

C  
C
C Output:
C      C = vector with discretized heat capacities (elements =sum(nSubLay))
C      R = vector with discretized thermal resistances (elements =sum(nSubLay)+1) 
 
C integers in
      INTEGER nLay,nCells, subMode
C pseudo integers in
      REAL*8 nSubLayR(nLay)
C local integers
      INTEGER i,j,k, ni
C parameters in
      DOUBLE PRECISION l(nLay),rho(nLay),cp(nLay),lambda(nLay),
     &  C(nCells),R(nCells+1)
C The point with most sparse cells
      DOUBLE PRECISION  xSparse
C Temp vars
      DOUBLE PRECISION a, b, qa, qb, qar, qbr, niR, cf, w,
     &       Rprel, Cprel, rPrev, xPrev, xj, wj, rj

c      WRITE(7,99) 'nlay,ncells,nofrs',nlay,ncells
c      WRITE(7,98) (nsublayR(i),i=1,nlay)
c      WRITE(7,98) (l(i),i=1,nlay)
c      WRITE(7,98) (lambda(i),i=1,nlay)
c      WRITE(7,98) (rho(i),i=1,nlay)
c      WRITE(7,98) (cp(i),i=1,nlay)
c   99 FORMAT(1x,a,3i5)
c   98 FORMAT(1x,5f8.2)



      k=0
      rPrev = 0.0

      IF (subMode .EQ. 0) THEN
        DO 200 i=1,nLay
          niR = nSubLayR(i)
          ni = NINT(niR)
          Rprel = (l(i)/niR)/(2*lambda(i))
          Cprel = (l(i)/niR)*rho(i)*cp(i)
          DO 300 j=1,ni
            k=k+1

            C(k)=Cprel

            IF (j.EQ.1) THEN
              R(k) = rPrev + Rprel
            ELSE
              R(k) = 2*Rprel
            END IF

300       CONTINUE
          rPrev = Rprel
200     CONTINUE

      ELSE
C find the thickness and the point with most sparse cells
        w = 0.0
        DO 100 i=1,nLay
          w = w + l(i)
100     CONTINUE

        IF (subMode .EQ. 1) THEN
          xSparse = w
        ELSE IF (subMode .EQ. 2) THEN
          xSparse  = 0.0
        ELSE
          xSparse = w/2
        END IF


        b = 0.0
        qb = 0.0
        qbr = SQRT(w)
        rPrev = 0.0

        DO 210 i=1,nLay
          niR = nSubLayR(i)
          ni = NINT(niR)
          Rprel = 1.0/(2*lambda(i))
          Cprel = rho(i)*cp(i)

          a = b
          qa = qb
          qar = qbr
          b = a + l(i)
          qb = SQRT(b)
cle>      qbr = SQRT(w-b)
          qbr = SQRT(MAX(0.0d0,w-b))

          IF (b .LE. xSparse) THEN
            cf = (qb - qa) / niR
          ELSE IF (a .GE. xSparse) THEN
            cf = (qar - qbr) / niR
          ELSE
            cf = (2*SQRT(xSparse) - qa - qbr) / niR
          END IF

          xPrev = a
          DO 310 j=1,ni
            k=k+1

            xj = (qa + cf*j)**2
            IF (xj .GT. xSparse) THEN
               xj = w - (qbr + cf*(ni-j))**2
            END IF
            wj = xj - xPrev

            C(k) = Cprel * wj
            rj = Rprel * wj
            R(k) = rj + rPrev

            rPrev = rj
            xPrev = xj
310       CONTINUE
210     CONTINUE

      END IF

      R(k+1) = rPrev

      RETURN

      END