	program makemarcs1
c
c a routine to generate an interpolated MARCS model
c for Snedens MOOG
c
c A.McWilliam

        real*8 teff,logg,feh,micro
        character*10 mtype

	write(6,'(a)')'Enter Teff, logg, [Fe/H] and micro.'
        read(5,*)teff,logg,feh,micro

c
c interpolate model of specified type
c
        write(6,*)
     %  'enter MARCSsp for standard comp. plane-parallel atmosphere'
        write(6,*)
     %  'enter MARCSss for standard comp. 1D spherical atmosphere'
        read(5,'(a)')mtype
        if (mtype.ne.'MARCSsp' .and.  mtype.ne.'MARCSSP' .and.
     %      mtype.ne.'MARCSss' .and.  mtype.ne.'MARCSSS' ) then
           write(6,*)
     %    'Invalid type: must be MARCSsp or MARCSss'
           stop
        endif
c

        call make_marcs(teff,logg,feh,micro,mtype)
	stop
	end





      SUBROUTINE MAKE_MARCS(TEFF,LOGG,FE,MICRO,MTYPEIN)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL*8 TEFF,LOGG,FE,MICRO
      LOGICAL ERROR
      CHARACTER*80 TITLE,MTYPEIN*10
      COMMON/MODTYPE/ MTYPE
      CHARACTER*10 MTYPE

      TITLE = '          '
      ERROR = .FALSE.
      MTYPE = MTYPEIN

c write parameters
      WRITE(TITLE,100)TEFF,LOGG,FE,MICRO
 100  FORMAT(F7.0,1h/,F7.2,1h/,F7.2,5x,7h mic = ,F6.4)


      CALL GET_MODEL_MARCS(TEFF,LOGG,FE,ERROR)
c
      IF (ERROR) THEN
          WRITE(6,'(a)')'ERROR IN MAKE_MARCS SUBROUTINE '
          STOP 
      ENDIF
      CALL SET_UP_LINES_INPUT_FILE(MICRO,FE,TITLE)
      CLOSE(UNIT=7)
      RETURN
      END





      
      SUBROUTINE GET_MODEL_MARCS(TEFF,LOGG,FE,ERROR)

      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 LOGG,TEFF,FE,METAL(15),M1,M2
      LOGICAL ERROR
      CHARACTER*8 M1FIL,M2FIL,MODEL
      DATA METAL/+1.0,+0.75,+0.5,+0.25,0.00,-0.25,-0.50,-0.75,-1.0,
     %   -1.5,-2.0,-2.5,-3.0,-4.0,-5.0/

c     DATA M1FIL,M2FIL,MODEL/'M1      ','M2      ','model   '/
      
      MODEL = 'model   '
      M1FIL = 'M1      '
      M2FIL = 'M2      '

      IF (FE .GT. +1.0D0) THEN   
         WRITE(6,'(a)')
     %      'STAR IS MORE METAL RICH THAN +1.0: +1.0 MODEL USED'
         FE = 1.0D0   
      ENDIF
      IF (FE .LT. -5.0D0)   THEN
         WRITE(6,'(a)')
     %      'STAR IS MORE METAL POOR THAN -5.0: -5.0 MODEL USED'
         FE = -5.0D0  
      ENDIF
      IF (FE .EQ. 1.00 .OR. FE .EQ. 0.75 .OR. FE .EQ.  0.50 .OR.
     %    FE .EQ. 0.25 .OR. FE .EQ. 0.00 .OR. FE .EQ. -0.25 .OR.
     %    FE .EQ.-0.50 .OR. FE .EQ.-0.75 .OR. FE .EQ. -1.00 .OR.
     %    FE .EQ.-1.50 .OR. FE .EQ.-2.00 .OR. FE .EQ. -2.50 .OR.
     %    FE .EQ.-3.00 .OR. FE .EQ.-4.00 .OR. FE .EQ. -5.00 ) THEN

         CALL COMPUTE_MODEL_ATMOSPHERE(TEFF,LOGG,FE,MODEL,ERROR)   
      ELSE
         DO I = 2, 15
            IF (FE .GT. METAL(I)) THEN
               M1 = METAL(I-1)
               M2 = METAL(I)
               GOTO 800
            ENDIF
         ENDDO
 800     CONTINUE
         CALL COMPUTE_MODEL_ATMOSPHERE(TEFF,LOGG,M1,M1FIL,ERROR)
         CALL COMPUTE_MODEL_ATMOSPHERE(TEFF,LOGG,M2,M2FIL,ERROR)
         IF (ERROR)  RETURN
         CALL KINTERP(M1,M2,FE,M1FIL,M2FIL,MODEL)
      ENDIF
      RETURN
      END







      SUBROUTINE COMPUTE_MODEL_ATMOSPHERE(TEFF,LOGG,FE,OUTFILE,ERROR)
c
c determines the flow of models to be interpolated
c

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*8 OUTFILE
      REAL*8 TEFF,TEFFS(2),LOGG,GRAVS(2,2),FE
      LOGICAL OKG,OKT,ERROR

      CALL PARAM_BOUNDS(FE,TEFF,LOGG,TEFFS,GRAVS)

      IF (TEFF.EQ.TEFFS(1)) THEN
         IF (LOGG.EQ.GRAVS(1,1) .AND. LOGG.EQ.GRAVS(1,2)) THEN
            CALL PRINT_MODEL(TEFF,LOGG,OUTFILE,FE,ERROR)
         ELSE
            CALL PRINT_MODEL(TEFF,GRAVS(1,1),'MOD1    ',FE,ERROR)
            CALL PRINT_MODEL(TEFF,GRAVS(1,2),'MOD2    ',FE,ERROR)
            IF (ERROR) RETURN
            CALL KINTERP(GRAVS(1,1),GRAVS(1,2),LOGG,'MOD1    ','MOD2    ',OUTFILE)
         ENDIF
      ELSE
         IF (LOGG.EQ.GRAVS(1,1)) THEN
            CALL PRINT_MODEL(TEFFS(1),LOGG,'MOD1    ',FE,ERROR)
            IF (LOGG.EQ.GRAVS(2,1)) THEN
               CALL PRINT_MODEL(TEFFS(2),LOGG,'MOD2    ',FE,ERROR)
               IF (ERROR) RETURN
               CALL KINTERP(TEFFS(1),TEFFS(2),TEFF,'MOD1    ','MOD2    ',OUTFILE)
            ELSE
               CALL PRINT_MODEL(TEFFS(2),GRAVS(2,1),'MOD2    ',FE,ERROR)
               CALL PRINT_MODEL(TEFFS(2),GRAVS(2,2),'MOD3    ',FE,ERROR)
               IF (ERROR) RETURN
               CALL KINTERP(GRAVS(2,1),GRAVS(2,2),LOGG,'MOD2    ','MOD3 ','MOD4    ')
               CALL KINTERP(TEFFS(1),TEFFS(2),TEFF,'MOD1    ','MOD4 ',OUTFILE)
            ENDIF
         ELSEIF (LOGG.EQ.GRAVS(2,1)) THEN
            CALL PRINT_MODEL(TEFFS(2),LOGG,'MOD1    ',FE,ERROR)
            CALL PRINT_MODEL(TEFFS(1),GRAVS(1,1),'MOD2    ',FE,ERROR)
            CALL PRINT_MODEL(TEFFS(1),GRAVS(1,2),'MOD3    ',FE,ERROR)
            IF (ERROR) RETURN
            CALL KINTERP(GRAVS(1,1),GRAVS(1,2),LOGG,'MOD2    ','MOD3    ','MOD4    ')
            CALL KINTERP(TEFFS(2),TEFFS(1),TEFF,'MOD1    ','MOD4    ',OUTFILE)
         ELSE
            CALL PRINT_MODEL(TEFFS(1),GRAVS(1,1),'MOD1    ',FE,ERROR)
            CALL PRINT_MODEL(TEFFS(1),GRAVS(1,2),'MOD2    ',FE,ERROR)
            CALL PRINT_MODEL(TEFFS(2),GRAVS(2,1),'MOD3    ',FE,ERROR)
            CALL PRINT_MODEL(TEFFS(2),GRAVS(2,2),'MOD4    ',FE,ERROR)
            IF (ERROR) RETURN
            CALL KINTERP(GRAVS(1,1),GRAVS(1,2),LOGG,'MOD1    ','MOD2    ','MOD5    ')
            CALL KINTERP(GRAVS(2,1),GRAVS(2,2),LOGG,'MOD3    ','MOD4    ','MOD6    ')
            CALL KINTERP(TEFFS(1),TEFFS(2),TEFF,'MOD5    ','MOD6    ',OUTFILE)
         ENDIF
      ENDIF

      RETURN
      END




        subroutine param_bounds(fe,tin,gin,teffs,gravs)
c
c a routine to find the nearest models to a specified Teff and logg
c in a given MARCS file
c
        implicit real*8 (a-h,o-z)

        real*8 t(1000),g(1000),gravs(2,2),teffs(2)
        real*8 tin,gin,glo,ghi,tlo,thi,fe
        logical err


        call read_params(fe,t,g,n)

c first find the closest temperatures
        call get_teffs(t,n,tin,teffs,err)

        if (err) then
           write(6,*)'Teff out of bounds'
           stop
        else
c          write(6,*)'Teff,  T_lo, T_hi'
c          write(6,*)tin,teffs(1),teffs(2)

c now find the closest gravities for given teff=tlo
           tlo = teffs(1)
           call get_gravs(t,g,n,tlo,gin,glo,ghi,err)
           gravs(1,1) = glo
           gravs(1,2) = ghi
           if (err) then
              write(6,*)'logg out of bounds for T_lo'
              stop
           else
c             write(6,*)'T_lo,  log g_lo, log g_hi'
c             write(6,*)tlo,gravs(1,1),gravs(1,2)
           endif

c now find the closest gravities for given teff=thi
           thi = teffs(2)
           call get_gravs(t,g,n,thi,gin,glo,ghi,err)
           gravs(2,1) = glo
           gravs(2,2) = ghi
           if (err) then
              write(6,*)'logg out of bounds for T_hi'
              stop
           else
c             write(6,*)'T_lo,  log g_lo, log g_hi'
c             write(6,*)thi,gravs(2,1),gravs(2,2)
           endif
        endif

        return
        end



        subroutine read_params(fe,t,g,n)
c
c simply read the Teffs and log g values in the MARCS files
c
        character*80 line,fname
        real*8 t(1000),g(1000),fe
        logical error
        integer n

        call access_model_file(fe,fname,error)
        n = 0
        do j = 1,10000000
           read(11,'(a)',end=900)line
           if (line(:5).eq.'MODEL') then
              n = n + 1
              read(11,*)t(n),g(n)
           endif
        enddo
        stop 'eeek! EOF not reached'

 900    continue
        rewind(unit=11)
        close(unit=11)

        return
        end




        subroutine get_teffs(t,n,tin,teffs,err)
c
c first find the closest temperatures
c teffs(1) = tlo, teffs(2)=thi
c
        real*8 t(1000),tin,dtlo,dthi
        real*8 teffs(2)
        logical err
        integer n

        teffs(1) = 0.0
        teffs(2) = 0.0
        dtlo = 10000.
        dthi = 10000.

        err = .false.

        do i = 1,n
           if (t(i).eq.tin) then
              teffs(1) = t(i)
              teffs(2) = t(i)
              return
           elseif (t(i).lt.tin) then
              if ( (tin-t(i)).le.dtlo) then
                 dtlo = tin-t(i)
                 teffs(1) = t(i)
              endif
           elseif (t(i).gt.tin) then
              if ( (t(i)-tin).le.dthi) then
                 dthi = t(i)-tin
                 teffs(2) = t(i)
              endif
           endif
        enddo

        if (dtlo.ge.10000. .or. dthi.ge.10000.) err = .true.
        

        return
        end




        subroutine get_gravs(t,g,n,teff,logg,glo,ghi,err)
c now find the closest gravities for given teff=thi
        real*8 teff,logg,glo,ghi,t(1000),g(1000)
        logical err

        err = .false.
           
        dglo = 10.
        dghi = 10.
        glo  = 10.
        ghi  = 10.

        do i = 1,n
           if (t(i).eq.teff) then
              if (g(i).eq.logg) then
                 glo = g(i)
                 ghi = g(i)
                 return
              elseif (g(i).lt.logg) then
                 if ( (logg-g(i)).le.dglo) then
                    dglo = logg-g(i)
                    glo = g(i)
                 endif
              elseif (g(i).gt.logg) then
                 if ( (g(i)-logg).le.dghi) then
                    dghi = g(i)-logg
                    ghi = g(i)
                 endif
              endif
           endif
        enddo

        if (glo.ge.10.0 .or. ghi.ge.10.) then
           err = .true.
        endif

        return
        end



      SUBROUTINE PRINT_MODEL(TEFF,LOGG,MODEL,FE,ERROR)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL ERROR
      CHARACTER*8 MODEL
      CHARACTER*(80) LINE
      CHARACTER*(80) FNAME
      REAL*8 TEMP,GRAV,LOGG,TEFF,FE

      ERROR = .TRUE.

      CALL ACCESS_MODEL_FILE(FE,FNAME,ERROR)

      OPEN(FILE=MODEL,UNIT=4)
      REWIND 4

      DO I = 1,1000000
         READ(11,'(A80)',END=900)LINE
         IF (LINE(:5).EQ.'MODEL') THEN
            READ(11,*)TEMP,GRAV
            IF (LOGG.EQ.GRAV  .AND. TEMP.EQ.TEFF) THEN
               ERROR = .FALSE.
               DO J = 1,1000000
                  READ(11,'(A80)',END=900)LINE
                  IF (LINE(:5).EQ.'MODEL') GOTO 900  
                  WRITE(4,'(A80)')LINE
               ENDDO
            ENDIF
         ENDIF
      ENDDO
 900  REWIND 4
      REWIND 11
      CLOSE(UNIT=11)
      CLOSE(UNIT=4)
      IF (ERROR) THEN
         WRITE(6,100)TEFF,LOGG,FNAME
         STOP
      ENDIF
 100  FORMAT('COULD NOT FIND FILE T = ',f6.0,'  LOGG = ',f6.2,' IN ',A80)
      RETURN
      END
      




      SUBROUTINE ACCESS_MODEL_FILE(FE,FNAME,ERROR)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/MODTYPE/ MTYPE
      CHARACTER*10 MTYPE,PREFIX*1
      LOGICAL ERROR
      CHARACTER*(5) MODEL*(8),FILES(15)
      CHARACTER*(80) DIR
      CHARACTER*(80) FNAME
      REAL*8 TEMP,GRAV,LOGG,TEFF,FE
      DATA DIR/'                                              '/
      DATA FILES /'sp100','sp075','sp050','sp025','sp000','sm025',
     %            'sm050','sm075','sm100','sm150','sm200','sm250',
     %            'sm300','sm400','sm500'/

      ERROR = .TRUE.
c
c Set the directory where models are located, currently only standard comp.
c also set prefix to indicate plane-parallel or 1D-spherical models 
c
      IF (MTYPE(:5).EQ.'MARCS') THEN
         DIR  = '/home/andy/fluff/marcs.models/standard/'
      ENDIF
      IF (MTYPE(6:7).EQ.'sp' .OR. MTYPE(6:7).EQ.'SP') THEN
        PREFIX ='p'
      ELSEIF (MTYPE(6:7).EQ.'ss' .OR. MTYPE(6:7).EQ.'SS') THEN
        PREFIX ='s'
      ENDIF
c
      IPOSN = INDEX(DIR,' ')
      IF (FE .EQ. 1.0D0) THEN  
         IFILE = 1
      ELSEIF (FE .EQ. 0.75D0) THEN
          IFILE = 2
      ELSEIF (FE .EQ.  0.50D0)  THEN
         IFILE = 3  
      ELSEIF (FE .EQ.  0.25D0)  THEN
         IFILE = 4  
      ELSEIF (FE .EQ.  0.00D0)  THEN
         IFILE = 5  
      ELSEIF (FE .EQ. -0.25D0)  THEN
         IFILE = 6  
      ELSEIF (FE .EQ. -0.50D0)  THEN  
         IFILE = 7  
      ELSEIF (FE .EQ. -0.75D0)  THEN  
         IFILE = 8  
      ELSEIF (FE .EQ. -1.0D0)  THEN  
         IFILE = 9  
      ELSEIF (FE .EQ. -1.5D0)  THEN  
         IFILE = 10  
      ELSEIF (FE .EQ. -2.0D0)  THEN  
         IFILE = 11  
      ELSEIF (FE .EQ. -2.5D0)  THEN  
         IFILE = 12  
      ELSEIF (FE .EQ. -3.0D0)  THEN  
         IFILE = 13  
      ELSEIF (FE .EQ. -4.0D0)  THEN  
         IFILE = 14  
      ELSEIF (FE .EQ. -5.0D0)  THEN  
         IFILE = 15
      ELSE  
         ERROR = .TRUE.
         WRITE(6,'a,f5.2,a')'MODEL METALLICITY [Fe/H] = ',FE,
     %                      ' NOT IN GRID'
         RETURN   
      ENDIF
      FNAME = DIR
      FNAME(IPOSN:) = FILES(IFILE)
      FNAME(IPOSN:IPOSN) = PREFIX
      OPEN(FILE=FNAME,UNIT=11,STATUS='OLD')
      REWIND 11

      RETURN
      END









      

      SUBROUTINE SET_UP_LINES_INPUT_FILE(MICRO,FE,TITLE)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/MODTYPE/ MTYPE
      CHARACTER*10 MTYPE
      CHARACTER*80 TITLE,LINE
      REAL*8 MICRO

c I must insert code for LINES to scale the SAD as required for RHOX models "

      OPEN(UNIT=4,FILE='model')
      OPEN(UNIT=7,FILE='MODEL')
      REWIND 4
      WRITE(7,'(A)')'KURTYPE'
      WRITE(7,'(A80)')TITLE

      IF (MTYPE.EQ.'KUROLD') THEN
         WRITE(7,'(13X,2H64)')
      ELSE
         WRITE(7,'(13X,2H56)')
      ENDIF

      WRITE(7,'(6H5000.0)')
      DO I = 1,1000000
         READ(4,'(A44)',END=900)LINE(:44)
         WRITE(7,'(A44)')LINE(:44)
      ENDDO
 900  CONTINUE
      WRITE(7,'(F13.2)')MICRO
      WRITE(7,'(6HNATOMS,8X,1H0,10X,F5.2)')FE
      WRITE(7,'(4HNMOL,10X,1H0)')
      REWIND 4
      CLOSE(UNIT=4,STATUS='DELETE')
      RETURN
      END


