SUBROUTINE lyman_forest(wl,fl,iw,z)

c ..... Computes the Lyman forest mean absorption of an 
c       input spectrum, according to <D_A> and <D_B> evolution 
c       from Madau (1995).

      INTEGER mxwl,nw
      PARAMETER (nw=100)
      REAL da,db,z,wtemp(nw),w1,w2,wstep
      INCLUDE 'dimension.dec'
      DIMENSION wl(mxwl),fl(mxwl),ptau(nw)

      IF (z.lt.0.) z=0.
      w1=1050.*(1.+z)
      w2=1170.*(1.+z)
      wstep=(w2-w1)/FLOAT(nw)
      
      DO i=1,nw
         wtemp(i)=w1+(i-1)*wstep
         ptau(i)=EXP(-3.6E-3*(wtemp(i)/1216.)**3.46)
      END DO

      da=(1./(120.*(1.+z)))*trapz1(wtemp,ptau,nw) ! da = 1-<D_A>

      w1=920.*(1.+z)
      w2=1015.*(1.+z)
      wstep=(w2-w1)/FLOAT(nw)

      DO i=1,nw
         wtemp(i)=w1+(i-1)*wstep
         ptau(i)=EXP(-1.7E-3*(wtemp(i)/1026.)**3.46
     .        -1.2E-3*(wtemp(i)/972.5)**3.46
     .        -9.3E-4*(wtemp(i)/950.)**3.46)
      END DO

      db=(1./(95.*(1.+z)))*trapz1(wtemp,ptau,nw) ! db = 1-<D_B>

c More Lyman blanketing
c      da=0.5*da
c      db=0.5*db
c Less Lyman blanketing
c      da=2.0*da
c      db=2.0*db

      IF (da.gt.1.) da=1.
      IF (db.gt.1.) db=1.
      IF (da.lt.0.) da=0.
      IF (db.lt.0.) db=0.

c        correcting spectrum:

      DO i=1,iw
         IF (wl(i).lt.1216.0.and.wl(i).gt.1026.0) THEN !between Lya and Lyb
            fl(i)=fl(i)*da
         END IF
         IF (wl(i).le.1026.0) THEN               !between Lyb and 912A break
            fl(i)=fl(i)*db
         END IF
      END DO

      RETURN
      END