SUBROUTINE red_spec(wl,fl,av,yob,ired,iwavl,law_rd)
c     Reads a spectrum and calculates the new one with internal
c     reddening. 
c     law_rd=1 : the standard variation of absorption with wavelength 
c                is used (Allen, 76) - MW
c     law_rd=2 : the analytic fit of Seaton (79) curve by 
c                Fitzpatrick (86) - MW 
c     law_rd=3 : Fitzpatrick (86) - LMC
c     law_rd=4 : Prevot et al. (1984) and Bouchet et al. (1985) - SMC
c     law_rd=5 : Calzetti (astro-ph/9911459) - SB
c    
      INCLUDE 'dimension.dec'

      INTEGER ired,iwavl,law_rd
      REAL yob,alambda,wl,fl,wab,amab,aab,av,ff,ff11,ff12,slope,rv
      DIMENSION yob(mxwl),wl(mxwl),fl(mxwl),wab(30),amab(30),aab(30)

c ..... 1. Allen (1976) : wab = lambda in A, amab = k(lambda)/R, 
c          A(lambda) = k(lambda)*E(B-V) = amab*A_V ; R = 3.1

      IF (law_rd.eq.1) THEN
         rv=3.1
         wab(1)=1000.
         amab(1)=4.2
         wab(2)=1110.
         amab(2)=3.7
         wab(3)=1250.
         amab(3)=3.3
         wab(4)=1430.
         amab(4)=3.
         wab(5)=1670.
         amab(5)=2.7
         wab(6)=2000.
         amab(6)=2.8
         wab(7)=2220.
         amab(7)=2.9
         wab(8)=2500.
         amab(8)=2.30
         wab(9)=2850.
         amab(9)=1.97
         wab(10)=3330.
         amab(10)=1.69
         wab(11)=3650.
         amab(11)=1.58
         wab(12)=4000.
         amab(12)=1.45
         wab(13)=4400.
         amab(13)=1.32
         wab(14)=5000.
         amab(14)=1.13
         wab(15)=5530.
         amab(15)=1.00
         wab(16)=6700.
         amab(16)=0.74
         wab(17)=9000.
         amab(17)=0.46
         wab(18)=10000.
         amab(18)=0.38
         wab(19)=20000.
         amab(19)=0.11
         wab(20)=100000.
         amab(20)=0.0
         DO i=1,19
            aab(i)=amab(i)*av
         END DO
         DO ik=1,iwavl
            IF (wl(ik).le.wab(1)) then 
               ccc=(aab(1)-aab(2))/(wab(2)-wab(1))
               alambda=aab(1)+(wab(1)-wl(ik))*ccc
            ELSE IF(wl(ik).ge.wab(19))then 
               ccc=(aab(19)-aab(20))/(wab(20)-wab(19))
               alambda=aab(19)+(wab(19)-wl(ik))*ccc
            ELSE IF (wl(ik).gt.wab(1).and.wl(ik).lt.wab(19)) then
               DO i=1,19
                  IF (wl(ik).ge.wab(i).and.wl(ik).lt.wab(i+1)) then
                     ccc=(aab(i)-aab(i+1))/(wab(i+1)-wab(i))
                     alambda=aab(i)+(wab(i)-wl(ik))*ccc
                  END IF
               END DO
            END IF
            IF (alambda.lt.0.) alambda=0.
            IF (ired.eq.1) yob(ik)=fl(ik)*(10**(-0.4*alambda))
            IF (ired.eq.2) yob(ik)=fl(ik)*(10**(0.4*alambda))
         END DO
         RETURN

c ..... 2. Seaton (1979) - Fitzpatrick fit : 
c                       for lambda < 1200 A slope from 1100-1200
c                       for lambda > 3650 A points from Allen
c                       ff = E(lambda-V)/E(B-V) = k(lambda)-R
c                       A(lambda)=A_V*(ff/R +1) ; R = 3.1

      ELSE IF (law_rd.eq.2) THEN
         rv=3.1
         al0=4.595
         ga=1.051
         c1=-0.38
         c2=0.74
         c3=3.96
         c4=0.26
         call red(1100.,al0,ga,c1,c2,c3,c4,ff)
         ff11=ff
         call red(1200.,al0,ga,c1,c2,c3,c4,ff)
         ff12=ff
         slope=(ff12-ff11)/100.
         wab(11)=3650.
         amab(11)=1.58
         wab(12)=4000.
         amab(12)=1.45
         wab(13)=4400.
         amab(13)=1.32
         wab(14)=5000.
         amab(14)=1.13
         wab(15)=5530.
         amab(15)=1.00
         wab(16)=6700.
         amab(16)=0.74
         wab(17)=9000.
         amab(17)=0.46
         wab(18)=10000.
         amab(18)=0.38
         wab(19)=20000.
         amab(19)=0.11
         wab(20)=100000.
         amab(20)=0.0
         DO i=11,19
            aab(i)=amab(i)*av
         END DO
         DO ik=1,iwavl
            IF (wl(ik).ge.1200..and.wl(ik).le.3650.) then
               call red(wl(ik),al0,ga,c1,c2,c3,c4,ff)
               alambda=av*(ff/rv+1.) 
            ELSE IF (wl(ik).lt.1200.) then
               ff=ff11+(wl(ik)-1100.)*slope
               alambda=av*(ff/rv+1.)
            ELSE IF(wl(ik).ge.wab(19))then !extrapol at high wavelengths
               ccc=(aab(19)-aab(20))/(wab(20)-wab(19))
               alambda=aab(19)+(wab(19)-wl(ik))*ccc
            ELSE IF (wl(ik).gt.3650..and.wl(ik).lt.wab(19)) then
               DO i=11,19
                  IF (wl(ik).gt.wab(i).and.wl(ik).lt.wab(i+1)) then
                     ccc=(aab(i)-aab(i+1))/(wab(i+1)-wab(i))
                     alambda=aab(i)+(wab(i)-wl(ik))*ccc
                  END IF
               END DO
            END IF
            IF (alambda.lt.0.) alambda=0.
            IF (ired.eq.1) yob(ik)=fl(ik)*(10**(-0.4*alambda))
            IF (ired.eq.2) yob(ik)=fl(ik)*(10**(0.4*alambda))
         END DO
         RETURN
      
c ..... 3. Fitzpatrick (1986) : 
c                       for lambda < 1200 A slope derived from 1100-1200
c                       for lambda > 3330 A points from Allen
c                       ff = E(lambda-V)/E(B-V) = k(lambda)-R
c                       A(lambda) = A_V*(ff/R +1) ; R = 3.1

      ELSE IF (law_rd.eq.3) THEN
         rv=3.1
         al0=4.608
         ga=0.994
         c1=-0.69
         c2=0.89
         c3=2.55
         c4=0.50
         call red(1100.,al0,ga,c1,c2,c3,c4,ff)
         ff11=ff
         call red(1200.,al0,ga,c1,c2,c3,c4,ff)
         ff12=ff
         slope=(ff12-ff11)/100.
         wab(10)=3330.
         amab(10)=1.682
         wab(11)=3650.
         amab(11)=1.58
         wab(12)=4000.
         amab(12)=1.45
         wab(13)=4400.
         amab(13)=1.32
         wab(14)=5000.
         amab(14)=1.13
         wab(15)=5530.
         amab(15)=1.00
         wab(16)=6700.
         amab(16)=0.74
         wab(17)=9000.
         amab(17)=0.46
         wab(18)=10000.
         amab(18)=0.38
         wab(19)=20000.
         amab(19)=0.11
         wab(20)=100000.
         amab(20)=0.0
         DO i=10,19
            aab(i)=amab(i)*av
         END DO
         DO ik=1,iwavl
            IF (wl(ik).ge.1200..and.wl(ik).le.3330.) then
               call red(wl(ik),al0,ga,c1,c2,c3,c4,ff)
               alambda=av*(ff/rv+1.)
            ELSE IF (wl(ik).lt.1200.) then
               ff=ff11+(wl(ik)-1100.)*slope
               alambda=av*(ff/rv+1.)
            ELSE IF(wl(ik).ge.wab(19))then !extrapol at high wavelengths
               ccc=(aab(19)-aab(20))/(wab(20)-wab(19))
               alambda=aab(19)+(wab(19)-wl(ik))*ccc
            ELSE IF (wl(ik).gt.3330..and.wl(ik).lt.wab(19)) then
               DO i=10,19
                  IF (wl(ik).gt.wab(i).and.wl(ik).lt.wab(i+1)) then
                     ccc=(aab(i)-aab(i+1))/(wab(i+1)-wab(i))
                     alambda=aab(i)+(wab(i)-wl(ik))*ccc
                  END IF
               END DO
            END IF
            IF (alambda.lt.0.) alambda=0.
            IF (ired.eq.1) yob(ik)=fl(ik)*(10**(-0.4*alambda))
            IF (ired.eq.2) yob(ik)=fl(ik)*(10**(0.4*alambda))
         END DO
         RETURN

c ..... 4. Prevot (1984) and Bouchet (1985) :
c                        wab = lambda in A,  
c                        amab = E(lambda-V)/E(B-V) = k(lambda)-R
c                        A(lambda) = A_V*(amab/R +1)
c                        N.B. R = 2.72

         ELSE IF (law_rd.eq.4) THEN
            rv=2.72
            wab(1)=1275.
            amab(1)=13.54
            wab(2)=1330.
            amab(2)=12.52
            wab(3)=1385.
            amab(3)=11.51
            wab(4)=1435.
            amab(4)=10.80
            wab(5)=1490.
            amab(5)=9.84
            wab(6)=1545.
            amab(6)=9.28
            wab(7)=1595.
            amab(7)=9.06
            wab(8)=1647.
            amab(8)=8.49
            wab(9)=1700.
            amab(9)=8.01
            wab(10)=1755.
            amab(10)=7.71
            wab(11)=1810.
            amab(11)=7.17
            wab(12)=1860.
            amab(12)=6.90
            wab(13)=1910.
            amab(13)=6.76
            wab(14)=2000.
            amab(14)=6.38
            wab(15)=2115.
            amab(15)=5.85
            wab(16)=2220.
            amab(16)=5.30
            wab(17)=2335.
            amab(17)=4.53
            wab(18)=2445.
            amab(18)=4.24
            wab(19)=2550.
            amab(19)=3.91
            wab(20)=2665.
            amab(20)=3.49
            wab(21)=2778.
            amab(21)=3.15
            wab(22)=2890.
            amab(22)=3.00
            wab(23)=2995.
            amab(23)=2.65
            wab(24)=3105.
            amab(24)=2.29
            wab(25)=3704.
            amab(25)=1.81
            wab(26)=4255.
            amab(26)=1.00
            wab(27)=5291.
            amab(27)=0.00
            wab(28)=12500.
            amab(28)=-2.02
            wab(29)=16500.
            amab(29)=-2.36
            wab(30)=22000.
            amab(30)=-2.47

            DO i=1,30
               aab(i)=av*(amab(i)/rv+1.)
            END DO

            DO ik=1,iwavl
               IF (wl(ik).le.wab(1)) THEN !extrapol at short wavelengths
                  ccc=(aab(1)-aab(2))/(wab(2)-wab(1))
                  alambda=aab(1)+(wab(1)-wl(ik))*ccc
               ELSE IF(wl(ik).ge.wab(30))then !extrapol at high wavelengths
                  ccc=(aab(29)-aab(30))/(wab(30)-wab(29))
                  alambda=aab(29)+(wab(29)-wl(ik))*ccc  
               ELSE IF (wl(ik).gt.wab(1).and.wl(ik).lt.wab(30)) THEN
                  DO i=1,30
                     IF (wl(ik).gt.wab(i).and.wl(ik).lt.wab(i+1)) THEN
                        ccc=(aab(i)-aab(i+1))/(wab(i+1)-wab(i))
                        alambda=aab(i)+(wab(i)-wl(ik))*ccc
                     END IF
                  END DO
               END IF   
               IF (alambda.lt.0.) alambda=0.
               IF (ired.eq.1) yob(ik)=fl(ik)*(10**(-0.4*alambda))
               IF (ired.eq.2) yob(ik)=fl(ik)*(10**(0.4*alambda))
            END DO
            RETURN


c ..... Calzetti (1999) : 
c                   for lambda < 1200 A slope derived from 1100-1200 A
c                   for lambda > 22000A slope derived from 21900-22000 A
c                   ff = k(lambda); A(lambda) = A_V*k(lambda)/R ; R = 4.05

      ELSE IF (law_rd.eq.5) THEN
         rv=4.05
         p11=1/0.11
         ff11=2.659*(-2.156+1.509*p11-0.198*p11**2
     .        +0.011*p11**3)+rv
         p12=1/0.12
         ff12=2.659*(-2.156+1.509*p12-0.198*p12**2
     .        +0.011*p12**3)+rv
         slope1=(ff12-ff11)/100.
         ff99=2.659*(-1.857+1.040/2.19)+rv
         ff100=2.659*(-1.857+1.040/2.2)+rv
         slope2=(ff100-ff99)/100.
         DO ik=1,iwavl
            ala=wl(ik)*1.E-4    !wavelength in microns
            p=1./ala
            IF (ala.ge.0.63.and.ala.le.2.2) then
               ff=2.659*(-1.857+1.040*p)+rv
            ELSE IF (ala.ge.0.12.and.ala.lt.0.63) then
               ff=2.659*(-2.156+1.509*p-0.198*p**2+0.011*p**3)+rv
            ELSE IF (ala.lt.0.12) THEN
               ff=ff11+(wl(ik)-1100.)*slope1
            ELSE IF (ala.gt.2.2) THEN
               ff=ff99+(wl(ik)-21900.)*slope2
            END IF
            alambda=av*ff/rv
            IF (alambda.lt.0.) alambda=0.
            IF (ired.eq.1) yob(ik)=fl(ik)*(10**(-0.4*alambda))
            IF (ired.eq.2) yob(ik)=fl(ik)*(10**(0.4*alambda))
         END DO
         RETURN
      END IF

      WRITE(*,*) ' No reddening law is given.....'
      RETURN
      END

c ----------------------------------------------------------------------
      SUBROUTINE red(alan,al0,ga,c1,c2,c3,c4,ff)

c.......According to Fitzpatrick, 86

      ala=alan*1.E-4           !wavelength in microns
      p=1.0/ala
      IF (p.lt.5.9) c4=0.0
      t1=c3/(((p-(al0**2/p))**2)+ga*ga)
      t2=c4*(0.539*((p-5.9)**2)+0.0564*((p-5.9)**3))
      ff=c1+c2*p+t1+t2
      RETURN
      END