C       OPTIONS(SAVE,DREAL)
C       OPTIONS(SAVE,DREAL,CHECK)
        OPTIONS(SAVE,DREAL,OPTIMISE)
C               A CODE FOR Teff-LOG G DETERMINATION
C               FROM UBV,uvbybeta AND Geneva PHOTOMETRY AS WELL AS
C               FROM H5 AND H6 BALMER LINE PROFILES
C
c       authors: J. Budaj
c       contact: budaj@ta3.sk
c       published: http://www.ta3.sk/~budaj/software
c       special requirements: DBOS environement
c       non original routins:  intep (Press.et.al. 1992, Numerical recipes)
c       input: keyboard + phot.dat+north.dat+listsub+
c                 h5.dat+h6.dat+h5-1.dat+h6-1.dat+h5-2.dat+h6-2.dat
c       output:  kurucz.out
c
C               NMOM1 - NUMBER OF MODELS(+20 INTERP.MODELS),
C               NSYN  - NUMBER OF TABULATED POINTS IN BALM. PROFIL
C               NEZ2 - DIMENZATION PARAMETER
        PARAMETER (NMOM1=103,NSYN=26,NOBS=10000,NEZ2=20)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        CHARACTER*12 INN1,INN2
        DIMENSION CM(NOBS),B(NOBS),TEP(NMOM1),GRAV(NMOM1)
     P  ,H(NMOM1,NSYN),CC(NSYN),CAF(NEZ2),CDF(NEZ2),OMM(NMOM1),L(NMOM1)
     P  ,D(NMOM1),POM1(NSYN),POM2(NSYN)
C               INPUT DIALOG
        WRITE(*,*)'                     BEFORE RUN THIS CODE:'
        WRITE(*,*)'                     DOS-COMMAND <GRAPHICS>'
        WRITE(*,*)'                     ALLOWS YOU TO PRINT FIGURES'
        WRITE(*,*)'                     BY MEANS OF [PRTSC]'
        OPEN(9,FILE='KURUCZ.OUT',STATUS='UNKNOWN')
10      WRITE(*,*)
        WRITE(*,'(A,$)')' INPUT DATA FILE => '
        READ(5,'(A12)')INN1
        OPEN(7,FILE=INN1,STATUS='OLD')
        READ(7,*)
        READ(7,*)
        READ(7,*)UB1,BV1,VV1,EBV1,UB0,BV0
        READ(7,*)
        READ(7,*)SUB,SBY,SM1,SC1,BETA,WGAMA,EBY
        READ(7,*)
        READ(7,*)GU,GV,GB1,GB2,GV1,GG
        READ(7,*)
12      WRITE(*,*)
        WRITE(*,*)'---------------------'
        WRITE(*,*)'  PHOTOMETRY                         -1-  '
        WRITE(*,*)'  SPECTROSCOPY                       -2-  '
        WRITE(*,*)'  END                                -3-  '
        WRITE(*,'(A,$)')' ---------------------                   => '
        READ(5,*)IAD
        IF (IAD.EQ.2)GOTO 14
        IF (IAD.EQ.3)GOTO 400
C------------------------------PHOTOMETRY-------------------------------
        IF(ABS(UB1-1.E6).GT.1.E-10.OR.ABS(BV1-1.E6).GT.1.E-10.OR.
     P  ABS(UB0-1.E6).GT.1.E-10.OR.ABS(BV0-1.E6).GT.1.E-10.OR.
     P  ABS(SUB-1.E6).GT.1.E-10.OR.ABS(SBY-1.E6).GT.1.E-10.OR.
     P  ABS(SM1-1.E6).GT.1.E-10.OR.ABS(SC1-1.E6).GT.1.E-10.OR.
     P  ABS(GU-1.E6).GT.1.E-10.OR.ABS(GV-1.E6).GT.1.E-10.OR.
     P  ABS(GB1-1.E6).GT.1.E-10.OR.ABS(GB2-1.E6).GT.1.E-10.OR.
     P  ABS(GV1-1.E6).GT.1.E-10.OR.ABS(GG-1.E6).GT.1.E-10)THEN
C     P  ABS(VV1-1.E6).GT.1.E-10.OR.
C     P  ABS(BETA-1.E6).GT.1.E-10.OR.ABS(WGAMA-1.E6).GT.1.E-10.OR.
          CALL PHOT(UB1,BV1,VV1,EBV1,UB0,BV0,SUB,SBY,SC1,SM1,BETA,
     P    WGAMA,EBY,GU,GV,GB1,GB2,GV1,GG)
          GOTO 12
        ENDIF
C------------------------------SPECTROSCOPY-----------------------------
C        WRITE(*,*)'                 TABULATED DELTA LAMBDA'
C        WRITE(*,*)'                 0,0.02,0.04,0.06,0.08,0.1,0.15,0.2,'
C        WRITE(*,*)'                 0.3,0.4,0.5,0.6,0.8,1,1.2,1.4,1.6,'
C        WRITE(*,*)'                 1.8,2,2.2,2.4,2.8,3.2,3.6,4,5 NM'
C        WRITE(*,*)
C               READS FROM THE INPUT FILE
14      I=1
        READ(7,*)IEPER,CENTER,UNIT
        READ(7,*)
20      READ(7,*,END=30)CM(I),B(I)
        CM(I)=(CM(I)-CENTER)*UNIT
        B(I)=B(I)*1.0D3
        I=I+1
        GOTO 20
30      NBOD=I-1
        WRITE(*,*)
        WRITE(*,'(30X,A15,I5)')' No. OF POINTS ',NBOD
        CLOSE(7)
        WRITE(*,*)
        WRITE(*,*)'SCATTERED LIGHT IN %,'
        WRITE(*,'(A,$)')' IF NO PUT         -0.-                    => '
        READ(5,*)AM
C               OBSERVED PROFIL WILL BE CORRECTED TO THE SCATTERED LIGHT
        IF (AM.GT.0.)THEN
          DO 60 I=1,NBOD
            B(I)=(B(I)-AM/100.*1000.)/(1.-AM/100.)
60        CONTINUE
        ENDIF
15      WRITE(*,*)
        WRITE(*,*)'TO NORMALIZE SYNTHETIC KURUCZ PROFILES GIVE AT WHICH'
        WRITE(*,*)'DELTA LAMBDA [NM] THEY SHOULD MERGE INTO COMTINUUM'
        WRITE(*,'(A,$)')' IF NO PUT         -0.-                    => '
        READ(5,*)AL
        WRITE(*,*)
        WRITE(*,*)'ROTATION V*SIN I [KM/S]'
        WRITE(*,'(A,$)')' IF NO PUT         -0.-                    => '
        READ(5,*)BV
        IF (BV.GT.0.) THEN
          WRITE(*,*)
          WRITE(*,*)'DOUBLE LETTER CODE OF EXPECTED SP. TYPE "XY"'
          WRITE(*,*)'(X: O=0,B=1,A=2,F=3,G=4,K=5, Y=No. OF SP.SUBCLASS'
          WRITE(*,'(A,$)')' E.G. A4=24)                             => '
          READ(5,*)IBA
        ENDIF
        WRITE(*,*)
        WRITE(*,*)'SOLAR ABUNDANCES                     -0- '
        WRITE(*,*)'METALS 10-TIMES LOWERED              -1- '
        WRITE(*,'(A,$)')' METALS 100-TIMES LOWERED             -2-  => '
        READ(5,*)IEPP
        WRITE(*,*)
        WRITE(*,*)'TEMPERATURE RANGE OF INTEREST T1<T2 [K] AND'
        WRITE(*,'(A,$)')' LOG G RANGE G1<G2 [CGS] <T1,T2,G1,G2>     => '
        READ(5,*)CA,CB,CD,CE
        WRITE(*,*)
        WRITE(*,*)'TO INTERPOLATE BETWEEN MODELS:  NO  -0- '
        WRITE(*,*)'                                YES -1- '
        WRITE(*,*)'TO INTERPOLATE+LIST OF MODELS       -3- '
        WRITE(*,*)'(YOU WILL NEED IT IF YOU WANT TO'
        WRITE(*,'(A,$)')' CONSIDER NON GRID MODELS)                 => '
        READ(5,*)INTRP
        IF (INTRP.EQ.3)THEN
          WRITE(*,*)
          WRITE(*,*)'-  S W I T C H   O N   T H E   P R I N T E R -'
          WRITE(*,*)'           AND PRESS [ENTER]'
          PAUSE
          CALL TLAC(IEPER,IEPP)
          INTRP=1
        ENDIF
C               NEPREPIS SI IVOL1, ABY NEINTERP. INY POCET MODELOV
        IVOL1=0
        IF (INTRP.EQ.1)THEN
          WRITE(*,*)
          WRITE(*,'(15X,A,A)')'PARAMETERS OF INTERP. MODELS (T,LOG G)',
     P    ' MUST LAY BETWEEN'
          WRITE(*,'(15X,A,A)')'PARAMETERS OF TWO COUPLES OF KURUCZ',
     P    ' DISCRETE MODELS'
          WRITE(*,'(15X,A)')'( EACH COUPLE MUST HAVE THE SAME LOG G)'
          WRITE(*,*)
          WRITE(*,*)'SEQUENCE No. OF KURUCZ MODELS WITH LARGER TEMPER.'
          WRITE(*,'(A,$)')' FIRST THAT WITH LOWER LOG G  <-,->      => '
          READ(5,*)I21,I2
          WRITE(*,*)'SEQUENCE No. OF KURUCZ MODELS WITH LOWER TEMPER.'
          WRITE(*,'(A,$)')' FIRST THAT WITH LOWER LOG G  <-,->      => '
          READ(5,*)I1,I11
          WRITE(*,*)
          WRITE(*,*)'HOW MANY MODELS DO YOU WANT TO INTERP. WITHIN'
          WRITE(*,'(A,$)')' ABOVE BOUNDARIES?       (MAXIM. -20-)   => '
          READ(5,*)IVOL1
          DO 17 I=1,IVOL1
          WRITE(*,'(A,$)')' INPUT <T,LOG G> OF INTERP. MODELS       => '
            READ(5,*)CAF(I),CDF(I)
17        CONTINUE
        ENDIF
        WRITE(*,*)
        WRITE(*,'(15X,A,A)')' NOTE: OBSERVED PROFIL IS CORRECTED FOR',
     P  ' SCATTERED LIGHT'
        WRITE(*,'(15X,A)')' ALL SYNTHETIC PROFILES ARE NORMALIZED'
        WRITE(*,'(15X,A,A)')' ONLY SYNTH. PROF. WITHIN GIVEN',
     P  ' INTERVAL ARE CORRECTED FOR ROTATION '
C       WRITE(*,*)'                       ZOZNAM MODELOV BUDE NA TERMINA
C    PLI A V SUBORE -LIST.SUB-'
        WRITE(*,*)
        WRITE(*,*)'---------------------'
        WRITE(*,*)'  REPEATE DIALOGUE FROM BEGINNING     -1-  '
        WRITE(*,*)'  REPEATE IT WITHOUT LOADING DATA     -2-  '
        WRITE(*,*)'  END                                 -3-  '
        WRITE(*,*)'  GO ON                               -4-  '
        WRITE(*,'(A,$)')' ---------------------                     => '
        READ(5,*)IAD
        IF (IAD.EQ.1)GOTO 10
        IF (IAD.EQ.2)GOTO 15
        IF (IAD.EQ.3)GOTO 400
C-----------------------------------------------------------------------
C               THIS READS SYNTHETIC PROFILES
        DATA CC/0.,0.02,0.04,0.06,0.08,0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.8
     P  ,1.,1.2,1.4,1.6,1.8,2.,2.2,2.4,2.8,3.2,3.6,4.,5./
        IF (IEPER.EQ.5)THEN
          IF (IEPP.EQ.0)OPEN(7,FILE='H5.DAT',STATUS='OLD')
          IF (IEPP.EQ.1)OPEN(7,FILE='H5-1.DAT',STATUS='OLD')
          IF (IEPP.EQ.2)OPEN(7,FILE='H5-2.DAT',STATUS='OLD')
        ELSE
          IF (IEPP.EQ.0)OPEN(7,FILE='H6.DAT',STATUS='OLD')
          IF (IEPP.EQ.1)OPEN(7,FILE='H6-1.DAT',STATUS='OLD')
          IF (IEPP.EQ.2)OPEN(7,FILE='H6-2.DAT',STATUS='OLD')
        ENDIF
        DO 70 I=1,NMOM1
          READ(7,*,END=75)TEP(I),GRAV(I),(H(I,J),J=1,NSYN)
70      CONTINUE
75      NMOD=I-1
        NMOM=NMOD+IVOL1
        CLOSE(7)
        DO 77 I=1,IVOL1
          TEP(NMOD+I)=CAF(I)
          GRAV(NMOD+I)=CDF(I)
77      CONTINUE
C------------INTERPOLATION OF ADDITIONAL SYNTH. KURUCZ MODELS-----------
        IF (INTRP.EQ.1)THEN
          DO 170 I=1,IVOL1
            DO 160 J=1,NSYN
C               INTERP. TO CORRECT TEMPER. AT HIGHER LOG G
              POM1(J)=(H(I2,J)-H(I11,J))*(CAF(I)-TEP(I11))
     P        /(TEP(I2)-TEP(I11))+H(I11,J)
C               INTERP. TO CORRECT TEMPER. AT LOWER LOG G
              POM2(J)=(H(I21,J)-H(I1,J))*(CAF(I)-TEP(I1))
     P        /(TEP(I21)-TEP(I1))+H(I1,J)
C               INTERP. TO CORRECT LOG G
              H(NMOD+I,J)=(POM1(J)-POM2(J))*(CDF(I)-GRAV(I1))
     P        /(GRAV(I2)-GRAV(I1))+POM2(J)
160         CONTINUE
170       CONTINUE
        ENDIF
C---------------ROTATIONAL BROADENING OF MODELS WITHIN GIVEN INTERVAL---
        IF (BV.GT.0.)THEN
          CALL STEM(IBA,IEPER,BV,UKOEF,DLL)
          DO 80 I=1,NMOM
            IF (TEP(I).GE.CA.AND.TEP(I).LE.CB.AND.GRAV(I).GE.CD.AND.
     P      GRAV(I).LE.CE)THEN
              DO 78 J=1,NSYN
                POM1(J)=H(I,J)
78            CONTINUE
              WRITE(*,'(F6.0,F6.2,A)')TEP(I),GRAV(I),
     P        '  IS BEING CORRECTED FOR ROTATION'
              CALL AKONV(NSYN,DLL,UKOEF,CC,POM1)
              DO 79 J=1,NSYN
                H(I,J)=POM1(J)
79            CONTINUE
            ENDIF
80        CONTINUE
        ENDIF
C---------------NORMALIZATION OF ALL SYNTH .PROFILES--------------------
        IF (AL.GT.0.)THEN
          DO 90 I=1,NMOM
            CALL NORM(I,NMOM1,NSYN,AL,CC,H)
90        CONTINUE
        ENDIF
C---------------COMPARING OF OBSERV. AND SYNTH. PROFILES----------------
        IKK=0
        DO 210 I=1,NMOM
          D(I)=0.
210     CONTINUE
        DO 250 I=1,NMOM
          IF (TEP(I).GE.CA.AND.TEP(I).LE.CB.AND.GRAV(I).GE.CD.AND.
     P    GRAV(I).LE.CE) THEN
            IKK=IKK+1
            IKK1=0
            DO 230 IL1=1,NSYN
              KRA=NSYN-IL1+1
              IF(-CC(KRA).GT.CM(1).AND.-CC(KRA).LT.CM(NBOD))THEN
                CALL INTEP(-CC(KRA),YDEL,CM,B,NBOD,IER)
                D(I)=(H(I,KRA)-YDEL)*(H(I,KRA)-YDEL)+D(I)
                IKK1=IKK1+1
              ENDIF
              IF(CC(IL1).GT.CM(1).AND.CC(IL1).LT.CM(NBOD))THEN
                CALL INTEP(CC(IL1),YDEL,CM,B,NBOD,IER)
                D(I)=(H(I,IL1)-YDEL)*(H(I,IL1)-YDEL)+D(I)
                IKK1=IKK1+1
              ENDIF
230         CONTINUE
            D(I)=SQRT(D(I))/DBLE(IKK1)
          ENDIF
250     CONTINUE
C-----------------------------------------------------------------------
        WRITE(*,*)
        WRITE(*,'(A,$)')' HOW MANY BEST MODELS AM I CHOOSE          => '
        READ(5,*)IVOL
        WRITE(*,*)
        IF (IKK.LE.IVOL) IVOL=IKK
        CALL ZORAD(D,OMM,L,NMOM,IVOL,NMOM1)
        WRITE(*,*)
        WRITE(9,*)
        WRITE(9,'(15X,A21,I1)')'   FROM BALMER LINE H',IEPER
        WRITE(9,*)
        DO 300 J=1,IVOL
          WRITE(*,310)L(J),TEP(L(J)),GRAV(L(J)),OMM(J)
          WRITE(9,310)L(J),TEP(L(J)),GRAV(L(J)),OMM(J)
300     CONTINUE
310     FORMAT(' No:',I2,'  T=',F6.0,'  LOG G=',F4.2,
     P  '  DEV.:',E11.4)
320     WRITE(*,*)
        WRITE(*,*)'------------------'
        WRITE(*,*)'VIEW                              -1- '
        WRITE(*,*)'PROCEED WITH NEW DATA             -2- '
        WRITE(*,*)'PROCEED WITH OLD DATA             -3- '
        WRITE(*,*)'SAVE DATA FROM VIEW TO DISC       -4- '
        WRITE(*,*)'END                               -5- '
        WRITE(*,'(A,$)')' ------------------                        => '
        READ(5,*)IAD
        IF (IAD.EQ.1)THEN
          WRITE(*,*)
          WRITE(*,*)'                        NOTE: [ENTER] OR [ESC] '
          WRITE(*,*)'                              TO ESCAPE GRAPHICS'
          WRITE(*,'(A,$)')' INPUT No. OF MODEL                      => '
          READ(5,*)IVOL
          CALL GRAF(CC,CM,B,H,TEP,GRAV,NSYN,NOBS,NMOM1,NBOD,IVOL
     P    ,AM,AL,BV,IEPER,IEPP,CA,CB,CD,CE)
          GOTO 320
        ENDIF
        IF (IAD.EQ.2) GOTO 10
        IF (IAD.EQ.3) GOTO 15
        IF (IAD.EQ.4)THEN
          WRITE(*,*)
          WRITE(*,'(A,$)')' No. OF MODEL                              =>
     P '
          READ(5,*)IVOL
          WRITE(*,'(A,$)')' FILE NAME WITH OBSERVED PROFIL            =>
     P '
          READ(5,'(A12)')INN1
          WRITE(*,'(A,$)')' FILE NAME WITH SYNTHETIC PROFIL           =>
     P '
          READ(5,'(A12)')INN2
          CALL NAHRAJ(CC,CM,B,H,IVOL,NSYN,NOBS,NBOD,NMOM1,INN1
     P    ,INN2)
          GOTO 320
        ENDIF
        IF (IAD.EQ.5) GOTO 400
        GOTO 10
400     CLOSE(9)
        STOP
        END
C-----------------------------------------------------------------------
C              PODPROGRAMY
C-----------------------------------------------------------------------
        SUBROUTINE NORM(I,NMOM1,NSYN,AL,CC,H)
C       IT WILL NORMALIZE SYNTHETIC PROFILES
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION CC(NSYN),H(NMOM1,NSYN)
        IF (AL.GT.CC(NSYN))THEN
          AP=(1000.-H(I,NSYN))*(AL-CC(NSYN))/(10.-CC(NSYN))+H(I,NSYN)
          GOTO 110
        ENDIF
        DO 100 IL=1,NSYN
          IF(CC(IL)-AL) 100,80,90
80          AP=H(I,IL)
            GOTO 110
90          AP=(H(I,IL)-H(I,IL-1))*(AL-CC(IL-1))/(CC(IL)-CC(IL-1))
     P      +H(I,IL-1)
            GOTO 110
100     CONTINUE
110     IF(ABS(AP-0.).LT.1.E-8) AP=1000.
        DO 120 J=1,NSYN
          H(I,J)=H(I,J)/AP*1000.
          IF (H(I,J).GT.1000.) H(I,J)=1000.
120     CONTINUE
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE TLAC(IEPER,IEPP)
C       IT WILL PRINT THE CURRENT LIST OF MODELS AVAILABLE
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        IF (IEPER.EQ.5)THEN
          IF (IEPP.EQ.0)OPEN(4,FILE='H5.DAT',STATUS='OLD')
          IF (IEPP.EQ.1)OPEN(4,FILE='H5-1.DAT',STATUS='OLD')
          IF (IEPP.EQ.2)OPEN(4,FILE='H5-2.DAT',STATUS='OLD')
        ELSE
          IF (IEPP.EQ.0)OPEN(4,FILE='H6.DAT',STATUS='OLD')
          IF (IEPP.EQ.1)OPEN(4,FILE='H6-1.DAT',STATUS='OLD')
          IF (IEPP.EQ.2)OPEN(4,FILE='H6-2.DAT',STATUS='OLD')
        ENDIF
        OPEN(3,FILE='LIST.SUB',STATUS='UNKNOWN')
        OPEN(6,FILE='PRN')
        WRITE(*,5)IEPER,10.0**FLOAT(-IEPP)
        WRITE(3,5)IEPER,10.0**FLOAT(-IEPP)
        WRITE(6,5)IEPER,10.0**FLOAT(-IEPP)
5       FORMAT(' MODELS FOR H',I1,' AND METAL.ABUND./SOLAR =',F4.2)
        ZMAZ=0.
10      READ(4,*,END=20)ZMAZ1,ZMAZ2
        ZMAZ=ZMAZ+1.
        WRITE(*,30)ZMAZ,ZMAZ1,ZMAZ2
        WRITE(3,30)ZMAZ,ZMAZ1,ZMAZ2
        WRITE(6,30)ZMAZ,ZMAZ1,ZMAZ2
        GOTO 10
20      CONTINUE
        WRITE(*,*)'IF YOU DID NOT CATCH IT IT WILL APPEAR IN FILE:'
        WRITE(*,*)' "LIST.SUB" AND ON YOUR PRINTER AS WELL'
        CLOSE(6)
        CLOSE(4)
        CLOSE(3)
30      FORMAT(F4.0,F7.0,F5.2)
        RETURN
        END
C-----------------------------------------------------------------------
        FUNCTION ROT(DL,DLL,UKOEF)
C               ROTATION PROFIL IN [NM]
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        ROT=(2.*(1.-UKOEF)*SQRT(1.-DL*DL/DLL/DLL)+0.5*3.1415*UKOEF
     P  *(1.-DL*DL/DLL/DLL))/(3.1415*DLL*(1.-UKOEF/3.))
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE STEM(IBA,IEPER,BV,UKOEF,DLL)
C       OUTPUT: UKOEF - LIMB DARK. COEFF.
C               DLL[NM] - DOPPLER WIDTH
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION G5(50),G6(50),G(50)
        DO 5 I=1,50
          IF (IEPER.EQ.5) G(I)=G5(I)
          IF (IEPER.EQ.6) G(I)=G6(I)
5       CONTINUE
        DO 10 I=0,10
          IF (ABS(G(IBA-I)-0.).LT.1.E-20)GOTO 10
          IBE=IBA-I
          GOTO 20
10      CONTINUE
20      DO 30 I=0,10
          IF (ABS(G(IBA+I)-0.).LT.1.E-20)GOTO 30
          IBF=IBA+I
          GOTO 40
30      CONTINUE
40      IF (IBA.EQ.IBE)THEN
          UKOEF=G(IBA)
        ELSE
          UKOEF=FLOAT(IBA-IBE)/FLOAT(IBF-IBE)*(G(IBF)-G(IBE))+G(IBE)
        ENDIF
        WRITE(*,*)
        WRITE(*,*)'LIMB DARKENING COEFFICIENT',UKOEF
        WRITE(*,*)
        IF (IEPER.EQ.5)THEN
          DLL=BV*434.0/299792.
        ELSE
          DLL=BV*410.2/299792.
        ENDIF
C               DEFINED BY EQUATION:  I=I0[1-UKOEF+UKOEF*COS(THETA)]
C       FROM GRYGAR ...
C               LIMB DARK. COEFF. -U, FOR H5------------
        DATA G5/0.,0.,0.,0.,0.,0.,0.,.35,.36,.38,0.,.4,0.,.44,0.,.48
     P  ,0.,.52,0.,.56,0.,0.,.69,0.,0.,.75,0.,0.,0.,.78,0.,0.,0.,0.
     P  ,0.,0.,.8,0.,0.,0.,0.,.77,0.,0.,.84,0.,0.,0.,0.,.92/
C               LIMB DARK. COEFF. -U, FOR H6------------
        DATA G6/0.,0.,0.,0.,0.,0.,0.,.36,.37,.39,0.,.42,0.,.46,0.,.5
     P  ,0.,.54,0.,.58,0.,0.,.74,0.,0.,.8,0.,0.,0.,.82,0.,0.,0.,0.
     P  ,0.,0.,.83,0.,0.,0.,0.,.83,0.,0.,.89,0.,0.,0.,0.,.97/
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE AKONV(NSYN,DLL,UKOEF,CC,POM1)
C       KONVOLUTION WITH ROTATION PROFIL
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION CC(NSYN),POM2(30),POM1(NSYN)
        DO 30 J=1,NSYN
          YDL0=0.
          POM2(J)=0.
          ROT0=0.
          DO 20 K=1,49
            DL=DLL*(K-25)/25.
            IF (CC(J)+DL.LE.0.)THEN
              DL1=-CC(J)-DL
            ELSE
              IF(CC(J)+DL.GE.CC(NSYN))THEN
                IF (CC(J)+DL.GE.10.)WRITE(*,*)'TOO SPEED ROTATION'
                YDL=(CC(J)+DL-CC(NSYN))/(10.-CC(NSYN))*(1.-POM1(J))
     P          +POM1(J)
                POM2(J)=ROT(DL,DLL,UKOEF)*(YDL+YDL0)/2.+POM2(J)
                YDL0=YDL
                GOTO 20
              ELSE
                DL1=CC(J)+DL
              ENDIF
            ENDIF
            CALL INTEP(DL1,YDL,CC,POM1,NSYN,IER)
            ROT1=ROT(DL,DLL,UKOEF)
C       WRITE(*,*)'ROT YDL YDL0 POM',ROT(DL,DLL,UKOEF),YDL,YDL0,POM2(J)
            POM2(J)=0.5*(ROT1*YDL+ROT0*YDL0)+POM2(J)
            YDL0=YDL
            ROT0=ROT1
20        CONTINUE
          POM2(J)=0.5*YDL*ROT1+POM2(J)
          POM2(J)=POM2(J)*2.*DLL/50.
30      CONTINUE
        DO 40 J=1,NSYN
          POM1(J)=POM2(J)
40      CONTINUE
        RETURN
        END
C-----------------------------------------------------------------------
C       SUBROUTINE   I N T E P  : SPLINE INTERPOLATION SCHEME BASED
C                                 ON HERMITE POLYNOMIALS
C        INPUT : XP  - THE CHOSEN ARGUMENT VALUE
C                X   - THE VECTOR OF INDEPENDENT VALUES
C                F   - THE VECTOR OF FUNCTION OR DEPENDENT VALUES
C                N   - THE NUMBER OF POINTS IN THE (X,F) VECTORS
C       OUTPUT : P   - THE RESULTANT INTERPOLATED VALUE
C                IER - THE RESULTANT ERROR PARAMETER
        SUBROUTINE INTEP(XP,P,X,F,N,IER)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        REAL LP1,LP2,L1,L2
        DIMENSION F(N),X(N)
        IER=1
        IO=1
        IUP=0
        IF(X(2).LT.X(1)) IUP=1
        N1=N-1
        IF((XP.GE.X(N).AND.IUP.EQ.0).OR.(XP.LE.X(N).AND.IUP.EQ.1)) THEN
         P=F(N)
         IER=2
         RETURN
        ELSE IF((XP.LE.X(1).AND.IUP.EQ.0).OR.
     *   (XP.GE.X(1).AND.IUP.EQ.1)) THEN
         P=F(1)
         IER=2
         RETURN
        END IF
        DO 1 I=IO,N
         IF(XP.LT.X(I).AND.IUP.EQ.0) GOTO 2
         IF(XP.GT.X(I).AND.IUP.EQ.1) GOTO 2
1       CONTINUE
        P=F(N)
        IER=2
        RETURN
2       I=I-1
        IF(I.EQ.IO-1) GOTO 4
        IO=I+1
        LP1=1./(X(I)-X(I+1))
        LP2=1./(X(I+1)-X(I))
        IF(I.EQ.1) FP1=(F(2)-F(1))/(X(2)-X(1))
        IF(I.EQ.1) GOTO 3
        FP1=(F(I+1)-F(I-1))/(X(I+1)-X(I-1))
3       IF(I.GE.N1) FP2=(F(N)-F(N-1))/(X(N)-X(N-1))
        IF(I.GE.N1) GOTO 4
        FP2=(F(I+2)-F(I))/(X(I+2)-X(I))
4       XPI1=XP-X(I+1)
        XPI=XP-X(I)
        L1=XPI1*LP1
        L2=XPI*LP2
        P=F(I)*(1.-2.*LP1*XPI)*L1*L1+F(I+1)*(1.-2.*LP2*XPI1)
     *  *L2*L2+FP2*XPI1*L2*L2+FP1*XPI*L1*L1
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE GRAF(CC,CM,B,H,TEP,GRAV,NSYN,NOBS,NMOM1,NBOD
     P  ,IVOL,AM,AL,BV,IEPER,IEPP,CA,CB,CD,CE)
C       GRAPHICS
        PARAMETER (NOBS1=10000,NSYN1=26)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION CC(NSYN),CM(NOBS),B(NOBS),H(NMOM1,NSYN),TEP(NMOM1)
     P  ,GRAV(NMOM1),IC(NSYN1),ICM(NOBS1),IB(NOBS1),IH(NSYN1)
        SCAL1=640./6.
        SCAL2=480./1200.
        IP1=40
        IP2=480-30
        DO 10 J=1,NSYN
          IC(J)=INT(CC(J)*SCAL1)+IP1
          IH(J)=-INT(H(IVOL,J)*SCAL2)+IP2
10      CONTINUE
        DO 20 J=1,NBOD
          ICM(J)=INT(DABS(CM(J))*SCAL1)+IP1
          IB(J)=-INT(B(J)*SCAL2)+IP2
20      CONTINUE
        CALL VGA@
C               X-AXIS
        CALL DRAW_LINE@(IP1,IP2,640,IP2,2)
        DO 30 J=1,5
          CALL DRAW_LINE@(INT(J*SCAL1)+IP1,IP2,INT(J*SCAL1)+IP1,IP2-5,2)
30      CONTINUE
C               Y-AXIS
        CALL DRAW_LINE@(IP1,IP2,IP1,0,2)
        DO 40 J=1,5
          CALL DRAW_LINE@(IP1,-INT(200*J*SCAL2)+IP2,IP1+3
     P    ,-INT(200*J*SCAL2)+IP2,2)
40      CONTINUE
C               OBSERVED (CORRECTED TO SCATTERED LIGHT) PROFIL
        DO 50 J=1,NBOD
          CALL CROSS(ICM(J),IB(J),2,5)
50      CONTINUE
C               SYNTHETIC (CORRECTED ...) PROFIL
        CALL POLYLINE@(IC,IH,NSYN,1)
C               DESCRIPTION OF THE GRAPH
        CALL SET_CURSOR_POS@(47,11)
        WRITE(*,60)TEP(IVOL)
60      FORMAT('T=',F7.0,' K')
        CALL SET_CURSOR_POS@(47,13)
        WRITE(*,70)GRAV(IVOL)
70      FORMAT('LOG G=',F5.2,' [CGS]')
        CALL SET_CURSOR_POS@(47,15)
        WRITE(*,80)IEPER
80      FORMAT('LINE: H',I1)
        CALL SET_CURSOR_POS@(47,17)
        WRITE(*,90)10.0**FLOAT(-IEPP)
90      FORMAT('MET.ABUND./SOLAR=',F4.2)
        BVDEL=0.
        IF (TEP(IVOL).GE.CA.AND.TEP(IVOL).LE.CB.AND.GRAV(IVOL).GE.CD
     P  .AND.GRAV(IVOL).LE.CE)BVDEL=BV
        CALL SET_CURSOR_POS@(47,19)
        WRITE(*,100)BVDEL
100     FORMAT('ROT.VEL.=',F4.0,' KM/S')
        CALL SET_CURSOR_POS@(47,21)
        WRITE(*,110)AM
110     FORMAT('SCAT.LIGHT:',F4.1,' %')
        CALL SET_CURSOR_POS@(47,23)
        WRITE(*,120)AL
120     FORMAT('NORM.:',F4.1,' NM')
130     CALL GET_KEY@(IVOL)
        IF (IVOL.NE.13.AND.IVOL.NE.27)GOTO 130
        CALL TEXT_MODE@
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE CROSS(IX,IY,ILENG,ICOLL)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        CALL DRAW_LINE@(IX-ILENG,IY,IX+ILENG,IY,ICOLL)
        CALL DRAW_LINE@(IX,IY-ILENG,IX,IY+ILENG,ICOLL)
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE NAHRAJ(CC,CM,B,H,IVOL,NSYN,NOBS,NBOD,NMOM1,INN1,INN2)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        CHARACTER*12 INN1,INN2
        DIMENSION CC(NSYN),CM(NOBS),B(NOBS),H(NMOM1,NSYN)
        OPEN(6,FILE=INN1,STATUS='UNKNOWN')
        OPEN(7,FILE=INN2,STATUS='UNKNOWN')
        DO 10 I=1,NBOD
          WRITE(6,30)CM(I),B(I)/1000.
10      CONTINUE
        DO 20 I=1,NSYN
          WRITE(7,40)CC(I),H(IVOL,I)/1000.
20      CONTINUE
30      FORMAT(F6.3,2F7.3)
40      FORMAT(F6.3,F7.3)
        CLOSE(6)
        CLOSE(7)
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE ZORAD(D,OMM,L,NMOM,IVOL,NMOM1)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION D(NMOM1),OMM(NMOM1),L(NMOM1)
        DO 260 I=1,IVOL
          OMM(I)=1.E20
          L(I)=0
260     CONTINUE
        DO 300 J=1,IVOL
          DO 290 I=1,NMOM
            DO 280 IL=1,IVOL
              IF (I.EQ.L(IL).OR.ABS(D(I)-0.).LT.1.E-20) GOTO 290
280         CONTINUE
            IF (D(I).LT.OMM(J)) THEN
              OMM(J)=D(I)
              L(J)=I
            ENDIF
290       CONTINUE
300     CONTINUE
        RETURN
        END
C------------------------------PHOTOMETRY-------------------------------
        SUBROUTINE PHOT(UB1,BV1,VV1,EBV1,UB0,BV0,SUB,SBY,SC1,SM1,BETA,
     P  WGAMA,EBY1,GU,GV,GB1,GB2,GV1,GG)
C               NMOM1 - POCET MODELOV DIMENZACIA
C               NMOM  - POCET MODELOV
C               UB1,BV1,V1=OBSERVED, UB0,BV0,V0=DEREDDENED
        PARAMETER (NMOM1=500)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION TEP(NMOM1),GRAV(NMOM1),AB(NMOM1)
     P  ,UB(NMOM1),BV(NMOM1),C1(NMOM1),C2(NMOM1),BY(NMOM1)
     P  ,D(NMOM1),OMM(NMOM1),L(NMOM1)
        OPEN(8,FILE='PHOT.DAT',STATUS='OLD')
        READ(8,*)
        I=1
10      READ(8,*,END=30)TEP(I),GRAV(I),AB(I),UB(I),BV(I),C1(I),BY(I)
        I=I+1
        GOTO 10
30      NMOM=I-1
        CLOSE(8)
C---------------U-B-V-----
        IF(ABS(BV1-1.E6).GT.1.E-10.OR.ABS(UB1-1.E6).GT.1.E-10.OR.
     P  ABS(BV0-1.E6).GT.1.E-10.OR.ABS(UB0-1.E6).GT.1.E-10)THEN
          WRITE(*,*)
          WRITE(9,*)
          WRITE(*,'(10X,A10)')'  U B V : '
          WRITE(9,'(10X,A10)')'  U B V : '
          UB2=1.E6
          BV2=1.E6
          EBV2=1.E6
          IF(ABS(EBV1-1.E6).LT.1.E-10.AND.ABS(UB0-1.E6).LT.1.E-10.AND.
     P    ABS(BV0-1.E6).LT.1.E-10.AND.ABS(UB1-1.E6).GT.1.E-10.AND.
     P    ABS(BV1-1.E6).GT.1.E-10) THEN
            CALL HEINTZE(UB1,BV1,EBV2,UB2,BV2)
          ENDIF
        ENDIF
C
        IF(ABS(BV0-1.E6).GT.1.E-10.AND.ABS(UB0-1.E6).GT.1.E-10)THEN
          CALL HYLAND(UB0,BV0)
          CALL SCHILD(UB0,BV0)
          PAUSE
          GOTO 40
        ENDIF
        IF(ABS(BV1-1.E6).GT.1.E-10.AND.ABS(UB1-1.E6).GT.1.E-10)THEN
          CALL HYLAND(UB1,BV1)
          CALL SCHILD(UB1,BV1)
          PAUSE
        ENDIF
C
40      IF(ABS(BV1-1.E6).GT.1.E-10.OR.ABS(BV0-1.E6).GT.1.E-10)THEN
          IF(ABS(BV0-1.E6).GT.1.E-10) THEN
            CALL UNDERHILL(BV0)
          ELSE
            IF(ABS(BV2-1.E6).GT.1.E-10) CALL UNDERHILL(BV2)
            IF(ABS(BV1-1.E6).GT.1.E-10) CALL UNDERHILL(BV1)
          ENDIF
          PAUSE
        ENDIF
C
        IF(ABS(BV1-1.E6).GT.1.E-10.OR.ABS(BV0-1.E6).GT.1.E-10)THEN
          IF(ABS(BV0-1.E6).GT.1.E-10) THEN
            CALL NAPI3(BV0)
          ELSE
            IF(ABS(BV2-1.E6).GT.1.E-10) CALL NAPI3(BV2)
            IF(ABS(BV1-1.E6).GT.1.E-10) CALL NAPI3(BV1)
          ENDIF
          PAUSE
        ENDIF
C
        IF(ABS(UB1-1.E6).GT.1.E-10.OR.ABS(UB0-1.E6).GT.1.E-10)THEN
          IF(ABS(UB0-1.E6).GT.1.E-10) THEN
            CALL HHSA(UB0)
          ELSE
            IF(ABS(UB2-1.E6).GT.1.E-10) CALL HHSA(UB2)
            IF(ABS(UB1-1.E6).GT.1.E-10) CALL HHSA(UB1)
          ENDIF
          PAUSE
        ENDIF
C
        IF(ABS(BV1-1.E6).GT.1.E-10)BV3=BV1
        IF(ABS(BV2-1.E6).GT.1.E-10)BV3=BV2
        IF(ABS(BV0-1.E6).GT.1.E-10)BV3=BV0
        IF(ABS(BV1-1.E6).GT.1.E-10.OR.ABS(BV0-1.E6).GT.1.E-10)THEN
          DO 250 I=1,NMOM
            D(I)=ABS(BV(I)-BV3)
            IF(ABS(D(I)-0.).LT.1.E-4) D(I)=0.0001
250       CONTINUE
          WRITE(*,*)
          WRITE(9,*)
          WRITE(*,'(5X,A20,F6.3)')' considering (B-V)o=',BV3
          WRITE(9,'(5X,A20,F6.3)')' considering (B-V)o=',BV3
          WRITE(*,*)
          WRITE(*,'(A,$)')' HOW MANY BEST MODELS TO FIT IT          => '
          READ(5,*)IVOL
          WRITE(*,*)
          IF (NMOM.LE.IVOL) IVOL=NMOM
          CALL ZORAD(D,OMM,L,NMOM,IVOL,NMOM1)
          WRITE(*,*)' chosen from Kurucz models to fit (B-V)o'
          WRITE(9,*)' chosen from Kurucz models to fit (B-V)o'
          DO 300 J=1,IVOL
            WRITE(*,310)L(J),TEP(L(J)),GRAV(L(J)),AB(L(J)),BV(L(J))
     P      ,UB(L(J)),OMM(J)
            WRITE(9,310)L(J),TEP(L(J)),GRAV(L(J)),AB(L(J)),BV(L(J))
     P      ,UB(L(J)),OMM(J)
300       CONTINUE
          PAUSE
310       FORMAT(' PC:',I2,'  T=',F6.0,'  LOG G=',F4.2,'  LOG A=',F4.2,
     P    ' B-V=',F6.3,' U-B=',F6.3,'  DEV.:',F6.3)
        ENDIF
C---------------u-v-b-y-beta----
        IF(ABS(EBY1-1.E6).LT.1.E-10.AND.ABS(SUB-1.E6).GT.1.E-10.AND.
     P  ABS(SBY-1.E6).GT.1.E-10)THEN
          EBY=(13.608*SBY-SUB+1.467)/(13.608-1.53)
C       CASTELLI 1991 AA,251,106-116 for B stars only
        ELSE
          IF(ABS(EBY1-1.E6).LT.1.E-10)THEN
            EBY=0.D0
          ELSE
            EBY=EBY1
          ENDIF
        ENDIF
C       Philip et al. (1976)
C       SC0=SC1-0.25*EBY
C       SM0=SM1+0.32*EBY
C       SUB0=SUB-1.61*EBY
C       SBY0=SBY-EBY
C
        SC0=SC1-0.19*EBY
        SM0=SM1+0.33*EBY
        SUB0=SUB-1.53*EBY
        SBY0=SBY-EBY
C       CASTELLI 1991 AA,251,106-116
        call SM(SBY0,SBY,SM0,SC0,SC1,BETA)
        WRITE(*,*)
        WRITE(9,*)
        WRITE(*,*)'                 u v b y : '
        WRITE(9,*)'                 u v b y : '
        IF(ABS(SBY-1.E6).GT.1.E-10)THEN
          CALL NAPI4(SBY0)
        ENDIF
        IF(ABS(SC1-1.E6).GT.1.E-10)THEN
          IF(ABS(BETA-1.E6).GT.1.E-10.OR.ABS(WGAMA-1.E6).GT.1.E-10)THEN
            CALL BALONA(BETA,WGAMA,SC0)
            CALL CASTEL(SC0,BETA)
          ENDIF
          IF(ABS(SBY-1.E6).GT.1.E-10)THEN
            CALL NAPI2(SC1,SBY)
          ENDIF
        ENDIF
        IF(ABS(SUB-1.E6).GT.1.E-10.AND.ABS(SBY-1.E6).GT.1.E-10)THEN
          CALL NAPI1(SUB,SBY)
        ENDIF
        WRITE(*,*)
        WRITE(9,*)
        WRITE(*,'(10X,A41)')'  u v b y :      CORRECTED FOR REDDENING'
        WRITE(*,'(10X,2(A9,F6.3),2(A6,F6.3),A9,F6.3)')
     P  '  (u-b)o=',SUB0,'  (b-y)o=',SBY0,'   C0=',SC0,'   M0=',SM0,
     P  '  E(b-y)=',EBY
        WRITE(9,'(10X,A41)')'  u v b y :      CORRECTED FOR REDDENING'
        WRITE(9,'(10X,2(A9,F6.3),2(A6,F6.3),A9,F6.3)')
     P  '  (u-b)o=',SUB0,'  (b-y)o=',SBY0,'   C0=',SC0,'   M0=',SM0,
     P  '  E(b-y)=',EBY
        IF(ABS(SC1-1.E6).GT.1.E-10)THEN
          DO 350 I=1,NMOM
            D(I)=ABS(C1(I)-SC0)
            IF(ABS(D(I)-0.).LT.1.E-4) D(I)=0.0001
350       CONTINUE
          WRITE(*,*)
          WRITE(*,'(A,$)')' HOW MANY BEST MODELS TO FIT IT          => '
          READ(5,*)IVOL
          WRITE(*,*)
          IF (NMOM.LE.IVOL) IVOL=NMOM
          CALL ZORAD(D,OMM,L,NMOM,IVOL,NMOM1)
          WRITE(*,*)' chosen from Kurucz models to fit C0'
          WRITE(9,*)' chosen from Kurucz models to fit C0'
          DO 400 J=1,IVOL
            WRITE(*,410)L(J),TEP(L(J)),GRAV(L(J)),AB(L(J)),C1(L(J))
     P      ,OMM(J)
            WRITE(9,410)L(J),TEP(L(J)),GRAV(L(J)),AB(L(J)),C1(L(J))
     P      ,OMM(J)
400       CONTINUE
          PAUSE
410       FORMAT(' PC:',I2,'  T=',F6.0,'  LOG G=',F4.2,'  LOG A=',F4.2,
     P    ' C0=',F6.3,'  DEV.:',F6.3)
        ENDIF
        IF(ABS(SC1-1.E6).GT.1.E-10.AND.ABS(SBY-1.E6).GT.1.E-10)THEN
C       Philip et al. (1976)
          SC2=SC1-0.20*SBY
          SM2=SM1+0.18*SBY
          WRITE(*,*)
          WRITE(9,*)
          WRITE(*,'(10X,A41)')'  u v b y :   [c1] reddening independent'
          WRITE(*,'(10X,2(A6,F6.3))')
     P    ' [c1]=',SC2,' [m1]=',SM2
          WRITE(9,'(10X,A41)')'  u v b y :   [c1] reddening independent'
          WRITE(9,'(10X,2(A6,F6.3))')
     P    ' [c1]=',SC2,' [m1]=',SM2
          DO 450 I=1,NMOM
            C2(I)=C1(I)-0.20*BY(I)
            D(I)=ABS(C2(I)-SC2)
            IF(ABS(D(I)-0.).LT.1.E-4) D(I)=0.0001
450       CONTINUE
          WRITE(*,*)
          WRITE(*,'(A,$)')' HOW MANY BEST MODELS TO FIT IT          => '
          READ(5,*)IVOL
          WRITE(*,*)
          IF (NMOM.LE.IVOL) IVOL=NMOM
          CALL ZORAD(D,OMM,L,NMOM,IVOL,NMOM1)
          WRITE(*,*)' chosen from Kurucz models to fit [c1]'
          WRITE(9,*)' chosen from Kurucz models to fit [c1]'
          DO 460 J=1,IVOL
            WRITE(*,470)L(J),TEP(L(J)),GRAV(L(J)),AB(L(J)),C2(L(J))
     P      ,OMM(J)
            WRITE(9,470)L(J),TEP(L(J)),GRAV(L(J)),AB(L(J)),C2(L(J))
     P      ,OMM(J)
460       CONTINUE
          PAUSE
470       FORMAT(' PC:',I2,'  T=',F6.0,'  LOG G=',F4.2,'  LOG A=',F4.2,
     P    ' [c1]=',F6.3,'  DEV.:',F6.3)
        ENDIF
C---------------U-V-B1-B2-V1-G-----
        IF((ABS(GB2-1.E6).GT.1.E-10.AND.ABS(GV1-1.E6).GT.1.E-10).OR.
     P  (ABS(GB2-1.E6).GT.1.E-10.AND.ABS(GG-1.E6).GT.1.E-10))THEN
          WRITE(*,*)
          WRITE(9,*)
          WRITE(*,'(10X,A41)')'           GENEVA PHOTOMETRY :          '
          WRITE(9,'(10X,A41)')'           GENEVA PHOTOMETRY :          '
          CALL HAUCK(GB2,GV1,GG)
        ENDIF
        IF(ABS(GU-1.E6).GT.1.E-10.AND.ABS(GV-1.E6).GT.1.E-10.AND.
     P  ABS(GB1-1.E6).GT.1.E-10.AND.ABS(GB2-1.E6).GT.1.E-10.AND.
     P  ABS(GV1-1.E6).GT.1.E-10.AND.ABS(GG-1.E6).GT.1.E-10)THEN
          CALL CRAMER(GU,GV,GB1,GB2,GV1,GG)
          CALL MEGESS(GU,GV,GB1,GB2,GV1,GG)
          CALL NORTH(GU,GV,GB1,GB2,GV1,GG)
        ENDIF
        RETURN
        END
C-----------------------------------------------------------------------
c       SUBROUTINE HYLAND(UB0,BV0)
c       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
c       WRITE(*,*)
c       WRITE(9,*)
c       WRITE(*,*)' Hylandov vztah  from Qeff-Q    ?'
c       WRITE(9,*)' Hylandov vztah  from Qeff-Q    ?'
c       Q=UB0-0.645*BV0
c       QEFF=0.406*Q+0.51
c       TEFF=5040./QEFF
c       WRITE(*,'(10X,A6,F7.0)')' Teff=',TEFF
c       WRITE(9,'(10X,A6,F7.0)')' Teff=',TEFF
c       RETURN
c       END
C-----------------------------------------------------------------------
        SUBROUTINE HEINTZE(UB1,BV1,EBV,UB0,BV0)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        Q=UB1-0.645*BV1
        BV2=0.277*Q-0.045
        IF(BV2.LT.-0.32.OR.BV2.GT.-0.05)RETURN
        BV0=BV2
        EBV=BV1-BV0
        UB0=4.246*BV0+0.161
        WRITE(*,*)
        WRITE(9,*)
        WRITE(*,'(A38,A23)')' Heintze J.R.W., 1973, in Problems of ',
     P'calibration of absolute'
        WRITE(*,*)'  magnitudes and temperatures of stars, '
        WRITE(*,*)'  eds. Hauck B., Westelund B.E., IAU, p.231'
        WRITE(*,*)'     from reddening indep. Q for B1V-A0V'
        WRITE(9,'(A38,A23)')' Heintze J.R.W., 1973, in Problems of ',
     P'calibration of absolute'
        WRITE(9,*)'  magnitudes and temperatures of stars, '
        WRITE(9,*)'  eds. Hauck B., Westelund B.E., IAU, p.231'
        WRITE(9,*)'     from reddening indep. Q for B1V-A0V'
        WRITE(*,'(10X,3(A8,F6.3))')' (U-B)o=',UB0,' (B-V)o=',BV0,
     P' E(B-V)=',EBV
        WRITE(9,'(10X,3(A8,F6.3))')' (U-B)o=',UB0,' (B-V)o=',BV0,
     P' E(B-V)=',EBV
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE HYLAND(UB1,BV1)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION Q0(12),QEFF0(12)
        DATA Q0/0.0,-0.078,-0.148,-0.235,-0.334,-0.389,-0.445,-0.566
     P,-0.687,-0.743,-0.808,-0.864/
        DATA QEFF0/0.520,0.485,0.455,0.417,0.375,0.350,0.330,0.290
     P,0.246,0.225,0.197,0.175/
        Q=UB1-0.72*BV1
        IF(Q.LT.-0.864.OR.Q.GT.0.0)RETURN
        CALL INTEP(Q,QEFF,Q0,QEFF0,12,IER)
        TEFF=5040./QEFF
        WRITE(*,*)
        WRITE(9,*)
        WRITE(*,*)' Hyland A.R., 1969, in Theory and Observations of Nor
     pmal Stellar Atmospheres'
        WRITE(*,*)' p.271, The M.I.T. Press'
        WRITE(*,*)' from redd. indep.  Qeff-Q relation for B stars'
        WRITE(9,*)' Hyland A.R., 1969, in Theory and Observations of Nor
     pmal Stellar Atmospheres'
        WRITE(9,*)' p.271, The M.I.T. Press'
        WRITE(9,*)' from redd. indep.  Qeff-Q relation for B stars'
        WRITE(*,'(10X,A6,F7.0,5X,A6,F6.3)')' Teff=',TEFF,'for Q=',Q
        WRITE(9,'(10X,A6,F7.0,5X,A6,F6.3)')' Teff=',TEFF,'for Q=',Q
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE SCHILD(UB,BV)
C       IT DOES NOT MATTER WEATHER UB, BV OR UB0, BV0
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        WRITE(*,*)
        WRITE(9,*)
        WRITE(*,*)' Schild, R., Peterson, D., M., Oke, J., B., 1971,',
     P  ' ApJ 166, 95'
        WRITE(*,*)'       from Qeff-Q (red. indep. param.) for (B-V)<0'
        WRITE(9,*)' Schild, R., Peterson, D., M., Oke, J., B., 1971,',
     P  ' ApJ 166, 95'
        WRITE(9,*)'       from Qeff-Q (red. indep. param.) for (B-V)<0'
        Q=UB-0.72*BV
        IF(Q.LT.-0.1)THEN
          QEFF=0.378*Q+0.500
        ELSE
          QEFF=0.625*Q+0.525
        ENDIF
        TEFF=5040./QEFF
        WRITE(*,'(10X,A6,F7.0,5X,A6,F6.3)')' Teff=',TEFF,'for Q=',Q
        WRITE(9,'(10X,A6,F7.0,5X,A6,F6.3)')' Teff=',TEFF,'for Q=',Q
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE UNDERHILL(BV0)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION BV00(14),TEFF0(14)
c       CHARACTER*4 SP(14)
        DATA BV00/-0.30,-0.28,-0.26,-0.25,-0.24,-0.22,-0.20,-0.18,-0.16,
     P-0.14,-0.13,-0.10,-0.07,-0.04/
        DATA TEFF0/3.078,2.927,2.690,2.572,2.282,2.038,1.853,1.634,1.517
     P,1.374,1.298,1.190,1.092,1.020/
c       DATA SP/'B0  ','B0.5','B1  ','B1.5','B2  ','B2.5','B3  ','B4  ',
c    P'B5  ','B6  ','B7  ','B8  ','B9  ','B9.5'/
        IF(BV0.LT.-0.30.OR.BV0.GT.-0.04)RETURN
        WRITE(*,*)
        WRITE(9,*)
        WRITE(*,*)' Underhill A., Doazan V., 1982, "B stars with and wit
     Phout emision lines, NASA'
        WRITE(*,*)' using intrinsic colors (B-V)o from FitzGerald(1970)'
        WRITE(*,*)' just for BIV-V stars'
        WRITE(9,*)' Underhill A., Doazan V., 1982, "B stars with and wit
     Phout emision lines, NASA'
        WRITE(9,*)' using intrinsic colors (B-V)o from FitzGerald(1970)'
        WRITE(9,*)' just for BIV-V stars'
        CALL INTEP(BV0,TEFF,BV00,TEFF0,14,IER)
        TEFF=TEFF*1.D4
        WRITE(*,'(10X,A6,F7.0,5X,A11,F6.3)')' Teff=',TEFF,'for (B-V)o='
     P,BV0
        WRITE(9,'(10X,A6,F7.0,5X,A11,F6.3)')' Teff=',TEFF,'for (B-V)o='
     P,BV0
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE HHSA(UB0)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        WRITE(*,*)
        WRITE(9,*)
        WRITE(*,'(A38,A23)')' Heintze J.R.W., 1973, in Problems of ',
     P'calibration of absolute'
        WRITE(*,*)'  magnitudes and temperatures of stars, '
        WRITE(*,*)'  eds. Hauck B., Westelund B.E., IAU, p.231'
        WRITE(*,*)'     from reddening indep. Q for B1V-A0V'
        WRITE(*,*)'  calibrated results of different authors on (U-B)o'
        WRITE(9,'(A38,A23)')' Heintze J.R.W., 1973, in Problems of ',
     P'calibration of absolute'
        WRITE(9,*)'  magnitudes and temperatures of stars, '
        WRITE(9,*)'  eds. Hauck B., Westelund B.E., IAU, p.231'
        WRITE(9,*)'     from reddening indep. Q for B1V-A0V'
        WRITE(9,*)'  calibrated results of different authors on (U-B)o'
        QEFF=0.3145*UB0+0.513
        TEFF=5040./QEFF
        WRITE(*,'(10X,A18,A6,F7.0)')' from Hyland(1969)',' Teff=',TEFF
        WRITE(9,'(10X,A18,A6,F7.0)')' from Hyland(1969)',' Teff=',TEFF
        QEFF=0.306*UB0+0.507
        TEFF=5040./QEFF
        WRITE(*,'(10X,A19,A6,F7.0)')' from Heintze(1969)',' Teff=',TEFF
        WRITE(9,'(10X,A19,A6,F7.0)')' from Heintze(1969)',' Teff=',TEFF
        QEFF=0.3045*UB0+0.502
        TEFF=5040./QEFF
        WRITE(*,'(10X,A31,A26,F7.0)')' from Auer & Mihalas(1972) and ',
     p'Schild et.al(1971)   Teff=',TEFF
        WRITE(9,'(10X,A31,A26,F7.0)')' from Auer & Mihalas(1972) and ',
     p'Schild et.al(1971)   Teff=',TEFF
        WRITE(*,'(A15,F6.3)')' all for (U-B)o=',UB0
        WRITE(9,'(A15,F6.3)')' all for (U-B)o=',UB0
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE HAUCK(GB2,GV1,GG)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        WRITE(*,*)
        WRITE(9,*)
        WRITE(*,*)' Hauck, B., North, P., 1982, AA 114, 23',
     P  '   for normal B and A stars'
        WRITE(9,*)' Hauck, B., North, P., 1982, AA 114, 23',
     P  '   for normal B and A stars'
        IF(ABS(GB2-1.E6).GT.1.E-10.AND.ABS(GV1-1.E6).GT.1.E-10)THEN
          IF((GB2-GV1).GT.-0.172)THEN
            TEFF=5040./(0.635*(GB2-GV1)+0.634)
            WRITE(*,'(10X,A6,F7.0,A12)')' Teff=',TEFF,'     z B2-V1'
            WRITE(9,'(10X,A6,F7.0,A12)')' Teff=',TEFF,'     z B2-V1'
          ENDIF
          IF((GB2-GV1).LT.-0.138)THEN
            TEFF=5040./(1.8*(GB2-GV1)+0.798)
            WRITE(*,'(10X,A6,F7.0,A12)')' Teff=',TEFF,'     z B2-V1'
            WRITE(9,'(10X,A6,F7.0,A12)')' Teff=',TEFF,'     z B2-V1'
          ENDIF
        ENDIF
        IF(ABS(GB2-1.E6).GT.1.E-10.AND.ABS(GG-1.E6).GT.1.E-10)THEN
          IF((GB2-GG).GT.-0.677)THEN
            TEFF=5040./(0.458*(GB2-GG)+0.835)
            WRITE(*,'(10X,A6,F7.0,A12)')' Teff=',TEFF,'      z B2-G'
            WRITE(9,'(10X,A6,F7.0,A12)')' Teff=',TEFF,'      z B2-G'
          ENDIF
          IF((GB2-GG).LT.-0.632)THEN
            TEFF=5040./(1.362*(GB2-GG)+1.411)
            WRITE(*,'(10X,A6,F7.0,A12)')' Teff=',TEFF,'      z B2-G'
            WRITE(9,'(10X,A6,F7.0,A12)')' Teff=',TEFF,'      z B2-G'
          ENDIF
        ENDIF
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE CRAMER(U,V,B1,B2,V1,G)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION XI(19),BCI(19),X0I(15),A0I(15),A1I(15),A2I(15)
        DATA X0I/0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,
     P  1.7,1.8/
        DATA A0I/-2.16,-1.65,-1.26,-0.94,-0.66,-0.40,-0.13,0.15,0.43,
     P  0.69,0.92,1.06,1.07,0.89,0.45/
        DATA A1I/-19.24,-16.62,-13.88,-11.31,-9.11,-7.40,-6.25,-5.63,
     P  -5.44,-5.50,-5.56,-5.30,-4.31,-2.10,1.87/
        DATA A2I/-424.34,-347.14,-327.90,-285.62,-247.29,-212.92,
     P  -182.51,-156.05,-133.55,-115.01,-100.43,-89.80,-83.13,-80.42,
     P  -81.67/
        DATA XI/-0.076,-0.020,0.042,0.110,0.185,0.269,0.364,0.473,
     P  0.600,0.752,0.840,0.940,1.054,1.188,1.349,1.552,1.678,1.836,
     P  2.044/
        DATA BCI/-3.21,-3.06,-2.91,-2.75,-2.58,-2.41,-2.23,-2.02,-1.80,
     P  -1.52,-1.35,-1.14,-0.90,-0.68,-0.48,-0.28,-0.18,-0.08,-0.02/
        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
        IF(X.GT.-0.076.AND.X.LT.2.044)THEN
          TEFF=10**(4.496-0.453*X+0.086*X*X)
          CALL INTEP(X,BC,XI,BCI,19,IER)
        ELSE
          RETURN
        ENDIF
        BETA=2.568+0.190*X-0.487*Y
        IF(X.GT.0.4.AND.X.LT.1.8.AND.Y.GT.-0.06)THEN
          CALL INTEP(X,A0,X0I,A0I,15,IER)
          CALL INTEP(X,A1,X0I,A1I,15,IER)
          CALL INTEP(X,A2,X0I,A2I,15,IER)
          AMV=A0+A1*Y+A2*Y*Y
        ELSE
          AMV=1.E6
        ENDIF
        WRITE(*,*)
        WRITE(9,*)
        WRITE(*,*)' Cramer, N., Maeder, A., 1979, AA 78, 305'
        WRITE(*,*)'      for B-type stars, from Geneva. phot.: redd.inde
     Pp. X-calib.'
        WRITE(9,*)' Cramer, N., Maeder, A., 1979, AA 78, 305'
        WRITE(9,*)'      for B-type stars, from Geneva. phot.: redd.inde
     Pp. X-calib.'
        WRITE(*,'(10X,A6,F7.0,2(A5,F7.3),A7,F5.3)')' Teff=',TEFF,
     P  '  Mv=',AMV,'  BC=',BC,'  BETA=',BETA
        WRITE(9,'(10X,A6,F7.0,2(A5,F7.3),A7,F5.3)')' Teff=',TEFF,
     P  '  Mv=',AMV,'  BC=',BC,'  BETA=',BETA
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE MEGESS(U,V,B1,B2,V1,G)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION XM(5),TM(5)
        DATA XM/1.55,1.33,1.15,0.96,0.78/
        DATA TM/10100.,11325.,12570.,14000.,15450./
        X=0.3788+1.3764*U-1.2162*B1-0.8498*B2-0.1554*V1+0.8450*G
        IF(X.GT.0.529.AND.X.LT.1.55)THEN
          WRITE(*,*)
          WRITE(9,*)
          WRITE(*,*)' Megessier C., 1988, AA 206, 74'
          WRITE(9,*)' Megessier C., 1988, AA 206, 74'
          WRITE(*,*)'        from Geneva phot., X-calib.'
          WRITE(9,*)'        from Geneva phot., X-calib.'
        ENDIF
        IF(X.GT.0.78.AND.X.LT.1.55)THEN
          CALL INTEP(X,TEFF,XM,TM,5,IER)
          WRITE(*,'(10X,A6,F7.0,A17)')' Teff=',TEFF,'  for norm. stars'
          WRITE(9,'(10X,A6,F7.0,A17)')' Teff=',TEFF,'  for norm. stars'
        ENDIF
        IF(X.GT.0.529.AND.X.LT.1.429)THEN
          TEFF=10**(4.3373-0.2516*X)
          WRITE(*,'(10X,A6,F7.0,A15)')' Teff=',TEFF,'  for CP2 stars'
          WRITE(9,'(10X,A6,F7.0,A15)')' Teff=',TEFF,'  for CP2 stars'
        ENDIF
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE NORTH(U,V,B1,B2,V1,G)
        PARAMETER (NX=13,NY=22)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION XI(NY),YI(NX),Q(NY,NX),GLOG(NY,NX)
        DATA XI/0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,
     P  1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1/
        DATA YI/-0.050,-0.025,0.0,0.025,0.050,0.075,0.1,0.125,0.15,
     P  0.175,0.2,0.225,0.25/
        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
        WRITE(*,*)
        WRITE(9,*)
        WRITE(*,*)' North, P., Nicolet, B., 1990, AA 228, 78'
        WRITE(*,*)'             for B stars, from X-Y plane'
        WRITE(9,*)' North, P., Nicolet, B., 1990, AA 228, 78'
        WRITE(9,*)'             for B stars, from X-Y plane'
        WRITE(*,'(10X,2(A5,F6.3))')'   X=',X,'   Y=',Y
        WRITE(9,'(10X,2(A5,F6.3))')'   X=',X,'   Y=',Y
        IF(X.LT.0.0.OR.X.GT.2.1.OR.Y.LT.-0.05.OR.Y.GT.0.25)THEN
          RETURN
        ENDIF
        OPEN(8,FILE='NORTH.DAT',STATUS='OLD')
        DO 10,I=1,NY
          READ(8,*)(Q(I,J),J=1,NX)
          READ(8,*)(GLOG(I,J),J=1,NX)
10      CONTINUE
        CLOSE(8)
        DO 20,I=1,NY
          IF(X.LT.XI(I))THEN
            KKX=I-1
            GOTO 25
          ENDIF
20      CONTINUE
25      DO 30,I=1,NX
          IF(Y.LT.YI(I))THEN
            KKY=I-1
            GOTO 35
          ENDIF
30      CONTINUE
35      IF(ABS(Q(KKX,KKY)-0.0).LT.1.E-10.OR.
     P  ABS(Q(KKX+1,KKY)-0.0).LT.1.E-10.OR.
     P  ABS(Q(KKX,KKY+1)-0.0).LT.1.E-10.OR.
     P  ABS(Q(KKX,KKY+1)-0.0).LT.1.E-10)THEN
          RETURN
        ENDIF
        POMQ1=(Y-YI(KKY))*(Q(KKX,KKY+1)-Q(KKX,KKY))/(YI(KKY+1)-YI(KKY))
     P  +Q(KKX,KKY)
        POMQ2=(Y-YI(KKY))*(Q(KKX+1,KKY+1)-Q(KKX+1,KKY))
     P  /(YI(KKY+1)-YI(KKY))+Q(KKX+1,KKY)
        TEFF=5040./((X-XI(KKX))*(POMQ2-POMQ1)/(XI(KKX+1)-XI(KKX))+POMQ1)
        POMG1=(Y-YI(KKY))*(GLOG(KKX,KKY+1)-GLOG(KKX,KKY))/
     P  (YI(KKY+1)-YI(KKY))+GLOG(KKX,KKY)
        POMG2=(Y-YI(KKY))*(GLOG(KKX+1,KKY+1)-GLOG(KKX+1,KKY))/
     P  (YI(KKY+1)-YI(KKY))+GLOG(KKX+1,KKY)
        GRAV=(X-XI(KKX))*(POMG2-POMG1)/(XI(KKX+1)-XI(KKX))+POMG1
        WRITE(*,'(10X,A6,F7.0,A8,F4.2,A23)')' Teff=',TEFF,'  log g=',
     P  GRAV,'  for solar metallicity'
        WRITE(9,'(10X,A6,F7.0,A8,F4.2,A23)')' Teff=',TEFF,'  log g=',
     P  GRAV,'  for solar metallicity'
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE BALONA(BETA1,WGAMA,C00)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        WRITE(*,*)
        WRITE(9,*)
        WRITE(*,*)' Balona, L., A., 1984, Mon. Not. R. astr. Soc.,',
     P  '211, 973'
        WRITE(*,*)'           for early-type stars   from beta and c0'
        WRITE(9,*)' Balona, L., A., 1984, Mon. Not. R. astr. Soc.,',
     P  '211, 973'
        WRITE(9,*)'           for early-type stars   from beta and c0'
        IF(ABS(BETA1-1.E6).LT.1.E-10)THEN
          BETA=0.02565*WGAMA+2.5123
          WRITE(*,'(A7,F5.3)')'  BETA=',BETA
          WRITE(9,'(A7,F5.3)')'  BETA=',BETA
        ELSE
          BETA=BETA1
        ENDIF
        BET=ALOG10(BETA-2.500)
        C00T=ALOG10(C00+0.200)
        TEFF=10**(3.9036-0.4816*C00T-0.5290*BET-
     P  0.1260*C00T*C00T+0.0924*BET*C00T-0.4013*BET*BET)
C       FROM KURUCZ MODELS AND BETA-WGAMA(KURUCZ-LTE) CALIB.
        GRAV=5.9046-3.2262*C00T+4.0883*BET-0.5383*C00T*C00T-
     P  0.2774*BET*C00T-0.0007*BET*BET
        EG=ALOG10(BETA-2.515)-1.60*ALOG10(C00+0.322)
        AMV=3.499+7.203*ALOG10(BETA-2.515)-2.319*EG+2.938*EG*EG*EG
        BC=0.2900+2.8467*C00T+2.8334*BET+0.6481*C00T*C00T
     P  -0.2997*BET*C00T+2.1487*BET*BET
        WRITE(*,'(10X,A6,F7.0,A8,F4.2,2(A5,F7.3))')' Teff=',TEFF,
     P  '  log g=',GRAV,'  Mv=',AMV,'  BC=',BC
        WRITE(9,'(10X,A6,F7.0,A8,F4.2,2(A5,F7.3))')' Teff=',TEFF,
     P  '  log g=',GRAV,'  Mv=',AMV,'  BC=',BC
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE CASTEL(C0,BETA)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION P(10,2),Q(7,2)
        DATA Q/9.372853735692513,0.1388984909726595,-0.3772751850721322,
     P  -3.441202247069613,-1.2919132663971867E-2,0.1563028804655211,
     P  0.5329196601502442,41.12891680060272,0.3074969330177897,
     P  0.1249648644440178,-25.75280684590345,-1.8620519285948986E-2,
     P  -0.1306123353098876,4.340621513119619/
        DATA P/-34699.13571279180,-9.389628024837989,-14.85477894733068,
     P  49675.89762128829,17.42254584910010,-26659.84653847129,
     P  -14.25166444541836,6364.673937260428,4.924571477326692,
     P  -570.2261422955868,-102313.9917581691,-4.047990456431145,
     P  -8.706895982165089,150275.8361010616,17.53795434313011,
     P  -82781.98328595359,-44.19718767241700,20273.39016383133,
     P  44.95052473400207,-1862.257533577166/
        WRITE(*,*)
        WRITE(*,*)' Castelli, F., 1991, AA 251, 106 from C0,BETA',
     P  ' for B-stars'
5       WRITE(*,'(A,$)')' Teff initial ><15000 (TO SKIP -0-)  => '
        READ(*,*)TEFF
        IF(DABS(TEFF-0.).LT.1.E-30)GOTO 20
        WRITE(9,*)
        WRITE(9,*)' Castelli, F., 1991, AA 251, 106 from C0,BETA',
     P  ' for B-stars'
        WRITE(*,'(A,$)')' No. OF ITERATION => '
        READ(*,*)N
        IF(TEFF.LT.15000)THEN
          J=1
        ELSE
          J=2
        ENDIF
        TEFFI=DLOG10(TEFF)
        DO 10 I=1,N
          GRAVI=P(1,J)+P(2,J)*TEFFI+P(3,J)*C0+P(4,J)*BETA+P(5,J)*C0*C0+
     P    P(6,J)*BETA*BETA+P(7,J)*C0*C0*C0+P(8,J)*BETA*BETA*BETA+
     P    P(9,J)*C0**4.+P(10,J)*BETA**4.
          TEFFI=Q(1,J)+Q(2,J)*GRAVI+Q(3,J)*C0+Q(4,J)*BETA+
     P    Q(5,J)*GRAVI*GRAVI+Q(6,J)*C0*C0+Q(7,J)*BETA*BETA
          WRITE(*,*)' TEFF, LOGG',TEFFI,GRAVI
10      CONTINUE
        WRITE(*,'(10X,A6,F7.0,A8,F4.2,A19,A10,F7.0)')' Teff=',
     P  DEXP(TEFFI*2.302585D0),
     P  '  log g=',GRAVI,'  for 11000<T<20000',' Tinitial=',TEFF
        WRITE(9,'(10X,A6,F7.0,A8,F4.2,A19,A10,F7.0)')' Teff=',
     P  DEXP(TEFFI*2.302585D0),
     P  '  log g=',GRAVI,'  for 11000<T<20000',' Tinitial=',TEFF
        GOTO 5
20      RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE NAPI1(UB1,BY1)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        WRITE(*,*)
        WRITE(9,*)
        WRITE(*,*)' Napiwotzky R., Schonberner D., Wenske V., 1993,',
     P  ' AA 268, 653'
        WRITE(*,*)'      for 9500-30000K'
        WRITE(9,*)' Napiwotzky R., Schonberner D., Wenske V., 1993,',
     P  ' AA 268, 653'
        WRITE(9,*)'      for 9500-30000K'
        UBH=UB1-1.53*BY1
        QEFF=0.1692+0.2828*UBH-0.0195*UBH*UBH
        TEFF=5040.0/QEFF
        WRITE(*,'(10X,A6,F7.0,A17,A11)')' Teff=',TEFF,' for norm. stars'
     P,' from [u-b]'
        WRITE(9,'(10X,A6,F7.0,A17,A11)')' Teff=',TEFF,' for norm. stars'
     P,' from [u-b]'
        QEFF=0.2162+0.2301*UBH
        TEFF=5040.0/QEFF
        WRITE(*,'(10X,A6,F7.0,A15,A11)')' Teff=',TEFF,'  for CP2 stars'
     P,' from [u-b]'
        WRITE(9,'(10X,A6,F7.0,A15,A11)')' Teff=',TEFF,'  for CP2 stars'
     P,' from [u-b]'
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE NAPI2(C1,BY1)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        WRITE(*,*)
        WRITE(9,*)
        WRITE(*,*)' Napiwotzky R., Schonberner D., Wenske V., 1993,',
     P  ' AA 268, 653'
        WRITE(9,*)' Napiwotzky R., Schonberner D., Wenske V., 1993,',
     P  ' AA 268, 653'
        C1H=C1-0.20*BY1
C       Philip et al. (1976)
        QEFF=0.2489+0.2698*C1H
        TEFF=5040.0/QEFF
        WRITE(*,'(10X,A6,F7.0,A15,A10)')' Teff=',TEFF,'  for CP2 stars'
     P,' from [C1]'
        WRITE(9,'(10X,A6,F7.0,A15,A10)')' Teff=',TEFF,'  for CP2 stars'
     P,' from [C1]'
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE NAPI4(SBY0)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        IF (SBY0.LT.-0.0215)THEN
          Q=(SBY0+0.175)/0.307
        ELSE
          IF (SBY0.GT.0.0494)THEN
            Q=(SBY0+0.745)/1.324
          ELSE
            Q=(SBY0+0.376)/0.709
          ENDIF
        ENDIF
        TEFF=5040.0/Q
        WRITE(*,*)
        WRITE(9,*)
        WRITE(*,*)' Napiwotzky R., Schonberner D., Wenske V., 1993,',
     P  ' AA 268, 653'
        WRITE(*,*)'      for B, A, F stars'
        WRITE(9,*)' Napiwotzky R., Schonberner D., Wenske V., 1993,',
     P  ' AA 268, 653'
        WRITE(9,*)'      for B, A, F stars'
        WRITE(*,'(10X,A6,F7.0,A17,A13,F6.3)')' Teff=',TEFF,
     P  ' for norm. stars',' from (b-y)o=',SBY0
        WRITE(9,'(10X,A6,F7.0,A17,A13,F6.3)')' Teff=',TEFF,
     P  ' for norm. stars',' from (b-y)o=',SBY0
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE NAPI3(BV0)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        IF (BV0.LT.-0.0505)THEN
          Q=(BV0+0.384)/0.667
        ELSE
          IF (BV0.GT.0.0842)THEN
            Q=(BV0+1.114)/1.997
          ELSE
            Q=(BV0+0.724)/1.347
          ENDIF
        ENDIF
        TEFF=5040.0/Q
        WRITE(*,*)
        WRITE(9,*)
        WRITE(*,*)' Napiwotzky R., Schonberner D., Wenske V., 1993,',
     P  ' AA 268, 653'
        WRITE(*,*)'      for B, A, F stars'
        WRITE(9,*)' Napiwotzky R., Schonberner D., Wenske V., 1993,',
     P  ' AA 268, 653'
        WRITE(9,*)'      for B, A, F stars'
        WRITE(*,'(10X,A6,F7.0,A17,A13,F6.3)')' Teff=',TEFF,
     P  ' for norm. stars',' from (B-V)o=',BV0
        WRITE(9,'(10X,A6,F7.0,A17,A13,F6.3)')' Teff=',TEFF,
     P  ' for norm. stars',' from (B-V)o=',BV0
        RETURN
        END
C-----------------------------------------------------------------------
        SUBROUTINE SM(SBY0,SBY,SM0,SC0,SC1,BETA)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        A0=1.36*SBY0+0.36*SM0+0.18*SC0-0.2448
        R=0.35*SC1-0.07*SBY-BETA+2.565
        IF (A0.LT.-0.09.OR.A0.GT.0.16.OR.R.LT.-0.03.OR.R.GT.0.34)THEN
          RETURN
        ENDIF
        WRITE(*,*)
        WRITE(9,*)
        WRITE(*,'(10X,A6,F6.3,A13,F6.3)')'   a0=',A0,
     P  '         r* =',R
        WRITE(9,'(10X,A6,F6.3,A13,F6.3)')'   a0=',A0,
     P  '         r* =',R
        WRITE(*,*)' to be used with intermediate A star grid e.g.'
        WRITE(*,'(A,A)')' Smalley B., Dworetsky M.M., 1994,',
     P  ' poster in GA22, Hague'
        WRITE(9,*)' to be used with intermediate A star grid e.g.'
        WRITE(9,'(A,A)')' Smalley B., Dworetsky M.M., 1994,',
     P  ' poster in GA22, Hague'
        RETURN
        END