program auto_ab4 c c This program is intended to take a list of atmosphere parameters c (Teff, logg, [Fe/H], and mic) with a set of EW input files, and performs c an abundance analysis using the Kurucz grid of model atmospheres. c c The list file format will be a30 for the EW input name, followed by a20 c output filename, folowed by free format Teff, logg, [Fe/H] and micro c implicit real*8 (a-h,o-z) character*30 listfile,ew_file,line*80 character*20 out_file common/mic/micro,omicro,slope,oslope,mic_step,sig_slope,slope1 real*8 micro,omicro,slope,oslope,mic_step,sig_slope,slope1 common/rinit/teff,logg,feh real*8 teff,logg,feh c c grab list name and open it c write(6,'(27h enter input list filename )') read(5,'(a30)')listfile open(unit=30,file=listfile,status='OLD') c c read EW filename in 1st 30 chars, and T, logg, [Fe/H] and micro c then call make kurucz to produce the model c do 100 i = 1,1000 read(30,'(a80)',end=999)line ew_file = line(:30) out_file = line(31:50) open(unit=31,file=ew_file,status='OLD') read(line(50:),*)teff,logg,feh,micro call make_kurucz(teff,logg,feh,micro) call auto(ew_file,sig_mic,sig_t) call save_output(micro,sig_mic,sig_t,out_file) 100 continue 999 continue stop end subroutine auto(ewfile,sig_mic,sig_t) implicit real*8 (a-h,o-z) c c This routine will converge on Teff and microturbulence, based on Fe I c line abundances, and then computes abundances of all lines. c c unit 7 = EW file c unit 8 = model file c common/rtemp/temp,otemp,tslope,otslope,t_step,sig_tslope,tslope1 real*8 temp,otemp,tslope,otslope,t_step,sig_tslope,tslope1 common/ltemp/firstt logical firstt common/rinit/teff,logg,feh real*8 teff,logg,feh common/mic/micro,omicro,slope,oslope,mic_step,sig_slope,slope1 real*8 micro,omicro,slope,oslope,mic_step,sig_slope,slope1 common/lmic/first logical first common/badfe/w(100),nbad real*8 w,sig_mic,sig_t integer nbad real*8 rew(500),ab(500),ep(500) integer nfe logical micro_ok,teff_ok character*30 ewfile call open_files(ewfile) call init_params call read_fe_list(ewfile) do 10 it = 1,20 call read_model mic_step = 0.5 first = .true. do 100 i = 1,10 call set_up_input(feh,micro) call system('~/bin/clines') call read_abunds(rew,ep,ab,nfe) c call fit_slope(rew,ab,nfe,slope) call line_fit(rew,ab,nfe,slope,sig_slope) if (micro_ok()) goto 9 100 continue call no_converge 9 continue call line_fit(ep,ab,nfe,tslope,sig_tslope) if (teff_ok()) then call make_kurucz(temp,logg,feh,micro) call setup_all(ewfile,micro) call system('~/bin/clines') call get_sigmic(sig_mic) call get_sigt(sig_t) return endif call make_kurucz(temp,logg,feh,micro) 10 continue call no_converge return end subroutine init_params c c a routine to initialise the Teff and microturbulence parameters c common/mic/micro,omicro,slope,oslope,mic_step,sig_slope,slope1 real*8 micro,omicro,slope,oslope,mic_step,sig_slope,slope1 common/lmic/first logical first common/rtemp/temp,otemp,tslope,otslope,t_step,sig_tslope,tslope1 real*8 temp,otemp,tslope,otslope,t_step,sig_tslope,tslope1 common/ltemp/firstt logical firstt common/rinit/teff,logg,feh real*8 teff,logg,feh first = .true. mic_step = 0.5 firstt = .true. t_step = 200.0 temp = teff return end subroutine open_files(ewfile) implicit real*8 (a-h,o-z) c c a routine to prompt for, and open input files c character*30 ewfile open(unit=13,file='auto.output') write(13,'(a30)')ewfile write(6,'(a30)')ewfile return end subroutine read_fe_list(ewfile) implicit real*8 (a-h,o-z) c c This routine will read through the EW line list and place in a special c file, 'autoab2.lines' only those lines to be used for microturbulence determination c c The routine could be added to in order to eliminate "bad" lines from c the Fe I list c character*80 line,id*5 character*30 ewfile real*8 ew,atom,wave logical bad_fe open(unit=9,file='autoab2.lines') open(unit=7,file=ewfile) c c read lines, and add to autoab2.lines file if Fe I less than log(RW)=-4.7 c write(9,'(a)')'Title Line' do 100 i = 1,10000 read(7,'(a80)',end=999)line read(line,'(a5,f10.3,f10.3,50x,f5.1)')id,wave,atom,ew if (atom.eq.26.0 .and. wave.ge.3850.) then rewlog = dlog10(ew/wave)-3.0 if (rewlog .le. -4.62 .and. c if (rewlog .le. -4.447 .and. rewlog .ge.-5.32 .and. c if (rewlog .le. -4.62 .and. rewlog .ge.-5.48 .and. c if (rewlog .le. -4.62 .and. rewlog .ge.-5.48) then c if (rewlog .le. -4.694 .and. rewlog .ge.-5.23 .and. % .not. bad_fe(wave) ) then write(9,'(a80)')line endif endif 100 continue 999 continue close (unit=9) rewind(unit=9) close (unit=7) return end subroutine read_model implicit real*8 (a-h,o-z) c c A routine to place the model in file 'autoab2.model', but first remove the IPR c control line c character*80 line open(unit=10,file='autoab2.model') open(unit=8,file='MODEL',status='old') do 100 i = 1,1000 read(8,'(a80)',end=999)line write(10,'(a80)')line 100 continue 999 continue rewind(unit=10) close(unit=10) close(unit=8) return end subroutine set_up_input(feh,micro) implicit real*8 (a-h,o-z) c c This routine will set up the INPUT file for LINES with the micro value c c head is set for a Kurucz atmosphere, with Kapref computed internally c real*8 micro integer imicro character*80 head,line data head/' '/ c c now copy the model and add microturbulence and [Fe/H] c open(unit=14,file='autoab2.model') open(unit=11,file='MODEL') do i = 1,68 read(14,'(a)')line write(11,'(a)')line enddo write(11,'(f13.2)')micro write(11,'(6HNATOMS,8x,1H0,5x,f10.2)')feh write(11,'(4HNMOL,10x,1H0)') close(unit=11) rewind(unit=11) close(unit=14) rewind(unit=14) c c now the lines c call system('cat autoab2.lines > INPUT') return end subroutine read_abunds(rew,ep,ab,nfe) implicit real*8 (a-h,o-z) c c reads the computed abundances from OUTPUT c real*8 rew(500),ab(500),ep(500) integer nfe open(unit=12,file='OUTPUT',status='OLD') rewind(unit=12) do 100 i = 1,500 read(12,'(5x,2f10.3,7x,2f8.3)',end=999)wave,ep(i),ew,ab(i) rew(i) = dlog10(ew/wave) - 3.00 100 continue write(6,'(45h more than 500 Fe I lines: change array size )') 999 continue nfe = i - 1 close(unit=12) rewind(unit=12) return end subroutine fit_slope(rew,ab,nfe,slope) c c This calls a minimum sum fitting routine, based on the simplex method. c As complex as this sounds, it doesn't give reasonable answers! c implicit real*8 (a-h,o-z) real*8 rew(500),ab(500),slope integer nfe real*8 q(502,4),toler,x(4),res(500),error,cu(2,502) integer k,l,m,n,klmd,klm2d,nklmd,n2d,iter,iu(2,502),s(500),kode k = 500 l = 0 m = 0 n = 2 klmd = 500 klm2d = 500 + 2 nklmd = 500 + 2 n2d = 4 toler = 1.0d-5 iter = 10 * 500 kode = 0 c c Now set up the q matrix c do 100 i = 1,nfe q(i,1) = rew(i) q(i,2) = 1.0 q(i,3) = ab(i) q(i,4) = 0.0 100 continue call cl1(k,l,m,n,klmd,klm2d,nklmd,n2d,q,kode,toler,iter,x,res, % error,cu,iu,s) if (kode .gt. 0) then write(5,'(29h minsum routine aborted code ,i2)')kode return endif slope = x(1) return end logical function micro_ok() c c A routine to see if the microturbulence has converged; if not, then set c micro to the next guess c implicit real*8 (a-h,o-z) common/mic/micro,omicro,slope,oslope,mic_step,sig_slope,slope1 real*8 micro,omicro,slope,oslope,mic_step,sig_slope,slope1 common/lmic/first logical first real*8 delta_mic micro_ok = .false. c c first time through need another abundance run to begin convergence on micro c write(13,'(9h micro = ,f10.3,15h Km/s; slope = ,g12.6, % 5h +/- ,f8.4)')micro,slope,sig_slope write(6,'(9h micro = ,f10.3,15h Km/s; slope = ,g12.6, % 5h +/- ,f8.4)')micro,slope,sig_slope if (first) then slope1 = slope first = .false. mic_step = slope/1.40 if (dabs(mic_step).lt.0.02) micro_ok = .true. omicro = micro oslope = slope micro = micro + mic_step else c c Test for convergence when slope changes sign; otherwise use delta_mic if c it is in range 0.04 to 0.50 Km/s c delta_mic = - slope * (omicro-micro)/(oslope-slope) if (slope/oslope .lt. 0.0 .or. dabs(delta_mic).le.0.02) then if (dabs(delta_mic).lt.0.02) then micro_ok = .true. endif else if (dabs(delta_mic) .gt. 0.50) then delta_mic = 0.50 * delta_mic/dabs(delta_mic) endif endif c c keep solution with smallest slope c if ( dabs(slope).lt.dabs(oslope) ) then omicro = micro oslope = slope endif c c set next micro to try c micro = micro + delta_mic endif if (micro_ok) then write(13,'(22h converged at micro = ,f10.3,6h Km/s)')micro write(6,'(22h converged at micro = ,f10.3,6h Km/s)')micro endif return end logical function teff_ok() c c A routine to see if the Teff has converged; if not, then set c Teff to the next guess c implicit real*8 (a-h,o-z) common/rtemp/temp,otemp,tslope,otslope,t_step,sig_tslope,tslope1 real*8 temp,otemp,tslope,otslope,t_step,sig_tslope,tslope1 common/ltemp/firstt logical firstt real*8 delta_t teff_ok = .false. c c first time through need another abundance run to begin convergence on teff c write(13,'(9h Teff = ,f10.1,12h K; slope = ,g12.6, % 5h +/- ,f8.4)')temp,tslope,sig_tslope write(6,'(9h Teff = ,f10.1,12h K; slope = ,g12.6, % 5h +/- ,f8.4)')temp,tslope,sig_tslope if (firstt) then tslope1 = tslope firstt = .false. t_step = tslope/0.00033 if (dabs(tslope).lt.0.003) teff_ok = .true. otemp = temp otslope = tslope temp = temp + t_step else c c Test for convergence when slope changes sign; otherwise use delta_t if c it is in range 30 to 300 K c delta_t = - tslope * (otemp-temp)/(otslope-tslope) if (tslope/otslope .lt. 0.0 .or. dabs(delta_t) .le. 10.0)then if (dabs(delta_t).lt.10.0) then teff_ok = .true. endif else if (dabs(delta_t) .gt. 300.0) then delta_mic = 300.0 * delta_t/dabs(delta_t) endif endif c c keep solution with smallest slope c if ( dabs(tslope).lt.dabs(otslope) ) then otemp = temp otslope = tslope endif c c set next Teff to try c temp = temp + delta_t endif if (teff_ok) then write(13,'(21h converged at Teff = ,f10.3,3h K)')temp write(6,'(21h converged at Teff = ,f10.3,3h K)')temp endif return end subroutine setup_all(ewfile,micro) implicit real*8 (a-h,o-z) c c This routine will set up the INPUT file for LINES with the micro value c for all lines in the input ewfile line list. c head is set for a Kurucz atmosphere, with Kapref computed internally c real*8 micro integer imicro character*80 head,command,ewfile*30 data head,command/' ',' '/ c c first set the header c head = ' 0 0 0 0 0 0 0 0 0 175 0 0 0' imicro = nint(100. * micro) write(head(38:40),'(i3)')imicro open(unit=11,file='INPUT') write(11,'(a80)')head close(unit=11) rewind(unit=11) c c now copy the model c c call system('cat autoab2.model > MODEL') c c now the lines c command = 'cat '//ewfile command(50:) = ' >> INPUT' call system(command) close(unit=11) rewind(unit=11) return end subroutine no_converge c c This routine just tells the user that convergence was not obtained c write(6,'(26h Convergence not obtained )') write(13,'(26h Convergence not obtained )') return end subroutine line_fit(rew,ab,nfe,slope,sig_slope) c c this is a general purpose least squares fitting program. The current set up c is for a polynomial fit of order specified by the user. To change the form c of the fit the user must change the subroutine "funcs". c implicit real*8 (a-h,o-z) real*8 x(1000),y(1000),sig(1000),a(50),rew(500),ab(500),slope real*8 covar(50,50),sig_slope,scale integer nfe c c set up data arrays, assuming equal weights c iord = 2 do 200 i = 1,nfe x(i) = rew(i) y(i) = ab(i) sig(i) = 1.0 200 continue call polylin(x,y,sig,iord,a,covar,chisq,nfe) scale = chisq/dble(nfe-2) sig_slope = dsqrt( scale * covar(2,2) ) slope = a(2) return end logical function bad_fe(wav) c c a logical function to identify bad Fe I lines, which should not be used c for abundance analysis c implicit real*8 (a-h,o-z) common/badfe/w(100),nbad real*8 w,wav integer nbad data nbad/40/ data w/ 3906.490, 3917.184, 3922.923, 4014.530, 4076.637, % 4132.908, 4134.685, 4147.675, 4174.917, 4181.764, 4182.387, % 4187.812, 4191.437, 4195.340, 4206.702, 4216.191, 4294.142, % 4307.912, 4337.055, 4427.317, 4430.622, 4459.100, 4461.660, % 4466.562, 4531.158, 4632.918, 4654.504, 4691.420, 4903.316, % 5041.076, 5041.763, 5051.640, 5079.745, 5110.435, 5254.953, % 5328.051, 5371.501, 5446.924, 5455.624, 6301.508, 60*0.0/ do 100 i = 1,nbad if (wav.eq.w(i)) then bad_fe = .true. return elseif (w(i).gt.wav) then bad_fe = .false. return endif 100 continue bad_fe = .false. return end subroutine save_output(micro,sig_mic,sig_t,out_file) implicit real*8 (a-h,o-z) real*8 teff,logg,feh,micro,sig_mic,sig_t character*20 out_file,command*80 common/rinit/teff,logg,feh common/rtemp/temp,otemp,tslope,otslope,t_step,sig_tslope,tslope1 real*8 temp,otemp,tslope,otslope,t_step,sig_tslope,tslope1 c open(unit=32,file=out_file,status='NEW') open(unit=32,file=out_file) write(32,'(a20)')out_file write(32,'(a,f9.3,f9.3,f6.3)')'Input Parameters:',teff,logg,feh write(6,'(7hTeff = ,f7.2,4h+/- ,f4.0)')temp,sig_t write(6,'(8hmicro = ,f7.3,4h+/- ,f7.2)')micro,sig_mic write(32,'(7hTeff = ,f7.2,4h+/- ,f4.0)')temp,sig_t write(32,'(8hmicro = ,f7.3,4h+/- ,f7.2)')micro,sig_mic close(unit=32) command = 'cat OUTPUT >> '//out_file call system(command) return end subroutine get_sigmic(sig_mic) c c a simple routine to compute the uncertainty of the microturbulent velocity c common/mic/micro,omicro,slope,oslope,mic_step,sig_slope,slope1 real*8 micro,omicro,slope,oslope,mic_step,sig_slope,slope1 real*8 sig_mic common/lmic/first logical first if (dabs(slope1).gt.0.20) then sig_mic = sig_slope * (micro-2.00)/(slope-slope1) else sig_mic = sig_slope / 1.40 endif return end subroutine get_sigt(sig_t) c c a simple routine to compute the uncertainty of the microturbulent velocity c common/rinit/teff,logg,feh real*8 teff,logg,feh common/rtemp/temp,otemp,tslope,otslope,t_step,sig_tslope,tslope1 real*8 temp,otemp,tslope,otslope,t_step,sig_tslope,tslope1 real*8 sig_t if (dabs(tslope1).gt.0.01) then sig_t = dabs(sig_tslope * (temp-teff)/(tslope-tslope1)) else sig_t = 0.00 endif return end c c These routines produce interpolated models in the format expected by c Chris Sneden's 1998 MOOG. c c When changing systems you need to change variable DIR in subroutine PRINT_MODEL c SUBROUTINE MAKE_KURUCZ(TEFF,LOGG,FE,MICRO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) REAL*8 TEFF,LOGG,FE,MICRO LOGICAL ERROR CHARACTER*80 TITLE TITLE = ' ' ERROR = .FALSE. WRITE(TITLE,100)TEFF,LOGG,FE,MICRO 100 FORMAT(F7.0,1h/,F6.1,1h/,F7.1,5x,7h mic = ,F6.4) CALL GET_MODEL(TEFF,LOGG,FE,ERROR) IF (ERROR) THEN WRITE(6,'(a)')'ERROR IN MAKE_KURUCZ SUBROUTINE ' STOP ENDIF CALL SET_UP_LINES_INPUT_FILE(MICRO,FE,TITLE) CLOSE(UNIT=7) RETURN END SUBROUTINE GET_MODEL(TEFF,LOGG,FE,ERROR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) REAL*8 LOGG,TEFF,FE,METAL(10),M1,M2 LOGICAL ERROR CHARACTER*8 M1FIL,M2FIL,MODEL DATA METAL/+1.0D0,+0.5D0,0.0D0,-0.5D0,-1.0D0,-1.5D0,-2.0D0,-2.5D0, % -3.0D0,-3.5D0/ DATA M1FIL,M2FIL,MODEL/'M1 ','M2 ','model '/ 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. -3.5D0) THEN WRITE(6,'(a)') % 'STAR IS MORE METAL POOR THAN -3.5: -3.5 MODEL USED' FE = -3.5D0 ENDIF IF (FE .EQ. +1.0D0 .OR. FE .EQ. +0.5D0 .OR. FE .EQ. 0.0D0 .OR. % FE .EQ. -0.5D0 .OR. FE .EQ. -1.0D0 .OR. FE .EQ. -1.5D0 .OR. % FE .EQ. -2.0D0 .OR. FE .EQ. -2.5D0 .OR. FE. EQ. -3.0D0 .OR. % FE .EQ. -3.5D0) THEN CALL COMPUTE_MODEL_ATMOSPHERE(TEFF,LOGG,FE,MODEL,ERROR) ELSE DO I = 2, 10 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) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*8 OUTFILE REAL*8 TEFF,TEFFS(2),LOGG,GRAVS(2),FE LOGICAL OKG,OKT,ERROR CALL CHECK_PARAMETERS(OKG,OKT,TEFF,LOGG,TEFFS,GRAVS,ERROR) IF (ERROR) RETURN IF (OKG .AND. OKT) THEN CALL PRINT_MODEL(TEFF,LOGG,OUTFILE,FE,ERROR) IF (ERROR) RETURN ELSEIF (OKG) THEN CALL PRINT_MODEL(TEFFS(1),LOGG,'MOD1 ',FE,ERROR) CALL PRINT_MODEL(TEFFS(2),LOGG,'MOD2 ',FE,ERROR) IF (ERROR) RETURN CALL KINTERP(TEFFS(1),TEFFS(2),TEFF,'MOD1 ','MOD2 ',OUTFILE) ELSEIF (OKT) THEN CALL PRINT_MODEL(TEFF,GRAVS(1),'MOD1 ',FE,ERROR) CALL PRINT_MODEL(TEFF,GRAVS(2),'MOD2 ',FE,ERROR) IF (ERROR) RETURN CALL KINTERP(GRAVS(1),GRAVS(2),LOGG,'MOD1 ','MOD2 ',OUTFILE) ELSE CALL PRINT_MODEL(TEFFS(1),GRAVS(1),'MOD1 ',FE,ERROR) CALL PRINT_MODEL(TEFFS(2),GRAVS(1),'MOD2 ',FE,ERROR) CALL PRINT_MODEL(TEFFS(1),GRAVS(2),'MOD3 ',FE,ERROR) CALL PRINT_MODEL(TEFFS(2),GRAVS(2),'MOD4 ',FE,ERROR) IF (ERROR) RETURN CALL KINTERP(GRAVS(1),GRAVS(2),LOGG,'MOD1 ','MOD3 ','MOD5 ') CALL KINTERP(GRAVS(1),GRAVS(2),LOGG,'MOD2 ','MOD4 ','MOD6 ') CALL KINTERP(TEFFS(1),TEFFS(2),TEFF,'MOD5 ','MOD6 ',OUTFILE) ENDIF RETURN END SUBROUTINE CHECK_PARAMETERS(OKG,OKT,TEFF,GRAV,TEFFS,GRAVS,ERROR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) LOGICAL OKG,OKT,ERROR INTEGER I REAL*8 TEFF,TEFFS(2),GRAV,GRAVS(2) IF (TEFF.LT.3500.0D0 .OR. TEFF.GT.10000.0D0) THEN WRITE(6,'(a)')'ERROR: TEMPERATURE IS OUT OF BOUNDS' STOP ENDIF c check to see if Teff is a multiple of 250K I = INT(TEFF/250.0D0) IF (DABS(TEFF - DBLE(I)*250.0D0) .LT. 0.5D0) THEN TEFFS(1) = DBLE(I)*250.0D0 TEFFS(2) = DBLE(I)*250.0D0 TEFF = DBLE(NINT(TEFF)) OKT = .TRUE. ELSE TEFFS(1) = DBLE(I)*250.0D0 TEFFS(2) = DBLE(I+1)*250.0D0 OKT = .FALSE. ENDIF c Now logg IF (GRAV.LT.0.0D0 .OR. GRAV.GT.5.0D0) THEN WRITE(6,'(a)')'ERROR: GRAVITY IS OUT OF BOUNDS' STOP ENDIF c check to see if logg is a multiple of 0.5 I = INT(GRAV/0.5D0) IF (DABS(GRAV - DBLE(I)*0.5D0) .LT. 0.005D0) THEN GRAVS(1) = DBLE(I)*0.5D0 GRAVS(2) = DBLE(I)*0.5D0 GRAV = DBLE(NINT(GRAV*2.0))*0.5 OKG = .TRUE. ELSE GRAVS(1) = DBLE(I)*0.5D0 GRAVS(2) = DBLE(I+1)*0.5D0 OKG = .FALSE. ENDIF RETURN END SUBROUTINE PRINT_MODEL(TEFF,LOGG,MODEL,FE,ERROR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) LOGICAL ERROR CHARACTER*(4) MODEL*(8),FILES(10) CHARACTER*(80) DIR,LINE CHARACTER*(80) FNAME REAL*8 TEMP,GRAV,LOGG,TEFF,FE DATA DIR/' '/ DATA FILES /'ap10','ap05','ap00','am05','am10','am15','am20', % 'am25', 'am30','am35'/ ERROR = .TRUE. DIR = '/home/fluff/andy/kurucz.models/' IPOSN = INDEX(DIR,' ') IF (FE .EQ. 1.0D0) THEN IFILE = 1 ELSEIF (FE .EQ. 0.5D0) THEN IFILE = 2 ELSEIF (FE .EQ. 0.0D0) THEN IFILE = 3 ELSEIF (FE .EQ. -0.5D0) THEN IFILE = 4 ELSEIF (FE .EQ. -1.0D0) THEN IFILE = 5 ELSEIF (FE .EQ. -1.5D0) THEN IFILE = 6 ELSEIF (FE .EQ. -2.0D0) THEN IFILE = 7 ELSEIF (FE .EQ. -2.5D0) THEN IFILE = 8 ELSEIF (FE .EQ. -3.0D0) THEN IFILE = 9 ELSEIF (FE .EQ. -3.5D0) THEN IFILE = 10 ELSE ERROR = .TRUE. WRITE(6,'(a)') 'FE/H NOT +1.0, +0.5, 0.0, -0.5, -1.0, -1.5, -2.0, -2.5, -3.0, OR -3.5' RETURN ENDIF FNAME = DIR FNAME(IPOSN:) = FILES(IFILE) OPEN(FILE=FNAME,UNIT=11,STATUS='OLD') OPEN(FILE=MODEL,UNIT=4) REWIND 11 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 STOP ENDIF 100 FORMAT('COULD NOT FIND FILE T = ',f6.0,' LOGG = ',f6.2) RETURN END SUBROUTINE SET_UP_LINES_INPUT_FILE(MICRO,FE,TITLE) IMPLICIT DOUBLE PRECISION (A-H,O-Z) 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 WRITE(7,'(13X,2H64)') 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