      PROGRAM zphoto
c     Version: 11/12/2000
      IMPLICIT NONE
c
      INTEGER i,ii,iii,j,k,l,ks,n_files,icomp,ncomp,kfil,iw_vega,iw,
     .     n_filter,iob,nob,id,iout_s,iout_m,law_rd,nfil,iroot,matrix,
     .     iw_ob,i_typ,i_red,i_age,nz,itype,iage,izmin,itypem,iagem,
     .     iredm,i_break,npoiz,n_sol,iwav,npoiz1,npoiz2,last,nred,ired,
     .     law_nd,i_lyman,mxz,izsol,mxage,mxtyp,mxfil,mxred,mxwl,ndeg,
     .     magtyp,iunit_sed,i1,i2,icat,iblank,input,nparam,ipass,
     .     izmin1,iagemax,iread,jread,nf,np,mxft,mxwf,nfilt_mabs,
     .     iz_clus,n_fil_wa,n_fil_w,itypea,ireda,iagea,iztyp,idum,
     .     izphot,i_sol,iout,iz,iforbid,n_head
c
      REAL w_vega,f_vega,T_red,f_red,f,w,z1,z2,amag,err_p,mag_nd,
     .     err_nd,err_min,avmin,avmax,averr,wl,surf,c_corr,val,p,
     .     flux,err_fl,dz,zmax,zmin,zz,z_min_array,filter,z,y,zop,av,
     .     q_pro,al_lyman,aw_limite,dz_f,zmin_f,zmax_f,age,fl,qq,
     .     poids,poids_z,q_pro1,q_pro2,q_proz,chi2_b1,chi2_b2,z_b,
     .     z_b1,z_b2,dd,z_poids,qq_p,w0,wmin1,wmax1,wsyn,ysyn,d1,d2,
     .     yob,zz1,zz2,t,vega,pthresh,conv_mag,bp1,cputime,time_tot,
     .     part_time,f_nu,f_jy,err_jy,deg,pchi2,qchi2,zero,chi2,
     .     r_chi2,amabs_1,amabs_2,sauce,amagb,bmag,dif_w,dif_w_min,
     .     fl_temp,mag_temp,func,prob,err,x1,x2,ww,h0,om,ov,ageuni,tz,
     .     r,zeropoint,d_mod,color,z_clus,amagb_clus,p_clus,dummy,
     .     dummyz,zphot,chi2sol,qq2,fl_nd,fl1,ebv_mw,av_mw,fl_ext,
     .     rlin,chi2min,bb_array,c,fnu_vega,rnu_vega,wrnu,rrnu,
     .     f1,f2,ffnu,trapz1
      REAL*8 bm,wt,sxx,sxy,b,bb,ddummy
c
      PARAMETER (err=1.E-8,nparam=35,mxft=300,mxwf=200)
      INCLUDE 'dimension.dec'

      LOGICAL header
      CHARACTER header_line*5000
      CHARACTER NAMECH(mxtyp)*160,sed_abc(mxtyp)*2,line*240
      CHARACTER*160 spectra,amatrix,filedef,filein,filt_file,
     .     cat_in,spectra1,spectra2,spectra3,spectra4,spectra5,
     .     spec_vega,filtresfile,line_sed,filterout,copyright1,
     .     copyright2,copyright3
      CHARACTER clog*1,ccat*1,cz*1,csed1*1,csed2*1,cage*1,cout_m*1,
     .     cout_s*1,root*8,date*24,vv*5
      CHARACTER fmt1*90,fmt2*90,fmt3*90,fmt4*90
      CHARACTER*40 cname,param(nparam)
c
      DIMENSION input(nparam),kfil(10)
      DIMENSION t_red(225),t(225),f_red(mxwl,225,mxtyp),ks(mxtyp)
      DIMENSION w_Vega(4000),f_vega(4000),vega(mxfil)
      DIMENSION bp1(0:mxfil),conv_mag(mxfil)
      DIMENSION c_corr(mxfil),fl_nd(mxfil),fl1(mxfil)
      DIMENSION f_jy(mxfil),err_jy(mxfil),fl_temp(mxfil),mag_temp(mxfil)
      DIMENSION amag(mxfil),err_p(mxfil),mag_nd(mxfil),err_nd(mxfil)
      DIMENSION law_nd(mxfil),aw_limite(mxfil)
      DIMENSION ysyn(mxwl),wsyn(mxwl),ww(mxwl),yob(mxwl),y(mxwl)
      DIMENSION z_min_array(mxz),i_typ(mxz),i_age(mxz),i_red(mxz)
      DIMENSION wl(0:mxfil),flux(0:mxfil,0:mxz,mxage,mxtyp,mxred)
      DIMENSION q_pro(mxz),q_pro1(mxz),q_pro2(mxz),q_proz(mxz)
      DIMENSION z(0:mxz),z_b(mxz),z_b1(mxz),z_b2(mxz)
      DIMENSION chi2_b1(mxz),chi2_b2(mxz),iforbid(mxz),bb_array(mxz)
      DIMENSION izphot(50),zphot(50),bb(50),chi2sol(50)
      DIMENSION val(3),zz1(3),zz2(3),p(3),zeropoint(mxfil),n_fil_w(mxz)
      DIMENSION fnu_vega(4000),rnu_vega(4000)
      DIMENSION wrnu(mxwf),rrnu(mxwf),ffnu(mxwf)
c
      EXTERNAL cputime,rlin,filter,pchi2,qchi2,zero,func,tz,trapz1
      EXTERNAL d_mod,color
      COMMON/et1/surf(0:mxfil)
      COMMON/et2/W(mxwl,mxtyp),F(mxwl,mxage,mxtyp),fl(mxfil),
     .     err_fl(mxfil)
      COMMON/et3/law_rd
      COMMON/et4/n_filter,nfil(0:mxfil),iw(mxtyp),i_lyman(mxfil)
      COMMON/et5/bm
      COMMON/compose/icomp
      COMMON/nfree/deg,prob
      COMMON/filterres/filtresfile
      COMMON/c_filter/iread,jread,nf,np(mxft),r(mxwf,mxft,2)
c_________________________________________________________________________

      vv='1.1'
      copyright1='    hyperz photometric redshift code - Version '//vv
     .     //' (December 11, 2000)' 
      copyright2=' by Micol Bolzonella, Roser Pello'','
     .     //' Joan-Marc Miralles'
      copyright3=' More information at '
     .     //'http://webast.ast.obs-mip.fr/hyperz/'
      WRITE(*,*)
      WRITE(*,'(a80)') copyright1
      WRITE(*,*)
      WRITE(*,'(a80)') copyright2
      WRITE(*,'(a80)') copyright3
      WRITE(*,*)

      jread=0                   ! Control on reading A0VSED
      iread=0                   ! Control on reading Filter Library
      c=2.9979E+18              ! speed of light in AA/s

c ..... Probabilities for the confidence intervals
      p(1)=0.99
      p(2)=0.90
      p(3)=0.68
      x1=0.
      x2=200.
      ndeg=1
      deg=REAL(ndeg)
      DO i=1,3
         prob=p(i)
         val(i)=zero(func,x1,x2,err)  ! chi2 for P = p(i) 
      END DO

c....... limits in absolute magnitude

      amabs_1=-100.0
      amabs_2=100.0
      z_clus=-1.

 11   filein='zphot.param'
      last=INDEX(filein,'  ')-1
      WRITE(*,*) 'Name of parameter file'
      WRITE(*,*) '[default = ',filein(1:last),']'
      READ(*,'(A)') line
      IF (line.ne.' ') THEN
         READ(line,'(A)',ERR=11) filein
      END IF
      time_tot=cputime()
c ------------------------------------------------------------------------
c                       Reading parameters file 
c ------------------------------------------------------------------------

      param(1)='A0VSED'
      param(2)='FILTERS_RES'
      param(3)='FILTERS_FILE'
      param(4)='TEMPLATES_FILE'
      param(5)='FILT_M_ABS'
      param(6)='ERR_MAG_MIN'
      param(7)='Z_MIN'
      param(8)='Z_MAX'
      param(9)='Z_STEP'
      param(10)='ZSTEP_TYPE'
      param(11)='CATALOG_TYPE'
      param(12)='CATALOG_FILE'
      param(13)='MAG_TYPE'
      param(14)='REDDENING_LAW'
      param(15)='AV_MIN'
      param(16)='AV_MAX'
      param(17)='AV_STEP'
      param(18)='AGE_CHECK'
      param(19)='PROB_THRESH'
      param(20)='OUTPUT_FILE'
      param(21)='OUTPUT_TYPE'
      param(22)='SED_OBS_FILE'
      param(23)='SED_TEMP_FILE'
      param(24)='LOGPHOT_FILE'
      param(25)='CATPHOT_FILE'
      param(26)='ZPHOT_FILE'
      param(27)='H0'
      param(28)='OMEGA_M'
      param(29)='OMEGA_V'
      param(30)='Z_CLUSTER'
      param(31)='M_ABS_MIN'
      param(32)='M_ABS_MAX'
      param(33)='MATRIX'
      param(34)='SPECTRUM'
      param(35)='EBV_MW'

      DO i=1,nparam
         input(i)=0
      END DO
      ebv_mw=0.

      OPEN(77,FILE=filein,STATUS='old',ERR=11)
      WRITE(*,*) 'Reading parameter file :'

 7    READ(77,'(A)',ERR=10,END=10) line
      DO i=1,80
         IF (line(i:i).ne.' ') THEN
            line=line(i:)
            GO TO 8
         END IF
      END DO
      GO TO 7  ! blank line
 8    IF (line(1:1).eq.'#') GO TO 7 ! comment line
      iblank=INDEX(line,' ')
      cname=line(1:iblank-1) ! variable name
      line=line(iblank:)
      DO i=1,80
         IF (line(i:i).ne.' ') THEN
            line=line(i:)
            IF (line(1:1).eq.'#') THEN
               WRITE(*,*) 'Missed value of variable ', cname
               STOP
            END IF
            GO TO 9
         END IF
      END DO
      WRITE(*,*) 'Missed value of variable ', cname
      STOP
 9    iblank=INDEX(line,' ')
      line=line(1:iblank-1) ! part of string containing value of variable
c ..... File containing Vega spectrum
      IF (cname.eq.param(1)) THEN
         READ(line,'(A)') spec_vega
         input(1)=1
         GO TO 7
c ..... File FILTER.RES
      ELSE IF (cname.eq.param(2)) THEN
         READ(line,'(A)') filtresfile
         input(2)=1
         GO TO 7
c ..... File containing filter parameters and instructions for 
c       non detected objects (nd) or out of file (oof) for each filter
c       0 -> flux = -1 (filter not considered) for either nd or oof
c       1 -> flux = -1 for oof, flux = 0, err_flux = flux_nd for nd
c       2 -> flux = -1 for oof, flux = flux_nd/2, err_flux = flux_nd/2 for nd
c       3 -> flux = -1 for oof, mag = mag_nd,  err = err_nd for nd
      ELSE IF (cname.eq.param(3)) THEN
         READ(line,'(A)') filt_file
         input(3)=1
         GO TO 7
c ..... File containing synthetic spectra filenames
      ELSE IF (cname.eq.param(4)) THEN
         READ(line,'(A)') filedef
         input(4)=1
         GO TO 7
c ..... Filter for the absolute magnitude
      ELSE IF (cname.eq.param(5)) THEN
         READ(line,*) nfilt_mabs
         input(5)=1
         GO TO 7
c ..... Minimum error in magnitude (if error < err_min => err = err_min)
      ELSE IF (cname.eq.param(6)) THEN
         READ(line,*) err_min
         input(6)=1
         GO TO 7
c ..... Redshift range (utilized only if icat=0) and step for scan
      ELSE IF (cname.eq.param(7)) THEN
         READ(line,*) zmin
         input(7)=1
         GO TO 7
      ELSE IF (cname.eq.param(8)) THEN
         READ(line,*) zmax
         input(8)=1
         GO TO 7
      ELSE IF (cname.eq.param(9)) THEN
         READ(line,*) dz
         input(9)=1
         GO TO 7
      ELSE IF (cname.eq.param(10)) THEN
         READ(line,*) iztyp
         IF (iztyp.ne.0.and.iztyp.ne.1) THEN
            WRITE(*,*) 'ZSTEP_TYPE must be 0 or 1'
            STOP
         END IF
         input(10)=1
         GO TO 7
c ..... Catalogue format 
c       Format of objects catalogue: 
c       N,[z_min,z_max,](mag(i),i=1,n),(err(i),i=1,n)
      ELSE IF (cname.eq.param(11)) THEN
         READ(line,*) icat
         IF (icat.ne.0.and.icat.ne.1) THEN
            WRITE(*,*) 'CATALOG_TYPE must be 0 or 1'
            STOP
         END IF
         input(11)=1
         GO TO 7
c ..... File containing photometric catalogue
      ELSE IF (cname.eq.param(12)) THEN
         READ(line,'(A)') cat_in
         input(12)=1
         GO TO 7
c ..... Magnitude type:
c       0 -> standard magnitudes
c       1 -> AB magnitudes
      ELSE IF (cname.eq.param(13)) THEN
         READ(line,*) magtyp
         IF (magtyp.ne.0.and.magtyp.ne.1) THEN
            WRITE(*,*) 'MAG_TYPE must be 0 or 1'
            STOP
         END IF
         input(13)=1
         GO TO 7
c ..... Reddening:
c       0 -> No reddening
c       1 -> Allen, 76, standard MW
c       2 -> Seaton, 79, standard MW
c       3 -> Fitzpatrick, 1986, LMC
c       4 -> Prevot, 1984 + Bouchet, 1985, SMC
c       5 -> Calzetti, astro-ph/9911459
      ELSE IF (cname.eq.param(14)) THEN
         READ(line,*) law_rd
         IF (law_rd.ne.0.and.law_rd.ne.1.and.law_rd.ne.2.and.
     .        law_rd.ne.3.and.law_rd.ne.4.and.law_rd.ne.5) THEN
            WRITE(*,*) 'REDDENING_LAW must be 0 -- 5'
            STOP
         END IF
         input(14)=1
         GO TO 7
c ..... Av_max, Av_min, Av_err
      ELSE IF (cname.eq.param(15)) THEN
         READ(line,*) avmin
         input(15)=1
         GO TO 7
      ELSE IF (cname.eq.param(16)) THEN
         READ(line,*) avmax
         input(16)=1
         GO TO 7
      ELSE IF (cname.eq.param(17)) THEN
         READ(line,*) averr
         input(17)=1
         GO TO 7
c ..... Check the age of the object < universe age
      ELSE IF (cname.eq.param(18)) THEN
         cage=line(1:1)
         input(18)=1
         GO TO 7
c ..... Probability threshold for the study of secondary maxima
      ELSE IF (cname.eq.param(19)) THEN
         READ(line,*) pthresh
         input(19)=1
         GO TO 7
c ..... Output files: no extension (file = *)
c       *.cat_phot -> catalogue
c       *.log_phot -> log file
c       *.sed      -> SED
c       *.z_phot   -> catalogue with only first maximum
      ELSE IF (cname.eq.param(20)) THEN
         READ(line,'(A)') spectra
         iblank=INDEX(spectra,' ')
         spectra1=spectra(1:iblank-1)//'.cat_phot'
         spectra2=spectra(1:iblank-1)//'.log_phot'
         spectra3=spectra(1:iblank-1)//'.obs_sed'
         spectra4=spectra(1:iblank-1)//'.z_phot'
         spectra5=spectra(1:iblank-1)//'.temp_sed'
         filterout=spectra(1:iblank-1)//'.out_phot'
         input(20)=1
         GO TO 7
c ..... Output SED in erg/s/cm^2/A (0), microJansky (1) or mag AB (2)
      ELSE IF (cname.eq.param(21)) THEN
         READ(line,*) iunit_sed
         IF (iunit_sed.ne.0.and.iunit_sed.ne.1.and.iunit_sed.ne.2) THEN
            WRITE(*,*) 'OUTPUT_TYPE must be 0, 1 or 2'
            STOP
         END IF
         input(21)=1
         GO TO 7
c ..... File .obs_sed
      ELSE IF (cname.eq.param(22)) THEN
         csed1=line(1:1)
         input(22)=1
         GO TO 7
c ..... File .temp_sed
      ELSE IF (cname.eq.param(23)) THEN
         csed2=line(1:1)
         input(23)=1
         GO TO 7
c ..... File .log_phot
      ELSE IF (cname.eq.param(24)) THEN
         clog=line(1:1)
         input(24)=1
         GO TO 7
c ..... File .cat_phot
      ELSE IF (cname.eq.param(25)) THEN
         ccat=line(1:1)
         input(25)=1
         GO TO 7
c ..... File .z_phot
      ELSE IF (cname.eq.param(26)) THEN
         cz=line(1:1)
         input(26)=1
         GO TO 7
c ..... Cosmological parameters to compute the universe age and M_abs 
      ELSE IF (cname.eq.param(27)) THEN
         READ(line,*) h0
         input(27)=1
         GO TO 7
      ELSE IF (cname.eq.param(28)) THEN
         READ(line,*) om
         input(28)=1
         GO TO 7
      ELSE IF (cname.eq.param(29)) THEN
         READ(line,*) ov
         input(29)=1
         GO TO 7
c ..... Cluster redshift (optional)
      ELSE IF (cname.eq.param(30)) THEN
         READ(line,*) z_clus
         input(30)=1
         GO TO 7
c ..... Limits for absolute magnitude
      ELSE IF (cname.eq.param(31)) THEN
         READ(line,*) amabs_1
         input(31)=1
         GO TO 7
      ELSE IF (cname.eq.param(32)) THEN
         READ(line,*) amabs_2
         input(32)=1
         GO TO 7
c ..... Matrix/obj (optional)
      ELSE IF (cname.eq.param(33)) THEN
         cout_m=line(1:1)
         input(33)=1
         GO TO 7
c ..... Spectrum/obj (optional)
      ELSE IF (cname.eq.param(34)) THEN
         cout_s=line(1:1)
         input(34)=1
         GO TO 7
c .... galactic dereddening
      ELSE IF (cname.eq.param(35)) THEN
         READ(line,*) ebv_mw
         input(35)=1
         GO TO 7
      ELSE 
         WRITE(*,'(a,1x,a)') 'Unknown parameter (ignored)',cname
         GO TO 7
      END IF

 10   CONTINUE
      CLOSE(77)
      DO i=1,nparam-6
         IF (input(i).eq.0) THEN
            WRITE(*,*) 'Missed parameter ', param(i)
            STOP
         END IF
      END DO
c ..... Control of optional parameters      
      IF (input(nparam-5).eq.1.and.icat.eq.0) THEN
         IF (z_clus.lt.zmin.or.z_clus.gt.zmax) THEN
            WRITE(*,*) 'Z_CLUSTER must be in the interval [Z_MIN,Z_MAX]'
            STOP
         END IF
      END IF
      iout_m=1
      iout_s=1
      IF (input(nparam-2).eq.1.and.
     .     cout_m.eq.'y'.or.cout_m.eq.'Y') iout_m=2
      IF (input(nparam-1).eq.1.and.
     .     cout_s.eq.'y'.or.cout_s.eq.'Y') iout_s=2

      av_mw=3.1*ebv_mw

      IF (law_rd.eq.0) THEN
         avmin=0.0
         avmax=0.0
         averr=1.0
      END IF
      nred=INT(ABS(avmax-avmin)/averr)+1
      IF (nred.gt.mxred) THEN
         WRITE(*,'(1x,a,i3,a)') 
     .        'Number of reddening steps exceeding the maximum [max=',
     .        mxred,']'
         STOP
      END IF

      OPEN(10,FILE=filedef,STATUS='old',ERR=78)
      iii=1
 2230 READ(10,'(A)',ERR=2220,END=2220) line_sed
      IF (line_sed(1:1).eq.'#'.or.line_sed.eq.' ') GO TO 2230
      last=INDEX(line_sed,' ')-1
      NAMECH(iii)=line_sed(1:last)
      line_sed=line_sed(last+1:)
      DO i=1,80
         IF (line_sed(i:i).ne.' ') THEN
            line_sed=line_sed(i:)
            GO TO 2221
         END IF
      END DO
      WRITE(*,*) 'Missed type of SED file ', NAMECH(iii)
      STOP
 2221 sed_abc(iii)=line_sed(1:2)
      iii=iii+1
      GO TO 2230
 2220 n_files=iii-1
      IF (n_files.gt.mxtyp) THEN
         WRITE(*,'(1x,a,i3,a)') 
     .        'Number of files exceeding the maximum [max=',mxtyp,']'
         STOP
      END IF
      WRITE(*,*) 'There are ',n_files,' synthetic spectra'
      CLOSE(10)

      CALL filter0
      OPEN(10,FILE=filt_file,STATUS='old',ERR=80)
      i=0
 134  CONTINUE
      READ(10,'(A)',ERR=135,END=135) line
      IF (line(1:1).eq.'#'.or.line.eq.' ') GO TO 134
      i=i+1
      ncomp=1
      DO k=1,80
         IF (line(k:k).eq.'&') THEN
            line(k:k)=' '
            ncomp=ncomp+1
         END IF
      END DO
      IF (ncomp.eq.1) THEN 
         READ(line,*,ERR=82) nfil(i),c_corr(i),law_nd(i),
     .        mag_nd(i),err_nd(i)
         IF (law_nd(i).ne.0.and.law_nd(i).ne.1.and.law_nd(i).ne.2.and.
     .        law_nd(i).ne.3) THEN
            WRITE(*,*) 'Non detection law must be 0, 1, 2 or 3'
            WRITE(*,*) 'Check filter ', nfil(i)
            STOP
         END IF
         IF (nfil(i).gt.icomp) THEN
            WRITE(*,1998) nfil(i),icomp
            STOP
         ELSE
c ..... simple filter
            kfil(1)=nfil(i)
            CALL compos_new(1,kfil,wl(i),surf(i),bp1(i))
            GO TO 134
         END IF
      ELSE 
         READ(line,*,ERR=82) (kfil(j),j=1,ncomp),c_corr(i),law_nd(i),
     .                 mag_nd(i),err_nd(i)
         IF (law_nd(i).ne.0.and.law_nd(i).ne.1.and.law_nd(i).ne.2.and.
     .        law_nd(i).ne.3) THEN
            WRITE(*,*) 'Non detection law must be 0, 1, 2 or 3'
            WRITE(*,*) 'Check filter ', nfil(i)
            STOP
         END IF
c ..... composite filter
         icomp=icomp+1
         nfil(i)=icomp
         CALL compos_new(ncomp,kfil,wl(i),surf(i),bp1(i))
         GO TO 134
      END IF
 1998 FORMAT('Filter number ',I4, ' out of range!',/,
     .       'Maximum = ',I4)

 135  n_filter=i
      IF (n_filter.gt.mxfil) THEN
         WRITE(*,'(1x,a,i3,a)') 
     .        'Number of filters exceeding the maximum [max=',mxfil,']'
         STOP
      END IF
      WRITE(*,*) 'There are ',n_filter,' filters'
      nfil(0)=nfilt_mabs ! filter for absolute magnitude estimates
      CALL compos_new(1,nfilt_mabs,wl(0),surf(0),bp1(0))
      CLOSE(10)

      OPEN(10,FILE=cat_in,STATUS='old',ERR=79)

c ..... Reading the header (thanks to V. Le Brun)
      header = .FALSE.
      READ(10,'(A)',ERR=79) header_line
      IF (header_line(1:1).eq.'#'.or.header_line.eq.' ') THEN
         header= .TRUE.
         i=1
 44      READ(10,'(A)',ERR=79) header_line
         IF (header_line(1:1).eq.'#'.or.header_line.eq.' ') THEN
            i=i+1
            GO TO 44
         ELSE
            n_head=i
         END IF
      END IF
      CLOSE(10)
      OPEN(10,FILE=cat_in,STATUS='old',ERR=79)
c ..... Skipping the header (thanks to V. Le Brun)
      IF (header) THEN
         DO i=1,n_head
            READ(10,'(A)') header_line
         END DO
      END IF

      IF (icat.eq.0) THEN
         iob=0
 1340    CONTINUE
         iob=iob+1
         READ(10,*,ERR=1350,END=1350) id,(amag(i),i=1,n_filter),
     .        (err_p(i),i=1,n_filter)
         z1=zmin
         z2=zmax
         GO TO 1340
      ELSE IF (icat.eq.1) THEN
         iob=0
 1341    CONTINUE
         iob=iob+1
         READ(10,*,ERR=1350,END=1350) id,z1,z2,(amag(i),i=1,n_filter),
     .        (err_p(i),i=1,n_filter)
         IF(z1.lt.zmin) zmin=z1
         IF(z2.gt.zmax) zmax=z2
         GO TO 1341
      ELSE
         GO TO 79
      END IF
 1350 nob=iob-1
      WRITE(*,*) 'There are ',nob,' objects in the catalogue'
      CLOSE(10)

      IF (iout_m.gt.1.and.nob.gt.100) THEN
         iout=0
         WRITE(*,'(a,i6,a)') ' You chose output matrix for all ',
     .        nob, ' objects' 
         WRITE(*,'(a)') ' Are you sure?'
         WRITE(*,*) 'Type 0 to continue with matrix, 1 without matrix'
         WRITE(*,*) '[<return>=1]'
         READ(*,'(a)') line 
         IF (line.ne.' ') THEN
            READ(line,*) iout
            IF (iout.eq.0) iout_m=2
         END IF
      END IF
c      z(0)=0.
c      z(1)=zmin
      z(0)=zmin
      z(1)=zmin+dz
      i=1
 88   i=i+1
      IF (iztyp.eq.0) THEN
         z(i)=z(i-1)+dz
      ELSE IF (iztyp.eq.1) THEN
         z(i)=z(i-1)+(1.+z(i-1))*dz
      END IF
      IF (zmax-z(i).gt.1.e-3) GO TO 88
      nz=i
      IF (nz.gt.mxz) THEN
         WRITE(*,'(1x,a,i4,a)') 
     .        'Number of redshift steps exceeding the maximum [max=',
     .        mxz,']'
         WRITE(*,*) 'Choose a larger Z_STEP or a smaller Z_MAX'
         STOP
      END IF

c -----------------------------------------------------------------------
c                       Reading synthetic spectra 
c -----------------------------------------------------------------------

c ..... Vega

      WRITE (*,'(1x,a)') 'Reading reference sed:'
      OPEN(10,FILE=spec_vega,STATUS='old')
      READ(10,'(A)') line
      WRITE(*,'(1x,a)') line (1:78)
      iii=1
 223  READ(10,*,ERR=222,END=222) w_vega(iii),f_vega(iii)
      iii=iii+1
      GO TO 223
 222  iw_vega=iii-1
      CLOSE(10)
      WRITE(*,*) 'There are ',iw_vega,' points in Vega'
      jread=iw_vega

c ..... SED

      WRITE(*,*) 'Reading synthetic spectra :'
      DO i=1,n_files
         WRITE(*,2005) i, NAMECH(i)
c GISSEL's templates
         IF (sed_abc(i).eq.'BC') THEN
            OPEN(1,FILE=NAMECH(i),STATUS='old',
     .           FORM='unformatted',ERR=25)
            READ (1) KS(i),(T_red(ii),ii=1,KS(i)) ! KS = n. temporal intervals
            READ (1) IW(i),(W(ii,i),ii=1,IW(i)) !Number of wavelengths
            DO j=1,KS(i)
               READ (1) IW(i),(F_red(ii,j,i),ii=1,IW(i))
            ENDDO
            CLOSE(1)
c extended CWW's templates
         ELSE IF (sed_abc(i).eq.'AS'.or.sed_abc(i).eq.'A ')THEN
            OPEN(1,FILE=NAMECH(i),STATUS='old',
     .           FORM='formatted',ERR=25)
            ii=1
 33         READ(1,*,ERR=34,END=34) W(ii,i),F_red(ii,2,i)
            ii=ii+1
            GO TO 33
 34         iw(i)=ii-1
            CLOSE(1)
         ELSE
            WRITE(*,*) 'Don''t know how to read :'
            WRITE(*,'(A)') NAMECH(i)
            STOP
         END IF
      END DO
 2005 FORMAT(' Type ', I2, ' : ', A60)
c -----------------------------------------------------------------------
c                           Rebinning in age
c -----------------------------------------------------------------------

      DO i=1,n_files
         last=INDEX(NAMECH(i),' ')-1
         NAMECH(i)=NAMECH(i)(1:last)
c templates with 221 ages -> 51 after the rebinning
         IF (sed_abc(i).eq.'BC') THEN
            DO k=1,7            ! age=0 to 0.01 Gyr
               j=k*10
               t(k)=t_red(j)
               DO ii=1,iw(i)
                  F(ii,k,i)=F_red(ii,j,i)
                  IF (w(ii,i).lt.912.) F(ii,k,i)=0.
               END DO
            END DO
            DO k=8,34           ! age=0.01 to 3.0 Gyr
               j=70+(k-7)*3
               t(k)=t_red(j)
               DO ii=1,iw(i)
                  F(ii,k,i)=F_red(ii,j,i)
                  IF (w(ii,i).lt.912.) F(ii,k,i)=0.
              END DO
            END DO
            DO k=35,51          ! age=3.0 to 20.0 Gyr
               j=151+(k-34)*4
               t(k)=t_red(j)
               DO ii=1,iw(i)
                  F(ii,k,i)=F_red(ii,j,i)
                  IF (w(ii,i).lt.912.) F(ii,k,i)=0.
               END DO
            END DO
            KS(i)=51
c templates with only 1 age 
         ELSE IF (sed_abc(i).eq.'AS'.or.sed_abc(i).eq.'A ')THEN
            DO ii=1,iw(i)
               F(ii,2,i)=F_red(ii,2,i)
               IF (w(ii,i).lt.912.) F(ii,2,i)=0.
            END DO
            KS(i)=2
         END IF
      END DO
c -----------------------------------------------------------------------
c                    Filters and .out_phot file
c -----------------------------------------------------------------------

      OPEN(18,FILE=filterout,STATUS='unknown')
      WRITE(18,'(a)')
     .     '##########################################################'
      WRITE(18,'(a,1x,a)')
     .     '# This is the output file of hyperz - Version '//vv
      WRITE(18,'(a)')
     .     '##########################################################'
      CALL fdate(date)
c      CALL fdate_(date)   ! IBM
      WRITE(18,'(a)') '# Date: '//date
      WRITE(18,'(a)') '# Photometric catalogue : '
      last=INDEX(cat_in,' ')-1
      WRITE(18,'(2x,a)')  cat_in(1:last)
      WRITE(18,'(a)') 
     .     '# Templates SEDs:'
      DO itype=1,n_files
         last=INDEX(NAMECH(itype),' ')-1
         WRITE(18,'(2x,i2,2x,a,2x,a)') itype,NAMECH(itype)(1:last),
     .     sed_abc(itype)
      END DO
      WRITE(18,'(a)') '# Characteristics of filters'
      WRITE(18,'(a)') 
     .     '#  n     wl_eff     surface   bandpass1  conv_AB mag_lim nd'
      zz=0.

c ..... Conversion from AB to standard magnitudes
c compute f_Vega(nu)
      DO i=1,iw_vega
c         fnu_vega(i)=f_vega(i)*2.E-17*w_vega(i)*w_vega(i)/c
         fnu_vega(i)=f_vega(i)*2.*w_vega(i)*w_vega(i)/2.9979
         rnu_vega(i)=c/w_vega(i)
      END DO
      DO k=1,n_filter
         vega(k)=filter(nfil(k),w_vega,f_vega,iw_vega,zz) ! unit 2.E-17
         zeropoint(k)=-color(nfil(k),nfilt_mabs,
     .                       w_vega,f_vega,iw_vega,zz)
c         f_nu=(vega(k)*2.E-17/surf(k))*wl(k)**2/(2.9979E+18)
         f_nu=(vega(k)*2./surf(k))*wl(k)**2/2.9979
c         conv_mag(k)=-2.5*(-35.+LOG10(f_nu))-48.60 ! approx value

c transmission function R(nu): R(nu) d(nu)= R(lambda) d(lambda)
         DO i=1,np(nfil(k))
            wrnu(i)=c/r(i,nfil(k),1) ! frequency
            rrnu(i)=r(i,nfil(k),2)*c/(wrnu(i)*wrnu(i)) ! R(nu)
         END DO
c R(nu)*f_vega(nu)
         DO i=1,np(nfil(k))
            l=0
            ffnu(i)=rrnu(i)*rlin(wrnu(i),rnu_vega,fnu_vega,iw_vega,l)
         END DO
c integrals with trapezoids
         f1=trapz1(wrnu,ffnu,np(nfil(k)))
         f2=trapz1(wrnu,rrnu,np(nfil(k)))
         conv_mag(k)=-2.5*(-35.+LOG10(f1/f2))-48.60
c ..... Write file 'filters.dat'
         WRITE(18,1999) nfil(k),wl(k),surf(k),bp1(k),conv_mag(k),
     .        mag_nd(k),law_nd(k)
c ..... Wavelength limit for the break in each filter 
         aw_limite(k)=r(np(nfil(k))-1,nfil(k),1)
c ..... Flux corresponding to the limiting magnitude (N.B. unit 2.E-17)
         IF (magtyp.eq.0) fl_nd(k)=vega(k)*
     .           10**(-0.4*(mag_nd(k)+c_corr(k)))/surf(k) 
         IF (magtyp.eq.1) fl_nd(k)=10**(-0.4*(mag_nd(k)+48.6))*
     .        (2.9979E+18/(wl(k)**2))/2.E-17  
         mag_nd(k)=-2.5*LOG10(fl_nd(k)*surf(k)/vega(k)) ! Vega limiting mag.
      END DO
 1999 FORMAT(2X,I3,1X,F11.3,1X,2(F10.3,1X),1X,F7.3,1X,F7.3,1X,I2)
      WRITE(18,'(a)') '#    zmin      zmax     zstep'
      WRITE(18,'(4X,3(F5.2,5X))') zmin,zmax,dz
      WRITE(18,'(a)') '#    AVmin     AVmax    AVstep  red_law'
      WRITE(18,'(4X,3(F5.2,5X),I2)') avmin,avmax,averr,law_rd
      WRITE(18,'(a)') '#    Mmin      Mmax     Mfilter'
      WRITE(18,'(4X,2(F7.2,2X),3X,I3)') amabs_1,amabs_2,nfilt_mabs
      WRITE(18,'(a)') '#    Omega_0   Omega_v  H_0'
c      WRITE(18,'(4X,3(F5.2,5X))') om,ov,h0
c     Ken's original
      WRITE(18,'(4X,3(F6.3,5X))') om,ov,h0
c -----------------------------------------------------------------------
c                          Building Hypercube! 
c -----------------------------------------------------------------------
      WRITE(*,'(1x,a)') 'Building hypercube ...'
      part_time=cputime()
      DO ired=1,nred
         av=avmin+(ired-1)*averr
         DO i=0,nz
            IF (cage.eq.'y'.or.cage.eq.'Y') ageuni=tz(h0,om,ov,z(i))
            DO itype=1,n_files
               IF (cage.eq.'y'.or.cage.eq.'Y') THEN
                  iagemax=2
                  IF (sed_abc(itype).eq.'BC') THEN
                     DO iage=2,KS(itype)
c                        IF (ageuni.gt.(T(iage-1)/1.e9).and.
c     .                       ageuni.le.(T(iage)/1.e9)) iagemax=iage
c                       Ken's original
                        IF (ageuni.ge.(T(iage-1)/1.e9).and.
     .                       ageuni.lt.(T(iage)/1.e9)) THEN 
                             iagemax=iage-1
                        END IF
                     END DO
                  END IF
               ELSE
                  iagemax=KS(itype)
               END IF
               DO iage=2,iagemax
                  DO k=1,iw(itype)
                     y(k)=f(k,iage,itype)
                     ww(k)=w(k,itype)
                  END DO
c                  IF (av.gt.0.0) THEN
c                Ken's original
                  IF (av.gt.-1.1) THEN
                     CALL red_spec(ww,y,av,yob,1,iw(itype),law_rd) 
                     DO k=1,iw(itype)
                        y(k)=yob(k)
                     END DO
                  END IF
c ..... Lyman alpha forest
                  CALL lyman_forest(ww,y,iw(itype),z(i))
c number of filters + the filter for absolute mag estimates
                  DO k=0,n_filter
                     flux(k,i,iage,itype,ired)=
     .                    filter(nfil(k),ww,y,iw(itype),z(i))/surf(k)
                  END DO
               END DO
            END DO
         END DO
      END DO

      WRITE(*,'(1x,a)') '                      ... done'
      part_time=cputime()-part_time
      WRITE(*,*) 'CPU time involved : ',part_time,' seconds'
c -----------------------------------------------------------------------
c                          Loop over the objects 
c -----------------------------------------------------------------------
      WRITE(*,'(1x,a)') 'I''m working ...'

      IF (ccat.eq.'y'.or.ccat.eq.'Y') 
     .     OPEN(10,FILE=spectra1,STATUS='unknown',ERR=81)
      IF (clog.eq.'y'.or.clog.eq.'Y')
     .     OPEN(14,FILE=spectra2,STATUS='unknown',ERR=81)
      IF (csed1.eq.'y'.or.csed1.eq.'Y')
     .     OPEN(15,FILE=spectra3,STATUS='unknown',ERR=81)
      IF (cz.eq.'y'.or.cz.eq.'Y')
     .     OPEN(17,FILE=spectra4,STATUS='unknown',ERR=81)
      IF (csed2.eq.'y'.or.csed2.eq.'Y')
     .     OPEN(19,FILE=spectra5,STATUS='unknown',ERR=81)

      OPEN(16,FILE=cat_in,STATUS='old',ERR=79)
c ..... Skipping the header (thanks to V. Le Brun)
      IF (header) THEN
         DO i=1,n_head
            READ(16,'(A)') header_line
         END DO
      END IF

      iob=0
13400 CONTINUE
      iob=iob+1
      IF (icat.eq.0) THEN
         READ(16,*,ERR=13500,END=13500) id,(amag(i),i=1,n_filter),
     .        (err_p(i),i=1,n_filter)
      ELSE IF (icat.eq.1) THEN
         READ(16,*,ERR=13500,END=13500) id,z1,z2,(amag(i),i=1,n_filter),
     .        (err_p(i),i=1,n_filter)
      END IF
     
      DO iii=1,n_filter
         IF (amag(iii).lt.90.) GO TO 13401
      END DO
      WRITE(*,*) 'WARNING: non detected object id ', id
      GO TO 13400
         
c ..... Computing the fluxes (unit 2.E-17) and errors in the filters
13401 iw_ob=0
      DO iii=1,n_filter
         i_lyman(iii)=0
         IF (err_p(iii).lt.err_min) err_p(iii)=err_min
         IF (amag(iii).gt.1..and.amag(iii).lt.90.) THEN
            iw_ob=iw_ob+1
            IF (magtyp.eq.0) fl(iii)=vega(iii)*
     .           10**(-0.4*(amag(iii)+c_corr(iii)))/surf(iii) 
            IF (magtyp.eq.1) fl(iii)=10**(-0.4*(amag(iii)+48.6))*
     .           (2.9979E+18/(wl(iii)**2))/2.E-17 
            err_fl(iii)=(2.3/2.5)*fl(iii)*err_p(iii)
c ..... Dereddening
            IF (ebv_mw.ne.0.) THEN
               CALL red_spec(wl(iii),fl(iii),av_mw,fl_ext,2,1,1)
               fl(iii)=fl_ext
            END IF
c amag -> Vega magnitude
            amag(iii)=-2.5*LOG10(fl(iii)*surf(iii)/vega(iii))
         ELSE IF (amag(iii).le.1.) THEN ! out of file objects
            fl(iii)=-1.0
         ELSE IF (amag(iii).ge.90.) THEN ! non detected objects
            IF (law_nd(iii).eq.0) THEN
               fl(iii)=-1.0
               err_fl(iii)=0.
            ELSE IF (law_nd(iii).eq.1) THEN
               fl(iii)=0.
               err_fl(iii)=fl_nd(iii)
               i_lyman(iii)=2
               iw_ob=iw_ob+1
            ELSE IF (law_nd(iii).eq.2) THEN
               fl(iii)=fl_nd(iii)/2.
               err_fl(iii)=fl(iii)
               i_lyman(iii)=1
               iw_ob=iw_ob+1
            ELSE IF (law_nd(iii).eq.3) THEN
               fl(iii)=fl_nd(iii)
               err_p(iii)=err_nd(iii)
               err_fl(iii)=(2.3/2.5)*fl(iii)*err_p(iii)
               i_lyman(iii)=1
               iw_ob=iw_ob+1
            END IF
         END IF
      END DO

      ndeg=iw_ob-1  ! degrees of freedom
      deg=REAL(ndeg)

      IF (ndeg.lt.1) THEN
         WRITE(*,*) 'WARNING: impossible to compute object ', iob
         GO TO 13400
      END IF
      
c -----------------------------------------------------------------------
c                         Computation 
c -----------------------------------------------------------------------

c ..... Open file for the matrix of results (if case)

      matrix=0
      WRITE(root,'(i8)') id
      DO k=1,8
         IF (root(k:k).ne.' ') THEN
            root=root(k:)
            GO TO 22
         END IF
      END DO
 22   iroot=INDEX(root,' ')-1
      IF (iout_m.gt.1) THEN
         amatrix=root(1:iroot)//'.m'
         matrix=1
         OPEN(11,FILE=amatrix,STATUS='unknown')
      END IF

      IF (icat.eq.1) THEN
c         DO i=1,nz-1
         DO i=0,nz-1
            IF (z1.ge.z(i).and.z1.lt.z(i+1)) i1=i
            IF (z2.gt.z(i).and.z2.le.z(i+1)) i2=i+1
         END DO
         IF (z2.gt.z(nz)) i2=nz
      ELSE IF (icat.eq.0) THEN
c         i1=1
         i1=0
         i2=nz
      END IF

      IF (z_clus.ge.0.) THEN
c         DO i=1,nz-1
         DO i=0,nz-1
            IF (z_clus.ge.z(i).and.z_clus.lt.z(i+1)) iz_clus=i
         END DO
         IF (z_clus.eq.z(nz)) iz_clus=nz
      END IF

c -----------------------------------------------------------------------
      DO i=i1,i2                                 !Scan in z
c -----------------------------------------------------------------------
         z_min_array(i)=1.E+30
         i_typ(i)=1
         i_age(i)=2
         i_red(i)=1
         IF (cage.eq.'y'.or.cage.eq.'Y') ageuni=tz(h0,om,ov,z(i))
         dif_w_min=1.e6
         n_fil_wa=1
         DO k=1,n_filter
            dif_w=ABS(wl(k)-wl(0)*(z_clus+1.))
            IF (amag(k).gt.1..and.amag(k).lt.90.) THEN
               IF (dif_w.lt.dif_w_min) THEN 
                  dif_w_min=dif_w
                  n_fil_wa=k
               END IF
            END IF
         END DO
c -----------------------------------------------------------------------
         DO itype=1,n_files     !Scan types
            IF (cage.eq.'y'.or.cage.eq.'Y') THEN
               iagemax=2
               IF (sed_abc(itype).eq.'BC') THEN
                  DO iage=2,KS(itype)
c                     IF (ageuni.gt.(T(iage-1)/1.e9).and.
c     .                   ageuni.le.(T(iage)/1.e9)) iagemax=iage
c                    Ken's original
                     IF (ageuni.ge.(T(iage-1)/1.e9).and.
     .                   ageuni.lt.(T(iage)/1.e9)) iagemax=iage-1
                  END DO
               END IF
            ELSE
               iagemax=KS(itype)
            END IF
c -----------------------------------------------------------------------
            DO iage=2,iagemax !Scan ages
c -----------------------------------------------------------------------
               DO ired=1,nred   !Scan reddening
c -----------------------------------------------------------------------
                  av=avmin+(ired-1)*averr
                  sxx=0.D0
                  sxy=0.D0
                  b=0.D0

c ..... Probability for z(i), iage, itype, ired

                  DO k=1,n_filter
                     IF (fl(k).ge.0.) THEN ! no oof objects
                        IF (flux(k,i,iage,itype,ired).eq.0..and.
     .                       i_lyman(k).eq.1) THEN
                           fl1(k)=0.
                        ELSE
                           fl1(k)=fl(k)
                           wt=1.D0/(err_fl(k)**2)
                           sxx=sxx+(flux(k,i,iage,itype,ired)**2)*wt
                           sxy=sxy+flux(k,i,iage,itype,ired)*fl1(k)*wt
                        END IF
                     ELSE
                        fl1(k)=-1.
                     END IF
                  END DO
                  b=sxy/sxx
                  chi2=0.
                  DO k=1,n_filter
                     IF (fl1(k).ge.0.) THEN ! no oof obj.
                        chi2=chi2+((fl1(k)-b*flux(k,i,iage,itype,ired))
     .                       /err_fl(k))**2
                     END IF
                  END DO
                  IF (chi2.lt.z_min_array(i)) THEN
                     z_min_array(i)=chi2
                     bb_array(i)=b
                     i_typ(i)=itype
                     i_age(i)=iage
                     i_red(i)=ired
                     IF (z_clus.ge.0..and.i.eq.iz_clus) THEN
                        itypea=itype
                        iagea=iage
                        ireda=ired
                     END IF
                  END IF

                  IF (sed_abc(itype).eq.'BC') THEN
                     age=T(iage)/1.e9
                  ELSE IF (sed_abc(itype).eq.'AS'.or.
     .                    sed_abc(itype).eq.'A ') THEN
                     age=19.9999
                  END IF

                  IF (matrix.gt.0) THEN
                     r_chi2=chi2/deg ! reduced chi2
                     WRITE(11,2012) id,z(i),itype,iage,age,av,r_chi2
                  END IF
               END DO
            END DO
         END DO

         q_pro(i)=100.*qchi2(0.5*deg,0.5*z_min_array(i))
         
         IF (sed_abc(i_typ(i)).eq.'BC') THEN
            age=T(i_age(i))/1.e9
         ELSE IF (sed_abc(i_typ(i)).eq.'AS'.or.
     .           sed_abc(i_typ(i)).eq.'A ') THEN
            age=19.9999
         END IF
         
         r_chi2=z_min_array(i)/deg ! reduced chi2

         IF(clog.eq.'y'.or.clog.eq.'Y') THEN
            WRITE(14,2002) id,z(i),r_chi2,q_pro(i),i_typ(i),age,
     .                     avmin+(i_red(i)-1)*averr
         END IF

      END DO
      CLOSE(11)                 !close output matrix

c -----------------------------------------------------------------------
c                     End scan
c -----------------------------------------------------------------------
c ..... Searching minima of chi2

      n_sol=0
      DO i=1,nz
         dif_w_min=1.e6
         n_fil_w(i)=1
         DO k=1,n_filter
            dif_w=ABS(wl(k)-wl(0)*(z(i)+1.))
            IF (amag(k).gt.1..and.amag(k).lt.90.
     .           .and.i_lyman(k).lt.1) THEN
               IF (dif_w.lt.dif_w_min) THEN   
                  dif_w_min=dif_w 
                  n_fil_w(i)=k  
               END IF
            END IF
         END DO
         dd=1.
         IF (i.eq.1) THEN
            dd=z_min_array(1)-z_min_array(2)
            d1=z_min_array(1)-z_min_array(2)
         ELSE IF (i.eq.nz) THEN
            dd=z_min_array(nz)-z_min_array(nz-1)
            d1=z_min_array(nz)-z_min_array(nz-1)
         ELSE
            d1=z_min_array(i)-z_min_array(i-1)
            d2=z_min_array(i+1)-z_min_array(i)
            IF (d1.ne.0.) dd=d2/d1
         END IF
         IF (dd.le.0..and.d1.lt.0.) THEN !dd < 0 extrema
            n_sol=n_sol+1
            izphot(n_sol)=i
         END IF
      END DO

c ..... Refinement around minima
      DO i=1,n_sol
         izsol=izphot(i)
         IF (izsol.eq.1) THEN
            dz_f=(z(izsol+1)-z(izsol))/10.
            zmax_f=z(izsol+1)
            zmin_f=z(izsol)
         ELSE IF (izsol.eq.nz) THEN
            dz_f=(z(izsol)-z(izsol-1))/10.
            zmax_f=z(izsol)
            zmin_f=z(izsol-1)
         ELSE
            dz_f=(z(izsol+1)-z(izsol-1))/20.
            zmax_f=z(izsol+1)
            zmin_f=z(izsol-1)
         END IF
c         IF (zmin_f.lt.z(1)) zmin_f=z(1)
         IF (zmin_f.lt.z(0)) zmin_f=z(0)
         IF (zmax_f.gt.z(nz)) zmax_f=z(nz)
         av=avmin+(i_red(izsol)-1)*averr
         CALL zfit_min(i_age(izsol),i_typ(izsol),av,
     .        zmin_f,zmax_f,dz_f,chi2,zop)
         IF (izsol.eq.1.or.izsol.eq.nz) THEN 
            IF (chi2.lt.z_min_array(izsol)) THEN
               zphot(i)=zop
               chi2sol(i)=chi2
               bb(i)=bm
            ELSE 
               zphot(i)=z(izsol)
               chi2sol(i)=z_min_array(izsol)
               bb(i)=bb_array(izsol)
            END IF
         ELSE
            zphot(i)=zop
            chi2sol(i)=chi2
            bb(i)=bm
         END IF
      END DO

c ..... Ordering solutions

      DO i=2,n_sol
         IF (chi2sol(i).lt.chi2sol(i-1)) THEN
            DO j=i-1,1,-1
               IF (chi2sol(j+1).lt.chi2sol(j)) THEN
                  dummy=chi2sol(j+1)
                  chi2sol(j+1)=chi2sol(j)
                  chi2sol(j)=dummy
                  dummy=zphot(j+1)
                  zphot(j+1)=zphot(j)
                  zphot(j)=dummy
                  idum=izphot(j+1)
                  izphot(j+1)=izphot(j)
                  izphot(j)=idum
                  ddummy=bb(j+1)
                  bb(j+1)=bb(j)
                  bb(j)=ddummy
               END IF
            END DO
         END IF
      END DO

c ..... Check absolute magnitude of the first solution
c       if amagb not in [amabs_1,amabs_2] search the next 
c       solution satisfying constraint
      DO i=1,nz
         iforbid(i)=0
      END DO
      iz=izphot(1)
c absolute magnitude
      amagb=0.0
      bmag=1000.0
      sauce=flux(n_fil_w(iz),iz,i_age(iz),i_typ(iz),i_red(iz))/
     .     flux(n_fil_w(iz),0,i_age(iz),i_typ(iz),i_red(iz))
      IF (sauce.gt.0.0) bmag=amag(n_fil_w(iz))+ ! absolute k magnitude
     .     2.5*LOG10(sauce)-d_mod(h0,om,ov,z(iz)) 
      sauce=flux(0,0,i_age(iz),i_typ(iz),i_red(iz))*surf(0)
      IF (sauce.gt.0.0) amagb=bmag-(zeropoint(n_fil_w(iz))-
     .     2.5*LOG10(flux(n_fil_w(iz),0,i_age(iz),
     .     i_typ(iz),i_red(iz))*surf(n_fil_w(iz))/sauce)) 
c constraint absolute magnitude
      IF (amagb.ge.amabs_1.and.amagb.le.amabs_2) GO TO 2000 ! solution ok
c first solution forbidden, search the next one
      iforbid(iz)=1
 200  chi2min=1.e6
      DO i=1,nz
         IF (z_min_array(i).lt.chi2min.and.iforbid(i).ne.1) THEN
            chi2min=z_min_array(i)
            iz=i
         END IF
      END DO
      DO i=1,n_sol
         IF (iz.eq.izphot(i)) THEN 
c check the absolute magnitude for the i-th solution
            amagb=0.0
            bmag=1000.0
            sauce=flux(n_fil_w(iz),iz,i_age(iz),i_typ(iz),i_red(iz))/
     .           flux(n_fil_w(iz),0,i_age(iz),i_typ(iz),i_red(iz))
            IF (sauce.gt.0.0) bmag=amag(n_fil_w(iz))+ ! absolute k magnitude
     .           2.5*LOG10(sauce)-d_mod(h0,om,ov,z(iz)) 
            sauce=flux(0,0,i_age(iz),i_typ(iz),i_red(iz))*surf(0)
            IF (sauce.gt.0.0) amagb=bmag-(zeropoint(n_fil_w(iz))-
     .           2.5*LOG10(flux(n_fil_w(iz),0,i_age(iz),
     .           i_typ(iz),i_red(iz))*surf(n_fil_w(iz))/sauce)) 
c constraint absolute magnitude
c i-th solution ok -> re-order solutions
            IF (amagb.ge.amabs_1.and.amagb.le.amabs_2) GO TO 3000 
c i-th solution forbidden
            iforbid(iz)=1
            GO TO 200 
         END IF
      END DO
c refinement if the new minimum is not a secondary solution
      IF (iz.eq.1) THEN
         dz_f=(z(iz+1)-z(iz))/10.
         zmax_f=z(iz+1)
         zmin_f=z(iz)
      ELSE IF (iz.eq.nz) THEN
         dz_f=(z(iz)-z(iz-1))/10.
         zmax_f=z(iz)
         zmin_f=z(iz-1)
      ELSE
         dz_f=(z(iz+1)-z(iz-1))/20.
         zmax_f=z(iz+1)
         zmin_f=z(iz-1)
      END IF
c      IF (zmin_f.lt.z(1)) zmin_f=z(1)
      IF (zmin_f.lt.z(0)) zmin_f=z(0)
      IF (zmax_f.gt.z(nz)) zmax_f=z(nz)
      av=avmin+(i_red(iz)-1)*averr
      CALL zfit_min(i_age(iz),i_typ(iz),av,
     .     zmin_f,zmax_f,dz_f,chi2,zop)
c check absolut magnitude
      amagb=0.0
      bmag=1000.0
      sauce=flux(n_fil_w(iz),iz,i_age(iz),i_typ(iz),i_red(iz))/
     .     flux(n_fil_w(iz),0,i_age(iz),i_typ(iz),i_red(iz))
      IF (sauce.gt.0.0) bmag=amag(n_fil_w(iz))+ 
     .     2.5*LOG10(sauce)-d_mod(h0,om,ov,z(iz)) 
      sauce=flux(0,0,i_age(iz),i_typ(iz),i_red(iz))*surf(0)
      IF (sauce.gt.0.0) amagb=bmag-(zeropoint(n_fil_w(iz))-
     .     2.5*LOG10(flux(n_fil_w(iz),0,i_age(iz),
     .     i_typ(iz),i_red(iz))*surf(n_fil_w(iz))/sauce)) 
      IF (amagb.ge.amabs_1.and.amagb.le.amabs_2) THEN
c found new first solution, add solution and shift the others
         n_sol=n_sol+1
         DO i=n_sol,2,-1
            chi2sol(i)=chi2sol(i-1)
            zphot(i)=zphot(i-1)
            izphot(i)=izphot(i-1)
            bb(i)=bb(i-1)
         END DO
         zphot(1)=zop
         chi2sol(1)=chi2
         bb(1)=bm
         izphot(1)=iz
         GO TO 2000
      ELSE
c forbidden solution, try again
         iforbid(iz)=1
         GO TO 200
      END IF

c re-ordering
 3000 izmin=izphot(i)
      dummy=chi2sol(i)
      dummyz=zphot(i)
      idum=izphot(i)
      ddummy=bb(i)
      i_sol=i
      DO i=i_sol,2,-1
         chi2sol(i)=chi2sol(i-1)
         zphot(i)=zphot(i-1)
         izphot(i)=izphot(i-1)
         bb(i)=bb(i-1)
      END DO
      chi2sol(1)=dummy
      zphot(1)=dummyz
      izphot(1)=idum
      bb(1)=ddummy

c ..... main solution 
 2000 i_sol=1
      izmin=izphot(1)
      itypem=i_typ(izmin)
      iagem=i_age(izmin)
      iredm=i_red(izmin)
      av=avmin+(iredm-1)*averr
      DO k=1,n_filter
         IF (flux(k,izmin,iagem,itypem,iredm).eq.0..and.
     .        i_lyman(k).eq.1) THEN
            fl1(k)=0.
         ELSE 
            fl1(k)=fl(k)
         END IF
      END DO
      IF (sed_abc(i_typ(izmin)).eq.'BC') THEN
         age=T(i_age(izmin))/1.e9
      ELSE IF (sed_abc(i_typ(izmin)).eq.'AS'.or.
     .        sed_abc(i_typ(izmin)).eq.'A ') THEN
         age=19.9999
      END IF
      r_chi2=chi2sol(1)/deg
      IF (r_chi2.ge.1.e4) r_chi2=9999.99
      qq=100.*qchi2(0.5*deg,0.5*chi2sol(1))      
      IF (z(izmin)-zphot(1).gt.dz_f/2.) THEN
         izmin1=izmin-1
      ELSE
         izmin1=izmin
      END IF

c ..... study of confidence intervals
      DO i=1,3
         zz1(i)=z1
         zz2(i)=z2
      END DO
      DO i=1,mxz
         q_pro1(i)=0.
         q_pro2(i)=0.
         q_proz(i)=0.
      END DO
      poids=0.0
      poids_z=0.0
      z_poids=0.
      qq_p=0.

c ..... Values after izmin
      z_b1(1)=zphot(1)
      chi2_b1(1)=chi2sol(1)
      q_pro1(1)=qq
      i=izmin1+1
      npoiz1=1
      IF (ABS(z(izmin1)-zphot(1)).lt.1.e-5.and.izmin1.eq.nz) THEN
         DO ii=1,3              ! z at P = 99 90 and 68%
            zz2(ii)=z2
         END DO
         npoiz1=1
      ELSE
 1001    CONTINUE
         k=i-izmin1+1
         q_pro1(k)=q_pro(i)
         z_b1(k)=z(i)
         chi2_b1(k)=z_min_array(i)

         IF (i.lt.nz.and.i.le.izmin1) THEN
            i=i+1
            GO TO 1001
         ELSE IF (i.lt.nz.and.(chi2_b1(k)-chi2sol(1)).le.val(1)) THEN
            i=i+1
            GO TO 1001
         END IF
         npoiz1=k

         DO ii=1,3              ! z at P = 99 90 and 68%
            l=0
            zz2(ii)=rlin(val(ii)+chi2sol(1),chi2_b1,z_b1,npoiz1,l)
            IF (zz2(ii).lt.z1) zz2(ii)=z1
            IF (zz2(ii).gt.z2) zz2(ii)=z2
            IF (zz2(ii).lt.zphot(1)) zz2(ii)=z2
         END DO
      END IF

c ..... Values before izmin

      z_b2(1)=zphot(1)
      chi2_b2(1)=chi2sol(1)
      i=izmin1
      ipass=0
      npoiz2=1
      IF (ABS(z(izmin1)-zphot(1)).lt.1.e-5) THEN
         q_pro2(1)=q_pro(izmin1)
         ipass=1
      END IF
      IF (ABS(z(izmin1)-zphot(1)).lt.1.e-5.and.izmin1.eq.1) THEN
         DO ii=1,3                 !z at P = 99 90 and 68%
            zz1(ii)=z1
         END DO
         npoiz2=1
      ELSE
 1003    CONTINUE
         k=izmin1-i+2
         q_pro2(k)=q_pro(i-ipass)
         z_b2(k)=z(i-ipass)
         chi2_b2(k)=z_min_array(i-ipass)
         IF (i-ipass.gt.1.and.i-ipass.ge.izmin1) THEN
            i=i-1
            GO TO 1003
         ELSE IF (i-ipass.gt.1.and.
     .           (chi2_b2(k)-chi2sol(1)).le.val(1)) THEN
            i=i-1
            GO TO 1003
         END IF
         npoiz2=k

         DO ii=1,3                 !z at P = 99 90 and 68%
         l=0
            zz1(ii)=rlin(val(ii)+chi2sol(1),chi2_b2,z_b2,npoiz2,l)
            IF (zz1(ii).lt.z1) zz1(ii)=z1
            IF (zz1(ii).gt.z2) zz1(ii)=z2
            IF (zz1(ii).gt.zphot(1)) zz1(ii)=zphot(1)
         END DO
      END IF

c ..... Ordering q_proz et z_b

      IF (ABS(z(izmin1)-zphot(1)).lt.1.e-5.and.izmin1.eq.1) THEN
         npoiz=npoiz1           ! points after only 
         DO ii=1,npoiz
            q_proz(ii)=q_pro1(ii)
            z_b(ii)=z_b1(ii)
         END DO         
      ELSE IF (ABS(z(izmin1)-zphot(1)).lt.1.e-5.and.izmin1.eq.nz) THEN
         npoiz=npoiz2           ! points before only  
         DO ii=1,npoiz
            q_proz(ii)=q_pro2(npoiz2-(ii-1))
            z_b(ii)=z_b2(npoiz2-(ii-1))
         END DO         
      ELSE IF (ABS(z(izmin1)-zphot(1)).lt.1.e-5.and.
     .        izmin1.ne.1.and.izmin1.ne.nz) THEN
         npoiz=npoiz1+npoiz2-1
         DO ii=1,npoiz2         ! points before
            q_proz(ii)=q_pro2(npoiz2-(ii-1))
            z_b(ii)=z_b2(npoiz2-(ii-1))
         END DO
         DO ii=1,npoiz1-1       ! points after (skip zphot(1))
            q_proz(npoiz2+ii)=q_pro1(ii+1)
            z_b(npoiz2+ii)=z_b1(ii+1)
         END DO
      ELSE
         npoiz=npoiz1+npoiz2-2
         DO ii=1,npoiz2-1       ! points before (skip zphot(1))
            q_proz(ii)=q_pro2(npoiz2-(ii-1))
            z_b(ii)=z_b2(npoiz2-(ii-1))
         END DO
         DO ii=1,npoiz1-1       ! points after (skip zphot(1))
            q_proz(npoiz2+ii-1)=q_pro1(ii+1)
            z_b(npoiz2+ii-1)=z_b1(ii+1)
         END DO
      END IF

      DO ii=1,npoiz
         poids=poids+q_proz(ii)
         poids_z=poids_z+q_proz(ii)*z_b(ii)
      END DO

c ..... z weighted mean

      IF (poids.gt.0.0) THEN
         z_poids=poids_z/poids
         l=0
         qq_p=rlin(z_poids,z_b,q_proz,npoiz,l)
      ELSE
         z_poids=zphot(1)
         qq_p=qq
      END IF

c ..... secondary solution
      IF (n_sol.lt.2) THEN      ! no secondary solution
         zphot(2)=9.99
         qq2=0.
      ELSE
         qq2=100.*qchi2(0.5*deg,0.5*chi2sol(2))
      END IF

c ..... exact absolute magnitude
c      amagb=0.0
c      bmag=1000.0
c      DO k=1,iw(itypem)
c         y(k)=f(k,iagem,itypem)
c         ww(k)=w(k,itypem)
c      END DO
c      IF (av.gt.0.0) THEN
c         CALL red_spec(ww,y,av,yob,1,iw(itypem),law_rd)
c         DO k=1,iw(itypem)
c            y(k)=yob(k)
c         END DO
c      END IF
c      CALL lyman_forest(ww,y,iw(itypem),zphot(1))
c      bmag=amag(n_fil_w(iz))-d_mod(h0,om,ov,zphot(1))+2.5*
c     .     LOG10(filter(nfil(n_fil_w(iz)),ww,y,iw(itypem),zphot(1))/
c     .     filter(nfil(n_fil_w(iz)),ww,y,iw(itypem),0.))
c      amagb=bmag-zeropoint(n_fil_w(iz))
c     .     +2.5*LOG10(filter(nfil(n_fil_w(iz)),ww,y,iw(itypem),0.)/
c     .     filter(nfilt_mabs,ww,y,iw(itypem),0.))

c ..... Probability to be in the cluster
      IF (z_clus.ge.0.) THEN
         IF (z_clus.eq.z(iz_clus)) THEN
            p_clus=q_pro(iz_clus)
         ELSE
            p_clus=q_pro(iz_clus)+(q_pro(iz_clus+1)-q_pro(iz_clus))
     .           *(z_clus-z(iz_clus))/(z(iz_clus+1)-z(iz_clus))
         END IF

c ..... Absolute magnitude in the cluster
         DO k=1,iw(itypea)
            y(k)=f(k,iagea,itypea)
            ww(k)=w(k,itypea)
         END DO
c         IF (avmin+(ireda-1)*averr.gt.0.0) THEN
c        Ken's original
         IF (avmin+(ireda-1)*averr.gt.-1.1) THEN
            CALL red_spec(ww,y,av,yob,1,iw(itypea),law_rd)
            DO k=1,iw(itypea)
               y(k)=yob(k)
            END DO
         END IF
         CALL lyman_forest(ww,y,iw(itypea),z_clus)
         bmag=amag(n_fil_wa)-d_mod(h0,om,ov,z_clus)+2.5*
     .        LOG10(filter(nfil(n_fil_wa),ww,y,iw(itypea),z_clus)/
     .        filter(nfil(n_fil_wa),ww,y,iw(itypea),0.))
         amagb_clus=bmag-zeropoint(n_fil_wa)
     .        +2.5*LOG10(filter(nfil(n_fil_wa),ww,y,iw(itypea),0.)/
     .        filter(nfilt_mabs,ww,y,iw(itypea),0.))

c ..... writing output files
         IF (cz.eq.'y'.or.cz.eq.'Y') WRITE(17,2009) id,zphot(1),r_chi2,
     .        qq,itypem,iagem,age,av,bb(1),zz1(1),zz2(1),zz1(2),
     .        zz2(2),zz1(3),zz2(3),z_poids,qq_p,amagb,zphot(2),qq2,
     .        amagb_clus,p_clus,amag(n_fil_wa)
      ELSE
         IF (cz.eq.'y'.or.cz.eq.'Y') WRITE(17,2008) id,zphot(1),r_chi2,
     .        qq,itypem,iagem,age,av,bb(1),zz1(1),zz2(1),zz1(2),
     .        zz2(2),zz1(3),zz2(3),z_poids,qq_p,amagb,zphot(2),qq2
      END IF

      IF (ccat.eq.'y'.or.ccat.eq.'Y') WRITE(10,2001) id,zphot(1),r_chi2,
     .     qq,itypem,iagem,age,av,bb(1),zz1(1),zz2(1),
     .     zz1(2),zz2(2),zz1(3),zz2(3),z_poids,qq_p,amagb,i_sol

c -------------------------------------------------------------------------

c ..... other solutions (if .cat_phot requested, check threshold)

      IF (ccat.eq.'y'.or.ccat.eq.'Y') THEN
         DO i_sol=2,n_sol
            qq=100.*qchi2(0.5*deg,0.5*chi2sol(i_sol))
            IF (qq.ge.pthresh) THEN
               izmin=izphot(i_sol)
               IF (sed_abc(i_typ(izmin)).eq.'BC') THEN
                  age=T(i_age(izmin))/1.e9
               ELSE IF (sed_abc(i_typ(izmin)).eq.'AS'.or.
     .                 sed_abc(i_typ(izmin)).eq.'A ') THEN
                  age=19.9999
               END IF
               r_chi2=chi2sol(i_sol)/deg
               IF (r_chi2.ge.1.e4) r_chi2=9999.99
               qq=100.*qchi2(0.5*deg,0.5*chi2sol(i_sol))      
               IF (z(izmin)-zphot(i_sol).gt.dz_f/2.) THEN
                  izmin1=izmin-1
               ELSE
                  izmin1=izmin
               END IF

c ..... study of confidence intervals
               DO i=1,3
                  zz1(i)=z1
                  zz2(i)=z2
               END DO
               DO i=1,mxz
                  q_pro1(i)=0.
                  q_pro2(i)=0.
                  q_proz(i)=0.
               END DO
               poids=0.0
               poids_z=0.0
               z_poids=0.
               qq_p=0.

c ..... Values after izmin
               z_b1(1)=zphot(i_sol)
               chi2_b1(1)=chi2sol(i_sol)
               q_pro1(1)=qq
               i=izmin1+1
               npoiz1=1
               IF (ABS(z(izmin1)-zphot(i_sol)).lt.1.e-5.and.
     .              izmin1.eq.nz) THEN
                  DO ii=1,3     ! z at P = 99 90 and 68%
                     zz2(ii)=z2
                  END DO
                  npoiz1=1
               ELSE
 1002             CONTINUE
                  k=i-izmin1+1
                  q_pro1(k)=q_pro(i)
                  z_b1(k)=z(i)
                  chi2_b1(k)=z_min_array(i)

                  IF (i.lt.nz.and.i.le.izmin1) THEN
                     i=i+1
                     GO TO 1002
                  ELSE IF (i.lt.nz.and.
     .                    (chi2_b1(k)-chi2sol(i_sol)).le.val(1)) THEN
                     i=i+1
                     GO TO 1002
                  END IF
                  npoiz1=k

                  DO ii=1,3     ! z at P = 99 90 and 68%
                     l=0
                     zz2(ii)=rlin(val(ii)+chi2sol(i_sol),
     .                    chi2_b1,z_b1,npoiz1,l)
                     IF (zz2(ii).lt.z1) zz2(ii)=z1
                     IF (zz2(ii).gt.z2) zz2(ii)=z2
                     IF (zz2(ii).lt.zphot(i_sol)) zz2(ii)=z2
                  END DO
               END IF

c ..... Values before izmin

               z_b2(1)=zphot(i_sol)
               chi2_b2(1)=chi2sol(i_sol)
               i=izmin1
               ipass=0
               npoiz2=1
               IF (ABS(z(izmin1)-zphot(i_sol)).lt.1.e-5) THEN
                  q_pro2(1)=q_pro(izmin1)
                  ipass=1
               END IF
               IF (ABS(z(izmin1)-zphot(i_sol)).lt.1.e-5.and.
     .              izmin1.eq.1) THEN
                  DO ii=1,3
                     zz1(ii)=z1
                  END DO
                  npoiz2=1
               ELSE
 1004             CONTINUE
                  k=izmin1-i+2
                  q_pro2(k)=q_pro(i-ipass)
                  z_b2(k)=z(i-ipass)
                  chi2_b2(k)=z_min_array(i-ipass)
                  IF (i-ipass.gt.1.and.i-ipass.ge.izmin1) THEN
                     i=i-1
                     GO TO 1004
                  ELSE IF (i-ipass.gt.1.and.
     .                    (chi2_b2(k)-chi2sol(i_sol)).le.val(1)) THEN
                     i=i-1
                     GO TO 1004
                  END IF
                  npoiz2=k

                  DO ii=1,3     !z at P = 99 90 and 68%
                     l=0
                     zz1(ii)=rlin(val(ii)+chi2sol(i_sol),
     .                    chi2_b2,z_b2,npoiz2,l)
                     IF (zz1(ii).lt.z1) zz1(ii)=z1
                     IF (zz1(ii).gt.z2) zz1(ii)=z2
                     IF (zz1(ii).gt.zphot(i_sol)) zz1(ii)=zphot(i_sol)
                  END DO
               END IF

c ..... Ordering q_proz et z_b

               IF (ABS(z(izmin1)-zphot(i_sol)).lt.1.e-5.
     .              and.izmin1.eq.1) THEN
                  npoiz=npoiz1  ! points after only 
                  DO ii=1,npoiz
                     q_proz(ii)=q_pro1(ii)
                     z_b(ii)=z_b1(ii)
                  END DO         
               ELSE IF (ABS(z(izmin1)-zphot(i_sol)).lt.1.e-5.and.
     .                 izmin1.eq.nz) THEN
                  npoiz=npoiz2  ! points before only  
                  DO ii=1,npoiz
                     q_proz(ii)=q_pro2(npoiz2-(ii-1))
                     z_b(ii)=z_b2(npoiz2-(ii-1))
                  END DO         
               ELSE IF (ABS(z(izmin1)-zphot(i_sol)).lt.1.e-5.and.
     .                 izmin1.ne.1.and.izmin1.ne.nz) THEN
                  npoiz=npoiz1+npoiz2-1
                  DO ii=1,npoiz2 ! points before
                     q_proz(ii)=q_pro2(npoiz2-(ii-1))
                     z_b(ii)=z_b2(npoiz2-(ii-1))
                  END DO
                  DO ii=1,npoiz1-1 ! points after (skip zphot(i_sol))
                     q_proz(npoiz2+ii)=q_pro1(ii+1)
                     z_b(npoiz2+ii)=z_b1(ii+1)
                  END DO
               ELSE
                  npoiz=npoiz1+npoiz2-2
                  DO ii=1,npoiz2-1 ! points before (skip zphot(i_sol))
                     q_proz(ii)=q_pro2(npoiz2-(ii-1))
                     z_b(ii)=z_b2(npoiz2-(ii-1))
                  END DO
                  DO ii=1,npoiz1-1 ! points after (skip zphot(i_sol))
                     q_proz(npoiz2+ii-1)=q_pro1(ii+1)
                     z_b(npoiz2+ii-1)=z_b1(ii+1)
                  END DO
               END IF

               DO ii=1,npoiz
                  poids=poids+q_proz(ii)
                  poids_z=poids_z+q_proz(ii)*z_b(ii)
               END DO

c ..... absolute magnitude
               itypem=i_typ(izmin)
               iagem=i_age(izmin)
               iredm=i_red(izmin)
               av=avmin+(iredm-1)*averr
               sauce=flux(n_fil_w(izmin),izmin,iagem,itypem,iredm)/
     .              flux(n_fil_w(izmin),0,iagem,itypem,iredm)
               IF (sauce.gt.0.0) bmag=amag(n_fil_w(izmin))+
     .              2.5*LOG10(sauce)-d_mod(h0,om,ov,zphot(i_sol)) 

               sauce=flux(0,0,iagem,itypem,iredm)*surf(0)
               IF (sauce.gt.0.0) amagb=bmag-(zeropoint(n_fil_w(izmin))-
     .              2.5*LOG10(flux(n_fil_w(izmin),0,iagem,itypem,iredm)
     .              *surf(n_fil_w(izmin))/sauce)) 

c               DO k=1,iw(itypem)
c                  y(k)=f(k,iagem,itypem)
c                  ww(k)=w(k,itypem)
c               END DO
c               IF (av.gt.0.0) THEN
c                  CALL red_spec(ww,y,av,yob,1,iw(itypem),law_rd)
c                  DO k=1,iw(itypem)
c                     y(k)=yob(k)
c                  END DO
c               END IF
c               CALL lyman_forest(ww,y,iw(itypem),zphot(i_sol))
c               bmag=amag(n_fil_w(izmin))-d_mod(h0,om,ov,zphot(i_sol))+
c     .              2.5*LOG10(filter(nfil(n_fil_w(izmin)),ww,y,
c     .              iw(itypem),zphot(i_sol))/
c     .              filter(nfil(n_fil_w(izmin)),ww,y,iw(itypem),0.))
c               amagb=bmag-zeropoint(n_fil_w(izmin))
c     .              +2.5*LOG10(filter(nfil(n_fil_w(izmin)),ww,y,
c     .              iw(itypem),0.)/
c     .              filter(nfilt_mabs,ww,y,iw(itypem),0.))

c ..... z weighted mean
               IF (poids.gt.0.0) THEN
                  z_poids=poids_z/poids
                  l=0
                  qq_p=rlin(z_poids,z_b,q_proz,npoiz,l)
               ELSE 
                  z_poids=zphot(i_sol)
                  qq_p=qq
               END IF

               WRITE(10,2001) id,zphot(i_sol),r_chi2,qq,itypem,iagem,
     .              age,av,bb(i_sol),zz1(1),zz2(1),zz1(2),
     .              zz2(2),zz1(3),zz2(3),z_poids,qq_p,amagb,i_sol
            END IF
         END DO
      END IF

c -----------------------------------------------------------------------
c                    Graphics: Synthetic spectrum at zphot
c -----------------------------------------------------------------------

c ..... Write SED of object 'id'

      izmin=izphot(1)
      itypem=i_typ(izmin)
      iagem=i_age(izmin)
      iredm=i_red(izmin)
      av=avmin+(iredm-1)*averr

      WRITE(fmt1,2006) n_filter,n_filter
      WRITE(fmt2,2007) n_filter,n_filter
 2006 FORMAT('(1x,i5,',i2,'(1x,F8.4),',i2,'(1x,F8.4))')
 2007 FORMAT('(1x,i5,',i2,'(1x,E9.3),',i2,'(1x,E9.3))')
      IF (iunit_sed.eq.1) THEN
         DO k=1,n_filter
            f_jy(k)=fl1(k)*2.E-17*(wl(k)**2/(2.9979E+18))*1.E+29
            err_jy(k)=(2.3/2.5)*f_jy(k)*err_p(k)
         END DO
      END IF
      IF (iunit_sed.eq.2) THEN
c ..... Vega mag -> AB mag
         DO k=1,n_filter
            amag(k)=amag(k)+conv_mag(k)
            IF (amag(k).gt.90.) amag(k)=mag_nd(k)+conv_mag(k)
         END DO
      END IF
      IF (csed1.eq.'y'.or.csed1.eq.'Y') THEN
         IF (iunit_sed.eq.0) THEN
            WRITE(15,fmt2) id,(fl1(i),i=1,n_filter),
     .           (err_fl(i),i=1,n_filter)
         ELSE IF (iunit_sed.eq.1) THEN
            WRITE(15,fmt2) id,(f_jy(i),i=1,n_filter),
     .           (err_jy(i),i=1,n_filter)
         ELSE IF (iunit_sed.eq.2) THEN
            WRITE(15,fmt1) id,(amag(i),i=1,n_filter),
     .           (err_p(i),i=1,n_filter)
         END IF
      END IF

      WRITE(fmt3,2010) n_filter
      WRITE(fmt4,2011) n_filter
 2010 FORMAT('(1x,i5,',i2,'(1x,F8.4))')
 2011 FORMAT('(1x,i5,',i2,'(1x,E9.3))')
      IF (csed2.eq.'y'.or.csed2.eq.'Y'.or.iout_s.gt.1) THEN 
         DO k=1,iw(itypem)
            y(k)=f(k,iagem,itypem)
            ww(k)=w(k,itypem)
         END DO
         av=avmin+(iredm-1)*averr
c         IF (av.gt.0.0) THEN
c        Ken's original
         IF (av.gt.-1.1) THEN
            CALL red_spec(ww,y,av,yob,1,iw(itypem),law_rd)
            DO k=1,iw(itypem)  
               y(k)=yob(k)
            END DO
         END IF
         CALL lyman_forest(ww,y,iw(itypem),zphot(1))
      END IF

      IF (csed2.eq.'y'.or.csed2.eq.'Y') THEN 
         DO k=1,n_filter
            fl_temp(k)=filter(nfil(k),ww,y,iw(itypem),zphot(1))/surf(k)
         END DO
         IF (iunit_sed.eq.0) THEN
            WRITE(19,fmt4) id,(bb(1)*fl_temp(k),k=1,n_filter)
         ELSE IF (iunit_sed.eq.1) THEN
            WRITE(19,fmt4) id,(bb(1)*2.E-17*fl_temp(k)*
     .     (wl(k)**2/(2.9979E+18))*1.E+29,k=1,n_filter)
         ELSE IF (iunit_sed.eq.2) THEN
            DO k=1,n_filter
               IF (fl_temp(k).gt.0.) THEN
                  mag_temp(k)=-2.5*LOG10(bb(1)*fl_temp(k)*
     .                 surf(k)/vega(k))+conv_mag(k)
               ELSE
                  mag_temp(k)=100.
               END IF
            END DO
            WRITE(19,fmt3) id,(mag_temp(k),k=1,n_filter)
         END IF
      END IF

      IF (iout_s.gt.1) THEN
         spectra=root(1:iroot)//'.spe'

         IF (iunit_sed.eq.0) THEN
            DO k=1,iw(itypem)
               y(k)=bb(1)*y(k)/(1.+zphot(1))
            END DO
         ELSE IF (iunit_sed.eq.1) THEN
            DO k=1,iw(itypem)
               y(k)=bb(1)*y(k)*2.E-17*(ww(k)**2/2.9979E+18)*1.E+29*
     .              (1.+zphot(1))
            END DO
         ELSE IF (iunit_sed.eq.2) THEN
            DO k=1,iw(itypem)
               y(k)=bb(1)*y(k)*2.E-17*(ww(k)**2/2.9979E+18)*
     .              (1.+zphot(1))
               IF (y(k).ne.0.) THEN
                  y(k)=-2.5*LOG10(y(k))-48.60
               ELSE
                  y(k)=100.
               END IF
            END DO
         END IF
         iwav=0
         DO k=1,iw(itypem)
            w0=ww(k)*(1.+zphot(1))
            wmin1=wl(1)-surf(1)/2.-1000.
            wmax1=wl(n_filter)+surf(n_filter)/2.+1000.
            IF (w0.ge.wmin1.and.w0.le.wmax1) THEN
               iwav=iwav+1
               wsyn(iwav)=ww(k)*(1.+zphot(1))
               ysyn(iwav)=y(k)
            END IF
         END DO
         OPEN(12,FILE=spectra,STATUS='unknown')
         DO k=1,iwav
            WRITE(12,*) wsyn(k),ysyn(k)
         END DO
         CLOSE(12)
      END IF

c.......................................................................
      GO TO 13400
13500 CONTINUE
c.......................................................................

      CLOSE(16)
      IF (ccat.eq.'y'.or.ccat.eq.'Y') CLOSE(10)
      IF (clog.eq.'y'.or.clog.eq.'Y') CLOSE(14)
      IF (csed1.eq.'y'.or.csed1.eq.'Y') CLOSE(15)
      IF (cz.eq.'y'.or.cz.eq.'Y') CLOSE(17)
      IF (csed2.eq.'y'.or.csed2.eq.'Y') CLOSE(19)

      time_tot=cputime()-time_tot
      WRITE(*,'(a)') ' Done !!!'
      WRITE(*,*) 'Total CPU time involved :', time_tot, ' seconds'
      WRITE(18,'(a,1x,f8.2,1x,a)') '#  Total CPU time involved :',  
     .     time_tot,' seconds'
      CLOSE(18)
      STOP
 2001 FORMAT(i10,1x,f6.3,1x,f8.3,1x,f6.2,1x,i2,1x,i3,1x,f8.4,1x,
     .     f5.2,1x,e9.3,1x,6(f6.3,1x),f5.3,1x,f6.2,1x,f7.2,1x,i3)
 2002 FORMAT(i10,1x,f5.2,1x,e9.3,1x,f6.2,1x,i2,1x,f8.4,1x,F4.2)
 2008 FORMAT(i10,1x,f6.3,1x,f8.3,1x,f6.2,1x,i2,1x,i3,1x,f8.4,1x,
     .     f5.2,1x,e9.3,1x,6(f6.3,1x),f5.3,1x,f6.2,1x,f7.2,1x,
     .     f6.3,1x,f6.2)
 2009 FORMAT(i10,1x,f6.3,1x,f8.3,1x,f6.2,1x,i2,1x,i3,1x,f8.4,1x,
     .     f5.2,1x,e9.3,1x,6(f6.3,1x),f5.3,1x,f6.2,1x,f7.2,1x,
     .     f6.3,1x,f6.2,1x,f6.2,1x,f6.2,1x,f6.2)
 2012 FORMAT(i10,1x,f5.2,1x,i2,1x,i3,1x,f8.4,1x,f5.2,1x,E9.3)
 25   STOP 'Error opening synthetic spectra files'
 78   STOP 'Error opening templates file'
 79   STOP 'Error opening or reading photometric catalog'
 80   STOP 'Error opening filters file'
 81   STOP 'Error opening output files'
 82   STOP 'Error reading filters file'
      END
c----------------------------------------------------------------------
      FUNCTION func(x)
      REAL func,pchi2,x,deg,prob
      EXTERNAL pchi2
      COMMON/nfree/deg,prob
c
      func=prob-pchi2(0.5*deg,0.5*x)
      RETURN
      END
c----------------------------------------------------------------------
c ..... Distance module: d_mod = 5*LOG10(d_L[pc]/10)
 
      FUNCTION d_mod(h0,om,ov,za)
      IMPLICIT NONE
      INTEGER nstep,i
      REAL d_mod,h0,om,ov,z,za,c,dmin,d_l,ee,s,s1,s2,dz,a,b,x
      PARAMETER (c=2.9979E+5,nstep=100)
      ee(x)=SQRT(om*(1.+x)**3+(1.-om-ov)*(1.+x)**2+ov)
c
      IF (za.lt.1.e-5) THEN
         z=1.e-5
      ELSE
         z=za
      END IF
      IF (om.lt.0.001.and.ov.eq.0.) THEN
         d_l=(1.+z)*LOG(1.+z)
      ELSE IF (om.eq.1..and.ov.eq.0.) THEN
         d_l=2.*(1.+z)*(1.-1./SQRT(1.+z))
      ELSE IF (om.lt.1..and.om.gt.0.01.and.ov.eq.0.) THEN
         d_l=(2./om**2)*(om*z-(2.-om)*(-1.+SQRT(om*z+1.))) ! Mattig
      ELSE
         dz=z/FLOAT(nstep)
         s=0.
         d_l=0.
         DO i=1,nstep+1
            a=(i-1)*dz
            b=i*dz
            s1=(1/ee(a))
            s2=(1/ee(b))
            s=dz*(s1+s2)/2.
            d_l=d_l+s
         END DO
         d_l=(1.+z)*d_l
      END IF
      d_mod=5.*LOG10((c/h0)*d_l)+25.
c   minimum d_mod: corresponding to d(~M31)=1000 Kpc
      dmin=25.
      IF (d_mod.lt.dmin) d_mod=dmin
      RETURN
      END
c ----------------------------------------------------------------------




