      PROGRAM CALIB  
C
C     This fortran code applies the calibrations proposed in the paper 
C     "A calibration of Geneva photometry for B to G stars in terms of Teff,
C     logg and [M/H]" (A&AS,1996, , ).
C
C     It allows to determine the physical parameters Teff and logg for 
C     unreddned B to mid-G stars of the main sequence or just above it.
C     Metallicity is also given for cool stars ( T<7500K ), but for hotter
C     stars, [M/H] has to be known a priori and can take value comprised
C     betweem -1 and 1.
C
C
C
C     We have defined three range of temperature:
C
C     - For hot stars (Teff>11000 K), we use the reddening-free parameters X
C       and Y defined by Cramer and Maeder (1979) by 
C
C            X= 0.3788+1.3784U-1.2162B1-0.8498B2-0.1554V1+0.8450G
C            Y=-0.8288+0.3250U-2.3228B1+2.3363B2+0.7495V1-1.0865G
C
C       They are sensitive to effective temperature and gravity respectively. 
C
C
C     - For the intermediate range (8000K<Teff<11000K), we use pT and pG which
C       are defined in a similar way as the a and r of the Stromgren photometry
C       by the following formula 
C       
C            pT = B2-V1 + 0.1X + 0.0635(Y-0.2d) - 0.006
C            pG =     Y - 0.2d + 1.1(B2-V1+0.1X) + 0.31
C
C       where X and Y are the parameters defined above, B2-V1 is the temperature
C       parameter for cool stars and d is the reddening-free parameter measuring
C       the Balmer jump. As pT and pG are not independent of interstellar 
C       extinction because of the presence of B2-V1 in these parameters, we
C       have to correct them by the following relations
C
C            pTo = pT - E(B2-V1)
C            pGo = pG - 1.1E(B2-V1)
C            
C
C     - For cool stars, B2-V1 colour index is used as a temperature indicator, 
C       the reddening-free parameter d as a gravity indicator and m2 as a 
C       metallicity indicator. I recall here for convenience the definiton of
C       d and m2
C 
C            d  = U-B1-1.43(B1-B2)
C            m2 = B1-B2-0.457(B2-V1)
C
C     
C
C     grids   : - The file bcdci contains numerical values 
C                 of the grids  { (B2-V1) , d } --> { Teff , log(g) } for 
C                 different metallicities ( +1 +0.5 +0.3 +0.2 +0.1 +0.0
C                 -0.1 -0.2 -0.3 -0.5 -1.0 -1.5 -2.0 -2.5 -3.0 -3.5 -4.0 -4.5
C                 -5.0) with in free format
C                         - [M/H]k  (B2-V1)i  dj  (log(g))ijk  (Te)ijk  
C               - The file bcmci contains numerical values 
C                 of the grids  { (B2-V1) , m2 } --> { [M/H] , Teff } for 
C                 different gravities ( 2.5 3 3.5 4 4.5 5), with in 
C                 free format
C                         - log(g)k  (B2-V1)i  m2j  [M/H]ijk  (Te)ijk   
C               - The files xcycp10i, xcycp00i and xcycm10i
C                 contain numerical values of the grids 
C                 {X,Y} --> {Teff,log(g)} for [M/H] = +1, 0 and -1
C                 respectively, with in free format
C                         - Xi  Yj  (log(g))ij  (theta)ij  (Teff)ij
C               - The files ptpgp10i, ptpgp00i and ptpgm10i
C                 contain numerical values of the grids 
C                 {pT,pG} --> {Teff,log(g)} for [M/H] = +1, 0 and -1
C                 respectively, with in free format
C                      - pTi  pGj  [M/H] (log(g))ij  (theta)ij  (Teff)ij
C
C 
C     the six colours of the Geneva photometric system 
      REAL U,V,B1,B2,V1,G
C
C     parameters for hot and intermediate stars   
      REAL X,Y,pT,pG
C     variables for XYTL and PTGTL
      REAL tef3(3),log3(3),covar(2,2),tef2(3),log2(3)
      REAL dte2(3),dl2(3),extr2(3)
      REAL dte3(3),dl3(3),extr3(3)
      REAL MH,xMH,log
C
C     declaration for calibration (B2-V1,d,m2)    
      REAL m2min,m2pas                                                     
C     Minimal values for B2-V1,d et m2
      PARAMETER (b2v1min=-0.175,dmin=+0.25,m2min=-0.66)  
C     steps in the inverted grids
      PARAMETER (b2v1pas=+0.025,dpas=+0.05,m2pas=+0.02)
C     Maximal number of values  B2-V1, d et m2. Nv > MAX(...)
      PARAMETER (Nb2v1=50,Nd=60,Nm2=80,Nv=80)   
C     Number of values for log(g) and [M/H]                          
      PARAMETER (Nm=19,Ng=6)
C     Maximal deviation allowed during the iterative process
      PARAMETER (dtmax=20,dgmax=0.02,dmmax=0.02,Nitermax=6)               
      PARAMETER (Nwarn=2)
      PARAMETER (NT1=Nd*Nb2v1*Nm, NT2=Nm2*Nb2v1*Ng) 
      REAL Vb2v1(Nb2v1),Vd(Nd),Vm2(Nm2)                                   
      REAL Vm(Nm),Vg(Ng)
C     Physical parameters and indices                                                  
      REAL teff,logg,metal,b2v1,d,m2                                
C     Scatters for Te,log(g),[M/H]
      REAL st,sg,sm                                                 
C     Deviations in Te, [M/H] and log(g) during the iterative process
      REAL dt,dm,dg                                                       
C     Temporary values during the iterative process
      REAL teff1,teff2,logg1,logg2,metal1,metal2                          
      REAL Mt(Nd,Nb2v1,Nm),Mg(Nd,Nb2v1,Nm),Mm(Nm2,Nb2v1,Ng)               
      LOGICAL exist1(Nd,Nb2v1,Nm),exist2(Nm2,Nb2v1,Ng)                  
      CHARACTER cas*1,choix*1,data*20,titre*80,fileres*20
C     Remarks NB and SSI          
      CHARACTER*3 warn(Nwarn)
      REAL tr(Nv,7),ext,exg,exm 
      DATA exist1/NT1*.FALSE./,exist2/NT2*.FALSE./
C     Values of [M/H] for the grids (d,B2-V1)
      DATA Vm/-5.,-4.5,-4.,-3.5,-3.,-2.5,-2.,-1.5,-1,-0.5,
     &-0.3,-0.2,-0.1,0.,0.1,0.2,0.3,0.5,1./
C     Values of log(g) for the grids (m2,B2-V1)
      DATA Vg/2.5,3.0,3.5,4.0,4.5,5./   
      DATA p/0/,val/4/                                     
C
C                
      OPEN(15,FILE='bcdci',STATUS='old')
      OPEN(16,FILE='bcmci',STATUS='old')      
C
C     Definition des vecteurs Vb2v1, Vd et Vm2

      DO 1 i=1,Nb2v1
         Vb2v1(i)=b2v1pas*(i-1)+b2v1min
    1 CONTINUE
      DO 2 i=1,Nd
         Vd(i)=dpas*(i-1)+dmin
    2 CONTINUE
      DO 3 i=1,Nm2
         Vm2(i)=m2pas*(i-1)+m2min
    3 CONTINUE
C
C
C     Reading of the file bcdci and bcmci
C
    4 READ(15,*,END=6)metal,b2v1,d,logg,teff
         i=NINT(1+(b2v1-b2v1min)/b2v1pas)  
         j=NINT(1+(d-dmin)/dpas)
         DO 5 n=1,Nm
            IF (metal.EQ.Vm(n)) k=n
    5    CONTINUE
         Mt(j,i,k)=teff
         Mg(j,i,k)=logg
         exist1(j,i,k)=.TRUE.
      GOTO 4
    6 CONTINUE
C
    7 READ(16,*,END=9)logg,b2v1,m2,metal
         i=NINT(1+(b2v1-b2v1min)/b2v1pas)  
         j=NINT(1+(m2-m2min)/m2pas)
         DO 8 n=1,Ng
            IF (logg.EQ.Vg(n)) k=n
    8    CONTINUE
         Mm(j,i,k)=metal
         exist2(j,i,k)=.TRUE.
      GOTO 7
    9 CONTINUE
C
C
C     Reading of the Geneva colours on screen for each star or in a file.
C     If you use a file, the first line of it has to contain a title.
C
C
      PRINT*
      WRITE(6,95),' Would you work star by star (y/n) ? '
  95  FORMAT(a,$)
      READ(*,200)cas
      IF (cas.NE.'y') THEN
         PRINT*
         PRINT*,'Name of the input file which contains on each line identification' 
         PRINT*,'of the star, U V B1 B2 V1 G, E(B2-V1) and [M/H]=-1,0,1 (in free format)'
         READ(*,201)data
         PRINT*,'Name of the output file'
         READ(*,201) fileres
         OPEN(30,FILE=data,STATUS='old')
         OPEN(31,FILE=fileres,STATUS='unknown')
         READ(30,98) titre
         WRITE(31,98) titre
         WRITE(31,*)
         WRITE(31,*) ' (1) : calibration in grids (B2-V1,d) or (B2-V1,m2)'
         WRITE(31,*) ' (2) : calibration in grids (pT,pG)'
         WRITE(31,*) ' (3) : calibration in grids (X,Y)'   
         WRITE(31,*)    
         WRITE(31,*) ' SSI : Physical parameters are determined  with only one'   
         WRITE(31,*) '       iteration in the grids for cool stars'   
         WRITE(31,*) ' NC  : Physical paramaters do not converge for cool stars'   
         WRITE(31,*) '       In this case, the process is stopped after 6 iterations'   
         WRITE(31,*)    
         WRITE(31,*) ' If there is no extrapolation in the grids (Te,logg), (Te,[M/H]), we display the value 0.0'    
         WRITE(31,*) ' If the value of extrapolation is larger than 10, the value of physical parameters become absurd'
         WRITE(31,*) ' In this case, we change absurd values to 0.'
         WRITE(31,*)    
  98     FORMAT(A80)   
         WRITE(31,*)
         WRITE(31,500)
  500    FORMAT('Star  ',9X,'Teff',18X,'Logg',18X,'[M/H]',
     &          13x,'Remarks  ',3x,'Extrapolation',/,89x,'(Te,logg)',
     &          x,'(Te,[M/H])')     
	  IF (p.EQ.0) then
              PRINT*,' '
              PRINT*,'Which grids would you want to use ?'
              PRINT*,' '
              PRINT*,'1) calibration in grids (b2v1,d,m2)'
              PRINT*,' '
              PRINT*,'2) calibration in grids (pT,pG)'
              PRINT*,' '
              PRINT*,'3) calibration in grids (X,Y)'
              PRINT*,' '
              PRINT*,'4) calibration in grids choosen automatically by the code'
              PRINT*,' '
              READ(*,*)val
          p=p+1
         ENDIF
      ENDIF
C      
   10 IF (cas.EQ.'y') THEN
         PRINT*
         PRINT*,'Could you give me the 6 colours of Geneva U,V,B1,B2,V1 and G (U=10 to finish)?'
         READ(*,*) U,V,B1,B2,V1,G
         IF (U.EQ.10) GOTO 90
         PRINT*,'Colour excess E(B2-V1) ?'
         READ(*,*) Eb2v1
         PRINT*,'Metallicity of the star ([M/H]=-1,0,1) ?'
         READ(*,*) MH
	 PRINT*,' '
      ELSE
         READ(30,*,END=90)ident,U,V,B1,B2,V1,G,Eb2v1,MH 
      ENDIF
C
C     B2-V1,d,m2,X,Y,pT,pG
C
      b2v1 = B2-V1
      d = U-B1-1.43*(B1-B2)
      m2 = B1-B2-0.457*(B2-V1)
      X = 0.3788+1.3764*U-1.2162*B1-0.8498*B2-0.1554*V1+0.8450*G
      Y = -0.8288+0.3235*U-2.3228*B1+2.3363*B2+0.7495*V1-1.0865*G
      pT = b2v1+0.1*X+0.0635*(Y-0.2*d)-0.006
      pG = Y-0.2*d+1.1*(b2v1+0.1*X)+0.31
C
C     parameters are corrected by a possible interstellar extinction
C

      b2v1 = b2v1 - Eb2v1
      pT = pT - Eb2v1
      pG = pG - 1.1*Eb2v1
C
      if (cas.EQ.'y') THEN
	write(*,27)b2v1,d,m2
   27 format('B2-V1 = ',f6.3,3x,'d = ',f6.3,3x,'m2 = ',f6.3,/)
	write(*,28)pT,pG
   28 format('pT = ',f6.3,3x,'pG = ',f6.3,/)
	write(*,29)X,Y
   29 format('X = ',f6.3,3x,'Y = ',f6.3,/)
      PRINT*,' '
      PRINT*,'Which grids would you want to use ?'
      PRINT*,' '
      PRINT*,'1) calibration in grids (b2v1,d,m2)'
      PRINT*,' '
      PRINT*,'2) calibration in grids (pT,pG)'
      PRINT*,' '
      PRINT*,'3) calibration in grids (X,Y)'
      PRINT*,' '
      PRINT*,'4) calibration in grids choosen automatically by the code'
      PRINT*,' '
      READ(*,*)val
      endif
C
C
C 
C     determination of the grids
C
       if (val.eq.4) then
          test1=X-1.4
          test2=Y+0.04
          test3=b2v1+0.029
          test4=X-1.3
          test5=Y+0.06
       else if (val.eq.3) then
          test1=-1
          test2=1
          test4=-1
          test5=1
          go to 72
       else if (val.eq.2) then
          test1=1
          test2=-1
          test3=-1
          test4=1
          test5=-1
          go to 81
       else if (val.eq.1) then
          test1=1
          test2=-1
          test3=1
          test4=1
          test5=-1
          goto 82
       end if

C     calibration (X,Y)
C
   72 IF (test1.le.0.and.test2.ge.0.or.test4.le.0.and.test5.ge.0) THEN

        if(MH.eq.0.or.MH.eq.-1.or.MH.eq.1) then
            CALL XYTL(X,Y,MH,thef,log,dt,dl,covar,extr)
               extr3(3)=extr
            if (extr3(3).le.10) then
               tef3(3)=5040./thef
               dte3(3)=5040.*dt/(thef*thef)
               log3(3)=log
	       dl3(3)=dl
               rho=0
            else
               tef3(3)=0
               dte3(3)=0
               log3(3)=0
               dl3(3)=0
               rho=0
            endif

        else

            xMH=0
            do 56 i=1,2
            CALL XYTL(X,Y,xMH,thef,log,dt,dl,covar,extr)
               extr3(i)=extr
               if (extr3(i).le.10) then
                  tef3(i)=5040./thef   
                  dte3(i)=5040.*dt/(thef*thef)
                  log3(i)=log
                  dl3(i)=dl
               else
                  tef3(i)=0
                  dte3(i)=0
                  log3(i)=0
                  dl3(i)=0
               endif
               if (MH.ge.0) then 
                  xMH=1
               else 
                  xMH=-1
               endif
   56     enddo
               if (extr3(1).le.10.and.extr3(2).le.10) then
                  tef3(3)= tef3(1)*(1-abs(MH))+tef3(2)*abs(MH)
                  dte3(3)=dte3(1)*(1-abs(MH))+dte3(2)*abs(MH)
                  log3(3)= log3(1)*(1-abs(MH))+log3(2)*abs(MH)
                  dl3(3)=dl3(1)*(1-abs(MH))+dl3(2)*abs(MH)
                  extr3(3)=extr3(1)*(1-abs(MH))+extr3(2)*abs(MH)
               else
                  tef3(3)=0
                  dte3(3)=0
                  log3(3)=0
                  dl3(3)=0
               endif
        endif
        goto 99
      ENDIF

 
C
C     calibration (pT,pG)
C
   81 IF (test3.le.0) THEN

        if(MH.eq.0.or.MH.eq.-1.or.MH.eq.1) then

            CALL PTGTL(pT,pG,MH,thef,log,dt,dl,covar,extr)
            extr2(3)=extr   
            if (extr2(3).le.10) then
               tef2(3)=5040./thef
               dte2(3)=5040.*dt/(thef*thef)
               log2(3)=log
	       dl2(3)=dl
               rho=0
            else
               tef2(3)=0
               dte2(3)=0
               log2(3)=0
               dl2(3)=0
               rho=0
            endif

        else

            xMH=0
            do 55 i=1,2
               CALL PTGTL(pT,pG,xMH,thef,log,dt,dl,covar,extr)
               extr2(i)=extr
               if (extr2(i).le.10) then
                  tef2(i)=5040./thef    
                  dte2(i)=5040.*dt/(thef*thef)
                  log2(i)=log
                  dl2(i)=dl
               else
                  tef2(i)=0
                  dte2(i)=0
                  log2(i)=0
                  dl2(i)=0
               endif
               if (MH.ge.0) then 
                  xMH=1
               else 
                  xMH=-1
               endif
   55     enddo
               if (extr2(1).le.10.and.extr2(2).le.10) then
                  tef2(3)= tef2(1)*(1-abs(MH))+tef2(2)*abs(MH)
                  dte2(3)=dte2(1)*(1-abs(MH))+dte2(2)*abs(MH)
                  log2(3)= log2(1)*(1-abs(MH))+log2(2)*abs(MH)
                  dl2(3)=dl2(1)*(1-abs(MH))+dl2(2)*abs(MH)
                  extr2(3)=extr2(1)*(1-abs(MH))+extr2(2)*abs(MH)
               else
                  tef2(3)=0
                  dte2(3)=0
                  log2(3)=0
                  dl2(3)=0
               endif
        endif
        goto 99
      ENDIF


C
C     calibration (B2-V1,d.m2)
C
C     First iteration : we suppose that [M/H]=0 for the determination of the 
C     temperature effective and gravity and that log(g)=4 for the determination
C     of [M/H]
C
   82 IF (test3.gt.0) THEN
C
C     Stars cooler than B2-V1>0.5 are not calibrate
C
      IF (b2v1.gt.0.5) THEN  
         IF (cas.EQ.'y') THEN  
            PRINT*,'This star is out of the range of validity :(B2-V1) > 0.5 '
	    PRINT*
         ELSE
            WRITE(31,87)ident
   87       FORMAT(I8,2X,'This star is out of the range of validity :(B2-V1) > 0.5')  
         ENDIF
       GOTO 10
      ENDIF
      DO 18 i=1,Nwarn
         warn(i)=' '
   18 CONTINUE
      iter=0
C     Teff
      CALL BISPLI(d,b2v1,Vd,Vb2v1,Mt(1,1,14),exist1(1,1,14),Nd,Nb2v1,       
     &   tr,Nv,teff,dfdy,dfdx,ext)
C     Sigma(teff)
      st = SIGMA(3,dfdx,dfdy,0,0,0)                                       
C     log(g)
      CALL BISPLI(d,b2v1,Vd,Vb2v1,Mg(1,1,14),exist1(1,1,14),Nd,Nb2v1,       
     &   tr,Nv,logg,dfdy,dfdx,exg)
C     Sigma(log(g))
      sg = SIGMA(3,dfdx,dfdy,0,0,0)                                       
C     [M/H]
      CALL BISPLI(m2,b2v1,Vm2,Vb2v1,Mm(1,1,4),exist2(1,1,4),Nm2,          
     &   Nb2v1,tr,Nv,metal,dfdy,dfdx,exm)
C     Sigma([M/H])
      sm = SIGMA(2,dfdx,dfdy,0,0,0) 
C
C     We stop here if the uncertainty of [M/H] is too large, because it can bend
C     the other parameters via the iterative process
C
      IF (b2v1.le.0.025) THEN
         warn(1)='SSI'
         GOTO 15
      ENDIF
C                                      
C     additional iteration
C
      DO 16 iter=1,100                                                     
         dt=teff
         dg=logg
         dm=metal
C
C     Determination of the grids with metallicity between the initial value of [M/H]
C
         IF (metal.LT.Vm(1)) THEN                                         
            imetal=1
            GOTO 12
         ELSE IF (metal.GT.Vm(Nm)) THEN
            imetal=Nm-1
            GOTO 12
         ENDIF
         DO 11 i=1,Nm-1
            IF (metal.GE.Vm(i).AND.metal.LE.Vm(i+1)) THEN
               imetal=i
               GOTO 12
            ENDIF
   11    CONTINUE     
C        
   12    fracm=(metal-Vm(imetal))/(Vm(imetal+1)-Vm(imetal))               
C
         CALL BISPLI(d,b2v1,Vd,Vb2v1,Mt(1,1,imetal),
     &    exist1(1,1,imetal),Nd,Nb2v1,tr,Nv,teff1,dfdy1,dfdx1,ext)
         CALL BISPLI(d,b2v1,Vd,Vb2v1,Mt(1,1,imetal+1),
     &    exist1(1,1,imetal+1),Nd,Nb2v1,tr,Nv,teff2,dfdy2,dfdx2,ext)
         teff = teff1 + fracm*(teff2-teff1)
         st = SIGMA(1,dfdx1,dfdy1,dfdx2,dfdy2,fracm)
C
         CALL BISPLI(d,b2v1,Vd,Vb2v1,Mg(1,1,imetal),
     &    exist1(1,1,imetal),Nd,Nb2v1,tr,Nv,logg1,dfdy1,dfdx1,exg)
         CALL BISPLI(d,b2v1,Vd,Vb2v1,Mg(1,1,imetal+1),
     &    exist1(1,1,imetal+1),Nd,Nb2v1,tr,Nv,logg2,dfdy2,dfdx2,exg)
         logg = logg1 + fracm*(logg2-logg1)
         sg = SIGMA(1,dfdx1,dfdy1,dfdx2,dfdy2,fracm)
C
C        Determination of the grids with gravity between the initial value of logg
C
         IF (logg.LT.Vg(1)) THEN
            ilogg=1
            GOTO 14
         ELSE IF (logg.GT.Vg(Ng)) THEN
            ilogg=Ng-1
            GOTO 14
         ENDIF
         DO 13 i=1,Ng-1
            IF (logg.GE.Vg(i).AND.logg.LE.Vg(i+1)) THEN
               ilogg=i
               GOTO 14
            ENDIF
   13    CONTINUE  
C   
   14    fracg=(logg-Vg(ilogg))/(Vg(ilogg+1)-Vg(ilogg))                   
C        Determination de [M/H], ainsi que du sigma par interpolation entre 
C        les deux grilles (m2,B2-V1) de log(g) differents
C
         CALL BISPLI(m2,b2v1,Vm2,Vb2v1,Mm(1,1,ilogg),
     &    exist2(1,1,ilogg),Nm2,Nb2v1,tr,Nv,metal1,dfdy1,dfdx1,exm)
         CALL BISPLI(m2,b2v1,Vm2,Vb2v1,Mm(1,1,ilogg+1),
     &    exist2(1,1,ilogg+1),Nm2,Nb2v1,tr,Nv,metal2,dfdy2,dfdx2,exm)
         metal=metal1+fracg*(metal2-metal1)
         sm = SIGMA(2,dfdx1,dfdy1,dfdx2,dfdy2,fracg)
C
C        If the differences between the last two iterations is less than 20 K for
C        temperature, 0.02 dex for gravity and metallicity, we stop the process here
C
         dt=ABS(dt-teff)
         dg=ABS(dg-logg)
         dm=ABS(dm-metal)
         IF (dt.LE.dtmax.AND.dg.LE.dgmax.AND.dm.LE.dmmax) GOTO 17
C        Do not converge after Nitermax iterations
         IF (iter.GE.Nitermax) THEN                                       
            IF (cas.EQ.'y') THEN
               IF (pT-0.1.gt.0) then
                   WRITE(*,103)dm,dt,dg
  103              FORMAT(' NC dm=',F4.2,' dt=',F5.0,' dg=',F4.2)
                   PRINT*,' On continue ? (O/N)'
                   READ(*,200)choix
               ENDIF
               IF (choix.NE.'o'.AND.choix.NE.'O') GOTO 17
            ELSE
               warn(2)='NC'                                               
               GOTO 17
            ENDIF
         ENDIF
C
   16 CONTINUE
   17 CONTINUE
   15 CONTINUE
C
C        if the value of the extrapolation (ext in the grid d vs B2-V1 and exm in
C        the grid m2 vs B2-V1) is larger than 10, we obtain absurd values for Teff,
C        logg and [M/H]. For this reason, I prefer changing absurd values to 0.
C 
	IF (ext.gt.10) THEN
           teff=0.
           st=0.
           logg=0.
           sg=0.
        ENDIF
	IF (exm.gt.10) THEN
           metal=0.
           sm=0.
        ENDIF
C        
C        
C
      ENDIF
C
   99 CONTINUE
C
C     Results
C
      IF (cas.EQ.'y') THEN      
         IF (test1.le.0.and.test2.ge.0.or.test4.le.0.and.test5.ge.0) then
            PRINT*,'temperature and gravity are estimated'
            PRINT*,'by X and Y'    
      	    WRITE(*,75) ifix(tef3(3)),ifix(dte3(3)),log3(3),dl3(3),extr3(3)
   75       FORMAT(/' Teff = ',I6,' +-'I5,/, 
     &             ' Log(g) = ',F6.2,' +-',F5.2,/,' extrap(X,Y)',F6.1,/)  
         ELSE IF (test3.lt.0) THEN
      	    PRINT*,'temperature and gravity are estimated' 
            PRINT*,'by pT and pG' 
      	    WRITE(*,73)ifix(tef2(3)),ifix(dte2(3)),log2(3),dl2(3),extr2(3)
   73       FORMAT(/' Teff = ',I6,' +-'I5,/,
     &          ' Log(g) = ',F6.2,' +-',F5.2,/,' extrap(pT,pG)',F6.1,/)
         ELSE IF (test3.ge.0) THEN       
            PRINT*,'temperature, gravity and metallicity are estimated'
            PRINT*,'by B2-V1, d and m2'
            WRITE(*,71) ifix(teff),ifix(st),metal,sm,logg,sg,ext,exm,warn
   71       FORMAT(/' Teff   = ',I6,' +-',I5,/,' [M/H]  = ',F6.2,' +-',F5.2
     &             ,/,' Log(g) = ',F6.2,' +-',F5.2,/,' extrap(B2V1,d)',F7.1,/,
     &             ' extrap(B2V1,m2)',F7.1,/,2A4,/)         
         ENDIF 
      ELSE  
         IF (test1.le.0.and.test2.ge.0.or.test4.le.0.and.test5.ge.0) THEN
      	    WRITE(31,80)ident,ifix(tef3(3)),ifix(dte3(3)),log3(3),dl3(3),MH,extr3(3)
   80       FORMAT(I8,5X,I6,' +-'I5,'(3)',5X, 
     &             F6.2,' +-',F5.2,'(3)',5X,F6.2,23x,F6.1)  
         ELSE IF (test3.lt.0) THEN
      	    WRITE(31,78)ident,ifix(tef2(3)),ifix(dte2(3)),log2(3),dl2(3),MH,extr2(3)
   78       FORMAT(I8,5X,I6,' +-',I5,'(2)',5X,
     &          F6.2,' +-',F5.2,'(2)',5X,F6.2,23x,F6.1)
C
         ELSE IF (test3.ge.0.) THEN       
            WRITE(31,76) ident,ifix(teff),ifix(st),logg,sg,
     &                   metal,sm,warn,ext,exm
   76       FORMAT(I8,5X,I6,' +-',I5,'(1)',5X,F6.2,' +-',
     &             F5.2,'(1)',5X,F6.2,' +-',F5.2,'(1)',2x,2A4,2x,f6.1,4x,f6.1) 
   83	    FORMAT(7f11.3)        
C
         ENDIF
      ENDIF                          
C                
      GOTO 10
C
  200 FORMAT(A1)
  201 FORMAT(A20)
   90 STOP
      END
C
C***************************************************************************
C
      FUNCTION SIGMA(I,dfdx1,dfdy1,dfdx2,dfdy2,frac)
C
C     Variance S = sigma^2 et Covariance observationnelle  
C     sur les couleurs pour un poids P=1 (en E-6 mag^2).
C     (corrigee le 4.3.1993)
      PARAMETER (Sb2v1=40.0,Sd=186.6,Sm2=60.8)                            
      PARAMETER (Sdb2v1=18.6,Sm2b2v1=-31.3)                               
C
C     Cas de Te et log(g)
      IF (I.EQ.1) THEN                                                    
         s1=SQRT(Sb2v1*dfdx1**2+Sd*dfdy1**2+2.*Sdb2v1*dfdx1*dfdy1)
         s2=SQRT(Sb2v1*dfdx2**2+Sd*dfdy2**2+2.*Sdb2v1*dfdx2*dfdy2)
         SIGMA=(s1 + frac*(s2-s1))/1000.0
C     Cas de [M/H]
      ELSE IF (I.EQ.2) THEN                                               
         s1=SQRT(Sb2v1*dfdx1**2+Sm2*dfdy1**2+2.*Sm2b2v1*dfdx1*dfdy1)
         s2=SQRT(Sb2v1*dfdx2**2+Sm2*dfdy2**2+2.*Sm2b2v1*dfdx2*dfdy2)
         SIGMA=(s1 + frac*(s2-s1))/1000.0
C     Cas de la masse, ou Iter0 pour Te et logg SEULS
      ELSE IF (I.EQ.3) THEN                                               
         s=SQRT(Sb2v1*dfdx1**2+Sd*dfdy1**2+2.*Sdb2v1*dfdx1*dfdy1)            
         SIGMA=s/1000.0
      ENDIF
      RETURN
      END
C
C***************************************************************************
C
      SUBROUTINE BISPLI(Tt,Uu,   Tm,Um,Xm,Loxm,Nt,Nu,  Vm,Nv                    
     &                 ,Xx,Dxt,Dxu,   Extrap)                                   
      save debut
C                                                                               
C        Routine d'interpolation pour fonction a deux variables definie         
C     par Xm(it,iu) en certains points d'un reseau [Tm(it), Um(iu)].            
C        La variable LOGICAL Loxm(it,iu) dit si Xm(it,iu) est                 
C     significatif. Si on se trouve a l'interieur de la grille des points       
C     significatifs, Extrap vaut 0. Si Extrap excede .5 ou 1. il faut           
C     se mefier.                                                                
C        Xx est la valeur de la fonction en (Tt,Uu).                            
C        Dxt et Dxu sont des derivees partielles.                               
C                                                                               
C     On donne:                                                                 
C     =========                                                                 
C        - Tt,     Uu      Les deux variables independantes.                    
C        - Tm(Nt), Um(Nu)  Grille des variables independantes                   
C        - Xm(Nt,Nu)       Valeurs (eventuellement bidon) de la                 
C                          fonction en (Tm., Um.)                               
C        - Loxm(Nt,Nu)*1   Matrice logique indiquant la signification           
C                          des Xm correspondants.                               
C        - Nt, Nu          Dimensions de Tm, Um, Xm et Loxm                     
C        - Vm(Nv,7)        Matrice de travail                                   
C        - Nv              Dimension de Vm. Nv doit etre .GE.Nt et .GE.Nv       
C                                                                               
C     On recoit:                                                                
C     ==========                                                                
C        - Xx              Valeur de la fonction en (Tt, Uu)                    
C        - Dxt             Derivee partielle p. r. a Tt en (Tt, Uu)             
C        - Dxu             Derivee partielle p. r. a Uu en (Tt, Uu)             
C        - Extrap = 0.     (Tt, Uu) se trouve a l'interieur du                  
C                          reseau de (Tm., Um.) significatifs.                  
C                          Sinon nombre de mailles d'extrapolation.             
C                                                                               
      REAL Tm(Nt),Um(Nu),Xm(Nt,Nu),Vm(Nv,7)                                     
      LOGICAL incid                                                             
      LOGICAL debut, Loxm(Nt,Nu)                                              
      COMMON/EXTREM/Icb,Ypb,Ich,Ysh                                             
      DATA debut/.TRUE./                                                        
C                                                                               
      IF(debut)THEN                                                             
         debut=.FALSE.                                                          
         Icb=0                                                                  
         Ich=0                                                                  
      END IF                                                                    
                                                                               
      mu=0                                                                      
      Extrap=0.                                                                 
      ju=NRCORD(Uu,Um,Nu,incid)                                                 
      DO 9 iu=1,Nu                                                              
C ----------- Succession d'interpolations sur Tt -----------                    
         mt=0                                                                   
         DO 8 it=1,Nt                                                           
            IF(loxm(it,iu))THEN                                                 
               mt=mt+1                                                          
               Vm(mt,1)=Tm(it)                                                  
               Vm(mt,2)=Xm(it,iu)                                               
            END IF                                                              
    8    CONTINUE                                                               
C ------------ Pas assez de points significatifs -----------------              
         IF(mt.LE.1)GO TO 9                                                     
C ------------ Extrapolation trop "violente" ------------------                 
         extrab=(Tt-Vm(1 ,1))/(Vm(1 ,1)-Vm(2   ,1))                             
         extrah=(Tt-Vm(mt,1))/(Vm(mt,1)-Vm(mt-1,1))                             
         extral=0.                                                              
         IF(iu.EQ.ju-1.OR.iu.EQ.ju+2)extral=1.                                  
         IF(iu.EQ.ju  .OR.iu.EQ.ju+1)extral=2.                                  
         IF(extrab.GE.extral)GO TO 9                                            
         Extrap=MAX(Extrap,extrab)                                              
         IF(extrah.GE.extral)GO TO 9                                            
         Extrap=MAX(Extrap,extrah)                                              
C ------------ Cas normal -----------------                                     
         mu=mu+1                                                                
         Vm(mu,3)=Um(iu)                                                        
         CALL SPLINF(Tt,Vm(1,1),Vm(1,2),mt,   Vm(1,6),Vm(1,7)                   
     &              ,Vm(mu,4),Vm(mu,5),xs)                                      
    9 CONTINUE                                                                  
      IF(mu.EQ.0)THEN                                                           
C ------------ Grille trop peu fournie --------------                           
         Extrap=1000.                                                           
         Xx =0.                                                                 
         Dxt=0.                                                                 
         Dxu=0.                                                                 
      ELSE IF(mu.EQ.1)THEN                                                      
         Extrap=100.                                                            
         Xx =Vm(1,4)                                                            
         Dxt=Vm(1,5)                                                            
         Dxu=0.                                                                 
      ELSE                                                                      
         Extrap=MAX(Extrap,(Uu-Vm(1 ,3))/(Vm(1 ,3)-Vm(2   ,3)))                 
         Extrap=MAX(Extrap,(Uu-Vm(mu,3))/(Vm(mu,3)-Vm(mu-1,3)))                 
         CALL SPLINF(Uu,Vm(1,3),Vm(1,4),mu,Vm(1,6),Vm(1,7),Xx,Dxu,xs)           
         CALL SPLINF(Uu,Vm(1,3),Vm(1,5),mu,Vm(1,6),Vm(1,7),Dxt,xp,xs)           
      END IF                                                                    
      END                                                                       
C
C***************************************************************************
C
      FUNCTION NRCORD(Rx,Rv,Nb,Incid)
C
C        NRCORD sert a trouver quelle composante du vecteur ordonne Rv
C     est egale (alors Incid est .TRUE.) ou suit immediatement Rx. 
C
C     On donne:
C     =========
C        - Rx  Nombre a situer
C        - Rv  Vecteur ordonne
C        - Nb  Nombre de composantes
C
C     On recoit:
C     ==========
C        - NRCORD 
C        - Incid  Verite de Rv(NRCORD)=Rx
C     Declarer LOGICAL Incid et DIMENSION Rv(nn) ou nn  >= Nb.
C
      LOGICAL Incid
      DIMENSION Rv(Nb)
C
      rdif=Rv(Nb)-Rv(1)
      rsig=SIGN(1.,rdif)
      NRCORD=(Nb+1)/2
      IF(rdif.EQ.0)GO TO 30
C --------- Cas ou Rx sort de l'intervalle definit par Rv(1), Rv(Nb) -------
      IF((Rx-Rv(Nb))*rsig.GE.0.)THEN
         NRCORD=Nb
      ELSE IF((Rx-Rv(1))*rsig.LT.0.)THEN
         NRCORD=0
      ELSE
C ------------------- Cas ou Rv a peu de composantes -----------------------
         IF(Nb.LE.30)THEN
            DO 19 ib=1,Nb
               IF((Rx-Rv(ib))*rsig)11,12,19
   11             NRCORD=ib-1
               GO TO 30
   12             NRCORD=ib
               GO TO 30
   19       CONTINUE
C --------------------------- Grand vecteur --------------------------------
         ELSE
            nbm=Nb-1
            nit=NRCORD
C
   20       nita=nit
            nit=(nit+1)/2
            IF(rsig*(Rv(NRCORD)-Rx))21,30,22
   21          IF(nita.LE.1)GO TO 30
               NRCORD=MIN(nbm,NRCORD+nit)
            GO TO 20
   22          NRCORD=MAX(  1,NRCORD-nit)
            IF(nita.GT.1)GO TO 20
         END IF
      END IF
   30 IF(NRCORD.GE.1.AND.NRCORD.LE.Nb)THEN
         Incid=Rx.EQ.Rv(NRCORD)
      ELSE
         Incid=.FALSE.
      END IF
      END
C 
C***************************************************************************
C
      SUBROUTINE SPLINF(Xx,Xd,Yd,Nn,   Ypd,Ysd,   Yy,Yp,Ys)
C
C        Procedure d'interpolation de type bicubique a travers les
C     points P(i) [Xd(i),Yd(i)], i=1,Nn ou les Xd(i) forment une suite
C     strictement monotone.
C        On a une fonction Yy = F(Xx). Evidemment F(Xd(i)) = Yd(i).
C        F(Xx) est formee d'arcs de fonctions du 3eme deg, elle est
C     partout continue et 2 fois derivable. La premiere et
C     la seconde  derivee sont partout continues. La valeur de la
C     premiere (resp. seconde) derivee en Xx est Yp (resp. Ys).
C        L'algorithme de base est le meme que SPLINE et SPLINT dans
C     les "NUMERICAL RECIPES".
C
C     On donne:
C     ---------
C        - Xx        Valeur de la variable independante
C        - Xd(Nn)    Abscisses des P(i)
C        - Yd(Nn)    Ordonnees  "   "
C        - Nn        Nombre de points P(i)
C
C     On fournit:
C     -----------
C        - Ypd(Nn)   De la memoire qui contiendra les F'(Xd(i))
C        - Ysd(Nn)   De la memoire qui contiendra les F"(Xd(i))
C
C        Si on fait plusieurs interpolations avec le meme jeu de P(i)
C     On gagne un temps de calcul considerable en utilisant SPLINU
C     des le deuxieme appel. SPLINF calcule Ypd et Ysd tandis que
C     SPLINU utilise les Xpd et Xsd obtenus par SPLINF.
C
C     On recoit:
C     ----------
C        - Yy        Valeur de la fonction interpolee en Xx
C        - Yp        Valeur de la premiere derivee en Xx
C        - Ys        Valeur de la seconde  derivee en Xx.
C
C     COMMON/EXTREM/Icb,Ypsb,Ich,Ypsh
C     -------------------------------
C        L'algorithme laisse deux degres de liberte au probleme.
C     En pratique il s'agira toujours de conditions aux limites portant
C     soit sur la premiere, soit sur la seconde derivee.
C        Icb (resp. Ich) indique si on veut agir sur la premiere ou
C     la seconde derivee en Xx(1) (resp.(Xx(Nn)). Ypsb/h donne la valeur
C     que l'on veut eventuellement imposer a ladite derivee.
C
C        - Icb/h = 0 Valeur par defaut. La seconde derivee est nulle
C                    a l'extremite consideree.
C                    Option recommandee si on veut extrapoler. Prudence alors!
C                = 1 Premiere derivee nulle.
C                    Option recommandee si la fonction est constante
C                    au dela ou en deca de l'extremite consideree
C                    Si Icb/h = 0 ou 1, Ipsb/h = 0.
C                = 2 On impose la valeur Ipsb/h a la seconde  derivee
C                    en Xx(1/Nn).
C                = 3 On impose la valeur Ipsb/h a la premiere derivee
C                    en Xx(1/Nn).
C
      REAL Xd(Nn),Yd(Nn),Ysd(Nn),Ypd(Nn)
      LOGICAL incid
      COMMON/EXTREM/Icb,Ypsb,Ich,Ypsh
C
      IF(Nn.LE.1)GO TO 20
      dxp=Xd(2)-Xd(1)
      IF(dxp.EQ.0.)THEN
         PRINT*,'SPLINF: Xd(1)=Xd(2)=',Xd(1)
         RETURN
      END IF
      ypp=(Yd(2)-Yd(1))/dxp
      IF(Nn.EQ.2)GO TO 20
C
C ---------------- Calcul des y"(Xd(i)) ------------------
C
      IF(Icb.LE.0.OR.Icb.GT.3)Icb=0
      IF(Ich.LE.0.OR.Ich.GT.3)Ich=0
      IF(Icb.LE.1)Ypsb=0.
      IF(Ich.LE.1)Ypsh=0.
C ----------- Condition initiale en Xx(1) ---------------
      IF(MOD(Icb,2).EQ.0)THEN
C ---------- y"(Xx(1)) est definie -----------------
         Ysd(1)=Ypsb
         Ypd(1)=0.
      ELSE
C ---------- y'(Xx(1)) est definie -----------------
         Ysd(1)=3.*(ypp-Ypsb)/dxp
         Ypd(1)=-.5
      END IF
C ---------- Recurrences "aller" -------------------
      DO 11 i=2,Nn-1
         dxm=dxp
         dxp=Xd(i+1)-Xd(i)
         IF(dxp*dxm.LE.0.)THEN
            PRINT 101,(j,Xd(j),j=i-1,i+1)
  101       FORMAT(' SPLINF',3(' Xd(',I2,') =',E12.5),' mal ordonnes')
            RETURN
         END IF
         dxx=dxp+dxm
         rap=dxm/dxx
         ypm=ypp
         ypp=(Yd(i+1)-Yd(i))/dxp
         ys =2.*(ypp-ypm)/dxx
         den=2.+rap*Ypd(i-1)
         Ysd(i)=(3.*ys-rap*Ysd(i-1))/den
         Ypd(i)=(rap-1.)/den
   11 CONTINUE
C ----------- Condition initiale en Xx(Nn) ---------------
         den=2.+Ypd(Nn-1)
      IF(MOD(Ich,2).EQ.0)THEN
C ---------- y"(Xx(Nn)) est definie -----------------
         Ysd(Nn)=Ypsh
         Ypd(Nn)=ypp+(den*Ypsh+Ysd(Nn-1))*dxp/6.
      ELSE
C ---------- y'(Xx(Nn)) est definie -----------------
         Ypd(Nn)=Ypsh
         ys =2.*(Ypsh-ypp)/dxp
         Ysd(Nn)=((Ypsh-ypp)*6./dxp-Ysd(Nn-1))/den
      END IF
C -------- Calcul des y"(i) par resolution d'un systeme tridiagonal
C               recurrence "retour" ---------------------------------
      DO 12 i=Nn-1,1,-1
         Ysd(i)=Ysd(i)+Ypd(i)*Ysd(i+1)
         Ypd(i)=Ypd(i+1)-(Xd(i+1)-Xd(i))*(Ysd(i)+Ysd(i+1))/2.
   12 CONTINUE
C
      ENTRY SPLINU(Xx,Xd,Yd,Nn,Ypd,Ysd,   Yy,Yp,Ys)
C
C ---------------- Cas avec peu de points --------------------
   20 IF(Nn.EQ.1)THEN
         Ys=0.
         Yp=0.
         Yy=Yd(1)
      ELSE IF(Nn.EQ.2)THEN
         Ys=0.
         Yp=(Yd(2)-Yd(1))/(Xd(2)-Xd(1))
         Yy=Yd(1)+(Xx-Xd(1))*Yp
      ELSE
         im=NRCORD(Xx,Xd,Nn,incid)
         im=MAX(1,MIN(Nn-1,im))
         in=im+1
         dx=Xd(in)-Xd(im)
         tim=(Xx-Xd(im))/dx
         tin=tim-1.
         ysm=Ysd(im)
         ysn=Ysd(in)
         Ys=ysn    *tim-ysm    *tin
         Yp=Ypd(in)*tim-Ypd(im)*tin+(ysn-ysm   )*tim*tin*dx   /2.
         Yy=Yd (in)*tim-Yd (im)*tin+(ysn+ysm+Ys)*tim*tin*dx*dx/6.
      END IF
      END
C
C****************************************************************************   
C                                                    
      SUBROUTINE PTGTL(Xx,Yy,Ic,  Thef,Logg,  Dthef,Dlogg,Covar,Extrap)
C
C        Calcul de Theta effective (Thef = 5040./tef1) et de Log(g)
C     a partir de pT et pG, parametres pour les etoiles de temperature
C     intermediaire.[Fe/H] est suppose connu, egal a 0 (solaire), -1 ou +1.
C     La routine donne une estimation des erreurs a partir de la
C     matrice de covariance sur les couleurs de Geneve. Ces erreurs
C     s'acroissent rapidement si on sort de la region ou la fonction
C     F: (Thef, Log(g) |----> (pT, pG) n'est plus bijective.
C        Dans ladite region Extrap est nulle. Ailleurs l'extrapolation
C     presente des risques, du moins si Extrap depasse 0.5 a 1.
C
C     On donne :  Xx,    Yy    Parametres de T intermediaire
C                 Ic           +1, 0, -1 [Fe/H]
C     ==========
C     On recoit:  Thef,  Logg  5040./Teff et log(g)
C     ========== Dthef, Dlogg  erreurs
C                Covar(2,2)    Matrice de covariance des erreurs
C                   Extrap     0. si on est dans la grille de definition
C                              nombre de mailles d'extrap. sinon.
C
      save debut
      PARAMETER (Nx=60,Ny=50,Nxy=Nx*Ny,Nv=60)
      REAL xm(Nx),ym(Ny),Logg,them(Ny,Nx),logm(Ny,Nx),vm(Nv,7)
     &,Covar(2,2)
      LOGICAL debut,lotl(Ny,Nx)
      DATA debut/.TRUE./,lua/0/,lum1,lu0,lup1/10,11,12/
     &,them/Nxy*0./,logm/Nxy*0./
C
      IF(debut)THEN
         debut=.FALSE.
         OPEN(lum1,FILE='ptpgm10i',STATUS='OLD')
         OPEN(lu0 ,FILE='ptpgp00i',STATUS='OLD')
         OPEN(lup1,FILE='ptpgp10i',STATUS='OLD')
      ENDIF
C
         DO 1 ix=1,Nx
         DO 2 iy=1,Ny
              IF (Ic.EQ.0.) THEN
                 xm(ix)=0.02*(ix-10)
                 ym(iy)=0.02*(iy-6)
              ELSE IF (Ic.LT.0) THEN
                 xm(ix)=0.02*(ix-9)
                 ym(iy)=0.02*(iy-6)
              ELSE 
                 xm(ix)=0.02*(ix-10)
                 ym(iy)=0.02*(iy-11)
              ENDIF
    2    CONTINUE
    1    CONTINUE
C
      IF(Ic.EQ.0)THEN
         lul=lu0
      ELSE IF(Ic.LT.0)THEN
         lul=lum1
      ELSE
         lul=lup1
      END IF
      IF(lul.NE.lua)THEN
         lua=lul
         DO 3 iy=1,Ny
         DO 3 ix=1,Nx
    3       lotl(iy,ix)=.FALSE.
    4    READ(lul,*,END=9)xd,yd,xmet,Logg,Thef,xtef
          IF (Ic.EQ.0.) THEN
             ix=NINT(50.*xd+10.)
             iy=NINT(50.*yd+6.)
          ELSE IF (Ic.LT.0.) THEN  
             ix=NINT(50.*xd+9.)
             iy=NINT(50.*yd+6.) 
          ELSE
            ix=NINT(50.*xd+10.)
            iy=NINT(50.*yd+11.)
          ENDIF
         IF(ix.LT.1.OR.ix.GT.Nx)GO TO 4
         IF(iy.LT.1.OR.iy.GT.Ny)GO TO 4
         lotl(iy,ix)=.TRUE.                               
         them(iy,ix)=Thef
         logm(iy,ix)=Logg
         GO TO 4
    9    REWIND lul
      END IF
C
      CALL BISPLI(Yy,Xx,  ym,xm,them,lotl,Ny,Nx,   vm,Nv
     &           ,Thef,dtdy,dtdx,Extrap)
      fact=1.+2.*Extrap*Extrap
      Covar(1,1)=1E-4*(.4116*dtdx*dtdx+1.5064*dtdx*dtdy
     &           +2.5896*dtdy*dtdy)
      Dthef=fact*SQRT(Covar(1,1))
      CALL BISPLI(Yy,Xx,  ym,xm,logm,lotl,Ny,Nx,   vm,Nv
     &           ,Logg,dldy,dldx,Extrap)
      Covar(2,2)=1E-4*(.4116*dldx*dldx+1.5064*dldx*dldy
     &           +2.5896*dldy*dldy)
      Covar(1,2)=1E-4*(.4116*dtdx*dldx+.7532*dtdx*dldy+.7532*dldx*dtdy
     &+2.5896*dtdy*dldy)
      Covar(2,1)=Covar(1,2)
      Dlogg=fact*SQRT(Covar(2,2))
      END
C
C*******************************************************************************
C
      SUBROUTINE XYTL(Xx,Yy,Ic,  Thef,Logg,  Dthef,Dlogg,Covar,Extrap)
C
C        Calcul de Theta effective (Thef = 5040./tef1) et de Log(g)
C     a partir de X et Y, parametres de Cramer & Maeder (1979).
C        [Fe/H] est suppose connu, egal a 0 (solaire), -1 ou +1.
C        La routine donne une estimation des erreurs a partir de la
C     matrice de covariance sur les couleurs de Geneve. Ces erreurs
C     s'acroissent rapidement si on sort de la region ou la fonction
C     F: (Thef, Log(g) |----> (X, Y) n'est plus bijective.
C        Dans ladite region Extrap est nulle. Ailleurs l'extrapolation
C     presente des risques, du moins si Extrap depasse 0.5 a 1.
C
C     On donne :  Xx,    Yy    Parametres de Cramer & Maeder
C                 Ic           +1, 0, -1 [Fe/H]
C     ==========
C     On recoit:  Thef,  Logg  5040./Teff et log(g)
C     ========== Dthef, Dlogg  erreurs
C                Covar(2,2)    Matrice de covariance des erreurs
C                   Extrap     0. si on est dans la grille de definition
C                              nombre de mailles d'extrap. sinon.
C
      save debut
      PARAMETER (Nx=100,Ny=150,Nxy=Nx*Ny,Nv=150)
      REAL xm(Nx),ym(Ny),Logg,them(Ny,Nx),logm(Ny,Nx),vm(Nv,7)
     &,Covar(2,2)
      LOGICAL debut,lotl(Ny,Nx)
      DATA debut/.TRUE./,lua/0/,lum1,lu0,lup1/20,21,22/
     &,them/Nxy*0./,logm/Nxy*0./
C
      IF(debut)THEN
         debut=.FALSE.
         OPEN(lum1,FILE='xcycm10i',STATUS='OLD')
         OPEN(lu0 ,FILE='xcycp00i',STATUS='OLD')
         OPEN(lup1,FILE='xcycp10i',STATUS='OLD')
	END IF
C
         DO 1 ix=1,Nx         
         DO 2 iy=1,Ny
             IF (Ic.EQ.0.) THEN
                 xm(ix)=0.05*(ix-1)
                 ym(iy)=0.025*(iy-11)
              ELSE IF (Ic.LT.0) THEN
                 xm(ix)=0.05*ix
                 ym(iy)=0.025*(iy-9)
              ELSE 
                 xm(ix)=0.05*ix
                 ym(iy)=0.025*(iy-15)
             ENDIF
    2    CONTINUE
    1    CONTINUE      
C
      IF(Ic.EQ.0)THEN
         lul=lu0
      ELSE IF(Ic.LT.0)THEN
         lul=lum1
      ELSE
         lul=lup1
      END IF
      IF(lul.NE.lua)THEN
         lua=lul
         DO 3 iy=1,Ny
         DO 3 ix=1,Nx
    3       lotl(iy,ix)=.FALSE.
    4    READ(lul,*,END=9)xd,yd,xmet,Logg,Thef,xtef
          IF (Ic.EQ.0.) THEN
             ix=NINT(20.*xd+1.)
             iy=NINT(40.*yd+11.)
          ELSE IF (Ic.LT.0.) THEN  
             ix=NINT(20.*xd)
             iy=NINT(40.*yd+9.) 
          ELSE
            ix=NINT(20.*xd)
            iy=NINT(40.*yd+15.)
          ENDIF   
         IF(ix.LT.1.OR.ix.GT.Nx)GO TO 4
         IF(iy.LT.1.OR.iy.GT.Ny)GO TO 4
         lotl(iy,ix)=.TRUE.                               
         them(iy,ix)=Thef
         logm(iy,ix)=Logg
         GO TO 4
    9    REWIND lul
      END IF
C
      CALL BISPLI(Yy,Xx,  ym,xm,them,lotl,Ny,Nx,   vm,Nv
     &           ,Thef,dtdy,dtdx,Extrap)
      fact=1.+2.*Extrap*Extrap
      Covar(1,1)=1E-4*(.8586*dtdx*dtdx+.2590*dtdx*dtdy+1.9330*dtdy*dtdy)

      Dthef=fact*SQRT(Covar(1,1))
      CALL BISPLI(Yy,Xx,  ym,xm,logm,lotl,Ny,Nx,   vm,Nv
     &           ,Logg,dldy,dldx,Extrap)
      Covar(2,2)=1E-4*(1.7969*dldx*dldx+.0234*dldx*dldy+2.2319*dldy*
     &dldy)
      Covar(1,2)=1E-4*(1.7969*dtdx*dldx+.0117*dtdx*dldy+.0117*dldx*dtdy
     &+2.2319*dtdy*dldy)
      Covar(2,1)=Covar(1,2)
      Dlogg=fact*SQRT(Covar(2,2))
      END
